feral 0.1.0

//! Global scaling for sparse symmetric indefinite matrices.
//!
//! Implements MC64-style matching-based scaling following
//! Duff & Koster 2001 and Duff & Pralet 2005, using a pure-Rust
//! Hungarian algorithm. The resulting scaling vector `s` is applied
//! symmetrically: `A ↦ diag(s) · A · diag(s)` before factorization.
//!
//! Design: see `dev/research/mc64-scaling.md`.
//! Plan:   see `dev/plans/mc64-scaling.md`.
//!
//! This module is Phase 2.2.1 work — closing the residual gap that
//! Phase 2.1.2's sanity check exposed on n > 500 matrices.
//!
//! ## Quick reference
//!
//! The caller computes scaling via `compute_scaling(matrix, strategy)`,
//! which returns `(Vec<f64>, ScalingInfo)`. The vector is in user-order
//! indexing (same numbering as the input CSC's row/column indices).
//! It is the responsibility of later symbolic-factorization code to
//! permute the vector into pivot-order before handing off to the
//! numeric phase.
//!
//! Once the scaling vector is available, three things must happen:
//!
//!   1. During frontal assembly in `numeric::factorize`, each original
//!      matrix entry `a[i,j]` is multiplied by `s[i] * s[j]` as it is
//!      scattered into the frontal matrix.
//!   2. In `numeric::solve`, the right-hand side `b` is pre-scaled by
//!      `b[i] *= s[i]` at the permutation boundary before the forward
//!      sweep.
//!   3. In `numeric::solve`, the solution `x` is post-scaled by
//!      `x[i] *= s[i]` at the un-permutation boundary after the
//!      backward sweep. **Same vector on both ends**, not its
//!      inverse — see the research note for the derivation.

use crate::error::FeralError;
use crate::sparse::csc::CscMatrix;

#[allow(dead_code)] // Real uses arrive in Step 3 of the implementation plan.
mod hungarian;
mod infnorm;
mod mc64;

/// Compute the MC64 symmetric matching on the lower-triangle CSC
/// matrix and return the column-to-row permutation (`perm[j]` is the
/// row matched to column `j`; `usize::MAX` marks unmatched columns).
///
/// Exposed for Phase 2.6.5 ordering-compression diagnostic work —
/// the matching cycle structure drives the MUMPS-style
/// `ICNTL(12)=2` quotient-graph compression. Internally this is
/// the same Hungarian call that `Mc64Symmetric` scaling uses, minus
/// the symmetric-average post-processing.
pub fn mc64_matching(matrix: &CscMatrix) -> Result<(Vec<usize>, usize), FeralError> {
    mc64::matching_perm(matrix)
}

/// Cached MC64 output: everything needed to both drive ordering
/// compression (`perm`) and derive the symmetric scaling vector
/// (`u`, `v`, `cmax`) without rerunning the expensive Hungarian
/// kernel. Produced by [`compute_mc64_cache`], consumed by
/// [`compute_scaling_with_cache`]. See Phase 2.4.4 compression
/// symbolic-speedup work.
pub(crate) use mc64::Mc64Cache;

/// Run the full MC64 pipeline once and return the cached output.
/// Used by the symbolic `LdltCompress` preprocessor so the numeric
/// phase can reuse the matching if its scaling strategy also
/// resolves to `Mc64Symmetric`.
pub(crate) fn compute_mc64_cache(matrix: &CscMatrix) -> Result<Mc64Cache, FeralError> {
    mc64::compute_matching(matrix)
}

