Nonlinear time-series analysis (NTSA) with ntsa + dynamodels

Given a single long trajectory of a dynamical system — a pressure probe in a combustor, one component of a chaotic attractor — how do we decide what kind of dynamics we are looking at: fixed point, limit cycle, quasiperiodic, frequency-locked, or chaotic? This tutorial walks through the ntsa package, which packages the classical toolbox of nonlinear time-series analysis (Kantz & Schreiber, 2004) for any dynamodels-style model. The package has three modules (models come from dynamodels):

Module	Contents
ntsa.tools	trajectory runners (respawn, run_long); Takens embedding diagnostics (delay_embed, AMI optimal_lag, false_nearest_neighbours); return maps, recurrence matrices, signal statistics; classical MDS; the classify_regime decision tree; generic bifurcation_sweep + plot_bifurcation
ntsa.lyapunov	Benettin QR lyapunov_spectrum (analytic Jacobians for Lorenz63/96, finite differences for everything else); perturbation-growth leading_lyapunov; a documented stub for covariant Lyapunov vectors
ntsa.characterize	the one-call characterize() pipeline and the plotting helpers (plot_row, plot_lyapunov_fit, plot_lyapunov_spectrum, plot_mds)

Why a data-assimilation group cares. Thermoacoustic systems (Rijke tube, annular combustor, azimuthal instabilities) bifurcate from steady operation to limit-cycle, quasiperiodic, frequency-locked and chaotic oscillations as operating conditions vary. Before we assimilate data or train a reduced-order model, we must characterize the regime — and the evidence we present for each regime call is exactly the 5-panel diagnostic row of our paper supplement: time series (+ zoom inset) | power spectral density | 3-D delay portrait | first-return map | recurrence plot, backed by Lyapunov exponents. characterize() reproduces that figure for any model in one call (Section 9).

Throughout we use the Lorenz systems as testbeds, because they are cheap and every number has a known reference value: for Lorenz-63 at $\rho = 28$ the Lyapunov spectrum is $(0.906,\ 0,\ -14.57)$, and for Lorenz-96 at $F = 8$, $N_x = 10$ the leading exponent is $\lambda_1 \approx 1.2$ with three positive exponents. Every stochastic step is seeded, so the notebook self-checks with asserts on re-runs.

# ruff: noqa: E402, I001  (the env pinning below must run before numpy is imported)
import os
# Pin the BLAS thread pools BEFORE importing numpy. Everything in this module is small or
# serial linear algebra (3x3 QR steps, one 2000x2000 eigendecomposition); on a busy
# many-core machine, letting OpenBLAS spawn one thread per core makes np.linalg.eigh
# 10-40x SLOWER through contention. One thread is the fast setting here.
for var in ('OMP_NUM_THREADS', 'OPENBLAS_NUM_THREADS', 'MKL_NUM_THREADS'):
    os.environ.setdefault(var, '1')

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import Image

from ntsa import tools as ntsa_tools
from ntsa import lyapunov as lyap
from ntsa import characterize as chz

from dynamodels.physical import Lorenz63, Lorenz96

l63 = Lorenz63(rho=28., dt=0.01)
print(f"{type(l63).__name__}: rho={l63.rho}, sigma={l63.sigma}, beta={l63.beta:.4g}, dt={l63.dt}")
print(f"t_lyap = {l63.t_lyap:.4f}  (1/0.9056),  t_CR = {l63.t_CR:.4f},  t_transient = {l63.t_transient:.4f}")
print(f"Nphi = {l63.Nphi}, Nq = {l63.Nq}, obs_labels = {l63.obs_labels}")

Lorenz63: rho=28.0, sigma=10.0, beta=2.667, dt=0.01
t_lyap = 1.0905  (1/0.9056),  t_CR = 4.4170,  t_transient = 11.0424
Nphi = 3, Nq = 3, obs_labels = ['$x$', '$y$', '$z$']

The model carries its own time scales: the Lyapunov time $t_\mathrm{lyap} = 1/\lambda_1 = 1/0.9056$ (the e-folding time of infinitesimal errors — the natural unit for "how long is long"), a characteristic recurrence time t_CR, and the transient t_transient to discard before the trajectory settles on the attractor.

2. Generating clean long trajectories

Two idioms are used everywhere in this module.

respawn(model, psi0=None, dt=None, **overrides) — dynamodels models freeze their parameters at construction (alpha0) and accumulate history as they are integrated, so the clean way to change a parameter — or simply to get a pristine copy with an empty history — is to build a fresh instance. respawn reads type(model), carries over alpha0, the fixed_params (e.g. Nx for Lorenz-96) and the observation configuration (observe_dims/observed_idx), then applies your overrides on top.

run_long(model, t_run, t_transient=None) — the transient-discard idiom: integrate for t_transient (default model.t_transient), reset the history clock to $t = 0$, then record t_run time units. It returns (t, y, psi) with shapes (Nt,), (Nt, Nq) and (Nt, Nphi): the time stamps, the observable history and the full-state history of ensemble member 0 (everything in this module is strictly serial, $m = 1$). Note that run_long mutates and closes the model — always feed it a respawned copy.

l63b = ntsa_tools.respawn(l63, rho=30.)   # the override wins; everything else is carried over
print(f"respawned copy: rho={l63b.rho} (overridden), sigma={l63b.sigma}, "
      f"beta={l63b.beta:.4g}, dt={l63b.dt}, observe_dims={l63b.observe_dims}")
l63b.close()

t, y, psi = ntsa_tools.run_long(ntsa_tools.respawn(l63), t_run=60 * l63.t_lyap)
dt = float(t[1] - t[0])
x = y[:, 0]
print(f"t: {t.shape},  y: {y.shape} = (Nt, Nq),  psi: {psi.shape} = (Nt, Nphi),  dt = {dt:.3g}")

respawned copy: rho=30.0 (overridden), sigma=10.0, beta=2.667, dt=0.01, observe_dims=[0, 1, 2]
t: (6544,),  y: (6544, 3) = (Nt, Nq),  psi: (6544, 3) = (Nt, Nphi),  dt = 0.01

fig, axs = plt.subplots(3, 1, figsize=(9, 5), sharex=True, layout='constrained')
for iq, ax in enumerate(axs):
    ax.plot(t, y[:, iq], lw=0.6, color='tab:red')
    ax.set_ylabel(l63.obs_labels[iq])
    ax.margins(x=0)
axs[-1].set_xlabel('$t$')
fig.suptitle(r'Lorenz-63 ($\rho = 28$): 60 Lyapunov times, transient discarded');

3. Signal statistics

signal_stats(x, dt) bundles the standard first look at a scalar signal into one dict: the first four moments (mean, standard deviation, skewness, kurtosis), the probability density function (histogram), the power spectral density (reusing ntsa.tools.fun_PSD, the same PSD routine used across the repo), and the normalized autocorrelation function.

