Skip to content

docs: fix continuous tutorial forecasts #506

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
MartinuzziFrancesco merged 1 commit into from
Aug 19, 2026
+57 −43
Merged

docs: fix continuous tutorial forecasts #506

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
51 changes: 29 additions & 22 deletions docs/src/tutorials/continuous_esn.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,30 +46,29 @@ test = data[:, (shift + train_len):(shift + train_len + predict_len - 1)]
## Constructing the `ContinuousESN`

```@example continuous-esn-lorenz
N_res = 100
N_res = 300
res_radius = 0.9
res_sparsity = 6 / N_res

# Float64 initialisers so the reservoir, the solve, and the input all
# share a numeric type. Without these the cell would default to
# Float32 via `scaled_rand` / `rand_sparse` / `zeros32`.
init_input_f64(rng, d...) = scaled_rand(rng, Float64, d...)
init_reservoir_f64(rng, d...) = rand_sparse(rng, Float64, d...)
init_bias_f64(rng, d...) = zeros(Float64, d...)
init_input_f64(rng, d...) = scaled_rand(rng, Float64, d...)
init_reservoir_f64(rng, d...) = rand_sparse(
rng, Float64, d...; radius = res_radius, sparsity = res_sparsity
)

esn_train = ContinuousESN(
3, N_res, 3, (0.0, Float64(train_len)), Tsit5();
use_bias = true,
init_input = init_input_f64,
init_reservoir = init_reservoir_f64,
init_bias = init_bias_f64,
state_modifiers = (NLAT2(),),
reltol = 1.0e-6, abstol = 1.0e-8
)
esn_pred = ContinuousESN(
3, N_res, 3, (0.0, Float64(predict_len)), Tsit5();
use_bias = true,
init_input = init_input_f64,
init_reservoir = init_reservoir_f64,
init_bias = init_bias_f64,
state_modifiers = (NLAT2(),),
reltol = 1.0e-6, abstol = 1.0e-8
)
Expand All @@ -80,28 +79,36 @@ ps, st = setup(rng, esn_train)
## Training

```@example continuous-esn-lorenz
ps, st = train(esn_train, input_data, target_data, ps, st)
ps, st = train(esn_train, input_data, target_data, ps, st;
objective = RidgeRegression(1.0e-6))
```

## Autoregressive rollout

```@example continuous-esn-lorenz
ps_pred, st_pred = setup(rng, esn_pred)
ps_pred = merge(ps_pred, (readout = ps.readout,))
st_pred = merge(st_pred, (readout = st.readout,))

output, _ = predict(
esn_pred, predict_len, ps_pred, st_pred; initialdata = test[:, 1]
esn_pred, predict_len, ps, st; initialdata = test[:, 1]
)
```

plot(
transpose(output)[:, 1], transpose(output)[:, 2],
transpose(output)[:, 3]; label = "predicted"
)
plot!(
transpose(test)[:, 1], transpose(test)[:, 2],
transpose(test)[:, 3]; label = "actual"
)
```@example continuous-esn-lorenz
using Plots.PlotMeasures

dt = 0.02
lorenz_maxlyap = 0.9056
lyap_time = (0:(predict_len - 1)) .* dt .* (1 / lorenz_maxlyap)

p1 = plot(lyap_time, [test[1, :] output[1, :]]; label = ["actual" "predicted"],

@chatgpt-codex-connector chatgpt-codex-connector Bot Aug 18, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P2 Badge Align actual series with one-step forecasts

When warm-starting from the st returned by train, test[:, 1] is the seed value already used as initialdata, and predict stores the output only after advancing one autoregressive step. This means output[:, 1] forecasts the sample after test[:, 1], but the plot compares it against test[:, 1] itself, shifting the actual curve by one step; the same pattern appears in the SciML reservoir Lorenz and Mackey-Glass plots. Use an actual slice that starts at the next sample, or compare against test[:, 2:end] with matching prediction length.

Useful? React with 👍 / 👎.

@Saswatsusmoy Saswatsusmoy Aug 18, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same overlay as lorenz_basic / getting started / the readme. predict always writes the first forecast after the seed, so those plots are also one sample off at dt=0.02. not visible on the lyapunov-time axes and not what made the old continuous plots collapse. leaving it so this page stays in lockstep with the discrete example.

ylabel = "x(t)", linewidth = 2.5, xticks = false, yticks = -15:15:15);
p2 = plot(lyap_time, [test[2, :] output[2, :]]; label = ["actual" "predicted"],
ylabel = "y(t)", linewidth = 2.5, xticks = false, yticks = -20:20:20);
p3 = plot(lyap_time, [test[3, :] output[3, :]]; label = ["actual" "predicted"],
ylabel = "z(t)", linewidth = 2.5, xlabel = "max(λ)*t", yticks = 10:15:40);

