# Module containing functions to compute the SPIKE directionality and the # spike train order profile # Copyright 2015, Mario Mulansky # 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_values ############################################################ def spike_directionality_values(*args, **kwargs): """ Computes the spike directionality value for each spike in each spike train. Returns a list containing an array of spike directionality values for every given spike train. Valid call structures:: spike_directionality_values(st1, st2) # returns the bi-variate profile spike_directionality_values(st1, st2, st3) # multi-variate profile of 3 # spike trains spike_trains = [st1, st2, st3, st4] # list of spike trains spike_directionality_values(spike_trains) # profile of the list of spike trains spike_directionality_values(spike_trains, indices=[0, 1]) # use only the spike trains # given by the indices Additonal arguments: :param max_tau: Upper bound for coincidence window (default=None). :param indices: list of indices defining which spike trains to use, if None all given spike trains are used (default=None) :returns: The spike directionality values :math:`D^n_i` as a list of arrays. """ if len(args) == 1: return _spike_directionality_values_impl(args[0], **kwargs) else: return _spike_directionality_values_impl(args, **kwargs) def _spike_directionality_values_impl(spike_trains, indices=None, interval=None, max_tau=None): """ Computes the multi-variate spike directionality profile of the given spike trains. :param spike_trains: List of spike trains. :type 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 spike-directionality values. """ if interval is not None: raise NotImplementedError("Parameter `interval` not supported.") 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 resulting asymmetry values asymmetry_list = [np.zeros_like(spike_trains[n].spikes) for n in indices] # 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_directionality ############################################################ def spike_directionality(spike_train1, spike_train2, normalize=True, interval=None, max_tau=None): """ Computes the overall spike directionality of the first spike train with respect to the second spike train. :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 normalize: Normalize by the number of spikes (multiplicity). :param max_tau: Maximum coincidence window size. If 0 or `None`, the coincidence window has no upper bound. :returns: The spike train order profile :math:`E(t)`. """ 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: 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 profile. d1, x = spike_directionality_values([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("Parameter `interval` not supported.") ############################################################ # spike_directionality_matrix ############################################################ def spike_directionality_matrix(spike_trains, normalize=True, indices=None, interval=None, max_tau=None): """ Computes the spike directionality matrix for the given spike trains. :param spike_trains: List of spike trains. :type spike_trains: List of :class:`pyspike.SpikeTrain` :param normalize: Normalize by the number of spikes (multiplicity). :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 spike-directionality values. """ 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_train_order_profile ############################################################ def spike_train_order_profile(*args, **kwargs): """ Computes the spike train order profile :math:`E(t)` of the given spike trains. Returns the profile as a DiscreteFunction object. Valid call structures:: spike_train_order_profile(st1, st2) # returns the bi-variate profile spike_train_order_profile(st1, st2, st3) # multi-variate profile of 3 # spike trains spike_trains = [st1, st2, st3, st4] # list of spike trains spike_train_order_profile(spike_trains) # profile of the list of spike trains spike_train_order_profile(spike_trains, indices=[0, 1]) # use only the spike trains # given by the indices Additonal arguments: :param max_tau: Upper bound for coincidence window, `default=None`. :param indices: list of indices defining which spike trains to use, if None all given spike trains are used (default=None) :returns: The spike train order profile :math:`E(t)` :rtype: :class:`.DiscreteFunction` """ if len(args) == 1: return spike_train_order_profile_multi(args[0], **kwargs) elif len(args) == 2: return spike_train_order_profile_bi(args[0], args[1], **kwargs) else: return spike_train_order_profile_multi(args, **kwargs) ############################################################ # spike_train_order_profile_bi ############################################################ def spike_train_order_profile_bi(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 train order profile :math:`E(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_profile_multi ############################################################ def spike_train_order_profile_multi(spike_trains, indices=None, max_tau=None): """ Computes the multi-variate spike train order 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:`(t)` :rtype: :class:`pyspike.function.DiscreteFunction` """ prof_func = partial(spike_train_order_profile_bi, max_tau=max_tau) average_prof, M = _generic_profile_multi(spike_trains, prof_func, indices) return average_prof ############################################################ # _spike_train_order_impl ############################################################ def _spike_train_order_impl(spike_train1, spike_train2, interval=None, max_tau=None): """ Implementation of bi-variatae spike train order value (Synfire Indicator). :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 train order value (Synfire Indicator) """ 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_func if max_tau is None: max_tau = 0.0 c, mp = spike_train_order_func(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=max_tau).integral(interval) return c, mp else: # some specific interval is provided: not yet implemented raise NotImplementedError("Parameter `interval` not supported.") ############################################################ # spike_train_order ############################################################ def spike_train_order(*args, **kwargs): """ Computes the spike train order (Synfire Indicator) of the given spike trains. Valid call structures:: spike_train_order(st1, st2, normalize=True) # normalized bi-variate # spike train order spike_train_order(st1, st2, st3) # multi-variate result of 3 spike trains spike_trains = [st1, st2, st3, st4] # list of spike trains spike_train_order(spike_trains) # result for the list of spike trains spike_train_order(spike_trains, indices=[0, 1]) # use only the spike trains # given by the indices Additonal arguments: - `max_tau` Upper bound for coincidence window, `default=None`. - `normalize` Flag indicating if the reslut should be normalized by the number of spikes , default=`False` :returns: The spike train order value (Synfire Indicator) """ if len(args) == 1: return spike_train_order_multi(args[0], **kwargs) elif len(args) == 2: return spike_train_order_bi(args[0], args[1], **kwargs) else: return spike_train_order_multi(args, **kwargs) ############################################################ # spike_train_order_bi ############################################################ def spike_train_order_bi(spike_train1, spike_train2, normalize=True, interval=None, max_tau=None): """ Computes the overall spike train order value (Synfire Indicator) for two spike trains. :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 normalize: Normalize by the number of spikes (multiplicity). :param max_tau: Maximum coincidence window size. If 0 or `None`, the coincidence window has no upper bound. :returns: The spike train order value (Synfire Indicator) """ c, mp = _spike_train_order_impl(spike_train1, spike_train2, interval, max_tau) if normalize: return 1.0*c/mp else: return c ############################################################ # spike_train_order_multi ############################################################ def spike_train_order_multi(spike_trains, indices=None, normalize=True, interval=None, max_tau=None): """ Computes the overall spike train order value (Synfire Indicator) for many spike trains. :param spike_trains: list of :class:`.SpikeTrain` :param indices: list of indices defining which spike trains to use, if None all given spike trains are used (default=None) :param normalize: Normalize by the number of spike (multiplicity). :param interval: averaging interval given as a pair of floats, if None the average over the whole function is computed. :type interval: Pair of floats or None. :param max_tau: Maximum coincidence window size. If 0 or `None`, the coincidence window has no upper bound. :returns: Spike train order values (Synfire Indicator) F for the given spike trains. :rtype: double """ 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:]] e_total = 0.0 m_total = 0.0 for (i, j) in pairs: e, m = _spike_train_order_impl(spike_trains[i], spike_trains[j], interval, max_tau) e_total += e m_total += m if m == 0.0: return 1.0 else: return e_total/m_total ############################################################ # optimal_spike_train_sorting_from_matrix ############################################################ def _optimal_spike_train_sorting_from_matrix(D, full_output=False): """ Finds the best sorting via simulated annealing. Returns the optimal permutation p and A value. Not for direct use, call :func:`.optimal_spike_train_sorting` instead. :param D: The directionality (Spike-ORDER) matrix. :param full_output: If true, then function will additionally return the number of performed iterations (default=False) :return: (p, F) - tuple with the optimal permutation and synfire indicator. if `full_output=True` , (p, F, iter) is returned. """ 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 try: from .cython.cython_simulated_annealing import sim_ann_cython as sim_ann except ImportError: raise NotImplementedError("PySpike with Cython required for computing spike train" " sorting!") p, A, total_iter = sim_ann(D, T_start, T_end, alpha) if full_output: return p, A, total_iter else: return p, A ############################################################ # optimal_spike_train_sorting ############################################################ def optimal_spike_train_sorting(spike_trains, indices=None, interval=None, max_tau=None, full_output=False): """ Finds the best sorting of the given spike trains by computing the spike directionality matrix and optimize the order using simulated annealing. For a detailed description of the algorithm see: `http://iopscience.iop.org/article/10.1088/1367-2630/aa68c3/meta` :param spike_trains: list of :class:`.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 interval: time interval filter given as a pair of floats, if None the full spike trains are used (default=None). :type interval: Pair of floats or None. :param max_tau: Maximum coincidence window size. If 0 or `None`, the coincidence window has no upper bound (default=None). :param full_output: If true, then function will additionally return the number of performed iterations (default=False) :return: (p, F) - tuple with the optimal permutation and synfire indicator. if `full_output=True` , (p, F, iter) is returned. """ D = spike_directionality_matrix(spike_trains, normalize=False, indices=indices, interval=interval, max_tau=max_tau) return _optimal_spike_train_sorting_from_matrix(D, full_output) ############################################################ # permutate_matrix ############################################################ def permutate_matrix(D, p): """ Helper function that applies the permutation p to the columns and rows of matrix D. Return the permutated matrix :math:`D'[n,m] = D[p[n], p[m]]`. :param D: The matrix. :param d: The permutation. :return: The permuated matrix D', ie :math:`D'[n,m] = D[p[n], p[m]]` """ 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