Metadata-Version: 2.4
Name: rc_prediction
Version: 0.1.1
Summary: Parameter-aware reservoir computing for critical transition and system collapse prediction
Project-URL: Homepage, https://github.com/jinchen7-cmd/Reservoir-Computing
Project-URL: Documentation, https://github.com/jinchen7-cmd/Reservoir-Computing#readme
Project-URL: Repository, https://github.com/jinchen7-cmd/Reservoir-Computing
Project-URL: Issues, https://github.com/jinchen7-cmd/Reservoir-Computing/issues
Author: Jincheng (Jeffery) Rao
License-Expression: MIT
License-File: LICENSE
Keywords: critical-transition,echo-state-network,esn,machine-learning,reservoir-computing,time-series,tipping-point,transient-chaos
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Topic :: Scientific/Engineering :: Artificial Intelligence
Requires-Python: >=3.10
Requires-Dist: numpy>=1.24
Requires-Dist: scikit-learn>=1.3
Requires-Dist: scipy>=1.10
Provides-Extra: dev
Requires-Dist: build>=1.0; extra == 'dev'
Requires-Dist: mypy>=1.8; extra == 'dev'
Requires-Dist: pytest-cov>=4.0; extra == 'dev'
Requires-Dist: pytest>=7.0; extra == 'dev'
Requires-Dist: ruff>=0.4; extra == 'dev'
Requires-Dist: twine>=5.0; extra == 'dev'
Provides-Extra: examples
Requires-Dist: ipykernel>=6.0; extra == 'examples'
Requires-Dist: jupyter>=1.0; extra == 'examples'
Requires-Dist: matplotlib>=3.7; extra == 'examples'
Description-Content-Type: text/markdown

# 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*](https://doi.org/10.1103/PhysRevResearch.3.013090), Phys. Rev. Research **3**, 013090 (2021).

---

## Contents

- [Overview](#overview)
- [Installation](#installation)
- [Quick start](#quick-start)
- [Method](#method)
- [Package structure](#package-structure)
- [API reference](#api-reference)
- [Benchmark systems](#benchmark-systems)
- [Tutorial and examples](#tutorial-and-examples)
- [Tests](#tests)
- [References](#references)
- [License](#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

```bash
pip install rc_prediction
```

Optional extras:

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

Editable install from source:

```bash
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

```python
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

```python
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

```python
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)

```python
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:

```python
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):

```python
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

```bash
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

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

---

## Tests

```bash
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](LICENSE).
