gam-sae 0.3.155

//! #1782 — `sae_manifold_fit` with `threshold_gate`/`softmax` assignments and
//! `euclidean`/`linear` topologies failed at "no candidate seeds passed outer
//! startup validation" on clean planted-circle data where `ordered_beta_bernoulli`+`circle`
//! converges. Root causes: (1) the Euclidean/Linear PCA seed read the SAME
//! leading principal-component scores for EVERY atom, so a K-atom dictionary
//! started as K identical atoms — a rank-deficient joint decoder whose undamped
//! Laplace factor is non-PD (the #1094 seed-startup refusal); (2) the
//! separable-gate (softmax / threshold_gate) seed dispersion-scaling WEAKENED
//! the decoder-smoothness / ARD seed toward zero on clean data, driving the
//! multi-atom joint Hessian indefinite at the seed. Both are recoverable
//! infeasible-ρ refusals that the single-seed EFS startup validation turned into
//! a fatal abort.
//!
//! These tests fit tiny (N=60, p=8, K=4) planted-circle dictionaries — the same
//! shape as the issue's repro but small enough to run in seconds / a few MB under
//! the RAM-tight shared build gate — through the real outer `OuterProblem::run`
//! ("SAE manifold") cascade, and assert each assignment/topology combination now
//! converges to a finite reconstruction EV instead of throwing
//! `RemlConvergenceError`.

#![cfg(test)]

use super::tests::{global_ev, planted_circle_embedded};
use super::*;
use crate::basis::{EuclideanPatchEvaluator, PeriodicHarmonicEvaluator, SaeBasisSecondJet};
use crate::sparse_dict::{SparseDictConfig, fit_sparse_dictionary};
use gam_linalg::faer_ndarray::{FaerCholesky, fast_atb};
use gam_solve::rho_optimizer::OuterObjective;
use ndarray::{Array1, Array2, ArrayView2, array, s};
use std::sync::Arc;

#[derive(Clone, Copy)]
pub(crate) enum Topo {
    Circle,
    Euclidean,
    Linear,
}

/// Build a K-atom, d=1 SAE term seeded exactly the way the production cold path
/// does (PCA-seed the per-atom coordinates, ridge-LSQ each per-atom decoder),
/// for the requested topology and assignment mode. Returns the term and the
/// seed reconstruction dispersion the outer cascade scales its ρ seed by.
pub(crate) fn build_term(
    z: ArrayView2<'_, f64>,
    k: usize,
    topo: Topo,
    mode: AssignmentMode,
) -> (SaeManifoldTerm, f64) {
    let n = z.nrows();
    let (basis_kind, dim, topo_name): (SaeAtomBasisKind, usize, &str) = match topo {
        Topo::Circle => (SaeAtomBasisKind::Periodic, 1, "circle"),
        Topo::Euclidean => (SaeAtomBasisKind::EuclideanPatch, 1, "euclidean"),
        Topo::Linear => (SaeAtomBasisKind::Linear, 1, "linear"),
    };
    let evaluator: Arc<dyn SaeBasisSecondJet> = match topo {
        Topo::Circle => Arc::new(PeriodicHarmonicEvaluator::new(3).unwrap()),
        Topo::Euclidean => Arc::new(EuclideanPatchEvaluator::new(dim, 2).unwrap()),
        Topo::Linear => Arc::new(EuclideanPatchEvaluator::new(dim, 1).unwrap()),
    };
    let basis_kinds = vec![basis_kind.clone(); k];
    let atom_dims = vec![dim; k];
    let seed_coords = sae_pca_seed_initial_coords(z, &basis_kinds, &atom_dims).unwrap();
    let mut atoms = Vec::with_capacity(k);
    let mut coords_blocks = Vec::with_capacity(k);
    let mut manifolds = Vec::with_capacity(k);
    let mut rss = 0.0_f64;
    for atom_idx in 0..k {
        let coords = seed_coords.slice(s![atom_idx, .., 0..dim]).to_owned();
        let (phi, jet) = evaluator.evaluate(coords.view()).unwrap();
        let mm = phi.ncols();
        let mut xtx = fast_atb(&phi, &phi);
        for i in 0..mm {
            xtx[[i, i]] += 1.0e-8;
        }
        let xtz = fast_atb(&phi, &z.to_owned());
        let decoder = xtx.cholesky(Side::Lower).unwrap().solve_mat(&xtz);
        let fitted = phi.dot(&decoder);
        for row in 0..n {
            for col in 0..z.ncols() {
                let r = z[[row, col]] - fitted[[row, col]];
                rss += r * r;
            }
        }
        let atom = SaeManifoldAtom::new_with_provided_function_gram(
            topo_name,
            basis_kind.clone(),
            dim,
            phi,
            jet,
            decoder,
            Array2::<f64>::eye(mm),
        )
        .unwrap()
        .with_basis_evaluator(evaluator.clone());
        atoms.push(atom);
        coords_blocks.push(coords);
        manifolds.push(match topo {
            Topo::Circle => LatentManifold::Circle { period: 1.0 },
            _ => LatentManifold::Euclidean,
        });
    }
    let seed_dispersion = (rss / (k * n * z.ncols()) as f64).max(1.0e-12);
    // Routing seed. ordered Beta--Bernoulli starts every gate on (the production cold seed). The
    // separable gates start from a round-robin row->atom assignment — a stand-in
    // for the deterministic alternating routing refine — so the routing is not degenerately
    // symmetric (every atom carries mass; no atom is a duplicate of another).
    let mut logits = Array2::<f64>::zeros((n, k));
    for row in 0..n {
        for atom in 0..k {
            logits[[row, atom]] = match mode {
                AssignmentMode::OrderedBetaBernoulli { .. } => 6.0,
                AssignmentMode::Softmax { .. } => {
                    if atom == row % k {
                        3.0
                    } else {
                        0.0
                    }
                }
                AssignmentMode::ThresholdGate { .. } => {
                    if atom == row % k {
                        3.0
                    } else {
                        -3.0
                    }
                }
                // TopK routes each row to its highest-logit atoms: the same
                // round-robin favored-atom seed keeps every atom carrying mass
                // (no degenerate symmetry, no duplicate atoms), and the margin
                // makes the per-row top-k selection deterministic.
                AssignmentMode::TopK { .. } => {
                    if atom == row % k {
                        3.0
                    } else {
                        0.0
                    }
                }
            };
        }
    }
    let assignment =
        SaeAssignment::from_blocks_with_mode_and_manifolds(logits, coords_blocks, manifolds, mode)
            .unwrap();
    (
        SaeManifoldTerm::new(atoms, assignment).unwrap(),
        seed_dispersion,
    )
}

