//! The tracker's peaked circuits: the heavy-hex `49x4020` and `49x5072`, and
//! the all-to-all `P9_Hqap_56x1917`.
//!
//! `cargo run --release --example tracker_peaked --features quantum_engine,quantum_receipt,quantum_bptn,quantum_mpsim,quantum_peaked -- <tracker checkout> <49x4020 | 49x5072 | P9_Hqap_56x1917> <chain | chainverify STRING | structure | verify STRING | junction | profile | run | gauged | mpo | scan> [...]`
//!
//! A peaked circuit prepares a state whose weight on one hidden bit string is
//! orders of magnitude above any other's. These were built as a circuit `R`, a
//! block `U` and its approximate inverse `U†` (re-angled and regrouped, and on
//! all-to-all hardware moved by hidden permutations), and a circuit `P`
//! trained to peak `R`.
//!
//! - `chain [cutoff=2e-3] [chi=8192] [split=…] [unswap=1000000] [rounds=20]
//!   [lookahead=20] [threads=…] [shots=1000] [seed=123]`: the published
//!   method for the all-to-all circuit (`quantum_peaked::chain`). The circuit
//!   is cut at half its blocks, the blocks before the cut absorbed on the
//!   operator's columns and those after it on its rows, each side routed onto
//!   a chain as it goes; the operator is a matrix-product operator truncated
//!   at the cutoff and cap, and when it holds more than `unswap` numbers,
//!   SWAPs that shrink a bond are kept as relabellings, stripping the hidden
//!   permutations. The peak is read on `|0…0⟩` by exact conditional decoding,
//!   checked by `shots` exact samples, and given with its amplitude and the
//!   published string's. `chainverify STRING [cutoffs=…] [splits=…]`: the
//!   string's amplitude in the operator grown at each cutoff from each cut.
//! - `structure [window=4] [chi=16] [cutoff=0.01] [read_chi=8] [unswap=T]`:
//!   reads the circuit through that structure (`quantum_peaked`). Every centre
//!   is probed by growing the operator `window` layers each way; the smallest
//!   marks where `U` meets `U†`. The whole circuit is grown outward from there
//!   as an operator on the graph its gates build, the peak read on `|0…0⟩` by
//!   exact contraction, and the string's amplitude given with the published
//!   string's, where one was printed. `unswap=T` keeps the SWAPs that shrink a
//!   bond of `T` or more, as relabellings.
//! - `verify STRING [chi=64] [cutoffs=…] [centre=…]`: the string's amplitude
//!   `⟨s|C|0…0⟩` in the whole circuit at each cutoff: no readout, no marginals,
//!   the network closed with the string and contracted exactly.
//! - `junction [window=2] [chi=8] [unswap=T]`: the scan alone, every centre's
//!   probe; `profile centre=H [steps=2,4,…]`: one centre's operator after each
//!   number of steps, printed as each finishes.
//!
//! The earlier readings, kept for comparison, swept over `χ`:
//!
//! - `run` and `gauged`: the state, as a tensor network on the coupling graph
//!   by simple update (`gauged` re-gauges after every layer), each qubit's
//!   `⟨Z⟩` through belief propagation. No peak at `χ ≤ 32`.
//! - `mpo`: the circuit as an operator on a chain (each qubit's row and column
//!   side by side, qubits in breadth-first order of the coupling graph), grown
//!   from the middle gate outward. No peak at `χ = 16`: a chain cut crosses
//!   several couplings, each carrying a bond, and their product is past `χ`.
//! - `scan`: the chain operator grown 150 gates a side from 39 junctions.
//!
//! Each run writes a signed receipt to `RECEIPTS_DIR` (default `receipts`).

#[path = "tracker/chain.rs"]
mod chain;
#[path = "tracker/common.rs"]
mod common;
#[path = "tracker/operator.rs"]
mod operator;
#[path = "tracker/structure.rs"]
mod structure;

use serde_json::json;
use wai_quantum::quantum_bptn::Tns;
use wai_quantum::quantum_spd::{Pauli, RotCircuit};

