spike train synchrony: new function called spike contrast (NeuralEnse…

…mble#354)

.gitignore

@@ -44,11 +44,9 @@ lib
 lib64
 # sphinx build directory
 doc/_build
-doc/reference/toctree/asset/elephant.asset.synchronous*
-doc/reference/toctree/parallel
-doc/reference/toctree/statistics
-doc/reference/toctree/unitary_event_analysis
-doc/reference/toctree/causality
+doc/reference/toctree/*
+!doc/reference/toctree/asset/elephant.asset.ASSET.rst
+!doc/reference/toctree/kernels
 *.h5
 # setup.py dist directory
 dist

doc/authors.rst

@@ -49,6 +49,8 @@ Do you want to contribute to Elephant? Please refer to the
 * Danylo Ulianych [1]
 * Anno Kurth [1]
 * Regimantas Jurkus [1]
+* Philipp Steigerwald [12]
+* Manuel Ciba [12]
 
 1. Institute of Neuroscience and Medicine (INM-6), Computational and Systems Neuroscience & Institute for Advanced Simulation (IAS-6), Theoretical Neuroscience, Jülich Research Centre and JARA, Jülich, Germany
 2. Unité de Neurosciences, Information et Complexité, CNRS UPR 3293, Gif-sur-Yvette, France
@@ -61,5 +63,6 @@ Do you want to contribute to Elephant? Please refer to the
 9. Nencki Institute of Experimental Biology, Warsaw, Poland
 10. Instituto de Neurobiología, Universidad Nacional Autónoma de México, Mexico City, Mexico
 11. Case Western Reserve University (CWRU), Cleveland, OH, USA
+12. BioMEMS Lab, TH Aschaffenburg University of applied sciences, Germany
 
 If we've somehow missed you off the list we're very sorry - please let us know.

doc/bib/elephant.bib

@@ -115,6 +115,16 @@ @article{Torre16_e1004939
   doi={10.1371/journal.pcbi.1004939}
 }
 
+@article{Ciba18_136,
+  title={Spike-contrast: A novel time scale independent and multivariate measure of spike train synchrony},
+  author={Ciba, M. and Isomura, T. and Jimbo, Y. and Bahmer, A. and Thielemann, C.},
+  year={2018},
+  journal={J. Neurosci. Meth.},
+  volume={293},
+  pages={136--143},
+  doi={10.1016/j.jneumeth.2017.09.008}
+}
+
 @article{Rostami17_3,
   title = {[{Re}] {S}pike synchronization and rate modulation differentially involved in motor cortical function},
   author = {Rostami, V. and Ito, J. and Denker, M. and Gr{\"u}n, S.},
@@ -125,4 +135,3 @@ @article{Rostami17_3
   pages = {3},
   doi = {10.5281/zenodo.583814}
 }
-

doc/modules.rst

@@ -21,11 +21,15 @@ Function Reference by Module
    reference/signal_processing
    reference/spade
    reference/spectral
-   reference/spike_train_correlation
-   reference/spike_train_dissimilarity
    reference/spike_train_generation
    reference/spike_train_surrogates
    reference/sta
    reference/statistics
    reference/unitary_event_analysis
    reference/waveform_features
+
+
+.. toctree::
+   :maxdepth: 2
+
+   reference/_spike_train_processing

doc/reference/_spike_train_processing.rst

@@ -0,0 +1,11 @@
+======================
+Spike train processing
+======================
+
+
+.. toctree::
+   :maxdepth: 1
+
+   spike_train_correlation
+   spike_train_dissimilarity
+   spike_train_synchrony

doc/reference/spike_train_dissimilarity.rst

@@ -1,6 +1,6 @@
-=======================================
-Spike train dissimilarity and synchrony
-=======================================
+=========================
+Spike train dissimilarity
+=========================
 
 
 .. automodule:: elephant.spike_train_dissimilarity

doc/reference/spike_train_synchrony.rst

@@ -0,0 +1,14 @@
+=====================
+Spike train synchrony
+=====================
+
+.. automodule:: elephant.spike_train_synchrony
+
+
+References
+----------
+
+.. bibliography:: ../bib/elephant.bib
+   :labelprefix: syn
+   :keyprefix: synchrony-
+   :style: unsrt

elephant/spike_train_synchrony.py

@@ -0,0 +1,218 @@
+# -*- coding: utf-8 -*-
+"""
+Functions to measure the synchrony of several spike trains.
+
+
+Synchrony Measures
+------------------
+
+.. autosummary::
+    :toctree: toctree/spike_train_synchrony/
+
+    spike_contrast
+
+
+:copyright: Copyright 2015-2020 by the Elephant team, see `doc/authors.rst`.
+:license: Modified BSD, see LICENSE.txt for details.
+"""
+from __future__ import division, print_function, unicode_literals
+
+from collections import namedtuple
+
+import neo
+import numpy as np
+import quantities as pq
+
+from elephant.utils import is_time_quantity
+
+SpikeContrastTrace = namedtuple("SpikeContrastTrace", (
+    "contrast", "active_spiketrains", "synchrony"))
+
+
+def _get_theta_and_n_per_bin(spiketrains, t_start, t_stop, bin_size):
+    """
+    Calculates theta (amount of spikes per bin) and the amount of active spike
+    trains per bin of one spike train.
+    """
+    bin_step = bin_size / 2
+    edges = np.arange(t_start, t_stop + bin_step, bin_step)
+    # Calculate histogram for every spike train
+    histogram = np.vstack([
+        _binning_half_overlap(st, edges=edges)
+        for st in spiketrains
+    ])
+    # Amount of spikes per bin
+    theta = histogram.sum(axis=0)
+    # Amount of active spike trains per bin
+    n_active_per_bin = np.count_nonzero(histogram, axis=0)
+
+    return theta, n_active_per_bin
+
+
+def _binning_half_overlap(spiketrain, edges):
+    """
+    Referring to [1] overlapping the bins creates a better result.
+    """
+    histogram, bin_edges = np.histogram(spiketrain, bins=edges)
+    histogram = histogram[:-1] + histogram[1:]
+    return histogram
+
+
+def spike_contrast(spiketrains, t_start=None, t_stop=None,
+                   min_bin=10 * pq.ms, bin_shrink_factor=0.9,
+                   return_trace=False):
+    """
+    Calculates the synchrony of spike trains, according to
+    :cite:`synchrony-Ciba18_136`. The spike trains can have different lengths.
+
+    Original implementation by: Philipp Steigerwald [s160857@th-ab.de]
+
+    Parameters
+    ----------
+    spiketrains : list of neo.SpikeTrain
+        A list of input spike trains to calculate the synchrony from.
+    t_start : pq.Quantity, optional
+        The beginning of the spike train. If None, it's taken as the minimum
+        value of `t_start`s of the input spike trains.
+        Default: None
+    t_stop : pq.Quantity, optional
+        The end of the spike train. If None, it's taken as the maximum value
+        of `t_stop` of the input spike trains.
+        Default: None
+    min_bin : pq.Quantity, optional
+        Sets the minimum value for the `bin_min` that is calculated by the
+        algorithm and defines the smallest bin size to compute the histogram
+        of the input `spiketrains`.
+        Default: 0.01 ms
+    bin_shrink_factor : float, optional
+        A multiplier to shrink the bin size on each iteration. The value must
+        be in range `(0, 1)`.
+        Default: 0.9
+    return_trace : bool, optional
+        If set to True, returns a history of spike-contrast synchrony, computed
+        for a range of different bin sizes, alongside with the maximum value of
+        the synchrony.
+        Default: False
+
+    Returns
+    -------
+    synchrony : float
+        Returns the synchrony of the input spike trains.
+    spike_contrast_trace : namedtuple
+        If `return_trace` is set to True, a `SpikeContrastTrace` namedtuple is
+        returned with the following attributes:
+          `.contrast` - the average sum of differences of the number of spikes
+          in subsuequent bins;
+
+          `.active_spiketrains` - the average number of spikes per bin,
+          weighted by the number of spike trains containing at least one spike
+          inside the bin;
+
+          `.synchrony` - the product of `contrast` and `active_spiketrains`.
+
+    Raises
+    ------
+    ValueError
+        If `bin_shrink_factor` is not in (0, 1) range.
+
+        If the input spike trains constist of a single spiketrain.
+
+        If all input spike trains contain no more than 1 spike.
+    TypeError
+        If the input spike trains is not a list of `neo.SpikeTrain` objects.
+
+        If `t_start`, `t_stop`, or `min_bin` are not time quantities.
+
+    Examples
+    --------
+    >>> import quantities as pq
+    >>> from elephant.spike_train_generation import homogeneous_poisson_process
+    >>> from elephant.spike_train_synchrony import spike_contrast
+    >>> spiketrain_1 = homogeneous_poisson_process(rate=20*pq.Hz,
+    ...     t_stop=1000*pq.ms)
+    >>> spiketrain_2 = homogeneous_poisson_process(rate=20*pq.Hz,
+    ...     t_stop=1000*pq.ms)
+    >>> spike_contrast([spiketrain_1, spiketrain_2])
+    0.4192546583850932
+
+    """
+    if not 0. < bin_shrink_factor < 1.:
+        raise ValueError("'bin_shrink_factor' ({}) must be in range (0, 1)."
+                         .format(bin_shrink_factor))
+    if not len(spiketrains) > 1:
+        raise ValueError("Spike contrast measure requires more than 1 input "
+                         "spiketrain.")
+    if not all(isinstance(st, neo.SpikeTrain) for st in spiketrains):
+        raise TypeError("Input spike trains must be a list of neo.SpikeTrain.")
+    if not is_time_quantity(t_start, allow_none=True) \
+            or not is_time_quantity(t_stop, allow_none=True):
+        raise TypeError("'t_start' and 't_stop' must be time quantities.")
+    if not is_time_quantity(min_bin):
+        raise TypeError("'min_bin' must be a time quantity.")
+
+    if t_start is None:
+        t_start = min(st.t_start for st in spiketrains)
+    if t_stop is None:
+        t_stop = max(st.t_stop for st in spiketrains)
+    spiketrains = [st.time_slice(t_start=t_start, t_stop=t_stop)
+                   for st in spiketrains]
+
+    # convert everything to seconds
+    spiketrains = [st.simplified.magnitude for st in spiketrains]
+    t_start = t_start.simplified.item()
+    t_stop = t_stop.simplified.item()
+    min_bin = min_bin.simplified.item()
+
+    n_spiketrains = len(spiketrains)
+    n_spikes_total = sum(map(len, spiketrains))
+
+    duration = t_stop - t_start
+    bin_max = duration / 2
+
+    try:
+        isi_min = min(np.diff(st).min() for st in spiketrains if len(st) > 1)
+    except TypeError:
+        raise ValueError("All input spiketrains contain no more than 1 spike.")
+    bin_min = max(isi_min / 2, min_bin)
+
+    contrast_list = []
+    active_spiketrains = []
+    synchrony_curve = []
+
+    bin_size = bin_max
+    while bin_size >= bin_min:
+        # Set new time boundaries
+        t_start = -isi_min
+        t_stop = duration + isi_min
+        # Calculate Theta and n
+        theta_k, n_k = _get_theta_and_n_per_bin(spiketrains,
+                                                t_start=t_start,
+                                                t_stop=t_stop,
+                                                bin_size=bin_size)
+
+        # calculate synchrony_curve = contrast * active_st
+        active_st = (np.sum(n_k * theta_k) / np.sum(theta_k) - 1) / (
+                    n_spiketrains - 1)
+        contrast = np.sum(np.abs(np.diff(theta_k))) / (2 * n_spikes_total)
+        # Contrast: sum(|derivation|) / (2*#Spikes)
+        synchrony = contrast * active_st
+
+        contrast_list.append(contrast)
+        active_spiketrains.append(active_st)
+        synchrony_curve.append(synchrony)
+
+        # New bin size
+        bin_size *= bin_shrink_factor
+
+    # Sync value is maximum of the cost function C
+    synchrony = max(synchrony_curve)
+
+    if return_trace:
+        spike_contrast_trace = SpikeContrastTrace(
+            contrast=contrast_list,
+            active_spiketrains=active_spiketrains,
+            synchrony=synchrony_curve
+        )
+        return synchrony, spike_contrast_trace
+
+    return synchrony