//! Cross-entropy benchmarking of the published random circuits from the //! processor's own bitstrings (Arute et al. 2019, doi:10.5061/dryad.k6t1rj8, //! CC0), on the dense state vector, with a signed receipt. //! //! ```text //! cargo run --release --example xeb_sv --features quantum_noisy,quantum_spd_receipt -- DIR [options] //! ``` //! //! `DIR` is one of the record's unpacked tarballs (`n12_m14/`, …). For every //! `circuit_.qsim` in it with a `measurements_.txt` beside it, this //! reads the circuit and the bitstrings, computes each region's ideal //! amplitudes at those bitstrings, and prints the linear and log XEB with their //! standard errors. Where the record has `amplitudes_.txt` (the authors' //! simulator's amplitudes of the same bitstrings) it compares them with ours //! and gives the XEB theirs make. With `--som FILE` (the record's //! `som_params_by_qubit.csv`) and `circuit_.py` beside the circuit (read as //! text for its qubits' device positions), it gives the digital error model's //! fidelity. At the end, per kind of circuit (full, elided, patch), the //! inverse-variance weighted mean, which is how the record's per-qubit-count //! figures were made, beside the published value (Fig. 4a or 4b). //! //! Options: //! //! - `--kind full|elided|patch`: only that kind (all by default); //! - `--seeds 0-9` or `--seeds 0,3`: only those seeds; //! - `--max-region K`: skip a circuit with a region of more than `K` qubits (the //! state vector's 32 by default), and say so; //! - `--width W`: the fusion width (3); //! - `--gpu`: the single-precision state vector on the platform's GPU (needs the //! `quantum_sv_metal` feature; macOS); //! - `--means E1,E2,EM`: the error model from mean errors instead of the //! record's table: the single-qubit, two-qubit and readout Pauli errors (for //! a record with no per-qubit table); //! - `--reading TAG`: how the circuit files' Z rotations are read //! ([`wai_quantum::quantum_xeb::ZReading`]): `as-written`, //! `divide-single-qubit-z-by-pi`, `drop-fsim-z` or //! `divide-single-qubit-z-by-pi+drop-fsim-z`. By default the 2024 record's //! phase-matched three-patch circuits (`…_patches3_…_phase_match`) are read //! with both corrections, the reading that reproduces them, and every other //! circuit as written. The reading is printed, a circuit outside the form a //! correction needs is refused, and a receipt of any reading but //! `as-written` carries it in its set name (` read `); //! - `--receipt PATH`: where the signed receipt goes; `--table PATH`: one //! tab-separated row a circuit; //! - the signing key is 32 bytes of hex in the file `XEB_SIGNING_KEY` names, or //! a fixed demonstration key. use std::cell::RefCell; use std::collections::HashMap; use std::path::{Path, PathBuf}; use std::rc::Rc; use std::time::Instant; use wai_quantum::ed25519_dalek::SigningKey; use wai_quantum::gates::{C, Gate}; use wai_quantum::quantum_otoc::process_cpu_energy_nj; use wai_quantum::quantum_xeb::{ Circuit, CircuitOutcome, Method, Prober, Shots, SvProber, XebCircuit, XebError, XebJob, ZReading, device_errors, digital_model, estimate, published, qubit_order, seal_metered, weighted_mean, }; /// The record's name for a circuit, parsed: `n12_m14_s0_e0_pEFGH`, /// `patch_n12_m14_s10_e18_pEFGH`, or the 2024 record's /// `n69_m14_s000_patches3_pABCD_num_sq3_[_phase_match]`. #[derive(Clone, Debug)] struct Name { full: String, patch: bool, n: u32, m: u32, s: u32, e: u32, /// The 2024 record's patch count, where the name gives one. patches: Option, /// The coupler activation sequence (`EFGH`, `ABCDCDAB`, `ABCD`). pattern: String, /// The 2024 record's `phase_match`. phase_matched: bool, } fn parse_name(stem: &str) -> Option { let (patch, rest) = match stem.strip_prefix("patch_") { Some(r) => (true, r), None => (false, stem), }; let num = |part: &str, key: &str| part.strip_prefix(key).and_then(|v| v.parse::().ok()); let (mut n, mut m, mut s, mut e, mut patches) = (None, None, None, None, None); let mut pattern = String::new(); let mut extra = String::new(); let mut phase_matched = false; for part in rest.split('_') { if let Some(p) = part.strip_prefix('p').filter(|p| !p.is_empty() && p.bytes().all(|b| b.is_ascii_uppercase())) { pattern = p.to_string(); } // The 2024 record's gate set (`num_sq3`, `num_sq8`) and phase matching. if part.starts_with("sq") && part[2..].parse::().is_ok() { extra += &format!("/{part}"); } if part == "match" { extra += "/phase-matched"; phase_matched = true; } n = n.or(num(part, "n")); m = m.or(num(part, "m")); s = s.or(num(part, "s")); e = e.or(num(part, "e")); patches = patches.or(num(part, "patches")); } Some(Name { full: stem.to_string(), patch, n: n?, m: m?, s: s?, e: e.unwrap_or(0), patches, pattern: pattern + &extra, phase_matched }) } /// The reading a circuit gets unless `--reading` says otherwise: both /// corrections for the 2024 record's phase-matched three-patch circuits, as /// written for every other. fn default_reading(n: &Name) -> ZReading { if n.phase_matched && n.patches == Some(3) { ZReading::DIVIDE_SINGLE_QUBIT_Z_BY_PI_AND_DROP_FSIM_Z } else { ZReading::AS_WRITTEN } } fn kind(n: &Name) -> &'static str { match n.patches { Some(1) => "full", Some(2) => "patches2", Some(3) => "patches3", Some(_) => "patchesK", None if n.patch => "patch", None if n.e == 0 => "full", None => "elided", } } /// A prober that keeps the last region's amplitudes where the job's callback /// reads them, for the comparison with the published ones. struct Keep<'a> { inner: &'a mut dyn Prober, last: Rc>>, } impl Prober for Keep<'_> { fn probe(&mut self, n: u32, gates: &[Gate], at: &[u64]) -> Result<(Vec, f64), XebError> { let (a, i) = self.inner.probe(n, gates, at)?; *self.last.borrow_mut() = a.clone(); Ok((a, i)) } fn method(&self) -> Method { self.inner.method() } } #[cfg(feature = "quantum_sv_metal")] struct GpuProber { width: u32, device: String, } #[cfg(feature = "quantum_sv_metal")] impl Prober for GpuProber { fn probe(&mut self, n: u32, gates: &[Gate], at: &[u64]) -> Result<(Vec, f64), XebError> { use wai_quantum::quantum_sv_metal::GpuState; let mut s = GpuState::zero(n).map_err(|e| XebError::Prober(format!("{e:?}")))?; s.run(gates, self.width as usize).map_err(|e| XebError::Prober(format!("{e:?}")))?; let amps = at.iter().map(|&x| s.amplitude(x as usize)).collect(); let d = (1u64 << n) as f64; let mut sum = 0.0f64; for a in s.amplitudes() { let p = f64::from(a[0]) * f64::from(a[0]) + f64::from(a[1]) * f64::from(a[1]); sum += p * p; } Ok((amps, d * sum - 1.0)) } fn method(&self) -> Method { Method::Gpu32 { width: self.width, device: self.device.clone() } } } /// The published amplitudes of a circuit's bitstrings, by bitstring. fn published_amplitudes(path: &Path, n: u32) -> Option> { let text = std::fs::read_to_string(path).ok()?; let mut map = HashMap::new(); for line in text.lines() { let mut f = line.split_whitespace(); let (b, re, im) = (f.next()?, f.next()?, f.next()?); let x = u64::from_str_radix(b, 2).ok()?; assert_eq!(b.len(), n as usize, "a bitstring of the circuit's width"); let a = C { re: re.parse().ok()?, im: im.parse().ok()? }; if let Some(old) = map.insert(x, a) { assert!(old == a, "one bitstring, one amplitude"); } } Some(map) } struct Row { name: Name, /// Each region's linear XEB and its standard error. regions: Vec<(f64, f64)>, lin: (f64, f64), log: (f64, f64), /// The error model with the simultaneous and with the isolated readout errors. model: Option<(f64, f64)>, theirs: Option<(f64, f64)>, } fn main() { let args: Vec = std::env::args().skip(1).collect(); let dir = PathBuf::from(args.first().expect("usage: xeb_sv DIR [--kind K] [--seeds A-B] [--max-region K] [--width W] [--gpu] [--som CSV] [--receipt PATH] [--table PATH]")); let opt = |k: &str| args.iter().position(|a| a == k).map(|i| args[i + 1].clone()); let width: u32 = opt("--width").map_or(3, |w| w.parse().expect("a width")); let max_region: usize = opt("--max-region").map_or(32, |k| k.parse().expect("a qubit count")); let seeds = opt("--seeds").map(|s| indices(&s)); let only = opt("--kind"); let som = opt("--som").map(|p| std::fs::read_to_string(p).expect("the table reads")); let forced: Option = opt("--reading").map(|t| ZReading::from_tag(&t).unwrap_or_else(|| panic!("--reading {t}: one of as-written, divide-single-qubit-z-by-pi, drop-fsim-z, divide-single-qubit-z-by-pi+drop-fsim-z"))); let means: Option<(f64, f64, f64)> = opt("--means").map(|v| { let x: Vec = v.split(',').map(|t| t.parse().expect("a number")).collect(); (x[0], x[1], x[2]) }); let mut names: Vec = std::fs::read_dir(&dir) .expect("the directory reads") .filter_map(|e| { let f = e.ok()?.file_name().into_string().ok()?; parse_name(f.strip_prefix("circuit_")?.strip_suffix(".qsim")?) }) .filter(|n| only.as_deref().is_none_or(|k| k == kind(n))) .filter(|n| seeds.as_ref().is_none_or(|s| s.contains(&n.s))) .collect(); names.sort_by_key(|n| (kind(n), n.pattern.clone(), n.s, n.e)); let mut circuits = Vec::new(); let mut meta = Vec::new(); // Each reading used, with how many circuits and Z rotations it changed. let mut readings: Vec<(ZReading, usize, u64, u64)> = Vec::new(); for n in names { let written = Circuit::parse(&std::fs::read_to_string(dir.join(format!("circuit_{}.qsim", n.full))).expect("the circuit reads")).expect("a circuit"); let reading = forced.unwrap_or_else(|| default_reading(&n)); let (circuit, changed) = written.read(reading).unwrap_or_else(|e| panic!("{}: the reading {} refuses it: {e} (--reading as-written reads it as written)", n.full, reading.tag())); let largest = circuit.components().iter().map(Vec::len).max().unwrap_or(0); if largest > max_region { println!("{:<34} not run: a region of {largest} qubits, more than --max-region {max_region}", n.full); continue; } let Ok(text) = std::fs::read_to_string(dir.join(format!("measurements_{}.txt", n.full))) else { println!("{:<34} not run: the record has no measurements for it", n.full); continue; }; let shots = Shots::parse_lines(&text, circuit.n).expect("the bitstrings read"); let from_means = means.map(|(e1, e2, em)| { let k = circuit.n as usize; let f = digital_model(&circuit, &vec![e1; k], e2, &vec![em; k]); (f, f) }); let model = from_means.or_else(|| som.as_ref().and_then(|s| { let py = std::fs::read_to_string(dir.join(format!("circuit_{}.py", n.full))).ok()?; let order = qubit_order(&py)?; assert_eq!(order.len(), circuit.n as usize, "a position for every qubit"); let e = device_errors(s, &order)?; let e2 = published::E2_SIMULTANEOUS_MEAN; Some((digital_model(&circuit, &e.single, e2, &e.readout_simultaneous), digital_model(&circuit, &e.single, e2, &e.readout_isolated))) })); match readings.iter_mut().find(|r| r.0 == reading) { Some(r) => { r.1 += 1; r.2 += changed.divided; r.3 += changed.dropped; } None => readings.push((reading, 1, changed.divided, changed.dropped)), } meta.push((n.clone(), model)); circuits.push(XebCircuit { name: n.full.clone(), circuit, shots }); } let set = dir.file_name().map_or("xeb".into(), |s| s.to_string_lossy().into_owned()); // A receipt names any reading but as-written in its set name. let read: String = readings.iter().filter(|r| r.0 != ZReading::AS_WRITTEN).map(|r| format!(" read {}", r.0.tag())).collect(); let job = XebJob { name: format!("{set}{}{read}", only.as_ref().map_or(String::new(), |k| format!("-{k}"))), circuits }; let (key, id) = signer(); let mut sv = SvProber { width }; #[cfg(feature = "quantum_sv_metal")] let mut gpu = GpuProber { width, // The receipt names the device before the first circuit runs. device: wai_quantum::quantum_sv_metal::GpuState::zero(1).map_or_else(|e| format!("no device: {e:?}"), |s| s.device()), }; let inner: &mut dyn Prober = if args.iter().any(|a| a == "--gpu") { #[cfg(feature = "quantum_sv_metal")] { &mut gpu } #[cfg(not(feature = "quantum_sv_metal"))] panic!("--gpu needs the quantum_sv_metal feature") } else { &mut sv }; let kept = Rc::new(RefCell::new(Vec::new())); let mut keep = Keep { inner, last: Rc::clone(&kept) }; println!("set {} · {} circuits · fusion width {width}", job.name, job.circuits.len()); for (r, k, divided, dropped) in &readings { println!("reading {}: {k} circuits · {divided} Z rotations divided by π · {dropped} dropped{}", r.tag(), if forced.is_some() { " (--reading)" } else { " (default for these circuits)" }); } println!("name kind regions shots linear XEB ± se log XEB ± se ideal model sim/iso theirs(linear) max|Δp|·D max|Δa|·√D phase secs CPU J"); let mut rows: Vec = Vec::new(); let mut table = String::from("name\tkind\tn\tm\ts\te\tregions\tshots\tlinear\tlinear_se\tlog\tlog_se\tideal\tmodel_sim\tmodel_iso\ttheirs_linear\ttheirs_linear_se\tmax_dp_D\tmax_da_sqrtD\tphase\tsecs\tcpu_j\n"); let mut last = Instant::now(); let mut energy = process_cpu_energy_nj(); let mut i = 0usize; let (receipt, json) = seal_metered(&key, &id, &job, &mut keep, &mut |c: &XebCircuit, o: &CircuitOutcome| { let secs = last.elapsed().as_secs_f64(); last = Instant::now(); let now = process_cpu_energy_nj(); let joules = match (energy, now) { (Some(a), Some(b)) if b > a => Some((b - a) as f64 / 1e9), _ => None, }; energy = now; let (name, model) = meta[i].clone(); i += 1; let regions: Vec = o.regions.iter().map(|e| e.qubits.len().to_string()).collect(); // The published amplitudes: the record gives them for one-region circuits. let theirs = if o.regions.len() == 1 { published_amplitudes(&dir.join(format!("amplitudes_{}.txt", name.full)), c.circuit.n) } else { None }; let cmp = theirs.map(|map| compare(&c.shots, &kept.borrow(), &map)); let lin = o.linear(); let log = o.log(); let ideal = o.regions.iter().map(|e| format!("{:.4}", e.ideal)).collect::>().join("×"); let f = |v: Option, p: usize| v.map_or("-".to_string(), |v| format!("{v:.p$}")); let g = |v: Option| v.map_or("-".to_string(), |v| format!("{v:.2e}")); println!( "{:<34} {:<7} {:<11} {:>8} {:+.5} ± {:.5} {:+.5} ± {:.5} {:<8} {:<18} {:<18} {:<10} {:<11} {:<9} {secs:5.1} {}", name.full, kind(&name), regions.join("+"), c.shots.total(), lin.0, lin.1, log.0, log.1, ideal, model.map_or("-".into(), |m| format!("{:.5}/{:.5}", m.0, m.1)), cmp.map_or("-".into(), |c| format!("{:+.5} ± {:.5}", c.lin.0, c.lin.1)), g(cmp.map(|c| c.max_dp)), g(cmp.map(|c| c.max_da)), f(cmp.map(|c| c.phase), 5), f(joules, 1), ); table += &format!( "{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{:e}\t{:e}\t{:e}\t{:e}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{secs:.3}\t{}\n", name.full, kind(&name), name.n, name.m, name.s, name.e, regions.join("+"), c.shots.total(), lin.0, lin.1, log.0, log.1, ideal, f(model.map(|m| m.0), 6), f(model.map(|m| m.1), 6), f(cmp.map(|c| c.lin.0), 6), f(cmp.map(|c| c.lin.1), 6), g(cmp.map(|c| c.max_dp)), g(cmp.map(|c| c.max_da)), f(cmp.map(|c| c.phase), 6), f(joules, 3), ); rows.push(Row { name, regions: o.regions.iter().map(|e| (e.linear, e.linear_se)).collect(), lin, log, model, theirs: cmp.map(|c| c.lin) }); }) .expect("the job runs"); println!(); let mut groups: Vec<(&'static str, String)> = Vec::new(); for r in &rows { let g = (kind(&r.name), r.name.pattern.clone()); if !groups.contains(&g) { groups.push(g); } } for (k, pattern) in groups { let these: Vec<&Row> = rows.iter().filter(|r| kind(&r.name) == k && r.name.pattern == pattern).collect(); let (n, m) = (these[0].name.n, these[0].name.m); let lin = weighted_mean(&these.iter().map(|r| r.lin).collect::>()); let log = weighted_mean(&these.iter().map(|r| r.log).collect::>()); let models: Vec<(f64, f64)> = these.iter().filter_map(|r| r.model).collect(); let k_m = models.len() as f64; let model = (!models.is_empty()).then(|| (models.iter().map(|m| m.0).sum::() / k_m, models.iter().map(|m| m.1).sum::() / k_m)); let theirs: Vec<(f64, f64)> = these.iter().filter_map(|r| r.theirs).collect(); // Fig. 4a is the EFGH circuits at 14 cycles; Fig. 4b the 53-qubit ABCDCDAB ones. let pubv: Option<(Option, f64)> = if m == 14 && pattern == "EFGH" { published::fig4a(n).map(|(p, e, f, pred)| { let v = match k { "full" => f, "elided" => Some(e), "patch" => Some(p), _ => None, }; (v, pred) }) } else if n == 53 && pattern == "ABCDCDAB" { published::fig4b(m).map(|(p, e, _, pred)| { let v = match k { "elided" => Some(e), "patch" => Some(p), _ => None, }; (v, pred) }) } else { None }; println!( "{k:<7} p{pattern} n {n} m {m}: {} circuits · weighted linear {:.5} ± {:.5} · weighted log {:.5} ± {:.5}{}{} · published {} · model mean (simultaneous/isolated readout) {} · published prediction {}", these.len(), lin.0, lin.1, log.0, log.1, if theirs.is_empty() { String::new() } else { let t = weighted_mean(&theirs); format!(" · theirs weighted linear {:.5} ± {:.5}", t.0, t.1) }, if k == "patch" { " (each a product over its regions)" } else { "" }, pubv.and_then(|p| p.0).map_or("-".into(), |v| format!("{v}")), model.map_or("-".into(), |v| format!("{:.5}/{:.5}", v.0, v.1)), pubv.map_or("-".into(), |p| format!("{}", p.1)), ); // A circuit of several regions: how many regions' linear XEB is positive, and the smallest in standard errors. let regions: Vec<(f64, f64)> = these.iter().flat_map(|r| r.regions.iter().copied()).collect(); if regions.len() > these.len() { let z_min = regions.iter().map(|&(v, s)| v / s).fold(f64::INFINITY, f64::min); println!("{k:<7} p{pattern} n {n} m {m}: regions with a positive linear XEB {} of {} · smallest linear XEB / se {z_min:+.1}", regions.iter().filter(|r| r.0 > 0.0).count(), regions.len()); } } let work: u64 = receipt.outcomes.iter().map(|o| o.work).sum(); println!("work {work} amplitude-gate updates · {:.1} s · {} µJ sealed (this process's CPU energy; 0 = unmetered)", receipt.elapsed_ms as f64 / 1000.0, receipt.joules_micro); if let Some(out) = opt("--receipt") { std::fs::write(&out, json + "\n").expect("the receipt writes"); println!("receipt {out} (receipt hash {})", hex(&receipt.receipt_hash())); } if let Some(out) = opt("--table") { std::fs::write(&out, table).expect("the table writes"); } } #[derive(Clone, Copy)] struct Comparison { /// The linear XEB the published amplitudes give, with its standard error. lin: (f64, f64), /// The largest difference of a probability, in units of `1/D`. max_dp: f64, /// The largest difference of an amplitude after the best global phase, in /// units of `1/√D`. max_da: f64, /// That global phase, in radians. phase: f64, } fn compare(shots: &Shots, ours: &[C], theirs: &HashMap) -> Comparison { let n = shots.n; let d = (1u64 << n) as f64; // A one-region circuit is within the state vector's 32 qubits, so its bitstrings fit a u64. let pairs: Vec<(C, C, u32)> = shots.counts.iter().zip(ours).map(|(&(x, c), &a)| (a, *theirs.get(&(x as u64)).expect("a published amplitude for every bitstring"), c)).collect(); // The global phase that best aligns ours to theirs: arg Σ conj(a)·b. let (sr, si) = pairs.iter().fold((0.0, 0.0), |(r, i), (a, b, _)| (r + a.re * b.re + a.im * b.im, i + a.re * b.im - a.im * b.re)); let norm = (sr * sr + si * si).sqrt(); let (cr, ci) = (sr / norm, si / norm); let mut max_dp = 0.0f64; let mut max_da = 0.0f64; for (a, b, _) in &pairs { let (pa, pb) = (a.re * a.re + a.im * a.im, b.re * b.re + b.im * b.im); max_dp = max_dp.max((pa - pb).abs() * d); let (rr, ri) = (a.re * cr - a.im * ci, a.re * ci + a.im * cr); max_da = max_da.max(((rr - b.re).powi(2) + (ri - b.im).powi(2)).sqrt() * d.sqrt()); } let probs: Vec<(f64, u32)> = pairs.iter().map(|(_, b, c)| (b.re * b.re + b.im * b.im, *c)).collect(); let e = estimate((0..n).collect(), &probs, f64::NAN).expect("shots"); Comparison { lin: (e.linear, e.linear_se), max_dp, max_da, phase: si.atan2(sr) } } fn indices(s: &str) -> Vec { s.split(',') .flat_map(|p| match p.split_once('-') { Some((a, b)) => (a.parse::().expect("an index")..=b.parse().expect("an index")).collect::>(), None => vec![p.parse().expect("an index")], }) .collect() } fn signer() -> (SigningKey, String) { match std::env::var("XEB_SIGNING_KEY") { Ok(p) => { let h = std::fs::read_to_string(&p).expect("the key file reads"); let h = h.trim(); let b: Vec = (0..32).map(|i| u8::from_str_radix(&h[2 * i..2 * i + 2], 16).expect("hex")).collect(); (SigningKey::from_bytes(&b.try_into().expect("32 bytes")), "wai-quantum xeb reproduction".into()) } Err(_) => (SigningKey::from_bytes(&[0x0b; 32]), "wai-quantum xeb example (demonstration key)".into()), } } fn hex(b: &[u8]) -> String { b.iter().map(|x| format!("{x:02x}")).collect() }