integral 0.1.0

//! Cross-engine validation: the OS/HGP and Rys ERI engines must be
//! **transparent** — same Coulomb integral, same `(a,b,c,d)` block layout — and
//! the `Auto` dispatch must equal whichever engine it picks.
//!
//! Pure Rust, no external library. These tests lock the
//! engine-to-engine agreement, the 8-fold symmetry and falsifiable layout from
//! independently built OS blocks, and the dispatch-boundary behaviour.

use integral::{select_engine, Basis, Engine, Shell};
use integral_math::am::n_cart;

#[test]
fn dispatch_policy_matches_5c_calibration() {
    // The 2D crossover (DESIGN_NOTES.md D13): OS/HGP
    // wins down to deg 1 for l_total ≤ 4, then once deg ≥ 16 up through l_total
    // 16; l_total ≥ 17 stays Rys (no measured OS win at high contraction there).
    assert_eq!(select_engine(0, 1), Engine::OsHgp); // L0: OS even uncontracted
    assert_eq!(select_engine(4, 1), Engine::OsHgp); // l_total 4, deg 1 → OS (was Rys)
    assert_eq!(select_engine(8, 1), Engine::Rys); // l_total 8, deg 1 → Rys
    assert_eq!(select_engine(8, 15), Engine::Rys); // below deg 16 → Rys
    assert_eq!(select_engine(8, 16), Engine::OsHgp); // deg 16 (K≥2) → OS (was 81)
    assert_eq!(select_engine(16, 16), Engine::OsHgp); // cap relaxed up to l_total 16
    assert_eq!(select_engine(16, 1), Engine::Rys); // …but deg 1 still Rys there
    assert_eq!(select_engine(17, usize::MAX), Engine::Rys); // l_total ≥ 17: always Rys
}

/// `|os − rys| ≤ atol + rtol·|rys|` with `atol = 1e-11`, `rtol = 1e-10`. The two
/// engines are independent f64 recurrences agreeing to ~1e-12 on dominant
/// elements; the atol floor absorbs benign near-cancellation on tiny components.
fn assert_close(os: &[f64], rys: &[f64], tag: &str) {
    assert_eq!(os.len(), rys.len(), "{tag}: length mismatch");
    const ATOL: f64 = 1e-11;
    const RTOL: f64 = 1e-10;
    let mut worst = 0.0_f64;
    for (i, (o, r)) in os.iter().zip(rys.iter()).enumerate() {
        let d = (o - r).abs();
        worst = worst.max(d);
        assert!(
            d <= ATOL + RTOL * r.abs(),
            "{tag}[{i}]: OS {o} vs Rys {r} (Δ={d:e})"
        );
    }
    assert!(worst.is_finite(), "{tag}: non-finite Δ");
}

/// s (contracted), p, d, f on four distinct centres — every (ij|kl) block is
/// genuinely distinct, and several shells are contracted.
fn mixed_basis() -> Basis {
    Basis::new(vec![
        Shell::new(0, [0.0, 0.0, 0.0], vec![1.4, 0.5], vec![0.6, 0.5]).unwrap(), // s, 2 prim
        Shell::new(1, [0.6, -0.3, 0.2], vec![0.9, 0.35], vec![0.55, 0.5]).unwrap(), // p, 2 prim
        Shell::new(2, [-0.4, 0.7, -0.1], vec![1.1], vec![1.0]).unwrap(),         // d
        Shell::new(3, [0.2, 0.5, 0.8], vec![0.7], vec![1.0]).unwrap(),           // f
    ])
}

#[test]
fn os_matches_rys_block_sweep() {
    // A spread of quartets with mixed L on four distinct centres, including the
    // mixed-high-L (dd|ff) case and contracted indices (shells 0 and 1).
    let basis = mixed_basis();
    let quartets = [
        (0, 0, 0, 0),
        (1, 0, 0, 0),
        (0, 1, 0, 0),
        (1, 1, 1, 1),
        (2, 0, 1, 0),
        (2, 1, 2, 1),
        (0, 1, 2, 3), // (sp|df) — four different L
        (3, 0, 2, 1),
        (2, 2, 3, 3), // (dd|ff) — mixed high-L
        (3, 3, 3, 3), // (ff|ff)
    ];
    for (i, j, k, l) in quartets {
        let os = basis.eri_block_with(Engine::OsHgp, i, j, k, l);
        let rys = basis.eri_block_with(Engine::Rys, i, j, k, l);
        let s = basis.shells();
        assert_close(
            &os,
            &rys,
            &format!("(l{} l{}|l{} l{})", s[i].l(), s[j].l(), s[k].l(), s[l].l()),
        );
    }
}

