#![cfg(sidereon_repo_tests)]
//! Reweighting-math correctness gate for the Huber/IRLS SPP path.
//!
//! Replays the hand-rolled explicit outer-loop Huber IRLS golden
//! (`fixtures/huber_irls_trace.json`, generated by `gen_huber_irls_trace.py`
//! with scipy 1.17.1 / numpy 2.4.6, NOT `scipy loss='huber'`) on the SAME
//! deterministic linear least-squares case with one injected outlier.
//!
//! Two assertions, mirroring the substrate's bit-exact-vs-agreement split:
//!   - the per-outer-iteration Huber WEIGHT vector matches the golden at
//!     0 ULP (pure `f64` abs/compare/divide via the crate `robust` primitives,
//!     no FMA);
//!   - the converged parameter vector at each iteration matches to the existing
//!     sub-micron AGREEMENT bound (`abs < 1e-6`, the linear-algebra step is
//!     BLAS-bound, not bit-reproducible).
//!
//! The recipe constants (k, MAD constant, scale floor, stop rule) are read from
//! the fixture metadata and asserted against the crate constants so the golden
//! and the implementation can never silently measure different algorithms.

use nalgebra::{DMatrix, DVector};
use serde_json::Value;
use sidereon_core::astro::math::robust::{huber_weight, mad_scale, HUBER_K, MAD_NORMAL_CONST};
use std::path::PathBuf;

/// Sub-micron solver-agreement bound, the same figure used at
/// `parity_least_squares.rs` for the BLAS-bound converged parameter vector.
const AGREEMENT_BOUND: f64 = 1e-6;

fn parse_hex_float(s: &str) -> f64 {
    let s = s.trim();
    match s {
        "-inf" => return f64::NEG_INFINITY,
        "+inf" | "inf" => return f64::INFINITY,
        _ => {}
    }
    let (neg, rest) = match s.strip_prefix('-') {
        Some(r) => (true, r),
        None => (false, s),
    };
    let rest = rest
        .strip_prefix("0x")
        .or_else(|| rest.strip_prefix("0X"))
        .unwrap_or_else(|| panic!("not a hex float (missing 0x): {s:?}"));
    let (mantissa, exp_str) = rest
        .split_once(['p', 'P'])
        .unwrap_or_else(|| panic!("not a hex float (missing p exponent): {s:?}"));
    let exp2: i32 = exp_str
        .parse()
        .unwrap_or_else(|_| panic!("bad exponent {s:?}"));
    let (int_part, frac_part) = match mantissa.split_once('.') {
        Some((i, f)) => (i, f),
        None => (mantissa, ""),
    };
    let int_val: f64 =
        i64::from_str_radix(int_part, 16).unwrap_or_else(|_| panic!("bad int hex {s:?}")) as f64;
    let mut frac_val = 0.0f64;
    let mut scale = 1.0f64 / 16.0;
    for c in frac_part.chars() {
        let d = c
            .to_digit(16)
            .unwrap_or_else(|| panic!("bad frac digit in {s:?}"));
        frac_val += (d as f64) * scale;
        scale /= 16.0;
    }
    let val = (int_val + frac_val) * 2.0f64.powi(exp2);
    if neg {
        -val
    } else {
        val
    }
}

fn float_hex(x: f64) -> String {
    if x == 0.0 {
        return if x.is_sign_negative() {
            "-0x0.0p+0".into()
        } else {
            "0x0.0p+0".into()
        };
    }
    let bits = x.to_bits();
    let sign = if (bits >> 63) & 1 == 1 { "-" } else { "" };
    let exp = ((bits >> 52) & 0x7ff) as i64;
    let mantissa = bits & 0x000f_ffff_ffff_ffff;
    let unbiased = exp - 1023;
    if unbiased >= 0 {
        format!("{sign}0x1.{mantissa:013x}p+{unbiased}")
    } else {
        format!("{sign}0x1.{mantissa:013x}p{unbiased}")
    }
}

fn ordered_i64(x: f64) -> i64 {
    let bits = x.to_bits() as i64;
    if bits < 0 {
        i64::MIN - bits
    } else {
        bits
    }
}

fn ulp_distance(a: f64, b: f64) -> u64 {
    if a.is_nan() || b.is_nan() {
        return u64::MAX;
    }
    ordered_i64(a).abs_diff(ordered_i64(b))
}

fn hex_vec(v: &Value) -> Vec<f64> {
    v.as_array()
        .expect("hex array")
        .iter()
        .map(|e| parse_hex_float(e.as_str().expect("hex string")))
        .collect()
}

fn fixture_path() -> PathBuf {
    PathBuf::from(env!("CARGO_MANIFEST_DIR"))
        .join("tests/fixtures/huber_irls_trace.json")
        .canonicalize()
        .expect("locate huber_irls_trace.json")
}

/// Solve the weighted normal equations `(A^T W A) x = A^T W y` with a dense
/// factorization, mirroring the numpy `linalg.solve` of the generator.
fn solve_weighted(a: &DMatrix<f64>, y: &DVector<f64>, w: &[f64]) -> DVector<f64> {
    let wdiag = DMatrix::from_diagonal(&DVector::from_row_slice(w));
    let lhs = a.transpose() * &wdiag * a;
    let rhs = a.transpose() * &wdiag * y;
    lhs.lu().solve(&rhs).expect("weighted normal eqs solvable")
}

