gam-sae 0.3.155

Sparse-autoencoder latent-manifold terms for the gam penalized-likelihood engine
//! #2015 — column-equilibration primitive test, and the revived planted
//! ~1e4-spread convergence gate now driven by the SOLVER-LEVEL fix.
//!
//! HISTORY: unit-RMS DATA equilibration was wired into the crosscoder fit
//! path (`run_auto_sae_crosscoder_fit`/`run_sae_crosscoder_fit`) and REVERTED
//! the same night: for a homoscedastic reconstruction objective, dividing
//! columns by their RMS is not a reparametrization — it changes the estimand
//! (noise-dominated columns are amplified to unit RMS and the fit spends
//! capacity explaining them). Measured on MSI 13021686 at the wired commit:
//! `zz2015_tiny_inner_crawl_terminates` refused at the co-collapse floor
//! (EV 0.4566 vs null 0.4583) and the planted transport-law verdicts
//! collapsed (phase R² 0.139, smooth R² 0.837). The `equilibrate_crosscoder_columns`
//! primitive below is kept (spec'd + unit-tested); the fit path passes a unit
//! `column_scale` and no longer calls it.
//!
//! The κ≈1e8 conditioning fix instead lives in the inner SOLVER's linear
//! algebra (design: issue 2015 comment 4949898801, landed in
//! `gam_solve::arrow_schur::reduced_solve::factor_dense_reduced_schur`):
//! Jacobi/Van der Sluis diagonal equilibration of the dense reduced-Schur
//! complement BEFORE its Cholesky factorization, undone exactly (`L = D·L̃`)
//! before the factor is returned — a pure numerical-conditioning aid with no
//! signature change and no estimand change (unlike the reverted data-frame
//! attempt). Because that fix is internal to the solver, the test below
//! exercises the UNMODIFIED `run_auto_sae_behavior_fit` front door directly on
//! an ill-conditioned planted fixture — no equilibration call anywhere in
//! this test — and the fit should converge because the solver now factors the
//! ill-conditioned Schur accurately.

use ndarray::Array2;

use crate::manifold::{SaeBehaviorAutoFitRequest, SaeCrosscoderAutoFitConfig, run_auto_sae_behavior_fit};

/// Row-aligned circle activation (`p_x = 4`) and behavior probabilities
/// (`vocab = 4`, so `p_y = 3`) at `n` positions. `residual_amp` sets the
/// amplitude of ONE tangent direction (`e3`) relative to the fixed
/// `radius = 0.25` circular signal carried by the other two tangent
/// directions (`e1`, `e2`) — the column-scale-spread lever. The activation
/// construction does not depend on `residual_amp` at all, so the two fixtures
/// the test below compares share a BYTE-IDENTICAL activation target; only the
/// behavior block's own column-scale spread differs.
fn planted_fixture(n: usize, residual_amp: f64) -> (Array2<f64>, Array2<f64>) {
    let inv_sqrt_two = std::f64::consts::FRAC_1_SQRT_2;
    let b = [0.5_f64; 4];
    let e1 = [inv_sqrt_two, -inv_sqrt_two, 0.0, 0.0];
    let e2 = [0.5, 0.5, -0.5, -0.5];
    let e3 = [0.0, 0.0, inv_sqrt_two, -inv_sqrt_two];
    let radius = 0.25_f64;
    let mut z = Array2::<f64>::zeros((n, 4));
    let mut probs = Array2::<f64>::zeros((n, 4));
    for i in 0..n {
        let t = i as f64 / n as f64;
        let theta = std::f64::consts::TAU * t;
        // A small deterministic high-harmonic wiggle keeps the activation
        // residual positive (so the fixed-λ_y fit sees a genuine two-block
        // residual) without perturbing the circle the fitted decoder chases.
        let wiggle = 0.05 * (9.0 * theta).sin();
        z[[i, 0]] = theta.cos();
        z[[i, 1]] = theta.sin();
        z[[i, 2]] = 0.4 * (2.0 * theta).cos() + wiggle;
        z[[i, 3]] = 0.4 * (2.0 * theta).sin();

        // The behavior half-density: the circular signal in e1/e2 at fixed
        // `radius`, plus a residual in e3 whose amplitude is the spread lever.
        let residual = residual_amp * (3.0 * theta).sin();
        let radial = (1.0 - radius * radius - residual * residual).sqrt();
        let (cos_t, sin_t) = (theta.cos(), theta.sin());
        for token in 0..4 {
            let q = radial * b[token]
                + radius * (cos_t * e1[token] + sin_t * e2[token])
                + residual * e3[token];
            probs[[i, token]] = q * q;
        }
    }
    (z, probs)
}

