Merge pull request #2 from mariomulansky/restructuring

Restructuring

16 files changed, 1130 insertions, 1113 deletions

diff --git a/pyspike/DiscreteFunc.py b/pyspike/DiscreteFunc.py
new file mode 100644
index 0000000..2283e03
--- /dev/null
+++ b/pyspike/DiscreteFunc.py
@@ -0,0 +1,242 @@
+"""
+Class representing discrete functions.
+
+Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net>
+
+Distributed under the BSD License
+
+"""
+from __future__ import print_function
+
+import numpy as np
+import collections
+
+
+##############################################################
+# DiscreteFunc
+##############################################################
+class DiscreteFunc(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:`DiscreteFunc`
+        """
+        return DiscreteFunc(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:`DiscreteFunc`
+        :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, averaging_window_size=0):
+        """ Returns two arrays containing x- and y-coordinates for plotting
+        the interval sequence. The optional parameter `averaging_window_size`
+        determines the size of an averaging window to smoothen the profile. If
+        this value is 0, no averaging is performed.
+
+        :param averaging_window_size: size of the averaging window, default=0.
+        :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 averaging_window_size > 0:
+            # for the averaged profile we have to take the multiplicity into
+            # account. values with higher multiplicity should be consider as if
+            # they appeared several times. Hence we can not know how many
+            # entries we have to consider to the left and right. Rather, we
+            # will iterate until some wanted multiplicity is reached.
+
+            # the first value in self.mp contains the number of averaged
+            # profiles without any possible extra multiplicities
+            # (by implementation)
+            expected_mp = (averaging_window_size+1) * int(self.mp[0])
+            y_plot = np.zeros_like(self.y)
+            # compute the values in a loop, could be done in cython if required
+            for i in xrange(len(y_plot)):
+
+                if self.mp[i] >= expected_mp:
+                    # the current value contains already all the wanted
+                    # multiplicity
+                    y_plot[i] = self.y[i]/self.mp[i]
+                    continue
+
+                # first look to the right
+                y = self.y[i]
+                mp_r = self.mp[i]
+                j = i+1
+                while j < len(y_plot):
+                    if mp_r+self.mp[j] < expected_mp:
+                        # if we still dont reach the required multiplicity
+                        # we take the whole value
+                        y += self.y[j]
+                        mp_r += self.mp[j]
+                    else:
+                        # otherwise, just some fraction
+                        y += self.y[j] * (expected_mp - mp_r)/self.mp[j]
+                        mp_r += (expected_mp - mp_r)
+                        break
+                    j += 1
+
+                # same story to the left
+                mp_l = self.mp[i]
+                j = i-1
+                while j >= 0:
+                    if mp_l+self.mp[j] < expected_mp:
+                        y += self.y[j]
+                        mp_l += self.mp[j]
+                    else:
+                        y += self.y[j] * (expected_mp - mp_l)/self.mp[j]
+                        mp_l += (expected_mp - mp_l)
+                        break
+                    j -= 1
+                y_plot[i] = y/(mp_l+mp_r-self.mp[i])
+            return 1.0*self.x, y_plot
+
+        else:  # k = 0
+
+            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 discrete
+        function, this amounts to the sum over all values divided by the total
+        multiplicity.
+
+        :param interval: integration interval given as a pair of floats, or a
+                         sequence of pairs in case of multiple intervals, if
+                         None the integral over the whole function is computed.
+        :type interval: Pair, sequence of pairs, or None.
+        :returns: the integral
+        :rtype: float
+        """
+
+        def get_indices(ival):
+            """ Retuns the indeces surrounding the given interval"""
+            start_ind = np.searchsorted(self.x, ival[0], side='right')
+            end_ind = np.searchsorted(self.x, ival[1], side='left')
+            assert start_ind > 0 and end_ind < len(self.x), \
+                "Invalid averaging interval"
+            return start_ind, end_ind
+
+        if interval is None:
+            # no interval given, integrate over the whole spike train
+            # don't count the first value, which is zero by definition
+            return 1.0 * np.sum(self.y[1:-1]) / np.sum(self.mp[1:-1])
+
+        # check if interval is as sequence
+        assert isinstance(interval, collections.Sequence), \
+            "Invalid value for `interval`. None, Sequence or Tuple expected."
+        # check if interval is a sequence of intervals
+        if not isinstance(interval[0], collections.Sequence):
+            # find the indices corresponding to the interval
+            start_ind, end_ind = get_indices(interval)
+            return (np.sum(self.y[start_ind:end_ind]) /
+                    np.sum(self.mp[start_ind:end_ind]))
+        else:
+            value = 0.0
+            multiplicity = 0.0
+            for ival in interval:
+                # find the indices corresponding to the interval
+                start_ind, end_ind = get_indices(ival)
+                value += np.sum(self.y[start_ind:end_ind])
+                multiplicity += np.sum(self.mp[start_ind:end_ind])
+        return value/multiplicity
+
+    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
+        """
+        return self.integral(interval)
+
+    def add(self, f):
+        """ Adds another `DiscreteFunc` function to this function.
+        Note: only functions defined on the same interval can be summed.
+
+        :param f: :class:`DiscreteFunc` 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.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 cython.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.
+
+    :param profiles: list of :class:`PieceWiseConstFunc` or
+                     :class:`PieceWiseLinFunc` representing ISI- or
+                     SPIKE-profiles to be averaged.
+    :returns: the averages profile :math:`<S_{isi}>` or :math:`<S_{spike}>`.
+    :rtype: :class:`PieceWiseConstFunc` or :class:`PieceWiseLinFunc`
+    """
+    assert len(profiles) > 1
+
+    avrg_profile = profiles[0].copy()
+    for i in xrange(1, len(profiles)):
+        avrg_profile.add(profiles[i])
+    avrg_profile.mul_scalar(1.0/len(profiles))  # normalize
+
+    return avrg_profile
diff --git a/pyspike/PieceWiseConstFunc.py b/pyspike/PieceWiseConstFunc.py
new file mode 100644
index 0000000..dc57ab1
--- /dev/null
+++ b/pyspike/PieceWiseConstFunc.py
@@ -0,0 +1,168 @@
+"""
+Class representing piece-wise constant functions.
+
+Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net>
+
+Distributed under the BSD License
+
+"""
+from __future__ import print_function
+
+import numpy as np
+import collections
+
+
+##############################################################
+# PieceWiseConstFunc
+##############################################################
+class PieceWiseConstFunc(object):
+    """ A class representing a piece-wise constant function. """
+
+    def __init__(self, x, y):
+        """ Constructs the piece-wise const function.
+
+        :param x: array of length N+1 defining the edges of the intervals of
+                  the pwc function.
+        :param y: array of length N defining the function values at the
+                  intervals.
+        """
+        # convert parameters to arrays, also ensures copying
+        self.x = np.array(x)
+        self.y = np.array(y)
+
+    def copy(self):
+        """ Returns a copy of itself
+
+        :rtype: :class:`PieceWiseConstFunc`
+        """
+        return PieceWiseConstFunc(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:`PieceWiseConstFunc`
+        :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 piece-wise function.
+
+        :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
+        y_plot[1::2] = self.y
+
+        return x_plot, y_plot
+
+    def integral(self, interval=None):
+        """ Returns the integral over the given interval.
+
+        :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
+            a = np.sum((self.x[1:]-self.x[:-1]) * self.y)
+        else:
+            # find the indices corresponding to the interval
+            start_ind = np.searchsorted(self.x, interval[0], side='right')
+            end_ind = np.searchsorted(self.x, interval[1], side='left')-1
+            assert start_ind > 0 and end_ind < len(self.x), \
+                "Invalid averaging interval"
+            # first the contribution from between the indices
+            a = np.sum((self.x[start_ind+1:end_ind+1] -
+                        self.x[start_ind:end_ind]) *
+                       self.y[start_ind:end_ind])
+            # correction from start to first index
+            a += (self.x[start_ind]-interval[0]) * self.y[start_ind-1]
+            # correction from last index to end
+            a += (interval[1]-self.x[end_ind]) * self.y[end_ind]
+        return a
+
+    def avrg(self, interval=None):
+        """ Computes the average of the piece-wise const function:
+        :math:`a = 1/T int_0^T f(x) dx` where T is the length of the interval.
+
+        :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
+            return self.integral() / (self.x[-1]-self.x[0])
+
+        # check if interval is as sequence
+        assert isinstance(interval, collections.Sequence), \
+            "Invalid value for `interval`. None, Sequence or Tuple expected."
+        # check if interval is a sequence of intervals
+        if not isinstance(interval[0], collections.Sequence):
+            # just one interval
+            a = self.integral(interval) / (interval[1]-interval[0])
+        else:
+            # several intervals
+            a = 0.0
+            int_length = 0.0
+            for ival in interval:
+                a += self.integral(ival)
+                int_length += ival[1] - ival[0]
+            a /= int_length
+        return a
+
+    def add(self, f):
+        """ Adds another PieceWiseConst function to this function.
+        Note: only functions defined on the same interval can be summed.
+
+        :param f: :class:`PieceWiseConstFunc` 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.cython_add import add_piece_wise_const_cython as \
+                add_piece_wise_const_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 cython.python_backend import add_piece_wise_const_python as \
+                add_piece_wise_const_impl
+
+        self.x, self.y = add_piece_wise_const_impl(self.x, self.y, f.x, f.y)
+
+    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
diff --git a/pyspike/PieceWiseLinFunc.py b/pyspike/PieceWiseLinFunc.py
new file mode 100644
index 0000000..bc0aa2a
--- /dev/null
+++ b/pyspike/PieceWiseLinFunc.py
@@ -0,0 +1,197 @@
+"""
+Class representing piece-wise linear functions.
+
+Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net>
+
+Distributed under the BSD License
+
+"""
+from __future__ import print_function
+
+import numpy as np
+import collections
+
+
+##############################################################
+# PieceWiseLinFunc
+##############################################################
+class PieceWiseLinFunc:
+    """ A class representing a piece-wise linear function. """
+
+    def __init__(self, x, y1, y2):
+        """ Constructs the piece-wise linear function.
+
+        :param x: array of length N+1 defining the edges of the intervals of
+                  the pwc function.
+        :param y1: array of length N defining the function values at the left
+                  of the intervals.
+        :param y2: array of length N defining the function values at the right
+                  of the intervals.
+        """
+        # convert to array, which also ensures copying
+        self.x = np.array(x)
+        self.y1 = np.array(y1)
+        self.y2 = np.array(y2)
+
+    def copy(self):
+        """ Returns a copy of itself
+
+        :rtype: :class:`PieceWiseLinFunc`
+        """
+        return PieceWiseLinFunc(self.x, self.y1, self.y2)
+
+    def almost_equal(self, other, decimal=14):
+        """ Checks if the function is equal to another function up to `decimal`
+        precision.
+
+        :param other: another :class:`PieceWiseLinFunc`
+        :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.y1, other.y1, atol=eps, rtol=0.0) and \
+            np.allclose(self.y2, other.y2, atol=eps, rtol=0.0)
+
+    def get_plottable_data(self):
+        """ Returns two arrays containing x- and y-coordinates for immeditate
+        plotting of the piece-wise function.
+
+        :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_like(x_plot)
+        y_plot[0::2] = self.y1
+        y_plot[1::2] = self.y2
+        return x_plot, y_plot
+
+    def integral(self, interval=None):
+        """ Returns the integral over the given interval.
+
+        :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
+        """
+
+        def intermediate_value(x0, x1, y0, y1, x):
+            """ computes the intermediate value of a linear function """
+            return y0 + (y1-y0)*(x-x0)/(x1-x0)
+
+        if interval is None:
+            # no interval given, integrate over the whole spike train
+            integral = np.sum((self.x[1:]-self.x[:-1]) * 0.5*(self.y1+self.y2))
+        else:
+            # find the indices corresponding to the interval
+            start_ind = np.searchsorted(self.x, interval[0], side='right')
+            end_ind = np.searchsorted(self.x, interval[1], side='left')-1
+            assert start_ind > 0 and end_ind < len(self.x), \
+                "Invalid averaging interval"
+            # first the contribution from between the indices
+            integral = np.sum((self.x[start_ind+1:end_ind+1] -
+                               self.x[start_ind:end_ind]) *
+                              0.5*(self.y1[start_ind:end_ind] +
+                                   self.y2[start_ind:end_ind]))
+            # correction from start to first index
+            integral += (self.x[start_ind]-interval[0]) * 0.5 * \
+                        (self.y2[start_ind-1] +
+                         intermediate_value(self.x[start_ind-1],
+                                            self.x[start_ind],
+                                            self.y1[start_ind-1],
+                                            self.y2[start_ind-1],
+                                            interval[0]
+                                            ))
+            # correction from last index to end
+            integral += (interval[1]-self.x[end_ind]) * 0.5 * \
+                        (self.y1[end_ind] +
+                         intermediate_value(self.x[end_ind], self.x[end_ind+1],
+                                            self.y1[end_ind], self.y2[end_ind],
+                                            interval[1]
+                                            ))
+        return integral
+
+    def avrg(self, interval=None):
+        """ Computes the average of the piece-wise linear function:
+        :math:`a = 1/T int_0^T f(x) dx` where T is the length of the interval.
+
+        :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
+            return self.integral() / (self.x[-1]-self.x[0])
+
+        # check if interval is as sequence
+        assert isinstance(interval, collections.Sequence), \
+            "Invalid value for `interval`. None, Sequence or Tuple expected."
+        # check if interval is a sequence of intervals
+        if not isinstance(interval[0], collections.Sequence):
+            # just one interval
+            a = self.integral(interval) / (interval[1]-interval[0])
+        else:
+            # several intervals
+            a = 0.0
+            int_length = 0.0
+            for ival in interval:
+                a += self.integral(ival)
+                int_length += ival[1] - ival[0]
+            a /= int_length
+        return a
+
+    def add(self, f):
+        """ Adds another PieceWiseLin function to this function.
+        Note: only functions defined on the same interval can be summed.
+
+        :param f: :class:`PieceWiseLinFunc` 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"
+
+        # python implementation
+        # from python_backend import add_piece_wise_lin_python
+        # self.x, self.y1, self.y2 = add_piece_wise_lin_python(
+        #     self.x, self.y1, self.y2, f.x, f.y1, f.y2)
+
+        # cython version
+        try:
+            from cython.cython_add import add_piece_wise_lin_cython as \
+                add_piece_wise_lin_impl
+        except ImportError:
+            print("Warning: add_piece_wise_lin_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 cython.python_backend import add_piece_wise_lin_python as \
+                add_piece_wise_lin_impl
+
+        self.x, self.y1, self.y2 = add_piece_wise_lin_impl(
+            self.x, self.y1, self.y2, f.x, f.y1, f.y2)
+
+    def mul_scalar(self, fac):
+        """ Multiplies the function with a scalar value
+
+        :param fac: Value to multiply
+        :type fac: double
+        :rtype: None
+        """
+        self.y1 *= fac
+        self.y2 *= fac
diff --git a/pyspike/__init__.py b/pyspike/__init__.py
index 55687e6..945dd4e 100644
--- a/pyspike/__init__.py
+++ b/pyspike/__init__.py
@@ -1,18 +1,23 @@
 """
-Copyright 2014, Mario Mulansky <mario.mulansky@gmx.net>
+Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net>
 
 Distributed under the BSD License
 """
 
-__all__ = ["function", "distances", "spikes"]
+__all__ = ["isi_distance", "spike_distance", "spike_sync",
+           "spikes", "PieceWiseConstFunc", "PieceWiseLinFunc",
+           "DiscreteFunc"]
 
-from function import PieceWiseConstFunc, PieceWiseLinFunc, \
-    DiscreteFunction, average_profile
-from distances import isi_profile, isi_distance, \
-    spike_profile, spike_distance, \
-    spike_sync_profile, spike_sync, \
-    isi_profile_multi, isi_distance_multi, isi_distance_matrix, \
-    spike_profile_multi, spike_distance_multi, spike_distance_matrix, \
+from PieceWiseConstFunc import PieceWiseConstFunc
+from PieceWiseLinFunc import PieceWiseLinFunc
+from DiscreteFunc import DiscreteFunc
+
+from isi_distance import isi_profile, isi_distance, isi_profile_multi,\
+    isi_distance_multi, isi_distance_matrix
+from spike_distance import spike_profile, spike_distance, spike_profile_multi,\
+    spike_distance_multi, spike_distance_matrix
+from spike_sync import spike_sync_profile, spike_sync,\
     spike_sync_profile_multi, spike_sync_multi, spike_sync_matrix
+
 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/cython/__init__.py b/pyspike/cython/__init__.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/pyspike/cython/__init__.py
diff --git a/pyspike/cython_add.pyx b/pyspike/cython/cython_add.pyx
index ac64005..ac64005 100644
--- a/pyspike/cython_add.pyx
+++ b/pyspike/cython/cython_add.pyx
diff --git a/pyspike/cython_distance.pyx b/pyspike/cython/cython_distance.pyx
index 489aab9..489aab9 100644
--- a/pyspike/cython_distance.pyx
+++ b/pyspike/cython/cython_distance.pyx
diff --git a/pyspike/python_backend.py b/pyspike/cython/python_backend.py
index 481daf9..481daf9 100644
--- a/pyspike/python_backend.py
+++ b/pyspike/cython/python_backend.py
diff --git a/pyspike/distances.py b/pyspike/distances.py
deleted file mode 100644
index 5476b6f..0000000
--- a/pyspike/distances.py
+++ /dev/null
@@ -1,504 +0,0 @@
-""" distances.py
-
-Module containing several functions to compute spike distances
-
-Copyright 2014, Mario Mulansky <mario.mulansky@gmx.net>
-
-Distributed under the BSD License
-"""
-
-import numpy as np
-import threading
-from functools import partial
-
-from pyspike import PieceWiseConstFunc, PieceWiseLinFunc, DiscreteFunction
-
-
-############################################################
-# isi_profile
-############################################################
-def isi_profile(spikes1, spikes2):
-    """ Computes the isi-distance profile :math:`S_{isi}(t)` of the two given
-    spike trains. Retruns the profile as a PieceWiseConstFunc object. The S_isi
-    values are defined positive S_isi(t)>=0.  The spike trains are expected
-    to have auxiliary spikes at the beginning and end of the interval. Use the
-    function add_auxiliary_spikes to add those spikes to the spike train.
-
-    :param spikes1: ordered array of spike times with auxiliary spikes.
-    :param spikes2: ordered array of spike times with auxiliary spikes.
-    :returns: The isi-distance profile :math:`S_{isi}(t)`
-    :rtype: :class:`pyspike.function.PieceWiseConstFunc`
-
-    """
-    # check for auxiliary spikes - first and last spikes should be identical
-    assert spikes1[0] == spikes2[0], \
-        "Given spike trains seems not to have auxiliary spikes!"
-    assert spikes1[-1] == spikes2[-1], \
-        "Given spike trains seems not to have auxiliary spikes!"
-
-    # load cython implementation
-    try:
-        from cython_distance import isi_distance_cython as isi_distance_impl
-    except ImportError:
-        print("Warning: isi_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 isi_distance_python as isi_distance_impl
-
-    times, values = isi_distance_impl(spikes1, spikes2)
-    return PieceWiseConstFunc(times, values)
-
-
-############################################################
-# isi_distance
-############################################################
-def isi_distance(spikes1, spikes2, interval=None):
-    """ Computes the isi-distance I of the given spike trains. The
-    isi-distance is the integral over the isi distance profile
-    :math:`S_{isi}(t)`:
-
-    .. math:: I = \int_{T_0}^{T_1} S_{isi}(t) dt.
-
-    :param spikes1: ordered array of spike times with auxiliary spikes.
-    :param spikes2: ordered array of spike times with auxiliary spikes.
-    :param interval: averaging interval given as a pair of floats (T0, T1),
-                     if None the average over the whole function is computed.
-    :type interval: Pair of floats or None.
-    :returns: The isi-distance I.
-    :rtype: double
-    """
-    return isi_profile(spikes1, spikes2).avrg(interval)
-
-
-############################################################
-# spike_profile
-############################################################
-def spike_profile(spikes1, spikes2):
-    """ Computes the spike-distance profile S_spike(t) of the two given spike
-    trains. Returns the profile as a PieceWiseLinFunc object. The S_spike
-    values are defined positive S_spike(t)>=0. The spike trains are expected to
-    have auxiliary spikes at the beginning and end of the interval. Use the
-    function add_auxiliary_spikes to add those spikes to the spike train.
-
-    :param spikes1: ordered array of spike times with auxiliary spikes.
-    :param spikes2: ordered array of spike times with auxiliary spikes.
-    :returns: The spike-distance profile :math:`S_{spike}(t)`.
-    :rtype: :class:`pyspike.function.PieceWiseLinFunc`
-
-    """
-    # check for auxiliary spikes - first and last spikes should be identical
-    assert spikes1[0] == spikes2[0], \
-        "Given spike trains seems not to have auxiliary spikes!"
-    assert spikes1[-1] == spikes2[-1], \
-        "Given spike trains seems not to have auxiliary spikes!"
-
-    # cython implementation
-    try:
-        from cython_distance import spike_distance_cython \
-            as spike_distance_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 spike_distance_python as spike_distance_impl
-
-    times, y_starts, y_ends = spike_distance_impl(spikes1, spikes2)
-    return PieceWiseLinFunc(times, y_starts, y_ends)
-
-
-############################################################
-# spike_distance
-############################################################
-def spike_distance(spikes1, spikes2, interval=None):
-    """ Computes the spike-distance S of the given spike trains. The
-    spike-distance is the integral over the isi distance profile S_spike(t):
-
-    .. math:: S = \int_{T_0}^{T_1} S_{spike}(t) dt.
-
-    :param spikes1: ordered array of spike times with auxiliary spikes.
-    :param spikes2: ordered array of spike times with auxiliary spikes.
-    :param interval: averaging interval given as a pair of floats (T0, T1),
-                     if None the average over the whole function is computed.
-    :type interval: Pair of floats or None.
-    :returns: The spike-distance.
-    :rtype: double
-
-    """
-    return spike_profile(spikes1, spikes2).avrg(interval)
-
-
-############################################################
-# spike_sync_profile
-############################################################
-def spike_sync_profile(spikes1, spikes2):
-    """ Computes the spike-synchronization profile S_sync(t) of the two given
-    spike trains. Returns the profile as a DiscreteFunction object. The S_sync
-    values are either 1 or 0, indicating the presence or absence of a
-    coincidence. The spike trains are expected to have auxiliary spikes at the
-    beginning and end of the interval. Use the function add_auxiliary_spikes to
-    add those spikes to the spike train.
-
-    :param spikes1: ordered array of spike times with auxiliary spikes.
-    :param spikes2: ordered array of spike times with auxiliary spikes.
-    :returns: The spike-distance profile :math:`S_{sync}(t)`.
-    :rtype: :class:`pyspike.function.DiscreteFunction`
-
-    """
-
-    # cython implementation
-    try:
-        from cython_distance import coincidence_cython \
-            as coincidence_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 coincidence_python \
-            as coincidence_impl
-
-    times, coincidences, multiplicity = coincidence_impl(spikes1, spikes2)
-
-    return DiscreteFunction(times, coincidences, multiplicity)
-
-
-############################################################
-# spike_sync
-############################################################
-def spike_sync(spikes1, spikes2, interval=None):
-    """ Computes the spike synchronization value SYNC of the given spike
-    trains. The spike synchronization value is the computed as the total number
-    of coincidences divided by the total number of spikes:
-
-    .. math:: SYNC = \sum_n C_n / N.
-
-    :param spikes1: ordered array of spike times with auxiliary spikes.
-    :param spikes2: ordered array of spike times with auxiliary spikes.
-    :param interval: averaging interval given as a pair of floats (T0, T1),
-                     if None the average over the whole function is computed.
-    :type interval: Pair of floats or None.
-    :returns: The spike synchronization value.
-    :rtype: double
-    """
-    return spike_sync_profile(spikes1, spikes2).avrg(interval)
-
-
-############################################################
-# _generic_profile_multi
-############################################################
-def _generic_profile_multi(spike_trains, pair_distance_func, indices=None):
-    """ Internal implementation detail, don't call this function directly,
-    use isi_profile_multi or spike_profile_multi instead.
-
-    Computes the multi-variate distance for a set of spike-trains using the
-    pair_dist_func to compute pair-wise distances. That is it computes the
-    average distance of all pairs of spike-trains:
-    :math:`S(t) = 2/((N(N-1)) sum_{<i,j>} S_{i,j}`,
-    where the sum goes over all pairs <i,j>.
-    Args:
-    - spike_trains: list of spike trains
-    - pair_distance_func: function computing the distance of two spike trains
-    - indices: list of indices defining which spike trains to use,
-    if None all given spike trains are used (default=None)
-    Returns:
-    - The averaged multi-variate distance of all pairs
-    """
-    if indices is None:
-        indices = np.arange(len(spike_trains))
-    indices = np.array(indices)
-    # check validity of indices
-    assert (indices < len(spike_trains)).all() and (indices >= 0).all(), \
-        "Invalid index list."
-    # generate a list of possible index pairs
-    pairs = [(indices[i], j) for i in range(len(indices)) for j in indices[i+1:]]
-    # start with first pair
-    (i, j) = pairs[0]
-    average_dist = pair_distance_func(spike_trains[i], spike_trains[j])
-    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
-    return average_dist, len(pairs)
-
-
-############################################################
-# multi_distance_par
-############################################################
-def _multi_distance_par(spike_trains, pair_distance_func, indices=None):
-    """ parallel implementation of the multi-distance. Not currently used as
-    it does not improve the performance.
-    """
-
-    num_threads = 2
-    lock = threading.Lock()
-
-    def run(spike_trains, index_pairs, average_dist):
-        (i, j) = index_pairs[0]
-        # print(i,j)
-        this_avrg = pair_distance_func(spike_trains[i], spike_trains[j])
-        for (i, j) in index_pairs[1:]:
-            # print(i,j)
-            current_dist = pair_distance_func(spike_trains[i], spike_trains[j])
-            this_avrg.add(current_dist)
-        with lock:
-            average_dist.add(this_avrg)
-
-    if indices is None:
-        indices = np.arange(len(spike_trains))
-    indices = np.array(indices)
-    # check validity of indices
-    assert (indices < len(spike_trains)).all() and (indices >= 0).all(), \
-        "Invalid index list."
-    # generate a list of possible index pairs
-    pairs = [(indices[i], j) for i in range(len(indices)) for j in indices[i+1:]]
-    num_pairs = len(pairs)
-
-    # start with first pair
-    (i, j) = pairs[0]
-    average_dist = pair_distance_func(spike_trains[i], spike_trains[j])
-    # remove the one we already computed
-    pairs = pairs[1:]
-    # distribute the rest into num_threads pieces
-    clustered_pairs = [pairs[n::num_threads] for n in xrange(num_threads)]
-
-    threads = []
-    for pairs in clustered_pairs:
-        t = threading.Thread(target=run, args=(spike_trains, pairs,
-                                               average_dist))
-        threads.append(t)
-        t.start()
-    for t in threads:
-        t.join()
-    average_dist.mul_scalar(1.0/num_pairs)  # normalize
-    return average_dist
-
-
-############################################################
-# isi_profile_multi
-############################################################
-def isi_profile_multi(spike_trains, indices=None):
-    """ computes the multi-variate isi distance profile for a set of spike
-    trains. That is the average isi-distance of all pairs of spike-trains:
-    S_isi(t) = 2/((N(N-1)) sum_{<i,j>} S_{isi}^{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 state: list or None
-    :returns: The averaged isi profile :math:`<S_{isi}>(t)`
-    :rtype: :class:`pyspike.function.PieceWiseConstFunc`
-    """
-    average_dist, M = _generic_profile_multi(spike_trains, isi_profile,
-                                             indices)
-    average_dist.mul_scalar(1.0/M)  # normalize
-    return average_dist
-
-
-############################################################
-# isi_distance_multi
-############################################################
-def isi_distance_multi(spike_trains, indices=None, interval=None):
-    """ computes the multi-variate isi-distance for a set of spike-trains.
-    That is the time average of the multi-variate spike profile:
-    I = \int_0^T 2/((N(N-1)) sum_{<i,j>} S_{isi}^{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)
-    :param interval: averaging interval given as a pair of floats, if None
-                     the average over the whole function is computed.
-    :type interval: Pair of floats or None.
-    :returns: The time-averaged isi distance :math:`I`
-    :rtype: double
-    """
-    return isi_profile_multi(spike_trains, indices).avrg(interval)
-
-
-############################################################
-# spike_profile_multi
-############################################################
-def spike_profile_multi(spike_trains, indices=None):
-    """ Computes the multi-variate spike distance profile for a set of spike
-    trains. That is the average spike-distance of all pairs of spike-trains:
-    :math:`S_spike(t) = 2/((N(N-1)) sum_{<i,j>} S_{spike}^{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_{spike}>(t)`
-    :rtype: :class:`pyspike.function.PieceWiseLinFunc`
-
-    """
-    average_dist, M = _generic_profile_multi(spike_trains, spike_profile,
-                                             indices)
-    average_dist.mul_scalar(1.0/M)  # normalize
-    return average_dist
-
-
-############################################################
-# spike_distance_multi
-############################################################
-def spike_distance_multi(spike_trains, indices=None, interval=None):
-    """ Computes the multi-variate spike distance for a set of spike trains.
-    That is the time average of the multi-variate spike profile:
-    S_{spike} = \int_0^T 2/((N(N-1)) sum_{<i,j>} S_{spike}^{i, j} dt
-    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
-    :param interval: averaging interval given as a pair of floats, if None
-                     the average over the whole function is computed.
-    :type interval: Pair of floats or None.
-    :returns: The averaged spike distance S.
-    :rtype: double
-    """
-    return spike_profile_multi(spike_trains, indices).avrg(interval)
-
-
-############################################################
-# spike_profile_multi
-############################################################
-def spike_sync_profile_multi(spike_trains, indices=None):
-    """ Computes the multi-variate spike synchronization profile for a set of
-    spike trains. For each spike in the set of spike trains, the multi-variate
-    profile is defined as the number of coincidences divided by the number of
-    spike trains pairs involving the spike train of containing this spike,
-    which is the number of spike trains minus one (N-1).
-
-    :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 multi-variate spike sync profile :math:`<S_{sync}>(t)`
-    :rtype: :class:`pyspike.function.DiscreteFunction`
-
-    """
-    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!
-    return average_dist
-
-
-############################################################
-# spike_distance_multi
-############################################################
-def spike_sync_multi(spike_trains, indices=None, interval=None):
-    """ Computes the multi-variate spike synchronization value for a set of
-    spike trains.
-
-    :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
-    :param interval: averaging interval given as a pair of floats, if None
-                     the average over the whole function is computed.
-    :type interval: Pair of floats or None.
-    :returns: The multi-variate spike synchronization value SYNC.
-    :rtype: double
-    """
-    return spike_sync_profile_multi(spike_trains, indices).avrg(interval)
-
-
-############################################################
-# generic_distance_matrix
-############################################################
-def _generic_distance_matrix(spike_trains, dist_function,
-                             indices=None, interval=None):
-    """ Internal implementation detail. Don't use this function directly.
-    Instead use isi_distance_matrix or spike_distance_matrix.
-    Computes the time averaged distance of all pairs of spike-trains.
-    Args:
-    - spike_trains: list of spike trains
-    - indices: list of indices defining which spike-trains to use
-    if None all given spike-trains are used (default=None)
-    Return:
-    - a 2D array of size len(indices)*len(indices) containing the average
-    pair-wise distance
-    """
-    if indices is None:
-        indices = np.arange(len(spike_trains))
-    indices = np.array(indices)
-    # check validity of indices
-    assert (indices < len(spike_trains)).all() and (indices >= 0).all(), \
-        "Invalid index list."
-    # generate a list of possible index pairs
-    pairs = [(indices[i], j) for i in range(len(indices)) for j in indices[i+1:]]
-
-    distance_matrix = np.zeros((len(indices), len(indices)))
-    for i, j in pairs:
-        d = dist_function(spike_trains[i], spike_trains[j], interval)
-        distance_matrix[i, j] = d
-        distance_matrix[j, i] = d
-    return distance_matrix
-
-
-############################################################
-# isi_distance_matrix
-############################################################
-def isi_distance_matrix(spike_trains, indices=None, interval=None):
-    """ Computes the time averaged isi-distance of all pairs of spike-trains.
-
-    :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
-    :param interval: averaging interval given as a pair of floats, if None
-                     the average over the whole function is computed.
-    :type interval: Pair of floats or None.
-    :returns: 2D array with the pair wise time average isi distances
-              :math:`I_{ij}`
-    :rtype: np.array
-    """
-    return _generic_distance_matrix(spike_trains, isi_distance,
-                                    indices, interval)
-
-
-############################################################
-# spike_distance_matrix
-############################################################
-def spike_distance_matrix(spike_trains, indices=None, interval=None):
-    """ Computes the time averaged spike-distance of all pairs of spike-trains.
-
-    :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
-    :param interval: averaging interval given as a pair of floats, if None
-                     the average over the whole function is computed.
-    :type interval: Pair of floats or None.
-    :returns: 2D array with the pair wise time average spike distances
-              :math:`S_{ij}`
-    :rtype: np.array
-    """
-    return _generic_distance_matrix(spike_trains, spike_distance,
-                                    indices, interval)
-
-
-############################################################
-# spike_sync_matrix
-############################################################
-def spike_sync_matrix(spike_trains, indices=None, interval=None):
-    """ Computes the overall spike-synchronization value of all pairs of
-    spike-trains.
-
-    :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
-    :param interval: averaging interval given as a pair of floats, if None
-                     the average over the whole function is computed.
-    :type interval: Pair of floats or None.
-    :returns: 2D array with the pair wise time spike synchronization values
-              :math:`SYNC_{ij}`
-    :rtype: np.array
-    """
-    return _generic_distance_matrix(spike_trains, spike_sync,
-                                    indices, interval)
diff --git a/pyspike/function.py b/pyspike/function.py
deleted file mode 100644
index 047c88a..0000000
--- a/pyspike/function.py
+++ /dev/null
@@ -1,585 +0,0 @@
-""" function.py
-
-Module containing classes representing piece-wise constant and piece-wise
-linear functions.
-
-Copyright 2014, Mario Mulansky <mario.mulansky@gmx.net>
-
-Distributed under the BSD License
-
-"""
-from __future__ import print_function
-
-import numpy as np
-import collections
-
-
-##############################################################
-# PieceWiseConstFunc
-##############################################################
-class PieceWiseConstFunc(object):
-    """ A class representing a piece-wise constant function. """
-
-    def __init__(self, x, y):
-        """ Constructs the piece-wise const function.
-
-        :param x: array of length N+1 defining the edges of the intervals of
-                  the pwc function.
-        :param y: array of length N defining the function values at the
-                  intervals.
-        """
-        # convert parameters to arrays, also ensures copying
-        self.x = np.array(x)
-        self.y = np.array(y)
-
-    def copy(self):
-        """ Returns a copy of itself
-
-        :rtype: :class:`PieceWiseConstFunc`
-        """
-        return PieceWiseConstFunc(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:`PieceWiseConstFunc`
-        :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 piece-wise function.
-
-        :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
-        y_plot[1::2] = self.y
-
-        return x_plot, y_plot
-
-    def integral(self, interval=None):
-        """ Returns the integral over the given interval.
-
-        :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
-            a = np.sum((self.x[1:]-self.x[:-1]) * self.y)
-        else:
-            # find the indices corresponding to the interval
-            start_ind = np.searchsorted(self.x, interval[0], side='right')
-            end_ind = np.searchsorted(self.x, interval[1], side='left')-1
-            assert start_ind > 0 and end_ind < len(self.x), \
-                "Invalid averaging interval"
-            # first the contribution from between the indices
-            a = np.sum((self.x[start_ind+1:end_ind+1] -
-                        self.x[start_ind:end_ind]) *
-                       self.y[start_ind:end_ind])
-            # correction from start to first index
-            a += (self.x[start_ind]-interval[0]) * self.y[start_ind-1]
-            # correction from last index to end
-            a += (interval[1]-self.x[end_ind]) * self.y[end_ind]
-        return a
-
-    def avrg(self, interval=None):
-        """ Computes the average of the piece-wise const function:
-        :math:`a = 1/T int_0^T f(x) dx` where T is the length of the interval.
-
-        :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
-            return self.integral() / (self.x[-1]-self.x[0])
-
-        # check if interval is as sequence
-        assert isinstance(interval, collections.Sequence), \
-            "Invalid value for `interval`. None, Sequence or Tuple expected."
-        # check if interval is a sequence of intervals
-        if not isinstance(interval[0], collections.Sequence):
-            # just one interval
-            a = self.integral(interval) / (interval[1]-interval[0])
-        else:
-            # several intervals
-            a = 0.0
-            int_length = 0.0
-            for ival in interval:
-                a += self.integral(ival)
-                int_length += ival[1] - ival[0]
-            a /= int_length
-        return a
-
-    def add(self, f):
-        """ Adds another PieceWiseConst function to this function.
-        Note: only functions defined on the same interval can be summed.
-
-        :param f: :class:`PieceWiseConstFunc` 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_piece_wise_const_cython as \
-                add_piece_wise_const_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_piece_wise_const_python as \
-                add_piece_wise_const_impl
-
-        self.x, self.y = add_piece_wise_const_impl(self.x, self.y, f.x, f.y)
-
-    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:
-    """ A class representing a piece-wise linear function. """
-
-    def __init__(self, x, y1, y2):
-        """ Constructs the piece-wise linear function.
-
-        :param x: array of length N+1 defining the edges of the intervals of
-                  the pwc function.
-        :param y1: array of length N defining the function values at the left
-                  of the intervals.
-        :param y2: array of length N defining the function values at the right
-                  of the intervals.
-        """
-        # convert to array, which also ensures copying
-        self.x = np.array(x)
-        self.y1 = np.array(y1)
-        self.y2 = np.array(y2)
-
-    def copy(self):
-        """ Returns a copy of itself
-
-        :rtype: :class:`PieceWiseLinFunc`
-        """
-        return PieceWiseLinFunc(self.x, self.y1, self.y2)
-
-    def almost_equal(self, other, decimal=14):
-        """ Checks if the function is equal to another function up to `decimal`
-        precision.
-
-        :param other: another :class:`PieceWiseLinFunc`
-        :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.y1, other.y1, atol=eps, rtol=0.0) and \
-            np.allclose(self.y2, other.y2, atol=eps, rtol=0.0)
-
-    def get_plottable_data(self):
-        """ Returns two arrays containing x- and y-coordinates for immeditate
-        plotting of the piece-wise function.
-
-        :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_like(x_plot)
-        y_plot[0::2] = self.y1
-        y_plot[1::2] = self.y2
-        return x_plot, y_plot
-
-    def integral(self, interval=None):
-        """ Returns the integral over the given interval.
-
-        :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
-        """
-
-        def intermediate_value(x0, x1, y0, y1, x):
-            """ computes the intermediate value of a linear function """
-            return y0 + (y1-y0)*(x-x0)/(x1-x0)
-
-        if interval is None:
-            # no interval given, integrate over the whole spike train
-            integral = np.sum((self.x[1:]-self.x[:-1]) * 0.5*(self.y1+self.y2))
-        else:
-            # find the indices corresponding to the interval
-            start_ind = np.searchsorted(self.x, interval[0], side='right')
-            end_ind = np.searchsorted(self.x, interval[1], side='left')-1
-            assert start_ind > 0 and end_ind < len(self.x), \
-                "Invalid averaging interval"
-            # first the contribution from between the indices
-            integral = np.sum((self.x[start_ind+1:end_ind+1] -
-                               self.x[start_ind:end_ind]) *
-                              0.5*(self.y1[start_ind:end_ind] +
-                                   self.y2[start_ind:end_ind]))
-            # correction from start to first index
-            integral += (self.x[start_ind]-interval[0]) * 0.5 * \
-                        (self.y2[start_ind-1] +
-                         intermediate_value(self.x[start_ind-1],
-                                            self.x[start_ind],
-                                            self.y1[start_ind-1],
-                                            self.y2[start_ind-1],
-                                            interval[0]
-                                            ))
-            # correction from last index to end
-            integral += (interval[1]-self.x[end_ind]) * 0.5 * \
-                        (self.y1[end_ind] +
-                         intermediate_value(self.x[end_ind], self.x[end_ind+1],
-                                            self.y1[end_ind], self.y2[end_ind],
-                                            interval[1]
-                                            ))
-        return integral
-
-    def avrg(self, interval=None):
-        """ Computes the average of the piece-wise linear function:
-        :math:`a = 1/T int_0^T f(x) dx` where T is the length of the interval.
-
-        :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
-            return self.integral() / (self.x[-1]-self.x[0])
-
-        # check if interval is as sequence
-        assert isinstance(interval, collections.Sequence), \
-            "Invalid value for `interval`. None, Sequence or Tuple expected."
-        # check if interval is a sequence of intervals
-        if not isinstance(interval[0], collections.Sequence):
-            # just one interval
-            a = self.integral(interval) / (interval[1]-interval[0])
-        else:
-            # several intervals
-            a = 0.0
-            int_length = 0.0
-            for ival in interval:
-                a += self.integral(ival)
-                int_length += ival[1] - ival[0]
-            a /= int_length
-        return a
-
-    def add(self, f):
-        """ Adds another PieceWiseLin function to this function.
-        Note: only functions defined on the same interval can be summed.
-
-        :param f: :class:`PieceWiseLinFunc` 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"
-
-        # python implementation
-        # from python_backend import add_piece_wise_lin_python
-        # self.x, self.y1, self.y2 = add_piece_wise_lin_python(
-        #     self.x, self.y1, self.y2, f.x, f.y1, f.y2)
-
-        # cython version
-        try:
-            from cython_add import add_piece_wise_lin_cython as \
-                add_piece_wise_lin_impl
-        except ImportError:
-            print("Warning: add_piece_wise_lin_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_piece_wise_lin_python as \
-                add_piece_wise_lin_impl
-
-        self.x, self.y1, self.y2 = add_piece_wise_lin_impl(
-            self.x, self.y1, self.y2, f.x, f.y1, f.y2)
-
-    def mul_scalar(self, fac):
-        """ Multiplies the function with a scalar value
-
-        :param fac: Value to multiply
-        :type fac: double
-        :rtype: None
-        """
-        self.y1 *= fac
-        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, averaging_window_size=0):
-        """ Returns two arrays containing x- and y-coordinates for plotting
-        the interval sequence. The optional parameter `averaging_window_size`
-        determines the size of an averaging window to smoothen the profile. If
-        this value is 0, no averaging is performed.
-
-        :param averaging_window_size: size of the averaging window, default=0.
-        :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 averaging_window_size > 0:
-            # for the averaged profile we have to take the multiplicity into
-            # account. values with higher multiplicity should be consider as if
-            # they appeared several times. Hence we can not know how many
-            # entries we have to consider to the left and right. Rather, we
-            # will iterate until some wanted multiplicity is reached.
-
-            # the first value in self.mp contains the number of averaged
-            # profiles without any possible extra multiplicities
-            # (by implementation)
-            expected_mp = (averaging_window_size+1) * int(self.mp[0])
-            y_plot = np.zeros_like(self.y)
-            # compute the values in a loop, could be done in cython if required
-            for i in xrange(len(y_plot)):
-
-                if self.mp[i] >= expected_mp:
-                    # the current value contains already all the wanted
-                    # multiplicity
-                    y_plot[i] = self.y[i]/self.mp[i]
-                    continue
-
-                # first look to the right
-                y = self.y[i]
-                mp_r = self.mp[i]
-                j = i+1
-                while j < len(y_plot):
-                    if mp_r+self.mp[j] < expected_mp:
-                        # if we still dont reach the required multiplicity
-                        # we take the whole value
-                        y += self.y[j]
-                        mp_r += self.mp[j]
-                    else:
-                        # otherwise, just some fraction
-                        y += self.y[j] * (expected_mp - mp_r)/self.mp[j]
-                        mp_r += (expected_mp - mp_r)
-                        break
-                    j += 1
-
-                # same story to the left
-                mp_l = self.mp[i]
-                j = i-1
-                while j >= 0:
-                    if mp_l+self.mp[j] < expected_mp:
-                        y += self.y[j]
-                        mp_l += self.mp[j]
-                    else:
-                        y += self.y[j] * (expected_mp - mp_l)/self.mp[j]
-                        mp_l += (expected_mp - mp_l)
-                        break
-                    j -= 1
-                y_plot[i] = y/(mp_l+mp_r-self.mp[i])
-            return 1.0*self.x, y_plot
-
-        else:  # k = 0
-
-            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 discrete
-        function, this amounts to the sum over all values divided by the total
-        multiplicity.
-
-        :param interval: integration interval given as a pair of floats, or a
-                         sequence of pairs in case of multiple intervals, if
-                         None the integral over the whole function is computed.
-        :type interval: Pair, sequence of pairs, or None.
-        :returns: the integral
-        :rtype: float
-        """
-
-        def get_indices(ival):
-            """ Retuns the indeces surrounding the given interval"""
-            start_ind = np.searchsorted(self.x, ival[0], side='right')
-            end_ind = np.searchsorted(self.x, ival[1], side='left')
-            assert start_ind > 0 and end_ind < len(self.x), \
-                "Invalid averaging interval"
-            return start_ind, end_ind
-
-        if interval is None:
-            # no interval given, integrate over the whole spike train
-            # don't count the first value, which is zero by definition
-            return 1.0 * np.sum(self.y[1:-1]) / np.sum(self.mp[1:-1])
-
-        # check if interval is as sequence
-        assert isinstance(interval, collections.Sequence), \
-            "Invalid value for `interval`. None, Sequence or Tuple expected."
-        # check if interval is a sequence of intervals
-        if not isinstance(interval[0], collections.Sequence):
-            # find the indices corresponding to the interval
-            start_ind, end_ind = get_indices(interval)
-            return (np.sum(self.y[start_ind:end_ind]) /
-                    np.sum(self.mp[start_ind:end_ind]))
-        else:
-            value = 0.0
-            multiplicity = 0.0
-            for ival in interval:
-                # find the indices corresponding to the interval
-                start_ind, end_ind = get_indices(ival)
-                value += np.sum(self.y[start_ind:end_ind])
-                multiplicity += np.sum(self.mp[start_ind:end_ind])
-        return value/multiplicity
-
-    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
-        """
-        return self.integral(interval)
-
-    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.
-
-    :param profiles: list of :class:`PieceWiseConstFunc` or
-                     :class:`PieceWiseLinFunc` representing ISI- or
-                     SPIKE-profiles to be averaged.
-    :returns: the averages profile :math:`<S_{isi}>` or :math:`<S_{spike}>`.
-    :rtype: :class:`PieceWiseConstFunc` or :class:`PieceWiseLinFunc`
-    """
-    assert len(profiles) > 1
-
-    avrg_profile = profiles[0].copy()
-    for i in xrange(1, len(profiles)):
-        avrg_profile.add(profiles[i])
-    avrg_profile.mul_scalar(1.0/len(profiles))  # normalize
-
-    return avrg_profile
diff --git a/pyspike/generic.py b/pyspike/generic.py
new file mode 100644
index 0000000..4f278d2
--- /dev/null
+++ b/pyspike/generic.py
@@ -0,0 +1,83 @@
+"""
+
+Generic functions to compute multi-variate profiles and distance matrices.
+
+Copyright 2015, Mario Mulansky <mario.mulansky@gmx.net>
+
+Distributed under the BSD License
+"""
+
+
+import numpy as np
+
+
+############################################################
+# _generic_profile_multi
+############################################################
+def _generic_profile_multi(spike_trains, pair_distance_func, indices=None):
+    """ Internal implementation detail, don't call this function directly,
+    use isi_profile_multi or spike_profile_multi instead.
+
+    Computes the multi-variate distance for a set of spike-trains using the
+    pair_dist_func to compute pair-wise distances. That is it computes the
+    average distance of all pairs of spike-trains:
+    :math:`S(t) = 2/((N(N-1)) sum_{<i,j>} S_{i,j}`,
+    where the sum goes over all pairs <i,j>.
+    Args:
+    - spike_trains: list of spike trains
+    - pair_distance_func: function computing the distance of two spike trains
+    - indices: list of indices defining which spike trains to use,
+    if None all given spike trains are used (default=None)
+    Returns:
+    - The averaged multi-variate distance of all pairs
+    """
+    if indices is None:
+        indices = np.arange(len(spike_trains))
+    indices = np.array(indices)
+    # check validity of indices
+    assert (indices < len(spike_trains)).all() and (indices >= 0).all(), \
+        "Invalid index list."
+    # generate a list of possible index pairs
+    pairs = [(indices[i], j) for i in range(len(indices))
+             for j in indices[i+1:]]
+    # start with first pair
+    (i, j) = pairs[0]
+    average_dist = pair_distance_func(spike_trains[i], spike_trains[j])
+    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
+    return average_dist, len(pairs)
+
+
+############################################################
+# generic_distance_matrix
+############################################################
+def _generic_distance_matrix(spike_trains, dist_function,
+                             indices=None, interval=None):
+    """ Internal implementation detail. Don't use this function directly.
+    Instead use isi_distance_matrix or spike_distance_matrix.
+    Computes the time averaged distance of all pairs of spike-trains.
+    Args:
+    - spike_trains: list of spike trains
+    - indices: list of indices defining which spike-trains to use
+    if None all given spike-trains are used (default=None)
+    Return:
+    - a 2D array of size len(indices)*len(indices) containing the average
+    pair-wise distance
+    """
+    if indices is None:
+        indices = np.arange(len(spike_trains))
+    indices = np.array(indices)
+    # check validity of indices
+    assert (indices < len(spike_trains)).all() and (indices >= 0).all(), \
+        "Invalid index list."
+    # generate a list of possible index pairs
+    pairs = [(indices[i], j) for i in range(len(indices))
+             for j in indices[i+1:]]
+
+    distance_matrix = np.zeros((len(indices), len(indices)))
+    for i, j in pairs:
+        d = dist_function(spike_trains[i], spike_trains[j], interval)
+        distance_matrix[i, j] = d
+        distance_matrix[j, i] = d
+    return distance_matrix
diff --git a/pyspike/isi_distance.py b/pyspike/isi_distance.py
new file mode 100644
index 0000000..c2ef8e8
--- /dev/null
+++ b/pyspike/isi_distance.py
@@ -0,0 +1,134 @@
+"""
+
+Module containing several functions to compute the ISI profiles and distances
+
+Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net>
+
+Distributed under the BSD License
+"""
+
+from pyspike import PieceWiseConstFunc
+from pyspike.generic import _generic_profile_multi, _generic_distance_matrix
+
+
+############################################################
+# isi_profile
+############################################################
+def isi_profile(spikes1, spikes2):
+    """ Computes the isi-distance profile :math:`S_{isi}(t)` of the two given
+    spike trains. Retruns the profile as a PieceWiseConstFunc object. The S_isi
+    values are defined positive S_isi(t)>=0.  The spike trains are expected
+    to have auxiliary spikes at the beginning and end of the interval. Use the
+    function add_auxiliary_spikes to add those spikes to the spike train.
+
+    :param spikes1: ordered array of spike times with auxiliary spikes.
+    :param spikes2: ordered array of spike times with auxiliary spikes.
+    :returns: The isi-distance profile :math:`S_{isi}(t)`
+    :rtype: :class:`pyspike.function.PieceWiseConstFunc`
+
+    """
+    # check for auxiliary spikes - first and last spikes should be identical
+    assert spikes1[0] == spikes2[0], \
+        "Given spike trains seems not to have auxiliary spikes!"
+    assert spikes1[-1] == spikes2[-1], \
+        "Given spike trains seems not to have auxiliary spikes!"
+
+    # load cython implementation
+    try:
+        from cython.cython_distance import isi_distance_cython \
+            as isi_distance_impl
+    except ImportError:
+        print("Warning: isi_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 cython.python_backend import isi_distance_python \
+            as isi_distance_impl
+
+    times, values = isi_distance_impl(spikes1, spikes2)
+    return PieceWiseConstFunc(times, values)
+
+
+############################################################
+# isi_distance
+############################################################
+def isi_distance(spikes1, spikes2, interval=None):
+    """ Computes the isi-distance I of the given spike trains. The
+    isi-distance is the integral over the isi distance profile
+    :math:`S_{isi}(t)`:
+
+    .. math:: I = \int_{T_0}^{T_1} S_{isi}(t) dt.
+
+    :param spikes1: ordered array of spike times with auxiliary spikes.
+    :param spikes2: ordered array of spike times with auxiliary spikes.
+    :param interval: averaging interval given as a pair of floats (T0, T1),
+                     if None the average over the whole function is computed.
+    :type interval: Pair of floats or None.
+    :returns: The isi-distance I.
+    :rtype: double
+    """
+    return isi_profile(spikes1, spikes2).avrg(interval)
+
+
+############################################################
+# isi_profile_multi
+############################################################
+def isi_profile_multi(spike_trains, indices=None):
+    """ computes the multi-variate isi distance profile for a set of spike
+    trains. That is the average isi-distance of all pairs of spike-trains:
+    S_isi(t) = 2/((N(N-1)) sum_{<i,j>} S_{isi}^{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 state: list or None
+    :returns: The averaged isi profile :math:`<S_{isi}>(t)`
+    :rtype: :class:`pyspike.function.PieceWiseConstFunc`
+    """
+    average_dist, M = _generic_profile_multi(spike_trains, isi_profile,
+                                             indices)
+    average_dist.mul_scalar(1.0/M)  # normalize
+    return average_dist
+
+
+############################################################
+# isi_distance_multi
+############################################################
+def isi_distance_multi(spike_trains, indices=None, interval=None):
+    """ computes the multi-variate isi-distance for a set of spike-trains.
+    That is the time average of the multi-variate spike profile:
+    I = \int_0^T 2/((N(N-1)) sum_{<i,j>} S_{isi}^{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)
+    :param interval: averaging interval given as a pair of floats, if None
+                     the average over the whole function is computed.
+    :type interval: Pair of floats or None.
+    :returns: The time-averaged isi distance :math:`I`
+    :rtype: double
+    """
+    return isi_profile_multi(spike_trains, indices).avrg(interval)
+
+
+############################################################
+# isi_distance_matrix
+############################################################
+def isi_distance_matrix(spike_trains, indices=None, interval=None):
+    """ Computes the time averaged isi-distance of all pairs of spike-trains.
+
+    :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
+    :param interval: averaging interval given as a pair of floats, if None
+                     the average over the whole function is computed.
+    :type interval: Pair of floats or None.
+    :returns: 2D array with the pair wise time average isi distances
+              :math:`I_{ij}`
+    :rtype: np.array
+    """
+    return _generic_distance_matrix(spike_trains, isi_distance,
+                                    indices, interval)
diff --git a/pyspike/spike_distance.py b/pyspike/spike_distance.py
new file mode 100644
index 0000000..f721c86
--- /dev/null
+++ b/pyspike/spike_distance.py
@@ -0,0 +1,136 @@
+"""
+
+Module containing several functions to compute SPIKE profiles and distances
+
+Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net>
+
+Distributed under the BSD License
+"""
+
+from pyspike import PieceWiseLinFunc
+from pyspike.generic import _generic_profile_multi, _generic_distance_matrix
+
+
+############################################################
+# spike_profile
+############################################################
+def spike_profile(spikes1, spikes2):
+    """ Computes the spike-distance profile S_spike(t) of the two given spike
+    trains. Returns the profile as a PieceWiseLinFunc object. The S_spike
+    values are defined positive S_spike(t)>=0. The spike trains are expected to
+    have auxiliary spikes at the beginning and end of the interval. Use the
+    function add_auxiliary_spikes to add those spikes to the spike train.
+
+    :param spikes1: ordered array of spike times with auxiliary spikes.
+    :param spikes2: ordered array of spike times with auxiliary spikes.
+    :returns: The spike-distance profile :math:`S_{spike}(t)`.
+    :rtype: :class:`pyspike.function.PieceWiseLinFunc`
+
+    """
+    # check for auxiliary spikes - first and last spikes should be identical
+    assert spikes1[0] == spikes2[0], \
+        "Given spike trains seems not to have auxiliary spikes!"
+    assert spikes1[-1] == spikes2[-1], \
+        "Given spike trains seems not to have auxiliary spikes!"
+
+    # cython implementation
+    try:
+        from cython.cython_distance import spike_distance_cython \
+            as spike_distance_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 cython.python_backend import spike_distance_python \
+            as spike_distance_impl
+
+    times, y_starts, y_ends = spike_distance_impl(spikes1, spikes2)
+    return PieceWiseLinFunc(times, y_starts, y_ends)
+
+
+############################################################
+# spike_distance
+############################################################
+def spike_distance(spikes1, spikes2, interval=None):
+    """ Computes the spike-distance S of the given spike trains. The
+    spike-distance is the integral over the isi distance profile S_spike(t):
+
+    .. math:: S = \int_{T_0}^{T_1} S_{spike}(t) dt.
+
+    :param spikes1: ordered array of spike times with auxiliary spikes.
+    :param spikes2: ordered array of spike times with auxiliary spikes.
+    :param interval: averaging interval given as a pair of floats (T0, T1),
+                     if None the average over the whole function is computed.
+    :type interval: Pair of floats or None.
+    :returns: The spike-distance.
+    :rtype: double
+
+    """
+    return spike_profile(spikes1, spikes2).avrg(interval)
+
+
+############################################################
+# spike_profile_multi
+############################################################
+def spike_profile_multi(spike_trains, indices=None):
+    """ Computes the multi-variate spike distance profile for a set of spike
+    trains. That is the average spike-distance of all pairs of spike-trains:
+    :math:`S_spike(t) = 2/((N(N-1)) sum_{<i,j>} S_{spike}^{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_{spike}>(t)`
+    :rtype: :class:`pyspike.function.PieceWiseLinFunc`
+
+    """
+    average_dist, M = _generic_profile_multi(spike_trains, spike_profile,
+                                             indices)
+    average_dist.mul_scalar(1.0/M)  # normalize
+    return average_dist
+
+
+############################################################
+# spike_distance_multi
+############################################################
+def spike_distance_multi(spike_trains, indices=None, interval=None):
+    """ Computes the multi-variate spike distance for a set of spike trains.
+    That is the time average of the multi-variate spike profile:
+    S_{spike} = \int_0^T 2/((N(N-1)) sum_{<i,j>} S_{spike}^{i, j} dt
+    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
+    :param interval: averaging interval given as a pair of floats, if None
+                     the average over the whole function is computed.
+    :type interval: Pair of floats or None.
+    :returns: The averaged spike distance S.
+    :rtype: double
+    """
+    return spike_profile_multi(spike_trains, indices).avrg(interval)
+
+
+############################################################
+# spike_distance_matrix
+############################################################
+def spike_distance_matrix(spike_trains, indices=None, interval=None):
+    """ Computes the time averaged spike-distance of all pairs of spike-trains.
+
+    :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
+    :param interval: averaging interval given as a pair of floats, if None
+                     the average over the whole function is computed.
+    :type interval: Pair of floats or None.
+    :returns: 2D array with the pair wise time average spike distances
+              :math:`S_{ij}`
+    :rtype: np.array
+    """
+    return _generic_distance_matrix(spike_trains, spike_distance,
+                                    indices, interval)
diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py
new file mode 100644
index 0000000..342bf69
--- /dev/null
+++ b/pyspike/spike_sync.py
@@ -0,0 +1,136 @@
+"""
+
+Module containing several functions to compute SPIKE-Synchronization profiles
+and distances
+
+Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net>
+
+Distributed under the BSD License
+"""
+
+from functools import partial
+from pyspike import DiscreteFunc
+from pyspike.generic import _generic_profile_multi, _generic_distance_matrix
+
+
+############################################################
+# spike_sync_profile
+############################################################
+def spike_sync_profile(spikes1, spikes2):
+    """ Computes the spike-synchronization profile S_sync(t) of the two given
+    spike trains. Returns the profile as a DiscreteFunction object. The S_sync
+    values are either 1 or 0, indicating the presence or absence of a
+    coincidence. The spike trains are expected to have auxiliary spikes at the
+    beginning and end of the interval. Use the function add_auxiliary_spikes to
+    add those spikes to the spike train.
+
+    :param spikes1: ordered array of spike times with auxiliary spikes.
+    :param spikes2: ordered array of spike times with auxiliary spikes.
+    :returns: The spike-distance profile :math:`S_{sync}(t)`.
+    :rtype: :class:`pyspike.function.DiscreteFunction`
+
+    """
+
+    # cython implementation
+    try:
+        from cython.cython_distance import coincidence_cython \
+            as coincidence_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 cython.python_backend import coincidence_python \
+            as coincidence_impl
+
+    times, coincidences, multiplicity = coincidence_impl(spikes1, spikes2)
+
+    return DiscreteFunc(times, coincidences, multiplicity)
+
+
+############################################################
+# spike_sync
+############################################################
+def spike_sync(spikes1, spikes2, interval=None):
+    """ Computes the spike synchronization value SYNC of the given spike
+    trains. The spike synchronization value is the computed as the total number
+    of coincidences divided by the total number of spikes:
+
+    .. math:: SYNC = \sum_n C_n / N.
+
+    :param spikes1: ordered array of spike times with auxiliary spikes.
+    :param spikes2: ordered array of spike times with auxiliary spikes.
+    :param interval: averaging interval given as a pair of floats (T0, T1),
+                     if None the average over the whole function is computed.
+    :type interval: Pair of floats or None.
+    :returns: The spike synchronization value.
+    :rtype: double
+    """
+    return spike_sync_profile(spikes1, spikes2).avrg(interval)
+
+
+############################################################
+# spike_sync_profile_multi
+############################################################
+def spike_sync_profile_multi(spike_trains, indices=None):
+    """ Computes the multi-variate spike synchronization profile for a set of
+    spike trains. For each spike in the set of spike trains, the multi-variate
+    profile is defined as the number of coincidences divided by the number of
+    spike trains pairs involving the spike train of containing this spike,
+    which is the number of spike trains minus one (N-1).
+
+    :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 multi-variate spike sync profile :math:`<S_{sync}>(t)`
+    :rtype: :class:`pyspike.function.DiscreteFunction`
+
+    """
+    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!
+    return average_dist
+
+
+############################################################
+# spike_sync_multi
+############################################################
+def spike_sync_multi(spike_trains, indices=None, interval=None):
+    """ Computes the multi-variate spike synchronization value for a set of
+    spike trains.
+
+    :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
+    :param interval: averaging interval given as a pair of floats, if None
+                     the average over the whole function is computed.
+    :type interval: Pair of floats or None.
+    :returns: The multi-variate spike synchronization value SYNC.
+    :rtype: double
+    """
+    return spike_sync_profile_multi(spike_trains, indices).avrg(interval)
+
+
+############################################################
+# spike_sync_matrix
+############################################################
+def spike_sync_matrix(spike_trains, indices=None, interval=None):
+    """ Computes the overall spike-synchronization value of all pairs of
+    spike-trains.
+
+    :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
+    :param interval: averaging interval given as a pair of floats, if None
+                     the average over the whole function is computed.
+    :type interval: Pair of floats or None.
+    :returns: 2D array with the pair wise time spike synchronization values
+              :math:`SYNC_{ij}`
+    :rtype: np.array
+    """
+    return _generic_distance_matrix(spike_trains, spike_sync,
+                                    indices, interval)
diff --git a/setup.py b/setup.py
index 9cab5de..289d521 100644
--- a/setup.py
+++ b/setup.py
@@ -21,8 +21,8 @@ except ImportError:
 else:
     use_cython = True
 
-if os.path.isfile("pyspike/cython_add.c") and \
-   os.path.isfile("pyspike/cython_distance.c"):
+if os.path.isfile("pyspike/cython/cython_add.c") and \
+   os.path.isfile("pyspike/cython/cython_distance.c"):
     use_c = True
 else:
     use_c = False
@@ -32,14 +32,14 @@ ext_modules = []
 
 if use_cython:  # Cython is available, compile .pyx -> .c
     ext_modules += [
-        Extension("pyspike.cython_add", ["pyspike/cython_add.pyx"]),
-        Extension("pyspike.cython_distance", ["pyspike/cython_distance.pyx"]),
+        Extension("pyspike.cython.cython_add", ["pyspike/cython/cython_add.pyx"]),
+        Extension("pyspike.cython.cython_distance", ["pyspike/cython/cython_distance.pyx"]),
     ]
     cmdclass.update({'build_ext': build_ext})
 elif use_c:  # c files are there, compile to binaries
     ext_modules += [
-        Extension("pyspike.cython_add", ["pyspike/cython_add.c"]),
-        Extension("pyspike.cython_distance", ["pyspike/cython_distance.c"]),
+        Extension("pyspike.cython.cython_add", ["pyspike/cython/cython_add.c"]),
+        Extension("pyspike.cython.cython_distance", ["pyspike/cython/cython_distance.c"]),
     ]
 # neither cython nor c files available -> automatic fall-back to python backend
 
@@ -78,7 +78,7 @@ train similarity',
         'Programming Language :: Python :: 2.7',
     ],
     package_data={
-        'pyspike': ['cython_add.c', 'cython_distance.c'],
+        'pyspike': ['cython/cython_add.c', 'cython/cython_distance.c'],
         'test': ['Spike_testdata.txt']
     }
 )
diff --git a/test/test_function.py b/test/test_function.py
index 933fd2e..d81b03a 100644
--- a/test/test_function.py
+++ b/test/test_function.py
@@ -92,11 +92,12 @@ def test_pwc_avrg():
     y = [0.5, 1.0, -0.25, 0.0, 1.5]
     f2 = spk.PieceWiseConstFunc(x, y)
 
-    f_avrg = spk.average_profile([f1, f2])
+    f1.add(f2)
+    f1.mul_scalar(0.5)
     x_expected = [0.0, 0.75, 1.0, 2.0, 2.5, 2.7, 4.0]
     y_expected = [0.75, 1.0, 0.25, 0.625, 0.375, 1.125]
-    assert_array_almost_equal(f_avrg.x, x_expected, decimal=16)
-    assert_array_almost_equal(f_avrg.y, y_expected, decimal=16)
+    assert_array_almost_equal(f1.x, x_expected, decimal=16)
+    assert_array_almost_equal(f1.y, y_expected, decimal=16)
 
 
 def test_pwl():
@@ -196,11 +197,12 @@ def test_pwl_avrg():
     y2_expected = np.array([0.8+1.0+0.5*0.75, 1.5+1.0-0.8*0.25/1.25, -0.4+0.2,
                             1.5-1.0, 0.75-0.5*0.2/1.5, 2.25]) / 2
 
-    f_avrg = spk.average_profile([f1, f2])
+    f1.add(f2)
+    f1.mul_scalar(0.5)
 
-    assert_array_almost_equal(f_avrg.x, x_expected, decimal=16)
-    assert_array_almost_equal(f_avrg.y1, y1_expected, decimal=16)
-    assert_array_almost_equal(f_avrg.y2, y2_expected, decimal=16)
+    assert_array_almost_equal(f1.x, x_expected, decimal=16)
+    assert_array_almost_equal(f1.y1, y1_expected, decimal=16)
+    assert_array_almost_equal(f1.y2, y2_expected, decimal=16)
 
 
 def test_df():
@@ -208,7 +210,7 @@ def test_df():
     x = [0.0, 1.0, 2.0, 2.5, 4.0]
     y = [0.0, 1.0, 1.0, 0.0, 1.0]
     mp = [1.0, 2.0, 1.0, 2.0, 1.0]
-    f = spk.DiscreteFunction(x, y, mp)
+    f = spk.DiscreteFunc(x, y, mp)
     xp, yp = f.get_plottable_data()
 
     xp_expected = [0.0, 1.0, 2.0, 2.5, 4.0]
@@ -237,6 +239,9 @@ if __name__ == "__main__":
     test_pwc()
     test_pwc_add()
     test_pwc_mul()
+    test_pwc_avrg()
     test_pwl()
     test_pwl_add()
     test_pwl_mul()
+    test_pwl_avrg()
+    test_df()