From 8c98f0043fa785b8352b3c685615da24b30e6149 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Fri, 26 Dec 2014 15:31:15 -0600 Subject: spike sync --- pyspike/__init__.py | 5 +- pyspike/distances.py | 57 ++++++++++++------ pyspike/function.py | 149 +++++++++++++++++++++++++++++++++++++++++++--- pyspike/python_backend.py | 82 ++++++++++++++++++++++++- test/test_distance.py | 54 +++++++++++++++++ 5 files changed, 316 insertions(+), 31 deletions(-) diff --git a/pyspike/__init__.py b/pyspike/__init__.py index dca5722..fa90d99 100644 --- a/pyspike/__init__.py +++ b/pyspike/__init__.py @@ -6,10 +6,11 @@ Distributed under the BSD License __all__ = ["function", "distances", "spikes"] -from function import PieceWiseConstFunc, PieceWiseLinFunc, average_profile +from function import PieceWiseConstFunc, PieceWiseLinFunc, IntervalSequence,\ + average_profile from distances import isi_profile, isi_distance, \ spike_profile, spike_distance, \ - spike_sync_profile, \ + spike_sync_profile, spike_sync_distance, \ isi_profile_multi, isi_distance_multi, isi_distance_matrix, \ spike_profile_multi, spike_distance_multi, spike_distance_matrix, \ spike_sync_profile_multi diff --git a/pyspike/distances.py b/pyspike/distances.py index c28fd7a..38c5cc2 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, IntervalSequence ############################################################ @@ -148,23 +148,34 @@ def spike_sync_profile(spikes1, spikes2, k=3): from python_backend import coincidence_python \ as coincidence_impl - st, c = coincidence_impl(spikes1, spikes2) + st, J = coincidence_impl(spikes1, spikes2) - dc = np.zeros(len(c)) - for i in xrange(2*k): - dc[k:-k] += c[i:-2*k+i] + N = len(J) - 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] + # 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 - dc *= 1.0/k + return IntervalSequence(st, J_w) - return PieceWiseConstFunc(st, dc) + +############################################################ +# spike_sync_distance +############################################################ +def spike_sync_distance(spikes1, spikes2, k=3): + return spike_sync_profile(spikes1, spikes2, k).avrg() ############################################################ @@ -201,8 +212,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 +283,10 @@ def isi_profile_multi(spike_trains, indices=None): :returns: The averaged isi profile :math:`(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,7 +327,10 @@ 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 ############################################################ @@ -336,7 +352,10 @@ def spike_sync_profile_multi(spike_trains, indices=None, k=3): """ prof_func = partial(spike_sync_profile, k=k) - return _generic_profile_multi(spike_trains, prof_func, indices) + average_dist, M = _generic_profile_multi(spike_trains, prof_func, + indices) + # average_dist.mul_scalar(1.0/M) # no normalization here! + return average_dist ############################################################ diff --git a/pyspike/function.py b/pyspike/function.py index 662606c..f5a1133 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. @@ -179,6 +170,146 @@ that PySpike is installed by running\n 'python setup.py build_ext --inplace'! \ self.y *= fac +############################################################## +# IntervalSequence +############################################################## +class IntervalSequence(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. + + :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. + """ + # convert parameters to arrays, also ensures copying + self.x = np.array(x) + self.y = np.array(y) + self.extra_zero_intervals = 0 + + def copy(self): + """ Returns a copy of itself + + :rtype: :class:`IntervalSequence` + """ + return IntervalSequence(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:`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) + + def get_plottable_data(self): + """ 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") + """ + + 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 + + return x_plot, y_plot * normalization + + 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 = np.sum(self.y) + 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() / (len(self.y)-1+self.extra_zero_intervals) + else: + raise NotImplementedError() + + def add(self, f): + """ Adds another `IntervalSequence` function to this function. + Note: only functions defined on the same interval can be summed. + + :param f: :class:`IntervalSequence` 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_interval_sequence_python as \ + add_interval_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 + + 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 ############################################################## diff --git a/pyspike/python_backend.py b/pyspike/python_backend.py index 7f8ea8c..154d250 100644 --- a/pyspike/python_backend.py +++ b/pyspike/python_backend.py @@ -289,7 +289,7 @@ def coincidence_python(spikes1, spikes2): n += 1 st[n] = spikes1[i] c[n] = 1 - c[0] = c[2] + #c[0] = c[2] st[0] = spikes1[0] st[-1] = spikes1[-1] @@ -340,6 +340,86 @@ def add_piece_wise_const_python(x1, y1, x2, y2): return x_new[:index+2], y_new[:index+1] +############################################################ +# add_interval_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 + + x_new = np.empty(len(x1) + len(x2)) + y_new = np.empty(len(x_new)-1) + 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 + 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] + # y_new[index] = y1[index1]*yscale1[index1] + \ + # y2[index2]*yscale2[index2] + index1 += 1 + # x_new[index] = x1[index1] + index2 += 1 + # index += 1 + x_new[index] = x1[index1] + additional_intervals += 1 + y_new[index] = y1[index1]*yscale1[index1] + y2[index2]*yscale2[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 + 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] + # 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 + + ############################################################ # add_piece_lin_const_python ############################################################ diff --git a/test/test_distance.py b/test/test_distance.py index 7be0d9b..d98069d 100644 --- a/test/test_distance.py +++ b/test/test_distance.py @@ -130,6 +130,60 @@ def test_spike(): decimal=16) +def test_spike_sync(): + spikes1 = np.array([1.0, 2.0, 3.0]) + spikes2 = np.array([2.1]) + spikes1 = spk.add_auxiliary_spikes(spikes1, 4.0) + spikes2 = spk.add_auxiliary_spikes(spikes2, 4.0) + for k in xrange(1, 3): + assert_almost_equal(spk.spike_sync_distance(spikes1, spikes2, k=k), + 0.5, decimal=16) + + spikes2 = np.array([3.1]) + spikes2 = spk.add_auxiliary_spikes(spikes2, 4.0) + for k in xrange(1, 3): + assert_almost_equal(spk.spike_sync_distance(spikes1, spikes2, k=k), + 0.5, decimal=16) + + spikes2 = np.array([1.1]) + spikes2 = spk.add_auxiliary_spikes(spikes2, 4.0) + for k in xrange(1, 3): + assert_almost_equal(spk.spike_sync_distance(spikes1, spikes2, k=k), + 0.5, decimal=16) + + spikes2 = np.array([0.9]) + spikes2 = spk.add_auxiliary_spikes(spikes2, 4.0) + for k in xrange(1, 3): + assert_almost_equal(spk.spike_sync_distance(spikes1, spikes2, k=k), + 0.5, decimal=16) + + spikes1 = np.array([100, 300, 400, 405, 410, 500, 700, 800, + 805, 810, 815, 900]) + spikes2 = np.array([100, 200, 205, 210, 295, 350, 400, 510, + 600, 605, 700, 910]) + spikes3 = np.array([100, 180, 198, 295, 412, 420, 510, 640, + 695, 795, 820, 920]) + spikes1 = spk.add_auxiliary_spikes(spikes1, 1000) + spikes2 = spk.add_auxiliary_spikes(spikes2, 1000) + spikes3 = spk.add_auxiliary_spikes(spikes3, 1000) + for k in xrange(1, 10): + assert_almost_equal(spk.spike_sync_distance(spikes1, spikes2, k=k), + 0.5, decimal=15) + assert_almost_equal(spk.spike_sync_distance(spikes1, spikes3, k=k), + 0.5, decimal=15) + assert_almost_equal(spk.spike_sync_distance(spikes2, spikes3, k=k), + 0.5, decimal=15) + + f1 = spk.spike_sync_profile(spikes1, spikes2, k=1) + f2 = spk.spike_sync_profile(spikes1, spikes3, k=1) + f3 = spk.spike_sync_profile(spikes2, spikes3, k=1) + f = spk.spike_sync_profile_multi([spikes1, spikes2, spikes3], k=1) + # hands on definition of the average multivariate spike synchronization + expected = (f1.integral() + f2.integral() + f3.integral()) / \ + (len(f1.y)+len(f2.y)+len(f3.y)-3) + assert_almost_equal(f.avrg(), expected, decimal=15) + + def check_multi_profile(profile_func, profile_func_multi): # generate spike trains: t1 = spk.add_auxiliary_spikes(np.array([0.2, 0.4, 0.6, 0.7]), 1.0) -- cgit v1.2.3