Stimulus

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

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:

Synthesized array

Return type:

ndarray

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

Maps

• henon_map

• logistic_map

• narma

Processes

• lorenz

• lorenz96

• mackey_glass

• multiscroll

• doublescroll

• rabinovich_fabrikant

• rossler

• kuramoto_sivashinsky

Elementary Functions

hetero.stims.elem_funcs.null(n_timesteps: int, **kwargs: dict) → ndarray

Generats a null (zero) stimulus.

Parameters:

n_timesteps (int) – Number of timesteps to generate.

Returns:

A zero 1-dimensional array.

Return type:

np.ndarray of shape (n_timesteps, 1)

hetero.stims.elem_funcs.sin(n_timesteps: int, h: float = 0.05, freq: float = 10.0, phase: float = 0.0) → ndarray

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 10.

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

Returns:

A 1-dimensional array.

Return type:

np.ndarray of shape (n_timesteps, 1)

hetero.stims.elem_funcs.sign(n_timesteps: int, h: float = 0.05, freq: float = 1.0, phase: float = 0.0) → ndarray

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 10.

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

Returns:

A 1-dimensional array.

Return type:

np.ndarray of shape (n_timesteps, 1)

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: int, a: float = 1.4, b: float = 0.3, x0: list | ndarray = [0.0, 0.0], **kwargs) → ndarray

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:

Hénon map discrete timeseries.

Return type:

array of shape (n_timesteps, 2)

hetero.stims.maps.logistic_map(n_timesteps: int, r: float = 3.9, x0: float = 0.5, **kwargs) → ndarray

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:

Logistic map discrete timeseries.

Return type:

array of shape (n_timesteps, 1)

hetero.stims.maps.narma(n_timesteps: int, order: int = 30, a1: float = 0.2, a2: float = 0.04, b: float = 1.5, c: float = 0.001, x0: list | ndarray = [0.0], seed: int | RandomState = None) → ndarray

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.

Returns:

NARMA timeseries.

Return type:

array of shape (n_timesteps, 1)

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: int, rho: float = 28.0, sigma: float = 10.0, beta: float = 2.6666666666666665, x0: list | ndarray = [-15.0577155, -19.28572772, 31.35055511], h: float = 0.03, **kwargs) → ndarray

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:

Lorenz attractor timeseries.

Return type:

array of shape (n_timesteps, 3)

hetero.stims.processes.mackey_glass(n_timesteps: int, tau: int = 17, a: float = 0.2, b: float = 0.1, n: int = 10, x0: float = 1.2, h: float = 1.0, seed: int | RandomState | Generator = None, **kwargs) → ndarray

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:

Mackey-Glass timeseries.

Return type:

array of shape (n_timesteps, 1)

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: int, a: float = 40.0, b: float = 3.0, c: float = 28.0, x0: list | ndarray = [-0.1, 0.5, -0.6], h: float = 0.01, **kwargs) → ndarray

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:

Multiscroll attractor timeseries.

Return type:

array of shape (n_timesteps, 3)

hetero.stims.processes.doublescroll(n_timesteps: int, r1: float = 1.2, r2: float = 3.44, r4: float = 0.193, ir: float = 4.5e-05, beta: float = 11.6, x0: list | ndarray = [0.37926545, 0.058339, -0.08167691], h: float = 0.25, **kwargs) → ndarray

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:

Multiscroll attractor timeseries.

Return type:

array of shape (n_timesteps, 3)

hetero.stims.processes.rabinovich_fabrikant(n_timesteps: int, alpha: float = 1.1, gamma: float = 0.89, x0: list | ndarray = [-1, 0, 0.5], h: float = 0.05, **kwargs) → ndarray

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:

Rabinovitch-Fabrikant system timeseries.

Return type:

array of shape (n_timesteps, 3)

hetero.stims.processes.lorenz96(n_timesteps: int, warmup: int = 0, N: int = 36, F: float = 8.0, dF: float = 0.01, h: float = 0.01, x0: list | ndarray = None, **kwargs) → ndarray

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:

Lorenz96 timeseries.

Return type:

array of shape (n_timesteps - warmup, N)

hetero.stims.processes.rossler(n_timesteps: int, a: float = 0.2, b: float = 0.2, c: float = 5.7, x0: list | ndarray = [-0.1, 0.0, 0.02], h: float = 0.1, **kwargs) → ndarray

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:

Rössler attractor timeseries.

Return type:

array of shape (n_timesteps, 3)

hetero.stims.processes.kuramoto_sivashinsky(n_timesteps: int, warmup: int = 0, N: int = 128, M: float = 16, x0: list | ndarray = None, h: float = 0.25) → ndarray

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:

Kuramoto-Sivashinsky equation solution.

Return type:

array of shape (n_timesteps - warmup, N)

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.