Entropyk/crates/solver/tests/jacobian_scaling.rs

//! CM1.5 — acceptance tests for Jacobian row/column equilibration (NFR1).
//!
//! These tests prove the equilibration requirement on a *multi-circuit,
//! mixed-unit* Jacobian (the kind produced by a two-circuit `(ṁ, P, h)` system,
//! where ṁ ≈ 1 kg/s, P ≈ 1e6 Pa, h ≈ 3e5 J/kg):
//!
//! 1. The condition number drops by ≥ 1e4 versus the unscaled matrix.
//! 2. The equilibrated solve returns the same Newton step as an unscaled
//!    reference solve, within tight relative tolerance (solution-preserving).
//!
//! A faithful synthetic stand-in is used so the test is deterministic and free
//! of any fluid-backend dependency: a well-conditioned base matrix `W` is framed
//! by physical magnitudes via `J = diag(mag) · W · diag(mag)`. This reproduces
//! exactly the ill-scaling that wrecks conditioning in the real assembled
//! Jacobian, while keeping the *intrinsic* problem (`W`) benign — so any κ blow-up
//! is purely a scaling artifact that equilibration must remove.

use entropyk_solver::{equilibrate, JacobianMatrix};
use nalgebra::{DMatrix, DVector};

/// Well-conditioned, diagonally-dominant base matrix for a 2-circuit layout.
///
/// Indices 0,1,2 = (ṁ, P, h) of circuit A; 3,4,5 = circuit B. The (0,5), (2,3),
/// (3,2), (5,0) entries model weak inter-circuit (thermal) coupling, so the
/// matrix is NOT block-diagonal — a realistic coupled system.
fn base_matrix() -> DMatrix<f64> {
    DMatrix::from_row_slice(
        6,
        6,
        &[
            2.0, 0.4, 0.1, 0.0, 0.0, 0.05, //
            0.3, 2.0, 0.5, 0.0, 0.0, 0.0, //
            0.1, 0.3, 2.0, 0.05, 0.0, 0.0, //
            0.0, 0.0, 0.05, 2.0, 0.4, 0.1, //
            0.0, 0.0, 0.0, 0.3, 2.0, 0.5, //
            0.05, 0.0, 0.0, 0.1, 0.3, 2.0, //
        ],
    )
}

/// Builds `J = diag(mag) · W · diag(mag)` and returns it as a `JacobianMatrix`.
fn scaled_system(mag: &[f64]) -> (DMatrix<f64>, JacobianMatrix) {
    let w = base_matrix();
    let n = w.nrows();
    let mut entries = Vec::with_capacity(n * n);
    let mut dense = DMatrix::zeros(n, n);
    for i in 0..n {
        for j in 0..n {
            let v = mag[i] * w[(i, j)] * mag[j];
            dense[(i, j)] = v;
            entries.push((i, j, v));
        }
    }
    (dense.clone(), JacobianMatrix::from_builder(&entries, n, n))
}

/// κ via SVD (σ_max / σ_min), skipping exact-zero singular values.
fn condition_number(m: &DMatrix<f64>) -> f64 {
    let svd = m.clone().svd(false, false);
    let sv = svd.singular_values;
    let sigma_max = sv.max();
    let sigma_min = sv
        .iter()
        .filter(|&&s| s > 0.0)
        .cloned()
        .fold(f64::INFINITY, f64::min);
    sigma_max / sigma_min
}

/// AC #3 (bullet 1+2): on a realistic mixed-unit (Pa + J/kg + kg/s) two-circuit
/// Jacobian, equilibration slashes the condition number by ≥ 1e4.
#[test]
fn test_equilibration_reduces_condition_number_realistic_magnitudes() {
    // ṁ ≈ 1, P ≈ 1e6 Pa, h ≈ 3e5 J/kg, repeated for two circuits.
    let mag = [1.0, 1.0e6, 3.0e5, 1.0, 1.0e6, 3.0e5];
    let (dense, _jac) = scaled_system(&mag);

    let cond_before = condition_number(&dense);
    // Sanity: the raw problem really is badly conditioned.
    assert!(
        cond_before > 1.0e8,
        "raw κ should be large for mixed units, got {:.3e}",
        cond_before
    );

    let (d_r, d_c) = equilibrate(&dense);
    let mut scaled = dense.clone();
    for i in 0..6 {
        for j in 0..6 {
            scaled[(i, j)] *= d_r[i] * d_c[j];
        }
    }
    let cond_after = condition_number(&scaled);

    assert!(
        cond_after <= cond_before / 1.0e4,
        "equilibration must cut κ by ≥1e4: before={:.3e}, after={:.3e} (ratio {:.3e})",
        cond_before,
        cond_after,
        cond_before / cond_after
    );
}

/// AC #3 (bullet 3) + AC #4: the equilibrated `JacobianMatrix::solve` returns the
/// same Newton step as an unscaled reference LU solve, within 1e-9 relative — the
/// scaling is solution-preserving. Uses a mixed-unit system whose conditioning
/// (κ ≈ 1e6) is still comfortably resolvable in f64, so the 1e-9 comparison is
/// meaningful while κ reduction (≥1e4) still holds.
#[test]
fn test_equilibrated_solve_matches_unscaled_reference() {
    // Mixed scales spanning 1e3 (kg/s vs reduced-pressure scale): κ_raw ≈ 1e6.
    let mag = [1.0, 1.0e3, 3.0e2, 1.0, 1.0e3, 3.0e2];
    let (dense, jac) = scaled_system(&mag);

    // Known step we want to recover.
    let x_true = DVector::from_row_slice(&[0.7, -1.3, 2.1, -0.4, 0.9, -1.1]);
    // b = J · x_true; we want J · Δx = b, i.e. solve() with r = -b → Δx = x_true.
    let b = &dense * &x_true;
    let r: Vec<f64> = b.iter().map(|v| -v).collect();

    // Equilibrated solve (the production path).
    let delta = jac.solve(&r).expect("non-singular");

    // Unscaled reference: direct LU on the raw matrix.
    let dx_ref = dense.clone().lu().solve(&b).expect("reference LU solves");

    for k in 0..6 {
        let scale = x_true[k].abs().max(1.0);
        assert!(
            (delta[k] - x_true[k]).abs() / scale < 1e-9,
            "equilibrated step differs from x_true at {}: got {}, want {}",
            k,
            delta[k],
            x_true[k]
        );
        assert!(
            (delta[k] - dx_ref[k]).abs() / scale < 1e-9,
            "equilibrated step differs from unscaled reference at {}: {} vs {}",
            k,
            delta[k],
            dx_ref[k]
        );
    }

    // κ reduction also holds for this system (≥1e4).
    let cond_before = condition_number(&dense);
    let (d_r, d_c) = equilibrate(&dense);
    let mut scaled = dense.clone();
    for i in 0..6 {
        for j in 0..6 {
            scaled[(i, j)] *= d_r[i] * d_c[j];
        }
    }
    let cond_after = condition_number(&scaled);
    assert!(
        cond_after <= cond_before / 1.0e4,
        "κ reduction ≥1e4 expected: before={:.3e}, after={:.3e}",
        cond_before,
        cond_after
    );
}