//! The 5-step magnetization of the heavy-hex kicked-Ising experiment, against
//! the exactly computed values published for it (θ_h = kπ/32).
//!
//! `cargo run --release --example utility_mz --features quantum_spd -- [threshold]`

use std::time::Instant;
use wai_quantum::quantum_spd::{Pauli, SpdConfig, expectation, heavy_hex_127, kicked_ising};

/// Exact Mz after 5 steps at θ_h = kπ/32, k = 1…16.
const EXACT: [f64; 16] = [
    0.995131, 0.979836, 0.951490, 0.904843, 0.833134, 0.731464, 0.601104, 0.452804, 0.305680, 0.180317, 0.090301,
    0.037263, 0.012260, 0.002936, 0.000324, 0.0,
];

fn main() {
    let threshold: f64 = std::env::args().nth(1).map_or(0.0, |a| a.parse().expect("a number"));
    let kmax: u32 = std::env::args().nth(2).map_or(16, |a| a.parse().expect("a number"));
    let edges = heavy_hex_127();
    let mz: Vec<(f64, Vec<(u32, Pauli)>)> = (0..127).map(|q| (1.0 / 127.0, vec![(q, Pauli::Z)])).collect();
    let mut worst: f64 = 0.0;
    for (k, exact) in (1..=kmax).zip(EXACT) {
        let theta_h = f64::from(k) * std::f64::consts::PI / 32.0;
        let started = Instant::now();
        let r = expectation(&kicked_ising(127, &edges, theta_h, 5), &mz, &[], &SpdConfig { threshold, max_weight: None })
            .expect("a valid circuit");
        let err = (r.value - exact).abs();
        worst = worst.max(err);
        println!(
            "k={k:2}  Mz={:.6}  exact={exact:.6}  |err|={err:.1e}  bound={:.1e}  peak={}  {:.2}s",
            r.value,
            r.truncation_bound,
            r.peak_terms,
            started.elapsed().as_secs_f64()
        );
    }
    println!("worst |err| = {worst:.2e} (the published values carry 6 decimals)");
}