gam-sae 0.3.151

use ndarray::{Array2, Array3, Array4, Array5, ArrayView2};
use std::sync::Arc;

pub trait SaeBasisEvaluator: Send + Sync + std::fmt::Debug {
    fn evaluate(&self, coords: ArrayView2<'_, f64>) -> Result<(Array2<f64>, Array3<f64>), String>;

    /// Evaluate `Φ` and its first-location jet DIRECTLY into caller-owned
    /// buffers, avoiding the fresh `(N, M)` Φ + `(N, M, d)` jet allocation
    /// [`Self::evaluate`] performs on every call.
    ///
    /// The inner Newton loop re-evaluates each atom's basis on every
    /// line-search trial ([`crate::manifold::atom::SaeManifoldAtom::refresh_basis`]);
    /// there the term already owns correctly-shaped Φ / jet arrays, so reusing
    /// them removes the per-trial allocation churn (multiple trials per Newton
    /// iteration under LM / backtracking). `phi` must be shaped
    /// `(coords.nrows(), M)` and `jet` `(coords.nrows(), M, coords.ncols())`
    /// for this evaluator's basis width `M`; implementations return a
    /// descriptive `Err` on any mismatch — the same shape guard `refresh_basis`
    /// previously applied to the freshly-allocated arrays.
    ///
    /// The default forwards to [`Self::evaluate`] and copies, so evaluators
    /// that have not specialized it stay correct at the cost of one extra copy;
    /// the hot evaluators override it to fill in place.
    fn evaluate_into(
        &self,
        phi: &mut Array2<f64>,
        jet: &mut Array3<f64>,
        coords: ArrayView2<'_, f64>,
    ) -> Result<(), String> {
        let (new_phi, new_jet) = self.evaluate(coords)?;
        if new_phi.dim() != phi.dim() {
            return Err(format!(
                "SaeBasisEvaluator::evaluate_into: evaluator returned Φ {:?}, target buffer {:?}",
                new_phi.dim(),
                phi.dim()
            ));
        }
        if new_jet.dim() != jet.dim() {
            return Err(format!(
                "SaeBasisEvaluator::evaluate_into: evaluator returned jet {:?}, target buffer {:?}",
                new_jet.dim(),
                jet.dim()
            ));
        }
        phi.assign(&new_phi);
        jet.assign(&new_jet);
        Ok(())
    }

    /// Return the same evaluator after the coordinate change
    /// `old_t = shift + scale * new_t`, when the basis family can transport the
    /// decoder coefficients exactly enough for the accepted-iterate gauge fix.
    fn affine_transformed_evaluator(
        &self,
        shift: &[f64],
        scale: &[f64],
        n_basis: usize,
    ) -> Result<Option<Arc<dyn SaeBasisSecondJet>>, String> {
        if shift.len() == usize::MAX || scale.len() == usize::MAX || n_basis == usize::MAX {
            return Err("SaeBasisEvaluator::affine_transformed_evaluator: unreachable affine metadata width".to_string());
        }
        Ok(None)
    }

    /// Column split for the curvature homotopy `Phi_eta = [base, eta*curved]`.
    ///
    /// `base` columns are held FIXED as `eta` walks `0 → 1`; `curved` columns are
    /// scaled by `eta`. The `eta = 0` endpoint is therefore the *base-topology*
    /// relaxation — the atom evaluated on its base (η-invariant) columns only —
    /// NOT a linear/affine model: for the harmonic and sphere-chart bases the
    /// base block already carries extrinsic curvature (a first-harmonic
    /// `[sin, cos]` pair traces a circle, the sphere chart's `[x, y, z]` block
    /// traces the unit sphere). "Base" names "does not scale with η", never
    /// "curvature-free". See [`PhiEtaSplit`].
    ///
    /// The default is a flat (genuinely affine) monomial-style basis where every
    /// column is a base column. Curved atom evaluators override this with their
    /// topology-specific split; callers pass `n_basis` so the split is checked
    /// against the concrete design width currently being evaluated.
    fn phi_eta_split(&self, n_basis: usize) -> Result<PhiEtaSplit, String> {
        Ok(PhiEtaSplit::all_base(n_basis))
    }

    /// Per-factor basis sizes `(M₁, M₂)` of a `d = 2` TENSOR-PRODUCT evaluator
    /// whose fused basis is the Kronecker product of two factor bases in
    /// row-major column order `flat = j·M₂ + k` (the within-atom carve / #993
    /// convention). `M₁·M₂` equals the fused basis width. `None` for evaluators
    /// that are not a two-factor product (single-axis, sphere chart, monomial
    /// patch) — the within-atom functional-ANOVA carve is only defined on a
    /// genuine product manifold, so a `None` here is the honest "no factor
    /// split" signal, never a guess.
    fn factor_basis_sizes(&self) -> Option<(usize, usize)> {
        None
    }

    /// Evaluate the basis at curvature scale `eta in [0, 1]` plus the analytic
    /// derivative with respect to eta.
    ///
    /// At `eta == 1.0` this leaves the existing basis and jet arrays untouched,
    /// so the returned `phi`/`jet` are exactly the same values as [`Self::evaluate`].
    fn evaluate_phi_eta(
        &self,
        coords: ArrayView2<'_, f64>,
        eta: f64,
    ) -> Result<PhiEtaEvaluation, String> {
        if !(eta.is_finite() && (0.0..=1.0).contains(&eta)) {
            return Err(format!(
                "SaeBasisEvaluator::evaluate_phi_eta: eta must be finite in [0, 1]; got {eta}"
            ));
        }
        let (mut phi, mut jet) = self.evaluate(coords)?;
        let split = self.phi_eta_split(phi.ncols())?;
        let mut dphi_deta = Array2::<f64>::zeros(phi.dim());
        let mut djet_deta = Array3::<f64>::zeros(jet.dim());
        for &col in &split.curved_cols {
            if col >= phi.ncols() {
                return Err(format!(
                    "SaeBasisEvaluator::evaluate_phi_eta: curved column {col} exceeds basis width {}",
                    phi.ncols()
                ));
            }
            for row in 0..phi.nrows() {
                dphi_deta[[row, col]] = phi[[row, col]];
                if eta != 1.0 {
                    phi[[row, col]] *= eta;
                }
                for axis in 0..jet.shape()[2] {
                    djet_deta[[row, col, axis]] = jet[[row, col, axis]];
                    if eta != 1.0 {
                        jet[[row, col, axis]] *= eta;
                    }
                }
            }
        }
        Ok(PhiEtaEvaluation {
            phi,
            jet,
            dphi_deta,
            djet_deta,
            split,
        })
    }

    /// Object-safe forwarder to [`SaeBasisSecondJet::second_jet`] for callers
    /// holding `&dyn SaeBasisEvaluator` / `Arc<dyn SaeBasisEvaluator>`.
    ///
    /// Implementations return `Some(result)` only when an analytic second jet
    /// exists for this evaluator. Returning `None` is an explicit capability
    /// declaration, not a default sentinel hidden in the trait.
    fn second_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array4<f64>, String>>;

    /// Object-safe forwarder to the basis third jet
    /// `T[n, m, a, c, e] = ∂³Φ_m / ∂t_a ∂t_c ∂t_e`, for callers holding
    /// `&dyn SaeBasisEvaluator` / `Arc<dyn SaeBasisSecondJet>`. The exact
    /// isometry Hessian (`IsometryPenalty::hvp`) needs the *decoder* third jet
    /// `K = Σ_m T[..,m,..]·B[m,:]` for its residual·curvature term; without it
    /// that exact Hessian silently drops the residual and collapses to
    /// Gauss-Newton (issue #458).
    ///
    /// Implementations return `Some(result)` only when an analytic third jet
    /// exists for this evaluator. Evaluators without one return `None`
    /// explicitly; there is no finite-difference fallback.
    fn third_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array5<f64>, String>>;
}

/// Curvature-homotopy column split `Phi_eta = [base, eta*curved]`.
///
/// `curved_cols` are the columns scaled by the dial `eta`; `base_cols` are the
/// columns held fixed across the whole walk. "Base" means "η-invariant", NOT
/// "linear/curvature-free": for the harmonic and sphere-chart bases the base
/// block already embeds curvature (first-harmonic `[sin, cos]`, the sphere
/// chart's `[x, y, z]`). The `eta = 0` endpoint is thus the base-topology
/// relaxation, not an affine/Eckart-Young linear model. The genuine
/// low-rank (Eckart-Young / PCA) certificate lives in
/// [`crate::manifold::outer_objective::linear_span_anchor`] and is a rank
/// ceiling that bounds every `eta`, independent of this split.
#[derive(Debug, Clone, PartialEq, Eq)]
pub struct PhiEtaSplit {
    pub base_cols: Vec<usize>,
    pub curved_cols: Vec<usize>,
}

impl PhiEtaSplit {
    pub fn all_base(n_basis: usize) -> Self {
        Self {
            base_cols: (0..n_basis).collect(),
            curved_cols: Vec::new(),
        }
    }

    fn from_curved_mask(mask: Vec<bool>) -> Self {
        let mut base_cols = Vec::new();
        let mut curved_cols = Vec::new();
        for (col, curved) in mask.into_iter().enumerate() {
            if curved {
                curved_cols.push(col);
            } else {
                base_cols.push(col);
            }
        }
        Self {
            base_cols,
            curved_cols,
        }
    }
}

#[derive(Debug, Clone)]
pub struct PhiEtaEvaluation {
    pub phi: Array2<f64>,
    pub jet: Array3<f64>,
    pub dphi_deta: Array2<f64>,
    pub djet_deta: Array3<f64>,
    pub split: PhiEtaSplit,
}

fn monomial_linear_mask(dimension: usize, max_total_degree: usize) -> Vec<bool> {
    gam_terms::basis::monomial_exponents(dimension, max_total_degree)
        .iter()
        .map(|alpha| alpha.iter().sum::<usize>() <= 1)
        .collect()
}

fn duchon_effective_order_for_eta(
    centers: ArrayView2<'_, f64>,
    order: gam_terms::basis::DuchonNullspaceOrder,
) -> gam_terms::basis::DuchonNullspaceOrder {
    let mut effective = order;
    while effective != gam_terms::basis::DuchonNullspaceOrder::Zero
        && centers.nrows() <= duchon_polynomial_column_count(centers.ncols(), effective)
    {
        effective = match effective {
            gam_terms::basis::DuchonNullspaceOrder::Zero => {
                gam_terms::basis::DuchonNullspaceOrder::Zero
            }
            gam_terms::basis::DuchonNullspaceOrder::Linear => {
                gam_terms::basis::DuchonNullspaceOrder::Zero
            }
            gam_terms::basis::DuchonNullspaceOrder::Degree(2) => {
                gam_terms::basis::DuchonNullspaceOrder::Linear
            }
            gam_terms::basis::DuchonNullspaceOrder::Degree(k) => {
                gam_terms::basis::DuchonNullspaceOrder::Degree(k - 1)
            }
        };
    }
    effective
}

fn duchon_polynomial_column_count(
    dimension: usize,
    order: gam_terms::basis::DuchonNullspaceOrder,
) -> usize {
    match order {
        gam_terms::basis::DuchonNullspaceOrder::Zero => 1,
        gam_terms::basis::DuchonNullspaceOrder::Linear => dimension + 1,
        gam_terms::basis::DuchonNullspaceOrder::Degree(degree) => {
            gam_terms::basis::monomial_exponents(dimension, degree).len()
        }
    }
}

/// Bases that expose an analytic second jet
/// `H[n, m, a, c] = ∂²Phi_k[n, m] / (∂t_{n,a} ∂t_{n,c})`,
/// shape `(n_rows, n_basis, latent_dim, latent_dim)`.
///
/// Implemented only by evaluators with a closed-form Hessian (periodic
/// harmonic, sphere chart, torus). Callers that need an analytic
/// `∂J/∂t` require this bound; evaluators without it must use a
/// derivative-free fallback. Replaces the previous `Option<Array4<f64>>`
/// return on the base trait so the "no second jet" case is encoded by
/// trait absence rather than a sentinel `None`, and shape mismatches
/// surface as descriptive errors instead of silently collapsing to
/// `None`.
pub trait SaeBasisSecondJet: SaeBasisEvaluator {
    fn second_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array4<f64>, String>;
}

/// Bases that expose an analytic third jet
/// `T[n, m, a, c, e] = ∂³Φ_m[n] / (∂t_{n,a} ∂t_{n,c} ∂t_{n,e})`,
/// shape `(n_rows, n_basis, latent_dim, latent_dim, latent_dim)`.
///
/// The exact isometry Hessian (`IsometryPenalty::hvp`) needs the third decoder
/// jet `K = ∂³φ/∂t³ = Σ_m T[..,m,..] · B[m, :]` for its residual·curvature term
/// `B_{ab,cd} = K_{a,cd}ᵀ W J_b + H_{a,c}ᵀ W H_{b,d} + H_{a,d}ᵀ W H_{b,c}
/// + J_aᵀ W K_{b,cd}`. Bases that supply a closed-form `H` (the
/// [`SaeBasisSecondJet`] super-bound) but not `K` leave that exact Hessian
/// silently dropping the residual term; this trait closes that gap for every
/// analytic basis: the curved bases (sphere chart, periodic harmonic, torus
/// harmonic), the Euclidean monomial patch, the trivially-zero affine basis,
/// and the Duchon basis (radial third-derivative kernel block + monomial
/// nullspace block, both in closed form). The full third jet is symmetric in
/// its three trailing axes.
pub trait SaeBasisThirdJet: SaeBasisSecondJet {
    fn third_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array5<f64>, String>;
}

/// Periodic harmonic basis evaluator for a single-dimensional circle latent.
///
/// Produces `M = 2*num_harmonics + 1` basis functions
/// `[1, sin(2π·1·t), cos(2π·1·t), …, sin(2π·H·t), cos(2π·H·t)]` where
/// `H = (M − 1) / 2`. The latent must have `latent_dim == 1`.
#[derive(Debug, Clone)]
pub struct PeriodicHarmonicEvaluator {
    pub num_basis: usize,
}

impl PeriodicHarmonicEvaluator {
    pub fn new(num_basis: usize) -> Result<Self, String> {
        if num_basis == 0 || num_basis % 2 == 0 {
            return Err(format!(
                "PeriodicHarmonicEvaluator requires odd num_basis >= 1; got {num_basis}"
            ));
        }
        Ok(Self { num_basis })
    }
}

impl SaeBasisEvaluator for PeriodicHarmonicEvaluator {
    fn phi_eta_split(&self, n_basis: usize) -> Result<PhiEtaSplit, String> {
        if n_basis != self.num_basis {
            return Err(format!(
                "PeriodicHarmonicEvaluator::phi_eta_split: n_basis {n_basis} != evaluator width {}",
                self.num_basis
            ));
        }
        // Base (η-invariant) columns: constant `1` plus the FIRST-harmonic
        // `[sin 2πt, cos 2πt]`. These are not linear — the first-harmonic pair
        // already traces the unit circle — but they define the base topology the
        // homotopy relaxes onto; only harmonics ≥ 2 scale with `eta`.
        let mut curved = vec![false; n_basis];
        for h in 2..=(n_basis - 1) / 2 {
            curved[2 * h - 1] = true;
            curved[2 * h] = true;
        }
        Ok(PhiEtaSplit::from_curved_mask(curved))
    }

    fn second_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array4<f64>, String>> {
        Some(<Self as SaeBasisSecondJet>::second_jet(self, coords))
    }

    fn third_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array5<f64>, String>> {
        Some(<Self as SaeBasisThirdJet>::third_jet(self, coords))
    }

    fn evaluate(&self, coords: ArrayView2<'_, f64>) -> Result<(Array2<f64>, Array3<f64>), String> {
        // Single source of truth: allocate the correctly-shaped buffers and let
        // the in-place path fill them, so `evaluate` and `evaluate_into` can
        // never numerically diverge.
        let n = coords.nrows();
        let m = self.num_basis;
        let mut phi = Array2::<f64>::zeros((n, m));
        let mut jet = Array3::<f64>::zeros((n, m, 1));
        self.evaluate_into(&mut phi, &mut jet, coords)?;
        Ok((phi, jet))
    }

    fn evaluate_into(
        &self,
        phi: &mut Array2<f64>,
        jet: &mut Array3<f64>,
        coords: ArrayView2<'_, f64>,
    ) -> Result<(), String> {
        let n = coords.nrows();
        let d = coords.ncols();
        if d != 1 {
            return Err(format!(
                "PeriodicHarmonicEvaluator: expected latent_dim == 1, got {d}"
            ));
        }
        let m = self.num_basis;
        if phi.dim() != (n, m) {
            return Err(format!(
                "PeriodicHarmonicEvaluator::evaluate_into: Φ buffer {:?} != ({n}, {m})",
                phi.dim()
            ));
        }
        if jet.dim() != (n, m, 1) {
            return Err(format!(
                "PeriodicHarmonicEvaluator::evaluate_into: jet buffer {:?} != ({n}, {m}, 1)",
                jet.dim()
            ));
        }
        let num_harmonics = (m - 1) / 2;
        let two_pi = 2.0 * std::f64::consts::PI;
        // The constant column carries a zero jet and is never written in the
        // harmonic loop below, so clear both buffers to erase any stale
        // (reused-workspace) contents before filling.
        phi.fill(0.0);
        jet.fill(0.0);
        for row in 0..n {
            let t = coords[[row, 0]];
            phi[[row, 0]] = 1.0;
            for h in 1..=num_harmonics {
                let angle = two_pi * (h as f64) * t;
                let s = angle.sin();
                let c = angle.cos();
                let s_idx = 2 * h - 1;
                let c_idx = 2 * h;
                phi[[row, s_idx]] = s;
                phi[[row, c_idx]] = c;
                jet[[row, s_idx, 0]] = two_pi * (h as f64) * c;
                jet[[row, c_idx, 0]] = -two_pi * (h as f64) * s;
            }
        }
        Ok(())
    }
}

impl SaeBasisSecondJet for PeriodicHarmonicEvaluator {
    /// Second derivative of the 1D Fourier basis on the unit circle.
    ///
    /// For `Phi = [1, sin(2π h t), cos(2π h t), ...]` we have
    /// `Phi'' = [0, -(2π h)² sin(...), -(2π h)² cos(...), ...]`, i.e.
    /// the second derivative is `-(2π h)² · phi(t)` on each harmonic pair.
    fn second_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array4<f64>, String> {
        let n = coords.nrows();
        let d = coords.ncols();
        if d != 1 {
            return Err(format!(
                "PeriodicHarmonicEvaluator::second_jet: expected latent_dim == 1, got {d}"
            ));
        }
        let m = self.num_basis;
        let num_harmonics = (m - 1) / 2;
        let two_pi = 2.0 * std::f64::consts::PI;
        let mut h = Array4::<f64>::zeros((n, m, 1, 1));
        for row in 0..n {
            let t = coords[[row, 0]];
            for k in 1..=num_harmonics {
                let freq = two_pi * (k as f64);
                let freq2 = freq * freq;
                let angle = freq * t;
                let s = angle.sin();
                let c = angle.cos();
                let s_idx = 2 * k - 1;
                let c_idx = 2 * k;
                h[[row, s_idx, 0, 0]] = -freq2 * s;
                h[[row, c_idx, 0, 0]] = -freq2 * c;
            }
        }
        Ok(h)
    }
}

impl SaeBasisThirdJet for PeriodicHarmonicEvaluator {
    /// Third derivative of the 1-D Fourier basis on the unit circle.
    ///
    /// For `Phi = [1, sin(2π h t), cos(2π h t), …]` the chain of derivatives is
    /// `sin → ωc → −ω²s → −ω³c` and `cos → −ωs → −ω²c → ω³s`, so the third
    /// derivative is `[0, −(2π h)³ cos(…), +(2π h)³ sin(…), …]`.
    fn third_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array5<f64>, String> {
        let n = coords.nrows();
        let d = coords.ncols();
        if d != 1 {
            return Err(format!(
                "PeriodicHarmonicEvaluator::third_jet: expected latent_dim == 1, got {d}"
            ));
        }
        let m = self.num_basis;
        let num_harmonics = (m - 1) / 2;
        let two_pi = 2.0 * std::f64::consts::PI;
        let mut t3 = Array5::<f64>::zeros((n, m, 1, 1, 1));
        for row in 0..n {
            let t = coords[[row, 0]];
            for k in 1..=num_harmonics {
                let freq = two_pi * (k as f64);
                let freq3 = freq * freq * freq;
                let angle = freq * t;
                let s = angle.sin();
                let c = angle.cos();
                let s_idx = 2 * k - 1;
                let c_idx = 2 * k;
                t3[[row, s_idx, 0, 0, 0]] = -freq3 * c;
                t3[[row, c_idx, 0, 0, 0]] = freq3 * s;
            }
        }
        Ok(t3)
    }
}

/// Raw-angle periodic evaluator for the minimal SAE-manifold front-end.
///
/// The basis is exactly `[cos(t), sin(t)]` with `t` measured in radians. If
/// the latent coordinate has more than one axis, the first axis carries the
/// circle phase and the remaining axes are left available to the optimizer but
/// do not enter this basis.
#[derive(Debug, Clone)]
pub struct RawPeriodicCircleEvaluator {
    pub latent_dim: usize,
}

impl RawPeriodicCircleEvaluator {
    pub fn new(latent_dim: usize) -> Result<Self, String> {
        if latent_dim == 0 {
            return Err("RawPeriodicCircleEvaluator requires latent_dim >= 1".to_string());
        }
        Ok(Self { latent_dim })
    }
}

impl SaeBasisEvaluator for RawPeriodicCircleEvaluator {
    fn phi_eta_split(&self, n_basis: usize) -> Result<PhiEtaSplit, String> {
        if n_basis != 2 {
            return Err(format!(
                "RawPeriodicCircleEvaluator::phi_eta_split: n_basis {n_basis} != 2"
            ));
        }
        // Both raw circle columns `[cos, sin]` are base (η-invariant) columns:
        // the homotopy has nothing to dial because the base topology IS the full
        // basis. These columns embed a curved unit circle, so this is a
        // base-topology relaxation, not a linear one.
        Ok(PhiEtaSplit::all_base(n_basis))
    }

    fn second_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array4<f64>, String>> {
        if coords.ncols() != self.latent_dim {
            return Some(Err(format!(
                "RawPeriodicCircleEvaluator::second_jet_dyn: expected latent_dim {}, got {}",
                self.latent_dim,
                coords.ncols()
            )));
        }
        None
    }

    fn third_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array5<f64>, String>> {
        if coords.ncols() != self.latent_dim {
            return Some(Err(format!(
                "RawPeriodicCircleEvaluator::third_jet_dyn: expected latent_dim {}, got {}",
                self.latent_dim,
                coords.ncols()
            )));
        }
        None
    }

