wai-quantum 0.3.38

//! Noisy Clifford circuits at QEC scale — `wai.quantum.frame`.
//!
//! An error-correction experiment is a Clifford circuit with noise sprinkled
//! through it, a set of **detectors** (parities of measurements that are
//! deterministic in the absence of noise) and **observables** (the logical
//! readout). Two questions are asked of it, millions of times:
//!
//! - **What does a run look like?** Which detectors fired, and did the
//!   observable flip? [`sample`] answers it by **Pauli-frame simulation**. Only
//!   the *difference* a noise realisation makes is tracked, as an X bit and a Z
//!   bit per qubit, and 64 shots are packed into each machine word. Gates are
//!   word-wide XORs. Noise is drawn by skipping geometrically between hits, so a
//!   rate-`p` channel costs `p` work per shot rather than one random draw.
//! - **How can it fail?** [`error_model`] answers it by propagating
//!   *sensitivities* backwards through the circuit. A sensitivity is the set of
//!   detectors and observables that an X or a Z error at that point would flip.
//!   Every noise channel then becomes a list of independent mechanisms, each with
//!   the detectors and observables it flips. Identical ones are merged. The
//!   result is the detector error model a decoder consumes, written in the
//!   ecosystem's text syntax by [`ErrorModel::to_text`].
//!
//! # The contract
//!
//! - **Sampling is reproducible bit for bit.** A seeded `splitmix64` stream is
//!   consumed in a fixed order, and the geometric skips take their logarithm
//!   in-crate. A `(circuit, shots, seed)` triple yields the same detection
//!   events on every machine and target.
//! - **The error model is exact under the stated independence.** A depolarizing
//!   channel's disjoint outcomes are converted to independent per-Pauli
//!   mechanisms with exactly the same output distribution: the per-channel rate
//!   is `(1 − √(1 − 4p/3))/2` for one qubit and `(1 − (1 − 16p/15)^{1/8})/2`
//!   for two. Mechanisms with the same symptoms are XOR-combined,
//!   `p ← p₁ + p₂ − 2p₁p₂`.
//! - **A non-deterministic detector is refused.** A detector or observable
//!   that anticommutes with a reset or a measurement would be random even
//!   without noise. The model reports such a circuit as an error rather than
//!   mis-scoring it.
//!
//! # Circuits
//!
//! [`parse`] reads the ecosystem's circuit text for the instructions an
//! error-correction experiment uses:
//! - resets and measurements in Z and X, with measurement flip noise;
//! - `H`, `S`, `√X`, Paulis, `CX`, `CZ`, `SWAP`;
//! - X/Y/Z flips and one- and two-qubit depolarizing noise;
//! - detectors, observables and `REPEAT` blocks.
//!
//! An instruction outside that set is an error, not a skipped line.
//!
//! [`surface_code_memory`] generates the rotated surface-code memory experiment.
//! It emits the same text, byte for byte, as the layout the ecosystem's
//! generators use, so models and decoders can be compared across tools on
//! identical circuits.

use crate::repro::ln;

/// A measurement or reset basis.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum Basis {
    X,
    Z,
}

/// A single-qubit noise channel.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum Noise1 {
    X,
    Y,
    Z,
    /// X, Y or Z, each with probability `p/3`, mutually exclusive.
    Depolarize,
}

/// One operation of a flattened circuit, on one qubit or one pair.
#[derive(Clone, Copy, Debug, PartialEq)]
pub enum Op {
    Reset { basis: Basis, q: u32 },
    /// Measure (and, if `reset`, reset) `q`, writing record `index`; the
    /// recorded result is flipped with probability `flip`.
    Measure { basis: Basis, reset: bool, flip: f64, q: u32, index: u32 },
    H(u32),
    /// `S` or `S†`: identical on Pauli frames.
    S(u32),
    /// `√X` or `√X†`.
    SqrtX(u32),
    Cx(u32, u32),
    Cz(u32, u32),
    Swap(u32, u32),
    Noise1 { kind: Noise1, p: f64, q: u32 },
    /// Each of the 15 non-identity two-qubit Paulis with probability `p/15`.
    Depolarize2 { p: f64, a: u32, b: u32 },
}

/// A circuit with its repeats unrolled and its measurement records resolved.
#[derive(Clone, Debug, PartialEq, Default)]
pub struct Circuit {
    pub qubits: u32,
    pub ops: Vec<Op>,
    pub measurements: u32,
    /// Each detector's measurement records (absolute indices).
    pub detectors: Vec<Vec<u32>>,
    /// Each observable's measurement records.
    pub observables: Vec<Vec<u32>>,
}

#[derive(Clone, Debug, PartialEq)]
pub enum FrameError {
    /// Text that is not a circuit this module reads, with its line.
    Parse { line: usize, message: String },
    /// A record reference before the first measurement.
    RecordOutOfRange { line: usize },
    /// A detector or observable that is random even without noise: it
    /// anticommutes with the reset or measurement at this flattened op.
    NonDeterministic { op: usize },
    /// A probability outside `[0, 1]` (or above `3/4` and `15/16` for the
    /// depolarizing channels, where they stop being mixtures).
    BadProbability(f64),
}

impl core::fmt::Display for FrameError {
    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
        match self {
            FrameError::Parse { line, message } => write!(f, "line {line}: {message}"),
            FrameError::RecordOutOfRange { line } => write!(f, "line {line}: measurement record out of range"),
            FrameError::NonDeterministic { op } => {
                write!(f, "a detector or observable is not deterministic (it anticommutes with op {op})")
            }
            FrameError::BadProbability(p) => write!(f, "probability {p} is out of range"),
        }
    }
}

// ---------------------------------------------------------------------------
// Parsing
// ---------------------------------------------------------------------------

/// Parse circuit text, unrolling `REPEAT` blocks and resolving `rec[-k]`.
pub fn parse(text: &str) -> Result<Circuit, FrameError> {
    let lines: Vec<(usize, String)> = text
        .lines()
        .enumerate()
        .map(|(i, l)| (i + 1, l.split('#').next().unwrap_or("").trim().to_string()))
        .filter(|(_, l)| !l.is_empty())
        .collect();
    let mut c = Circuit::default();
    let mut pos = 0;
    parse_block(&lines, &mut pos, &mut c, false)?;
    Ok(c)
}

