1 files changed, 64 insertions, 0 deletions
diff --git a/pyspike/cython/cython_distances.pyx b/pyspike/cython/cython_distances.pyx
index 65c2872..bf90638 100644
--- a/pyspike/cython/cython_distances.pyx
+++ b/pyspike/cython/cython_distances.pyx
@@ -326,3 +326,67 @@ def spike_distance_cython(double[:] t1, double[:] t2,
     # use only the data added above 
     # could be less than original length due to equal spike times
     return spike_value / (t_end-t_start)
+
+
+
+############################################################
+# get_tau
+############################################################
+cdef inline double get_tau(double[:] spikes1, double[:] spikes2,
+                           int i, int j, double max_tau):
+    cdef double m = 1E100   # some huge number
+    cdef int N1 = spikes1.shape[0]-1  # len(spikes1)-1
+    cdef int N2 = spikes2.shape[0]-1  # len(spikes2)-1
+    if i < N1 and i > -1:
+        m = fmin(m, spikes1[i+1]-spikes1[i])
+    if j < N2 and j > -1:
+        m = fmin(m, spikes2[j+1]-spikes2[j])
+    if i > 0:
+        m = fmin(m, spikes1[i]-spikes1[i-1])
+    if j > 0:
+        m = fmin(m, spikes2[j]-spikes2[j-1])
+    m *= 0.5
+    if max_tau > 0.0:
+        m = fmin(m, max_tau)
+    return m
+    
+
+############################################################
+# coincidence_value_cython
+############################################################
+def coincidence_value_cython(double[:] spikes1, double[:] spikes2,
+                             double t_start, double t_end, double max_tau):
+
+    cdef int N1 = len(spikes1)
+    cdef int N2 = len(spikes2)
+    cdef int i = -1
+    cdef int j = -1
+    cdef double coinc = 0.0
+    cdef double mp = 0.0
+    cdef double tau
+    while i + j < N1 + N2 - 2:
+        if (i < N1-1) and (j == N2-1 or spikes1[i+1] < spikes2[j+1]):
+            i += 1
+            mp += 1
+            tau = get_tau(spikes1, spikes2, i, j, max_tau)
+            if j > -1 and spikes1[i]-spikes2[j] < tau:
+                # coincidence between the current spike and the previous spike
+                # both get marked with 1
+                coinc += 2
+        elif (j < N2-1) and (i == N1-1 or spikes1[i+1] > spikes2[j+1]):
+            j += 1
+            mp += 1
+            tau = get_tau(spikes1, spikes2, i, j, max_tau)
+            if i > -1 and spikes2[j]-spikes1[i] < tau:
+                # coincidence between the current spike and the previous spike
+                # both get marked with 1
+                coinc += 2
+        else:   # spikes1[i+1] = spikes2[j+1]
+            # advance in both spike trains
+            j += 1
+            i += 1
+            # add only one event, but with coincidence 2 and multiplicity 2
+            mp += 2
+            coinc += 2
+
+    return coinc, mp