From ae90eb6bdd6afd716cb90a4aeeec8e6e77635d10 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Wed, 28 Jan 2015 17:42:19 +0100 Subject: each function class in separate source file --- pyspike/__init__.py | 9 +- pyspike/distances.py | 13 +- pyspike/function.py | 585 -------------------------------------------------- test/test_function.py | 21 +- 4 files changed, 27 insertions(+), 601 deletions(-) delete mode 100644 pyspike/function.py diff --git a/pyspike/__init__.py b/pyspike/__init__.py index 55687e6..f480964 100644 --- a/pyspike/__init__.py +++ b/pyspike/__init__.py @@ -4,10 +4,13 @@ Copyright 2014, Mario Mulansky Distributed under the BSD License """ -__all__ = ["function", "distances", "spikes"] +__all__ = ["distances", "spikes", "PieceWiseConstFunc", "PieceWiseLinFunc", + "DiscreteFunc"] + +from PieceWiseConstFunc import PieceWiseConstFunc +from PieceWiseLinFunc import PieceWiseLinFunc +from DiscreteFunc import 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, \ diff --git a/pyspike/distances.py b/pyspike/distances.py index 5476b6f..9077871 100644 --- a/pyspike/distances.py +++ b/pyspike/distances.py @@ -11,7 +11,7 @@ import numpy as np import threading from functools import partial -from pyspike import PieceWiseConstFunc, PieceWiseLinFunc, DiscreteFunction +from pyspike import PieceWiseConstFunc, PieceWiseLinFunc, DiscreteFunc ############################################################ @@ -161,7 +161,7 @@ Falling back to slow python backend.") times, coincidences, multiplicity = coincidence_impl(spikes1, spikes2) - return DiscreteFunction(times, coincidences, multiplicity) + return DiscreteFunc(times, coincidences, multiplicity) ############################################################ @@ -212,7 +212,8 @@ def _generic_profile_multi(spike_trains, pair_distance_func, indices=None): 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 = [(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]) @@ -251,7 +252,8 @@ def _multi_distance_par(spike_trains, pair_distance_func, indices=None): 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 = [(indices[i], j) for i in range(len(indices)) + for j in indices[i+1:]] num_pairs = len(pairs) # start with first pair @@ -430,7 +432,8 @@ def _generic_distance_matrix(spike_trains, dist_function, assert (indices < len(spike_trains)).all() and (indices >= 0).all(), \ "Invalid index list." # generate a list of possible index pairs - pairs = [(indices[i], j) for i in range(len(indices)) for j in indices[i+1:]] + 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: 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 - -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:`` or :math:``. - :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/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() -- cgit v1.2.3 From 01fa64bd8a8931544a1734de104e2fe72018694f Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Wed, 28 Jan 2015 19:22:02 +0100 Subject: added missing py files --- pyspike/DiscreteFunc.py | 242 ++++++++++++++++++++++++++++++++++++++++++ pyspike/PieceWiseConstFunc.py | 168 +++++++++++++++++++++++++++++ pyspike/PieceWiseLinFunc.py | 197 ++++++++++++++++++++++++++++++++++ 3 files changed, 607 insertions(+) create mode 100644 pyspike/DiscreteFunc.py create mode 100644 pyspike/PieceWiseConstFunc.py create mode 100644 pyspike/PieceWiseLinFunc.py diff --git a/pyspike/DiscreteFunc.py b/pyspike/DiscreteFunc.py new file mode 100644 index 0000000..3e24284 --- /dev/null +++ b/pyspike/DiscreteFunc.py @@ -0,0 +1,242 @@ +""" +Class representing discrete functions. + +Copyright 2014-2015, Mario Mulansky + +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_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:`` or :math:``. + :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..e639dfc --- /dev/null +++ b/pyspike/PieceWiseConstFunc.py @@ -0,0 +1,168 @@ +""" +Class representing piece-wise constant functions. + +Copyright 2014-2015, Mario Mulansky + +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 diff --git a/pyspike/PieceWiseLinFunc.py b/pyspike/PieceWiseLinFunc.py new file mode 100644 index 0000000..58a20e5 --- /dev/null +++ b/pyspike/PieceWiseLinFunc.py @@ -0,0 +1,197 @@ +""" +Class representing piece-wise linear functions. + +Copyright 2014-2015, Mario Mulansky + +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_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 -- cgit v1.2.3 From be9eeb2e48115134b93ec0cd9035d97117bd019e Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Tue, 3 Feb 2015 10:45:32 +0100 Subject: split distance.py into 3 separate modules --- pyspike/__init__.py | 14 +- pyspike/distances.py | 507 ---------------------------------------------- pyspike/generic.py | 83 ++++++++ pyspike/isi_distance.py | 132 ++++++++++++ pyspike/spike_distance.py | 135 ++++++++++++ pyspike/spike_sync.py | 136 +++++++++++++ 6 files changed, 494 insertions(+), 513 deletions(-) delete mode 100644 pyspike/distances.py create mode 100644 pyspike/generic.py create mode 100644 pyspike/isi_distance.py create mode 100644 pyspike/spike_distance.py create mode 100644 pyspike/spike_sync.py diff --git a/pyspike/__init__.py b/pyspike/__init__.py index f480964..1c2efbc 100644 --- a/pyspike/__init__.py +++ b/pyspike/__init__.py @@ -4,18 +4,20 @@ Copyright 2014, Mario Mulansky Distributed under the BSD License """ -__all__ = ["distances", "spikes", "PieceWiseConstFunc", "PieceWiseLinFunc", +__all__ = ["isi_distance", "spike_distance", "spike_sync", + "spikes", "PieceWiseConstFunc", "PieceWiseLinFunc", "DiscreteFunc"] from PieceWiseConstFunc import PieceWiseConstFunc from PieceWiseLinFunc import PieceWiseLinFunc from DiscreteFunc import DiscreteFunc -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 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/distances.py b/pyspike/distances.py deleted file mode 100644 index 9077871..0000000 --- a/pyspike/distances.py +++ /dev/null @@ -1,507 +0,0 @@ -""" distances.py - -Module containing several functions to compute spike distances - -Copyright 2014, Mario Mulansky - -Distributed under the BSD License -""" - -import numpy as np -import threading -from functools import partial - -from pyspike import PieceWiseConstFunc, PieceWiseLinFunc, DiscreteFunc - - -############################################################ -# 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 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) - - -############################################################ -# _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_{} S_{i,j}`, - where the sum goes over all pairs . - Args: - - spike_trains: list of spike trains - - pair_distance_func: function computing the distance of two spike trains - - indices: list of indices defining which spike trains to use, - if None all given spike trains are used (default=None) - Returns: - - The averaged multi-variate distance of all pairs - """ - if indices is None: - indices = np.arange(len(spike_trains)) - indices = np.array(indices) - # 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_{} S_{isi}^{i,j}, - where the sum goes over all pairs - - :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:`(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_{} S_{isi}^{i,j}, - where the sum goes over all pairs - - :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_{} S_{spike}^{i, j}`, - where the sum goes over all pairs - - :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:`(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_{} S_{spike}^{i, j} dt - where the sum goes over all pairs - - :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:`(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/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 + +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_{} S_{i,j}`, + where the sum goes over all pairs . + Args: + - spike_trains: list of spike trains + - pair_distance_func: function computing the distance of two spike trains + - indices: list of indices defining which spike trains to use, + if None all given spike trains are used (default=None) + Returns: + - The averaged multi-variate distance of all pairs + """ + if indices is None: + indices = np.arange(len(spike_trains)) + indices = np.array(indices) + # 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..745d280 --- /dev/null +++ b/pyspike/isi_distance.py @@ -0,0 +1,132 @@ +""" + +Module containing several functions to compute the ISI profiles and distances + +Copyright 2014-2015, Mario Mulansky + +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_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) + + +############################################################ +# 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_{} S_{isi}^{i,j}, + where the sum goes over all pairs + + :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:`(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_{} S_{isi}^{i,j}, + where the sum goes over all pairs + + :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..2c989a4 --- /dev/null +++ b/pyspike/spike_distance.py @@ -0,0 +1,135 @@ +""" + +Module containing several functions to compute SPIKE profiles and distances + +Copyright 2014-2015, Mario Mulansky + +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_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_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_{} S_{spike}^{i, j}`, + where the sum goes over all pairs + + :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:`(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_{} S_{spike}^{i, j} dt + where the sum goes over all pairs + + :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..f7cbe1d --- /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 + +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_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 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:`(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) + + +############################################################ +# 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) -- cgit v1.2.3 From 7f20d9a8076326c1800373a7f95f4871873f14b0 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Tue, 3 Feb 2015 10:46:51 +0100 Subject: copyright, docs --- pyspike/__init__.py | 2 +- pyspike/spike_sync.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/pyspike/__init__.py b/pyspike/__init__.py index 1c2efbc..945dd4e 100644 --- a/pyspike/__init__.py +++ b/pyspike/__init__.py @@ -1,5 +1,5 @@ """ -Copyright 2014, Mario Mulansky +Copyright 2014-2015, Mario Mulansky Distributed under the BSD License """ diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py index f7cbe1d..bded8da 100644 --- a/pyspike/spike_sync.py +++ b/pyspike/spike_sync.py @@ -95,7 +95,7 @@ def spike_sync_profile_multi(spike_trains, indices=None): ############################################################ -# spike_distance_multi +# spike_sync_multi ############################################################ def spike_sync_multi(spike_trains, indices=None, interval=None): """ Computes the multi-variate spike synchronization value for a set of -- cgit v1.2.3 From 6eb6bc486027d3d5304a94cfb417a2257f2b6fd9 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Tue, 3 Feb 2015 12:19:53 +0100 Subject: moved cython functions to subdirectory --- pyspike/DiscreteFunc.py | 4 +- pyspike/PieceWiseConstFunc.py | 4 +- pyspike/PieceWiseLinFunc.py | 4 +- pyspike/cython/cython_add.pyx | 235 ++++++++++++++++++ pyspike/cython/cython_distance.pyx | 312 ++++++++++++++++++++++++ pyspike/cython/python_backend.py | 485 +++++++++++++++++++++++++++++++++++++ pyspike/cython_add.pyx | 235 ------------------ pyspike/cython_distance.pyx | 312 ------------------------ pyspike/isi_distance.py | 6 +- pyspike/python_backend.py | 485 ------------------------------------- pyspike/spike_distance.py | 5 +- pyspike/spike_sync.py | 4 +- setup.py | 14 +- 13 files changed, 1054 insertions(+), 1051 deletions(-) create mode 100644 pyspike/cython/cython_add.pyx create mode 100644 pyspike/cython/cython_distance.pyx create mode 100644 pyspike/cython/python_backend.py delete mode 100644 pyspike/cython_add.pyx delete mode 100644 pyspike/cython_distance.pyx delete mode 100644 pyspike/python_backend.py diff --git a/pyspike/DiscreteFunc.py b/pyspike/DiscreteFunc.py index 3e24284..2283e03 100644 --- a/pyspike/DiscreteFunc.py +++ b/pyspike/DiscreteFunc.py @@ -198,7 +198,7 @@ class DiscreteFunc(object): # cython version try: - from 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 \ @@ -206,7 +206,7 @@ sure that PySpike is installed by running\n\ 'python setup.py build_ext --inplace'! \ \n Falling back to slow python backend.") # use python backend - from python_backend import add_discrete_function_python as \ + from cython.python_backend import add_discrete_function_python as \ add_discrete_function_impl self.x, self.y, self.mp = \ diff --git a/pyspike/PieceWiseConstFunc.py b/pyspike/PieceWiseConstFunc.py index e639dfc..dc57ab1 100644 --- a/pyspike/PieceWiseConstFunc.py +++ b/pyspike/PieceWiseConstFunc.py @@ -146,14 +146,14 @@ class PieceWiseConstFunc(object): # cython version try: - from 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'! \ \n Falling back to slow python backend.") # use python backend - from 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 58a20e5..bc0aa2a 100644 --- a/pyspike/PieceWiseLinFunc.py +++ b/pyspike/PieceWiseLinFunc.py @@ -173,14 +173,14 @@ class PieceWiseLinFunc: # cython version try: - from 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.") # use python backend - from 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/cython/cython_add.pyx b/pyspike/cython/cython_add.pyx new file mode 100644 index 0000000..ac64005 --- /dev/null +++ b/pyspike/cython/cython_add.pyx @@ -0,0 +1,235 @@ +#cython: boundscheck=False +#cython: wraparound=False +#cython: cdivision=True + +""" +cython_add.pyx + +cython implementation of the add function for piece-wise const and +piece-wise linear functions + +Note: using cython memoryviews (e.g. double[:]) instead of ndarray objects +improves the performance of spike_distance by a factor of 10! + +Copyright 2014, Mario Mulansky + +Distributed under the BSD License + +""" + +""" +To test whether things can be optimized: remove all yellow stuff +in the html output:: + + cython -a cython_add.pyx + +which gives:: + + cython_add.html + +""" + +import numpy as np +cimport numpy as np + +from libc.math cimport fabs + +DTYPE = np.float +ctypedef np.float_t DTYPE_t + + +############################################################ +# add_piece_wise_const_cython +############################################################ +def add_piece_wise_const_cython(double[:] x1, double[:] y1, + double[:] x2, double[:] y2): + + cdef int N1 = len(x1) + cdef int N2 = len(x2) + cdef double[:] x_new = np.empty(N1+N2) + cdef double[:] y_new = np.empty(N1+N2-1) + cdef int index1 = 0 + cdef int index2 = 0 + cdef int index = 0 + cdef int i + with nogil: # release the interpreter lock to allow multi-threading + x_new[0] = x1[0] + y_new[0] = y1[0] + y2[0] + while (index1+1 < N1-1) and (index2+1 < N2-1): + index += 1 + # print(index1+1, x1[index1+1], y1[index1+1], x_new[index]) + if x1[index1+1] < x2[index2+1]: + index1 += 1 + x_new[index] = x1[index1] + elif x1[index1+1] > x2[index2+1]: + index2 += 1 + x_new[index] = x2[index2] + else: # x1[index1+1] == x2[index2+1]: + index1 += 1 + index2 += 1 + x_new[index] = x1[index1] + y_new[index] = y1[index1] + y2[index2] + # one array reached the end -> copy the contents of the other to the end + if index1+1 < N1-1: + x_new[index+1:index+1+N1-index1-1] = x1[index1+1:] + for i in xrange(N1-index1-2): + y_new[index+1+i] = y1[index1+1+i] + y2[N2-2] + index += N1-index1-2 + elif index2+1 < N2-1: + x_new[index+1:index+1+N2-index2-1] = x2[index2+1:] + for i in xrange(N2-index2-2): + y_new[index+1+i] = y2[index2+1+i] + y1[N1-2] + index += N2-index2-2 + else: # both arrays reached the end simultaneously + # only the last x-value missing + x_new[index+1] = x1[N1-1] + # the last value is again the end of the interval + # x_new[index+1] = x1[-1] + # only use the data that was actually filled + x1 = x_new[:index+2] + y1 = y_new[:index+1] + # end nogil + return np.array(x_new[:index+2]), np.array(y_new[:index+1]) + + +############################################################ +# add_piece_wise_lin_cython +############################################################ +def add_piece_wise_lin_cython(double[:] x1, double[:] y11, double[:] y12, + double[:] x2, double[:] y21, double[:] y22): + cdef int N1 = len(x1) + cdef int N2 = len(x2) + cdef double[:] x_new = np.empty(N1+N2) + cdef double[:] y1_new = np.empty(N1+N2-1) + cdef double[:] y2_new = np.empty_like(y1_new) + cdef int index1 = 0 # index for self + cdef int index2 = 0 # index for f + cdef int index = 0 # index for new + cdef int i + cdef double y + with nogil: # release the interpreter lock to allow multi-threading + x_new[0] = x1[0] + y1_new[0] = y11[0] + y21[0] + while (index1+1 < N1-1) and (index2+1 < N2-1): + # print(index1+1, x1[index1+1], self.y[index1+1], x_new[index]) + if x1[index1+1] < x2[index2+1]: + # first compute the end value of the previous interval + # linear interpolation of the interval + y = y21[index2] + (y22[index2]-y21[index2]) * \ + (x1[index1+1]-x2[index2]) / (x2[index2+1]-x2[index2]) + y2_new[index] = y12[index1] + y + index1 += 1 + index += 1 + x_new[index] = x1[index1] + # and the starting value for the next interval + y1_new[index] = y11[index1] + y + elif x1[index1+1] > x2[index2+1]: + # first compute the end value of the previous interval + # linear interpolation of the interval + y = y11[index1] + (y12[index1]-y11[index1]) * \ + (x2[index2+1]-x1[index1]) / \ + (x1[index1+1]-x1[index1]) + y2_new[index] = y22[index2] + y + index2 += 1 + index += 1 + x_new[index] = x2[index2] + # and the starting value for the next interval + y1_new[index] = y21[index2] + y + else: # x1[index1+1] == x2[index2+1]: + y2_new[index] = y12[index1] + y22[index2] + index1 += 1 + index2 += 1 + index += 1 + x_new[index] = x1[index1] + y1_new[index] = y11[index1] + y21[index2] + # one array reached the end -> copy the contents of the other to the end + if index1+1 < N1-1: + x_new[index+1:index+1+N1-index1-1] = x1[index1+1:] + for i in xrange(N1-index1-2): + # compute the linear interpolations value + y = y21[index2] + (y22[index2]-y21[index2]) * \ + (x1[index1+1+i]-x2[index2]) / (x2[index2+1]-x2[index2]) + y1_new[index+1+i] = y11[index1+1+i] + y + y2_new[index+i] = y12[index1+i] + y + index += N1-index1-2 + elif index2+1 < N2-1: + x_new[index+1:index+1+N2-index2-1] = x2[index2+1:] + # compute the linear interpolations values + for i in xrange(N2-index2-2): + y = y11[index1] + (y12[index1]-y11[index1]) * \ + (x2[index2+1+i]-x1[index1]) / \ + (x1[index1+1]-x1[index1]) + y1_new[index+1+i] = y21[index2+1+i] + y + y2_new[index+i] = y22[index2+i] + y + index += N2-index2-2 + else: # both arrays reached the end simultaneously + # only the last x-value missing + x_new[index+1] = x1[N1-1] + # finally, the end value for the last interval + y2_new[index] = y12[N1-2]+y22[N2-2] + # only use the data that was actually filled + # end nogil + return (np.array(x_new[:index+2]), + np.array(y1_new[:index+1]), + np.array(y2_new[:index+1])) + + +############################################################ +# add_discrete_function_cython +############################################################ +def add_discrete_function_cython(double[:] x1, double[:] y1, double[:] mp1, + double[:] x2, double[:] y2, double[:] mp2): + + cdef double[:] x_new = np.empty(len(x1) + len(x2)) + cdef double[:] y_new = np.empty_like(x_new) + cdef double[:] mp_new = np.empty_like(x_new) + cdef int index1 = 0 + cdef int index2 = 0 + cdef int index = 0 + cdef int N1 = len(y1) + cdef int N2 = len(y2) + x_new[0] = x1[0] + while (index1+1 < N1) and (index2+1 < N2): + if x1[index1+1] < x2[index2+1]: + index1 += 1 + index += 1 + x_new[index] = x1[index1] + y_new[index] = y1[index1] + mp_new[index] = mp1[index1] + elif x1[index1+1] > x2[index2+1]: + index2 += 1 + index += 1 + x_new[index] = x2[index2] + y_new[index] = y2[index2] + mp_new[index] = mp2[index2] + else: # x1[index1+1] == x2[index2+1] + index1 += 1 + index2 += 1 + index += 1 + x_new[index] = x1[index1] + y_new[index] = y1[index1] + y2[index2] + mp_new[index] = mp1[index1] + mp2[index2] + # one array reached the end -> copy the contents of the other to the end + if index1+1 < len(y1): + x_new[index+1:index+1+len(x1)-index1-1] = x1[index1+1:] + y_new[index+1:index+1+len(x1)-index1-1] = y1[index1+1:] + mp_new[index+1:index+1+len(x1)-index1-1] = mp1[index1+1:] + index += len(x1)-index1-1 + elif index2+1 < len(y2): + x_new[index+1:index+1+len(x2)-index2-1] = x2[index2+1:] + y_new[index+1:index+1+len(x2)-index2-1] = y2[index2+1:] + mp_new[index+1:index+1+len(x2)-index2-1] = mp2[index2+1:] + index += len(x2)-index2-1 + # else: # both arrays reached the end simultaneously + # x_new[index+1] = x1[-1] + # y_new[index+1] = y1[-1] + y2[-1] + # mp_new[index+1] = mp1[-1] + mp2[-1] + + y_new[0] = y_new[1] + mp_new[0] = mp_new[1] + + # the last value is again the end of the interval + # only use the data that was actually filled + return (np.array(x_new[:index+1]), + np.array(y_new[:index+1]), + np.array(mp_new[:index+1])) diff --git a/pyspike/cython/cython_distance.pyx b/pyspike/cython/cython_distance.pyx new file mode 100644 index 0000000..489aab9 --- /dev/null +++ b/pyspike/cython/cython_distance.pyx @@ -0,0 +1,312 @@ +#cython: boundscheck=False +#cython: wraparound=False +#cython: cdivision=True + +""" +cython_distances.pyx + +cython implementation of the isi- and spike-distance + +Note: using cython memoryviews (e.g. double[:]) instead of ndarray objects +improves the performance of spike_distance by a factor of 10! + +Copyright 2014, Mario Mulansky + +Distributed under the BSD License + +""" + +""" +To test whether things can be optimized: remove all yellow stuff +in the html output:: + + cython -a cython_distance.pyx + +which gives:: + + cython_distance.html + +""" + +import numpy as np +cimport numpy as np + +from libc.math cimport fabs +from libc.math cimport fmax +from libc.math cimport fmin + +DTYPE = np.float +ctypedef np.float_t DTYPE_t + + +############################################################ +# isi_distance_cython +############################################################ +def isi_distance_cython(double[:] s1, + double[:] s2): + + cdef double[:] spike_events + cdef double[:] isi_values + cdef int index1, index2, index + cdef int N1, N2 + cdef double nu1, nu2 + N1 = len(s1)-1 + N2 = len(s2)-1 + + nu1 = s1[1]-s1[0] + nu2 = s2[1]-s2[0] + spike_events = np.empty(N1+N2) + spike_events[0] = s1[0] + # the values have one entry less - the number of intervals between events + isi_values = np.empty(N1+N2-1) + + with nogil: # release the interpreter to allow multithreading + isi_values[0] = fabs(nu1-nu2)/fmax(nu1, nu2) + index1 = 0 + index2 = 0 + index = 1 + while True: + # check which spike is next - from s1 or s2 + if s1[index1+1] < s2[index2+1]: + index1 += 1 + # break condition relies on existence of spikes at T_end + if index1 >= N1: + break + spike_events[index] = s1[index1] + nu1 = s1[index1+1]-s1[index1] + elif s1[index1+1] > s2[index2+1]: + index2 += 1 + if index2 >= N2: + break + spike_events[index] = s2[index2] + nu2 = s2[index2+1]-s2[index2] + else: # s1[index1+1] == s2[index2+1] + index1 += 1 + index2 += 1 + if (index1 >= N1) or (index2 >= N2): + break + spike_events[index] = s1[index1] + nu1 = s1[index1+1]-s1[index1] + nu2 = s2[index2+1]-s2[index2] + # compute the corresponding isi-distance + isi_values[index] = fabs(nu1 - nu2) / fmax(nu1, nu2) + index += 1 + # the last event is the interval end + spike_events[index] = s1[N1] + # end nogil + + return spike_events[:index+1], isi_values[:index] + + +############################################################ +# get_min_dist_cython +############################################################ +cdef inline double get_min_dist_cython(double spike_time, + double[:] spike_train, + # use memory view to ensure inlining + # np.ndarray[DTYPE_t,ndim=1] spike_train, + int N, + int start_index=0) nogil: + """ Returns the minimal distance |spike_time - spike_train[i]| + with i>=start_index. + """ + cdef double d, d_temp + d = fabs(spike_time - spike_train[start_index]) + start_index += 1 + while start_index < N: + d_temp = fabs(spike_time - spike_train[start_index]) + if d_temp > d: + break + else: + d = d_temp + start_index += 1 + return d + + +############################################################ +# isi_avrg_cython +############################################################ +cdef inline double isi_avrg_cython(double isi1, double isi2) nogil: + return 0.5*(isi1+isi2)*(isi1+isi2) + # alternative definition to obtain ~ 0.5 for Poisson spikes + # return 0.5*(isi1*isi1+isi2*isi2) + + +############################################################ +# spike_distance_cython +############################################################ +def spike_distance_cython(double[:] t1, + double[:] t2): + + cdef double[:] spike_events + cdef double[:] y_starts + cdef double[:] y_ends + + cdef int N1, N2, index1, index2, index + cdef double dt_p1, dt_p2, dt_f1, dt_f2, isi1, isi2, s1, s2 + + N1 = len(t1) + N2 = len(t2) + + spike_events = np.empty(N1+N2-2) + spike_events[0] = t1[0] + y_starts = np.empty(len(spike_events)-1) + y_ends = np.empty(len(spike_events)-1) + + with nogil: # release the interpreter to allow multithreading + index1 = 0 + index2 = 0 + index = 1 + dt_p1 = 0.0 + dt_f1 = get_min_dist_cython(t1[1], t2, N2, 0) + dt_p2 = 0.0 + dt_f2 = get_min_dist_cython(t2[1], t1, N1, 0) + isi1 = max(t1[1]-t1[0], t1[2]-t1[1]) + isi2 = max(t2[1]-t2[0], t2[2]-t2[1]) + s1 = dt_f1*(t1[1]-t1[0])/isi1 + s2 = dt_f2*(t2[1]-t2[0])/isi2 + y_starts[0] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2) + while True: + # print(index, index1, index2) + if t1[index1+1] < t2[index2+1]: + index1 += 1 + # break condition relies on existence of spikes at T_end + if index1+1 >= N1: + break + spike_events[index] = t1[index1] + # first calculate the previous interval end value + dt_p1 = dt_f1 # the previous time now was the following time before + s1 = dt_p1 + s2 = (dt_p2*(t2[index2+1]-t1[index1]) + + dt_f2*(t1[index1]-t2[index2])) / isi2 + y_ends[index-1] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2) + # now the next interval start value + dt_f1 = get_min_dist_cython(t1[index1+1], t2, N2, index2) + isi1 = t1[index1+1]-t1[index1] + # s2 is the same as above, thus we can compute y2 immediately + y_starts[index] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2) + elif t1[index1+1] > t2[index2+1]: + index2 += 1 + if index2+1 >= N2: + break + spike_events[index] = t2[index2] + # first calculate the previous interval end value + dt_p2 = dt_f2 # the previous time now was the following time before + s1 = (dt_p1*(t1[index1+1]-t2[index2]) + + dt_f1*(t2[index2]-t1[index1])) / isi1 + s2 = dt_p2 + y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2) + # now the next interval start value + dt_f2 = get_min_dist_cython(t2[index2+1], t1, N1, index1) + #s2 = dt_f2 + isi2 = t2[index2+1]-t2[index2] + # 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: # t1[index1+1] == t2[index2+1] - generate only one event + index1 += 1 + index2 += 1 + if (index1+1 >= N1) or (index2+1 >= N2): + break + spike_events[index] = t1[index1] + y_ends[index-1] = 0.0 + y_starts[index] = 0.0 + dt_p1 = 0.0 + dt_p2 = 0.0 + dt_f1 = get_min_dist_cython(t1[index1+1], t2, N2, index2) + dt_f2 = get_min_dist_cython(t2[index2+1], t1, N1, index1) + isi1 = t1[index1+1]-t1[index1] + isi2 = t2[index2+1]-t2[index2] + index += 1 + # the last event is the interval end + spike_events[index] = t1[N1-1] + # the ending value of the last interval + isi1 = max(t1[N1-1]-t1[N1-2], t1[N1-2]-t1[N1-3]) + isi2 = max(t2[N2-1]-t2[N2-2], t2[N2-2]-t2[N2-3]) + s1 = dt_p1*(t1[N1-1]-t1[N1-2])/isi1 + s2 = dt_p2*(t2[N2-1]-t2[N2-2])/isi2 + y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2) + # end nogil + + # use only the data added above + # could be less than original length due to equal spike times + return spike_events[:index+1], y_starts[:index], y_ends[:index] + + + +############################################################ +# coincidence_python +############################################################ +cdef inline double get_tau(double[:] spikes1, double[:] spikes2, int i, int j): + cdef double m = 1E100 # some huge number + cdef int N1 = len(spikes1)-2 + cdef int N2 = len(spikes2)-2 + if i < N1: + m = fmin(m, spikes1[i+1]-spikes1[i]) + if j < N2: + m = fmin(m, spikes2[j+1]-spikes2[j]) + if i > 1: + m = fmin(m, spikes1[i]-spikes1[i-1]) + if j > 1: + m = fmin(m, spikes2[j]-spikes2[j-1]) + return 0.5*m + + +############################################################ +# coincidence_cython +############################################################ +def coincidence_cython(double[:] spikes1, double[:] spikes2): + + cdef int N1 = len(spikes1) + cdef int N2 = len(spikes2) + cdef int i = 0 + cdef int j = 0 + cdef int n = 0 + cdef double[:] st = np.zeros(N1 + N2 - 2) # spike times + cdef double[:] c = np.zeros(N1 + N2 - 2) # coincidences + cdef double[:] mp = np.ones(N1 + N2 - 2) # multiplicity + cdef double tau + while n < N1 + N2 - 2: + if spikes1[i+1] < spikes2[j+1]: + i += 1 + n += 1 + tau = get_tau(spikes1, spikes2, i, j) + st[n] = spikes1[i] + if j > 0 and spikes1[i]-spikes2[j] < tau: + # coincidence between the current spike and the previous spike + # both get marked with 1 + c[n] = 1 + c[n-1] = 1 + elif spikes1[i+1] > spikes2[j+1]: + j += 1 + n += 1 + tau = get_tau(spikes1, spikes2, i, j) + st[n] = spikes2[j] + if i > 0 and spikes2[j]-spikes1[i] < tau: + # coincidence between the current spike and the previous spike + # both get marked with 1 + c[n] = 1 + c[n-1] = 1 + else: # spikes1[i+1] = spikes2[j+1] + # advance in both spike trains + j += 1 + i += 1 + if i == N1-1 or j == N2-1: + break + n += 1 + # add only one event, but with coincidence 2 and multiplicity 2 + st[n] = spikes1[i] + c[n] = 2 + mp[n] = 2 + + st = st[:n+2] + c = c[:n+2] + mp = mp[:n+2] + + st[0] = spikes1[0] + st[len(st)-1] = spikes1[len(spikes1)-1] + c[0] = c[1] + c[len(c)-1] = c[len(c)-2] + mp[0] = mp[1] + mp[len(mp)-1] = mp[len(mp)-2] + + return st, c, mp diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py new file mode 100644 index 0000000..481daf9 --- /dev/null +++ b/pyspike/cython/python_backend.py @@ -0,0 +1,485 @@ +""" python_backend.py + +Collection of python functions that can be used instead of the cython +implementation. + +Copyright 2014, Mario Mulansky + +Distributed under the BSD License + +""" + +import numpy as np + + +############################################################ +# isi_distance_python +############################################################ +def isi_distance_python(s1, s2): + """ Plain Python implementation of the isi distance. + """ + # compute the interspike interval + nu1 = s1[1:] - s1[:-1] + nu2 = s2[1:] - s2[:-1] + + # compute the isi-distance + spike_events = np.empty(len(nu1) + len(nu2)) + spike_events[0] = s1[0] + # the values have one entry less - the number of intervals between events + isi_values = np.empty(len(spike_events) - 1) + # add the distance of the first events + # isi_values[0] = nu1[0]/nu2[0] - 1.0 if nu1[0] <= nu2[0] \ + # else 1.0 - nu2[0]/nu1[0] + isi_values[0] = abs(nu1[0] - nu2[0]) / max(nu1[0], nu2[0]) + index1 = 0 + index2 = 0 + index = 1 + while True: + # check which spike is next - from s1 or s2 + if s1[index1+1] < s2[index2+1]: + index1 += 1 + # break condition relies on existence of spikes at T_end + if index1 >= len(nu1): + break + spike_events[index] = s1[index1] + elif s1[index1+1] > s2[index2+1]: + index2 += 1 + if index2 >= len(nu2): + break + spike_events[index] = s2[index2] + else: # s1[index1 + 1] == s2[index2 + 1] + index1 += 1 + index2 += 1 + if (index1 >= len(nu1)) or (index2 >= len(nu2)): + break + spike_events[index] = s1[index1] + # compute the corresponding isi-distance + isi_values[index] = abs(nu1[index1] - nu2[index2]) / \ + max(nu1[index1], nu2[index2]) + index += 1 + # the last event is the interval end + spike_events[index] = s1[-1] + # use only the data added above + # could be less than original length due to equal spike times + return spike_events[:index + 1], isi_values[:index] + + +############################################################ +# get_min_dist +############################################################ +def get_min_dist(spike_time, spike_train, start_index=0): + """ Returns the minimal distance |spike_time - spike_train[i]| + with i>=start_index. + """ + d = abs(spike_time - spike_train[start_index]) + start_index += 1 + while start_index < len(spike_train): + d_temp = abs(spike_time - spike_train[start_index]) + if d_temp > d: + break + else: + d = d_temp + start_index += 1 + return d + + +############################################################ +# spike_distance_python +############################################################ +def spike_distance_python(spikes1, spikes2): + """ Computes the instantaneous spike-distance S_spike (t) of the two given + spike trains. 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. + Args: + - spikes1, spikes2: ordered arrays of spike times with auxiliary spikes. + Returns: + - PieceWiseLinFunc describing the spike-distance. + """ + # 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!" + # shorter variables + t1 = spikes1 + t2 = spikes2 + + spike_events = np.empty(len(t1) + len(t2) - 2) + spike_events[0] = t1[0] + y_starts = np.empty(len(spike_events) - 1) + y_ends = np.empty(len(spike_events) - 1) + + index1 = 0 + index2 = 0 + index = 1 + dt_p1 = 0.0 + dt_f1 = get_min_dist(t1[1], t2, 0) + dt_p2 = 0.0 + dt_f2 = get_min_dist(t2[1], t1, 0) + isi1 = max(t1[1]-t1[0], t1[2]-t1[1]) + isi2 = max(t2[1]-t2[0], t2[2]-t2[1]) + s1 = dt_f1*(t1[1]-t1[0])/isi1 + s2 = dt_f2*(t2[1]-t2[0])/isi2 + y_starts[0] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2) + while True: + # print(index, index1, index2) + if t1[index1+1] < t2[index2+1]: + index1 += 1 + # break condition relies on existence of spikes at T_end + if index1+1 >= len(t1): + break + spike_events[index] = t1[index1] + # first calculate the previous interval end value + dt_p1 = dt_f1 # the previous time was the following time before + s1 = dt_p1 + s2 = (dt_p2*(t2[index2+1]-t1[index1]) + + dt_f2*(t1[index1]-t2[index2])) / isi2 + y_ends[index-1] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2) + # now the next interval start value + dt_f1 = get_min_dist(t1[index1+1], t2, index2) + isi1 = t1[index1+1]-t1[index1] + # s2 is the same as above, thus we can compute y2 immediately + y_starts[index] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2) + elif t1[index1+1] > t2[index2+1]: + index2 += 1 + if index2+1 >= len(t2): + break + spike_events[index] = t2[index2] + # first calculate the previous interval end value + dt_p2 = dt_f2 # the previous time was the following time before + s1 = (dt_p1*(t1[index1+1]-t2[index2]) + + dt_f1*(t2[index2]-t1[index1])) / isi1 + s2 = dt_p2 + y_ends[index-1] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2) + # now the next interval start value + dt_f2 = get_min_dist(t2[index2+1], t1, index1) + #s2 = dt_f2 + isi2 = t2[index2+1]-t2[index2] + # s2 is the same as above, thus we can compute y2 immediately + y_starts[index] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2) + else: # t1[index1+1] == t2[index2+1] - generate only one event + index1 += 1 + index2 += 1 + if (index1+1 >= len(t1)) or (index2+1 >= len(t2)): + break + assert dt_f2 == 0.0 + assert dt_f1 == 0.0 + spike_events[index] = t1[index1] + y_ends[index-1] = 0.0 + y_starts[index] = 0.0 + dt_p1 = 0.0 + dt_p2 = 0.0 + dt_f1 = get_min_dist(t1[index1+1], t2, index2) + dt_f2 = get_min_dist(t2[index2+1], t1, index1) + isi1 = t1[index1+1]-t1[index1] + isi2 = t2[index2+1]-t2[index2] + index += 1 + # the last event is the interval end + spike_events[index] = t1[-1] + # the ending value of the last interval + isi1 = max(t1[-1]-t1[-2], t1[-2]-t1[-3]) + isi2 = max(t2[-1]-t2[-2], t2[-2]-t2[-3]) + s1 = dt_p1*(t1[-1]-t1[-2])/isi1 + s2 = dt_p2*(t2[-1]-t2[-2])/isi2 + y_ends[index-1] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2) + # use only the data added above + # could be less than original length due to equal spike times + return spike_events[:index+1], y_starts[:index], y_ends[:index] + + +############################################################ +# cumulative_sync_python +############################################################ +def cumulative_sync_python(spikes1, spikes2): + + def get_tau(spikes1, spikes2, i, j): + return 0.5*min([spikes1[i]-spikes1[i-1], spikes1[i+1]-spikes1[i], + spikes2[j]-spikes2[j-1], spikes2[j+1]-spikes2[j]]) + N1 = len(spikes1) + N2 = len(spikes2) + i = 0 + j = 0 + n = 0 + st = np.zeros(N1 + N2 - 2) + c = np.zeros(N1 + N2 - 3) + c[0] = 0 + st[0] = 0 + while n < N1 + N2: + if spikes1[i+1] < spikes2[j+1]: + i += 1 + n += 1 + tau = get_tau(spikes1, spikes2, i, j) + st[n] = spikes1[i] + if spikes1[i]-spikes2[j] > tau: + c[n] = c[n-1] + else: + c[n] = c[n-1]+1 + elif spikes1[i+1] > spikes2[j+1]: + j += 1 + n += 1 + tau = get_tau(spikes1, spikes2, i, j) + st[n] = spikes2[j] + if spikes2[j]-spikes1[i] > tau: + c[n] = c[n-1] + else: + c[n] = c[n-1]+1 + else: # spikes1[i+1] = spikes2[j+1] + j += 1 + i += 1 + if i == N1-1 or j == N2-1: + break + n += 1 + st[n] = spikes1[i] + c[n] = c[n-1] + n += 1 + st[n] = spikes1[i] + c[n] = c[n-1]+1 + c[0] = 0 + st[0] = spikes1[0] + st[-1] = spikes1[-1] + + return st, c + + +############################################################ +# coincidence_python +############################################################ +def coincidence_python(spikes1, spikes2): + + def get_tau(spikes1, spikes2, i, j): + m = 1E100 # some huge number + if i < len(spikes1)-2: + m = min(m, spikes1[i+1]-spikes1[i]) + if j < len(spikes2)-2: + m = min(m, spikes2[j+1]-spikes2[j]) + if i > 1: + m = min(m, spikes1[i]-spikes1[i-1]) + if j > 1: + m = min(m, spikes2[j]-spikes2[j-1]) + return 0.5*m + N1 = len(spikes1) + N2 = len(spikes2) + i = 0 + j = 0 + n = 0 + st = np.zeros(N1 + N2 - 2) # spike times + c = np.zeros(N1 + N2 - 2) # coincidences + mp = np.ones(N1 + N2 - 2) # multiplicity + while n < N1 + N2 - 2: + if spikes1[i+1] < spikes2[j+1]: + i += 1 + n += 1 + tau = get_tau(spikes1, spikes2, i, j) + st[n] = spikes1[i] + if j > 0 and spikes1[i]-spikes2[j] < tau: + # coincidence between the current spike and the previous spike + # both get marked with 1 + c[n] = 1 + c[n-1] = 1 + elif spikes1[i+1] > spikes2[j+1]: + j += 1 + n += 1 + tau = get_tau(spikes1, spikes2, i, j) + st[n] = spikes2[j] + if i > 0 and spikes2[j]-spikes1[i] < tau: + # coincidence between the current spike and the previous spike + # both get marked with 1 + c[n] = 1 + c[n-1] = 1 + else: # spikes1[i+1] = spikes2[j+1] + # advance in both spike trains + j += 1 + i += 1 + if i == N1-1 or j == N2-1: + break + n += 1 + # add only one event, but with coincidence 2 and multiplicity 2 + st[n] = spikes1[i] + c[n] = 2 + mp[n] = 2 + + st = st[:n+2] + c = c[:n+2] + mp = mp[:n+2] + + st[0] = spikes1[0] + st[-1] = spikes1[-1] + c[0] = c[1] + c[-1] = c[-2] + mp[0] = mp[1] + mp[-1] = mp[-2] + + return st, c, mp + + +############################################################ +# add_piece_wise_const_python +############################################################ +def add_piece_wise_const_python(x1, y1, x2, y2): + x_new = np.empty(len(x1) + len(x2)) + y_new = np.empty(len(x_new)-1) + x_new[0] = x1[0] + y_new[0] = y1[0] + y2[0] + index1 = 0 + index2 = 0 + index = 0 + while (index1+1 < len(y1)) and (index2+1 < len(y2)): + index += 1 + # print(index1+1, x1[index1+1], y1[index1+1], x_new[index]) + if x1[index1+1] < x2[index2+1]: + index1 += 1 + x_new[index] = x1[index1] + elif x1[index1+1] > x2[index2+1]: + index2 += 1 + x_new[index] = x2[index2] + else: # x1[index1+1] == x2[index2+1]: + index1 += 1 + index2 += 1 + x_new[index] = x1[index1] + y_new[index] = y1[index1] + y2[index2] + # one array reached the end -> copy the contents of the other to the end + if index1+1 < len(y1): + x_new[index+1:index+1+len(x1)-index1-1] = x1[index1+1:] + y_new[index+1:index+1+len(y1)-index1-1] = y1[index1+1:] + y2[-1] + index += len(x1)-index1-2 + elif index2+1 < len(y2): + x_new[index+1:index+1+len(x2)-index2-1] = x2[index2+1:] + y_new[index+1:index+1+len(y2)-index2-1] = y2[index2+1:] + y1[-1] + index += len(x2)-index2-2 + else: # both arrays reached the end simultaneously + # only the last x-value missing + x_new[index+1] = x1[-1] + # the last value is again the end of the interval + # x_new[index+1] = x1[-1] + # only use the data that was actually filled + + return x_new[:index+2], y_new[:index+1] + + +############################################################ +# add_piece_lin_const_python +############################################################ +def add_piece_wise_lin_python(x1, y11, y12, x2, y21, y22): + x_new = np.empty(len(x1) + len(x2)) + y1_new = np.empty(len(x_new)-1) + y2_new = np.empty_like(y1_new) + x_new[0] = x1[0] + y1_new[0] = y11[0] + y21[0] + index1 = 0 # index for self + index2 = 0 # index for f + index = 0 # index for new + while (index1+1 < len(y11)) and (index2+1 < len(y21)): + # print(index1+1, x1[index1+1], self.y[index1+1], x_new[index]) + if x1[index1+1] < x2[index2+1]: + # first compute the end value of the previous interval + # linear interpolation of the interval + y = y21[index2] + (y22[index2]-y21[index2]) * \ + (x1[index1+1]-x2[index2]) / (x2[index2+1]-x2[index2]) + y2_new[index] = y12[index1] + y + index1 += 1 + index += 1 + x_new[index] = x1[index1] + # and the starting value for the next interval + y1_new[index] = y11[index1] + y + elif x1[index1+1] > x2[index2+1]: + # first compute the end value of the previous interval + # linear interpolation of the interval + y = y11[index1] + (y12[index1]-y11[index1]) * \ + (x2[index2+1]-x1[index1]) / \ + (x1[index1+1]-x1[index1]) + y2_new[index] = y22[index2] + y + index2 += 1 + index += 1 + x_new[index] = x2[index2] + # and the starting value for the next interval + y1_new[index] = y21[index2] + y + else: # x1[index1+1] == x2[index2+1]: + y2_new[index] = y12[index1] + y22[index2] + index1 += 1 + index2 += 1 + index += 1 + x_new[index] = x1[index1] + y1_new[index] = y11[index1] + y21[index2] + # one array reached the end -> copy the contents of the other to the end + if index1+1 < len(y11): + # compute the linear interpolations values + y = y21[index2] + (y22[index2]-y21[index2]) * \ + (x1[index1+1:-1]-x2[index2]) / (x2[index2+1]-x2[index2]) + x_new[index+1:index+1+len(x1)-index1-1] = x1[index1+1:] + y1_new[index+1:index+1+len(y11)-index1-1] = y11[index1+1:]+y + y2_new[index:index+len(y12)-index1-1] = y12[index1:-1] + y + index += len(x1)-index1-2 + elif index2+1 < len(y21): + # compute the linear interpolations values + y = y11[index1] + (y12[index1]-y11[index1]) * \ + (x2[index2+1:-1]-x1[index1]) / \ + (x1[index1+1]-x1[index1]) + x_new[index+1:index+1+len(x2)-index2-1] = x2[index2+1:] + y1_new[index+1:index+1+len(y21)-index2-1] = y21[index2+1:] + y + y2_new[index:index+len(y22)-index2-1] = y22[index2:-1] + y + index += len(x2)-index2-2 + else: # both arrays reached the end simultaneously + # only the last x-value missing + x_new[index+1] = x1[-1] + # finally, the end value for the last interval + y2_new[index] = y12[-1]+y22[-1] + # only use the data that was actually filled + return x_new[:index+2], y1_new[:index+1], y2_new[:index+1] + + +############################################################ +# add_discrete_function_python +############################################################ +def add_discrete_function_python(x1, y1, mp1, x2, y2, mp2): + + x_new = np.empty(len(x1) + len(x2)) + y_new = np.empty_like(x_new) + mp_new = np.empty_like(x_new) + x_new[0] = x1[0] + index1 = 0 + index2 = 0 + index = 0 + while (index1+1 < len(y1)) and (index2+1 < len(y2)): + if x1[index1+1] < x2[index2+1]: + index1 += 1 + index += 1 + x_new[index] = x1[index1] + y_new[index] = y1[index1] + mp_new[index] = mp1[index1] + elif x1[index1+1] > x2[index2+1]: + index2 += 1 + index += 1 + x_new[index] = x2[index2] + y_new[index] = y2[index2] + mp_new[index] = mp2[index2] + else: # x1[index1+1] == x2[index2+1] + index1 += 1 + index2 += 1 + index += 1 + x_new[index] = x1[index1] + y_new[index] = y1[index1] + y2[index2] + mp_new[index] = mp1[index1] + mp2[index2] + # one array reached the end -> copy the contents of the other to the end + if index1+1 < len(y1): + x_new[index+1:index+1+len(x1)-index1-1] = x1[index1+1:] + y_new[index+1:index+1+len(x1)-index1-1] = y1[index1+1:] + mp_new[index+1:index+1+len(x1)-index1-1] = mp1[index1+1:] + index += len(x1)-index1-1 + elif index2+1 < len(y2): + x_new[index+1:index+1+len(x2)-index2-1] = x2[index2+1:] + y_new[index+1:index+1+len(x2)-index2-1] = y2[index2+1:] + mp_new[index+1:index+1+len(x2)-index2-1] = mp2[index2+1:] + index += len(x2)-index2-1 + # else: # both arrays reached the end simultaneously + # x_new[index+1] = x1[-1] + # y_new[index+1] = y1[-1] + y2[-1] + # mp_new[index+1] = mp1[-1] + mp2[-1] + + y_new[0] = y_new[1] + mp_new[0] = mp_new[1] + + # the last value is again the end of the interval + # only use the data that was actually filled + return x_new[:index+1], y_new[:index+1], mp_new[:index+1] + diff --git a/pyspike/cython_add.pyx b/pyspike/cython_add.pyx deleted file mode 100644 index ac64005..0000000 --- a/pyspike/cython_add.pyx +++ /dev/null @@ -1,235 +0,0 @@ -#cython: boundscheck=False -#cython: wraparound=False -#cython: cdivision=True - -""" -cython_add.pyx - -cython implementation of the add function for piece-wise const and -piece-wise linear functions - -Note: using cython memoryviews (e.g. double[:]) instead of ndarray objects -improves the performance of spike_distance by a factor of 10! - -Copyright 2014, Mario Mulansky - -Distributed under the BSD License - -""" - -""" -To test whether things can be optimized: remove all yellow stuff -in the html output:: - - cython -a cython_add.pyx - -which gives:: - - cython_add.html - -""" - -import numpy as np -cimport numpy as np - -from libc.math cimport fabs - -DTYPE = np.float -ctypedef np.float_t DTYPE_t - - -############################################################ -# add_piece_wise_const_cython -############################################################ -def add_piece_wise_const_cython(double[:] x1, double[:] y1, - double[:] x2, double[:] y2): - - cdef int N1 = len(x1) - cdef int N2 = len(x2) - cdef double[:] x_new = np.empty(N1+N2) - cdef double[:] y_new = np.empty(N1+N2-1) - cdef int index1 = 0 - cdef int index2 = 0 - cdef int index = 0 - cdef int i - with nogil: # release the interpreter lock to allow multi-threading - x_new[0] = x1[0] - y_new[0] = y1[0] + y2[0] - while (index1+1 < N1-1) and (index2+1 < N2-1): - index += 1 - # print(index1+1, x1[index1+1], y1[index1+1], x_new[index]) - if x1[index1+1] < x2[index2+1]: - index1 += 1 - x_new[index] = x1[index1] - elif x1[index1+1] > x2[index2+1]: - index2 += 1 - x_new[index] = x2[index2] - else: # x1[index1+1] == x2[index2+1]: - index1 += 1 - index2 += 1 - x_new[index] = x1[index1] - y_new[index] = y1[index1] + y2[index2] - # one array reached the end -> copy the contents of the other to the end - if index1+1 < N1-1: - x_new[index+1:index+1+N1-index1-1] = x1[index1+1:] - for i in xrange(N1-index1-2): - y_new[index+1+i] = y1[index1+1+i] + y2[N2-2] - index += N1-index1-2 - elif index2+1 < N2-1: - x_new[index+1:index+1+N2-index2-1] = x2[index2+1:] - for i in xrange(N2-index2-2): - y_new[index+1+i] = y2[index2+1+i] + y1[N1-2] - index += N2-index2-2 - else: # both arrays reached the end simultaneously - # only the last x-value missing - x_new[index+1] = x1[N1-1] - # the last value is again the end of the interval - # x_new[index+1] = x1[-1] - # only use the data that was actually filled - x1 = x_new[:index+2] - y1 = y_new[:index+1] - # end nogil - return np.array(x_new[:index+2]), np.array(y_new[:index+1]) - - -############################################################ -# add_piece_wise_lin_cython -############################################################ -def add_piece_wise_lin_cython(double[:] x1, double[:] y11, double[:] y12, - double[:] x2, double[:] y21, double[:] y22): - cdef int N1 = len(x1) - cdef int N2 = len(x2) - cdef double[:] x_new = np.empty(N1+N2) - cdef double[:] y1_new = np.empty(N1+N2-1) - cdef double[:] y2_new = np.empty_like(y1_new) - cdef int index1 = 0 # index for self - cdef int index2 = 0 # index for f - cdef int index = 0 # index for new - cdef int i - cdef double y - with nogil: # release the interpreter lock to allow multi-threading - x_new[0] = x1[0] - y1_new[0] = y11[0] + y21[0] - while (index1+1 < N1-1) and (index2+1 < N2-1): - # print(index1+1, x1[index1+1], self.y[index1+1], x_new[index]) - if x1[index1+1] < x2[index2+1]: - # first compute the end value of the previous interval - # linear interpolation of the interval - y = y21[index2] + (y22[index2]-y21[index2]) * \ - (x1[index1+1]-x2[index2]) / (x2[index2+1]-x2[index2]) - y2_new[index] = y12[index1] + y - index1 += 1 - index += 1 - x_new[index] = x1[index1] - # and the starting value for the next interval - y1_new[index] = y11[index1] + y - elif x1[index1+1] > x2[index2+1]: - # first compute the end value of the previous interval - # linear interpolation of the interval - y = y11[index1] + (y12[index1]-y11[index1]) * \ - (x2[index2+1]-x1[index1]) / \ - (x1[index1+1]-x1[index1]) - y2_new[index] = y22[index2] + y - index2 += 1 - index += 1 - x_new[index] = x2[index2] - # and the starting value for the next interval - y1_new[index] = y21[index2] + y - else: # x1[index1+1] == x2[index2+1]: - y2_new[index] = y12[index1] + y22[index2] - index1 += 1 - index2 += 1 - index += 1 - x_new[index] = x1[index1] - y1_new[index] = y11[index1] + y21[index2] - # one array reached the end -> copy the contents of the other to the end - if index1+1 < N1-1: - x_new[index+1:index+1+N1-index1-1] = x1[index1+1:] - for i in xrange(N1-index1-2): - # compute the linear interpolations value - y = y21[index2] + (y22[index2]-y21[index2]) * \ - (x1[index1+1+i]-x2[index2]) / (x2[index2+1]-x2[index2]) - y1_new[index+1+i] = y11[index1+1+i] + y - y2_new[index+i] = y12[index1+i] + y - index += N1-index1-2 - elif index2+1 < N2-1: - x_new[index+1:index+1+N2-index2-1] = x2[index2+1:] - # compute the linear interpolations values - for i in xrange(N2-index2-2): - y = y11[index1] + (y12[index1]-y11[index1]) * \ - (x2[index2+1+i]-x1[index1]) / \ - (x1[index1+1]-x1[index1]) - y1_new[index+1+i] = y21[index2+1+i] + y - y2_new[index+i] = y22[index2+i] + y - index += N2-index2-2 - else: # both arrays reached the end simultaneously - # only the last x-value missing - x_new[index+1] = x1[N1-1] - # finally, the end value for the last interval - y2_new[index] = y12[N1-2]+y22[N2-2] - # only use the data that was actually filled - # end nogil - return (np.array(x_new[:index+2]), - np.array(y1_new[:index+1]), - np.array(y2_new[:index+1])) - - -############################################################ -# add_discrete_function_cython -############################################################ -def add_discrete_function_cython(double[:] x1, double[:] y1, double[:] mp1, - double[:] x2, double[:] y2, double[:] mp2): - - cdef double[:] x_new = np.empty(len(x1) + len(x2)) - cdef double[:] y_new = np.empty_like(x_new) - cdef double[:] mp_new = np.empty_like(x_new) - cdef int index1 = 0 - cdef int index2 = 0 - cdef int index = 0 - cdef int N1 = len(y1) - cdef int N2 = len(y2) - x_new[0] = x1[0] - while (index1+1 < N1) and (index2+1 < N2): - if x1[index1+1] < x2[index2+1]: - index1 += 1 - index += 1 - x_new[index] = x1[index1] - y_new[index] = y1[index1] - mp_new[index] = mp1[index1] - elif x1[index1+1] > x2[index2+1]: - index2 += 1 - index += 1 - x_new[index] = x2[index2] - y_new[index] = y2[index2] - mp_new[index] = mp2[index2] - else: # x1[index1+1] == x2[index2+1] - index1 += 1 - index2 += 1 - index += 1 - x_new[index] = x1[index1] - y_new[index] = y1[index1] + y2[index2] - mp_new[index] = mp1[index1] + mp2[index2] - # one array reached the end -> copy the contents of the other to the end - if index1+1 < len(y1): - x_new[index+1:index+1+len(x1)-index1-1] = x1[index1+1:] - y_new[index+1:index+1+len(x1)-index1-1] = y1[index1+1:] - mp_new[index+1:index+1+len(x1)-index1-1] = mp1[index1+1:] - index += len(x1)-index1-1 - elif index2+1 < len(y2): - x_new[index+1:index+1+len(x2)-index2-1] = x2[index2+1:] - y_new[index+1:index+1+len(x2)-index2-1] = y2[index2+1:] - mp_new[index+1:index+1+len(x2)-index2-1] = mp2[index2+1:] - index += len(x2)-index2-1 - # else: # both arrays reached the end simultaneously - # x_new[index+1] = x1[-1] - # y_new[index+1] = y1[-1] + y2[-1] - # mp_new[index+1] = mp1[-1] + mp2[-1] - - y_new[0] = y_new[1] - mp_new[0] = mp_new[1] - - # the last value is again the end of the interval - # only use the data that was actually filled - return (np.array(x_new[:index+1]), - np.array(y_new[:index+1]), - np.array(mp_new[:index+1])) diff --git a/pyspike/cython_distance.pyx b/pyspike/cython_distance.pyx deleted file mode 100644 index 489aab9..0000000 --- a/pyspike/cython_distance.pyx +++ /dev/null @@ -1,312 +0,0 @@ -#cython: boundscheck=False -#cython: wraparound=False -#cython: cdivision=True - -""" -cython_distances.pyx - -cython implementation of the isi- and spike-distance - -Note: using cython memoryviews (e.g. double[:]) instead of ndarray objects -improves the performance of spike_distance by a factor of 10! - -Copyright 2014, Mario Mulansky - -Distributed under the BSD License - -""" - -""" -To test whether things can be optimized: remove all yellow stuff -in the html output:: - - cython -a cython_distance.pyx - -which gives:: - - cython_distance.html - -""" - -import numpy as np -cimport numpy as np - -from libc.math cimport fabs -from libc.math cimport fmax -from libc.math cimport fmin - -DTYPE = np.float -ctypedef np.float_t DTYPE_t - - -############################################################ -# isi_distance_cython -############################################################ -def isi_distance_cython(double[:] s1, - double[:] s2): - - cdef double[:] spike_events - cdef double[:] isi_values - cdef int index1, index2, index - cdef int N1, N2 - cdef double nu1, nu2 - N1 = len(s1)-1 - N2 = len(s2)-1 - - nu1 = s1[1]-s1[0] - nu2 = s2[1]-s2[0] - spike_events = np.empty(N1+N2) - spike_events[0] = s1[0] - # the values have one entry less - the number of intervals between events - isi_values = np.empty(N1+N2-1) - - with nogil: # release the interpreter to allow multithreading - isi_values[0] = fabs(nu1-nu2)/fmax(nu1, nu2) - index1 = 0 - index2 = 0 - index = 1 - while True: - # check which spike is next - from s1 or s2 - if s1[index1+1] < s2[index2+1]: - index1 += 1 - # break condition relies on existence of spikes at T_end - if index1 >= N1: - break - spike_events[index] = s1[index1] - nu1 = s1[index1+1]-s1[index1] - elif s1[index1+1] > s2[index2+1]: - index2 += 1 - if index2 >= N2: - break - spike_events[index] = s2[index2] - nu2 = s2[index2+1]-s2[index2] - else: # s1[index1+1] == s2[index2+1] - index1 += 1 - index2 += 1 - if (index1 >= N1) or (index2 >= N2): - break - spike_events[index] = s1[index1] - nu1 = s1[index1+1]-s1[index1] - nu2 = s2[index2+1]-s2[index2] - # compute the corresponding isi-distance - isi_values[index] = fabs(nu1 - nu2) / fmax(nu1, nu2) - index += 1 - # the last event is the interval end - spike_events[index] = s1[N1] - # end nogil - - return spike_events[:index+1], isi_values[:index] - - -############################################################ -# get_min_dist_cython -############################################################ -cdef inline double get_min_dist_cython(double spike_time, - double[:] spike_train, - # use memory view to ensure inlining - # np.ndarray[DTYPE_t,ndim=1] spike_train, - int N, - int start_index=0) nogil: - """ Returns the minimal distance |spike_time - spike_train[i]| - with i>=start_index. - """ - cdef double d, d_temp - d = fabs(spike_time - spike_train[start_index]) - start_index += 1 - while start_index < N: - d_temp = fabs(spike_time - spike_train[start_index]) - if d_temp > d: - break - else: - d = d_temp - start_index += 1 - return d - - -############################################################ -# isi_avrg_cython -############################################################ -cdef inline double isi_avrg_cython(double isi1, double isi2) nogil: - return 0.5*(isi1+isi2)*(isi1+isi2) - # alternative definition to obtain ~ 0.5 for Poisson spikes - # return 0.5*(isi1*isi1+isi2*isi2) - - -############################################################ -# spike_distance_cython -############################################################ -def spike_distance_cython(double[:] t1, - double[:] t2): - - cdef double[:] spike_events - cdef double[:] y_starts - cdef double[:] y_ends - - cdef int N1, N2, index1, index2, index - cdef double dt_p1, dt_p2, dt_f1, dt_f2, isi1, isi2, s1, s2 - - N1 = len(t1) - N2 = len(t2) - - spike_events = np.empty(N1+N2-2) - spike_events[0] = t1[0] - y_starts = np.empty(len(spike_events)-1) - y_ends = np.empty(len(spike_events)-1) - - with nogil: # release the interpreter to allow multithreading - index1 = 0 - index2 = 0 - index = 1 - dt_p1 = 0.0 - dt_f1 = get_min_dist_cython(t1[1], t2, N2, 0) - dt_p2 = 0.0 - dt_f2 = get_min_dist_cython(t2[1], t1, N1, 0) - isi1 = max(t1[1]-t1[0], t1[2]-t1[1]) - isi2 = max(t2[1]-t2[0], t2[2]-t2[1]) - s1 = dt_f1*(t1[1]-t1[0])/isi1 - s2 = dt_f2*(t2[1]-t2[0])/isi2 - y_starts[0] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2) - while True: - # print(index, index1, index2) - if t1[index1+1] < t2[index2+1]: - index1 += 1 - # break condition relies on existence of spikes at T_end - if index1+1 >= N1: - break - spike_events[index] = t1[index1] - # first calculate the previous interval end value - dt_p1 = dt_f1 # the previous time now was the following time before - s1 = dt_p1 - s2 = (dt_p2*(t2[index2+1]-t1[index1]) + - dt_f2*(t1[index1]-t2[index2])) / isi2 - y_ends[index-1] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2) - # now the next interval start value - dt_f1 = get_min_dist_cython(t1[index1+1], t2, N2, index2) - isi1 = t1[index1+1]-t1[index1] - # s2 is the same as above, thus we can compute y2 immediately - y_starts[index] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2) - elif t1[index1+1] > t2[index2+1]: - index2 += 1 - if index2+1 >= N2: - break - spike_events[index] = t2[index2] - # first calculate the previous interval end value - dt_p2 = dt_f2 # the previous time now was the following time before - s1 = (dt_p1*(t1[index1+1]-t2[index2]) + - dt_f1*(t2[index2]-t1[index1])) / isi1 - s2 = dt_p2 - y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2) - # now the next interval start value - dt_f2 = get_min_dist_cython(t2[index2+1], t1, N1, index1) - #s2 = dt_f2 - isi2 = t2[index2+1]-t2[index2] - # 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: # t1[index1+1] == t2[index2+1] - generate only one event - index1 += 1 - index2 += 1 - if (index1+1 >= N1) or (index2+1 >= N2): - break - spike_events[index] = t1[index1] - y_ends[index-1] = 0.0 - y_starts[index] = 0.0 - dt_p1 = 0.0 - dt_p2 = 0.0 - dt_f1 = get_min_dist_cython(t1[index1+1], t2, N2, index2) - dt_f2 = get_min_dist_cython(t2[index2+1], t1, N1, index1) - isi1 = t1[index1+1]-t1[index1] - isi2 = t2[index2+1]-t2[index2] - index += 1 - # the last event is the interval end - spike_events[index] = t1[N1-1] - # the ending value of the last interval - isi1 = max(t1[N1-1]-t1[N1-2], t1[N1-2]-t1[N1-3]) - isi2 = max(t2[N2-1]-t2[N2-2], t2[N2-2]-t2[N2-3]) - s1 = dt_p1*(t1[N1-1]-t1[N1-2])/isi1 - s2 = dt_p2*(t2[N2-1]-t2[N2-2])/isi2 - y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2) - # end nogil - - # use only the data added above - # could be less than original length due to equal spike times - return spike_events[:index+1], y_starts[:index], y_ends[:index] - - - -############################################################ -# coincidence_python -############################################################ -cdef inline double get_tau(double[:] spikes1, double[:] spikes2, int i, int j): - cdef double m = 1E100 # some huge number - cdef int N1 = len(spikes1)-2 - cdef int N2 = len(spikes2)-2 - if i < N1: - m = fmin(m, spikes1[i+1]-spikes1[i]) - if j < N2: - m = fmin(m, spikes2[j+1]-spikes2[j]) - if i > 1: - m = fmin(m, spikes1[i]-spikes1[i-1]) - if j > 1: - m = fmin(m, spikes2[j]-spikes2[j-1]) - return 0.5*m - - -############################################################ -# coincidence_cython -############################################################ -def coincidence_cython(double[:] spikes1, double[:] spikes2): - - cdef int N1 = len(spikes1) - cdef int N2 = len(spikes2) - cdef int i = 0 - cdef int j = 0 - cdef int n = 0 - cdef double[:] st = np.zeros(N1 + N2 - 2) # spike times - cdef double[:] c = np.zeros(N1 + N2 - 2) # coincidences - cdef double[:] mp = np.ones(N1 + N2 - 2) # multiplicity - cdef double tau - while n < N1 + N2 - 2: - if spikes1[i+1] < spikes2[j+1]: - i += 1 - n += 1 - tau = get_tau(spikes1, spikes2, i, j) - st[n] = spikes1[i] - if j > 0 and spikes1[i]-spikes2[j] < tau: - # coincidence between the current spike and the previous spike - # both get marked with 1 - c[n] = 1 - c[n-1] = 1 - elif spikes1[i+1] > spikes2[j+1]: - j += 1 - n += 1 - tau = get_tau(spikes1, spikes2, i, j) - st[n] = spikes2[j] - if i > 0 and spikes2[j]-spikes1[i] < tau: - # coincidence between the current spike and the previous spike - # both get marked with 1 - c[n] = 1 - c[n-1] = 1 - else: # spikes1[i+1] = spikes2[j+1] - # advance in both spike trains - j += 1 - i += 1 - if i == N1-1 or j == N2-1: - break - n += 1 - # add only one event, but with coincidence 2 and multiplicity 2 - st[n] = spikes1[i] - c[n] = 2 - mp[n] = 2 - - st = st[:n+2] - c = c[:n+2] - mp = mp[:n+2] - - st[0] = spikes1[0] - st[len(st)-1] = spikes1[len(spikes1)-1] - c[0] = c[1] - c[len(c)-1] = c[len(c)-2] - mp[0] = mp[1] - mp[len(mp)-1] = mp[len(mp)-2] - - return st, c, mp diff --git a/pyspike/isi_distance.py b/pyspike/isi_distance.py index 745d280..c2ef8e8 100644 --- a/pyspike/isi_distance.py +++ b/pyspike/isi_distance.py @@ -35,13 +35,15 @@ def isi_profile(spikes1, spikes2): # load cython implementation try: - from cython_distance import isi_distance_cython as isi_distance_impl + 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 python_backend import isi_distance_python as isi_distance_impl + from cython.python_backend import isi_distance_python \ + as isi_distance_impl times, values = isi_distance_impl(spikes1, spikes2) return PieceWiseConstFunc(times, values) diff --git a/pyspike/python_backend.py b/pyspike/python_backend.py deleted file mode 100644 index 481daf9..0000000 --- a/pyspike/python_backend.py +++ /dev/null @@ -1,485 +0,0 @@ -""" python_backend.py - -Collection of python functions that can be used instead of the cython -implementation. - -Copyright 2014, Mario Mulansky - -Distributed under the BSD License - -""" - -import numpy as np - - -############################################################ -# isi_distance_python -############################################################ -def isi_distance_python(s1, s2): - """ Plain Python implementation of the isi distance. - """ - # compute the interspike interval - nu1 = s1[1:] - s1[:-1] - nu2 = s2[1:] - s2[:-1] - - # compute the isi-distance - spike_events = np.empty(len(nu1) + len(nu2)) - spike_events[0] = s1[0] - # the values have one entry less - the number of intervals between events - isi_values = np.empty(len(spike_events) - 1) - # add the distance of the first events - # isi_values[0] = nu1[0]/nu2[0] - 1.0 if nu1[0] <= nu2[0] \ - # else 1.0 - nu2[0]/nu1[0] - isi_values[0] = abs(nu1[0] - nu2[0]) / max(nu1[0], nu2[0]) - index1 = 0 - index2 = 0 - index = 1 - while True: - # check which spike is next - from s1 or s2 - if s1[index1+1] < s2[index2+1]: - index1 += 1 - # break condition relies on existence of spikes at T_end - if index1 >= len(nu1): - break - spike_events[index] = s1[index1] - elif s1[index1+1] > s2[index2+1]: - index2 += 1 - if index2 >= len(nu2): - break - spike_events[index] = s2[index2] - else: # s1[index1 + 1] == s2[index2 + 1] - index1 += 1 - index2 += 1 - if (index1 >= len(nu1)) or (index2 >= len(nu2)): - break - spike_events[index] = s1[index1] - # compute the corresponding isi-distance - isi_values[index] = abs(nu1[index1] - nu2[index2]) / \ - max(nu1[index1], nu2[index2]) - index += 1 - # the last event is the interval end - spike_events[index] = s1[-1] - # use only the data added above - # could be less than original length due to equal spike times - return spike_events[:index + 1], isi_values[:index] - - -############################################################ -# get_min_dist -############################################################ -def get_min_dist(spike_time, spike_train, start_index=0): - """ Returns the minimal distance |spike_time - spike_train[i]| - with i>=start_index. - """ - d = abs(spike_time - spike_train[start_index]) - start_index += 1 - while start_index < len(spike_train): - d_temp = abs(spike_time - spike_train[start_index]) - if d_temp > d: - break - else: - d = d_temp - start_index += 1 - return d - - -############################################################ -# spike_distance_python -############################################################ -def spike_distance_python(spikes1, spikes2): - """ Computes the instantaneous spike-distance S_spike (t) of the two given - spike trains. 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. - Args: - - spikes1, spikes2: ordered arrays of spike times with auxiliary spikes. - Returns: - - PieceWiseLinFunc describing the spike-distance. - """ - # 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!" - # shorter variables - t1 = spikes1 - t2 = spikes2 - - spike_events = np.empty(len(t1) + len(t2) - 2) - spike_events[0] = t1[0] - y_starts = np.empty(len(spike_events) - 1) - y_ends = np.empty(len(spike_events) - 1) - - index1 = 0 - index2 = 0 - index = 1 - dt_p1 = 0.0 - dt_f1 = get_min_dist(t1[1], t2, 0) - dt_p2 = 0.0 - dt_f2 = get_min_dist(t2[1], t1, 0) - isi1 = max(t1[1]-t1[0], t1[2]-t1[1]) - isi2 = max(t2[1]-t2[0], t2[2]-t2[1]) - s1 = dt_f1*(t1[1]-t1[0])/isi1 - s2 = dt_f2*(t2[1]-t2[0])/isi2 - y_starts[0] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2) - while True: - # print(index, index1, index2) - if t1[index1+1] < t2[index2+1]: - index1 += 1 - # break condition relies on existence of spikes at T_end - if index1+1 >= len(t1): - break - spike_events[index] = t1[index1] - # first calculate the previous interval end value - dt_p1 = dt_f1 # the previous time was the following time before - s1 = dt_p1 - s2 = (dt_p2*(t2[index2+1]-t1[index1]) + - dt_f2*(t1[index1]-t2[index2])) / isi2 - y_ends[index-1] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2) - # now the next interval start value - dt_f1 = get_min_dist(t1[index1+1], t2, index2) - isi1 = t1[index1+1]-t1[index1] - # s2 is the same as above, thus we can compute y2 immediately - y_starts[index] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2) - elif t1[index1+1] > t2[index2+1]: - index2 += 1 - if index2+1 >= len(t2): - break - spike_events[index] = t2[index2] - # first calculate the previous interval end value - dt_p2 = dt_f2 # the previous time was the following time before - s1 = (dt_p1*(t1[index1+1]-t2[index2]) + - dt_f1*(t2[index2]-t1[index1])) / isi1 - s2 = dt_p2 - y_ends[index-1] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2) - # now the next interval start value - dt_f2 = get_min_dist(t2[index2+1], t1, index1) - #s2 = dt_f2 - isi2 = t2[index2+1]-t2[index2] - # s2 is the same as above, thus we can compute y2 immediately - y_starts[index] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2) - else: # t1[index1+1] == t2[index2+1] - generate only one event - index1 += 1 - index2 += 1 - if (index1+1 >= len(t1)) or (index2+1 >= len(t2)): - break - assert dt_f2 == 0.0 - assert dt_f1 == 0.0 - spike_events[index] = t1[index1] - y_ends[index-1] = 0.0 - y_starts[index] = 0.0 - dt_p1 = 0.0 - dt_p2 = 0.0 - dt_f1 = get_min_dist(t1[index1+1], t2, index2) - dt_f2 = get_min_dist(t2[index2+1], t1, index1) - isi1 = t1[index1+1]-t1[index1] - isi2 = t2[index2+1]-t2[index2] - index += 1 - # the last event is the interval end - spike_events[index] = t1[-1] - # the ending value of the last interval - isi1 = max(t1[-1]-t1[-2], t1[-2]-t1[-3]) - isi2 = max(t2[-1]-t2[-2], t2[-2]-t2[-3]) - s1 = dt_p1*(t1[-1]-t1[-2])/isi1 - s2 = dt_p2*(t2[-1]-t2[-2])/isi2 - y_ends[index-1] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2) - # use only the data added above - # could be less than original length due to equal spike times - return spike_events[:index+1], y_starts[:index], y_ends[:index] - - -############################################################ -# cumulative_sync_python -############################################################ -def cumulative_sync_python(spikes1, spikes2): - - def get_tau(spikes1, spikes2, i, j): - return 0.5*min([spikes1[i]-spikes1[i-1], spikes1[i+1]-spikes1[i], - spikes2[j]-spikes2[j-1], spikes2[j+1]-spikes2[j]]) - N1 = len(spikes1) - N2 = len(spikes2) - i = 0 - j = 0 - n = 0 - st = np.zeros(N1 + N2 - 2) - c = np.zeros(N1 + N2 - 3) - c[0] = 0 - st[0] = 0 - while n < N1 + N2: - if spikes1[i+1] < spikes2[j+1]: - i += 1 - n += 1 - tau = get_tau(spikes1, spikes2, i, j) - st[n] = spikes1[i] - if spikes1[i]-spikes2[j] > tau: - c[n] = c[n-1] - else: - c[n] = c[n-1]+1 - elif spikes1[i+1] > spikes2[j+1]: - j += 1 - n += 1 - tau = get_tau(spikes1, spikes2, i, j) - st[n] = spikes2[j] - if spikes2[j]-spikes1[i] > tau: - c[n] = c[n-1] - else: - c[n] = c[n-1]+1 - else: # spikes1[i+1] = spikes2[j+1] - j += 1 - i += 1 - if i == N1-1 or j == N2-1: - break - n += 1 - st[n] = spikes1[i] - c[n] = c[n-1] - n += 1 - st[n] = spikes1[i] - c[n] = c[n-1]+1 - c[0] = 0 - st[0] = spikes1[0] - st[-1] = spikes1[-1] - - return st, c - - -############################################################ -# coincidence_python -############################################################ -def coincidence_python(spikes1, spikes2): - - def get_tau(spikes1, spikes2, i, j): - m = 1E100 # some huge number - if i < len(spikes1)-2: - m = min(m, spikes1[i+1]-spikes1[i]) - if j < len(spikes2)-2: - m = min(m, spikes2[j+1]-spikes2[j]) - if i > 1: - m = min(m, spikes1[i]-spikes1[i-1]) - if j > 1: - m = min(m, spikes2[j]-spikes2[j-1]) - return 0.5*m - N1 = len(spikes1) - N2 = len(spikes2) - i = 0 - j = 0 - n = 0 - st = np.zeros(N1 + N2 - 2) # spike times - c = np.zeros(N1 + N2 - 2) # coincidences - mp = np.ones(N1 + N2 - 2) # multiplicity - while n < N1 + N2 - 2: - if spikes1[i+1] < spikes2[j+1]: - i += 1 - n += 1 - tau = get_tau(spikes1, spikes2, i, j) - st[n] = spikes1[i] - if j > 0 and spikes1[i]-spikes2[j] < tau: - # coincidence between the current spike and the previous spike - # both get marked with 1 - c[n] = 1 - c[n-1] = 1 - elif spikes1[i+1] > spikes2[j+1]: - j += 1 - n += 1 - tau = get_tau(spikes1, spikes2, i, j) - st[n] = spikes2[j] - if i > 0 and spikes2[j]-spikes1[i] < tau: - # coincidence between the current spike and the previous spike - # both get marked with 1 - c[n] = 1 - c[n-1] = 1 - else: # spikes1[i+1] = spikes2[j+1] - # advance in both spike trains - j += 1 - i += 1 - if i == N1-1 or j == N2-1: - break - n += 1 - # add only one event, but with coincidence 2 and multiplicity 2 - st[n] = spikes1[i] - c[n] = 2 - mp[n] = 2 - - st = st[:n+2] - c = c[:n+2] - mp = mp[:n+2] - - st[0] = spikes1[0] - st[-1] = spikes1[-1] - c[0] = c[1] - c[-1] = c[-2] - mp[0] = mp[1] - mp[-1] = mp[-2] - - return st, c, mp - - -############################################################ -# add_piece_wise_const_python -############################################################ -def add_piece_wise_const_python(x1, y1, x2, y2): - x_new = np.empty(len(x1) + len(x2)) - y_new = np.empty(len(x_new)-1) - x_new[0] = x1[0] - y_new[0] = y1[0] + y2[0] - index1 = 0 - index2 = 0 - index = 0 - while (index1+1 < len(y1)) and (index2+1 < len(y2)): - index += 1 - # print(index1+1, x1[index1+1], y1[index1+1], x_new[index]) - if x1[index1+1] < x2[index2+1]: - index1 += 1 - x_new[index] = x1[index1] - elif x1[index1+1] > x2[index2+1]: - index2 += 1 - x_new[index] = x2[index2] - else: # x1[index1+1] == x2[index2+1]: - index1 += 1 - index2 += 1 - x_new[index] = x1[index1] - y_new[index] = y1[index1] + y2[index2] - # one array reached the end -> copy the contents of the other to the end - if index1+1 < len(y1): - x_new[index+1:index+1+len(x1)-index1-1] = x1[index1+1:] - y_new[index+1:index+1+len(y1)-index1-1] = y1[index1+1:] + y2[-1] - index += len(x1)-index1-2 - elif index2+1 < len(y2): - x_new[index+1:index+1+len(x2)-index2-1] = x2[index2+1:] - y_new[index+1:index+1+len(y2)-index2-1] = y2[index2+1:] + y1[-1] - index += len(x2)-index2-2 - else: # both arrays reached the end simultaneously - # only the last x-value missing - x_new[index+1] = x1[-1] - # the last value is again the end of the interval - # x_new[index+1] = x1[-1] - # only use the data that was actually filled - - return x_new[:index+2], y_new[:index+1] - - -############################################################ -# add_piece_lin_const_python -############################################################ -def add_piece_wise_lin_python(x1, y11, y12, x2, y21, y22): - x_new = np.empty(len(x1) + len(x2)) - y1_new = np.empty(len(x_new)-1) - y2_new = np.empty_like(y1_new) - x_new[0] = x1[0] - y1_new[0] = y11[0] + y21[0] - index1 = 0 # index for self - index2 = 0 # index for f - index = 0 # index for new - while (index1+1 < len(y11)) and (index2+1 < len(y21)): - # print(index1+1, x1[index1+1], self.y[index1+1], x_new[index]) - if x1[index1+1] < x2[index2+1]: - # first compute the end value of the previous interval - # linear interpolation of the interval - y = y21[index2] + (y22[index2]-y21[index2]) * \ - (x1[index1+1]-x2[index2]) / (x2[index2+1]-x2[index2]) - y2_new[index] = y12[index1] + y - index1 += 1 - index += 1 - x_new[index] = x1[index1] - # and the starting value for the next interval - y1_new[index] = y11[index1] + y - elif x1[index1+1] > x2[index2+1]: - # first compute the end value of the previous interval - # linear interpolation of the interval - y = y11[index1] + (y12[index1]-y11[index1]) * \ - (x2[index2+1]-x1[index1]) / \ - (x1[index1+1]-x1[index1]) - y2_new[index] = y22[index2] + y - index2 += 1 - index += 1 - x_new[index] = x2[index2] - # and the starting value for the next interval - y1_new[index] = y21[index2] + y - else: # x1[index1+1] == x2[index2+1]: - y2_new[index] = y12[index1] + y22[index2] - index1 += 1 - index2 += 1 - index += 1 - x_new[index] = x1[index1] - y1_new[index] = y11[index1] + y21[index2] - # one array reached the end -> copy the contents of the other to the end - if index1+1 < len(y11): - # compute the linear interpolations values - y = y21[index2] + (y22[index2]-y21[index2]) * \ - (x1[index1+1:-1]-x2[index2]) / (x2[index2+1]-x2[index2]) - x_new[index+1:index+1+len(x1)-index1-1] = x1[index1+1:] - y1_new[index+1:index+1+len(y11)-index1-1] = y11[index1+1:]+y - y2_new[index:index+len(y12)-index1-1] = y12[index1:-1] + y - index += len(x1)-index1-2 - elif index2+1 < len(y21): - # compute the linear interpolations values - y = y11[index1] + (y12[index1]-y11[index1]) * \ - (x2[index2+1:-1]-x1[index1]) / \ - (x1[index1+1]-x1[index1]) - x_new[index+1:index+1+len(x2)-index2-1] = x2[index2+1:] - y1_new[index+1:index+1+len(y21)-index2-1] = y21[index2+1:] + y - y2_new[index:index+len(y22)-index2-1] = y22[index2:-1] + y - index += len(x2)-index2-2 - else: # both arrays reached the end simultaneously - # only the last x-value missing - x_new[index+1] = x1[-1] - # finally, the end value for the last interval - y2_new[index] = y12[-1]+y22[-1] - # only use the data that was actually filled - return x_new[:index+2], y1_new[:index+1], y2_new[:index+1] - - -############################################################ -# add_discrete_function_python -############################################################ -def add_discrete_function_python(x1, y1, mp1, x2, y2, mp2): - - x_new = np.empty(len(x1) + len(x2)) - y_new = np.empty_like(x_new) - mp_new = np.empty_like(x_new) - x_new[0] = x1[0] - index1 = 0 - index2 = 0 - index = 0 - while (index1+1 < len(y1)) and (index2+1 < len(y2)): - if x1[index1+1] < x2[index2+1]: - index1 += 1 - index += 1 - x_new[index] = x1[index1] - y_new[index] = y1[index1] - mp_new[index] = mp1[index1] - elif x1[index1+1] > x2[index2+1]: - index2 += 1 - index += 1 - x_new[index] = x2[index2] - y_new[index] = y2[index2] - mp_new[index] = mp2[index2] - else: # x1[index1+1] == x2[index2+1] - index1 += 1 - index2 += 1 - index += 1 - x_new[index] = x1[index1] - y_new[index] = y1[index1] + y2[index2] - mp_new[index] = mp1[index1] + mp2[index2] - # one array reached the end -> copy the contents of the other to the end - if index1+1 < len(y1): - x_new[index+1:index+1+len(x1)-index1-1] = x1[index1+1:] - y_new[index+1:index+1+len(x1)-index1-1] = y1[index1+1:] - mp_new[index+1:index+1+len(x1)-index1-1] = mp1[index1+1:] - index += len(x1)-index1-1 - elif index2+1 < len(y2): - x_new[index+1:index+1+len(x2)-index2-1] = x2[index2+1:] - y_new[index+1:index+1+len(x2)-index2-1] = y2[index2+1:] - mp_new[index+1:index+1+len(x2)-index2-1] = mp2[index2+1:] - index += len(x2)-index2-1 - # else: # both arrays reached the end simultaneously - # x_new[index+1] = x1[-1] - # y_new[index+1] = y1[-1] + y2[-1] - # mp_new[index+1] = mp1[-1] + mp2[-1] - - y_new[0] = y_new[1] - mp_new[0] = mp_new[1] - - # the last value is again the end of the interval - # only use the data that was actually filled - return x_new[:index+1], y_new[:index+1], mp_new[:index+1] - diff --git a/pyspike/spike_distance.py b/pyspike/spike_distance.py index 2c989a4..f721c86 100644 --- a/pyspike/spike_distance.py +++ b/pyspike/spike_distance.py @@ -35,14 +35,15 @@ def spike_profile(spikes1, spikes2): # cython implementation try: - from cython_distance import spike_distance_cython \ + 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 python_backend import spike_distance_python as spike_distance_impl + 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) diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py index bded8da..342bf69 100644 --- a/pyspike/spike_sync.py +++ b/pyspike/spike_sync.py @@ -33,14 +33,14 @@ def spike_sync_profile(spikes1, spikes2): # cython implementation try: - from cython_distance import coincidence_cython \ + 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 python_backend import coincidence_python \ + from cython.python_backend import coincidence_python \ as coincidence_impl times, coincidences, multiplicity = coincidence_impl(spikes1, spikes2) diff --git a/setup.py b/setup.py index 9cab5de..289d521 100644 --- a/setup.py +++ b/setup.py @@ -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'] } ) -- cgit v1.2.3 From bfee34894f21212bcef5fe271e701b1cf8c86fde Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Tue, 3 Feb 2015 12:26:02 +0100 Subject: +init py --- pyspike/cython/__init__.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 pyspike/cython/__init__.py diff --git a/pyspike/cython/__init__.py b/pyspike/cython/__init__.py new file mode 100644 index 0000000..e69de29 -- cgit v1.2.3