#[test]
fn os_matches_rys_diffuse_tight() {
    // Extreme exponent ratio (1e3 / 1e-3) on a (pp|pp) quartet — the conditioning
    // stress case. Both engines must still agree.
    let basis = Basis::new(vec![
        Shell::new(1, [0.0, 0.0, 0.0], vec![1.0e3], vec![1.0]).unwrap(),
        Shell::new(1, [0.0, 0.0, 1.2], vec![1.0e-3], vec![1.0]).unwrap(),
    ]);
    let os = basis.eri_block_with(Engine::OsHgp, 0, 1, 0, 1);
    let rys = basis.eri_block_with(Engine::Rys, 0, 1, 0, 1);
    assert_close(&os, &rys, "diffuse×tight (pp|pp)");
}

#[test]
fn os_matches_rys_high_contraction() {
    // High contraction degree on every centre (the regime OS/HGP targets): a
    // 6-primitive s and a 4-primitive p, so the (sp|sp) quartet contracts
    // 6·4·6·4 = 576 primitive quartets behind a single early-contraction HRR.
    let basis = Basis::new(vec![
        Shell::new(
            0,
            [0.0, 0.0, 0.0],
            vec![30.0, 10.0, 3.0, 1.0, 0.4, 0.15],
            vec![0.02, 0.08, 0.2, 0.35, 0.3, 0.12],
        )
        .unwrap(),
        Shell::new(
            1,
            [0.3, 0.2, 1.0],
            vec![1.2, 0.6, 0.3, 0.12],
            vec![0.3, 0.4, 0.25, 0.1],
        )
        .unwrap(),
    ]);
    let os = basis.eri_block_with(Engine::OsHgp, 0, 1, 0, 1);
    let rys = basis.eri_block_with(Engine::Rys, 0, 1, 0, 1);
    assert_close(&os, &rys, "high-contraction (sp|sp)");
}

#[test]
fn eight_fold_symmetry_from_os_blocks() {
    // The OS/HGP tensor, built block-by-block with no symmetrization, must obey
    // the 8-fold permutational symmetry of real Cartesian GTOs.
    let basis = mixed_basis();
    let n = basis.nao_cart();
    let eri = basis.eri_with(Engine::OsHgp);
    let at = |i: usize, j: usize, k: usize, l: usize| eri[((i * n + j) * n + k) * n + l];

    let mut worst = 0.0_f64;
    for i in 0..n {
        for j in 0..n {
            for k in 0..n {
                for l in 0..n {
                    let base = at(i, j, k, l);
                    let perms = [
                        at(j, i, k, l),
                        at(i, j, l, k),
                        at(j, i, l, k),
                        at(k, l, i, j),
                        at(l, k, i, j),
                        at(k, l, j, i),
                        at(l, k, j, i),
                    ];
                    for &p in &perms {
                        let d = (base - p).abs();
                        worst = worst.max(d / base.abs().max(1e-300));
                        // atol floor: the same integral evaluated as (ij|kl) vs
                        // (kl|ij) differs only by f64 recurrence-ordering round-off
                        // (~1e-17 abs on near-cancellation elements); rtol catches
                        // any genuine asymmetry.
                        assert!(
                            d <= 1e-12 + 1e-10 * base.abs(),
                            "OS 8-fold symmetry broken at ({i}{j}|{k}{l}): {base} vs {p} (Δ={d:e})"
                        );
                    }
                }
            }
        }
    }
    assert!(worst < 1e-10, "worst OS symmetry residual {worst:e}");
}