/// Build the objective and dispersion-scaled seed ρ exactly the way the FFI
/// does (single seed, `inner_max_iter` short — the seed decoder is already
/// LSQ-fit, so its inner solve starts near-optimal and converges quickly).
pub(crate) fn objective_and_seed(
    z: ArrayView2<'_, f64>,
    k: usize,
    topo: Topo,
    mode: AssignmentMode,
) -> (SaeManifoldOuterObjective, Array1<f64>) {
    let (term, seed_dispersion) = build_term(z, k, topo, mode);
    let init_rho = SaeManifoldRho::new(0.02_f64.ln(), 1.0_f64.ln(), vec![array![0.0]; k])
        .seed_scaled_by_dispersion_for_assignment(seed_dispersion, mode)
        .unwrap();
    let init_rho_flat = init_rho.to_flat();
    let objective =
        SaeManifoldOuterObjective::new(term, z.to_owned(), None, init_rho, 8, 0.04, 1.0e-6, 1.0e-6);
    (objective, init_rho_flat)
}

/// Assert the seed passes the EFS OUTER STARTUP VALIDATION — the exact gate the
/// issue reports failing. For an all-penalty-like objective with `n_params > 8`
/// the outer planner selects the EFS solver, whose seed validation is a single
/// `eval_efs(seed)` (`run_fixed_point_outer_solver`): it must return a finite
/// cost and finite steps. Before the fix a recoverable non-PD-seed / did-not-
/// converge refusal `?`-propagated out of `efs_step` as a fatal error and — with
/// the single SAE seed — surfaced as `RemlConvergenceError`: "no candidate seeds
/// passed outer startup validation (SAE manifold)".
fn seed_passes_startup_validation(
    z: ArrayView2<'_, f64>,
    k: usize,
    topo: Topo,
    mode: AssignmentMode,
) -> Result<f64, String> {
    let (mut objective, seed) = objective_and_seed(z, k, topo, mode);
    // n_params = 1 (sparse) + K (smooth) + K (ARD) = 1 + 2K; K = 4 -> 9 > 8, so
    // the production planner routes this through the EFS lane, whose startup
    // validation is exactly this call.
    assert!(
        seed.len() > 8,
        "test must exercise the EFS lane (n_params={} must exceed 8)",
        seed.len()
    );
    let eval = objective.eval_efs(&seed).map_err(|e| e.to_string())?;
    if eval.cost == f64::INFINITY {
        // `+inf` is not a numerical accident here: `OuterEval::infeasible` uses it
        // as the conventional representation of an infeasible trial point. Name
        // the EFS startup lane that produced the verdict; #2609 demonstrated that
        // attributing it to the criterion was false when a different basin-entry
        // path at the same rho had a finite criterion.
        return Err(
            "seed rejected: the EFS startup lane reports this rho INFEASIBLE \
             (cost = +inf, the conventional infeasible encoding)"
                .to_string(),
        );
    }
    if !eval.cost.is_finite() {
        // NaN or -inf, which the infeasible convention does not cover and which
        // would be a genuine numerical failure.
        return Err(format!(
            "EFS seed cost is numerically invalid ({}) — this is NOT the infeasible              convention and does indicate a divergence",
            eval.cost
        ));
    }
    if let Some((idx, v)) = eval.steps.iter().enumerate().find(|(_, v)| !v.is_finite()) {
        return Err(format!("EFS seed step[{idx}] is non-finite ({v})"));
    }
    Ok(eval.cost)
}

/// #2609 sequential-lane sanity check. The ordinary value lane selects and
/// installs the authoritative basin before EFS evaluates the same seed, so both
/// calls must remain finite. This ordering cannot diagnose the cold EFS entry
/// path: `seed_verdict_depends_on_lane_call_order_2609` below uses independent
/// objectives for that purpose.
///
/// `circle/ordered_beta_bernoulli` is the cell that reports
/// "EFS seed cost is non-finite (inf)" in the matrix below — notable because
/// #1782 records it as the one configuration that used to be the ONLY survivor.
/// This narrows where to look before anyone instruments the criterion.
#[test]
fn efs_and_value_lanes_agree_on_finiteness_at_the_seed_2609() {
    let z = planted_circle_embedded(48, 6, 0.03);
    let (mut objective, seed) = objective_and_seed(
        z.view(),
        4,
        Topo::Circle,
        AssignmentMode::ordered_beta_bernoulli(1.0, 1.0, false),
    );
    let value_lane = objective.eval(&seed).map(|evaluation| evaluation.cost);
    let efs_lane = objective.eval_efs(&seed).map(|evaluation| evaluation.cost);
    eprintln!("2609 value lane = {value_lane:?}");
    eprintln!("2609 efs   lane = {efs_lane:?}");

    let value_cost = value_lane.expect("the value lane must evaluate the seed");
    let efs_cost = efs_lane.expect("the EFS lane must evaluate the seed");
    assert!(
        value_cost.is_finite(),
        "the value lane returned a non-finite seed cost ({value_cost}); the defect is \
         then in the criterion itself, not in the EFS step"
    );
    assert!(
        efs_cost.is_finite(),
        "the value lane installed a finite authoritative basin ({value_cost}) but the \
         following EFS evaluation rejected that same seed ({efs_cost})"
    );
}

