From fed0ceec753fc1a7e5a1e20632de5a9800fe4fb1 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Mon, 19 Jan 2015 16:39:17 +0100 Subject: final version for spike sync --- pyspike/__init__.py | 4 +- pyspike/distances.py | 36 +++---------- pyspike/function.py | 58 +++++++++----------- pyspike/python_backend.py | 135 ++++++++++++++++++++++------------------------ pyspike/spikes.py | 4 +- 5 files changed, 102 insertions(+), 135 deletions(-) (limited to 'pyspike') diff --git a/pyspike/__init__.py b/pyspike/__init__.py index fa90d99..74d52c5 100644 --- a/pyspike/__init__.py +++ b/pyspike/__init__.py @@ -6,8 +6,8 @@ Distributed under the BSD License __all__ = ["function", "distances", "spikes"] -from function import PieceWiseConstFunc, PieceWiseLinFunc, IntervalSequence,\ - average_profile +from function import PieceWiseConstFunc, PieceWiseLinFunc, \ + MultipleValueSequence, average_profile from distances import isi_profile, isi_distance, \ spike_profile, spike_distance, \ spike_sync_profile, spike_sync_distance, \ diff --git a/pyspike/distances.py b/pyspike/distances.py index 38c5cc2..5ee8261 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, IntervalSequence +from pyspike import PieceWiseConstFunc, PieceWiseLinFunc, MultipleValueSequence ############################################################ @@ -132,9 +132,7 @@ def spike_distance(spikes1, spikes2, interval=None): ############################################################ # spike_sync_profile ############################################################ -def spike_sync_profile(spikes1, spikes2, k=3): - - assert k > 0 +def spike_sync_profile(spikes1, spikes2): # cython implementation try: @@ -148,34 +146,16 @@ def spike_sync_profile(spikes1, spikes2, k=3): from python_backend import coincidence_python \ as coincidence_impl - st, J = coincidence_impl(spikes1, spikes2) - - N = len(J) - - # compute the cumulative sum, include some extra values for boundary - # conditions - c = np.zeros(N + 2*k) - c[k:-k] = np.cumsum(J) - # set the boundary values - # on the left: c_0 = -c_1, c_{-1} = -c_2, ..., c{-k+1} = c_k - # on the right: c_{N+1} = c_N, c_{N+2} = 2*c_N - c_{N-1}, - # c_{N+2} = 2*c_N - c_{N-2}, ..., c_{N+k} = 2*c_N - c_{N-k+1} - for n in xrange(k): - c[k-n-1] = -c[k+n] - c[-k+n] = 2*c[-k-1] - c[-k-1-n] - # with the right boundary values, the differences become trivial - J_w = c[2*k:] - c[:-2*k] - # normalize to half the interval width - J_w *= 1.0/k + times, coincidences, multiplicity = coincidence_impl(spikes1, spikes2) - return IntervalSequence(st, J_w) + return MultipleValueSequence(times, coincidences, multiplicity) ############################################################ # spike_sync_distance ############################################################ -def spike_sync_distance(spikes1, spikes2, k=3): - return spike_sync_profile(spikes1, spikes2, k).avrg() +def spike_sync_distance(spikes1, spikes2): + return spike_sync_profile(spikes1, spikes2).avrg() ############################################################ @@ -336,7 +316,7 @@ def spike_profile_multi(spike_trains, indices=None): ############################################################ # 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: @@ -351,7 +331,7 @@ def spike_sync_profile_multi(spike_trains, indices=None, k=3): :rtype: :class:`pyspike.function.PieceWiseConstFunc` """ - prof_func = partial(spike_sync_profile, k=k) + 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! diff --git a/pyspike/function.py b/pyspike/function.py index f5a1133..f10c136 100644 --- a/pyspike/function.py +++ b/pyspike/function.py @@ -171,32 +171,32 @@ that PySpike is installed by running\n 'python setup.py build_ext --inplace'! \ ############################################################## -# IntervalSequence +# MultipleValueSequence ############################################################## -class IntervalSequence(object): +class MultipleValueSequence(object): """ A class representing a sequence of values defined in some interval. - This is very similar to a `PieceWiseConstFunc`, but with different - averaging and addition. """ - def __init__(self, x, y): - """ Constructs the interval sequence. + def __init__(self, x, y, multiplicity): + """ Constructs the value sequence. - :param x: array of length N+1 defining the edges of the intervals of - the intervals. - :param y: array of length N defining the values at the intervals. + :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.extra_zero_intervals = 0 + self.mp = np.array(multiplicity) def copy(self): """ Returns a copy of itself :rtype: :class:`IntervalSequence` """ - return IntervalSequence(self.x, self.y) + 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` @@ -209,9 +209,10 @@ class IntervalSequence(object): """ 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) + 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): + def get_plottable_data(self, k=0): """ Returns two arrays containing x- and y-coordinates for immeditate plotting of the interval sequence. @@ -224,17 +225,10 @@ class IntervalSequence(object): 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 - normalization = 1.0 * (len(self.y)-1) / (len(self.y) + - self.extra_zero_intervals-1) - y_plot[1::2] = self.y + if k > 0: + raise NotImplementedError() - return x_plot, y_plot * normalization + 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 @@ -250,7 +244,7 @@ class IntervalSequence(object): if interval is None: # no interval given, integrate over the whole spike train # don't count the first value, which is zero by definition - a = np.sum(self.y) + a = 1.0*np.sum(self.y[1:-1]) else: raise NotImplementedError() return a @@ -270,15 +264,15 @@ class IntervalSequence(object): if interval is None: # no interval given, average over the whole spike train # don't count the first interval for normalization - return self.integral() / (len(self.y)-1+self.extra_zero_intervals) + return self.integral() / np.sum(self.mp[1:-1]) else: raise NotImplementedError() def add(self, f): - """ Adds another `IntervalSequence` function to this function. + """ Adds another `MultipleValueSequence` function to this function. Note: only functions defined on the same interval can be summed. - :param f: :class:`IntervalSequence` function to be added. + :param f: :class:`MultipleValueSequence` function to be added. :rtype: None """ assert self.x[0] == f.x[0], "The functions have different intervals" @@ -293,12 +287,12 @@ class IntervalSequence(object): # 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_interval_sequence_python as \ - add_interval_sequence_impl + from python_backend import add_multiple_value_sequence_python as \ + add_multiple_value_sequence_impl - self.x, self.y, extra_intervals = \ - add_interval_sequence_impl(self.x, self.y, f.x, f.y) - self.extra_zero_intervals += extra_intervals + 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 diff --git a/pyspike/python_backend.py b/pyspike/python_backend.py index 154d250..bbbd572 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 ############################################################ @@ -341,83 +358,59 @@ def add_piece_wise_const_python(x1, y1, x2, y2): ############################################################ -# add_interval_sequence_python +# add_multiple_value_sequence_python ############################################################ -def add_interval_sequence_python(x1, y1, x2, y2): - yscale1 = np.empty_like(y1) - index2 = 1 - # s1 = (len(y1)+len(y2)-2.0) / (len(y1)-1.0) - # s2 = (len(y1)+len(y2)-2.0) / (len(y2)-1.0) - s1 = 1.0 - s2 = 1.0 - for i in xrange(len(y1)): - c = 1 - while index2 < len(x2)-1 and x2[index2] < x1[i+1]: - index2 += 1 - c += 1 - if index2 < len(x2)-1 and x2[index2] == x1[i+1]: - index2 += 1 - # c += 1 - yscale1[i] = s1/c - - yscale2 = np.empty_like(y2) - index1 = 1 - for i in xrange(len(y2)): - c = 1 - while index1 < len(x1)-1 and x1[index1] < x2[i+1]: - index1 += 1 - c += 1 - if index1 < len(x1)-1 and x1[index1] == x2[i+1]: - index1 += 1 - # c += 1 - yscale2[i] = s2/c +def add_multiple_value_sequence_python(x1, y1, mp1, x2, y2, mp2): x_new = np.empty(len(x1) + len(x2)) - y_new = np.empty(len(x_new)-1) + 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 - additional_intervals = 0 while (index1+1 < len(y1)) and (index2+1 < len(y2)): - y_new[index] = y1[index1]*yscale1[index1] + y2[index2]*yscale2[index2] - index += 1 - # print(index1+1, x1[index1+1], y1[index1+1], x_new[index]) 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] - # y_new[index] = y1[index1]*yscale1[index1] + \ - # y2[index2]*yscale2[index2] index1 += 1 - # x_new[index] = x1[index1] index2 += 1 - # index += 1 + index += 1 x_new[index] = x1[index1] - additional_intervals += 1 - y_new[index] = y1[index1]*yscale1[index1] + y2[index2]*yscale2[index2] + 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(y1)-index1-1] = \ - y1[index1+1:]*yscale1[index1+1:] + y2[-1]*yscale2[-1] - index += len(x1)-index1-2 + 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(y2)-index2-1] = \ - y2[index2+1:]*yscale2[index2+1:] + y1[-1]*yscale1[-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] + 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 - # x_new[index+1] = x1[-1] # only use the data that was actually filled - - return x_new[:index+2], y_new[:index+1], additional_intervals + 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) -- cgit v1.2.3