treatment of empty spike trains in spike sync

5 files changed, 58 insertions, 9 deletions

diff --git a/pyspike/DiscreteFunc.py b/pyspike/DiscreteFunc.py
index dfe2cab..6ade87e 100644
--- a/pyspike/DiscreteFunc.py
+++ b/pyspike/DiscreteFunc.py
@@ -136,6 +136,10 @@ class DiscreteFunc(object):
         :rtype: pair of float
         """
 
+        if len(self.y) <= 2:
+            # no actual values in the profile, return spike sync of 0
+            return 0.0, 1.0
+
         def get_indices(ival):
             """ Retuns the indeces surrounding the given interval"""
             start_ind = np.searchsorted(self.x, ival[0], side='right')
diff --git a/pyspike/cython/cython_distances.pyx b/pyspike/cython/cython_distances.pyx
index bf90638..16780f2 100644
--- a/pyspike/cython/cython_distances.pyx
+++ b/pyspike/cython/cython_distances.pyx
@@ -333,8 +333,8 @@ def spike_distance_cython(double[:] t1, double[:] t2,
 # 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
+                           int i, int j, double interval, double max_tau):
+    cdef double m = interval   # use interval length as initial tau
     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:
@@ -363,12 +363,13 @@ def coincidence_value_cython(double[:] spikes1, double[:] spikes2,
     cdef int j = -1
     cdef double coinc = 0.0
     cdef double mp = 0.0
+    cdef double interval = t_end - t_start
     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)
+            tau = get_tau(spikes1, spikes2, i, j, interval, 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
@@ -376,7 +377,7 @@ def coincidence_value_cython(double[:] spikes1, double[:] spikes2,
         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)
+            tau = get_tau(spikes1, spikes2, i, j, interval, 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
@@ -389,4 +390,8 @@ def coincidence_value_cython(double[:] spikes1, double[:] spikes2,
             mp += 2
             coinc += 2
 
+    if coinc == 0 and mp == 0:
+        # empty spike trains -> set mp to one to avoid 0/0
+        mp = 1
+
     return coinc, mp
diff --git a/pyspike/cython/cython_profiles.pyx b/pyspike/cython/cython_profiles.pyx
index 3690091..d937a02 100644
--- a/pyspike/cython/cython_profiles.pyx
+++ b/pyspike/cython/cython_profiles.pyx
@@ -345,8 +345,8 @@ def spike_profile_cython(double[:] t1, double[:] t2,
 # 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
+                           int i, int j, double interval, double max_tau):
+    cdef double m = interval   # use interval as initial tau
     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:
@@ -377,12 +377,13 @@ def coincidence_profile_cython(double[:] spikes1, double[:] spikes2,
     cdef double[:] st = np.zeros(N1 + N2 + 2)  # spike times
     cdef double[:] c = np.zeros(N1 + N2 + 2)   # coincidences
     cdef double[:] mp = np.ones(N1 + N2 + 2)   # multiplicity
+    cdef double interval = t_end - t_start
     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
             n += 1
-            tau = get_tau(spikes1, spikes2, i, j, max_tau)
+            tau = get_tau(spikes1, spikes2, i, j, interval, max_tau)
             st[n] = spikes1[i]
             if j > -1 and spikes1[i]-spikes2[j] < tau:
                 # coincidence between the current spike and the previous spike
@@ -392,7 +393,7 @@ def coincidence_profile_cython(double[:] spikes1, double[:] spikes2,
         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, interval, max_tau)
             st[n] = spikes2[j]
             if i > -1 and spikes2[j]-spikes1[i] < tau:
                 # coincidence between the current spike and the previous spike
diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py
index 1fd8c42..830dc69 100644
--- a/pyspike/cython/python_backend.py
+++ b/pyspike/cython/python_backend.py
@@ -340,7 +340,7 @@ def cumulative_sync_python(spikes1, spikes2):
 def coincidence_python(spikes1, spikes2, t_start, t_end, max_tau):
 
     def get_tau(spikes1, spikes2, i, j, max_tau):
-        m = 1E100   # some huge number
+        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:
diff --git a/test/test_empty.py b/test/test_empty.py
index 42c3716..b31d8a4 100644
--- a/test/test_empty.py
+++ b/test/test_empty.py
@@ -102,7 +102,46 @@ def test_spike_empty():
     assert_array_almost_equal(prof.y2, expected_y2, decimal=15)
 
 
+def test_spike_sync_empty():
+    st1 = SpikeTrain([], edges=(0.0, 1.0))
+    st2 = SpikeTrain([], edges=(0.0, 1.0))
+    d = spk.spike_sync(st1, st2)
+    assert_equal(d, 0.0)
+    prof = spk.spike_sync_profile(st1, st2)
+    assert_equal(d, prof.avrg())
+    assert_array_equal(prof.x, [0.0, 1.0])
+    assert_array_equal(prof.y, [0.0, 0.0])
+
+    st1 = SpikeTrain([], edges=(0.0, 1.0))
+    st2 = SpikeTrain([0.4, ], edges=(0.0, 1.0))
+    d = spk.spike_sync(st1, st2)
+    assert_equal(d, 0.0)
+    prof = spk.spike_sync_profile(st1, st2)
+    assert_equal(d, prof.avrg())
+    assert_array_equal(prof.x, [0.0, 0.4, 1.0])
+    assert_array_equal(prof.y, [0.0, 0.0, 0.0])
+
+    st1 = SpikeTrain([0.6, ], edges=(0.0, 1.0))
+    st2 = SpikeTrain([0.4, ], edges=(0.0, 1.0))
+    d = spk.spike_sync(st1, st2)
+    assert_almost_equal(d, 1.0, decimal=15)
+    prof = spk.spike_sync_profile(st1, st2)
+    assert_equal(d, prof.avrg())
+    assert_array_almost_equal(prof.x, [0.0, 0.4, 0.6, 1.0], decimal=15)
+    assert_array_almost_equal(prof.y, [1.0, 1.0, 1.0, 1.0], decimal=15)
+
+    st1 = SpikeTrain([0.2, ], edges=(0.0, 1.0))
+    st2 = SpikeTrain([0.8, ], edges=(0.0, 1.0))
+    d = spk.spike_sync(st1, st2)
+    assert_almost_equal(d, 0.0, decimal=15)
+    prof = spk.spike_sync_profile(st1, st2)
+    assert_equal(d, prof.avrg())
+    assert_array_almost_equal(prof.x, [0.0, 0.2, 0.8, 1.0], decimal=15)
+    assert_array_almost_equal(prof.y, [0.0, 0.0, 0.0, 0.0], decimal=15)
+
+
 if __name__ == "__main__":
     test_get_non_empty()
     test_isi_empty()
     test_spike_empty()
+    test_spike_sync_empty()