    fn evaluate(&self, coords: ArrayView2<'_, f64>) -> Result<(Array2<f64>, Array3<f64>), String> {
        if coords.ncols() != self.latent_dim {
            return Err(format!(
                "RawPeriodicCircleEvaluator: expected latent_dim {}, got {}",
                self.latent_dim,
                coords.ncols()
            ));
        }
        let n = coords.nrows();
        let mut phi = Array2::<f64>::zeros((n, 2));
        let mut jet = Array3::<f64>::zeros((n, 2, self.latent_dim));
        for row in 0..n {
            let t = coords[[row, 0]];
            phi[[row, 0]] = t.cos();
            phi[[row, 1]] = t.sin();
            jet[[row, 0, 0]] = -t.sin();
            jet[[row, 1, 0]] = t.cos();
        }
        Ok((phi, jet))
    }
}

/// Diagonal of the chart-local seven-column sphere basis penalty.
///
/// The columns are `[1, x, y, z, xy, yz, xz]`; the constant column carries a
/// numerically-negligible ridge (`1e-8`) so the penalty stays positive
/// definite, the three degree-one columns `[x, y, z]` are penalized at unit
/// weight, and the three bilinear columns at weight `4` (their second-order
/// angular content).
/// This is the single source of truth for the chart penalty shared between the
/// core SAE path and the PyFFI `sphere_chart_basis_with_jet` helper.
pub const SPHERE_CHART_PENALTY_DIAGONAL: [f64; 7] = [1e-8, 1.0, 1.0, 1.0, 4.0, 4.0, 4.0];

/// Shared single source of truth for the chart-local sphere basis and its
/// analytic first-derivative (lat/lon) jet.
///
/// `coords` is an `(N, 2)` array of latitude/longitude pairs in radians. The
/// returned `phi` has shape `(N, 7)` with columns `[1, x, y, z, xy, yz, xz]`
/// for the unit-sphere embedding `x = cos(lat)cos(lon)`, `y = cos(lat)sin(lon)`,
/// `z = sin(lat)`; the returned `jet` has shape `(N, 7, 2)` with the last axis
/// indexing `[∂/∂lat, ∂/∂lon]`.
///
/// NOTE: this is a lat/lon *product* chart, not an intrinsic / rotation-
/// invariant `S²` parametrization. Two caveats follow from the chart:
///   * **Pole gauge degeneracy.** At the poles `cos(lat) = 0`, so
///     `x = y = xy = yz = xz = 0` and `z = ±1` for *every* longitude — all
///     longitudes collapse to the same physical point (longitude is a gauge
///     coordinate there). The longitude jet vanishes at the poles
///     (`∂x/∂lon = ∂y/∂lon = 0`), so the Hessian carries a longitude gauge
///     degeneracy while the ARD von-Mises prior can still prefer a longitude —
///     i.e. the fit is not invariant under pole equivalence.
///   * **Not rotationally invariant.** The quadratic block `[xy, yz, xz]` is
///     not a rotation-invariant degree-2 spherical-harmonic basis: rotating
///     `xy` produces `x² − y²`, which lies outside its span. A proper
///     atlas / Wahba-style basis is a follow-up; do not read this as intrinsic.
///
/// Within a fixed chart the map and its jet are everywhere `C^∞` in
/// `(lat, lon)`: every column is a
/// polynomial in `cos`/`sin` of the two coordinates, and `cos`/`sin` are entire,
/// so the exact analytic derivatives `∂x/∂lat = -sin(lat)cos(lon)`, … are
/// globally smooth. Latitude is therefore **not** clamped and the latitude
/// derivatives are **not** gated here.
///
/// The physical `lat ∈ [-π/2, π/2]` box that pins a canonical latitude range is
/// enforced where it belongs — in the latent retraction / tangent projection
/// ([`gam_terms::latent::LatentManifold::Interval`]), which clamps the
/// coordinate after each step and zeroes only the *outward-normal* component of
/// the tangent velocity at an active bound (a correct KKT projection). The old
/// binary `chain_lat` gate instead zeroed the *entire* latitude jet at the
/// boundary, making the basis nonsmooth there: an atom whose latitude reached
/// `±π/2` saw a zero latitude gradient and froze, even for the tangential
/// (in-box) direction along which the loss does decrease. Computing the exact
/// jet here and letting the retraction handle the bound restores a smooth
/// objective and the correct boundary behaviour. Both the core path
/// ([`SphereChartEvaluator`]) and the PyFFI helper route through this function.
pub fn sphere_chart_basis_jet(
    coords: ArrayView2<'_, f64>,
) -> Result<(Array2<f64>, Array3<f64>), String> {
    if coords.ncols() != 2 {
        return Err(format!(
            "sphere_chart_basis_jet expects latent_dim == 2, got {}",
            coords.ncols()
        ));
    }
    let n = coords.nrows();
    let mut phi = Array2::<f64>::zeros((n, 7));
    let mut jet = Array3::<f64>::zeros((n, 7, 2));
    for row in 0..n {
        let lat = coords[[row, 0]];
        let lon = coords[[row, 1]];
        let clat = lat.cos();
        let slat = lat.sin();
        let clon = lon.cos();
        let slon = lon.sin();
        let x = clat * clon;
        let y = clat * slon;
        let z = slat;
        phi[[row, 0]] = 1.0;
        phi[[row, 1]] = x;
        phi[[row, 2]] = y;
        phi[[row, 3]] = z;
        phi[[row, 4]] = x * y;
        phi[[row, 5]] = y * z;
        phi[[row, 6]] = x * z;

        let dx_dlat = -slat * clon;
        let dx_dlon = -clat * slon;
        let dy_dlat = -slat * slon;
        let dy_dlon = clat * clon;
        let dz_dlat = clat;
        jet[[row, 1, 0]] = dx_dlat;
        jet[[row, 1, 1]] = dx_dlon;
        jet[[row, 2, 0]] = dy_dlat;
        jet[[row, 2, 1]] = dy_dlon;
        jet[[row, 3, 0]] = dz_dlat;
        jet[[row, 4, 0]] = dx_dlat * y + x * dy_dlat;
        jet[[row, 4, 1]] = dx_dlon * y + x * dy_dlon;
        jet[[row, 5, 0]] = dy_dlat * z + y * dz_dlat;
        jet[[row, 5, 1]] = dy_dlon * z;
        jet[[row, 6, 0]] = dx_dlat * z + x * dz_dlat;
        jet[[row, 6, 1]] = dx_dlon * z;
    }
    Ok((phi, jet))
}

/// Lat/lon sphere chart evaluator used by the Rust-owned minimal SAE path.
#[derive(Debug, Clone)]
pub struct SphereChartEvaluator;

impl SaeBasisEvaluator for SphereChartEvaluator {
    fn phi_eta_split(&self, n_basis: usize) -> Result<PhiEtaSplit, String> {
        if n_basis != 7 {
            return Err(format!(
                "SphereChartEvaluator::phi_eta_split: n_basis {n_basis} != 7"
            ));
        }
        // Base (η-invariant) columns `[1, x, y, z]`; the quadratic cross-terms
        // `[xy, yz, xz]` scale with `eta`. The `[x, y, z]` block is NOT linear in
        // the chart — evaluated on the sphere it traces the unit sphere (the base
        // topology). So `eta = 0` is the base-topology relaxation, not an affine
        // model.
        let mut curved = vec![false; n_basis];
        for col in 4..7 {
            curved[col] = true;
        }
        Ok(PhiEtaSplit::from_curved_mask(curved))
    }

    fn second_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array4<f64>, String>> {
        Some(<Self as SaeBasisSecondJet>::second_jet(self, coords))
    }

    fn third_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array5<f64>, String>> {
        Some(<Self as SaeBasisThirdJet>::third_jet(self, coords))
    }

    fn evaluate(&self, coords: ArrayView2<'_, f64>) -> Result<(Array2<f64>, Array3<f64>), String> {
        sphere_chart_basis_jet(coords)
    }
}

impl SaeBasisSecondJet for SphereChartEvaluator {
    /// Analytic Hessian of the 7-column lat/lon sphere chart basis.
    ///
    /// With `x = cos(lat) cos(lon)`, `y = cos(lat) sin(lon)`, `z = sin(lat)`
    /// the non-trivial second derivatives are
    ///
    /// ```text
    /// x_{lat,lat} = -x,     x_{lon,lon} = -x,     x_{lat,lon} = sin(lat)·sin(lon)
    /// y_{lat,lat} = -y,     y_{lon,lon} = -y,     y_{lat,lon} = -sin(lat)·cos(lon)
    /// z_{lat,lat} = -z,     z_{lon,lon} =  0,     z_{lat,lon} =  0
    /// ```
    ///
    /// Bilinear basis entries `xy, yz, xz` follow the product rule
    /// `(fg)_{αβ} = f_{αβ} g + f_α g_β + f_β g_α + f g_{αβ}`. The map is `C^∞`
    /// in `(lat, lon)`, so the Hessian is the exact analytic one with no clamp
    /// or boundary gating; the `lat ∈ [-π/2, π/2]` box is enforced by the
    /// retraction, not by truncating derivatives (see [`sphere_chart_basis_jet`]).
    fn second_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array4<f64>, String> {
        if coords.ncols() != 2 {
            return Err(format!(
                "SphereChartEvaluator::second_jet expects latent_dim == 2, got {}",
                coords.ncols()
            ));
        }
        let n = coords.nrows();
        let mut h = Array4::<f64>::zeros((n, 7, 2, 2));
        for row in 0..n {
            let lat = coords[[row, 0]];
            let lon = coords[[row, 1]];
            let clat = lat.cos();
            let slat = lat.sin();
            let clon = lon.cos();
            let slon = lon.sin();
            let x = clat * clon;
            let y = clat * slon;
            let z = slat;
            let dx = [-slat * clon, -clat * slon];
            let dy = [-slat * slon, clat * clon];
            let dz = [clat, 0.0];
            let hx = [[-x, slat * slon], [slat * slon, -x]];
            let hy = [[-y, -slat * clon], [-slat * clon, -y]];
            let hz = [[-z, 0.0], [0.0, 0.0]];
            for axis_a in 0..2 {
                for axis_b in 0..2 {
                    h[[row, 1, axis_a, axis_b]] = hx[axis_a][axis_b];
                    h[[row, 2, axis_a, axis_b]] = hy[axis_a][axis_b];
                    h[[row, 3, axis_a, axis_b]] = hz[axis_a][axis_b];
                }
            }
            let pair = |hf: [[f64; 2]; 2],
                        df: [f64; 2],
                        f: f64,
                        hg: [[f64; 2]; 2],
                        dg: [f64; 2],
                        g: f64|
             -> [[f64; 2]; 2] {
                let mut out = [[0.0; 2]; 2];
                for axis_a in 0..2 {
                    for axis_b in 0..2 {
                        out[axis_a][axis_b] = hf[axis_a][axis_b] * g
                            + df[axis_a] * dg[axis_b]
                            + df[axis_b] * dg[axis_a]
                            + f * hg[axis_a][axis_b];
                    }
                }
                out
            };
            let hxy = pair(hx, dx, x, hy, dy, y);
            let hyz = pair(hy, dy, y, hz, dz, z);
            let hxz = pair(hx, dx, x, hz, dz, z);
            for axis_a in 0..2 {
                for axis_b in 0..2 {
                    h[[row, 4, axis_a, axis_b]] = hxy[axis_a][axis_b];
                    h[[row, 5, axis_a, axis_b]] = hyz[axis_a][axis_b];
                    h[[row, 6, axis_a, axis_b]] = hxz[axis_a][axis_b];
                }
            }
        }
        Ok(h)
    }
}

impl SaeBasisThirdJet for SphereChartEvaluator {
    /// Third derivative of the 7-column lat/lon sphere chart basis
    /// `[1, x, y, z, xy, yz, xz]`.
    ///
    /// Each Cartesian coordinate is *separable* in (lat, lon):
    /// `x = cos(lat) cos(lon)`, `y = cos(lat) sin(lon)`, `z = sin(lat)·1`. A
    /// separable coordinate's mixed derivative is the product of the per-axis
    /// derivative of the right order, so it is fully described by two
    /// length-4 derivative tables (orders 0..3) — one per axis. The map is
    /// `C^∞` in `(lat, lon)`; the tables are the exact analytic derivatives
    /// with no clamp or boundary gating (the `lat ∈ [-π/2, π/2]` box is
    /// enforced by the retraction, see [`sphere_chart_basis_jet`]).
    ///
    /// The bilinear columns `xy, yz, xz` are products of two separable
    /// coordinates; their third derivative is the symmetric triple-Leibniz sum
    /// over the `2³` ways to route the three derivative operators to the two
    /// factors. This is the order-3 generalization of the `pair` Leibniz used
    /// in [`SaeBasisSecondJet::second_jet`], so the two stay structurally
    /// identical and a finite difference of `second_jet` pins it.
    fn third_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array5<f64>, String> {
        if coords.ncols() != 2 {
            return Err(format!(
                "SphereChartEvaluator::third_jet expects latent_dim == 2, got {}",
                coords.ncols()
            ));
        }
        let n = coords.nrows();
        let mut t3 = Array5::<f64>::zeros((n, 7, 2, 2, 2));
        // Derivative of a separable coordinate along axes `ax` (each 0 = lat,
        // 1 = lon): product of the lat table at order `#lat` and the lon table
        // at order `#lon`.
        let single = |lat: &[f64; 4], lon: &[f64; 4], ax: [usize; 3]| -> f64 {
            let n_lat = ax.iter().filter(|&&q| q == 0).count();
            lat[n_lat] * lon[3 - n_lat]
        };
        // Third derivative of a product of two separable coordinates: sum over
        // all 2³ routings of the three operators to factor f vs g (Leibniz).
        let product = |f_lat: &[f64; 4],
                       f_lon: &[f64; 4],
                       g_lat: &[f64; 4],
                       g_lon: &[f64; 4],
                       ax: [usize; 3]|
         -> f64 {
            let mut acc = 0.0;
            for mask in 0u8..8 {
                let (mut f_lat_n, mut f_lon_n, mut g_lat_n, mut g_lon_n) = (0, 0, 0, 0);
                for (i, &axis) in ax.iter().enumerate() {
                    let to_f = (mask >> i) & 1 == 1;
                    match (to_f, axis == 0) {
                        (true, true) => f_lat_n += 1,
                        (true, false) => f_lon_n += 1,
                        (false, true) => g_lat_n += 1,
                        (false, false) => g_lon_n += 1,
                    }
                }
                acc += f_lat[f_lat_n] * f_lon[f_lon_n] * g_lat[g_lat_n] * g_lon[g_lon_n];
            }
            acc
        };
        for row in 0..n {
            let lat = coords[[row, 0]];
            let lon = coords[[row, 1]];
            let clat = lat.cos();
            let slat = lat.sin();
            let clon = lon.cos();
            let slon = lon.sin();
            // Per-axis derivative tables, orders 0..3 (exact analytic, no clamp).
            let cos_lat = [clat, -slat, -clat, slat];
            let sin_lat = [slat, clat, -slat, -clat];
            let cos_lon = [clon, -slon, -clon, slon];
            let sin_lon = [slon, clon, -slon, -clon];
            let const_lon = [1.0, 0.0, 0.0, 0.0];
            // x = cos(lat)cos(lon), y = cos(lat)sin(lon), z = sin(lat).
            let (x_lat, x_lon) = (&cos_lat, &cos_lon);
            let (y_lat, y_lon) = (&cos_lat, &sin_lon);
            let (z_lat, z_lon) = (&sin_lat, &const_lon);
            for axis_a in 0..2 {
                for axis_b in 0..2 {
                    for axis_c in 0..2 {
                        let ax = [axis_a, axis_b, axis_c];
                        t3[[row, 1, axis_a, axis_b, axis_c]] = single(x_lat, x_lon, ax);
                        t3[[row, 2, axis_a, axis_b, axis_c]] = single(y_lat, y_lon, ax);
                        t3[[row, 3, axis_a, axis_b, axis_c]] = single(z_lat, z_lon, ax);
                        t3[[row, 4, axis_a, axis_b, axis_c]] =
                            product(x_lat, x_lon, y_lat, y_lon, ax);
                        t3[[row, 5, axis_a, axis_b, axis_c]] =
                            product(y_lat, y_lon, z_lat, z_lon, ax);
                        t3[[row, 6, axis_a, axis_b, axis_c]] =
                            product(x_lat, x_lon, z_lat, z_lon, ax);
                    }
                }
            }
        }
        Ok(t3)
    }
}

/// One real spherical-harmonic column `Y_l^m`, fixed at construction.
///
/// Each column factors as `R_{l,m}(lat) · T_m(lon)` — a pure-latitude
/// amplitude times a pure-longitude phase — so its full jet in `(lat, lon)` is
/// the separable outer product of the two per-axis derivative tables. The
/// latitude amplitude is `R_{l,m}(lat) = N_{l,m} · cos(lat)^{|m|} · Q(sin lat)`,
/// where `Q = d^{|m|}/du^{|m|} P_l(u)` is the `|m|`-th `u`-derivative of the
/// Legendre polynomial (the associated-Legendre function's polynomial part) and
/// `N_{l,m}` is the orthonormalization constant. The `(-1)^m` Condon–Shortley
/// phase is intentionally dropped: it is a per-column sign the decoder absorbs
/// and it has no effect on the span, orthonormality, or reconstruction.
#[derive(Debug, Clone)]
struct SphHarmonicColumn {
    /// Spherical-harmonic degree `l`.
    degree: usize,
    /// Signed order `m ∈ [-l, l]`: `m ≥ 0` uses `cos(m·lon)`, `m < 0` uses
    /// `sin(|m|·lon)`.
    m: i64,
    /// `|m|`, the `cos(lat)` power and the Legendre `u`-derivative order.
    am: usize,
    /// Orthonormalization constant `N_{l,m}` (includes the `√2` for `m ≠ 0`).
    norm: f64,
    /// `Q = d^{|m|}/du^{|m|} P_l(u)`, ascending powers of `u = sin(lat)`.
    assoc: Vec<f64>,
    /// `true` iff `l ≥ 2` — the curved (η-dialed) refinement above the base
    /// monopole+dipole sphere embedding.
    curved: bool,
}

/// Real orthonormal spherical-harmonic evaluator on `S^2`, charted by
/// `(lat, lon)` with `x = cos(lat)cos(lon)`, `y = cos(lat)sin(lon)`,
/// `z = sin(lat)` (matching [`SphereChartEvaluator`]).
///
/// This is the rotation-covariant basis the fixed 7-column
/// [`SphereChartEvaluator`] is not: its columns are the `(degree+1)²` real
/// spherical harmonics `Y_l^m`, `l = 0..=degree`, `m = -l..=l`, an orthonormal
/// basis for every band-limited field on the sphere. The degree-2 chart spans
/// only `[1, x, y, z, xy, yz, xz]` — the monopole, the dipole, and *three of the
/// five* quadrupoles — so it cannot represent `x²−y²`, `3z²−1`, or any `l ≥ 3`
/// content; a genuinely sphere-class field with higher-degree structure is
/// out of its span. This evaluator closes that resolution gap, with the working
/// degree chosen by [`select_spherical_harmonic_degree`] under the same
/// spectral-noise-floor bandwidth doctrine the torus/circle already use.
///
/// **Smoothness / poles.** Every column is a polynomial in `sin(lat)`,
/// `cos(lat)`, `sin(m·lon)`, `cos(m·lon)`, all entire, so the map and its jets
/// are globally `C^∞` in `(lat, lon)` — latitude is never clamped here. The
/// pole gauge degeneracy is intrinsic to the lat/lon chart, not this basis: at
/// `cos(lat) = 0` every `m ≠ 0` harmonic and its longitude derivative vanish, so
/// all longitudes collapse to one physical point. That is handled exactly where
/// [`SphereChartEvaluator`] handles it — by the retraction / tangent projection
/// enforcing the `lat ∈ [-π/2, π/2]` box and the pole-seam seeding (issue
/// #1890) — never by truncating a derivative.
#[derive(Debug, Clone)]
pub struct SphericalHarmonicEvaluator {
    degree: usize,
    columns: Vec<SphHarmonicColumn>,
}

/// Spectral identity of one real spherical-harmonic column.
///
/// The evaluator's public column order is triangular in `(degree, order)`:
/// degree `l` contributes orders `-l..=l`.  Carrying this identity beside the
/// analytic evaluator lets quotient constructions restrict the cover's exact
/// eigenspaces without re-deriving indexing conventions.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct SphericalHarmonicMode {
    pub degree: usize,
    pub order: i64,
    /// Positive-Laplacian eigenvalue `l(l+1)` on the unit sphere.
    pub laplace_eigenvalue: f64,
    /// Diagonal `L²(S²)` Gram entry.  The real harmonics are orthonormal.
    pub l2_gram_weight: f64,
}

/// Legendre polynomials `P_0..=P_degree` as ascending-power coefficient vectors,
/// via the exact three-term recurrence `l·P_l = (2l-1)·u·P_{l-1} − (l-1)·P_{l-2}`.
fn legendre_polynomials(degree: usize) -> Vec<Vec<f64>> {
    let mut polys: Vec<Vec<f64>> = Vec::with_capacity(degree + 1);
    polys.push(vec![1.0]);
    if degree >= 1 {
        polys.push(vec![0.0, 1.0]);
    }
    for l in 2..=degree {
        let prev = &polys[l - 1];
        let prev2 = &polys[l - 2];
        let mut next = vec![0.0_f64; l + 1];
        // (2l-1)·u·P_{l-1}
        for (k, &coeff) in prev.iter().enumerate() {
            next[k + 1] += (2 * l - 1) as f64 * coeff;
        }
        // −(l-1)·P_{l-2}
        for (k, &coeff) in prev2.iter().enumerate() {
            next[k] -= (l - 1) as f64 * coeff;
        }
        let inv = 1.0 / l as f64;
        for c in next.iter_mut() {
            *c *= inv;
        }
        polys.push(next);
    }
    polys
}

/// `order`-th derivative of an ascending-power polynomial (exact).
fn polynomial_derivative(coeffs: &[f64], order: usize) -> Vec<f64> {
    let mut c = coeffs.to_vec();
    for _ in 0..order {
        if c.len() <= 1 {
            return vec![0.0];
        }
        c = (1..c.len()).map(|k| k as f64 * c[k]).collect();
    }
    c
}

/// Orthonormalization constant `N_{l,m} = √((2l+1)/(4π) · (l−|m|)!/(l+|m|)!)`,
/// times `√2` for `m ≠ 0` (the real-harmonic combination of `±m`).
fn spherical_harmonic_norm(l: usize, m: i64) -> f64 {
    let am = m.unsigned_abs() as usize;
    let mut ratio = 1.0_f64; // (l-am)! / (l+am)!
    for k in (l - am + 1)..=(l + am) {
        ratio /= k as f64;
    }
    let base = ((2 * l + 1) as f64 / (4.0 * std::f64::consts::PI) * ratio).sqrt();
    if m == 0 {
        base
    } else {
        base * std::f64::consts::SQRT_2
    }
}

/// One-dimensional order-3 jet `[f, f', f'', f''']` product (Leibniz).
#[inline]
fn sph_jet_mul(a: [f64; 4], b: [f64; 4]) -> [f64; 4] {
    [
        a[0] * b[0],
        a[1] * b[0] + a[0] * b[1],
        a[2] * b[0] + 2.0 * a[1] * b[1] + a[0] * b[2],
        a[3] * b[0] + 3.0 * a[2] * b[1] + 3.0 * a[1] * b[2] + a[0] * b[3],
    ]
}

