Tutorial 5: 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, qlOpinf or qlESN 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$.

Contents

  • 1. Charts
    • 1.1. A chart is an affine space
    • 1.2. Chart to chart: the pairwise transition
    • 1.3 ... and why it does not come back
  • 2. The atlas: one affine space that contains all of them
    • 2.1. Transitions in the atlas
      • a. The lossless local -> atlas -> local transition
      • b. Going through the atlas is the same switch, in a shared frame
    • 2.2. Exact distances
    • 2.3. When it stops being lossless: truncation
    • 2.4. A non-identity inner product changes nothing
  • 3. The same identity for operators
  • 4. The API in one table
  • 5. Figure
In [1]:
import time

import matplotlib.pyplot as plt
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

t_start = time.time()
torch.manual_seed(0)

C_ATLAS, C_PAIR, C_PHYS = "#2a78d6", "#eb6834", "#1baf7a"

# every figure below also becomes a page of figs/04_atlas_and_transitions.pdf
start_pdf("04_atlas_and_transitions")
Out[1]:
PosixPath('../outputs/04_atlas_and_transitions.pdf')

1. Charts¶

The 1-D Kuramoto-Sivashinsky equation in its chaotic regime, sampled coarsely (dt = 0.25) since we only need 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. It returns an Atlas, which owns no dynamics at all.

In [2]:
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.

In [3]:
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.

In [4]:
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. On real KS coordinates this is a tens-of-percent effect, not a rounding artefact.

In [5]:
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 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}$ 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 -- the charts share directions, so the atlas is far smaller than the naive column count $Kr+K$.

In [6]:
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$: since $\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.

In [7]:
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.

In [8]:
tm = atlas.transitions          # method='auto' -> the exact atlas, since 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 an ensemble Kalman update needs when its members sit in different charts (see the qlromda data-assimilation tutorial), and what makes the distances in section 2.2 exact.

In [9]:
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 free and exact: since $\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, the fallback is the in-chart surrogate $\|\bm{d}_{ki}-\bm{a}_k\|$ (tshift[:, k]), which measures the centroid offsets after projecting them onto chart $k$ -- it silently drops whatever separates the centroids in directions chart $k$ cannot see, so it under-estimates distances and can rank them wrongly.

TransitionMaps.distances implements all three; method='physical' lifts the state to $\mathbb{R}^{N_h}$ and is the $O(N_h K)$ ground truth (it uses the plain Euclidean norm, so it coincides with the atlas metric exactly when $\mathbf{M}=\mathbf{I}$, as here).

In [10]:
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 |in-chart distance - physical distance| = {np.abs(d_x - d_p).max():.3f}   (biased low)")
print(f"nearest centroid identified correctly : atlas {n_near_a / n_probe:.1%} | in-chart {n_near_x / n_probe:.1%}")
print(f"FULL ordering of all K centroids right: atlas {n_ord_a / n_probe:.1%} | in-chart {n_ord_x / n_probe:.1%}")

max |atlas distance - physical distance| = 1.07e-14   (EXACT)
max |in-chart distance - physical distance| = 0.545   (biased low)
nearest centroid identified correctly : atlas 100.0% | in-chart 99.0%
FULL ordering of all K centroids right: atlas 100.0% | in-chart 83.5%

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 smoothly 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: the identity is broken for every $m<r_g$ and machine-exact at $m=r_g$. Dimension alone is not the criterion -- inclusion is. (A global POD basis of the same size as one chart, $m=r$, is off by ~20%: it cannot contain six $r$-dimensional flats.)

