Merge branch 'new_directionality' into develop

12 files changed, 1353 insertions, 27 deletions

diff --git a/pyspike/__init__.py b/pyspike/__init__.py
index 08253fb..3cf416f 100644
--- a/pyspike/__init__.py
+++ b/pyspike/__init__.py
@@ -7,8 +7,8 @@ Distributed under the BSD License
 from __future__ import absolute_import
 
 __all__ = ["isi_distance", "spike_distance", "spike_sync", "psth",
-           "spikes", "SpikeTrain", "PieceWiseConstFunc", "PieceWiseLinFunc",
-           "DiscreteFunc", "directionality"]
+           "spikes", "spike_directionality", "SpikeTrain",
+           "PieceWiseConstFunc", "PieceWiseLinFunc", "DiscreteFunc"]
 
 from .PieceWiseConstFunc import PieceWiseConstFunc
 from .PieceWiseLinFunc import PieceWiseLinFunc
@@ -20,13 +20,26 @@ from .isi_distance import isi_profile, isi_distance, isi_profile_multi,\
 from .spike_distance import spike_profile, spike_distance, spike_profile_multi,\
     spike_distance_multi, spike_distance_matrix
 from .spike_sync import spike_sync_profile, spike_sync,\
-    spike_sync_profile_multi, spike_sync_multi, spike_sync_matrix
+    spike_sync_profile_multi, spike_sync_multi, spike_sync_matrix,\
+    filter_by_spike_sync
 from .psth import psth
 
 from .spikes import load_spike_trains_from_txt, save_spike_trains_to_txt, \
     spike_train_from_string, import_spike_trains_from_time_series, \
     merge_spike_trains, generate_poisson_spikes
 
+from .spike_directionality import spike_directionality, \
+    spike_directionality_profiles, spike_directionality_matrix, \
+    spike_train_order_profile, spike_train_order, \
+    spike_train_order_profile_multi, optimal_spike_train_order_from_matrix, \
+    optimal_spike_train_order, permutate_matrix
+
+from .spike_directionality import spike_directionality, \
+    spike_directionality_profiles, spike_directionality_matrix, \
+    spike_train_order_profile, spike_train_order, \
+    spike_train_order_profile_multi, optimal_spike_train_order_from_matrix, \
+    optimal_spike_train_order, permutate_matrix
+
 # define the __version__ following
 # http://stackoverflow.com/questions/17583443
 from pkg_resources import get_distribution, DistributionNotFound
diff --git a/pyspike/cython/cython_directionality.pyx b/pyspike/cython/cython_directionality.pyx
new file mode 100644
index 0000000..ac37690
--- /dev/null
+++ b/pyspike/cython/cython_directionality.pyx
@@ -0,0 +1,262 @@
+#cython: boundscheck=False
+#cython: wraparound=False
+#cython: cdivision=True
+
+"""
+cython_directionality.pyx
+
+cython implementation of the spike delay asymmetry measures
+
+Copyright 2015, Mario Mulansky <mario.mulansky@gmx.net>
+
+Distributed under the BSD License
+
+"""
+
+"""
+To test whether things can be optimized: remove all yellow stuff
+in the html output::
+
+  cython -a cython_directionality.pyx
+
+which gives::
+
+  cython_directionality.html
+
+"""
+
+import numpy as np
+cimport numpy as np
+
+from libc.math cimport fabs
+from libc.math cimport fmax
+from libc.math cimport fmin
+
+# from pyspike.cython.cython_distances cimport get_tau
+
+DTYPE = np.float
+ctypedef np.float_t DTYPE_t
+
+
+############################################################
+# get_tau
+############################################################
+cdef inline double get_tau(double[:] spikes1, double[:] spikes2,
+                           int i, int j, double interval, double max_tau):
+    cdef double m = interval   # use interval length as initial tau
+    cdef int N1 = spikes1.shape[0]-1  # len(spikes1)-1
+    cdef int N2 = spikes2.shape[0]-1  # len(spikes2)-1
+    if i < N1 and i > -1:
+        m = fmin(m, spikes1[i+1]-spikes1[i])
+    if j < N2 and j > -1:
+        m = fmin(m, spikes2[j+1]-spikes2[j])
+    if i > 0:
+        m = fmin(m, spikes1[i]-spikes1[i-1])
+    if j > 0:
+        m = fmin(m, spikes2[j]-spikes2[j-1])
+    m *= 0.5
+    if max_tau > 0.0:
+        m = fmin(m, max_tau)
+    return m
+
+
+############################################################
+# spike_train_order_profile_cython
+############################################################
+def spike_train_order_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 i = -1
+    cdef int j = -1
+    cdef int n = 0
+    cdef double[:] st = np.zeros(N1 + N2 + 2)  # spike times
+    cdef double[:] a = np.zeros(N1 + N2 + 2)   # asymmetry values
+    cdef double[:] mp = np.ones(N1 + N2 + 2)   # multiplicity
+    cdef double interval = t_end - t_start
+    cdef double tau
+    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, interval, max_tau)
+            st[n] = spikes1[i]
+            if j > -1 and spikes1[i]-spikes2[j] < tau:
+                # coincidence between the current spike and the previous spike
+                # spike from spike train 1 after spike train 2
+                # 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, interval, max_tau)
+            st[n] = spikes2[j]
+            if i > -1 and spikes2[j]-spikes1[i] < tau:
+                # coincidence between the current spike and the previous spike
+                # spike from spike train 1 before spike train 2
+                # 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
+
+
+############################################################
+# spike_train_order_cython
+############################################################
+def spike_train_order_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 i = -1
+    cdef int j = -1
+    cdef int d = 0
+    cdef int mp = 0
+    cdef double interval = t_end - t_start
+    cdef double tau
+    while i + j < N1 + N2 - 2:
+        if (i < N1-1) and (j == N2-1 or spikes1[i+1] < spikes2[j+1]):
+            i += 1
+            mp += 1
+            tau = get_tau(spikes1, spikes2, i, j, interval, max_tau)
+            if j > -1 and spikes1[i]-spikes2[j] < tau:
+                # coincidence between the current spike and the previous spike
+                # spike in spike train 2 appeared before spike in spike train 1
+                # mark with -1
+                d -= 2
+        elif (j < N2-1) and (i == N1-1 or spikes1[i+1] > spikes2[j+1]):
+            j += 1
+            mp += 1
+            tau = get_tau(spikes1, spikes2, i, j, interval, max_tau)
+            if i > -1 and spikes2[j]-spikes1[i] < tau:
+                # coincidence between the current spike and the previous spike
+                # spike in spike train 1 appeared before spike in spike train 2
+                # mark with +1
+                d += 2
+        else:   # spikes1[i+1] = spikes2[j+1]
+            # advance in both spike trains
+            j += 1
+            i += 1
+            # add only one event with multiplicity 2, but no asymmetry counting
+            mp += 2
+
+    if d == 0 and mp == 0:
+        # empty spike trains -> spike sync = 1 by definition
+        d = 1
+        mp = 1
+
+    return d, mp
+
+
+############################################################
+# spike_directionality_profiles_cython
+############################################################
+def spike_directionality_profiles_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 i = -1
+    cdef int j = -1
+    cdef double[:] d1 = np.zeros(N1)  # directionality values
+    cdef double[:] d2 = np.zeros(N2)  # directionality values
+    cdef double interval = t_end - t_start
+    cdef double tau
+    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, interval, max_tau)
+            if j > -1 and spikes1[i]-spikes2[j] < tau:
+                # coincidence between the current spike and the previous spike
+                # spike from spike train 1 after spike train 2
+                # leading spike gets +1, following spike -1
+                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, interval, max_tau)
+            if i > -1 and spikes2[j]-spikes1[i] < tau:
+                # coincidence between the current spike and the previous spike
+                # spike from spike train 1 before spike train 2
+                # leading spike gets +1, following spike -1
+                d1[i] = +1
+                d2[j] = -1
+        else:   # spikes1[i+1] = spikes2[j+1]
+            # advance in both spike trains
+            j += 1
+            i += 1
+            # equal spike times: zero asymmetry value
+            d1[i] = 0
+            d2[j] = 0
+
+    return d1, d2
+
+
+############################################################
+# spike_directionality_cython
+############################################################
+def spike_directionality_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 i = -1
+    cdef int j = -1
+    cdef int d = 0  # directionality value
+    cdef double interval = t_end - t_start
+    cdef double tau
+    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, interval, max_tau)
+            if j > -1 and spikes1[i]-spikes2[j] < tau:
+                # coincidence between the current spike and the previous spike
+                # spike from spike train 1 after spike train 2
+                # leading spike gets +1, following spike -1
+                d -= 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, interval, max_tau)
+            if i > -1 and spikes2[j]-spikes1[i] < tau:
+                # coincidence between the current spike and the previous spike
+                # spike from spike train 1 before spike train 2
+                # leading spike gets +1, following spike -1
+                d += 1
+        else:   # spikes1[i+1] = spikes2[j+1]
+            # advance in both spike trains
+            j += 1
+            i += 1
+
+    return d
diff --git a/pyspike/cython/cython_distances.pyx b/pyspike/cython/cython_distances.pyx
index ac5f226..d4070ae 100644
--- a/pyspike/cython/cython_distances.pyx
+++ b/pyspike/cython/cython_distances.pyx
@@ -178,6 +178,8 @@ cdef inline double isi_avrg_cython(double isi1, double isi2) nogil:
     return 0.5*(isi1+isi2)*(isi1+isi2)
     # alternative definition to obtain <S> ~ 0.5 for Poisson spikes
     # return 0.5*(isi1*isi1+isi2*isi2)
