From b09561705ab9c67c93a384248f7c3bc9ad5bdd32 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Wed, 3 Feb 2016 13:07:26 +0100 Subject: fixed spike-sync bug fixed ugly bugs in code for computing multi-variate spike sync profile and multi-variate spike sync value. --- pyspike/DiscreteFunc.py | 9 ++++---- pyspike/cython/cython_add.pyx | 33 ++++++++++++++------------- pyspike/cython/cython_distances.pyx | 5 ---- pyspike/cython/python_backend.py | 33 +++++++++++++++------------ pyspike/spike_sync.py | 10 ++++++-- test/numeric/test_regression_random_spikes.py | 27 +++++++++++++++------- test/test_distance.py | 33 +++++++++++++++++++++++++++ 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: -- cgit v1.2.3