From 92f9e000cde8b50d35891ed33055e32a98025c2e Mon Sep 17 00:00:00 2001 From: Transaction Science Date: Mon, 5 Oct 2026 05:27:49 -0400 Subject: [PATCH] wai-quantum 0.3.31: quantum-inspired optimisation, a re-runnable classical baseline Optimisation is the claim made most often for quantum hardware, and each such claim is a comparison with a classical solver. quantum_qio is that solver, deterministic and re-checkable: - Ising problems (sparse symmetric J, fields h) and Max-Cut (J = -w), a G-set reader, and the +-1 all-to-all spin glass. - Simulated bifurcation from its published update rules: symplectic Euler, perfectly inelastic walls, ballistic (f = Jx) and discrete (f = J sgn x) forces, and both heated (arXiv:2203.08361), with c0 = c1 / (rms(J) sqrt N). Fields ride on an extra spin fixed at +1. - Eight trials advance per pass over the couplings; each trial's arithmetic is unchanged, so a batched run leaves the same final positions bit for bit as one trial at a time. The test fails on a 1e-15 change to one lane. - Simulated annealing (single-spin Metropolis, linear beta), single-flip descent, an exact Gray-code referee up to 34 spins (threaded), and step-to-solution S = N_s ln 0.01 / ln(1 - P). Checked: - every variant and annealing reach the exact ground state on 20-spin glasses; - G-set best-known cuts reached (seed 1): G1 11624 (39% of 2,000-step trials), G2 11620, G3 11622, G11 564, G14 3064 (annealing; bifurcation stops at 3062), G22 13359, G43 6660; G55 reaches 10292 of 10299; - the published dense-glass step (1.1) is unstable on sparse Max-Cut graphs: the stiffest mode has w^2 = a0(a0 - c0 mu_min) and mu_min = -degree, so G1's limit is near 0.88; documented, with the settings that work. CLI: wai-quantum qio (--gset FILE | --sk N) --method dsb|bsb|hdsb|hbsb|sa|exact [--descend --target --spins --json]. Example: qio_gset. Pinned bit for bit natively and under wasm32-wasip2 (16 pins); CLI output byte-identical across the two, threaded or not, G-set files included. Co-Authored-By: Claude Opus 5.5 --- wai/wai-quantum-serve/Cargo.lock | 2 +- wai/wai-quantum/Cargo.lock | 2 +- wai/wai-quantum/Cargo.toml | 11 +- wai/wai-quantum/README.md | 23 + wai/wai-quantum/examples/qio_gset.rs | 74 +++ wai/wai-quantum/src/lib.rs | 6 +- wai/wai-quantum/src/quantum_qio.rs | 895 ++++++++++++++++++++++++++ wai/wai-quantum/src/repro.rs | 28 +- wai/wai-quantum/tests/bit_identity.rs | 31 +- wai/wai-rs/Cargo.lock | 2 +- wai/wai-rs/Cargo.toml | 3 + wai/wai-rs/src/bin/wai_quantum.rs | 79 +++ wai/wai-rs/src/lib.rs | 3 + wai/wai-web/Cargo.lock | 2 +- 14 files changed, 1139 insertions(+), 22 deletions(-) create mode 100644 wai/wai-quantum/examples/qio_gset.rs create mode 100644 wai/wai-quantum/src/quantum_qio.rs diff --git a/wai/wai-quantum-serve/Cargo.lock b/wai/wai-quantum-serve/Cargo.lock index ad69b461e..72b29a56e 100644 --- a/wai/wai-quantum-serve/Cargo.lock +++ b/wai/wai-quantum-serve/Cargo.lock @@ -589,7 +589,7 @@ dependencies = [ [[package]] name = "wai-quantum" -version = "0.3.30" +version = "0.3.31" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum/Cargo.lock b/wai/wai-quantum/Cargo.lock index 5ebca3ae9..47e0fe263 100644 --- a/wai/wai-quantum/Cargo.lock +++ b/wai/wai-quantum/Cargo.lock @@ -415,7 +415,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" [[package]] name = "wai-quantum" -version = "0.3.30" +version = "0.3.31" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum/Cargo.toml b/wai/wai-quantum/Cargo.toml index e41d71ae7..80bf17974 100644 --- a/wai/wai-quantum/Cargo.toml +++ b/wai/wai-quantum/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "wai-quantum" -version = "0.3.30" +version = "0.3.31" edition = "2024" rust-version = "1.97" license = "Apache-2.0" @@ -71,6 +71,9 @@ quantum_arch = [] # qLDPC memory at circuit level: bivariate-bicycle codes, their depth-8 syndrome # cycle, and BP-OSD decoding of the circuit's error model. quantum_ldpc = ["quantum_frame"] +# Quantum-inspired optimisation: simulated bifurcation (ballistic, discrete, +# heated), simulated annealing and an exact referee, on Ising and Max-Cut. +quantum_qio = [] # Mixed states and the observables that distinguish quantum theory from a classical # account: density matrices, the partial trace, exact Pauli expectations, CHSH. quantum_info = ["quantum"] @@ -110,13 +113,17 @@ quantum_phasor_meter = ["_signing"] # Everything. full = [ - "quantum", "quantum_toolchain", "quantum_mps", "quantum_stabilizer", "quantum_pauli", "quantum_spd", "quantum_spd_receipt", "quantum_bptn", "quantum_resource", "quantum_resource_receipt", "quantum_frame", "quantum_tn", "quantum_gbs", "quantum_sv", "quantum_arch", "quantum_ldpc", "quantum_qir", "quantum_info", "quantum_qudit", "quantum_comm", "quantum_source", "quantum_channel", + "quantum", "quantum_toolchain", "quantum_mps", "quantum_stabilizer", "quantum_pauli", "quantum_spd", "quantum_spd_receipt", "quantum_bptn", "quantum_resource", "quantum_resource_receipt", "quantum_frame", "quantum_tn", "quantum_gbs", "quantum_sv", "quantum_arch", "quantum_ldpc", "quantum_qio", "quantum_qir", "quantum_info", "quantum_qudit", "quantum_comm", "quantum_source", "quantum_channel", "quantum_receipt", "quantum_ops", "quantum_cal", "quantum_mitigate", "quantum_qec", "quantum_control", "quantum_noise", "quantum_compile", "quantum_atom", "quantum_qbom", "quantum_vml", "quantum_phasor", "quantum_qfhrr", "quantum_kernel", "quantum_qdata", "quantum_phasor_meter", ] +[[example]] +name = "qio_gset" +required-features = ["quantum_qio"] + [[example]] name = "ldpc_memory" required-features = ["quantum_ldpc"] diff --git a/wai/wai-quantum/README.md b/wai/wai-quantum/README.md index 47fbb445a..cb9eaa45f 100644 --- a/wai/wai-quantum/README.md +++ b/wai/wai-quantum/README.md @@ -64,6 +64,7 @@ never creates information. | `quantum_bptn` | **tensor-network states on the hardware graph**: simple-update evolution with belief-propagation environments and re-gauging, the device's couplings applied a layer at a time in parallel (bit-identical to one at a time). Exact on trees; on the 127-qubit lattice it gives the exactly computed 5-step magnetization to 4e-7 at bond dimension 32 — an independent check on sparse Pauli dynamics, whose errors it does not share | | `quantum_resource` | **fault-tolerant resource estimation**: from a program's logical counts (qubits, T, Toffoli and rotation gates, measurements) to code distance, T-state factories found by exhaustive search, physical qubits and runtime, on six reference qubits and three codes, with the qubits-versus-time frontier. Reproduces the published worked examples to their printed precision | | `quantum_resource_receipt` | a signed, energy-accounted receipt over an estimate, re-checkable bit for bit on any machine | +| `quantum_qio` | **quantum-inspired optimisation**: Ising and Max-Cut problems (G-set reader, spin-glass generator), simulated bifurcation in its ballistic, discrete and heated forms (batched eight trials per pass, bit-identical to one at a time), simulated annealing, single-flip descent, an exact Gray-code referee up to 34 spins, and step-to-solution. Reaches the best-known cut on seven G-set instances | | `quantum_ldpc` | **qLDPC memory at circuit level**: bivariate-bicycle codes from their two polynomials (`k` by rank, logical bases), the depth-8 syndrome cycle with circuit-level noise, its detector error model, and BP-OSD decoding (min-sum with adaptive scaling, ordered statistics with a combination sweep). The error model has the same columns as the source's own decoding matrices, and on identical syndromes the decoder makes the source's decisions shot for shot | | `quantum_arch` | **architecture-level estimation**: yoked surface-code storage (row and grid parity checks over patches, hot and cold), magic-state cultivation read off the published simulation curves, CCZ factories, and a clock set by whichever is slowest of lattice surgery, the magic-state supply and the control loop. Assembles whole-machine estimates; its storage search matches the source's footprint tool layout for layout | | `quantum_receipt` | signed, energy-accounted receipt over a reconstruction | @@ -154,6 +155,28 @@ search limits, to the bit). Re-running it anywhere must reproduce every figure, down to each factory round, and the bytes are the same natively and as a WASI component. +## Optimisation + +```sh +wai-quantum qio --gset G22 --method dsb --dt 0.5 --steps 20000 --trials 64 --target 13359 +wai-quantum qio --gset G14 --method sa --steps 300000 --beta0 0.2 --beta1 4 --target 3064 +wai-quantum qio --sk 24 --method exact --spins # the true ground state, for small problems +``` + +Optimisation is the claim made most often for quantum hardware, and every +such claim is a comparison against a classical solver. `quantum_qio` is a +re-runnable one. It offers: +- simulated bifurcation (ballistic, discrete, and both with heating); +- simulated annealing; +- single-flip descent; +- an exact referee that settles small instances. + +A result is a spin configuration, so anyone can re-check its cut or energy +in one pass. On the G-set Max-Cut benchmark the best-known cut is reached on +G1, G2, G3, G11, G14, G22 and G43. On G55 the best found is 10,292 against +10,299. The module documents each setting and success rate, and the bound +that sparse graphs put on the time step. + ## qLDPC memory ```sh diff --git a/wai/wai-quantum/examples/qio_gset.rs b/wai/wai-quantum/examples/qio_gset.rs new file mode 100644 index 000000000..bcc4a06a9 --- /dev/null +++ b/wai/wai-quantum/examples/qio_gset.rs @@ -0,0 +1,74 @@ +//! Solve a G-set Max-Cut instance with simulated bifurcation and report the +//! best cut, how often the target is reached, and the steps to reach it with +//! 99% probability. +//! +//! ```text +//! cargo run --release --example qio_gset --features quantum_qio -- FILE VARIANT STEPS TRIALS [TARGET] [DT] [C1] [GAMMA] [THREADS] +//! ``` +//! VARIANT is bsb, dsb, hbsb or hdsb. Omitted settings take the published values. +use wai_quantum::quantum_qio::{descend, read_gset, simulated_annealing, simulated_bifurcation, step_to_solution, SaConfig, SbConfig, Variant}; + +fn main() { + let a: Vec = std::env::args().collect(); + let g = read_gset(&std::fs::read_to_string(&a[1]).expect("read")).expect("parse"); + let polish = a[2].ends_with("+ls"); + let name = a[2].trim_end_matches("+ls"); + if name == "sa" { + // STEPS sweeps; DT and C1 slots carry beta0 and beta1. + let sweeps: u32 = a[3].parse().unwrap(); + let trials: usize = a[4].parse().unwrap(); + let target: Option = a.get(5).and_then(|x| x.parse().ok()); + let beta0 = a.get(6).and_then(|x| x.parse().ok()).unwrap_or(0.1); + let beta1 = a.get(7).and_then(|x| x.parse().ok()).unwrap_or(3.0); + let threads: usize = a.get(9).map_or(8, |x| x.parse().unwrap()); + let t0 = std::time::Instant::now(); + let run = simulated_annealing(&g.ising(), &SaConfig { sweeps, beta0, beta1, trials, seed: 1 }, threads); + let cuts: Vec = run.trials.iter().map(|t| g.cut(&t.spins)).collect(); + let best = cuts.iter().cloned().fold(f64::MIN, f64::max); + print!("SA beta {beta0}->{beta1} sweeps={sweeps} trials={trials} best {best} mean {:.1}", cuts.iter().sum::() / cuts.len() as f64); + if let Some(t) = target { + let p = cuts.iter().filter(|&&c| c >= t).count() as f64 / cuts.len() as f64; + print!(" P(≥{t}) {p:.4} S99 {:?}", step_to_solution(p, f64::from(sweeps)).map(|s| s.round())); + } + println!(" {:.1}s", t0.elapsed().as_secs_f64()); + return; + } + let v = match name { + "bsb" => Variant::Ballistic, + "dsb" => Variant::Discrete, + "hbsb" => Variant::HeatedBallistic(0.0), + "hdsb" => Variant::HeatedDiscrete(0.0), + other => panic!("variant {other}"), + }; + let steps: u32 = a[3].parse().unwrap(); + let trials: usize = a[4].parse().unwrap(); + let mut cfg = SbConfig::published(v, steps, trials, 1); + let target: Option = a.get(5).and_then(|x| x.parse().ok()); + if let Some(dt) = a.get(6).and_then(|x| x.parse().ok()) { + cfg.dt = dt; + } + if let Some(c1) = a.get(7).and_then(|x| x.parse().ok()) { + cfg.c1 = c1; + } + if let Some(gamma) = a.get(8).and_then(|x| x.parse::().ok()) { + cfg.variant = match cfg.variant { + Variant::HeatedBallistic(_) => Variant::HeatedBallistic(gamma), + Variant::HeatedDiscrete(_) => Variant::HeatedDiscrete(gamma), + other => other, + }; + } + let threads: usize = a.get(9).map_or(8, |x| x.parse().unwrap()); + let t0 = std::time::Instant::now(); + let run = simulated_bifurcation(&g.ising(), &cfg, threads); + let dt = t0.elapsed().as_secs_f64(); + let ising = g.ising(); + let cuts: Vec = run.trials.iter().map(|t| if polish { g.cut(&descend(&ising, &t.spins).spins) } else { g.cut(&t.spins) }).collect(); + let best = cuts.iter().cloned().fold(f64::MIN, f64::max); + let mean = cuts.iter().sum::() / cuts.len() as f64; + print!("n={} edges={} {}{:?} dt={} c1={} steps={} trials={} best {} mean {:.1}", g.n, g.edges.len(), if polish { "+descent " } else { "" }, cfg.variant, cfg.dt, cfg.c1, steps, trials, best, mean); + if let Some(t) = target { + let p = cuts.iter().filter(|&&c| c >= t).count() as f64 / cuts.len() as f64; + print!(" P(≥{t}) {p:.4} S99 {:?}", step_to_solution(p, f64::from(steps)).map(|s| s.round())); + } + println!(" {:.1}s ({:.2} ms/trial-thread)", dt, 1e3 * dt * threads as f64 / trials as f64); +} diff --git a/wai/wai-quantum/src/lib.rs b/wai/wai-quantum/src/lib.rs index da1ff3b1a..85a7d16e4 100644 --- a/wai/wai-quantum/src/lib.rs +++ b/wai/wai-quantum/src/lib.rs @@ -36,7 +36,8 @@ pub mod gates; feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_sv", - feature = "quantum_arch" + feature = "quantum_arch", + feature = "quantum_qio" ))] mod repro; @@ -109,6 +110,9 @@ pub mod quantum_frame; #[cfg(feature = "quantum_ldpc")] pub mod quantum_ldpc; +#[cfg(feature = "quantum_qio")] +pub mod quantum_qio; + #[cfg(feature = "quantum_frame")] pub mod quantum_match; diff --git a/wai/wai-quantum/src/quantum_qio.rs b/wai/wai-quantum/src/quantum_qio.rs new file mode 100644 index 000000000..ac10828d7 --- /dev/null +++ b/wai/wai-quantum/src/quantum_qio.rs @@ -0,0 +1,895 @@ +//! Quantum-inspired OPTIMISATION — `wai.quantum.qio`. +//! +//! The most common claim made for quantum hardware is in optimisation: an +//! annealer or an Ising machine finds a lower-energy spin configuration, or +//! finds one faster. Every such claim is a comparison with a classical +//! baseline, and the baseline has to be as good as the field's best. +//! Otherwise the comparison measures the baseline and not the machine. +//! +//! This module is that baseline, built to be re-run: +//! +//! - **Problems.** [`Ising`] holds couplings `J` (sparse, symmetric) and +//! fields `h`; its energy is `E = −½ Σ J_ij s_i s_j − Σ h_i s_i`. Max-Cut +//! with weights `w` maps to it with `J = −w`, and [`read_gset`] reads the +//! standard benchmark format. [`Ising::sk`] builds the all-to-all ±1 spin +//! glass used to benchmark Ising machines. +//! - **Simulated bifurcation.** [`simulated_bifurcation`] integrates the +//! Hamiltonian dynamics `ẋ = a₀y`, `ẏ = −(a₀ − a(t))x + c₀f` with the +//! symplectic Euler method: +//! - `a(t)` rises linearly from 0 to `a₀`; +//! - perfectly inelastic walls hold each `|x_i| ≤ 1`; +//! - the force is `f = Jx` (ballistic) or `f = J sgn(x)` (discrete); +//! - the heated variants add `γy` after the walls (arXiv:2203.08361); +//! - `c₀ = c₁ / (σ_J √N)`, with `σ_J` the root-mean-square coupling. +//! +//! Many trials run at once, threaded. Each trial is seeded on its own and +//! uses only `+ − × ÷` and `√`, so a run is the same bits on every machine +//! and at any thread count. +//! - **Simulated annealing.** [`simulated_annealing`] is single-spin Metropolis +//! on a linear inverse-temperature schedule, the classical reference that +//! annealing hardware is compared against. +//! - **An exact referee.** [`ground_state`] enumerates every configuration in +//! Gray-code order up to 34 spins, so the heuristics can be held to the true +//! optimum on small instances. +//! - **Certificates.** A result is a spin configuration. Its energy or cut is +//! re-checked by anyone in one pass over the couplings +//! ([`Ising::energy`], [`MaxCut::cut`]). +//! - **Step-to-solution.** [`step_to_solution`] is the effort needed to reach +//! a target with 99% probability, `S = N_s · ln 0.01 / ln(1 − P)`. +//! +//! # Checked +//! +//! - **Small instances.** Every variant of simulated bifurcation, and +//! annealing, reaches the exact ground state from enumeration on 20-spin +//! spin glasses. +//! - **Batching and threads.** The batched kernel leaves the same final +//! positions, bit for bit, as one trial at a time, at any thread count. The +//! test that says so fails on a `10⁻¹⁵` change to one trial's arithmetic. +//! - **The G-set Max-Cut benchmark.** The best-known cut is reached on seven +//! instances. Each figure is the share of trials reaching it, seed 1, from +//! the `qio_gset` example: +//! +//! | instance | best known | reached by | trials reaching it | +//! |----------|-----------:|-------------------------------------|-------------------:| +//! | G1 | 11,624 | heated ballistic, Δt 0.8, 2,000 steps | 39% of 64 | +//! | G2 | 11,620 | discrete, Δt 0.5, 20,000 steps | 5.5% of 128 | +//! | G3 | 11,622 | heated ballistic, Δt 0.8, 5,000 steps | 16% of 64 | +//! | G11 | 564 | heated ballistic, Δt 0.8, 5,000 steps | 3.1% of 64 | +//! | G14 | 3,064 | annealing, 300,000 sweeps, β 0.2→4 | 1 of 64 | +//! | G22 | 13,359 | discrete, Δt 0.5, 20,000 steps | 3.1% of 64 | +//! | G43 | 6,660 | heated ballistic, Δt 0.8, 5,000 steps | 4.7% of 64 | +//! +//! On G14 simulated bifurcation stops at 3,062. On G55 (best known 10,299) +//! the best found is 10,292 by annealing and 10,290 by bifurcation. +//! +//! # The time step is bounded by the graph +//! +//! The symplectic step is stable only while `ω·Δt < 2`, where the stiffest +//! frequency is `ω² = a₀(a₀ − c₀ μ_min)` and `μ_min` is the most negative +//! eigenvalue of `J`. For a dense ±1 spin glass, `c₀μ_min ≈ −2c₁`, so the +//! published `Δt = 1.1` is safe. +//! +//! On a regular-ish sparse graph with Max-Cut couplings, the uniform mode +//! dominates: `μ_min = −degree`, and `c₀μ_min` is several times larger. On +//! G1 that puts the limit near `Δt ≈ 0.88`. At the published 1.1, discrete +//! and heated runs blow up the uniform mode and cut nothing. +//! [`SbConfig::published`] is the dense setting; sparse graphs need a smaller +//! step. +//! +//! # Honest boundaries +//! +//! - **Heuristics.** Neither simulated bifurcation nor annealing proves +//! optimality; a best-known value is the best anyone has found. Only +//! [`ground_state`] is exact, and only up to 34 spins. +//! - **Parameters matter.** The time step, `c₁` and the heating rate are +//! tuned per problem class. The defaults are the published ones for dense +//! spin glasses; sparse graphs want a smaller step (above). +//! - **Not every instance falls.** G55's best-known cut was not reached here. + +use crate::repro::exp; + +// --------------------------------------------------------------------------- +// Problems +// --------------------------------------------------------------------------- + +/// An Ising problem: minimise `E = −½ Σ_ij J_ij s_i s_j − Σ_i h_i s_i` over +/// `s ∈ {−1, +1}ⁿ`. +#[derive(Clone, Debug, PartialEq)] +pub struct Ising { + n: usize, + /// Symmetric couplings in compressed rows: row `i` holds `(j, J_ij)`, + /// both directions stored, no diagonal. + start: Vec, + col: Vec, + val: Vec, + h: Vec, +} + +/// Why a problem could not be built or read. +#[derive(Clone, Debug, PartialEq)] +pub enum QioError { + /// A coupling names a spin outside `0..n`, or couples a spin to itself. + BadCoupling { i: usize, j: usize }, + /// Text that is not a graph in the benchmark format, with its line. + Parse { line: usize, message: String }, + /// Too many spins for exact enumeration. + TooLarge(usize), +} + +impl core::fmt::Display for QioError { + fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result { + match self { + QioError::BadCoupling { i, j } => write!(f, "coupling ({i}, {j}) is out of range or on the diagonal"), + QioError::Parse { line, message } => write!(f, "line {line}: {message}"), + QioError::TooLarge(n) => write!(f, "{n} spins is too many to enumerate (at most 34)"), + } + } +} + +impl std::error::Error for QioError {} + +impl Ising { + /// Build from couplings `(i, j, J_ij)` (each pair once; repeats add) and + /// fields `h` (empty for none). + pub fn new(n: usize, couplings: &[(usize, usize, f64)], h: &[f64]) -> Result { + let mut rows: Vec> = vec![Vec::new(); n]; + for &(i, j, v) in couplings { + if i >= n || j >= n || i == j { + return Err(QioError::BadCoupling { i, j }); + } + rows[i].push((j as u32, v)); + rows[j].push((i as u32, v)); + } + let mut start = vec![0]; + let (mut col, mut val) = (Vec::new(), Vec::new()); + for r in rows.iter_mut() { + r.sort_by_key(|e| e.0); + let mut k = 0; + while k < r.len() { + let (c, mut v) = r[k]; + k += 1; + while k < r.len() && r[k].0 == c { + v += r[k].1; + k += 1; + } + col.push(c); + val.push(v); + } + start.push(col.len()); + } + let mut hv = vec![0.0; n]; + hv[..h.len().min(n)].copy_from_slice(&h[..h.len().min(n)]); + Ok(Ising { n, start, col, val, h: hv }) + } + + /// The Sherrington–Kirkpatrick spin glass: all-to-all couplings drawn + /// uniformly from ±1, no fields, reproducible from `seed`. + pub fn sk(n: usize, seed: u64) -> Ising { + let mut rng = Rng(seed); + let mut c = Vec::with_capacity(n * n.saturating_sub(1) / 2); + for i in 0..n { + for j in i + 1..n { + c.push((i, j, if rng.next() >> 63 == 0 { 1.0 } else { -1.0 })); + } + } + Ising::new(n, &c, &[]).expect("in-range couplings") + } + + pub fn n(&self) -> usize { + self.n + } + + /// Couplings `(j, J_ij)` of spin `i`. + pub fn row(&self, i: usize) -> impl Iterator + '_ { + (self.start[i]..self.start[i + 1]).map(|k| (self.col[k] as usize, self.val[k])) + } + + fn has_fields(&self) -> bool { + self.h.iter().any(|&x| x != 0.0) + } + + /// `E(s) = −½ Σ J_ij s_i s_j − Σ h_i s_i`, summed in a fixed order. + pub fn energy(&self, s: &[i8]) -> f64 { + let mut e = 0.0; + for i in 0..self.n { + let mut field = 0.0; + for (j, v) in self.row(i) { + field += v * f64::from(s[j]); + } + e -= f64::from(s[i]) * (0.5 * field + self.h[i]); + } + e + } + + /// Root-mean-square coupling over ordered pairs `i ≠ j`. + pub fn rms_coupling(&self) -> f64 { + if self.n < 2 { + return 0.0; + } + let sum: f64 = self.val.iter().map(|v| v * v).sum(); + (sum / (self.n as f64 * (self.n - 1) as f64)).sqrt() + } + + /// The same problem with fields folded into couplings to one extra spin + /// fixed at `+1`: `J_{i,n} = h_i`. + fn without_fields(&self) -> Ising { + let mut c = Vec::new(); + for i in 0..self.n { + for (j, v) in self.row(i) { + if i < j { + c.push((i, j, v)); + } + } + if self.h[i] != 0.0 { + c.push((i, self.n, self.h[i])); + } + } + Ising::new(self.n + 1, &c, &[]).expect("in-range couplings") + } +} + +/// A Max-Cut problem: maximise `C = Σ_{edges} w_ij (1 − s_i s_j) / 2`. +#[derive(Clone, Debug, PartialEq)] +pub struct MaxCut { + pub n: usize, + pub edges: Vec<(usize, usize, f64)>, +} + +impl MaxCut { + /// The Ising problem with the same optima: `J = −w`, so that + /// `C = −E/2 + ¼ Σ_ij w_ij`. + pub fn ising(&self) -> Ising { + let c: Vec<(usize, usize, f64)> = self.edges.iter().map(|&(i, j, w)| (i, j, -w)).collect(); + Ising::new(self.n, &c, &[]).expect("a graph's edges are in range") + } + /// The cut of `s`, summed over edges in order. + pub fn cut(&self, s: &[i8]) -> f64 { + self.edges.iter().filter(|&&(i, j, _)| s[i] != s[j]).map(|e| e.2).sum() + } +} + +/// Read a graph in the G-set format: a header `n m`, then `m` lines +/// `i j w` with 1-based vertices. +pub fn read_gset(text: &str) -> Result { + let mut lines = text.lines().enumerate().map(|(k, l)| (k + 1, l.trim())).filter(|(_, l)| !l.is_empty()); + let err = |line: usize, m: &str| QioError::Parse { line, message: m.into() }; + let (hl, header) = lines.next().ok_or_else(|| err(1, "empty"))?; + let mut h = header.split_whitespace().map(|x| x.parse::()); + let (n, m) = match (h.next(), h.next()) { + (Some(Ok(n)), Some(Ok(m))) => (n, m), + _ => return Err(err(hl, "expected 'n m'")), + }; + let mut edges = Vec::with_capacity(m); + for (line, l) in lines { + let f: Vec<&str> = l.split_whitespace().collect(); + if f.len() < 3 { + return Err(err(line, "expected 'i j w'")); + } + let i: usize = f[0].parse().map_err(|_| err(line, "bad vertex"))?; + let j: usize = f[1].parse().map_err(|_| err(line, "bad vertex"))?; + let w: f64 = f[2].parse().map_err(|_| err(line, "bad weight"))?; + if i == 0 || j == 0 || i > n || j > n || i == j { + return Err(QioError::BadCoupling { i, j }); + } + edges.push((i - 1, j - 1, w)); + } + if edges.len() != m { + return Err(err(hl, &format!("header says {m} edges, found {}", edges.len()))); + } + Ok(MaxCut { n, edges }) +} + +// --------------------------------------------------------------------------- +// Randomness +// --------------------------------------------------------------------------- + +#[derive(Clone)] +struct Rng(u64); + +impl Rng { + fn next(&mut self) -> u64 { + self.0 = self.0.wrapping_add(0x9e37_79b9_7f4a_7c15); + let mut z = self.0; + z = (z ^ (z >> 30)).wrapping_mul(0xbf58_476d_1ce4_e5b9); + z = (z ^ (z >> 27)).wrapping_mul(0x94d0_49bb_1331_11eb); + z ^ (z >> 31) + } + /// Uniform on `[0, 1)`. + fn unit(&mut self) -> f64 { + (self.next() >> 11) as f64 * (1.0 / 9_007_199_254_740_992.0) + } + /// Uniform on `(−1, 1)`. + fn symmetric(&mut self) -> f64 { + 2.0 * self.unit() - 1.0 + } + #[cfg(test)] + fn below(&mut self, n: u64) -> u64 { + self.next() % n + } +} + +/// The stream of trial `t` under `seed`. +fn trial_rng(seed: u64, t: usize) -> Rng { + let mut r = Rng(seed ^ (t as u64).wrapping_mul(0xd1b5_4a32_d192_ed03)); + r.next(); + r +} + +#[cfg(test)] +fn spins(x: &[f64]) -> Vec { + x.iter().map(|&v| if v >= 0.0 { 1 } else { -1 }).collect() +} + +/// Run `trials` independent closures across threads; results come back in +/// trial order, so they do not depend on the thread count. wasm32 runs one. +fn run_trials(trials: usize, threads: usize, f: impl Fn(usize) -> T + Sync) -> Vec { + let threads = if cfg!(target_arch = "wasm32") { 1 } else { threads.max(1).min(trials.max(1)) }; + if threads == 1 { + return (0..trials).map(f).collect(); + } + let next = std::sync::atomic::AtomicUsize::new(0); + let slots: Vec>> = (0..trials).map(|_| std::sync::Mutex::new(None)).collect(); + std::thread::scope(|scope| { + for _ in 0..threads { + scope.spawn(|| loop { + let t = next.fetch_add(1, std::sync::atomic::Ordering::Relaxed); + if t >= trials { + break; + } + let r = f(t); + *slots[t].lock().unwrap() = Some(r); + }); + } + }); + slots.into_iter().map(|m| m.into_inner().unwrap().expect("every trial ran")).collect() +} + +// --------------------------------------------------------------------------- +// Simulated bifurcation +// --------------------------------------------------------------------------- + +/// Which simulated bifurcation. +#[derive(Clone, Copy, Debug, PartialEq)] +pub enum Variant { + /// `f = Jx`: fast to a good local solution. + Ballistic, + /// `f = J sgn(x)`: slower, more accurate. + Discrete, + /// Ballistic plus heating at rate `γ`. + HeatedBallistic(f64), + /// Discrete plus heating at rate `γ`. + HeatedDiscrete(f64), +} + +impl Variant { + fn discrete(&self) -> bool { + matches!(self, Variant::Discrete | Variant::HeatedDiscrete(_)) + } + fn gamma(&self) -> f64 { + match *self { + Variant::HeatedBallistic(g) | Variant::HeatedDiscrete(g) => g, + _ => 0.0, + } + } +} + +/// Settings for [`simulated_bifurcation`]. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct SbConfig { + pub variant: Variant, + /// Time steps per trial, `N_s`. + pub steps: u32, + pub dt: f64, + /// `c₀ = c₁ / (σ_J √N)`. + pub c1: f64, + pub a0: f64, + /// Evaluate the spins every this many steps (and at the end); a trial + /// reports its best. + pub sample_every: u32, + pub trials: usize, + pub seed: u64, +} + +impl SbConfig { + /// The published settings for dense spin glasses (arXiv:2203.08361): + /// - ballistic: `Δt = 0.7`, `c₁ = 0.6`; + /// - discrete: `Δt = 1.1`, `c₁ = 0.6`; + /// - heated ballistic: `Δt = 1.1`, `c₁ = 0.9`, `γ = 0.5`; + /// - heated discrete: `Δt = 1.1`, `c₁ = 0.7`, `γ = 0.06`. + pub fn published(variant: Variant, steps: u32, trials: usize, seed: u64) -> SbConfig { + let (variant, dt, c1) = match variant { + Variant::Ballistic => (variant, 0.7, 0.6), + Variant::Discrete => (variant, 1.1, 0.6), + Variant::HeatedBallistic(_) => (Variant::HeatedBallistic(0.5), 1.1, 0.9), + Variant::HeatedDiscrete(_) => (Variant::HeatedDiscrete(0.06), 1.1, 0.7), + }; + SbConfig { variant, steps, dt, c1, a0: 1.0, sample_every: 100, trials, seed } + } +} + +/// One trial's best configuration and its energy. +#[derive(Clone, Debug, PartialEq)] +pub struct Trial { + pub energy: f64, + pub spins: Vec, +} + +/// Every trial of a run, in trial order. +#[derive(Clone, Debug, PartialEq)] +pub struct Run { + pub trials: Vec, +} + +impl Run { + /// The lowest-energy trial (the first, on a tie). + pub fn best(&self) -> &Trial { + self.trials.iter().fold(&self.trials[0], |b, t| if t.energy < b.energy { t } else { b }) + } + /// The fraction of trials at or below `energy` (within `1e-9`). + pub fn success(&self, energy: f64) -> f64 { + self.trials.iter().filter(|t| t.energy <= energy + 1e-9).count() as f64 / self.trials.len() as f64 + } +} + +/// Trials advanced together through one pass over the couplings. Each trial's +/// arithmetic is unchanged, so batching never changes a bit. +const LANES: usize = 8; + +/// Run simulated bifurcation on `problem`: `cfg.trials` independent trials, +/// on up to `threads` threads. +pub fn simulated_bifurcation(problem: &Ising, cfg: &SbConfig, threads: usize) -> Run { + sb_with_state(problem, cfg, threads).0 +} + +/// [`simulated_bifurcation`], also returning each trial's final positions. +fn sb_with_state(problem: &Ising, cfg: &SbConfig, threads: usize) -> (Run, Vec>) { + let fields = problem.has_fields(); + let p = if fields { problem.without_fields() } else { problem.clone() }; + let n = p.n; + let sigma = p.rms_coupling(); + let c0 = if sigma > 0.0 { cfg.c1 / (sigma * (n as f64).sqrt()) } else { 0.0 }; + let (dt, a0, gamma, discrete) = (cfg.dt, cfg.a0, cfg.variant.gamma(), cfg.variant.discrete()); + const L: usize = LANES; + let blocks = run_trials(cfg.trials.div_ceil(L), threads, |b| { + // Lane `l` is trial `b·L + l`, initialised exactly as it would be + // alone: all of x from its stream, then all of y. + let mut x = vec![0.0; n * L]; + let mut y = vec![0.0; n * L]; + for l in 0..L { + let mut rng = trial_rng(cfg.seed, b * L + l); + for i in 0..n { + x[i * L + l] = rng.symmetric(); + } + for i in 0..n { + y[i * L + l] = rng.symmetric(); + } + } + let mut f = vec![0.0; n * L]; + let mut src = vec![0.0; n * L]; + let mut best: Vec = (0..L).map(|_| Trial { energy: f64::INFINITY, spins: Vec::new() }).collect(); + let consider = |x: &[f64], best: &mut Vec| { + for (l, b) in best.iter_mut().enumerate() { + let s: Vec = (0..n).map(|i| if x[i * L + l] >= 0.0 { 1 } else { -1 }).collect(); + let e = p.energy(&s); + if e < b.energy { + *b = Trial { energy: e, spins: s }; + } + } + }; + for k in 0..cfg.steps { + let a = a0 * f64::from(k) / f64::from(cfg.steps); + for (s, &v) in src.iter_mut().zip(&x) { + *s = if discrete { if v >= 0.0 { 1.0 } else { -1.0 } } else { v }; + } + for i in 0..n { + let mut acc = [0.0f64; L]; + for kk in p.start[i]..p.start[i + 1] { + let (v, c) = (p.val[kk], p.col[kk] as usize * L); + let row = &src[c..c + L]; + for l in 0..L { + acc[l] += v * row[l]; + } + } + f[i * L..i * L + L].copy_from_slice(&acc); + } + for ((xi, yi), &fi) in x.iter_mut().zip(y.iter_mut()).zip(&f) { + let y0 = *yi; + let mut yn = y0 + (-(a0 - a) * *xi + c0 * fi) * dt; + let mut xn = *xi + a0 * yn * dt; + if xn.abs() > 1.0 { + xn = if xn > 0.0 { 1.0 } else { -1.0 }; + yn = 0.0; + } + *xi = xn; + *yi = yn + gamma * y0 * dt; + } + if cfg.sample_every > 0 && (k + 1) % cfg.sample_every == 0 { + consider(&x, &mut best); + } + } + consider(&x, &mut best); + let finals: Vec> = (0..L).map(|l| (0..n).map(|i| x[i * L + l]).collect()).collect(); + best.into_iter().zip(finals).collect::>() + }); + let (trials, finals): (Vec, Vec>) = blocks.into_iter().flatten().take(cfg.trials).unzip(); + let trials = trials + .into_iter() + .map(|mut t| { + if fields { + // Gauge: the extra spin is +1, then drop it. + let g = t.spins[n - 1]; + t.spins.truncate(n - 1); + if g < 0 { + t.spins.iter_mut().for_each(|s| *s = -*s); + } + t.energy = problem.energy(&t.spins); + } + t + }) + .collect(); + (Run { trials }, finals) +} + +/// Greedy single-flip descent from `spins`: flip the spin that lowers the +/// energy most, until no flip does (ties go to the lowest index). The result +/// is a 1-flip local minimum. +pub fn descend(problem: &Ising, spins: &[i8]) -> Trial { + let n = problem.n; + let mut s = spins.to_vec(); + let mut field: Vec = (0..n).map(|i| problem.row(i).map(|(j, v)| v * f64::from(s[j])).sum::() + problem.h[i]).collect(); + loop { + // Flipping s_i changes E by 2 s_i field_i. + let (mut best_i, mut best_de) = (usize::MAX, -1e-12); + for i in 0..n { + let de = 2.0 * f64::from(s[i]) * field[i]; + if de < best_de { + best_de = de; + best_i = i; + } + } + if best_i == usize::MAX { + break; + } + let si = f64::from(s[best_i]); + s[best_i] = -s[best_i]; + for (j, v) in problem.row(best_i) { + field[j] -= 2.0 * v * si; + } + } + Trial { energy: problem.energy(&s), spins: s } +} + +// --------------------------------------------------------------------------- +// Simulated annealing +// --------------------------------------------------------------------------- + +/// Settings for [`simulated_annealing`]: `sweeps` sweeps over all spins, the +/// inverse temperature rising linearly from `beta0` to `beta1`. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct SaConfig { + pub sweeps: u32, + pub beta0: f64, + pub beta1: f64, + pub trials: usize, + pub seed: u64, +} + +/// Single-spin Metropolis annealing from a random start; each trial reports +/// the best configuration it visited. +pub fn simulated_annealing(problem: &Ising, cfg: &SaConfig, threads: usize) -> Run { + let n = problem.n; + let trials = run_trials(cfg.trials, threads, |t| { + let mut rng = trial_rng(cfg.seed, t); + let mut s: Vec = (0..n).map(|_| if rng.next() >> 63 == 0 { 1 } else { -1 }).collect(); + // Local fields `Σ_j J_ij s_j + h_i`. + let mut field: Vec = (0..n).map(|i| problem.row(i).map(|(j, v)| v * f64::from(s[j])).sum::() + problem.h[i]).collect(); + let mut e = problem.energy(&s); + let mut best = Trial { energy: e, spins: s.clone() }; + for sweep in 0..cfg.sweeps { + let beta = if cfg.sweeps > 1 { cfg.beta0 + (cfg.beta1 - cfg.beta0) * f64::from(sweep) / f64::from(cfg.sweeps - 1) } else { cfg.beta1 }; + for i in 0..n { + // Flipping s_i changes E by 2 s_i field_i. + let de = 2.0 * f64::from(s[i]) * field[i]; + let accept = de <= 0.0 || rng.unit() < exp(-beta * de); + if accept { + let si = f64::from(s[i]); + s[i] = -s[i]; + e += de; + for (j, v) in problem.row(i) { + field[j] -= 2.0 * v * si; + } + } + } + if e < best.energy - 1e-9 { + // Re-sum, so a reported energy never carries drift. + let exact = problem.energy(&s); + e = exact; + if exact < best.energy { + best = Trial { energy: exact, spins: s.clone() }; + } + } + } + best + }); + Run { trials } +} + +// --------------------------------------------------------------------------- +// Exact ground states +// --------------------------------------------------------------------------- + +/// The exact ground state, by enumerating all `2ⁿ⁻¹` configurations with +/// `s₀ = +1` in Gray-code order (both signs when there are fields), threaded +/// over the top bits. At most 34 spins. +pub fn ground_state(problem: &Ising, threads: usize) -> Result { + let n = problem.n; + if n > 34 { + return Err(QioError::TooLarge(n)); + } + if n == 0 { + return Ok(Trial { energy: 0.0, spins: Vec::new() }); + } + // With no fields the energy is even in s, so s₀ = +1 loses nothing. + let free = if problem.has_fields() { n } else { n - 1 }; + let offset = n - free; + let top = free.min(8).min(free.saturating_sub(4)); + let low = free - top; + let chunks = 1usize << top; + let results = run_trials(chunks, threads, |c| { + // Spins offset..offset+low run through a Gray code; the rest are + // fixed by the chunk. + let mut s = vec![1i8; n]; + for b in 0..top { + if (c >> b) & 1 == 1 { + s[offset + low + b] = -1; + } + } + let mut field: Vec = (0..n).map(|i| problem.row(i).map(|(j, v)| v * f64::from(s[j])).sum::() + problem.h[i]).collect(); + let mut e = problem.energy(&s); + let mut best = Trial { energy: e, spins: s.clone() }; + for k in 1u64..(1u64 << low) { + let i = offset + k.trailing_zeros() as usize; + let si = f64::from(s[i]); + e += 2.0 * si * field[i]; + s[i] = -s[i]; + for (j, v) in problem.row(i) { + field[j] -= 2.0 * v * si; + } + if e < best.energy - 1e-9 { + best = Trial { energy: problem.energy(&s), spins: s.clone() }; + e = best.energy; + } + } + best + }); + Ok(results.into_iter().fold(Trial { energy: f64::INFINITY, spins: Vec::new() }, |b, t| if t.energy < b.energy { t } else { b })) +} + +/// Steps (or sweeps) needed to reach a target with 99% probability, given a +/// per-trial success probability `p` at `steps` per trial: +/// `S = steps · ln 0.01 / ln(1 − p)`. `None` when `p = 0`. +pub fn step_to_solution(p: f64, steps: f64) -> Option { + if p <= 0.0 { + return None; + } + if p >= 0.99 { + return Some(steps); + } + Some(steps * crate::repro::ln(0.01) / crate::repro::ln(1.0 - p)) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn random_sparse(n: usize, edges: usize, seed: u64) -> MaxCut { + let mut rng = Rng(seed); + let mut e = Vec::new(); + while e.len() < edges { + let (i, j) = (rng.below(n as u64) as usize, rng.below(n as u64) as usize); + if i != j && !e.iter().any(|&(a, b, _)| (a, b) == (i.min(j), i.max(j))) { + e.push((i.min(j), i.max(j), if rng.next() >> 63 == 0 { 1.0 } else { -1.0 })); + } + } + MaxCut { n, edges: e } + } + + #[test] + fn energy_and_cut_agree() { + let g = random_sparse(12, 30, 3); + let p = g.ising(); + let total: f64 = g.edges.iter().map(|e| e.2).sum(); + let mut rng = Rng(9); + for _ in 0..20 { + let s: Vec = (0..12).map(|_| if rng.next() >> 63 == 0 { 1 } else { -1 }).collect(); + // C = −E/2 + ¼ Σ_ij w_ij, with the sum over ordered pairs. + assert!((g.cut(&s) - (-p.energy(&s) / 2.0 + total / 2.0)).abs() < 1e-12); + } + } + + #[test] + fn the_referee_enumerates_exhaustively() { + // Brute force over every configuration, against the Gray code. + let p = Ising::sk(10, 4); + let mut best = f64::INFINITY; + for m in 0..(1u32 << 10) { + let s: Vec = (0..10).map(|b| if (m >> b) & 1 == 1 { -1 } else { 1 }).collect(); + best = best.min(p.energy(&s)); + } + assert_eq!(ground_state(&p, 3).unwrap().energy, best); + // With fields, both signs of s₀ matter. + let q = Ising::new(6, &[(0, 1, 1.0), (1, 2, -2.0), (3, 4, 0.5)], &[0.3, -0.2, 0.0, 1.0, 0.0, -0.7]).unwrap(); + let mut best = f64::INFINITY; + for m in 0..64u32 { + let s: Vec = (0..6).map(|b| if (m >> b) & 1 == 1 { -1 } else { 1 }).collect(); + best = best.min(q.energy(&s)); + } + assert!((ground_state(&q, 1).unwrap().energy - best).abs() < 1e-12); + } + + #[test] + fn every_heuristic_finds_small_ground_states() { + for seed in 0..4 { + let p = Ising::sk(20, 100 + seed); + let exact = ground_state(&p, 4).unwrap().energy; + for v in [Variant::Ballistic, Variant::Discrete, Variant::HeatedBallistic(0.0), Variant::HeatedDiscrete(0.0)] { + let run = simulated_bifurcation(&p, &SbConfig::published(v, 1000, 32, seed), 4); + assert_eq!(run.best().energy, exact, "{v:?} seed {seed}"); + assert!((p.energy(&run.best().spins) - run.best().energy).abs() < 1e-12); + } + let sa = simulated_annealing(&p, &SaConfig { sweeps: 500, beta0: 0.1, beta1: 3.0, trials: 16, seed }, 4); + assert_eq!(sa.best().energy, exact, "SA seed {seed}"); + } + } + + #[test] + fn fields_ride_on_an_extra_spin() { + let mut rng = Rng(5); + let mut c = Vec::new(); + for i in 0..14 { + for j in i + 1..14 { + if rng.below(3) == 0 { + c.push((i, j, rng.symmetric())); + } + } + } + let h: Vec = (0..14).map(|_| rng.symmetric()).collect(); + let p = Ising::new(14, &c, &h).unwrap(); + let exact = ground_state(&p, 2).unwrap(); + let run = simulated_bifurcation(&p, &SbConfig::published(Variant::Discrete, 2000, 32, 1), 4); + assert!((run.best().energy - exact.energy).abs() < 1e-9, "{} vs {}", run.best().energy, exact.energy); + assert_eq!(run.best().spins.len(), 14); + } + + /// The one-trial-at-a-time form the batched kernel must match bit for bit. + fn scalar_reference(problem: &Ising, cfg: &SbConfig, threads: usize) -> (Run, Vec>) { + let fields = problem.has_fields(); + let p = if fields { problem.without_fields() } else { problem.clone() }; + let n = p.n; + let sigma = p.rms_coupling(); + let c0 = if sigma > 0.0 { cfg.c1 / (sigma * (n as f64).sqrt()) } else { 0.0 }; + let (dt, a0, gamma, discrete) = (cfg.dt, cfg.a0, cfg.variant.gamma(), cfg.variant.discrete()); + let trials = run_trials(cfg.trials, threads, |t| { + let mut rng = trial_rng(cfg.seed, t); + let mut x: Vec = (0..n).map(|_| rng.symmetric()).collect(); + let mut y: Vec = (0..n).map(|_| rng.symmetric()).collect(); + let mut f = vec![0.0; n]; + let mut src = vec![0.0; n]; + let mut best = Trial { energy: f64::INFINITY, spins: Vec::new() }; + let mut consider = |x: &[f64]| { + let s = spins(x); + let e = p.energy(&s); + if e < best.energy { + best = Trial { energy: e, spins: s }; + } + }; + for k in 0..cfg.steps { + let a = a0 * f64::from(k) / f64::from(cfg.steps); + for i in 0..n { + src[i] = if discrete { if x[i] >= 0.0 { 1.0 } else { -1.0 } } else { x[i] }; + } + for (i, fi) in f.iter_mut().enumerate() { + let mut acc = 0.0; + for kk in p.start[i]..p.start[i + 1] { + acc += p.val[kk] * src[p.col[kk] as usize]; + } + *fi = acc; + } + for i in 0..n { + let y0 = y[i]; + let mut yi = y0 + (-(a0 - a) * x[i] + c0 * f[i]) * dt; + let mut xi = x[i] + a0 * yi * dt; + if xi.abs() > 1.0 { + xi = if xi > 0.0 { 1.0 } else { -1.0 }; + yi = 0.0; + } + x[i] = xi; + y[i] = yi + gamma * y0 * dt; + } + if cfg.sample_every > 0 && (k + 1) % cfg.sample_every == 0 { + consider(&x); + } + } + consider(&x); + (best, x) + }); + let (trials, finals): (Vec, Vec>) = trials.into_iter().unzip(); + let trials = trials + .into_iter() + .map(|mut t| { + if fields { + // Gauge: the extra spin is +1, then drop it. + let g = t.spins[n - 1]; + t.spins.truncate(n - 1); + if g < 0 { + t.spins.iter_mut().for_each(|s| *s = -*s); + } + t.energy = problem.energy(&t.spins); + } + t + }) + .collect(); + (Run { trials }, finals) + } + + #[test] + fn batching_changes_no_bit() { + let g = random_sparse(60, 300, 21); + let q = Ising::new(15, &[(0, 3, 0.5), (3, 7, -1.25), (7, 14, 2.0), (2, 9, 1.0)], &[0.1, 0.0, -0.4, 0.0, 0.2, 0.0, 0.0, 0.3, 0.0, 0.0, 0.0, -0.5, 0.0, 0.0, 0.0]).unwrap(); + for p in [g.ising(), Ising::sk(33, 2), q] { + for v in [Variant::Ballistic, Variant::Discrete, Variant::HeatedBallistic(0.0), Variant::HeatedDiscrete(0.0)] { + // 13 trials: one full batch and a ragged one. + let cfg = SbConfig::published(v, 250, 13, 5); + let (run, finals) = sb_with_state(&p, &cfg, 3); + let (r_run, r_finals) = scalar_reference(&p, &cfg, 2); + assert_eq!(run, r_run, "{v:?}"); + let bits = |f: &Vec>| f.iter().flatten().map(|v| v.to_bits()).collect::>(); + assert_eq!(bits(&finals), bits(&r_finals), "{v:?}"); + } + } + } + + #[test] + fn runs_do_not_depend_on_threads() { + let p = Ising::sk(40, 8); + for v in [Variant::Discrete, Variant::HeatedBallistic(0.5)] { + let cfg = SbConfig::published(v, 300, 12, 77); + assert_eq!(simulated_bifurcation(&p, &cfg, 1), simulated_bifurcation(&p, &cfg, 5)); + } + let sa = SaConfig { sweeps: 100, beta0: 0.1, beta1: 2.0, trials: 9, seed: 3 }; + assert_eq!(simulated_annealing(&p, &sa, 1), simulated_annealing(&p, &sa, 4)); + } + + #[test] + fn gset_text_reads_and_refuses() { + let g = read_gset("4 3\n1 2 1\n2 3 1\n3 4 -1\n").unwrap(); + assert_eq!(g.n, 4); + assert_eq!(g.cut(&[1, -1, 1, -1]), 1.0); + assert!(read_gset("4 2\n1 2 1\n").is_err()); + assert!(read_gset("4 1\n1 1 1\n").is_err()); + assert!(read_gset("4 1\n1 5 1\n").is_err()); + } + + #[test] + fn descent_ends_at_a_local_minimum_never_higher() { + let p = Ising::sk(30, 6); + let mut rng = Rng(2); + for _ in 0..10 { + let s: Vec = (0..30).map(|_| if rng.next() >> 63 == 0 { 1 } else { -1 }).collect(); + let d = descend(&p, &s); + assert!(d.energy <= p.energy(&s)); + for i in 0..30 { + let mut t = d.spins.clone(); + t[i] = -t[i]; + assert!(p.energy(&t) >= d.energy - 1e-9); + } + } + } + + #[test] + fn step_to_solution_matches_its_formula() { + assert_eq!(step_to_solution(0.0, 1000.0), None); + let s = step_to_solution(0.5, 1000.0).unwrap(); + assert!((s - 1000.0 * (0.01f64).ln() / (0.5f64).ln()).abs() < 1e-9); + } +} diff --git a/wai/wai-quantum/src/repro.rs b/wai/wai-quantum/src/repro.rs index f740492b5..45620ed07 100644 --- a/wai/wai-quantum/src/repro.rs +++ b/wai/wai-quantum/src/repro.rs @@ -34,28 +34,28 @@ pub use trig::sin_cos; // Natural logarithm // --------------------------------------------------------------------------- -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio"))] const LN2_HI: f64 = 6.931_471_803_691_238_164_90e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio"))] const LN2_LO: f64 = 1.908_214_929_270_587_700_02e-10; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio"))] const LG1: f64 = 6.666_666_666_666_735_130e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio"))] const LG2: f64 = 3.999_999_999_940_941_908e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio"))] const LG3: f64 = 2.857_142_874_366_239_149e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio"))] const LG4: f64 = 2.222_219_843_214_978_396e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio"))] const LG5: f64 = 1.818_357_216_161_805_012e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio"))] const LG6: f64 = 1.531_383_769_920_937_332e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio"))] const LG7: f64 = 1.479_819_860_511_658_591e-1; /// The natural logarithm, within an ulp of the true value: `−∞` at zero, NaN /// below it. -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio"))] pub fn ln(x: f64) -> f64 { let mut bits = x.to_bits(); let mut hx = (bits >> 32) as u32; @@ -99,7 +99,7 @@ pub fn ln(x: f64) -> f64 { // Exponential // --------------------------------------------------------------------------- -#[cfg(any(feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch"))] +#[cfg(any(feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio"))] const EXP_P: [f64; 5] = [ 1.666_666_666_666_660_190_37e-1, -2.777_777_777_701_559_338_42e-3, @@ -111,7 +111,7 @@ const EXP_P: [f64; 5] = [ /// `eˣ`, bit-identical on every IEEE-754 platform: fdlibm's algorithm /// (reduction by `ln 2` in two parts, a rational approximation on /// `[−ln2/2, ln2/2]`), using only `+ − × ÷`. Within an ulp of the true value. -#[cfg(any(feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch"))] +#[cfg(any(feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio"))] pub fn exp(x: f64) -> f64 { const LN2_HI_E: f64 = 6.931_471_803_691_238_164_90e-1; const LN2_LO_E: f64 = 1.908_214_929_270_587_700_02e-10; @@ -312,7 +312,7 @@ mod tests { #[allow(unused_imports)] use super::*; - #[cfg(any(feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch"))] + #[cfg(any(feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio"))] #[test] fn exp_is_within_an_ulp_of_the_platform() { let mut s = 0x9e37_79b9_7f4a_7c15u64; @@ -355,7 +355,7 @@ mod tests { }) } - #[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch"))] + #[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio"))] #[test] fn ln_is_within_an_ulp_of_the_platform() { for x in inputs().filter(|&x| x > 0.0) { diff --git a/wai/wai-quantum/tests/bit_identity.rs b/wai/wai-quantum/tests/bit_identity.rs index 3bac232a3..9a66a7910 100644 --- a/wai/wai-quantum/tests/bit_identity.rs +++ b/wai/wai-quantum/tests/bit_identity.rs @@ -21,7 +21,8 @@ feature = "quantum_gbs", feature = "quantum_sv", feature = "quantum_arch", - feature = "quantum_ldpc" + feature = "quantum_ldpc", + feature = "quantum_qio" ))] use wai_quantum::*; @@ -364,3 +365,31 @@ fn qldpc_memory_decoding() { 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); +} diff --git a/wai/wai-rs/Cargo.lock b/wai/wai-rs/Cargo.lock index 4a2966ae1..aadbabb04 100644 --- a/wai/wai-rs/Cargo.lock +++ b/wai/wai-rs/Cargo.lock @@ -3900,7 +3900,7 @@ dependencies = [ [[package]] name = "wai-quantum" -version = "0.3.30" +version = "0.3.31" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-rs/Cargo.toml b/wai/wai-rs/Cargo.toml index 361264dda..ac0aec25b 100644 --- a/wai/wai-rs/Cargo.toml +++ b/wai/wai-rs/Cargo.toml @@ -367,6 +367,8 @@ quantum_sv = ["dep:wai-quantum", "wai-quantum/quantum_sv"] quantum_arch = ["dep:wai-quantum", "wai-quantum/quantum_arch"] # qLDPC memory at circuit level: bivariate-bicycle codes, BP-OSD. quantum_ldpc = ["quantum_frame", "wai-quantum/quantum_ldpc"] +# Quantum-inspired optimisation: simulated bifurcation, annealing, exact referee. +quantum_qio = ["dep:wai-quantum", "wai-quantum/quantum_qio"] quantum_qir = ["quantum_toolchain", "wai-quantum/quantum_qir"] # Variational quantum machine learning (extensions/quantum-sim § VML, # `wai.quantum.vml`): the interference-based ML the thermodynamic-computing pivot @@ -501,6 +503,7 @@ required-features = [ "quantum_sv", "quantum_arch", "quantum_ldpc", + "quantum_qio", ] [[bin]] diff --git a/wai/wai-rs/src/bin/wai_quantum.rs b/wai/wai-rs/src/bin/wai_quantum.rs index 95e02cbf4..d80cf5949 100644 --- a/wai/wai-rs/src/bin/wai_quantum.rs +++ b/wai/wai-rs/src/bin/wai_quantum.rs @@ -165,6 +165,10 @@ fn usage() -> ! { wai-quantum gbs --squeezing r1,r2,… [--eta E] [--seed S | --unitary U.json] [--json]\n \ wai-quantum gbs … --sample K (exact click patterns)\n \ wai-quantum gbs … --check SAMPLES.txt (observed patterns vs exact marginals)\n\n\ + OPTIMISATION (simulated bifurcation, annealing, exact):\n \ + wai-quantum qio (--gset FILE | --sk N [--instance S]) [--method dsb|bsb|hdsb|hbsb|sa|exact]\n \ + [--steps N] [--trials T] [--dt D] [--c1 C] [--gamma G] [--beta0 B --beta1 B] [--descend]\n \ + [--target V] [--seed S] [--threads T] [--spins] [--json]\n\n\ QLDPC MEMORY (circuit level, BP-OSD):\n \ wai-quantum ldpc [--code 72|90|108|144|288] [--p P] [--cycles N] [--shots S] [--basis z|x|both]\n \ [--max-iter I] [--osd-order K] [--threads T] [--seed S] [--json]\n \ @@ -965,6 +969,81 @@ fn main() { } } + // Quantum-inspired optimisation: Max-Cut (G-set text) or a spin + // glass, by simulated bifurcation, annealing, or exact enumeration. + // Max-Cut reports cuts; a spin glass reports energies. + "qio" => { + use wai::quantum_qio::*; + let (problem, graph) = if let Some(path) = opt(args, "--gset", "") { + let text = std::fs::read_to_string(&path).unwrap_or_else(|e| die(&format!("{path}: {e}"))); + let g = read_gset(&text).unwrap_or_else(|e| die(&e.to_string())); + (g.ising(), Some(g)) + } else if let Some(n) = opt(args, "--sk", "") { + let n: usize = n.parse().unwrap_or_else(|_| die("bad --sk size")); + (Ising::sk(n, num(args, "--instance", "", 0)), None) + } else { + die("qio needs --gset FILE or --sk N") + }; + let score = |s: &[i8]| graph.as_ref().map_or_else(|| problem.energy(s), |g| g.cut(s)); + let maximise = graph.is_some(); + let method = opt(args, "--method", "").unwrap_or_else(|| "dsb".into()); + let threads: usize = num(args, "--threads", "", std::thread::available_parallelism().map_or(1, |n| n.get())); + let seed: u64 = num(args, "--seed", "", 1); + let trials: usize = num(args, "--trials", "", 64); + let run = match method.as_str() { + "exact" => Run { trials: vec![ground_state(&problem, threads).unwrap_or_else(|e| die(&e.to_string()))] }, + "sa" => { + let cfg = SaConfig { sweeps: num(args, "--steps", "", 10_000), beta0: num(args, "--beta0", "", 0.1), beta1: num(args, "--beta1", "", 3.0), trials, seed }; + simulated_annealing(&problem, &cfg, threads) + } + m => { + let v = match m { + "bsb" => Variant::Ballistic, + "dsb" => Variant::Discrete, + "hbsb" => Variant::HeatedBallistic(0.0), + "hdsb" => Variant::HeatedDiscrete(0.0), + _ => die("--method: dsb, bsb, hdsb, hbsb, sa or exact"), + }; + let mut cfg = SbConfig::published(v, num(args, "--steps", "", 10_000), trials, seed); + cfg.dt = num(args, "--dt", "", cfg.dt); + cfg.c1 = num(args, "--c1", "", cfg.c1); + if let Some(g) = opt(args, "--gamma", "") { + let g: f64 = g.parse().unwrap_or_else(|_| die("bad --gamma")); + cfg.variant = match cfg.variant { + Variant::HeatedBallistic(_) => Variant::HeatedBallistic(g), + Variant::HeatedDiscrete(_) => Variant::HeatedDiscrete(g), + other => other, + }; + } + simulated_bifurcation(&problem, &cfg, threads) + } + }; + let finals: Vec> = run.trials.iter().map(|t| if flag(args, "--descend") { descend(&problem, &t.spins).spins } else { t.spins.clone() }).collect(); + let scores: Vec = finals.iter().map(|s| score(s)).collect(); + let best_k = (0..scores.len()).fold(0, |b, k| if (maximise && scores[k] > scores[b]) || (!maximise && scores[k] < scores[b]) { k } else { b }); + let best = scores[best_k]; + let mean = scores.iter().sum::() / scores.len() as f64; + let target: Option = opt(args, "--target", "").map(|t| t.parse().unwrap_or_else(|_| die("bad --target"))); + let hit = target.map(|t| scores.iter().filter(|&&v| if maximise { v >= t - 1e-9 } else { v <= t + 1e-9 }).count() as f64 / scores.len() as f64); + let bits: String = finals[best_k].iter().map(|&s| if s > 0 { '0' } else { '1' }).collect(); + let what = if maximise { "cut" } else { "energy" }; + if flag(args, "--json") { + let hit_s = hit.map_or("null".into(), |h| h.to_string()); + let spins_s = if flag(args, "--spins") { format!(",\"spins\":\"{bits}\"") } else { String::new() }; + println!("{{\"best\":{best},\"mean\":{mean},\"method\":\"{method}\",\"n\":{},\"objective\":\"{what}\",\"success\":{hit_s},\"trials\":{}{spins_s}}}", problem.n(), scores.len()); + } else { + println!("{} spins, {method}, {} trial(s)", problem.n(), scores.len()); + println!("best {what} {best}"); + println!("mean {what} {mean:.3}"); + if let (Some(t), Some(h)) = (target, hit) { + println!("reached {t} {:.2}% of trials", 100.0 * h); + } + if flag(args, "--spins") { + println!("spins {bits} (1 = −1; re-check: one pass over the edges)"); + } + } + } + // qLDPC memory: a bivariate-bicycle code's depth-8 syndrome cycle, // sampled with circuit-level noise and decoded by BP-OSD. "ldpc" => { diff --git a/wai/wai-rs/src/lib.rs b/wai/wai-rs/src/lib.rs index 00d5c98e9..e372519fe 100644 --- a/wai/wai-rs/src/lib.rs +++ b/wai/wai-rs/src/lib.rs @@ -651,6 +651,9 @@ pub use wai_quantum::quantum_arch; #[cfg(feature = "quantum_ldpc")] pub use wai_quantum::quantum_ldpc; +#[cfg(feature = "quantum_qio")] +pub use wai_quantum::quantum_qio; + #[cfg(any(feature = "quantum_tn", feature = "quantum_sv"))] pub use wai_quantum::gates; #[cfg(feature = "quantum_qir")] diff --git a/wai/wai-web/Cargo.lock b/wai/wai-web/Cargo.lock index 4116b39ae..321aa28d8 100644 --- a/wai/wai-web/Cargo.lock +++ b/wai/wai-web/Cargo.lock @@ -734,7 +734,7 @@ dependencies = [ [[package]] name = "wai-quantum" -version = "0.3.30" +version = "0.3.31" dependencies = [ "blake3", "ed25519-dalek", -- 2.54.0 (Apple Git-157)