impl SphericalHarmonicEvaluator {
    pub fn new(degree: usize) -> Result<Self, String> {
        let side = degree.checked_add(1).ok_or_else(|| {
            "SphericalHarmonicEvaluator: basis width overflowed usize".to_string()
        })?;
        let basis_size = side.checked_mul(side).ok_or_else(|| {
            "SphericalHarmonicEvaluator: basis width overflowed usize".to_string()
        })?;
        let polys = legendre_polynomials(degree);
        let mut columns = Vec::with_capacity(basis_size);
        for l in 0..=degree {
            for m in -(l as i64)..=(l as i64) {
                let am = m.unsigned_abs() as usize;
                columns.push(SphHarmonicColumn {
                    degree: l,
                    m,
                    am,
                    norm: spherical_harmonic_norm(l, m),
                    assoc: polynomial_derivative(&polys[l], am),
                    curved: l >= 2,
                });
            }
        }
        Ok(Self { degree, columns })
    }

    pub fn degree(&self) -> usize {
        self.degree
    }

    pub fn basis_size(&self) -> usize {
        self.columns.len()
    }

    /// Exact `(l,m)` and Laplace-eigenvalue metadata in evaluator column order.
    pub fn spectral_modes(&self) -> Vec<SphericalHarmonicMode> {
        self.columns
            .iter()
            .map(|column| {
                let degree = column.degree as f64;
                SphericalHarmonicMode {
                    degree: column.degree,
                    order: column.m,
                    laplace_eigenvalue: degree * (degree + 1.0),
                    l2_gram_weight: 1.0,
                }
            })
            .collect()
    }

    /// Latitude amplitude jet `[R, R', R'', R''']` of a column via order-3 jet
    /// arithmetic in `lat`: `R = N · cos(lat)^{|m|} · Q(sin lat)`.
    fn lat_table(&self, col: &SphHarmonicColumn, lat: f64) -> [f64; 4] {
        let (s, c) = lat.sin_cos();
        let slat = [s, c, -s, -c];
        let clat = [c, -s, -c, s];
        let mut pow = [1.0, 0.0, 0.0, 0.0];
        for _ in 0..col.am {
            pow = sph_jet_mul(pow, clat);
        }
        // Horner evaluation of Q at slat.
        let mut acc = [0.0, 0.0, 0.0, 0.0];
        for &coeff in col.assoc.iter().rev() {
            acc = sph_jet_mul(acc, slat);
            acc[0] += coeff;
        }
        let mut r = sph_jet_mul(pow, acc);
        for value in r.iter_mut() {
            *value *= col.norm;
        }
        r
    }

    /// Longitude phase jet `[T, T', T'', T''']`: `cos(m·lon)` for `m ≥ 0`,
    /// `sin(|m|·lon)` for `m < 0`.
    fn lon_table(m: i64, lon: f64) -> [f64; 4] {
        if m == 0 {
            return [1.0, 0.0, 0.0, 0.0];
        }
        let mf = m.unsigned_abs() as f64;
        let (s, c) = (mf * lon).sin_cos();
        if m > 0 {
            [c, -mf * s, -mf * mf * c, mf * mf * mf * s]
        } else {
            [s, mf * c, -mf * mf * s, -mf * mf * mf * c]
        }
    }

    fn check_coords(&self, coords: ArrayView2<'_, f64>, what: &str) -> Result<(), String> {
        if coords.ncols() != 2 {
            return Err(format!(
                "SphericalHarmonicEvaluator::{what}: expected latent_dim == 2 (lat, lon), got {}",
                coords.ncols()
            ));
        }
        Ok(())
    }
}

impl SaeBasisEvaluator for SphericalHarmonicEvaluator {
    fn phi_eta_split(&self, n_basis: usize) -> Result<PhiEtaSplit, String> {
        let expected = self.basis_size();
        if n_basis != expected {
            return Err(format!(
                "SphericalHarmonicEvaluator::phi_eta_split: n_basis {n_basis} != evaluator width {expected}"
            ));
        }
        // Base (η-invariant) block: monopole + dipole `l ≤ 1` (the base sphere
        // embedding). The quadrupole and higher harmonics `l ≥ 2` are the
        // η-dialed curvature refinement, mirroring the fixed sphere chart.
        let curved = self
            .columns
            .iter()
            .map(|col| col.curved)
            .collect::<Vec<_>>();
        Ok(PhiEtaSplit::from_curved_mask(curved))
    }

    /// The `(l, m)` triangular indexing is not a clean per-axis tensor product.
    fn factor_basis_sizes(&self) -> Option<(usize, usize)> {
        None
    }

    fn second_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array4<f64>, String>> {
        Some(<Self as SaeBasisSecondJet>::second_jet(self, coords))
    }

    fn third_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array5<f64>, String>> {
        Some(<Self as SaeBasisThirdJet>::third_jet(self, coords))
    }

    fn evaluate(&self, coords: ArrayView2<'_, f64>) -> Result<(Array2<f64>, Array3<f64>), String> {
        self.check_coords(coords, "evaluate")?;
        let n = coords.nrows();
        let m = self.basis_size();
        let mut phi = Array2::<f64>::zeros((n, m));
        let mut jet = Array3::<f64>::zeros((n, m, 2));
        for row in 0..n {
            let lat = coords[[row, 0]];
            let lon = coords[[row, 1]];
            for (col_idx, col) in self.columns.iter().enumerate() {
                let r = self.lat_table(col, lat);
                let t = Self::lon_table(col.m, lon);
                phi[[row, col_idx]] = r[0] * t[0];
                jet[[row, col_idx, 0]] = r[1] * t[0];
                jet[[row, col_idx, 1]] = r[0] * t[1];
            }
        }
        Ok((phi, jet))
    }
}

impl SaeBasisSecondJet for SphericalHarmonicEvaluator {
    fn second_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array4<f64>, String> {
        self.check_coords(coords, "second_jet")?;
        let n = coords.nrows();
        let m = self.basis_size();
        let mut h = Array4::<f64>::zeros((n, m, 2, 2));
        for row in 0..n {
            let lat = coords[[row, 0]];
            let lon = coords[[row, 1]];
            for (col_idx, col) in self.columns.iter().enumerate() {
                let r = self.lat_table(col, lat);
                let t = Self::lon_table(col.m, lon);
                // Separable column: ∂^{p+q}/∂lat^p ∂lon^q = R^{(p)} · T^{(q)}.
                h[[row, col_idx, 0, 0]] = r[2] * t[0];
                h[[row, col_idx, 0, 1]] = r[1] * t[1];
                h[[row, col_idx, 1, 0]] = r[1] * t[1];
                h[[row, col_idx, 1, 1]] = r[0] * t[2];
            }
        }
        Ok(h)
    }
}

impl SaeBasisThirdJet for SphericalHarmonicEvaluator {
    fn third_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array5<f64>, String> {
        self.check_coords(coords, "third_jet")?;
        let n = coords.nrows();
        let m = self.basis_size();
        let mut t3 = Array5::<f64>::zeros((n, m, 2, 2, 2));
        for row in 0..n {
            let lat = coords[[row, 0]];
            let lon = coords[[row, 1]];
            for (col_idx, col) in self.columns.iter().enumerate() {
                let r = self.lat_table(col, lat);
                let t = Self::lon_table(col.m, lon);
                for axis_a in 0..2 {
                    for axis_b in 0..2 {
                        for axis_c in 0..2 {
                            // Separable: the mixed derivative depends only on how
                            // many of the three operators are latitude (axis 0).
                            let n_lat = (axis_a == 0) as usize
                                + (axis_b == 0) as usize
                                + (axis_c == 0) as usize;
                            t3[[row, col_idx, axis_a, axis_b, axis_c]] = r[n_lat] * t[3 - n_lat];
                        }
                    }
                }
            }
        }
        Ok(t3)
    }
}

/// Select the working spherical-harmonic degree by the spectral-noise-floor
/// bandwidth doctrine: fit the full `max_degree` basis on the chart, form the
/// per-degree mean coefficient energy `E_l = ‖coefficients at degree l‖² /
/// (2l+1)`, and cut at the largest cliff in that spectrum — the degree above
/// which the energy collapses to the noise floor. This is the sphere analogue
/// of the torus/circle bandwidth selection and the Duchon-sheet resolution
/// ladder: no tuned resolution constant, the data's own spectrum sets the band.
///
/// Returns the smallest sufficient degree in `0..=max_degree`. If the spectrum
/// has no pronounced floor (every degree carries signal), the full `max_degree`
/// is kept.
pub fn select_spherical_harmonic_degree(
    coords: ArrayView2<'_, f64>,
    target: ArrayView2<'_, f64>,
    max_degree: usize,
) -> Result<usize, String> {
    use faer::Side;
    use gam_linalg::faer_ndarray::{FaerCholesky, fast_ata, fast_atb};

    if coords.nrows() != target.nrows() {
        return Err(format!(
            "select_spherical_harmonic_degree: coords rows {} != target rows {}",
            coords.nrows(),
            target.nrows()
        ));
    }
    if max_degree == 0 {
        return Ok(0);
    }
    let evaluator = SphericalHarmonicEvaluator::new(max_degree)?;
    let (phi, _) = evaluator.evaluate(coords)?;
    // Ridge-stabilized normal-equations fit of the full basis to the target.
    let mut gram = fast_ata(&phi);
    let scale = gram.diag().iter().copied().fold(0.0_f64, f64::max);
    let ridge = if scale > 0.0 { scale } else { 1.0 } * 1e-9;
    for d in gram.diag_mut().iter_mut() {
        *d += ridge;
    }
    let rhs = fast_atb(&phi, &target.to_owned());
    let decoder = gram
        .cholesky(Side::Lower)
        .map_err(|err| {
            format!("select_spherical_harmonic_degree: gram factorization failed: {err:?}")
        })?
        .solve_mat(&rhs);

    // Per-degree mean coefficient energy over the (2l+1) orders and all outputs.
    let mut energy = vec![0.0_f64; max_degree + 1];
    let mut col = 0usize;
    for (l, e) in energy.iter_mut().enumerate() {
        let width = 2 * l + 1;
        let mut acc = 0.0_f64;
        for _ in 0..width {
            for out in 0..decoder.ncols() {
                acc += decoder[[col, out]] * decoder[[col, out]];
            }
            col += 1;
        }
        *e = acc / width as f64;
    }

    // Largest downward cliff `E_l / E_{l+1}`: the boundary between signal and
    // the noise floor. Require a decisive factor (100×) before trusting it as a
    // floor; otherwise keep the full band.
    let mut cut = max_degree;
    let mut best_ratio = 1.0_f64;
    for l in 0..max_degree {
        let hi = energy[l];
        let lo = energy[l + 1].max(f64::MIN_POSITIVE);
        let ratio = hi / lo;
        if ratio > best_ratio {
            best_ratio = ratio;
            cut = l;
        }
    }
    if best_ratio < 100.0 {
        cut = max_degree;
    }
    Ok(cut)
}

/// Tensor-product periodic harmonic evaluator for a `d`-dimensional torus
/// `T^d = (S^1)^d`. The basis is the tensor product over each axis of the
/// 1-D circle basis, stored **sine-first** within each harmonic (matching the
/// `evaluate` layout below, `s_idx = 2h-1`, `c_idx = 2h`, and the
/// `PeriodicHarmonicEvaluator` convention):
/// `[1, sin(2π·1·t), cos(2π·1·t), …, sin(2π·H·t), cos(2π·H·t)]`
/// (each axis contributes `2H+1` factors, so the total basis size is
/// `(2H+1)^d`). The latent coords are angular phases in `[0, 1)` (consistent
/// with the periodic 1-D atoms).
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum RealHarmonicComponent {
    Constant,
    Sine { harmonic: usize },
    Cosine { harmonic: usize },
}

impl RealHarmonicComponent {
    pub fn harmonic(self) -> usize {
        match self {
            Self::Constant => 0,
            Self::Sine { harmonic } | Self::Cosine { harmonic } => harmonic,
        }
    }

    /// Sign under reflection of the corresponding circle coordinate.
    pub fn reflection_sign(self) -> i8 {
        match self {
            Self::Sine { .. } => -1,
            Self::Constant | Self::Cosine { .. } => 1,
        }
    }
}

/// Spectral identity of one tensor-product real Fourier column.
#[derive(Debug, Clone, PartialEq)]
pub struct TorusHarmonicMode {
    pub components: Vec<RealHarmonicComponent>,
    /// Positive-Laplacian eigenvalue `sum_a k_a^2`, up to the common `(2π)^2`
    /// chart-unit scale absorbed by the smoothing strength.
    pub laplace_eigenvalue: f64,
    /// Diagonal Gram entry under normalized Haar measure on the torus.  Each
    /// sine/cosine factor contributes `1/2`; constant factors contribute one.
    pub l2_gram_weight: f64,
}

#[derive(Debug, Clone)]
pub struct TorusHarmonicEvaluator {
    latent_dim: usize,
    num_harmonics: usize,
    axis_basis_size: usize,
    basis_size: usize,
}

impl TorusHarmonicEvaluator {
    pub fn new(latent_dim: usize, num_harmonics: usize) -> Result<Self, String> {
        if latent_dim == 0 {
            return Err("TorusHarmonicEvaluator requires latent_dim >= 1".to_string());
        }
        if num_harmonics == 0 {
            return Err("TorusHarmonicEvaluator requires num_harmonics >= 1".to_string());
        }
        let axis_width = num_harmonics
            .checked_mul(2)
            .and_then(|twice| twice.checked_add(1))
            .ok_or_else(|| {
                "TorusHarmonicEvaluator: per-axis basis width overflowed usize".to_string()
            })?;
        let basis_size = (0..latent_dim)
            .try_fold(1usize, |width, _| width.checked_mul(axis_width))
            .ok_or_else(|| {
                "TorusHarmonicEvaluator: tensor basis width overflowed usize".to_string()
            })?;
        Ok(Self {
            latent_dim,
            num_harmonics,
            axis_basis_size: axis_width,
            basis_size,
        })
    }

    pub fn latent_dim(&self) -> usize {
        self.latent_dim
    }

    pub fn num_harmonics(&self) -> usize {
        self.num_harmonics
    }

    pub fn axis_basis_size(&self) -> usize {
        self.axis_basis_size
    }

    pub fn basis_size(&self) -> usize {
        self.basis_size
    }

    /// Real Fourier component represented by one per-axis column index.
    pub fn axis_component(axis_column: usize) -> RealHarmonicComponent {
        if axis_column == 0 {
            RealHarmonicComponent::Constant
        } else {
            let harmonic = axis_column.div_ceil(2);
            if axis_column % 2 == 1 {
                RealHarmonicComponent::Sine { harmonic }
            } else {
                RealHarmonicComponent::Cosine { harmonic }
            }
        }
    }

    /// Exact tensor Fourier modes and Laplace eigenvalues in evaluator column
    /// order (last axis fastest).
    pub fn spectral_modes(&self) -> Vec<TorusHarmonicMode> {
        let axis_m = self.axis_basis_size();
        let mut index = vec![0usize; self.latent_dim];
        let mut modes = Vec::with_capacity(self.basis_size());
        for _ in 0..self.basis_size() {
            let components: Vec<RealHarmonicComponent> =
                index.iter().copied().map(Self::axis_component).collect();
            let laplace_eigenvalue = components
                .iter()
                .map(|component| {
                    let harmonic = component.harmonic() as f64;
                    harmonic * harmonic
                })
                .sum();
            let l2_gram_weight = components.iter().fold(1.0, |weight, component| {
                if *component == RealHarmonicComponent::Constant {
                    weight
                } else {
                    0.5 * weight
                }
            });
            modes.push(TorusHarmonicMode {
                components,
                laplace_eigenvalue,
                l2_gram_weight,
            });
            for axis in (0..self.latent_dim).rev() {
                index[axis] += 1;
                if index[axis] < axis_m {
                    break;
                }
                index[axis] = 0;
            }
        }
        modes
    }
}

impl SaeBasisEvaluator for TorusHarmonicEvaluator {
    fn phi_eta_split(&self, n_basis: usize) -> Result<PhiEtaSplit, String> {
        let expected = self.basis_size();
        if n_basis != expected {
            return Err(format!(
                "TorusHarmonicEvaluator::phi_eta_split: n_basis {n_basis} != evaluator width {expected}"
            ));
        }
        // Base (η-invariant) columns: the constant plus each single-axis
        // first-harmonic `[sin, cos]` (one non-constant axis, harmonic ≤ 1).
        // These embed the per-axis circles — the torus base topology — so they
        // are curved, not linear; only higher harmonics and cross-axis products
        // scale with `eta`. `eta = 0` is thus the base-topology relaxation.
        let d = self.latent_dim;
        let axis_m = self.axis_basis_size();
        let mut curved = Vec::with_capacity(n_basis);
        let mut idx = vec![0usize; d];
        for _flat in 0..n_basis {
            let mut nonconstant_axes = 0usize;
            let mut has_higher_harmonic = false;
            for &axis_col in &idx {
                if axis_col > 0 {
                    nonconstant_axes += 1;
                    if axis_col > 2 {
                        has_higher_harmonic = true;
                    }
                }
            }
            curved.push(has_higher_harmonic || nonconstant_axes > 1);
            for axis in (0..d).rev() {
                idx[axis] += 1;
                if idx[axis] < axis_m {
                    break;
                }
                idx[axis] = 0;
            }
        }
        Ok(PhiEtaSplit::from_curved_mask(curved))
    }

    /// A `d = 2` torus is the tensor product of two equal circle factors, each
    /// of width `axis_basis_size()` (constant-leading: column 0 ≡ 1). Higher-
    /// dimensional tori are not a TWO-factor product, so they report `None`.
    fn factor_basis_sizes(&self) -> Option<(usize, usize)> {
        if self.latent_dim == 2 {
            let m = self.axis_basis_size();
            Some((m, m))
        } else {
            None
        }
    }

    fn second_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array4<f64>, String>> {
        Some(<Self as SaeBasisSecondJet>::second_jet(self, coords))
    }

    fn third_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array5<f64>, String>> {
        Some(<Self as SaeBasisThirdJet>::third_jet(self, coords))
    }

    fn evaluate(&self, coords: ArrayView2<'_, f64>) -> Result<(Array2<f64>, Array3<f64>), String> {
        // Single source of truth: allocate correctly-shaped buffers and fill
        // them through the in-place path (see `evaluate_into`).
        let n = coords.nrows();
        let m = self.basis_size();
        let d = self.latent_dim;
        let mut phi = Array2::<f64>::zeros((n, m));
        let mut jet = Array3::<f64>::zeros((n, m, d));
        self.evaluate_into(&mut phi, &mut jet, coords)?;
        Ok((phi, jet))
    }

    fn evaluate_into(
        &self,
        phi: &mut Array2<f64>,
        jet: &mut Array3<f64>,
        coords: ArrayView2<'_, f64>,
    ) -> Result<(), String> {
        let d = self.latent_dim;
        if coords.ncols() != d {
            return Err(format!(
                "TorusHarmonicEvaluator: expected latent_dim {d}, got {}",
                coords.ncols()
            ));
        }
        let n = coords.nrows();
        let axis_m = self.axis_basis_size();
        let m = self.basis_size();
        if phi.dim() != (n, m) {
            return Err(format!(
                "TorusHarmonicEvaluator::evaluate_into: Φ buffer {:?} != ({n}, {m})",
                phi.dim()
            ));
        }
        if jet.dim() != (n, m, d) {
            return Err(format!(
                "TorusHarmonicEvaluator::evaluate_into: jet buffer {:?} != ({n}, {m}, {d})",
                jet.dim()
            ));
        }
        let h_max = self.num_harmonics;
        let two_pi = 2.0 * std::f64::consts::PI;
        // Every `(row, flat)` Φ entry and every `(row, flat, axis)` jet entry is
        // written unconditionally in the product loops below, so no pre-clear is
        // needed — stale workspace contents cannot survive.
        // Per-axis evaluation buffer: phi_axis[axis][col] and dphi_axis[axis][col].
        let mut phi_axis = vec![vec![0.0_f64; axis_m]; d];
        let mut dphi_axis = vec![vec![0.0_f64; axis_m]; d];
        for row in 0..n {
            for axis in 0..d {
                let t = coords[[row, axis]];
                phi_axis[axis][0] = 1.0;
                dphi_axis[axis][0] = 0.0;
                for h in 1..=h_max {
                    let freq = two_pi * (h as f64);
                    let angle = freq * t;
                    let s = angle.sin();
                    let c = angle.cos();
                    let s_idx = 2 * h - 1;
                    let c_idx = 2 * h;
                    phi_axis[axis][s_idx] = s;
                    phi_axis[axis][c_idx] = c;
                    dphi_axis[axis][s_idx] = freq * c;
                    dphi_axis[axis][c_idx] = -freq * s;
                }
            }
            // Enumerate the Cartesian product of per-axis indices in
            // lexicographic order (axis 0 is the slowest).
            let mut idx = vec![0usize; d];
            for flat in 0..m {
                let mut val = 1.0_f64;
                for axis in 0..d {
                    val *= phi_axis[axis][idx[axis]];
                }
                phi[[row, flat]] = val;
                // ∂/∂coords[row, axis_target] = product over axes, replacing
                // phi_axis[axis_target] with its derivative.
                for axis_target in 0..d {
                    let mut deriv = 1.0_f64;
                    for axis in 0..d {
                        deriv *= if axis == axis_target {
                            dphi_axis[axis][idx[axis]]
                        } else {
                            phi_axis[axis][idx[axis]]
                        };
                    }
                    jet[[row, flat, axis_target]] = deriv;
                }
                // Increment lexicographic index (last axis fastest).
                for axis in (0..d).rev() {
                    idx[axis] += 1;
                    if idx[axis] < axis_m {
                        break;
                    }
                    idx[axis] = 0;
                }
            }
        }
        Ok(())
    }
}

