+isi distance matrix with examples

4 files changed, 52 insertions, 4 deletions

diff --git a/examples/isi_matrix.py b/examples/isi_matrix.py
new file mode 100644
index 0000000..a149cd6
--- /dev/null
+++ b/examples/isi_matrix.py
@@ -0,0 +1,21 @@
+from __future__ import print_function
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+import pyspike as spk
+
+# first load the data
+spike_trains = []
+spike_file = open("SPIKY_testdata.txt", 'r')
+for line in spike_file:
+    spike_trains.append(spk.add_auxiliary_spikes(
+        spk.spike_train_from_string(line), 4000))
+
+print(len(spike_trains))
+
+m = spk.isi_distance_matrix(spike_trains)
+
+plt.imshow(m, interpolation='none')
+plt.show()
+
diff --git a/pyspike/__init__.py b/pyspike/__init__.py
index 2143bdc..21005e9 100644
--- a/pyspike/__init__.py
+++ b/pyspike/__init__.py
@@ -2,5 +2,5 @@ __all__ = ["function", "distances", "spikes"]
 
 from function import PieceWiseConstFunc, PieceWiseLinFunc
 from distances import add_auxiliary_spikes, isi_distance, spike_distance, \
-    isi_distance_multi, spike_distance_multi
+    isi_distance_multi, spike_distance_multi, isi_distance_matrix
 from spikes import spike_train_from_string, merge_spike_trains
diff --git a/pyspike/cython_distance.pyx b/pyspike/cython_distance.pyx
index 23ffc37..2be8525 100644
--- a/pyspike/cython_distance.pyx
+++ b/pyspike/cython_distance.pyx
@@ -29,6 +29,7 @@ import numpy as np
 cimport numpy as np
 
 from libc.math cimport fabs
+from libc.math cimport fmax
 
 DTYPE = np.float
 ctypedef np.float_t DTYPE_t
@@ -56,7 +57,7 @@ def isi_distance_cython(double[:] s1,
     isi_values = np.empty(N1+N2-1)
 
     with nogil: # release the interpreter to allow multithreading
-        isi_values[0] = (nu1-nu2)/max(nu1,nu2)
+        isi_values[0] = (nu1-nu2)/fmax(nu1,nu2)
         index1 = 0
         index2 = 0
         index = 1
@@ -84,7 +85,7 @@ def isi_distance_cython(double[:] s1,
                 nu1 = s1[index1+1]-s1[index1]
                 nu2 = s2[index2+1]-s2[index2]            
             # compute the corresponding isi-distance
-            isi_values[index] = (nu1 - nu2) / max(nu1, nu2)
+            isi_values[index] = (nu1 - nu2) / fmax(nu1, nu2)
             index += 1
         # the last event is the interval end
         spike_events[index] = s1[N1]
diff --git a/pyspike/distances.py b/pyspike/distances.py
index 35650f7..f78c0d4 100644
--- a/pyspike/distances.py
+++ b/pyspike/distances.py
@@ -206,10 +206,36 @@ def spike_distance_multi(spike_trains, indices=None):
     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, 
+    - 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)
 
+
+def isi_distance_matrix(spike_trains, indices=None):
+    """ Computes the average isi-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 isi-distance
+    """
+    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:]]
+    
+    distance_matrix = np.zeros((len(indices), len(indices)))
+    for i,j in pairs:
+        d = isi_distance(spike_trains[i], spike_trains[j]).abs_avrg()
+        distance_matrix[i,j] = d
+        distance_matrix[j,i] = d
+    return distance_matrix