bugfix for edge correction

fixed bug within new edge correction (auxiliary spike was ignored in some
cases)
added regression test with 10000 random spike train sets

2 files changed, 26 insertions, 5 deletions

diff --git a/pyspike/cython/cython_distances.pyx b/pyspike/cython/cython_distances.pyx
index 41baa60..6c6e7e5 100644
--- a/pyspike/cython/cython_distances.pyx
+++ b/pyspike/cython/cython_distances.pyx
@@ -184,7 +184,8 @@ def spike_distance_cython(double[:] t1, double[:] t2,
     N1 = len(t1)
     N2 = len(t2)
 
-    with nogil: # release the interpreter to allow multithreading
+    # with nogil: # release the interpreter to allow multithreading
+    if True:
         t_last = t_start
         # t_p1 = t_start
         # t_p2 = t_start
@@ -193,6 +194,7 @@ def spike_distance_cython(double[:] t1, double[:] t2,
         t_aux1[1] = fmax(t_end, t1[N1-1]+(t1[N1-1]-t1[N1-2]))
         t_aux2[0] = fmin(t_start, t2[0]-(t2[1]-t2[0]))
         t_aux2[1] = fmax(t_end, t2[N2-1]+(t2[N2-1]-t2[N2-2]))
+        # print "aux spikes %.15f, %.15f ; %.15f, %.15f" % (t_aux1[0], t_aux1[1], t_aux2[0], t_aux2[1])
         t_p1 = t_start if (t1[0] == t_start) else t_aux1[0]
         t_p2 = t_start if (t2[0] == t_start) else t_aux2[0]
         if t1[0] > t_start:
@@ -207,7 +209,8 @@ def spike_distance_cython(double[:] t1, double[:] t2,
         else:
             t_f1 = t1[1]
             dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_aux2[0], t_aux2[1])
-            dt_p1 = 0.0
+            # dt_p1 = t_start-t_p2  # 0.0
+            dt_p1 = get_min_dist_cython(t_p1, t2, N2, 0, t_aux2[0], t_aux2[1])
             isi1 = t1[1]-t1[0]
             s1 = dt_p1
             index1 = 0
@@ -223,7 +226,8 @@ def spike_distance_cython(double[:] t1, double[:] t2,
         else:
             t_f2 = t2[1]
             dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_aux1[0], t_aux1[1])
-            dt_p2 = 0.0
+            # dt_p2 = t_start-t_p1  # 0.0
+            dt_p2 = get_min_dist_cython(t_p2, t1, N1, 0, t_aux1[0], t_aux1[1])
             isi2 = t2[1]-t2[0]
             s2 = dt_p2
             index2 = 0
diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py
index a5f7ae4..a007a7f 100644
--- a/pyspike/cython/python_backend.py
+++ b/pyspike/cython/python_backend.py
@@ -153,6 +153,8 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
     t_p1 = t_start if (t1[0] == t_start) else t_aux1[0]
     t_p2 = t_start if (t2[0] == t_start) else t_aux2[0]
 
+    print "t_aux1", t_aux1, ", t_aux2:", t_aux2
+
     spike_events[0] = t_start
     if t1[0] > t_start:
         t_f1 = t1[0]
@@ -163,7 +165,8 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
         s1 = dt_p1
         index1 = -1
     else:
-        dt_p1 = 0.0
+        # dt_p1 = t_start-t_p2
+        dt_p1 = get_min_dist(t_p1, t2, 0, t_aux2[0], t_aux2[1])
         t_f1 = t1[1]
         dt_f1 = get_min_dist(t_f1, t2, 0, t_aux2[0], t_aux2[1])
         isi1 = t1[1]-t1[0]
@@ -179,13 +182,18 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
         s2 = dt_p2
         index2 = -1
     else:
-        dt_p2 = 0.0
+        dt_p2 = t_start-t_p1
+        dt_p2 = get_min_dist(t_p2, t1, 0, t_aux1[0], t_aux1[1])
         t_f2 = t2[1]
         dt_f2 = get_min_dist(t_f2, t1, 0, t_aux1[0], t_aux1[1])
         isi2 = t2[1]-t2[0]
         s2 = dt_p2
         index2 = 0
 
+    print "t_p1:", repr(t_p1), ", t_f1:", repr(t_f1), ", dt_p1:", repr(dt_p1), ", dt_f1:", repr(dt_f1)
+    print "t_p2:", repr(t_p2), ", t_f2:", repr(t_f2), ", dt_p2:", repr(dt_p2), ", dt_f2:", repr(dt_f2)
+    print "s1: ", repr(s1), ", s2:", repr(s2)
+
     y_starts[0] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2)
     index = 1
 
@@ -276,6 +284,11 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
                 dt_f2 = dt_p2
                 isi2 = max(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
         index += 1
+
+        print "t_p1:", repr(t_p1), ", t_f1:", repr(t_f1), ", dt_p1:", repr(dt_p1), ", dt_f1:", repr(dt_f1)
+        print "t_p2:", repr(t_p2), ", t_f2:", repr(t_f2), ", dt_p2:", repr(dt_p2), ", dt_f2:", repr(dt_f2)
+        print "s1: ", repr(s1), ", s2:", repr(s2)
+
     # the last event is the interval end
     if spike_events[index-1] == t_end:
         index -= 1
@@ -288,6 +301,10 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
         s2 = dt_f2  # *(t_end-t2[N2-1])/isi2
         y_ends[index-1] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2)
 
+    print "t_p1:", repr(t_p1), ", t_f1:", repr(t_f1), ", dt_p1:", repr(dt_p1), ", dt_f1:", repr(dt_f1)
+    print "t_p2:", repr(t_p2), ", t_f2:", repr(t_f2), ", dt_p2:", repr(dt_p2), ", dt_f2:", repr(dt_f2)
+    print "s1: ", repr(s1), ", s2:", repr(s2)
+
     # 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]