//! The crate promises identical results natively, in the browser and at the
//! edge. Before every transcendental moved in-crate, four of these layers gave
//! different bits natively and under `wasm32-wasip2`. These digests were taken
//! natively and confirmed under `wasm32-wasip2`; any platform whose math
//! differs fails here, not in a user's diff.
//!
//! `CARGO_TARGET_WASM32_WASIP2_RUNNER=wasmtime cargo test --release --target wasm32-wasip2 --features full`
#![cfg(all(
    feature = "quantum",
    feature = "quantum_vml",
    feature = "quantum_phasor",
    feature = "quantum_qfhrr",
    feature = "quantum_kernel",
    feature = "quantum_mps",
    feature = "quantum_pauli",
    feature = "quantum_spd",
    feature = "quantum_bptn",
    feature = "quantum_resource",
    feature = "quantum_frame",
    feature = "quantum_tn",
    feature = "quantum_gbs",
    feature = "quantum_sv",
    feature = "quantum_arch",
    feature = "quantum_ldpc",
    feature = "quantum_qio",
    feature = "quantum_uf",
    feature = "quantum_noisy",
    feature = "quantum_dmrg",
    feature = "quantum_tdvp",
    feature = "quantum_chem",
    feature = "quantum_nqs"
))]

use wai_quantum::*;

/// FNV-1a over every value's IEEE-754 bits.
fn digest(xs: &[f64]) -> u64 {
    let mut h: u64 = 0xcbf2_9ce4_8422_2325;
    for x in xs {
        for b in x.to_bits().to_le_bytes() {
            h = (h ^ u64::from(b)).wrapping_mul(0x100_0000_01b3);
        }
    }
    h
}

fn trained() -> quantum_vml::BornFit {
    quantum_vml::train_born(&quantum_vml::parity_target(3), 3, 2, 40, 0.1, 7)
}

#[test]
fn born_machine_training() {
    let fit = trained();
    let mut v = fit.params.clone();
    v.extend(&fit.loss_curve);
    v.extend(&fit.final_probs);
    v.push(fit.final_kl);
    assert_eq!(digest(&v), 0xd8e5_bd02_2201_2e5d);
}

#[test]
fn phasor_bridge() {
    let p = quantum_phasor::random_params(3, 5);
    let (a, b, c, d) = quantum_phasor::compare(3, &p, 64);
    let mut v = a;
    v.extend(b);
    v.extend(c);
    v.push(d);
    assert_eq!(digest(&v), 0x08be_30dd_7ce1_9561);
}

#[test]
fn qfhrr_born_machine() {
    let fit = trained();
    let mut v = quantum_qfhrr::born_probs(3, 2, &fit.params, 4);
    v.push(quantum_qfhrr::tv_vs_exact(3, 2, &fit.params, 4));
    assert_eq!(digest(&v), 0x5407_7c4d_1be0_d985);
}

#[test]
fn phasor_kernel() {
    let (x, y) = quantum_kernel::circles(200, 3);
    let rff = quantum_kernel::sample_rff(2, 64, 1.0, 9);
    let w = quantum_kernel::fit_ridge(&rff, &x, &y, 1e-3);
    let mut v: Vec<f64> = x.iter().flatten().copied().collect();
    v.extend(&w);
    assert_eq!(digest(&v), 0x58ba_1666_5ef6_fd9e);
}

fn brickwork(gs: quantum::Gateset) -> quantum::Circuit {
    let mut c = quantum::Circuit::with_gateset(8, gs);
    for l in 0..4u8 {
        for q in 0..8 {
            c.h(q);
            c.t(q);
        }
        for q in (l % 2..7).step_by(2) {
            c.cx(q, q + 1);
        }
    }
    c
}

#[test]
fn mps_truncation() {
    let c = brickwork(quantum::Gateset::V1);
    let (mps, _) = quantum_mps::run(&c, 4).unwrap();
    assert_eq!(digest(&[quantum_mps::fidelity(&mps, &c.simulate().unwrap())]), 0xa807_44ca_a97f_bff8);
}

