gam-sae 0.3.155

#![cfg(test)]
//! Stationary-cache `∂log|H|/∂θ` adjoint regression tests (#1416),
//! split verbatim out of `tests.rs` to keep that tracked file under the #780
//! 10k-line gate. Declared as a sibling `#[cfg(test)] mod` in `mod.rs`; shared
//! `gamma_fd_tiny_fixture` / `fixed_state_logdet_sample` are sourced from the sibling
//! `tests` module.

#![cfg(test)]

use super::*;

#[derive(Clone, Copy)]
struct TinyComplex {
    re: f64,
    im: f64,
}

impl TinyComplex {
    fn real(re: f64) -> Self {
        Self { re, im: 0.0 }
    }

    fn add(self, other: Self) -> Self {
        Self {
            re: self.re + other.re,
            im: self.im + other.im,
        }
    }

    fn mul(self, other: Self) -> Self {
        Self {
            re: self.re * other.re - self.im * other.im,
            im: self.re * other.im + self.im * other.re,
        }
    }

    fn div(self, other: Self) -> Self {
        let denom = other.re * other.re + other.im * other.im;
        Self {
            re: (self.re * other.re + self.im * other.im) / denom,
            im: (self.im * other.re - self.re * other.im) / denom,
        }
    }

    fn exp(self) -> Self {
        let e = self.re.exp();
        Self {
            re: e * self.im.cos(),
            im: e * self.im.sin(),
        }
    }
}

fn real_softmax(logits: &[f64], tau: f64) -> Vec<f64> {
    let max_logit = logits.iter().copied().fold(f64::NEG_INFINITY, f64::max);
    let mut weights: Vec<f64> = logits
        .iter()
        .map(|&z| ((z - max_logit) / tau).exp())
        .collect();
    let sum: f64 = weights.iter().sum();
    for weight in weights.iter_mut() {
        *weight /= sum;
    }
    weights
}

fn complex_softmax_weight_product_derivative(
    logits: &[f64],
    tau: f64,
    atom_a: usize,
    atom_b: usize,
    atom_w: usize,
    block_inner: f64,
) -> f64 {
    let h = 1.0e-30;
    let max_logit = logits.iter().copied().fold(f64::NEG_INFINITY, f64::max);
    let mut denom = TinyComplex::real(0.0);
    let mut numer_a = TinyComplex::real(0.0);
    let mut numer_b = TinyComplex::real(0.0);
    for (atom, &logit) in logits.iter().enumerate() {
        let z = TinyComplex {
            re: (logit - max_logit) / tau,
            im: if atom == atom_w { h / tau } else { 0.0 },
        };
        let exp_z = z.exp();
        denom = denom.add(exp_z);
        if atom == atom_a {
            numer_a = exp_z;
        }
        if atom == atom_b {
            numer_b = exp_z;
        }
    }
    let a = numer_a.div(denom);
    let b = numer_b.div(denom);
    a.mul(b).mul(TinyComplex::real(block_inner)).im / h
}

#[test]
pub(crate) fn softmax_tt_weight_product_logit_adjoint_hits_both_factors_2156() {
    let logits = [0.31_f64, -0.27, 0.14, -0.08];
    let tau = 0.73_f64;
    let inv_tau = 1.0 / tau;
    let assignments = real_softmax(&logits, tau);
    let block_inner = 1.417_f64;

    for (atom_a, atom_b, atom_w) in [(0usize, 2usize, 1usize), (2usize, 2usize, 2usize)] {
        let h_ab = assignments[atom_a] * assignments[atom_b] * block_inner;
        let one_factor =
            h_ab * (if atom_w == atom_a { 1.0 } else { 0.0 } - assignments[atom_w]) * inv_tau;
        let fixed = h_ab
            * SaeManifoldTerm::softmax_data_weight_product_logit_factor(
                &assignments,
                atom_a,
                atom_b,
                atom_w,
                inv_tau,
            );
        let complex_step = complex_softmax_weight_product_derivative(
            &logits,
            tau,
            atom_a,
            atom_b,
            atom_w,
            block_inner,
        );
        let ratio = fixed / one_factor;
        assert!(
            (ratio - 2.0).abs() <= 1.0e-12,
            "one-factor softmax product derivative must be 2x low: got ratio {ratio:.12}"
        );
        assert!(
            (fixed - complex_step).abs() <= 1.0e-6 * (1.0 + complex_step.abs()),
            "fixed softmax product derivative must match complex-step: fixed={fixed:.12e}, complex={complex_step:.12e}"
        );
    }
}