fn parse_block(lines: &[(usize, String)], pos: &mut usize, c: &mut Circuit, nested: bool) -> Result<(), FrameError> {
    while *pos < lines.len() {
        let (line, text) = &lines[*pos];
        let line = *line;
        *pos += 1;
        if text == "}" {
            if nested {
                return Ok(());
            }
            return Err(FrameError::Parse { line, message: "unmatched '}'".into() });
        }
        if let Some(rest) = text.strip_prefix("REPEAT") {
            let rest = rest.trim();
            let count = rest
                .strip_suffix('{')
                .map(str::trim)
                .and_then(|n| n.parse::<u64>().ok())
                .ok_or(FrameError::Parse { line, message: "expected 'REPEAT <n> {'".into() })?;
            let start = *pos;
            // Parse the body once to find its end, then replay it.
            let mut probe = Circuit { qubits: c.qubits, measurements: c.measurements, ..Default::default() };
            parse_block(lines, pos, &mut probe, true)?;
            let end = *pos;
            for _ in 0..count {
                let mut p = start;
                parse_block(&lines[..end], &mut p, c, true)?;
            }
            continue;
        }
        parse_instruction(line, text, c)?;
    }
    if nested {
        return Err(FrameError::Parse { line: lines.last().map_or(0, |l| l.0), message: "unterminated REPEAT".into() });
    }
    Ok(())
}

fn parse_instruction(line: usize, text: &str, c: &mut Circuit) -> Result<(), FrameError> {
    let err = |m: &str| FrameError::Parse { line, message: m.into() };
    // `NAME(args) targets` or `NAME targets`; arguments may contain spaces.
    let (name, args, targets) = match text.find('(') {
        Some(open) if text[..open].chars().all(|ch| ch.is_ascii_alphanumeric() || ch == '_') => {
            let close = text.find(')').ok_or_else(|| err("unclosed '('"))?;
            (&text[..open], &text[open + 1..close], text[close + 1..].trim())
        }
        _ => match text.find(char::is_whitespace) {
            Some(i) => (&text[..i], "", text[i..].trim()),
            None => (text, "", ""),
        },
    };
    let nums: Vec<f64> = if args.trim().is_empty() {
        vec![]
    } else {
        args.split(',').map(|a| a.trim().parse::<f64>().map_err(|_| err("bad argument"))).collect::<Result<_, _>>()?
    };
    let prob = |i: usize| -> Result<f64, FrameError> {
        let p = *nums.get(i).ok_or_else(|| err("missing probability"))?;
        if !(0.0..=1.0).contains(&p) {
            return Err(FrameError::BadProbability(p));
        }
        Ok(p)
    };
    let qs = || -> Result<Vec<u32>, FrameError> {
        targets.split_whitespace().map(|t| t.parse::<u32>().map_err(|_| err("bad qubit target"))).collect()
    };
    let pairs = || -> Result<Vec<(u32, u32)>, FrameError> {
        let q = qs()?;
        if q.len() % 2 != 0 {
            return Err(err("two-qubit gate needs an even number of targets"));
        }
        if q.chunks(2).any(|p| p[0] == p[1]) {
            return Err(err("two-qubit gate on one qubit"));
        }
        Ok(q.chunks(2).map(|p| (p[0], p[1])).collect())
    };
    let recs = |c: &Circuit| -> Result<Vec<u32>, FrameError> {
        targets
            .split_whitespace()
            .map(|t| {
                let k = t
                    .strip_prefix("rec[-")
                    .and_then(|r| r.strip_suffix(']'))
                    .and_then(|r| r.parse::<u32>().ok())
                    .ok_or_else(|| err("expected rec[-k]"))?;
                if k == 0 || k > c.measurements {
                    return Err(FrameError::RecordOutOfRange { line });
                }
                Ok(c.measurements - k)
            })
            .collect()
    };
    let touch = |c: &mut Circuit, q: &[u32]| {
        if let Some(m) = q.iter().max() {
            c.qubits = c.qubits.max(m + 1);
        }
    };
    match name {
        "TICK" | "QUBIT_COORDS" | "SHIFT_COORDS" => {}
        "I" | "X" | "Y" | "Z" => touch(c, &qs()?),
        "R" | "RZ" | "RX" => {
            let basis = if name == "RX" { Basis::X } else { Basis::Z };
            let q = qs()?;
            touch(c, &q);
            c.ops.extend(q.into_iter().map(|q| Op::Reset { basis, q }));
        }
        "M" | "MZ" | "MX" | "MR" | "MRZ" | "MRX" => {
            let basis = if name.ends_with('X') { Basis::X } else { Basis::Z };
            let reset = name.starts_with("MR");
            let flip = if nums.is_empty() { 0.0 } else { prob(0)? };
            let q = qs()?;
            touch(c, &q);
            for q in q {
                c.ops.push(Op::Measure { basis, reset, flip, q, index: c.measurements });
                c.measurements += 1;
            }
        }
        "H" | "S" | "S_DAG" | "SQRT_X" | "SQRT_X_DAG" => {
            let q = qs()?;
            touch(c, &q);
            c.ops.extend(q.into_iter().map(|q| match name {
                "H" => Op::H(q),
                "S" | "S_DAG" => Op::S(q),
                _ => Op::SqrtX(q),
            }));
        }
        "CX" | "CNOT" | "ZCX" | "CZ" | "ZCZ" | "SWAP" => {
            let p = pairs()?;
            for &(a, b) in &p {
                touch(c, &[a, b]);
            }
            c.ops.extend(p.into_iter().map(|(a, b)| match name {
                "CZ" | "ZCZ" => Op::Cz(a, b),
                "SWAP" => Op::Swap(a, b),
                _ => Op::Cx(a, b),
            }));
        }
        "X_ERROR" | "Y_ERROR" | "Z_ERROR" | "DEPOLARIZE1" => {
            let p = prob(0)?;
            if name == "DEPOLARIZE1" && p > 0.75 {
                return Err(FrameError::BadProbability(p));
            }
            let kind = match name {
                "X_ERROR" => Noise1::X,
                "Y_ERROR" => Noise1::Y,
                "Z_ERROR" => Noise1::Z,
                _ => Noise1::Depolarize,
            };
            let q = qs()?;
            touch(c, &q);
            c.ops.extend(q.into_iter().map(|q| Op::Noise1 { kind, p, q }));
        }
        "DEPOLARIZE2" => {
            let p = prob(0)?;
            if p > 15.0 / 16.0 {
                return Err(FrameError::BadProbability(p));
            }
            let pr = pairs()?;
            for &(a, b) in &pr {
                touch(c, &[a, b]);
            }
            c.ops.extend(pr.into_iter().map(|(a, b)| Op::Depolarize2 { p, a, b }));
        }
        "DETECTOR" => {
            let r = recs(c)?;
            c.detectors.push(r);
        }
        "OBSERVABLE_INCLUDE" => {
            let k = *nums.first().ok_or_else(|| err("missing observable index"))?;
            if k < 0.0 || k.fract() != 0.0 || k > 63.0 {
                return Err(err("observable index must be 0..=63"));
            }
            let k = k as usize;
            let r = recs(c)?;
            if c.observables.len() <= k {
                c.observables.resize(k + 1, Vec::new());
            }
            c.observables[k].extend(r);
        }
        other => return Err(err(&format!("unsupported instruction '{other}'"))),
    }
    Ok(())
}

