From d6462d271aeaf1be635cbc7c4317ae6a3b30b63f Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Fri, 12 Jun 2015 14:55:07 +0200 Subject: implement __getitem__ and __len__ for SpikeTrain This allows to use SpikeTrain objects to be used in many applications as if they were arrays with spike times. --- pyspike/SpikeTrain.py | 15 +++++++++++++++ pyspike/spikes.py | 3 ++- 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/pyspike/SpikeTrain.py b/pyspike/SpikeTrain.py index 9127b60..4b59a5d 100644 --- a/pyspike/SpikeTrain.py +++ b/pyspike/SpikeTrain.py @@ -32,6 +32,21 @@ class SpikeTrain(object): self.t_start = 0.0 self.t_end = float(edges) + def __getitem__(self, index): + """ Returns the time of the spike given by index. + + :param index: Index of the spike. + :return: spike time. + """ + return self.spikes[index] + + def __len__(self): + """ Returns the number of spikes. + + :return: Number of spikes. + """ + return len(self.spikes) + def sort(self): """ Sorts the spike times of this spike train using `np.sort` """ diff --git a/pyspike/spikes.py b/pyspike/spikes.py index 35d8533..b18d7eb 100644 --- a/pyspike/spikes.py +++ b/pyspike/spikes.py @@ -28,7 +28,8 @@ def spike_train_from_string(s, edges, sep=' ', is_sorted=False): # load_spike_trains_txt ############################################################ def load_spike_trains_from_txt(file_name, edges, - separator=' ', comment='#', is_sorted=False): + separator=' ', comment='#', is_sorted=False, + ignore_empty_lines=True): """ 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 -- cgit v1.2.3 From 7984d32e767e5833f1aaee06b6aeda8cc3f4500d Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Fri, 12 Jun 2015 14:57:50 +0200 Subject: update example to use new SpikeTrain capability Make use of __getitem__ and __len__ of SpikeTrains in some examples. --- examples/multivariate.py | 2 +- examples/plot.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/multivariate.py b/examples/multivariate.py index 53dbf0f..9a44758 100644 --- a/examples/multivariate.py +++ b/examples/multivariate.py @@ -23,7 +23,7 @@ spike_trains = spk.load_spike_trains_from_txt("PySpike_testdata.txt", t_loading = time.clock() print("Number of spike trains: %d" % len(spike_trains)) -num_of_spikes = sum([len(spike_trains[i].spikes) +num_of_spikes = sum([len(spike_trains[i]) for i in xrange(len(spike_trains))]) print("Number of spikes: %d" % num_of_spikes) diff --git a/examples/plot.py b/examples/plot.py index 9670286..5841baf 100644 --- a/examples/plot.py +++ b/examples/plot.py @@ -19,9 +19,9 @@ import pyspike as spk spike_trains = spk.load_spike_trains_from_txt("PySpike_testdata.txt", edges=(0, 4000)) -# plot the spike time +# plot the spike times for (i, spike_train) in enumerate(spike_trains): - plt.plot(spike_train.spikes, i*np.ones_like(spike_train.spikes), 'o') + plt.plot(spike_train, i*np.ones_like(spike_train), 'o') f = spk.isi_profile(spike_trains[0], spike_trains[1]) x, y = f.get_plottable_data() -- cgit v1.2.3 From 65883934996239907856b6e51641560dc7e1b7c7 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Wed, 24 Jun 2015 11:09:22 +0200 Subject: fixed typo --- pyspike/spike_distance.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyspike/spike_distance.py b/pyspike/spike_distance.py index ac2d260..dd6d4f8 100644 --- a/pyspike/spike_distance.py +++ b/pyspike/spike_distance.py @@ -54,7 +54,8 @@ Falling back to slow python backend.") ############################################################ def spike_distance(spike_train1, spike_train2, interval=None): """ Computes the spike-distance :math:`D_S` of the given spike trains. The - spike-distance is the integral over the isi distance profile :math:`S(t)`: + spike-distance is the integral over the spike distance profile + :math:`S(t)`: .. math:: D_S = \int_{T_0}^{T_1} S(t) dt. -- cgit v1.2.3 From 0ece782e1579660cdb71e077cbaaf9f76e97bef4 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Tue, 7 Jul 2015 18:06:12 +0200 Subject: better spike train plot (scatter) in plot.py --- examples/plot.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/plot.py b/examples/plot.py index 5841baf..c44afd1 100644 --- a/examples/plot.py +++ b/examples/plot.py @@ -21,7 +21,7 @@ spike_trains = spk.load_spike_trains_from_txt("PySpike_testdata.txt", # plot the spike times for (i, spike_train) in enumerate(spike_trains): - plt.plot(spike_train, i*np.ones_like(spike_train), 'o') + plt.scatter(spike_train, i*np.ones_like(spike_train), marker='|') f = spk.isi_profile(spike_trains[0], spike_trains[1]) x, y = f.get_plottable_data() -- cgit v1.2.3 From 5119d47d0f00c3f7203cf94460730b59a7e473ec Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Tue, 7 Jul 2015 18:55:32 +0200 Subject: add disable_backend_warning property Users can now disable the warning messages produced when the cython backend is not available by writing spk.disable_backend_warning = True in the beginning --- examples/performance.py | 3 +++ examples/plot.py | 1 + pyspike/DiscreteFunc.py | 4 +++- pyspike/PieceWiseConstFunc.py | 7 +++++-- pyspike/PieceWiseLinFunc.py | 9 ++++++--- pyspike/__init__.py | 2 ++ pyspike/isi_distance.py | 6 ++++-- pyspike/spike_distance.py | 4 +++- pyspike/spike_sync.py | 4 +++- 9 files changed, 30 insertions(+), 10 deletions(-) diff --git a/examples/performance.py b/examples/performance.py index 1c31e8f..d0c3b91 100644 --- a/examples/performance.py +++ b/examples/performance.py @@ -14,6 +14,9 @@ from datetime import datetime import cProfile import pstats +# in case you dont have the cython backends, disable the warnings as follows: +# spk.disable_backend_warning = True + M = 100 # number of spike trains r = 1.0 # rate of Poisson spike times T = 1E3 # length of spike trains diff --git a/examples/plot.py b/examples/plot.py index c44afd1..1922939 100644 --- a/examples/plot.py +++ b/examples/plot.py @@ -16,6 +16,7 @@ import matplotlib.pyplot as plt import pyspike as spk + spike_trains = spk.load_spike_trains_from_txt("PySpike_testdata.txt", edges=(0, 4000)) diff --git a/pyspike/DiscreteFunc.py b/pyspike/DiscreteFunc.py index 17153ee..a8c054e 100644 --- a/pyspike/DiscreteFunc.py +++ b/pyspike/DiscreteFunc.py @@ -6,6 +6,7 @@ from __future__ import print_function import numpy as np import collections +import pyspike ############################################################## @@ -202,7 +203,8 @@ class DiscreteFunc(object): from cython.cython_add import add_discrete_function_cython as \ add_discrete_function_impl except ImportError: - print("Warning: add_discrete_function_cython not found. Make \ + if not(pyspike.disable_backend_warning): + print("Warning: add_discrete_function_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.") diff --git a/pyspike/PieceWiseConstFunc.py b/pyspike/PieceWiseConstFunc.py index 2705443..23ff536 100644 --- a/pyspike/PieceWiseConstFunc.py +++ b/pyspike/PieceWiseConstFunc.py @@ -6,6 +6,7 @@ from __future__ import print_function import numpy as np import collections +import pyspike ############################################################## @@ -191,8 +192,10 @@ class PieceWiseConstFunc(object): from cython.cython_add import add_piece_wise_const_cython as \ add_piece_wise_const_impl except ImportError: - print("Warning: add_piece_wise_const_cython not found. Make sure \ -that PySpike is installed by running\n 'python setup.py build_ext --inplace'! \ + if not(pyspike.disable_backend_warning): + print("Warning: add_piece_wise_const_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 add_piece_wise_const_python as \ diff --git a/pyspike/PieceWiseLinFunc.py b/pyspike/PieceWiseLinFunc.py index c0dd475..0d51c76 100644 --- a/pyspike/PieceWiseLinFunc.py +++ b/pyspike/PieceWiseLinFunc.py @@ -6,6 +6,7 @@ from __future__ import print_function import numpy as np import collections +import pyspike ############################################################## @@ -230,9 +231,11 @@ class PieceWiseLinFunc: from cython.cython_add import add_piece_wise_lin_cython as \ add_piece_wise_lin_impl except ImportError: - print("Warning: add_piece_wise_lin_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.") + if not(pyspike.disable_backend_warning): + print("Warning: add_piece_wise_lin_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 add_piece_wise_lin_python as \ add_piece_wise_lin_impl diff --git a/pyspike/__init__.py b/pyspike/__init__.py index 3e836bd..2060f73 100644 --- a/pyspike/__init__.py +++ b/pyspike/__init__.py @@ -42,3 +42,5 @@ except DistributionNotFound: __version__ = 'Please install this project with setup.py' else: __version__ = _dist.version + +disable_backend_warning = False diff --git a/pyspike/isi_distance.py b/pyspike/isi_distance.py index 5ea555d..e50f203 100644 --- a/pyspike/isi_distance.py +++ b/pyspike/isi_distance.py @@ -2,6 +2,7 @@ # Copyright 2014-2015, Mario Mulansky # Distributed under the BSD License +import pyspike from pyspike import PieceWiseConstFunc from pyspike.generic import _generic_profile_multi, _generic_distance_multi, \ _generic_distance_matrix @@ -34,8 +35,9 @@ def isi_profile(spike_train1, spike_train2): from cython.cython_profiles import isi_profile_cython \ as isi_profile_impl except ImportError: - print("Warning: isi_distance_cython not found. Make sure that PySpike \ -is installed by running\n 'python setup.py build_ext --inplace'!\n \ + if not(pyspike.disable_backend_warning): + print("Warning: isi_profile_cython not found. Make sure that \ +PySpike is installed by running\n 'python setup.py build_ext --inplace'!\n \ Falling back to slow python backend.") # use python backend from cython.python_backend import isi_distance_python \ diff --git a/pyspike/spike_distance.py b/pyspike/spike_distance.py index dd6d4f8..feea0c1 100644 --- a/pyspike/spike_distance.py +++ b/pyspike/spike_distance.py @@ -2,6 +2,7 @@ # Copyright 2014-2015, Mario Mulansky # Distributed under the BSD License +import pyspike from pyspike import PieceWiseLinFunc from pyspike.generic import _generic_profile_multi, _generic_distance_multi, \ _generic_distance_matrix @@ -34,7 +35,8 @@ def spike_profile(spike_train1, spike_train2): from cython.cython_profiles import spike_profile_cython \ as spike_profile_impl except ImportError: - print("Warning: spike_profile_cython not found. Make sure that \ + if not(pyspike.disable_backend_warning): + print("Warning: spike_profile_cython not found. Make sure that \ PySpike is installed by running\n 'python setup.py build_ext --inplace'!\n \ Falling back to slow python backend.") # use python backend diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py index 40d98d2..10ebdc7 100644 --- a/pyspike/spike_sync.py +++ b/pyspike/spike_sync.py @@ -5,6 +5,7 @@ import numpy as np from functools import partial +import pyspike from pyspike import DiscreteFunc from pyspike.generic import _generic_profile_multi, _generic_distance_matrix @@ -39,7 +40,8 @@ def spike_sync_profile(spike_train1, spike_train2, max_tau=None): from cython.cython_profiles import coincidence_profile_cython \ as coincidence_profile_impl except ImportError: - print("Warning: spike_distance_cython not found. Make sure that \ + 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 -- cgit v1.2.3 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(-) 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(-) 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(-) 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(-) 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 2d9e1ea84162c8f964221df6851e4a4fe5e955e3 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Tue, 11 Aug 2015 15:38:08 +0200 Subject: add regression #11 test case --- test/regression/test_regression_11.py | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) create mode 100644 test/regression/test_regression_11.py diff --git a/test/regression/test_regression_11.py b/test/regression/test_regression_11.py new file mode 100644 index 0000000..7e66bc2 --- /dev/null +++ b/test/regression/test_regression_11.py @@ -0,0 +1,28 @@ +""" 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) -- 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 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 5731dad11d6d7129864fa6273d780000b34fd8a9 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Thu, 13 Aug 2015 17:45:02 +0200 Subject: more directionality functions + spike delay asymmetry profile + spike delay dual profile --- pyspike/directionality/__init__.py | 5 +- .../cython/cython_directionality.pyx | 46 ++++++++++++++++++ pyspike/directionality/spike_delay_asymmetry.py | 55 +++++++++++++++++++++- 3 files changed, 102 insertions(+), 4 deletions(-) diff --git a/pyspike/directionality/__init__.py b/pyspike/directionality/__init__.py index e6de1de..9a45b54 100644 --- a/pyspike/directionality/__init__.py +++ b/pyspike/directionality/__init__.py @@ -8,5 +8,6 @@ __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 + spike_delay_asymmetry_matrix, spike_delay_asymmetry_full, \ + optimal_asymmetry_order, optimal_asymmetry_order_from_D, \ + permutate_asymmetry_matrix diff --git a/pyspike/directionality/cython/cython_directionality.pyx b/pyspike/directionality/cython/cython_directionality.pyx index f5ea752..00e3b86 100644 --- a/pyspike/directionality/cython/cython_directionality.pyx +++ b/pyspike/directionality/cython/cython_directionality.pyx @@ -129,6 +129,52 @@ def spike_delay_asymmetry_profile_cython(double[:] spikes1, double[:] spikes2, +############################################################ +# spike_delay_dual_profile_cython +############################################################ +def spike_delay_dual_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 double[:] a1 = np.zeros(N1) # asymmetry values + cdef double[:] a2 = np.zeros(N2) # asymmetry values + 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 + 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 from spike train 1 after spike train 2 + # leading spike gets +1, following spike -1 + a1[i] = -1 + a2[j] = +1 + elif (j < N2-1) and (i == N1-1 or spikes1[i+1] > spikes2[j+1]): + j += 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 from spike train 1 before spike train 2 + # leading spike gets +1, following spike -1 + a1[i] = +1 + a2[j] = -1 + else: # spikes1[i+1] = spikes2[j+1] + # advance in both spike trains + j += 1 + i += 1 + # equal spike times: zero asymmetry value + a1[i] = 0 + a2[j] = 0 + + return a1, a2 + + ############################################################ # spike_delay_asymmetry_cython ############################################################ diff --git a/pyspike/directionality/spike_delay_asymmetry.py b/pyspike/directionality/spike_delay_asymmetry.py index 7d59601..0676f6d 100644 --- a/pyspike/directionality/spike_delay_asymmetry.py +++ b/pyspike/directionality/spike_delay_asymmetry.py @@ -146,6 +146,57 @@ def spike_delay_asymmetry_matrix(spike_trains, indices=None, return distance_matrix +############################################################ +# spike_delay_asymmetry_full +############################################################ +def spike_delay_asymmetry_full(spike_trains, indices=None, + interval=None, max_tau=None): + """ Computes the spike train symmetry value for each spike in each spike + train. + """ + 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." + # list of arrays for reulting asymmetry values + asymmetry_list = [np.zeros_like(st.spikes) for st in spike_trains] + # generate a list of possible index pairs + pairs = [(indices[i], j) for i in range(len(indices)) + for j in indices[i+1:]] + + # cython implementation + try: + from cython.cython_directionality import \ + spike_delay_dual_profile_cython as \ + sda_dual_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 + + for i, j in pairs: + a1, a2 = sda_dual_profile_impl(spike_trains[i].spikes, + spike_trains[j].spikes, + spike_trains[i].t_start, + spike_trains[i].t_end, + max_tau) + asymmetry_list[i] += a1 + asymmetry_list[j] += a2 + for a in asymmetry_list: + a /= len(spike_trains)-1 + return asymmetry_list + + ############################################################ # optimal_asymmetry_order_from_D ############################################################ @@ -188,7 +239,7 @@ def optimal_asymmetry_order_from_D(D, full_output=False): ############################################################ -# _optimal_asymmetry_order +# optimal_asymmetry_order ############################################################ def optimal_asymmetry_order(spike_trains, indices=None, interval=None, max_tau=None, full_output=False): @@ -203,7 +254,7 @@ def optimal_asymmetry_order(spike_trains, indices=None, interval=None, ############################################################ # reorder_asymmetry_matrix ############################################################ -def reorder_asymmetry_matrix(D, p): +def permutate_asymmetry_matrix(D, p): N = len(D) D_p = np.empty_like(D) for n in xrange(N): -- 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 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 0d7255c1fb7398720de10653efee617075c30892 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Mon, 17 Aug 2015 15:05:08 +0200 Subject: fix for #14 test case and fix for Issue #14. Spike-Sync function now correctly deal with empty intervals as well. --- pyspike/DiscreteFunc.py | 43 +++++++++++++++++++++++-------------------- test/test_empty.py | 12 ++++++++++++ 2 files changed, 35 insertions(+), 20 deletions(-) diff --git a/pyspike/DiscreteFunc.py b/pyspike/DiscreteFunc.py index a8c054e..4fd496d 100644 --- a/pyspike/DiscreteFunc.py +++ b/pyspike/DiscreteFunc.py @@ -137,9 +137,8 @@ class DiscreteFunc(object): :rtype: pair of float """ - if len(self.y) <= 2: - # no actual values in the profile, return spike sync of 1 - return 1.0, 1.0 + value = 0.0 + multiplicity = 0.0 def get_indices(ival): """ Retuns the indeces surrounding the given interval""" @@ -152,25 +151,29 @@ class DiscreteFunc(object): if interval is None: # no interval given, integrate over the whole spike train # don't count the first value, which is zero by definition - return (1.0 * np.sum(self.y[1:-1]), np.sum(self.mp[1:-1])) - - # check if interval is as sequence - assert isinstance(interval, collections.Sequence), \ - "Invalid value for `interval`. None, Sequence or Tuple expected." - # check if interval is a sequence of intervals - if not isinstance(interval[0], collections.Sequence): - # find the indices corresponding to the interval - start_ind, end_ind = get_indices(interval) - return (np.sum(self.y[start_ind:end_ind]), - np.sum(self.mp[start_ind:end_ind])) + value = 1.0 * np.sum(self.y[1:-1]) + multiplicity = np.sum(self.mp[1:-1]) else: - value = 0.0 - multiplicity = 0.0 - for ival in interval: + # check if interval is as sequence + assert isinstance(interval, collections.Sequence), \ + "Invalid value for `interval`. None, Sequence or Tuple \ +expected." + # check if interval is a sequence of intervals + if not isinstance(interval[0], collections.Sequence): # find the indices corresponding to the interval - 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]) + start_ind, end_ind = get_indices(interval) + value = np.sum(self.y[start_ind:end_ind]) + multiplicity = np.sum(self.mp[start_ind:end_ind]) + else: + for ival in interval: + # find the indices corresponding to the interval + 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): diff --git a/test/test_empty.py b/test/test_empty.py index 48be25d..af7fb36 100644 --- a/test/test_empty.py +++ b/test/test_empty.py @@ -139,6 +139,18 @@ def test_spike_sync_empty(): assert_array_almost_equal(prof.x, [0.0, 0.2, 0.8, 1.0], decimal=15) assert_array_almost_equal(prof.y, [0.0, 0.0, 0.0, 0.0], decimal=15) + # test with empty intervals + st1 = SpikeTrain([2.0, 5.0], [0, 10.0]) + st2 = SpikeTrain([2.1, 7.0], [0, 10.0]) + st3 = SpikeTrain([5.1, 6.0], [0, 10.0]) + res = spk.spike_sync_profile(st1, st2).avrg(interval=[3.0, 4.0]) + assert_equal(res, 1.0) + res = spk.spike_sync(st1, st2, interval=[3.0, 4.0]) + assert_equal(res, 1.0) + + sync_matrix = spk.spike_sync_matrix([st1, st2, st3], interval=[3.0, 4.0]) + assert_array_equal(sync_matrix, np.ones((3, 3)) - np.diag(np.ones(3))) + if __name__ == "__main__": test_get_non_empty() -- cgit v1.2.3 From 84f3333317e119cc87c4cb5f9c444ff66f1f7a23 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Mon, 17 Aug 2015 17:19:18 +0200 Subject: first unit tests for directionality --- .../cython/cython_directionality.pyx | 4 +-- pyspike/directionality/spike_delay_asymmetry.py | 11 +++--- test/test_spike_delay_asymmetry.py | 41 ++++++++++++++++++++++ 3 files changed, 50 insertions(+), 6 deletions(-) create mode 100644 test/test_spike_delay_asymmetry.py diff --git a/pyspike/directionality/cython/cython_directionality.pyx b/pyspike/directionality/cython/cython_directionality.pyx index 00e3b86..3936372 100644 --- a/pyspike/directionality/cython/cython_directionality.pyx +++ b/pyspike/directionality/cython/cython_directionality.pyx @@ -198,7 +198,7 @@ def spike_delay_asymmetry_cython(double[:] spikes1, double[:] spikes2, # 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 + asym -= 2 elif (j < N2-1) and (i == N1-1 or spikes1[i+1] > spikes2[j+1]): j += 1 mp += 1 @@ -207,7 +207,7 @@ def spike_delay_asymmetry_cython(double[:] spikes1, double[:] spikes2, # 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 + asym += 2 else: # spikes1[i+1] = spikes2[j+1] # advance in both spike trains j += 1 diff --git a/pyspike/directionality/spike_delay_asymmetry.py b/pyspike/directionality/spike_delay_asymmetry.py index 0676f6d..7da49ee 100644 --- a/pyspike/directionality/spike_delay_asymmetry.py +++ b/pyspike/directionality/spike_delay_asymmetry.py @@ -64,7 +64,7 @@ def spike_delay_asymmetry_profile(spike_train1, spike_train2, max_tau=None): ############################################################ # spike_delay_asymmetry ############################################################ -def spike_delay_asymmetry(spike_train1, spike_train2, +def spike_delay_asymmetry(spike_train1, spike_train2, normalize=True, interval=None, max_tau=None): """ Computes the overall spike delay asymmetry value for two spike trains. """ @@ -81,7 +81,10 @@ def spike_delay_asymmetry(spike_train1, spike_train2, spike_train1.t_start, spike_train1.t_end, max_tau) - return c + if normalize: + return 1.0*c/mp + else: + return c except ImportError: # Cython backend not available: fall back to profile averaging raise NotImplementedError() @@ -123,7 +126,7 @@ def spike_delay_asymmetry_profile_multi(spike_trains, indices=None, ############################################################ # spike_delay_asymmetry_matrix ############################################################ -def spike_delay_asymmetry_matrix(spike_trains, indices=None, +def spike_delay_asymmetry_matrix(spike_trains, normalize=True, indices=None, interval=None, max_tau=None): """ Computes the spike delay asymmetry matrix for the given spike trains. """ @@ -139,7 +142,7 @@ def spike_delay_asymmetry_matrix(spike_trains, indices=None, distance_matrix = np.zeros((len(indices), len(indices))) for i, j in pairs: - d = spike_delay_asymmetry(spike_trains[i], spike_trains[j], + d = spike_delay_asymmetry(spike_trains[i], spike_trains[j], normalize, interval, max_tau=max_tau) distance_matrix[i, j] = d distance_matrix[j, i] = -d diff --git a/test/test_spike_delay_asymmetry.py b/test/test_spike_delay_asymmetry.py new file mode 100644 index 0000000..93c18de --- /dev/null +++ b/test/test_spike_delay_asymmetry.py @@ -0,0 +1,41 @@ +""" test_spike_delay_asymmetry.py + +Tests the asymmetry functions + +Copyright 2015, Mario Mulansky + +Distributed under the BSD License + +""" + +import numpy as np +from numpy.testing import assert_equal, assert_almost_equal, \ + assert_array_equal + +import pyspike as spk +from pyspike import SpikeTrain, DiscreteFunc + + +def test_profile(): + st1 = SpikeTrain([100, 200, 300], [0, 1000]) + st2 = SpikeTrain([105, 205, 300], [0, 1000]) + expected_x = np.array([0, 100, 105, 200, 205, 300, 1000]) + expected_y = np.array([1, 1, 1, 1, 1, 0, 0]) + expected_mp = np.array([1, 1, 1, 1, 1, 2, 2]) + + f = spk.drct.spike_delay_asymmetry_profile(st1, st2) + + assert f.almost_equal(DiscreteFunc(expected_x, expected_y, expected_mp)) + assert_almost_equal(f.avrg(), 2.0/3.0) + assert_almost_equal(spk.drct.spike_delay_asymmetry(st1, st2), 2.0/3.0) + assert_almost_equal(spk.drct.spike_delay_asymmetry(st1, st2, + normalize=False), + 4.0) + + st3 = SpikeTrain([105, 195, 500], [0, 1000]) + expected_x = np.array([0, 100, 105, 195, 200, 300, 500, 1000]) + expected_y = np.array([1, 1, 1, -1, -1, 0, 0, 0]) + expected_mp = np.array([1, 1, 1, 1, 1, 1, 1, 1]) + + f = spk.drct.spike_delay_asymmetry_profile(st1, st3) + assert f.almost_equal(DiscreteFunc(expected_x, expected_y, expected_mp)) -- cgit v1.2.3 From 083a75679a2b15cceb13fa80ff1bb9277bd729d5 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Mon, 17 Aug 2015 17:48:06 +0200 Subject: major renaming of spike train order functions --- pyspike/DiscreteFunc.py | 7 +- pyspike/directionality/__init__.py | 12 +- .../cython/cython_directionality.pyx | 24 +- pyspike/directionality/spike_delay_asymmetry.py | 266 --------------------- pyspike/directionality/spike_train_order.py | 265 ++++++++++++++++++++ test/test_spike_delay_asymmetry.py | 9 +- 6 files changed, 292 insertions(+), 291 deletions(-) delete mode 100644 pyspike/directionality/spike_delay_asymmetry.py create mode 100644 pyspike/directionality/spike_train_order.py diff --git a/pyspike/DiscreteFunc.py b/pyspike/DiscreteFunc.py index 4fd496d..9cc7bd5 100644 --- a/pyspike/DiscreteFunc.py +++ b/pyspike/DiscreteFunc.py @@ -176,7 +176,7 @@ expected." multiplicity = 1.0 return (value, multiplicity) - def avrg(self, interval=None): + def avrg(self, interval=None, normalize=True): """ Computes the average of the interval sequence: :math:`a = 1/N \\sum f_n` where N is the number of intervals. @@ -189,7 +189,10 @@ expected." :rtype: float """ val, mp = self.integral(interval) - return val/mp + if normalize: + return val/mp + else: + return val def add(self, f): """ Adds another `DiscreteFunc` function to this function. diff --git a/pyspike/directionality/__init__.py b/pyspike/directionality/__init__.py index 9a45b54..6f74c50 100644 --- a/pyspike/directionality/__init__.py +++ b/pyspike/directionality/__init__.py @@ -4,10 +4,10 @@ Copyright 2015, Mario Mulansky Distributed under the BSD License """ -__all__ = ["spike_delay_asymmetry"] +__all__ = ["spike_train_order"] -from spike_delay_asymmetry import spike_delay_asymmetry_profile, \ - spike_delay_asymmetry, spike_delay_asymmetry_profile_multi, \ - spike_delay_asymmetry_matrix, spike_delay_asymmetry_full, \ - optimal_asymmetry_order, optimal_asymmetry_order_from_D, \ - permutate_asymmetry_matrix +from spike_train_order import spike_train_order_profile, \ + spike_train_order, spike_train_order_profile_multi, \ + spike_train_order_matrix, spike_order_values, \ + optimal_spike_train_order, optimal_spike_train_order_from_matrix, \ + permutate_matrix diff --git a/pyspike/directionality/cython/cython_directionality.pyx b/pyspike/directionality/cython/cython_directionality.pyx index 3936372..e1f63c4 100644 --- a/pyspike/directionality/cython/cython_directionality.pyx +++ b/pyspike/directionality/cython/cython_directionality.pyx @@ -61,11 +61,11 @@ cdef inline double get_tau(double[:] spikes1, double[:] spikes2, ############################################################ -# spike_delay_asymmetry_profile_cython +# spike_train_order_profile_cython ############################################################ -def spike_delay_asymmetry_profile_cython(double[:] spikes1, double[:] spikes2, - double t_start, double t_end, - double max_tau): +def spike_train_order_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) @@ -130,12 +130,12 @@ def spike_delay_asymmetry_profile_cython(double[:] spikes1, double[:] spikes2, ############################################################ -# spike_delay_dual_profile_cython +# spike_order_values_cython ############################################################ -def spike_delay_dual_profile_cython(double[:] spikes1, - double[:] spikes2, - double t_start, double t_end, - double max_tau): +def spike_order_values_cython(double[:] spikes1, + double[:] spikes2, + double t_start, double t_end, + double max_tau): cdef int N1 = len(spikes1) cdef int N2 = len(spikes2) @@ -176,10 +176,10 @@ def spike_delay_dual_profile_cython(double[:] spikes1, ############################################################ -# spike_delay_asymmetry_cython +# spike_train_order_cython ############################################################ -def spike_delay_asymmetry_cython(double[:] spikes1, double[:] spikes2, - double t_start, double t_end, double max_tau): +def spike_train_order_cython(double[:] spikes1, double[:] spikes2, + double t_start, double t_end, double max_tau): cdef int N1 = len(spikes1) cdef int N2 = len(spikes2) diff --git a/pyspike/directionality/spike_delay_asymmetry.py b/pyspike/directionality/spike_delay_asymmetry.py deleted file mode 100644 index 7da49ee..0000000 --- a/pyspike/directionality/spike_delay_asymmetry.py +++ /dev/null @@ -1,266 +0,0 @@ -# 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, normalize=True, - 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) - if normalize: - return 1.0*c/mp - else: - 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, normalize=True, 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], normalize, - interval, max_tau=max_tau) - distance_matrix[i, j] = d - distance_matrix[j, i] = -d - return distance_matrix - - -############################################################ -# spike_delay_asymmetry_full -############################################################ -def spike_delay_asymmetry_full(spike_trains, indices=None, - interval=None, max_tau=None): - """ Computes the spike train symmetry value for each spike in each spike - train. - """ - 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." - # list of arrays for reulting asymmetry values - asymmetry_list = [np.zeros_like(st.spikes) for st in spike_trains] - # generate a list of possible index pairs - pairs = [(indices[i], j) for i in range(len(indices)) - for j in indices[i+1:]] - - # cython implementation - try: - from cython.cython_directionality import \ - spike_delay_dual_profile_cython as \ - sda_dual_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 - - for i, j in pairs: - a1, a2 = sda_dual_profile_impl(spike_trains[i].spikes, - spike_trains[j].spikes, - spike_trains[i].t_start, - spike_trains[i].t_end, - max_tau) - asymmetry_list[i] += a1 - asymmetry_list[j] += a2 - for a in asymmetry_list: - a /= len(spike_trains)-1 - return asymmetry_list - - -############################################################ -# 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 permutate_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/directionality/spike_train_order.py b/pyspike/directionality/spike_train_order.py new file mode 100644 index 0000000..f8c8615 --- /dev/null +++ b/pyspike/directionality/spike_train_order.py @@ -0,0 +1,265 @@ +# Module containing functions to compute multivariate spike train order +# 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_train_order_profile +############################################################ +def spike_train_order_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_train_order_profile_cython as \ + spike_train_order_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_train_order_profile_impl(spike_train1.spikes, + spike_train2.spikes, + spike_train1.t_start, + spike_train1.t_end, + max_tau) + + return DiscreteFunc(times, coincidences, multiplicity) + + +############################################################ +# spike_train_order +############################################################ +def spike_train_order(spike_train1, spike_train2, normalize=True, + 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_train_order_cython as spike_train_order_impl + if max_tau is None: + max_tau = 0.0 + c, mp = spike_train_order_impl(spike_train1.spikes, + spike_train2.spikes, + spike_train1.t_start, + spike_train1.t_end, + max_tau) + if normalize: + return 1.0*c/mp + else: + 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_train_order_profile_multi +############################################################ +def spike_train_order_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_train_order_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_train_order_matrix +############################################################ +def spike_train_order_matrix(spike_trains, normalize=True, 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_train_order(spike_trains[i], spike_trains[j], normalize, + interval, max_tau=max_tau) + distance_matrix[i, j] = d + distance_matrix[j, i] = -d + return distance_matrix + + +############################################################ +# spike_order_values +############################################################ +def spike_order_values(spike_trains, indices=None, + interval=None, max_tau=None): + """ Computes the spike train symmetry value for each spike in each spike + train. + """ + 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." + # list of arrays for reulting asymmetry values + asymmetry_list = [np.zeros_like(st.spikes) for st in spike_trains] + # generate a list of possible index pairs + pairs = [(indices[i], j) for i in range(len(indices)) + for j in indices[i+1:]] + + # cython implementation + try: + from cython.cython_directionality import \ + spike_order_values_cython as spike_order_values_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 + + for i, j in pairs: + a1, a2 = spike_order_values_impl(spike_trains[i].spikes, + spike_trains[j].spikes, + spike_trains[i].t_start, + spike_trains[i].t_end, + max_tau) + asymmetry_list[i] += a1 + asymmetry_list[j] += a2 + for a in asymmetry_list: + a /= len(spike_trains)-1 + return asymmetry_list + + +############################################################ +# optimal_asymmetry_order_from_matrix +############################################################ +def optimal_spike_train_order_from_matrix(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_spike_train_order +############################################################ +def optimal_spike_train_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_train_order_matrix(spike_trains, indices, interval, max_tau) + return optimal_asymmetry_order_from_matrix(D, full_output) + + +############################################################ +# permutate_matrix +############################################################ +def permutate_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/test/test_spike_delay_asymmetry.py b/test/test_spike_delay_asymmetry.py index 93c18de..9de16e5 100644 --- a/test/test_spike_delay_asymmetry.py +++ b/test/test_spike_delay_asymmetry.py @@ -23,13 +23,12 @@ def test_profile(): expected_y = np.array([1, 1, 1, 1, 1, 0, 0]) expected_mp = np.array([1, 1, 1, 1, 1, 2, 2]) - f = spk.drct.spike_delay_asymmetry_profile(st1, st2) + f = spk.drct.spike_train_order_profile(st1, st2) assert f.almost_equal(DiscreteFunc(expected_x, expected_y, expected_mp)) assert_almost_equal(f.avrg(), 2.0/3.0) - assert_almost_equal(spk.drct.spike_delay_asymmetry(st1, st2), 2.0/3.0) - assert_almost_equal(spk.drct.spike_delay_asymmetry(st1, st2, - normalize=False), + assert_almost_equal(spk.drct.spike_train_order(st1, st2), 2.0/3.0) + assert_almost_equal(spk.drct.spike_train_order(st1, st2, normalize=False), 4.0) st3 = SpikeTrain([105, 195, 500], [0, 1000]) @@ -37,5 +36,5 @@ def test_profile(): expected_y = np.array([1, 1, 1, -1, -1, 0, 0, 0]) expected_mp = np.array([1, 1, 1, 1, 1, 1, 1, 1]) - f = spk.drct.spike_delay_asymmetry_profile(st1, st3) + f = spk.drct.spike_train_order_profile(st1, st3) assert f.almost_equal(DiscreteFunc(expected_x, expected_y, expected_mp)) -- cgit v1.2.3 From 691bf73c06322a2c47c37a5c48d085b789c8e8bf Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Mon, 17 Aug 2015 18:02:24 +0200 Subject: add python backend for directionality --- .../cython/directionality_python_backend.py | 89 ++++++++++++++++++++++ pyspike/directionality/spike_train_order.py | 33 ++++---- 2 files changed, 105 insertions(+), 17 deletions(-) create mode 100644 pyspike/directionality/cython/directionality_python_backend.py diff --git a/pyspike/directionality/cython/directionality_python_backend.py b/pyspike/directionality/cython/directionality_python_backend.py new file mode 100644 index 0000000..e14238f --- /dev/null +++ b/pyspike/directionality/cython/directionality_python_backend.py @@ -0,0 +1,89 @@ +""" directionality_python_backend.py + +Collection of python functions that can be used instead of the cython +implementation. + +Copyright 2015, Mario Mulansky + +Distributed under the BSD License + +""" + +import numpy as np + + +############################################################ +# spike_train_order_python +############################################################ +def spike_train_order_python(spikes1, spikes2, t_start, t_end, max_tau): + + def get_tau(spikes1, spikes2, i, j, max_tau): + m = t_end - t_start # use interval as initial tau + if i < len(spikes1)-1 and i > -1: + m = min(m, spikes1[i+1]-spikes1[i]) + if j < len(spikes2)-1 and j > -1: + m = min(m, spikes2[j+1]-spikes2[j]) + if i > 0: + m = min(m, spikes1[i]-spikes1[i-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 = -1 + j = -1 + n = 0 + st = np.zeros(N1 + N2 + 2) # spike times + a = 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 > -1 and spikes1[i]-spikes2[j] < tau: + # coincidence between the current spike and the previous spike + # 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, max_tau) + st[n] = spikes2[j] + if i > -1 and spikes2[j]-spikes1[i] < tau: + # coincidence between the current spike and the previous spike + # both get marked with 1 + 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 diff --git a/pyspike/directionality/spike_train_order.py b/pyspike/directionality/spike_train_order.py index f8c8615..892ffd0 100644 --- a/pyspike/directionality/spike_train_order.py +++ b/pyspike/directionality/spike_train_order.py @@ -5,7 +5,7 @@ import numpy as np from math import exp from functools import partial -# import pyspike +import pyspike from pyspike import DiscreteFunc from pyspike.generic import _generic_profile_multi @@ -39,14 +39,14 @@ def spike_train_order_profile(spike_train1, spike_train2, max_tau=None): spike_train_order_profile_cython as \ spike_train_order_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 + # 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.directionality_python_backend import \ + spike_train_order_python as spike_train_order_profile_impl if max_tau is None: max_tau = 0.0 @@ -81,15 +81,14 @@ def spike_train_order(spike_train1, spike_train2, normalize=True, spike_train1.t_start, spike_train1.t_end, max_tau) - if normalize: - return 1.0*c/mp - else: - 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) + c, mp = spike_train_order_profile(spike_train1, spike_train2, + max_tau).integral(interval) + if normalize: + return 1.0*c/mp + else: + return c else: # some specific interval is provided: not yet implemented raise NotImplementedError() @@ -250,7 +249,7 @@ def optimal_spike_train_order(spike_trains, indices=None, interval=None, Returns the optimal permutation p and A value. """ D = spike_train_order_matrix(spike_trains, indices, interval, max_tau) - return optimal_asymmetry_order_from_matrix(D, full_output) + return optimal_spike_train_order_from_matrix(D, full_output) ############################################################ -- cgit v1.2.3 From eeb4918ec2181f136e85bce976ec46a35a74b8f1 Mon Sep 17 00:00:00 2001 From: Igor Gnatenko Date: Sun, 13 Dec 2015 10:55:30 +0100 Subject: py3: absolute_import Signed-off-by: Igor Gnatenko --- pyspike/DiscreteFunc.py | 6 +++--- pyspike/PieceWiseConstFunc.py | 6 +++--- pyspike/PieceWiseLinFunc.py | 8 ++++---- pyspike/__init__.py | 22 ++++++++++++---------- pyspike/directionality/__init__.py | 4 +++- pyspike/directionality/spike_train_order.py | 12 +++++++----- pyspike/isi_distance.py | 8 +++++--- pyspike/spike_distance.py | 8 +++++--- pyspike/spike_sync.py | 8 +++++--- 9 files changed, 47 insertions(+), 35 deletions(-) diff --git a/pyspike/DiscreteFunc.py b/pyspike/DiscreteFunc.py index 9cc7bd5..55c0bc8 100644 --- a/pyspike/DiscreteFunc.py +++ b/pyspike/DiscreteFunc.py @@ -2,7 +2,7 @@ # Copyright 2014-2015, Mario Mulansky # Distributed under the BSD License -from __future__ import print_function +from __future__ import absolute_import, print_function import numpy as np import collections @@ -206,7 +206,7 @@ expected." # cython version try: - from cython.cython_add import add_discrete_function_cython as \ + from .cython.cython_add import add_discrete_function_cython as \ add_discrete_function_impl except ImportError: if not(pyspike.disable_backend_warning): @@ -215,7 +215,7 @@ 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 add_discrete_function_python as \ + from .cython.python_backend import add_discrete_function_python as \ add_discrete_function_impl self.x, self.y, self.mp = \ diff --git a/pyspike/PieceWiseConstFunc.py b/pyspike/PieceWiseConstFunc.py index 23ff536..5ce5f27 100644 --- a/pyspike/PieceWiseConstFunc.py +++ b/pyspike/PieceWiseConstFunc.py @@ -2,7 +2,7 @@ # Copyright 2014-2015, Mario Mulansky # Distributed under the BSD License -from __future__ import print_function +from __future__ import absolute_import, print_function import numpy as np import collections @@ -189,7 +189,7 @@ class PieceWiseConstFunc(object): # cython version try: - from cython.cython_add import add_piece_wise_const_cython as \ + from .cython.cython_add import add_piece_wise_const_cython as \ add_piece_wise_const_impl except ImportError: if not(pyspike.disable_backend_warning): @@ -198,7 +198,7 @@ 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 add_piece_wise_const_python as \ + from .cython.python_backend import add_piece_wise_const_python as \ add_piece_wise_const_impl self.x, self.y = add_piece_wise_const_impl(self.x, self.y, f.x, f.y) diff --git a/pyspike/PieceWiseLinFunc.py b/pyspike/PieceWiseLinFunc.py index 0d51c76..8145e63 100644 --- a/pyspike/PieceWiseLinFunc.py +++ b/pyspike/PieceWiseLinFunc.py @@ -2,7 +2,7 @@ # Copyright 2014-2015, Mario Mulansky # Distributed under the BSD License -from __future__ import print_function +from __future__ import absolute_import, print_function import numpy as np import collections @@ -222,13 +222,13 @@ class PieceWiseLinFunc: assert self.x[-1] == f.x[-1], "The functions have different intervals" # python implementation - # from python_backend import add_piece_wise_lin_python + # from .python_backend import add_piece_wise_lin_python # self.x, self.y1, self.y2 = add_piece_wise_lin_python( # self.x, self.y1, self.y2, f.x, f.y1, f.y2) # cython version try: - from cython.cython_add import add_piece_wise_lin_cython as \ + from .cython.cython_add import add_piece_wise_lin_cython as \ add_piece_wise_lin_impl except ImportError: if not(pyspike.disable_backend_warning): @@ -237,7 +237,7 @@ 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 add_piece_wise_lin_python as \ + from .cython.python_backend import add_piece_wise_lin_python as \ add_piece_wise_lin_impl self.x, self.y1, self.y2 = add_piece_wise_lin_impl( diff --git a/pyspike/__init__.py b/pyspike/__init__.py index 8d92ea4..335b1d3 100644 --- a/pyspike/__init__.py +++ b/pyspike/__init__.py @@ -4,27 +4,29 @@ Copyright 2014-2015, Mario Mulansky Distributed under the BSD License """ +from __future__ import absolute_import + __all__ = ["isi_distance", "spike_distance", "spike_sync", "psth", "spikes", "SpikeTrain", "PieceWiseConstFunc", "PieceWiseLinFunc", "DiscreteFunc", "directionality"] -from PieceWiseConstFunc import PieceWiseConstFunc -from PieceWiseLinFunc import PieceWiseLinFunc -from DiscreteFunc import DiscreteFunc -from SpikeTrain import SpikeTrain +from .PieceWiseConstFunc import PieceWiseConstFunc +from .PieceWiseLinFunc import PieceWiseLinFunc +from .DiscreteFunc import DiscreteFunc +from .SpikeTrain import SpikeTrain -from isi_distance import isi_profile, isi_distance, isi_profile_multi,\ +from .isi_distance import isi_profile, isi_distance, isi_profile_multi,\ isi_distance_multi, isi_distance_matrix -from spike_distance import spike_profile, spike_distance, spike_profile_multi,\ +from .spike_distance import spike_profile, spike_distance, spike_profile_multi,\ spike_distance_multi, spike_distance_matrix -from spike_sync import spike_sync_profile, spike_sync,\ +from .spike_sync import spike_sync_profile, spike_sync,\ spike_sync_profile_multi, spike_sync_multi, spike_sync_matrix -from psth import psth +from .psth import psth -from spikes import load_spike_trains_from_txt, spike_train_from_string, \ +from .spikes import load_spike_trains_from_txt, spike_train_from_string, \ merge_spike_trains, generate_poisson_spikes -import directionality as drct +from . 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 index 6f74c50..6ea38b2 100644 --- a/pyspike/directionality/__init__.py +++ b/pyspike/directionality/__init__.py @@ -4,9 +4,11 @@ Copyright 2015, Mario Mulansky Distributed under the BSD License """ +from __future__ import absolute_import + __all__ = ["spike_train_order"] -from spike_train_order import spike_train_order_profile, \ +from .spike_train_order import spike_train_order_profile, \ spike_train_order, spike_train_order_profile_multi, \ spike_train_order_matrix, spike_order_values, \ optimal_spike_train_order, optimal_spike_train_order_from_matrix, \ diff --git a/pyspike/directionality/spike_train_order.py b/pyspike/directionality/spike_train_order.py index 892ffd0..44d931d 100644 --- a/pyspike/directionality/spike_train_order.py +++ b/pyspike/directionality/spike_train_order.py @@ -2,6 +2,8 @@ # Copyright 2015, Mario Mulansky # Distributed under the BSD License +from __future__ import absolute_import + import numpy as np from math import exp from functools import partial @@ -35,7 +37,7 @@ def spike_train_order_profile(spike_train1, spike_train2, max_tau=None): # cython implementation try: - from cython.cython_directionality import \ + from .cython.cython_directionality import \ spike_train_order_profile_cython as \ spike_train_order_profile_impl except ImportError: @@ -45,7 +47,7 @@ def spike_train_order_profile(spike_train1, spike_train2, max_tau=None): PySpike is installed by running\n 'python setup.py build_ext --inplace'!\n \ Falling back to slow python backend.") # use python backend - from cython.directionality_python_backend import \ + from .cython.directionality_python_backend import \ spike_train_order_python as spike_train_order_profile_impl if max_tau is None: @@ -72,7 +74,7 @@ def spike_train_order(spike_train1, spike_train2, normalize=True, # distance over the whole interval is requested: use specific function # for optimal performance try: - from cython.cython_directionality import \ + from .cython.cython_directionality import \ spike_train_order_cython as spike_train_order_impl if max_tau is None: max_tau = 0.0 @@ -170,7 +172,7 @@ def spike_order_values(spike_trains, indices=None, # cython implementation try: - from cython.cython_directionality import \ + from .cython.cython_directionality import \ spike_order_values_cython as spike_order_values_impl except ImportError: raise NotImplementedError() @@ -179,7 +181,7 @@ def spike_order_values(spike_trains, indices=None, # 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 \ +# from .cython.python_backend import coincidence_python \ # as coincidence_profile_impl if max_tau is None: diff --git a/pyspike/isi_distance.py b/pyspike/isi_distance.py index e50f203..0ae7393 100644 --- a/pyspike/isi_distance.py +++ b/pyspike/isi_distance.py @@ -2,6 +2,8 @@ # Copyright 2014-2015, Mario Mulansky # Distributed under the BSD License +from __future__ import absolute_import + import pyspike from pyspike import PieceWiseConstFunc from pyspike.generic import _generic_profile_multi, _generic_distance_multi, \ @@ -32,7 +34,7 @@ def isi_profile(spike_train1, spike_train2): # load cython implementation try: - from cython.cython_profiles import isi_profile_cython \ + from .cython.cython_profiles import isi_profile_cython \ as isi_profile_impl except ImportError: if not(pyspike.disable_backend_warning): @@ -40,7 +42,7 @@ def isi_profile(spike_train1, spike_train2): PySpike is installed by running\n 'python setup.py build_ext --inplace'!\n \ Falling back to slow python backend.") # use python backend - from cython.python_backend import isi_distance_python \ + from .cython.python_backend import isi_distance_python \ as isi_profile_impl times, values = isi_profile_impl(spike_train1.get_spikes_non_empty(), @@ -74,7 +76,7 @@ def isi_distance(spike_train1, spike_train2, interval=None): # distance over the whole interval is requested: use specific function # for optimal performance try: - from cython.cython_distances import isi_distance_cython \ + from .cython.cython_distances import isi_distance_cython \ as isi_distance_impl return isi_distance_impl(spike_train1.get_spikes_non_empty(), diff --git a/pyspike/spike_distance.py b/pyspike/spike_distance.py index feea0c1..e418283 100644 --- a/pyspike/spike_distance.py +++ b/pyspike/spike_distance.py @@ -2,6 +2,8 @@ # Copyright 2014-2015, Mario Mulansky # Distributed under the BSD License +from __future__ import absolute_import + import pyspike from pyspike import PieceWiseLinFunc from pyspike.generic import _generic_profile_multi, _generic_distance_multi, \ @@ -32,7 +34,7 @@ def spike_profile(spike_train1, spike_train2): # cython implementation try: - from cython.cython_profiles import spike_profile_cython \ + from .cython.cython_profiles import spike_profile_cython \ as spike_profile_impl except ImportError: if not(pyspike.disable_backend_warning): @@ -40,7 +42,7 @@ def spike_profile(spike_train1, spike_train2): PySpike is installed by running\n 'python setup.py build_ext --inplace'!\n \ Falling back to slow python backend.") # use python backend - from cython.python_backend import spike_distance_python \ + from .cython.python_backend import spike_distance_python \ as spike_profile_impl times, y_starts, y_ends = spike_profile_impl( @@ -76,7 +78,7 @@ def spike_distance(spike_train1, spike_train2, interval=None): # distance over the whole interval is requested: use specific function # for optimal performance try: - from cython.cython_distances import spike_distance_cython \ + from .cython.cython_distances import spike_distance_cython \ as spike_distance_impl return spike_distance_impl(spike_train1.get_spikes_non_empty(), spike_train2.get_spikes_non_empty(), diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py index 10ebdc7..3dc29ff 100644 --- a/pyspike/spike_sync.py +++ b/pyspike/spike_sync.py @@ -3,6 +3,8 @@ # Copyright 2014-2015, Mario Mulansky # Distributed under the BSD License +from __future__ import absolute_import + import numpy as np from functools import partial import pyspike @@ -37,7 +39,7 @@ def spike_sync_profile(spike_train1, spike_train2, max_tau=None): # cython implementation try: - from cython.cython_profiles import coincidence_profile_cython \ + from .cython.cython_profiles import coincidence_profile_cython \ as coincidence_profile_impl except ImportError: if not(pyspike.disable_backend_warning): @@ -45,7 +47,7 @@ def spike_sync_profile(spike_train1, spike_train2, max_tau=None): 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 \ + from .cython.python_backend import coincidence_python \ as coincidence_profile_impl if max_tau is None: @@ -73,7 +75,7 @@ def _spike_sync_values(spike_train1, spike_train2, interval, max_tau): # distance over the whole interval is requested: use specific function # for optimal performance try: - from cython.cython_distances import coincidence_value_cython \ + from .cython.cython_distances import coincidence_value_cython \ as coincidence_value_impl if max_tau is None: max_tau = 0.0 -- 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(-) 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(-) 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 From 61a9785822d402d5269bdff608f6259c1b6bffa4 Mon Sep 17 00:00:00 2001 From: Igor Gnatenko Date: Sun, 13 Dec 2015 11:04:33 +0100 Subject: add metadata and testing under py3 Signed-off-by: Igor Gnatenko --- .travis.yml | 3 +++ setup.py | 5 +++++ 2 files changed, 8 insertions(+) diff --git a/.travis.yml b/.travis.yml index 1035775..1562d21 100644 --- a/.travis.yml +++ b/.travis.yml @@ -2,6 +2,9 @@ language: python python: - "2.6" - "2.7" + - "3.3" + - "3.4" + - "3.5" env: - CYTHON_INSTALL="pip install -q cython" - CYTHON_INSTALL="" diff --git a/setup.py b/setup.py index 960c684..0bd9173 100644 --- a/setup.py +++ b/setup.py @@ -89,6 +89,11 @@ train similarity', 'Programming Language :: Python :: 2', 'Programming Language :: Python :: 2.6', 'Programming Language :: Python :: 2.7', + + 'Programming Language :: Python :: 3', + 'Programming Language :: Python :: 3.3', + 'Programming Language :: Python :: 3.4', + 'Programming Language :: Python :: 3.5', ], package_data={ 'pyspike': ['cython/cython_add.c', 'cython/cython_profiles.c', -- cgit v1.2.3 From 692a5773c85e08a690bc5406147e58c4f10b6628 Mon Sep 17 00:00:00 2001 From: Igor Gnatenko Date: Sun, 13 Dec 2015 11:17:37 +0100 Subject: stop using package_data Signed-off-by: Igor Gnatenko --- MANIFEST.in | 7 +++++++ setup.py | 8 +------- 2 files changed, 8 insertions(+), 7 deletions(-) create mode 100644 MANIFEST.in diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000..aed0ae0 --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1,7 @@ +include *.rst +include *.txt +include pyspike/cython/*.c +include directionality/cython/*.c +recursive-include examples *.py *.txt +recursive-include test *.py *.txt +recursive-include doc * diff --git a/setup.py b/setup.py index 0bd9173..62c1951 100644 --- a/setup.py +++ b/setup.py @@ -94,11 +94,5 @@ train similarity', 'Programming Language :: Python :: 3.3', 'Programming Language :: Python :: 3.4', 'Programming Language :: Python :: 3.5', - ], - package_data={ - 'pyspike': ['cython/cython_add.c', 'cython/cython_profiles.c', - 'cython/cython_distances.c', - 'directionality/cython/cython_directionality.c'], - 'test': ['Spike_testdata.txt'] - } + ] ) -- cgit v1.2.3 From 45d4fac7b207bb5ce98f17433a443837fce0455a Mon Sep 17 00:00:00 2001 From: Igor Gnatenko Date: Mon, 14 Dec 2015 01:53:21 +0100 Subject: tests: allow to run out-of-tree Signed-off-by: Igor Gnatenko --- test/test_distance.py | 14 ++++++++------ test/test_regression/test_regression_15.py | 12 ++++++------ test/test_spikes.py | 9 +++++---- 3 files changed, 19 insertions(+), 16 deletions(-) diff --git a/test/test_distance.py b/test/test_distance.py index d5bce30..d405839 100644 --- a/test/test_distance.py +++ b/test/test_distance.py @@ -17,6 +17,8 @@ from numpy.testing import assert_equal, assert_almost_equal, \ import pyspike as spk from pyspike import SpikeTrain +import os +TEST_PATH = os.path.dirname(os.path.realpath(__file__)) def test_isi(): # generate two spike trains: @@ -275,8 +277,8 @@ def test_multi_spike_sync(): expected, decimal=15) # multivariate regression test - spike_trains = spk.load_spike_trains_from_txt("test/SPIKE_Sync_Test.txt", - edges=[0, 4000]) + spike_trains = spk.load_spike_trains_from_txt( + os.path.join(TEST_PATH, "SPIKE_Sync_Test.txt"), edges=[0, 4000]) # extract all spike times spike_times = np.array([]) for st in spike_trains: @@ -352,8 +354,8 @@ def test_regression_spiky(): # multivariate check - spike_trains = spk.load_spike_trains_from_txt("test/PySpike_testdata.txt", - (0.0, 4000.0)) + spike_trains = spk.load_spike_trains_from_txt( + os.path.join(TEST_PATH, "PySpike_testdata.txt"), (0.0, 4000.0)) isi_dist = spk.isi_distance_multi(spike_trains) # get the full precision from SPIKY assert_almost_equal(isi_dist, 0.17051816816999129656, decimal=15) @@ -371,8 +373,8 @@ def test_regression_spiky(): def test_multi_variate_subsets(): - spike_trains = spk.load_spike_trains_from_txt("test/PySpike_testdata.txt", - (0.0, 4000.0)) + spike_trains = spk.load_spike_trains_from_txt( + os.path.join(TEST_PATH, "PySpike_testdata.txt"), (0.0, 4000.0)) sub_set = [1, 3, 5, 7] spike_trains_sub_set = [spike_trains[i] for i in sub_set] diff --git a/test/test_regression/test_regression_15.py b/test/test_regression/test_regression_15.py index 42a39ea..dcacae2 100644 --- a/test/test_regression/test_regression_15.py +++ b/test/test_regression/test_regression_15.py @@ -16,11 +16,13 @@ from numpy.testing import assert_equal, assert_almost_equal, \ import pyspike as spk +import os +TEST_PATH = os.path.dirname(os.path.realpath(__file__)) +TEST_DATA = os.path.join(TEST_PATH, "..", "SPIKE_Sync_Test.txt") def test_regression_15_isi(): # load spike trains - spike_trains = spk.load_spike_trains_from_txt("test/SPIKE_Sync_Test.txt", - edges=[0, 4000]) + spike_trains = spk.load_spike_trains_from_txt(TEST_DATA, edges=[0, 4000]) N = len(spike_trains) @@ -38,8 +40,7 @@ def test_regression_15_isi(): def test_regression_15_spike(): # load spike trains - spike_trains = spk.load_spike_trains_from_txt("test/SPIKE_Sync_Test.txt", - edges=[0, 4000]) + spike_trains = spk.load_spike_trains_from_txt(TEST_DATA, edges=[0, 4000]) N = len(spike_trains) @@ -57,8 +58,7 @@ def test_regression_15_spike(): def test_regression_15_sync(): # load spike trains - spike_trains = spk.load_spike_trains_from_txt("test/SPIKE_Sync_Test.txt", - edges=[0, 4000]) + spike_trains = spk.load_spike_trains_from_txt(TEST_DATA, edges=[0, 4000]) N = len(spike_trains) diff --git a/test/test_spikes.py b/test/test_spikes.py index d4eb131..609a819 100644 --- a/test/test_spikes.py +++ b/test/test_spikes.py @@ -13,10 +13,12 @@ from numpy.testing import assert_equal import pyspike as spk +import os +TEST_PATH = os.path.dirname(os.path.realpath(__file__)) +TEST_DATA = os.path.join(TEST_PATH, "PySpike_testdata.txt") def test_load_from_txt(): - spike_trains = spk.load_spike_trains_from_txt("test/PySpike_testdata.txt", - edges=(0, 4000)) + spike_trains = spk.load_spike_trains_from_txt(TEST_DATA, edges=(0, 4000)) assert len(spike_trains) == 40 # check the first spike train @@ -48,8 +50,7 @@ def check_merged_spikes(merged_spikes, spike_trains): def test_merge_spike_trains(): # first load the data - spike_trains = spk.load_spike_trains_from_txt("test/PySpike_testdata.txt", - edges=(0, 4000)) + spike_trains = spk.load_spike_trains_from_txt(TEST_DATA, edges=(0, 4000)) merged_spikes = spk.merge_spike_trains([spike_trains[0], spike_trains[1]]) # test if result is sorted -- cgit v1.2.3 From 8baf994af3c5c8f7c76ca5a2c06e5220ba33604f Mon Sep 17 00:00:00 2001 From: Igor Gnatenko Date: Mon, 14 Dec 2015 01:54:11 +0100 Subject: contributors: add myself Signed-off-by: Igor Gnatenko --- Contributors.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/Contributors.txt b/Contributors.txt index 83563c7..512aa7f 100644 --- a/Contributors.txt +++ b/Contributors.txt @@ -1,5 +1,6 @@ Python/C Programming: - Mario Mulansky +- Igor Gnatenko Scientific Methods: - Thomas Kreuz -- cgit v1.2.3 From b970055641b215d30b671ee810e29c6a55e6214a Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Mon, 14 Dec 2015 14:23:02 +0100 Subject: improved edge correction for spike distance Improvement following Eeros suggestions to use auxiliary spike at the edges consistently with the corresponding corrected ISI intervals. --- pyspike/cython/cython_distances.pyx | 55 +++++++++++++++++++++------------ pyspike/cython/cython_profiles.pyx | 61 +++++++++++++++++++++++++------------ test/test_distance.py | 50 ++++++++++++++++++++++++++---- test/test_empty.py | 4 +-- 4 files changed, 122 insertions(+), 48 deletions(-) diff --git a/pyspike/cython/cython_distances.pyx b/pyspike/cython/cython_distances.pyx index c4f2349..41baa60 100644 --- a/pyspike/cython/cython_distances.pyx +++ b/pyspike/cython/cython_distances.pyx @@ -176,6 +176,8 @@ def spike_distance_cython(double[:] t1, double[:] t2, cdef double t_p1, t_f1, t_p2, t_f2, dt_p1, dt_p2, dt_f1, dt_f2 cdef double isi1, isi2, s1, s2 cdef double y_start, y_end, t_last, t_current, spike_value + cdef double[:] t_aux1 = np.empty(2) + cdef double[:] t_aux2 = np.empty(2) spike_value = 0.0 @@ -184,19 +186,27 @@ def spike_distance_cython(double[:] t1, double[:] t2, with nogil: # release the interpreter to allow multithreading t_last = t_start - t_p1 = t_start - t_p2 = t_start + # t_p1 = t_start + # t_p2 = t_start + # auxiliary spikes for edge correction - consistent with first/last ISI + t_aux1[0] = fmin(t_start, t1[0]-(t1[1]-t1[0])) + t_aux1[1] = fmax(t_end, t1[N1-1]+(t1[N1-1]-t1[N1-2])) + t_aux2[0] = fmin(t_start, t2[0]-(t2[1]-t2[0])) + t_aux2[1] = fmax(t_end, t2[N2-1]+(t2[N2-1]-t2[N2-2])) + t_p1 = t_start if (t1[0] == t_start) else t_aux1[0] + t_p2 = t_start if (t2[0] == t_start) else t_aux2[0] if t1[0] > t_start: # dt_p1 = t2[0]-t_start t_f1 = t1[0] - dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_start, t_end) + dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_aux2[0], t_aux2[1]) isi1 = fmax(t_f1-t_start, t1[1]-t1[0]) dt_p1 = dt_f1 - s1 = dt_p1*(t_f1-t_start)/isi1 + # s1 = dt_p1*(t_f1-t_start)/isi1 + s1 = dt_p1 index1 = -1 else: t_f1 = t1[1] - dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_start, t_end) + dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_aux2[0], t_aux2[1]) dt_p1 = 0.0 isi1 = t1[1]-t1[0] s1 = dt_p1 @@ -204,14 +214,15 @@ def spike_distance_cython(double[:] t1, double[:] t2, if t2[0] > t_start: # dt_p1 = t2[0]-t_start t_f2 = t2[0] - dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_start, t_end) + dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_aux1[0], t_aux1[1]) dt_p2 = dt_f2 isi2 = fmax(t_f2-t_start, t2[1]-t2[0]) - s2 = dt_p2*(t_f2-t_start)/isi2 + # s2 = dt_p2*(t_f2-t_start)/isi2 + s2 = dt_p2 index2 = -1 else: t_f2 = t2[1] - dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_start, t_end) + dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_aux1[0], t_aux1[1]) dt_p2 = 0.0 isi2 = t2[1]-t2[0] s2 = dt_p2 @@ -233,7 +244,7 @@ def spike_distance_cython(double[:] t1, double[:] t2, if index1 < N1-1: t_f1 = t1[index1+1] else: - t_f1 = t_end + t_f1 = t_aux1[1] t_curr = t_p1 s2 = (dt_p2*(t_f2-t_p1) + dt_f2*(t_p1-t_p2)) / isi2 y_end = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2) @@ -243,14 +254,16 @@ def spike_distance_cython(double[:] t1, double[:] t2, # now the next interval start value if index1 < N1-1: dt_f1 = get_min_dist_cython(t_f1, t2, N2, index2, - t_start, t_end) + t_aux2[0], t_aux2[1]) isi1 = t_f1-t_p1 s1 = dt_p1 else: dt_f1 = dt_p1 isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2]) # s1 needs adjustment due to change of isi1 - s1 = dt_p1*(t_end-t1[N1-1])/isi1 + # s1 = dt_p1*(t_end-t1[N1-1])/isi1 + # Eero's correction: no adjustment + s1 = dt_p1 # s2 is the same as above, thus we can compute y2 immediately y_start = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2) elif (index2 < N2-1) and (t_f1 > t_f2 or index1 == N1-1): @@ -264,7 +277,7 @@ def spike_distance_cython(double[:] t1, double[:] t2, if index2 < N2-1: t_f2 = t2[index2+1] else: - t_f2 = t_end + t_f2 = t_aux2[1] t_curr = t_p2 s1 = (dt_p1*(t_f1-t_p2) + dt_f1*(t_p2-t_p1)) / isi1 y_end = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2) @@ -274,14 +287,16 @@ def spike_distance_cython(double[:] t1, double[:] t2, # now the next interval start value if index2 < N2-1: dt_f2 = get_min_dist_cython(t_f2, t1, N1, index1, - t_start, t_end) + t_aux1[0], t_aux1[1]) isi2 = t_f2-t_p2 s2 = dt_p2 else: dt_f2 = dt_p2 isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2]) # s2 needs adjustment due to change of isi2 - s2 = dt_p2*(t_end-t2[N2-1])/isi2 + # s2 = dt_p2*(t_end-t2[N2-1])/isi2 + # Eero's correction: no adjustment + s2 = dt_p2 # s1 is the same as above, thus we can compute y2 immediately y_start = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2) else: # t_f1 == t_f2 - generate only one event @@ -298,27 +313,27 @@ def spike_distance_cython(double[:] t1, double[:] t2, if index1 < N1-1: t_f1 = t1[index1+1] dt_f1 = get_min_dist_cython(t_f1, t2, N2, index2, - t_start, t_end) + t_aux2[0], t_aux2[1]) isi1 = t_f1 - t_p1 else: - t_f1 = t_end + t_f1 = t_aux1[1] 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) + t_aux1[0], t_aux1[1]) isi2 = t_f2 - t_p2 else: - t_f2 = t_end + t_f2 = t_aux2[1] dt_f2 = dt_p2 isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2]) index += 1 t_last = t_curr # isi1 = max(t_end-t1[N1-1], t1[N1-1]-t1[N1-2]) # isi2 = max(t_end-t2[N2-1], t2[N2-1]-t2[N2-2]) - s1 = dt_f1*(t_end-t1[N1-1])/isi1 - s2 = dt_f2*(t_end-t2[N2-1])/isi2 + s1 = dt_f1 # *(t_end-t1[N1-1])/isi1 + s2 = dt_f2 # *(t_end-t2[N2-1])/isi2 y_end = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2) spike_value += 0.5*(y_start + y_end) * (t_end - t_last) # end nogil diff --git a/pyspike/cython/cython_profiles.pyx b/pyspike/cython/cython_profiles.pyx index f9893eb..f08de6e 100644 --- a/pyspike/cython/cython_profiles.pyx +++ b/pyspike/cython/cython_profiles.pyx @@ -181,6 +181,8 @@ def spike_profile_cython(double[:] t1, double[:] t2, cdef double[:] spike_events cdef double[:] y_starts cdef double[:] y_ends + cdef double[:] t_aux1 = np.empty(2) + cdef double[:] t_aux2 = np.empty(2) cdef int N1, N2, index1, index2, index cdef double t_p1, t_f1, t_p2, t_f2, dt_p1, dt_p2, dt_f1, dt_f2 @@ -189,6 +191,10 @@ def spike_profile_cython(double[:] t1, double[:] t2, N1 = len(t1) N2 = len(t2) + # we can assume at least two spikes per spike train + assert N1 > 1 + assert N2 > 1 + spike_events = np.empty(N1+N2+2) y_starts = np.empty(len(spike_events)-1) @@ -196,19 +202,27 @@ def spike_profile_cython(double[:] t1, double[:] t2, with nogil: # release the interpreter to allow multithreading spike_events[0] = t_start - t_p1 = t_start - t_p2 = t_start + # t_p1 = t_start + # t_p2 = t_start + # auxiliary spikes for edge correction - consistent with first/last ISI + t_aux1[0] = fmin(t_start, t1[0]-(t1[1]-t1[0])) + t_aux1[1] = fmax(t_end, t1[N1-1]+(t1[N1-1]-t1[N1-2])) + t_aux2[0] = fmin(t_start, t2[0]-(t2[1]-t2[0])) + t_aux2[1] = fmax(t_end, t2[N2-1]+(t2[N2-1]-t2[N2-2])) + t_p1 = t_start if (t1[0] == t_start) else t_aux1[0] + t_p2 = t_start if (t2[0] == t_start) else t_aux2[0] if t1[0] > t_start: # dt_p1 = t2[0]-t_start t_f1 = t1[0] - dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_start, t_end) + dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_aux2[0], t_aux2[1]) isi1 = fmax(t_f1-t_start, t1[1]-t1[0]) dt_p1 = dt_f1 - s1 = dt_p1*(t_f1-t_start)/isi1 + # s1 = dt_p1*(t_f1-t_start)/isi1 + s1 = dt_p1 index1 = -1 else: t_f1 = t1[1] - dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_start, t_end) + dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_aux2[0], t_aux2[1]) dt_p1 = 0.0 isi1 = t1[1]-t1[0] s1 = dt_p1 @@ -216,14 +230,15 @@ def spike_profile_cython(double[:] t1, double[:] t2, if t2[0] > t_start: # dt_p1 = t2[0]-t_start t_f2 = t2[0] - dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_start, t_end) + dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_aux1[0], t_aux1[1]) dt_p2 = dt_f2 isi2 = fmax(t_f2-t_start, t2[1]-t2[0]) - s2 = dt_p2*(t_f2-t_start)/isi2 + # s2 = dt_p2*(t_f2-t_start)/isi2 + s2 = dt_p2 index2 = -1 else: t_f2 = t2[1] - dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_start, t_end) + dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_aux1[0], t_aux1[1]) dt_p2 = 0.0 isi2 = t2[1]-t2[0] s2 = dt_p2 @@ -245,7 +260,7 @@ def spike_profile_cython(double[:] t1, double[:] t2, if index1 < N1-1: t_f1 = t1[index1+1] else: - t_f1 = t_end + t_f1 = t_aux1[1] spike_events[index] = t_p1 s2 = (dt_p2*(t_f2-t_p1) + dt_f2*(t_p1-t_p2)) / isi2 y_ends[index-1] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, @@ -253,14 +268,16 @@ def spike_profile_cython(double[:] t1, double[:] t2, # now the next interval start value if index1 < N1-1: dt_f1 = get_min_dist_cython(t_f1, t2, N2, index2, - t_start, t_end) + t_aux2[0], t_aux2[1]) isi1 = t_f1-t_p1 s1 = dt_p1 else: dt_f1 = dt_p1 isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2]) # s1 needs adjustment due to change of isi1 - s1 = dt_p1*(t_end-t1[N1-1])/isi1 + # s1 = dt_p1*(t_end-t1[N1-1])/isi1 + # Eero's correction: no adjustment + s1 = dt_p1 # s2 is the same as above, thus we can compute y2 immediately y_starts[index] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2) @@ -275,7 +292,7 @@ def spike_profile_cython(double[:] t1, double[:] t2, if index2 < N2-1: t_f2 = t2[index2+1] else: - t_f2 = t_end + t_f2 = t_aux2[1] spike_events[index] = t_p2 s1 = (dt_p1*(t_f1-t_p2) + dt_f1*(t_p2-t_p1)) / isi1 y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, @@ -283,14 +300,16 @@ def spike_profile_cython(double[:] t1, double[:] t2, # now the next interval start value if index2 < N2-1: dt_f2 = get_min_dist_cython(t_f2, t1, N1, index1, - t_start, t_end) + t_aux1[0], t_aux1[1]) isi2 = t_f2-t_p2 s2 = dt_p2 else: dt_f2 = dt_p2 isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2]) # s2 needs adjustment due to change of isi2 - s2 = dt_p2*(t_end-t2[N2-1])/isi2 + # s2 = dt_p2*(t_end-t2[N2-1])/isi2 + # Eero's correction: no adjustment + s2 = dt_p2 # s2 is the same as above, thus we can compute y2 immediately y_starts[index] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2) else: # t_f1 == t_f2 - generate only one event @@ -306,19 +325,19 @@ def spike_profile_cython(double[:] t1, double[:] t2, if index1 < N1-1: t_f1 = t1[index1+1] dt_f1 = get_min_dist_cython(t_f1, t2, N2, index2, - t_start, t_end) + t_aux2[0], t_aux2[1]) isi1 = t_f1 - t_p1 else: - t_f1 = t_end + t_f1 = t_aux1[1] 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) + t_aux1[0], t_aux1[1]) isi2 = t_f2 - t_p2 else: - t_f2 = t_end + t_f2 = t_aux2[1] dt_f2 = dt_p2 isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2]) index += 1 @@ -330,8 +349,10 @@ def spike_profile_cython(double[:] t1, double[:] t2, # 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 + # s1 = dt_f1*(t_end-t1[N1-1])/isi1 + # s2 = dt_f2*(t_end-t2[N2-1])/isi2 + s1 = dt_f1 + s2 = dt_f2 y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2) # end nogil diff --git a/test/test_distance.py b/test/test_distance.py index e45ac16..626b438 100644 --- a/test/test_distance.py +++ b/test/test_distance.py @@ -36,6 +36,7 @@ def test_isi(): f = spk.isi_profile(t1, t2) # print("ISI: ", f.y) + print("ISI value:", expected_isi_val) assert_equal(f.x, expected_times) assert_array_almost_equal(f.y, expected_isi, decimal=15) @@ -73,8 +74,19 @@ def test_spike(): assert_equal(f.x, expected_times) - assert_almost_equal(f.avrg(), 1.6624149659863946e-01, decimal=15) - assert_almost_equal(f.y2[-1], 0.1394558, decimal=6) + # from SPIKY: + y_all = np.array([0.000000000000000000, 0.555555555555555580, + 0.222222222222222210, 0.305555555555555580, + 0.255102040816326536, 0.000000000000000000, + 0.000000000000000000, 0.255102040816326536, + 0.255102040816326536, 0.285714285714285698, + 0.285714285714285698, 0.285714285714285698]) + + #assert_array_almost_equal(f.y1, y_all[::2]) + assert_array_almost_equal(f.y2, y_all[1::2]) + + assert_almost_equal(f.avrg(), 0.186309523809523814, decimal=15) + assert_equal(spk.spike_distance(t1, t2), f.avrg()) 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) @@ -99,6 +111,8 @@ def test_spike(): (expected_y1+expected_y2)/2) expected_spike_val /= (expected_times[-1]-expected_times[0]) + print("SPIKE value:", expected_spike_val) + f = spk.spike_profile(t1, t2) assert_equal(f.x, expected_times) @@ -117,9 +131,14 @@ def test_spike(): # 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, + # 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]) + # eero's edge correction: + s2_r = np.array([0.1, 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, + s2_l = np.array([0.1, 0.1*0.3/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]) @@ -345,7 +364,7 @@ def test_regression_spiky(): 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) + assert_equal(spike_dist, 0.211058782487353908) spike_sync = spk.spike_sync(st1, st2) assert_equal(spike_sync, 8.6956521739130432e-01) @@ -363,12 +382,31 @@ def test_regression_spiky(): spike_dist = spk.spike_distance_multi(spike_trains) # get the full precision from SPIKY - assert_almost_equal(spike_dist, 2.4432433330596512e-01, decimal=15) + assert_almost_equal(spike_dist, 0.25188056475463755, decimal=15) spike_sync = spk.spike_sync_multi(spike_trains) # get the full precision from SPIKY assert_equal(spike_sync, 0.7183531505298066) + # Eero's edge correction example + st1 = SpikeTrain([0.5, 1.5, 2.5], 6.0) + st2 = SpikeTrain([3.5, 4.5, 5.5], 6.0) + + f = spk.spike_profile(st1, st2) + + expected_times = np.array([0.0, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.0]) + y_all = np.array([0.271604938271605, 0.271604938271605, 0.271604938271605, + 0.617283950617284, 0.617283950617284, 0.444444444444444, + 0.285714285714286, 0.285714285714286, 0.444444444444444, + 0.617283950617284, 0.617283950617284, 0.271604938271605, + 0.271604938271605, 0.271604938271605]) + expected_y1 = y_all[::2] + expected_y2 = y_all[1::2] + + assert_equal(f.x, expected_times) + assert_array_almost_equal(f.y1, expected_y1, decimal=14) + assert_array_almost_equal(f.y2, expected_y2, decimal=14) + def test_multi_variate_subsets(): spike_trains = spk.load_spike_trains_from_txt("test/PySpike_testdata.txt", diff --git a/test/test_empty.py b/test/test_empty.py index af7fb36..5a0042f 100644 --- a/test/test_empty.py +++ b/test/test_empty.py @@ -70,8 +70,8 @@ def test_spike_empty(): st1 = SpikeTrain([], edges=(0.0, 1.0)) st2 = SpikeTrain([0.4, ], edges=(0.0, 1.0)) d = spk.spike_distance(st1, st2) - assert_almost_equal(d, 0.4*0.4*1.0/(0.4+1.0)**2 + 0.6*0.4*1.0/(0.6+1.0)**2, - decimal=15) + d_expect = 0.4*0.4*1.0/(0.4+1.0)**2 + 0.6*0.4*1.0/(0.6+1.0)**2 + assert_almost_equal(d, d_expect, decimal=15) prof = spk.spike_profile(st1, st2) assert_equal(d, prof.avrg()) assert_array_equal(prof.x, [0.0, 0.4, 1.0]) -- cgit v1.2.3 From a37dc02a0f6c11a8773cba8a6137a2604222fa26 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Mon, 14 Dec 2015 16:17:28 +0100 Subject: python backend updated for new edge correction Eero's improved edge correction now also implemented in the python backend. --- pyspike/cython/python_backend.py | 53 +++++++++++++++++++++++++--------------- 1 file changed, 33 insertions(+), 20 deletions(-) diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py index 69a420f..a5f7ae4 100644 --- a/pyspike/cython/python_backend.py +++ b/pyspike/cython/python_backend.py @@ -144,35 +144,44 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end): y_starts = np.empty(len(spike_events)-1) y_ends = np.empty(len(spike_events)-1) + t_aux1 = np.zeros(2) + t_aux2 = np.zeros(2) + t_aux1[0] = min(t_start, t1[0]-(t1[1]-t1[0])) + t_aux1[1] = max(t_end, t1[N1-1]+(t1[N1-1]-t1[N1-2])) + t_aux2[0] = min(t_start, t2[0]-(t2[1]-t2[0])) + t_aux2[1] = max(t_end, t2[N2-1]+(t2[N2-1]-t2[N2-2])) + t_p1 = t_start if (t1[0] == t_start) else t_aux1[0] + t_p2 = t_start if (t2[0] == t_start) else t_aux2[0] + spike_events[0] = t_start - t_p1 = t_start - t_p2 = t_start if t1[0] > t_start: t_f1 = t1[0] - dt_f1 = get_min_dist(t_f1, t2, 0, t_start, t_end) + dt_f1 = get_min_dist(t_f1, t2, 0, t_aux2[0], t_aux2[1]) dt_p1 = dt_f1 isi1 = max(t_f1-t_start, t1[1]-t1[0]) - s1 = dt_p1*(t_f1-t_start)/isi1 + # s1 = dt_p1*(t_f1-t_start)/isi1 + s1 = dt_p1 index1 = -1 else: dt_p1 = 0.0 t_f1 = t1[1] - dt_f1 = get_min_dist(t_f1, t2, 0, t_start, t_end) + dt_f1 = get_min_dist(t_f1, t2, 0, t_aux2[0], t_aux2[1]) isi1 = t1[1]-t1[0] s1 = dt_p1 index1 = 0 if t2[0] > t_start: # dt_p1 = t2[0]-t_start t_f2 = t2[0] - dt_f2 = get_min_dist(t_f2, t1, 0, t_start, t_end) + dt_f2 = get_min_dist(t_f2, t1, 0, t_aux1[0], t_aux1[1]) dt_p2 = dt_f2 isi2 = max(t_f2-t_start, t2[1]-t2[0]) - s2 = dt_p2*(t_f2-t_start)/isi2 + # s2 = dt_p2*(t_f2-t_start)/isi2 + s2 = dt_p2 index2 = -1 else: dt_p2 = 0.0 t_f2 = t2[1] - dt_f2 = get_min_dist(t_f2, t1, 0, t_start, t_end) + dt_f2 = get_min_dist(t_f2, t1, 0, t_aux1[0], t_aux1[1]) isi2 = t2[1]-t2[0] s2 = dt_p2 index2 = 0 @@ -193,20 +202,22 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end): if index1 < N1-1: t_f1 = t1[index1+1] else: - t_f1 = t_end + t_f1 = t_aux1[1] spike_events[index] = t_p1 s2 = (dt_p2*(t_f2-t_p1) + dt_f2*(t_p1-t_p2)) / isi2 y_ends[index-1] = (s1*isi2 + s2*isi1) / (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) + dt_f1 = get_min_dist(t_f1, t2, index2, t_aux2[0], t_aux2[1]) 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]) # s1 needs adjustment due to change of isi1 - s1 = dt_p1*(t_end-t1[N1-1])/isi1 + # s1 = dt_p1*(t_end-t1[N1-1])/isi1 + # Eero's correction: no adjustment + s1 = dt_p1 # s2 is the same as above, thus we can compute y2 immediately 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): @@ -220,20 +231,22 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end): if index2 < N2-1: t_f2 = t2[index2+1] else: - t_f2 = t_end + t_f2 = t_aux2[1] spike_events[index] = t_p2 s1 = (dt_p1*(t_f1-t_p2) + dt_f1*(t_p2-t_p1)) / isi1 y_ends[index-1] = (s1*isi2 + s2*isi1) / (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) + dt_f2 = get_min_dist(t_f2, t1, index1, t_aux1[0], t_aux1[1]) 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]) # s2 needs adjustment due to change of isi2 - s2 = dt_p2*(t_end-t2[N2-1])/isi2 + # s2 = dt_p2*(t_end-t2[N2-1])/isi2 + # Eero's adjustment: no correction + s2 = dt_p2 # s2 is the same as above, thus we can compute y2 immediately y_starts[index] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2) else: # t_f1 == t_f2 - generate only one event @@ -248,18 +261,18 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end): 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) + dt_f1 = get_min_dist(t_f1, t2, index2, t_aux2[0], t_aux2[1]) isi1 = t_f1 - t_p1 else: - t_f1 = t_end + t_f1 = t_aux1[1] 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) + dt_f2 = get_min_dist(t_f2, t1, index1, t_aux1[0], t_aux1[1]) isi2 = t_f2 - t_p2 else: - t_f2 = t_end + t_f2 = t_aux2[1] dt_f2 = dt_p2 isi2 = max(t_end-t2[N2-1], t2[N2-1]-t2[N2-2]) index += 1 @@ -271,8 +284,8 @@ def spike_distance_python(spikes1, spikes2, t_start, 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 + 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 -- cgit v1.2.3 From 2958fac1f12dfc743214615a3e8252dae54ea784 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Mon, 14 Dec 2015 16:38:48 +0100 Subject: remove directionality from develop directionality development happens in separate branch. --- pyspike/directionality/__init__.py | 15 -- pyspike/directionality/cython/__init__.py | 0 .../cython/cython_directionality.pyx | 223 ----------------- .../cython/directionality_python_backend.py | 89 ------- pyspike/directionality/spike_train_order.py | 266 --------------------- test/test_spike_delay_asymmetry.py | 40 ---- 6 files changed, 633 deletions(-) delete mode 100644 pyspike/directionality/__init__.py delete mode 100644 pyspike/directionality/cython/__init__.py delete mode 100644 pyspike/directionality/cython/cython_directionality.pyx delete mode 100644 pyspike/directionality/cython/directionality_python_backend.py delete mode 100644 pyspike/directionality/spike_train_order.py delete mode 100644 test/test_spike_delay_asymmetry.py diff --git a/pyspike/directionality/__init__.py b/pyspike/directionality/__init__.py deleted file mode 100644 index 6ea38b2..0000000 --- a/pyspike/directionality/__init__.py +++ /dev/null @@ -1,15 +0,0 @@ -""" -Copyright 2015, Mario Mulansky - -Distributed under the BSD License -""" - -from __future__ import absolute_import - -__all__ = ["spike_train_order"] - -from .spike_train_order import spike_train_order_profile, \ - spike_train_order, spike_train_order_profile_multi, \ - spike_train_order_matrix, spike_order_values, \ - optimal_spike_train_order, optimal_spike_train_order_from_matrix, \ - permutate_matrix diff --git a/pyspike/directionality/cython/__init__.py b/pyspike/directionality/cython/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/pyspike/directionality/cython/cython_directionality.pyx b/pyspike/directionality/cython/cython_directionality.pyx deleted file mode 100644 index e1f63c4..0000000 --- a/pyspike/directionality/cython/cython_directionality.pyx +++ /dev/null @@ -1,223 +0,0 @@ -#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_train_order_profile_cython -############################################################ -def spike_train_order_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_order_values_cython -############################################################ -def spike_order_values_cython(double[:] spikes1, - double[:] spikes2, - double t_start, double t_end, - double max_tau): - - cdef int N1 = len(spikes1) - cdef int N2 = len(spikes2) - cdef int i = -1 - cdef int j = -1 - cdef double[:] a1 = np.zeros(N1) # asymmetry values - cdef double[:] a2 = np.zeros(N2) # asymmetry values - 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 - 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 from spike train 1 after spike train 2 - # leading spike gets +1, following spike -1 - a1[i] = -1 - a2[j] = +1 - elif (j < N2-1) and (i == N1-1 or spikes1[i+1] > spikes2[j+1]): - j += 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 from spike train 1 before spike train 2 - # leading spike gets +1, following spike -1 - a1[i] = +1 - a2[j] = -1 - else: # spikes1[i+1] = spikes2[j+1] - # advance in both spike trains - j += 1 - i += 1 - # equal spike times: zero asymmetry value - a1[i] = 0 - a2[j] = 0 - - return a1, a2 - - -############################################################ -# spike_train_order_cython -############################################################ -def spike_train_order_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 -= 2 - elif (j < N2-1) and (i == N1-1 or spikes1[i+1] > spikes2[j+1]): - j += 1 - mp += 1 - tau = get_tau(spikes1, spikes2, i, j, 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 += 2 - 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/cython/directionality_python_backend.py b/pyspike/directionality/cython/directionality_python_backend.py deleted file mode 100644 index e14238f..0000000 --- a/pyspike/directionality/cython/directionality_python_backend.py +++ /dev/null @@ -1,89 +0,0 @@ -""" directionality_python_backend.py - -Collection of python functions that can be used instead of the cython -implementation. - -Copyright 2015, Mario Mulansky - -Distributed under the BSD License - -""" - -import numpy as np - - -############################################################ -# spike_train_order_python -############################################################ -def spike_train_order_python(spikes1, spikes2, t_start, t_end, max_tau): - - def get_tau(spikes1, spikes2, i, j, max_tau): - m = t_end - t_start # use interval as initial tau - if i < len(spikes1)-1 and i > -1: - m = min(m, spikes1[i+1]-spikes1[i]) - if j < len(spikes2)-1 and j > -1: - m = min(m, spikes2[j+1]-spikes2[j]) - if i > 0: - m = min(m, spikes1[i]-spikes1[i-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 = -1 - j = -1 - n = 0 - st = np.zeros(N1 + N2 + 2) # spike times - a = 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 > -1 and spikes1[i]-spikes2[j] < tau: - # coincidence between the current spike and the previous spike - # 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, max_tau) - st[n] = spikes2[j] - if i > -1 and spikes2[j]-spikes1[i] < tau: - # coincidence between the current spike and the previous spike - # both get marked with 1 - 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 diff --git a/pyspike/directionality/spike_train_order.py b/pyspike/directionality/spike_train_order.py deleted file mode 100644 index e6c9830..0000000 --- a/pyspike/directionality/spike_train_order.py +++ /dev/null @@ -1,266 +0,0 @@ -# Module containing functions to compute multivariate spike train order -# Copyright 2015, Mario Mulansky -# Distributed under the BSD License - -from __future__ import absolute_import - -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_train_order_profile -############################################################ -def spike_train_order_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_train_order_profile_cython as \ - spike_train_order_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.directionality_python_backend import \ - spike_train_order_python as spike_train_order_profile_impl - - if max_tau is None: - max_tau = 0.0 - - times, coincidences, multiplicity \ - = spike_train_order_profile_impl(spike_train1.spikes, - spike_train2.spikes, - spike_train1.t_start, - spike_train1.t_end, - max_tau) - - return DiscreteFunc(times, coincidences, multiplicity) - - -############################################################ -# spike_train_order -############################################################ -def spike_train_order(spike_train1, spike_train2, normalize=True, - 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_train_order_cython as spike_train_order_impl - if max_tau is None: - max_tau = 0.0 - c, mp = spike_train_order_impl(spike_train1.spikes, - spike_train2.spikes, - spike_train1.t_start, - spike_train1.t_end, - max_tau) - except ImportError: - # Cython backend not available: fall back to profile averaging - c, mp = spike_train_order_profile(spike_train1, spike_train2, - max_tau).integral(interval) - if normalize: - return 1.0*c/mp - else: - return c - else: - # some specific interval is provided: not yet implemented - raise NotImplementedError() - - -############################################################ -# spike_train_order_profile_multi -############################################################ -def spike_train_order_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_train_order_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_train_order_matrix -############################################################ -def spike_train_order_matrix(spike_trains, normalize=True, 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_train_order(spike_trains[i], spike_trains[j], normalize, - interval, max_tau=max_tau) - distance_matrix[i, j] = d - distance_matrix[j, i] = -d - return distance_matrix - - -############################################################ -# spike_order_values -############################################################ -def spike_order_values(spike_trains, indices=None, - interval=None, max_tau=None): - """ Computes the spike train symmetry value for each spike in each spike - train. - """ - 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." - # list of arrays for reulting asymmetry values - asymmetry_list = [np.zeros_like(st.spikes) for st in spike_trains] - # generate a list of possible index pairs - pairs = [(indices[i], j) for i in range(len(indices)) - for j in indices[i+1:]] - - # cython implementation - try: - from .cython.cython_directionality import \ - spike_order_values_cython as spike_order_values_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 - - for i, j in pairs: - a1, a2 = spike_order_values_impl(spike_trains[i].spikes, - spike_trains[j].spikes, - spike_trains[i].t_start, - spike_trains[i].t_end, - max_tau) - asymmetry_list[i] += a1 - asymmetry_list[j] += a2 - for a in asymmetry_list: - a /= len(spike_trains)-1 - return asymmetry_list - - -############################################################ -# optimal_asymmetry_order_from_matrix -############################################################ -def optimal_spike_train_order_from_matrix(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_spike_train_order -############################################################ -def optimal_spike_train_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_train_order_matrix(spike_trains, indices, interval, max_tau) - return optimal_spike_train_order_from_matrix(D, full_output) - - -############################################################ -# permutate_matrix -############################################################ -def permutate_matrix(D, p): - N = len(D) - D_p = np.empty_like(D) - 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/test/test_spike_delay_asymmetry.py b/test/test_spike_delay_asymmetry.py deleted file mode 100644 index 9de16e5..0000000 --- a/test/test_spike_delay_asymmetry.py +++ /dev/null @@ -1,40 +0,0 @@ -""" test_spike_delay_asymmetry.py - -Tests the asymmetry functions - -Copyright 2015, Mario Mulansky - -Distributed under the BSD License - -""" - -import numpy as np -from numpy.testing import assert_equal, assert_almost_equal, \ - assert_array_equal - -import pyspike as spk -from pyspike import SpikeTrain, DiscreteFunc - - -def test_profile(): - st1 = SpikeTrain([100, 200, 300], [0, 1000]) - st2 = SpikeTrain([105, 205, 300], [0, 1000]) - expected_x = np.array([0, 100, 105, 200, 205, 300, 1000]) - expected_y = np.array([1, 1, 1, 1, 1, 0, 0]) - expected_mp = np.array([1, 1, 1, 1, 1, 2, 2]) - - f = spk.drct.spike_train_order_profile(st1, st2) - - assert f.almost_equal(DiscreteFunc(expected_x, expected_y, expected_mp)) - assert_almost_equal(f.avrg(), 2.0/3.0) - assert_almost_equal(spk.drct.spike_train_order(st1, st2), 2.0/3.0) - assert_almost_equal(spk.drct.spike_train_order(st1, st2, normalize=False), - 4.0) - - st3 = SpikeTrain([105, 195, 500], [0, 1000]) - expected_x = np.array([0, 100, 105, 195, 200, 300, 500, 1000]) - expected_y = np.array([1, 1, 1, -1, -1, 0, 0, 0]) - expected_mp = np.array([1, 1, 1, 1, 1, 1, 1, 1]) - - f = spk.drct.spike_train_order_profile(st1, st3) - assert f.almost_equal(DiscreteFunc(expected_x, expected_y, expected_mp)) -- cgit v1.2.3 From 64b680238ac5bcafbf3ce67ed2a4d7487098b43d Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Mon, 14 Dec 2015 16:40:49 +0100 Subject: update version strings and Changelog --- Changelog | 6 ++++++ doc/conf.py | 4 ++-- setup.py | 2 +- 3 files changed, 9 insertions(+), 3 deletions(-) diff --git a/Changelog b/Changelog index 519dd3b..2be5e52 100644 --- a/Changelog +++ b/Changelog @@ -1,3 +1,9 @@ +PySpike v0.4: + * Python 3 support (thanks to Igor Gnatenko) + * list interface to SpikeTrain class + * disable_backend_warning property + * several bugfixes + PySpike v0.3: * addition of __version__ attribute * restructured docs, Readme now only contains basic examples diff --git a/doc/conf.py b/doc/conf.py index 8011ea9..7e13eae 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -64,9 +64,9 @@ copyright = u'2014-2015, Mario Mulansky' # built documents. # # The short X.Y version. -version = '0.3' +version = '0.4' # The full version, including alpha/beta/rc tags. -release = '0.3.0' +release = '0.4.0' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/setup.py b/setup.py index 62c1951..a1ab122 100644 --- a/setup.py +++ b/setup.py @@ -60,7 +60,7 @@ elif use_c: # c files are there, compile to binaries setup( name='pyspike', packages=find_packages(exclude=['doc']), - version='0.3.0', + version='0.4.0', cmdclass=cmdclass, ext_modules=ext_modules, include_dirs=[numpy.get_include()], -- cgit v1.2.3 From 9061f2a0c13134e53f937d730295a421fd671ea3 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Mon, 14 Dec 2015 16:50:55 +0100 Subject: removed directionality from __init__ and setup.py --- pyspike/__init__.py | 2 -- setup.py | 11 +++-------- 2 files changed, 3 insertions(+), 10 deletions(-) diff --git a/pyspike/__init__.py b/pyspike/__init__.py index 335b1d3..069090b 100644 --- a/pyspike/__init__.py +++ b/pyspike/__init__.py @@ -26,8 +26,6 @@ from .psth import psth from .spikes import load_spike_trains_from_txt, spike_train_from_string, \ merge_spike_trains, generate_poisson_spikes -from . import directionality as drct - # define the __version__ following # http://stackoverflow.com/questions/17583443 from pkg_resources import get_distribution, DistributionNotFound diff --git a/setup.py b/setup.py index a1ab122..8ef431a 100644 --- a/setup.py +++ b/setup.py @@ -23,8 +23,7 @@ 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") and \ - os.path.isfile("pyspike/directionality/cython/cython_directionality.c"): + os.path.isfile("pyspike/cython/cython_distances.c"): use_c = True else: use_c = False @@ -39,9 +38,7 @@ if use_cython: # Cython is available, compile .pyx -> .c 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"]) + ["pyspike/cython/cython_distances.pyx"]) ] cmdclass.update({'build_ext': build_ext}) elif use_c: # c files are there, compile to binaries @@ -51,9 +48,7 @@ elif use_c: # c files are there, compile to binaries 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"]) + ["pyspike/cython/cython_distances.c"]) ] # neither cython nor c files available -> automatic fall-back to python backend -- cgit v1.2.3