// #2330 Patch D — fixed-θ EXACT-A logdet for the θ-adjoint FD arbiter: rebuild
// the fixed-θ̂ cache at the (perturbed) state and return log|A| (not log|B|).
// `None` when the criterion refuses or A is indefinite there (so the FD probe
// can report that instead of panicking).
fn fixed_state_exact_a_logdet(
    mut term: SaeManifoldTerm,
    target: &Array2<f64>,
    rho: &SaeManifoldRho,
) -> Option<f64> {
    let (_v, _l, cache) = term
        .penalized_quasi_laplace_criterion_with_cache(
            target.view(),
            rho,
            None,
            0,
            0.4,
            1.0e-6,
            1.0e-6,
        )
        .ok()?;
    term.exact_observed_information_log_dets(rho, target.view(), &cache)
        .ok()
        .map(|(log_a, _log_a_tt)| log_a)
}

// #2330 Patch D — an ordered-Beta--Bernoulli fixture whose target is generated
// with the SAME independent-logistic gates the model applies. The shared
// `gamma_fd_tiny_fixture` builds its target from NORMALIZED softmax weights, so
// simply flipping that fixture's mode to ordered Beta--Bernoulli leaves a target
// the model cannot reach: the resulting large residual drives the dropped
// residual curvature `ΔC = ⟨error_metric, ∂²f⟩` big enough to push the exact
// `A = B + ΔC` indefinite, and the Phase-2a criterion then refuses at
// construction. `residual_scale` adds a deterministic model-unreachable
// component on top of the reachable target, so `ΔC` — the object Patch D
// differentiates — is nonzero and tunable rather than either zero (a fixture
// that would false-green the arbiter) or saddle-inducing.
pub(crate) fn obb_patchd_fixture(
    residual_scale: f64,
    log_lambda_sparse: f64,
) -> (SaeManifoldTerm, Array2<f64>, SaeManifoldRho) {
    let n = 10usize;
    let p = 3usize;
    let k_atoms = 2usize;
    let m = 3usize;
    let tau = 0.7_f64;
    let evaluator = Arc::new(PeriodicHarmonicEvaluator::new(m).unwrap());
    let mut logits = Array2::<f64>::zeros((n, k_atoms));
    let mut coords = vec![Array2::<f64>::zeros((n, 1)), Array2::<f64>::zeros((n, 1))];
    let weights = [
        [
            [0.10, -0.05, 0.03],
            [0.35, -0.20, 0.12],
            [-0.16, 0.18, 0.08],
        ],
        [
            [-0.08, 0.04, 0.06],
            [0.22, 0.10, -0.18],
            [0.11, -0.24, 0.15],
        ],
    ];
    let mut target = Array2::<f64>::zeros((n, p));
    for row in 0..n {
        let phase = (row as f64 + 0.35) / n as f64;
        coords[0][[row, 0]] = phase;
        coords[1][[row, 0]] = (phase + 0.21).fract();
        logits[[row, 0]] = if row % 2 == 0 { 0.8 } else { -0.6 };
        logits[[row, 1]] = if row % 3 == 0 { -0.4 } else { 0.5 };
        for atom in 0..k_atoms {
            // Ordered Beta--Bernoulli gate: independent per-atom logistic, NOT a
            // normalized simplex weight.
            let gate = 1.0 / (1.0 + (-logits[[row, atom]] / tau).exp());
            let theta = std::f64::consts::TAU * coords[atom][[row, 0]];
            let basis = [1.0, theta.sin(), theta.cos()];
            for out_col in 0..p {
                for basis_col in 0..m {
                    target[[row, out_col]] +=
                        gate * basis[basis_col] * weights[atom][basis_col][out_col];
                }
            }
        }
        for out_col in 0..p {
            target[[row, out_col]] +=
                residual_scale * (((row * 7 + out_col * 3) as f64) * 0.7).sin();
        }
    }
    let mut atoms = Vec::with_capacity(k_atoms);
    for atom in 0..k_atoms {
        let (phi, jet) = evaluator.evaluate(coords[atom].view()).unwrap();
        let decoder = Array2::from_shape_fn((m, p), |(basis_col, out_col)| {
            weights[atom][basis_col][out_col]
        });
        atoms.push(
            SaeManifoldAtom::new_with_provided_function_gram(
                format!("patchd_{atom}"),
                SaeAtomBasisKind::Periodic,
                1,
                phi,
                jet,
                decoder,
                Array2::<f64>::eye(m),
            )
            .unwrap()
            .with_basis_second_jet(evaluator.clone()),
        );
    }
    let assignment = SaeAssignment::from_blocks_with_mode_and_manifolds(
        logits,
        coords,
        vec![LatentManifold::Circle { period: 1.0 }; k_atoms],
        AssignmentMode::ordered_beta_bernoulli(tau, 0.9, false),
    )
    .unwrap();
    let term = SaeManifoldTerm::new(atoms, assignment).unwrap();
    let rho = SaeManifoldRho::new(
        log_lambda_sparse,
        -6.0,
        vec![Array1::from_vec(vec![-6.0]), Array1::from_vec(vec![-6.0])],
    );
    (term, target, rho)
}

