renamed IntervalSequence -> DiscreteFunction, cython implementation of add

5 files changed, 255 insertions, 194 deletions

diff --git a/pyspike/__init__.py b/pyspike/__init__.py
index 74d52c5..1b18569 100644
--- a/pyspike/__init__.py
+++ b/pyspike/__init__.py
@@ -7,7 +7,7 @@ Distributed under the BSD License
 __all__ = ["function", "distances", "spikes"]
 
 from function import PieceWiseConstFunc, PieceWiseLinFunc, \
-    MultipleValueSequence, average_profile
+    DiscreteFunction, average_profile
 from distances import isi_profile, isi_distance, \
     spike_profile, spike_distance, \
     spike_sync_profile, spike_sync_distance, \
diff --git a/pyspike/cython_add.pyx b/pyspike/cython_add.pyx
index bfbe208..817799e 100644
--- a/pyspike/cython_add.pyx
+++ b/pyspike/cython_add.pyx
@@ -172,3 +172,62 @@ def add_piece_wise_lin_cython(double[:] x1, double[:] y11, double[:] y12,
     return (np.array(x_new[:index+2]),
             np.array(y1_new[:index+1]), 
             np.array(y2_new[:index+1]))
+
+
+############################################################
+# add_discrete_function_cython
+############################################################
+def add_discrete_function_cython(double[:] x1, double[:] y1, double[:] mp1,
+                                 double[:] x2, double[:] y2, double[:] mp2):
+
+    cdef double[:] x_new = np.empty(len(x1) + len(x2))
+    cdef double[:] y_new = np.empty_like(x_new)
+    cdef double[:] mp_new = np.empty_like(x_new)
+    cdef int index1 = 0
+    cdef int index2 = 0
+    cdef int index = 0
+    cdef int N1 = len(y1)
+    cdef int N2 = len(y2)
+    x_new[0] = x1[0]
+    while (index1+1 < N1) and (index2+1 < N2):
+        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]
+            index1 += 1
+            index2 += 1
+            index += 1
+            x_new[index] = x1[index1]
+            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(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(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
+    # only use the data that was actually filled
+    return x_new[:index+1], y_new[:index+1], mp_new[:index+1]
diff --git a/pyspike/distances.py b/pyspike/distances.py
index 8bde724..0f0efa9 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, MultipleValueSequence
+from pyspike import PieceWiseConstFunc, PieceWiseLinFunc, DiscreteFunction
 
 
 ############################################################
@@ -148,7 +148,7 @@ Falling back to slow python backend.")
 
     times, coincidences, multiplicity = coincidence_impl(spikes1, spikes2)
 
-    return MultipleValueSequence(times, coincidences, multiplicity)
+    return DiscreteFunction(times, coincidences, multiplicity)
 
 
 ############################################################
@@ -328,7 +328,7 @@ def spike_sync_profile_multi(spike_trains, indices=None):
                     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`
+    :rtype: :class:`pyspike.function.DiscreteFunction`
 
     """
     prof_func = partial(spike_sync_profile)
diff --git a/pyspike/function.py b/pyspike/function.py
index f10c136..6fb7537 100644
--- a/pyspike/function.py
+++ b/pyspike/function.py
@@ -171,140 +171,6 @@ that PySpike is installed by running\n 'python setup.py build_ext --inplace'! \
 
 
 ##############################################################
-# MultipleValueSequence
-##############################################################
-class MultipleValueSequence(object):
-    """ A class representing a sequence of values defined in some interval.
-    """
-
-    def __init__(self, x, y, multiplicity):
-        """ Constructs the value sequence.
-
-        :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.mp = np.array(multiplicity)
-
-    def copy(self):
-        """ Returns a copy of itself
-
-        :rtype: :class:`IntervalSequence`
-        """
-        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`
-        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) and \
-            np.allclose(self.mp, other.mp, atol=eps, rtol=0.0)
-
-    def get_plottable_data(self, k=0):
-        """ 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")
-        """
-
-        if k > 0:
-            raise NotImplementedError()
-
-        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
-        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 = 1.0*np.sum(self.y[1:-1])
-        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() / np.sum(self.mp[1:-1])
-        else:
-            raise NotImplementedError()
-
-    def add(self, f):
-        """ Adds another `MultipleValueSequence` function to this function.
-        Note: only functions defined on the same interval can be summed.
-
-        :param f: :class:`MultipleValueSequence` 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_multiple_value_sequence_python as \
-                add_multiple_value_sequence_impl
-
-        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
-
-        :param fac: Value to multiply
-        :type fac: double
-        :rtype: None
-        """
-        self.y *= fac
-
-
-##############################################################
 # PieceWiseLinFunc
 ##############################################################
 class PieceWiseLinFunc:
@@ -489,6 +355,141 @@ that PySpike is installed by running\n 'python setup.py build_ext --inplace'! \
         self.y2 *= fac
 
 
+##############################################################
+# DiscreteFunction
+##############################################################
+class DiscreteFunction(object):
+    """ A class representing values defined on a discrete set of points.
+    """
+
+    def __init__(self, x, y, multiplicity):
+        """ Constructs the discrete function.
+
+        :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.mp = np.array(multiplicity)
+
+    def copy(self):
+        """ Returns a copy of itself
+
+        :rtype: :class:`DiscreteFunction`
+        """
+        return DiscreteFunction(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`
+        precision.
+
+        :param other: another :class:`DiscreteFunction`
+        :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) and \
+            np.allclose(self.mp, other.mp, atol=eps, rtol=0.0)
+
+    def get_plottable_data(self, k=0):
+        """ 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="Discrete function")
+        """
+
+        if k > 0:
+            raise NotImplementedError()
+
+        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
+        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 = 1.0*np.sum(self.y[1:-1])
+        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() / np.sum(self.mp[1:-1])
+        else:
+            raise NotImplementedError()
+
+    def add(self, f):
+        """ Adds another `DiscreteFunction` function to this function.
+        Note: only functions defined on the same interval can be summed.
+
+        :param f: :class:`DiscreteFunction` 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_discrete_function_cython as \
+                add_discrete_function_impl
+        except ImportError:
+            print("Warning: add_discrete_function_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_discrete_function_python as \
+                add_discrete_function_impl
+
+        self.x, self.y, self.mp = \
+            add_discrete_function_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
+
+        :param fac: Value to multiply
+        :type fac: double
+        :rtype: None
+        """
+        self.y *= fac
+
+
 def average_profile(profiles):
     """ Computes the average profile from the given ISI- or SPIKE-profiles.
 
diff --git a/pyspike/python_backend.py b/pyspike/python_backend.py
index bbbd572..481daf9 100644
--- a/pyspike/python_backend.py
+++ b/pyspike/python_backend.py
@@ -358,62 +358,6 @@ def add_piece_wise_const_python(x1, y1, x2, y2):
 
 
 ############################################################
-# add_multiple_value_sequence_python
-############################################################
-def add_multiple_value_sequence_python(x1, y1, mp1, x2, y2, mp2):
-
-    x_new = np.empty(len(x1) + len(x2))
-    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
-    while (index1+1 < len(y1)) and (index2+1 < len(y2)):
-        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]
-            index1 += 1
-            index2 += 1
-            index += 1
-            x_new[index] = x1[index1]
-            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(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(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
-    # only use the data that was actually filled
-    return x_new[:index+1], y_new[:index+1], mp_new[:index+1]
-
-
-############################################################
 # add_piece_lin_const_python
 ############################################################
 def add_piece_wise_lin_python(x1, y11, y12, x2, y21, y22):
@@ -482,3 +426,60 @@ def add_piece_wise_lin_python(x1, y11, y12, x2, y21, y22):
     y2_new[index] = y12[-1]+y22[-1]
     # only use the data that was actually filled
     return x_new[:index+2], y1_new[:index+1], y2_new[:index+1]
+
+
+############################################################
+# add_discrete_function_python
+############################################################
+def add_discrete_function_python(x1, y1, mp1, x2, y2, mp2):
+
+    x_new = np.empty(len(x1) + len(x2))
+    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
+    while (index1+1 < len(y1)) and (index2+1 < len(y2)):
+        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]
+            index1 += 1
+            index2 += 1
+            index += 1
+            x_new[index] = x1[index1]
+            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(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(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
+    # only use the data that was actually filled
+    return x_new[:index+1], y_new[:index+1], mp_new[:index+1]
+