+cython version for isi and performance tests

2 files changed, 89 insertions, 3 deletions

diff --git a/pyspike/cython_distance.pyx b/pyspike/cython_distance.pyx
new file mode 100644
index 0000000..330eea4
--- /dev/null
+++ b/pyspike/cython_distance.pyx
@@ -0,0 +1,75 @@
+#cython: boundscheck=False
+#cython: wraparound=False
+#cython: cdivision=True
+
+"""
+Doc
+
+To test whether things can be optimized: remove all yellow stuff
+in the html output::
+
+  cython -a cython_distance.pyx
+
+which gives::
+
+  cython_distance.html
+
+
+"""
+
+import numpy as np
+cimport numpy as np
+
+DTYPE = np.float
+ctypedef np.float_t DTYPE_t
+
+def isi_distance_cython(np.ndarray[DTYPE_t, ndim=1] s1, np.ndarray[DTYPE_t, ndim=1] s2):
+
+    cdef np.ndarray[DTYPE_t, ndim=1] spike_events
+    # the values have one entry less - the number of intervals between events
+    cdef np.ndarray[DTYPE_t, ndim=1] isi_values
+    cdef int index1, index2, index
+    cdef int N1, N2
+    cdef double nu1, nu2
+    N1 = len(s1)-1
+    N2 = len(s2)-1
+
+    nu1 = s1[1]-s1[0]
+    nu2 = s2[1]-s2[0]
+    spike_events = np.empty(N1+N2)
+    spike_events[0] = s1[0]
+    isi_values = np.empty(N1+N2-1)
+    isi_values[0] = (nu1-nu2)/max(nu1,nu2)
+    index1 = 0
+    index2 = 0
+    index = 1
+    while True:
+        # check which spike is next - from s1 or s2
+        if s1[index1+1] < s2[index2+1]:
+            index1 += 1
+            # break condition relies on existence of spikes at T_end
+            if index1 >= N1:
+                break
+            spike_events[index] = s1[index1]
+            nu1 = s1[index1+1]-s1[index1]
+        elif s1[index1+1] > s2[index2+1]:
+            index2 += 1
+            if index2 >= N2:
+                break
+            spike_events[index] = s2[index2]
+            nu2 = s2[index2+1]-s2[index2]
+        else: # s1[index1+1] == s2[index2+1]
+            index1 += 1
+            index2 += 1
+            if (index1 >= N1) or (index2 >= N2):
+                break
+            spike_events[index] = s1[index1]
+            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)
+        index += 1
+    # the last event is the interval end
+    spike_events[index] = s1[N1]
+
+    return spike_events[:index+1], isi_values[:index]
diff --git a/pyspike/distances.py b/pyspike/distances.py
index 52c6640..76dcd83 100644
--- a/pyspike/distances.py
+++ b/pyspike/distances.py
@@ -53,10 +53,21 @@ def isi_distance(spikes1, spikes2):
     assert spikes1[-1]==spikes2[-1], \
         "Given spike trains seems not to have auxiliary spikes!"
 
-    # shorter names
-    s1 = spikes1
-    s2 = spikes2
+    # compile and load cython implementation
+    import pyximport
+    pyximport.install(setup_args={'include_dirs': [np.get_include()]})
+    from cython_distance import isi_distance_cython
 
+    times, values = isi_distance_cython(spikes1, spikes2)
+    return PieceWiseConstFunc(times, values)
+
+
+############################################################
+# isi_distance_python
+############################################################
+def isi_distance_python(s1, s2):
+    """ Plain Python implementation of the isi distance.
+    """
     # compute the interspike interval
     nu1 = s1[1:]-s1[:-1]
     nu2 = s2[1:]-s2[:-1]