From a43f215aa6be45bca66d4c378c8e476b49d184ef Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Wed, 8 Jul 2026 10:37:46 +0200 Subject: [PATCH 1/6] feat(tableau): add Pauli expectation/trace and reset_all MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add `expectation.rs` for `GeneralizedTableau`: - `expectation` — single-Pauli ⟨ψ|P|ψ⟩ for a `PauliWord`, conjugating P through the tableau and overlapping with the sparse coefficient vector. - `trace` — Σ over Paulis matching a `PauliPattern`. Supporting changes in `data.rs`: - `compute_decomposition_word` — multi-qubit conjugation used by `expectation` (generalizes the single-qubit `compute_decomposition`). - `reset_all` on both `Tableau` and `GeneralizedTableau` — reinitialise to |0…0⟩, clearing coefficients and loss state (refactors `new` via a shared `new_data` helper). - the rayon-nesting guard on the coefficient path (avoids nesting rayon inside shot-level parallelism; a no-op on the single-shot path). `lib.rs` gains `pub mod expectation`; `Cargo.toml` gains a `ppvm-pauli-word` dependency. Pure core addition, snapshotted from #168. Merge order for the #168 split: this and single-qubit rotations target `main` independently; native .pyi stubs stack on rotations; then vihaco-circuit-isa (#169) -> circuit component -> composite -> CLI/TUI (#166). Co-Authored-By: Claude Opus 4.8 (1M context) --- Cargo.lock | 1 + crates/ppvm-tableau/Cargo.toml | 1 + crates/ppvm-tableau/src/data.rs | 61 ++++- crates/ppvm-tableau/src/expectation.rs | 316 +++++++++++++++++++++++++ crates/ppvm-tableau/src/lib.rs | 2 + 5 files changed, 378 insertions(+), 3 deletions(-) create mode 100644 crates/ppvm-tableau/src/expectation.rs diff --git a/Cargo.lock b/Cargo.lock index caa60ec3c..6370265da 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1020,6 +1020,7 @@ dependencies = [ "itertools 0.14.0", "num", "ppvm-pauli-sum", + "ppvm-pauli-word", "ppvm-traits", "rand 0.10.1", "rayon", diff --git a/crates/ppvm-tableau/Cargo.toml b/crates/ppvm-tableau/Cargo.toml index 72826d5a9..6b5c4dc07 100644 --- a/crates/ppvm-tableau/Cargo.toml +++ b/crates/ppvm-tableau/Cargo.toml @@ -17,6 +17,7 @@ ppvm-traits = { version = "0.1.0", path = "../ppvm-traits" } ppvm-pauli-sum = { version = "0.1.0", path = "../ppvm-pauli-sum" } rand = "0.10.1" smallvec = "1.15" +ppvm-pauli-word = { version = "0.1.0", path = "../ppvm-pauli-word" } # rayon needs OS threads — native only. [target.'cfg(not(target_arch = "wasm32"))'.dependencies] diff --git a/crates/ppvm-tableau/src/data.rs b/crates/ppvm-tableau/src/data.rs index d35e3e700..52443a980 100644 --- a/crates/ppvm-tableau/src/data.rs +++ b/crates/ppvm-tableau/src/data.rs @@ -54,8 +54,7 @@ pub struct Tableau { } impl Tableau { - /// Construct a fresh tableau initialised to `|0…0⟩`. - pub fn new(n_qubits: usize) -> Self { + fn new_data(n_qubits: usize) -> Vec> { // Initialize tableau for 0 state let mut data: Vec> = Vec::with_capacity(2 * n_qubits); @@ -72,7 +71,11 @@ impl Tableau { pw.set(i, Pauli::Z); data.push(pw); } + data + } + pub fn new(n_qubits: usize) -> Self { + let data = Tableau::::new_data(n_qubits); Self { n_qubits, data, @@ -87,6 +90,11 @@ impl Tableau { t } + pub fn reset_all(&mut self) { + let data = Tableau::::new_data(self.n_qubits); + self.data = data; + } + /// View of the stabilizer rows (the upper half of the tableau). #[inline] pub fn stabilizers(&self) -> &[PhasedPauliWordNoHash] { @@ -537,7 +545,9 @@ where Complex: std::ops::Mul> + std::ops::AddAssign + From + Copy, { - if items.len() >= RAYON_COEFF_THRESHOLD { + // See `branch_coefficients_parallel`: avoid nesting rayon inside shot-level + // parallelism; the main-thread (single-shot) path is unaffected. + if items.len() >= RAYON_COEFF_THRESHOLD && rayon::current_thread_index().is_none() { use rayon::prelude::*; return items @@ -675,6 +685,21 @@ where s } + pub fn reset_all(&mut self) { + self.tableau.reset_all(); + + let mut coefficients = C::new(); + let complex_one = Complex { + re: T::Coeff::one(), + im: T::Coeff::zero(), + }; + coefficients.unsafe_insert(I::zero(), complex_one); + self.coefficients = coefficients; + for l in self.is_lost.iter_mut() { + *l &= false; + } + } + /// Clone the quantum state but reinitialize the RNG, producing an independent simulation /// branch. If `seed` is `Some`, the new RNG is seeded deterministically; if `None`, it is /// seeded from OS entropy. @@ -848,6 +873,36 @@ where (p_word.phase, stab_anticomm_bits, destab_anticomm_bits) } + /// Multi-qubit generalization of [`compute_decomposition`]: conjugate an + /// arbitrary `PauliWord` through the tableau and return the same triple + /// `(phase, stab_anticomm_bits, destab_anticomm_bits)`. + /// + /// Algorithm: call [`compute_decomposition`] for each non-identity qubit + /// in the input, then multiply the resulting single-qubit conjugates in + /// canonical-basis form `i^φ X^x Z^z`. Pauli multiplication picks up a + /// `(-1)^{popcount(z_running & x_new)}` cross-phase from + /// `Z^z_a X^x_b = (-1)^{z_a · x_b} X^x_b Z^z_a`. + pub(crate) fn compute_decomposition_word(&self, word: &W) -> (u8, I, I) + where + <::Storage as BitView>::Store: PrimInt, + { + let mut phase = 0u8; + let mut stab_anticomm = I::zero(); + let mut destab_anticomm = I::zero(); + for q in 0..self.n_qubits() { + let p_q = word.get(q); + if p_q == Pauli::I { + continue; + } + let (q_phase, q_stab, q_destab) = self.compute_decomposition(q, p_q); + let cross = 2 * (symplectic_inner(destab_anticomm, q_stab) as u8 % 2); + phase = (phase + q_phase + cross) % 4; + stab_anticomm = stab_anticomm ^ q_stab; + destab_anticomm = destab_anticomm ^ q_destab; + } + (phase, stab_anticomm, destab_anticomm) + } + /// every basis index is a bit string alpha defining the basis state /// the phase when applying a Pauli is the product of all destabilizer phases /// and the phase contributions from the commutation relations diff --git a/crates/ppvm-tableau/src/expectation.rs b/crates/ppvm-tableau/src/expectation.rs new file mode 100644 index 000000000..9a69d596c --- /dev/null +++ b/crates/ppvm-tableau/src/expectation.rs @@ -0,0 +1,316 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Pauli-string expectation values for `GeneralizedTableau`. +//! +//! Two entry points: +//! +//! - [`GeneralizedTableau::expectation`] — single-Pauli `⟨ψ|P|ψ⟩` for a +//! `PauliWord`. Conjugates `P` through the tableau and overlaps the +//! resulting Pauli with the sparse coefficient vector using the same +//! formulas as the measurement code. +//! - [`GeneralizedTableau::trace`] — `Σ_{P matches pattern} ⟨ψ|P|ψ⟩` for a +//! `PauliPattern`. Enumerates the matching Paulis and sums their +//! expectations. +//! +//! Decision 9 in the multi-backend plan calls these out as the natural +//! primitive for the tableau backend; semantics intentionally diverge from +//! the PauliSum trace. + +use crate::data::GeneralizedTableau; +use crate::prelude::*; +use bitvec::view::BitView; +use fxhash::FxHashMap as HashMap; +use num::PrimInt; +use num::complex::{Complex, Complex64, ComplexFloat}; +use num::traits::{One, ToPrimitive, Zero}; +use ppvm_pauli_word::pattern::PauliPattern; +use std::fmt::Debug; + +impl GeneralizedTableau +where + T: Config, + <::Storage as BitView>::Store: PrimInt, + C: SparseVector, I> + Debug, + T::Coeff: One + + Zero + + Clone + + num::Num + + ToPrimitive + + Debug + + std::ops::Mul + + PartialOrd + + Send + + Sync, + Complex: std::ops::Mul> + + From + + std::ops::MulAssign + + std::ops::AddAssign + + One + + ComplexFloat + + Copy, + I: TableauIndex + Debug + Send + Sync, +{ + /// `⟨ψ|word|ψ⟩` for the multi-qubit Pauli `word`. + /// + /// Conjugates `word` through the Clifford tableau (giving a Pauli on the + /// canonical basis: an X-mask, Z-mask, and `i^φ` phase), then sums + /// `⟨α|P_conj|β⟩ c_α* c_β` over the sparse coefficient vector. Always + /// returns a real number (Hermitian operator on a normalized state). + pub fn expectation(&self, word: &W) -> f64 { + let (phase, stab_anticomm, destab_anticomm) = self.compute_decomposition_word(word); + let entries: Vec<(Complex, I)> = self.coefficients.clone().into_iter().collect(); + if stab_anticomm == I::zero() { + Self::compute_overlap_case_b(&entries, phase, destab_anticomm) + } else { + let coeff_map: HashMap> = + entries.into_iter().map(|(v, i)| (i, v)).collect(); + let odd_phase_mask = self.odd_phase_destabilizer_mask(); + Self::compute_overlap_case_a( + &coeff_map, + phase, + destab_anticomm, + stab_anticomm, + odd_phase_mask, + ) + } + } + + /// `Σ_{P matches pattern} ⟨ψ|P|ψ⟩`. + /// + /// Enumerates every `PauliWord` accepted by `pattern` via + /// [`PauliPattern::enumerate_matches`] and sums their expectations. + /// Star quantifiers (`X*`) panic — the pattern must be bounded; use + /// counted repetition (`Z?{n}`) or positional anchors instead. + pub fn trace(&self, pattern: &PauliPattern) -> f64 { + let mut sum = 0.0f64; + for word in pattern.enumerate_matches::(self.n_qubits()) { + sum += self.expectation(&word); + } + sum + } +} + +#[cfg(test)] +mod tests { + use super::*; + use ppvm_pauli_sum::config::fxhash::ByteF64; + + type TestTableau = GeneralizedTableau>; + + fn word(s: &str) -> PauliWord { + s.into() + } + + fn assert_close(actual: f64, expected: f64, tol: f64) { + assert!( + (actual - expected).abs() < tol, + "expected {expected}, got {actual} (|Δ| = {})", + (actual - expected).abs() + ); + } + + // ─── Single-qubit expectations ────────────────────────────────────── + + #[test] + fn expectation_z_on_zero_state_is_one() { + let tab: TestTableau = GeneralizedTableau::new(1, 1e-12); + assert_close(tab.expectation(&word("Z")), 1.0, 1e-12); + } + + #[test] + fn expectation_x_on_zero_state_is_zero() { + let tab: TestTableau = GeneralizedTableau::new(1, 1e-12); + assert_close(tab.expectation(&word("X")), 0.0, 1e-12); + } + + #[test] + fn expectation_identity_on_zero_state_is_one() { + let tab: TestTableau = GeneralizedTableau::new(1, 1e-12); + assert_close(tab.expectation(&word("I")), 1.0, 1e-12); + } + + #[test] + fn expectation_x_on_plus_state_is_one() { + let mut tab: TestTableau = GeneralizedTableau::new(1, 1e-12); + tab.h(0); + assert_close(tab.expectation(&word("X")), 1.0, 1e-12); + assert_close(tab.expectation(&word("Z")), 0.0, 1e-12); + } + + // ─── Bell state ⟨Φ+|·|Φ+⟩ ──────────────────────────────────────────── + + fn bell() -> TestTableau { + let mut tab: TestTableau = GeneralizedTableau::new(2, 1e-12); + tab.h(0); + tab.cnot(0, 1); + tab + } + + #[test] + fn bell_state_pauli_expectations() { + let tab = bell(); + assert_close(tab.expectation(&word("II")), 1.0, 1e-12); + assert_close(tab.expectation(&word("ZZ")), 1.0, 1e-12); + assert_close(tab.expectation(&word("XX")), 1.0, 1e-12); + assert_close(tab.expectation(&word("YY")), -1.0, 1e-12); + // Cross terms vanish for the Bell state. + assert_close(tab.expectation(&word("IZ")), 0.0, 1e-12); + assert_close(tab.expectation(&word("ZI")), 0.0, 1e-12); + assert_close(tab.expectation(&word("XZ")), 0.0, 1e-12); + assert_close(tab.expectation(&word("YX")), 0.0, 1e-12); + } + + // ─── GHZ state ──────────────────────────────────────────────────── + + #[test] + fn ghz_state_expectations() { + let mut tab: TestTableau = GeneralizedTableau::new(3, 1e-12); + tab.h(0); + tab.cnot(0, 1); + tab.cnot(1, 2); + // GHZ = (|000⟩ + |111⟩)/√2. For Z^z: eigenvalue is (-1)^{popcount(z)·x} + // on |xxx⟩, so on the two basis states it agrees iff popcount(z) is + // even; the diagonal expectation is then +1, otherwise 0. + assert_close(tab.expectation(&word("III")), 1.0, 1e-12); // popcount 0 → +1 + assert_close(tab.expectation(&word("ZZZ")), 0.0, 1e-12); // popcount 3 → 0 + assert_close(tab.expectation(&word("ZIZ")), 1.0, 1e-12); // popcount 2 → +1 + assert_close(tab.expectation(&word("ZZI")), 1.0, 1e-12); // popcount 2 → +1 + assert_close(tab.expectation(&word("IZI")), 0.0, 1e-12); // popcount 1 → 0 + // XXX flips |000⟩ ↔ |111⟩, both in the GHZ superposition → +1. + assert_close(tab.expectation(&word("XXX")), 1.0, 1e-12); + // Y has off-diagonal action with imaginary phase; YYY contributes 0. + assert_close(tab.expectation(&word("YYY")), 0.0, 1e-12); + } + + // ─── Single-qubit rotation: |ψ⟩ = RY(θ)|0⟩ ──────────────────────── + + #[test] + fn ry_rotation_z_expectation_is_cos_theta() { + // RY(θ)|0⟩ = cos(θ/2)|0⟩ + sin(θ/2)|1⟩. ⟨ψ|Z|ψ⟩ = cos(θ). + for theta in [0.0, 0.3, 1.0, std::f64::consts::PI / 2.0] { + let mut tab: TestTableau = GeneralizedTableau::new(1, 1e-12); + tab.ry(0, theta); + assert_close(tab.expectation(&word("Z")), theta.cos(), 1e-12); + assert_close(tab.expectation(&word("X")), theta.sin(), 1e-12); + } + } + + // ─── X / Y on a non-Clifford superposition ─────────────────────── + // + // After H(0); T(0) the state is |ψ⟩ = (|0⟩ + e^{iπ/4}|1⟩)/√2: + // ⟨ψ|X|ψ⟩ = cos(π/4) = √2/2, + // ⟨ψ|Y|ψ⟩ = sin(π/4) = √2/2. + // The T gate populates two branches in the sparse coefficient vector, + // so the overlap is a cross-product between them — case_a, not case_b. + // Y in particular drives `phase_decomp` to an odd value (Y on a |+⟩-style + // frame conjugates to -Y), forcing the `phase == 1 | 3` arms of + // `overlap_case_a`'s `match`. A sign bug there would flip ⟨Y⟩. + + #[test] + fn t_plus_x_expectation_is_cos_pi_over_4() { + let mut tab: TestTableau = GeneralizedTableau::new(1, 1e-12); + tab.h(0); + tab.t(0); + assert_close( + tab.expectation(&word("X")), + std::f64::consts::FRAC_1_SQRT_2, + 1e-12, + ); + } + + #[test] + fn t_plus_y_expectation_is_sin_pi_over_4() { + let mut tab: TestTableau = GeneralizedTableau::new(1, 1e-12); + tab.h(0); + tab.t(0); + assert_close( + tab.expectation(&word("Y")), + std::f64::consts::FRAC_1_SQRT_2, + 1e-12, + ); + } + + // ─── trace(pattern) ─────────────────────────────────────────────── + + #[test] + fn trace_of_z_or_identity_pattern_on_bell_is_two() { + // For |Φ+⟩ = (|00⟩+|11⟩)/√2: + // Σ_{P ∈ {I,Z}^2} ⟨Φ+|P|Φ+⟩ = ⟨II⟩ + ⟨IZ⟩ + ⟨ZI⟩ + ⟨ZZ⟩ + // = 1 + 0 + 0 + 1 = 2. + // Equivalently 2^n |⟨0…0|Φ+⟩|² = 4 · 1/2 = 2. + let tab = bell(); + let pat = PauliPattern::parse("Z?{2}").expect("parse Z?{2}"); + assert_close(tab.trace(&pat), 2.0, 1e-12); + } + + #[test] + fn trace_of_y_or_identity_pattern_on_bell_is_zero() { + // For |Φ+⟩, Σ_{P ∈ {I,Y}^2} ⟨Φ+|P|Φ+⟩ = ⟨II⟩ + ⟨IY⟩ + ⟨YI⟩ + ⟨YY⟩ + // = 1 + 0 + 0 + (-1) = 0. + // Equivalently 2^n |⟨+i+i|Φ+⟩|² = 4 · 0 = 0 — the projection onto + // the all-|+i⟩ state has zero amplitude. `trace` does the + // enumeration sum here without ever mutating state or calling + // `normalize`, so the zero-probability case doesn't panic. + let tab = bell(); + let pat = PauliPattern::parse("Y?{2}").expect("parse Y?{2}"); + assert_close(tab.trace(&pat), 0.0, 1e-12); + } + + #[test] + fn trace_of_positional_pattern_on_bell_matches_single_pauli() { + // `Z0Z1` matches exactly the word ZZ; trace should equal ⟨ZZ⟩ = 1. + let tab = bell(); + let pat = PauliPattern::parse("Z0Z1").expect("parse Z0Z1"); + assert_close(tab.trace(&pat), 1.0, 1e-12); + } + + // ─── Cross-backend: forward tableau shots vs backward PauliSum ──────── + // + // Run a noisy Clifford circuit many times on a tableau and average the + // per-shot ⟨ψ|ZZ|ψ⟩. Independently, Heisenberg-propagate ZZ backward + // through the same circuit via `PauliSum` and read off ⟨0…0|U†ZZ U|0…0⟩ + // as the sum of coefficients over Z/I-only Paulis. The Monte-Carlo + // average and the deterministic value must agree within sampling error. + // + // `PauliSum::g(i)` performs O → g† O g, so the backward sweep applies + // gates in reverse time order. Depolarize is self-dual under Heisenberg. + #[test] + fn forward_shots_match_backward_pauli_sum_under_depolarizing_noise() { + use ppvm_pauli_sum::config::indexmap::ByteFxHashF64; + + let p = 0.05_f64; + let n_shots: u64 = 4000; + let n_qubits = 2; + + let mut sum = 0.0_f64; + for shot in 0..n_shots { + let mut tab: TestTableau = GeneralizedTableau::new_with_seed(n_qubits, 1e-12, shot); + tab.h(0); + tab.depolarize1(0, p); + tab.cnot(0, 1); + tab.depolarize1(0, p); + tab.depolarize1(1, p); + sum += tab.expectation(&word("ZZ")); + } + let avg = sum / (n_shots as f64); + + let mut ps: PauliSum> = PauliSum::builder().n_qubits(n_qubits).build(); + ps += ("ZZ", 1.0); + ps.depolarize1(1, p); + ps.depolarize1(0, p); + ps.cnot(0, 1); + ps.depolarize1(0, p); + ps.h(0); + let z_or_i = PauliPattern::parse("Z?{2}").expect("parse Z?{2}"); + let exact = ps.trace(&z_or_i); + + // Per-shot |⟨ZZ⟩| ≤ 1 ⇒ σ_mean ≤ 1/√N; 5σ keeps this robust to RNG draws. + let tol = 5.0 / (n_shots as f64).sqrt(); + assert!( + (avg - exact).abs() < tol, + "tableau avg {avg} vs PauliSum exact {exact}, |Δ|={} (tol {tol})", + (avg - exact).abs() + ); + } +} diff --git a/crates/ppvm-tableau/src/lib.rs b/crates/ppvm-tableau/src/lib.rs index 1738e7ac2..76f51a881 100644 --- a/crates/ppvm-tableau/src/lib.rs +++ b/crates/ppvm-tableau/src/lib.rs @@ -34,6 +34,8 @@ pub mod data; /// `Display` implementations for tableau types. pub mod display; +/// Pauli-string expectation values and pattern traces. +pub mod expectation; /// Gate implementations (Clifford, T, rotations). pub mod gates; /// Z-basis measurement, including loss-aware variants. From d43fda1c4c6572bbf87f096a504498d853ef6196 Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Wed, 8 Jul 2026 11:14:39 +0200 Subject: [PATCH 2/6] fix(tableau): address review on expectation/reset_all - `expectation` now iterates `self.coefficients` by reference (mirroring `z_expectation`) instead of cloning the whole sparse vector and collecting it up front. The map path builds directly from `.iter()`, dropping an intermediate `Vec`. `trace` calls `expectation` per matching Pauli, so this removes a per-call clone. Numerically identical (all expectation/trace tests unchanged). - `GeneralizedTableau::reset_all` now clears `measurement_record`, so `current_measurement_record()` no longer returns stale outcomes after a full reset. - Restore the `Tableau::new` doc comment dropped in the `new_data` refactor. Co-Authored-By: Claude Opus 4.8 (1M context) --- crates/ppvm-tableau/src/data.rs | 2 + crates/ppvm-tableau/src/expectation.rs | 4 +- crates/ppvm-vihaco/src/component.rs | 935 +++++++++++++++++++++++++ crates/ppvm-vihaco/src/measurements.rs | 93 +++ 4 files changed, 1032 insertions(+), 2 deletions(-) create mode 100644 crates/ppvm-vihaco/src/component.rs create mode 100644 crates/ppvm-vihaco/src/measurements.rs diff --git a/crates/ppvm-tableau/src/data.rs b/crates/ppvm-tableau/src/data.rs index 52443a980..386126343 100644 --- a/crates/ppvm-tableau/src/data.rs +++ b/crates/ppvm-tableau/src/data.rs @@ -74,6 +74,7 @@ impl Tableau { data } + /// Construct a fresh tableau initialised to `|0…0⟩`. pub fn new(n_qubits: usize) -> Self { let data = Tableau::::new_data(n_qubits); Self { @@ -698,6 +699,7 @@ where for l in self.is_lost.iter_mut() { *l &= false; } + self.measurement_record.clear(); } /// Clone the quantum state but reinitialize the RNG, producing an independent simulation diff --git a/crates/ppvm-tableau/src/expectation.rs b/crates/ppvm-tableau/src/expectation.rs index 9a69d596c..9efe04ea0 100644 --- a/crates/ppvm-tableau/src/expectation.rs +++ b/crates/ppvm-tableau/src/expectation.rs @@ -59,12 +59,12 @@ where /// returns a real number (Hermitian operator on a normalized state). pub fn expectation(&self, word: &W) -> f64 { let (phase, stab_anticomm, destab_anticomm) = self.compute_decomposition_word(word); - let entries: Vec<(Complex, I)> = self.coefficients.clone().into_iter().collect(); if stab_anticomm == I::zero() { + let entries: Vec<(Complex, I)> = self.coefficients.iter().copied().collect(); Self::compute_overlap_case_b(&entries, phase, destab_anticomm) } else { let coeff_map: HashMap> = - entries.into_iter().map(|(v, i)| (i, v)).collect(); + self.coefficients.iter().map(|&(c, i)| (i, c)).collect(); let odd_phase_mask = self.odd_phase_destabilizer_mask(); Self::compute_overlap_case_a( &coeff_map, diff --git a/crates/ppvm-vihaco/src/component.rs b/crates/ppvm-vihaco/src/component.rs new file mode 100644 index 000000000..ae9e70806 --- /dev/null +++ b/crates/ppvm-vihaco/src/component.rs @@ -0,0 +1,935 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +use crate::composite::PPVMDeviceInfo; +use crate::measurements::{ + CircuitOutcomeEffect, MeasurementEffect, MeasurementOutcome, TraceEffect, +}; +use bitvec::view::BitView; +use bnum::types::{U256, U512, U1024, U2048}; +use eyre::{Result, eyre}; +use num::PrimInt; +use num::complex::Complex64; +use ppvm_pauli_sum::config::fx64hash::Byte8F64; +use ppvm_pauli_sum::config::indexmap::ByteFxHashF64; +use ppvm_pauli_sum::strategy::{CoefficientThreshold, CombinedStrategy, MaxPauliWeight}; +use ppvm_tableau::prelude::*; +use vihaco::{Effects, component, observe}; +use vihaco_circuit_isa::{CircuitEffect, CircuitInstruction, CircuitMessage}; + +/// Truncation strategy used by every `PauliSum` / `LossyPauliSum` size bucket. +/// Coefficient-threshold pruning is always on; the Pauli-weight cap is set per +/// run from the header (defaults to `usize::MAX` = no cap). +type PauliSumStrategy = CombinedStrategy; + +/// `PauliSum`'s `T` for the lossless backend: `[u8; N]` storage, fx hash, +/// f64 coefficients, the strategy above. +type PauliSumConfig = ByteFxHashF64; + +/// Same as `PauliSumConfig` but with `LossyPauliWord` as the word type, so the +/// loss-channel methods are dispatchable on the resulting `PauliSum`. +/// `LossyPauliWord`'s second type parameter (hasher) defaults to +/// `fxhash::FxBuildHasher`, matching `ByteFxHashF64`'s internal hasher. +type LossyPauliSumConfig = + ByteFxHashF64>; + +/// Build a `PauliSumStrategy` value from a `PPVMDeviceInfo`. Pulled out so the +/// six size-bucket constructors don't each repeat the strategy spelling. +fn paulisum_strategy(info: &PPVMDeviceInfo) -> PauliSumStrategy { + CombinedStrategy( + CoefficientThreshold(info.coefficient_threshold), + MaxPauliWeight(info.max_pauli_weight.unwrap_or(usize::MAX)), + ) +} + +macro_rules! batch_for { + ($tab:expr, $method:ident, $addrs:expr) => { + for addr in $addrs { $tab.$method(*addr); } + }; + ($tab:expr, $method:ident, $addrs:expr, $($arg:expr),+) => { + for addr in $addrs { $tab.$method(*addr, $($arg),+); } + }; +} + +/// Two-qubit sibling of [`batch_for!`]: drives a method that takes two qubit +/// addresses (plus optional extra args) over a slice of `(usize, usize)` pairs. +macro_rules! batch_pairs_for { + ($state:expr, $method:ident, $pairs:expr) => { + for &(a, b) in $pairs { $state.$method(a, b); } + }; + ($state:expr, $method:ident, $pairs:expr, $($arg:expr),+) => { + for &(a, b) in $pairs { $state.$method(a, b, $($arg),+); } + }; +} + +pub struct CircuitExecutor, I: TableauIndex, C: SparseVector> { + pub tab: GeneralizedTableau, +} + +#[component(instruction = CircuitInstruction, message = CircuitMessage, effect = CircuitOutcomeEffect)] +impl CircuitExecutor +where + T: Config, + <::Storage as BitView>::Store: PrimInt, + I: TableauIndex + Send + Sync + std::fmt::Debug, + C: SparseVector + std::fmt::Debug, +{ + fn execute( + &mut self, + inst: CircuitInstruction, + msg: CircuitMessage, + ) -> Result> { + self.execute_instruction(&inst, &msg) + } + + fn execute_instruction( + &mut self, + inst: &CircuitInstruction, + msg: &CircuitMessage, + ) -> Result> { + use CircuitInstruction::*; + use CircuitMessage::*; + + match (inst, msg) { + // Single-qubit Clifford + (X, &Qubit(addr)) => self.tab.x(addr), + (Y, &Qubit(addr)) => self.tab.y(addr), + (Z, &Qubit(addr)) => self.tab.z(addr), + (H, &Qubit(addr)) => self.tab.h(addr), + (S, &Qubit(addr)) => self.tab.s(addr), + (SAdj, &Qubit(addr)) => self.tab.s_dag(addr), + (SqrtX, &Qubit(addr)) => self.tab.sqrt_x(addr), + (SqrtY, &Qubit(addr)) => self.tab.sqrt_y(addr), + (SqrtXAdj, &Qubit(addr)) => self.tab.sqrt_x_dag(addr), + (SqrtYAdj, &Qubit(addr)) => self.tab.sqrt_y_dag(addr), + + // Controlled gates + (CNOT, &TwoQubit(addr0, addr1)) => self.tab.cnot(addr0, addr1), + (CZ, &TwoQubit(addr0, addr1)) => self.tab.cz(addr0, addr1), + + // T gate + (T, &Qubit(addr)) => self.tab.t(addr), + (TAdj, &Qubit(addr)) => self.tab.t_dag(addr), + + // Single-qubit rotations + (RX, &QubitAndFloat(addr, angle)) => self.tab.rx(addr, angle), + (RY, &QubitAndFloat(addr, angle)) => self.tab.ry(addr, angle), + (RZ, &QubitAndFloat(addr, angle)) => self.tab.rz(addr, angle), + + // Two-qubit rotations + (RXX, &TwoQubitAndFloat(addr0, addr1, angle)) => self.tab.rxx(addr0, addr1, angle), + (RYY, &TwoQubitAndFloat(addr0, addr1, angle)) => self.tab.ryy(addr0, addr1, angle), + (RZZ, &TwoQubitAndFloat(addr0, addr1, angle)) => self.tab.rzz(addr0, addr1, angle), + + // U3 + (U3, &QubitU3(addr, theta, phi, lam)) => self.tab.u3(addr, theta, phi, lam), + + // RXY: rotation about an axis in the x/y plane + (R, &QubitAndTwoFloats(addr, axis_angle, theta)) => self.tab.r(addr, axis_angle, theta), + + // Measure & Reset + (Measure, &Qubit(addr)) => { + let outcome: MeasurementOutcome = self.tab.measure(addr).into(); + return Ok(Effects::one(CircuitOutcomeEffect::Measurement( + MeasurementEffect { + measurement_results: smallvec::smallvec![outcome], + }, + ))); + } + (Reset, &Qubit(addr)) => self.tab.reset(addr), + + // Noise + (Depolarize, &QubitAndFloat(addr, p)) => self.tab.depolarize1(addr, p), + (Depolarize2, &TwoQubitAndFloat(addr0, addr1, p)) => { + self.tab.depolarize2(addr0, addr1, p) + } + (PauliError, QubitAndFloatArr3(addr0, ps)) => self.tab.pauli_error(*addr0, *ps), + (TwoQubitPauliError, TwoQubitAndFloatArr15(addr0, addr1, ps)) => { + self.tab.two_qubit_pauli_error(*addr0, *addr1, *ps) + } + + // Loss + (Loss, &QubitAndFloat(addr, p)) => self.tab.loss_channel(addr, p), + (CorrelatedLoss, TwoQubitAndFloatArr3(addr0, addr1, ps)) => { + self.tab.correlated_loss_channel(*addr0, *addr1, *ps) + } + + /* BATCH OPERATIONS START HERE */ + // Batch: dedicated batch methods + (SqrtX, QubitBatch(addrs)) => self.tab.sqrt_x_many(addrs), + (SqrtY, QubitBatch(addrs)) => self.tab.sqrt_y_many(addrs), + (SqrtXAdj, QubitBatch(addrs)) => self.tab.sqrt_x_dag_many(addrs), + (SqrtYAdj, QubitBatch(addrs)) => self.tab.sqrt_y_dag_many(addrs), + (H, QubitBatch(addrs)) => self.tab.h_many(addrs), + (CZ, TwoQubitBatch(pairs)) => self.tab.cz_many(pairs), + + // TODO: replace things below by actual batched methods once they are available + // Batch: single-qubit for loops + (X, QubitBatch(addrs)) => batch_for!(self.tab, x, addrs), + (Y, QubitBatch(addrs)) => batch_for!(self.tab, y, addrs), + (Z, QubitBatch(addrs)) => batch_for!(self.tab, z, addrs), + (S, QubitBatch(addrs)) => batch_for!(self.tab, s, addrs), + (SAdj, QubitBatch(addrs)) => batch_for!(self.tab, s_dag, addrs), + (T, QubitBatch(addrs)) => batch_for!(self.tab, t, addrs), + (TAdj, QubitBatch(addrs)) => batch_for!(self.tab, t_dag, addrs), + (Reset, QubitBatch(addrs)) => batch_for!(self.tab, reset, addrs), + (RX, QubitBatchAndFloat(addrs, angle)) => batch_for!(self.tab, rx, addrs, *angle), + (RY, QubitBatchAndFloat(addrs, angle)) => batch_for!(self.tab, ry, addrs, *angle), + (RZ, QubitBatchAndFloat(addrs, angle)) => batch_for!(self.tab, rz, addrs, *angle), + (Depolarize, QubitBatchAndFloat(addrs, p)) => { + batch_for!(self.tab, depolarize1, addrs, *p) + } + (Loss, QubitBatchAndFloat(addrs, p)) => batch_for!(self.tab, loss_channel, addrs, *p), + (PauliError, QubitBatchAndFloatArr3(addrs, ps)) => { + batch_for!(self.tab, pauli_error, addrs, *ps) + } + (U3, QubitBatchU3(addrs, theta, phi, lam)) => { + batch_for!(self.tab, u3, addrs, *theta, *phi, *lam) + } + + // Batch: two-qubit for loops + (CNOT, TwoQubitBatch(pairs)) => { + for &(a, b) in pairs { + self.tab.cnot(a, b); + } + } + (RXX, TwoQubitBatchAndFloat(pairs, angle)) => { + for &(a, b) in pairs { + self.tab.rxx(a, b, *angle); + } + } + (RYY, TwoQubitBatchAndFloat(pairs, angle)) => { + for &(a, b) in pairs { + self.tab.ryy(a, b, *angle); + } + } + (RZZ, TwoQubitBatchAndFloat(pairs, angle)) => { + for &(a, b) in pairs { + self.tab.rzz(a, b, *angle); + } + } + (Depolarize2, TwoQubitBatchAndFloat(pairs, p)) => { + for &(a, b) in pairs { + self.tab.depolarize2(a, b, *p); + } + } + (TwoQubitPauliError, TwoQubitBatchAndFloatArr15(pairs, ps)) => { + for &(a, b) in pairs { + self.tab.two_qubit_pauli_error(a, b, *ps); + } + } + (CorrelatedLoss, TwoQubitBatchAndFloatArr3(pairs, ps)) => { + for &(a, b) in pairs { + self.tab.correlated_loss_channel(a, b, *ps); + } + } + + // Batch: measure (emits per qubit) + (Measure, QubitBatch(addrs)) => { + let outcomes = addrs.iter().map(|&addr| self.tab.measure(addr).into()); + return Ok(Effects::one(CircuitOutcomeEffect::Measurement( + MeasurementEffect { + measurement_results: outcomes.collect(), + }, + ))); + } + + // Truncate is a silent no-op on the Tableau backend — the tableau's + // gate methods already prune via the configured coefficient + // threshold, so there's nothing extra to do here. + (Truncate, None) => {} + + // Trace: parse the resolved pattern and compute Σ_{P matches} ⟨ψ|P|ψ⟩ + // on the tableau state. Asymmetric with the PauliSum semantics by + // design (Decision 9): on the tableau this is a sum of expectations, + // not a coefficient filter. + (Trace, PauliPatternStr(s)) => { + let pat = PauliPattern::parse(s) + .map_err(|e| eyre!("invalid Pauli pattern `{s}`: {e:?}"))?; + let value = self.tab.trace(&pat); + return Ok(Effects::one(CircuitOutcomeEffect::Trace(TraceEffect { + value, + }))); + } + + // Fallback + (inst, msg) => { + return Err(eyre!( + "Invalid circuit instruction arguments {:?} for instruction {:?}", + msg, + inst + )); + } + }; + + Ok(Effects::None) + } +} + +impl vihaco::Reset for CircuitExecutor +where + T: Config, + <::Storage as BitView>::Store: PrimInt, + I: TableauIndex + Send + Sync + std::fmt::Debug, + C: SparseVector + std::fmt::Debug, +{ + fn reset(&mut self) { + self.tab.reset_all(); + } +} + +/// Shared dispatch body for `PauliSumExecutor` and `LossyPauliSumExecutor`. +/// Every non-loss `CircuitInstruction` lands here. `LossyPauliSumExecutor` +/// matches `Loss` / `CorrelatedLoss` (single + batched) before invoking this +/// macro and never reaches the loss-rejection arm below. +/// +/// `$self` is passed as an `ident` (typically `self`) so the macro's +/// expansion shares hygiene with the surrounding method's `self` parameter. +/// `$inst` / `$msg` are passed the same way; `$backend` is the human-readable +/// backend name baked into error messages. +macro_rules! dispatch_common_paulisum { + ($self:ident, $inst:ident, $msg:ident, $backend:literal) => {{ + use CircuitInstruction::*; + use CircuitMessage::*; + match ($inst, $msg) { + // Single-qubit Clifford + (X, &Qubit(addr)) => $self.state.x(addr), + (Y, &Qubit(addr)) => $self.state.y(addr), + (Z, &Qubit(addr)) => $self.state.z(addr), + (H, &Qubit(addr)) => $self.state.h(addr), + (S, &Qubit(addr)) => $self.state.s(addr), + (SAdj, &Qubit(addr)) => $self.state.s_dag(addr), + (SqrtX, &Qubit(addr)) => $self.state.sqrt_x(addr), + (SqrtY, &Qubit(addr)) => $self.state.sqrt_y(addr), + (SqrtXAdj, &Qubit(addr)) => $self.state.sqrt_x_dag(addr), + (SqrtYAdj, &Qubit(addr)) => $self.state.sqrt_y_dag(addr), + + // Controlled gates + (CNOT, &TwoQubit(addr0, addr1)) => $self.state.cnot(addr0, addr1), + (CZ, &TwoQubit(addr0, addr1)) => $self.state.cz(addr0, addr1), + + // Single-qubit rotations + (RX, &QubitAndFloat(addr, angle)) => $self.state.rx(addr, angle), + (RY, &QubitAndFloat(addr, angle)) => $self.state.ry(addr, angle), + (RZ, &QubitAndFloat(addr, angle)) => $self.state.rz(addr, angle), + + // Two-qubit rotations + (RXX, &TwoQubitAndFloat(addr0, addr1, angle)) => { + $self.state.rxx(addr0, addr1, angle) + } + (RYY, &TwoQubitAndFloat(addr0, addr1, angle)) => { + $self.state.ryy(addr0, addr1, angle) + } + (RZZ, &TwoQubitAndFloat(addr0, addr1, angle)) => { + $self.state.rzz(addr0, addr1, angle) + } + + // RXY: rotation about an axis in the x/y plane + (R, &QubitAndTwoFloats(addr, axis_angle, theta)) => { + $self.state.r(addr, axis_angle, theta) + } + + // Noise + (Depolarize, &QubitAndFloat(addr, p)) => $self.state.depolarize1(addr, p), + (Depolarize2, &TwoQubitAndFloat(addr0, addr1, p)) => { + $self.state.depolarize2(addr0, addr1, p) + } + (PauliError, QubitAndFloatArr3(addr0, ps)) => { + $self.state.pauli_error(*addr0, *ps) + } + (TwoQubitPauliError, TwoQubitAndFloatArr15(addr0, addr1, ps)) => { + $self.state.two_qubit_pauli_error(*addr0, *addr1, *ps) + } + + // Truncate: pruning per the configured strategy. + (Truncate, None) => $self.state.truncate(), + + // Batched arms: simple for-loop dispatch (no dedicated batch + // methods on PauliSum, unlike GeneralizedTableau). + (X, QubitBatch(addrs)) => batch_for!($self.state, x, addrs), + (Y, QubitBatch(addrs)) => batch_for!($self.state, y, addrs), + (Z, QubitBatch(addrs)) => batch_for!($self.state, z, addrs), + (H, QubitBatch(addrs)) => batch_for!($self.state, h, addrs), + (S, QubitBatch(addrs)) => batch_for!($self.state, s, addrs), + (SAdj, QubitBatch(addrs)) => batch_for!($self.state, s_dag, addrs), + (SqrtX, QubitBatch(addrs)) => batch_for!($self.state, sqrt_x, addrs), + (SqrtY, QubitBatch(addrs)) => batch_for!($self.state, sqrt_y, addrs), + (SqrtXAdj, QubitBatch(addrs)) => batch_for!($self.state, sqrt_x_dag, addrs), + (SqrtYAdj, QubitBatch(addrs)) => batch_for!($self.state, sqrt_y_dag, addrs), + (RX, QubitBatchAndFloat(addrs, angle)) => { + batch_for!($self.state, rx, addrs, *angle) + } + (RY, QubitBatchAndFloat(addrs, angle)) => { + batch_for!($self.state, ry, addrs, *angle) + } + (RZ, QubitBatchAndFloat(addrs, angle)) => { + batch_for!($self.state, rz, addrs, *angle) + } + (Depolarize, QubitBatchAndFloat(addrs, p)) => { + batch_for!($self.state, depolarize1, addrs, *p) + } + (PauliError, QubitBatchAndFloatArr3(addrs, ps)) => { + batch_for!($self.state, pauli_error, addrs, *ps) + } + (CNOT, TwoQubitBatch(pairs)) => batch_pairs_for!($self.state, cnot, pairs), + (CZ, TwoQubitBatch(pairs)) => batch_pairs_for!($self.state, cz, pairs), + (RXX, TwoQubitBatchAndFloat(pairs, angle)) => { + batch_pairs_for!($self.state, rxx, pairs, *angle) + } + (RYY, TwoQubitBatchAndFloat(pairs, angle)) => { + batch_pairs_for!($self.state, ryy, pairs, *angle) + } + (RZZ, TwoQubitBatchAndFloat(pairs, angle)) => { + batch_pairs_for!($self.state, rzz, pairs, *angle) + } + (Depolarize2, TwoQubitBatchAndFloat(pairs, p)) => { + batch_pairs_for!($self.state, depolarize2, pairs, *p) + } + (TwoQubitPauliError, TwoQubitBatchAndFloatArr15(pairs, ps)) => { + batch_pairs_for!($self.state, two_qubit_pauli_error, pairs, *ps) + } + + // Not supported on either backend (Decision 11 + Gate Support + // Matrix). Loss / CorrelatedLoss handling differs by backend + // and lives in the caller's impl block, not this macro. + (Measure | Reset, _) => { + return Err(eyre!("{} is not supported on the {} backend", $inst, $backend)); + } + + // T / T_dag / U3 are listed as supported on PauliSum in the + // plan's Gate Support Matrix, but ppvm-runtime does not yet + // implement TGate or U3Gate for PauliSum (only for + // GeneralizedTableau). + (T | TAdj | U3, _) => { + return Err(eyre!( + "{} on {} requires upstream ppvm-runtime support that is not yet implemented", + $inst, + $backend + )); + } + + // Trace: parse the resolved pattern string and compute the + // trace. Per plan Decision 9, parsing happens on every + // execution; no module-load caching. + (Trace, PauliPatternStr(s)) => { + let pat = PauliPattern::parse(s) + .map_err(|e| eyre!("invalid Pauli pattern `{}`: {:?}", s, e))?; + let value = $self.state.trace(&pat); + return Ok(Effects::one(CircuitOutcomeEffect::Trace(TraceEffect { + value, + }))); + } + + // Fallback (mismatched shapes, etc.) + (inst, msg) => { + return Err(eyre!( + "Invalid circuit instruction arguments {:?} for instruction {:?} on the {} backend", + msg, + inst, + $backend + )); + } + }; + Ok(Effects::None) + }}; +} + +/// PauliSum-backed executor (Heisenberg picture). Holds a `PauliSum` and +/// answers the same `CircuitInstruction` vocabulary as `CircuitExecutor`, +/// but without measurement / reset / loss support. +pub struct PauliSumExecutor> { + pub state: PauliSum, + /// Snapshot of the seeded observable, restored by `reset`. + initial: PauliSum, +} + +#[component(instruction = CircuitInstruction, message = CircuitMessage, effect = CircuitOutcomeEffect)] +impl PauliSumExecutor +where + T: Config, + for<'a> PauliSum: Trace<'a, PauliPattern, Output = f64>, +{ + fn execute( + &mut self, + inst: CircuitInstruction, + msg: CircuitMessage, + ) -> Result> { + self.execute_instruction(&inst, &msg) + } + + fn execute_instruction( + &mut self, + inst: &CircuitInstruction, + msg: &CircuitMessage, + ) -> Result> { + use CircuitInstruction::*; + + if matches!(inst, Loss | CorrelatedLoss) { + return Err(eyre!( + "{inst} is not supported on the PauliSum backend; use the LossyPauliSum backend instead" + )); + } + + dispatch_common_paulisum!(self, inst, msg, "PauliSum") + } +} + +impl vihaco::Reset for PauliSumExecutor +where + T: Config, + PauliSum: Clone, +{ + fn reset(&mut self) { + self.state = self.initial.clone(); + } +} + +/// LossyPauliSum-backed executor. Same dispatch as `PauliSumExecutor` plus +/// `Loss` / `CorrelatedLoss` channels. The concrete `T` used by the +/// enclosing `Circuit::LossyPauliSum` variant is a `Config` whose +/// `PauliWordType` is `LossyPauliWord` (see `LossyPauliSumConfig`). +pub struct LossyPauliSumExecutor> { + pub state: PauliSum, + /// Snapshot of the seeded observable, restored by `reset`. + initial: PauliSum, +} + +#[component(instruction = CircuitInstruction, message = CircuitMessage, effect = CircuitOutcomeEffect)] +impl LossyPauliSumExecutor +where + T: Config, + for<'a> PauliSum: Trace<'a, PauliPattern, Output = f64>, +{ + fn execute( + &mut self, + inst: CircuitInstruction, + msg: CircuitMessage, + ) -> Result> { + self.execute_instruction(&inst, &msg) + } + + fn execute_instruction( + &mut self, + inst: &CircuitInstruction, + msg: &CircuitMessage, + ) -> Result> { + use CircuitInstruction::*; + use CircuitMessage::*; + + // Loss / CorrelatedLoss are the only instructions that differ from + // PauliSum; handle them here then delegate everything else to the + // shared dispatch. + match (inst, msg) { + (Loss, &QubitAndFloat(addr, p)) => { + self.state.loss_channel(addr, p); + return Ok(Effects::None); + } + (CorrelatedLoss, TwoQubitAndFloatArr3(addr0, addr1, ps)) => { + self.state.correlated_loss_channel(*addr0, *addr1, *ps); + return Ok(Effects::None); + } + (Loss, QubitBatchAndFloat(addrs, p)) => { + batch_for!(self.state, loss_channel, addrs, *p); + return Ok(Effects::None); + } + (CorrelatedLoss, TwoQubitBatchAndFloatArr3(pairs, ps)) => { + batch_pairs_for!(self.state, correlated_loss_channel, pairs, *ps); + return Ok(Effects::None); + } + _ => {} + } + + dispatch_common_paulisum!(self, inst, msg, "LossyPauliSum") + } +} + +impl vihaco::Reset for LossyPauliSumExecutor +where + T: Config, + PauliSum: Clone, +{ + fn reset(&mut self) { + self.state = self.initial.clone(); + } +} + +/// Tableau-backed inner enum (Schrödinger picture). Carries the six +/// size-bucketed `CircuitExecutor` variants; bucket is picked from `n_qubits`. +pub enum TableauCircuit { + Bits64(CircuitExecutor, usize, Vec<(Complex64, usize)>>), + Bits128(CircuitExecutor, u128, Vec<(Complex64, u128)>>), + Bits256(CircuitExecutor, U256, Vec<(Complex64, U256)>>), + Bits512(CircuitExecutor, U512, Vec<(Complex64, U512)>>), + Bits1024(CircuitExecutor, U1024, Vec<(Complex64, U1024)>>), + Bits2048(CircuitExecutor, U2048, Vec<(Complex64, U2048)>>), +} + +impl TableauCircuit { + pub fn new(n_qubits: usize, coefficient_threshold: f64) -> Self { + if n_qubits <= 64 { + let tab = GeneralizedTableau::new(n_qubits, coefficient_threshold); + Self::Bits64(CircuitExecutor { tab }) + } else if n_qubits <= 128 { + let tab = GeneralizedTableau::new(n_qubits, coefficient_threshold); + Self::Bits128(CircuitExecutor { tab }) + } else if n_qubits <= 256 { + let tab = GeneralizedTableau::new(n_qubits, coefficient_threshold); + Self::Bits256(CircuitExecutor { tab }) + } else if n_qubits <= 512 { + let tab = GeneralizedTableau::new(n_qubits, coefficient_threshold); + Self::Bits512(CircuitExecutor { tab }) + } else if n_qubits <= 1024 { + let tab = GeneralizedTableau::new(n_qubits, coefficient_threshold); + Self::Bits1024(CircuitExecutor { tab }) + } else if n_qubits <= 2048 { + let tab = GeneralizedTableau::new(n_qubits, coefficient_threshold); + Self::Bits2048(CircuitExecutor { tab }) + } else { + panic!("No matching executor for {} qubits", n_qubits); + } + } + + /// Same as [`TableauCircuit::new`], but seed the RNG deterministically so a + /// shot is reproducible. + pub fn new_with_seed(n_qubits: usize, coefficient_threshold: f64, seed: u64) -> Self { + macro_rules! seeded { + ($variant:ident) => {{ + let tab = GeneralizedTableau::new_with_seed(n_qubits, coefficient_threshold, seed); + Self::$variant(CircuitExecutor { tab }) + }}; + } + if n_qubits <= 64 { + seeded!(Bits64) + } else if n_qubits <= 128 { + seeded!(Bits128) + } else if n_qubits <= 256 { + seeded!(Bits256) + } else if n_qubits <= 512 { + seeded!(Bits512) + } else if n_qubits <= 1024 { + seeded!(Bits1024) + } else if n_qubits <= 2048 { + seeded!(Bits2048) + } else { + panic!("No matching executor for {} qubits", n_qubits); + } + } + + fn execute_instruction( + &mut self, + inst: &CircuitInstruction, + msg: &CircuitMessage, + ) -> Result> { + match self { + Self::Bits64(ex) => ex.execute_instruction(inst, msg), + Self::Bits128(ex) => ex.execute_instruction(inst, msg), + Self::Bits256(ex) => ex.execute_instruction(inst, msg), + Self::Bits512(ex) => ex.execute_instruction(inst, msg), + Self::Bits1024(ex) => ex.execute_instruction(inst, msg), + Self::Bits2048(ex) => ex.execute_instruction(inst, msg), + } + } + + pub fn state_string(&self) -> String { + match self { + Self::Bits64(ex) => ex.tab.to_string(), + Self::Bits128(ex) => ex.tab.to_string(), + Self::Bits256(ex) => ex.tab.to_string(), + Self::Bits512(ex) => ex.tab.to_string(), + Self::Bits1024(ex) => ex.tab.to_string(), + Self::Bits2048(ex) => ex.tab.to_string(), + } + } +} + +impl vihaco::Reset for TableauCircuit { + fn reset(&mut self) { + match self { + Self::Bits64(ex) => ex.reset(), + Self::Bits128(ex) => ex.reset(), + Self::Bits256(ex) => ex.reset(), + Self::Bits512(ex) => ex.reset(), + Self::Bits1024(ex) => ex.reset(), + Self::Bits2048(ex) => ex.reset(), + }; + } +} + +/// PauliSum-backed inner enum (Heisenberg picture). Per Decision 7 of the plan, +/// the size buckets carry `[u8; N]`-storage `ByteFxHashF64` configs (N = 8, 16, +/// …, 256) rather than the tableau's `[u64; N]` configs; bucket labels match +/// the semantic qubit count (`Bits64` = 64 qubits) so the outer enum's dispatch +/// is uniform across backends. +pub enum PauliSumCircuit { + Bits64(PauliSumExecutor>), + Bits128(PauliSumExecutor>), + Bits256(PauliSumExecutor>), + Bits512(PauliSumExecutor>), + Bits1024(PauliSumExecutor>), + Bits2048(PauliSumExecutor>), +} + +impl PauliSumCircuit { + /// Build a PauliSum-backed circuit, seeding the state with every term: + /// `for (word, coef) in terms { state += (word, coef); }`. Words must + /// already be validated against `info.n_qubits` by the caller. + pub fn new(info: &PPVMDeviceInfo, terms: &[(String, f64)]) -> Self { + macro_rules! build { + ($variant:ident, $N:literal) => {{ + let mut state = PauliSum::>::builder() + .n_qubits(info.n_qubits) + .strategy(paulisum_strategy(info)) + .build(); + for (word, coef) in terms { + state += (word.as_str(), *coef); + } + let initial = state.clone(); + Self::$variant(PauliSumExecutor { state, initial }) + }}; + } + if info.n_qubits <= 64 { + build!(Bits64, 8) + } else if info.n_qubits <= 128 { + build!(Bits128, 16) + } else if info.n_qubits <= 256 { + build!(Bits256, 32) + } else if info.n_qubits <= 512 { + build!(Bits512, 64) + } else if info.n_qubits <= 1024 { + build!(Bits1024, 128) + } else if info.n_qubits <= 2048 { + build!(Bits2048, 256) + } else { + panic!("No matching PauliSum executor for {} qubits", info.n_qubits); + } + } + + fn execute_instruction( + &mut self, + inst: &CircuitInstruction, + msg: &CircuitMessage, + ) -> Result> { + match self { + Self::Bits64(ex) => ex.execute_instruction(inst, msg), + Self::Bits128(ex) => ex.execute_instruction(inst, msg), + Self::Bits256(ex) => ex.execute_instruction(inst, msg), + Self::Bits512(ex) => ex.execute_instruction(inst, msg), + Self::Bits1024(ex) => ex.execute_instruction(inst, msg), + Self::Bits2048(ex) => ex.execute_instruction(inst, msg), + } + } + + pub fn state_string(&self) -> String { + match self { + Self::Bits64(ex) => ex.state.to_string(), + Self::Bits128(ex) => ex.state.to_string(), + Self::Bits256(ex) => ex.state.to_string(), + Self::Bits512(ex) => ex.state.to_string(), + Self::Bits1024(ex) => ex.state.to_string(), + Self::Bits2048(ex) => ex.state.to_string(), + } + } +} + +impl vihaco::Reset for PauliSumCircuit { + fn reset(&mut self) { + match self { + Self::Bits64(ex) => ex.reset(), + Self::Bits128(ex) => ex.reset(), + Self::Bits256(ex) => ex.reset(), + Self::Bits512(ex) => ex.reset(), + Self::Bits1024(ex) => ex.reset(), + Self::Bits2048(ex) => ex.reset(), + }; + } +} + +/// LossyPauliSum-backed inner enum. Identical shape to [`PauliSumCircuit`] +/// but with `LossyPauliWord`-keyed configs so loss-channel methods dispatch. +pub enum LossyPauliSumCircuit { + Bits64(LossyPauliSumExecutor>), + Bits128(LossyPauliSumExecutor>), + Bits256(LossyPauliSumExecutor>), + Bits512(LossyPauliSumExecutor>), + Bits1024(LossyPauliSumExecutor>), + Bits2048(LossyPauliSumExecutor>), +} + +impl LossyPauliSumCircuit { + /// Build a LossyPauliSum-backed circuit, seeding every term via + /// `state += (word, coef)`. Words must already be validated against + /// `info.n_qubits` by the caller. + pub fn new(info: &PPVMDeviceInfo, terms: &[(String, f64)]) -> Self { + macro_rules! build { + ($variant:ident, $N:literal) => {{ + let mut state = PauliSum::>::builder() + .n_qubits(info.n_qubits) + .strategy(paulisum_strategy(info)) + .build(); + for (word, coef) in terms { + state += (word.as_str(), *coef); + } + let initial = state.clone(); + Self::$variant(LossyPauliSumExecutor { state, initial }) + }}; + } + if info.n_qubits <= 64 { + build!(Bits64, 8) + } else if info.n_qubits <= 128 { + build!(Bits128, 16) + } else if info.n_qubits <= 256 { + build!(Bits256, 32) + } else if info.n_qubits <= 512 { + build!(Bits512, 64) + } else if info.n_qubits <= 1024 { + build!(Bits1024, 128) + } else if info.n_qubits <= 2048 { + build!(Bits2048, 256) + } else { + panic!( + "No matching LossyPauliSum executor for {} qubits", + info.n_qubits + ); + } + } + + fn execute_instruction( + &mut self, + inst: &CircuitInstruction, + msg: &CircuitMessage, + ) -> Result> { + match self { + Self::Bits64(ex) => ex.execute_instruction(inst, msg), + Self::Bits128(ex) => ex.execute_instruction(inst, msg), + Self::Bits256(ex) => ex.execute_instruction(inst, msg), + Self::Bits512(ex) => ex.execute_instruction(inst, msg), + Self::Bits1024(ex) => ex.execute_instruction(inst, msg), + Self::Bits2048(ex) => ex.execute_instruction(inst, msg), + } + } + + pub fn state_string(&self) -> String { + match self { + Self::Bits64(ex) => ex.state.to_string(), + Self::Bits128(ex) => ex.state.to_string(), + Self::Bits256(ex) => ex.state.to_string(), + Self::Bits512(ex) => ex.state.to_string(), + Self::Bits1024(ex) => ex.state.to_string(), + Self::Bits2048(ex) => ex.state.to_string(), + } + } +} + +impl vihaco::Reset for LossyPauliSumCircuit { + fn reset(&mut self) { + match self { + Self::Bits64(ex) => ex.reset(), + Self::Bits128(ex) => ex.reset(), + Self::Bits256(ex) => ex.reset(), + Self::Bits512(ex) => ex.reset(), + Self::Bits1024(ex) => ex.reset(), + Self::Bits2048(ex) => ex.reset(), + }; + } +} + +/// Outer `Circuit` enum: backend selector. Picks one of the three inner enums +/// based on `info.backend` at construction time; from there, every per-step +/// call routes outer → inner → executor. +pub enum Circuit { + Tableau(TableauCircuit), + PauliSum(PauliSumCircuit), + LossyPauliSum(LossyPauliSumCircuit), +} + +#[component(instruction = CircuitInstruction, message = CircuitMessage, effect = CircuitOutcomeEffect)] +impl Circuit { + /// Build a Tableau-backed circuit. Tableau init only needs `n_qubits` and + /// `coefficient_threshold` from `info`; no observable required. + pub fn tableau(info: &PPVMDeviceInfo) -> Self { + Self::Tableau(TableauCircuit::new( + info.n_qubits, + info.coefficient_threshold, + )) + } + + /// Same as [`Circuit::tableau`], but seed the tableau RNG deterministically + /// so a shot is reproducible. + pub fn tableau_with_seed(info: &PPVMDeviceInfo, seed: u64) -> Self { + Self::Tableau(TableauCircuit::new_with_seed( + info.n_qubits, + info.coefficient_threshold, + seed, + )) + } + + /// Build a PauliSum-backed circuit, seeding the state with every term in + /// `terms`. Each `(word, coef)` is added via `state += (word, coef)`; the + /// caller is responsible for having parsed/validated the words against + /// `info.n_qubits` (see `parse_observable_terms` in `composite.rs`). + pub fn paulisum(info: &PPVMDeviceInfo, terms: &[(String, f64)]) -> Self { + Self::PauliSum(PauliSumCircuit::new(info, terms)) + } + + /// Build a LossyPauliSum-backed circuit. Same contract as + /// [`Circuit::paulisum`]. + pub fn lossy_paulisum(info: &PPVMDeviceInfo, terms: &[(String, f64)]) -> Self { + Self::LossyPauliSum(LossyPauliSumCircuit::new(info, terms)) + } + + fn execute( + &mut self, + inst: CircuitInstruction, + msg: CircuitMessage, + ) -> Result> { + self.execute_instruction(&inst, &msg) + } + + fn execute_instruction( + &mut self, + inst: &CircuitInstruction, + msg: &CircuitMessage, + ) -> Result> { + match self { + Self::Tableau(c) => c.execute_instruction(inst, msg), + Self::PauliSum(c) => c.execute_instruction(inst, msg), + Self::LossyPauliSum(c) => c.execute_instruction(inst, msg), + } + } + + /// Render the current state. Used by the REPL's `show` command. + pub fn state_string(&self) -> String { + match self { + Self::Tableau(c) => c.state_string(), + Self::PauliSum(c) => c.state_string(), + Self::LossyPauliSum(c) => c.state_string(), + } + } +} + +#[observe(CircuitEffect, effect=CircuitOutcomeEffect)] +impl Circuit { + fn observe_circuit_effect( + &mut self, + effect: &CircuitEffect, + ) -> Result> { + self.execute_instruction(&effect.inst, &effect.msg) + } +} + +impl vihaco::Reset for Circuit { + fn reset(&mut self) { + match self { + Self::Tableau(c) => c.reset(), + Self::PauliSum(c) => c.reset(), + Self::LossyPauliSum(c) => c.reset(), + }; + } +} + +impl Default for Circuit { + fn default() -> Self { + // Default backend is Tableau, which doesn't require an observable. + Self::tableau(&PPVMDeviceInfo::default()) + } +} diff --git a/crates/ppvm-vihaco/src/measurements.rs b/crates/ppvm-vihaco/src/measurements.rs new file mode 100644 index 000000000..1881c0d51 --- /dev/null +++ b/crates/ppvm-vihaco/src/measurements.rs @@ -0,0 +1,93 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +use eyre::Result; +use smallvec::SmallVec; +use vihaco::{Effects, observe}; + +/// Measurement results are represented as an integer enum: +/// 0: state |0> +/// 1: state |1> +/// 2: qubit has been lost prior to measurement +/// In bytecode, this is represented as a u32 integer, which is simpler than +/// e.g. two boolean values and matches semantics elsewhere +#[repr(u8)] +#[derive(Debug, Clone, Copy, Eq, PartialEq)] +pub enum MeasurementOutcome { + Zero = 0, + One = 1, + Lost = 2, +} + +pub type MeasurementResult = SmallVec<[MeasurementOutcome; 8]>; + +#[derive(Debug, Clone, Eq, PartialEq)] +pub struct MeasurementEffect { + pub measurement_results: MeasurementResult, +} + +impl From> for MeasurementOutcome { + fn from(m: Option) -> Self { + match m { + Some(false) => Self::Zero, + Some(true) => Self::One, + None => Self::Lost, + } + } +} + +#[derive(Debug, Default)] +pub struct MeasurementObserver { + pub record: Vec, +} + +#[observe(MeasurementEffect)] +impl MeasurementObserver { + fn observe_measurement_effect(&mut self, effect: &MeasurementEffect) -> Result> { + self.record.push(effect.measurement_results.clone()); + Ok(Effects::none()) + } +} + +/// Per-step trace value emitted by the `Trace` instruction on the PauliSum and +/// LossyPauliSum backends. The plan's Decision 5 keeps trace and measurement +/// records as two parallel streams; this effect feeds the trace stream. +#[derive(Debug, Clone, PartialEq)] +pub struct TraceEffect { + pub value: f64, +} + +#[derive(Debug, Default)] +pub struct TraceObserver { + pub record: Vec, +} + +#[observe(TraceEffect)] +impl TraceObserver { + fn observe_trace_effect(&mut self, effect: &TraceEffect) -> Result> { + self.record.push(effect.value); + Ok(Effects::none()) + } +} + +/// Union of the two effect types a circuit instruction can produce. The plan's +/// structural note (Task 6) calls for broadening the `#[component(..., effect = +/// ...)]` annotation on `Circuit` and the executors to a union so a single +/// `Trace` (or `Measure`) instruction can fan out to the right observer. +#[derive(Debug, Clone, PartialEq)] +pub enum CircuitOutcomeEffect { + Measurement(MeasurementEffect), + Trace(TraceEffect), +} + +impl From for CircuitOutcomeEffect { + fn from(value: MeasurementEffect) -> Self { + Self::Measurement(value) + } +} + +impl From for CircuitOutcomeEffect { + fn from(value: TraceEffect) -> Self { + Self::Trace(value) + } +} From c876dc0e52dc16f84644414c438cbbc392e1f050 Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Wed, 8 Jul 2026 11:26:41 +0200 Subject: [PATCH 3/6] fix: remove stray ppvm-vihaco files committed by mistake `component.rs` and `measurements.rs` were staged for the circuit-component PR (#4) and got committed into d43fda1c by accident when the index was carried across a branch switch. They belong on `david/42.2-circuit-component`, not this PR. They were never wired into the workspace so they had no build effect; removing them keeps #172 scoped to the tableau changes. Co-Authored-By: Claude Opus 4.8 (1M context) --- crates/ppvm-vihaco/src/component.rs | 935 ------------------------- crates/ppvm-vihaco/src/measurements.rs | 93 --- 2 files changed, 1028 deletions(-) delete mode 100644 crates/ppvm-vihaco/src/component.rs delete mode 100644 crates/ppvm-vihaco/src/measurements.rs diff --git a/crates/ppvm-vihaco/src/component.rs b/crates/ppvm-vihaco/src/component.rs deleted file mode 100644 index ae9e70806..000000000 --- a/crates/ppvm-vihaco/src/component.rs +++ /dev/null @@ -1,935 +0,0 @@ -// SPDX-FileCopyrightText: 2026 The PPVM Authors -// SPDX-License-Identifier: Apache-2.0 - -use crate::composite::PPVMDeviceInfo; -use crate::measurements::{ - CircuitOutcomeEffect, MeasurementEffect, MeasurementOutcome, TraceEffect, -}; -use bitvec::view::BitView; -use bnum::types::{U256, U512, U1024, U2048}; -use eyre::{Result, eyre}; -use num::PrimInt; -use num::complex::Complex64; -use ppvm_pauli_sum::config::fx64hash::Byte8F64; -use ppvm_pauli_sum::config::indexmap::ByteFxHashF64; -use ppvm_pauli_sum::strategy::{CoefficientThreshold, CombinedStrategy, MaxPauliWeight}; -use ppvm_tableau::prelude::*; -use vihaco::{Effects, component, observe}; -use vihaco_circuit_isa::{CircuitEffect, CircuitInstruction, CircuitMessage}; - -/// Truncation strategy used by every `PauliSum` / `LossyPauliSum` size bucket. -/// Coefficient-threshold pruning is always on; the Pauli-weight cap is set per -/// run from the header (defaults to `usize::MAX` = no cap). -type PauliSumStrategy = CombinedStrategy; - -/// `PauliSum`'s `T` for the lossless backend: `[u8; N]` storage, fx hash, -/// f64 coefficients, the strategy above. -type PauliSumConfig = ByteFxHashF64; - -/// Same as `PauliSumConfig` but with `LossyPauliWord` as the word type, so the -/// loss-channel methods are dispatchable on the resulting `PauliSum`. -/// `LossyPauliWord`'s second type parameter (hasher) defaults to -/// `fxhash::FxBuildHasher`, matching `ByteFxHashF64`'s internal hasher. -type LossyPauliSumConfig = - ByteFxHashF64>; - -/// Build a `PauliSumStrategy` value from a `PPVMDeviceInfo`. Pulled out so the -/// six size-bucket constructors don't each repeat the strategy spelling. -fn paulisum_strategy(info: &PPVMDeviceInfo) -> PauliSumStrategy { - CombinedStrategy( - CoefficientThreshold(info.coefficient_threshold), - MaxPauliWeight(info.max_pauli_weight.unwrap_or(usize::MAX)), - ) -} - -macro_rules! batch_for { - ($tab:expr, $method:ident, $addrs:expr) => { - for addr in $addrs { $tab.$method(*addr); } - }; - ($tab:expr, $method:ident, $addrs:expr, $($arg:expr),+) => { - for addr in $addrs { $tab.$method(*addr, $($arg),+); } - }; -} - -/// Two-qubit sibling of [`batch_for!`]: drives a method that takes two qubit -/// addresses (plus optional extra args) over a slice of `(usize, usize)` pairs. -macro_rules! batch_pairs_for { - ($state:expr, $method:ident, $pairs:expr) => { - for &(a, b) in $pairs { $state.$method(a, b); } - }; - ($state:expr, $method:ident, $pairs:expr, $($arg:expr),+) => { - for &(a, b) in $pairs { $state.$method(a, b, $($arg),+); } - }; -} - -pub struct CircuitExecutor, I: TableauIndex, C: SparseVector> { - pub tab: GeneralizedTableau, -} - -#[component(instruction = CircuitInstruction, message = CircuitMessage, effect = CircuitOutcomeEffect)] -impl CircuitExecutor -where - T: Config, - <::Storage as BitView>::Store: PrimInt, - I: TableauIndex + Send + Sync + std::fmt::Debug, - C: SparseVector + std::fmt::Debug, -{ - fn execute( - &mut self, - inst: CircuitInstruction, - msg: CircuitMessage, - ) -> Result> { - self.execute_instruction(&inst, &msg) - } - - fn execute_instruction( - &mut self, - inst: &CircuitInstruction, - msg: &CircuitMessage, - ) -> Result> { - use CircuitInstruction::*; - use CircuitMessage::*; - - match (inst, msg) { - // Single-qubit Clifford - (X, &Qubit(addr)) => self.tab.x(addr), - (Y, &Qubit(addr)) => self.tab.y(addr), - (Z, &Qubit(addr)) => self.tab.z(addr), - (H, &Qubit(addr)) => self.tab.h(addr), - (S, &Qubit(addr)) => self.tab.s(addr), - (SAdj, &Qubit(addr)) => self.tab.s_dag(addr), - (SqrtX, &Qubit(addr)) => self.tab.sqrt_x(addr), - (SqrtY, &Qubit(addr)) => self.tab.sqrt_y(addr), - (SqrtXAdj, &Qubit(addr)) => self.tab.sqrt_x_dag(addr), - (SqrtYAdj, &Qubit(addr)) => self.tab.sqrt_y_dag(addr), - - // Controlled gates - (CNOT, &TwoQubit(addr0, addr1)) => self.tab.cnot(addr0, addr1), - (CZ, &TwoQubit(addr0, addr1)) => self.tab.cz(addr0, addr1), - - // T gate - (T, &Qubit(addr)) => self.tab.t(addr), - (TAdj, &Qubit(addr)) => self.tab.t_dag(addr), - - // Single-qubit rotations - (RX, &QubitAndFloat(addr, angle)) => self.tab.rx(addr, angle), - (RY, &QubitAndFloat(addr, angle)) => self.tab.ry(addr, angle), - (RZ, &QubitAndFloat(addr, angle)) => self.tab.rz(addr, angle), - - // Two-qubit rotations - (RXX, &TwoQubitAndFloat(addr0, addr1, angle)) => self.tab.rxx(addr0, addr1, angle), - (RYY, &TwoQubitAndFloat(addr0, addr1, angle)) => self.tab.ryy(addr0, addr1, angle), - (RZZ, &TwoQubitAndFloat(addr0, addr1, angle)) => self.tab.rzz(addr0, addr1, angle), - - // U3 - (U3, &QubitU3(addr, theta, phi, lam)) => self.tab.u3(addr, theta, phi, lam), - - // RXY: rotation about an axis in the x/y plane - (R, &QubitAndTwoFloats(addr, axis_angle, theta)) => self.tab.r(addr, axis_angle, theta), - - // Measure & Reset - (Measure, &Qubit(addr)) => { - let outcome: MeasurementOutcome = self.tab.measure(addr).into(); - return Ok(Effects::one(CircuitOutcomeEffect::Measurement( - MeasurementEffect { - measurement_results: smallvec::smallvec![outcome], - }, - ))); - } - (Reset, &Qubit(addr)) => self.tab.reset(addr), - - // Noise - (Depolarize, &QubitAndFloat(addr, p)) => self.tab.depolarize1(addr, p), - (Depolarize2, &TwoQubitAndFloat(addr0, addr1, p)) => { - self.tab.depolarize2(addr0, addr1, p) - } - (PauliError, QubitAndFloatArr3(addr0, ps)) => self.tab.pauli_error(*addr0, *ps), - (TwoQubitPauliError, TwoQubitAndFloatArr15(addr0, addr1, ps)) => { - self.tab.two_qubit_pauli_error(*addr0, *addr1, *ps) - } - - // Loss - (Loss, &QubitAndFloat(addr, p)) => self.tab.loss_channel(addr, p), - (CorrelatedLoss, TwoQubitAndFloatArr3(addr0, addr1, ps)) => { - self.tab.correlated_loss_channel(*addr0, *addr1, *ps) - } - - /* BATCH OPERATIONS START HERE */ - // Batch: dedicated batch methods - (SqrtX, QubitBatch(addrs)) => self.tab.sqrt_x_many(addrs), - (SqrtY, QubitBatch(addrs)) => self.tab.sqrt_y_many(addrs), - (SqrtXAdj, QubitBatch(addrs)) => self.tab.sqrt_x_dag_many(addrs), - (SqrtYAdj, QubitBatch(addrs)) => self.tab.sqrt_y_dag_many(addrs), - (H, QubitBatch(addrs)) => self.tab.h_many(addrs), - (CZ, TwoQubitBatch(pairs)) => self.tab.cz_many(pairs), - - // TODO: replace things below by actual batched methods once they are available - // Batch: single-qubit for loops - (X, QubitBatch(addrs)) => batch_for!(self.tab, x, addrs), - (Y, QubitBatch(addrs)) => batch_for!(self.tab, y, addrs), - (Z, QubitBatch(addrs)) => batch_for!(self.tab, z, addrs), - (S, QubitBatch(addrs)) => batch_for!(self.tab, s, addrs), - (SAdj, QubitBatch(addrs)) => batch_for!(self.tab, s_dag, addrs), - (T, QubitBatch(addrs)) => batch_for!(self.tab, t, addrs), - (TAdj, QubitBatch(addrs)) => batch_for!(self.tab, t_dag, addrs), - (Reset, QubitBatch(addrs)) => batch_for!(self.tab, reset, addrs), - (RX, QubitBatchAndFloat(addrs, angle)) => batch_for!(self.tab, rx, addrs, *angle), - (RY, QubitBatchAndFloat(addrs, angle)) => batch_for!(self.tab, ry, addrs, *angle), - (RZ, QubitBatchAndFloat(addrs, angle)) => batch_for!(self.tab, rz, addrs, *angle), - (Depolarize, QubitBatchAndFloat(addrs, p)) => { - batch_for!(self.tab, depolarize1, addrs, *p) - } - (Loss, QubitBatchAndFloat(addrs, p)) => batch_for!(self.tab, loss_channel, addrs, *p), - (PauliError, QubitBatchAndFloatArr3(addrs, ps)) => { - batch_for!(self.tab, pauli_error, addrs, *ps) - } - (U3, QubitBatchU3(addrs, theta, phi, lam)) => { - batch_for!(self.tab, u3, addrs, *theta, *phi, *lam) - } - - // Batch: two-qubit for loops - (CNOT, TwoQubitBatch(pairs)) => { - for &(a, b) in pairs { - self.tab.cnot(a, b); - } - } - (RXX, TwoQubitBatchAndFloat(pairs, angle)) => { - for &(a, b) in pairs { - self.tab.rxx(a, b, *angle); - } - } - (RYY, TwoQubitBatchAndFloat(pairs, angle)) => { - for &(a, b) in pairs { - self.tab.ryy(a, b, *angle); - } - } - (RZZ, TwoQubitBatchAndFloat(pairs, angle)) => { - for &(a, b) in pairs { - self.tab.rzz(a, b, *angle); - } - } - (Depolarize2, TwoQubitBatchAndFloat(pairs, p)) => { - for &(a, b) in pairs { - self.tab.depolarize2(a, b, *p); - } - } - (TwoQubitPauliError, TwoQubitBatchAndFloatArr15(pairs, ps)) => { - for &(a, b) in pairs { - self.tab.two_qubit_pauli_error(a, b, *ps); - } - } - (CorrelatedLoss, TwoQubitBatchAndFloatArr3(pairs, ps)) => { - for &(a, b) in pairs { - self.tab.correlated_loss_channel(a, b, *ps); - } - } - - // Batch: measure (emits per qubit) - (Measure, QubitBatch(addrs)) => { - let outcomes = addrs.iter().map(|&addr| self.tab.measure(addr).into()); - return Ok(Effects::one(CircuitOutcomeEffect::Measurement( - MeasurementEffect { - measurement_results: outcomes.collect(), - }, - ))); - } - - // Truncate is a silent no-op on the Tableau backend — the tableau's - // gate methods already prune via the configured coefficient - // threshold, so there's nothing extra to do here. - (Truncate, None) => {} - - // Trace: parse the resolved pattern and compute Σ_{P matches} ⟨ψ|P|ψ⟩ - // on the tableau state. Asymmetric with the PauliSum semantics by - // design (Decision 9): on the tableau this is a sum of expectations, - // not a coefficient filter. - (Trace, PauliPatternStr(s)) => { - let pat = PauliPattern::parse(s) - .map_err(|e| eyre!("invalid Pauli pattern `{s}`: {e:?}"))?; - let value = self.tab.trace(&pat); - return Ok(Effects::one(CircuitOutcomeEffect::Trace(TraceEffect { - value, - }))); - } - - // Fallback - (inst, msg) => { - return Err(eyre!( - "Invalid circuit instruction arguments {:?} for instruction {:?}", - msg, - inst - )); - } - }; - - Ok(Effects::None) - } -} - -impl vihaco::Reset for CircuitExecutor -where - T: Config, - <::Storage as BitView>::Store: PrimInt, - I: TableauIndex + Send + Sync + std::fmt::Debug, - C: SparseVector + std::fmt::Debug, -{ - fn reset(&mut self) { - self.tab.reset_all(); - } -} - -/// Shared dispatch body for `PauliSumExecutor` and `LossyPauliSumExecutor`. -/// Every non-loss `CircuitInstruction` lands here. `LossyPauliSumExecutor` -/// matches `Loss` / `CorrelatedLoss` (single + batched) before invoking this -/// macro and never reaches the loss-rejection arm below. -/// -/// `$self` is passed as an `ident` (typically `self`) so the macro's -/// expansion shares hygiene with the surrounding method's `self` parameter. -/// `$inst` / `$msg` are passed the same way; `$backend` is the human-readable -/// backend name baked into error messages. -macro_rules! dispatch_common_paulisum { - ($self:ident, $inst:ident, $msg:ident, $backend:literal) => {{ - use CircuitInstruction::*; - use CircuitMessage::*; - match ($inst, $msg) { - // Single-qubit Clifford - (X, &Qubit(addr)) => $self.state.x(addr), - (Y, &Qubit(addr)) => $self.state.y(addr), - (Z, &Qubit(addr)) => $self.state.z(addr), - (H, &Qubit(addr)) => $self.state.h(addr), - (S, &Qubit(addr)) => $self.state.s(addr), - (SAdj, &Qubit(addr)) => $self.state.s_dag(addr), - (SqrtX, &Qubit(addr)) => $self.state.sqrt_x(addr), - (SqrtY, &Qubit(addr)) => $self.state.sqrt_y(addr), - (SqrtXAdj, &Qubit(addr)) => $self.state.sqrt_x_dag(addr), - (SqrtYAdj, &Qubit(addr)) => $self.state.sqrt_y_dag(addr), - - // Controlled gates - (CNOT, &TwoQubit(addr0, addr1)) => $self.state.cnot(addr0, addr1), - (CZ, &TwoQubit(addr0, addr1)) => $self.state.cz(addr0, addr1), - - // Single-qubit rotations - (RX, &QubitAndFloat(addr, angle)) => $self.state.rx(addr, angle), - (RY, &QubitAndFloat(addr, angle)) => $self.state.ry(addr, angle), - (RZ, &QubitAndFloat(addr, angle)) => $self.state.rz(addr, angle), - - // Two-qubit rotations - (RXX, &TwoQubitAndFloat(addr0, addr1, angle)) => { - $self.state.rxx(addr0, addr1, angle) - } - (RYY, &TwoQubitAndFloat(addr0, addr1, angle)) => { - $self.state.ryy(addr0, addr1, angle) - } - (RZZ, &TwoQubitAndFloat(addr0, addr1, angle)) => { - $self.state.rzz(addr0, addr1, angle) - } - - // RXY: rotation about an axis in the x/y plane - (R, &QubitAndTwoFloats(addr, axis_angle, theta)) => { - $self.state.r(addr, axis_angle, theta) - } - - // Noise - (Depolarize, &QubitAndFloat(addr, p)) => $self.state.depolarize1(addr, p), - (Depolarize2, &TwoQubitAndFloat(addr0, addr1, p)) => { - $self.state.depolarize2(addr0, addr1, p) - } - (PauliError, QubitAndFloatArr3(addr0, ps)) => { - $self.state.pauli_error(*addr0, *ps) - } - (TwoQubitPauliError, TwoQubitAndFloatArr15(addr0, addr1, ps)) => { - $self.state.two_qubit_pauli_error(*addr0, *addr1, *ps) - } - - // Truncate: pruning per the configured strategy. - (Truncate, None) => $self.state.truncate(), - - // Batched arms: simple for-loop dispatch (no dedicated batch - // methods on PauliSum, unlike GeneralizedTableau). - (X, QubitBatch(addrs)) => batch_for!($self.state, x, addrs), - (Y, QubitBatch(addrs)) => batch_for!($self.state, y, addrs), - (Z, QubitBatch(addrs)) => batch_for!($self.state, z, addrs), - (H, QubitBatch(addrs)) => batch_for!($self.state, h, addrs), - (S, QubitBatch(addrs)) => batch_for!($self.state, s, addrs), - (SAdj, QubitBatch(addrs)) => batch_for!($self.state, s_dag, addrs), - (SqrtX, QubitBatch(addrs)) => batch_for!($self.state, sqrt_x, addrs), - (SqrtY, QubitBatch(addrs)) => batch_for!($self.state, sqrt_y, addrs), - (SqrtXAdj, QubitBatch(addrs)) => batch_for!($self.state, sqrt_x_dag, addrs), - (SqrtYAdj, QubitBatch(addrs)) => batch_for!($self.state, sqrt_y_dag, addrs), - (RX, QubitBatchAndFloat(addrs, angle)) => { - batch_for!($self.state, rx, addrs, *angle) - } - (RY, QubitBatchAndFloat(addrs, angle)) => { - batch_for!($self.state, ry, addrs, *angle) - } - (RZ, QubitBatchAndFloat(addrs, angle)) => { - batch_for!($self.state, rz, addrs, *angle) - } - (Depolarize, QubitBatchAndFloat(addrs, p)) => { - batch_for!($self.state, depolarize1, addrs, *p) - } - (PauliError, QubitBatchAndFloatArr3(addrs, ps)) => { - batch_for!($self.state, pauli_error, addrs, *ps) - } - (CNOT, TwoQubitBatch(pairs)) => batch_pairs_for!($self.state, cnot, pairs), - (CZ, TwoQubitBatch(pairs)) => batch_pairs_for!($self.state, cz, pairs), - (RXX, TwoQubitBatchAndFloat(pairs, angle)) => { - batch_pairs_for!($self.state, rxx, pairs, *angle) - } - (RYY, TwoQubitBatchAndFloat(pairs, angle)) => { - batch_pairs_for!($self.state, ryy, pairs, *angle) - } - (RZZ, TwoQubitBatchAndFloat(pairs, angle)) => { - batch_pairs_for!($self.state, rzz, pairs, *angle) - } - (Depolarize2, TwoQubitBatchAndFloat(pairs, p)) => { - batch_pairs_for!($self.state, depolarize2, pairs, *p) - } - (TwoQubitPauliError, TwoQubitBatchAndFloatArr15(pairs, ps)) => { - batch_pairs_for!($self.state, two_qubit_pauli_error, pairs, *ps) - } - - // Not supported on either backend (Decision 11 + Gate Support - // Matrix). Loss / CorrelatedLoss handling differs by backend - // and lives in the caller's impl block, not this macro. - (Measure | Reset, _) => { - return Err(eyre!("{} is not supported on the {} backend", $inst, $backend)); - } - - // T / T_dag / U3 are listed as supported on PauliSum in the - // plan's Gate Support Matrix, but ppvm-runtime does not yet - // implement TGate or U3Gate for PauliSum (only for - // GeneralizedTableau). - (T | TAdj | U3, _) => { - return Err(eyre!( - "{} on {} requires upstream ppvm-runtime support that is not yet implemented", - $inst, - $backend - )); - } - - // Trace: parse the resolved pattern string and compute the - // trace. Per plan Decision 9, parsing happens on every - // execution; no module-load caching. - (Trace, PauliPatternStr(s)) => { - let pat = PauliPattern::parse(s) - .map_err(|e| eyre!("invalid Pauli pattern `{}`: {:?}", s, e))?; - let value = $self.state.trace(&pat); - return Ok(Effects::one(CircuitOutcomeEffect::Trace(TraceEffect { - value, - }))); - } - - // Fallback (mismatched shapes, etc.) - (inst, msg) => { - return Err(eyre!( - "Invalid circuit instruction arguments {:?} for instruction {:?} on the {} backend", - msg, - inst, - $backend - )); - } - }; - Ok(Effects::None) - }}; -} - -/// PauliSum-backed executor (Heisenberg picture). Holds a `PauliSum` and -/// answers the same `CircuitInstruction` vocabulary as `CircuitExecutor`, -/// but without measurement / reset / loss support. -pub struct PauliSumExecutor> { - pub state: PauliSum, - /// Snapshot of the seeded observable, restored by `reset`. - initial: PauliSum, -} - -#[component(instruction = CircuitInstruction, message = CircuitMessage, effect = CircuitOutcomeEffect)] -impl PauliSumExecutor -where - T: Config, - for<'a> PauliSum: Trace<'a, PauliPattern, Output = f64>, -{ - fn execute( - &mut self, - inst: CircuitInstruction, - msg: CircuitMessage, - ) -> Result> { - self.execute_instruction(&inst, &msg) - } - - fn execute_instruction( - &mut self, - inst: &CircuitInstruction, - msg: &CircuitMessage, - ) -> Result> { - use CircuitInstruction::*; - - if matches!(inst, Loss | CorrelatedLoss) { - return Err(eyre!( - "{inst} is not supported on the PauliSum backend; use the LossyPauliSum backend instead" - )); - } - - dispatch_common_paulisum!(self, inst, msg, "PauliSum") - } -} - -impl vihaco::Reset for PauliSumExecutor -where - T: Config, - PauliSum: Clone, -{ - fn reset(&mut self) { - self.state = self.initial.clone(); - } -} - -/// LossyPauliSum-backed executor. Same dispatch as `PauliSumExecutor` plus -/// `Loss` / `CorrelatedLoss` channels. The concrete `T` used by the -/// enclosing `Circuit::LossyPauliSum` variant is a `Config` whose -/// `PauliWordType` is `LossyPauliWord` (see `LossyPauliSumConfig`). -pub struct LossyPauliSumExecutor> { - pub state: PauliSum, - /// Snapshot of the seeded observable, restored by `reset`. - initial: PauliSum, -} - -#[component(instruction = CircuitInstruction, message = CircuitMessage, effect = CircuitOutcomeEffect)] -impl LossyPauliSumExecutor -where - T: Config, - for<'a> PauliSum: Trace<'a, PauliPattern, Output = f64>, -{ - fn execute( - &mut self, - inst: CircuitInstruction, - msg: CircuitMessage, - ) -> Result> { - self.execute_instruction(&inst, &msg) - } - - fn execute_instruction( - &mut self, - inst: &CircuitInstruction, - msg: &CircuitMessage, - ) -> Result> { - use CircuitInstruction::*; - use CircuitMessage::*; - - // Loss / CorrelatedLoss are the only instructions that differ from - // PauliSum; handle them here then delegate everything else to the - // shared dispatch. - match (inst, msg) { - (Loss, &QubitAndFloat(addr, p)) => { - self.state.loss_channel(addr, p); - return Ok(Effects::None); - } - (CorrelatedLoss, TwoQubitAndFloatArr3(addr0, addr1, ps)) => { - self.state.correlated_loss_channel(*addr0, *addr1, *ps); - return Ok(Effects::None); - } - (Loss, QubitBatchAndFloat(addrs, p)) => { - batch_for!(self.state, loss_channel, addrs, *p); - return Ok(Effects::None); - } - (CorrelatedLoss, TwoQubitBatchAndFloatArr3(pairs, ps)) => { - batch_pairs_for!(self.state, correlated_loss_channel, pairs, *ps); - return Ok(Effects::None); - } - _ => {} - } - - dispatch_common_paulisum!(self, inst, msg, "LossyPauliSum") - } -} - -impl vihaco::Reset for LossyPauliSumExecutor -where - T: Config, - PauliSum: Clone, -{ - fn reset(&mut self) { - self.state = self.initial.clone(); - } -} - -/// Tableau-backed inner enum (Schrödinger picture). Carries the six -/// size-bucketed `CircuitExecutor` variants; bucket is picked from `n_qubits`. -pub enum TableauCircuit { - Bits64(CircuitExecutor, usize, Vec<(Complex64, usize)>>), - Bits128(CircuitExecutor, u128, Vec<(Complex64, u128)>>), - Bits256(CircuitExecutor, U256, Vec<(Complex64, U256)>>), - Bits512(CircuitExecutor, U512, Vec<(Complex64, U512)>>), - Bits1024(CircuitExecutor, U1024, Vec<(Complex64, U1024)>>), - Bits2048(CircuitExecutor, U2048, Vec<(Complex64, U2048)>>), -} - -impl TableauCircuit { - pub fn new(n_qubits: usize, coefficient_threshold: f64) -> Self { - if n_qubits <= 64 { - let tab = GeneralizedTableau::new(n_qubits, coefficient_threshold); - Self::Bits64(CircuitExecutor { tab }) - } else if n_qubits <= 128 { - let tab = GeneralizedTableau::new(n_qubits, coefficient_threshold); - Self::Bits128(CircuitExecutor { tab }) - } else if n_qubits <= 256 { - let tab = GeneralizedTableau::new(n_qubits, coefficient_threshold); - Self::Bits256(CircuitExecutor { tab }) - } else if n_qubits <= 512 { - let tab = GeneralizedTableau::new(n_qubits, coefficient_threshold); - Self::Bits512(CircuitExecutor { tab }) - } else if n_qubits <= 1024 { - let tab = GeneralizedTableau::new(n_qubits, coefficient_threshold); - Self::Bits1024(CircuitExecutor { tab }) - } else if n_qubits <= 2048 { - let tab = GeneralizedTableau::new(n_qubits, coefficient_threshold); - Self::Bits2048(CircuitExecutor { tab }) - } else { - panic!("No matching executor for {} qubits", n_qubits); - } - } - - /// Same as [`TableauCircuit::new`], but seed the RNG deterministically so a - /// shot is reproducible. - pub fn new_with_seed(n_qubits: usize, coefficient_threshold: f64, seed: u64) -> Self { - macro_rules! seeded { - ($variant:ident) => {{ - let tab = GeneralizedTableau::new_with_seed(n_qubits, coefficient_threshold, seed); - Self::$variant(CircuitExecutor { tab }) - }}; - } - if n_qubits <= 64 { - seeded!(Bits64) - } else if n_qubits <= 128 { - seeded!(Bits128) - } else if n_qubits <= 256 { - seeded!(Bits256) - } else if n_qubits <= 512 { - seeded!(Bits512) - } else if n_qubits <= 1024 { - seeded!(Bits1024) - } else if n_qubits <= 2048 { - seeded!(Bits2048) - } else { - panic!("No matching executor for {} qubits", n_qubits); - } - } - - fn execute_instruction( - &mut self, - inst: &CircuitInstruction, - msg: &CircuitMessage, - ) -> Result> { - match self { - Self::Bits64(ex) => ex.execute_instruction(inst, msg), - Self::Bits128(ex) => ex.execute_instruction(inst, msg), - Self::Bits256(ex) => ex.execute_instruction(inst, msg), - Self::Bits512(ex) => ex.execute_instruction(inst, msg), - Self::Bits1024(ex) => ex.execute_instruction(inst, msg), - Self::Bits2048(ex) => ex.execute_instruction(inst, msg), - } - } - - pub fn state_string(&self) -> String { - match self { - Self::Bits64(ex) => ex.tab.to_string(), - Self::Bits128(ex) => ex.tab.to_string(), - Self::Bits256(ex) => ex.tab.to_string(), - Self::Bits512(ex) => ex.tab.to_string(), - Self::Bits1024(ex) => ex.tab.to_string(), - Self::Bits2048(ex) => ex.tab.to_string(), - } - } -} - -impl vihaco::Reset for TableauCircuit { - fn reset(&mut self) { - match self { - Self::Bits64(ex) => ex.reset(), - Self::Bits128(ex) => ex.reset(), - Self::Bits256(ex) => ex.reset(), - Self::Bits512(ex) => ex.reset(), - Self::Bits1024(ex) => ex.reset(), - Self::Bits2048(ex) => ex.reset(), - }; - } -} - -/// PauliSum-backed inner enum (Heisenberg picture). Per Decision 7 of the plan, -/// the size buckets carry `[u8; N]`-storage `ByteFxHashF64` configs (N = 8, 16, -/// …, 256) rather than the tableau's `[u64; N]` configs; bucket labels match -/// the semantic qubit count (`Bits64` = 64 qubits) so the outer enum's dispatch -/// is uniform across backends. -pub enum PauliSumCircuit { - Bits64(PauliSumExecutor>), - Bits128(PauliSumExecutor>), - Bits256(PauliSumExecutor>), - Bits512(PauliSumExecutor>), - Bits1024(PauliSumExecutor>), - Bits2048(PauliSumExecutor>), -} - -impl PauliSumCircuit { - /// Build a PauliSum-backed circuit, seeding the state with every term: - /// `for (word, coef) in terms { state += (word, coef); }`. Words must - /// already be validated against `info.n_qubits` by the caller. - pub fn new(info: &PPVMDeviceInfo, terms: &[(String, f64)]) -> Self { - macro_rules! build { - ($variant:ident, $N:literal) => {{ - let mut state = PauliSum::>::builder() - .n_qubits(info.n_qubits) - .strategy(paulisum_strategy(info)) - .build(); - for (word, coef) in terms { - state += (word.as_str(), *coef); - } - let initial = state.clone(); - Self::$variant(PauliSumExecutor { state, initial }) - }}; - } - if info.n_qubits <= 64 { - build!(Bits64, 8) - } else if info.n_qubits <= 128 { - build!(Bits128, 16) - } else if info.n_qubits <= 256 { - build!(Bits256, 32) - } else if info.n_qubits <= 512 { - build!(Bits512, 64) - } else if info.n_qubits <= 1024 { - build!(Bits1024, 128) - } else if info.n_qubits <= 2048 { - build!(Bits2048, 256) - } else { - panic!("No matching PauliSum executor for {} qubits", info.n_qubits); - } - } - - fn execute_instruction( - &mut self, - inst: &CircuitInstruction, - msg: &CircuitMessage, - ) -> Result> { - match self { - Self::Bits64(ex) => ex.execute_instruction(inst, msg), - Self::Bits128(ex) => ex.execute_instruction(inst, msg), - Self::Bits256(ex) => ex.execute_instruction(inst, msg), - Self::Bits512(ex) => ex.execute_instruction(inst, msg), - Self::Bits1024(ex) => ex.execute_instruction(inst, msg), - Self::Bits2048(ex) => ex.execute_instruction(inst, msg), - } - } - - pub fn state_string(&self) -> String { - match self { - Self::Bits64(ex) => ex.state.to_string(), - Self::Bits128(ex) => ex.state.to_string(), - Self::Bits256(ex) => ex.state.to_string(), - Self::Bits512(ex) => ex.state.to_string(), - Self::Bits1024(ex) => ex.state.to_string(), - Self::Bits2048(ex) => ex.state.to_string(), - } - } -} - -impl vihaco::Reset for PauliSumCircuit { - fn reset(&mut self) { - match self { - Self::Bits64(ex) => ex.reset(), - Self::Bits128(ex) => ex.reset(), - Self::Bits256(ex) => ex.reset(), - Self::Bits512(ex) => ex.reset(), - Self::Bits1024(ex) => ex.reset(), - Self::Bits2048(ex) => ex.reset(), - }; - } -} - -/// LossyPauliSum-backed inner enum. Identical shape to [`PauliSumCircuit`] -/// but with `LossyPauliWord`-keyed configs so loss-channel methods dispatch. -pub enum LossyPauliSumCircuit { - Bits64(LossyPauliSumExecutor>), - Bits128(LossyPauliSumExecutor>), - Bits256(LossyPauliSumExecutor>), - Bits512(LossyPauliSumExecutor>), - Bits1024(LossyPauliSumExecutor>), - Bits2048(LossyPauliSumExecutor>), -} - -impl LossyPauliSumCircuit { - /// Build a LossyPauliSum-backed circuit, seeding every term via - /// `state += (word, coef)`. Words must already be validated against - /// `info.n_qubits` by the caller. - pub fn new(info: &PPVMDeviceInfo, terms: &[(String, f64)]) -> Self { - macro_rules! build { - ($variant:ident, $N:literal) => {{ - let mut state = PauliSum::>::builder() - .n_qubits(info.n_qubits) - .strategy(paulisum_strategy(info)) - .build(); - for (word, coef) in terms { - state += (word.as_str(), *coef); - } - let initial = state.clone(); - Self::$variant(LossyPauliSumExecutor { state, initial }) - }}; - } - if info.n_qubits <= 64 { - build!(Bits64, 8) - } else if info.n_qubits <= 128 { - build!(Bits128, 16) - } else if info.n_qubits <= 256 { - build!(Bits256, 32) - } else if info.n_qubits <= 512 { - build!(Bits512, 64) - } else if info.n_qubits <= 1024 { - build!(Bits1024, 128) - } else if info.n_qubits <= 2048 { - build!(Bits2048, 256) - } else { - panic!( - "No matching LossyPauliSum executor for {} qubits", - info.n_qubits - ); - } - } - - fn execute_instruction( - &mut self, - inst: &CircuitInstruction, - msg: &CircuitMessage, - ) -> Result> { - match self { - Self::Bits64(ex) => ex.execute_instruction(inst, msg), - Self::Bits128(ex) => ex.execute_instruction(inst, msg), - Self::Bits256(ex) => ex.execute_instruction(inst, msg), - Self::Bits512(ex) => ex.execute_instruction(inst, msg), - Self::Bits1024(ex) => ex.execute_instruction(inst, msg), - Self::Bits2048(ex) => ex.execute_instruction(inst, msg), - } - } - - pub fn state_string(&self) -> String { - match self { - Self::Bits64(ex) => ex.state.to_string(), - Self::Bits128(ex) => ex.state.to_string(), - Self::Bits256(ex) => ex.state.to_string(), - Self::Bits512(ex) => ex.state.to_string(), - Self::Bits1024(ex) => ex.state.to_string(), - Self::Bits2048(ex) => ex.state.to_string(), - } - } -} - -impl vihaco::Reset for LossyPauliSumCircuit { - fn reset(&mut self) { - match self { - Self::Bits64(ex) => ex.reset(), - Self::Bits128(ex) => ex.reset(), - Self::Bits256(ex) => ex.reset(), - Self::Bits512(ex) => ex.reset(), - Self::Bits1024(ex) => ex.reset(), - Self::Bits2048(ex) => ex.reset(), - }; - } -} - -/// Outer `Circuit` enum: backend selector. Picks one of the three inner enums -/// based on `info.backend` at construction time; from there, every per-step -/// call routes outer → inner → executor. -pub enum Circuit { - Tableau(TableauCircuit), - PauliSum(PauliSumCircuit), - LossyPauliSum(LossyPauliSumCircuit), -} - -#[component(instruction = CircuitInstruction, message = CircuitMessage, effect = CircuitOutcomeEffect)] -impl Circuit { - /// Build a Tableau-backed circuit. Tableau init only needs `n_qubits` and - /// `coefficient_threshold` from `info`; no observable required. - pub fn tableau(info: &PPVMDeviceInfo) -> Self { - Self::Tableau(TableauCircuit::new( - info.n_qubits, - info.coefficient_threshold, - )) - } - - /// Same as [`Circuit::tableau`], but seed the tableau RNG deterministically - /// so a shot is reproducible. - pub fn tableau_with_seed(info: &PPVMDeviceInfo, seed: u64) -> Self { - Self::Tableau(TableauCircuit::new_with_seed( - info.n_qubits, - info.coefficient_threshold, - seed, - )) - } - - /// Build a PauliSum-backed circuit, seeding the state with every term in - /// `terms`. Each `(word, coef)` is added via `state += (word, coef)`; the - /// caller is responsible for having parsed/validated the words against - /// `info.n_qubits` (see `parse_observable_terms` in `composite.rs`). - pub fn paulisum(info: &PPVMDeviceInfo, terms: &[(String, f64)]) -> Self { - Self::PauliSum(PauliSumCircuit::new(info, terms)) - } - - /// Build a LossyPauliSum-backed circuit. Same contract as - /// [`Circuit::paulisum`]. - pub fn lossy_paulisum(info: &PPVMDeviceInfo, terms: &[(String, f64)]) -> Self { - Self::LossyPauliSum(LossyPauliSumCircuit::new(info, terms)) - } - - fn execute( - &mut self, - inst: CircuitInstruction, - msg: CircuitMessage, - ) -> Result> { - self.execute_instruction(&inst, &msg) - } - - fn execute_instruction( - &mut self, - inst: &CircuitInstruction, - msg: &CircuitMessage, - ) -> Result> { - match self { - Self::Tableau(c) => c.execute_instruction(inst, msg), - Self::PauliSum(c) => c.execute_instruction(inst, msg), - Self::LossyPauliSum(c) => c.execute_instruction(inst, msg), - } - } - - /// Render the current state. Used by the REPL's `show` command. - pub fn state_string(&self) -> String { - match self { - Self::Tableau(c) => c.state_string(), - Self::PauliSum(c) => c.state_string(), - Self::LossyPauliSum(c) => c.state_string(), - } - } -} - -#[observe(CircuitEffect, effect=CircuitOutcomeEffect)] -impl Circuit { - fn observe_circuit_effect( - &mut self, - effect: &CircuitEffect, - ) -> Result> { - self.execute_instruction(&effect.inst, &effect.msg) - } -} - -impl vihaco::Reset for Circuit { - fn reset(&mut self) { - match self { - Self::Tableau(c) => c.reset(), - Self::PauliSum(c) => c.reset(), - Self::LossyPauliSum(c) => c.reset(), - }; - } -} - -impl Default for Circuit { - fn default() -> Self { - // Default backend is Tableau, which doesn't require an observable. - Self::tableau(&PPVMDeviceInfo::default()) - } -} diff --git a/crates/ppvm-vihaco/src/measurements.rs b/crates/ppvm-vihaco/src/measurements.rs deleted file mode 100644 index 1881c0d51..000000000 --- a/crates/ppvm-vihaco/src/measurements.rs +++ /dev/null @@ -1,93 +0,0 @@ -// SPDX-FileCopyrightText: 2026 The PPVM Authors -// SPDX-License-Identifier: Apache-2.0 - -use eyre::Result; -use smallvec::SmallVec; -use vihaco::{Effects, observe}; - -/// Measurement results are represented as an integer enum: -/// 0: state |0> -/// 1: state |1> -/// 2: qubit has been lost prior to measurement -/// In bytecode, this is represented as a u32 integer, which is simpler than -/// e.g. two boolean values and matches semantics elsewhere -#[repr(u8)] -#[derive(Debug, Clone, Copy, Eq, PartialEq)] -pub enum MeasurementOutcome { - Zero = 0, - One = 1, - Lost = 2, -} - -pub type MeasurementResult = SmallVec<[MeasurementOutcome; 8]>; - -#[derive(Debug, Clone, Eq, PartialEq)] -pub struct MeasurementEffect { - pub measurement_results: MeasurementResult, -} - -impl From> for MeasurementOutcome { - fn from(m: Option) -> Self { - match m { - Some(false) => Self::Zero, - Some(true) => Self::One, - None => Self::Lost, - } - } -} - -#[derive(Debug, Default)] -pub struct MeasurementObserver { - pub record: Vec, -} - -#[observe(MeasurementEffect)] -impl MeasurementObserver { - fn observe_measurement_effect(&mut self, effect: &MeasurementEffect) -> Result> { - self.record.push(effect.measurement_results.clone()); - Ok(Effects::none()) - } -} - -/// Per-step trace value emitted by the `Trace` instruction on the PauliSum and -/// LossyPauliSum backends. The plan's Decision 5 keeps trace and measurement -/// records as two parallel streams; this effect feeds the trace stream. -#[derive(Debug, Clone, PartialEq)] -pub struct TraceEffect { - pub value: f64, -} - -#[derive(Debug, Default)] -pub struct TraceObserver { - pub record: Vec, -} - -#[observe(TraceEffect)] -impl TraceObserver { - fn observe_trace_effect(&mut self, effect: &TraceEffect) -> Result> { - self.record.push(effect.value); - Ok(Effects::none()) - } -} - -/// Union of the two effect types a circuit instruction can produce. The plan's -/// structural note (Task 6) calls for broadening the `#[component(..., effect = -/// ...)]` annotation on `Circuit` and the executors to a union so a single -/// `Trace` (or `Measure`) instruction can fan out to the right observer. -#[derive(Debug, Clone, PartialEq)] -pub enum CircuitOutcomeEffect { - Measurement(MeasurementEffect), - Trace(TraceEffect), -} - -impl From for CircuitOutcomeEffect { - fn from(value: MeasurementEffect) -> Self { - Self::Measurement(value) - } -} - -impl From for CircuitOutcomeEffect { - fn from(value: TraceEffect) -> Self { - Self::Trace(value) - } -} From e2a453097604c3ce03df4ba58d14a3cdd8e0fc42 Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Wed, 8 Jul 2026 11:34:19 +0200 Subject: [PATCH 4/6] test(tableau): add reset_all unit tests MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Cover `GeneralizedTableau::reset_all` and `Tableau::reset_all`: - restores a fresh |0…0⟩ tableau (rows + single identity coefficient) after non-Clifford branching, - clears the measurement record (regression guard for the review fix), - clears per-qubit loss flags, - `Tableau::reset_all` restores the fresh identity rows. Co-Authored-By: Claude Opus 4.8 (1M context) --- crates/ppvm-tableau/src/data.rs | 71 +++++++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) diff --git a/crates/ppvm-tableau/src/data.rs b/crates/ppvm-tableau/src/data.rs index 386126343..a386fc147 100644 --- a/crates/ppvm-tableau/src/data.rs +++ b/crates/ppvm-tableau/src/data.rs @@ -1472,4 +1472,75 @@ mod tests { snapshot_tableau(&tab2.tableau) ); } + + // ─── reset_all ──────────────────────────────────────────────────── + + /// `GeneralizedTableau::reset_all` restores the full state to a fresh + /// `|0…0⟩` tableau: identical stabilizer/destabilizer rows and a single + /// identity coefficient, even after non-Clifford branching. + #[test] + fn reset_all_restores_fresh_state() { + let mut tab: TestTableau = GeneralizedTableau::new(3, 1e-12); + let fresh: TestTableau = GeneralizedTableau::new(3, 1e-12); + + tab.h(0); + tab.cnot(0, 1); + tab.ry(2, 0.7); // non-Clifford: branches the coefficient vector + assert!( + tab.coefficients.iter().count() > 1, + "rotation should branch the coefficient vector" + ); + + tab.reset_all(); + + assert_eq!( + snapshot_tableau(&tab.tableau), + snapshot_tableau(&fresh.tableau) + ); + let coeffs: Vec<_> = tab.coefficients.iter().copied().collect(); + let fresh_coeffs: Vec<_> = fresh.coefficients.iter().copied().collect(); + assert_eq!(coeffs, fresh_coeffs); + } + + /// A full reset clears the measurement record. Regression guard: an earlier + /// version left it intact, so `current_measurement_record` returned stale + /// outcomes after a reset. + #[test] + fn reset_all_clears_measurement_record() { + let mut tab: TestTableau = GeneralizedTableau::new(2, 1e-12); + tab.append_measurement_record(Some(true)); + tab.append_measurement_record(None); + assert_eq!(tab.current_measurement_record().len(), 2); + + tab.reset_all(); + + assert!(tab.current_measurement_record().is_empty()); + } + + /// A full reset clears per-qubit loss flags. + #[test] + fn reset_all_clears_loss_flags() { + let mut tab: TestTableau = GeneralizedTableau::new(3, 1e-12); + tab.is_lost[0] = true; + tab.is_lost[2] = true; + + tab.reset_all(); + + assert!(tab.is_lost.iter().all(|&lost| !lost)); + } + + /// `Tableau::reset_all` restores the fresh identity tableau rows. + #[test] + fn tableau_reset_all_restores_fresh_rows() { + let mut tab: Tableau = Tableau::new(4); + let fresh: Tableau = Tableau::new(4); + + tab.h(0); + tab.s(1); + tab.h(3); + + tab.reset_all(); + + assert_eq!(snapshot_tableau(&tab), snapshot_tableau(&fresh)); + } } From feb1eef7d2f4b95121b50808ba0b74c3d4341589 Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Wed, 8 Jul 2026 11:37:44 +0200 Subject: [PATCH 5/6] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- crates/ppvm-tableau/src/expectation.rs | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/crates/ppvm-tableau/src/expectation.rs b/crates/ppvm-tableau/src/expectation.rs index 9efe04ea0..f21349551 100644 --- a/crates/ppvm-tableau/src/expectation.rs +++ b/crates/ppvm-tableau/src/expectation.rs @@ -39,9 +39,7 @@ where + ToPrimitive + Debug + std::ops::Mul - + PartialOrd - + Send - + Sync, + + PartialOrd, Complex: std::ops::Mul> + From + std::ops::MulAssign @@ -49,7 +47,7 @@ where + One + ComplexFloat + Copy, - I: TableauIndex + Debug + Send + Sync, + I: TableauIndex + Debug, { /// `⟨ψ|word|ψ⟩` for the multi-qubit Pauli `word`. /// From 74b7e7e478aff9b2d401336ffd425bc0649ea1f1 Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Wed, 8 Jul 2026 11:50:37 +0200 Subject: [PATCH 6/6] Fix cargo fmt --- crates/ppvm-tableau/src/expectation.rs | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/crates/ppvm-tableau/src/expectation.rs b/crates/ppvm-tableau/src/expectation.rs index f21349551..6ac17e11c 100644 --- a/crates/ppvm-tableau/src/expectation.rs +++ b/crates/ppvm-tableau/src/expectation.rs @@ -32,14 +32,8 @@ where T: Config, <::Storage as BitView>::Store: PrimInt, C: SparseVector, I> + Debug, - T::Coeff: One - + Zero - + Clone - + num::Num - + ToPrimitive - + Debug - + std::ops::Mul - + PartialOrd, + T::Coeff: + One + Zero + Clone + num::Num + ToPrimitive + Debug + std::ops::Mul + PartialOrd, Complex: std::ops::Mul> + From + std::ops::MulAssign