fn fit_config() -> SaeCrosscoderAutoFitConfig {
    let mut config = SaeCrosscoderAutoFitConfig::standard(1, 3);
    config.max_iter = 40;
    // Fixed-rho (λ_y held at its seed value) isolates the INNER solve's
    // conditioning from the outer λ_y selection, matching the fast
    // `zz2015_tiny_inner_crawl_terminates` convergence gate this test mirrors.
    config.run_outer_rho_search = false;
    config
}

/// #2015 — a behavior block whose own tangent-coordinate columns carry a
/// genuine ~1e4 within-block scale spread must still converge through the
/// UNMODIFIED `run_auto_sae_behavior_fit` front door (no equilibration call
/// anywhere in this test — the fix lives entirely inside the solver's dense
/// Schur factorization), and must recover essentially the SAME joint
/// reconstruction as an otherwise-identical well-scaled (spread ~1) fixture —
/// proof that the solver-level Jacobi/Van der Sluis fix neutralizes the
/// spread's effect on the inner solve's conditioning without changing what
/// gets fitted (the estimand is untouched this time: nothing in the fit
/// pipeline's data frame changed, only the solver's internal factorization
/// path — see `factor_dense_reduced_schur_reconstructs_original_illconditioned_matrix_2015`
/// in `gam-solve` for the primitive-level identity this end-to-end gate
/// depends on).
#[test]
fn planted_1e4_column_spread_behavior_block_converges_and_matches_well_scaled() {
    let n = 96usize;
    // Well-scaled: the residual direction's amplitude is the same order as the
    // signal directions'.
    let (z_well, probs_well) = planted_fixture(n, 0.05);
    // Ill-conditioned: the residual direction's amplitude is 1e-4x the
    // well-scaled fixture's — a deliberate ~1e4 within-behavior-block
    // column-scale spread. The activation target is untouched (identical to
    // `z_well`), so any degradation in ITS reconstruction quality is entirely
    // a symptom of the shared joint solve's conditioning, not new activation
    // structure.
    let (z_ill, probs_ill) = planted_fixture(n, 0.05e-4);
    assert_eq!(
        z_well, z_ill,
        "the activation target must be untouched by the spread lever"
    );

    let well = run_auto_sae_behavior_fit(SaeBehaviorAutoFitRequest {
        activation: z_well,
        probabilities: probs_well,
        config: fit_config(),
        cancel: None,
    })
    .expect("well-scaled fixture must fit");

    let ill = run_auto_sae_behavior_fit(SaeBehaviorAutoFitRequest {
        activation: z_ill,
        probabilities: probs_ill,
        config: fit_config(),
        cancel: None,
    })
    .expect(
        "the solver's Jacobi-conditioned dense Schur factorization must let the ~1e4-spread \
         behavior block converge with NO data-frame equilibration in the fit path",
    );

    assert_eq!(ill.crosscoder.layers.len(), 2);
    for layer in &ill.crosscoder.layers {
        assert!(
            layer.reconstruction_r2.is_finite() && layer.reconstruction_r2 > 0.9,
            "{}: ill-conditioned fit must still reconstruct well, got {}",
            layer.label,
            layer.reconstruction_r2
        );
    }
    assert_eq!(
        ill.kl.infinite_rows, 0,
        "no fitted row may decode off-simplex"
    );

    // The recovered ACTIVATION reconstruction (gauge-invariant: it is decoded
    // output, not an internal coordinate) must match the well-scaled fixture's
    // closely — the shared latent/gate the behavior block's spread could have
    // destabilized did not, in fact, destabilize the activation block either.
    let activation_diff = &ill.crosscoder.layers[0].fitted - &well.crosscoder.layers[0].fitted;
    let activation_max_abs_diff = activation_diff
        .iter()
        .cloned()
        .fold(0.0_f64, |m, v| m.max(v.abs()));
    assert!(
        activation_max_abs_diff < 0.1,
        "activation reconstruction must match the well-scaled fixture's to tolerance, got max diff {activation_max_abs_diff}"
    );

    // The recovered BEHAVIOR (decoded distributions) must likewise match.
    let mut behavior_max_abs_diff = 0.0_f64;
    for (a, b) in ill
        .fitted_probabilities
        .iter()
        .zip(well.fitted_probabilities.iter())
    {
        behavior_max_abs_diff = behavior_max_abs_diff.max((a - b).abs());
    }
    assert!(
        behavior_max_abs_diff < 0.1,
        "fitted behavior probabilities must match the well-scaled fixture's to tolerance, got max diff {behavior_max_abs_diff}"
    );
}