#[test]
fn huber_irls_trace_replay() {
    let raw = std::fs::read_to_string(fixture_path()).expect("read fixture");
    let doc: Value = serde_json::from_str(&raw).expect("parse fixture");

    // Parser self-check.
    let probe = "0x1.5851eb851eb85p+0";
    assert_eq!(
        float_hex(parse_hex_float(probe)),
        probe,
        "hex-float parser round-trip broken"
    );

    // Recipe constants must match the crate's, or the two gates measure
    // different algorithms.
    let meta = &doc["metadata"];
    let k = parse_hex_float(meta["huber_k"].as_str().unwrap());
    let mad_c = parse_hex_float(meta["mad_normal_const"].as_str().unwrap());
    let scale_floor = parse_hex_float(meta["scale_floor_m"].as_str().unwrap());
    let outer_tol = parse_hex_float(meta["outer_tol"].as_str().unwrap());
    let max_outer = meta["max_outer"].as_u64().unwrap() as usize;
    assert_eq!(k, HUBER_K, "fixture k diverges from crate HUBER_K");
    assert_eq!(
        mad_c, MAD_NORMAL_CONST,
        "fixture MAD constant diverges from crate"
    );

    // Rebuild the problem.
    let a_rows: Vec<Vec<f64>> = doc["problem"]["A"]
        .as_array()
        .unwrap()
        .iter()
        .map(hex_vec)
        .collect();
    let n = a_rows.len();
    let ncol = a_rows[0].len();
    let a = DMatrix::from_fn(n, ncol, |i, j| a_rows[i][j]);
    let y = DVector::from_row_slice(&hex_vec(&doc["problem"]["y"]));

    let iters = doc["iterations"].as_array().unwrap();

    // Iteration 0: unweighted warm start (identity weights).
    let w0 = vec![1.0f64; n];
    let mut x = solve_weighted(&a, &y, &w0);
    let x0_exp = hex_vec(&iters[0]["x"]);
    for (i, xe) in x0_exp.iter().enumerate() {
        assert!(
            (x[i] - xe).abs() < AGREEMENT_BOUND,
            "warm-start x[{i}] abs={:e} exceeds {AGREEMENT_BOUND:e}",
            (x[i] - xe).abs()
        );
    }

    let mut weight_ulp_failures: Vec<String> = Vec::new();
    let mut outer_done = 0usize;

    for (it, iter) in iters.iter().enumerate().skip(1).take(max_outer) {
        // The 0-ULP reweighting-math gate is on the pure transform from a GIVEN
        // residual vector to the scale and weights. Feed the golden's recorded
        // residual (the residual itself is BLAS-bound and is checked at the
        // sub-micron x level below), so the weight/scale comparison isolates the
        // pure f64 abs/compare/divide path with no solver disagreement leaking
        // in. Cross-check that our own re-solved residual agrees with it to the
        // sub-micron bound, proving the recorded residual is the same problem.
        let r_golden = hex_vec(&iter["residual"]);
        let r_rust: Vec<f64> = (&y - &a * &x).iter().copied().collect();
        for (i, rg) in r_golden.iter().enumerate() {
            assert!(
                (r_rust[i] - rg).abs() < AGREEMENT_BOUND,
                "outer {it} residual[{i}] agreement abs={:e} exceeds {AGREEMENT_BOUND:e}",
                (r_rust[i] - rg).abs()
            );
        }
        let r = r_golden;
        let s = mad_scale(&r, scale_floor).expect("valid robust scale");
        let w: Vec<f64> = r.iter().map(|&ri| huber_weight(ri / s, k)).collect();

        // ---- 0-ULP Huber WEIGHT vector against the golden ----
        let w_exp = hex_vec(&iter["huber_weight"]);
        for (i, we) in w_exp.iter().enumerate() {
            let ulp = ulp_distance(w[i], *we);
            if ulp != 0 {
                weight_ulp_failures.push(format!(
                    "outer {it} w[{i}]: {ulp} ULP (rust={} ref={})",
                    float_hex(w[i]),
                    float_hex(*we)
                ));
            }
        }
        // The recorded scale must match at 0 ULP too (it drives the weights).
        let s_exp = parse_hex_float(iter["scale_m"].as_str().unwrap());
        let s_ulp = ulp_distance(s, s_exp);
        if s_ulp != 0 {
            weight_ulp_failures.push(format!(
                "outer {it} scale: {s_ulp} ULP (rust={} ref={})",
                float_hex(s),
                float_hex(s_exp)
            ));
        }

        let x_prev = x.clone();
        x = solve_weighted(&a, &y, &w);

        // ---- sub-micron converged x against the golden ----
        let x_exp = hex_vec(&iter["x"]);
        for (i, xe) in x_exp.iter().enumerate() {
            assert!(
                (x[i] - xe).abs() < AGREEMENT_BOUND,
                "outer {it} x[{i}] abs={:e} exceeds {AGREEMENT_BOUND:e} (rust={}, ref={xe})",
                (x[i] - xe).abs(),
                x[i]
            );
        }

        outer_done += 1;
        let dpos =
            ((x[0] - x_prev[0]).powi(2) + (x[1] - x_prev[1]).powi(2) + (x[2] - x_prev[2]).powi(2))
                .sqrt();
        if dpos < outer_tol {
            break;
        }
    }

    assert!(
        weight_ulp_failures.is_empty(),
        "Huber weight/scale not 0-ULP:\n{}",
        weight_ulp_failures.join("\n")
    );
    assert_eq!(
        outer_done,
        iters.len() - 1,
        "outer iteration count diverged from the golden"
    );
}
