gam-sae 0.3.154

//! ═══════════════════════════════════════════════════════════════════════════
//! Structure-theory grounding (memo Part VI, "binding = transport", Prop H):
//!
//! A fitted atom's decoder `γ = Φ(t)·B` and its latent chart `t` are related by
//! an intrinsic (gauge-invariant) smoothness penalty that scores every
//! reparameterization `t ↦ φ(t)`, `φ ∈ Diff(M)`, identically — the penalty sees
//! only the pullback metric / curve up to diffeomorphism, never the coordinate
//! itself. So the fit does not pick a point in coordinate space; it picks an
//! entire `Diff(M)` ORBIT, and that orbit IS a fiber of the memo's GAUGE
//! GROUPOID: two charts `t`, `φ(t)` on the same orbit are gauge-equivalent
//! representations of one physical atom, and reconstruction EV cannot see the
//! difference (the planted-circle witness: a full loop compressed into ~1 rad
//! of chart span still scores EV 0.9979 — image perfect, chart arbitrary).
//!
//! This module is a SLICE of that groupoid: for each reference topology it
//! picks out ONE canonical representative per orbit —
//!   * `d = 1` (circle / interval): the arc-length (unit-speed) representative
//!     ([`unit_speed_reparameterization`]);
//!   * `d = 2` torus / free patch: the minimum-isometry-defect flow chart
//!     ([`torus_isometry_flow_reparameterization`],
//!     [`patch_isometry_flow_reparameterization`]);
//!   * `d = 2` sphere: the round-sphere conformal-boost chart
//!     ([`sphere_isometry_flow_reparameterization`]).
//! Each slice is EXACT and IMAGE-FROZEN: the recomposed decoder reproduces the
//! original decoded curve/surface pointwise (`recompose_decoder_exact_ls`) and
//! the transport is a congruence that preserves the smoothness quadratic form
//! (`B̃ᵀ S̃ B̃ = Bᵀ S B`). So picking the slice changes nothing the objective (data
//! fit + intrinsic smoothness) can see — it is a gauge RETRACTION along the
//! orbit, the concrete realization of the memo's "binding = transport" slogan:
//! moving the *representation* to a canonical gauge without moving the physics.
//!
//! After the slice is taken, what remains of the chart freedom is no longer the
//! full continuous `Diff(M)` — it collapses to the finite LINEAR STABILIZER of
//! the reference manifold under its own isometry group: `O(2)` (rotation +
//! reflection) for the circle, reflection + translation for the interval,
//! `Isom(T², flat) = U(1)² ⋊ D₄` for the torus, `O(2) ⋉ ℝ²` for the flat patch,
//! `O(3)` for the round sphere. THIS residual group is exactly Proposition H's
//! uniqueness condition for canonical layer transport: two fitted atoms (or two
//! layers of the same atom across refits) are honestly bindable/comparable iff
//! they agree up to this residual stabilizer — nothing coarser is guaranteed,
//! and nothing finer is left to compare. The `PinnedByCanonicalization`
//! provenance ([`crate::identifiability::VerdictProvenance`]) is the
//! certificate that this collapse — continuous `Diff(M)` freedom sliced down to
//! the finite residual — has actually been performed for a given atom.
//!
//! The honesty gate threaded through every arm (`CHART_RECOMPOSITION_REL_TOL`,
//! `recompose_decoder_exact_ls` returning `Ok(None)`) is the module's
//! realization of the memo's honesty doctrine: a basis too poor to carry the
//! reparameterized image through exactly REFUSES the slice rather than taking
//! it lossily — canonicalization is a retraction or it does not happen. The
//! diffeomorphism guards (`SAE_FLOW_DIFFEO_MIN_DET`,
//! `certified_min_jacobian_det`) are the complementary safety rail: they keep
//! the candidate slice INSIDE `Diff(M)`, because a folded map leaves the gauge
//! orbit entirely rather than moving along it.
//!
//! ⚠ `det Dφ > 0` everywhere does NOT by itself imply the map is injective — a
//! polynomial map can have a positive Jacobian determinant everywhere and still
//! be many-to-one (Pinchuk 1994). Each arm supplies the extra hypothesis that
//! actually closes the inference, and it is the STRUCTURE of the flow family,
//! never the guard, that closes it:
//!
//!  * torus — the flow is a degree-`(1,1)` self-map of `T²` (the displacement is
//!    periodic), and a degree-`(1,1)` self-map with `det Dφ > 0` is a global
//!    diffeomorphism;
//!  * free patch — `PATCH_FLOW_MAX_DEGREE = 1`, so the flow is AFFINE and
//!    nonzero determinant gives global invertibility by linear algebra;
//!  * sphere — the conformal Möbius boosts are diffeomorphisms by construction;
//!  * `d = 1` arc length — monotone by construction, hence degree 1.
//!
//! So the guard's job is narrow and should be read narrowly: it refuses folds
//! (`det Dφ ≤ δ`) within a family whose injectivity is already established by
//! degree. Widening any family's degree invalidates the injectivity argument and
//! requires a new one — no amount of determinant checking substitutes for it.
//! ═══════════════════════════════════════════════════════════════════════════
//!
//! #1019 stage 1 / #2022 — arc-length (unit-speed) chart canonicalization for `d = 1`.
//! The arc-length reparameterization is exact and IMAGE-FROZEN (reconstruction
//! unchanged; the transport preserves `BᵀSB`), so it is a gauge RETRACTION along
//! the `Diff(M)` orbit that leaves the data-fit AND intrinsic-smoothness objective
//! invariant. #2022 promotes it from a post-fit pass to IN-LOOP enforcement (see
//! [`unit_speed_retraction`]), applied at chart-refresh boundaries — never inside a
//! line search, because it re-gauges `t` and thus changes the ARD *coordinate*
//! prior energy (the term that pins the residual gauge to `t → ±t + c`). Post-fit
//! canonicalization then reduces to a verification no-op.
//! #1019 stage 1 — arc-length (unit-speed) chart canonicalization for `d = 1`
//! manifold atoms (circle and interval topologies).
//!
//! A fitted atom's latent chart is gauge-arbitrary: the gauge-invariant
//! (intrinsic) smoothness penalty makes every reparameterization of the latent
//! coordinate equal-cost BY DESIGN, so nothing in the likelihood prefers the
//! honest chart and reconstruction metrics cannot detect chart dishonesty
//! (the planted-circle failure that motivated #1019 compressed the full loop
//! into ~1 rad of coordinate span at reconstruction EV 0.9979 — image
//! perfect, chart arbitrary).
//!
//! For `d = 1` the canonical representative of the `Diff(S¹)` /
//! `Diff([0, 1])` orbit is mathematically solved: the **arc-length
//! reparameterization**. Post-fit and image-frozen, compute the cumulative
//! arc length `s(t) = ∫_lo^t ‖γ'(u)‖ du` along the fitted decoder curve
//! `γ(t) = Φ(t) B`, normalize by the total length `L` to the chart's native
//! span (the basis period for a circle, the unit interval for a line
//! segment), and reparameterize: new coordinates `t̃_i = s(t_i)`, new decoder
//! `B̃` refit by exact least squares of the ORIGINAL decoded curve on a fine
//! grid against the basis at the new coordinates. The refit is linear and
//! exact up to basis expressiveness; the recomposition residual is recorded
//! and the canonicalization is REFUSED when it exceeds a small tolerance
//! relative to the curve scale — an honest fallback, never a lossy silent
//! swap.
//!
//! After canonicalization the atom's residual chart freedom downgrades from
//! the full diffeomorphism group to the finite isometry group of the
//! reference manifold: rotation + reflection (`O(2)`) for the circle,
//! reflection + translation for the interval. The certificate records this
//! with the `PinnedByCanonicalization` provenance
//! ([`crate::identifiability::VerdictProvenance`]).
//!
//! #1019 stage 2 (`d = 2`, torus): the canonical representative of the
//! `Diff(T²)` orbit is the **minimum-isometry-defect flow** chart. The chart
//! map is parameterized as `φ_θ(t) = t + Σ_k θ_k v_k(t)` with `v_k` a fixed
//! truncated Fourier vector-field basis on `T²` (orders ≤ 2 per axis, both
//! components — a few tens of coefficients, wrap-around respected by
//! construction), and `θ` minimizes the discretized isometry defect over the
//! fitted rows with an exact analytic Gauss–Newton (see
//! [`torus_isometry_flow_reparameterization`] for the full derivation). A
//! hard diffeomorphism guard `det Dφ_θ > δ` on a check grid means a folded
//! chart is REFUSED, never produced. The decoder transport is the same
//! exact-LS recomposition — and the same honesty gate — as the `d = 1` path
//! (shared helper `recompose_decoder_exact_ls`).
//!
//! #1019 free-chart arm (`d = 2`, free/patch): a contractible Euclidean-patch
//! atom (Duchon / EuclideanPatch basis kind) admits a **global** truncated
//! flow basis — polynomial vector fields `v_{c,(a,b)}(t) = e_c · u₀^a u₁^b`
//! on the normalized patch box ([`FreePatchFlowBasis`]) — with no hairy-ball
//! obstruction (a contractible patch carries nowhere-vanishing vector fields).
//! So unlike the sphere it is genuinely MINIMIZED, not merely measured: the
//! isometry defect is descended against the flat reference `g_ref = I`
//! ([`patch_isometry_flow_reparameterization`]), pinning the uniform-speed
//! (minimum-anisotropy) chart, with the residual chart freedom downgraded to
//! the flat isometry group `O(2) ⋉ ℝ²` and provenance
//! `PinnedByCanonicalization`. The torus, free-patch, and sphere arms share one
//! pullback-metric extraction (`extract_pullback_metric_d2`) and the two
//! flow-pinned arms share one exact Gauss–Newton core
//! (`minimize_isometry_defect_flow`).
//!
//! `S²` (sphere atoms): the hairy-ball theorem rules out a single global
//! pole-free flow basis the way the torus and free-patch paths use, so the
//! sphere representative is the **round-sphere conformal-boost flow** — the
//! gradients of the three degree-1 harmonics (`K_z` zonal, pole-free in its one
//! latitude component; `K_x`, `K_y` carrying the longitudinal pole the theorem
//! forces somewhere). Minimizing the isometry defect over these three boosts
//! breaks the conformal (Möbius) moduli down to the round sphere's isometry
//! group `O(3)` — the chart pathology a bare harmonic energy cannot pin
//! ([`sphere_isometry_flow_reparameterization`], on a pole-margin band so the
//! `1/cos lat` boost stays well-conditioned). The same exact-image-frozen LS
//! decoder transport gates the commit: a boosted image only freezes to the
//! recomposition floor inside a basis rich enough to absorb it, so a too-poor
//! basis honestly refuses rather than silently altering the image. The
//! round-sphere isometry DEFECT remains available as a standalone measurement
//! (`sphere_chart_isometry_defect`).

use faer::Side as FaerSide;
use ndarray::{Array1, Array2, ArrayView1, ArrayView2};

use crate::manifold::{SaeBasisEvaluator, solve_design_least_squares};
use gam_linalg::faer_ndarray::{FaerCholesky, fast_ab, fast_ata, fast_atb};

/// Number of integration cells for the fitted-turning quadrature (#1026).
const TURNING_QUADRATURE_CELLS: usize = 256;

/// Number of grid CELLS for the arc-length quadrature and the decoder
/// recomposition least squares. Each cell carries a node, a midpoint, and the
/// next node (composite Simpson), so the speed field is sampled at
/// `2·ARC_LENGTH_GRID_CELLS + 1` points; the per-cell quadrature error is
/// `O(Δu⁴)`, far below the recomposition tolerance.
pub const ARC_LENGTH_GRID_CELLS: usize = 2048;

/// Relative image-recomposition tolerance: the canonicalization is refused
/// (honest fallback to the fitted chart) when the max-abs difference between
/// the original decoded curve and its recomposition through the new chart
/// exceeds this fraction of the curve scale — on the audit grid OR on the
/// fitted rows. Matched to the image-invariance contract the certificate and
/// the tests assert (reconstruction unchanged within 1e-8).
pub const CHART_RECOMPOSITION_REL_TOL: f64 = 1.0e-9;

/// The `d = 1` reference topology the canonical chart lives on.
#[derive(Debug, Clone, PartialEq)]
pub enum CanonicalChartTopology {
    /// `S¹` with the basis evaluator's native period (`1.0` for the
    /// fraction-of-period harmonic evaluators). Arc length is accumulated
    /// over one full period and rescaled so the canonical chart spans
    /// exactly one period; the residual chart freedom is `O(2)` (base-point
    /// rotation + orientation reflection).
    Circle { period: f64 },
    /// A line-segment chart. Arc length is accumulated over the fitted
    /// coordinate range and rescaled to `[0, 1]`; the residual chart freedom
    /// is reflection + translation.
    Interval,
}

/// The exact, image-frozen arc-length reparameterization of one `d = 1` atom.
#[derive(Debug, Clone)]
pub struct UnitSpeedReparameterization {
    /// Canonical per-row coordinates `t̃_i = span · s(t_i) / L`, length `n`.
    pub new_row_coords: Array1<f64>,
    /// Recomposed decoder coefficients `B̃ = T · B`, shape `(m, p)`: the exact
    /// least-squares refit of the original decoded curve (on the audit grid)
    /// against the basis at the canonical coordinates.
    pub new_decoder: Array2<f64>,
    /// The `(m, m)` basis transport `T` with `Φ(t̃) · T ≈ Φ(t)` on the grid —
    /// the same object the affine gauge canonicalization uses to transport
    /// the smoothness Gram (`S̃ = T⁻ᵀ S T⁻¹` preserves `B̃ᵀ S̃ B̃ = Bᵀ S B`).
    pub decoder_transport: Array2<f64>,
    /// Total arc length `L` of the decoder curve over the canonical domain.
    pub total_arc_length: f64,
    /// Max-abs recomposition error on the audit grid, relative to the curve
    /// scale. Always `≤ CHART_RECOMPOSITION_REL_TOL` when `Some(..)` is
    /// returned.
    pub recomposition_residual: f64,
}

/// Decoder-curve speed `‖Φ'(u) B‖₂` for each evaluated coordinate row, from
/// the basis jet `(rows, m, 1)` and the decoder `(m, p)`.
pub(crate) fn curve_speeds(
    jet: &ndarray::Array3<f64>,
    decoder: ArrayView2<'_, f64>,
) -> Result<Vec<f64>, String> {
    let (rows, m, d) = jet.dim();
    if d != 1 {
        return Err(format!(
            "sae_chart_canonicalization: expected a 1-D latent jet, got latent_dim {d}"
        ));
    }
    if decoder.nrows() != m {
        return Err(format!(
            "sae_chart_canonicalization: jet basis width {m} != decoder rows {}",
            decoder.nrows()
        ));
    }
    let p = decoder.ncols();
    let mut speeds = Vec::with_capacity(rows);
    let mut tangent = vec![0.0_f64; p];
    for row in 0..rows {
        for slot in tangent.iter_mut() {
            *slot = 0.0;
        }
        for bm in 0..m {
            let dphi = jet[[row, bm, 0]];
            if dphi == 0.0 {
                continue;
            }
            for (j, slot) in tangent.iter_mut().enumerate() {
                *slot += dphi * decoder[[bm, j]];
            }
        }
        speeds.push(tangent.iter().map(|v| v * v).sum::<f64>().sqrt());
    }
    Ok(speeds)
}

/// Exact integral of the cell-local quadratic speed interpolant (through the
/// node, midpoint, and next-node speeds) over `[0, x]`, `x ∈ [0, h]`. At
/// `x = h` this is exactly the Simpson cell weight `h(f0 + 4fm + f1)/6`.
fn partial_cell_arc(f0: f64, fm: f64, f1: f64, h: f64, x: f64) -> f64 {
    if h <= 0.0 {
        return 0.0;
    }
    let a = (2.0 * f0 - 4.0 * fm + 2.0 * f1) / (h * h);
    let b = (-3.0 * f0 + 4.0 * fm - f1) / h;
    let x2 = x * x;
    a * x2 * x / 3.0 + b * x2 / 2.0 + f0 * x
}

/// Compute the arc-length (unit-speed) reparameterization of a fitted `d = 1`
/// atom: the canonical per-row coordinates and the exactly-recomposed decoder.
///
/// Image-frozen: the decoder curve is never refit against data — only
/// re-expressed in the canonical chart. Returns `Ok(None)` (honest skip,
/// never a lossy swap) when:
/// * the chart is degenerate (no rows, empty basis, zero/non-finite total
///   arc length, collapsed interval range), or
/// * the basis family cannot absorb the reparameterized curve within
///   [`CHART_RECOMPOSITION_REL_TOL`] of the curve scale on the audit grid.
// Gauge-theory reading: `row_coords` ranges over one point of the `Diff(S¹)` /
// `Diff([0,1])` orbit the fit actually picked; this function computes the
// SLICE representative of that same orbit (arc length is the unique
// unit-speed member) and the transport `T` that carries the decoder there.
// Everything downstream of a successful `Some(..)` return is gauge-equivalent
// to the input under the residual `O(2)` / reflection+translation stabilizer
// — never a different physical fit.
pub fn unit_speed_reparameterization(
    evaluator: &dyn SaeBasisEvaluator,
    decoder: ArrayView2<'_, f64>,
    row_coords: ArrayView1<'_, f64>,
    topology: &CanonicalChartTopology,
) -> Result<Option<UnitSpeedReparameterization>, String> {
    let n = row_coords.len();
    let m = decoder.nrows();
    let p = decoder.ncols();
    if n == 0 || m == 0 || p == 0 {
        return Ok(None);
    }
    for &t in row_coords.iter() {
        if !t.is_finite() {
            return Ok(None);
        }
    }

    // ── Canonical quadrature domain `[lo, hi]` and target span ──────────────
    let (lo, hi, span) = match topology {
        CanonicalChartTopology::Circle { period } => {
            if !(period.is_finite() && *period > 0.0) {
                return Err(format!(
                    "sae_chart_canonicalization: circle period must be finite and positive; got {period}"
                ));
            }
            // The decoder curve is defined over the WHOLE period (the chart
            // dishonesty being cured is exactly rows compressed into a sliver
            // of it), so the arc length integrates the full circle.
            (0.0, *period, *period)
        }
        CanonicalChartTopology::Interval => {
            let mut t_min = f64::INFINITY;
            let mut t_max = f64::NEG_INFINITY;
            for &t in row_coords.iter() {
                t_min = t_min.min(t);
                t_max = t_max.max(t);
            }
            let scale = t_min.abs().max(t_max.abs()).max(1.0);
            if !(t_max - t_min > 1.0e-12 * scale) {
                // Collapsed chart: every row at one point — arc length cannot
                // define a chart there.
                return Ok(None);
            }
            (t_min, t_max, 1.0)
        }
    };

    // ── Speed field on the Simpson grid (nodes + midpoints in one call) ─────
    let cells = ARC_LENGTH_GRID_CELLS;
    let h = (hi - lo) / cells as f64;
    let mut quad_coords = Array2::<f64>::zeros((2 * cells + 1, 1));
    for j in 0..=cells {
        quad_coords[[2 * j, 0]] = lo + j as f64 * h;
        if j < cells {
            quad_coords[[2 * j + 1, 0]] = lo + (j as f64 + 0.5) * h;
        }
    }
    let (grid_phi_all, grid_jet_all) = evaluator.evaluate(quad_coords.view())?;
    if grid_phi_all.ncols() != m {
        return Err(format!(
            "sae_chart_canonicalization: evaluator basis width {} != decoder rows {m}",
            grid_phi_all.ncols()
        ));
    }
    let speeds = curve_speeds(&grid_jet_all, decoder)?;
    if speeds.iter().any(|s| !s.is_finite()) {
        return Ok(None);
    }

    // ── Composite-Simpson cumulative arc length at the nodes ────────────────
    let mut cumulative = vec![0.0_f64; cells + 1];
    for j in 0..cells {
        let f0 = speeds[2 * j];
        let fm = speeds[2 * j + 1];
        let f1 = speeds[2 * j + 2];
        cumulative[j + 1] = cumulative[j] + h * (f0 + 4.0 * fm + f1) / 6.0;
    }
    let total = cumulative[cells];
    if !(total.is_finite() && total > 0.0) {
        return Ok(None);
    }
    let rescale = span / total;

    // ── The canonical chart map `t ↦ span · s(t) / L` ───────────────────────
    let map_coord = |t: f64| -> f64 {
        let local = match topology {
            CanonicalChartTopology::Circle { period } => (t - lo).rem_euclid(*period),
            CanonicalChartTopology::Interval => (t - lo).clamp(0.0, hi - lo),
        };
        let cell = ((local / h).floor() as usize).min(cells - 1);
        let x = local - cell as f64 * h;
        let s = cumulative[cell]
            + partial_cell_arc(
                speeds[2 * cell],
                speeds[2 * cell + 1],
                speeds[2 * cell + 2],
                h,
                x,
            );
        let mapped = rescale * s;
        match topology {
            CanonicalChartTopology::Circle { period } => mapped.rem_euclid(*period),
            CanonicalChartTopology::Interval => mapped.clamp(0.0, span),
        }
    };

    let new_row_coords = Array1::from_iter(row_coords.iter().map(|&t| map_coord(t)));

    // ── Decoder recomposition: exact LS of the original curve on the grid ───
    // Audit grid = the quadrature nodes. `Φ_new · T ≈ Φ_old` (row j of Φ_old
    // is the basis at node u_j, row j of Φ_new is the basis at the node's
    // canonical image s̃(u_j)), so `B̃ = T · B` reproduces the original curve
    // values `Φ_old · B` at the canonical coordinates — the image is frozen.
    let mut node_new_coords = Array2::<f64>::zeros((cells + 1, 1));
    let mut old_phi = Array2::<f64>::zeros((cells + 1, m));
    for j in 0..=cells {
        node_new_coords[[j, 0]] = map_coord(lo + j as f64 * h);
        for bm in 0..m {
            old_phi[[j, bm]] = grid_phi_all[[2 * j, bm]];
        }
    }
    let Some(recomposition) =
        recompose_decoder_exact_ls(evaluator, decoder, old_phi.view(), node_new_coords.view())?
    else {
        // The basis family is not expressive enough to carry the
        // arc-length-reparameterized curve: refuse rather than swap lossily.
        return Ok(None);
    };

    Ok(Some(UnitSpeedReparameterization {
        new_row_coords,
        new_decoder: recomposition.new_decoder,
        decoder_transport: recomposition.transport,
        total_arc_length: total,
        recomposition_residual: recomposition.recomposition_residual,
    }))
}

