//! The 127-qubit heavy-hex kicked-Ising experiment against published values.
//!
//! - The 5-step magnetization `Mz = (1/127) Σ ⟨Z_q⟩` at `θ_h = kπ/32` was
//!   computed exactly and published as `exact.csv` with arXiv:2308.05077.
//!   With the truncation threshold at 0 this crate reproduces `k = 1` to
//!   3.5e-7 (0.995131 against 0.995131, 357M terms); here a threshold of
//!   1e-7 keeps the test fast, and the error is checked against the bound.
//! - ⟨Z₆₂⟩ after 20 steps at `θ_h = 2π/32`: the converged classical values
//!   published with the same paper (sparse Pauli dynamics, a mixed method,
//!   PEPS, PEPO) span 0.9776–0.9796.
#![cfg(feature = "quantum_spd")]

use wai_quantum::quantum_spd::{Pauli, SpdConfig, expect_z, expectation, heavy_hex_127, kicked_ising};

/// The same magnetization by belief-propagation tensor networks at χ = 16, an
/// independent method: the light cone of five steps fits, and the values are
/// exact to rounding of the published six decimals.
#[cfg(feature = "quantum_bptn")]
#[test]
fn five_step_magnetization_by_tensor_networks_matches_the_exact_values() {
    use wai_quantum::quantum_bptn::Tns;
    let edges = heavy_hex_127();
    for (k, exact) in [(2, 0.979836), (4, 0.904843)] {
        let mut t = Tns::new(127, &edges, 16).unwrap();
        t.run(&kicked_ising(127, &edges, theta(k), 5)).unwrap();
        let mz: f64 = (0..127).map(|q| t.expect1(q, Pauli::Z).unwrap()).sum::<f64>() / 127.0;
        assert!((mz - exact).abs() < 1e-6, "k = {k}: Mz = {mz} against {exact}");
    }
}

fn theta(k: u32) -> f64 {
    f64::from(k) * std::f64::consts::PI / 32.0
}

#[test]
fn five_step_magnetization_matches_the_exact_values() {
    let mz: Vec<(f64, Vec<(u32, Pauli)>)> = (0..127).map(|q| (1.0 / 127.0, vec![(q, Pauli::Z)])).collect();
    for (k, exact) in [(1, 0.995131), (2, 0.979836)] {
        let c = kicked_ising(127, &heavy_hex_127(), theta(k), 5);
        let r = expectation(&c, &mz, &[], &SpdConfig { threshold: 1e-7, max_weight: None }).unwrap();
        let err = (r.value - exact).abs();
        // The published values carry six decimals: allow their rounding.
        assert!(err <= r.truncation_bound + 5e-7, "k = {k}: error {err:e} exceeds the bound {:e}", r.truncation_bound);
        assert!(err < 3e-5, "k = {k}: Mz = {} against {exact}", r.value);
    }
}

#[test]
fn twenty_step_z62_lies_among_the_converged_classical_values() {
    let c = kicked_ising(127, &heavy_hex_127(), theta(2), 20);
    let r = expect_z(&c, 62, &SpdConfig { threshold: 1e-4, max_weight: None }).unwrap();
    assert!((0.9776..=0.9796).contains(&r.value), "⟨Z₆₂⟩ = {}", r.value);
}