rc-prediction 0.1.0

Parameter-aware reservoir computing for critical transition and system collapse prediction

rc_prediction

Model-free prediction of critical transitions and system collapse using parameter-aware reservoir computing.

Python implementation of Kong, Fan, Grebogi & Lai, Machine learning prediction of critical transition and system collapse, Phys. Rev. Research 3, 013090 (2021).

---

Contents

• Overview
• Installation
• Quick start
• Method
• Package structure
• API reference
• Benchmark systems
• Tutorial and examples
• Tests
• References
• License

---

Overview

When a control parameter $p$ crosses a critical value $p_c$, many nonlinear systems switch from sustained chaos to transient escape:

Regime	Condition	Typical behavior
Pre-critical	$p < p_c$	Trajectory stays on the attractor
Post-critical	$p > p_c$	Transient chaos, then collapse / escape

This package trains on time series from at least three pre-critical parameter values (equations unknown) and predicts:

1. Whether a test parameter produces collapse
2. An estimated critical point $p_c^{*}$
3. Transient lifetimes beyond $p_c$

Main entry point: ParameterAwareRC

---

Installation

pip install rc_prediction

Optional extras:

pip install "rc_prediction[examples]"   # matplotlib, Jupyter, ipykernel
pip install "rc_prediction[dev]"        # pytest, ruff, mypy, build, twine

Editable install from source:

git clone https://github.com/jinchen7-cmd/Reservoir-Computing.git
cd Reservoir-Computing
pip install -e ".[dev]"

• Python: $\geq$ 3.10
• Dependencies: NumPy, SciPy, scikit-learn
• PyPI: https://pypi.org/project/rc-prediction/

---

Quick start

Ikeda map

from rc_prediction import IkedaMap, ParameterAwareRC
from rc_prediction.systems.ikeda import DEFAULT_TRAINING_MU, MU_CRITICAL

ikeda = IkedaMap()
data = ikeda.simulate_training_set(
    DEFAULT_TRAINING_MU, n_steps=1500, burn_in=500, random_state=42
)

model = ParameterAwareRC(n_units=300, random_state=7)
model.fit(data, parameter_name="mu", train_length=800, collapse_bound=6.0)

safe = model.predict_closed_loop(0.99, n_steps=2000)
collapse = model.predict_closed_loop(1.01, n_steps=2000)
scan = model.scan_critical_point((0.98, 1.02), n_points=15, n_steps=1500)

print(MU_CRITICAL, scan.p_critical, safe.collapsed, collapse.collapsed)

Food chain

from rc_prediction import FoodChain, ParameterAwareRC
from rc_prediction.systems.food_chain import DEFAULT_TRAINING_K, K_CRITICAL

food = FoodChain()
data = food.simulate_training_set(
    DEFAULT_TRAINING_K, t_max=1200.0, dt=1.0, burn_in=600.0, random_state=42
)

model = ParameterAwareRC(n_units=400, random_state=3)
model.fit(
    data,
    parameter_name="K",
    train_length=400,
    predator_index=food.predator_index,
)

result = model.predict_closed_loop(1.01, n_steps=1500)
print(K_CRITICAL, result.collapsed)

Standard ESN

from rc_prediction import ESN
from rc_prediction.utils import lorenz_system, train_test_split_sequence

series = lorenz_system(5000)
X, y = series[:-1], series[1:]
X_train, X_test, y_train, y_test = train_test_split_sequence(X, y)

model = ESN(n_units=500, spectral_radius=0.9, leaking_rate=0.3)
model.fit(X_train, y_train, warmup=500)
print(model.score(X_test, y_test, warmup=0))

---

Method

Reservoir dynamics

State $\mathbf{u}(t) \in \mathbb{R}^d$ and bifurcation parameter $p$ are combined into one input vector (see arc/core.py):

