12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157 |
- # -*- coding: utf-8 -*-
- """
- Statistical measures of spike trains (e.g., Fano factor) and functions to estimate firing rates.
- :copyright: Copyright 2014-2016 by the Elephant team, see AUTHORS.txt.
- :license: Modified BSD, see LICENSE.txt for details.
- """
- from __future__ import division, print_function
- import numpy as np
- import quantities as pq
- import scipy.stats
- import scipy.signal
- import neo
- from neo.core import SpikeTrain
- import elephant.conversion as conv
- import elephant.kernels as kernels
- import warnings
- # warnings.simplefilter('always', DeprecationWarning)
- def isi(spiketrain, axis=-1):
- """
- Return an array containing the inter-spike intervals of the SpikeTrain.
- Accepts a Neo SpikeTrain, a Quantity array, or a plain NumPy array.
- If either a SpikeTrain or Quantity array is provided, the return value will
- be a quantities array, otherwise a plain NumPy array. The units of
- the quantities array will be the same as spiketrain.
- Parameters
- ----------
- spiketrain : Neo SpikeTrain or Quantity array or NumPy ndarray
- The spike times.
- axis : int, optional
- The axis along which the difference is taken.
- Default is the last axis.
- Returns
- -------
- NumPy array or quantities array.
- """
- if axis is None:
- axis = -1
- if isinstance(spiketrain, neo.SpikeTrain):
- intervals = np.diff(
- np.sort(spiketrain.times.view(pq.Quantity)), axis=axis)
- else:
- intervals = np.diff(np.sort(spiketrain), axis=axis)
- return intervals
- def mean_firing_rate(spiketrain, t_start=None, t_stop=None, axis=None):
- """
- Return the firing rate of the SpikeTrain.
- Accepts a Neo SpikeTrain, a Quantity array, or a plain NumPy array.
- If either a SpikeTrain or Quantity array is provided, the return value will
- be a quantities array, otherwise a plain NumPy array. The units of
- the quantities array will be the inverse of the spiketrain.
- The interval over which the firing rate is calculated can be optionally
- controlled with `t_start` and `t_stop`
- Parameters
- ----------
- spiketrain : Neo SpikeTrain or Quantity array or NumPy ndarray
- The spike times.
- t_start : float or Quantity scalar, optional
- The start time to use for the interval.
- If not specified, retrieved from the``t_start`
- attribute of `spiketrain`. If that is not present, default to
- `0`. Any value from `spiketrain` below this value is ignored.
- t_stop : float or Quantity scalar, optional
- The stop time to use for the time points.
- If not specified, retrieved from the `t_stop`
- attribute of `spiketrain`. If that is not present, default to
- the maximum value of `spiketrain`. Any value from
- `spiketrain` above this value is ignored.
- axis : int, optional
- The axis over which to do the calculation.
- Default is `None`, do the calculation over the flattened array.
- Returns
- -------
- float, quantities scalar, NumPy array or quantities array.
- Notes
- -----
- If `spiketrain` is a Quantity or Neo SpikeTrain and `t_start` or `t_stop`
- are not, `t_start` and `t_stop` are assumed to have the same units as
- `spiketrain`.
- Raises
- ------
- TypeError
- If `spiketrain` is a NumPy array and `t_start` or `t_stop`
- is a quantity scalar.
- """
- if t_start is None:
- t_start = getattr(spiketrain, 't_start', 0)
- found_t_start = False
- if t_stop is None:
- if hasattr(spiketrain, 't_stop'):
- t_stop = spiketrain.t_stop
- else:
- t_stop = np.max(spiketrain, axis=axis)
- found_t_start = True
- # figure out what units, if any, we are dealing with
- if hasattr(spiketrain, 'units'):
- units = spiketrain.units
- else:
- units = None
- # convert everything to the same units
- if hasattr(t_start, 'units'):
- if units is None:
- raise TypeError('t_start cannot be a Quantity if '
- 'spiketrain is not a quantity')
- t_start = t_start.rescale(units)
- elif units is not None:
- t_start = pq.Quantity(t_start, units=units)
- if hasattr(t_stop, 'units'):
- if units is None:
- raise TypeError('t_stop cannot be a Quantity if '
- 'spiketrain is not a quantity')
- t_stop = t_stop.rescale(units)
- elif units is not None:
- t_stop = pq.Quantity(t_stop, units=units)
- if not axis or not found_t_start:
- return np.sum((spiketrain >= t_start) & (spiketrain <= t_stop),
- axis=axis) / (t_stop - t_start)
- else:
- # this is needed to handle broadcasting between spiketrain and t_stop
- t_stop_test = np.expand_dims(t_stop, axis)
- return np.sum((spiketrain >= t_start) & (spiketrain <= t_stop_test),
- axis=axis) / (t_stop - t_start)
- # we make `cv` an alias for scipy.stats.variation for the convenience
- # of former NeuroTools users
- cv = scipy.stats.variation
- def fanofactor(spiketrains):
- """
- Evaluates the empirical Fano factor F of the spike counts of
- a list of `neo.core.SpikeTrain` objects.
- Given the vector v containing the observed spike counts (one per
- spike train) in the time window [t0, t1], F is defined as:
- F := var(v)/mean(v).
- The Fano factor is typically computed for spike trains representing the
- activity of the same neuron over different trials. The higher F, the larger
- the cross-trial non-stationarity. In theory for a time-stationary Poisson
- process, F=1.
- Parameters
- ----------
- spiketrains : list of neo.SpikeTrain objects, quantity arrays, numpy arrays or lists
- Spike trains for which to compute the Fano factor of spike counts.
- Returns
- -------
- fano : float or nan
- The Fano factor of the spike counts of the input spike trains. If an
- empty list is specified, or if all spike trains are empty, F:=nan.
- """
- # Build array of spike counts (one per spike train)
- spike_counts = np.array([len(t) for t in spiketrains])
- # Compute FF
- if all([count == 0 for count in spike_counts]):
- fano = np.nan
- else:
- fano = spike_counts.var() / spike_counts.mean()
- return fano
- def lv(v):
- """
- Calculate the measure of local variation LV for
- a sequence of time intervals between events.
- Given a vector v containing a sequence of intervals, the LV is
- defined as:
- .math $$ LV := \\frac{1}{N}\\sum_{i=1}^{N-1}
- \\frac{3(isi_i-isi_{i+1})^2}
- {(isi_i+isi_{i+1})^2} $$
- The LV is typically computed as a substitute for the classical
- coefficient of variation for sequences of events which include
- some (relatively slow) rate fluctuation. As with the CV, LV=1 for
- a sequence of intervals generated by a Poisson process.
- Parameters
- ----------
- v : quantity array, numpy array or list
- Vector of consecutive time intervals
- Returns
- -------
- lvar : float
- The LV of the inter-spike interval of the input sequence.
- Raises
- ------
- AttributeError :
- If an empty list is specified, or if the sequence has less
- than two entries, an AttributeError will be raised.
- ValueError :
- Only vector inputs are supported. If a matrix is passed to the
- function a ValueError will be raised.
- References
- ----------
- ..[1] Shinomoto, S., Shima, K., & Tanji, J. (2003). Differences in spiking
- patterns among cortical neurons. Neural Computation, 15, 2823–2842.
- """
- # convert to array, cast to float
- v = np.asarray(v)
- # ensure we have enough entries
- if v.size < 2:
- raise AttributeError("Input size is too small. Please provide "
- "an input with more than 1 entry.")
- # calculate LV and return result
- # raise error if input is multi-dimensional
- return 3. * np.mean(np.power(np.diff(v) / (v[:-1] + v[1:]), 2))
- # sigma2kw and kw2sigma only needed for oldfct_instantaneous_rate!
- # to finally be taken out of Elephant
- def sigma2kw(form):
- warnings.simplefilter('always', DeprecationWarning)
- warnings.warn("deprecated", DeprecationWarning, stacklevel=2)
- if form.upper() == 'BOX':
- coeff = 2.0 * np.sqrt(3)
- elif form.upper() == 'TRI':
- coeff = 2.0 * np.sqrt(6)
- elif form.upper() == 'EPA':
- coeff = 2.0 * np.sqrt(5)
- elif form.upper() == 'GAU':
- coeff = 2.0 * 2.7 # > 99% of distribution weight
- elif form.upper() == 'ALP':
- coeff = 5.0
- elif form.upper() == 'EXP':
- coeff = 5.0
- return coeff
- def kw2sigma(form):
- warnings.simplefilter('always', DeprecationWarning)
- warnings.warn("deprecated", DeprecationWarning, stacklevel=2)
- return 1/sigma2kw(form)
- # to finally be taken out of Elephant
- def make_kernel(form, sigma, sampling_period, direction=1):
- """
- Creates kernel functions for convolution.
- Constructs a numeric linear convolution kernel of basic shape to be used
- for data smoothing (linear low pass filtering) and firing rate estimation
- from single trial or trial-averaged spike trains.
- Exponential and alpha kernels may also be used to represent postynaptic
- currents / potentials in a linear (current-based) model.
- Parameters
- ----------
- form : {'BOX', 'TRI', 'GAU', 'EPA', 'EXP', 'ALP'}
- Kernel form. Currently implemented forms are BOX (boxcar),
- TRI (triangle), GAU (gaussian), EPA (epanechnikov), EXP (exponential),
- ALP (alpha function). EXP and ALP are asymmetric kernel forms and
- assume optional parameter `direction`.
- sigma : Quantity
- Standard deviation of the distribution associated with kernel shape.
- This parameter defines the time resolution of the kernel estimate
- and makes different kernels comparable (cf. [1] for symmetric kernels).
- This is used here as an alternative definition to the cut-off
- frequency of the associated linear filter.
- sampling_period : float
- Temporal resolution of input and output.
- direction : {-1, 1}
- Asymmetric kernels have two possible directions.
- The values are -1 or 1, default is 1. The
- definition here is that for direction = 1 the
- kernel represents the impulse response function
- of the linear filter. Default value is 1.
- Returns
- -------
- kernel : numpy.ndarray
- Array of kernel. The length of this array is always an odd
- number to represent symmetric kernels such that the center bin
- coincides with the median of the numeric array, i.e for a
- triangle, the maximum will be at the center bin with equal
- number of bins to the right and to the left.
- norm : float
- For rate estimates. The kernel vector is normalized such that
- the sum of all entries equals unity sum(kernel)=1. When
- estimating rate functions from discrete spike data (0/1) the
- additional parameter `norm` allows for the normalization to
- rate in spikes per second.
- For example:
- ``rate = norm * scipy.signal.lfilter(kernel, 1, spike_data)``
- m_idx : int
- Index of the numerically determined median (center of gravity)
- of the kernel function.
- Examples
- --------
- To obtain single trial rate function of trial one should use::
- r = norm * scipy.signal.fftconvolve(sua, kernel)
- To obtain trial-averaged spike train one should use::
- r_avg = norm * scipy.signal.fftconvolve(sua, np.mean(X,1))
- where `X` is an array of shape `(l,n)`, `n` is the number of trials and
- `l` is the length of each trial.
- See also
- --------
- elephant.statistics.instantaneous_rate
- References
- ----------
- .. [1] Meier R, Egert U, Aertsen A, Nawrot MP, "FIND - a unified framework
- for neural data analysis"; Neural Netw. 2008 Oct; 21(8):1085-93.
- .. [2] Nawrot M, Aertsen A, Rotter S, "Single-trial estimation of neuronal
- firing rates - from single neuron spike trains to population activity";
- J. Neurosci Meth 94: 81-92; 1999.
- """
- warnings.simplefilter('always', DeprecationWarning)
- warnings.warn("deprecated", DeprecationWarning, stacklevel=2)
- forms_abbreviated = np.array(['BOX', 'TRI', 'GAU', 'EPA', 'EXP', 'ALP'])
- forms_verbose = np.array(['boxcar', 'triangle', 'gaussian', 'epanechnikov',
- 'exponential', 'alpha'])
- if form in forms_verbose:
- form = forms_abbreviated[forms_verbose == form][0]
- assert form.upper() in ('BOX', 'TRI', 'GAU', 'EPA', 'EXP', 'ALP'), \
- "form must be one of either 'BOX','TRI','GAU','EPA','EXP' or 'ALP'!"
- assert direction in (1, -1), "direction must be either 1 or -1"
- # conversion to SI units (s)
- if sigma < 0:
- raise ValueError('sigma must be positive!')
- SI_sigma = sigma.rescale('s').magnitude
- SI_time_stamp_resolution = sampling_period.rescale('s').magnitude
- norm = 1. / SI_time_stamp_resolution
- if form.upper() == 'BOX':
- w = 2.0 * SI_sigma * np.sqrt(3)
- # always odd number of bins
- width = 2 * np.floor(w / 2.0 / SI_time_stamp_resolution) + 1
- height = 1. / width
- kernel = np.ones((1, width)) * height # area = 1
- elif form.upper() == 'TRI':
- w = 2 * SI_sigma * np.sqrt(6)
- halfwidth = np.floor(w / 2.0 / SI_time_stamp_resolution)
- trileft = np.arange(1, halfwidth + 2)
- triright = np.arange(halfwidth, 0, -1) # odd number of bins
- triangle = np.append(trileft, triright)
- kernel = triangle / triangle.sum() # area = 1
- elif form.upper() == 'EPA':
- w = 2.0 * SI_sigma * np.sqrt(5)
- halfwidth = np.floor(w / 2.0 / SI_time_stamp_resolution)
- base = np.arange(-halfwidth, halfwidth + 1)
- parabula = base**2
- epanech = parabula.max() - parabula # inverse parabula
- kernel = epanech / epanech.sum() # area = 1
- elif form.upper() == 'GAU':
- w = 2.0 * SI_sigma * 2.7 # > 99% of distribution weight
- halfwidth = np.floor(w / 2.0 / SI_time_stamp_resolution) # always odd
- base = np.arange(-halfwidth, halfwidth + 1) * SI_time_stamp_resolution
- g = np.exp(
- -(base**2) / 2.0 / SI_sigma**2) / SI_sigma / np.sqrt(2.0 * np.pi)
- kernel = g / g.sum() # area = 1
- elif form.upper() == 'ALP':
- w = 5.0 * SI_sigma
- alpha = np.arange(
- 1, (
- 2.0 * np.floor(w / SI_time_stamp_resolution / 2.0) + 1) +
- 1) * SI_time_stamp_resolution
- alpha = (2.0 / SI_sigma**2) * alpha * np.exp(
- -alpha * np.sqrt(2) / SI_sigma)
- kernel = alpha / alpha.sum() # normalization
- if direction == -1:
- kernel = np.flipud(kernel)
- elif form.upper() == 'EXP':
- w = 5.0 * SI_sigma
- expo = np.arange(
- 1, (
- 2.0 * np.floor(w / SI_time_stamp_resolution / 2.0) + 1) +
- 1) * SI_time_stamp_resolution
- expo = np.exp(-expo / SI_sigma)
- kernel = expo / expo.sum()
- if direction == -1:
- kernel = np.flipud(kernel)
- kernel = kernel.ravel()
- m_idx = np.nonzero(kernel.cumsum() >= 0.5)[0].min()
- return kernel, norm, m_idx
- # to finally be taken out of Elephant
- def oldfct_instantaneous_rate(spiketrain, sampling_period, form,
- sigma='auto', t_start=None, t_stop=None,
- acausal=True, trim=False):
- """
- Estimate instantaneous firing rate by kernel convolution.
- Parameters
- -----------
- spiketrain: 'neo.SpikeTrain'
- Neo object that contains spike times, the unit of the time stamps
- and t_start and t_stop of the spike train.
- sampling_period : Quantity
- time stamp resolution of the spike times. the same resolution will
- be assumed for the kernel
- form : {'BOX', 'TRI', 'GAU', 'EPA', 'EXP', 'ALP'}
- Kernel form. Currently implemented forms are BOX (boxcar),
- TRI (triangle), GAU (gaussian), EPA (epanechnikov), EXP (exponential),
- ALP (alpha function). EXP and ALP are asymmetric kernel forms and
- assume optional parameter `direction`.
- sigma : string or Quantity
- Standard deviation of the distribution associated with kernel shape.
- This parameter defines the time resolution of the kernel estimate
- and makes different kernels comparable (cf. [1] for symmetric kernels).
- This is used here as an alternative definition to the cut-off
- frequency of the associated linear filter.
- Default value is 'auto'. In this case, the optimized kernel width for
- the rate estimation is calculated according to [1]. Note that the
- automatized calculation of the kernel width ONLY works for gaussian
- kernel shapes!
- t_start : Quantity (Optional)
- start time of the interval used to compute the firing rate, if None
- assumed equal to spiketrain.t_start
- Default:None
- t_stop : Qunatity
- End time of the interval used to compute the firing rate (included).
- If none assumed equal to spiketrain.t_stop
- Default:None
- acausal : bool
- if True, acausal filtering is used, i.e., the gravity center of the
- filter function is aligned with the spike to convolve
- Default:None
- m_idx : int
- index of the value in the kernel function vector that corresponds
- to its gravity center. this parameter is not mandatory for
- symmetrical kernels but it is required when asymmetrical kernels
- are to be aligned at their gravity center with the event times if None
- is assumed to be the median value of the kernel support
- Default : None
- trim : bool
- if True, only the 'valid' region of the convolved
- signal are returned, i.e., the points where there
- isn't complete overlap between kernel and spike train
- are discarded
- NOTE: if True and an asymmetrical kernel is provided
- the output will not be aligned with [t_start, t_stop]
- Returns
- -------
- rate : neo.AnalogSignal
- Contains the rate estimation in unit hertz (Hz).
- Has a property 'rate.times' which contains the time axis of the rate
- estimate. The unit of this property is the same as the resolution that
- is given as an argument to the function.
- Raises
- ------
- TypeError:
- If argument value for the parameter `sigma` is not a quantity object
- or string 'auto'.
- See also
- --------
- elephant.statistics.make_kernel
- References
- ----------
- ..[1] H. Shimazaki, S. Shinomoto, J Comput Neurosci (2010) 29:171–182.
- """
- warnings.simplefilter('always', DeprecationWarning)
- warnings.warn("deprecated", DeprecationWarning, stacklevel=2)
- if sigma == 'auto':
- form = 'GAU'
- unit = spiketrain.units
- kernel_width = sskernel(spiketrain.magnitude, tin=None,
- bootstrap=True)['optw']
- sigma = kw2sigma(form) * kernel_width * unit
- elif not isinstance(sigma, pq.Quantity):
- raise TypeError('sigma must be either a quantities object or "auto".'
- ' Found: %s, value %s' % (type(sigma), str(sigma)))
- kernel, norm, m_idx = make_kernel(form=form, sigma=sigma,
- sampling_period=sampling_period)
- units = pq.CompoundUnit(
- "%s*s" % str(sampling_period.rescale('s').magnitude))
- spiketrain = spiketrain.rescale(units)
- if t_start is None:
- t_start = spiketrain.t_start
- else:
- t_start = t_start.rescale(spiketrain.units)
- if t_stop is None:
- t_stop = spiketrain.t_stop
- else:
- t_stop = t_stop.rescale(spiketrain.units)
- time_vector = np.zeros(int((t_stop - t_start)) + 1)
- spikes_slice = spiketrain.time_slice(t_start, t_stop) \
- if len(spiketrain) else np.array([])
- for spike in spikes_slice:
- index = int((spike - t_start))
- time_vector[index] += 1
- r = norm * scipy.signal.fftconvolve(time_vector, kernel, 'full')
- if np.any(r < 0):
- warnings.warn('Instantaneous firing rate approximation contains '
- 'negative values, possibly caused due to machine '
- 'precision errors')
- if acausal:
- if not trim:
- r = r[m_idx:-(kernel.size - m_idx)]
- elif trim:
- r = r[2 * m_idx:-2 * (kernel.size - m_idx)]
- t_start = t_start + m_idx * spiketrain.units
- t_stop = t_stop - ((kernel.size) - m_idx) * spiketrain.units
- else:
- if not trim:
- r = r[m_idx:-(kernel.size - m_idx)]
- elif trim:
- r = r[2 * m_idx:-2 * (kernel.size - m_idx)]
- t_start = t_start + m_idx * spiketrain.units
- t_stop = t_stop - ((kernel.size) - m_idx) * spiketrain.units
- rate = neo.AnalogSignal(signal=r.reshape(r.size, 1),
- sampling_period=sampling_period,
- units=pq.Hz, t_start=t_start)
- return rate, sigma
- def instantaneous_rate(spiketrain, sampling_period, kernel='auto',
- cutoff=5.0, t_start=None, t_stop=None, trim=False):
- """
- Estimates instantaneous firing rate by kernel convolution.
- Parameters
- -----------
- spiketrain : 'neo.SpikeTrain'
- Neo object that contains spike times, the unit of the time stamps
- and t_start and t_stop of the spike train.
- sampling_period : Time Quantity
- Time stamp resolution of the spike times. The same resolution will
- be assumed for the kernel
- kernel : string 'auto' or callable object of :class:`Kernel` from module
- 'kernels.py'. Currently implemented kernel forms are rectangular,
- triangular, epanechnikovlike, gaussian, laplacian, exponential,
- and alpha function.
- Example: kernel = kernels.RectangularKernel(sigma=10*ms, invert=False)
- The kernel is used for convolution with the spike train and its
- standard deviation determines the time resolution of the instantaneous
- rate estimation.
- Default: 'auto'. In this case, the optimized kernel width for the
- rate estimation is calculated according to [1] and with this width
- a gaussian kernel is constructed. Automatized calculation of the
- kernel width is not available for other than gaussian kernel shapes.
- cutoff : float
- This factor determines the cutoff of the probability distribution of
- the kernel, i.e., the considered width of the kernel in terms of
- multiples of the standard deviation sigma.
- Default: 5.0
- t_start : Time Quantity (optional)
- Start time of the interval used to compute the firing rate. If None
- assumed equal to spiketrain.t_start
- Default: None
- t_stop : Time Quantity (optional)
- End time of the interval used to compute the firing rate (included).
- If None assumed equal to spiketrain.t_stop
- Default: None
- trim : bool
- if False, the output of the Fast Fourier Transformation being a longer
- vector than the input vector by the size of the kernel is reduced back
- to the original size of the considered time interval of the spiketrain
- using the median of the kernel.
- if True, only the region of the convolved signal is returned, where
- there is complete overlap between kernel and spike train. This is
- achieved by reducing the length of the output of the Fast Fourier
- Transformation by a total of two times the size of the kernel, and
- t_start and t_stop are adjusted.
- Default: False
- Returns
- -------
- rate : neo.AnalogSignal
- Contains the rate estimation in unit hertz (Hz).
- Has a property 'rate.times' which contains the time axis of the rate
- estimate. The unit of this property is the same as the resolution that
- is given via the argument 'sampling_period' to the function.
- Raises
- ------
- TypeError:
- If `spiketrain` is not an instance of :class:`SpikeTrain` of Neo.
- If `sampling_period` is not a time quantity.
- If `kernel` is neither instance of :class:`Kernel` or string 'auto'.
- If `cutoff` is neither float nor int.
- If `t_start` and `t_stop` are neither None nor a time quantity.
- If `trim` is not bool.
- ValueError:
- If `sampling_period` is smaller than zero.
- Example
- --------
- kernel = kernels.AlphaKernel(sigma = 0.05*s, invert = True)
- rate = instantaneous_rate(spiketrain, sampling_period = 2*ms, kernel)
- References
- ----------
- ..[1] H. Shimazaki, S. Shinomoto, J Comput Neurosci (2010) 29:171–182.
- """
- # Checks of input variables:
- if not isinstance(spiketrain, SpikeTrain):
- raise TypeError(
- "spiketrain must be instance of :class:`SpikeTrain` of Neo!\n"
- " Found: %s, value %s" % (type(spiketrain), str(spiketrain)))
- if not (isinstance(sampling_period, pq.Quantity) and
- sampling_period.dimensionality.simplified ==
- pq.Quantity(1, "s").dimensionality):
- raise TypeError(
- "The sampling period must be a time quantity!\n"
- " Found: %s, value %s" % (type(sampling_period), str(sampling_period)))
- if sampling_period.magnitude < 0:
- raise ValueError("The sampling period must be larger than zero.")
- if kernel == 'auto':
- kernel_width = sskernel(spiketrain.magnitude, tin=None,
- bootstrap=True)['optw']
- unit = spiketrain.units
- sigma = 1/(2.0 * 2.7) * kernel_width * unit
- # factor 2.0 connects kernel width with its half width,
- # factor 2.7 connects half width of Gaussian distribution with
- # 99% probability mass with its standard deviation.
- kernel = kernels.GaussianKernel(sigma)
- elif not isinstance(kernel, kernels.Kernel):
- raise TypeError(
- "kernel must be either instance of :class:`Kernel` "
- "or the string 'auto'!\n"
- " Found: %s, value %s" % (type(kernel), str(kernel)))
- if not (isinstance(cutoff, float) or isinstance(cutoff, int)):
- raise TypeError("cutoff must be float or integer!")
- if not (t_start is None or (isinstance(t_start, pq.Quantity) and
- t_start.dimensionality.simplified ==
- pq.Quantity(1, "s").dimensionality)):
- raise TypeError("t_start must be a time quantity!")
- if not (t_stop is None or (isinstance(t_stop, pq.Quantity) and
- t_stop.dimensionality.simplified ==
- pq.Quantity(1, "s").dimensionality)):
- raise TypeError("t_stop must be a time quantity!")
- if not (isinstance(trim, bool)):
- raise TypeError("trim must be bool!")
- # main function:
- units = pq.CompoundUnit("%s*s" % str(sampling_period.rescale('s').magnitude))
- spiketrain = spiketrain.rescale(units)
- if t_start is None:
- t_start = spiketrain.t_start
- else:
- t_start = t_start.rescale(spiketrain.units)
- if t_stop is None:
- t_stop = spiketrain.t_stop
- else:
- t_stop = t_stop.rescale(spiketrain.units)
- time_vector = np.zeros(int((t_stop - t_start)) + 1)
- spikes_slice = spiketrain.time_slice(t_start, t_stop) \
- if len(spiketrain) else np.array([])
- for spike in spikes_slice:
- index = int((spike - t_start))
- time_vector[index] += 1
- if cutoff < kernel.min_cutoff:
- cutoff = kernel.min_cutoff
- warnings.warn("The width of the kernel was adjusted to a minimally "
- "allowed width.")
- t_arr = np.arange(-cutoff * kernel.sigma.rescale(units).magnitude,
- cutoff * kernel.sigma.rescale(units).magnitude +
- sampling_period.rescale(units).magnitude,
- sampling_period.rescale(units).magnitude) * units
- r = scipy.signal.fftconvolve(time_vector,
- kernel(t_arr).rescale(pq.Hz).magnitude, 'full')
- if np.any(r < 0):
- warnings.warn("Instantaneous firing rate approximation contains "
- "negative values, possibly caused due to machine "
- "precision errors.")
- if not trim:
- r = r[kernel.median_index(t_arr):-(kernel(t_arr).size -
- kernel.median_index(t_arr))]
- elif trim:
- r = r[2 * kernel.median_index(t_arr):-2 * (kernel(t_arr).size -
- kernel.median_index(t_arr))]
- t_start += kernel.median_index(t_arr) * spiketrain.units
- t_stop -= (kernel(t_arr).size -
- kernel.median_index(t_arr)) * spiketrain.units
- rate = neo.AnalogSignal(signal=r.reshape(r.size, 1),
- sampling_period=sampling_period,
- units=pq.Hz, t_start=t_start, t_stop=t_stop)
- return rate
- def time_histogram(spiketrains, binsize, t_start=None, t_stop=None,
- output='counts', binary=False):
- """
- Time Histogram of a list of :attr:`neo.SpikeTrain` objects.
- Parameters
- ----------
- spiketrains : List of neo.SpikeTrain objects
- Spiketrains with a common time axis (same `t_start` and `t_stop`)
- binsize : quantities.Quantity
- Width of the histogram's time bins.
- t_start, t_stop : Quantity (optional)
- Start and stop time of the histogram. Only events in the input
- `spiketrains` falling between `t_start` and `t_stop` (both included)
- are considered in the histogram. If `t_start` and/or `t_stop` are not
- specified, the maximum `t_start` of all :attr:spiketrains is used as
- `t_start`, and the minimum `t_stop` is used as `t_stop`.
- Default: t_start = t_stop = None
- output : str (optional)
- Normalization of the histogram. Can be one of:
- * `counts`'`: spike counts at each bin (as integer numbers)
- * `mean`: mean spike counts per spike train
- * `rate`: mean spike rate per spike train. Like 'mean', but the
- counts are additionally normalized by the bin width.
- binary : bool (optional)
- If **True**, indicates whether all spiketrain objects should first
- binned to a binary representation (using the `BinnedSpikeTrain` class
- in the `conversion` module) and the calculation of the histogram is
- based on this representation.
- Note that the output is not binary, but a histogram of the converted,
- binary representation.
- Default: False
- Returns
- -------
- time_hist : neo.AnalogSignal
- A neo.AnalogSignal object containing the histogram values.
- `AnalogSignal[j]` is the histogram computed between
- `t_start + j * binsize` and `t_start + (j + 1) * binsize`.
- See also
- --------
- elephant.conversion.BinnedSpikeTrain
- """
- min_tstop = 0
- if t_start is None:
- # Find the internal range for t_start, where all spike trains are
- # defined; cut all spike trains taking that time range only
- max_tstart, min_tstop = conv._get_start_stop_from_input(spiketrains)
- t_start = max_tstart
- if not all([max_tstart == t.t_start for t in spiketrains]):
- warnings.warn(
- "Spiketrains have different t_start values -- "
- "using maximum t_start as t_start.")
- if t_stop is None:
- # Find the internal range for t_stop
- if min_tstop:
- t_stop = min_tstop
- if not all([min_tstop == t.t_stop for t in spiketrains]):
- warnings.warn(
- "Spiketrains have different t_stop values -- "
- "using minimum t_stop as t_stop.")
- else:
- min_tstop = conv._get_start_stop_from_input(spiketrains)[1]
- t_stop = min_tstop
- if not all([min_tstop == t.t_stop for t in spiketrains]):
- warnings.warn(
- "Spiketrains have different t_stop values -- "
- "using minimum t_stop as t_stop.")
- sts_cut = [st.time_slice(t_start=t_start, t_stop=t_stop) for st in
- spiketrains]
- # Bin the spike trains and sum across columns
- bs = conv.BinnedSpikeTrain(sts_cut, t_start=t_start, t_stop=t_stop,
- binsize=binsize)
- if binary:
- bin_hist = bs.to_sparse_bool_array().sum(axis=0)
- else:
- bin_hist = bs.to_sparse_array().sum(axis=0)
- # Flatten array
- bin_hist = np.ravel(bin_hist)
- # Renormalise the histogram
- if output == 'counts':
- # Raw
- bin_hist = bin_hist * pq.dimensionless
- elif output == 'mean':
- # Divide by number of input spike trains
- bin_hist = bin_hist * 1. / len(spiketrains) * pq.dimensionless
- elif output == 'rate':
- # Divide by number of input spike trains and bin width
- bin_hist = bin_hist * 1. / len(spiketrains) / binsize
- else:
- raise ValueError('Parameter output is not valid.')
- return neo.AnalogSignal(signal=bin_hist.reshape(bin_hist.size, 1),
- sampling_period=binsize, units=bin_hist.units,
- t_start=t_start)
- def complexity_pdf(spiketrains, binsize):
- """
- Complexity Distribution [1] of a list of :attr:`neo.SpikeTrain` objects.
- Probability density computed from the complexity histogram which is the
- histogram of the entries of the population histogram of clipped (binary)
- spike trains computed with a bin width of binsize.
- It provides for each complexity (== number of active neurons per bin) the
- number of occurrences. The normalization of that histogram to 1 is the
- probability density.
- Parameters
- ----------
- spiketrains : List of neo.SpikeTrain objects
- Spiketrains with a common time axis (same `t_start` and `t_stop`)
- binsize : quantities.Quantity
- Width of the histogram's time bins.
- Returns
- -------
- time_hist : neo.AnalogSignal
- A neo.AnalogSignal object containing the histogram values.
- `AnalogSignal[j]` is the histogram computed between .
- See also
- --------
- elephant.conversion.BinnedSpikeTrain
- References
- ----------
- [1]Gruen, S., Abeles, M., & Diesmann, M. (2008). Impact of higher-order
- correlations on coincidence distributions of massively parallel data.
- In Dynamic Brain-from Neural Spikes to Behaviors (pp. 96-114).
- Springer Berlin Heidelberg.
- """
- # Computing the population histogram with parameter binary=True to clip the
- # spike trains before summing
- pophist = time_histogram(spiketrains, binsize, binary=True)
- # Computing the histogram of the entries of pophist (=Complexity histogram)
- complexity_hist = np.histogram(
- pophist.magnitude, bins=range(0, len(spiketrains) + 2))[0]
- # Normalization of the Complexity Histogram to 1 (probabilty distribution)
- complexity_hist = complexity_hist / complexity_hist.sum()
- # Convert the Complexity pdf to an neo.AnalogSignal
- complexity_distribution = neo.AnalogSignal(
- np.array(complexity_hist).reshape(len(complexity_hist), 1) *
- pq.dimensionless, t_start=0 * pq.dimensionless,
- sampling_period=1 * pq.dimensionless)
- return complexity_distribution
- """Kernel Bandwidth Optimization.
- Python implementation by Subhasis Ray.
- Original matlab code (sskernel.m) here:
- http://2000.jukuin.keio.ac.jp/shimazaki/res/kernel.html
- This was translated into Python by Subhasis Ray, NCBS. Tue Jun 10
- 23:01:43 IST 2014
- """
- def nextpow2(x):
- """ Return the smallest integral power of 2 that >= x """
- n = 2
- while n < x:
- n = 2 * n
- return n
- def fftkernel(x, w):
- """
- y = fftkernel(x,w)
- Function `fftkernel' applies the Gauss kernel smoother to an input
- signal using FFT algorithm.
- Input argument
- x: Sample signal vector.
- w: Kernel bandwidth (the standard deviation) in unit of
- the sampling resolution of x.
- Output argument
- y: Smoothed signal.
- MAY 5/23, 2012 Author Hideaki Shimazaki
- RIKEN Brain Science Insitute
- http://2000.jukuin.keio.ac.jp/shimazaki
- Ported to Python: Subhasis Ray, NCBS. Tue Jun 10 10:42:38 IST 2014
- """
- L = len(x)
- Lmax = L + 3 * w
- n = nextpow2(Lmax)
- X = np.fft.fft(x, n)
- f = np.arange(0, n, 1.0) / n
- f = np.concatenate((-f[:int(n / 2)], f[int(n / 2):0:-1]))
- K = np.exp(-0.5 * (w * 2 * np.pi * f)**2)
- y = np.fft.ifft(X * K, n)
- y = y[:L].copy()
- return y
- def logexp(x):
- if x < 1e2:
- y = np.log(1 + np.exp(x))
- else:
- y = x
- return y
- def ilogexp(x):
- if x < 1e2:
- y = np.log(np.exp(x) - 1)
- else:
- y = x
- return y
- def cost_function(x, N, w, dt):
- """
- The cost function
- Cn(w) = sum_{i,j} int k(x - x_i) k(x - x_j) dx - 2 sum_{i~=j} k(x_i - x_j)
- """
- yh = np.abs(fftkernel(x, w / dt)) # density
- # formula for density
- C = np.sum(yh ** 2) * dt - 2 * np.sum(yh * x) * \
- dt + 2 / np.sqrt(2 * np.pi) / w / N
- C = C * N * N
- # formula for rate
- # C = dt*sum( yh.^2 - 2*yh.*y_hist + 2/sqrt(2*pi)/w*y_hist )
- return C, yh
- def sskernel(spiketimes, tin=None, w=None, bootstrap=False):
- """
- Calculates optimal fixed kernel bandwidth.
- spiketimes: sequence of spike times (sorted to be ascending).
- tin: (optional) time points at which the kernel bandwidth is to be estimated.
- w: (optional) vector of kernel bandwidths. If specified, optimal
- bandwidth is selected from this.
- bootstrap (optional): whether to calculate the 95% confidence
- interval. (default False)
- Returns
- A dictionary containing the following key value pairs:
- 'y': estimated density,
- 't': points at which estimation was computed,
- 'optw': optimal kernel bandwidth,
- 'w': kernel bandwidths examined,
- 'C': cost functions of w,
- 'confb95': (lower bootstrap confidence level, upper bootstrap confidence level),
- 'yb': bootstrap samples.
- Ref: Shimazaki, Hideaki, and Shigeru Shinomoto. 2010. Kernel
- Bandwidth Optimization in Spike Rate Estimation. Journal of
- Computational Neuroscience 29 (1-2):
- 171-82. doi:10.1007/s10827-009-0180-4.
- """
- if tin is None:
- time = np.max(spiketimes) - np.min(spiketimes)
- isi = np.diff(spiketimes)
- isi = isi[isi > 0].copy()
- dt = np.min(isi)
- tin = np.linspace(np.min(spiketimes),
- np.max(spiketimes),
- min(int(time / dt + 0.5), 1000)) # The 1000 seems somewhat arbitrary
- t = tin
- else:
- time = np.max(tin) - np.min(tin)
- spiketimes = spiketimes[(spiketimes >= np.min(tin)) &
- (spiketimes <= np.max(tin))].copy()
- isi = np.diff(spiketimes)
- isi = isi[isi > 0].copy()
- dt = np.min(isi)
- if dt > np.min(np.diff(tin)):
- t = np.linspace(np.min(tin), np.max(tin),
- min(int(time / dt + 0.5), 1000))
- else:
- t = tin
- dt = np.min(np.diff(tin))
- yhist, bins = np.histogram(spiketimes, np.r_[t - dt / 2, t[-1] + dt / 2])
- N = np.sum(yhist)
- yhist = yhist / (N * dt) # density
- optw = None
- y = None
- if w is not None:
- C = np.zeros(len(w))
- Cmin = np.inf
- for k, w_ in enumerate(w):
- C[k], yh = cost_function(yhist, N, w_, dt)
- if C[k] < Cmin:
- Cmin = C[k]
- optw = w_
- y = yh
- else:
- # Golden section search on a log-exp scale
- wmin = 2 * dt
- wmax = max(spiketimes) - min(spiketimes)
- imax = 20 # max iterations
- w = np.zeros(imax)
- C = np.zeros(imax)
- tolerance = 1e-5
- phi = 0.5 * (np.sqrt(5) + 1) # The Golden ratio
- a = ilogexp(wmin)
- b = ilogexp(wmax)
- c1 = (phi - 1) * a + (2 - phi) * b
- c2 = (2 - phi) * a + (phi - 1) * b
- f1, y1 = cost_function(yhist, N, logexp(c1), dt)
- f2, y2 = cost_function(yhist, N, logexp(c2), dt)
- k = 0
- while (np.abs(b - a) > (tolerance * (np.abs(c1) + np.abs(c2))))\
- and (k < imax):
- if f1 < f2:
- b = c2
- c2 = c1
- c1 = (phi - 1) * a + (2 - phi) * b
- f2 = f1
- f1, y1 = cost_function(yhist, N, logexp(c1), dt)
- w[k] = logexp(c1)
- C[k] = f1
- optw = logexp(c1)
- y = y1 / (np.sum(y1 * dt))
- else:
- a = c1
- c1 = c2
- c2 = (2 - phi) * a + (phi - 1) * b
- f1 = f2
- f2, y2 = cost_function(yhist, N, logexp(c2), dt)
- w[k] = logexp(c2)
- C[k] = f2
- optw = logexp(c2)
- y = y2 / np.sum(y2 * dt)
- k = k + 1
- # Bootstrap confidence intervals
- confb95 = None
- yb = None
- if bootstrap:
- nbs = 1000
- yb = np.zeros((nbs, len(tin)))
- for ii in range(nbs):
- idx = np.floor(np.random.rand(N) * N).astype(int)
- xb = spiketimes[idx]
- y_histb, bins = np.histogram(
- xb, np.r_[t - dt / 2, t[-1] + dt / 2]) / dt / N
- yb_buf = fftkernel(y_histb, optw / dt).real
- yb_buf = yb_buf / np.sum(yb_buf * dt)
- yb[ii, :] = np.interp(tin, t, yb_buf)
- ybsort = np.sort(yb, axis=0)
- y95b = ybsort[np.floor(0.05 * nbs).astype(int), :]
- y95u = ybsort[np.floor(0.95 * nbs).astype(int), :]
- confb95 = (y95b, y95u)
- ret = np.interp(tin, t, y)
- return {'y': ret,
- 't': tin,
- 'optw': optw,
- 'w': w,
- 'C': C,
- 'confb95': confb95,
- 'yb': yb}
|