For a chaotic signal these carry characteristic fingerprints: a broadband PSD (continuous spectrum riding under any dominant peaks — contrast the line spectrum of a limit cycle), an ACF that decays within a few Lyapunov times (loss of memory), and a generally non-Gaussian PDF shaped by the attractor geometry — for Lorenz-63's $x$ it is bimodal, one lobe per wing of the butterfly.

stats = ntsa_tools.signal_stats(x, dt)
print({k: round(stats[k], 3) for k in ('mean', 'std', 'skew', 'kurtosis')})

fig, axs = plt.subplots(1, 3, figsize=(12, 3), layout='constrained')
axs[0].plot(stats['pdf_centers'], stats['pdf'], color='k')
axs[0].set(xlabel='$x$', ylabel='PDF')
axs[1].semilogy(stats['f'][1:], stats['psd'][1:], color='purple', lw=0.6)
axs[1].set(xlabel='frequency', ylabel='PSD', xlim=(0, 10))
axs[2].plot(stats['acf_lags'], stats['acf'], color='tab:blue', lw=0.8)
axs[2].axhline(0, color='grey', lw=0.6)
axs[2].set(xlabel='lag [t.u.]', ylabel='ACF', xlim=(0, 10 * l63.t_lyap));

{'mean': 0.332, 'std': 7.891, 'skew': -0.082, 'kurtosis': -0.577}

4. Takens delay embedding

In experiments we rarely measure the full state — often just one scalar $x(t)$. The embedding theorem (Takens 1981; see Kantz & Schreiber ch. 3 and 9) states that the delay-coordinate map

$$\mathbf{Y}(t) = \big(x(t),\ x(t+\zeta),\ \dots,\ x(t+(d-1)\zeta)\big)$$

is generically a diffeomorphic copy of the underlying attractor, provided the embedding dimension satisfies $d \ge 2D + 1$ where $D$ is the attractor's box-counting dimension. In practice a smaller, attractor-adapted $d$ usually suffices (unfolding, not the worst-case bound), and two knobs must be chosen from the data:

• the delay $\zeta$ — too small and consecutive coordinates are redundant (the portrait stretches along the diagonal); too large and they decorrelate completely (the portrait folds into a featureless blob). Fraser & Swinney (1986) proposed the first local minimum of the average mutual information (AMI) between $x(t)$ and $x(t+\zeta)$ — a nonlinear generalization of the first zero of the autocorrelation. optimal_lag implements exactly this (with an ACF fallback for noiseless periodic signals, where the AMI is flat).
• the dimension $d$ — Kennel, Brown & Abarbanel (1992): increase $d$ and count false nearest neighbours, points that are close in $d$ dimensions only because the projection folds the attractor onto itself; they fly apart when the $(d{+}1)$-th coordinate is added. The smallest $d$ with (almost) no false neighbours unfolds the attractor. false_nearest_neighbours returns that $d$ and the FNN fraction per dimension.

max_lag = max(2, int(l63.t_CR / l63.dt))
ami = ntsa_tools.average_mutual_information(x, max_lag)
zeta = ntsa_tools.optimal_lag(x, max_lag=max_lag)
print(f"zeta = {zeta} samples = {zeta * dt:.3f} t.u. = {zeta * dt / l63.t_lyap:.2f} t_lyap")

fig, ax = plt.subplots(figsize=(7, 3), layout='constrained')
lags = np.arange(1, max_lag + 1) * dt
ax.plot(lags, ami, color='k', lw=1)
ax.axvline(zeta * dt, color='tab:red', ls='--', label=rf'first local minimum: $\zeta$ = {zeta} samples')
ax.set(xlabel=r'lag $\zeta\,\Delta t$ [t.u.]', ylabel='AMI [nats]')
ax.legend();

zeta = 16 samples = 0.160 t.u. = 0.15 t_lyap

d, fractions = ntsa_tools.false_nearest_neighbours(x, zeta)
print(f"FNN embedding dimension d = {d};  fractions = {np.array2string(fractions, precision=4)}")

fig, ax = plt.subplots(figsize=(7, 3), layout='constrained')
dims = np.arange(1, len(fractions) + 1)
ax.plot(dims, fractions, 'o-', color='k')
ax.axvline(d, color='tab:red', ls='--', label=f'$d$ = {d}')
ax.axhline(0.01, color='grey', lw=0.8, ls=':', label='1% threshold')
ax.set(xlabel='embedding dimension', ylabel='FNN fraction')
ax.legend();

FNN embedding dimension d = 3;  fractions = [0.9876 0.0756 0.        nan    nan    nan    nan    nan    nan    nan]

The payoff: from the single scalar $x(t)$, the 3-D delay portrait reconstructs a diffeomorphic copy of the Lorenz butterfly — compare it with the true $(x, y, z)$ attractor from psi. Wing shapes and the fold are all recovered; only a smooth change of coordinates separates the two pictures.

Y = ntsa_tools.delay_embed(x, dim=3, lag=zeta)
print(f"embedded trajectory: {Y.shape} = (Nt - (dim-1)*lag, dim)")

fig = plt.figure(figsize=(11, 5), layout='constrained')
ax1 = fig.add_subplot(1, 2, 1, projection='3d')
ax1.plot(psi[:, 0], psi[:, 1], psi[:, 2], lw=0.4, color='k')
ax1.set(xlabel='$x$', ylabel='$y$', zlabel='$z$', title='true state space $(x, y, z)$')
ax2 = fig.add_subplot(1, 2, 2, projection='3d')
ax2.plot(Y[:, 0], Y[:, 1], Y[:, 2], lw=0.4, color='green')
ax2.set(xlabel='$x(t)$', ylabel=r'$x(t+\zeta)$', zlabel=r'$x(t+2\zeta)$',
        title='delay embedding of $x$ alone');

embedded trajectory: (6512, 3) = (Nt - (dim-1)*lag, dim)

5. First-return map and recurrence plot

First-return map. Lorenz (1963) noticed that plotting each local maximum of $z$ against the next one collapses the chaotic flow onto a nearly one-dimensional curve — the famous cusp map. That unimodal, tent-like structure is the smoking gun of low-dimensional deterministic chaos: the maxima do not scatter randomly (stochasticity) and do not collapse onto $k$ isolated points (a period-$k$ limit cycle) — they fill a continuum lying on a deterministic curve. first_return_map extracts the successive maxima, and count_peak_clusters counts distinct peak levels; when it detects a continuum band it returns the number of peaks itself as a sentinel, so a chaotic signal is never mistaken for a small-period cycle.

z = y[:, 2]
zm_i, zm_ip1, pk_idx = ntsa_tools.first_return_map(z)
k = ntsa_tools.count_peak_clusters(z[pk_idx], x_range=np.ptp(z))
print(f"{len(pk_idx)} local maxima;  count_peak_clusters -> {k} "
      f"(= n_peaks sentinel: a continuum band, not a period-k cycle)")

