open-standards/wai/wai-quantum/tests/bit_identity.rs
Transaction Science a84f5d704f wai-quantum 0.6.12: cheaper contraction trees, from V-cycles and tree annealing
Each bisection now runs two V-cycles after the first projection. The graph is
coarsened again with every match inside one side, so the split carries to the coarse
graph with the same cut. It is refined there and on the way back down, which never
raises the cut.

The best few trees are annealed after reconfiguration (quantum_tn::anneal). Random
rotations at each node are accepted on the change in log total cost at a falling
temperature; the cheapest tree seen is kept, then reconfigured again. A rotation keeps
the node's index set, so it changes the cost of two contractions only.

On the 61-qubit Nighthawk logical circuits (split, reduced, 1024 trials, 8 restarts),
the median of 4 seeds, in log10 complex multiply-adds for one amplitude:

  depth   before   after
  8       5.30     5.29
  12      8.10     8.09
  16      11.12    10.97
  20      14.02    13.68

At depth 20 the search takes about 1.7 times as long.

Also tried, and left out: 24 initial partitions per bisection (no gain at depth 16 at
equal time), a coarsening rating that keeps coarse nodes balanced (mixed), more
refinement passes (no change), and 200 annealing sweeps (no better than 50).

The tensor-network plan pinned in tests/bit_identity.rs changes, by design. The native
and WASI builds compute the same new digest. A test checks that annealing never costs
more, keeps every tensor contracted once and after it exists, repeats its tree for a
seed, and gives the dense amplitude.

The changelog also records what 0.6.11 gained after its entry was written: the
circuit's layout as data, path-search restarts and the command line's release builds.

Co-Authored-By: Claude Opus 5.5 <[email protected]>
(cherry picked from commit 5c253c907402ce4a8e6140a676e3e58d86e01398)
2026-10-20 17:34:17 -04:00

582 lines
23 KiB
Rust
Raw Permalink Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

//! 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);
}