From 6eb6bc486027d3d5304a94cfb417a2257f2b6fd9 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Tue, 3 Feb 2015 12:19:53 +0100 Subject: moved cython functions to subdirectory --- pyspike/cython/cython_add.pyx | 235 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 235 insertions(+) create mode 100644 pyspike/cython/cython_add.pyx (limited to 'pyspike/cython/cython_add.pyx') diff --git a/pyspike/cython/cython_add.pyx b/pyspike/cython/cython_add.pyx new file mode 100644 index 0000000..ac64005 --- /dev/null +++ b/pyspike/cython/cython_add.pyx @@ -0,0 +1,235 @@ +#cython: boundscheck=False +#cython: wraparound=False +#cython: cdivision=True + +""" +cython_add.pyx + +cython implementation of the add function for piece-wise const and +piece-wise linear functions + +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_add.pyx + +which gives:: + + cython_add.html + +""" + +import numpy as np +cimport numpy as np + +from libc.math cimport fabs + +DTYPE = np.float +ctypedef np.float_t DTYPE_t + + +############################################################ +# add_piece_wise_const_cython +############################################################ +def add_piece_wise_const_cython(double[:] x1, double[:] y1, + double[:] x2, double[:] y2): + + cdef int N1 = len(x1) + cdef int N2 = len(x2) + cdef double[:] x_new = np.empty(N1+N2) + cdef double[:] y_new = np.empty(N1+N2-1) + cdef int index1 = 0 + cdef int index2 = 0 + cdef int index = 0 + cdef int i + with nogil: # release the interpreter lock to allow multi-threading + x_new[0] = x1[0] + y_new[0] = y1[0] + y2[0] + while (index1+1 < N1-1) and (index2+1 < N2-1): + index += 1 + # print(index1+1, x1[index1+1], y1[index1+1], x_new[index]) + if x1[index1+1] < x2[index2+1]: + index1 += 1 + x_new[index] = x1[index1] + elif x1[index1+1] > x2[index2+1]: + index2 += 1 + x_new[index] = x2[index2] + else: # x1[index1+1] == x2[index2+1]: + index1 += 1 + index2 += 1 + x_new[index] = x1[index1] + y_new[index] = y1[index1] + y2[index2] + # one array reached the end -> copy the contents of the other to the end + if index1+1 < N1-1: + x_new[index+1:index+1+N1-index1-1] = x1[index1+1:] + for i in xrange(N1-index1-2): + y_new[index+1+i] = y1[index1+1+i] + y2[N2-2] + index += N1-index1-2 + elif index2+1 < N2-1: + x_new[index+1:index+1+N2-index2-1] = x2[index2+1:] + for i in xrange(N2-index2-2): + y_new[index+1+i] = y2[index2+1+i] + y1[N1-2] + index += N2-index2-2 + else: # both arrays reached the end simultaneously + # only the last x-value missing + x_new[index+1] = x1[N1-1] + # the last value is again the end of the interval + # x_new[index+1] = x1[-1] + # only use the data that was actually filled + x1 = x_new[:index+2] + y1 = y_new[:index+1] + # end nogil + return np.array(x_new[:index+2]), np.array(y_new[:index+1]) + + +############################################################ +# add_piece_wise_lin_cython +############################################################ +def add_piece_wise_lin_cython(double[:] x1, double[:] y11, double[:] y12, + double[:] x2, double[:] y21, double[:] y22): + cdef int N1 = len(x1) + cdef int N2 = len(x2) + cdef double[:] x_new = np.empty(N1+N2) + cdef double[:] y1_new = np.empty(N1+N2-1) + cdef double[:] y2_new = np.empty_like(y1_new) + cdef int index1 = 0 # index for self + cdef int index2 = 0 # index for f + cdef int index = 0 # index for new + cdef int i + cdef double y + with nogil: # release the interpreter lock to allow multi-threading + x_new[0] = x1[0] + y1_new[0] = y11[0] + y21[0] + while (index1+1 < N1-1) and (index2+1 < N2-1): + # print(index1+1, x1[index1+1], self.y[index1+1], x_new[index]) + if x1[index1+1] < x2[index2+1]: + # first compute the end value of the previous interval + # linear interpolation of the interval + y = y21[index2] + (y22[index2]-y21[index2]) * \ + (x1[index1+1]-x2[index2]) / (x2[index2+1]-x2[index2]) + y2_new[index] = y12[index1] + y + index1 += 1 + index += 1 + x_new[index] = x1[index1] + # and the starting value for the next interval + y1_new[index] = y11[index1] + y + elif x1[index1+1] > x2[index2+1]: + # first compute the end value of the previous interval + # linear interpolation of the interval + y = y11[index1] + (y12[index1]-y11[index1]) * \ + (x2[index2+1]-x1[index1]) / \ + (x1[index1+1]-x1[index1]) + y2_new[index] = y22[index2] + y + index2 += 1 + index += 1 + x_new[index] = x2[index2] + # and the starting value for the next interval + y1_new[index] = y21[index2] + y + else: # x1[index1+1] == x2[index2+1]: + y2_new[index] = y12[index1] + y22[index2] + index1 += 1 + index2 += 1 + index += 1 + x_new[index] = x1[index1] + y1_new[index] = y11[index1] + y21[index2] + # one array reached the end -> copy the contents of the other to the end + if index1+1 < N1-1: + x_new[index+1:index+1+N1-index1-1] = x1[index1+1:] + for i in xrange(N1-index1-2): + # compute the linear interpolations value + y = y21[index2] + (y22[index2]-y21[index2]) * \ + (x1[index1+1+i]-x2[index2]) / (x2[index2+1]-x2[index2]) + y1_new[index+1+i] = y11[index1+1+i] + y + y2_new[index+i] = y12[index1+i] + y + index += N1-index1-2 + elif index2+1 < N2-1: + x_new[index+1:index+1+N2-index2-1] = x2[index2+1:] + # compute the linear interpolations values + for i in xrange(N2-index2-2): + y = y11[index1] + (y12[index1]-y11[index1]) * \ + (x2[index2+1+i]-x1[index1]) / \ + (x1[index1+1]-x1[index1]) + y1_new[index+1+i] = y21[index2+1+i] + y + y2_new[index+i] = y22[index2+i] + y + index += N2-index2-2 + else: # both arrays reached the end simultaneously + # only the last x-value missing + x_new[index+1] = x1[N1-1] + # finally, the end value for the last interval + y2_new[index] = y12[N1-2]+y22[N2-2] + # only use the data that was actually filled + # end nogil + return (np.array(x_new[:index+2]), + np.array(y1_new[:index+1]), + np.array(y2_new[:index+1])) + + +############################################################ +# add_discrete_function_cython +############################################################ +def add_discrete_function_cython(double[:] x1, double[:] y1, double[:] mp1, + double[:] x2, double[:] y2, double[:] mp2): + + cdef double[:] x_new = np.empty(len(x1) + len(x2)) + cdef double[:] y_new = np.empty_like(x_new) + cdef double[:] mp_new = np.empty_like(x_new) + cdef int index1 = 0 + cdef int index2 = 0 + cdef int index = 0 + cdef int N1 = len(y1) + cdef int N2 = len(y2) + x_new[0] = x1[0] + while (index1+1 < N1) and (index2+1 < N2): + if x1[index1+1] < x2[index2+1]: + index1 += 1 + index += 1 + x_new[index] = x1[index1] + y_new[index] = y1[index1] + mp_new[index] = mp1[index1] + elif x1[index1+1] > x2[index2+1]: + index2 += 1 + index += 1 + x_new[index] = x2[index2] + y_new[index] = y2[index2] + mp_new[index] = mp2[index2] + else: # x1[index1+1] == x2[index2+1] + index1 += 1 + index2 += 1 + index += 1 + x_new[index] = x1[index1] + y_new[index] = y1[index1] + y2[index2] + mp_new[index] = mp1[index1] + mp2[index2] + # one array reached the end -> copy the contents of the other to the end + if index1+1 < len(y1): + x_new[index+1:index+1+len(x1)-index1-1] = x1[index1+1:] + y_new[index+1:index+1+len(x1)-index1-1] = y1[index1+1:] + mp_new[index+1:index+1+len(x1)-index1-1] = mp1[index1+1:] + index += len(x1)-index1-1 + elif index2+1 < len(y2): + x_new[index+1:index+1+len(x2)-index2-1] = x2[index2+1:] + y_new[index+1:index+1+len(x2)-index2-1] = y2[index2+1:] + mp_new[index+1:index+1+len(x2)-index2-1] = mp2[index2+1:] + index += len(x2)-index2-1 + # else: # both arrays reached the end simultaneously + # x_new[index+1] = x1[-1] + # y_new[index+1] = y1[-1] + y2[-1] + # mp_new[index+1] = mp1[-1] + mp2[-1] + + y_new[0] = y_new[1] + mp_new[0] = mp_new[1] + + # the last value is again the end of the interval + # only use the data that was actually filled + return (np.array(x_new[:index+1]), + np.array(y_new[:index+1]), + np.array(mp_new[:index+1])) -- cgit v1.2.3 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. --- examples/performance.py | 42 +++++++++++++++++++++++ pyspike/cython/cython_add.pyx | 20 +++++------ pyspike/generic.py | 77 +++++++++++++++++++++++++++++++++++++++---- pyspike/isi_distance.py | 6 ++-- pyspike/spike_distance.py | 6 ++-- test/test_distance.py | 14 ++++++-- 6 files changed, 139 insertions(+), 26 deletions(-) create mode 100644 examples/performance.py (limited to 'pyspike/cython/cython_add.pyx') diff --git a/examples/performance.py b/examples/performance.py new file mode 100644 index 0000000..469b5ab --- /dev/null +++ b/examples/performance.py @@ -0,0 +1,42 @@ +""" Compute distances of large sets of spike trains for performance tests +""" + +from __future__ import print_function + +import pyspike as spk +from datetime import datetime +import cProfile + +M = 100 # number of spike trains +r = 1.0 # rate of Poisson spike times +T = 1E3 # length of spike trains + +print("%d spike trains with %d spikes" % (M, int(r*T))) + +spike_trains = [] + +t_start = datetime.now() +for i in xrange(M): + spike_trains.append(spk.generate_poisson_spikes(r, T)) +t_end = datetime.now() +runtime = (t_end-t_start).total_seconds() + +print("Spike generation runtime: %.3fs" % runtime) + +print("================ ISI COMPUTATIONS ================") +print(" MULTIVARIATE DISTANCE") +cProfile.run('spk.isi_distance_multi(spike_trains)') +print(" MULTIVARIATE PROFILE") +cProfile.run('spk.isi_profile_multi(spike_trains)') + +print("================ SPIKE COMPUTATIONS ================") +print(" MULTIVARIATE DISTANCE") +cProfile.run('spk.spike_distance_multi(spike_trains)') +print(" MULTIVARIATE PROFILE") +cProfile.run('spk.spike_profile_multi(spike_trains)') + +print("================ SPIKE-SYNC COMPUTATIONS ================") +print(" MULTIVARIATE DISTANCE") +cProfile.run('spk.spike_sync_multi(spike_trains)') +print(" MULTIVARIATE PROFILE") +cProfile.run('spk.spike_sync_profile_multi(spike_trains)') diff --git a/pyspike/cython/cython_add.pyx b/pyspike/cython/cython_add.pyx index ac64005..8da1e53 100644 --- a/pyspike/cython/cython_add.pyx +++ b/pyspike/cython/cython_add.pyx @@ -83,13 +83,9 @@ def add_piece_wise_const_cython(double[:] x1, double[:] y1, else: # both arrays reached the end simultaneously # only the last x-value missing x_new[index+1] = x1[N1-1] - # the last value is again the end of the interval - # x_new[index+1] = x1[-1] - # only use the data that was actually filled - x1 = x_new[:index+2] - y1 = y_new[:index+1] # end nogil - return np.array(x_new[:index+2]), np.array(y_new[:index+1]) + # return np.asarray(x_new[:index+2]), np.asarray(y_new[:index+1]) + return np.asarray(x_new[:index+2]), np.asarray(y_new[:index+1]) ############################################################ @@ -169,9 +165,9 @@ def add_piece_wise_lin_cython(double[:] x1, double[:] y11, double[:] y12, y2_new[index] = y12[N1-2]+y22[N2-2] # only use the data that was actually filled # end nogil - return (np.array(x_new[:index+2]), - np.array(y1_new[:index+1]), - np.array(y2_new[:index+1])) + return (np.asarray(x_new[:index+2]), + np.asarray(y1_new[:index+1]), + np.asarray(y2_new[:index+1])) ############################################################ @@ -230,6 +226,6 @@ def add_discrete_function_cython(double[:] x1, double[:] y1, double[:] mp1, # the last value is again the end of the interval # only use the data that was actually filled - return (np.array(x_new[:index+1]), - np.array(y_new[:index+1]), - np.array(mp_new[:index+1])) + return (np.asarray(x_new[:index+1]), + np.asarray(y_new[:index+1]), + np.asarray(mp_new[:index+1])) diff --git a/pyspike/generic.py b/pyspike/generic.py index 4f278d2..41affcb 100644 --- a/pyspike/generic.py +++ b/pyspike/generic.py @@ -31,6 +31,69 @@ def _generic_profile_multi(spike_trains, pair_distance_func, indices=None): Returns: - The averaged multi-variate distance of all pairs """ + + def divide_and_conquer(pairs1, pairs2): + """ recursive calls by splitting the two lists in half. + """ + L1 = len(pairs1) + if L1 > 1: + dist_prof1 = divide_and_conquer(pairs1[:L1/2], pairs1[L1/2:]) + else: + dist_prof1 = pair_distance_func(spike_trains[pairs1[0][0]], + spike_trains[pairs1[0][1]]) + L2 = len(pairs2) + if L2 > 1: + dist_prof2 = divide_and_conquer(pairs2[:L2/2], pairs2[L2/2:]) + else: + dist_prof2 = pair_distance_func(spike_trains[pairs2[0][0]], + spike_trains[pairs2[0][1]]) + dist_prof1.add(dist_prof2) + return dist_prof1 + + if indices is None: + indices = np.arange(len(spike_trains)) + indices = np.array(indices) + # check validity of indices + assert (indices < len(spike_trains)).all() and (indices >= 0).all(), \ + "Invalid index list." + # generate a list of possible index pairs + pairs = [(indices[i], j) for i in range(len(indices)) + for j in indices[i+1:]] + + L = len(pairs) + if L > 1: + # recursive iteration through the list of pairs to get average profile + avrg_dist = divide_and_conquer(pairs[:len(pairs)/2], + pairs[len(pairs)/2:]) + else: + avrg_dist = pair_distance_func(spike_trains[pairs[0][0]], + spike_trains[pairs[0][1]]) + + return avrg_dist, L + + +############################################################ +# _generic_distance_multi +############################################################ +def _generic_distance_multi(spike_trains, pair_distance_func, + indices=None, interval=None): + """ Internal implementation detail, don't call this function directly, + use isi_distance_multi or spike_distance_multi instead. + + Computes the multi-variate distance for a set of spike-trains using the + pair_dist_func to compute pair-wise distances. That is it computes the + average distance of all pairs of spike-trains: + :math:`S(t) = 2/((N(N-1)) sum_{} D_{i,j}`, + where the sum goes over all pairs . + Args: + - spike_trains: list of spike trains + - pair_distance_func: function computing the distance of two spike trains + - indices: list of indices defining which spike trains to use, + if None all given spike trains are used (default=None) + Returns: + - The averaged multi-variate distance of all pairs + """ + if indices is None: indices = np.arange(len(spike_trains)) indices = np.array(indices) @@ -40,13 +103,13 @@ def _generic_profile_multi(spike_trains, pair_distance_func, indices=None): # generate a list of possible index pairs pairs = [(indices[i], j) for i in range(len(indices)) for j in indices[i+1:]] - # start with first pair - (i, j) = pairs[0] - average_dist = pair_distance_func(spike_trains[i], spike_trains[j]) - for (i, j) in pairs[1:]: - current_dist = pair_distance_func(spike_trains[i], spike_trains[j]) - average_dist.add(current_dist) # add to the average - return average_dist, len(pairs) + + avrg_dist = 0.0 + for (i, j) in pairs: + avrg_dist += pair_distance_func(spike_trains[i], spike_trains[j], + interval) + + return avrg_dist/len(pairs) ############################################################ diff --git a/pyspike/isi_distance.py b/pyspike/isi_distance.py index aeab0df..21df561 100644 --- a/pyspike/isi_distance.py +++ b/pyspike/isi_distance.py @@ -3,7 +3,8 @@ # Distributed under the BSD License from pyspike import PieceWiseConstFunc -from pyspike.generic import _generic_profile_multi, _generic_distance_matrix +from pyspike.generic import _generic_profile_multi, _generic_distance_multi, \ + _generic_distance_matrix ############################################################ @@ -112,7 +113,8 @@ def isi_distance_multi(spike_trains, indices=None, interval=None): :returns: The time-averaged multivariate ISI distance :math:`D_I` :rtype: double """ - return isi_profile_multi(spike_trains, indices).avrg(interval) + return _generic_distance_multi(spike_trains, isi_distance, indices, + interval) ############################################################ 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) ############################################################ diff --git a/test/test_distance.py b/test/test_distance.py index 19da35f..e45ac16 100644 --- a/test/test_distance.py +++ b/test/test_distance.py @@ -196,7 +196,7 @@ def test_spike_sync(): 0.4, decimal=16) -def check_multi_profile(profile_func, profile_func_multi): +def check_multi_profile(profile_func, profile_func_multi, dist_func_multi): # generate spike trains: 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) @@ -213,10 +213,14 @@ def check_multi_profile(profile_func, profile_func_multi): f_multi = profile_func_multi(spike_trains, [0, 1]) assert f_multi.almost_equal(f12, decimal=14) + d = dist_func_multi(spike_trains, [0, 1]) + assert_equal(f_multi.avrg(), d) f_multi1 = profile_func_multi(spike_trains, [1, 2, 3]) f_multi2 = profile_func_multi(spike_trains[1:]) assert f_multi1.almost_equal(f_multi2, decimal=14) + d = dist_func_multi(spike_trains, [1, 2, 3]) + assert_almost_equal(f_multi1.avrg(), d, decimal=14) f = copy(f12) f.add(f13) @@ -224,6 +228,8 @@ def check_multi_profile(profile_func, profile_func_multi): f.mul_scalar(1.0/3) f_multi = profile_func_multi(spike_trains, [0, 1, 2]) assert f_multi.almost_equal(f, decimal=14) + d = dist_func_multi(spike_trains, [0, 1, 2]) + assert_almost_equal(f_multi.avrg(), d, decimal=14) f.mul_scalar(3) # revert above normalization f.add(f14) @@ -235,11 +241,13 @@ def check_multi_profile(profile_func, profile_func_multi): def test_multi_isi(): - check_multi_profile(spk.isi_profile, spk.isi_profile_multi) + check_multi_profile(spk.isi_profile, spk.isi_profile_multi, + spk.isi_distance_multi) def test_multi_spike(): - check_multi_profile(spk.spike_profile, spk.spike_profile_multi) + check_multi_profile(spk.spike_profile, spk.spike_profile_multi, + spk.spike_distance_multi) def test_multi_spike_sync(): -- cgit v1.2.3 From b09561705ab9c67c93a384248f7c3bc9ad5bdd32 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Wed, 3 Feb 2016 13:07:26 +0100 Subject: fixed spike-sync bug fixed ugly bugs in code for computing multi-variate spike sync profile and multi-variate spike sync value. --- pyspike/DiscreteFunc.py | 9 ++++---- pyspike/cython/cython_add.pyx | 33 ++++++++++++++------------- pyspike/cython/cython_distances.pyx | 5 ---- pyspike/cython/python_backend.py | 33 +++++++++++++++------------ pyspike/spike_sync.py | 10 ++++++-- test/numeric/test_regression_random_spikes.py | 27 +++++++++++++++------- test/test_distance.py | 33 +++++++++++++++++++++++++++ 7 files changed, 99 insertions(+), 51 deletions(-) (limited to 'pyspike/cython/cython_add.pyx') diff --git a/pyspike/DiscreteFunc.py b/pyspike/DiscreteFunc.py index fe97bc2..caad290 100644 --- a/pyspike/DiscreteFunc.py +++ b/pyspike/DiscreteFunc.py @@ -170,10 +170,6 @@ expected." start_ind, end_ind = get_indices(ival) value += np.sum(self.y[start_ind:end_ind]) multiplicity += np.sum(self.mp[start_ind:end_ind]) - if multiplicity == 0.0: - # empty profile, return spike sync of 1 - value = 1.0 - multiplicity = 1.0 return (value, multiplicity) def avrg(self, interval=None, normalize=True): @@ -190,7 +186,10 @@ expected." """ val, mp = self.integral(interval) if normalize: - return val/mp + if mp > 0: + return val/mp + else: + return 1.0 else: return val diff --git a/pyspike/cython/cython_add.pyx b/pyspike/cython/cython_add.pyx index 8da1e53..25f1181 100644 --- a/pyspike/cython/cython_add.pyx +++ b/pyspike/cython/cython_add.pyx @@ -182,8 +182,8 @@ def add_discrete_function_cython(double[:] x1, double[:] y1, double[:] mp1, cdef int index1 = 0 cdef int index2 = 0 cdef int index = 0 - cdef int N1 = len(y1) - cdef int N2 = len(y2) + cdef int N1 = len(y1)-1 + cdef int N2 = len(y2)-1 x_new[0] = x1[0] while (index1+1 < N1) and (index2+1 < N2): if x1[index1+1] < x2[index2+1]: @@ -206,20 +206,21 @@ def add_discrete_function_cython(double[:] x1, double[:] y1, double[:] mp1, y_new[index] = y1[index1] + y2[index2] mp_new[index] = mp1[index1] + mp2[index2] # one array reached the end -> copy the contents of the other to the end - if index1+1 < len(y1): - x_new[index+1:index+1+len(x1)-index1-1] = x1[index1+1:] - y_new[index+1:index+1+len(x1)-index1-1] = y1[index1+1:] - mp_new[index+1:index+1+len(x1)-index1-1] = mp1[index1+1:] - index += len(x1)-index1-1 - elif index2+1 < len(y2): - x_new[index+1:index+1+len(x2)-index2-1] = x2[index2+1:] - y_new[index+1:index+1+len(x2)-index2-1] = y2[index2+1:] - mp_new[index+1:index+1+len(x2)-index2-1] = mp2[index2+1:] - index += len(x2)-index2-1 - # else: # both arrays reached the end simultaneously - # x_new[index+1] = x1[-1] - # y_new[index+1] = y1[-1] + y2[-1] - # mp_new[index+1] = mp1[-1] + mp2[-1] + if index1+1 < N1: + x_new[index+1:index+1+N1-index1] = x1[index1+1:] + y_new[index+1:index+1+N1-index1] = y1[index1+1:] + mp_new[index+1:index+1+N1-index1] = mp1[index1+1:] + index += N1-index1 + elif index2+1 < N2: + x_new[index+1:index+1+N2-index2] = x2[index2+1:] + y_new[index+1:index+1+N2-index2] = y2[index2+1:] + mp_new[index+1:index+1+N2-index2] = mp2[index2+1:] + index += N2-index2 + else: # both arrays reached the end simultaneously + x_new[index+1] = x1[index1+1] + y_new[index+1] = y1[index1+1] + y2[index2+1] + mp_new[index+1] = mp1[index1+1] + mp2[index2+1] + index += 1 y_new[0] = y_new[1] mp_new[0] = mp_new[1] diff --git a/pyspike/cython/cython_distances.pyx b/pyspike/cython/cython_distances.pyx index 7dc9cb9..ac5f226 100644 --- a/pyspike/cython/cython_distances.pyx +++ b/pyspike/cython/cython_distances.pyx @@ -428,9 +428,4 @@ def coincidence_value_cython(double[:] spikes1, double[:] spikes2, mp += 2 coinc += 2 - if coinc == 0 and mp == 0: - # empty spike trains -> spike sync = 1 by definition - coinc = 1 - mp = 1 - return coinc, mp diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py index a1e3a65..6b7209a 100644 --- a/pyspike/cython/python_backend.py +++ b/pyspike/cython/python_backend.py @@ -560,7 +560,9 @@ def add_discrete_function_python(x1, y1, mp1, x2, y2, mp2): index1 = 0 index2 = 0 index = 0 - while (index1+1 < len(y1)) and (index2+1 < len(y2)): + N1 = len(x1)-1 + N2 = len(x2)-1 + while (index1+1 < N1) and (index2+1 < N2): if x1[index1+1] < x2[index2+1]: index1 += 1 index += 1 @@ -581,20 +583,21 @@ def add_discrete_function_python(x1, y1, mp1, x2, y2, mp2): y_new[index] = y1[index1] + y2[index2] mp_new[index] = mp1[index1] + mp2[index2] # one array reached the end -> copy the contents of the other to the end - if index1+1 < len(y1): - x_new[index+1:index+1+len(x1)-index1-1] = x1[index1+1:] - y_new[index+1:index+1+len(x1)-index1-1] = y1[index1+1:] - mp_new[index+1:index+1+len(x1)-index1-1] = mp1[index1+1:] - index += len(x1)-index1-1 - elif index2+1 < len(y2): - x_new[index+1:index+1+len(x2)-index2-1] = x2[index2+1:] - y_new[index+1:index+1+len(x2)-index2-1] = y2[index2+1:] - mp_new[index+1:index+1+len(x2)-index2-1] = mp2[index2+1:] - index += len(x2)-index2-1 - # else: # both arrays reached the end simultaneously - # x_new[index+1] = x1[-1] - # y_new[index+1] = y1[-1] + y2[-1] - # mp_new[index+1] = mp1[-1] + mp2[-1] + if index1+1 < N1: + x_new[index+1:index+1+N1-index1] = x1[index1+1:] + y_new[index+1:index+1+N1-index1] = y1[index1+1:] + mp_new[index+1:index+1+N1-index1] = mp1[index1+1:] + index += N1-index1 + elif index2+1 < N2: + x_new[index+1:index+1+N2-index2] = x2[index2+1:] + y_new[index+1:index+1+N2-index2] = y2[index2+1:] + mp_new[index+1:index+1+N2-index2] = mp2[index2+1:] + index += N2-index2 + else: # both arrays reached the end simultaneously + x_new[index+1] = x1[-1] + y_new[index+1] = y1[-1] + y2[-1] + mp_new[index+1] = mp1[-1] + mp2[-1] + index += 1 y_new[0] = y_new[1] mp_new[0] = mp_new[1] diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py index 3dc29ff..802b98c 100644 --- a/pyspike/spike_sync.py +++ b/pyspike/spike_sync.py @@ -119,7 +119,10 @@ def spike_sync(spike_train1, spike_train2, interval=None, max_tau=None): """ c, mp = _spike_sync_values(spike_train1, spike_train2, interval, max_tau) - return 1.0*c/mp + if mp == 0: + return 1.0 + else: + return 1.0*c/mp ############################################################ @@ -187,7 +190,10 @@ def spike_sync_multi(spike_trains, indices=None, interval=None, max_tau=None): coincidence += c mp += m - return coincidence/mp + if mp == 0.0: + return 1.0 + else: + return coincidence/mp ############################################################ diff --git a/test/numeric/test_regression_random_spikes.py b/test/numeric/test_regression_random_spikes.py index 6156bb4..73d39a6 100644 --- a/test/numeric/test_regression_random_spikes.py +++ b/test/numeric/test_regression_random_spikes.py @@ -35,7 +35,9 @@ def test_regression_random(): spike = spk.spike_distance_multi(spike_trains) spike_prof = spk.spike_profile_multi(spike_trains).avrg() - # spike_sync = spk.spike_sync_multi(spike_trains) + + spike_sync = spk.spike_sync_multi(spike_trains) + spike_sync_prof = spk.spike_sync_profile_multi(spike_trains).avrg() assert_almost_equal(isi, results_cSPIKY[i][0], decimal=14, err_msg="Index: %d, ISI" % i) @@ -47,6 +49,9 @@ def test_regression_random(): assert_almost_equal(spike_prof, results_cSPIKY[i][1], decimal=14, err_msg="Index: %d, SPIKE" % i) + assert_almost_equal(spike_sync, spike_sync_prof, decimal=14, + err_msg="Index: %d, SPIKE-Sync" % i) + def check_regression_dataset(spike_file="benchmark.mat", spikes_name="spikes", @@ -109,19 +114,25 @@ def check_single_spike_train_set(index): spike_train_data = spike_train_sets[index] spike_trains = [] + N = 0 for spikes in spike_train_data[0]: - print("Spikes:", spikes.flatten()) - spike_trains.append(spk.SpikeTrain(spikes.flatten(), 100.0)) + N += len(spikes.flatten()) + print("Spikes:", len(spikes.flatten())) + spikes_array = spikes.flatten() + if len(spikes_array > 0) and (spikes_array[-1] > 100.0): + spikes_array[-1] = 100.0 + spike_trains.append(spk.SpikeTrain(spikes_array, 100.0)) + print(spike_trains[-1].spikes) - print(spk.spike_distance_multi(spike_trains)) + print(N) - print(results_cSPIKY[index][1]) + print(spk.spike_sync_multi(spike_trains)) - print(spike_trains[1].spikes) + print(spk.spike_sync_profile_multi(spike_trains).integral()) if __name__ == "__main__": - test_regression_random() + # test_regression_random() # check_regression_dataset() - # check_single_spike_train_set(7633) + check_single_spike_train_set(4) diff --git a/test/test_distance.py b/test/test_distance.py index 083d8a3..fe09f34 100644 --- a/test/test_distance.py +++ b/test/test_distance.py @@ -217,6 +217,28 @@ def test_spike_sync(): assert_almost_equal(spk.spike_sync(spikes1, spikes2), 0.4, decimal=16) + spikes1 = SpikeTrain([1.0, 2.0, 4.0], 4.0) + spikes2 = SpikeTrain([3.8], 4.0) + spikes3 = SpikeTrain([3.9, ], 4.0) + + expected_x = np.array([0.0, 1.0, 2.0, 3.8, 4.0, 4.0]) + expected_y = np.array([0.0, 0.0, 0.0, 1.0, 1.0, 1.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) + + f2 = spk.spike_sync_profile(spikes2, spikes3) + + i1 = f.integral() + i2 = f2.integral() + f.add(f2) + i12 = f.integral() + + assert_equal(i1[0]+i2[0], i12[0]) + assert_equal(i1[1]+i2[1], i12[1]) + def check_multi_profile(profile_func, profile_func_multi, dist_func_multi): # generate spike trains: @@ -313,6 +335,17 @@ def test_multi_spike_sync(): assert_equal(np.sum(f.y[1:-1]), 39932) assert_equal(np.sum(f.mp[1:-1]), 85554) + # example with 2 empty spike trains + sts = [] + sts.append(SpikeTrain([1, 9], [0, 10])) + sts.append(SpikeTrain([1, 3], [0, 10])) + sts.append(SpikeTrain([], [0, 10])) + sts.append(SpikeTrain([], [0, 10])) + + assert_almost_equal(spk.spike_sync_multi(sts), 1.0/6.0, decimal=15) + assert_almost_equal(spk.spike_sync_profile_multi(sts).avrg(), 1.0/6.0, + decimal=15) + def check_dist_matrix(dist_func, dist_matrix_func): # generate spike trains: -- cgit v1.2.3