+    # another alternative definition without second normalization
+    # return 0.5*(isi1+isi2)
 
 
 ############################################################
@@ -248,6 +250,8 @@ def spike_distance_cython(double[:] t1, double[:] t2,
             index2 = 0
 
         y_start = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2)
+        # alternative definition without second normalization
+        # y_start = (s1 + s2) / isi_avrg_cython(isi1, isi2)
         index = 1
 
         while index1+index2 < N1+N2-2:
@@ -267,6 +271,8 @@ def spike_distance_cython(double[:] t1, double[:] t2,
                 t_curr =  t_p1
                 s2 = (dt_p2*(t_f2-t_p1) + dt_f2*(t_p1-t_p2)) / isi2
                 y_end = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2)
+                # alternative definition without second normalization
+                # y_end = (s1 + s2) / isi_avrg_cython(isi1, isi2)
 
                 spike_value += 0.5*(y_start + y_end) * (t_curr - t_last)
 
@@ -286,6 +292,8 @@ def spike_distance_cython(double[:] t1, double[:] t2,
                     s1 = dt_p1
                 # s2 is the same as above, thus we can compute y2 immediately
                 y_start = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2)
+                # alternative definition without second normalization
+                # y_start = (s1 + s2) / isi_avrg_cython(isi1, isi2)
             elif (index2 < N2-1) and (t_f1 > t_f2 or index1 == N1-1):
                 index2 += 1
                 # first calculate the previous interval end value
@@ -301,6 +309,8 @@ def spike_distance_cython(double[:] t1, double[:] t2,
                 t_curr = t_p2
                 s1 = (dt_p1*(t_f1-t_p2) + dt_f1*(t_p2-t_p1)) / isi1
                 y_end = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2)
+                # alternative definition without second normalization
+                # y_end = (s1 + s2) / isi_avrg_cython(isi1, isi2)
 
                 spike_value += 0.5*(y_start + y_end) * (t_curr - t_last)
 
@@ -320,6 +330,9 @@ def spike_distance_cython(double[:] t1, double[:] t2,
                     s2 = dt_p2
                 # s1 is the same as above, thus we can compute y2 immediately
                 y_start = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2)
+                # alternative definition without second normalization
+                # y_start = (s1 + s2) / isi_avrg_cython(isi1, isi2)
+
             else: # t_f1 == t_f2 - generate only one event
                 index1 += 1
                 index2 += 1
@@ -358,6 +371,193 @@ def spike_distance_cython(double[:] t1, double[:] t2,
         s1 = dt_f1 # *(t_end-t1[N1-1])/isi1
         s2 = dt_f2 # *(t_end-t2[N2-1])/isi2
         y_end = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2)