/// User-facing scaling strategy selector.
///
/// Default is `Auto` — adaptive shape-based routing that picks
/// `Mc64Symmetric` for matrices with the arrow-KKT signature
/// (`diag_only / n >= 0.30`) and `InfNorm` everywhere else. Flipped
/// from the prior `InfNorm` default on 2026-04-19 after the
/// per-matrix residual-set diff confirmed the trade: 8× tail
/// compression on factor/MUMPS (worst case 83× → 10×) and material
/// wins on the VESUVIO/CRESC IPM corpus, against a net −9 change
/// in the residual_pass count out of 154 588. Of the 21 regressions,
/// 14 are oracle-`numerically_intractable` and 1 is `excluded`
/// (boundary flicker on already-hard matrices); 5 of the remaining
/// 6 `definitive` regressions are tolerance-edge effects (residuals
/// 1e-10 → 1e-9 around the `n·ε·1e6` threshold). The lone material
/// residual regression is MSS1_0009 (6e-12 → 1e-6, inertia preserved).
/// Inertia hard rule is satisfied on every regression. See
/// `dev/research/lever-c-residual-diff-2026-04-19.md`.
///
/// `InfNorm` (Knight-Ruiz iterative ∞-norm equilibration) is still
/// available as an opt-in; it is the only choice that solves
/// MSS1_0009 to working precision today and is the right pick for
/// pipelines that cannot tolerate the MSS1-class residual loss
/// pending Policy 4 (post-scaling trial-residual diagnostic).
///
/// `Mc64Symmetric` is also opt-in; it is useful on matrices where
/// matching provides better conditioning than ∞-norm balancing
/// (e.g. SSINE_2529, VESUVIA_0000 in the parity panel) but pays the
/// MC64 symbolic overhead unconditionally.
#[derive(Debug, Clone, PartialEq, Default)]
pub enum ScalingStrategy {
    /// Knight-Ruiz ∞-norm iterative equilibration. Matches the
    /// scaling algorithm used by the dense BK path. Was the default
    /// from Phase 2.2.3 through the 2026-04-19 lever-C residual diff
    /// (now opt-in).
    InfNorm,
    /// MC64-style symmetric matching-based scaling. Matches the
    /// default behavior of MUMPS (SYM=2) and SSIDS
    /// (options%scaling=1). Useful on matrices where matching
    /// provides better conditioning than ∞-norm balancing.
    Mc64Symmetric,
    /// Identity scaling (no-op). Use for regression testing and for
    /// inputs where any scaling is inappropriate.
    Identity,
    /// User-supplied pre-computed scaling vector in user-order
    /// indexing. Length must equal the matrix dimension.
    External(Vec<f64>),
    /// Adaptive shape-based routing: `Mc64Symmetric` when the matrix
    /// has the arrow-KKT signature (many degree-1 "constraint slack"
    /// columns), else `InfNorm`. The routing rule is documented at
    /// [`pick_scaling_strategy`]; threshold is `diag_only / n >= 0.3`.
    /// Default since 2026-04-19. See
    /// `dev/research/lever-c-residual-diff-2026-04-19.md`.
    #[default]
    Auto,
}

/// Diagnostic information about how the scaling was computed.
#[derive(Debug, Clone, PartialEq)]
pub enum ScalingInfo {
    /// MC64 matching ran to completion on a non-singular matrix.
    Applied,
    /// MC64 matching found a partial solution; unmatched rows and
    /// columns fall back to identity scaling. `n_unmatched` is the
    /// number of variables that could not be matched. The returned
    /// scaling vector has `1.0` at the unmatched positions.
    PartialSingular { n_unmatched: usize },
    /// No matching-based scaling was applied (e.g., the caller
    /// requested `Identity` or `External`).
    NotApplied,
}

/// Compute the symmetric scaling vector for a sparse symmetric
/// matrix stored in CSC with only the lower triangle, following
/// `strategy`.
///
/// Returns a vector of length `n` in **user-order** indexing such
/// that applying `D = diag(scaling)` as the congruence transform
/// `D · A · D` produces a matrix whose largest-magnitude entries lie
/// on the diagonal. The off-diagonals are bounded by 1 in absolute
/// value when MC64 succeeds on a non-singular matrix.
///
/// Users of the result must permute the vector into pivot-order
/// indexing before the numeric phase looks it up.
pub fn compute_scaling(
    matrix: &CscMatrix,
    strategy: &ScalingStrategy,
) -> Result<(Vec<f64>, ScalingInfo), FeralError> {
    compute_scaling_with_cache(matrix, strategy, None)
}

