nsim package¶

Submodules¶

nsim.nsim module¶

Implements the core functionality of nsim.

Classes:¶

Model base class for different kinds of dynamical model ODEModel system of ordinary differential equations SDEModel system of stochastic differential equations DelaySDEModel system of delay stochastic differential equations

Simulation single simulation run of a model, with simulation results

MultipleSim set of simulations, distributed. RepeatedSim repeated simulations of the same model (to get statistics) ParameterSim multiple simulations of a model exploring parameter space NetworkSim simulate many instances of a model coupled in a network

RemoteTimeseries Local proxy representing a Timeseries on a remote engine DistTimeseries Timeseries with one axis distributed onto multiple engines

functions:¶

newmodel() Create a new model class dynamically at runtime newsim() Create a new simulation dynamically at runtime

class nsim.nsim.DelaySDEModel¶

Bases: nsim.nsim.Model

Model defined by a system of stochastic delay differential equations

class nsim.nsim.DistTimeseries(subts, axis=None, axislabels=None)¶

Bases: distob.arrays.DistArray

a Timeseries with one axis distributed onto multiple computing engines.

For example, a multi-channel timeseries could be distributed so that each channel is held on a different computer (for parallel computation)

Currently only a non-time axis can be distributed.

__array_prepare__(out_arr, context=None)¶: Fetch underlying data to user’s computer and apply ufunc locally. Only used as a fallback, for numpy versions < 1.10.0 which lack support for the __numpy_ufunc__ mechanism.

__getitem__(index)¶: Slice the distributed timeseries

abs()¶: Calculate the absolute value element-wise.

absolute()¶

Calculate the absolute value element-wise.

Returns:	absolute – Absolute value. For complex input (a + bj) gives sqrt(aa + b*2)
Return type:	Timeseries

classmethod add_analyses(source, vectorize=False)¶

Dynamically add new analysis methods to the DistTimeseries class. :param source: Can be a function, module or the filename of a python file.

If a filename or a module is given, then all functions defined inside not starting with _ will be added as methods.

Parameters:	vectorize (bool) – Whether to apply distob.vectorize() to the function before making it a method of DistTimeseries.

The only restriction on the functions is that they can accept a Timeseries as their first argument.

angle(deg=False)¶

Return the angle of a complex Timeseries

Parameters:	deg (bool, optional) – Return angle in degrees if True, radians if False (default).
Returns:	angle – The counterclockwise angle from the positive real axis on the complex plane, with dtype as numpy.float64.
Return type:	Timeseries

ap_entropy(obj, *args, **kwargs)¶

Apply ap_entropy() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for ap_entropy() detailed below:

Computer approximate entropy (ApEN) of series X, specified by M and R.

Suppose given time series is X = [x(1), x(2), ... , x(N)]. We first build embedding matrix Em, of dimension (N-M+1)-by-M, such that the i-th row of Em is x(i),x(i+1), ... , x(i+M-1). Hence, the embedding lag and dimension are 1 and M-1 respectively. Such a matrix can be built by calling pyeeg function as Em = embed_seq(X, 1, M). Then we build matrix Emp, whose only difference with Em is that the length of each embedding sequence is M + 1

Denote the i-th and j-th row of Em as Em[i] and Em[j]. Their k-th elements are Em[i][k] and Em[j][k] respectively. The distance between Em[i] and Em[j] is defined as 1) the maximum difference of their corresponding scalar components, thus, max(Em[i]-Em[j]), or 2) Euclidean distance. We say two 1-D vectors Em[i] and Em[j] match in tolerance R, if the distance between them is no greater than R, thus, max(Em[i]-Em[j]) <= R. Mostly, the value of R is defined as 20% - 30% of standard deviation of X.

Pick Em[i] as a template, for all j such that 0 < j < N - M + 1, we can check whether Em[j] matches with Em[i]. Denote the number of Em[j], which is in the range of Em[i], as k[i], which is the i-th element of the vector k. The probability that a random row in Em matches Em[i] is simga_1^{N-M+1} k[i] / (N - M + 1), thus sum(k)/ (N - M + 1), denoted as Cm[i].

We repeat the same process on Emp and obtained Cmp[i], but here 0<i<N-M since the length of each sequence in Emp is M + 1.

The probability that any two embedding sequences in Em match is then sum(Cm)/ (N - M +1 ). We define Phi_m = sum(log(Cm)) / (N - M + 1) and Phi_mp = sum(log(Cmp)) / (N - M ).

And the ApEn is defined as Phi_m - Phi_mp.

Please be aware that self-match is also counted in ApEn.

This function now runs very slow. We are still trying to speed it up.

Costa M, Goldberger AL, Peng CK, Multiscale entropy analysis of biological signals, Physical Review E, 71:021906, 2005

samp_entropy: sample entropy of a time series

Extremely slow implementation. Do NOT use if your dataset is not small.

autocorrelation(obj, *args, **kwargs)¶

Apply autocorrelation() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for autocorrelation() detailed below:

Returns the discrete, linear convolution of a time series with itself, optionally using unbiased normalization.

N.B. Autocorrelation estimates are necessarily inaccurate for longer lags, as there are less pairs of points to convolve separated by that lag. Therefore best to throw out the results except for shorter lags, e.g. keep lags from tau=0 up to one quarter of the total time series length.

Args:

normalized (boolean): If True, the time series will first be normalized

to a mean of 0 and variance of 1. This gives autocorrelation 1 at zero lag.

unbiased (boolean): If True, the result at each lag m will be scaled by

1/(N-m). This gives an unbiased estimation of the autocorrelation of a stationary process from a finite length sample.

Ref: S. J. Orfanidis (1996) “Optimum Signal Processing”, 2nd Ed.

bandpass(obj, *args, **kwargs)¶

Apply bandpass() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for bandpass() detailed below:

