Fix incorrect integrals in PieceWiseConstFunc (#36)

* Add (some currently failing) tests for PieceWiseConstFunc.integral

* Fix implementation of PieceWiseConstFunc.integral

Just by adding a special condition for when we are only taking an
integral "between" two edges of a PieceWiseConstFunc
All tests now pass.
Fixes #33.

* Add PieceWiseConstFunc.integral tests for ValueError

* Add testing bounds of integral
* Raise ValueError in function implementation

2 files changed, 72 insertions, 10 deletions

diff --git a/pyspike/PieceWiseConstFunc.py b/pyspike/PieceWiseConstFunc.py
index 5ce5f27..17fdd3f 100644
--- a/pyspike/PieceWiseConstFunc.py
+++ b/pyspike/PieceWiseConstFunc.py
@@ -129,19 +129,31 @@ class PieceWiseConstFunc(object):
             # no interval given, integrate over the whole spike train
             a = np.sum((self.x[1:]-self.x[:-1]) * self.y)
         else:
+            if interval[0]>interval[1]:
+                raise ValueError("Invalid averaging interval: interval[0]>=interval[1]")
+            if interval[0]<self.x[0]:
+                raise ValueError("Invalid averaging interval: interval[0]<self.x[0]")
+            if interval[1]>self.x[-1]:
+                raise ValueError("Invalid averaging interval: interval[0]<self.x[-1]")
             # find the indices corresponding to the interval
             start_ind = np.searchsorted(self.x, interval[0], side='right')
             end_ind = np.searchsorted(self.x, interval[1], side='left')-1
-            assert start_ind > 0 and end_ind < len(self.x), \
-                "Invalid averaging interval"
-            # first the contribution from between the indices
-            a = np.sum((self.x[start_ind+1:end_ind+1] -
-                        self.x[start_ind:end_ind]) *
-                       self.y[start_ind:end_ind])
-            # correction from start to first index
-            a += (self.x[start_ind]-interval[0]) * self.y[start_ind-1]
-            # correction from last index to end
-            a += (interval[1]-self.x[end_ind]) * self.y[end_ind]
+            if start_ind > end_ind:
+                # contribution from between two closest edges
+                a = (self.x[start_ind]-self.x[end_ind]) * self.y[end_ind]
+                # minus the part that is not within the interval
+                a -= ((interval[0]-self.x[end_ind])+(self.x[start_ind]-interval[1])) * self.y[end_ind]
+            else:
+                assert start_ind > 0 and end_ind < len(self.x), \
+                    "Invalid averaging interval"
+                # first the contribution from between the indices
+                a = np.sum((self.x[start_ind+1:end_ind+1] -
+                            self.x[start_ind:end_ind]) *
+                           self.y[start_ind:end_ind])
+                # correction from start to first index
+                a += (self.x[start_ind]-interval[0]) * self.y[start_ind-1]
+                # correction from last index to end
+                a += (interval[1]-self.x[end_ind]) * self.y[end_ind]
         return a
 
     def avrg(self, interval=None):
diff --git a/test/test_function.py b/test/test_function.py
index 92d378d..ddeb4b1 100644
--- a/test/test_function.py
+++ b/test/test_function.py
@@ -10,6 +10,7 @@ Distributed under the BSD License
 from __future__ import print_function
 import numpy as np
 from copy import copy
+from nose.tools import raises
 from numpy.testing import assert_equal, assert_almost_equal, \
     assert_array_equal, assert_array_almost_equal
 
@@ -49,6 +50,8 @@ def test_pwc():
     assert_almost_equal(a, (0.5-0.5+0.5*1.5+1.0*0.75)/3.0, decimal=16)
     a = f.avrg([1.5, 3.5])
     assert_almost_equal(a, (-0.5*0.5+0.5*1.5+1.0*0.75)/2.0, decimal=16)
+    a = f.avrg([1.0, 2.0])
+    assert_almost_equal(a, (1.0*-0.5)/1.0, decimal=16)
     a = f.avrg([1.0, 3.5])
     assert_almost_equal(a, (-0.5*1.0+0.5*1.5+1.0*0.75)/2.5, decimal=16)
     a = f.avrg([1.0, 4.0])
@@ -120,6 +123,53 @@ def test_pwc_avrg():
     assert_array_almost_equal(f1.x, x_expected, decimal=16)
     assert_array_almost_equal(f1.y, y_expected, decimal=16)
 
+def test_pwc_integral():
+    # some random data
+    x = [0.0, 1.0, 2.0, 2.5, 4.0]
+    y = [1.0, -0.5, 1.5, 0.75]
+    f1 = spk.PieceWiseConstFunc(x, y)
+
+    # test full interval
+    full = 1.0*1.0 + 1.0*-0.5 + 0.5*1.5 + 1.5*0.75;
+    assert_equal(f1.integral(), full)
+    assert_equal(f1.integral((np.min(x),np.max(x))), full)
+    # test part interval, spanning an edge
+    assert_equal(f1.integral((0.5,1.5)), 0.5*1.0 + 0.5*-0.5)
+    # test part interval, just over two edges
+    assert_almost_equal(f1.integral((1.0-1e-16,2+1e-16)), 1.0*-0.5, decimal=16)
+    # test part interval, between two edges
+    assert_equal(f1.integral((1.0,2.0)), 1.0*-0.5)
+    assert_equal(f1.integral((1.2,1.7)), (1.7-1.2)*-0.5)
+    # test part interval, start to before and after edge
+    assert_equal(f1.integral((0.0,0.7)), 0.7*1.0)
+    assert_equal(f1.integral((0.0,1.1)), 1.0*1.0+0.1*-0.5)
+    # test part interval, before and after edge till end
+    assert_equal(f1.integral((2.6,4.0)), (4.0-2.6)*0.75)
+    assert_equal(f1.integral((2.4,4.0)), (2.5-2.4)*1.5+(4-2.5)*0.75)
+
+@raises(ValueError)
+def test_pwc_integral_bad_bounds_inv():
+    # some random data
+    x = [0.0, 1.0, 2.0, 2.5, 4.0]
+    y = [1.0, -0.5, 1.5, 0.75]
+    f1 = spk.PieceWiseConstFunc(x, y)
+    f1.integral((3,2))
+
+@raises(ValueError)
+def test_pwc_integral_bad_bounds_oob_1():
+    # some random data
+    x = [0.0, 1.0, 2.0, 2.5, 4.0]
+    y = [1.0, -0.5, 1.5, 0.75]
+    f1 = spk.PieceWiseConstFunc(x, y)
+    f1.integral((1,6))
+
+@raises(ValueError)
+def test_pwc_integral_bad_bounds_oob_2():
+    # some random data
+    x = [0.0, 1.0, 2.0, 2.5, 4.0]
+    y = [1.0, -0.5, 1.5, 0.75]
+    f1 = spk.PieceWiseConstFunc(x, y)
+    f1.integral((-1,3))
 
 def test_pwl():
     x = [0.0, 1.0, 2.0, 2.5, 4.0]