/// #2609 value-lane localisation sweep. This varies the ordered-Beta assignment
/// knobs on identical data and records where the authoritative envelope has a
/// finite value; it deliberately does not infer anything about cold EFS entry.
///
/// `alpha` sets the column shapes `a_k = 1/expm1((k+1)·ln(1+1/α))`, `tau` sets
/// the concrete-logit sharpness, and `k` sets how many columns exist — so a
/// divergence that tracks `alpha` implicates the marginal, one that tracks
/// `tau` implicates the logits saturating, and one that tracks `k` implicates
/// the column ladder. Printed as a table rather than asserted, because the
/// point is to say where the boundary IS before anyone claims a cause.
#[test]
fn ordered_beta_finiteness_sweep_2609() {
    let z = planted_circle_embedded(48, 6, 0.03);
    let mut any_finite = false;
    for &k in &[2usize, 4, 8] {
        for &tau in &[0.25_f64, 1.0, 4.0] {
            for &alpha in &[0.25_f64, 1.0, 4.0] {
                let (mut objective, seed) = objective_and_seed(
                    z.view(),
                    k,
                    Topo::Circle,
                    AssignmentMode::ordered_beta_bernoulli(tau, alpha, false),
                );
                let cost = objective.eval(&seed).map(|evaluation| evaluation.cost);
                let verdict = match &cost {
                    Ok(value) if value.is_finite() => {
                        any_finite = true;
                        format!("finite {value:.6e}")
                    }
                    Ok(value) => format!("NON-FINITE {value}"),
                    Err(error) => format!("Err {error}"),
                };
                eprintln!("2609sweep k={k} tau={tau} alpha={alpha}: {verdict}");
            }
        }
    }
    assert!(
        any_finite,
        "every ordered-Beta configuration diverged, so the sweep localises nothing; \
         widen it before concluding the assignment kind itself is at fault"
    );
}

/// The #1782 startup-validation matrix: on identical clean planted-circle data
/// every assignment kind (ordered_beta_bernoulli, softmax, threshold_gate) and every
/// atom topology (circle, euclidean, linear) must PASS outer startup validation.
/// Before the fix only circle/ordered_beta_bernoulli survived; the rest threw "no candidate
/// seeds passed outer startup validation (SAE manifold)". Fast: one inner solve
/// per config from the near-optimal LSQ seed.
#[test]
fn all_assignment_topology_combinations_pass_startup_validation_1782() {
    let z = planted_circle_embedded(48, 6, 0.03);
    let k = 4usize;
    let cases: Vec<(&str, Topo, AssignmentMode)> = vec![
        (
            "circle/ordered_beta_bernoulli",
            Topo::Circle,
            AssignmentMode::ordered_beta_bernoulli(1.0, 1.0, false),
        ),
        ("circle/softmax", Topo::Circle, AssignmentMode::softmax(1.0)),
        (
            "circle/threshold_gate",
            Topo::Circle,
            AssignmentMode::threshold_gate(1.0, 0.0),
        ),
        (
            "euclidean/ordered_beta_bernoulli",
            Topo::Euclidean,
            AssignmentMode::ordered_beta_bernoulli(1.0, 1.0, false),
        ),
        (
            "linear/ordered_beta_bernoulli",
            Topo::Linear,
            AssignmentMode::ordered_beta_bernoulli(1.0, 1.0, false),
        ),
    ];
    // Measure every cell before judging any of them. Unwrapping inside the loop
    // made the FIRST failing cell decide what got measured: `circle/ordered_beta_bernoulli`
    // is case 1, so while #2609 was open the other four cells never ran, and the
    // claim "only that cell fails" was not something this test could support.
    let mut verdicts: Vec<(&str, Result<f64, String>)> = Vec::new();
    for (label, topo, mode) in cases {
        let result = seed_passes_startup_validation(z.view(), k, topo, mode);
        match &result {
            Ok(cost) => eprintln!("REPRO1782 {label}: startup OK (cost={cost:.4e})"),
            Err(e) => eprintln!("REPRO1782 {label}: startup ERR={e}"),
        }
        verdicts.push((label, result));
    }
    let failures: Vec<String> = verdicts
        .iter()
        .filter_map(|(label, result)| {
            result.as_ref().err().map(|e| format!("{label}: {e}"))
        })
        .collect();
    assert!(
        failures.is_empty(),
        "#1782: {} of {} assignment/topology cells failed outer startup validation:\n  {}",
        failures.len(),
        verdicts.len(),
        failures.join("\n  ")
    );
}