/// The same truncation against the circuit2 reference statevector. The MPS
/// itself is an ideal-angle simulation and does not read the phase table, so
/// only the reference it is compared with changes.
#[test]
fn mps_truncation_v2() {
    let c = brickwork(quantum::Gateset::V2);
    let (mps, _) = quantum_mps::run(&c, 4).unwrap();
    assert_eq!(digest(&[quantum_mps::fidelity(&mps, &c.simulate().unwrap())]), 0xe408_42e9_2c89_20ed);
}

#[test]
fn sparse_pauli_propagation() {
    let (e, _) = quantum_pauli::expect_z(&brickwork(quantum::Gateset::V1), 3, 1e-6).unwrap();
    assert_eq!(digest(&[e]), 0x8cc6_5829_1fb1_8909);
}

#[test]
fn sparse_pauli_dynamics() {
    let chain: Vec<(u32, u32)> = (0..29).map(|q| (q, q + 1)).collect();
    let circuit = quantum_spd::kicked_ising(30, &chain, 0.5, 6);
    let cfg = quantum_spd::SpdConfig { threshold: 1e-6, max_weight: None };
    let r = quantum_spd::expect_z(&circuit, 7, &cfg).unwrap();
    assert_eq!(digest(&[r.value, r.truncation_bound]), 0xfc44_f4b6_0e70_68b2);
}

#[test]
fn tensor_network_evolution() {
    // Every qubit's ⟨Z⟩ and the discarded weight after three steps on the
    // 127-qubit lattice at χ = 8: QR, the converged SVD and the parallel
    // layers all feed these bits.
    let edges = quantum_spd::heavy_hex_127();
    let mut t = quantum_bptn::Tns::new(127, &edges, 8).unwrap();
    t.run(&quantum_spd::kicked_ising(127, &edges, 0.6, 3)).unwrap();
    let mut v: Vec<f64> = (0..127).map(|q| t.expect1(q, quantum_spd::Pauli::Z).unwrap()).collect();
    v.push(t.discarded());
    assert_eq!(digest(&v), 0x0e51_33ff_7987_aa26);
}

/// Every figure of the factoring and chemistry estimates on all six reference
/// qubits: distances, factory rounds, qubit counts, times and error rates.
/// The rotation T-count goes through the in-crate logarithm.
#[test]
fn resource_estimates() {
    use quantum_resource::*;
    let programs = [
        (LogicalCounts { qubits: 12_581, measurements: 1_080_000_000, rotations: 12, rotation_depth: 12, t_gates: 12, ccz: 3_730_000_000 }, 1.0 / 3.0),
        (LogicalCounts { qubits: 1318, measurements: 1_370_000_000, rotations: 206_000_000, rotation_depth: 205_000_000, t_gates: 55_300_000, ccz: 135_000_000_000 }, 0.01),
    ];
    let mut v = Vec::new();
    for (c, eps) in programs {
        for q in PhysicalQubit::presets() {
            let e = estimate_best(&c, &q, &EstimateConfig::new(eps)).unwrap();
            v.extend([e.distance as f64, e.steps as f64, e.runtime_ns, e.logical_error, e.t_target]);
            v.extend([e.factories as f64, e.physical_qubits as f64, e.isa.t_per_rotation as f64]);
            for r in &e.factory.unwrap().rounds {
                v.extend([r.distance as f64, r.copies as f64, r.duration_ns, r.failure, r.output_error]);
            }
        }
    }
    assert_eq!(digest(&v), 0xce9d_43b2_d498_f379);
}

/// Detection counts of a sampled surface-code run, the failures its decoding
/// leaves, and the code model fitted from a small grid: the sampler's random
/// stream, its logarithm-driven skips, the matching and the fit's `exp`.
#[test]
fn frame_sampling_and_matching() {
    use quantum_frame::{parse, sample, surface_code_memory};
    use quantum_match::{fit_code_model, memory_experiment};
    let c = parse(&surface_code_memory(3, 3, 0.01)).unwrap();
    let s = sample(&c, 10_000, 5);
    let mut v: Vec<f64> = (0..s.detectors).map(|d| s.count(d) as f64).collect();
    let mut pts = Vec::new();
    for d in [3, 5] {
        for p in [0.004, 0.008] {
            let t = memory_experiment(d, d, p, 4_000, 9).unwrap();
            v.push(t.failures as f64);
            v.push(t.per_round());
            pts.push(t);
        }
    }
    let f = fit_code_model(&pts).unwrap();
    v.extend([f.prefactor, f.threshold]);
    assert_eq!(digest(&v), 0xcbb3_99f9_6d5c_e8f0);
}