forward-backward butterworth band-pass filter

bin_power(obj, *args, **kwargs)¶

Apply bin_power() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for bin_power() detailed below:

Compute power in each frequency bin specified by Band from FFT result of

By default, X is a real signal.

A real signal can be synthesized, thus not real.

Band

list

boundary frequencies (in Hz) of bins. They can be unequal bins, e.g. [0.5,4,7,12,30] which are delta, theta, alpha and beta respectively. You can also use range() function of Python to generate equal bins and pass the generated list to this function.

Each element of Band is a physical frequency and shall not exceed the Nyquist frequency, i.e., half of sampling frequency.

X

list

a 1-D real time series.

Fs

integer

the sampling rate in physical frequency

Power

list

spectral power in each frequency bin.

Power_ratio

list

spectral power in each frequency bin normalized by total power in ALL frequency bins.

circmean(dts, axis=2)¶: Circular mean phase

circstd(dts, axis=2)¶: Circular standard deviation

crossing_indices(obj, *args, **kwargs)¶

Apply crossing_indices() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for crossing_indices() detailed below:

For a single variable timeseries, find the time indices each time the

value crosses c from above or below. Can optionally set a non-zero d to impose the condition that the value must wander at least d units away from c between crossings.

If the timeseries begins (or ends) exactly at c, then time zero (or the ending time) is also included as a crossing event, so that the boundaries of the first and last excursions are included.

Args:

ts: Timeseries (single variable)

c (float): Critical value at which to report crossings.

d (float): Optional min distance from c to be attained between crossings.

Returns:

array of indices

cwt(obj, *args, **kwargs)¶

Apply cwt() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for cwt() detailed below:

Continuous wavelet transform

Note the full results can use a huge amount of memory at 64-bit precision

Args:

ts: Timeseries of m variables, shape (n, m). Assumed constant timestep. freqs: list of frequencies (in Hz) to use for the tranform.

(default is 50 frequency bins logarithmic from 1Hz to 100Hz)

wavelet: the wavelet to use. may be complex. see scipy.signal.wavelets plot: whether to plot time-resolved power spectrum

Returns:

coefs: Continuous wavelet transform output array, shape (n,len(freqs),m)

cwt_distributed(obj, *args, **kwargs)¶

Apply cwt_distributed() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for cwt_distributed() detailed below:

Continuous wavelet transform using distributed computation.

(Currently just splits the data by channel. TODO split it further.) Note: this function requires an IPython cluster to be started first.

Args:

ts: Timeseries of m variables, shape (n, m). Assumed constant timestep. freqs: list of frequencies (in Hz) to use for the tranform.

(default is 50 frequency bins logarithmic from 1Hz to 100Hz)

wavelet: the wavelet to use. may be complex. see scipy.signal.wavelets plot: whether to plot time-resolved power spectrum

Returns:

coefs: Continuous wavelet transform output array, shape (n,len(freqs),m)

dfa(obj, *args, **kwargs)¶

Apply dfa() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for dfa() detailed below:

Compute Detrended Fluctuation Analysis from a time series X and length of

boxes L.

The first step to compute DFA is to integrate the signal. Let original series be X= [x(1), x(2), ..., x(N)].

The integrated signal Y = [y(1), y(2), ..., y(N)] is obtained as follows y(k) = sum_{i=1}^{k}{x(i)-Ave} where Ave is the mean of X.

The second step is to partition/slice/segment the integrated sequence Y into boxes. At least two boxes are needed for computing DFA. Box sizes are specified by the L argument of this function. By default, it is from 1/5 of signal length to one (x-5)-th of the signal length, where x is the nearest power of 2 from the length of the signal, i.e., 1/16, 1/32, 1/64, 1/128, ...

In each box, a linear least square fitting is employed on data in the box. Denote the series on fitted line as Yn. Its k-th elements, yn(k), corresponds to y(k).

For fitting in each box, there is a residue, the sum of squares of all offsets, difference between actual points and points on fitted line.

F(n) denotes the square root of average total residue in all boxes when box length is n, thus Total_Residue = sum_{k=1}^{N}{(y(k)-yn(k))} F(n) = sqrt(Total_Residue/N)

The computing to F(n) is carried out for every box length n. Therefore, a relationship between n and F(n) can be obtained. In general, F(n) increases when n increases.

Finally, the relationship between F(n) and n is analyzed. A least square fitting is performed between log(F(n)) and log(n). The slope of the fitting line is the DFA value, denoted as Alpha. To white noise, Alpha should be 0.5. Higher level of signal complexity is related to higher Alpha.

X:: 1-D Python list or numpy array a time series
Ave:: integer, optional The average value of the time series
L:: 1-D Python list of integers A list of box size, integers in ascending order

Alpha:: integer the result of DFA analysis, thus the slope of fitting line of log(F(n)) vs. log(n). where n is the

>>> import pyeeg
>>> from numpy.random import randn
>>> print pyeeg.dfa(randn(4096))
0.490035110345

Peng C-K, Havlin S, Stanley HE, Goldberger AL. Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series. _Chaos_ 1995;5:82-87

This value depends on the box sizes very much. When the input is a white noise, this value should be 0.5. But, some choices on box sizes can lead to the value lower or higher than 0.5, e.g. 0.38 or 0.58.

Based on many test, I set the box sizes from 1/5 of signal length to one (x-5)-th of the signal length, where x is the nearest power of 2 from the length of the signal, i.e., 1/16, 1/32, 1/64, 1/128, ...

You may generate a list of box sizes and pass in such a list as a parameter.

embed_seq(obj, *args, **kwargs)¶

Apply embed_seq() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for embed_seq() detailed below:

Build a set of embedding sequences from given time series X with lag Tau