$$ \tilde{\mathbf{u}}(t) = \begin{bmatrix} \mathbf{u}(t) \ k_b (p + b_0) \end{bmatrix} $$

Leaky echo-state update:

$$ \mathbf{r}(t+\Delta t) = (1-\alpha)\mathbf{r}(t) + \alpha \tanh\left(\mathbf{W}\mathbf{r}(t) + \mathbf{W}_{in}\tilde{\mathbf{u}}(t)\right) $$

Readout with squared even-index reservoir units (0-based indices $1, 3, 5, \ldots$):

$$ \phi(\mathbf{r})_i = \begin{cases} r_i^2 & i \text{ odd (0-based)} \ r_i & \text{otherwise} \end{cases} $$

$$ \mathbf{v}(t) = \mathbf{W}_{out},\phi(\mathbf{r}(t)) $$

Component	Role
$\mathbf{W}_{out}$	Trained by ridge regression (one-step: $\mathbf{v}(t) \approx \mathbf{u}(t+\Delta t)$)
$\mathbf{W}_{in}$, $\mathbf{W}$	Fixed after random initialization
$k_b$, $b_0$, $\alpha$, $\lambda$	Hyperparameters

Training data

ParameterAwareRC.fit expects a dictionary mapping parameter values to arrays of shape (n_steps, n_dims):

Rule	Requirement
Minimum parameters	$\geq 3$ distinct values (InsufficientParameterValues otherwise)
Regime	All training values pre-critical ($p < p_c$)
Series length	$\geq 20$ timesteps per parameter

Training concatenates teacher-forced segments from all parameters and fits a single $\mathbf{W}_{out}$.

Closed-loop prediction

At test parameter $p_{\mathrm{test}}$:

$$ \mathbf{u}(t+\Delta t) = \mathbf{v}(t) $$

Collapse detection (arc/predictor.py):

• Ikeda / bounded systems: collapse_bound — any state component exceeds the bound
• Food chain: predator_index — predator density drops below a threshold

Primary methods on ParameterAwareRC:

Method	Returns	Purpose
predict_closed_loop(p_test, ...)	ClosedLoopResult	Autonomous rollout at one $p$
scan_critical_point(p_range, ...)	ScanResult	Sweep $p$ and estimate $p_c^{*}$
ensemble_predict(p_test, ...)	EnsembleResult	Average over reservoir seeds

---

Package structure

Reservoir-Computing/
├── pyproject.toml
├── LICENSE
├── README.md
├── .github/
│   └── workflows/
│       └── publish-pypi.yml
├── src/
│   └── rc_prediction/
│       ├── __init__.py           # top-level public API
│       ├── base.py
│       ├── reservoir.py
│       ├── readout.py
│       ├── esn.py
│       ├── topology.py
│       ├── metrics.py
│       ├── utils.py
│       ├── arc/
│       │   ├── __init__.py
│       │   ├── core.py             # ParameterAwareReservoir
│       │   ├── parameter_aware_rc.py
│       │   ├── trainer.py
│       │   ├── predictor.py
│       │   ├── results.py
│       │   └── exceptions.py
│       ├── systems/
│       │   ├── __init__.py
│       │   ├── ikeda.py
│       │   ├── food_chain.py
│       │   └── kuramoto_sivashinsky.py
│       ├── analysis/
│       │   ├── __init__.py
│       │   ├── critical_point.py
│       │   ├── transient_lifetime.py
│       │   └── ensemble.py
│       └── hpo/
│           ├── __init__.py
│           └── bayesian_opt.py
├── examples/
│   ├── ikeda_prediction.py
│   ├── food_chain_prediction.py
│   ├── food_chain_simulation.py
│   └── lorenz_prediction.py
├── tutorials/
│   └── rc_prediction_tutorial.ipynb
└── tests/
    ├── __init__.py
    ├── conftest.py
    ├── test_parameter_aware_rc.py
    ├── test_ikeda.py
    ├── test_food_chain.py
    ├── test_reservoirkit.py
    └── test_hyperparameter_tuning.py