/// A contraction plan found by the seeded hyper-optimiser, and the amplitude
/// its sliced, threaded contraction gives: the search's random stream and
/// logarithms, and the fixed-order complex arithmetic.
#[test]
fn tensor_network_plan_and_amplitude() {
    use quantum_tn::{complex, contract_plan, find_path, Families, Gate, Network, PathOptions};
    // A 3 × 4 grid, ten layers of fixed two-qubit gates and phased rotations.
    let (w, h) = (4u32, 3u32);
    let mut gates = Vec::new();
    let mut k = 0.0f64;
    for layer in 0..10u32 {
        for q in 0..w * h {
            k += 0.37;
            let (c, s) = (1.0 / (1.0 + k * k).sqrt(), k / (1.0 + k * k).sqrt());
            gates.push(Gate { qubits: vec![q], matrix: vec![complex(c, 0.0), complex(0.0, -s), complex(0.0, -s), complex(c, 0.0)] });
        }
        let horizontal = layer % 2 == 0;
        for y in 0..h {
            for x in 0..w {
                let (a, b) = if horizontal { (y * w + x, y * w + x + 1) } else { (y * w + x, (y + 1) * w + x) };
                let ok = if horizontal { x + 1 < w && (x + layer / 2) % 2 == 0 } else { y + 1 < h && (y + layer / 2) % 2 == 0 };
                if ok {
                    // An iSWAP-like gate with a controlled phase.
                    let z = complex(0.0, 0.0);
                    let one = complex(1.0, 0.0);
                    let mi = complex(0.0, -1.0);
                    let ph = complex(0.5, -0.866_025_403_784_438_6);
                    gates.push(Gate { qubits: vec![a, b], matrix: vec![one, z, z, z, z, z, mi, z, z, mi, z, z, z, z, z, ph] });
                }
            }
        }
    }
    let n = w * h;
    let bits: Vec<u8> = (0..n).map(|q| (q % 3 == 1) as u8).collect();
    let mut net = Network::amplitude(n, &gates, &bits, &[]).unwrap();
    net.simplify();
    let shape = net.shape();
    let free = find_path(&shape, &net.dims, &net.open, &PathOptions { trials: 12, seed: 5, max_size: None, reconfigure: 6, families: Families::Both });
    let cap = free.per_slice.max_size / 8.0;
    let plan = find_path(&shape, &net.dims, &net.open, &PathOptions { trials: 12, seed: 5, max_size: Some(cap), reconfigure: 6, families: Families::Both });
    let amp = contract_plan(&net, &plan).unwrap().data[0];
    let steps: Vec<f64> = plan.path.steps.iter().flat_map(|&(a, b)| [a as f64, b as f64]).collect();
    let mut v = vec![amp.re, amp.im, plan.per_slice.flops, plan.per_slice.max_size, plan.slices];
    v.extend(steps);
    assert_eq!(digest(&v), 0xf39b_7370_3ae9_ea0c);
}

/// A 24 × 24 hafnian and a 14-click torontonian (both summed in fixed chunks,
/// threaded natively and sequential under WASI), lossy-state probabilities,
/// exact samples and marginals.
#[test]
fn gaussian_boson_sampling() {
    use quantum_gbs::{hafnian, random_unitary, GaussianState, Mat, C};
    let u = random_unitary(24, 5);
    let mut m = Mat::zeros(24);
    for i in 0..24 {
        for j in 0..24 {
            let (x, y) = (u.at(i, j), u.at(j, i));
            m.a[i * 24 + j] = C { re: x.re + y.re, im: x.im + y.im };
        }
    }
    let h = hafnian(&m);
    let r: Vec<f64> = (0..16).map(|i| 0.3 + 0.04 * i as f64).collect();
    let st = GaussianState::squeezed(&r, &random_unitary(16, 6), 0.7).unwrap();
    let pattern: Vec<bool> = (0..16).map(|i| i % 8 != 3).collect();
    let mut v = vec![h.re, h.im, st.click_probability(&pattern).unwrap()];
    let small = GaussianState::squeezed(&[0.8, 0.5, 0.6, 0.4, 0.7], &random_unitary(5, 7), 0.6).unwrap();
    v.push(small.photon_probability(&[2, 0, 1, 1, 0]).unwrap());
    for s in small.sample_clicks(300, 9) {
        v.push(s.iter().enumerate().map(|(i, &b)| if b { (1u32 << i) as f64 } else { 0.0 }).sum());
    }
    let (p1, p2) = small.click_marginals();
    v.extend(p1);
    v.extend(p2);
    assert_eq!(digest(&v), 0x039d_e851_06c3_567d);
}

