wai-quantum 0.3.31: quantum-inspired optimisation, a re-runnable classical baseline #151
14 changed files with 1139 additions and 22 deletions
2
wai/wai-quantum-serve/Cargo.lock
generated
2
wai/wai-quantum-serve/Cargo.lock
generated
|
|
@ -589,7 +589,7 @@ dependencies = [
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.3.30"
|
||||
version = "0.3.31"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"ed25519-dalek",
|
||||
|
|
|
|||
2
wai/wai-quantum/Cargo.lock
generated
2
wai/wai-quantum/Cargo.lock
generated
|
|
@ -415,7 +415,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.3.30"
|
||||
version = "0.3.31"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"ed25519-dalek",
|
||||
|
|
|
|||
|
|
@ -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"]
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
74
wai/wai-quantum/examples/qio_gset.rs
Normal file
74
wai/wai-quantum/examples/qio_gset.rs
Normal file
|
|
@ -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<String> = 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<f64> = 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<f64> = 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::<f64>() / 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<f64> = 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::<f64>().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<f64> = 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::<f64>() / 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);
|
||||
}
|
||||
|
|
@ -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;
|
||||
|
||||
|
|
|
|||
895
wai/wai-quantum/src/quantum_qio.rs
Normal file
895
wai/wai-quantum/src/quantum_qio.rs
Normal file
|
|
@ -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<usize>,
|
||||
col: Vec<u32>,
|
||||
val: Vec<f64>,
|
||||
h: Vec<f64>,
|
||||
}
|
||||
|
||||
/// 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<Ising, QioError> {
|
||||
let mut rows: Vec<Vec<(u32, f64)>> = 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<Item = (usize, f64)> + '_ {
|
||||
(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<MaxCut, QioError> {
|
||||
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::<usize>());
|
||||
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<i8> {
|
||||
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<T: Send>(trials: usize, threads: usize, f: impl Fn(usize) -> T + Sync) -> Vec<T> {
|
||||
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<std::sync::Mutex<Option<T>>> = (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<i8>,
|
||||
}
|
||||
|
||||
/// Every trial of a run, in trial order.
|
||||
#[derive(Clone, Debug, PartialEq)]
|
||||
pub struct Run {
|
||||
pub trials: Vec<Trial>,
|
||||
}
|
||||
|
||||
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<Vec<f64>>) {
|
||||
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<Trial> = (0..L).map(|_| Trial { energy: f64::INFINITY, spins: Vec::new() }).collect();
|
||||
let consider = |x: &[f64], best: &mut Vec<Trial>| {
|
||||
for (l, b) in best.iter_mut().enumerate() {
|
||||
let s: Vec<i8> = (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<Vec<f64>> = (0..L).map(|l| (0..n).map(|i| x[i * L + l]).collect()).collect();
|
||||
best.into_iter().zip(finals).collect::<Vec<_>>()
|
||||
});
|
||||
let (trials, finals): (Vec<Trial>, Vec<Vec<f64>>) = 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<f64> = (0..n).map(|i| problem.row(i).map(|(j, v)| v * f64::from(s[j])).sum::<f64>() + 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<i8> = (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<f64> = (0..n).map(|i| problem.row(i).map(|(j, v)| v * f64::from(s[j])).sum::<f64>() + 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<Trial, QioError> {
|
||||
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<f64> = (0..n).map(|i| problem.row(i).map(|(j, v)| v * f64::from(s[j])).sum::<f64>() + 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<f64> {
|
||||
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<i8> = (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<i8> = (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<i8> = (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<f64> = (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<Vec<f64>>) {
|
||||
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<f64> = (0..n).map(|_| rng.symmetric()).collect();
|
||||
let mut y: Vec<f64> = (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<Trial>, Vec<Vec<f64>>) = 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<Vec<f64>>| f.iter().flatten().map(|v| v.to_bits()).collect::<Vec<u64>>();
|
||||
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<i8> = (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);
|
||||
}
|
||||
}
|
||||
|
|
@ -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) {
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
}
|
||||
|
|
|
|||
2
wai/wai-rs/Cargo.lock
generated
2
wai/wai-rs/Cargo.lock
generated
|
|
@ -3900,7 +3900,7 @@ dependencies = [
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.3.30"
|
||||
version = "0.3.31"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"ed25519-dalek",
|
||||
|
|
|
|||
|
|
@ -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]]
|
||||
|
|
|
|||
|
|
@ -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<Vec<i8>> = run.trials.iter().map(|t| if flag(args, "--descend") { descend(&problem, &t.spins).spins } else { t.spins.clone() }).collect();
|
||||
let scores: Vec<f64> = 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::<f64>() / scores.len() as f64;
|
||||
let target: Option<f64> = 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" => {
|
||||
|
|
|
|||
|
|
@ -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")]
|
||||
|
|
|
|||
2
wai/wai-web/Cargo.lock
generated
2
wai/wai-web/Cargo.lock
generated
|
|
@ -734,7 +734,7 @@ dependencies = [
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.3.30"
|
||||
version = "0.3.31"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"ed25519-dalek",
|
||||
|
|
|
|||
Loading…
Reference in a new issue