fig, ax = plt.subplots(figsize=(4.5, 4.5), layout='constrained')
lo, hi = zm_i.min(), zm_i.max()
ax.plot([lo, hi], [lo, hi], color='grey', lw=0.8)
ax.scatter(zm_i, zm_ip1, s=6, color='tab:blue')
ax.set(xlabel=r'$z_{\max}(i)$', ylabel=r'$z_{\max}(i+1)$', title='Lorenz cusp map')
ax.set_aspect('equal');

86 local maxima;  count_peak_clusters -> 86 (= n_peaks sentinel: a continuum band, not a period-k cycle)

Recurrence plot (Eckmann, Kamphorst & Ruelle 1987; review: Marwan et al. 2007). Using the embedded trajectory, mark all pairs of times whose states are close:

$$R_{ij} = \Theta\big(\varepsilon - \lVert \mathbf{Y}_i - \mathbf{Y}_j \rVert\big),$$

here with $\varepsilon$ = 10% of the maximum pairwise distance (recurrence_matrix). The texture of the black dots reads like a fingerprint of the regime:

• periodic — unbroken diagonal lines with evenly spaced offsets (the orbit returns exactly every period);
• quasiperiodic — diagonals with unevenly spaced offsets (two incommensurate periods never quite line up);
• chaotic — short, broken diagonals: trajectories shadow each other for a while (the determinism) and then diverge (the positive Lyapunov exponent). Isolated points would indicate noise.

n_win = 1500                       # trailing window of 15 t.u. (~3.4 recurrence times)
Y_win = ntsa_tools.delay_embed(x[-n_win:], dim=d, lag=zeta)
R = ntsa_tools.recurrence_matrix(Y_win, eps_frac=0.10)
t_w = t[-n_win:][:R.shape[0]]
print(f"recurrence matrix: {R.shape}, density = {R.mean():.3f}")

fig, ax = plt.subplots(figsize=(5, 5), layout='constrained')
ax.imshow(R, cmap='binary', origin='lower', extent=[t_w[0], t_w[-1], t_w[0], t_w[-1]])
ax.set(xlabel='$t$', ylabel='$t$', title='recurrence plot (chaotic: short, broken diagonals)');

recurrence matrix: (1468, 1468), density = 0.055

6. Lyapunov exponents

6a. Leading exponent from perturbation growth

The leading Lyapunov exponent measures the mean exponential divergence of nearby trajectories,

$$\lambda_1 = \lim_{t\to\infty} \frac{1}{t} \ln \frac{\lVert \delta(t) \rVert}{\lVert \delta(0) \rVert},$$

and $\lambda_1 > 0$ is the operational definition of chaos. leading_lyapunov implements the most model-agnostic estimator: integrate one reference trajectory and n_pert copies perturbed by $\varepsilon \sim 10^{-6}$ (serially, $m = 1$ each — it works for any model, including discrete maps like the ESN), and fit a line to the mean log-separation. The fit window is chosen automatically — from the separation minimum up to 70% of the rise towards the saturation level (once the separation reaches the attractor size, growth must stop) — and fits with $R^2 < 0.5$ are rejected (lam1 = nan), which is what happens for non-chaotic signals where there is no exponential growth to fit. Reference value: $\lambda_1 \approx 0.906$, i.e. $t_\mathrm{lyap} = 1/0.9056 \approx 1.10$ time units.

lam1, lam1_std, res = lyap.leading_lyapunov(ntsa_tools.respawn(l63), n_pert=8)
print(f"lambda_1 = {lam1:.4f} +/- {lam1_std:.4f}   (R^2 = {res['r2']:.3f}; reference 0.906)")
assert 0.6 < lam1 < 1.2, f"lam1 = {lam1} outside (0.6, 1.2)"
chz.plot_lyapunov_fit(res);

lambda_1 = 0.9469 +/- 0.0896   (R^2 = 0.889; reference 0.906)

Grey curves are the individual perturbations, black their mean, and the red dashed line the automatic fit window — note how the curve saturates at $\ln(\text{attractor size})$, which is exactly why the fit must stop before saturation.

6b. The full spectrum: Benettin's QR algorithm

A $d$-dimensional system has $d$ Lyapunov exponents, describing the growth rates of infinitesimal volumes: a $k$-dimensional parallelepiped of perturbations grows like $\exp\!\big(t \sum_{j\le k} \lambda_j\big)$. They are computed from the variational (tangent-linear) dynamics $\dot{\mathbf{U}} = J(\mathbf{u}(t))\,\mathbf{U}$ propagated alongside the state. Numerically all tangent columns collapse onto the leading direction, so Benettin et al. (1980) periodically re-orthonormalize with a QR decomposition and accumulate the logs of the diagonal of $R$:

$$\lambda_j = \frac{1}{T} \sum_{\mathrm{QR\ steps}} \ln \lvert R_{jj} \rvert .$$

lyapunov_spectrum integrates state + tangent matrix with a joint fixed-step RK4 (sharing the RK stages), QR every N_gs = 4 steps, dt = min(model.dt, 0.01). It needs the Jacobian: get_jacobian returns the analytic one for Lorenz-63/96 and falls back to central finite differences (fd_jacobian) for any other model. For Lorenz-63,

$$J(\mathbf{u}) = \begin{pmatrix} -\sigma & \sigma & 0 \\ \rho - z & -1 & -x \\ y & x & -\beta \end{pmatrix},$$

and since the RHS is quadratic, central differences should agree to near machine precision:

f = lyap.get_rhs(l63)               # f(u) -> du/dt from model.time_derivative
jac = lyap.get_jacobian(l63)        # analytic for Lorenz63/96, FD fallback otherwise
u_test = np.array([-6.0, 8.0, 27.0])
J_an, J_fd = jac(u_test), lyap.fd_jacobian(f, u_test)
rel_err = np.abs(J_fd - J_an).max() / np.abs(J_an).max()
print("analytic J(u_test) =\n", J_an)
print(f"max relative error, finite differences vs analytic: {rel_err:.2e}")
assert rel_err < 1e-6, rel_err

analytic J(u_test) =
 [[-10.          10.           0.        ]
 [  1.          -1.           6.        ]
 [  8.          -6.          -2.66666667]]
max relative error, finite differences vs analytic: 1.56e-10

exps = lyap.lyapunov_spectrum(l63)   # floor 200 t_lyap, auto-extends until converged  (~5 s)
print(f"\nBenettin spectrum: {np.array2string(exps, precision=4)}   reference: (0.906, 0, -14.57)")
assert abs(exps[0] - 0.906) < 0.1, f"lambda1 = {exps[0]}"
assert abs(exps[1]) < 2e-3, f"lambda2 = {exps[1]}"

  GS 8724/436205  lam1=0.9173

  GS 17448/436205  lam1=0.8959
  GS 26172/436205  lam1=0.9028
  GS 34896/436205  lam1=0.9062
  GS 43620/436205  lam1=0.9025
Lyapunov spectrum (T=501.4, converged): [ 9.0572e-01  7.1831e-04 -1.4573e+01]

