diff options
author | Mario Mulansky <mario.mulansky@gmx.net> | 2015-02-03 12:43:21 +0100 |
---|---|---|
committer | Mario Mulansky <mario.mulansky@gmx.net> | 2015-02-03 12:43:21 +0100 |
commit | 57b282e117c73913eed1f62dbe080544aadf4bbe (patch) | |
tree | 544ced40a0ef6d030af5f95f0f65dd975e8058c8 | |
parent | 877480b65c94019d25e44cec0e67876b96b4f418 (diff) | |
parent | bfee34894f21212bcef5fe271e701b1cf8c86fde (diff) |
Merge pull request #2 from mariomulansky/restructuring
Restructuring
-rw-r--r-- | pyspike/DiscreteFunc.py | 242 | ||||
-rw-r--r-- | pyspike/PieceWiseConstFunc.py | 168 | ||||
-rw-r--r-- | pyspike/PieceWiseLinFunc.py | 197 | ||||
-rw-r--r-- | pyspike/__init__.py | 23 | ||||
-rw-r--r-- | pyspike/cython/__init__.py | 0 | ||||
-rw-r--r-- | pyspike/cython/cython_add.pyx (renamed from pyspike/cython_add.pyx) | 0 | ||||
-rw-r--r-- | pyspike/cython/cython_distance.pyx (renamed from pyspike/cython_distance.pyx) | 0 | ||||
-rw-r--r-- | pyspike/cython/python_backend.py (renamed from pyspike/python_backend.py) | 0 | ||||
-rw-r--r-- | pyspike/distances.py | 504 | ||||
-rw-r--r-- | pyspike/function.py | 585 | ||||
-rw-r--r-- | pyspike/generic.py | 83 | ||||
-rw-r--r-- | pyspike/isi_distance.py | 134 | ||||
-rw-r--r-- | pyspike/spike_distance.py | 136 | ||||
-rw-r--r-- | pyspike/spike_sync.py | 136 | ||||
-rw-r--r-- | setup.py | 14 | ||||
-rw-r--r-- | test/test_function.py | 21 |
16 files changed, 1130 insertions, 1113 deletions
diff --git a/pyspike/DiscreteFunc.py b/pyspike/DiscreteFunc.py new file mode 100644 index 0000000..2283e03 --- /dev/null +++ b/pyspike/DiscreteFunc.py @@ -0,0 +1,242 @@ +""" +Class representing discrete functions. + +Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net> + +Distributed under the BSD License + +""" +from __future__ import print_function + +import numpy as np +import collections + + +############################################################## +# DiscreteFunc +############################################################## +class DiscreteFunc(object): + """ A class representing values defined on a discrete set of points. + """ + + def __init__(self, x, y, multiplicity): + """ Constructs the discrete function. + + :param x: array of length N defining the points at which the values are + defined. + :param y: array of length N degining the values at the points x. + :param multiplicity: array of length N defining the multiplicity of the + values. + """ + # convert parameters to arrays, also ensures copying + self.x = np.array(x) + self.y = np.array(y) + self.mp = np.array(multiplicity) + + def copy(self): + """ Returns a copy of itself + + :rtype: :class:`DiscreteFunc` + """ + return DiscreteFunc(self.x, self.y, self.mp) + + def almost_equal(self, other, decimal=14): + """ Checks if the function is equal to another function up to `decimal` + precision. + + :param other: another :class:`DiscreteFunc` + :returns: True if the two functions are equal up to `decimal` decimals, + False otherwise + :rtype: bool + """ + eps = 10.0**(-decimal) + return np.allclose(self.x, other.x, atol=eps, rtol=0.0) and \ + np.allclose(self.y, other.y, atol=eps, rtol=0.0) and \ + np.allclose(self.mp, other.mp, atol=eps, rtol=0.0) + + def get_plottable_data(self, averaging_window_size=0): + """ Returns two arrays containing x- and y-coordinates for plotting + the interval sequence. The optional parameter `averaging_window_size` + determines the size of an averaging window to smoothen the profile. If + this value is 0, no averaging is performed. + + :param averaging_window_size: size of the averaging window, default=0. + :returns: (x_plot, y_plot) containing plottable data + :rtype: pair of np.array + + Example:: + + x, y = f.get_plottable_data() + plt.plot(x, y, '-o', label="Discrete function") + """ + + if averaging_window_size > 0: + # for the averaged profile we have to take the multiplicity into + # account. values with higher multiplicity should be consider as if + # they appeared several times. Hence we can not know how many + # entries we have to consider to the left and right. Rather, we + # will iterate until some wanted multiplicity is reached. + + # the first value in self.mp contains the number of averaged + # profiles without any possible extra multiplicities + # (by implementation) + 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)): + + if self.mp[i] >= expected_mp: + # the current value contains already all the wanted + # multiplicity + y_plot[i] = self.y[i]/self.mp[i] + continue + + # first look to the right + y = self.y[i] + mp_r = self.mp[i] + j = i+1 + while j < len(y_plot): + if mp_r+self.mp[j] < expected_mp: + # if we still dont reach the required multiplicity + # we take the whole value + y += self.y[j] + mp_r += self.mp[j] + else: + # otherwise, just some fraction + y += self.y[j] * (expected_mp - mp_r)/self.mp[j] + mp_r += (expected_mp - mp_r) + break + j += 1 + + # same story to the left + mp_l = self.mp[i] + j = i-1 + while j >= 0: + if mp_l+self.mp[j] < expected_mp: + y += self.y[j] + mp_l += self.mp[j] + else: + y += self.y[j] * (expected_mp - mp_l)/self.mp[j] + mp_l += (expected_mp - mp_l) + break + j -= 1 + y_plot[i] = y/(mp_l+mp_r-self.mp[i]) + return 1.0*self.x, y_plot + + else: # k = 0 + + return 1.0*self.x, 1.0*self.y/self.mp + + def integral(self, interval=None): + """ Returns the integral over the given interval. For the discrete + function, this amounts to the sum over all values divided by the total + multiplicity. + + :param interval: integration interval given as a pair of floats, or a + sequence of pairs in case of multiple intervals, if + None the integral over the whole function is computed. + :type interval: Pair, sequence of pairs, or None. + :returns: the integral + :rtype: float + """ + + def get_indices(ival): + """ Retuns the indeces surrounding the given interval""" + start_ind = np.searchsorted(self.x, ival[0], side='right') + end_ind = np.searchsorted(self.x, ival[1], side='left') + assert start_ind > 0 and end_ind < len(self.x), \ + "Invalid averaging interval" + return start_ind, end_ind + + 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])) + else: + value = 0.0 + multiplicity = 0.0 + 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]) + return value/multiplicity + + def avrg(self, interval=None): + """ Computes the average of the interval sequence: + :math:`a = 1/N sum f_n ` where N is the number of intervals. + + :param interval: averaging interval given as a pair of floats, a + sequence of pairs for averaging multiple intervals, or + None, if None the average over the whole function is + computed. + :type interval: Pair, sequence of pairs, or None. + :returns: the average a. + :rtype: float + """ + return self.integral(interval) + + def add(self, f): + """ Adds another `DiscreteFunc` function to this function. + Note: only functions defined on the same interval can be summed. + + :param f: :class:`DiscreteFunc` function to be added. + :rtype: None + """ + assert self.x[0] == f.x[0], "The functions have different intervals" + assert self.x[-1] == f.x[-1], "The functions have different intervals" + + # cython version + try: + 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 \ +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 \ + add_discrete_function_impl + + self.x, self.y, self.mp = \ + add_discrete_function_impl(self.x, self.y, self.mp, + f.x, f.y, f.mp) + + def mul_scalar(self, fac): + """ Multiplies the function with a scalar value + + :param fac: Value to multiply + :type fac: double + :rtype: None + """ + self.y *= fac + + +def average_profile(profiles): + """ Computes the average profile from the given ISI- or SPIKE-profiles. + + :param profiles: list of :class:`PieceWiseConstFunc` or + :class:`PieceWiseLinFunc` representing ISI- or + SPIKE-profiles to be averaged. + :returns: the averages profile :math:`<S_{isi}>` or :math:`<S_{spike}>`. + :rtype: :class:`PieceWiseConstFunc` or :class:`PieceWiseLinFunc` + """ + assert len(profiles) > 1 + + avrg_profile = profiles[0].copy() + for i in xrange(1, len(profiles)): + avrg_profile.add(profiles[i]) + avrg_profile.mul_scalar(1.0/len(profiles)) # normalize + + return avrg_profile diff --git a/pyspike/PieceWiseConstFunc.py b/pyspike/PieceWiseConstFunc.py new file mode 100644 index 0000000..dc57ab1 --- /dev/null +++ b/pyspike/PieceWiseConstFunc.py @@ -0,0 +1,168 @@ +""" +Class representing piece-wise constant functions. + +Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net> + +Distributed under the BSD License + +""" +from __future__ import print_function + +import numpy as np +import collections + + +############################################################## +# PieceWiseConstFunc +############################################################## +class PieceWiseConstFunc(object): + """ A class representing a piece-wise constant function. """ + + def __init__(self, x, y): + """ Constructs the piece-wise const function. + + :param x: array of length N+1 defining the edges of the intervals of + the pwc function. + :param y: array of length N defining the function values at the + intervals. + """ + # convert parameters to arrays, also ensures copying + self.x = np.array(x) + self.y = np.array(y) + + def copy(self): + """ Returns a copy of itself + + :rtype: :class:`PieceWiseConstFunc` + """ + return PieceWiseConstFunc(self.x, self.y) + + def almost_equal(self, other, decimal=14): + """ Checks if the function is equal to another function up to `decimal` + precision. + + :param other: another :class:`PieceWiseConstFunc` + :returns: True if the two functions are equal up to `decimal` decimals, + False otherwise + :rtype: bool + """ + eps = 10.0**(-decimal) + return np.allclose(self.x, other.x, atol=eps, rtol=0.0) and \ + np.allclose(self.y, other.y, atol=eps, rtol=0.0) + + def get_plottable_data(self): + """ Returns two arrays containing x- and y-coordinates for immeditate + plotting of the piece-wise function. + + :returns: (x_plot, y_plot) containing plottable data + :rtype: pair of np.array + + Example:: + + x, y = f.get_plottable_data() + plt.plot(x, y, '-o', label="Piece-wise const function") + """ + + x_plot = np.empty(2*len(self.x)-2) + x_plot[0] = self.x[0] + x_plot[1::2] = self.x[1:] + x_plot[2::2] = self.x[1:-1] + y_plot = np.empty(2*len(self.y)) + y_plot[::2] = self.y + y_plot[1::2] = self.y + + return x_plot, y_plot + + def integral(self, interval=None): + """ Returns the integral over the given interval. + + :param interval: integration interval given as a pair of floats, if + None the integral over the whole function is computed. + :type interval: Pair of floats or None. + :returns: the integral + :rtype: float + """ + if interval is None: + # no interval given, integrate over the whole spike train + a = np.sum((self.x[1:]-self.x[:-1]) * self.y) + else: + # find the indices corresponding to the interval + start_ind = np.searchsorted(self.x, interval[0], side='right') + end_ind = np.searchsorted(self.x, interval[1], side='left')-1 + assert start_ind > 0 and end_ind < len(self.x), \ + "Invalid averaging interval" + # first the contribution from between the indices + a = np.sum((self.x[start_ind+1:end_ind+1] - + self.x[start_ind:end_ind]) * + self.y[start_ind:end_ind]) + # correction from start to first index + a += (self.x[start_ind]-interval[0]) * self.y[start_ind-1] + # correction from last index to end + a += (interval[1]-self.x[end_ind]) * self.y[end_ind] + return a + + def avrg(self, interval=None): + """ Computes the average of the piece-wise const function: + :math:`a = 1/T int_0^T f(x) dx` where T is the length of the interval. + + :param interval: averaging interval given as a pair of floats, a + sequence of pairs for averaging multiple intervals, or + None, if None the average over the whole function is + computed. + :type interval: Pair, sequence of pairs, or None. + :returns: the average a. + :rtype: float + """ + if interval is None: + # no interval given, average over the whole spike train + return self.integral() / (self.x[-1]-self.x[0]) + + # 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): + # just one interval + a = self.integral(interval) / (interval[1]-interval[0]) + else: + # several intervals + a = 0.0 + int_length = 0.0 + for ival in interval: + a += self.integral(ival) + int_length += ival[1] - ival[0] + a /= int_length + return a + + def add(self, f): + """ Adds another PieceWiseConst function to this function. + Note: only functions defined on the same interval can be summed. + + :param f: :class:`PieceWiseConstFunc` function to be added. + :rtype: None + """ + assert self.x[0] == f.x[0], "The functions have different intervals" + assert self.x[-1] == f.x[-1], "The functions have different intervals" + + # cython version + try: + 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'! \ +\n Falling back to slow python backend.") + # use python backend + 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) + + def mul_scalar(self, fac): + """ Multiplies the function with a scalar value + + :param fac: Value to multiply + :type fac: double + :rtype: None + """ + self.y *= fac diff --git a/pyspike/PieceWiseLinFunc.py b/pyspike/PieceWiseLinFunc.py new file mode 100644 index 0000000..bc0aa2a --- /dev/null +++ b/pyspike/PieceWiseLinFunc.py @@ -0,0 +1,197 @@ +""" +Class representing piece-wise linear functions. + +Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net> + +Distributed under the BSD License + +""" +from __future__ import print_function + +import numpy as np +import collections + + +############################################################## +# PieceWiseLinFunc +############################################################## +class PieceWiseLinFunc: + """ A class representing a piece-wise linear function. """ + + def __init__(self, x, y1, y2): + """ Constructs the piece-wise linear function. + + :param x: array of length N+1 defining the edges of the intervals of + the pwc function. + :param y1: array of length N defining the function values at the left + of the intervals. + :param y2: array of length N defining the function values at the right + of the intervals. + """ + # convert to array, which also ensures copying + self.x = np.array(x) + self.y1 = np.array(y1) + self.y2 = np.array(y2) + + def copy(self): + """ Returns a copy of itself + + :rtype: :class:`PieceWiseLinFunc` + """ + return PieceWiseLinFunc(self.x, self.y1, self.y2) + + def almost_equal(self, other, decimal=14): + """ Checks if the function is equal to another function up to `decimal` + precision. + + :param other: another :class:`PieceWiseLinFunc` + :returns: True if the two functions are equal up to `decimal` decimals, + False otherwise + :rtype: bool + """ + eps = 10.0**(-decimal) + return np.allclose(self.x, other.x, atol=eps, rtol=0.0) and \ + np.allclose(self.y1, other.y1, atol=eps, rtol=0.0) and \ + np.allclose(self.y2, other.y2, atol=eps, rtol=0.0) + + def get_plottable_data(self): + """ Returns two arrays containing x- and y-coordinates for immeditate + plotting of the piece-wise function. + + :returns: (x_plot, y_plot) containing plottable data + :rtype: pair of np.array + + Example:: + + x, y = f.get_plottable_data() + plt.plot(x, y, '-o', label="Piece-wise const function") + """ + x_plot = np.empty(2*len(self.x)-2) + x_plot[0] = self.x[0] + x_plot[1::2] = self.x[1:] + x_plot[2::2] = self.x[1:-1] + y_plot = np.empty_like(x_plot) + y_plot[0::2] = self.y1 + y_plot[1::2] = self.y2 + return x_plot, y_plot + + def integral(self, interval=None): + """ Returns the integral over the given interval. + + :param interval: integration interval given as a pair of floats, if + None the integral over the whole function is computed. + :type interval: Pair of floats or None. + :returns: the integral + :rtype: float + """ + + def intermediate_value(x0, x1, y0, y1, x): + """ computes the intermediate value of a linear function """ + return y0 + (y1-y0)*(x-x0)/(x1-x0) + + if interval is None: + # no interval given, integrate over the whole spike train + integral = np.sum((self.x[1:]-self.x[:-1]) * 0.5*(self.y1+self.y2)) + else: + # find the indices corresponding to the interval + start_ind = np.searchsorted(self.x, interval[0], side='right') + end_ind = np.searchsorted(self.x, interval[1], side='left')-1 + assert start_ind > 0 and end_ind < len(self.x), \ + "Invalid averaging interval" + # first the contribution from between the indices + integral = np.sum((self.x[start_ind+1:end_ind+1] - + self.x[start_ind:end_ind]) * + 0.5*(self.y1[start_ind:end_ind] + + self.y2[start_ind:end_ind])) + # correction from start to first index + integral += (self.x[start_ind]-interval[0]) * 0.5 * \ + (self.y2[start_ind-1] + + intermediate_value(self.x[start_ind-1], + self.x[start_ind], + self.y1[start_ind-1], + self.y2[start_ind-1], + interval[0] + )) + # correction from last index to end + integral += (interval[1]-self.x[end_ind]) * 0.5 * \ + (self.y1[end_ind] + + intermediate_value(self.x[end_ind], self.x[end_ind+1], + self.y1[end_ind], self.y2[end_ind], + interval[1] + )) + return integral + + def avrg(self, interval=None): + """ Computes the average of the piece-wise linear function: + :math:`a = 1/T int_0^T f(x) dx` where T is the length of the interval. + + :param interval: averaging interval given as a pair of floats, a + sequence of pairs for averaging multiple intervals, or + None, if None the average over the whole function is + computed. + :type interval: Pair, sequence of pairs, or None. + :returns: the average a. + :rtype: float + + """ + + if interval is None: + # no interval given, average over the whole spike train + return self.integral() / (self.x[-1]-self.x[0]) + + # 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): + # just one interval + a = self.integral(interval) / (interval[1]-interval[0]) + else: + # several intervals + a = 0.0 + int_length = 0.0 + for ival in interval: + a += self.integral(ival) + int_length += ival[1] - ival[0] + a /= int_length + return a + + def add(self, f): + """ Adds another PieceWiseLin function to this function. + Note: only functions defined on the same interval can be summed. + + :param f: :class:`PieceWiseLinFunc` function to be added. + :rtype: None + """ + assert self.x[0] == f.x[0], "The functions have different intervals" + assert self.x[-1] == f.x[-1], "The functions have different intervals" + + # python implementation + # 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 \ + 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.") + # use python backend + 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( + self.x, self.y1, self.y2, f.x, f.y1, f.y2) + + def mul_scalar(self, fac): + """ Multiplies the function with a scalar value + + :param fac: Value to multiply + :type fac: double + :rtype: None + """ + self.y1 *= fac + self.y2 *= fac diff --git a/pyspike/__init__.py b/pyspike/__init__.py index 55687e6..945dd4e 100644 --- a/pyspike/__init__.py +++ b/pyspike/__init__.py @@ -1,18 +1,23 @@ """ -Copyright 2014, Mario Mulansky <mario.mulansky@gmx.net> +Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net> Distributed under the BSD License """ -__all__ = ["function", "distances", "spikes"] +__all__ = ["isi_distance", "spike_distance", "spike_sync", + "spikes", "PieceWiseConstFunc", "PieceWiseLinFunc", + "DiscreteFunc"] -from function import PieceWiseConstFunc, PieceWiseLinFunc, \ - DiscreteFunction, average_profile -from distances import isi_profile, isi_distance, \ - spike_profile, spike_distance, \ - spike_sync_profile, spike_sync, \ - isi_profile_multi, isi_distance_multi, isi_distance_matrix, \ - spike_profile_multi, spike_distance_multi, spike_distance_matrix, \ +from PieceWiseConstFunc import PieceWiseConstFunc +from PieceWiseLinFunc import PieceWiseLinFunc +from DiscreteFunc import DiscreteFunc + +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,\ + spike_distance_multi, spike_distance_matrix +from spike_sync import spike_sync_profile, spike_sync,\ spike_sync_profile_multi, spike_sync_multi, spike_sync_matrix + from spikes import add_auxiliary_spikes, load_spike_trains_from_txt, \ spike_train_from_string, merge_spike_trains, generate_poisson_spikes diff --git a/pyspike/cython/__init__.py b/pyspike/cython/__init__.py new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/pyspike/cython/__init__.py diff --git a/pyspike/cython_add.pyx b/pyspike/cython/cython_add.pyx index ac64005..ac64005 100644 --- a/pyspike/cython_add.pyx +++ b/pyspike/cython/cython_add.pyx diff --git a/pyspike/cython_distance.pyx b/pyspike/cython/cython_distance.pyx index 489aab9..489aab9 100644 --- a/pyspike/cython_distance.pyx +++ b/pyspike/cython/cython_distance.pyx diff --git a/pyspike/python_backend.py b/pyspike/cython/python_backend.py index 481daf9..481daf9 100644 --- a/pyspike/python_backend.py +++ b/pyspike/cython/python_backend.py diff --git a/pyspike/distances.py b/pyspike/distances.py deleted file mode 100644 index 5476b6f..0000000 --- a/pyspike/distances.py +++ /dev/null @@ -1,504 +0,0 @@ -""" distances.py - -Module containing several functions to compute spike distances - -Copyright 2014, Mario Mulansky <mario.mulansky@gmx.net> - -Distributed under the BSD License -""" - -import numpy as np -import threading -from functools import partial - -from pyspike import PieceWiseConstFunc, PieceWiseLinFunc, DiscreteFunction - - -############################################################ -# isi_profile -############################################################ -def isi_profile(spikes1, spikes2): - """ Computes the isi-distance profile :math:`S_{isi}(t)` of the two given - spike trains. Retruns the profile as a PieceWiseConstFunc object. The S_isi - values are defined positive S_isi(t)>=0. The spike trains are expected - to have auxiliary spikes at the beginning and end of the interval. Use the - function add_auxiliary_spikes to add those spikes to the spike train. - - :param spikes1: ordered array of spike times with auxiliary spikes. - :param spikes2: ordered array of spike times with auxiliary spikes. - :returns: The isi-distance profile :math:`S_{isi}(t)` - :rtype: :class:`pyspike.function.PieceWiseConstFunc` - - """ - # check for auxiliary spikes - first and last spikes should be identical - assert spikes1[0] == spikes2[0], \ - "Given spike trains seems not to have auxiliary spikes!" - assert spikes1[-1] == spikes2[-1], \ - "Given spike trains seems not to have auxiliary spikes!" - - # load cython implementation - try: - from cython_distance import isi_distance_cython as isi_distance_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 \ -Falling back to slow python backend.") - # use python backend - from python_backend import isi_distance_python as isi_distance_impl - - times, values = isi_distance_impl(spikes1, spikes2) - return PieceWiseConstFunc(times, values) - - -############################################################ -# isi_distance -############################################################ -def isi_distance(spikes1, spikes2, interval=None): - """ Computes the isi-distance I of the given spike trains. The - isi-distance is the integral over the isi distance profile - :math:`S_{isi}(t)`: - - .. math:: I = \int_{T_0}^{T_1} S_{isi}(t) dt. - - :param spikes1: ordered array of spike times with auxiliary spikes. - :param spikes2: ordered array of spike times with auxiliary spikes. - :param interval: averaging interval given as a pair of floats (T0, T1), - if None the average over the whole function is computed. - :type interval: Pair of floats or None. - :returns: The isi-distance I. - :rtype: double - """ - return isi_profile(spikes1, spikes2).avrg(interval) - - -############################################################ -# spike_profile -############################################################ -def spike_profile(spikes1, spikes2): - """ Computes the spike-distance profile S_spike(t) of the two given spike - trains. Returns the profile as a PieceWiseLinFunc object. The S_spike - values are defined positive S_spike(t)>=0. The spike trains are expected to - have auxiliary spikes at the beginning and end of the interval. Use the - function add_auxiliary_spikes to add those spikes to the spike train. - - :param spikes1: ordered array of spike times with auxiliary spikes. - :param spikes2: ordered array of spike times with auxiliary spikes. - :returns: The spike-distance profile :math:`S_{spike}(t)`. - :rtype: :class:`pyspike.function.PieceWiseLinFunc` - - """ - # check for auxiliary spikes - first and last spikes should be identical - assert spikes1[0] == spikes2[0], \ - "Given spike trains seems not to have auxiliary spikes!" - assert spikes1[-1] == spikes2[-1], \ - "Given spike trains seems not to have auxiliary spikes!" - - # cython implementation - try: - from cython_distance import spike_distance_cython \ - as spike_distance_impl - except ImportError: - 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 python_backend import spike_distance_python as spike_distance_impl - - times, y_starts, y_ends = spike_distance_impl(spikes1, spikes2) - return PieceWiseLinFunc(times, y_starts, y_ends) - - -############################################################ -# spike_distance -############################################################ -def spike_distance(spikes1, spikes2, interval=None): - """ Computes the spike-distance S of the given spike trains. The - spike-distance is the integral over the isi distance profile S_spike(t): - - .. math:: S = \int_{T_0}^{T_1} S_{spike}(t) dt. - - :param spikes1: ordered array of spike times with auxiliary spikes. - :param spikes2: ordered array of spike times with auxiliary spikes. - :param interval: averaging interval given as a pair of floats (T0, T1), - if None the average over the whole function is computed. - :type interval: Pair of floats or None. - :returns: The spike-distance. - :rtype: double - - """ - return spike_profile(spikes1, spikes2).avrg(interval) - - -############################################################ -# spike_sync_profile -############################################################ -def spike_sync_profile(spikes1, spikes2): - """ Computes the spike-synchronization profile S_sync(t) of the two given - spike trains. Returns the profile as a DiscreteFunction object. The S_sync - values are either 1 or 0, indicating the presence or absence of a - coincidence. The spike trains are expected to have auxiliary spikes at the - beginning and end of the interval. Use the function add_auxiliary_spikes to - add those spikes to the spike train. - - :param spikes1: ordered array of spike times with auxiliary spikes. - :param spikes2: ordered array of spike times with auxiliary spikes. - :returns: The spike-distance profile :math:`S_{sync}(t)`. - :rtype: :class:`pyspike.function.DiscreteFunction` - - """ - - # cython implementation - try: - from cython_distance import coincidence_cython \ - as coincidence_impl - except ImportError: - 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 python_backend import coincidence_python \ - as coincidence_impl - - times, coincidences, multiplicity = coincidence_impl(spikes1, spikes2) - - return DiscreteFunction(times, coincidences, multiplicity) - - -############################################################ -# spike_sync -############################################################ -def spike_sync(spikes1, spikes2, interval=None): - """ Computes the spike synchronization value SYNC of the given spike - trains. The spike synchronization value is the computed as the total number - of coincidences divided by the total number of spikes: - - .. math:: SYNC = \sum_n C_n / N. - - :param spikes1: ordered array of spike times with auxiliary spikes. - :param spikes2: ordered array of spike times with auxiliary spikes. - :param interval: averaging interval given as a pair of floats (T0, T1), - if None the average over the whole function is computed. - :type interval: Pair of floats or None. - :returns: The spike synchronization value. - :rtype: double - """ - return spike_sync_profile(spikes1, spikes2).avrg(interval) - - -############################################################ -# _generic_profile_multi -############################################################ -def _generic_profile_multi(spike_trains, pair_distance_func, indices=None): - """ Internal implementation detail, don't call this function directly, - use isi_profile_multi or spike_profile_multi instead. - - Computes the multi-variate distance for a set of spike-trains using the - pair_dist_func to compute pair-wise distances. That is it computes the - average distance of all pairs of spike-trains: - :math:`S(t) = 2/((N(N-1)) sum_{<i,j>} S_{i,j}`, - where the sum goes over all pairs <i,j>. - Args: - - spike_trains: list of spike trains - - pair_distance_func: function computing the distance of two spike trains - - indices: list of indices defining which spike trains to use, - if None all given spike trains are used (default=None) - Returns: - - The averaged multi-variate distance of all pairs - """ - if indices is None: - indices = np.arange(len(spike_trains)) - indices = np.array(indices) - # 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:]] - # start with first pair - (i, j) = pairs[0] - average_dist = pair_distance_func(spike_trains[i], spike_trains[j]) - for (i, j) in pairs[1:]: - current_dist = pair_distance_func(spike_trains[i], spike_trains[j]) - average_dist.add(current_dist) # add to the average - return average_dist, len(pairs) - - -############################################################ -# multi_distance_par -############################################################ -def _multi_distance_par(spike_trains, pair_distance_func, indices=None): - """ parallel implementation of the multi-distance. Not currently used as - it does not improve the performance. - """ - - num_threads = 2 - lock = threading.Lock() - - def run(spike_trains, index_pairs, average_dist): - (i, j) = index_pairs[0] - # print(i,j) - this_avrg = pair_distance_func(spike_trains[i], spike_trains[j]) - for (i, j) in index_pairs[1:]: - # print(i,j) - current_dist = pair_distance_func(spike_trains[i], spike_trains[j]) - this_avrg.add(current_dist) - with lock: - average_dist.add(this_avrg) - - 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:]] - num_pairs = len(pairs) - - # start with first pair - (i, j) = pairs[0] - average_dist = pair_distance_func(spike_trains[i], spike_trains[j]) - # remove the one we already computed - pairs = pairs[1:] - # distribute the rest into num_threads pieces - clustered_pairs = [pairs[n::num_threads] for n in xrange(num_threads)] - - threads = [] - for pairs in clustered_pairs: - t = threading.Thread(target=run, args=(spike_trains, pairs, - average_dist)) - threads.append(t) - t.start() - for t in threads: - t.join() - average_dist.mul_scalar(1.0/num_pairs) # normalize - return average_dist - - -############################################################ -# isi_profile_multi -############################################################ -def isi_profile_multi(spike_trains, indices=None): - """ computes the multi-variate isi distance profile for a set of spike - trains. That is the average isi-distance of all pairs of spike-trains: - S_isi(t) = 2/((N(N-1)) sum_{<i,j>} S_{isi}^{i,j}, - where the sum goes over all pairs <i,j> - - :param spike_trains: list of spike trains - :param indices: list of indices defining which spike trains to use, - if None all given spike trains are used (default=None) - :type state: list or None - :returns: The averaged isi profile :math:`<S_{isi}>(t)` - :rtype: :class:`pyspike.function.PieceWiseConstFunc` - """ - average_dist, M = _generic_profile_multi(spike_trains, isi_profile, - indices) - average_dist.mul_scalar(1.0/M) # normalize - return average_dist - - -############################################################ -# isi_distance_multi -############################################################ -def isi_distance_multi(spike_trains, indices=None, interval=None): - """ computes the multi-variate isi-distance for a set of spike-trains. - That is the time average of the multi-variate spike profile: - I = \int_0^T 2/((N(N-1)) sum_{<i,j>} S_{isi}^{i,j}, - where the sum goes over all pairs <i,j> - - :param spike_trains: list of spike trains - :param indices: list of indices defining which spike trains to use, - if None all given spike trains are used (default=None) - :param interval: averaging interval given as a pair of floats, if None - the average over the whole function is computed. - :type interval: Pair of floats or None. - :returns: The time-averaged isi distance :math:`I` - :rtype: double - """ - return isi_profile_multi(spike_trains, indices).avrg(interval) - - -############################################################ -# spike_profile_multi -############################################################ -def spike_profile_multi(spike_trains, indices=None): - """ Computes the multi-variate spike distance profile for a set of spike - trains. That is the average spike-distance of all pairs of spike-trains: - :math:`S_spike(t) = 2/((N(N-1)) sum_{<i,j>} S_{spike}^{i, j}`, - where the sum goes over all pairs <i,j> - - :param spike_trains: list of spike trains - :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 - :returns: The averaged spike profile :math:`<S_{spike}>(t)` - :rtype: :class:`pyspike.function.PieceWiseLinFunc` - - """ - average_dist, M = _generic_profile_multi(spike_trains, spike_profile, - indices) - average_dist.mul_scalar(1.0/M) # normalize - return average_dist - - -############################################################ -# spike_distance_multi -############################################################ -def spike_distance_multi(spike_trains, indices=None, interval=None): - """ Computes the multi-variate spike distance for a set of spike trains. - That is the time average of the multi-variate spike profile: - S_{spike} = \int_0^T 2/((N(N-1)) sum_{<i,j>} S_{spike}^{i, j} dt - where the sum goes over all pairs <i,j> - - :param spike_trains: list of spike trains - :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 interval: averaging interval given as a pair of floats, if None - the average over the whole function is computed. - :type interval: Pair of floats or None. - :returns: The averaged spike distance S. - :rtype: double - """ - return spike_profile_multi(spike_trains, indices).avrg(interval) - - -############################################################ -# spike_profile_multi -############################################################ -def spike_sync_profile_multi(spike_trains, indices=None): - """ Computes the multi-variate spike synchronization profile for a set of - spike trains. For each spike in the set of spike trains, the multi-variate - profile is defined as the number of coincidences 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 spike trains - :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 - :returns: The multi-variate spike sync profile :math:`<S_{sync}>(t)` - :rtype: :class:`pyspike.function.DiscreteFunction` - - """ - prof_func = partial(spike_sync_profile) - average_dist, M = _generic_profile_multi(spike_trains, prof_func, - indices) - # average_dist.mul_scalar(1.0/M) # no normalization here! - return average_dist - - -############################################################ -# spike_distance_multi -############################################################ -def spike_sync_multi(spike_trains, indices=None, interval=None): - """ Computes the multi-variate spike synchronization value for a set of - spike trains. - - :param spike_trains: list of spike trains - :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 interval: averaging interval given as a pair of floats, if None - the average over the whole function is computed. - :type interval: Pair of floats or None. - :returns: The multi-variate spike synchronization value SYNC. - :rtype: double - """ - return spike_sync_profile_multi(spike_trains, indices).avrg(interval) - - -############################################################ -# generic_distance_matrix -############################################################ -def _generic_distance_matrix(spike_trains, dist_function, - indices=None, interval=None): - """ Internal implementation detail. Don't use this function directly. - Instead use isi_distance_matrix or spike_distance_matrix. - Computes the time averaged distance of all pairs of spike-trains. - Args: - - spike_trains: list of spike trains - - indices: list of indices defining which spike-trains to use - if None all given spike-trains are used (default=None) - Return: - - a 2D array of size len(indices)*len(indices) containing the average - pair-wise distance - """ - 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 = dist_function(spike_trains[i], spike_trains[j], interval) - distance_matrix[i, j] = d - distance_matrix[j, i] = d - return distance_matrix - - -############################################################ -# isi_distance_matrix -############################################################ -def isi_distance_matrix(spike_trains, indices=None, interval=None): - """ Computes the time averaged isi-distance of all pairs of spike-trains. - - :param spike_trains: list of spike trains - :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 interval: averaging interval given as a pair of floats, if None - the average over the whole function is computed. - :type interval: Pair of floats or None. - :returns: 2D array with the pair wise time average isi distances - :math:`I_{ij}` - :rtype: np.array - """ - return _generic_distance_matrix(spike_trains, isi_distance, - indices, interval) - - -############################################################ -# spike_distance_matrix -############################################################ -def spike_distance_matrix(spike_trains, indices=None, interval=None): - """ Computes the time averaged spike-distance of all pairs of spike-trains. - - :param spike_trains: list of spike trains - :param 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 interval: averaging interval given as a pair of floats, if None - the average over the whole function is computed. - :type interval: Pair of floats or None. - :returns: 2D array with the pair wise time average spike distances - :math:`S_{ij}` - :rtype: np.array - """ - return _generic_distance_matrix(spike_trains, spike_distance, - indices, interval) - - -############################################################ -# spike_sync_matrix -############################################################ -def spike_sync_matrix(spike_trains, indices=None, interval=None): - """ Computes the overall spike-synchronization value of all pairs of - spike-trains. - - :param spike_trains: list of spike trains - :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 interval: averaging interval given as a pair of floats, if None - the average over the whole function is computed. - :type interval: Pair of floats or None. - :returns: 2D array with the pair wise time spike synchronization values - :math:`SYNC_{ij}` - :rtype: np.array - """ - return _generic_distance_matrix(spike_trains, spike_sync, - indices, interval) diff --git a/pyspike/function.py b/pyspike/function.py deleted file mode 100644 index 047c88a..0000000 --- a/pyspike/function.py +++ /dev/null @@ -1,585 +0,0 @@ -""" function.py - -Module containing classes representing piece-wise constant and piece-wise -linear functions. - -Copyright 2014, Mario Mulansky <mario.mulansky@gmx.net> - -Distributed under the BSD License - -""" -from __future__ import print_function - -import numpy as np -import collections - - -############################################################## -# PieceWiseConstFunc -############################################################## -class PieceWiseConstFunc(object): - """ A class representing a piece-wise constant function. """ - - def __init__(self, x, y): - """ Constructs the piece-wise const function. - - :param x: array of length N+1 defining the edges of the intervals of - the pwc function. - :param y: array of length N defining the function values at the - intervals. - """ - # convert parameters to arrays, also ensures copying - self.x = np.array(x) - self.y = np.array(y) - - def copy(self): - """ Returns a copy of itself - - :rtype: :class:`PieceWiseConstFunc` - """ - return PieceWiseConstFunc(self.x, self.y) - - def almost_equal(self, other, decimal=14): - """ Checks if the function is equal to another function up to `decimal` - precision. - - :param other: another :class:`PieceWiseConstFunc` - :returns: True if the two functions are equal up to `decimal` decimals, - False otherwise - :rtype: bool - """ - eps = 10.0**(-decimal) - return np.allclose(self.x, other.x, atol=eps, rtol=0.0) and \ - np.allclose(self.y, other.y, atol=eps, rtol=0.0) - - def get_plottable_data(self): - """ Returns two arrays containing x- and y-coordinates for immeditate - plotting of the piece-wise function. - - :returns: (x_plot, y_plot) containing plottable data - :rtype: pair of np.array - - Example:: - - x, y = f.get_plottable_data() - plt.plot(x, y, '-o', label="Piece-wise const function") - """ - - x_plot = np.empty(2*len(self.x)-2) - x_plot[0] = self.x[0] - x_plot[1::2] = self.x[1:] - x_plot[2::2] = self.x[1:-1] - y_plot = np.empty(2*len(self.y)) - y_plot[::2] = self.y - y_plot[1::2] = self.y - - return x_plot, y_plot - - def integral(self, interval=None): - """ Returns the integral over the given interval. - - :param interval: integration interval given as a pair of floats, if - None the integral over the whole function is computed. - :type interval: Pair of floats or None. - :returns: the integral - :rtype: float - """ - if interval is None: - # no interval given, integrate over the whole spike train - a = np.sum((self.x[1:]-self.x[:-1]) * self.y) - else: - # find the indices corresponding to the interval - start_ind = np.searchsorted(self.x, interval[0], side='right') - end_ind = np.searchsorted(self.x, interval[1], side='left')-1 - assert start_ind > 0 and end_ind < len(self.x), \ - "Invalid averaging interval" - # first the contribution from between the indices - a = np.sum((self.x[start_ind+1:end_ind+1] - - self.x[start_ind:end_ind]) * - self.y[start_ind:end_ind]) - # correction from start to first index - a += (self.x[start_ind]-interval[0]) * self.y[start_ind-1] - # correction from last index to end - a += (interval[1]-self.x[end_ind]) * self.y[end_ind] - return a - - def avrg(self, interval=None): - """ Computes the average of the piece-wise const function: - :math:`a = 1/T int_0^T f(x) dx` where T is the length of the interval. - - :param interval: averaging interval given as a pair of floats, a - sequence of pairs for averaging multiple intervals, or - None, if None the average over the whole function is - computed. - :type interval: Pair, sequence of pairs, or None. - :returns: the average a. - :rtype: float - """ - if interval is None: - # no interval given, average over the whole spike train - return self.integral() / (self.x[-1]-self.x[0]) - - # 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): - # just one interval - a = self.integral(interval) / (interval[1]-interval[0]) - else: - # several intervals - a = 0.0 - int_length = 0.0 - for ival in interval: - a += self.integral(ival) - int_length += ival[1] - ival[0] - a /= int_length - return a - - def add(self, f): - """ Adds another PieceWiseConst function to this function. - Note: only functions defined on the same interval can be summed. - - :param f: :class:`PieceWiseConstFunc` function to be added. - :rtype: None - """ - assert self.x[0] == f.x[0], "The functions have different intervals" - assert self.x[-1] == f.x[-1], "The functions have different intervals" - - # cython version - try: - from 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'! \ -\n Falling back to slow python backend.") - # use python backend - from 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) - - def mul_scalar(self, fac): - """ Multiplies the function with a scalar value - - :param fac: Value to multiply - :type fac: double - :rtype: None - """ - self.y *= fac - - -############################################################## -# PieceWiseLinFunc -############################################################## -class PieceWiseLinFunc: - """ A class representing a piece-wise linear function. """ - - def __init__(self, x, y1, y2): - """ Constructs the piece-wise linear function. - - :param x: array of length N+1 defining the edges of the intervals of - the pwc function. - :param y1: array of length N defining the function values at the left - of the intervals. - :param y2: array of length N defining the function values at the right - of the intervals. - """ - # convert to array, which also ensures copying - self.x = np.array(x) - self.y1 = np.array(y1) - self.y2 = np.array(y2) - - def copy(self): - """ Returns a copy of itself - - :rtype: :class:`PieceWiseLinFunc` - """ - return PieceWiseLinFunc(self.x, self.y1, self.y2) - - def almost_equal(self, other, decimal=14): - """ Checks if the function is equal to another function up to `decimal` - precision. - - :param other: another :class:`PieceWiseLinFunc` - :returns: True if the two functions are equal up to `decimal` decimals, - False otherwise - :rtype: bool - """ - eps = 10.0**(-decimal) - return np.allclose(self.x, other.x, atol=eps, rtol=0.0) and \ - np.allclose(self.y1, other.y1, atol=eps, rtol=0.0) and \ - np.allclose(self.y2, other.y2, atol=eps, rtol=0.0) - - def get_plottable_data(self): - """ Returns two arrays containing x- and y-coordinates for immeditate - plotting of the piece-wise function. - - :returns: (x_plot, y_plot) containing plottable data - :rtype: pair of np.array - - Example:: - - x, y = f.get_plottable_data() - plt.plot(x, y, '-o', label="Piece-wise const function") - """ - x_plot = np.empty(2*len(self.x)-2) - x_plot[0] = self.x[0] - x_plot[1::2] = self.x[1:] - x_plot[2::2] = self.x[1:-1] - y_plot = np.empty_like(x_plot) - y_plot[0::2] = self.y1 - y_plot[1::2] = self.y2 - return x_plot, y_plot - - def integral(self, interval=None): - """ Returns the integral over the given interval. - - :param interval: integration interval given as a pair of floats, if - None the integral over the whole function is computed. - :type interval: Pair of floats or None. - :returns: the integral - :rtype: float - """ - - def intermediate_value(x0, x1, y0, y1, x): - """ computes the intermediate value of a linear function """ - return y0 + (y1-y0)*(x-x0)/(x1-x0) - - if interval is None: - # no interval given, integrate over the whole spike train - integral = np.sum((self.x[1:]-self.x[:-1]) * 0.5*(self.y1+self.y2)) - else: - # find the indices corresponding to the interval - start_ind = np.searchsorted(self.x, interval[0], side='right') - end_ind = np.searchsorted(self.x, interval[1], side='left')-1 - assert start_ind > 0 and end_ind < len(self.x), \ - "Invalid averaging interval" - # first the contribution from between the indices - integral = np.sum((self.x[start_ind+1:end_ind+1] - - self.x[start_ind:end_ind]) * - 0.5*(self.y1[start_ind:end_ind] + - self.y2[start_ind:end_ind])) - # correction from start to first index - integral += (self.x[start_ind]-interval[0]) * 0.5 * \ - (self.y2[start_ind-1] + - intermediate_value(self.x[start_ind-1], - self.x[start_ind], - self.y1[start_ind-1], - self.y2[start_ind-1], - interval[0] - )) - # correction from last index to end - integral += (interval[1]-self.x[end_ind]) * 0.5 * \ - (self.y1[end_ind] + - intermediate_value(self.x[end_ind], self.x[end_ind+1], - self.y1[end_ind], self.y2[end_ind], - interval[1] - )) - return integral - - def avrg(self, interval=None): - """ Computes the average of the piece-wise linear function: - :math:`a = 1/T int_0^T f(x) dx` where T is the length of the interval. - - :param interval: averaging interval given as a pair of floats, a - sequence of pairs for averaging multiple intervals, or - None, if None the average over the whole function is - computed. - :type interval: Pair, sequence of pairs, or None. - :returns: the average a. - :rtype: float - - """ - - if interval is None: - # no interval given, average over the whole spike train - return self.integral() / (self.x[-1]-self.x[0]) - - # 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): - # just one interval - a = self.integral(interval) / (interval[1]-interval[0]) - else: - # several intervals - a = 0.0 - int_length = 0.0 - for ival in interval: - a += self.integral(ival) - int_length += ival[1] - ival[0] - a /= int_length - return a - - def add(self, f): - """ Adds another PieceWiseLin function to this function. - Note: only functions defined on the same interval can be summed. - - :param f: :class:`PieceWiseLinFunc` function to be added. - :rtype: None - """ - assert self.x[0] == f.x[0], "The functions have different intervals" - assert self.x[-1] == f.x[-1], "The functions have different intervals" - - # python implementation - # 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_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.") - # use python backend - from 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( - self.x, self.y1, self.y2, f.x, f.y1, f.y2) - - def mul_scalar(self, fac): - """ Multiplies the function with a scalar value - - :param fac: Value to multiply - :type fac: double - :rtype: None - """ - self.y1 *= fac - self.y2 *= fac - - -############################################################## -# DiscreteFunction -############################################################## -class DiscreteFunction(object): - """ A class representing values defined on a discrete set of points. - """ - - def __init__(self, x, y, multiplicity): - """ Constructs the discrete function. - - :param x: array of length N defining the points at which the values are - defined. - :param y: array of length N degining the values at the points x. - :param multiplicity: array of length N defining the multiplicity of the - values. - """ - # convert parameters to arrays, also ensures copying - self.x = np.array(x) - self.y = np.array(y) - self.mp = np.array(multiplicity) - - def copy(self): - """ Returns a copy of itself - - :rtype: :class:`DiscreteFunction` - """ - return DiscreteFunction(self.x, self.y, self.mp) - - def almost_equal(self, other, decimal=14): - """ Checks if the function is equal to another function up to `decimal` - precision. - - :param other: another :class:`DiscreteFunction` - :returns: True if the two functions are equal up to `decimal` decimals, - False otherwise - :rtype: bool - """ - eps = 10.0**(-decimal) - return np.allclose(self.x, other.x, atol=eps, rtol=0.0) and \ - np.allclose(self.y, other.y, atol=eps, rtol=0.0) and \ - np.allclose(self.mp, other.mp, atol=eps, rtol=0.0) - - def get_plottable_data(self, averaging_window_size=0): - """ Returns two arrays containing x- and y-coordinates for plotting - the interval sequence. The optional parameter `averaging_window_size` - determines the size of an averaging window to smoothen the profile. If - this value is 0, no averaging is performed. - - :param averaging_window_size: size of the averaging window, default=0. - :returns: (x_plot, y_plot) containing plottable data - :rtype: pair of np.array - - Example:: - - x, y = f.get_plottable_data() - plt.plot(x, y, '-o', label="Discrete function") - """ - - if averaging_window_size > 0: - # for the averaged profile we have to take the multiplicity into - # account. values with higher multiplicity should be consider as if - # they appeared several times. Hence we can not know how many - # entries we have to consider to the left and right. Rather, we - # will iterate until some wanted multiplicity is reached. - - # the first value in self.mp contains the number of averaged - # profiles without any possible extra multiplicities - # (by implementation) - 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)): - - if self.mp[i] >= expected_mp: - # the current value contains already all the wanted - # multiplicity - y_plot[i] = self.y[i]/self.mp[i] - continue - - # first look to the right - y = self.y[i] - mp_r = self.mp[i] - j = i+1 - while j < len(y_plot): - if mp_r+self.mp[j] < expected_mp: - # if we still dont reach the required multiplicity - # we take the whole value - y += self.y[j] - mp_r += self.mp[j] - else: - # otherwise, just some fraction - y += self.y[j] * (expected_mp - mp_r)/self.mp[j] - mp_r += (expected_mp - mp_r) - break - j += 1 - - # same story to the left - mp_l = self.mp[i] - j = i-1 - while j >= 0: - if mp_l+self.mp[j] < expected_mp: - y += self.y[j] - mp_l += self.mp[j] - else: - y += self.y[j] * (expected_mp - mp_l)/self.mp[j] - mp_l += (expected_mp - mp_l) - break - j -= 1 - y_plot[i] = y/(mp_l+mp_r-self.mp[i]) - return 1.0*self.x, y_plot - - else: # k = 0 - - return 1.0*self.x, 1.0*self.y/self.mp - - def integral(self, interval=None): - """ Returns the integral over the given interval. For the discrete - function, this amounts to the sum over all values divided by the total - multiplicity. - - :param interval: integration interval given as a pair of floats, or a - sequence of pairs in case of multiple intervals, if - None the integral over the whole function is computed. - :type interval: Pair, sequence of pairs, or None. - :returns: the integral - :rtype: float - """ - - def get_indices(ival): - """ Retuns the indeces surrounding the given interval""" - start_ind = np.searchsorted(self.x, ival[0], side='right') - end_ind = np.searchsorted(self.x, ival[1], side='left') - assert start_ind > 0 and end_ind < len(self.x), \ - "Invalid averaging interval" - return start_ind, end_ind - - 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])) - else: - value = 0.0 - multiplicity = 0.0 - 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]) - return value/multiplicity - - def avrg(self, interval=None): - """ Computes the average of the interval sequence: - :math:`a = 1/N sum f_n ` where N is the number of intervals. - - :param interval: averaging interval given as a pair of floats, a - sequence of pairs for averaging multiple intervals, or - None, if None the average over the whole function is - computed. - :type interval: Pair, sequence of pairs, or None. - :returns: the average a. - :rtype: float - """ - return self.integral(interval) - - def add(self, f): - """ Adds another `DiscreteFunction` function to this function. - Note: only functions defined on the same interval can be summed. - - :param f: :class:`DiscreteFunction` function to be added. - :rtype: None - """ - assert self.x[0] == f.x[0], "The functions have different intervals" - assert self.x[-1] == f.x[-1], "The functions have different intervals" - - # cython version - try: - from cython_add import add_discrete_function_cython as \ - add_discrete_function_impl - except ImportError: - 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 python_backend import add_discrete_function_python as \ - add_discrete_function_impl - - self.x, self.y, self.mp = \ - add_discrete_function_impl(self.x, self.y, self.mp, - f.x, f.y, f.mp) - - def mul_scalar(self, fac): - """ Multiplies the function with a scalar value - - :param fac: Value to multiply - :type fac: double - :rtype: None - """ - self.y *= fac - - -def average_profile(profiles): - """ Computes the average profile from the given ISI- or SPIKE-profiles. - - :param profiles: list of :class:`PieceWiseConstFunc` or - :class:`PieceWiseLinFunc` representing ISI- or - SPIKE-profiles to be averaged. - :returns: the averages profile :math:`<S_{isi}>` or :math:`<S_{spike}>`. - :rtype: :class:`PieceWiseConstFunc` or :class:`PieceWiseLinFunc` - """ - assert len(profiles) > 1 - - avrg_profile = profiles[0].copy() - for i in xrange(1, len(profiles)): - avrg_profile.add(profiles[i]) - avrg_profile.mul_scalar(1.0/len(profiles)) # normalize - - return avrg_profile diff --git a/pyspike/generic.py b/pyspike/generic.py new file mode 100644 index 0000000..4f278d2 --- /dev/null +++ b/pyspike/generic.py @@ -0,0 +1,83 @@ +""" + +Generic functions to compute multi-variate profiles and distance matrices. + +Copyright 2015, Mario Mulansky <mario.mulansky@gmx.net> + +Distributed under the BSD License +""" + + +import numpy as np + + +############################################################ +# _generic_profile_multi +############################################################ +def _generic_profile_multi(spike_trains, pair_distance_func, indices=None): + """ Internal implementation detail, don't call this function directly, + use isi_profile_multi or spike_profile_multi instead. + + Computes the multi-variate distance for a set of spike-trains using the + pair_dist_func to compute pair-wise distances. That is it computes the + average distance of all pairs of spike-trains: + :math:`S(t) = 2/((N(N-1)) sum_{<i,j>} S_{i,j}`, + where the sum goes over all pairs <i,j>. + Args: + - spike_trains: list of spike trains + - pair_distance_func: function computing the distance of two spike trains + - indices: list of indices defining which spike trains to use, + if None all given spike trains are used (default=None) + Returns: + - The averaged multi-variate distance of all pairs + """ + if indices is None: + indices = np.arange(len(spike_trains)) + indices = np.array(indices) + # 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:]] + # start with first pair + (i, j) = pairs[0] + average_dist = pair_distance_func(spike_trains[i], spike_trains[j]) + for (i, j) in pairs[1:]: + current_dist = pair_distance_func(spike_trains[i], spike_trains[j]) + average_dist.add(current_dist) # add to the average + return average_dist, len(pairs) + + +############################################################ +# generic_distance_matrix +############################################################ +def _generic_distance_matrix(spike_trains, dist_function, + indices=None, interval=None): + """ Internal implementation detail. Don't use this function directly. + Instead use isi_distance_matrix or spike_distance_matrix. + Computes the time averaged distance of all pairs of spike-trains. + Args: + - spike_trains: list of spike trains + - indices: list of indices defining which spike-trains to use + if None all given spike-trains are used (default=None) + Return: + - a 2D array of size len(indices)*len(indices) containing the average + pair-wise distance + """ + 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 = dist_function(spike_trains[i], spike_trains[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 new file mode 100644 index 0000000..c2ef8e8 --- /dev/null +++ b/pyspike/isi_distance.py @@ -0,0 +1,134 @@ +""" + +Module containing several functions to compute the ISI profiles and distances + +Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net> + +Distributed under the BSD License +""" + +from pyspike import PieceWiseConstFunc +from pyspike.generic import _generic_profile_multi, _generic_distance_matrix + + +############################################################ +# isi_profile +############################################################ +def isi_profile(spikes1, spikes2): + """ Computes the isi-distance profile :math:`S_{isi}(t)` of the two given + spike trains. Retruns the profile as a PieceWiseConstFunc object. The S_isi + values are defined positive S_isi(t)>=0. The spike trains are expected + to have auxiliary spikes at the beginning and end of the interval. Use the + function add_auxiliary_spikes to add those spikes to the spike train. + + :param spikes1: ordered array of spike times with auxiliary spikes. + :param spikes2: ordered array of spike times with auxiliary spikes. + :returns: The isi-distance profile :math:`S_{isi}(t)` + :rtype: :class:`pyspike.function.PieceWiseConstFunc` + + """ + # check for auxiliary spikes - first and last spikes should be identical + assert spikes1[0] == spikes2[0], \ + "Given spike trains seems not to have auxiliary spikes!" + assert spikes1[-1] == spikes2[-1], \ + "Given spike trains seems not to have auxiliary spikes!" + + # load cython implementation + try: + from cython.cython_distance import isi_distance_cython \ + as isi_distance_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 \ +Falling back to slow python backend.") + # use python backend + from cython.python_backend import isi_distance_python \ + as isi_distance_impl + + times, values = isi_distance_impl(spikes1, spikes2) + return PieceWiseConstFunc(times, values) + + +############################################################ +# isi_distance +############################################################ +def isi_distance(spikes1, spikes2, interval=None): + """ Computes the isi-distance I of the given spike trains. The + isi-distance is the integral over the isi distance profile + :math:`S_{isi}(t)`: + + .. math:: I = \int_{T_0}^{T_1} S_{isi}(t) dt. + + :param spikes1: ordered array of spike times with auxiliary spikes. + :param spikes2: ordered array of spike times with auxiliary spikes. + :param interval: averaging interval given as a pair of floats (T0, T1), + if None the average over the whole function is computed. + :type interval: Pair of floats or None. + :returns: The isi-distance I. + :rtype: double + """ + return isi_profile(spikes1, spikes2).avrg(interval) + + +############################################################ +# isi_profile_multi +############################################################ +def isi_profile_multi(spike_trains, indices=None): + """ computes the multi-variate isi distance profile for a set of spike + trains. That is the average isi-distance of all pairs of spike-trains: + S_isi(t) = 2/((N(N-1)) sum_{<i,j>} S_{isi}^{i,j}, + where the sum goes over all pairs <i,j> + + :param spike_trains: list of spike trains + :param indices: list of indices defining which spike trains to use, + if None all given spike trains are used (default=None) + :type state: list or None + :returns: The averaged isi profile :math:`<S_{isi}>(t)` + :rtype: :class:`pyspike.function.PieceWiseConstFunc` + """ + average_dist, M = _generic_profile_multi(spike_trains, isi_profile, + indices) + average_dist.mul_scalar(1.0/M) # normalize + return average_dist + + +############################################################ +# isi_distance_multi +############################################################ +def isi_distance_multi(spike_trains, indices=None, interval=None): + """ computes the multi-variate isi-distance for a set of spike-trains. + That is the time average of the multi-variate spike profile: + I = \int_0^T 2/((N(N-1)) sum_{<i,j>} S_{isi}^{i,j}, + where the sum goes over all pairs <i,j> + + :param spike_trains: list of spike trains + :param indices: list of indices defining which spike trains to use, + if None all given spike trains are used (default=None) + :param interval: averaging interval given as a pair of floats, if None + the average over the whole function is computed. + :type interval: Pair of floats or None. + :returns: The time-averaged isi distance :math:`I` + :rtype: double + """ + return isi_profile_multi(spike_trains, indices).avrg(interval) + + +############################################################ +# isi_distance_matrix +############################################################ +def isi_distance_matrix(spike_trains, indices=None, interval=None): + """ Computes the time averaged isi-distance of all pairs of spike-trains. + + :param spike_trains: list of spike trains + :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 interval: averaging interval given as a pair of floats, if None + the average over the whole function is computed. + :type interval: Pair of floats or None. + :returns: 2D array with the pair wise time average isi distances + :math:`I_{ij}` + :rtype: np.array + """ + return _generic_distance_matrix(spike_trains, isi_distance, + indices, interval) diff --git a/pyspike/spike_distance.py b/pyspike/spike_distance.py new file mode 100644 index 0000000..f721c86 --- /dev/null +++ b/pyspike/spike_distance.py @@ -0,0 +1,136 @@ +""" + +Module containing several functions to compute SPIKE profiles and distances + +Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net> + +Distributed under the BSD License +""" + +from pyspike import PieceWiseLinFunc +from pyspike.generic import _generic_profile_multi, _generic_distance_matrix + + +############################################################ +# spike_profile +############################################################ +def spike_profile(spikes1, spikes2): + """ Computes the spike-distance profile S_spike(t) of the two given spike + trains. Returns the profile as a PieceWiseLinFunc object. The S_spike + values are defined positive S_spike(t)>=0. The spike trains are expected to + have auxiliary spikes at the beginning and end of the interval. Use the + function add_auxiliary_spikes to add those spikes to the spike train. + + :param spikes1: ordered array of spike times with auxiliary spikes. + :param spikes2: ordered array of spike times with auxiliary spikes. + :returns: The spike-distance profile :math:`S_{spike}(t)`. + :rtype: :class:`pyspike.function.PieceWiseLinFunc` + + """ + # check for auxiliary spikes - first and last spikes should be identical + assert spikes1[0] == spikes2[0], \ + "Given spike trains seems not to have auxiliary spikes!" + assert spikes1[-1] == spikes2[-1], \ + "Given spike trains seems not to have auxiliary spikes!" + + # cython implementation + try: + from cython.cython_distance import spike_distance_cython \ + as spike_distance_impl + except ImportError: + 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 spike_distance_python \ + as spike_distance_impl + + times, y_starts, y_ends = spike_distance_impl(spikes1, spikes2) + return PieceWiseLinFunc(times, y_starts, y_ends) + + +############################################################ +# spike_distance +############################################################ +def spike_distance(spikes1, spikes2, interval=None): + """ Computes the spike-distance S of the given spike trains. The + spike-distance is the integral over the isi distance profile S_spike(t): + + .. math:: S = \int_{T_0}^{T_1} S_{spike}(t) dt. + + :param spikes1: ordered array of spike times with auxiliary spikes. + :param spikes2: ordered array of spike times with auxiliary spikes. + :param interval: averaging interval given as a pair of floats (T0, T1), + if None the average over the whole function is computed. + :type interval: Pair of floats or None. + :returns: The spike-distance. + :rtype: double + + """ + return spike_profile(spikes1, spikes2).avrg(interval) + + +############################################################ +# spike_profile_multi +############################################################ +def spike_profile_multi(spike_trains, indices=None): + """ Computes the multi-variate spike distance profile for a set of spike + trains. That is the average spike-distance of all pairs of spike-trains: + :math:`S_spike(t) = 2/((N(N-1)) sum_{<i,j>} S_{spike}^{i, j}`, + where the sum goes over all pairs <i,j> + + :param spike_trains: list of spike trains + :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 + :returns: The averaged spike profile :math:`<S_{spike}>(t)` + :rtype: :class:`pyspike.function.PieceWiseLinFunc` + + """ + average_dist, M = _generic_profile_multi(spike_trains, spike_profile, + indices) + average_dist.mul_scalar(1.0/M) # normalize + return average_dist + + +############################################################ +# spike_distance_multi +############################################################ +def spike_distance_multi(spike_trains, indices=None, interval=None): + """ Computes the multi-variate spike distance for a set of spike trains. + That is the time average of the multi-variate spike profile: + S_{spike} = \int_0^T 2/((N(N-1)) sum_{<i,j>} S_{spike}^{i, j} dt + where the sum goes over all pairs <i,j> + + :param spike_trains: list of spike trains + :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 interval: averaging interval given as a pair of floats, if None + the average over the whole function is computed. + :type interval: Pair of floats or None. + :returns: The averaged spike distance S. + :rtype: double + """ + return spike_profile_multi(spike_trains, indices).avrg(interval) + + +############################################################ +# spike_distance_matrix +############################################################ +def spike_distance_matrix(spike_trains, indices=None, interval=None): + """ Computes the time averaged spike-distance of all pairs of spike-trains. + + :param spike_trains: list of spike trains + :param 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 interval: averaging interval given as a pair of floats, if None + the average over the whole function is computed. + :type interval: Pair of floats or None. + :returns: 2D array with the pair wise time average spike distances + :math:`S_{ij}` + :rtype: np.array + """ + return _generic_distance_matrix(spike_trains, spike_distance, + indices, interval) diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py new file mode 100644 index 0000000..342bf69 --- /dev/null +++ b/pyspike/spike_sync.py @@ -0,0 +1,136 @@ +""" + +Module containing several functions to compute SPIKE-Synchronization profiles +and distances + +Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net> + +Distributed under the BSD License +""" + +from functools import partial +from pyspike import DiscreteFunc +from pyspike.generic import _generic_profile_multi, _generic_distance_matrix + + +############################################################ +# spike_sync_profile +############################################################ +def spike_sync_profile(spikes1, spikes2): + """ Computes the spike-synchronization profile S_sync(t) of the two given + spike trains. Returns the profile as a DiscreteFunction object. The S_sync + values are either 1 or 0, indicating the presence or absence of a + coincidence. The spike trains are expected to have auxiliary spikes at the + beginning and end of the interval. Use the function add_auxiliary_spikes to + add those spikes to the spike train. + + :param spikes1: ordered array of spike times with auxiliary spikes. + :param spikes2: ordered array of spike times with auxiliary spikes. + :returns: The spike-distance profile :math:`S_{sync}(t)`. + :rtype: :class:`pyspike.function.DiscreteFunction` + + """ + + # cython implementation + try: + from cython.cython_distance import coincidence_cython \ + as coincidence_impl + except ImportError: + 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_impl + + times, coincidences, multiplicity = coincidence_impl(spikes1, spikes2) + + return DiscreteFunc(times, coincidences, multiplicity) + + +############################################################ +# spike_sync +############################################################ +def spike_sync(spikes1, spikes2, interval=None): + """ Computes the spike synchronization value SYNC of the given spike + trains. The spike synchronization value is the computed as the total number + of coincidences divided by the total number of spikes: + + .. math:: SYNC = \sum_n C_n / N. + + :param spikes1: ordered array of spike times with auxiliary spikes. + :param spikes2: ordered array of spike times with auxiliary spikes. + :param interval: averaging interval given as a pair of floats (T0, T1), + if None the average over the whole function is computed. + :type interval: Pair of floats or None. + :returns: The spike synchronization value. + :rtype: double + """ + return spike_sync_profile(spikes1, spikes2).avrg(interval) + + +############################################################ +# spike_sync_profile_multi +############################################################ +def spike_sync_profile_multi(spike_trains, indices=None): + """ Computes the multi-variate spike synchronization profile for a set of + spike trains. For each spike in the set of spike trains, the multi-variate + profile is defined as the number of coincidences 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 spike trains + :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 + :returns: The multi-variate spike sync profile :math:`<S_{sync}>(t)` + :rtype: :class:`pyspike.function.DiscreteFunction` + + """ + prof_func = partial(spike_sync_profile) + average_dist, M = _generic_profile_multi(spike_trains, prof_func, + indices) + # average_dist.mul_scalar(1.0/M) # no normalization here! + return average_dist + + +############################################################ +# spike_sync_multi +############################################################ +def spike_sync_multi(spike_trains, indices=None, interval=None): + """ Computes the multi-variate spike synchronization value for a set of + spike trains. + + :param spike_trains: list of spike trains + :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 interval: averaging interval given as a pair of floats, if None + the average over the whole function is computed. + :type interval: Pair of floats or None. + :returns: The multi-variate spike synchronization value SYNC. + :rtype: double + """ + return spike_sync_profile_multi(spike_trains, indices).avrg(interval) + + +############################################################ +# spike_sync_matrix +############################################################ +def spike_sync_matrix(spike_trains, indices=None, interval=None): + """ Computes the overall spike-synchronization value of all pairs of + spike-trains. + + :param spike_trains: list of spike trains + :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 interval: averaging interval given as a pair of floats, if None + the average over the whole function is computed. + :type interval: Pair of floats or None. + :returns: 2D array with the pair wise time spike synchronization values + :math:`SYNC_{ij}` + :rtype: np.array + """ + return _generic_distance_matrix(spike_trains, spike_sync, + indices, interval) @@ -21,8 +21,8 @@ except ImportError: else: use_cython = True -if os.path.isfile("pyspike/cython_add.c") and \ - os.path.isfile("pyspike/cython_distance.c"): +if os.path.isfile("pyspike/cython/cython_add.c") and \ + os.path.isfile("pyspike/cython/cython_distance.c"): use_c = True else: use_c = False @@ -32,14 +32,14 @@ ext_modules = [] if use_cython: # Cython is available, compile .pyx -> .c ext_modules += [ - Extension("pyspike.cython_add", ["pyspike/cython_add.pyx"]), - Extension("pyspike.cython_distance", ["pyspike/cython_distance.pyx"]), + Extension("pyspike.cython.cython_add", ["pyspike/cython/cython_add.pyx"]), + Extension("pyspike.cython.cython_distance", ["pyspike/cython/cython_distance.pyx"]), ] cmdclass.update({'build_ext': build_ext}) elif use_c: # c files are there, compile to binaries ext_modules += [ - Extension("pyspike.cython_add", ["pyspike/cython_add.c"]), - Extension("pyspike.cython_distance", ["pyspike/cython_distance.c"]), + Extension("pyspike.cython.cython_add", ["pyspike/cython/cython_add.c"]), + Extension("pyspike.cython.cython_distance", ["pyspike/cython/cython_distance.c"]), ] # neither cython nor c files available -> automatic fall-back to python backend @@ -78,7 +78,7 @@ train similarity', 'Programming Language :: Python :: 2.7', ], package_data={ - 'pyspike': ['cython_add.c', 'cython_distance.c'], + 'pyspike': ['cython/cython_add.c', 'cython/cython_distance.c'], 'test': ['Spike_testdata.txt'] } ) diff --git a/test/test_function.py b/test/test_function.py index 933fd2e..d81b03a 100644 --- a/test/test_function.py +++ b/test/test_function.py @@ -92,11 +92,12 @@ def test_pwc_avrg(): y = [0.5, 1.0, -0.25, 0.0, 1.5] f2 = spk.PieceWiseConstFunc(x, y) - f_avrg = spk.average_profile([f1, f2]) + f1.add(f2) + f1.mul_scalar(0.5) x_expected = [0.0, 0.75, 1.0, 2.0, 2.5, 2.7, 4.0] y_expected = [0.75, 1.0, 0.25, 0.625, 0.375, 1.125] - assert_array_almost_equal(f_avrg.x, x_expected, decimal=16) - assert_array_almost_equal(f_avrg.y, y_expected, decimal=16) + assert_array_almost_equal(f1.x, x_expected, decimal=16) + assert_array_almost_equal(f1.y, y_expected, decimal=16) def test_pwl(): @@ -196,11 +197,12 @@ def test_pwl_avrg(): y2_expected = np.array([0.8+1.0+0.5*0.75, 1.5+1.0-0.8*0.25/1.25, -0.4+0.2, 1.5-1.0, 0.75-0.5*0.2/1.5, 2.25]) / 2 - f_avrg = spk.average_profile([f1, f2]) + f1.add(f2) + f1.mul_scalar(0.5) - assert_array_almost_equal(f_avrg.x, x_expected, decimal=16) - assert_array_almost_equal(f_avrg.y1, y1_expected, decimal=16) - assert_array_almost_equal(f_avrg.y2, y2_expected, decimal=16) + assert_array_almost_equal(f1.x, x_expected, decimal=16) + assert_array_almost_equal(f1.y1, y1_expected, decimal=16) + assert_array_almost_equal(f1.y2, y2_expected, decimal=16) def test_df(): @@ -208,7 +210,7 @@ def test_df(): x = [0.0, 1.0, 2.0, 2.5, 4.0] y = [0.0, 1.0, 1.0, 0.0, 1.0] mp = [1.0, 2.0, 1.0, 2.0, 1.0] - f = spk.DiscreteFunction(x, y, mp) + f = spk.DiscreteFunc(x, y, mp) xp, yp = f.get_plottable_data() xp_expected = [0.0, 1.0, 2.0, 2.5, 4.0] @@ -237,6 +239,9 @@ if __name__ == "__main__": test_pwc() test_pwc_add() test_pwc_mul() + test_pwc_avrg() test_pwl() test_pwl_add() test_pwl_mul() + test_pwl_avrg() + test_df() |