/// A 16-qubit random circuit at fusion widths 0, 2 and 3: every kernel path
/// (block-parallel, split by the top bit, split into quarters, general fused,
/// diagonal), threaded natively and sequential under WASI, and the chunked
/// reductions (norm, expectation, samples).
#[test]
fn dense_state_vector() {
    use gates::Gate;
    use quantum_sv::StateVector;
    let n = 16u32;
    let mut g = Vec::new();
    let mut a = 0.1f64;
    for layer in 0..6u32 {
        for q in 0..n {
            a += 0.731;
            g.push(match (q + layer) % 5 {
                0 => Gate::rx(q, a),
                1 => Gate::ry(q, a),
                2 => Gate::rz(q, a),
                3 => Gate::sqrt_w(q),
                _ => Gate::t(q),
            });
        }
        let mut q = layer % 2;
        while q + 1 < n {
            g.push(if q % 3 == 0 { Gate::fsim(q, q + 1, a, 0.5 * a) } else { Gate::cz(q, q + 1) });
            q += 2;
        }
        g.push(Gate::iswap(0, n - 1));
    }
    let mut v = Vec::new();
    for width in [0usize, 2, 3] {
        let mut sv = StateVector::zero(n).unwrap();
        sv.run(&g, width).unwrap();
        v.extend(sv.re.iter().step_by(97).copied());
        v.extend(sv.im.iter().step_by(89).copied());
        v.push(sv.norm2());
        v.push(sv.expectation(&[(0.5, vec![(0, 'Z'), (7, 'X')]), (1.5, vec![(3, 'Y'), (15, 'Y')])]));
        v.extend(sv.sample(64, width as u64).into_iter().map(|x| x as f64));
    }
    assert_eq!(digest(&v), 0xd0df_e437_5b28_903f);
}


/// Storage layouts under every fit set, cultivation interpolation (in-crate
/// `ln` and `exp`), the threshold-to-storage conversion, and the whole factoring
/// machine at three key sizes.
#[test]
fn architecture_estimates() {
    use quantum_arch::*;
    let mut v = Vec::new();
    for model in [StorageModel::si1000(), StorageModel::uniform(), StorageModel::uniform_refit()] {
        for target in [1e-9, 3.7e-13, 1e-15] {
            for b in [
                best_block(&model, Storage::Cold, Yoke::Rows, target, 250),
                best_block(&model, Storage::Hot, Yoke::Rows, target, 250),
                grid_block(&model, 16, target),
            ]
            .into_iter()
            .flatten()
            {
                v.push(f64::from(b.distance));
                v.push(b.qubits_per_logical());
                v.push(b.error_per_logical_round);
            }
        }
    }
    for (p, t) in [(5e-4, 3.3e-9), (1e-3, 1e-7), (1e-3, 4.4e-6), (2e-3, 2.2e-5)] {
        let c = cultivate(p, t).unwrap();
        v.push(c.volume);
        v.push(c.error);
    }
    let f = StorageFit::flat_from_threshold(0.023, 0.0100, 1e-3);
    v.push(f.suppression);
    v.push(f.per_patch_round(17, 1, 1));
    for bits in [1024, 2048, 4096] {
        let e = factoring_estimate(&FactoringParams::published(bits).unwrap(), &Machine::published()).unwrap();
        v.push(e.physical_qubits() as f64);
        v.push(e.ccz_rounds);
        v.push(e.hours_per_shot);
        v.push(e.hours_per_shot_rounded);
        v.push(e.shot_survival);
        v.push(e.days);
        v.push(e.toffolis_per_shot);
    }
    assert_eq!(digest(&v), 0xd166_b4c6_c731_7a52);
}

