From 8f5f7f4c07ec2475bf37680e8b4d4b89aea8a99f Mon Sep 17 00:00:00 2001 From: Jonathan Wurtz Date: Tue, 6 Oct 2026 12:59:44 -0400 Subject: [PATCH 1/6] feat(tableau): add post-selected project measurement (#239) Add GeneralizedTableau::project(addr0, outcome) -> Result, a forced-outcome Z measurement that projects the state onto the requested outcome, appends it to the measurement record, and returns its probability. Built on the existing compute_overlap_case_{a,b} and project_case_{a,b} helpers. Lost qubits and zero-probability outcomes return an error and leave the state unchanged. Expose it in Python as GeneralizedTableau.project(addr0, value) -> float, raising NotImplementedError for loss and ValueError for zero probability. Co-Authored-By: Claude Opus 5.5 --- .../src/interface_tableau.rs | 12 + crates/ppvm-tableau/src/lib.rs | 1 + crates/ppvm-tableau/src/measure.rs | 134 ++++++++++ crates/ppvm-tableau/tests/project.rs | 245 ++++++++++++++++++ ppvm-python/src/ppvm/_core.pyi | 1 + ppvm-python/src/ppvm/generalized_tableau.py | 31 +++ .../test/generalized_tableau/test_project.py | 130 ++++++++++ 7 files changed, 554 insertions(+) create mode 100644 crates/ppvm-tableau/tests/project.rs create mode 100644 ppvm-python/test/generalized_tableau/test_project.py diff --git a/crates/ppvm-python-native/src/interface_tableau.rs b/crates/ppvm-python-native/src/interface_tableau.rs index 52706c802..e8ca3cef7 100644 --- a/crates/ppvm-python-native/src/interface_tableau.rs +++ b/crates/ppvm-python-native/src/interface_tableau.rs @@ -47,6 +47,18 @@ macro_rules! create_interface { measurement_to_u8(self.inner.measure(addr0)) as i64 } + /// Post-select qubit `addr0` onto `outcome` and return its probability. + pub fn project(&mut self, addr0: usize, outcome: bool) -> PyResult { + self.inner.project(addr0, outcome).map_err(|e| match e { + ProjectError::QubitLost(_) => { + pyo3::exceptions::PyNotImplementedError::new_err(e.to_string()) + } + ProjectError::ZeroProbability { .. } => { + pyo3::exceptions::PyValueError::new_err(e.to_string()) + } + }) + } + pub fn measure_many(&mut self, targets: Vec) -> Vec { self.inner .measure_many(targets.as_slice()) diff --git a/crates/ppvm-tableau/src/lib.rs b/crates/ppvm-tableau/src/lib.rs index 0133fb2fc..d0c916a1a 100644 --- a/crates/ppvm-tableau/src/lib.rs +++ b/crates/ppvm-tableau/src/lib.rs @@ -96,6 +96,7 @@ pub mod tableau_like; /// Convenience re-exports for downstream code. pub mod prelude { pub use crate::data::{GeneralizedTableau, Tableau}; + pub use crate::measure::{PROJECT_ZERO_TOL, ProjectError}; pub use crate::sparsevec::SparseVector; pub use crate::tableau_index::TableauIndex; pub use crate::tableau_like::TableauLike; diff --git a/crates/ppvm-tableau/src/measure.rs b/crates/ppvm-tableau/src/measure.rs index e386192c8..a21997458 100644 --- a/crates/ppvm-tableau/src/measure.rs +++ b/crates/ppvm-tableau/src/measure.rs @@ -83,6 +83,45 @@ impl Default for MeasureScratch { } } +/// Outcome probabilities below this are treated as zero by +/// [`project`](GeneralizedTableau::project): the projected state cannot be +/// normalized, so the projection is refused. +pub const PROJECT_ZERO_TOL: f64 = 1e-12; + +/// Error returned by [`project`](GeneralizedTableau::project). +#[derive(Debug, Clone, PartialEq)] +pub enum ProjectError { + /// The target qubit is lost; projecting a lost qubit is not implemented. + QubitLost(usize), + /// The requested outcome has (numerically) zero probability. + ZeroProbability { + addr0: usize, + outcome: bool, + prob: f64, + }, +} + +impl std::fmt::Display for ProjectError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + ProjectError::QubitLost(addr0) => { + write!(f, "cannot project qubit {addr0}: qubit is lost") + } + ProjectError::ZeroProbability { + addr0, + outcome, + prob, + } => write!( + f, + "cannot project qubit {addr0} onto outcome {}: probability {prob:e} is zero", + *outcome as u8 + ), + } + } +} + +impl std::error::Error for ProjectError {} + impl, I>> LossyMeasure for GeneralizedTableau where @@ -374,6 +413,101 @@ where } } + /// Post-select qubit `addr0` onto the Z-basis `outcome` (`false` = |0⟩, + /// `true` = |1⟩) and return the probability of that outcome. + /// + /// This is a non-physical, forced-outcome version of + /// [`measure`](LossyMeasure::measure): the state is projected onto + /// `outcome` and renormalized, and the outcome is appended to the + /// measurement record, but no RNG draw is made. Chaining it over every + /// qubit gives `P(z) = |⟨z|ψ⟩|²` as the product of the returned + /// probabilities. + /// + /// Acts on the current pure state of this trajectory: any noise or loss + /// channels applied earlier have already been sampled. + /// + /// # Errors + /// + /// - [`ProjectError::QubitLost`] if `addr0` is lost. + /// - [`ProjectError::ZeroProbability`] if the outcome's probability is + /// below [`PROJECT_ZERO_TOL`]. + /// + /// The state is left unchanged when an error is returned. + pub fn project(&mut self, addr0: usize, outcome: bool) -> Result { + if self.is_lost[addr0] { + return Err(ProjectError::QubitLost(addr0)); + } + + let (phase_decomp, stab_anticomm_bits, destab_anticomm_bits) = + self.compute_decomposition(addr0, Pauli::Z); + + if stab_anticomm_bits == I::zero() { + // Case b: Z is a stabilizer — self-pairing overlap, then filter. + let entries: Vec<(Complex, I)> = self.coefficients.iter().copied().collect(); + let z = Self::compute_overlap_case_b(&entries, phase_decomp, destab_anticomm_bits); + let prob = Self::outcome_probability(z, outcome); + if prob < PROJECT_ZERO_TOL { + return Err(ProjectError::ZeroProbability { + addr0, + outcome, + prob, + }); + } + // `project_case_b` refills `self.coefficients` from `entries`. + self.coefficients = C::new(); + self.project_case_b(&entries, outcome, phase_decomp, destab_anticomm_bits); + self.measurement_record.push(Some(outcome)); + Ok(prob) + } else { + // Case a: Z is not a stabilizer — cross-index pairing via HashMap. + let odd_phase_mask = self.odd_phase_destabilizer_mask(); + let coeff_map: HashMap> = + self.coefficients.iter().map(|&(c, i)| (i, c)).collect(); + let z = Self::compute_overlap_case_a( + &coeff_map, + phase_decomp, + destab_anticomm_bits, + stab_anticomm_bits, + odd_phase_mask, + ); + let prob = Self::outcome_probability(z, outcome); + if prob < PROJECT_ZERO_TOL { + return Err(ProjectError::ZeroProbability { + addr0, + outcome, + prob, + }); + } + + // `project_case_a` expects the coefficients drained into + // `scratch.coeff_map` and refills `self.coefficients` from it. + self.coefficients = C::new(); + let mut scratch = MeasureScratch::new(); + scratch.coeff_map = coeff_map; + scratch.odd_phase_mask = Some(odd_phase_mask); + self.project_case_a( + outcome, + &mut scratch, + phase_decomp, + stab_anticomm_bits, + destab_anticomm_bits, + addr0, + ); + self.measurement_record.push(Some(outcome)); + Ok(prob) + } + } + + /// Probability of `outcome` given `⟨Z⟩ = z`, clamped to `[0, 1]`. + fn outcome_probability(z: f64, outcome: bool) -> f64 { + let prob = if outcome { + 0.5 - 0.5 * z + } else { + 0.5 + 0.5 * z + }; + prob.clamp(0.0, 1.0) + } + pub fn project_case_a( &mut self, outcome: bool, diff --git a/crates/ppvm-tableau/tests/project.rs b/crates/ppvm-tableau/tests/project.rs new file mode 100644 index 000000000..a01746f26 --- /dev/null +++ b/crates/ppvm-tableau/tests/project.rs @@ -0,0 +1,245 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +use ppvm_pauli_sum::config::fxhash::ByteF64; +use ppvm_tableau::prelude::*; + +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() + ); +} + +/// All `4^n` Pauli strings on `n` qubits. +fn all_paulis(n: usize) -> Vec { + let mut out = vec![String::new()]; + for _ in 0..n { + out = out + .into_iter() + .flat_map(|s| ["I", "X", "Y", "Z"].map(|p| format!("{s}{p}"))) + .collect(); + } + out +} + +/// `P(b) = 2^-n Σ_T (-1)^{b·T} ⟨Z_T⟩`, computed without collapsing the state. +fn bitstring_prob_from_expectations(tab: &TestTableau, bits: &[bool]) -> f64 { + let n = bits.len(); + let mut sum = 0.0; + for mask in 0..(1usize << n) { + let w: String = (0..n) + .map(|q| if mask >> q & 1 == 1 { 'Z' } else { 'I' }) + .collect(); + let sign = (0..n).filter(|&q| mask >> q & 1 == 1 && bits[q]).count() % 2; + let e = tab.expectation(&word(&w)); + sum += if sign == 0 { e } else { -e }; + } + sum / (1usize << n) as f64 +} + +/// Chain-rule bitstring probability via successive projections. +fn bitstring_prob_by_projection(tab: &TestTableau, bits: &[bool]) -> f64 { + let mut t = tab.fork(Some(0)); + let mut prob = 1.0; + for (q, &b) in bits.iter().enumerate() { + match t.project(q, b) { + Ok(p) => prob *= p, + Err(ProjectError::ZeroProbability { .. }) => return 0.0, + Err(e) => panic!("unexpected error: {e}"), + } + } + prob +} + +fn all_bitstrings(n: usize) -> Vec> { + (0..(1usize << n)) + .map(|m| (0..n).map(|q| m >> q & 1 == 1).collect()) + .collect() +} + +/// A small entangled, non-Clifford circuit exercising both measurement cases. +fn magic_circuit() -> TestTableau { + let mut tab = TestTableau::new(3, 1e-12); + tab.h(0); + tab.t(0); + tab.cnot(0, 1); + tab.h(2); + tab.t(2); + tab.cnot(1, 2); + tab.ry(1, 0.7); + tab.t(1); + tab.h(1); + tab.rx(2, 0.3); + tab +} + +#[test] +fn project_zero_state() { + let mut tab = TestTableau::new(1, 1e-12); + assert_close(tab.project(0, false).unwrap(), 1.0, 1e-12); + assert_eq!(tab.current_measurement_record(), &[Some(false)]); +} + +#[test] +fn project_zero_probability_leaves_state_unchanged() { + let mut tab = TestTableau::new(1, 1e-12); + let err = tab.project(0, true).unwrap_err(); + assert!(matches!( + err, + ProjectError::ZeroProbability { + addr0: 0, + outcome: true, + .. + } + )); + assert_close(tab.expectation(&word("Z")), 1.0, 1e-12); + assert!(tab.current_measurement_record().is_empty()); + + // Same for a case-a (random-outcome) qubit whose outcome is still impossible: + // |+⟩ with a T phase has P(1) = 0.5, but after projecting onto 0, P(1) = 0. + let mut tab = TestTableau::new(1, 1e-12); + tab.h(0); + tab.t(0); + tab.project(0, false).unwrap(); + assert!(tab.project(0, true).is_err()); + assert_close(tab.expectation(&word("Z")), 1.0, 1e-12); +} + +#[test] +fn project_plus_state() { + for outcome in [false, true] { + let mut tab = TestTableau::new(1, 1e-12); + tab.h(0); + assert_close(tab.project(0, outcome).unwrap(), 0.5, 1e-12); + let z = if outcome { -1.0 } else { 1.0 }; + assert_close(tab.expectation(&word("Z")), z, 1e-12); + // The outcome is now deterministic. + assert_close(tab.project(0, outcome).unwrap(), 1.0, 1e-12); + assert_eq!( + tab.current_measurement_record(), + &[Some(outcome), Some(outcome)] + ); + } +} + +#[test] +fn project_ry_matches_cos_squared() { + for theta in [0.1, 0.7, 1.3, 2.0, 2.9] { + for outcome in [false, true] { + let mut tab = TestTableau::new(1, 1e-12); + tab.ry(0, theta); + let p0 = (theta / 2.0_f64).cos().powi(2); + let expected = if outcome { 1.0 - p0 } else { p0 }; + assert_close(tab.project(0, outcome).unwrap(), expected, 1e-10); + let z = if outcome { -1.0 } else { 1.0 }; + assert_close(tab.expectation(&word("Z")), z, 1e-10); + } + } +} + +#[test] +fn project_matches_measure() { + let base = magic_circuit(); + let paulis = all_paulis(3); + for q in 0..3 { + for outcome in [false, true] { + let mut projected = base.fork(Some(0)); + if projected.project(q, outcome).is_err() { + continue; + } + let mut measured = (0..1000u64) + .map(|s| { + let mut t = base.fork(Some(s)); + let m = t.measure(q); + (t, m) + }) + .find(|(_, m)| *m == Some(outcome)) + .expect("outcome with nonzero probability should be sampled") + .0; + for p in &paulis { + assert_close( + projected.expectation(&word(p)), + measured.expectation(&word(p)), + 1e-9, + ); + } + // Subsequent measurements on the two states agree in distribution; + // check that the next qubit's probabilities agree exactly. + let next = (q + 1) % 3; + let p_proj = projected.fork(Some(1)).project(next, false); + let p_meas = measured.project(next, false); + match (p_proj, p_meas) { + (Ok(a), Ok(b)) => assert_close(a, b, 1e-9), + (Err(_), Err(_)) => {} + other => panic!("mismatch: {other:?}"), + } + } + } +} + +#[test] +fn project_chain_rule_matches_expectations() { + let tab = magic_circuit(); + let mut total = 0.0; + for bits in all_bitstrings(3) { + let by_projection = bitstring_prob_by_projection(&tab, &bits); + let by_expectation = bitstring_prob_from_expectations(&tab, &bits); + assert_close(by_projection, by_expectation, 1e-9); + total += by_projection; + } + assert_close(total, 1.0, 1e-9); +} + +#[test] +fn project_chain_rule_order_independent() { + let tab = magic_circuit(); + for bits in all_bitstrings(3) { + let mut t = tab.fork(Some(0)); + let mut prob = 1.0; + for q in [2, 0, 1] { + match t.project(q, bits[q]) { + Ok(p) => prob *= p, + Err(_) => { + prob = 0.0; + break; + } + } + } + assert_close(prob, bitstring_prob_from_expectations(&tab, &bits), 1e-9); + } +} + +#[test] +fn project_ghz() { + let mut tab = TestTableau::new(3, 1e-12); + tab.h(0); + tab.cnot(0, 1); + tab.cnot(1, 2); + for bits in all_bitstrings(3) { + let expected = if bits.iter().all(|&b| b) || bits.iter().all(|&b| !b) { + 0.5 + } else { + 0.0 + }; + assert_close(bitstring_prob_by_projection(&tab, &bits), expected, 1e-12); + } +} + +#[test] +fn project_lost_qubit_errors() { + let mut tab = TestTableau::new(2, 1e-12); + tab.h(0); + tab.loss_channel(0, 1.0); + assert_eq!(tab.project(0, false), Err(ProjectError::QubitLost(0))); + assert!(tab.current_measurement_record().is_empty()); + // Other qubits are unaffected. + assert_close(tab.project(1, false).unwrap(), 1.0, 1e-12); +} diff --git a/ppvm-python/src/ppvm/_core.pyi b/ppvm-python/src/ppvm/_core.pyi index c4017324c..0926eb468 100644 --- a/ppvm-python/src/ppvm/_core.pyi +++ b/ppvm-python/src/ppvm/_core.pyi @@ -113,6 +113,7 @@ class _GeneralizedTableauBase: def __copy__(self) -> _GeneralizedTableauBase: ... def __deepcopy__(self, memo: object) -> _GeneralizedTableauBase: ... def measure(self, addr0: int) -> int: ... + def project(self, addr0: int, outcome: bool) -> float: ... def measure_many(self, targets: Sequence[int]) -> list[int]: ... def current_measurement_record(self) -> list[int]: ... def coefficients(self) -> dict[int, complex]: ... diff --git a/ppvm-python/src/ppvm/generalized_tableau.py b/ppvm-python/src/ppvm/generalized_tableau.py index f127c63f6..0c7399bf4 100644 --- a/ppvm-python/src/ppvm/generalized_tableau.py +++ b/ppvm-python/src/ppvm/generalized_tableau.py @@ -191,6 +191,37 @@ def measure(self, addr0: int) -> MeasurementResult: """ return _BY_VALUE[self._interface.measure(addr0)] + def project(self, addr0: int, value: MeasurementResult) -> float: + """Post-select the specified qubit onto a Z-basis outcome. + + A non-physical, forced-outcome version of `measure`: the state is + projected onto ``value`` and renormalized, and ``value`` is appended + to the measurement record, but no random outcome is drawn. Chaining + it over every qubit gives the bitstring probability + ``P(z) = |⟨z|ψ⟩|²`` as the product of the returned probabilities. + + This acts on the current pure state of this trajectory: any noise or + loss channels applied earlier have already been sampled. + + Args: + addr0: The index of the target qubit. + value: The outcome to project onto, ``ZERO`` or ``ONE``. + + Returns: + The probability of measuring ``value`` on the state before projection. + + Raises: + NotImplementedError: If ``value`` is ``LOST`` or the qubit has + been lost; projection with loss is not implemented. + ValueError: If ``value`` has zero probability (the state cannot + be normalized). The state is left unchanged. + """ + if value == MeasurementResult.LOST: + raise NotImplementedError("Projecting onto a LOST outcome is not implemented.") + if value not in (MeasurementResult.ZERO, MeasurementResult.ONE): + raise ValueError(f"value must be MeasurementResult.ZERO or ONE (got {value!r}).") + return self._interface.project(addr0, value == MeasurementResult.ONE) + def measure_many(self, *targets: int | Iterable[int]) -> list[MeasurementResult]: """Measure several qubits in the Z basis. diff --git a/ppvm-python/test/generalized_tableau/test_project.py b/ppvm-python/test/generalized_tableau/test_project.py new file mode 100644 index 000000000..9fed27bef --- /dev/null +++ b/ppvm-python/test/generalized_tableau/test_project.py @@ -0,0 +1,130 @@ +import itertools +import math + +import pytest + +from ppvm import GeneralizedTableau +from ppvm.generalized_tableau import MeasurementResult + +ZERO, ONE, LOST = MeasurementResult.ZERO, MeasurementResult.ONE, MeasurementResult.LOST + + +def _magic_circuit() -> GeneralizedTableau: + tab = GeneralizedTableau(3, min_abs_coeff=1e-12) + tab.h(0) + tab.t(0) + tab.cnot(0, 1) + tab.h(2) + tab.t(2) + tab.cnot(1, 2) + tab.ry(1, 0.7) + tab.t(1) + tab.h(1) + tab.rx(2, 0.3) + return tab + + +def _prob_by_projection(tab: GeneralizedTableau, bits: tuple[int, ...]) -> float: + t = tab.fork(seed=0) + prob = 1.0 + for q, b in enumerate(bits): + try: + prob *= t.project(q, MeasurementResult(b)) + except ValueError: + return 0.0 + return prob + + +def _prob_by_expectation(tab: GeneralizedTableau, bits: tuple[int, ...]) -> float: + # P(b) = 2^-n Σ_T (-1)^{b·T} ⟨Z_T⟩ + n = len(bits) + total = 0.0 + for mask in itertools.product((0, 1), repeat=n): + word = "".join("Z" if m else "I" for m in mask) + sign = (-1) ** sum(m * b for m, b in zip(mask, bits)) + total += sign * tab.expectation(word) + return total / 2**n + + +def test_project_zero_state(): + tab = GeneralizedTableau(1) + assert tab.project(0, ZERO) == pytest.approx(1.0, abs=1e-12) + assert tab.current_measurement_record() == [ZERO] + + +def test_project_zero_probability_raises_and_preserves_state(): + tab = GeneralizedTableau(1) + with pytest.raises(ValueError, match="probability"): + tab.project(0, ONE) + assert tab.expectation("Z") == pytest.approx(1.0, abs=1e-12) + assert tab.current_measurement_record() == [] + + +@pytest.mark.parametrize("value", [ZERO, ONE]) +def test_project_plus_state(value): + tab = GeneralizedTableau(1) + tab.h(0) + assert tab.project(0, value) == pytest.approx(0.5, abs=1e-12) + expected_z = 1.0 if value == ZERO else -1.0 + assert tab.expectation("Z") == pytest.approx(expected_z, abs=1e-12) + assert tab.project(0, value) == pytest.approx(1.0, abs=1e-12) + assert tab.current_measurement_record() == [value, value] + + +@pytest.mark.parametrize("theta", [0.1, 0.7, 1.3, 2.0, 2.9]) +def test_project_ry_matches_cos_squared(theta): + p0 = math.cos(theta / 2) ** 2 + for value, expected in [(ZERO, p0), (ONE, 1 - p0)]: + tab = GeneralizedTableau(1) + tab.ry(0, theta) + assert tab.project(0, value) == pytest.approx(expected, abs=1e-10) + + +def test_project_matches_measure(): + base = _magic_circuit() + paulis = ["".join(p) for p in itertools.product("IXYZ", repeat=3)] + for q in range(3): + for value in (ZERO, ONE): + projected = base.fork(seed=0) + try: + projected.project(q, value) + except ValueError: + continue + measured = next(t for s in range(1000) if (t := base.fork(seed=s)).measure(q) == value) + for p in paulis: + assert projected.expectation(p) == pytest.approx(measured.expectation(p), abs=1e-9) + + +def test_project_chain_rule_gives_bitstring_probabilities(): + tab = _magic_circuit() + total = 0.0 + for bits in itertools.product((0, 1), repeat=3): + p = _prob_by_projection(tab, bits) + assert p == pytest.approx(_prob_by_expectation(tab, bits), abs=1e-9) + total += p + assert total == pytest.approx(1.0, abs=1e-9) + + +def test_project_ghz(): + tab = GeneralizedTableau(3) + tab.h(0) + tab.cnot(0, 1) + tab.cnot(1, 2) + for bits in itertools.product((0, 1), repeat=3): + expected = 0.5 if len(set(bits)) == 1 else 0.0 + assert _prob_by_projection(tab, bits) == pytest.approx(expected, abs=1e-12) + + +def test_project_lost_value_not_implemented(): + tab = GeneralizedTableau(1) + with pytest.raises(NotImplementedError): + tab.project(0, LOST) + + +def test_project_lost_qubit_not_implemented(): + tab = GeneralizedTableau(2) + tab.loss_channel(0, 1.0) + with pytest.raises(NotImplementedError, match="lost"): + tab.project(0, ZERO) + assert tab.current_measurement_record() == [] + assert tab.project(1, ZERO) == pytest.approx(1.0, abs=1e-12) From 028beacf8b75da25fb1cf18957a8b4ffdb4130fd Mon Sep 17 00:00:00 2001 From: Jonathan Wurtz Date: Tue, 6 Oct 2026 13:21:59 -0400 Subject: [PATCH 2/6] test(tableau): cover zero-probability projection in both cases Split the zero-probability project test into case-b and case-a variants and assert which measurement path each projection takes via compute_decomposition. The case-a test uses H then RY(-pi/2), which returns to |0> while the stabilizer frame still holds X, so P(1) = 0 on the case-a path. Both tests check that the tableau, coefficients, all Pauli expectations, and the measurement record are unchanged after the error. Mirror both scenarios in the Python tests. Co-Authored-By: Claude Opus 5.5 --- crates/ppvm-tableau/tests/project.rs | 84 +++++++++++++++---- .../test/generalized_tableau/test_project.py | 49 +++++++++-- 2 files changed, 113 insertions(+), 20 deletions(-) diff --git a/crates/ppvm-tableau/tests/project.rs b/crates/ppvm-tableau/tests/project.rs index a01746f26..73f2cb5aa 100644 --- a/crates/ppvm-tableau/tests/project.rs +++ b/crates/ppvm-tableau/tests/project.rs @@ -1,6 +1,8 @@ // SPDX-FileCopyrightText: 2026 The PPVM Authors // SPDX-License-Identifier: Apache-2.0 +use std::f64::consts::FRAC_PI_2; + use ppvm_pauli_sum::config::fxhash::ByteF64; use ppvm_tableau::prelude::*; @@ -88,29 +90,81 @@ fn project_zero_state() { assert_eq!(tab.current_measurement_record(), &[Some(false)]); } +/// True if projecting `addr0` takes the case-a path (Z anticommutes with some +/// stabilizer), false for case b (Z is a stabilizer up to sign). +fn is_case_a(tab: &TestTableau, addr0: usize) -> bool { + let (_, stab_anticomm_bits, _) = tab.compute_decomposition(addr0, Pauli::Z); + stab_anticomm_bits != 0 +} + +/// Projects `addr0` onto an outcome that must have zero probability, and checks +/// that the error is returned and the tableau, coefficients, Pauli expectations, +/// and measurement record are all unchanged. +fn assert_zero_probability_preserves_state(tab: &mut TestTableau, addr0: usize, outcome: bool) { + let n = tab.n_qubits(); + let paulis = all_paulis(n); + let before_display = format!("{tab}"); + let before_expectations: Vec = paulis.iter().map(|p| tab.expectation(&word(p))).collect(); + let before_record = tab.current_measurement_record().to_vec(); + + let err = tab.project(addr0, outcome).unwrap_err(); + assert!( + matches!(err, ProjectError::ZeroProbability { addr0: a, outcome: o, .. } if a == addr0 && o == outcome), + "unexpected error: {err:?}" + ); + + assert_eq!(format!("{tab}"), before_display); + for (p, before) in paulis.iter().zip(before_expectations) { + assert_close(tab.expectation(&word(p)), before, 1e-12); + } + assert_eq!(tab.current_measurement_record(), before_record.as_slice()); +} + #[test] -fn project_zero_probability_leaves_state_unchanged() { +fn project_zero_probability_case_b() { + // |0⟩: Z is a stabilizer, P(1) = 0. let mut tab = TestTableau::new(1, 1e-12); - let err = tab.project(0, true).unwrap_err(); - assert!(matches!( - err, - ProjectError::ZeroProbability { - addr0: 0, - outcome: true, - .. - } - )); - assert_close(tab.expectation(&word("Z")), 1.0, 1e-12); - assert!(tab.current_measurement_record().is_empty()); + assert!(!is_case_a(&tab, 0)); + assert_zero_probability_preserves_state(&mut tab, 0, true); - // Same for a case-a (random-outcome) qubit whose outcome is still impossible: - // |+⟩ with a T phase has P(1) = 0.5, but after projecting onto 0, P(1) = 0. + // |+⟩ with a T phase, projected onto 0: the projection makes Z a + // stabilizer, so the impossible second projection also goes through case b. let mut tab = TestTableau::new(1, 1e-12); tab.h(0); tab.t(0); + assert!(is_case_a(&tab, 0)); tab.project(0, false).unwrap(); - assert!(tab.project(0, true).is_err()); + assert!(!is_case_a(&tab, 0)); + assert_zero_probability_preserves_state(&mut tab, 0, true); +} + +#[test] +fn project_zero_probability_case_a() { + // H then RY(-π/2) returns to |0⟩, but the rotation only reweights the + // coefficients: the stabilizer frame still holds X, so Z is not a stabilizer + // and the projection takes case a even though P(1) = 0. + let mut tab = TestTableau::new(1, 1e-12); + tab.h(0); + tab.ry(0, -FRAC_PI_2); + assert!(is_case_a(&tab, 0)); assert_close(tab.expectation(&word("Z")), 1.0, 1e-12); + assert_zero_probability_preserves_state(&mut tab, 0, true); + // The possible outcome still projects, with probability 1. + assert_close(tab.project(0, false).unwrap(), 1.0, 1e-12); + + // Same on 3 qubits, with qubit 0 in |0⟩ (case a) next to an entangled, + // non-Clifford pair on qubits 1 and 2. + let mut tab = TestTableau::new(3, 1e-12); + tab.h(0); + tab.h(1); + tab.h(2); + tab.t(1); + tab.cz(1, 2); + tab.ry(0, -FRAC_PI_2); + assert!(is_case_a(&tab, 0)); + assert!(tab.coefficients.len() > 1); + assert_zero_probability_preserves_state(&mut tab, 0, true); + assert_close(tab.project(0, false).unwrap(), 1.0, 1e-12); } #[test] diff --git a/ppvm-python/test/generalized_tableau/test_project.py b/ppvm-python/test/generalized_tableau/test_project.py index 9fed27bef..d838716bd 100644 --- a/ppvm-python/test/generalized_tableau/test_project.py +++ b/ppvm-python/test/generalized_tableau/test_project.py @@ -52,12 +52,51 @@ def test_project_zero_state(): assert tab.current_measurement_record() == [ZERO] -def test_project_zero_probability_raises_and_preserves_state(): - tab = GeneralizedTableau(1) +def _assert_zero_probability_preserves_state(tab: GeneralizedTableau, addr0: int, value): + paulis = ["".join(p) for p in itertools.product("IXYZ", repeat=tab.n_qubits)] + before_str = str(tab) + before = [tab.expectation(p) for p in paulis] + before_coeffs = tab.coefficients() + before_record = tab.current_measurement_record() with pytest.raises(ValueError, match="probability"): - tab.project(0, ONE) - assert tab.expectation("Z") == pytest.approx(1.0, abs=1e-12) - assert tab.current_measurement_record() == [] + tab.project(addr0, value) + assert str(tab) == before_str + assert tab.coefficients() == before_coeffs + assert [tab.expectation(p) for p in paulis] == pytest.approx(before, abs=1e-12) + assert tab.current_measurement_record() == before_record + + +def test_project_zero_probability_z_stabilizer(): + # |0⟩: Z is a stabilizer (case b). + tab = GeneralizedTableau(1) + _assert_zero_probability_preserves_state(tab, 0, ONE) + + # Projecting |+⟩·T onto 0 makes Z a stabilizer; projecting onto 1 is then impossible. + tab = GeneralizedTableau(1) + tab.h(0) + tab.t(0) + tab.project(0, ZERO) + _assert_zero_probability_preserves_state(tab, 0, ONE) + + +def test_project_zero_probability_z_not_stabilizer(): + # H then RY(-π/2) is |0⟩, but the stabilizer frame still holds X, so the + # projection takes the case-a path with P(1) = 0. + tab = GeneralizedTableau(1) + tab.h(0) + tab.ry(0, -math.pi / 2) + _assert_zero_probability_preserves_state(tab, 0, ONE) + assert tab.project(0, ZERO) == pytest.approx(1.0, abs=1e-12) + + # 3 qubits: qubit 0 as above, next to an entangled non-Clifford pair. + tab = GeneralizedTableau(3, min_abs_coeff=1e-12) + tab.h(0, 1, 2) + tab.t(1) + tab.cz(1, 2) + tab.ry(0, -math.pi / 2) + assert tab.num_coefficients() > 1 + _assert_zero_probability_preserves_state(tab, 0, ONE) + assert tab.project(0, ZERO) == pytest.approx(1.0, abs=1e-12) @pytest.mark.parametrize("value", [ZERO, ONE]) From 87ecd89a7f71ceb166f9724b19473c96a806ce74 Mon Sep 17 00:00:00 2001 From: Jonathan Wurtz Date: Tue, 6 Oct 2026 13:34:02 -0400 Subject: [PATCH 3/6] feat(tableau): add project_many and probability Add GeneralizedTableau::project_many(targets, outcomes), which post-selects each target in order and returns the joint probability. It is atomic: on any error the state and measurement record are restored. Mismatched lengths return the new ProjectError::LengthMismatch. Add GeneralizedTableau::probability(targets, outcomes), which forks the state and calls project_many on the fork, leaving the state unchanged. Impossible outcomes return 0.0 instead of an error. Expose both in Python with the signature (targets, values) and share the argument validation through _projection_args. Co-Authored-By: Claude Opus 5.5 --- .../src/interface_tableau.rs | 40 ++++++++-- crates/ppvm-tableau/src/measure.rs | 67 +++++++++++++++- crates/ppvm-tableau/tests/project.rs | 79 +++++++++++++++++++ ppvm-python/src/ppvm/_core.pyi | 2 + ppvm-python/src/ppvm/generalized_tableau.py | 78 ++++++++++++++++-- .../test/generalized_tableau/test_project.py | 56 +++++++++++++ 6 files changed, 308 insertions(+), 14 deletions(-) diff --git a/crates/ppvm-python-native/src/interface_tableau.rs b/crates/ppvm-python-native/src/interface_tableau.rs index e8ca3cef7..38dede5ba 100644 --- a/crates/ppvm-python-native/src/interface_tableau.rs +++ b/crates/ppvm-python-native/src/interface_tableau.rs @@ -16,6 +16,17 @@ pub(crate) fn measurement_to_u8(m: Option) -> u8 { } } +fn project_error_to_py(e: ProjectError) -> PyErr { + match e { + ProjectError::QubitLost(_) => { + pyo3::exceptions::PyNotImplementedError::new_err(e.to_string()) + } + ProjectError::ZeroProbability { .. } | ProjectError::LengthMismatch { .. } => { + pyo3::exceptions::PyValueError::new_err(e.to_string()) + } + } +} + macro_rules! create_interface { ($name: ident, $type: ident, $indexType: ident) => { #[pyclass] @@ -49,14 +60,27 @@ macro_rules! create_interface { /// Post-select qubit `addr0` onto `outcome` and return its probability. pub fn project(&mut self, addr0: usize, outcome: bool) -> PyResult { - self.inner.project(addr0, outcome).map_err(|e| match e { - ProjectError::QubitLost(_) => { - pyo3::exceptions::PyNotImplementedError::new_err(e.to_string()) - } - ProjectError::ZeroProbability { .. } => { - pyo3::exceptions::PyValueError::new_err(e.to_string()) - } - }) + self.inner + .project(addr0, outcome) + .map_err(project_error_to_py) + } + + /// Joint probability of `outcomes` on `targets`, without modifying the state. + pub fn probability(&self, targets: Vec, outcomes: Vec) -> PyResult { + self.inner + .probability(&targets, &outcomes) + .map_err(project_error_to_py) + } + + /// Post-select each target onto its outcome and return the joint probability. + pub fn project_many( + &mut self, + targets: Vec, + outcomes: Vec, + ) -> PyResult { + self.inner + .project_many(&targets, &outcomes) + .map_err(project_error_to_py) } pub fn measure_many(&mut self, targets: Vec) -> Vec { diff --git a/crates/ppvm-tableau/src/measure.rs b/crates/ppvm-tableau/src/measure.rs index a21997458..ed8601161 100644 --- a/crates/ppvm-tableau/src/measure.rs +++ b/crates/ppvm-tableau/src/measure.rs @@ -88,7 +88,8 @@ impl Default for MeasureScratch { /// normalized, so the projection is refused. pub const PROJECT_ZERO_TOL: f64 = 1e-12; -/// Error returned by [`project`](GeneralizedTableau::project). +/// Error returned by [`project`](GeneralizedTableau::project) and +/// [`project_many`](GeneralizedTableau::project_many). #[derive(Debug, Clone, PartialEq)] pub enum ProjectError { /// The target qubit is lost; projecting a lost qubit is not implemented. @@ -99,6 +100,8 @@ pub enum ProjectError { outcome: bool, prob: f64, }, + /// `project_many` was given different numbers of targets and outcomes. + LengthMismatch { targets: usize, outcomes: usize }, } impl std::fmt::Display for ProjectError { @@ -116,6 +119,10 @@ impl std::fmt::Display for ProjectError { "cannot project qubit {addr0} onto outcome {}: probability {prob:e} is zero", *outcome as u8 ), + ProjectError::LengthMismatch { targets, outcomes } => write!( + f, + "got {targets} targets but {outcomes} outcomes; they must have the same length" + ), } } } @@ -498,6 +505,64 @@ where } } + /// Post-select each `targets[k]` onto `outcomes[k]`, in order, and return + /// the joint probability of those outcomes. + /// + /// Equivalent to calling [`project`](Self::project) on each pair and + /// multiplying the returned conditional probabilities, but atomic: if any + /// projection fails, the state (including the measurement record) is + /// restored to what it was before the call. + /// + /// # Errors + /// + /// - [`ProjectError::LengthMismatch`] if `targets` and `outcomes` differ in + /// length. + /// - Any error from [`project`](Self::project) on one of the targets. + pub fn project_many( + &mut self, + targets: &[usize], + outcomes: &[bool], + ) -> Result { + if targets.len() != outcomes.len() { + return Err(ProjectError::LengthMismatch { + targets: targets.len(), + outcomes: outcomes.len(), + }); + } + let backup = self.clone(); + let mut prob = 1.0; + for (&addr0, &outcome) in targets.iter().zip(outcomes) { + match self.project(addr0, outcome) { + Ok(p) => prob *= p, + Err(e) => { + *self = backup; + return Err(e); + } + } + } + Ok(prob) + } + + /// Joint probability of measuring `outcomes` on `targets`, without + /// modifying this state. + /// + /// Forks the state and calls [`project_many`](Self::project_many) on the + /// fork. Outcomes with zero joint probability return `Ok(0.0)` rather than + /// an error. + /// + /// # Errors + /// + /// - [`ProjectError::LengthMismatch`] if `targets` and `outcomes` differ in + /// length. + /// - [`ProjectError::QubitLost`] if any target qubit is lost. + pub fn probability(&self, targets: &[usize], outcomes: &[bool]) -> Result { + // Projection never draws from the RNG, so the fork's seed is irrelevant. + match self.fork(Some(0)).project_many(targets, outcomes) { + Err(ProjectError::ZeroProbability { .. }) => Ok(0.0), + result => result, + } + } + /// Probability of `outcome` given `⟨Z⟩ = z`, clamped to `[0, 1]`. fn outcome_probability(z: f64, outcome: bool) -> f64 { let prob = if outcome { diff --git a/crates/ppvm-tableau/tests/project.rs b/crates/ppvm-tableau/tests/project.rs index 73f2cb5aa..5d308ce3a 100644 --- a/crates/ppvm-tableau/tests/project.rs +++ b/crates/ppvm-tableau/tests/project.rs @@ -297,3 +297,82 @@ fn project_lost_qubit_errors() { // Other qubits are unaffected. assert_close(tab.project(1, false).unwrap(), 1.0, 1e-12); } + +#[test] +fn project_many_gives_joint_probability() { + let tab = magic_circuit(); + for bits in all_bitstrings(3) { + let mut t = tab.fork(Some(0)); + let expected = bitstring_prob_from_expectations(&tab, &bits); + match t.project_many(&[0, 1, 2], &bits) { + Ok(p) => { + assert_close(p, expected, 1e-9); + let record: Vec> = bits.iter().map(|&b| Some(b)).collect(); + assert_eq!(t.current_measurement_record(), record.as_slice()); + } + Err(ProjectError::ZeroProbability { .. }) => assert_close(expected, 0.0, 1e-9), + Err(e) => panic!("unexpected error: {e}"), + } + } +} + +#[test] +fn project_many_rolls_back_on_error() { + let mut tab = TestTableau::new(3, 1e-12); + tab.h(0); + tab.cnot(0, 1); + tab.cnot(1, 2); + let before = format!("{tab}"); + + // The first projection succeeds; the second is impossible for GHZ. + let err = tab.project_many(&[0, 1], &[false, true]).unwrap_err(); + assert!(matches!( + err, + ProjectError::ZeroProbability { + addr0: 1, + outcome: true, + .. + } + )); + assert_eq!(format!("{tab}"), before); + assert!(tab.current_measurement_record().is_empty()); + + assert_eq!( + tab.project_many(&[0, 1], &[false]), + Err(ProjectError::LengthMismatch { + targets: 2, + outcomes: 1 + }) + ); + assert_eq!(format!("{tab}"), before); +} + +#[test] +fn probability_does_not_modify_state() { + let tab = magic_circuit(); + let before = format!("{tab}"); + for bits in all_bitstrings(3) { + let p = tab.probability(&[0, 1, 2], &bits).unwrap(); + assert_close(p, bitstring_prob_from_expectations(&tab, &bits), 1e-9); + } + // Marginal over a subset: P(q0 = 0) = (1 + ⟨Z₀⟩) / 2. + let z0 = tab.expectation(&word("ZII")); + assert_close( + tab.probability(&[0], &[false]).unwrap(), + 0.5 + 0.5 * z0, + 1e-9, + ); + assert_eq!(format!("{tab}"), before); + assert!(tab.current_measurement_record().is_empty()); + + // Zero-probability outcomes return 0.0; other errors propagate. + let mut ghz = TestTableau::new(3, 1e-12); + ghz.h(0); + ghz.cnot(0, 1); + ghz.cnot(1, 2); + assert_eq!(ghz.probability(&[0, 1], &[false, true]), Ok(0.0)); + assert!(matches!( + ghz.probability(&[0, 1], &[false]), + Err(ProjectError::LengthMismatch { .. }) + )); +} diff --git a/ppvm-python/src/ppvm/_core.pyi b/ppvm-python/src/ppvm/_core.pyi index 0926eb468..51ce31020 100644 --- a/ppvm-python/src/ppvm/_core.pyi +++ b/ppvm-python/src/ppvm/_core.pyi @@ -114,6 +114,8 @@ class _GeneralizedTableauBase: def __deepcopy__(self, memo: object) -> _GeneralizedTableauBase: ... def measure(self, addr0: int) -> int: ... def project(self, addr0: int, outcome: bool) -> float: ... + def project_many(self, targets: Sequence[int], outcomes: Sequence[bool]) -> float: ... + def probability(self, targets: Sequence[int], outcomes: Sequence[bool]) -> float: ... def measure_many(self, targets: Sequence[int]) -> list[int]: ... def current_measurement_record(self) -> list[int]: ... def coefficients(self) -> dict[int, complex]: ... diff --git a/ppvm-python/src/ppvm/generalized_tableau.py b/ppvm-python/src/ppvm/generalized_tableau.py index 0c7399bf4..55848f0a0 100644 --- a/ppvm-python/src/ppvm/generalized_tableau.py +++ b/ppvm-python/src/ppvm/generalized_tableau.py @@ -55,6 +55,29 @@ class MeasurementResult(enum.IntEnum): _BY_VALUE = (MeasurementResult.ZERO, MeasurementResult.ONE, MeasurementResult.LOST) +def _projection_outcome(value: MeasurementResult) -> bool: + """Convert a projection target value to the native bool outcome.""" + if value == MeasurementResult.LOST: + raise NotImplementedError("Projecting onto a LOST outcome is not implemented.") + if value not in (MeasurementResult.ZERO, MeasurementResult.ONE): + raise ValueError(f"value must be MeasurementResult.ZERO or ONE (got {value!r}).") + return value == MeasurementResult.ONE + + +def _projection_args( + targets: Iterable[int], values: Iterable[MeasurementResult] +) -> tuple[list[int], list[bool]]: + """Validate ``project_many`` / ``probability`` arguments for the native call.""" + targets = list(targets) + outcomes = [_projection_outcome(v) for v in values] + if len(targets) != len(outcomes): + raise ValueError( + f"got {len(targets)} targets but {len(outcomes)} values; " + "they must have the same length." + ) + return targets, outcomes + + @dataclass(frozen=True) class GeneralizedTableau( CliffordMixin, @@ -216,11 +239,56 @@ def project(self, addr0: int, value: MeasurementResult) -> float: ValueError: If ``value`` has zero probability (the state cannot be normalized). The state is left unchanged. """ - if value == MeasurementResult.LOST: - raise NotImplementedError("Projecting onto a LOST outcome is not implemented.") - if value not in (MeasurementResult.ZERO, MeasurementResult.ONE): - raise ValueError(f"value must be MeasurementResult.ZERO or ONE (got {value!r}).") - return self._interface.project(addr0, value == MeasurementResult.ONE) + return self._interface.project(addr0, _projection_outcome(value)) + + def project_many(self, targets: Iterable[int], values: Iterable[MeasurementResult]) -> float: + """Post-select several qubits onto Z-basis outcomes, in order. + + Equivalent to calling `project` on each ``(target, value)`` pair and + multiplying the returned probabilities, but atomic: if any projection + fails, the state and measurement record are restored to what they + were before the call. + + Args: + targets: The indices of the target qubits. + values: The outcome to project each target onto, ``ZERO`` or + ``ONE``; must have the same length as ``targets``. + + Returns: + The joint probability of measuring ``values`` on ``targets``. + + Raises: + NotImplementedError: If any value is ``LOST`` or any target qubit + has been lost. + ValueError: If ``targets`` and ``values`` differ in length, or the + outcomes have zero joint probability. The state is left + unchanged. + """ + return self._interface.project_many(*_projection_args(targets, values)) + + def probability(self, targets: Iterable[int], values: Iterable[MeasurementResult]) -> float: + """Return the joint probability of measuring ``values`` on ``targets``. + + Unlike `project_many`, this does not modify the state: it forks the + state and calls `project_many` on the fork. Outcomes with zero joint + probability return ``0.0`` rather than raising. Passing every qubit + gives the bitstring probability ``P(z) = |⟨z|ψ⟩|²``; passing a subset + gives the marginal probability over those qubits. + + Args: + targets: The indices of the target qubits. + values: The outcome for each target, ``ZERO`` or ``ONE``; must + have the same length as ``targets``. + + Returns: + The joint probability of measuring ``values`` on ``targets``. + + Raises: + NotImplementedError: If any value is ``LOST`` or any target qubit + has been lost. + ValueError: If ``targets`` and ``values`` differ in length. + """ + return self._interface.probability(*_projection_args(targets, values)) def measure_many(self, *targets: int | Iterable[int]) -> list[MeasurementResult]: """Measure several qubits in the Z basis. diff --git a/ppvm-python/test/generalized_tableau/test_project.py b/ppvm-python/test/generalized_tableau/test_project.py index d838716bd..09f352586 100644 --- a/ppvm-python/test/generalized_tableau/test_project.py +++ b/ppvm-python/test/generalized_tableau/test_project.py @@ -167,3 +167,59 @@ def test_project_lost_qubit_not_implemented(): tab.project(0, ZERO) assert tab.current_measurement_record() == [] assert tab.project(1, ZERO) == pytest.approx(1.0, abs=1e-12) + + +def test_project_many_gives_joint_probability(): + tab = _magic_circuit() + for bits in itertools.product((0, 1), repeat=3): + t = tab.fork(seed=0) + expected = _prob_by_expectation(tab, bits) + values = [MeasurementResult(b) for b in bits] + try: + p = t.project_many([0, 1, 2], values) + except ValueError: + assert expected == pytest.approx(0.0, abs=1e-9) + continue + assert p == pytest.approx(expected, abs=1e-9) + assert t.current_measurement_record() == values + + +def test_project_many_rolls_back_on_error(): + tab = GeneralizedTableau(3) + tab.h(0) + tab.cnot(0, 1) + tab.cnot(1, 2) + before = str(tab) + with pytest.raises(ValueError, match="probability"): + tab.project_many([0, 1], [ZERO, ONE]) + assert str(tab) == before + assert tab.current_measurement_record() == [] + + with pytest.raises(ValueError, match="same length"): + tab.project_many([0, 1], [ZERO]) + with pytest.raises(NotImplementedError): + tab.project_many([0], [LOST]) + assert str(tab) == before + + +def test_probability_does_not_modify_state(): + tab = _magic_circuit() + before = str(tab) + for bits in itertools.product((0, 1), repeat=3): + p = tab.probability([0, 1, 2], [MeasurementResult(b) for b in bits]) + assert p == pytest.approx(_prob_by_expectation(tab, bits), abs=1e-9) + # Marginal over a subset: P(q0 = 0) = (1 + ⟨Z₀⟩) / 2. + z0 = tab.expectation("ZII") + assert tab.probability([0], [ZERO]) == pytest.approx(0.5 + 0.5 * z0, abs=1e-9) + assert str(tab) == before + assert tab.current_measurement_record() == [] + + ghz = GeneralizedTableau(3) + ghz.h(0) + ghz.cnot(0, 1) + ghz.cnot(1, 2) + assert ghz.probability([0, 1], [ZERO, ONE]) == 0.0 + with pytest.raises(ValueError, match="same length"): + ghz.probability([0, 1], [ZERO]) + with pytest.raises(NotImplementedError): + ghz.probability([0], [LOST]) From ceea75d05eeb7f9f9e98d324cd9ff8709439ba69 Mon Sep 17 00:00:00 2001 From: Jonathan Wurtz Date: Tue, 6 Oct 2026 14:13:16 -0400 Subject: [PATCH 4/6] fix(tableau): validate project targets up front and compute case-b probability exactly Address review findings on project / project_many / probability: - Validate every target (index range, then loss) before any projection, clone, or fork. An out-of-range index no longer panics mid-batch and leaves project_many half-applied, and a lost target now errors regardless of its position in the list. Adds ProjectError::QubitOutOfRange, raised as IndexError in Python; the Python wrapper also rejects negative indices with IndexError. - In case b (Z is a stabilizer), compute the probability as the kept share of the coefficient norm instead of (1 +/- )/2, so small probabilities are accurate and only exact zeros are refused. Case a keeps PROJECT_ZERO_TOL, now documented. - probability() forks once and projects without a rollback backup. - Restore the phase debug_assert from measure in the case-b path. - Drop the duplicate Python length check; the native LengthMismatch surfaces as ValueError. Document the new API in skills/ppvm-usage/SKILL.md and the Python tableau quickstart, including error behavior, precision limits, and the pure-trajectory scope. Co-Authored-By: Claude Opus 5.5 --- .../src/interface_tableau.rs | 3 + crates/ppvm-tableau/src/measure.rs | 213 ++++++++++++------ crates/ppvm-tableau/tests/project.rs | 49 ++++ .../src/pages/quickstart/python/tableau.astro | 56 +++++ ppvm-python/src/ppvm/generalized_tableau.py | 32 ++- .../test/generalized_tableau/test_project.py | 39 ++++ skills/ppvm-usage/SKILL.md | 22 +- 7 files changed, 330 insertions(+), 84 deletions(-) diff --git a/crates/ppvm-python-native/src/interface_tableau.rs b/crates/ppvm-python-native/src/interface_tableau.rs index 38dede5ba..653dc7ad7 100644 --- a/crates/ppvm-python-native/src/interface_tableau.rs +++ b/crates/ppvm-python-native/src/interface_tableau.rs @@ -18,6 +18,9 @@ pub(crate) fn measurement_to_u8(m: Option) -> u8 { fn project_error_to_py(e: ProjectError) -> PyErr { match e { + ProjectError::QubitOutOfRange { .. } => { + pyo3::exceptions::PyIndexError::new_err(e.to_string()) + } ProjectError::QubitLost(_) => { pyo3::exceptions::PyNotImplementedError::new_err(e.to_string()) } diff --git a/crates/ppvm-tableau/src/measure.rs b/crates/ppvm-tableau/src/measure.rs index ed8601161..5e566542b 100644 --- a/crates/ppvm-tableau/src/measure.rs +++ b/crates/ppvm-tableau/src/measure.rs @@ -83,15 +83,22 @@ impl Default for MeasureScratch { } } -/// Outcome probabilities below this are treated as zero by -/// [`project`](GeneralizedTableau::project): the projected state cannot be -/// normalized, so the projection is refused. +/// When `Z` on the target is not a stabilizer of the frame, outcome +/// probabilities below this are treated as zero by +/// [`project`](GeneralizedTableau::project) and its batched variants. There the +/// probability comes from `(1 ± ⟨Z⟩)/2`, whose cancellation leaves noise of +/// order machine epsilon, and normalizing a noise-sized projection would +/// produce garbage amplitudes. (When `Z` is a stabilizer the probability is +/// computed exactly and only exact zeros are refused.) pub const PROJECT_ZERO_TOL: f64 = 1e-12; -/// Error returned by [`project`](GeneralizedTableau::project) and -/// [`project_many`](GeneralizedTableau::project_many). +/// Error returned by [`project`](GeneralizedTableau::project), +/// [`project_many`](GeneralizedTableau::project_many), and +/// [`probability`](GeneralizedTableau::probability). #[derive(Debug, Clone, PartialEq)] pub enum ProjectError { + /// The target qubit index is not less than the number of qubits. + QubitOutOfRange { addr0: usize, n_qubits: usize }, /// The target qubit is lost; projecting a lost qubit is not implemented. QubitLost(usize), /// The requested outcome has (numerically) zero probability. @@ -100,13 +107,18 @@ pub enum ProjectError { outcome: bool, prob: f64, }, - /// `project_many` was given different numbers of targets and outcomes. + /// `project_many` / `probability` got different numbers of targets and + /// outcomes. LengthMismatch { targets: usize, outcomes: usize }, } impl std::fmt::Display for ProjectError { fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { match self { + ProjectError::QubitOutOfRange { addr0, n_qubits } => write!( + f, + "cannot project qubit {addr0}: index out of range for {n_qubits} qubits" + ), ProjectError::QubitLost(addr0) => { write!(f, "cannot project qubit {addr0}: qubit is lost") } @@ -435,25 +447,137 @@ where /// /// # Errors /// + /// - [`ProjectError::QubitOutOfRange`] if `addr0 >= n_qubits`. /// - [`ProjectError::QubitLost`] if `addr0` is lost. - /// - [`ProjectError::ZeroProbability`] if the outcome's probability is - /// below [`PROJECT_ZERO_TOL`]. + /// - [`ProjectError::ZeroProbability`] if the outcome has zero probability. + /// When `Z` on `addr0` is a stabilizer of the frame, the probability is + /// computed exactly and only an exact zero is refused; otherwise it comes + /// from `⟨Z⟩` and anything below [`PROJECT_ZERO_TOL`] is refused. /// /// The state is left unchanged when an error is returned. pub fn project(&mut self, addr0: usize, outcome: bool) -> Result { + self.check_project_target(addr0)?; + self.project_unchecked(addr0, outcome) + } + + /// Post-select each `targets[k]` onto `outcomes[k]`, in order, and return + /// the joint probability of those outcomes. + /// + /// Equivalent to calling [`project`](Self::project) on each pair and + /// multiplying the returned conditional probabilities, but atomic: every + /// target is validated before any projection, and if a later outcome has + /// zero probability the state (including the measurement record) is + /// restored to what it was before the call. + /// + /// # Errors + /// + /// - [`ProjectError::LengthMismatch`] if `targets` and `outcomes` differ in + /// length. + /// - Any error from [`project`](Self::project) on one of the targets. + pub fn project_many( + &mut self, + targets: &[usize], + outcomes: &[bool], + ) -> Result { + Self::check_lengths(targets, outcomes)?; + for &addr0 in targets { + self.check_project_target(addr0)?; + } + let backup = self.clone(); + let result = self.project_sequence(targets, outcomes); + if result.is_err() { + *self = backup; + } + result + } + + /// Joint probability of measuring `outcomes` on `targets`, without + /// modifying this state. + /// + /// Forks the state and projects the fork as in + /// [`project_many`](Self::project_many). Outcomes with zero joint + /// probability return `Ok(0.0)` rather than an error. + /// + /// # Errors + /// + /// - [`ProjectError::LengthMismatch`] if `targets` and `outcomes` differ in + /// length. + /// - [`ProjectError::QubitOutOfRange`] or [`ProjectError::QubitLost`] if + /// any target is out of range or lost, regardless of target order. + pub fn probability(&self, targets: &[usize], outcomes: &[bool]) -> Result { + Self::check_lengths(targets, outcomes)?; + for &addr0 in targets { + self.check_project_target(addr0)?; + } + // Projection never draws from the RNG, so the fork's seed is irrelevant. + // The fork is discarded, so no rollback backup is needed. + match self.fork(Some(0)).project_sequence(targets, outcomes) { + Err(ProjectError::ZeroProbability { .. }) => Ok(0.0), + result => result, + } + } + + fn check_lengths(targets: &[usize], outcomes: &[bool]) -> Result<(), ProjectError> { + if targets.len() != outcomes.len() { + return Err(ProjectError::LengthMismatch { + targets: targets.len(), + outcomes: outcomes.len(), + }); + } + Ok(()) + } + + fn check_project_target(&self, addr0: usize) -> Result<(), ProjectError> { + let n_qubits = self.is_lost.len(); + if addr0 >= n_qubits { + return Err(ProjectError::QubitOutOfRange { addr0, n_qubits }); + } if self.is_lost[addr0] { return Err(ProjectError::QubitLost(addr0)); } + Ok(()) + } + + /// Project each target in order and multiply the probabilities. Targets + /// must already be validated; on error the state may be partially projected. + fn project_sequence( + &mut self, + targets: &[usize], + outcomes: &[bool], + ) -> Result { + let mut prob = 1.0; + for (&addr0, &outcome) in targets.iter().zip(outcomes) { + prob *= self.project_unchecked(addr0, outcome)?; + } + Ok(prob) + } + /// [`project`](Self::project) without the range and loss checks. Returns + /// `ZeroProbability` before mutating anything. + fn project_unchecked(&mut self, addr0: usize, outcome: bool) -> Result { let (phase_decomp, stab_anticomm_bits, destab_anticomm_bits) = self.compute_decomposition(addr0, Pauli::Z); if stab_anticomm_bits == I::zero() { - // Case b: Z is a stabilizer — self-pairing overlap, then filter. - let entries: Vec<(Complex, I)> = self.coefficients.iter().copied().collect(); - let z = Self::compute_overlap_case_b(&entries, phase_decomp, destab_anticomm_bits); - let prob = Self::outcome_probability(z, outcome); - if prob < PROJECT_ZERO_TOL { + // Case b: Z is a stabilizer, so each entry lies wholly in one + // outcome. The probability is the kept share of the norm, computed + // directly to avoid the cancellation in (1 ± ⟨Z⟩)/2. + debug_assert!( + phase_decomp == 0 || phase_decomp == 2, + "Measurement result cannot be imaginary!" + ); + let z_sign = phase_decomp == 2; + let (mut kept, mut total) = (0.0f64, 0.0f64); + for &(coeff, alpha) in self.coefficients.iter() { + let norm_sq = coeff.norm_sqr().to_f64().unwrap_or(0.0); + total += norm_sq; + let parity = symplectic_inner(alpha, destab_anticomm_bits) % 2 != 0; + if (parity ^ outcome) == z_sign { + kept += norm_sq; + } + } + let prob = if total > 0.0 { kept / total } else { 0.0 }; + if prob <= 0.0 { return Err(ProjectError::ZeroProbability { addr0, outcome, @@ -461,7 +585,10 @@ where }); } // `project_case_b` refills `self.coefficients` from `entries`. - self.coefficients = C::new(); + let entries: Vec<(Complex, I)> = + std::mem::replace(&mut self.coefficients, C::new()) + .into_iter() + .collect(); self.project_case_b(&entries, outcome, phase_decomp, destab_anticomm_bits); self.measurement_record.push(Some(outcome)); Ok(prob) @@ -505,64 +632,6 @@ where } } - /// Post-select each `targets[k]` onto `outcomes[k]`, in order, and return - /// the joint probability of those outcomes. - /// - /// Equivalent to calling [`project`](Self::project) on each pair and - /// multiplying the returned conditional probabilities, but atomic: if any - /// projection fails, the state (including the measurement record) is - /// restored to what it was before the call. - /// - /// # Errors - /// - /// - [`ProjectError::LengthMismatch`] if `targets` and `outcomes` differ in - /// length. - /// - Any error from [`project`](Self::project) on one of the targets. - pub fn project_many( - &mut self, - targets: &[usize], - outcomes: &[bool], - ) -> Result { - if targets.len() != outcomes.len() { - return Err(ProjectError::LengthMismatch { - targets: targets.len(), - outcomes: outcomes.len(), - }); - } - let backup = self.clone(); - let mut prob = 1.0; - for (&addr0, &outcome) in targets.iter().zip(outcomes) { - match self.project(addr0, outcome) { - Ok(p) => prob *= p, - Err(e) => { - *self = backup; - return Err(e); - } - } - } - Ok(prob) - } - - /// Joint probability of measuring `outcomes` on `targets`, without - /// modifying this state. - /// - /// Forks the state and calls [`project_many`](Self::project_many) on the - /// fork. Outcomes with zero joint probability return `Ok(0.0)` rather than - /// an error. - /// - /// # Errors - /// - /// - [`ProjectError::LengthMismatch`] if `targets` and `outcomes` differ in - /// length. - /// - [`ProjectError::QubitLost`] if any target qubit is lost. - pub fn probability(&self, targets: &[usize], outcomes: &[bool]) -> Result { - // Projection never draws from the RNG, so the fork's seed is irrelevant. - match self.fork(Some(0)).project_many(targets, outcomes) { - Err(ProjectError::ZeroProbability { .. }) => Ok(0.0), - result => result, - } - } - /// Probability of `outcome` given `⟨Z⟩ = z`, clamped to `[0, 1]`. fn outcome_probability(z: f64, outcome: bool) -> f64 { let prob = if outcome { diff --git a/crates/ppvm-tableau/tests/project.rs b/crates/ppvm-tableau/tests/project.rs index 5d308ce3a..5121ea9b1 100644 --- a/crates/ppvm-tableau/tests/project.rs +++ b/crates/ppvm-tableau/tests/project.rs @@ -376,3 +376,52 @@ fn probability_does_not_modify_state() { Err(ProjectError::LengthMismatch { .. }) )); } + +#[test] +fn project_out_of_range_is_atomic() { + let mut tab = TestTableau::new(2, 1e-12); + tab.h(0); + let before = format!("{tab}"); + let oor = Err(ProjectError::QubitOutOfRange { + addr0: 9, + n_qubits: 2, + }); + assert_eq!(tab.project(9, false), oor); + // Target 0 is valid, but nothing is projected because target 9 is checked first. + assert_eq!(tab.project_many(&[0, 9], &[false, false]), oor); + assert_eq!(tab.probability(&[0, 9], &[false, false]), oor); + assert_eq!(format!("{tab}"), before); + assert!(tab.current_measurement_record().is_empty()); +} + +#[test] +fn project_lost_target_errors_regardless_of_order() { + // Qubit 0 is |0⟩ (so outcome 1 is impossible) and qubit 1 is lost. + let mut tab = TestTableau::new(2, 1e-12); + tab.loss_channel(1, 1.0); + let lost = Err(ProjectError::QubitLost(1)); + assert_eq!(tab.probability(&[0, 1], &[true, false]), lost); + assert_eq!(tab.probability(&[1, 0], &[false, true]), lost); + assert_eq!(tab.project_many(&[0, 1], &[true, false]), lost); +} + +#[test] +fn project_small_probability_is_exact_when_z_is_stabilizer() { + // RX(θ)|0⟩ keeps Z as a stabilizer of the frame (case b), and + // P(1) = sin²(θ/2) ≈ 2.5e-13 is below PROJECT_ZERO_TOL. + let theta: f64 = 1e-6; + let expected = (theta / 2.0).sin().powi(2); + assert!(expected < PROJECT_ZERO_TOL); + + let mut tab = TestTableau::new(1, 1e-12); + tab.rx(0, theta); + assert!(!is_case_a(&tab, 0)); + let p = tab.probability(&[0], &[true]).unwrap(); + assert!( + (p - expected).abs() < 1e-9 * expected, + "expected {expected}, got {p}" + ); + let p = tab.project(0, true).unwrap(); + assert!((p - expected).abs() < 1e-9 * expected); + assert_close(tab.expectation(&word("Z")), -1.0, 1e-9); +} diff --git a/docs/src/pages/quickstart/python/tableau.astro b/docs/src/pages/quickstart/python/tableau.astro index 8c1e667be..8177bf7c5 100644 --- a/docs/src/pages/quickstart/python/tableau.astro +++ b/docs/src/pages/quickstart/python/tableau.astro @@ -8,6 +8,7 @@ const toc = [ { label: "Overview", href: "#overview" }, { label: "Drive it from Python", href: "#basic-usage" }, { label: "Forking for multi-shot", href: "#forking" }, + { label: "Post-selection", href: "#postselect" }, { label: "Stim circuit support", href: "#stim" }, { label: "Noise and loss", href: "#noise" }, { label: "Coefficient pruning", href: "#pruning" }, @@ -24,6 +25,19 @@ for shot in range(100): t = tab.fork(seed=shot) print(t.measure(0), t.measure(1))`; +const postselect = `from ppvm.generalized_tableau import MeasurementResult as M + +tab = GeneralizedTableau(n_qubits=2) +tab.h(0) +tab.cnot(0, 1) + +print(tab.probability([0, 1], [M.ZERO, M.ZERO])) # 0.5, state unchanged +print(tab.probability([0, 1], [M.ZERO, M.ONE])) # 0.0 +print(tab.probability([1], [M.ONE])) # 0.5, marginal on qubit 1 + +p = tab.project(0, M.ONE) # 0.5; qubit 0 is now |1> +print(tab.measure(1) == M.ONE) # True (Bell correlation)`; + const stimUsage = `from ppvm import GeneralizedTableau, StimProgram, sample_stim prog = StimProgram.parse(""" @@ -133,6 +147,48 @@ const sumPruning = `tab = GeneralizedTableauSum(n_qubits=4, min_abs_coeff=1e-10,

+
+

Post-selection and bitstring probabilities

+

+ measure draws its outcome from the RNG. To force the + outcome instead, use project. It projects the qubit onto + the requested Z-basis value, records that value in the measurement + record, and returns the outcome's probability. This is non-physical, + but by the chain rule it computes bitstring and marginal + probabilities, P(z) = |⟨z|ψ⟩|², directly instead of + estimating them from samples: +

+
{postselect}
+

+ project_many projects several qubits in order and returns + their joint probability. If any projection fails, the state is rolled + back to what it was before the call. probability takes the + same arguments but runs on a fork, so the state is never modified. + project and project_many raise + ValueError for an outcome with zero probability, leaving + the state unchanged; probability returns 0.0 + instead. All three raise NotImplementedError for + LOST values or lost qubits, and IndexError + for out-of-range targets. Every target is checked before any is + projected. +

+

+ When Z on a target is a stabilizer of the tableau frame, + the probability is computed exactly from the coefficient weights. + Otherwise it comes from (1 ± ⟨Z⟩)/2, whose rounding error + grows with the number of coefficients, and values below about + 1e-12 are reported as zero. In both cases the result + reflects any coefficients already pruned by + min_abs_coeff. +

+

+ These methods act on the current trajectory's pure state: noise + channels applied earlier have already sampled their outcomes, so the + result is the probability for that trajectory. Average over + trajectories to estimate it under noise. +

+
+

Stim circuit support

diff --git a/ppvm-python/src/ppvm/generalized_tableau.py b/ppvm-python/src/ppvm/generalized_tableau.py index 55848f0a0..00fb83018 100644 --- a/ppvm-python/src/ppvm/generalized_tableau.py +++ b/ppvm-python/src/ppvm/generalized_tableau.py @@ -67,15 +67,21 @@ def _projection_outcome(value: MeasurementResult) -> bool: def _projection_args( targets: Iterable[int], values: Iterable[MeasurementResult] ) -> tuple[list[int], list[bool]]: - """Validate ``project_many`` / ``probability`` arguments for the native call.""" + """Convert ``project_many`` / ``probability`` arguments for the native call. + + Length mismatches and indices ``>= n_qubits`` are checked natively. + """ targets = list(targets) - outcomes = [_projection_outcome(v) for v in values] - if len(targets) != len(outcomes): - raise ValueError( - f"got {len(targets)} targets but {len(outcomes)} values; " - "they must have the same length." - ) - return targets, outcomes + for addr0 in targets: + _check_nonnegative_index(addr0) + return targets, [_projection_outcome(v) for v in values] + + +def _check_nonnegative_index(addr0: int) -> None: + # The native layer takes unsigned indices; a negative int would surface as + # OverflowError instead of the documented IndexError. + if addr0 < 0: + raise IndexError(f"qubit index must be non-negative (got {addr0}).") @dataclass(frozen=True) @@ -238,16 +244,18 @@ def project(self, addr0: int, value: MeasurementResult) -> float: been lost; projection with loss is not implemented. ValueError: If ``value`` has zero probability (the state cannot be normalized). The state is left unchanged. + IndexError: If ``addr0`` is not a valid qubit index. """ + _check_nonnegative_index(addr0) return self._interface.project(addr0, _projection_outcome(value)) def project_many(self, targets: Iterable[int], values: Iterable[MeasurementResult]) -> float: """Post-select several qubits onto Z-basis outcomes, in order. Equivalent to calling `project` on each ``(target, value)`` pair and - multiplying the returned probabilities, but atomic: if any projection - fails, the state and measurement record are restored to what they - were before the call. + multiplying the returned probabilities, but atomic: every target is + validated first, and if any projection fails, the state and + measurement record are restored to what they were before the call. Args: targets: The indices of the target qubits. @@ -263,6 +271,7 @@ def project_many(self, targets: Iterable[int], values: Iterable[MeasurementResul ValueError: If ``targets`` and ``values`` differ in length, or the outcomes have zero joint probability. The state is left unchanged. + IndexError: If any target is not a valid qubit index. """ return self._interface.project_many(*_projection_args(targets, values)) @@ -287,6 +296,7 @@ def probability(self, targets: Iterable[int], values: Iterable[MeasurementResult NotImplementedError: If any value is ``LOST`` or any target qubit has been lost. ValueError: If ``targets`` and ``values`` differ in length. + IndexError: If any target is not a valid qubit index. """ return self._interface.probability(*_projection_args(targets, values)) diff --git a/ppvm-python/test/generalized_tableau/test_project.py b/ppvm-python/test/generalized_tableau/test_project.py index 09f352586..09ff73382 100644 --- a/ppvm-python/test/generalized_tableau/test_project.py +++ b/ppvm-python/test/generalized_tableau/test_project.py @@ -223,3 +223,42 @@ def test_probability_does_not_modify_state(): ghz.probability([0, 1], [ZERO]) with pytest.raises(NotImplementedError): ghz.probability([0], [LOST]) + + +def test_project_out_of_range_is_atomic(): + tab = GeneralizedTableau(2) + tab.h(0) + before = str(tab) + with pytest.raises(IndexError): + tab.project(9, ZERO) + with pytest.raises(IndexError): + tab.project_many([0, 9], [ZERO, ZERO]) + with pytest.raises(IndexError): + tab.probability([0, 9], [ZERO, ZERO]) + with pytest.raises(IndexError): + tab.project(-1, ZERO) + with pytest.raises(IndexError): + tab.project_many([0, -1], [ZERO, ZERO]) + with pytest.raises(IndexError): + tab.probability([-1], [ZERO]) + assert str(tab) == before + assert tab.current_measurement_record() == [] + + +def test_lost_target_errors_regardless_of_order(): + tab = GeneralizedTableau(2) + tab.loss_channel(1, 1.0) + for targets, values in [([0, 1], [ONE, ZERO]), ([1, 0], [ZERO, ONE])]: + with pytest.raises(NotImplementedError): + tab.probability(targets, values) + with pytest.raises(NotImplementedError): + tab.project_many(targets, values) + + +def test_small_probability_is_exact(): + theta = 1e-6 + expected = math.sin(theta / 2) ** 2 # ≈ 2.5e-13 + tab = GeneralizedTableau(1, min_abs_coeff=1e-12) + tab.rx(0, theta) + assert tab.probability([0], [ONE]) == pytest.approx(expected, rel=1e-9) + assert tab.project(0, ONE) == pytest.approx(expected, rel=1e-9) diff --git a/skills/ppvm-usage/SKILL.md b/skills/ppvm-usage/SKILL.md index 4c290a72f..e7a315827 100644 --- a/skills/ppvm-usage/SKILL.md +++ b/skills/ppvm-usage/SKILL.md @@ -134,6 +134,21 @@ r1 = tab.measure(1) # correlated with r0 (Bell state) For throughput on tableau circuits, batch same-gate layers in one Python call. Single-qubit gates accept variadic targets or a sequence: `tab.h(0, 2, 4)` and `tab.h([0, 2, 4])` are equivalent. Two-qubit gates consume a flat target list as consecutive pairs: `tab.cnot([0, 1, 2, 3])` applies `(0, 1)` and `(2, 3)`. Rotations and Pauli/depolarizing noise use the same convention with `theta=...` or `p=...`. This avoids one Python→Rust call per target and forwards to fused Rust tableau kernels internally. `measure` stays scalar; use `tab.measure_many([0, 1, 2])` for readout layers. +To post-select instead of sampling, use `project`, `project_many`, or `probability`. They force the Z-basis outcome rather than drawing it from the RNG, and return its probability. This is non-physical, but by the chain rule it computes bitstring and marginal probabilities, P(z) = |⟨z|ψ⟩|², directly instead of estimating them from samples: + +```python +from ppvm.generalized_tableau import MeasurementResult as M + +tab = GeneralizedTableau(n_qubits=2) +tab.h(0) +tab.cnot(0, 1) +p = tab.probability([0, 1], [M.ZERO, M.ZERO]) # P(00) = 0.5; the state is unchanged +p = tab.project_many([0, 1], [M.ONE, M.ONE]) # P(11) = 0.5; collapses to |11>, rolls back on error +p = tab.project(0, M.ONE) # 1.0: projects one qubit, returns its probability +``` + +`project` / `project_many` raise `ValueError` for an outcome with zero probability, leaving the state unchanged; `probability` returns `0.0` instead. All three raise `NotImplementedError` for `LOST` values or lost qubits, and `IndexError` for out-of-range targets, checking every target before projecting any. Probabilities below about 1e-12 can be reported as zero when `Z` on the target is not a stabilizer of the tableau frame. They act on the current trajectory's pure state: noise channels applied earlier have already been sampled. + Non-Clifford gates and Stim programs: ```python @@ -212,6 +227,8 @@ let outcome = tab.measure(0); For layer-style tableau circuits in Rust, prefer explicit batch methods instead of per-target loops: `tab.h_many(&[0, 2, 4])`, `tab.cnot_many(&[(0, 1), (2, 3)])`, `tab.rx_many(&targets, theta)`, `tab.rzz_many(&pairs, theta)`, `tab.depolarize1_many(&targets, p)`, `tab.measure_many(&targets)`, `tab.reset_many(&targets)`, and the analogous `*_many` forms. `GeneralizedTableau` specializes these into fused bit operations. Other backends may expose trait-default `*_many` methods too, but the fused speedup is tableau-specific. +For post-selection in Rust, use `tab.project(q, outcome)`, `tab.project_many(&targets, &outcomes)`, and `tab.probability(&targets, &outcomes)`, with `bool` outcomes. Each returns `Result`. These are not fused kernels: they project one qubit at a time, and `project_many` / `probability` copy the state once per call. + Pick `IndexType` by qubit count: `usize` up to ~64, `u128` up to 128, `bnum::types::U256` / `U512` / `U1024` beyond. **Using `usize` past 64 qubits silently overflows** — this is the second-most-common bug after Heisenberg-order mistakes. ### Running Stim programs (Rust) @@ -238,7 +255,7 @@ With the `rayon` feature, `sample` fans shots across the global thread pool (ser Availability varies by backend and language binding. Don't trust intuition: the Python `PauliSum` exposes a deliberately narrower surface than the Rust `PauliSum` or the `GeneralizedTableau`. Names: everything is `snake_case`; daggers are `_dag` (e.g. `s_dag`, `sqrt_x_dag`, `t_dag`) on both Rust and Python. -In the tables below: **R** = Rust on both backends, **P-S** = Python `PauliSum` / `LossyPauliSum`, **P-T** = Python `GeneralizedTableau`. A check means the method is exposed there. +In the tables below: **R** = Rust on both backends, **P-S** = Python `PauliSum` / `LossyPauliSum`, **P-T** = Python `GeneralizedTableau`. A check means the method is exposed there. "tableau" in the **R** column means Rust `GeneralizedTableau` only, not Rust `PauliSum`. ### Clifford gates @@ -267,6 +284,9 @@ Important: the six off-diagonal two-qubit rotations (`rxy`, `rxz`, `ryx`, `ryz`, | Method | R | P-S | P-T | |-----------------------------------------------------------------------------|---|-----|-----| | `measure(q)` → `MeasurementResult` (`ZERO`/`ONE`/`LOST` on tableau) | ✓ | — | ✓ | +| `project(q, value)` → probability (post-select; mutates the state) | tableau | — | ✓ | +| `project_many(targets, values)` → joint probability (atomic on error) | tableau | — | ✓ | +| `probability(targets, values)` → joint probability (state unchanged) | tableau | — | ✓ | | `reset(q)` | ✓ | — | ✓ | | `depolarize1(q, p=...)` | ✓ | ✓ | ✓ | | `depolarize2(q0, q1, p=...)` | ✓ | — | — | From 5a4859dac109a64b82383525403cabbd76a87991 Mon Sep 17 00:00:00 2001 From: Jonathan Wurtz Date: Tue, 6 Oct 2026 14:55:48 -0400 Subject: [PATCH 5/6] fix(tableau): normalize case-a projection probability for pruned states Gate branching prunes coefficients below the threshold without renormalizing, so the coefficient norm can drop below 1. The case-a path of project computed (1 +/- z)/2 from the raw overlap z = N, giving wrong conditional probabilities (a deterministic outcome could come back below 1). Divide the overlap by the coefficient norm first. Case b already computed kept/total and was unaffected. Add a 3-qubit regression test with aggressive pruning (min_abs_coeff = 0.25, norm^2 ~ 0.96) where qubit 0 is |0> on the case-a path; it reported P(0) = 0.980 before this fix. Reported by Copilot review on #240. Co-Authored-By: Claude Opus 5.5 --- crates/ppvm-tableau/src/measure.rs | 14 ++++++++++++- crates/ppvm-tableau/tests/project.rs | 30 ++++++++++++++++++++++++++++ 2 files changed, 43 insertions(+), 1 deletion(-) diff --git a/crates/ppvm-tableau/src/measure.rs b/crates/ppvm-tableau/src/measure.rs index 5e566542b..d8f943592 100644 --- a/crates/ppvm-tableau/src/measure.rs +++ b/crates/ppvm-tableau/src/measure.rs @@ -597,13 +597,25 @@ where let odd_phase_mask = self.odd_phase_destabilizer_mask(); let coeff_map: HashMap> = self.coefficients.iter().map(|&(c, i)| (i, c)).collect(); - let z = Self::compute_overlap_case_a( + let overlap = Self::compute_overlap_case_a( &coeff_map, phase_decomp, destab_anticomm_bits, stab_anticomm_bits, odd_phase_mask, ); + // Gate branching prunes small coefficients without renormalizing, + // so the overlap is ⟨ψ|Z|ψ⟩ for an unnormalized ψ. Divide by the + // norm to get ⟨Z⟩ before taking (1 ± ⟨Z⟩)/2. + let norm_sq: f64 = coeff_map + .values() + .map(|c| c.norm_sqr().to_f64().unwrap_or(0.0)) + .sum(); + let z = if norm_sq > 0.0 { + overlap / norm_sq + } else { + 0.0 + }; let prob = Self::outcome_probability(z, outcome); if prob < PROJECT_ZERO_TOL { return Err(ProjectError::ZeroProbability { diff --git a/crates/ppvm-tableau/tests/project.rs b/crates/ppvm-tableau/tests/project.rs index 5121ea9b1..61b8e8340 100644 --- a/crates/ppvm-tableau/tests/project.rs +++ b/crates/ppvm-tableau/tests/project.rs @@ -425,3 +425,33 @@ fn project_small_probability_is_exact_when_z_is_stabilizer() { assert!((p - expected).abs() < 1e-9 * expected); assert_close(tab.expectation(&word("Z")), -1.0, 1e-9); } + +#[test] +fn project_case_a_normalizes_pruned_state() { + // Gate branching prunes coefficients below the threshold without + // renormalizing, so the coefficient norm can drop below 1. Here RY(0.4) on + // qubit 1 prunes its sin(0.2) ≈ 0.199 branch (threshold 0.25), leaving + // norm² = cos²(0.2) ≈ 0.96. Qubit 0 is then |0⟩ on the case-a path, so + // P(0) must be exactly 1 and P(1) exactly 0 regardless of that norm. + let mut tab = TestTableau::new(3, 0.25); + tab.ry(1, 0.4); + tab.h(2); + tab.cz(1, 2); + tab.h(0); + tab.ry(0, -FRAC_PI_2); + let norm_sq: f64 = tab.coefficients.iter().map(|(c, _)| c.norm_sqr()).sum(); + assert!( + norm_sq < 0.99, + "test needs a pruned, unnormalized state: {norm_sq}" + ); + assert!(is_case_a(&tab, 0)); + + assert_close(tab.probability(&[0], &[false]).unwrap(), 1.0, 1e-12); + assert_eq!(tab.probability(&[0], &[true]), Ok(0.0)); + let total: f64 = all_bitstrings(3) + .iter() + .map(|bits| tab.probability(&[0, 1, 2], bits).unwrap()) + .sum(); + assert_close(total, 1.0, 1e-12); + assert_close(tab.project(0, false).unwrap(), 1.0, 1e-12); +} From 0b61e8693ba8bf0840e83a5ba9730130f4c0e786 Mon Sep 17 00:00:00 2001 From: Jonathan Wurtz Date: Tue, 6 Oct 2026 15:03:06 -0400 Subject: [PATCH 6/6] fix(python): raise IndexError for any out-of-range projection target A target larger than usize::MAX passed the Python non-negativity check and failed during PyO3 conversion with OverflowError, contradicting the documented IndexError. Validate 0 <= target < n_qubits in the Python wrapper before calling the native layer, covering negative, oversized, and arbitrarily large integers. The native range check remains for Rust callers. Reported by Copilot review on #240. Co-Authored-By: Claude Opus 5.5 --- ppvm-python/src/ppvm/generalized_tableau.py | 25 ++++++++++--------- .../test/generalized_tableau/test_project.py | 5 ++++ 2 files changed, 18 insertions(+), 12 deletions(-) diff --git a/ppvm-python/src/ppvm/generalized_tableau.py b/ppvm-python/src/ppvm/generalized_tableau.py index 00fb83018..c6815a584 100644 --- a/ppvm-python/src/ppvm/generalized_tableau.py +++ b/ppvm-python/src/ppvm/generalized_tableau.py @@ -65,23 +65,24 @@ def _projection_outcome(value: MeasurementResult) -> bool: def _projection_args( - targets: Iterable[int], values: Iterable[MeasurementResult] + targets: Iterable[int], values: Iterable[MeasurementResult], n_qubits: int ) -> tuple[list[int], list[bool]]: - """Convert ``project_many`` / ``probability`` arguments for the native call. + """Validate and convert ``project_many`` / ``probability`` arguments. - Length mismatches and indices ``>= n_qubits`` are checked natively. + Length mismatches are checked natively. """ targets = list(targets) for addr0 in targets: - _check_nonnegative_index(addr0) + _check_qubit_index(addr0, n_qubits) return targets, [_projection_outcome(v) for v in values] -def _check_nonnegative_index(addr0: int) -> None: - # The native layer takes unsigned indices; a negative int would surface as - # OverflowError instead of the documented IndexError. - if addr0 < 0: - raise IndexError(f"qubit index must be non-negative (got {addr0}).") +def _check_qubit_index(addr0: int, n_qubits: int) -> None: + # Check here rather than relying on the native check: the native layer + # takes unsigned machine-sized indices, so a negative or very large int + # would surface as OverflowError instead of the documented IndexError. + if not 0 <= addr0 < n_qubits: + raise IndexError(f"qubit index {addr0} is out of range for {n_qubits} qubits.") @dataclass(frozen=True) @@ -246,7 +247,7 @@ def project(self, addr0: int, value: MeasurementResult) -> float: be normalized). The state is left unchanged. IndexError: If ``addr0`` is not a valid qubit index. """ - _check_nonnegative_index(addr0) + _check_qubit_index(addr0, self.n_qubits) return self._interface.project(addr0, _projection_outcome(value)) def project_many(self, targets: Iterable[int], values: Iterable[MeasurementResult]) -> float: @@ -273,7 +274,7 @@ def project_many(self, targets: Iterable[int], values: Iterable[MeasurementResul unchanged. IndexError: If any target is not a valid qubit index. """ - return self._interface.project_many(*_projection_args(targets, values)) + return self._interface.project_many(*_projection_args(targets, values, self.n_qubits)) def probability(self, targets: Iterable[int], values: Iterable[MeasurementResult]) -> float: """Return the joint probability of measuring ``values`` on ``targets``. @@ -298,7 +299,7 @@ def probability(self, targets: Iterable[int], values: Iterable[MeasurementResult ValueError: If ``targets`` and ``values`` differ in length. IndexError: If any target is not a valid qubit index. """ - return self._interface.probability(*_projection_args(targets, values)) + return self._interface.probability(*_projection_args(targets, values, self.n_qubits)) def measure_many(self, *targets: int | Iterable[int]) -> list[MeasurementResult]: """Measure several qubits in the Z basis. diff --git a/ppvm-python/test/generalized_tableau/test_project.py b/ppvm-python/test/generalized_tableau/test_project.py index 09ff73382..e0ac3a128 100644 --- a/ppvm-python/test/generalized_tableau/test_project.py +++ b/ppvm-python/test/generalized_tableau/test_project.py @@ -241,6 +241,11 @@ def test_project_out_of_range_is_atomic(): tab.project_many([0, -1], [ZERO, ZERO]) with pytest.raises(IndexError): tab.probability([-1], [ZERO]) + # Larger than usize::MAX: must still be IndexError, not OverflowError. + with pytest.raises(IndexError): + tab.project(2**64, ZERO) + with pytest.raises(IndexError): + tab.probability([0, 2**64], [ZERO, ZERO]) assert str(tab) == before assert tab.current_measurement_record() == []