// ---------------------------------------------------------------------------
// The rotated surface-code memory experiment
// ---------------------------------------------------------------------------

/// Text of a rotated surface-code memory experiment (Z basis): distance `d`,
/// `rounds` syndrome rounds, every noise knob at `p`. The four knobs are:
/// - depolarizing noise after every Clifford;
/// - depolarizing noise on the data before each round;
/// - flips before measurements;
/// - flips after resets.
///
/// The layout, gate order and detector order follow the ecosystem's reference
/// generator, so the text matches it byte for byte.
pub fn surface_code_memory(d: u32, rounds: u32, p: f64) -> String {
    assert!(d >= 2 && rounds >= 1, "distance ≥ 2 and at least one round");
    let w = 2 * d + 1;
    let index = |x: u32, y: u32| x + (y / 2) * w;
    let mut data: Vec<(u32, u32)> = Vec::new();
    for x in 0..d {
        for y in 0..d {
            data.push((2 * x + 1, 2 * y + 1));
        }
    }
    // Measure qubits, in generation order (x outer, y inner), tagged X-type.
    let mut meas: Vec<((u32, u32), bool)> = Vec::new();
    for x in 0..=d {
        for y in 0..=d {
            let edge1 = x == 0 || x == d;
            let edge2 = y == 0 || y == d;
            let parity = (x % 2) != (y % 2);
            if (edge1 && parity) || (edge2 && !parity) {
                continue;
            }
            meas.push(((2 * x, 2 * y), parity));
        }
    }
    let is_data = |x: i64, y: i64| x > 0 && y > 0 && x < 2 * d as i64 && y < 2 * d as i64 && x % 2 == 1 && y % 2 == 1;
    let mut all: Vec<(u32, (u32, u32))> = data.iter().chain(meas.iter().map(|(c, _)| c)).map(|&(x, y)| (index(x, y), (x, y))).collect();
    all.sort();
    let mut data_ix: Vec<u32> = data.iter().map(|&(x, y)| index(x, y)).collect();
    data_ix.sort();
    let mut meas_ix: Vec<u32> = meas.iter().map(|&((x, y), _)| index(x, y)).collect();
    meas_ix.sort();
    let mut x_ix: Vec<u32> = meas.iter().filter(|m| m.1).map(|&((x, y), _)| index(x, y)).collect();
    x_ix.sort();
    let x_order: [(i64, i64); 4] = [(1, 1), (-1, 1), (1, -1), (-1, -1)];
    let z_order: [(i64, i64); 4] = [(1, 1), (1, -1), (-1, 1), (-1, -1)];
    let mut layers: Vec<Vec<u32>> = Vec::new();
    for k in 0..4 {
        let mut targets = Vec::new();
        for want_x in [true, false] {
            for &((mx, my), is_x) in &meas {
                if is_x != want_x {
                    continue;
                }
                let (dx, dy) = if is_x { x_order[k] } else { z_order[k] };
                let (tx, ty) = (mx as i64 + dx, my as i64 + dy);
                if !is_data(tx, ty) {
                    continue;
                }
                let (m, t) = (index(mx, my), index(tx as u32, ty as u32));
                if is_x {
                    targets.extend([m, t]);
                } else {
                    targets.extend([t, m]);
                }
            }
        }
        layers.push(targets);
    }
    let list = |v: &[u32]| v.iter().map(|q| q.to_string()).collect::<Vec<_>>().join(" ");
    let fp = fmt_prob(p);
    let mut o = String::new();
    for (q, (x, y)) in &all {
        o += &format!("QUBIT_COORDS({x}, {y}) {q}\n");
    }
    let rec_of = |q: u32, list: &[u32]| -> i64 { list.iter().position(|&m| m == q).unwrap() as i64 - list.len() as i64 };
    let round = |o: &mut String, indent: &str| {
        *o += &format!("{indent}DEPOLARIZE1({fp}) {}\n", list(&data_ix));
        *o += &format!("{indent}H {}\n", list(&x_ix));
        *o += &format!("{indent}DEPOLARIZE1({fp}) {}\n", list(&x_ix));
        *o += &format!("{indent}TICK\n");
        for l in &layers {
            *o += &format!("{indent}CX {}\n", list(l));
            *o += &format!("{indent}DEPOLARIZE2({fp}) {}\n", list(l));
            *o += &format!("{indent}TICK\n");
        }
        *o += &format!("{indent}H {}\n", list(&x_ix));
        *o += &format!("{indent}DEPOLARIZE1({fp}) {}\n", list(&x_ix));
        *o += &format!("{indent}TICK\n");
        *o += &format!("{indent}X_ERROR({fp}) {}\n", list(&meas_ix));
        *o += &format!("{indent}MR {}\n", list(&meas_ix));
        *o += &format!("{indent}X_ERROR({fp}) {}\n", list(&meas_ix));
    };
    o += &format!("R {}\nX_ERROR({fp}) {}\n", list(&data_ix), list(&data_ix));
    o += &format!("R {}\nX_ERROR({fp}) {}\nTICK\n", list(&meas_ix), list(&meas_ix));
    round(&mut o, "");
    for &((x, y), is_x) in &meas {
        if !is_x {
            o += &format!("DETECTOR({x}, {y}, 0) rec[{}]\n", rec_of(index(x, y), &meas_ix));
        }
    }
    if rounds > 1 {
        // A single repeated round is written inline, not as `REPEAT 1`.
        let indent = if rounds > 2 { "    " } else { "" };
        if rounds > 2 {
            o += &format!("REPEAT {} {{\n", rounds - 1);
        }
        o += &format!("{indent}TICK\n");
        round(&mut o, indent);
        o += &format!("{indent}SHIFT_COORDS(0, 0, 1)\n");
        let n = meas_ix.len() as i64;
        for &q in &meas_ix {
            let (x, y) = all.iter().find(|a| a.0 == q).unwrap().1;
            let r = rec_of(q, &meas_ix);
            o += &format!("{indent}DETECTOR({x}, {y}, 0) rec[{r}] rec[{}]\n", r - n);
        }
        if rounds > 2 {
            o += "}\n";
        }
    }
    o += &format!("X_ERROR({fp}) {}\nM {}\n", list(&data_ix), list(&data_ix));
    let nd = data_ix.len() as i64;
    for &((x, y), is_x) in &meas {
        if is_x {
            continue;
        }
        let mut r: Vec<i64> = [(1, 1), (1, -1), (-1, 1), (-1, -1)]
            .iter()
            .map(|&(dx, dy)| (x as i64 + dx, y as i64 + dy))
            .filter(|&(tx, ty)| is_data(tx, ty))
            .map(|(tx, ty)| rec_of(index(tx as u32, ty as u32), &data_ix))
            .collect();
        r.sort_by(|a, b| b.cmp(a));
        let anc = rec_of(index(x, y), &meas_ix) - nd;
        let recs: Vec<String> = r.iter().chain(core::iter::once(&anc)).map(|k| format!("rec[{k}]")).collect();
        o += &format!("DETECTOR({x}, {y}, 1) {}\n", recs.join(" "));
    }
    let mut obs: Vec<i64> = (0..d).map(|x| rec_of(index(2 * x + 1, 1), &data_ix)).collect();
    obs.sort_by(|a, b| b.cmp(a));
    o += &format!("OBSERVABLE_INCLUDE(0) {}\n", obs.iter().map(|k| format!("rec[{k}]")).collect::<Vec<_>>().join(" "));
    o
}