/// A bivariate-bicycle memory: the circuit's error model (probabilities from
/// the independent-channel conversion), BP over it with in-crate logarithms,
/// and OSD on the shots BP leaves unexplained.
#[test]
fn qldpc_memory_decoding() {
    use quantum_frame::{error_model, parse, sample, Basis};
    use quantum_ldpc::*;
    let c = parse(&memory_circuit(&BbCode::bb72(), 3, 0.005, Basis::Z)).unwrap();
    let model = error_model(&c).unwrap();
    let dec = BpOsd::new(&model, BpOsdConfig { max_iter: 40, ..BpOsdConfig::default() });
    let det = sample(&c, 192, 11);
    let mut v: Vec<f64> = model.mechanisms.iter().step_by(37).map(|m| m.probability).collect();
    let mut osd = 0;
    for s in 0..192 {
        let d = dec.decode(&det.fired(s));
        osd += u32::from(!d.converged);
        v.push(d.observables as f64);
        v.push(f64::from(d.iterations));
        v.push(det.flips(s) as f64);
    }
    assert!(osd > 0, "the sample must exercise OSD");
    assert_eq!(digest(&v), 0x5479_ebe5_f648_b739);
}

/// Optimisation: every simulated-bifurcation variant (batched, with fields
/// riding on an extra spin), annealing with in-crate `exp`, single-flip
/// descent, and the exact referee.
#[test]
fn quantum_inspired_optimisation() {
    use quantum_qio::*;
    let p = Ising::sk(40, 3);
    let mut v = Vec::new();
    for var in [Variant::Ballistic, Variant::Discrete, Variant::HeatedBallistic(0.0), Variant::HeatedDiscrete(0.0)] {
        for t in simulated_bifurcation(&p, &SbConfig::published(var, 300, 11, 9), 3).trials {
            v.push(t.energy);
            v.extend(t.spins.iter().map(|&s| f64::from(s)));
        }
    }
    let q = Ising::new(9, &[(0, 4, 0.75), (4, 8, -1.5), (1, 2, 1.0), (2, 6, 0.25)], &[0.5, 0.0, -0.3, 0.0, 0.1, 0.0, 0.0, 0.0, -0.9]).unwrap();
    for t in simulated_bifurcation(&q, &SbConfig::published(Variant::Discrete, 200, 5, 2), 2).trials {
        v.push(t.energy);
    }
    for t in simulated_annealing(&p, &SaConfig { sweeps: 200, beta0: 0.1, beta1: 2.5, trials: 6, seed: 4 }, 3).trials {
        v.push(t.energy);
        v.push(descend(&p, &t.spins).energy);
    }
    let g = ground_state(&Ising::sk(18, 5), 4).unwrap();
    v.push(g.energy);
    v.extend(g.spins.iter().map(|&s| f64::from(s)));
    assert_eq!(digest(&v), 0xc01a_7afa_39c0_c487);
}

/// Union-find decoding: event-driven growth on the radix heap, merges and
/// peeling, over a sampled surface-code memory, with a reused workspace.
#[test]
fn union_find_decoding() {
    use quantum_frame::{error_model, parse, sample, surface_code_memory};
    let c = parse(&surface_code_memory(7, 7, 0.006)).unwrap();
    let g = quantum_match::MatchingGraph::from_model(&error_model(&c).unwrap()).unwrap();
    let uf = quantum_uf::UnionFind::new(&g);
    let mut ws = uf.workspace();
    let det = sample(&c, 512, 21);
    let v: Vec<f64> = (0..512).map(|s| uf.decode(&det.fired(s), &mut ws).unwrap() as f64).collect();
    assert!(v.iter().any(|&x| x != 0.0), "the sample must include flips");
    assert_eq!(digest(&v), 0xe43c_2ad4_9244_3118);
}

