new spike delay asymmetry measures

added first version of spike delay asymmetry functions. still incomplete and
untested.

7 files changed, 427 insertions, 11 deletions

diff --git a/pyspike/__init__.py b/pyspike/__init__.py
index 2060f73..8d92ea4 100644
--- a/pyspike/__init__.py
+++ b/pyspike/__init__.py
@@ -6,7 +6,7 @@ Distributed under the BSD License
 
 __all__ = ["isi_distance", "spike_distance", "spike_sync", "psth",
            "spikes", "SpikeTrain", "PieceWiseConstFunc", "PieceWiseLinFunc",
-           "DiscreteFunc"]
+           "DiscreteFunc", "directionality"]
 
 from PieceWiseConstFunc import PieceWiseConstFunc
 from PieceWiseLinFunc import PieceWiseLinFunc
@@ -24,6 +24,7 @@ from psth import psth
 from spikes import load_spike_trains_from_txt, spike_train_from_string, \
     merge_spike_trains, generate_poisson_spikes
 
+import directionality as drct
 
 # define the __version__ following
 # http://stackoverflow.com/questions/17583443
diff --git a/pyspike/directionality/__init__.py b/pyspike/directionality/__init__.py
new file mode 100644
index 0000000..e6de1de
--- /dev/null
+++ b/pyspike/directionality/__init__.py
@@ -0,0 +1,12 @@
+"""
+Copyright 2015, Mario Mulansky <mario.mulansky@gmx.net>
+
+Distributed under the BSD License
+"""
+
+__all__ = ["spike_delay_asymmetry"]
+
+from spike_delay_asymmetry import spike_delay_asymmetry_profile, \
+    spike_delay_asymmetry, spike_delay_asymmetry_profile_multi, \
+    spike_delay_asymmetry_matrix, optimal_asymmetry_order, \
+    optimal_asymmetry_order_from_D, reorder_asymmetry_matrix
diff --git a/pyspike/directionality/cython/__init__.py b/pyspike/directionality/cython/__init__.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/pyspike/directionality/cython/__init__.py
diff --git a/pyspike/directionality/cython/cython_directionality.pyx b/pyspike/directionality/cython/cython_directionality.pyx
new file mode 100644
index 0000000..f5ea752
--- /dev/null
+++ b/pyspike/directionality/cython/cython_directionality.pyx
@@ -0,0 +1,177 @@
+#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_delay_asymmetry_profile_cython
+############################################################
+def spike_delay_asymmetry_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_delay_asymmetry_cython
+############################################################
+def spike_delay_asymmetry_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 asym = 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
+                asym -= 1
+        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
+                asym += 1
+        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 asym == 0 and mp == 0:
+        # empty spike trains -> spike sync = 1 by definition
+        asym = 1
+        mp = 1
+
+    return asym, mp
diff --git a/pyspike/directionality/spike_delay_asymmetry.py b/pyspike/directionality/spike_delay_asymmetry.py
new file mode 100644
index 0000000..7d59601
--- /dev/null
+++ b/pyspike/directionality/spike_delay_asymmetry.py
@@ -0,0 +1,212 @@
+# Module containing functions to compute multivariate spike delay asymmetry
+# Copyright 2015, Mario Mulansky <mario.mulansky@gmx.net>
+# Distributed under the BSD License
+
+import numpy as np
+from math import exp
+from functools import partial
+# import pyspike
+from pyspike import DiscreteFunc
+from pyspike.generic import _generic_profile_multi
+
+
+############################################################
+# spike_delay_asymmetry_profile
+############################################################
+def spike_delay_asymmetry_profile(spike_train1, spike_train2, max_tau=None):
+    """ Computes the spike delay asymmetry profile A(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_delay_asymmetry_profile_cython as \
+            spike_delay_asymmetry_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 max_tau is None:
+        max_tau = 0.0
+
+    times, coincidences, multiplicity \
+        = spike_delay_asymmetry_profile_impl(spike_train1.spikes,
+                                             spike_train2.spikes,
+                                             spike_train1.t_start,
+                                             spike_train1.t_end,
+                                             max_tau)
+
+    return DiscreteFunc(times, coincidences, multiplicity)
+
+
+############################################################
+# spike_delay_asymmetry
+############################################################
+def spike_delay_asymmetry(spike_train1, spike_train2,
+                          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_delay_asymmetry_cython as spike_delay_impl
+            if max_tau is None:
+                max_tau = 0.0
+            c, mp = spike_delay_impl(spike_train1.spikes,
+                                     spike_train2.spikes,
+                                     spike_train1.t_start,
+                                     spike_train1.t_end,
+                                     max_tau)
+            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)
+    else:
+        # some specific interval is provided: not yet implemented
+        raise NotImplementedError()
+
+
+############################################################
+# spike_delay_asymmetry_profile_multi
+############################################################
+def spike_delay_asymmetry_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_delay_asymmetry_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
+
+
+############################################################
+# spike_delay_asymmetry_matrix
+############################################################
+def spike_delay_asymmetry_matrix(spike_trains, indices=None,
+                                 interval=None, max_tau=None):
+    """ Computes the spike delay asymmetry 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_delay_asymmetry(spike_trains[i], spike_trains[j],
+                                  interval, max_tau=max_tau)
+        distance_matrix[i, j] = d
+        distance_matrix[j, i] = -d
+    return distance_matrix
+
+
+############################################################
+# optimal_asymmetry_order_from_D
+############################################################
+def optimal_asymmetry_order_from_D(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 = 2*np.max(D)    # starting temperature
+    T_end = 1E-5 * T   # final temperature
+    alpha = 0.9        # cooling factor
+    total_iter = 0
+    while T > T_end:
+        iterations = 0
+        succ_iter = 0
+        while iterations < 100*N and succ_iter < 10*N:
+            # exchange two rows and cols
+            ind1 = np.random.randint(N-1)
+            delta_A = -2*D[p[ind1], p[ind1+1]]
+            if delta_A > 0.0 or exp(delta_A/T) > np.random.random():
+                # swap indices
+                p[ind1], p[ind1+1] = p[ind1+1], 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_asymmetry_order
+############################################################
+def optimal_asymmetry_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_delay_asymmetry_matrix(spike_trains, indices, interval, max_tau)
+    return optimal_asymmetry_order_from_D(D, full_output)
+
+
+############################################################
+# reorder_asymmetry_matrix
+############################################################
+def reorder_asymmetry_matrix(D, p):
+    N = len(D)
+    D_p = np.empty_like(D)
+    for n in xrange(N):
+        for m in xrange(N):
+            D_p[n, m] = D[p[n], p[m]]
+    return D_p
diff --git a/pyspike/generic.py b/pyspike/generic.py
index 2df34f1..515cbf4 100644
--- a/pyspike/generic.py
+++ b/pyspike/generic.py
@@ -37,13 +37,15 @@ def _generic_profile_multi(spike_trains, pair_distance_func, indices=None):
         """
         L1 = len(pairs1)
         if L1 > 1:
-            dist_prof1 = divide_and_conquer(pairs1[:L1//2], pairs1[int(L1//2):])
+            dist_prof1 = divide_and_conquer(pairs1[:L1//2],
+                                            pairs1[int(L1//2):])
         else:
             dist_prof1 = pair_distance_func(spike_trains[pairs1[0][0]],
                                             spike_trains[pairs1[0][1]])
         L2 = len(pairs2)
         if L2 > 1:
-            dist_prof2 = divide_and_conquer(pairs2[:L2//2], pairs2[int(L2//2):])
+            dist_prof2 = divide_and_conquer(pairs2[:L2//2],
+                                            pairs2[int(L2//2):])
         else:
             dist_prof2 = pair_distance_func(spike_trains[pairs2[0][0]],
                                             spike_trains[pairs2[0][1]])
diff --git a/setup.py b/setup.py
index d853cdf..960c684 100644
--- a/setup.py
+++ b/setup.py
@@ -23,7 +23,8 @@ else:
 
 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/directionality/cython/cython_directionality.c"):
     use_c = True
 else:
     use_c = False
@@ -33,16 +34,26 @@ ext_modules = []
 
 if use_cython:  # Cython is available, compile .pyx -> .c
     ext_modules += [
-        Extension("pyspike.cython.cython_add", ["pyspike/cython/cython_add.pyx"]),
-        Extension("pyspike.cython.cython_profiles", ["pyspike/cython/cython_profiles.pyx"]),
-        Extension("pyspike.cython.cython_distances", ["pyspike/cython/cython_distances.pyx"]),
+        Extension("pyspike.cython.cython_add",
+                  ["pyspike/cython/cython_add.pyx"]),
+        Extension("pyspike.cython.cython_profiles",
+                  ["pyspike/cython/cython_profiles.pyx"]),
+        Extension("pyspike.cython.cython_distances",
+                  ["pyspike/cython/cython_distances.pyx"]),
+        Extension("pyspike.directionality.cython.cython_directionality",
+                  ["pyspike/directionality/cython/cython_directionality.pyx"])
     ]
     cmdclass.update({'build_ext': build_ext})
 elif use_c:  # c files are there, compile to binaries
     ext_modules += [
-        Extension("pyspike.cython.cython_add", ["pyspike/cython/cython_add.c"]),
-        Extension("pyspike.cython.cython_profiles", ["pyspike/cython/cython_profiles.c"]),
-        Extension("pyspike.cython.cython_distances", ["pyspike/cython/cython_distances.c"]),
+        Extension("pyspike.cython.cython_add",
+                  ["pyspike/cython/cython_add.c"]),
+        Extension("pyspike.cython.cython_profiles",
+                  ["pyspike/cython/cython_profiles.c"]),
+        Extension("pyspike.cython.cython_distances",
+                  ["pyspike/cython/cython_distances.c"]),
+        Extension("pyspike.directionality.cython.cython_directionality",
+                  ["pyspike/directionality/cython/cython_directionality.c"])
     ]
 # neither cython nor c files available -> automatic fall-back to python backend
 
@@ -81,7 +92,8 @@ train similarity',
     ],
     package_data={
         'pyspike': ['cython/cython_add.c', 'cython/cython_profiles.c',
-                    'cython_distances.c'],
+                    'cython/cython_distances.c',
+                    'directionality/cython/cython_directionality.c'],
         'test': ['Spike_testdata.txt']
     }
 )