/// Variant of [`compute_scaling`] that accepts a precomputed MC64
/// cache. When the strategy resolves to `Mc64Symmetric` (including
/// via `Auto` routing), the cache is consumed in O(n) — no Hungarian
/// rerun. When the strategy does not end up running MC64 (Identity,
/// External, InfNorm, or Auto resolving to InfNorm with Policy 4
/// fallback), the cache is ignored and the regular path runs.
///
/// `cache` must be `compute_mc64_cache(matrix)` on the same matrix,
/// else the produced scaling is wrong.
pub(crate) fn compute_scaling_with_cache(
    matrix: &CscMatrix,
    strategy: &ScalingStrategy,
    cache: Option<&Mc64Cache>,
) -> Result<(Vec<f64>, ScalingInfo), FeralError> {
    match strategy {
        ScalingStrategy::Identity => Ok((vec![1.0; matrix.n], ScalingInfo::NotApplied)),
        ScalingStrategy::External(s) => {
            if s.len() != matrix.n {
                return Err(FeralError::InvalidInput(format!(
                    "external scaling has length {} but matrix has n={}",
                    s.len(),
                    matrix.n,
                )));
            }
            Ok((s.clone(), ScalingInfo::NotApplied))
        }
        ScalingStrategy::InfNorm => Ok(infnorm::compute_infnorm(matrix)),
        ScalingStrategy::Mc64Symmetric => match cache {
            Some(c) => Ok(mc64::scaling_from_cache(c)),
            None => mc64::compute_symmetric(matrix),
        },
        ScalingStrategy::Auto => compute_scaling_auto_with_cache(matrix, cache),
    }
}

/// Resolve `ScalingStrategy::Auto` with a Policy 4 fallback rule:
/// when `pick_scaling_strategy` would pick `Mc64Symmetric`, check
/// whether MC64 has produced a scaling that is catastrophically
/// worse than InfNorm on a matrix where InfNorm would have done
/// fine. If so, fall back to InfNorm.
///
/// Rule (all three must fire):
/// 1. `raw_diag_range < RAW_GUARD` — the raw matrix's diagonal
///    spans only a few orders of magnitude. MC64 has nothing
///    to recover from raw ill-conditioning here, so any huge
///    scaled off/diag ratio it produces is pure artifact, not
///    reflection of inherent matrix difficulty.
/// 2. `mc_off > MC_OFF_GUARD` — MC64's scaled `max(|off|/|diag|)`
///    is large in absolute terms.
/// 3. `mc_off / in_off > RATIO_GUARD` — and is much larger
///    than what InfNorm produces.
///
/// The first guard is the critical one: it lets matrices like
/// MEYER3NE_0220 (raw_drng=4.77e19, but MC64 actually works) keep
/// MC64, while still catching MSS1_0009 (raw_drng=51, where MC64
/// produces noise).
///
/// Validated on a 17-matrix panel: MSS1_0009 falls back (recovers
/// the 6e-12 InfNorm residual instead of the 1e-6 MC64 residual);
/// VESUVIA / VESUVIO / VESUVIOU / MUONSINE / CRESC132 / HS75 /
/// MEYER3NE all keep MC64 (preserving the 84× → 9.4× factor
/// speedup, the 4-order HS75 residual win, and the MEYER3NE parity
/// tests). See `dev/research/policy-4-scaling-fallback.md`.
fn compute_scaling_auto_with_cache(
    matrix: &CscMatrix,
    cache: Option<&Mc64Cache>,
) -> Result<(Vec<f64>, ScalingInfo), FeralError> {
    const RAW_GUARD: f64 = 1e6;
    const MC_OFF_GUARD: f64 = 1e6;
    const RATIO_GUARD: f64 = 1e5;

    let picked = pick_scaling_strategy(matrix);
    if !matches!(picked, ScalingStrategy::Mc64Symmetric) {
        // Auto picked InfNorm-class — no fallback needed. Cache is
        // unused; MC64 was speculative work for compression and has
        // no payoff on this branch.
        return compute_scaling(matrix, &picked);
    }

    let mc64_from_cache = |matrix: &CscMatrix| -> Result<(Vec<f64>, ScalingInfo), FeralError> {
        match cache {
            Some(c) => Ok(mc64::scaling_from_cache(c)),
            None => mc64::compute_symmetric(matrix),
        }
    };

    // Cheap pre-filter: a wide raw |diag| range means MC64 has
    // genuine work to do. Skip the diagnostic and use MC64.
    if raw_diag_range(matrix) >= RAW_GUARD {
        return mc64_from_cache(matrix);
    }

    let (mc_vec, mc_info) = mc64_from_cache(matrix)?;
    let mc_off = max_off_diag_ratio(matrix, &mc_vec);
    if mc_off <= MC_OFF_GUARD {
        // MC64 produced a well-conditioned scaled matrix.
        return Ok((mc_vec, mc_info));
    }
    let (in_vec, in_info) = infnorm::compute_infnorm(matrix);
    let in_off = max_off_diag_ratio(matrix, &in_vec);
    let ratio = if in_off > 0.0 {
        mc_off / in_off
    } else {
        f64::INFINITY
    };
    if ratio > RATIO_GUARD {
        // MC64 is catastrophically worse than InfNorm AND the raw
        // matrix is already well-behaved — fall back to InfNorm.
        // Forward `in_info` (`ScalingInfo::Applied`) so the solve
        // path actually uses the scaling vector instead of treating
        // it as identity.
        Ok((in_vec, in_info))
    } else {
        Ok((mc_vec, mc_info))
    }
}