fn main() {
    let mut args = std::env::args().skip(1);
    let root = common::tracker_root(args.next());
    let name = args.next().unwrap_or_else(|| "49x4020".into());
    let how = args.next().unwrap_or_else(|| "run".into());
    let rest: Vec<String> = args.collect();
    let path = format!("data/classically-verifiable-problems/circuit-models/peaked_circuit/peaked_circuit_{}.qasm", structure::instance(&name));
    let src = String::from_utf8(common::read_pinned(&root, &path)).expect("UTF-8");
    let (circuit, map) = common::compact(&common::rotations(&src));
    let edges = common::coupling(&circuit);
    println!("instance        peaked_circuit_{}", structure::instance(&name));
    println!("qubits          {}, {} gates, {} couplings", circuit.n, circuit.ops.len(), edges.len());
    if how == "structure" {
        return structure::structure(&circuit, &map, &name, &path, &rest);
    }
    if how == "junction" {
        return structure::junction(&circuit, &rest);
    }
    if how == "profile" {
        return structure::profile(&circuit, &rest);
    }
    if how == "chain" {
        return chain::chain(&circuit, &map, &name, &path, &rest);
    }
    if how == "chainverify" {
        let string = rest.first().expect("the string to verify, bit q the qubit q");
        return chain::verify(&circuit, &map, &name, &path, string, &rest[1..]);
    }
    if how == "verify" {
        let string = rest.first().expect("the string to verify, bit q the qubit q");
        structure::verify(&circuit, &map, &name, &path, string, &rest[1..]);
        return;
    }
    let chis: Vec<usize> = rest.iter().map(|a| a.parse().expect("a bond dimension")).collect();
    let chis = if chis.is_empty() { vec![8, 16, 32] } else { chis };
    if how == "mpo" {
        return mpo(&circuit, &map, &edges, &name, &path, &chis);
    }
    if how == "scan" {
        // Where do the circuit's halves undo each other? The operator grown
        // 150 two-qubit gates a side from each of 39 junctions, at the first χ.
        let chi = chis[0];
        for k in 1..40 {
            let mid = circuit.ops.len() * k / 40;
            let g = operator::grow_from(&circuit, &edges, chi, mid, 150).expect("a window");
            println!("scan chi {chi} junction {mid:>6} ({:>4.1}%)  bond {}  discarded {:.3e}", 100.0 * k as f64 / 40.0, g.max_bond, g.discarded);
        }
        return;
    }
    for chi in chis {
        let ((z, discarded, bond), secs, energy) = common::metered(|| {
            let mut t = Tns::new(circuit.n, &edges, chi).expect("the graph");
            match how.as_str() {
                "run" => t.run(&circuit).expect("gates on the graph"),
                "gauged" => t.run_gauged(&circuit, 20, 1e-10).expect("gates on the graph"),
                other => panic!("run or gauged, not {other}"),
            }
            let env = t.belief_propagation(500, 1e-12);
            let z: Vec<f64> = (0..circuit.n).map(|q| t.expect1_bp(&env, q, Pauli::Z).expect("a qubit")).collect();
            (z, t.discarded(), t.max_bond())
        });
        // Bit q of the string is qubit q's (device order), 1 where ⟨Z⟩ < 0.
        let bits: String = z.iter().map(|&v| if v < 0.0 { '1' } else { '0' }).collect();
        let weakest = z.iter().map(|v| v.abs()).fold(f64::INFINITY, f64::min);
        println!("chi {chi:<4} string {bits}  min |<Z>| {weakest:.3}  discarded {discarded:.3e}  bond {bond}  {secs:.1} s");
        let request = json!({
            "instance": format!("peaked_circuit_heavy_hex_{name}"),
            "source": common::source(&path),
            "method": if how == "run" {
                "tensor network on the coupling graph by simple update (quantum_bptn::Tns::run), each qubit's ⟨Z⟩ through belief propagation"
            } else {
                "tensor network on the coupling graph by simple update re-gauged after every layer (quantum_bptn::Tns::run_gauged), each qubit's ⟨Z⟩ through belief propagation"
            },
            "chi": chi,
            "qubits": map,
        });
        let result = json!({
            "z": z,
            "string": bits,
            "min_abs_z": weakest,
            "discarded": discarded,
            "max_bond": bond,
            "wall_seconds": secs,
            "energy_figure": energy,
        });
        let out = common::receipt_path(&format!("peaked_{name}_{how}_chi{chi}.json"));
        assert!(common::verifies(&common::seal(request, result, Some(&out))));
        println!("receipt         {}", out.display());
    }
}

/// The circuit as an operator ([`operator::grow_from_middle`]), at each `χ`.
fn mpo(circuit: &RotCircuit, map: &[u32], edges: &[(u32, u32)], name: &str, path: &str, chis: &[usize]) {
    for &chi in chis {
        let (g, secs, energy) = common::metered(|| operator::grow_from_middle(circuit, edges, chi));
        let bits: String = g.string.iter().map(|&b| if b == 1 { '1' } else { '0' }).collect();
        let weakest = g.z.iter().map(|v| v.abs()).fold(f64::INFINITY, f64::min);
        println!(
            "mpo chi {chi:<4} string {bits}  p {:.4e}  min |<Z>| {weakest:.3}  bond {}  discarded {:.3e}  fidelity estimate {:.3e}  {secs:.1} s",
            g.probability, g.max_bond, g.discarded, g.fidelity_estimate
        );
        let request = json!({
            "instance": format!("peaked_circuit_heavy_hex_{name}"),
            "source": common::source(path),
            "method": "the circuit as an operator, a matrix-product state on row and column sites grown from the middle outward (quantum_mpsim::Chain), its columns projected on |0⟩; each qubit's ⟨Z⟩ and the string's probability",
            "chi": chi,
            "chain_order": g.order.iter().map(|&q| map[q as usize]).collect::<Vec<_>>(),
            "qubits": map,
        });
        let result = json!({
            "z": g.z,
            "string": bits,
            "probability": g.probability,
            "min_abs_z": weakest,
            "max_bond": g.max_bond,
            "discarded": g.discarded,
            "fidelity_estimate": g.fidelity_estimate,
            "wall_seconds": secs,
            "energy_figure": energy,
        });
        let out = common::receipt_path(&format!("peaked_{name}_mpo_chi{chi}.json"));
        assert!(common::verifies(&common::seal(request, result, Some(&out))));
        println!("receipt         {}", out.display());
    }
}