added SpikeTrain class, changed isi_distance

spike trains are now represented as SpikeTrain objects consisting of the
spike times and the interval edges.
The implementation of the ISI-distance has been modified accordingly.
The SPIKE-distance and SPIKE-Synchronization are still to be done.

6 files changed, 156 insertions, 74 deletions

diff --git a/pyspike/SpikeTrain.py b/pyspike/SpikeTrain.py
new file mode 100644
index 0000000..4760014
--- /dev/null
+++ b/pyspike/SpikeTrain.py
@@ -0,0 +1,34 @@
+""" Module containing the class representing spike trains for PySpike.
+
+Copyright 2015, Mario Mulansky <mario.mulansky@gmx.net>
+
+Distributed under the BSD License
+"""
+
+import numpy as np
+import collections
+
+
+class SpikeTrain:
+    """ Class representing spike trains for the PySpike Module."""
+
+    def __init__(self, spike_times, interval):
+        """ Constructs the SpikeTrain
+        :param spike_times: ordered array of spike times.
+        :param interval: interval defining the edges of the spike train.
+        Given as a pair of floats (T0, T1) or a single float T1, where T0=0 is
+        assumed.
+        """
+
+        # TODO: sanity checks
+        self.spikes = np.array(spike_times)
+
+        # check if interval is as sequence
+        if not isinstance(interval, collections.Sequence):
+            # treat value as end time and assume t_start = 0
+            self.t_start = 0.0
+            self.t_end = interval
+        else:
+            # extract times from sequence
+            self.t_start = interval[0]
+            self.t_end = interval[1]
diff --git a/pyspike/__init__.py b/pyspike/__init__.py
index 4d3f9f6..76e58a1 100644
--- a/pyspike/__init__.py
+++ b/pyspike/__init__.py
@@ -5,12 +5,13 @@ Distributed under the BSD License
 """
 
 __all__ = ["isi_distance", "spike_distance", "spike_sync", "psth",
-           "spikes", "PieceWiseConstFunc", "PieceWiseLinFunc",
+           "spikes", "SpikeTrain", "PieceWiseConstFunc", "PieceWiseLinFunc",
            "DiscreteFunc"]
 
 from PieceWiseConstFunc import PieceWiseConstFunc
 from PieceWiseLinFunc import PieceWiseLinFunc
 from DiscreteFunc import DiscreteFunc
+from SpikeTrain import SpikeTrain
 
 from isi_distance import isi_profile, isi_distance, isi_profile_multi,\
     isi_distance_multi, isi_distance_matrix
diff --git a/pyspike/cython/cython_distance.pyx b/pyspike/cython/cython_distance.pyx
index 2834ca5..1d652ee 100644
--- a/pyspike/cython/cython_distance.pyx
+++ b/pyspike/cython/cython_distance.pyx
@@ -42,57 +42,82 @@ ctypedef np.float_t DTYPE_t
 ############################################################
 # isi_distance_cython
 ############################################################
-def isi_distance_cython(double[:] s1, 
-                        double[:] s2):
+def isi_distance_cython(double[:] s1, double[:] s2,
+                        double t_start, double t_end):
 
     cdef double[:] spike_events
     cdef double[:] isi_values
     cdef int index1, index2, index
     cdef int N1, N2
     cdef double nu1, nu2
-    N1 = len(s1)-1
-    N2 = len(s2)-1
+    N1 = len(s1)
+    N2 = len(s2)
+
+    spike_events = np.empty(N1+N2+2)
+    # the values have one entry less as they are defined at the intervals
+    isi_values = np.empty(N1+N2+1)
+
+    # first x-value of the profile
+    spike_events[0] = t_start
+
+    # first interspike interval - check if a spike exists at the start time
+    if s1[0] > t_start:
+        nu1 = s1[0] - t_start
+        index1 = -1
+    else:
+        nu1 = s1[1]-s1[0]
+        index1 = 0
+
+    if s2[0] > t_start:
+        nu2 = s2[0] - t_start
+        index2 = -1
+    else:
+        nu2 = s2[1]-s2[0]
+        index2 = 0
 
-    nu1 = s1[1]-s1[0]
-    nu2 = s2[1]-s2[0]
-    spike_events = np.empty(N1+N2)
-    spike_events[0] = s1[0]
-    # the values have one entry less - the number of intervals between events
-    isi_values = np.empty(N1+N2-1)
+    isi_values[0] = fabs(nu1-nu2)/fmax(nu1, nu2)
+    index = 1
 
     with nogil: # release the interpreter to allow multithreading
-        isi_values[0] = fabs(nu1-nu2)/fmax(nu1, nu2)
-        index1 = 0
-        index2 = 0
-        index = 1
-        while True:
-            # check which spike is next - from s1 or s2
-            if s1[index1+1] < s2[index2+1]:
+        while index1+index2 < N1+N2-2:
+            # check which spike is next, only if there are spikes left in 1
+            # next spike in 1 is earlier, or there are no spikes left in 2
+            if (index1 < N1-1) and ((index2 == N2-1) or
+                                    (s1[index1+1] < s2[index2+1])):
                 index1 += 1
-                # break condition relies on existence of spikes at T_end
-                if index1 >= N1:
-                    break
                 spike_events[index] = s1[index1]
-                nu1 = s1[index1+1]-s1[index1]
-            elif s1[index1+1] > s2[index2+1]:
+                if index1 < N1-1:
+                    nu1 = s1[index1+1]-s1[index1]
+                else:
+                    nu1 = t_end-s1[index1]
+            elif (index2 < N2-1) and ((index1 == N1-1) or
+                                      (s1[index1+1] > s2[index2+1])):
                 index2 += 1
-                if index2 >= N2:
-                    break
                 spike_events[index] = s2[index2]
-                nu2 = s2[index2+1]-s2[index2]
+                if index2 < N2-1:
+                    nu2 = s2[index2+1]-s2[index2]
+                else:
+                    nu2 = t_end-s2[index2]
             else: # s1[index1+1] == s2[index2+1]
                 index1 += 1
                 index2 += 1
-                if (index1 >= N1) or (index2 >= N2):
-                    break
                 spike_events[index] = s1[index1]
-                nu1 = s1[index1+1]-s1[index1]
-                nu2 = s2[index2+1]-s2[index2]            
+                if index1 < N1-1:
+                    nu1 = s1[index1+1]-s1[index1]
+                else:
+                    nu1 = t_end-s1[index1]
+                if index2 < N2-1:
+                    nu2 = s2[index2+1]-s2[index2]
+                else:
+                    nu2 = t_end-s2[index2]
             # compute the corresponding isi-distance
             isi_values[index] = fabs(nu1 - nu2) / fmax(nu1, nu2)
             index += 1
         # the last event is the interval end
-        spike_events[index] = s1[N1]
+        if spike_events[index-1] == t_end:
+            index -= 1
+        else:
+            spike_events[index] = t_end
     # end nogil
 
     return spike_events[:index+1], isi_values[:index]
diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py
index 749507a..4c37236 100644
--- a/pyspike/cython/python_backend.py
+++ b/pyspike/cython/python_backend.py
@@ -15,50 +15,72 @@ import numpy as np
 ############################################################
 # isi_distance_python
 ############################################################
-def isi_distance_python(s1, s2):
+def isi_distance_python(s1, s2, t_start, t_end):
     """ Plain Python implementation of the isi distance.
     """
-    # compute the interspike interval
-    nu1 = s1[1:] - s1[:-1]
-    nu2 = s2[1:] - s2[:-1]
+    N1 = len(s1)
+    N2 = len(s2)
 
     # compute the isi-distance
-    spike_events = np.empty(len(nu1) + len(nu2))
-    spike_events[0] = s1[0]
+    spike_events = np.empty(N1+N2+2)
+    spike_events[0] = t_start
     # the values have one entry less - the number of intervals between events
     isi_values = np.empty(len(spike_events) - 1)
-    # add the distance of the first events
-    # isi_values[0] = nu1[0]/nu2[0] - 1.0 if nu1[0] <= nu2[0] \
-    #                 else 1.0 - nu2[0]/nu1[0]
-    isi_values[0] = abs(nu1[0] - nu2[0]) / max(nu1[0], nu2[0])
-    index1 = 0
-    index2 = 0
+    if s1[0] > t_start:
+        nu1 = s1[0] - t_start
+        index1 = -1
+    else:
+        nu1 = s1[1] - s1[0]
+        index1 = 0
+    if s2[0] > t_start:
+        nu2 = s2[0] - t_start
+        index2 = -1
+    else:
+        nu2 = s2[1] - s2[0]
+        index2 = 0
+
+    isi_values[0] = abs(nu1 - nu2) / max(nu1, nu2)
     index = 1
-    while True:
+    while index1+index2 < N1+N2-2:
         # check which spike is next - from s1 or s2
-        if s1[index1+1] < s2[index2+1]:
+        if (index1 < N1-1) and (index2 == N2-1 or s1[index1+1] < s2[index2+1]):
             index1 += 1
-            # break condition relies on existence of spikes at T_end
-            if index1 >= len(nu1):
-                break
             spike_events[index] = s1[index1]
-        elif s1[index1+1] > s2[index2+1]:
+            if index1 < N1-1:
+                nu1 = s1[index1+1]-s1[index1]
+            else:
+                nu1 = t_end-s1[index1]
+
+        elif (index2 < N2-1) and (index1 == N1-1 or
+                                  s1[index1+1] > s2[index2+1]):
             index2 += 1
-            if index2 >= len(nu2):
-                break
             spike_events[index] = s2[index2]
+            if index2 < N2-1:
+                nu2 = s2[index2+1]-s2[index2]
+            else:
+                nu2 = t_end-s2[index2]
+
         else:  # s1[index1 + 1] == s2[index2 + 1]
             index1 += 1
             index2 += 1
-            if (index1 >= len(nu1)) or (index2 >= len(nu2)):
-                break
             spike_events[index] = s1[index1]
+            if index1 < N1-1:
+                nu1 = s1[index1+1]-s1[index1]
+            else:
+                nu1 = t_end-s1[index1]
+            if index2 < N2-1:
+                nu2 = s2[index2+1]-s2[index2]
+            else:
+                nu2 = t_end-s2[index2]
         # compute the corresponding isi-distance
-        isi_values[index] = abs(nu1[index1] - nu2[index2]) / \
-            max(nu1[index1], nu2[index2])
+        isi_values[index] = abs(nu1 - nu2) / \
+            max(nu1, nu2)
         index += 1
     # the last event is the interval end
-    spike_events[index] = s1[-1]
+    if spike_events[index-1] == t_end:
+        index -= 1
+    else:
+        spike_events[index] = t_end
     # use only the data added above
     # could be less than original length due to equal spike times
     return spike_events[:index + 1], isi_values[:index]
diff --git a/pyspike/isi_distance.py b/pyspike/isi_distance.py
index c2ef8e8..a34e135 100644
--- a/pyspike/isi_distance.py
+++ b/pyspike/isi_distance.py
@@ -14,23 +14,25 @@ from pyspike.generic import _generic_profile_multi, _generic_distance_matrix
 ############################################################
 # isi_profile
 ############################################################
-def isi_profile(spikes1, spikes2):
+def isi_profile(spike_train1, spike_train2):
     """ Computes the isi-distance profile :math:`S_{isi}(t)` of the two given
     spike trains. Retruns the profile as a PieceWiseConstFunc object. The S_isi
     values are defined positive S_isi(t)>=0.  The spike trains are expected
     to have auxiliary spikes at the beginning and end of the interval. Use the
     function add_auxiliary_spikes to add those spikes to the spike train.
 