/// Noisy circuits: trajectories under every channel (Pauli branches, damping
/// jumps with renormalisation, readout flips), the exact density matrix and
/// both cross-entropy figures.
#[test]
fn noisy_circuits() {
    use gates::Gate;
    use quantum_noisy::*;
    let n = 6u32;
    let mut g = Vec::new();
    let mut a = 0.3f64;
    for layer in 0..5u32 {
        for q in 0..n {
            a += 0.917;
            g.push(if (q + layer) % 2 == 0 { Gate::rx(q, a) } else { Gate::ry(q, a) });
        }
        let mut q = layer % 2;
        while q + 1 < n {
            g.push(Gate::cz(q, q + 1));
            q += 2;
        }
    }
    let noise = NoiseModel { depol1: 0.01, depol2: 0.03, damping: 0.02, dephasing: 0.01, readout: 0.02 };
    let obs = vec![vec![(1.0, vec![(0, 'Z')])], vec![(1.0, vec![(1, 'X'), (4, 'Y')])]];
    let run = trajectories(&g, n, &noise, &obs, 40, 6, 3).unwrap();
    let mut v: Vec<f64> = run.expectations.iter().flat_map(|e| [e.mean, e.stderr]).collect();
    v.extend(run.samples.iter().map(|&x| x as f64));
    let rho = DensityMatrix::run(&g, n, &noise).unwrap();
    v.extend(obs.iter().map(|o| rho.expectation(o)));
    v.push(rho.trace());
    let mut ideal = quantum_sv::StateVector::zero(n).unwrap();
    ideal.run(&g, 0).unwrap();
    v.push(exact_xeb(&ideal, &rho, noise.readout));
    let x = linear_xeb(&ideal, &run.samples);
    v.push(x.mean);
    v.push(x.stderr);
    v.push(digital_fidelity(&g, n, &noise));
    assert_eq!(digest(&v), 0x78b2_ddfb_e7c3_6bea);
}

/// DMRG: the symmetric eigensolver, Lanczos with re-orthogonalisation,
/// density-matrix truncation and entropies (in-crate `ln`), on two chains;
/// the free-fermion Ising energy; and exact diagonalisation.
#[test]
fn dmrg_ground_states() {
    use quantum_dmrg::*;
    let mut v = Vec::new();
    let cfg = DmrgConfig { max_bond: 16, sweeps: 2, ..DmrgConfig::default() };
    for chain in [Chain::heisenberg(10, 1.0, 0.8, 0.1), Chain::ising(12, 1.0, 0.9)] {
        let r = dmrg(&chain, &cfg);
        v.extend(r.sweep_energies);
        v.push(r.discarded);
        v.extend(r.entropies);
        v.push(exact_ground_energy(&chain).unwrap());
    }
    v.push(ising_exact_energy(40, 1.0, 0.7));
    let a: Vec<f64> = (0..36).map(|k| ((k * 7 % 11) as f64 - 5.0) / (1.0 + (k % 6) as f64)).collect();
    let sym: Vec<f64> = (0..36).map(|k| a[k] + a[(k % 6) * 6 + k / 6]).collect();
    let (vals, vecs) = symmetric_eigen(&sym, 6);
    v.extend(vals);
    v.extend(vecs.iter().map(|x| x.abs()));
    assert_eq!(digest(&v), 0xcdd9_d9fb_5ba6_d71c);
}

#[test]
fn tdvp_time_evolution() {
    use quantum_dmrg::{Chain, Mps, Op, Scalar, C};
    use quantum_tdvp::*;
    let mut v = Vec::new();
    // A truncated Heisenberg quench from the Néel state, and the Ising quench
    // with its free-fermion referee.
    let n = 10;
    let neel = Mps::basis(&(0..n).map(|i| i % 2 == 1).collect::<Vec<_>>());
    let mut ev = Evolution::new(&Chain::heisenberg(n, 1.0, 0.7, 0.1), &neel, TdvpConfig { dt: 0.1, max_bond: 6, threads: 1, ..TdvpConfig::default() });
    for _ in 0..6 {
        v.push(ev.step().discarded);
    }
    v.extend(ev.expect(Op::z()));
    v.extend(ev.bond_expect(Op::sz(), Op::sz()));
    v.extend(ev.entropies());
    v.push(ev.energy());
    let o = ev.overlap(&neel);
    v.extend([o.re, o.im]);
    let plus = Mps::product(&vec![[C::ONE, C::ONE]; 8]);
    let mut ev = Evolution::new(&Chain::ising(8, 1.0, 0.9), &plus, TdvpConfig { dt: 0.2, max_bond: 16, cutoff: 0.0, threads: 1, ..TdvpConfig::default() });
    for _ in 0..3 {
        ev.step();
    }
    v.extend(ev.expect(Op::x()));
    let ff = ising_quench(8, 1.0, 0.9, ev.time());
    v.extend(ff.x);
    v.extend(ff.zz);
    v.push(ff.echo);
    assert_eq!(digest(&v), 0xe1b4_fd28_af23_50cd);
}