Every subpackage (arc, systems, analysis, hpo) is a proper Python package with __init__.py.

---

API reference

Exported from rc_prediction (top level)

import rc_prediction
rc_prediction.__version__   # "0.2.0"

Category	Names
Main model	ParameterAwareRC
Core RC	ESN, Reservoir, RidgeReadout, BaseEstimator
Systems	IkedaMap, FoodChain, FoodChainParams, KuramotoSivashinsky
Results	ClosedLoopResult, ScanResult, EnsembleResult
Analysis	scan_critical_point, ensemble_predict, transient_lifetime
HPO	optimize_hyperparameters, HyperparameterSpace
Metrics	rmse, nrmse, mae, memory_capacity
Utils	lorenz_system, standardize, apply_standardize, add_noise, train_test_split_sequence
Exceptions	InsufficientParameterValues, ModelNotFittedError

lifetime_distribution is exported from rc_prediction.analysis only:

from rc_prediction.analysis import lifetime_distribution

Submodule exports

Submodule	Key exports
rc_prediction.arc	ParameterAwareRC, ClosedLoopResult, ScanResult, EnsembleResult, exceptions
rc_prediction.systems	IkedaMap, FoodChain, KuramotoSivashinsky, DEFAULT_TRAINING_MU, MU_CRITICAL, DEFAULT_TRAINING_ALPHA
rc_prediction.analysis	scan_critical_point, ensemble_predict, transient_lifetime, lifetime_distribution
rc_prediction.hpo	optimize_hyperparameters, HyperparameterSpace

System constants not re-exported at systems level (import from module):

from rc_prediction.systems.food_chain import DEFAULT_TRAINING_K, K_CRITICAL
from rc_prediction.systems.ikeda import DEFAULT_TRAINING_MU, MU_CRITICAL

Result dataclasses

ClosedLoopResult

Field	Type	Description
trajectory	ndarray	Closed-loop states, shape (n_steps, n_dims)
p_test	float	Test parameter
collapsed	bool	Whether collapse was detected
collapse_step	int | None	Step index of collapse
parameter_name	str	Name set in fit (default "p")

ScanResult

Field	Type	Description
p_values	ndarray	Scanned parameter grid
collapsed	ndarray	Boolean collapse flags
collapse_steps	ndarray	Collapse step per grid point ($-1$ if none)
p_critical	float	Estimated $p_c^{*}$
parameter_name	str	Parameter name

EnsembleResult

Field	Type	Description
p_test	float	Test parameter
collapsed_fraction	float	Fraction of realizations that collapsed
mean_lifetime	float	Mean transient lifetime
lifetimes	ndarray	Per-realization lifetimes
trajectories	list[ndarray]	Optional stored trajectories

ParameterAwareRC hyperparameters

Argument	Paper symbol	Default	Role
n_units	—	500	Reservoir size
average_degree	—	4.0	Mean degree of sparse $\mathbf{W}$
spectral_radius	$\rho(\mathbf{W})$	0.9	Spectral radius scaling
input_scaling	—	1.0	Scale of $\mathbf{W}_{in}$
param_gain	$k_b$	0.5	Parameter-channel gain
param_bias	$b_0$	0.0	Parameter-channel bias
leaking_rate	$\alpha$	0.3	Leak rate
ridge	$\lambda$	1e-8	Readout regularization
washout	—	10	Discarded initial training steps
random_state	—	None	RNG seed for weight init

optimize_hyperparameters in hpo/bayesian_opt.py performs random search over these ranges (the filename is historical; it is not full Gaussian-process Bayesian optimization).

---

Benchmark systems

Ikeda map (systems/ikeda.py)

Complex map with bifurcation parameter mu ($\mu_c = 1.0027$):

$$ z_{n+1} = \mu + \gamma z_n \exp\left(i\left(\kappa - \frac{\eta}{1+|z_n|^2}\right)\right) $$

