From 017893e782d1747e1f031131077a40a48f882e86 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Mon, 11 May 2015 14:32:48 +0200 Subject: treatment of empty spike trains in isi functions --- test/test_empty.py | 62 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) create mode 100644 test/test_empty.py (limited to 'test/test_empty.py') diff --git a/test/test_empty.py b/test/test_empty.py new file mode 100644 index 0000000..22982c7 --- /dev/null +++ b/test/test_empty.py @@ -0,0 +1,62 @@ +""" test_empty.py + +Tests the distance measure for empty spike trains + +Copyright 2014, Mario Mulansky + +Distributed under the BSD License + +""" + +from __future__ import print_function +import numpy as np +from copy import copy +from numpy.testing import assert_equal, assert_almost_equal, \ + assert_array_equal, assert_array_almost_equal + +import pyspike as spk +from pyspike import SpikeTrain + + +def test_get_non_empty(): + st = SpikeTrain([], edges=(0.0, 1.0)) + spikes = st.get_spikes_non_empty() + assert_array_equal(spikes, [0.0, 1.0]) + + st = SpikeTrain([0.5, ], edges=(0.0, 1.0)) + spikes = st.get_spikes_non_empty() + assert_array_equal(spikes, [0.0, 0.5, 1.0]) + + +def test_isi_empty(): + st1 = SpikeTrain([], edges=(0.0, 1.0)) + st2 = SpikeTrain([], edges=(0.0, 1.0)) + d = spk.isi_distance(st1, st2) + assert_equal(d, 0.0) + prof = spk.isi_profile(st1, st2) + assert_equal(d, prof.avrg()) + assert_array_equal(prof.x, [0.0, 1.0]) + assert_array_equal(prof.y, [0.0, ]) + + st1 = SpikeTrain([], edges=(0.0, 1.0)) + st2 = SpikeTrain([0.4, ], edges=(0.0, 1.0)) + d = spk.isi_distance(st1, st2) + assert_equal(d, 0.6*0.4+0.4*0.6) + prof = spk.isi_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.6, 0.4]) + + st1 = SpikeTrain([0.6, ], edges=(0.0, 1.0)) + st2 = SpikeTrain([0.4, ], edges=(0.0, 1.0)) + d = spk.isi_distance(st1, st2) + assert_almost_equal(d, 0.2/0.6*0.4 + 0.0 + 0.2/0.6*0.4, decimal=15) + prof = spk.isi_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, [0.2/0.6, 0.0, 0.2/0.6], decimal=15) + + +if __name__ == "__main__": + test_get_non_empty() + test_isi_empty() -- cgit v1.2.3 From b6521869cad89eae119391557bfa57e818dc9894 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Mon, 11 May 2015 15:54:19 +0200 Subject: treatment of empty spike trains in spike distance --- pyspike/spike_distance.py | 13 +++++++------ test/test_empty.py | 48 ++++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 54 insertions(+), 7 deletions(-) (limited to 'test/test_empty.py') diff --git a/pyspike/spike_distance.py b/pyspike/spike_distance.py index d727fa2..ac2d260 100644 --- a/pyspike/spike_distance.py +++ b/pyspike/spike_distance.py @@ -41,10 +41,11 @@ Falling back to slow python backend.") from cython.python_backend import spike_distance_python \ as spike_profile_impl - times, y_starts, y_ends = spike_profile_impl(spike_train1.spikes, - spike_train2.spikes, - spike_train1.t_start, - spike_train1.t_end) + times, y_starts, y_ends = spike_profile_impl( + spike_train1.get_spikes_non_empty(), + spike_train2.get_spikes_non_empty(), + spike_train1.t_start, spike_train1.t_end) + return PieceWiseLinFunc(times, y_starts, y_ends) @@ -74,8 +75,8 @@ def spike_distance(spike_train1, spike_train2, interval=None): try: from cython.cython_distances import spike_distance_cython \ as spike_distance_impl - return spike_distance_impl(spike_train1.spikes, - spike_train2.spikes, + return spike_distance_impl(spike_train1.get_spikes_non_empty(), + spike_train2.get_spikes_non_empty(), spike_train1.t_start, spike_train1.t_end) except ImportError: diff --git a/test/test_empty.py b/test/test_empty.py index 22982c7..e08bd1a 100644 --- a/test/test_empty.py +++ b/test/test_empty.py @@ -10,7 +10,6 @@ Distributed under the BSD License from __future__ import print_function import numpy as np -from copy import copy from numpy.testing import assert_equal, assert_almost_equal, \ assert_array_equal, assert_array_almost_equal @@ -57,6 +56,53 @@ def test_isi_empty(): assert_array_almost_equal(prof.y, [0.2/0.6, 0.0, 0.2/0.6], decimal=15) +def test_spike_empty(): + st1 = SpikeTrain([], edges=(0.0, 1.0)) + st2 = SpikeTrain([], edges=(0.0, 1.0)) + d = spk.spike_distance(st1, st2) + assert_equal(d, 0.0) + prof = spk.spike_profile(st1, st2) + assert_equal(d, prof.avrg()) + assert_array_equal(prof.x, [0.0, 1.0]) + assert_array_equal(prof.y1, [0.0, ]) + assert_array_equal(prof.y2, [0.0, ]) + + st1 = SpikeTrain([], edges=(0.0, 1.0)) + st2 = SpikeTrain([0.4, ], edges=(0.0, 1.0)) + d = spk.spike_distance(st1, st2) + assert_almost_equal(d, 0.4*0.4*1.0/(0.4+1.0)**2 + 0.6*0.4*1.0/(0.6+1.0)**2, + decimal=15) + prof = spk.spike_profile(st1, st2) + assert_equal(d, prof.avrg()) + assert_array_equal(prof.x, [0.0, 0.4, 1.0]) + assert_array_almost_equal(prof.y1, [0.0, 2*0.4*1.0/(0.6+1.0)**2], + decimal=15) + assert_array_almost_equal(prof.y2, [2*0.4*1.0/(0.4+1.0)**2, 0.0], + decimal=15) + + st1 = SpikeTrain([0.6, ], edges=(0.0, 1.0)) + st2 = SpikeTrain([0.4, ], edges=(0.0, 1.0)) + d = spk.spike_distance(st1, st2) + s1 = np.array([0.0, 0.4*0.2/0.6, 0.2, 0.0]) + s2 = np.array([0.0, 0.2, 0.2*0.4/0.6, 0.0]) + isi1 = np.array([0.6, 0.6, 0.4]) + isi2 = np.array([0.4, 0.6, 0.6]) + expected_y1 = (s1[:-1]*isi2+s2[:-1]*isi1) / (0.5*(isi1+isi2)**2) + expected_y2 = (s1[1:]*isi2+s2[1:]*isi1) / (0.5*(isi1+isi2)**2) + expected_times = np.array([0.0, 0.4, 0.6, 1.0]) + expected_spike_val = sum((expected_times[1:] - expected_times[:-1]) * + (expected_y1+expected_y2)/2) + expected_spike_val /= (expected_times[-1]-expected_times[0]) + + assert_almost_equal(d, expected_spike_val, decimal=15) + prof = spk.spike_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.y1, expected_y1, decimal=15) + assert_array_almost_equal(prof.y2, expected_y2, decimal=15) + + if __name__ == "__main__": test_get_non_empty() test_isi_empty() + test_spike_empty() -- cgit v1.2.3 From 3b10b416940ae674df6d9dff8cdddc31085d8cf5 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Mon, 11 May 2015 17:03:29 +0200 Subject: updated copyright year --- test/test_empty.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'test/test_empty.py') diff --git a/test/test_empty.py b/test/test_empty.py index e08bd1a..42c3716 100644 --- a/test/test_empty.py +++ b/test/test_empty.py @@ -2,7 +2,7 @@ Tests the distance measure for empty spike trains -Copyright 2014, Mario Mulansky +Copyright 2015, Mario Mulansky Distributed under the BSD License -- cgit v1.2.3 From bec2529367f1bdd5dac6d6fbaec560a30feec3c7 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Mon, 11 May 2015 17:41:08 +0200 Subject: treatment of empty spike trains in spike sync --- pyspike/DiscreteFunc.py | 4 ++++ pyspike/cython/cython_distances.pyx | 13 +++++++++---- pyspike/cython/cython_profiles.pyx | 9 +++++---- pyspike/cython/python_backend.py | 2 +- test/test_empty.py | 39 +++++++++++++++++++++++++++++++++++++ 5 files changed, 58 insertions(+), 9 deletions(-) (limited to 'test/test_empty.py') 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() -- cgit v1.2.3 From a35402c208bd0ad31e5e60b6ddc55a3470e7bdde Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Mon, 11 May 2015 17:54:02 +0200 Subject: bugfix: spike_sync=1 for empty spike trains --- pyspike/DiscreteFunc.py | 4 ++-- pyspike/cython/cython_distances.pyx | 3 ++- pyspike/cython/cython_profiles.pyx | 12 ++++++++---- pyspike/cython/python_backend.py | 12 ++++++++---- test/test_empty.py | 4 ++-- 5 files changed, 22 insertions(+), 13 deletions(-) (limited to 'test/test_empty.py') diff --git a/pyspike/DiscreteFunc.py b/pyspike/DiscreteFunc.py index 6ade87e..17153ee 100644 --- a/pyspike/DiscreteFunc.py +++ b/pyspike/DiscreteFunc.py @@ -137,8 +137,8 @@ class DiscreteFunc(object): """ if len(self.y) <= 2: - # no actual values in the profile, return spike sync of 0 - return 0.0, 1.0 + # no actual values in the profile, return spike sync of 1 + return 1.0, 1.0 def get_indices(ival): """ Retuns the indeces surrounding the given interval""" diff --git a/pyspike/cython/cython_distances.pyx b/pyspike/cython/cython_distances.pyx index 16780f2..c4f2349 100644 --- a/pyspike/cython/cython_distances.pyx +++ b/pyspike/cython/cython_distances.pyx @@ -391,7 +391,8 @@ def coincidence_value_cython(double[:] spikes1, double[:] spikes2, coinc += 2 if coinc == 0 and mp == 0: - # empty spike trains -> set mp to one to avoid 0/0 + # empty spike trains -> spike sync = 1 by definition + coinc = 1 mp = 1 return coinc, mp diff --git a/pyspike/cython/cython_profiles.pyx b/pyspike/cython/cython_profiles.pyx index d937a02..f9893eb 100644 --- a/pyspike/cython/cython_profiles.pyx +++ b/pyspike/cython/cython_profiles.pyx @@ -416,9 +416,13 @@ def coincidence_profile_cython(double[:] spikes1, double[:] spikes2, st[0] = t_start st[len(st)-1] = t_end - c[0] = c[1] - c[len(c)-1] = c[len(c)-2] - mp[0] = mp[1] - mp[len(mp)-1] = mp[len(mp)-2] + if N1 + N2 > 0: + c[0] = c[1] + c[len(c)-1] = c[len(c)-2] + mp[0] = mp[1] + mp[len(mp)-1] = mp[len(mp)-2] + else: + c[0] = 1 + c[1] = 1 return st, c, mp diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py index 830dc69..69a420f 100644 --- a/pyspike/cython/python_backend.py +++ b/pyspike/cython/python_backend.py @@ -399,10 +399,14 @@ def coincidence_python(spikes1, spikes2, t_start, t_end, max_tau): st[0] = t_start st[len(st)-1] = t_end - c[0] = c[1] - c[len(c)-1] = c[len(c)-2] - mp[0] = mp[1] - mp[len(mp)-1] = mp[len(mp)-2] + if N1 + N2 > 0: + c[0] = c[1] + c[len(c)-1] = c[len(c)-2] + mp[0] = mp[1] + mp[len(mp)-1] = mp[len(mp)-2] + else: + c[0] = 1 + c[1] = 1 return st, c, mp diff --git a/test/test_empty.py b/test/test_empty.py index b31d8a4..48be25d 100644 --- a/test/test_empty.py +++ b/test/test_empty.py @@ -106,11 +106,11 @@ 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) + assert_equal(d, 1.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]) + assert_array_equal(prof.y, [1.0, 1.0]) st1 = SpikeTrain([], edges=(0.0, 1.0)) st2 = SpikeTrain([0.4, ], edges=(0.0, 1.0)) -- cgit v1.2.3