spike sync

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:`<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,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.
@@ -180,6 +171,146 @@ that PySpike is installed by running\n 'python setup.py build_ext --inplace'! \
 
 
 ##############################################################
+# 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
 ##############################################################
 class 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]
 
@@ -341,6 +341,86 @@ def add_piece_wise_const_python(x1, y1, x2, y2):
 
 
 ############################################################
+# 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
 ############################################################
 def add_piece_wise_lin_python(x1, y11, y12, x2, y21, y22):
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)