diff --git a/benchmarks/continuous_warmup/.gitignore b/benchmarks/continuous_warmup/.gitignore
new file mode 100644
index 000000000..acd10262f
--- /dev/null
+++ b/benchmarks/continuous_warmup/.gitignore
@@ -0,0 +1,7 @@
+Manifest.toml
+# Latest-run convenience copies (tagged files are committed when useful)
+results/summary.json
+results/summary.md
+# Keep FINDINGS.md + summary_full.* for the investigation PR
+*.png
+*.jld2
diff --git a/benchmarks/continuous_warmup/Project.toml b/benchmarks/continuous_warmup/Project.toml
new file mode 100644
index 000000000..8667f0a59
--- /dev/null
+++ b/benchmarks/continuous_warmup/Project.toml
@@ -0,0 +1,21 @@
+# Investigation environment only — not a Julia package.
+# Continuous autoregressive predict warmup (cold vs warm u0).
+
+[deps]
+DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
+JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+ReservoirComputing = "7c2d2b1e-3dd4-11ea-355a-8f6a8116e294"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[compat]
+DataInterpolations = "6, 7, 8, 9"
+JSON = "0.21, 1"
+OrdinaryDiffEq = "6"
+OrdinaryDiffEqTsit5 = "1"
+SciMLBase = "2, 3"
+julia = "1.10"
diff --git a/benchmarks/continuous_warmup/README.md b/benchmarks/continuous_warmup/README.md
new file mode 100644
index 000000000..68ba0550e
--- /dev/null
+++ b/benchmarks/continuous_warmup/README.md
@@ -0,0 +1,89 @@
+# Continuous AR predict warmup investigation
+
+**Branch:** `investigate/continuous-predict-warmup`  
+**Context:** Deferred from [#456](https://github.com/SciML/ReservoirComputing.jl/pull/456)  
+(Francesco: *“add the warmup as part of a separate PR… investigate that a bit further”*).
+
+This harness **does not modify package `src/` or root `Project.toml`**.  
+It only measures cold vs warm reservoir state at the start of continuous
+autoregressive `predict`, and probes candidate API shapes **locally**.
+
+## Problem
+
+| Path | Initial reservoir state for AR `predict` |
+|------|------------------------------------------|
+| Discrete `ESN` | Carries through `StatefulLayer` / sequential `apply` |
+| `SciMLProblemReservoir` | `res.prob.u0` (often zeros) |
+| `ContinuousESN` / `ContinuousESNCell` | Always `zeros(out_dims)` |
+
+On Lorenz, PR #456 measured roughly **NRMSE ~1.5 (cold)** vs **~0.11 (warmed)**.
+
+## Experiments
+
+| ID | Question |
+|----|----------|
+| E1 | Reproduce cold vs terminal-train warm on `ContinuousESN` (Lorenz) |
+| E2 | Same for hand-rolled `SciMLProblemReservoir` (eq. 5) |
+| E3 | Warmup **length** sweep (last `K` teacher-forced inputs before AR) |
+| E4 | Seed variants: zeros / random / train terminal / test-prefix warm / wrong state |
+| E5 | Horizon curve: NRMSE vs Lyapunov time (cold vs warm) |
+| E6 | Does `st` after `train!` / `collectstates` hold usable continuous state? |
+| E7 | Discrete `ESN` control (same data split) |
+| E8 | Washout interaction (train with washout, warm from post-washout tail) |
+
+## Setup
+
+```bash
+cd benchmarks/continuous_warmup
+julia --project=. -e 'using Pkg; Pkg.develop(path="../.."); Pkg.instantiate()'
+```
+
+## Run
+
+Smoke (plumbing only — small N, short series, ~20s):
+
+```bash
+julia --project=. run.jl --smoke
+```
+
+**Full matrix** (forecast analysis — matches #456 Lorenz scale):
+
+| knob | full default |
+|------|----------------|
+| `n_res` | 300 |
+| train / predict | 5000 / 1250 |
+| Wr spectral radius | 0.9 |
+| Win scale / bias scale | 0.1 / 0.05 |
+| solver | Tsit5 (`reltol=1e-6`) |
+
+```bash
+julia --project=. run.jl            # full
+julia --project=. run.jl --n-res=300
+julia --project=. run.jl --only=E1,E2,E5,E7
+```
+
+Expect minutes–tens of minutes for continuous `collectstates` at N=300 / T=5000 (wall times are recorded per experiment).
+
+Outputs land in `results/`:
+
+- `results/summary_full.json` / `summary_smoke.json` — machine-readable  
+- `results/summary_full.md` / `summary_smoke.md` — tables  
+- `results/summary.json` — latest run copy  
+- `results/FINDINGS.md` — interpretation
+
+## Design options under test (no API change yet)
+
+Local experimental helpers in `src/predict_variants.jl`:
+
+1. **`initial_state`** — pass reservoir `u0` into a mirrored AR loop  
+2. **`warmup_data`** — run teacher-forced `collectstates` on a prefix, take last column as `u0`, then AR  
+3. **`remake(prob; u0=…)`** — only for `SciMLProblemReservoir` public path  
+4. **Read `st` after train** — check whether continuous `st` already carries terminal state (expected: **no**)
+
+Package `predict` is left untouched so results stay comparable to merged `master`.
+
+## Relation to other work
+
+- Not PR4 perf ([#467](https://github.com/SciML/ReservoirComputing.jl/issues/467))  
+- Not a user-facing API PR until API choice is agreed  
+- Feeds a future tracking issue + design PR
diff --git a/benchmarks/continuous_warmup/results/FINDINGS.md b/benchmarks/continuous_warmup/results/FINDINGS.md
new file mode 100644
index 000000000..2336fa549
--- /dev/null
+++ b/benchmarks/continuous_warmup/results/FINDINGS.md
@@ -0,0 +1,176 @@
+# Findings — continuous AR predict warmup
+
+Branch: `investigate/continuous-predict-warmup`  
+PR: https://github.com/SciML/ReservoirComputing.jl/pull/476  
+Related: #456 (merged), #397  
+
+## Reproduce
+
+```bash
+cd benchmarks/continuous_warmup
+julia --project=. -e 'using Pkg; Pkg.develop(path="../.."); Pkg.instantiate()'
+
+# plumbing only (~20s)
+julia --project=. run.jl --smoke
+
+# forecast analysis — #456-scale (~5 min on M-series laptop)
+julia --project=. run.jl
+# → results/summary_full.{json,md}
+```
+
+### Full-run config
+
+| knob | value |
+|------|-------|
+| mode | **full** (not smoke) |
+| n_res | 300 |
+| train / predict | 5000 / 1250 samples |
+| dt / λ_max | 0.02 / 0.9056 → predict span ≈ **22.6 t_λ** |
+| Wr spectral radius | 0.9 |
+| Win scale / bias scale | 0.1 / 0.05 |
+| state modifiers | `NLAT2()` |
+| solver | Tsit5, `reltol=1e-6`, `abstol=1e-8` |
+| ridge | 1e-6 |
+| seed | 17 |
+| suite wall | **~275 s** (this machine) |
+
+---
+
+## 1. Structural (independent of HPs)
+
+| Check | Result |
+|-------|--------|
+| Seeded `u0=0` vs package `predict` | **max abs diff = 0** |
+| `remake(prob; u0=…)` vs seeded warm | **identical NRMSE** |
+| Continuous `st` after `train!` / `collectstates` | **no carry** |
+| Discrete `st` after `train!` | **has carry** |
+| Discrete rewarm via `collectstates` | ≡ post-train `st` |
+
+**Conclusion:** continuous AR always cold-starts today. There is nowhere in `st` to put a terminal ODE state. Public warm workaround for generic SciML reservoirs: `remake(prob; u0=terminal)`. `ContinuousESN` has no public `u0` hook without a custom AR loop.
+
+---
+
+## 2. Forecast quality — full matrix (the real analysis)
+
+### Headline: short-horizon NRMSE (ContinuousESN, E5 / E1 horizons)
+
+| Horizon | steps | t_λ | cold NRMSE | warm (train-terminal) NRMSE |
+|---------|------:|----:|-----------:|----------------------------:|
+| short | 28 | 0.5 | 0.75 | **0.21** |
+| 1 t_λ | 55 | 1.0 | 0.69 | **0.15** |
+| 2 t_λ | 110 | 2.0 | 1.34 | **0.085** |
+| 3 t_λ | 166 | 3.0 | 1.51 | **0.093** |
+| 4 t_λ | 221 | 4.0 | 1.66 | **0.14** |
+| 6 t_λ | 331 | 6.0 | 1.51 | **0.67** |
+
+Warm recovers the #456-class short-horizon gap: at 2–3 Lyapunov times, warm is ~**0.09** vs cold ~**1.3–1.5**.
+
+Full-horizon NRMSE over all 1250 steps (≈22 t_λ) is high for both (cold **1.48**, warm **1.23**) — expected once chaos has diverged; the interesting regime is the first few t_λ.
+
+### Valid prediction time (E3, threshold 0.5)
+
+| Warmup K | VPT (t_λ) | full-horizon NRMSE |
+|---------:|----------:|-------------------:|
+| 0 (cold) | **0.22** | 1.48 |
+| 10 | **4.87** | 1.14 |
+| 50–2000 | ~4.13 | ~1.18–1.28 |
+
+Even a short teacher-forced warmup (K=10) lifts VPT from ~0.2 → ~5 t_λ.
+
+### Seed ranking (E4, full-horizon NRMSE)
+
+| seed | NRMSE |
+|------|------:|
+| train_terminal | **1.23** |
+| zeros (cold) | 1.48 |
+| oracle test prefix (K=100) | 1.45 |
+| randn / shuffled terminal | ~13.5 (catastrophic) |
+
+Wrong non-zero seeds are much worse than cold. Warm is not “any u0” — it has to be a dynamically consistent terminal state.
+
+### SciMLProblemReservoir eq.5 (E2)
+
+| variant | full-horizon NRMSE |
+|---------|-------------------:|
+| cold | 1.47 |
+| warm train-terminal | **1.31** |
+| `remake(prob; u0=…)` | **1.31** (matches seeded) |
+
+### Discrete control (E7)
+
+| variant | full-horizon NRMSE |
+|---------|-------------------:|
+| cold (fresh `st`) | 1.55 |
+| warm (post-train `st`) | **1.14** |
+| rewarm via `collectstates` | **1.14** |
+
+### Washout (E8, washout=200)
+
+| variant | NRMSE |
+|---------|------:|
+| cold | 1.48 |
+| warm full train | 1.18 |
+| warm post-washout tail | **1.13** |
+
+Post-washout tail is slightly better than full-train terminal.
+
+---
+
+## 3. Wall times (full, not a perf PR — order-of-magnitude only)
+
+| stage | wall |
+|-------|------|
+| ContinuousESN train (`collectstates` + ridge), N=300, T=5000 | ~18 s |
+| AR cold predict, 1250 steps | ~2 s |
+| AR warm predict, 1250 steps | ~1 s |
+| Warmup collect (full train pass for terminal u0) | ~train-scale |
+| Full E1–E8 suite | ~275 s |
+
+Warm vs cold AR cost is comparable; the extra cost of warm is mainly the one-time teacher-forced collect for `u0`. Not a throughput regression story — a correctness/usability one. Perf work stays on #467.
+
+---
+
+## 4. Smoke vs full (why smoke lied)
+
+| | smoke (n=80, T=800) | full (n=300, T=5000) |
+|--|---------------------|----------------------|
+| ContinuousESN warm vs cold | warm **hurt** | warm **helps** (esp. short horizon) |
+| SciML warm vs cold | helps a little | helps |
+| Structural checks | same | same |
+
+Smoke is fine for plumbing (`match_ok`, carry inspection). **Do not** use smoke for API justification.
+
+---
+
+## 5. API recommendation (updated after full run)
+
+| Option | Verdict after full data |
+|--------|-------------------------|
+| `predict(...; initial_state=u0)` | **Yes — primitive.** Matches remake; E4 shows wrong seeds are dangerous so it should be explicit. |
+| `predict(...; warmup_data=W)` | **Yes — sugar.** K≥10 already unlocks ~5 t_λ VPT; must use unit windows (`Δt=1`). |
+| Persist terminal in `st` | Still attractive for discrete-like ergonomics; needs continuous `st` design — separate decision |
+| Docs-only | Insufficient: cold looks like a broken model at 1–4 t_λ |
+
+**Proposed first code PR:**  
+`initial_state` on continuous AR `predict` + unit tests that seeded zeros ≡ today’s cold path and that a known `u0` is actually used. Optional `warmup_data` in the same PR or immediately after. Lorenz short-horizon eye-test optional but strong.
+
+---
+
+## 6. Questions for @MartinuzziFrancesco
+
+1. OK to treat short-horizon (1–4 t_λ) NRMSE / VPT as the acceptance metric for warmup, not full-horizon NRMSE?
+2. Ship `initial_state` first, or `warmup_data` first, or both?
+3. Should `train!` / `collectstates` start writing a continuous terminal into `st` (discrete-like), or keep state external via kwargs?
+4. Keep this harness under `benchmarks/` long-term?
+
+---
+
+## File map
+
+| Path | Role |
+|------|------|
+| `run.jl` | `--smoke` / full / `--only=` / `--n-res=` |
+| `src/predict_variants.jl` | experimental seeded AR (no package API) |
+| `src/experiments.jl` | E1–E8 |
+| `results/summary_full.md` | latest full table |
+| `results/FINDINGS.md` | this document |
diff --git a/benchmarks/continuous_warmup/results/summary_full.json b/benchmarks/continuous_warmup/results/summary_full.json
new file mode 100644
index 000000000..d2dc2c2c0
--- /dev/null
+++ b/benchmarks/continuous_warmup/results/summary_full.json
@@ -0,0 +1,625 @@
+[
+  {
+    "experiment": "E1",
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res": 300,
+    "n_res_cfg": 300,
+    "nrmse": 1.4782807405733622,
+    "nrmse_global": 0.8686309038147971,
+    "predict_len": 1250,
+    "timestamp": "2026-07-18T02:19:39.220",
+    "train_len": 5000,
+    "variant": "cold_package_predict",
+    "vpt_lyap": 0.21734399999999998,
+    "wall_ar_s": 2.060436125,
+    "wall_train_s": 18.117923542
+  },
+  {
+    "experiment": "E1",
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res": 300,
+    "n_res_cfg": 300,
+    "nrmse": 1.2292196160902893,
+    "nrmse_global": 0.7208150983122181,
+    "predict_len": 1250,
+    "timestamp": "2026-07-18T02:19:39.229",
+    "train_len": 5000,
+    "u0_norm": 13.724870901821692,
+    "variant": "warm_train_terminal_u0",
+    "vpt_lyap": 4.129536,
+    "wall_ar_s": 0.952867458,
+    "wall_train_s": 18.117923542,
+    "wall_warmup_collect_s": 16.314653625
+  },
+  {
+    "experiment": "E1",
+    "match_ok": true,
+    "max_abs_diff": 0.0,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res": 300,
+    "n_res_cfg": 300,
+    "timestamp": "2026-07-18T02:19:39.237",
+    "variant": "seeded_zero_matches_cold_maxabs"
+  },
+  {
+    "experiment": "E1",
+    "horizon_lyap": 0.507136,
+    "horizon_steps": 28,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 0.7535572065523274,
+    "timestamp": "2026-07-18T02:19:39.372",
+    "variant": "cold_horizon"
+  },
+  {
+    "experiment": "E1",
+    "horizon_lyap": 0.507136,
+    "horizon_steps": 28,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 0.21463110240826375,
+    "timestamp": "2026-07-18T02:19:39.372",
+    "variant": "warm_horizon"
+  },
+  {
+    "experiment": "E1",
+    "horizon_lyap": 0.99616,
+    "horizon_steps": 55,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 0.6853049713043267,
+    "timestamp": "2026-07-18T02:19:39.372",
+    "variant": "cold_horizon"
+  },
+  {
+    "experiment": "E1",
+    "horizon_lyap": 0.99616,
+    "horizon_steps": 55,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 0.14675113554137212,
+    "timestamp": "2026-07-18T02:19:39.372",
+    "variant": "warm_horizon"
+  },
+  {
+    "experiment": "E1",
+    "horizon_lyap": 1.99232,
+    "horizon_steps": 110,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.3393780477247423,
+    "timestamp": "2026-07-18T02:19:39.372",
+    "variant": "cold_horizon"
+  },
+  {
+    "experiment": "E1",
+    "horizon_lyap": 1.99232,
+    "horizon_steps": 110,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 0.08497214024585852,
+    "timestamp": "2026-07-18T02:19:39.372",
+    "variant": "warm_horizon"
+  },
+  {
+    "experiment": "E1",
+    "horizon_lyap": 3.006592,
+    "horizon_steps": 166,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.506054216154576,
+    "timestamp": "2026-07-18T02:19:39.372",
+    "variant": "cold_horizon"
+  },
+  {
+    "experiment": "E1",
+    "horizon_lyap": 3.006592,
+    "horizon_steps": 166,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 0.093388739272057,
+    "timestamp": "2026-07-18T02:19:39.372",
+    "variant": "warm_horizon"
+  },
+  {
+    "experiment": "E1",
+    "horizon_lyap": 4.002752,
+    "horizon_steps": 221,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.6554797621233028,
+    "timestamp": "2026-07-18T02:19:39.372",
+    "variant": "cold_horizon"
+  },
+  {
+    "experiment": "E1",
+    "horizon_lyap": 4.002752,
+    "horizon_steps": 221,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 0.13715717405603614,
+    "timestamp": "2026-07-18T02:19:39.372",
+    "variant": "warm_horizon"
+  },
+  {
+    "experiment": "E1",
+    "horizon_lyap": 5.9950719999999995,
+    "horizon_steps": 331,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.5117492481722878,
+    "timestamp": "2026-07-18T02:19:39.372",
+    "variant": "cold_horizon"
+  },
+  {
+    "experiment": "E1",
+    "horizon_lyap": 5.9950719999999995,
+    "horizon_steps": 331,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 0.6692777566540231,
+    "timestamp": "2026-07-18T02:19:39.372",
+    "variant": "warm_horizon"
+  },
+  {
+    "experiment": "E2",
+    "mode": "full",
+    "model": "SciMLProblemReservoir",
+    "n_res": 300,
+    "n_res_cfg": 300,
+    "nrmse": 1.470005354266804,
+    "nrmse_global": 0.8610435641262816,
+    "timestamp": "2026-07-18T02:20:19.420",
+    "variant": "cold_package_predict",
+    "vpt_lyap": 0.018112
+  },
+  {
+    "experiment": "E2",
+    "mode": "full",
+    "model": "SciMLProblemReservoir",
+    "n_res": 300,
+    "n_res_cfg": 300,
+    "nrmse": 1.312672661950845,
+    "nrmse_global": 0.7709100283139908,
+    "timestamp": "2026-07-18T02:20:19.429",
+    "u0_norm": 13.604509626860247,
+    "variant": "warm_train_terminal_u0",
+    "vpt_lyap": 2.843584
+  },
+  {
+    "experiment": "E2",
+    "max_abs_diff_vs_seeded": 0.0,
+    "mode": "full",
+    "model": "SciMLProblemReservoir",
+    "n_res": 300,
+    "n_res_cfg": 300,
+    "nrmse": 1.312672661950845,
+    "nrmse_global": 0.7709100283139908,
+    "timestamp": "2026-07-18T02:20:19.437",
+    "variant": "remake_prob_u0_package_predict",
+    "vpt_lyap": 2.843584
+  },
+  {
+    "experiment": "E3",
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.4782807405733622,
+    "nrmse_global": 0.8686309038147971,
+    "timestamp": "2026-07-18T02:20:36.899",
+    "variant": "K=0_cold",
+    "vpt_lyap": 0.21734399999999998,
+    "warmup_len": 0
+  },
+  {
+    "experiment": "E3",
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.1420303622679457,
+    "nrmse_global": 0.6695443570238542,
+    "timestamp": "2026-07-18T02:20:37.882",
+    "variant": "K=10",
+    "vpt_lyap": 4.872128,
+    "warmup_len": 10
+  },
+  {
+    "experiment": "E3",
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.272365316049519,
+    "nrmse_global": 0.7447275884595839,
+    "timestamp": "2026-07-18T02:20:38.974",
+    "variant": "K=50",
+    "vpt_lyap": 4.129536,
+    "warmup_len": 50
+  },
+  {
+    "experiment": "E3",
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.2309111911851687,
+    "nrmse_global": 0.7277626070587372,
+    "timestamp": "2026-07-18T02:20:40.254",
+    "variant": "K=100",
+    "vpt_lyap": 4.129536,
+    "warmup_len": 100
+  },
+  {
+    "experiment": "E3",
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.2163625936067015,
+    "nrmse_global": 0.7141959888660526,
+    "timestamp": "2026-07-18T02:20:42.025",
+    "variant": "K=250",
+    "vpt_lyap": 4.129536,
+    "warmup_len": 250
+  },
+  {
+    "experiment": "E3",
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.2843248580360316,
+    "nrmse_global": 0.7517730578199671,
+    "timestamp": "2026-07-18T02:20:44.597",
+    "variant": "K=500",
+    "vpt_lyap": 4.129536,
+    "warmup_len": 500
+  },
+  {
+    "experiment": "E3",
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.1789022279431667,
+    "nrmse_global": 0.6993686867334836,
+    "timestamp": "2026-07-18T02:20:48.850",
+    "variant": "K=1000",
+    "vpt_lyap": 4.129536,
+    "warmup_len": 1000
+  },
+  {
+    "experiment": "E3",
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.1835926940405561,
+    "nrmse_global": 0.692584966096159,
+    "timestamp": "2026-07-18T02:20:56.329",
+    "variant": "K=2000",
+    "vpt_lyap": 4.129536,
+    "warmup_len": 2000
+  },
+  {
+    "experiment": "E4",
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 13.54940655959303,
+    "nrmse_global": 9.826713305822715,
+    "timestamp": "2026-07-18T02:21:29.835",
+    "u0_norm": 13.724870901821683,
+    "variant": "shuffled_train_terminal",
+    "vpt_lyap": 0.0
+  },
+  {
+    "experiment": "E4",
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 13.558170985119219,
+    "nrmse_global": 9.833322066308458,
+    "timestamp": "2026-07-18T02:21:30.115",
+    "u0_norm": 17.851901593047117,
+    "variant": "randn",
+    "vpt_lyap": 0.0
+  },
+  {
+    "experiment": "E4",
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.4782807405733622,
+    "nrmse_global": 0.8686309038147971,
+    "timestamp": "2026-07-18T02:21:31.073",
+    "u0_norm": 0.0,
+    "variant": "zeros",
+    "vpt_lyap": 0.21734399999999998
+  },
+  {
+    "experiment": "E4",
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.2292196160902893,
+    "nrmse_global": 0.7208150983122181,
+    "timestamp": "2026-07-18T02:21:32.024",
+    "u0_norm": 13.724870901821692,
+    "variant": "train_terminal",
+    "vpt_lyap": 4.129536
+  },
+  {
+    "experiment": "E4",
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.4509837707185447,
+    "nrmse_global": 0.8508405432103581,
+    "timestamp": "2026-07-18T02:21:32.968",
+    "u0_norm": 15.374881292468654,
+    "variant": "oracle_test_prefix",
+    "vpt_lyap": 0.959936
+  },
+  {
+    "experiment": "E5",
+    "horizon_lyap": 0.507136,
+    "horizon_steps": 28,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 0.7535572065523274,
+    "timestamp": "2026-07-18T02:22:07.810",
+    "variant": "cold"
+  },
+  {
+    "experiment": "E5",
+    "horizon_lyap": 0.507136,
+    "horizon_steps": 28,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 0.21463110240826375,
+    "timestamp": "2026-07-18T02:22:07.810",
+    "variant": "warm_train_terminal"
+  },
+  {
+    "experiment": "E5",
+    "horizon_lyap": 0.99616,
+    "horizon_steps": 55,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 0.6853049713043267,
+    "timestamp": "2026-07-18T02:22:07.821",
+    "variant": "cold"
+  },
+  {
+    "experiment": "E5",
+    "horizon_lyap": 0.99616,
+    "horizon_steps": 55,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 0.14675113554137212,
+    "timestamp": "2026-07-18T02:22:07.821",
+    "variant": "warm_train_terminal"
+  },
+  {
+    "experiment": "E5",
+    "horizon_lyap": 1.99232,
+    "horizon_steps": 110,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.3393780477247423,
+    "timestamp": "2026-07-18T02:22:07.821",
+    "variant": "cold"
+  },
+  {
+    "experiment": "E5",
+    "horizon_lyap": 1.99232,
+    "horizon_steps": 110,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 0.08497214024585852,
+    "timestamp": "2026-07-18T02:22:07.821",
+    "variant": "warm_train_terminal"
+  },
+  {
+    "experiment": "E5",
+    "horizon_lyap": 3.006592,
+    "horizon_steps": 166,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.506054216154576,
+    "timestamp": "2026-07-18T02:22:07.821",
+    "variant": "cold"
+  },
+  {
+    "experiment": "E5",
+    "horizon_lyap": 3.006592,
+    "horizon_steps": 166,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 0.093388739272057,
+    "timestamp": "2026-07-18T02:22:07.821",
+    "variant": "warm_train_terminal"
+  },
+  {
+    "experiment": "E5",
+    "horizon_lyap": 4.002752,
+    "horizon_steps": 221,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.6554797621233028,
+    "timestamp": "2026-07-18T02:22:07.821",
+    "variant": "cold"
+  },
+  {
+    "experiment": "E5",
+    "horizon_lyap": 4.002752,
+    "horizon_steps": 221,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 0.13715717405603614,
+    "timestamp": "2026-07-18T02:22:07.821",
+    "variant": "warm_train_terminal"
+  },
+  {
+    "experiment": "E5",
+    "horizon_lyap": 5.9950719999999995,
+    "horizon_steps": 331,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.5117492481722878,
+    "timestamp": "2026-07-18T02:22:07.821",
+    "variant": "cold"
+  },
+  {
+    "experiment": "E5",
+    "horizon_lyap": 5.9950719999999995,
+    "horizon_steps": 331,
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 0.6692777566540231,
+    "timestamp": "2026-07-18T02:22:07.821",
+    "variant": "warm_train_terminal"
+  },
+  {
+    "experiment": "E6",
+    "has_carry": false,
+    "keys": {},
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "st_type": "@NamedTuple{}",
+    "timestamp": "2026-07-18T02:22:41.403",
+    "variant": "continuous_st0"
+  },
+  {
+    "experiment": "E6",
+    "has_carry": false,
+    "keys": {},
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "st_type": "@NamedTuple{}",
+    "summary": "NamedTuple()",
+    "timestamp": "2026-07-18T02:22:41.412",
+    "variant": "continuous_st_after_train"
+  },
+  {
+    "experiment": "E6",
+    "has_carry": false,
+    "keys": {},
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "st_type": "@NamedTuple{}",
+    "summary": "NamedTuple()",
+    "timestamp": "2026-07-18T02:22:41.412",
+    "variant": "continuous_st_after_collectstates"
+  },
+  {
+    "experiment": "E6",
+    "has_carry": true,
+    "keys": [
+      "cell",
+      "carry"
+    ],
+    "mode": "full",
+    "model": "ESN",
+    "n_res_cfg": 300,
+    "st_type": "@NamedTuple{cell::@NamedTuple{rng::MersenneTwister}, carry::Tuple{Matrix{Float64}}}",
+    "summary": "(cell = (rng = MersenneTwister(17, (0, 222606, 221156, 931)),), carry = ([-0.8459772455802655; -0.48366492020831586; -0.8995900043994268; 0.9292338400401382; 0.7345940463048126; 0.842570735371275; 0.1602216455975732; 0.796943209513598; -0.8512128203797109; 0.046498299310999684; 0.8838683222507832; -0.9223242369221242; -0.6519113630471485; -0.9505135159408366; -0.5897116160716003; 0.5127238699590044; 0.5354341653853366; -0.2167309691344319; 0.985432447331775; 0.9344514907704993; 0.8145426255832817; 0.9051898459753677; 0.8962051796538499; 0.137759426856421; -0.7519732623630886; -0.20909311108898415; 0.9812807430448122; 0.6410289905641998; -0.6728006780332255; 0.664373535093157; 0.9766478047976146; 0.914544128529556; 0.27372974609337847; -0.9559208773872372; -0.9426315592580121; 0.6328265155719525; -0.25541348400825364; 0.31589005395925657; 0.5874013583401602; 0.6302850655211003; -0.883167140502558; -0.6756226790046802; 0.22586646200224147; -0.3033582102461979; 0.9477418131446771; 0.767710699638399; 0.4201556162825191; -0.9758012319887187; 0.5844813589539288; 0.648496264993395; 0.9609372370844307; 0.990302950955413; 0.846382953995659; -0.42228945910258475; 0.00523671371802411; -0.9733695922037601; -0.820489249569788; -0.6535195035919952; 0.6675207884781951; -0.7028988827929498; 0.5556076878653934; -0.6715772240494514; 0.6947637693463453; -0.9899517772404052; 0.17144218006467082; 0.4077795504510223; -0.652425209957331; -0.9419907702395707; -0.4519756441618818; -0.7038513519182468; 0.4977929719966706; -0.7014394265798016; 0.7541835358747554; 0.9301383003521996; -0.5362117226565815; 0.18873336485259964; 0.8875076536413823; 0.9822517972762308; 0.669011868522317; 0.07788935520668964; -0.04441664047121452; 0.6510802194786907; -0.5812102464471012; -0.8260085665152631; -0.9848503272424689; -0.41692469602444604; -0.7248076857310832; -0.22042986226192968; 0.40678548336583736; -0.11686465193849968; 0.6211276960123986; 0.8632711135608314; 0.6798899693921168; 0.9417805911702982; 0.8284028339230047; -0.8632014606518237; 0.7010704444408193; -0.19409747469218605; 0.32605845995723015; 0.11137635744171862; 0.7232349295110461; 0.5288191261418995; -0.8602099147912512; 0.6659429370642505; 0.29856432492754276; -0.4132264031475149; 0.2714976646117236; -0.9506295300151116; 0.9945690440130633; 0.09487869328434272; 0.6364383039092928; -0.3100450897775931; -0.44618462491572336; -0.9205760393779723; -0.5817518095393643; -0.3569748261011347; -0.12139133175997965; 0.8037471454674672; 0.6205796968890384; -0.13333512145187593; 0.03277340603932864; 0.6421086675274126; 0.4423466831147324; -0.7321815375377059; 0.4795760808984699; 0.8820010330669477; -0.4826418160753365; -0.9548985854368519; -0.9440353040153748; 0.8578351379955963; -0.06160592298640819; 0.6432160612960786; 0.45402726603339727; -0.09758076073355423; 0.897917733227205; 0.9647144656540431; -0.890701400588892; 0.22978498293575111; 0.6193533139534437; -0.882104563130712; -0.7983985116140266; -0.5595327390757587; -0.8280436905363562; -0.5323732914443493; -0.6265092358657431; 0.8106471165139244; -0.8641473730933273; -0.5353575698176427; 0.9624022974474143; -0.06989477178555709; -0.7738741988302407; 0.3714708206847453; 0.28149824861190464; -0.9692691302592675; -0.1375259877946844; -0.6610357991578912; 0.8146785524843005; -0.04658208591897321; 0.9172088122035855; 0.9327247467448289; -0.946682514426566; 0.7938365153278566; -0.23760288054436338; -0.530511316355213; 0.19048954905470655; -0.9965366112635179; 0.21895259034872067; 0.9524371308203091; -0.44910563641119483; -0.8633037434084547; 0.9262871500728647; -0.2730127550963002; -0.8129254360492008; 0.2970538365982147; -0.7222122162909113; 0.1740655153073884; 0.4905189458865581; -0.17621971001845615; -0.9362712422178356; 0.2873605816875309; -0.2152210876610816; -0.8163408836427131; -0.625571435144257; 0.3338799038713975; 0.9693167693763403; -0.9464384305050012; 0.5527996802044771; 0.6666161741354575; -0.9875056025662363; 0.2933523329367735; 0.7945637952096477; 0.6848699298402934; 0.6490173585859442; 0.7121880769822362; -0.2050215316213654; -0.9464917146560637; -0.9344099699156242; -0.8547797484659176; 0.9259890956140381; 0.7738671507833101; 0.9642585389306354; -0.8186651425910476; 0.7689150170866851; -0.01517593411380655; 0.9029068492878274; 0.7700865111798508; 0.9391347492348087; 0.7035858495850119; 0.22658454521285898; -0.8945644915691096; 0.8981745857254391; -0.025559124195450252; -0.8745995559836761; -0.3085204916509622; -0.7356657931574957; 0.5211330987819445; 0.6450410140833305; 0.8525806634644354; -0.46619714775539356; -0.8590800702345783; -0.7286830676134923; -0.20594607384239044; 0.791115315345083; -0.11052696748761318; 0.9486810127502264; 0.4136600811193309; 0.37265775729587985; -0.05600314969338066; -0.973570671173324; -0.671441985829284; -0.45569634189587643; -0.7211615609080226; 0.014704909810522037; 0.8048112315075959; -0.6481325205587654; 0.9534192047643347; -0.8930635000865577; 0.4084394383801676; -0.8000347487930202; -0.9073078266439747; -0.4997137739097308; 0.7339105813275539; -0.8472063620600302; -0.8386002726899132; -0.7665100536580777; 0.8931869354544422; 0.9323123145306803; 0.7163078260624565; -0.16517049635126113; -0.7498589887611965; -0.8172868254049097; -0.9270744867029621; -0.9148573034409024; 0.6942638232334406; 0.800730588744462; -0.5815811636183894; -0.8783919868476993; 0.5764684951965523; 0.4277939916454236; 0.9850259752563585; 0.31314596679856016; -0.011097033801206316; -0.19022005049022067; 0.8047063192550168; 0.9981176614360656; 0.9907908065742174; -0.312655175073292; 0.7181730407132259; -0.4798611404298997; 0.9145611970952096; 0.2494115455251776; -0.7500082441544813; 0.8741729659723434; -0.15941779767815806; -0.857814506950908; -0.8377440624455936; -0.39018199425008676; 0.8607088655003254; 0.21156308876872038; -0.5390233551686354; 0.973741803980379; -0.4872638494293887; 0.7985475451363055; -0.7275849054361054; -0.8663567279481859; -0.4106374635999801; -0.9668857343152594; -0.6291793682394918; -0.7743726475162771; -0.47922129105278854; -0.9790459410762601; -0.3139641883961854; 0.8548847063182738; -0.02583371326313954; 0.5138527998796617; -0.7600583280671034; -0.8977722082061435; 0.8071661649346491; -0.36194127746207366; 0.013380833786579114;;],))",
+    "timestamp": "2026-07-18T02:22:41.421",
+    "variant": "discrete_st_after_train"
+  },
+  {
+    "experiment": "E7",
+    "mode": "full",
+    "model": "ESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.55491741812422,
+    "timestamp": "2026-07-18T02:22:42.181",
+    "variant": "cold_fresh_st",
+    "vpt_lyap": 0.0
+  },
+  {
+    "experiment": "E7",
+    "mode": "full",
+    "model": "ESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.1363799866710107,
+    "timestamp": "2026-07-18T02:22:42.181",
+    "variant": "warm_st_after_train",
+    "vpt_lyap": 4.1114239999999995
+  },
+  {
+    "experiment": "E7",
+    "mode": "full",
+    "model": "ESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.1363799866710107,
+    "timestamp": "2026-07-18T02:22:42.181",
+    "variant": "rewarm_collectstates_then_ar",
+    "vpt_lyap": 4.1114239999999995
+  },
+  {
+    "experiment": "E8",
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.4840906014114033,
+    "timestamp": "2026-07-18T02:23:34.772",
+    "variant": "cold",
+    "washout": 200
+  },
+  {
+    "experiment": "E8",
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.1839060794299072,
+    "timestamp": "2026-07-18T02:23:34.772",
+    "variant": "warm_full_train",
+    "washout": 200
+  },
+  {
+    "experiment": "E8",
+    "mode": "full",
+    "model": "ContinuousESN",
+    "n_res_cfg": 300,
+    "nrmse": 1.1280934497932662,
+    "timestamp": "2026-07-18T02:23:34.772",
+    "variant": "warm_post_washout_tail",
+    "washout": 200
+  },
+  {
+    "experiment": "META",
+    "mode": "full",
+    "model": "harness",
+    "n_res": 300,
+    "predict_len": 1250,
+    "train_len": 5000,
+    "variant": "suite_wall_s",
+    "wall_s": 275.113651291
+  }
+]
\ No newline at end of file
diff --git a/benchmarks/continuous_warmup/results/summary_full.md b/benchmarks/continuous_warmup/results/summary_full.md
new file mode 100644
index 000000000..6e3869bdd
--- /dev/null
+++ b/benchmarks/continuous_warmup/results/summary_full.md
@@ -0,0 +1,60 @@
+# Continuous warmup investigation results
+
+Generated: 2026-07-18T02:23:35.043
+
+| exp | variant | model | nrmse | nrmse_global | vpt_lyap | notes |
+|-----|---------|-------|-------|--------------|----------|-------|
+| E1 | cold_package_predict | ContinuousESN | 1.4782807405733622 | 0.8686309038147971 | 0.21734399999999998 | wall_train_s=18.117923542; wall_ar_s=2.060436125; mode=full |
+| E1 | warm_train_terminal_u0 | ContinuousESN | 1.2292196160902893 | 0.7208150983122181 | 4.129536 | u0_norm=13.724870901821692; wall_train_s=18.117923542; wall_ar_s=0.952867458; wall_warmup_collect_s=16.314653625; mode=full |
+| E1 | seeded_zero_matches_cold_maxabs | ContinuousESN |  |  |  | max_abs_diff=0.0; match_ok=true; mode=full |
+| E1 | cold_horizon | ContinuousESN | 0.7535572065523274 |  |  | horizon_steps=28; horizon_lyap=0.507136; mode=full |
+| E1 | warm_horizon | ContinuousESN | 0.21463110240826375 |  |  | horizon_steps=28; horizon_lyap=0.507136; mode=full |
+| E1 | cold_horizon | ContinuousESN | 0.6853049713043267 |  |  | horizon_steps=55; horizon_lyap=0.99616; mode=full |
+| E1 | warm_horizon | ContinuousESN | 0.14675113554137212 |  |  | horizon_steps=55; horizon_lyap=0.99616; mode=full |
+| E1 | cold_horizon | ContinuousESN | 1.3393780477247423 |  |  | horizon_steps=110; horizon_lyap=1.99232; mode=full |
+| E1 | warm_horizon | ContinuousESN | 0.08497214024585852 |  |  | horizon_steps=110; horizon_lyap=1.99232; mode=full |
+| E1 | cold_horizon | ContinuousESN | 1.506054216154576 |  |  | horizon_steps=166; horizon_lyap=3.006592; mode=full |
+| E1 | warm_horizon | ContinuousESN | 0.093388739272057 |  |  | horizon_steps=166; horizon_lyap=3.006592; mode=full |
+| E1 | cold_horizon | ContinuousESN | 1.6554797621233028 |  |  | horizon_steps=221; horizon_lyap=4.002752; mode=full |
+| E1 | warm_horizon | ContinuousESN | 0.13715717405603614 |  |  | horizon_steps=221; horizon_lyap=4.002752; mode=full |
+| E1 | cold_horizon | ContinuousESN | 1.5117492481722878 |  |  | horizon_steps=331; horizon_lyap=5.9950719999999995; mode=full |
+| E1 | warm_horizon | ContinuousESN | 0.6692777566540231 |  |  | horizon_steps=331; horizon_lyap=5.9950719999999995; mode=full |
+| E2 | cold_package_predict | SciMLProblemReservoir | 1.470005354266804 | 0.8610435641262816 | 0.018112 | mode=full |
+| E2 | warm_train_terminal_u0 | SciMLProblemReservoir | 1.312672661950845 | 0.7709100283139908 | 2.843584 | u0_norm=13.604509626860247; mode=full |
+| E2 | remake_prob_u0_package_predict | SciMLProblemReservoir | 1.312672661950845 | 0.7709100283139908 | 2.843584 | mode=full |
+| E3 | K=0_cold | ContinuousESN | 1.4782807405733622 | 0.8686309038147971 | 0.21734399999999998 | warmup_len=0; mode=full |
+| E3 | K=10 | ContinuousESN | 1.1420303622679457 | 0.6695443570238542 | 4.872128 | warmup_len=10; mode=full |
+| E3 | K=50 | ContinuousESN | 1.272365316049519 | 0.7447275884595839 | 4.129536 | warmup_len=50; mode=full |
+| E3 | K=100 | ContinuousESN | 1.2309111911851687 | 0.7277626070587372 | 4.129536 | warmup_len=100; mode=full |
+| E3 | K=250 | ContinuousESN | 1.2163625936067015 | 0.7141959888660526 | 4.129536 | warmup_len=250; mode=full |
+| E3 | K=500 | ContinuousESN | 1.2843248580360316 | 0.7517730578199671 | 4.129536 | warmup_len=500; mode=full |
+| E3 | K=1000 | ContinuousESN | 1.1789022279431667 | 0.6993686867334836 | 4.129536 | warmup_len=1000; mode=full |
+| E3 | K=2000 | ContinuousESN | 1.1835926940405561 | 0.692584966096159 | 4.129536 | warmup_len=2000; mode=full |
+| E4 | shuffled_train_terminal | ContinuousESN | 13.54940655959303 | 9.826713305822715 | 0.0 | u0_norm=13.724870901821683; mode=full |
+| E4 | randn | ContinuousESN | 13.558170985119219 | 9.833322066308458 | 0.0 | u0_norm=17.851901593047117; mode=full |
+| E4 | zeros | ContinuousESN | 1.4782807405733622 | 0.8686309038147971 | 0.21734399999999998 | u0_norm=0.0; mode=full |
+| E4 | train_terminal | ContinuousESN | 1.2292196160902893 | 0.7208150983122181 | 4.129536 | u0_norm=13.724870901821692; mode=full |
+| E4 | oracle_test_prefix | ContinuousESN | 1.4509837707185447 | 0.8508405432103581 | 0.959936 | u0_norm=15.374881292468654; mode=full |
+| E5 | cold | ContinuousESN | 0.7535572065523274 |  |  | horizon_steps=28; horizon_lyap=0.507136; mode=full |
+| E5 | warm_train_terminal | ContinuousESN | 0.21463110240826375 |  |  | horizon_steps=28; horizon_lyap=0.507136; mode=full |
+| E5 | cold | ContinuousESN | 0.6853049713043267 |  |  | horizon_steps=55; horizon_lyap=0.99616; mode=full |
+| E5 | warm_train_terminal | ContinuousESN | 0.14675113554137212 |  |  | horizon_steps=55; horizon_lyap=0.99616; mode=full |
+| E5 | cold | ContinuousESN | 1.3393780477247423 |  |  | horizon_steps=110; horizon_lyap=1.99232; mode=full |
+| E5 | warm_train_terminal | ContinuousESN | 0.08497214024585852 |  |  | horizon_steps=110; horizon_lyap=1.99232; mode=full |
+| E5 | cold | ContinuousESN | 1.506054216154576 |  |  | horizon_steps=166; horizon_lyap=3.006592; mode=full |
+| E5 | warm_train_terminal | ContinuousESN | 0.093388739272057 |  |  | horizon_steps=166; horizon_lyap=3.006592; mode=full |
+| E5 | cold | ContinuousESN | 1.6554797621233028 |  |  | horizon_steps=221; horizon_lyap=4.002752; mode=full |
+| E5 | warm_train_terminal | ContinuousESN | 0.13715717405603614 |  |  | horizon_steps=221; horizon_lyap=4.002752; mode=full |
+| E5 | cold | ContinuousESN | 1.5117492481722878 |  |  | horizon_steps=331; horizon_lyap=5.9950719999999995; mode=full |
+| E5 | warm_train_terminal | ContinuousESN | 0.6692777566540231 |  |  | horizon_steps=331; horizon_lyap=5.9950719999999995; mode=full |
+| E6 | continuous_st0 | ContinuousESN |  |  |  | has_carry=false; mode=full |
+| E6 | continuous_st_after_train | ContinuousESN |  |  |  | has_carry=false; mode=full |
+| E6 | continuous_st_after_collectstates | ContinuousESN |  |  |  | has_carry=false; mode=full |
+| E6 | discrete_st_after_train | ESN |  |  |  | has_carry=true; mode=full |
+| E7 | cold_fresh_st | ESN | 1.55491741812422 |  | 0.0 | mode=full |
+| E7 | warm_st_after_train | ESN | 1.1363799866710107 |  | 4.1114239999999995 | mode=full |
+| E7 | rewarm_collectstates_then_ar | ESN | 1.1363799866710107 |  | 4.1114239999999995 | mode=full |
+| E8 | cold | ContinuousESN | 1.4840906014114033 |  |  | washout=200; mode=full |
+| E8 | warm_full_train | ContinuousESN | 1.1839060794299072 |  |  | washout=200; mode=full |
+| E8 | warm_post_washout_tail | ContinuousESN | 1.1280934497932662 |  |  | washout=200; mode=full |
+| META | suite_wall_s | harness |  |  |  | wall_s=275.113651291; mode=full |
diff --git a/benchmarks/continuous_warmup/run.jl b/benchmarks/continuous_warmup/run.jl
new file mode 100644
index 000000000..bf42133af
--- /dev/null
+++ b/benchmarks/continuous_warmup/run.jl
@@ -0,0 +1,137 @@
+#!/usr/bin/env julia
+# Continuous AR predict warmup investigation entrypoint.
+#
+# Usage:
+#   julia --project=. run.jl --smoke
+#   julia --project=. run.jl
+#   julia --project=. run.jl --only=E1,E6
+
+using Pkg
+Pkg.instantiate()
+
+using Random
+using JSON
+using Dates
+using SciMLBase
+using DataInterpolations
+using OrdinaryDiffEqTsit5
+using OrdinaryDiffEq
+using ReservoirComputing
+
+const ROOT = @__DIR__
+include(joinpath(ROOT, "src", "metrics.jl"))
+include(joinpath(ROOT, "src", "data.jl"))
+include(joinpath(ROOT, "src", "models.jl"))
+include(joinpath(ROOT, "src", "predict_variants.jl"))
+include(joinpath(ROOT, "src", "experiments.jl"))
+
+function parse_args(args)
+    smoke = "--smoke" in args
+    # Alias kept for clarity in docs/CI; default without --smoke is full.
+    full = "--full" in args || !smoke
+    only = nothing
+    n_res = nothing
+    for a in args
+        if startswith(a, "--only=")
+            only = split(a[8:end], ','; keepempty = false)
+        elseif startswith(a, "--n-res=")
+            n_res = parse(Int, a[9:end])
+        end
+    end
+    return (; smoke, full, only, n_res)
+end
+
+function rows_to_markdown(rows)
+    io = IOBuffer()
+    println(io, "# Continuous warmup investigation results")
+    println(io)
+    println(io, "Generated: $(Dates.now())")
+    println(io)
+    println(io, "| exp | variant | model | nrmse | nrmse_global | vpt_lyap | notes |")
+    println(io, "|-----|---------|-------|-------|--------------|----------|-------|")
+    for r in rows
+        nrmse = get(r, "nrmse", "")
+        nrmse_g = get(r, "nrmse_global", "")
+        vpt = get(r, "vpt_lyap", "")
+        notes = String[]
+        for k in ("warmup_len", "horizon_steps", "horizon_lyap", "washout",
+            "max_abs_diff", "match_ok", "u0_norm", "has_carry",
+            "wall_train_s", "wall_ar_s", "wall_warmup_collect_s", "wall_s", "mode")
+            haskey(r, k) && push!(notes, "$k=$(r[k])")
+        end
+        println(
+            io,
+            "| $(r["experiment"]) | $(r["variant"]) | $(r["model"]) | ",
+            "$(nrmse) | $(nrmse_g) | $(vpt) | $(join(notes, "; ")) |",
+        )
+    end
+    return String(take!(io))
+end
+
+function main(args)
+    opts = parse_args(args)
+    results_dir = joinpath(ROOT, "results")
+    mkpath(results_dir)
+
+    data = opts.smoke ? make_lorenz_smoke() : make_lorenz()
+    n_res = something(opts.n_res, opts.smoke ? 80 : 300)
+    mode = opts.smoke ? "smoke" : "full"
+    cfg = (
+        smoke = opts.smoke,
+        mode = mode,
+        seed = 17,
+        n_res = n_res,
+        ridge = 1.0e-6,
+        data = data,
+    )
+
+    println("Continuous warmup investigation")
+    println("  mode=$(mode) n_res=$(cfg.n_res) train=$(data.train_len) ",
+        "predict=$(data.predict_len)  (t_λ span ≈ ",
+        round(data.predict_len * data.dt * data.λ_max; digits = 1), ")")
+    println("  only=$(opts.only === nothing ? "all" : join(opts.only, ","))")
+    println("  HPs: radius=0.9, input_scale=0.1, bias_scale=0.05, ridge=$(cfg.ridge)")
+
+    t_all0 = time_ns()
+    rows = run_experiments(cfg; only = opts.only)
+    wall_all = (time_ns() - t_all0) / 1.0e9
+    for r in rows
+        r["mode"] = mode
+        r["n_res_cfg"] = n_res
+    end
+    push!(
+        rows,
+        Dict{String, Any}(
+            "experiment" => "META",
+            "variant" => "suite_wall_s",
+            "model" => "harness",
+            "mode" => mode,
+            "wall_s" => wall_all,
+            "n_res" => n_res,
+            "train_len" => data.train_len,
+            "predict_len" => data.predict_len,
+        ),
+    )
+
+    tag = mode
+    json_path = joinpath(results_dir, "summary_$(tag).json")
+    open(json_path, "w") do io
+        JSON.print(io, rows, 2)
+    end
+    # Keep untagged names as the latest run for convenience.
+    cp(json_path, joinpath(results_dir, "summary.json"); force = true)
+
+    md_path = joinpath(results_dir, "summary_$(tag).md")
+    write(md_path, rows_to_markdown(rows))
+    cp(md_path, joinpath(results_dir, "summary.md"); force = true)
+
+    println("\nSuite wall time: $(round(wall_all; digits=1)) s")
+    println("Wrote $json_path")
+    println("Wrote $md_path")
+    println("Update results/FINDINGS.md with interpretation after review.")
+    return rows
+end
+
+if abspath(PROGRAM_FILE) == @__FILE__
+    main(ARGS)
+end
diff --git a/benchmarks/continuous_warmup/src/data.jl b/benchmarks/continuous_warmup/src/data.jl
new file mode 100644
index 000000000..18f16a9b5
--- /dev/null
+++ b/benchmarks/continuous_warmup/src/data.jl
@@ -0,0 +1,66 @@
+using Random
+using SciMLBase
+using OrdinaryDiffEqTsit5
+
+const LORENZ_λ_MAX = 0.9056  # typical max Lyapunov for classical Lorenz-63
+
+function lorenz!(du, u, p, t)
+    du[1] = p[1] * (u[2] - u[1])
+    du[2] = u[1] * (p[2] - u[3]) - u[2]
+    du[3] = u[1] * u[2] - p[3] * u[3]
+    return nothing
+end
+
+"""
+    make_lorenz(; dt, tspan, shift, train_len, predict_len, u0, p)
+
+Returns named tuple with `input_data`, `target_data`, `test_data`,
+`dt`, and full `data`.
+"""
+function make_lorenz(;
+        dt::Float64 = 0.02,
+        t_end::Float64 = 200.0,
+        shift::Int = 300,
+        train_len::Int = 5000,
+        # Match the #456 Lorenz probe (predict ≈ 1250 samples ≈ 22.6 t_λ at dt=0.02).
+        predict_len::Int = 1250,
+        u0 = [1.0, 0.0, 0.0],
+        p = [10.0, 28.0, 8 / 3],
+    )
+    data_prob = ODEProblem(lorenz!, u0, (0.0, t_end), p)
+    data = Array(solve(data_prob, Tsit5(); saveat = dt))
+
+    need = shift + train_len + predict_len
+    size(data, 2) ≥ need || throw(
+        ArgumentError(
+            "Lorenz series too short: need $need samples, got $(size(data, 2)). " *
+                "Increase t_end."
+        )
+    )
+
+    input_data = data[:, shift:(shift + train_len - 1)]
+    target_data = data[:, (shift + 1):(shift + train_len)]
+    test_data = data[:, (shift + train_len):(shift + train_len + predict_len - 1)]
+
+    return (
+        data = data,
+        input_data = input_data,
+        target_data = target_data,
+        test_data = test_data,
+        dt = dt,
+        shift = shift,
+        train_len = train_len,
+        predict_len = predict_len,
+        λ_max = LORENZ_λ_MAX,
+    )
+end
+
+"""Shorter Lorenz split for --smoke runs."""
+function make_lorenz_smoke()
+    return make_lorenz(;
+        t_end = 80.0,
+        shift = 100,
+        train_len = 800,
+        predict_len = 200,
+    )
+end
diff --git a/benchmarks/continuous_warmup/src/experiments.jl b/benchmarks/continuous_warmup/src/experiments.jl
new file mode 100644
index 000000000..345ed6567
--- /dev/null
+++ b/benchmarks/continuous_warmup/src/experiments.jl
@@ -0,0 +1,678 @@
+using Dates
+using JSON
+using Random
+using Statistics
+using LinearAlgebra
+
+"""Accumulate one experiment row."""
+function _row(; experiment, variant, model, kwargs...)
+    d = Dict{String, Any}(
+        "experiment" => string(experiment),
+        "variant" => string(variant),
+        "model" => string(model),
+        "timestamp" => string(Dates.now()),
+    )
+    for (k, v) in pairs(kwargs)
+        d[string(k)] = v
+    end
+    return d
+end
+
+function _score(pred, truth, data)
+    return (
+        nrmse = nrmse(pred, truth),
+        nrmse_global = nrmse_global(pred, truth),
+        vpt = valid_prediction_time(
+            pred, truth; dt = data.dt, λ_max = data.λ_max, threshold = 0.5
+        ),
+    )
+end
+
+function _horizons(predict_len)
+    # steps ≈ Lyapunov times * λ_max / dt  →  t_λ = steps * dt * λ_max
+    # invert: steps = t_λ / (dt * λ_max)
+    return nothing  # filled per-call with data
+end
+
+function lyap_horizons(data; t_λs = (0.5, 1.0, 2.0, 3.0, 4.0, 6.0))
+    steps = Int[]
+    for tλ in t_λs
+        h = max(1, round(Int, tλ / (data.dt * data.λ_max)))
+        h ≤ data.predict_len && push!(steps, h)
+    end
+    return unique(steps)
+end
+
+# ---------------------------------------------------------------------------
+# E1 — ContinuousESN cold vs train-terminal warm
+# ---------------------------------------------------------------------------
+
+function run_E1(cfg)
+    println("\n=== E1: ContinuousESN cold vs train-terminal warm ===")
+    rng = MersenneTwister(cfg.seed)
+    data = cfg.data
+    n_res = cfg.n_res
+
+    m_train = build_continuous_esn(n_res, data.train_len)
+    m_pred = build_continuous_esn(n_res, data.predict_len)
+    (tr, train_s) = timed() do
+        train_pair(m_train, m_pred, data; ridge = cfg.ridge, rng = rng)
+    end
+    ps_ar, st_ar = align_pred_params(
+        tr.ps_train, tr.st_train, tr.ps_pred, tr.st_pred
+    )
+
+    (cold, cold_s) = timed() do
+        predict_ar_cold(
+            tr.model_pred, data.predict_len, ps_ar, st_ar;
+            initialdata = data.test_data[:, 1],
+        )
+    end
+    cold_out = first(cold)
+    sc_cold = _score(cold_out, data.test_data, data)
+
+    (u0, warm_collect_s) = timed() do
+        raw_terminal_ode_state(
+            tr.model_train, data.input_data, tr.ps_train, tr.st_train
+        )
+    end
+    (warm, warm_ar_s) = timed() do
+        predict_ar_seeded(
+            tr.model_pred, data.predict_len, ps_ar, st_ar;
+            initialdata = data.test_data[:, 1],
+            initial_state = u0,
+        )
+    end
+    warm_out = first(warm)
+    sc_warm = _score(warm_out, data.test_data, data)
+
+    # Sanity: seeded with zeros should match package cold path
+    zero_u0 = zeros(eltype(u0), length(u0))
+    seeded_zero, _ = predict_ar_seeded(
+        tr.model_pred, data.predict_len, ps_ar, st_ar;
+        initialdata = data.test_data[:, 1],
+        initial_state = zero_u0,
+    )
+    match_cold = maximum(abs.(seeded_zero .- cold_out))
+
+    hs = lyap_horizons(data)
+    rows = [
+        _row(;
+            experiment = "E1",
+            variant = "cold_package_predict",
+            model = "ContinuousESN",
+            n_res = n_res,
+            nrmse = sc_cold.nrmse,
+            nrmse_global = sc_cold.nrmse_global,
+            vpt_lyap = sc_cold.vpt,
+            train_len = data.train_len,
+            predict_len = data.predict_len,
+            wall_train_s = train_s,
+            wall_ar_s = cold_s,
+        ),
+        _row(;
+            experiment = "E1",
+            variant = "warm_train_terminal_u0",
+            model = "ContinuousESN",
+            n_res = n_res,
+            nrmse = sc_warm.nrmse,
+            nrmse_global = sc_warm.nrmse_global,
+            vpt_lyap = sc_warm.vpt,
+            train_len = data.train_len,
+            predict_len = data.predict_len,
+            u0_norm = norm(u0),
+            wall_train_s = train_s,
+            wall_warmup_collect_s = warm_collect_s,
+            wall_ar_s = warm_ar_s,
+        ),
+        _row(;
+            experiment = "E1",
+            variant = "seeded_zero_matches_cold_maxabs",
+            model = "ContinuousESN",
+            n_res = n_res,
+            max_abs_diff = match_cold,
+            match_ok = match_cold < 1.0e-8,
+        ),
+    ]
+    for h in hs
+        tλ = h * data.dt * data.λ_max
+        push!(
+            rows,
+            _row(;
+                experiment = "E1",
+                variant = "cold_horizon",
+                model = "ContinuousESN",
+                horizon_steps = h,
+                horizon_lyap = tλ,
+                nrmse = nrmse(@view(cold_out[:, 1:h]), @view(data.test_data[:, 1:h])),
+            )
+        )
+        push!(
+            rows,
+            _row(;
+                experiment = "E1",
+                variant = "warm_horizon",
+                model = "ContinuousESN",
+                horizon_steps = h,
+                horizon_lyap = tλ,
+                nrmse = nrmse(@view(warm_out[:, 1:h]), @view(data.test_data[:, 1:h])),
+            )
+        )
+    end
+
+    @info "E1 cold nrmse=$(round(sc_cold.nrmse; digits=4)) warm=$(round(sc_warm.nrmse; digits=4)) " *
+        "Δ=$(round(sc_cold.nrmse - sc_warm.nrmse; digits=4)) zero_match=$(match_cold) " *
+        "train=$(round(train_s; digits=1))s ar_cold=$(round(cold_s; digits=1))s " *
+        "ar_warm=$(round(warm_ar_s; digits=1))s"
+    return rows
+end
+
+# ---------------------------------------------------------------------------
+# E2 — SciMLProblemReservoir eq.5 cold vs warm
+# ---------------------------------------------------------------------------
+
+function run_E2(cfg)
+    println("\n=== E2: SciMLProblemReservoir (eq.5) cold vs warm ===")
+    rng = MersenneTwister(cfg.seed)
+    data = cfg.data
+    n_res = cfg.n_res
+
+    rng2 = MersenneTwister(cfg.seed)
+    m_train = build_sciml_eq5(rng2, n_res, data.train_len)
+    # Same Wr/Win/b as train; only tspan / u0 sized for predict length.
+    res_tr = m_train.reservoir
+    tspan_p = (0.0, Float64(data.predict_len))
+    prob_p = remake(res_tr.prob; u0 = zeros(n_res), tspan = tspan_p)
+    res_p = SciMLProblemReservoir(
+        prob_p, TerminalStateSampling(), tspan_p, Tsit5();
+        reltol = RELTOL, abstol = ABSTOL,
+    )
+    m_pred = ReservoirComputer(res_p, m_train.states_modifiers, LinearReadout(n_res => 3))
+
+    tr = train_pair(m_train, m_pred, data; ridge = cfg.ridge, rng = MersenneTwister(cfg.seed))
+
+    cold, _ = predict_ar_cold(
+        tr.model_pred, data.predict_len, tr.ps_pred, tr.st_pred;
+        initialdata = data.test_data[:, 1],
+    )
+    sc_cold = _score(cold, data.test_data, data)
+
+    u0 = raw_terminal_ode_state(
+        tr.model_train, data.input_data, tr.ps_train, tr.st_train
+    )
+    warm, _ = predict_ar_seeded(
+        tr.model_pred, data.predict_len, tr.ps_pred, tr.st_pred;
+        initialdata = data.test_data[:, 1],
+        initial_state = u0,
+    )
+    sc_warm = _score(warm, data.test_data, data)
+
+    # remake public-style: set pred.prob.u0 then package predict
+    res_pred = tr.model_pred.reservoir
+    res_warm = SciMLProblemReservoir(
+        remake(res_pred.prob; u0 = u0),
+        res_pred.sampler,
+        res_pred.tspan,
+        res_pred.args,
+        res_pred.kwargs,
+    )
+    m_remake = ReservoirComputer(
+        res_warm, tr.model_pred.states_modifiers, tr.model_pred.readout
+    )
+    remake_out, _ = predict(
+        m_remake, data.predict_len, tr.ps_pred, tr.st_pred;
+        initialdata = data.test_data[:, 1],
+    )
+    sc_remake = _score(remake_out, data.test_data, data)
+    match_warm = maximum(abs.(remake_out .- warm))
+
+    rows = [
+        _row(;
+            experiment = "E2",
+            variant = "cold_package_predict",
+            model = "SciMLProblemReservoir",
+            n_res = n_res,
+            nrmse = sc_cold.nrmse,
+            nrmse_global = sc_cold.nrmse_global,
+            vpt_lyap = sc_cold.vpt,
+        ),
+        _row(;
+            experiment = "E2",
+            variant = "warm_train_terminal_u0",
+            model = "SciMLProblemReservoir",
+            n_res = n_res,
+            nrmse = sc_warm.nrmse,
+            nrmse_global = sc_warm.nrmse_global,
+            vpt_lyap = sc_warm.vpt,
+            u0_norm = norm(u0),
+        ),
+        _row(;
+            experiment = "E2",
+            variant = "remake_prob_u0_package_predict",
+            model = "SciMLProblemReservoir",
+            n_res = n_res,
+            nrmse = sc_remake.nrmse,
+            nrmse_global = sc_remake.nrmse_global,
+            vpt_lyap = sc_remake.vpt,
+            max_abs_diff_vs_seeded = match_warm,
+        ),
+    ]
+
+    @info "E2 cold=$(round(sc_cold.nrmse; digits=4)) warm=$(round(sc_warm.nrmse; digits=4)) " *
+        "remake=$(round(sc_remake.nrmse; digits=4))"
+    return rows
+end
+
+# ---------------------------------------------------------------------------
+# E3 — warmup length sweep
+# ---------------------------------------------------------------------------
+
+function run_E3(cfg)
+    println("\n=== E3: warmup length sweep (ContinuousESN) ===")
+    rng = MersenneTwister(cfg.seed)
+    data = cfg.data
+    n_res = cfg.n_res
+
+    m_train = build_continuous_esn(n_res, data.train_len)
+    m_pred = build_continuous_esn(n_res, data.predict_len)
+    tr = train_pair(m_train, m_pred, data; ridge = cfg.ridge, rng = rng)
+
+    ps_ar, st_ar = align_pred_params(
+        tr.ps_train, tr.st_train, tr.ps_pred, tr.st_pred
+    )
+
+    Ks = cfg.smoke ? [0, 10, 50, 100, 200] : [0, 10, 50, 100, 250, 500, 1000, 2000]
+    Ks = filter(k -> k == 0 || k ≤ data.train_len, Ks)
+
+    rows = Dict{String, Any}[]
+    for K in Ks
+        if K == 0
+            pred, _ = predict_ar_cold(
+                tr.model_pred, data.predict_len, ps_ar, st_ar;
+                initialdata = data.test_data[:, 1],
+            )
+            variant = "K=0_cold"
+        else
+            warm_data = data.input_data[:, (end - K + 1):end]
+            # collectstates on pred model needs tspan matching warm_data length
+            # ContinuousESN tspan is fixed to predict_len — length mismatch!
+            # Use train model for warmup collect (tspan = train_len), then AR on pred.
+            u0 = raw_terminal_ode_state(
+                tr.model_train, warm_data, tr.ps_train, tr.st_train
+            )
+            pred, _ = predict_ar_seeded(
+                tr.model_pred, data.predict_len, ps_ar, st_ar;
+                initialdata = data.test_data[:, 1],
+                initial_state = u0,
+            )
+            variant = "K=$(K)"
+        end
+        sc = _score(pred, data.test_data, data)
+        push!(
+            rows,
+            _row(;
+                experiment = "E3",
+                variant = variant,
+                model = "ContinuousESN",
+                warmup_len = K,
+                nrmse = sc.nrmse,
+                nrmse_global = sc.nrmse_global,
+                vpt_lyap = sc.vpt,
+            )
+        )
+        @info "E3 K=$K nrmse=$(round(sc.nrmse; digits=4)) vpt=$(round(sc.vpt; digits=3))"
+    end
+    return rows
+end
+
+# ---------------------------------------------------------------------------
+# E4 — seed variants
+# ---------------------------------------------------------------------------
+
+function run_E4(cfg)
+    println("\n=== E4: seed variants (ContinuousESN) ===")
+    rng = MersenneTwister(cfg.seed)
+    data = cfg.data
+    n_res = cfg.n_res
+
+    m_train = build_continuous_esn(n_res, data.train_len)
+    m_pred = build_continuous_esn(n_res, data.predict_len)
+    tr = train_pair(m_train, m_pred, data; ridge = cfg.ridge, rng = rng)
+    ps_ar, st_ar = align_pred_params(
+        tr.ps_train, tr.st_train, tr.ps_pred, tr.st_pred
+    )
+
+    u_train = raw_terminal_ode_state(
+        tr.model_train, data.input_data, tr.ps_train, tr.st_train
+    )
+    # test-prefix warm: first min(100, predict_len÷2) of test via teacher force
+    # but test is the target for AR — use teacher force on true test prefix
+    # as "oracle warm" upper bound
+    K_prefix = min(100, max(10, data.predict_len ÷ 5))
+    # For raw terminal on train model, use train-length tspan with a short
+    # series of test inputs placed as if they were a mini drive.
+    # Use train model with test prefix as data (grid scales to cell.tspan).
+    u_oracle = raw_terminal_ode_state(
+        tr.model_train,
+        data.test_data[:, 1:K_prefix],
+        tr.ps_train,
+        tr.st_train,
+    )
+
+    seeds = Dict(
+        "zeros" => zeros(n_res),
+        "randn" => randn(MersenneTwister(0), n_res),
+        "train_terminal" => u_train,
+        "oracle_test_prefix" => u_oracle,
+        "shuffled_train_terminal" => u_train[randperm(MersenneTwister(1), n_res)],
+    )
+
+    rows = Dict{String, Any}[]
+    for (name, u0) in seeds
+        pred, _ = predict_ar_seeded(
+            tr.model_pred, data.predict_len, ps_ar, st_ar;
+            initialdata = data.test_data[:, 1],
+            initial_state = Vector{Float64}(u0),
+        )
+        sc = _score(pred, data.test_data, data)
+        push!(
+            rows,
+            _row(;
+                experiment = "E4",
+                variant = name,
+                model = "ContinuousESN",
+                nrmse = sc.nrmse,
+                nrmse_global = sc.nrmse_global,
+                vpt_lyap = sc.vpt,
+                u0_norm = norm(u0),
+            )
+        )
+        @info "E4 $name nrmse=$(round(sc.nrmse; digits=4))"
+    end
+    return rows
+end
+
+# ---------------------------------------------------------------------------
+# E5 — horizon curves cold vs warm
+# ---------------------------------------------------------------------------
+
+function run_E5(cfg)
+    println("\n=== E5: horizon NRMSE curves ===")
+    rng = MersenneTwister(cfg.seed)
+    data = cfg.data
+    n_res = cfg.n_res
+
+    m_train = build_continuous_esn(n_res, data.train_len)
+    m_pred = build_continuous_esn(n_res, data.predict_len)
+    tr = train_pair(m_train, m_pred, data; ridge = cfg.ridge, rng = rng)
+    ps_ar, st_ar = align_pred_params(
+        tr.ps_train, tr.st_train, tr.ps_pred, tr.st_pred
+    )
+
+    cold, _ = predict_ar_cold(
+        tr.model_pred, data.predict_len, ps_ar, st_ar;
+        initialdata = data.test_data[:, 1],
+    )
+    u0 = raw_terminal_ode_state(
+        tr.model_train, data.input_data, tr.ps_train, tr.st_train
+    )
+    warm, _ = predict_ar_seeded(
+        tr.model_pred, data.predict_len, ps_ar, st_ar;
+        initialdata = data.test_data[:, 1],
+        initial_state = u0,
+    )
+
+    hs = lyap_horizons(data)
+    rows = Dict{String, Any}[]
+    for h in hs
+        tλ = h * data.dt * data.λ_max
+        for (variant, pred) in (("cold", cold), ("warm_train_terminal", warm))
+            sc = nrmse(@view(pred[:, 1:h]), @view(data.test_data[:, 1:h]))
+            push!(
+                rows,
+                _row(;
+                    experiment = "E5",
+                    variant = variant,
+                    model = "ContinuousESN",
+                    horizon_steps = h,
+                    horizon_lyap = tλ,
+                    nrmse = sc,
+                )
+            )
+        end
+        @info "E5 h=$h (tλ=$(round(tλ; digits=2))) cold=$(round(nrmse(cold[:,1:h], data.test_data[:,1:h]); digits=4)) " *
+            "warm=$(round(nrmse(warm[:,1:h], data.test_data[:,1:h]); digits=4))"
+    end
+    return rows
+end
+
+# ---------------------------------------------------------------------------
+# E6 — does st after train hold continuous carry?
+# ---------------------------------------------------------------------------
+
+function run_E6(cfg)
+    println("\n=== E6: inspect st after train! / collectstates ===")
+    rng = MersenneTwister(cfg.seed)
+    data = cfg.data
+    n_res = cfg.n_res
+
+    m = build_continuous_esn(n_res, data.train_len)
+    ps, st0 = setup(rng, m)
+    ps, st_train = train!(
+        m, data.input_data, data.target_data, ps, st0, StandardRidge(cfg.ridge)
+    )
+    _, st_collect = collectstates(m, data.input_data, ps, st_train)
+
+    d0 = inspect_st_reservoir(st0)
+    d1 = inspect_st_reservoir(st_train)
+    d2 = inspect_st_reservoir(st_collect)
+
+    m_disc = build_discrete_esn(n_res)
+    ps_d, st_d0 = setup(MersenneTwister(cfg.seed), m_disc)
+    ps_d, st_d = train!(
+        m_disc, data.input_data, data.target_data, ps_d, st_d0, StandardRidge(cfg.ridge)
+    )
+    d_disc = inspect_st_reservoir(st_d)
+
+    rows = [
+        _row(;
+            experiment = "E6",
+            variant = "continuous_st0",
+            model = "ContinuousESN",
+            st_type = d0.type,
+            keys = something(d0.keys, []),
+            has_carry = d0.has_carry,
+        ),
+        _row(;
+            experiment = "E6",
+            variant = "continuous_st_after_train",
+            model = "ContinuousESN",
+            st_type = d1.type,
+            keys = something(d1.keys, []),
+            has_carry = d1.has_carry,
+            summary = d1.summary,
+        ),
+        _row(;
+            experiment = "E6",
+            variant = "continuous_st_after_collectstates",
+            model = "ContinuousESN",
+            st_type = d2.type,
+            keys = something(d2.keys, []),
+            has_carry = d2.has_carry,
+            summary = d2.summary,
+        ),
+        _row(;
+            experiment = "E6",
+            variant = "discrete_st_after_train",
+            model = "ESN",
+            st_type = d_disc.type,
+            keys = something(d_disc.keys, []),
+            has_carry = d_disc.has_carry,
+            summary = d_disc.summary,
+        ),
+    ]
+
+    @info "E6 continuous has_carry after train=$(d1.has_carry); discrete has_carry=$(d_disc.has_carry)"
+    return rows
+end
+
+# ---------------------------------------------------------------------------
+# E7 — discrete ESN control
+# ---------------------------------------------------------------------------
+
+function run_E7(cfg)
+    println("\n=== E7: discrete ESN cold vs warm st ===")
+    rng = MersenneTwister(cfg.seed)
+    data = cfg.data
+    n_res = cfg.n_res
+
+    m = build_discrete_esn(n_res)
+    ps, st0 = setup(rng, m)
+    ps, st_train = train!(
+        m, data.input_data, data.target_data, ps, st0, StandardRidge(cfg.ridge)
+    )
+
+    # cold: fresh states, only readout trained
+    _, st_cold = setup(MersenneTwister(cfg.seed), m)
+    st_cold = merge(st_cold, (readout = st_train.readout,))
+    cold, _ = predict(
+        m, data.predict_len, ps, st_cold; initialdata = data.test_data[:, 1]
+    )
+
+    # warm: continue from post-train states
+    warm, _ = predict(
+        m, data.predict_len, ps, st_train; initialdata = data.test_data[:, 1]
+    )
+
+    # re-drive last train window then AR (like warmup_data)
+    _, st_rewarm = collectstates(m, data.input_data, ps, st0)
+    st_rewarm = merge(st_rewarm, (readout = st_train.readout,))
+    rewarm, _ = predict(
+        m, data.predict_len, ps, st_rewarm; initialdata = data.test_data[:, 1]
+    )
+
+    sc_c = _score(cold, data.test_data, data)
+    sc_w = _score(warm, data.test_data, data)
+    sc_r = _score(rewarm, data.test_data, data)
+
+    rows = [
+        _row(;
+            experiment = "E7",
+            variant = "cold_fresh_st",
+            model = "ESN",
+            nrmse = sc_c.nrmse,
+            vpt_lyap = sc_c.vpt,
+        ),
+        _row(;
+            experiment = "E7",
+            variant = "warm_st_after_train",
+            model = "ESN",
+            nrmse = sc_w.nrmse,
+            vpt_lyap = sc_w.vpt,
+        ),
+        _row(;
+            experiment = "E7",
+            variant = "rewarm_collectstates_then_ar",
+            model = "ESN",
+            nrmse = sc_r.nrmse,
+            vpt_lyap = sc_r.vpt,
+        ),
+    ]
+    @info "E7 cold=$(round(sc_c.nrmse; digits=4)) warm_train_st=$(round(sc_w.nrmse; digits=4)) " *
+        "rewarm=$(round(sc_r.nrmse; digits=4))"
+    return rows
+end
+
+# ---------------------------------------------------------------------------
+# E8 — washout interaction
+# ---------------------------------------------------------------------------
+
+function run_E8(cfg)
+    println("\n=== E8: washout + warm (ContinuousESN) ===")
+    rng = MersenneTwister(cfg.seed)
+    data = cfg.data
+    n_res = cfg.n_res
+    washout = cfg.smoke ? 50 : 200
+
+    m_train = build_continuous_esn(n_res, data.train_len)
+    m_pred = build_continuous_esn(n_res, data.predict_len)
+    tr = train_pair(
+        m_train, m_pred, data; ridge = cfg.ridge, rng = rng, washout = washout
+    )
+    ps_ar, st_ar = align_pred_params(
+        tr.ps_train, tr.st_train, tr.ps_pred, tr.st_pred
+    )
+
+    cold, _ = predict_ar_cold(
+        tr.model_pred, data.predict_len, ps_ar, st_ar;
+        initialdata = data.test_data[:, 1],
+    )
+    u_full = raw_terminal_ode_state(
+        tr.model_train, data.input_data, tr.ps_train, tr.st_train
+    )
+    # post-washout tail only
+    tail = data.input_data[:, (washout + 1):end]
+    u_tail = if size(tail, 2) ≥ 2
+        raw_terminal_ode_state(tr.model_train, tail, tr.ps_train, tr.st_train)
+    else
+        u_full
+    end
+
+    warm_full, _ = predict_ar_seeded(
+        tr.model_pred, data.predict_len, ps_ar, st_ar;
+        initialdata = data.test_data[:, 1], initial_state = u_full,
+    )
+    warm_tail, _ = predict_ar_seeded(
+        tr.model_pred, data.predict_len, ps_ar, st_ar;
+        initialdata = data.test_data[:, 1], initial_state = u_tail,
+    )
+
+    rows = [
+        _row(;
+            experiment = "E8",
+            variant = "cold",
+            model = "ContinuousESN",
+            washout = washout,
+            nrmse = nrmse(cold, data.test_data),
+        ),
+        _row(;
+            experiment = "E8",
+            variant = "warm_full_train",
+            model = "ContinuousESN",
+            washout = washout,
+            nrmse = nrmse(warm_full, data.test_data),
+        ),
+        _row(;
+            experiment = "E8",
+            variant = "warm_post_washout_tail",
+            model = "ContinuousESN",
+            washout = washout,
+            nrmse = nrmse(warm_tail, data.test_data),
+        ),
+    ]
+    @info "E8 washout=$washout cold=$(round(rows[1]["nrmse"]; digits=4)) " *
+        "full=$(round(rows[2]["nrmse"]; digits=4)) tail=$(round(rows[3]["nrmse"]; digits=4))"
+    return rows
+end
+
+const EXPERIMENT_FUNS = Dict(
+    "E1" => run_E1,
+    "E2" => run_E2,
+    "E3" => run_E3,
+    "E4" => run_E4,
+    "E5" => run_E5,
+    "E6" => run_E6,
+    "E7" => run_E7,
+    "E8" => run_E8,
+)
+
+function run_experiments(cfg; only = nothing)
+    ids = only === nothing ? collect(keys(EXPERIMENT_FUNS)) : only
+    sort!(ids)
+    all_rows = Dict{String, Any}[]
+    for id in ids
+        haskey(EXPERIMENT_FUNS, id) || throw(ArgumentError("unknown experiment $id"))
+        append!(all_rows, EXPERIMENT_FUNS[id](cfg))
+    end
+    return all_rows
+end
diff --git a/benchmarks/continuous_warmup/src/metrics.jl b/benchmarks/continuous_warmup/src/metrics.jl
new file mode 100644
index 000000000..ddd6dc5e8
--- /dev/null
+++ b/benchmarks/continuous_warmup/src/metrics.jl
@@ -0,0 +1,78 @@
+using Statistics
+
+"""
+    nrmse(pred, truth)
+
+Per-channel NRMSE (Lukoševičius-style scale by channel std), then average
+over channels. `pred` and `truth` are `(n_channels, T)`.
+"""
+function nrmse(pred::AbstractMatrix, truth::AbstractMatrix)
+    size(pred) == size(truth) ||
+        throw(DimensionMismatch("pred $(size(pred)) vs truth $(size(truth))"))
+    n_ch = size(pred, 1)
+    acc = 0.0
+    for c in 1:n_ch
+        y = @view truth[c, :]
+        ŷ = @view pred[c, :]
+        σ = std(y)
+        σ = σ > 0 ? σ : one(σ)
+        acc += sqrt(mean((ŷ .- y) .^ 2)) / σ
+    end
+    return acc / n_ch
+end
+
+"""
+    nrmse_global(pred, truth)
+
+Single scalar NRMSE using global std of `truth` (matches some #456 notes).
+"""
+function nrmse_global(pred::AbstractMatrix, truth::AbstractMatrix)
+    return sqrt(mean((pred .- truth) .^ 2)) / std(truth)
+end
+
+"""
+    valid_prediction_time(pred, truth; dt, λ_max, threshold=0.5)
+
+First time (in Lyapunov units) where per-step relative error exceeds
+`threshold`. Relative error at step k is
+`‖pred[:,k] - truth[:,k]‖ / √(mean(truth.^2) + ε)`.
+Returns `Inf` if never exceeded.
+"""
+function valid_prediction_time(
+        pred::AbstractMatrix,
+        truth::AbstractMatrix;
+        dt::Real,
+        λ_max::Real,
+        threshold::Real = 0.5
+    )
+    scale = sqrt(mean(truth .^ 2)) + eps(Float64)
+    T = size(pred, 2)
+    for k in 1:T
+        err = sqrt(sum(abs2, @view(pred[:, k]) .- @view(truth[:, k]))) / scale
+        if err > threshold
+            return (k - 1) * dt * λ_max
+        end
+    end
+    return Inf
+end
+
+"""
+    horizon_nrmse(pred, truth; horizons)
+
+NRMSE on prefixes of length `h` for each `h` in `horizons`.
+"""
+function horizon_nrmse(pred::AbstractMatrix, truth::AbstractMatrix; horizons)
+    return Dict(string(h) => nrmse(@view(pred[:, 1:h]), @view(truth[:, 1:h])) for h in horizons)
+end
+
+"""
+    timed(f) -> (result, wall_s)
+
+Wall-clock seconds for `f()` via `time_ns` (includes compile on first call —
+callers should warm once if they need steady-state timings).
+"""
+function timed(f)
+    t0 = time_ns()
+    result = f()
+    return result, (time_ns() - t0) / 1.0e9
+end
diff --git a/benchmarks/continuous_warmup/src/models.jl b/benchmarks/continuous_warmup/src/models.jl
new file mode 100644
index 000000000..18d126417
--- /dev/null
+++ b/benchmarks/continuous_warmup/src/models.jl
@@ -0,0 +1,190 @@
+using Random
+using LinearAlgebra
+using ReservoirComputing
+using SciMLBase
+using OrdinaryDiffEqTsit5
+
+const RELTOL = 1.0e-6
+const ABSTOL = 1.0e-8
+
+init_input_f64(rng, dims...) = scaled_rand(rng, Float64, dims...)
+init_reservoir_f64(rng, dims...) = rand_sparse(rng, Float64, dims...)
+init_bias_f64(rng, dims...) = zeros(Float64, dims...)
+init_state_f64(rng, dims...) = zeros(Float64, dims...)
+
+"""
+#456-style continuous inits: spectral radius 0.9, modest input drive (~0.1),
+small bias (~0.05). Package `scaled_rand` defaults are too aggressive for
+continuous eq. (5) forecasting and wash out the cold/warm signal.
+"""
+function continuous_inits(;
+        use_bias::Bool = true,
+        radius::Float64 = 0.9,
+        input_scale::Float64 = 0.1,
+        bias_scale::Float64 = 0.05,
+    )
+    init_res = (rng, dims...) -> rand_sparse(rng, Float64, dims...; radius = radius)
+    init_in = (rng, dims...) -> input_scale .* randn(rng, Float64, dims...)
+    init_b = if use_bias
+        (rng, dims...) -> bias_scale .* randn(rng, Float64, dims...)
+    else
+        init_bias_f64
+    end
+    return (
+        use_bias = use_bias,
+        init_input = init_in,
+        init_reservoir = init_res,
+        init_bias = init_b,
+        init_state = init_state_f64,
+        reltol = RELTOL,
+        abstol = ABSTOL,
+    )
+end
+
+function f64_inits(; use_bias = true)
+    return (
+        use_bias = use_bias,
+        init_input = init_input_f64,
+        init_reservoir = init_reservoir_f64,
+        init_bias = init_bias_f64,
+        init_state = init_state_f64,
+        reltol = RELTOL,
+        abstol = ABSTOL,
+    )
+end
+
+"""
+    build_continuous_esn(n_res, n_steps; radius=0.9, kwargs...)
+
+`tspan = (0, n_steps)` so one sample ≈ one unit of reservoir time
+(tutorial / #456 convention). Defaults match the #456 Lorenz probe.
+"""
+function build_continuous_esn(
+        n_res::Integer,
+        n_steps::Integer;
+        state_modifiers = (NLAT2(),),
+        use_bias::Bool = true,
+        solver = Tsit5(),
+        radius::Float64 = 0.9,
+        input_scale::Float64 = 0.1,
+        bias_scale::Float64 = 0.05,
+        extra...,
+    )
+    return ContinuousESN(
+        3, n_res, 3, (0.0, Float64(n_steps)), solver;
+        continuous_inits(; use_bias, radius, input_scale, bias_scale)...,
+        state_modifiers = state_modifiers,
+        extra...,
+    )
+end
+
+"""
+Copy trained reservoir + readout parameters into a predict-length model
+setup so cold/warm AR use the **same** `W_in` / `W_r` / bias / `W_out`.
+"""
+function align_pred_params(ps_train, st_train, ps_pred, st_pred)
+    ps = merge(ps_pred, (reservoir = ps_train.reservoir, readout = ps_train.readout))
+    st = merge(st_pred, (readout = st_train.readout,))
+    return ps, st
+end
+
+"""
+    eq5_rhs!(dx, x, p, t)
+
+Lukoševičius §3.2.6 eq. (5) for hand-rolled `SciMLProblemReservoir`.
+"""
+function eq5_rhs!(dx, x, p, t)
+    input_t = p.input(t)
+    # Match hand-rolled #456 / test style: Wr, Win, b in `prob.p`.
+    dx .= .-x .+ tanh.(p.Wr * x .+ p.Win * input_t .+ p.b)
+    return nothing
+end
+
+"""
+    build_sciml_eq5(rng, n_res, n_steps; radius=0.9, …)
+
+Sparse random Wr rescaled to spectral radius `radius` (continuous ESP-ish).
+"""
+function build_sciml_eq5(
+        rng::AbstractRNG,
+        n_res::Integer,
+        n_steps::Integer;
+        radius::Float64 = 0.9,
+        input_scale::Float64 = 0.1,
+        use_bias::Bool = true,
+        sparsity::Float64 = 6 / n_res,
+        state_modifiers = (NLAT2(),),
+        solver = Tsit5(),
+    )
+    Wr_raw = randn(rng, n_res, n_res)
+    mask = rand(rng, n_res, n_res) .< sparsity
+    Wr_sparse = Wr_raw .* mask
+    ρ = maximum(abs.(eigvals(Matrix(Wr_sparse))))
+    ρ = ρ > 0 ? ρ : 1.0
+    Wr = (radius / ρ) .* Wr_sparse
+    Win = input_scale .* randn(rng, n_res, 3)
+    bias = use_bias ? (0.05 .* randn(rng, n_res)) : zeros(n_res)
+    p0 = (Wr = Wr, Win = Win, b = bias)
+
+    tspan = (0.0, Float64(n_steps))
+    u0 = zeros(n_res)
+    prob = ODEProblem(eq5_rhs!, u0, tspan, p0)
+    res = SciMLProblemReservoir(
+        prob, TerminalStateSampling(), tspan, solver;
+        reltol = RELTOL, abstol = ABSTOL,
+    )
+    return ReservoirComputer(res, state_modifiers, LinearReadout(n_res => 3))
+end
+
+function build_discrete_esn(
+        n_res::Integer;
+        state_modifiers = (NLAT2(),),
+        use_bias::Bool = true,
+        radius::Float64 = 0.9,
+    )
+    return ESN(
+        3, n_res, 3;
+        use_bias = use_bias,
+        init_input = init_input_f64,
+        init_reservoir = (rng, dims...) -> rand_sparse(
+            rng, Float64, dims...; radius = radius
+        ),
+        init_bias = init_bias_f64,
+        init_state = init_state_f64,
+        state_modifiers = state_modifiers,
+    )
+end
+
+"""
+    train_pair(model_train, model_pred, data; washout, ridge, rng)
+
+Train on `model_train`, copy readout into a fresh `model_pred` parameter
+set (same pattern as continuous tutorial).
+"""
+function train_pair(
+        model_train,
+        model_pred,
+        data;
+        washout::Int = 0,
+        ridge::Float64 = 1.0e-6,
+        rng = MersenneTwister(17),
+    )
+    ps, st = setup(rng, model_train)
+    ps, st = train!(
+        model_train, data.input_data, data.target_data, ps, st,
+        StandardRidge(ridge); washout = washout
+    )
+
+    ps_pred, st_pred = setup(rng, model_pred)
+    ps_pred = merge(ps_pred, (readout = ps.readout,))
+    st_pred = merge(st_pred, (readout = st.readout,))
+
+    return (
+        ps_train = ps,
+        st_train = st,
+        ps_pred = ps_pred,
+        st_pred = st_pred,
+        model_train = model_train,
+        model_pred = model_pred,
+    )
+end
diff --git a/benchmarks/continuous_warmup/src/predict_variants.jl b/benchmarks/continuous_warmup/src/predict_variants.jl
new file mode 100644
index 000000000..171456e65
--- /dev/null
+++ b/benchmarks/continuous_warmup/src/predict_variants.jl
@@ -0,0 +1,278 @@
+# Experimental AR predict variants — mirror extension logic without
+# changing package API. Used only for investigation.
+
+using ReservoirComputing  # reexports LuxCore.apply
+using SciMLBase
+using OrdinaryDiffEq
+
+function _require_ext()
+    ext = Base.get_extension(ReservoirComputing, :RCODEReservoirExt)
+    ext === nothing && error(
+        "RCODEReservoirExt not loaded. Import SciMLBase, DataInterpolations, " *
+            "and an ODE solver package first."
+    )
+    return ext
+end
+
+"""
+    terminal_state_from_collect(rc, data, ps, st)
+
+Teacher-forced `collectstates` then last state column (raw continuous
+state after modifiers are applied — same matrix the readout sees).
+
+Note: for continuous models this is the **sampled** reservoir feature
+after modifiers, not necessarily the raw ODE `u`. For AR seeding we need
+the raw ODE state. Prefer [`raw_terminal_ode_state`](@ref) when seeding `u0`.
+"""
+function terminal_state_from_collect(rc, data::AbstractMatrix, ps, st)
+    states, new_st = collectstates(rc, data, ps, st)
+    return states[:, end], new_st
+end
+
+"""
+    raw_terminal_ode_state(rc, data, ps, st)
+
+Re-run the continuous collect path and return the last **unmodified**
+ODE state (before state_modifiers). This is the correct seed for `u0`.
+"""
+function raw_terminal_ode_state(rc, data::AbstractMatrix, ps, st)
+    ext = _require_ext()
+    res = rc.reservoir
+
+    if res isa ContinuousESNCell
+        return _raw_terminal_continuous_esn(ext, res, rc, data, ps, st)
+    elseif res isa AbstractSciMLProblemReservoir
+        return _raw_terminal_sciml(ext, res, rc, data, ps, st)
+    else
+        # Discrete: collectstates last column is the recurrent carry
+        # after modifiers — for ESN without modifiers this is fine.
+        states, _ = collectstates(rc, data, ps, st)
+        return copy(states[:, end])
+    end
+end
+
+function _raw_terminal_continuous_esn(ext, cell::ContinuousESNCell, rc, data, ps, st)
+    n_samples = size(data, 2)
+    # Unit window width (Δt = 1), independent of the model’s train/predict
+    # `tspan`. Using `cell.tspan` with a short warmup of length K would stretch
+    # those K samples across the full train interval and change the ODE pacing.
+    t0, t1 = 0.0, Float64(n_samples)
+    Δt = (t1 - t0) / n_samples  # == 1
+    input_ts = collect(range(t0, t1 - Δt; length = n_samples))
+    sample_ts = collect(range(t0 + Δt, t1; length = n_samples))
+
+    input_interp = ext._make_input_fn(data, input_ts)
+    solve_p = ext._build_solve_params(nothing, ps.reservoir, input_interp)
+    u0 = zeros(eltype(ps.reservoir.input_matrix), cell.out_dims)
+    prob = ODEProblem(cell.equations, u0, (t0, t1), solve_p)
+    sol = solve(
+        prob, cell.args...;
+        saveat = sample_ts,
+        save_everystep = false,
+        dense = false,
+        cell.kwargs...,
+    )
+    return copy(sol.u[end])
+end
+
+function _raw_terminal_sciml(ext, res, rc, data, ps, st)
+    n_samples = size(data, 2)
+    t0, t1 = 0.0, Float64(n_samples)
+    Δt = (t1 - t0) / n_samples
+    input_ts = collect(range(t0, t1 - Δt; length = n_samples))
+    sample_ts = collect(range(t0 + Δt, t1; length = n_samples))
+
+    input_interp = ext._make_input_fn(data, input_ts)
+    solve_p = ext._build_solve_params(res.prob.p, ps.reservoir, input_interp)
+    remade = remake(res.prob; tspan = (t0, t1), p = solve_p)
+    sol = solve(
+        remade, res.args...;
+        saveat = sample_ts,
+        save_everystep = false,
+        dense = false,
+        res.kwargs...,
+    )
+    return copy(sol.u[end])
+end
+
+"""
+    predict_ar_seeded(rc, steps, ps, st; initialdata, initial_state)
+
+Autoregressive rollout with an explicit reservoir seed `initial_state`.
+Mirrors `RCODEReservoirExt` continuous AR `predict` for both
+`ContinuousESNCell` and generic `SciMLProblemReservoir`.
+"""
+function predict_ar_seeded(
+        rc,
+        steps::Integer,
+        ps,
+        st;
+        initialdata::AbstractVector,
+        initial_state::AbstractVector,
+    )
+    steps ≥ 1 || throw(ArgumentError("steps must be ≥ 1, got $steps"))
+    ext = _require_ext()
+    res = rc.reservoir
+
+    if res isa ContinuousESNCell
+        return _ar_continuous_esn(ext, res, rc, steps, ps, st; initialdata, initial_state)
+    elseif res isa AbstractSciMLProblemReservoir
+        return _ar_sciml(ext, res, rc, steps, ps, st; initialdata, initial_state)
+    else
+        throw(ArgumentError("predict_ar_seeded only for continuous reservoirs, got $(typeof(res))"))
+    end
+end
+
+function _ar_continuous_esn(
+        ext, cell, rc, steps, ps, st;
+        initialdata, initial_state,
+    )
+    t0, t1 = cell.tspan
+    ts = collect(range(t0, t1; length = steps + 1))
+    window_starts = @view ts[1:(end - 1)]
+    window_ends = @view ts[2:end]
+
+    current_state = copy(initial_state)
+    current_input = initialdata
+    st_mods = st.states_modifiers
+    st_ro = st.readout
+
+    local outputs
+    for (step_idx, (t_lo, t_hi)) in enumerate(zip(window_starts, window_ends))
+        input_fn = ext._make_const_input_fn(current_input, t_lo, t_hi)
+        solve_p = ext._build_solve_params(nothing, ps.reservoir, input_fn)
+        sub_prob = ODEProblem(cell.equations, current_state, (t_lo, t_hi), solve_p)
+        sol = solve(
+            sub_prob, cell.args...;
+            saveat = [t_hi],
+            save_everystep = false,
+            dense = false,
+            cell.kwargs...,
+        )
+        current_state = sol.u[end]
+
+        if !isempty(rc.states_modifiers)
+            state_after_mods, st_mods = ReservoirComputing._apply_seq(
+                rc.states_modifiers, current_state, ps.states_modifiers, st_mods
+            )
+        else
+            state_after_mods = current_state
+        end
+
+        current_output, st_ro = apply(rc.readout, state_after_mods, ps.readout, st_ro)
+        if step_idx == 1
+            outputs = similar(current_output, length(current_output), steps)
+        end
+        outputs[:, step_idx] .= current_output
+        current_input = current_output
+    end
+
+    newst = (
+        reservoir = st.reservoir,
+        states_modifiers = st_mods,
+        readout = st_ro,
+    )
+    return outputs, newst
+end
+
+function _ar_sciml(
+        ext, res, rc, steps, ps, st;
+        initialdata, initial_state,
+    )
+    t0, t1 = res.tspan
+    ts = collect(range(t0, t1; length = steps + 1))
+    window_starts = @view ts[1:(end - 1)]
+    window_ends = @view ts[2:end]
+
+    current_state = copy(initial_state)
+    current_input = initialdata
+    st_mods = st.states_modifiers
+    st_ro = st.readout
+
+    local outputs
+    for (step_idx, (t_lo, t_hi)) in enumerate(zip(window_starts, window_ends))
+        input_fn = ext._make_const_input_fn(current_input, t_lo, t_hi)
+        solve_p = ext._build_solve_params(res.prob.p, ps.reservoir, input_fn)
+        sub_prob = remake(
+            res.prob;
+            tspan = (t_lo, t_hi),
+            p = solve_p,
+            u0 = current_state,
+        )
+        sol = solve(
+            sub_prob, res.args...;
+            saveat = [t_hi],
+            save_everystep = false,
+            dense = false,
+            res.kwargs...,
+        )
+        current_state = sol.u[end]
+
+        if !isempty(rc.states_modifiers)
+            state_after_mods, st_mods = ReservoirComputing._apply_seq(
+                rc.states_modifiers, current_state, ps.states_modifiers, st_mods
+            )
+        else
+            state_after_mods = current_state
+        end
+
+        current_output, st_ro = apply(rc.readout, state_after_mods, ps.readout, st_ro)
+        if step_idx == 1
+            outputs = similar(current_output, length(current_output), steps)
+        end
+        outputs[:, step_idx] .= current_output
+        current_input = current_output
+    end
+
+    newst = (
+        reservoir = st.reservoir,
+        states_modifiers = st_mods,
+        readout = st_ro,
+    )
+    return outputs, newst
+end
+
+"""
+    predict_ar_warmup(rc, steps, ps, st; initialdata, warmup_data)
+
+Teacher-force on `warmup_data`, seed AR from raw terminal ODE state,
+first AR input = `initialdata` (default: last column of warmup or
+first test input — caller passes explicitly).
+"""
+function predict_ar_warmup(
+        rc,
+        steps::Integer,
+        ps,
+        st;
+        initialdata::AbstractVector,
+        warmup_data::AbstractMatrix,
+    )
+    u0 = raw_terminal_ode_state(rc, warmup_data, ps, st)
+    return predict_ar_seeded(
+        rc, steps, ps, st;
+        initialdata = initialdata,
+        initial_state = u0,
+    )
+end
+
+"""Package AR predict — cold start as implemented on master."""
+function predict_ar_cold(rc, steps, ps, st; initialdata)
+    return predict(rc, steps, ps, st; initialdata = initialdata)
+end
+
+"""
+Inspect whether `st` after continuous train/collect carries a hidden
+state usable as `u0`. Returns a NamedTuple of diagnostics.
+"""
+function inspect_st_reservoir(st)
+    r = st.reservoir
+    return (
+        type = string(typeof(r)),
+        is_namedtuple = r isa NamedTuple,
+        keys = r isa NamedTuple ? collect(keys(r)) : nothing,
+        has_carry = r isa NamedTuple && haskey(r, :carry),
+        has_cell = r isa NamedTuple && haskey(r, :cell),
+        # ContinuousESNCell initialstates is typically (rng = …) only
+        summary = sprint(show, r),
+    )
+end