//! The 127-qubit heavy-hex kicked-Ising experiment by belief-propagation
//! tensor networks, against published values.
//!
//! `cargo run --release --example utility_bptn --features quantum_bptn -- [chi ...]`
//!
//! For each bond dimension χ: the 5-step magnetization against its exactly
//! computed values (arXiv:2308.05077), and ⟨Z₆₂⟩ after 20 steps against the
//! published belief-propagation tensor-network values at the same angles
//! (arXiv:2306.14887, extrapolated in χ).

use std::time::Instant;
use wai_quantum::quantum_bptn::Tns;
use wai_quantum::quantum_spd::{Pauli, heavy_hex_127, kicked_ising};

/// (k, exact 5-step Mz at θ_h = kπ/32)
const MZ5: [(u32, f64); 4] = [(2, 0.979836), (4, 0.904843), (6, 0.731464), (8, 0.452804)];
/// (θ_h, published 20-step ⟨Z₆₂⟩, extrapolated in χ)
const Z62: [(f64, f64); 4] = [(0.2, 0.9786), (0.4, 0.8498), (0.6, 0.6686), (0.8, 0.2267)];

fn main() {
    let chis: Vec<usize> = std::env::args().skip(1).map(|a| a.parse().expect("a bond dimension")).collect();
    let chis = if chis.is_empty() { vec![8, 16, 32] } else { chis };
    let edges = heavy_hex_127();
    for chi in chis {
        println!("== χ = {chi}");
        for (k, exact) in MZ5 {
            let started = Instant::now();
            let theta = f64::from(k) * std::f64::consts::PI / 32.0;
            let mut t = Tns::new(127, &edges, chi).unwrap();
            t.run(&kicked_ising(127, &edges, theta, 5)).unwrap();
            let mz: f64 = (0..127).map(|q| t.expect1(q, Pauli::Z).unwrap()).sum::<f64>() / 127.0;
            println!(
                "  5-step Mz   θ_h=kπ/32 k={k}  {mz:.6}  exact {exact:.6}  |err| {:.1e}  discarded {:.1e}  {:.1}s",
                (mz - exact).abs(),
                t.discarded(),
                started.elapsed().as_secs_f64()
            );
        }
        for (theta, published) in Z62 {
            let started = Instant::now();
            let mut t = Tns::new(127, &edges, chi).unwrap();
            t.run(&kicked_ising(127, &edges, theta, 20)).unwrap();
            let z = t.expect1(62, Pauli::Z).unwrap();
            println!(
                "  20-step Z62 θ_h={theta}  {z:.4}  published {published:.4}  |diff| {:.1e}  discarded {:.1e}  {:.1}s",
                (z - published).abs(),
                t.discarded(),
                started.elapsed().as_secs_f64()
            );
        }
    }
}