From 76a4bbcc733bdd24bb61072a341c43a14b7f83d1 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Fri, 8 May 2015 11:57:00 +0200 Subject: performance improvement for multivar spike sync dont compute the average profile in the function spike_sync_multi, but rather compute the overall average distance directly --- pyspike/spike_sync.py | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) (limited to 'pyspike/spike_sync.py') diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py index 9d2e363..0c78228 100644 --- a/pyspike/spike_sync.py +++ b/pyspike/spike_sync.py @@ -3,6 +3,7 @@ # Copyright 2014-2015, Mario Mulansky # Distributed under the BSD License +import numpy as np from functools import partial from pyspike import DiscreteFunc from pyspike.generic import _generic_profile_multi, _generic_distance_matrix @@ -131,8 +132,25 @@ def spike_sync_multi(spike_trains, indices=None, interval=None, max_tau=None): :rtype: double """ - return spike_sync_profile_multi(spike_trains, indices, - max_tau).avrg(interval) + if indices is None: + indices = np.arange(len(spike_trains)) + indices = np.array(indices) + # check validity of indices + assert (indices < len(spike_trains)).all() and (indices >= 0).all(), \ + "Invalid index list." + # generate a list of possible index pairs + pairs = [(indices[i], j) for i in range(len(indices)) + for j in indices[i+1:]] + + coincidence = 0.0 + mp = 0.0 + for (i, j) in pairs: + profile = spike_sync_profile(spike_trains[i], spike_trains[j]) + summed_vals = profile.integral(interval) + coincidence += summed_vals[0] + mp += summed_vals[1] + + return coincidence/mp ############################################################ -- cgit v1.2.3 From f44d78a5f0b4ab25bb443accdcb9fc1bd8ff57da Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Fri, 8 May 2015 12:08:48 +0200 Subject: improved error message --- pyspike/isi_distance.py | 4 ++-- pyspike/spike_distance.py | 4 ++-- pyspike/spike_sync.py | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) (limited to 'pyspike/spike_sync.py') diff --git a/pyspike/isi_distance.py b/pyspike/isi_distance.py index 21df561..8b1e9bf 100644 --- a/pyspike/isi_distance.py +++ b/pyspike/isi_distance.py @@ -25,9 +25,9 @@ def isi_profile(spike_train1, spike_train2): """ # check whether the spike trains are defined for the same interval assert spike_train1.t_start == spike_train2.t_start, \ - "Given spike trains seems not to have auxiliary spikes!" + "Given spike trains are not defined on the same interval!" assert spike_train1.t_end == spike_train2.t_end, \ - "Given spike trains seems not to have auxiliary spikes!" + "Given spike trains are not defined on the same interval!" # load cython implementation try: diff --git a/pyspike/spike_distance.py b/pyspike/spike_distance.py index 75b3b0e..499ab77 100644 --- a/pyspike/spike_distance.py +++ b/pyspike/spike_distance.py @@ -25,9 +25,9 @@ def spike_profile(spike_train1, spike_train2): """ # check whether the spike trains are defined for the same interval assert spike_train1.t_start == spike_train2.t_start, \ - "Given spike trains seems not to have auxiliary spikes!" + "Given spike trains are not defined on the same interval!" assert spike_train1.t_end == spike_train2.t_end, \ - "Given spike trains seems not to have auxiliary spikes!" + "Given spike trains are not defined on the same interval!" # cython implementation try: diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py index 0c78228..7d429e4 100644 --- a/pyspike/spike_sync.py +++ b/pyspike/spike_sync.py @@ -30,9 +30,9 @@ def spike_sync_profile(spike_train1, spike_train2, max_tau=None): """ # check whether the spike trains are defined for the same interval assert spike_train1.t_start == spike_train2.t_start, \ - "Given spike trains seems not to have auxiliary spikes!" + "Given spike trains are not defined on the same interval!" assert spike_train1.t_end == spike_train2.t_end, \ - "Given spike trains seems not to have auxiliary spikes!" + "Given spike trains are not defined on the same interval!" # cython implementation try: -- cgit v1.2.3 From 619ffd7105203938a26075c79a77d63960da9922 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Fri, 8 May 2015 12:29:47 +0200 Subject: renamed cython_distance module -> cython_profiles --- pyspike/cython/cython_distance.pyx | 423 ------------------------------------- pyspike/cython/cython_profiles.pyx | 423 +++++++++++++++++++++++++++++++++++++ pyspike/isi_distance.py | 10 +- pyspike/spike_distance.py | 16 +- pyspike/spike_sync.py | 15 +- setup.py | 8 +- 6 files changed, 447 insertions(+), 448 deletions(-) delete mode 100644 pyspike/cython/cython_distance.pyx create mode 100644 pyspike/cython/cython_profiles.pyx (limited to 'pyspike/spike_sync.py') diff --git a/pyspike/cython/cython_distance.pyx b/pyspike/cython/cython_distance.pyx deleted file mode 100644 index 6ee0181..0000000 --- a/pyspike/cython/cython_distance.pyx +++ /dev/null @@ -1,423 +0,0 @@ -#cython: boundscheck=False -#cython: wraparound=False -#cython: cdivision=True - -""" -cython_distances.pyx - -cython implementation of the isi- and spike-distance - -Note: using cython memoryviews (e.g. double[:]) instead of ndarray objects -improves the performance of spike_distance by a factor of 10! - -Copyright 2014, Mario Mulansky - -Distributed under the BSD License - -""" - -""" -To test whether things can be optimized: remove all yellow stuff -in the html output:: - - cython -a cython_distance.pyx - -which gives:: - - cython_distance.html - -""" - -import numpy as np -cimport numpy as np - -from libc.math cimport fabs -from libc.math cimport fmax -from libc.math cimport fmin - -DTYPE = np.float -ctypedef np.float_t DTYPE_t - - -############################################################ -# isi_distance_cython -############################################################ -def isi_distance_cython(double[:] s1, double[:] s2, - double t_start, double t_end): - - cdef double[:] spike_events - cdef double[:] isi_values - cdef int index1, index2, index - cdef int N1, N2 - cdef double nu1, nu2 - N1 = len(s1) - N2 = len(s2) - - spike_events = np.empty(N1+N2+2) - # the values have one entry less as they are defined at the intervals - isi_values = np.empty(N1+N2+1) - - # first x-value of the profile - spike_events[0] = t_start - - # first interspike interval - check if a spike exists at the start time - if 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: - # edge correction - nu2 = fmax(s2[0]-t_start, s2[1]-s2[0]) - index2 = -1 - else: - nu2 = s2[1]-s2[0] - index2 = 0 - - isi_values[0] = fabs(nu1-nu2)/fmax(nu1, nu2) - index = 1 - - with nogil: # release the interpreter to allow multithreading - while index1+index2 < N1+N2-2: - # check which spike is next, only if there are spikes left in 1 - # next spike in 1 is earlier, or there are no spikes left in 2 - if (index1 < N1-1) and ((index2 == N2-1) or - (s1[index1+1] < s2[index2+1])): - index1 += 1 - spike_events[index] = s1[index1] - if index1 < N1-1: - nu1 = s1[index1+1]-s1[index1] - else: - # 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 - spike_events[index] = s2[index2] - if index2 < N2-1: - nu2 = s2[index2+1]-s2[index2] - else: - # edge correction - nu2 = fmax(t_end-s2[index2], nu2) - else: # s1[index1+1] == s2[index2+1] - index1 += 1 - index2 += 1 - spike_events[index] = s1[index1] - if index1 < N1-1: - nu1 = s1[index1+1]-s1[index1] - else: - # edge correction - nu1 = fmax(t_end-s1[index1], nu1) - if index2 < N2-1: - nu2 = s2[index2+1]-s2[index2] - else: - # 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 - # the last event is the interval end - if spike_events[index-1] == t_end: - index -= 1 - else: - spike_events[index] = t_end - # end nogil - - return spike_events[:index+1], isi_values[:index] - - -############################################################ -# get_min_dist_cython -############################################################ -cdef inline double get_min_dist_cython(double spike_time, - double[:] spike_train, - # use memory view to ensure inlining - # np.ndarray[DTYPE_t,ndim=1] spike_train, - int N, - int start_index, - double t_start, double t_end) nogil: - """ Returns the minimal distance |spike_time - spike_train[i]| - with i>=start_index. - """ - cdef double d, d_temp - # start with the distance to the start time - d = fabs(spike_time - t_start) - if start_index < 0: - start_index = 0 - while start_index < N: - d_temp = fabs(spike_time - spike_train[start_index]) - if d_temp > d: - return d - else: - d = d_temp - start_index += 1 - - # finally, check the distance to end time - d_temp = fabs(t_end - spike_time) - if d_temp > d: - return d - else: - return d_temp - - -############################################################ -# isi_avrg_cython -############################################################ -cdef inline double isi_avrg_cython(double isi1, double isi2) nogil: - return 0.5*(isi1+isi2)*(isi1+isi2) - # alternative definition to obtain ~ 0.5 for Poisson spikes - # return 0.5*(isi1*isi1+isi2*isi2) - - -############################################################ -# spike_distance_cython -############################################################ -def spike_distance_cython(double[:] t1, double[:] t2, - double t_start, double t_end): - - cdef double[:] spike_events - cdef double[:] y_starts - cdef double[:] y_ends - - cdef int N1, N2, index1, index2, index - cdef double t_p1, t_f1, t_p2, t_f2, dt_p1, dt_p2, dt_f1, dt_f2 - cdef double isi1, isi2, s1, s2 - - N1 = len(t1) - N2 = len(t2) - - spike_events = np.empty(N1+N2+2) - - y_starts = np.empty(len(spike_events)-1) - y_ends = np.empty(len(spike_events)-1) - - with nogil: # release the interpreter to allow multithreading - spike_events[0] = t_start - t_p1 = t_start - t_p2 = t_start - if t1[0] > t_start: - # dt_p1 = t2[0]-t_start - t_f1 = t1[0] - dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_start, t_end) - isi1 = fmax(t_f1-t_start, t1[1]-t1[0]) - dt_p1 = dt_f1 - s1 = dt_p1*(t_f1-t_start)/isi1 - index1 = -1 - else: - t_f1 = t1[1] - dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_start, t_end) - dt_p1 = 0.0 - isi1 = t1[1]-t1[0] - s1 = dt_p1 - index1 = 0 - if t2[0] > t_start: - # dt_p1 = t2[0]-t_start - t_f2 = t2[0] - dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_start, t_end) - dt_p2 = dt_f2 - isi2 = fmax(t_f2-t_start, t2[1]-t2[0]) - s2 = dt_p2*(t_f2-t_start)/isi2 - index2 = -1 - else: - t_f2 = t2[1] - dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_start, t_end) - dt_p2 = 0.0 - isi2 = t2[1]-t2[0] - s2 = dt_p2 - index2 = 0 - - y_starts[0] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2) - index = 1 - - while index1+index2 < N1+N2-2: - # print(index, index1, index2) - if (index1 < N1-1) and (t_f1 < t_f2 or index2 == N2-1): - index1 += 1 - # first calculate the previous interval end value - s1 = dt_f1*(t_f1-t_p1) / isi1 - # the previous time now was the following time before: - dt_p1 = dt_f1 - t_p1 = t_f1 # t_p1 contains the current time point - # get the next time - if index1 < N1-1: - t_f1 = t1[index1+1] - else: - t_f1 = t_end - spike_events[index] = t_p1 - s2 = (dt_p2*(t_f2-t_p1) + dt_f2*(t_p1-t_p2)) / isi2 - y_ends[index-1] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, - isi2) - # now the next interval start value - if index1 < N1-1: - dt_f1 = get_min_dist_cython(t_f1, t2, N2, index2, - t_start, t_end) - isi1 = t_f1-t_p1 - s1 = dt_p1 - else: - dt_f1 = dt_p1 - isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2]) - # s1 needs adjustment due to change of isi1 - s1 = dt_p1*(t_end-t1[N1-1])/isi1 - # s2 is the same as above, thus we can compute y2 immediately - y_starts[index] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, - isi2) - elif (index2 < N2-1) and (t_f1 > t_f2 or index1 == N1-1): - index2 += 1 - # first calculate the previous interval end value - s2 = dt_f2*(t_f2-t_p2) / isi2 - # the previous time now was the following time before: - dt_p2 = dt_f2 - t_p2 = t_f2 # t_p2 contains the current time point - # get the next time - if index2 < N2-1: - t_f2 = t2[index2+1] - else: - t_f2 = t_end - spike_events[index] = t_p2 - s1 = (dt_p1*(t_f1-t_p2) + dt_f1*(t_p2-t_p1)) / isi1 - y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, - isi2) - # now the next interval start value - if index2 < N2-1: - dt_f2 = get_min_dist_cython(t_f2, t1, N1, index1, - t_start, t_end) - isi2 = t_f2-t_p2 - s2 = dt_p2 - else: - dt_f2 = dt_p2 - isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2]) - # s2 needs adjustment due to change of isi2 - s2 = dt_p2*(t_end-t2[N2-1])/isi2 - # s2 is the same as above, thus we can compute y2 immediately - y_starts[index] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2) - else: # t_f1 == t_f2 - generate only one event - index1 += 1 - index2 += 1 - t_p1 = t_f1 - t_p2 = t_f2 - dt_p1 = 0.0 - dt_p2 = 0.0 - spike_events[index] = t_f1 - y_ends[index-1] = 0.0 - y_starts[index] = 0.0 - if index1 < N1-1: - t_f1 = t1[index1+1] - dt_f1 = get_min_dist_cython(t_f1, t2, N2, index2, - t_start, t_end) - isi1 = t_f1 - t_p1 - else: - t_f1 = t_end - dt_f1 = dt_p1 - isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2]) - if index2 < N2-1: - t_f2 = t2[index2+1] - dt_f2 = get_min_dist_cython(t_f2, t1, N1, index1, - t_start, t_end) - isi2 = t_f2 - t_p2 - else: - t_f2 = t_end - dt_f2 = dt_p2 - isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2]) - index += 1 - # the last event is the interval end - if spike_events[index-1] == t_end: - index -= 1 - else: - spike_events[index] = t_end - # the ending value of the last interval - isi1 = max(t_end-t1[N1-1], t1[N1-1]-t1[N1-2]) - isi2 = max(t_end-t2[N2-1], t2[N2-1]-t2[N2-2]) - s1 = dt_f1*(t_end-t1[N1-1])/isi1 - s2 = dt_f2*(t_end-t2[N2-1])/isi2 - y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2) - # end nogil - - # use only the data added above - # could be less than original length due to equal spike times - return spike_events[:index+1], y_starts[:index], y_ends[:index] - - - -############################################################ -# coincidence_python -############################################################ -cdef inline double get_tau(double[:] spikes1, double[:] spikes2, - int i, int j, double max_tau): - cdef double m = 1E100 # some huge number - 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: - m = fmin(m, spikes1[i+1]-spikes1[i]) - if j < N2 and j > -1: - m = fmin(m, spikes2[j+1]-spikes2[j]) - if i > 0: - m = fmin(m, spikes1[i]-spikes1[i-1]) - if j > 0: - m = fmin(m, spikes2[j]-spikes2[j-1]) - m *= 0.5 - if max_tau > 0.0: - m = fmin(m, max_tau) - return m - - -############################################################ -# coincidence_cython -############################################################ -def coincidence_cython(double[:] spikes1, double[:] spikes2, - double t_start, double t_end, double max_tau): - - cdef int N1 = len(spikes1) - cdef int N2 = len(spikes2) - cdef int i = -1 - cdef int j = -1 - cdef int n = 0 - 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 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) - st[n] = spikes1[i] - if j > -1 and spikes1[i]-spikes2[j] < tau: - # coincidence between the current spike and the previous spike - # both get marked with 1 - c[n] = 1 - c[n-1] = 1 - 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) - st[n] = spikes2[j] - if i > -1 and spikes2[j]-spikes1[i] < tau: - # coincidence between the current spike and the previous spike - # both get marked with 1 - c[n] = 1 - c[n-1] = 1 - else: # spikes1[i+1] = spikes2[j+1] - # advance in both spike trains - j += 1 - i += 1 - n += 1 - # add only one event, but with coincidence 2 and multiplicity 2 - st[n] = spikes1[i] - c[n] = 2 - mp[n] = 2 - - st = st[:n+2] - c = c[:n+2] - mp = mp[:n+2] - - 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] - - return st, c, mp diff --git a/pyspike/cython/cython_profiles.pyx b/pyspike/cython/cython_profiles.pyx new file mode 100644 index 0000000..59a8d30 --- /dev/null +++ b/pyspike/cython/cython_profiles.pyx @@ -0,0 +1,423 @@ +#cython: boundscheck=False +#cython: wraparound=False +#cython: cdivision=True + +""" +cython_profiles.pyx + +cython implementation of the isi-, spike- and spike-sync profiles + +Note: using cython memoryviews (e.g. double[:]) instead of ndarray objects +improves the performance of spike_distance by a factor of 10! + +Copyright 2014, Mario Mulansky + +Distributed under the BSD License + +""" + +""" +To test whether things can be optimized: remove all yellow stuff +in the html output:: + + cython -a cython_profiles.pyx + +which gives:: + + cython_profiles.html + +""" + +import numpy as np +cimport numpy as np + +from libc.math cimport fabs +from libc.math cimport fmax +from libc.math cimport fmin + +DTYPE = np.float +ctypedef np.float_t DTYPE_t + + +############################################################ +# isi_profile_cython +############################################################ +def isi_profile_cython(double[:] s1, double[:] s2, + double t_start, double t_end): + + cdef double[:] spike_events + cdef double[:] isi_values + cdef int index1, index2, index + cdef int N1, N2 + cdef double nu1, nu2 + N1 = len(s1) + N2 = len(s2) + + spike_events = np.empty(N1+N2+2) + # the values have one entry less as they are defined at the intervals + isi_values = np.empty(N1+N2+1) + + # first x-value of the profile + spike_events[0] = t_start + + # first interspike interval - check if a spike exists at the start time + if 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: + # edge correction + nu2 = fmax(s2[0]-t_start, s2[1]-s2[0]) + index2 = -1 + else: + nu2 = s2[1]-s2[0] + index2 = 0 + + isi_values[0] = fabs(nu1-nu2)/fmax(nu1, nu2) + index = 1 + + with nogil: # release the interpreter to allow multithreading + while index1+index2 < N1+N2-2: + # check which spike is next, only if there are spikes left in 1 + # next spike in 1 is earlier, or there are no spikes left in 2 + if (index1 < N1-1) and ((index2 == N2-1) or + (s1[index1+1] < s2[index2+1])): + index1 += 1 + spike_events[index] = s1[index1] + if index1 < N1-1: + nu1 = s1[index1+1]-s1[index1] + else: + # 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 + spike_events[index] = s2[index2] + if index2 < N2-1: + nu2 = s2[index2+1]-s2[index2] + else: + # edge correction + nu2 = fmax(t_end-s2[index2], nu2) + else: # s1[index1+1] == s2[index2+1] + index1 += 1 + index2 += 1 + spike_events[index] = s1[index1] + if index1 < N1-1: + nu1 = s1[index1+1]-s1[index1] + else: + # edge correction + nu1 = fmax(t_end-s1[index1], nu1) + if index2 < N2-1: + nu2 = s2[index2+1]-s2[index2] + else: + # 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 + # the last event is the interval end + if spike_events[index-1] == t_end: + index -= 1 + else: + spike_events[index] = t_end + # end nogil + + return spike_events[:index+1], isi_values[:index] + + +############################################################ +# get_min_dist_cython +############################################################ +cdef inline double get_min_dist_cython(double spike_time, + double[:] spike_train, + # use memory view to ensure inlining + # np.ndarray[DTYPE_t,ndim=1] spike_train, + int N, + int start_index, + double t_start, double t_end) nogil: + """ Returns the minimal distance |spike_time - spike_train[i]| + with i>=start_index. + """ + cdef double d, d_temp + # start with the distance to the start time + d = fabs(spike_time - t_start) + if start_index < 0: + start_index = 0 + while start_index < N: + d_temp = fabs(spike_time - spike_train[start_index]) + if d_temp > d: + return d + else: + d = d_temp + start_index += 1 + + # finally, check the distance to end time + d_temp = fabs(t_end - spike_time) + if d_temp > d: + return d + else: + return d_temp + + +############################################################ +# isi_avrg_cython +############################################################ +cdef inline double isi_avrg_cython(double isi1, double isi2) nogil: + return 0.5*(isi1+isi2)*(isi1+isi2) + # alternative definition to obtain ~ 0.5 for Poisson spikes + # return 0.5*(isi1*isi1+isi2*isi2) + + +############################################################ +# spike_profile_cython +############################################################ +def spike_profile_cython(double[:] t1, double[:] t2, + double t_start, double t_end): + + cdef double[:] spike_events + cdef double[:] y_starts + cdef double[:] y_ends + + cdef int N1, N2, index1, index2, index + cdef double t_p1, t_f1, t_p2, t_f2, dt_p1, dt_p2, dt_f1, dt_f2 + cdef double isi1, isi2, s1, s2 + + N1 = len(t1) + N2 = len(t2) + + spike_events = np.empty(N1+N2+2) + + y_starts = np.empty(len(spike_events)-1) + y_ends = np.empty(len(spike_events)-1) + + with nogil: # release the interpreter to allow multithreading + spike_events[0] = t_start + t_p1 = t_start + t_p2 = t_start + if t1[0] > t_start: + # dt_p1 = t2[0]-t_start + t_f1 = t1[0] + dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_start, t_end) + isi1 = fmax(t_f1-t_start, t1[1]-t1[0]) + dt_p1 = dt_f1 + s1 = dt_p1*(t_f1-t_start)/isi1 + index1 = -1 + else: + t_f1 = t1[1] + dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_start, t_end) + dt_p1 = 0.0 + isi1 = t1[1]-t1[0] + s1 = dt_p1 + index1 = 0 + if t2[0] > t_start: + # dt_p1 = t2[0]-t_start + t_f2 = t2[0] + dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_start, t_end) + dt_p2 = dt_f2 + isi2 = fmax(t_f2-t_start, t2[1]-t2[0]) + s2 = dt_p2*(t_f2-t_start)/isi2 + index2 = -1 + else: + t_f2 = t2[1] + dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_start, t_end) + dt_p2 = 0.0 + isi2 = t2[1]-t2[0] + s2 = dt_p2 + index2 = 0 + + y_starts[0] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2) + index = 1 + + while index1+index2 < N1+N2-2: + # print(index, index1, index2) + if (index1 < N1-1) and (t_f1 < t_f2 or index2 == N2-1): + index1 += 1 + # first calculate the previous interval end value + s1 = dt_f1*(t_f1-t_p1) / isi1 + # the previous time now was the following time before: + dt_p1 = dt_f1 + t_p1 = t_f1 # t_p1 contains the current time point + # get the next time + if index1 < N1-1: + t_f1 = t1[index1+1] + else: + t_f1 = t_end + spike_events[index] = t_p1 + s2 = (dt_p2*(t_f2-t_p1) + dt_f2*(t_p1-t_p2)) / isi2 + y_ends[index-1] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, + isi2) + # now the next interval start value + if index1 < N1-1: + dt_f1 = get_min_dist_cython(t_f1, t2, N2, index2, + t_start, t_end) + isi1 = t_f1-t_p1 + s1 = dt_p1 + else: + dt_f1 = dt_p1 + isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2]) + # s1 needs adjustment due to change of isi1 + s1 = dt_p1*(t_end-t1[N1-1])/isi1 + # s2 is the same as above, thus we can compute y2 immediately + y_starts[index] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, + isi2) + elif (index2 < N2-1) and (t_f1 > t_f2 or index1 == N1-1): + index2 += 1 + # first calculate the previous interval end value + s2 = dt_f2*(t_f2-t_p2) / isi2 + # the previous time now was the following time before: + dt_p2 = dt_f2 + t_p2 = t_f2 # t_p2 contains the current time point + # get the next time + if index2 < N2-1: + t_f2 = t2[index2+1] + else: + t_f2 = t_end + spike_events[index] = t_p2 + s1 = (dt_p1*(t_f1-t_p2) + dt_f1*(t_p2-t_p1)) / isi1 + y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, + isi2) + # now the next interval start value + if index2 < N2-1: + dt_f2 = get_min_dist_cython(t_f2, t1, N1, index1, + t_start, t_end) + isi2 = t_f2-t_p2 + s2 = dt_p2 + else: + dt_f2 = dt_p2 + isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2]) + # s2 needs adjustment due to change of isi2 + s2 = dt_p2*(t_end-t2[N2-1])/isi2 + # s2 is the same as above, thus we can compute y2 immediately + y_starts[index] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2) + else: # t_f1 == t_f2 - generate only one event + index1 += 1 + index2 += 1 + t_p1 = t_f1 + t_p2 = t_f2 + dt_p1 = 0.0 + dt_p2 = 0.0 + spike_events[index] = t_f1 + y_ends[index-1] = 0.0 + y_starts[index] = 0.0 + if index1 < N1-1: + t_f1 = t1[index1+1] + dt_f1 = get_min_dist_cython(t_f1, t2, N2, index2, + t_start, t_end) + isi1 = t_f1 - t_p1 + else: + t_f1 = t_end + dt_f1 = dt_p1 + isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2]) + if index2 < N2-1: + t_f2 = t2[index2+1] + dt_f2 = get_min_dist_cython(t_f2, t1, N1, index1, + t_start, t_end) + isi2 = t_f2 - t_p2 + else: + t_f2 = t_end + dt_f2 = dt_p2 + isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2]) + index += 1 + # the last event is the interval end + if spike_events[index-1] == t_end: + index -= 1 + else: + spike_events[index] = t_end + # the ending value of the last interval + isi1 = max(t_end-t1[N1-1], t1[N1-1]-t1[N1-2]) + isi2 = max(t_end-t2[N2-1], t2[N2-1]-t2[N2-2]) + s1 = dt_f1*(t_end-t1[N1-1])/isi1 + s2 = dt_f2*(t_end-t2[N2-1])/isi2 + y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2) + # end nogil + + # use only the data added above + # could be less than original length due to equal spike times + return spike_events[:index+1], y_starts[:index], y_ends[:index] + + + +############################################################ +# 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 + 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: + m = fmin(m, spikes1[i+1]-spikes1[i]) + if j < N2 and j > -1: + m = fmin(m, spikes2[j+1]-spikes2[j]) + if i > 0: + m = fmin(m, spikes1[i]-spikes1[i-1]) + if j > 0: + m = fmin(m, spikes2[j]-spikes2[j-1]) + m *= 0.5 + if max_tau > 0.0: + m = fmin(m, max_tau) + return m + + +############################################################ +# coincidence_profile_cython +############################################################ +def coincidence_profile_cython(double[:] spikes1, double[:] spikes2, + double t_start, double t_end, double max_tau): + + cdef int N1 = len(spikes1) + cdef int N2 = len(spikes2) + cdef int i = -1 + cdef int j = -1 + cdef int n = 0 + 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 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) + st[n] = spikes1[i] + if j > -1 and spikes1[i]-spikes2[j] < tau: + # coincidence between the current spike and the previous spike + # both get marked with 1 + c[n] = 1 + c[n-1] = 1 + 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) + st[n] = spikes2[j] + if i > -1 and spikes2[j]-spikes1[i] < tau: + # coincidence between the current spike and the previous spike + # both get marked with 1 + c[n] = 1 + c[n-1] = 1 + else: # spikes1[i+1] = spikes2[j+1] + # advance in both spike trains + j += 1 + i += 1 + n += 1 + # add only one event, but with coincidence 2 and multiplicity 2 + st[n] = spikes1[i] + c[n] = 2 + mp[n] = 2 + + st = st[:n+2] + c = c[:n+2] + mp = mp[:n+2] + + 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] + + return st, c, mp diff --git a/pyspike/isi_distance.py b/pyspike/isi_distance.py index 8b1e9bf..164378d 100644 --- a/pyspike/isi_distance.py +++ b/pyspike/isi_distance.py @@ -31,18 +31,18 @@ def isi_profile(spike_train1, spike_train2): # load cython implementation try: - from cython.cython_distance import isi_distance_cython \ - as isi_distance_impl + from cython.cython_profiles import isi_profile_cython \ + as isi_profile_impl except ImportError: print("Warning: isi_distance_cython not found. Make 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 isi_distance_python \ - as isi_distance_impl + as isi_profile_impl - times, values = isi_distance_impl(spike_train1.spikes, spike_train2.spikes, - spike_train1.t_start, spike_train1.t_end) + times, values = isi_profile_impl(spike_train1.spikes, spike_train2.spikes, + spike_train1.t_start, spike_train1.t_end) return PieceWiseConstFunc(times, values) diff --git a/pyspike/spike_distance.py b/pyspike/spike_distance.py index 499ab77..3567585 100644 --- a/pyspike/spike_distance.py +++ b/pyspike/spike_distance.py @@ -31,20 +31,20 @@ def spike_profile(spike_train1, spike_train2): # cython implementation try: - from cython.cython_distance import spike_distance_cython \ - as spike_distance_impl + from cython.cython_profiles import spike_profile_cython \ + as spike_profile_impl except ImportError: - print("Warning: spike_distance_cython not found. Make sure that \ + print("Warning: spike_profile_cython not found. Make 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 spike_distance_python \ - as spike_distance_impl + as spike_profile_impl - times, y_starts, y_ends = spike_distance_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.spikes, + spike_train2.spikes, + spike_train1.t_start, + spike_train1.t_end) return PieceWiseLinFunc(times, y_starts, y_ends) diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py index 7d429e4..107734d 100644 --- a/pyspike/spike_sync.py +++ b/pyspike/spike_sync.py @@ -36,24 +36,23 @@ def spike_sync_profile(spike_train1, spike_train2, max_tau=None): # cython implementation try: - from cython.cython_distance import coincidence_cython \ - as coincidence_impl + from cython.cython_profiles import coincidence_profile_cython \ + as coincidence_profile_impl except ImportError: print("Warning: spike_distance_cython not found. Make 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_python \ - as coincidence_impl + as coincidence_profile_impl if max_tau is None: max_tau = 0.0 - times, coincidences, multiplicity = coincidence_impl(spike_train1.spikes, - spike_train2.spikes, - spike_train1.t_start, - spike_train1.t_end, - max_tau) + times, coincidences, multiplicity \ + = coincidence_profile_impl(spike_train1.spikes, spike_train2.spikes, + spike_train1.t_start, spike_train1.t_end, + max_tau) return DiscreteFunc(times, coincidences, multiplicity) diff --git a/setup.py b/setup.py index 289d521..d687240 100644 --- a/setup.py +++ b/setup.py @@ -22,7 +22,7 @@ else: use_cython = True if os.path.isfile("pyspike/cython/cython_add.c") and \ - os.path.isfile("pyspike/cython/cython_distance.c"): + os.path.isfile("pyspike/cython/cython_profiles.c"): use_c = True else: use_c = False @@ -33,13 +33,13 @@ ext_modules = [] if use_cython: # Cython is available, compile .pyx -> .c ext_modules += [ Extension("pyspike.cython.cython_add", ["pyspike/cython/cython_add.pyx"]), - Extension("pyspike.cython.cython_distance", ["pyspike/cython/cython_distance.pyx"]), + Extension("pyspike.cython.cython_profiles", ["pyspike/cython/cython_profiles.pyx"]), ] cmdclass.update({'build_ext': build_ext}) elif use_c: # c files are there, compile to binaries ext_modules += [ Extension("pyspike.cython.cython_add", ["pyspike/cython/cython_add.c"]), - Extension("pyspike.cython.cython_distance", ["pyspike/cython/cython_distance.c"]), + Extension("pyspike.cython.cython_profiles", ["pyspike/cython/cython_profiles.c"]), ] # neither cython nor c files available -> automatic fall-back to python backend @@ -78,7 +78,7 @@ train similarity', 'Programming Language :: Python :: 2.7', ], package_data={ - 'pyspike': ['cython/cython_add.c', 'cython/cython_distance.c'], + 'pyspike': ['cython/cython_add.c', 'cython/cython_profiles.c'], 'test': ['Spike_testdata.txt'] } ) -- cgit v1.2.3 From ad7961c9f73d77b87ee23349169618d128418a51 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Mon, 11 May 2015 12:06:19 +0200 Subject: performance improvement for spike sync Additional cython implementation for overall spike sync values. It is not necessary to compute the profile anymore if only the spike sync value is required. 3x performance gain. --- pyspike/cython/cython_distances.pyx | 64 +++++++++++++++++++++++++++++++++++++ pyspike/spike_sync.py | 46 ++++++++++++++++++++++---- 2 files changed, 104 insertions(+), 6 deletions(-) (limited to 'pyspike/spike_sync.py') diff --git a/pyspike/cython/cython_distances.pyx b/pyspike/cython/cython_distances.pyx index 65c2872..bf90638 100644 --- a/pyspike/cython/cython_distances.pyx +++ b/pyspike/cython/cython_distances.pyx @@ -326,3 +326,67 @@ def spike_distance_cython(double[:] t1, double[:] t2, # use only the data added above # could be less than original length due to equal spike times return spike_value / (t_end-t_start) + + + +############################################################ +# 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 + 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: + m = fmin(m, spikes1[i+1]-spikes1[i]) + if j < N2 and j > -1: + m = fmin(m, spikes2[j+1]-spikes2[j]) + if i > 0: + m = fmin(m, spikes1[i]-spikes1[i-1]) + if j > 0: + m = fmin(m, spikes2[j]-spikes2[j-1]) + m *= 0.5 + if max_tau > 0.0: + m = fmin(m, max_tau) + return m + + +############################################################ +# coincidence_value_cython +############################################################ +def coincidence_value_cython(double[:] spikes1, double[:] spikes2, + double t_start, double t_end, double max_tau): + + cdef int N1 = len(spikes1) + cdef int N2 = len(spikes2) + cdef int i = -1 + cdef int j = -1 + cdef double coinc = 0.0 + cdef double mp = 0.0 + 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) + if j > -1 and spikes1[i]-spikes2[j] < tau: + # coincidence between the current spike and the previous spike + # both get marked with 1 + coinc += 2 + 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) + if i > -1 and spikes2[j]-spikes1[i] < tau: + # coincidence between the current spike and the previous spike + # both get marked with 1 + coinc += 2 + else: # spikes1[i+1] = spikes2[j+1] + # advance in both spike trains + j += 1 + i += 1 + # add only one event, but with coincidence 2 and multiplicity 2 + mp += 2 + coinc += 2 + + return coinc, mp diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py index 107734d..40d98d2 100644 --- a/pyspike/spike_sync.py +++ b/pyspike/spike_sync.py @@ -57,6 +57,40 @@ Falling back to slow python backend.") return DiscreteFunc(times, coincidences, multiplicity) +############################################################ +# _spike_sync_values +############################################################ +def _spike_sync_values(spike_train1, spike_train2, interval, max_tau): + """" Internal function. Computes the summed coincidences and multiplicity + for spike synchronization of the two given spike trains. + + Do not call this function directly, use `spike_sync` or `spike_sync_multi` + instead. + """ + if interval is None: + # distance over the whole interval is requested: use specific function + # for optimal performance + try: + from cython.cython_distances import coincidence_value_cython \ + as coincidence_value_impl + if max_tau is None: + max_tau = 0.0 + c, mp = coincidence_value_impl(spike_train1.spikes, + spike_train2.spikes, + spike_train1.t_start, + spike_train1.t_end, + max_tau) + return c, mp + except ImportError: + # Cython backend not available: fall back to profile averaging + return spike_sync_profile(spike_train1, spike_train2, + max_tau).integral(interval) + else: + # some specific interval is provided: use profile + return spike_sync_profile(spike_train1, spike_train2, + max_tau).integral(interval) + + ############################################################ # spike_sync ############################################################ @@ -80,8 +114,8 @@ def spike_sync(spike_train1, spike_train2, interval=None, max_tau=None): :rtype: `double` """ - return spike_sync_profile(spike_train1, spike_train2, - max_tau).avrg(interval) + c, mp = _spike_sync_values(spike_train1, spike_train2, interval, max_tau) + return 1.0*c/mp ############################################################ @@ -144,10 +178,10 @@ def spike_sync_multi(spike_trains, indices=None, interval=None, max_tau=None): coincidence = 0.0 mp = 0.0 for (i, j) in pairs: - profile = spike_sync_profile(spike_trains[i], spike_trains[j]) - summed_vals = profile.integral(interval) - coincidence += summed_vals[0] - mp += summed_vals[1] + c, m = _spike_sync_values(spike_trains[i], spike_trains[j], + interval, max_tau) + coincidence += c + mp += m return coincidence/mp -- cgit v1.2.3