From 27aa30d1fdb830a04b608c702cf7b616115eeb50 Mon Sep 17 00:00:00 2001
From: Mario Mulansky <mario.mulansky@gmx.net>
Date: Wed, 22 Apr 2015 18:18:30 +0200
Subject: 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.
---
 pyspike/cython/python_backend.py | 72 ++++++++++++++++++++++++++--------------
 1 file changed, 47 insertions(+), 25 deletions(-)

(limited to 'pyspike/cython/python_backend.py')

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]
-- 
cgit v1.2.3


From ed85a9b72edcb7bba6ae1105e213b3b0a2f78d3a Mon Sep 17 00:00:00 2001
From: Mario Mulansky <mario.mulansky@gmx.net>
Date: Fri, 24 Apr 2015 00:49:16 +0200
Subject: changed spike distance to use new SpikeTrain class

---
 pyspike/cython/cython_distance.pyx | 199 +++++++++++++++++++++++++-----------
 pyspike/cython/python_backend.py   | 203 ++++++++++++++++++++++++-------------
 pyspike/spike_distance.py          |  37 ++++---
 test/test_distance.py              |  34 ++++---
 4 files changed, 313 insertions(+), 160 deletions(-)

(limited to 'pyspike/cython/python_backend.py')

diff --git a/pyspike/cython/cython_distance.pyx b/pyspike/cython/cython_distance.pyx
index 1d652ee..7999e0a 100644
--- a/pyspike/cython/cython_distance.pyx
+++ b/pyspike/cython/cython_distance.pyx
@@ -131,21 +131,30 @@ cdef inline double get_min_dist_cython(double spike_time,
                                        # use memory view to ensure inlining
                                        # np.ndarray[DTYPE_t,ndim=1] spike_train,
                                        int N,
-                                       int start_index=0) nogil:
+                                       int start_index,
+                                       double t_start, double t_end) nogil:
     """ Returns the minimal distance |spike_time - spike_train[i]| 
     with i>=start_index.
     """
     cdef double d, d_temp
-    d = fabs(spike_time - spike_train[start_index])
-    start_index += 1
+    # start with the distance to the start time
+    d = fabs(spike_time - t_start)
+    if start_index < 0:
+        start_index = 0
     while start_index < N:
         d_temp = fabs(spike_time - spike_train[start_index])
         if d_temp > d:
-            break
+            return d
         else:
             d = d_temp
         start_index += 1
-    return d
+
+    # finally, check the distance to end time
+    d_temp = fabs(t_end - spike_time)
+    if d_temp > d:
+        return d
+    else:
+        return d_temp
 
 
 ############################################################
@@ -160,96 +169,162 @@ cdef inline double isi_avrg_cython(double isi1, double isi2) nogil:
 ############################################################
 # spike_distance_cython
 ############################################################
-def spike_distance_cython(double[:] t1, 
-                          double[:] t2):
+def spike_distance_cython(double[:] t1, double[:] t2,
+                          double t_start, double t_end):
 
     cdef double[:] spike_events
     cdef double[:] y_starts
     cdef double[:] y_ends
 
     cdef int N1, N2, index1, index2, index
-    cdef double dt_p1, dt_p2, dt_f1, dt_f2, isi1, isi2, s1, s2
+    cdef double t_p1, t_f1, t_p2, t_f2, dt_p1, dt_p2, dt_f1, dt_f2
+    cdef double isi1, isi2, s1, s2
 
     N1 = len(t1)
     N2 = len(t2)
 
-    spike_events = np.empty(N1+N2-2)
-    spike_events[0] = t1[0]
+    spike_events = np.empty(N1+N2+2)
+
     y_starts = np.empty(len(spike_events)-1)
     y_ends = np.empty(len(spike_events)-1)
 
     with nogil: # release the interpreter to allow multithreading
-        index1 = 0
-        index2 = 0
-        index = 1
-        dt_p1 = 0.0
-        dt_f1 = get_min_dist_cython(t1[1], t2, N2, 0)
-        dt_p2 = 0.0
-        dt_f2 = get_min_dist_cython(t2[1], t1, N1, 0)
-        isi1 = max(t1[1]-t1[0], t1[2]-t1[1])
-        isi2 = max(t2[1]-t2[0], t2[2]-t2[1])
-        s1 = dt_f1*(t1[1]-t1[0])/isi1
-        s2 = dt_f2*(t2[1]-t2[0])/isi2
+        spike_events[0] = t_start
+        t_p1 = t_start
+        t_p2 = t_start
+        if t1[0] > t_start:
+            # dt_p1 = t2[0]-t_start
+            dt_p1 = 0.0
+            t_f1 = t1[0]
+            dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_start, t_end)
+            isi1 = fmax(t_f1-t_start, t1[1]-t1[0])
+            s1 = dt_f1*(t_f1-t_start)/isi1
+            index1 = -1
+        else:
+            dt_p1 = 0.0
+            t_f1 = t1[1]
+            dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_start, t_end)
+            isi1 = t1[1]-t1[0]
+            s1 = dt_p1
+            index1 = 0
+        if t2[0] > t_start:
+            # dt_p1 = t2[0]-t_start
+            dt_p2 = 0.0
+            t_f2 = t2[0]
+            dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_start, t_end)
+            isi2 = fmax(t_f2-t_start, t2[1]-t2[0])
+            s2 = dt_f2*(t_f2-t_start)/isi2
+            index2 = -1
+        else:
+            dt_p2 = 0.0
+            t_f2 = t2[1]
+            dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_start, t_end)
+            isi2 = t2[1]-t2[0]
+            s2 = dt_p2
+            index2 = 0
+
         y_starts[0] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2)
-        while True:
+        index = 1
+
+        while index1+index2 < N1+N2-2:
             # print(index, index1, index2)
-            if t1[index1+1] < t2[index2+1]:
+            if (index1 < N1-1) and (t_f1 < t_f2 or index2 == N2-1):
                 index1 += 1
-                # break condition relies on existence of spikes at T_end
-                if index1+1 >= N1:
-                    break
-                spike_events[index] = t1[index1]
                 # first calculate the previous interval end value
-                dt_p1 = dt_f1 # the previous time now was the following time before
+                # the previous time now was the following time before:
+                dt_p1 = dt_f1  
+                t_p1 = t_f1    # t_p1 contains the current time point
+                # get the next time 
+                if index1 < N1-1:
+                    t_f1 = t1[index1+1]
+                else:
+                    t_f1 = t_end
+                spike_events[index] = t_p1
                 s1 = dt_p1
-                s2 = (dt_p2*(t2[index2+1]-t1[index1]) + 
-                      dt_f2*(t1[index1]-t2[index2])) / isi2
-                y_ends[index-1] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2)
+                s2 = (dt_p2*(t_f2-t_p1) + dt_f2*(t_p1-t_p2)) / isi2
+                y_ends[index-1] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1,
+                                                                      isi2)
                 # now the next interval start value
-                dt_f1 = get_min_dist_cython(t1[index1+1], t2, N2, index2)
-                isi1 = t1[index1+1]-t1[index1]
+                if index1 < N1-1:
+                    dt_f1 = get_min_dist_cython(t_f1, t2, N2, index2,
+                                                t_start, t_end)
+                    isi1 = t_f1-t_p1
+                else:
+                    dt_f1 = dt_p1
+                    isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
+                    # s1 needs adjustment due to change of isi1
+                    s1 = dt_p1*(t_end-t1[N1-1])/isi1
                 # s2 is the same as above, thus we can compute y2 immediately