/// Exact-LS decoder recomposition shared by the `d = 1` (arc-length) and
/// `d = 2` (torus isometry-flow) canonicalizations.
///
/// `old_phi` is the ORIGINAL basis at the audit grid (so `old_phi · B` is the
/// original decoded image there) and `new_coords` are the same grid points'
/// canonical images. Solves the basis transport `Φ(new) · T ≈ Φ(old)` by
/// exact least squares, recomposes `B̃ = T · B` (so the decoded image is
/// reproduced at the transported coordinates — image-frozen), and applies the
/// honesty gate: returns `Ok(None)` (refuse, never a lossy silent swap) when
/// the max-abs image drift on the audit grid exceeds
/// [`CHART_RECOMPOSITION_REL_TOL`] of the image scale.
pub(crate) struct DecoderRecomposition {
    /// `(m, m)` basis transport `T` with `Φ(new) · T ≈ Φ(old)` on the grid.
    pub transport: Array2<f64>,
    /// Recomposed decoder `B̃ = T · B`, shape `(m, p)`.
    pub new_decoder: Array2<f64>,
    /// Max-abs recomposition error on the audit grid, relative to the image
    /// scale. Always `≤ CHART_RECOMPOSITION_REL_TOL` when returned.
    pub recomposition_residual: f64,
}

// This is the module's one TRANSPORT primitive in the "binding = transport"
// sense: `T` is the exact congruence carrying the decoder from the fitted
// gauge to the canonical-slice gauge with the decoded image held fixed
// pointwise (image-frozen — no new information enters, only the
// representation moves). The honesty gate below is what makes the slice
// trustworthy: if the basis cannot carry the reparameterized curve/surface to
// within `CHART_RECOMPOSITION_REL_TOL`, the function REFUSES (`Ok(None)`)
// instead of committing a lossy transport — the slice is taken only when it
// is exact, never approximately.
pub(crate) fn recompose_decoder_exact_ls(
    evaluator: &dyn SaeBasisEvaluator,
    decoder: ArrayView2<'_, f64>,
    old_phi: ArrayView2<'_, f64>,
    new_coords: ArrayView2<'_, f64>,
) -> Result<Option<DecoderRecomposition>, String> {
    let m = decoder.nrows();
    let (new_phi, new_jet) = evaluator.evaluate(new_coords)?;
    if new_phi.ncols() != m
        || new_phi.nrows() != old_phi.nrows()
        || new_jet.dim() != (new_coords.nrows(), m, new_coords.ncols())
    {
        return Err(format!(
            "sae_chart_canonicalization: evaluator returned basis {:?} / jet {:?} at the canonical grid; expected ({}, {m}) with latent_dim {}",
            new_phi.dim(),
            new_jet.dim(),
            old_phi.nrows(),
            new_coords.ncols()
        ));
    }
    let transport = solve_design_least_squares(new_phi.view(), old_phi)?;
    let new_decoder = fast_ab(&transport, &decoder);

    // ── Honest gate: max-abs recomposition error relative to image scale ────
    let old_fit = fast_ab(&old_phi, &decoder);
    let new_fit = fast_ab(&new_phi, &new_decoder);
    let mut fit_scale = 0.0_f64;
    let mut max_abs = 0.0_f64;
    for (a, b) in old_fit.iter().zip(new_fit.iter()) {
        fit_scale = fit_scale.max(a.abs()).max(b.abs());
        max_abs = max_abs.max((a - b).abs());
    }
    if !(fit_scale.is_finite() && fit_scale > 0.0 && max_abs.is_finite()) {
        return Ok(None);
    }
    let recomposition_residual = max_abs / fit_scale;
    if recomposition_residual > CHART_RECOMPOSITION_REL_TOL {
        // The basis family cannot absorb the reparameterized image: refuse.
        return Ok(None);
    }
    Ok(Some(DecoderRecomposition {
        transport,
        new_decoder,
        recomposition_residual,
    }))
}

// ════════════════════════════════════════════════════════════════════════════
// #1019 stage 2 — d = 2 torus isometry-flow chart canonicalization
//
// Gauge-theory reading: the `Diff(T²)` orbit here is infinite-dimensional; this
// arm slices it via a finite-dimensional PARAMETERIZED sub-family (the
// truncated Fourier flow `φ_θ`) and descends the isometry defect over `θ` to
// find the flat-isometric representative within that family. The slice is
// still exact (image-frozen recomposition, honesty-gated) even though the
// search over `θ` is a numerical minimization rather than a closed form — the
// residual gauge left over is `Isom(T², flat) = U(1)² ⋊ D₄` (pure torus
// translations, the axis-swap/reflection dihedral symmetry), which is why
// translations are excluded from the flow basis below: they are already in the
// residual stabilizer, not something the slice needs to remove.
// ════════════════════════════════════════════════════════════════════════════

/// Highest Fourier order per axis of the truncated flow basis on `T²`:
/// frequency vectors `(a, b)` with `|a|, |b| ≤ 2`. One representative per
/// antipodal pair `±(a, b)` (12 of them), `sin` + `cos` phases, both vector
/// components ⇒ `2 · 2 · 12 = 48` flow coefficients — the "dim ~ tens"
/// unconstrained smooth problem of #1019 stage 2. Constants (the pure torus
/// translations) are EXCLUDED on purpose: translations are exact isometries
/// of `(T², g_ref)`, so they leave the defect invariant and would only insert
/// null directions into the Gauss–Newton system.
pub const TORUS_FLOW_MAX_HARMONIC: i32 = 2;

/// Diffeomorphism floor `δ`: a candidate flow is REJECTED (the line search
/// treats it as a failed step; the final chart is never produced) when
/// `det Dφ_θ ≤ δ` anywhere on the check grid. `θ = 0` has `det Dφ = 1`
/// everywhere and only guarded steps are ever accepted, so the optimizer can
/// never walk through a fold. The floor is one-directional: raising it only
/// ever rejects MORE near-fold candidates, never accepts a fold, so it is a
/// conservative guard rather than a tuned operating point. The torus, free
/// patch, and sphere flows all enforce the IDENTICAL contract, so the value
/// lives here once and the per-topology floors below alias it.
///
/// Gauge-theory role: `det Dφ > 0` is precisely the condition that `φ` stays
/// inside `Diff(M)` — a fold (`det Dφ ≤ 0` somewhere) is not a reparameteriza-
/// tion at all, it is a map that identifies distinct points or reverses
/// orientation locally, so it would leave the gauge orbit rather than move
/// along it. This constant is what keeps every flow-based slice a genuine
/// gauge retraction instead of an information-destroying collapse.
pub const SAE_FLOW_DIFFEO_MIN_DET: f64 = 0.1;

/// Torus-flow diffeomorphism floor — see [`SAE_FLOW_DIFFEO_MIN_DET`].
pub const TORUS_FLOW_DIFFEO_MIN_DET: f64 = SAE_FLOW_DIFFEO_MIN_DET;

/// Per-axis node count of the diffeomorphism-guard check grid. The flow basis
/// is band-limited to `TORUS_FLOW_MAX_HARMONIC` (≤ 2 oscillations per axis), so
/// `det Dφ_θ` is band-limited to ≤ 4 per axis and 64 nodes oversample it by 16×.
///
/// Oversampling is not a proof. "Band-limited and sampled densely" bounds
/// nothing about the value BETWEEN nodes without a derivative bound, and the
/// guard's contract is `det Dφ_θ > δ` *everywhere*. The node count therefore
/// no longer carries the argument on its own:
/// [`TorusFlowBasis::certified_min_jacobian_det`] subtracts a closed-form
/// inter-node correction from this grid's minimum and is what the gate consumes.
/// `N` now sets only how TIGHT that certificate is — the correction scales as
/// `1/N` — never whether it is valid.
pub const TORUS_FLOW_GUARD_NODES_PER_AXIS: usize = 64;

/// Outer iteration cap for the damped Gauss–Newton flow optimization. The
/// problem is a 48-dimensional smooth nonlinear least squares; quadratic
/// local convergence makes this cap generous (termination is normally by the
/// relative step / improvement tolerances below).
pub const TORUS_FLOW_GN_MAX_ITERS: usize = 80;

/// Consecutive damping escalations before the Gauss–Newton declares the
/// current iterate a local minimum and stops.
pub const TORUS_FLOW_GN_MAX_REJECTS: usize = 12;

/// Levenberg damping growth factor applied on a rejected trial step (a Cholesky
/// factorization failure, or a folded / non-improving candidate) inside the
/// shared `d = 2` flow trust loop [`lm_damped_accept_sweep`]: `λ ←
/// LM_LAMBDA_GROWTH · λ` tightens the trust region toward the (safe) gradient
/// direction. Shared verbatim by the torus and sphere-boost Gauss–Newton cores.
const LM_LAMBDA_GROWTH: f64 = 10.0;

/// Levenberg damping decay factor applied on an accepted step: `λ ←
/// max(λ / LM_LAMBDA_DECAY, LM_LAMBDA_FLOOR)` relaxes the trust region back
/// toward the (fast) Gauss–Newton direction. Used by [`lm_damped_accept_sweep`].
const LM_LAMBDA_DECAY: f64 = 10.0;

/// Floor on the Levenberg damping so an accepted-step relaxation can never drive
/// `λ` to zero (which would leave the next damped normal-equation solve singular
/// on a rank-deficient `JᵀJ`). Used by [`lm_damped_accept_sweep`].
const LM_LAMBDA_FLOOR: f64 = 1.0e-12;

/// Relative improvement cut declaring the flow converged: once an accepted step
/// lowers the isometry defect by `≤ LM_IMPROVEMENT_REL_TOL · (1 + defect)` the
/// iterate is treated as a stationary point of the flow objective. Applied by
/// [`lm_damped_accept_sweep`].
const LM_IMPROVEMENT_REL_TOL: f64 = 1.0e-14;

/// Relative step-stall floor for the outer Gauss–Newton loop: after an accepted
/// step, if the squared step norm is `≤ LM_STEP_STALL_REL_FLOOR · (1 + ‖θ‖²)`
/// the iterate has stopped moving and the flow terminates. Applied by both
/// `d = 2` flow cores (torus and sphere-boost).
const LM_STEP_STALL_REL_FLOOR: f64 = 1.0e-24;

/// Minimum per-axis node count of the decoder-recomposition audit grid. The
/// actual count also scales with the basis width (`3·√m` per axis) so the
/// tensor harmonic basis is always Nyquist-oversampled on the audit grid.
pub const TORUS_TRANSPORT_MIN_NODES_PER_AXIS: usize = 48;

/// Identity of one flow mode (for tests and diagnostics): which coordinate
/// component the vector field moves, its integer frequency vector, and its
/// phase (`cos` vs `sin`).
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub struct TorusFlowModeKey {
    pub component: usize,
    pub freq: (i32, i32),
    pub is_cos: bool,
}

/// Per-mode sample of a `d = 2` flow basis at one chart point `t`: the scalar
/// field value `f(t)` (the displacement this mode adds to coordinate
/// `component`) and its gradient `∇f(t)` (the mode's contribution to row
/// `component` of the flow Jacobian `Dφ`). Shared by the torus
/// ([`TorusFlowBasis`]) and the free-patch ([`FreePatchFlowBasis`]) flow
/// families — the isometry-defect Gauss–Newton core
/// (`minimize_isometry_defect_flow`) consumes only this `(component, grad)`
/// contract, so both families descend the same exact optimizer.
#[derive(Debug, Clone, Copy)]
pub struct FlowModeSample {
    pub component: usize,
    pub value: f64,
    pub grad: [f64; 2],
}

/// Truncated Fourier vector-field basis on `T²` with period `period` per
/// axis: `v_{c,(a,b),trig}(t) = e_c · trig(2π(a·t₀ + b·t₁)/period)` for
/// `trig ∈ {sin, cos}`, `c ∈ {0, 1}`, and one frequency representative per
/// antipodal pair (`sin` of `−ω` is `−sin` of `ω`, so both signs would be
/// redundant). The flow map `φ_θ(t) = t + Σ_k θ_k v_k(t)` is automatically a
/// degree-(1,1) torus self-map (`φ(t + period·e_c) = φ(t) + period·e_c` —
/// wrap-around respected by periodicity of the displacement), and any such
/// map with `det Dφ > 0` everywhere is a global diffeomorphism of `T²`.
#[derive(Debug, Clone)]
pub struct TorusFlowBasis {
    pub period: f64,
    /// Canonical frequency representatives: `a > 0`, or `a == 0 && b > 0`.
    freqs: Vec<(i32, i32)>,
}

impl TorusFlowBasis {
    pub fn new(period: f64) -> Result<Self, String> {
        if !(period.is_finite() && period > 0.0) {
            return Err(format!(
                "TorusFlowBasis: period must be finite and positive; got {period}"
            ));
        }
        let h = TORUS_FLOW_MAX_HARMONIC;
        let mut freqs = Vec::new();
        for a in -h..=h {
            for b in -h..=h {
                if a > 0 || (a == 0 && b > 0) {
                    freqs.push((a, b));
                }
            }
        }
        Ok(Self { period, freqs })
    }

    /// Number of flow coefficients `θ`: 2 components × 2 phases × 12
    /// frequency representatives = 48 at the default order.
    pub fn dim(&self) -> usize {
        4 * self.freqs.len()
    }

    /// Sample every mode (value + gradient) at chart point `t`, in `θ` order.
    pub fn mode_samples(&self, t: [f64; 2]) -> Vec<FlowModeSample> {
        let tau = std::f64::consts::TAU;
        let mut out = Vec::with_capacity(self.dim());
        for component in 0..2 {
            for &(a, b) in &self.freqs {
                let w0 = tau * a as f64 / self.period;
                let w1 = tau * b as f64 / self.period;
                let angle = w0 * t[0] + w1 * t[1];
                let s = angle.sin();
                let c = angle.cos();
                out.push(FlowModeSample {
                    component,
                    value: s,
                    grad: [w0 * c, w1 * c],
                });
                out.push(FlowModeSample {
                    component,
                    value: c,
                    grad: [-w0 * s, -w1 * s],
                });
            }
        }
        out
    }

    /// `φ_θ(t)`, wrapped into `[0, period)` per axis.
    pub fn map_point(&self, theta: &[f64], t: [f64; 2]) -> [f64; 2] {
        assert_eq!(theta.len(), self.dim(), "TorusFlowBasis: theta length");
        let mut out = t;
        for (coef, sample) in theta.iter().zip(self.mode_samples(t)) {
            out[sample.component] += coef * sample.value;
        }
        [
            out[0].rem_euclid(self.period),
            out[1].rem_euclid(self.period),
        ]
    }

    /// Flow Jacobian `Dφ_θ(t) = I + Σ_k θ_k Dv_k(t)`, row-major.
    pub fn flow_jacobian(&self, theta: &[f64], t: [f64; 2]) -> [[f64; 2]; 2] {
        assert_eq!(theta.len(), self.dim(), "TorusFlowBasis: theta length");
        let mut jac = [[1.0, 0.0], [0.0, 1.0]];
        for (coef, sample) in theta.iter().zip(self.mode_samples(t)) {
            jac[sample.component][0] += coef * sample.grad[0];
            jac[sample.component][1] += coef * sample.grad[1];
        }
        jac
    }

    /// Minimum of `det Dφ_θ` over the
    /// [`TORUS_FLOW_GUARD_NODES_PER_AXIS`]² check grid — a SAMPLE, not the
    /// continuum minimum. The diffeomorphism guard consumes
    /// [`Self::certified_min_jacobian_det`], which corrects this to a rigorous
    /// everywhere-bound; this stays public because the pair (sample, certified
    /// bound) is what makes the correction's size observable.
    pub fn min_jacobian_det_on_grid(&self, theta: &[f64]) -> f64 {
        let nodes = TORUS_FLOW_GUARD_NODES_PER_AXIS;
        let mut min_det = f64::INFINITY;
        for i in 0..nodes {
            for j in 0..nodes {
                let t = [
                    self.period * i as f64 / nodes as f64,
                    self.period * j as f64 / nodes as f64,
                ];
                let jac = self.flow_jacobian(theta, t);
                let det = jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0];
                min_det = min_det.min(det);
            }
        }
        min_det
    }

    /// Half the guard grid's per-axis spacing — the largest distance from an
    /// arbitrary torus point to the nearest node along one axis. The grid wraps
    /// with the torus, so this is exact, not an edge-effect approximation.
    fn guard_grid_half_spacing(&self) -> f64 {
        0.5 * self.period / TORUS_FLOW_GUARD_NODES_PER_AXIS as f64
    }

    /// The inter-node correction `(h/2)·(L₀ + L₁)` with `L_m = ‖∂_m det Dφ_θ‖_∞`
    /// — the amount by which `det Dφ_θ` can dip below the grid minimum between
    /// nodes. Exactly `0` at `θ = 0`, where `det Dφ ≡ 1`.
    ///
    /// Every `L_m` is closed form in `θ` and the mode frequencies, with no
    /// sampling. Mode `k` is `e_{c(k)}·trig(w₀t₀ + w₁t₁)` with
    /// `w_j = 2π·freq_j/period`, so for EITHER phase
    /// `|∂_j v| ≤ |w_j|` and `|∂_m ∂_j v| ≤ |w_j w_m|`. Writing
    ///
    /// ```text
    /// A[c][j]    = Σ_{k: c(k)=c} |θ_k|·|w_{k,j}|
    /// D[c][j][m] = Σ_{k: c(k)=c} |θ_k|·|w_{k,j}|·|w_{k,m}|
    /// ```
    ///
    /// gives `‖J[c][j]‖_∞ ≤ δ_{cj} + A[c][j]` and `‖∂_m J[c][j]‖_∞ ≤ D[c][j][m]`
    /// for `J = Dφ_θ = I + Σ_k θ_k Dv_k`, and the product rule on
    /// `det = J₀₀J₁₁ − J₀₁J₁₀` gives `L_m`.
    pub fn jacobian_det_grid_correction(&self, theta: &[f64]) -> f64 {
        assert_eq!(theta.len(), self.dim(), "TorusFlowBasis: theta length");
        let tau = std::f64::consts::TAU;
        let mut a = [[0.0_f64; 2]; 2];
        let mut d = [[[0.0_f64; 2]; 2]; 2];
        let mut idx = 0usize;
        for component in 0..2 {
            for &(fa, fb) in &self.freqs {
                let w = [
                    (tau * fa as f64 / self.period).abs(),
                    (tau * fb as f64 / self.period).abs(),
                ];
                // The sin and cos phases of one frequency share `|w|`, so they
                // contribute through the same bounds; walk both in θ order.
                for _phase in 0..2 {
                    let coef = theta[idx].abs();
                    idx += 1;
                    for j in 0..2 {
                        a[component][j] += coef * w[j];
                        for m in 0..2 {
                            d[component][j][m] += coef * w[j] * w[m];
                        }
                    }
                }
            }
        }
        let jmax = |c: usize, j: usize| -> f64 { (if c == j { 1.0 } else { 0.0 }) + a[c][j] };
        let half = self.guard_grid_half_spacing();
        let mut correction = 0.0_f64;
        for m in 0..2 {
            let l_m = d[0][0][m] * jmax(1, 1)
                + jmax(0, 0) * d[1][1][m]
                + d[0][1][m] * jmax(1, 0)
                + jmax(0, 1) * d[1][0][m];
            correction += l_m * half;
        }
        correction
    }

    /// Rigorous LOWER BOUND on `det Dφ_θ` over the WHOLE torus — the quantity
    /// the diffeomorphism guard's contract actually needs, and the one the gate
    /// consumes.
    ///
    /// [`Self::min_jacobian_det_on_grid`] is a sample. `det Dφ_θ` is a trig
    /// polynomial of degree ≤ `2·TORUS_FLOW_MAX_HARMONIC` per axis (each entry
    /// of `Dφ_θ` has degree ≤ `H`, and `det` is a sum of pairwise products), so
    /// between two nodes it can dip below every value the grid saw. The grid is
    /// uniform with spacing `h = period/N` and wraps with the torus, so every
    /// point lies within `h/2` of a node along each axis and
    ///
    /// ```text
    /// min_{T²} det  ≥  min_grid det  −  (h/2)·(L₀ + L₁),   L_m = ‖∂_m det‖_∞,
    /// ```
    ///
    /// with `L_m` closed form in `θ` ([`Self::jacobian_det_grid_correction`]).
    ///
    /// # Why the coefficient bound rather than Bernstein
    ///
    /// Bernstein's inequality gives `L_m ≤ (2π·2H/period)·‖det‖_∞`, which is
    /// correct but loose exactly where this guard operates. Accepted flows are
    /// near-identity: `‖det‖_∞ ≈ 1` while the true `L_m` is `O(|θ|)`. Charging
    /// the correction against `‖det‖_∞` would make it `O(1)` and reject
    /// near-identity flows the guard has no reason to touch. Bernstein is the
    /// right tool when only `‖det‖_∞` is known; here the coefficients ARE the
    /// input, so the product-rule bound is both rigorous and `O(|θ|)`.
    pub fn certified_min_jacobian_det(&self, theta: &[f64]) -> f64 {
        self.min_jacobian_det_on_grid(theta) - self.jacobian_det_grid_correction(theta)
    }
}