and embedding dimension DE. Let X = [x(1), x(2), ... , x(N)], then for each i such that 1 < i < N - (D - 1) * Tau, we build an embedding sequence, Y(i) = [x(i), x(i + Tau), ... , x(i + (D - 1) * Tau)]. All embedding sequence are placed in a matrix Y.

X

list

a time series

Tau

integer

the lag or delay when building embedding sequence

D

integer

the embedding dimension

Y

2-D list

embedding matrix built

>>> import pyeeg
>>> a=range(0,9)
>>> pyeeg.embed_seq(a,1,4)
array([[ 0.,  1.,  2.,  3.],
       [ 1.,  2.,  3.,  4.],
       [ 2.,  3.,  4.,  5.],
       [ 3.,  4.,  5.,  6.],
       [ 4.,  5.,  6.,  7.],
       [ 5.,  6.,  7.,  8.]])
>>> pyeeg.embed_seq(a,2,3)
array([[ 0.,  2.,  4.],
       [ 1.,  3.,  5.],
       [ 2.,  4.,  6.],
       [ 3.,  5.,  7.],
       [ 4.,  6.,  8.]])
>>> pyeeg.embed_seq(a,4,1)
array([[ 0.],
       [ 1.],
       [ 2.],
       [ 3.],
       [ 4.],
       [ 5.],
       [ 6.],
       [ 7.],
       [ 8.]])

epochs(ts, variability=None, threshold=0.0, minlength=1.0, plot=True)¶

Identify “stationary” epochs within a time series, based on a continuous measure of variability. Epochs are defined to contain the points of minimal variability, and to extend as wide as possible with variability not exceeding the threshold.

Parameters:

Timeseries of m variables, shape (n, m). (ts) –
(optional) Timeseries of shape (n, m, q), giving q scalar (variability) – measures of the variability of timeseries ts near each point in time. (if None, we will use variability_fp()) Epochs require the mean of these to be below the threshold.
The maximum variability permitted in stationary epochs. (threshold) –
Shortest acceptable epoch length (in seconds) (minlength) –
bool Whether to display the output (plot) –

Returns: (variability, allchannels_epochs): variability: as above allchannels_epochs: (list of) list of tuples For each variable, a list of tuples (start, end) that give the starting and ending indices of stationary epochs. (epochs are inclusive of start point but not the end point)

epochs_distributed(ts, variability=None, threshold=0.0, minlength=1.0, plot=True)¶

Same as epochs(), but computes channels in parallel for speed.

(Note: This requires an IPython cluster to be started first,: e.g. on a workstation type ‘ipcluster start’)

Identify “stationary” epochs within a time series, based on a continuous measure of variability. Epochs are defined to contain the points of minimal variability, and to extend as wide as possible with variability not exceeding the threshold.

Parameters:

Timeseries of m variables, shape (n, m). (ts) –
(optional) Timeseries of shape (n, m, q), giving q scalar (variability) – measures of the variability of timeseries ts near each point in time. (if None, we will use variability_fp()) Epochs require the mean of these to be below the threshold.
The maximum variability permitted in stationary epochs. (threshold) –
Shortest acceptable epoch length (in seconds) (minlength) –
bool Whether to display the output (plot) –

Returns: (variability, allchannels_epochs): variability: as above allchannels_epochs: (list of) list of tuples For each variable, a list of tuples (start, end) that give the starting and ending indices of stationary epochs. (epochs are inclusive of start point but not the end point)

epochs_joint(ts, variability=None, threshold=0.0, minlength=1.0, proportion=0.75, plot=True)¶

Identify epochs within a multivariate time series where at least a certain proportion of channels are “stationary”, based on a previously computed variability measure.

(Note: This requires an IPython cluster to be started first,: e.g. on a workstation type ‘ipcluster start’)

Parameters:

Timeseries of m variables, shape (n, m). (ts) –
(optional) Timeseries of shape (n, m), giving a scalar (variability) – measure of the variability of timeseries ts near each point in time. (if None, we will use variability_fp())
The maximum variability permitted in stationary epochs. (threshold) –
Shortest acceptable epoch length (in seconds) (minlength) –
Require at least this fraction of channels to be "stationary" (proportion) –
bool Whether to display the output (plot) –

Returns: (variability, joint_epochs): joint_epochs: list of tuples A list of tuples (start, end) that give the starting and ending indices of time epochs that are stationary for at least proportion of channels. (epochs are inclusive of start point but not the end point)

expand_dims(axis)¶

Insert a new axis, at a given position in the array shape :param axis: Position (amongst axes) where new axis is to be inserted.

first_return_times(dts, c=None, d=0.0)¶

For an ensemble of time series, return the set of all time intervals between successive returns to value c for all instances in the ensemble. If c is not given, the default is the mean across all times and across all time series in the ensemble.

Parameters:	(DistTimeseries) (dts) – c (float) – Optional target value (default is the ensemble mean value) d (float) – Optional min distance from c to be attained between returns
Returns:	array of time intervals (Can take the mean of these to estimate the expected first return time for the whole ensemble)

fisher_info(obj, *args, **kwargs)¶

Apply fisher_info() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for fisher_info() detailed below:

Compute SVD Entropy from either two cases below:

a time series X, with lag tau and embedding dimension dE (default)

2. a list, W, of normalized singular values of a matrix (if W is provided, recommend to speed up.)

If W is None, the function will do as follows to prepare singular spectrum:

First, computer an embedding matrix from X, Tau and DE using pyeeg function embed_seq():

M = embed_seq(X, Tau, DE)

Second, use scipy.linalg function svd to decompose the embedding matrix M and obtain a list of singular values:

W = svd(M, compute_uv=0)

At last, normalize W:

W /= sum(W)

To speed up, it is recommended to compute W before calling this function because W may also be used by other functions whereas computing it here again will slow down.

hfd(obj, *args, **kwargs)¶

Apply hfd() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for hfd() detailed below:

