fixed spike-sync bug

fixed ugly bugs in code for computing multi-variate spike sync profile and
multi-variate spike sync value.

7 files changed, 99 insertions, 51 deletions

diff --git a/pyspike/DiscreteFunc.py b/pyspike/DiscreteFunc.py
index fe97bc2..caad290 100644
--- a/pyspike/DiscreteFunc.py
+++ b/pyspike/DiscreteFunc.py
@@ -170,10 +170,6 @@ expected."
                     start_ind, end_ind = get_indices(ival)
                     value += np.sum(self.y[start_ind:end_ind])
                     multiplicity += np.sum(self.mp[start_ind:end_ind])
-        if multiplicity == 0.0:
-            # empty profile, return spike sync of 1
-            value = 1.0
-            multiplicity = 1.0
         return (value, multiplicity)
 
     def avrg(self, interval=None, normalize=True):
@@ -190,7 +186,10 @@ expected."
         """
         val, mp = self.integral(interval)
         if normalize:
-            return val/mp
+            if mp > 0:
+                return val/mp
+            else:
+                return 1.0
         else:
             return val
 
diff --git a/pyspike/cython/cython_add.pyx b/pyspike/cython/cython_add.pyx
index 8da1e53..25f1181 100644
--- a/pyspike/cython/cython_add.pyx
+++ b/pyspike/cython/cython_add.pyx
@@ -182,8 +182,8 @@ def add_discrete_function_cython(double[:] x1, double[:] y1, double[:] mp1,
     cdef int index1 = 0
     cdef int index2 = 0
     cdef int index = 0
-    cdef int N1 = len(y1)
-    cdef int N2 = len(y2)
+    cdef int N1 = len(y1)-1
+    cdef int N2 = len(y2)-1
     x_new[0] = x1[0]
     while (index1+1 < N1) and (index2+1 < N2):
         if x1[index1+1] < x2[index2+1]:
@@ -206,20 +206,21 @@ def add_discrete_function_cython(double[:] x1, double[:] y1, double[:] mp1,
             y_new[index] = y1[index1] + y2[index2]
             mp_new[index] = mp1[index1] + mp2[index2]
     # one array reached the end -> copy the contents of the other to the end
-    if index1+1 < len(y1):
-        x_new[index+1:index+1+len(x1)-index1-1] = x1[index1+1:]
-        y_new[index+1:index+1+len(x1)-index1-1] = y1[index1+1:]
-        mp_new[index+1:index+1+len(x1)-index1-1] = mp1[index1+1:]
-        index += len(x1)-index1-1
-    elif index2+1 < len(y2):
-        x_new[index+1:index+1+len(x2)-index2-1] = x2[index2+1:]
-        y_new[index+1:index+1+len(x2)-index2-1] = y2[index2+1:]
-        mp_new[index+1:index+1+len(x2)-index2-1] = mp2[index2+1:]
-        index += len(x2)-index2-1
-    # else:  # both arrays reached the end simultaneously
-    #     x_new[index+1] = x1[-1]
-    #     y_new[index+1] = y1[-1] + y2[-1]
-    #     mp_new[index+1] = mp1[-1] + mp2[-1]
+    if index1+1 < N1:
+        x_new[index+1:index+1+N1-index1] = x1[index1+1:]
+        y_new[index+1:index+1+N1-index1] = y1[index1+1:]
+        mp_new[index+1:index+1+N1-index1] = mp1[index1+1:]
+        index += N1-index1
+    elif index2+1 < N2:
+        x_new[index+1:index+1+N2-index2] = x2[index2+1:]
+        y_new[index+1:index+1+N2-index2] = y2[index2+1:]
+        mp_new[index+1:index+1+N2-index2] = mp2[index2+1:]
+        index += N2-index2
+    else:  # both arrays reached the end simultaneously
+        x_new[index+1] = x1[index1+1]
+        y_new[index+1] = y1[index1+1] + y2[index2+1]
+        mp_new[index+1] = mp1[index1+1] + mp2[index2+1]
+        index += 1
 
     y_new[0] = y_new[1]
     mp_new[0] = mp_new[1]
diff --git a/pyspike/cython/cython_distances.pyx b/pyspike/cython/cython_distances.pyx
index 7dc9cb9..ac5f226 100644
--- a/pyspike/cython/cython_distances.pyx
+++ b/pyspike/cython/cython_distances.pyx
@@ -428,9 +428,4 @@ def coincidence_value_cython(double[:] spikes1, double[:] spikes2,
             mp += 2
             coinc += 2
 
-    if coinc == 0 and mp == 0:
-        # empty spike trains -> spike sync = 1 by definition
-        coinc = 1
-        mp = 1
-
     return coinc, mp
diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py
index a1e3a65..6b7209a 100644
--- a/pyspike/cython/python_backend.py
+++ b/pyspike/cython/python_backend.py
@@ -560,7 +560,9 @@ def add_discrete_function_python(x1, y1, mp1, x2, y2, mp2):
     index1 = 0
     index2 = 0
     index = 0
-    while (index1+1 < len(y1)) and (index2+1 < len(y2)):
+    N1 = len(x1)-1
+    N2 = len(x2)-1
+    while (index1+1 < N1) and (index2+1 < N2):
         if x1[index1+1] < x2[index2+1]:
             index1 += 1
             index += 1
@@ -581,20 +583,21 @@ def add_discrete_function_python(x1, y1, mp1, x2, y2, mp2):
             y_new[index] = y1[index1] + y2[index2]
             mp_new[index] = mp1[index1] + mp2[index2]
     # one array reached the end -> copy the contents of the other to the end
-    if index1+1 < len(y1):
-        x_new[index+1:index+1+len(x1)-index1-1] = x1[index1+1:]
-        y_new[index+1:index+1+len(x1)-index1-1] = y1[index1+1:]
-        mp_new[index+1:index+1+len(x1)-index1-1] = mp1[index1+1:]
-        index += len(x1)-index1-1
-    elif index2+1 < len(y2):
-        x_new[index+1:index+1+len(x2)-index2-1] = x2[index2+1:]
-        y_new[index+1:index+1+len(x2)-index2-1] = y2[index2+1:]
-        mp_new[index+1:index+1+len(x2)-index2-1] = mp2[index2+1:]
-        index += len(x2)-index2-1
-    # else:  # both arrays reached the end simultaneously
-    #     x_new[index+1] = x1[-1]
-    #     y_new[index+1] = y1[-1] + y2[-1]
-    #     mp_new[index+1] = mp1[-1] + mp2[-1]
+    if index1+1 < N1:
+        x_new[index+1:index+1+N1-index1] = x1[index1+1:]
+        y_new[index+1:index+1+N1-index1] = y1[index1+1:]
+        mp_new[index+1:index+1+N1-index1] = mp1[index1+1:]
+        index += N1-index1
+    elif index2+1 < N2:
+        x_new[index+1:index+1+N2-index2] = x2[index2+1:]
+        y_new[index+1:index+1+N2-index2] = y2[index2+1:]
+        mp_new[index+1:index+1+N2-index2] = mp2[index2+1:]
+        index += N2-index2
+    else:  # both arrays reached the end simultaneously
+        x_new[index+1] = x1[-1]
+        y_new[index+1] = y1[-1] + y2[-1]
+        mp_new[index+1] = mp1[-1] + mp2[-1]
+        index += 1
 
     y_new[0] = y_new[1]
     mp_new[0] = mp_new[1]
diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py
index 3dc29ff..802b98c 100644
--- a/pyspike/spike_sync.py
+++ b/pyspike/spike_sync.py
@@ -119,7 +119,10 @@ def spike_sync(spike_train1, spike_train2, interval=None, max_tau=None):
 
     """
     c, mp = _spike_sync_values(spike_train1, spike_train2, interval, max_tau)