Defaults: $\gamma=0.9$, $\kappa=0.4$, $\eta=6.0$.
simulate / simulate_training_set return shape (n_steps, 2) (real and imaginary parts).

Constant	Value
MU_CRITICAL	1.0027
DEFAULT_TRAINING_MU	(0.91, 0.94, 0.97)

Food chain (systems/food_chain.py)

Three-species Hastings–Powell / McCann–Yodzis model with carrying capacity K as bifurcation parameter ($K_c \approx 0.99976$).

State: resource R, consumer C, predator P — shape (n_steps, 3).

Constant	Value
K_CRITICAL	0.99976
DEFAULT_TRAINING_K	(0.97, 0.98, 0.99)

Post-critical simulations use warmup_K to settle on the pre-critical attractor first.

Kuramoto–Sivashinsky (systems/kuramoto_sivashinsky.py)

1D KS equation with bifurcation parameter alpha (paper supplementary material).
ETDRK4 spectral solver; output shape (n_steps, n_grid).

Constant	Value
DEFAULT_TRAINING_ALPHA	(196.0, 197.0, 198.0)
DEFAULT_N_GRID	32

---

Tutorial and examples

Notebook

pip install "rc_prediction[examples]"
jupyter notebook tutorials/rc_prediction_tutorial.ipynb

Sections: Ikeda collapse, food chain, ensemble prediction, transient lifetimes, optional HPO.
The first code cell adds src/ to sys.path when the package is not installed.

Scripts

python examples/ikeda_prediction.py
python examples/food_chain_prediction.py
python examples/food_chain_simulation.py
python examples/lorenz_prediction.py

---

Tests

pip install -e ".[dev]"
pytest

29 tests across five files. Use pytest, not python tests/test_*.py.

File	Coverage
test_parameter_aware_rc.py	ParameterAwareRC fit, predict, scan, ensemble
test_ikeda.py	Ikeda map, Kuramoto–Sivashinsky
test_food_chain.py	Food chain simulator and training set
test_reservoirkit.py	ESN, Reservoir, RidgeReadout, metrics
test_hyperparameter_tuning.py	optimize_hyperparameters

tests/conftest.py adds src/ to the import path and prints pass/fail summaries.

---

References

Primary method

1. Kong, L.-W., Fan, H.-W., Grebogi, C. and Lai, Y.-C. (2021) ‘Machine learning prediction of critical transition and system collapse’, Physical Review Research, 3(1), 013090. Available at: https://doi.org/10.1103/PhysRevResearch.3.013090

Reservoir computing

2. Jaeger, H. (2001) ‘The “echo state” approach to analysing and training recurrent neural networks’, GMD Report 148, German National Research Center for Information Technology. Available at: https://api.semanticscholar.org/CorpusID:15467150

Benchmark systems

3. Hastings, A. and Powell, T. (1991) ‘Chaos in a three-species food chain’, Ecology, 72(3), pp. 896–903. Available at: https://doi.org/10.2307/1940591.

4. McCann, K. and Yodzis, P. (1995) ‘Bifurcation structure of a three-species food-chain model’, Theoretical Population Biology, 48(2), pp. 93–125. Available at: https://doi.org/10.1006/tpbi.1995.1023.

5. Dhamala, M. and Lai, Y.-C. (1999) ‘Controlling transient chaos in deterministic flows with applications to electrical power systems and ecology’, Physical Review E, 59(2), pp. 1646–1655. Available at: https://doi.org/10.1103/PhysRevE.59.1646.

Related work

6. Panahi, S & Lai, YC 2024, 'Adaptable reservoir computing: A paradigm for model-free data-driven prediction of critical transitions in nonlinear dynamical systems', Chaos, vol. 34, no. 5, 051501. https://doi.org/10.1063/5.0200898

---

License

MIT — see LICENSE.