performance improvements

use recursive approach to compute average profile
for average multivariate distances, dont compute average multivariate profile,
but average distances directly.

4 files changed, 86 insertions, 23 deletions

diff --git a/pyspike/cython/cython_add.pyx b/pyspike/cython/cython_add.pyx
index ac64005..8da1e53 100644
--- a/pyspike/cython/cython_add.pyx
+++ b/pyspike/cython/cython_add.pyx
@@ -83,13 +83,9 @@ def add_piece_wise_const_cython(double[:] x1, double[:] y1,
         else: # both arrays reached the end simultaneously
             # only the last x-value missing
             x_new[index+1] = x1[N1-1]
-        # the last value is again the end of the interval
-        # x_new[index+1] = x1[-1]
-        # only use the data that was actually filled
-        x1 = x_new[:index+2]
-        y1 = y_new[:index+1]
     # end nogil
-    return np.array(x_new[:index+2]), np.array(y_new[:index+1])
+    # return np.asarray(x_new[:index+2]), np.asarray(y_new[:index+1])
+    return np.asarray(x_new[:index+2]), np.asarray(y_new[:index+1])
 
 
 ############################################################
@@ -169,9 +165,9 @@ def add_piece_wise_lin_cython(double[:] x1, double[:] y11, double[:] y12,
         y2_new[index] = y12[N1-2]+y22[N2-2]
         # only use the data that was actually filled
     # end nogil
-    return (np.array(x_new[:index+2]),
-            np.array(y1_new[:index+1]), 
-            np.array(y2_new[:index+1]))
+    return (np.asarray(x_new[:index+2]),
+            np.asarray(y1_new[:index+1]), 
+            np.asarray(y2_new[:index+1]))
 
 
 ############################################################
@@ -230,6 +226,6 @@ def add_discrete_function_cython(double[:] x1, double[:] y1, double[:] mp1,
 
     # the last value is again the end of the interval
     # only use the data that was actually filled
-    return (np.array(x_new[:index+1]), 
-            np.array(y_new[:index+1]), 
-            np.array(mp_new[:index+1]))
+    return (np.asarray(x_new[:index+1]), 
+            np.asarray(y_new[:index+1]), 
+            np.asarray(mp_new[:index+1]))
diff --git a/pyspike/generic.py b/pyspike/generic.py
index 4f278d2..41affcb 100644
--- a/pyspike/generic.py
+++ b/pyspike/generic.py
@@ -31,6 +31,69 @@ def _generic_profile_multi(spike_trains, pair_distance_func, indices=None):
     Returns:
     - The averaged multi-variate distance of all pairs
     """
+
+    def divide_and_conquer(pairs1, pairs2):
+        """ recursive calls by splitting the two lists in half.
+        """
+        L1 = len(pairs1)
+        if L1 > 1:
+            dist_prof1 = divide_and_conquer(pairs1[:L1/2], pairs1[L1/2:])
+        else:
+            dist_prof1 = pair_distance_func(spike_trains[pairs1[0][0]],
+                                            spike_trains[pairs1[0][1]])
+        L2 = len(pairs2)
+        if L2 > 1:
+            dist_prof2 = divide_and_conquer(pairs2[:L2/2], pairs2[L2/2:])
+        else:
+            dist_prof2 = pair_distance_func(spike_trains[pairs2[0][0]],
+                                            spike_trains[pairs2[0][1]])
+        dist_prof1.add(dist_prof2)
+        return dist_prof1
+
+    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:]]
+
+    L = len(pairs)
+    if L > 1:
+        # recursive iteration through the list of pairs to get average profile
+        avrg_dist = divide_and_conquer(pairs[:len(pairs)/2],
+                                       pairs[len(pairs)/2:])
+    else:
+        avrg_dist = pair_distance_func(spike_trains[pairs[0][0]],
+                                       spike_trains[pairs[0][1]])
+
+    return avrg_dist, L
+
+
+############################################################
+# _generic_distance_multi
+############################################################
+def _generic_distance_multi(spike_trains, pair_distance_func,
+                            indices=None, interval=None):
+    """ Internal implementation detail, don't call this function directly,
+    use isi_distance_multi or spike_distance_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>} D_{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)
@@ -40,13 +103,13 @@ def _generic_profile_multi(spike_trains, pair_distance_func, indices=None):
     # 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)
+
+    avrg_dist = 0.0
+    for (i, j) in pairs:
+        avrg_dist += pair_distance_func(spike_trains[i], spike_trains[j],
+                                        interval)
+
+    return avrg_dist/len(pairs)
 
 
 ############################################################
diff --git a/pyspike/isi_distance.py b/pyspike/isi_distance.py
index aeab0df..21df561 100644
--- a/pyspike/isi_distance.py
+++ b/pyspike/isi_distance.py
@@ -3,7 +3,8 @@
 # Distributed under the BSD License
 
 from pyspike import PieceWiseConstFunc
-from pyspike.generic import _generic_profile_multi, _generic_distance_matrix
+from pyspike.generic import _generic_profile_multi, _generic_distance_multi, \
+    _generic_distance_matrix
 
 
 ############################################################
@@ -112,7 +113,8 @@ def isi_distance_multi(spike_trains, indices=None, interval=None):
     :returns: The time-averaged multivariate ISI distance :math:`D_I`
     :rtype: double
     """
-    return isi_profile_multi(spike_trains, indices).avrg(interval)
+    return _generic_distance_multi(spike_trains, isi_distance, indices,
+                                   interval)
 
 
 ############################################################
diff --git a/pyspike/spike_distance.py b/pyspike/spike_distance.py
index cc620d4..75b3b0e 100644
--- a/pyspike/spike_distance.py
+++ b/pyspike/spike_distance.py
@@ -3,7 +3,8 @@
 # Distributed under the BSD License
 
 from pyspike import PieceWiseLinFunc
-from pyspike.generic import _generic_profile_multi, _generic_distance_matrix
+from pyspike.generic import _generic_profile_multi, _generic_distance_multi, \
+    _generic_distance_matrix
 
 
 ############################################################
@@ -117,7 +118,8 @@ def spike_distance_multi(spike_trains, indices=None, interval=None):
     :returns: The averaged multi-variate spike distance :math:`D_S`.
     :rtype: double
     """
-    return spike_profile_multi(spike_trains, indices).avrg(interval)
+    return _generic_distance_multi(spike_trains, spike_distance, indices,
+                                   interval)
 
 
 ############################################################