Compute Hjorth Fractal Dimension of a time series X, kmax

is an HFD parameter

highpass(obj, *args, **kwargs)¶

Apply highpass() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for highpass() detailed below:

forward-backward butterworth high-pass filter

hilbert(obj, *args, **kwargs)¶

Apply hilbert() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for hilbert() detailed below:

Analytic signal, using the Hilbert transform

hilbert_amplitude(obj, *args, **kwargs)¶

Apply hilbert_amplitude() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for hilbert_amplitude() detailed below:

Amplitude of the analytic signal, using the Hilbert transform

hilbert_phase(obj, *args, **kwargs)¶

Apply hilbert_phase() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for hilbert_phase() detailed below:

Phase of the analytic signal, using the Hilbert transform

hjorth(obj, *args, **kwargs)¶

Apply hjorth() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for hjorth() detailed below:

Compute Hjorth mobility and complexity of a time series from either two

cases below:

X, the time series of type list (default)

D, a first order differential sequence of X (if D is provided, recommended to speed up)

In case 1, D is computed using Numpy’s Difference function.

To speed up, it is recommended to compute D before calling this function because D may also be used by other functions whereas computing it here again will slow down.

X

list

a time series

D

list

first order differential sequence of a time series

As indicated in return line

Hjorth mobility and complexity

hurst(obj, *args, **kwargs)¶

Apply hurst() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for hurst() detailed below:

Compute the Hurst exponent of X. If the output H=0.5,the behavior
of the time-series is similar to random walk. If H<0.5, the time-series cover less “distance” than a random walk, vice verse.

X

list

a time series

H

float

Hurst exponent

Author of this function is Xin Liu
>>> import pyeeg
>>> from numpy.random import randn
>>> a = randn(4096)
>>> pyeeg.hurst(a)
0.5057444

in_range(obj, *args, **kwargs)¶

Apply in_range() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for in_range() detailed below:

Determines whether one vector is the the range of another vector.

The two vectors should have equal length.

Template

list The template vector, one of two vectors being compared

Scroll

list The scroll vector, one of the two vectors being compared

Distance

float Two vectors match if their distance is less than D

The distance between two vectors can be defined as Euclidean distance according to some publications.

The two vector should of equal length

lowpass(obj, *args, **kwargs)¶

Apply lowpass() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for lowpass() detailed below:

forward-backward butterworth low-pass filter

mod2pi(obj, *args, **kwargs)¶

Apply mod2pi() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for mod2pi() detailed below:

For a timeseries where all variables represent phases (in radians),: return an equivalent timeseries where all values are in the range (-pi, pi]

notch(obj, *args, **kwargs)¶

Apply notch() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for notch() detailed below:

notch filter to remove remove a particular frequency: Adapted from code by Sturla Molden

periods(dts, phi=0.0)¶

For an ensemble of oscillators, return the set of periods lengths of all successive oscillations of all oscillators.

An individual oscillation is defined to start and end when the phase passes phi (by default zero) after completing a full cycle.

If the timeseries of an oscillator phase begins (or ends) exactly at phi, then the first (or last) oscillation will be included.

Parameters:	dts (DistTimeseries) – where dts.shape[1] is 1 (single output variable representing phase) and axis 2 ranges over multiple realizations of the oscillator. phi=0.0 – float A single oscillation starts and ends at phase phi (by default zero).

pfd(obj, *args, **kwargs)¶

Apply pfd() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for pfd() detailed below:

Compute Petrosian Fractal Dimension of a time series from either two

cases below:

X, the time series of type list (default)
D, the first order differential sequence of X (if D is provided, recommended to speed up)

In case 1, D is computed using Numpy’s difference function.

To speed up, it is recommended to compute D before calling this function because D may also be used by other functions whereas computing it here again will slow down.

phase_crossings(obj, *args, **kwargs)¶

Apply phase_crossings() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for phase_crossings() detailed below:

For a single variable timeseries representing the phase of an oscillator,

find the time indices each time the phase crosses angle phi, with the condition that the phase must visit phi+pi between crossings.

(Thus if noise causes the phase to wander back and forth across angle phi without the oscillator doing a full revolution, then this is recorded as a single crossing event, giving the index of the earliest arrival.)

If the timeseries begins (or ends) exactly at phi, then time zero (or the ending time) is also included as a crossing event, so that the boundaries of the first and last oscillations are included.

Arguments:

ts: Timeseries (single variable): The timeseries of an angle variable (radians)

phi (float): Critical phase angle (radians) at which to report crossings.

Returns:

array of indices

phase_histogram(dts, times, nbins=30, colormap=<matplotlib.colors.LinearSegmentedColormap object>)¶

Plot a polar histogram of a phase variable’s probability distribution :param dts: DistTimeseries with axis 2 ranging over separate instances of an

oscillator (time series values are assumed to represent an angle)

Parameters:	times (float or sequence of floats) – The target times at which to plot the distribution nbins (int) – number of histogram bins colormap –

plot(dts, title=None, points=None, show=True)¶

Plot a distributed timeseries :param dts (DistTimeseries): :param title (str, optional): :param points: Limit the number of time points plotted.

If specified, will downsample to use this total number of time points, and only fetch back the necessary points to the client for plotting.

Returns:	fig

psd(obj, *args, **kwargs)¶

Apply psd() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for psd() detailed below:

plot Welch estimate of power spectral density

samp_entropy(obj, *args, **kwargs)¶

Apply samp_entropy() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for samp_entropy() detailed below:

Computer sample entropy (SampEn) of series X, specified by M and R.

SampEn is very close to ApEn.

Suppose given time series is X = [x(1), x(2), ... , x(N)]. We first build embedding matrix Em, of dimension (N-M+1)-by-M, such that the i-th row of Em is x(i),x(i+1), ... , x(i+M-1). Hence, the embedding lag and dimension are 1 and M-1 respectively. Such a matrix can be built by calling pyeeg function as Em = embed_seq(X, 1, M). Then we build matrix Emp, whose only difference with Em is that the length of each embedding sequence is M + 1