Benettin spectrum: [ 9.0572e-01  7.1831e-04 -1.4573e+01]   reference: (0.906, 0, -14.57)

How long is long enough? Finite-time Lyapunov estimates converge slowly — the neutral (along-flow) exponent, exactly 0 for any attractor of a flow, decays only as $O(1/T)$ — and a fixed horizon cannot know when it is done. lyapunov_spectrum therefore treats t_run (default $200\,t_\mathrm{ref}$) as a floor and keeps integrating (up to t_max = 20*t_run) until every running exponent passes a halving test at two consecutive checks,

$$\lvert \lambda_j(T) - \lambda_j(T/2) \rvert < \max(\mathrm{atol},\ \mathrm{rtol}\,\lvert \lambda_j(T) \rvert), \qquad \mathrm{atol} = 2\times10^{-3},\ \mathrm{rtol} = 5\%,$$

warning if the cap is hit unconverged; full_output=True returns the running-estimate history. Two guards matter especially for limit cycles, where finite-time $\lambda_1$ is easily overestimated: (i) the warmup (_settled_state) doubles the discarded transient until a probe run shows no residual amplitude drift — a still-settling state reads as tangent growth, and near a Hopf point no fixed transient multiple is safe — and (ii) the halving test keeps extending until the neutral exponent lands within $\sim$atol of zero (VdP: $+0.005$ at the auto-chosen $T \approx 31$, vs $+0.018$ at the old fixed $T = 8$).

A free consistency check: the sum of all exponents equals the time-averaged divergence of the flow, $\sum_j \lambda_j = \langle \nabla \cdot \mathbf{f} \rangle$. For Lorenz-63 the divergence is constant, $\nabla \cdot \mathbf{f} = -(\sigma + 1 + \beta) = -13.667$ — the trace identity the spectrum must reproduce:

trace_j = -(l63.sigma + 1.0 + l63.beta)
print(f"sum(lambda) = {exps.sum():.4f}    vs    -(sigma+1+beta) = {trace_j:.4f}")
assert abs(exps.sum() - trace_j) < 0.5, f"trace identity violated: {exps.sum()} vs {trace_j}"
chz.plot_lyapunov_spectrum(exps);

sum(lambda) = -13.6666    vs    -(sigma+1+beta) = -13.6667

6c. Outlook: covariant Lyapunov vectors

The Gram–Schmidt vectors of the Benettin pass are orthogonal by construction and not covariant with the dynamics. The covariant Lyapunov vectors of Ginelli et al. (2007) — obtained by storing the $Q, R$ factors of the forward pass and back-iterating the upper-triangular expansion coefficients — give the true intrinsic stable/unstable directions (useful e.g. for projected data assimilation). The module documents this as a stub:

try:
    lyap.covariant_lyapunov_vectors()
except NotImplementedError as err:
    print(f"NotImplementedError: {err}")

NotImplementedError: CLVs: Ginelli et al. 2007 — store Q,R in forward Benettin pass, backward-iterate upper-triangular coeffs. Not implemented yet.

6d. Kaplan–Yorke dimension

The spectrum also yields a fractal-dimension estimate. With $j$ the largest index such that $\sum_{i \le j} \lambda_i \ge 0$, the Kaplan–Yorke (Lyapunov) dimension is

$$D_{KY} = j + \frac{\sum_{i \le j} \lambda_i}{\lvert \lambda_{j+1} \rvert}.$$

For Lorenz-63, $j = 2$ (since $\lambda_2 = 0$) and $D_{KY} = 2 + \lambda_1 / \lvert \lambda_3 \rvert \approx 2.06$ — the attractor is "slightly more than a surface".

D_KY = lyap.kaplan_yorke(exps)   # module implementation (tolerance-guarded neutral jitter)
print(f"D_KY(Lorenz-63) = {D_KY:.3f}   (expect ~2.06)")
assert 1.9 < D_KY < 2.2, D_KY

D_KY(Lorenz-63) = 2.062   (expect ~2.06)

6e. Torus dimension: neutral exponents, correlation dimension, Poincaré section

Along a Ruelle–Takens–Newhouse route (periodic → 2-torus → breakdown → chaos) the regime is pinned down by three independent witnesses: the number of neutral exponents in the converged spectrum ($k$ zeros ⇔ attractor $T^k$; a positive one ⇔ chaos), the Grassberger–Procaccia correlation dimension (correlation_dimension, $D_2 \approx k$ on a $k$-torus, fractal for chaos), and the plane-crossing Poincaré section (poincare_section: a section of $T^k$ has $D_2 \approx k-1$ — dots → loop → band; chaos → fractal scatter). On Lorenz-96 ($N_x = 10$) these resolve the breakdown region: $F = 4.3$ is a clean 2-torus (two zeros, $D_2 \approx 1.8$, loop section) while $F = 4.6$ is a chaotic wrinkled torus ($\lambda_1 = +0.039$ converged, $D_2 \approx 2.15$, section still loop-like) — the geometry stays toroidal after the dynamics turn chaotic, so the spectrum, not the eye, makes the call.

# correlation dimension: L63 attractor (fractal ~2.06) vs a synthetic 2-torus signal
D2_l63, _ = ntsa_tools.correlation_dimension(psi)
tt = np.arange(0, 400, 0.05)
x2 = np.sin(tt) + 0.5 * np.sin(np.sqrt(2) * tt)          # two incommensurate tones -> T^2
zeta2 = ntsa_tools.optimal_lag(x2)
D2_torus, _ = ntsa_tools.correlation_dimension(ntsa_tools.delay_embed(x2, 3, zeta2))

# plane-crossing Poincare sections: chaos -> fractal scatter; 2-torus -> closed loop
P_l63 = ntsa_tools.poincare_section(x, zeta)
P_torus = ntsa_tools.poincare_section(x2, zeta2)
D2_sec, _ = ntsa_tools.correlation_dimension(P_torus)

fig, axs = plt.subplots(1, 2, figsize=(7, 3), layout='constrained')
axs[0].scatter(P_l63[:, 0], P_l63[:, 1], s=4, color='darkorange')
axs[0].set(title=f'L63 section (attractor $D_2$={D2_l63:.2f})', xlabel='$x(t)$', ylabel=r'$x(t+\zeta)$')
axs[1].scatter(P_torus[:, 0], P_torus[:, 1], s=4, color='darkorange')
axs[1].set(title=f'2-torus section: loop ($D_2$={D2_torus:.2f}, section {D2_sec:.2f})', xlabel='$x(t)$')
print(f"D2: L63 {D2_l63:.2f} (fractal ~2.06) | two-tone torus {D2_torus:.2f} (~2) | its section {D2_sec:.2f} (~1)")
assert 1.6 < D2_torus < 2.4 and 0.6 < D2_sec < 1.4

D2: L63 1.66 (fractal ~2.06) | two-tone torus 1.98 (~2) | its section 1.02 (~1)

7. Classical MDS attractor portraits

