Tutorial 3: diagnosing a ql-ROM¶

Two questions, one notebook. First, is the model right? -- the two errors a ql-ROM can make, where they come from, and how to choose $(K, r)$ against them. Then, is it the right kind of dynamics? -- because a chaotic forecast stops shadowing the true state long before it stops being a useful model, and a pointwise error cannot tell a good statistical reproduction from a model that has collapsed.

The second half runs on ntsa: the full-order models through their dynamodels interface, and the fitted ql-ROMs from their own closed-loop records through ntsa's equation-free DataSeries -- signal-level analysis, attractor geometry, regime classification, a data-only Lyapunov exponent and the eight-panel diagnostic row. The cases are those of tutorial 0.

In [1]:
import os
import time
from pathlib import Path

os.environ.setdefault("QLROM_DATA", str(Path.home() / ".cache" / "qlrom"))

import matplotlib.pyplot as plt

from qlroms.intrusive_qlroms import ks1d
from qlroms.intrusive_qlroms.build_ks import load_case
from qlroms.utils.diagnosis import (
    compare_local_global_reconstruction,
    compute_projection_mse,
    default_forecast,
    sweep_qlrom_diagnosis,
)
from qlroms.utils.diagnosis_plots import (
    plot_cluster_occupancy,
    plot_error_timeseries,
    plot_fom_vs_roms,
    plot_kr_error_map,
    plot_parameter_selection,
    plot_transition_matrix,
)
from qlroms.utils.plots import close_pdf, start_pdf
# the figures only this tutorial draws (tutorials/plots_03.py)
from plots_03 import plot_local_global_errors

plt.rcParams["figure.dpi"] = 160

CASE = "chaotic"          # 1-D KS, nu = 1, L = 20 pi, Nx = 128, dt = 0.05
NTRAIN, NTEST = 40_000, 4_000
K_LIST = [2, 5, 10, 20]
R_LIST = [10, 20, 30, 40]
SEED = 0
t_start = time.time()

# cached trajectory, split into training and test windows (build_ks.load_case)
case_1d = load_case(ks1d, CASE, NTRAIN, NTEST)
fom, dt, save_dir = case_1d["fom"], case_1d["dt"], case_1d["save_dir"]
Xtrain, Xtest = case_1d["Xtrain"], case_1d["Xtest"]

print(f"1-D KS [{CASE}]: Nx={fom.Nx}, dt={dt}, train {tuple(Xtrain.shape)}, "
      f"test {tuple(Xtest.shape)} ({NTEST * dt:.0f} t.u.)")
print(f"models and diagnosis caches live in\n  {save_dir}")

# every figure below also becomes a page of figs/03_diagnosing_a_qlrom.pdf
start_pdf("03_diagnosing_a_qlrom")

/home/anovoama/.conda/envs/qlrom/lib/python3.11/site-packages/torch/cuda/__init__.py:188: UserWarning: CUDA initialization: The NVIDIA driver on your system is too old (found version 12080). Please update your GPU driver by downloading and installing a new version from the URL: http://www.nvidia.com/Download/index.aspx Alternatively, go to: https://pytorch.org to install a PyTorch version that has been compiled with your version of the CUDA driver. (Triggered internally at /__w/pytorch/pytorch/c10/cuda/CUDAFunctions.cpp:119.)
  return torch._C._cuda_getDeviceCount() > 0

Loading full trajectory from cache: /home/anovoama/.cache/qlrom/ks1d/nu_1.0000_Lx62.8319_Nx128/dt0.0500/full_trajectory.pth.
1-D KS [chaotic]: Nx=128, dt=0.05, train (128, 40000), test (128, 4000) (200 t.u.)
models and diagnosis caches live in
  /home/anovoama/.cache/qlrom/ks1d/nu_1.0000_Lx62.8319_Nx128/dt0.0500/Ntrain40000
Out[1]:
PosixPath('figs/03_diagnosing_a_qlrom.pdf')

1. The two errors a ql-ROM can make¶

A chart $k$ is an affine subspace: centroid $\bm{c}_k$, $M$-orthonormal basis $\mathsf{\Phi}_k$, and the projector

$$ \Pi_k \bm{u} \;=\; \bm{c}_k + \mathsf{\Phi}_k\,\mathsf{\Phi}_k^{\top} \mathsf{M} (\bm{u} - \bm{c}_k), \qquad \bm{a} \;=\; \mathsf{\Phi}_k^{\top}\mathsf{M}(\bm{u}-\bm{c}_k) . $$

A priori (representation). Project every test snapshot into its nearest chart and come straight back. No dynamics are involved, so the mean-squared error

$$ \mathrm{MSE}(K,r) \;=\; \frac{1}{M}\sum_{m=1}^{M} \big\| \bm{u}_m - \Pi_{k(m)}\bm{u}_m \big\|_{\mathsf{M}}^{2}, \qquad k(m) = \arg\min_k \|\bm{u}_m - \bm{c}_k\| , $$

measures the charts alone. For $K=1$ this is the ordinary POD projection error. A forecast lives on the charts, so a chart set that cannot represent the data will not be rescued by any dynamics.

A posteriori (forecast). Advance the reduced state in closed loop from a single initial condition and compare with the true state,

$$ \varepsilon(t) \;=\; \frac{\|\bm{u}(t) - \tilde{\bm{u}}(t)\|}{\|\bm{u}(t)\|}\times 100\% . $$

The gap between the two is the whole subject: the first is bounded by the geometry, the second adds the model error and, on a chaotic case, the exponential divergence that no ROM can avoid.

In [2]:
K, r = 5, 40
local_rom = ks1d.build_local_model(Xtrain, fom, r=r, K=K, save_dir=save_dir,
                                   clustering_kwargs={"random_state": SEED})
global_rom = ks1d.build_local_model(Xtrain, fom, r=r, K=1, save_dir=save_dir,
                                    clustering_kwargs={"random_state": SEED})

mse_local = compute_projection_mse(Xtest, local_rom)
mse_global = compute_projection_mse(Xtest, global_rom)
print(f"projection MSE on the test set (r={r})")
print(f"  global (K=1) : {mse_global:.4e}")
print(f"  local  (K={K})  : {mse_local:.4e}   ->  ratio global/local = "
      f"{mse_global / mse_local:.3f}  ({'local better' if mse_local < mse_global else 'global better'})")

Cached model at /home/anovoama/.cache/qlrom/ks1d/nu_1.0000_Lx62.8319_Nx128/dt0.0500/Ntrain40000/local_model_K5.pth was not fitted on these training snapshots (same count but other steps, or saved before this check). Rebuilding model.

Fitted K=5 clusters with sizes: [7748, 8395, 8037, 8730, 7090].

Saved qlROM cache: K=5, r_store=40, file=dt0.0500/Ntrain40000/local_model_K5.pth.
Cached model at /home/anovoama/.cache/qlrom/ks1d/nu_1.0000_Lx62.8319_Nx128/dt0.0500/Ntrain40000/local_model_K1.pth was not fitted on these training snapshots (same count but other steps, or saved before this check). Rebuilding model.
Fitted K=1 clusters with sizes: [40000].