impl SaeBasisSecondJet for TorusHarmonicEvaluator {
    /// Hessian of the tensor-product torus basis.
    ///
    /// Each basis function factors as `Φ_flat = Π_axis f_axis(t_axis)`, so
    ///
    /// * `∂² Φ / ∂t_a ∂t_b = (Π_{k ∉ {a, b}} f_k) · f_a'(t_a) · f_b'(t_b)`
    ///   when `a ≠ b`,
    /// * `∂² Φ / ∂t_a²    = (Π_{k ≠ a} f_k) · f_a''(t_a)` on the diagonal.
    ///
    /// Per-axis the basis is `[1, sin(2π h t), cos(2π h t), …]`, so
    /// `f_axis''(t) = -(2π h)² · f_axis(t)` on the harmonic columns and 0 on
    /// the constant column.
    fn second_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array4<f64>, String> {
        let d = self.latent_dim;
        if coords.ncols() != d {
            return Err(format!(
                "TorusHarmonicEvaluator::second_jet expects latent_dim == {d}, got {}",
                coords.ncols()
            ));
        }
        let n = coords.nrows();
        let axis_m = self.axis_basis_size();
        let m = self.basis_size();
        let h_max = self.num_harmonics;
        let two_pi = 2.0 * std::f64::consts::PI;
        let mut hess = Array4::<f64>::zeros((n, m, d, d));
        let mut phi_axis = vec![vec![0.0_f64; axis_m]; d];
        let mut dphi_axis = vec![vec![0.0_f64; axis_m]; d];
        let mut d2phi_axis = vec![vec![0.0_f64; axis_m]; d];
        for row in 0..n {
            for axis in 0..d {
                let t = coords[[row, axis]];
                phi_axis[axis][0] = 1.0;
                dphi_axis[axis][0] = 0.0;
                d2phi_axis[axis][0] = 0.0;
                for k in 1..=h_max {
                    let freq = two_pi * (k as f64);
                    let freq2 = freq * freq;
                    let angle = freq * t;
                    let s = angle.sin();
                    let c = angle.cos();
                    let s_idx = 2 * k - 1;
                    let c_idx = 2 * k;
                    phi_axis[axis][s_idx] = s;
                    phi_axis[axis][c_idx] = c;
                    dphi_axis[axis][s_idx] = freq * c;
                    dphi_axis[axis][c_idx] = -freq * s;
                    d2phi_axis[axis][s_idx] = -freq2 * s;
                    d2phi_axis[axis][c_idx] = -freq2 * c;
                }
            }
            let mut idx = vec![0usize; d];
            for flat in 0..m {
                for axis_a in 0..d {
                    for axis_b in 0..d {
                        let mut prod = 1.0_f64;
                        for axis in 0..d {
                            let factor = if axis == axis_a && axis == axis_b {
                                d2phi_axis[axis][idx[axis]]
                            } else if axis == axis_a || axis == axis_b {
                                dphi_axis[axis][idx[axis]]
                            } else {
                                phi_axis[axis][idx[axis]]
                            };
                            prod *= factor;
                        }
                        hess[[row, flat, axis_a, axis_b]] = prod;
                    }
                }
                for axis in (0..d).rev() {
                    idx[axis] += 1;
                    if idx[axis] < axis_m {
                        break;
                    }
                    idx[axis] = 0;
                }
            }
        }
        Ok(hess)
    }
}

impl SaeBasisThirdJet for TorusHarmonicEvaluator {
    /// Third derivative of the tensor-product torus basis.
    ///
    /// Each basis function factors as `Φ_flat = Π_axis f_axis(t_axis)`, so its
    /// third derivative `∂³Φ / ∂t_a ∂t_b ∂t_c` is the product, over every
    /// axis, of `f_axis` differentiated as many times as that axis appears in
    /// `{a, b, c}` (0..3). Per axis the basis is `[1, sin(2π h t),
    /// cos(2π h t), …]`, whose order-3 derivative is `[0, −(2π h)³ cos(…),
    /// +(2π h)³ sin(…), …]`. This is the order-3 sibling of
    /// [`SaeBasisSecondJet::second_jet`].
    fn third_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array5<f64>, String> {
        let d = self.latent_dim;
        if coords.ncols() != d {
            return Err(format!(
                "TorusHarmonicEvaluator::third_jet expects latent_dim == {d}, got {}",
                coords.ncols()
            ));
        }
        let n = coords.nrows();
        let axis_m = self.axis_basis_size();
        let m = self.basis_size();
        let h_max = self.num_harmonics;
        let two_pi = 2.0 * std::f64::consts::PI;
        let mut t3 = Array5::<f64>::zeros((n, m, d, d, d));
        // Per-axis derivative tables indexed [axis][order 0..3][column].
        let mut deriv_axis = vec![vec![vec![0.0_f64; axis_m]; 4]; d];
        for row in 0..n {
            for axis in 0..d {
                let t = coords[[row, axis]];
                for order in 0..4 {
                    deriv_axis[axis][order][0] = 0.0;
                }
                deriv_axis[axis][0][0] = 1.0;
                for k in 1..=h_max {
                    let freq = two_pi * (k as f64);
                    let freq2 = freq * freq;
                    let freq3 = freq2 * freq;
                    let angle = freq * t;
                    let s = angle.sin();
                    let c = angle.cos();
                    let s_idx = 2 * k - 1;
                    let c_idx = 2 * k;
                    deriv_axis[axis][0][s_idx] = s;
                    deriv_axis[axis][0][c_idx] = c;
                    deriv_axis[axis][1][s_idx] = freq * c;
                    deriv_axis[axis][1][c_idx] = -freq * s;
                    deriv_axis[axis][2][s_idx] = -freq2 * s;
                    deriv_axis[axis][2][c_idx] = -freq2 * c;
                    deriv_axis[axis][3][s_idx] = -freq3 * c;
                    deriv_axis[axis][3][c_idx] = freq3 * s;
                }
            }
            let mut idx = vec![0usize; d];
            for flat in 0..m {
                for axis_a in 0..d {
                    for axis_b in 0..d {
                        for axis_c in 0..d {
                            let mut prod = 1.0_f64;
                            for axis in 0..d {
                                let order = (axis == axis_a) as usize
                                    + (axis == axis_b) as usize
                                    + (axis == axis_c) as usize;
                                prod *= deriv_axis[axis][order][idx[axis]];
                            }
                            t3[[row, flat, axis_a, axis_b, axis_c]] = prod;
                        }
                    }
                }
                for axis in (0..d).rev() {
                    idx[axis] += 1;
                    if idx[axis] < axis_m {
                        break;
                    }
                    idx[axis] = 0;
                }
            }
        }
        Ok(t3)
    }
}

/// One-dimensional real character of an involutive deck transformation.
///
/// This type is deliberately narrower than a signed-permutation action.  It
/// certifies that one cover column `f` is itself an eigenvector of pullback by
/// the non-identity deck map `g`: `g*f = f` or `g*f = -f`.  The exact Reynolds
/// projector is therefore
///
/// `P_Gamma f = (f + g*f) / 2 = f` for `Trivial`, and `0` for `Sign`.
///
/// A basis that is merely permuted by `g` cannot use this representation: for
/// example, swapping two columns retains their orbit sum, not either column.
/// Such an action requires an explicit restriction matrix `Q`.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum DeckInvolutionCharacter {
    Trivial,
    Sign,
}

impl DeckInvolutionCharacter {
    /// Eigenvalue of the exact group average `(I + g*) / 2` on this column.
    fn reynolds_projector_eigenvalue(self) -> u8 {
        match self {
            Self::Trivial => 1,
            Self::Sign => 0,
        }
    }
}

/// Spectral metadata for every cover column, in cover evaluator order.
///
/// Requiring a full character table, rather than accepting a caller-selected
/// mask, makes the group average the single authority for quotient width,
/// column restriction, Gram entries, penalty eigenvalues, and nullity.
#[derive(Debug, Clone, Copy, PartialEq)]
struct CoverSpectralCharacter {
    deck_character: DeckInvolutionCharacter,
    laplace_eigenvalue: f64,
    l2_gram_weight: f64,
    curved_if_retained: bool,
}

/// Exact diagonal-character restriction of an analytic cover evaluator.
///
/// Static column selection is exact only because each supported cover column
/// is a simultaneous one-dimensional character of the deck involution.  Mere
/// signed permutation of a basis is insufficient; non-diagonal actions require
/// orbit sums or a dense restriction matrix `Q`.  Values and analytic
/// first/second/third jets are restricted from one cover evaluator, while the
/// same complete character table determines width, Gram, spectral penalties,
/// and null space.
#[derive(Debug, Clone)]
pub struct QuotientSpectralEvaluator {
    quotient_name: String,
    cover: Arc<dyn SaeBasisThirdJet>,
    cover_width: usize,
    cover_columns: Vec<usize>,
    laplace_eigenvalues: Vec<f64>,
    l2_gram_weights: Vec<f64>,
    curved_columns: Vec<bool>,
}

/// Number of invariant real spherical harmonics on `RP²` through even cover
/// degree `2H`: `sum_{r=0}^H (4r+1) = (H+1)(2H+1)`.
pub fn projective_plane_basis_size(harmonic_order: usize) -> Result<usize, String> {
    if harmonic_order == 0 {
        return Err("projective_plane_basis_size requires harmonic_order >= 1".to_string());
    }
    harmonic_order
        .checked_add(1)
        .and_then(|left| {
            harmonic_order
                .checked_mul(2)
                .and_then(|twice| twice.checked_add(1))
                .and_then(|right| left.checked_mul(right))
        })
        .ok_or_else(|| "projective_plane_basis_size overflowed usize".to_string())
}

/// Number of invariant real torus harmonics on the flat Klein quotient through
/// equal per-axis order `H`: `1 + 2 floor(H/2) + H + 2H²`.
pub fn klein_bottle_basis_size(num_harmonics: usize) -> Result<usize, String> {
    if num_harmonics == 0 {
        return Err("klein_bottle_basis_size requires num_harmonics >= 1".to_string());
    }
    let cross_width = num_harmonics
        .checked_mul(num_harmonics)
        .and_then(|square| square.checked_mul(2))
        .ok_or_else(|| "klein_bottle_basis_size overflowed usize".to_string())?;
    1usize
        .checked_add(2 * (num_harmonics / 2))
        .and_then(|width| width.checked_add(num_harmonics))
        .and_then(|width| width.checked_add(cross_width))
        .ok_or_else(|| "klein_bottle_basis_size overflowed usize".to_string())
}

impl QuotientSpectralEvaluator {
    /// Build a quotient from a cover and its complete diagonal character table.
    ///
    /// `characters[cover_column]` must certify the pullback eigencharacter of
    /// that exact cover column.  The Reynolds projector derives the retained
    /// columns; callers cannot supply an already-filtered mask.  A connected
    /// closed cover and quotient have one zero-eigenvalue mode (the constant),
    /// whose character is trivial and which remains in the homotopy base.
    fn from_diagonal_involution_characters(
        quotient_name: impl Into<String>,
        cover: Arc<dyn SaeBasisThirdJet>,
        cover_width: usize,
        characters: Vec<CoverSpectralCharacter>,
    ) -> Result<Self, String> {
        let quotient_name = quotient_name.into();
        if quotient_name.trim().is_empty() {
            return Err("QuotientSpectralEvaluator requires a non-empty quotient name".to_string());
        }
        if cover_width == 0 {
            return Err(format!(
                "QuotientSpectralEvaluator[{quotient_name}]: cover width must be positive"
            ));
        }
        if characters.len() != cover_width {
            return Err(format!(
                "QuotientSpectralEvaluator[{quotient_name}]: diagonal character table width {} != cover width {cover_width}",
                characters.len()
            ));
        }

        let mut cover_columns = Vec::with_capacity(cover_width);
        let mut laplace_eigenvalues = Vec::with_capacity(cover_width);
        let mut l2_gram_weights = Vec::with_capacity(cover_width);
        let mut curved_columns = Vec::with_capacity(cover_width);
        let mut nullity = 0usize;
        for (cover_column, mode) in characters.into_iter().enumerate() {
            if !(mode.laplace_eigenvalue.is_finite() && mode.laplace_eigenvalue >= 0.0) {
                return Err(format!(
                    "QuotientSpectralEvaluator[{quotient_name}]: cover column {cover_column} has invalid Laplace eigenvalue {}",
                    mode.laplace_eigenvalue
                ));
            }
            if !(mode.l2_gram_weight.is_finite() && mode.l2_gram_weight > 0.0) {
                return Err(format!(
                    "QuotientSpectralEvaluator[{quotient_name}]: cover column {cover_column} has invalid L2 Gram weight {}",
                    mode.l2_gram_weight
                ));
            }
            if mode.laplace_eigenvalue == 0.0 {
                if mode.deck_character != DeckInvolutionCharacter::Trivial {
                    return Err(format!(
                        "QuotientSpectralEvaluator[{quotient_name}]: the constant cover mode at column {cover_column} must have trivial deck character"
                    ));
                }
                nullity += 1;
                if mode.curved_if_retained {
                    return Err(format!(
                        "QuotientSpectralEvaluator[{quotient_name}]: the constant null mode cannot be curvature-scaled"
                    ));
                }
            }
            if mode.deck_character.reynolds_projector_eigenvalue() == 1 {
                cover_columns.push(cover_column);
                laplace_eigenvalues.push(mode.laplace_eigenvalue);
                l2_gram_weights.push(mode.l2_gram_weight);
                curved_columns.push(mode.curved_if_retained);
            }
        }
        if nullity != 1 {
            return Err(format!(
                "QuotientSpectralEvaluator[{quotient_name}]: connected closed quotient requires exactly one constant null mode; found {nullity}"
            ));
        }

        Ok(Self {
            quotient_name,
            cover,
            cover_width,
            cover_columns,
            laplace_eigenvalues,
            l2_gram_weights,
            curved_columns,
        })
    }

    /// Real harmonics on `RP² = S²/{u ~ -u}` through quotient order `H`.
    ///
    /// Antipodal parity is `Y_lm(-u) = (-1)^l Y_lm(u)`, hence precisely the
    /// even degrees `l = 0, 2, ..., 2H` survive.
    pub fn projective_plane(harmonic_order: usize) -> Result<Self, String> {
        if harmonic_order == 0 {
            return Err(
                "QuotientSpectralEvaluator::projective_plane requires harmonic_order >= 1"
                    .to_string(),
            );
        }
        let max_degree = harmonic_order.checked_mul(2).ok_or_else(|| {
            "QuotientSpectralEvaluator::projective_plane: maximum cover degree overflowed usize"
                .to_string()
        })?;
        let cover_side = max_degree.checked_add(1).ok_or_else(|| {
            "QuotientSpectralEvaluator::projective_plane: cover width overflowed usize".to_string()
        })?;
        cover_side.checked_mul(cover_side).ok_or_else(|| {
            "QuotientSpectralEvaluator::projective_plane: cover width overflowed usize".to_string()
        })?;
        let expected_width = projective_plane_basis_size(harmonic_order)?;

        let cover = SphericalHarmonicEvaluator::new(max_degree)?;
        let cover_width = cover.basis_size();
        let characters = cover
            .spectral_modes()
            .into_iter()
            .map(|mode| CoverSpectralCharacter {
                deck_character: if mode.degree % 2 == 0 {
                    DeckInvolutionCharacter::Trivial
                } else {
                    DeckInvolutionCharacter::Sign
                },
                laplace_eigenvalue: mode.laplace_eigenvalue,
                l2_gram_weight: mode.l2_gram_weight,
                // `l = 0,2` contains the constant and Veronese embedding; only
                // higher even degrees are the curvature refinement.
                curved_if_retained: mode.degree > 2,
            })
            .collect::<Vec<_>>();
        let evaluator = Self::from_diagonal_involution_characters(
            "projective-plane",
            Arc::new(cover),
            cover_width,
            characters,
        )?;
        if evaluator.basis_size() != expected_width {
            return Err(format!(
                "QuotientSpectralEvaluator::projective_plane: group average produced width {}, expected {expected_width}",
                evaluator.basis_size()
            ));
        }
        Ok(evaluator)
    }

    /// Real harmonics on the flat Klein bottle
    /// `T²/{(theta, phi) ~ (theta + 1/2, -phi)}` through order `H` per axis.
    /// `H >= 2` is required because the standard smooth `R⁴` Klein embedding
    /// uses theta harmonics one and two; its constant plus six coordinate modes
    /// form the seven-column homotopy base.
    pub fn klein_bottle(num_harmonics: usize) -> Result<Self, String> {
        if num_harmonics < 2 {
            return Err(
                "QuotientSpectralEvaluator::klein_bottle requires num_harmonics >= 2 for the standard R4 embedding"
                    .to_string(),
            );
        }
        let expected_width = klein_bottle_basis_size(num_harmonics)?;

        let cover = TorusHarmonicEvaluator::new(2, num_harmonics)?;
        let cover_width = cover.basis_size();
        let mut characters = Vec::with_capacity(cover_width);
        for (cover_column, mode) in cover.spectral_modes().into_iter().enumerate() {
            let [theta_component, phi_component] = mode.components.as_slice() else {
                return Err(format!(
                    "QuotientSpectralEvaluator::klein_bottle: torus cover mode {cover_column} did not have two factors"
                ));
            };
            let theta_harmonic = theta_component.harmonic();
            let phi_harmonic = phi_component.harmonic();
            let half_turn_sign = if theta_harmonic % 2 == 0 { 1 } else { -1 };
            let deck_character = if half_turn_sign * phi_component.reflection_sign() == 1 {
                DeckInvolutionCharacter::Trivial
            } else {
                DeckInvolutionCharacter::Sign
            };
            // Standard R4 Klein embedding in this unit-period chart:
            //   {1,
            //    sin/cos(4πθ),
            //    sin/cos(4πθ) cos(2πφ),
            //    sin/cos(2πθ) sin(2πφ)}.
            // These are exactly one constant plus six coordinate columns.
            let embedding_mode = (theta_harmonic == 0 && phi_harmonic == 0)
                || (theta_harmonic == 2 && phi_harmonic == 0)
                || (theta_harmonic == 2
                    && phi_harmonic == 1
                    && phi_component.reflection_sign() == 1)
                || (theta_harmonic == 1
                    && phi_harmonic == 1
                    && phi_component.reflection_sign() == -1);
            characters.push(CoverSpectralCharacter {
                deck_character,
                laplace_eigenvalue: mode.laplace_eigenvalue,
                l2_gram_weight: mode.l2_gram_weight,
                curved_if_retained: !embedding_mode,
            });
        }
        let evaluator = Self::from_diagonal_involution_characters(
            "klein-bottle",
            Arc::new(cover),
            cover_width,
            characters,
        )?;
        if evaluator.basis_size() != expected_width {
            return Err(format!(
                "QuotientSpectralEvaluator::klein_bottle: group average produced width {}, expected {expected_width}",
                evaluator.basis_size()
            ));
        }
        Ok(evaluator)
    }

    pub fn quotient_name(&self) -> &str {
        &self.quotient_name
    }

    pub fn basis_size(&self) -> usize {
        self.cover_columns.len()
    }

    pub fn cover_width(&self) -> usize {
        self.cover_width
    }

    pub fn cover_columns(&self) -> &[usize] {
        &self.cover_columns
    }

    pub fn laplace_eigenvalues(&self) -> &[f64] {
        &self.laplace_eigenvalues
    }

    pub fn l2_gram_weights(&self) -> &[f64] {
        &self.l2_gram_weights
    }

    /// Exact diagonal `L²` function-space Gram inherited from the cover.
    pub fn function_space_gram(&self) -> Array2<f64> {
        let mut gram = Array2::<f64>::zeros((self.basis_size(), self.basis_size()));
        for (column, &weight) in self.l2_gram_weights.iter().enumerate() {
            gram[[column, column]] = weight;
        }
        gram
    }

    /// Exact spectral penalty `G · diag(lambda_laplace^power)`.
    ///
    /// The Gram factor matters for the unnormalized real Fourier columns: it
    /// makes the generalized penalty eigenvalues exactly the cover Laplacian's
    /// eigenvalues rather than silently changing them by column convention.
    pub fn spectral_penalty(&self, power: u32) -> Result<Array2<f64>, String> {
        if power == 0 {
            return Err(format!(
                "QuotientSpectralEvaluator[{}]::spectral_penalty requires power >= 1",
                self.quotient_name
            ));
        }
        let exponent = i32::try_from(power).map_err(|_| {
            format!(
                "QuotientSpectralEvaluator[{}]::spectral_penalty power {power} exceeds i32::MAX",
                self.quotient_name
            )
        })?;
        let mut penalty = Array2::<f64>::zeros((self.basis_size(), self.basis_size()));
        for column in 0..self.basis_size() {
            let value =
                self.l2_gram_weights[column] * self.laplace_eigenvalues[column].powi(exponent);
            if !value.is_finite() {
                return Err(format!(
                    "QuotientSpectralEvaluator[{}]::spectral_penalty overflowed at quotient column {column}",
                    self.quotient_name
                ));
            }
            penalty[[column, column]] = value;
        }
        Ok(penalty)
    }

    pub fn nullspace_dimension(&self) -> usize {
        self.laplace_eigenvalues
            .iter()
            .filter(|&&eigenvalue| eigenvalue == 0.0)
            .count()
    }

    fn validate_cover_value_jet(
        &self,
        phi: &Array2<f64>,
        jet: &Array3<f64>,
        n_rows: usize,
        latent_dim: usize,
    ) -> Result<(), String> {
        if phi.dim() != (n_rows, self.cover_width) {
            return Err(format!(
                "QuotientSpectralEvaluator[{}]: cover Phi shape {:?} != ({n_rows}, {})",
                self.quotient_name,
                phi.dim(),
                self.cover_width
            ));
        }
        if jet.dim() != (n_rows, self.cover_width, latent_dim) {
            return Err(format!(
                "QuotientSpectralEvaluator[{}]: cover jet shape {:?} != ({n_rows}, {}, {latent_dim})",
                self.quotient_name,
                jet.dim(),
                self.cover_width
            ));
        }
        Ok(())
    }
}

impl SaeBasisEvaluator for QuotientSpectralEvaluator {
    fn evaluate(&self, coords: ArrayView2<'_, f64>) -> Result<(Array2<f64>, Array3<f64>), String> {
        let n_rows = coords.nrows();
        let latent_dim = coords.ncols();
        let mut phi = Array2::<f64>::zeros((n_rows, self.basis_size()));
        let mut jet = Array3::<f64>::zeros((n_rows, self.basis_size(), latent_dim));
        self.evaluate_into(&mut phi, &mut jet, coords)?;
        Ok((phi, jet))
    }

    fn evaluate_into(
        &self,
        phi: &mut Array2<f64>,
        jet: &mut Array3<f64>,
        coords: ArrayView2<'_, f64>,
    ) -> Result<(), String> {
        let n_rows = coords.nrows();
        let latent_dim = coords.ncols();
        let quotient_width = self.basis_size();
        if phi.dim() != (n_rows, quotient_width) {
            return Err(format!(
                "QuotientSpectralEvaluator[{}]::evaluate_into: Phi buffer {:?} != ({n_rows}, {quotient_width})",
                self.quotient_name,
                phi.dim()
            ));
        }
        if jet.dim() != (n_rows, quotient_width, latent_dim) {
            return Err(format!(
                "QuotientSpectralEvaluator[{}]::evaluate_into: jet buffer {:?} != ({n_rows}, {quotient_width}, {latent_dim})",
                self.quotient_name,
                jet.dim()
            ));
        }

        let (cover_phi, cover_jet) = self.cover.evaluate(coords)?;
        self.validate_cover_value_jet(&cover_phi, &cover_jet, n_rows, latent_dim)?;
        for (quotient_column, &cover_column) in self.cover_columns.iter().enumerate() {
            for row in 0..n_rows {
                phi[[row, quotient_column]] = cover_phi[[row, cover_column]];
                for axis in 0..latent_dim {
                    jet[[row, quotient_column, axis]] = cover_jet[[row, cover_column, axis]];
                }
            }
        }
        Ok(())
    }

    fn phi_eta_split(&self, n_basis: usize) -> Result<PhiEtaSplit, String> {
        if n_basis != self.basis_size() {
            return Err(format!(
                "QuotientSpectralEvaluator[{}]::phi_eta_split: n_basis {n_basis} != evaluator width {}",
                self.quotient_name,
                self.basis_size()
            ));
        }
        Ok(PhiEtaSplit::from_curved_mask(self.curved_columns.clone()))
    }

    /// A quotient mask couples cover factors and is not a tensor-product basis.
    fn factor_basis_sizes(&self) -> Option<(usize, usize)> {
        None
    }

    fn second_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array4<f64>, String>> {
        Some(<Self as SaeBasisSecondJet>::second_jet(self, coords))
    }