/// Compute `max |A_{j,j}| / min(|A_{j,j}|)` over diagonal entries
/// that are present and nonzero. Returns `+∞` if no nonzero
/// diagonal is present. O(nnz), no allocations.
fn raw_diag_range(matrix: &CscMatrix) -> f64 {
    let n = matrix.n;
    if n == 0 {
        return 0.0;
    }
    let mut lo = f64::INFINITY;
    let mut hi = 0.0_f64;
    for j in 0..n {
        for k in matrix.col_ptr[j]..matrix.col_ptr[j + 1] {
            if matrix.row_idx[k] == j {
                let a = matrix.values[k].abs();
                if a > 0.0 {
                    if a < lo {
                        lo = a;
                    }
                    if a > hi {
                        hi = a;
                    }
                }
            }
        }
    }
    if lo.is_finite() && lo > 0.0 {
        hi / lo
    } else {
        f64::INFINITY
    }
}

/// Compute `max_j (max_{i ≠ j} |s_i · A_{i,j} · s_j|) / |s_j · A_{j,j} · s_j|`
/// over all columns of the symmetrically-scaled matrix `D · A · D`.
/// Diagonal columns with zero diagonal contribute `+∞` to the max.
/// O(nnz), no allocations.
fn max_off_diag_ratio(matrix: &CscMatrix, scaling: &[f64]) -> f64 {
    let n = matrix.n;
    if n == 0 {
        return 0.0;
    }
    let mut diag_abs = vec![0.0_f64; n];
    let mut max_off = vec![0.0_f64; n];
    for j in 0..n {
        for k in matrix.col_ptr[j]..matrix.col_ptr[j + 1] {
            let i = matrix.row_idx[k];
            let v = (matrix.values[k] * scaling[i] * scaling[j]).abs();
            if i == j {
                diag_abs[j] = v;
            } else {
                if v > max_off[i] {
                    max_off[i] = v;
                }
                if v > max_off[j] {
                    max_off[j] = v;
                }
            }
        }
    }
    let mut worst = 0.0_f64;
    for j in 0..n {
        let r = if diag_abs[j] > 0.0 {
            max_off[j] / diag_abs[j]
        } else if max_off[j] > 0.0 {
            f64::INFINITY
        } else {
            0.0
        };
        if r > worst {
            worst = r;
        }
    }
    worst
}

