From 6c68690b992a6dfabf8259e265cf427fc6f193eb Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Mon, 12 Oct 2015 11:20:03 +0200 Subject: added coincidence single profile to python backend missing function in python backend added, identified and fixed a bug in the implementation as well --- pyspike/cython/cython_profiles.pyx | 2 +- pyspike/cython/python_backend.py | 67 ++++++++++++++++++++++++++++---------- pyspike/spike_sync.py | 2 +- test/test_sync_filter.py | 3 +- 4 files changed, 53 insertions(+), 21 deletions(-) diff --git a/pyspike/cython/cython_profiles.pyx b/pyspike/cython/cython_profiles.pyx index eb4d157..aa24db4 100644 --- a/pyspike/cython/cython_profiles.pyx +++ b/pyspike/cython/cython_profiles.pyx @@ -473,7 +473,7 @@ def coincidence_single_profile_cython(double[:] spikes1, double[:] spikes2, if j > -1 and fabs(spikes1[i]-spikes2[j]) < tau: # current spike in st1 is coincident c[i] = 1 - if j < N2-1 and spikes2[j] < spikes1[i]: + if j < N2-1 and (j < 0 or spikes2[j] < spikes1[i]): # in case spikes2[j] is before spikes1[i] it has to be the one # right before (see above), hence we move one forward and also # check the next spike diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py index 6b7209a..a4a0c9e 100644 --- a/pyspike/cython/python_backend.py +++ b/pyspike/cython/python_backend.py @@ -3,7 +3,7 @@ Collection of python functions that can be used instead of the cython implementation. -Copyright 2014, Mario Mulansky +Copyright 2014-2015, Mario Mulansky Distributed under the BSD License @@ -356,26 +356,27 @@ def cumulative_sync_python(spikes1, spikes2): return st, c +def get_tau(spikes1, spikes2, i, j, max_tau, init_tau): + m = init_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: + m = min(m, spikes2[j+1]-spikes2[j]) + if i > 0: + m = min(m, spikes1[i]-spikes1[i-1]) + if j > 0: + m = min(m, spikes2[j]-spikes2[j-1]) + m *= 0.5 + if max_tau > 0.0: + m = min(m, max_tau) + return m + + ############################################################ # coincidence_python ############################################################ def coincidence_python(spikes1, spikes2, t_start, t_end, max_tau): - def get_tau(spikes1, spikes2, i, j, max_tau): - 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: - m = min(m, spikes2[j+1]-spikes2[j]) - if i > 0: - m = min(m, spikes1[i]-spikes1[i-1]) - if j > 0: - m = min(m, spikes2[j]-spikes2[j-1]) - m *= 0.5 - if max_tau > 0.0: - m = min(m, max_tau) - return m - N1 = len(spikes1) N2 = len(spikes2) i = -1 @@ -388,7 +389,7 @@ def coincidence_python(spikes1, spikes2, t_start, t_end, max_tau): 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, max_tau, t_end-t_start) st[n] = spikes1[i] if j > -1 and spikes1[i]-spikes2[j] < tau: # coincidence between the current spike and the previous spike @@ -398,7 +399,7 @@ def coincidence_python(spikes1, spikes2, t_start, t_end, max_tau): 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, max_tau, t_end-t_start) st[n] = spikes2[j] if i > -1 and spikes2[j]-spikes1[i] < tau: # coincidence between the current spike and the previous spike @@ -433,6 +434,36 @@ def coincidence_python(spikes1, spikes2, t_start, t_end, max_tau): return st, c, mp +############################################################ +# coincidence_single_profile_cython +############################################################ +def coincidence_single_python(spikes1, spikes2, t_start, t_end, max_tau): + + N1 = len(spikes1) + N2 = len(spikes2) + j = -1 + c = np.zeros(N1) # coincidences + for i in xrange(N1): + while j < N2-1 and spikes2[j+1] < spikes1[i]: + # move forward until spikes2[j] is the last spike before spikes1[i] + # note that if spikes2[j] is after spikes1[i] we dont do anything + j += 1 + tau = get_tau(spikes1, spikes2, i, j, max_tau, t_end-t_start) + if j > -1 and abs(spikes1[i]-spikes2[j]) < tau: + # current spike in st1 is coincident + c[i] = 1 + if j < N2-1 and (j < 0 or spikes2[j] < spikes1[i]): + # in case spikes2[j] is before spikes1[i] it has to be the first or + # the one right before (see above), hence we move one forward and + # also check the next spike + j += 1 + tau = get_tau(spikes1, spikes2, i, j, max_tau, t_end-t_start) + if abs(spikes2[j]-spikes1[i]) < tau: + # current spike in st1 is coincident + c[i] = 1 + return c + + ############################################################ # add_piece_wise_const_python ############################################################ diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py index 1d2ecdb..66aa906 100644 --- a/pyspike/spike_sync.py +++ b/pyspike/spike_sync.py @@ -315,7 +315,7 @@ sure that PySpike is installed by running\n \ 'python setup.py build_ext --inplace'!\n \ Falling back to slow python backend.") # use python backend - from cython.python_backend import coincidence_single_profile_python \ + from cython.python_backend import coincidence_single_python \ as coincidence_impl if max_tau is None: diff --git a/test/test_sync_filter.py b/test/test_sync_filter.py index 66ffcb6..e259903 100644 --- a/test/test_sync_filter.py +++ b/test/test_sync_filter.py @@ -28,12 +28,13 @@ def test_single_prof(): coincidence_single_profile_cython as coincidence_impl except ImportError: from pyspike.cython.python_backend import \ - coincidence_single_profile_python as coincidence_impl + coincidence_single_python as coincidence_impl sync_prof = spk.spike_sync_profile(SpikeTrain(st1, 5.0), SpikeTrain(st2, 5.0)) coincidences = np.array(coincidence_impl(st1, st2, 0, 5.0, 0.0)) + print(coincidences) for i, t in enumerate(st1): assert_equal(coincidences[i], sync_prof.y[sync_prof.x == t], "At index %d" % i) -- cgit v1.2.3