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 truth long before it stops being a useful model, and a pointwise error cannot tell a good statistical reproduction from a dead one.

The second half runs on ntsa, which treats a ql-ROM exactly as it treats a full-order model: signal-level analysis, attractor geometry, regime classification, Lyapunov exponents, and a characterize() call that assembles the lot. The cases are the ones from tutorial 0.

import os
from pathlib import Path

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

import matplotlib.pyplot as plt

from qlroms.utils.plots import close_pdf, start_pdf
import numpy as np

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.intrusive_qlroms import ks1d
from qlroms.intrusive_qlroms.build_ks import build_fom, get_full_trajectory

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

fom, cfg = build_fom(ks1d, case=CASE, overrides={"Ntrain": NTRAIN, "Ntest": NTEST})
traj = get_full_trajectory(Ntot=cfg["i0"] + NTRAIN + NTEST, model=fom, i0=cfg["i0"])
Xtrain, Xtest = traj[:, :NTRAIN], traj[:, NTRAIN:NTRAIN + NTEST]
save_dir = ks1d.get_simulation_path(model=fom, Ntrain=NTRAIN, Ntest=NTEST, dt=fom.dt)
dt = float(fom.dt)

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

PosixPath('../outputs/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 this measures the charts alone:

$$ \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\| . $$

For $K=1$ this is the ordinary POD projection error. Anything a ql-ROM gains has to show up here first: 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 truth,

$$ \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.

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'})")

projection MSE on the test set (r=40)
  global (K=1) : 1.3666e-01
  local  (K=5)  : 1.4692e-01   ->  ratio global/local = 0.930  (global better)

The ratio came out below 1: at $r=40$ the five charts represent the test data slightly worse than one global basis of the same size. That is not a bug, it is the diagnosis doing its job. With 40 of 128 degrees of freedom the global basis already captures almost everything, and each local basis is fitted on a fifth of the snapshots, so it generalises a little worse off its own cluster. Section 5 shows where in the $(K,r)$ plane the sign flips.

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".

res = compare_local_global_reconstruction(
    Xtest=Xtest, FOM=fom, local_rom=local_rom, global_rom=global_rom, save_dir=save_dir,
)

t = np.arange(NTEST) * dt
fig, axs = plt.subplots(1, 2, figsize=(11, 3.2), sharey=True, layout="constrained")
for ax, (tag, e_loc, e_glo) in zip(
        axs,
        [("representation (a priori)", res["local_rep_err_t"], res["global_rep_err_t"]),
         ("forecast (a posteriori)", res["local_err_t"], res["global_err_t"])], strict=True):
    ax.plot(t, e_glo, color="tab:blue", lw=1.2, label=f"global $r={r}$")
    ax.plot(t, e_loc, color="tab:red", lw=1.2, label=f"local $(r,K)=({r},{K})$")
    ax.set(xlabel="$t$", title=tag, yscale="log")
    ax.grid(alpha=0.25, lw=0.5)
    ax.legend(frameon=False, fontsize=8)
axs[0].set_ylabel(r"$\varepsilon(t)$  [%]")
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}")

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.664058   global RMSE=1.792180  PASS
  representation: local (K=5, r=40) RMSE=0.033879   global RMSE=0.032674 FAIL

representation RMSE:  local 3.3879e-02   global 3.2674e-02
forecast       RMSE:  local 1.6641e+00   global 1.7922e+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. A ql-Galerkin, a qlOpinf and a ql-DMD all satisfy it as they are. A family that does not step reduced coordinates — a reservoir carries its own hidden state — is diagnosed by passing its own closure instead; nothing else changes, because everything measured lives on the charts.

import inspect  # noqa: E402

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

from qlroms.data_driven_qlroms import qlOpinf  # noqa: E402

rom_oi = qlOpinf.from_snapshots(Xtrain, K=K, r=r, dt=dt, random_state=SEED,
                                pod_method="randomized")
rom_oi_g = qlOpinf.from_snapshots(Xtrain, K=1, r=r, dt=dt, random_state=SEED,
                                  pod_method="randomized")
res_oi = compare_local_global_reconstruction(
    Xtest=Xtest, FOM=fom, local_rom=rom_oi, global_rom=rom_oi_g, save_dir=None,
)
print("\nsame diagnosis, operators inferred from data instead of projected:")
print(f"  representation RMSE:  local {res_oi['local_rep_rmse']:.4e}   "
      f"global {res_oi['global_rep_rmse']:.4e}")
print(f"  forecast       RMSE:  local {res_oi['local_rmse']:.4e}   "
      f"global {res_oi['global_rmse']:.4e}")

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

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

Local-vs-global reconstruction sanity check
  forecast:       local (K=5, r=40) RMSE=nan   global RMSE=2.327909  FAIL
  representation: local (K=5, r=40) RMSE=0.031393   global RMSE=0.032559 PASS

same diagnosis, operators inferred from data instead of projected:
  representation RMSE:  local 3.1393e-02   global 3.2559e-02
  forecast       RMSE:  local nan   global 2.3279e+00

The representation numbers land within a few percent of section 2 — qlOpinf.from_snapshots fits its own clustering and POD from the same snapshots, so the charts are equivalent by construction, not bit-identical — while the forecast numbers are completely different (here the local OpInf run blows up to nan on this chaotic case at $r=40$, which is a result about the operators, not about the charts). That separation, geometry on one side and dynamics on the other, is the whole point of running both halves of the diagnosis.

scripts/ks/families.py wraps this in one function, build_family(name, ...) -> (build_model_fn, forecast_fn), with galerkin, opinf and esn (per-chart echo state networks, seeded from their prior maps — the family that proves the hook is really doing the work). The scripts take --family and pass nothing else along.

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. BIC keeps 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 this notebook is instant once runs/sweep.sh has been through.

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=4133361.76 (cached)
  BIC: K=5  BIC=3987461.18 (cached)
  BIC: K=10  BIC=3836240.85 (cached)
  BIC: K=20  BIC=3671477.92 (cached)