Saved qlROM cache: K=1, r_store=40, file=dt0.0500/Ntrain40000/local_model_K1.pth.
projection MSE on the test set (r=40)
  global (K=1) : 1.3569e-01
  local  (K=5)  : 1.2615e-01   ->  ratio global/local = 1.076  (local better)

The ratio global/local is the verdict of the a priori test: above 1, the $K$ charts represent unseen snapshots better than one global basis of the same size $r$; below 1, worse. A ratio below 1 is not a bug: the diagnosis is doing its job. When $r$ is a large fraction of the $N_h$ degrees of freedom, the global basis already captures most of the variance, while each local basis is fitted on about $1/K$ of the snapshots and generalizes worse off its own cluster. Section 5 maps the ratio over the $(K,r)$ plane.

2. Both at once: compare_local_global_reconstruction¶

One call runs the a priori and the a posteriori comparison for a given $(K, r)$ and returns both trajectories, so they can be looked at as fields rather than as scalars. It prints a PASS/FAIL line, which is only shorthand for "is the local RMSE below the global one". The figure shows both errors in time on a shared axis, the representation on the left and the forecast on the right, so the difference between the two panels is the error that the dynamics add.

In [3]:
res = compare_local_global_reconstruction(
    Xtest=Xtest, FOM=fom, local_rom=local_rom, global_rom=global_rom, save_dir=save_dir,
)

plot_local_global_errors(res, dt, K=K, r=r)
plt.show()

print(f"representation RMSE:  local {res['local_rep_rmse']:.4e}   global {res['global_rep_rmse']:.4e}")
print(f"forecast       RMSE:  local {res['local_rmse']:.4e}   global {res['global_rmse']:.4e}")
# part-one numbers for the SUMMARY: `res`, `K` and `r` are reused as loop names in part two
part1 = dict(K=K, r=r, ratio=mse_global / mse_local, rmse_local=res["local_rmse"], rmse_global=res["global_rmse"])

Loaded compare cache (forecast only): /home/anovoama/.cache/qlrom/ks1d/nu_1.0000_Lx62.8319_Nx128/dt0.0500/Ntrain40000/compare_K5_r40_Ntrain40000_Ntest4000.npz


 Computing representation (project-recover every snapshot)...

Local-vs-global reconstruction sanity check
  forecast:       local (K=5, r=40) RMSE=1.664310   global RMSE=1.792657  PASS
  representation: local (K=5, r=40) RMSE=0.031393   global RMSE=0.032559 PASS
No description has been provided for this image 
representation RMSE:  local 3.1393e-02   global 3.2559e-02
forecast       RMSE:  local 1.6643e+00   global 1.7927e+00

3. The diagnosis does not know the model family¶

Nothing above stepped anything itself. compare_local_global_reconstruction and sweep_qlrom_diagnosis reach a model through one hook,

forecast_fn(model, x0, n_steps) -> (Ndof, n_steps)      # physical trajectory

whose default is qlroms.utils.diagnosis.default_forecast: project, step, recover, through the common model interface, which every qlROM compilation satisfies as it is. A model that does not step reduced coordinates — one that carries its own hidden state — is diagnosed by passing its own closure instead; nothing else changes, because everything measured lives on the charts.

In [4]:
import inspect  # noqa: E402

print(inspect.getsource(default_forecast).split('"""')[0].strip())

def default_forecast(model, x0, n_steps: int) -> torch.Tensor:

That separation — geometry on one side (the charts), dynamics on the other (whatever forecast_fn steps) — is the whole point of running both halves of the diagnosis.

4. Choosing $K$: BIC over the clustering¶

The clustering is a model too, so it can be scored like one. For hard $k$-means under an isotropic Gaussian,

$$ \mathrm{BIC}(K) \;=\; \nu \log M \;-\; 2\hat{\ell}, \qquad \nu = K \cdot N_h, \qquad \hat{\sigma}^2 = \frac{J}{N_h M}, $$

with $J$ the within-cluster sum of squares. On snapshot data, BIC tends to keep falling as $K$ grows, so the useful statistic is where it stops falling fast — the elbow of $\Delta\mathrm{BIC}/\Delta K$, which is what sweep_k_bic and the sweep below report as $K_{\mathrm{opt}}$.

5. The $(K, r)$ sweep¶

sweep_qlrom_diagnosis runs the whole grid: for each $(K,r)$ it builds the model, measures the projection MSE on train and test, runs the closed loop, and scores the cluster occupancy of that run against the data. Results are cached in save_dir, so re-running with more $K$ or $r$ extends the file instead of redoing it, and the sweep cell is instant once runs/sweep.sh has filled the cache for the same case and windows.

In [5]:
diag = sweep_qlrom_diagnosis(
    Xtrain=Xtrain, Xtest=Xtest, FOM=fom, K_list=K_LIST, r_list=R_LIST,
    build_model_fn=ks1d.build_local_model, save_dir=save_dir,
)
print(f"\nK_opt = {diag['K_opt']}, r_check = {diag['r_check']}")

Loaded diagnosis cache: /home/anovoama/.cache/qlrom/ks1d/nu_1.0000_Lx62.8319_Nx128/dt0.0500/Ntrain40000/sweep_diagnosis_Ntrain40000_Ntest4000_dt0.05.npz
  BIC: K=2  BIC=4273483.26 (cached)
  BIC: K=5  BIC=4175442.78 (cached)
  BIC: K=10  BIC=4036326.64 (cached)
  BIC: K=20  BIC=3974837.28 (cached)
Cached model at /home/anovoama/.cache/qlrom/ks1d/nu_1.0000_Lx62.8319_Nx128/dt0.0500/Ntrain40000/local_model_K20.pth was not fitted on these training snapshots (same count but other steps, or saved before this check). Rebuilding model.

Fitted K=20 clusters with sizes: [1938, 1835, 2585, 3356, 2718, 1444, 1918, 1578, 1473, 1797, 1919, 1865, 2032, 1471, 2051, 2350, 1288, 2117, 1906, 2359].

Saved qlROM cache: K=20, r_store=40, file=dt0.0500/Ntrain40000/local_model_K20.pth.


Sweeping global (K=1) ROM...
  K=1, r= 40  err=1.7927  proj_mse_te=1.3666e-01  proj_mse_tr=1.1161e-01  kld=6.3196e+00  (cached)
  K=1, r= 30  err=2.3958  proj_mse_te=2.8140e+00  proj_mse_tr=2.0137e+00  kld=4.9437e-01  (cached)
  K=1, r= 20  err=15.3751  proj_mse_te=2.2019e+01  proj_mse_tr=1.8259e+01  kld=2.0344e-01  (cached)

  K=1, r= 10  err=nan  proj_mse_te=6.2516e+01  proj_mse_tr=5.4419e+01  kld=inf
Sweeping local ROM K=2...
  K=2, r= 40  err=1.6688  proj_mse_te=1.3043e-01  proj_mse_tr=1.1588e-01  kld=2.4999e+00  (cached)
  K=2, r= 30  err=2.2223  proj_mse_te=2.6952e+00  proj_mse_tr=1.9823e+00  kld=5.4096e-01  (cached)
  K=2, r= 20  err=3.4390  proj_mse_te=1.9655e+01  proj_mse_tr=1.6997e+01  kld=7.7441e-01  (cached)
  K=2, r= 10  err=7.6025  proj_mse_te=6.2545e+01  proj_mse_tr=5.4924e+01  kld=1.7463e+01  (cached)