-    :param spikes1: ordered array of spike times with auxiliary spikes.
-    :param spikes2: ordered array of spike times with auxiliary spikes.
+    :param spike_train1: First spike train.
+    :type spike_train1: :class:`pyspike.SpikeTrain`
+    :param spike_train2: Second spike train.
+    :type spike_train2: `SpikeTrain`
     :returns: The isi-distance profile :math:`S_{isi}(t)`
     :rtype: :class:`pyspike.function.PieceWiseConstFunc`
 
     """
-    # check for auxiliary spikes - first and last spikes should be identical
-    assert spikes1[0] == spikes2[0], \
+    # check whether the spike trains are defined for the same interval
+    assert spike_train1.t_start == spike_train2.t_start, \
         "Given spike trains seems not to have auxiliary spikes!"
-    assert spikes1[-1] == spikes2[-1], \
+    assert spike_train1.t_end == spike_train2.t_end, \
         "Given spike trains seems not to have auxiliary spikes!"
 
     # load cython implementation
@@ -45,7 +47,8 @@ Falling back to slow python backend.")
         from cython.python_backend import isi_distance_python \
             as isi_distance_impl
 
-    times, values = isi_distance_impl(spikes1, spikes2)
+    times, values = isi_distance_impl(spike_train1.spikes, spike_train2.spikes,
+                                      spike_train1.t_start, spike_train1.t_end)
     return PieceWiseConstFunc(times, values)
 
 
diff --git a/test/test_distance.py b/test/test_distance.py
index ba19f5e..b54e908 100644
--- a/test/test_distance.py
+++ b/test/test_distance.py
@@ -15,12 +15,13 @@ from numpy.testing import assert_equal, assert_almost_equal, \
     assert_array_almost_equal
 
 import pyspike as spk
+from pyspike import SpikeTrain
 
 
 def test_isi():
     # generate two spike trains:
-    t1 = np.array([0.2, 0.4, 0.6, 0.7])
-    t2 = np.array([0.3, 0.45, 0.8, 0.9, 0.95])
+    t1 = SpikeTrain([0.2, 0.4, 0.6, 0.7], 1.0)
+    t2 = SpikeTrain([0.3, 0.45, 0.8, 0.9, 0.95], 1.0)
 
     # pen&paper calculation of the isi distance
     expected_times = [0.0, 0.2, 0.3, 0.4, 0.45, 0.6, 0.7, 0.8, 0.9, 0.95, 1.0]
@@ -32,8 +33,6 @@ def test_isi():
     expected_isi_val = sum((expected_times[1:] - expected_times[:-1]) *
                            expected_isi)/(expected_times[-1]-expected_times[0])
 
-    t1 = spk.add_auxiliary_spikes(t1, 1.0)
-    t2 = spk.add_auxiliary_spikes(t2, 1.0)
     f = spk.isi_profile(t1, t2)
 
     # print("ISI: ", f.y)
@@ -44,8 +43,8 @@ def test_isi():
     assert_equal(spk.isi_distance(t1, t2), expected_isi_val)
 
     # check with some equal spike times
-    t1 = np.array([0.2, 0.4, 0.6])
-    t2 = np.array([0.1, 0.4, 0.5, 0.6])
+    t1 = SpikeTrain([0.2, 0.4, 0.6], [0.0, 1.0])
+    t2 = SpikeTrain([0.1, 0.4, 0.5, 0.6], [0.0, 1.0])
 
     expected_times = [0.0, 0.1, 0.2, 0.4, 0.5, 0.6, 1.0]
     expected_isi = [0.1/0.2, 0.1/0.3, 0.1/0.3, 0.1/0.2, 0.1/0.2, 0.0/0.5]
@@ -55,8 +54,6 @@ def test_isi():
     expected_isi_val = sum((expected_times[1:] - expected_times[:-1]) *
                            expected_isi)/(expected_times[-1]-expected_times[0])
 
-    t1 = spk.add_auxiliary_spikes(t1, 1.0)
-    t2 = spk.add_auxiliary_spikes(t2, 1.0)
     f = spk.isi_profile(t1, t2)
 
     assert_equal(f.x, expected_times)
@@ -318,6 +315,6 @@ def test_multi_variate_subsets():
 
 if __name__ == "__main__":
     test_isi()
-    test_spike()
-    test_multi_isi()
-    test_multi_spike()
+    # test_spike()
+    # test_multi_isi()
+    # test_multi_spike()