/// Run the real outer `OuterProblem::run` ("SAE manifold") cascade — the exact
/// FFI entry — for one topology/assignment pair on the tiny planted-circle
/// fixture, with the single-PCA-seed budget the production `sae_manifold_fit`
/// FFI uses, and return the reconstruction EV. A non-converged best-so-far
/// iterate is still returned as `Ok`, so a returned EV means the fit RAN to a
/// real reconstruction rather than aborting at startup / in the outer solver.
fn run_full_fit(
    z: ArrayView2<'_, f64>,
    k: usize,
    topo: Topo,
    mode: AssignmentMode,
    label: &str,
) -> f64 {
    let (mut objective, seed) = objective_and_seed(z, k, topo, mode);
    let n_params = seed.len();
    let result = gam_solve::rho_optimizer::OuterProblem::new(n_params)
        .with_initial_rho(seed)
        .with_max_iter(4)
        .with_seed_config(gam_problem::SeedConfig {
            max_seeds: 1,
            seed_budget: 1,
            ..Default::default()
        })
        .run(&mut objective, "SAE manifold")
        .unwrap_or_else(|e| {
            // The two #1782 failure surfaces both land here: the threshold-gate / euclidean
            // "no candidate seeds passed outer startup validation" abort, and the
            // softmax "BFGS aborted: globally infeasible neighbourhood at seed
            // (probe-refusal guard)" abort — both are the emptied / globally-refused
            // seed cascade the fit must avoid by entering a basin with defined
            // quasi-Laplace score; infeasible probes remain `+∞` and cannot certify.
            panic!("#1782 {label} fit must not abort at startup / in the outer solver, got: {e}")
        });
    objective
        .certify_outer_result(&result)
        .expect("#1782 outer result must certify the installed state");
    let fitted = objective.into_fitted().expect("outer fit was evaluated");
    let ev = global_ev(z, fitted.term.fitted().view());
    eprintln!("REPRO1782 {label} fit: ev={ev:.4}");
    assert!(
        ev.is_finite(),
        "#1782 {label} produced a non-finite reconstruction EV ({ev})"
    );
    ev
}

/// The assignment axis (the issue's headline: threshold-gate/softmax) must not just
/// pass validation but actually FIT: run the real outer `OuterProblem::run`
/// ("SAE manifold") cascade — the exact FFI entry — on circle atoms for each
/// assignment kind and require a finite reconstruction EV. Circle atoms are
/// well-conditioned, so a low outer-iteration cap keeps this fast; a
/// non-converged best-so-far iterate is still returned as `Ok`, so this asserts
/// the fit RUNS to a real reconstruction rather than aborting at startup.
///
/// `softmax` is the SECOND #1782 failure surface: its seed and its whole
/// neighbourhood land in the recoverable infeasible-ρ refusal class, so the
/// outer BFGS lane previously returned `+∞` for every probe, never accepted a
/// step, and the bridge's non-termination guard escalated the globally-refused
/// neighbourhood to a FATAL seed rejection ("BFGS aborted: globally infeasible
/// neighbourhood at seed (probe-refusal guard)"). `ordered_beta_bernoulli`+`circle` lands in
/// the PD region and never trips it — RED before the fix on `softmax`, GREEN
/// after (the entry path now reaches a basin with defined quasi-Laplace score).
#[test]
fn assignment_kinds_fit_on_circle_1782() {
    let z = planted_circle_embedded(48, 6, 0.03);
    let k = 4usize;
    for (label, mode) in [
        (
            "circle/ordered_beta_bernoulli",
            AssignmentMode::ordered_beta_bernoulli(1.0, 1.0, false),
        ),
        ("circle/softmax", AssignmentMode::softmax(1.0)),
        (
            "circle/threshold_gate",
            AssignmentMode::threshold_gate(1.0, 0.0),
        ),
    ] {
        run_full_fit(z.view(), k, Topo::Circle, mode, label);
    }
}

/// The topology axis of #1782: on identical clean planted-circle data the
/// `euclidean` and `linear` atom topologies (whose rank-deficient PCA seed lands
/// in the recoverable infeasible-ρ refusal class) must also FIT through the real
/// outer cascade, not abort with an emptied / globally-refused seed cascade. Same
/// single-PCA-seed budget and low outer-iteration cap as the assignment-axis
/// test, so it stays fast. RED before the fix (`euclidean`/`linear` aborted at
/// "no candidate seeds passed outer startup validation"); GREEN after.
#[test]
fn topologies_fit_on_circle_data_1782() {
    let z = planted_circle_embedded(48, 6, 0.03);
    let k = 4usize;
    for (label, topo) in [
        ("euclidean/ordered_beta_bernoulli", Topo::Euclidean),
        ("linear/ordered_beta_bernoulli", Topo::Linear),
    ] {
        run_full_fit(
            z.view(),
            k,
            topo,
            AssignmentMode::ordered_beta_bernoulli(1.0, 1.0, false),
            label,
        );
    }
}