Sweeping local ROM K=5...
  K=5, r= 40  err=1.6643  proj_mse_te=1.4692e-01  proj_mse_tr=1.2000e-01  kld=1.0653e+01  (cached)
  K=5, r= 30  err=1.9936  proj_mse_te=2.3514e+00  proj_mse_tr=2.5109e+00  kld=4.2959e+00  (cached)
  K=5, r= 20  err=2.4005  proj_mse_te=2.0048e+01  proj_mse_tr=1.6217e+01  kld=9.5648e-01  (cached)
  K=5, r= 10  err=2.0473  proj_mse_te=5.7335e+01  proj_mse_tr=5.2667e+01  kld=1.3073e+01  (cached)
Sweeping local ROM K=10...
  K=10, r= 40  err=1.7638  proj_mse_te=1.9034e-01  proj_mse_tr=1.9447e-01  kld=1.9348e+00  (cached)
  K=10, r= 30  err=1.9045  proj_mse_te=3.4288e+00  proj_mse_tr=2.7397e+00  kld=4.4501e-01  (cached)
  K=10, r= 20  err=1.8814  proj_mse_te=2.0059e+01  proj_mse_tr=1.6426e+01  kld=3.6247e+00  (cached)
  K=10, r= 10  err=2.0569  proj_mse_te=6.6800e+01  proj_mse_tr=5.5256e+01  kld=7.0193e+00  (cached)
Sweeping local ROM K=20...
  K=20, r= 40  err=1.8261  proj_mse_te=5.5806e-01  proj_mse_tr=5.2871e-01  kld=1.3045e+01  (cached)
  K=20, r= 30  err=1.8016  proj_mse_te=5.9067e+00  proj_mse_tr=4.7753e+00  kld=5.6895e+00  (cached)
  K=20, r= 20  err=1.8994  proj_mse_te=2.3460e+01  proj_mse_tr=2.0148e+01  kld=5.1060e+00  (cached)
  K=20, r= 10  err=1.9304  proj_mse_te=6.4922e+01  proj_mse_tr=5.9272e+01  kld=9.3152e+00  (cached)
  Loaded cluster probs/transition matrices from cache (K_opt=5).
  KL(P_test  || P_train)      = 6.51e-02  (cached)
  KL(P_test  || P_qlrom/test) = 2.77e-01  (cached)
  KL(P_test  || P_grom /test) = 1.27e-01  (cached)
  KL(P_test  || P_qlrom/repr) = 3.63e-06  (cached)
  KL(P_test  || P_grom /repr) = 0.00e+00  (cached)

K_opt = 5, r_check = 40
In [6]:
fig, _ = plot_kr_error_map(diag)
plt.show()
No description has been provided for this image

Read the maps left to right; the bottom row, below the white line, is the global ROM ($K=1$). The forecast error is the a posteriori score of each closed loop, and the star marks the lowest. The projection MSE is the a priori score: how it changes along $r$ and along $K$ tells whether the modes or the charts do the work. The ratio global/local is the a priori verdict: green is a ql-ROM win, red a loss, and the $K=1$ row is 1 by construction. Where it turns red, the local bases represent unseen snapshots worse than one global basis of the same size, and the next figure shows whether the training fit degrades too. The KLD map scores the cluster occupancy of the closed-loop forecast against the data; it is the only panel that looks at where the trajectory goes rather than how far it is at each instant.

In [7]:
fig, _ = plot_parameter_selection(diag)
plt.show()
No description has been provided for this image

Left: the BIC elbow. Middle and right: projection MSE against $K$ (at $r=r_{\rm check}$) and against $r$ (at $K=K_{\rm opt}$), with train as crosses and test as circles. The vertical gap between them measures how well the charts generalize: a test error that stops improving while the training error keeps improving is overfitting.

6. Statistics: does the model live where the data lives?¶

Error curves die with the trajectory: on a chaotic case every ROM decorrelates, and after that $\varepsilon(t)$ says nothing except "different". Two statistics survive. The cluster occupancy $P(\bm{c}_k)$ — the fraction of time spent in each chart — compared with the data's own, through

$$ D_{\mathrm{KL}}(P_{\mathrm{data}} \,\|\, P_{\mathrm{ROM}}) \;=\; \sum_k P_{\mathrm{data}}(\bm{c}_k)\, \log \frac{P_{\mathrm{data}}(\bm{c}_k)}{P_{\mathrm{ROM}}(\bm{c}_k)} , $$

and the Markov transition matrix $T_{ij} = P(\bm{c}_j \mid \bm{c}_i)$ of the chart sequence, which adds the order in which charts are visited.

In [8]:
fig, _ = plot_cluster_occupancy(diag)
plt.show()
fig, _ = plot_transition_matrix(diag)
plt.show()
No description has been provided for this image
No description has been provided for this image

The left occupancy panel is the a priori one. Projecting and recovering a true snapshot moves it by its representation error, so only snapshots near a chart boundary can change chart, and both ROMs should reproduce the data histogram closely, with a $D_{\mathrm{KL}}$ close to zero (the legend quotes it). The right panel is the closed-loop forecast: a long forecast on a chaotic attractor can redistribute its time among the charts, and $D_{\mathrm{KL}}$ measures by how much. The transition matrices (top row a priori, bottom row forecast) are diagonally dominant whenever the dwell time is long compared with $\Delta t$, because a chart is left after many steps, not at the next one; the off-diagonal entries carry the order in which the charts are visited.

7. Error time series across $K$¶

In [9]:
fig, _ = plot_error_timeseries(diag, dt=dt)
plt.show()
No description has been provided for this image

8. The physical page¶

The scalars above are summaries; the two figures below are the thing itself. Snapshots at four instants, the spatiotemporal map, four probe time series and the pointwise PDF $p(u)$ — FOM, global ROM and ql-ROM on the same axes, for the representation (first figure) and the forecast (second figure). This is the page runs/diagnosis.sh writes to a PDF for a single $(K,r)$.

In [10]:
fig = plot_fom_vs_roms(Xtest[:, :1500], res["Xlocal_rep"][:, :1500], res["Xglobal_rep"][:, :1500],
                       fom, K=K, r=r, figsize=(12, 12))
fig.suptitle(f"representation: FOM vs ROMs  (K={K}, r={r})", fontsize=11)
plt.show()
No description has been provided for this image
In [11]:
fig = plot_fom_vs_roms(Xtest[:, :1500], res["Xlocal_rec"][:, :1500], res["Xglobal_rec"][:, :1500],
                       fom, K=K, r=r, figsize=(12, 12))
fig.suptitle(f"forecast: FOM vs ROMs  (K={K}, r={r})", fontsize=11)
plt.show()
No description has been provided for this image

9. Running it offline¶

The two jobs above are also the two scripts, so a long sweep does not have to live in a notebook:

bash runs/diagnosis.sh    # one (K, r): scripts/ks/run_diagnosis.py -> diagnosis_<family>_K<K>_r<r>.pdf
bash runs/sweep.sh        # the grid:   scripts/ks/run_sweep.py     -> diagnosis_<family>.pdf + .npz

Both take --1d/--2d, --case and --Ntrain/--Ntest, and run_sweep.py also takes the clustering flags. runs/sweep.sh runs under nohup and tails its log, because a full grid is a long job; everything it computes lands in the .npz this notebook reads, so the figures above are a cache hit away.