/// Resolve `ScalingStrategy::Auto` to a concrete strategy based on
/// matrix shape.
///
/// Routes to `Mc64Symmetric` when the matrix has the arrow-KKT
/// signature — many degree-1 "constraint slack" columns whose only
/// stored row is the diagonal. Else routes to `InfNorm`.
///
/// Threshold: `diag_only / n >= 0.3`. Selected from the `vesuvio_diag`
/// shape distribution: VESUVIOU/VESUVIO/VESUVIA/MUONSINE/CRESC132 all
/// have ratios above 0.3 and benefit from MC64 (delays drop to zero,
/// 6×–229× factor speedup); HYDCAR20/METHANL8/SWOPF/HATFLDG (the
/// matrices that motivated the InfNorm default) have ratios below
/// 0.3. See `dev/research/lever-c-adaptive-scaling.md`.
///
/// One O(n) pass over the column pointers and one O(nnz) pass over
/// the row indices. No allocations.
pub fn pick_scaling_strategy(matrix: &CscMatrix) -> ScalingStrategy {
    let n = matrix.n;
    if n == 0 {
        return ScalingStrategy::InfNorm;
    }
    let mut diag_only = 0usize;
    for j in 0..n {
        let start = matrix.col_ptr[j];
        let end = matrix.col_ptr[j + 1];
        if end - start != 1 {
            continue;
        }
        if matrix.row_idx[start] == j {
            diag_only += 1;
        }
    }
    if diag_only as f64 / n as f64 >= 0.3 {
        ScalingStrategy::Mc64Symmetric
    } else {
        ScalingStrategy::InfNorm
    }
}

// Hungarian types are used by the `mc64` module once Step 3 lands.
// Not part of the public API.
#[allow(unused_imports)]
pub(crate) use hungarian::{hungarian_match, CostGraph, Matching};

#[cfg(test)]
mod tests {
    use super::*;
    use crate::sparse::csc::CscMatrix;

    /// Build a CSC with `n` columns where the first `diag_only`
    /// columns are degree-1 (just the diagonal), and the remaining
    /// `n - diag_only` columns each store the diagonal plus one
    /// off-diagonal row at column 0. Lower-triangular only — no
    /// validity beyond the column-degree pattern is required for
    /// `pick_scaling_strategy`, which only inspects col_ptr and
    /// row_idx.
    fn shape_csc(n: usize, diag_only: usize) -> CscMatrix {
        assert!(diag_only <= n);
        let mut col_ptr = Vec::with_capacity(n + 1);
        let mut row_idx: Vec<usize> = Vec::new();
        let mut values: Vec<f64> = Vec::new();
        col_ptr.push(0);
        for j in 0..n {
            row_idx.push(j);
            values.push(1.0);
            if j >= diag_only && j != 0 {
                row_idx.push(j.max(1) - 1);
                values.push(0.1);
            }
            col_ptr.push(row_idx.len());
        }
        CscMatrix {
            n,
            col_ptr,
            row_idx,
            values,
        }
    }

    #[test]
    fn pick_scaling_strategy_picks_mc64_for_arrow_kkt() {
        // 10 of 20 columns are diag-only → ratio = 0.5 ≥ 0.3.
        let csc = shape_csc(20, 10);
        assert_eq!(pick_scaling_strategy(&csc), ScalingStrategy::Mc64Symmetric);
    }

    #[test]
    fn pick_scaling_strategy_picks_infnorm_for_dense() {
        // 0 of 20 columns are diag-only → ratio = 0.0 < 0.3.
        let csc = shape_csc(20, 0);
        assert_eq!(pick_scaling_strategy(&csc), ScalingStrategy::InfNorm);
    }

    #[test]
    fn pick_scaling_strategy_threshold_boundary() {
        // 29 of 100 → 0.29 < 0.30 → InfNorm.
        let below = shape_csc(100, 29);
        assert_eq!(pick_scaling_strategy(&below), ScalingStrategy::InfNorm);
        // 30 of 100 → 0.30 ≥ 0.30 → MC64.
        let at = shape_csc(100, 30);
        assert_eq!(pick_scaling_strategy(&at), ScalingStrategy::Mc64Symmetric);
    }

    #[test]
    fn pick_scaling_strategy_empty_matrix_picks_infnorm() {
        let csc = CscMatrix {
            n: 0,
            col_ptr: vec![0],
            row_idx: vec![],
            values: vec![],
        };
        assert_eq!(pick_scaling_strategy(&csc), ScalingStrategy::InfNorm);
    }

