Tutorial 4: charts, transitions, and the exact atlas¶
This tutorial is about the geometry layer of qlroms -- Chart, Atlas,
TransitionMaps -- and about one specific claim those classes make:
mapping a reduced state local -> atlas -> local is lossless, i.e. the exact identity up to floating point (~$10^{-15}$), whereas the naive local -> local round trip through a neighbouring chart is not.
No dynamics appear anywhere below. A qlGalerkin adds a
vector field on top of this geometry; everything here is true before a single
time step is taken.
Notation follows docs/theory/atlas.txt:
| symbol | meaning | code |
|---|---|---|
| $\bm{q}\in\mathbb{R}^{N_h}$ | physical state | X[:, n] |
| $\bar{\bm{q}}_k$, $\bm{\Phi}_k$ | chart $k$ centroid and local basis | atlas.centroids[k], atlas.Phi_all[:, :, k] |
| $\bm{a}_k\in\mathbb{R}^{r}$ | local reduced coordinates | a |
| $\mathbf{M}$ | inner-product weight | atlas.Mw (None = identity) |
| $\bar{\bm{q}}_g$, $\bm{\Phi}_g$ | atlas reference state and basis | atlas.qbar_g, atlas.Phi_g |
| $\bm{z}\in\mathbb{R}^{r_g}$ | atlas (global) coordinates | z |
| $\mathbf{T}_{ji},\bm{d}_{ji}$ | chart $i\to j$ affine map | tmap[i, j], tshift[i, j] |
| $\mathbf{T}_{gk},\bm{d}_{gk}$ | chart $k\to$ atlas | atlas.Tgk[k], atlas.dgk[k] |
| $\mathbf{T}_{kg},\bm{d}_{kg}$ | atlas $\to$ chart $k$ | atlas.Tkg[k], atlas.dkg[k] |
The index convention is target first, source second: $\mathbf{T}_{ji}$ maps from $i$ to $j$.
import time
import numpy as np
import torch
from dynamodels.physical import KS
from qlroms import Chart, compute_pod_basis, fit_charts
from qlroms.atlas import build_global_assimilation_basis
from qlroms.utils.plots import close_pdf, save_figure, start_pdf
from qlroms.transitions import TransitionMaps
# the summary figure of section 5 (tutorials/plots_04.py)
from plots_04 import plot_summary
t_start = time.time()
torch.manual_seed(0)
# every figure below also becomes a page of figs/04_atlas_and_transitions.pdf
start_pdf("04_atlas_and_transitions")
PosixPath('figs/04_atlas_and_transitions.pdf')
1. Charts¶
The 1-D Kuramoto-Sivashinsky equation in its chaotic regime, sampled with a coarse
dt, because the geometry needs only snapshots that visit several regimes, not a
time-resolved trajectory.
fit_charts is the one-call offline geometry build: k-means on the snapshots,
a POD basis per cluster, and (by default, build_atlas=True) the exact shared
atlas. The result is an Atlas, which owns no dynamics at all.
fom = KS(Nx=128, dt=0.25) # nu = 1 on L = 2*pi/sqrt(0.08) -- chaotic
psi, t = fom.time_integrate(Nt=int(fom.t_transient / fom.dt))
fom.update_history(psi, t) # discard the transient
Nt = 6000
psi, t = fom.time_integrate(Nt=Nt)
fom.update_history(psi, t)
u = fom.get_observable_hist(Nt=Nt, loc="all") # (Nt, Nx, m)
X = torch.as_tensor(u[:, :, 0].T, dtype=torch.float64) # (N, Nt), columns = time
Ntrain = 5000
Xtrain, Xtest = X[:, :Ntrain], X[:, Ntrain:]
K, r = 6, 10
atlas = fit_charts(Xtrain, K=K, r=r, dt=fom.dt, random_state=0)
print(f"Atlas: N = {atlas.N} dofs, K = {atlas.K} charts, r = {atlas.r} modes/chart, "
f"r_g = {atlas.Phi_g.shape[1]} atlas dims, Mw = {atlas.Mw}")
Atlas: N = 128 dofs, K = 6 charts, r = 10 modes/chart, r_g = 20 atlas dims, Mw = None
1.1. A chart is an affine space¶
Chart $k$ approximates the state as
$$\bm{q}\;\approx\;\bar{\bm{q}}_k+\bm{\Phi}_k\bm{a}_k ,\qquad \bm{a}_k=\bm{\Phi}_k^\top\mathbf{M}\left(\bm{q}-\bar{\bm{q}}_k\right),$$
with $\bm{\Phi}_k^\top\mathbf{M}\bm{\Phi}_k=\mathbf{I}_r$. The set
$$\mathcal{A}_k=\bar{\bm{q}}_k+\operatorname{span}(\bm{\Phi}_k)$$
is chart $k$'s affine trial space -- an $r$-dimensional flat living in
$\mathbb{R}^{N_h}$, offset from the origin by its centroid. Chart is exactly
this pair $(\bar{\bm{q}}_k,\bm{\Phi}_k)$ plus the two maps.
Projection is a lossy operation: $\bm{P}_k=\bm{\Phi}_k\bm{\Phi}_k^\top\mathbf{M}$ is the $\mathbf{M}$-orthogonal projector onto $\operatorname{span}(\bm{\Phi}_k)$, so recovering a projected state discards whatever was orthogonal to the chart. Keep that number in mind -- it is the only lossy step in this entire tutorial.
k0 = 2
chart = Chart(atlas.Phi_all[:, :, k0], atlas.centroids[k0], Mw=atlas.Mw)
q = Xtest[:, 0]
a = chart.project_state(q) # (r, 1)
q_rec = chart.recover_state(a).reshape(-1)
proj_err = (q_rec - q).norm() / q.norm()
# idempotence: projecting an already-projected state changes nothing
a_again = chart.project_state(q_rec)
print(f"physical -> local -> physical, chart {k0}: relative error = {proj_err:.2e} (LOSSY: truncation)")
print(f"P_k is a projector, P_k^2 = P_k: |a(P q) - a(q)| = {float((a_again - a).abs().max()):.2e}")
physical -> local -> physical, chart 2: relative error = 5.03e-02 (LOSSY: truncation) P_k is a projector, P_k^2 = P_k: |a(P q) - a(q)| = 8.88e-15
1.2. Chart to chart: the pairwise transition¶
When a trajectory leaves chart $i$ and enters chart $j$, the reduced coordinates must be re-expressed. Substituting chart $i$'s reconstruction into chart $j$'s projection gives
$$\bm{a}_j=\bm{\Phi}_j^\top\mathbf{M}\bm{\Phi}_i\,\bm{a}_i +\bm{\Phi}_j^\top\mathbf{M}\left(\bar{\bm{q}}_i-\bar{\bm{q}}_j\right) \;=\;\mathbf{T}_{ji}\bm{a}_i+\bm{d}_{ji}.$$
Atlas precomputes all $K^2$ pairs once, offline, into tmap $(K,K,r,r)$ and
tshift $(K,K,r)$. At run time a switch is a single $(r,r)\!\times\!(r,)$ product --
no $N_h$-dimensional vector is ever formed. That is the entire point of caching
them.
pair = TransitionMaps.from_model(atlas, method="pairwise")
i, j = 0, 3
a_i = atlas.project_state(Xtest[:, 0], cluster_id=i)[:-1, 0] # drop the id row
chart_i = Chart(atlas.Phi_all[:, :, i], atlas.centroids[i], Mw=atlas.Mw)
chart_j = Chart(atlas.Phi_all[:, :, j], atlas.centroids[j], Mw=atlas.Mw)
a_j_cached = pair.map(a_i, i, j) # small: (r,r) @ (r,)
a_j_physical = chart_j.project_state(chart_i.recover_state(a_i)).reshape(-1) # the N_h round trip
print(f"tmap/tshift reproduce the physical round trip exactly: "
f"|cached - physical| = {float((a_j_cached - a_j_physical).abs().max()):.2e}")
tmap/tshift reproduce the physical round trip exactly: |cached - physical| = 3.55e-15
1.3 ... and why it does not come back¶
The pairwise map is exact as a definition -- it is by construction chart $j$'s projection of chart $i$'s own reconstruction. But it is a projection, so going there and back composes two of them:
$$\mathbf{T}_{ij}\mathbf{T}_{ji} =\bm{\Phi}_i^\top\mathbf{M}\bm{\Phi}_j\bm{\Phi}_j^\top\mathbf{M}\bm{\Phi}_i \;\neq\;\mathbf{I}_r \qquad\text{unless }\operatorname{span}(\bm{\Phi}_i)\subseteq\operatorname{span}(\bm{\Phi}_j).$$
Whatever of chart $i$'s space is invisible to chart $j$ is gone for good. The cell below measures this loss on the KS test states, for every ordered pair of charts; compare it with round-off.
apod = atlas.project_state(Xtest)
A_test, ids = apod[:-1], apod[-1].long() # (r, Nt_test) coords + their charts
rt_pair = np.full((K, K), np.nan)
for kk in range(K):
sel = ids == kk
a_k = A_test[:, sel]
for jj in range(K):
if jj == kk:
continue
back = pair.map(pair.map(a_k, kk, jj), jj, kk)
rt_pair[kk, jj] = float(((back - a_k).norm(dim=0) / a_k.norm(dim=0)).mean())
print("pairwise k -> j -> k relative error, mean over test states in chart k:")
print(f" min {np.nanmin(rt_pair):.3f} | mean {np.nanmean(rt_pair):.3f} | max {np.nanmax(rt_pair):.3f}")
pairwise k -> j -> k relative error, mean over test states in chart k: min 0.027 | mean 0.138 | max 0.376
2. The atlas: one affine space that contains all of them¶
The fix is not a better pair of charts, it is a third space large enough to hold every chart exactly. Take the global reference state $\bar{\bm{q}}_g=\frac{1}{M}\sum_i\bm{q}_i$ (the training-snapshot mean) and assemble the dictionary
$$\bm{U}=\left[\bm{\Phi}_1,\dots,\bm{\Phi}_K,\; \bar{\bm{q}}_1-\bar{\bm{q}}_g,\dots,\bar{\bm{q}}_K-\bar{\bm{q}}_g\right],$$
then take an $\mathbf{M}$-orthonormal basis of its range,
$$\bm{\Phi}_g=\operatorname{orth}_{\mathbf{M}}(\bm{U}),\qquad \bm{\Phi}_g^\top\mathbf{M}\bm{\Phi}_g=\mathbf{I}_{r_g},$$
computed, for a general $\mathbf{M}$, as $\bm{G}=\bm{U}^\top\mathbf{M}\bm{U}=\bm{V}\bm{\Lambda}\bm{V}^\top$,
$\bm{\Phi}_g=\bm{U}\bm{V}_r\bm{\Lambda}_r^{-1/2}$ (a QR factorization of $\bm{U}$ when
$\mathbf{M}$ is the identity or diagonal), after dropping the numerically negligible
directions (tol). That is
build_global_assimilation_basis, and $\mathcal{A}_g=\bar{\bm{q}}_g+
\operatorname{span}(\bm{\Phi}_g)$ is the atlas.
Inclusion is by construction. For any state on chart $k$,
$$\bm{q}-\bar{\bm{q}}_g=\underbrace{(\bar{\bm{q}}_k-\bar{\bm{q}}_g)}_{\text{column of }\bm{U}} +\underbrace{\bm{\Phi}_k\bm{a}_k}_{\in\,\operatorname{span}(\bm{\Phi}_k)\subseteq\operatorname{range}(\bm{U})},$$
so $\mathcal{A}_k\subseteq\mathcal{A}_g$ for every $k$: both the centroid offsets and the local bases were put in the dictionary. Note that $r_g$ is discovered, not chosen: when the charts share directions, the atlas is smaller than the naive column count $Kr+K$. The cell below prints both.
rg = atlas.Phi_g.shape[1]
G = atlas.Phi_g.T @ atlas.Phi_g
resid = (atlas.centroids - atlas.qbar_g) - atlas.dgk @ atlas.Phi_g.T # c_k on the atlas?
print(f"dictionary columns K*r + K = {K * r + K} -> numerical rank r_g = {rg} "
f"({(K * r + K) / rg:.1f}x compression: the charts overlap)")
print(f"M-orthonormality |Phi_g^T M Phi_g - I| = {float((G - torch.eye(rg, dtype=G.dtype)).abs().max()):.2e}")
print(f"centroids lie on the atlas: max relative residual = "
f"{float((resid.norm(dim=1) / (atlas.centroids - atlas.qbar_g).norm(dim=1)).max()):.2e}")
dictionary columns K*r + K = 66 -> numerical rank r_g = 20 (3.3x compression: the charts overlap) M-orthonormality |Phi_g^T M Phi_g - I| = 8.88e-16 centroids lie on the atlas: max relative residual = 2.91e-12
2.1. Transitions in the atlas¶
Two views of the same switch: first the round trip that proves it is lossless, then the switch written chart-to-chart through the shared frame.
a. The lossless local -> atlas -> local transition¶
With the atlas in place, chart $k$ is just another pair of affine maps -- same $\mathbf{T},\bm{d}$ notation, with $g$ as the target/source index:
$$\bm{z}=\mathbf{T}_{gk}\bm{a}_k+\bm{d}_{gk},\qquad \bm{a}_k=\mathbf{T}_{kg}\bm{z}+\bm{d}_{kg},$$
$$\mathbf{T}_{gk}=\bm{\Phi}_g^\top\mathbf{M}\bm{\Phi}_k,\quad \bm{d}_{gk}=\bm{\Phi}_g^\top\mathbf{M}(\bar{\bm{q}}_k-\bar{\bm{q}}_g),\quad \mathbf{T}_{kg}=\bm{\Phi}_k^\top\mathbf{M}\bm{\Phi}_g,\quad \bm{d}_{kg}=\bm{\Phi}_k^\top\mathbf{M}(\bar{\bm{q}}_g-\bar{\bm{q}}_k).$$
Composing them,
$$\widehat{\bm{a}}_k=\mathbf{T}_{kg}\left(\mathbf{T}_{gk}\bm{a}_k+\bm{d}_{gk}\right)+\bm{d}_{kg} =\mathbf{T}_{kg}\mathbf{T}_{gk}\,\bm{a}_k+\left(\mathbf{T}_{kg}\bm{d}_{gk}+\bm{d}_{kg}\right),$$
so $\widehat{\bm{a}}_k=\bm{a}_k$ for all $\bm{a}_k$ if and only if the two conditions
$$\boxed{\;\mathbf{T}_{kg}\mathbf{T}_{gk}=\mathbf{I}_r\;},\qquad \boxed{\;\mathbf{T}_{kg}\bm{d}_{gk}+\bm{d}_{kg}=\bm{0}\;}$$
hold. The first says the linear parts undo each other, the second that the two offsets cancel. Both follow from $\mathcal{A}_k\subseteq\mathcal{A}_g$: because $\bm{\Phi}_g\bm{\Phi}_g^\top\mathbf{M}$ acts as the identity on $\operatorname{span}(\bm{\Phi}_k)$,
$$\mathbf{T}_{kg}\mathbf{T}_{gk} =\bm{\Phi}_k^\top\mathbf{M}\underbrace{\bm{\Phi}_g\bm{\Phi}_g^\top\mathbf{M}\bm{\Phi}_k}_{=\;\bm{\Phi}_k} =\bm{\Phi}_k^\top\mathbf{M}\bm{\Phi}_k=\mathbf{I}_r,$$
and the same substitution turns $\mathbf{T}_{kg}\bm{d}_{gk}$ into $\bm{\Phi}_k^\top\mathbf{M}(\bar{\bm{q}}_k-\bar{\bm{q}}_g)=-\bm{d}_{kg}$. Let us watch both hold to machine precision.
I_r = torch.eye(r, dtype=atlas.Tgk.dtype)
print("chart | ||T_kg T_gk - I||_inf ||T_kg d_gk + d_kg||_inf")
for k in range(K):
e_lin = float((atlas.Tkg[k] @ atlas.Tgk[k] - I_r).abs().max())
e_off = float((atlas.Tkg[k] @ atlas.dgk[k] + atlas.dkg[k]).abs().max())
print(f" {k} | {e_lin:.2e} {e_off:.2e}")
chart | ||T_kg T_gk - I||_inf ||T_kg d_gk + d_kg||_inf 0 | 2.22e-15 8.88e-15 1 | 1.33e-15 2.66e-15 2 | 1.78e-15 3.55e-15 3 | 2.22e-15 2.55e-15 4 | 2.00e-15 2.66e-15 5 | 1.78e-15 3.55e-15
The class-level surface for the same statement: TransitionMaps.to_atlas /
from_atlas. Run every test snapshot through its own chart -> atlas -> back and
compare with the pairwise round trip of section 1.3.
tm = atlas.transitions # method='auto' -> the exact atlas, because it exists
print(f"transitions method = {tm.method!r}, has_atlas = {tm.has_atlas}, has_pairwise = {tm.has_pairwise}")
rt_atlas = np.zeros(K)
for k in range(K):
a_k = A_test[:, ids == k]
z = tm.to_atlas(a_k, k) # (r_g, m)
back = tm.from_atlas(z, k) # (r, m)
rt_atlas[k] = float(((back - a_k).norm(dim=0) / a_k.norm(dim=0)).mean())
print(f"local -> atlas -> local : relative error {rt_atlas.max():.2e} (LOSSLESS)")
print(f"local -> local -> local : relative error {np.nanmean(rt_pair):.3f} "
f"({np.nanmean(rt_pair) / rt_atlas.max():.0e}x larger)")
transitions method = 'auto', has_atlas = True, has_pairwise = True local -> atlas -> local : relative error 1.75e-15 (LOSSLESS) local -> local -> local : relative error 0.138 (8e+13x larger)
b. Going through the atlas is the same switch, in a shared frame¶
A fair question: if the atlas round trip is exact, does routing a chart switch $i\to g\to j$ give something different from the pairwise $i\to j$? No -- and that is the reassuring part. Using $\bm{\Phi}_g\bm{\Phi}_g^\top\mathbf{M}\bm{\Phi}_i=\bm{\Phi}_i$ again,
$$\mathbf{T}_{jg}\mathbf{T}_{gi}=\bm{\Phi}_j^\top\mathbf{M}\bm{\Phi}_g\bm{\Phi}_g^\top\mathbf{M}\bm{\Phi}_i =\bm{\Phi}_j^\top\mathbf{M}\bm{\Phi}_i=\mathbf{T}_{ji},$$
and likewise for the shifts. TransitionMaps.map therefore returns the same numbers
whether method='atlas' or method='pairwise' -- exactness of the atlas does not
change the switch. What it adds is the intermediate object $\bm{z}$: one frame in
which all $K$ charts' states are simultaneously comparable, which is what any
linear combination of states sitting in different charts needs, and in which the
distances of section 2.2 are plain Euclidean norms.
a_probe = A_test[:, :200]
worst = max(float((tm.map(a_probe, ii, jj) - pair.map(a_probe, ii, jj)).abs().max())
for ii in range(K) for jj in range(K))
print(f"max |atlas-routed switch - pairwise switch| over all K^2 pairs = {worst:.2e}")
max |atlas-routed switch - pairwise switch| over all K^2 pairs = 1.87e-14
2.2. Exact distances¶
Chart affiliation is decided by distance to the centroids, so a ROM switches charts using whatever metric it can afford. In atlas coordinates the metric is cheap and exact: because $\bm{q}-\bar{\bm{q}}_g=\bm{\Phi}_g\bm{z}$ and $\bar{\bm{q}}_i-\bar{\bm{q}}_g=\bm{\Phi}_g\bm{d}_{gi}$ with $\bm{\Phi}_g$ $\mathbf{M}$-orthonormal,
$$\left\|\bm{z}-\bm{d}_{gi}\right\|_2=\left\|\bm{q}-\bar{\bm{q}}_i\right\|_{\mathbf{M}} \qquad\text{for every chart }i,$$
i.e. an $r_g$-dimensional Euclidean norm equals the physical distance.
Without the atlas, chart $k$ sees only the in-chart part $\|\bm{d}_{ki}-\bm{a}_k\|$
(tshift[:, k]), which measures the centroid offsets after projecting them onto
chart $k$. On its own, this naive in-chart distance drops whatever separates the
centroids in directions that chart $k$ cannot see, so it under-estimates the distances
and can rank them wrongly. The dropped part depends on the geometry only, not on the
state. For a state $\bm{q}=\bar{\bm{q}}_k+\bm{\Phi}_k\bm{a}_k$ on chart $k$,
$$\left\|\bm{q}-\bar{\bm{q}}_i\right\|_{\mathbf{M}}^2=\left\|\bm{d}_{ki}-\bm{a}_k\right\|_2^2 +\left\|\left(\mathbf{I}-\bm{P}_k\right)\left(\bar{\bm{q}}_i-\bar{\bm{q}}_k\right)\right\|_{\mathbf{M}}^2,$$
in which the second term is the $(i,k)$ entry of oos2 (pairwise_oos2), computed
once, offline, next to tshift. The pairwise metric adds it back, so it is exact too,
and the two strategies rank the centroids identically. What separates them is the round
trip (lossy in section 1.3, lossless in section 2.1a), not the distance.
TransitionMaps.distances implements all three metrics; method='physical' lifts the
state to $\mathbb{R}^{N_h}$ and computes the $\mathbf{M}$-norm there, the $O(N_h K)$
reference that the other two reproduce.
phys = TransitionMaps.from_model(atlas, method="physical")
d_p, d_a, d_x = [], [], []
n_near_a = n_near_x = n_ord_a = n_ord_x = 0
probe = range(0, A_test.shape[1], 5)
for n in probe:
k = int(ids[n])
dp, da, dx = phys.distances(A_test[:, n], k), tm.distances(A_test[:, n], k), pair.distances(A_test[:, n], k)
d_p.append(dp), d_a.append(da), d_x.append(dx)
n_near_a += int(da.argmin() == dp.argmin())
n_near_x += int(dx.argmin() == dp.argmin())
n_ord_a += int(torch.equal(da.argsort(), dp.argsort()))
n_ord_x += int(torch.equal(dx.argsort(), dp.argsort()))
d_p, d_a, d_x = (torch.stack(v).numpy() for v in (d_p, d_a, d_x))
n_probe = len(d_p)
print(f"max |atlas distance - physical distance| = {np.abs(d_a - d_p).max():.2e} (EXACT)")
print(f"max |pairwise distance - physical distance| = {np.abs(d_x - d_p).max():.2e} (in-chart part + oos2)")
print(f"nearest centroid identified correctly : atlas {n_near_a / n_probe:.1%} | pairwise {n_near_x / n_probe:.1%}")
print(f"FULL ordering of all K centroids right: atlas {n_ord_a / n_probe:.1%} | pairwise {n_ord_x / n_probe:.1%}")
max |atlas distance - physical distance| = 1.07e-14 (EXACT) max |pairwise distance - physical distance| = 1.42e-14 (in-chart part + oos2) nearest centroid identified correctly : atlas 100.0% | pairwise 100.0% FULL ordering of all K centroids right: atlas 100.0% | pairwise 100.0%
2.3. When it stops being lossless: truncation¶
The identities of section 2.1a are conditional: they need the atlas to contain $\mathcal{A}_k$ exactly. Truncate $\bm{\Phi}_g$ and they degrade into a projected approximation -- the maps still run, they are just no longer inverses. We keep the same formulas but feed them the first $m\le r_g$ atlas modes.
The cliff is sharp: for $m<r_g$ some chart or centroid lies outside the truncated atlas and the identities break; at $m=r_g$ they return to machine precision. Dimension alone is not the criterion -- inclusion is. (A global POD basis of the same size as one chart, $m=r$, cannot contain $K$ distinct $r$-dimensional flats; the last line of the output measures how far the linear identity fails.)
def maps_from(Phi_g):
"""The section-2.1a formulas for an arbitrary atlas basis (M = I here)."""
Tgk = torch.stack([Phi_g.T @ atlas.Phi_all[:, :, k] for k in range(K)])
Tkg = torch.stack([atlas.Phi_all[:, :, k].T @ Phi_g for k in range(K)])
dgk = torch.stack([Phi_g.T @ (atlas.centroids[k] - atlas.qbar_g) for k in range(K)])
dkg = torch.stack([atlas.Phi_all[:, :, k].T @ (atlas.qbar_g - atlas.centroids[k]) for k in range(K)])
e_lin = max(float((Tkg[k] @ Tgk[k] - I_r).abs().max()) for k in range(K))
e_off = max(float((Tkg[k] @ dgk[k] + dkg[k]).abs().max()) for k in range(K))
return e_lin, e_off
m_sweep = list(range(2, rg + 1, 2))
err_sweep = np.array([maps_from(atlas.Phi_g[:, :m]) for m in m_sweep]) # (m, 2)
for m, (e_lin, e_off) in zip(m_sweep, err_sweep, strict=True):
tag = " <- exact inclusion" if m == rg else ""
print(f" atlas dims m = {m:3d}: |T_kg T_gk - I| = {e_lin:.2e} |T_kg d_gk + d_kg| = {e_off:.2e}{tag}")
Phi_pod_r = compute_pod_basis(Xtrain - atlas.qbar_g[:, None], r)
print(f"global POD basis of size r = {r} (one chart's worth): "
f"|T_kg T_gk - I| = {maps_from(Phi_pod_r)[0]:.2e}")
atlas dims m = 2: |T_kg T_gk - I| = 1.00e+00 |T_kg d_gk + d_kg| = 8.95e+00 atlas dims m = 4: |T_kg T_gk - I| = 1.00e+00 |T_kg d_gk + d_kg| = 4.98e+00 atlas dims m = 6: |T_kg T_gk - I| = 1.00e+00 |T_kg d_gk + d_kg| = 3.74e+00 atlas dims m = 8: |T_kg T_gk - I| = 1.00e+00 |T_kg d_gk + d_kg| = 2.34e+00 atlas dims m = 10: |T_kg T_gk - I| = 2.88e-01 |T_kg d_gk + d_kg| = 9.51e-01 atlas dims m = 12: |T_kg T_gk - I| = 2.64e-01 |T_kg d_gk + d_kg| = 4.57e-01 atlas dims m = 14: |T_kg T_gk - I| = 8.10e-02 |T_kg d_gk + d_kg| = 1.00e-01 atlas dims m = 16: |T_kg T_gk - I| = 1.24e-02 |T_kg d_gk + d_kg| = 1.77e-02 atlas dims m = 18: |T_kg T_gk - I| = 9.59e-03 |T_kg d_gk + d_kg| = 5.25e-03 atlas dims m = 20: |T_kg T_gk - I| = 2.22e-15 |T_kg d_gk + d_kg| = 8.88e-15 <- exact inclusion
global POD basis of size r = 10 (one chart's worth): |T_kg T_gk - I| = 2.12e-01
The builder polices this itself. build_global_assimilation_basis asserts that every
centroid is reproduced by the atlas ($\bar{\bm{q}}_k-\bar{\bm{q}}_g=\bm{\Phi}_g\bm{d}_{gk}$)
and refuses to return a basis that has thrown away real structure, rather than
quietly handing back maps that are not inverses. Its tol is the relative cutoff that
drops the numerically negligible directions of $\bm{U}$. The default, $10^{-11}$, keeps
the structure that the centroid offsets need, and stays above $10^{-12}$, below which
the orthonormality of $\bm{\Phi}_g$ starts to degrade.
for tol in [1e-11, 1e-4, 1e-2]:
try:
Phi_g_t, *_ = build_global_assimilation_basis(atlas.Phi_all, atlas.centroids, atlas.qbar_g, tol=tol)
print(f" tol = {tol:.0e}: r_g = {Phi_g_t.shape[1]}, accepted")
except AssertionError as exc:
print(f" tol = {tol:.0e}: rejected -- {exc}")
tol = 1e-11: r_g = 20, accepted tol = 1e-04: r_g = 20, accepted tol = 1e-02: rejected -- centroids off the atlas (max rel err 1.94e-02); lower tol.
2.4. A non-identity inner product changes nothing¶
Everything above was written with $\mathbf{M}$ in place, and the code carries it
through: pass Mw (a $(N_h,)$ diagonal of quadrature weights, or a dense/sparse FEM
mass matrix) and the local POD becomes $\mathbf{M}$-orthonormal, the transitions use
$\bm{\Phi}_j^\top\mathbf{M}\bm{\Phi}_i$, and the atlas is built with
$\operatorname{orth}_{\mathbf{M}}$. The lossless identities hold verbatim, and atlas
distances now equal the $\mathbf{M}$-weighted physical norm, which is the physically
meaningful norm on a non-uniform mesh.
w = torch.full((atlas.N,), float(fom.L) / atlas.N, dtype=torch.float64) # trapezoidal weights
atlas_w = fit_charts(Xtrain, K=K, r=r, dt=fom.dt, random_state=0, Mw=w)
Phi0 = atlas_w.Phi_all[:, :, 0]
e_lin = max(float((atlas_w.Tkg[k] @ atlas_w.Tgk[k] - I_r).abs().max()) for k in range(K))
e_off = max(float((atlas_w.Tkg[k] @ atlas_w.dgk[k] + atlas_w.dkg[k]).abs().max()) for k in range(K))
print(f"M-orthonormal local basis |Phi^T M Phi - I| = {float((Phi0.T @ (w[:, None] * Phi0) - I_r).abs().max()):.2e}")
print(f"weighted atlas (r_g = {atlas_w.Phi_g.shape[1]}): |T_kg T_gk - I| = {e_lin:.2e}, "
f"|T_kg d_gk + d_kg| = {e_off:.2e}")
apod_w = atlas_w.project_state(Xtest[:, 0])
a_w, k_w = apod_w[:-1, 0], int(apod_w[-1, 0])
z_w = atlas_w.transitions.to_atlas(a_w, k_w)
q_w = atlas_w.centroids[k_w] + atlas_w.Phi_all[:, :, k_w] @ a_w
diff = atlas_w.centroids - q_w[None, :]
d_M = torch.sqrt((diff * (w[None, :] * diff)).sum(dim=1))
print(f"atlas distance == M-weighted physical distance: "
f"{float(((atlas_w.dgk - z_w).norm(dim=1) - d_M).abs().max()):.2e}")
M-orthonormal local basis |Phi^T M Phi - I| = 2.33e-15 weighted atlas (r_g = 20): |T_kg T_gk - I| = 3.55e-15, |T_kg d_gk + d_kg| = 2.44e-15 atlas distance == M-weighted physical distance: 8.88e-16
3. The same identity for operators¶
States are not the only thing that moves between frames. A local quadratic ql-ROM field $\bm{g}_k(\bm{a})=\bm{b}_k+\mathbf{A}_k\bm{a}+\mathsf{B}_k(\bm{a},\bm{a})$ induces one in atlas coordinates, because $\dot{\bm{z}}=\mathbf{T}_{gk}\dot{\bm{a}}_k$ and $\bm{a}_k=\mathbf{T}_{kg}\bm{z}+\bm{d}_{kg}$:
$$\widetilde{\bm{b}}_k^g=\mathbf{T}_{gk}\left[\bm{b}_k+\mathbf{A}_k\bm{d}_{kg} +\mathsf{B}_k(\bm{d}_{kg},\bm{d}_{kg})\right],$$ $$\widetilde{\mathbf{A}}_k^g=\mathbf{T}_{gk}\left[\mathbf{A}_k\mathbf{T}_{kg} +\mathsf{B}_k(\bm{d}_{kg},\mathbf{T}_{kg}\cdot)+\mathsf{B}_k(\mathbf{T}_{kg}\cdot,\bm{d}_{kg})\right], \qquad \widetilde{\mathsf{B}}_k^g(\cdot,\cdot)=\mathbf{T}_{gk}\mathsf{B}_k(\mathbf{T}_{kg}\cdot,\mathbf{T}_{kg}\cdot).$$
Note the shift $\bm{d}_{kg}$ leaks into both the induced constant and the induced linear block. Mapping back is the mirror image, and the composition is again the identity under exactly the two boxed conditions of section 2.1a. This is what lets a correction of $(\bm{b},\mathbf{A})$ made in the shared frame be pushed back into each chart exactly.
g = torch.Generator().manual_seed(3)
b = torch.randn(r, generator=g, dtype=torch.float64)
A = torch.randn(r, r, generator=g, dtype=torch.float64)
B = torch.randn(r, r, r, generator=g, dtype=torch.float64)
k = 2
Tgk, Tkg, dgk, dkg = atlas.Tgk[k], atlas.Tkg[k], atlas.dgk[k], atlas.dkg[k]
b_g = Tgk @ (b + A @ dkg + torch.einsum("ijk,j,k->i", B, dkg, dkg))
A_g = Tgk @ (A @ Tkg + torch.einsum("ijk,j,kl->il", B, dkg, Tkg) + torch.einsum("ijk,jl,k->il", B, Tkg, dkg))
B_g = torch.einsum("ai,ijk,jm,kn->amn", Tgk, B, Tkg, Tkg)
a = torch.randn(r, generator=g, dtype=torch.float64)
z = Tgk @ a + dgk
g_local = b + A @ a + torch.einsum("ijk,j,k->i", B, a, a)
g_global = b_g + A_g @ z + torch.einsum("aij,i,j->a", B_g, z, z)
print(f"induced field agrees: |T_gk g_k(a) - g_g(z)| = {float((Tgk @ g_local - g_global).abs().max()):.2e}")
print(f"operator round trip exact: |T_kg g_g(z) - g_k(a)| = {float((Tkg @ g_global - g_local).abs().max()):.2e}")
induced field agrees: |T_gk g_k(a) - g_g(z)| = 8.88e-14 operator round trip exact: |T_kg g_g(z) - g_k(a)| = 1.19e-13
4. The API in one table¶
| you want | call |
|---|---|
| build charts + atlas from snapshots | fit_charts(X, K, r, dt, Mw=..., build_atlas=True) |
| attach an atlas to existing charts | atlas.build_global_atlas(X) (or set atlas.method = 'auto') |
| one chart's physical round trip | Chart(Phi, centroid, Mw=...).project_state / .recover_state |
| lift to / drop from the shared frame | atlas.transitions.to_atlas(a, k) / .from_atlas(z, k) |
| switch charts | atlas.transitions.map(a, i, j) |
| distance to every centroid | atlas.transitions.distances(a, k) |
| wrap an existing qlROM's charts | Atlas.from_rom(rom) |
| assemble from dolfinx/FEM charts | Atlas.from_charts([...]), Chart.from_fenics(...) |
TransitionMaps.method picks the strategy: 'atlas' (lossless round trip, needs
Tgk/dgk/Tkg/dkg), 'pairwise' (tmap/tshift, plus oos2 for exact distances),
'auto' (atlas when built, else pairwise -- the default), and 'physical', which
changes only distances to the $O(N_h K)$ reference while the maps stay exact. Setting atlas.method = 'auto' on an atlas-less
Atlas builds the shared basis on the spot.
5. Figure¶
# (a) round trips of sections 1.3 and 2.1, (b) distances of section 2.2, (c) truncation of section 2.3
save_figure(plot_summary(rt_pair, rt_atlas, d_p, d_x, d_a, n_ord_x, n_ord_a, m_sweep, err_sweep, rg))
print("SUMMARY:\n"
f" - K={K} charts, r={r} modes, dictionary {K * r + K} cols -> r_g={rg} atlas dims\n"
f" - local->atlas->local {rt_atlas.max():.1e} vs local->local->local {np.nanmean(rt_pair):.2f}\n"
f" - centroid ordering: atlas {n_ord_a / n_probe:.0%} vs pairwise {n_ord_x / n_probe:.0%}\n"
f" - runtime {time.time() - t_start:.0f}s")
SUMMARY: - K=6 charts, r=10 modes, dictionary 66 cols -> r_g=20 atlas dims - local->atlas->local 1.8e-15 vs local->local->local 0.14 - centroid ordering: atlas 100% vs pairwise 100% - runtime 5s
close_pdf()
1 figures written to figs/04_atlas_and_transitions.pdf
PosixPath('figs/04_atlas_and_transitions.pdf')