Summary of part one: is the model right?¶

question metric where it lives
can the charts represent the data? projection MSE, ratio vs $K=1$ compute_projection_mse, map panels 2-3
how many charts? $\Delta\mathrm{BIC}/\Delta K$ elbow sweep_k_bic, selection panel 1
how many modes? MSE vs $r$, train and test selection panel 3
does the model predict? $\varepsilon(t)$, forecast RMSE compare_local_global_reconstruction, map panel 1
does it live in the right place? occupancy $P(\bm{c}_k)$, $D_{KL}$, $T_{ij}$ occupancy / transition figures

Every row is computed from charts and trajectories only. The model enters once, as forecast_fn, which is why the same diagnosis serves any qlROM without a line of model-specific code.

10. Does it keep the character of the system?¶

Everything above measures how far the model is from the true state. Past the predictability time of a chaotic system, that distance saturates and stops discriminating: a model that has drifted onto the attractor and a model that has collapsed to a fixed point can score the same.

Tutorial 0 characterizes the cases -- spectra, attractor geometry, regime, leading Lyapunov exponent. Here the same measurements are run on the fitted ql-ROMs and laid next to them; the ql-ROMs are measured from their closed-loop records, as if they were data. The cells below rebuild the case baselines with the helpers of tutorial 0, whose method is not repeated here.

In [12]:
import os
import time
from pathlib import Path

os.environ.setdefault("QLROM_DATA", str(Path.home() / ".cache" / "qlrom"))

import matplotlib.pyplot as plt
import numpy as np
from ntsa import lyapunov as lyap
from ntsa import tools as nt
from ntsa.data import DataSeries

from qlroms import free_run
from qlroms.utils.characterization import (
    load_characterization_case,
    lyapunov_verdict,
    saturation_rate,
    signal_report,
    spectrum_text,
)
from qlroms.utils.diagnosis import closed_loop_forecast
from qlroms.utils.plots import outputs_dir, save_figure
# the figures only this tutorial draws (tutorials/plots_03.py)
from plots_03 import plot_attractor_overlay, plot_delay_portraits, plot_psd_comparison

plt.rcParams["figure.dpi"] = 160

SEED = 0

N_SIG = 12_000        # samples fed to every signal-level diagnostic
N_FC = 2_000          # closed-loop steps used only for the prediction-horizon contrast
EXP_WINDOW_MIN = 0.15  # a fitted lambda_1 counts only if its window spans this much of the record

# Each case brings its own (K, r) and its own characterization timescales from
# its module's TEST_CASES -- t_lam is the FOM Lyapunov horizon, which must cover the growth
# phase and stop near saturation. Only the WINDOWS are this notebook's own: the cases are
# generated far longer than a notebook of about 15 min can characterize.
WINDOWS = {"ks1d qp":    dict(Ntrain=10_000, Ntest=2_000),
           "ks1d chaos": dict(Ntrain=20_000, Ntest=2_000),
           "ks2d tw":    dict(Ntrain=15_000, Ntest=2_000),
           "ks2d chaos": dict(Ntrain=40_000, Ntest=2_000)}

_OUT = outputs_dir()          # tutorials/figs, from the repo root or from tutorials/
In [13]:
# cached trajectory, sensor signal and respawned FOM of every case (load_characterization_case)
cases = {name: load_characterization_case(name, n_sig=N_SIG, **window) for name, window in WINDOWS.items()}
for c in cases.values():
    n_h = c["traj"].shape[0]
    print(f"{c['name']:11s} ({c['P']['case']:14s}): N_h = {n_h:5d}, dt = {c['dt']:.3f}, "
          f"{c['traj'].shape[1]} cached snapshots | sensor row {c['obs_idx'][c['i_obs']]} "
          f"| analysis window {N_SIG} samples = {N_SIG * c['dt']:.0f} t.u. | recorded lamb1 = {c['lamb1_ref']}")

Loading full trajectory from cache: /home/anovoama/.cache/qlrom/ks1d/nu_0.2254_Lx6.2832_Nx128/dt0.1000/full_trajectory.pth.
Loading full trajectory from cache: /home/anovoama/.cache/qlrom/ks1d/nu_1.0000_Lx62.8319_Nx128/dt0.0500/full_trajectory.pth.
Loading full trajectory from cache: /home/anovoama/.cache/qlrom/ks2d/nu1_0.50_nu2_0.35_Nx32_Ny32/dt0.1000/full_trajectory.pth.

Loading full trajectory from cache: /home/anovoama/.cache/qlrom/ks2d/nu1_0.30_nu2_0.10_Nx64_Ny64/dt0.0100/full_trajectory.pth.

ks1d qp     (quasi-periodic): N_h =   128, dt = 0.100, 27000 cached snapshots | sensor row 42 | analysis window 12000 samples = 1200 t.u. | recorded lamb1 = 0.0
ks1d chaos  (chaotic       ): N_h =   128, dt = 0.050, 52000 cached snapshots | sensor row 42 | analysis window 12000 samples = 600 t.u. | recorded lamb1 = 0.062
ks2d tw     (travelling    ): N_h =  1024, dt = 0.100, 27000 cached snapshots | sensor row 341 | analysis window 12000 samples = 1200 t.u. | recorded lamb1 = 0.0
ks2d chaos  (chaotic_B     ): N_h =  4096, dt = 0.010, 52000 cached snapshots | sensor row 1365 | analysis window 12000 samples = 120 t.u. | recorded lamb1 = 1.8
In [14]:
# baselines: the tutorial-0 measurements for the FOM side, recomputed here with the same
# qlroms.utils.characterization helpers so this notebook stands alone -- see tutorial 0
# for what they mean
t0 = time.time()
for c in cases.values():
    c["sig"] = signal_report(c["x_fom"], c["dt"], f"{c['name']} FOM")
    s = c["sig"]
    print(f"{c['name']:11s}: PSD {spectrum_text(s)}; tallest f = {s['f_top']:.4f}, centroid "
          f"{s['f_c']:.4f} | ACF first zero {s['t_dec']:.1f} t.u., tail |rho| {s['acf_tail']:.2f} "
          f"| optimal lag zeta = {s['zeta']} samples = {s['zeta'] * c['dt']:.2f} t.u. | FNN "
          f"embedding dim = {s['dim']} | D2 = {s['D2']:.2f} (small-r {s['D2_lo']:.2f}, "
          f"large-r {s['D2_hi']:.2f})")
print(f"  (all {len(cases)} signal reports in {time.time() - t0:.1f} s)")

ks1d qp    : PSD  6 prominent lines, the strongest 6 on 2 generator(s) [0.0208 0.0658]; tallest f = 0.0208, centroid 0.3295 | ACF first zero 12.1 t.u., tail |rho| 0.69 | optimal lag zeta = 119 samples = 11.90 t.u. | FNN embedding dim = 10 | D2 = 0.50 (small-r 0.72, large-r 0.51)

ks1d chaos : PSD broadband, 18 periodogram maxima, no line lattice; tallest f = 0.0067, centroid 0.3996 | ACF first zero 18.2 t.u., tail |rho| 0.23 | optimal lag zeta = 77 samples = 3.85 t.u. | FNN embedding dim = 3 | D2 = 2.01 (small-r 2.11, large-r 1.71)