#[test]
fn chemistry_by_dmrg() {
    use quantum_chem::*;
    let mut v = Vec::new();
    let fd = read_fcidump(include_str!("data/lih.fcidump")).unwrap();
    v.push(fci(&fd).unwrap().energy);
    v.push(hartree_fock_energy(&fd));
    let r = ground_state(&fd, &ChemConfig { max_bond: 10, sweeps: 1, threads: 1, ..ChemConfig::default() });
    v.extend(r.sweep_energies);
    v.push(r.discarded);
    v.extend(r.entropies);
    let mpo = molecular_mpo(&fd, &ChemConfig::default());
    v.extend(mpo.dims.iter().map(|&d| d as f64));
    assert_eq!(digest(&v), 0xe955_a4ff_7eda_9960);
}

#[test]
fn neural_quantum_states() {
    use quantum_nqs::*;
    let mut v = Vec::new();
    let n = 8;
    let model = SpinModel::heisenberg(n, &chain_edges(n, false), 1.0);
    v.push(model.exact_ground_energy(Some(0)).unwrap());
    let rbm = Rbm::new(n, 2, 0.05, 1).with_marshall_sign(checkerboard(n, 1));
    // Exact expectations, then sampled ones, a few steps each.
    let mut exact = Vmc::new(model.clone(), rbm.clone(), VmcConfig { exact: true, magnetisation: Some(0), lr: 0.1, threads: 1, ..VmcConfig::default() });
    for _ in 0..5 {
        let e = exact.step();
        v.extend([e.energy, e.variance]);
    }
    let mut sampled = Vmc::new(SpinModel::ising(n, &chain_edges(n, false), 1.0, 0.8), Rbm::new(n, 1, 0.1, 2), VmcConfig { samples: 200, chains: 4, threads: 1, ..VmcConfig::default() });
    for _ in 0..3 {
        let e = sampled.step();
        v.extend([e.energy, e.variance, e.error, e.acceptance]);
    }
    v.extend(&sampled.rbm.w);
    // A complex machine on the Majumdar–Ghosh ring, which has no sign rule:
    // complex exponentials, logarithms and arctangents in every step.
    let mg = SpinModel::j1j2(n, &chain_edges(n, true), &chain_next_nearest(n, true), 1.0, 0.5);
    let mut cplx = Vmc::new(mg, Rbm::complex(n, 2, 0.05, 3), VmcConfig { exact: true, magnetisation: Some(0), lr: 0.1, threads: 1, ..VmcConfig::default() });
    for _ in 0..4 {
        let e = cplx.step();
        v.extend([e.energy, e.variance]);
    }
    v.extend(cplx.rbm.w.iter().flat_map(|c| [c.re, c.im]));
    // A projected machine on a 10-site ring (momentum π, odd), started near
    // even: the metric goes singular and the eigendecomposition takes over.
    let ring = SpinModel::heisenberg(10, &chain_edges(10, true), 1.0);
    let rbm = Rbm::complex(10, 1, 0.01, 7).with_symmetry(Symmetry::ring(10, true).with_spin_flip(true));
    let mut proj = Vmc::new(ring, rbm, VmcConfig { exact: true, magnetisation: Some(0), lr: 0.1, threads: 1, ..VmcConfig::default() });
    for k in 0..60 {
        let e = proj.step();
        if k % 10 == 9 {
            v.extend([e.energy, e.variance]);
        }
    }
    v.extend(proj.rbm.b.iter().flat_map(|c| [c.re, c.im]));
    assert_eq!(digest(&v), 0x834c_988e_b34e_814f);
}