    #[test]
    fn compute_scaling_auto_routes_to_mc64_on_arrow_kkt() {
        // Build a small symmetric arrow KKT: 4 diag-only "slack"
        // columns + 2 dense "linking" columns. Lower-triangular CSC.
        // Ratio diag_only / n = 4/6 = 0.67 → Auto resolves to MC64.
        let n = 6;
        let mut col_ptr = vec![0usize];
        let mut row_idx = Vec::new();
        let mut values = Vec::new();
        // 4 diag-only columns.
        for j in 0..4 {
            row_idx.push(j);
            values.push(2.0);
            col_ptr.push(row_idx.len());
        }
        // 2 dense columns (diagonal + all earlier rows).
        for j in 4..n {
            row_idx.push(j);
            values.push(2.0);
            for i in (j + 1)..n {
                row_idx.push(i);
                values.push(0.1);
            }
            col_ptr.push(row_idx.len());
        }
        let csc = CscMatrix {
            n,
            col_ptr,
            row_idx,
            values,
        };
        assert_eq!(pick_scaling_strategy(&csc), ScalingStrategy::Mc64Symmetric);
        // Auto and explicit Mc64Symmetric must produce the same vector
        // here — this is a well-conditioned shape, so the Policy 4
        // fallback rule (mc_off > 1e6 ∧ mc_off/in_off > 1e5) never fires.
        let (auto_s, _) =
            compute_scaling(&csc, &ScalingStrategy::Auto).expect("Auto routing should succeed");
        let (mc64_s, _) =
            compute_scaling(&csc, &ScalingStrategy::Mc64Symmetric).expect("MC64 should succeed");
        assert_eq!(auto_s, mc64_s);
    }

    #[test]
    fn max_off_diag_ratio_basic_well_conditioned() {
        // 3x3 well-conditioned matrix:
        //   [ 4  1  0 ]
        //   [ 1  3  1 ]
        //   [ 0  1  2 ]
        // With identity scaling, max ratio = 1/2 = 0.5.
        let csc = CscMatrix {
            n: 3,
            col_ptr: vec![0, 2, 4, 5],
            row_idx: vec![0, 1, 1, 2, 2],
            values: vec![4.0, 1.0, 3.0, 1.0, 2.0],
        };
        let s = vec![1.0; 3];
        let r = max_off_diag_ratio(&csc, &s);
        assert!((r - 0.5).abs() < 1e-12, "got {r}");
    }

    #[test]
    fn max_off_diag_ratio_zero_diag_gives_infinity() {
        // 2x2 with zero diagonal on column 0:
        //   [ 0  1 ]
        //   [ 1  1 ]
        // Column 0 has off=1, diag=0 → +inf. Column 1 has off=1,
        // diag=1 → 1.0. max = +inf.
        let csc = CscMatrix {
            n: 2,
            col_ptr: vec![0, 2, 3],
            row_idx: vec![0, 1, 1],
            values: vec![0.0, 1.0, 1.0],
        };
        let s = vec![1.0; 2];
        let r = max_off_diag_ratio(&csc, &s);
        assert!(r.is_infinite(), "got {r}");
    }

    /// Policy 4 fallback regression test — MSS1_0009 should resolve
    /// to InfNorm under Auto despite the diag_only/n=0.45 ratio
    /// triggering the MC64 routing rule. The fallback fires because
    /// MC64 produces a scaled `max(|off|/|diag|) ≈ 7.8e14` while
    /// InfNorm gets ≈ 2.0e8 — ratio 3.9e6 is well above the
    /// 1e5 RATIO_GUARD. See `dev/research/policy-4-scaling-fallback.md`
    /// table for the full numbers.
    #[test]
    fn auto_falls_back_to_infnorm_on_mss1_0009() {
        let path = std::path::Path::new("data/matrices/kkt/MSS1/MSS1_0009.mtx");
        let mtx = match crate::io::mtx::read_mtx(path) {
            Ok(m) => m,
            Err(_) => return, // fixture not present — skip
        };
        let csc = mtx.to_csc().expect("MSS1_0009 CSC build");

        // pick_scaling_strategy still picks MC64 — the routing rule
        // hasn't changed.
        assert_eq!(pick_scaling_strategy(&csc), ScalingStrategy::Mc64Symmetric);

        // But Auto should resolve to the InfNorm scaling because of
        // the Policy 4 fallback.
        let (auto_s, _) = compute_scaling(&csc, &ScalingStrategy::Auto)
            .expect("Auto on MSS1_0009 should succeed");
        let (in_s, _) = compute_scaling(&csc, &ScalingStrategy::InfNorm)
            .expect("InfNorm on MSS1_0009 should succeed");
        let (mc_s, _) = compute_scaling(&csc, &ScalingStrategy::Mc64Symmetric)
            .expect("MC64 on MSS1_0009 should succeed");
        assert_eq!(auto_s, in_s, "Auto must fall back to InfNorm on MSS1_0009");
        assert_ne!(
            auto_s, mc_s,
            "Auto must NOT use MC64 on MSS1_0009 (would regress residual to 1e-6)"
        );
    }