// #2330 Patch D prerequisite — map the residual scale at which the converged
// exact `A` stops being positive definite, and how big the residual-curvature
// block `ΔC` is inside that window. This decides whether the Patch-D FD arbiter
// can be anchored on a PD fixture at all, and separates "the shared fixture
// manufactured a saddle" from "every converged mode is an A-saddle" (the latter
// would gate #2330 behind #2336's saddle escape rather than behind Patch D).
#[test]
fn sae_exact_a_pd_window_scan_2330_patchd() {
    for &scale in &[0.0_f64, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.4] {
        let (mut term, target, rho) = obb_patchd_fixture(scale, -6.0);
        let built = term.penalized_quasi_laplace_criterion_with_cache(
            target.view(),
            &rho,
            None,
            200,
            0.4,
            1.0e-6,
            1.0e-6,
        );
        match built {
            Ok((_value, _loss, cache)) => {
                match term.exact_a_spectrum_summary(&rho, target.view(), &cache) {
                    Ok((min_eig, max_eig, n_neg, dc_frob, a_frob)) => {
                        eprintln!(
                            "PATCHD_WINDOW scale={scale:.4} PD_OK min_eig={min_eig:.6e} \
                             max_eig={max_eig:.6e} n_neg={n_neg} dc_frob={dc_frob:.6e} \
                             a_frob={a_frob:.6e} dc_rel={:.6e}",
                            dc_frob / a_frob.max(1.0e-300)
                        );
                        // This window scan decides where the Patch-D arbiter may be
                        // anchored, so each summary row must be a real spectrum
                        // summary: ordered finite extremes, Frobenius norms that are
                        // finite non-negative magnitudes, and a negative count that
                        // agrees with the reported minimum. A silently-NaN row would
                        // otherwise read as "PD_OK".
                        assert!(
                            min_eig.is_finite() && max_eig.is_finite() && min_eig <= max_eig,
                            "scale={scale}: the exact-A spectrum must be finite and ordered \
                             (min_eig={min_eig}, max_eig={max_eig})"
                        );
                        assert!(
                            dc_frob.is_finite()
                                && dc_frob >= 0.0
                                && a_frob.is_finite()
                                && a_frob >= 0.0,
                            "scale={scale}: Frobenius norms must be finite non-negative \
                             magnitudes (dc_frob={dc_frob}, a_frob={a_frob})"
                        );
                        // `n_neg` counts eigenvalues below the relative PD floor, so
                        // it may be 0 while `min_eig` is a hair negative — but it
                        // can never be positive unless the minimum is genuinely
                        // negative. That one-way implication is exact.
                        assert!(
                            n_neg == 0 || min_eig < 0.0,
                            "scale={scale}: {n_neg} eigenvalue(s) below the PD floor but the \
                             reported minimum is min_eig={min_eig} ≥ 0 — the count and the \
                             minimum are two readouts of one spectrum"
                        );
                    }
                    Err(e) => eprintln!("PATCHD_WINDOW scale={scale:.4} SPECTRUM_ERR {e}"),
                }
            }
            Err(e) => eprintln!("PATCHD_WINDOW scale={scale:.4} CRITERION_REFUSED {e:?}"),
        }
    }
}

