improved implementation of spike_sync

2 files changed, 72 insertions, 16 deletions

diff --git a/pyspike/distances.py b/pyspike/distances.py
index fbedce5..c28fd7a 100644
--- a/pyspike/distances.py
+++ b/pyspike/distances.py
@@ -138,30 +138,32 @@ def spike_sync_profile(spikes1, spikes2, k=3):
 
     # cython implementation
     try:
-        from cython_distance import cumulative_sync_cython \
-            as cumulative_sync_impl
+        from cython_distance import coincidence_cython \
+            as coincidence_impl
     except ImportError:
 #         print("Warning: spike_distance_cython not found. Make sure that \
 # PySpike is installed by running\n 'python setup.py build_ext --inplace'!\n \
 # Falling back to slow python backend.")
         # use python backend
-        from python_backend import cumulative_sync_python \
-            as cumulative_sync_impl
+        from python_backend import coincidence_python \
+            as coincidence_impl
 
-    st, c = cumulative_sync_impl(spikes1, spikes2)
-
-    # print c
-    # print 2*(c[-1]-c[0])/(len(spikes1)+len(spikes2)-2)
+    st, c = coincidence_impl(spikes1, spikes2)
 
     dc = np.zeros(len(c))
-    dc[k:-k] = (c[2*k:] - c[:-2*k]) / k
-    for n in xrange(1, k):
-        dc[n] = (c[2*n] - c[0]) / k
-        dc[-n-1] = (c[-1]-c[-2*n-1]) / k
-    dc[0] = dc[1]
-    dc[-1] = dc[-2]
-    # dc[-1] = (c[-1]-c[-2])/k
-    # print dc
+    for i in xrange(2*k):
+        dc[k:-k] += c[i:-2*k+i]
+
+    for n in xrange(0, k):
+        for i in xrange(n+k):
+            dc[n] += c[i]
+            dc[-n-1] += c[-i-1]
+        for i in xrange(k-n-1):
+            dc[n] += c[i]
+            dc[-n-1] += c[-i-1]
+
+    dc *= 1.0/k
+
     return PieceWiseConstFunc(st, dc)
 
 
diff --git a/pyspike/python_backend.py b/pyspike/python_backend.py
index b85262d..7f8ea8c 100644
--- a/pyspike/python_backend.py
+++ b/pyspike/python_backend.py
@@ -243,6 +243,60 @@ def cumulative_sync_python(spikes1, spikes2):
 
 
 ############################################################
+# coincidence_python
+############################################################
+def coincidence_python(spikes1, spikes2):
+
+    def get_tau(spikes1, spikes2, i, j):
+        return 0.5*min([spikes1[i]-spikes1[i-1], spikes1[i+1]-spikes1[i],
+                        spikes2[j]-spikes2[j-1], spikes2[j+1]-spikes2[j]])
+    N1 = len(spikes1)
+    N2 = len(spikes2)
+    i = 0
+    j = 0
+    n = 0
+    st = np.zeros(N1 + N2 - 2)
+    c = np.zeros(N1 + N2 - 3)
+    c[0] = 0
+    st[0] = 0
+    while n < N1 + N2:
+        if spikes1[i+1] < spikes2[j+1]:
+            i += 1
+            n += 1
+            tau = get_tau(spikes1, spikes2, i, j)
+            st[n] = spikes1[i]
+            if spikes1[i]-spikes2[j] > tau:
+                c[n] = 0
+            else:
+                c[n] = 1
+        elif spikes1[i+1] > spikes2[j+1]:
+            j += 1
+            n += 1
+            tau = get_tau(spikes1, spikes2, i, j)
+            st[n] = spikes2[j]
+            if spikes2[j]-spikes1[i] > tau:
+                c[n] = 0
+            else:
+                c[n] = 1
+        else:   # spikes1[i+1] = spikes2[j+1]
+            j += 1
+            i += 1
+            if i == N1-1 or j == N2-1:
+                break
+            n += 1
+            st[n] = spikes1[i]
+            c[n] = 0
+            n += 1
+            st[n] = spikes1[i]
+            c[n] = 1
+    c[0] = c[2]
+    st[0] = spikes1[0]
+    st[-1] = spikes1[-1]
+
+    return st, c
+
+
+############################################################
 # add_piece_wise_const_python
 ############################################################
 def add_piece_wise_const_python(x1, y1, x2, y2):