edge correction for the ISI-distance

2 files changed, 24 insertions, 12 deletions

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)