diff options
Diffstat (limited to 'pyspike')
-rw-r--r-- | pyspike/DiscreteFunc.py | 64 | ||||
-rw-r--r-- | pyspike/PieceWiseConstFunc.py | 13 | ||||
-rw-r--r-- | pyspike/PieceWiseLinFunc.py | 17 | ||||
-rw-r--r-- | pyspike/SpikeTrain.py | 15 | ||||
-rw-r--r-- | pyspike/__init__.py | 25 | ||||
-rw-r--r-- | pyspike/cython/cython_distances.pyx | 55 | ||||
-rw-r--r-- | pyspike/cython/cython_profiles.pyx | 61 | ||||
-rw-r--r-- | pyspike/cython/python_backend.py | 53 | ||||
-rw-r--r-- | pyspike/generic.py | 18 | ||||
-rw-r--r-- | pyspike/isi_distance.py | 14 | ||||
-rw-r--r-- | pyspike/psth.py | 2 | ||||
-rw-r--r-- | pyspike/spike_distance.py | 15 | ||||
-rw-r--r-- | pyspike/spike_sync.py | 12 | ||||
-rw-r--r-- | pyspike/spikes.py | 3 |
14 files changed, 233 insertions, 134 deletions
diff --git a/pyspike/DiscreteFunc.py b/pyspike/DiscreteFunc.py index 17153ee..fe97bc2 100644 --- a/pyspike/DiscreteFunc.py +++ b/pyspike/DiscreteFunc.py @@ -2,10 +2,11 @@ # Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net> # Distributed under the BSD License -from __future__ import print_function +from __future__ import absolute_import, print_function import numpy as np import collections +import pyspike ############################################################## @@ -79,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 @@ -136,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""" @@ -151,28 +151,32 @@ 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): + 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. @@ -185,7 +189,10 @@ class DiscreteFunc(object): :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. @@ -199,15 +206,16 @@ class DiscreteFunc(object): # 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: - 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.") # 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 = \ @@ -236,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/PieceWiseConstFunc.py b/pyspike/PieceWiseConstFunc.py index 2705443..5ce5f27 100644 --- a/pyspike/PieceWiseConstFunc.py +++ b/pyspike/PieceWiseConstFunc.py @@ -2,10 +2,11 @@ # Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net> # Distributed under the BSD License -from __future__ import print_function +from __future__ import absolute_import, print_function import numpy as np import collections +import pyspike ############################################################## @@ -188,14 +189,16 @@ 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: - 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 \ + 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 c0dd475..8145e63 100644 --- a/pyspike/PieceWiseLinFunc.py +++ b/pyspike/PieceWiseLinFunc.py @@ -2,10 +2,11 @@ # Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net> # Distributed under the BSD License -from __future__ import print_function +from __future__ import absolute_import, print_function import numpy as np import collections +import pyspike ############################################################## @@ -221,20 +222,22 @@ 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: - 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 \ + 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/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/__init__.py b/pyspike/__init__.py index 3e836bd..069090b 100644 --- a/pyspike/__init__.py +++ b/pyspike/__init__.py @@ -4,27 +4,28 @@ Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net> Distributed under the BSD License """ +from __future__ import absolute_import + __all__ = ["isi_distance", "spike_distance", "spike_sync", "psth", "spikes", "SpikeTrain", "PieceWiseConstFunc", "PieceWiseLinFunc", - "DiscreteFunc"] + "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 - # define the __version__ following # http://stackoverflow.com/questions/17583443 from pkg_resources import get_distribution, DistributionNotFound @@ -42,3 +43,5 @@ except DistributionNotFound: __version__ = 'Please install this project with setup.py' else: __version__ = _dist.version + +disable_backend_warning = False 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/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 diff --git a/pyspike/generic.py b/pyspike/generic.py index 41affcb..5ad06f1 100644 --- a/pyspike/generic.py +++ b/pyspike/generic.py @@ -7,6 +7,7 @@ Copyright 2015, Mario Mulansky <mario.mulansky@gmx.net> Distributed under the BSD License """ +from __future__ import division import numpy as np @@ -37,13 +38,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[L1/2:]) + dist_prof1 = divide_and_conquer(pairs1[:L1//2], + pairs1[L1//2:]) else: dist_prof1 = pair_distance_func(spike_trains[pairs1[0][0]], spike_trains[pairs1[0][1]]) L2 = len(pairs2) if L2 > 1: - dist_prof2 = divide_and_conquer(pairs2[:L2/2], pairs2[L2/2:]) + dist_prof2 = divide_and_conquer(pairs2[:L2//2], + pairs2[L2//2:]) else: dist_prof2 = pair_distance_func(spike_trains[pairs2[0][0]], spike_trains[pairs2[0][1]]) @@ -63,8 +66,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[:len(pairs)//2], + pairs[len(pairs)//2:]) else: avrg_dist = pair_distance_func(spike_trains[pairs[0][0]], spike_trains[pairs[0][1]]) @@ -135,12 +138,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 range(len(indices)) + for j in range(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/pyspike/isi_distance.py b/pyspike/isi_distance.py index 5ea555d..0ae7393 100644 --- a/pyspike/isi_distance.py +++ b/pyspike/isi_distance.py @@ -2,6 +2,9 @@ # Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net> # 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, \ _generic_distance_matrix @@ -31,14 +34,15 @@ 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: - 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 \ + from .cython.python_backend import isi_distance_python \ as isi_profile_impl times, values = isi_profile_impl(spike_train1.get_spikes_non_empty(), @@ -72,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/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/pyspike/spike_distance.py b/pyspike/spike_distance.py index ac2d260..e418283 100644 --- a/pyspike/spike_distance.py +++ b/pyspike/spike_distance.py @@ -2,6 +2,9 @@ # Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net> # 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, \ _generic_distance_matrix @@ -31,14 +34,15 @@ 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: - 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 - 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( @@ -54,7 +58,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. @@ -73,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 40d98d2..3dc29ff 100644 --- a/pyspike/spike_sync.py +++ b/pyspike/spike_sync.py @@ -3,8 +3,11 @@ # Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net> # Distributed under the BSD License +from __future__ import absolute_import + import numpy as np from functools import partial +import pyspike from pyspike import DiscreteFunc from pyspike.generic import _generic_profile_multi, _generic_distance_matrix @@ -36,14 +39,15 @@ 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: - 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 - from cython.python_backend import coincidence_python \ + from .cython.python_backend import coincidence_python \ as coincidence_profile_impl if max_tau is None: @@ -71,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 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 |