// #2330 Patch D FD ARBITER — the per-coordinate gap between the analytic
// exact-A θ-adjoint `Γ_A,w = tr(A⁺ ∂A/∂θ_w)` and a CENTRAL DIFFERENCE of
// `exact_observed_information_log_dets(...).0 = log|A|` over frozen θ̂ with the
// cache REBUILT at each perturbed state (a frozen cache would false-green the
// gate). At baseline — before the Patch-D `∂ΔC/∂θ` legs land — the residual
// here IS the missing term, coordinate by coordinate.
//
// Anchored on `obb_patchd_fixture`, whose exact A is positive definite at the
// converged mode (see `sae_exact_a_pd_window_scan_2330_patchd`); the shared
// softmax fixture is not OBB-reachable and lands on an A-saddle where the
// criterion refuses outright.
#[test]
fn sae_exact_a_theta_adjoint_gap_measure_2330_patchd() {
    let (mut term, target, rho) = obb_patchd_fixture(0.0, -6.0);
    let (_value, _loss, cache) = term
        .penalized_quasi_laplace_criterion_with_cache(
            target.view(),
            &rho,
            None,
            200,
            0.4,
            1.0e-6,
            1.0e-6,
        )
        .expect("PD converged cache");
    let (log_a, log_a_tt) = term
        .exact_observed_information_log_dets(&rho, target.view(), &cache)
        .expect("exact-A log dets at the converged mode");
    eprintln!("PATCHD base log|A|={log_a:.9e} log|A_tt|={log_a_tt:.9e}");
    let gamma = term
        .exact_a_theta_adjoint_joint(&rho, target.view(), &cache)
        .expect("analytic exact-A joint theta adjoint");

    // Probe slots read off the ACTUAL cache layout rather than hardcoded, so a
    // layout change cannot silently repoint the probes at the wrong variables.
    let mut probes: Vec<(usize, usize, SaeLocalRowVar)> = Vec::new();
    for row in 0..3usize {
        let vars = term
            .row_vars_for_cache_row(row, &cache)
            .expect("row vars for probe layout");
        for (local, var) in vars.iter().enumerate() {
            probes.push((row, local, *var));
        }
    }
    probes.truncate(8);

    // This state has a narrow admitted basin: a 1e-4 coordinate perturbation
    // can leave it, and at 1e-5 central-difference truncation was almost 1%.
    // Verify convergence of the scalar oracle itself before comparing it with
    // the analytic derivative. Richardson removes its leading O(h²) term.
    const COORD_TOL: f64 = 1.0e-3;
    const LOGIT_TOL: f64 = 3.0e-2;
    let mut max_coord_rel = 0.0_f64;
    let mut max_logit_rel = 0.0_f64;
    for &(row, local, var) in &probes {
        // Logits move on the gate-temperature scale. Their response is O(1),
        // so tiny coordinate-sized steps amplify the scalar eigensolve's
        // rounding error. Coordinate responses here reach O(1e4), and need
        // the smaller step to remain inside the admitted basin.
        let coarse_step = match var {
            SaeLocalRowVar::Logit { .. } => match term.assignment.mode {
                AssignmentMode::OrderedBetaBernoulli { temperature, .. } => 1.0e-3 * temperature,
                _ => unreachable!("OBB fixture"),
            },
            SaeLocalRowVar::Coord { .. } => 1.0e-5,
        };
        let mut estimates = [0.0_f64; 3];
        for (step_index, divisor) in [1.0_f64, 2.0, 4.0].into_iter().enumerate() {
            let h = coarse_step / divisor;
            let mut plus = term.clone();
            let mut minus = term.clone();
            // Clone drops these lagged operands. Preserve the same frozen
            // Newton objective whose Hessian/third derivative is contracted.
            for endpoint in [&mut plus, &mut minus] {
                endpoint.decoder_repulsion_gate = term.decoder_repulsion_gate.clone();
                endpoint.barrier_coactivation_gate = term.barrier_coactivation_gate.clone();
                endpoint.amplitude_barrier_gate = term.amplitude_barrier_gate;
                endpoint.streaming_gates_frozen = true;
            }
            match var {
                SaeLocalRowVar::Logit { atom } => {
                    plus.assignment.logits[[row, atom]] += h;
                    minus.assignment.logits[[row, atom]] -= h;
                }
                SaeLocalRowVar::Coord { atom, axis } => {
                    let mut fp = plus.assignment.coords[atom].as_flat().clone();
                    let mut fm = minus.assignment.coords[atom].as_flat().clone();
                    let idx = row * plus.assignment.coords[atom].latent_dim() + axis;
                    fp[idx] += h;
                    fm[idx] -= h;
                    plus.assignment.coords[atom].set_flat(fp.view());
                    minus.assignment.coords[atom].set_flat(fm.view());
                }
            }
            let a = fixed_state_exact_a_logdet(plus, &target, &rho)
                .expect("every positive finite-difference endpoint must be admitted");
            let b = fixed_state_exact_a_logdet(minus, &target, &rho)
                .expect("every negative finite-difference endpoint must be admitted");
            estimates[step_index] = (a - b) / (2.0 * h);
        }
        let coarse = (4.0 * estimates[1] - estimates[0]) / 3.0;
        let fine = (4.0 * estimates[2] - estimates[1]) / 3.0;
        let tolerance = match var {
            SaeLocalRowVar::Coord { .. } => COORD_TOL,
            SaeLocalRowVar::Logit { .. } => LOGIT_TOL,
        };
        let oracle_error = (fine - coarse).abs() / (1.0 + fine.abs().max(coarse.abs()));
        assert!(
            oracle_error < tolerance,
            "row={row} var={var:?}: scalar FD oracle has not converged: estimates={estimates:?}, Richardson error={oracle_error:e}"
        );
        let analytic = gamma.t[cache.row_offsets[row] + local];
        let rel = (fine - analytic).abs() / (1.0 + fine.abs().max(analytic.abs()));
        // Charge the oracle's remaining uncertainty to the same error budget
        // instead of allowing it in addition to the derivative tolerance.
        let accounted_error = rel + oracle_error;
        match var {
            SaeLocalRowVar::Coord { .. } => max_coord_rel = max_coord_rel.max(accounted_error),
            SaeLocalRowVar::Logit { .. } => max_logit_rel = max_logit_rel.max(accounted_error),
        }
        eprintln!(
            "PATCHD_GAP row={row} local={local} var={var:?} fd={fine:.9e} analytic={analytic:.9e} rel={rel:.3e} oracle_error={oracle_error:.3e}"
        );
    }
    eprintln!("PATCHD_ARBITER max_coord_rel={max_coord_rel:.3e} max_logit_rel={max_logit_rel:.3e}");
    assert!(
        max_coord_rel < COORD_TOL,
        "exact-A theta-adjoint coordinate channel must match FD: max_coord_rel={max_coord_rel:.3e} >= {COORD_TOL:.1e}"
    );
    assert!(
        max_logit_rel < LOGIT_TOL,
        "exact-A theta-adjoint logit channel must match FD: max_logit_rel={max_logit_rel:.3e} >= {LOGIT_TOL:.1e}"
    );
}