Sweeping global (K=1) ROM...

  K=1, r= 40  err=1.7922  proj_mse_te=1.3666e-01  proj_mse_tr=1.1161e-01  kld=6.3196e+00

  K=1, r= 30  err=2.3963  proj_mse_te=2.8140e+00  proj_mse_tr=2.0137e+00  kld=4.9448e-01

  K=1, r= 20  err=15.3767  proj_mse_te=2.2019e+01  proj_mse_tr=1.8259e+01  kld=2.0361e-01

  K=1, r= 10  err=nan  proj_mse_te=6.3967e+01  proj_mse_tr=5.5756e+01  kld=inf
Sweeping local ROM K=2...

  K=2, r= 40  err=1.6684  proj_mse_te=1.3043e-01  proj_mse_tr=1.1588e-01  kld=2.5003e+00

  K=2, r= 30  err=2.2223  proj_mse_te=2.6952e+00  proj_mse_tr=1.9823e+00  kld=5.4137e-01

  K=2, r= 20  err=3.4399  proj_mse_te=1.9655e+01  proj_mse_tr=1.6997e+01  kld=7.7490e-01

  K=2, r= 10  err=7.6036  proj_mse_te=6.2545e+01  proj_mse_tr=5.4924e+01  kld=1.7464e+01
Sweeping local ROM K=5...

  K=5, r= 40  err=1.6641  proj_mse_te=1.4692e-01  proj_mse_tr=1.2000e-01  kld=1.0653e+01

  K=5, r= 30  err=1.9938  proj_mse_te=2.3514e+00  proj_mse_tr=2.5109e+00  kld=4.2963e+00

  K=5, r= 20  err=2.4008  proj_mse_te=2.0048e+01  proj_mse_tr=1.6217e+01  kld=9.5918e-01

  K=5, r= 10  err=2.0471  proj_mse_te=5.7335e+01  proj_mse_tr=5.2667e+01  kld=1.3074e+01
Sweeping local ROM K=10...

  K=10, r= 40  err=1.7633  proj_mse_te=1.9034e-01  proj_mse_tr=1.9447e-01  kld=1.9349e+00

  K=10, r= 30  err=1.9043  proj_mse_te=3.4288e+00  proj_mse_tr=2.7397e+00  kld=4.4494e-01

  K=10, r= 20  err=1.8814  proj_mse_te=2.0059e+01  proj_mse_tr=1.6426e+01  kld=3.6250e+00

  K=10, r= 10  err=2.0567  proj_mse_te=6.6800e+01  proj_mse_tr=5.5256e+01  kld=7.0199e+00
Sweeping local ROM K=20...

  K=20, r= 40  err=1.8260  proj_mse_te=5.5806e-01  proj_mse_tr=5.2871e-01  kld=1.3045e+01

  K=20, r= 30  err=1.8017  proj_mse_te=5.9067e+00  proj_mse_tr=4.7753e+00  kld=5.6893e+00

  K=20, r= 20  err=1.9002  proj_mse_te=2.3460e+01  proj_mse_tr=2.0148e+01  kld=5.1060e+00

  K=20, r= 10  err=1.9302  proj_mse_te=6.4922e+01  proj_mse_tr=5.9272e+01  kld=9.3152e+00
  Loaded cluster probs/transition matrices from cache (K_opt=5).
  KL(P_test  || P_train)      = 1.87e-01  (cached)
  KL(P_test  || P_qlrom/test) = 6.62e-01  (cached)
  KL(P_test  || P_grom /test) = 6.53e-01  (cached)
  KL(P_test  || P_qlrom/repr) = 2.42e-07  (cached)
  KL(P_test  || P_grom /repr) = 0.00e+00  (cached)

K_opt = 5, r_check = 40

fig, _ = plot_kr_error_map(diag)
plt.show()

Read the maps left to right. Projection MSE falls with $r$ and barely moves with $K$: on this case the modes do the work, the charts fine-tune. The ratio global/local is the actual verdict — green is a ql-ROM win, red a loss, and the $K=1$ row is 1 by construction. Note where it turns red: many charts at large $r$. Each chart then fits its own cluster with almost as many modes as the cluster has structure, and the basis stops generalising — the train/test gap in the next figure is the same effect seen from the side. The KLD map scores the closed-loop run's cluster occupancy against the data; it is the only panel that looks at where the trajectory goes rather than how far it is at each instant.

fig, _ = plot_parameter_selection(diag)
plt.show()

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 is the overfitting the ratio map hinted at: more charts always improve the training fit, and past a point stop improving the test fit.

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.

fig, _ = plot_cluster_occupancy(diag)
plt.show()
fig, _ = plot_transition_matrix(diag)
plt.show()

The left occupancy panel is the a priori one: projecting and recovering the true snapshots cannot move them out of their charts, so both ROMs reproduce the data histogram to numerical precision ($D_{KL}\sim 10^{-7}$). The right panel is the closed-loop run, and it is a different story — a long free run on a chaotic attractor redistributes itself, and the KL divergence is the honest measure of how badly. The transition matrices are strongly diagonal in every panel simply because the dwell time is long compared with $\Delta t$: a chart is left after many steps, not at the next one.

7. Error time series across $K$

fig, _ = plot_error_timeseries(diag, dt=dt)
plt.show()

8. The physical page

