added SPIKE-Sync based filtering

new function filter_by_spike_sync removes spikes that have a multi-variate
Spike Sync value below some threshold

not yet fully tested, python backend missing.

2 files changed, 70 insertions, 1 deletions

diff --git a/pyspike/cython/cython_profiles.pyx b/pyspike/cython/cython_profiles.pyx
index f9893eb..3b33c5a 100644
--- a/pyspike/cython/cython_profiles.pyx
+++ b/pyspike/cython/cython_profiles.pyx
@@ -426,3 +426,34 @@ def coincidence_profile_cython(double[:] spikes1, double[:] spikes2,
         c[1] = 1
 
     return st, c, mp
+
+
+############################################################
+# coincidence_single_profile_cython
+############################################################
+def coincidence_single_profile_cython(double[:] spikes1, double[:] spikes2,
+                                      double t_start, double t_end, double max_tau):
+
+    cdef int N1 = len(spikes1)
+    cdef int N2 = len(spikes2)
+    cdef int j = -1
+    cdef double[:] c = np.zeros(N1)   # coincidences
+    cdef double interval = t_end - t_start
+    cdef double tau
+    for i in xrange(N1):
+        while j < N2-1 and spikes2[j+1] < spikes1[i]:
+            j += 1
+        tau = get_tau(spikes1, spikes2, i, j, interval, max_tau)
+        print i, j, spikes1[i], spikes2[j], tau
+        if j > -1 and spikes1[i]-spikes2[j] < tau:
+            # current spike in st1 is coincident
+            c[i] = 1
+        if j < N2-1:
+            j += 1
+            tau = get_tau(spikes1, spikes2, i, j, interval, max_tau)
+            print i, j, spikes1[i], spikes2[j], tau
+            if spikes2[j]-spikes1[i] < tau:
+                # current spike in st1 is coincident
+                c[i] = 1
+
+    return c
diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py
index 10ebdc7..f99018a 100644
--- a/pyspike/spike_sync.py
+++ b/pyspike/spike_sync.py
@@ -6,7 +6,7 @@
 import numpy as np
 from functools import partial
 import pyspike
-from pyspike import DiscreteFunc
+from pyspike import DiscreteFunc, SpikeTrain
 from pyspike.generic import _generic_profile_multi, _generic_distance_matrix
 
 
@@ -212,3 +212,41 @@ def spike_sync_matrix(spike_trains, indices=None, interval=None, max_tau=None):
     dist_func = partial(spike_sync, max_tau=max_tau)
     return _generic_distance_matrix(spike_trains, dist_func,
                                     indices, interval)
+
+
+############################################################
+# filter_by_spike_sync
+############################################################
+def filter_by_spike_sync(spike_trains, threshold, indices=None, max_tau=None):
+    """ Removes the spikes with a multi-variate spike_sync value below
+    threshold.
+    """
+    N = len(spike_trains)
+    filtered_spike_trains = []
+
+    # cython implementation
+    try:
+        from cython.cython_profiles import coincidence_single_profile_cython \
+            as coincidence_impl
+    except ImportError:
+        if not(pyspike.disable_backend_warning):
+            print("Warning: coincidence_single_profile_cytho 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_single_profile_python \
+            as coincidence_impl
+
+    if max_tau is None:
+        max_tau = 0.0
+
+    for i, st in enumerate(spike_trains):
+        coincidences = np.zeros_like(st)
+        for j in range(N).remove(i):
+            coincidences += coincidence_impl(st.spikes, spike_trains[j].spikes,
+                                             st.t_start, st.t_end, max_tau)
+        filtered_spikes = st[coincidences > threshold*(N-1)]
+        filtered_spike_trains.append(SpikeTrain(filtered_spikes,
+                                                [st.t_start, st.t_end]))
+    return filtered_spike_trains