From 2e7351393927ba9e9e0c3b7b59d05e8aeeb41d1f Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Tue, 28 Apr 2015 16:11:11 +0200 Subject: edge correction for the ISI-distance --- pyspike/cython/cython_distance.pyx | 18 ++++++++++++------ pyspike/cython/python_backend.py | 18 ++++++++++++------ test/test_distance.py | 8 +++++--- 3 files changed, 29 insertions(+), 15 deletions(-) diff --git a/pyspike/cython/cython_distance.pyx b/pyspike/cython/cython_distance.pyx index a41d8e8..6ee0181 100644 --- a/pyspike/cython/cython_distance.pyx +++ b/pyspike/cython/cython_distance.pyx @@ -62,14 +62,16 @@ def isi_distance_cython(double[:] s1, double[:] s2, # first interspike interval - check if a spike exists at the start time if s1[0] > t_start: - nu1 = s1[0] - t_start + # edge correction + nu1 = fmax(s1[0]-t_start, s1[1]-s1[0]) index1 = -1 else: nu1 = s1[1]-s1[0] index1 = 0 if s2[0] > t_start: - nu2 = s2[0] - t_start + # edge correction + nu2 = fmax(s2[0]-t_start, s2[1]-s2[0]) index2 = -1 else: nu2 = s2[1]-s2[0] @@ -89,7 +91,8 @@ def isi_distance_cython(double[:] s1, double[:] s2, if index1 < N1-1: nu1 = s1[index1+1]-s1[index1] else: - nu1 = t_end-s1[index1] + # edge correction + nu1 = fmax(t_end-s1[index1], nu1) elif (index2 < N2-1) and ((index1 == N1-1) or (s1[index1+1] > s2[index2+1])): index2 += 1 @@ -97,7 +100,8 @@ def isi_distance_cython(double[:] s1, double[:] s2, if index2 < N2-1: nu2 = s2[index2+1]-s2[index2] else: - nu2 = t_end-s2[index2] + # edge correction + nu2 = fmax(t_end-s2[index2], nu2) else: # s1[index1+1] == s2[index2+1] index1 += 1 index2 += 1 @@ -105,11 +109,13 @@ def isi_distance_cython(double[:] s1, double[:] s2, if index1 < N1-1: nu1 = s1[index1+1]-s1[index1] else: - nu1 = t_end-s1[index1] + # edge correction + nu1 = fmax(t_end-s1[index1], nu1) if index2 < N2-1: nu2 = s2[index2+1]-s2[index2] else: - nu2 = t_end-s2[index2] + # edge correction + nu2 = fmax(t_end-s2[index2], nu2) # compute the corresponding isi-distance isi_values[index] = fabs(nu1 - nu2) / fmax(nu1, nu2) index += 1 diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py index 317b568..1fd8c42 100644 --- a/pyspike/cython/python_backend.py +++ b/pyspike/cython/python_backend.py @@ -27,13 +27,15 @@ def isi_distance_python(s1, s2, t_start, t_end): # the values have one entry less - the number of intervals between events isi_values = np.empty(len(spike_events) - 1) if s1[0] > t_start: - nu1 = s1[0] - t_start + # edge correction + nu1 = max(s1[0] - t_start, s1[1] - s1[0]) index1 = -1 else: nu1 = s1[1] - s1[0] index1 = 0 if s2[0] > t_start: - nu2 = s2[0] - t_start + # edge correction + nu2 = max(s2[0] - t_start, s2[1] - s2[0]) index2 = -1 else: nu2 = s2[1] - s2[0] @@ -49,7 +51,8 @@ def isi_distance_python(s1, s2, t_start, t_end): if index1 < N1-1: nu1 = s1[index1+1]-s1[index1] else: - nu1 = t_end-s1[index1] + # edge correction + nu1 = max(t_end-s1[N1-1], s1[N1-1]-s1[N1-2]) elif (index2 < N2-1) and (index1 == N1-1 or s1[index1+1] > s2[index2+1]): @@ -58,7 +61,8 @@ def isi_distance_python(s1, s2, t_start, t_end): if index2 < N2-1: nu2 = s2[index2+1]-s2[index2] else: - nu2 = t_end-s2[index2] + # edge correction + nu2 = max(t_end-s2[N2-1], s2[N2-1]-s2[N2-2]) else: # s1[index1 + 1] == s2[index2 + 1] index1 += 1 @@ -67,11 +71,13 @@ def isi_distance_python(s1, s2, t_start, t_end): if index1 < N1-1: nu1 = s1[index1+1]-s1[index1] else: - nu1 = t_end-s1[index1] + # edge correction + nu1 = max(t_end-s1[N1-1], s1[N1-1]-s1[N1-2]) if index2 < N2-1: nu2 = s2[index2+1]-s2[index2] else: - nu2 = t_end-s2[index2] + # edge correction + nu2 = max(t_end-s2[N2-1], s2[N2-1]-s2[N2-2]) # compute the corresponding isi-distance isi_values[index] = abs(nu1 - nu2) / \ max(nu1, nu2) diff --git a/test/test_distance.py b/test/test_distance.py index 20b52e8..19da35f 100644 --- a/test/test_distance.py +++ b/test/test_distance.py @@ -47,7 +47,7 @@ def test_isi(): t2 = SpikeTrain([0.1, 0.4, 0.5, 0.6], [0.0, 1.0]) expected_times = [0.0, 0.1, 0.2, 0.4, 0.5, 0.6, 1.0] - expected_isi = [0.1/0.2, 0.1/0.3, 0.1/0.3, 0.1/0.2, 0.1/0.2, 0.0/0.5] + expected_isi = [0.1/0.3, 0.1/0.3, 0.1/0.3, 0.1/0.2, 0.1/0.2, 0.0/0.5] expected_times = np.array(expected_times) expected_isi = np.array(expected_isi) @@ -332,7 +332,9 @@ def test_regression_spiky(): st2 = SpikeTrain(np.arange(100, 1201, 110), 1300) isi_dist = spk.isi_distance(st1, st2) - assert_almost_equal(isi_dist, 7.6923076923076941e-02, decimal=15) + assert_almost_equal(isi_dist, 9.0909090909090939e-02, decimal=15) + isi_profile = spk.isi_profile(st1, st2) + assert_equal(isi_profile.y, 0.1/1.1 * np.ones_like(isi_profile.y)) spike_dist = spk.spike_distance(st1, st2) assert_equal(spike_dist, 2.1105878248735391e-01) @@ -346,7 +348,7 @@ def test_regression_spiky(): (0.0, 4000.0)) isi_dist = spk.isi_distance_multi(spike_trains) # get the full precision from SPIKY - assert_almost_equal(isi_dist, 1.8318789829845508e-01, decimal=15) + assert_almost_equal(isi_dist, 0.17051816816999129656, decimal=15) spike_profile = spk.spike_profile_multi(spike_trains) assert_equal(len(spike_profile.y1)+len(spike_profile.y2), 1252) -- cgit v1.2.3