/// The exact, image-frozen minimum-isometry-defect flow reparameterization of
/// one `d = 2` torus atom.
#[derive(Debug, Clone)]
pub struct TorusIsometryFlowReparameterization {
    /// Canonical per-row coordinates `t̃_i = φ_θ(t_i)`, shape `(n, 2)`,
    /// wrapped into `[0, period)` per axis.
    pub new_row_coords: Array2<f64>,
    /// Recomposed decoder `B̃ = T · B`, shape `(m, p)` — the exact LS refit of
    /// the original decoded image on the audit grid against the basis at the
    /// transported grid (so `γ̃ = γ ∘ φ⁻¹` without ever forming `φ⁻¹`).
    pub new_decoder: Array2<f64>,
    /// The `(m, m)` basis transport `T` with `Φ(φ(u)) · T ≈ Φ(u)` on the
    /// audit grid — the same congruence object the `d = 1` path and the
    /// affine gauge canonicalization use to transport the smoothness Gram.
    pub decoder_transport: Array2<f64>,
    /// Optimal flow coefficients `θ` (layout per
    /// `TorusFlowBasis::mode_layout`).
    pub flow_theta: Vec<f64>,
    /// Isometry defect `E(0)` of the fitted chart (identity flow).
    pub defect_initial: f64,
    /// Isometry defect `E(θ)` of the canonical chart. Strictly below
    /// `defect_initial` (the pass refuses no-improvement flows).
    pub defect_final: f64,
    /// The profiled global metric scale `c` at the optimum (the canonical
    /// chart's pullback metric is `≈ c·ḡ·I`).
    pub profiled_metric_scale: f64,
    /// Certified lower bound on `det Dφ_θ` over the whole torus
    /// ([`TorusFlowBasis::certified_min_jacobian_det`]) — NOT the raw grid
    /// minimum. Always `> TORUS_FLOW_DIFFEO_MIN_DET` when `Some(..)` is
    /// returned, so a folded chart is refused upstream on a bound that holds
    /// everywhere rather than on the nodes it happened to sample.
    pub min_flow_jacobian_det: f64,
    /// Max-abs recomposition error on the audit grid, relative to the image
    /// scale. Always `≤ CHART_RECOMPOSITION_REL_TOL` when `Some(..)`.
    pub recomposition_residual: f64,
}

/// Speed-uniformity defect above which the in-loop unit-speed retraction fires.
/// Below it the chart is already ~arc-length and re-gauging is skipped (keeps the
/// per-refresh cost negligible and makes the retraction idempotent at a boundary).
///
/// PRICED (#2071): the defect is a coefficient of variation
/// `stddev(‖γ'‖)/mean(‖γ'‖)` — dimensionless, `0` ⟺ exactly unit-speed — so
/// this is a numerical-idempotency FLOOR, not a modeling knob. The retraction is
/// objective-invariant for data-fit and intrinsic smoothness (image-frozen,
/// `BᵀSB` transported); ONLY the ARD coordinate prior re-evaluates at the
/// canonical `t̃`, so the threshold governs a single question: is the chart far
/// enough from arc-length that re-gauging would move the prior's evaluation point
/// measurably? The reparameterization is built by integrating the sampled speed
/// field, so it cannot itself drive the CV below its own quadrature/interpolation
/// resolution — many orders above machine eps (`~2e-16`) yet far below any
/// non-uniformity that perturbs the ARD prior. `1e-6` sits in that gap (same
/// numerical genus as the sibling [`CHART_RECOMPOSITION_REL_TOL`] `1e-9`
/// image-scale honesty tolerance, loosened to CV scale because the defect is a
/// first-derivative statistic): below it the retraction is a numerical no-op (its
/// own residual-CV floor), above it there is genuine, correctable non-uniformity.
/// What breaks at 10×: at `1e-5` a chart with a real, removable `~1e-5` speed CV
/// is skipped, so the ARD coordinate prior is evaluated on a slightly
/// non-arc-length gauge; at `1e-7` the early-out almost never fires and the
/// reparam is attempted on charts already uniform to its own resolution (a
/// wasted, idempotent re-gauge each refresh boundary).
pub const UNIT_SPEED_INLOOP_DEFECT_TOL: f64 = 1.0e-6;

/// Coefficient of variation `stddev(‖γ'‖)/mean(‖γ'‖)` of the decoder-curve speed
/// field. `0` ⇔ already unit-speed (arc-length). Cheap in-loop early-out signal;
/// returns `+∞` on a collapsed/non-positive mean so the caller defers to the full
/// reparameterization (which honest-skips a degenerate chart).
fn speed_uniformity_defect(speeds: &[f64]) -> f64 {
    let n = speeds.len();
    if n == 0 {
        return 0.0;
    }
    let mean = speeds.iter().sum::<f64>() / n as f64;
    if !(mean > 0.0) {
        return f64::INFINITY;
    }
    let var = speeds.iter().map(|&s| (s - mean) * (s - mean)).sum::<f64>() / n as f64;
    var.sqrt() / mean
}

/// #2022 — in-loop unit-speed RETRACTION for one `d = 1` atom, the primitive the
/// fit loop calls at each chart-refresh boundary. Returns `Ok(None)` (honest skip:
/// already ~unit-speed, degenerate chart, basis not closed under the reparam, or
/// image drift above tolerance) or `Ok(Some(repar))` to apply. Objective-invariant
/// for data-fit + intrinsic smoothness (image-frozen; transport preserves `BᵀSB`);
/// the ARD coordinate prior re-evaluates at the canonical `t̃`, so this MUST be
/// applied only at refresh boundaries, never inside a line search.
pub fn unit_speed_retraction(
    evaluator: &dyn SaeBasisEvaluator,
    decoder: ArrayView2<'_, f64>,
    row_coords: ArrayView1<'_, f64>,
    topology: &CanonicalChartTopology,
) -> Result<Option<UnitSpeedReparameterization>, String> {
    let n = row_coords.len();
    if n == 0 {
        return Ok(None);
    }
    // Cheap early-out: evaluate the jet at the CURRENT coords and skip the reparam
    // when the chart is already ~arc-length. `evaluate` returns `(phi, jet)`; only
    // the jet is needed for the speed field.
    let mut coords2 = Array2::<f64>::zeros((n, 1));
    for i in 0..n {
        coords2[[i, 0]] = row_coords[i];
    }
    let (_phi, jet) = evaluator.evaluate(coords2.view())?;
    let speeds = curve_speeds(&jet, decoder)?;
    if speeds.iter().all(|s| s.is_finite())
        && speed_uniformity_defect(&speeds) < UNIT_SPEED_INLOOP_DEFECT_TOL
    {
        return Ok(None);
    }
    unit_speed_reparameterization(evaluator, decoder, row_coords, topology)
}

/// #2081 — the unit-speed (arc-length) DEFECT of a fitted `d = 1` chart: the
/// coefficient of variation `stddev(‖γ'(t)‖) / mean(‖γ'(t)‖)` of the decoder
/// curve speed over a uniform grid of the chart's canonical domain. `0` ⟺ the
/// parameterization is exactly arc-length (unit-speed); a positive value is the
/// same gauge the in-loop retraction's early-out reads
/// (`speed_uniformity_defect`), promoted here from a control signal to a
/// REPORTED per-atom certificate quantity.
///
/// It is a pure property of `(decoder, basis)` — measured on a UNIFORM latent
/// grid, not on the data rows — so it certifies the CHART parameterization
/// independent of where the data falls (the complement of the coordinate
/// UNIFORMITY statistic, which certifies the coordinate distribution). This is
/// exactly the #2081 lever: an `n_basis = 4` Fourier chart can reconstruct its
/// ring at high EV while its parameterization is far from arc-length, and EV
/// cannot see that — this defect can.
///
/// The domain is the circle's full period `[0, period)` or, for an interval
/// chart, the fitted coordinate range `[min t_i, max t_i]` — the SAME `(lo, hi)`
/// [`unit_speed_reparameterization`] integrates over. Returns `Ok(None)` (honest
/// skip, never a fabricated number) on a degenerate chart: empty basis/decoder,
/// a non-finite / collapsed domain, or a non-finite / zero-mean speed field.
pub fn chart_unit_speed_defect(
    evaluator: &dyn SaeBasisEvaluator,
    decoder: ArrayView2<'_, f64>,
    row_coords: ArrayView1<'_, f64>,
    topology: &CanonicalChartTopology,
) -> Result<Option<f64>, String> {
    let m = decoder.nrows();
    let p = decoder.ncols();
    if m == 0 || p == 0 {
        return Ok(None);
    }
    // Canonical domain `[lo, hi]` — identical to the reparameterization's.
    let (lo, hi) = match topology {
        CanonicalChartTopology::Circle { period } => {
            if !(period.is_finite() && *period > 0.0) {
                return Ok(None);
            }
            (0.0, *period)
        }
        CanonicalChartTopology::Interval => {
            let mut t_min = f64::INFINITY;
            let mut t_max = f64::NEG_INFINITY;
            for &t in row_coords.iter() {
                if !t.is_finite() {
                    return Ok(None);
                }
                t_min = t_min.min(t);
                t_max = t_max.max(t);
            }
            let scale = t_min.abs().max(t_max.abs()).max(1.0);
            if !(t_max - t_min > 1.0e-12 * scale) {
                return Ok(None);
            }
            (t_min, t_max)
        }
    };
    // Uniform latent grid over the domain. The speed field is band-limited by
    // the basis, so reusing the arc-length quadrature cell count oversamples it
    // — no new resolution constant. A circle omits the wrap-duplicate endpoint;
    // an interval samples both ends.
    let cells = ARC_LENGTH_GRID_CELLS;
    let step = (hi - lo) / cells as f64;
    let nodes = match topology {
        CanonicalChartTopology::Circle { .. } => cells,
        CanonicalChartTopology::Interval => cells + 1,
    };
    let mut grid = Array2::<f64>::zeros((nodes, 1));
    for j in 0..nodes {
        grid[[j, 0]] = lo + j as f64 * step;
    }
    let (grid_phi, jet) = evaluator.evaluate(grid.view())?;
    if grid_phi.ncols() != m {
        return Err(format!(
            "chart_unit_speed_defect: evaluator basis width {} != decoder rows {m}",
            grid_phi.ncols()
        ));
    }
    let speeds = curve_speeds(&jet, decoder)?;
    if speeds.iter().any(|s| !s.is_finite()) {
        return Ok(None);
    }
    let defect = speed_uniformity_defect(&speeds);
    if !defect.is_finite() {
        return Ok(None);
    }
    Ok(Some(defect))
}

/// #2081 — the honest arc-length (unit-speed) coordinate of a fitted `d = 1`
/// chart, computed as a PURE READ: the per-row canonical coordinate
/// `u_i = s(t_i)/L` plus the decoder-curve speed profile, WITHOUT the decoder
/// mutation or the [`CHART_RECOMPOSITION_REL_TOL`] recomposition gate that
/// [`unit_speed_reparameterization`] applies.
///
/// The recomposition gate rightly protects a decoder REWRITE — a finite Fourier
/// basis generically cannot re-express a warped circle within `1e-9`, so the
/// mutating canonicalization honestly refuses. But the arc-length COORDINATE is
/// a property of the fitted curve alone (`s(t) = ∫‖γ'‖`), so it is always
/// computable and must always be reportable: this is the coordinate every
/// downstream angle/dose/adjacency claim should read instead of the gauge-
/// arbitrary raw `t` (#2081). The map is the SAME `t ↦ span·s(t)/L` the mutating
/// path uses, so `coords_u_arc` equals the reparameterized coordinates whenever
/// that path DID commit — this is its ungated sibling.
///
/// Returns `Ok(None)` (honest skip, never a fabricated coordinate) on a
/// degenerate chart: empty basis/decoder, a non-finite / collapsed domain, or a
/// non-finite / zero total arc length.
#[derive(Debug, Clone)]
pub struct ChartArcLengthReading {
    /// Per-row canonical coordinate `u_i = s(t_i)/L ∈ [0, 1)` (circle, wrapped)
    /// or `∈ [0, 1]` (interval) — the arc-length reparameterization normalized
    /// to the unit interval. The honest angle/position every coordinate consumer
    /// should read in place of the raw `t_i`.
    pub coords_u_arc: Array1<f64>,
    /// Speed coefficient of variation `stddev(‖γ'‖)/mean(‖γ'‖)` on the uniform
    /// quadrature grid. `0` ⟺ the raw chart is already exactly arc-length.
    pub speed_cv: f64,
    /// RMS of `log(‖γ'‖/mean)` over the grid — the scale-invariant log-speed
    /// spread (a lognormal-flavoured complement to `speed_cv`, robust to a few
    /// fast cells dominating the variance).
    pub log_speed_rms: f64,
    /// `min ‖γ'‖ / mean ‖γ'‖` over the grid. Approaching `0` means the chart
    /// nearly collapses somewhere (a flat spot in `u`), so the arc-length
    /// reparameterization stops being a well-conditioned diffeomorphism.
    pub min_speed_over_mean: f64,
    /// `max ‖γ'‖ / mean ‖γ'‖` over the grid.
    pub max_speed_over_mean: f64,
    /// Total arc length `L` of the decoder curve over the canonical domain.
    pub total_arc_length: f64,
}

// Gauge-theory reading: `coords_u_arc` is the gauge-invariant OBSERVABLE this
// module ultimately exists to expose. It is computed by the same map as the
// mutating slice ([`unit_speed_reparameterization`]) but without ever touching
// the decoder or requiring the recomposition gate to pass — the arc-length
// position along a fitted curve is well-defined from `(decoder, basis)` alone,
// independent of which point of the `Diff(M)` orbit the fit happened to land
// on. Any downstream consumer reading angle/dose/adjacency off the raw `t_i`
// is reading a gauge-ARBITRARY quantity; reading `coords_u_arc` instead reads
// the coordinate that is actually pinned by the physics.
pub fn chart_arclength_coordinates(
    evaluator: &dyn SaeBasisEvaluator,
    decoder: ArrayView2<'_, f64>,
    row_coords: ArrayView1<'_, f64>,
    topology: &CanonicalChartTopology,
) -> Result<Option<ChartArcLengthReading>, String> {
    let n = row_coords.len();
    let m = decoder.nrows();
    let p = decoder.ncols();
    if n == 0 || m == 0 || p == 0 {
        return Ok(None);
    }
    for &t in row_coords.iter() {
        if !t.is_finite() {
            return Ok(None);
        }
    }

    // Canonical quadrature domain `[lo, hi]` and unit target span — IDENTICAL to
    // `unit_speed_reparameterization`, but normalized to `[0, 1)` (we report the
    // unit-interval coordinate, not the native-period one).
    let (lo, hi) = match topology {
        CanonicalChartTopology::Circle { period } => {
            if !(period.is_finite() && *period > 0.0) {
                return Ok(None);
            }
            (0.0, *period)
        }
        CanonicalChartTopology::Interval => {
            let mut t_min = f64::INFINITY;
            let mut t_max = f64::NEG_INFINITY;
            for &t in row_coords.iter() {
                t_min = t_min.min(t);
                t_max = t_max.max(t);
            }
            let scale = t_min.abs().max(t_max.abs()).max(1.0);
            if !(t_max - t_min > 1.0e-12 * scale) {
                return Ok(None);
            }
            (t_min, t_max)
        }
    };

    // Speed field on the composite-Simpson grid (nodes + midpoints).
    let cells = ARC_LENGTH_GRID_CELLS;
    let h = (hi - lo) / cells as f64;
    let mut quad_coords = Array2::<f64>::zeros((2 * cells + 1, 1));
    for j in 0..=cells {
        quad_coords[[2 * j, 0]] = lo + j as f64 * h;
        if j < cells {
            quad_coords[[2 * j + 1, 0]] = lo + (j as f64 + 0.5) * h;
        }
    }
    let (_grid_phi, grid_jet) = evaluator.evaluate(quad_coords.view())?;
    let speeds = curve_speeds(&grid_jet, decoder)?;
    if speeds.iter().any(|s| !s.is_finite()) {
        return Ok(None);
    }

    // Cumulative arc length at the cell nodes (composite Simpson).
    let mut cumulative = vec![0.0_f64; cells + 1];
    for j in 0..cells {
        let f0 = speeds[2 * j];
        let fm = speeds[2 * j + 1];
        let f1 = speeds[2 * j + 2];
        cumulative[j + 1] = cumulative[j] + h * (f0 + 4.0 * fm + f1) / 6.0;
    }
    let total = cumulative[cells];
    if !(total.is_finite() && total > 0.0) {
        return Ok(None);
    }
    let inv_total = 1.0 / total;

    // The canonical chart map `t ↦ s(t)/L`, normalized to `[0, 1)`.
    let map_unit = |t: f64| -> f64 {
        let local = match topology {
            CanonicalChartTopology::Circle { period } => (t - lo).rem_euclid(*period),
            CanonicalChartTopology::Interval => (t - lo).clamp(0.0, hi - lo),
        };
        let cell = ((local / h).floor() as usize).min(cells - 1);
        let x = local - cell as f64 * h;
        let s = cumulative[cell]
            + partial_cell_arc(
                speeds[2 * cell],
                speeds[2 * cell + 1],
                speeds[2 * cell + 2],
                h,
                x,
            );
        let u = s * inv_total;
        match topology {
            CanonicalChartTopology::Circle { .. } => u.rem_euclid(1.0),
            CanonicalChartTopology::Interval => u.clamp(0.0, 1.0),
        }
    };
    let coords_u_arc = Array1::from_iter(row_coords.iter().map(|&t| map_unit(t)));

    // Speed profile over the same grid samples (nodes + midpoints).
    let mean = speeds.iter().sum::<f64>() / speeds.len() as f64;
    if !(mean > 0.0) {
        return Ok(None);
    }
    let speed_cv = speed_uniformity_defect(&speeds);
    let mut log_sq = 0.0_f64;
    let mut min_over_mean = f64::INFINITY;
    let mut max_over_mean = f64::NEG_INFINITY;
    for &s in &speeds {
        let ratio = s / mean;
        min_over_mean = min_over_mean.min(ratio);
        max_over_mean = max_over_mean.max(ratio);
        // `s` is finite and non-negative; a genuine zero-speed cell drives
        // `min_over_mean → 0` (caught by the diffeomorphism-floor verdict) and is
        // excluded from the log spread so it does not poison it with `-inf`.
        if ratio > 0.0 {
            let l = ratio.ln();
            log_sq += l * l;
        }
    }
    let log_speed_rms = (log_sq / speeds.len() as f64).sqrt();
    if !(speed_cv.is_finite() && log_speed_rms.is_finite() && min_over_mean.is_finite()) {
        return Ok(None);
    }

    Ok(Some(ChartArcLengthReading {
        coords_u_arc,
        speed_cv,
        log_speed_rms,
        min_speed_over_mean: min_over_mean,
        max_speed_over_mean: max_over_mean,
        total_arc_length: total,
    }))
}

/// State of the flow objective at one `θ`: the defect, the profiled scale,
/// and the per-row (possibly reference-whitened) flow Jacobians `A_i`
/// (row-major `[a00, a01, a10, a11]`) the Gauss--Newton rows are built from.
struct FlowObjectiveState {
    defect: f64,
    scale: f64,
    a_rows: Vec<[f64; 4]>,
}

/// Outcome of one [`lm_damped_accept_sweep`] trust step, reported back to the
/// outer Gauss–Newton loop so it can decide whether to continue.
struct LmTrustStep {
    /// A strict-descent, fold-free candidate was accepted this sweep.
    accepted: bool,
    /// The accepted step improved the defect by less than the relative
    /// convergence cut ([`LM_IMPROVEMENT_REL_TOL`]) — treat as stationary.
    converged: bool,
    /// Squared Euclidean norm of the accepted step (0 if none was accepted);
    /// the outer loop compares it against the relative step-stall floor.
    step_norm_sq: f64,
}

/// The shared Levenberg–Marquardt damped-accept trust step for the analytic
/// `d = 2` flow cores. Given the current normal-equation blocks (`JᵀJ`, `Jᵀr`) it
/// escalates the damping `λ` until it finds a fold-free strict-descent
/// candidate or exhausts [`TORUS_FLOW_GN_MAX_REJECTS`] rejections:
///
/// * on a Cholesky failure or a folded / non-improving candidate it grows the
///   damping (`λ ← LM_LAMBDA_GROWTH · λ`) and counts a rejection;
/// * on an accepted candidate it commits `theta`/`state`, relaxes the damping
///   (`λ ← max(λ / LM_LAMBDA_DECAY, LM_LAMBDA_FLOOR)`), and flags convergence
///   when the relative improvement drops below [`LM_IMPROVEMENT_REL_TOL`].
///
/// The `λ` here is decayed on accept and carried across outer Gauss–Newton
/// iterations (a stateful trust-region damping), which is why this is NOT the
/// generic ridge-escalation optimizer primitive. The only flow-specific piece
/// is `eval_candidate`, which folds the family's diffeomorphism guard and its
/// defect evaluation into one closure returning `None` for a folded or
/// out-of-band candidate. The numerics are bit-for-bit identical to the two
/// former hand-rolled loops it replaces.
fn lm_damped_accept_sweep(
    jtj: &Array2<f64>,
    jtr: &Array2<f64>,
    q: usize,
    lambda: &mut f64,
    theta: &mut Vec<f64>,
    state: &mut FlowObjectiveState,
    eval_candidate: impl Fn(&[f64]) -> Option<FlowObjectiveState>,
) -> LmTrustStep {
    let mut rejects = 0usize;
    let mut accepted = false;
    let mut converged = false;
    let mut step_norm_sq = 0.0_f64;
    while rejects < TORUS_FLOW_GN_MAX_REJECTS {
        let mut damped = jtj.clone();
        for d in 0..q {
            damped[[d, d]] += *lambda * (1.0 + jtj[[d, d]]);
        }
        let factor = match damped.cholesky(FaerSide::Lower) {
            Ok(factor) => factor,
            Err(_) => {
                *lambda *= LM_LAMBDA_GROWTH;
                rejects += 1;
                continue;
            }
        };
        let mut neg_jtr = jtr.clone();
        neg_jtr.mapv_inplace(|v| -v);
        let delta = factor.solve_mat(&neg_jtr);
        let mut candidate = theta.clone();
        step_norm_sq = 0.0;
        for k in 0..q {
            candidate[k] += delta[[k, 0]];
            step_norm_sq += delta[[k, 0]] * delta[[k, 0]];
        }
        match eval_candidate(&candidate) {
            Some(next) if next.defect < state.defect => {
                let improvement = state.defect - next.defect;
                *theta = candidate;
                *state = next;
                accepted = true;
                *lambda = (*lambda / LM_LAMBDA_DECAY).max(LM_LAMBDA_FLOOR);
                if improvement <= LM_IMPROVEMENT_REL_TOL * (1.0 + state.defect) {
                    // Converged: the accepted step no longer moves E.
                    converged = true;
                }
                break;
            }
            Some(..) | None => {
                *lambda *= LM_LAMBDA_GROWTH;
                rejects += 1;
            }
        }
    }
    LmTrustStep {
        accepted,
        converged,
        step_norm_sq,
    }
}

