kshana 0.22.0

// SPDX-License-Identifier: AGPL-3.0-only
//! RAIM / ARAIM **statistical-kernel** reference test (external oracle: SciPy).
//!
//! kshana's snapshot-RAIM / solution-separation / ARAIM stack builds its detection
//! thresholds and protection-level multipliers from a dependency-free regularized
//! incomplete-gamma function (`chi2_cdf`/`chi2_quantile`, `normal_cdf`/`normal_quantile`,
//! `noncentral_chi2_cdf`) plus a bisection for the non-centrality `pbias`. This test pins
//! every one of those primitives against **SciPy** (`scipy.stats.chi2`, `.norm`, `.ncx2`
//! and `scipy.optimize.brentq`) — a wholly independent implementation (Cephes/Boost under
//! the hood, a different algorithm from Kshana's Numerical-Recipes series/continued-
//! fraction), evaluated across the false-alarm, missed-detection and redundancy ranges
//! the RAIM stack actually uses. Reproducing SciPy's chi-squared / non-central / normal
//! laws to the tolerances below makes the integrity kernel **externally validated**, not
//! merely self-consistent.
//!
//! Oracle: SciPy 1.17.0 (NumPy 2.4.1). Reference values are the SciPy outputs to full f64.
//! The geometry that wraps this kernel (the `(GᵀG)⁻¹` slope/covariance projection) is
//! validated separately against gnss_lib_py (DOP) and on real IGS SP3 geometry; the
//! ARAIM MHSS *integrity-risk budget allocation* has no published numeric oracle and is
//! tracked honestly as founder-gated in `docs/ARAIM_REFERENCE.md` — not asserted here.

use kshana::raim::{
    chi2_cdf, chi2_quantile, noncentral_chi2_cdf, normal_cdf, normal_quantile, pbias,
};

/// Pass when `got` is within `abstol + reltol·|want|` of `want` (a combined absolute /
/// relative band, so both tail probabilities near 0 and O(10) quantiles are checked
/// meaningfully).
fn close(got: f64, want: f64, abstol: f64, reltol: f64, what: &str) {
    let d = (got - want).abs();
    assert!(
        d <= abstol + reltol * want.abs(),
        "{what}: kshana {got:.12e} vs SciPy {want:.12e} (|Δ|={d:.2e} > {abstol:.0e}+{reltol:.0e}·|want|)"
    );
}

#[test]
fn chi2_cdf_matches_scipy() {
    // (x, k, scipy.stats.chi2.cdf(x, k))
    let cases: &[(f64, f64, f64)] = &[
        (0.5, 1.0, 0.5204998778130466),
        (1.0, 1.0, 0.6826894921370859),
        (2.0, 1.0, 0.8427007929497151),
        (3.841459, 1.0, 0.9500000053468042),
        (5.0, 1.0, 0.9746526813225318),
        (9.21034, 1.0, 0.9975934800522014),
        (15.0, 1.0, 0.9998924888232705),
        (30.0, 1.0, 0.9999999567953695),
        (45.0, 1.0, 0.9999999999802965),
        (0.5, 2.0, 0.22119921692859512),
        (1.0, 2.0, 0.3934693402873665),
        (2.0, 2.0, 0.6321205588285577),
        (3.841459, 2.0, 0.8534999486480763),
        (5.0, 2.0, 0.9179150013761012),
        (9.21034, 2.0, 0.989999998140119),
        (15.0, 2.0, 0.9994469156298522),
        (30.0, 2.0, 0.9999996940976795),
        (45.0, 2.0, 0.9999999998308102),
        (0.5, 3.0, 0.08110858834532418),
        (1.0, 3.0, 0.19874804309879915),
        (2.0, 3.0, 0.42759329552912023),
        (3.841459, 3.0, 0.7208995567559311),
        (5.0, 3.0, 0.8282028557032668),
        (9.21034, 3.0, 0.973378842463947),
        (15.0, 3.0, 0.9981833510334277),
        (30.0, 3.0, 0.9999986199429687),
        (45.0, 3.0, 0.9999999990747298),
        (0.5, 5.0, 0.007876706767370404),
        (1.0, 5.0, 0.03743422675270362),
        (2.0, 5.0, 0.15085496391539038),
        (3.841459, 5.0, 0.42753956337478244),
        (5.0, 5.0, 0.5841198130044919),
        (9.21034, 5.0, 0.8990371607424125),
        (15.0, 5.0, 0.9896376620842136),
        (30.0, 5.0, 0.9999852514189616),
        (45.0, 5.0, 0.9999999854912283),
        (0.5, 8.0, 0.00013336965051406237),
        (1.0, 8.0, 0.001751622556290824),
        (2.0, 8.0, 0.01898815687615381),
        (3.841459, 8.0, 0.1288632972898739),
        (5.0, 8.0, 0.2424238668669339),
        (9.21034, 8.0, 0.675136018272993),
        (15.0, 8.0, 0.9408545401673161),
        (30.0, 8.0, 0.9997886214965332),
        (45.0, 8.0, 0.9999996320016274),
        (0.5, 12.0, 2.7381356338284025e-07),
        (1.0, 12.0, 1.4164937322342495e-05),
        (2.0, 12.0, 0.000594184817581693),
        (3.841459, 12.0, 0.013869813216517646),
        (5.0, 12.0, 0.04202103819530614),
        (9.21034, 12.0, 0.31513265527650813),
        (15.0, 12.0, 0.7585635490297244),
        (30.0, 12.0, 0.997207570667299),
        (45.0, 12.0, 0.9999896949927592),
        (0.5, 20.0, 2.0942485399973623e-13),
        (1.0, 20.0, 1.709670029348906e-10),
        (2.0, 20.0, 1.1142547833872071e-07),
        (3.841459, 20.0, 3.330932704410381e-05),
        (5.0, 20.0, 0.0002773520946208362),
        (9.21034, 20.0, 0.019659538269975428),
        (15.0, 20.0, 0.22359238698028533),
        (30.0, 20.0, 0.9301463393005901),
        (45.0, 20.0, 0.9988965307569716),
    ];
    for &(x, k, want) in cases {
        close(
            chi2_cdf(x, k),
            want,
            1e-9,
            1e-7,
            &format!("chi2_cdf({x},{k})"),
        );
    }
}