/// FRONTIER PROBE (#1026 co-collapse). Measure the largest OVERCOMPLETE
/// dictionary size `K` whose PCA-diversified cold seed still passes the EFS
/// outer startup validation on a single planted circle embedded in `p = 24`
/// dims — the regime where every extra atom competes for the same rank-2 signal
/// and used to co-collapse into a rank-deficient, non-PD seed. The seed
/// diversification (disjoint per-atom PC windows, `pca_seed.rs`) plus the `K > 1`
/// non-PD dispersion floor (`rho.rs`) are exactly the fixes this probes. Prints
/// the pass/fail frontier so the effect is visible in `--nocapture`. Asserts
/// only the known-safe `K = 4` so it can never red the shared tree; the printed
/// frontier drives the co-collapse-saddle work. K is capped small (the dense
/// per-seed arrow-Schur criterion is `O((K·b·p)³)`) so the sweep stays fast.
#[test]
fn cocollapse_startup_frontier_1026() {
    let z = planted_circle_embedded(96, 10, 0.03);
    let ks = [4usize, 8];
    // Compare the assignment modes: ordered Beta--Bernoulli couples all rows through a cross-row
    // Woodbury evidence with NO matrix-free log-det route (so large-K refuses on
    // the dense reduced Schur), whereas the smooth logistic threshold gate is
    // per-row independent and streams. This measures which mode
    // extends the startup frontier, decoupling the routing wall from seed
    // co-collapse.
    let modes: [(&str, fn() -> AssignmentMode); 3] = [
        ("ordered_beta_bernoulli    ", || {
            AssignmentMode::ordered_beta_bernoulli(1.0, 1.0, false)
        }),
        ("thresh_gate", || AssignmentMode::threshold_gate(1.0, 0.5)),
        ("softmax    ", || AssignmentMode::softmax(1.0)),
    ];
    let mut ordered_beta_bernoulli_frontier = 0usize;
    for (label, mk) in modes {
        let mut frontier = 0usize;
        for &k in &ks {
            match seed_passes_startup_validation(z.view(), k, Topo::Circle, mk()) {
                Ok(cost) => {
                    eprintln!("FRONTIER1026 {label} K={k:>3}: startup PASS (cost={cost:.4e})");
                    frontier = k;
                }
                Err(e) => {
                    eprintln!("FRONTIER1026 {label} K={k:>3}: startup FAIL ({e})");
                    break;
                }
            }
        }
        eprintln!("FRONTIER1026 {label}: largest passing K = {frontier}");
        if label.trim() == "ordered_beta_bernoulli" {
            ordered_beta_bernoulli_frontier = frontier;
        }
    }
    assert!(
        ordered_beta_bernoulli_frontier >= 4,
        "startup validation must hold at least to K=4 (got frontier {ordered_beta_bernoulli_frontier})"
    );
}

/// WIN artifact (#1026 / #1610). A PRINCIPLED joint manifold SAE — curved 1-D
/// circle fibers, smooth logistic gate (`threshold_gate`, the per-row
/// streaming assignment whose criterion log-det takes the matrix-free SLQ route)
/// — fit end-to-end by the real outer penalized quasi-Laplace
/// cascade must MATCH-OR-BEAT a traditional linear SAE (`fit_sparse_dictionary`,
/// the "large linear SAE" of #1026) at matched, OVERCOMPLETE dictionary size K on
/// genuinely curved data. On a planted circle a linear dictionary is rank-capped
/// while curved atoms bend to the ring, so the manifold decisively wins. This is
/// the joint-solve WIN (no alternating-minimization searcher, no Python): the
/// coupled inner arrow-Schur Newton exercises the landed disjoint-PC seed
/// diversification and the spectral Schur PD-floor that keep the overcomplete
/// (K > true-rank) joint block PD instead of co-collapsing. K is kept box-safe
/// here; the per-row work is `top_k`-bounded, so it is the same solve at larger
/// K (the streaming matrix-free criterion log-det, exercised by the outer penalized quasi-Laplace
/// cascade, is what carries it to K=32,000).
#[test]
fn manifold_beats_linear_joint_streaming_1026() {
    let z = planted_circle_embedded(120, 10, 0.03);
    for &k in &[8usize] {
        // Traditional linear SAE baseline at matched K (the sparse-dict lane is f32).
        let z32 = z.mapv(|v| v as f32);
        let lin = fit_sparse_dictionary(z32.view(), &SparseDictConfig::new(k))
            .expect("linear SAE baseline fits");
        let ev_linear = lin.explained_variance;

        // Principled joint manifold SAE: curved circle fibers, smooth logistic threshold gate,
        // solved directly by the coupled inner arrow-Schur joint Newton over the
        // (coords t, decoders β) block — the exact joint solve the outer penalized quasi-Laplace
        // cascade drives, run here at a fixed penalty seed so the comparison is a
        // fast, deterministic reconstruction check (no per-step criterion log-det).
        let mode = AssignmentMode::threshold_gate(1.0, 0.0);
        let (mut term, _disp) = build_term(z.view(), k, Topo::Circle, mode);
        let mut rho = SaeManifoldRho::new(
            1.0e-3_f64.ln(),
            1.0e-3_f64.ln(),
            vec![array![1.0e-3_f64.ln()]; k],
        );
        term.run_joint_fit_arrow_schur(z.view(), &mut rho, None, 24, 1.0, 1.0e-6, 1.0e-6)
            .unwrap_or_else(|e| {
                panic!("#1026 manifold K={k} joint inner fit must run e2e, got: {e}")
            });
        let fitted = term.try_fitted().expect("manifold fitted");
        let ev_manifold = global_ev(z.view(), fitted.view());

        eprintln!(
            "WIN1026 K={k:>3}: manifold EV={ev_manifold:.4}  linear EV={ev_linear:.4}  \
             margin={:+.4}",
            ev_manifold - ev_linear
        );
        assert!(
            ev_manifold.is_finite() && ev_linear.is_finite(),
            "#1026 K={k}: both EVs must be finite (manifold={ev_manifold}, linear={ev_linear})"
        );
        // Match-or-beat contract (#1026 strict generalization): the curved
        // dictionary generalizes the linear one, so it must never do worse.
        assert!(
            ev_manifold + 5.0e-2 >= ev_linear,
            "#1026 K={k}: principled manifold SAE must match-or-beat linear \
             (manifold={ev_manifold:.4} vs linear={ev_linear:.4})"
        );
    }
}

