wai-quantum 0.3.41

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