ks2d tw    : PSD  7 prominent lines, the strongest 7 on 2 generator(s) [0.0175 0.03  ]; tallest f = 0.0300, centroid 0.1777 | ACF first zero 6.7 t.u., tail |rho| 0.81 | optimal lag zeta = 57 samples = 5.70 t.u. | FNN embedding dim = 4 | D2 = 2.30 (small-r 2.10, large-r 2.48)

ks2d chaos : PSD broadband, 18 periodogram maxima, no line lattice; tallest f = 0.0750, centroid 3.0566 | ACF first zero 2.9 t.u., tail |rho| 0.16 | optimal lag zeta = 85 samples = 0.85 t.u. | FNN embedding dim = 3 | D2 = 2.02 (small-r 2.28, large-r 1.77)
  (all 4 signal reports in 3.6 s)
In [15]:
t0 = time.time()
for c in cases.values():
    P = c["P"]
    lam1, lam1_std, res = lyap.leading_lyapunov(c["fom_run"], n_pert=3, t_run=P["t_lam"])
    lam_ok, lam_str, frac = lyapunov_verdict(lam1, res, min_window=EXP_WINDOW_MIN)
    climb, span, rate = saturation_rate(res)
    c["lyap_fom"] = dict(lam1=lam1, std=lam1_std, res=res, ok=lam_ok, text=lam_str, frac=frac,
                         climb=climb, span=span, rate=rate)
    ref = c["lamb1_ref"]
    sat_txt = (f"climbs {climb:4.1f} nats in {span:5.1f} t.u., mean rate {rate:.3f}/t.u."
               if lam_ok is not None else
               f"plateaus at {100 * res['sat'] / res['diam']:5.2f}% of the attractor diameter")
    # the recorded exponent as a TIME: how long exp(lamb1 t) needs to cover the same climb
    ref_txt = (f" (at the recorded {ref}, that climb takes {climb / ref:.1f} t.u.)"
               if lam_ok is not None and ref > 0 else "")
    print(f"{c['name']:11s}: lam1 = {lam_str:22s} (raw fit {lam1:.3f}, R2 = {res['r2']:.2f}, "
          f"window {100 * frac:4.1f}% of {P['t_lam']:.0f} t.u.) | {sat_txt}{ref_txt} | "
          f"config records {ref}")
print(f"  ({len(cases)} leading_lyapunov runs in {time.time() - t0:.0f} s)")

ks1d qp    : lam1 = transient (9% win)     (raw fit 0.299, R2 = 0.72, window  9.1% of 200 t.u.) | plateaus at 13.40% of the attractor diameter | config records 0.0

ks1d chaos : lam1 = 0.067                  (raw fit 0.067, R2 = 0.98, window 45.6% of 250 t.u.) | climbs 13.0 nats in 197.5 t.u., mean rate 0.066/t.u. (at the recorded 0.062, that climb takes 210.1 t.u.) | config records 0.062

ks2d tw    : lam1 = rejected (no growth)   (raw fit nan, R2 = 0.00, window  1.7% of 120 t.u.) | plateaus at  0.00% of the attractor diameter | config records 0.0

ks2d chaos : lam1 = 0.303                  (raw fit 0.303, R2 = 0.95, window 39.2% of 90 t.u.) | climbs 15.1 nats in  53.0 t.u., mean rate 0.286/t.u. (at the recorded 1.8, that climb takes 8.4 t.u.) | config records 1.8
  (4 leading_lyapunov runs in 62 s)
In [16]:
rows = []
for c in cases.values():
    L = c["lyap_fom"]
    label, ev = nt.classify_regime(c["x_fom"], c["dt"], lam1=L["ok"],
                                   lam1_std=(L["std"] if L["ok"] is not None else 0.0),
                                   t_total=N_SIG * c["dt"])
    c["regime_fom"], c["evidence_fom"] = label, ev
    ratio = f"{ev['f2'] / ev['f1']:.3f}" if (ev["f1"] and ev["f2"]) else "-"
    rows.append((c["name"], L["text"], c["lamb1_ref"], label, ev["n_clusters"], c["sig"]["dim"],
                 c["sig"]["D2"], len(c["sig"]["gen"]) if c["sig"]["discrete"] else 0,
                 ratio, ev["rational_match"] or "-"))

# the last three columns are the two spectral verdicts side by side: OUR generator count
# (0 = broadband, no lattice quoted) and the f2/f1 ratio classify_regime forms internally
hdr = (f"{'case':<12} {'lam1 measured':>22} {'lamb1 cfg':>10} {'classify_regime':>18} "
       f"{'clusters':>9} {'d_FNN':>6} {'D2':>6} {'n_gen':>6} {'its f2/f1':>10} {'rational':>9}")
print(hdr)
print("-" * len(hdr))
for name, lam_s, lam_ref, label, ncl, dim, D2, ngen, ratio, rat in rows:
    print(f"{name:<12} {lam_s:>22} {lam_ref:>10.3f} {label:>18} {ncl:>9d} {dim:>6d} {D2:>6.2f} "
          f"{ngen:>6d} {ratio:>10} {rat:>9}")

case                  lam1 measured  lamb1 cfg    classify_regime  clusters  d_FNN     D2  n_gen  its f2/f1  rational
---------------------------------------------------------------------------------------------------------------------
ks1d qp          transient (9% win)      0.000 limit_cycle_period_3         3     10   0.50      2      3.160         -
ks1d chaos                    0.067      0.062            chaotic        35      3   2.01      0      2.500       5/2
ks2d tw        rejected (no growth)      0.000      quasiperiodic        59      4   2.30      2      1.583         -
ks2d chaos                    0.303      1.800            chaotic        32      3   2.02      0      1.333       4/3

11. The payoff: the same analysis on a fitted qlROM¶

Now build one quantized-local ROM per case with the intrusive (Galerkin) route of tutorial 2 and run it in closed loop. The qlROM side is characterized from its closed-loop record, through ntsa's equation-free path: free_run produces the physical trajectory, the FOM's sensor row is read off it, and every diagnostic below treats that series exactly as it would treat measured data. The only measurement that changes is the Lyapunov exponent: the FOM's comes from leading_lyapunov, which perturbs the model and follows the separation, while the qlROM's is Rosenstein's data-only estimate from its record. A model-based exponent for the qlROM needs it wrapped as a dynamodels model, and that adapter lives with the data-assimilation code in the qlROM-DA repository, not here. To keep the comparison like for like, the Rosenstein estimate is also computed on the FOM record.

First, the contrast that motivates the whole notebook. A closed-loop forecast of N_FC steps is scored by its predictability time, $T_{\mathrm{ph}}$, the time at which the relative error first exceeds 0.5 -- a trajectory measure. Then the ROM is integrated for N_SIG steps with no reference at all, and the signal diagnostics of section 10 are recomputed on its sensor signal: an attractor measure.

One reporting rule. prediction_horizon returns the end of the window when the threshold is never crossed, and that number is the window length, not a horizon. Those entries are printed as >= T (censored), with $T$ the window length: the run was still tracking when we stopped looking, and the only honest statement is a lower bound.