    fn third_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array5<f64>, String>> {
        Some(<Self as SaeBasisThirdJet>::third_jet(self, coords))
    }
}

impl SaeBasisSecondJet for QuotientSpectralEvaluator {
    fn second_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array4<f64>, String> {
        let n_rows = coords.nrows();
        let latent_dim = coords.ncols();
        let cover_hessian = self.cover.second_jet(coords)?;
        if cover_hessian.dim() != (n_rows, self.cover_width, latent_dim, latent_dim) {
            return Err(format!(
                "QuotientSpectralEvaluator[{}]: cover second-jet shape {:?} != ({n_rows}, {}, {latent_dim}, {latent_dim})",
                self.quotient_name,
                cover_hessian.dim(),
                self.cover_width
            ));
        }
        let mut hessian = Array4::<f64>::zeros((n_rows, self.basis_size(), latent_dim, latent_dim));
        for (quotient_column, &cover_column) in self.cover_columns.iter().enumerate() {
            for row in 0..n_rows {
                for axis_a in 0..latent_dim {
                    for axis_b in 0..latent_dim {
                        hessian[[row, quotient_column, axis_a, axis_b]] =
                            cover_hessian[[row, cover_column, axis_a, axis_b]];
                    }
                }
            }
        }
        Ok(hessian)
    }
}

impl SaeBasisThirdJet for QuotientSpectralEvaluator {
    fn third_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array5<f64>, String> {
        let n_rows = coords.nrows();
        let latent_dim = coords.ncols();
        let cover_third = self.cover.third_jet(coords)?;
        if cover_third.dim() != (n_rows, self.cover_width, latent_dim, latent_dim, latent_dim) {
            return Err(format!(
                "QuotientSpectralEvaluator[{}]: cover third-jet shape {:?} != ({n_rows}, {}, {latent_dim}, {latent_dim}, {latent_dim})",
                self.quotient_name,
                cover_third.dim(),
                self.cover_width
            ));
        }
        let mut third = Array5::<f64>::zeros((
            n_rows,
            self.basis_size(),
            latent_dim,
            latent_dim,
            latent_dim,
        ));
        for (quotient_column, &cover_column) in self.cover_columns.iter().enumerate() {
            for row in 0..n_rows {
                for axis_a in 0..latent_dim {
                    for axis_b in 0..latent_dim {
                        for axis_c in 0..latent_dim {
                            third[[row, quotient_column, axis_a, axis_b, axis_c]] =
                                cover_third[[row, cover_column, axis_a, axis_b, axis_c]];
                        }
                    }
                }
            }
        }
        Ok(third)
    }
}

/// Affine Euclidean/Duchon fallback for the minimal fit entrypoint.
#[derive(Debug, Clone)]
pub struct AffineCoordinateEvaluator {
    pub latent_dim: usize,
}

impl AffineCoordinateEvaluator {
    pub fn new(latent_dim: usize) -> Self {
        Self { latent_dim }
    }
}

impl SaeBasisEvaluator for AffineCoordinateEvaluator {
    fn phi_eta_split(&self, n_basis: usize) -> Result<PhiEtaSplit, String> {
        let expected = self.latent_dim + 1;
        if n_basis != expected {
            return Err(format!(
                "AffineCoordinateEvaluator::phi_eta_split: n_basis {n_basis} != {expected}"
            ));
        }
        // The affine `[1, x, y, z]` basis is genuinely linear, so every column is
        // a base column and the base-topology relaxation coincides with a true
        // linear model here (the one basis where the two notions agree).
        Ok(PhiEtaSplit::all_base(n_basis))
    }

    fn second_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array4<f64>, String>> {
        Some(<Self as SaeBasisSecondJet>::second_jet(self, coords))
    }

    fn third_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array5<f64>, String>> {
        Some(<Self as SaeBasisThirdJet>::third_jet(self, coords))
    }

    fn evaluate(&self, coords: ArrayView2<'_, f64>) -> Result<(Array2<f64>, Array3<f64>), String> {
        if coords.ncols() != self.latent_dim {
            return Err(format!(
                "AffineCoordinateEvaluator: expected latent_dim {}, got {}",
                self.latent_dim,
                coords.ncols()
            ));
        }
        let n = coords.nrows();
        let m = self.latent_dim + 1;
        let mut phi = Array2::<f64>::zeros((n, m));
        let mut jet = Array3::<f64>::zeros((n, m, self.latent_dim));
        phi.column_mut(0).fill(1.0);
        for row in 0..n {
            for axis in 0..self.latent_dim {
                phi[[row, axis + 1]] = coords[[row, axis]];
                jet[[row, axis + 1, axis]] = 1.0;
            }
        }
        Ok((phi, jet))
    }
}

impl SaeBasisSecondJet for AffineCoordinateEvaluator {
    /// Second derivative of the affine basis `[1, t_1, ..., t_d]`.
    ///
    /// Every basis function is at most linear in `t`, so all second derivatives
    /// are identically zero. Returns the all-zeros tensor of shape
    /// `(n_obs, d+1, d, d)`.
    fn second_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array4<f64>, String> {
        if coords.ncols() != self.latent_dim {
            return Err(format!(
                "AffineCoordinateEvaluator::second_jet: expected latent_dim {}, got {}",
                self.latent_dim,
                coords.ncols()
            ));
        }
        let n = coords.nrows();
        let m = self.latent_dim + 1;
        let d = self.latent_dim;
        Ok(Array4::<f64>::zeros((n, m, d, d)))
    }
}

impl SaeBasisThirdJet for AffineCoordinateEvaluator {
    /// Third derivative of the affine basis `[1, t_1, …, t_d]`. Every column is
    /// at most linear, so all third derivatives vanish identically.
    fn third_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array5<f64>, String> {
        if coords.ncols() != self.latent_dim {
            return Err(format!(
                "AffineCoordinateEvaluator::third_jet: expected latent_dim {}, got {}",
                self.latent_dim,
                coords.ncols()
            ));
        }
        let n = coords.nrows();
        let m = self.latent_dim + 1;
        let d = self.latent_dim;
        Ok(Array5::<f64>::zeros((n, m, d, d, d)))
    }
}

/// Scale-free Duchon atom evaluator for the SAE-manifold Newton loop.
///
/// Recomputes the radial+polynomial design `Φ(t)` and its first/second
/// input-location jets at arbitrary latent coordinates against a fixed set of
/// `centers` and Duchon null-space `order`. The column layout — the
/// kernel block `Φ_radial(t)·Z` followed by the polynomial block `P(t)`,
/// both carrying the same scalar kernel amplification `α` — matches
/// [`gam_terms::basis::build_duchon_basis`] under the SAE atom's spec
/// (`length_scale = None`, `power = 0`, no identifiability transform). The
/// forward design and the jet are produced from a single core entry point
/// ([`gam_terms::basis::duchon_sae_atom_basis_with_jet`]) so they always agree on
/// column count and scaling — the exact contract issue #247 pinned.
#[derive(Debug, Clone)]
pub struct DuchonCoordinateEvaluator {
    pub centers: Array2<f64>,
    pub order: gam_terms::basis::DuchonNullspaceOrder,
}

impl DuchonCoordinateEvaluator {
    /// Build from the atom's centers and Duchon `m` (`m = 1` → constant
    /// null space, `m = 2` → constant+linear, `m = k+1` → degree-`k`).
    pub fn new(centers: Array2<f64>, m: usize) -> Result<Self, String> {
        if centers.ncols() == 0 {
            return Err("DuchonCoordinateEvaluator: centers must have at least one column".into());
        }
        if m == 0 {
            return Err("DuchonCoordinateEvaluator: Duchon m must be at least 1".into());
        }
        let order = match m {
            1 => gam_terms::basis::DuchonNullspaceOrder::Zero,
            2 => gam_terms::basis::DuchonNullspaceOrder::Linear,
            other => gam_terms::basis::DuchonNullspaceOrder::Degree(other - 1),
        };
        Ok(Self { centers, order })
    }
}

impl SaeBasisEvaluator for DuchonCoordinateEvaluator {
    fn affine_transformed_evaluator(
        &self,
        shift: &[f64],
        scale: &[f64],
        n_basis: usize,
    ) -> Result<Option<Arc<dyn SaeBasisSecondJet>>, String> {
        let dim = self.centers.ncols();
        if shift.len() != dim || scale.len() != dim {
            return Err(format!(
                "DuchonCoordinateEvaluator::affine_transformed_evaluator: affine vectors must have length {dim}; got shift={} scale={}",
                shift.len(),
                scale.len()
            ));
        }
        if n_basis == usize::MAX {
            return Err(
                "DuchonCoordinateEvaluator::affine_transformed_evaluator: unreachable basis width"
                    .to_string(),
            );
        }
        if dim != 1 {
            return Ok(None);
        }
        if !(scale[0].is_finite() && scale[0] > 0.0 && shift[0].is_finite()) {
            return Ok(None);
        }
        let mut centers = self.centers.clone();
        for row in 0..centers.nrows() {
            centers[[row, 0]] = (centers[[row, 0]] - shift[0]) / scale[0];
        }
        Ok(Some(Arc::new(Self {
            centers,
            order: self.order,
        })))
    }

    fn phi_eta_split(&self, n_basis: usize) -> Result<PhiEtaSplit, String> {
        let dim = self.centers.ncols();
        let effective = duchon_effective_order_for_eta(self.centers.view(), self.order);
        let n_poly = duchon_polynomial_column_count(dim, effective);
        if n_basis < n_poly {
            return Err(format!(
                "DuchonCoordinateEvaluator::phi_eta_split: n_basis {n_basis} smaller than polynomial block {n_poly}"
            ));
        }
        let n_kernel = n_basis - n_poly;
        let mut curved = vec![false; n_basis];
        for col in 0..n_kernel {
            curved[col] = true;
        }
        if let gam_terms::basis::DuchonNullspaceOrder::Degree(degree) = effective {
            let linear_mask = monomial_linear_mask(dim, degree);
            if linear_mask.len() != n_poly {
                return Err(format!(
                    "DuchonCoordinateEvaluator::phi_eta_split: polynomial mask width {} != {n_poly}",
                    linear_mask.len()
                ));
            }
            for (local_col, linear) in linear_mask.into_iter().enumerate() {
                if !linear {
                    curved[n_kernel + local_col] = true;
                }
            }
        }
        Ok(PhiEtaSplit::from_curved_mask(curved))
    }

    fn second_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array4<f64>, String>> {
        Some(<Self as SaeBasisSecondJet>::second_jet(self, coords))
    }

    fn third_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array5<f64>, String>> {
        Some(<Self as SaeBasisThirdJet>::third_jet(self, coords))
    }

    fn evaluate(&self, coords: ArrayView2<'_, f64>) -> Result<(Array2<f64>, Array3<f64>), String> {
        if coords.ncols() != self.centers.ncols() {
            return Err(format!(
                "DuchonCoordinateEvaluator: expected latent_dim {}, got {}",
                self.centers.ncols(),
                coords.ncols()
            ));
        }
        gam_terms::basis::duchon_sae_atom_basis_with_jet(coords, self.centers.view(), self.order)
            .map_err(|err| err.to_string())
    }
}

impl SaeBasisSecondJet for DuchonCoordinateEvaluator {
    fn second_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array4<f64>, String> {
        if coords.ncols() != self.centers.ncols() {
            return Err(format!(
                "DuchonCoordinateEvaluator::second_jet: expected latent_dim {}, got {}",
                self.centers.ncols(),
                coords.ncols()
            ));
        }
        gam_terms::basis::duchon_sae_atom_second_jet(coords, self.centers.view(), self.order)
            .map_err(|err| err.to_string())
    }
}

impl SaeBasisThirdJet for DuchonCoordinateEvaluator {
    fn third_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array5<f64>, String> {
        if coords.ncols() != self.centers.ncols() {
            return Err(format!(
                "DuchonCoordinateEvaluator::third_jet: expected latent_dim {}, got {}",
                self.centers.ncols(),
                coords.ncols()
            ));
        }
        gam_terms::basis::duchon_sae_atom_third_jet(coords, self.centers.view(), self.order)
            .map_err(|err| err.to_string())
    }
}

/// Flat Euclidean tangent-patch evaluator for the SAE-manifold Newton loop.
///
/// The basis is the set of monomials of total degree ≤ `max_degree` in the
/// atom's latent coordinates (a zero-curvature polynomial expansion, distinct
/// from the thin-plate Duchon kernel). It recomputes the monomial design and
/// its first/second derivatives at arbitrary coordinates, so the inner Newton
/// latent update stays consistent with the deployed design.
#[derive(Debug, Clone)]
pub struct EuclideanPatchEvaluator {
    pub latent_dim: usize,
    pub max_degree: usize,
}

impl EuclideanPatchEvaluator {
    pub fn new(latent_dim: usize, max_degree: usize) -> Result<Self, String> {
        if latent_dim == 0 {
            return Err("EuclideanPatchEvaluator: latent_dim must be positive".into());
        }
        Ok(Self {
            latent_dim,
            max_degree,
        })
    }

    pub fn basis_size(&self) -> usize {
        gam_terms::basis::monomial_exponents(self.latent_dim, self.max_degree).len()
    }
}

impl SaeBasisEvaluator for EuclideanPatchEvaluator {
    fn affine_transformed_evaluator(
        &self,
        shift: &[f64],
        scale: &[f64],
        n_basis: usize,
    ) -> Result<Option<Arc<dyn SaeBasisSecondJet>>, String> {
        if shift.len() != self.latent_dim || scale.len() != self.latent_dim {
            return Err(format!(
                "EuclideanPatchEvaluator::affine_transformed_evaluator: affine vectors must have length {}; got shift={} scale={}",
                self.latent_dim,
                shift.len(),
                scale.len()
            ));
        }
        if n_basis != self.basis_size() {
            return Err(format!(
                "EuclideanPatchEvaluator::affine_transformed_evaluator: n_basis {n_basis} != evaluator width {}",
                self.basis_size()
            ));
        }
        if shift.iter().chain(scale.iter()).any(|v| !v.is_finite())
            || scale.iter().any(|&v| v <= 0.0)
        {
            return Ok(None);
        }
        Ok(Some(Arc::new(Self {
            latent_dim: self.latent_dim,
            max_degree: self.max_degree,
        })))
    }

    fn phi_eta_split(&self, n_basis: usize) -> Result<PhiEtaSplit, String> {
        let linear_mask = monomial_linear_mask(self.latent_dim, self.max_degree);
        if linear_mask.len() != n_basis {
            return Err(format!(
                "EuclideanPatchEvaluator::phi_eta_split: polynomial mask width {} != n_basis {n_basis}",
                linear_mask.len()
            ));
        }
        Ok(PhiEtaSplit::from_curved_mask(
            linear_mask.into_iter().map(|linear| !linear).collect(),
        ))
    }

    fn second_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array4<f64>, String>> {
        Some(<Self as SaeBasisSecondJet>::second_jet(self, coords))
    }

    fn third_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array5<f64>, String>> {
        Some(<Self as SaeBasisThirdJet>::third_jet(self, coords))
    }

    fn evaluate(&self, coords: ArrayView2<'_, f64>) -> Result<(Array2<f64>, Array3<f64>), String> {
        // Single source of truth: allocate correctly-shaped buffers and fill
        // them through the in-place path (see `evaluate_into`).
        let n = coords.nrows();
        let m = self.basis_size();
        let d = self.latent_dim;
        let mut phi = Array2::<f64>::zeros((n, m));
        let mut jet = Array3::<f64>::zeros((n, m, d));
        self.evaluate_into(&mut phi, &mut jet, coords)?;
        Ok((phi, jet))
    }

    fn evaluate_into(
        &self,
        phi: &mut Array2<f64>,
        jet: &mut Array3<f64>,
        coords: ArrayView2<'_, f64>,
    ) -> Result<(), String> {
        let d = self.latent_dim;
        if coords.ncols() != d {
            return Err(format!(
                "EuclideanPatchEvaluator: expected latent_dim {}, got {}",
                self.latent_dim,
                coords.ncols()
            ));
        }
        let exponents = gam_terms::basis::monomial_exponents(self.latent_dim, self.max_degree);
        let n = coords.nrows();
        let m = exponents.len();
        if phi.dim() != (n, m) {
            return Err(format!(
                "EuclideanPatchEvaluator::evaluate_into: Φ buffer {:?} != ({n}, {m})",
                phi.dim()
            ));
        }
        if jet.dim() != (n, m, d) {
            return Err(format!(
                "EuclideanPatchEvaluator::evaluate_into: jet buffer {:?} != ({n}, {m}, {d})",
                jet.dim()
            ));
        }
        // The jet is nonzero only where a monomial's axis exponent is positive,
        // so most entries are left untouched by the loops below; clear both
        // buffers to erase any stale (reused-workspace) contents first.
        phi.fill(0.0);
        jet.fill(0.0);
        for (col, alpha) in exponents.iter().enumerate() {
            for row in 0..n {
                let mut value = 1.0_f64;
                for (axis, &exp) in alpha.iter().enumerate() {
                    if exp != 0 {
                        value *= coords[[row, axis]].powi(exp as i32);
                    }
                }
                phi[[row, col]] = value;
            }
            // Monomial first derivative, written to match
            // `gam_terms::basis::duchon_polynomial_first_derivative_nd` operation
            // for operation (same factor ordering) so the value is bit-identical.
            for axis in 0..d {
                let a_axis = alpha[axis];
                if a_axis == 0 {
                    continue;
                }
                for row in 0..n {
                    let mut value = a_axis as f64;
                    for a in 0..d {
                        let exp_a = if a == axis { a_axis - 1 } else { alpha[a] };
                        if exp_a != 0 {
                            value *= coords[[row, a]].powi(exp_a as i32);
                        }
                    }
                    jet[[row, col, axis]] = value;
                }
            }
        }
        Ok(())
    }
}

impl SaeBasisSecondJet for EuclideanPatchEvaluator {
    fn second_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array4<f64>, String> {
        if coords.ncols() != self.latent_dim {
            return Err(format!(
                "EuclideanPatchEvaluator::second_jet: expected latent_dim {}, got {}",
                self.latent_dim,
                coords.ncols()
            ));
        }
        let exponents = gam_terms::basis::monomial_exponents(self.latent_dim, self.max_degree);
        let n = coords.nrows();
        let m = exponents.len();
        let d = self.latent_dim;
        let mut hess = Array4::<f64>::zeros((n, m, d, d));
        for (col, alpha) in exponents.iter().enumerate() {
            for a in 0..d {
                if alpha[a] == 0 {
                    continue;
                }
                for c in 0..d {
                    if a != c && alpha[c] == 0 {
                        continue;
                    }
                    let lead = if a == c {
                        (alpha[a] as f64) * (alpha[a].saturating_sub(1) as f64)
                    } else {
                        (alpha[a] as f64) * (alpha[c] as f64)
                    };
                    if lead == 0.0 {
                        continue;
                    }
                    for row in 0..n {
                        let mut value = lead;
                        for axis in 0..d {
                            let mut exp = alpha[axis];
                            if axis == a {
                                exp = exp.saturating_sub(1);
                            }
                            if axis == c {
                                exp = exp.saturating_sub(1);
                            }
                            if exp != 0 {
                                value *= coords[[row, axis]].powi(exp as i32);
                            }
                        }
                        hess[[row, col, a, c]] = value;
                    }
                }
            }
        }
        Ok(hess)
    }
}

impl SaeBasisThirdJet for EuclideanPatchEvaluator {
    /// Third derivative of the monomial basis `Φ_α = Π_axis t_axis^{α_axis}`.
    ///
    /// Differentiating axis `j` a total of `k_j` times (where `k_j` is how
    /// often axis `j` appears in `{a, b, c}`) contracts that factor to
    /// `falling(α_j, k_j) · t_j^{α_j − k_j}`, with `falling(α, k) = α(α−1)…
    /// (α−k+1)` and the term vanishing whenever `α_j < k_j`.
    fn third_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array5<f64>, String> {
        if coords.ncols() != self.latent_dim {
            return Err(format!(
                "EuclideanPatchEvaluator::third_jet: expected latent_dim {}, got {}",
                self.latent_dim,
                coords.ncols()
            ));
        }
        let exponents = gam_terms::basis::monomial_exponents(self.latent_dim, self.max_degree);
        let n = coords.nrows();
        let m = exponents.len();
        let d = self.latent_dim;
        let mut t3 = Array5::<f64>::zeros((n, m, d, d, d));
        let falling = |alpha: usize, k: usize| -> f64 {
            let mut acc = 1.0_f64;
            for j in 0..k {
                acc *= (alpha as f64) - (j as f64);
            }
            acc
        };
        for (col, alpha) in exponents.iter().enumerate() {
            for a in 0..d {
                if alpha[a] == 0 {
                    continue;
                }
                for b in 0..d {
                    for c in 0..d {
                        // Per-axis differentiation order in this (a, b, c) cell.
                        let mut order = vec![0usize; d];
                        order[a] += 1;
                        order[b] += 1;
                        order[c] += 1;
                        if (0..d).any(|axis| order[axis] > alpha[axis]) {
                            continue;
                        }
                        let mut lead = 1.0_f64;
                        for axis in 0..d {
                            lead *= falling(alpha[axis], order[axis]);
                        }
                        if lead == 0.0 {
                            continue;
                        }
                        for row in 0..n {
                            let mut value = lead;
                            for axis in 0..d {
                                let exp = alpha[axis] - order[axis];
                                if exp != 0 {
                                    value *= coords[[row, axis]].powi(exp as i32);
                                }
                            }
                            t3[[row, col, a, b, c]] = value;
                        }
                    }
                }
            }
        }
        Ok(t3)
    }
}

/// Tensor-product harmonic × polynomial evaluator for the cylinder `S¹ × ℝ`
/// (`d = 2`): a periodic circle axis crossed with a flat (Duchon-polynomial)
/// line axis.
///
/// This is the missing geometry in the `d = 2` topology race (torus vs sphere
/// vs euclidean-patch vs **cylinder**). A feature whose latent structure is
/// *periodic along one axis and unbounded-linear along the other* (e.g. a
/// phase-times-magnitude direction in the residual stream) lives on `S¹ × ℝ`,
/// not on `T²` (two circles) or `S²`; before this evaluator it was forced into
/// a torus stand-in (wrapping the linear axis spuriously) or a flat patch
/// (losing the periodicity). With it, the cylinder is a first-class candidate
/// adjudicated under the same TK-normalized evidence gate.
///
/// # Basis
///
/// Axis 0 is the circle, in the `PeriodicHarmonicEvaluator` form
/// `c(t₀) = [1, sin(2π·1·t₀), cos(2π·1·t₀), …, sin(2π·H·t₀), cos(2π·H·t₀)]`
/// (`Mc = 2H + 1` columns; `t₀` is a fraction of one period, matching the
/// periodic/torus convention). Axis 1 is the flat line, in the
/// `EuclideanPatchEvaluator` form `l(t₁) = [1, t₁, t₁², …, t₁^D]`
/// (`Ml = D + 1` columns). The product basis is
/// `Φ_{c,l}(t) = c_c(t₀) · l_l(t₁)`, `M = Mc · Ml` columns, enumerated
/// lexicographically with the circle index slowest (`col = c·Ml + l`), so the
/// constant `[c=0, l=0]` is column 0.
///
/// # Jets (exact product rule across the two axes)
///
/// Because the two factors depend on disjoint coordinates, the value/derivative
/// of `Φ_{c,l}` along a multi-index over `{axis 0, axis 1}` is simply the
/// product of the circle factor differentiated as many times as axis 0 appears
/// and the line factor differentiated as many times as axis 1 appears. Each
/// per-axis derivative table (orders 0..3) is closed form (circle: the Fourier
/// `sin → ωcos → −ω²sin → −ω³cos` chain; line: the falling-factorial monomial
/// chain), so the value, first, second and third jets are all exact and the
/// second/third jets pin against a finite difference of the level below.
#[derive(Debug, Clone)]
pub struct CylinderHarmonicEvaluator {
    /// Number of circle harmonics `H ≥ 1` (axis-0 width is `2H + 1`).
    pub circle_harmonics: usize,
    /// Polynomial degree `D ≥ 0` of the flat line axis (axis-1 width is `D + 1`).
    pub line_degree: usize,
}