    /// Policy 4 fallback must NOT fire on the VESUVIO/CRESC class —
    /// these are the matrices the lever-C win is built on. MC64
    /// produces a scaled `mc_off ≈ 4.84e12` for VESUVIA_0000 with
    /// `mc/in ≈ 40` — well below the 1e5 RATIO_GUARD.
    #[test]
    fn auto_keeps_mc64_on_vesuvia_0000() {
        let path = std::path::Path::new("data/matrices/kkt/VESUVIA/VESUVIA_0000.mtx");
        let mtx = match crate::io::mtx::read_mtx(path) {
            Ok(m) => m,
            Err(_) => return, // fixture not present — skip
        };
        let csc = mtx.to_csc().expect("VESUVIA_0000 CSC build");
        let (auto_s, _) =
            compute_scaling(&csc, &ScalingStrategy::Auto).expect("Auto on VESUVIA_0000");
        let (mc_s, _) =
            compute_scaling(&csc, &ScalingStrategy::Mc64Symmetric).expect("MC64 on VESUVIA_0000");
        assert_eq!(auto_s, mc_s, "Auto must keep MC64 on VESUVIA_0000");
    }

    /// Same shape as `auto_keeps_mc64_on_vesuvia_0000` for the
    /// VESUVIOU subfamily — the highest mc/in ratio in the
    /// validation panel (1.05e4) is on this matrix; the threshold
    /// has 10× margin.
    #[test]
    fn auto_keeps_mc64_on_vesuviou_0000() {
        let path = std::path::Path::new("data/matrices/kkt/VESUVIOU/VESUVIOU_0000.mtx");
        let mtx = match crate::io::mtx::read_mtx(path) {
            Ok(m) => m,
            Err(_) => return, // fixture not present — skip
        };
        let csc = mtx.to_csc().expect("VESUVIOU_0000 CSC build");
        let (auto_s, _) =
            compute_scaling(&csc, &ScalingStrategy::Auto).expect("Auto on VESUVIOU_0000");
        let (mc_s, _) =
            compute_scaling(&csc, &ScalingStrategy::Mc64Symmetric).expect("MC64 on VESUVIOU_0000");
        assert_eq!(auto_s, mc_s, "Auto must keep MC64 on VESUVIOU_0000");
    }

    /// HS75_0000 is one of the two "material wins" — MC64 scaling
    /// gives a 4-order residual improvement. The fallback rule
    /// must not fire here. mc_off=1.08e9, in_off=1e10, ratio=0.108
    /// — well below 1.0, far from the 1e5 RATIO_GUARD.
    #[test]
    fn auto_keeps_mc64_on_hs75_0000() {
        let path = std::path::Path::new("data/matrices/kkt/HS75/HS75_0000.mtx");
        let mtx = match crate::io::mtx::read_mtx(path) {
            Ok(m) => m,
            Err(_) => return, // fixture not present — skip
        };
        let csc = mtx.to_csc().expect("HS75_0000 CSC build");
        let (auto_s, _) = compute_scaling(&csc, &ScalingStrategy::Auto).expect("Auto on HS75_0000");
        let (mc_s, _) =
            compute_scaling(&csc, &ScalingStrategy::Mc64Symmetric).expect("MC64 on HS75_0000");
        assert_eq!(auto_s, mc_s, "Auto must keep MC64 on HS75_0000 (the win)");
    }
}