/// #2502 — the d>=2 arm at MEASURED parameter parity, against a REAL sphere.
///
/// The published "d>=2 loses decisively" number is unsupported for two reasons
/// with one root: it ASSERTED a parameter count instead of measuring one. The
/// sphere atom was built as a cylinder when the arm ran, and parity was derived
/// from an assumed 7-column `(lat, lon)` chart width. The pole-free ambient
/// basis carries the full `l <= 2` harmonics -- NINE columns -- so the arm
/// labelled "param parity" at K=1818 in fact carried ~18% more decoder
/// parameters than the linear arm it lost to.
///
/// So this measures. Decoder width is read off `phi.ncols()` per atom, which is
/// the realized column count of the basis that was actually built; the mixed
/// arm's `K` is then derived from that measurement. No width is assumed.
///
/// Reconstruction is at the ridge-LSQ SEED (the same construction
/// `objective_and_seed` scores), evaluated on a DISJOINT split: coordinates are
/// re-seeded on held-out rows and the TRAIN decoder is frozen. That is a
/// generalization measure, but of the seed rather than of the converged outer
/// fit -- it does not reproduce the campaign's 250-cycle EV and must not be
/// compared against those numbers.
///
/// KNOWN CONFOUND -- read the printed EVs with this in mind. The aggregation
/// here AVERAGES all `K` atoms' full-data fits, whereas a real SAE routes each
/// row to its top-`s` active atoms. Averaging rewards HOMOGENEITY as such: many
/// similar linear atoms average coherently, while a heterogeneous portfolio
/// whose atoms carry different coordinate systems averages incoherently. So the
/// linear arm's margin here is inflated by an unknown amount, and the mixed
/// arm's EV can even go NEGATIVE (worse than the mean) for a reason that is
/// about the aggregation, not the parameterization.
///
/// What this test therefore establishes is the PARITY MACHINERY -- a measured,
/// skew-guarded, real-ambient-sphere comparison harness -- plus the measured
/// per-atom cost ratio, which is a property of the bases alone and is NOT
/// affected by the confound. The EV ORDERING needs routed reconstruction before
/// it can carry the weight the published table put on it.
#[test]
fn d2_portfolio_loses_at_measured_parameter_parity_2502() {
    // Arm construction: per-atom basis kinds/dims, so the mixed portfolio is a
    // genuine mixture rather than a homogeneous fit under a mixed label.
    fn arm(
        z: ArrayView2<'_, f64>,
        z_test: ArrayView2<'_, f64>,
        kinds: &[SaeAtomBasisKind],
        dims: &[usize],
        evals: &[Arc<dyn SaeBasisSecondJet>],
    ) -> (usize, f64) {
        let k = kinds.len();
        let p = z.ncols();
        let seed = sae_pca_seed_initial_coords(z, kinds, dims).expect("train seed");
        let seed_test = sae_pca_seed_initial_coords(z_test, kinds, dims).expect("test seed");
        let mut params = 0usize;
        let n_test = z_test.nrows();
        let mut recon = Array2::<f64>::zeros((n_test, p));
        for a in 0..k {
            let d = dims[a];
            let coords = seed.slice(s![a, .., 0..d]).to_owned();
            let (phi, _) = evals[a].evaluate(coords.view()).expect("train phi");
            let mm = phi.ncols();
            // MEASURED cost of this atom: its realized decoder is mm x p.
            params += mm * p;
            let mut xtx = fast_atb(&phi, &phi);
            for i in 0..mm {
                xtx[[i, i]] += 1.0e-8;
            }
            let xtz = fast_atb(&phi, &z.to_owned());
            let decoder = xtx
                .cholesky(Side::Lower)
                .expect("seed gram is PD")
                .solve_mat(&xtz);
            // Held-out: re-seed coordinates, FREEZE the train decoder.
            let coords_t = seed_test.slice(s![a, .., 0..d]).to_owned();
            let (phi_t, _) = evals[a].evaluate(coords_t.view()).expect("test phi");
            let fitted = phi_t.dot(&decoder);
            for row in 0..n_test {
                for col in 0..p {
                    recon[[row, col]] += fitted[[row, col]] / k as f64;
                }
            }
        }
        let mut resid = 0.0_f64;
        let mut total = 0.0_f64;
        let mean = z_test.mean_axis(ndarray::Axis(0)).expect("test mean");
        for row in 0..n_test {
            for col in 0..p {
                let r = z_test[[row, col]] - recon[[row, col]];
                resid += r * r;
                let c = z_test[[row, col]] - mean[col];
                total += c * c;
            }
        }
        (params, 1.0 - resid / total.max(1.0e-30))
    }

    let z = planted_circle_embedded(600, 8, 0.03);
    let z_test = planted_circle_embedded(400, 8, 0.03);

    let linear_eval: Arc<dyn SaeBasisSecondJet> =
        Arc::new(EuclideanPatchEvaluator::new(1, 1).expect("linear evaluator"));
    let euclid_eval: Arc<dyn SaeBasisSecondJet> =
        Arc::new(EuclideanPatchEvaluator::new(2, 2).expect("euclidean evaluator"));
    let periodic_eval: Arc<dyn SaeBasisSecondJet> =
        Arc::new(PeriodicHarmonicEvaluator::new(3).expect("periodic evaluator"));
    let sphere_eval: Arc<dyn SaeBasisSecondJet> =
        Arc::new(crate::basis::AmbientSphereHarmonicEvaluator::new(2).expect("sphere evaluator"));

    // Stage 1 -- MEASURE per-atom cost of each arm on a small probe.
    let probe = 10usize;
    let lin_kinds = vec![SaeAtomBasisKind::Linear; probe];
    let lin_dims = vec![1usize; probe];
    let lin_evals = vec![linear_eval.clone(); probe];
    let (lin_probe_params, _) = arm(z.view(), z_test.view(), &lin_kinds, &lin_dims, &lin_evals);

    // The published portfolio, `atom % 5`: linear, euclidean, periodic, sphere,
    // sphere -- whose sphere entries NOW resolve to the ambient harmonic basis.
    let cycle: [(SaeAtomBasisKind, usize, usize); 5] = [
        (SaeAtomBasisKind::Linear, 1, 0),
        (SaeAtomBasisKind::EuclideanPatch, 2, 1),
        (SaeAtomBasisKind::Periodic, 1, 2),
        (SaeAtomBasisKind::Sphere, 3, 3),
        (SaeAtomBasisKind::Sphere, 3, 3),
    ];
    let evals_by_slot = [
        linear_eval.clone(),
        euclid_eval.clone(),
        periodic_eval.clone(),
        sphere_eval.clone(),
    ];
    let mix = |k: usize| {
        let mut kinds = Vec::with_capacity(k);
        let mut dims = Vec::with_capacity(k);
        let mut evals: Vec<Arc<dyn SaeBasisSecondJet>> = Vec::with_capacity(k);
        for a in 0..k {
            let (kind, d, slot) = cycle[a % cycle.len()].clone();
            kinds.push(kind);
            dims.push(d);
            evals.push(evals_by_slot[slot].clone());
        }
        (kinds, dims, evals)
    };
    let (pk, pd, pe) = mix(probe);
    let (mix_probe_params, _) = arm(z.view(), z_test.view(), &pk, &pd, &pe);

    let lin_per_atom = lin_probe_params as f64 / probe as f64;
    let mix_per_atom = mix_probe_params as f64 / probe as f64;
    let cost_ratio = mix_per_atom / lin_per_atom;
    println!(
        "[2502-parity] measured per-atom decoder params: linear={lin_per_atom:.1} \
         mixed={mix_per_atom:.1}  ratio={cost_ratio:.3}x"
    );

    // Stage 2 -- the comparison at the K the MEASUREMENT chose.
    let k_lin = 60usize;
    let k_mix = ((k_lin as f64) / cost_ratio).round().max(1.0) as usize;
    let lk = vec![SaeAtomBasisKind::Linear; k_lin];
    let ld = vec![1usize; k_lin];
    let le = vec![linear_eval.clone(); k_lin];
    let (lin_params, lin_ev) = arm(z.view(), z_test.view(), &lk, &ld, &le);
    let (mk, md, me) = mix(k_mix);
    let (mix_params, mix_ev) = arm(z.view(), z_test.view(), &mk, &md, &me);

    let skew = (mix_params as f64 - lin_params as f64).abs() / (lin_params as f64);
    println!(
        "[2502-parity] linear K={k_lin} params={lin_params} heldout_EV={lin_ev:.4}"
    );
    println!(
        "[2502-parity] mixed  K={k_mix} params={mix_params} heldout_EV={mix_ev:.4}"
    );
    println!(
        "[2502-parity] realized skew={:.1}%  linear_EV - mixed_EV = {:+.4}",
        skew * 100.0,
        lin_ev - mix_ev
    );

    // The guard the original comparison lacked: refuse to call this parity if it
    // is not. A skew this large is what made the K=1818 arm not-a-parity-arm.
    assert!(
        skew <= 0.05,
        "this is NOT a parity comparison: realized skew {:.1}% (linear {} vs mixed {} \
         decoder params). Do not report an EV difference across it.",
        skew * 100.0,
        lin_params,
        mix_params
    );
    assert!(
        lin_ev.is_finite() && mix_ev.is_finite(),
        "both arms must produce a finite held-out EV; got linear={lin_ev} mixed={mix_ev}"
    );
}


