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

//! Deterministic quantum LDPC DECODING — `wai.quantum.decode`
//! (extensions/quantum-ops § Decoding).
//!
//! A fault-tolerant quantum computer runs a classical decoder in its innermost
//! loop: every syndrome-measurement round produces a syndrome, and a decoder must
//! turn it into a correction fast enough and cheaply enough to keep up with the
//! code cycle. For quantum LDPC (qLDPC) codes — the high-rate codes that make
//! fault tolerance affordable — the decoder of record is **belief propagation**,
//! and plain BP *fails* on quantum codes: the stabilizer degeneracy and the short
//! cycles in the Tanner graph trap it in symmetric non-solutions. A published fix is
//! **Relay-BP**: an ensemble of BP passes with disordered per-node memory
//! strengths, relayed (warm-started) from one another, which breaks the symmetry.
//!
//! This module is Relay-BP done the WAI way: **normalized min-sum** message
//! passing in pure fixed-point integer arithmetic — no float, no `tanh`/`atanh` —
//! so a decode is **byte-identical on every machine** and the emitted correction
//! is a portable hash. Every decode seals a signed, joule-metered
//! [`crate::quantum_ops::DecodeReceipt`] carrying **joules-per-decode**.
//!
//! Min-sum on a single uniform channel is scale-invariant (check node = min,
//! variable node = sum, both homogeneous), so the channel prior is a fixed unit
//! and the exact physical error rate never enters the decoder — only the
//! deterministic error *sampler* uses it. That keeps the decode fully integer.
//!
//! # The codes verify their own claims
//!
//! A code is named `[[n, k, d]]`, and `d` is the half that matters: below it,
//! correction is guaranteed rather than likely. [`x_distance`] computes it exactly
//! by exhaustion, so `surface(3)` is *shown* to have distance 3 and `surface(5)`
//! distance 5 — every lighter logical operator ruled out — rather than being
//! assumed from the constructor's name. The decoder is held to the same standard:
//! every error of weight below `d` is decoded and checked, exhaustively.
//!
//! # Circuit-level noise: faults that move
//!
//! Phenomenological noise treats a measurement as a black box that sometimes
//! lies. A real extraction round is a *circuit*, and a fault inside it does not
//! stay put: a fault on a stabilizer's ancilla rides the remaining CNOTs and
//! lands as X errors on every data qubit that circuit has not yet touched — a
//! **hook error**. One fault, several correlated data errors, and no edge in the
//! matching graph that represents them.
//!
//! [`decode_circuit_level`] adds that channel. Surface code, `q = hook = p`,
//! `d` rounds, and the same run with hooks switched off:
//!
//! | `p`   | d=3    | d=5    | d=7    | d=5, hooks off |
//! |-------|--------|--------|--------|----------------|
//! | 0.002 | 0.0008 | 0.0000 | 0.0000 | 0.0000         |
//! | 0.005 | 0.0058 | 0.0033 | 0.0008 | 0.0008         |
//!
//! Hooks cost roughly a factor of four at `d = 5`, which is the point: they are
//! why circuit-level thresholds sit below phenomenological ones. With `hook = 0`
//! the model reduces to [`decode_phenomenological`] **bit for bit** — same
//! stream, same failures — which is asserted, not merely hoped.
//!
//! # Handing the model to other tools
//!
//! [`to_dem`] exports the noise model as a **detector error model** in Stim's DEM
//! syntax. A DEM says nothing about qubits or gates — it says here are the
//! independent ways this experiment can fault, here is which detectors each one
//! lights, and here is whether it flips the logical observable. That is exactly
//! the interface every matching decoder wants, which is why it is the QEC
//! ecosystem's interchange format, and it is to decoders what QIR is to compilers.
//!
//! The observable comes from [`logical_z_support`]: the Z checks detect X errors,
//! so failure is not "an error remained" but "the residual chain crossed the
//! code", and a minimum-weight logical Z is the ruler that measures it.
//!
//! The export is checked against the code's own algebra rather than eyeballed:
//! the mechanisms of a minimum-weight logical operator must light **no** detectors
//! yet flip the observable an **odd** number of times, a stabilizer's must light
//! none and flip it an **even** number, and no mechanism anywhere may be
//! undetectable while flipping the observable — that last one would be a logical
//! error no decoder could see coming, and would make every external consumer
//! silently mis-score.
//!
//! Ingest works too: [`from_dem`] parses a model — ours or somebody else's — and
//! [`decode_dem`] decodes a round of detector data against it, answering the only
//! question a QEC pipeline asks: did the logical observable flip? Decoding the
//! same shots through the exported DEM matches the native decoder:
//!
//! | d | `p`   | native | via DEM |
//! |---|-------|--------|---------|
//! | 3 | 0.01  | 0.0127 | 0.0113  |
//! | 5 | 0.005 | 0.0007 | 0.0007  |
//! | 5 | 0.02  | 0.0440 | 0.0440  |
//!
//! (The two judge success slightly differently — stabilizer membership against
//! observable parity — which coincide whenever the decoder fully explains the
//! syndrome, and account for the small gap at `d = 3`.)
//!
//! Ingest refuses rather than approximates: an unimplemented instruction is an
//! error, not a skipped line, because a DEM whose `repeat` block was ignored
//! describes a different experiment; and a hyperedge mechanism is refused by the
//! matching decoder instead of being dropped.
//!
//! # Two decoders, and which code each is for
//!
//! **Relay-BP** is a belief-propagation decoder and earns its name on the qLDPC
//! family it was designed for: on the bivariate-bicycle code the relay legs beat
//! plain BP (at `p = 0.06`, 0.30 → 0.27 logical failure), and more legs help.
//!
//! On **topological** codes it does not merely underperform — it inverts the
//! property that makes a code family worth having. Measured across `p = 0.001`
//! to `0.02`, a distance-5 surface code decodes *worse* than distance-3 under
//! Relay-BP, and the toric code shows the same inversion. That is the known
//! weakness of BP on degenerate codes with short cycles, not a defect in the code
//! constructions (both verify as valid CSS codes with the expected `k`).
//!
//! [`union_find_correct`] is the answer for those codes — a matching decoder that
//! pairs syndromes on the code's own graph. It restores the property BP loses
//! (surface code, logical failure rate, 2000 trials):
//!
//! | `p`   | d=3 BP | d=3 UF | d=5 BP | d=5 UF | d=7 BP | d=7 UF |
//! |-------|--------|--------|--------|--------|--------|--------|
//! | 0.005 | 0.0205 | 0.0010 | 0.0410 | 0.0000 | 0.0605 | 0.0000 |
//! | 0.01  | 0.0370 | 0.0020 | 0.0785 | 0.0000 | 0.1170 | 0.0000 |
//! | 0.02  | 0.0715 | 0.0090 | 0.1535 | 0.0020 | 0.2365 | 0.0005 |
//!
//! # Faulty measurement: decoding in space-time
//!
//! Everything above assumes the syndrome is read perfectly, which no device does.
//! [`decode_phenomenological`] runs `rounds` noisy measurement rounds plus a final
//! perfect readout and decodes the block as a **space-time** matching problem:
//! detectors are *differences* between consecutive rounds, so a data error lights
//! two checks at the same time while a lying measurement lights one check at two
//! adjacent times. Surface code, `q = p`, `d` rounds:
//!
//! | `p`   | d=3    | d=5    | d=7    |
//! |-------|--------|--------|--------|
//! | 0.003 | 0.0025 | 0.0000 | 0.0000 |
//! | 0.005 | 0.0050 | 0.0008 | 0.0008 |
//! | 0.01  | 0.0142 | 0.0042 | 0.0067 |
//! | 0.02  | 0.0475 | 0.0458 | 0.0417 |
//!
//! Distance pays below `p ≈ 1%` and the curves converge by `2%` — a threshold in
//! the region the literature puts the phenomenological surface code, which is the
//! outcome that would have exposed a broken implementation.
//!
//! Edges are weighted by `−ln(rate)`, so the decoder blames the likelier fault
//! rather than the nearer one. That matters more than it sounds: with data noise
//! off and measurement noise on, an unweighted decoder happily "explains" a lying
//! measurement by flipping data qubits, and enough of those in a line is a logical
//! error invented out of nothing.
//!
//! The near-linear cost is measured, not just claimed: single-shot decoding runs
//! about 1.2 µs per shot at `d = 3` and 7.5 µs at `d = 11`, so roughly six times
//! the work for thirteen times the qubits. Threshold studies are Monte Carlo, and
//! at those rates they are a coffee break rather than an overnight job.
//!
//! Distance now buys what it is supposed to buy, monotonically. Use union-find for
//! surface and toric codes, Relay-BP for the higher-degree qLDPC families where a
//! qubit touches more than two checks and the matching graph does not exist.
//!
//! Codes built in: the **rotated surface code** (`[[d², 1, d]]`, the planar code
//! hardware actually targets), the **toric code** (its periodic-boundary cousin)
//! and a
//! **bivariate-bicycle** qLDPC code (the high-rate family; the \[\[72,12\]\] preset is
//! the "gross"-style construction). Logical success — knowable only in simulation,
//! since it needs the true error — is reported by the benchmark harness, never put
//! in the receipt (a real decoder does not know it).

use crate::quantum_ops::{content_hash, DecodeReceipt, GrantRef};
use ed25519_dalek::SigningKey;

/// Fixed-point fractional bits for messages/damping. `1.0` = `1 << QEC_FRAC`.
pub const QEC_FRAC: u32 = 12;
const SCALE: i64 = 1 << QEC_FRAC;
/// Channel prior (uniform): min-sum is scale-invariant, so any positive unit works.
const L0: i64 = SCALE;
/// Message saturation, keeps fixed-point sums bounded without changing decisions.
const MSG_MAX: i64 = 64 * SCALE;
/// Normalized-min-sum scaling α = 7/8 (the standard hardware value).
const ALPHA_NUM: i64 = 7;
/// Per-node relay memory strength range. Relay-BP draws *disordered* memory,
/// including negative values (over-relaxation) that aggressively break the
/// symmetric trapping sets plain BP gets stuck in. γ ∈ [-0.6, 0.9].
const GAMMA_MIN: i64 = -(6 * SCALE) / 10;
const GAMMA_MAX: i64 = (9 * SCALE) / 10;

fn splitmix64(state: &mut u64) -> u64 {
    *state = state.wrapping_add(0x9E37_79B9_7F4A_7C15);
    let mut z = *state;
    z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
    z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
    z ^ (z >> 31)
}

// ===========================================================================
// GF(2) linear algebra — row space membership (logical-triviality test)
// ===========================================================================

/// A reduced GF(2) basis of a set of binary row vectors, supporting membership
/// testing. Rows are packed into `u64` words. Used to test whether a residual
/// error is a product of stabilizer generators (i.e. logically trivial).
#[derive(Clone, Debug)]
pub struct Gf2Basis {
    words: usize,
    /// Reduced rows, each with a distinct pivot column.
    rows: Vec<Vec<u64>>,
    /// Pivot column of each reduced row (parallel to `rows`).
    pivots: Vec<usize>,
}

fn bit_get(row: &[u64], i: usize) -> bool {
    (row[i >> 6] >> (i & 63)) & 1 == 1
}
fn bit_set(row: &mut [u64], i: usize) {
    row[i >> 6] |= 1u64 << (i & 63);
}
fn bit_xor(dst: &mut [u64], src: &[u64]) {
    for (d, s) in dst.iter_mut().zip(src) {
        *d ^= *s;
    }
}
fn is_zero(row: &[u64]) -> bool {
    row.iter().all(|&w| w == 0)
}

/// Fully reduced row-echelon form over GF(2), returned as `(pivot column, row)`.
///
/// [`Gf2Basis`] stops at row-echelon form — a pivot may still appear in earlier
/// rows — which is enough to test membership but not to read a row off as an
/// independent generator. Encoding needs the stronger form.
#[cfg(feature = "quantum")]
fn rref(rows_cols: &[Vec<usize>], n: usize) -> Vec<(usize, Vec<usize>)> {
    let mut m: Vec<Vec<bool>> = rows_cols
        .iter()
        .map(|cols| {
            let mut r = vec![false; n];
            for &c in cols {
                r[c] = true;
            }
            r
        })
        .collect();
    let mut pivots: Vec<usize> = Vec::new();
    let mut r = 0usize;
    for c in 0..n {
        let Some(sel) = (r..m.len()).find(|&i| m[i][c]) else { continue };
        m.swap(r, sel);
        for i in 0..m.len() {
            if i != r && m[i][c] {
                for j in 0..n {
                    let v = m[r][j];
                    m[i][j] ^= v;
                }
            }
        }
        pivots.push(c);
        r += 1;
        if r == m.len() {
            break;
        }
    }
    m.truncate(pivots.len());
    pivots
        .into_iter()
        .zip(m)
        .map(|(p, row)| (p, (0..n).filter(|&j| row[j]).collect()))
        .collect()
}

impl Gf2Basis {
    /// Whether a row lies in this basis's span — it reduces to zero against it.
    pub fn reduces_to_zero(&self, cols: &[usize], ncols: usize) -> bool {
        let mut v = vec![0u64; ncols.div_ceil(64).max(1)];
        for &c in cols {
            bit_set(&mut v, c);
        }
        for (bi, br) in self.rows.iter().enumerate() {
            if bit_get(&v, self.pivots[bi]) {
                bit_xor(&mut v, br);
            }
        }
        is_zero(&v)
    }

    /// Build from a list of rows, each given as the sorted set of its 1-columns.
    pub fn from_rows(rows_cols: &[Vec<usize>], ncols: usize) -> Gf2Basis {
        let words = ncols.div_ceil(64).max(1);
        let mut basis: Vec<Vec<u64>> = Vec::new();
        let mut pivots: Vec<usize> = Vec::new();
        for cols in rows_cols {
            let mut v = vec![0u64; words];
            for &c in cols {
                bit_set(&mut v, c);
            }
            // reduce against current basis
            for (bi, br) in basis.iter().enumerate() {
                if bit_get(&v, pivots[bi]) {
                    bit_xor(&mut v, br);
                }
            }
            if let Some(p) = Self::lowest_set(&v) {
                // eliminate this pivot from existing basis rows for a full RREF
                for br in basis.iter_mut() {
                    if bit_get(br, p) {
                        bit_xor(br, &v);
                    }
                }
                basis.push(v);
                pivots.push(p);
            }
        }
        Gf2Basis { words, rows: basis, pivots }
    }

    fn lowest_set(v: &[u64]) -> Option<usize> {
        for (w, &word) in v.iter().enumerate() {
            if word != 0 {
                return Some(w * 64 + word.trailing_zeros() as usize);
            }
        }
        None
    }

    /// Rank of the row space.
    pub fn rank(&self) -> usize {
        self.rows.len()
    }

    /// Is the vector (given as its 1-columns) in the row space? Reduce it against
    /// the basis and check it vanishes.
    pub fn contains(&self, cols: &[usize]) -> bool {
        let mut v = vec![0u64; self.words];
        for &c in cols {
            bit_set(&mut v, c);
        }
        for (bi, br) in self.rows.iter().enumerate() {
            if bit_get(&v, self.pivots[bi]) {
                bit_xor(&mut v, br);
            }
        }
        is_zero(&v)
    }
}

// ===========================================================================
// CSS code — check matrix (for decoding) + stabilizer basis (for success test)
// ===========================================================================

/// A CSS code prepared for decoding one error sector. `checks` are the parity
/// checks whose syndrome we decode (rows as sorted qubit lists); `stab_basis` is
/// the GF(2) row space of the complementary stabilizers, against which a residual
/// error is tested for logical triviality.
#[derive(Clone, Debug)]
pub struct CssCode {
    pub id: String,
    pub n: usize,
    /// Parity checks (syndrome producers), each a sorted list of qubit indices.
    pub checks: Vec<Vec<usize>>,
    /// The complementary stabilizer family, kept so callers can inspect the code
    /// (a CSS code is defined by both families, not just the one that produces the
    /// syndrome being decoded).
    pub stabs: Vec<Vec<usize>>,
    stab_basis: Gf2Basis,
    /// Number of encoded logical qubits, `k = n − rank(checks) − rank(stab)`.
    pub k: usize,
}

impl CssCode {
    fn new(id: String, n: usize, checks: Vec<Vec<usize>>, stab: Vec<Vec<usize>>) -> CssCode {
        let stab_basis = Gf2Basis::from_rows(&stab, n);
        let check_basis = Gf2Basis::from_rows(&checks, n);
        let k = n.saturating_sub(check_basis.rank() + stab_basis.rank());
        CssCode { id, n, checks, stabs: stab, stab_basis, k }
    }

    /// Canonical bytes of the check matrix (content-addresses the code).
    pub fn check_matrix_bytes(&self) -> Vec<u8> {
        let mut b = Vec::new();
        b.extend_from_slice(b"wai:qec-checks\x01");
        b.extend_from_slice(&(self.n as u64).to_le_bytes());
        for row in &self.checks {
            b.extend_from_slice(&(row.len() as u32).to_le_bytes());
            for &c in row {
                b.extend_from_slice(&(c as u32).to_le_bytes());
            }
        }
        b
    }

    /// A residual error (as its 1-qubits) is logically trivial iff it lies in the
    /// stabilizer row space.
    pub fn is_trivial(&self, residual_cols: &[usize]) -> bool {
        self.stab_basis.contains(residual_cols)
    }

    // --- Toric code -------------------------------------------------------

    /// The `L×L` toric code: `n = 2L²` qubits on edges, `L²` star (X) checks and
    /// `L²` plaquette (Z) checks, `k = 2`. We decode X-errors: `checks = Z`
    /// (plaquettes), triviality against the X-stabilizers (stars).
    pub fn toric(l: usize) -> CssCode {
        let n = 2 * l * l;
        let h = |i: usize, j: usize| (i % l) * l + (j % l); // horizontal edge id
        let v = |i: usize, j: usize| l * l + (i % l) * l + (j % l); // vertical edge id
        let mut plaq = Vec::with_capacity(l * l); // Z checks
        let mut star = Vec::with_capacity(l * l); // X stabilizers
        for i in 0..l {
            for j in 0..l {
                // plaquette (i,j): top/bottom horizontals, left/right verticals
                let mut p = vec![h(i, j), h(i + 1, j), v(i, j), v(i, j + 1)];
                p.sort_unstable();
                p.dedup();
                plaq.push(p);
                // star at vertex (i,j): incident horizontals + verticals
                let mut s = vec![h(i, j), h(i, j + l - 1), v(i, j), v(i + l - 1, j)];
                s.sort_unstable();
                s.dedup();
                star.push(s);
            }
        }
        CssCode::new(format!("toric:L{l}"), n, plaq, star)
    }

    // --- Rotated surface code ---------------------------------------------

    /// The distance-`d` **rotated surface code**, `[[d², 1, d]]` (`d` odd).
    ///
    /// The toric code above is this code's periodic-boundary cousin: elegant, but
    /// it needs a torus. The rotated surface code is the planar one — nearest-
    /// neighbour on a `d × d` patch with open boundaries — and it is what
    /// essentially every hardware roadmap actually targets.
    ///
    /// Data qubit `(r, c)` is index `r·d + c`. Bulk plaquettes cover each 2×2
    /// square and alternate in a checkerboard; the boundaries carry weight-2
    /// stabilizers, X along the top and bottom, Z along the left and right. The
    /// parities are the ones that make the two families commute — the neighbouring
    /// choice does not, which the `css_commutation` test would catch.
    pub fn surface(d: usize) -> CssCode {
        assert!(d >= 3 && d % 2 == 1, "rotated surface code needs an odd distance >= 3");
        let n = d * d;
        let q = |r: usize, c: usize| r * d + c;
        let mut z_checks: Vec<Vec<usize>> = Vec::new();
        let mut x_stabs: Vec<Vec<usize>> = Vec::new();

        // Bulk: every 2x2 square, checkerboarded.
        for r in 0..d - 1 {
            for c in 0..d - 1 {
                let face = vec![q(r, c), q(r, c + 1), q(r + 1, c), q(r + 1, c + 1)];
                if (r + c) % 2 == 0 {
                    z_checks.push(face);
                } else {
                    x_stabs.push(face);
                }
            }
        }
        // X boundary: top row on even columns, bottom row on odd columns.
        for c in (0..d - 1).step_by(2) {
            x_stabs.push(vec![q(0, c), q(0, c + 1)]);
        }
        for c in (1..d - 1).step_by(2) {
            x_stabs.push(vec![q(d - 1, c), q(d - 1, c + 1)]);
        }
        // Z boundary: left column on odd rows, right column on even rows.
        for r in (1..d - 1).step_by(2) {
            z_checks.push(vec![q(r, 0), q(r + 1, 0)]);
        }
        for r in (0..d - 1).step_by(2) {
            z_checks.push(vec![q(r, d - 1), q(r + 1, d - 1)]);
        }
        for v in z_checks.iter_mut().chain(x_stabs.iter_mut()) {
            v.sort_unstable();
            v.dedup();
        }
        CssCode::new(format!("surface:d{d}"), n, z_checks, x_stabs)
    }

    // --- Colour codes ------------------------------------------------------

    /// The `[[7,1,3]]` Steane code — the smallest 2D **colour code**.
    ///
    /// Its three checks are the rows of the Hamming `[7,4,3]` parity-check matrix,
    /// used for *both* the X and the Z family. That coincidence is the whole point
    /// of a colour code: the two stabilizer families sit on the same supports, so
    /// the code is **self-dual** as a CSS code.
    ///
    /// The payoff is [`CssCode::has_transversal_hadamard`]. A surface code has to
    /// implement a logical Hadamard with lattice surgery or code deformation; a
    /// colour code gets it by applying `H` to every physical qubit, one round, no
    /// ancillas. The trade is a lower threshold — this buys gates, not noise
    /// tolerance.
    pub fn color_steane() -> CssCode {
        let faces = vec![vec![3, 4, 5, 6], vec![1, 2, 5, 6], vec![0, 2, 4, 6]];
        CssCode::new("color:steane".into(), 7, faces.clone(), faces)
    }

    /// Whether the X and Z stabilizer families generate the **same** group.
    ///
    /// This is the algebraic content of "colour code". `H⊗n` exchanges X-type and
    /// Z-type Paulis, so it maps the stabilizer group to itself exactly when the
    /// two families coincide — and a map that preserves the stabilizer group
    /// preserves the code space.
    pub fn is_self_dual(&self) -> bool {
        let a = Gf2Basis::from_rows(&self.checks, self.n);
        let b = Gf2Basis::from_rows(&self.stabs, self.n);
        if a.rank() != b.rank() {
            return false;
        }
        // same rank and each family reduces to nothing against the other ⇒ same span
        self.checks.iter().all(|r| b.reduces_to_zero(r, self.n))
            && self.stabs.iter().all(|r| a.reduces_to_zero(r, self.n))
    }

    /// Whether applying `H` to every physical qubit is a logical operation.
    ///
    /// True exactly for a self-dual CSS code. Stated as its own question because
    /// it is the property people actually want from a colour code.
    pub fn has_transversal_hadamard(&self) -> bool {
        self.is_self_dual()
    }

    /// A circuit preparing the logical `|0⟩` of this CSS code from `|0…0⟩`.
    ///
    /// `|0⟩_L` is the uniform superposition over the row space of the X-stabilizer
    /// matrix. Reduce that matrix to *fully* reduced row-echelon form, put an `H`
    /// on each pivot qubit, and fan a `CNOT` out to the rest of its row: each row
    /// then contributes an independent branch, which is exactly that superposition.
    ///
    /// The full reduction matters — with pivots shared between rows the fan-out
    /// would interfere with itself and prepare something else.
    ///
    /// Needs the `quantum` feature: the decoding layer above is deliberately
    /// independent of the simulator, and this is the one place that borrows it.
    #[cfg(feature = "quantum")]
    pub fn encode_zero_circuit(&self) -> crate::quantum::Circuit {
        let mut c = crate::quantum::Circuit::new(self.n as u8);
        for (pivot, cols) in rref(&self.stabs, self.n) {
            c.h(pivot as u8);
            for col in cols {
                if col != pivot {
                    c.cx(pivot as u8, col as u8);
                }
            }
        }
        c
    }

    // --- Bivariate-bicycle qLDPC code -------------------------------------

    /// A bivariate-bicycle code on `Z_l × Z_m`. `a` and `b` are lists of monomials
    /// `(px, py)`; `A = Σ x^px y^py`, `B` likewise (mod-2 sums of shift matrices).
    /// `H_X = [A | B]`, `H_Z = [Bᵀ | Aᵀ]`, `n = 2·l·m`. We decode X-errors:
    /// `checks = H_Z`, triviality against `H_X`.
    pub fn bivariate_bicycle(l: usize, m: usize, a: &[(usize, usize)], b: &[(usize, usize)]) -> CssCode {
        let lm = l * m;
        let n = 2 * lm;
        let cell = |i: usize, j: usize| (i % l) * m + (j % m);
        // A[row][col]: monomial (px,py) sends col-cell (i,j) → row-cell (i+px, j+py).
        // Build row-wise supports for A and B (over the lm cells).
        let build = |monos: &[(usize, usize)]| -> Vec<Vec<usize>> {
            let mut rows = vec![Vec::new(); lm];
            for i in 0..l {
                for j in 0..m {
                    let col = cell(i, j);
                    for &(px, py) in monos {
                        let row = cell(i + px, j + py);
                        rows[row].push(col);
                    }
                }
            }
            for r in rows.iter_mut() {
                r.sort_unstable();
                r.dedup();
            }
            rows
        };
        let a_rows = build(a);
        let b_rows = build(b);
        // transpose helper (for Bᵀ, Aᵀ)
        let transpose = |rows: &[Vec<usize>]| -> Vec<Vec<usize>> {
            let mut t = vec![Vec::new(); lm];
            for (r, cols) in rows.iter().enumerate() {
                for &c in cols {
                    t[c].push(r);
                }
            }
            for x in t.iter_mut() {
                x.sort_unstable();
            }
            t
        };
        let at = transpose(&a_rows);
        let bt = transpose(&b_rows);

        // H_X = [A | B]  (columns 0..lm = block L, lm..2lm = block R)
        let hx: Vec<Vec<usize>> = (0..lm)
            .map(|r| {
                let mut row: Vec<usize> = a_rows[r].clone();
                row.extend(b_rows[r].iter().map(|&c| c + lm));
                row.sort_unstable();
                row
            })
            .collect();
        // H_Z = [Bᵀ | Aᵀ]
        let hz: Vec<Vec<usize>> = (0..lm)
            .map(|r| {
                let mut row: Vec<usize> = bt[r].clone();
                row.extend(at[r].iter().map(|&c| c + lm));
                row.sort_unstable();
                row
            })
            .collect();

        // decode X-errors: checks = H_Z, triviality against H_X
        let code = CssCode::new(String::new(), n, hz, hx);
        CssCode { id: format!("bb:[[{},{}]]", code.n, code.k), ..code }
    }

    /// The `[[72,12]]` bivariate-bicycle preset (l=m=6, A=x³+y+y², B=y³+x+x²).
    pub fn gross() -> CssCode {
        CssCode::bivariate_bicycle(6, 6, &[(3, 0), (0, 1), (0, 2)], &[(0, 3), (1, 0), (2, 0)])
    }
}

// ===========================================================================
// Tanner graph + normalized min-sum belief propagation
// ===========================================================================

struct Tanner {
    n: usize,
    n_edges: usize,
    e_var: Vec<u32>,          // variable (qubit) of each edge
    chk_range: Vec<(u32, u32)>, // contiguous edge range per check
    var_edges: Vec<Vec<u32>>, // edges incident to each variable
}

impl Tanner {
    fn build(checks: &[Vec<usize>], n: usize) -> Tanner {
        let mut e_var = Vec::new();
        let mut chk_range = Vec::with_capacity(checks.len());
        let mut var_edges = vec![Vec::new(); n];
        for cols in checks {
            let start = e_var.len() as u32;
            for &v in cols {
                var_edges[v].push(e_var.len() as u32);
                e_var.push(v as u32);
            }
            chk_range.push((start, e_var.len() as u32));
        }
        let n_edges = e_var.len();
        Tanner { n, n_edges, e_var, chk_range, var_edges }
    }
}

/// One min-sum BP leg over the current messages. `gamma` is the per-variable relay
/// memory strength (fixed-point, 0 = no memory). Returns `(hard_decision,
/// converged, iters, messages)`. Messages are warm-started from `mu_v2c`.
#[allow(clippy::too_many_arguments)]
fn bp_leg(
    t: &Tanner,
    syndrome: &[bool],
    gamma: &[i64],
    max_iter: u32,
    mu_v2c: &mut [i64],
    mu_c2v: &mut [i64],
) -> (Vec<bool>, bool, u32, u64) {
    let mut hard = vec![false; t.n];
    let mut messages = 0u64;
    let mut iters = 0;
    for _ in 0..max_iter {
        iters += 1;
        // --- check update (min-sum) ---
        for (&(a, b), &s) in t.chk_range.iter().zip(syndrome) {
            let (mut min1, mut min2) = (i64::MAX, i64::MAX);
            let mut arg = a;
            let mut neg_parity = false;
            for e in a..b {
                let x = mu_v2c[e as usize];
                if x < 0 {
                    neg_parity = !neg_parity;
                }
                let mag = x.abs();
                if mag < min1 {
                    min2 = min1;
                    min1 = mag;
                    arg = e;
                } else if mag < min2 {
                    min2 = mag;
                }
            }
            for e in a..b {
                let x = mu_v2c[e as usize];
                let base = if e == arg { min2 } else { min1 };
                let mag = (base * ALPHA_NUM) >> 3;
                let is_neg = neg_parity ^ (x < 0) ^ s;
                mu_c2v[e as usize] = if is_neg { -mag } else { mag };
            }
        }
        // --- variable update + hard decision ---
        for v in 0..t.n {
            let mut total = L0;
            for &e in &t.var_edges[v] {
                total += mu_c2v[e as usize];
            }
            hard[v] = total < 0;
            let g = gamma[v];
            for &e in &t.var_edges[v] {
                let target = total - mu_c2v[e as usize];
                let nv = if g == 0 {
                    target
                } else {
                    (g * mu_v2c[e as usize] + (SCALE - g) * target) >> QEC_FRAC
                };
                mu_v2c[e as usize] = nv.clamp(-MSG_MAX, MSG_MAX);
            }
        }
        messages += 2 * t.n_edges as u64;
        // --- convergence: H·hard == syndrome ---
        let mut ok = true;
        for (&(a, b), &s) in t.chk_range.iter().zip(syndrome) {
            let mut par = false;
            for e in a..b {
                par ^= hard[t.e_var[e as usize] as usize];
            }
            if par != s {
                ok = false;
                break;
            }
        }
        if ok {
            return (hard, true, iters, messages);
        }
    }
    (hard, false, iters, messages)
}

/// Relay-BP: leg 0 is plain min-sum; each subsequent leg draws disordered
/// per-variable memory strengths (deterministic from `seed`) and warm-starts
/// (relays) from the previous leg. Stops at the first leg that satisfies the
/// syndrome. `max_legs = 1` recovers plain BP.
struct DecodeCore {
    correction: Vec<bool>,
    converged: bool,
    legs: u32,
    iters: u32,
    messages: u64,
}

fn relay_bp(
    t: &Tanner,
    syndrome: &[bool],
    seed: u64,
    max_iter: u32,
    max_legs: u32,
) -> DecodeCore {
    let mut mu_v2c = vec![L0; t.n_edges];
    let mut mu_c2v = vec![0i64; t.n_edges];
    let gamma0 = vec![0i64; t.n];
    let (mut hard, mut conv, mut it, mut msg) =
        bp_leg(t, syndrome, &gamma0, max_iter, &mut mu_v2c, &mut mu_c2v);
    let mut legs = 1u32;
    let (mut tot_it, mut tot_msg) = (it, msg);
    while !conv && legs < max_legs {
        let mut gamma = vec![0i64; t.n];
        let mut st = seed
            .wrapping_mul(0x100_0001)
            .wrapping_add(legs as u64)
            .wrapping_add(0xD15EA5E);
        let span = (GAMMA_MAX - GAMMA_MIN) as u64 + 1;
        for gv in gamma.iter_mut() {
            let r = splitmix64(&mut st);
            *gv = GAMMA_MIN + (r % span) as i64;
        }
        let out = bp_leg(t, syndrome, &gamma, max_iter, &mut mu_v2c, &mut mu_c2v);
        hard = out.0;
        conv = out.1;
        it = out.2;
        msg = out.3;
        tot_it += it;
        tot_msg += msg;
        legs += 1;
    }
    let _ = it;
    let _ = msg;
    DecodeCore { correction: hard, converged: conv, legs, iters: tot_it, messages: tot_msg }
}

// ===========================================================================
// Simulation front-end + decode result
// ===========================================================================

/// Deterministic decoder configuration.
#[derive(Clone, Copy, Debug)]
pub struct DecodeConfig {
    /// Physical error rate per qubit (fixed-point, scale 2^12).
    pub p_fx: i64,
    /// Max BP iterations per leg.
    pub max_iter: u32,
    /// Max relay legs (1 = plain BP).
    pub max_legs: u32,
}

impl DecodeConfig {
    pub fn from_p(p: f64) -> DecodeConfig {
        DecodeConfig { p_fx: (p * SCALE as f64) as i64, max_iter: 32, max_legs: 12 }
    }
}

/// One decode: sample an error at rate `p`, decode its syndrome, and check
/// logical success against the true error. Deterministic given `(code, cfg, seed)`.
#[derive(Clone, Debug)]
pub struct DecodeResult {
    pub code_id: String,
    pub n: usize,
    pub error: Vec<u32>,      // true error qubits (simulation only)
    pub correction: Vec<u32>, // decoder output qubits
    pub syndrome: Vec<u32>,   // unsatisfied check indices
    pub converged: bool,
    pub logical_success: bool,
    pub legs: u32,
    pub iters: u32,
    pub messages: u64,
    syndrome_bits: Vec<bool>,
    correction_bits: Vec<bool>,
    check_matrix_hash: [u8; 32],
}

/// Sample a bit-flip error deterministically at rate `p_fx`.
fn sample_error(n: usize, p_fx: i64, seed: u64) -> Vec<bool> {
    let mut e = vec![false; n];
    let mut st = seed.wrapping_mul(0x2545_F491_4F6C_DD1D).wrapping_add(1);
    for ei in e.iter_mut() {
        let r = (splitmix64(&mut st) >> (64 - QEC_FRAC)) as i64; // [0, SCALE)
        if r < p_fx {
            *ei = true;
        }
    }
    e
}

fn syndrome_of(checks: &[Vec<usize>], e: &[bool]) -> Vec<bool> {
    checks
        .iter()
        .map(|row| row.iter().fold(false, |acc, &q| acc ^ e[q]))
        .collect()
}

fn bits_to_ids(bits: &[bool]) -> Vec<u32> {
    bits.iter()
        .enumerate()
        .filter(|&(_, &b)| b)
        .map(|(i, _)| i as u32)
        .collect()
}

fn canonical_bits(tag: &[u8], bits: &[bool]) -> Vec<u8> {
    let mut out = tag.to_vec();
    out.extend_from_slice(&(bits.len() as u64).to_le_bytes());
    for chunk in bits.chunks(8) {
        let mut byte = 0u8;
        for (i, &b) in chunk.iter().enumerate() {
            if b {
                byte |= 1 << i;
            }
        }
        out.push(byte);
    }
    out
}

/// Run one full decode on `code` with `cfg`, sampling error stream `seed`.
pub fn decode_once(code: &CssCode, cfg: &DecodeConfig, seed: u64) -> DecodeResult {
    let t = Tanner::build(&code.checks, code.n);
    let error = sample_error(code.n, cfg.p_fx, seed);
    let syndrome_bits = syndrome_of(&code.checks, &error);
    let core = relay_bp(&t, &syndrome_bits, seed, cfg.max_iter, cfg.max_legs);
    // residual = error XOR correction
    let residual: Vec<usize> = (0..code.n)
        .filter(|&i| error[i] ^ core.correction[i])
        .collect();
    let logical_success = core.converged && code.is_trivial(&residual);
    DecodeResult {
        code_id: code.id.clone(),
        n: code.n,
        error: bits_to_ids(&error),
        correction: bits_to_ids(&core.correction),
        syndrome: bits_to_ids(&syndrome_bits),
        converged: core.converged,
        logical_success,
        legs: core.legs,
        iters: core.iters,
        messages: core.messages,
        syndrome_bits,
        correction_bits: core.correction.clone(),
        check_matrix_hash: content_hash(&code.check_matrix_bytes()),
    }
}

impl DecodeResult {
    /// Seal this decode into a signed, joule-metered receipt. `joules_micro` is the
    /// decode's energy: on real hardware a RAPL/NVML meter; here a caller-supplied
    /// modeled figure (the `messages` work is the exact, deterministic basis).
    pub fn seal(
        &self,
        signer: &SigningKey,
        signer_id: impl Into<String>,
        joules_micro: u64,
        grant: GrantRef,
    ) -> DecodeReceipt {
        DecodeReceipt::seal(
            signer,
            signer_id,
            self.code_id.clone(),
            "relay-bp.min-sum",
            self.check_matrix_hash,
            content_hash(&canonical_bits(b"wai:qec-syndrome\x01", &self.syndrome_bits)),
            content_hash(&canonical_bits(b"wai:qec-correction\x01", &self.correction_bits)),
            self.converged,
            self.legs,
            self.messages,
            joules_micro,
            grant,
            None,
        )
    }

    /// Canonical correction bytes (what the receipt's `correction_hash` binds).
    pub fn correction_bytes(&self) -> Vec<u8> {
        canonical_bits(b"wai:qec-correction\x01", &self.correction_bits)
    }
}

/// A Monte-Carlo comparison of plain BP vs Relay-BP on the same error samples —
/// the demonstration that relay is what makes BP work on quantum codes.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct BenchResult {
    pub trials: u32,
    pub plain_logical_failures: u32,
    pub relay_logical_failures: u32,
    pub relay_total_messages: u64,
}

/// Run `trials` decodes at rate `p`, each with plain BP (`max_legs=1`) and
/// Relay-BP, over identical error samples, counting logical failures of each.
/// A matching-style **union-find decoder** (Delfosse–Nickerson).
///
/// Belief propagation struggles on topological codes: they are degenerate and
/// their Tanner graphs are full of short cycles, and the measurements in this
/// module's tests show the consequence — under Relay-BP a distance-5 surface code
/// decodes *worse* than distance-3. Union-find decodes the syndrome the way the
/// code's structure actually wants: as a matching problem on a graph whose nodes
/// are checks and whose edges are qubits.
///
/// Grow every odd cluster until no odd cluster remains, then peel each cluster's
/// spanning forest from the leaves inward to read off a correction. Almost-linear
/// time, and unlike BP it makes the distance count.
///
/// # Applicability
///
/// This only makes sense when a qubit touches at most two checks — true for
/// surface and toric codes, false for the higher-degree qLDPC families (a
/// bivariate-bicycle qubit sits in three). Returns `None` rather than guessing
/// when the code is not a matching graph; use Relay-BP for those.
/// A minimum-weight **logical Z** representative: the observable that decides
/// whether an X-error chain flipped the encoded qubit.
///
/// The Z checks detect X errors, so failure is not "an error remained" but "the
/// residual X chain crossed the code". A logical Z operator is exactly the ruler
/// that measures that: it commutes with every X stabilizer, is not itself a
/// product of Z checks, and its overlap with the residual has odd parity precisely
/// when the logical qubit flipped.
///
/// Returned as the qubit support. Exhaustive by increasing weight, so it is a
/// minimum-weight representative, not just some representative.
pub fn logical_z_support(code: &CssCode, max_weight: usize) -> Option<Vec<usize>> {
    let n = code.n;
    let check_basis = Gf2Basis::from_rows(&code.checks, n);
    for w in 1..=max_weight.min(n) {
        let mut combo: Vec<usize> = (0..w).collect();
        loop {
            let mut bits = vec![false; n];
            for &q in &combo {
                bits[q] = true;
            }
            // commutes with every X stabilizer, and is not a product of Z checks
            if syndrome_of(&code.stabs, &bits).iter().all(|b| !b) && !check_basis.contains(&combo) {
                return Some(combo);
            }
            if !next_combination(&mut combo, n) {
                break;
            }
        }
    }
    None
}

/// Advance `c` to the next `w`-combination of `0..n` in lexicographic order.
/// Returns `false` once the last combination has been passed.
fn next_combination(c: &mut [usize], n: usize) -> bool {
    let w = c.len();
    let mut i = w;
    while i > 0 {
        i -= 1;
        if c[i] < n - (w - i) {
            c[i] += 1;
            for j in i + 1..w {
                c[j] = c[j - 1] + 1;
            }
            return true;
        }
    }
    false
}

/// The code's **X-distance**: the least weight of an undetectable error that is
/// not a stabilizer — a logical operator.
///
/// This is the `d` in `[[n, k, d]]`, and it is the claim a code family actually
/// rests on: any error of weight below `d` is correctable, because no error that
/// small can be logical. Asserting `k` while leaving `d` unchecked verifies the
/// less interesting half of the name.
///
/// Searched exactly by increasing weight, so a returned `d` is a proof: every
/// lighter operator has been ruled out by exhaustion. `None` means no logical of
/// weight `≤ max_weight` exists. The search is exponential in the cap — `d ≤ 5`
/// is instant, `d = 7` is not.
pub fn x_distance(code: &CssCode, max_weight: usize) -> Option<usize> {
    let n = code.n;
    for w in 1..=max_weight.min(n) {
        let mut combo: Vec<usize> = (0..w).collect();
        loop {
            let mut bits = vec![false; n];
            for &q in &combo {
                bits[q] = true;
            }
            // undetectable (zero syndrome) and not a stabilizer => logical
            if syndrome_of(&code.checks, &bits).iter().all(|b| !b) && !code.is_trivial(&combo) {
                return Some(w);
            }
            if !next_combination(&mut combo, n) {
                break;
            }
        }
    }
    None
}

/// The union-find core, on an abstract graph.
///
/// `edges` are `(u, v)` node pairs; `boundary` is the node that may absorb any
/// leftover parity (open boundaries, and the time boundary of a space-time graph).
/// Returns which edges the decoder selected.
///
/// Grow every odd cluster until none remain, then peel each cluster's spanning
/// forest from the leaves inward. Growth adds a **half-edge from each odd
/// endpoint**, so an edge joining two syndromes closes in one round while an edge
/// to a neutral node needs two — that asymmetry is what makes the decoder pair
/// syndromes with each other rather than draining them all to the boundary.
///
/// `weights` are per-edge costs: an edge closes at `2 × weight`, so cheap edges
/// close first. Uniform weights recover plain union-find; deriving them from the
/// error rates is what stops the decoder explaining a lying measurement with a
/// data correction when data errors are the less likely cause.
fn uf_decode_graph(
    nodes: usize,
    boundary: usize,
    edges: &[(usize, usize)],
    weights: &[u32],
    syndrome: &[bool],
) -> Vec<bool> {
    let mut parent: Vec<usize> = (0..nodes).collect();
    fn find(parent: &mut Vec<usize>, mut x: usize) -> usize {
        while parent[x] != x {
            parent[x] = parent[parent[x]];
            x = parent[x];
        }
        x
    }
    let mut syn: Vec<bool> = (0..nodes).map(|i| syndrome.get(i).copied().unwrap_or(false)).collect();
    syn[boundary] = false;
    let mut grown: Vec<u32> = vec![0; edges.len()];
    let maxw = weights.iter().copied().max().unwrap_or(1).max(1) as usize;

    for _ in 0..(4 * nodes * maxw + 8) {
        let mut parity = vec![false; nodes];
        let mut touches = vec![false; nodes];
        for v in 0..nodes {
            let r = find(&mut parent, v);
            if syn[v] {
                parity[r] ^= true;
            }
            if v == boundary {
                touches[r] = true;
            }
        }
        if !(0..nodes).any(|v| find(&mut parent, v) == v && parity[v] && !touches[v]) {
            break;
        }
        let mut to_union: Vec<usize> = Vec::new();
        for (ei, &(a, b)) in edges.iter().enumerate() {
            if grown[ei] >= 2 * weights.get(ei).copied().unwrap_or(1).max(1) {
                continue;
            }
            let (ra, rb) = (find(&mut parent, a), find(&mut parent, b));
            let odd_a = parity[ra] && !touches[ra];
            let odd_b = parity[rb] && !touches[rb];
            let inc = if ra == rb { u8::from(odd_a) } else { u8::from(odd_a) + u8::from(odd_b) };
            if inc > 0 {
                let cap = 2 * weights.get(ei).copied().unwrap_or(1).max(1);
                grown[ei] = grown[ei].saturating_add(inc as u32);
                if grown[ei] >= cap {
                    to_union.push(ei);
                }
            }
        }
        for ei in to_union {
            let (a, b) = edges[ei];
            let (ra, rb) = (find(&mut parent, a), find(&mut parent, b));
            if ra != rb {
                parent[ra] = rb;
            }
        }
    }

    let mut adj: Vec<Vec<(usize, usize)>> = vec![Vec::new(); nodes];
    for (ei, &(a, b)) in edges.iter().enumerate() {
        if grown[ei] >= 2 * weights.get(ei).copied().unwrap_or(1).max(1) {
            adj[a].push((b, ei));
            adj[b].push((a, ei));
        }
    }
    let mut picked = vec![false; edges.len()];
    let mut seen = vec![false; nodes];
    for root in std::iter::once(boundary).chain(0..nodes) {
        if seen[root] {
            continue;
        }
        let mut order = vec![root];
        let mut parent_edge: Vec<Option<(usize, usize)>> = vec![None; nodes];
        seen[root] = true;
        let mut i = 0;
        while i < order.len() {
            let v = order[i];
            i += 1;
            for &(w, ei) in &adj[v] {
                if !seen[w] {
                    seen[w] = true;
                    parent_edge[w] = Some((v, ei));
                    order.push(w);
                }
            }
        }
        for &v in order.iter().rev() {
            if v == root || !syn[v] {
                continue;
            }
            if let Some((u, ei)) = parent_edge[v] {
                picked[ei] ^= true;
                syn[v] = false;
                syn[u] ^= true;
            }
        }
    }
    picked
}

/// A matching-style **union-find decoder** (Delfosse–Nickerson) for one shot of
/// syndrome, assuming perfect measurement.
///
/// Belief propagation struggles on topological codes: they are degenerate and
/// their Tanner graphs are full of short cycles. Union-find decodes the syndrome
/// the way the code's structure wants — as matching on a graph whose nodes are
/// checks and whose edges are qubits — and unlike BP it makes the distance count.
///
/// # Applicability
///
/// Only meaningful when a qubit touches at most two checks: true for surface and
/// toric codes, false for higher-degree qLDPC (a bivariate-bicycle qubit sits in
/// three). Returns `None` rather than guessing; use Relay-BP for those.
pub fn union_find_correct(code: &CssCode, syndrome: &[bool]) -> Option<Vec<bool>> {
    let m = code.checks.len();
    let boundary = m;
    let mut owners: Vec<Vec<usize>> = vec![Vec::new(); code.n];
    for (ci, chk) in code.checks.iter().enumerate() {
        for &q in chk {
            owners[q].push(ci);
        }
    }
    let mut edges: Vec<(usize, usize)> = Vec::new();
    let mut edge_qubit: Vec<usize> = Vec::new();
    for (q, own) in owners.iter().enumerate() {
        match own.len() {
            0 => {}
            1 => {
                edges.push((own[0], boundary));
                edge_qubit.push(q);
            }
            2 => {
                edges.push((own[0], own[1]));
                edge_qubit.push(q);
            }
            _ => return None,
        }
    }
    let weights = vec![1u32; edges.len()];
    let picked = uf_decode_graph(m + 1, boundary, &edges, &weights, syndrome);
    let mut correction = vec![false; code.n];
    for (ei, &on) in picked.iter().enumerate() {
        if on {
            correction[edge_qubit[ei]] ^= true;
        }
    }
    Some(correction)
}

/// Build the space-time matching graph for a block and decode it to a data-qubit
/// correction. Shared by the phenomenological and circuit-level models so both
/// decode through identical code — only the fault simulation differs.
fn phenom_correction(code: &CssCode, cfg: &PhenomConfig, detector_bits: &[bool]) -> Option<Vec<bool>> {
    let m = code.checks.len();
    let t_rounds = cfg.rounds.max(1) as usize;
    let layers = t_rounds + 1;
    let mut owners: Vec<Vec<usize>> = vec![Vec::new(); code.n];
    for (ci, chk) in code.checks.iter().enumerate() {
        for &q in chk {
            owners[q].push(ci);
        }
    }
    if owners.iter().any(|o| o.len() > 2) {
        return None;
    }
    // --- space-time matching graph ---
    let boundary = layers * m;
    let node = |c: usize, t: usize| t * m + c;
    // Cost of an edge is -ln(rate): the less likely a fault, the more the decoder
    // must pay to blame it. A rate of zero means the fault cannot happen, so the
    // edge is omitted rather than merely made expensive.
    let cost = |rate_fx: i64| -> u32 {
        let r = rate_fx as f64 / 4096.0;
        if r <= 0.0 { 0 } else { (-r.ln()).round().clamp(1.0, 12.0) as u32 }
    };
    let w_space = cost(cfg.p_data_fx);
    let w_time = cost(cfg.q_meas_fx);
    let mut edges: Vec<(usize, usize)> = Vec::new();
    let mut weights: Vec<u32> = Vec::new();
    // space edges: a data error at layer t lights its two checks at layer t
    let mut space_qubit: Vec<Option<usize>> = Vec::new();
    if w_space > 0 {
        for t in 0..layers {
            for (q, own) in owners.iter().enumerate() {
                match own.len() {
                    1 => {
                        edges.push((node(own[0], t), boundary));
                        weights.push(w_space);
                        space_qubit.push(Some(q));
                    }
                    2 => {
                        edges.push((node(own[0], t), node(own[1], t)));
                        weights.push(w_space);
                        space_qubit.push(Some(q));
                    }
                    _ => {}
                }
            }
        }
    }
    // time edges: a measurement error at round t lights the same check at t and t+1
    if w_time > 0 {
        for t in 0..layers - 1 {
            for c in 0..m {
                edges.push((node(c, t), node(c, t + 1)));
                weights.push(w_time);
                space_qubit.push(None);
            }
        }
    }

    let picked = uf_decode_graph(boundary + 1, boundary, &edges, &weights, detector_bits);
    let mut correction = vec![false; code.n];
    for (ei, &on) in picked.iter().enumerate() {
        if on && let Some(q) = space_qubit[ei] {
            correction[q] ^= true;
        }
    }

    Some(correction)
}

/// Repeated-round noise with **faulty syndrome measurement** — the
/// phenomenological model.
///
/// Single-shot decoding assumes the syndrome itself is read perfectly, which no
/// device does. Here each of `rounds` measurement rounds can misreport a check,
/// so a lone flipped bit is ambiguous: a data error, or a lying measurement?
/// The resolution is to decode in **space-time**. Detectors are *differences*
/// between consecutive rounds, so a measurement error lights two detectors in the
/// same check at adjacent times, while a data error lights two checks at the same
/// time. A final perfect round closes the block.
#[derive(Clone, Copy, Debug)]
pub struct PhenomConfig {
    /// Noisy measurement rounds (a perfect readout round is appended).
    pub rounds: u32,
    /// Data error per qubit per round (fixed-point, scale 2^12).
    pub p_data_fx: i64,
    /// Measurement error per check per noisy round (fixed-point, scale 2^12).
    pub q_meas_fx: i64,
}

impl PhenomConfig {
    pub fn from_p(rounds: u32, p_data: f64, q_meas: f64) -> PhenomConfig {
        PhenomConfig {
            rounds,
            p_data_fx: (p_data * 4096.0) as i64,
            q_meas_fx: (q_meas * 4096.0) as i64,
        }
    }
}

/// One space-time decode.
#[derive(Clone, Debug)]
pub struct PhenomResult {
    pub rounds: u32,
    /// Detectors that fired, as `round * n_checks + check`.
    pub detectors: Vec<u32>,
    /// Accumulated data error over the whole block.
    pub data_error: Vec<u32>,
    pub correction: Vec<u32>,
    pub logical_success: bool,
}

/// Decode a phenomenological-noise block with union-find over the space-time
/// graph. Returns `None` when the code is not a matching graph.
pub fn decode_phenomenological(
    code: &CssCode,
    cfg: &PhenomConfig,
    seed: u64,
) -> Option<PhenomResult> {
    let m = code.checks.len();
    let t_rounds = cfg.rounds.max(1) as usize;
    let layers = t_rounds + 1; // + the final perfect readout

    // qubit -> its checks (matching graph precondition)
    let mut owners: Vec<Vec<usize>> = vec![Vec::new(); code.n];
    for (ci, chk) in code.checks.iter().enumerate() {
        for &q in chk {
            owners[q].push(ci);
        }
    }
    if owners.iter().any(|o| o.len() > 2) {
        return None;
    }

    // --- simulate the block ---
    let mut st = seed ^ 0x51ED_2701_A55A_1234;
    let draw = |rate: i64, st: &mut u64| -> bool {
        *st = st
            .wrapping_mul(6364136223846793005)
            .wrapping_add(1442695040888963407);
        (((*st >> 40) & 0xFFF) as i64) < rate
    };
    let mut cum_err = vec![false; code.n];
    let mut s_prev = vec![false; m];
    let mut detector_bits = vec![false; layers * m];
    for t in 0..layers {
        if t < t_rounds {
            for q in 0..code.n {
                if draw(cfg.p_data_fx, &mut st) {
                    cum_err[q] ^= true;
                }
            }
        }
        let s_true = syndrome_of(&code.checks, &cum_err);
        let mut s_obs = s_true.clone();
        if t < t_rounds {
            for c in 0..m {
                if draw(cfg.q_meas_fx, &mut st) {
                    s_obs[c] ^= true;
                }
            }
        }
        for c in 0..m {
            detector_bits[t * m + c] = s_obs[c] ^ s_prev[c];
        }
        s_prev = s_obs;
    }

    let correction = phenom_correction(code, cfg, &detector_bits)?;
    let residual: Vec<usize> = (0..code.n).filter(|&i| cum_err[i] ^ correction[i]).collect();
    Some(PhenomResult {
        rounds: cfg.rounds,
        detectors: bits_to_ids(&detector_bits),
        data_error: bits_to_ids(&cum_err),
        correction: bits_to_ids(&correction),
        logical_success: code.is_trivial(&residual),
    })
}

/// Circuit-level noise: syndrome extraction that can fault mid-circuit.
///
/// Phenomenological noise treats a measurement as a black box that occasionally
/// lies. A real extraction round is a *circuit* — an ancilla and a sequence of
/// CNOTs — and a fault inside it does not stay put. A fault on a stabilizer's
/// ancilla propagates through every CNOT still to come, landing as X errors on
/// all the data qubits that circuit has not touched yet: a **hook error**.
///
/// Hook errors are why circuit-level thresholds sit well below phenomenological
/// ones. They are *correlated* — one fault, several data qubits, always adjacent
/// in the code — and a matching graph has no edge for them, so the decoder must
/// explain one fault with several, and a weight-`⌈d/2⌉` hook can already reach
/// across a distance-`d` code.
///
/// # Scope, stated plainly
///
/// This models the fault mechanism that distinguishes circuit-level noise —
/// ancilla faults propagating along a fixed CNOT order — on top of the data and
/// measurement faults of [`decode_phenomenological`]. It is not a full
/// gate-by-gate simulation of both stabilizer families with every propagation
/// path; `hook_fx = 0` reduces it exactly to the phenomenological model.
#[derive(Clone, Copy, Debug)]
pub struct CircuitConfig {
    pub rounds: u32,
    /// Data error per qubit per round.
    pub p_data_fx: i64,
    /// Ancilla measurement flip per check per round.
    pub q_meas_fx: i64,
    /// Ancilla fault per stabilizer per round, propagating as a hook error.
    pub hook_fx: i64,
}

impl CircuitConfig {
    pub fn from_p(rounds: u32, p_data: f64, q_meas: f64, hook: f64) -> CircuitConfig {
        CircuitConfig {
            rounds,
            p_data_fx: (p_data * 4096.0) as i64,
            q_meas_fx: (q_meas * 4096.0) as i64,
            hook_fx: (hook * 4096.0) as i64,
        }
    }
    fn phenom(&self) -> PhenomConfig {
        PhenomConfig { rounds: self.rounds, p_data_fx: self.p_data_fx, q_meas_fx: self.q_meas_fx }
    }
}

/// One circuit-level decode. `hook_events` records how many ancilla faults fired
/// and the data weight each landed as — the mechanism, visible rather than buried.
#[derive(Clone, Debug)]
pub struct CircuitResult {
    pub rounds: u32,
    pub data_error: Vec<u32>,
    pub correction: Vec<u32>,
    pub hook_events: Vec<u32>,
    pub logical_success: bool,
}

/// Decode a block under circuit-level noise. `None` when the code is not a
/// matching graph.
pub fn decode_circuit_level(
    code: &CssCode,
    cfg: &CircuitConfig,
    seed: u64,
) -> Option<CircuitResult> {
    let m = code.checks.len();
    let t_rounds = cfg.rounds.max(1) as usize;
    let layers = t_rounds + 1;

    let mut owners: Vec<Vec<usize>> = vec![Vec::new(); code.n];
    for (ci, chk) in code.checks.iter().enumerate() {
        for &q in chk {
            owners[q].push(ci);
        }
    }
    if owners.iter().any(|o| o.len() > 2) {
        return None;
    }

    // Same stream as decode_phenomenological, and the hook draws below are
    // guarded, so hook_fx = 0 reproduces that model bit-for-bit rather than merely
    // resembling it — a claim worth being able to test.
    let mut st = seed ^ 0x51ED_2701_A55A_1234;
    let draw = |rate: i64, st: &mut u64| -> bool {
        *st = st.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
        (((*st >> 40) & 0xFFF) as i64) < rate
    };
    let pick = |bound: usize, st: &mut u64| -> usize {
        *st = st.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
        if bound == 0 { 0 } else { ((*st >> 33) as usize) % bound }
    };

    let mut cum_err = vec![false; code.n];
    let mut s_prev = vec![false; m];
    let mut detector_bits = vec![false; layers * m];
    let mut hook_events: Vec<u32> = Vec::new();

    for t in 0..layers {
        if t < t_rounds {
            for q in 0..code.n {
                if draw(cfg.p_data_fx, &mut st) {
                    cum_err[q] ^= true;
                }
            }
            // Ancilla faults: propagate along the remaining CNOTs of that
            // stabilizer's extraction circuit.
            if cfg.hook_fx > 0 {
                for stab in &code.stabs {
                    if !stab.is_empty() && draw(cfg.hook_fx, &mut st) {
                        let j = pick(stab.len(), &mut st);
                        for &q in &stab[j..] {
                            cum_err[q] ^= true;
                        }
                        hook_events.push((stab.len() - j) as u32);
                    }
                }
            }
        }
        let s_true = syndrome_of(&code.checks, &cum_err);
        let mut s_obs = s_true.clone();
        if t < t_rounds {
            for c in 0..m {
                if draw(cfg.q_meas_fx, &mut st) {
                    s_obs[c] ^= true;
                }
            }
        }
        for c in 0..m {
            detector_bits[t * m + c] = s_obs[c] ^ s_prev[c];
        }
        s_prev = s_obs;
    }

    let correction = phenom_correction(code, &cfg.phenom(), &detector_bits)?;
    let residual: Vec<usize> = (0..code.n).filter(|&i| cum_err[i] ^ correction[i]).collect();
    Some(CircuitResult {
        rounds: cfg.rounds,
        data_error: bits_to_ids(&cum_err),
        correction: bits_to_ids(&correction),
        hook_events,
        logical_success: code.is_trivial(&residual),
    })
}

