+spike synchronization, python impl only

4 files changed, 125 insertions, 1 deletions

diff --git a/pyspike/__init__.py b/pyspike/__init__.py
index d700e7a..dca5722 100644
--- a/pyspike/__init__.py
+++ b/pyspike/__init__.py
@@ -9,7 +9,9 @@ __all__ = ["function", "distances", "spikes"]
 from function import PieceWiseConstFunc, PieceWiseLinFunc, average_profile
 from distances import isi_profile, isi_distance, \
     spike_profile, spike_distance, \
+    spike_sync_profile, \
     isi_profile_multi, isi_distance_multi, isi_distance_matrix, \
-    spike_profile_multi, spike_distance_multi, spike_distance_matrix
+    spike_profile_multi, spike_distance_multi, spike_distance_matrix, \
+    spike_sync_profile_multi
 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/distances.py b/pyspike/distances.py
index 34f7d78..fbedce5 100644
--- a/pyspike/distances.py
+++ b/pyspike/distances.py
@@ -9,6 +9,7 @@ Distributed under the BSD License
 
 import numpy as np
 import threading
+from functools import partial
 
 from pyspike import PieceWiseConstFunc, PieceWiseLinFunc
 
@@ -129,6 +130,42 @@ def spike_distance(spikes1, spikes2, interval=None):
 
 
 ############################################################
+# spike_sync_profile
+############################################################
+def spike_sync_profile(spikes1, spikes2, k=3):
+
+    assert k > 0
+
+    # cython implementation
+    try:
+        from cython_distance import cumulative_sync_cython \
+            as cumulative_sync_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.")
+        # use python backend
+        from python_backend import cumulative_sync_python \
+            as cumulative_sync_impl
+
+    st, c = cumulative_sync_impl(spikes1, spikes2)
+
+    # print c
+    # print 2*(c[-1]-c[0])/(len(spikes1)+len(spikes2)-2)
+
+    dc = np.zeros(len(c))
+    dc[k:-k] = (c[2*k:] - c[:-2*k]) / k
+    for n in xrange(1, k):
+        dc[n] = (c[2*n] - c[0]) / k
+        dc[-n-1] = (c[-1]-c[-2*n-1]) / k
+    dc[0] = dc[1]
+    dc[-1] = dc[-2]
+    # dc[-1] = (c[-1]-c[-2])/k
+    # print dc
+    return PieceWiseConstFunc(st, dc)
+
+
+############################################################
 # _generic_profile_multi
 ############################################################
 def _generic_profile_multi(spike_trains, pair_distance_func, indices=None):
@@ -279,6 +316,28 @@ def spike_profile_multi(spike_trains, indices=None):
 
 
 ############################################################
+# spike_profile_multi
+############################################################
+def spike_sync_profile_multi(spike_trains, indices=None, k=3):
+    """ 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>
+
+    :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`
+
+    """
+    prof_func = partial(spike_sync_profile, k=k)
+    return _generic_profile_multi(spike_trains, prof_func, indices)
+
+
+############################################################
 # spike_distance_multi
 ############################################################
 def spike_distance_multi(spike_trains, indices=None, interval=None):
diff --git a/pyspike/function.py b/pyspike/function.py
index 0c0a391..662606c 100644
--- a/pyspike/function.py
+++ b/pyspike/function.py
@@ -136,6 +136,15 @@ 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.
diff --git a/pyspike/python_backend.py b/pyspike/python_backend.py
index 874c689..b85262d 100644
--- a/pyspike/python_backend.py
+++ b/pyspike/python_backend.py
@@ -189,6 +189,60 @@ def spike_distance_python(spikes1, spikes2):
 
 
 ############################################################
+# cumulative_sync_python
+############################################################
+def cumulative_sync_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]])
+    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:
+        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] = c[n-1]
+            else:
+                c[n] = 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] = c[n-1]
+            else:
+                c[n] = c[n-1]+1
+        else:   # spikes1[i+1] = spikes2[j+1]
+            j += 1
+            i += 1
+            if i == N1-1 or j == N2-1:
+                break
+            n += 1
+            st[n] = spikes1[i]
+            c[n] = c[n-1]
+            n += 1
+            st[n] = spikes1[i]
+            c[n] = c[n-1]+1
+    c[0] = 0
+    st[0] = spikes1[0]
+    st[-1] = spikes1[-1]
+
+    return st, c
+
+
+############################################################
 # add_piece_wise_const_python
 ############################################################
 def add_piece_wise_const_python(x1, y1, x2, y2):