#[test]
fn chi2_quantile_matches_scipy() {
    // (p, k, scipy.stats.chi2.ppf(p, k))
    let cases: &[(f64, f64, f64)] = &[
        (0.9, 1.0, 2.705543454095418),
        (0.95, 1.0, 3.8414588206941205),
        (0.99, 1.0, 6.634896601021215),
        (0.999, 1.0, 10.827566170662733),
        (0.99999, 1.0, 19.511420964666268),
        (0.9999999, 1.0, 28.373987362798125),
        (0.9, 2.0, 4.605170185988092),
        (0.95, 2.0, 5.99146454710798),
        (0.99, 2.0, 9.21034037197618),
        (0.999, 2.0, 13.815510557964274),
        (0.99999, 2.0, 23.02585092994956),
        (0.9999999, 2.0, 32.23619130296935),
        (0.9, 3.0, 6.251388631170322),
        (0.95, 3.0, 7.8147279032511765),
        (0.99, 3.0, 11.34486673014437),
        (0.999, 3.0, 16.26623619623813),
        (0.99999, 3.0, 25.90174974567149),
        (0.9999999, 3.0, 35.40575158101893),
        (0.9, 5.0, 9.236356899781121),
        (0.95, 5.0, 11.070497693516351),
        (0.99, 5.0, 15.086272469388987),
        (0.999, 5.0, 20.515005652432876),
        (0.99999, 5.0, 30.85618994044592),
        (0.9999999, 5.0, 40.86302108923746),
        (0.9, 8.0, 13.361566136511728),
        (0.95, 8.0, 15.507313055865453),
        (0.99, 8.0, 20.090235029663233),
        (0.999, 8.0, 26.12448155837614),
        (0.99999, 8.0, 37.33159364443294),
        (0.9999999, 8.0, 47.97246494259998),
        (0.9, 12.0, 18.549347786703244),
        (0.95, 12.0, 21.02606981748307),
        (0.99, 12.0, 26.21696730553585),
        (0.999, 12.0, 32.90949040736021),
        (0.99999, 12.0, 45.076146524171676),
        (0.9999999, 12.0, 56.43366123067437),
        (0.9, 20.0, 28.411980584305635),
        (0.95, 20.0, 31.410432844230925),
        (0.99, 20.0, 37.56623478662506),
        (0.999, 20.0, 45.31474661812586),
        (0.99999, 20.0, 59.0445503868145),
        (0.9999999, 20.0, 71.58930867948122),
    ];
    for &(p, k, want) in cases {
        close(
            chi2_quantile(p, k),
            want,
            1e-6,
            1e-6,
            &format!("chi2_quantile({p},{k})"),
        );
    }
}

#[test]
fn normal_cdf_matches_scipy() {
    // (z, scipy.stats.norm.cdf(z))
    let cases: &[(f64, f64)] = &[
        (-3.0, 0.001349898031630093),
        (-1.0, 0.15865525393145707),
        (0.0, 0.5),
        (1.0, 0.8413447460685429),
        (1.959964, 0.9750000009035577),
        (2.0, 0.9772498680518208),
        (3.0, 0.9986501019683699),
        (4.0, 0.9999683287581669),
        (5.0, 0.9999997133484281),
        (5.3267, 0.9999999499934271),
    ];
    for &(z, want) in cases {
        close(normal_cdf(z), want, 1e-9, 1e-7, &format!("normal_cdf({z})"));
    }
}

