From ed592d76f9f672f4dc9a4afae7e95f064b66f05e Mon Sep 17 00:00:00 2001 From: Transaction Science Date: Sun, 4 Oct 2026 14:21:56 -0400 Subject: [PATCH] wai-quantum 0.3.25: error correction at circuit level, sampled, modelled and decoded exactly quantum_frame takes an error-correction experiment as what it is: a noisy Clifford circuit with detectors and observables. - parse reads the ecosystem's circuit text (resets and measurements in Z and X, H/S/sqrt(X)/Paulis/CX/CZ/SWAP, flip and depolarizing noise, detectors, observables, REPEAT), and refuses anything else rather than skipping it. - sample is bit-packed Pauli-frame simulation, 64 shots per word, with geometric skipping between noise hits; reproducible from (circuit, shots, seed). - error_model derives the detector error model by backward sensitivity propagation. It converts depolarizing channels to equivalent independent mechanisms exactly, XOR-merges identical ones, and refuses non-deterministic detectors. A graphlike mechanism stays one matching edge; otherwise it splits into X and Z halves, then into per-qubit parts. - surface_code_memory generates the rotated surface-code memory experiment. quantum_match decodes by exact minimum-weight perfect matching: Dijkstra on integer weights, with boundary paths computed once and dominated pairs cut, then Edmonds' blossom. memory_experiment and fit_code_model measure P(d) = a(p/p*)^((d+1)/2): a = 0.023, p* = 0.0100 under uniform circuit noise, and CodeFit::surface_code feeds that into a resource estimate. Checked against an independent implementation: - generated circuits identical byte for byte (135 configurations, d 2-15, 1-8 rounds); - error models identical mechanism for mechanism, to 1e-12 in probability (41 circuits: rotated and unrotated surface codes in both bases, repetition codes); - decoded logical error rates within 1.4 standard errors on all 12 grid points. tests/frame_reference.rs keeps a circuit and a model from that implementation as fixtures. Blossom equals brute force on 3000 random graphs. A mirrored random Clifford circuit over every gate type ties the sampler's forward rules to the model's backward ones. Mutating the depolarizing conversion, any gate rule on either side, the determinism check or the decomposition fails a named test. bit_identity pins sampling, decoding and the fit natively and under wasm32-wasip2; repro gains an in-crate exp, within an ulp. CLI: wai-quantum qec [-d --rounds --p --shots --seed --json] | --file | --circuit | --dem | --fit, and estimate --prefactor/--threshold for a measured code. Byte-identical natively and under wasm32-wasip2. Co-Authored-By: Claude Opus 5.5 --- wai/wai-quantum-serve/Cargo.lock | 2 +- wai/wai-quantum/Cargo.lock | 2 +- wai/wai-quantum/Cargo.toml | 15 +- wai/wai-quantum/README.md | 27 + wai/wai-quantum/examples/qec_frame.rs | 65 + wai/wai-quantum/examples/qec_threshold.rs | 26 + wai/wai-quantum/src/lib.rs | 9 +- wai/wai-quantum/src/quantum_frame.rs | 1155 +++++++++++++++++ wai/wai-quantum/src/quantum_match.rs | 1008 ++++++++++++++ wai/wai-quantum/src/repro.rs | 111 +- wai/wai-quantum/tests/bit_identity.rs | 27 +- .../tests/data/memory_x_d3.circuit | 89 ++ wai/wai-quantum/tests/data/memory_x_d3.dem | 221 ++++ .../tests/data/memory_z_d5.circuit | 153 +++ wai/wai-quantum/tests/frame_reference.rs | 46 + wai/wai-rs/Cargo.lock | 2 +- wai/wai-rs/Cargo.toml | 4 + wai/wai-rs/src/bin/wai_quantum.rs | 114 +- wai/wai-rs/src/lib.rs | 6 + wai/wai-web/Cargo.lock | 2 +- 20 files changed, 3062 insertions(+), 22 deletions(-) create mode 100644 wai/wai-quantum/examples/qec_frame.rs create mode 100644 wai/wai-quantum/examples/qec_threshold.rs create mode 100644 wai/wai-quantum/src/quantum_frame.rs create mode 100644 wai/wai-quantum/src/quantum_match.rs create mode 100644 wai/wai-quantum/tests/data/memory_x_d3.circuit create mode 100644 wai/wai-quantum/tests/data/memory_x_d3.dem create mode 100644 wai/wai-quantum/tests/data/memory_z_d5.circuit create mode 100644 wai/wai-quantum/tests/frame_reference.rs 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", -- 2.54.0 (Apple Git-157)