+        # alternative definition without second normalization
+        # y_end = (s1 + s2) / isi_avrg_cython(isi1, isi2)
+
+        spike_value += 0.5*(y_start + y_end) * (t_end - t_last)
+    # end nogil
+
+    # use only the data added above 
+    # could be less than original length due to equal spike times
+    return spike_value / (t_end-t_start)
+
+
+############################################################
+# isi_avrg_rf_cython
+############################################################
+cdef inline double isi_avrg_rf_cython(double isi1, double isi2) nogil:
+    # rate free version
+    return (isi1+isi2)
+
+
+############################################################
+# spike_distance_rf_cython
+############################################################
+def spike_distance_rf_cython(double[:] t1, double[:] t2,
+                             double t_start, double t_end):
+
+    cdef int N1, N2, index1, index2, index
+    cdef double t_p1, t_f1, t_p2, t_f2, dt_p1, dt_p2, dt_f1, dt_f2
+    cdef double isi1, isi2, s1, s2
+    cdef double y_start, y_end, t_last, t_current, spike_value
+    
+    spike_value = 0.0
+
+    N1 = len(t1)
+    N2 = len(t2)
+
+    with nogil: # release the interpreter to allow multithreading
+        t_last = t_start
+        t_p1 = t_start
+        t_p2 = t_start
+        if t1[0] > t_start:
+            # dt_p1 = t2[0]-t_start
+            t_f1 = t1[0]
+            dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_start, t_end)
+            isi1 = fmax(t_f1-t_start, t1[1]-t1[0])
+            dt_p1 = dt_f1
+            s1 = dt_p1*(t_f1-t_start)/isi1
+            index1 = -1
+        else:
+            t_f1 = t1[1]
+            dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_start, t_end)
+            dt_p1 = 0.0
+            isi1 = t1[1]-t1[0]
+            s1 = dt_p1
+            index1 = 0
+        if t2[0] > t_start:
+            # dt_p1 = t2[0]-t_start
+            t_f2 = t2[0]
+            dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_start, t_end)
+            dt_p2 = dt_f2
+            isi2 = fmax(t_f2-t_start, t2[1]-t2[0])
+            s2 = dt_p2*(t_f2-t_start)/isi2
+            index2 = -1
+        else:
+            t_f2 = t2[1]
+            dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_start, t_end)
+            dt_p2 = 0.0
+            isi2 = t2[1]-t2[0]
+            s2 = dt_p2
+            index2 = 0
+
+        # y_start = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2)
+        # alternative definition without second normalization
+        y_start = (s1 + s2) / isi_avrg_rf_cython(isi1, isi2)
+        index = 1
+
+        while index1+index2 < N1+N2-2:
+            # print(index, index1, index2)
+            if (index1 < N1-1) and (t_f1 < t_f2 or index2 == N2-1):
+                index1 += 1
+                # first calculate the previous interval end value
+                s1 = dt_f1*(t_f1-t_p1) / isi1
+                # the previous time now was the following time before:
+                dt_p1 = dt_f1
+                t_p1 = t_f1    # t_p1 contains the current time point
+                # get the next time
+                if index1 < N1-1:
+                    t_f1 = t1[index1+1]
+                else:
+                    t_f1 = t_end
+                t_curr =  t_p1
+                s2 = (dt_p2*(t_f2-t_p1) + dt_f2*(t_p1-t_p2)) / isi2
+                # y_end = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2)
+                # alternative definition without second normalization
+                y_end = (s1 + s2) / isi_avrg_rf_cython(isi1, isi2)
+
+                spike_value += 0.5*(y_start + y_end) * (t_curr - t_last)
+
+                # now the next interval start value
+                if index1 < N1-1:
+                    dt_f1 = get_min_dist_cython(t_f1, t2, N2, index2,
+                                                t_start, t_end)
+                    isi1 = t_f1-t_p1
+                    s1 = dt_p1
+                else:
+                    dt_f1 = dt_p1
+                    isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
+                    # s1 needs adjustment due to change of isi1
+                    s1 = dt_p1*(t_end-t1[N1-1])/isi1
+                # s2 is the same as above, thus we can compute y2 immediately
+                # y_start = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2)
+                # alternative definition without second normalization
+                y_start = (s1 + s2) / isi_avrg_rf_cython(isi1, isi2)
+            elif (index2 < N2-1) and (t_f1 > t_f2 or index1 == N1-1):
+                index2 += 1
+                # first calculate the previous interval end value
+                s2 = dt_f2*(t_f2-t_p2) / isi2
+                # the previous time now was the following time before:
+                dt_p2 = dt_f2
+                t_p2 = t_f2    # t_p2 contains the current time point
+                # get the next time
+                if index2 < N2-1:
+                    t_f2 = t2[index2+1]
+                else:
+                    t_f2 = t_end
+                t_curr = t_p2
+                s1 = (dt_p1*(t_f1-t_p2) + dt_f1*(t_p2-t_p1)) / isi1
+                # y_end = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2)
+                # alternative definition without second normalization
+                y_end = (s1 + s2) / isi_avrg_rf_cython(isi1, isi2)
+
+                spike_value += 0.5*(y_start + y_end) * (t_curr - t_last)
+
+                # now the next interval start value
+                if index2 < N2-1:
+                    dt_f2 = get_min_dist_cython(t_f2, t1, N1, index1,
+                                                t_start, t_end)
+                    isi2 = t_f2-t_p2
+                    s2 = dt_p2
+                else:
+                    dt_f2 = dt_p2
+                    isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
+                    # s2 needs adjustment due to change of isi2
+                    s2 = dt_p2*(t_end-t2[N2-1])/isi2
+                # s1 is the same as above, thus we can compute y2 immediately
+                # y_start = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2)
+                # alternative definition without second normalization
+                y_start = (s1 + s2) / isi_avrg_rf_cython(isi1, isi2)
+
+            else: # t_f1 == t_f2 - generate only one event
+                index1 += 1
+                index2 += 1
+                t_p1 = t_f1
+                t_p2 = t_f2
+                dt_p1 = 0.0
+                dt_p2 = 0.0
+                t_curr = t_f1
+                y_end = 0.0
+                spike_value += 0.5*(y_start + y_end) * (t_curr - t_last)
+                y_start = 0.0
+                if index1 < N1-1:
+                    t_f1 = t1[index1+1]
+                    dt_f1 = get_min_dist_cython(t_f1, t2, N2, index2,
+                                                t_start, t_end)
+                    isi1 = t_f1 - t_p1
+                else:
+                    t_f1 = t_end
+                    dt_f1 = dt_p1
+                    isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
+                if index2 < N2-1:
+                    t_f2 = t2[index2+1]
+                    dt_f2 = get_min_dist_cython(t_f2, t1, N1, index1,
+                                                t_start, t_end)
+                    isi2 = t_f2 - t_p2
+                else:
+                    t_f2 = t_end
+                    dt_f2 = dt_p2
+                    isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
+            index += 1
+            t_last = t_curr
+        # isi1 = max(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
+        # isi2 = max(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
+        s1 = dt_f1*(t_end-t1[N1-1])/isi1
+        s2 = dt_f2*(t_end-t2[N2-1])/isi2
+        # y_end = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2)
+        # alternative definition without second normalization
+        y_end = (s1 + s2) / isi_avrg_rf_cython(isi1, isi2)
+
         spike_value += 0.5*(y_start + y_end) * (t_end - t_last)
     # end nogil
 