/// The shortest decimal that reads back as `p`.
fn fmt_prob(p: f64) -> String {
    format!("{p}")
}

// ---------------------------------------------------------------------------
// Pauli-frame sampling
// ---------------------------------------------------------------------------

struct Rng(u64);

impl Rng {
    fn next(&mut self) -> u64 {
        self.0 = self.0.wrapping_add(0x9e37_79b9_7f4a_7c15);
        let mut z = self.0;
        z = (z ^ (z >> 30)).wrapping_mul(0xbf58_476d_1ce4_e5b9);
        z = (z ^ (z >> 27)).wrapping_mul(0x94d0_49bb_1331_11eb);
        z ^ (z >> 31)
    }
    /// Uniform on `(0, 1]`.
    fn unit(&mut self) -> f64 {
        ((self.next() >> 11) + 1) as f64 * (1.0 / 9_007_199_254_740_992.0)
    }
    /// Uniform on `0..n` (n small; the bias is below 2⁻⁶⁰).
    fn below(&mut self, n: u64) -> u64 {
        self.next() % n
    }
}

/// Call `hit` for every position in `0..n` that a rate-`p` Bernoulli trial
/// selects, skipping geometrically between hits.
fn bernoulli(rng: &mut Rng, p: f64, n: usize, mut hit: impl FnMut(usize, &mut Rng)) {
    if p <= 0.0 || n == 0 {
        return;
    }
    if p >= 1.0 {
        for i in 0..n {
            hit(i, rng);
        }
        return;
    }
    let lq = ln(1.0 - p);
    let mut i = 0usize;
    loop {
        let skip = ln(rng.unit()) / lq;
        if skip >= (n - i) as f64 {
            return;
        }
        i += skip as usize;
        hit(i, rng);
        i += 1;
        if i >= n {
            return;
        }
    }
}

/// Shots per batch, in 64-bit words. Part of the reproducibility contract:
/// the random stream is consumed batch by batch.
const BATCH_WORDS: usize = 8;

/// Detection events and observable flips for a number of shots, bit-packed.
#[derive(Clone, Debug, PartialEq)]
pub struct Detections {
    pub shots: usize,
    pub detectors: usize,
    pub observables: usize,
    words: usize,
    /// `det[d * words + w]`: bit `s` of word `w` is shot `64w + s`.
    det: Vec<u64>,
    obs: Vec<u64>,
}

impl Detections {
    /// The detectors that fired in shot `s`, ascending.
    pub fn fired(&self, s: usize) -> Vec<u32> {
        let (w, b) = (s / 64, s % 64);
        (0..self.detectors).filter(|d| (self.det[d * self.words + w] >> b) & 1 == 1).map(|d| d as u32).collect()
    }
    /// The observables that flipped in shot `s`, as a bitmask.
    pub fn flips(&self, s: usize) -> u64 {
        let (w, b) = (s / 64, s % 64);
        (0..self.observables).fold(0, |m, k| m | (((self.obs[k * self.words + w] >> b) & 1) << k))
    }
    /// How many shots fired detector `d`.
    pub fn count(&self, d: usize) -> u64 {
        let row = &self.det[d * self.words..(d + 1) * self.words];
        let full = self.shots / 64;
        let mut n: u64 = row[..full].iter().map(|w| w.count_ones() as u64).sum();
        if !self.shots.is_multiple_of(64) {
            n += (row[full] & ((1u64 << (self.shots % 64)) - 1)).count_ones() as u64;
        }
        n
    }
}