/// Evaluate the isometry-defect objective at `θ` (see
/// [`torus_isometry_flow_reparameterization`] for the derivation). Returns
/// `None` when the profiled scale degenerates (`c ≤ 0` or non-finite).
///
/// `row_base` is the per-row base `A0_i` (row-major
/// `[a00, a01, a10, a11]`) the flow modes are added on top of: the effective
/// Jacobian is `A_i = A0_i + Σ_k θ_k W_{ik}`. The torus/free-patch families
/// pass `A0_i = I`. The sphere has moved-point reference whitening and therefore
/// uses its own exact evaluator/Jacobian below.
fn evaluate_flow_defect(
    theta: &[f64],
    row_modes: &[Vec<FlowModeSample>],
    row_base: &[[f64; 4]],
    ghat: &[[f64; 3]],
    ghat_norm_sq: f64,
) -> Option<FlowObjectiveState> {
    let n = row_modes.len();
    let mut a_rows = Vec::with_capacity(n);
    let mut cross = 0.0_f64;
    for (modes, base) in row_modes.iter().zip(row_base.iter()) {
        let mut a = *base;
        for (coef, sample) in theta.iter().zip(modes.iter()) {
            a[2 * sample.component] += coef * sample.grad[0];
            a[2 * sample.component + 1] += coef * sample.grad[1];
        }
        a_rows.push(a);
    }
    for (a, g) in a_rows.iter().zip(ghat.iter()) {
        // AᵀA in symmetric storage [m00, m11, m01].
        let m00 = a[0] * a[0] + a[2] * a[2];
        let m11 = a[1] * a[1] + a[3] * a[3];
        let m01 = a[0] * a[1] + a[2] * a[3];
        cross += m00 * g[0] + m11 * g[1] + 2.0 * m01 * g[2];
    }
    let scale = cross / ghat_norm_sq;
    if !(scale.is_finite() && scale > 0.0) {
        return None;
    }
    let mut defect = 0.0_f64;
    for (a, g) in a_rows.iter().zip(ghat.iter()) {
        let m00 = a[0] * a[0] + a[2] * a[2];
        let m11 = a[1] * a[1] + a[3] * a[3];
        let m01 = a[0] * a[1] + a[2] * a[3];
        let r00 = m00 - scale * g[0];
        let r11 = m11 - scale * g[1];
        let r01 = m01 - scale * g[2];
        defect += r00 * r00 + r11 * r11 + 2.0 * r01 * r01;
    }
    if !defect.is_finite() {
        return None;
    }
    Some(FlowObjectiveState {
        defect,
        scale,
        a_rows,
    })
}

/// Outcome of the shared isometry-defect flow minimizer: the optimal flow
/// coefficients and the bracket of defects/scale they achieve. Returned only
/// when a strict, fold-free improvement over the identity flow was found.
struct FlowMinimization {
    theta: Vec<f64>,
    defect_initial: f64,
    defect_final: f64,
    profiled_scale: f64,
}

/// Exact damped Gauss–Newton for the `d = 2` isometry defect
/// `E(θ) = Σ_i ‖A_iᵀA_i − c·Ĝ_i‖²_F` over the flow coefficients `θ`, shared by
/// the torus and free-patch flow families (see
/// [`torus_isometry_flow_reparameterization`] for the full derivation of the
/// residual, the profiled scale `c`, and the analytic Gauss–Newton Jacobian).
///
/// The flow family enters ONLY through `row_modes` (the per-row mode samples
/// `W_{ik} = Dv_k(t_i)` and displacements) and the `min_det_on_grid` guard
/// closure, so the two families descend the identical optimizer with the
/// identical strict-descent + diffeomorphism accept test. The minimization
/// starts at `θ = 0` (`det Dφ = 1` everywhere) and never accepts a candidate
/// whose `min det Dφ_θ ≤ min_det` on the guard grid, so the iterate can never
/// walk through a fold. Returns `None` (honest skip — no lossy or folded swap)
/// when the identity chart is already isometric, the profiled scale
/// degenerates, or no strict improvement is reachable within the family.
fn minimize_isometry_defect_flow(
    row_modes: &[Vec<FlowModeSample>],
    row_base: &[[f64; 4]],
    ghat: &[[f64; 3]],
    ghat_norm_sq: f64,
    q: usize,
    min_det: f64,
    min_det_on_grid: &dyn Fn(&[f64]) -> f64,
) -> Option<FlowMinimization> {
    let n = row_modes.len();
    let mut theta = vec![0.0_f64; q];
    let mut state = evaluate_flow_defect(&theta, row_modes, row_base, ghat, ghat_norm_sq)?;
    let defect_initial = state.defect;
    if !(defect_initial > 0.0) {
        // Already exactly isometric — nothing to canonicalize.
        return None;
    }
    let sqrt2 = std::f64::consts::SQRT_2;
    let mut lambda = 1.0e-4_f64;
    let mut any_accepted = false;
    for iteration in 0..TORUS_FLOW_GN_MAX_ITERS {
        if iteration + 1 == TORUS_FLOW_GN_MAX_ITERS {
            break;
        }
        // Residual r and Gauss–Newton Jacobian J at the current θ.
        let mut jmat = Array2::<f64>::zeros((3 * n, q));
        let mut rcol = Array2::<f64>::zeros((3 * n, 1));
        for (i, (a, g)) in state.a_rows.iter().zip(ghat.iter()).enumerate() {
            let m00 = a[0] * a[0] + a[2] * a[2];
            let m11 = a[1] * a[1] + a[3] * a[3];
            let m01 = a[0] * a[1] + a[2] * a[3];
            rcol[[3 * i, 0]] = m00 - state.scale * g[0];
            rcol[[3 * i + 1, 0]] = m11 - state.scale * g[1];
            rcol[[3 * i + 2, 0]] = sqrt2 * (m01 - state.scale * g[2]);
            for (k, sample) in row_modes[i].iter().enumerate() {
                // W_{ik} has single nonzero row `component` = grad, so
                // M = W_{ik}ᵀ A_i has entries M_{ab} = grad[a]·A[component, b]
                // and S = M + Mᵀ.
                let ac0 = a[2 * sample.component];
                let ac1 = a[2 * sample.component + 1];
                let s00 = 2.0 * sample.grad[0] * ac0;
                let s11 = 2.0 * sample.grad[1] * ac1;
                let s01 = sample.grad[0] * ac1 + sample.grad[1] * ac0;
                jmat[[3 * i, k]] = s00;
                jmat[[3 * i + 1, k]] = s11;
                jmat[[3 * i + 2, k]] = sqrt2 * s01;
            }
        }
        let jtj = fast_ata(&jmat);
        let jtr = fast_atb(&jmat, &rcol);

        // Levenberg-damped step with the diffeomorphism guard in the accept
        // test: only strict-descent, fold-free candidates are ever taken.
        let step = lm_damped_accept_sweep(
            &jtj,
            &jtr,
            q,
            &mut lambda,
            &mut theta,
            &mut state,
            |candidate| {
                if min_det_on_grid(candidate) <= min_det {
                    None
                } else {
                    evaluate_flow_defect(candidate, row_modes, row_base, ghat, ghat_norm_sq)
                }
            },
        );
        any_accepted |= step.accepted;
        if !step.accepted {
            break;
        }
        if step.converged {
            break;
        }
        let theta_norm_sq: f64 = theta.iter().map(|v| v * v).sum();
        if step.step_norm_sq <= LM_STEP_STALL_REL_FLOOR * (1.0 + theta_norm_sq) {
            break;
        }
    }
    if !any_accepted || !(state.defect < defect_initial) {
        // No strict improvement within the flow family: the fitted chart is
        // already the canonical representative — honest skip.
        return None;
    }
    Some(FlowMinimization {
        theta,
        defect_initial,
        defect_final: state.defect,
        profiled_scale: state.scale,
    })
}

/// Extract the fitted pullback metric `G_i = J(t_i)ᵀ J(t_i)` (symmetric storage
/// `[g00, g11, g01]`) at every row of a `d = 2` atom from the exact decoder jet,
/// together with the geometric-mean metric scale `ḡ = exp(mean_i ½ log det G_i)`
/// used by every `d = 2` defect for its scale-invariant normalization. Shared by
/// the torus, free-patch, and sphere defect paths — the single source of truth
/// for the pullback-metric extraction.
///
/// Returns `Ok(None)` (honest refusal) on a degenerate chart: empty rows/basis,
/// a rank-deficient pullback metric (`det G_i ≤ 0`) anywhere — the chart is
/// collapsed along some direction there, so no isometric representative exists —
/// or a non-finite geometric-mean scale.
fn extract_pullback_metric_d2(
    label: &str,
    evaluator: &dyn SaeBasisEvaluator,
    decoder: ArrayView2<'_, f64>,
    row_coords: ArrayView2<'_, f64>,
) -> Result<Option<(Vec<[f64; 3]>, f64)>, String> {
    let n = row_coords.nrows();
    let m = decoder.nrows();
    let p = decoder.ncols();
    let (row_phi, row_jet) = evaluator.evaluate(row_coords)?;
    if row_phi.ncols() != m || row_jet.dim() != (n, m, 2) {
        return Err(format!(
            "{label}: evaluator returned basis {:?} / jet {:?}; expected width {m}, latent_dim 2",
            row_phi.dim(),
            row_jet.dim()
        ));
    }
    let mut g_rows: Vec<[f64; 3]> = Vec::with_capacity(n);
    let mut log_det_sum = 0.0_f64;
    let mut tangent0 = vec![0.0_f64; p];
    let mut tangent1 = vec![0.0_f64; p];
    for row in 0..n {
        for slot in tangent0.iter_mut() {
            *slot = 0.0;
        }
        for slot in tangent1.iter_mut() {
            *slot = 0.0;
        }
        for bm in 0..m {
            let d0 = row_jet[[row, bm, 0]];
            let d1 = row_jet[[row, bm, 1]];
            if d0 == 0.0 && d1 == 0.0 {
                continue;
            }
            for j in 0..p {
                let b = decoder[[bm, j]];
                tangent0[j] += d0 * b;
                tangent1[j] += d1 * b;
            }
        }
        let mut g00 = 0.0_f64;
        let mut g11 = 0.0_f64;
        let mut g01 = 0.0_f64;
        for j in 0..p {
            g00 += tangent0[j] * tangent0[j];
            g11 += tangent1[j] * tangent1[j];
            g01 += tangent0[j] * tangent1[j];
        }
        let det = g00 * g11 - g01 * g01;
        if !(det.is_finite() && det > 0.0) {
            return Ok(None);
        }
        log_det_sum += 0.5 * det.ln();
        g_rows.push([g00, g11, g01]);
    }
    let g_bar = (log_det_sum / n as f64).exp();
    if !(g_bar.is_finite() && g_bar > 0.0) {
        return Ok(None);
    }
    Ok(Some((g_rows, g_bar)))
}

/// Normalize the fitted metric rows against the **flat** reference `g_ref = I`:
/// `Ĝ_i = G_i / ḡ` plus the reference-norm² `Σ_i ‖Ĝ_i‖²_F` the isometry-defect
/// Gauss–Newton profiles its global scale against. Shared by the torus and
/// free-patch families (both pin to a flat uniform-speed reference); the sphere
/// uses the `diag(1, cos²lat)` reference instead and normalizes inline.
fn flat_normalized_metric(g_rows: &[[f64; 3]], g_bar: f64) -> Option<(Vec<[f64; 3]>, f64)> {
    let mut ghat: Vec<[f64; 3]> = Vec::with_capacity(g_rows.len());
    let mut ghat_norm_sq = 0.0_f64;
    for g in g_rows {
        let h = [g[0] / g_bar, g[1] / g_bar, g[2] / g_bar];
        ghat_norm_sq += h[0] * h[0] + h[1] * h[1] + 2.0 * h[2] * h[2];
        ghat.push(h);
    }
    if !(ghat_norm_sq.is_finite() && ghat_norm_sq > 0.0) {
        return None;
    }
    Some((ghat, ghat_norm_sq))
}

/// Compute the minimum-isometry-defect flow reparameterization of a fitted
/// `d = 2` torus atom: the canonical per-row coordinates `t̃_i = φ_θ(t_i)`
/// and the exactly-recomposed decoder.
///
/// # The defect functional (and why it is exactly the issue's isometry defect)
///
/// The new chart is `t̃ = φ(t)` with new decoded map `γ̃ = γ ∘ φ⁻¹` (image
/// frozen), so the pullback metric in the canonical chart at `φ(t)` is
/// `G̃(φ(t)) = Dφ(t)⁻ᵀ G(t) Dφ(t)⁻¹` where `G(t) = J(t)ᵀJ(t)` is the fitted
/// pullback metric (`J` = decoder Jacobian, from the exact `(Φ, ∂Φ)` jet).
/// The canonical chart is isometric to the flat reference torus up to a
/// global scale `s` iff `G̃ ≡ s·I`, i.e. iff `Dφᵀ Dφ ≡ G / s`. Measuring the
/// defect on THIS side of the equivalence,
///
/// ```text
/// E(θ) = Σ_i ‖ A_iᵀ A_i − c · Ĝ_i ‖²_F ,   A_i = Dφ_θ(t_i) = I + Σ_k θ_k W_{ik} ,
/// ```
///
/// keeps the residual polynomial (quadratic) in `θ` — no `Dφ⁻¹` anywhere.
/// Here `W_{ik} = Dv_k(t_i)` are the constant per-row mode Jacobians,
/// `Ĝ_i = G_i / ḡ` with `ḡ = exp( mean_i ½·log det G_i )` the geometric-mean
/// metric scale of the fitted rows (the scale-invariant normalization — the
/// `d = 2` analogue of the `d = 1` module's rescale-by-total-arc-length), and
/// `c = c(θ)` the analytically profiled residual global scale
///
/// ```text
/// c(θ) = Σ_i ⟨A_iᵀA_i, Ĝ_i⟩_F / Σ_i ‖Ĝ_i‖²_F   (the exact argmin over c),
/// ```
///
/// which absorbs the (second-order) arithmetic-vs-geometric mean mismatch so
/// the defect is exactly scale-invariant: a chart isometric up to ANY global
/// scale has `E = 0`.
///
/// # The analytic gradient / Gauss–Newton (FD-free)
///
/// With `R_i = A_iᵀA_i − c·Ĝ_i` (symmetric) and `c` profiled,
/// `∂E/∂c = 0` at `c(θ)` (envelope theorem), so the exact gradient treats `c`
/// as fixed:
///
/// ```text
/// ∂R_i/∂θ_k = W_{ik}ᵀ A_i + A_iᵀ W_{ik}
/// ∂E/∂θ_k   = 2 Σ_i ⟨R_i, W_{ik}ᵀA_i + A_iᵀW_{ik}⟩_F = 4 Σ_i ⟨R_i, A_iᵀ W_{ik}⟩_F .
/// ```
///
/// The Gauss–Newton residual vector stacks the norm-preserving symmetric
/// vectorization `svec(R_i) = (R_00, R_11, √2·R_01)` and its Jacobian stacks
/// `svec(∂R_i/∂θ_k)`, so `JᵀJ δ = −Jᵀr` is the exact Gauss–Newton system for
/// `E`; Levenberg damping plus the `det Dφ > δ` guard make every accepted
/// step a strict-descent diffeomorphism. Each `v_k` moves a single component,
/// so `W_{ik}` has one nonzero row and every row/mode contraction is a
/// handful of scalar ops.
///
/// # Honest refusals (`Ok(None)`, never a lossy or folded swap)
///
/// * degenerate chart: empty rows/basis, non-finite coordinates, or a
///   rank-deficient pullback metric (`det G_i ≤ 0`) anywhere;
/// * the optimizer finds no strict improvement over the identity flow (the
///   fitted chart is already minimum-defect within the flow family);
/// * every improving candidate violates the diffeomorphism guard;
/// * the basis cannot absorb `γ ∘ φ⁻¹` within
///   [`CHART_RECOMPOSITION_REL_TOL`] on the audit grid (shared gate with the
///   `d = 1` path).
pub fn torus_isometry_flow_reparameterization(
    evaluator: &dyn SaeBasisEvaluator,
    decoder: ArrayView2<'_, f64>,
    row_coords: ArrayView2<'_, f64>,
    period: f64,
) -> Result<Option<TorusIsometryFlowReparameterization>, String> {
    let n = row_coords.nrows();
    let m = decoder.nrows();
    let p = decoder.ncols();
    if row_coords.ncols() != 2 {
        return Err(format!(
            "torus_isometry_flow_reparameterization: expected (n, 2) row coordinates; got {:?}",
            row_coords.dim()
        ));
    }
    if n == 0 || m == 0 || p == 0 {
        return Ok(None);
    }
    for &t in row_coords.iter() {
        if !t.is_finite() {
            return Ok(None);
        }
    }

    // ── Fitted pullback metric G_i = J(t_i)ᵀ J(t_i) from the exact jet, then
    //    normalize against the flat reference g_ref = I (shared helpers) ──────
    let Some((g_rows, g_bar)) = extract_pullback_metric_d2(
        "torus_isometry_flow_reparameterization",
        evaluator,
        decoder,
        row_coords,
    )?
    else {
        return Ok(None);
    };
    let Some((ghat, ghat_norm_sq)) = flat_normalized_metric(&g_rows, g_bar) else {
        return Ok(None);
    };

    // ── Flow basis + per-row mode samples (W_{ik} and the displacements) ────
    let basis = TorusFlowBasis::new(period)?;
    let q = basis.dim();
    let mut row_modes: Vec<Vec<FlowModeSample>> = Vec::with_capacity(n);
    for row in 0..n {
        row_modes.push(basis.mode_samples([row_coords[[row, 0]], row_coords[[row, 1]]]));
    }
    // Flat reference ⇒ the whitened base is the identity at every row.
    let row_base = vec![[1.0_f64, 0.0, 0.0, 1.0]; n];

    // ── Damped Gauss–Newton on θ (shared exact core; derivation above) ──────
    let Some(minimization) = minimize_isometry_defect_flow(
        &row_modes,
        &row_base,
        &ghat,
        ghat_norm_sq,
        q,
        TORUS_FLOW_DIFFEO_MIN_DET,
        &|candidate: &[f64]| basis.certified_min_jacobian_det(candidate),
    ) else {
        return Ok(None);
    };
    let theta = minimization.theta;
    let defect_initial = minimization.defect_initial;
    let min_flow_jacobian_det = basis.certified_min_jacobian_det(&theta);
    if !(min_flow_jacobian_det > TORUS_FLOW_DIFFEO_MIN_DET) {
        // Unreachable through the guarded accept path; refuse defensively
        // rather than ever committing a folded chart.
        return Ok(None);
    }

    // ── Decoder transport on the Nyquist-oversampled audit grid ─────────────
    let axis_nodes = TORUS_TRANSPORT_MIN_NODES_PER_AXIS.max(3 * (m as f64).sqrt().ceil() as usize);
    let grid_rows = axis_nodes * axis_nodes;
    let mut grid = Array2::<f64>::zeros((grid_rows, 2));
    let mut new_grid = Array2::<f64>::zeros((grid_rows, 2));
    for i in 0..axis_nodes {
        for j in 0..axis_nodes {
            let idx = i * axis_nodes + j;
            let u = [
                period * i as f64 / axis_nodes as f64,
                period * j as f64 / axis_nodes as f64,
            ];
            grid[[idx, 0]] = u[0];
            grid[[idx, 1]] = u[1];
            let mapped = basis.map_point(&theta, u);
            new_grid[[idx, 0]] = mapped[0];
            new_grid[[idx, 1]] = mapped[1];
        }
    }
    let (grid_phi, grid_jet) = evaluator.evaluate(grid.view())?;
    if grid_phi.ncols() != m || grid_jet.dim() != (grid_rows, m, 2) {
        return Err(format!(
            "torus_isometry_flow_reparameterization: evaluator returned basis {:?} / jet {:?} on the audit grid; expected width {m}, latent_dim 2",
            grid_phi.dim(),
            grid_jet.dim()
        ));
    }
    let Some(recomposition) =
        recompose_decoder_exact_ls(evaluator, decoder, grid_phi.view(), new_grid.view())?
    else {
        return Ok(None);
    };

    // ── Canonical per-row coordinates t̃_i = φ_θ(t_i) ────────────────────────
    let mut new_row_coords = Array2::<f64>::zeros((n, 2));
    for row in 0..n {
        let mapped = basis.map_point(&theta, [row_coords[[row, 0]], row_coords[[row, 1]]]);
        new_row_coords[[row, 0]] = mapped[0];
        new_row_coords[[row, 1]] = mapped[1];
    }

    Ok(Some(TorusIsometryFlowReparameterization {
        new_row_coords,
        new_decoder: recomposition.new_decoder,
        decoder_transport: recomposition.transport,
        flow_theta: theta,
        defect_initial,
        defect_final: minimization.defect_final,
        profiled_metric_scale: minimization.profiled_scale,
        min_flow_jacobian_det,
        recomposition_residual: recomposition.recomposition_residual,
    }))
}