impl CylinderHarmonicEvaluator {
    pub fn new(circle_harmonics: usize, line_degree: usize) -> Result<Self, String> {
        if circle_harmonics == 0 {
            return Err(
                "CylinderHarmonicEvaluator requires circle_harmonics >= 1 (S¹ needs at least one \
                 harmonic pair)"
                    .to_string(),
            );
        }
        Ok(Self {
            circle_harmonics,
            line_degree,
        })
    }

    /// Circle-axis width `Mc = 2H + 1`.
    pub fn circle_basis_size(&self) -> usize {
        2 * self.circle_harmonics + 1
    }

    /// Line-axis width `Ml = D + 1`.
    pub fn line_basis_size(&self) -> usize {
        self.line_degree + 1
    }

    /// Product basis width `M = Mc · Ml`.
    pub fn basis_size(&self) -> usize {
        self.circle_basis_size() * self.line_basis_size()
    }

    /// Per-axis circle derivative tables, orders 0..=3, indexed `[order][col]`,
    /// length `Mc` each. Column 0 is the constant; columns `2h-1`/`2h` carry the
    /// `h`-th sin/cos pair at frequency `ω = 2π h`.
    fn circle_tables(&self, t: f64) -> [Vec<f64>; 4] {
        let mc = self.circle_basis_size();
        let two_pi = 2.0 * std::f64::consts::PI;
        let mut table = [
            vec![0.0_f64; mc],
            vec![0.0_f64; mc],
            vec![0.0_f64; mc],
            vec![0.0_f64; mc],
        ];
        // Constant column: value 1, all derivatives 0.
        table[0][0] = 1.0;
        for h in 1..=self.circle_harmonics {
            let omega = two_pi * (h as f64);
            let w2 = omega * omega;
            let w3 = w2 * omega;
            let angle = omega * t;
            let s = angle.sin();
            let c = angle.cos();
            let s_idx = 2 * h - 1;
            let c_idx = 2 * h;
            // sin chain: sin → ω cos → −ω² sin → −ω³ cos.
            table[0][s_idx] = s;
            table[1][s_idx] = omega * c;
            table[2][s_idx] = -w2 * s;
            table[3][s_idx] = -w3 * c;
            // cos chain: cos → −ω sin → −ω² cos → ω³ sin.
            table[0][c_idx] = c;
            table[1][c_idx] = -omega * s;
            table[2][c_idx] = -w2 * c;
            table[3][c_idx] = w3 * s;
        }
        table
    }

    /// Per-axis line (monomial) derivative tables, orders 0..=3, indexed
    /// `[order][col]`, length `Ml` each. Column `j` is `t^j`; its `k`-th
    /// derivative is `falling(j, k) · t^{j-k}` (zero once `k > j`).
    fn line_tables(&self, t: f64) -> [Vec<f64>; 4] {
        let ml = self.line_basis_size();
        let mut table = [
            vec![0.0_f64; ml],
            vec![0.0_f64; ml],
            vec![0.0_f64; ml],
            vec![0.0_f64; ml],
        ];
        for j in 0..ml {
            for k in 0..4 {
                if k > j {
                    // falling(j, k) = 0: the monomial is exhausted.
                    table[k][j] = 0.0;
                    continue;
                }
                let mut coeff = 1.0_f64;
                for q in 0..k {
                    coeff *= (j - q) as f64;
                }
                let residual = j - k;
                let pow = if residual == 0 {
                    1.0
                } else {
                    t.powi(residual as i32)
                };
                table[k][j] = coeff * pow;
            }
        }
        table
    }

    /// Analytic seed roughness Gram `S = ∫ (LΦ)ᵀ (LΦ)` for the cylinder, built
    /// as the tensor sum of a per-axis curvature energy: it penalizes the
    /// second derivative along the circle (the bending energy of the periodic
    /// factor) plus the second derivative along the line (the thin-plate energy
    /// of the flat factor). Because the basis is a clean tensor product and the
    /// two coordinate measures are independent on `[0,1) × ℝ`, the cross terms
    /// factor through the per-axis Grams.
    ///
    /// Concretely, with `Sc` the circle second-derivative Gram and `Gc` the
    /// circle value Gram (both `Mc × Mc`), `Sl`/`Gl` their line counterparts
    /// (`Ml × Ml`), the roughness operator
    /// `‖∂²_{t₀}Φ‖² + ‖∂²_{t₁}Φ‖²` integrates to
    /// `S = Sc ⊗ Gl + Gc ⊗ Sl` in the same lexicographic column order as the
    /// design. This is gauge-invariant: it depends only on the basis functions,
    /// not on any chart-specific normalization, and the constant column (zero in
    /// both `Sc` and `Sl`) sits in the null space exactly as the smooth-penalty
    /// nullity recovery expects.
    ///
    /// The circle blocks use the closed-form Fourier integrals on `[0,1)`:
    /// `∫₀¹ 1 dt = 1`, `∫₀¹ sin²(2πht) = ∫₀¹ cos²(2πht) = ½`, all distinct-mode
    /// and sin·cos cross integrals vanish, so `Gc` is diagonal
    /// `diag(1, ½, ½, …)` and `Sc = diag(0, (2πh)⁴·½, …)` (the second derivative
    /// of a mode scales its value by `(2πh)²`, squared and integrated → `(2πh)⁴`
    /// times the value integral). The line blocks use the monomial moments on a
    /// canonical unit interval `[0,1)` so the energy is finite and scale-fixed:
    /// `Gl[i,j] = ∫₀¹ tⁱ⁺ʲ dt = 1/(i+j+1)` and
    /// `Sl[i,j] = ∫₀¹ (i(i-1)t^{i-2})(j(j-1)t^{j-2}) dt`.
    pub fn roughness_gram(&self) -> Array2<f64> {
        let mc = self.circle_basis_size();
        let ml = self.line_basis_size();
        let two_pi = 2.0 * std::f64::consts::PI;

        // Circle value Gram Gc and second-derivative Gram Sc (both diagonal).
        let mut gc = Array2::<f64>::zeros((mc, mc));
        let mut sc = Array2::<f64>::zeros((mc, mc));
        gc[[0, 0]] = 1.0; // ∫₀¹ 1 dt
        for h in 1..=self.circle_harmonics {
            let omega = two_pi * (h as f64);
            let w4 = omega.powi(4);
            let s_idx = 2 * h - 1;
            let c_idx = 2 * h;
            // ∫₀¹ sin² = ∫₀¹ cos² = ½.
            gc[[s_idx, s_idx]] = 0.5;
            gc[[c_idx, c_idx]] = 0.5;
            // Second derivative scales the mode by (2πh)²; squared·integrated.
            sc[[s_idx, s_idx]] = w4 * 0.5;
            sc[[c_idx, c_idx]] = w4 * 0.5;
        }

        // Line value Gram Gl and second-derivative Gram Sl on the canonical
        // reference interval [0,1).
        //
        // NOTE: despite the cylinder being declared `S¹ × ℝ` (unbounded line
        // axis) in the latent-manifold spec, this line-factor roughness is a
        // *compact* [0,1) reference-domain measure, NOT an intrinsic roughness
        // integrated over all of ℝ (a uniform ∫_ℝ tⁱ⁺ʲ dt would diverge for
        // polynomials). It is a canonical reference penalty: the line
        // coordinate's scale and origin therefore matter through this reference
        // interval — the same monomial coefficients carry a different penalty
        // if the line axis is rescaled or shifted relative to [0,1).
        let mut gl = Array2::<f64>::zeros((ml, ml));
        let mut sl = Array2::<f64>::zeros((ml, ml));
        for i in 0..ml {
            for j in 0..ml {
                // Gl[i,j] = ∫₀¹ t^{i+j} dt = 1/(i+j+1).
                gl[[i, j]] = 1.0 / ((i + j + 1) as f64);
                // Sl[i,j] = ∫₀¹ (i(i-1) t^{i-2})(j(j-1) t^{j-2}) dt.
                if i >= 2 && j >= 2 {
                    let ci = (i * (i - 1)) as f64;
                    let cj = (j * (j - 1)) as f64;
                    let exp = (i - 2) + (j - 2);
                    sl[[i, j]] = ci * cj / ((exp + 1) as f64);
                }
            }
        }

        // S = Sc ⊗ Gl + Gc ⊗ Sl in lexicographic (circle-slow, line-fast) order.
        let m = mc * ml;
        let mut s = Array2::<f64>::zeros((m, m));
        for ca in 0..mc {
            for la in 0..ml {
                let row = ca * ml + la;
                for cb in 0..mc {
                    for lb in 0..ml {
                        let col = cb * ml + lb;
                        s[[row, col]] = sc[[ca, cb]] * gl[[la, lb]] + gc[[ca, cb]] * sl[[la, lb]];
                    }
                }
            }
        }
        s
    }

    fn check_coords(&self, coords: ArrayView2<'_, f64>, what: &str) -> Result<(), String> {
        if coords.ncols() != 2 {
            return Err(format!(
                "CylinderHarmonicEvaluator::{what}: expected latent_dim == 2 (S¹ × ℝ), got {}",
                coords.ncols()
            ));
        }
        Ok(())
    }
}

impl SaeBasisEvaluator for CylinderHarmonicEvaluator {
    fn phi_eta_split(&self, n_basis: usize) -> Result<PhiEtaSplit, String> {
        let expected = self.basis_size();
        if n_basis != expected {
            return Err(format!(
                "CylinderHarmonicEvaluator::phi_eta_split: n_basis {n_basis} != evaluator width {expected}"
            ));
        }
        let ml = self.line_basis_size();
        // A product column `[c, l]` is curved (η-dialed) iff either factor
        // carries dialed curvature: the circle factor on any harmonic above the
        // first (`c > 2`, the 2nd-and-higher sin/cos), the line factor on any
        // monomial of degree ≥ 2 (`l > 1`). The first-harmonic circle columns
        // crossed with the affine line columns are the base (η-invariant) block
        // and span the base-topology relaxation — the base circle columns still
        // embed the unit circle, so this block is not linear — matching the
        // `eta`-homotopy split of the two parent evaluators.
        let mut curved = vec![false; expected];
        for c in 0..self.circle_basis_size() {
            for l in 0..ml {
                let circle_curved = c > 2;
                let line_curved = l > 1;
                curved[c * ml + l] = circle_curved || line_curved;
            }
        }
        Ok(PhiEtaSplit::from_curved_mask(curved))
    }

    /// A cylinder `S¹ × ℝ` is the tensor product of the circle factor
    /// (`circle_basis_size()`, the slow axis 0) and the line factor
    /// (`line_basis_size()`, the fast axis 1), both constant-leading.
    fn factor_basis_sizes(&self) -> Option<(usize, usize)> {
        Some((self.circle_basis_size(), self.line_basis_size()))
    }

    fn second_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array4<f64>, String>> {
        Some(<Self as SaeBasisSecondJet>::second_jet(self, coords))
    }

    fn third_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array5<f64>, String>> {
        Some(<Self as SaeBasisThirdJet>::third_jet(self, coords))
    }

    fn evaluate(&self, coords: ArrayView2<'_, f64>) -> Result<(Array2<f64>, Array3<f64>), String> {
        self.check_coords(coords, "evaluate")?;
        let n = coords.nrows();
        let mc = self.circle_basis_size();
        let ml = self.line_basis_size();
        let m = mc * ml;
        let mut phi = Array2::<f64>::zeros((n, m));
        let mut jet = Array3::<f64>::zeros((n, m, 2));
        for row in 0..n {
            let t0 = coords[[row, 0]];
            let t1 = coords[[row, 1]];
            let circ = self.circle_tables(t0);
            let line = self.line_tables(t1);
            for c in 0..mc {
                for l in 0..ml {
                    let col = c * ml + l;
                    // Value: c·l. ∂/∂t₀ = c'·l. ∂/∂t₁ = c·l'.
                    phi[[row, col]] = circ[0][c] * line[0][l];
                    jet[[row, col, 0]] = circ[1][c] * line[0][l];
                    jet[[row, col, 1]] = circ[0][c] * line[1][l];
                }
            }
        }
        Ok((phi, jet))
    }
}

impl SaeBasisSecondJet for CylinderHarmonicEvaluator {
    /// Hessian of the cylinder product basis. With `Φ_{c,l} = c(t₀)·l(t₁)` and
    /// the two factors on disjoint coordinates:
    /// `∂²/∂t₀² = c''·l`, `∂²/∂t₁² = c·l''`, `∂²/∂t₀∂t₁ = c'·l'`.
    fn second_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array4<f64>, String> {
        self.check_coords(coords, "second_jet")?;
        let n = coords.nrows();
        let mc = self.circle_basis_size();
        let ml = self.line_basis_size();
        let m = mc * ml;
        let mut h = Array4::<f64>::zeros((n, m, 2, 2));
        for row in 0..n {
            let t0 = coords[[row, 0]];
            let t1 = coords[[row, 1]];
            let circ = self.circle_tables(t0);
            let line = self.line_tables(t1);
            for c in 0..mc {
                for l in 0..ml {
                    let col = c * ml + l;
                    h[[row, col, 0, 0]] = circ[2][c] * line[0][l];
                    h[[row, col, 1, 1]] = circ[0][c] * line[2][l];
                    let mixed = circ[1][c] * line[1][l];
                    h[[row, col, 0, 1]] = mixed;
                    h[[row, col, 1, 0]] = mixed;
                }
            }
        }
        Ok(h)
    }
}

impl SaeBasisThirdJet for CylinderHarmonicEvaluator {
    /// Third derivative of the cylinder product basis. The number of axis-0
    /// derivative operators `k₀` (and axis-1 `k₁ = 3 − k₀`) in a cell `(a,b,e)`
    /// routes that many derivatives to the circle factor and the rest to the
    /// line factor: `∂³Φ = c^{(k₀)}(t₀) · l^{(k₁)}(t₁)`. This is the order-3
    /// sibling of [`SaeBasisSecondJet::second_jet`].
    fn third_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array5<f64>, String> {
        self.check_coords(coords, "third_jet")?;
        let n = coords.nrows();
        let mc = self.circle_basis_size();
        let ml = self.line_basis_size();
        let m = mc * ml;
        let mut t3 = Array5::<f64>::zeros((n, m, 2, 2, 2));
        for row in 0..n {
            let t0 = coords[[row, 0]];
            let t1 = coords[[row, 1]];
            let circ = self.circle_tables(t0);
            let line = self.line_tables(t1);
            for c in 0..mc {
                for l in 0..ml {
                    let col = c * ml + l;
                    for a in 0..2 {
                        for b in 0..2 {
                            for e in 0..2 {
                                // k0 = number of axis-0 operators among (a,b,e).
                                let k0 = (a == 0) as usize + (b == 0) as usize + (e == 0) as usize;
                                let k1 = 3 - k0;
                                t3[[row, col, a, b, e]] = circ[k0][c] * line[k1][l];
                            }
                        }
                    }
                }
            }
        }
        Ok(t3)
    }
}

/// Möbius-band harmonic basis on the DOUBLE-COVER chart (#2240).
///
/// The band is charted by `(s, w)` with `s ∈ [0, 2)` an angle on the
/// double-cover circle (period 2) and `w ∈ [-1, 1]` the width coordinate. The
/// deck transformation of the double cover is `σ(s, w) = (s + 1, -w)`; the
/// Möbius band is the quotient, and functions on the band are exactly the
/// σ-invariant functions on the cylinder cover. With the tensor harmonics
/// `T_k(s)·w^m` (`T_k ∈ {cos(πks), sin(πks)}`, period-2 modes), invariance is
/// `(-1)^{k+m} = 1`, so the basis keeps exactly the columns with `k + m`
/// EVEN: `(k even, m even) ∪ (k odd, m odd)`.
///
/// This gives the band its defining behavior with NO new manifold machinery:
/// the optimizer retracts on an ordinary smooth cylinder (`Circle{period 2} ×
/// Interval[-1, 1]`), every basis column is C^∞ there, and a point and its
/// deck-twin `(s+1, -w)` produce IDENTICAL basis rows — the half-twist lives
/// in the parity culling, not in a seam. Odd-`m` (width-odd) structure is
/// forced to carry a half-period angular factor, which is precisely the
/// non-orientability a torus or flat-patch chart cannot express.
#[derive(Debug, Clone)]
pub struct MobiusHarmonicEvaluator {
    /// Circle harmonics `H ≥ 1` on the double-cover angle (mode `k ≤ H`).
    pub circle_harmonics: usize,
    /// Width monomial degree `D ≥ 1` (`w^m`, `m ≤ D`).
    pub width_degree: usize,
    /// Admitted `(circle_col, width_power)` pairs (deck-invariant columns),
    /// circle-column-slow / width-power-fast, fixed at construction.
    columns: Vec<(usize, usize)>,
}

impl MobiusHarmonicEvaluator {
    pub fn new(circle_harmonics: usize, width_degree: usize) -> Result<Self, String> {
        if circle_harmonics == 0 {
            return Err(
                "MobiusHarmonicEvaluator requires circle_harmonics >= 1 (the band core needs \
                 at least the half-period harmonic pair)"
                    .to_string(),
            );
        }
        if width_degree == 0 {
            return Err(
                "MobiusHarmonicEvaluator requires width_degree >= 1: with no width-odd \
                 columns the deck-invariant basis degenerates to a plain circle"
                    .to_string(),
            );
        }
        let mc = 2 * circle_harmonics + 1;
        let mut columns = Vec::new();
        for c in 0..mc {
            let k = Self::circle_mode(c);
            for m in 0..=width_degree {
                if (k + m) % 2 == 0 {
                    columns.push((c, m));
                }
            }
        }
        Ok(Self {
            circle_harmonics,
            width_degree,
            columns,
        })
    }

    /// Angular mode `k` of circle column `c` (col 0 = constant, cols
    /// `2h-1`/`2h` = the sin/cos pair of mode `h`).
    fn circle_mode(c: usize) -> usize {
        c.div_ceil(2)
    }

    pub fn basis_size(&self) -> usize {
        self.columns.len()
    }

    /// Circle derivative tables on the DOUBLE-COVER angle, orders 0..=3:
    /// mode `k` has frequency `ω = πk` (period 2), so odd modes are the
    /// half-period harmonics the quotient demands.
    fn circle_tables(&self, s: f64) -> [Vec<f64>; 4] {
        let mc = 2 * self.circle_harmonics + 1;
        let pi = std::f64::consts::PI;
        let mut table = [
            vec![0.0_f64; mc],
            vec![0.0_f64; mc],
            vec![0.0_f64; mc],
            vec![0.0_f64; mc],
        ];
        table[0][0] = 1.0;
        for h in 1..=self.circle_harmonics {
            let omega = pi * (h as f64);
            let w2 = omega * omega;
            let w3 = w2 * omega;
            let angle = omega * s;
            let sv = angle.sin();
            let cv = angle.cos();
            let s_idx = 2 * h - 1;
            let c_idx = 2 * h;
            table[0][s_idx] = sv;
            table[1][s_idx] = omega * cv;
            table[2][s_idx] = -w2 * sv;
            table[3][s_idx] = -w3 * cv;
            table[0][c_idx] = cv;
            table[1][c_idx] = -omega * sv;
            table[2][c_idx] = -w2 * cv;
            table[3][c_idx] = w3 * sv;
        }
        table
    }

    /// Width (monomial) derivative tables on `[-1, 1]`, orders 0..=3.
    fn width_tables(&self, w: f64) -> [Vec<f64>; 4] {
        let mw = self.width_degree + 1;
        let mut table = [
            vec![0.0_f64; mw],
            vec![0.0_f64; mw],
            vec![0.0_f64; mw],
            vec![0.0_f64; mw],
        ];
        for j in 0..mw {
            for k in 0..4 {
                if k > j {
                    table[k][j] = 0.0;
                    continue;
                }
                let mut coeff = 1.0_f64;
                for q in 0..k {
                    coeff *= (j - q) as f64;
                }
                let residual = j - k;
                let pow = if residual == 0 {
                    1.0
                } else {
                    w.powi(residual as i32)
                };
                table[k][j] = coeff * pow;
            }
        }
        table
    }

    /// Analytic seed roughness Gram `S = Sc ⊗ Gw + Gc ⊗ Sw` restricted to the
    /// admitted deck-invariant columns. On the double-cover measure
    /// `[0, 2) × [-1, 1]` the circle Grams are diagonal (`∫₀² 1 = 2`,
    /// `∫₀² sin²(πks) = ∫₀² cos²(πks) = 1`, all cross terms vanish over the
    /// full period; the second derivative scales a mode by `(πk)²`), and the
    /// width Grams are the even-moment tables `∫₋₁¹ w^{i+j} dw`
    /// (`= 2/(i+j+1)` for `i+j` even, `0` odd). The constant column sits in
    /// the null space exactly as the smooth-penalty nullity recovery expects.
    pub fn roughness_gram(&self) -> Array2<f64> {
        let pi = std::f64::consts::PI;
        let m = self.columns.len();
        let moment = |exp: usize| -> f64 {
            if exp % 2 == 0 {
                2.0 / ((exp + 1) as f64)
            } else {
                0.0
            }
        };
        let mut s = Array2::<f64>::zeros((m, m));
        for (row, &(c_a, m_a)) in self.columns.iter().enumerate() {
            for (col, &(c_b, m_b)) in self.columns.iter().enumerate() {
                if c_a != c_b {
                    // Distinct circle columns are orthogonal in BOTH circle
                    // Grams over the full double-cover period.
                    continue;
                }
                let k = Self::circle_mode(c_a) as f64;
                let gc = if c_a == 0 { 2.0 } else { 1.0 };
                let sc = if c_a == 0 {
                    0.0
                } else {
                    (pi * k).powi(4) * 1.0
                };
                let gw = moment(m_a + m_b);
                let sw = if m_a >= 2 && m_b >= 2 {
                    ((m_a * (m_a - 1)) as f64) * ((m_b * (m_b - 1)) as f64) * moment(m_a + m_b - 4)
                } else {
                    0.0
                };
                s[[row, col]] = sc * gw + gc * sw;
            }
        }
        s
    }

    fn check_coords(&self, coords: ArrayView2<'_, f64>, what: &str) -> Result<(), String> {
        if coords.ncols() != 2 {
            return Err(format!(
                "MobiusHarmonicEvaluator::{what}: expected latent_dim == 2 (double-cover \
                 angle × width), got {}",
                coords.ncols()
            ));
        }
        Ok(())
    }
}

