From 748b4c5fcc2f41d8d414e3714b66dd6c237c8ff7 Mon Sep 17 00:00:00 2001 From: ImmanuelSamuel Date: Mon, 10 Aug 2015 17:11:32 -0400 Subject: Indices for matrix slicing needs to be integers Dividing the indices by 2 gives a float which cannot be used for slicing the matrix. Changed the type of the indices to integer to ensure that matrix slicing does not return error. --- pyspike/generic.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) (limited to 'pyspike/generic.py') diff --git a/pyspike/generic.py b/pyspike/generic.py index 41affcb..4d70e65 100644 --- a/pyspike/generic.py +++ b/pyspike/generic.py @@ -37,13 +37,13 @@ def _generic_profile_multi(spike_trains, pair_distance_func, indices=None): """ L1 = len(pairs1) if L1 > 1: - dist_prof1 = divide_and_conquer(pairs1[:L1/2], pairs1[L1/2:]) + dist_prof1 = divide_and_conquer(pairs1[:int(L1/2)], pairs1[int(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:]) + dist_prof2 = divide_and_conquer(pairs2[:int(L2/2)], pairs2[int(L2/2):]) else: dist_prof2 = pair_distance_func(spike_trains[pairs2[0][0]], spike_trains[pairs2[0][1]]) @@ -63,8 +63,8 @@ def _generic_profile_multi(spike_trains, pair_distance_func, indices=None): 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:]) + avrg_dist = divide_and_conquer(pairs[:int(len(pairs)/2]), + pairs[int(len(pairs)/2):]) else: avrg_dist = pair_distance_func(spike_trains[pairs[0][0]], spike_trains[pairs[0][1]]) -- cgit v1.2.3 From d95d021f91f165b80be94b11c7bb526a055a94ef Mon Sep 17 00:00:00 2001 From: ImmanuelSamuel Date: Mon, 10 Aug 2015 17:29:08 -0400 Subject: Syntax correction --- pyspike/generic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'pyspike/generic.py') diff --git a/pyspike/generic.py b/pyspike/generic.py index 4d70e65..9acba26 100644 --- a/pyspike/generic.py +++ b/pyspike/generic.py @@ -63,7 +63,7 @@ def _generic_profile_multi(spike_trains, pair_distance_func, indices=None): L = len(pairs) if L > 1: # recursive iteration through the list of pairs to get average profile - avrg_dist = divide_and_conquer(pairs[:int(len(pairs)/2]), + avrg_dist = divide_and_conquer(pairs[:int(len(pairs)/2)], pairs[int(len(pairs)/2):]) else: avrg_dist = pair_distance_func(spike_trains[pairs[0][0]], -- cgit v1.2.3 From d97d1f798ea1b4fe90c0e3b6828a219c90996271 Mon Sep 17 00:00:00 2001 From: ImmanuelSamuel Date: Mon, 10 Aug 2015 18:39:06 -0400 Subject: Used // instead of type conversion --- pyspike/generic.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) (limited to 'pyspike/generic.py') diff --git a/pyspike/generic.py b/pyspike/generic.py index 9acba26..723c90d 100644 --- a/pyspike/generic.py +++ b/pyspike/generic.py @@ -37,13 +37,13 @@ def _generic_profile_multi(spike_trains, pair_distance_func, indices=None): """ L1 = len(pairs1) if L1 > 1: - dist_prof1 = divide_and_conquer(pairs1[:int(L1/2)], pairs1[int(L1/2):]) + dist_prof1 = divide_and_conquer(pairs1[:L1//2], pairs1[int(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[:int(L2/2)], pairs2[int(L2/2):]) + dist_prof2 = divide_and_conquer(pairs2[:L2//2], pairs2[int(L2/2):]) else: dist_prof2 = pair_distance_func(spike_trains[pairs2[0][0]], spike_trains[pairs2[0][1]]) @@ -63,8 +63,8 @@ def _generic_profile_multi(spike_trains, pair_distance_func, indices=None): L = len(pairs) if L > 1: # recursive iteration through the list of pairs to get average profile - avrg_dist = divide_and_conquer(pairs[:int(len(pairs)/2)], - pairs[int(len(pairs)/2):]) + 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]]) -- cgit v1.2.3 From 3c6cdfe943e54880325c5ee81ab0d96a6bb101de Mon Sep 17 00:00:00 2001 From: ImmanuelSamuel Date: Mon, 10 Aug 2015 19:06:50 -0400 Subject: Added more // --- pyspike/generic.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'pyspike/generic.py') diff --git a/pyspike/generic.py b/pyspike/generic.py index 723c90d..2df34f1 100644 --- a/pyspike/generic.py +++ b/pyspike/generic.py @@ -37,13 +37,13 @@ def _generic_profile_multi(spike_trains, pair_distance_func, indices=None): """ L1 = len(pairs1) if L1 > 1: - dist_prof1 = divide_and_conquer(pairs1[:L1//2], pairs1[int(L1/2):]) + dist_prof1 = divide_and_conquer(pairs1[:L1//2], pairs1[int(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[int(L2/2):]) + dist_prof2 = divide_and_conquer(pairs2[:L2//2], pairs2[int(L2//2):]) else: dist_prof2 = pair_distance_func(spike_trains[pairs2[0][0]], spike_trains[pairs2[0][1]]) -- cgit v1.2.3 From 29e50478dcfc31ce04c4343fa585463abe96caae Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Wed, 12 Aug 2015 18:42:54 +0200 Subject: new spike delay asymmetry measures added first version of spike delay asymmetry functions. still incomplete and untested. --- pyspike/__init__.py | 3 +- pyspike/directionality/__init__.py | 12 ++ pyspike/directionality/cython/__init__.py | 0 .../cython/cython_directionality.pyx | 177 +++++++++++++++++ pyspike/directionality/spike_delay_asymmetry.py | 212 +++++++++++++++++++++ pyspike/generic.py | 6 +- setup.py | 28 ++- 7 files changed, 427 insertions(+), 11 deletions(-) create mode 100644 pyspike/directionality/__init__.py create mode 100644 pyspike/directionality/cython/__init__.py create mode 100644 pyspike/directionality/cython/cython_directionality.pyx create mode 100644 pyspike/directionality/spike_delay_asymmetry.py (limited to 'pyspike/generic.py') diff --git a/pyspike/__init__.py b/pyspike/__init__.py index 2060f73..8d92ea4 100644 --- a/pyspike/__init__.py +++ b/pyspike/__init__.py @@ -6,7 +6,7 @@ Distributed under the BSD License __all__ = ["isi_distance", "spike_distance", "spike_sync", "psth", "spikes", "SpikeTrain", "PieceWiseConstFunc", "PieceWiseLinFunc", - "DiscreteFunc"] + "DiscreteFunc", "directionality"] from PieceWiseConstFunc import PieceWiseConstFunc from PieceWiseLinFunc import PieceWiseLinFunc @@ -24,6 +24,7 @@ from psth import psth from spikes import load_spike_trains_from_txt, spike_train_from_string, \ merge_spike_trains, generate_poisson_spikes +import directionality as drct # define the __version__ following # http://stackoverflow.com/questions/17583443 diff --git a/pyspike/directionality/__init__.py b/pyspike/directionality/__init__.py new file mode 100644 index 0000000..e6de1de --- /dev/null +++ b/pyspike/directionality/__init__.py @@ -0,0 +1,12 @@ +""" +Copyright 2015, Mario Mulansky + +Distributed under the BSD License +""" + +__all__ = ["spike_delay_asymmetry"] + +from spike_delay_asymmetry import spike_delay_asymmetry_profile, \ + spike_delay_asymmetry, spike_delay_asymmetry_profile_multi, \ + spike_delay_asymmetry_matrix, optimal_asymmetry_order, \ + optimal_asymmetry_order_from_D, reorder_asymmetry_matrix diff --git a/pyspike/directionality/cython/__init__.py b/pyspike/directionality/cython/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/pyspike/directionality/cython/cython_directionality.pyx b/pyspike/directionality/cython/cython_directionality.pyx new file mode 100644 index 0000000..f5ea752 --- /dev/null +++ b/pyspike/directionality/cython/cython_directionality.pyx @@ -0,0 +1,177 @@ +#cython: boundscheck=False +#cython: wraparound=False +#cython: cdivision=True + +""" +cython_directionality.pyx + +cython implementation of the spike delay asymmetry measures + +Copyright 2015, Mario Mulansky + +Distributed under the BSD License + +""" + +""" +To test whether things can be optimized: remove all yellow stuff +in the html output:: + + cython -a cython_directionality.pyx + +which gives:: + + cython_directionality.html + +""" + +import numpy as np +cimport numpy as np + +from libc.math cimport fabs +from libc.math cimport fmax +from libc.math cimport fmin + +# from pyspike.cython.cython_distances cimport get_tau + +DTYPE = np.float +ctypedef np.float_t DTYPE_t + + +############################################################ +# get_tau +############################################################ +cdef inline double get_tau(double[:] spikes1, double[:] spikes2, + int i, int j, double interval, double max_tau): + cdef double m = interval # use interval length as initial tau + cdef int N1 = spikes1.shape[0]-1 # len(spikes1)-1 + cdef int N2 = spikes2.shape[0]-1 # len(spikes2)-1 + if i < N1 and i > -1: + m = fmin(m, spikes1[i+1]-spikes1[i]) + if j < N2 and j > -1: + m = fmin(m, spikes2[j+1]-spikes2[j]) + if i > 0: + m = fmin(m, spikes1[i]-spikes1[i-1]) + if j > 0: + m = fmin(m, spikes2[j]-spikes2[j-1]) + m *= 0.5 + if max_tau > 0.0: + m = fmin(m, max_tau) + return m + + +############################################################ +# spike_delay_asymmetry_profile_cython +############################################################ +def spike_delay_asymmetry_profile_cython(double[:] spikes1, double[:] spikes2, + double t_start, double t_end, + double max_tau): + + cdef int N1 = len(spikes1) + cdef int N2 = len(spikes2) + cdef int i = -1 + cdef int j = -1 + cdef int n = 0 + cdef double[:] st = np.zeros(N1 + N2 + 2) # spike times + cdef double[:] a = np.zeros(N1 + N2 + 2) # asymmetry values + cdef double[:] mp = np.ones(N1 + N2 + 2) # multiplicity + cdef double interval = t_end - t_start + cdef double tau + while i + j < N1 + N2 - 2: + if (i < N1-1) and (j == N2-1 or spikes1[i+1] < spikes2[j+1]): + i += 1 + n += 1 + tau = get_tau(spikes1, spikes2, i, j, interval, max_tau) + st[n] = spikes1[i] + if j > -1 and spikes1[i]-spikes2[j] < tau: + # coincidence between the current spike and the previous spike + # spike from spike train 1 after spike train 2 + # both get marked with -1 + a[n] = -1 + a[n-1] = -1 + elif (j < N2-1) and (i == N1-1 or spikes1[i+1] > spikes2[j+1]): + j += 1 + n += 1 + tau = get_tau(spikes1, spikes2, i, j, interval, max_tau) + st[n] = spikes2[j] + if i > -1 and spikes2[j]-spikes1[i] < tau: + # coincidence between the current spike and the previous spike + # spike from spike train 1 before spike train 2 + # both get marked with 1 + a[n] = 1 + a[n-1] = 1 + else: # spikes1[i+1] = spikes2[j+1] + # advance in both spike trains + j += 1 + i += 1 + n += 1 + # add only one event with zero asymmetry value and multiplicity 2 + st[n] = spikes1[i] + a[n] = 0 + mp[n] = 2 + + st = st[:n+2] + a = a[:n+2] + mp = mp[:n+2] + + st[0] = t_start + st[len(st)-1] = t_end + if N1 + N2 > 0: + a[0] = a[1] + a[len(a)-1] = a[len(a)-2] + mp[0] = mp[1] + mp[len(mp)-1] = mp[len(mp)-2] + else: + a[0] = 1 + a[1] = 1 + + return st, a, mp + + + +############################################################ +# spike_delay_asymmetry_cython +############################################################ +def spike_delay_asymmetry_cython(double[:] spikes1, double[:] spikes2, + double t_start, double t_end, double max_tau): + + cdef int N1 = len(spikes1) + cdef int N2 = len(spikes2) + cdef int i = -1 + cdef int j = -1 + cdef int asym = 0 + cdef int mp = 0 + cdef double interval = t_end - t_start + cdef double tau + while i + j < N1 + N2 - 2: + if (i < N1-1) and (j == N2-1 or spikes1[i+1] < spikes2[j+1]): + i += 1 + mp += 1 + tau = get_tau(spikes1, spikes2, i, j, interval, max_tau) + if j > -1 and spikes1[i]-spikes2[j] < tau: + # coincidence between the current spike and the previous spike + # spike in spike train 2 appeared before spike in spike train 1 + # mark with -1 + asym -= 1 + elif (j < N2-1) and (i == N1-1 or spikes1[i+1] > spikes2[j+1]): + j += 1 + mp += 1 + tau = get_tau(spikes1, spikes2, i, j, interval, max_tau) + if i > -1 and spikes2[j]-spikes1[i] < tau: + # coincidence between the current spike and the previous spike + # spike in spike train 1 appeared before spike in spike train 2 + # mark with +1 + asym += 1 + else: # spikes1[i+1] = spikes2[j+1] + # advance in both spike trains + j += 1 + i += 1 + # add only one event with multiplicity 2, but no asymmetry counting + mp += 2 + + if asym == 0 and mp == 0: + # empty spike trains -> spike sync = 1 by definition + asym = 1 + mp = 1 + + return asym, mp diff --git a/pyspike/directionality/spike_delay_asymmetry.py b/pyspike/directionality/spike_delay_asymmetry.py new file mode 100644 index 0000000..7d59601 --- /dev/null +++ b/pyspike/directionality/spike_delay_asymmetry.py @@ -0,0 +1,212 @@ +# Module containing functions to compute multivariate spike delay asymmetry +# Copyright 2015, Mario Mulansky +# Distributed under the BSD License + +import numpy as np +from math import exp +from functools import partial +# import pyspike +from pyspike import DiscreteFunc +from pyspike.generic import _generic_profile_multi + + +############################################################ +# spike_delay_asymmetry_profile +############################################################ +def spike_delay_asymmetry_profile(spike_train1, spike_train2, max_tau=None): + """ Computes the spike delay asymmetry profile A(t) of the two given + spike trains. Returns the profile as a DiscreteFunction object. + + :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 are not defined on the same interval!" + assert spike_train1.t_end == spike_train2.t_end, \ + "Given spike trains are not defined on the same interval!" + + # cython implementation + try: + from cython.cython_directionality import \ + spike_delay_asymmetry_profile_cython as \ + spike_delay_asymmetry_profile_impl + except ImportError: + raise NotImplementedError() +# if not(pyspike.disable_backend_warning): +# print("Warning: spike_distance_cython not found. Make sure that \ +# PySpike is installed by running\n 'python setup.py build_ext --inplace'!\n \ +# Falling back to slow python backend.") +# # use python backend +# from cython.python_backend import coincidence_python \ +# as coincidence_profile_impl + + if max_tau is None: + max_tau = 0.0 + + times, coincidences, multiplicity \ + = spike_delay_asymmetry_profile_impl(spike_train1.spikes, + spike_train2.spikes, + spike_train1.t_start, + spike_train1.t_end, + max_tau) + + return DiscreteFunc(times, coincidences, multiplicity) + + +############################################################ +# spike_delay_asymmetry +############################################################ +def spike_delay_asymmetry(spike_train1, spike_train2, + interval=None, max_tau=None): + """ Computes the overall spike delay asymmetry value for two spike trains. + """ + if interval is None: + # distance over the whole interval is requested: use specific function + # for optimal performance + try: + from cython.cython_directionality import \ + spike_delay_asymmetry_cython as spike_delay_impl + if max_tau is None: + max_tau = 0.0 + c, mp = spike_delay_impl(spike_train1.spikes, + spike_train2.spikes, + spike_train1.t_start, + spike_train1.t_end, + max_tau) + return c + except ImportError: + # Cython backend not available: fall back to profile averaging + raise NotImplementedError() + # return spike_sync_profile(spike_train1, spike_train2, + # max_tau).integral(interval) + else: + # some specific interval is provided: not yet implemented + raise NotImplementedError() + + +############################################################ +# spike_delay_asymmetry_profile_multi +############################################################ +def spike_delay_asymmetry_profile_multi(spike_trains, indices=None, + max_tau=None): + """ Computes the multi-variate spike delay asymmetry profile for a set of + spike trains. For each spike in the set of spike trains, the multi-variate + profile is defined as the sum of asymmetry values divided by the number of + 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 :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 + :param max_tau: Maximum coincidence window size. If 0 or `None`, the + coincidence window has no upper bound. + :returns: The multi-variate spike sync profile :math:`(t)` + :rtype: :class:`pyspike.function.DiscreteFunction` + + """ + prof_func = partial(spike_delay_asymmetry_profile, max_tau=max_tau) + average_prof, M = _generic_profile_multi(spike_trains, prof_func, + indices) + # average_dist.mul_scalar(1.0/M) # no normalization here! + return average_prof + + +############################################################ +# spike_delay_asymmetry_matrix +############################################################ +def spike_delay_asymmetry_matrix(spike_trains, indices=None, + interval=None, max_tau=None): + """ Computes the spike delay asymmetry matrix for the given spike trains. + """ + 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:]] + + distance_matrix = np.zeros((len(indices), len(indices))) + for i, j in pairs: + d = spike_delay_asymmetry(spike_trains[i], spike_trains[j], + interval, max_tau=max_tau) + distance_matrix[i, j] = d + distance_matrix[j, i] = -d + return distance_matrix + + +############################################################ +# optimal_asymmetry_order_from_D +############################################################ +def optimal_asymmetry_order_from_D(D, full_output=False): + """ finds the best sorting via simulated annealing. + Returns the optimal permutation p and A value. + Internal function, don't call directly! Use optimal_asymmetry_order + instead. + """ + N = len(D) + A = np.sum(np.triu(D, 0)) + + p = np.arange(N) + + T = 2*np.max(D) # starting temperature + T_end = 1E-5 * T # final temperature + alpha = 0.9 # cooling factor + total_iter = 0 + while T > T_end: + iterations = 0 + succ_iter = 0 + while iterations < 100*N and succ_iter < 10*N: + # exchange two rows and cols + ind1 = np.random.randint(N-1) + delta_A = -2*D[p[ind1], p[ind1+1]] + if delta_A > 0.0 or exp(delta_A/T) > np.random.random(): + # swap indices + p[ind1], p[ind1+1] = p[ind1+1], p[ind1] + A += delta_A + succ_iter += 1 + iterations += 1 + total_iter += iterations + T *= alpha # cool down + if succ_iter == 0: + break + if full_output: + return p, A, total_iter + else: + return p, A + + +############################################################ +# _optimal_asymmetry_order +############################################################ +def optimal_asymmetry_order(spike_trains, indices=None, interval=None, + max_tau=None, full_output=False): + """ finds the best sorting of the given spike trains via simulated + annealing. + Returns the optimal permutation p and A value. + """ + D = spike_delay_asymmetry_matrix(spike_trains, indices, interval, max_tau) + return optimal_asymmetry_order_from_D(D, full_output) + + +############################################################ +# reorder_asymmetry_matrix +############################################################ +def reorder_asymmetry_matrix(D, p): + N = len(D) + D_p = np.empty_like(D) + for n in xrange(N): + for m in xrange(N): + D_p[n, m] = D[p[n], p[m]] + return D_p diff --git a/pyspike/generic.py b/pyspike/generic.py index 2df34f1..515cbf4 100644 --- a/pyspike/generic.py +++ b/pyspike/generic.py @@ -37,13 +37,15 @@ def _generic_profile_multi(spike_trains, pair_distance_func, indices=None): """ L1 = len(pairs1) if L1 > 1: - dist_prof1 = divide_and_conquer(pairs1[:L1//2], pairs1[int(L1//2):]) + dist_prof1 = divide_and_conquer(pairs1[:L1//2], + pairs1[int(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[int(L2//2):]) + dist_prof2 = divide_and_conquer(pairs2[:L2//2], + pairs2[int(L2//2):]) else: dist_prof2 = pair_distance_func(spike_trains[pairs2[0][0]], spike_trains[pairs2[0][1]]) diff --git a/setup.py b/setup.py index d853cdf..960c684 100644 --- a/setup.py +++ b/setup.py @@ -23,7 +23,8 @@ else: if os.path.isfile("pyspike/cython/cython_add.c") and \ os.path.isfile("pyspike/cython/cython_profiles.c") and \ - os.path.isfile("pyspike/cython/cython_distances.c"): + os.path.isfile("pyspike/cython/cython_distances.c") and \ + os.path.isfile("pyspike/directionality/cython/cython_directionality.c"): use_c = True else: use_c = False @@ -33,16 +34,26 @@ ext_modules = [] if use_cython: # Cython is available, compile .pyx -> .c ext_modules += [ - Extension("pyspike.cython.cython_add", ["pyspike/cython/cython_add.pyx"]), - Extension("pyspike.cython.cython_profiles", ["pyspike/cython/cython_profiles.pyx"]), - Extension("pyspike.cython.cython_distances", ["pyspike/cython/cython_distances.pyx"]), + Extension("pyspike.cython.cython_add", + ["pyspike/cython/cython_add.pyx"]), + Extension("pyspike.cython.cython_profiles", + ["pyspike/cython/cython_profiles.pyx"]), + Extension("pyspike.cython.cython_distances", + ["pyspike/cython/cython_distances.pyx"]), + Extension("pyspike.directionality.cython.cython_directionality", + ["pyspike/directionality/cython/cython_directionality.pyx"]) ] cmdclass.update({'build_ext': build_ext}) elif use_c: # c files are there, compile to binaries ext_modules += [ - Extension("pyspike.cython.cython_add", ["pyspike/cython/cython_add.c"]), - Extension("pyspike.cython.cython_profiles", ["pyspike/cython/cython_profiles.c"]), - Extension("pyspike.cython.cython_distances", ["pyspike/cython/cython_distances.c"]), + Extension("pyspike.cython.cython_add", + ["pyspike/cython/cython_add.c"]), + Extension("pyspike.cython.cython_profiles", + ["pyspike/cython/cython_profiles.c"]), + Extension("pyspike.cython.cython_distances", + ["pyspike/cython/cython_distances.c"]), + Extension("pyspike.directionality.cython.cython_directionality", + ["pyspike/directionality/cython/cython_directionality.c"]) ] # neither cython nor c files available -> automatic fall-back to python backend @@ -81,7 +92,8 @@ train similarity', ], package_data={ 'pyspike': ['cython/cython_add.c', 'cython/cython_profiles.c', - 'cython_distances.c'], + 'cython/cython_distances.c', + 'directionality/cython/cython_directionality.c'], 'test': ['Spike_testdata.txt'] } ) -- cgit v1.2.3 From be74318a2b269ec0c1e16981e7286679746f1a49 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Mon, 17 Aug 2015 11:57:53 +0200 Subject: fix #15 add test case and fix for Issue #15 closes #15 --- pyspike/generic.py | 7 +-- test/regression/test_regression_11.py | 28 ----------- test/test_regression/test_regression_15.py | 78 ++++++++++++++++++++++++++++++ 3 files changed, 82 insertions(+), 31 deletions(-) delete mode 100644 test/regression/test_regression_11.py create mode 100644 test/test_regression/test_regression_15.py (limited to 'pyspike/generic.py') diff --git a/pyspike/generic.py b/pyspike/generic.py index 515cbf4..904c3c2 100644 --- a/pyspike/generic.py +++ b/pyspike/generic.py @@ -137,12 +137,13 @@ def _generic_distance_matrix(spike_trains, dist_function, 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:]] + pairs = [(i, j) for i in xrange(len(indices)) + for j in xrange(i+1, len(indices))] distance_matrix = np.zeros((len(indices), len(indices))) for i, j in pairs: - d = dist_function(spike_trains[i], spike_trains[j], interval) + d = dist_function(spike_trains[indices[i]], spike_trains[indices[j]], + interval) distance_matrix[i, j] = d distance_matrix[j, i] = d return distance_matrix diff --git a/test/regression/test_regression_11.py b/test/regression/test_regression_11.py deleted file mode 100644 index 7e66bc2..0000000 --- a/test/regression/test_regression_11.py +++ /dev/null @@ -1,28 +0,0 @@ -""" test_regression_11.py - -Regression test for Issue 11 -https://github.com/mariomulansky/PySpike/issues/11 - -Copyright 2015, Mario Mulansky - -Distributed under the BSD License -""" - -# run as -# python -Qnew test_regression_11.py -# to check correct behavior with new division - -# import division to see if everythin works with new division operator -from __future__ import division - -import pyspike as spk - -M = 19 # uneven number of spike trains - -spike_trains = [spk.generate_poisson_spikes(1.0, [0, 100]) for i in xrange(M)] - -isi_prof = spk.isi_profile_multi(spike_trains) -isi = spk.isi_distance_multi(spike_trains) - -spike_prof = spk.spike_profile_multi(spike_trains) -spike = spk.spike_distance_multi(spike_trains) diff --git a/test/test_regression/test_regression_15.py b/test/test_regression/test_regression_15.py new file mode 100644 index 0000000..1ce1290 --- /dev/null +++ b/test/test_regression/test_regression_15.py @@ -0,0 +1,78 @@ +""" test_regression_15.py + +Regression test for Issue #15 + +Copyright 2015, Mario Mulansky + +Distributed under the BSD License + +""" + +import numpy as np +from numpy.testing import assert_equal, assert_almost_equal, \ + assert_array_almost_equal + +import pyspike as spk + + +def test_regression_15_isi(): + # load spike trains + spike_trains = spk.load_spike_trains_from_txt("test/SPIKE_Sync_Test.txt", + edges=[0, 4000]) + + N = len(spike_trains) + + dist_mat = spk.isi_distance_matrix(spike_trains) + assert_equal(dist_mat.shape, (N, N)) + + ind = np.arange(N/2) + dist_mat = spk.isi_distance_matrix(spike_trains, ind) + assert_equal(dist_mat.shape, (N/2, N/2)) + + ind = np.arange(N/2, N) + dist_mat = spk.isi_distance_matrix(spike_trains, ind) + assert_equal(dist_mat.shape, (N/2, N/2)) + + +def test_regression_15_spike(): + # load spike trains + spike_trains = spk.load_spike_trains_from_txt("test/SPIKE_Sync_Test.txt", + edges=[0, 4000]) + + N = len(spike_trains) + + dist_mat = spk.spike_distance_matrix(spike_trains) + assert_equal(dist_mat.shape, (N, N)) + + ind = np.arange(N/2) + dist_mat = spk.spike_distance_matrix(spike_trains, ind) + assert_equal(dist_mat.shape, (N/2, N/2)) + + ind = np.arange(N/2, N) + dist_mat = spk.spike_distance_matrix(spike_trains, ind) + assert_equal(dist_mat.shape, (N/2, N/2)) + + +def test_regression_15_sync(): + # load spike trains + spike_trains = spk.load_spike_trains_from_txt("test/SPIKE_Sync_Test.txt", + edges=[0, 4000]) + + N = len(spike_trains) + + dist_mat = spk.spike_sync_matrix(spike_trains) + assert_equal(dist_mat.shape, (N, N)) + + ind = np.arange(N/2) + dist_mat = spk.spike_sync_matrix(spike_trains, ind) + assert_equal(dist_mat.shape, (N/2, N/2)) + + ind = np.arange(N/2, N) + dist_mat = spk.spike_sync_matrix(spike_trains, ind) + assert_equal(dist_mat.shape, (N/2, N/2)) + + +if __name__ == "__main__": + test_regression_15_isi() + test_regression_15_spike() + test_regression_15_sync() -- cgit v1.2.3 From 2c42e59e5097d3b9745e6eae2bee8f1ff27f7e09 Mon Sep 17 00:00:00 2001 From: Igor Gnatenko Date: Sun, 13 Dec 2015 10:57:12 +0100 Subject: py3: xrange() -> range() Signed-off-by: Igor Gnatenko --- examples/multivariate.py | 2 +- examples/performance.py | 2 +- pyspike/DiscreteFunc.py | 4 ++-- pyspike/directionality/spike_train_order.py | 4 ++-- pyspike/generic.py | 4 ++-- pyspike/psth.py | 2 +- test/test_distance.py | 6 +++--- 7 files changed, 12 insertions(+), 12 deletions(-) (limited to 'pyspike/generic.py') diff --git a/examples/multivariate.py b/examples/multivariate.py index 9a44758..93f8516 100644 --- a/examples/multivariate.py +++ b/examples/multivariate.py @@ -24,7 +24,7 @@ t_loading = time.clock() print("Number of spike trains: %d" % len(spike_trains)) num_of_spikes = sum([len(spike_trains[i]) - for i in xrange(len(spike_trains))]) + for i in range(len(spike_trains))]) print("Number of spikes: %d" % num_of_spikes) # calculate the multivariate spike distance diff --git a/examples/performance.py b/examples/performance.py index d0c3b91..ec6c830 100644 --- a/examples/performance.py +++ b/examples/performance.py @@ -26,7 +26,7 @@ print("%d spike trains with %d spikes" % (M, int(r*T))) spike_trains = [] t_start = datetime.now() -for i in xrange(M): +for i in range(M): spike_trains.append(spk.generate_poisson_spikes(r, T)) t_end = datetime.now() runtime = (t_end-t_start).total_seconds() diff --git a/pyspike/DiscreteFunc.py b/pyspike/DiscreteFunc.py index 55c0bc8..fe97bc2 100644 --- a/pyspike/DiscreteFunc.py +++ b/pyspike/DiscreteFunc.py @@ -80,7 +80,7 @@ class DiscreteFunc(object): expected_mp = (averaging_window_size+1) * int(self.mp[0]) y_plot = np.zeros_like(self.y) # compute the values in a loop, could be done in cython if required - for i in xrange(len(y_plot)): + for i in range(len(y_plot)): if self.mp[i] >= expected_mp: # the current value contains already all the wanted @@ -244,7 +244,7 @@ def average_profile(profiles): assert len(profiles) > 1 avrg_profile = profiles[0].copy() - for i in xrange(1, len(profiles)): + for i in range(1, len(profiles)): avrg_profile.add(profiles[i]) avrg_profile.mul_scalar(1.0/len(profiles)) # normalize diff --git a/pyspike/directionality/spike_train_order.py b/pyspike/directionality/spike_train_order.py index 44d931d..e6c9830 100644 --- a/pyspike/directionality/spike_train_order.py +++ b/pyspike/directionality/spike_train_order.py @@ -260,7 +260,7 @@ def optimal_spike_train_order(spike_trains, indices=None, interval=None, def permutate_matrix(D, p): N = len(D) D_p = np.empty_like(D) - for n in xrange(N): - for m in xrange(N): + for n in range(N): + for m in range(N): D_p[n, m] = D[p[n], p[m]] return D_p diff --git a/pyspike/generic.py b/pyspike/generic.py index 904c3c2..81ae660 100644 --- a/pyspike/generic.py +++ b/pyspike/generic.py @@ -137,8 +137,8 @@ def _generic_distance_matrix(spike_trains, dist_function, assert (indices < len(spike_trains)).all() and (indices >= 0).all(), \ "Invalid index list." # generate a list of possible index pairs - pairs = [(i, j) for i in xrange(len(indices)) - for j in xrange(i+1, len(indices))] + pairs = [(i, j) for i in range(len(indices)) + for j in range(i+1, len(indices))] distance_matrix = np.zeros((len(indices), len(indices))) for i, j in pairs: diff --git a/pyspike/psth.py b/pyspike/psth.py index 4027215..7cf1140 100644 --- a/pyspike/psth.py +++ b/pyspike/psth.py @@ -24,7 +24,7 @@ def psth(spike_trains, bin_size): # N = len(spike_trains) combined_spike_train = spike_trains[0].spikes - for i in xrange(1, len(spike_trains)): + for i in range(1, len(spike_trains)): combined_spike_train = np.append(combined_spike_train, spike_trains[i].spikes) diff --git a/test/test_distance.py b/test/test_distance.py index e45ac16..d5bce30 100644 --- a/test/test_distance.py +++ b/test/test_distance.py @@ -309,10 +309,10 @@ def check_dist_matrix(dist_func, dist_matrix_func): f_matrix = dist_matrix_func(spike_trains) # check zero diagonal - for i in xrange(4): + for i in range(4): assert_equal(0.0, f_matrix[i, i]) - for i in xrange(4): - for j in xrange(i+1, 4): + for i in range(4): + for j in range(i+1, 4): assert_equal(f_matrix[i, j], f_matrix[j, i]) assert_equal(f12, f_matrix[1, 0]) assert_equal(f13, f_matrix[2, 0]) -- cgit v1.2.3 From fe1f6179cce645df2511bfedae3af90167308f5f Mon Sep 17 00:00:00 2001 From: Igor Gnatenko Date: Sun, 13 Dec 2015 11:01:34 +0100 Subject: py3: division Signed-off-by: Igor Gnatenko --- pyspike/generic.py | 5 +++-- test/test_regression/test_regression_15.py | 26 ++++++++++++++------------ 2 files changed, 17 insertions(+), 14 deletions(-) (limited to 'pyspike/generic.py') diff --git a/pyspike/generic.py b/pyspike/generic.py index 81ae660..5ad06f1 100644 --- a/pyspike/generic.py +++ b/pyspike/generic.py @@ -7,6 +7,7 @@ Copyright 2015, Mario Mulansky Distributed under the BSD License """ +from __future__ import division import numpy as np @@ -38,14 +39,14 @@ def _generic_profile_multi(spike_trains, pair_distance_func, indices=None): L1 = len(pairs1) if L1 > 1: dist_prof1 = divide_and_conquer(pairs1[:L1//2], - pairs1[int(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[int(L2//2):]) + pairs2[L2//2:]) else: dist_prof2 = pair_distance_func(spike_trains[pairs2[0][0]], spike_trains[pairs2[0][1]]) diff --git a/test/test_regression/test_regression_15.py b/test/test_regression/test_regression_15.py index 1ce1290..42a39ea 100644 --- a/test/test_regression/test_regression_15.py +++ b/test/test_regression/test_regression_15.py @@ -8,6 +8,8 @@ Distributed under the BSD License """ +from __future__ import division + import numpy as np from numpy.testing import assert_equal, assert_almost_equal, \ assert_array_almost_equal @@ -25,13 +27,13 @@ def test_regression_15_isi(): dist_mat = spk.isi_distance_matrix(spike_trains) assert_equal(dist_mat.shape, (N, N)) - ind = np.arange(N/2) + ind = np.arange(N//2) dist_mat = spk.isi_distance_matrix(spike_trains, ind) - assert_equal(dist_mat.shape, (N/2, N/2)) + assert_equal(dist_mat.shape, (N//2, N//2)) - ind = np.arange(N/2, N) + ind = np.arange(N//2, N) dist_mat = spk.isi_distance_matrix(spike_trains, ind) - assert_equal(dist_mat.shape, (N/2, N/2)) + assert_equal(dist_mat.shape, (N//2, N//2)) def test_regression_15_spike(): @@ -44,13 +46,13 @@ def test_regression_15_spike(): dist_mat = spk.spike_distance_matrix(spike_trains) assert_equal(dist_mat.shape, (N, N)) - ind = np.arange(N/2) + ind = np.arange(N//2) dist_mat = spk.spike_distance_matrix(spike_trains, ind) - assert_equal(dist_mat.shape, (N/2, N/2)) + assert_equal(dist_mat.shape, (N//2, N//2)) - ind = np.arange(N/2, N) + ind = np.arange(N//2, N) dist_mat = spk.spike_distance_matrix(spike_trains, ind) - assert_equal(dist_mat.shape, (N/2, N/2)) + assert_equal(dist_mat.shape, (N//2, N//2)) def test_regression_15_sync(): @@ -63,13 +65,13 @@ def test_regression_15_sync(): dist_mat = spk.spike_sync_matrix(spike_trains) assert_equal(dist_mat.shape, (N, N)) - ind = np.arange(N/2) + ind = np.arange(N//2) dist_mat = spk.spike_sync_matrix(spike_trains, ind) - assert_equal(dist_mat.shape, (N/2, N/2)) + assert_equal(dist_mat.shape, (N//2, N//2)) - ind = np.arange(N/2, N) + ind = np.arange(N//2, N) dist_mat = spk.spike_sync_matrix(spike_trains, ind) - assert_equal(dist_mat.shape, (N/2, N/2)) + assert_equal(dist_mat.shape, (N//2, N//2)) if __name__ == "__main__": -- cgit v1.2.3