diff --git a/wai/wai-quantum-serve/Cargo.lock b/wai/wai-quantum-serve/Cargo.lock index 5e9fe6ed9..34b20486b 100644 --- a/wai/wai-quantum-serve/Cargo.lock +++ b/wai/wai-quantum-serve/Cargo.lock @@ -589,7 +589,7 @@ dependencies = [ [[package]] name = "wai-quantum" -version = "0.3.24" +version = "0.3.25" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum/Cargo.lock b/wai/wai-quantum/Cargo.lock index e382e5e30..c2286cbf0 100644 --- a/wai/wai-quantum/Cargo.lock +++ b/wai/wai-quantum/Cargo.lock @@ -415,7 +415,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" [[package]] name = "wai-quantum" -version = "0.3.24" +version = "0.3.25" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum/Cargo.toml b/wai/wai-quantum/Cargo.toml index 12b401b8e..cfc38cdda 100644 --- a/wai/wai-quantum/Cargo.toml +++ b/wai/wai-quantum/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "wai-quantum" -version = "0.3.24" +version = "0.3.25" edition = "2024" rust-version = "1.97" license = "Apache-2.0" @@ -53,6 +53,9 @@ quantum_bptn = ["quantum_spd"] quantum_resource = [] # Signed, energy-accounted receipts for it, re-checkable bit for bit anywhere. quantum_resource_receipt = ["quantum_resource", "_signing"] +# Noisy Clifford circuits at QEC scale: bit-packed Pauli-frame sampling, detector +# error models by backward sensitivity, and exact minimum-weight matching. +quantum_frame = [] # Mixed states and the observables that distinguish quantum theory from a classical # account: density matrices, the partial trace, exact Pauli expectations, CHSH. quantum_info = ["quantum"] @@ -92,13 +95,21 @@ quantum_phasor_meter = ["_signing"] # Everything. full = [ - "quantum", "quantum_toolchain", "quantum_mps", "quantum_stabilizer", "quantum_pauli", "quantum_spd", "quantum_spd_receipt", "quantum_bptn", "quantum_resource", "quantum_resource_receipt", "quantum_qir", "quantum_info", "quantum_qudit", "quantum_comm", "quantum_source", "quantum_channel", + "quantum", "quantum_toolchain", "quantum_mps", "quantum_stabilizer", "quantum_pauli", "quantum_spd", "quantum_spd_receipt", "quantum_bptn", "quantum_resource", "quantum_resource_receipt", "quantum_frame", "quantum_qir", "quantum_info", "quantum_qudit", "quantum_comm", "quantum_source", "quantum_channel", "quantum_receipt", "quantum_ops", "quantum_cal", "quantum_mitigate", "quantum_qec", "quantum_control", "quantum_noise", "quantum_compile", "quantum_atom", "quantum_qbom", "quantum_vml", "quantum_phasor", "quantum_qfhrr", "quantum_kernel", "quantum_qdata", "quantum_phasor_meter", ] +[[example]] +name = "qec_threshold" +required-features = ["quantum_frame"] + +[[example]] +name = "qec_frame" +required-features = ["quantum_frame"] + [[example]] name = "utility" required-features = ["quantum_spd"] diff --git a/wai/wai-quantum/README.md b/wai/wai-quantum/README.md index 9ec24cdb5..f1682622b 100644 --- a/wai/wai-quantum/README.md +++ b/wai/wai-quantum/README.md @@ -69,6 +69,7 @@ never creates information. | `quantum_cal`, `quantum_control`, `quantum_noise` | calibration, filter-function robust control, DD noise spectroscopy, Cycle-Benchmarking noise learning | | `quantum_mitigate` | zero-noise extrapolation, readout M3, classical shadows | | `quantum_qec` | CSS codes (distance verified by exhaustive search, not asserted) — **rotated surface** `[[d²,1,d]]`, toric, bivariate-bicycle, and the **`[[7,1,3]]` colour code**, whose self-duality buys a *transversal* logical Hadamard — one round, no ancillas, no surgery — verified on the simulator, with the surface code shown failing the same move — with **two decoders**: union-find (weighted matching, single-shot, space-time with faulty measurement, **and circuit-level with hook errors**; **imports and exports** Stim-syntax detector error models, so other tools' noise models decode here and ours decode there) and Relay-BP (for qLDPC) | +| `quantum_frame` | **noisy Clifford circuits at QEC scale**: reads the ecosystem's circuit text; samples detection events by bit-packed Pauli-frame simulation (64 shots per word); derives the detector error model by backward sensitivity, refusing non-deterministic detectors; generates the rotated surface-code memory experiment byte for byte as the reference layout; and decodes by **exact minimum-weight perfect matching** (Edmonds' blossom on integer weights). Measures logical error rates and fits the code model resource estimates assume | | `quantum_compile` | Clifford routing + stabilizer-tableau equivalence proof | | `quantum_atom` | neutral-atom register preparation (Hungarian / LSAP) | | `quantum_qir` | **QIR export + ingest** — emit and parse QIR (the QIR Alliance's LLVM-based IR); emitted modules validated with `llvm-as`, and 9 of 10 catalog algorithms re-ingest to a bit-identical statevector | @@ -148,6 +149,32 @@ search limits, to the bit). Re-running it anywhere must reproduce every figure, down to each factory round, and the bytes are the same natively and as a WASI component. +## Error correction at circuit level + +An error-correction experiment is a noisy Clifford circuit, its detectors and +its logical observable. `quantum_frame` samples it, derives its detector error +model, and decodes it exactly: + +```sh +wai-quantum qec -d 5 --p 0.003 --shots 100000 # surface-code memory, per-round logical error +wai-quantum qec --file circuit.txt --shots 100000 # any circuit in the ecosystem's text +wai-quantum qec --dem -d 3 # its detector error model +wai-quantum qec --fit # measure a and p* of P(d) = a·(p/p*)^((d+1)/2) +``` + +The model is checked against an independent implementation: +- the generated circuits are identical, byte for byte, across distances 2–15 + and 1–8 rounds; +- the derived error models match mechanism for mechanism, to 10⁻¹² in + probability, across surface codes in both bases (rotated and unrotated) and + repetition codes; +- decoded logical error rates agree point by point within sampling error. + +`tests/frame_reference.rs` keeps two of those fixtures in the repository. The +fitted code model comes out at `a ≈ 0.023`, `p* ≈ 0.0100`. That is a measured +value for the constants a resource estimate (`quantum_resource`) otherwise takes +on trust, and it can be fed straight into one. + ## Determinism Simulation is **byte-exact**: amplitudes are dyadic fixed-point, so a circuit diff --git a/wai/wai-quantum/examples/qec_frame.rs b/wai/wai-quantum/examples/qec_frame.rs new file mode 100644 index 000000000..50b2658aa --- /dev/null +++ b/wai/wai-quantum/examples/qec_frame.rs @@ -0,0 +1,65 @@ +//! Noisy Clifford circuits from the command line. +//! +//! ```text +//! cargo run --release --example qec_frame --features quantum_frame -- gen 5 5 0.001 +//! cargo run --release --example qec_frame --features quantum_frame -- dem circuit.txt +//! cargo run --release --example qec_frame --features quantum_frame -- sample circuit.txt 100000 7 +//! cargo run --release --example qec_frame --features quantum_frame -- decode circuit.txt events.txt +//! cargo run --release --example qec_frame --features quantum_frame -- rate circuit.txt 100000 7 +//! ``` +use wai_quantum::quantum_frame::{error_model, parse, sample, surface_code_memory}; +use wai_quantum::quantum_match::MatchingGraph; + +fn main() { + let a: Vec = std::env::args().collect(); + let read = |p: &str| std::fs::read_to_string(p).expect("read circuit"); + match a.get(1).map(String::as_str) { + Some("gen") => { + let d: u32 = a[2].parse().unwrap(); + let r: u32 = a[3].parse().unwrap(); + let p: f64 = a[4].parse().unwrap(); + print!("{}", surface_code_memory(d, r, p)); + } + Some("dem") => { + let c = parse(&read(&a[2])).expect("parse"); + print!("{}", error_model(&c).expect("model").to_text()); + } + Some("sample") => { + let c = parse(&read(&a[2])).expect("parse"); + let shots: usize = a[3].parse().unwrap(); + let seed: u64 = a[4].parse().unwrap(); + let s = sample(&c, shots, seed); + for d in 0..s.detectors { + println!("D{d} {}", s.count(d)); + } + let flips = (0..shots).filter(|&i| s.flips(i) & 1 == 1).count(); + println!("L0 {flips}"); + } + // One line per shot of fired detector indices; prints one predicted + // observable mask per line. + Some("decode") => { + let c = parse(&read(&a[2])).expect("parse"); + let g = MatchingGraph::from_model(&error_model(&c).expect("model")).expect("graphlike"); + let mut out = String::new(); + for line in read(&a[3]).lines() { + let fired: Vec = line.split_whitespace().map(|t| t.parse().unwrap()).collect(); + out += &format!("{}\n", g.decode(&fired).expect("decode")); + } + print!("{out}"); + } + // Sample and decode: logical failures out of shots, and the time spent. + Some("rate") => { + let c = parse(&read(&a[2])).expect("parse"); + let shots: usize = a[3].parse().unwrap(); + let seed: u64 = a[4].parse().unwrap(); + let g = MatchingGraph::from_model(&error_model(&c).expect("model")).expect("graphlike"); + let t0 = std::time::Instant::now(); + let s = sample(&c, shots, seed); + let t1 = t0.elapsed(); + let f = g.failures(&s).expect("decode"); + let t2 = t0.elapsed() - t1; + println!("{f} {shots} sample_s={:.3} decode_s={:.3}", t1.as_secs_f64(), t2.as_secs_f64()); + } + _ => eprintln!("usage: qec_frame gen D ROUNDS P | dem FILE | sample FILE SHOTS SEED | decode FILE EVENTS | rate FILE SHOTS SEED"), + } +} diff --git a/wai/wai-quantum/examples/qec_threshold.rs b/wai/wai-quantum/examples/qec_threshold.rs new file mode 100644 index 000000000..91e202609 --- /dev/null +++ b/wai/wai-quantum/examples/qec_threshold.rs @@ -0,0 +1,26 @@ +//! Measure the rotated surface code's logical error rate per round, and fit +//! the code model `P(d) = a·(p/p*)^((d+1)/2)` that resource estimates assume. +//! +//! ```text +//! cargo run --release --example qec_threshold --features quantum_frame -- [SHOTS] [SEED] +//! ``` +use wai_quantum::quantum_match::{fit_code_model, memory_experiment}; + +fn main() { + let a: Vec = std::env::args().collect(); + let shots: usize = a.get(1).map_or(200_000, |s| s.parse().unwrap()); + let seed: u64 = a.get(2).map_or(1, |s| s.parse().unwrap()); + let mut pts = Vec::new(); + println!("{:>3} {:>7} {:>9} {:>9} {:>12}", "d", "p", "shots", "failures", "per round"); + for d in [3u32, 5, 7] { + for p in [0.002, 0.003, 0.004, 0.005] { + // Fewer shots where decoding costs more and failures are plentiful. + let n = if d == 7 { shots / 4 } else { shots }; + let t = memory_experiment(d, d, p, n, seed).expect("decode"); + println!("{:>3} {:>7} {:>9} {:>9} {:>12.4e}", d, p, t.shots, t.failures, t.per_round()); + pts.push(t); + } + } + let f = fit_code_model(&pts).expect("fit"); + println!("fit: a = {:.4}, p* = {:.5}", f.prefactor, f.threshold); +} diff --git a/wai/wai-quantum/src/lib.rs b/wai/wai-quantum/src/lib.rs index 5a868bd7c..e2628802f 100644 --- a/wai/wai-quantum/src/lib.rs +++ b/wai/wai-quantum/src/lib.rs @@ -27,7 +27,8 @@ mod linalg; feature = "quantum_vml", feature = "quantum_phasor", feature = "quantum_kernel", - feature = "quantum_resource" + feature = "quantum_resource", + feature = "quantum_frame" ))] mod repro; @@ -91,6 +92,12 @@ pub mod quantum_resource; #[cfg(feature = "quantum_resource_receipt")] pub mod quantum_resource_receipt; +#[cfg(feature = "quantum_frame")] +pub mod quantum_frame; + +#[cfg(feature = "quantum_frame")] +pub mod quantum_match; + #[cfg(feature = "quantum_info")] pub mod quantum_info; diff --git a/wai/wai-quantum/src/quantum_frame.rs b/wai/wai-quantum/src/quantum_frame.rs new file mode 100644 index 000000000..283504b47 --- /dev/null +++ b/wai/wai-quantum/src/quantum_frame.rs @@ -0,0 +1,1155 @@ +//! 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, + pub measurements: u32, + /// Each detector's measurement records (absolute indices). + pub detectors: Vec>, + /// Each observable's measurement records. + pub observables: Vec>, +} + +#[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 { + 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::().ok()) + .ok_or(FrameError::Parse { line, message: "expected 'REPEAT {'".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 = if args.trim().is_empty() { + vec![] + } else { + args.split(',').map(|a| a.trim().parse::().map_err(|_| err("bad argument"))).collect::>()? + }; + let prob = |i: usize| -> Result { + 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, FrameError> { + targets.split_whitespace().map(|t| t.parse::().map_err(|_| err("bad qubit target"))).collect() + }; + let pairs = || -> Result, 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, FrameError> { + targets + .split_whitespace() + .map(|t| { + let k = t + .strip_prefix("rec[-") + .and_then(|r| r.strip_suffix(']')) + .and_then(|r| r.parse::().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 = data.iter().map(|&(x, y)| index(x, y)).collect(); + data_ix.sort(); + let mut meas_ix: Vec = meas.iter().map(|&((x, y), _)| index(x, y)).collect(); + meas_ix.sort(); + let mut x_ix: Vec = 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::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::>().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 = [(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 = 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 = (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::>().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, + obs: Vec, +} + +impl Detections { + /// The detectors that fired in shot `s`, ascending. + pub fn fired(&self, s: usize) -> Vec { + 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, + /// 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, 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, +} + +/// Symmetric difference of two ascending sets, into `a`. +fn xor_into(a: &mut Vec, 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 { + let nd = c.detectors.len() as u32; + // Symptom ids: detectors 0..nd, observable k as nd + k. + let mut of_meas: Vec> = 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![Vec::new(); nq]; + let mut zs: Vec> = 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, f64, Vec>)> = 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, u64) { + let dets: Vec = 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 = 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, f64, Vec>)>; + +/// An error's X components and Z components, each as the symptom set it flips. +type Halves = (Vec>, Vec>); + +/// 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], zs: &[Vec]) -> 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]| { + 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, 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 = (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}"); + } + } +} diff --git a/wai/wai-quantum/src/quantum_match.rs b/wai/wai-quantum/src/quantum_match.rs new file mode 100644 index 000000000..cde7fdcda --- /dev/null +++ b/wai/wai-quantum/src/quantum_match.rs @@ -0,0 +1,1008 @@ +//! Exact minimum-weight perfect-matching decoding — `wai.quantum.match`. +//! +//! A graphlike detector error model is a graph. Each detector is a node; one +//! more node is the boundary. Each error mechanism is an edge between the +//! detectors it flips, or between one detector and the boundary. An edge +//! weighs `ln((1 − p)/p)`: the less likely the error, the heavier the edge. +//! +//! The fired detectors of a shot must be explained by a set of errors. The +//! most likely set is a minimum-weight perfect matching: +//! - each fired detector pairs with another, or with the boundary; +//! - a pair costs the shortest path between its two ends; +//! - the predicted observable flip is the parity the matched paths carry. +//! +//! This module finds that matching **exactly**. +//! - Shortest paths come from Dijkstra on integer edge weights, so path costs +//! are exact sums. +//! - The matching itself is Edmonds' weighted blossom algorithm, maximising +//! `C − cost` over maximum-cardinality matchings. Fired detectors get +//! boundary twins joined at zero cost, so a perfect matching always exists +//! when the boundary is reachable. +//! +//! Two cuts keep the searches local without giving up exactness. Every +//! detector's path to the boundary is computed once per graph. A pair costing +//! at least both ends' boundary paths is never needed, because the zero-cost +//! twins absorb the swap; so each search stops at that cap. +//! +//! Everything is integer and ordered, so a decode is the same on every machine. +//! [`MatchingGraph::from_model`] refuses a mechanism with a part touching more +//! than two detectors, rather than dropping it: such a model is not graphlike, +//! and a matching decoder would silently mis-score it. +//! +//! # What it measures +//! +//! [`memory_experiment`] samples and decodes the rotated surface-code memory +//! experiment (`d` rounds, every noise knob at `p`). Logical error per round: +//! +//! | d | p = 0.002 | 0.003 | 0.004 | 0.005 | +//! |---|-----------|-------|-------|-------| +//! | 3 | 1.0e-3 | 2.2e-3 | 3.7e-3 | 5.6e-3 | +//! | 5 | 2.1e-4 | 6.6e-4 | 1.5e-3 | 2.8e-3 | +//! | 7 | 4.0e-5 | 1.9e-4 | 6.6e-4 | 1.4e-3 | +//! +//! Distance pays at every rate here, and [`fit_code_model`] turns the grid into +//! the model resource estimates assume, `P(d) = a·(p/p*)^((d+1)/2)` per round: +//! `a ≈ 0.023`, `p* ≈ 0.0100`. Every point was checked against an independent +//! sampler and decoder on the same circuits, and none differs by more than 1.4 +//! standard errors. +//! +//! # Cost +//! +//! The decoder is exact first and fast second. At small distances it decodes +//! in microseconds. At `d = 9` a fired detector's boundary cap still spans +//! much of the lattice, so a shot costs about half a millisecond. A +//! region-growing matcher is the step after this one. Sampling is not the +//! bottleneck: the frame simulator draws over a million `d = 9` shots per second. + +use crate::quantum_frame::{Detections, ErrorModel}; +use crate::repro::ln; +use std::cmp::Reverse; +use std::collections::BinaryHeap; + +/// Integer units per unit of log-likelihood weight. +const SCALE: f64 = 16_384.0; + +#[derive(Clone, Debug, PartialEq, Eq)] +pub enum MatchError { + /// A mechanism part touches this many detectors (more than two). + Hyperedge(usize), + /// No perfect matching: a fired detector cannot reach a partner or the + /// boundary. + Unmatchable, +} + +/// The matching graph of a detector error model. +#[derive(Clone, Debug, PartialEq)] +pub struct MatchingGraph { + detectors: usize, + /// Neighbours of each node (detectors, then the boundary at index + /// `detectors`): `(node, integer weight, observable mask)`. + adj: Vec>, + /// Probability mass of parts that flip observables but no detector: + /// logical errors no decoder can see. + pub undetectable: f64, + /// Each detector's shortest path to the boundary: (cost, observable + /// parity), or `None` if it cannot reach it. + to_boundary: Vec>, +} + +impl MatchingGraph { + /// Build the graph from the model's per-component parts. Parallel edges + /// with the same observable effect are XOR-combined; otherwise the more + /// likely one is kept. + pub fn from_model(m: &ErrorModel) -> Result { + let n = m.detectors as usize; + let boundary = n as u32; + // (u, v) -> (p, obs), u < v. + let mut edges: std::collections::BTreeMap<(u32, u32), (f64, u64)> = std::collections::BTreeMap::new(); + let mut undetectable = 0.0; + for mech in &m.mechanisms { + let p = mech.probability; + for (dets, obs) in &mech.parts { + let key = match dets.as_slice() { + [] => { + if *obs != 0 { + undetectable = undetectable + p - 2.0 * undetectable * p; + } + continue; + } + [a] => (*a, boundary), + [a, b] => (*a.min(b), *a.max(b)), + more => return Err(MatchError::Hyperedge(more.len())), + }; + edges + .entry(key) + .and_modify(|e| { + if e.1 == *obs { + e.0 = e.0 + p - 2.0 * e.0 * p; + } else if p > e.0 { + *e = (p, *obs); + } + }) + .or_insert((p, *obs)); + } + } + let mut adj = vec![Vec::new(); n + 1]; + for ((u, v), (p, obs)) in edges { + let w = if p >= 0.5 { 0 } else { (ln((1.0 - p) / p) * SCALE).round() as i64 }; + adj[u as usize].push((v, w, obs)); + adj[v as usize].push((u, w, obs)); + } + let mut g = MatchingGraph { detectors: n, adj, undetectable, to_boundary: Vec::new() }; + g.to_boundary = g.paths_from(boundary, None)[..n].to_vec(); + Ok(g) + } + + /// Shortest paths from `src`, as (cost, observable parity), never passing + /// through the boundary. With `targets` (each with a cost cap), the search + /// stops once every target is settled or no unsettled target can still be + /// reached within its cap. + fn paths_from(&self, src: u32, targets: Option<&[(u32, i64)]>) -> Vec> { + let n = self.detectors + 1; + let boundary = self.detectors as u32; + let mut best: Vec> = vec![None; n]; + let mut done = vec![false; n]; + let mut open: Vec<(u32, i64)> = targets.map_or(Vec::new(), |t| t.to_vec()); + let mut limit = open.iter().map(|t| t.1).max().unwrap_or(i64::MAX); + let mut heap = BinaryHeap::new(); + best[src as usize] = Some((0, 0)); + heap.push(Reverse((0i64, src, 0u64))); + while let Some(Reverse((c, u, obs))) = heap.pop() { + if done[u as usize] { + continue; + } + if targets.is_some() && c >= limit { + break; + } + done[u as usize] = true; + if let Some(i) = open.iter().position(|t| t.0 == u) { + open.swap_remove(i); + if open.is_empty() { + break; + } + limit = open.iter().map(|t| t.1).max().unwrap(); + } + if u == boundary && u != src { + continue; + } + for &(v, w, o) in &self.adj[u as usize] { + // Per-shot searches leave the boundary to the precomputed table. + if targets.is_some() && v == boundary { + continue; + } + let nc = c + w; + if !done[v as usize] && best[v as usize].is_none_or(|(bc, _)| nc < bc) { + best[v as usize] = Some((nc, obs ^ o)); + heap.push(Reverse((nc, v, obs ^ o))); + } + } + } + best + } + + /// Predict which observables flipped, from the detectors that fired. + pub fn decode(&self, fired: &[u32]) -> Result { + let k = fired.len(); + if k == 0 { + return Ok(0); + } + let reach = |i: usize| self.to_boundary[fired[i] as usize]; + let mut edges: Vec<(usize, usize, i64)> = Vec::new(); + let mut parity: Vec = Vec::new(); + let mut max_cost = 0; + for i in 0..k { + // A pair costing at least both ends' boundary paths never helps: + // the zero-cost boundary twins absorb the swap. So the search from + // `i` looks for each later `j` only within that cap. + let cap = |j: usize| match (reach(i), reach(j)) { + (Some((a, _)), Some((b, _))) => a + b, + _ => i64::MAX, + }; + let targets: Vec<(u32, i64)> = (i + 1..k).map(|j| (fired[j], cap(j))).collect(); + let paths = if targets.is_empty() { Vec::new() } else { self.paths_from(fired[i], Some(&targets)) }; + for j in i + 1..k { + if let Some((c, o)) = paths[fired[j] as usize] { + let dominated = match (reach(i), reach(j)) { + (Some((a, _)), Some((b, _))) => c >= a + b, + _ => false, + }; + if !dominated { + edges.push((i, j, c)); + parity.push(o); + max_cost = max_cost.max(c); + } + } + } + if let Some((c, o)) = reach(i) { + edges.push((i, k + i, c)); + parity.push(o); + max_cost = max_cost.max(c); + } + for j in i + 1..k { + edges.push((k + i, k + j, 0)); + parity.push(0); + } + } + // Maximise C − cost over maximum-cardinality matchings. + let big = max_cost + 1; + let weighted: Vec<(usize, usize, i64)> = edges.iter().map(|&(a, b, c)| (a, b, big - c)).collect(); + let mate = max_weight_matching(2 * k, &weighted, true); + if mate.iter().any(|m| m.is_none()) { + return Err(MatchError::Unmatchable); + } + let mut out = 0; + for (e, &(a, b, _)) in edges.iter().enumerate() { + if mate[a] == Some(b) { + out ^= parity[e]; + } + } + Ok(out) + } + + /// Decode every shot; returns how many predictions missed the actual + /// observable flips. + pub fn failures(&self, shots: &Detections) -> Result { + let mut n = 0; + for s in 0..shots.shots { + if self.decode(&shots.fired(s))? != shots.flips(s) { + n += 1; + } + } + Ok(n) + } +} + +// --------------------------------------------------------------------------- +// Memory experiments and the code model they measure +// --------------------------------------------------------------------------- + +/// The outcome of one surface-code memory experiment. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct MemoryPoint { + pub distance: u32, + pub rounds: u32, + /// Physical error rate (every noise knob). + pub p: f64, + pub shots: u64, + /// Shots whose decoded observable was wrong. + pub failures: u64, +} + +impl MemoryPoint { + /// Logical error per round, `(1 − (1 − 2P)^{1/r}) / 2` for a per-shot + /// failure rate `P` over `r` rounds. + pub fn per_round(&self) -> f64 { + let shot = (self.failures as f64 / self.shots as f64).min(0.499_999); + (1.0 - crate::repro::exp(ln(1.0 - 2.0 * shot) / self.rounds as f64)) / 2.0 + } +} + +/// Sample and decode `shots` runs of the rotated surface-code memory +/// experiment (distance `d`, `rounds` rounds, every noise knob at `p`). +pub fn memory_experiment(d: u32, rounds: u32, p: f64, shots: usize, seed: u64) -> Result { + use crate::quantum_frame::{error_model, parse, sample, surface_code_memory}; + let c = parse(&surface_code_memory(d, rounds, p)).expect("the generated circuit parses"); + let g = MatchingGraph::from_model(&error_model(&c).expect("the generated circuit is deterministic"))?; + let failures = g.failures(&sample(&c, shots, seed))?; + Ok(MemoryPoint { distance: d, rounds, p, shots: shots as u64, failures }) +} + +/// A code model `P(d) = a·(p/p*)^((d+1)/2)` per round, as measured. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct CodeFit { + /// The prefactor `a`. + pub prefactor: f64, + /// The threshold `p*`. + pub threshold: f64, +} + +#[cfg(feature = "quantum_resource")] +impl CodeFit { + /// The gate-based surface-code model of a resource estimate, with this + /// fit's constants in place of assumed ones. + pub fn surface_code(&self) -> crate::quantum_resource::QecModel { + crate::quantum_resource::QecModel { + name: "surface code (measured)".into(), + prefactor: self.prefactor, + threshold: self.threshold, + ..crate::quantum_resource::QecModel::surface_gate() + } + } +} + +/// Fit the code model to measured points by least squares on +/// `ln P − x·ln p = ln a − x·ln p*`, with `x = (d+1)/2`. Each point is +/// weighted by its failure count, the inverse variance of `ln P` to first +/// order. Points without failures carry no information and are skipped. +/// `None` without two distinct distances. +pub fn fit_code_model(points: &[MemoryPoint]) -> Option { + let pts: Vec<(f64, f64, f64)> = points + .iter() + .filter(|q| q.failures > 0) + .map(|q| { + let x = f64::from(q.distance + 1) / 2.0; + (x, ln(q.per_round()) - x * ln(q.p), q.failures as f64) + }) + .collect(); + let wsum: f64 = pts.iter().map(|t| t.2).sum(); + if pts.is_empty() || wsum == 0.0 { + return None; + } + let xm = pts.iter().map(|t| t.2 * t.0).sum::() / wsum; + let ym = pts.iter().map(|t| t.2 * t.1).sum::() / wsum; + let sxx: f64 = pts.iter().map(|t| t.2 * (t.0 - xm) * (t.0 - xm)).sum(); + if sxx == 0.0 { + return None; + } + let sxy: f64 = pts.iter().map(|t| t.2 * (t.0 - xm) * (t.1 - ym)).sum(); + let slope = sxy / sxx; + let intercept = ym - slope * xm; + Some(CodeFit { prefactor: crate::repro::exp(intercept), threshold: crate::repro::exp(-slope) }) +} + +// --------------------------------------------------------------------------- +// Edmonds' weighted blossom algorithm, O(n³) +// --------------------------------------------------------------------------- + +const NONE: usize = usize::MAX; + +/// A maximum-weight matching of an undirected graph on `n` vertices with +/// integer edge weights. With `max_cardinality`, the maximum weight among +/// maximum-cardinality matchings. Returns each vertex's mate. +pub fn max_weight_matching(n: usize, edges: &[(usize, usize, i64)], max_cardinality: bool) -> Vec> { + if edges.is_empty() || n == 0 { + return vec![None; n]; + } + Blossom::new(n, edges).solve(max_cardinality) +} + +struct Blossom<'a> { + n: usize, + edges: &'a [(usize, usize, i64)], + endpoint: Vec, + neighbend: Vec>, + mate: Vec, + label: Vec, + labelend: Vec, + inblossom: Vec, + blossomparent: Vec, + blossomchilds: Vec>, + blossombase: Vec, + blossomendps: Vec>, + bestedge: Vec, + blossombestedges: Vec>>, + unusedblossoms: Vec, + dualvar: Vec, + allowedge: Vec, + queue: Vec, +} + +impl<'a> Blossom<'a> { + fn new(n: usize, edges: &'a [(usize, usize, i64)]) -> Self { + let ne = edges.len(); + let maxweight = edges.iter().map(|e| e.2).max().unwrap_or(0).max(0); + let mut endpoint = Vec::with_capacity(2 * ne); + let mut neighbend = vec![Vec::new(); n]; + for (k, &(i, j, _)) in edges.iter().enumerate() { + endpoint.push(i); + endpoint.push(j); + neighbend[i].push(2 * k + 1); + neighbend[j].push(2 * k); + } + let mut dualvar = vec![maxweight; n]; + dualvar.extend(std::iter::repeat_n(0, n)); + Blossom { + n, + edges, + endpoint, + neighbend, + mate: vec![NONE; n], + label: vec![0; 2 * n], + labelend: vec![NONE; 2 * n], + inblossom: (0..n).collect(), + blossomparent: vec![NONE; 2 * n], + blossomchilds: vec![Vec::new(); 2 * n], + blossombase: (0..n).chain(std::iter::repeat_n(NONE, n)).collect(), + blossomendps: vec![Vec::new(); 2 * n], + bestedge: vec![NONE; 2 * n], + blossombestedges: vec![None; 2 * n], + unusedblossoms: (n..2 * n).collect(), + dualvar, + allowedge: vec![false; ne], + queue: Vec::new(), + } + } + + fn slack(&self, k: usize) -> i64 { + let (i, j, w) = self.edges[k]; + self.dualvar[i] + self.dualvar[j] - 2 * w + } + + fn leaves(&self, b: usize, out: &mut Vec) { + if b < self.n { + out.push(b); + } else { + for &t in &self.blossomchilds[b] { + self.leaves(t, out); + } + } + } + + fn leaves_of(&self, b: usize) -> Vec { + let mut v = Vec::new(); + self.leaves(b, &mut v); + v + } + + fn assign_label(&mut self, w: usize, t: u8, p: usize) { + let b = self.inblossom[w]; + self.label[w] = t; + self.label[b] = t; + self.labelend[w] = p; + self.labelend[b] = p; + self.bestedge[w] = NONE; + self.bestedge[b] = NONE; + if t == 1 { + let l = self.leaves_of(b); + self.queue.extend(l); + } else if t == 2 { + let base = self.blossombase[b]; + let mb = self.mate[base]; + self.assign_label(self.endpoint[mb], 1, mb ^ 1); + } + } + + fn scan_blossom(&mut self, mut v: usize, mut w: usize) -> usize { + let mut path = Vec::new(); + let mut base = NONE; + while v != NONE || w != NONE { + let mut b = self.inblossom[v]; + if self.label[b] & 4 != 0 { + base = self.blossombase[b]; + break; + } + path.push(b); + self.label[b] = 5; + if self.labelend[b] == NONE { + v = NONE; + } else { + v = self.endpoint[self.labelend[b]]; + b = self.inblossom[v]; + v = self.endpoint[self.labelend[b]]; + } + if w != NONE { + core::mem::swap(&mut v, &mut w); + } + } + for b in path { + self.label[b] = 1; + } + base + } + + fn add_blossom(&mut self, base: usize, k: usize) { + let (mut v, mut w, _) = self.edges[k]; + let bb = self.inblossom[base]; + let mut bv = self.inblossom[v]; + let mut bw = self.inblossom[w]; + let b = self.unusedblossoms.pop().expect("a free blossom slot"); + self.blossombase[b] = base; + self.blossomparent[b] = NONE; + self.blossomparent[bb] = b; + let mut path = Vec::new(); + let mut endps = Vec::new(); + while bv != bb { + self.blossomparent[bv] = b; + path.push(bv); + endps.push(self.labelend[bv]); + v = self.endpoint[self.labelend[bv]]; + bv = self.inblossom[v]; + } + path.push(bb); + path.reverse(); + endps.reverse(); + endps.push(2 * k); + while bw != bb { + self.blossomparent[bw] = b; + path.push(bw); + endps.push(self.labelend[bw] ^ 1); + w = self.endpoint[self.labelend[bw]]; + bw = self.inblossom[w]; + } + self.label[b] = 1; + self.labelend[b] = self.labelend[bb]; + self.dualvar[b] = 0; + self.blossomchilds[b] = path.clone(); + self.blossomendps[b] = endps; + for v in self.leaves_of(b) { + if self.label[self.inblossom[v]] == 2 { + self.queue.push(v); + } + self.inblossom[v] = b; + } + let mut bestedgeto = vec![NONE; 2 * self.n]; + for &bv in &path { + let lists: Vec> = match self.blossombestedges[bv].take() { + None => self.leaves_of(bv).iter().map(|&v| self.neighbend[v].iter().map(|p| p / 2).collect()).collect(), + Some(l) => vec![l], + }; + for list in lists { + for k in list { + let (mut i, mut j, _) = self.edges[k]; + if self.inblossom[j] == b { + core::mem::swap(&mut i, &mut j); + } + let _ = i; + let bj = self.inblossom[j]; + if bj != b + && self.label[bj] == 1 + && (bestedgeto[bj] == NONE || self.slack(k) < self.slack(bestedgeto[bj])) + { + bestedgeto[bj] = k; + } + } + } + self.bestedge[bv] = NONE; + } + let best: Vec = bestedgeto.into_iter().filter(|&k| k != NONE).collect(); + self.bestedge[b] = NONE; + for &k in &best { + if self.bestedge[b] == NONE || self.slack(k) < self.slack(self.bestedge[b]) { + self.bestedge[b] = k; + } + } + self.blossombestedges[b] = Some(best); + } + + fn expand_blossom(&mut self, b: usize, endstage: bool) { + let childs = self.blossomchilds[b].clone(); + for &s in &childs { + self.blossomparent[s] = NONE; + if s < self.n { + self.inblossom[s] = s; + } else if endstage && self.dualvar[s] == 0 { + self.expand_blossom(s, endstage); + } else { + for v in self.leaves_of(s) { + self.inblossom[v] = s; + } + } + } + if !endstage && self.label[b] == 2 { + let entrychild = self.inblossom[self.endpoint[self.labelend[b] ^ 1]]; + let len = childs.len() as isize; + let mut j = childs.iter().position(|&c| c == entrychild).unwrap() as isize; + let (jstep, endptrick): (isize, usize) = if j & 1 == 1 { + j -= len; + (1, 0) + } else { + (-1, 1) + }; + let at = |j: isize| -> usize { (j.rem_euclid(len)) as usize }; + let endps = self.blossomendps[b].clone(); + let mut p = self.labelend[b]; + while j != 0 { + let ep = self.endpoint[p ^ 1]; + self.label[ep] = 0; + let q = endps[at(j - endptrick as isize)] ^ endptrick ^ 1; + self.label[self.endpoint[q]] = 0; + self.assign_label(ep, 2, p); + self.allowedge[endps[at(j - endptrick as isize)] / 2] = true; + j += jstep; + p = endps[at(j - endptrick as isize)] ^ endptrick; + self.allowedge[p / 2] = true; + j += jstep; + } + let bv = childs[at(j)]; + let ep = self.endpoint[p ^ 1]; + self.label[ep] = 2; + self.label[bv] = 2; + self.labelend[ep] = p; + self.labelend[bv] = p; + self.bestedge[bv] = NONE; + j += jstep; + while childs[at(j)] != entrychild { + let bv = childs[at(j)]; + if self.label[bv] == 1 { + j += jstep; + continue; + } + let mut found = NONE; + for v in self.leaves_of(bv) { + if self.label[v] != 0 { + found = v; + break; + } + } + if found != NONE { + let v = found; + self.label[v] = 0; + let mb = self.mate[self.blossombase[bv]]; + self.label[self.endpoint[mb]] = 0; + let le = self.labelend[v]; + self.assign_label(v, 2, le); + } + j += jstep; + } + } + self.label[b] = 0; + self.labelend[b] = NONE; + self.blossomchilds[b].clear(); + self.blossomendps[b].clear(); + self.blossombase[b] = NONE; + self.blossombestedges[b] = None; + self.bestedge[b] = NONE; + self.unusedblossoms.push(b); + } + + fn augment_blossom(&mut self, b: usize, v: usize) { + let mut t = v; + while self.blossomparent[t] != b { + t = self.blossomparent[t]; + } + if t >= self.n { + self.augment_blossom(t, v); + } + let childs = self.blossomchilds[b].clone(); + let endps = self.blossomendps[b].clone(); + let len = childs.len() as isize; + let i = childs.iter().position(|&c| c == t).unwrap(); + let mut j = i as isize; + let (jstep, endptrick): (isize, usize) = if i & 1 == 1 { + j -= len; + (1, 0) + } else { + (-1, 1) + }; + let at = |j: isize| -> usize { (j.rem_euclid(len)) as usize }; + while j != 0 { + j += jstep; + let t = childs[at(j)]; + let p = endps[at(j - endptrick as isize)] ^ endptrick; + if t >= self.n { + self.augment_blossom(t, self.endpoint[p]); + } + j += jstep; + let t = childs[at(j)]; + if t >= self.n { + self.augment_blossom(t, self.endpoint[p ^ 1]); + } + self.mate[self.endpoint[p]] = p ^ 1; + self.mate[self.endpoint[p ^ 1]] = p; + } + let mut c = childs[i..].to_vec(); + c.extend_from_slice(&childs[..i]); + let mut e = endps[i..].to_vec(); + e.extend_from_slice(&endps[..i]); + self.blossombase[b] = self.blossombase[c[0]]; + self.blossomchilds[b] = c; + self.blossomendps[b] = e; + } + + fn augment_matching(&mut self, k: usize) { + let (v, w, _) = self.edges[k]; + for (mut s, mut p) in [(v, 2 * k + 1), (w, 2 * k)] { + loop { + let bs = self.inblossom[s]; + if bs >= self.n { + self.augment_blossom(bs, s); + } + self.mate[s] = p; + if self.labelend[bs] == NONE { + break; + } + let t = self.endpoint[self.labelend[bs]]; + let bt = self.inblossom[t]; + s = self.endpoint[self.labelend[bt]]; + let j = self.endpoint[self.labelend[bt] ^ 1]; + if bt >= self.n { + self.augment_blossom(bt, j); + } + self.mate[j] = self.labelend[bt]; + p = self.labelend[bt] ^ 1; + } + } + } + + fn solve(mut self, max_cardinality: bool) -> Vec> { + let n = self.n; + for _ in 0..n { + self.label.iter_mut().for_each(|l| *l = 0); + self.bestedge.iter_mut().for_each(|e| *e = NONE); + for b in n..2 * n { + self.blossombestedges[b] = None; + } + self.allowedge.iter_mut().for_each(|a| *a = false); + self.queue.clear(); + for v in 0..n { + if self.mate[v] == NONE && self.label[self.inblossom[v]] == 0 { + self.assign_label(v, 1, NONE); + } + } + let mut augmented = false; + loop { + while let Some(v) = self.queue.pop() { + if augmented { + break; + } + for idx in 0..self.neighbend[v].len() { + let p = self.neighbend[v][idx]; + let k = p / 2; + let w = self.endpoint[p]; + if self.inblossom[v] == self.inblossom[w] { + continue; + } + let mut kslack = 0; + if !self.allowedge[k] { + kslack = self.slack(k); + if kslack <= 0 { + self.allowedge[k] = true; + } + } + if self.allowedge[k] { + if self.label[self.inblossom[w]] == 0 { + self.assign_label(w, 2, p ^ 1); + } else if self.label[self.inblossom[w]] == 1 { + let base = self.scan_blossom(v, w); + if base != NONE { + self.add_blossom(base, k); + } else { + self.augment_matching(k); + augmented = true; + break; + } + } else if self.label[w] == 0 { + self.label[w] = 2; + self.labelend[w] = p ^ 1; + } + } else if self.label[self.inblossom[w]] == 1 { + let b = self.inblossom[v]; + if self.bestedge[b] == NONE || kslack < self.slack(self.bestedge[b]) { + self.bestedge[b] = k; + } + } else if self.label[w] == 0 + && (self.bestedge[w] == NONE || kslack < self.slack(self.bestedge[w])) + { + self.bestedge[w] = k; + } + } + } + if augmented { + break; + } + // Choose the dual adjustment. + let mut deltatype = 0u8; + let mut delta = 0i64; + let mut deltaedge = NONE; + let mut deltablossom = NONE; + if !max_cardinality { + deltatype = 1; + delta = self.dualvar[..n].iter().copied().min().unwrap(); + } + for v in 0..n { + if self.label[self.inblossom[v]] == 0 && self.bestedge[v] != NONE { + let d = self.slack(self.bestedge[v]); + if deltatype == 0 || d < delta { + delta = d; + deltatype = 2; + deltaedge = self.bestedge[v]; + } + } + } + for b in 0..2 * n { + if self.blossomparent[b] == NONE && self.label[b] == 1 && self.bestedge[b] != NONE { + let kslack = self.slack(self.bestedge[b]); + debug_assert_eq!(kslack % 2, 0); + let d = kslack / 2; + if deltatype == 0 || d < delta { + delta = d; + deltatype = 3; + deltaedge = self.bestedge[b]; + } + } + } + for b in n..2 * n { + if self.blossombase[b] != NONE + && self.blossomparent[b] == NONE + && self.label[b] == 2 + && (deltatype == 0 || self.dualvar[b] < delta) + { + delta = self.dualvar[b]; + deltatype = 4; + deltablossom = b; + } + } + if deltatype == 0 { + deltatype = 1; + delta = self.dualvar[..n].iter().copied().min().unwrap().max(0); + } + for v in 0..n { + match self.label[self.inblossom[v]] { + 1 => self.dualvar[v] -= delta, + 2 => self.dualvar[v] += delta, + _ => {} + } + } + for b in n..2 * n { + if self.blossombase[b] != NONE && self.blossomparent[b] == NONE { + match self.label[b] { + 1 => self.dualvar[b] += delta, + 2 => self.dualvar[b] -= delta, + _ => {} + } + } + } + match deltatype { + 1 => break, + 2 => { + self.allowedge[deltaedge] = true; + let (mut i, mut j, _) = self.edges[deltaedge]; + if self.label[self.inblossom[i]] == 0 { + core::mem::swap(&mut i, &mut j); + } + let _ = j; + self.queue.push(i); + } + 3 => { + self.allowedge[deltaedge] = true; + let (i, _, _) = self.edges[deltaedge]; + self.queue.push(i); + } + _ => self.expand_blossom(deltablossom, false), + } + } + if !augmented { + break; + } + for b in n..2 * n { + if self.blossomparent[b] == NONE + && self.blossombase[b] != NONE + && self.label[b] == 1 + && self.dualvar[b] == 0 + { + self.expand_blossom(b, true); + } + } + } + self.mate.iter().map(|&p| if p == NONE { None } else { Some(self.endpoint[p]) }).collect() + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::quantum_frame::{error_model, parse}; + + /// Minimum-cost perfect matching by exhaustion (None if there is none). + fn brute(n: usize, cost: &[Vec>]) -> Option { + fn go(used: &mut Vec, cost: &[Vec>]) -> Option { + let Some(i) = used.iter().position(|u| !u) else { return Some(0) }; + used[i] = true; + let mut best = None; + for j in i + 1..used.len() { + if used[j] { + continue; + } + if let Some(c) = cost[i][j] { + used[j] = true; + if let Some(rest) = go(used, cost) { + best = Some(best.map_or(c + rest, |b: i64| b.min(c + rest))); + } + used[j] = false; + } + } + used[i] = false; + best + } + go(&mut vec![false; n], cost) + } + + #[test] + #[allow(clippy::needless_range_loop)] + fn blossom_finds_the_minimum_perfect_matching() { + let mut state = 12345u64; + let mut next = |m: u64| { + state = state.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + (state >> 33) % m + }; + for trial in 0..3000 { + let n = 2 * (1 + next(5) as usize); // 2..10 vertices + let density = 40 + next(61); // percent of pairs present + let mut cost = vec![vec![None; n]; n]; + let mut edges = Vec::new(); + for i in 0..n { + for j in i + 1..n { + if next(100) < density { + let c = next(30) as i64; + cost[i][j] = Some(c); + cost[j][i] = Some(c); + edges.push((i, j, c)); + } + } + } + let want = brute(n, &cost); + let big = 1000; + let w: Vec<(usize, usize, i64)> = edges.iter().map(|&(a, b, c)| (a, b, big - c)).collect(); + let mate = max_weight_matching(n, &w, true); + match want { + None => assert!(mate.iter().any(|m| m.is_none()), "trial {trial}: matched with no perfect matching"), + Some(c) => { + assert!(mate.iter().all(|m| m.is_some()), "trial {trial}: not perfect"); + let mut total = 0; + for (i, m) in mate.iter().enumerate() { + let j = m.unwrap(); + assert_eq!(mate[j], Some(i)); + if i < j { + total += cost[i][j].expect("a real edge"); + } + } + assert_eq!(total, c, "trial {trial}"); + } + } + } + } + + #[test] + fn the_repetition_code_decodes_every_correctable_error() { + // Distance 5, one round of perfect checks: every pattern of at most 2 + // data flips is corrected. + let d = 5; + let mut text = String::from("R 0 1 2 3 4 5 6 7 8\n"); + text += "X_ERROR(0.01) 0 1 2 3 4\n"; + for i in 0..d - 1 { + text += &format!("CX {} {} {} {}\n", i, d + i, i + 1, d + i); + } + text += "MR 5 6 7 8\nM 0 1 2 3 4\n"; + for i in 0..d - 1 { + text += &format!("DETECTOR rec[-{}]\n", 9 - i); + } + text += "OBSERVABLE_INCLUDE(0) rec[-1]\n"; + let c = parse(&text).unwrap(); + let g = MatchingGraph::from_model(&error_model(&c).unwrap()).unwrap(); + for mask in 0u32..32 { + if mask.count_ones() > 2 { + continue; + } + let flipped = |q: u32| mask >> q & 1; + let fired: Vec = (0..d - 1).filter(|&i| flipped(i) ^ flipped(i + 1) == 1).collect(); + assert_eq!(g.decode(&fired).unwrap(), flipped(d - 1) as u64, "mask {mask:05b}"); + } + } + + #[test] + fn the_fit_recovers_a_known_model() { + // Exact per-round rates from a = 0.07, p* = 0.009, turned into failure + // counts over a huge number of shots (so rounding is negligible). + let (a, ps): (f64, f64) = (0.07, 0.009); + let mut pts = Vec::new(); + for d in [3u32, 5, 7] { + for p in [0.002, 0.004] { + let per_round: f64 = a * (p / ps).powi(d.div_ceil(2) as i32); + let rounds = d; + let shot = (1.0 - (1.0 - 2.0 * per_round).powi(rounds as i32)) / 2.0; + let shots = 1u64 << 50; + pts.push(MemoryPoint { distance: d, rounds, p, shots, failures: (shot * shots as f64).round() as u64 }); + } + } + let f = fit_code_model(&pts).unwrap(); + assert!((f.prefactor / a - 1.0).abs() < 1e-6, "{f:?}"); + assert!((f.threshold / ps - 1.0).abs() < 1e-6, "{f:?}"); + // One distance cannot separate a from p*. + assert!(fit_code_model(&pts[..2]).is_none()); + #[cfg(feature = "quantum_resource")] + { + let m = f.surface_code(); + assert_eq!((m.prefactor, m.threshold), (f.prefactor, f.threshold)); + assert_eq!(m.tile_qubits(13), crate::quantum_resource::QecModel::surface_gate().tile_qubits(13)); + } + } + + #[test] + fn distance_pays_on_the_surface_code() { + // Below threshold, d = 5 beats d = 3 at the same physical error rate. + let p = 0.003; + let rate = |d: u32| memory_experiment(d, d, p, 20_000, 11).unwrap().per_round(); + let (r3, r5) = (rate(3), rate(5)); + assert!(r5 < r3, "d=3 {r3}, d=5 {r5}"); + } +} diff --git a/wai/wai-quantum/src/repro.rs b/wai/wai-quantum/src/repro.rs index 99902a644..d95b88f82 100644 --- a/wai/wai-quantum/src/repro.rs +++ b/wai/wai-quantum/src/repro.rs @@ -30,28 +30,28 @@ pub use trig::sin_cos; // Natural logarithm // --------------------------------------------------------------------------- -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame"))] const LN2_HI: f64 = 6.931_471_803_691_238_164_90e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame"))] const LN2_LO: f64 = 1.908_214_929_270_587_700_02e-10; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame"))] const LG1: f64 = 6.666_666_666_666_735_130e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame"))] const LG2: f64 = 3.999_999_999_940_941_908e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame"))] const LG3: f64 = 2.857_142_874_366_239_149e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame"))] const LG4: f64 = 2.222_219_843_214_978_396e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame"))] const LG5: f64 = 1.818_357_216_161_805_012e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame"))] const LG6: f64 = 1.531_383_769_920_937_332e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame"))] const LG7: f64 = 1.479_819_860_511_658_591e-1; /// The natural logarithm, within an ulp of the true value: `−∞` at zero, NaN /// below it. -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame"))] pub fn ln(x: f64) -> f64 { let mut bits = x.to_bits(); let mut hx = (bits >> 32) as u32; @@ -91,6 +91,74 @@ pub fn ln(x: f64) -> f64 { s * (hfsq + r) + dk * LN2_LO - hfsq + f + dk * LN2_HI } +// --------------------------------------------------------------------------- +// Exponential +// --------------------------------------------------------------------------- + +#[cfg(feature = "quantum_frame")] +const EXP_P: [f64; 5] = [ + 1.666_666_666_666_660_190_37e-1, + -2.777_777_777_701_559_338_42e-3, + 6.613_756_321_437_934_361_17e-5, + -1.653_390_220_546_525_153_90e-6, + 4.138_136_797_057_238_460_39e-8, +]; + +/// `eˣ`, bit-identical on every IEEE-754 platform: fdlibm's algorithm +/// (reduction by `ln 2` in two parts, a rational approximation on +/// `[−ln2/2, ln2/2]`), using only `+ − × ÷`. Within an ulp of the true value. +#[cfg(feature = "quantum_frame")] +pub fn exp(x: f64) -> f64 { + const LN2_HI_E: f64 = 6.931_471_803_691_238_164_90e-1; + const LN2_LO_E: f64 = 1.908_214_929_270_587_700_02e-10; + const INV_LN2: f64 = 1.442_695_040_888_963_387_00; + if x.is_nan() { + return x; + } + if x > 7.097_827_128_933_839_730_96e2 { + return f64::INFINITY; + } + if x < -7.451_332_191_019_411_084_20e2 { + return 0.0; + } + let ax = x.abs(); + let (hi, lo, k): (f64, f64, i32); + if ax > 0.5 * core::f64::consts::LN_2 { + let kk = if ax < 1.5 * core::f64::consts::LN_2 { + if x < 0.0 { -1 } else { 1 } + } else { + (INV_LN2 * x + if x < 0.0 { -0.5 } else { 0.5 }) as i32 + }; + let t = f64::from(kk); + hi = x - t * LN2_HI_E; + lo = t * LN2_LO_E; + k = kk; + } else if ax < 3.725_290_298_461_914e-9 { + // |x| < 2⁻²⁸: 1 + x to the last bit. + return 1.0 + x; + } else { + hi = x; + lo = 0.0; + k = 0; + } + let r = hi - lo; + let t = r * r; + let c = r - t * (EXP_P[0] + t * (EXP_P[1] + t * (EXP_P[2] + t * (EXP_P[3] + t * EXP_P[4])))); + if k == 0 { + return 1.0 - ((r * c) / (c - 2.0) - r); + } + let y = 1.0 - ((lo - (r * c) / (2.0 - c)) - hi); + // y · 2ᵏ, in two steps where 2ᵏ alone would leave the normal range. + let pow2 = |e: i32| f64::from_bits(((e + 1023) as u64) << 52); + if k > 1023 { + y * pow2(1023) * pow2(k - 1023) + } else if k < -1021 { + y * pow2(k + 1000) * pow2(-1000) + } else { + y * pow2(k) + } +} + // --------------------------------------------------------------------------- // Arctangent // --------------------------------------------------------------------------- @@ -240,6 +308,27 @@ mod tests { #[allow(unused_imports)] use super::*; + #[cfg(feature = "quantum_frame")] + #[test] + fn exp_is_within_an_ulp_of_the_platform() { + let mut s = 0x9e37_79b9_7f4a_7c15u64; + for i in 0..400_000 { + s ^= s << 13; + s ^= s >> 7; + s ^= s << 17; + let u = (s >> 11) as f64 / (1u64 << 53) as f64; + let span = [1e-10, 0.3, 0.7, 2.0, 20.0, 700.0, 740.0][i % 7]; + let x = (2.0 * u - 1.0) * span; + assert!(ulps(exp(x), x.exp()) <= 1, "exp({x:e}) = {:e} vs {:e}", exp(x), x.exp()); + } + assert_eq!(exp(0.0), 1.0); + assert_eq!(exp(f64::INFINITY), f64::INFINITY); + assert_eq!(exp(f64::NEG_INFINITY), 0.0); + assert!(exp(f64::NAN).is_nan()); + assert_eq!(exp(800.0), f64::INFINITY); + assert_eq!(exp(-800.0), 0.0); + } + #[allow(dead_code)] fn ulps(a: f64, b: f64) -> u64 { if a == b { @@ -262,7 +351,7 @@ mod tests { }) } - #[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] + #[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame"))] #[test] fn ln_is_within_an_ulp_of_the_platform() { for x in inputs().filter(|&x| x > 0.0) { diff --git a/wai/wai-quantum/tests/bit_identity.rs b/wai/wai-quantum/tests/bit_identity.rs index 5edf8663a..e28da2d84 100644 --- a/wai/wai-quantum/tests/bit_identity.rs +++ b/wai/wai-quantum/tests/bit_identity.rs @@ -15,7 +15,8 @@ feature = "quantum_pauli", feature = "quantum_spd", feature = "quantum_bptn", - feature = "quantum_resource" + feature = "quantum_resource", + feature = "quantum_frame" ))] use wai_quantum::*; @@ -147,3 +148,27 @@ fn resource_estimates() { assert_eq!(digest(&v), 0xce9d_43b2_d498_f379); } +/// Detection counts of a sampled surface-code run, the failures its decoding +/// leaves, and the code model fitted from a small grid: the sampler's random +/// stream, its logarithm-driven skips, the matching and the fit's `exp`. +#[test] +fn frame_sampling_and_matching() { + use quantum_frame::{parse, sample, surface_code_memory}; + use quantum_match::{fit_code_model, memory_experiment}; + let c = parse(&surface_code_memory(3, 3, 0.01)).unwrap(); + let s = sample(&c, 10_000, 5); + let mut v: Vec = (0..s.detectors).map(|d| s.count(d) as f64).collect(); + let mut pts = Vec::new(); + for d in [3, 5] { + for p in [0.004, 0.008] { + let t = memory_experiment(d, d, p, 4_000, 9).unwrap(); + v.push(t.failures as f64); + v.push(t.per_round()); + pts.push(t); + } + } + let f = fit_code_model(&pts).unwrap(); + v.extend([f.prefactor, f.threshold]); + assert_eq!(digest(&v), 0xcbb3_99f9_6d5c_e8f0); +} + diff --git a/wai/wai-quantum/tests/data/memory_x_d3.circuit b/wai/wai-quantum/tests/data/memory_x_d3.circuit new file mode 100644 index 000000000..ffec8cf7e --- /dev/null +++ b/wai/wai-quantum/tests/data/memory_x_d3.circuit @@ -0,0 +1,89 @@ +QUBIT_COORDS(1, 1) 1 +QUBIT_COORDS(2, 0) 2 +QUBIT_COORDS(3, 1) 3 +QUBIT_COORDS(5, 1) 5 +QUBIT_COORDS(1, 3) 8 +QUBIT_COORDS(2, 2) 9 +QUBIT_COORDS(3, 3) 10 +QUBIT_COORDS(4, 2) 11 +QUBIT_COORDS(5, 3) 12 +QUBIT_COORDS(6, 2) 13 +QUBIT_COORDS(0, 4) 14 +QUBIT_COORDS(1, 5) 15 +QUBIT_COORDS(2, 4) 16 +QUBIT_COORDS(3, 5) 17 +QUBIT_COORDS(4, 4) 18 +QUBIT_COORDS(5, 5) 19 +QUBIT_COORDS(4, 6) 25 +RX 1 3 5 8 10 12 15 17 19 +Z_ERROR(0.002) 1 3 5 8 10 12 15 17 19 +R 2 9 11 13 14 16 18 25 +X_ERROR(0.002) 2 9 11 13 14 16 18 25 +TICK +DEPOLARIZE1(0.002) 1 3 5 8 10 12 15 17 19 +H 2 11 16 25 +DEPOLARIZE1(0.002) 2 11 16 25 +TICK +CX 2 3 16 17 11 12 15 14 10 9 19 18 +DEPOLARIZE2(0.002) 2 3 16 17 11 12 15 14 10 9 19 18 +TICK +CX 2 1 16 15 11 10 8 14 3 9 12 18 +DEPOLARIZE2(0.002) 2 1 16 15 11 10 8 14 3 9 12 18 +TICK +CX 16 10 11 5 25 19 8 9 17 18 12 13 +DEPOLARIZE2(0.002) 16 10 11 5 25 19 8 9 17 18 12 13 +TICK +CX 16 8 11 3 25 17 1 9 10 18 5 13 +DEPOLARIZE2(0.002) 16 8 11 3 25 17 1 9 10 18 5 13 +TICK +H 2 11 16 25 +DEPOLARIZE1(0.002) 2 11 16 25 +TICK +X_ERROR(0.002) 2 9 11 13 14 16 18 25 +MR 2 9 11 13 14 16 18 25 +X_ERROR(0.002) 2 9 11 13 14 16 18 25 +DETECTOR(2, 0, 0) rec[-8] +DETECTOR(2, 4, 0) rec[-3] +DETECTOR(4, 2, 0) rec[-6] +DETECTOR(4, 6, 0) rec[-1] +REPEAT 2 { + TICK + DEPOLARIZE1(0.002) 1 3 5 8 10 12 15 17 19 + H 2 11 16 25 + DEPOLARIZE1(0.002) 2 11 16 25 + TICK + CX 2 3 16 17 11 12 15 14 10 9 19 18 + DEPOLARIZE2(0.002) 2 3 16 17 11 12 15 14 10 9 19 18 + TICK + CX 2 1 16 15 11 10 8 14 3 9 12 18 + DEPOLARIZE2(0.002) 2 1 16 15 11 10 8 14 3 9 12 18 + TICK + CX 16 10 11 5 25 19 8 9 17 18 12 13 + DEPOLARIZE2(0.002) 16 10 11 5 25 19 8 9 17 18 12 13 + TICK + CX 16 8 11 3 25 17 1 9 10 18 5 13 + DEPOLARIZE2(0.002) 16 8 11 3 25 17 1 9 10 18 5 13 + TICK + H 2 11 16 25 + DEPOLARIZE1(0.002) 2 11 16 25 + TICK + X_ERROR(0.002) 2 9 11 13 14 16 18 25 + MR 2 9 11 13 14 16 18 25 + X_ERROR(0.002) 2 9 11 13 14 16 18 25 + SHIFT_COORDS(0, 0, 1) + DETECTOR(2, 0, 0) rec[-8] rec[-16] + DETECTOR(2, 2, 0) rec[-7] rec[-15] + DETECTOR(4, 2, 0) rec[-6] rec[-14] + DETECTOR(6, 2, 0) rec[-5] rec[-13] + DETECTOR(0, 4, 0) rec[-4] rec[-12] + DETECTOR(2, 4, 0) rec[-3] rec[-11] + DETECTOR(4, 4, 0) rec[-2] rec[-10] + DETECTOR(4, 6, 0) rec[-1] rec[-9] +} +Z_ERROR(0.002) 1 3 5 8 10 12 15 17 19 +MX 1 3 5 8 10 12 15 17 19 +DETECTOR(2, 0, 1) rec[-8] rec[-9] rec[-17] +DETECTOR(2, 4, 1) rec[-2] rec[-3] rec[-5] rec[-6] rec[-12] +DETECTOR(4, 2, 1) rec[-4] rec[-5] rec[-7] rec[-8] rec[-15] +DETECTOR(4, 6, 1) rec[-1] rec[-2] rec[-10] +OBSERVABLE_INCLUDE(0) rec[-3] rec[-6] rec[-9] diff --git a/wai/wai-quantum/tests/data/memory_x_d3.dem b/wai/wai-quantum/tests/data/memory_x_d3.dem new file mode 100644 index 000000000..efb96edda --- /dev/null +++ b/wai/wai-quantum/tests/data/memory_x_d3.dem @@ -0,0 +1,221 @@ +error(0.003858066319639164) D0 D2 +error(0.007686363287824861) D0 D4 +error(0.003858066319639164) D0 L0 +error(0.004387566933333155) D1 D2 +error(0.001066666666666516) D1 D2 D5 +error(0.003858066319639164) D1 D3 +error(0.000400266975246834) D1 D4 +error(0.000400266975246834) D1 D4 D5 +error(0.0001334579439585259) D1 D4 D5 D8 +error(0.0001334579439585259) D1 D4 D8 +error(0.000400266975246834) D1 D5 D6 +error(0.0001334579439585259) D1 D5 D6 D8 +error(0.000400266975246834) D1 D5 D8 D9 +error(0.000400266975246834) D1 D5 D8 L0 +error(0.0001334579439585259) D1 D5 D9 +error(0.0001334579439585259) D1 D5 L0 +error(0.000400266975246834) D1 D6 +error(0.0001334579439585259) D1 D6 D8 +error(0.000400266975246834) D1 D8 D9 +error(0.001466079739032824) D1 D8 L0 +error(0.007817769619301565) D1 D9 +error(0.008343044294438688) D1 L0 +error(0.00821142736622363) D10 +error(0.001200017820547319) D10 D11 +error(0.0001334579439585259) D10 D11 D14 +error(0.0002668802658714406) D10 D11 D14 D18 +error(0.000400266975246834) D10 D11 D17 +error(0.0001334579439585259) D10 D11 D17 D18 +error(0.000400266975246834) D10 D11 D18 +error(0.0008002135231907215) D10 D13 +error(0.000400266975246834) D10 D13 D14 D17 +error(0.0001334579439585259) D10 D14 D16 D17 +error(0.0001334579439585259) D10 D14 D17 +error(0.000400266975246834) D10 D14 D17 D18 +error(0.0001334579439585259) D10 D14 D18 +error(0.0002668802658714406) D10 D16 +error(0.0001334579439585259) D10 D16 D17 L0 +error(0.0001334579439585259) D10 D17 L0 +error(0.004785824618421568) D10 D18 +error(0.001599324140787608) D11 +error(0.0002668802658714406) D11 D14 +error(0.0001334579439585259) D11 D14 D18 +error(0.0006669335944047278) D11 D17 +error(0.000400266975246834) D11 D17 D18 +error(0.000533618081590258) D11 D18 +error(0.000533618081590258) D11 D18 D19 +error(0.007160387030215617) D11 D19 +error(0.001200017820547319) D12 D13 D14 +error(0.0001334579439585259) D12 D13 D16 L0 +error(0.000533618081590258) D12 D13 D20 +error(0.001200017820547319) D12 D13 L0 +error(0.001200017820547319) D12 D14 +error(0.007160387030215617) D12 D20 +error(0.001599324140787608) D12 L0 +error(0.008211427366223628) D13 +error(0.0001334579439585259) D13 D14 +error(0.0001334579439585259) D13 D14 D16 D17 +error(0.000533618081590258) D13 D14 D17 +error(0.001200017820547319) D13 D14 D17 D18 +error(0.000533618081590258) D13 D14 D20 +error(0.001466257564268764) D13 D16 +error(0.001200017820547319) D13 D16 D17 L0 +error(0.000533618081590258) D13 D17 D20 +error(0.000533618081590258) D13 D17 L0 +error(0.001200017820547319) D13 D18 +error(0.001066666666666516) D13 D20 L0 +error(0.002928215306047978) D14 +error(0.001200017820547319) D14 D15 +error(0.0001334579439585259) D14 D15 D17 +error(0.0001334579439585259) D14 D15 D17 D18 +error(0.001200017820547319) D14 D15 D18 +error(0.000400266975246834) D14 D15 D18 D22 +error(0.000400266975246834) D14 D15 D22 +error(0.001865350750555233) D14 D17 +error(0.000533618081590258) D14 D17 D18 +error(0.0001334579439585259) D14 D18 D22 +error(0.001066666666666516) D14 D20 +error(0.007817769619301565) D14 D22 +error(0.006501598321819228) D15 +error(0.0001334579439585259) D15 D17 D18 D22 +error(0.0001334579439585259) D15 D17 D22 +error(0.001466257564268764) D15 D18 +error(0.0001334579439585259) D15 D18 D19 +error(0.0001334579439585259) D15 D18 D19 D22 +error(0.000400266975246834) D15 D18 D22 +error(0.0001334579439585259) D15 D19 +error(0.0001334579439585259) D15 D19 D22 +error(0.001466079739032824) D15 D22 +error(0.005842457790162801) D16 +error(0.001599324140787608) D16 D17 L0 +error(0.001200017820547319) D17 D18 D19 +error(0.0008002135231907215) D17 D18 D21 +error(0.000400266975246834) D17 D18 D22 +error(0.0002668802658714406) D17 D18 L0 +error(0.001200017820547319) D17 D19 +error(0.000533618081590258) D17 D20 +error(0.007949140876340624) D17 D21 +error(0.000400266975246834) D17 D22 +error(0.003591243050134198) D17 L0 +error(0.008867408007091262) D18 +error(0.001599324140787608) D18 D19 +error(0.000533618081590258) D18 D19 D21 +error(0.000400266975246834) D18 D19 D22 +error(0.001066666666666516) D18 D21 D22 +error(0.0002668802658714406) D18 D21 L0 +error(0.0001334579439585259) D18 D22 +error(0.00213123536581817) D19 +error(0.001066666666666516) D19 D21 +error(0.000400266975246834) D19 D22 +error(0.007686363287824861) D19 D23 +error(0.007423445364098325) D2 +error(0.001066666666666516) D2 D4 +error(0.000533618081590258) D2 D4 D5 +error(0.0002668802658714406) D2 D5 +error(0.0008002135231907215) D2 D5 D6 +error(0.007949140876340624) D2 D6 +error(0.002531483609263898) D20 D22 +error(0.003062399999999851) D20 L0 +error(0.003062399999999851) D21 D22 +error(0.002531483609263898) D21 D23 +error(0.004785824618421568) D21 L0 +error(0.005710347764487468) D22 +error(0.002531483609263898) D23 +error(0.004387566933333155) D3 +error(0.001066666666666516) D3 D10 +error(0.000533618081590258) D3 D10 D11 +error(0.007160387030215617) D3 D11 +error(0.000533618081590258) D3 D6 +error(0.000533618081590258) D3 D6 D10 +error(0.001066666666666516) D3 D9 +error(0.000533618081590258) D3 D9 D10 +error(0.007160387030215617) D4 D12 +error(0.000533618081590258) D4 D5 D12 +error(0.001200017820547319) D4 D5 D6 +error(0.0001334579439585259) D4 D5 D8 L0 +error(0.001599324140787608) D4 D5 L0 +error(0.001200017820547319) D4 D6 +error(0.0001334579439585259) D4 D8 L0 +error(0.00213123536581817) D4 L0 +error(0.008867408007091262) D5 +error(0.001200017820547319) D5 D10 +error(0.000400266975246834) D5 D12 D13 L0 +error(0.0001334579439585259) D5 D12 D16 L0 +error(0.000533618081590258) D5 D12 L0 +error(0.004785824618421568) D5 D13 +error(0.000400266975246834) D5 D16 +error(0.0002668802658714406) D5 D6 +error(0.000400266975246834) D5 D6 D12 +error(0.0001334579439585259) D5 D6 D12 D13 +error(0.0001334579439585259) D5 D6 D8 D9 +error(0.000533618081590258) D5 D6 D9 +error(0.001200017820547319) D5 D6 D9 D10 +error(0.000400266975246834) D5 D6 D9 D13 +error(0.001466257564268764) D5 D8 +error(0.001200017820547319) D5 D8 D9 L0 +error(0.0001334579439585259) D5 D9 D12 +error(0.0002668802658714406) D5 D9 D12 D13 +error(0.0001334579439585259) D5 D9 D12 D16 +error(0.0001334579439585259) D5 D9 D13 L0 +error(0.000400266975246834) D5 D9 D16 L0 +error(0.0035912430501342) D6 +error(0.000533618081590258) D6 D10 +error(0.0001334579439585259) D6 D10 D13 D14 +error(0.0006669335944047278) D6 D12 +error(0.000400266975246834) D6 D12 D13 +error(0.0001334579439585259) D6 D13 +error(0.000400266975246834) D6 D13 D14 +error(0.007423445364098326) D6 D14 +error(0.001599324140787608) D6 D7 +error(0.001200017820547319) D6 D7 D10 +error(0.000400266975246834) D6 D7 D10 D14 +error(0.0001334579439585259) D6 D7 D13 +error(0.0002668802658714406) D6 D7 D13 D14 +error(0.0001334579439585259) D6 D7 D14 +error(0.0001334579439585259) D6 D7 D9 D10 +error(0.0001334579439585259) D6 D7 D9 D13 +error(0.0001334579439585259) D6 D8 D9 +error(0.001865350750555233) D6 D9 +error(0.000533618081590258) D6 D9 D10 +error(0.000400266975246834) D6 D9 D10 D13 +error(0.0001334579439585259) D6 D9 D13 +error(0.005842457790162799) D7 +error(0.001466257564268764) D7 D10 +error(0.0001334579439585259) D7 D10 D11 +error(0.0001334579439585259) D7 D10 D11 D14 +error(0.000400266975246834) D7 D10 D14 +error(0.0001334579439585259) D7 D11 D14 D18 +error(0.0001334579439585259) D7 D11 D18 +error(0.0002668802658714406) D7 D13 +error(0.0001334579439585259) D7 D13 D14 +error(0.000400266975246834) D7 D14 +error(0.000533618081590258) D7 D14 D15 +error(0.000400266975246834) D7 D14 D18 +error(0.004521357674826842) D7 D15 +error(0.000400266975246834) D7 D18 +error(0.0001334579439585259) D7 D9 D10 D14 +error(0.0001334579439585259) D7 D9 D13 D14 +error(0.006501598321819227) D8 +error(0.004521357674826842) D8 D16 +error(0.000533618081590258) D8 D9 D16 L0 +error(0.001200017820547319) D8 D9 L0 +error(0.001200017820547319) D9 D10 D11 +error(0.0002668802658714406) D9 D10 D13 D14 +error(0.0001334579439585259) D9 D10 D13 D17 +error(0.0001334579439585259) D9 D10 D14 D16 +error(0.0002668802658714406) D9 D10 D16 D17 +error(0.0001334579439585259) D9 D10 D16 L0 +error(0.000400266975246834) D9 D10 D17 +error(0.0001334579439585259) D9 D10 L0 +error(0.001200017820547319) D9 D11 +error(0.0002668802658714406) D9 D12 +error(0.0001334579439585259) D9 D12 D13 +error(0.0001334579439585259) D9 D12 D13 D16 +error(0.0001334579439585259) D9 D13 D14 D16 +error(0.000400266975246834) D9 D13 D16 D17 +error(0.000400266975246834) D9 D13 D16 L0 +error(0.0002668802658714406) D9 D14 +error(0.0001334579439585259) D9 D16 D17 +error(0.000400266975246834) D9 D16 L0 +error(0.007423445364098326) D9 D17 +error(0.002928215306047977) D9 L0 diff --git a/wai/wai-quantum/tests/data/memory_z_d5.circuit b/wai/wai-quantum/tests/data/memory_z_d5.circuit new file mode 100644 index 000000000..dc4bd5570 --- /dev/null +++ b/wai/wai-quantum/tests/data/memory_z_d5.circuit @@ -0,0 +1,153 @@ +QUBIT_COORDS(1, 1) 1 +QUBIT_COORDS(2, 0) 2 +QUBIT_COORDS(3, 1) 3 +QUBIT_COORDS(5, 1) 5 +QUBIT_COORDS(6, 0) 6 +QUBIT_COORDS(7, 1) 7 +QUBIT_COORDS(9, 1) 9 +QUBIT_COORDS(1, 3) 12 +QUBIT_COORDS(2, 2) 13 +QUBIT_COORDS(3, 3) 14 +QUBIT_COORDS(4, 2) 15 +QUBIT_COORDS(5, 3) 16 +QUBIT_COORDS(6, 2) 17 +QUBIT_COORDS(7, 3) 18 +QUBIT_COORDS(8, 2) 19 +QUBIT_COORDS(9, 3) 20 +QUBIT_COORDS(10, 2) 21 +QUBIT_COORDS(0, 4) 22 +QUBIT_COORDS(1, 5) 23 +QUBIT_COORDS(2, 4) 24 +QUBIT_COORDS(3, 5) 25 +QUBIT_COORDS(4, 4) 26 +QUBIT_COORDS(5, 5) 27 +QUBIT_COORDS(6, 4) 28 +QUBIT_COORDS(7, 5) 29 +QUBIT_COORDS(8, 4) 30 +QUBIT_COORDS(9, 5) 31 +QUBIT_COORDS(1, 7) 34 +QUBIT_COORDS(2, 6) 35 +QUBIT_COORDS(3, 7) 36 +QUBIT_COORDS(4, 6) 37 +QUBIT_COORDS(5, 7) 38 +QUBIT_COORDS(6, 6) 39 +QUBIT_COORDS(7, 7) 40 +QUBIT_COORDS(8, 6) 41 +QUBIT_COORDS(9, 7) 42 +QUBIT_COORDS(10, 6) 43 +QUBIT_COORDS(0, 8) 44 +QUBIT_COORDS(1, 9) 45 +QUBIT_COORDS(2, 8) 46 +QUBIT_COORDS(3, 9) 47 +QUBIT_COORDS(4, 8) 48 +QUBIT_COORDS(5, 9) 49 +QUBIT_COORDS(6, 8) 50 +QUBIT_COORDS(7, 9) 51 +QUBIT_COORDS(8, 8) 52 +QUBIT_COORDS(9, 9) 53 +QUBIT_COORDS(4, 10) 59 +QUBIT_COORDS(8, 10) 63 +R 1 3 5 7 9 12 14 16 18 20 23 25 27 29 31 34 36 38 40 42 45 47 49 51 53 +X_ERROR(0.001) 1 3 5 7 9 12 14 16 18 20 23 25 27 29 31 34 36 38 40 42 45 47 49 51 53 +R 2 6 13 15 17 19 21 22 24 26 28 30 35 37 39 41 43 44 46 48 50 52 59 63 +X_ERROR(0.001) 2 6 13 15 17 19 21 22 24 26 28 30 35 37 39 41 43 44 46 48 50 52 59 63 +TICK +DEPOLARIZE1(0.001) 1 3 5 7 9 12 14 16 18 20 23 25 27 29 31 34 36 38 40 42 45 47 49 51 53 +H 2 6 15 19 24 28 37 41 46 50 59 63 +DEPOLARIZE1(0.001) 2 6 15 19 24 28 37 41 46 50 59 63 +TICK +CX 2 3 24 25 46 47 15 16 37 38 6 7 28 29 50 51 19 20 41 42 23 22 45 44 14 13 36 35 27 26 49 48 18 17 40 39 31 30 53 52 +DEPOLARIZE2(0.001) 2 3 24 25 46 47 15 16 37 38 6 7 28 29 50 51 19 20 41 42 23 22 45 44 14 13 36 35 27 26 49 48 18 17 40 39 31 30 53 52 +TICK +CX 2 1 24 23 46 45 15 14 37 36 6 5 28 27 50 49 19 18 41 40 12 22 34 44 3 13 25 35 16 26 38 48 7 17 29 39 20 30 42 52 +DEPOLARIZE2(0.001) 2 1 24 23 46 45 15 14 37 36 6 5 28 27 50 49 19 18 41 40 12 22 34 44 3 13 25 35 16 26 38 48 7 17 29 39 20 30 42 52 +TICK +CX 24 14 46 36 15 5 37 27 59 49 28 18 50 40 19 9 41 31 63 53 12 13 34 35 25 26 47 48 16 17 38 39 29 30 51 52 20 21 42 43 +DEPOLARIZE2(0.001) 24 14 46 36 15 5 37 27 59 49 28 18 50 40 19 9 41 31 63 53 12 13 34 35 25 26 47 48 16 17 38 39 29 30 51 52 20 21 42 43 +TICK +CX 24 12 46 34 15 3 37 25 59 47 28 16 50 38 19 7 41 29 63 51 1 13 23 35 14 26 36 48 5 17 27 39 18 30 40 52 9 21 31 43 +DEPOLARIZE2(0.001) 24 12 46 34 15 3 37 25 59 47 28 16 50 38 19 7 41 29 63 51 1 13 23 35 14 26 36 48 5 17 27 39 18 30 40 52 9 21 31 43 +TICK +H 2 6 15 19 24 28 37 41 46 50 59 63 +DEPOLARIZE1(0.001) 2 6 15 19 24 28 37 41 46 50 59 63 +TICK +X_ERROR(0.001) 2 6 13 15 17 19 21 22 24 26 28 30 35 37 39 41 43 44 46 48 50 52 59 63 +MR 2 6 13 15 17 19 21 22 24 26 28 30 35 37 39 41 43 44 46 48 50 52 59 63 +X_ERROR(0.001) 2 6 13 15 17 19 21 22 24 26 28 30 35 37 39 41 43 44 46 48 50 52 59 63 +DETECTOR(0, 4, 0) rec[-17] +DETECTOR(0, 8, 0) rec[-7] +DETECTOR(2, 2, 0) rec[-22] +DETECTOR(2, 6, 0) rec[-12] +DETECTOR(4, 4, 0) rec[-15] +DETECTOR(4, 8, 0) rec[-5] +DETECTOR(6, 2, 0) rec[-20] +DETECTOR(6, 6, 0) rec[-10] +DETECTOR(8, 4, 0) rec[-13] +DETECTOR(8, 8, 0) rec[-3] +DETECTOR(10, 2, 0) rec[-18] +DETECTOR(10, 6, 0) rec[-8] +REPEAT 3 { + TICK + DEPOLARIZE1(0.001) 1 3 5 7 9 12 14 16 18 20 23 25 27 29 31 34 36 38 40 42 45 47 49 51 53 + H 2 6 15 19 24 28 37 41 46 50 59 63 + DEPOLARIZE1(0.001) 2 6 15 19 24 28 37 41 46 50 59 63 + TICK + CX 2 3 24 25 46 47 15 16 37 38 6 7 28 29 50 51 19 20 41 42 23 22 45 44 14 13 36 35 27 26 49 48 18 17 40 39 31 30 53 52 + DEPOLARIZE2(0.001) 2 3 24 25 46 47 15 16 37 38 6 7 28 29 50 51 19 20 41 42 23 22 45 44 14 13 36 35 27 26 49 48 18 17 40 39 31 30 53 52 + TICK + CX 2 1 24 23 46 45 15 14 37 36 6 5 28 27 50 49 19 18 41 40 12 22 34 44 3 13 25 35 16 26 38 48 7 17 29 39 20 30 42 52 + DEPOLARIZE2(0.001) 2 1 24 23 46 45 15 14 37 36 6 5 28 27 50 49 19 18 41 40 12 22 34 44 3 13 25 35 16 26 38 48 7 17 29 39 20 30 42 52 + TICK + CX 24 14 46 36 15 5 37 27 59 49 28 18 50 40 19 9 41 31 63 53 12 13 34 35 25 26 47 48 16 17 38 39 29 30 51 52 20 21 42 43 + DEPOLARIZE2(0.001) 24 14 46 36 15 5 37 27 59 49 28 18 50 40 19 9 41 31 63 53 12 13 34 35 25 26 47 48 16 17 38 39 29 30 51 52 20 21 42 43 + TICK + CX 24 12 46 34 15 3 37 25 59 47 28 16 50 38 19 7 41 29 63 51 1 13 23 35 14 26 36 48 5 17 27 39 18 30 40 52 9 21 31 43 + DEPOLARIZE2(0.001) 24 12 46 34 15 3 37 25 59 47 28 16 50 38 19 7 41 29 63 51 1 13 23 35 14 26 36 48 5 17 27 39 18 30 40 52 9 21 31 43 + TICK + H 2 6 15 19 24 28 37 41 46 50 59 63 + DEPOLARIZE1(0.001) 2 6 15 19 24 28 37 41 46 50 59 63 + TICK + X_ERROR(0.001) 2 6 13 15 17 19 21 22 24 26 28 30 35 37 39 41 43 44 46 48 50 52 59 63 + MR 2 6 13 15 17 19 21 22 24 26 28 30 35 37 39 41 43 44 46 48 50 52 59 63 + X_ERROR(0.001) 2 6 13 15 17 19 21 22 24 26 28 30 35 37 39 41 43 44 46 48 50 52 59 63 + SHIFT_COORDS(0, 0, 1) + DETECTOR(2, 0, 0) rec[-24] rec[-48] + DETECTOR(6, 0, 0) rec[-23] rec[-47] + DETECTOR(2, 2, 0) rec[-22] rec[-46] + DETECTOR(4, 2, 0) rec[-21] rec[-45] + DETECTOR(6, 2, 0) rec[-20] rec[-44] + DETECTOR(8, 2, 0) rec[-19] rec[-43] + DETECTOR(10, 2, 0) rec[-18] rec[-42] + DETECTOR(0, 4, 0) rec[-17] rec[-41] + DETECTOR(2, 4, 0) rec[-16] rec[-40] + DETECTOR(4, 4, 0) rec[-15] rec[-39] + DETECTOR(6, 4, 0) rec[-14] rec[-38] + DETECTOR(8, 4, 0) rec[-13] rec[-37] + DETECTOR(2, 6, 0) rec[-12] rec[-36] + DETECTOR(4, 6, 0) rec[-11] rec[-35] + DETECTOR(6, 6, 0) rec[-10] rec[-34] + DETECTOR(8, 6, 0) rec[-9] rec[-33] + DETECTOR(10, 6, 0) rec[-8] rec[-32] + DETECTOR(0, 8, 0) rec[-7] rec[-31] + DETECTOR(2, 8, 0) rec[-6] rec[-30] + DETECTOR(4, 8, 0) rec[-5] rec[-29] + DETECTOR(6, 8, 0) rec[-4] rec[-28] + DETECTOR(8, 8, 0) rec[-3] rec[-27] + DETECTOR(4, 10, 0) rec[-2] rec[-26] + DETECTOR(8, 10, 0) rec[-1] rec[-25] +} +X_ERROR(0.001) 1 3 5 7 9 12 14 16 18 20 23 25 27 29 31 34 36 38 40 42 45 47 49 51 53 +M 1 3 5 7 9 12 14 16 18 20 23 25 27 29 31 34 36 38 40 42 45 47 49 51 53 +DETECTOR(0, 4, 1) rec[-15] rec[-20] rec[-42] +DETECTOR(0, 8, 1) rec[-5] rec[-10] rec[-32] +DETECTOR(2, 2, 1) rec[-19] rec[-20] rec[-24] rec[-25] rec[-47] +DETECTOR(2, 6, 1) rec[-9] rec[-10] rec[-14] rec[-15] rec[-37] +DETECTOR(4, 4, 1) rec[-13] rec[-14] rec[-18] rec[-19] rec[-40] +DETECTOR(4, 8, 1) rec[-3] rec[-4] rec[-8] rec[-9] rec[-30] +DETECTOR(6, 2, 1) rec[-17] rec[-18] rec[-22] rec[-23] rec[-45] +DETECTOR(6, 6, 1) rec[-7] rec[-8] rec[-12] rec[-13] rec[-35] +DETECTOR(8, 4, 1) rec[-11] rec[-12] rec[-16] rec[-17] rec[-38] +DETECTOR(8, 8, 1) rec[-1] rec[-2] rec[-6] rec[-7] rec[-28] +DETECTOR(10, 2, 1) rec[-16] rec[-21] rec[-43] +DETECTOR(10, 6, 1) rec[-6] rec[-11] rec[-33] +OBSERVABLE_INCLUDE(0) rec[-21] rec[-22] rec[-23] rec[-24] rec[-25] diff --git a/wai/wai-quantum/tests/frame_reference.rs b/wai/wai-quantum/tests/frame_reference.rs new file mode 100644 index 000000000..7dc797807 --- /dev/null +++ b/wai/wai-quantum/tests/frame_reference.rs @@ -0,0 +1,46 @@ +//! Cross-tool fixtures for noisy Clifford circuits. +//! +//! `data/memory_z_d5.circuit` and `data/memory_x_d3.circuit` were written by +//! an independent circuit generator. `data/memory_x_d3.dem` is that tool's +//! detector error model of the second circuit, flattened, with identical +//! mechanisms XOR-combined. Matching them checks three things against an +//! implementation that shares no code with this one: +//! - the generator, byte for byte; +//! - the parser, including the X-basis resets and measurements our generator +//! never emits; +//! - the backward-sensitivity model, mechanism for mechanism. +#![cfg(feature = "quantum_frame")] + +use std::collections::BTreeMap; +use wai_quantum::quantum_frame::{error_model, parse, surface_code_memory}; + +#[test] +fn the_generator_matches_the_reference_text() { + let want = include_str!("data/memory_z_d5.circuit"); + assert_eq!(surface_code_memory(5, 4, 0.001), want); +} + +/// `"D0 D4 L0"` → probability. +fn model_map(text: &str) -> BTreeMap { + text.lines() + .filter_map(|l| { + let rest = l.trim().strip_prefix("error(")?; + let (p, syms) = rest.split_once(')')?; + let mut s: Vec<&str> = syms.split_whitespace().collect(); + s.sort_by_key(|t| (t.as_bytes()[0], t[1..].parse::().unwrap())); + Some((s.join(" "), p.parse().unwrap())) + }) + .collect() +} + +#[test] +fn the_model_matches_the_reference_model() { + let c = parse(include_str!("data/memory_x_d3.circuit")).unwrap(); + let ours = model_map(&error_model(&c).unwrap().to_text()); + let theirs = model_map(include_str!("data/memory_x_d3.dem")); + assert_eq!(ours.len(), theirs.len()); + for (k, p) in &theirs { + let q = ours.get(k).unwrap_or_else(|| panic!("missing mechanism {k}")); + assert!((q / p - 1.0).abs() < 1e-9, "{k}: {q} vs {p}"); + } +} diff --git a/wai/wai-rs/Cargo.lock b/wai/wai-rs/Cargo.lock index 283b86ca3..44da8e1a9 100644 --- a/wai/wai-rs/Cargo.lock +++ b/wai/wai-rs/Cargo.lock @@ -3900,7 +3900,7 @@ dependencies = [ [[package]] name = "wai-quantum" -version = "0.3.24" +version = "0.3.25" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-rs/Cargo.toml b/wai/wai-rs/Cargo.toml index 6fe930a16..4bf6b742f 100644 --- a/wai/wai-rs/Cargo.toml +++ b/wai/wai-rs/Cargo.toml @@ -348,6 +348,9 @@ quantum_bptn = ["quantum_spd", "wai-quantum/quantum_bptn"] # and its signed receipts. quantum_resource = ["dep:wai-quantum", "wai-quantum/quantum_resource"] quantum_resource_receipt = ["quantum_resource", "wai-quantum/quantum_resource_receipt"] +# Noisy Clifford circuits at QEC scale: frame sampling, detector error models, +# exact matching, measured code models. +quantum_frame = ["dep:wai-quantum", "wai-quantum/quantum_frame"] quantum_qir = ["quantum_toolchain", "wai-quantum/quantum_qir"] # Variational quantum machine learning (extensions/quantum-sim § VML, # `wai.quantum.vml`): the interference-based ML the thermodynamic-computing pivot @@ -476,6 +479,7 @@ required-features = [ "quantum_spd_receipt", "quantum_bptn", "quantum_resource_receipt", + "quantum_frame", ] [[bin]] diff --git a/wai/wai-rs/src/bin/wai_quantum.rs b/wai/wai-rs/src/bin/wai_quantum.rs index 059a29314..bf2935145 100644 --- a/wai/wai-rs/src/bin/wai_quantum.rs +++ b/wai/wai-rs/src/bin/wai_quantum.rs @@ -47,6 +47,16 @@ //! wai-quantum estimate factoring --qubit ns-e4 --receipt | wai-quantum verify //! ``` //! +//! Error correction itself, at circuit level: sample a noisy Clifford circuit, +//! derive its detector error model, decode it exactly, and measure the code +//! model that estimates otherwise take on trust: +//! +//! ```text +//! wai-quantum qec -d 5 --p 0.003 --shots 100000 +//! wai-quantum qec --file circuit.txt +//! wai-quantum qec --fit # then: estimate --prefactor A --threshold P +//! ``` +//! //! The `channel` command prints the Holevo bound before and after, which can only //! fall: noise never creates information. Two results there are easy to miss — //! fully damping `|1⟩` yields a *pure* output that still carries nothing (both @@ -145,7 +155,12 @@ fn usage() -> ! { wai-quantum estimate [factoring|chemistry] [--qubit NAME|all] [--code surface|surface-meas|floquet]\n \ [--budget E] [--slowdown N] [--frontier] [--json] [--receipt]\n \ wai-quantum estimate --qubits N [--ccz N] [--t N] [--rotations N] [--depth N] [--measurements N] ...\n \ - Qubits: us-e3, us-e4, ns-e3, ns-e4, meas-e4, meas-e6.\n\n\ + Qubits: us-e3, us-e4, ns-e3, ns-e4, meas-e4, meas-e6. A measured code: --prefactor A --threshold P.\n\n\ + ERROR CORRECTION AT CIRCUIT LEVEL:\n \ + wai-quantum qec [-d D] [--rounds R] [--p P] [--shots N] [--seed S] [--json]\n \ + wai-quantum qec --file CIRCUIT [--shots N] (the ecosystem's circuit text)\n \ + wai-quantum qec [--circuit | --dem] [-d D | --file CIRCUIT]\n \ + wai-quantum qec --fit [--shots N] (measure a and p* of the code model)\n\n\ Channels: identity, bitflip, dephasing, depolarizing, damping.\n\ Run `wai-quantum list` for the algorithm catalog." ); @@ -560,6 +575,20 @@ fn main() { "floquet" => QecModel::floquet(), other => die(&format!("unknown code '{other}' (surface, surface-meas, floquet)")), }); + // A code model with given constants, e.g. as `qec --fit` measured them. + let a = opt(args, "--prefactor", "").map(|v| v.parse::().unwrap_or_else(|_| die("bad --prefactor"))); + let t = opt(args, "--threshold", "").map(|v| v.parse::().unwrap_or_else(|_| die("bad --threshold"))); + let code = if a.is_some() || t.is_some() { + let base = code.unwrap_or_else(QecModel::surface_gate); + Some(QecModel { + name: format!("{} (given)", base.name), + prefactor: a.unwrap_or(base.prefactor), + threshold: t.unwrap_or(base.threshold), + ..base + }) + } else { + code + }; let which = opt(args, "--qubit", "").unwrap_or_else(|| { if flag(args, "--receipt") || flag(args, "--frontier") { "ns-e4".into() } else { "all".into() } }); @@ -631,13 +660,13 @@ fn main() { let isa = qr::layout(&counts, budget); println!("program {} tiles, {} logical steps, {} T states (budget {budget})", isa.tiles, isa.min_steps, isa.t_states); println!( - "{:<12} {:<13} {:>3} {:>9} {:>8} {:>14} {}", + "{:<12} {:<22} {:>3} {:>9} {:>8} {:>14} {}", "qubit", "code", "d", "factories", "share", "physical", "runtime" ); for (q, e) in &estimates { match e { Ok(e) => println!( - "{:<12} {:<13} {:>3} {:>9} {:>7.1}% {:>14} {}", + "{:<12} {:<22} {:>3} {:>9} {:>7.1}% {:>14} {}", e.qubit, e.code, e.distance, @@ -651,6 +680,85 @@ fn main() { } } + // Error correction at circuit level: sample a noisy Clifford circuit, + // derive its detector error model, decode it exactly. + "qec" => { + use wai::quantum_frame::{error_model, parse, sample, surface_code_memory}; + use wai::quantum_match::{fit_code_model, memory_experiment, MatchingGraph}; + let shots: usize = num(args, "--shots", "", 100_000); + let seed: u64 = num(args, "--seed", "", 1); + let json = flag(args, "--json"); + if flag(args, "--fit") { + // The surface code's logical error per round on a grid, and + // the code model a resource estimate would otherwise assume. + let mut pts = Vec::new(); + if !json { + println!("{:>3} {:>7} {:>9} {:>9} {:>11}", "d", "p", "shots", "failures", "per round"); + } + for d in [3u32, 5, 7] { + for p in [0.002, 0.003, 0.004, 0.005] { + let n = if d == 7 { shots / 4 } else { shots }; + let t = memory_experiment(d, d, p, n, seed).unwrap_or_else(|e| die(&format!("{e:?}"))); + if !json { + println!("{:>3} {:>7} {:>9} {:>9} {:>11.3e}", d, p, t.shots, t.failures, t.per_round()); + } + pts.push(t); + } + } + let f = fit_code_model(&pts).unwrap_or_else(|| die("too few failures to fit; raise --shots")); + if json { + println!("{{\"prefactor\":{},\"threshold\":{}}}", f.prefactor, f.threshold); + } else { + println!("fit P(d) = a·(p/p*)^((d+1)/2) per round, a = {:.4}, p* = {:.5}", f.prefactor, f.threshold); + println!("estimate wai-quantum estimate factoring --prefactor {:.4} --threshold {:.5}", f.prefactor, f.threshold); + } + return; + } + let d: u32 = num(args, "--distance", "-d", 3); + let rounds: u32 = num(args, "--rounds", "", d); + let p: f64 = num(args, "--p", "", 0.001); + let text = match opt(args, "--file", "") { + Some(path) => std::fs::read_to_string(&path).unwrap_or_else(|e| die(&format!("{path}: {e}"))), + None => { + if d < 2 || rounds < 1 || !(0.0..=0.5).contains(&p) { + die("need -d ≥ 2, --rounds ≥ 1 and 0 ≤ --p ≤ 0.5"); + } + surface_code_memory(d, rounds, p) + } + }; + if flag(args, "--circuit") { + print!("{text}"); + return; + } + let c = parse(&text).unwrap_or_else(|e| die(&format!("circuit: {e}"))); + let model = error_model(&c).unwrap_or_else(|e| die(&format!("model: {e}"))); + if flag(args, "--dem") { + print!("{}", model.to_text()); + return; + } + let g = MatchingGraph::from_model(&model).unwrap_or_else(|e| die(&format!("matching graph: {e:?}"))); + let s = sample(&c, shots, seed); + let failures = g.failures(&s).unwrap_or_else(|e| die(&format!("decode: {e:?}"))); + let rate = failures as f64 / shots as f64; + if json { + println!( + "{{\"detectors\":{},\"failures\":{failures},\"mechanisms\":{},\"observables\":{},\"qubits\":{},\"seed\":{seed},\"shots\":{shots}}}", + model.detectors, + model.mechanisms.len(), + model.observables, + c.qubits + ); + } else { + println!("circuit {} qubits, {} detectors, {} observables", c.qubits, model.detectors, model.observables); + println!("model {} error mechanisms", model.mechanisms.len()); + println!("failures {failures} of {shots} shots ({rate:.3e} per shot)"); + if opt(args, "--file", "").is_none() { + let per_round = wai::quantum_match::MemoryPoint { distance: d, rounds, p, shots: shots as u64, failures }.per_round(); + println!("per round {per_round:.3e} (d = {d}, {rounds} rounds, p = {p})"); + } + } + } + "entropy" => { let c = circuit_from(args, algo); let sv = c.simulate().unwrap_or_else(|e| die(&format!("simulate: {e:?}"))); diff --git a/wai/wai-rs/src/lib.rs b/wai/wai-rs/src/lib.rs index 5782eddfd..2ee536b01 100644 --- a/wai/wai-rs/src/lib.rs +++ b/wai/wai-rs/src/lib.rs @@ -618,6 +618,12 @@ pub use wai_quantum::quantum_resource; #[cfg(feature = "quantum_resource_receipt")] pub use wai_quantum::quantum_resource_receipt; + +#[cfg(feature = "quantum_frame")] +pub use wai_quantum::quantum_frame; + +#[cfg(feature = "quantum_frame")] +pub use wai_quantum::quantum_match; #[cfg(feature = "quantum_qir")] pub use wai_quantum::quantum_qir; diff --git a/wai/wai-web/Cargo.lock b/wai/wai-web/Cargo.lock index 0da1aeda9..4453ab014 100644 --- a/wai/wai-web/Cargo.lock +++ b/wai/wai-web/Cargo.lock @@ -734,7 +734,7 @@ dependencies = [ [[package]] name = "wai-quantum" -version = "0.3.24" +version = "0.3.25" dependencies = [ "blake3", "ed25519-dalek",