Stimulus

This module provides convenient functions for synthesizing timeseries used as the input to the network.

hetero.stimulus is the dispatcher: it looks the label given on the command line (--stim) up in the stim block of hetero.config.CFG, calls the matching generator with that block’s parameters, and concatenates the requested timeseries into the single input u the network is driven with. The generators themselves are grouped into three families, one module each:

hetero.stims.elementary

Trivial shapes – a null signal and a constant offset.

hetero.stims.maps

Iterated one- and two-dimensional discrete maps.

hetero.stims.processes

Continuous-time dynamical systems integrated onto a discrete time series.

Note

Because the map and process generators are the workhorses of the input benchmark, they follow the conventions of the ReservoirPy (https://github.com/reservoirpy/reservoirpy) package, with some adaptation for cross-modular compatibility.

hetero.stimulus.synthesize(name, n_timesteps, dt_scaler, **stim_params)

Synthesizes timeseries from the specified process.

Parameters:

• name (str) – The name of the underlying process. See notes for valid names.

• n_timesteps (int) – Number of timesteps.

• dt_scaler (float) – A scaling factor to up/down scale the integrator’s timestep (only relevant for processes.)

Returns:

ndarray – Synthesized array

Notes

All supported generative processes can be categorized into three categories:

• elementary functions: simple functions of time

• maps: functions of previous step (without explicit notion of physical time)

• processes: (system of) ODEs with explicit notion of time

The valid processes are listed below. Also see their respective documentation for potential arguments.

Elementary functions

• sin

• sign

• null

• constant

Maps

• henon_map

• logistic_map

• narma

Processes

• lorenz

• lorenz96

• mackey_glass

• multiscroll

• doublescroll

• rabinovich_fabrikant

• rossler

• kuramoto_sivashinsky

hetero.stimulus.is_constant_stim(stim_params)

Checks whether the stimulus configuration is a pure constant (DC) input.

A constant (DC) input carries no characteristic timescale. Hence, its temporal characterization (base-lag estimation) and the subsequent time rescaling are irrelevant and can be bypassed.

Parameters:

stim_params (dict) – Stimulus parameters, e.g. {'constant': {'amplitude': 1.}} or a combination thereof, e.g. {'constant': ..., 'constant-4': ...}.

Returns:

bool – True if every (compound) component is a constant input.

Return type:

bool

Elementary functions

hetero.stims.elementary.null(n_timesteps, **kwargs)

Generats a null (zero) stimulus.

Parameters:

• n_timesteps (int) – Number of timesteps to generate.

• kwargs (dict)

Returns:

np.ndarray of shape (n_timesteps, 1) – A zero 1-dimensional array.

Return type:

ndarray

hetero.stims.elementary.sin(n_timesteps, h=0.05, freq=1., phase=0.)

Generates a sinusoidal input of given frequency.

Parameters:

• n_timesteps (int) – Number of timesteps to generate.

• h (float, optional) – Time delta between two discrete timesteps, by default 0.05.

• freq (float, optional) – Oscillation frequency, by default 1.

• phase (float, optional) – Oscillation starting phase in radian, by default 0.

Returns:

np.ndarray of shape (n_timesteps, 1) – A 1-dimensional array.

Return type:

ndarray

hetero.stims.elementary.sign(n_timesteps, h=0.05, freq=1., phase=0.)

Generates a signal that is the sign of a sinusoidal input of given frequency.

Parameters:

• n_timesteps (int) – Number of timesteps to generate.

• h (float, optional) – Time delta between two discrete timesteps, by default 0.05.

• freq (float, optional) – Oscillation frequency, by default 1.

• phase (float, optional) – Oscillation starting phase in radian, by default 0.

Returns:

np.ndarray of shape (n_timesteps, 1) – A 1-dimensional array.

Return type:

ndarray

hetero.stims.elementary.sinabs(n_timesteps, h=0.05, freq=1., phase=0.)

Generates the abs of a sinusoidal input of given frequency.

Parameters:

• n_timesteps (int) – Number of timesteps to generate.

• h (float, optional) – Time delta between two discrete timesteps, by default 0.05.

• freq (float, optional) – Oscillation frequency, by default 1.

• phase (float, optional) – Oscillation starting phase in radian, by default 0.

Returns:

np.ndarray of shape (n_timesteps, 1) – A 1-dimensional array.

Return type:

ndarray

hetero.stims.elementary.constant(n_timesteps, amplitude=1., h=0.05, **kwargs)

Generates a constant amplitude input

Parameters:

• n_timesteps (int) – Number of timesteps to generate.

• amplitude (float, optional) – Input amplitude.

• h (float, optional) – Time delta between two discrete timesteps, by default 0.05. (unused: a constant input does not depend on time)

• kwargs (dict)

Returns:

np.ndarray of shape (n_timesteps, 1) – A 1-dimensional array.

Return type:

ndarray

Maps

Note

This section is directly reproduced from the ReservoirPy (https://github.com/reservoirpy/reservoirpy) package, with some adaptation for cross-modular compatibility.

hetero.stims.maps.henon_map(n_timesteps, a=1.4, b=0.3, x0=[0.0, 0.0], **kwargs)

Hénon map discrete timeseries [2] [3].

\[\begin{split}x(n+1) &= 1 - ax(n)^2 + y(n)\\ y(n+1) &= bx(n)\end{split}\]

Parameters:

• n_timesteps (int) – Number of timesteps to generate.

• a (float, default to 1.4) – \(a\) parameter of the system.

• b (float, default to 0.3) – \(b\) parameter of the system.

• x0 (array-like of shape (2,), default to [0.0, 0.0]) – Initial conditions of the system.

Returns:

array of shape (n_timesteps, 2) – Hénon map discrete timeseries.

Return type:

ndarray

hetero.stims.maps.logistic_map(n_timesteps, r=3.9, x0=0.5, **kwargs)

Logistic map discrete timeseries [4] [5].

\[x(n+1) = rx(n)(1-x(n))\]

Parameters:

• n_timesteps (int) – Number of timesteps to generate.

• r (float, default to 3.9) – \(r\) parameter of the system.

• x0 (float, default to 0.5) – Initial condition of the system.

Returns:

array of shape (n_timesteps, 1) – Logistic map discrete timeseries.

Return type:

ndarray

hetero.stims.maps.narma(n_timesteps, order=30, a1=0.2, a2=0.04, b=1.5, c=0.001, x0=[0.0], seed=None, h=1.0)

Non-linear Autoregressive Moving Average (NARMA) timeseries, as first defined in [14], and as used in [15].

NARMA n-th order dynamical system is defined by the recurrent relation:

\[y[t+1] = a_1 y[t] + a_2 y[t] (\sum_{i=0}^{n-1} y[t-i]) + b u[t-( n-1)]u[t] + c\]

where \(u[t]\) are sampled following a uniform distribution in \([0, 0.5]\).

Parameters:

• n_timesteps (int) – Number of timesteps to generate.

• order (int, default to 30) – Order of the system.

• a1 (float, default to 0.2) – \(a_1\) parameter of the system.

• a2 (float, default to 0.04) – \(a_2\) parameter of the system.

• b (float, default to 1.5) – \(b\) parameter of the system.

• c (float, default to 0.001) – \(c\) parameter of the system.

• x0 (array-like of shape (init_steps,), default to [0.0]) – Initial conditions of the system.

• seed (int or numpy.random.Generator, optional) – Random state seed for reproducibility.

• h (float)

Returns:

array of shape (n_timesteps, 1) – NARMA timeseries.

Return type:

ndarray

Processes

Note

This section is directly reproduced from the ReservoirPy (https://github.com/reservoirpy/reservoirpy) package, with some adaptation for cross-modular compatibility.

hetero.stims.processes.lorenz(n_timesteps, rho=28.0, sigma=10.0, beta=8.0 / 3.0, x0=[-15.0577155, -19.28572772, 31.35055511], h=0.03, **kwargs)

Lorenz attractor timeseries as defined by Lorenz in 1963 [6] [7].

\[\begin{split}\frac{\mathrm{d}x}{\mathrm{d}t} &= \sigma (y-x) \\ \frac{\mathrm{d}y}{\mathrm{d}t} &= x(\rho - z) - y \\ \frac{\mathrm{d}z}{\mathrm{d}t} &= xy - \beta z\end{split}\]

Parameters:

• n_timesteps (int) – Number of timesteps to generate.

• rho (float, default to 28.0) – \(\rho\) parameter of the system.

• sigma (float, default to 10.0) – \(\sigma\) parameter of the system.

• beta (float, default to 8/3) – \(\beta\) parameter of the system.

• x0 (array-like of shape (3,), default to [1.0, 1.0, 1.0]) – Initial conditions of the system.

• h (float, default to 0.03) – Time delta between two discrete timesteps.

• **kwargs – Other parameters to pass to the scipy.integrate.solve_ivp solver.

Returns:

array of shape (n_timesteps, 3) – Lorenz attractor timeseries.

Return type:

ndarray

hetero.stims.processes.mackey_glass(n_timesteps, tau=17, a=0.2, b=0.1, n=10, x0=1.2, h=1.0, seed=None, **kwargs)

Mackey-Glass timeseries [8] [9], computed from the Mackey-Glass delayed differential equation.

\[\frac{x}{t} = \frac{ax(t-\tau)}{1+x(t-\tau)^n} - bx(t)\]

Parameters:

• n_timesteps (int) – Number of timesteps to compute.

• tau (int, default to 17) – Time delay \(\tau\) of Mackey-Glass equation. By defaults, equals to 17. Other values can change the choatic behaviour of the timeseries.

• a (float, default to 0.2) – \(a\) parameter of the equation.

• b (float, default to 0.1) – \(b\) parameter of the equation.

• n (int, default to 10) – \(n\) parameter of the equation.

• x0 (float, optional, default to 1.2) – Initial condition of the timeseries.

• h (float, default to 1.0) – Time delta between two discrete timesteps.

• seed (int or numpy.random.Generator, optional) – Random state seed for reproducibility.

Returns:

array of shape (n_timesteps, 1) – Mackey-Glass timeseries.

Return type:

ndarray

Note

As Mackey-Glass is defined by delayed time differential equations, the first timesteps of the timeseries can’t be initialized at 0 (otherwise, the first steps of computation involving these not-computed-yet-timesteps would yield inconsistent results). A random number generator is therefore used to produce random initial timesteps based on the value of the initial condition passed as parameter. A default seed is hard-coded to ensure reproducibility in any case. It can be changed with the set_seed() function.

hetero.stims.processes.multiscroll(n_timesteps, a=40.0, b=3.0, c=28.0, x0=[-0.1, 0.5, -0.6], h=0.01, **kwargs)

Double scroll attractor timeseries [10] [11], a particular case of multiscroll attractor timeseries.

\[\begin{split}\frac{\mathrm{d}x}{\mathrm{d}t} &= a(y - x) \\ \frac{\mathrm{d}y}{\mathrm{d}t} &= (c - a)x - xz + cy \\ \frac{\mathrm{d}z}{\mathrm{d}t} &= xy - bz\end{split}\]

Parameters:

• n_timesteps (int) – Number of timesteps to generate.

• a (float, default to 40.0) – \(a\) parameter of the system.

• b (float, default to 3.0) – \(b\) parameter of the system.

• c (float, default to 28.0) – \(c\) parameter of the system.

• x0 (array-like of shape (3,), default to [-0.1, 0.5, -0.6]) – Initial conditions of the system.

• h (float, default to 0.01) – Time delta between two discrete timesteps.

Returns:

array of shape (n_timesteps, 3) – Multiscroll attractor timeseries.

Return type:

ndarray

hetero.stims.processes.doublescroll(n_timesteps, r1=1.2, r2=3.44, r4=0.193, ir=2 * 2.25e-5, beta=11.6, x0=[0.37926545, 0.058339, -0.08167691], h=0.25, **kwargs)

Double scroll attractor timeseries [10] [11], a particular case of multiscroll attractor timeseries.

\[\begin{split}\frac{\mathrm{d}V_1}{\mathrm{d}t} &= \frac{V_1}{R_1} - \frac{\Delta V}{R_2} - 2I_r \sinh(\beta\Delta V) \\ \frac{\mathrm{d}V_2}{\mathrm{d}t} &= \frac{\Delta V}{R_2} +2I_r \sinh(\beta\Delta V) - I\\ \frac{\mathrm{d}I}{\mathrm{d}t} &= V_2 - R_4 I\end{split}\]

where \(\Delta V = V_1 - V_2\).

Parameters:

• n_timesteps (int) – Number of timesteps to generate.

• r1 (float, default to 1.2) – \(R_1\) parameter of the system.

• r2 (float, default to 3.44) – \(R_2\) parameter of the system.

• r4 (float, default to 0.193) – \(R_4\) parameter of the system.

• ir (float, default to 2*2e.25e-5) – \(I_r\) parameter of the system.

• beta (float, default to 11.6) – \(\beta\) parameter of the system.

• x0 (array-like of shape (3,), default to [0.37926545, 0.058339, -0.08167691]) – Initial conditions of the system.

• h (float, default to 0.01) – Time delta between two discrete timesteps.

Returns:

array of shape (n_timesteps, 3) – Multiscroll attractor timeseries.

Return type:

ndarray

hetero.stims.processes.rabinovich_fabrikant(n_timesteps, alpha=1.1, gamma=0.89, x0=[-1, 0, 0.5], h=0.05, **kwargs)

Rabinovitch-Fabrikant system [12] [13] timeseries.

\[\begin{split}\frac{\mathrm{d}x}{\mathrm{d}t} &= y(z - 1 + x^2) + \gamma x \\ \frac{\mathrm{d}y}{\mathrm{d}t} &= x(3z + 1 - x^2) + \gamma y \\ \frac{\mathrm{d}z}{\mathrm{d}t} &= -2z(\alpha + xy)\end{split}\]

Parameters:

• n_timesteps (int) – Number of timesteps to generate.

• alpha (float, default to 1.1) – \(\alpha\) parameter of the system.

• gamma (float, default to 0.89) – \(\gamma\) parameter of the system.

• x0 (array-like of shape (3,), default to [-1, 0, 0.5]) – Initial conditions of the system.

• h (float, default to 0.05) – Time delta between two discrete timesteps.

• **kwargs – Other parameters to pass to the scipy.integrate.solve_ivp solver.

Returns:

array of shape (n_timesteps, 3) – Rabinovitch-Fabrikant system timeseries.

Return type:

ndarray

hetero.stims.processes.lorenz96(n_timesteps, warmup=0, N=36, F=8.0, dF=0.01, h=0.01, x0=None, **kwargs)

Lorenz96 attractor timeseries as defined by Lorenz in 1996 [17].

\[\frac{\mathrm{d}x_i}{\mathrm{d} t} = (x_{i+1} - x_{i-2}) x_{i-1} - x_i + F\]

where \(i = 1, \dots, N\) and \(x_{-1} = x_{N-1}\) and \(x_{N+1} = x_1\) and \(N \geq 4\).

Parameters:

• n_timesteps (int) – Number of timesteps to generate.

• warmup (int, default to 0) – Number of timesteps to discard at the begining of the signal, to remove transient states.

• N (int, default to 36) – Dimension of the system.

• F (float, default to F) – \(F\) parameter of the system.

• dF (float, default to 0.01) – Pertubation applied to initial condition if x0 is None.

• h (float, default to 0.01) – Time delta between two discrete timesteps.

• x0 (array-like of shape (N,), default to None) – Initial conditions of the system. If None, the array is initialized to an array of shape (N, ) with value F, except for the first value of the array that takes the value F + dF.

• **kwargs – Other parameters to pass to the scipy.integrate.solve_ivp solver.

Returns:

array of shape (n_timesteps - warmup, N) – Lorenz96 timeseries.

Return type:

ndarray

hetero.stims.processes.rossler(n_timesteps, a=0.2, b=0.2, c=5.7, x0=[-0.1, 0.0, 0.02], h=0.1, **kwargs)

Rössler attractor timeseries [18].

\[\begin{split}\frac{\mathrm{d}x}{\mathrm{d}t} &= -y - z \\ \frac{\mathrm{d}y}{\mathrm{d}t} &= x + a y \\ \frac{\mathrm{d}z}{\mathrm{d}t} &= b + z (x - c)\end{split}\]

Parameters:

• n_timesteps (int) – Number of timesteps to generate.

• a (float, default to 0.2) – \(a\) parameter of the system.

• b (float, default to 0.2) – \(b\) parameter of the system.

• c (float, default to 5.7) – \(c\) parameter of the system.

• x0 (array-like of shape (3,), default to [-0.1, 0.0, 0.02]) – Initial conditions of the system.

• h (float, default to 0.1) – Time delta between two discrete timesteps.

• **kwargs – Other parameters to pass to the scipy.integrate.solve_ivp solver.

Returns:

array of shape (n_timesteps, 3) – Rössler attractor timeseries.

Return type:

ndarray

hetero.stims.processes.kuramoto_sivashinsky(n_timesteps, warmup=0, N=128, M=16, x0=None, h=0.25)

Kuramoto-Sivashinsky oscillators [19] [20] [21].

\[y_t = -yy_x - y_{xx} - y_{xxxx}, ~~ x \in [0, 32\pi]\]

This 1D partial differential equation is solved using ETDRK4 (Exponential Time-Differencing 4th order Runge-Kutta) method, as described in [22].

Parameters:

• n_timesteps (int) – Number of timesteps to generate.

• warmup (int, default to 0) – Number of timesteps to discard at the begining of the signal, to remove transient states.

• N (int, default to 128) – Dimension of the system.

• M (float, default to 0.2) – Number of points for complex means. Modify beahviour of the resulting multivariate timeseries.

• x0 (array-like of shape (N,), default to None.) – Initial conditions of the system. If None, x0 is equal to \(\cos (\frac{y}{M}) * (1 + \sin(\frac{y}{M}))\) with \(y = 2M\pi x / N, ~~ x \in [1, N]\).

• h (float, default to 0.25) – Time delta between two discrete timesteps.

Returns:

array of shape (n_timesteps - warmup, N) – Kuramoto-Sivashinsky equation solution.

Return type:

ndarray

Seeding helpers

hetero.stims.utils.get_seed(seed=None)

hetero.stims.utils.rand_generator(seed=None)

Parameters:

seed (int | Generator | RandomState)

Return type:

Generator

hetero.stims.utils.check_vector(array, allow_reshape=True, allow_timespans=True, caller=None)

References

[2]

M. Hénon, ‘A two-dimensional mapping with a strange attractor’, Comm. Math. Phys., vol. 50, no. 1, pp. 69–77, 1976.

[3]

Hénon map on Wikipédia

[4]

R. M. May, ‘Simple mathematical models with very complicated dynamics’, Nature, vol. 261, no. 5560, Art. no. 5560, Jun. 1976, doi: 10.1038/261459a0.

[5]

Logistic map on Wikipédia

[6]

E. N. Lorenz, ‘Deterministic Nonperiodic Flow’, Journal of the Atmospheric Sciences, vol. 20, no. 2, pp. 130–141, Mar. 1963, doi: 10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2.

[7]

Lorenz system on Wikipedia.

[8]

M. C. Mackey and L. Glass, ‘Oscillation and chaos in physiological control systems’, Science, vol. 197, no. 4300, pp. 287–289, Jul. 1977, doi: 10.1126/science.267326.

[9]

Mackey-Glass equations on Wikipedia.

[10] (1,2)

G. Chen and T. Ueta, ‘Yet another chaotic attractor’, Int. J. Bifurcation Chaos, vol. 09, no. 07, pp. 1465–1466, Jul. 1999, doi: 10.1142/S0218127499001024.

[11] (1,2)

Chen double scroll attractor on Wikipedia.

[12]

M. I. Rabinovich and A. L. Fabrikant, ‘Stochastic self-modulation of waves in nonequilibrium media’, p. 8, 1979.

[13]

Rabinovich-Fabrikant equations on Wikipedia.

[14]

A. F. Atiya and A. G. Parlos, ‘New results on recurrent network training: unifying the algorithms and accelerating convergence,’ in IEEE Transactions on Neural Networks, vol. 11, no. 3, pp. 697-709, May 2000, doi: 10.1109/72.846741.

[15]

B. Schrauwen, M. Wardermann, D. Verstraeten, J. Steil, D. Stroobandt, ‘Improving reservoirs using intrinsic plasticity’, Neurocomputing, 71. 1159-1171, 2008, doi: 10.1016/j.neucom.2007.12.020.

[17]

Lorenz, E. N. (1996, September). Predictability: A problem partly solved. In Proc. Seminar on predictability (Vol. 1, No. 1).

[18]

O.E. Rössler, “An equation for continuous chaos”, Physics Letters A, vol 57, Issue 5, Pages 397-398, ISSN 0375-9601, 1976, https://doi.org/10.1016/0375-9601(76)90101-8.

[19]

Kuramoto, Y. (1978). Diffusion-Induced Chaos in Reaction Systems. Progress of Theoretical Physics Supplement, 64, 346–367. https://doi.org/10.1143/PTPS.64.346

[20]

Sivashinsky, G. I. (1977). Nonlinear analysis of hydrodynamic instability in laminar flames—I. Derivation of basic equations. Acta Astronautica, 4(11), 1177–1206. https://doi.org/10.1016/0094-5765(77)90096-0

[21]

Sivashinsky, G. I. (1980). On Flame Propagation Under Conditions of Stoichiometry. SIAM Journal on Applied Mathematics, 39(1), 67–82. https://doi.org/10.1137/0139007

[22]

Kassam, A. K., & Trefethen, L. N. (2005). Fourth-order time-stepping for stiff PDEs. SIAM Journal on Scientific Computing, 26(4), 1214-1233.