+ multivariate distances

3 files changed, 112 insertions, 7 deletions

diff --git a/pyspike/__init__.py b/pyspike/__init__.py
index 1784037..2143bdc 100644
--- a/pyspike/__init__.py
+++ b/pyspike/__init__.py
@@ -1,5 +1,6 @@
 __all__ = ["function", "distances", "spikes"]
 
 from function import PieceWiseConstFunc, PieceWiseLinFunc
-from distances import add_auxiliary_spikes, isi_distance, spike_distance
+from distances import add_auxiliary_spikes, isi_distance, spike_distance, \
+    isi_distance_multi, spike_distance_multi
 from spikes import spike_train_from_string, merge_spike_trains
diff --git a/pyspike/distances.py b/pyspike/distances.py
index f4be625..52c6640 100644
--- a/pyspike/distances.py
+++ b/pyspike/distances.py
@@ -224,3 +224,79 @@ def spike_distance(spikes1, spikes2):
     # could be less than original length due to equal spike times
     return PieceWiseLinFunc(spike_events[:index+1], 
                             y_starts[:index], y_ends[:index])
+
+
+
+
+############################################################
+# multi_distance
+############################################################
+def multi_distance(spike_trains, pair_distance_func, indices=None):
+    """ Internal implementation detail, use isi_distance_multi or 
+    spike_distance_multi.
+
+    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:
+    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==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 = [(i,j) for i in 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
+    average_dist.mul_scalar(1.0/len(pairs)) # normalize
+    return average_dist
+
+
+############################################################
+# isi_distance_multi
+############################################################
+def isi_distance_multi(spike_trains, indices=None):
+    """ computes the multi-variate isi-distance for a set of spike-trains. That
+    is the average isi-distance of all pairs of spike-trains:
+    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
+    - indices: list of indices defining which spike trains to use, 
+    if None all given spike trains are used (default=None)
+    Returns:
+    - A PieceWiseConstFunc representing the averaged isi distance S
+    """
+    return multi_distance(spike_trains, isi_distance, indices)
+
+
+############################################################
+# spike_distance_multi
+############################################################
+def spike_distance_multi(spike_trains, indices=None):
+    """ computes the multi-variate spike-distance for a set of spike-trains. 
+    That is the average spike-distance of all pairs of spike-trains:
+    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
+    - indices: list of indices defining which spike trains to use, 
+    if None all given spike trains are used (default=None)
+    Returns:
+    - A PieceWiseLinFunc representing the averaged spike distance S
+    """
+    return multi_distance(spike_trains, spike_distance, indices)
diff --git a/pyspike/function.py b/pyspike/function.py
index 3a5a01c..26ca4b2 100644
--- a/pyspike/function.py
+++ b/pyspike/function.py
@@ -10,6 +10,7 @@ from __future__ import print_function
 
 import numpy as np
 
+
 ##############################################################
 # PieceWiseConstFunc
 ##############################################################
@@ -18,7 +19,7 @@ class PieceWiseConstFunc:
     
     def __init__(self, x, y):
         """ Constructs the piece-wise const function.
-        Params:
+        Args:
         - x: array of length N+1 defining the edges of the intervals of the pwc
         function.
         - y: array of length N defining the function values at the intervals.
@@ -26,6 +27,19 @@ class PieceWiseConstFunc:
         self.x = np.array(x)
         self.y = np.array(y)
 
+    def almost_equal(self, other, decimal=14):
+        """ Checks if the function is equal to another function up to `decimal`
+        precision.
+        Args:
+        - other: another PieceWiseConstFunc object
+        Returns:
+        True if the two functions are equal up to `decimal` decimals, 
+        False otherwise
+        """
+        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.
@@ -63,7 +77,7 @@ class PieceWiseConstFunc:
     def add(self, f):
         """ Adds another PieceWiseConst function to this function. 
         Note: only functions defined on the same interval can be summed.
-        Params:
+        Args:
         - f: PieceWiseConst function to be added.
         """
         assert self.x[0] == f.x[0], "The functions have different intervals"
@@ -111,7 +125,7 @@ class PieceWiseConstFunc:
 
     def mul_scalar(self, fac):
         """ Multiplies the function with a scalar value
-        Params:
+        Args:
         - fac: Value to multiply
         """
         self.y *= fac
@@ -125,7 +139,7 @@ class PieceWiseLinFunc:
     
     def __init__(self, x, y1, y2):
         """ Constructs the piece-wise linear function.
-        Params:
+        Args:
         - x: array of length N+1 defining the edges of the intervals of the pwc
         function.
         - y1: array of length N defining the function values at the left of the 
@@ -137,6 +151,20 @@ class PieceWiseLinFunc:
         self.y1 = np.array(y1)
         self.y2 = np.array(y2)
 
+    def almost_equal(self, other, decimal=14):
+        """ Checks if the function is equal to another function up to `decimal`
+        precision.
+        Args:
+        - other: another PieceWiseLinFunc object
+        Returns:
+        True if the two functions are equal up to `decimal` decimals, 
+        False otherwise
+        """
+        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.
@@ -171,7 +199,7 @@ class PieceWiseLinFunc:
     def add(self, f):
         """ Adds another PieceWiseLin function to this function. 
         Note: only functions defined on the same interval can be summed.
-        Params:
+        Args:
         - f: PieceWiseLin function to be added.
         """
         assert self.x[0] == f.x[0], "The functions have different intervals"
@@ -246,7 +274,7 @@ class PieceWiseLinFunc:
 
     def mul_scalar(self, fac):
         """ Multiplies the function with a scalar value
-        Params:
+        Args:
         - fac: Value to multiply
         """
         self.y1 *= fac