/// #2609 attribution. Both lanes were already shown to agree on `Ok(+inf)`, and
/// `+inf` is the documented infeasible encoding — so the open question is WHICH
/// of the four channels that can emit it actually fires here. `eval` collapses
/// all four into one value, so reading the value cannot answer it, and
/// `log::debug!` cannot either: a `--lib` test binary installs no logger backend.
///
/// This calls the criterion the objective calls, on a term this test builds
/// itself, and matches the TYPED outcome:
///   * `Err(VanishedAtoms)`          — fixed-K structural boundary
///   * `Err(IndefiniteObservedInformation)` — `½log|A|` undefined at the mode
///   * `Err(Numerical)`              — a genuine defect (never mapped to `+inf`)
///   * `Ok(non-finite)`              — the assembled-value class
/// and prints the same readout for the two assignments that PASS on identical
/// data, so the discriminating channel is named rather than inferred.
#[test]
fn seed_infeasibility_channel_is_named_2609() {
    let z = planted_circle_embedded(48, 6, 0.03);
    let cases: [(&str, Topo, AssignmentMode); 5] = [
        (
            "circle/ordered_beta_bernoulli",
            Topo::Circle,
            AssignmentMode::ordered_beta_bernoulli(1.0, 1.0, false),
        ),
        ("circle/softmax", Topo::Circle, AssignmentMode::softmax(1.0)),
        (
            "circle/threshold_gate",
            Topo::Circle,
            AssignmentMode::threshold_gate(1.0, 0.0),
        ),
        (
            "euclidean/ordered_beta_bernoulli",
            Topo::Euclidean,
            AssignmentMode::ordered_beta_bernoulli(1.0, 1.0, false),
        ),
        (
            "linear/ordered_beta_bernoulli",
            Topo::Linear,
            AssignmentMode::ordered_beta_bernoulli(1.0, 1.0, false),
        ),
    ];
    let mut fatal_channels = Vec::new();
    for (label, topo, mode) in cases {
        let (mut term, seed_dispersion) = build_term(z.view(), 4, topo, mode);
        let rho = SaeManifoldRho::new(0.02_f64.ln(), 1.0_f64.ln(), vec![array![0.0]; 4])
            .seed_scaled_by_dispersion_for_assignment(seed_dispersion, mode)
            .unwrap();
        let warm = term.warm_start_latents_from_amortized_encoder(z.view(), &rho);
        let outcome = term.penalized_quasi_laplace_criterion_with_refine_policy_and_lane(
            z.view(),
            &rho,
            None,
            8,
            0.04,
            1.0e-6,
            1.0e-6,
            true,
            None,
        );
        let channel = match &outcome {
            Ok((cost, _)) if cost.is_finite() => format!("FINITE cost={cost:.6e}"),
            Ok((cost, _)) => format!("ASSEMBLED-NON-FINITE cost={cost}"),
            Err(SaeCriterionError::VanishedAtoms(atoms)) => {
                format!("VANISHED-ATOMS {atoms}")
            }
            Err(SaeCriterionError::IndefiniteObservedInformation { block }) => {
                format!("INDEFINITE-OBSERVED-INFORMATION block={block}")
            }
            Err(SaeCriterionError::Numerical(message)) => {
                format!("NUMERICAL(fatal, never mapped to +inf) {message}")
            }
        };
        println!("[2609-channel] {label}: {channel}  warm_start={warm:?}");
        if matches!(&outcome, Err(SaeCriterionError::Numerical(_))) {
            fatal_channels.push((label, outcome.err()));
        }
    }
    // The invariant this run must not violate: the criterion is being asked at a
    // seed the production planner would hand it, so a `Numerical` refusal here
    // would be a defect that never reaches the `+inf` convention at all. Every
    // other outcome is a legitimate verdict and is reported above.
    assert!(
        fatal_channels.is_empty(),
        "seed criteria returned fatal Numerical refusals, which the +inf infeasible \
         convention does not cover: {fatal_channels:?}"
    );
}

