From a0262fc04e4b084f4dd270a75938d4ad029783d4 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Thu, 7 May 2015 23:08:12 +0200 Subject: performance improvements use recursive approach to compute average profile for average multivariate distances, dont compute average multivariate profile, but average distances directly. --- pyspike/spike_distance.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) (limited to 'pyspike/spike_distance.py') diff --git a/pyspike/spike_distance.py b/pyspike/spike_distance.py index cc620d4..75b3b0e 100644 --- a/pyspike/spike_distance.py +++ b/pyspike/spike_distance.py @@ -3,7 +3,8 @@ # Distributed under the BSD License from pyspike import PieceWiseLinFunc -from pyspike.generic import _generic_profile_multi, _generic_distance_matrix +from pyspike.generic import _generic_profile_multi, _generic_distance_multi, \ + _generic_distance_matrix ############################################################ @@ -117,7 +118,8 @@ def spike_distance_multi(spike_trains, indices=None, interval=None): :returns: The averaged multi-variate spike distance :math:`D_S`. :rtype: double """ - return spike_profile_multi(spike_trains, indices).avrg(interval) + return _generic_distance_multi(spike_trains, spike_distance, indices, + interval) ############################################################ -- 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_distance.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_distance.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 f688dc2e8616f914040746de845646abb158125d Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Fri, 8 May 2015 18:06:59 +0200 Subject: introduce backend for distance function isi- and spike distances over complete intervals are now computed without obtaining the profile first. This gives more than x2 performance improvements. --- examples/performance.py | 8 +- pyspike/cython/cython_distances.pyx | 328 ++++++++++++++++++++++++++++++++++++ pyspike/cython/cython_profiles.pyx | 2 +- pyspike/isi_distance.py | 16 +- pyspike/spike_distance.py | 17 +- setup.py | 5 +- 6 files changed, 371 insertions(+), 5 deletions(-) create mode 100644 pyspike/cython/cython_distances.pyx (limited to 'pyspike/spike_distance.py') diff --git a/examples/performance.py b/examples/performance.py index 94dae25..1c31e8f 100644 --- a/examples/performance.py +++ b/examples/performance.py @@ -1,4 +1,10 @@ -""" Compute distances of large sets of spike trains for performance tests +""" +Compute distances of large sets of spike trains for performance tests + +Copyright 2015, Mario Mulansky + +Distributed under the BSD License + """ from __future__ import print_function diff --git a/pyspike/cython/cython_distances.pyx b/pyspike/cython/cython_distances.pyx new file mode 100644 index 0000000..65c2872 --- /dev/null +++ b/pyspike/cython/cython_distances.pyx @@ -0,0 +1,328 @@ +#cython: boundscheck=False +#cython: wraparound=False +#cython: cdivision=True + +""" +cython_distances.pyx + +cython implementation of the isi-, spike- and spike-sync distances + +Note: using cython memoryviews (e.g. double[:]) instead of ndarray objects +improves the performance of spike_distance by a factor of 10! + +Copyright 2015, 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_distances.pyx + +which gives:: + + cython_distances.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 isi_value + cdef int index1, index2, index + cdef int N1, N2 + cdef double nu1, nu2 + cdef double last_t, curr_t, curr_isi + isi_value = 0.0 + N1 = len(s1) + N2 = len(s2) + + # 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 + + last_t = t_start + curr_isi = 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 + curr_t = 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 + curr_t = 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 + curr_t = 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_value += curr_isi * (curr_t - last_t) + curr_isi = fabs(nu1 - nu2) / fmax(nu1, nu2) + last_t = curr_t + index += 1 + + isi_value += curr_isi * (t_end - last_t) + # end nogil + + return isi_value / (t_end-t_start) + + +############################################################ +# 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 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 + cdef double y_start, y_end, t_last, t_current, spike_value + + spike_value = 0.0 + + N1 = len(t1) + N2 = len(t2) + + with nogil: # release the interpreter to allow multithreading + t_last = 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_start = (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 + t_curr = t_p1 + s2 = (dt_p2*(t_f2-t_p1) + dt_f2*(t_p1-t_p2)) / isi2 + y_end = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2) + + spike_value += 0.5*(y_start + y_end) * (t_curr - t_last) + + # 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_start = (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 + t_curr = t_p2 + s1 = (dt_p1*(t_f1-t_p2) + dt_f1*(t_p2-t_p1)) / isi1 + y_end = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2) + + spike_value += 0.5*(y_start + y_end) * (t_curr - t_last) + + # 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 + # s1 is the same as above, thus we can compute y2 immediately + y_start = (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 + t_curr = t_f1 + y_end = 0.0 + spike_value += 0.5*(y_start + y_end) * (t_curr - t_last) + y_start = 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 + t_last = t_curr + # 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_end = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2) + spike_value += 0.5*(y_start + y_end) * (t_end - t_last) + # end nogil + + # use only the data added above + # could be less than original length due to equal spike times + return spike_value / (t_end-t_start) diff --git a/pyspike/cython/cython_profiles.pyx b/pyspike/cython/cython_profiles.pyx index 59a8d30..3690091 100644 --- a/pyspike/cython/cython_profiles.pyx +++ b/pyspike/cython/cython_profiles.pyx @@ -10,7 +10,7 @@ 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 +Copyright 2014-2015, Mario Mulansky Distributed under the BSD License diff --git a/pyspike/isi_distance.py b/pyspike/isi_distance.py index 164378d..2a1ed3a 100644 --- a/pyspike/isi_distance.py +++ b/pyspike/isi_distance.py @@ -66,7 +66,21 @@ def isi_distance(spike_train1, spike_train2, interval=None): :returns: The isi-distance :math:`D_I`. :rtype: double """ - return isi_profile(spike_train1, spike_train2).avrg(interval) + + if interval is None: + # distance over the whole interval is requested: use specific function + # for optimal performance + try: + from cython.cython_distances import isi_distance_cython \ + as isi_distance_impl + return isi_distance_impl(spike_train1.spikes, spike_train2.spikes, + spike_train1.t_start, spike_train1.t_end) + except ImportError: + # Cython backend not available: fall back to profile averaging + return isi_profile(spike_train1, spike_train2).avrg(interval) + else: + # some specific interval is provided: use profile + return isi_profile(spike_train1, spike_train2).avrg(interval) ############################################################ diff --git a/pyspike/spike_distance.py b/pyspike/spike_distance.py index 3567585..d727fa2 100644 --- a/pyspike/spike_distance.py +++ b/pyspike/spike_distance.py @@ -68,7 +68,22 @@ def spike_distance(spike_train1, spike_train2, interval=None): :rtype: double """ - return spike_profile(spike_train1, spike_train2).avrg(interval) + if interval is None: + # distance over the whole interval is requested: use specific function + # for optimal performance + try: + from cython.cython_distances import spike_distance_cython \ + as spike_distance_impl + return spike_distance_impl(spike_train1.spikes, + spike_train2.spikes, + spike_train1.t_start, + spike_train1.t_end) + except ImportError: + # Cython backend not available: fall back to average profile + return spike_profile(spike_train1, spike_train2).avrg(interval) + else: + # some specific interval is provided: compute the whole profile + return spike_profile(spike_train1, spike_train2).avrg(interval) ############################################################ diff --git a/setup.py b/setup.py index d687240..7902066 100644 --- a/setup.py +++ b/setup.py @@ -22,7 +22,8 @@ else: use_cython = True if os.path.isfile("pyspike/cython/cython_add.c") and \ - os.path.isfile("pyspike/cython/cython_profiles.c"): + os.path.isfile("pyspike/cython/cython_profiles.c") and \ + os.path.isfile("pyspike/cython/cython_distances.c"): use_c = True else: use_c = False @@ -34,12 +35,14 @@ 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_profiles", ["pyspike/cython/cython_profiles.pyx"]), + Extension("pyspike.cython.cython_distances", ["pyspike/cython/cython_distances.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_profiles", ["pyspike/cython/cython_profiles.c"]), + Extension("pyspike.cython.cython_distances", ["pyspike/cython/cython_distances.c"]), ] # neither cython nor c files available -> automatic fall-back to python backend -- 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 'pyspike/spike_distance.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