When the state dimension is large (Lorenz-96, POD coefficients of a flow, ...) we cannot plot the attractor directly. Classical multidimensional scaling (Torgerson 1952) finds the low-dimensional coordinates that best preserve all pairwise distances: square the Euclidean distance matrix $D$, double-centre it with $C = I - \tfrac{1}{T}\mathbf{1}\mathbf{1}^\top$ to recover the Gram matrix

$$A = -\tfrac{1}{2}\, C D^{(2)} C = V \Lambda V^\top,$$

and map each snapshot to $\gamma = V_{1:k} \sqrt{\Lambda_{1:k}}$ using the top eigenpairs. For Euclidean input this is equivalent to PCA on the centred coordinates, so the $\gamma_i$ are ordered by variance. classical_mds subsamples the trajectory to at most 2000 snapshots (the Gram matrix is $T \times T$) and returns both $\gamma$ and the snapshot indices, so you can colour the portrait by time.

gamma, idx = ntsa_tools.classical_mds(psi)       # (T, 3) coordinates + subsample indices
print(f"gamma: {gamma.shape}, coordinate variances (descending): {np.round(gamma.var(axis=0), 2)}")
chz.plot_mds(gamma, t[idx]);

gamma: (2000, 3), coordinate variances (descending): [135.51  76.83   9.1 ]

8. Regime classification

classify_regime(x, dt, lam1=None, ...) distils the evidence above into a label, via a decision tree (Kantz & Schreiber ch. 1, 3, 5):

1. fixed point — the trailing fifth of the signal is flat (oscillation died out), or the full Lyapunov spectrum is entirely negative (decisive when integrator noise keeps the tail wiggling at a stable focus — e.g. Lorenz-63 at $\rho = 10$: $\lambda = [-0.6, -0.6, -12.5]$);
2. limit cycle, period-$k$ — the local maxima form $k \le k_\mathrm{max}$ tight clusters (count_peak_clusters). Checked before the Lyapunov step: $k$ tight maxima levels over a long record are incompatible with chaos, whereas a perturbation-growth $\lambda_1$ can read large and positive on a stable orbit through non-normal transient amplification (thermoacoustic limit cycles otherwise read $\lambda_1 \sim 200$);
3. chaotic — the leading Lyapunov exponent is positive beyond a trust tolerance lam_tol = max(3*lam1_std, 10/T) (a finite-time estimate on a horizon $T$ cannot resolve rates below ~$10/T$);
4. otherwise: with a full spectrum, $\ge 2$ neutral exponents ($|\lambda| <$ neutral_tol) $\to$ quasiperiodic — a 2-torus has two zero exponents, a locked periodic orbit exactly one, which is more robust than the PSD rational-ratio test (limit_denominator(10) rationals are dense: L96 at $F = 4.4$ has $f_2/f_1 = 0.552 \approx 5/9$ yet $\lambda = [0.001, 0.000, -0.003]$ is a torus). Else PSD peak analysis: after removing harmonics of the dominant peak $f_1$, is there a second incommensurate peak $f_2$? A rational ratio $f_2/f_1$ (matched to a fraction with denominator $\le 10$) means frequency-locked; an irrational ratio means quasiperiodic; no $f_2$ at all means a plain limit cycle.

It returns (label, evidence) — the evidence dict is what you show a reviewer. We test the tree on three signals: chaotic Lorenz-63 ($\rho = 28$, with the $\lambda_1$ measured in Section 6), the $\rho = 350$ periodic window of Lorenz-63 (a fast period-1 orbit, hence the finer dt = 0.005), and a synthetic two-tone signal $\sin(t) + \sin(\sqrt{2}\, t)$ whose frequency ratio is irrational by construction.

def show_evidence(ev):
    for key, val in ev.items():
        if isinstance(val, np.ndarray):
            val = (f"array of {len(val)} peaks: " + np.array2string(val[:5], precision=3)
                   + (' ...' if len(val) > 5 else ''))
        elif isinstance(val, float):
            val = round(val, 4)
        print(f"  {key:>16}: {val}")

label_ch, ev_ch = ntsa_tools.classify_regime(x, dt, lam1=lam1, lam1_std=lam1_std,
                                             t_total=t[-1] - t[0])
print(f"Lorenz-63, rho=28   ->   {label_ch}")
show_evidence(ev_ch)
assert label_ch == 'chaotic', label_ch

Lorenz-63, rho=28   ->   chaotic
           lambda1: 0.9469
       lambda1_std: 0.0896
           lam_tol: 0.2687
        n_clusters: 61
    psd_peak_freqs: array of 54 peaks: [0.031 0.076 0.153 0.199 0.29 ] ...
                f1: 0.0306
                f2: 0.5656
    rational_match: 37/2
         tail_flat: False
         n_neutral: None
             label: chaotic

l63_350 = ntsa_tools.respawn(l63, rho=350., dt=0.005)
t3, y3, _ = ntsa_tools.run_long(l63_350, t_run=10.)
label_lc, ev_lc = ntsa_tools.classify_regime(y3[:, 0], float(t3[1] - t3[0]))
print(f"Lorenz-63, rho=350   ->   {label_lc}")
show_evidence(ev_lc)
assert label_lc == 'limit_cycle_period_1', label_lc

Lorenz-63, rho=350   ->   limit_cycle_period_1
           lambda1: None
       lambda1_std: 0.0
           lam_tol: None
        n_clusters: 1
    psd_peak_freqs: array of 2 peaks: [2.603 7.708]
                f1: 2.6026
                f2: None
    rational_match: None
         tail_flat: False
         n_neutral: None
             label: limit_cycle_period_1

tt = np.arange(0, 500, 0.05)
xq = np.sin(tt) + np.sin(np.sqrt(2.0) * tt)      # two incommensurate tones
label_qp, ev_qp = ntsa_tools.classify_regime(xq, dt=0.05)
print(f"two-tone sin(t) + sin(sqrt(2) t)   ->   {label_qp}")
show_evidence(ev_qp)
assert label_qp == 'quasiperiodic', label_qp

two-tone sin(t) + sin(sqrt(2) t)   ->   quasiperiodic
           lambda1: None
       lambda1_std: 0.0
           lam_tol: 0.02
        n_clusters: 113
    psd_peak_freqs: array of 2 peaks: [0.16  0.226]
                f1: 0.16
                f2: 0.226
    rational_match: None
         tail_flat: False
         n_neutral: None
             label: quasiperiodic

fig, axs = plt.subplots(1, 3, figsize=(12, 2.6), layout='constrained')
for ax, (tv, sig, ttl) in zip(axs, [(t, x, f'rho=28: {label_ch}'),
                                    (t3, y3[:, 0], f'rho=350: {label_lc}'),
                                    (tt, xq, f'two tones: {label_qp}')]):
    n = min(len(tv), 2000)
    ax.plot(tv[-n:], sig[-n:], lw=0.6, color='tab:red')
    ax.set(title=ttl, xlabel='$t$')
    ax.margins(x=0)