In [11]:
def maps_from(Phi_g):
    """The section-5 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 ($\bm{c}_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 eigenvalue cutoff on $\bm{\Lambda}$; the 1e-11 default is deliberately far below the $10^{-12}$ regime where $\bm{\Phi}_g$'s own orthonormality would start to suffer.

In [12]:
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 -- the physically meaningful one on a non-uniform mesh.

In [13]:
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 data-assimilation step correct $(\bm{b},\mathbf{A})$ in the shared frame and push the correction back into each chart (qlroms.model.QLModel, est_theta).

In [14]:
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' (exact, needs Tgk/dgk/Tkg/dkg), 'pairwise' (tmap/tshift only), 'auto' (atlas when built, else pairwise -- the default), and 'physical', which changes only distances to the $O(N_h K)$ ground truth while the maps stay exact. Setting atlas.method = 'auto' on an atlas-less Atlas builds the shared basis on the spot.

5. Figure¶

In [15]:
fig, axs = plt.subplot_mosaic([["a", "b", "c"]], figsize=(14, 4.2), constrained_layout=True)

# (a) round-trip error: pairwise vs atlas, per chart
xs = np.arange(K)
axs["a"].bar(xs - 0.2, np.nanmean(rt_pair, axis=1), width=0.4, color=C_PAIR, label=r"pairwise $k\to j\to k$")
axs["a"].bar(xs + 0.2, np.maximum(rt_atlas, 1e-17), width=0.4, color=C_ATLAS, label=r"atlas $k\to g\to k$")
axs["a"].axhline(np.finfo(float).eps, color="0.4", lw=1.0, ls=":")
axs["a"].text(-0.45, 3.5e-16, "machine eps", fontsize=8, color="0.4", ha="left")
axs["a"].set(yscale="log", xlabel="chart $k$", ylabel="relative round-trip error",
             xticks=xs, ylim=(1e-17, 1e3), title="(a) Round trips: lossy vs lossless")
axs["a"].legend(loc="upper center", frameon=False, fontsize=9, ncols=2)

# (b) how wrong is the metric each strategy switches on?
axs["b"].axhline(0.0, color="0.4", lw=1.0, ls="--", zorder=0)
axs["b"].scatter(d_p.ravel(), (d_x - d_p).ravel(), s=6, alpha=0.35, color=C_PAIR,
                 label=f"in-chart ({n_ord_x / n_probe:.0%} of orderings right)")
axs["b"].scatter(d_p.ravel(), (d_a - d_p).ravel(), s=6, alpha=0.35, color=C_ATLAS,
                 label=f"atlas ({n_ord_a / n_probe:.0%} of orderings right)")
axs["b"].set(xlabel=r"physical distance $\|q - \bar q_i\|$",
             ylabel="metric $-$ physical distance",
             title="(b) Chart-switching metric vs truth")
axs["b"].legend(loc="lower left", frameon=False, fontsize=9)

# (c) truncation sweep
axs["c"].semilogy(m_sweep, np.maximum(err_sweep[:, 0], 1e-17), "o-", color=C_ATLAS,
                  label=r"$\|T_{kg}T_{gk} - I\|_\infty$")
axs["c"].semilogy(m_sweep, np.maximum(err_sweep[:, 1], 1e-17), "s-", color=C_PHYS,
                  label=r"$\|T_{kg}d_{gk} + d_{kg}\|_\infty$")
axs["c"].axvline(rg, color="0.4", lw=1.0, ls="--")
axs["c"].text(rg - 0.4, 1e-8, rf"$r_g = {rg}$: exact inclusion", rotation=90, fontsize=9, ha="right")
axs["c"].set(xlabel="atlas dimension $m$ retained", ylabel="identity residual",
             ylim=(1e-17, 3e1), title="(c) Losslessness needs exact inclusion")
axs["c"].legend(loc="lower left", frameon=False, fontsize=9)

for ax in axs.values():
    ax.grid(alpha=0.25, lw=0.5)

# outputs land next to this tutorial whether it runs from the repo root or from tutorials/
save_figure(fig)
No description has been provided for this image
In [16]:
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 in-chart {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 in-chart 84%
  - runtime 4s
In [17]:
close_pdf()

1 figures written to ../outputs/04_atlas_and_transitions.pdf
Out[17]:
PosixPath('../outputs/04_atlas_and_transitions.pdf')