/// Sample `shots` runs of `c`: which detectors fired, which observables
/// flipped. Reproducible from `(c, shots, seed)`.
pub fn sample(c: &Circuit, shots: usize, seed: u64) -> Detections {
    let words = shots.div_ceil(64);
    let (nd, no) = (c.detectors.len(), c.observables.len());
    let mut det = vec![0u64; nd * words];
    let mut obs = vec![0u64; no * words];
    let mut rng = Rng(seed);
    let nq = c.qubits as usize;
    let nm = c.measurements as usize;
    let mut w0 = 0;
    while w0 < words {
        let bw = BATCH_WORDS.min(words - w0);
        let nshots = (shots - 64 * w0).min(64 * bw);
        let mut x = vec![0u64; nq * bw];
        let mut z = vec![0u64; nq * bw];
        let mut m = vec![0u64; nm * bw];
        for op in &c.ops {
            run_op(op, &mut x, &mut z, &mut m, bw, nshots, &mut rng);
        }
        for (d, recs) in c.detectors.iter().enumerate() {
            for w in 0..bw {
                det[d * words + w0 + w] = recs.iter().fold(0, |a, &r| a ^ m[r as usize * bw + w]);
            }
        }
        for (k, recs) in c.observables.iter().enumerate() {
            for w in 0..bw {
                obs[k * words + w0 + w] = recs.iter().fold(0, |a, &r| a ^ m[r as usize * bw + w]);
            }
        }
        w0 += bw;
    }
    Detections { shots, detectors: nd, observables: no, words, det, obs }
}

fn run_op(op: &Op, x: &mut [u64], z: &mut [u64], m: &mut [u64], bw: usize, nshots: usize, rng: &mut Rng) {
    let row = |q: u32| q as usize * bw..(q as usize + 1) * bw;
    let flip_bit = |v: &mut [u64], q: u32, s: usize| v[q as usize * bw + s / 64] ^= 1u64 << (s % 64);
    match *op {
        Op::Reset { basis, q } => {
            let (keep, rand) = match basis {
                Basis::Z => (&mut *x, &mut *z),
                Basis::X => (&mut *z, &mut *x),
            };
            for i in row(q) {
                keep[i] = 0;
                rand[i] = rng.next();
            }
        }
        Op::Measure { basis, reset, flip, q, index } => {
            let mrow = index as usize * bw;
            let (seen, gauge) = match basis {
                Basis::Z => (&mut *x, &mut *z),
                Basis::X => (&mut *z, &mut *x),
            };
            for (k, i) in row(q).enumerate() {
                m[mrow + k] = seen[i];
                if reset {
                    seen[i] = 0;
                }
                gauge[i] = rng.next();
            }
            bernoulli(rng, flip, nshots, |s, _| m[mrow + s / 64] ^= 1u64 << (s % 64));
        }
        Op::H(q) => {
            for i in row(q) {
                core::mem::swap(&mut x[i], &mut z[i]);
            }
        }
        Op::S(q) => {
            for i in row(q) {
                z[i] ^= x[i];
            }
        }
        Op::SqrtX(q) => {
            for i in row(q) {
                x[i] ^= z[i];
            }
        }
        Op::Cx(a, b) => {
            for k in 0..bw {
                let (ia, ib) = (a as usize * bw + k, b as usize * bw + k);
                x[ib] ^= x[ia];
                z[ia] ^= z[ib];
            }
        }
        Op::Cz(a, b) => {
            for k in 0..bw {
                let (ia, ib) = (a as usize * bw + k, b as usize * bw + k);
                z[ia] ^= x[ib];
                z[ib] ^= x[ia];
            }
        }
        Op::Swap(a, b) => {
            for k in 0..bw {
                let (ia, ib) = (a as usize * bw + k, b as usize * bw + k);
                x.swap(ia, ib);
                z.swap(ia, ib);
            }
        }
        Op::Noise1 { kind, p, q } => bernoulli(rng, p, nshots, |s, r| {
            let pauli = match kind {
                Noise1::X => 1,
                Noise1::Y => 2,
                Noise1::Z => 3,
                Noise1::Depolarize => 1 + r.below(3),
            };
            if pauli != 3 {
                flip_bit(x, q, s);
            }
            if pauli != 1 {
                flip_bit(z, q, s);
            }
        }),
        Op::Depolarize2 { p, a, b } => bernoulli(rng, p, nshots, |s, r| {
            let k = 1 + r.below(15);
            for (q, pauli) in [(a, k >> 2), (b, k & 3)] {
                if pauli == 1 || pauli == 2 {
                    flip_bit(x, q, s);
                }
                if pauli == 2 || pauli == 3 {
                    flip_bit(z, q, s);
                }
            }
        }),
    }
}

// ---------------------------------------------------------------------------
// The detector error model
// ---------------------------------------------------------------------------

/// One independent way the experiment can fail.
#[derive(Clone, Debug, PartialEq)]
pub struct Mechanism {
    pub probability: f64,
    /// Detectors it flips, ascending.
    pub detectors: Vec<u32>,
    /// Observables it flips, as a bitmask.
    pub observables: u64,
    /// How a matching decoder sees it: the whole symptom set when it touches
    /// at most two detectors, else its X and Z halves, else its per-qubit
    /// components. A matching decoder needs every part to touch at most two
    /// detectors.
    pub parts: Vec<(Vec<u32>, u64)>,
}

/// The circuit's detector error model.
#[derive(Clone, Debug, PartialEq)]
pub struct ErrorModel {
    pub detectors: u32,
    pub observables: u32,
    /// Merged by symptoms, sorted by (detectors, observables).
    pub mechanisms: Vec<Mechanism>,
}

/// Symmetric difference of two ascending sets, into `a`.
fn xor_into(a: &mut Vec<u32>, b: &[u32]) {
    if b.is_empty() {
        return;
    }
    let mut out = Vec::with_capacity(a.len() + b.len());
    let (mut i, mut j) = (0, 0);
    while i < a.len() && j < b.len() {
        match a[i].cmp(&b[j]) {
            core::cmp::Ordering::Less => {
                out.push(a[i]);
                i += 1;
            }
            core::cmp::Ordering::Greater => {
                out.push(b[j]);
                j += 1;
            }
            core::cmp::Ordering::Equal => {
                i += 1;
                j += 1;
            }
        }
    }
    out.extend_from_slice(&a[i..]);
    out.extend_from_slice(&b[j..]);
    *a = out;
}

