added coincidence single profile to python backend

missing function in python backend added, identified and fixed a bug in the
implementation as well

4 files changed, 53 insertions, 21 deletions

diff --git a/pyspike/cython/cython_profiles.pyx b/pyspike/cython/cython_profiles.pyx
index eb4d157..aa24db4 100644
--- a/pyspike/cython/cython_profiles.pyx
+++ b/pyspike/cython/cython_profiles.pyx
@@ -473,7 +473,7 @@ def coincidence_single_profile_cython(double[:] spikes1, double[:] spikes2,
         if j > -1 and fabs(spikes1[i]-spikes2[j]) < tau:
             # current spike in st1 is coincident
             c[i] = 1
-        if j < N2-1 and spikes2[j] < spikes1[i]:
+        if j < N2-1 and (j < 0 or spikes2[j] < spikes1[i]):
             # in case spikes2[j] is before spikes1[i] it has to be the one 
             # right before (see above), hence we move one forward and also 
             # check the next spike
diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py
index 6b7209a..a4a0c9e 100644
--- a/pyspike/cython/python_backend.py
+++ b/pyspike/cython/python_backend.py
@@ -3,7 +3,7 @@
 Collection of python functions that can be used instead of the cython
 implementation.
 
-Copyright 2014, Mario Mulansky <mario.mulansky@gmx.net>
+Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net>
 
 Distributed under the BSD License
 
@@ -356,26 +356,27 @@ def cumulative_sync_python(spikes1, spikes2):
     return st, c
 
 
+def get_tau(spikes1, spikes2, i, j, max_tau, init_tau):
+    m = init_tau
+    if i < len(spikes1)-1 and i > -1:
+        m = min(m, spikes1[i+1]-spikes1[i])
+    if j < len(spikes2)-1 and j > -1:
+        m = min(m, spikes2[j+1]-spikes2[j])
+    if i > 0:
+        m = min(m, spikes1[i]-spikes1[i-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
+
+
 ############################################################
 # coincidence_python
 ############################################################
 def coincidence_python(spikes1, spikes2, t_start, t_end, max_tau):
 
-    def get_tau(spikes1, spikes2, i, j, max_tau):
-        m = t_end - t_start   # use interval as initial tau
-        if i < len(spikes1)-1 and i > -1:
-            m = min(m, spikes1[i+1]-spikes1[i])
-        if j < len(spikes2)-1 and j > -1:
-            m = min(m, spikes2[j+1]-spikes2[j])
-        if i > 0:
-            m = min(m, spikes1[i]-spikes1[i-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 = -1
@@ -388,7 +389,7 @@ def coincidence_python(spikes1, spikes2, t_start, t_end, max_tau):
         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)
+            tau = get_tau(spikes1, spikes2, i, j, max_tau, t_end-t_start)
             st[n] = spikes1[i]
             if j > -1 and spikes1[i]-spikes2[j] < tau:
                 # coincidence between the current spike and the previous spike
@@ -398,7 +399,7 @@ def coincidence_python(spikes1, spikes2, t_start, t_end, max_tau):
         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)
+            tau = get_tau(spikes1, spikes2, i, j, max_tau, t_end-t_start)
             st[n] = spikes2[j]
             if i > -1 and spikes2[j]-spikes1[i] < tau:
                 # coincidence between the current spike and the previous spike
@@ -434,6 +435,36 @@ def coincidence_python(spikes1, spikes2, t_start, t_end, max_tau):
 
 
 ############################################################
+# coincidence_single_profile_cython
+############################################################
+def coincidence_single_python(spikes1, spikes2, t_start, t_end, max_tau):
+
+    N1 = len(spikes1)
+    N2 = len(spikes2)
+    j = -1
+    c = np.zeros(N1)   # coincidences
+    for i in xrange(N1):
+        while j < N2-1 and spikes2[j+1] < spikes1[i]:
+            # move forward until spikes2[j] is the last spike before spikes1[i]
+            # note that if spikes2[j] is after spikes1[i] we dont do anything
+            j += 1
+        tau = get_tau(spikes1, spikes2, i, j, max_tau, t_end-t_start)
+        if j > -1 and abs(spikes1[i]-spikes2[j]) < tau:
+            # current spike in st1 is coincident
+            c[i] = 1
+        if j < N2-1 and (j < 0 or spikes2[j] < spikes1[i]):
+            # in case spikes2[j] is before spikes1[i] it has to be the first or
+            # the one right before (see above), hence we move one forward and
+            # also check the next spike
+            j += 1
+            tau = get_tau(spikes1, spikes2, i, j, max_tau, t_end-t_start)
+            if abs(spikes2[j]-spikes1[i]) < tau:
+                # current spike in st1 is coincident
+                c[i] = 1
+    return c
+
+
+############################################################
 # add_piece_wise_const_python
 ############################################################
 def add_piece_wise_const_python(x1, y1, x2, y2):
diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py
index 1d2ecdb..66aa906 100644
--- a/pyspike/spike_sync.py
+++ b/pyspike/spike_sync.py
@@ -315,7 +315,7 @@ 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 cython.python_backend import coincidence_single_profile_python \
+        from cython.python_backend import coincidence_single_python \
             as coincidence_impl
 
     if max_tau is None:
diff --git a/test/test_sync_filter.py b/test/test_sync_filter.py
index 66ffcb6..e259903 100644
--- a/test/test_sync_filter.py
+++ b/test/test_sync_filter.py
@@ -28,12 +28,13 @@ def test_single_prof():
             coincidence_single_profile_cython as coincidence_impl
     except ImportError:
         from pyspike.cython.python_backend import \
-            coincidence_single_profile_python as coincidence_impl
+            coincidence_single_python as coincidence_impl
 
     sync_prof = spk.spike_sync_profile(SpikeTrain(st1, 5.0),
                                        SpikeTrain(st2, 5.0))
 
     coincidences = np.array(coincidence_impl(st1, st2, 0, 5.0, 0.0))
+    print(coincidences)
     for i, t in enumerate(st1):
         assert_equal(coincidences[i], sync_prof.y[sync_prof.x == t],
                      "At index %d" % i)