/// One independent error mechanism of a detector error model.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct DemMechanism {
    /// Detectors flipped, as `round * n_checks + check`, sorted.
    pub detectors: Vec<u32>,
    /// Whether the mechanism flips the logical observable.
    pub flips_observable: bool,
    /// Probability, in parts per million (fixed point keeps the model portable).
    pub probability_ppm: u32,
}

/// Export the noise model as a **detector error model** in Stim's DEM syntax.
///
/// A DEM says nothing about qubits or gates. It says: here are the independent
/// ways this experiment can fault, here is which detectors each one lights, and
/// here is whether it flips the logical observable. That is precisely the
/// interface every matching decoder wants, which is why it is the interchange
/// format of the QEC ecosystem — the analogue of what QIR is for circuits.
///
/// Emits `detector(check, round) D#` coordinates and one `error(p) D… L0` line per
/// mechanism: data errors, measurement errors, and — when `hook_fx > 0` — one line
/// per ancilla-fault split.
///
/// Whether a hook becomes a *hyperedge* depends on the CNOT order, which is real
/// physics rather than an implementation detail: scheduling is how surface-code
/// circuits control hook damage. With the index order used here, a partial hook on
/// the rotated surface code lights at most two detectors, so the emitted model is
/// an ordinary graph. The hook still costs accuracy — it flips several *data*
/// qubits while presenting as one edge, and that mismatch is the damage — but a
/// consumer of this DEM should not assume hyperedges appear.
///
/// Returns `None` when the code is not a matching graph.
pub fn to_dem(code: &CssCode, cfg: &CircuitConfig, max_logical_weight: usize) -> Option<String> {
    let mechs = dem_mechanisms(code, cfg, max_logical_weight)?;
    let m = code.checks.len();
    let layers = cfg.rounds.max(1) as usize + 1;
    let mut s = String::new();
    s.push_str("# detector error model emitted by wai-quantum\n");
    s.push_str(&format!("# code={} rounds={} detectors={}\n", code.id, cfg.rounds, layers * m));
    for t in 0..layers {
        for c in 0..m {
            s.push_str(&format!("detector({c}, {t}) D{}\n", t * m + c));
        }
    }
    for mech in &mechs {
        let p = mech.probability_ppm as f64 / 1.0e6;
        let mut line = format!("error({p:.9})");
        for d in &mech.detectors {
            line.push_str(&format!(" D{d}"));
        }
        if mech.flips_observable {
            line.push_str(" L0");
        }
        s.push('\n');
        s.truncate(s.len() - 1);
        s.push_str(&line);
        s.push('\n');
    }
    Some(s)
}

/// The mechanisms behind [`to_dem`], as data — so they can be compared against
/// the simulator rather than eyeballed as text.
pub fn dem_mechanisms(
    code: &CssCode,
    cfg: &CircuitConfig,
    max_logical_weight: usize,
) -> Option<Vec<DemMechanism>> {
    let m = code.checks.len();
    let t_rounds = cfg.rounds.max(1) as usize;

    let mut owners: Vec<Vec<usize>> = vec![Vec::new(); code.n];
    for (ci, chk) in code.checks.iter().enumerate() {
        for &q in chk {
            owners[q].push(ci);
        }
    }
    if owners.iter().any(|o| o.len() > 2) {
        return None;
    }
    let logical = logical_z_support(code, max_logical_weight)?;
    let flips = |qs: &[usize]| qs.iter().filter(|q| logical.contains(q)).count() % 2 == 1;
    let ppm = |fx: i64| ((fx as f64 / 4096.0) * 1.0e6).round().max(0.0) as u32;

    let mut out: Vec<DemMechanism> = Vec::new();
    // data errors: light this qubit's checks, at that round only
    if cfg.p_data_fx > 0 {
        for t in 0..t_rounds {
            for q in 0..code.n {
                let mut ds: Vec<u32> = owners[q].iter().map(|&c| (t * m + c) as u32).collect();
                ds.sort_unstable();
                out.push(DemMechanism {
                    detectors: ds,
                    flips_observable: flips(&[q]),
                    probability_ppm: ppm(cfg.p_data_fx),
                });
            }
        }
    }
    // measurement errors: same check, adjacent rounds; never touch the observable
    if cfg.q_meas_fx > 0 {
        for t in 0..t_rounds {
            for c in 0..m {
                out.push(DemMechanism {
                    detectors: vec![(t * m + c) as u32, ((t + 1) * m + c) as u32],
                    flips_observable: false,
                    probability_ppm: ppm(cfg.q_meas_fx),
                });
            }
        }
    }
    // hook errors: one ancilla fault, several data qubits — a hyperedge
    if cfg.hook_fx > 0 {
        for t in 0..t_rounds {
            for stab in &code.stabs {
                if stab.is_empty() {
                    continue;
                }
                let each = ppm(cfg.hook_fx) / stab.len() as u32;
                for j in 0..stab.len() {
                    let qs = &stab[j..];
                    let mut set: Vec<u32> = Vec::new();
                    for &q in qs {
                        for &c in &owners[q] {
                            let d = (t * m + c) as u32;
                            if let Some(pos) = set.iter().position(|&x| x == d) {
                                set.remove(pos); // XOR: two flips cancel
                            } else {
                                set.push(d);
                            }
                        }
                    }
                    set.sort_unstable();
                    // A hook that covers the whole stabilizer IS that stabilizer:
                    // undetectable and harmless. Emitting it would put a target-less
                    // `error(p)` in the DEM, which is not a mechanism at all.
                    let obs = flips(qs);
                    if set.is_empty() && !obs {
                        continue;
                    }
                    out.push(DemMechanism {
                        detectors: set,
                        flips_observable: obs,
                        probability_ppm: each,
                    });
                }
            }
        }
    }
    Some(out)
}

/// A parsed detector error model.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct DemModel {
    pub num_detectors: usize,
    pub mechanisms: Vec<DemMechanism>,
}

/// Parse a detector error model in Stim's DEM syntax.
///
/// Accepts what [`to_dem`] emits and, more usefully, a DEM produced elsewhere —
/// the point of a format is that both directions work. `detector(...)` lines carry
/// only coordinates and are read for the detector count; `error(p) D… L…` lines
/// are the mechanisms. Unknown instructions are refused rather than skipped: a DEM
/// whose `repeat` block was silently ignored would decode a different experiment
/// from the one it describes.
pub fn from_dem(src: &str) -> Result<DemModel, String> {
    let mut mechanisms = Vec::new();
    let mut max_det: i64 = -1;
    for (lineno, raw) in src.lines().enumerate() {
        let line = raw.split('#').next().unwrap_or("").trim();
        if line.is_empty() {
            continue;
        }
        let head = line.split(['(', ' ']).next().unwrap_or("");
        match head {
            "detector" | "shift_detectors" | "logical_observable" | "observable_include" => {
                for tok in line.split_whitespace() {
                    if let Some(d) = tok.strip_prefix('D')
                        && let Ok(i) = d.trim_end_matches(',').parse::<i64>()
                    {
                        max_det = max_det.max(i);
                    }
                }
            }
            "error" => {
                let open = line.find('(').ok_or_else(|| format!("line {}: error without probability", lineno + 1))?;
                let close = line[open..].find(')').ok_or_else(|| format!("line {}: unclosed probability", lineno + 1))? + open;
                let p: f64 = line[open + 1..close]
                    .trim()
                    .parse()
                    .map_err(|_| format!("line {}: bad probability", lineno + 1))?;
                if !(0.0..=1.0).contains(&p) {
                    return Err(format!("line {}: probability {p} outside [0,1]", lineno + 1));
                }
                let mut detectors = Vec::new();
                let mut flips = false;
                for tok in line[close + 1..].split_whitespace() {
                    if let Some(d) = tok.strip_prefix('D') {
                        let i: u32 = d.parse().map_err(|_| format!("line {}: bad detector {tok}", lineno + 1))?;
                        max_det = max_det.max(i as i64);
                        detectors.push(i);
                    } else if let Some(l) = tok.strip_prefix('L') {
                        let _: u32 = l.parse().map_err(|_| format!("line {}: bad observable {tok}", lineno + 1))?;
                        flips = !flips;
                    } else {
                        return Err(format!("line {}: unexpected target {tok}", lineno + 1));
                    }
                }
                detectors.sort_unstable();
                if detectors.is_empty() && !flips {
                    continue; // a mechanism that does nothing
                }
                mechanisms.push(DemMechanism {
                    detectors,
                    flips_observable: flips,
                    probability_ppm: (p * 1.0e6).round().clamp(0.0, 1.0e6) as u32,
                });
            }
            other => return Err(format!("line {}: unsupported DEM instruction `{other}`", lineno + 1)),
        }
    }
    Ok(DemModel { num_detectors: (max_det + 1).max(0) as usize, mechanisms })
}

/// What a DEM decode concluded.
#[derive(Clone, Debug)]
pub struct DemDecode {
    /// Indices of the mechanisms the decoder blamed.
    pub mechanisms: Vec<u32>,
    /// Whether the logical observable is predicted to have flipped.
    pub observable_flipped: bool,
}

/// Decode a syndrome against a detector error model.
///
/// This is the whole point of the format: given somebody else's noise model and a
/// round of detector data, say whether the logical observable flipped. The
/// mechanisms become graph edges — one detector means an edge to the boundary, two
/// means an ordinary edge — weighted by `−ln p`, and union-find matches them.
///
/// Mechanisms touching three or more detectors are hyperedges, which matching
/// cannot represent. Rather than drop them and quietly decode a different model,
/// this refuses.
pub fn decode_dem(model: &DemModel, detectors: &[bool]) -> Result<DemDecode, String> {
    let nd = model.num_detectors;
    let boundary = nd;
    let mut edges: Vec<(usize, usize)> = Vec::new();
    let mut weights: Vec<u32> = Vec::new();
    let mut index: Vec<usize> = Vec::new();
    for (i, m) in model.mechanisms.iter().enumerate() {
        let w = {
            let p = m.probability_ppm as f64 / 1.0e6;
            if p <= 0.0 { continue } else { (-p.ln()).round().clamp(1.0, 12.0) as u32 }
        };
        match m.detectors.len() {
            0 => {}
            1 => {
                edges.push((m.detectors[0] as usize, boundary));
                weights.push(w);
                index.push(i);
            }
            2 => {
                edges.push((m.detectors[0] as usize, m.detectors[1] as usize));
                weights.push(w);
                index.push(i);
            }
            k => {
                return Err(format!(
                    "mechanism {i} touches {k} detectors; matching needs at most two — \
                     decompose the hyperedge or use a decoder that accepts one"
                ));
            }
        }
    }
    let picked = uf_decode_graph(nd + 1, boundary, &edges, &weights, detectors);
    let mut used = Vec::new();
    let mut flipped = false;
    for (ei, &on) in picked.iter().enumerate() {
        if on {
            let mi = index[ei];
            used.push(mi as u32);
            flipped ^= model.mechanisms[mi].flips_observable;
        }
    }
    Ok(DemDecode { mechanisms: used, observable_flipped: flipped })
}

/// Logical failure rate under circuit-level noise.
pub fn benchmark_circuit_level(
    code: &CssCode,
    cfg: &CircuitConfig,
    trials: u32,
    seed: u64,
) -> (u32, u32) {
    let mut fails = 0;
    for tix in 0..trials {
        let s = seed.wrapping_add(tix as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15);
        match decode_circuit_level(code, cfg, s) {
            Some(r) if r.logical_success => {}
            _ => fails += 1,
        }
    }
    (trials, fails)
}

/// Logical failure rate under phenomenological noise.
pub fn benchmark_phenomenological(
    code: &CssCode,
    cfg: &PhenomConfig,
    trials: u32,
    seed: u64,
) -> (u32, u32) {
    let mut fails = 0;
    for tix in 0..trials {
        let s = seed.wrapping_add(tix as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15);
        match decode_phenomenological(code, cfg, s) {
            Some(r) if r.logical_success => {}
            _ => fails += 1,
        }
    }
    (trials, fails)
}

/// One union-find decode of a sampled error — the matching counterpart to
/// [`decode_once`].
pub fn decode_once_uf(code: &CssCode, cfg: &DecodeConfig, seed: u64) -> DecodeResult {
    let error = sample_error(code.n, cfg.p_fx, seed);
    let syndrome_bits = syndrome_of(&code.checks, &error);
    let (correction, converged) = match union_find_correct(code, &syndrome_bits) {
        Some(c) => (c, true),
        None => (vec![false; code.n], false),
    };
    let residual: Vec<usize> = (0..code.n).filter(|&i| error[i] ^ correction[i]).collect();
    let logical_success = converged && code.is_trivial(&residual);
    DecodeResult {
        code_id: code.id.clone(),
        n: code.n,
        error: bits_to_ids(&error),
        correction: bits_to_ids(&correction),
        syndrome: bits_to_ids(&syndrome_bits),
        converged,
        logical_success,
        legs: 0,
        iters: 0,
        messages: 0,
        syndrome_bits,
        correction_bits: correction,
        check_matrix_hash: content_hash(&code.check_matrix_bytes()),
    }
}

/// Union-find logical failure rate over `trials` sampled errors.
pub fn benchmark_uf(code: &CssCode, p: f64, trials: u32, seed: u64) -> (u32, u32) {
    let cfg = DecodeConfig::from_p(p);
    let mut fails = 0;
    for tix in 0..trials {
        let s = seed.wrapping_add(tix as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15);
        if !decode_once_uf(code, &cfg, s).logical_success {
            fails += 1;
        }
    }
    (trials, fails)
}

pub fn benchmark(code: &CssCode, p: f64, trials: u32, seed: u64, max_legs: u32) -> BenchResult {
    let base = DecodeConfig::from_p(p);
    let plain = DecodeConfig { max_legs: 1, ..base };
    let relay = DecodeConfig { max_legs, ..base };
    let mut pf = 0;
    let mut rf = 0;
    let mut msg = 0u64;
    for tix in 0..trials {
        let s = seed.wrapping_add(tix as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15);
        if !decode_once(code, &plain, s).logical_success {
            pf += 1;
        }
        let r = decode_once(code, &relay, s);
        if !r.logical_success {
            rf += 1;
        }
        msg += r.messages;
    }
    BenchResult { trials, plain_logical_failures: pf, relay_logical_failures: rf, relay_total_messages: msg }
}

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

    fn key(s: u8) -> SigningKey {
        SigningKey::from_bytes(&[s; 32])
    }

    /// The `[[7,1,3]]` claim, with the distance found by exhaustive search rather
    /// than asserted — the same standard the other codes here are held to.
    #[test]
    fn steane_is_a_distance_three_colour_code() {
        let c = CssCode::color_steane();
        assert_eq!(c.n, 7);
        assert_eq!(c.k, 1, "one logical qubit");
        assert_eq!(x_distance(&c, 5), Some(3), "distance must be found, not claimed");
        assert!(c.is_self_dual(), "a colour code's two stabilizer families coincide");
        assert!(c.has_transversal_hadamard());
    }

    /// The contrast that makes the point. A surface code is *not* self-dual, so the
    /// one-round transversal Hadamard is unavailable to it — that is why it needs
    /// lattice surgery or code deformation for a logical H.
    #[test]
    fn the_surface_code_has_no_transversal_hadamard() {
        let s = CssCode::surface(3);
        assert!(!s.is_self_dual(), "surface X and Z families sit on different supports");
        assert!(!s.has_transversal_hadamard());
    }

    /// Encoding really does land in the code space: every stabilizer of the code
    /// is satisfied by the prepared state, checked by syndrome rather than by
    /// construction.
    #[cfg(feature = "quantum")]
    #[test]
    fn the_encoder_prepares_a_codeword() {
        for code in [CssCode::color_steane(), CssCode::surface(3)] {
            let sv = code.encode_zero_circuit().simulate().unwrap();
            // |0>_L is a superposition of X-stabilizer row-space vectors, so every
            // basis state with nonzero amplitude must satisfy every Z check.
            for (i, a) in sv.amps.iter().enumerate() {
                if a.norm2() == 0 {
                    continue;
                }
                for chk in &code.checks {
                    let par = chk.iter().filter(|&&q| (i >> q) & 1 == 1).count() % 2;
                    assert_eq!(par, 0, "{} left basis state {i:b} outside the code", code.id);
                }
            }
        }
    }

    /// The headline, verified end-to-end on the byte-exact simulator rather than
    /// argued: applying `H` to all seven physical qubits turns the logical `|0⟩`
    /// into the logical `|+⟩`. That is a logical Hadamard, done in one round with
    /// no ancillas and no surgery.
    #[cfg(feature = "quantum_info")]
    #[test]
    fn transversal_hadamard_implements_a_logical_hadamard() {
        use crate::quantum::ONE;
        use crate::quantum_info::{expect_pauli, Pauli};

        let code = CssCode::color_steane();
        let zl: Vec<(u8, Pauli)> = (0..7).map(|q| (q, Pauli::Z)).collect();
        let xl: Vec<(u8, Pauli)> = (0..7).map(|q| (q, Pauli::X)).collect();
        let near = |got: i64, want: i64, what: &str| {
            assert!((got - want).abs() < 1 << 13, "{what}: got {got}, want {want}");
        };

        // |0>_L is the +1 eigenstate of the logical Z and carries no logical X.
        let mut c = code.encode_zero_circuit();
        let sv = c.simulate().unwrap();
        near(expect_pauli(&sv, &zl), ONE, "<Z_L> on |0>_L");
        near(expect_pauli(&sv, &xl), 0, "<X_L> on |0>_L");

        // Transversal H: one H on every physical qubit, nothing else.
        for q in 0..7 {
            c.h(q);
        }
        let out = c.simulate().unwrap();

        // The logical operators have swapped — the definition of a logical H.
        near(expect_pauli(&out, &xl), ONE, "<X_L> after transversal H");
        near(expect_pauli(&out, &zl), 0, "<Z_L> after transversal H");

        // And it never left the code space: every stabilizer still reads +1.
        for chk in &code.checks {
            let zs: Vec<(u8, Pauli)> = chk.iter().map(|&q| (q as u8, Pauli::Z)).collect();
            let xs: Vec<(u8, Pauli)> = chk.iter().map(|&q| (q as u8, Pauli::X)).collect();
            near(expect_pauli(&out, &zs), ONE, "Z stabilizer after transversal H");
            near(expect_pauli(&out, &xs), ONE, "X stabilizer after transversal H");
        }
    }

    /// The same move on a surface code does not work, so the previous test is
    /// measuring the colour code's structure and not an artefact of the check.
    #[cfg(feature = "quantum_info")]
    #[test]
    fn transversal_hadamard_breaks_the_surface_code() {
        use crate::quantum::ONE;
        use crate::quantum_info::{expect_pauli, Pauli};

        let code = CssCode::surface(3);
        let mut c = code.encode_zero_circuit();
        for q in 0..code.n as u8 {
            c.h(q);
        }
        let out = c.simulate().unwrap();
        let violated = code.checks.iter().any(|chk| {
            let zs: Vec<(u8, Pauli)> = chk.iter().map(|&q| (q as u8, Pauli::Z)).collect();
            (expect_pauli(&out, &zs) - ONE).abs() >= 1 << 13
        });
        assert!(violated, "transversal H must knock the surface code out of its code space");
    }

    #[test]
    fn gf2_membership() {
        // rows {0,1}, {1,2} span {0,1},{1,2},{0,2}; not {0},{2}.
        let basis = Gf2Basis::from_rows(&[vec![0, 1], vec![1, 2]], 3);
        assert_eq!(basis.rank(), 2);
        assert!(basis.contains(&[0, 1]));
        assert!(basis.contains(&[1, 2]));
        assert!(basis.contains(&[0, 2]));
        assert!(basis.contains(&[])); // zero vector always in the space
        assert!(!basis.contains(&[0]));
        assert!(!basis.contains(&[2]));
    }


    /// A CSS code is only a code if the two families commute: every X stabilizer
    /// must meet every Z check in an even number of qubits. This is the test that
    /// decides whether a boundary convention is right — the first parities I tried
    /// for the surface code failed it.
    #[test]
    fn surface_code_is_a_valid_css_code() {
        for d in [3usize, 5, 7] {
            let code = CssCode::surface(d);
            assert_eq!(code.n, d * d, "d={d} qubit count");
            assert_eq!(code.k, 1, "d={d} must encode exactly one logical qubit");
            let x: Vec<&Vec<usize>> = code.stabs.iter().collect();
            for zc in &code.checks {
                for xs in &x {
                    let overlap = zc.iter().filter(|q| xs.contains(q)).count();
                    assert_eq!(overlap % 2, 0, "d={d}: X{xs:?} and Z{zc:?} anticommute");
                }
            }
            // stabilizer counts: (d^2-1)/2 of each family
            assert_eq!(code.checks.len(), (d * d - 1) / 2, "d={d} Z-check count");
            assert_eq!(x.len(), (d * d - 1) / 2, "d={d} X-stabilizer count");
            // planarity: every stabilizer is weight 4 (bulk) or weight 2 (boundary)
            for v in code.checks.iter().chain(x.iter().copied()) {
                assert!(v.len() == 4 || v.len() == 2, "d={d}: weight {} stabilizer", v.len());
            }
        }
    }

    /// The decoder must actually protect the code, and Relay-BP must earn its
    /// name: plain belief propagation degrades on quantum codes (degeneracy plus
    /// short cycles), which is the whole reason the relay legs exist.
    #[test]
    fn surface_code_decodes() {
        let d3 = CssCode::surface(3);
        let r = benchmark(&d3, 0.01, 500, 7, 4);
        let relay = r.relay_logical_failures as f64 / r.trials as f64;
        assert!(relay <= 0.05, "d=3 at p=0.01 should rarely fail: {relay}");
        assert!(
            r.relay_logical_failures <= r.plain_logical_failures,
            "relay {} must not do worse than plain BP {}",
            r.relay_logical_failures,
            r.plain_logical_failures
        );
    }

    /// Relay-BP's actual claim, on the code family it was designed for: the relay
    /// legs must beat plain belief propagation on a qLDPC (bivariate-bicycle) code.
    ///
    /// Note what is deliberately NOT asserted here: that a larger *topological*
    /// code decodes better. Measured across p = 0.001…0.02, distance-5 loses to
    /// distance-3 on both the surface code and the toric code — see
    /// `probe_surface_threshold`. That is a property of BP-family decoding on
    /// degenerate topological codes, not of the code constructions, and inventing
    /// a passing threshold for it would be asserting something untrue.
    #[test]
    fn relay_beats_plain_bp_on_qldpc() {
        let bb = CssCode::gross();
        let r = benchmark(&bb, 0.06, 400, 11, 16);
        assert!(
            r.relay_logical_failures < r.plain_logical_failures,
            "relay {} should beat plain BP {} on a qLDPC code",
            r.relay_logical_failures,
            r.plain_logical_failures
        );
    }

    #[test]
    #[ignore]
    fn probe_surface_threshold() {
        println!("\n-- relay legs on the code Relay-BP was designed for: BB [[72,12]] --");
        let bb = CssCode::gross();
        println!(" p       legs=1(plain)  legs=4    legs=16");
        for &pp in &[0.02f64, 0.04, 0.06] {
            let f = |legs: u32| {
                let r = benchmark(&bb, pp, 400, 11, legs);
                r.relay_logical_failures as f64 / r.trials as f64
            };
            let plain = {
                let r = benchmark(&bb, pp, 400, 11, 4);
                r.plain_logical_failures as f64 / r.trials as f64
            };
            println!(" {pp:<8} {plain:<13.4} {:<9.4} {:<9.4}", f(4), f(16));
        }
        println!("\n-- the same sweep on topological codes --");
        println!(" p       surf d=3  surf d=5  toric L=3 toric L=5");
        for &pp in &[0.001f64, 0.005, 0.01] {
            let sf = |d: usize| {
                let r = benchmark(&CssCode::surface(d), pp, 1500, 11, 16);
                r.relay_logical_failures as f64 / r.trials as f64
            };
            let tf = |l: usize| {
                let r = benchmark(&CssCode::toric(l), pp, 1500, 11, 16);
                r.relay_logical_failures as f64 / r.trials as f64
            };
            println!(" {pp:<8} {:<9.4} {:<9.4} {:<9.4} {:<9.4}", sf(3), sf(5), tf(3), tf(5));
        }
    }



    /// The property BP could not deliver: with a matching decoder, distance counts.
    /// Under Relay-BP a distance-5 surface code decoded WORSE than distance-3;
    /// union-find inverts that decisively.
    #[test]
    fn union_find_makes_distance_count() {
        let p = 0.01;
        let rate = |d: usize| {
            let (t, f) = benchmark_uf(&CssCode::surface(d), p, 1500, 11);
            f as f64 / t as f64
        };
        let (r3, r5) = (rate(3), rate(5));
        assert!(r5 < r3, "d=5 ({r5}) must beat d=3 ({r3}) below threshold");
        assert!(r5 <= 0.005, "d=5 at p={p} should almost never fail: {r5}");
    }

    /// And it must beat belief propagation on the codes BP is bad at.
    #[test]
    fn union_find_beats_bp_on_topological_codes() {
        let p = 0.02;
        for d in [5usize, 7] {
            let code = CssCode::surface(d);
            let bp = {
                let r = benchmark(&code, p, 1000, 11, 16);
                r.relay_logical_failures as f64 / r.trials as f64
            };
            let (t, f) = benchmark_uf(&code, p, 1000, 11);
            let uf = f as f64 / t as f64;
            assert!(uf < bp, "d={d}: union-find {uf} should beat Relay-BP {bp}");
        }
    }

    /// A qubit in three checks is not a matching graph. Refuse, do not guess —
    /// Relay-BP is the decoder for those families.
    #[test]
    fn union_find_refuses_non_matching_codes() {
        assert!(
            union_find_correct(&CssCode::gross(), &vec![false; CssCode::gross().checks.len()]).is_none(),
            "a bivariate-bicycle qubit sits in three checks; UF must decline"
        );
        assert!(union_find_correct(&CssCode::surface(5), &vec![false; 12]).is_some());
    }

    /// Ground truth: every single-qubit error must be corrected exactly.
    #[test]
    fn union_find_corrects_every_single_qubit_error() {
        for d in [3usize, 5] {
            let code = CssCode::surface(d);
            for q in 0..code.n {
                let mut err = vec![false; code.n];
                err[q] = true;
                let syn = syndrome_of(&code.checks, &err);
                let corr = union_find_correct(&code, &syn).expect("matching graph");
                let residual: Vec<usize> = (0..code.n).filter(|&i| err[i] ^ corr[i]).collect();
                assert!(
                    code.is_trivial(&residual),
                    "d={d}: single error on qubit {q} left a logical residual {residual:?}"
                );
            }
        }
    }



    /// Measurement errors alone must never become a logical error. With no data
    /// noise the true error is empty, so any space edge the decoder picks IS the
    /// residual — this asserts it explains lying measurements with time edges
    /// instead of inventing data corrections.
    #[test]
    fn measurement_noise_alone_never_causes_a_logical_error() {
        for d in [3usize, 5] {
            let code = CssCode::surface(d);
            let cfg = PhenomConfig::from_p(d as u32, 0.0, 0.05);
            for t in 0..300u64 {
                let r = decode_phenomenological(&code, &cfg, t.wrapping_mul(0x9E37_79B9))
                    .expect("matching graph");
                assert!(
                    r.logical_success,
                    "d={d} seed={t}: pure measurement noise produced a logical failure \
                     (correction {:?} on an empty data error)",
                    r.correction
                );
            }
        }
    }

    /// Below threshold, distance must still buy protection when the syndrome
    /// itself is unreliable — the property that matters for a real device.
    #[test]
    fn phenomenological_distance_helps() {
        let p = 0.005;
        let rate = |d: usize| {
            let cfg = PhenomConfig::from_p(d as u32, p, p);
            let (t, f) = benchmark_phenomenological(&CssCode::surface(d), &cfg, 1200, 11);
            f as f64 / t as f64
        };
        let (r3, r5) = (rate(3), rate(5));
        assert!(r5 < r3, "d=5 ({r5}) must beat d=3 ({r3}) under noisy measurement");
    }

    /// With perfect measurement the space-time decoder must reduce to the
    /// single-shot one.
    #[test]
    fn perfect_measurement_reduces_to_single_shot() {
        for d in [3usize, 5] {
            let code = CssCode::surface(d);
            let cfg = PhenomConfig::from_p(1, 0.02, 0.0);
            let (t1, f1) = benchmark_phenomenological(&code, &cfg, 1200, 11);
            let (t2, f2) = benchmark_uf(&code, 0.02, 1200, 11);
            let (a, b) = (f1 as f64 / t1 as f64, f2 as f64 / t2 as f64);
            assert!((a - b).abs() < 0.02, "d={d}: spacetime {a} vs single-shot {b}");
        }
    }

    #[test]
    fn phenomenological_refuses_non_matching_codes() {
        let cfg = PhenomConfig::from_p(3, 0.01, 0.01);
        assert!(decode_phenomenological(&CssCode::gross(), &cfg, 1).is_none());
    }


    /// The `d` in `[[d², 1, d]]` — verified, not asserted by naming. A code whose
    /// distance is smaller than advertised silently corrects less than it claims.
    #[test]
    fn surface_code_distance_is_d() {
        for d in [3usize, 5] {
            let code = CssCode::surface(d);
            assert_eq!(
                x_distance(&code, d + 1),
                Some(d),
                "surface d={d} should have X-distance exactly {d}"
            );
        }
    }

    /// The toric code's distance is L, and its k is 2 — a different point in the
    /// same family, so this also guards the distance search itself.
    #[test]
    fn toric_code_distance_is_l() {
        let code = CssCode::toric(3);
        assert_eq!(code.k, 2, "toric encodes two logical qubits");
        assert_eq!(x_distance(&code, 4), Some(3));
    }

    /// Below the distance, correction must be guaranteed — not merely likely.
    /// Every error of weight < d, exhaustively, must decode to a trivial residual.
    #[test]
    fn union_find_corrects_every_error_below_the_distance() {
        for d in [3usize, 5] {
            let code = CssCode::surface(d);
            let w = (d - 1) / 2; // guaranteed-correctable weight
            let mut combo: Vec<usize> = (0..w).collect();
            let mut checked = 0u32;
            loop {
                let mut err = vec![false; code.n];
                for &q in &combo {
                    err[q] = true;
                }
                let syn = syndrome_of(&code.checks, &err);
                let corr = union_find_correct(&code, &syn).expect("matching graph");
                let residual: Vec<usize> =
                    (0..code.n).filter(|&i| err[i] ^ corr[i]).collect();
                assert!(
                    code.is_trivial(&residual),
                    "d={d}: weight-{w} error {combo:?} decoded to a logical residual {residual:?}"
                );
                checked += 1;
                if !next_combination(&mut combo, code.n) {
                    break;
                }
            }
            assert!(checked > 0, "d={d}: nothing checked");
        }
    }



    /// With no ancilla faults, circuit-level noise IS phenomenological noise. Not
    /// "statistically similar" — the same stream, the same failures, bit for bit.
    /// A model that only approximately reduces to its own limiting case is hiding
    /// a difference somewhere.
    #[test]
    fn hooks_off_is_exactly_phenomenological() {
        for &pp in &[0.005f64, 0.02] {
            for d in [3usize, 5] {
                let code = CssCode::surface(d);
                let ph = PhenomConfig::from_p(d as u32, pp, pp);
                let cl = CircuitConfig::from_p(d as u32, pp, pp, 0.0);
                let (_, a) = benchmark_phenomenological(&code, &ph, 600, 11);
                let (_, b) = benchmark_circuit_level(&code, &cl, 600, 11);
                assert_eq!(a, b, "p={pp} d={d}: hook-free circuit level must equal phenomenological");
            }
        }
    }

    /// Hook errors must cost something — that is the entire point of modelling
    /// them. One ancilla fault becomes several correlated data errors the matching
    /// graph has no edge for.
    #[test]
    fn hook_errors_cost_accuracy() {
        let (p, d) = (0.005, 5usize);
        let code = CssCode::surface(d);
        let off = CircuitConfig::from_p(d as u32, p, p, 0.0);
        let on = CircuitConfig::from_p(d as u32, p, p, p);
        let (_, f_off) = benchmark_circuit_level(&code, &off, 1500, 11);
        let (_, f_on) = benchmark_circuit_level(&code, &on, 1500, 11);
        assert!(f_on > f_off, "hooks on ({f_on}) should fail more than hooks off ({f_off})");
    }

    /// A hook lands on the CNOTs a circuit has not reached yet, so it can never be
    /// heavier than the stabilizer it came from — and never empty.
    #[test]
    fn hook_weight_is_bounded_by_the_stabilizer() {
        let code = CssCode::surface(5);
        let maxw = code.stabs.iter().map(|s| s.len()).max().unwrap();
        let cfg = CircuitConfig::from_p(5, 0.0, 0.0, 0.2);
        let mut seen = 0;
        for t in 0..200u64 {
            let r = decode_circuit_level(&code, &cfg, t.wrapping_mul(0x9E37_79B9)).unwrap();
            for &w in &r.hook_events {
                assert!(w >= 1 && w as usize <= maxw, "hook weight {w} outside 1..={maxw}");
                seen += 1;
            }
        }
        assert!(seen > 0, "the probe never fired a hook");
    }

    /// A noiseless circuit must never fail. Trivial, and exactly the kind of thing
    /// that quietly breaks when a model grows a new fault channel.
    #[test]
    fn noiseless_circuit_never_fails() {
        let code = CssCode::surface(5);
        let cfg = CircuitConfig::from_p(5, 0.0, 0.0, 0.0);
        for t in 0..200u64 {
            let r = decode_circuit_level(&code, &cfg, t).unwrap();
            assert!(r.logical_success && r.data_error.is_empty() && r.correction.is_empty());
        }
    }

    #[test]
    fn circuit_level_refuses_non_matching_codes() {
        let cfg = CircuitConfig::from_p(3, 0.01, 0.01, 0.01);
        assert!(decode_circuit_level(&CssCode::gross(), &cfg, 1).is_none());
    }


    /// The logical Z must actually be the ruler it claims: it commutes with every
    /// X stabilizer, is no product of Z checks, and — the property that makes it
    /// usable — meets a minimum-weight logical X in an ODD number of qubits. Even
    /// overlap would mean it cannot detect that logical at all.
    #[test]
    fn logical_z_is_a_valid_observable() {
        for d in [3usize, 5] {
            let code = CssCode::surface(d);
            let lz = logical_z_support(&code, d + 1).expect("a logical Z exists");
            assert_eq!(lz.len(), d, "d={d}: minimum-weight logical Z should have weight {d}");
            // commutes with all X stabilizers
            let mut bits = vec![false; code.n];
            for &q in &lz { bits[q] = true; }
            assert!(syndrome_of(&code.stabs, &bits).iter().all(|b| !b));
            // and anticommutes with a minimum-weight logical X
            let dx = x_distance(&code, d + 1).expect("a logical X exists");
            assert_eq!(dx, d);
            // find that logical X and check odd overlap
            let mut combo: Vec<usize> = (0..d).collect();
            let mut found = false;
            loop {
                let mut e = vec![false; code.n];
                for &q in &combo { e[q] = true; }
                if syndrome_of(&code.checks, &e).iter().all(|b| !b) && !code.is_trivial(&combo) {
                    let overlap = combo.iter().filter(|q| lz.contains(q)).count();
                    assert_eq!(overlap % 2, 1, "d={d}: logical Z {lz:?} misses logical X {combo:?}");
                    found = true;
                    break;
                }
                if !next_combination(&mut combo, code.n) { break; }
            }
            assert!(found, "d={d}: no minimum-weight logical X found");
        }
    }


    /// The defining property of a logical error, stated purely in DEM terms: take
    /// the mechanisms of a minimum-weight logical operator, and their combined
    /// effect must light NO detectors while flipping the observable an ODD number
    /// of times. Undetectable, yet it changes the answer.
    ///
    /// A DEM that gets this wrong is worse than useless — an external decoder
    /// would silently mis-score every shot — and it checks the detector sets and
    /// the observable flags against each other rather than against my own formula.
    #[test]
    fn dem_encodes_the_codes_logical_structure() {
        for d in [3usize, 5] {
            let code = CssCode::surface(d);
            let cfg = CircuitConfig::from_p(2, 0.001, 0.001, 0.0);
            let mechs = dem_mechanisms(&code, &cfg, d + 1).expect("matching graph");

            // a minimum-weight logical X
            let mut combo: Vec<usize> = (0..d).collect();
            let logical_x = loop {
                let mut e = vec![false; code.n];
                for &q in &combo { e[q] = true; }
                if syndrome_of(&code.checks, &e).iter().all(|b| !b) && !code.is_trivial(&combo) {
                    break combo.clone();
                }
                assert!(next_combination(&mut combo, code.n), "d={d}: no logical X of weight {d}");
            };

            // its data mechanisms at round 0 are the first `n` entries
            let mut detectors: Vec<u32> = Vec::new();
            let mut obs = false;
            for &q in &logical_x {
                let mech = &mechs[q]; // round 0, qubit q
                obs ^= mech.flips_observable;
                for &dt in &mech.detectors {
                    if let Some(pos) = detectors.iter().position(|&x| x == dt) {
                        detectors.remove(pos);
                    } else {
                        detectors.push(dt);
                    }
                }
            }
            assert!(detectors.is_empty(), "d={d}: logical X {logical_x:?} lit detectors {detectors:?}");
            assert!(obs, "d={d}: logical X {logical_x:?} did not flip the observable");
        }
    }

    /// A stabilizer is the converse case: also undetectable, but it must NOT flip
    /// the observable. Getting only the first half right would still be wrong.
    #[test]
    fn dem_says_stabilizers_are_harmless() {
        let d = 5usize;
        let code = CssCode::surface(d);
        let cfg = CircuitConfig::from_p(2, 0.001, 0.001, 0.0);
        let mechs = dem_mechanisms(&code, &cfg, d + 1).unwrap();
        for stab in &code.stabs {
            let mut detectors: Vec<u32> = Vec::new();
            let mut obs = false;
            for &q in stab {
                let mech = &mechs[q];
                obs ^= mech.flips_observable;
                for &dt in &mech.detectors {
                    if let Some(pos) = detectors.iter().position(|&x| x == dt) {
                        detectors.remove(pos);
                    } else {
                        detectors.push(dt);
                    }
                }
            }
            assert!(detectors.is_empty(), "stabilizer {stab:?} lit detectors {detectors:?}");
            assert!(!obs, "stabilizer {stab:?} flipped the observable");
        }
    }

    /// The invariant that actually matters for a DEM: no mechanism may be
    /// undetectable while flipping the observable. Such a fault is a logical error
    /// no decoder could ever see coming, and a DEM containing one would make every
    /// external decoder silently mis-score.
    ///
    /// Note what is NOT asserted: that hooks appear as hyperedges. Whether they do
    /// depends on the CNOT order — with the index order used here a partial hook on
    /// the rotated surface code lights at most two detectors, so this model is an
    /// ordinary graph. Asserting hyperedges would have been asserting an assumption.
    #[test]
    fn no_dem_mechanism_is_an_invisible_logical() {
        for d in [3usize, 5] {
            let code = CssCode::surface(d);
            for hook in [0.0, 0.002] {
                let cfg = CircuitConfig::from_p(2, 0.001, 0.001, hook);
                for m in dem_mechanisms(&code, &cfg, d + 1).unwrap() {
                    assert!(
                        !(m.detectors.is_empty() && m.flips_observable),
                        "d={d} hook={hook}: undetectable mechanism flips the observable"
                    );
                    assert!(
                        !m.detectors.is_empty(),
                        "d={d} hook={hook}: a mechanism with no detectors should have been dropped"
                    );
                }
            }
        }
    }

    #[test]
    fn dem_text_is_wellformed() {
        let code = CssCode::surface(3);
        let cfg = CircuitConfig::from_p(2, 0.002, 0.002, 0.001);
        let dem = to_dem(&code, &cfg, 4).unwrap();
        let n_det = 3 * code.checks.len(); // rounds+1 layers
        assert_eq!(dem.lines().filter(|l| l.starts_with("detector(")).count(), n_det);
        let errs: Vec<&str> = dem.lines().filter(|l| l.starts_with("error(")).collect();
        assert!(!errs.is_empty());
        for l in &errs {
            let p: f64 = l[6..l.find(')').unwrap()].parse().expect("probability parses");
            assert!(p > 0.0 && p < 1.0, "probability {p} out of range in {l}");
            assert!(l.contains(" D"), "mechanism with no detectors: {l}");
            for tok in l.split_whitespace().skip(1) {
                let idx: u32 = tok.trim_start_matches(['D', 'L']).parse().expect("target index");
                if tok.starts_with('D') {
                    assert!((idx as usize) < n_det, "detector {idx} out of range");
                } else {
                    assert_eq!(idx, 0, "only one observable is defined");
                }
            }
        }
        assert!(dem.contains("logical") || dem.contains("L0"), "no observable referenced");
    }




    /// A format that only writes is half a format.
    #[test]
    fn dem_round_trips() {
        for d in [3usize, 5] {
            let code = CssCode::surface(d);
            let cfg = CircuitConfig::from_p(2, 0.002, 0.003, 0.0);
            let text = to_dem(&code, &cfg, d + 1).unwrap();
            let parsed = from_dem(&text).expect("our own output must parse");
            let native = dem_mechanisms(&code, &cfg, d + 1).unwrap();
            assert_eq!(parsed.mechanisms.len(), native.len(), "d={d}: mechanism count");
            assert_eq!(parsed.mechanisms, native, "d={d}: mechanisms changed through the text");
            assert_eq!(parsed.num_detectors, (cfg.rounds as usize + 1) * code.checks.len());
        }
    }

    /// The test that decides whether the export is worth anything: decode the SAME
    /// shots through the DEM and natively, and the DEM must do as well. If it does
    /// not, the exported model is missing information the decoder needs, and every
    /// external consumer would inherit that loss.
    #[test]
    fn decoding_through_the_dem_matches_the_native_decoder() {
        for d in [3usize, 5] {
            let code = CssCode::surface(d);
            let cfg = CircuitConfig::from_p(d as u32, 0.01, 0.01, 0.0);
            let phen = PhenomConfig { rounds: cfg.rounds, p_data_fx: cfg.p_data_fx, q_meas_fx: cfg.q_meas_fx };
            let model = from_dem(&to_dem(&code, &cfg, d + 1).unwrap()).unwrap();
            let logical = logical_z_support(&code, d + 1).unwrap();
            let m = code.checks.len();
            let n_det = (cfg.rounds as usize + 1) * m;

            let (mut native_fail, mut dem_fail) = (0u32, 0u32);
            let trials = 600u32;
            for tix in 0..trials {
                let seed = (tix as u64).wrapping_add(7).wrapping_mul(0x9E37_79B9_7F4A_7C15);
                let r = decode_phenomenological(&code, &phen, seed).unwrap();
                if !r.logical_success {
                    native_fail += 1;
                }
                let mut bits = vec![false; n_det];
                for &dd in &r.detectors {
                    bits[dd as usize] = true;
                }
                let truth = r.data_error.iter().filter(|q| logical.contains(&(**q as usize))).count() % 2 == 1;
                let pred = decode_dem(&model, &bits).unwrap().observable_flipped;
                if pred != truth {
                    dem_fail += 1;
                }
            }
            let (a, b) = (native_fail as f64 / trials as f64, dem_fail as f64 / trials as f64);
            assert!(
                b <= a + 0.02,
                "d={d}: decoding through the DEM ({b}) lost accuracy vs native ({a})"
            );
        }
    }

    /// Refuse rather than quietly decode a different experiment.
    #[test]
    fn dem_ingest_refuses_what_it_cannot_honour() {
        // a hyperedge cannot be matched
        let hyper = "error(0.001) D0 D1 D2\n";
        let model = from_dem(hyper).unwrap();
        let e = decode_dem(&model, &[false; 3]).unwrap_err();
        assert!(e.contains("three or more") || e.contains("at most two"), "{e}");

        // an instruction we do not implement must not be skipped
        assert!(from_dem("repeat 5 {\n").unwrap_err().contains("unsupported"));
        assert!(from_dem("error(2.0) D0\n").unwrap_err().contains("outside"));
        assert!(from_dem("error(0.1) X3\n").unwrap_err().contains("unexpected target"));
    }

    /// A foreign DEM — not our text — must decode. Two detectors, one boundary
    /// edge, one observable-flipping mechanism.
    #[test]
    fn ingests_a_hand_written_dem() {
        let src = "\
# hand written
detector(0, 0) D0
detector(1, 0) D1
error(0.01) D0 D1
error(0.02) D0 L0
error(0.02) D1 L0
";
        let model = from_dem(src).expect("parses");
        assert_eq!(model.num_detectors, 2);
        assert_eq!(model.mechanisms.len(), 3);
        // lighting both detectors is explained by the single cheap D0 D1 edge,
        // which does not touch the observable
        let r = decode_dem(&model, &[true, true]).unwrap();
        assert!(!r.observable_flipped, "the D0-D1 mechanism should explain it");
        // lighting one detector must go to the boundary, flipping the observable
        let r = decode_dem(&model, &[true, false]).unwrap();
        assert!(r.observable_flipped);
    }



    #[test]
    #[ignore]
    fn probe_perf() {
        use std::time::Instant;
        println!("\n  decoder throughput (surface code, p=0.01)");
        for d in [3usize, 5, 7, 9, 11] {
            let code = CssCode::surface(d);
            let t0 = Instant::now();
            let (trials, _) = benchmark_uf(&code, 0.01, 500, 11);
            let el = t0.elapsed();
            println!("   d={d:<3} n={:<4} checks={:<4} {:>8.1} shots/s   ({:?}/shot)",
                code.n, code.checks.len(),
                trials as f64 / el.as_secs_f64(), el / trials);
        }
        println!("\n  space-time decoder (d rounds, p=q=0.01)");
        for d in [3usize, 5, 7] {
            let code = CssCode::surface(d);
            let cfg = PhenomConfig::from_p(d as u32, 0.01, 0.01);
            let t0 = Instant::now();
            let (trials, _) = benchmark_phenomenological(&code, &cfg, 200, 11);
            let el = t0.elapsed();
            println!("   d={d:<3} detectors={:<5} {:>8.1} shots/s", (d + 1) * code.checks.len(),
                trials as f64 / el.as_secs_f64());
        }
        println!("\n  code construction + distance verification");
        for d in [3usize, 5] {
            let t0 = Instant::now(); let c = CssCode::surface(d); let build = t0.elapsed();
            let t1 = Instant::now(); let dd = x_distance(&c, d + 1); let dist = t1.elapsed();
            println!("   d={d}: build {build:?}, distance {dd:?} in {dist:?}");
        }
    }

    #[test]
    #[ignore]
    fn probe_dem_equivalence() {
        println!("\n  same shots, decoded natively vs through the exported DEM");
        println!("  d   p       native   via DEM");
        for d in [3usize, 5] {
            for &pp in &[0.005f64, 0.01, 0.02] {
                let code = CssCode::surface(d);
                let cfg = CircuitConfig::from_p(d as u32, pp, pp, 0.0);
                let phen = PhenomConfig { rounds: cfg.rounds, p_data_fx: cfg.p_data_fx, q_meas_fx: cfg.q_meas_fx };
                let model = from_dem(&to_dem(&code, &cfg, d + 1).unwrap()).unwrap();
                let logical = logical_z_support(&code, d + 1).unwrap();
                let n_det = (cfg.rounds as usize + 1) * code.checks.len();
                let (mut a, mut b) = (0u32, 0u32);
                let trials = 1500u32;
                for tix in 0..trials {
                    let seed = (tix as u64).wrapping_add(7).wrapping_mul(0x9E37_79B9_7F4A_7C15);
                    let r = decode_phenomenological(&code, &phen, seed).unwrap();
                    if !r.logical_success { a += 1; }
                    let mut bits = vec![false; n_det];
                    for &dd in &r.detectors { bits[dd as usize] = true; }
                    let truth = r.data_error.iter().filter(|q| logical.contains(&(**q as usize))).count() % 2 == 1;
                    if decode_dem(&model, &bits).unwrap().observable_flipped != truth { b += 1; }
                }
                println!("  {d}   {pp:<7} {:<8.4} {:<8.4}", a as f64/trials as f64, b as f64/trials as f64);
            }
        }
    }

    #[test]
    #[ignore]
    fn probe_print_dem() {
        let code = CssCode::surface(3);
        let cfg = CircuitConfig::from_p(1, 0.002, 0.002, 0.001);
        let dem = to_dem(&code, &cfg, 4).unwrap();
        for l in dem.lines().take(22) { println!("{l}"); }
        println!("... ({} lines total)", dem.lines().count());
    }

    #[test]
    #[ignore]
    fn probe_dem_shapes() {
        let code = CssCode::surface(5);
        let cfg = CircuitConfig::from_p(2, 0.001, 0.001, 0.001);
        let mechs = dem_mechanisms(&code, &cfg, 6).unwrap();
        let mut hist = std::collections::BTreeMap::new();
        for m in &mechs { *hist.entry(m.detectors.len()).or_insert(0u32) += 1; }
        println!("\n  detector-set size histogram: {hist:?}");
        let empty = mechs.iter().filter(|m| m.detectors.is_empty()).count();
        println!("  mechanisms with NO detectors: {empty} (of {})", mechs.len());
        let empty_obs = mechs.iter().filter(|m| m.detectors.is_empty() && m.flips_observable).count();
        println!("  empty AND flips observable (would be a silent logical): {empty_obs}");
        println!("  X stabilizer weights: {:?}", {
            let mut w: Vec<usize> = code.stabs.iter().map(|s| s.len()).collect(); w.sort_unstable(); w.dedup(); w });
        // what do partial hooks look like?
        for stab in code.stabs.iter().take(3) {
            print!("  stab {stab:?}: ");
            for j in 0..stab.len() {
                let qs = &stab[j..];
                let mut set: Vec<usize> = Vec::new();
                for &q in qs {
                    for (ci, chk) in code.checks.iter().enumerate() {
                        if chk.contains(&q) {
                            if let Some(pos) = set.iter().position(|&x| x == ci) { set.remove(pos); } else { set.push(ci); }
                        }
                    }
                }
                print!("j={j}->{} ", set.len());
            }
            println!();
        }
    }

    #[test]
    #[ignore]
    fn probe_circuit_level() {
        println!("\n  hook=0 must reproduce phenomenological exactly");
        for &pp in &[0.005f64, 0.02] {
            for d in [3usize, 5] {
                let ph = PhenomConfig::from_p(d as u32, pp, pp);
                let cl = CircuitConfig::from_p(d as u32, pp, pp, 0.0);
                let (_, a) = benchmark_phenomenological(&CssCode::surface(d), &ph, 800, 11);
                let (_, b) = benchmark_circuit_level(&CssCode::surface(d), &cl, 800, 11);
                println!("   p={pp} d={d}: phenom {a}  circuit(hook=0) {b}  {}",
                    if a == b { "IDENTICAL" } else { "DIFFER" });
            }
        }
        println!("\n  hook errors lower the threshold (surface code, q=p, hook=p)");
        println!("  p        d=3       d=5       d=7");
        for &pp in &[0.001f64, 0.002, 0.005, 0.01] {
            let r = |d: usize| {
                let cfg = CircuitConfig::from_p(d as u32, pp, pp, pp);
                let (t, f) = benchmark_circuit_level(&CssCode::surface(d), &cfg, 1200, 11);
                f as f64 / t as f64
            };
            println!("  {pp:<8} {:<9.4} {:<9.4} {:<9.4}", r(3), r(5), r(7));
        }
        println!("\n  same p, hooks off vs on (d=5)");
        for &pp in &[0.002f64, 0.005] {
            let off = { let c=CircuitConfig::from_p(5,pp,pp,0.0); let (t,f)=benchmark_circuit_level(&CssCode::surface(5),&c,1200,11); f as f64/t as f64 };
            let on  = { let c=CircuitConfig::from_p(5,pp,pp,pp);  let (t,f)=benchmark_circuit_level(&CssCode::surface(5),&c,1200,11); f as f64/t as f64 };
            println!("   p={pp:<7} hooks off {:<8.4} hooks on {:<8.4}", off, on);
        }
    }

    #[test]
    #[ignore]
    fn probe_phenomenological() {
        println!("\n  surface code, phenomenological noise (data p, measurement q=p), d rounds");
        println!("  p        d=3       d=5       d=7");
        for &pp in &[0.001f64, 0.003, 0.005, 0.01, 0.02] {
            let r = |d: usize| {
                let cfg = PhenomConfig::from_p(d as u32, pp, pp);
                let (t, f) = benchmark_phenomenological(&CssCode::surface(d), &cfg, 1200, 11);
                f as f64 / t as f64
            };
            println!("  {pp:<8} {:<9.4} {:<9.4} {:<9.4}", r(3), r(5), r(7));
        }
        println!("\n  sanity: q=0 (perfect measurement) should track the single-shot decoder");
        for &pp in &[0.005f64, 0.02] {
            let st = |d: usize| {
                let cfg = PhenomConfig::from_p(1, pp, 0.0);
                let (t, f) = benchmark_phenomenological(&CssCode::surface(d), &cfg, 1200, 11);
                f as f64 / t as f64
            };
            let ss = |d: usize| { let (t,f)=benchmark_uf(&CssCode::surface(d), pp, 1200, 11); f as f64/t as f64 };
            println!("  p={pp:<7} d=3 spacetime {:<8.4} single-shot {:<8.4} | d=5 spacetime {:<8.4} single-shot {:<8.4}",
                st(3), ss(3), st(5), ss(5));
        }
    }

    #[test]
    #[ignore]
    fn probe_union_find_vs_bp() {
        println!("\n  surface code — logical failure rate");
        println!("  p        d=3 BP    d=3 UF    d=5 BP    d=5 UF    d=7 BP    d=7 UF");
        for &pp in &[0.001f64, 0.005, 0.01, 0.02, 0.05] {
            let bp = |d: usize| {
                let r = benchmark(&CssCode::surface(d), pp, 2000, 11, 16);
                r.relay_logical_failures as f64 / r.trials as f64
            };
            let uf = |d: usize| {
                let (t, f) = benchmark_uf(&CssCode::surface(d), pp, 2000, 11);
                f as f64 / t as f64
            };
            println!("  {pp:<8} {:<9.4} {:<9.4} {:<9.4} {:<9.4} {:<9.4} {:<9.4}",
                bp(3), uf(3), bp(5), uf(5), bp(7), uf(7));
        }
        println!("\n  toric code — UF");
        for &pp in &[0.005f64, 0.02] {
            let uf = |l: usize| { let (t,f)=benchmark_uf(&CssCode::toric(l), pp, 2000, 11); f as f64/t as f64 };
            println!("  p={pp:<7} L=3 {:<8.4} L=5 {:<8.4}", uf(3), uf(5));
        }
    }

    #[test]
    fn toric_code_shape() {
        let c = CssCode::toric(5);
        assert_eq!(c.n, 50);
        assert_eq!(c.checks.len(), 25);
        assert_eq!(c.k, 2, "toric code encodes 2 logical qubits");
        // every check has weight 4
        assert!(c.checks.iter().all(|r| r.len() == 4));
    }

    #[test]
    fn gross_code_is_72_12() {
        let c = CssCode::gross();
        assert_eq!(c.n, 72);
        assert_eq!(c.k, 12, "the gross-style BB code encodes 12 logical qubits");
        // qLDPC: bounded check weight (6 here), independent of n
        assert!(c.checks.iter().all(|r| r.len() == 6));
    }

    #[test]
    fn css_orthogonality_gross() {
        // H_X · H_Zᵀ = 0: every Z-check commutes with every X-stabilizer, i.e.
        // every decode check overlaps every stabilizer in an even number of qubits.
        let c = CssCode::gross();
        // stab_basis holds H_X's row space; reconstruct raw H_X rows to test.
        let hx = CssCode::bivariate_bicycle(6, 6, &[(3, 0), (0, 1), (0, 2)], &[(0, 3), (1, 0), (2, 0)]);
        for zc in &c.checks {
            let zset: std::collections::HashSet<usize> = zc.iter().copied().collect();
            for xs in raw_hx(&hx) {
                let overlap = xs.iter().filter(|q| zset.contains(q)).count();
                assert_eq!(overlap % 2, 0, "CSS orthogonality violated");
            }
        }
    }

    // helper: recover raw H_X rows (stab side) for the orthogonality test
    fn raw_hx(_c: &CssCode) -> Vec<Vec<usize>> {
        // rebuild H_X = [A|B] directly
        let (l, m) = (6usize, 6usize);
        let a = [(3usize, 0usize), (0, 1), (0, 2)];
        let b = [(0usize, 3usize), (1, 0), (2, 0)];
        let lm = l * m;
        let cell = |i: usize, j: usize| (i % l) * m + (j % m);
        let build = |monos: &[(usize, usize)]| -> Vec<Vec<usize>> {
            let mut rows = vec![Vec::new(); lm];
            for i in 0..l {
                for j in 0..m {
                    let col = cell(i, j);
                    for &(px, py) in monos {
                        rows[cell(i + px, j + py)].push(col);
                    }
                }
            }
            rows
        };
        let ar = build(&a);
        let br = build(&b);
        (0..lm)
            .map(|r| {
                let mut row = ar[r].clone();
                row.extend(br[r].iter().map(|&c| c + lm));
                row.sort_unstable();
                row
            })
            .collect()
    }

    #[test]
    fn decode_recovers_low_weight_errors() {
        // At a low error rate, decoding should almost always succeed logically.
        let c = CssCode::toric(5);
        let cfg = DecodeConfig::from_p(0.02);
        let mut fails = 0;
        for s in 0..200u64 {
            if !decode_once(&c, &cfg, s.wrapping_mul(0x9E37_79B9)).logical_success {
                fails += 1;
            }
        }
        assert!(fails < 20, "too many logical failures at p=0.02: {fails}/200");
    }

    #[test]
    fn decode_is_deterministic() {
        let c = CssCode::gross();
        let cfg = DecodeConfig::from_p(0.03);
        let a = decode_once(&c, &cfg, 12345);
        let b = decode_once(&c, &cfg, 12345);
        assert_eq!(a.correction, b.correction);
        assert_eq!(a.messages, b.messages);
        assert_eq!(a.legs, b.legs);
    }

    #[test]
    fn relay_beats_plain_bp() {
        // The whole point of Relay-BP: on a quantum code, relay legs recover
        // decodes that plain BP leaves failed. Relay failures ≤ plain failures.
        let c = CssCode::gross();
        let bench = benchmark(&c, 0.04, 150, 0xBEEF, 12);
        assert!(
            bench.relay_logical_failures <= bench.plain_logical_failures,
            "relay {} must not exceed plain {}",
            bench.relay_logical_failures, bench.plain_logical_failures
        );
    }

    #[test]
    fn decode_seals_verifying_receipt() {
        let c = CssCode::gross();
        let cfg = DecodeConfig::from_p(0.03);
        let r = decode_once(&c, &cfg, 7);
        let rec = r.seal(&key(1), "did:key:lab", 500_000, GrantRef::unbounded("quantum.decode"));
        assert!(rec.verify());
        assert!(rec.correction_matches(&r.correction_bytes()));
        assert_eq!(rec.code_id, "bb:[[72,12]]");
        assert_eq!(rec.work_messages, r.messages);
    }
}