// #2330 Patch D — channel-2 exercise gate. The main arbiter fixture sets
// log_lambda_sparse=-6 (OBB prior weight e^{-6}≈0.0025), so the ordered-BB prior
// curvature channel-2 (∂ΔC_obb/∂logit) is nearly inert there — correct-in-form
// but numerically ~0, which would let a channel-2 SIGN error ship silently.
// This variant raises the prior weight so channel-2 carries measurable weight;
// the logit slots staying FD-consistent here is what actually exercises its sign.
#[test]
fn sae_exact_a_theta_adjoint_gap_measure_2330_patchd_weighted() {
    // Scan a few sparse weights at residual_scale 0; report PD + the logit gaps so
    // a channel-2 sign error shows up as a blown logit slot.
    for &(rs, lls) in &[
        (0.005_f64, -4.0_f64),
        (0.005, -3.0),
        (0.01, -3.0),
        (0.02, -2.0),
    ] {
        let (mut term, target, rho) = obb_patchd_fixture(rs, lls);
        let built = term.penalized_quasi_laplace_criterion_with_cache(
            target.view(),
            &rho,
            None,
            200,
            0.4,
            1.0e-6,
            1.0e-6,
        );
        let cache = match built {
            Ok((_v, _l, c)) => c,
            Err(e) => {
                eprintln!("PATCHD_W lls={lls:.2} CRITERION_REFUSED {e:?}");
                continue;
            }
        };
        let gamma = match term.exact_a_theta_adjoint_joint(&rho, target.view(), &cache) {
            Ok(g) => g,
            Err(e) => {
                eprintln!("PATCHD_W lls={lls:.2} GAMMA_ERR {e}");
                continue;
            }
        };
        let h = 1.0e-5;
        for row in 0..2usize {
            let vars = term.row_vars_for_cache_row(row, &cache).expect("vars");
            for (local, var) in vars.iter().enumerate() {
                let mut plus = term.clone();
                let mut minus = term.clone();
                match *var {
                    SaeLocalRowVar::Logit { atom } => {
                        plus.assignment.logits[[row, atom]] += h;
                        minus.assignment.logits[[row, atom]] -= h;
                    }
                    SaeLocalRowVar::Coord { atom, axis } => {
                        let mut fp = plus.assignment.coords[atom].as_flat().clone();
                        let mut fm = minus.assignment.coords[atom].as_flat().clone();
                        let idx = row * plus.assignment.coords[atom].latent_dim() + axis;
                        fp[idx] += h;
                        fm[idx] -= h;
                        plus.assignment.coords[atom].set_flat(fp.view());
                        minus.assignment.coords[atom].set_flat(fm.view());
                    }
                }
                let analytic = gamma.t[cache.row_offsets[row] + local];
                match (
                    fixed_state_exact_a_logdet(plus, &target, &rho),
                    fixed_state_exact_a_logdet(minus, &target, &rho),
                ) {
                    (Some(a), Some(b)) => {
                        let fd = (a - b) / (2.0 * h);
                        let rel = (fd - analytic).abs() / (1.0 + fd.abs().max(analytic.abs()));
                        eprintln!(
                            "PATCHD_W rs={rs:.3} lls={lls:.2} row={row} var={var:?} fd={fd:.6e} \
                             analytic={analytic:.6e} rel={rel:.3e}"
                        );
                    }
                    _ => {
                        eprintln!("PATCHD_W rs={rs:.3} lls={lls:.2} row={row} var={var:?} refused")
                    }
                }
            }
        }
    }
}