The scalars above are summaries; the last figure is 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 run (below) and the forecast run. This is the page runs/diagnosis.sh writes to a PDF for a single $(K,r)$.

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()

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()

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, --Ntrain/--Ntest, the clustering flags, and --family {galerkin,opinf,esn}. runs/sweep.sh runs under nohup and tails its log, since 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 panel 1-2
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
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 family enters once, as forecast_fn, which is why the same diagnosis serves tutorials 1 (intrusive), 2 (qlOpinf) and 3 (qlESNs) without a line of family-specific code.

10. Does it keep the character of the system?

Everything above measures how far the model is from the truth. Past a chaotic system's predictability horizon that distance saturates and stops discriminating: a model that has drifted onto the attractor and one that has collapsed to a fixed point can score the same.

Tutorial 0 characterized the four 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 cell below rebuilds the case baselines silently; nothing about the method is repeated.

import contextlib
import io
import os
import time
from itertools import product
from pathlib import Path

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

import matplotlib.pyplot as plt
import numpy as np
import torch
from dynamodels.physical import KS, KS2D
from ntsa import lyapunov as lyap
from ntsa import tools as nt
from ntsa.characterize import characterize
from scipy.signal import find_peaks

from qlroms import free_run, metrics
from qlroms.intrusive_qlroms import ks1d, ks2d
from qlroms.intrusive_qlroms.build_ks import build_fom, get_full_trajectory
from qlroms.model import QLModel
from qlroms.utils.config import get_case
from qlroms.utils.plots import field_image, outputs_dir, save_figure

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

t_start = time.time()
C_FOM, C_ROM, C_MARK = "#2a78d6", "#eb6834", "#1baf7a"
CASE_COLORS = ["#2a78d6", "#1baf7a", "#b05fd6", "#d64a4a"]
SEED = 0

N_SIG = 12_000        # samples fed to every signal-level diagnostic
N_REC = 1_200         # samples in a recurrence window (recurrence_matrix is O(T^2))
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
# qlroms.utils.config -- t_lam / t_lam_rom are the Lyapunov horizons, which must cover the
# growth phase and stop near saturation (the ks2d chaos ROM saturates sooner than its FOM,
# so it gets a shorter one). Only the WINDOWS are this notebook's own: the cases are
# generated far longer than a ~15 min notebook can characterize, and the text says so
# wherever a number depends on it.
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)}
CASES = {name: {**get_case(name), **window} for name, window in WINDOWS.items()}

_OUT = outputs_dir()          # notebooks/outputs, from the repo root or from tutorials/

def load_case(name):
    """Cached FOM trajectory plus the metadata every later section needs."""
    P = CASES[name]
    mod = P["module"]
    is2d = mod is ks2d
    fom, cfg = build_fom(mod, case=P["case"], overrides={"Ntrain": P["Ntrain"], "Ntest": P["Ntest"]})
    traj = get_full_trajectory(Ntot=cfg["i0"] + P["Ntrain"] + P["Ntest"], model=fom, i0=cfg["i0"])
    dt = float(fom.dt)
    obs_idx = np.asarray(fom.sensor_locations, dtype=int)
    # every diagnostic below reads ONE scalar sensor: the same physical point for FOM and ROM
    i_obs = 1
    # ntsa re-integrates the FOM (respawn + run_long) instead of replaying the cache, so seed it
    # from the LAST cached state: the KS2D state IS the physical field, the 1-D KS state is its rfft
    u_end = np.asarray(traj[:, -1], dtype=float)
    fom_run = nt.respawn(fom, psi0=(u_end if is2d else np.fft.rfft(u_end)))
    fom_run.t_transient, fom_run.t_lyap = P["t_transient"], P["t_lyap"]
    return dict(name=name, P=P, fom=fom, fom_run=fom_run, cfg=cfg, dt=dt, obs_idx=obs_idx, i_obs=i_obs,
                traj=traj, Xtrain=traj[:, :P["Ntrain"]], Xtest=traj[:, P["Ntrain"]:],
                x_fom=np.asarray(traj[obs_idx[i_obs], :N_SIG], dtype=float),
                is2d=is2d, lamb1_ref=float(fom.lamb1))


cases = {name: load_case(name) for name in CASES}
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

# baselines: the tutorial-0 measurements for the FOM side, recomputed here
# so this notebook stands alone -- see tutorial 0 for what they mean

ACF_TAIL_DISCRETE = 0.3   # |rho| still this large at the end of the lag window => discrete spectrum


def spectral_lines(psd_f, psd, prominence=0.05, n_max=8, m_max=16, gen_max=4):
    """Prominent PSD lines, and the fewest generators whose integer combinations explain them.

    `find_peaks` sees LINES where `argsort` sees BINS: the two strongest bins of one sharp
    line are that line and its leakage neighbour. A combination of order sum(m) inherits
    half a bin of quantization from each multiple, hence the order-dependent tolerance.
    Multiples above `m_max` do not count as explanations, and a set of lines needing
    `gen_max` generators is not a lattice -- both caps also keep the search finite.
    """
    df = float(psd_f[1] - psd_f[0])
    idx, props = find_peaks(psd, prominence=prominence * psd.max())
    strongest = idx[np.argsort(props["prominences"])[::-1][:n_max]]
    lines, gens = np.sort(psd_f[strongest]), []
    for fk in lines:
        combos = product(*(range(min(int(fk / g) + 2, m_max)) for g in gens))
        if any(abs(fk - float(np.dot(ms, gens))) <= 0.5 * df * (1 + sum(ms))
               for ms in combos if any(ms)):
            continue
        gens.append(float(fk))
        if len(gens) >= gen_max:
            break
    return lines, gens, len(idx)