Reading the evidence dicts:

• lambda1, lambda1_std, lam_tol — the Lyapunov estimate, its spread across perturbations, and the trust tolerance it must exceed for a 'chaotic' call. For $\rho = 28$, $\lambda_1 \approx 0.9 \gg$ lam_tol, so the tree stops at step 2. Where no estimate is supplied (None), the tree relies on the return map and PSD alone.
• n_clusters — the number of distinct peak levels. Exactly 1 for the $\rho = 350$ orbit (every maximum identical → period-1); for the two-tone signal the maxima form a continuum, the sentinel keeps it out of the limit-cycle branch, and the decision falls through to the PSD.
• f1, f2, psd_peak_freqs, rational_match — the dominant PSD peak, the strongest non-harmonic second peak, and whether $f_2/f_1$ matched a small rational. For the two tones $f_2/f_1 \approx \sqrt{2}$, no rational match → 'quasiperiodic'. A match like 2/3 would have meant 'frequency_locked'.
• tail_flat — the fixed-point test from step 1.

9. The one-call characterization: characterize()

Everything above is bundled into a single call. For each model, characterize: respawns it (your instance is untouched) → runs 100*t_CR past the transient (residual drift trimmed by stationary_start) → picks $\zeta$ (AMI) and $d$ (FNN) → runs leading_lyapunov and, for continuous models (spectrum='auto'), the Benettin spectrum → classify_regime → classical MDS → and writes a multi-page PDF: the 7-panel diagnostic rows (time series, PSD, delay portrait, return map, recurrence plot, 3-D MDS, Lyapunov spectrum) first, then per-case Lyapunov-fit, Lyapunov-spectrum and MDS pages. The first row-grid page is also saved as a PNG with the same name — handy because the PDF writer closes the figures. It returns one result dict per case with all the numbers (zeta, dim, lambda1, spectrum, regime, evidence, stats, gamma, ...).

We run it on the chaotic/periodic Lorenz-63 pair (the busiest cell of the notebook: two full pipelines, spectra included).

cases = [ntsa_tools.respawn(l63),                       # rho = 28 : chaotic
         ntsa_tools.respawn(l63, rho=350., dt=0.005)]   # rho = 350: period-1 limit cycle
results = chz.characterize(cases, pdf_name='figs/tutorial_l63_pair.pdf')

for r in results:
    print(f"{r['label']:<45} -> {r['regime']:<22} zeta={r['zeta']}, d={r['dim']}")
assert results[0]['regime'] == 'chaotic', results[0]['regime']
assert results[1]['regime'] == 'limit_cycle_period_1', results[1]['regime']

-- characterizing: Lorenz63 (rho=28, sigma=10, beta=2.667)

  GS 8724/436205  lam1=0.9173
  GS 17448/436205  lam1=0.8959
  GS 26172/436205  lam1=0.9028
  GS 34896/436205  lam1=0.9062
  GS 43620/436205  lam1=0.9025
Lyapunov spectrum (T=501.4, converged): [ 9.0572e-01  7.1831e-04 -1.4573e+01]

-- characterizing: Lorenz63 (rho=350, sigma=10, beta=2.667)

  GS 4416/883392  lam1=0.0100
  GS 8832/883392  lam1=0.0272
  GS 13248/883392  lam1=0.0008
  GS 17664/883392  lam1=0.0021
  GS 22080/883392  lam1=0.0108
  GS 26496/883392  lam1=0.0007
  GS 30912/883392  lam1=0.0010
  GS 35328/883392  lam1=0.0067
  GS 39744/883392  lam1=0.0007
  GS 44160/883392  lam1=0.0005
  GS 48576/883392  lam1=0.0049
  GS 52992/883392  lam1=0.0007

  GS 57408/883392  lam1=0.0003
  GS 61824/883392  lam1=0.0038
  GS 66240/883392  lam1=0.0007
  GS 70656/883392  lam1=0.0001
  GS 75072/883392  lam1=0.0031
  GS 79488/883392  lam1=0.0007
Lyapunov spectrum (T=397.44, converged): [ 6.8799e-04 -7.3134e-01 -1.2936e+01]

Saved figures --> figs/tutorial_l63_pair.pdf
Lorenz63 (rho=28, sigma=10, beta=2.667)       -> chaotic                zeta=16, d=3
Lorenz63 (rho=350, sigma=10, beta=2.667)      -> limit_cycle_period_1   zeta=20, d=2

Image('figs/tutorial_l63_pair.png')    # first page of figs/tutorial_l63_pair.pdf

How to read a row, left to right (this is the layout of the paper-supplement figure):

1. time series (red) with a zoom inset spanning ~5 dominant periods — the row title carries the regime call, $\zeta$, $d$ and $\lambda_1$;
2. PSD (purple, semilogy) — broadband floor for chaos vs a clean line spectrum for the limit cycle. The black triangles are the peaks found by classify_regime, and the text box lists $f_1$, $f_2$ and their ratio — precisely the numbers that drive the LC / frequency-locked / QP / chaotic call, so the classification can be audited by eye;
3. 3-D delay portrait (green) — butterfly vs closed loop;
4. first-return map (blue) — continuum on the cusp curve vs a single point on the identity line;
5. recurrence plot — broken short diagonals vs evenly spaced unbroken diagonals.

10. A higher-dimensional test: Lorenz-96

The Lorenz-96 lattice,

$$\dot{x}_i = (x_{i+1} - x_{i-2})\, x_{i-1} - x_i + F, \qquad i = 0, \dots, N_x - 1 \ \ (\text{cyclic}),$$

is the standard data-assimilation testbed: advection-like quadratic coupling, linear damping, constant forcing $F$. With $N_x = 10$ and $F = 8$ it is chaotic with several positive Lyapunov exponents ($\lambda_1 \approx 1.2$, three positive) — genuinely higher-dimensional chaos, unlike Lorenz-63. All tools above apply unchanged.

l96 = Lorenz96(Nx=10, F=8., dt=0.01)
t96, y96, psi96 = ntsa_tools.run_long(ntsa_tools.respawn(l96), t_run=60 * l96.t_lyap)
dt96 = float(t96[1] - t96[0])
x96 = y96[:, 0]
print(f"t: {t96.shape}, y: {y96.shape}, psi: {psi96.shape};  observed sites: {l96.observed_idx}")

fig, (ax0, ax1) = plt.subplots(2, 1, figsize=(9, 5), sharex=True, layout='constrained',
                               gridspec_kw=dict(height_ratios=[1, 1.6]))
ax0.plot(t96, x96, lw=0.6, color='tab:red')
ax0.set_ylabel(l96.obs_labels[0])
ax0.margins(x=0)
im = ax1.imshow(psi96.T, aspect='auto', origin='lower', cmap='RdBu_r',
                extent=[t96[0], t96[-1], -0.5, l96.Nx - 0.5])
ax1.set(xlabel='$t$', ylabel='site $i$')
fig.colorbar(im, ax=ax1, label='$x_i$');