#[test]
fn normal_quantile_matches_scipy() {
    // (p, scipy.stats.norm.ppf(p)) — includes the DO-229E/Rayleigh K-factor tails.
    let cases: &[(f64, f64)] = &[
        (0.5, 0.0),
        (0.9, 1.2815515655446004),
        (0.975, 1.959963984540054),
        (0.99, 2.3263478740408408),
        (0.999, 3.090232306167813),
        (0.99999, 4.264890793923841),
        (0.9999999, 5.199337582290661),
        (0.99999995, 5.326723886278398),
        (0.999999995, 5.73072886926708),
        (0.9999999995, 6.109410191663287),
    ];
    for &(p, want) in cases {
        close(
            normal_quantile(p),
            want,
            1e-4,
            1e-6,
            &format!("normal_quantile({p})"),
        );
    }
}

#[test]
fn noncentral_chi2_cdf_matches_scipy() {
    // (x, k, lambda, scipy.stats.ncx2.cdf(x, k, lambda))
    let cases: &[(f64, f64, f64, f64)] = &[
        (1.0, 1.0, 1.0, 0.4772498680518209),
        (5.0, 1.0, 1.0, 0.8911774090184011),
        (10.0, 1.0, 1.0, 0.9846858552504759),
        (20.0, 1.0, 1.0, 0.9997418103789223),
        (40.0, 1.0, 1.0, 0.99999994939973),
        (1.0, 2.0, 4.0, 0.08189230363059402),
        (5.0, 2.0, 4.0, 0.4958760288617339),
        (10.0, 2.0, 4.0, 0.8314310864698685),
        (20.0, 2.0, 4.0, 0.9894594878972105),
        (40.0, 2.0, 4.0, 0.9999860484700465),
        (1.0, 3.0, 9.0, 0.004766081943871699),
        (5.0, 3.0, 9.0, 0.12312784130847292),
        (10.0, 3.0, 9.0, 0.43321512043956073),
        (20.0, 3.0, 9.0, 0.8845101735411985),
        (40.0, 3.0, 9.0, 0.9990278644693452),
        (1.0, 5.0, 16.0, 3.416191268071441e-05),
        (5.0, 5.0, 16.0, 0.00737337689392745),
        (10.0, 5.0, 16.0, 0.07974831517569392),
        (20.0, 5.0, 16.0, 0.4981975393611857),
        (40.0, 5.0, 16.0, 0.973100244008288),
        (1.0, 3.0, 25.0, 4.9054252690677785e-06),
        (5.0, 3.0, 25.0, 0.0011052916318189403),
        (10.0, 3.0, 25.0, 0.018308683698619183),
        (20.0, 3.0, 25.0, 0.22938493093002565),
        (40.0, 3.0, 25.0, 0.8741537512138643),
        (1.0, 8.0, 40.0, 2.0378777781848078e-11),
        (5.0, 8.0, 40.0, 2.922341625994304e-07),
        (10.0, 8.0, 40.0, 4.093966478244554e-05),
        (20.0, 8.0, 40.0, 0.0060038094400554215),
        (40.0, 8.0, 40.0, 0.28925765033413664),
    ];
    for &(x, k, lam, want) in cases {
        close(
            noncentral_chi2_cdf(x, k, lam),
            want,
            1e-9,
            1e-6,
            &format!("ncx2_cdf({x},{k},{lam})"),
        );
    }
}

#[test]
fn pbias_matches_scipy_noncentrality() {
    // (t2, dof, p_md, sqrt(lambda*)) where lambda* solves
    // scipy.stats.ncx2.cdf(t2, dof, lambda*) = p_md (scipy.optimize.brentq).
    // t2 = scipy.stats.chi2.ppf(1 - p_fa, dof) for p_fa in {1e-2, 1e-5}.
    let cases: &[(f64, f64, f64, f64)] = &[
        (6.634896601021215, 1.0, 0.001, 5.66606160971669),
        (6.634896601021215, 1.0, 1e-07, 7.775166885741718),
        (9.21034037197618, 2.0, 0.001, 6.009091854391071),
        (9.21034037197618, 2.0, 1e-07, 8.136814116719645),
        (11.34486673014437, 3.0, 0.001, 6.241795405377199),
        (11.34486673014437, 3.0, 1e-07, 8.384444168227217),
        (15.086272469388987, 5.0, 0.001, 6.580104500822473),
        (15.086272469388987, 5.0, 1e-07, 8.747325265907646),
        (19.511420964666268, 1.0, 0.001, 7.50740571963782),
        (19.511420964666268, 1.0, 1e-07, 9.616510995662832),
        (23.02585092994956, 2.0, 0.001, 7.807486418391467),
        (23.02585092994956, 2.0, 1e-07, 9.926719718434315),
        (25.90174974567149, 3.0, 0.001, 8.023757521044631),
        (25.90174974567149, 3.0, 1e-07, 10.151806144388557),
        (30.85618994044592, 5.0, 0.001, 8.352219555274482),
        (30.85618994044592, 5.0, 1e-07, 10.495652642699305),
    ];
    for &(t2, dof, p_md, want) in cases {
        close(
            pbias(t2, dof, p_md),
            want,
            1e-3,
            1e-4,
            &format!("pbias({t2},{dof},{p_md})"),
        );
    }
}