diff --git a/pyspike/cython/cython_profiles.pyx b/pyspike/cython/cython_profiles.pyx
index 4a42cdb..aa24db4 100644
--- a/pyspike/cython/cython_profiles.pyx
+++ b/pyspike/cython/cython_profiles.pyx
@@ -450,3 +450,36 @@ 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]:
+            # move forward until spikes2[j] is the last spike before spikes1[i]
+            # note that if spikes2[j] is after spikes1[i] we dont do anything
+            j += 1
+        tau = get_tau(spikes1, spikes2, i, j, interval, max_tau)
+        if j > -1 and fabs(spikes1[i]-spikes2[j]) < tau:
+            # current spike in st1 is coincident
+            c[i] = 1
+        if j < N2-1 and (j < 0 or spikes2[j] < spikes1[i]):
+            # in case spikes2[j] is before spikes1[i] it has to be the one 
+            # right before (see above), hence we move one forward and also 
+            # check the next spike
+            j += 1
+            tau = get_tau(spikes1, spikes2, i, j, interval, max_tau)
+            if fabs(spikes2[j]-spikes1[i]) < tau:
+                # current spike in st1 is coincident
+                c[i] = 1
+    return c
diff --git a/pyspike/cython/cython_simulated_annealing.pyx b/pyspike/cython/cython_simulated_annealing.pyx
new file mode 100644
index 0000000..be9423c
--- /dev/null
+++ b/pyspike/cython/cython_simulated_annealing.pyx
@@ -0,0 +1,82 @@
+#cython: boundscheck=False
+#cython: wraparound=False
+#cython: cdivision=True
+
+"""
+cython_simulated_annealing.pyx
+
+cython implementation of a simulated annealing algorithm to find the optimal
+spike train order
+
+Note: using cython memoryviews (e.g. double[:]) instead of ndarray objects
+improves the performance of spike_distance by a factor of 10!
+
+Copyright 2015, Mario Mulansky <mario.mulansky@gmx.net>
+
+Distributed under the BSD License
+
+"""
+
+"""
+To test whether things can be optimized: remove all yellow stuff
+in the html output::
+
+  cython -a cython_simulated_annealing.pyx
+
+which gives:
+
+  cython_simulated_annealing.html
+
+"""
+
+import numpy as np
+cimport numpy as np
+
+from libc.math cimport exp
+from libc.math cimport fmod
+from libc.stdlib cimport rand
+from libc.stdlib cimport RAND_MAX
+
+DTYPE = np.float
+ctypedef np.float_t DTYPE_t
+
+
+def sim_ann_cython(double[:, :] D, double T_start, double T_end, double alpha):
+
+    cdef long N = len(D)
+    cdef double A = np.sum(np.triu(D, 0))
+    cdef long[:] p = np.arange(N)
+    cdef double T = T_start
+    cdef long iterations
+    cdef long succ_iter
+    cdef long total_iter = 0
+    cdef double delta_A
+    cdef long ind1
+    cdef long ind2
+
+    while T > T_end:
+        iterations = 0
+        succ_iter = 0
+        # equilibrate for 100*N steps or 10*N successful steps
+        while iterations < 100*N and succ_iter < 10*N:
+            # exchange two rows and cols
+            # ind1 = np.random.randint(N-1)
+            ind1 = rand() % (N-1)
+            if ind1 < N-1:
+                ind2 = ind1+1
+            else:  # this can never happen!
+                ind2 = 0
+            delta_A = -2*D[p[ind1], p[ind2]]
+            if delta_A > 0.0 or exp(delta_A/T) > ((1.0*rand()) / RAND_MAX):
+                # swap indices
+                p[ind1], p[ind2] = p[ind2], p[ind1]
+                A += delta_A
+                succ_iter += 1
+            iterations += 1
+        total_iter += iterations
+        T *= alpha   # cool down
+        if succ_iter == 0:
+            # no successful step -> we believe we have converged
+            break
+
+    return p, A, total_iter
diff --git a/pyspike/cython/directionality_python_backend.py b/pyspike/cython/directionality_python_backend.py
new file mode 100644
index 0000000..c1d820b
--- /dev/null
+++ b/pyspike/cython/directionality_python_backend.py
@@ -0,0 +1,144 @@
+""" 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_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
+        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/cython/python_backend.py b/pyspike/cython/python_backend.py
index 6b7209a..e75f181 100644
--- a/pyspike/cython/python_backend.py
+++ b/pyspike/cython/python_backend.py
@@ -3,7 +3,7 @@
 Collection of python functions that can be used instead of the cython
 implementation.
 
-Copyright 2014, Mario Mulansky <mario.mulansky@gmx.net>
+Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net>
 
 Distributed under the BSD License
 
@@ -356,26 +356,27 @@ def cumulative_sync_python(spikes1, spikes2):
     return st, c
 
 
+def get_tau(spikes1, spikes2, i, j, max_tau, init_tau):
+    m = init_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
+
+
 ############################################################
 # coincidence_python
 ############################################################
 def coincidence_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
@@ -388,7 +389,7 @@ def coincidence_python(spikes1, spikes2, t_start, t_end, max_tau):
         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)
+            tau = get_tau(spikes1, spikes2, i, j, max_tau, t_end-t_start)
             st[n] = spikes1[i]
             if j > -1 and spikes1[i]-spikes2[j] < tau:
                 # coincidence between the current spike and the previous spike
@@ -398,7 +399,7 @@ def coincidence_python(spikes1, spikes2, t_start, t_end, max_tau):
         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)
+            tau = get_tau(spikes1, spikes2, i, j, max_tau, t_end-t_start)
             st[n] = spikes2[j]
             if i > -1 and spikes2[j]-spikes1[i] < tau:
                 # coincidence between the current spike and the previous spike
@@ -434,6 +435,36 @@ def coincidence_python(spikes1, spikes2, t_start, t_end, max_tau):
 
 
 ############################################################
+# coincidence_single_profile_cython
+############################################################
+def coincidence_single_python(spikes1, spikes2, t_start, t_end, max_tau):
+
+    N1 = len(spikes1)
+    N2 = len(spikes2)
+    j = -1
+    c = np.zeros(N1)   # coincidences
+    for i in range(N1):
+        while j < N2-1 and spikes2[j+1] < spikes1[i]:
+            # move forward until spikes2[j] is the last spike before spikes1[i]
+            # note that if spikes2[j] is after spikes1[i] we dont do anything
+            j += 1
+        tau = get_tau(spikes1, spikes2, i, j, max_tau, t_end-t_start)
+        if j > -1 and abs(spikes1[i]-spikes2[j]) < tau:
+            # current spike in st1 is coincident
+            c[i] = 1
+        if j < N2-1 and (j < 0 or spikes2[j] < spikes1[i]):
+            # in case spikes2[j] is before spikes1[i] it has to be the first or
+            # the one right before (see above), hence we move one forward and
+            # also check the next spike
+            j += 1
+            tau = get_tau(spikes1, spikes2, i, j, max_tau, t_end-t_start)
+            if abs(spikes2[j]-spikes1[i]) < tau:
+                # current spike in st1 is coincident
+                c[i] = 1
+    return c
+
+
+############################################################
 # add_piece_wise_const_python
 ############################################################
 def add_piece_wise_const_python(x1, y1, x2, y2):
diff --git a/pyspike/spike_directionality.py b/pyspike/spike_directionality.py
new file mode 100644
index 0000000..d1a525e
--- /dev/null
+++ b/pyspike/spike_directionality.py
@@ -0,0 +1,312 @@
+# Module containing functions to compute the SPIKE directionality and the
+# spike train order profile
+# Copyright 2015, Mario Mulansky <mario.mulansky@gmx.net>
+# Distributed under the BSD License
+
+from __future__ import absolute_import
+
+import numpy as np
+import pyspike
+from pyspike import DiscreteFunc
+from functools import partial
+from pyspike.generic import _generic_profile_multi
+
+
+############################################################
+# spike_directionality
+############################################################
+def spike_directionality(spike_train1, spike_train2, normalize=True,
+                         interval=None, max_tau=None):
+    """ Computes the overall spike directionality for two spike trains.
+    """
+    if interval is None:
+        # distance over the whole interval is requested: use specific function
+        # for optimal performance
+        try:
+            from .cython.cython_directionality import \
+                spike_directionality_cython as spike_directionality_impl
+            if max_tau is None:
+                max_tau = 0.0
+            d = spike_directionality_impl(spike_train1.spikes,
+                                          spike_train2.spikes,
+                                          spike_train1.t_start,
+                                          spike_train1.t_end,
+                                          max_tau)
+            c = len(spike_train1.spikes)
+        except ImportError:
+            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:
+            return d
+    else:
+        # some specific interval is provided: not yet implemented
+        raise NotImplementedError()
+
+
+############################################################
+# spike_directionality_matrix
+############################################################
+def spike_directionality_matrix(spike_trains, normalize=True, indices=None,
+                                interval=None, max_tau=None):
+    """ Computes the spike directionaity matrix for the given spike trains.
+    """
+    if indices is None:
+        indices = np.arange(len(spike_trains))
+    indices = np.array(indices)
+    # check validity of indices
+    assert (indices < len(spike_trains)).all() and (indices >= 0).all(), \
+        "Invalid index list."
+    # generate a list of possible index pairs
+    pairs = [(indices[i], j) for i in range(len(indices))
+             for j in indices[i+1:]]
+
+    distance_matrix = np.zeros((len(indices), len(indices)))
+    for i, j in pairs:
+        d = spike_directionality(spike_trains[i], spike_trains[j], normalize,
+                                 interval, max_tau=max_tau)
+        distance_matrix[i, j] = d
+        distance_matrix[j, i] = -d
+    return distance_matrix
+
+
+############################################################
+# spike_directionality_profiles
+############################################################
+def spike_directionality_profiles(spike_trains, indices=None,
+                                  interval=None, max_tau=None):
+    """ Computes the spike directionality value for each spike in each spike
+    train.
+    """
+    if indices is None:
+        indices = np.arange(len(spike_trains))
+    indices = np.array(indices)
+    # check validity of indices
+    assert (indices < len(spike_trains)).all() and (indices >= 0).all(), \
+        "Invalid index list."
+    # list of arrays for reulting asymmetry values
+    asymmetry_list = [np.zeros_like(st.spikes) for st in spike_trains]
+    # generate a list of possible index pairs
+    pairs = [(indices[i], j) for i in range(len(indices))
+             for j in indices[i+1:]]
+
+    # cython implementation
+    try:
+        from .cython.cython_directionality import \
+            spike_directionality_profiles_cython as profile_impl
+    except ImportError:
+        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:
+        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] += d1
+        asymmetry_list[j] += d2
+    for a in asymmetry_list:
+        a /= len(spike_trains)-1
+    return asymmetry_list
+
+
+############################################################
+# spike_train_order_profile
+############################################################
+def spike_train_order_profile(spike_train1, spike_train2, max_tau=None):
+    """ Computes the spike train order profile P(t) of the two given
+    spike trains. Returns the profile as a DiscreteFunction object.
+    :param spike_train1: First spike train.
+    :type spike_train1: :class:`pyspike.SpikeTrain`
+    :param spike_train2: Second spike train.
+    :type spike_train2: :class:`pyspike.SpikeTrain`
+    :param max_tau: Maximum coincidence window size. If 0 or `None`, the
+                    coincidence window has no upper bound.
+    :returns: The spike-distance profile :math:`S_{sync}(t)`.
+    :rtype: :class:`pyspike.function.DiscreteFunction`
+    """
+    # check whether the spike trains are defined for the same interval
+    assert spike_train1.t_start == spike_train2.t_start, \
+        "Given spike trains are not defined on the same interval!"
+    assert spike_train1.t_end == spike_train2.t_end, \
+        "Given spike trains are not defined on the same interval!"
+
+    # cython implementation
+    try:
+        from .cython.cython_directionality import \
+            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.directionality_python_backend import \
+            spike_train_order_profile_python as spike_train_order_profile_impl
+
+    if max_tau is None:
+        max_tau = 0.0
+
+    times, coincidences, multiplicity \
+        = spike_train_order_profile_impl(spike_train1.spikes,
+                                         spike_train2.spikes,
+                                         spike_train1.t_start,
+                                         spike_train1.t_end,
+                                         max_tau)
+
+    return DiscreteFunc(times, coincidences, multiplicity)
+
+
+############################################################
+# spike_train_order
+############################################################
+def spike_train_order(spike_train1, spike_train2, normalize=True,
+                      interval=None, max_tau=None):
+    """ Computes the overall spike delay asymmetry value for two spike trains.
+    """
+    if interval is None:
+        # distance over the whole interval is requested: use specific function
+        # for optimal performance
+        try:
+            from .cython.cython_directionality import \
+                spike_train_order_cython as spike_train_order_impl
+            if max_tau is None:
+                max_tau = 0.0
+            c, mp = spike_train_order_impl(spike_train1.spikes,
+                                           spike_train2.spikes,
+                                           spike_train1.t_start,
+                                           spike_train1.t_end,
+                                           max_tau)
+        except ImportError:
+            # Cython backend not available: fall back to profile averaging
+            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()
+
+
+############################################################
+# spike_train_order_profile_multi
+############################################################
+def spike_train_order_profile_multi(spike_trains, indices=None,
+                                    max_tau=None):
+    """ Computes the multi-variate spike delay asymmetry profile for a set of
+    spike trains. For each spike in the set of spike trains, the multi-variate
+    profile is defined as the sum of asymmetry values divided by the number of
+    spike trains pairs involving the spike train of containing this spike,
+    which is the number of spike trains minus one (N-1).
+    :param spike_trains: list of :class:`pyspike.SpikeTrain`
+    :param indices: list of indices defining which spike trains to use,
+                    if None all given spike trains are used (default=None)
+    :type indices: list or None
+    :param max_tau: Maximum coincidence window size. If 0 or `None`, the
+                    coincidence window has no upper bound.
+    :returns: The multi-variate spike sync profile :math:`<S_{sync}>(t)`
+    :rtype: :class:`pyspike.function.DiscreteFunction`
+    """
+    prof_func = partial(spike_train_order_profile, max_tau=max_tau)
+    average_prof, M = _generic_profile_multi(spike_trains, prof_func,
+                                             indices)
+    # average_dist.mul_scalar(1.0/M)  # no normalization here!
+    return average_prof
+
+
+############################################################
+# optimal_spike_train_order_from_matrix
+############################################################
+def optimal_spike_train_order_from_matrix(D, full_output=False):
+    """ finds the best sorting via simulated annealing.
+    Returns the optimal permutation p and A value.
+    Internal function, don't call directly! Use optimal_asymmetry_order
+    instead.
+    """
+    N = len(D)
+    A = np.sum(np.triu(D, 0))
+
+    p = np.arange(N)
+
+    T_start = 2*np.max(D)    # starting temperature
+    T_end = 1E-5 * T_start   # final temperature
+    alpha = 0.9              # cooling factor
+
+    from .cython.cython_simulated_annealing import sim_ann_cython as sim_ann
+
+    p, A, total_iter = sim_ann(D, T_start, T_end, alpha)
+
+    # T = T_start
+    # total_iter = 0
+    # while T > T_end:
+    #     iterations = 0
+    #     succ_iter = 0
+    #     # equilibrate for 100*N steps or 10*N successful steps
+    #     while iterations < 100*N and succ_iter < 10*N:
+    #         # exchange two rows and cols
+    #         ind1 = np.random.randint(N-1)
+    #         if ind1 < N-1:
+    #             ind2 = ind1+1
+    #         else:   # this can never happend
+    #             ind2 = 0
+    #         delta_A = -2*D[p[ind1], p[ind2]]
+    #         if delta_A > 0.0 or exp(delta_A/T) > np.random.random():
+    #             # swap indices
+    #             p[ind1], p[ind2] = p[ind2], p[ind1]
+    #             A += delta_A
+    #             succ_iter += 1
+    #         iterations += 1
+    #     total_iter += iterations
+    #     T *= alpha   # cool down
+    #     if succ_iter == 0:
+    #         break
+
+    if full_output:
+        return p, A, total_iter
+    else:
+        return p, A
+
+
+############################################################
+# optimal_spike_train_order
+############################################################
+def optimal_spike_train_order(spike_trains,  indices=None, interval=None,
+                              max_tau=None, full_output=False):
+    """ finds the best sorting of the given spike trains via simulated
+    annealing.
+    Returns the optimal permutation p and A value.
+    """
+    D = spike_directionality_matrix(spike_trains, normalize=False,
+                                    indices=indices, interval=interval,
+                                    max_tau=max_tau)
+    return optimal_spike_train_order_from_matrix(D, full_output)
+
+
+############################################################
+# permutate_matrix
+############################################################
+def permutate_matrix(D, p):
+    """ Applies the permutation p to the columns and rows of matrix D.
+    Return the new permutated matrix.
+    """
+    N = len(D)
+    D_p = np.empty_like(D)
+    for n in range(N):
+        for m in range(N):
+            D_p[n, m] = D[p[n], p[m]]
+    return D_p
diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py
index 80f7805..32e6bf4 100644
--- a/pyspike/spike_sync.py
+++ b/pyspike/spike_sync.py
@@ -8,7 +8,7 @@ from __future__ import absolute_import
 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
 
 
@@ -290,3 +290,52 @@ def spike_sync_matrix(spike_trains, indices=None, interval=None, max_tau=None):
     dist_func = partial(spike_sync_bi, 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,
+                         return_removed_spikes=False):
+    """ Removes the spikes with a multi-variate spike_sync value below
+    threshold.
+    """
+    N = len(spike_trains)
+    filtered_spike_trains = []
+    removed_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_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_single_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):
+            if i == j:
+                continue
+            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]))
+        if return_removed_spikes:
+            removed_spikes = st[coincidences <= threshold*(N-1)]
+            removed_spike_trains.append(SpikeTrain(removed_spikes,
+                                                   [st.t_start, st.t_end]))
+    if return_removed_spikes:
+        return [filtered_spike_trains, removed_spike_trains]
+    else:
+        return filtered_spike_trains
diff --git a/setup.py b/setup.py
index 5b9e677..52bebd9 100644
--- a/setup.py
+++ b/setup.py
@@ -30,7 +30,9 @@ class numpy_include(object):
 
 if os.path.isfile("pyspike/cython/cython_add.c") and \
    os.path.isfile("pyspike/cython/cython_profiles.c") and \
-   os.path.isfile("pyspike/cython/cython_distances.c"):
+   os.path.isfile("pyspike/cython/cython_distances.c") and \
+   os.path.isfile("pyspike/cython/cython_directionality.c") and \
+   os.path.isfile("pyspike/cython/cython_simulated_annealing.c"):
     use_c = True
 else:
     use_c = False
@@ -45,7 +47,11 @@ if use_cython:  # Cython is available, compile .pyx -> .c
         Extension("pyspike.cython.cython_profiles",
                   ["pyspike/cython/cython_profiles.pyx"]),
         Extension("pyspike.cython.cython_distances",
-                  ["pyspike/cython/cython_distances.pyx"])
+                  ["pyspike/cython/cython_distances.pyx"]),
+        Extension("pyspike.cython.cython_directionality",
+                  ["pyspike/cython/cython_directionality.pyx"]),
+        Extension("pyspike.cython.cython_simulated_annealing",
+                  ["pyspike/cython/cython_simulated_annealing.pyx"])
     ]
     cmdclass.update({'build_ext': build_ext})
 elif use_c:  # c files are there, compile to binaries
@@ -55,7 +61,11 @@ elif use_c:  # c files are there, compile to binaries
         Extension("pyspike.cython.cython_profiles",
                   ["pyspike/cython/cython_profiles.c"]),
         Extension("pyspike.cython.cython_distances",
-                  ["pyspike/cython/cython_distances.c"])
+                  ["pyspike/cython/cython_distances.c"]),
+        Extension("pyspike.cython.cython_directionality",
+                  ["pyspike/cython/cython_directionality.c"]),
+        Extension("pyspike.cython.cython_simulated_annealing",
+                  ["pyspike/cython/cython_simulated_annealing.c"])
     ]
 # neither cython nor c files available -> automatic fall-back to python backend
 
@@ -90,11 +100,17 @@ train similarity',
 
         'Programming Language :: Python :: 2',
         'Programming Language :: Python :: 2.7',
-
         'Programming Language :: Python :: 3',
         'Programming Language :: Python :: 3.3',
         'Programming Language :: Python :: 3.4',
         'Programming Language :: Python :: 3.5',
         'Programming Language :: Python :: 3.6'
-    ]
+    ],
+    package_data={
+        'pyspike': ['cython/cython_add.c', 'cython/cython_profiles.c',
+                    'cython/cython_distances.c',
+                    'cython/cython_directionality.c',
+                    'cython/cython_simulated_annealing.c'],
+        'test': ['Spike_testdata.txt']
+    }
 )
diff --git a/test/test_directionality.py b/test/test_directionality.py
new file mode 100644
index 0000000..63865cc
--- /dev/null
+++ b/test/test_directionality.py
@@ -0,0 +1,89 @@
+""" test_directionality.py
+
+Tests the directionality functions
+
+Copyright 2015, Mario Mulansky <mario.mulansky@gmx.net>
+
+Distributed under the BSD License
+
+"""
+
+import numpy as np
+from numpy.testing import assert_equal, assert_almost_equal, \
+    assert_array_equal
+
+import pyspike as spk
+from pyspike import SpikeTrain, DiscreteFunc
+# from pyspike.spike_directionality import _spike_directionality_profile
+
+
+def test_spike_directionality():
+    st1 = SpikeTrain([100, 200, 300], [0, 1000])
+    st2 = SpikeTrain([105, 205, 300], [0, 1000])
+    assert_almost_equal(spk.spike_directionality(st1, st2), 2.0/3.0)
+    assert_almost_equal(spk.spike_directionality(st1, st2, normalize=False),
+                        2.0)
+
+    # exchange order of spike trains should give exact negative profile
+    assert_almost_equal(spk.spike_directionality(st2, st1), -2.0/3.0)
+    assert_almost_equal(spk.spike_directionality(st2, st1, normalize=False),
+                        -2.0)
+
+    st3 = SpikeTrain([105, 195, 500], [0, 1000])
+    assert_almost_equal(spk.spike_directionality(st1, st3), 0.0)
+    assert_almost_equal(spk.spike_directionality(st1, st3, normalize=False),
+                        0.0)
+    assert_almost_equal(spk.spike_directionality(st3, st1), 0.0)
+
+    D = spk.spike_directionality_matrix([st1, st2, st3], normalize=False)
+    D_expected = np.array([[0, 2.0, 0.0], [-2.0, 0.0, -1.0], [0.0, 1.0, 0.0]])
+    assert_array_equal(D, D_expected)
+
+    dir_profs = spk.spike_directionality_profiles([st1, st2, st3])
+    assert_array_equal(dir_profs[0], [1.0, 0.0, 0.0])
+    assert_array_equal(dir_profs[1], [-0.5, -1.0, 0.0])
+
+
+def test_spike_train_order():
+    st1 = SpikeTrain([100, 200, 300], [0, 1000])
+    st2 = SpikeTrain([105, 205, 300], [0, 1000])
+    st3 = SpikeTrain([105, 195, 500], [0, 1000])
+
+    expected_x12 = np.array([0, 100, 105, 200, 205, 300, 1000])
+    expected_y12 = np.array([1, 1, 1, 1, 1, 0, 0])
+    expected_mp12 = np.array([1, 1, 1, 1, 1, 2, 2])
+
+    f = spk.spike_train_order_profile(st1, st2)
+
+    assert f.almost_equal(DiscreteFunc(expected_x12, expected_y12,
+                                       expected_mp12))
+    assert_almost_equal(f.avrg(), 2.0/3.0)
+    assert_almost_equal(f.avrg(normalize=False), 4.0)
+    assert_almost_equal(spk.spike_train_order(st1, st2), 2.0/3.0)
+    assert_almost_equal(spk.spike_train_order(st1, st2, normalize=False), 4.0)
+
+    expected_x23 = np.array([0, 105, 195, 205, 300, 500, 1000])
+    expected_y23 = np.array([0, 0, -1, -1, 0, 0, 0])
+    expected_mp23 = np.array([2, 2, 1, 1, 1, 1, 1])
+
+    f = spk.spike_train_order_profile(st2, st3)
+
+    assert_array_equal(f.x, expected_x23)
+    assert_array_equal(f.y, expected_y23)
+    assert_array_equal(f.mp, expected_mp23)
+    assert f.almost_equal(DiscreteFunc(expected_x23, expected_y23,
+                                       expected_mp23))
+    assert_almost_equal(f.avrg(), -1.0/3.0)
+    assert_almost_equal(f.avrg(normalize=False), -2.0)
+    assert_almost_equal(spk.spike_train_order(st2, st3), -1.0/3.0)
+    assert_almost_equal(spk.spike_train_order(st2, st3, normalize=False), -2.0)
+
+    f = spk.spike_train_order_profile_multi([st1, st2, st3])
+
+    expected_x = np.array([0, 100, 105, 195, 200, 205, 300, 500, 1000])
+    expected_y = np.array([2, 2, 2, -2, 0, 0, 0, 0, 0])
+    expected_mp = np.array([2, 2, 4, 2, 2, 2, 4, 2, 2])
+
+    assert_array_equal(f.x, expected_x)
+    assert_array_equal(f.y, expected_y)
+    assert_array_equal(f.mp, expected_mp)
diff --git a/test/test_sync_filter.py b/test/test_sync_filter.py
new file mode 100644
index 0000000..e259903
--- /dev/null
+++ b/test/test_sync_filter.py
@@ -0,0 +1,95 @@
+""" test_sync_filter.py
+
+Tests the spike sync based filtering
+
+Copyright 2015, Mario Mulansky <mario.mulansky@gmx.net>
+
+Distributed under the BSD License
+
+"""
+
+from __future__ import print_function
+import numpy as np
+from numpy.testing import assert_equal, assert_almost_equal, \
+    assert_array_almost_equal
+
+import pyspike as spk
+from pyspike import SpikeTrain
+
+
+def test_single_prof():
+    st1 = np.array([1.0, 2.0, 3.0, 4.0])
+    st2 = np.array([1.1, 2.1, 3.8])
+    st3 = np.array([0.9, 3.1, 4.1])
+
+    # cython implementation
+    try:
+        from pyspike.cython.cython_profiles import \
+            coincidence_single_profile_cython as coincidence_impl
+    except ImportError:
+        from pyspike.cython.python_backend import \
+            coincidence_single_python as coincidence_impl
+
+    sync_prof = spk.spike_sync_profile(SpikeTrain(st1, 5.0),
+                                       SpikeTrain(st2, 5.0))
+
+    coincidences = np.array(coincidence_impl(st1, st2, 0, 5.0, 0.0))
+    print(coincidences)
+    for i, t in enumerate(st1):
+        assert_equal(coincidences[i], sync_prof.y[sync_prof.x == t],
+                     "At index %d" % i)
+
+    coincidences = np.array(coincidence_impl(st2, st1, 0, 5.0, 0.0))
+    for i, t in enumerate(st2):
+        assert_equal(coincidences[i], sync_prof.y[sync_prof.x == t],
+                     "At index %d" % i)
+
+    sync_prof = spk.spike_sync_profile(SpikeTrain(st1, 5.0),
+                                       SpikeTrain(st3, 5.0))
+
+    coincidences = np.array(coincidence_impl(st1, st3, 0, 5.0, 0.0))
+    for i, t in enumerate(st1):
+        assert_equal(coincidences[i], sync_prof.y[sync_prof.x == t],
+                     "At index %d" % i)
+
+    st1 = np.array([1.0, 2.0, 3.0, 4.0])
+    st2 = np.array([1.0, 2.0, 4.0])
+
+    sync_prof = spk.spike_sync_profile(SpikeTrain(st1, 5.0),
+                                       SpikeTrain(st2, 5.0))
+
+    coincidences = np.array(coincidence_impl(st1, st2, 0, 5.0, 0.0))
+    for i, t in enumerate(st1):
+        expected = sync_prof.y[sync_prof.x == t]/sync_prof.mp[sync_prof.x == t]
+        assert_equal(coincidences[i], expected,
+                     "At index %d" % i)
+
+
+def test_filter():
+    st1 = SpikeTrain(np.array([1.0, 2.0, 3.0, 4.0]), 5.0)
+    st2 = SpikeTrain(np.array([1.1, 2.1, 3.8]), 5.0)
+    st3 = SpikeTrain(np.array([0.9, 3.1, 4.1]), 5.0)
+
+    # filtered_spike_trains = spk.filter_by_spike_sync([st1, st2], 0.5)
+
+    # assert_equal(filtered_spike_trains[0].spikes, [1.0, 2.0, 4.0])
+    # assert_equal(filtered_spike_trains[1].spikes, [1.1, 2.1, 3.8])
+
+    # filtered_spike_trains = spk.filter_by_spike_sync([st2, st1], 0.5)
+
+    # assert_equal(filtered_spike_trains[0].spikes, [1.1, 2.1, 3.8])
+    # assert_equal(filtered_spike_trains[1].spikes, [1.0, 2.0, 4.0])
+
+    filtered_spike_trains = spk.filter_by_spike_sync([st1, st2, st3], 0.75)
+
+    for st in filtered_spike_trains:
+        print(st.spikes)
+
+    assert_equal(filtered_spike_trains[0].spikes, [1.0, 4.0])
+    assert_equal(filtered_spike_trains[1].spikes, [1.1, 3.8])
+    assert_equal(filtered_spike_trains[2].spikes, [0.9, 4.1])
+
+
+if __name__ == "main":
+    test_single_prof()
+    test_filter()