t: (5069,), y: (5069, 3), psi: (5069, 10);  observed sites: [0, 5, 9]

zeta96 = ntsa_tools.optimal_lag(x96, max_lag=max(2, int(l96.t_CR / l96.dt)))
d96, fr96 = ntsa_tools.false_nearest_neighbours(x96, zeta96)
print(f"zeta = {zeta96} samples ({zeta96 * dt96:.2f} t.u.);  FNN dimension d = {d96}")
print(f"FNN fractions: {np.array2string(fr96, precision=3)}")

zeta = 26 samples (0.26 t.u.);  FNN dimension d = 4
FNN fractions: [0.996 0.521 0.069 0.004   nan   nan   nan   nan   nan   nan]

The full Benettin spectrum with n_exp = 10 (all of them — the analytic cyclic-banded Jacobian makes this cheap). Two things to verify: $\lambda_1 \approx 1.2$ with three positive exponents, and the trace identity — for Lorenz-96 every diagonal Jacobian entry is $-1$, so $\sum_j \lambda_j = -N_x = -10$ exactly.

exps96 = lyap.lyapunov_spectrum(l96, n_exp=10)     # ~40 s (auto-extends until converged)
n_pos = int((exps96 > 0.05).sum())
print(f"\nlambda_1 = {exps96[0]:.3f} (expect ~1.2);  {n_pos} positive exponents (expect 3)")
print(f"sum(lambda) = {exps96.sum():.3f}  vs  -Nx = {-l96.Nx}")
assert abs(exps96[0] - 1.2) < 0.4, exps96[0]
assert 2 <= n_pos <= 4, n_pos
assert abs(exps96.sum() + l96.Nx) < 0.5, exps96.sum()
D_KY96 = lyap.kaplan_yorke(exps96)
print(f"Kaplan-Yorke dimension = {D_KY96:.2f}  (of Nphi = {l96.Nphi})")
chz.plot_lyapunov_spectrum(exps96);

  GS 6756/337837  lam1=1.3522

  GS 13512/337837  lam1=1.2854

  GS 20268/337837  lam1=1.2232

  GS 27024/337837  lam1=1.1832

  GS 33780/337837  lam1=1.2034

  GS 40536/337837  lam1=1.2085

  GS 47292/337837  lam1=1.2012

  GS 54048/337837  lam1=1.1979

  GS 60804/337837  lam1=1.2098

  GS 67560/337837  lam1=1.2008

  GS 74316/337837  lam1=1.2158

  GS 81072/337837  lam1=1.2004

  GS 87828/337837  lam1=1.1910

  GS 94584/337837  lam1=1.1963

  GS 101340/337837  lam1=1.2018

  GS 108096/337837  lam1=1.2009

  GS 114852/337837  lam1=1.2064

  GS 121608/337837  lam1=1.2024

  GS 128364/337837  lam1=1.2008

  GS 135120/337837  lam1=1.2005

  GS 141876/337837  lam1=1.1966

  GS 148632/337837  lam1=1.1875

  GS 155388/337837  lam1=1.1864

  GS 162144/337837  lam1=1.1891

  GS 168900/337837  lam1=1.1891

  GS 175656/337837  lam1=1.1867

  GS 182412/337837  lam1=1.1900

  GS 189168/337837  lam1=1.1839

  GS 195924/337837  lam1=1.1884

  GS 202680/337837  lam1=1.1913

  GS 209436/337837  lam1=1.1845

  GS 216192/337837  lam1=1.1855

  GS 222948/337837  lam1=1.1791

  GS 229704/337837  lam1=1.1777

  GS 236460/337837  lam1=1.1806

  GS 243216/337837  lam1=1.1797

  GS 249972/337837  lam1=1.1809

  GS 256728/337837  lam1=1.1841

  GS 263484/337837  lam1=1.1888

  GS 270240/337837  lam1=1.1901

  GS 276996/337837  lam1=1.1919

  GS 283752/337837  lam1=1.1869

  GS 290508/337837  lam1=1.1819

  GS 297264/337837  lam1=1.1846

  GS 304020/337837  lam1=1.1848

  GS 310776/337837  lam1=1.1861

  GS 317532/337837  lam1=1.1843

  GS 324288/337837  lam1=1.1833

Lyapunov spectrum (T=3274.72, converged): [ 1.1809e+00  6.9042e-01  8.4933e-02 -4.5108e-03 -4.5208e-01 -8.7725e-01
 -1.3333e+00 -1.9145e+00 -2.7762e+00 -4.5984e+00]

lambda_1 = 1.181 (expect ~1.2);  3 positive exponents (expect 3)
sum(lambda) = -10.000  vs  -Nx = -10
Kaplan-Yorke dimension = 6.47  (of Nphi = 10)

res96 = chz.characterize(ntsa_tools.respawn(l96), spectrum=False,   # spectrum computed above
                         pdf_name='figs/tutorial_l96.pdf')
r = res96[0]
print(f"{r['label']} -> {r['regime']}  (zeta={r['zeta']}, d={r['dim']}, lambda1={r['lambda1']:.3f})")
assert r['regime'] == 'chaotic', r['regime']
Image('figs/tutorial_l96.png')

-- characterizing: Lorenz96 (F=8)

Saved figures --> figs/tutorial_l96.pdf
Lorenz96 (F=8) -> chaotic  (zeta=39, d=5, lambda1=1.018)

A coarse bifurcation diagram in $F$

bifurcation_sweep is the module's generic route-to-chaos tool: for each parameter value it integrates past a transient and collects the local extrema of every observable. By default the whole sweep runs as one ensemble forecast — the swept parameter is augmented into the state (init_ensemble(est_alpha=[param])), so member $k$ integrates at values[k] from a shared initial state and the parameter values run in parallel. Plotted against the parameter, a single point per column means a period-1 orbit, $k$ points a period-$k$ orbit, and a filled band a chaotic (or quasiperiodic) attractor. With continuation=True the sweep falls back to a serial loop where each run starts from the previous endpoint (plus a tiny perturbation), which keeps the sweep on the same attractor branch — the standard numeric continuation trick, useful across hysteresis, but inherently sequential.

We sweep $F \in [2, 18]$ coarsely (28 values, short samples — the whole sweep takes a few seconds). The sweep starts just past the Hopf bifurcation of the fixed point $x_i = F$: at $F \approx 2$–$4$ the lattice carries a periodic travelling wave — a few isolated peak levels per column — and by $F \gtrsim 5$–$6$ the columns fill out into the bands of chaos, widening as the forcing grows.

F_values = np.linspace(2., 18., 28)
vals, peaks = ntsa_tools.bifurcation_sweep(l96, 'F', F_values,
                                           t_transient=10., t_sample=30.,
                                           extrema=('max', 'min'))
fig, axs = ntsa_tools.plot_bifurcation(vals, peaks, l96.alpha_labels['F'], l96.obs_labels)