-                y_starts[index] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2)
-            elif t1[index1+1] > t2[index2+1]:
+                y_starts[index] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1,
+                                                                      isi2)
+            elif (index2 < N2-1) and (t_f1 > t_f2 or index1 == N1-1):
                 index2 += 1
-                if index2+1 >= N2:
-                    break
-                spike_events[index] = t2[index2]
                 # first calculate the previous interval end value
-                dt_p2 = dt_f2 # the previous time now was the following time before
-                s1 = (dt_p1*(t1[index1+1]-t2[index2]) + 
-                      dt_f1*(t2[index2]-t1[index1])) / isi1
+                # the previous time now was the following time before:
+                dt_p2 = dt_f2
+                t_p2 = t_f2    # t_p1 contains the current time point
+                # get the next time
+                if index2 < N2-1:
+                    t_f2 = t2[index2+1]
+                else:
+                    t_f2 = t_end
+                spike_events[index] = t_p2
+                s1 = (dt_p1*(t_f1-t_p2) + dt_f1*(t_p2-t_p1)) / isi1
                 s2 = dt_p2
-                y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2)
+                y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1,
+                                                                        isi2)
                 # now the next interval start value
-                dt_f2 = get_min_dist_cython(t2[index2+1], t1, N1, index1)
-                #s2 = dt_f2
-                isi2 = t2[index2+1]-t2[index2]
+                if index2 < N2-1:
+                    dt_f2 = get_min_dist_cython(t_f2, t1, N1, index1,
+                                                t_start, t_end)
+                    isi2 = t_f2-t_p2
+                else:
+                    dt_f2 = dt_p2
+                    isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
+                    # s2 needs adjustment due to change of isi2
+                    s2 = dt_p2*(t_end-t2[N2-1])/isi2
                 # s2 is the same as above, thus we can compute y2 immediately
                 y_starts[index] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2)
-            else: # t1[index1+1] == t2[index2+1] - generate only one event
+            else: # t_f1 == t_f2 - generate only one event
                 index1 += 1
                 index2 += 1
-                if (index1+1 >= N1) or (index2+1 >= N2):
-                    break
-                spike_events[index] = t1[index1]
-                y_ends[index-1] = 0.0
-                y_starts[index] = 0.0
+                t_p1 = t_f1
+                t_p2 = t_f2
                 dt_p1 = 0.0
                 dt_p2 = 0.0
-                dt_f1 = get_min_dist_cython(t1[index1+1], t2, N2, index2)
-                dt_f2 = get_min_dist_cython(t2[index2+1], t1, N1, index1)
-                isi1 = t1[index1+1]-t1[index1]
-                isi2 = t2[index2+1]-t2[index2]
+                spike_events[index] = t_f1
+                y_ends[index-1] = 0.0
+                y_starts[index] = 0.0
+                if index1 < N1-1:
+                    t_f1 = t1[index1+1]
+                    dt_f1 = get_min_dist_cython(t_f1, t2, N2, index2,
+                                                t_start, t_end)
+                    isi1 = t_f1 - t_p1
+                else:
+                    t_f1 = t_end
+                    dt_f1 = dt_p1
+                    isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
+                if index2 < N2-1:
+                    t_f2 = t2[index2+1]
+                    dt_f2 = get_min_dist_cython(t_f2, t1, N1, index1,
+                                                t_start, t_end)
+                    isi2 = t_f2 - t_p2
+                else:
+                    t_f2 = t_end
+                    dt_f2 = dt_p2
+                    isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
             index += 1
         # the last event is the interval end
-        spike_events[index] = t1[N1-1]
-        # the ending value of the last interval
-        isi1 = max(t1[N1-1]-t1[N1-2], t1[N1-2]-t1[N1-3])
-        isi2 = max(t2[N2-1]-t2[N2-2], t2[N2-2]-t2[N2-3])
-        s1 = dt_p1*(t1[N1-1]-t1[N1-2])/isi1
-        s2 = dt_p2*(t2[N2-1]-t2[N2-2])/isi2
-        y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2)
+        if spike_events[index-1] == t_end:
+            index -= 1
+        else:
+            spike_events[index] = t_end
+            # the ending value of the last interval
+            isi1 = max(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
+            isi2 = max(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
+            s1 = dt_f1*(t_end-t1[N1-1])/isi1
+            s2 = dt_f2*(t_end-t2[N2-1])/isi2
+            y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2)
     # end nogil
 
     # use only the data added above 
diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py
index 4c37236..bcf9c30 100644
--- a/pyspike/cython/python_backend.py
+++ b/pyspike/cython/python_backend.py
@@ -89,122 +89,185 @@ def isi_distance_python(s1, s2, t_start, t_end):
 ############################################################
 # get_min_dist
 ############################################################
-def get_min_dist(spike_time, spike_train, start_index=0):
+def get_min_dist(spike_time, spike_train, start_index, t_start, t_end):
     """ Returns the minimal distance |spike_time - spike_train[i]|
     with i>=start_index.
     """
-    d = abs(spike_time - spike_train[start_index])
-    start_index += 1
+    d = abs(spike_time - t_start)
+    if start_index < 0:
+        start_index = 0
     while start_index < len(spike_train):
         d_temp = abs(spike_time - spike_train[start_index])
         if d_temp > d:
-            break
+            return d
         else:
             d = d_temp
         start_index += 1
-    return d
+    # finally, check the distance to end time
+    d_temp = abs(t_end - spike_time)
+    if d_temp > d:
+        return d
+    else:
+        return d_temp
 
 
 ############################################################
 # spike_distance_python
 ############################################################
-def spike_distance_python(spikes1, spikes2):
+def spike_distance_python(spikes1, spikes2, t_start, t_end):
     """ Computes the instantaneous spike-distance S_spike (t) of the two given
     spike trains. 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.
     Args:
     - spikes1, spikes2: ordered arrays of spike times with auxiliary spikes.
+    - t_start, t_end: edges of the spike train
     Returns:
     - PieceWiseLinFunc describing the spike-distance.
     """
-    # check for auxiliary spikes - first and last spikes should be identical
-    assert spikes1[0] == spikes2[0], \
-        "Given spike trains seems not to have auxiliary spikes!"
-    assert spikes1[-1] == spikes2[-1], \
-        "Given spike trains seems not to have auxiliary spikes!"
+
     # shorter variables
     t1 = spikes1
     t2 = spikes2
 
-    spike_events = np.empty(len(t1) + len(t2) - 2)
-    spike_events[0] = t1[0]
-    y_starts = np.empty(len(spike_events) - 1)
-    y_ends = np.empty(len(spike_events) - 1)
+    N1 = len(t1)
+    N2 = len(t2)
 
-    index1 = 0
-    index2 = 0
+    spike_events = np.empty(N1+N2+2)
+
+    y_starts = np.empty(len(spike_events)-1)
+    y_ends = np.empty(len(spike_events)-1)
+
+    spike_events[0] = t_start
+    t_p1 = t_start
+    t_p2 = t_start
+    if t1[0] > t_start:
+        # dt_p1 = t2[0]-t_start
+        dt_p1 = 0.0
+        t_f1 = t1[0]
+        dt_f1 = get_min_dist(t_f1, t2, 0, t_start, t_end)
+        isi1 = max(t_f1-t_start, t1[1]-t1[0])
+        s1 = dt_f1*(t_f1-t_start)/isi1
+        index1 = -1
+    else:
+        dt_p1 = 0.0
+        t_f1 = t1[1]
+        dt_f1 = get_min_dist(t_f1, t2, 0, t_start, t_end)
+        isi1 = t1[1]-t1[0]
+        s1 = dt_p1
+        index1 = 0
+    if t2[0] > t_start:
+        # dt_p1 = t2[0]-t_start
+        dt_p2 = 0.0
+        t_f2 = t2[0]
+        dt_f2 = get_min_dist(t_f2, t1, 0, t_start, t_end)
+        isi2 = max(t_f2-t_start, t2[1]-t2[0])
+        s2 = dt_f2*(t_f2-t_start)/isi2
+        index2 = -1
+    else:
+        dt_p2 = 0.0
+        t_f2 = t2[1]
+        dt_f2 = get_min_dist(t_f2, t1, 0, t_start, t_end)
+        isi2 = t2[1]-t2[0]
+        s2 = dt_p2
+        index2 = 0
+
+    y_starts[0] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2)
     index = 1