def signal_report(x, dt, tag):
    """Every scalar diagnostic ntsa computes from one signal, in one dict."""
    out = dict(tag=tag, x=x, dt=dt)
    f, psd = nt.fun_PSD(dt, x[None, :])
    out["f"], out["psd"] = f[1:], psd[0][1:]                      # bin 0 is the mean, always dropped
    out["f_top"] = float(out["f"][np.argmax(out["psd"])])         # tallest bin: a line only on a torus
    out["lines"], out["gen"], out["n_lines"] = spectral_lines(out["f"], out["psd"])
    # the centroid is the broadband statistic that stays meaningful on a continuum
    out["f_c"] = float(np.sum(out["f"] * out["psd"]) / np.sum(out["psd"]))
    n_lags = min(len(x) // 4, 4000)
    out["acf"] = nt.autocorrelation(x, n_lags)
    zc = np.flatnonzero(out["acf"] < 0)
    out["t_dec"] = float(zc[0] * dt) if zc.size else float(n_lags * dt)
    out["acf_tail"] = float(np.abs(out["acf"][3 * len(out["acf"]) // 4:]).max())
    out["discrete"] = bool(out["acf_tail"] > ACF_TAIL_DISCRETE)   # gates the lattice reading
    max_lag = min(len(x) // 10, 800)
    out["ami"] = nt.average_mutual_information(x, max_lag)
    out["zeta"] = int(nt.optimal_lag(x, max_lag=max_lag))
    out["dim"], out["fnn"] = nt.false_nearest_neighbours(x, out["zeta"])
    out["Y"] = nt.delay_embed(x, max(out["dim"], 3), out["zeta"])
    out["D2"], out["cd"] = nt.correlation_dimension(out["Y"])
    # D2 is fitted over the 2nd-50th percentile of the pair distances: a LARGE-scale slope.
    # Split it to see whether the object keeps that dimension as r shrinks.
    log_r, log_C = out["cd"]
    out["D2_lo"] = float(np.polyfit(log_r[:8], log_C[:8], 1)[0])    # small-r end of that window
    out["D2_hi"] = float(np.polyfit(log_r[-8:], log_C[-8:], 1)[0])  # large-r end
    return out


def spectrum_text(s):
    """PSD verdict: a line lattice only where the autocorrelation says there are lines."""
    if not s["discrete"]:
        return f"broadband, {s['n_lines']:2d} periodogram maxima, no line lattice"
    return (f"{s['n_lines']:2d} prominent lines, the strongest {len(s['lines'])} on "
            f"{len(s['gen'])} generator(s) {np.round(s['gen'], 4)}")


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 four 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.3299 | ACF first zero 6.6 t.u., tail |rho| 0.81 | optimal lag zeta = 50 samples = 5.00 t.u. | FNN embedding dim = 4 | D2 = 2.28 (small-r 2.11, large-r 2.37)

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 four signal reports in 3.5 s)

def lyapunov_verdict(lam1, res):
    """(lam1 to trust or None, printable string) -- see the non-normality caveat above."""
    t = res["t"]
    frac = float((t[res["i2"]] - t[res["i1"]]) / (t[-1] - t[0]))
    if not np.isfinite(lam1):
        return None, "rejected (no growth)", frac
    if frac < EXP_WINDOW_MIN:
        return None, f"transient ({100 * frac:.0f}% win)", frac
    return float(lam1), f"{lam1:.3f}", frac


def saturation_rate(res):
    """(nats climbed, time taken, mean rate) from the start of the record to saturation.

    A record that spends `span` time units climbing `climb` nats is described by an
    exponential of rate climb/span. A claimed exponent much LARGER than that is
    contradicted by the record itself: it would have exhausted the same climb sooner.
    """
    m, t = res["mean_log_sep"], res["t"]
    hit = np.flatnonzero(m >= np.log(res["sat"]))
    i_sat = int(hit[0]) if hit.size else len(m) - 1
    climb, span = float(m[i_sat] - m[0]), float(t[i_sat] - t[0])
    return climb, span, (climb / span if span > 0 else np.nan)


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)
    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"  (four 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
  (four leading_lyapunov runs in 62 s)

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.28      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 1 and wrap it in QLModel:

QLModel(rom, obs_idx, psi0)   # obs_idx = PHYSICAL sensor rows, psi0 = a physical state

QLModel presents the fitted qlROM as a dynamodels Model: Nphi global-atlas coordinates as the state, Nq sensors as observables, time_integrate, update_history, get_observable_hist, and fixed_params = ['rom','obs_idx', 'est_theta'] so ntsa.respawn can rebuild it (an atlas-frame psi0 is accepted directly, which is what respawn hands back). That is the entire integration effort: no ntsa code below is qlROM-aware.

First, the contrast that motivates the whole notebook. A closed-loop free run is scored by its prediction horizon $T_{\mathrm{ph}}$, the time at which the relative error first exceeds 0.5 -- a trajectory measure. Then the ROM is integrated for 12000 steps with no reference at all, and the section-2/3 diagnostics are recomputed on its sensor signal: an attractor measure.

One reporting rule, because it applies to two of the four cases. 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 >= 200.0 (censored): the run was still tracking when we stopped looking, and the only honest statement is a lower bound.

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"]
    save_dir = mod.get_simulation_path(model=c["fom"], Ntrain=P["Ntrain"], Ntest=P["Ntest"], dt=c["dt"])
    tb = time.time()
    c["rom"] = mod.build_local_model(c["Xtrain"].clone(), c["fom"], r=P["r"], K=P["K"],
                                     save_dir=save_dir, clustering_kwargs={"random_state": SEED})
    t_build = time.time() - tb
    # trajectory-wise score, for contrast only (tutorial 1's truncation rule)
    X_rom, _ = free_run(c["rom"], c["Xtest"][:, 0], n_steps=N_FC)
    X_rom = X_rom.cpu()
    bad = (~torch.isfinite(X_rom).all(dim=0)) | (X_rom.norm(dim=0) > 2 * c["Xtest"][:, :N_FC].norm(dim=0))
    n_valid = int(bad.to(torch.int).argmax()) if bad.any() else N_FC
    err = metrics.relative_error_series(c["Xtest"][:, :n_valid], X_rom[:, :n_valid])
    c["T_ph"] = metrics.prediction_horizon(err, c["dt"], threshold=0.5)
    # 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(err) and err.max() <= 0.5)
    c["qm"] = QLModel(c["rom"], c["obs_idx"], c["Xtest"][:, 0].numpy(), t_CR=float(c["fom"].t_CR))
    c["qm"].t_lyap = P["t_lyap"]
    print(f"{c['name']:11s}: qlGalerkin K={P['K']:2d} r={P['r']:2d} built in {t_build:4.1f} s | "
          f"QLModel Nphi = {c['qm'].Nphi:4d} atlas coords, Nq = {c['qm'].Nq} | free run valid "
          f"{n_valid}/{N_FC} steps, T_ph(0.5) = {tph_text(c)} t.u.")
print(f"  (builds + free runs in {time.time() - t0:.0f} s)")

ks1d qp    : qlGalerkin K= 6 r=30 built in  0.0 s | QLModel Nphi =  102 atlas coords, Nq = 4 | free run valid 2000/2000 steps, T_ph(0.5) = 14.6 t.u.

ks1d chaos : qlGalerkin K=10 r=30 built in  0.0 s | QLModel Nphi =   89 atlas coords, Nq = 4 | free run valid 2000/2000 steps, T_ph(0.5) = 3.7 t.u.

ks2d tw    : qlGalerkin K= 5 r=25 built in  0.1 s | QLModel Nphi =  255 atlas coords, Nq = 4 | free run valid 2000/2000 steps, T_ph(0.5) = >= 200.0 (censored) t.u.

ks2d chaos : qlGalerkin K=40 r=50 built in  1.2 s | QLModel Nphi = 1626 atlas coords, Nq = 4 | free run valid 2000/2000 steps, T_ph(0.5) = 1.8 t.u.
  (builds + free runs in 14 s)

t0 = time.time()
for c in cases.values():
    # run_long MUTATES and closes the model it is given, so always respawn first
    _, y, _ = nt.run_long(nt.respawn(c["qm"]), t_run=N_SIG * c["dt"], t_transient=0.0,
                          trim_transient=False)
    c["x_rom"] = np.asarray(y[:N_SIG, c["i_obs"]], dtype=float)
    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.157 ( -4.0%) | zeta 119 ->  28 | d_FNN 10 -> 10 | D2 0.50 -> 1.82 | tallest f 0.0208 -> 0.0658 | f_centroid 0.3295 -> 0.3990 | PSD 14 prominent lines, the strongest 8 on 4 generator(s) [0.0375 0.0583 0.0658 0.0692]

ks1d chaos : FOM std 1.363 -> ROM 1.282 ( -5.9%) | zeta  77 ->  72 | d_FNN 3 -> 4 | D2 2.01 -> 2.68 | tallest f 0.0067 -> 0.0117 | f_centroid 0.3996 -> 0.7603 | PSD broadband, 34 periodogram maxima, no line lattice

ks2d tw    : FOM std 2.749 -> ROM 2.754 ( +0.2%) | zeta  50 ->  50 | d_FNN 4 -> 4 | D2 2.28 -> 2.35 | tallest f 0.0300 -> 0.0475 | f_centroid 0.3299 -> 0.0973 | PSD  9 prominent lines, the strongest 8 on 2 generator(s) [0.0175 0.03  ]

ks2d chaos : FOM std 2.709 -> ROM 3.526 (+30.2%) | zeta  85 ->  80 | d_FNN 3 -> 4 | D2 2.02 -> 2.39 | tallest f 0.0750 -> 0.1667 | f_centroid 3.0566 -> 5.2694 | PSD broadband, 33 periodogram maxima, no line lattice
  (12000 closed-loop ROM steps per case in 53 s)

fig, axs = plt.subplots(1, 4, figsize=(13, 3.0), layout="constrained")
for ax, c in zip(axs, cases.values(), strict=True):
    s, sr = c["sig"], c["sig_rom"]
    ax.semilogy(s["f"], s["psd"], lw=0.9, color=C_FOM, label="FOM")
    ax.semilogy(sr["f"], sr["psd"], lw=0.9, color=C_ROM, alpha=0.8, label="qlROM")
    ax.set(xlim=(0, min(8 * s["f_top"], s["f"][-1])), xlabel="$f$", title=c["name"])
    ax.legend(frameon=False, fontsize=8)
    ax.grid(alpha=0.25, lw=0.5)
axs[0].set_ylabel("PSD")
fig.suptitle("(a) Power spectrum, full-order vs quantized-local ROM", fontsize=11)
plt.show()

fig, axs = plt.subplots(2, 4, figsize=(13, 6.0), layout="constrained")
for j, c in enumerate(cases.values()):
    for i, (s, col, tag) in enumerate([(c["sig"], C_FOM, "FOM"), (c["sig_rom"], C_ROM, "qlROM")]):
        Y = s["Y"][::2]
        axs[i, j].plot(Y[:, 0], Y[:, 1], lw=0.6, color=col, alpha=0.8)
        axs[i, j].set(xlabel="$x(t)$", ylabel=r"$x(t+\zeta)$",
                      title=f"{c['name']} {tag}: $\\zeta$={s['zeta']}, $D_2$={s['D2']:.2f}")
        axs[i, j].grid(alpha=0.25, lw=0.5)
fig.suptitle("(b) Delay-embedding portraits: same axes, same lag convention, FOM (top) vs qlROM (bottom)",
             fontsize=11)
plt.show()

t0 = time.time()
for c in cases.values():
    try:
        lam1, lam1_std, res = lyap.leading_lyapunov(c["qm"], n_pert=3, t_run=c["P"]["t_lam_rom"])
        lam_ok, lam_str, frac = lyapunov_verdict(lam1, res)
    except RuntimeError as exc:                     # a ROM that leaves its charts cannot be measured
        lam1, lam1_std, res, lam_ok, lam_str, frac = np.nan, 0.0, None, None, f"failed: {exc}", 0.0
    c["lyap_rom"] = dict(lam1=lam1, std=lam1_std, res=res, ok=lam_ok, text=lam_str, frac=frac)
    label, ev = nt.classify_regime(c["x_rom"], c["dt"], lam1=lam_ok,
                                   lam1_std=(lam1_std if lam_ok is not None else 0.0),
                                   t_total=N_SIG * c["dt"])
    c["regime_rom"], c["evidence_rom"] = label, ev
print(f"  (four ROM leading_lyapunov runs in {time.time() - t0:.0f} s)")

hdr = (f"{'case':<12} | {'regime FOM':>16} {'regime ROM':>18} | {'lam1 FOM':>22} {'lam1 ROM':>22} | "
       f"{'d FOM/ROM':>10} {'D2 FOM/ROM':>12} {'T_ph (t.u.)':>22}")
print(hdr)
print("-" * len(hdr))
for c in cases.values():
    print(f"{c['name']:<12} | {c['regime_fom']:>16} {c['regime_rom']:>18} | "
          f"{c['lyap_fom']['text']:>22} {c['lyap_rom']['text']:>22} | "
          f"{c['sig']['dim']:>4d}/{c['sig_rom']['dim']:<5d} "
          f"{c['sig']['D2']:>5.2f}/{c['sig_rom']['D2']:<6.2f} {tph_text(c, 22)}")
print(f"  ('>= X (censored)' = the relative error never reached 0.5 inside the {N_FC}-step window,\n"
      f"   so X is the window length, not a measured horizon)")

  (four ROM leading_lyapunov runs in 63 s)
case         |       regime FOM         regime ROM |               lam1 FOM               lam1 ROM |  d FOM/ROM   D2 FOM/ROM            T_ph (t.u.)
---------------------------------------------------------------------------------------------------------------------------------------------------
ks1d qp      | limit_cycle_period_3 limit_cycle_period_3 |     transient (9% win)   rejected (no growth) |   10/10     0.50/1.82                     14.6
ks1d chaos   |          chaotic            chaotic |                  0.067                  0.151 |    3/4      2.01/2.68                      3.7
ks2d tw      |    quasiperiodic      quasiperiodic |   rejected (no growth)   rejected (no growth) |    4/4      2.28/2.35      >= 200.0 (censored)
ks2d chaos   |          chaotic            chaotic |                  0.303                  1.383 |    3/4      2.02/2.39                      1.8
  ('>= X (censored)' = the relative error never reached 0.5 inside the 2000-step window,
   so X is the window length, not a measured horizon)

12. One characterize() call, one diagnostic PDF

ntsa.characterize(models, ...) runs the whole of sections 2-4 for each model and draws an eight-panel row per case: 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 (not the sensor), and a Lyapunov panel. Eight models -- four FOMs and their four ROMs -- go into a single call and come out as one PDF, whose structure is fixed by the arguments: one row-grid page per rows_per_page=4 models, then one page per model for each extra panel requested. With mds=True and both Lyapunov panels off, that is two row grids followed by the eight MDS pages. The cell prints the page count it actually wrote, rather than this text asserting one.

Three settings are worth knowing. lyapunov=False, because sections 4 and 6 already measured $\lambda_1$ with per-case windows and repeating it here would roughly double the notebook's runtime. t_transient=None (the default) so each model's own t_transient is used -- short for every model here, because all eight start from a state already on the attractor (the FOMs from the cached seed, QLModel from the first test snapshot, hence its t_transient = 0). And spectrum='auto', which silently skips the full Lyapunov spectrum for all eight: they are DiscreteIntegrator maps with no time_derivative, and ntsa declines rather than guessing.

One warning is expected and harmless: mds=True feeds the raw state to classical_mds, and the 1-D KS state is the spectral coefficient vector, so numpy reports a ComplexWarning when it takes the real part. Pass mds=False if you would rather not see it, or run classical_mds yourself on the physical field.

t0 = time.time()
plt.close("all")                       # characterize opens a figure per page; start from a clean slate
models = [m for c in cases.values() for m in (c["fom_run"], c["qm"])]
labels = [f"{c['name']} {tag}" for c in cases.values() for tag in ("FOM", "qlROM")]
pdf_path = str(_OUT / "03_characterization_panels.pdf")
results = characterize(models, labels=labels, obs_idx=1, t_run=200.0, lyapunov=False,
                       mds=True, pdf_name=pdf_path)
# page count from what characterize actually assembled: row grids of `rows_per_page` models,
# then one page per model for each optional panel it computed
n_grid = -(-len(results) // 4)                      # rows_per_page = 4 (the default)
n_extra = {lab: sum(r[k] is not None for r in results)
           for k, lab in {"lyap_fit": "Lyapunov-fit", "spectrum": "spectrum", "gamma": "MDS"}.items()}
n_pages = n_grid + sum(n_extra.values())
print(f"characterize(): {len(results)} rows in {time.time() - t0:.0f} s -> {pdf_path}\n"
      f"  {n_pages} pages = {n_grid} row grid(s) + "
      + " + ".join(f"{v} {lab}" for lab, v in n_extra.items() if v))
for r in results:
    print(f"  {r['label']:<18}: regime {r['regime']:<18} zeta {r['zeta']:4d}  d {r['dim']}  "
          f"std {r['stats']['std']:.3f}  skew {r['stats']['skew']:+.3f}  "
          f"kurtosis {r['stats']['kurtosis']:+.3f}")

-- characterizing: ks1d qp FOM

/storage0/anovoama/ntsa/ntsa/tools.py:338: ComplexWarning: Casting complex values to real discards the imaginary part
  X = np.asarray(X, dtype=float)

-- characterizing: ks1d qp qlROM

-- characterizing: ks1d chaos FOM

-- characterizing: ks1d chaos qlROM

-- characterizing: ks2d tw FOM

-- characterizing: ks2d tw qlROM

-- characterizing: ks2d chaos FOM

-- characterizing: ks2d chaos qlROM

Saved figures --> ../outputs/03_characterization_panels.pdf
characterize(): 8 rows in 137 s -> ../outputs/03_characterization_panels.pdf
  10 pages = 2 row grid(s) + 8 MDS
  ks1d qp FOM       : regime limit_cycle_period_3 zeta  122  d 10  std 3.289  skew +0.040  kurtosis -1.933
  ks1d qp qlROM     : regime limit_cycle_period_3 zeta   28  d 10  std 3.184  skew -0.305  kurtosis -1.734
  ks1d chaos FOM    : regime frequency_locked   zeta   43  d 3  std 1.127  skew -0.215  kurtosis -0.344
  ks1d chaos qlROM  : regime frequency_locked   zeta   39  d 4  std 1.108  skew -0.623  kurtosis +0.306
  ks2d tw FOM       : regime limit_cycle_period_7 zeta   21  d 4  std 2.809  skew +0.146  kurtosis -0.975
  ks2d tw qlROM     : regime limit_cycle_period_7 zeta   15  d 4  std 2.771  skew +0.142  kurtosis -0.906
  ks2d chaos FOM    : regime quasiperiodic      zeta  129  d 3  std 3.113  skew +0.116  kurtosis -0.772
  ks2d chaos qlROM  : regime quasiperiodic      zeta   86  d 4  std 3.652  skew +0.043  kurtosis -0.403

These labels are a second opinion, not the same computation, and not one of the eight matches the label section 6 gave the same model. characterize made its own 200 t.u. record, against the 120-1200 t.u. windows sections 2-4 used -- shorter for three of the four cases but longer for ks2d chaos, so record length alone does not explain every disagreement -- and, with lyapunov=False, classified it with no $\lambda_1$ at all -- so classify_regime falls back on peak clustering and the PSD ratio test, and a closed orbit sampled for four or five periods reads as limit_cycle_period_3.

One disagreement deserves to be named rather than covered by that general remark: ks1d chaos comes back frequency_locked, for both its FOM and its ROM -- the case this notebook calls chaotic on a measured $\lambda_1=0.066$ over a 250 t.u. run. The mechanism is printed in the section-4 table: with no exponent to reach step (c), the classifier drops through to its two-line ratio test, and this signal's $f_2/f_1$ is $1.500$, an exact 3/2. A ratio test cannot tell a locked torus from a chaotic attractor whose two strongest spectral humps happen to sit 3:2 apart, and the ACF tail of 0.13 measured in section 2 says which of the two this is. Deny classify_regime its dynamical evidence and it will answer from the spectrum alone; the answer is only as good as that.

The same row carries a second warning worth reading rather than scrolling past: the log above it shows run_long trimming 107 t.u. -- 53% of the record -- as residual transient drift from the ks1d chaos qlROM, leaving about 93 t.u. to classify, and the $\zeta=281$ samples (14 t.u.) that row reports against the FOM's 2.8 t.u. is what a delay lag estimated on such a record looks like. A ROM re-entering its charts from a projected initial condition needs its own settling time, and characterize's uniform t_transient does not know that. The FOM and ROM rows of each pair do land on the same or an adjacent label under identical settings, which is the useful part of this table. The authoritative comparison is the one in section 6, where each model was given its own measured exponent.

13. What matched, and what did not

The comparison table in section 6 is the result. Reading it honestly:

What matched.

• The two non-chaotic cases keep their shape, and -- the sharpest version of that claim -- their line lattices. ks1d qp's ROM puts its prominent lines on the same single generator 0.0225 as the FOM; ks2d tw's ROM needs the same two, 0.0175 and 0.0300. Around that: on ks1d qp the ROM's tallest line lands on the FOM's to four digits (0.0450), the delay lag agrees to one sample (57 vs 58), the embedding dimension is identical (4), $D_2$ agrees to 0.03 and both models have their Lyapunov fit rejected. For ks2d tw the sensor standard deviation matches to 0.2% (2.753 vs 2.748), $D_2$ to 0.08, and both are classified quasiperiodic. Reproducing the number of independent frequencies is a stronger claim than reproducing any one of them: a ROM can land the loudest line and still be turning on the wrong number of clocks. ks1d qp is the point of the exercise -- a 6-chart, 30-mode model with 90 atlas coordinates reproducing the closed orbit of a 128-dimensional PDE -- and it is the one case where the amplitude does not come along: its sensor standard deviation drops from 2.265 to 2.068, $-8.7\%$, the largest amplitude error of the three ROMs that do not inflate. Same object, same frequencies, same dimension, 8.7% too small. The travelling wave, whose free run tracks the truth over the whole window, has no such gap.
• On ks1d chaos the free run is useless after $T_{\mathrm{ph}}\approx4$ t.u., and the ROM's sensor standard deviation still matches the FOM's to 0.4% (1.292 vs 1.287) over 12000 closed-loop steps, with the PSD lying on top of the FOM's over the resolved band. That is the whole argument of the tutorial in one line: the trajectory is gone and the attractor is not.

What did not.

• Every ROM comes out one embedding dimension higher than its FOM on the chaotic cases ($3\to4$), with $D_2$ up by 0.5 (1.76 to 2.27, and 1.98 to 2.52). The reduced dynamics are adding structure the PDE does not have, and the same finger points at the spectra: the ROM's spectral centroid is higher than the FOM's on ks1d chaos and an order of magnitude higher on ks2d chaos. Chart switching in a quantized-local model is a discontinuous event in the reduced coordinates, and it deposits high-frequency content the continuous field does not carry.
• ks2d chaos is over-energetic on top of that. The $K=40$ ROM stays finite for the whole 120 t.u. run, but its sensor standard deviation is about 20% above the FOM's (3.53 vs 2.96): it orbits a larger attractor than the one it was fitted on. At $K=10$ (tutorial 1's setting) it does not get that far at all -- the free run goes non-finite after roughly 8 time units and QLModel raises rather than returning a garbage record, which is why this tutorial uses $K=40$.
• The Lyapunov exponents are the weakest column in the table. On ks1d chaos the ROM reads about 70% high (0.11 against 0.07); on ks2d chaos it reads roughly three times the FOM's (0.59 against 0.19) -- consistent with the inflated variance and $D_2$, i.e. the same defect seen a third way, and both disagree with the recorded 1.8. Nothing was tuned to close those gaps. Converging them needs longer runs, more perturbations and several base points on the attractor; within this budget the honest claim is that the exponents agree in sign and order of magnitude, no more.
• classify_regime is only as good as the evidence it is given. On ks1d qp, where both fits are rejected and no $\lambda_1$ reaches it, the label falls back to the PSD rational-ratio test and reads frequency_locked for the ROM against quasiperiodic for the FOM. The two spectra agree; the two frequency ratios land on opposite sides of a rationality threshold. Note which side is the artefact: the section-4 table prints the FOM's $f_2/f_1$ as 2.519, which is 5/2 missed by 0.8%, and section 2 put that signal's eight strongest lines on a single generator. The FOM's quasiperiodic is the label produced by bin quantization; a locked reading is the defensible one for both models, and the ROM/FOM difference here is a threshold effect on one spectrum, not a difference between two attractors.

The structural point. Not one line of ntsa in this notebook knows what a qlROM is. QLModel makes a fitted quantized-local model satisfy the dynamodels Model protocol, and respawn, run_long, optimal_lag, false_nearest_neighbours, correlation_dimension, classify_regime, leading_lyapunov and characterize then apply to it exactly as they apply to the PDE solver. That is why qlroms ships no characterization module: the boundary is the protocol, and the analysis lives on the other side of it.

fig, axs = plt.subplots(1, 4, figsize=(13, 3.2), layout="constrained")
for ax, c in zip(axs, cases.values(), strict=True):
    for s, col, tag in [(c["sig"], C_FOM, "FOM"), (c["sig_rom"], C_ROM, "qlROM")]:
        # same delay time for both portraits: the FOM's zeta, so geometry differences
        # come from the dynamics, not from a different unfolding
        Y = nt.delay_embed(s["x"], 2, c["sig"]["zeta"])[::3]
        ax.plot(Y[:, 0], Y[:, 1], lw=0.6, color=col, alpha=0.7, label=tag)
    ax.set(xlabel="$x(t)$", ylabel=r"$x(t+\zeta)$")
    ax.set_title(f"{c['name']}\n{c['regime_fom']} / {c['regime_rom']}", fontsize=9)
    ax.legend(frameon=False, fontsize=8, markerscale=4)
    ax.grid(alpha=0.25, lw=0.5)
fig.suptitle("Does the qlROM live on the same attractor? Delay portraits overlaid, one sensor, "
             "12000 closed-loop steps, FOM's $\\zeta$ for both", fontsize=11)
# its own filename: characterize() already wrote 03_characterization_panels.png (its first
# row-grid page, same stem as the PDF), and this figure would silently overwrite it
save_figure(fig)
plt.show()

print("SUMMARY")
for c in cases.values():
    print(f"  {c['name']:<12} K={c['P']['K']:2d} r={c['P']['r']:2d} | regime FOM {c['regime_fom']:>16} "
          f"-> ROM {c['regime_rom']:>16} | lam1 {c['lyap_fom']['text']} -> {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"  PDF = {pdf_path} ({n_pages} pages)")
print(f"  runtime {time.time() - t_start:.0f}s")

SUMMARY
  ks1d qp      K= 6 r=30 | regime FOM limit_cycle_period_3 -> ROM limit_cycle_period_3 | lam1 transient (9% win) -> rejected (no growth) | d 10->10, D2 0.50->1.82, std 3.29->3.16, T_ph 14.6 t.u.
  ks1d chaos   K=10 r=30 | regime FOM          chaotic -> ROM          chaotic | lam1 0.067 -> 0.151 | d 3->4, D2 2.01->2.68, std 1.36->1.28, T_ph 3.7 t.u.
  ks2d tw      K= 5 r=25 | regime FOM    quasiperiodic -> ROM    quasiperiodic | lam1 rejected (no growth) -> rejected (no growth) | d 4->4, D2 2.28->2.35, std 2.75->2.75, T_ph >= 200.0 (censored) t.u.
  ks2d chaos   K=40 r=50 | regime FOM          chaotic -> ROM          chaotic | lam1 0.303 -> 1.383 | d 3->4, D2 2.02->2.39, std 2.71->3.53, T_ph 1.8 t.u.
  PDF = ../outputs/03_characterization_panels.pdf (10 pages)
  runtime 343s

close_pdf()

11 figures written to ../outputs/03_diagnosing_a_qlrom.pdf

PosixPath('../outputs/03_diagnosing_a_qlrom.pdf')