// ─── #2080 / #2712 from-probes θ-adjoint parity, restored (#2818) ───────────
//
// `c0a21b554` deleted this gate because it no longer compiled, and it no longer
// compiled because `d484a091a` had deleted the FD-anchor scaffolding it stood
// on — `certified_fd_anchor`, `FdAnchorRegime`, `FdAnchorCandidate`,
// `rho_ladder_family`, `sparse_lift_ladder`, `deflation_blind_cache`. Every one
// of those was `#[cfg(test)]`, where the sweep's criterion ("no production
// artifact links this function") is true of everything by construction.
//
// The three production entry points the gate actually grades —
// `SaeManifoldTerm::logdet_theta_adjoint`,
// `SaeManifoldTerm::logdet_theta_adjoint_from_probes`, and
// `ArrowFactorCache::schur_inverse_apply` — were untouched, so this is a
// rebuild against them directly. Everything the anchor machinery did for the
// `any_maximum()` regime this gate declared is inlined as closures: walk the
// declared `log λ_sparse` ladder, converge each member's own inner mode, freeze
// it at `inner_max_iter = 0`, and accept the first member the criterion prices
// finitely. Nothing in that acceptance can see the finite difference or the
// analytic value, so it still cannot converge on "whatever agrees".

/// #2080 θ-adjoint from-probes — SOFTMAX fixture. Exercises the softmax entropy
/// dense off-diagonal channel + the core t–t / t–β / β–β selected-inverse folds.
///
/// The matrix-free θ-adjoint reconstructed from the FULL-BASIS probe bundle
/// (`z_j = √k·e_j`, exact dense `S⁻¹` via `cache.schur_inverse_apply`) must
/// reproduce the dense selected-inverse θ-adjoint. This isolates the from-probes
/// reconstruction; the dense adjoint is already FD-validated against `log|H|`
/// elsewhere.
///
/// On a DEFLATED cache the two have to be told apart from the deflation-blind
/// operator before agreement means anything: the deflated and undeflated
/// θ-adjoints coincide wherever the deflation is inactive, so machine-precision
/// agreement is ALSO what a port that ignored deflation would produce. The gate
/// measures `‖Γ_dense − Γ_deflation-blind‖∞` first and refuses to read parity as
/// evidence unless the two provably separate.
#[test]
fn sae_logdet_theta_adjoint_from_probes_matches_dense_softmax_2080() {
    // The declared `log λ_sparse` ladder. The assignment-strength penalty is the
    // dial that moves a state between the deflating and non-deflating regimes,
    // so it is the natural declared axis for a regime the gate needs but cannot
    // control directly. Ordered by lift; the accepted member is reported.
    const SPARSE_LIFTS: [f64; 9] = [0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5];

    let (base_term, target, base_rho) =
        crate::manifold::tests_recovery_split_780::gamma_fd_tiny_fixture();

    // The `any_maximum()` anchor regime, inlined: the frozen state must merely
    // BE a maximum the criterion will price. Everything the gate then
    // differentiates is defined there; nothing further is asserted about it.
    let mut rejections: Vec<String> = Vec::with_capacity(SPARSE_LIFTS.len());
    let mut certified: Option<(String, SaeManifoldTerm, SaeManifoldRho, ArrowFactorCache)> = None;
    for &lift in &SPARSE_LIFTS {
        let description = format!("log_lambda_sparse={lift:.2}");
        let mut rho = base_rho.clone();
        rho.log_lambda_sparse = lift;
        let mut term = base_term.clone();
        // Reach this member's own mode first. A solve that refuses is a
        // rejection of this member, not a panic.
        if let Err(error) = term.penalized_quasi_laplace_criterion_with_cache(
            target.view(),
            &rho,
            None,
            200,
            0.4,
            1.0e-6,
            1.0e-6,
        ) {
            rejections.push(format!("  {description}: inner solve refused: {error}"));
            continue;
        }
        // `inner_max_iter = 0` freezes θ̂ where the ladder put it: the anchor is
        // the point the gate declares, not wherever the solve would wander.
        match term.penalized_quasi_laplace_criterion_with_cache(
            target.view(),
            &rho,
            None,
            0,
            0.4,
            1.0e-6,
            1.0e-6,
        ) {
            Ok((value, loss, cache)) if value.is_finite() && loss.total().is_finite() => {
                certified = Some((description, term, rho, cache));
                break;
            }
            Ok((value, loss, _)) => rejections.push(format!(
                "  {description}: frozen state priced non-finitely (value={value}, loss={})",
                loss.total()
            )),
            Err(error) => rejections.push(format!(
                "  {description}: criterion refused the frozen state: {error}"
            )),
        }
    }
    let (accepted, term, rho, cache) = certified.unwrap_or_else(|| {
        panic!(
            "#2080 from-probes softmax parity: no member of the declared ladder is a maximum \
             the criterion will price. Widening the ladder is a fixture decision and dropping \
             the regime would change what is proved. Rejections:\n{}",
            rejections.join("\n")
        )
    });
    eprintln!("#2080 from-probes softmax parity: anchor certified at {accepted}");

    let solver = DeflatedArrowSolver::plain(&cache);
    let dense = term
        .logdet_theta_adjoint(&rho, &cache, &solver)
        .expect("dense theta-adjoint");

    let deflated_rows = cache
        .deflated_row_directions
        .iter()
        .filter(|d| !d.is_empty())
        .count();
    // The deflation-blind operator: the production dense adjoint against the
    // same cache with ONLY the deflation metadata stripped — the per-row
    // Cholesky factors and the reduced Schur are untouched. That is exactly what
    // a from-probes port which silently dropped the Daleckii–Krein correction
    // would return, so the distance to it is the resolution this gate has. It is
    // a REFERENCE, never a route.
    let separation = if deflated_rows == 0 {
        0.0
    } else {
        let mut blind = cache.clone();
        let rows = cache.deflated_row_directions.len();
        blind.deflated_row_directions = std::sync::Arc::from(vec![Vec::new(); rows]);
        blind.deflation_row_spectra = std::sync::Arc::from(vec![None; rows]);
        let blind_solver = DeflatedArrowSolver::plain(&blind);
        let blind_gamma = term
            .logdet_theta_adjoint(&rho, &blind, &blind_solver)
            .expect("deflation-blind dense theta-adjoint");
        dense
            .t
            .iter()
            .zip(blind_gamma.t.iter())
            .chain(dense.beta.iter().zip(blind_gamma.beta.iter()))
            .map(|(a, b)| (a - b).abs())
            .fold(0.0_f64, f64::max)
    };

    let k = cache.k;
    assert!(
        k > 0,
        "fixture must have a non-empty border to exercise S⁻¹ folds"
    );
    let sqrt_k = (k as f64).sqrt();
    let probes: Vec<Array1<f64>> = (0..k)
        .map(|j| {
            let mut v = Array1::<f64>::zeros(k);
            v[j] = sqrt_k;
            v
        })
        .collect();
    let sinv: Vec<Array1<f64>> = probes
        .iter()
        .map(|v| {
            cache
                .schur_inverse_apply(v.view())
                .expect("schur_inverse_apply")
        })
        .collect();
    let mf = term
        .logdet_theta_adjoint_from_probes(
            &rho,
            &cache,
            &probes,
            &sinv,
            EvidenceOperator::Majorizer,
            None,
        )
        .expect("matrix-free theta-adjoint");

    assert_eq!(dense.t.len(), mf.t.len());
    assert_eq!(dense.beta.len(), mf.beta.len());
    let mut max_abs = 0.0_f64;
    let mut parity = 0.0_f64;
    for (d, m) in dense
        .t
        .iter()
        .zip(mf.t.iter())
        .chain(dense.beta.iter().zip(mf.beta.iter()))
    {
        parity = parity.max((d - m).abs());
        max_abs = max_abs.max(d.abs());
    }
    eprintln!(
        "#2080/#2712 from-probes θ-adjoint gate: {deflated_rows} deflated row(s), \
         ‖Γ_dense‖∞ = {max_abs:.6e}, ‖Γ_dense − Γ_from-probes‖∞ = {parity:.6e}, \
         ‖Γ_dense − Γ_deflation-blind‖∞ = {separation:.6e}"
    );
    // #2712: `1e-8` per entry is the historical undeflated bar; on a deflated
    // cache it has to be finer than the correction it is supposed to be
    // sensitive to, which the assertion at the end checks against the measured
    // separation.
    let relative_tolerance = if deflated_rows > 0 { 1.0e-10 } else { 1.0e-8 };
    for (i, (d, m)) in dense.t.iter().zip(mf.t.iter()).enumerate() {
        assert!(
            (d - m).abs() <= relative_tolerance * (1.0 + d.abs()),
            "theta-adjoint gamma_t[{i}] mismatch: dense={d:.10e}, from_probes={m:.10e}"
        );
    }
    for (i, (d, m)) in dense.beta.iter().zip(mf.beta.iter()).enumerate() {
        assert!(
            (d - m).abs() <= relative_tolerance * (1.0 + d.abs()),
            "theta-adjoint gamma_beta[{i}] mismatch: dense={d:.10e}, from_probes={m:.10e}"
        );
    }
    assert!(
        max_abs > 0.0 && max_abs.is_finite(),
        "the theta-adjoint must be non-trivial to make the parity check meaningful"
    );
    if deflated_rows > 0 {
        // Non-vacuity, stated as a RATIO against the MEASURED separation rather
        // than as an absolute threshold. On the historical fixture the
        // correction moves Γ by 8.5e-8 against ‖Γ‖∞ = 98.9 — the deflated
        // direction is a near-null the raw derivative barely touches — so an
        // absolute floor would reject an honest fixture, while the per-entry
        // `1e-8·(1+|Γ|)` parity tolerance ALONE would admit a deflation-blind
        // port (8.5e-8 < 1e-6). The ratio is the margin by which such a port is
        // actually caught.
        assert!(
            separation > 0.0 && parity * 1.0e3 <= separation,
            "the deflated and deflation-blind θ-adjoints must SEPARATE before parity \
             is evidence of anything: a port that dropped the Daleckii–Krein \
             correction would also agree here. Measured separation {separation:.6e} \
             against parity error {parity:.6e} on {deflated_rows} deflated row(s)."
        );
        // ...and the parity tolerance the loops above applied must itself be
        // finer than the separation, or a deflation-blind port would slip
        // through them even though the ratio above holds.
        let loop_tolerance = relative_tolerance * (1.0 + max_abs);
        assert!(
            loop_tolerance < separation,
            "the per-entry parity tolerance {loop_tolerance:.6e} is coarser than the \
             {separation:.6e} distance to the deflation-blind operator, so the \
             element-wise assertions above would pass a port that dropped the \
             Daleckii–Krein correction. Tighten them or pick a fixture on which the \
             correction is larger."
        );
    }
}