-    return 1.0*c/mp
+    if mp == 0:
+        return 1.0
+    else:
+        return 1.0*c/mp
 
 
 ############################################################
@@ -187,7 +190,10 @@ def spike_sync_multi(spike_trains, indices=None, interval=None, max_tau=None):
         coincidence += c
         mp += m
 
-    return coincidence/mp
+    if mp == 0.0:
+        return 1.0
+    else:
+        return coincidence/mp
 
 
 ############################################################
diff --git a/test/numeric/test_regression_random_spikes.py b/test/numeric/test_regression_random_spikes.py
index 6156bb4..73d39a6 100644
--- a/test/numeric/test_regression_random_spikes.py
+++ b/test/numeric/test_regression_random_spikes.py
@@ -35,7 +35,9 @@ def test_regression_random():
 
         spike = spk.spike_distance_multi(spike_trains)
         spike_prof = spk.spike_profile_multi(spike_trains).avrg()
-        # spike_sync = spk.spike_sync_multi(spike_trains)
+
+        spike_sync = spk.spike_sync_multi(spike_trains)
+        spike_sync_prof = spk.spike_sync_profile_multi(spike_trains).avrg()
 
         assert_almost_equal(isi, results_cSPIKY[i][0], decimal=14,
                             err_msg="Index: %d, ISI" % i)