plot(p1, p2, p3; plot_title = "Lorenz System Coordinates",
layout = (3, 1), xtickfontsize = 12, ytickfontsize = 12, xguidefontsize = 15,
yguidefontsize = 15,
legendfontsize = 12, titlefontsize = 20)
```

The two trajectories agree on the early portion of the rollout before
Expand Down
49 changes: 28 additions & 21 deletions docs/src/tutorials/sciml_reservoir.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -133,9 +133,9 @@ predict(rc, steps, ps, st; initialdata) # autoregressive rollout
the next sub-interval, and stitches the per-window readouts into
the returned output matrix.

In both cases the reservoir's initial state is `prob.u0`. To continue
from a previously computed trajectory, `remake(prob; u0 = …)` before
constructing the reservoir.
After `train`, pass the same `st` into `predict` so the rollout
continues from the trained reservoir state. A fresh `st` starts from
`prob.u0`.

## Eye test: Lorenz chaos forecasting with a continuous ESN

Expand All @@ -159,6 +159,7 @@ using SciMLBase
using DataInterpolations
using OrdinaryDiffEqTsit5
using Plots
using Plots.PlotMeasures
using Random

Random.seed!(42)
Expand All @@ -179,10 +180,10 @@ target_data = data[:, (shift + 1):(shift + train_len)]
test = data[:, (shift + train_len):(shift + train_len + predict_len - 1)]

# 2. Continuous ESN reservoir parameters
N_res = 100
Wr = 0.3 .* randn(rng, N_res, N_res) ./ sqrt(N_res)
Win = 0.5 .* randn(rng, N_res, 3)
bias = 0.05 .* randn(rng, N_res)
N_res = 300
Wr = rand_sparse(rng, Float64, N_res, N_res; radius = 0.9, sparsity = 6 / N_res)
Win = scaled_rand(rng, Float64, N_res, 3)
bias = zeros(N_res)
initial_state = zeros(N_res)

# 3. Raw ODE equations — leaky-integrator continuous ESN
Expand Down Expand Up @@ -210,18 +211,27 @@ rc_predict = build_rc(predict_len)
ps, st = setup(rng, rc_train)

# 5. Fit the linear readout on the collected continuous states
ps, st = train(rc_train, input_data, target_data, ps, st)
ps, st = train(rc_train, input_data, target_data, ps, st;
objective = RidgeRegression(1.0e-6))

# 6. Autoregressive rollout under the same continuous dynamics
ps_pred, st_pred = setup(rng, rc_predict)
ps_pred = merge(ps_pred, (readout = ps.readout,))
st_pred = merge(st_pred, (readout = st.readout,))
output, _ = predict(rc_predict, predict_len, ps_pred, st_pred; initialdata = test[:, 1])

plot(transpose(output)[:, 1], transpose(output)[:, 2], transpose(output)[:, 3];
label = "predicted")
plot!(transpose(test)[:, 1], transpose(test)[:, 2], transpose(test)[:, 3];
label = "actual")
output, _ = predict(rc_predict, predict_len, ps, st; initialdata = test[:, 1])

dt = 0.02
lorenz_maxlyap = 0.9056
lyap_time = (0:(predict_len - 1)) .* dt .* (1 / lorenz_maxlyap)

p1 = plot(lyap_time, [test[1, :] output[1, :]]; label = ["actual" "predicted"],
ylabel = "x(t)", linewidth = 2.5, xticks = false, yticks = -15:15:15);
p2 = plot(lyap_time, [test[2, :] output[2, :]]; label = ["actual" "predicted"],
ylabel = "y(t)", linewidth = 2.5, xticks = false, yticks = -20:20:20);
p3 = plot(lyap_time, [test[3, :] output[3, :]]; label = ["actual" "predicted"],
ylabel = "z(t)", linewidth = 2.5, xlabel = "max(λ)*t", yticks = 10:15:40);

plot(p1, p2, p3; plot_title = "Lorenz System Coordinates",
layout = (3, 1), xtickfontsize = 12, ytickfontsize = 12, xguidefontsize = 15,
yguidefontsize = 15,
legendfontsize = 12, titlefontsize = 20)
```

The two trajectories should agree on the early portion of the rollout
Expand Down Expand Up @@ -311,10 +321,7 @@ ps_mg, st_mg = setup(rng, rc_mg_train)
ps_mg, st_mg = train(rc_mg_train, input_data, target_data, ps_mg, st_mg;
objective = RidgeRegression(1.0e-6), washout = 0)

ps_pred, st_pred = setup(rng, rc_mg_predict)
ps_pred = merge(ps_pred, (readout = ps_mg.readout,))
st_pred = merge(st_pred, (readout = st_mg.readout,))
mg_output, _ = predict(rc_mg_predict, predict_len, ps_pred, st_pred;
mg_output, _ = predict(rc_mg_predict, predict_len, ps_mg, st_mg;
initialdata = test_data[:, 1])

plot([test_data[1, :], mg_output[1, :]];
Expand Down