// ════════════════════════════════════════════════════════════════════════════
// #1019 stage 2 — d = 2 free/patch (Euclidean-patch) isometry-flow chart
// canonicalization
//
// Gauge-theory reading: a contractible patch carries no hairy-ball obstruction,
// so (unlike the sphere arm below) this slice is a genuine global minimization
// over an affine flow family, pinning the uniform-speed representative of the
// `Diff(patch)` orbit. The residual gauge left after the slice is the flat
// isometry group `O(2) ⋉ ℝ²` (rotation + reflection + translation) — reported
// via `PinnedByCanonicalization`, exactly the Prop-H stabilizer two patch fits
// must agree up to before they can be called the same bound feature.
// ════════════════════════════════════════════════════════════════════════════

/// Highest total polynomial degree of the truncated vector-field flow basis on
/// a `d = 2` Euclidean patch: monomial modes `u₀^a · u₁^b` with
/// `1 ≤ a+b ≤ PATCH_FLOW_MAX_DEGREE`.
///
/// The degree is **1** (affine flows) by exact-image-freezing necessity, not
/// timidity. The decoder lives in a polynomial basis of some fixed degree `d_b`
/// (the Duchon / EuclideanPatch monomial families), and the image is frozen by
/// re-expressing `γ̃ = γ ∘ φ⁻¹` in that SAME basis. Composing a degree-`d_b`
/// decoder image with a degree-`d_f` flow yields a degree-`d_b·d_f` polynomial,
/// which is back in the decoder basis **iff `d_f = 1`** (an affine flow leaves
/// the polynomial degree unchanged). A higher-degree flow would push `γ ∘ φ⁻¹`
/// out of the basis, the exact-LS recomposition gate
/// (`recompose_decoder_exact_ls`) would see a large image drift, and the
/// canonicalization would be honestly REFUSED — so a degree > 1 flow basis
/// could never actually commit a chart here. The affine family is exactly the
/// one that keeps the image-freezing exact, and it captures the dominant patch
/// chart dishonesty: an anisotropic / sheared / rotated affine reparameteriza-
/// tion (the `d = 2` analogue of the `d = 1` non-uniform-speed pathology).
///
/// The monomials are `(1, 0)` and `(0, 1)` per component, both components ⇒
/// `2 · 2 = 4` flow coefficients. The constant `(0, 0)` is EXCLUDED: a uniform
/// translation is an exact isometry of the flat patch, so it leaves the defect
/// invariant and would only insert a null direction into the Gauss–Newton
/// system. The 4 linear modes span the full `GL(2)` action on the chart
/// (rotation / shear / anisotropic scale); the rotation and global-scale
/// sub-directions are defect-null and the Levenberg damping absorbs them
/// harmlessly.
pub const PATCH_FLOW_MAX_DEGREE: usize = 1;

/// Diffeomorphism floor `δ` for the free-patch flow — identical contract to
/// [`TORUS_FLOW_DIFFEO_MIN_DET`]: a candidate with `det Dφ_θ ≤ δ` anywhere on
/// the check grid is rejected, so the optimizer can never walk through a fold.
pub const PATCH_FLOW_DIFFEO_MIN_DET: f64 = SAE_FLOW_DIFFEO_MIN_DET;

/// Per-axis node count of the free-patch diffeomorphism-guard check grid. With
/// the affine flow basis (`PATCH_FLOW_MAX_DEGREE = 1`) the Jacobian `Dφ_θ` is
/// CONSTANT over the patch, so `det Dφ_θ` is a single value and any grid — this
/// one, or a single node — is exact.
///
/// ⚠ The node count is therefore NOT a safety margin against a future degree
/// rise, and must not be read as one. TWO separate things rest on
/// `PATCH_FLOW_MAX_DEGREE = 1`, and a finer grid rescues only the first:
///
///  1. Grid exactness. At degree `> 1` the Jacobian varies over the patch and a
///     sampled minimum stops being the minimum. That is repairable, but not by
///     nodes alone — it needs a closed-form inter-node correction, exactly as
///     [`TorusFlowBasis::certified_min_jacobian_det`] does for the torus.
///  2. Injectivity, which NO grid can rescue. The guard exists to keep the
///     candidate inside `Diff(patch)`, and the step `det Dφ > 0` ⇒ injective is
///     a DEGREE-1 fact: an affine map with nonzero determinant is globally
///     invertible by linear algebra. At degree `> 1` it is simply false — a real
///     polynomial map can have `det Dφ > 0` everywhere and still be
///     non-injective (Pinchuk 1994). A degree-`> 1` patch flow needs a NEW
///     injectivity argument; more nodes are not one.
///
/// Raising `PATCH_FLOW_MAX_DEGREE` therefore invalidates this constant's
/// justification outright and must force a re-derivation of both points above.
/// (`PATCH_FLOW_MAX_DEGREE`'s own doc records a second, independent reason the
/// degree is pinned at 1: exact image-freezing needs the composed image to stay
/// inside the decoder's polynomial basis.)
///
/// The guard grid spans the normalized patch box `[-1, 1]²` slightly widened to
/// `[-1.1, 1.1]²` so a fold just outside the data hull is still refused.
pub const PATCH_FLOW_GUARD_NODES_PER_AXIS: usize = 48;

/// Minimum per-axis node count of the free-patch decoder-recomposition audit
/// grid (scaled up with the basis width like the torus path).
pub const PATCH_TRANSPORT_MIN_NODES_PER_AXIS: usize = 48;

/// Identity of one free-patch flow mode (for tests / diagnostics): which
/// coordinate component the vector field moves and its monomial exponents.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub struct PatchFlowModeKey {
    pub component: usize,
    pub exps: (usize, usize),
}

/// Truncated polynomial vector-field basis on a `d = 2` Euclidean patch,
/// `v_{c,(a,b)}(t) = e_c · u₀^a · u₁^b` where `u = (t − center) ⊙ inv_half ∈
/// [−1, 1]²` is the affinely-normalized patch coordinate (conditioning: the raw
/// `t` box can be far from `[−1, 1]`). Because the patch is **contractible**
/// there is no hairy-ball obstruction — these are smooth global vector fields,
/// so the flow `φ_θ(t) = t + Σ_k θ_k v_k(t)` parameterizes a neighborhood of
/// the identity in `Diff(patch)` exactly the way the Fourier basis does on the
/// torus, and any such map with `det Dφ > 0` on the patch is a diffeomorphism
/// onto its image. The reference metric is the flat `g_ref = I` (uniform-speed
/// coordinates), so the canonical chart is the minimum-anisotropy-defect one.
#[derive(Debug, Clone)]
pub struct FreePatchFlowBasis {
    center: [f64; 2],
    /// `2 / span` per axis: the Jacobian `∂u/∂t` of the normalization.
    inv_half: [f64; 2],
    /// Monomial exponents `(a, b)` with `1 ≤ a+b ≤ PATCH_FLOW_MAX_DEGREE`.
    exps: Vec<(usize, usize)>,
}

impl FreePatchFlowBasis {
    /// Build the basis for the patch box `[lo, hi]` per axis. Returns an error
    /// when an axis has collapsed (`hi ≤ lo`) — a patch with no extent in some
    /// direction has no honest flat chart.
    pub fn new(lo: [f64; 2], hi: [f64; 2]) -> Result<Self, String> {
        let mut center = [0.0_f64; 2];
        let mut inv_half = [0.0_f64; 2];
        for axis in 0..2 {
            let span = hi[axis] - lo[axis];
            let scale = lo[axis].abs().max(hi[axis].abs()).max(1.0);
            if !(span.is_finite() && span > 1.0e-12 * scale) {
                return Err(format!(
                    "FreePatchFlowBasis: patch axis {axis} has collapsed extent [{}, {}]",
                    lo[axis], hi[axis]
                ));
            }
            center[axis] = 0.5 * (lo[axis] + hi[axis]);
            inv_half[axis] = 2.0 / span;
        }
        let mut exps = Vec::new();
        for total in 1..=PATCH_FLOW_MAX_DEGREE {
            for a in (0..=total).rev() {
                let b = total - a;
                exps.push((a, b));
            }
        }
        Ok(Self {
            center,
            inv_half,
            exps,
        })
    }

    /// Number of flow coefficients `θ`: 2 components × #monomials.
    pub fn dim(&self) -> usize {
        2 * self.exps.len()
    }

    /// Normalized patch coordinate `u = (t − center) ⊙ inv_half`.
    fn normalize(&self, t: [f64; 2]) -> [f64; 2] {
        [
            (t[0] - self.center[0]) * self.inv_half[0],
            (t[1] - self.center[1]) * self.inv_half[1],
        ]
    }

    /// Sample every mode (value + gradient **in `t`**) at chart point `t`, in
    /// `θ` order. The monomial `f(u) = u₀^a · u₁^b` has `∂f/∂t_d =
    /// (∂f/∂u_d)·(∂u_d/∂t_d) = (∂f/∂u_d)·inv_half[d]` by the chain rule, so the
    /// returned gradient is already in the chart coordinate the flow Jacobian
    /// `Dφ = I + Σ θ_k ∂v_k/∂t` lives in.
    pub fn mode_samples(&self, t: [f64; 2]) -> Vec<FlowModeSample> {
        let u = self.normalize(t);
        let mut out = Vec::with_capacity(self.dim());
        for component in 0..2 {
            for &(a, b) in &self.exps {
                let value = pow_u(u[0], a) * pow_u(u[1], b);
                // ∂/∂u₀ (u₀^a u₁^b) = a·u₀^{a−1}·u₁^b ; ∂/∂u₁ = b·u₀^a·u₁^{b−1}.
                let du0 = if a == 0 {
                    0.0
                } else {
                    a as f64 * pow_u(u[0], a - 1) * pow_u(u[1], b)
                };
                let du1 = if b == 0 {
                    0.0
                } else {
                    b as f64 * pow_u(u[0], a) * pow_u(u[1], b - 1)
                };
                out.push(FlowModeSample {
                    component,
                    value,
                    grad: [du0 * self.inv_half[0], du1 * self.inv_half[1]],
                });
            }
        }
        out
    }

    /// `φ_θ(t) = t + Σ_k θ_k v_k(t)` (no wrap — the patch is not periodic).
    pub fn map_point(&self, theta: &[f64], t: [f64; 2]) -> [f64; 2] {
        assert_eq!(theta.len(), self.dim(), "FreePatchFlowBasis: theta length");
        let mut out = t;
        for (coef, sample) in theta.iter().zip(self.mode_samples(t)) {
            out[sample.component] += coef * sample.value;
        }
        out
    }

    /// Flow Jacobian `Dφ_θ(t) = I + Σ_k θ_k Dv_k(t)`, row-major.
    pub fn flow_jacobian(&self, theta: &[f64], t: [f64; 2]) -> [[f64; 2]; 2] {
        assert_eq!(theta.len(), self.dim(), "FreePatchFlowBasis: theta length");
        let mut jac = [[1.0, 0.0], [0.0, 1.0]];
        for (coef, sample) in theta.iter().zip(self.mode_samples(t)) {
            jac[sample.component][0] += coef * sample.grad[0];
            jac[sample.component][1] += coef * sample.grad[1];
        }
        jac
    }

    /// Minimum of `det Dφ_θ` over the widened normalized check grid `[−1.1,
    /// 1.1]²` (in `t` coordinates) — the diffeomorphism guard's decision
    /// quantity. The grid is built in `t` by inverting the normalization.
    pub fn min_jacobian_det_on_grid(&self, theta: &[f64]) -> f64 {
        let nodes = PATCH_FLOW_GUARD_NODES_PER_AXIS;
        let mut min_det = f64::INFINITY;
        for i in 0..nodes {
            for j in 0..nodes {
                // u ∈ [−1.1, 1.1] per axis; t = center + u / inv_half.
                let u0 = -1.1 + 2.2 * i as f64 / (nodes - 1) as f64;
                let u1 = -1.1 + 2.2 * j as f64 / (nodes - 1) as f64;
                let t = [
                    self.center[0] + u0 / self.inv_half[0],
                    self.center[1] + u1 / self.inv_half[1],
                ];
                let jac = self.flow_jacobian(theta, t);
                let det = jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0];
                min_det = min_det.min(det);
            }
        }
        min_det
    }
}

/// `u^k` for small non-negative integer `k` (avoids `f64::powi`'s branch on the
/// hot mode-sampling path; `k ≤ PATCH_FLOW_MAX_DEGREE`).
fn pow_u(u: f64, k: usize) -> f64 {
    let mut acc = 1.0_f64;
    for _ in 0..k {
        acc *= u;
    }
    acc
}

/// The exact, image-frozen minimum-isometry-defect flow reparameterization of
/// one `d = 2` **free/patch** (Euclidean-patch) atom.
#[derive(Debug, Clone)]
pub struct PatchIsometryFlowReparameterization {
    /// Canonical per-row coordinates `t̃_i = φ_θ(t_i)`, shape `(n, 2)`.
    pub new_row_coords: Array2<f64>,
    /// Recomposed decoder `B̃ = T · B`, shape `(m, p)` — exact LS refit of the
    /// original decoded image on the audit grid against the basis at the
    /// transported grid (so `γ̃ = γ ∘ φ⁻¹` without ever forming `φ⁻¹`).
    pub new_decoder: Array2<f64>,
    /// The `(m, m)` basis transport `T` with `Φ(φ(u)) · T ≈ Φ(u)` on the audit
    /// grid — the congruence object that transports the smoothness Gram.
    pub decoder_transport: Array2<f64>,
    /// Optimal flow coefficients `θ` (layout per
    /// `FreePatchFlowBasis::mode_layout`).
    pub flow_theta: Vec<f64>,
    /// Isometry defect `E(0)` of the fitted chart (identity flow).
    pub defect_initial: f64,
    /// Isometry defect `E(θ)` of the canonical chart. Strictly below
    /// `defect_initial` (the pass refuses no-improvement flows).
    pub defect_final: f64,
    /// The profiled global metric scale `c` at the optimum (the canonical
    /// chart's pullback metric is `≈ c·ḡ·I` — flat, uniform-speed).
    pub profiled_metric_scale: f64,
    /// `min det Dφ_θ` on the guard grid. Always `> PATCH_FLOW_DIFFEO_MIN_DET`
    /// when `Some(..)` is returned — a folded chart is refused upstream.
    pub min_flow_jacobian_det: f64,
    /// Max-abs recomposition error on the audit grid, relative to the image
    /// scale. Always `≤ CHART_RECOMPOSITION_REL_TOL` when `Some(..)`.
    pub recomposition_residual: f64,
}

/// Compute the minimum-isometry-defect flow reparameterization of a fitted
/// `d = 2` **free/patch** atom: the canonical per-row coordinates
/// `t̃_i = φ_θ(t_i)` and the exactly-recomposed decoder.
///
/// This is the #1019 unblocked free-chart case the issue charter calls out: for
/// a manifold patch (a contractible Euclidean-patch atom) a **global** truncated
/// flow basis DOES exist (no hairy-ball obstruction — [`FreePatchFlowBasis`]),
/// so the defect is genuinely MINIMIZED (not merely measured as on the sphere).
/// The reference metric is the flat `g_ref = I`, so the canonical chart is the
/// uniform-speed (minimum-anisotropy) one. Everything else — the scale-invariant
/// isometry defect `E(θ) = Σ_i ‖A_iᵀA_i − c·Ĝ_i‖²_F`, the analytic profiled scale
/// `c`, the exact Gauss–Newton, the `det Dφ > δ` diffeomorphism guard, and the
/// exact-LS decoder transport with the [`CHART_RECOMPOSITION_REL_TOL`] honesty
/// gate — is the SHARED machinery the torus path uses
/// (`minimize_isometry_defect_flow`, `recompose_decoder_exact_ls`); see
/// [`torus_isometry_flow_reparameterization`] for the full derivation.
///
/// The residual chart freedom after pinning is the finite isometry group of the
/// flat patch with the reference uniform metric: `O(2) ⋉ ℝ²` (rotation +
/// reflection + translation) — reported on the certificate as the
/// `PinnedByCanonicalization` residual gauge.
///
/// # Honest refusals (`Ok(None)`, never a lossy or folded swap)
///
/// * degenerate chart: empty rows/basis, non-finite coordinates, a collapsed
///   coordinate box on some axis, or a rank-deficient pullback metric anywhere;
/// * the optimizer finds no strict improvement over the identity flow (the
///   fitted chart is already minimum-defect within the flow family);
/// * every improving candidate violates the diffeomorphism guard;
/// * the basis cannot absorb `γ ∘ φ⁻¹` within
///   [`CHART_RECOMPOSITION_REL_TOL`] on the audit grid (shared gate).
pub fn patch_isometry_flow_reparameterization(
    evaluator: &dyn SaeBasisEvaluator,
    decoder: ArrayView2<'_, f64>,
    row_coords: ArrayView2<'_, f64>,
) -> Result<Option<PatchIsometryFlowReparameterization>, String> {
    let n = row_coords.nrows();
    let m = decoder.nrows();
    let p = decoder.ncols();
    if row_coords.ncols() != 2 {
        return Err(format!(
            "patch_isometry_flow_reparameterization: expected (n, 2) row coordinates; got {:?}",
            row_coords.dim()
        ));
    }
    if n == 0 || m == 0 || p == 0 {
        return Ok(None);
    }
    let mut lo = [f64::INFINITY; 2];
    let mut hi = [f64::NEG_INFINITY; 2];
    for row in 0..n {
        for axis in 0..2 {
            let t = row_coords[[row, axis]];
            if !t.is_finite() {
                return Ok(None);
            }
            lo[axis] = lo[axis].min(t);
            hi[axis] = hi[axis].max(t);
        }
    }

    // ── Fitted pullback metric, normalized against the flat reference I ──────
    let Some((g_rows, g_bar)) = extract_pullback_metric_d2(
        "patch_isometry_flow_reparameterization",
        evaluator,
        decoder,
        row_coords,
    )?
    else {
        return Ok(None);
    };
    let Some((ghat, ghat_norm_sq)) = flat_normalized_metric(&g_rows, g_bar) else {
        return Ok(None);
    };

    // ── Flow basis + per-row mode samples ───────────────────────────────────
    let basis = match FreePatchFlowBasis::new(lo, hi) {
        Ok(basis) => basis,
        // A collapsed patch axis: no honest flat chart — refuse, don't error.
        Err(_) => return Ok(None),
    };
    let q = basis.dim();
    let mut row_modes: Vec<Vec<FlowModeSample>> = Vec::with_capacity(n);
    for row in 0..n {
        row_modes.push(basis.mode_samples([row_coords[[row, 0]], row_coords[[row, 1]]]));
    }
    // Flat reference ⇒ the whitened base is the identity at every row.
    let row_base = vec![[1.0_f64, 0.0, 0.0, 1.0]; n];

    // ── Damped Gauss–Newton on θ (shared exact core) ────────────────────────
    let Some(minimization) = minimize_isometry_defect_flow(
        &row_modes,
        &row_base,
        &ghat,
        ghat_norm_sq,
        q,
        PATCH_FLOW_DIFFEO_MIN_DET,
        &|candidate: &[f64]| basis.min_jacobian_det_on_grid(candidate),
    ) else {
        return Ok(None);
    };
    let theta = minimization.theta;
    let defect_initial = minimization.defect_initial;
    let min_flow_jacobian_det = basis.min_jacobian_det_on_grid(&theta);
    if !(min_flow_jacobian_det > PATCH_FLOW_DIFFEO_MIN_DET) {
        return Ok(None);
    }

    // ── Decoder transport on the audit grid spanning the patch box ──────────
    let axis_nodes = PATCH_TRANSPORT_MIN_NODES_PER_AXIS.max(3 * (m as f64).sqrt().ceil() as usize);
    let grid_rows = axis_nodes * axis_nodes;
    let mut grid = Array2::<f64>::zeros((grid_rows, 2));
    let mut new_grid = Array2::<f64>::zeros((grid_rows, 2));
    for i in 0..axis_nodes {
        for j in 0..axis_nodes {
            let idx = i * axis_nodes + j;
            let u = [
                lo[0] + (hi[0] - lo[0]) * i as f64 / (axis_nodes - 1) as f64,
                lo[1] + (hi[1] - lo[1]) * j as f64 / (axis_nodes - 1) as f64,
            ];
            grid[[idx, 0]] = u[0];
            grid[[idx, 1]] = u[1];
            let mapped = basis.map_point(&theta, u);
            new_grid[[idx, 0]] = mapped[0];
            new_grid[[idx, 1]] = mapped[1];
        }
    }
    let (grid_phi, grid_jet) = evaluator.evaluate(grid.view())?;
    if grid_phi.ncols() != m || grid_jet.dim() != (grid_rows, m, 2) {
        return Err(format!(
            "patch_isometry_flow_reparameterization: evaluator returned basis {:?} / jet {:?} on the audit grid; expected width {m}, latent_dim 2",
            grid_phi.dim(),
            grid_jet.dim()
        ));
    }
    let Some(recomposition) =
        recompose_decoder_exact_ls(evaluator, decoder, grid_phi.view(), new_grid.view())?
    else {
        return Ok(None);
    };

    // ── Canonical per-row coordinates t̃_i = φ_θ(t_i) ────────────────────────
    let mut new_row_coords = Array2::<f64>::zeros((n, 2));
    for row in 0..n {
        let mapped = basis.map_point(&theta, [row_coords[[row, 0]], row_coords[[row, 1]]]);
        new_row_coords[[row, 0]] = mapped[0];
        new_row_coords[[row, 1]] = mapped[1];
    }

    Ok(Some(PatchIsometryFlowReparameterization {
        new_row_coords,
        new_decoder: recomposition.new_decoder,
        decoder_transport: recomposition.transport,
        flow_theta: theta,
        defect_initial,
        defect_final: minimization.defect_final,
        profiled_metric_scale: minimization.profiled_scale,
        min_flow_jacobian_det,
        recomposition_residual: recomposition.recomposition_residual,
    }))
}

// ════════════════════════════════════════════════════════════════════════════
// #1019 sphere arm — d = 2 sphere (S²) isometry-flow chart canonicalization
//
// Gauge-theory reading: on `S²` the hairy-ball theorem forbids a pole-free
// global flow basis, so the residual chart freedom after the gauge-invariant
// fit is not the finite isometry group directly — it is the full conformal
// (Möbius) group `PSL(2,ℂ)`, of which `O(3)` (the round sphere's isometries)
// is only a subgroup. This arm's slice therefore targets the non-isometric
// complement of that subgroup: the three conformal BOOSTS (gradients of the
// degree-1 harmonics). Minimizing the isometry defect over just the boosts
// breaks the conformal moduli down to `O(3)` — the residual stabilizer
// `PinnedByCanonicalization` reports — while leaving the true rotational
// isometries alone (they are already defect-null, i.e. already in the
// stabilizer, so including them would only add null directions to the solve).
// ════════════════════════════════════════════════════════════════════════════