@@ -47,6 +49,9 @@ def test_regression_random():
         assert_almost_equal(spike_prof, results_cSPIKY[i][1], decimal=14,
                             err_msg="Index: %d, SPIKE" % i)
 
+        assert_almost_equal(spike_sync, spike_sync_prof, decimal=14,
+                            err_msg="Index: %d, SPIKE-Sync" % i)
+
 
 def check_regression_dataset(spike_file="benchmark.mat",
                              spikes_name="spikes",
@@ -109,19 +114,25 @@ def check_single_spike_train_set(index):
     spike_train_data = spike_train_sets[index]
 
     spike_trains = []
+    N = 0
     for spikes in spike_train_data[0]:
-        print("Spikes:", spikes.flatten())
-        spike_trains.append(spk.SpikeTrain(spikes.flatten(), 100.0))
+        N += len(spikes.flatten())
+        print("Spikes:", len(spikes.flatten()))
+        spikes_array = spikes.flatten()
+        if len(spikes_array > 0) and (spikes_array[-1] > 100.0):
+            spikes_array[-1] = 100.0
+        spike_trains.append(spk.SpikeTrain(spikes_array, 100.0))
+        print(spike_trains[-1].spikes)
 
-    print(spk.spike_distance_multi(spike_trains))
+    print(N)
 
-    print(results_cSPIKY[index][1])
+    print(spk.spike_sync_multi(spike_trains))
 
-    print(spike_trains[1].spikes)
+    print(spk.spike_sync_profile_multi(spike_trains).integral())
 
 
 if __name__ == "__main__":
 
-    test_regression_random()
+    # test_regression_random()
     # check_regression_dataset()
-    # check_single_spike_train_set(7633)
+    check_single_spike_train_set(4)
diff --git a/test/test_distance.py b/test/test_distance.py
index 083d8a3..fe09f34 100644
--- a/test/test_distance.py
+++ b/test/test_distance.py
@@ -217,6 +217,28 @@ def test_spike_sync():
     assert_almost_equal(spk.spike_sync(spikes1, spikes2),
                         0.4, decimal=16)
 
+    spikes1 = SpikeTrain([1.0, 2.0, 4.0], 4.0)
+    spikes2 = SpikeTrain([3.8], 4.0)
+    spikes3 = SpikeTrain([3.9, ], 4.0)
+
+    expected_x = np.array([0.0, 1.0, 2.0, 3.8, 4.0, 4.0])
+    expected_y = np.array([0.0, 0.0, 0.0, 1.0, 1.0, 1.0])
+
+    f = spk.spike_sync_profile(spikes1, spikes2)
+
+    assert_array_almost_equal(f.x, expected_x, decimal=16)
+    assert_array_almost_equal(f.y, expected_y, decimal=16)
+
+    f2 = spk.spike_sync_profile(spikes2, spikes3)
+
+    i1 = f.integral()
+    i2 = f2.integral()
+    f.add(f2)
+    i12 = f.integral()
+
+    assert_equal(i1[0]+i2[0], i12[0])
+    assert_equal(i1[1]+i2[1], i12[1])
+
 
 def check_multi_profile(profile_func, profile_func_multi, dist_func_multi):
     # generate spike trains:
@@ -313,6 +335,17 @@ def test_multi_spike_sync():
     assert_equal(np.sum(f.y[1:-1]), 39932)
     assert_equal(np.sum(f.mp[1:-1]), 85554)
 
+    # example with 2 empty spike trains
+    sts = []
+    sts.append(SpikeTrain([1, 9], [0, 10]))
+    sts.append(SpikeTrain([1, 3], [0, 10]))
+    sts.append(SpikeTrain([], [0, 10]))
+    sts.append(SpikeTrain([], [0, 10]))
+
+    assert_almost_equal(spk.spike_sync_multi(sts), 1.0/6.0, decimal=15)
+    assert_almost_equal(spk.spike_sync_profile_multi(sts).avrg(), 1.0/6.0,
+                        decimal=15)
+
 
 def check_dist_matrix(dist_func, dist_matrix_func):
     # generate spike trains: