open-standards/wai/wai-quantum/src/quantum_ldpc.rs

//! qLDPC MEMORY at circuit level — `wai.quantum.ldpc`.
//!
//! High-rate quantum LDPC codes promise many logical qubits per physical
//! qubit. The bivariate-bicycle family (arXiv:2308.07915) is the
//! best-documented case: 12 logical qubits in 144 data and 144 check qubits at
//! distance 12. The claim that matters is how often such a memory fails per
//! syndrome cycle when every gate, idle, preparation and measurement is noisy.
//! That number comes from a circuit simulation and a decoder.
//!
//! This module builds both, so the claim can be re-derived rather than quoted:
//!
//! - **The codes.** [`BbCode`] builds the check matrices from the two
//!   polynomials and counts the logical qubits by rank. It also finds a basis
//!   of logical operators. Five codes from the source are presets.
//! - **The circuit.** [`memory_circuit`] writes the source's depth-8 syndrome
//!   cycle in the ecosystem's circuit text:
//!   - each check qubit is coupled to its six neighbours in the published
//!     order;
//!   - X checks are measured in the X basis and Z checks in the Z basis, all
//!     interleaved;
//!   - circuit-level depolarizing noise of strength `p` follows every CNOT,
//!     idle, preparation and measurement;
//!   - two noiseless cycles at the end read the final syndrome.
//!
//!   `quantum_frame` samples it and derives its detector error model.
//! - **The decoder.** [`BpOsd`] is belief propagation with ordered-statistics
//!   post-processing:
//!   - min-sum belief propagation in the log domain, with the adaptive
//!     scaling `α = 1 − 2^(−t)`;
//!   - when it fails to converge, ordered statistics: columns are ranked by
//!     the posterior, an information set is taken in that order, and the
//!     combination sweep is scored by `Σ ln(1/p)`.
//!
//!   It runs in `f64` with in-crate logarithms, so a decode is the same bits on
//!   every machine.
//! - **The experiment.** [`memory_experiment`] runs Z-basis or X-basis memory,
//!   sampled in fixed chunks across threads, and counts logical failures.
//!   [`LdpcPoint::per_cycle`] turns them into a rate per syndrome cycle.
//!
//! # Checked against the source's own software
//!
//! The source's simulation software is Apache-2.0. Running it gives the
//! oracle:
//! - **The error model.** Its decoding matrices for `[[72,12,6]]` over six
//!   cycles at `p = 0.003` have the same columns as the ones derived here
//!   from the circuit, detector for detector: 2,268 for X faults and 2,232
//!   for Z faults. Half the probabilities agree to `10⁻¹³`. The largest gap
//!   is 0.94%, on columns where several faults merge: the source adds their
//!   probabilities, while this crate combines them as independent events.
//! - **Logical failure rates.** These agree with the source's decoder on the
//!   same codes and noise. On `[[72,12,6]]` over six cycles, the failure
//!   rates per shot are:
//!
//!   | `p`   | here (8,000 shots) | source's software (2,000 trials) |
//!   |-------|--------------------|----------------------------------|
//!   | 0.003 | 0.028              | 0.034                            |
//!   | 0.004 | 0.090              | 0.086                            |
//!   | 0.005 | 0.230              | 0.209                            |
//!
//!   "Here" combines the two halves as independent. The source decodes both
//!   halves of each run, and a Y fault can fail both together: in its software
//!   at `p = 0.005` both fail in 2.9% of runs, against 1.3% if they were
//!   independent. That accounts for most of the gap in the last row.
//!
//!   Taken half by half at `p = 0.005`, the source's software (9,000 runs for
//!   the X-type half, 3,000 for the Z-type) gives 0.111 and 0.118. This
//!   crate gives 0.121 (20,000 shots) and 0.122 (4,000). The Z-type half
//!   agrees. The X-type half sits 2.3σ above.
//!
//!   Each piece of the chain was checked separately:
//!   - the error models match column for column;
//!   - the source's decoder, given this crate's syndromes and error model,
//!     makes the same decision as [`BpOsd`] on every one of 1,500 shots;
//!   - decoding with the source's priors instead of these moves 476
//!     failures out of 4,000 to 478;
//!   - both samplers' detector rates match the exact expectation of the
//!     error model.
//!
//!   The residual lies in sampling, at a level not yet resolved.
//! - **The 144-qubit code.** On `[[144,12,12]]` over twelve cycles, the
//!   failure rate per run is 0.034 at `p = 0.004` and 0.180 at `p = 0.005`
//!   (8,000 and 2,000 shots). The source's software gives 0.039 (1,522 runs)
//!   and 0.193 (519 runs). Per cycle, the rates land on the source's
//!   published curve: 2.90·10⁻³ at `p = 0.004` (fit 3.15·10⁻³) and
//!   1.64·10⁻² at `p = 0.005` (fit 1.54·10⁻²).
//!
//! # Honest boundaries
//!
//! - **Memory only.** Logical gates on these codes (automorphisms, the
//!   ZX-duality, measurements through an ancilla system) are not simulated.
//! - **Two experiments, not one.** A Z-basis memory sees X-type failures and an
//!   X-basis memory sees Z-type failures. Their rates are combined as
//!   independent. The source decodes both halves of one noisy run instead,
//!   which differs at second order.
//! - **The decoder is the source's.** It is not the best decoder known for
//!   these codes. A better one would make the code look better. The point here
//!   is to reproduce the claim on its own terms.

use crate::quantum_frame::{error_model, parse, sample, Basis, ErrorModel, FrameError};
use crate::repro::ln;

// ---------------------------------------------------------------------------
// GF(2)
// ---------------------------------------------------------------------------

/// A set of GF(2) rows over `cols` columns, kept in reduced echelon form as
/// they arrive, so membership and rank are cheap.
#[derive(Clone, Debug)]
struct Span {
    /// Reduced rows with the column of their leading one.
    rows: Vec<(usize, Vec<u64>)>,
}

impl Span {
    fn new(cols: usize) -> Span {
        let _ = cols;
        Span { rows: Vec::new() }
    }
    fn reduce(&self, v: &mut [u64]) {
        for (lead, r) in &self.rows {
            if (v[lead / 64] >> (lead % 64)) & 1 == 1 {
                for (a, b) in v.iter_mut().zip(r) {
                    *a ^= b;
                }
            }
        }
    }
    /// Add `v` unless it is already in the span; whether it was new.
    fn insert(&mut self, v: &[u64]) -> bool {
        let mut v = v.to_vec();
        self.reduce(&mut v);
        let Some(w) = v.iter().position(|&x| x != 0) else { return false };
        let lead = w * 64 + v[w].trailing_zeros() as usize;
        // Keep every row reduced at the new leading column.
        for (_, r) in self.rows.iter_mut() {
            if (r[lead / 64] >> (lead % 64)) & 1 == 1 {
                for (a, b) in r.iter_mut().zip(&v) {
                    *a ^= b;
                }
            }
        }
        self.rows.push((lead, v));
        true
    }
    fn rank(&self) -> usize {
        self.rows.len()
    }
}

fn bits(cols: usize, ones: impl IntoIterator<Item = usize>) -> Vec<u64> {
    let mut v = vec![0u64; cols.div_ceil(64)];
    for c in ones {
        v[c / 64] ^= 1 << (c % 64);
    }
    v
}

fn ones(v: &[u64]) -> Vec<usize> {
    let mut out = Vec::new();
    for (w, &x) in v.iter().enumerate() {
        let mut x = x;
        while x != 0 {
            out.push(w * 64 + x.trailing_zeros() as usize);
            x &= x - 1;
        }
    }
    out
}

/// A basis of the null space of `rows` (each a list of columns) over `cols`
/// columns.
fn kernel(rows: &[Vec<usize>], cols: usize) -> Vec<Vec<u64>> {
    // Row-reduce, then read the kernel off the free columns.
    let mut span = Span::new(cols);
    for r in rows {
        span.insert(&bits(cols, r.iter().copied()));
    }
    let leads: Vec<usize> = span.rows.iter().map(|(l, _)| *l).collect();
    let mut is_lead = vec![false; cols];
    for &l in &leads {
        is_lead[l] = true;
    }
    let mut out = Vec::new();
    for free in (0..cols).filter(|&c| !is_lead[c]) {
        let mut v = bits(cols, [free]);
        for (lead, r) in &span.rows {
            if (r[free / 64] >> (free % 64)) & 1 == 1 {
                v[lead / 64] ^= 1 << (lead % 64);
            }
        }
        out.push(v);
    }
    out
}

// ---------------------------------------------------------------------------
// Bivariate-bicycle codes
// ---------------------------------------------------------------------------

/// A bivariate-bicycle code: `A = x^a₁ + y^a₂ + y^a₃` and
/// `B = y^b₁ + x^b₂ + x^b₃` over `x^ℓ = y^m = 1`, with `H_X = [A | B]` and
/// `H_Z = [Bᵀ | Aᵀ]`.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct BbCode {
    pub l: u32,
    pub m: u32,
    pub a: [u32; 3],
    pub b: [u32; 3],
}

/// Which monomial of `A` or `B` a check's neighbour comes from.
#[derive(Clone, Copy)]
enum Mono {
    X(u32),
    Y(u32),
}

impl BbCode {
    /// `[[72,12,6]]`.
    pub fn bb72() -> BbCode {
        BbCode { l: 6, m: 6, a: [3, 1, 2], b: [3, 1, 2] }
    }
    /// `[[90,8,10]]`.
    pub fn bb90() -> BbCode {
        BbCode { l: 15, m: 3, a: [9, 1, 2], b: [0, 2, 7] }
    }
    /// `[[108,8,10]]`.
    pub fn bb108() -> BbCode {
        BbCode { l: 9, m: 6, a: [3, 1, 2], b: [3, 1, 2] }
    }
    /// `[[144,12,12]]`.
    pub fn bb144() -> BbCode {
        BbCode { l: 12, m: 6, a: [3, 1, 2], b: [3, 1, 2] }
    }
    /// `[[288,12,18]]`.
    pub fn bb288() -> BbCode {
        BbCode { l: 12, m: 12, a: [3, 2, 7], b: [3, 1, 2] }
    }

    /// Checks of each type, `ℓm`.
    pub fn half(&self) -> usize {
        (self.l * self.m) as usize
    }
    /// Data qubits, `2ℓm`.
    pub fn n(&self) -> usize {
        2 * self.half()
    }

    fn apply(&self, mono: Mono, i: usize) -> usize {
        let (l, m) = (self.l as usize, self.m as usize);
        let (r, c) = (i / m, i % m);
        match mono {
            Mono::X(s) => ((r + s as usize) % l) * m + c,
            Mono::Y(s) => r * m + (c + s as usize) % m,
        }
    }
    fn apply_t(&self, mono: Mono, i: usize) -> usize {
        let (l, m) = (self.l as usize, self.m as usize);
        let (r, c) = (i / m, i % m);
        match mono {
            Mono::X(s) => ((r + l - s as usize % l) % l) * m + c,
            Mono::Y(s) => r * m + (c + m - s as usize % m) % m,
        }
    }
    fn a_mono(&self) -> [Mono; 3] {
        [Mono::X(self.a[0]), Mono::Y(self.a[1]), Mono::Y(self.a[2])]
    }
    fn b_mono(&self) -> [Mono; 3] {
        [Mono::Y(self.b[0]), Mono::X(self.b[1]), Mono::X(self.b[2])]
    }

    /// The six data qubits of X check `i`, by direction 0 to 5: `A₁(i), A₂(i),
    /// A₃(i)` on the left half, then `B₁(i), B₂(i), B₃(i)` on the right.
    pub fn x_neighbours(&self, i: usize) -> [usize; 6] {
        let (a, b, h) = (self.a_mono(), self.b_mono(), self.half());
        [self.apply(a[0], i), self.apply(a[1], i), self.apply(a[2], i), h + self.apply(b[0], i), h + self.apply(b[1], i), h + self.apply(b[2], i)]
    }
    /// The six data qubits of Z check `i`: `B₁ᵀ(i), B₂ᵀ(i), B₃ᵀ(i)` on the
    /// left, then `A₁ᵀ(i), A₂ᵀ(i), A₃ᵀ(i)` on the right.
    pub fn z_neighbours(&self, i: usize) -> [usize; 6] {
        let (a, b, h) = (self.a_mono(), self.b_mono(), self.half());
        [self.apply_t(b[0], i), self.apply_t(b[1], i), self.apply_t(b[2], i), h + self.apply_t(a[0], i), h + self.apply_t(a[1], i), h + self.apply_t(a[2], i)]
    }

    /// X checks as lists of data qubits.
    pub fn hx(&self) -> Vec<Vec<usize>> {
        (0..self.half()).map(|i| self.x_neighbours(i).to_vec()).collect()
    }
    /// Z checks as lists of data qubits.
    pub fn hz(&self) -> Vec<Vec<usize>> {
        (0..self.half()).map(|i| self.z_neighbours(i).to_vec()).collect()
    }

    /// Logical qubits, `n − rank H_X − rank H_Z`.
    pub fn k(&self) -> usize {
        let rank = |rows: Vec<Vec<usize>>| {
            let mut s = Span::new(self.n());
            for r in rows {
                s.insert(&bits(self.n(), r));
            }
            s.rank()
        };
        self.n() - rank(self.hx()) - rank(self.hz())
    }

    /// Whether every X check commutes with every Z check.
    pub fn is_valid(&self) -> bool {
        let hz = self.hz();
        self.hx().iter().all(|x| hz.iter().all(|z| x.iter().filter(|q| z.contains(q)).count() % 2 == 0))
    }

    /// A basis of `k` logical operators of one type. `Basis::Z` gives Z-type
    /// operators, which commute with the X checks and are not products of Z
    /// checks; `Basis::X` gives the reverse. Each is a list of data qubits.
    pub fn logicals(&self, basis: Basis) -> Vec<Vec<usize>> {
        let (commute, stabilizers) = match basis {
            Basis::Z => (self.hx(), self.hz()),
            Basis::X => (self.hz(), self.hx()),
        };
        let n = self.n();
        let mut span = Span::new(n);
        for s in stabilizers {
            span.insert(&bits(n, s));
        }
        let mut out = Vec::new();
        for v in kernel(&commute, n) {
            if span.insert(&v) {
                out.push(ones(&v));
            }
        }
        out
    }
}

// ---------------------------------------------------------------------------
// The syndrome-measurement circuit
// ---------------------------------------------------------------------------

/// The order in which X checks (`SX`) and Z checks (`SZ`) visit their six
/// neighbours over the seven CNOT rounds of a cycle (arXiv:2308.07915, its
/// table of the depth-8 cycle). `None` is a round without a CNOT.
const SX: [Option<usize>; 7] = [None, Some(1), Some(4), Some(3), Some(5), Some(0), Some(2)];
const SZ: [Option<usize>; 7] = [Some(3), Some(5), Some(0), Some(1), Some(2), Some(4), None];

fn line(out: &mut String, name: &str, arg: Option<f64>, targets: &[usize]) {
    if targets.is_empty() {
        return;
    }
    out.push_str(name);
    if let Some(a) = arg {
        out.push_str(&format!("({a})"));
    }
    for t in targets {
        out.push_str(&format!(" {t}"));
    }
    out.push('\n');
}

/// Text of a bivariate-bicycle memory experiment in `basis`.
///
/// - **Layout.** Qubits are X checks `0..ℓm`, data `ℓm..3ℓm` (left half, then
///   right) and Z checks `3ℓm..4ℓm`.
/// - **Start.** Data and Z checks start in `basis`, noiselessly.
/// - **Cycles.** `cycles` noisy syndrome cycles at strength `p`, then two
///   noiseless ones.
/// - **Detectors.** Each check of `basis`, cycle by cycle: its first outcome,
///   then each change.
/// - **Observables.** The logical operators of [`BbCode::logicals`], read from
///   a final noiseless measurement of the data.
///
/// Noise:
/// - two-qubit depolarizing noise after each CNOT;
/// - one-qubit depolarizing noise on idle data;
/// - a flip after each preparation;
/// - a flipped result on each measurement.
pub fn memory_circuit(code: &BbCode, cycles: u32, p: f64, basis: Basis) -> String {
    let h = code.half();
    let xc: Vec<usize> = (0..h).collect();
    let data: Vec<usize> = (h..3 * h).collect();
    let zc: Vec<usize> = (3 * h..4 * h).collect();
    let xn: Vec<[usize; 6]> = (0..h).map(|i| code.x_neighbours(i).map(|q| q + h)).collect();
    let zn: Vec<[usize; 6]> = (0..h).map(|i| code.z_neighbours(i).map(|q| q + h)).collect();
    let mut out = String::new();
    match basis {
        Basis::Z => line(&mut out, "R", None, &data),
        Basis::X => line(&mut out, "RX", None, &data),
    }
    line(&mut out, "R", None, &zc);
    let mut measured = 0usize;
    // Absolute record of each check's outcome in the previous cycle.
    let mut previous: Option<Vec<usize>> = None;
    for cycle in 0..cycles + 2 {
        let noise = if cycle < cycles { Some(p) } else { None };
        let idle_cnots = |out: &mut String, t: usize, x_rounds: bool, z_rounds: bool| {
            let mut pairs: Vec<usize> = Vec::new();
            let mut busy = vec![false; 4 * h];
            if let (true, Some(dir)) = (z_rounds, SZ[t]) {
                for i in 0..h {
                    pairs.extend([zn[i][dir], zc[i]]);
                    busy[zn[i][dir]] = true;
                }
            }
            if let (true, Some(dir)) = (x_rounds, SX[t]) {
                for i in 0..h {
                    pairs.extend([xc[i], xn[i][dir]]);
                    busy[xn[i][dir]] = true;
                }
            }
            line(out, "CX", None, &pairs);
            if let Some(p) = noise {
                line(out, "DEPOLARIZE2", Some(p), &pairs);
                let idle: Vec<usize> = data.iter().copied().filter(|&q| !busy[q]).collect();
                line(out, "DEPOLARIZE1", Some(p), &idle);
            }
        };
        // Round 1: prepare the X checks; the Z checks take their first CNOT.
        line(&mut out, "RX", None, &xc);
        if let Some(p) = noise {
            line(&mut out, "Z_ERROR", Some(p), &xc);
        }
        idle_cnots(&mut out, 0, false, true);
        out.push_str("TICK\n");
        // Rounds 2 to 6: both kinds of check.
        for t in 1..6 {
            idle_cnots(&mut out, t, true, true);
            out.push_str("TICK\n");
        }
        // Round 7: measure the Z checks; the X checks take their last CNOT.
        line(&mut out, "M", noise, &zc);
        let z_records: Vec<usize> = (measured..measured + h).collect();
        measured += h;
        idle_cnots(&mut out, 6, true, false);
        out.push_str("TICK\n");
        // Round 8: data idle; measure the X checks; prepare the Z checks.
        if let Some(p) = noise {
            line(&mut out, "DEPOLARIZE1", Some(p), &data);
        }
        line(&mut out, "MX", noise, &xc);
        let x_records: Vec<usize> = (measured..measured + h).collect();
        measured += h;
        line(&mut out, "R", None, &zc);
        if let Some(p) = noise {
            line(&mut out, "X_ERROR", Some(p), &zc);
        }
        out.push_str("TICK\n");
        let current = match basis {
            Basis::Z => z_records,
            Basis::X => x_records,
        };
        for (i, &r) in current.iter().enumerate() {
            match &previous {
                None => out.push_str(&format!("DETECTOR rec[-{}]\n", measured - r)),
                Some(prev) => out.push_str(&format!("DETECTOR rec[-{}] rec[-{}]\n", measured - r, measured - prev[i])),
            }
        }
        previous = Some(current);
    }
    match basis {
        Basis::Z => line(&mut out, "M", None, &data),
        Basis::X => line(&mut out, "MX", None, &data),
    }
    let first = measured;
    measured += data.len();
    for (k, op) in code.logicals(basis).iter().enumerate() {
        out.push_str(&format!("OBSERVABLE_INCLUDE({k})"));
        for &q in op {
            out.push_str(&format!(" rec[-{}]", measured - (first + q)));
        }
        out.push('\n');
    }
    out
}

// ---------------------------------------------------------------------------
// BP-OSD
// ---------------------------------------------------------------------------

/// How belief propagation scales its min-sum messages.
#[derive(Clone, Copy, Debug, PartialEq)]
pub enum Scaling {
    /// `α = 1 − 2^(−t)` at iteration `t`.
    Adaptive,
    Fixed(f64),
}

/// Decoder settings. The default is the source's: min-sum, adaptive scaling,
/// up to 10,000 iterations, combination sweep of order 7.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct BpOsdConfig {
    pub max_iter: u32,
    pub scaling: Scaling,
    /// Order of the combination sweep; 0 is OSD-0 alone.
    pub osd_order: u32,
}

impl Default for BpOsdConfig {
    fn default() -> BpOsdConfig {
        BpOsdConfig { max_iter: 10_000, scaling: Scaling::Adaptive, osd_order: 7 }
    }
}

/// A decode: the observables the correction flips, and how it was found.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct Decoded {
    pub observables: u64,
    /// Whether belief propagation alone explained the syndrome.
    pub converged: bool,
    pub iterations: u32,
}

/// Belief propagation with ordered-statistics post-processing over a
/// detector error model.
#[derive(Clone, Debug)]
pub struct BpOsd {
    cfg: BpOsdConfig,
    checks: usize,
    vars: usize,
    /// Edges by check: `(var)`, CSR.
    row_start: Vec<usize>,
    row_var: Vec<usize>,
    /// Edges by variable: indices into the check-ordered edge list, CSR.
    col_start: Vec<usize>,
    col_edge: Vec<usize>,
    prior: Vec<f64>,
    /// `ln(1/p)`, the OSD cost of each fault.
    cost: Vec<f64>,
    observables: Vec<u64>,
    rank: usize,
}

impl BpOsd {
    /// Build the decoder for `model`. Faults with probability 0 are dropped.
    pub fn new(model: &ErrorModel, cfg: BpOsdConfig) -> BpOsd {
        let mechs: Vec<_> = model.mechanisms.iter().filter(|m| m.probability > 0.0 && !m.detectors.is_empty()).collect();
        let checks = model.detectors as usize;
        let vars = mechs.len();
        let mut by_row: Vec<Vec<usize>> = vec![Vec::new(); checks];
        for (j, m) in mechs.iter().enumerate() {
            for &d in &m.detectors {
                by_row[d as usize].push(j);
            }
        }
        let mut row_start = vec![0];
        let mut row_var: Vec<usize> = Vec::new();
        for r in &by_row {
            row_var.extend(r);
            row_start.push(row_var.len());
        }
        let mut by_col: Vec<Vec<usize>> = vec![Vec::new(); vars];
        for (e, &j) in row_var.iter().enumerate() {
            by_col[j].push(e);
        }
        let mut col_start = vec![0];
        let mut col_edge: Vec<usize> = Vec::new();
        for c in &by_col {
            col_edge.extend(c);
            col_start.push(col_edge.len());
        }
        let prior: Vec<f64> = mechs.iter().map(|m| ln((1.0 - m.probability) / m.probability)).collect();
        let cost: Vec<f64> = mechs.iter().map(|m| ln(1.0 / m.probability)).collect();
        let observables = mechs.iter().map(|m| m.observables).collect();
        let mut span = Span::new(vars);
        // Rank of H: row-reduce the checks over the faults.
        for r in &by_row {
            span.insert(&bits(vars, r.iter().copied()));
        }
        let rank = span.rank();
        BpOsd { cfg, checks, vars, row_start, row_var, col_start, col_edge, prior, cost, observables, rank }
    }

    /// Faults (columns) in the model the decoder uses.
    pub fn faults(&self) -> usize {
        self.vars
    }

    /// Decode the detectors that fired.
    pub fn decode(&self, fired: &[u32]) -> Decoded {
        let mut synd = vec![false; self.checks];
        for &d in fired {
            synd[d as usize] = true;
        }
        let edges = self.row_var.len();
        let mut v2c = vec![0.0f64; edges];
        let mut c2v = vec![0.0f64; edges];
        let mut sgn = vec![0u32; edges];
        let mut post = vec![0.0f64; self.vars];
        let mut hard = vec![false; self.vars];
        for j in 0..self.vars {
            for &e in &self.col_edge[self.col_start[j]..self.col_start[j + 1]] {
                v2c[e] = self.prior[j];
            }
        }
        let mut iterations = 0;
        for t in 1..=self.cfg.max_iter {
            iterations = t;
            let alpha = match self.cfg.scaling {
                Scaling::Adaptive => 1.0 - pow2_neg(t),
                Scaling::Fixed(a) => a,
            };
            // Check to variable: the smallest other magnitude, the parity of
            // the other signs and the syndrome; zero counts as negative.
            for (i, &fired) in synd.iter().enumerate() {
                let (s, e_end) = (self.row_start[i], self.row_start[i + 1]);
                let mut min = 1e308f64;
                let mut negatives = u32::from(fired);
                for e in s..e_end {
                    c2v[e] = min;
                    sgn[e] = negatives;
                    if v2c[e].abs() < min {
                        min = v2c[e].abs();
                    }
                    if v2c[e] <= 0.0 {
                        negatives += 1;
                    }
                }
                let (mut min, mut negatives) = (1e308f64, 0u32);
                for e in (s..e_end).rev() {
                    if min < c2v[e] {
                        c2v[e] = min;
                    }
                    sgn[e] += negatives;
                    c2v[e] *= if sgn[e] % 2 == 1 { -alpha } else { alpha };
                    if v2c[e].abs() < min {
                        min = v2c[e].abs();
                    }
                    if v2c[e] <= 0.0 {
                        negatives += 1;
                    }
                }
            }
            // Variable to check: the prior plus every other incoming message,
            // summed forward then backward.
            for j in 0..self.vars {
                let col = &self.col_edge[self.col_start[j]..self.col_start[j + 1]];
                let mut acc = self.prior[j];
                for &e in col {
                    v2c[e] = acc;
                    acc += c2v[e];
                }
                post[j] = acc;
                hard[j] = acc <= 0.0;
                let mut acc = 0.0;
                for &e in col.iter().rev() {
                    v2c[e] += acc;
                    acc += c2v[e];
                }
            }
            if self.explains(&hard, &synd) {
                return Decoded { observables: self.flips(&hard), converged: true, iterations };
            }
        }
        let x = self.osd(&post, &synd);
        Decoded { observables: self.flips(&x), converged: false, iterations }
    }

    fn explains(&self, x: &[bool], synd: &[bool]) -> bool {
        (0..self.checks).all(|i| {
            let parity = self.row_var[self.row_start[i]..self.row_start[i + 1]].iter().filter(|&&j| x[j]).count() % 2 == 1;
            parity == synd[i]
        })
    }

    fn flips(&self, x: &[bool]) -> u64 {
        x.iter().zip(&self.observables).filter(|(b, _)| **b).fold(0, |m, (_, o)| m ^ o)
    }

    /// Ordered-statistics decoding: rank faults by posterior (likeliest
    /// first), take the first independent columns as the information set, and
    /// solve; then sweep single flips of every other column and pairs among
    /// the first `osd_order`, keeping the cheapest solution.
    fn osd(&self, post: &[f64], synd: &[bool]) -> Vec<bool> {
        let mut order: Vec<usize> = (0..self.vars).collect();
        order.sort_by(|&a, &b| post[a].total_cmp(&post[b]).then(a.cmp(&b)));
        // Row-major matrix over the sorted columns, the syndrome appended.
        let cols = self.vars + 1;
        let words = cols.div_ceil(64);
        let mut m = vec![0u64; self.checks * words];
        let mut position = vec![0usize; self.vars];
        for (p, &j) in order.iter().enumerate() {
            position[j] = p;
        }
        for i in 0..self.checks {
            for &j in &self.row_var[self.row_start[i]..self.row_start[i + 1]] {
                let p = position[j];
                m[i * words + p / 64] ^= 1 << (p % 64);
            }
            if synd[i] {
                let p = self.vars;
                m[i * words + p / 64] ^= 1 << (p % 64);
            }
        }
        // Reduced echelon form, pivots in sorted order.
        let mut pivots: Vec<usize> = Vec::with_capacity(self.rank);
        let mut row = 0;
        for p in 0..self.vars {
            if row == self.checks || pivots.len() == self.rank {
                break;
            }
            let (w, b) = (p / 64, p % 64);
            let Some(r) = (row..self.checks).find(|&r| (m[r * words + w] >> b) & 1 == 1) else { continue };
            if r != row {
                for k in 0..words {
                    m.swap(r * words + k, row * words + k);
                }
            }
            let (head, tail) = m.split_at_mut(row * words);
            let (pivot, rest) = tail.split_at_mut(words);
            for other in head.chunks_exact_mut(words).chain(rest.chunks_exact_mut(words)) {
                // Whole rows: the free columns before `p` are read later.
                if (other[w] >> b) & 1 == 1 {
                    for (a, c) in other.iter_mut().zip(pivot.iter()) {
                        *a ^= c;
                    }
                }
            }
            pivots.push(p);
            row += 1;
        }
        let bit = |r: usize, p: usize| (m[r * words + p / 64] >> (p % 64)) & 1 == 1;
        let s = self.vars;
        // OSD-0: the pivots take the reduced syndrome; everything else is 0.
        let mut base = vec![false; self.vars];
        for (r, &p) in pivots.iter().enumerate() {
            base[p] = bit(r, s);
        }
        let cost_at = |p: usize| self.cost[order[p]];
        let base_cost: f64 = (0..self.vars).filter(|&p| base[p]).map(cost_at).sum();
        let mut best = (base_cost, Vec::new());
        if self.cfg.osd_order > 0 {
            let mut is_pivot = vec![false; self.vars];
            for &p in &pivots {
                is_pivot[p] = true;
            }
            let free: Vec<usize> = (0..self.vars).filter(|&p| !is_pivot[p]).collect();
            // The change in cost from flipping a set of free columns.
            let delta = |set: &[usize]| {
                let mut d: f64 = set.iter().map(|&p| cost_at(p)).sum();
                for (r, &pv) in pivots.iter().enumerate() {
                    let flip = set.iter().fold(false, |a, &p| a ^ bit(r, p));
                    if flip {
                        d += if base[pv] { -cost_at(pv) } else { cost_at(pv) };
                    }
                }
                d
            };
            let order_k = (self.cfg.osd_order as usize).min(free.len());
            let mut candidates: Vec<Vec<usize>> = free.iter().map(|&p| vec![p]).collect();
            for i in 0..order_k {
                for j in i + 1..order_k {
                    candidates.push(vec![free[i], free[j]]);
                }
            }
            for set in candidates {
                let c = base_cost + delta(&set);
                if c < best.0 {
                    best = (c, set);
                }
            }
        }
        let mut x = base;
        for &p in &best.1 {
            x[p] = true;
            for (r, &pv) in pivots.iter().enumerate() {
                if bit(r, p) {
                    x[pv] = !x[pv];
                }
            }
        }
        let mut out = vec![false; self.vars];
        for (p, &j) in order.iter().enumerate() {
            out[j] = x[p];
        }
        out
    }
}

fn pow2_neg(t: u32) -> f64 {
    if t >= 1075 {
        return 0.0;
    }
    let mut x = 1.0;
    for _ in 0..t {
        x /= 2.0;
    }
    x
}

// ---------------------------------------------------------------------------
// Memory experiments
// ---------------------------------------------------------------------------

/// A memory experiment's tally.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct LdpcPoint {
    pub p: f64,
    pub cycles: u32,
    pub shots: u64,
    /// Shots whose decoded observables were wrong (any logical flipped).
    pub failures: u64,
    /// Shots belief propagation alone decoded.
    pub converged: u64,
}

impl LdpcPoint {
    /// Failure rate per shot.
    pub fn per_shot(&self) -> f64 {
        self.failures as f64 / self.shots as f64
    }
    /// Failure rate per syndrome cycle, `1 − (1 − P)^(1/N)`.
    pub fn per_cycle(&self) -> f64 {
        1.0 - crate::repro::exp(ln(1.0 - self.per_shot()) / f64::from(self.cycles))
    }
}

/// Shots per independently seeded chunk: results do not depend on how many
/// threads ran them.
const CHUNK: usize = 64;

/// Run `shots` shots of the memory experiment in `basis` and decode each.
/// Chunk `c` is sampled with seed `seed + c`. The tally is the same at any
/// `threads` (one on wasm32, which has none).
#[allow(clippy::too_many_arguments)]
pub fn memory_experiment(code: &BbCode, cycles: u32, p: f64, basis: Basis, shots: usize, seed: u64, cfg: BpOsdConfig, threads: usize) -> Result<LdpcPoint, FrameError> {
    let circuit = parse(&memory_circuit(code, cycles, p, basis))?;
    let model = error_model(&circuit)?;
    let decoder = BpOsd::new(&model, cfg);
    let chunks = shots.div_ceil(CHUNK);
    let run = |c: usize| {
        let n = CHUNK.min(shots - c * CHUNK);
        let det = sample(&circuit, n, seed.wrapping_add(c as u64));
        let (mut fail, mut conv) = (0u64, 0u64);
        for s in 0..n {
            let d = decoder.decode(&det.fired(s));
            fail += u64::from(d.observables != det.flips(s));
            conv += u64::from(d.converged);
        }
        (fail, conv)
    };
    // wasm32 has no threads; the tally is the same either way.
    let threads = if cfg!(target_arch = "wasm32") { 1 } else { threads.max(1).min(chunks.max(1)) };
    let mut tallies = vec![(0u64, 0u64); chunks];
    if threads == 1 {
        for (c, t) in tallies.iter_mut().enumerate() {
            *t = run(c);
        }
    } else {
        let next = std::sync::atomic::AtomicUsize::new(0);
        let results = std::sync::Mutex::new(&mut tallies);
        std::thread::scope(|scope| {
            for _ in 0..threads {
                scope.spawn(|| loop {
                    let c = next.fetch_add(1, std::sync::atomic::Ordering::Relaxed);
                    if c >= chunks {
                        break;
                    }
                    let r = run(c);
                    results.lock().unwrap()[c] = r;
                });
            }
        });
    }
    let (failures, converged) = tallies.iter().fold((0, 0), |(f, c), t| (f + t.0, c + t.1));
    Ok(LdpcPoint { p, cycles, shots: shots as u64, failures, converged })
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn the_presets_have_their_published_parameters() {
        for (code, n, k) in [(BbCode::bb72(), 72, 12), (BbCode::bb90(), 90, 8), (BbCode::bb108(), 108, 8), (BbCode::bb144(), 144, 12), (BbCode::bb288(), 288, 12)] {
            assert!(code.is_valid());
            assert_eq!((code.n(), code.k()), (n, k));
            for basis in [Basis::X, Basis::Z] {
                assert_eq!(code.logicals(basis).len(), k);
            }
        }
    }

    #[test]
    fn logical_operators_commute_with_checks_and_pair_up() {
        let code = BbCode::bb72();
        let (lx, lz) = (code.logicals(Basis::X), code.logicals(Basis::Z));
        let overlap = |a: &[usize], b: &[usize]| a.iter().filter(|q| b.contains(q)).count() % 2;
        for z in &lz {
            assert!(code.hx().iter().all(|c| overlap(c, z) == 0));
        }
        for x in &lx {
            assert!(code.hz().iter().all(|c| overlap(c, x) == 0));
        }
        // The pairing matrix between the two bases is invertible.
        let mut span = Span::new(lz.len());
        for x in &lx {
            span.insert(&bits(lz.len(), lz.iter().enumerate().filter(|(_, z)| overlap(x, z) == 1).map(|(i, _)| i)));
        }
        assert_eq!(span.rank(), lz.len());
    }

    #[test]
    fn the_noiseless_cycle_is_deterministic_and_silent() {
        for basis in [Basis::Z, Basis::X] {
            let c = parse(&memory_circuit(&BbCode::bb72(), 3, 0.0, basis)).unwrap();
            let model = error_model(&c).unwrap();
            assert!(model.mechanisms.iter().all(|m| m.probability == 0.0) || model.mechanisms.is_empty());
            let det = sample(&c, 64, 1);
            assert!((0..64).all(|s| det.fired(s).is_empty() && det.flips(s) == 0));
        }
    }

    #[test]
    fn the_error_model_has_the_reference_columns() {
        // The source's decoding matrices for [[72,12,6]], six cycles, p = 0.003.
        let reference = include_str!("../tests/data/bb72_p003_c6.dem");
        type Column = (f64, Vec<u32>);
        let mut sections: Vec<(char, Vec<Column>)> = Vec::new();
        for l in reference.lines() {
            if let Some(h) = l.strip_prefix("# ") {
                sections.push((h.chars().next().unwrap(), Vec::new()));
                continue;
            }
            let mut f = l.split_whitespace();
            let p: f64 = f.next().unwrap().parse().unwrap();
            let rows: Vec<u32> = f.map(|x| x.parse().unwrap()).collect();
            // Its matrix drops the logical rows, so columns that differ only
            // there repeat; their faults add.
            let cols = &mut sections.last_mut().unwrap().1;
            match cols.iter_mut().find(|(_, r)| *r == rows) {
                Some((q, _)) => *q += p,
                None => cols.push((p, rows)),
            }
        }
        for (kind, columns) in sections {
            // X faults are seen by Z checks: the Z-basis memory.
            let basis = if kind == 'X' { Basis::Z } else { Basis::X };
            let c = parse(&memory_circuit(&BbCode::bb72(), 6, 0.003, basis)).unwrap();
            let model = error_model(&c).unwrap();
            let mut ours: Vec<(Vec<u32>, f64)> = Vec::new();
            for m in model.mechanisms.iter().filter(|m| !m.detectors.is_empty()) {
                match ours.iter_mut().find(|(d, _)| *d == m.detectors) {
                    // The same detectors with different observables: the
                    // source keeps logical rows out of its decoding matrix.
                    Some((_, q)) => *q = *q + m.probability - 2.0 * *q * m.probability,
                    None => ours.push((m.detectors.clone(), m.probability)),
                }
            }
            // The source keeps one column of faults that light nothing (its
            // probabilities summed); here they leave the model. None of them
            // flips a logical.
            let columns: Vec<_> = columns.into_iter().filter(|(_, rows)| !rows.is_empty()).collect();
            assert!(model.mechanisms.iter().all(|m| !m.detectors.is_empty() || m.observables == 0));
            assert_eq!(ours.len(), columns.len(), "{kind}");
            for (p, rows) in &columns {
                let (_, q) = ours.iter().find(|(d, _)| d == rows).unwrap_or_else(|| panic!("{kind}: no column {rows:?}"));
                assert!((q / p - 1.0).abs() < 0.01, "{kind} {rows:?}: {q} vs {p}");
            }
        }
    }

    #[test]
    fn bp_osd_corrects_every_single_fault() {
        let c = parse(&memory_circuit(&BbCode::bb72(), 2, 0.003, Basis::Z)).unwrap();
        let model = error_model(&c).unwrap();
        let dec = BpOsd::new(&model, BpOsdConfig::default());
        for m in model.mechanisms.iter().filter(|m| !m.detectors.is_empty()) {
            assert_eq!(dec.decode(&m.detectors).observables, m.observables, "{:?}", m.detectors);
        }
    }

    #[test]
    fn the_tally_does_not_depend_on_threads() {
        let cfg = BpOsdConfig { max_iter: 50, ..BpOsdConfig::default() };
        let a = memory_experiment(&BbCode::bb72(), 2, 0.004, Basis::Z, 600, 9, cfg, 1).unwrap();
        let b = memory_experiment(&BbCode::bb72(), 2, 0.004, Basis::Z, 600, 9, cfg, 3).unwrap();
        assert_eq!(a, b);
        assert!(a.failures > 0 && a.converged > 0);
    }
}