/// Diffeomorphism floor `δ` for the sphere conformal-boost flow: a candidate
/// flow with `det Dφ_θ ≤ δ` anywhere on the data band is rejected, so the
/// optimizer never walks through a fold (identical contract to the torus /
/// patch floors).
pub const SPHERE_FLOW_DIFFEO_MIN_DET: f64 = SAE_FLOW_DIFFEO_MIN_DET;

/// Latitude band margin (radians) from each pole inside which the sphere
/// conformal-boost flow is well-conditioned. The off-meridian boost generators
/// carry a `1/cos lat` longitudinal component, so a chart whose data reaches
/// within this margin of a pole (`|lat| > π/2 − margin`) is honestly REFUSED
/// rather than pinned with a near-singular flow — the genuine residue of the
/// hairy-ball obstruction, scoped to exactly where it bites.
pub const SPHERE_FLOW_POLE_MARGIN: f64 = 0.20;

/// Identity of one sphere conformal-boost flow mode (for tests / diagnostics):
/// which round-sphere axis the boost points along.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum SphereBoostAxis {
    /// Zonal boost toward the poles, `cos(lat) ∂_lat` — pole-free in its single
    /// (latitude) component, the dominant sphere chart pathology.
    Z,
    /// Boost toward the `x = (lat 0, lon 0)` point.
    X,
    /// Boost toward the `y = (lat 0, lon π/2)` point.
    Y,
}

/// The three conformal-boost vector fields on the round `S²` in `(lat, lon)`
/// coordinates — the **non-isometric** part of the conformal group, i.e. the
/// gradient fields of the degree-1 spherical harmonics `z, x, y`.
///
/// # Why these three modes, and why they are the right (and complete) pin
///
/// After the gauge-invariant fit the sphere chart's residual freedom is the
/// conformal group of `(S², g_round)`, which is `PSL(2, ℂ)` (6-real-dim): the
/// three **rotation** Killing fields (the isometries `O(3)`) plus the three
/// **boosts**. The isometry defect `‖g_φ − g_ref‖²` is null on the rotations
/// (true isometries) and strictly positive on the boosts (a boost is a
/// non-isometric conformal map: it scales the metric by a non-constant
/// conformal factor `μ(x) ≠ 1`, compressing one hemisphere into a cap — the
/// `d = 2` analogue of the `d = 1` circle-into-`1.0`-rad pathology). So
/// minimizing the defect over the **boost** subspace breaks the conformal
/// moduli down to the residual isometry group `O(3)` — exactly the finite-dim
/// residual the certificate reports. The rotation modes are deliberately
/// EXCLUDED: they are defect-null and would only insert null directions into
/// the Gauss–Newton system (the same reason torus translations / patch
/// rotations are excluded from their flow bases).
///
/// In `(lat, lon)` the boost fields are (with `s = sin lat`, `c = cos lat`):
///
/// ```text
/// K_z = ( c,                0                 )      (gradient of z = sin lat)
/// K_x = ( s·cos lon,       −sin lon / c       )      (gradient of x)
/// K_y = ( s·sin lon,        cos lon / c       )      (gradient of y)
/// ```
///
/// `K_z` is globally smooth (pole-free); `K_x, K_y` carry the `1/c`
/// longitudinal pole the hairy-ball theorem guarantees somewhere. The flow is
/// the same Euler convention as the torus / patch families,
/// `φ_θ(t) = t + Σ_k θ_k v_k(t)`, and the band guard ([`SPHERE_FLOW_POLE_MARGIN`])
/// keeps every evaluated row away from the `1/c` singularity, so the three
/// fields are smooth and bounded on the data and the pin is well-posed.
#[derive(Debug, Clone)]
pub struct SphereBoostFlowBasis;

impl SphereBoostFlowBasis {
    /// The three boost modes are always present; the dimension is fixed at 3.
    pub fn dim(&self) -> usize {
        3
    }

    /// The displacement `v_k(t)` of each boost mode at chart point `t`.
    fn mode_displacements(t: [f64; 2]) -> [[f64; 2]; 3] {
        let (lat, lon) = (t[0], t[1]);
        let (s, c) = (lat.sin(), lat.cos());
        let (cl, sl) = (lon.cos(), lon.sin());
        [
            [c, 0.0],          // K_z
            [s * cl, -sl / c], // K_x
            [s * sl, cl / c],  // K_y
        ]
    }

    /// Jacobian `Dv_k(t)` (`[[∂v0/∂lat, ∂v0/∂lon], [∂v1/∂lat, ∂v1/∂lon]]`) of
    /// each boost mode at `t`. With `s = sin lat`, `c = cos lat`,
    /// `∂(1/c)/∂lat = s/c²`.
    fn mode_jacobians(t: [f64; 2]) -> [[[f64; 2]; 2]; 3] {
        let (lat, lon) = (t[0], t[1]);
        let (s, c) = (lat.sin(), lat.cos());
        let (cl, sl) = (lon.cos(), lon.sin());
        let c2 = c * c;
        [
            // K_z = (c, 0)
            [[-s, 0.0], [0.0, 0.0]],
            // K_x = (s·cl, −sl/c)
            [[c * cl, -s * sl], [-sl * s / c2, -cl / c]],
            // K_y = (s·sl, cl/c)
            [[c * sl, s * cl], [cl * s / c2, -sl / c]],
        ]
    }

    /// Min `det Dφ_θ` over a band-restricted check grid for the
    /// diffeomorphism guard. `φ_θ(t) = t + Σ_k θ_k v_k(t)`, so
    /// `Dφ_θ = I + Σ_k θ_k Dv_k`. The grid spans the data latitude band
    /// `[lat_lo, lat_hi]` (kept `SPHERE_FLOW_POLE_MARGIN` off each pole) and a
    /// full longitude turn.
    fn min_jacobian_det_on_band(theta: &[f64], lat_lo: f64, lat_hi: f64) -> f64 {
        let nodes = 48usize;
        let mut min_det = f64::INFINITY;
        for i in 0..nodes {
            let lat = lat_lo + (lat_hi - lat_lo) * i as f64 / (nodes - 1) as f64;
            for j in 0..nodes {
                let lon = -std::f64::consts::PI + std::f64::consts::TAU * j as f64 / nodes as f64;
                let jac = Self::mode_jacobians([lat, lon]);
                let mut a = [[1.0_f64, 0.0], [0.0, 1.0]];
                for (k, dv) in jac.iter().enumerate() {
                    a[0][0] += theta[k] * dv[0][0];
                    a[0][1] += theta[k] * dv[0][1];
                    a[1][0] += theta[k] * dv[1][0];
                    a[1][1] += theta[k] * dv[1][1];
                }
                let det = a[0][0] * a[1][1] - a[0][1] * a[1][0];
                min_det = min_det.min(det);
            }
        }
        min_det
    }

    /// Apply the flow `φ_θ(t) = t + Σ_k θ_k v_k(t)` to one chart point.
    fn map_point(theta: &[f64], t: [f64; 2]) -> [f64; 2] {
        let disp = Self::mode_displacements(t);
        let mut out = t;
        for (k, v) in disp.iter().enumerate() {
            out[0] += theta[k] * v[0];
            out[1] += theta[k] * v[1];
        }
        out
    }
}

/// The exact, image-frozen conformal-boost flow reparameterization of one
/// `d = 2` **sphere** atom: the canonical per-row `(lat, lon)` and the exactly
/// recomposed decoder.
#[derive(Debug, Clone)]
pub struct SphereIsometryFlowReparameterization {
    /// Canonical per-row coordinates `t̃_i = φ_θ(t_i)`, shape `(n, 2)`.
    pub new_row_coords: Array2<f64>,
    /// Recomposed decoder `B̃ = T · B`, shape `(m, p)` — the exact LS refit of
    /// the original decoded image on the audit grid against the basis at the
    /// transported grid (so `γ̃ = γ ∘ φ⁻¹` without ever forming `φ⁻¹`).
    pub new_decoder: Array2<f64>,
    /// The `(m, m)` basis transport `T` with `Φ(φ(u)) · T ≈ Φ(u)` on the audit
    /// grid — the congruence object the smoothness Gram is transported by.
    pub decoder_transport: Array2<f64>,
    /// Optimal boost coefficients `θ` (layout `[Z, X, Y]`).
    pub flow_theta: Vec<f64>,
    /// Round-sphere isometry defect `E(0)` of the fitted chart (identity flow).
    pub defect_initial: f64,
    /// Round-sphere isometry defect `E(θ)` of the canonical chart. Strictly
    /// below `defect_initial` (the pass refuses no-improvement flows).
    pub defect_final: f64,
    /// `min det Dφ_θ` on the band guard grid. Always `>
    /// SPHERE_FLOW_DIFFEO_MIN_DET` when `Some(..)` is returned.
    pub min_flow_jacobian_det: f64,
    /// Max-abs recomposition error on the audit grid, relative to the image
    /// scale. Always `≤ CHART_RECOMPOSITION_REL_TOL` when `Some(..)`.
    pub recomposition_residual: f64,
}

/// Compute the minimum-isometry-defect conformal-boost flow reparameterization
/// of a fitted `d = 2` sphere atom against the round-sphere reference metric.
///
/// # The defect, the reference-whitening reduction, and why the shared core fits
///
/// The canonical sphere chart is `t̃ = φ(t)` with `γ̃ = γ ∘ φ⁻¹` (image frozen),
/// and is round-isometric up to a global scale iff its pullback metric equals
/// `c · g_ref(lat) = c · diag(1, cos²lat)`. Measuring on the `φ` side (no
/// `Dφ⁻¹`), the chart is canonical iff `Dφ(t)ᵀ g_ref(φ(t)) Dφ(t) ≡ s · G(t)`
/// for some global `s` (`G = JᵀJ` the fitted pullback metric). The round
/// reference is diagonal, `g_ref(lat) = L(lat)ᵀ L(lat)` with
/// `L(lat) = diag(1, cos lat)`, so writing the **whitened** flow Jacobian
/// `Ã_i = L(φ(t_i)) · Dφ(t_i)`,
///
/// ```text
/// Dφ_iᵀ g_ref,i Dφ_i = Ã_iᵀ Ã_i ,
/// ```
///
/// the residual `R_i = Ã_iᵀ Ã_i − s · Ĝ_i` is the same symmetric-metric
/// residual used by the flat-reference arms. At the identity flow `Dφ = I`,
/// `Ã_i = L(lat_i) = diag(1, cos lat_i)`. Away from identity, both factors in
/// `Ã_i = L(φ_θ(t_i))Dφ_θ(t_i)` are evaluated at the current `θ`; the sphere
/// therefore assembles its own exact residual and analytic Gauss--Newton
/// Jacobian (the boost modes move both coordinate components, unlike the
/// single-component `FlowModeSample` contract).
///
/// For boost `k`, write `V_ik = Dv_k(t_i)`, `u_ik = v_k(t_i)_lat`, moved
/// latitude `ℓ̃_i = φ_θ(t_i)_lat`, and `D_i = Dφ_θ(t_i)`. The complete chain is
///
/// ```text
/// ∂D_i/∂θ_k = V_ik,
/// ∂L_i/∂θ_k = diag(0, −sin(ℓ̃_i) u_ik),
/// ∂Ã_i/∂θ_k = (∂L_i/∂θ_k) D_i + L_i V_ik,
/// ∂(Ã_iᵀÃ_i)/∂θ_k = (∂Ã_i/∂θ_k)ᵀÃ_i + Ã_iᵀ(∂Ã_i/∂θ_k).
/// ```
///
/// The global scale `s(θ)` is profiled exactly. Its derivative projects each
/// metric differential off the stacked `Ĝ` direction, so the residual Jacobian
/// is `J = (I − ĝ ĝᵀ / ‖ĝ‖²) D_θ svec(ÃᵀÃ)`. Consequently the moved-latitude
/// whitening, metric product, and profiled scale all participate analytically
/// in the LM system; only the tests use finite differences as an audit oracle.
///
/// # Honest refusals (`Ok(None)`)
///
/// * degenerate chart (empty / non-finite, rank-deficient pullback metric);
/// * data within [`SPHERE_FLOW_POLE_MARGIN`] of a pole (the `1/cos lat` boost
///   singularity — the scoped hairy-ball residue);
/// * the fitted chart is already minimum-defect (no strict improvement);
/// * every improving candidate violates the diffeomorphism guard;
/// * the basis cannot absorb `γ ∘ φ⁻¹` within [`CHART_RECOMPOSITION_REL_TOL`].
pub fn sphere_isometry_flow_reparameterization(
    evaluator: &dyn SaeBasisEvaluator,
    decoder: ArrayView2<'_, f64>,
    row_coords: ArrayView2<'_, f64>,
) -> Result<Option<SphereIsometryFlowReparameterization>, String> {
    let n = row_coords.nrows();
    let m = decoder.nrows();
    let p = decoder.ncols();
    if row_coords.ncols() != 2 {
        return Err(format!(
            "sphere_isometry_flow_reparameterization: expected (n, 2) row coordinates; got {:?}",
            row_coords.dim()
        ));
    }
    if n == 0 || m == 0 || p == 0 {
        return Ok(None);
    }
    for &t in row_coords.iter() {
        if !t.is_finite() {
            return Ok(None);
        }
    }

    // ── Pole-band guard: refuse charts reaching into the 1/cos lat singularity ─
    let mut lat_lo = f64::INFINITY;
    let mut lat_hi = f64::NEG_INFINITY;
    for row in 0..n {
        let lat = row_coords[[row, 0]];
        lat_lo = lat_lo.min(lat);
        lat_hi = lat_hi.max(lat);
    }
    let pole = std::f64::consts::FRAC_PI_2 - SPHERE_FLOW_POLE_MARGIN;
    if !(lat_lo > -pole && lat_hi < pole) {
        return Ok(None);
    }

    // ── Fitted pullback metric G_i = J(t_i)ᵀ J(t_i) (raw, NOT flat-normalized:
    //    the sphere target is g_ref, not I), with the geometric-mean scale ḡ ───
    let Some((g_rows, g_bar)) = extract_pullback_metric_d2(
        "sphere_isometry_flow_reparameterization",
        evaluator,
        decoder,
        row_coords,
    )?
    else {
        return Ok(None);
    };
    // ĝ_i = G_i / ḡ; the round reference target enters per row via the whitened
    // base below, so the GN `ghat` here is the scale-normalized fitted metric.
    let mut ghat: Vec<[f64; 3]> = Vec::with_capacity(n);
    let mut ghat_norm_sq = 0.0_f64;
    for g in &g_rows {
        let h = [g[0] / g_bar, g[1] / g_bar, g[2] / g_bar];
        ghat_norm_sq += h[0] * h[0] + h[1] * h[1] + 2.0 * h[2] * h[2];
        ghat.push(h);
    }
    if !(ghat_norm_sq.is_finite() && ghat_norm_sq > 0.0) {
        return Ok(None);
    }

    // ── Damped Gauss–Newton over the 3 boost coefficients ───────────────────
    let Some(minimization) =
        sphere_minimize_boost_defect(&ghat, ghat_norm_sq, row_coords, lat_lo, lat_hi)
    else {
        return Ok(None);
    };
    let theta = minimization.theta;
    let defect_initial = minimization.defect_initial;
    let min_flow_jacobian_det =
        SphereBoostFlowBasis::min_jacobian_det_on_band(&theta, lat_lo, lat_hi);
    if !(min_flow_jacobian_det > SPHERE_FLOW_DIFFEO_MIN_DET) {
        return Ok(None);
    }

    // ── Decoder transport on an audit grid spanning the data band ───────────
    let axis_nodes = TORUS_TRANSPORT_MIN_NODES_PER_AXIS.max(3 * (m as f64).sqrt().ceil() as usize);
    let grid_rows = axis_nodes * axis_nodes;
    let mut grid = Array2::<f64>::zeros((grid_rows, 2));
    let mut new_grid = Array2::<f64>::zeros((grid_rows, 2));
    for i in 0..axis_nodes {
        for j in 0..axis_nodes {
            let idx = i * axis_nodes + j;
            let lat = lat_lo + (lat_hi - lat_lo) * i as f64 / (axis_nodes - 1) as f64;
            let lon = -std::f64::consts::PI + std::f64::consts::TAU * j as f64 / axis_nodes as f64;
            grid[[idx, 0]] = lat;
            grid[[idx, 1]] = lon;
            let mapped = SphereBoostFlowBasis::map_point(&theta, [lat, lon]);
            new_grid[[idx, 0]] = mapped[0];
            new_grid[[idx, 1]] = mapped[1];
        }
    }
    let (grid_phi, grid_jet) = evaluator.evaluate(grid.view())?;
    if grid_phi.ncols() != m || grid_jet.dim() != (grid_rows, m, 2) {
        return Err(format!(
            "sphere_isometry_flow_reparameterization: evaluator returned basis {:?} / jet {:?} on the audit grid; expected width {m}, latent_dim 2",
            grid_phi.dim(),
            grid_jet.dim()
        ));
    }
    let Some(recomposition) =
        recompose_decoder_exact_ls(evaluator, decoder, grid_phi.view(), new_grid.view())?
    else {
        return Ok(None);
    };

    // ── Canonical per-row coordinates t̃_i = φ_θ(t_i) ────────────────────────
    let mut new_row_coords = Array2::<f64>::zeros((n, 2));
    for row in 0..n {
        let mapped =
            SphereBoostFlowBasis::map_point(&theta, [row_coords[[row, 0]], row_coords[[row, 1]]]);
        new_row_coords[[row, 0]] = mapped[0];
        new_row_coords[[row, 1]] = mapped[1];
    }

    Ok(Some(SphereIsometryFlowReparameterization {
        new_row_coords,
        new_decoder: recomposition.new_decoder,
        decoder_transport: recomposition.transport,
        flow_theta: theta,
        defect_initial,
        defect_final: minimization.defect_final,
        min_flow_jacobian_det,
        recomposition_residual: recomposition.recomposition_residual,
    }))
}

/// Outcome of the sphere boost-defect minimizer.
struct SphereFlowMinimization {
    theta: Vec<f64>,
    defect_initial: f64,
    defect_final: f64,
}

/// Whitened sphere-flow Jacobian `Ã = L(φ_θ(t))Dφ_θ(t)` and its exact
/// derivatives with respect to the `[Z, X, Y]` boost coefficients.
///
/// The derivative includes both the linear flow-Jacobian term and the
/// moved-latitude whitening term:
/// `∂Ã_k = L Dv_k + diag(0, -sin(lat̃) v_k,lat) Dφ`.
fn sphere_whitened_boost_row(theta: &[f64], t: [f64; 2]) -> Option<([f64; 4], [[f64; 4]; 3])> {
    let q = SphereBoostFlowBasis.dim();
    if theta.len() != q {
        return None;
    }
    let displacements = SphereBoostFlowBasis::mode_displacements(t);
    let mode_jacobians = SphereBoostFlowBasis::mode_jacobians(t);
    let mut dphi = [[1.0_f64, 0.0], [0.0, 1.0]];
    let mut moved_lat = t[0];
    for k in 0..q {
        moved_lat += theta[k] * displacements[k][0];
        let dv = mode_jacobians[k];
        dphi[0][0] += theta[k] * dv[0][0];
        dphi[0][1] += theta[k] * dv[0][1];
        dphi[1][0] += theta[k] * dv[1][0];
        dphi[1][1] += theta[k] * dv[1][1];
    }

    let (sin_lat_new, cos_lat_new) = moved_lat.sin_cos();
    // Mirror the pole floor in `sphere_chart_isometry_defect`: at smaller
    // `|cos(lat̃)|`, the longitudinal whitening row has no reliable metric
    // content and the chart is honestly outside this coordinate patch.
    const SPHERE_EVAL_COS_FLOOR: f64 = 1.0e-6;
    if !(cos_lat_new.is_finite() && cos_lat_new > SPHERE_EVAL_COS_FLOOR) {
        return None;
    }

    let a = [
        dphi[0][0],
        dphi[0][1],
        cos_lat_new * dphi[1][0],
        cos_lat_new * dphi[1][1],
    ];
    let mut da = [[0.0_f64; 4]; 3];
    for k in 0..q {
        let dv = mode_jacobians[k];
        let dcos = -sin_lat_new * displacements[k][0];
        da[k] = [
            dv[0][0],
            dv[0][1],
            dcos * dphi[1][0] + cos_lat_new * dv[1][0],
            dcos * dphi[1][1] + cos_lat_new * dv[1][1],
        ];
    }
    if !a.iter().all(|value| value.is_finite())
        || !da.iter().flatten().all(|value| value.is_finite())
    {
        return None;
    }
    Some((a, da))
}