In [17]:
def tph_text(c, width=0):
    """T_ph, or the bound it really is when the error never crossed the threshold."""
    txt = f">= {c['T_ph']:.1f} (censored)" if c["T_ph_cens"] else f"{c['T_ph']:.1f}"
    return f"{txt:>{width}}" if width else txt


t0 = time.time()
for c in cases.values():
    P, mod = c["P"], c["P"]["module"]
    tb = time.time()
    c["rom"] = mod.build_local_model(c["Xtrain"].clone(), c["fom"], r=P["r"], K=P["K"],
                                     save_dir=c["save_dir"], clustering_kwargs={"random_state": SEED})
    t_build = time.time() - tb
    # trajectory-wise score, for contrast only (tutorial 2's truncation rule): column n is
    # time n, from x0 as the model represents it, cut at the first non-finite state or the
    # first state above twice the norm of the true state
    run = closed_loop_forecast(c["rom"], c["Xtest"], N_FC)
    c["T_ph"] = run["horizon"]
    # prediction_horizon returns the END of the window when the threshold is never crossed:
    # that is a censored observation (T_ph >= window), not a measured horizon. Flag it.
    c["T_ph_cens"] = bool(len(run["err"]) and run["err"].max() <= 0.5)
    print(f"{c['name']:11s}: qlGalerkin K={P['K']:2d} r={P['r']:2d} built in {t_build:4.1f} s | "
          f"forecast valid "
          f"{run['n_valid']}/{N_FC} steps, T_ph(0.5) = {tph_text(c)} t.u.")
print(f"  (builds + forecasts in {time.time() - t0:.0f} s)")

Cached model at /home/anovoama/.cache/qlrom/ks1d/nu_0.2254_Lx6.2832_Nx128/dt0.1000/Ntrain10000/local_model_K6.pth was not fitted on these training snapshots (same count but other steps, or saved before this check). Rebuilding model.

Fitted K=6 clusters with sizes: [4697, 4772, 102, 136, 162, 131].
Saved qlROM cache: K=6, r_store=30, file=dt0.1000/Ntrain10000/local_model_K6.pth.

ks1d qp    : qlGalerkin K= 6 r=30 built in  0.4 s | forecast valid 2000/2000 steps, T_ph(0.5) = 14.7 t.u.
Cached model at /home/anovoama/.cache/qlrom/ks1d/nu_1.0000_Lx62.8319_Nx128/dt0.0500/Ntrain20000/local_model_K10.pth was not fitted on these training snapshots (same count but other steps, or saved before this check). Rebuilding model.

Fitted K=10 clusters with sizes: [2361, 1483, 1546, 2414, 1382, 3541, 1904, 1719, 1756, 1894].

Saved qlROM cache: K=10, r_store=30, file=dt0.0500/Ntrain20000/local_model_K10.pth.

ks1d chaos : qlGalerkin K=10 r=30 built in  0.9 s | forecast valid 2000/2000 steps, T_ph(0.5) = 3.7 t.u.
Failed to load cached model (will rebuild): Cached model at /home/anovoama/.cache/qlrom/ks2d/nu1_0.50_nu2_0.35_Nx32_Ny32/dt0.1000/Ntrain15000/local_model_K5.pth was not fitted on these training snapshots (same count but other steps, or saved before this check).
Building qlROM with K=5 clusters... 
 and r=50 (requested r=25).

Saved qlROM cache: K=5, r_store=50, file=dt0.1000/Ntrain15000/local_model_K5.pth.

ks2d tw    : qlGalerkin K= 5 r=25 built in  1.3 s | forecast valid 2000/2000 steps, T_ph(0.5) = >= 200.0 (censored) t.u.

Failed to load cached model (will rebuild): Cached model at /home/anovoama/.cache/qlrom/ks2d/nu1_0.30_nu2_0.10_Nx64_Ny64/dt0.0100/Ntrain40000/local_model_K40.pth was not fitted on these training snapshots (same count but other steps, or saved before this check).
Building qlROM with K=40 clusters... 
 and r=50 (requested r=50).

Saved qlROM cache: K=40, r_store=50, file=dt0.0100/Ntrain40000/local_model_K40.pth.

ks2d chaos : qlGalerkin K=40 r=50 built in 30.1 s | forecast valid 1616/2000 steps, T_ph(0.5) = 1.8 t.u.
  (builds + forecasts in 48 s)
In [18]:
t0 = time.time()
for c in cases.values():
    # the qlROM's own closed-loop record: N_SIG steps from the first test snapshot, read at
    # the FOM's sensor row -- from here on the ROM is a data series, not a model
    X_rom, _ = free_run(c["rom"], c["Xtest"][:, 0], n_steps=N_SIG)
    c["X_rom"] = X_rom.cpu().numpy()
    if not np.isfinite(c["X_rom"]).all():
        raise RuntimeError(f"{c['name']}: the qlROM closed-loop forecast went non-finite")
    c["x_rom"] = c["X_rom"][c["obs_idx"][c["i_obs"]]]
    c["sig_rom"] = signal_report(c["x_rom"], c["dt"], f"{c['name']} ROM")
    s, sr = c["sig"], c["sig_rom"]
    print(f"{c['name']:11s}: FOM std {s['x'].std():.3f} -> ROM {sr['x'].std():.3f} "
          f"({100 * (sr['x'].std() / s['x'].std() - 1):+5.1f}%) | zeta {s['zeta']:3d} -> "
          f"{sr['zeta']:3d} | d_FNN {s['dim']} -> {sr['dim']} | D2 {s['D2']:.2f} -> "
          f"{sr['D2']:.2f} | tallest f {s['f_top']:.4f} -> {sr['f_top']:.4f} | f_centroid "
          f"{s['f_c']:.4f} -> {sr['f_c']:.4f} | PSD {spectrum_text(sr)}")
print(f"  ({N_SIG} closed-loop ROM steps per case in {time.time() - t0:.0f} s)")

ks1d qp    : FOM std 3.287 -> ROM 3.178 ( -3.3%) | zeta 119 ->  28 | d_FNN 10 -> 10 | D2 0.50 -> 1.73 | tallest f 0.0208 -> 0.0650 | f_centroid 0.3295 -> 0.4042 | PSD 20 prominent lines, the strongest 8 on 4 generator(s) [0.0583 0.0617 0.0633 0.065 ]

ks1d chaos : FOM std 1.363 -> ROM 1.267 ( -7.1%) | zeta  77 ->  69 | d_FNN 3 -> 4 | D2 2.01 -> 2.51 | tallest f 0.0067 -> 0.0017 | f_centroid 0.3996 -> 0.5548 | PSD broadband, 28 periodogram maxima, no line lattice

ks2d tw    : FOM std 2.747 -> ROM 2.731 ( -0.6%) | zeta  57 ->  50 | d_FNN 4 -> 4 | D2 2.30 -> 2.34 | tallest f 0.0300 -> 0.0475 | f_centroid 0.1777 -> 0.4472 | PSD  9 prominent lines, the strongest 8 on 2 generator(s) [0.0175 0.03  ]

ks2d chaos : FOM std 2.709 -> ROM 3.512 (+29.7%) | zeta  85 ->  79 | d_FNN 3 -> 4 | D2 2.02 -> 2.73 | tallest f 0.0750 -> 0.0583 | f_centroid 3.0566 -> 4.3320 | PSD broadband, 29 periodogram maxima, no line lattice
  (12000 closed-loop ROM steps per case in 126 s)
