From 66bf08417c651d1c0d96cb980571efb2043d7f1a Mon Sep 17 00:00:00 2001
From: Mario Mulansky <mario.mulansky@gmx.net>
Date: Wed, 26 Aug 2015 18:26:58 +0200
Subject: completed python directionality backend

---
 pyspike/cython/directionality_python_backend.py | 57 ++++++++++++++++++++++++-
 pyspike/spike_directionality.py                 | 29 +++++++------
 2 files changed, 72 insertions(+), 14 deletions(-)

diff --git a/pyspike/cython/directionality_python_backend.py b/pyspike/cython/directionality_python_backend.py
index e14238f..c1d820b 100644
--- a/pyspike/cython/directionality_python_backend.py
+++ b/pyspike/cython/directionality_python_backend.py
@@ -15,7 +15,62 @@ import numpy as np
 ############################################################
 # spike_train_order_python
 ############################################################
-def spike_train_order_python(spikes1, spikes2, t_start, t_end, max_tau):
+def spike_directionality_profile_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
+    d1 = np.zeros(N1)   # directionality values
+    d2 = np.zeros(N2)   # directionality values
+    while i + j < N1 + N2 - 2:
+        if (i < N1-1) and (j == N2-1 or spikes1[i+1] < spikes2[j+1]):
+            i += 1
+            tau = get_tau(spikes1, spikes2, i, j, max_tau)
+            if j > -1 and spikes1[i]-spikes2[j] < tau:
+                # coincidence between the current spike and the previous spike
+                # spike in first spike train occurs after second
+                d1[i] = -1
+                d2[j] = +1
+        elif (j < N2-1) and (i == N1-1 or spikes1[i+1] > spikes2[j+1]):
+            j += 1
+            tau = get_tau(spikes1, spikes2, i, j, max_tau)
+            if i > -1 and spikes2[j]-spikes1[i] < tau:
+                # coincidence between the current spike and the previous spike
+                # spike in second spike train occurs after first
+                d1[i] = +1
+                d2[j] = -1
+        else:   # spikes1[i+1] = spikes2[j+1]
+            # advance in both spike trains
+            j += 1
+            i += 1
+            d1[i] = 0
+            d2[j] = 0
+
+    return d1, d2
+
+
+############################################################
+# spike_train_order_python
+############################################################
+def spike_train_order_profile_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
diff --git a/pyspike/spike_directionality.py b/pyspike/spike_directionality.py
index f608ecc..cda7fe3 100644
--- a/pyspike/spike_directionality.py
+++ b/pyspike/spike_directionality.py
@@ -33,7 +33,11 @@ def spike_directionality(spike_train1, spike_train2, normalize=True,
                                           max_tau)
             c = len(spike_train1.spikes)
         except ImportError:
-            raise NotImplementedError()
+            d1, x = spike_directionality_profiles([spike_train1, spike_train2],
+                                                  interval=interval,
+                                                  max_tau=max_tau)
+            d = np.sum(d1)
+            c = len(spike_train1.spikes)
         if normalize:
             return 1.0*d/c
         else:
@@ -94,24 +98,23 @@ def spike_directionality_profiles(spike_trains, indices=None,
         from cython.cython_directionality import \
             spike_directionality_profiles_cython as 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
+        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_directionality_profile_python as profile_impl
 
     if max_tau is None:
         max_tau = 0.0
 
     for i, j in pairs:
-        a1, a2 = profile_impl(spike_trains[i].spikes, spike_trains[j].spikes,
+        d1, d2 = profile_impl(spike_trains[i].spikes, spike_trains[j].spikes,
                               spike_trains[i].t_start, spike_trains[i].t_end,
                               max_tau)
-        asymmetry_list[i] += a1
-        asymmetry_list[j] += a2
+        asymmetry_list[i] += d1
+        asymmetry_list[j] += d2
     for a in asymmetry_list:
         a /= len(spike_trains)-1
     return asymmetry_list
@@ -151,7 +154,7 @@ 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
+            spike_train_order_profile_python as spike_train_order_profile_impl
 
     if max_tau is None:
         max_tau = 0.0
-- 
cgit v1.2.3