/// Per-row whitened flow Jacobian `Ã_i = L(φ(t_i)) · Dφ_θ(t_i)` (row-major
/// `[a00, a01, a10, a11]`) and the round-sphere defect at `θ`.
///
/// `Dφ_θ(t) = I + Σ_k θ_k Dv_k(t)`; whitening by `L(lat̃) = diag(1, cos lat̃)`
/// at the moved latitude `lat̃ = φ_θ(t)[0]` realizes the reference metric (see
/// [`sphere_isometry_flow_reparameterization`] for the reduction). Returns
/// `None` when the profiled scale degenerates or any whitened latitude leaves
/// the valid band (`cos lat̃ ≤ 0`).
fn sphere_eval_boost_defect(
    theta: &[f64],
    row_coords: ArrayView2<'_, f64>,
    ghat: &[[f64; 3]],
    ghat_norm_sq: f64,
) -> Option<FlowObjectiveState> {
    let n = row_coords.nrows();
    let mut a_rows: Vec<[f64; 4]> = Vec::with_capacity(n);
    let mut cross = 0.0_f64;
    for row in 0..n {
        let t = [row_coords[[row, 0]], row_coords[[row, 1]]];
        let (a, _da) = sphere_whitened_boost_row(theta, t)?;
        a_rows.push(a);
    }
    for (a, g) in a_rows.iter().zip(ghat.iter()) {
        let m00 = a[0] * a[0] + a[2] * a[2];
        let m11 = a[1] * a[1] + a[3] * a[3];
        let m01 = a[0] * a[1] + a[2] * a[3];
        cross += m00 * g[0] + m11 * g[1] + 2.0 * m01 * g[2];
    }
    let scale = cross / ghat_norm_sq;
    if !(scale.is_finite() && scale > 0.0) {
        return None;
    }
    let mut defect = 0.0_f64;
    for (a, g) in a_rows.iter().zip(ghat.iter()) {
        let m00 = a[0] * a[0] + a[2] * a[2];
        let m11 = a[1] * a[1] + a[3] * a[3];
        let m01 = a[0] * a[1] + a[2] * a[3];
        let r00 = m00 - scale * g[0];
        let r11 = m11 - scale * g[1];
        let r01 = m01 - scale * g[2];
        defect += r00 * r00 + r11 * r11 + 2.0 * r01 * r01;
    }
    if !defect.is_finite() {
        return None;
    }
    Some(FlowObjectiveState {
        defect,
        scale,
        a_rows,
    })
}

/// Exact profiled residual and analytic Jacobian for the sphere boost flow.
///
/// If `m_i = svec(Ã_iᵀÃ_i)` and `g_i = svec(Ĝ_i)`, the profiled scale is
/// `s = (Σ g_iᵀm_i) / (Σ g_iᵀg_i)`. Therefore
///
/// ```text
/// ds/dθ_k = Σ_i g_iᵀ dm_i/dθ_k / ‖g‖²,
/// dr_i/dθ_k = dm_i/dθ_k − (ds/dθ_k) g_i.
/// ```
///
/// The subtraction is load-bearing even though it cancels from `Jᵀr` by the
/// envelope condition: it changes `JᵀJ`, and hence the LM chart selected by the
/// profiled objective.
fn sphere_boost_residual_jacobian(
    theta: &[f64],
    row_coords: ArrayView2<'_, f64>,
    ghat: &[[f64; 3]],
    ghat_norm_sq: f64,
    state: &FlowObjectiveState,
) -> Option<(Array2<f64>, Array2<f64>)> {
    let n = row_coords.nrows();
    let q = SphereBoostFlowBasis.dim();
    if theta.len() != q
        || row_coords.ncols() != 2
        || ghat.len() != n
        || state.a_rows.len() != n
        || !(ghat_norm_sq.is_finite() && ghat_norm_sq > 0.0 && state.scale.is_finite())
    {
        return None;
    }

    let sqrt2 = std::f64::consts::SQRT_2;
    let mut residual = Array2::<f64>::zeros((3 * n, 1));
    let mut jacobian = Array2::<f64>::zeros((3 * n, q));
    let mut scale_gradient = vec![0.0_f64; q];

    for row in 0..n {
        let t = [row_coords[[row, 0]], row_coords[[row, 1]]];
        let (a, da) = sphere_whitened_boost_row(theta, t)?;
        let g = ghat[row];
        let m00 = a[0] * a[0] + a[2] * a[2];
        let m11 = a[1] * a[1] + a[3] * a[3];
        let m01 = a[0] * a[1] + a[2] * a[3];
        residual[[3 * row, 0]] = m00 - state.scale * g[0];
        residual[[3 * row + 1, 0]] = m11 - state.scale * g[1];
        residual[[3 * row + 2, 0]] = sqrt2 * (m01 - state.scale * g[2]);

        for k in 0..q {
            let b = da[k];
            let dm00 = 2.0 * (a[0] * b[0] + a[2] * b[2]);
            let dm11 = 2.0 * (a[1] * b[1] + a[3] * b[3]);
            let dm01 = b[0] * a[1] + a[0] * b[1] + b[2] * a[3] + a[2] * b[3];
            jacobian[[3 * row, k]] = dm00;
            jacobian[[3 * row + 1, k]] = dm11;
            jacobian[[3 * row + 2, k]] = sqrt2 * dm01;
            scale_gradient[k] += dm00 * g[0] + dm11 * g[1] + 2.0 * dm01 * g[2];
        }
    }

    for k in 0..q {
        scale_gradient[k] /= ghat_norm_sq;
        for row in 0..n {
            let g = ghat[row];
            jacobian[[3 * row, k]] -= scale_gradient[k] * g[0];
            jacobian[[3 * row + 1, k]] -= scale_gradient[k] * g[1];
            jacobian[[3 * row + 2, k]] -= sqrt2 * scale_gradient[k] * g[2];
        }
    }

    if residual.iter().all(|value| value.is_finite())
        && jacobian.iter().all(|value| value.is_finite())
    {
        Some((residual, jacobian))
    } else {
        None
    }
}

/// Damped Gauss--Newton on the three sphere conformal-boost coefficients. The
/// residual `svec(Ã_iᵀÃ_i − s·Ĝ_i)` and its `θ`-Jacobian are analytic through
/// the boost, moved-latitude whitening, metric product, and profiled scale (see
/// [`sphere_boost_residual_jacobian`]). The Levenberg damping, strict-descent
/// gate, and diffeomorphism accept test match the torus / patch core. Starts at
/// `θ = 0` (`Dφ = I`); only fold-free strict-descent candidates are accepted.
/// Returns `None` (honest skip) when the identity chart is already
/// round-isometric or no strict improvement is reachable.
fn sphere_minimize_boost_defect(
    ghat: &[[f64; 3]],
    ghat_norm_sq: f64,
    row_coords: ArrayView2<'_, f64>,
    lat_lo: f64,
    lat_hi: f64,
) -> Option<SphereFlowMinimization> {
    let q = SphereBoostFlowBasis.dim();
    let mut theta = vec![0.0_f64; q];
    let mut state = sphere_eval_boost_defect(&theta, row_coords, ghat, ghat_norm_sq)?;
    let defect_initial = state.defect;
    if !(defect_initial > 0.0) {
        return None;
    }
    let mut lambda = 1.0e-4_f64;
    let mut any_accepted = false;
    for iteration in 0..TORUS_FLOW_GN_MAX_ITERS {
        if iteration + 1 == TORUS_FLOW_GN_MAX_ITERS {
            break;
        }
        let (rcol, jmat) =
            sphere_boost_residual_jacobian(&theta, row_coords, ghat, ghat_norm_sq, &state)?;
        let jtj = fast_ata(&jmat);
        let jtr = fast_atb(&jmat, &rcol);

        // Same Levenberg-damped strict-descent + diffeomorphism accept test as
        // the torus core; the sphere-specific fold guard and defect evaluation
        // enter only through the `eval_candidate` closure.
        let step = lm_damped_accept_sweep(
            &jtj,
            &jtr,
            q,
            &mut lambda,
            &mut theta,
            &mut state,
            |candidate| {
                if SphereBoostFlowBasis::min_jacobian_det_on_band(candidate, lat_lo, lat_hi)
                    <= SPHERE_FLOW_DIFFEO_MIN_DET
                {
                    None
                } else {
                    sphere_eval_boost_defect(candidate, row_coords, ghat, ghat_norm_sq)
                }
            },
        );
        any_accepted |= step.accepted;
        if !step.accepted || step.converged {
            break;
        }
        let theta_norm_sq: f64 = theta.iter().map(|v| v * v).sum();
        if step.step_norm_sq <= LM_STEP_STALL_REL_FLOOR * (1.0 + theta_norm_sq) {
            break;
        }
    }
    if !any_accepted || !(state.defect < defect_initial) {
        return None;
    }
    Some(SphereFlowMinimization {
        theta,
        defect_initial,
        defect_final: state.defect,
    })
}

/// Total fitted turning `Θ = ∫ κ ds` of a `d = 1` atom's decoded curve (#1026).
///
/// # Why this is the discriminating measurement
///
/// The hybrid-vs-shatter question — does a curved atom genuinely earn its
/// curvature, or is it just more linear directions in disguise — is answered by
/// pairing each atom's EV contribution with its fitted **turning** `Θ`
/// (integrated curvature). A linear SAE shatters a curved feature of total
/// turning `Θ` into `N(ε) ≈ Θ/(2√(2ε))` rank-1 atoms at relative error `ε`
/// (radius cancels — relative error is scale-free), so the curved win is
/// concentrated on high-`Θ` features and vanishes as `Θ → 0`. A near-linear
/// atom (`Θ ≈ 0`) contributing EV is a linear direction wearing a curved basis;
/// a high-`Θ` atom contributing EV is a genuine curved family. Reporting EV-per-
/// atom vs `Θ` (not EV vs K) directly shows which.
///
/// # The functional
///
/// For the decoded curve `γ(t) = Φ(t)·B` the unsigned total curvature is
///
/// ```text
/// Θ = ∫ κ(t) ‖γ'(t)‖ dt ,   κ = ‖γ'(t) ∧ γ''(t)‖ / ‖γ'(t)‖³ ,
/// ```
///
/// so `Θ = ∫ ‖γ' ∧ γ''‖ / ‖γ'‖² dt`, where the wedge norm in `ℝ^p` is the
/// parallelogram area `‖γ'∧γ''‖ = √(‖γ'‖²‖γ''‖² − ⟨γ',γ''⟩²)` (Lagrange). `Θ`
/// is reparameterization-invariant (it is an integral of `κ ds`), so it is the
/// honest chart-free geometric content. Integrated by Simpson's rule over a
/// uniform grid spanning the fitted coordinate range from the exact `(γ', γ'')`
/// jets. Units: radians of total turning (a full hue circle ≈ `2π`).
///
/// Returns `None` on a degenerate atom (no rows, no second jet, a collapsed
/// coordinate range, or a non-finite integrand) — an honest refusal, never a
/// fabricated number.
pub fn d1_atom_fitted_turning(
    evaluator: &dyn SaeBasisEvaluator,
    decoder: ArrayView2<'_, f64>,
    row_coords: ArrayView1<'_, f64>,
) -> Result<Option<f64>, String> {
    let m = decoder.nrows();
    let p = decoder.ncols();
    if m == 0 || p == 0 || row_coords.is_empty() {
        return Ok(None);
    }
    let mut lo = f64::INFINITY;
    let mut hi = f64::NEG_INFINITY;
    for &t in row_coords.iter() {
        if !t.is_finite() {
            return Ok(None);
        }
        lo = lo.min(t);
        hi = hi.max(t);
    }
    if !(hi > lo) {
        return Ok(None);
    }
    let cells = TURNING_QUADRATURE_CELLS;
    let nodes = 2 * cells + 1; // node, midpoint, node per cell (Simpson)
    let h = (hi - lo) / (nodes - 1) as f64;
    let mut grid = Array2::<f64>::zeros((nodes, 1));
    for (i, mut row) in grid.outer_iter_mut().enumerate() {
        row[0] = lo + h * i as f64;
    }
    let (_phi, jet) = evaluator.evaluate(grid.view())?;
    if jet.dim() != (nodes, m, 1) {
        return Err(format!(
            "d1_atom_fitted_turning: evaluator returned jet {:?}; expected ({nodes}, {m}, 1)",
            jet.dim()
        ));
    }
    let Some(hess_result) = evaluator.second_jet_dyn(grid.view()) else {
        // No analytic second jet for this basis family → no honest turning.
        return Ok(None);
    };
    let hess = hess_result?;
    if hess.dim() != (nodes, m, 1, 1) {
        return Err(format!(
            "d1_atom_fitted_turning: second_jet returned {:?}; expected ({nodes}, {m}, 1, 1)",
            hess.dim()
        ));
    }
    // Per-node curvature integrand f(t) = ‖γ'∧γ''‖/‖γ'‖² (= κ·‖γ'‖).
    let mut integrand = vec![0.0_f64; nodes];
    // Distinguish a GLOBALLY constant image (every node stationary → a degenerate
    // point with total turning exactly 0) from a curve that is stationary at only
    // SOME nodes (a genuine cusp whose turning is ill-defined there, refused).
    let mut any_moving = false;
    let mut any_stationary = false;
    let mut g1 = vec![0.0_f64; p];
    let mut g2 = vec![0.0_f64; p];
    for node in 0..nodes {
        for slot in g1.iter_mut() {
            *slot = 0.0;
        }
        for slot in g2.iter_mut() {
            *slot = 0.0;
        }
        for bm in 0..m {
            let d1 = jet[[node, bm, 0]];
            let d2 = hess[[node, bm, 0, 0]];
            if d1 == 0.0 && d2 == 0.0 {
                continue;
            }
            for j in 0..p {
                let b = decoder[[bm, j]];
                g1[j] += d1 * b;
                g2[j] += d2 * b;
            }
        }
        let mut n1 = 0.0_f64; // ‖γ'‖²
        let mut n2 = 0.0_f64; // ‖γ''‖²
        let mut dot = 0.0_f64; // ⟨γ',γ''⟩
        for j in 0..p {
            n1 += g1[j] * g1[j];
            n2 += g2[j] * g2[j];
            dot += g1[j] * g2[j];
        }
        if !(n1 > 0.0) {
            // Zero speed at this node: the arc-length measure `ds = ‖γ'‖dt`
            // vanishes here, so this node contributes zero turning regardless of
            // κ. Record it and continue; the all-stationary (constant image) vs
            // mixed (cusp) cases are resolved after the loop.
            any_stationary = true;
            integrand[node] = 0.0;
            continue;
        }
        any_moving = true;
        // Wedge norm² = ‖γ'‖²‖γ''‖² − ⟨γ',γ''⟩² (Lagrange identity); clamp tiny
        // negative round-off to 0 before the sqrt.
        let raw_wedge_sq = n1 * n2 - dot * dot;
        let roundoff_floor = 64.0 * f64::EPSILON * (n1 * n2).abs().max(dot.abs() * dot.abs());
        let wedge_sq = if raw_wedge_sq <= roundoff_floor {
            0.0
        } else {
            raw_wedge_sq
        };
        integrand[node] = wedge_sq.sqrt() / n1;
        if !integrand[node].is_finite() {
            return Ok(None);
        }
    }
    if !any_moving {
        // The decoded image never moves over the coordinate span: a degenerate
        // single point (e.g. a `d = 1` atom straightened to its constant DC
        // component). A point has no arc to turn through, so its total turning is
        // exactly 0 — the ultimate linear-tail signature — not "undefined".
        return Ok(Some(0.0));
    }
    if any_stationary {
        // Stationary at SOME nodes but moving at others: a cusp where the turning
        // integrand is genuinely ill-defined. Refuse rather than under-count the
        // sharp turn (the historical conservative behavior for a partial cusp).
        return Ok(None);
    }
    // Composite Simpson over `cells` cells: each cell [2i, 2i+1, 2i+2] gets
    // (h/3)(f0 + 4 f_mid + f1).
    let mut theta = 0.0_f64;
    for cell in 0..cells {
        let f0 = integrand[2 * cell];
        let fm = integrand[2 * cell + 1];
        let f1 = integrand[2 * cell + 2];
        theta += h / 3.0 * (f0 + 4.0 * fm + f1);
    }
    if !(theta.is_finite() && theta >= 0.0) {
        return Ok(None);
    }
    Ok(Some(theta))
}

#[cfg(test)]
mod patch_flow_tests {
    use super::*;
    use ndarray::{Array2, Array3, ArrayView2};

    /// Mock free-patch evaluator with the affine basis `Φ(t) = [1, t₀, t₁]`
    /// (`m = 3`) and its exact jet. The decoder plants the affine decoded image
    /// `γ(t) = M·t` (a constant anisotropic warp), so the fitted pullback metric
    /// is the constant `G = MᵀM` everywhere — an anisotropically warped flat
    /// chart. The affine basis is closed under affine reparameterization, so the
    /// linear flow modes can transport the image exactly (the recomposition gate
    /// passes), and the minimizer recovers `Dφ ≈ M⁻¹` (up to a flat O(2)
    /// isometry), driving the anisotropy defect to ≈ 0.
    #[derive(Debug)]
    struct MockPatchEvaluator;

    impl SaeBasisEvaluator for MockPatchEvaluator {
        fn evaluate(
            &self,
            coords: ArrayView2<'_, f64>,
        ) -> Result<(Array2<f64>, Array3<f64>), String> {
            let n = coords.nrows();
            let mut phi = Array2::<f64>::zeros((n, 3));
            let mut jet = Array3::<f64>::zeros((n, 3, 2));
            for row in 0..n {
                let t0 = coords[[row, 0]];
                let t1 = coords[[row, 1]];
                phi[[row, 0]] = 1.0;
                phi[[row, 1]] = t0;
                phi[[row, 2]] = t1;
                // ∂Φ/∂t₀ = [0, 1, 0]; ∂Φ/∂t₁ = [0, 0, 1].
                jet[[row, 1, 0]] = 1.0;
                jet[[row, 2, 1]] = 1.0;
            }
            Ok((phi, jet))
        }

        fn second_jet_dyn(
            &self,
            coords: ArrayView2<'_, f64>,
        ) -> Option<Result<ndarray::Array4<f64>, String>> {
            // Affine basis: second jet is exactly zero. Return it honestly
            // (shape (n, 3, 2, 2)) after validating the coordinate width.
            if coords.ncols() != 2 {
                return Some(Err(format!(
                    "MockPatchEvaluator::second_jet_dyn: expected 2 cols, got {}",
                    coords.ncols()
                )));
            }
            Some(Ok(ndarray::Array4::<f64>::zeros((coords.nrows(), 3, 2, 2))))
        }

        fn third_jet_dyn(
            &self,
            coords: ArrayView2<'_, f64>,
        ) -> Option<Result<ndarray::Array5<f64>, String>> {
            if coords.ncols() != 2 {
                return Some(Err(format!(
                    "MockPatchEvaluator::third_jet_dyn: expected 2 cols, got {}",
                    coords.ncols()
                )));
            }
            Some(Ok(ndarray::Array5::<f64>::zeros((
                coords.nrows(),
                3,
                2,
                2,
                2,
            ))))
        }
    }

    /// Decoder for `γ(t) = M·t`: basis `[1, t₀, t₁]` ↦ outputs `(γ₀, γ₁)`.
    /// Row 0 (const) = 0; row 1 (t₀) = (M₀₀, M₁₀); row 2 (t₁) = (M₀₁, M₁₁).
    fn warp_decoder(m: [[f64; 2]; 2]) -> Array2<f64> {
        let mut b = Array2::<f64>::zeros((3, 2));
        b[[1, 0]] = m[0][0];
        b[[1, 1]] = m[1][0];
        b[[2, 0]] = m[0][1];
        b[[2, 1]] = m[1][1];
        b
    }

    /// A grid of fitted coordinates over the patch box `[0, 1]²`.
    fn patch_coords() -> Array2<f64> {
        let g = 9usize;
        let mut c = Array2::<f64>::zeros((g * g, 2));
        for i in 0..g {
            for j in 0..g {
                let row = i * g + j;
                c[[row, 0]] = i as f64 / (g - 1) as f64;
                c[[row, 1]] = j as f64 / (g - 1) as f64;
            }
        }
        c
    }

    /// Scale-invariant anisotropy defect of a constant pullback metric `G`
    /// against the flat reference — the optimum the minimizer chases is 0.
    fn flat_defect_of_constant_metric(g: [f64; 3], n: usize) -> f64 {
        // ḡ = sqrt(det G); Ĝ = G/ḡ; profile c against I; sum over n identical
        // rows. With one constant metric the per-row defect is identical.
        let det = g[0] * g[1] - g[2] * g[2];
        let g_bar = det.sqrt();
        let h = [g[0] / g_bar, g[1] / g_bar, g[2] / g_bar];
        // c = <Ĝ, I> / <I, I> = (h00 + h11) / 2 ; residual = Ĝ − c·I.
        let c = 0.5 * (h[0] + h[1]);
        let r00 = h[0] - c;
        let r11 = h[1] - c;
        let r01 = h[2];
        n as f64 * (r00 * r00 + r11 * r11 + 2.0 * r01 * r01)
    }

    #[test]
    fn planted_warped_patch_recovers_uniform_speed_coords() {
        // A deliberately anisotropic + sheared affine warp M: the fitted chart
        // stretches axis 0 by 1.6, axis 1 by 0.8, with a 0.5 shear.
        let m = [[1.6, 0.5], [0.0, 0.8]];
        let ev = MockPatchEvaluator;
        let decoder = warp_decoder(m);
        let coords = patch_coords();
        let n = coords.nrows();

        // Pullback metric G = MᵀM (constant over the patch).
        let g00 = m[0][0] * m[0][0] + m[1][0] * m[1][0];
        let g11 = m[0][1] * m[0][1] + m[1][1] * m[1][1];
        let g01 = m[0][0] * m[0][1] + m[1][0] * m[1][1];
        let defect_initial = flat_defect_of_constant_metric([g00, g11, g01], n);
        assert!(
            defect_initial > 1e-2,
            "the planted anisotropic warp must start with a sizeable defect; got {defect_initial:.3e}"
        );

        let repar = patch_isometry_flow_reparameterization(&ev, decoder.view(), coords.view())
            .expect("patch reparameterization must evaluate")
            .expect("a warped patch with a global flow basis must canonicalize");

        // Acceptance (#1019): the canonical chart recovers uniform-speed coords
        // — the residual defect is within 10% of the optimum (0). The optimum is
        // exactly 0 (a flat patch IS isometric to itself), so the bar is a small
        // absolute fraction of the initial defect.
        assert!(
            repar.defect_final <= 0.10 * defect_initial,
            "canonicalization must drive the anisotropy defect to within 10% of the optimum; \
             initial {defect_initial:.3e}, final {:.3e}",
            repar.defect_final
        );
        // The defect-final the struct reports must match what the optimizer
        // descended to (strict improvement over the identity flow).
        assert!(
            repar.defect_final < repar.defect_initial,
            "the pass must report a strict improvement; initial {:.3e}, final {:.3e}",
            repar.defect_initial,
            repar.defect_final
        );
        // The canonical chart is a genuine diffeomorphism (guard cleared).
        assert!(
            repar.min_flow_jacobian_det > PATCH_FLOW_DIFFEO_MIN_DET,
            "the canonical flow must be fold-free; min det {:.3e}",
            repar.min_flow_jacobian_det
        );
        // Image frozen: the recomposition residual is within the honest gate.
        assert!(
            repar.recomposition_residual <= CHART_RECOMPOSITION_REL_TOL,
            "the decoded image must be reproduced within the recomposition tolerance; got {:.3e}",
            repar.recomposition_residual
        );
    }