Denote the i-th and j-th row of Em as Em[i] and Em[j]. Their k-th elements are Em[i][k] and Em[j][k] respectively. The distance between Em[i] and Em[j] is defined as 1) the maximum difference of their corresponding scalar components, thus, max(Em[i]-Em[j]), or 2) Euclidean distance. We say two 1-D vectors Em[i] and Em[j] match in tolerance R, if the distance between them is no greater than R, thus, max(Em[i]-Em[j]) <= R. Mostly, the value of R is defined as 20% - 30% of standard deviation of X.

Pick Em[i] as a template, for all j such that 0 < j < N - M , we can check whether Em[j] matches with Em[i]. Denote the number of Em[j], which is in the range of Em[i], as k[i], which is the i-th element of the vector k.

We repeat the same process on Emp and obtained Cmp[i], 0 < i < N - M.

The SampEn is defined as log(sum(Cm)/sum(Cmp))

Costa M, Goldberger AL, Peng C-K, Multiscale entropy analysis of biological signals, Physical Review E, 71:021906, 2005

ap_entropy: approximate entropy of a time series

Extremely slow computation. Do NOT use if your dataset is not small and you are not patient enough.

spectral_entropy(obj, *args, **kwargs)¶

Apply spectral_entropy() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for spectral_entropy() detailed below:

Compute spectral entropy of a time series from either two cases below:

X, the time series (default)

2. Power_Ratio, a list of normalized signal power in a set of frequency bins defined in Band (if Power_Ratio is provided, recommended to speed up)

In case 1, Power_Ratio is computed by bin_power() function.

To speed up, it is recommended to compute Power_Ratio before calling this function because it may also be used by other functions whereas computing it here again will slow down.

Band

list

boundary frequencies (in Hz) of bins. They can be unequal bins, e.g. [0.5,4,7,12,30] which are delta, theta, alpha and beta respectively. You can also use range() function of Python to generate equal bins and pass the generated list to this function.

Each element of Band is a physical frequency and shall not exceed the Nyquist frequency, i.e., half of sampling frequency.

X

list

a 1-D real time series.

Fs

integer

the sampling rate in physical frequency

As indicated in return line

bin_power: pyeeg function that computes spectral power in frequency bins

svd_entropy(obj, *args, **kwargs)¶

Apply svd_entropy() in parallel to a list or array

Args:

obj (Sequence of objects or an array) other args are the same as for svd_entropy() detailed below:

Compute SVD Entropy from either two cases below:

a time series X, with lag tau and embedding dimension dE (default)

2. a list, W, of normalized singular values of a matrix (if W is provided, recommend to speed up.)

If W is None, the function will do as follows to prepare singular spectrum:

First, computer an embedding matrix from X, Tau and DE using pyeeg function embed_seq():

M = embed_seq(X, Tau, DE)

Second, use scipy.linalg function svd to decompose the embedding matrix M and obtain a list of singular values:

W = svd(M, compute_uv=0)

At last, normalize W:

W /= sum(W)

To speed up, it is recommended to compute W before calling this function because W may also be used by other functions whereas computing it here again will slow down.

variability_fp(ts, freqs=None, ncycles=6, plot=True)¶

Example variability function. Gives two continuous, time-resolved measures of the variability of a time series, ranging between -1 and 1. The two measures are based on variance of the centroid frequency and variance of the height of the spectral peak, respectively. (Centroid frequency meaning the power-weighted average frequency) These measures are calculated over sliding time windows of variable size. See also: Blenkinsop et al. (2012) The dynamic evolution of focal-onset

epilepsies - combining theoretical and clinical observations

Parameters:	Timeseries of m variables, shape (n, m). Assumed constant timestep. (ts) – (optional) List of frequencies to examine. If None, defaults to (freqs) – 50 frequency bands ranging 1Hz to 60Hz, logarithmically spaced. Window size, in number of cycles of the centroid frequency. (ncycles) – bool Whether to display the output (plot) –
Returns:	variability Timeseries of shape (n, m, 2) variability[:, :, 0] gives a measure of variability between -1 and 1 based on variance of centroid frequency. variability[:, :, 1] gives a measure of variability between -1 and 1 based on variance of maximum power.

exception nsim.nsim.Error¶: Bases: exceptions.Exception

class nsim.nsim.Model¶

Bases: object

Base class for different kinds of dynamical systems

integrate(tspan)¶: numerical integration function to use

class nsim.nsim.MultipleSim(systems, T=60.0, dt=0.005)¶

Bases: object

Represents multiple simulations, possibly running on different hosts

Like a list, indexing with [i] gives access to the ith simulation

timeseries¶: resulting timeseries: all variables of all simulations

output¶: resulting timeseries: output variables of all simulations

output: Rank 3 array representing output time series. 1st axis is time, 2nd axis ranges across output variables of a single simulation, 3rd axis ranges across different simulation instances.

timeseries: Rank 3 array representing multiple time series. 1st axis is time, 2nd axis ranges across all dynamical variables in a single simulation, 3rd axis ranges across different simulation instances.

class nsim.nsim.NetworkSim(systems, T=60.0, dt=0.005)¶

Bases: nsim.nsim.MultipleSim

Simulation of many coupled instances of a model connected in a network

class nsim.nsim.ODEModel¶

Bases: nsim.nsim.Model

Model defined by a system of ordinary differential equations

dimension¶

integer

Dimension of the state space

output_vars¶

list of integers

If i is in this list then y[i] is considered an output variable

f¶

y, t

right hand side of the ODE system dy/dt = f(y, t)

Instance attributes:: y0 (array of shape (ndim,)): Initial state vector

