final version for spike sync

5 files changed, 102 insertions, 135 deletions

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)