impl SaeBasisEvaluator for MobiusHarmonicEvaluator {
    fn phi_eta_split(&self, n_basis: usize) -> Result<PhiEtaSplit, String> {
        let expected = self.basis_size();
        if n_basis != expected {
            return Err(format!(
                "MobiusHarmonicEvaluator::phi_eta_split: n_basis {n_basis} != evaluator width {expected}"
            ));
        }
        // Base (η-invariant) block: the constant, the half-period first
        // harmonic pair crossed with the affine width — the band's core
        // embedding. Higher angular modes (k ≥ 2) or width curvature (m ≥ 2)
        // are the η-dialed refinement, mirroring the cylinder's split.
        let curved = self
            .columns
            .iter()
            .map(|&(c, m)| Self::circle_mode(c) >= 2 || m >= 2)
            .collect::<Vec<_>>();
        Ok(PhiEtaSplit::from_curved_mask(curved))
    }

    /// The parity culling breaks the clean tensor-product factorization, so
    /// the Möbius basis does not expose per-axis factor sizes.
    fn factor_basis_sizes(&self) -> Option<(usize, usize)> {
        None
    }

    fn second_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array4<f64>, String>> {
        Some(<Self as SaeBasisSecondJet>::second_jet(self, coords))
    }

    fn third_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array5<f64>, String>> {
        Some(<Self as SaeBasisThirdJet>::third_jet(self, coords))
    }

    fn evaluate(&self, coords: ArrayView2<'_, f64>) -> Result<(Array2<f64>, Array3<f64>), String> {
        self.check_coords(coords, "evaluate")?;
        let n = coords.nrows();
        let m = self.basis_size();
        let mut phi = Array2::<f64>::zeros((n, m));
        let mut jet = Array3::<f64>::zeros((n, m, 2));
        for row in 0..n {
            let circ = self.circle_tables(coords[[row, 0]]);
            let width = self.width_tables(coords[[row, 1]]);
            for (col, &(c, wm)) in self.columns.iter().enumerate() {
                phi[[row, col]] = circ[0][c] * width[0][wm];
                jet[[row, col, 0]] = circ[1][c] * width[0][wm];
                jet[[row, col, 1]] = circ[0][c] * width[1][wm];
            }
        }
        Ok((phi, jet))
    }
}

impl SaeBasisSecondJet for MobiusHarmonicEvaluator {
    fn second_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array4<f64>, String> {
        self.check_coords(coords, "second_jet")?;
        let n = coords.nrows();
        let m = self.basis_size();
        let mut h = Array4::<f64>::zeros((n, m, 2, 2));
        for row in 0..n {
            let circ = self.circle_tables(coords[[row, 0]]);
            let width = self.width_tables(coords[[row, 1]]);
            for (col, &(c, wm)) in self.columns.iter().enumerate() {
                h[[row, col, 0, 0]] = circ[2][c] * width[0][wm];
                h[[row, col, 1, 1]] = circ[0][c] * width[2][wm];
                let mixed = circ[1][c] * width[1][wm];
                h[[row, col, 0, 1]] = mixed;
                h[[row, col, 1, 0]] = mixed;
            }
        }
        Ok(h)
    }
}

impl SaeBasisThirdJet for MobiusHarmonicEvaluator {
    fn third_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array5<f64>, String> {
        self.check_coords(coords, "third_jet")?;
        let n = coords.nrows();
        let m = self.basis_size();
        let mut t3 = Array5::<f64>::zeros((n, m, 2, 2, 2));
        for row in 0..n {
            let circ = self.circle_tables(coords[[row, 0]]);
            let width = self.width_tables(coords[[row, 1]]);
            for (col, &(c, wm)) in self.columns.iter().enumerate() {
                for a in 0..2 {
                    for b in 0..2 {
                        for e in 0..2 {
                            let k0 = (a == 0) as usize + (b == 0) as usize + (e == 0) as usize;
                            let k1 = 3 - k0;
                            t3[[row, col, a, b, e]] = circ[k0][c] * width[k1][wm];
                        }
                    }
                }
            }
        }
        Ok(t3)
    }
}

/// Rank-revealing subspace reparametrization of an inner basis evaluator.
///
/// Issue #1117: a decoder basis (e.g. [`PeriodicHarmonicEvaluator`]) emits a
/// *fixed* number of columns `M` independent of the data, so on a
/// near-degenerate checkpoint the higher columns are unexcited and the
/// decoder design `Φ` is rank-deficient by construction (the OLMo
/// `stage1-step0` PCA-32 circle: data Gram rank `3/5`, the 2nd-harmonic pair
/// dead). A rank-deficient design leaves the inner solve conditioned only by
/// ridges/deflation and flattens the outer REML surface, stalling BFGS.
///
/// This wrapper makes the design **full-rank by construction**: the
/// data-supported subspace of the inner basis is discovered ONCE at fit entry
/// (the eigenvectors of the weighted data Gram `G = Φᵀ W Φ` whose eigenvalue
/// clears the relative spectral cutoff) and frozen into an orthonormal column
/// map `Q ∈ ℝ^{M × r}` (`r = rank(G) ≤ M`). The wrapped evaluator then emits
/// the reduced design `Φ̃ = Φ Q` (and its jets `∂Φ̃ = (∂Φ) Q`, …) on every
/// refresh, so the reduction *survives* re-evaluation — unlike a step-time
/// projector, which the evaluator's next `evaluate` overwrites.
///
/// Because `Q` is a fixed linear remix of the inner columns, the reduced basis
/// is exactly as smooth as the inner one and every derivative composes by the
/// same right-multiply: `∂^g Φ̃ = (∂^g Φ) Q`, contracting only the basis
/// (column) axis. The retained columns span exactly the data-identified part of
/// the inner basis; the smooth/REML penalty then shrinks within that span and
/// is never asked to identify a direction the data cannot see.
///
/// When the inner Gram is full rank (`r == M`, the `base`/`step_2300` case),
/// the fit-entry installer skips the wrap entirely and the inner evaluator is
/// used unchanged, so the well-conditioned path is byte-identical.
#[derive(Debug, Clone)]
pub struct SubspaceReducedEvaluator {
    inner: Arc<dyn SaeBasisSecondJet>,
    /// `(M × r)` orthonormal column map onto the data-supported subspace.
    q: Array2<f64>,
}

impl SubspaceReducedEvaluator {
    /// Wrap `inner` with the column map `q` (`M_inner × r`, `r ≤ M_inner`). The
    /// retained width is `q.ncols()`. The columns of `q` are expected to be
    /// orthonormal (the eigenvectors of a symmetric data Gram); orthonormality
    /// is not re-checked here — it is the caller's contract at fit entry.
    pub fn new(inner: Arc<dyn SaeBasisSecondJet>, q: Array2<f64>) -> Result<Self, String> {
        if q.nrows() == 0 || q.ncols() == 0 {
            return Err(format!(
                "SubspaceReducedEvaluator: column map must be non-empty; got {:?}",
                q.dim()
            ));
        }
        if q.ncols() > q.nrows() {
            return Err(format!(
                "SubspaceReducedEvaluator: retained rank {} exceeds inner basis width {}",
                q.ncols(),
                q.nrows()
            ));
        }
        Ok(Self { inner, q })
    }

    /// Inner basis width `M` (the wrapped evaluator's column count).
    pub fn inner_width(&self) -> usize {
        self.q.nrows()
    }

    /// Retained (data-supported) width `r`.
    pub fn reduced_width(&self) -> usize {
        self.q.ncols()
    }

    fn check_inner_width(&self, got: usize, what: &str) -> Result<(), String> {
        if got != self.q.nrows() {
            return Err(format!(
                "SubspaceReducedEvaluator::{what}: inner evaluator returned width {got}, \
                 column map expects {}",
                self.q.nrows()
            ));
        }
        Ok(())
    }
}

/// Right-multiply the basis (column) axis of a per-row value matrix
/// `phi` (`n × M`) by `q` (`M × r`), returning `(n × r)`.
fn remix_cols_2(phi: &Array2<f64>, q: &Array2<f64>) -> Array2<f64> {
    phi.dot(q)
}

/// Right-multiply the basis axis of a jet `jet[n, M, ..]` by `q` (`M × r`),
/// returning the same trailing shape with the basis axis reduced to `r`. The
/// trailing derivative axes are flattened, the `(M)`→`(r)` remix applied as one
/// matmul, then reshaped back; this is the exact `∂^g Φ̃ = (∂^g Φ) Q` contract.
fn remix_cols_along_basis(
    jet: ndarray::ArrayViewD<'_, f64>,
    q: &Array2<f64>,
) -> Result<ndarray::ArrayD<f64>, String> {
    let shape = jet.shape().to_vec();
    if shape.len() < 2 {
        return Err(format!(
            "SubspaceReducedEvaluator: jet must have at least (n, M) axes; got {shape:?}"
        ));
    }
    let n = shape[0];
    let m = shape[1];
    if m != q.nrows() {
        return Err(format!(
            "SubspaceReducedEvaluator: jet basis axis {m} != column-map rows {}",
            q.nrows()
        ));
    }
    let r = q.ncols();
    let trailing: usize = shape[2..].iter().product::<usize>().max(1);
    let mut out_shape = shape.clone();
    out_shape[1] = r;
    // Flatten the trailing derivative axes so the remix is a single
    // `(M)→(r)` contraction over the basis axis for every (row, trailing) fiber.
    // `to_owned()` produces a standard (row-major contiguous) layout, so the
    // flatten and the final reshape back to `out_shape` are exact.
    let jet_std = jet.to_owned();
    let jet_flat = jet_std
        .to_shape((n, m, trailing))
        .map_err(|err| format!("SubspaceReducedEvaluator: jet reshape failed: {err}"))?;
    let mut out_flat = Array3::<f64>::zeros((n, r, trailing));
    for row in 0..n {
        for t in 0..trailing {
            for rc in 0..r {
                let mut acc = 0.0_f64;
                for mc in 0..m {
                    acc += jet_flat[[row, mc, t]] * q[[mc, rc]];
                }
                out_flat[[row, rc, t]] = acc;
            }
        }
    }
    let out = out_flat
        .into_shape_with_order(ndarray::IxDyn(&out_shape))
        .map_err(|err| format!("SubspaceReducedEvaluator: out reshape failed: {err}"))?;
    Ok(out)
}

impl SaeBasisEvaluator for SubspaceReducedEvaluator {
    fn phi_eta_split(&self, n_basis: usize) -> Result<PhiEtaSplit, String> {
        if n_basis != self.q.ncols() {
            return Err(format!(
                "SubspaceReducedEvaluator::phi_eta_split: n_basis {n_basis} != reduced width {}",
                self.q.ncols()
            ));
        }
        // A reduced column is "curved" iff its data-supported direction draws on
        // any curved inner column. `Q[:, rc]` mixes inner columns; the reduced
        // column carries curvature when `Q[curved_inner, rc]` is non-zero.
        let inner_split = self.inner.phi_eta_split(self.q.nrows())?;
        let mut inner_curved = vec![false; self.q.nrows()];
        for &col in &inner_split.curved_cols {
            if col < inner_curved.len() {
                inner_curved[col] = true;
            }
        }
        let mut curved = vec![false; self.q.ncols()];
        for rc in 0..self.q.ncols() {
            for mc in 0..self.q.nrows() {
                if inner_curved[mc] && self.q[[mc, rc]] != 0.0 {
                    curved[rc] = true;
                    break;
                }
            }
        }
        Ok(PhiEtaSplit::from_curved_mask(curved))
    }

    fn second_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array4<f64>, String>> {
        Some(<Self as SaeBasisSecondJet>::second_jet(self, coords))
    }

    fn third_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array5<f64>, String>> {
        match self.inner.third_jet_dyn(coords) {
            Some(Ok(t3)) => {
                if let Err(err) = self.check_inner_width(t3.shape()[1], "third_jet_dyn") {
                    return Some(Err(err));
                }
                Some(
                    remix_cols_along_basis(t3.view().into_dyn(), &self.q).and_then(|out| {
                        out.into_dimensionality::<ndarray::Ix5>().map_err(|err| {
                            format!("SubspaceReducedEvaluator: third jet dim: {err}")
                        })
                    }),
                )
            }
            Some(Err(err)) => Some(Err(err)),
            None => None,
        }
    }

    fn evaluate(&self, coords: ArrayView2<'_, f64>) -> Result<(Array2<f64>, Array3<f64>), String> {
        let (phi, jet) = self.inner.evaluate(coords)?;
        self.check_inner_width(phi.ncols(), "evaluate")?;
        let phi_red = remix_cols_2(&phi, &self.q);
        let jet_red = remix_cols_along_basis(jet.view().into_dyn(), &self.q)?
            .into_dimensionality::<ndarray::Ix3>()
            .map_err(|err| format!("SubspaceReducedEvaluator: jet dim: {err}"))?;
        Ok((phi_red, jet_red))
    }
}

impl SaeBasisSecondJet for SubspaceReducedEvaluator {
    fn second_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array4<f64>, String> {
        let h = self.inner.second_jet(coords)?;
        self.check_inner_width(h.shape()[1], "second_jet")?;
        remix_cols_along_basis(h.view().into_dyn(), &self.q)?
            .into_dimensionality::<ndarray::Ix4>()
            .map_err(|err| format!("SubspaceReducedEvaluator: second jet dim: {err}"))
    }
}

/// F2 finite-set (discrete anchor) basis: the INDICATOR / one-hot design over a
/// fixed set of `anchors` anchors. The latent coordinate `t` (a single flat axis)
/// is read as a categorical assignment — `t` is snapped to its nearest anchor
/// index `round(t)` clamped to `[0, anchors)` — and the design row is the one-hot
/// vector selecting that anchor. This is the honest model for cluster-like
/// structure (weekdays as a finite point set, not an occupied circle): unlike
/// every continuous evaluator here, the decoded map is piecewise CONSTANT in `t`,
/// so its first, second, and third jets are identically zero (the derivative of a
/// step is zero a.e.; the anchor assignment moves by re-labelling, not by a
/// tangent step). The design width equals `anchors`; the rank charge the race
/// prices is `anchors − 1` ([`crate::manifold::finite_set_rank_charge`]), one
/// anchor being the reference contrast.
#[derive(Debug, Clone)]
pub struct AnchorIndicatorEvaluator {
    pub anchors: usize,
}

impl AnchorIndicatorEvaluator {
    pub fn new(anchors: usize) -> Result<Self, String> {
        if anchors < 2 {
            return Err(format!(
                "AnchorIndicatorEvaluator requires anchors >= 2 (a finite set of at \
                 least two points); got {anchors}"
            ));
        }
        Ok(Self { anchors })
    }
}

impl SaeBasisEvaluator for AnchorIndicatorEvaluator {
    fn second_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array4<f64>, String>> {
        Some(<Self as SaeBasisSecondJet>::second_jet(self, coords))
    }

    fn third_jet_dyn(&self, coords: ArrayView2<'_, f64>) -> Option<Result<Array5<f64>, String>> {
        // The indicator design is piecewise constant, so every jet order is zero.
        let n = coords.nrows();
        Some(Ok(Array5::<f64>::zeros((n, self.anchors, 1, 1, 1))))
    }

    fn evaluate(&self, coords: ArrayView2<'_, f64>) -> Result<(Array2<f64>, Array3<f64>), String> {
        let n = coords.nrows();
        let d = coords.ncols();
        if d != 1 {
            return Err(format!(
                "AnchorIndicatorEvaluator: expected latent_dim == 1 (a single \
                 categorical axis), got {d}"
            ));
        }
        let m = self.anchors;
        let mut phi = Array2::<f64>::zeros((n, m));
        // The jet is identically zero: a one-hot indicator is constant between
        // anchors, so ∂Φ/∂t = 0 a.e.
        let jet = Array3::<f64>::zeros((n, m, 1));
        for row in 0..n {
            let t = coords[[row, 0]];
            if !t.is_finite() {
                return Err("AnchorIndicatorEvaluator: non-finite coordinate".to_string());
            }
            // Snap to the nearest anchor index, clamped into range.
            let idx = t.round().clamp(0.0, (m - 1) as f64) as usize;
            phi[[row, idx]] = 1.0;
        }
        Ok((phi, jet))
    }
}

impl SaeBasisSecondJet for AnchorIndicatorEvaluator {
    fn second_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array4<f64>, String> {
        let n = coords.nrows();
        Ok(Array4::<f64>::zeros((n, self.anchors, 1, 1)))
    }
}