#[test]
fn os_block_layout_matches_dense_strides_and_is_falsifiable() {
    // The forced-OS `eri_block` must equal the forced-OS dense slice under the
    // documented row-major (a,b,c,d) layout; a deliberate (c,d) transpose must
    // make the comparison fail (so the layout test has teeth).
    let basis = mixed_basis();
    let n = basis.nao_cart();
    let dense = basis.eri_with(Engine::OsHgp);
    let shells = basis.shells();
    let offs: Vec<usize> = {
        let mut acc = 0;
        let mut v = Vec::new();
        for s in shells {
            v.push(acc);
            acc += s.n_cart();
        }
        v
    };

    let (qi, qj, qk, ql) = (1usize, 0usize, 2usize, 3usize); // (ps|df), all distinct L
    let block = basis.eri_block_with(Engine::OsHgp, qi, qj, qk, ql);
    let (na, nb, nc, nd) = (
        n_cart(shells[qi].l()),
        n_cart(shells[qj].l()),
        n_cart(shells[qk].l()),
        n_cart(shells[ql].l()),
    );
    assert_eq!(block.len(), na * nb * nc * nd);

    let mut transpose_would_differ = false;
    for a in 0..na {
        for b in 0..nb {
            for c in 0..nc {
                for d in 0..nd {
                    let src = ((a * nb + b) * nc + c) * nd + d;
                    let (i, j, k, l) = (offs[qi] + a, offs[qj] + b, offs[qk] + c, offs[ql] + d);
                    let dval = dense[((i * n + j) * n + k) * n + l];
                    assert!(
                        (block[src] - dval).abs() <= 1e-12 * dval.abs().max(1e-12),
                        "OS layout mismatch at {src} ({i}{j}|{k}{l}): {} vs {}",
                        block[src],
                        dval
                    );
                    // The deliberately wrong (c,d)-transposed index: if it differs
                    // from the correct value anywhere, the layout test is sensitive.
                    if c != d && nc == nd {
                        let wrong = ((a * nb + b) * nd + d) * nc + c;
                        if (block[wrong] - dval).abs() > 1e-9 * dval.abs().max(1e-12) {
                            transpose_would_differ = true;
                        }
                    }
                }
            }
        }
    }
    // For (ps|df), nc (d=6) ≠ nd (f=10), so the transpose is a different shape;
    // assert the block is non-trivial instead (the dense-slice equality above is
    // the falsifiable layout check — a transposed engine output would fail it).
    let _ = transpose_would_differ;
    assert!(
        block.iter().any(|&x| x.abs() > 1e-6),
        "OS block unexpectedly ~zero"
    );
}

#[test]
fn auto_dispatch_equals_both_forced_engines() {
    // Across the crossover boundary, `Auto` must (a) equal the engine it selects
    // and (b) since both engines agree, equal *both* forced results. These quartets
    // straddle the policy (DESIGN_NOTES.md D13): the contracted low-L ones
    // dispatch to OS/HGP, the low-contraction higher-L ones to Rys — yet `Auto`
    // matches both. (The exhaustive per-boundary fingerprint is in `dispatch.rs`.)
    let basis = mixed_basis();
    let quartets = [
        (0, 0, 0, 0), // l_total 0
        (1, 1, 0, 0), // 2
        (1, 1, 1, 0), // 3
        (1, 1, 1, 1), // 4  (boundary, OS side)
        (2, 1, 1, 1), // 5  (Rys side)
        (2, 2, 3, 3), // 10 (Rys side)
    ];
    for (i, j, k, l) in quartets {
        let auto = basis.eri_block_with(Engine::Auto, i, j, k, l);
        let os = basis.eri_block_with(Engine::OsHgp, i, j, k, l);
        let rys = basis.eri_block_with(Engine::Rys, i, j, k, l);
        let s = basis.shells();
        let tag = format!(
            "auto-vs-os (l{} l{}|l{} l{})",
            s[i].l(),
            s[j].l(),
            s[k].l(),
            s[l].l()
        );
        assert_close(&auto, &os, &tag);
        let tag2 = format!(
            "auto-vs-rys (l{} l{}|l{} l{})",
            s[i].l(),
            s[j].l(),
            s[k].l(),
            s[l].l()
        );
        assert_close(&auto, &rys, &tag2);
    }
}