From 691bf73c06322a2c47c37a5c48d085b789c8e8bf Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Mon, 17 Aug 2015 18:02:24 +0200 Subject: add python backend for directionality --- .../cython/directionality_python_backend.py | 89 ++++++++++++++++++++++ pyspike/directionality/spike_train_order.py | 33 ++++---- 2 files changed, 105 insertions(+), 17 deletions(-) create mode 100644 pyspike/directionality/cython/directionality_python_backend.py 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 + +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) ############################################################ -- cgit v1.2.3