dimension = 1

f(y, t)

integrate(tspan)¶

output_vars = [0]

class nsim.nsim.ParameterSim(systems, T=60.0, dt=0.005)¶

Bases: nsim.nsim.MultipleSim

Independent simulations of a model exploring different parameters

class nsim.nsim.RepeatedSim(model, T=60.0, dt=0.005, repeat=1, identical=True)¶

Bases: nsim.nsim.MultipleSim

Independent simulations of the same model multiple times, with results.

Like a list, indexing the object with [i] gives access to the ith simulation

modelclass¶: the Model class common to all the simulations

timeseries¶: resulting timeseries: all variables of all simulations

output¶: resulting timeseries: output variables of all simulations

class nsim.nsim.SDEModel¶

Bases: nsim.nsim.Model

Model defined by system of (ordinary) stochastic differential equations dy = f(y, t) dt + G(y, t) dW (N.B. currently must be in Stratonovich form)

dimension¶

integer

Dimension of the state space

output_vars¶

list of integers

If i is in this list then y[i] is considered an output variable

f¶

y, t

deterministic part of Stratonovich SDE system

G¶

y, t

noise coefficient matrix of Stratonovich SDE system

Instance attributes:: y0 (array of shape (ndim,)): Initial state vector

G(y, t)

dimension = 1

f(y, t)

integrate(tspan)¶

output_vars = [0]

exception nsim.nsim.SimTypeError¶: Bases: nsim.nsim.Error

exception nsim.nsim.SimValueError¶: Bases: nsim.nsim.Error

class nsim.nsim.Simulation(system, T=60.0, dt=0.005)¶

Bases: object

Represents simulation of a single system and the resulting time series.

system¶

Model

The dynamical system being simulated. (Can provide either a Model subclass or Model instance)

tspan¶

array

The sequence of time points simulated

timeseries¶

array of shape (len(tspan), len(y0))

Multivariate time series of full simulation results.

output¶: Some function of the simulated timeseries, for example a univariate time series of a single output variable.

compute()¶

output: Simulated model output

timeseries: Simulated time series

nsim.nsim.newmodel(f, G, y0, name='NewModel')¶

Use the functions f and G to define a new Model class for simulations.

It will take functions f and G from global scope and make a new Model class out of them. It will automatically gather any globals used in the definition of f and G and turn them into attributes of the new Model.

Parameters:	f – callable(y, t) (defined in global scope) returning (n,) array Scalar or vector-valued function to define the deterministic part G – callable(y, t) (defined in global scope) returning (n,m) array Optional scalar or matrix-valued function to define noise coefficients y0 (Number or array) – Initial condition name (str) – Optional class name for the new model
Returns:	new class (subclass of Model)
Raises:	SimValueError, SimTypeError –

nsim.nsim.newsim(f, G, y0, name='NewModel', T=60.0, dt=0.005, repeat=1, identical=True)¶

Make a simulation of the system defined by functions f and G.

dy = f(y,t)dt + G(y,t).dW with initial condition y0 This helper function is for convenience, making it easy to define one-off simulations interactively in ipython.

Parameters:	f – callable(y, t) (defined in global scope) returning (n,) array Vector-valued function to define the deterministic part of the system G – callable(y, t) (defined in global scope) returning (n,m) array Optional matrix-valued function to define noise coefficients y0 (array) – Initial condition name (str) – Optional class name for the new model T – Total length of time to simulate, in seconds. dt – Timestep for numerical integration. (int, optional) (repeat) – (bool, optional) (identical) –
Returns:	Simulation
Raises:	SimValueError, SimTypeError –

nsim.readfile module¶

functions:: timeseries_from_mat() load a multi-channel Timeseries from a MATLAB .mat file timeseries_from_file() load a multi-channel Timeseries from many file types

nsim.readfile.timeseries_from_file(filename)¶

Load a multi-channel Timeseries from any file type supported by biosig

Supported file formats include EDF/EDF+, BDF/BDF+, EEG, CNT and GDF. Full list is here: http://pub.ist.ac.at/~schloegl/biosig/TESTED

For EDF, EDF+, BDF and BDF+ files, we will use python-edf if it is installed, otherwise will fall back to python-biosig.

Parameters:	filename –
Returns:	Timeseries

nsim.readfile.timeseries_from_mat(filename, varname=None, fs=1.0)¶

load a multi-channel Timeseries from a MATLAB .mat file

Parameters:	filename (str) – .mat file to load varname (str) – variable name. only needed if there is more than one variable saved in the .mat file fs (scalar) – sample rate of timeseries in Hz. (constant timestep assumed)
Returns:	Timeseries

nsim.sde module¶

Provides functions to integrate systems of stochastic differential equations

sodeint integrates Stratonovich SDE systems using the Heun algorithm

exception nsim.sde.Error¶: Bases: exceptions.Exception

exception nsim.sde.SDEValueError¶

Bases: nsim.sde.Error

Thrown if integration arguments fail some basic sanity checks

nsim.sde.sodeint(f, G, y0, tspan)¶

Integrate an (ordinary) stochastic differential equation dy = f(y,t)dt + G(y,t).dW(t) (which must be given in Stratonovich form)

where y is the n-dimensional state vector, f is a vector-valued function, G is an nxm matrix-valued function giving the noise coefficients and dW(t) = (dW_1, dW_2, ... dW_m) is a vector of independent Weiner increments.

Integration is currently done using the Heun algorithm.

Parameters:	f – callable(y, t) returning (n,) array Vector-valued function to define the deterministic part of the system G – callable(y, t) returning (n,m) array Matrix-valued function to define the noise coefficients of the system y0 – array of shape (n,) giving the initial state vector y(t==0) tspan (array) – The sequence of time points for which to solve for y. These must be equally spaced, e.g. np.arange(0,10,0.005) tspan[0] is the intial time corresponding to the initial state y0.
Returns:	array, with shape (len(tspan), len(y0)) With the initial value y0 in the first row
Return type:	y
Raises:	`SDEValueError` –

