add python backend for directionality

2 files changed, 105 insertions, 17 deletions

diff --git a/pyspike/directionality/cython/directionality_python_backend.py b/pyspike/directionality/cython/directionality_python_backend.py
new file mode 100644
index 0000000..e14238f
--- /dev/null
+++ b/pyspike/directionality/cython/directionality_python_backend.py
@@ -0,0 +1,89 @@
+""" directionality_python_backend.py
+
+Collection of python functions that can be used instead of the cython
+implementation.
+
+Copyright 2015, Mario Mulansky <mario.mulansky@gmx.net>
+
+Distributed under the BSD License
+
+"""
+
+import numpy as np
+
+
+############################################################
+# spike_train_order_python
+############################################################
+def spike_train_order_python(spikes1, spikes2, t_start, t_end, max_tau):
+
+    def get_tau(spikes1, spikes2, i, j, max_tau):
+        m = t_end - t_start   # use interval as initial tau
+        if i < len(spikes1)-1 and i > -1:
+            m = min(m, spikes1[i+1]-spikes1[i])
+        if j < len(spikes2)-1 and j > -1:
+            m = min(m, spikes2[j+1]-spikes2[j])
+        if i > 0:
+            m = min(m, spikes1[i]-spikes1[i-1])
+        if j > 0:
+            m = min(m, spikes2[j]-spikes2[j-1])
+        m *= 0.5
+        if max_tau > 0.0:
+            m = min(m, max_tau)
+        return m
+
+    N1 = len(spikes1)
+    N2 = len(spikes2)
+    i = -1
+    j = -1
+    n = 0
+    st = np.zeros(N1 + N2 + 2)  # spike times
+    a = np.zeros(N1 + N2 + 2)   # coincidences
+    mp = np.ones(N1 + N2 + 2)   # multiplicity
+    while i + j < N1 + N2 - 2:
+        if (i < N1-1) and (j == N2-1 or spikes1[i+1] < spikes2[j+1]):
+            i += 1
+            n += 1
+            tau = get_tau(spikes1, spikes2, i, j, max_tau)
+            st[n] = spikes1[i]
+            if j > -1 and spikes1[i]-spikes2[j] < tau:
+                # coincidence between the current spike and the previous spike
+                # both get marked with 1
+                a[n] = -1
+                a[n-1] = -1
+        elif (j < N2-1) and (i == N1-1 or spikes1[i+1] > spikes2[j+1]):
+            j += 1
+            n += 1
+            tau = get_tau(spikes1, spikes2, i, j, max_tau)
+            st[n] = spikes2[j]
+            if i > -1 and spikes2[j]-spikes1[i] < tau:
+                # coincidence between the current spike and the previous spike
+                # both get marked with 1
+                a[n] = 1
+                a[n-1] = 1
+        else:   # spikes1[i+1] = spikes2[j+1]
+            # advance in both spike trains
+            j += 1
+            i += 1
+            n += 1
+            # add only one event with zero asymmetry value and multiplicity 2
+            st[n] = spikes1[i]
+            a[n] = 0
+            mp[n] = 2
+
+    st = st[:n+2]
+    a = a[:n+2]
+    mp = mp[:n+2]
+
+    st[0] = t_start
+    st[len(st)-1] = t_end
+    if N1 + N2 > 0:
+        a[0] = a[1]
+        a[len(a)-1] = a[len(a)-2]
+        mp[0] = mp[1]
+        mp[len(mp)-1] = mp[len(mp)-2]
+    else:
+        a[0] = 1
+        a[1] = 1
+
+    return st, a, mp
diff --git a/pyspike/directionality/spike_train_order.py b/pyspike/directionality/spike_train_order.py
index f8c8615..892ffd0 100644
--- a/pyspike/directionality/spike_train_order.py
+++ b/pyspike/directionality/spike_train_order.py
@@ -5,7 +5,7 @@
 import numpy as np
 from math import exp
 from functools import partial
-# import pyspike
+import pyspike
 from pyspike import DiscreteFunc
 from pyspike.generic import _generic_profile_multi
 
@@ -39,14 +39,14 @@ def spike_train_order_profile(spike_train1, spike_train2, max_tau=None):
             spike_train_order_profile_cython as \
             spike_train_order_profile_impl
     except ImportError:
-        raise NotImplementedError()
-#         if not(pyspike.disable_backend_warning):
-#             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_profile_impl
+        # raise NotImplementedError()
+        if not(pyspike.disable_backend_warning):
+            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.directionality_python_backend import \
+            spike_train_order_python as spike_train_order_profile_impl
 
     if max_tau is None:
         max_tau = 0.0
@@ -81,15 +81,14 @@ def spike_train_order(spike_train1, spike_train2, normalize=True,
                                            spike_train1.t_start,
                                            spike_train1.t_end,
                                            max_tau)
-            if normalize:
-                return 1.0*c/mp
-            else:
-                return c
         except ImportError:
             # Cython backend not available: fall back to profile averaging
-            raise NotImplementedError()
-            # return spike_sync_profile(spike_train1, spike_train2,
-            #                           max_tau).integral(interval)
+            c, mp = spike_train_order_profile(spike_train1, spike_train2,
+                                              max_tau).integral(interval)
+        if normalize:
+            return 1.0*c/mp
+        else:
+            return c
     else:
         # some specific interval is provided: not yet implemented
         raise NotImplementedError()
@@ -250,7 +249,7 @@ def optimal_spike_train_order(spike_trains,  indices=None, interval=None,
     Returns the optimal permutation p and A value.
     """
     D = spike_train_order_matrix(spike_trains, indices, interval, max_tau)
-    return optimal_asymmetry_order_from_matrix(D, full_output)
+    return optimal_spike_train_order_from_matrix(D, full_output)
 
 
 ############################################################