diff options
Diffstat (limited to 'pyspike')
| -rw-r--r-- | pyspike/__init__.py | 7 | ||||
| -rw-r--r-- | pyspike/cython_add.pyx | 59 | ||||
| -rw-r--r-- | pyspike/cython_distance.pyx | 81 | ||||
| -rw-r--r-- | pyspike/distances.py | 146 | ||||
| -rw-r--r-- | pyspike/function.py | 165 | ||||
| -rw-r--r-- | pyspike/python_backend.py | 112 | ||||
| -rw-r--r-- | pyspike/spikes.py | 4 |
7 files changed, 502 insertions, 72 deletions
diff --git a/pyspike/__init__.py b/pyspike/__init__.py index dca5722..55687e6 100644 --- a/pyspike/__init__.py +++ b/pyspike/__init__.py @@ -6,12 +6,13 @@ Distributed under the BSD License __all__ = ["function", "distances", "spikes"] -from function import PieceWiseConstFunc, PieceWiseLinFunc, average_profile +from function import PieceWiseConstFunc, PieceWiseLinFunc, \ + DiscreteFunction, average_profile from distances import isi_profile, isi_distance, \ spike_profile, spike_distance, \ - spike_sync_profile, \ + spike_sync_profile, spike_sync, \ isi_profile_multi, isi_distance_multi, isi_distance_matrix, \ spike_profile_multi, spike_distance_multi, spike_distance_matrix, \ - spike_sync_profile_multi + spike_sync_profile_multi, spike_sync_multi, spike_sync_matrix from spikes import add_auxiliary_spikes, load_spike_trains_from_txt, \ spike_train_from_string, merge_spike_trains, generate_poisson_spikes diff --git a/pyspike/cython_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/cython_distance.pyx b/pyspike/cython_distance.pyx index 779ff94..489aab9 100644 --- a/pyspike/cython_distance.pyx +++ b/pyspike/cython_distance.pyx @@ -33,6 +33,7 @@ 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 @@ -229,3 +230,83 @@ def spike_distance_cython(double[:] t1, # 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/distances.py b/pyspike/distances.py index 04ea77b..5476b6f 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 +from pyspike import PieceWiseConstFunc, PieceWiseLinFunc, DiscreteFunction ############################################################ @@ -132,39 +132,57 @@ def spike_distance(spikes1, spikes2, interval=None): ############################################################ # spike_sync_profile ############################################################ -def spike_sync_profile(spikes1, spikes2, k=3): +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. - assert k > 0 + :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.") + 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 - st, c = coincidence_impl(spikes1, spikes2) + times, coincidences, multiplicity = coincidence_impl(spikes1, spikes2) - dc = np.zeros(len(c)) - for i in xrange(2*k): - dc[k:-k] += c[i:-2*k+i] + return DiscreteFunction(times, coincidences, multiplicity) - for n in xrange(0, k): - for i in xrange(n+k): - dc[n] += c[i] - dc[-n-1] += c[-i-1] - for i in xrange(k-n-1): - dc[n] += c[i] - dc[-n-1] += c[-i-1] - dc *= 1.0/k +############################################################ +# 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: - return PieceWiseConstFunc(st, dc) + .. 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) ############################################################ @@ -201,8 +219,7 @@ def _generic_profile_multi(spike_trains, pair_distance_func, indices=None): 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 - average_dist.mul_scalar(1.0/len(pairs)) # normalize - return average_dist + return average_dist, len(pairs) ############################################################ @@ -273,7 +290,10 @@ def isi_profile_multi(spike_trains, indices=None): :returns: The averaged isi profile :math:`<S_{isi}>(t)` :rtype: :class:`pyspike.function.PieceWiseConstFunc` """ - return _generic_profile_multi(spike_trains, isi_profile, indices) + average_dist, M = _generic_profile_multi(spike_trains, isi_profile, + indices) + average_dist.mul_scalar(1.0/M) # normalize + return average_dist ############################################################ @@ -314,39 +334,65 @@ def spike_profile_multi(spike_trains, indices=None): :rtype: :class:`pyspike.function.PieceWiseLinFunc` """ - return _generic_profile_multi(spike_trains, spike_profile, indices) + average_dist, M = _generic_profile_multi(spike_trains, spike_profile, + indices) + average_dist.mul_scalar(1.0/M) # normalize + return average_dist + + +############################################################ +# spike_distance_multi +############################################################ +def spike_distance_multi(spike_trains, indices=None, interval=None): + """ Computes the multi-variate spike distance for a set of spike trains. + That is the time average of the multi-variate spike profile: + S_{spike} = \int_0^T 2/((N(N-1)) sum_{<i,j>} S_{spike}^{i, j} dt + where the sum goes over all pairs <i,j> + + :param spike_trains: list of spike trains + :param indices: list of indices defining which spike trains to use, + if None all given spike trains are used (default=None) + :type indices: list or None + :param interval: averaging interval given as a pair of floats, if None + the average over the whole function is computed. + :type interval: Pair of floats or None. + :returns: The averaged spike distance S. + :rtype: double + """ + return spike_profile_multi(spike_trains, indices).avrg(interval) ############################################################ # spike_profile_multi ############################################################ -def spike_sync_profile_multi(spike_trains, indices=None, k=3): +def spike_sync_profile_multi(spike_trains, indices=None): """ Computes the multi-variate spike synchronization profile for a set of - spike trains. That is the average spike-distance of all pairs of spike - trains: - :math:`S_ss(t) = 2/((N(N-1)) sum_{<i,j>} S_{ss}^{i, j}`, - where the sum goes over all pairs <i,j> + 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 averaged spike profile :math:`<S_{ss}>(t)` - :rtype: :class:`pyspike.function.PieceWiseConstFunc` + :returns: The multi-variate spike sync profile :math:`<S_{sync}>(t)` + :rtype: :class:`pyspike.function.DiscreteFunction` """ - prof_func = partial(spike_sync_profile, k=k) - return _generic_profile_multi(spike_trains, prof_func, indices) + 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_distance_multi(spike_trains, indices=None, interval=None): - """ Computes the multi-variate spike distance for a set of spike trains. - That is the time average of the multi-variate spike profile: - S_{spike} = \int_0^T 2/((N(N-1)) sum_{<i,j>} S_{spike}^{i, j} dt - where the sum goes over all pairs <i,j> +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, @@ -355,10 +401,10 @@ def spike_distance_multi(spike_trains, indices=None, interval=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. + :returns: The multi-variate spike synchronization value SYNC. :rtype: double """ - return spike_profile_multi(spike_trains, indices).avrg(interval) + return spike_sync_profile_multi(spike_trains, indices).avrg(interval) ############################################################ @@ -434,3 +480,25 @@ def spike_distance_matrix(spike_trains, indices=None, interval=None): """ return _generic_distance_matrix(spike_trains, spike_distance, indices, interval) + + +############################################################ +# spike_sync_matrix +############################################################ +def spike_sync_matrix(spike_trains, indices=None, interval=None): + """ Computes the overall spike-synchronization value of all pairs of + spike-trains. + + :param spike_trains: list of spike trains + :param indices: list of indices defining which spike trains to use, + if None all given spike trains are used (default=None) + :type indices: list or None + :param interval: averaging interval given as a pair of floats, if None + the average over the whole function is computed. + :type interval: Pair of floats or None. + :returns: 2D array with the pair wise time spike synchronization values + :math:`SYNC_{ij}` + :rtype: np.array + """ + return _generic_distance_matrix(spike_trains, spike_sync, + indices, interval) diff --git a/pyspike/function.py b/pyspike/function.py index 662606c..e0dadf6 100644 --- a/pyspike/function.py +++ b/pyspike/function.py @@ -136,15 +136,6 @@ class PieceWiseConstFunc(object): a /= int_length return a - def avrg_function_value(self): - """ Computes the average function value of the piece-wise const - function: :math:`a = 1/N sum_i f_i` where N is the number of intervals. - - :returns: the average a. - :rtype: float - """ - return sum(self.y)/(len(self.y)) - def add(self, f): """ Adds another PieceWiseConst function to this function. Note: only functions defined on the same interval can be summed. @@ -364,6 +355,162 @@ 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 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. diff --git a/pyspike/python_backend.py b/pyspike/python_backend.py index 7f8ea8c..481daf9 100644 --- a/pyspike/python_backend.py +++ b/pyspike/python_backend.py @@ -248,52 +248,69 @@ def cumulative_sync_python(spikes1, spikes2): def coincidence_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]]) + 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) - c = np.zeros(N1 + N2 - 3) - c[0] = 0 - st[0] = 0 - while n < N1 + N2: + 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 spikes1[i]-spikes2[j] > tau: - c[n] = 0 - else: + 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 spikes2[j]-spikes1[i] > tau: - c[n] = 0 - else: + 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] = 0 - n += 1 - st[n] = spikes1[i] - c[n] = 1 - c[0] = c[2] + 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 + return st, c, mp ############################################################ @@ -409,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] + diff --git a/pyspike/spikes.py b/pyspike/spikes.py index aa25c48..6a3353e 100644 --- a/pyspike/spikes.py +++ b/pyspike/spikes.py @@ -67,7 +67,7 @@ def spike_train_from_string(s, sep=' ', is_sorted=False): # load_spike_trains_txt ############################################################ def load_spike_trains_from_txt(file_name, time_interval=None, - separator=' ', comment='#', sort=True): + separator=' ', comment='#', is_sorted=False): """ Loads a number of spike trains from a text file. Each line of the text file should contain one spike train as a sequence of spike times separated by `separator`. Empty lines as well as lines starting with `comment` are @@ -94,7 +94,7 @@ def load_spike_trains_from_txt(file_name, time_interval=None, for line in spike_file: if len(line) > 1 and not line.startswith(comment): # use only the lines with actual data and not commented - spike_train = spike_train_from_string(line, separator, sort) + spike_train = spike_train_from_string(line, separator, is_sorted) if time_interval is not None: # add auxil. spikes if times given spike_train = add_auxiliary_spikes(spike_train, time_interval) spike_trains.append(spike_train) |