Add measure_batch and adopt it for STIM measurement

crates/ppvm-stim/src/executor.rs

@@ -10,6 +10,7 @@ use num::{Complex, One, ToPrimitive, Zero};
 use std::fmt::Debug;
 
 use ppvm_pauli_sum::prelude::*;
+use ppvm_tableau::measure_all::LossyMeasureAll;
 use ppvm_tableau::prelude::*;
 use stim_parser::ast::{GateName, MeasureName, NoiseName};
 use stim_parser::extended::{Axis, ExtendedInstruction, ExtendedProgram, RawPassthrough};
@@ -592,8 +593,17 @@ pub fn execute_validated<T, I, C>(
                 let noise = args.first().copied().unwrap_or(0.0);
                 match name {
                     MeasureName::M | MeasureName::MZ => {
-                        for &q in targets {
-                            results.push(tab.measure_noisy(q, noise));
+                        // Noiseless readout: batch the measurement so the shared
+                        // scratch is reused across targets, mirroring the batched
+                        // gate paths above. With readout noise, keep the per-qubit
+                        // loop so each measure/flip pair's RNG draws stay
+                        // interleaved exactly as before.
+                        if noise > 0.0 {
+                            for &q in targets {
+                                results.push(tab.measure_noisy(q, noise));
+                            }
+                        } else {
+                            results.extend(tab.measure_batch(targets));
                         }
                     }
                     // MR cannot delegate to `measure_noisy` because the reset must use

crates/ppvm-stim/tests/executor.rs

@@ -14,6 +14,12 @@ fn run(src: &str, n_qubits: usize) -> (Vec<Option<bool>>, Tab) {
     (results, tab)
 }
 
+fn run_seeded(src: &str, n_qubits: usize, seed: u64) -> Vec<Option<bool>> {
+    let prog = parse_extended(src).expect("parse_extended");
+    let mut tab: Tab = GeneralizedTableau::new_with_seed(n_qubits, 1e-10, seed);
+    execute(&prog, &mut tab).expect("execute")
+}
+
 #[test]
 fn x_then_measure_returns_one() {
     let (results, _) = run("X 0\nM 0", 1);
@@ -97,6 +103,22 @@ fn comments_and_annotations_are_no_ops() {
     assert_eq!(results, vec![Some(true)]);
 }
 
+#[test]
+fn grouped_measurement_matches_separate_measurements() {
+    // A single multi-target `M` must produce identical results to the same
+    // qubits measured by separate single-target `M` instructions, for every
+    // seed. Grouping only changes internal batching (shared measurement
+    // scratch); it must never change outcomes or the RNG-draw order.
+    let prelude = "H 0\nCX 0 1\nH 2\nCX 2 3\nS 1\nCZ 1 2\n";
+    let grouped = format!("{prelude}M 0 1 2 3");
+    let separate = format!("{prelude}M 0\nM 1\nM 2\nM 3");
+    for seed in 0..50 {
+        let g = run_seeded(&grouped, 4, seed);
+        let s = run_seeded(&separate, 4, seed);
+        assert_eq!(g, s, "grouped vs separate mismatch at seed={seed}");
+    }
+}
+
 #[test]
 fn measurement_buffer_is_pre_sized() {
     let prog = parse_extended("X 0\nM 0 1 2 3 4").unwrap();

crates/ppvm-tableau/benches/measure-all.rs

@@ -181,6 +181,16 @@ pub fn benchmark_suite_msd(c: &mut Criterion, name: impl AsRef<str>) {
         );
     });
 
+    // measure_batch over every qubit must match measure_all (no regression).
+    let all_indices: Vec<usize> = (0..tab.n_qubits()).collect();
+    group.bench_function("measure-batch-all", |b| {
+        b.iter_batched_ref(
+            || tab.clone(),
+            |tab| tab.measure_batch(&all_indices),
+            criterion::BatchSize::SmallInput,
+        );
+    });
+
     group.finish();
 }
 

crates/ppvm-tableau/benches/measure-scaling.rs

@@ -5,6 +5,7 @@ use std::time::Duration;
 
 use criterion::{Criterion, criterion_group, criterion_main};
 use ppvm_pauli_sum::config::fx64hash::Byte8F64;
+use ppvm_tableau::measure_all::LossyMeasureAll;
 use ppvm_tableau::prelude::*;
 
 type Tab = GeneralizedTableau<Byte8F64<2>, u128>;
@@ -25,6 +26,7 @@ fn bench_measure_scaling(c: &mut Criterion) {
     let n_qubits = 85;
     let mut group = c.benchmark_group("measure-scaling");
 
+    let all: Vec<usize> = (0..n_qubits).collect();
     for n_t in [8, 10, 12, 14] {
         let tab = build_circuit(n_qubits, n_t);
         group.bench_function(format!("t{n_t}-{n_qubits}q"), |b| {
@@ -38,6 +40,16 @@ fn bench_measure_scaling(c: &mut Criterion) {
                 criterion::BatchSize::SmallInput,
             );
         });
+        // Batching every qubit should be no slower than the per-qubit loop.
+        group.bench_function(format!("batch-t{n_t}-{n_qubits}q"), |b| {
+            b.iter_batched_ref(
+                || tab.fork(Some(42)),
+                |t| {
+                    let _ = t.measure_batch(&all);
+                },
+                criterion::BatchSize::SmallInput,
+            );
+        });
     }
 
     group.finish();

crates/ppvm-tableau/benches/measure_batch_profiling.md

@@ -0,0 +1,70 @@
+<!--
+SPDX-FileCopyrightText: 2026 The PPVM Authors
+SPDX-License-Identifier: Apache-2.0
+-->
+
+# `measure_batch` profiling: where the time goes (2026-06-22)
+
+Question: beyond reusing a common `MeasureScratch`, is there an algorithmic win
+for `measure_batch`? This note records the profiling that answers it. **No
+optimization has been applied yet** — these are findings to act on later.
+
+## Harness
+
+- `examples/profile_measure_batch.rs` — builds a full-width (85-qubit) entangled
+  Clifford state so every measurement takes the case-a (random) path, with `n_t`
+  T-gates controlling the coefficient count. Three workloads bracket the regimes:
+  `few` (~1 coefficient), `mid` (~128), `large` (~1024). `quick` mode prints
+  timings + the achieved coefficient count; `flame` mode runs a sustained loop.
+- `scripts/profile_measure_batch.sh` — records all three with **samply** (no
+  sudo on macOS, unlike dtrace/flamegraph) into `target/profiles/`, then prints
+  a top-functions summary.
+- `scripts/samply_top.py` — turns a samply `profile.json.gz` (+ presymbolicated
+  `.syms.json` sidecar) into a top-functions-by-self/inclusive-time table,
+  headless (no browser UI needed).
+
+Reproduce: `./scripts/profile_measure_batch.sh`
+
+## Results (self time; per-shot from `quick` mode)
+
+Per-shot `measure_batch` cost: few 17.7 µs, mid 21 µs, large 48 µs (85 qubits).
+
+| function (self)                          | few (1) | mid (128) | large (1024) |
+|------------------------------------------|---------|-----------|--------------|
+| `compute_decomposition` (tableau scan)   | 39%     | 33%       | 14%          |
+| `measure_one_with_scratch` *(see note)*  | 56%     | 54%       | 49%          |
+| `insert` (build case-a HashMap)          | ~1%     | 7%        | 23%          |
+| `retain` (drain Vec↔HashMap)             | ~1%     | 3%        | 9%           |
+| memset / memmove                         | 1.4%    | 1.3%      | 4%           |
+
+*Note: `measure_with_scratch`, `project_case_a`, and
+`update_tableau_according_to_outcome` are inlined into `measure_one_with_scratch`.*
+
+## Conclusions
+
+The cost is two distinct regimes:
+
+1. **Tableau-bound floor (dominates few/mid).** `compute_decomposition` costs a
+   *constant* ~6.8 µs regardless of coefficient count (39% at 1 coeff → 14% at
+   1024, same absolute µs). With the inlined `update_tableau_according_to_outcome`,
+   the tableau work is a fixed **~16 µs O(n²)-per-batch floor** (each of the n
+   measured qubits scans all 2n stabilizer/destabilizer rows). This is ~95% of
+   the cost on Clifford-heavy / low-coefficient states.
+
+2. **Coefficient-bound `Vec`↔`HashMap` round-trip (dominates large).** Each
+   case-a measurement drains `self.coefficients` (`Vec`) into the scratch
+   HashMap, works, then drains back; the next qubit rebuilds it. At 1024
+   coefficients `insert` + `retain` + their memset/memmove churn is **~35%
+   (~17 µs)** of the cost, and ~0% at low coefficient counts. The shared scratch
+   keeps the *allocation*, but not the *representation*, across qubits.
+
+## Recommended next steps (not yet implemented)
+
+1. **Keep case-a coefficients HashMap-resident across the batch** — convert
+   `Vec`→`HashMap` once and back once *per batch* instead of per qubit.
+   Estimated ~30–35% on coefficient-heavy states; well-scoped refactor of
+   `measure_with_scratch` / `project_case_a` that leaves the RNG draw order
+   unchanged (so the existing characterization tests still apply).
+2. **Batch the O(n²) tableau floor** — a blocked/incremental measurement that
+   avoids rescanning all 2n rows per qubit. Larger, riskier algorithmic change;
+   payoff concentrated on cheap (low-coefficient) states. Treat separately.

crates/ppvm-tableau/examples/profile_measure_batch.rs

@@ -0,0 +1,139 @@
+// SPDX-FileCopyrightText: 2026 The PPVM Authors
+// SPDX-License-Identifier: Apache-2.0
+
+//! Profiling driver for `measure_batch` across coefficient-count regimes.
+//!
+//! `measure_batch` has two cost components per measured qubit:
+//!   - tableau work  — `compute_decomposition` (scans all stabilizers +
+//!     destabilizers) and, on case-a, `update_tableau_according_to_outcome`.
+//!     This is O(n_qubits) per qubit and *independent of the coefficient count*.
+//!   - coefficient work — the case-a HashMap build / overlap / merge / normalize
+//!     passes, which scale with the number of stored coefficients.
+//!
+//! The three workloads bracket the crossover between the two regimes. Every
+//! state is a full-width entangled Clifford state (so measurements take the
+//! case-a path), with `n_t` T-gates controlling the coefficient count:
+//!   - `few`   — ~1 coefficient    (tableau-bound)
+//!   - `mid`   — ~100 coefficients
+//!   - `large` — ~1000 coefficients (coefficient-bound)
+//!
+//! Usage:
+//!   # quick timing + achieved coefficient count (use this to calibrate n_t)
+//!   cargo run -p ppvm-tableau --example profile_measure_batch --release -- mid
+//!
+//!   # sustained hot loop (~FLAME_SECS seconds, default 6) for a sampling
+//!   # profiler. scripts/profile_measure_batch.sh drives this with samply
+//!   # (no sudo on macOS) and summarizes via scripts/samply_top.py:
+//!   samply record --save-only --unstable-presymbolicate -o p.json.gz -- \
+//!     target/release/examples/profile_measure_batch mid flame
+
+use ppvm_pauli_sum::config::fx64hash::Byte8F64;
+use ppvm_tableau::measure_all::LossyMeasureAll;
+use ppvm_tableau::prelude::*;
+use std::time::{Duration, Instant};
+
+type Tab = GeneralizedTableau<Byte8F64<2>, u128>;
+
+const N_QUBITS: usize = 85;
+
+/// Number of T-gates per workload, tuned so the post-circuit coefficient count
+/// lands near the target. Recalibrate with the quick-timing mode if the state
+/// builder changes.
+fn n_t_for(workload: &str) -> usize {
+    match workload {
+        "few" => 0,
+        "mid" => 7,
+        "large" => 10,
+        other => panic!("unknown workload {other:?}; use few|mid|large"),
+    }
+}
+
+/// Full-width entangled state: H on every qubit (Clifford superposition, so
+/// measurements hit the case-a path), `n_t` T-gates to branch the coefficient
+/// vector, then a CZ chain to entangle.
+fn build_state(n_t: usize) -> Tab {
+    let mut tab: Tab = GeneralizedTableau::new(N_QUBITS, 1e-10);
+    for q in 0..N_QUBITS {
+        tab.h(q);
+    }
+    for q in 0..n_t {
+        tab.t(q);
+    }
+    for q in 0..N_QUBITS - 1 {
+        tab.cz(q, q + 1);
+    }
+    tab
+}
+
+fn median(mut xs: Vec<Duration>) -> Duration {
+    xs.sort();
+    xs[xs.len() / 2]
+}
+
+fn main() {
+    let mut args = std::env::args().skip(1);
+    let workload = args.next().unwrap_or_else(|| "mid".to_string());
+    let mode = args.next().unwrap_or_else(|| "quick".to_string());
+
+    let n_t = n_t_for(&workload);
+    let base = build_state(n_t);
+    let n = base.n_qubits();
+    let n_coeffs = base.coefficients.len();
+    let all: Vec<usize> = (0..n).collect();
+
+    // Setup info goes to stderr so it never pollutes a flamegraph's sampled window.
+    eprintln!("workload={workload}  n_qubits={n}  n_t={n_t}  coefficients={n_coeffs}  mode={mode}");
+
+    match mode.as_str() {
+        "quick" => {
+            let n_runs = 200;
+            let mut shot = Vec::with_capacity(n_runs);
+            for _ in 0..n_runs {
+                let start = Instant::now();
+                let mut t = base.fork(Some(42));
+                std::hint::black_box(t.measure_batch(&all));
+                shot.push(start.elapsed());
+            }
+            let mut forks = Vec::with_capacity(n_runs);
+            for _ in 0..n_runs {
+                let start = Instant::now();
+                let t = base.fork(Some(42));
+                forks.push(start.elapsed());
+                std::hint::black_box(t);
+            }
+            let shot_m = median(shot);
+            let fork_m = median(forks);
+            eprintln!("per-shot (fork + measure_batch) median over {n_runs}: {shot_m:?}");
+            eprintln!("    fork only:     {fork_m:?}");
+            eprintln!("    measure_batch: {:?}", shot_m.saturating_sub(fork_m));
+            eprintln!(
+                "    per qubit:     {:?}",
+                shot_m.saturating_sub(fork_m) / n as u32
+            );
+        }
+        "flame" => {
+            // Size the hot loop to ~10s of sampling regardless of workload cost.
+            let mut probe = Vec::with_capacity(50);
+            for _ in 0..50 {
+                let start = Instant::now();
+                let mut t = base.fork(Some(42));
+                std::hint::black_box(t.measure_batch(&all));
+                probe.push(start.elapsed());
+            }
+            let per = median(probe).max(Duration::from_nanos(1));
+            let secs: u64 = std::env::var("FLAME_SECS")
+                .ok()
+                .and_then(|s| s.parse().ok())
+                .unwrap_or(6);
+            let iters = (Duration::from_secs(secs).as_nanos() / per.as_nanos()).max(2_000) as u64;
+            eprintln!("flame: ~{per:?}/shot -> {iters} iterations (~{secs}s)");
+            // fork feeds measure_batch and the result is black-boxed, so neither
+            // call can be hoisted out of the loop.
+            for _ in 0..iters {
+                let mut t = std::hint::black_box(base.fork(Some(42)));
+                std::hint::black_box(t.measure_batch(&all));
+            }
+        }
+        other => panic!("unknown mode {other:?}; use quick|flame"),
+    }
+}