/// #2609 — the seed verdict depends on the ORDER of the two lane calls, which is
/// why the lane bisect could not see the failure it was built to attribute.
///
/// `efs_and_value_lanes_agree_on_finiteness_at_the_seed_2609` calls `eval` and
/// then `eval_efs` on ONE objective. A finite `eval` selects the authoritative
/// lower-envelope basin and leaves its converged term installed, so the next
/// evaluation at the same ρ opens AT the inner KKT optimum instead of tracing
/// from the cold LSQ seed. Before #2609, the first call could therefore heal the
/// second, and "both lanes agree" described mutable call history rather than the
/// point being evaluated.
///
/// Startup validation (`seed_passes_startup_validation`) calls `eval_efs` on a
/// FRESH objective with nothing parked. This measures both orders on separate
/// objectives so the two are not the same experiment:
///   * `A` — `eval_efs` alone, cold. This is what startup validation does.
///   * `B` — `eval` then `eval_efs`, one objective. This is what the bisect does.
/// Both orders must now select the same authoritative basin before EFS commits
/// its accepted state; clearing the pre-commit bundle is legal only after that
/// selection.
#[test]
fn seed_verdict_depends_on_lane_call_order_2609() {
    let z = planted_circle_embedded(48, 6, 0.03);
    let mode = AssignmentMode::ordered_beta_bernoulli(1.0, 1.0, false);

    let (mut cold, seed_a) = objective_and_seed(z.view(), 4, Topo::Circle, mode);
    let a_efs = cold.eval_efs(&seed_a).map(|evaluation| evaluation.cost);
    eprintln!("[2609-order] A cold eval_efs FIRST      = {a_efs:?}");

    let (mut warm, seed_b) = objective_and_seed(z.view(), 4, Topo::Circle, mode);
    let b_value = warm.eval(&seed_b).map(|evaluation| evaluation.cost);
    let b_efs = warm.eval_efs(&seed_b).map(|evaluation| evaluation.cost);
    eprintln!("[2609-order] B eval THEN eval_efs       = {b_value:?} then {b_efs:?}");

    let a_cost = a_efs.expect("the cold EFS lane must evaluate the seed");
    let b_value_cost = b_value.expect("the value lane must evaluate the seed");
    let b_efs_cost = b_efs.expect("the warmed EFS lane must evaluate the seed");
    eprintln!(
        "[2609-order] cold_efs_finite={} value_finite={} warmed_efs_finite={}",
        a_cost.is_finite(),
        b_value_cost.is_finite(),
        b_efs_cost.is_finite()
    );

    // This planted seed has a measured finite authoritative envelope. Both EFS
    // call orders must therefore install that basin and return a finite value;
    // merely making both orders agree on +inf would preserve the startup defect.
    assert!(
        b_value_cost.is_finite(),
        "the authoritative value lane must establish that this seed is feasible, got {b_value_cost}"
    );
    assert!(
        a_cost.is_finite(),
        "cold EFS rejected a seed whose authoritative envelope is finite ({b_value_cost}); \
         startup validation must select the basin before committing it"
    );
    assert!(
        b_efs_cost.is_finite(),
        "EFS rejected the same seed after a finite value evaluation: {b_efs_cost}"
    );
}