F sweep:   0%|          | 0/28 [00:00<?, ?it/s]

F sweep:  11%|█         | 3/28 [00:00<00:01, 23.54it/s]

F sweep:  21%|██▏       | 6/28 [00:00<00:01, 18.17it/s]

F sweep:  29%|██▊       | 8/28 [00:00<00:01, 14.46it/s]

F sweep:  36%|███▌      | 10/28 [00:00<00:01, 11.90it/s]

F sweep:  43%|████▎     | 12/28 [00:01<00:01,  9.89it/s]

F sweep:  50%|█████     | 14/28 [00:01<00:01,  8.75it/s]

F sweep:  54%|█████▎    | 15/28 [00:01<00:01,  8.23it/s]

F sweep:  57%|█████▋    | 16/28 [00:01<00:01,  7.78it/s]

F sweep:  61%|██████    | 17/28 [00:01<00:01,  7.37it/s]

F sweep:  64%|██████▍   | 18/28 [00:01<00:01,  6.91it/s]

F sweep:  68%|██████▊   | 19/28 [00:02<00:01,  6.61it/s]

F sweep:  71%|███████▏  | 20/28 [00:02<00:01,  6.34it/s]

F sweep:  75%|███████▌  | 21/28 [00:02<00:01,  6.04it/s]

F sweep:  79%|███████▊  | 22/28 [00:02<00:01,  5.74it/s]

F sweep:  82%|████████▏ | 23/28 [00:02<00:00,  5.52it/s]

F sweep:  86%|████████▌ | 24/28 [00:03<00:00,  5.36it/s]

F sweep:  89%|████████▉ | 25/28 [00:03<00:00,  5.21it/s]

F sweep:  93%|█████████▎| 26/28 [00:03<00:00,  4.94it/s]

F sweep:  96%|█████████▋| 27/28 [00:03<00:00,  4.81it/s]

F sweep: 100%|██████████| 28/28 [00:03<00:00,  4.73it/s]

F sweep: 100%|██████████| 28/28 [00:03<00:00,  7.12it/s]

A closing note on scope: the sweep is completely generic — it only needs model.params, respawn and the ensemble machinery every dynamodels model already carries, so the same two lines produce bifurcation diagrams for VdP (e.g. over $\beta$), Rijke, Annular, or any other model. For IVP models the ensemble mode uses a multiprocessing pool, so plain scripts (not notebooks) should guard the call with if __name__ == '__main__': under forkserver/spawn.

11. Summary

Question about a signal	Tool	Key reference
clean trajectory on the attractor	respawn + run_long	—
moments / PDF / PSD / ACF	signal_stats	Kantz & Schreiber (2004) ch. 2
embedding delay $\zeta$	optimal_lag (AMI minimum)	Fraser & Swinney (1986)
embedding dimension $d$	false_nearest_neighbours	Kennel et al. (1992)
phase portrait from one scalar	delay_embed	Takens (1981)
periodicity of the orbit	first_return_map + count_peak_clusters	Lorenz (1963)
determinism / regime texture	recurrence_matrix	Eckmann et al. (1987); Marwan et al. (2007)
is it chaotic? ($\lambda_1$)	leading_lyapunov	Kantz & Schreiber (2004) ch. 5
full spectrum, $D_{KY}$, trace check	lyapunov_spectrum (+ get_rhs/get_jacobian)	Benettin et al. (1980)
stable/unstable directions	covariant_lyapunov_vectors (stub)	Ginelli et al. (2007)
portrait of a high-dim state	classical_mds	Torgerson (1952)
the regime label + evidence	classify_regime	Kantz & Schreiber (2004)
route to chaos in a parameter	bifurcation_sweep + plot_bifurcation	—
all of the above, one PDF	characterize	this module

print('Measured in this notebook (all seeded, all assert-checked):')
print(f'  L63 rho=28  leading exponent (perturbation growth): {lam1:.3f} +/- {lam1_std:.3f}')
print(f'  L63 rho=28  Benettin spectrum: {np.array2string(exps, precision=3)}'
      f'   [sum = {exps.sum():.3f} vs -(sigma+1+beta) = {trace_j:.3f}]')
print(f'  L63 rho=28  Kaplan-Yorke dimension: {D_KY:.3f}')
print(f'  L96 F=8, Nx=10  spectrum: lambda1 = {exps96[0]:.3f}, {n_pos} positive, '
      f'sum = {exps96.sum():.2f} vs -Nx = {-l96.Nx}, D_KY = {D_KY96:.2f}')
print(f'  regimes: rho=28 -> {label_ch};  rho=350 -> {label_lc};  two-tone -> {label_qp}')

Measured in this notebook (all seeded, all assert-checked):
  L63 rho=28  leading exponent (perturbation growth): 0.947 +/- 0.090
  L63 rho=28  Benettin spectrum: [ 9.057e-01  7.183e-04 -1.457e+01]   [sum = -13.667 vs -(sigma+1+beta) = -13.667]
  L63 rho=28  Kaplan-Yorke dimension: 2.062
  L96 F=8, Nx=10  spectrum: lambda1 = 1.181, 3 positive, sum = -10.00 vs -Nx = -10, D_KY = 6.47
  regimes: rho=28 -> chaotic;  rho=350 -> limit_cycle_period_1;  two-tone -> quasiperiodic

References

• Kantz, H. & Schreiber, T. (2004). Nonlinear Time Series Analysis, 2nd ed., Cambridge University Press.
• Lorenz, E. N. (1963). Deterministic nonperiodic flow. J. Atmos. Sci. 20, 130–141.
• Takens, F. (1981). Detecting strange attractors in turbulence. Lecture Notes in Math. 898.
• Fraser, A. M. & Swinney, H. L. (1986). Independent coordinates for strange attractors from mutual information. Phys. Rev. A 33, 1134.
• Kennel, M. B., Brown, R. & Abarbanel, H. D. I. (1992). Determining embedding dimension for phase-space reconstruction using a geometrical construction. Phys. Rev. A 45, 3403.
• Benettin, G., Galgani, L., Giorgilli, A. & Strelcyn, J.-M. (1980). Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems. Meccanica 15, 9–30.
• Eckmann, J.-P., Kamphorst, S. O. & Ruelle, D. (1987). Recurrence plots of dynamical systems. Europhys. Lett. 4, 973.
• Marwan, N., Romano, M. C., Thiel, M. & Kurths, J. (2007). Recurrence plots for the analysis of complex systems. Phys. Rep. 438, 237–329.
• Torgerson, W. S. (1952). Multidimensional scaling: I. Theory and method. Psychometrika 17, 401–419.
• Ginelli, F., Poggi, P., Turchi, A., Chaté, H., Livi, R. & Politi, A. (2007). Characterizing dynamics with covariant Lyapunov vectors. Phys. Rev. Lett. 99, 130601.
• Lorenz, E. N. (1996). Predictability: a problem partly solved. Proc. ECMWF Seminar on Predictability, Vol. 1, 1–18.