diff --git a/crates/ppvm-python-native/src/interface_tableau.rs b/crates/ppvm-python-native/src/interface_tableau.rs index 268986c67..f80705052 100644 --- a/crates/ppvm-python-native/src/interface_tableau.rs +++ b/crates/ppvm-python-native/src/interface_tableau.rs @@ -55,6 +55,14 @@ macro_rules! create_interface { .collect() } + /// Post-selected Z measurement; `ProjectError` becomes `ValueError`. + #[pyo3(signature = (addr0, outcome, tol = 1e-12))] + pub fn project(&mut self, addr0: usize, outcome: bool, tol: f64) -> PyResult { + self.inner + .project(addr0, outcome, tol) + .map_err(|e| pyo3::exceptions::PyValueError::new_err(e.to_string())) + } + pub fn current_measurement_record(&self) -> Vec { self.inner .current_measurement_record() diff --git a/crates/ppvm-tableau/src/measure.rs b/crates/ppvm-tableau/src/measure.rs index 79e5c4d38..30a8104b5 100644 --- a/crates/ppvm-tableau/src/measure.rs +++ b/crates/ppvm-tableau/src/measure.rs @@ -83,6 +83,43 @@ impl Default for MeasureScratch { } } +/// How [`GeneralizedTableau::collapse_with_scratch`] picks the outcome. +#[derive(Clone, Copy, Debug)] +enum OutcomeChoice { + /// Draw the outcome from the tableau's RNG (ordinary measurement). + Sample, + /// Use the given outcome (post-selection). Rejected if its probability + /// is below `tol`. + Forced { outcome: bool, tol: f64 }, +} + +/// Why [`GeneralizedTableau::project`] refused to project. The tableau is +/// left as it was before the call. +#[derive(Clone, Copy, Debug, PartialEq)] +pub enum ProjectError { + /// The qubit has been lost, so there is no Z outcome to project onto. + LostQubit { addr0: usize }, + /// The requested outcome has probability below `tol`. + ImpossibleOutcome { outcome: bool, prob: f64, tol: f64 }, +} + +impl std::fmt::Display for ProjectError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Self::LostQubit { addr0 } => { + write!(f, "cannot project qubit {addr0}: the qubit has been lost") + } + Self::ImpossibleOutcome { outcome, prob, tol } => write!( + f, + "cannot project onto outcome {}: its probability {prob:e} is below tol = {tol:e}", + *outcome as u8 + ), + } + } +} + +impl std::error::Error for ProjectError {} + impl, I>> LossyMeasure for GeneralizedTableau where @@ -169,6 +206,31 @@ where stab_anticomm_bits: I, destab_anticomm_bits: I, ) -> Option { + let (outcome, _) = self + .collapse_with_scratch( + addr0, + scratch, + phase_decomp, + stab_anticomm_bits, + destab_anticomm_bits, + OutcomeChoice::Sample, + ) + .expect("sampling an outcome cannot fail"); + Some(outcome) + } + + /// Shared body of `measure` and [`project`](Self::project): collapse + /// `addr0` onto the outcome picked by `choice`, record it, and return it + /// with its probability. A rejected `Forced` outcome mutates nothing. + fn collapse_with_scratch( + &mut self, + addr0: usize, + scratch: &mut MeasureScratch, + phase_decomp: u8, + stab_anticomm_bits: I, + destab_anticomm_bits: I, + choice: OutcomeChoice, + ) -> Result<(bool, f64), ProjectError> { if stab_anticomm_bits == I::zero() { // Case b (fast path): Z is already a stabilizer. Overlap + filter in place. @@ -194,7 +256,7 @@ where // Gate-time truncation drops branches without renormalizing, so // divide by ⟨ψ|ψ⟩ rather than assuming it is 1. let prob_1 = 0.5 - 0.5 * z_overlap_re / norm_sqr; - let outcome = self.tableau.rng.random::() < prob_1; + let (outcome, prob) = self.choose_outcome(prob_1, choice)?; debug_assert!( phase_decomp == 0 || phase_decomp == 2, @@ -213,7 +275,7 @@ where // Case-b doesn't mutate destabilizers, so the cached mask remains valid. self.measurement_record.push(Some(outcome)); - Some(outcome) + Ok((outcome, prob)) } else { // Case a: Z is not a stabilizer — sort-merge instead of HashMap. let mut by_idx: Vec<(I, Complex)> = @@ -301,7 +363,16 @@ where .map(|(_, c)| c.norm_sqr().to_f64().unwrap_or(0.0)) .sum(); let prob_1 = 0.5 - 0.5 * z_overlap_re / norm_sqr; - let outcome = self.tableau.rng.random::() < prob_1; + let (outcome, prob) = match self.choose_outcome(prob_1, choice) { + Ok(chosen) => chosen, + Err(err) => { + // Restore the coefficients drained into `by_idx` (possibly re-sorted). + for (idx, coeff) in by_idx { + self.coefficients.unsafe_insert(idx, coeff); + } + return Err(err); + } + }; // PROJECTION: partition A (k-bit=0) and B (k-bit=1), transform B, merge. let q_idx = stab_anticomm_bits.trailing_zeros() as usize; @@ -379,10 +450,29 @@ where .update_tableau_according_to_outcome(addr0, q_idx, outcome); scratch.odd_phase_mask = None; self.measurement_record.push(Some(outcome)); - Some(outcome) + Ok((outcome, prob)) } } + /// Pick the outcome given `prob_1 = P(1)`; returns it with its probability. + fn choose_outcome( + &mut self, + prob_1: f64, + choice: OutcomeChoice, + ) -> Result<(bool, f64), ProjectError> { + let outcome = match choice { + OutcomeChoice::Sample => self.tableau.rng.random::() < prob_1, + OutcomeChoice::Forced { outcome, .. } => outcome, + }; + let prob = if outcome { prob_1 } else { 1.0 - prob_1 }; + if let OutcomeChoice::Forced { tol, .. } = choice + && (prob.is_nan() || prob < tol) + { + return Err(ProjectError::ImpossibleOutcome { outcome, prob, tol }); + } + Ok((outcome, prob)) + } + pub fn project_case_a( &mut self, outcome: bool, @@ -551,6 +641,34 @@ where bit } } + + /// Post-selected Z measurement: collapse qubit `addr0` onto `outcome`, + /// record it like [`measure`](LossyMeasure::measure), and return its + /// probability. Chaining it gives `P(z) = Π_i P(z_i | z_ Result { + debug_assert!( + tol > 0.0 && tol.is_finite(), + "tol must be positive and finite, got {tol}" + ); + if self.qubit_status[addr0] == QubitStatus::Lost { + return Err(ProjectError::LostQubit { addr0 }); + } + let (phase_decomp, stab_anticomm_bits, destab_anticomm_bits) = + self.compute_decomposition(addr0, Pauli::Z); + let mut scratch = MeasureScratch::new(); + let (_, prob) = self.collapse_with_scratch( + addr0, + &mut scratch, + phase_decomp, + stab_anticomm_bits, + destab_anticomm_bits, + OutcomeChoice::Forced { outcome, tol }, + )?; + Ok(prob) + } } /// Measurement overlap helper functions, with optional rayon parallelism. diff --git a/crates/ppvm-tableau/tests/project.rs b/crates/ppvm-tableau/tests/project.rs new file mode 100644 index 000000000..c13cfeff4 --- /dev/null +++ b/crates/ppvm-tableau/tests/project.rs @@ -0,0 +1,214 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Tests for `GeneralizedTableau::project` (post-selected Z measurement). + +use std::f64::consts::FRAC_PI_2; + +use ppvm_pauli_sum::config::dashmap::ByteFxHashF64; +use ppvm_tableau::measure::ProjectError; +use ppvm_tableau::prelude::*; + +type Tab = GeneralizedTableau, usize>; + +const TOL: f64 = 1e-12; + +/// Full observable snapshot of a tableau: stabilizer frame (via `Display`), +/// coefficients sorted by index, and the measurement record. +type Snapshot = (String, Vec<(usize, f64, f64)>, Vec>); + +fn snapshot(tab: &Tab) -> Snapshot { + let mut coeffs: Vec<_> = tab + .coefficients + .iter() + .map(|(c, i)| (*i, c.re, c.im)) + .collect(); + coeffs.sort_by_key(|e| e.0); + ( + format!("{}", tab.tableau), + coeffs, + tab.measurement_record.clone(), + ) +} + +/// `ry(θ)|0⟩`, with `Z_0` either in the stabilizer group or anticommuting +/// with a stabilizer (via `H`); measurement collapses these differently. +fn ry_state(theta: f64, z_anticommutes: bool) -> Tab { + let mut tab = Tab::new_with_seed(1, 1e-12, 0); + if z_anticommutes { + // H|0⟩ = ry(π/2)|0⟩, so H then ry(θ − π/2) is ry(θ)|0⟩. + tab.h(0); + tab.ry(0, theta - FRAC_PI_2); + } else { + tab.ry(0, theta); + } + let (_, stab_anticomm_bits, _) = tab.compute_decomposition(0, Pauli::Z); + assert_eq!(stab_anticomm_bits != 0, z_anticommutes); + tab +} + +#[test] +fn project_returns_born_probability_and_collapses() { + for z_anticommutes in [false, true] { + for theta in [0.3f64, 1.0, 2.0, 2.9] { + let p1 = (theta / 2.0).sin().powi(2); + + let mut tab = ry_state(theta, z_anticommutes); + let p = tab.project(0, true, TOL).unwrap(); + assert!((p - p1).abs() < 1e-10, "θ={theta}: {p} vs {p1}"); + assert!((tab.z_expectation(0) + 1.0).abs() < 1e-10); + + let mut tab = ry_state(theta, z_anticommutes); + let p = tab.project(0, false, TOL).unwrap(); + assert!((p - (1.0 - p1)).abs() < 1e-10); + assert!((tab.z_expectation(0) - 1.0).abs() < 1e-10); + } + } +} + +#[test] +fn project_matches_measure_that_landed_on_same_outcome() { + // A post-selected projection must leave exactly the state that a sampled + // measurement with the same outcome leaves. + let mut base = Tab::new_with_seed(3, 1e-12, 7); + base.h(0); + base.t(0); + base.cnot(0, 1); + base.ry(1, 0.7); + base.h(2); + base.t(2); + base.cz(1, 2); + + for q in 0..3 { + for seed in 0..20 { + let mut measured = base.fork(Some(seed)); + let outcome = measured.measure(q).unwrap(); + let mut projected = base.fork(Some(seed)); + projected.project(q, outcome, TOL).unwrap(); + assert_eq!( + snapshot(&measured), + snapshot(&projected), + "q={q} seed={seed}" + ); + } + } +} + +#[test] +fn ghz_correlations_and_impossible_outcome() { + let mut tab = Tab::new_with_seed(3, 1e-12, 0); + tab.h(0); + tab.t(0); + tab.cnot(0, 1); + tab.cnot(1, 2); + + assert!((tab.project(0, true, TOL).unwrap() - 0.5).abs() < 1e-12); + assert!((tab.project(1, true, TOL).unwrap() - 1.0).abs() < 1e-12); + + assert!(matches!( + tab.project(2, false, TOL), + Err(ProjectError::ImpossibleOutcome { .. }) + )); + assert!((tab.project(2, true, TOL).unwrap() - 1.0).abs() < 1e-12); + assert_eq!( + tab.measurement_record, + vec![Some(true), Some(true), Some(true)] + ); +} + +#[test] +fn impossible_outcome_is_rejected_without_mutation() { + // |0⟩ in an `H`-rotated frame, where collapse rebuilds the coefficient + // vector; a rejection must hand back the original one. + let mut tab = ry_state(0.0, true); + tab.t(0); + let before = snapshot(&tab); + assert!(matches!( + tab.project(0, true, TOL), + Err(ProjectError::ImpossibleOutcome { .. }) + )); + assert_eq!(snapshot(&tab), before); + assert!((tab.project(0, false, TOL).unwrap() - 1.0).abs() < 1e-10); +} + +#[test] +fn tol_threshold_is_respected() { + // P(1) = sin²(θ/2) = 1e-8: accepted with tol = 1e-12, rejected with 1e-6. + let theta = 2.0 * 1e-4f64.asin(); + let mut tab = ry_state(theta, false); + assert!(matches!( + tab.project(0, true, 1e-6), + Err(ProjectError::ImpossibleOutcome { .. }) + )); + let p = tab.project(0, true, TOL).unwrap(); + assert!((p - 1e-8).abs() < 1e-14); +} + +#[test] +fn lost_qubit_is_rejected() { + let mut tab = Tab::new_with_seed(1, 1e-12, 0); + tab.loss_channel(0, 1.0); + assert_eq!( + tab.project(0, false, TOL), + Err(ProjectError::LostQubit { addr0: 0 }) + ); +} + +#[test] +fn project_after_truncation_divides_by_norm() { + // Same setup as `test_measure_after_truncation_divides_by_norm`: rx drops + // a branch without renormalizing, but qubit 0 is still exactly |0⟩. + let theta = 2.0 * 0.25f64.asin(); + let mut tab = Tab::new_with_seed(1, 0.3, 0); + tab.rx(0, theta); + assert!((tab.project(0, false, TOL).unwrap() - 1.0).abs() < 1e-12); +} + +#[test] +fn chain_rule_matches_z_string_reference() { + let n = 3; + let mut base = Tab::new_with_seed(n, 1e-12, 0); + base.h(0); + base.t(0); + base.cnot(0, 1); + base.ry(1, 0.9); + base.rx(2, 1.3); + base.cz(1, 2); + base.h(2); + base.t(2); + base.cnot(2, 0); + + let mut total = 0.0; + for z in 0..1usize << n { + // P(z) = Π_q P(z_q | z_> q) & 1 == 1, TOL).unwrap_or(0.0)) + .product(); + + // Reference: ⟨ψ|z⟩⟨z|ψ⟩ = ⟨Π_q (I + (-1)^{z_q} Z_q)/2⟩, expanded into + // Z-string expectations. + let reference: f64 = (0..1usize << n) + .map(|subset| { + let word: String = (0..n) + .map(|q| if (subset >> q) & 1 == 1 { 'Z' } else { 'I' }) + .collect(); + let word: PauliWord< as Config>::Storage> = word.into(); + let sign = if (z & subset).count_ones() % 2 == 0 { + 1.0 + } else { + -1.0 + }; + sign * base.expectation(&word) + }) + .sum::() + / (1usize << n) as f64; + + assert!( + (p - reference).abs() < 1e-10, + "z={z:03b}: {p} vs {reference}" + ); + total += p; + } + assert!((total - 1.0).abs() < 1e-10, "Σ P(z) = {total}"); +} diff --git a/ppvm-python/src/ppvm/_core.pyi b/ppvm-python/src/ppvm/_core.pyi index d6da337f2..9c024c88d 100644 --- a/ppvm-python/src/ppvm/_core.pyi +++ b/ppvm-python/src/ppvm/_core.pyi @@ -127,6 +127,7 @@ class _GeneralizedTableauBase: def __deepcopy__(self, memo: object) -> _GeneralizedTableauBase: ... def measure(self, addr0: int) -> int: ... def measure_many(self, targets: Sequence[int]) -> list[int]: ... + def project(self, addr0: int, outcome: bool, tol: float = 1e-12) -> float: ... def current_measurement_record(self) -> list[int]: ... def coefficients(self) -> dict[int, complex]: ... def num_coefficients(self) -> int: ... diff --git a/ppvm-python/src/ppvm/generalized_tableau.py b/ppvm-python/src/ppvm/generalized_tableau.py index 528e5139d..52935e7a2 100644 --- a/ppvm-python/src/ppvm/generalized_tableau.py +++ b/ppvm-python/src/ppvm/generalized_tableau.py @@ -2,6 +2,7 @@ # SPDX-License-Identifier: Apache-2.0 import enum +import math from collections.abc import Iterable from dataclasses import InitVar, dataclass, field from typing import Literal, overload @@ -202,6 +203,38 @@ def measure_many(self, *targets: int | Iterable[int]) -> list[MeasurementResult] """ return [_BY_VALUE[v] for v in self._interface.measure_many(_normalize_targets(targets))] + def project(self, addr0: int, outcome: MeasurementResult, *, tol: float = 1e-12) -> float: + """Post-select a Z measurement of the specified qubit onto ``outcome``. + + Collapses the state as `measure` would have on ``outcome`` and records + it. Chaining it gives exact bitstring probabilities: + + ```python + p = 1.0 + for q, bit in enumerate(bits): + p *= tab.project(q, MeasurementResult(bit)) + ``` + + Args: + addr0: The index of the target qubit. + outcome: ``ZERO`` or ``ONE``. + tol: Smallest probability accepted for ``outcome``; must be + positive and finite. + + Returns: + The probability of measuring ``outcome`` before the projection. + + Raises: + ValueError: If ``outcome`` is ``LOST`` or ``tol`` is invalid, or, + leaving the state unchanged, if the qubit is lost or the + probability of ``outcome`` is below ``tol``. + """ + if outcome not in (MeasurementResult.ZERO, MeasurementResult.ONE): + raise ValueError(f"can only project onto ZERO or ONE, got {outcome!r}") + if not (tol > 0 and math.isfinite(tol)): + raise ValueError(f"tol must be positive and finite, got {tol}") + return self._interface.project(addr0, bool(outcome), tol) + def current_measurement_record(self) -> list[MeasurementResult]: """Return all measurement outcomes recorded so far. 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..59b5c923b --- /dev/null +++ b/ppvm-python/test/generalized_tableau/test_project.py @@ -0,0 +1,46 @@ +import math + +import pytest + +from ppvm import GeneralizedTableau +from ppvm.generalized_tableau import MeasurementResult + +ZERO, ONE, LOST = MeasurementResult.ZERO, MeasurementResult.ONE, MeasurementResult.LOST + + +def test_project_returns_probability_and_records(): + theta = 1.1 + tab = GeneralizedTableau(1) + tab.ry(0, theta=theta) + assert tab.project(0, ONE) == pytest.approx(math.sin(theta / 2) ** 2, abs=1e-10) + assert tab.expectation("Z") == pytest.approx(-1.0, abs=1e-10) + assert tab.current_measurement_record() == [ONE] + + +def test_impossible_outcome_raises(): + tab = GeneralizedTableau(1) + with pytest.raises(ValueError, match="probability"): + tab.project(0, ONE) + + +def test_tol_keyword(): + tab = GeneralizedTableau(1) + tab.ry(0, theta=2 * math.asin(1e-4)) # P(1) = 1e-8 + with pytest.raises(ValueError, match="probability"): + tab.project(0, ONE, tol=1e-6) + assert tab.project(0, ONE) == pytest.approx(1e-8, rel=1e-6) + + +@pytest.mark.parametrize("tol", [0.0, -1.0, math.inf, math.nan]) +def test_invalid_tol_raises(tol): + with pytest.raises(ValueError, match="tol"): + GeneralizedTableau(1).project(0, ZERO, tol=tol) + + +def test_lost_raises(): + tab = GeneralizedTableau(1) + with pytest.raises(ValueError, match="ZERO or ONE"): + tab.project(0, LOST) + tab.loss_channel(0, 1.0) + with pytest.raises(ValueError, match="lost"): + tab.project(0, ZERO) diff --git a/skills/ppvm-usage/SKILL.md b/skills/ppvm-usage/SKILL.md index 5e0913b71..1f0cc297c 100644 --- a/skills/ppvm-usage/SKILL.md +++ b/skills/ppvm-usage/SKILL.md @@ -150,6 +150,8 @@ shots = sample_stim(prog, shots=1000, n_qubits=5) `MeasurementResult` is an `IntEnum` (`ZERO`, `ONE`, `LOST`). Loss is first-class — a lost qubit measures `LOST`. +`tab.project(q, MeasurementResult.ONE)` post-selects instead of sampling: it collapses onto that outcome, records it, and returns its probability, so chaining it over all qubits gives the exact bitstring probability. It raises `ValueError` if the outcome has probability below `tol` (default `1e-12`) or the qubit is lost. + Stim leakage is `I_ERROR[leakage](p0, p1)` — same `p0`/`p1` as `leakage_channel`. A leaked qubit measures the pinned `ZERO`/`ONE`, not `LOST`. `I_ERROR[reset_loss] q…` / `I_ERROR[reset_leakage] q…` / `I_ERROR[reset_state] q…` (no args) map to `reset_loss_channel` / `reset_leakage_channel` / `reset_state`. `R`/`RX`/`RY` and the reset half of `MR`/`MRX`/`MRY` are full resets: lost and leaked qubits come back live (`MR` reports `LOST` / the pinned bit first). ## Rust API @@ -269,6 +271,7 @@ 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, outcome)` → probability of `outcome` (post-selection) | ✓ | — | ✓ | | `reset(q)` (also un-loses / un-leaks) | ✓ | — | ✓ | | `reset_state(q)` (live qubits only; lost / leaked untouched) | ✓ | — | ✓ | | `depolarize1(q, p=...)` | ✓ | ✓ | ✓ |