See also

Mannella (2002) Integration of Stochastic Differential Equations on a Computer

Rumelin (1982) Numerical Treatment of Stochastic Differential Equations

nsim.timeseries module¶

classes:

Timeseries numpy array with extra methods for time series analyses

class nsim.timeseries.Timeseries(input_array, tspan=None, labels=None, fs=None)¶

Bases: numpy.ndarray

A numpy array with extra methods for applying time series analyses. It is an array of up to 3 dimensions.

axis 0 ranges across values of the system at different points in time.: For a single variable time series this is the only axis needed.
axis 1, if present, ranges across the different variables/channels of a: multivariate time series, at a single node.
axis 2, if present, ranges across different nodes of a network simulation: or repeated simulation.

Thus the shape of a Timeseries array is (N, n, m) where N is the total number of time points, n is the number of variables or channels of a single node and m is the number of nodes.

Slice by array index: timeseries[103280:104800] Slice by time in seconds: timeseries.t[10.0:20.5]

All functions defined in the package nsim.analyses1 are available as

methods of the Timeseries (as well as the usual methods of numpy arrays)

Your own analysis functions can also be included by calling

add_analyses('file.py')¶

tspan¶: An array of shape (N,) giving the time in seconds at each time point along axis 0 of the array. tspan will always remain sorted. The time points are not necessarily evenly spaced.

t¶: Allows slicing by time: timeseries.t[starttime:endtime, ...]

labels¶

list

For each axis>0, optionally gives names to label points along the axis. (labels[0] is always None as axis 0 is the time axis, labeled by numbers in tspan rather than by strings) For i>0, labels[i] is a list of str of length shape(timeseries)[i] giving names to each position along axis i. For example labels[1] can hold names for each variable or channel in a multivariate timeseries.

__distob_scatter__(axis=-1, destination=None)¶: Turn a Timeseries into a distributed timeseries

__getitem__(index)¶: When a Timeseries is sliced, tspan will be sliced in the same way as axis 0 of the Timeseries, and labels[i] will be sliced in the same way as axis i of the Timeseries. If the resulting array is not a Timeseries (for example if the time axis is sliced to a scalar value) then returns an ndarray.

__reduce__()¶: Support pickling Timeseries instances by saving __dict__

__setstate__(tsstate)¶: Support unpickling Timeseries instances by loading __dict__

abs()¶: Calculate the absolute value element-wise.

absolute()¶

Calculate the absolute value element-wise.

Returns:	absolute – Absolute value. For complex input (a + bj) gives sqrt(aa + b*2)
Return type:	Timeseries

classmethod add_analyses(source)

Dynamically add new analysis methods to the Timeseries class. :param source: Can be a function, module or the filename of a python file.

If a filename or a module is given, then all functions defined inside not starting with _ will be added as methods.

The only restriction on the functions is that they can accept a Timeseries as their first argument. So existing functions that take a ndarray or array or even a list will usually also work.

angle(deg=False)¶

Return the angle of the complex argument.

Parameters:	deg (bool, optional) – Return angle in degrees if True, radians if False (default).
Returns:	angle – The counterclockwise angle from the positive real axis on the complex plane, with dtype as numpy.float64.
Return type:	Timeseries

ap_entropy(X, M, R)¶

Computer approximate entropy (ApEN) of series X, specified by M and R.

Suppose given time series is X = [x(1), x(2), ... , x(N)]. We first build embedding matrix Em, of dimension (N-M+1)-by-M, such that the i-th row of Em is x(i),x(i+1), ... , x(i+M-1). Hence, the embedding lag and dimension are 1 and M-1 respectively. Such a matrix can be built by calling pyeeg function as Em = embed_seq(X, 1, M). Then we build matrix Emp, whose only difference with Em is that the length of each embedding sequence is M + 1

Denote the i-th and j-th row of Em as Em[i] and Em[j]. Their k-th elements are Em[i][k] and Em[j][k] respectively. The distance between Em[i] and Em[j] is defined as 1) the maximum difference of their corresponding scalar components, thus, max(Em[i]-Em[j]), or 2) Euclidean distance. We say two 1-D vectors Em[i] and Em[j] match in tolerance R, if the distance between them is no greater than R, thus, max(Em[i]-Em[j]) <= R. Mostly, the value of R is defined as 20% - 30% of standard deviation of X.

Pick Em[i] as a template, for all j such that 0 < j < N - M + 1, we can check whether Em[j] matches with Em[i]. Denote the number of Em[j], which is in the range of Em[i], as k[i], which is the i-th element of the vector k. The probability that a random row in Em matches Em[i] is simga_1^{N-M+1} k[i] / (N - M + 1), thus sum(k)/ (N - M + 1), denoted as Cm[i].

We repeat the same process on Emp and obtained Cmp[i], but here 0<i<N-M since the length of each sequence in Emp is M + 1.

The probability that any two embedding sequences in Em match is then sum(Cm)/ (N - M +1 ). We define Phi_m = sum(log(Cm)) / (N - M + 1) and Phi_mp = sum(log(Cmp)) / (N - M ).

And the ApEn is defined as Phi_m - Phi_mp.

Notes

Please be aware that self-match is also counted in ApEn.
This function now runs very slow. We are still trying to speed it up.

References

Costa M, Goldberger AL, Peng CK, Multiscale entropy analysis of biological signals, Physical Review E, 71:021906, 2005

See also

samp_entropy(): sample entropy of a time series

Notes

Extremely slow implementation. Do NOT use if your dataset is not small.

argmax(axis=None, out=None)¶

argmin(axis=None, out=None)¶