-    dt_p1 = 0.0
-    dt_f1 = get_min_dist(t1[1], t2, 0)
-    dt_p2 = 0.0
-    dt_f2 = get_min_dist(t2[1], t1, 0)
-    isi1 = max(t1[1]-t1[0], t1[2]-t1[1])
-    isi2 = max(t2[1]-t2[0], t2[2]-t2[1])
-    s1 = dt_f1*(t1[1]-t1[0])/isi1
-    s2 = dt_f2*(t2[1]-t2[0])/isi2
-    y_starts[0] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2)
-    while True:
+
+    while index1+index2 < N1+N2-2:
         # print(index, index1, index2)
-        if t1[index1+1] < t2[index2+1]:
+        if (index1 < N1-1) and (t_f1 < t_f2 or index2 == N2-1):
             index1 += 1
-            # break condition relies on existence of spikes at T_end
-            if index1+1 >= len(t1):
-                break
-            spike_events[index] = t1[index1]
             # first calculate the previous interval end value
-            dt_p1 = dt_f1  # the previous time was the following time before
+            # the previous time now was the following time before:
+            dt_p1 = dt_f1
+            t_p1 = t_f1    # t_p1 contains the current time point
+            # get the next time
+            if index1 < N1-1:
+                t_f1 = t1[index1+1]
+            else:
+                t_f1 = t_end
+            spike_events[index] = t_p1
             s1 = dt_p1
-            s2 = (dt_p2*(t2[index2+1]-t1[index1]) +
-                  dt_f2*(t1[index1]-t2[index2])) / isi2
-            y_ends[index-1] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2)
+            s2 = (dt_p2*(t_f2-t_p1) + dt_f2*(t_p1-t_p2)) / isi2
+            y_ends[index-1] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2)
             # now the next interval start value
-            dt_f1 = get_min_dist(t1[index1+1], t2, index2)
-            isi1 = t1[index1+1]-t1[index1]
+            if index1 < N1-1:
+                dt_f1 = get_min_dist(t_f1, t2, index2, t_start, t_end)
+                isi1 = t_f1-t_p1
+            else:
+                dt_f1 = dt_p1
+                isi1 = max(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
+                # s1 needs adjustment due to change of isi1
+                s1 = dt_p1*(t_end-t1[N1-1])/isi1
             # s2 is the same as above, thus we can compute y2 immediately
-            y_starts[index] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2)
-        elif t1[index1+1] > t2[index2+1]:
+            y_starts[index] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2)
+        elif (index2 < N2-1) and (t_f1 > t_f2 or index1 == N1-1):
             index2 += 1
-            if index2+1 >= len(t2):
-                break
-            spike_events[index] = t2[index2]
             # first calculate the previous interval end value
-            dt_p2 = dt_f2  # the previous time was the following time before
-            s1 = (dt_p1*(t1[index1+1]-t2[index2]) +
-                  dt_f1*(t2[index2]-t1[index1])) / isi1
+            # the previous time now was the following time before:
+            dt_p2 = dt_f2
+            t_p2 = t_f2    # t_p1 contains the current time point
+            # get the next time
+            if index2 < N2-1:
+                t_f2 = t2[index2+1]
+            else:
+                t_f2 = t_end
+            spike_events[index] = t_p2
+            s1 = (dt_p1*(t_f1-t_p2) + dt_f1*(t_p2-t_p1)) / isi1
             s2 = dt_p2
-            y_ends[index-1] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2)
+            y_ends[index-1] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2)
             # now the next interval start value
-            dt_f2 = get_min_dist(t2[index2+1], t1, index1)
-            #s2 = dt_f2
-            isi2 = t2[index2+1]-t2[index2]
+            if index2 < N2-1:
+                dt_f2 = get_min_dist(t_f2, t1, index1, t_start, t_end)
+                isi2 = t_f2-t_p2
+            else:
+                dt_f2 = dt_p2
+                isi2 = max(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
+                # s2 needs adjustment due to change of isi2
+                s2 = dt_p2*(t_end-t2[N2-1])/isi2
             # s2 is the same as above, thus we can compute y2 immediately
-            y_starts[index] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2)
-        else:  # t1[index1+1] == t2[index2+1] - generate only one event
+            y_starts[index] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2)
+        else:  # t_f1 == t_f2 - generate only one event
             index1 += 1
             index2 += 1
-            if (index1+1 >= len(t1)) or (index2+1 >= len(t2)):
-                break
-            assert dt_f2 == 0.0
-            assert dt_f1 == 0.0
-            spike_events[index] = t1[index1]
-            y_ends[index-1] = 0.0
-            y_starts[index] = 0.0
+            t_p1 = t_f1
+            t_p2 = t_f2
             dt_p1 = 0.0
             dt_p2 = 0.0