/// Per-Pauli rate of independent channels equivalent to a depolarizing one.
fn independent_rate(p: f64, qubits: u32) -> f64 {
    if qubits == 1 {
        (1.0 - (1.0 - 4.0 * p / 3.0).sqrt()) / 2.0
    } else {
        (1.0 - (1.0 - 16.0 * p / 15.0).sqrt().sqrt().sqrt()) / 2.0
    }
}

/// Derive the detector error model by backward sensitivity propagation.
pub fn error_model(c: &Circuit) -> Result<ErrorModel, FrameError> {
    let nd = c.detectors.len() as u32;
    // Symptom ids: detectors 0..nd, observable k as nd + k.
    let mut of_meas: Vec<Vec<u32>> = vec![Vec::new(); c.measurements as usize];
    for (d, recs) in c.detectors.iter().enumerate() {
        for &r in recs {
            xor_into(&mut of_meas[r as usize], &[d as u32]);
        }
    }
    for (k, recs) in c.observables.iter().enumerate() {
        for &r in recs {
            xor_into(&mut of_meas[r as usize], &[nd + k as u32]);
        }
    }
    let nq = c.qubits as usize;
    let mut xs: Vec<Vec<u32>> = vec![Vec::new(); nq];
    let mut zs: Vec<Vec<u32>> = vec![Vec::new(); nq];
    // (symptom set, probability, parts) in discovery order.
    let mut found: Found = Vec::new();
    for (i, op) in c.ops.iter().enumerate().rev() {
        match *op {
            Op::Reset { basis, q } => {
                let q = q as usize;
                let gauge = if basis == Basis::Z { &zs[q] } else { &xs[q] };
                if !gauge.is_empty() {
                    return Err(FrameError::NonDeterministic { op: i });
                }
                xs[q].clear();
                zs[q].clear();
            }
            Op::Measure { basis, reset, flip, q, index } => {
                let q = q as usize;
                if reset {
                    let gauge = if basis == Basis::Z { &zs[q] } else { &xs[q] };
                    if !gauge.is_empty() {
                        return Err(FrameError::NonDeterministic { op: i });
                    }
                    xs[q].clear();
                    zs[q].clear();
                }
                let (seen, gauge) = if basis == Basis::Z { (&mut xs, &zs) } else { (&mut zs, &xs) };
                if !gauge[q].is_empty() {
                    return Err(FrameError::NonDeterministic { op: i });
                }
                let s = &of_meas[index as usize];
                if flip > 0.0 && !s.is_empty() {
                    found.push((s.clone(), flip, vec![s.clone()]));
                }
                xor_into(&mut seen[q], s);
            }
            Op::H(q) => {
                let q = q as usize;
                core::mem::swap(&mut xs[q], &mut zs[q]);
            }
            Op::S(q) => {
                let q = q as usize;
                let t = zs[q].clone();
                xor_into(&mut xs[q], &t);
            }
            Op::SqrtX(q) => {
                let q = q as usize;
                let t = xs[q].clone();
                xor_into(&mut zs[q], &t);
            }
            Op::Cx(a, b) => {
                let (a, b) = (a as usize, b as usize);
                let t = xs[b].clone();
                xor_into(&mut xs[a], &t);
                let t = zs[a].clone();
                xor_into(&mut zs[b], &t);
            }
            Op::Cz(a, b) => {
                let (a, b) = (a as usize, b as usize);
                let (za, zb) = (zs[a].clone(), zs[b].clone());
                xor_into(&mut xs[a], &zb);
                xor_into(&mut xs[b], &za);
            }
            Op::Swap(a, b) => {
                let (a, b) = (a as usize, b as usize);
                xs.swap(a, b);
                zs.swap(a, b);
            }
            Op::Noise1 { kind, p, q } => {
                if p == 0.0 {
                    continue;
                }
                let q = q as usize;
                let (terms, n) = match kind {
                    Noise1::X => ([(1, p), (0, 0.0), (0, 0.0)], 1),
                    Noise1::Y => ([(2, p), (0, 0.0), (0, 0.0)], 1),
                    Noise1::Z => ([(3, p), (0, 0.0), (0, 0.0)], 1),
                    Noise1::Depolarize => {
                        let r = independent_rate(p, 1);
                        ([(1, r), (2, r), (3, r)], 3)
                    }
                };
                for &(pauli, r) in &terms[..n] {
                    push_found(&mut found, pauli_parts(&[(q, pauli)], &xs, &zs), r, nd);
                }
            }
            Op::Depolarize2 { p, a, b } => {
                if p == 0.0 {
                    continue;
                }
                let r = independent_rate(p, 2);
                for k in 1..16u8 {
                    push_found(&mut found, pauli_parts(&[(a as usize, k >> 2), (b as usize, k & 3)], &xs, &zs), r, nd);
                }
            }
        }
    }
    // Merge by full symptom set (XOR-combine), keeping the first parts seen.
    found.sort_by(|a, b| a.0.cmp(&b.0));
    let mut merged: Vec<(Vec<u32>, f64, Vec<Vec<u32>>)> = Vec::new();
    for (s, p, parts) in found {
        match merged.last_mut() {
            Some(last) if last.0 == s => last.1 = last.1 + p - 2.0 * last.1 * p,
            _ => merged.push((s, p, parts)),
        }
    }
    let split = |s: &[u32]| -> (Vec<u32>, u64) {
        let dets: Vec<u32> = s.iter().copied().filter(|&x| x < nd).collect();
        let obs = s.iter().filter(|&&x| x >= nd).fold(0u64, |m, &x| m | 1 << (x - nd));
        (dets, obs)
    };
    let mut mechanisms: Vec<Mechanism> = merged
        .into_iter()
        .map(|(s, p, parts)| {
            let (detectors, observables) = split(&s);
            Mechanism { probability: p, detectors, observables, parts: parts.iter().map(|x| split(x)).collect() }
        })
        .collect();
    mechanisms.sort_by(|a, b| (&a.detectors, a.observables).cmp(&(&b.detectors, b.observables)));
    Ok(ErrorModel { detectors: nd, observables: c.observables.len() as u32, mechanisms })
}

/// Errors found so far: (symptom set, probability, parts for matching).
type Found = Vec<(Vec<u32>, f64, Vec<Vec<u32>>)>;

/// An error's X components and Z components, each as the symptom set it flips.
type Halves = (Vec<Vec<u32>>, Vec<Vec<u32>>);

/// The X and Z components of a Pauli on some qubits, each as the symptom set
/// it flips: `(x components, z components)`, empty ones dropped.
fn pauli_parts(terms: &[(usize, u8)], xs: &[Vec<u32>], zs: &[Vec<u32>]) -> Halves {
    let (mut xp, mut zp) = (Vec::new(), Vec::new());
    for &(q, pauli) in terms {
        if (pauli == 1 || pauli == 2) && !xs[q].is_empty() {
            xp.push(xs[q].clone());
        }
        if (pauli == 2 || pauli == 3) && !zs[q].is_empty() {
            zp.push(zs[q].clone());
        }
    }
    (xp, zp)
}

/// Record an error and how a matching decoder should see it: the whole
/// symptom if it touches at most two detectors; otherwise its X and Z halves
/// if each does; otherwise its per-qubit components.
fn push_found(found: &mut Found, (xp, zp): Halves, p: f64, nd: u32) {
    let fold = |parts: &[Vec<u32>]| {
        let mut s = Vec::new();
        for part in parts {
            xor_into(&mut s, part);
        }
        s
    };
    let (xt, zt) = (fold(&xp), fold(&zp));
    let mut full = xt.clone();
    xor_into(&mut full, &zt);
    if full.is_empty() {
        return;
    }
    let dets = |s: &[u32]| s.iter().filter(|&&x| x < nd).count();
    let parts = if dets(&full) <= 2 {
        vec![full.clone()]
    } else if dets(&xt) <= 2 && dets(&zt) <= 2 {
        [xt, zt].into_iter().filter(|s| !s.is_empty()).collect()
    } else {
        xp.into_iter().chain(zp).collect()
    };
    found.push((full, p, parts));
}

impl ErrorModel {
    /// The model in the ecosystem's text syntax: one `error(p) D… L…` line per
    /// mechanism, with `p` printed as the shortest decimal that reads back
    /// exactly.
    pub fn to_text(&self) -> String {
        let mut o = String::new();
        for m in &self.mechanisms {
            o += &format!("error({})", m.probability);
            for d in &m.detectors {
                o += &format!(" D{d}");
            }
            for k in 0..64 {
                if m.observables >> k & 1 == 1 {
                    o += &format!(" L{k}");
                }
            }
            o.push('\n');
        }
        o
    }

    /// Mechanisms with a part that touches more than two detectors: the
    /// errors a matching decoder cannot represent.
    pub fn hyperedges(&self) -> usize {
        self.mechanisms.iter().filter(|m| m.parts.iter().any(|(d, _)| d.len() > 2)).count()
    }
}

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

    #[test]
    fn xor_is_symmetric_difference() {
        let mut a = vec![1, 3, 5, 9];
        xor_into(&mut a, &[0, 3, 9, 10]);
        assert_eq!(a, vec![0, 1, 5, 10]);
        xor_into(&mut a, &[0, 1, 5, 10]);
        assert!(a.is_empty());
    }

    #[test]
    fn the_independent_rates_reproduce_the_depolarizing_channels() {
        // One qubit: P(net X) = q(1 - q) for independent X, Y, Z at q.
        let p = 0.01;
        let q = independent_rate(p, 1);
        assert!((q * (1.0 - q) - p / 3.0).abs() < 1e-15);
        // Two qubits: P(net P) = (1 - (1 - 2q)^8) / 16 for each of 15.
        let q = independent_rate(p, 2);
        let net = (1.0 - (1.0 - 2.0 * q).powi(8)) / 16.0;
        assert!((net - p / 15.0).abs() < 1e-15);
    }

    #[test]
    fn a_repetition_code_has_the_expected_model() {
        // Three data qubits, two Z-parity checks, one round, X flips on data.
        let text = "R 0 1 2 3 4\nX_ERROR(0.1) 0 1 2\nCX 0 3 1 3 1 4 2 4\nMR 3 4\nM 0 1 2\n\
                    DETECTOR rec[-5]\nDETECTOR rec[-4]\nDETECTOR rec[-3] rec[-2] rec[-5]\n\
                    DETECTOR rec[-2] rec[-1] rec[-4]\nOBSERVABLE_INCLUDE(0) rec[-1]\n";
        let c = parse(text).unwrap();
        assert_eq!((c.qubits, c.measurements, c.detectors.len()), (5, 5, 4));
        let m = error_model(&c).unwrap();
        let got: Vec<(Vec<u32>, u64)> = m.mechanisms.iter().map(|m| (m.detectors.clone(), m.observables)).collect();
        // Qubit 0 flips check 0; qubit 1 both checks; qubit 2 check 1 and L0.
        assert_eq!(got, vec![(vec![0], 0), (vec![0, 1], 0), (vec![1], 1)]);
        assert!(m.mechanisms.iter().all(|x| x.probability == 0.1));
    }

    #[test]
    fn a_non_deterministic_detector_is_refused() {
        // An X-basis reset measured in Z: random without any noise.
        let c = parse("RX 0\nM 0\nDETECTOR rec[-1]\n").unwrap();
        assert!(matches!(error_model(&c), Err(FrameError::NonDeterministic { .. })));
    }

    #[test]
    fn unsupported_text_is_an_error_not_a_skip() {
        assert!(matches!(parse("T 0\n"), Err(FrameError::Parse { .. })));
        assert!(matches!(parse("M 0\nDETECTOR rec[-2]\n"), Err(FrameError::RecordOutOfRange { .. })));
        assert!(matches!(parse("REPEAT 2 {\nH 0\n"), Err(FrameError::Parse { .. })));
        assert!(matches!(parse("X_ERROR(1.5) 0\n"), Err(FrameError::BadProbability(_))));
    }

    #[test]
    fn repeat_blocks_unroll() {
        let c = parse("R 0\nREPEAT 3 {\n  X_ERROR(0.5) 0\n  MR 0\n  DETECTOR rec[-1]\n}\n").unwrap();
        assert_eq!((c.measurements, c.detectors.len()), (3, 3));
        assert_eq!(c.detectors, vec![vec![0], vec![1], vec![2]]);
    }

    #[test]
    fn sampling_matches_the_model_marginals_and_is_reproducible() {
        let c = parse(&surface_code_memory(3, 3, 0.01)).unwrap();
        let m = error_model(&c).unwrap();
        let shots = 200_000;
        let a = sample(&c, shots, 7);
        assert_eq!(a, sample(&c, shots, 7));
        assert_ne!(a, sample(&c, shots, 8));
        // A detector fires iff an odd number of its mechanisms fire.
        for d in 0..m.detectors {
            let mut q = 0.0;
            for mech in m.mechanisms.iter().filter(|x| x.detectors.contains(&d)) {
                q = q + mech.probability - 2.0 * q * mech.probability;
            }
            let rate = a.count(d as usize) as f64 / shots as f64;
            let sigma = (q * (1.0 - q) / shots as f64).sqrt();
            assert!((rate - q).abs() < 5.0 * sigma, "D{d}: sampled {rate} vs model {q}");
        }
    }

    /// Random Clifford layers over every gate type, noise in the middle, then
    /// the layers undone: every measurement is deterministic, so every one is a
    /// detector. The sampler applies the forward rules and the model the
    /// backward ones, so their marginals agree only if both are right.
    #[test]
    fn every_gate_agrees_between_sampler_and_model() {
        let mut state = 0x5eed_u64;
        let mut next = |m: u32| {
            state = state.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
            ((state >> 33) % m as u64) as u32
        };
        let n = 6;
        for trial in 0..4 {
            let mut layers: Vec<(String, String)> = Vec::new();
            for _ in 0..12 {
                let (a, b) = (next(n), next(n));
                let b = if a == b { (b + 1) % n } else { b };
                let (fwd, inv) = match next(7) {
                    0 => (format!("H {a}"), format!("H {a}")),
                    1 => (format!("S {a}"), format!("S_DAG {a}")),
                    2 => (format!("SQRT_X {a}"), format!("SQRT_X_DAG {a}")),
                    3 => (format!("CX {a} {b}"), format!("CX {a} {b}")),
                    4 => (format!("CZ {a} {b}"), format!("CZ {a} {b}")),
                    5 => (format!("SWAP {a} {b}"), format!("SWAP {a} {b}")),
                    _ => (format!("S_DAG {a}"), format!("S {a}")),
                };
                layers.push((fwd, inv));
            }
            let qs: Vec<String> = (0..n).map(|q| q.to_string()).collect();
            let mut t = format!("R {}\n", qs.join(" "));
            for (f, _) in &layers {
                t += &format!("{f}\nDEPOLARIZE1(0.01) {}\n", qs.join(" "));
            }
            t += "DEPOLARIZE2(0.02) 0 1 2 3 4 5\nY_ERROR(0.03) 2\nZ_ERROR(0.03) 4\n";
            for (_, i) in layers.iter().rev() {
                t += &format!("{i}\nX_ERROR(0.005) {}\n", qs.join(" "));
            }
            t += &format!("M(0.01) {}\n", qs.join(" "));
            for k in 1..=n {
                t += &format!("DETECTOR rec[-{k}]\n");
            }
            t += "OBSERVABLE_INCLUDE(0) rec[-1] rec[-2]\n";
            let c = parse(&t).unwrap();
            let m = error_model(&c).unwrap();
            let shots = 100_000;
            let smp = sample(&c, shots, 3 + trial);
            for d in 0..m.detectors {
                let mut q = 0.0;
                for mech in m.mechanisms.iter().filter(|x| x.detectors.contains(&d)) {
                    q = q + mech.probability - 2.0 * q * mech.probability;
                }
                let rate = smp.count(d as usize) as f64 / shots as f64;
                let sigma = (q * (1.0 - q) / shots as f64).sqrt().max(1e-4);
                assert!((rate - q).abs() < 5.0 * sigma, "trial {trial} D{d}: sampled {rate} vs model {q}\n{t}");
            }
            // The observable too.
            let mut q = 0.0;
            for mech in m.mechanisms.iter().filter(|x| x.observables & 1 == 1) {
                q = q + mech.probability - 2.0 * q * mech.probability;
            }
            let rate = (0..shots).filter(|&s| smp.flips(s) & 1 == 1).count() as f64 / shots as f64;
            assert!((rate - q).abs() < 5.0 * (q * (1.0 - q) / shots as f64).sqrt(), "trial {trial} L0: {rate} vs {q}");
        }
    }

    #[test]
    fn a_graphlike_error_stays_one_edge() {
        // A Y error inside a Bell pair: its Z half flips D0 and its X half
        // flips D1. Two detectors in all, so a matching decoder must see one
        // D0–D1 edge, not two boundary edges that lose the correlation.
        let c = parse("R 0 1\nH 0\nCX 0 1\nY_ERROR(0.1) 0\nCX 0 1\nH 0\nM 0 1\nDETECTOR rec[-2]\nDETECTOR rec[-1]\n").unwrap();
        let m = error_model(&c).unwrap();
        assert_eq!(m.mechanisms.len(), 1);
        assert_eq!(m.mechanisms[0].detectors, vec![0, 1]);
        assert_eq!(m.mechanisms[0].parts, vec![(vec![0, 1], 0)]);
    }

    #[test]
    fn noiseless_runs_fire_nothing() {
        let c = parse(&surface_code_memory(5, 4, 0.0)).unwrap();
        let s = sample(&c, 1000, 1);
        assert!((0..s.shots).all(|i| s.fired(i).is_empty() && s.flips(i) == 0));
        assert!(error_model(&c).unwrap().mechanisms.is_empty());
    }

    #[test]
    fn the_surface_code_model_is_graphlike() {
        for d in [3, 5] {
            let m = error_model(&parse(&surface_code_memory(d, d, 0.001)).unwrap()).unwrap();
            assert_eq!(m.hyperedges(), 0, "d = {d}");
        }
    }
}