In [19]:
plot_psd_comparison(cases)
plt.show()
No description has been provided for this image
In [20]:
plot_delay_portraits(cases)
plt.show()
No description has been provided for this image
In [21]:
def lam_text(lam):
    """A Rosenstein estimate, or the rejection it really is (nan: the fit guards refused)."""
    return f"{lam:.3f}" if lam is not None and np.isfinite(lam) else "rejected"


t0 = time.time()
for c in cases.values():
    # the same equation-free pipeline on both records: Rosenstein lambda_1, D2, classify_regime,
    # signal statistics, and MDS of the full state (the cached FOM field / the recovered ROM field)
    c["ds_fom"] = DataSeries(c["x_fom"], c["dt"], label=f"{c['name']} FOM",
                             Y=c["traj"][:, :N_SIG].T.numpy())
    c["ds_rom"] = DataSeries(c["x_rom"], c["dt"], label=f"{c['name']} qlROM", Y=c["X_rom"].T)
    c["data_fom"], c["data_rom"] = c["ds_fom"].analyze(), c["ds_rom"].analyze()
    c["regime_rom"] = c["data_rom"]["regime"]
    c["lyap_rom"] = dict(text=lam_text(c["data_rom"]["lambda1"]))
print(f"  ({2 * len(cases)} equation-free analyses in {time.time() - t0:.0f} s)")

hdr = (f"{'case':<11} | {'regime FOM':>20} {'FOM record':>20} {'qlROM record':>20} | "
       f"{'lam1 FOM model':>21} {'FOM Rosenstein':>14} {'qlROM Rosenstein':>16} | "
       f"{'d FOM/ROM':>9} {'D2 FOM/ROM':>11} {'T_ph (t.u.)':>21}")
print(hdr)
print("-" * len(hdr))
for c in cases.values():
    print(f"{c['name']:<11} | {c['regime_fom']:>20} {c['data_fom']['regime']:>20} {c['regime_rom']:>20} | "
          f"{c['lyap_fom']['text']:>21} {lam_text(c['data_fom']['lambda1']):>14} "
          f"{c['lyap_rom']['text']:>16} | {c['sig']['dim']:>4d}/{c['sig_rom']['dim']:<4d} "
          f"{c['sig']['D2']:>5.2f}/{c['sig_rom']['D2']:<5.2f} {tph_text(c, 21)}")
print("  regime FOM = classify_regime with the model-based lambda_1 (section 10); FOM / qlROM record =\n"
      "  DataSeries with Rosenstein's lambda_1; '>= X (censored)' = the error never reached 0.5 inside\n"
      f"  the {N_FC}-step window, so X is the window length, not a measured horizon")

  (8 equation-free analyses in 42 s)
case        |           regime FOM           FOM record         qlROM record |        lam1 FOM model FOM Rosenstein qlROM Rosenstein | d FOM/ROM  D2 FOM/ROM           T_ph (t.u.)
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
ks1d qp     | limit_cycle_period_3 limit_cycle_period_3 limit_cycle_period_3 |    transient (9% win)          0.014            0.134 |   10/10    0.50/1.73                   14.7
ks1d chaos  |              chaotic              chaotic     frequency_locked |                 0.067          0.232         rejected |    3/4     2.01/2.51                    3.7
ks2d tw     |        quasiperiodic        quasiperiodic        quasiperiodic |  rejected (no growth)       rejected         rejected |    4/4     2.30/2.34    >= 200.0 (censored)
ks2d chaos  |              chaotic     frequency_locked     frequency_locked |                 0.303       rejected         rejected |    3/4     2.02/2.73                    1.8
  regime FOM = classify_regime with the model-based lambda_1 (section 10); FOM / qlROM record =
  DataSeries with Rosenstein's lambda_1; '>= X (censored)' = the error never reached 0.5 inside
  the 2000-step window, so X is the window length, not a measured horizon

12. The eight-panel diagnostic row, for FOM and qlROM alike¶

DataSeries.characterize(pdf_name=...) draws ntsa's diagnostic row for one record: time series with a zoom inset, PSD, 3-D delay portrait, first return map, Poincare section, recurrence plot, classical-MDS embedding of the full state passed as Y (subsampled to 2000 snapshots), and the Lyapunov panel with the Rosenstein marker -- plus, when the Rosenstein fit is accepted, a second page with its log-divergence fit. It reuses the analysis cached by analyze() above, so nothing is recomputed. Every record goes through the same call, so each FOM page and its qlROM page are directly comparable; they land in figs/03_characterization/, one PDF (and same-name PNG) per record.

In [22]:
t0 = time.time()
plt.close("all")                       # characterize opens its own figures; they live in the PDFs
pdf_dir = _OUT / "03_characterization"
for c in cases.values():
    for side in ("fom", "rom"):
        c[f"ds_{side}"].characterize(pdf_name=str(pdf_dir / f"{c['name'].replace(' ', '_')}_{side}.pdf"))
    plt.close("all")
print(f"characterize(): {2 * len(cases)} records in {time.time() - t0:.0f} s -> {pdf_dir}")
for c in cases.values():
    for rec in (c["data_fom"], c["data_rom"]):
        print(f"  {rec['label']:<17}: regime {rec['regime']:<20} zeta {rec['zeta']:4d}  d {rec['dim']}  "
              f"D2 {rec['D2']:.2f}  std {rec['stats']['std']:.3f}  skew {rec['stats']['skew']:+.3f}  "
              f"kurtosis {rec['stats']['kurtosis']:+.3f}")

Saved figures --> figs/03_characterization/ks1d_qp_fom.pdf

Saved figures --> figs/03_characterization/ks1d_qp_rom.pdf

Saved figures --> figs/03_characterization/ks1d_chaos_fom.pdf

Saved figures --> figs/03_characterization/ks1d_chaos_rom.pdf

Saved figures --> figs/03_characterization/ks2d_tw_fom.pdf

Saved figures --> figs/03_characterization/ks2d_tw_rom.pdf

Saved figures --> figs/03_characterization/ks2d_chaos_fom.pdf

Saved figures --> figs/03_characterization/ks2d_chaos_rom.pdf
characterize(): 8 records in 19 s -> figs/03_characterization
  ks1d qp FOM      : regime limit_cycle_period_3 zeta  119  d 10  D2 0.50  std 3.287  skew +0.016  kurtosis -1.935
  ks1d qp qlROM    : regime limit_cycle_period_3 zeta   28  d 10  D2 1.73  std 3.178  skew -0.319  kurtosis -1.717
  ks1d chaos FOM   : regime chaotic              zeta   77  d 3  D2 2.01  std 1.363  skew -0.153  kurtosis -0.998
  ks1d chaos qlROM : regime frequency_locked     zeta   69  d 4  D2 2.51  std 1.267  skew -0.366  kurtosis -0.125
  ks2d tw FOM      : regime quasiperiodic        zeta   57  d 4  D2 2.30  std 2.747  skew +0.220  kurtosis -0.820
  ks2d tw qlROM    : regime quasiperiodic        zeta   50  d 4  D2 2.34  std 2.731  skew +0.221  kurtosis -0.799
  ks2d chaos FOM   : regime frequency_locked     zeta   85  d 3  D2 2.02  std 2.709  skew -0.018  kurtosis -0.199
  ks2d chaos qlROM : regime frequency_locked     zeta   79  d 4  D2 2.73  std 3.512  skew +0.349  kurtosis +0.613

