From 27aa30d1fdb830a04b608c702cf7b616115eeb50 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Wed, 22 Apr 2015 18:18:30 +0200 Subject: added SpikeTrain class, changed isi_distance spike trains are now represented as SpikeTrain objects consisting of the spike times and the interval edges. The implementation of the ISI-distance has been modified accordingly. The SPIKE-distance and SPIKE-Synchronization are still to be done. --- pyspike/cython/cython_distance.pyx | 85 ++++++++++++++++++++++++-------------- 1 file changed, 55 insertions(+), 30 deletions(-) (limited to 'pyspike/cython/cython_distance.pyx') diff --git a/pyspike/cython/cython_distance.pyx b/pyspike/cython/cython_distance.pyx index 2834ca5..1d652ee 100644 --- a/pyspike/cython/cython_distance.pyx +++ b/pyspike/cython/cython_distance.pyx @@ -42,57 +42,82 @@ ctypedef np.float_t DTYPE_t ############################################################ # isi_distance_cython ############################################################ -def isi_distance_cython(double[:] s1, - double[:] s2): +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)-1 - N2 = len(s2)-1 + 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: + nu1 = s1[0] - t_start + index1 = -1 + else: + nu1 = s1[1]-s1[0] + index1 = 0 + + if s2[0] > t_start: + nu2 = s2[0] - t_start + index2 = -1 + else: + nu2 = s2[1]-s2[0] + index2 = 0 - nu1 = s1[1]-s1[0] - nu2 = s2[1]-s2[0] - spike_events = np.empty(N1+N2) - spike_events[0] = s1[0] - # the values have one entry less - the number of intervals between events - isi_values = np.empty(N1+N2-1) + isi_values[0] = fabs(nu1-nu2)/fmax(nu1, nu2) + index = 1 with nogil: # release the interpreter to allow multithreading - isi_values[0] = fabs(nu1-nu2)/fmax(nu1, nu2) - index1 = 0 - index2 = 0 - index = 1 - while True: - # check which spike is next - from s1 or s2 - if s1[index1+1] < s2[index2+1]: + 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 - # break condition relies on existence of spikes at T_end - if index1 >= N1: - break spike_events[index] = s1[index1] - nu1 = s1[index1+1]-s1[index1] - elif s1[index1+1] > s2[index2+1]: + if index1 < N1-1: + nu1 = s1[index1+1]-s1[index1] + else: + nu1 = t_end-s1[index1] + elif (index2 < N2-1) and ((index1 == N1-1) or + (s1[index1+1] > s2[index2+1])): index2 += 1 - if index2 >= N2: - break spike_events[index] = s2[index2] - nu2 = s2[index2+1]-s2[index2] + if index2 < N2-1: + nu2 = s2[index2+1]-s2[index2] + else: + nu2 = t_end-s2[index2] else: # s1[index1+1] == s2[index2+1] index1 += 1 index2 += 1 - if (index1 >= N1) or (index2 >= N2): - break spike_events[index] = s1[index1] - nu1 = s1[index1+1]-s1[index1] - nu2 = s2[index2+1]-s2[index2] + if index1 < N1-1: + nu1 = s1[index1+1]-s1[index1] + else: + nu1 = t_end-s1[index1] + if index2 < N2-1: + nu2 = s2[index2+1]-s2[index2] + else: + nu2 = t_end-s2[index2] # compute the corresponding isi-distance isi_values[index] = fabs(nu1 - nu2) / fmax(nu1, nu2) index += 1 # the last event is the interval end - spike_events[index] = s1[N1] + 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] -- cgit v1.2.3 From ed85a9b72edcb7bba6ae1105e213b3b0a2f78d3a Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Fri, 24 Apr 2015 00:49:16 +0200 Subject: changed spike distance to use new SpikeTrain class --- pyspike/cython/cython_distance.pyx | 199 +++++++++++++++++++++++++----------- pyspike/cython/python_backend.py | 203 ++++++++++++++++++++++++------------- pyspike/spike_distance.py | 37 ++++--- test/test_distance.py | 34 ++++--- 4 files changed, 313 insertions(+), 160 deletions(-) (limited to 'pyspike/cython/cython_distance.pyx') diff --git a/pyspike/cython/cython_distance.pyx b/pyspike/cython/cython_distance.pyx index 1d652ee..7999e0a 100644 --- a/pyspike/cython/cython_distance.pyx +++ b/pyspike/cython/cython_distance.pyx @@ -131,21 +131,30 @@ cdef inline double get_min_dist_cython(double spike_time, # use memory view to ensure inlining # np.ndarray[DTYPE_t,ndim=1] spike_train, int N, - int start_index=0) nogil: + 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 - d = fabs(spike_time - spike_train[start_index]) - start_index += 1 + # 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: - break + return d else: d = d_temp start_index += 1 - return d + + # finally, check the distance to end time + d_temp = fabs(t_end - spike_time) + if d_temp > d: + return d + else: + return d_temp ############################################################ @@ -160,96 +169,162 @@ cdef inline double isi_avrg_cython(double isi1, double isi2) nogil: ############################################################ # spike_distance_cython ############################################################ -def spike_distance_cython(double[:] t1, - double[:] t2): +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 dt_p1, dt_p2, dt_f1, dt_f2, isi1, isi2, s1, s2 + 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) - spike_events[0] = t1[0] + 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 - index1 = 0 - index2 = 0 - index = 1 - dt_p1 = 0.0 - dt_f1 = get_min_dist_cython(t1[1], t2, N2, 0) - dt_p2 = 0.0 - dt_f2 = get_min_dist_cython(t2[1], t1, N1, 0) - isi1 = max(t1[1]-t1[0], t1[2]-t1[1]) - isi2 = max(t2[1]-t2[0], t2[2]-t2[1]) - s1 = dt_f1*(t1[1]-t1[0])/isi1 - s2 = dt_f2*(t2[1]-t2[0])/isi2 + spike_events[0] = t_start + t_p1 = t_start + t_p2 = t_start + if t1[0] > t_start: + # dt_p1 = t2[0]-t_start + dt_p1 = 0.0 + 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]) + s1 = dt_f1*(t_f1-t_start)/isi1 + index1 = -1 + else: + dt_p1 = 0.0 + t_f1 = t1[1] + dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_start, t_end) + isi1 = t1[1]-t1[0] + s1 = dt_p1 + index1 = 0 + if t2[0] > t_start: + # dt_p1 = t2[0]-t_start + dt_p2 = 0.0 + t_f2 = t2[0] + dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_start, t_end) + isi2 = fmax(t_f2-t_start, t2[1]-t2[0]) + s2 = dt_f2*(t_f2-t_start)/isi2 + index2 = -1 + else: + dt_p2 = 0.0 + t_f2 = t2[1] + dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_start, t_end) + isi2 = t2[1]-t2[0] + s2 = dt_p2 + index2 = 0 + y_starts[0] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2) - while True: + index = 1 + + while index1+index2 < N1+N2-2: # print(index, index1, index2) - if t1[index1+1] < t2[index2+1]: + if (index1 < N1-1) and (t_f1 < t_f2 or index2 == N2-1): index1 += 1 - # break condition relies on existence of spikes at T_end - if index1+1 >= N1: - break - spike_events[index] = t1[index1] # first calculate the previous interval end value - dt_p1 = dt_f1 # the previous time now was the following time before + # 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 s1 = dt_p1 - s2 = (dt_p2*(t2[index2+1]-t1[index1]) + - dt_f2*(t1[index1]-t2[index2])) / isi2 - y_ends[index-1] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2) + 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 - dt_f1 = get_min_dist_cython(t1[index1+1], t2, N2, index2) - isi1 = t1[index1+1]-t1[index1] + if index1 < N1-1: + dt_f1 = get_min_dist_cython(t_f1, t2, N2, index2, + t_start, t_end) + isi1 = t_f1-t_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 t1[index1+1] > t2[index2+1]: + 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 - if index2+1 >= N2: - break - spike_events[index] = t2[index2] # first calculate the previous interval end value - dt_p2 = dt_f2 # the previous time now was the following time before - s1 = (dt_p1*(t1[index1+1]-t2[index2]) + - dt_f1*(t2[index2]-t1[index1])) / isi1 + # the previous time now was the following time before: + dt_p2 = dt_f2 + t_p2 = t_f2 # t_p1 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 s2 = dt_p2 - y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2) + y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, + isi2) # now the next interval start value - dt_f2 = get_min_dist_cython(t2[index2+1], t1, N1, index1) - #s2 = dt_f2 - isi2 = t2[index2+1]-t2[index2] + if index2 < N2-1: + dt_f2 = get_min_dist_cython(t_f2, t1, N1, index1, + t_start, t_end) + isi2 = t_f2-t_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: # t1[index1+1] == t2[index2+1] - generate only one event + else: # t_f1 == t_f2 - generate only one event index1 += 1 index2 += 1 - if (index1+1 >= N1) or (index2+1 >= N2): - break - spike_events[index] = t1[index1] - y_ends[index-1] = 0.0 - y_starts[index] = 0.0 + t_p1 = t_f1 + t_p2 = t_f2 dt_p1 = 0.0 dt_p2 = 0.0 - dt_f1 = get_min_dist_cython(t1[index1+1], t2, N2, index2) - dt_f2 = get_min_dist_cython(t2[index2+1], t1, N1, index1) - isi1 = t1[index1+1]-t1[index1] - isi2 = t2[index2+1]-t2[index2] + 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 - spike_events[index] = t1[N1-1] - # the ending value of the last interval - isi1 = max(t1[N1-1]-t1[N1-2], t1[N1-2]-t1[N1-3]) - isi2 = max(t2[N2-1]-t2[N2-2], t2[N2-2]-t2[N2-3]) - s1 = dt_p1*(t1[N1-1]-t1[N1-2])/isi1 - s2 = dt_p2*(t2[N2-1]-t2[N2-2])/isi2 - y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2) + 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 diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py index 4c37236..bcf9c30 100644 --- a/pyspike/cython/python_backend.py +++ b/pyspike/cython/python_backend.py @@ -89,122 +89,185 @@ def isi_distance_python(s1, s2, t_start, t_end): ############################################################ # get_min_dist ############################################################ -def get_min_dist(spike_time, spike_train, start_index=0): +def get_min_dist(spike_time, spike_train, start_index, t_start, t_end): """ Returns the minimal distance |spike_time - spike_train[i]| with i>=start_index. """ - d = abs(spike_time - spike_train[start_index]) - start_index += 1 + d = abs(spike_time - t_start) + if start_index < 0: + start_index = 0 while start_index < len(spike_train): d_temp = abs(spike_time - spike_train[start_index]) if d_temp > d: - break + return d else: d = d_temp start_index += 1 - return d + # finally, check the distance to end time + d_temp = abs(t_end - spike_time) + if d_temp > d: + return d + else: + return d_temp ############################################################ # spike_distance_python ############################################################ -def spike_distance_python(spikes1, spikes2): +def spike_distance_python(spikes1, spikes2, t_start, t_end): """ Computes the instantaneous spike-distance S_spike (t) of the two given spike trains. The spike trains are expected to have auxiliary spikes at the beginning and end of the interval. Use the function add_auxiliary_spikes to add those spikes to the spike train. Args: - spikes1, spikes2: ordered arrays of spike times with auxiliary spikes. + - t_start, t_end: edges of the spike train Returns: - PieceWiseLinFunc describing the spike-distance. """ - # check for auxiliary spikes - first and last spikes should be identical - assert spikes1[0] == spikes2[0], \ - "Given spike trains seems not to have auxiliary spikes!" - assert spikes1[-1] == spikes2[-1], \ - "Given spike trains seems not to have auxiliary spikes!" + # shorter variables t1 = spikes1 t2 = spikes2 - spike_events = np.empty(len(t1) + len(t2) - 2) - spike_events[0] = t1[0] - y_starts = np.empty(len(spike_events) - 1) - y_ends = np.empty(len(spike_events) - 1) + N1 = len(t1) + N2 = len(t2) - index1 = 0 - index2 = 0 + spike_events = np.empty(N1+N2+2) + + y_starts = np.empty(len(spike_events)-1) + y_ends = np.empty(len(spike_events)-1) + + spike_events[0] = t_start + t_p1 = t_start + t_p2 = t_start + if t1[0] > t_start: + # dt_p1 = t2[0]-t_start + dt_p1 = 0.0 + t_f1 = t1[0] + dt_f1 = get_min_dist(t_f1, t2, 0, t_start, t_end) + isi1 = max(t_f1-t_start, t1[1]-t1[0]) + s1 = dt_f1*(t_f1-t_start)/isi1 + index1 = -1 + else: + dt_p1 = 0.0 + t_f1 = t1[1] + dt_f1 = get_min_dist(t_f1, t2, 0, t_start, t_end) + isi1 = t1[1]-t1[0] + s1 = dt_p1 + index1 = 0 + if t2[0] > t_start: + # dt_p1 = t2[0]-t_start + dt_p2 = 0.0 + t_f2 = t2[0] + dt_f2 = get_min_dist(t_f2, t1, 0, t_start, t_end) + isi2 = max(t_f2-t_start, t2[1]-t2[0]) + s2 = dt_f2*(t_f2-t_start)/isi2 + index2 = -1 + else: + dt_p2 = 0.0 + t_f2 = t2[1] + dt_f2 = get_min_dist(t_f2, t1, 0, t_start, t_end) + isi2 = t2[1]-t2[0] + s2 = dt_p2 + index2 = 0 + + y_starts[0] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2) index = 1 - dt_p1 = 0.0 - dt_f1 = get_min_dist(t1[1], t2, 0) - dt_p2 = 0.0 - dt_f2 = get_min_dist(t2[1], t1, 0) - isi1 = max(t1[1]-t1[0], t1[2]-t1[1]) - isi2 = max(t2[1]-t2[0], t2[2]-t2[1]) - s1 = dt_f1*(t1[1]-t1[0])/isi1 - s2 = dt_f2*(t2[1]-t2[0])/isi2 - y_starts[0] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2) - while True: + + while index1+index2 < N1+N2-2: # print(index, index1, index2) - if t1[index1+1] < t2[index2+1]: + if (index1 < N1-1) and (t_f1 < t_f2 or index2 == N2-1): index1 += 1 - # break condition relies on existence of spikes at T_end - if index1+1 >= len(t1): - break - spike_events[index] = t1[index1] # first calculate the previous interval end value - dt_p1 = dt_f1 # the previous time was the following time before + # 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 s1 = dt_p1 - s2 = (dt_p2*(t2[index2+1]-t1[index1]) + - dt_f2*(t1[index1]-t2[index2])) / isi2 - y_ends[index-1] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2) + s2 = (dt_p2*(t_f2-t_p1) + dt_f2*(t_p1-t_p2)) / isi2 + y_ends[index-1] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2) # now the next interval start value - dt_f1 = get_min_dist(t1[index1+1], t2, index2) - isi1 = t1[index1+1]-t1[index1] + if index1 < N1-1: + dt_f1 = get_min_dist(t_f1, t2, index2, t_start, t_end) + isi1 = t_f1-t_p1 + else: + dt_f1 = dt_p1 + isi1 = max(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) / ((isi1+isi2)**2/2) - elif t1[index1+1] > t2[index2+1]: + y_starts[index] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2) + elif (index2 < N2-1) and (t_f1 > t_f2 or index1 == N1-1): index2 += 1 - if index2+1 >= len(t2): - break - spike_events[index] = t2[index2] # first calculate the previous interval end value - dt_p2 = dt_f2 # the previous time was the following time before - s1 = (dt_p1*(t1[index1+1]-t2[index2]) + - dt_f1*(t2[index2]-t1[index1])) / isi1 + # the previous time now was the following time before: + dt_p2 = dt_f2 + t_p2 = t_f2 # t_p1 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 s2 = dt_p2 - y_ends[index-1] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2) + y_ends[index-1] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2) # now the next interval start value - dt_f2 = get_min_dist(t2[index2+1], t1, index1) - #s2 = dt_f2 - isi2 = t2[index2+1]-t2[index2] + if index2 < N2-1: + dt_f2 = get_min_dist(t_f2, t1, index1, t_start, t_end) + isi2 = t_f2-t_p2 + else: + dt_f2 = dt_p2 + isi2 = max(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) / ((isi1+isi2)**2/2) - else: # t1[index1+1] == t2[index2+1] - generate only one event + y_starts[index] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2) + else: # t_f1 == t_f2 - generate only one event index1 += 1 index2 += 1 - if (index1+1 >= len(t1)) or (index2+1 >= len(t2)): - break - assert dt_f2 == 0.0 - assert dt_f1 == 0.0 - spike_events[index] = t1[index1] - y_ends[index-1] = 0.0 - y_starts[index] = 0.0 + t_p1 = t_f1 + t_p2 = t_f2 dt_p1 = 0.0 dt_p2 = 0.0 - dt_f1 = get_min_dist(t1[index1+1], t2, index2) - dt_f2 = get_min_dist(t2[index2+1], t1, index1) - isi1 = t1[index1+1]-t1[index1] - isi2 = t2[index2+1]-t2[index2] + 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(t_f1, t2, index2, t_start, t_end) + isi1 = t_f1 - t_p1 + else: + t_f1 = t_end + dt_f1 = dt_p1 + isi1 = max(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(t_f2, t1, index1, t_start, t_end) + isi2 = t_f2 - t_p2 + else: + t_f2 = t_end + dt_f2 = dt_p2 + isi2 = max(t_end-t2[N2-1], t2[N2-1]-t2[N2-2]) index += 1 # the last event is the interval end - spike_events[index] = t1[-1] - # the ending value of the last interval - isi1 = max(t1[-1]-t1[-2], t1[-2]-t1[-3]) - isi2 = max(t2[-1]-t2[-2], t2[-2]-t2[-3]) - s1 = dt_p1*(t1[-1]-t1[-2])/isi1 - s2 = dt_p2*(t2[-1]-t2[-2])/isi2 - y_ends[index-1] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2) + 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) / (0.5*(isi1+isi2)**2) + # 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] diff --git a/pyspike/spike_distance.py b/pyspike/spike_distance.py index f721c86..8d03d70 100644 --- a/pyspike/spike_distance.py +++ b/pyspike/spike_distance.py @@ -14,23 +14,23 @@ from pyspike.generic import _generic_profile_multi, _generic_distance_matrix ############################################################ # spike_profile ############################################################ -def spike_profile(spikes1, spikes2): +def spike_profile(spike_train1, spike_train2): """ Computes the spike-distance profile S_spike(t) of the two given spike trains. Returns the profile as a PieceWiseLinFunc object. The S_spike - values are defined positive S_spike(t)>=0. The spike trains are expected to - have auxiliary spikes at the beginning and end of the interval. Use the - function add_auxiliary_spikes to add those spikes to the spike train. + values are defined positive S_spike(t)>=0. - :param spikes1: ordered array of spike times with auxiliary spikes. - :param spikes2: ordered array of spike times with auxiliary spikes. + :param spike_train1: First spike train. + :type spike_train1: :class:`pyspike.SpikeTrain` + :param spike_train2: Second spike train. + :type spike_train2: :class:`pyspike.SpikeTrain` :returns: The spike-distance profile :math:`S_{spike}(t)`. :rtype: :class:`pyspike.function.PieceWiseLinFunc` """ - # check for auxiliary spikes - first and last spikes should be identical - assert spikes1[0] == spikes2[0], \ + # 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!" - assert spikes1[-1] == spikes2[-1], \ + assert spike_train1.t_end == spike_train2.t_end, \ "Given spike trains seems not to have auxiliary spikes!" # cython implementation @@ -45,21 +45,26 @@ Falling back to slow python backend.") from cython.python_backend import spike_distance_python \ as spike_distance_impl - times, y_starts, y_ends = spike_distance_impl(spikes1, spikes2) + times, y_starts, y_ends = spike_distance_impl(spike_train1.spikes, + spike_train2.spikes, + spike_train1.t_start, + spike_train1.t_end) return PieceWiseLinFunc(times, y_starts, y_ends) ############################################################ # spike_distance ############################################################ -def spike_distance(spikes1, spikes2, interval=None): +def spike_distance(spike_train1, spike_train2, interval=None): """ Computes the spike-distance S of the given spike trains. The spike-distance is the integral over the isi distance profile S_spike(t): .. math:: S = \int_{T_0}^{T_1} S_{spike}(t) dt. - :param spikes1: ordered array of spike times with auxiliary spikes. - :param spikes2: ordered array of spike times with auxiliary spikes. + :param spike_train1: First spike train. + :type spike_train1: :class:`pyspike.SpikeTrain` + :param spike_train2: Second spike train. + :type spike_train2: :class:`pyspike.SpikeTrain` :param interval: averaging interval given as a pair of floats (T0, T1), if None the average over the whole function is computed. :type interval: Pair of floats or None. @@ -67,7 +72,7 @@ def spike_distance(spikes1, spikes2, interval=None): :rtype: double """ - return spike_profile(spikes1, spikes2).avrg(interval) + return spike_profile(spike_train1, spike_train2).avrg(interval) ############################################################ @@ -102,7 +107,7 @@ def spike_distance_multi(spike_trains, indices=None, interval=None): S_{spike} = \int_0^T 2/((N(N-1)) sum_{} S_{spike}^{i, j} dt where the sum goes over all pairs - :param spike_trains: list of spike trains + :param spike_trains: list of :class:`pyspike.SpikeTrain` :param indices: list of indices defining which spike trains to use, if None all given spike trains are used (default=None) :type indices: list or None @@ -121,7 +126,7 @@ def spike_distance_multi(spike_trains, indices=None, interval=None): def spike_distance_matrix(spike_trains, indices=None, interval=None): """ Computes the time averaged spike-distance of all pairs of spike-trains. - :param spike_trains: list of spike trains + :param spike_trains: list of :class:`pyspike.SpikeTrain` :param indices: list of indices defining which spike trains to use, if None all given spike trains are used (default=None) :type indices: list or None diff --git a/test/test_distance.py b/test/test_distance.py index b54e908..4af0e63 100644 --- a/test/test_distance.py +++ b/test/test_distance.py @@ -64,13 +64,26 @@ def test_isi(): def test_spike(): # generate two spike trains: - t1 = np.array([0.2, 0.4, 0.6, 0.7]) - t2 = np.array([0.3, 0.45, 0.8, 0.9, 0.95]) + t1 = SpikeTrain([0.0, 2.0, 5.0, 8.0], 10.0) + t2 = SpikeTrain([0.0, 1.0, 5.0, 9.0], 10.0) + + expected_times = np.array([0.0, 1.0, 2.0, 5.0, 8.0, 9.0, 10.0]) + + f = spk.spike_profile(t1, t2) + + assert_equal(f.x, expected_times) + + assert_almost_equal(f.avrg(), 0.1662415, decimal=6) + assert_almost_equal(f.y2[-1], 0.1394558, decimal=6) + + t1 = SpikeTrain([0.2, 0.4, 0.6, 0.7], 1.0) + t2 = SpikeTrain([0.3, 0.45, 0.8, 0.9, 0.95], 1.0) # pen&paper calculation of the spike distance expected_times = [0.0, 0.2, 0.3, 0.4, 0.45, 0.6, 0.7, 0.8, 0.9, 0.95, 1.0] s1 = np.array([0.1, 0.1, (0.1*0.1+0.05*0.1)/0.2, 0.05, (0.05*0.15 * 2)/0.2, - 0.15, 0.1, 0.1*0.2/0.3, 0.1**2/0.3, 0.1*0.05/0.3, 0.1]) + 0.15, 0.1, (0.1*0.1+0.1*0.2)/0.3, (0.1*0.2+0.1*0.1)/0.3, + (0.1*0.05+0.1*0.25)/0.3, 0.1]) s2 = np.array([0.1, 0.1*0.2/0.3, 0.1, (0.1*0.05 * 2)/.15, 0.05, (0.05*0.2+0.1*0.15)/0.35, (0.05*0.1+0.1*0.25)/0.35, 0.1, 0.1, 0.05, 0.05]) @@ -86,19 +99,18 @@ def test_spike(): (expected_y1+expected_y2)/2) expected_spike_val /= (expected_times[-1]-expected_times[0]) - t1 = spk.add_auxiliary_spikes(t1, 1.0) - t2 = spk.add_auxiliary_spikes(t2, 1.0) f = spk.spike_profile(t1, t2) assert_equal(f.x, expected_times) assert_array_almost_equal(f.y1, expected_y1, decimal=15) assert_array_almost_equal(f.y2, expected_y2, decimal=15) - assert_equal(f.avrg(), expected_spike_val) - assert_equal(spk.spike_distance(t1, t2), expected_spike_val) + assert_almost_equal(f.avrg(), expected_spike_val, decimal=15) + assert_almost_equal(spk.spike_distance(t1, t2), expected_spike_val, + decimal=15) # check with some equal spike times - t1 = np.array([0.2, 0.4, 0.6]) - t2 = np.array([0.1, 0.4, 0.5, 0.6]) + t1 = SpikeTrain([0.2, 0.4, 0.6], [0.0, 1.0]) + t2 = SpikeTrain([0.1, 0.4, 0.5, 0.6], [0.0, 1.0]) expected_times = [0.0, 0.1, 0.2, 0.4, 0.5, 0.6, 1.0] s1 = np.array([0.1, 0.1*0.1/0.2, 0.1, 0.0, 0.0, 0.0, 0.0]) @@ -115,8 +127,6 @@ def test_spike(): (expected_y1+expected_y2)/2) expected_spike_val /= (expected_times[-1]-expected_times[0]) - t1 = spk.add_auxiliary_spikes(t1, 1.0) - t2 = spk.add_auxiliary_spikes(t2, 1.0) f = spk.spike_profile(t1, t2) assert_equal(f.x, expected_times) @@ -315,6 +325,6 @@ def test_multi_variate_subsets(): if __name__ == "__main__": test_isi() - # test_spike() + test_spike() # test_multi_isi() # test_multi_spike() -- cgit v1.2.3 From 7da6da8533f9f76a99b959c9de37138377119ffc Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Fri, 24 Apr 2015 12:08:05 +0200 Subject: changed spike sync implementation to SpikeTrain --- pyspike/SpikeTrain.py | 8 ++--- pyspike/cython/cython_distance.pyx | 37 +++++++++++------------ pyspike/spike_sync.py | 37 +++++++++++++++-------- test/test_distance.py | 62 +++++++++++++++++++++++--------------- 4 files changed, 83 insertions(+), 61 deletions(-) (limited to 'pyspike/cython/cython_distance.pyx') diff --git a/pyspike/SpikeTrain.py b/pyspike/SpikeTrain.py index 4760014..041f897 100644 --- a/pyspike/SpikeTrain.py +++ b/pyspike/SpikeTrain.py @@ -21,14 +21,14 @@ class SpikeTrain: """ # TODO: sanity checks - self.spikes = np.array(spike_times) + self.spikes = np.array(spike_times, dtype=float) # check if interval is as sequence if not isinstance(interval, collections.Sequence): # treat value as end time and assume t_start = 0 self.t_start = 0.0 - self.t_end = interval + self.t_end = float(interval) else: # extract times from sequence - self.t_start = interval[0] - self.t_end = interval[1] + self.t_start = float(interval[0]) + self.t_end = float(interval[1]) diff --git a/pyspike/cython/cython_distance.pyx b/pyspike/cython/cython_distance.pyx index 7999e0a..6d998b9 100644 --- a/pyspike/cython/cython_distance.pyx +++ b/pyspike/cython/cython_distance.pyx @@ -339,11 +339,11 @@ def spike_distance_cython(double[:] t1, double[:] t2, cdef inline double get_tau(double[:] spikes1, double[:] spikes2, int i, int j, max_tau): cdef double m = 1E100 # some huge number - cdef int N1 = len(spikes1)-2 - cdef int N2 = len(spikes2)-2 - if i < N1: + cdef int N1 = len(spikes1)-1 + cdef int N2 = len(spikes2)-1 + if i < N1 and i > -1: m = fmin(m, spikes1[i+1]-spikes1[i]) - if j < N2: + if j < N2 and j > -1: m = fmin(m, spikes2[j+1]-spikes2[j]) if i > 1: m = fmin(m, spikes1[i]-spikes1[i-1]) @@ -358,34 +358,35 @@ cdef inline double get_tau(double[:] spikes1, double[:] spikes2, ############################################################ # coincidence_cython ############################################################ -def coincidence_cython(double[:] spikes1, double[:] spikes2, double max_tau): +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 = 0 - cdef int j = 0 + 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[:] 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 n < N1 + N2 - 2: - if spikes1[i+1] < spikes2[j+1]: + while i + j < N1 + N2 - 2: + if (i < N1-1) and (spikes1[i+1] < spikes2[j+1] or j == N2-1): i += 1 n += 1 tau = get_tau(spikes1, spikes2, i, j, max_tau) st[n] = spikes1[i] - if j > 0 and spikes1[i]-spikes2[j] < tau: + 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 spikes1[i+1] > spikes2[j+1]: + elif (j < N2-1) and (spikes1[i+1] > spikes2[j+1] or i == N1-1): j += 1 n += 1 tau = get_tau(spikes1, spikes2, i, j, max_tau) st[n] = spikes2[j] - if i > 0 and spikes2[j]-spikes1[i] < tau: + 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 @@ -394,8 +395,6 @@ def coincidence_cython(double[:] spikes1, double[:] spikes2, double max_tau): # advance in both spike trains j += 1 i += 1 - if i == N1-1 or j == N2-1: - break n += 1 # add only one event, but with coincidence 2 and multiplicity 2 st[n] = spikes1[i] @@ -406,8 +405,8 @@ def coincidence_cython(double[:] spikes1, double[:] spikes2, double max_tau): c = c[:n+2] mp = mp[:n+2] - st[0] = spikes1[0] - st[len(st)-1] = spikes1[len(spikes1)-1] + 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] diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py index e12ebb8..bca6f73 100644 --- a/pyspike/spike_sync.py +++ b/pyspike/spike_sync.py @@ -16,22 +16,27 @@ from pyspike.generic import _generic_profile_multi, _generic_distance_matrix ############################################################ # spike_sync_profile ############################################################ -def spike_sync_profile(spikes1, spikes2, max_tau=None): +def spike_sync_profile(spike_train1, spike_train2, max_tau=None): """ Computes the spike-synchronization profile S_sync(t) of the two given spike trains. Returns the profile as a DiscreteFunction object. The S_sync values are either 1 or 0, indicating the presence or absence of a - coincidence. The spike trains are expected to have auxiliary spikes at the - beginning and end of the interval. Use the function add_auxiliary_spikes to - add those spikes to the spike train. + coincidence. - :param spikes1: ordered array of spike times with auxiliary spikes. - :param spikes2: ordered array of spike times with auxiliary spikes. + :param spike_train1: First spike train. + :type spike_train1: :class:`pyspike.SpikeTrain` + :param spike_train2: Second spike train. + :type spike_train2: :class:`pyspike.SpikeTrain` :param max_tau: Maximum coincidence window size. If 0 or `None`, the coincidence window has no upper bound. :returns: The spike-distance profile :math:`S_{sync}(t)`. :rtype: :class:`pyspike.function.DiscreteFunction` """ + # 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!" + assert spike_train1.t_end == spike_train2.t_end, \ + "Given spike trains seems not to have auxiliary spikes!" # cython implementation try: @@ -48,7 +53,10 @@ Falling back to slow python backend.") if max_tau is None: max_tau = 0.0 - times, coincidences, multiplicity = coincidence_impl(spikes1, spikes2, + times, coincidences, multiplicity = coincidence_impl(spike_train1.spikes, + spike_train2.spikes, + spike_train1.t_start, + spike_train1.t_end, max_tau) return DiscreteFunc(times, coincidences, multiplicity) @@ -57,15 +65,17 @@ Falling back to slow python backend.") ############################################################ # spike_sync ############################################################ -def spike_sync(spikes1, spikes2, interval=None, max_tau=None): +def spike_sync(spike_train1, spike_train2, interval=None, max_tau=None): """ Computes the spike synchronization value SYNC of the given spike trains. The spike synchronization value is the computed as the total number of coincidences divided by the total number of spikes: .. math:: SYNC = \sum_n C_n / N. - :param spikes1: ordered array of spike times with auxiliary spikes. - :param spikes2: ordered array of spike times with auxiliary spikes. + :param spike_train1: First spike train. + :type spike_train1: :class:`pyspike.SpikeTrain` + :param spike_train2: Second spike train. + :type spike_train2: :class:`pyspike.SpikeTrain` :param interval: averaging interval given as a pair of floats (T0, T1), if None the average over the whole function is computed. :type interval: Pair of floats or None. @@ -74,7 +84,8 @@ def spike_sync(spikes1, spikes2, interval=None, max_tau=None): :returns: The spike synchronization value. :rtype: double """ - return spike_sync_profile(spikes1, spikes2, max_tau).avrg(interval) + return spike_sync_profile(spike_train1, spike_train2, + max_tau).avrg(interval) ############################################################ @@ -87,7 +98,7 @@ def spike_sync_profile_multi(spike_trains, indices=None, max_tau=None): spike trains pairs involving the spike train of containing this spike, which is the number of spike trains minus one (N-1). - :param spike_trains: list of spike trains + :param spike_trains: list of :class:`pyspike.SpikeTrain` :param indices: list of indices defining which spike trains to use, if None all given spike trains are used (default=None) :type indices: list or None @@ -134,7 +145,7 @@ def spike_sync_matrix(spike_trains, indices=None, interval=None, max_tau=None): """ Computes the overall spike-synchronization value of all pairs of spike-trains. - :param spike_trains: list of spike trains + :param spike_trains: list of :class:`pyspike.SpikeTrain` :param indices: list of indices defining which spike trains to use, if None all given spike trains are used (default=None) :type indices: list or None diff --git a/test/test_distance.py b/test/test_distance.py index 4af0e63..dbb72f1 100644 --- a/test/test_distance.py +++ b/test/test_distance.py @@ -138,10 +138,17 @@ def test_spike(): def test_spike_sync(): - spikes1 = np.array([1.0, 2.0, 3.0]) - spikes2 = np.array([2.1]) - spikes1 = spk.add_auxiliary_spikes(spikes1, 4.0) - spikes2 = spk.add_auxiliary_spikes(spikes2, 4.0) + spikes1 = SpikeTrain([1.0, 2.0, 3.0], 4.0) + spikes2 = SpikeTrain([2.1], 4.0) + + expected_x = np.array([0.0, 1.0, 2.0, 2.1, 3.0, 4.0]) + expected_y = np.array([0.0, 0.0, 1.0, 1.0, 0.0, 0.0]) + + f = spk.spike_sync_profile(spikes1, spikes2) + + assert_array_almost_equal(f.x, expected_x, decimal=16) + assert_array_almost_equal(f.y, expected_y, decimal=16) + assert_almost_equal(spk.spike_sync(spikes1, spikes2), 0.5, decimal=16) @@ -149,28 +156,34 @@ def test_spike_sync(): assert_almost_equal(spk.spike_sync(spikes1, spikes2, max_tau=0.05), 0.0, decimal=16) - spikes2 = np.array([3.1]) - spikes2 = spk.add_auxiliary_spikes(spikes2, 4.0) + spikes2 = SpikeTrain([3.1], 4.0) assert_almost_equal(spk.spike_sync(spikes1, spikes2), 0.5, decimal=16) - spikes2 = np.array([1.1]) - spikes2 = spk.add_auxiliary_spikes(spikes2, 4.0) + spikes2 = SpikeTrain([1.1], 4.0) + + expected_x = np.array([0.0, 1.0, 1.1, 2.0, 3.0, 4.0]) + expected_y = np.array([1.0, 1.0, 1.0, 0.0, 0.0, 0.0]) + + f = spk.spike_sync_profile(spikes1, spikes2) + + assert_array_almost_equal(f.x, expected_x, decimal=16) + assert_array_almost_equal(f.y, expected_y, decimal=16) + assert_almost_equal(spk.spike_sync(spikes1, spikes2), 0.5, decimal=16) - spikes2 = np.array([0.9]) - spikes2 = spk.add_auxiliary_spikes(spikes2, 4.0) + spikes2 = SpikeTrain([0.9], 4.0) assert_almost_equal(spk.spike_sync(spikes1, spikes2), 0.5, decimal=16) def check_multi_profile(profile_func, profile_func_multi): # generate spike trains: - t1 = spk.add_auxiliary_spikes(np.array([0.2, 0.4, 0.6, 0.7]), 1.0) - t2 = spk.add_auxiliary_spikes(np.array([0.3, 0.45, 0.8, 0.9, 0.95]), 1.0) - t3 = spk.add_auxiliary_spikes(np.array([0.2, 0.4, 0.6]), 1.0) - t4 = spk.add_auxiliary_spikes(np.array([0.1, 0.4, 0.5, 0.6]), 1.0) + t1 = SpikeTrain([0.2, 0.4, 0.6, 0.7], 1.0) + t2 = SpikeTrain([0.3, 0.45, 0.8, 0.9, 0.95], 1.0) + t3 = SpikeTrain([0.2, 0.4, 0.6], 1.0) + t4 = SpikeTrain([0.1, 0.4, 0.5, 0.6], 1.0) spike_trains = [t1, t2, t3, t4] f12 = profile_func(t1, t2) @@ -213,15 +226,12 @@ def test_multi_spike(): def test_multi_spike_sync(): # some basic multivariate check - spikes1 = np.array([100, 300, 400, 405, 410, 500, 700, 800, - 805, 810, 815, 900], dtype=float) - spikes2 = np.array([100, 200, 205, 210, 295, 350, 400, 510, - 600, 605, 700, 910], dtype=float) - spikes3 = np.array([100, 180, 198, 295, 412, 420, 510, 640, - 695, 795, 820, 920], dtype=float) - spikes1 = spk.add_auxiliary_spikes(spikes1, 1000) - spikes2 = spk.add_auxiliary_spikes(spikes2, 1000) - spikes3 = spk.add_auxiliary_spikes(spikes3, 1000) + spikes1 = SpikeTrain([100, 300, 400, 405, 410, 500, 700, 800, + 805, 810, 815, 900], 1000) + spikes2 = SpikeTrain([100, 200, 205, 210, 295, 350, 400, 510, + 600, 605, 700, 910], 1000) + spikes3 = SpikeTrain([100, 180, 198, 295, 412, 420, 510, 640, + 695, 795, 820, 920], 1000) assert_almost_equal(spk.spike_sync(spikes1, spikes2), 0.5, decimal=15) assert_almost_equal(spk.spike_sync(spikes1, spikes3), @@ -326,5 +336,7 @@ def test_multi_variate_subsets(): if __name__ == "__main__": test_isi() test_spike() - # test_multi_isi() - # test_multi_spike() + test_spike_sync() + test_multi_isi() + test_multi_spike() + test_multi_spike_sync() -- cgit v1.2.3 From 3bf9e12e6b5667fb1ea72c969848dacaff3cb470 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Fri, 24 Apr 2015 14:58:39 +0200 Subject: further adjustments in spike sync --- pyspike/SpikeTrain.py | 2 +- pyspike/__init__.py | 4 +- pyspike/cython/cython_distance.pyx | 6 +- pyspike/spike_sync.py | 6 +- pyspike/spikes.py | 113 ++++++++++++------------------------- test/test_distance.py | 5 +- 6 files changed, 49 insertions(+), 87 deletions(-) (limited to 'pyspike/cython/cython_distance.pyx') diff --git a/pyspike/SpikeTrain.py b/pyspike/SpikeTrain.py index 041f897..89520c9 100644 --- a/pyspike/SpikeTrain.py +++ b/pyspike/SpikeTrain.py @@ -9,7 +9,7 @@ import numpy as np import collections -class SpikeTrain: +class SpikeTrain(object): """ Class representing spike trains for the PySpike Module.""" def __init__(self, spike_times, interval): diff --git a/pyspike/__init__.py b/pyspike/__init__.py index 76e58a1..a5f9f0a 100644 --- a/pyspike/__init__.py +++ b/pyspike/__init__.py @@ -21,5 +21,5 @@ from spike_sync import spike_sync_profile, spike_sync,\ spike_sync_profile_multi, spike_sync_multi, spike_sync_matrix from psth import psth -from spikes import add_auxiliary_spikes, load_spike_trains_from_txt, \ - spike_train_from_string, merge_spike_trains, generate_poisson_spikes +from spikes import load_spike_trains_from_txt, spike_train_from_string, \ + merge_spike_trains, generate_poisson_spikes diff --git a/pyspike/cython/cython_distance.pyx b/pyspike/cython/cython_distance.pyx index 6d998b9..2841da8 100644 --- a/pyspike/cython/cython_distance.pyx +++ b/pyspike/cython/cython_distance.pyx @@ -337,10 +337,10 @@ def spike_distance_cython(double[:] t1, double[:] t2, # coincidence_python ############################################################ cdef inline double get_tau(double[:] spikes1, double[:] spikes2, - int i, int j, max_tau): + int i, int j, double max_tau): cdef double m = 1E100 # some huge number - cdef int N1 = len(spikes1)-1 - cdef int N2 = len(spikes2)-1 + 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: diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py index bca6f73..8ddd32c 100644 --- a/pyspike/spike_sync.py +++ b/pyspike/spike_sync.py @@ -109,10 +109,10 @@ def spike_sync_profile_multi(spike_trains, indices=None, max_tau=None): """ prof_func = partial(spike_sync_profile, max_tau=max_tau) - average_dist, M = _generic_profile_multi(spike_trains, prof_func, + average_prof, M = _generic_profile_multi(spike_trains, prof_func, indices) # average_dist.mul_scalar(1.0/M) # no normalization here! - return average_dist + return average_prof ############################################################ @@ -122,7 +122,7 @@ def spike_sync_multi(spike_trains, indices=None, interval=None, max_tau=None): """ Computes the multi-variate spike synchronization value for a set of spike trains. - :param spike_trains: list of spike trains + :param spike_trains: list of :class:`pyspike.SpikeTrain` :param indices: list of indices defining which spike trains to use, if None all given spike trains are used (default=None) :type indices: list or None diff --git a/pyspike/spikes.py b/pyspike/spikes.py index 9d7d6f4..128873d 100644 --- a/pyspike/spikes.py +++ b/pyspike/spikes.py @@ -8,82 +8,46 @@ Distributed under the BSD License """ import numpy as np - - -############################################################ -# add_auxiliary_spikes -############################################################ -def add_auxiliary_spikes(spike_train, time_interval): - """ Adds spikes at the beginning and end of the given time interval. - - :param spike_train: ordered array of spike times - :param time_interval: A pair (T_start, T_end) of values representing the - start and end time of the spike train measurement or - a single value representing the end time, the T_start - is then assuemd as 0. Auxiliary spikes will be added - to the spike train at the beginning and end of this - interval, if they are not yet present. - :type time_interval: pair of doubles or double - :returns: spike train with additional spikes at T_start and T_end. - - """ - try: - T_start = time_interval[0] - T_end = time_interval[1] - except: - T_start = 0 - T_end = time_interval - - assert spike_train[0] >= T_start, \ - "Spike train has events before the given start time" - assert spike_train[-1] <= T_end, \ - "Spike train has events after the given end time" - if spike_train[0] != T_start: - spike_train = np.insert(spike_train, 0, T_start) - if spike_train[-1] != T_end: - spike_train = np.append(spike_train, T_end) - return spike_train +from pyspike import SpikeTrain ############################################################ # spike_train_from_string ############################################################ -def spike_train_from_string(s, sep=' ', is_sorted=False): - """ Converts a string of times into an array of spike times. +def spike_train_from_string(s, interval, sep=' ', is_sorted=False): + """ Converts a string of times into a :class:`pyspike.SpikeTrain`. - :param s: the string with (ordered) spike times + :param s: the string with (ordered) spike times. + :param interval: interval defining the edges of the spike train. + Given as a pair of floats (T0, T1) or a single float T1, where T0=0 is + assumed. :param sep: The separator between the time numbers, default=' '. :param is_sorted: if True, the spike times are not sorted after loading, if False, spike times are sorted with `np.sort` - :returns: array of spike times + :returns: :class:`pyspike.SpikeTrain` """ if not(is_sorted): - return np.sort(np.fromstring(s, sep=sep)) + return SpikeTrain(np.sort(np.fromstring(s, sep=sep)), interval) else: - return np.fromstring(s, sep=sep) + return SpikeTrain(np.fromstring(s, sep=sep), interval) ############################################################ # load_spike_trains_txt ############################################################ -def load_spike_trains_from_txt(file_name, time_interval=None, +def load_spike_trains_from_txt(file_name, interval=None, separator=' ', comment='#', is_sorted=False): """ Loads a number of spike trains from a text file. Each line of the text file should contain one spike train as a sequence of spike times separated by `separator`. Empty lines as well as lines starting with `comment` are - neglected. The `time_interval` represents the start and the end of the - spike trains and it is used to add auxiliary spikes at the beginning and - end of each spike train. However, if `time_interval == None`, no auxiliary - spikes are added, but note that the Spike and ISI distance both require - auxiliary spikes. + neglected. The `interval` represents the start and the end of the + spike trains. :param file_name: The name of the text file. - :param time_interval: A pair (T_start, T_end) of values representing the - start and end time of the spike train measurement - or a single value representing the end time, the - T_start is then assuemd as 0. Auxiliary spikes will - be added to the spike train at the beginning and end - of this interval. + :param interval: A pair (T_start, T_end) of values representing the + start and end time of the spike train measurement + or a single value representing the end time, the + T_start is then assuemd as 0. :param separator: The character used to seprate the values in the text file :param comment: Lines starting with this character are ignored. :param sort: If true, the spike times are order via `np.sort`, default=True @@ -94,9 +58,8 @@ def load_spike_trains_from_txt(file_name, time_interval=None, for line in spike_file: if len(line) > 1 and not line.startswith(comment): # use only the lines with actual data and not commented - spike_train = spike_train_from_string(line, separator, is_sorted) - if time_interval is not None: # add auxil. spikes if times given - spike_train = add_auxiliary_spikes(spike_train, time_interval) + spike_train = spike_train_from_string(line, interval, + separator, is_sorted) spike_trains.append(spike_train) return spike_trains @@ -111,14 +74,14 @@ def merge_spike_trains(spike_trains): :returns: spike train with the merged spike times """ # get the lengths of the spike trains - lens = np.array([len(st) for st in spike_trains]) + lens = np.array([len(st.spikes) for st in spike_trains]) merged_spikes = np.empty(np.sum(lens)) index = 0 # the index for merged_spikes indices = np.zeros_like(lens) # indices of the spike trains index_list = np.arange(len(indices)) # indices of indices of spike trains # that have not yet reached the end # list of the possible events in the spike trains - vals = [spike_trains[i][indices[i]] for i in index_list] + vals = [spike_trains[i].spikes[indices[i]] for i in index_list] while len(index_list) > 0: i = np.argmin(vals) # the next spike is the minimum merged_spikes[index] = vals[i] # put it to the merged spike train @@ -127,33 +90,34 @@ def merge_spike_trains(spike_trains): indices[i] += 1 # next index for the chosen spike train if indices[i] >= lens[i]: # remove spike train index if ended index_list = index_list[index_list != i] - vals = [spike_trains[n][indices[n]] for n in index_list] - return merged_spikes + vals = [spike_trains[n].spikes[indices[n]] for n in index_list] + return SpikeTrain(merged_spikes, [spike_trains[0].t_start, + spike_trains[0].t_end]) ############################################################ # generate_poisson_spikes ############################################################ -def generate_poisson_spikes(rate, time_interval, add_aux_spikes=True): +def generate_poisson_spikes(rate, interval): """ Generates a Poisson spike train with the given rate in the given time interval :param rate: The rate of the spike trains - :param time_interval: A pair (T_start, T_end) of values representing the - start and end time of the spike train measurement or - a single value representing the end time, the T_start - is then assuemd as 0. Auxiliary spikes will be added - to the spike train at the beginning and end of this - interval, if they are not yet present. - :type time_interval: pair of doubles or double - :returns: Poisson spike train + :param interval: A pair (T_start, T_end) of values representing the + start and end time of the spike train measurement or + a single value representing the end time, the T_start + is then assuemd as 0. Auxiliary spikes will be added + to the spike train at the beginning and end of this + interval, if they are not yet present. + :type interval: pair of doubles or double + :returns: Poisson spike train as a :class:`pyspike.SpikeTrain` """ try: - T_start = time_interval[0] - T_end = time_interval[1] + T_start = interval[0] + T_end = interval[1] except: T_start = 0 - T_end = time_interval + T_end = interval # roughly how many spikes are required to fill the interval N = max(1, int(1.2 * rate * (T_end-T_start))) N_append = max(1, int(0.1 * rate * (T_end-T_start))) @@ -165,7 +129,4 @@ def generate_poisson_spikes(rate, time_interval, add_aux_spikes=True): np.random.exponential(1.0/rate, N_append)) spikes = T_start + np.cumsum(intervals) spikes = spikes[spikes < T_end] - if add_aux_spikes: - return add_auxiliary_spikes(spikes, time_interval) - else: - return spikes + return SpikeTrain(spikes, interval) diff --git a/test/test_distance.py b/test/test_distance.py index dbb72f1..0fff840 100644 --- a/test/test_distance.py +++ b/test/test_distance.py @@ -250,7 +250,8 @@ def test_multi_spike_sync(): # multivariate regression test spike_trains = spk.load_spike_trains_from_txt("test/SPIKE_Sync_Test.txt", - time_interval=(0, 4000)) + interval=(0, 4000)) + print(spike_trains[0].spikes) f = spk.spike_sync_profile_multi(spike_trains) assert_equal(np.sum(f.y[1:-1]), 39932) assert_equal(np.sum(f.mp[1:-1]), 85554) @@ -339,4 +340,4 @@ if __name__ == "__main__": test_spike_sync() test_multi_isi() test_multi_spike() - test_multi_spike_sync() + # test_multi_spike_sync() -- cgit v1.2.3 From 36d80c9ec1d28488f9b5c97cd202c196efff694e Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Fri, 24 Apr 2015 15:58:35 +0200 Subject: distance tests now pass with new spike trains --- pyspike/cython/cython_distance.pyx | 8 ++++---- pyspike/cython/python_backend.py | 41 +++++++++++++++++++------------------- test/test_distance.py | 40 ++++++++++++++++++++++++++++++------- 3 files changed, 57 insertions(+), 32 deletions(-) (limited to 'pyspike/cython/cython_distance.pyx') diff --git a/pyspike/cython/cython_distance.pyx b/pyspike/cython/cython_distance.pyx index 2841da8..dc2557f 100644 --- a/pyspike/cython/cython_distance.pyx +++ b/pyspike/cython/cython_distance.pyx @@ -345,9 +345,9 @@ cdef inline double get_tau(double[:] spikes1, double[:] spikes2, m = fmin(m, spikes1[i+1]-spikes1[i]) if j < N2 and j > -1: m = fmin(m, spikes2[j+1]-spikes2[j]) - if i > 1: + if i > 0: m = fmin(m, spikes1[i]-spikes1[i-1]) - if j > 1: + if j > 0: m = fmin(m, spikes2[j]-spikes2[j-1]) m *= 0.5 if max_tau > 0.0: @@ -371,7 +371,7 @@ def coincidence_cython(double[:] spikes1, double[:] spikes2, cdef double[:] mp = np.ones(N1 + N2 + 2) # multiplicity cdef double tau while i + j < N1 + N2 - 2: - if (i < N1-1) and (spikes1[i+1] < spikes2[j+1] or j == N2-1): + 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) @@ -381,7 +381,7 @@ def coincidence_cython(double[:] spikes1, double[:] spikes2, # both get marked with 1 c[n] = 1 c[n-1] = 1 - elif (j < N2-1) and (spikes1[i+1] > spikes2[j+1] or i == N1-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) diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py index bcf9c30..c65bfb0 100644 --- a/pyspike/cython/python_backend.py +++ b/pyspike/cython/python_backend.py @@ -330,47 +330,48 @@ def cumulative_sync_python(spikes1, spikes2): ############################################################ # coincidence_python ############################################################ -def coincidence_python(spikes1, spikes2, max_tau): +def coincidence_python(spikes1, spikes2, t_start, t_end, max_tau): def get_tau(spikes1, spikes2, i, j, max_tau): m = 1E100 # some huge number - if i < len(spikes1)-2: + if i < len(spikes1)-1 and i > -1: m = min(m, spikes1[i+1]-spikes1[i]) - if j < len(spikes2)-2: + if j < len(spikes2)-1 and j > -1: m = min(m, spikes2[j+1]-spikes2[j]) - if i > 1: + if i > 0: m = min(m, spikes1[i]-spikes1[i-1]) - if j > 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 = 0 - j = 0 + i = -1 + j = -1 n = 0 - st = np.zeros(N1 + N2 - 2) # spike times - c = np.zeros(N1 + N2 - 2) # coincidences - mp = np.ones(N1 + N2 - 2) # multiplicity - while n < N1 + N2 - 2: - if spikes1[i+1] < spikes2[j+1]: + st = np.zeros(N1 + N2 + 2) # spike times + c = np.zeros(N1 + N2 + 2) # coincidences + mp = np.ones(N1 + N2 + 2) # multiplicity + 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 > 0 and spikes1[i]-spikes2[j] < tau: + 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 spikes1[i+1] > spikes2[j+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 > 0 and spikes2[j]-spikes1[i] < tau: + 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 @@ -379,8 +380,6 @@ def coincidence_python(spikes1, spikes2, max_tau): # advance in both spike trains j += 1 i += 1 - if i == N1-1 or j == N2-1: - break n += 1 # add only one event, but with coincidence 2 and multiplicity 2 st[n] = spikes1[i] @@ -391,12 +390,12 @@ def coincidence_python(spikes1, spikes2, max_tau): c = c[:n+2] mp = mp[:n+2] - st[0] = spikes1[0] - st[-1] = spikes1[-1] + st[0] = t_start + st[len(st)-1] = t_end c[0] = c[1] - c[-1] = c[-2] + c[len(c)-1] = c[len(c)-2] mp[0] = mp[1] - mp[-1] = mp[-2] + mp[len(mp)-1] = mp[len(mp)-2] return st, c, mp diff --git a/test/test_distance.py b/test/test_distance.py index 0fff840..88cf40e 100644 --- a/test/test_distance.py +++ b/test/test_distance.py @@ -177,6 +177,18 @@ def test_spike_sync(): assert_almost_equal(spk.spike_sync(spikes1, spikes2), 0.5, decimal=16) + spikes2 = SpikeTrain([3.0], 4.0) + assert_almost_equal(spk.spike_sync(spikes1, spikes2), + 0.5, decimal=16) + + spikes2 = SpikeTrain([1.0], 4.0) + assert_almost_equal(spk.spike_sync(spikes1, spikes2), + 0.5, decimal=16) + + spikes2 = SpikeTrain([1.5, 3.0], 4.0) + assert_almost_equal(spk.spike_sync(spikes1, spikes2), + 0.4, decimal=16) + def check_multi_profile(profile_func, profile_func_multi): # generate spike trains: @@ -250,19 +262,28 @@ def test_multi_spike_sync(): # multivariate regression test spike_trains = spk.load_spike_trains_from_txt("test/SPIKE_Sync_Test.txt", - interval=(0, 4000)) - print(spike_trains[0].spikes) + interval=[0, 4000]) + # extract all spike times + spike_times = np.array([]) + for st in spike_trains: + spike_times = np.append(spike_times, st.spikes) + spike_times = np.unique(np.sort(spike_times)) + f = spk.spike_sync_profile_multi(spike_trains) + + assert_equal(spike_times, f.x[1:-1]) + assert_equal(len(f.x), len(f.y)) + assert_equal(np.sum(f.y[1:-1]), 39932) assert_equal(np.sum(f.mp[1:-1]), 85554) def check_dist_matrix(dist_func, dist_matrix_func): # generate spike trains: - t1 = spk.add_auxiliary_spikes(np.array([0.2, 0.4, 0.6, 0.7]), 1.0) - t2 = spk.add_auxiliary_spikes(np.array([0.3, 0.45, 0.8, 0.9, 0.95]), 1.0) - t3 = spk.add_auxiliary_spikes(np.array([0.2, 0.4, 0.6]), 1.0) - t4 = spk.add_auxiliary_spikes(np.array([0.1, 0.4, 0.5, 0.6]), 1.0) + t1 = SpikeTrain([0.2, 0.4, 0.6, 0.7], 1.0) + t2 = SpikeTrain([0.3, 0.45, 0.8, 0.9, 0.95], 1.0) + t3 = SpikeTrain([0.2, 0.4, 0.6], 1.0) + t4 = SpikeTrain([0.1, 0.4, 0.5, 0.6], 1.0) spike_trains = [t1, t2, t3, t4] f12 = dist_func(t1, t2) @@ -340,4 +361,9 @@ if __name__ == "__main__": test_spike_sync() test_multi_isi() test_multi_spike() - # test_multi_spike_sync() + test_multi_spike_sync() + test_isi_matrix() + test_spike_matrix() + test_spike_sync_matrix() + test_regression_spiky() + test_multi_variate_subsets() -- cgit v1.2.3 From 795e16ffe7afb469ef07a548c1f6a31d924196b3 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Fri, 24 Apr 2015 23:29:05 +0200 Subject: bugfixes for spike distance --- pyspike/cython/cython_distance.pyx | 24 +++++++++--------- pyspike/cython/python_backend.py | 15 ++++++------ test/test_distance.py | 50 ++++++++++++++++++++++++++++---------- 3 files changed, 58 insertions(+), 31 deletions(-) (limited to 'pyspike/cython/cython_distance.pyx') diff --git a/pyspike/cython/cython_distance.pyx b/pyspike/cython/cython_distance.pyx index dc2557f..a41d8e8 100644 --- a/pyspike/cython/cython_distance.pyx +++ b/pyspike/cython/cython_distance.pyx @@ -194,31 +194,31 @@ def spike_distance_cython(double[:] t1, double[:] t2, t_p2 = t_start if t1[0] > t_start: # dt_p1 = t2[0]-t_start - dt_p1 = 0.0 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]) - s1 = dt_f1*(t_f1-t_start)/isi1 + dt_p1 = dt_f1 + s1 = dt_p1*(t_f1-t_start)/isi1 index1 = -1 else: - dt_p1 = 0.0 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 - dt_p2 = 0.0 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_f2*(t_f2-t_start)/isi2 + s2 = dt_p2*(t_f2-t_start)/isi2 index2 = -1 else: - dt_p2 = 0.0 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 @@ -231,16 +231,16 @@ def spike_distance_cython(double[:] t1, double[:] t2, 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 + dt_p1 = dt_f1 t_p1 = t_f1 # t_p1 contains the current time point - # get the next time + # get the next time if index1 < N1-1: t_f1 = t1[index1+1] else: t_f1 = t_end spike_events[index] = t_p1 - s1 = dt_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) @@ -249,6 +249,7 @@ def spike_distance_cython(double[:] t1, double[:] t2, 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]) @@ -260,9 +261,10 @@ def spike_distance_cython(double[:] t1, double[:] t2, 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_p1 contains the current time point + t_p2 = t_f2 # t_p2 contains the current time point # get the next time if index2 < N2-1: t_f2 = t2[index2+1] @@ -270,7 +272,6 @@ def spike_distance_cython(double[:] t1, double[:] t2, t_f2 = t_end spike_events[index] = t_p2 s1 = (dt_p1*(t_f1-t_p2) + dt_f1*(t_p2-t_p1)) / isi1 - s2 = dt_p2 y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2) # now the next interval start value @@ -278,6 +279,7 @@ def spike_distance_cython(double[:] t1, double[:] t2, 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]) diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py index c65bfb0..317b568 100644 --- a/pyspike/cython/python_backend.py +++ b/pyspike/cython/python_backend.py @@ -142,12 +142,11 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end): t_p1 = t_start t_p2 = t_start if t1[0] > t_start: - # dt_p1 = t2[0]-t_start - dt_p1 = 0.0 t_f1 = t1[0] dt_f1 = get_min_dist(t_f1, t2, 0, t_start, t_end) + dt_p1 = dt_f1 isi1 = max(t_f1-t_start, t1[1]-t1[0]) - s1 = dt_f1*(t_f1-t_start)/isi1 + s1 = dt_p1*(t_f1-t_start)/isi1 index1 = -1 else: dt_p1 = 0.0 @@ -158,11 +157,11 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end): index1 = 0 if t2[0] > t_start: # dt_p1 = t2[0]-t_start - dt_p2 = 0.0 t_f2 = t2[0] dt_f2 = get_min_dist(t_f2, t1, 0, t_start, t_end) + dt_p2 = dt_f2 isi2 = max(t_f2-t_start, t2[1]-t2[0]) - s2 = dt_f2*(t_f2-t_start)/isi2 + s2 = dt_p2*(t_f2-t_start)/isi2 index2 = -1 else: dt_p2 = 0.0 @@ -180,6 +179,7 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end): 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 @@ -189,13 +189,13 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end): else: t_f1 = t_end spike_events[index] = t_p1 - s1 = dt_p1 s2 = (dt_p2*(t_f2-t_p1) + dt_f2*(t_p1-t_p2)) / isi2 y_ends[index-1] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2) # now the next interval start value if index1 < N1-1: dt_f1 = get_min_dist(t_f1, t2, index2, t_start, t_end) isi1 = t_f1-t_p1 + s1 = dt_p1 else: dt_f1 = dt_p1 isi1 = max(t_end-t1[N1-1], t1[N1-1]-t1[N1-2]) @@ -206,6 +206,7 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end): 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_p1 contains the current time point @@ -216,12 +217,12 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end): t_f2 = t_end spike_events[index] = t_p2 s1 = (dt_p1*(t_f1-t_p2) + dt_f1*(t_p2-t_p1)) / isi1 - s2 = dt_p2 y_ends[index-1] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2) # now the next interval start value if index2 < N2-1: dt_f2 = get_min_dist(t_f2, t1, index1, t_start, t_end) isi2 = t_f2-t_p2 + s2 = dt_p2 else: dt_f2 = dt_p2 isi2 = max(t_end-t2[N2-1], t2[N2-1]-t2[N2-2]) diff --git a/test/test_distance.py b/test/test_distance.py index 0059001..20b52e8 100644 --- a/test/test_distance.py +++ b/test/test_distance.py @@ -73,7 +73,7 @@ def test_spike(): assert_equal(f.x, expected_times) - assert_almost_equal(f.avrg(), 0.1662415, decimal=6) + assert_almost_equal(f.avrg(), 1.6624149659863946e-01, decimal=15) assert_almost_equal(f.y2[-1], 0.1394558, decimal=6) t1 = SpikeTrain([0.2, 0.4, 0.6, 0.7], 1.0) @@ -84,7 +84,7 @@ def test_spike(): s1 = np.array([0.1, 0.1, (0.1*0.1+0.05*0.1)/0.2, 0.05, (0.05*0.15 * 2)/0.2, 0.15, 0.1, (0.1*0.1+0.1*0.2)/0.3, (0.1*0.2+0.1*0.1)/0.3, (0.1*0.05+0.1*0.25)/0.3, 0.1]) - s2 = np.array([0.1, 0.1*0.2/0.3, 0.1, (0.1*0.05 * 2)/.15, 0.05, + s2 = np.array([0.1, (0.1*0.2+0.1*0.1)/0.3, 0.1, (0.1*0.05 * 2)/.15, 0.05, (0.05*0.2+0.1*0.15)/0.35, (0.05*0.1+0.1*0.25)/0.35, 0.1, 0.1, 0.05, 0.05]) isi1 = np.array([0.2, 0.2, 0.2, 0.2, 0.2, 0.1, 0.3, 0.3, 0.3, 0.3]) @@ -113,12 +113,18 @@ def test_spike(): t2 = SpikeTrain([0.1, 0.4, 0.5, 0.6], [0.0, 1.0]) expected_times = [0.0, 0.1, 0.2, 0.4, 0.5, 0.6, 1.0] - s1 = np.array([0.1, 0.1*0.1/0.2, 0.1, 0.0, 0.0, 0.0, 0.0]) - s2 = np.array([0.1*0.1/0.3, 0.1, 0.1*0.2/0.3, 0.0, 0.1, 0.0, 0.0]) + # due to the edge correction in the beginning, s1 and s2 are different + # for left and right values + s1_r = np.array([0.1, (0.1*0.1+0.1*0.1)/0.2, 0.1, 0.0, 0.0, 0.0, 0.0]) + s1_l = np.array([0.1, (0.1*0.1+0.1*0.1)/0.2, 0.1, 0.0, 0.0, 0.0, 0.0]) + s2_r = np.array([0.1*0.1/0.3, 0.1*0.3/0.3, 0.1*0.2/0.3, + 0.0, 0.1, 0.0, 0.0]) + s2_l = np.array([0.1*0.1/0.3, 0.1*0.1/0.3, 0.1*0.2/0.3, 0.0, + 0.1, 0.0, 0.0]) isi1 = np.array([0.2, 0.2, 0.2, 0.2, 0.2, 0.4]) isi2 = np.array([0.3, 0.3, 0.3, 0.1, 0.1, 0.4]) - 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_y1 = (s1_r[:-1]*isi2+s2_r[:-1]*isi1) / (0.5*(isi1+isi2)**2) + expected_y2 = (s1_l[1:]*isi2+s2_l[1:]*isi1) / (0.5*(isi1+isi2)**2) expected_times = np.array(expected_times) expected_y1 = np.array(expected_y1) @@ -321,19 +327,37 @@ def test_spike_sync_matrix(): def test_regression_spiky(): + # standard example + st1 = SpikeTrain(np.arange(100, 1201, 100), 1300) + st2 = SpikeTrain(np.arange(100, 1201, 110), 1300) + + isi_dist = spk.isi_distance(st1, st2) + assert_almost_equal(isi_dist, 7.6923076923076941e-02, decimal=15) + + spike_dist = spk.spike_distance(st1, st2) + assert_equal(spike_dist, 2.1105878248735391e-01) + + spike_sync = spk.spike_sync(st1, st2) + assert_equal(spike_sync, 8.6956521739130432e-01) + + # multivariate check + spike_trains = spk.load_spike_trains_from_txt("test/PySpike_testdata.txt", (0.0, 4000.0)) - isi_profile = spk.isi_profile_multi(spike_trains) - isi_dist = isi_profile.avrg() - print(isi_dist) + isi_dist = spk.isi_distance_multi(spike_trains) # get the full precision from SPIKY - # assert_equal(isi_dist, 0.1832) + assert_almost_equal(isi_dist, 1.8318789829845508e-01, decimal=15) spike_profile = spk.spike_profile_multi(spike_trains) - spike_dist = spike_profile.avrg() - print(spike_dist) + assert_equal(len(spike_profile.y1)+len(spike_profile.y2), 1252) + + spike_dist = spk.spike_distance_multi(spike_trains) + # get the full precision from SPIKY + assert_almost_equal(spike_dist, 2.4432433330596512e-01, decimal=15) + + spike_sync = spk.spike_sync_multi(spike_trains) # get the full precision from SPIKY - # assert_equal(spike_dist, 0.2445) + assert_equal(spike_sync, 0.7183531505298066) def test_multi_variate_subsets(): -- cgit v1.2.3 From 2e7351393927ba9e9e0c3b7b59d05e8aeeb41d1f Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Tue, 28 Apr 2015 16:11:11 +0200 Subject: edge correction for the ISI-distance --- pyspike/cython/cython_distance.pyx | 18 ++++++++++++------ pyspike/cython/python_backend.py | 18 ++++++++++++------ test/test_distance.py | 8 +++++--- 3 files changed, 29 insertions(+), 15 deletions(-) (limited to 'pyspike/cython/cython_distance.pyx') diff --git a/pyspike/cython/cython_distance.pyx b/pyspike/cython/cython_distance.pyx index a41d8e8..6ee0181 100644 --- a/pyspike/cython/cython_distance.pyx +++ b/pyspike/cython/cython_distance.pyx @@ -62,14 +62,16 @@ def isi_distance_cython(double[:] s1, double[:] s2, # first interspike interval - check if a spike exists at the start time if s1[0] > t_start: - nu1 = s1[0] - t_start + # edge correction + nu1 = fmax(s1[0]-t_start, s1[1]-s1[0]) index1 = -1 else: nu1 = s1[1]-s1[0] index1 = 0 if s2[0] > t_start: - nu2 = s2[0] - t_start + # edge correction + nu2 = fmax(s2[0]-t_start, s2[1]-s2[0]) index2 = -1 else: nu2 = s2[1]-s2[0] @@ -89,7 +91,8 @@ def isi_distance_cython(double[:] s1, double[:] s2, if index1 < N1-1: nu1 = s1[index1+1]-s1[index1] else: - nu1 = t_end-s1[index1] + # edge correction + nu1 = fmax(t_end-s1[index1], nu1) elif (index2 < N2-1) and ((index1 == N1-1) or (s1[index1+1] > s2[index2+1])): index2 += 1 @@ -97,7 +100,8 @@ def isi_distance_cython(double[:] s1, double[:] s2, if index2 < N2-1: nu2 = s2[index2+1]-s2[index2] else: - nu2 = t_end-s2[index2] + # edge correction + nu2 = fmax(t_end-s2[index2], nu2) else: # s1[index1+1] == s2[index2+1] index1 += 1 index2 += 1 @@ -105,11 +109,13 @@ def isi_distance_cython(double[:] s1, double[:] s2, if index1 < N1-1: nu1 = s1[index1+1]-s1[index1] else: - nu1 = t_end-s1[index1] + # edge correction + nu1 = fmax(t_end-s1[index1], nu1) if index2 < N2-1: nu2 = s2[index2+1]-s2[index2] else: - nu2 = t_end-s2[index2] + # edge correction + nu2 = fmax(t_end-s2[index2], nu2) # compute the corresponding isi-distance isi_values[index] = fabs(nu1 - nu2) / fmax(nu1, nu2) index += 1 diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py index 317b568..1fd8c42 100644 --- a/pyspike/cython/python_backend.py +++ b/pyspike/cython/python_backend.py @@ -27,13 +27,15 @@ def isi_distance_python(s1, s2, t_start, t_end): # the values have one entry less - the number of intervals between events isi_values = np.empty(len(spike_events) - 1) if s1[0] > t_start: - nu1 = s1[0] - t_start + # edge correction + nu1 = max(s1[0] - t_start, s1[1] - s1[0]) index1 = -1 else: nu1 = s1[1] - s1[0] index1 = 0 if s2[0] > t_start: - nu2 = s2[0] - t_start + # edge correction + nu2 = max(s2[0] - t_start, s2[1] - s2[0]) index2 = -1 else: nu2 = s2[1] - s2[0] @@ -49,7 +51,8 @@ def isi_distance_python(s1, s2, t_start, t_end): if index1 < N1-1: nu1 = s1[index1+1]-s1[index1] else: - nu1 = t_end-s1[index1] + # edge correction + nu1 = max(t_end-s1[N1-1], s1[N1-1]-s1[N1-2]) elif (index2 < N2-1) and (index1 == N1-1 or s1[index1+1] > s2[index2+1]): @@ -58,7 +61,8 @@ def isi_distance_python(s1, s2, t_start, t_end): if index2 < N2-1: nu2 = s2[index2+1]-s2[index2] else: - nu2 = t_end-s2[index2] + # edge correction + nu2 = max(t_end-s2[N2-1], s2[N2-1]-s2[N2-2]) else: # s1[index1 + 1] == s2[index2 + 1] index1 += 1 @@ -67,11 +71,13 @@ def isi_distance_python(s1, s2, t_start, t_end): if index1 < N1-1: nu1 = s1[index1+1]-s1[index1] else: - nu1 = t_end-s1[index1] + # edge correction + nu1 = max(t_end-s1[N1-1], s1[N1-1]-s1[N1-2]) if index2 < N2-1: nu2 = s2[index2+1]-s2[index2] else: - nu2 = t_end-s2[index2] + # edge correction + nu2 = max(t_end-s2[N2-1], s2[N2-1]-s2[N2-2]) # compute the corresponding isi-distance isi_values[index] = abs(nu1 - nu2) / \ max(nu1, nu2) diff --git a/test/test_distance.py b/test/test_distance.py index 20b52e8..19da35f 100644 --- a/test/test_distance.py +++ b/test/test_distance.py @@ -47,7 +47,7 @@ def test_isi(): t2 = SpikeTrain([0.1, 0.4, 0.5, 0.6], [0.0, 1.0]) expected_times = [0.0, 0.1, 0.2, 0.4, 0.5, 0.6, 1.0] - expected_isi = [0.1/0.2, 0.1/0.3, 0.1/0.3, 0.1/0.2, 0.1/0.2, 0.0/0.5] + expected_isi = [0.1/0.3, 0.1/0.3, 0.1/0.3, 0.1/0.2, 0.1/0.2, 0.0/0.5] expected_times = np.array(expected_times) expected_isi = np.array(expected_isi) @@ -332,7 +332,9 @@ def test_regression_spiky(): st2 = SpikeTrain(np.arange(100, 1201, 110), 1300) isi_dist = spk.isi_distance(st1, st2) - assert_almost_equal(isi_dist, 7.6923076923076941e-02, decimal=15) + assert_almost_equal(isi_dist, 9.0909090909090939e-02, decimal=15) + isi_profile = spk.isi_profile(st1, st2) + assert_equal(isi_profile.y, 0.1/1.1 * np.ones_like(isi_profile.y)) spike_dist = spk.spike_distance(st1, st2) assert_equal(spike_dist, 2.1105878248735391e-01) @@ -346,7 +348,7 @@ def test_regression_spiky(): (0.0, 4000.0)) isi_dist = spk.isi_distance_multi(spike_trains) # get the full precision from SPIKY - assert_almost_equal(isi_dist, 1.8318789829845508e-01, decimal=15) + assert_almost_equal(isi_dist, 0.17051816816999129656, decimal=15) spike_profile = spk.spike_profile_multi(spike_trains) assert_equal(len(spike_profile.y1)+len(spike_profile.y2), 1252) -- cgit v1.2.3