    #[test]
    fn already_uniform_patch_is_left_as_fitted() {
        // M = I: the fitted chart is already uniform-speed (defect 0), so the
        // minimizer finds no strict improvement and honestly skips.
        let ev = MockPatchEvaluator;
        let decoder = warp_decoder([[1.0, 0.0], [0.0, 1.0]]);
        let coords = patch_coords();
        let out = patch_isometry_flow_reparameterization(&ev, decoder.view(), coords.view())
            .expect("patch reparameterization must evaluate");
        assert!(
            out.is_none(),
            "an already-uniform patch chart must be left as fitted (honest skip), got Some"
        );
    }

    #[test]
    fn collapsed_patch_axis_is_refused() {
        // All rows share a single t₁ value: the patch has no extent on axis 1,
        // so there is no honest flat chart — the basis builder refuses and the
        // function returns None rather than erroring.
        let ev = MockPatchEvaluator;
        let decoder = warp_decoder([[1.3, 0.0], [0.0, 0.9]]);
        let mut coords = patch_coords();
        for row in 0..coords.nrows() {
            coords[[row, 1]] = 0.5;
        }
        let out = patch_isometry_flow_reparameterization(&ev, decoder.view(), coords.view())
            .expect("patch reparameterization must evaluate");
        assert!(
            out.is_none(),
            "a patch collapsed along one axis must be refused (None), got Some"
        );
    }

    /// Finite-difference check of the analytic mode gradients (the `grad` the
    /// Gauss–Newton Jacobian is built from) against the monomial values.
    #[test]
    fn free_patch_mode_gradients_match_finite_difference() {
        let basis = FreePatchFlowBasis::new([-0.5, 0.2], [1.5, 2.2]).expect("patch basis");
        let t = [0.3, 1.1];
        let eps = 1e-6;
        let base = basis.mode_samples(t);
        let plus0 = basis.mode_samples([t[0] + eps, t[1]]);
        let plus1 = basis.mode_samples([t[0], t[1] + eps]);
        for k in 0..base.len() {
            let fd0 = (plus0[k].value - base[k].value) / eps;
            let fd1 = (plus1[k].value - base[k].value) / eps;
            assert!(
                (fd0 - base[k].grad[0]).abs() < 1e-4,
                "mode {k} ∂/∂t₀ FD {fd0} vs analytic {}",
                base[k].grad[0]
            );
            assert!(
                (fd1 - base[k].grad[1]).abs() < 1e-4,
                "mode {k} ∂/∂t₁ FD {fd1} vs analytic {}",
                base[k].grad[1]
            );
        }
    }
}

#[cfg(test)]
mod sphere_defect_tests {
    use super::*;
    use ndarray::Array2;

    /// FD-gate the analytic conformal-boost mode Jacobians `Dv_k(t)` against
    /// central differences of the displacements `v_k(t)` — the exact object the
    /// sphere Gauss–Newton residual derivative is built from.
    #[test]
    fn sphere_boost_mode_jacobians_match_finite_difference() {
        let eps = 1e-6;
        // Two interior points well off the poles (cos lat bounded away from 0).
        for t in [[0.3, 0.7], [-0.6, -1.2]] {
            let jac = SphereBoostFlowBasis::mode_jacobians(t);
            let vp0 = SphereBoostFlowBasis::mode_displacements([t[0] + eps, t[1]]);
            let vm0 = SphereBoostFlowBasis::mode_displacements([t[0] - eps, t[1]]);
            let vp1 = SphereBoostFlowBasis::mode_displacements([t[0], t[1] + eps]);
            let vm1 = SphereBoostFlowBasis::mode_displacements([t[0], t[1] - eps]);
            for k in 0..3 {
                for comp in 0..2 {
                    let fd_dlat = (vp0[k][comp] - vm0[k][comp]) / (2.0 * eps);
                    let fd_dlon = (vp1[k][comp] - vm1[k][comp]) / (2.0 * eps);
                    assert!(
                        (fd_dlat - jac[k][comp][0]).abs() < 1e-5,
                        "boost {k} comp {comp} ∂/∂lat FD {fd_dlat} vs analytic {} at {t:?}",
                        jac[k][comp][0]
                    );
                    assert!(
                        (fd_dlon - jac[k][comp][1]).abs() < 1e-5,
                        "boost {k} comp {comp} ∂/∂lon FD {fd_dlon} vs analytic {} at {t:?}",
                        jac[k][comp][1]
                    );
                }
            }
        }
    }

    /// Test-only finite-difference oracle for the complete production Jacobian:
    /// boost flow, moved-latitude whitening, `AᵀA`, and the profiled-scale
    /// projection. The perturbed residual uses each perturbation's own profiled
    /// scale, so this catches the former fixed-scale Jacobian as well as a missing
    /// whitening chain-rule term.
    #[test]
    fn sphere_profiled_residual_jacobian_matches_finite_difference_oracle() {
        let row_coords = Array2::from_shape_vec(
            (5, 2),
            vec![-0.72, -1.1, -0.31, 0.8, 0.04, -0.2, 0.43, 1.7, 0.77, -2.0],
        )
        .expect("five sphere rows");
        // Positive-definite, non-diagonal fitted metrics exercise all three
        // symmetric entries and make the scale projection observable.
        let ghat = vec![
            [1.20, 0.82, 0.11],
            [0.91, 1.31, -0.14],
            [1.43, 0.73, 0.19],
            [0.79, 1.18, 0.08],
            [1.08, 0.96, -0.17],
        ];
        let ghat_norm_sq = ghat
            .iter()
            .map(|g| g[0] * g[0] + g[1] * g[1] + 2.0 * g[2] * g[2])
            .sum::<f64>();
        let theta = vec![0.07, -0.05, 0.04];
        let state = sphere_eval_boost_defect(&theta, row_coords.view(), &ghat, ghat_norm_sq)
            .expect("interior boost state");
        let (residual, jacobian) =
            sphere_boost_residual_jacobian(&theta, row_coords.view(), &ghat, ghat_norm_sq, &state)
                .expect("analytic profiled residual Jacobian");

        let profiled_residual = |trial: &[f64]| {
            let trial_state =
                sphere_eval_boost_defect(trial, row_coords.view(), &ghat, ghat_norm_sq)
                    .expect("test perturbation remains in the sphere chart");
            let mut out = Array2::<f64>::zeros((3 * ghat.len(), 1));
            for (row, (a, g)) in trial_state.a_rows.iter().zip(ghat.iter()).enumerate() {
                let m00 = a[0] * a[0] + a[2] * a[2];
                let m11 = a[1] * a[1] + a[3] * a[3];
                let m01 = a[0] * a[1] + a[2] * a[3];
                out[[3 * row, 0]] = m00 - trial_state.scale * g[0];
                out[[3 * row + 1, 0]] = m11 - trial_state.scale * g[1];
                out[[3 * row + 2, 0]] = std::f64::consts::SQRT_2 * (m01 - trial_state.scale * g[2]);
            }
            out
        };
        let direct_residual = profiled_residual(&theta);
        for row in 0..residual.nrows() {
            assert_eq!(residual[[row, 0]], direct_residual[[row, 0]]);
        }

        let eps = 1.0e-6;
        for k in 0..theta.len() {
            let mut plus = theta.clone();
            let mut minus = theta.clone();
            plus[k] += eps;
            minus[k] -= eps;
            let r_plus = profiled_residual(&plus);
            let r_minus = profiled_residual(&minus);
            for row in 0..jacobian.nrows() {
                let fd = (r_plus[[row, 0]] - r_minus[[row, 0]]) / (2.0 * eps);
                let exact = jacobian[[row, k]];
                let tolerance = 2.0e-7 * (1.0 + exact.abs());
                assert!(
                    (fd - exact).abs() <= tolerance,
                    "profiled sphere residual row {row}, boost {k}: FD {fd:.12e}, analytic {exact:.12e}"
                );
            }

            // Variable projection makes every residual-Jacobian column
            // orthogonal to the stacked profiled metric direction.
            let mut projection = 0.0_f64;
            for (row, g) in ghat.iter().enumerate() {
                projection += g[0] * jacobian[[3 * row, k]]
                    + g[1] * jacobian[[3 * row + 1, k]]
                    + std::f64::consts::SQRT_2 * g[2] * jacobian[[3 * row + 2, k]];
            }
            assert!(
                projection.abs() <= 1.0e-10,
                "profile projection failed for boost {k}: gᵀJ={projection:.12e}"
            );
        }
    }

}

#[cfg(test)]
mod turning_tests {
    use super::*;
    use crate::basis::PeriodicHarmonicEvaluator;
    use ndarray::{Array1, Array2};
    use std::f64::consts::TAU;

    /// A unit circle `γ(t) = (cos 2πt, sin 2πt)` traversed once over `t ∈ [0,1]`
    /// has constant curvature κ = 1 and speed 2π, so the total turning is
    /// `Θ = ∫₀¹ κ·‖γ'‖ dt = 2π`.
    #[test]
    fn full_circle_turning_is_two_pi() {
        let ev = PeriodicHarmonicEvaluator::new(3).expect("3-basis circle");
        // Basis [1, sin2πt, cos2πt]; B maps cos→x (output 0), sin→y (output 1).
        let mut decoder = Array2::<f64>::zeros((3, 2));
        decoder[[2, 0]] = 1.0; // cos -> x
        decoder[[1, 1]] = 1.0; // sin -> y
        // Span the full period [0, 1]; the integral runs over [min, max] of the
        // supplied coordinates, so the endpoint t = 1 must be present to close
        // the loop (the circle's speed 2π is nonzero everywhere, no stationary
        // point at the seam).
        let coords = Array1::from_iter((0..=50).map(|i| i as f64 / 50.0));
        let theta = d1_atom_fitted_turning(&ev, decoder.view(), coords.view())
            .expect("turning must evaluate")
            .expect("a non-degenerate circle must return Some");
        assert!(
            (theta - TAU).abs() < 1e-6,
            "a full unit circle has total turning 2π; got {theta:.9}"
        );
    }

    /// A half circle (`t ∈ [0, 0.5]`) turns through π.
    #[test]
    fn half_circle_turning_is_pi() {
        let ev = PeriodicHarmonicEvaluator::new(3).expect("3-basis circle");
        let mut decoder = Array2::<f64>::zeros((3, 2));
        decoder[[2, 0]] = 1.0;
        decoder[[1, 1]] = 1.0;
        let coords = Array1::from_iter((0..=25).map(|i| 0.5 * i as f64 / 25.0));
        let theta = d1_atom_fitted_turning(&ev, decoder.view(), coords.view())
            .expect("turning must evaluate")
            .expect("a non-degenerate half-circle must return Some");
        assert!(
            (theta - std::f64::consts::PI).abs() < 1e-6,
            "a half circle turns through π; got {theta:.9}"
        );
    }

    /// A straight line (only the constant + one linear-in-image direction, no
    /// genuine curvature) has zero turning — exactly the `Θ → 0` linear-tail
    /// signature the EV-vs-Θ measurement uses to flag a curved atom that is
    /// really just a linear direction. Here the decoder uses a SINGLE harmonic
    /// pair scaled so the image is a 1-D segment (x and y both ∝ cos 2πt), i.e.
    /// the decoded curve lies on a line through the origin → wedge ≡ 0.
    #[test]
    fn straight_line_image_has_zero_turning() {
        let ev = PeriodicHarmonicEvaluator::new(3).expect("3-basis circle");
        let mut decoder = Array2::<f64>::zeros((3, 2));
        // Both outputs ∝ cos 2πt → the image is the line y = 2x, γ collinear with
        // γ', γ'' at every t, so the wedge norm vanishes identically.
        decoder[[2, 0]] = 1.0;
        decoder[[2, 1]] = 2.0;
        // Span a coordinate range strictly inside (0, 0.25) so the speed
        // `‖γ'‖ ∝ |sin 2πt|` never hits the stationary zero at t = 0 (where the
        // turning integrand is genuinely undefined and the function refuses).
        let coords = Array1::from_iter((0..=20).map(|i| 0.05 + 0.15 * i as f64 / 20.0));
        let theta = d1_atom_fitted_turning(&ev, decoder.view(), coords.view())
            .expect("turning must evaluate")
            .expect("a non-degenerate segment must return Some");
        assert!(
            theta < 1e-9,
            "a straight-line image has zero turning (the linear-tail signature); got {theta:.3e}"
        );
    }

    /// #1610/#1026 — a GLOBALLY constant image (only the DC basis row is
    /// nonzero) decodes to a single fixed point: `γ'(t) ≡ 0` at every node, so
    /// the curve is stationary everywhere. A point has no arc to turn through,
    /// so its total turning is exactly `0` — the ultimate linear-tail signature.
    /// This is the path the hybrid-collapse witness relies on (a `d = 1` atom
    /// straightened to its DC component must read `Some(0.0)`, not the historical
    /// `None`, so its slot can collapse to a legitimate linear tail).
    #[test]
    fn constant_image_turning_is_some_zero() {
        let ev = PeriodicHarmonicEvaluator::new(3).expect("3-basis circle");
        let mut decoder = Array2::<f64>::zeros((3, 2));
        // Only the DC (constant) basis row → γ(t) is a fixed point for all t.
        decoder[[0, 0]] = 0.7;
        decoder[[0, 1]] = -0.3;
        let coords = Array1::from_iter((0..=20).map(|i| 0.05 + 0.15 * i as f64 / 20.0));
        let theta = d1_atom_fitted_turning(&ev, decoder.view(), coords.view())
            .expect("turning must evaluate")
            .expect("a globally constant image returns Some(0.0), not None");
        assert_eq!(
            theta, 0.0,
            "a globally constant (all-stationary) image has exactly zero turning"
        );
    }

    /// #1610/#1026 — a curve that is stationary at SOME nodes but moving at
    /// others is a cusp: the unsigned curvature integrand is genuinely
    /// ill-defined at the stationary node, so the function REFUSES (`None`)
    /// rather than under-count the sharp turn. Built as an image `∝ cos 2πt` (a
    /// line through the origin) over a span whose lower endpoint is exactly
    /// `t = 0`, where `γ' ∝ -sin 2πt` vanishes: that endpoint lands on a grid
    /// node and is stationary while every interior node moves — the mixed case
    /// the constant-vs-cusp split must keep conservative.
    #[test]
    fn partial_cusp_turning_is_none() {
        let ev = PeriodicHarmonicEvaluator::new(3).expect("3-basis circle");
        let mut decoder = Array2::<f64>::zeros((3, 2));
        decoder[[2, 0]] = 1.0;
        decoder[[2, 1]] = 2.0;
        // Span [0, 0.2]: the Simpson grid's lower node sits exactly on t = 0,
        // where the speed `‖γ'‖ ∝ |sin 2πt|` is zero (stationary), while every
        // interior node is moving.
        let coords = Array1::from_iter((0..=20).map(|i| 0.2 * i as f64 / 20.0));
        let result = d1_atom_fitted_turning(&ev, decoder.view(), coords.view())
            .expect("turning must evaluate");
        assert!(
            result.is_none(),
            "a curve stationary at some nodes but moving at others (a cusp) must \
             refuse with None; got {result:?}"
        );
    }
}

#[cfg(test)]
mod certified_diffeo_guard_2518_tests {
    use super::*;

    /// #2518 item 3 — the guard's contract is `det Dφ_θ > δ` EVERYWHERE, so its
    /// decision quantity must be a lower bound on the continuum minimum, not a
    /// sample. Resample `det Dφ_θ` on a much finer grid deliberately OFFSET from
    /// the guard nodes — so it sees exactly the inter-node values the guard
    /// cannot — and require the certificate to sit below everything it finds.
    #[test]
    fn certified_torus_jacobian_det_lower_bounds_the_offgrid_minimum_2518() {
        let period = 1.0_f64;
        let basis = TorusFlowBasis::new(period).expect("torus flow basis");
        let dim = basis.dim();

        // Identity flow: `det Dφ ≡ 1`, the grid is already exact, and the
        // correction must be EXACTLY zero — a certificate that charged a
        // near-identity flow would be unusable.
        let zero = vec![0.0_f64; dim];
        assert_eq!(
            basis.jacobian_det_grid_correction(&zero),
            0.0,
            "the identity flow has a constant determinant; its correction is exactly 0"
        );
        assert!(
            (basis.certified_min_jacobian_det(&zero) - 1.0).abs() < 1.0e-12,
            "identity flow must certify det = 1"
        );

        // Off-grid dense resample helper: 4x the guard resolution, shifted by
        // half the dense spacing so no dense node coincides with a guard node.
        let offgrid_min = |theta: &[f64]| -> f64 {
            let dense = 256usize;
            let shift = 0.5 * period / dense as f64;
            let mut worst = f64::INFINITY;
            for i in 0..dense {
                for j in 0..dense {
                    let t = [
                        period * i as f64 / dense as f64 + shift,
                        period * j as f64 / dense as f64 + shift,
                    ];
                    let jac = basis.flow_jacobian(theta, t);
                    worst = worst.min(jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0]);
                }
            }
            worst
        };
        let mixed_signs = |scale: f64| -> Vec<f64> {
            (0..dim)
                .map(|k| {
                    let sign = if k % 2 == 0 { 1.0 } else { -1.0 };
                    scale * ((k % 7) as f64 - 3.0) * sign
                })
                .collect::<Vec<f64>>()
        };

        // SOUNDNESS, at two amplitudes an order of magnitude apart. The bound is
        // the guard's whole contract, so it must hold everywhere the guard could
        // be asked, not at one tuned point.
        for scale in [2.0e-5_f64, 4.0e-4] {
            let theta = mixed_signs(scale);
            let grid_min = basis.min_jacobian_det_on_grid(&theta);
            let correction = basis.jacobian_det_grid_correction(&theta);
            let certified = basis.certified_min_jacobian_det(&theta);
            let dense_min = offgrid_min(&theta);
            println!(
                "[#2518 torus det certificate] scale={scale:.1e} grid_min={grid_min:.9} \
                 correction={correction:.9} certified={certified:.9} offgrid_min={dense_min:.9}"
            );
            assert!(
                correction > 0.0,
                "a flow with a varying Jacobian must pay a positive inter-node correction"
            );
            assert!(
                certified < grid_min,
                "the certificate must sit strictly below the sample it corrects: \
                 {certified} vs {grid_min}"
            );
            assert!(
                certified <= dense_min,
                "certified bound {certified} must not exceed the off-grid minimum \
                 {dense_min} (correction {correction}, grid minimum {grid_min})"
            );
        }

        // NON-VACUITY. The correction is `O(|θ|)` precisely so a near-identity
        // flow still clears the fold floor — a certificate that refused every
        // flow would satisfy the soundness arm above and be worthless.
        let near_identity = mixed_signs(2.0e-5);
        let certified_near = basis.certified_min_jacobian_det(&near_identity);
        assert!(
            certified_near > TORUS_FLOW_DIFFEO_MIN_DET,
            "a near-identity flow must still certify above the fold floor δ={}; got \
             {certified_near}",
            TORUS_FLOW_DIFFEO_MIN_DET
        );
    }

    /// The correction is a bound on how far `det` can move between nodes, so it
    /// must vanish with the flow amplitude and grow with it — otherwise it is a
    /// constant tax rather than a derived quantity.
    #[test]
    fn certified_torus_correction_tracks_the_flow_amplitude_2518() {
        let basis = TorusFlowBasis::new(1.0).expect("torus flow basis");
        let dim = basis.dim();
        let base: Vec<f64> = (0..dim)
            .map(|k| 0.001 * ((k % 5) as f64 - 2.0))
            .collect();
        let small: Vec<f64> = base.iter().map(|v| v * 0.5).collect();
        let large: Vec<f64> = base.iter().map(|v| v * 2.0).collect();
        let c_small = basis.jacobian_det_grid_correction(&small);
        let c_base = basis.jacobian_det_grid_correction(&base);
        let c_large = basis.jacobian_det_grid_correction(&large);
        assert!(
            c_small < c_base && c_base < c_large,
            "correction must be strictly increasing in the flow amplitude: \
             {c_small} < {c_base} < {c_large}"
        );
        assert!(
            basis.jacobian_det_grid_correction(&vec![0.0_f64; dim]) == 0.0,
            "zero amplitude pays nothing"
        );
    }

    /// #2518 item 0 — the free-patch guard's injectivity argument is DEGREE-1
    /// ONLY: `det Dφ > 0` ⇒ injective holds for an affine map by linear algebra
    /// and is FALSE for higher-degree polynomial maps (Pinchuk 1994 constructs a
    /// real polynomial counterexample). `PATCH_FLOW_GUARD_NODES_PER_AXIS` is
    /// documented against that proof, not against grid resolution, so pin the
    /// degree here: raising it must break a test rather than silently invalidate
    /// the guard's justification.
    #[test]
    fn patch_flow_degree_is_pinned_to_its_injectivity_proof_2518() {
        assert_eq!(
            PATCH_FLOW_MAX_DEGREE, 1,
            "PATCH_FLOW_GUARD_NODES_PER_AXIS = {} is justified by the flow being AFFINE: \
             Dφ_θ is then CONSTANT (so any grid is exact) and det ≠ 0 gives global \
             injectivity by linear algebra. At degree > 1 BOTH facts fail and no node count \
             restores either — the guard needs a new injectivity argument first, and the \
             inter-node correction of TorusFlowBasis::certified_min_jacobian_det for the \
             sampling half",
            PATCH_FLOW_GUARD_NODES_PER_AXIS
        );
    }
}