impl SaeBasisThirdJet for AnchorIndicatorEvaluator {
    fn third_jet(&self, coords: ArrayView2<'_, f64>) -> Result<Array5<f64>, String> {
        let n = coords.nrows();
        Ok(Array5::<f64>::zeros((n, self.anchors, 1, 1, 1)))
    }
}

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

    /// Deterministic splitmix64 stream in `[0, 1)` for reproducible fixtures.
    fn uniform_stream(seed: u64) -> impl FnMut() -> f64 {
        let mut state = seed;
        move || {
            state = state.wrapping_add(0x9E3779B97F4A7C15);
            let mut z = state;
            z = (z ^ (z >> 30)).wrapping_mul(0xBF58476D1CE4E5B9);
            z = (z ^ (z >> 27)).wrapping_mul(0x94D049BB133111EB);
            z ^= z >> 31;
            (z >> 11) as f64 / (1u64 << 53) as f64
        }
    }

    /// Least-squares decoder on `train_rows`, held-out R² on `test_rows`.
    fn heldout_r2(
        phi: &Array2<f64>,
        target: &Array2<f64>,
        train_rows: &[usize],
        test_rows: &[usize],
    ) -> f64 {
        use faer::Side;
        use gam_linalg::faer_ndarray::FaerCholesky;
        let p = phi.ncols();
        let q = target.ncols();
        let mut gram = Array2::<f64>::zeros((p, p));
        let mut rhs = Array2::<f64>::zeros((p, q));
        for &row in train_rows {
            for a in 0..p {
                for b in 0..p {
                    gram[[a, b]] += phi[[row, a]] * phi[[row, b]];
                }
                for c in 0..q {
                    rhs[[a, c]] += phi[[row, a]] * target[[row, c]];
                }
            }
        }
        let scale = gram.diag().iter().copied().fold(0.0_f64, f64::max);
        for d in gram.diag_mut().iter_mut() {
            *d += scale * 64.0 * f64::EPSILON;
        }
        let decoder = gram.cholesky(Side::Lower).unwrap().solve_mat(&rhs);
        let mut mean = vec![0.0_f64; q];
        for &row in test_rows {
            for c in 0..q {
                mean[c] += target[[row, c]];
            }
        }
        for v in mean.iter_mut() {
            *v /= test_rows.len() as f64;
        }
        let (mut residual, mut total) = (0.0_f64, 0.0_f64);
        for &row in test_rows {
            for c in 0..q {
                let mut fit = 0.0_f64;
                for a in 0..p {
                    fit += phi[[row, a]] * decoder[[a, c]];
                }
                residual += (target[[row, c]] - fit).powi(2);
                total += (target[[row, c]] - mean[c]).powi(2);
            }
        }
        1.0 - residual / total
    }

    /// A band-limited field on `S²` (degree ≤ 3) lies exactly in the real
    /// spherical-harmonic span, so the `Y_l^m` basis reconstructs held-out rows
    /// to near-exactness — while the fixed degree-2 lat/lon chart, which spans
    /// only `[1, x, y, z, xy, yz, xz]` (missing two quadrupoles and every `l ≥ 3`
    /// mode), cannot. That gap is the whole point of #2243's sphere-resolution
    /// follow-up. The bandwidth selector must also recover the planted degree 3
    /// from a noisy over-fit at degree 6. No reference tool is asserted against;
    /// the fixed chart is a baseline we beat on held-out reconstruction.
    #[test]
    fn spherical_harmonic_recovers_bandlimited_field_and_beats_fixed_chart() {
        let mut rng = uniform_stream(0xACE1);
        let n = 3000usize;
        let mut coords = Array2::<f64>::zeros((n, 2));
        for row in 0..n {
            // Uniform on the sphere: lat = asin(U(-1,1)), lon = U(0, 2π).
            coords[[row, 0]] = (2.0 * rng() - 1.0).asin();
            coords[[row, 1]] = std::f64::consts::TAU * rng();
        }
        // Plant a degree-3 field: random decoder over the 16 real harmonics.
        let planted = SphericalHarmonicEvaluator::new(3).unwrap();
        let (basis3, _) = planted.evaluate(coords.view()).unwrap();
        let p_out = 6usize;
        let mut decoder = Array2::<f64>::zeros((basis3.ncols(), p_out));
        for a in 0..basis3.ncols() {
            for c in 0..p_out {
                decoder[[a, c]] = 2.0 * rng() - 1.0;
            }
        }
        let target = basis3.dot(&decoder);

        let test_rows: Vec<usize> = (0..n).filter(|r| r % 4 == 0).collect();
        let train_rows: Vec<usize> = (0..n).filter(|r| r % 4 != 0).collect();

        let sh_r2 = heldout_r2(&basis3, &target, &train_rows, &test_rows);
        let (fixed_phi, _) = SphereChartEvaluator.evaluate(coords.view()).unwrap();
        let fixed_r2 = heldout_r2(&fixed_phi, &target, &train_rows, &test_rows);

        assert!(
            sh_r2 > 0.999,
            "the spherical-harmonic basis must reconstruct a band-limited sphere \
             field to near-exactness; SH held-out R²={sh_r2}, fixed-chart R²={fixed_r2}"
        );
        assert!(
            fixed_r2 < 0.9,
            "the fixed degree-2 chart must NOT reach the higher-degree field \
             (its span omits two quadrupoles and every l≥3 mode); \
             SH held-out R²={sh_r2}, fixed-chart R²={fixed_r2}"
        );
        assert!(
            sh_r2 > fixed_r2,
            "the spherical-harmonic basis must beat the fixed chart on held-out \
             reconstruction; SH held-out R²={sh_r2}, fixed-chart R²={fixed_r2}"
        );

        // Bandwidth selection: a noisy fit at max_degree=6 must cut back to the
        // planted degree 3 at the spectral noise floor.
        let mut noisy = target.clone();
        for row in 0..n {
            for c in 0..p_out {
                noisy[[row, c]] += 0.01 * (2.0 * rng() - 1.0);
            }
        }
        let selected = select_spherical_harmonic_degree(coords.view(), noisy.view(), 6).unwrap();
        assert_eq!(
            selected, 3,
            "spectral-noise-floor bandwidth selection must recover the planted degree 3"
        );
    }

    /// The analytic Hessian and third jet must match finite differences of the
    /// lower jet (the Newton solver consumes them). Order-1 FD of a lower
    /// analytic derivative is well-conditioned, unlike a raw high-order FD of the
    /// value, so this pins every entry: `second_jet` against a central
    /// difference of the `evaluate` jet, and `third_jet` against a central
    /// difference of `second_jet`.
    #[test]
    fn spherical_harmonic_jets_match_finite_differences() {
        let evaluator = SphericalHarmonicEvaluator::new(4).unwrap();
        let coords =
            Array2::from_shape_vec((4, 2), vec![0.3, 0.7, -0.9, 2.1, 1.2, -1.3, -0.1, 4.0])
                .unwrap();
        let h = 1e-5;
        let m = evaluator.basis_size();
        let hess = evaluator.second_jet(coords.view()).unwrap();
        let third = evaluator.third_jet(coords.view()).unwrap();

        let shifted = |axis: usize, step: f64| -> Array2<f64> {
            let mut c = coords.clone();
            for row in 0..c.nrows() {
                c[[row, axis]] += step;
            }
            c
        };

        // second_jet vs central FD of the analytic first jet.
        let mut max_h_err = 0.0_f64;
        for axis in 0..2 {
            let (_, jp) = evaluator.evaluate(shifted(axis, h).view()).unwrap();
            let (_, jm) = evaluator.evaluate(shifted(axis, -h).view()).unwrap();
            for row in 0..coords.nrows() {
                for col in 0..m {
                    for other in 0..2 {
                        let fd = (jp[[row, col, other]] - jm[[row, col, other]]) / (2.0 * h);
                        max_h_err = max_h_err.max((hess[[row, col, other, axis]] - fd).abs());
                    }
                }
            }
        }
        assert!(
            max_h_err < 1e-4,
            "spherical-harmonic Hessian must match FD of the first jet; max err {max_h_err}"
        );

        // third_jet vs central FD of the analytic Hessian.
        let mut max_t_err = 0.0_f64;
        for axis in 0..2 {
            let hp = evaluator.second_jet(shifted(axis, h).view()).unwrap();
            let hm = evaluator.second_jet(shifted(axis, -h).view()).unwrap();
            for row in 0..coords.nrows() {
                for col in 0..m {
                    for a in 0..2 {
                        for b in 0..2 {
                            let fd = (hp[[row, col, a, b]] - hm[[row, col, a, b]]) / (2.0 * h);
                            max_t_err = max_t_err.max((third[[row, col, a, b, axis]] - fd).abs());
                        }
                    }
                }
            }
        }
        assert!(
            max_t_err < 1e-4,
            "spherical-harmonic third jet must match FD of the Hessian; max err {max_t_err}"
        );
    }

    fn assert_machine_close(actual: f64, expected: f64, label: &str) {
        let tolerance = 4096.0 * f64::EPSILON * (1.0 + actual.abs().max(expected.abs()));
        assert!(
            (actual - expected).abs() <= tolerance,
            "{label}: actual={actual:.17e}, expected={expected:.17e}, tolerance={tolerance:.3e}"
        );
    }

    /// Normwise roundoff certificate for a deck-covariant spectral derivative.
    ///
    /// The quotient stores Laplace eigenvalues without the common torus chart
    /// factor `(2*pi)^2`.  A retained mode with eigenvalue `lambda` therefore
    /// has per-coordinate angular frequency at most `2*pi*sqrt(lambda)`, and an
    /// order-`r` derivative has absolute envelope at most that frequency to the
    /// `r`th power.  This envelope is the relevant scale when the exact result
    /// vanishes: evaluating the deck twin forms a phase shifted by an integer
    /// multiple of pi, whose libm argument-reduction error is then amplified by
    /// the derivative frequency even though the rounded output itself is near
    /// zero.  Spherical quotient frequencies are smaller, so the same bound is
    /// conservative for the projective-plane cover.
    fn assert_spectral_deck_close(
        actual: f64,
        expected: f64,
        laplace_eigenvalue: f64,
        derivative_order: i32,
        label: &str,
    ) {
        let angular_frequency = std::f64::consts::TAU * laplace_eigenvalue.sqrt();
        let derivative_envelope = if derivative_order == 0 {
            1.0
        } else {
            angular_frequency.powi(derivative_order)
        };
        // 64 rounded-operation units cover the deck-coordinate map, phase
        // construction/reduction, tensor products, and sign application.  The
        // analytic envelope—not the observed near-zero result—sets the scale.
        let tolerance =
            64.0 * f64::EPSILON * (1.0 + derivative_envelope + actual.abs().max(expected.abs()));
        assert!(
            (actual - expected).abs() <= tolerance,
            "{label}: actual={actual:.17e}, expected={expected:.17e}, spectral envelope={derivative_envelope:.6e}, tolerance={tolerance:.3e}"
        );
    }

    /// Independent oracle for the restriction constructor: evaluate the exact
    /// Reynolds group average `(f(x) + f(gx)) / 2` on every cover column.
    /// Retained columns must be fixed by the average and omitted columns must
    /// be annihilated.  This catches both an incomplete character table and the
    /// invalid use of column selection for a deck action that mixes columns.
    fn assert_group_average_equals_diagonal_restriction(
        evaluator: &QuotientSpectralEvaluator,
        cover: &dyn SaeBasisEvaluator,
        cover_eigenvalues: &[f64],
        coords: &Array2<f64>,
        twins: &Array2<f64>,
    ) {
        let (cover_phi, _) = cover.evaluate(coords.view()).unwrap();
        let (twin_cover_phi, _) = cover.evaluate(twins.view()).unwrap();
        let (quotient_phi, _) = evaluator.evaluate(coords.view()).unwrap();
        assert_eq!(cover_phi.dim(), twin_cover_phi.dim());
        assert_eq!(cover_phi.ncols(), evaluator.cover_width());
        assert_eq!(cover_eigenvalues.len(), evaluator.cover_width());
        assert_eq!(quotient_phi.ncols(), evaluator.basis_size());

        let mut quotient_column = 0usize;
        for cover_column in 0..cover_phi.ncols() {
            let retained =
                evaluator.cover_columns().get(quotient_column).copied() == Some(cover_column);
            for row in 0..cover_phi.nrows() {
                let group_average =
                    0.5 * (cover_phi[[row, cover_column]] + twin_cover_phi[[row, cover_column]]);
                if retained {
                    assert_spectral_deck_close(
                        group_average,
                        cover_phi[[row, cover_column]],
                        cover_eigenvalues[cover_column],
                        0,
                        "trivial-character column fixed by exact group average",
                    );
                    assert_spectral_deck_close(
                        quotient_phi[[row, quotient_column]],
                        group_average,
                        cover_eigenvalues[cover_column],
                        0,
                        "quotient restriction equals exact group average",
                    );
                } else {
                    assert_spectral_deck_close(
                        group_average,
                        0.0,
                        cover_eigenvalues[cover_column],
                        0,
                        "sign-character column annihilated by exact group average",
                    );
                }
            }
            if retained {
                quotient_column += 1;
            }
        }
        assert_eq!(quotient_column, evaluator.basis_size());
    }

    /// If `q(twin(t)) = q(t)`, jets in cover coordinates obey the chain rule
    /// `J(twin)D = J(t)`, `D' H(twin)D = H(t)`, and the analogous order-3
    /// identity.  Checking all four orders pins both invariance and the exact
    /// analytic quotient jets without finite differencing.
    fn assert_quotient_deck_covariance(
        evaluator: &QuotientSpectralEvaluator,
        coords: &Array2<f64>,
        twins: &Array2<f64>,
        deck_jacobian_diagonal: [f64; 2],
    ) {
        let (phi, jet) = evaluator.evaluate(coords.view()).unwrap();
        let (twin_phi, twin_jet) = evaluator.evaluate(twins.view()).unwrap();
        let hessian = evaluator.second_jet(coords.view()).unwrap();
        let twin_hessian = evaluator.second_jet(twins.view()).unwrap();
        let third = evaluator.third_jet(coords.view()).unwrap();
        let twin_third = evaluator.third_jet(twins.view()).unwrap();
        for row in 0..coords.nrows() {
            for column in 0..evaluator.basis_size() {
                assert_spectral_deck_close(
                    twin_phi[[row, column]],
                    phi[[row, column]],
                    evaluator.laplace_eigenvalues()[column],
                    0,
                    "deck-invariant quotient value",
                );
                for axis_a in 0..2 {
                    assert_spectral_deck_close(
                        twin_jet[[row, column, axis_a]] * deck_jacobian_diagonal[axis_a],
                        jet[[row, column, axis_a]],
                        evaluator.laplace_eigenvalues()[column],
                        1,
                        "deck-covariant quotient first jet",
                    );
                    for axis_b in 0..2 {
                        assert_spectral_deck_close(
                            twin_hessian[[row, column, axis_a, axis_b]]
                                * deck_jacobian_diagonal[axis_a]
                                * deck_jacobian_diagonal[axis_b],
                            hessian[[row, column, axis_a, axis_b]],
                            evaluator.laplace_eigenvalues()[column],
                            2,
                            "deck-covariant quotient second jet",
                        );
                        for axis_c in 0..2 {
                            assert_spectral_deck_close(
                                twin_third[[row, column, axis_a, axis_b, axis_c]]
                                    * deck_jacobian_diagonal[axis_a]
                                    * deck_jacobian_diagonal[axis_b]
                                    * deck_jacobian_diagonal[axis_c],
                                third[[row, column, axis_a, axis_b, axis_c]],
                                evaluator.laplace_eigenvalues()[column],
                                3,
                                "deck-covariant quotient third jet",
                            );
                        }
                    }
                }
            }
        }
    }

    fn assert_exact_quotient_spectral_contract(evaluator: &QuotientSpectralEvaluator) {
        let gram = evaluator.function_space_gram();
        let penalty = evaluator.spectral_penalty(2).unwrap();
        assert_eq!(gram.dim(), (evaluator.basis_size(), evaluator.basis_size()));
        assert_eq!(penalty.dim(), gram.dim());
        assert_eq!(evaluator.nullspace_dimension(), 1);
        let mut positive_penalty_diagonal = 0usize;
        for row in 0..evaluator.basis_size() {
            assert!(gram[[row, row]] > 0.0);
            assert_machine_close(
                penalty[[row, row]],
                gram[[row, row]] * evaluator.laplace_eigenvalues()[row].powi(2),
                "restricted spectral penalty",
            );
            if penalty[[row, row]] > 0.0 {
                positive_penalty_diagonal += 1;
            }
            for column in 0..evaluator.basis_size() {
                if row != column {
                    assert_eq!(gram[[row, column]], 0.0);
                    assert_eq!(penalty[[row, column]], 0.0);
                }
            }
        }
        assert_eq!(
            positive_penalty_diagonal,
            evaluator.basis_size() - 1,
            "the exact restricted penalty null space must contain only constants"
        );
    }

    #[test]
    fn quotient_spectral_projective_plane_group_average_jets_penalty_and_null_are_exact() {
        let harmonic_order = 3;
        let evaluator = QuotientSpectralEvaluator::projective_plane(harmonic_order).unwrap();
        assert_eq!(evaluator.quotient_name(), "projective-plane");
        assert_eq!(evaluator.basis_size(), 28);

        let cover = SphericalHarmonicEvaluator::new(2 * harmonic_order).unwrap();
        let cover_modes = cover.spectral_modes();
        assert_eq!(evaluator.cover_width(), cover.basis_size());
        for (quotient_column, &cover_column) in evaluator.cover_columns().iter().enumerate() {
            assert_eq!(
                evaluator.laplace_eigenvalues()[quotient_column],
                cover_modes[cover_column].laplace_eigenvalue
            );
            assert_eq!(
                evaluator.l2_gram_weights()[quotient_column],
                cover_modes[cover_column].l2_gram_weight
            );
        }
        let split = evaluator.phi_eta_split(evaluator.basis_size()).unwrap();
        assert_eq!(split.base_cols.len(), 6, "l=0 and l=2 form the RP2 core");
        assert_exact_quotient_spectral_contract(&evaluator);

        let coords = Array2::from_shape_vec(
            (5, 2),
            vec![-1.2, -2.8, -0.47, -0.31, 0.0, 0.91, 0.63, 2.17, 1.31, 5.4],
        )
        .unwrap();
        let mut twins = coords.clone();
        for row in 0..twins.nrows() {
            twins[[row, 0]] = -twins[[row, 0]];
            twins[[row, 1]] += std::f64::consts::PI;
        }
        let cover_eigenvalues = cover_modes
            .iter()
            .map(|mode| mode.laplace_eigenvalue)
            .collect::<Vec<_>>();
        assert_group_average_equals_diagonal_restriction(
            &evaluator,
            &cover,
            &cover_eigenvalues,
            &coords,
            &twins,
        );
        assert_quotient_deck_covariance(&evaluator, &coords, &twins, [-1.0, 1.0]);
    }

    #[test]
    fn quotient_spectral_klein_group_average_jets_penalty_and_null_are_exact() {
        let num_harmonics = 3;
        let evaluator = QuotientSpectralEvaluator::klein_bottle(num_harmonics).unwrap();
        assert_eq!(evaluator.quotient_name(), "klein-bottle");
        assert_eq!(evaluator.basis_size(), 24);

        let cover = TorusHarmonicEvaluator::new(2, num_harmonics).unwrap();
        let cover_modes = cover.spectral_modes();
        assert_eq!(evaluator.cover_width(), cover.basis_size());
        for (quotient_column, &cover_column) in evaluator.cover_columns().iter().enumerate() {
            assert_eq!(
                evaluator.laplace_eigenvalues()[quotient_column],
                cover_modes[cover_column].laplace_eigenvalue
            );
            assert_eq!(
                evaluator.l2_gram_weights()[quotient_column],
                cover_modes[cover_column].l2_gram_weight
            );
        }
        let split = evaluator.phi_eta_split(evaluator.basis_size()).unwrap();
        assert_eq!(
            split.base_cols.len(),
            7,
            "the standard R4 Klein embedding is one constant plus six coordinates"
        );
        assert_eq!(split.curved_cols.len(), evaluator.basis_size() - 7);
        let mut embedding_category_counts = [0usize; 4];
        for quotient_column in split.base_cols {
            let cover_mode = &cover_modes[evaluator.cover_columns()[quotient_column]];
            let theta = cover_mode.components[0];
            let phi = cover_mode.components[1];
            match (theta.harmonic(), phi.harmonic(), phi.reflection_sign()) {
                (0, 0, 1) => embedding_category_counts[0] += 1,
                (2, 0, 1) => embedding_category_counts[1] += 1,
                (2, 1, 1) => embedding_category_counts[2] += 1,
                (1, 1, -1) => embedding_category_counts[3] += 1,
                identity => panic!("unexpected Klein embedding-base character {identity:?}"),
            }
        }
        assert_eq!(embedding_category_counts, [1, 2, 2, 2]);
        assert_exact_quotient_spectral_contract(&evaluator);

        let coords = Array2::from_shape_vec(
            (5, 2),
            vec![-0.3, -0.41, 0.0, 0.0, 0.17, 0.29, 0.73, -0.88, 1.21, 1.7],
        )
        .unwrap();
        let mut twins = coords.clone();
        for row in 0..twins.nrows() {
            twins[[row, 0]] += 0.5;
            twins[[row, 1]] = -twins[[row, 1]];
        }
        let cover_eigenvalues = cover_modes
            .iter()
            .map(|mode| mode.laplace_eigenvalue)
            .collect::<Vec<_>>();
        assert_group_average_equals_diagonal_restriction(
            &evaluator,
            &cover,
            &cover_eigenvalues,
            &coords,
            &twins,
        );
        assert_quotient_deck_covariance(&evaluator, &coords, &twins, [1.0, -1.0]);
    }

    #[test]
    fn quotient_and_cover_constructors_enforce_minimum_order_and_checked_widths() {
        assert!(projective_plane_basis_size(0).is_err());
        assert!(klein_bottle_basis_size(0).is_err());
        assert!(QuotientSpectralEvaluator::projective_plane(0).is_err());
        assert!(QuotientSpectralEvaluator::klein_bottle(0).is_err());
        assert!(QuotientSpectralEvaluator::klein_bottle(1).is_err());

        let projective_plane = QuotientSpectralEvaluator::projective_plane(1).unwrap();
        assert_eq!(projective_plane.basis_size(), 6);
        assert_eq!(
            projective_plane.basis_size(),
            projective_plane_basis_size(1).unwrap()
        );
        let klein_bottle = QuotientSpectralEvaluator::klein_bottle(2).unwrap();
        assert_eq!(klein_bottle.basis_size(), 13);
        assert_eq!(
            klein_bottle.basis_size(),
            klein_bottle_basis_size(2).unwrap()
        );
        assert_eq!(
            klein_bottle
                .phi_eta_split(klein_bottle.basis_size())
                .unwrap()
                .base_cols
                .len(),
            7
        );

        assert!(projective_plane_basis_size(usize::MAX).is_err());
        assert!(klein_bottle_basis_size(usize::MAX).is_err());
        assert!(QuotientSpectralEvaluator::projective_plane(usize::MAX).is_err());
        assert!(QuotientSpectralEvaluator::klein_bottle(usize::MAX).is_err());
        assert!(SphericalHarmonicEvaluator::new(usize::MAX).is_err());
        assert!(TorusHarmonicEvaluator::new(1, usize::MAX).is_err());
        assert!(TorusHarmonicEvaluator::new(usize::BITS as usize, 1).is_err());
    }

    #[test]
    fn mobius_basis_is_invariant_under_deck_transform() {
        let evaluator = MobiusHarmonicEvaluator::new(3, 2).unwrap();
        let coords = Array2::from_shape_vec(
            (5, 2),
            vec![0.0, -0.8, 0.17, -0.3, 0.51, 0.0, 0.88, 0.4, 1.41, 0.9],
        )
        .unwrap();
        let mut twins = coords.clone();
        for row in 0..twins.nrows() {
            twins[[row, 0]] += 1.0;
            twins[[row, 1]] = -twins[[row, 1]];
        }
        let (phi, _) = evaluator.evaluate(coords.view()).unwrap();
        let (phi_twin, _) = evaluator.evaluate(twins.view()).unwrap();
        let max_error = (&phi - &phi_twin)
            .iter()
            .map(|value| value.abs())
            .fold(0.0_f64, f64::max);
        assert!(
            max_error <= 32.0 * f64::EPSILON,
            "every basis column must descend to the Möbius quotient; max deck error {max_error}"
        );
    }

    /// The in-place `evaluate_into` must reproduce `evaluate` EXACTLY, even when
    /// the target workspace arrives pre-loaded with garbage — the hot Newton loop
    /// reuses one buffer across trials, so any entry the in-place path forgets to
    /// (re)write would carry a stale value forward. Prefilling with a sentinel and
    /// asserting bit-equality against a fresh `evaluate` pins that.
    fn assert_into_matches_evaluate(eval: &dyn SaeBasisEvaluator, coords: &Array2<f64>) {
        let (phi_ref, jet_ref) = eval.evaluate(coords.view()).expect("evaluate");
        let mut phi = Array2::<f64>::from_elem(phi_ref.dim(), 999.0);
        let mut jet = Array3::<f64>::from_elem(jet_ref.dim(), 999.0);
        eval.evaluate_into(&mut phi, &mut jet, coords.view())
            .expect("evaluate_into");
        assert_eq!(
            phi, phi_ref,
            "evaluate_into Φ must equal evaluate Φ exactly"
        );
        assert_eq!(
            jet, jet_ref,
            "evaluate_into jet must equal evaluate jet exactly"
        );
    }

    /// Reusing ONE workspace across two different coordinate sets (the line-search
    /// cadence) must leave no contamination: the second fill must match a fresh
    /// `evaluate` on the second coordinates.
    fn assert_workspace_reuse(
        eval: &dyn SaeBasisEvaluator,
        coords_a: &Array2<f64>,
        coords_b: &Array2<f64>,
    ) {
        let (phi_a, jet_a) = eval.evaluate(coords_a.view()).expect("evaluate a");
        let mut phi = Array2::<f64>::zeros(phi_a.dim());
        let mut jet = Array3::<f64>::zeros(jet_a.dim());
        eval.evaluate_into(&mut phi, &mut jet, coords_a.view())
            .expect("into a");
        assert_eq!(phi, phi_a);
        assert_eq!(jet, jet_a);
        let (phi_b_ref, jet_b_ref) = eval.evaluate(coords_b.view()).expect("evaluate b");
        eval.evaluate_into(&mut phi, &mut jet, coords_b.view())
            .expect("into b (reused workspace)");
        assert_eq!(
            phi, phi_b_ref,
            "reused workspace Φ must not carry stale data"
        );
        assert_eq!(
            jet, jet_b_ref,
            "reused workspace jet must not carry stale data"
        );
    }

    #[test]
    fn periodic_harmonic_evaluate_into_matches() {
        let eval = PeriodicHarmonicEvaluator::new(5).unwrap();
        let coords_a = Array2::from_shape_vec((4, 1), vec![0.10, 0.35, 0.60, 0.85]).unwrap();
        let coords_b = Array2::from_shape_vec((4, 1), vec![0.20, 0.45, 0.70, 0.05]).unwrap();
        assert_into_matches_evaluate(&eval, &coords_a);
        assert_workspace_reuse(&eval, &coords_a, &coords_b);
    }

    #[test]
    fn euclidean_patch_evaluate_into_matches() {
        let eval = EuclideanPatchEvaluator::new(2, 2).unwrap();
        let coords_a =
            Array2::from_shape_vec((4, 2), vec![0.1, -0.2, 0.3, 0.4, -0.5, 0.6, 0.7, -0.8])
                .unwrap();
        let coords_b =
            Array2::from_shape_vec((4, 2), vec![-0.3, 0.9, 0.2, -0.1, 0.5, 0.5, -0.7, 0.3])
                .unwrap();
        assert_into_matches_evaluate(&eval, &coords_a);
        assert_workspace_reuse(&eval, &coords_a, &coords_b);
    }

    #[test]
    fn torus_harmonic_evaluate_into_matches() {
        let eval = TorusHarmonicEvaluator::new(2, 2).unwrap();
        let coords_a =
            Array2::from_shape_vec((4, 2), vec![0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8]).unwrap();
        let coords_b =
            Array2::from_shape_vec((4, 2), vec![0.9, 0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65])
                .unwrap();
        assert_into_matches_evaluate(&eval, &coords_a);
        assert_workspace_reuse(&eval, &coords_a, &coords_b);
    }

    #[test]
    fn default_evaluate_into_matches_for_unspecialized_evaluator() {
        // SphereChartEvaluator does not override `evaluate_into`, so this pins the
        // allocate-and-copy DEFAULT trait method against `evaluate`.
        let eval = SphereChartEvaluator;
        let coords_a =
            Array2::from_shape_vec((3, 2), vec![0.2, 0.5, -0.4, 1.1, 0.9, -0.7]).unwrap();
        let coords_b =
            Array2::from_shape_vec((3, 2), vec![-0.1, 0.3, 0.6, -0.9, -0.5, 0.8]).unwrap();
        assert_into_matches_evaluate(&eval, &coords_a);
        assert_workspace_reuse(&eval, &coords_a, &coords_b);
    }

    #[test]
    fn evaluate_into_rejects_mismatched_buffer() {
        let eval = PeriodicHarmonicEvaluator::new(5).unwrap();
        let coords = Array2::from_shape_vec((4, 1), vec![0.1, 0.2, 0.3, 0.4]).unwrap();
        // Wrong Φ width (4 columns instead of 5) must be rejected, not silently
        // written past — the shape guard `refresh_basis` relied on.
        let mut phi = Array2::<f64>::zeros((4, 4));
        let mut jet = Array3::<f64>::zeros((4, 5, 1));
        assert!(
            eval.evaluate_into(&mut phi, &mut jet, coords.view())
                .is_err()
        );
    }
}