diff options
-rw-r--r-- | pyspike/__init__.py | 2 | ||||
-rw-r--r-- | pyspike/cython_add.pyx | 59 | ||||
-rw-r--r-- | pyspike/distances.py | 6 | ||||
-rw-r--r-- | pyspike/function.py | 269 | ||||
-rw-r--r-- | pyspike/python_backend.py | 113 |
5 files changed, 255 insertions, 194 deletions
diff --git a/pyspike/__init__.py b/pyspike/__init__.py index 74d52c5..1b18569 100644 --- a/pyspike/__init__.py +++ b/pyspike/__init__.py @@ -7,7 +7,7 @@ Distributed under the BSD License __all__ = ["function", "distances", "spikes"] from function import PieceWiseConstFunc, PieceWiseLinFunc, \ - MultipleValueSequence, average_profile + DiscreteFunction, average_profile from distances import isi_profile, isi_distance, \ spike_profile, spike_distance, \ spike_sync_profile, spike_sync_distance, \ diff --git a/pyspike/cython_add.pyx b/pyspike/cython_add.pyx index bfbe208..817799e 100644 --- a/pyspike/cython_add.pyx +++ b/pyspike/cython_add.pyx @@ -172,3 +172,62 @@ def add_piece_wise_lin_cython(double[:] x1, double[:] y11, double[:] y12, 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 x_new[:index+1], y_new[:index+1], mp_new[:index+1] diff --git a/pyspike/distances.py b/pyspike/distances.py index 8bde724..0f0efa9 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, MultipleValueSequence +from pyspike import PieceWiseConstFunc, PieceWiseLinFunc, DiscreteFunction ############################################################ @@ -148,7 +148,7 @@ Falling back to slow python backend.") times, coincidences, multiplicity = coincidence_impl(spikes1, spikes2) - return MultipleValueSequence(times, coincidences, multiplicity) + return DiscreteFunction(times, coincidences, multiplicity) ############################################################ @@ -328,7 +328,7 @@ def spike_sync_profile_multi(spike_trains, indices=None): if None all given spike trains are used (default=None) :type indices: list or None :returns: The averaged spike profile :math:`<S_{ss}>(t)` - :rtype: :class:`pyspike.function.PieceWiseConstFunc` + :rtype: :class:`pyspike.function.DiscreteFunction` """ prof_func = partial(spike_sync_profile) diff --git a/pyspike/function.py b/pyspike/function.py index f10c136..6fb7537 100644 --- a/pyspike/function.py +++ b/pyspike/function.py @@ -171,140 +171,6 @@ that PySpike is installed by running\n 'python setup.py build_ext --inplace'! \ ############################################################## -# MultipleValueSequence -############################################################## -class MultipleValueSequence(object): - """ A class representing a sequence of values defined in some interval. - """ - - def __init__(self, x, y, multiplicity): - """ Constructs the value sequence. - - :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:`IntervalSequence` - """ - return MultipleValueSequence(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:`IntervalSequence` - :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, k=0): - """ Returns two arrays containing x- and y-coordinates for immeditate - plotting of the interval sequence. - - :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") - """ - - if k > 0: - raise NotImplementedError() - - 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 interval - sequence this amounts to the sum over all values divided by the count - of intervals. - - :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 - # don't count the first value, which is zero by definition - a = 1.0*np.sum(self.y[1:-1]) - else: - raise NotImplementedError() - return a - - 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 - """ - if interval is None: - # no interval given, average over the whole spike train - # don't count the first interval for normalization - return self.integral() / np.sum(self.mp[1:-1]) - else: - raise NotImplementedError() - - def add(self, f): - """ Adds another `MultipleValueSequence` function to this function. - Note: only functions defined on the same interval can be summed. - - :param f: :class:`MultipleValueSequence` 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_interval_sequence_cython as \ - add_interval_sequence_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_multiple_value_sequence_python as \ - add_multiple_value_sequence_impl - - self.x, self.y, self.mp = \ - add_multiple_value_sequence_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 - - -############################################################## # PieceWiseLinFunc ############################################################## class PieceWiseLinFunc: @@ -489,6 +355,141 @@ that PySpike is installed by running\n 'python setup.py build_ext --inplace'! \ 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, k=0): + """ Returns two arrays containing x- and y-coordinates for immeditate + plotting of the interval sequence. + + :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 k > 0: + raise NotImplementedError() + + 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 interval + sequence this amounts to the sum over all values divided by the count + of intervals. + + :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 + # don't count the first value, which is zero by definition + a = 1.0*np.sum(self.y[1:-1]) + else: + raise NotImplementedError() + return a + + 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 + """ + if interval is None: + # no interval given, average over the whole spike train + # don't count the first interval for normalization + return self.integral() / np.sum(self.mp[1:-1]) + else: + raise NotImplementedError() + + 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. diff --git a/pyspike/python_backend.py b/pyspike/python_backend.py index bbbd572..481daf9 100644 --- a/pyspike/python_backend.py +++ b/pyspike/python_backend.py @@ -358,62 +358,6 @@ def add_piece_wise_const_python(x1, y1, x2, y2): ############################################################ -# add_multiple_value_sequence_python -############################################################ -def add_multiple_value_sequence_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] - - -############################################################ # add_piece_lin_const_python ############################################################ def add_piece_wise_lin_python(x1, y11, y12, x2, y21, y22): @@ -482,3 +426,60 @@ def add_piece_wise_lin_python(x1, y11, y12, x2, y21, y22): 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] + |