Reading the table. The two regime columns built from records -- FOM and qlROM, each classified with its own Rosenstein exponent -- are a second opinion, not the computation behind the regime column of section 10, which employs the model-based $\lambda_1$.

  • Where the verdict does not hinge on an exponent, the three columns should agree: tight clusters of maxima outrank $\lambda_1$ in classify_regime, so a clean period-$k$ cycle keeps its label whatever the exponent.
  • Where it does, the columns can split. With no accepted exponent, classify_regime drops through to its two-line ratio test. On a quasi-periodic record this test is the right one; on a broadband record whose two strongest humps sit near a small rational ratio, it reads frequency_locked or quasiperiodic. Before trusting such a label, check the spectrum of the record (sections 10 and 11) and, for the FOM, the tail of its autocorrelation: a broadband spectrum with a decaying autocorrelation marks a chaotic record, and the label is then an artefact of the missing exponent, the same failure that the guard of section 10 exists for.

Rosenstein is a weak witness. Compare each Rosenstein column with the model-based FOM column. A positive estimate on a case whose model-based fit shows no expansion is not chaos, a positive estimate on a chaotic case can still be off by a large factor, and a rejected estimate says nothing either way. The model-based FOM column therefore stays the reference. A model-based exponent for the qlROM -- perturb the reduced state, follow the separation -- needs the qlROM wrapped as a dynamodels model, the adapter kept in the qlROM-DA repository.

13. What to compare, and how to read it¶

Each row of the tables above pairs a FOM record with the qlROM record of the same sensor. Each comparison tests a different property of the attractor.

  • Amplitude. The sensor standard deviation tests whether the qlROM orbits an attractor of the right size. A standard deviation well above the FOM's marks an over-energetic model, which orbits a larger attractor than the attractor it was fitted on; a standard deviation well below marks a model that loses energy and, in the limit, collapses to a fixed point.
  • Clocks. The tallest line, the generators of the line lattice and the spectral centroid (figure (a) of section 11) test whether the qlROM oscillates at the right frequencies. A qlROM can keep the size of its orbit and still lose a generator, or crowd its lines onto new generators; reproducing the number of independent frequencies is a stronger test than landing any single line. A spectral centroid above the FOM's marks high-frequency content that the PDE does not carry, which the discontinuous chart switches of a quantized-local model can deposit.
  • Geometry. The optimal lag $\zeta$, the embedding dimension $d_{\mathrm{FNN}}$ and the correlation dimension $D_2$ (figure (b) of section 11) test whether the reconstructed attractor has the same shape. An extra embedding dimension or a larger $D_2$ means that the reduced dynamics add structure that the PDE does not carry.
  • Regime and exponent. The regime and exponent columns test the kind of dynamics; read them with the caveats of section 12.
  • Trajectory against attractor. Set $T_{\mathrm{ph}}$ against the rest of the row. On a chaotic case, a short $T_{\mathrm{ph}}$ is expected however accurate the model; if the statistics still match, the trajectory is gone but the attractor is of the right size and kind. On a case without exponential separation, a long or censored $T_{\mathrm{ph}}$ should come with statistics that match on every count.

The structural point. Not one line of ntsa in this notebook knows what a qlROM is. The FOMs go through the dynamodels Model protocol (respawn, run_long, leading_lyapunov); the qlROMs go through the equation-free path any measured record takes (free_run, then DataSeries: optimal_lag, false_nearest_neighbours, correlation_dimension, Rosenstein, classify_regime, the diagnostic row). That is why qlroms.utils.characterization, which holds the signal-level helpers shared with tutorial 0, never touches a qlROM: the analysis lives on the other side of a data boundary.

The figure below overlays the delay portraits of the FOM and of the qlROM, both unfolded with the lag of the FOM, so that a difference in geometry comes from the dynamics and not from a different unfolding.

In [23]:
save_figure(plot_attractor_overlay(cases))
plt.show()
No description has been provided for this image
In [24]:
print("SUMMARY")
print(f"  1-D KS [{CASE}], train {NTRAIN}, test {NTEST}: K_opt = {diag['K_opt']}, "
      f"r_check = {diag['r_check']}")
print(f"  (K, r) = ({part1['K']}, {part1['r']}): projection MSE global/local = {part1['ratio']:.3f}, forecast "
      f"RMSE local {part1['rmse_local']:.3e} vs global {part1['rmse_global']:.3e}")
for c in cases.values():
    print(f"  {c['name']:<11} K={c['P']['K']:2d} r={c['P']['r']:2d} | regime FOM {c['regime_fom']:>20} "
          f"-> qlROM record {c['regime_rom']:>20} | lam1 FOM model {c['lyap_fom']['text']}, "
          f"Rosenstein FOM {lam_text(c['data_fom']['lambda1'])} -> qlROM {c['lyap_rom']['text']} | "
          f"d {c['sig']['dim']}->{c['sig_rom']['dim']}, D2 {c['sig']['D2']:.2f}->{c['sig_rom']['D2']:.2f}, "
          f"std {c['sig']['x'].std():.2f}->{c['sig_rom']['x'].std():.2f}, T_ph {tph_text(c)} t.u.")
print(f"  characterization PDFs in {pdf_dir}")
print(f"  runtime {time.time() - t_start:.0f}s")

SUMMARY
  1-D KS [chaotic], train 40000, test 4000: K_opt = 5, r_check = 40
  (K, r) = (5, 40): projection MSE global/local = 1.076, forecast RMSE local 1.664e+00 vs global 1.793e+00
  ks1d qp     K= 6 r=30 | regime FOM limit_cycle_period_3 -> qlROM record limit_cycle_period_3 | lam1 FOM model transient (9% win), Rosenstein FOM 0.014 -> qlROM 0.134 | d 10->10, D2 0.50->1.73, std 3.29->3.18, T_ph 14.7 t.u.
  ks1d chaos  K=10 r=30 | regime FOM              chaotic -> qlROM record     frequency_locked | lam1 FOM model 0.067, Rosenstein FOM 0.232 -> qlROM rejected | d 3->4, D2 2.01->2.51, std 1.36->1.27, T_ph 3.7 t.u.
  ks2d tw     K= 5 r=25 | regime FOM        quasiperiodic -> qlROM record        quasiperiodic | lam1 FOM model rejected (no growth), Rosenstein FOM rejected -> qlROM rejected | d 4->4, D2 2.30->2.34, std 2.75->2.73, T_ph >= 200.0 (censored) t.u.
  ks2d chaos  K=40 r=50 | regime FOM              chaotic -> qlROM record     frequency_locked | lam1 FOM model 0.303, Rosenstein FOM rejected -> qlROM rejected | d 3->4, D2 2.02->2.73, std 2.71->3.51, T_ph 1.8 t.u.
  characterization PDFs in figs/03_characterization
  runtime 328s
In [25]:
close_pdf()

11 figures written to figs/03_diagnosing_a_qlrom.pdf
Out[25]:
PosixPath('figs/03_diagnosing_a_qlrom.pdf')