-            dt_f1 = get_min_dist(t1[index1+1], t2, index2)
-            dt_f2 = get_min_dist(t2[index2+1], t1, index1)
-            isi1 = t1[index1+1]-t1[index1]
-            isi2 = t2[index2+1]-t2[index2]
+            spike_events[index] = t_f1
+            y_ends[index-1] = 0.0
+            y_starts[index] = 0.0
+            if index1 < N1-1:
+                t_f1 = t1[index1+1]
+                dt_f1 = get_min_dist(t_f1, t2, index2, t_start, t_end)
+                isi1 = t_f1 - t_p1
+            else:
+                t_f1 = t_end
+                dt_f1 = dt_p1
+                isi1 = max(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
+            if index2 < N2-1:
+                t_f2 = t2[index2+1]
+                dt_f2 = get_min_dist(t_f2, t1, index1, t_start, t_end)
+                isi2 = t_f2 - t_p2
+            else:
+                t_f2 = t_end
+                dt_f2 = dt_p2
+                isi2 = max(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
         index += 1
     # the last event is the interval end
-    spike_events[index] = t1[-1]
-    # the ending value of the last interval
-    isi1 = max(t1[-1]-t1[-2], t1[-2]-t1[-3])
-    isi2 = max(t2[-1]-t2[-2], t2[-2]-t2[-3])
-    s1 = dt_p1*(t1[-1]-t1[-2])/isi1
-    s2 = dt_p2*(t2[-1]-t2[-2])/isi2
-    y_ends[index-1] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2)
+    if spike_events[index-1] == t_end:
+        index -= 1
+    else:
+        spike_events[index] = t_end
+        # the ending value of the last interval
+        isi1 = max(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
+        isi2 = max(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
+        s1 = dt_f1*(t_end-t1[N1-1])/isi1
+        s2 = dt_f2*(t_end-t2[N2-1])/isi2
+        y_ends[index-1] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2)
+
     # use only the data added above
     # could be less than original length due to equal spike times
     return spike_events[:index+1], y_starts[:index], y_ends[:index]
diff --git a/pyspike/spike_distance.py b/pyspike/spike_distance.py
index f721c86..8d03d70 100644
--- a/pyspike/spike_distance.py
+++ b/pyspike/spike_distance.py
@@ -14,23 +14,23 @@ from pyspike.generic import _generic_profile_multi, _generic_distance_matrix
 ############################################################
 # spike_profile
 ############################################################
-def spike_profile(spikes1, spikes2):
+def spike_profile(spike_train1, spike_train2):
     """ Computes the spike-distance profile S_spike(t) of the two given spike
     trains. Returns the profile as a PieceWiseLinFunc object. The S_spike
-    values are defined positive S_spike(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.
+    values are defined positive S_spike(t)>=0.
 
-    :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: :class:`pyspike.SpikeTrain`
     :returns: The spike-distance profile :math:`S_{spike}(t)`.
     :rtype: :class:`pyspike.function.PieceWiseLinFunc`
 
     """
-    # 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!"
 
     # cython implementation
@@ -45,21 +45,26 @@ Falling back to slow python backend.")
         from cython.python_backend import spike_distance_python \
             as spike_distance_impl
 
-    times, y_starts, y_ends = spike_distance_impl(spikes1, spikes2)
+    times, y_starts, y_ends = spike_distance_impl(spike_train1.spikes,
+                                                  spike_train2.spikes,
+                                                  spike_train1.t_start,
+                                                  spike_train1.t_end)
     return PieceWiseLinFunc(times, y_starts, y_ends)
 
 
 ############################################################
 # spike_distance
 ############################################################
-def spike_distance(spikes1, spikes2, interval=None):
+def spike_distance(spike_train1, spike_train2, interval=None):
     """ Computes the spike-distance S of the given spike trains. The
     spike-distance is the integral over the isi distance profile S_spike(t):
 
     .. math:: S = \int_{T_0}^{T_1} S_{spike}(t) dt.
 
-    :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: :class:`pyspike.SpikeTrain`
     :param interval: averaging interval given as a pair of floats (T0, T1),
                      if None the average over the whole function is computed.
     :type interval: Pair of floats or None.
@@ -67,7 +72,7 @@ def spike_distance(spikes1, spikes2, interval=None):
     :rtype: double
 
     """
-    return spike_profile(spikes1, spikes2).avrg(interval)
+    return spike_profile(spike_train1, spike_train2).avrg(interval)
 
 
 ############################################################
@@ -102,7 +107,7 @@ def spike_distance_multi(spike_trains, indices=None, interval=None):
     S_{spike} = \int_0^T 2/((N(N-1)) sum_{<i,j>} S_{spike}^{i, j} dt
     where the sum goes over all pairs <i,j>
 
-    :param spike_trains: list of spike trains
+    :param spike_trains: list of :class:`pyspike.SpikeTrain`
     :param indices: list of indices defining which spike trains to use,
                     if None all given spike trains are used (default=None)
     :type indices: list or None
@@ -121,7 +126,7 @@ def spike_distance_multi(spike_trains, indices=None, interval=None):
 def spike_distance_matrix(spike_trains, indices=None, interval=None):
     """ Computes the time averaged spike-distance of all pairs of spike-trains.
 
-    :param spike_trains: list of spike trains
+    :param spike_trains: list of :class:`pyspike.SpikeTrain`
     :param indices: list of indices defining which spike trains to use,
                     if None all given spike trains are used (default=None)
     :type indices: list or None
diff --git a/test/test_distance.py b/test/test_distance.py
index b54e908..4af0e63 100644
--- a/test/test_distance.py
+++ b/test/test_distance.py
@@ -64,13 +64,26 @@ def test_isi():
 
 def test_spike():
     # 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.0, 2.0, 5.0, 8.0], 10.0)
+    t2 = SpikeTrain([0.0, 1.0, 5.0, 9.0], 10.0)
+
+    expected_times = np.array([0.0, 1.0, 2.0, 5.0, 8.0, 9.0, 10.0])
+
+    f = spk.spike_profile(t1, t2)
+
+    assert_equal(f.x, expected_times)
+
+    assert_almost_equal(f.avrg(), 0.1662415, decimal=6)
+    assert_almost_equal(f.y2[-1], 0.1394558, decimal=6)
+
+    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 spike 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]
     s1 = np.array([0.1, 0.1, (0.1*0.1+0.05*0.1)/0.2, 0.05, (0.05*0.15 * 2)/0.2,
-                   0.15, 0.1, 0.1*0.2/0.3, 0.1**2/0.3, 0.1*0.05/0.3, 0.1])
+                   0.15, 0.1, (0.1*0.1+0.1*0.2)/0.3, (0.1*0.2+0.1*0.1)/0.3,
+                   (0.1*0.05+0.1*0.25)/0.3, 0.1])
     s2 = np.array([0.1, 0.1*0.2/0.3, 0.1, (0.1*0.05 * 2)/.15, 0.05,
                    (0.05*0.2+0.1*0.15)/0.35, (0.05*0.1+0.1*0.25)/0.35,
                    0.1, 0.1, 0.05, 0.05])
@@ -86,19 +99,18 @@ def test_spike():
                              (expected_y1+expected_y2)/2)
     expected_spike_val /= (expected_times[-1]-expected_times[0])
 
-    t1 = spk.add_auxiliary_spikes(t1, 1.0)
-    t2 = spk.add_auxiliary_spikes(t2, 1.0)
     f = spk.spike_profile(t1, t2)
 
     assert_equal(f.x, expected_times)
     assert_array_almost_equal(f.y1, expected_y1, decimal=15)
     assert_array_almost_equal(f.y2, expected_y2, decimal=15)
-    assert_equal(f.avrg(), expected_spike_val)
-    assert_equal(spk.spike_distance(t1, t2), expected_spike_val)
+    assert_almost_equal(f.avrg(), expected_spike_val, decimal=15)
+    assert_almost_equal(spk.spike_distance(t1, t2), expected_spike_val,
+                        decimal=15)
 
     # 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]
     s1 = np.array([0.1, 0.1*0.1/0.2, 0.1, 0.0, 0.0, 0.0, 0.0])
@@ -115,8 +127,6 @@ def test_spike():
                              (expected_y1+expected_y2)/2)
     expected_spike_val /= (expected_times[-1]-expected_times[0])
 
-    t1 = spk.add_auxiliary_spikes(t1, 1.0)
-    t2 = spk.add_auxiliary_spikes(t2, 1.0)
     f = spk.spike_profile(t1, t2)
 
     assert_equal(f.x, expected_times)
@@ -315,6 +325,6 @@ def test_multi_variate_subsets():
 
 if __name__ == "__main__":
     test_isi()
-    # test_spike()
+    test_spike()
     # test_multi_isi()
     # test_multi_spike()
-- 
cgit v1.2.3


From 36d80c9ec1d28488f9b5c97cd202c196efff694e Mon Sep 17 00:00:00 2001
From: Mario Mulansky <mario.mulansky@gmx.net>
Date: Fri, 24 Apr 2015 15:58:35 +0200
Subject: distance tests now pass with new spike trains

---
 pyspike/cython/cython_distance.pyx |  8 ++++----
 pyspike/cython/python_backend.py   | 41 +++++++++++++++++++-------------------
 test/test_distance.py              | 40 ++++++++++++++++++++++++++++++-------
 3 files changed, 57 insertions(+), 32 deletions(-)

(limited to 'pyspike/cython/python_backend.py')

diff --git a/pyspike/cython/cython_distance.pyx b/pyspike/cython/cython_distance.pyx
index 2841da8..dc2557f 100644
--- a/pyspike/cython/cython_distance.pyx
+++ b/pyspike/cython/cython_distance.pyx
@@ -345,9 +345,9 @@ cdef inline double get_tau(double[:] spikes1, double[:] spikes2,
         m = fmin(m, spikes1[i+1]-spikes1[i])
     if j < N2 and j > -1:
         m = fmin(m, spikes2[j+1]-spikes2[j])
-    if i > 1:
+    if i > 0:
         m = fmin(m, spikes1[i]-spikes1[i-1])
-    if j > 1:
+    if j > 0:
         m = fmin(m, spikes2[j]-spikes2[j-1])
     m *= 0.5
     if max_tau > 0.0:
@@ -371,7 +371,7 @@ def coincidence_cython(double[:] spikes1, double[:] spikes2,
     cdef double[:] mp = np.ones(N1 + N2 + 2)   # multiplicity
     cdef double tau
     while i + j < N1 + N2 - 2:
-        if (i < N1-1) and (spikes1[i+1] < spikes2[j+1] or j == N2-1):
+        if (i < N1-1) and (j == N2-1 or spikes1[i+1] < spikes2[j+1]):
             i += 1
             n += 1
             tau = get_tau(spikes1, spikes2, i, j, max_tau)
@@ -381,7 +381,7 @@ def coincidence_cython(double[:] spikes1, double[:] spikes2,
                 # both get marked with 1
                 c[n] = 1
                 c[n-1] = 1
-        elif (j < N2-1) and (spikes1[i+1] > spikes2[j+1] or i == N1-1):
+        elif (j < N2-1) and (i == N1-1 or spikes1[i+1] > spikes2[j+1]):
             j += 1
             n += 1
             tau = get_tau(spikes1, spikes2, i, j, max_tau)
diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py
index bcf9c30..c65bfb0 100644
--- a/pyspike/cython/python_backend.py
+++ b/pyspike/cython/python_backend.py
@@ -330,47 +330,48 @@ def cumulative_sync_python(spikes1, spikes2):
 ############################################################
 # coincidence_python
 ############################################################
-def coincidence_python(spikes1, spikes2, max_tau):
+def coincidence_python(spikes1, spikes2, t_start, t_end, max_tau):
 
     def get_tau(spikes1, spikes2, i, j, max_tau):
         m = 1E100   # some huge number
-        if i < len(spikes1)-2:
+        if i < len(spikes1)-1 and i > -1:
             m = min(m, spikes1[i+1]-spikes1[i])
-        if j < len(spikes2)-2:
+        if j < len(spikes2)-1 and j > -1:
             m = min(m, spikes2[j+1]-spikes2[j])
-        if i > 1:
+        if i > 0:
             m = min(m, spikes1[i]-spikes1[i-1])
-        if j > 1:
+        if j > 0:
             m = min(m, spikes2[j]-spikes2[j-1])
         m *= 0.5
         if max_tau > 0.0:
             m = min(m, max_tau)
         return m
+
     N1 = len(spikes1)
     N2 = len(spikes2)
-    i = 0
-    j = 0
+    i = -1
+    j = -1
     n = 0
-    st = np.zeros(N1 + N2 - 2)  # spike times
-    c = np.zeros(N1 + N2 - 2)   # coincidences
-    mp = np.ones(N1 + N2 - 2)   # multiplicity
-    while n < N1 + N2 - 2:
-        if spikes1[i+1] < spikes2[j+1]:
+    st = np.zeros(N1 + N2 + 2)  # spike times
+    c = np.zeros(N1 + N2 + 2)   # coincidences
+    mp = np.ones(N1 + N2 + 2)   # multiplicity
+    while i + j < N1 + N2 - 2:
+        if (i < N1-1) and (j == N2-1 or spikes1[i+1] < spikes2[j+1]):
             i += 1
             n += 1
             tau = get_tau(spikes1, spikes2, i, j, max_tau)
             st[n] = spikes1[i]
-            if j > 0 and spikes1[i]-spikes2[j] < tau:
+            if j > -1 and spikes1[i]-spikes2[j] < tau:
                 # coincidence between the current spike and the previous spike
                 # both get marked with 1
                 c[n] = 1
                 c[n-1] = 1
-        elif spikes1[i+1] > spikes2[j+1]:
+        elif (j < N2-1) and (i == N1-1 or spikes1[i+1] > spikes2[j+1]):
             j += 1
             n += 1
             tau = get_tau(spikes1, spikes2, i, j, max_tau)
             st[n] = spikes2[j]
-            if i > 0 and spikes2[j]-spikes1[i] < tau:
+            if i > -1 and spikes2[j]-spikes1[i] < tau:
                 # coincidence between the current spike and the previous spike
                 # both get marked with 1
                 c[n] = 1
@@ -379,8 +380,6 @@ def coincidence_python(spikes1, spikes2, max_tau):
             # advance in both spike trains
             j += 1
             i += 1
-            if i == N1-1 or j == N2-1:
-                break
             n += 1
             # add only one event, but with coincidence 2 and multiplicity 2
             st[n] = spikes1[i]
@@ -391,12 +390,12 @@ def coincidence_python(spikes1, spikes2, max_tau):
     c = c[:n+2]
     mp = mp[:n+2]
 
-    st[0] = spikes1[0]
-    st[-1] = spikes1[-1]
+    st[0] = t_start
+    st[len(st)-1] = t_end
     c[0] = c[1]
-    c[-1] = c[-2]
+    c[len(c)-1] = c[len(c)-2]
     mp[0] = mp[1]
-    mp[-1] = mp[-2]
+    mp[len(mp)-1] = mp[len(mp)-2]
 
     return st, c, mp
 
diff --git a/test/test_distance.py b/test/test_distance.py
index 0fff840..88cf40e 100644
--- a/test/test_distance.py
+++ b/test/test_distance.py
@@ -177,6 +177,18 @@ def test_spike_sync():
     assert_almost_equal(spk.spike_sync(spikes1, spikes2),
                         0.5, decimal=16)
 
+    spikes2 = SpikeTrain([3.0], 4.0)
+    assert_almost_equal(spk.spike_sync(spikes1, spikes2),
+                        0.5, decimal=16)
+
+    spikes2 = SpikeTrain([1.0], 4.0)
+    assert_almost_equal(spk.spike_sync(spikes1, spikes2),
+                        0.5, decimal=16)
+
+    spikes2 = SpikeTrain([1.5, 3.0], 4.0)
+    assert_almost_equal(spk.spike_sync(spikes1, spikes2),
+                        0.4, decimal=16)
+
 
 def check_multi_profile(profile_func, profile_func_multi):
     # generate spike trains:
@@ -250,19 +262,28 @@ def test_multi_spike_sync():
 
     # multivariate regression test
     spike_trains = spk.load_spike_trains_from_txt("test/SPIKE_Sync_Test.txt",
-                                                  interval=(0, 4000))
-    print(spike_trains[0].spikes)
+                                                  interval=[0, 4000])
+    # extract all spike times
+    spike_times = np.array([])
+    for st in spike_trains:
+        spike_times = np.append(spike_times, st.spikes)
+    spike_times = np.unique(np.sort(spike_times))
+
     f = spk.spike_sync_profile_multi(spike_trains)
+
+    assert_equal(spike_times, f.x[1:-1])
+    assert_equal(len(f.x), len(f.y))
+
     assert_equal(np.sum(f.y[1:-1]), 39932)
     assert_equal(np.sum(f.mp[1:-1]), 85554)
 
 
 def check_dist_matrix(dist_func, dist_matrix_func):
     # generate spike trains:
-    t1 = spk.add_auxiliary_spikes(np.array([0.2, 0.4, 0.6, 0.7]), 1.0)
-    t2 = spk.add_auxiliary_spikes(np.array([0.3, 0.45, 0.8, 0.9, 0.95]), 1.0)
-    t3 = spk.add_auxiliary_spikes(np.array([0.2, 0.4, 0.6]), 1.0)
-    t4 = spk.add_auxiliary_spikes(np.array([0.1, 0.4, 0.5, 0.6]), 1.0)
+    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)
+    t3 = SpikeTrain([0.2, 0.4, 0.6], 1.0)
+    t4 = SpikeTrain([0.1, 0.4, 0.5, 0.6], 1.0)
     spike_trains = [t1, t2, t3, t4]
 
     f12 = dist_func(t1, t2)
@@ -340,4 +361,9 @@ if __name__ == "__main__":
     test_spike_sync()
     test_multi_isi()
     test_multi_spike()
-    # test_multi_spike_sync()
+    test_multi_spike_sync()
+    test_isi_matrix()
+    test_spike_matrix()
+    test_spike_sync_matrix()
+    test_regression_spiky()
+    test_multi_variate_subsets()
-- 
cgit v1.2.3


From 795e16ffe7afb469ef07a548c1f6a31d924196b3 Mon Sep 17 00:00:00 2001
From: Mario Mulansky <mario.mulansky@gmx.net>
Date: Fri, 24 Apr 2015 23:29:05 +0200
Subject: bugfixes for spike distance

---
 pyspike/cython/cython_distance.pyx | 24 +++++++++---------
 pyspike/cython/python_backend.py   | 15 ++++++------
 test/test_distance.py              | 50 ++++++++++++++++++++++++++++----------
 3 files changed, 58 insertions(+), 31 deletions(-)

(limited to 'pyspike/cython/python_backend.py')

diff --git a/pyspike/cython/cython_distance.pyx b/pyspike/cython/cython_distance.pyx
index dc2557f..a41d8e8 100644
--- a/pyspike/cython/cython_distance.pyx
+++ b/pyspike/cython/cython_distance.pyx
@@ -194,31 +194,31 @@ def spike_distance_cython(double[:] t1, double[:] t2,
         t_p2 = t_start
         if t1[0] > t_start:
             # dt_p1 = t2[0]-t_start
-            dt_p1 = 0.0
             t_f1 = t1[0]
             dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_start, t_end)
             isi1 = fmax(t_f1-t_start, t1[1]-t1[0])
-            s1 = dt_f1*(t_f1-t_start)/isi1
+            dt_p1 = dt_f1
+            s1 = dt_p1*(t_f1-t_start)/isi1
             index1 = -1
         else:
-            dt_p1 = 0.0
             t_f1 = t1[1]
             dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_start, t_end)
+            dt_p1 = 0.0
             isi1 = t1[1]-t1[0]
             s1 = dt_p1
             index1 = 0
         if t2[0] > t_start:
             # dt_p1 = t2[0]-t_start
-            dt_p2 = 0.0
             t_f2 = t2[0]
             dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_start, t_end)
+            dt_p2 = dt_f2
             isi2 = fmax(t_f2-t_start, t2[1]-t2[0])
-            s2 = dt_f2*(t_f2-t_start)/isi2
+            s2 = dt_p2*(t_f2-t_start)/isi2
             index2 = -1
         else:
-            dt_p2 = 0.0
             t_f2 = t2[1]
             dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_start, t_end)
+            dt_p2 = 0.0
             isi2 = t2[1]-t2[0]
             s2 = dt_p2
             index2 = 0
@@ -231,16 +231,16 @@ def spike_distance_cython(double[:] t1, double[:] t2,
             if (index1 < N1-1) and (t_f1 < t_f2 or index2 == N2-1):
                 index1 += 1
                 # first calculate the previous interval end value
+                s1 = dt_f1*(t_f1-t_p1) / isi1
                 # the previous time now was the following time before:
-                dt_p1 = dt_f1  
+                dt_p1 = dt_f1
                 t_p1 = t_f1    # t_p1 contains the current time point
-                # get the next time 
+                # get the next time
                 if index1 < N1-1:
                     t_f1 = t1[index1+1]
                 else:
                     t_f1 = t_end
                 spike_events[index] = t_p1
-                s1 = dt_p1
                 s2 = (dt_p2*(t_f2-t_p1) + dt_f2*(t_p1-t_p2)) / isi2
                 y_ends[index-1] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1,
                                                                       isi2)
@@ -249,6 +249,7 @@ def spike_distance_cython(double[:] t1, double[:] t2,
                     dt_f1 = get_min_dist_cython(t_f1, t2, N2, index2,
                                                 t_start, t_end)
                     isi1 = t_f1-t_p1
+                    s1 = dt_p1
                 else:
                     dt_f1 = dt_p1
                     isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
@@ -260,9 +261,10 @@ def spike_distance_cython(double[:] t1, double[:] t2,
             elif (index2 < N2-1) and (t_f1 > t_f2 or index1 == N1-1):
                 index2 += 1
                 # first calculate the previous interval end value
+                s2 = dt_f2*(t_f2-t_p2) / isi2
                 # the previous time now was the following time before:
                 dt_p2 = dt_f2
-                t_p2 = t_f2    # t_p1 contains the current time point
+                t_p2 = t_f2    # t_p2 contains the current time point
                 # get the next time
                 if index2 < N2-1:
                     t_f2 = t2[index2+1]
@@ -270,7 +272,6 @@ def spike_distance_cython(double[:] t1, double[:] t2,
                     t_f2 = t_end
                 spike_events[index] = t_p2
                 s1 = (dt_p1*(t_f1-t_p2) + dt_f1*(t_p2-t_p1)) / isi1
-                s2 = dt_p2
                 y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1,
                                                                         isi2)
                 # now the next interval start value
@@ -278,6 +279,7 @@ def spike_distance_cython(double[:] t1, double[:] t2,
                     dt_f2 = get_min_dist_cython(t_f2, t1, N1, index1,
                                                 t_start, t_end)
                     isi2 = t_f2-t_p2
+                    s2 = dt_p2
                 else:
                     dt_f2 = dt_p2
                     isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py
index c65bfb0..317b568 100644
--- a/pyspike/cython/python_backend.py
+++ b/pyspike/cython/python_backend.py
@@ -142,12 +142,11 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
     t_p1 = t_start
     t_p2 = t_start
     if t1[0] > t_start:
-        # dt_p1 = t2[0]-t_start
-        dt_p1 = 0.0
         t_f1 = t1[0]
         dt_f1 = get_min_dist(t_f1, t2, 0, t_start, t_end)
+        dt_p1 = dt_f1
         isi1 = max(t_f1-t_start, t1[1]-t1[0])
-        s1 = dt_f1*(t_f1-t_start)/isi1
+        s1 = dt_p1*(t_f1-t_start)/isi1
         index1 = -1
     else:
         dt_p1 = 0.0
@@ -158,11 +157,11 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
         index1 = 0
     if t2[0] > t_start:
         # dt_p1 = t2[0]-t_start
-        dt_p2 = 0.0
         t_f2 = t2[0]
         dt_f2 = get_min_dist(t_f2, t1, 0, t_start, t_end)
+        dt_p2 = dt_f2
         isi2 = max(t_f2-t_start, t2[1]-t2[0])
-        s2 = dt_f2*(t_f2-t_start)/isi2
+        s2 = dt_p2*(t_f2-t_start)/isi2
         index2 = -1
     else:
         dt_p2 = 0.0
@@ -180,6 +179,7 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
         if (index1 < N1-1) and (t_f1 < t_f2 or index2 == N2-1):
             index1 += 1
             # first calculate the previous interval end value
+            s1 = dt_f1*(t_f1-t_p1) / isi1
             # the previous time now was the following time before:
             dt_p1 = dt_f1
             t_p1 = t_f1    # t_p1 contains the current time point
@@ -189,13 +189,13 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
             else:
                 t_f1 = t_end
             spike_events[index] = t_p1
-            s1 = dt_p1
             s2 = (dt_p2*(t_f2-t_p1) + dt_f2*(t_p1-t_p2)) / isi2
             y_ends[index-1] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2)
             # now the next interval start value
             if index1 < N1-1:
                 dt_f1 = get_min_dist(t_f1, t2, index2, t_start, t_end)
                 isi1 = t_f1-t_p1
+                s1 = dt_p1
             else:
                 dt_f1 = dt_p1
                 isi1 = max(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
@@ -206,6 +206,7 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
         elif (index2 < N2-1) and (t_f1 > t_f2 or index1 == N1-1):
             index2 += 1
             # first calculate the previous interval end value
+            s2 = dt_f2*(t_f2-t_p2) / isi2
             # the previous time now was the following time before:
             dt_p2 = dt_f2
             t_p2 = t_f2    # t_p1 contains the current time point
@@ -216,12 +217,12 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
                 t_f2 = t_end
             spike_events[index] = t_p2
             s1 = (dt_p1*(t_f1-t_p2) + dt_f1*(t_p2-t_p1)) / isi1
-            s2 = dt_p2
             y_ends[index-1] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2)
             # now the next interval start value
             if index2 < N2-1:
                 dt_f2 = get_min_dist(t_f2, t1, index1, t_start, t_end)
                 isi2 = t_f2-t_p2
+                s2 = dt_p2
             else:
                 dt_f2 = dt_p2
                 isi2 = max(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
diff --git a/test/test_distance.py b/test/test_distance.py
index 0059001..20b52e8 100644
--- a/test/test_distance.py
+++ b/test/test_distance.py
@@ -73,7 +73,7 @@ def test_spike():
 
     assert_equal(f.x, expected_times)
 
-    assert_almost_equal(f.avrg(), 0.1662415, decimal=6)
+    assert_almost_equal(f.avrg(), 1.6624149659863946e-01, decimal=15)
     assert_almost_equal(f.y2[-1], 0.1394558, decimal=6)
 
     t1 = SpikeTrain([0.2, 0.4, 0.6, 0.7], 1.0)
@@ -84,7 +84,7 @@ def test_spike():
     s1 = np.array([0.1, 0.1, (0.1*0.1+0.05*0.1)/0.2, 0.05, (0.05*0.15 * 2)/0.2,
                    0.15, 0.1, (0.1*0.1+0.1*0.2)/0.3, (0.1*0.2+0.1*0.1)/0.3,
                    (0.1*0.05+0.1*0.25)/0.3, 0.1])
-    s2 = np.array([0.1, 0.1*0.2/0.3, 0.1, (0.1*0.05 * 2)/.15, 0.05,
+    s2 = np.array([0.1, (0.1*0.2+0.1*0.1)/0.3, 0.1, (0.1*0.05 * 2)/.15, 0.05,
                    (0.05*0.2+0.1*0.15)/0.35, (0.05*0.1+0.1*0.25)/0.35,
                    0.1, 0.1, 0.05, 0.05])
     isi1 = np.array([0.2, 0.2, 0.2, 0.2, 0.2, 0.1, 0.3, 0.3, 0.3, 0.3])
@@ -113,12 +113,18 @@ def test_spike():
     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]
-    s1 = np.array([0.1, 0.1*0.1/0.2, 0.1, 0.0, 0.0, 0.0, 0.0])
-    s2 = np.array([0.1*0.1/0.3, 0.1, 0.1*0.2/0.3, 0.0, 0.1, 0.0, 0.0])
+    # due to the edge correction in the beginning, s1 and s2 are different
+    # for left and right values
+    s1_r = np.array([0.1, (0.1*0.1+0.1*0.1)/0.2, 0.1, 0.0, 0.0, 0.0, 0.0])
+    s1_l = np.array([0.1, (0.1*0.1+0.1*0.1)/0.2, 0.1, 0.0, 0.0, 0.0, 0.0])
+    s2_r = np.array([0.1*0.1/0.3, 0.1*0.3/0.3, 0.1*0.2/0.3,
+                     0.0, 0.1, 0.0, 0.0])
+    s2_l = np.array([0.1*0.1/0.3, 0.1*0.1/0.3, 0.1*0.2/0.3, 0.0,
+                     0.1, 0.0, 0.0])
     isi1 = np.array([0.2, 0.2, 0.2, 0.2, 0.2, 0.4])
     isi2 = np.array([0.3, 0.3, 0.3, 0.1, 0.1, 0.4])
-    expected_y1 = (s1[:-1]*isi2+s2[:-1]*isi1) / (0.5*(isi1+isi2)**2)
-    expected_y2 = (s1[1:]*isi2+s2[1:]*isi1) / (0.5*(isi1+isi2)**2)
+    expected_y1 = (s1_r[:-1]*isi2+s2_r[:-1]*isi1) / (0.5*(isi1+isi2)**2)
+    expected_y2 = (s1_l[1:]*isi2+s2_l[1:]*isi1) / (0.5*(isi1+isi2)**2)
 
     expected_times = np.array(expected_times)
     expected_y1 = np.array(expected_y1)
@@ -321,19 +327,37 @@ def test_spike_sync_matrix():
 
 
 def test_regression_spiky():
+    # standard example
+    st1 = SpikeTrain(np.arange(100, 1201, 100), 1300)
+    st2 = SpikeTrain(np.arange(100, 1201, 110), 1300)
+
+    isi_dist = spk.isi_distance(st1, st2)
+    assert_almost_equal(isi_dist, 7.6923076923076941e-02, decimal=15)
+
+    spike_dist = spk.spike_distance(st1, st2)
+    assert_equal(spike_dist, 2.1105878248735391e-01)
+
+    spike_sync = spk.spike_sync(st1, st2)
+    assert_equal(spike_sync, 8.6956521739130432e-01)
+
+    # multivariate check
+
     spike_trains = spk.load_spike_trains_from_txt("test/PySpike_testdata.txt",
                                                   (0.0, 4000.0))
-    isi_profile = spk.isi_profile_multi(spike_trains)
-    isi_dist = isi_profile.avrg()
-    print(isi_dist)
+    isi_dist = spk.isi_distance_multi(spike_trains)
     # get the full precision from SPIKY
-    # assert_equal(isi_dist, 0.1832)
+    assert_almost_equal(isi_dist, 1.8318789829845508e-01, decimal=15)
 
     spike_profile = spk.spike_profile_multi(spike_trains)
-    spike_dist = spike_profile.avrg()
-    print(spike_dist)
+    assert_equal(len(spike_profile.y1)+len(spike_profile.y2), 1252)
+
+    spike_dist = spk.spike_distance_multi(spike_trains)
+    # get the full precision from SPIKY
+    assert_almost_equal(spike_dist, 2.4432433330596512e-01, decimal=15)
+
+    spike_sync = spk.spike_sync_multi(spike_trains)
     # get the full precision from SPIKY
-    # assert_equal(spike_dist, 0.2445)
+    assert_equal(spike_sync, 0.7183531505298066)
 
 
 def test_multi_variate_subsets():
-- 
cgit v1.2.3


From 2e7351393927ba9e9e0c3b7b59d05e8aeeb41d1f Mon Sep 17 00:00:00 2001
From: Mario Mulansky <mario.mulansky@gmx.net>
Date: Tue, 28 Apr 2015 16:11:11 +0200
Subject: edge correction for the ISI-distance

---
 pyspike/cython/cython_distance.pyx | 18 ++++++++++++------
 pyspike/cython/python_backend.py   | 18 ++++++++++++------
 test/test_distance.py              |  8 +++++---
 3 files changed, 29 insertions(+), 15 deletions(-)

(limited to 'pyspike/cython/python_backend.py')

diff --git a/pyspike/cython/cython_distance.pyx b/pyspike/cython/cython_distance.pyx
index a41d8e8..6ee0181 100644
--- a/pyspike/cython/cython_distance.pyx
+++ b/pyspike/cython/cython_distance.pyx
@@ -62,14 +62,16 @@ def isi_distance_cython(double[:] s1, double[:] s2,
 
     # first interspike interval - check if a spike exists at the start time
     if s1[0] > t_start:
-        nu1 = s1[0] - t_start
+        # edge correction
+        nu1 = fmax(s1[0]-t_start, s1[1]-s1[0])
         index1 = -1
     else:
         nu1 = s1[1]-s1[0]
         index1 = 0
 
     if s2[0] > t_start:
-        nu2 = s2[0] - t_start
+        # edge correction
+        nu2 = fmax(s2[0]-t_start, s2[1]-s2[0])
         index2 = -1
     else:
         nu2 = s2[1]-s2[0]
@@ -89,7 +91,8 @@ def isi_distance_cython(double[:] s1, double[:] s2,
                 if index1 < N1-1:
                     nu1 = s1[index1+1]-s1[index1]
                 else:
-                    nu1 = t_end-s1[index1]
+                    # edge correction
+                    nu1 = fmax(t_end-s1[index1], nu1)
             elif (index2 < N2-1) and ((index1 == N1-1) or
                                       (s1[index1+1] > s2[index2+1])):
                 index2 += 1
@@ -97,7 +100,8 @@ def isi_distance_cython(double[:] s1, double[:] s2,
                 if index2 < N2-1:
                     nu2 = s2[index2+1]-s2[index2]
                 else:
-                    nu2 = t_end-s2[index2]
+                    # edge correction
+                    nu2 = fmax(t_end-s2[index2], nu2)
             else: # s1[index1+1] == s2[index2+1]
                 index1 += 1
                 index2 += 1
@@ -105,11 +109,13 @@ def isi_distance_cython(double[:] s1, double[:] s2,
                 if index1 < N1-1:
                     nu1 = s1[index1+1]-s1[index1]
                 else:
-                    nu1 = t_end-s1[index1]
+                    # edge correction
+                    nu1 = fmax(t_end-s1[index1], nu1)
                 if index2 < N2-1:
                     nu2 = s2[index2+1]-s2[index2]
                 else:
-                    nu2 = t_end-s2[index2]
+                    # edge correction
+                    nu2 = fmax(t_end-s2[index2], nu2)
             # compute the corresponding isi-distance
             isi_values[index] = fabs(nu1 - nu2) / fmax(nu1, nu2)
             index += 1
diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py
index 317b568..1fd8c42 100644
--- a/pyspike/cython/python_backend.py
+++ b/pyspike/cython/python_backend.py
@@ -27,13 +27,15 @@ def isi_distance_python(s1, s2, t_start, t_end):
     # the values have one entry less - the number of intervals between events
     isi_values = np.empty(len(spike_events) - 1)
     if s1[0] > t_start:
-        nu1 = s1[0] - t_start
+        # edge correction
+        nu1 = max(s1[0] - t_start, s1[1] - s1[0])
         index1 = -1
     else:
         nu1 = s1[1] - s1[0]
         index1 = 0
     if s2[0] > t_start:
-        nu2 = s2[0] - t_start
+        # edge correction
+        nu2 = max(s2[0] - t_start, s2[1] - s2[0])
         index2 = -1
     else:
         nu2 = s2[1] - s2[0]
@@ -49,7 +51,8 @@ def isi_distance_python(s1, s2, t_start, t_end):
             if index1 < N1-1:
                 nu1 = s1[index1+1]-s1[index1]
             else:
-                nu1 = t_end-s1[index1]
+                # edge correction
+                nu1 = max(t_end-s1[N1-1], s1[N1-1]-s1[N1-2])
 
         elif (index2 < N2-1) and (index1 == N1-1 or
                                   s1[index1+1] > s2[index2+1]):
@@ -58,7 +61,8 @@ def isi_distance_python(s1, s2, t_start, t_end):
             if index2 < N2-1:
                 nu2 = s2[index2+1]-s2[index2]
             else:
-                nu2 = t_end-s2[index2]
+                # edge correction
+                nu2 = max(t_end-s2[N2-1], s2[N2-1]-s2[N2-2])
 
         else:  # s1[index1 + 1] == s2[index2 + 1]
             index1 += 1
@@ -67,11 +71,13 @@ def isi_distance_python(s1, s2, t_start, t_end):
             if index1 < N1-1:
                 nu1 = s1[index1+1]-s1[index1]
             else:
-                nu1 = t_end-s1[index1]
+                # edge correction
+                nu1 = max(t_end-s1[N1-1], s1[N1-1]-s1[N1-2])
             if index2 < N2-1:
                 nu2 = s2[index2+1]-s2[index2]
             else:
-                nu2 = t_end-s2[index2]
+                # edge correction
+                nu2 = max(t_end-s2[N2-1], s2[N2-1]-s2[N2-2])
         # compute the corresponding isi-distance
         isi_values[index] = abs(nu1 - nu2) / \
             max(nu1, nu2)
diff --git a/test/test_distance.py b/test/test_distance.py
index 20b52e8..19da35f 100644
--- a/test/test_distance.py
+++ b/test/test_distance.py
@@ -47,7 +47,7 @@ def test_isi():
     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]
+    expected_isi = [0.1/0.3, 0.1/0.3, 0.1/0.3, 0.1/0.2, 0.1/0.2, 0.0/0.5]
     expected_times = np.array(expected_times)
     expected_isi = np.array(expected_isi)
 
@@ -332,7 +332,9 @@ def test_regression_spiky():
     st2 = SpikeTrain(np.arange(100, 1201, 110), 1300)
 
     isi_dist = spk.isi_distance(st1, st2)
-    assert_almost_equal(isi_dist, 7.6923076923076941e-02, decimal=15)
+    assert_almost_equal(isi_dist, 9.0909090909090939e-02, decimal=15)
+    isi_profile = spk.isi_profile(st1, st2)
+    assert_equal(isi_profile.y, 0.1/1.1 * np.ones_like(isi_profile.y))
 
     spike_dist = spk.spike_distance(st1, st2)
     assert_equal(spike_dist, 2.1105878248735391e-01)
@@ -346,7 +348,7 @@ def test_regression_spiky():
                                                   (0.0, 4000.0))
     isi_dist = spk.isi_distance_multi(spike_trains)
     # get the full precision from SPIKY
-    assert_almost_equal(isi_dist, 1.8318789829845508e-01, decimal=15)
+    assert_almost_equal(isi_dist, 0.17051816816999129656, decimal=15)
 
     spike_profile = spk.spike_profile_multi(spike_trains)
     assert_equal(len(spike_profile.y1)+len(spike_profile.y2), 1252)
-- 
cgit v1.2.3