autocorrelation(ts, normalized=False, unbiased=False)¶

Returns the discrete, linear convolution of a time series with itself, optionally using unbiased normalization.

N.B. Autocorrelation estimates are necessarily inaccurate for longer lags, as there are less pairs of points to convolve separated by that lag. Therefore best to throw out the results except for shorter lags, e.g. keep lags from tau=0 up to one quarter of the total time series length.

Parameters:	normalized (boolean) – If True, the time series will first be normalized to a mean of 0 and variance of 1. This gives autocorrelation 1 at zero lag. unbiased (boolean) – If True, the result at each lag m will be scaled by 1/(N-m). This gives an unbiased estimation of the autocorrelation of a stationary process from a finite length sample.

Ref: S. J. Orfanidis (1996) “Optimum Signal Processing”, 2nd Ed.

bandpass(ts, low_hz, high_hz, order=3)¶: forward-backward butterworth band-pass filter

bin_power(X, Band, Fs)¶

Compute power in each frequency bin specified by Band from FFT result of X. By default, X is a real signal.

A real signal can be synthesized, thus not real.

Band

list

boundary frequencies (in Hz) of bins. They can be unequal bins, e.g. [0.5,4,7,12,30] which are delta, theta, alpha and beta respectively. You can also use range() function of Python to generate equal bins and pass the generated list to this function.

Each element of Band is a physical frequency and shall not exceed the Nyquist frequency, i.e., half of sampling frequency.

X

list

a 1-D real time series.

Fs

integer

the sampling rate in physical frequency

Returns:	Power – list spectral power in each frequency bin. Power_ratio – list spectral power in each frequency bin normalized by total power in ALL frequency bins.

circmean(ts, axis=2)¶: Circular mean phase

circstd(ts, axis=2)¶: Circular standard deviation

concatenate(tup, axis=0)¶

Join a sequence of Timeseries to this one :param tup: timeseries to be joined with this one.

They must have the same shape as this Timeseries, except in the dimension corresponding to axis.

Parameters:	axis (int, optional) – The axis along which timeseries will be joined.
Returns:	res (Timeseries or ndarray)

crossing_indices(ts, c=0.0, d=0.0)¶

For a single variable timeseries, find the time indices each time the value crosses c from above or below. Can optionally set a non-zero d to impose the condition that the value must wander at least d units away from c between crossings.

If the timeseries begins (or ends) exactly at c, then time zero (or the ending time) is also included as a crossing event, so that the boundaries of the first and last excursions are included.

Parameters:	ts – Timeseries (single variable) c (float) – Critical value at which to report crossings. d (float) – Optional min distance from c to be attained between crossings.
Returns:	array of indices

cwt(ts, freqs=array([ 1., 1.09854114, 1.20679264, 1.32571137, 1.45634848, 1.59985872, 1.75751062, 1.93069773, 2.12095089, 2.32995181, 2.55954792, 2.8117687, 3.0888436, 3.39322177, 3.72759372, 4.09491506, 4.49843267, 4.94171336, 5.42867544, 5.96362332, 6.55128557, 7.19685673, 7.90604321, 8.68511374, 9.54095476, 10.48113134, 11.51395399, 12.64855217, 13.89495494, 15.26417967, 16.76832937, 18.42069969, 20.23589648, 22.22996483, 24.42053095, 26.82695795, 29.47051703, 32.37457543, 35.56480306, 39.06939937, 42.9193426, 47.14866363, 51.79474679, 56.89866029, 62.50551925, 68.6648845, 75.43120063, 82.86427729, 91.0298178, 100. ]), wavelet=<function cwtmorlet>, plot=True)¶

Continuous wavelet transform Note the full results can use a huge amount of memory at 64-bit precision

Parameters:	ts – Timeseries of m variables, shape (n, m). Assumed constant timestep. freqs – list of frequencies (in Hz) to use for the tranform. (default is 50 frequency bins logarithmic from 1Hz to 100Hz) wavelet – the wavelet to use. may be complex. see scipy.signal.wavelets plot – whether to plot time-resolved power spectrum
Returns:	Continuous wavelet transform output array, shape (n,len(freqs),m)
Return type:	coefs

cwt_distributed(ts, freqs=array([ 1., 1.09854114, 1.20679264, 1.32571137, 1.45634848, 1.59985872, 1.75751062, 1.93069773, 2.12095089, 2.32995181, 2.55954792, 2.8117687, 3.0888436, 3.39322177, 3.72759372, 4.09491506, 4.49843267, 4.94171336, 5.42867544, 5.96362332, 6.55128557, 7.19685673, 7.90604321, 8.68511374, 9.54095476, 10.48113134, 11.51395399, 12.64855217, 13.89495494, 15.26417967, 16.76832937, 18.42069969, 20.23589648, 22.22996483, 24.42053095, 26.82695795, 29.47051703, 32.37457543, 35.56480306, 39.06939937, 42.9193426, 47.14866363, 51.79474679, 56.89866029, 62.50551925, 68.6648845, 75.43120063, 82.86427729, 91.0298178, 100. ]), wavelet=<function cwtmorlet>, plot=True)¶

Continuous wavelet transform using distributed computation. (Currently just splits the data by channel. TODO split it further.) Note: this function requires an IPython cluster to be started first.

Parameters:	ts – Timeseries of m variables, shape (n, m). Assumed constant timestep. freqs – list of frequencies (in Hz) to use for the tranform. (default is 50 frequency bins logarithmic from 1Hz to 100Hz) wavelet – the wavelet to use. may be complex. see scipy.signal.wavelets plot – whether to plot time-resolved power spectrum
Returns:	Continuous wavelet transform output array, shape (n,len(freqs),m)
Return type:	coefs

dfa(X, Ave=None, L=None)¶

Compute Detrended Fluctuation Analysis from a time series X and length of boxes L.