From fdaa6a7a47c9cae6600454d9f0d3387d44aed728 Mon Sep 17 00:00:00 2001
From: Ciaran Ryan-Anderson <ciaranra@gmail.com>
Date: Thu, 6 Aug 2026 21:38:39 -0600
Subject: [PATCH 1/2] Consolidate DEM batch sampling onto SampleBatch with bulk
 accessors (#448)

---
 docs/user-guide/dem-from-guppy.md             |  46 +-
 examples/surface/brickwork_sweep.py           |   4 +-
 examples/surface/decoder_comparison.py        |   2 +-
 examples/surface/dem_comparison.py            |   2 +-
 examples/surface/dem_method_ler_comparison.py |   2 +-
 examples/surface/dem_tutorial.py              |   2 +-
 examples/surface/dem_vs_stabilizer.py         |   2 +-
 examples/surface/eeg_vs_statevec.py           |   4 +-
 examples/surface/generate_data.py             |   2 +-
 examples/surface/inner_decoder_study.py       |   2 +-
 examples/surface/ml_lookup_decoder.py         |   4 +-
 examples/surface/validate_dem_correlations.py |   2 +-
 examples/surface/validate_dem_generators.py   |   4 +-
 python/pecos-rslib/pecos_rslib.pyi            | 489 ++++++++++++++++++
 .../src/fault_tolerance_bindings.rs           | 248 ++++++---
 .../src/pecos/qec/surface/decode.py           |  13 +-
 .../tests/qec/surface/test_circuit_fuzz.py    |  16 +-
 ...test_logical_subgraph_region_comparison.py |  16 +-
 .../qec/surface/test_pauli_twirl_handoff.py   |  53 ++
 .../tests/qec/test_dem_equivalence.py         |  30 +-
 .../tests/qec/test_dem_sampler.py             |  12 +-
 .../tests/qec/test_dem_sampler_modes.py       |  23 +-
 .../tests/qec/test_dem_sampler_vs_stim.py     |   2 +-
 .../tests/qec/test_from_guppy_dem.py          |   2 +-
 .../tests/qec/test_parsed_dem_sampler.py      |   8 +-
 .../tests/qec/test_sample_batch.py            |  52 +-
 .../qec/test_traced_qis_slow_integration.py   |   2 +-
 .../tests/qec/test_wide_observables.py        |  16 +-
 scripts/bench_raw_meas_sampling.py            |   8 +-
 scripts/compare_meas_sampling_pipeline.py     |   4 +-
 30 files changed, 926 insertions(+), 146 deletions(-)

diff --git a/docs/user-guide/dem-from-guppy.md b/docs/user-guide/dem-from-guppy.md
index 1430202fe..f04552091 100644
--- a/docs/user-guide/dem-from-guppy.md
+++ b/docs/user-guide/dem-from-guppy.md
@@ -161,7 +161,7 @@ assert dem.num_observables == 1
 
 # Sample syndromes/observables and decode them, all PECOS-native.
 sampler = dem.to_sampler()
-batch = sampler.generate_samples(1000, 0)
+batch = sampler.sample_batch(1000, 0)
 assert batch.num_shots == 1000
 
 decoder = PyMatchingDecoder.from_dem(dem.to_string_decomposed())
@@ -177,6 +177,50 @@ The DEM built this way is identical to the reference DEM produced by the
 surface traced-QIS pipeline — the abstract builder's metadata and the
 traced Guppy program agree on measurement order.
 
+## Sampling and Comparing Decoders with `SampleBatch`
+
+`DemSampler.sample_batch` and `ParsedDem.sample_batch` return a
+`SampleBatch`. The detector events and observable flips remain in Rust memory,
+so the same shots can be passed to several decoders without copying them
+through Python. When Python data is needed, `detector_events()` and
+`observable_flips()` return shots-major `list[list[bool]]` values.
+
+<!--test-name: dem_sample_batch_workflow-->
+```python
+from pecos_rslib.qec import ParsedDem, SampleBatch
+
+dem_text = "error(0.1) D0 L0"
+sampler = ParsedDem.from_string(dem_text).to_dem_sampler()
+batch = sampler.sample_batch(32, seed=42)
+
+assert isinstance(batch, SampleBatch)
+assert batch.num_shots == 32
+
+# Materialize shots-major Python data only when it is needed.
+detector_events = batch.detector_events()
+observable_flips = batch.observable_flips()
+assert len(detector_events) == len(observable_flips) == 32
+assert all(len(shot) == 1 for shot in detector_events)
+assert all(len(shot) == 1 for shot in observable_flips)
+
+# Aggregate errors, inspect individual predictions, or collect timings.
+error_count = batch.decode_count(dem_text, "pymatching")
+predictions = batch.decode_each(dem_text, "pymatching")
+stats = batch.decode_stats(dem_text, "pymatching")
+
+assert 0 <= error_count <= batch.num_shots
+assert len(predictions) == batch.num_shots
+assert stats.num_shots == batch.num_shots
+```
+
+Use `decode_count` for a logical-error total, `decode_each` to inspect the
+prediction for every shot, and `decode_stats` for error counts plus per-shot
+timing statistics. The parallel `decode_count_parallel` and
+`decode_stats_parallel` variants distribute slow decoder work across multiple
+workers. A former raw-list call such as
+`detectors, observables = sampler.sample_batch(...)` becomes a batch call
+followed by the two bulk accessors shown above.
+
 ## Choosing the Selene Runtime
 
 `from_guppy(..., runtime=...)` forwards to `pecos.selene_engine(runtime)`,
diff --git a/examples/surface/brickwork_sweep.py b/examples/surface/brickwork_sweep.py
index 55d501e3a..3782d50ba 100644
--- a/examples/surface/brickwork_sweep.py
+++ b/examples/surface/brickwork_sweep.py
@@ -196,7 +196,7 @@ def run_sweep(
                     t0 = time.perf_counter()
                     parsed = ParsedDem.from_string(dem_str)
                     rust_sampler = parsed.to_dem_sampler()
-                    batch = rust_sampler.generate_samples(shots, seed=circuit_seed + cell_idx)
+                    batch = rust_sampler.sample_batch(shots, seed=circuit_seed + cell_idx)
                     sample_sec = time.perf_counter() - t0
 
                     point = BrickworkPoint(
@@ -690,7 +690,7 @@ def main():
                 sc = b.stab_coords()
                 dem_str = b.build_dem(p1=p, p2=p, p_meas=p, p_prep=p)
                 parsed = ParsedDem.from_string(dem_str)
-                batch = parsed.to_dem_sampler().generate_samples(args.shots, seed=args.seed)
+                batch = parsed.to_dem_sampler().sample_batch(args.shots, seed=args.seed)
 
                 point = BrickworkPoint(
                     distance=d,
diff --git a/examples/surface/decoder_comparison.py b/examples/surface/decoder_comparison.py
index 926b5de7b..a53efc71f 100644
--- a/examples/surface/decoder_comparison.py
+++ b/examples/surface/decoder_comparison.py
@@ -201,7 +201,7 @@ def run_comparison(
 
             # Generate samples once
             t0 = time.perf_counter()
-            sample_batch = sampler.sampler.generate_samples(shots, seed=seed + config_idx)
+            sample_batch = sampler.sampler.sample_batch(shots, seed=seed + config_idx)
             sample_seconds = time.perf_counter() - t0
 
             results: list[DecoderResult] = []
diff --git a/examples/surface/dem_comparison.py b/examples/surface/dem_comparison.py
index d9b178981..1d103118f 100644
--- a/examples/surface/dem_comparison.py
+++ b/examples/surface/dem_comparison.py
@@ -92,7 +92,7 @@ def extract_det_rates(results):
     # 2. DemSampler.from_circuit
     t0 = time.perf_counter()
     sampler_fc = DemSampler.from_circuit(dag, p1=p, p2=p, p_meas=p, p_prep=p)
-    batch_fc = sampler_fc.generate_samples(num_shots=shots, seed=seed)
+    batch_fc = sampler_fc.sample_batch(num_shots=shots, seed=seed)
     dem_fc = [0.0] * num_dets
     for i in range(shots):
         syn = batch_fc.get_syndrome(i)
diff --git a/examples/surface/dem_method_ler_comparison.py b/examples/surface/dem_method_ler_comparison.py
index 84dd7ef28..e313764fd 100644
--- a/examples/surface/dem_method_ler_comparison.py
+++ b/examples/surface/dem_method_ler_comparison.py
@@ -313,7 +313,7 @@ def run_comparison(
                     **sampler_params,
                     idle_rz=idle_rz if idle_rz > 0 else None,
                 )
-                batch = sampler.generate_samples(shots, seed=seed)
+                batch = sampler.sample_batch(shots, seed=seed)
             t_sample = time.perf_counter() - t0
             print(f"  Sampled {shots} shots in {t_sample:.2f}s")
 
diff --git a/examples/surface/dem_tutorial.py b/examples/surface/dem_tutorial.py
index 545d66f18..92b91b66c 100644
--- a/examples/surface/dem_tutorial.py
+++ b/examples/surface/dem_tutorial.py
@@ -81,7 +81,7 @@ def main():
     # ================================================================
     shots = 100_000
     sampler = DemSampler.from_circuit(tc, p1=p, p2=p, p_meas=p, p_prep=p)
-    batch = sampler.generate_samples(num_shots=shots, seed=42)
+    batch = sampler.sample_batch(num_shots=shots, seed=42)
 
     # Compute per-detector firing rates from DEM sampling
     num_dets = len(dets)
diff --git a/examples/surface/dem_vs_stabilizer.py b/examples/surface/dem_vs_stabilizer.py
index 5a48f938a..b791a5579 100644
--- a/examples/surface/dem_vs_stabilizer.py
+++ b/examples/surface/dem_vs_stabilizer.py
@@ -69,7 +69,7 @@ def run_comparison(*, distance, rounds, basis, p, shots, seed):
     # 2. DemSampler.from_circuit
     t0 = time.perf_counter()
     sampler_fc = DemSampler.from_circuit(dag, p1=p, p2=p, p_meas=p, p_prep=p)
-    batch_fc = sampler_fc.generate_samples(num_shots=shots, seed=seed)
+    batch_fc = sampler_fc.sample_batch(num_shots=shots, seed=seed)
     dem_fc = [0.0] * num_dets
     for i in range(shots):
         syn = batch_fc.get_syndrome(i)
diff --git a/examples/surface/eeg_vs_statevec.py b/examples/surface/eeg_vs_statevec.py
index 82fd38915..d3e07322f 100644
--- a/examples/surface/eeg_vs_statevec.py
+++ b/examples/surface/eeg_vs_statevec.py
@@ -93,14 +93,14 @@ def run_comparison(
             t0 = time.perf_counter()
             dem_taylor_str = perturbative_dem(tc, idle_rz=theta)
             sampler_taylor = DemSampler.from_dem_string(dem_taylor_str)
-            batch_taylor = sampler_taylor.generate_samples(num_shots=shots, seed=seed)
+            batch_taylor = sampler_taylor.sample_batch(num_shots=shots, seed=seed)
             taylor_sample_time = time.perf_counter() - t0
 
             # Heisenberg DEM → sampler
             t0 = time.perf_counter()
             dem_heis_str = coherent_dem_exact(tc, idle_rz=theta)
             sampler_heis = DemSampler.from_dem_string(dem_heis_str)
-            batch_heis = sampler_heis.generate_samples(num_shots=shots, seed=seed)
+            batch_heis = sampler_heis.sample_batch(num_shots=shots, seed=seed)
             heis_sample_time = time.perf_counter() - t0
 
             # Compute per-detector rates from DEM samples
diff --git a/examples/surface/generate_data.py b/examples/surface/generate_data.py
index 1eb99df7b..9a2ca6dc7 100644
--- a/examples/surface/generate_data.py
+++ b/examples/surface/generate_data.py
@@ -221,7 +221,7 @@ def generate(
 
                 # Sample once
                 t0 = time.perf_counter()
-                batch = sampler.sampler.generate_samples(shots, seed=seed + cell_idx)
+                batch = sampler.sampler.sample_batch(shots, seed=seed + cell_idx)
                 sample_seconds = time.perf_counter() - t0
 
                 point = DataPoint(
diff --git a/examples/surface/inner_decoder_study.py b/examples/surface/inner_decoder_study.py
index 9e6ba9a84..9c04acbbd 100644
--- a/examples/surface/inner_decoder_study.py
+++ b/examples/surface/inner_decoder_study.py
@@ -164,7 +164,7 @@ def measure_cell(
     else:
         dem = builder.build_dem(p1=p, p2=p, p_meas=p)
     sc = builder.stab_coords()
-    batch = ParsedDem.from_string(dem).to_dem_sampler().generate_samples(n, seed=seed)
+    batch = ParsedDem.from_string(dem).to_dem_sampler().sample_batch(n, seed=seed)
 
     cells: list[Cell] = []
     for inner in inners:
diff --git a/examples/surface/ml_lookup_decoder.py b/examples/surface/ml_lookup_decoder.py
index 6fefa9ca2..8329a0d1c 100644
--- a/examples/surface/ml_lookup_decoder.py
+++ b/examples/surface/ml_lookup_decoder.py
@@ -156,7 +156,7 @@ def main():
 
         sampler_params = {k: v for k, v in noise_params.items() if k in ("p1", "p2", "p_meas", "p_prep")}
         sampler = DemSampler.from_circuit(tc, **sampler_params)
-        train_batch = sampler.generate_samples(args.shots, seed=args.seed)
+        train_batch = sampler.sample_batch(args.shots, seed=args.seed)
 
     t_sample = time.perf_counter() - t0
     print(f"  Sampled in {t_sample:.2f}s")
@@ -198,7 +198,7 @@ def main():
             observable_masks2.append(obs_mask)
         test_batch = SampleBatch(detection_events2, observable_masks2)
     else:
-        test_batch = sampler.generate_samples(test_shots, seed=args.seed + 1000)
+        test_batch = sampler.sample_batch(test_shots, seed=args.seed + 1000)
 
     # Decode with lookup
     errors_lookup, n = decode_with_lookup(test_batch, table, num_dets)
diff --git a/examples/surface/validate_dem_correlations.py b/examples/surface/validate_dem_correlations.py
index 78f84ed9e..cf3abd026 100644
--- a/examples/surface/validate_dem_correlations.py
+++ b/examples/surface/validate_dem_correlations.py
@@ -88,7 +88,7 @@ def dem_detector_events(tc, noise_kw, shots, seed):
 
     full_kw = {k: noise_kw.get(k, 0.0) for k in ["p1", "p2", "p_meas", "p_prep"]}
     sampler = DemSampler.from_circuit(tc, **full_kw)
-    batch = sampler.generate_samples(num_shots=shots, seed=seed)
+    batch = sampler.sample_batch(num_shots=shots, seed=seed)
     num_dets = len(json.loads(tc.get_meta("detectors")))
 
     events = []
diff --git a/examples/surface/validate_dem_generators.py b/examples/surface/validate_dem_generators.py
index 66efd020a..af8fb6765 100644
--- a/examples/surface/validate_dem_generators.py
+++ b/examples/surface/validate_dem_generators.py
@@ -136,7 +136,7 @@ def dem_sampler_rates(tc, noise_kw, shots, seed, num_dets):
     from pecos_rslib.qec import DemSampler
 
     sampler = DemSampler.from_circuit(tc, **full_noise_kw(noise_kw))
-    batch = sampler.generate_samples(num_shots=shots, seed=seed)
+    batch = sampler.sample_batch(num_shots=shots, seed=seed)
     rates = [0.0] * num_dets
     for i in range(shots):
         syn = batch.get_syndrome(i)
@@ -162,7 +162,7 @@ def dem_builder_rates(tc, noise_kw, shots, seed, num_dets):
         .build()
     )
     sampler = dem.to_sampler()
-    batch = sampler.generate_samples(num_shots=shots, seed=seed)
+    batch = sampler.sample_batch(num_shots=shots, seed=seed)
     rates = [0.0] * num_dets
     for i in range(shots):
         syn = batch.get_syndrome(i)
diff --git a/python/pecos-rslib/pecos_rslib.pyi b/python/pecos-rslib/pecos_rslib.pyi
index f5a792e9f..619a48ae4 100644
--- a/python/pecos-rslib/pecos_rslib.pyi
+++ b/python/pecos-rslib/pecos_rslib.pyi
@@ -2030,6 +2030,495 @@ class WasmForeignObject:
         """Support for pickle deserialization."""
         ...
 
+# =============================================================================
+# Quantum Error Correction Types
+# =============================================================================
+
+class qec:
+    """Fault-tolerance and detector-error-model submodule."""
+
+    PAULI_I: int
+    PAULI_X: int
+    PAULI_Y: int
+    PAULI_Z: int
+
+    class FaultLocation:
+        @property
+        def node(self) -> int: ...
+        @property
+        def qubits(self) -> list[int]: ...
+        @property
+        def before(self) -> bool: ...
+        @property
+        def gate_type(self) -> str: ...
+
+    class DagFaultInfluenceMap:
+        @property
+        def num_locations(self) -> int: ...
+        @property
+        def num_detectors(self) -> int: ...
+        @property
+        def num_dem_outputs(self) -> int: ...
+        @property
+        def num_observables(self) -> int: ...
+        @property
+        def num_tracked_paulis(self) -> int: ...
+        def get_locations(self) -> list[qec.FaultLocation]: ...
+        def get_location(self, loc_idx: int) -> qec.FaultLocation | None: ...
+        def classify_fault(self, loc_idx: int, pauli: int) -> tuple[bool, bool]: ...
+        def get_detector_indices(self, loc_idx: int, pauli: int) -> list[int]: ...
+        def get_dem_output_indices(self, loc_idx: int, pauli: int) -> list[int]: ...
+        def get_internal_dem_output_indices(self, loc_idx: int, pauli: int) -> list[int]: ...
+        def get_tracked_pauli_indices(self, loc_idx: int, pauli: int) -> list[int]: ...
+        def get_observable_indices(self, loc_idx: int, pauli: int) -> list[int]: ...
+        def has_detector_flips(self, loc_idx: int, pauli: int) -> bool: ...
+        def has_dem_output_flips(self, loc_idx: int, pauli: int) -> bool: ...
+        def has_observable_flips(self, loc_idx: int, pauli: int) -> bool: ...
+        def has_tracked_pauli_flips(self, loc_idx: int, pauli: int) -> bool: ...
+        def merge_dem_outputs_from(self, other: qec.DagFaultInfluenceMap) -> None: ...
+        def memory_stats(self) -> dict[str, Any]: ...
+        def export_csr(self) -> dict[str, Any]: ...
+        def measurements(self) -> list[tuple[int, int, int]]: ...
+        def __len__(self) -> int: ...
+
+    class DagFaultAnalyzer:
+        def __init__(self, dag: DagCircuit) -> None: ...
+        def build_influence_map(self) -> qec.DagFaultInfluenceMap: ...
+        @property
+        def max_node(self) -> int: ...
+        @property
+        def max_qubit(self) -> int: ...
+
+    class InfluenceBuilder:
+        def __init__(self, dag: DagCircuit) -> None: ...
+        def with_tracked_x(self, qubits: Sequence[int]) -> qec.InfluenceBuilder: ...
+        def with_tracked_z(self, qubits: Sequence[int]) -> qec.InfluenceBuilder: ...
+        def with_tracked_pauli(self, entries: Sequence[tuple[int, str]]) -> qec.InfluenceBuilder: ...
+        def with_circuit_annotations(self) -> qec.InfluenceBuilder: ...
+        def build(self) -> qec.DagFaultInfluenceMap: ...
+
+    class PauliFrameLookup:
+        @staticmethod
+        def from_circuit(
+            dag: DagCircuit,
+            detectors: Sequence[Sequence[int]],
+            observables: Sequence[Sequence[int]],
+        ) -> qec.PauliFrameLookup: ...
+        @property
+        def num_pauli_sites(self) -> int: ...
+        @property
+        def num_tracked_paulis(self) -> int: ...
+        @property
+        def num_detectors(self) -> int: ...
+        @property
+        def num_observables(self) -> int: ...
+        def row(self, tracked_idx: int) -> tuple[list[int], list[int]]: ...
+        def mask_firings(self, pauli_masks: Any) -> list[list[bool]]: ...
+        def compute_mask_xor(self, pauli_masks: Any) -> tuple[list[list[bool]], list[list[bool]]]: ...
+
+    class DetectorErrorModel:
+        @staticmethod
+        def from_circuit(
+            circuit: DagCircuit | TickCircuit,
+            p1: float = ...,
+            p2: float = ...,
+            p_meas: float = ...,
+            p_prep: float = ...,
+            **noise_options: Any,
+        ) -> qec.DetectorErrorModel: ...
+        @staticmethod
+        def from_pecos_metadata_json(json: str) -> qec.DetectorErrorModel: ...
+        @property
+        def num_detectors(self) -> int: ...
+        @property
+        def num_dem_outputs(self) -> int: ...
+        @property
+        def num_observables(self) -> int: ...
+        @property
+        def num_tracked_paulis(self) -> int: ...
+        @property
+        def num_contributions(self) -> int: ...
+        def to_string(self) -> str: ...
+        def to_string_decomposed(self) -> str: ...
+        def to_string_decomposed_maximally(self) -> str: ...
+        def to_string_decomposed_with_two_detector_direct_policy(self, policy: str) -> str: ...
+        def to_string_source_decomposed(self) -> str: ...
+        def to_string_source_graphlike_decomposed(self) -> str: ...
+        def to_string_terminal_graphlike_decomposed(self) -> str: ...
+        def to_string_graphlike_search_decomposed(self) -> str: ...
+        def all_contribution_effects(self) -> list[Any]: ...
+        def contribution_effect_summaries(self) -> list[Any]: ...
+        def contribution_render_records(self) -> list[Any]: ...
+        def contribution_render_records_with_two_detector_direct_policy(self, policy: str) -> list[Any]: ...
+        def contribution_render_summaries(self) -> list[Any]: ...
+        def contribution_render_summaries_with_two_detector_direct_policy(self, policy: str) -> list[Any]: ...
+        def contribution_source_graphlike_render_records(self) -> list[Any]: ...
+        def contributions_for_effect(self, detectors: Sequence[int], dem_outputs: Sequence[int]) -> list[Any]: ...
+        def contributions_for_mechanism(self, detectors: Sequence[int]) -> list[Any]: ...
+        def to_sampler(self) -> qec.DemSampler: ...
+
+    class DemBuilder:
+        def __init__(self, influence_map: qec.DagFaultInfluenceMap) -> None: ...
+        def with_noise(
+            self,
+            p1: float,
+            p2: float,
+            p_meas: float,
+            p_prep: float,
+            **noise_options: Any,
+        ) -> qec.DemBuilder: ...
+        def with_detectors_json(self, json: str) -> qec.DemBuilder: ...
+        def with_observables_json(self, json: str) -> qec.DemBuilder: ...
+        def with_measurement_order(self, order: Sequence[int]) -> qec.DemBuilder: ...
+        def with_num_measurements(self, num: int) -> qec.DemBuilder: ...
+        def with_exact_branch_replay_circuit(self, circuit: DagCircuit) -> qec.DemBuilder: ...
+        def build(self) -> qec.DetectorErrorModel: ...
+        def build_with_source_tracking(self) -> qec.DetectorErrorModel: ...
+
+    class SampleBatch:
+        def __init__(
+            self,
+            detection_events: Sequence[Sequence[int | bool]],
+            observable_masks: Sequence[int],
+            *,
+            num_observables: int | None = ...,
+        ) -> None: ...
+        @property
+        def num_shots(self) -> int: ...
+        def get_syndrome(self, i: int) -> list[int]: ...
+        def get_observable_mask(self, i: int) -> int: ...
+        def get_observable_mask_wide(self, i: int) -> int: ...
+        def detector_events(self) -> list[list[bool]]: ...
+        def observable_flips(self) -> list[list[bool]]: ...
+        def decode_count(self, dem: str, decoder_type: str = ...) -> int: ...
+        def decode_each(self, dem: str, decoder_type: str = ...) -> list[int]: ...
+        def decode_count_parallel(
+            self,
+            dem: str,
+            decoder_type: str = ...,
+            num_workers: int | None = ...,
+        ) -> int: ...
+        def decode_count_batch(self, dem: str) -> int: ...
+        def decode_stats(self, dem: str, decoder_type: str = ...) -> qec.DecodeStats: ...
+        def decode_stats_parallel(
+            self,
+            dem: str,
+            decoder_type: str = ...,
+            num_workers: int | None = ...,
+        ) -> qec.DecodeStats: ...
+
+    class DecodeStats:
+        num_shots: int
+        num_errors: int
+        logical_error_rate: float
+        total_seconds: float
+        per_shot_mean: float
+        per_shot_median: float
+        per_shot_p99: float
+        per_shot_min: float
+        per_shot_max: float
+        quantiles: list[tuple[float, float]]
+
+    class DemSampler:
+        @staticmethod
+        def from_circuit(
+            circuit: DagCircuit | TickCircuit,
+            p1: float = ...,
+            p2: float = ...,
+            p_meas: float = ...,
+            p_prep: float = ...,
+            **noise_options: Any,
+        ) -> qec.DemSampler: ...
+        @staticmethod
+        def from_dem_string(dem_string: str) -> qec.DemSampler: ...
+        @staticmethod
+        def raw_uniform(influence_map: qec.DagFaultInfluenceMap, p_error: float) -> qec.DemSampler: ...
+        @staticmethod
+        def raw(
+            influence_map: qec.DagFaultInfluenceMap,
+            p1: float,
+            p2: float,
+            p_meas: float,
+            p_prep: float,
+        ) -> qec.DemSampler: ...
+        @staticmethod
+        def with_detectors(
+            influence_map: qec.DagFaultInfluenceMap,
+            detectors: Sequence[Sequence[int]],
+            observables: Sequence[Sequence[int]],
+            p1: float,
+            p2: float,
+            p_meas: float,
+            p_prep: float,
+            **noise_options: Any,
+        ) -> qec.DemSampler: ...
+        @staticmethod
+        def from_influence_map(
+            influence_map: qec.DagFaultInfluenceMap,
+            p_error: float,
+        ) -> qec.DemSampler: ...
+        @staticmethod
+        def from_influence_map_circuit_noise(
+            influence_map: qec.DagFaultInfluenceMap,
+            p1: float,
+            p2: float,
+            p_meas: float,
+            p_prep: float,
+        ) -> qec.DemSampler: ...
+        @property
+        def num_mechanisms(self) -> int: ...
+        @property
+        def num_outputs(self) -> int: ...
+        @property
+        def num_detectors(self) -> int: ...
+        @property
+        def num_observables(self) -> int: ...
+        @property
+        def num_dem_outputs(self) -> int: ...
+        @property
+        def num_tracked_paulis(self) -> int: ...
+        def sample(self, seed: int | None = ...) -> tuple[list[bool], list[bool]]: ...
+        def sample_batch(self, num_shots: int, seed: int | None = ...) -> qec.SampleBatch: ...
+        def sample_batch_with_pauli_masks(
+            self,
+            num_shots: int,
+            lookup: qec.PauliFrameLookup,
+            pauli_masks: Any,
+            seed: int | None = ...,
+        ) -> qec.SampleBatch: ...
+        def sample_tracked_paulis(self, seed: int | None = ...) -> list[bool]: ...
+        def sample_tracked_pauli_batch(self, num_shots: int, seed: int | None = ...) -> list[list[bool]]: ...
+        def sample_statistics(self, num_shots: int, seed: int | None = ...) -> dict[str, Any]: ...
+        def sample_decode_count(
+            self,
+            dem: str,
+            num_shots: int,
+            decoder_type: str = ...,
+            seed: int | None = ...,
+        ) -> int: ...
+        def sample_decode_count_parallel(
+            self,
+            dem: str,
+            num_shots: int,
+            decoder_type: str = ...,
+            seed: int | None = ...,
+            num_workers: int | None = ...,
+        ) -> int: ...
+        def labels(self) -> dict[str, Any]: ...
+
+    class DemSamplerBuilder:
+        def __init__(self, influence_map: qec.DagFaultInfluenceMap) -> None: ...
+        def with_noise(
+            self,
+            p1: float,
+            p2: float,
+            p_meas: float,
+            p_prep: float,
+            **noise_options: Any,
+        ) -> qec.DemSamplerBuilder: ...
+        def with_detectors_json(self, json: str) -> qec.DemSamplerBuilder: ...
+        def with_observables_json(self, json: str) -> qec.DemSamplerBuilder: ...
+        def with_measurement_order(self, order: Sequence[int]) -> qec.DemSamplerBuilder: ...
+        def build(self) -> qec.DemSampler: ...
+
+    class EquivalenceResult:
+        equivalent: bool
+        max_rate_difference: float
+        max_relative_difference: float
+        correlation: float
+        syndrome_rate_correlation: float
+        detector_rate_differences: list[float]
+        observable_rate_differences: list[float]
+        dem1_mechanism_count: int
+        dem2_mechanism_count: int
+        only_in_dem1: list[str]
+        only_in_dem2: list[str]
+        def details(self) -> dict[str, Any]: ...
+
+    class ParsedDem:
+        @staticmethod
+        def from_string(dem_str: str) -> qec.ParsedDem: ...
+        @property
+        def num_mechanisms(self) -> int: ...
+        @property
+        def num_detectors(self) -> int: ...
+        @property
+        def num_observables(self) -> int: ...
+        @property
+        def num_dem_outputs(self) -> int: ...
+        @property
+        def num_tracked_paulis(self) -> int: ...
+        def to_string_decomposed(self) -> str: ...
+        def aggregate(self) -> dict[tuple[tuple[int, ...], tuple[int, ...]], float]: ...
+        def sample(self, seed: int | None = ...) -> tuple[list[bool], list[bool]]: ...
+        def sample_batch(self, num_shots: int, seed: int | None = ...) -> qec.SampleBatch: ...
+        def to_dem_sampler(self) -> qec.DemSampler: ...
+
+    class CssUfDecoder:
+        def __init__(self, x_dem: str, z_dem: str) -> None: ...
+        @property
+        def num_qubit_pairs(self) -> int: ...
+        def decode_css(self, x_syndrome: Sequence[int], z_syndrome: Sequence[int]) -> tuple[int, int]: ...
+        def count_erasures(self, x_syndrome: Sequence[int], z_syndrome: Sequence[int]) -> int: ...
+        def decode_count_batch(
+            self,
+            syndromes: Sequence[Sequence[int]],
+            true_obs_masks: Sequence[int],
+            x_num_detectors: int,
+        ) -> int: ...
+
+    class LogicalSubgraphDecoder:
+        def __init__(
+            self,
+            dem: str,
+            stab_coords: Sequence[Mapping[str, Any]],
+            inner_decoder: str = ...,
+            max_time_radius: int | None = ...,
+        ) -> None: ...
+        @staticmethod
+        def from_membership(
+            dem: str,
+            membership: Sequence[Sequence[int]],
+            inner_decoder: str = ...,
+        ) -> qec.LogicalSubgraphDecoder: ...
+        @staticmethod
+        def count_ghost_edges(dem: str, stab_coords: Sequence[Mapping[str, Any]]) -> tuple[int, int]: ...
+        @property
+        def inner_decoder(self) -> str: ...
+        def decode(self, syndrome: Sequence[int]) -> int: ...
+        def decode_batch(self, syndromes: Sequence[Sequence[int]]) -> list[int]: ...
+        def decode_count(self, batch: qec.SampleBatch) -> int: ...
+        def decode_count_parallel(
+            self,
+            batch: qec.SampleBatch,
+            dem: str,
+            stab_coords: Sequence[Mapping[str, Any]],
+            inner_decoder: str | None = ...,
+            num_workers: int | None = ...,
+            max_time_radius: int | None = ...,
+        ) -> int: ...
+        def num_observables(self) -> int: ...
+        def observing_regions(self) -> list[Any]: ...
+        def subgraph_dems(self) -> list[str]: ...
+        def subgraph_detector_maps(self) -> list[list[int]]: ...
+        def subgraph_diagnostics(self) -> list[Any]: ...
+        def subgraph_sizes(self) -> list[int]: ...
+
+    class WindowedLogicalSubgraphDecoder:
+        def __init__(
+            self,
+            dem: str,
+            stab_coords: Sequence[Mapping[str, Any]],
+            step: int = ...,
+            buffer: int = ...,
+        ) -> None: ...
+        def decode(self, syndrome: Sequence[int]) -> int: ...
+        def decode_count(self, batch: qec.SampleBatch) -> int: ...
+        def num_windows(self) -> int: ...
+
+    class LogicalAlgorithmDecoder:
+        def __init__(self, descriptor: Mapping[str, Any], inner_decoder: str = ...) -> None: ...
+        def feed_dense(self, syndrome: Sequence[int]) -> None: ...
+        def feed_sparse(self, detectors: Sequence[int]) -> None: ...
+        def flush(self) -> int: ...
+        def decode(self, syndrome: Sequence[int]) -> int: ...
+        def decode_count(self, batch: qec.SampleBatch) -> int: ...
+        def reset(self) -> None: ...
+        def accumulated_obs(self) -> list[bool]: ...
+        def accumulated_obs_mask(self) -> int: ...
+        def rounds_fed(self) -> int: ...
+        def num_segments(self) -> int: ...
+
+    class LogicalCircuitDecoder:
+        def __init__(
+            self,
+            descriptor: Mapping[str, Any],
+            budget: str = ...,
+            inner_decoder: str = ...,
+            strict: bool = ...,
+        ) -> None: ...
+        def decode(self, syndrome: Sequence[int]) -> int: ...
+        def decode_count(self, batch: qec.SampleBatch) -> int: ...
+        def reset(self) -> None: ...
+        def num_segments(self) -> int: ...
+        def can_window(self) -> bool: ...
+        def actual_num_windows(self) -> int: ...
+        def effective_windowing(self) -> bool: ...
+        def has_decision_points(self) -> bool: ...
+        def num_decision_points(self) -> int: ...
+        def total_detectors(self) -> int: ...
+
+    @staticmethod
+    def compare_dems_exact(
+        dem1: str | qec.ParsedDem, dem2: str | qec.ParsedDem, prob_tolerance: float = ...
+    ) -> qec.EquivalenceResult: ...
+    @staticmethod
+    def compare_dems_statistical(
+        dem1: str | qec.ParsedDem,
+        dem2: str | qec.ParsedDem,
+        num_shots: int = ...,
+        seed: int = ...,
+        tolerance: float = ...,
+    ) -> qec.EquivalenceResult: ...
+    @staticmethod
+    def verify_dem_equivalence(
+        dem1: str | qec.ParsedDem,
+        dem2: str | qec.ParsedDem,
+        method: str = ...,
+        prob_tolerance: float = ...,
+        num_shots: int = ...,
+        tolerance: float = ...,
+        seed: int = ...,
+    ) -> qec.EquivalenceResult: ...
+    @staticmethod
+    def assert_dems_equivalent(
+        dem1: str | qec.ParsedDem,
+        dem2: str | qec.ParsedDem,
+        method: str = ...,
+        prob_tolerance: float = ...,
+        num_shots: int = ...,
+        tolerance: float = ...,
+        seed: int = ...,
+    ) -> None: ...
+    @staticmethod
+    def detector_flip_matrix(fired_per_shot: Sequence[Sequence[int]], num_detectors: int) -> list[list[float]]: ...
+    @staticmethod
+    def detector_flip_matrices_by_round(
+        fired_per_shot: Sequence[Sequence[int]],
+        num_detectors: int,
+        dets_per_round: int,
+    ) -> list[Any]: ...
+    @staticmethod
+    def detector_k_body_rates(
+        fired_per_shot: Sequence[Sequence[int]],
+        num_detectors: int,
+        max_order: int = ...,
+    ) -> dict[Any, float]: ...
+    @staticmethod
+    def detector_k_body_rates_by_round(
+        fired_per_shot: Sequence[Sequence[int]],
+        num_detectors: int,
+        dets_per_round: int,
+        max_order: int = ...,
+    ) -> list[Any]: ...
+    @staticmethod
+    def compare_flip_matrices_rs(sim: Any, dem: Any, num_detectors: int, min_rate: float = ...) -> dict[str, Any]: ...
+    @staticmethod
+    def compare_k_body_rates_rs(sim: Any, dem: Any, min_rate: float = ...) -> dict[str, Any]: ...
+    @staticmethod
+    def fit_dem_to_marginals(
+        mechanisms: Sequence[tuple[float, Sequence[int], Sequence[int]]],
+        target_marginals: Sequence[float],
+        max_iterations: int = ...,
+        tolerance: float = ...,
+    ) -> tuple[list[float], list[float]]: ...
+    @staticmethod
+    def mechanisms_to_dem_string(mechanisms: Sequence[tuple[float, Sequence[int], Sequence[int]]]) -> str: ...
+    @staticmethod
+    def decoder_dem_requirement(decoder_type: str) -> str: ...
+
+DemSampler = qec.DemSampler
+
 # =============================================================================
 # Decoder Types
 # =============================================================================
diff --git a/python/pecos-rslib/src/fault_tolerance_bindings.rs b/python/pecos-rslib/src/fault_tolerance_bindings.rs
index fbb855d8b..0b9f4a899 100644
--- a/python/pecos-rslib/src/fault_tolerance_bindings.rs
+++ b/python/pecos-rslib/src/fault_tolerance_bindings.rs
@@ -55,8 +55,8 @@ use pecos_qec::fault_tolerance::dem_builder::{
     DetectorErrorModel as RustDetectorErrorModel, DirectSourceFamily as RustDirectSourceFamily,
     EquivalenceResult as RustEquivalenceResult, FaultContribution as RustFaultContribution,
     FaultSourceType as RustFaultSourceType, MeasurementCrosstalkDemMode,
-    MeasurementCrosstalkTransitionModel, NoiseConfig, PAULI_2Q_ORDER, ParsedDem as RustParsedDem,
-    PauliWeights, ReplacementBranchApproximation,
+    MeasurementCrosstalkTransitionModel, NoiseConfig, OutputMode, PAULI_2Q_ORDER,
+    ParsedDem as RustParsedDem, PauliWeights, ReplacementBranchApproximation,
     TwoDetectorDirectRenderPolicy as RustTwoDetectorDirectRenderPolicy,
     compare_dems_exact as rust_compare_dems_exact,
     compare_dems_statistical as rust_compare_dems_statistical,
@@ -3348,14 +3348,19 @@ fn create_observable_decoder(
 
 /// Pre-generated sample batch held in Rust memory.
 ///
-/// Created by `DemSampler.generate_samples()`. Can be decoded by multiple
+/// Created by `DemSampler.sample_batch()`. Can be decoded by multiple
 /// decoders without re-sampling, and without crossing the Rust/Python boundary
 /// per shot.
 ///
+/// A batch produced by a raw-measurement `DemSampler` uses the same container,
+/// but its detector columns contain raw measurements rather than detector
+/// events. Data accessors remain available for those batches; decode methods
+/// reject them because raw measurements are not decoder syndromes.
+///
 /// # Example
 ///
 /// ```python
-/// samples = sampler.generate_samples(10000, seed=42)
+/// samples = sampler.sample_batch(10000, seed=42)
 /// pm_errors = samples.decode_count(dem, "pymatching")
 /// ts_errors = samples.decode_count(dem, "tesseract")
 /// # Both decoders ran on the exact same samples.
@@ -3368,6 +3373,7 @@ pub struct PySampleBatch {
     obs_columns: Vec<Vec<u64>>,
     num_detectors: usize,
     num_shots: usize,
+    raw_measurements: bool,
 }
 
 impl PySampleBatch {
@@ -3400,6 +3406,16 @@ impl PySampleBatch {
         Ok(())
     }
 
+    /// Reject raw-measurement batches before treating their rows as syndromes.
+    fn ensure_detector_events(&self) -> PyResult<()> {
+        if self.raw_measurements {
+            return Err(pyo3::exceptions::PyValueError::new_err(
+                "raw-measurement SampleBatch rows carry measurements, not detector events, and cannot be decoded",
+            ));
+        }
+        Ok(())
+    }
+
     /// Extract observable mask for one shot (`u64`; observables 0..=63 only).
     ///
     /// The caller must have rejected wide batches via
@@ -3435,7 +3451,7 @@ impl PySampleBatch {
         mask
     }
 
-    /// Build from columnar data (from generate_samples).
+    /// Build from columnar sampling data.
     fn from_columnar(
         det_columns: Vec<Vec<u64>>,
         obs_columns: Vec<Vec<u64>>,
@@ -3447,7 +3463,73 @@ impl PySampleBatch {
             obs_columns,
             num_detectors,
             num_shots,
+            raw_measurements: false,
+        }
+    }
+
+    /// Build from rectangular row-major boolean detector and observable data.
+    ///
+    /// Both outer lists must have equal length, and each list's rows must have
+    /// the same width as its row 0.
+    fn from_bool_rows(
+        detection_events: Vec<Vec<bool>>,
+        observable_flips: Vec<Vec<bool>>,
+        raw_measurements: bool,
+    ) -> Self {
+        debug_assert_eq!(observable_flips.len(), detection_events.len());
+        let num_shots = detection_events.len();
+        let num_detectors = detection_events.first().map_or(0, Vec::len);
+        let num_observables = observable_flips.first().map_or(0, Vec::len);
+        debug_assert!(
+            detection_events
+                .iter()
+                .all(|row| row.len() == num_detectors)
+        );
+        debug_assert!(
+            observable_flips
+                .iter()
+                .all(|row| row.len() == num_observables)
+        );
+        let num_words = num_shots.div_ceil(64);
+        let mut det_columns = vec![vec![0u64; num_words]; num_detectors];
+        let mut obs_columns = vec![vec![0u64; num_words]; num_observables];
+
+        for (shot, row) in detection_events.iter().enumerate() {
+            let word_idx = shot / 64;
+            let bit_mask = 1u64 << (shot % 64);
+            for (det_idx, &value) in row.iter().enumerate() {
+                if value {
+                    det_columns[det_idx][word_idx] |= bit_mask;
+                }
+            }
         }
+        for (shot, row) in observable_flips.iter().enumerate() {
+            let word_idx = shot / 64;
+            let bit_mask = 1u64 << (shot % 64);
+            for (obs_idx, &value) in row.iter().enumerate() {
+                if value {
+                    obs_columns[obs_idx][word_idx] |= bit_mask;
+                }
+            }
+        }
+
+        let mut batch = Self::from_columnar(det_columns, obs_columns, num_shots);
+        batch.raw_measurements = raw_measurements;
+        batch
+    }
+
+    /// Materialize bit-packed columns as shots-major boolean rows.
+    fn columns_as_rows(columns: &[Vec<u64>], num_shots: usize) -> Vec<Vec<bool>> {
+        (0..num_shots)
+            .map(|shot| {
+                let word_idx = shot / 64;
+                let bit_mask = 1u64 << (shot % 64);
+                columns
+                    .iter()
+                    .map(|column| column[word_idx] & bit_mask != 0)
+                    .collect()
+            })
+            .collect()
     }
 
     /// Build from row-major data (from Python constructor). Observable masks are
@@ -3455,6 +3537,7 @@ impl PySampleBatch {
     fn from_row_major(
         detection_events: Vec<Vec<u8>>,
         observable_masks: &[pecos_decoder_core::obs_mask::ObsMask],
+        num_observables: usize,
     ) -> Self {
         let num_shots = detection_events.len();
         let num_detectors = detection_events.first().map_or(0, Vec::len);
@@ -3472,14 +3555,7 @@ impl PySampleBatch {
             }
         }
 
-        // One observable column per observable index; sized to the highest set
-        // bit across all shots (supports >64 observables).
-        let max_obs = observable_masks
-            .iter()
-            .filter_map(|m| m.iter_set_bits().max())
-            .max()
-            .map_or(0, |b| b + 1);
-        let mut obs_columns = vec![vec![0u64; num_words]; max_obs];
+        let mut obs_columns = vec![vec![0u64; num_words]; num_observables];
         for (shot, mask) in observable_masks.iter().enumerate() {
             let word_idx = shot / 64;
             let bit_mask = 1u64 << (shot % 64);
@@ -3493,6 +3569,7 @@ impl PySampleBatch {
             obs_columns,
             num_detectors,
             num_shots,
+            raw_measurements: false,
         }
     }
 }
@@ -3506,11 +3583,16 @@ impl PySampleBatch {
     ///     observable_masks: List of true observable flip masks as Python ints
     ///         (arbitrary precision; bit ``i`` = observable ``i``, so more than 64
     ///         observables are supported).
+    ///     num_observables: Optional exact observable-column width. Every set
+    ///         mask bit must be below this width. When omitted, the width is
+    ///         inferred as one greater than the highest set bit across all masks;
+    ///         consequently, all-zero masks infer zero observable columns.
     #[new]
-    #[pyo3(signature = (detection_events, observable_masks))]
+    #[pyo3(signature = (detection_events, observable_masks, *, num_observables=None))]
     fn new(
         detection_events: Vec<Vec<u8>>,
         observable_masks: Vec<pyo3::Bound<'_, pyo3::PyAny>>,
+        num_observables: Option<usize>,
     ) -> PyResult<Self> {
         if detection_events.len() != observable_masks.len() {
             return Err(pyo3::exceptions::PyValueError::new_err(format!(
@@ -3533,7 +3615,29 @@ impl PySampleBatch {
             .iter()
             .map(py_to_obsmask)
             .collect::<PyResult<_>>()?;
-        Ok(Self::from_row_major(detection_events, &masks))
+        let observable_width = if let Some(width) = num_observables {
+            if let Some(bit) = masks
+                .iter()
+                .flat_map(pecos_decoder_core::obs_mask::ObsMask::iter_set_bits)
+                .find(|&bit| bit >= width)
+            {
+                return Err(pyo3::exceptions::PyValueError::new_err(format!(
+                    "observable mask bit {bit} is outside num_observables={width}",
+                )));
+            }
+            width
+        } else {
+            masks
+                .iter()
+                .filter_map(|mask| mask.iter_set_bits().max())
+                .max()
+                .map_or(0, |bit| bit + 1)
+        };
+        Ok(Self::from_row_major(
+            detection_events,
+            &masks,
+            observable_width,
+        ))
     }
 
     /// Number of shots in this batch.
@@ -3579,6 +3683,25 @@ impl PySampleBatch {
         obsmask_to_py(py, &self.extract_obs_mask_wide(i))
     }
 
+    /// Return all detector events as shots-major boolean lists.
+    ///
+    /// The result has shape (`num_shots`, `num_detectors`).
+    fn detector_events(&self) -> Vec<Vec<bool>> {
+        Self::columns_as_rows(&self.det_columns, self.num_shots)
+    }
+
+    /// Return all observable flips as shots-major boolean lists.
+    ///
+    /// The result has shape (`num_shots`, stored observable-column width) and
+    /// does not truncate batches containing more than 64 observables. For the
+    /// Python constructor, the width is `num_observables` when supplied and is
+    /// otherwise inferred from the highest set mask bit (all-zero masks infer
+    /// width zero). Sampler-produced columns contain all DEM outputs, which can
+    /// be a superset of the logical observables.
+    fn observable_flips(&self) -> Vec<Vec<bool>> {
+        Self::columns_as_rows(&self.obs_columns, self.num_shots)
+    }
+
     /// Decode all samples with the given decoder type and return the error count.
     ///
     /// This runs entirely in Rust -- no per-shot Python crossing.
@@ -3593,6 +3716,7 @@ impl PySampleBatch {
     ///     Number of logical errors.
     #[pyo3(signature = (dem, decoder_type="pymatching"))]
     fn decode_count(&self, dem: &str, decoder_type: &str) -> PyResult<usize> {
+        self.ensure_detector_events()?;
         let mut decoder = create_observable_decoder(dem, decoder_type)?;
         let mut errors = 0usize;
         let mut syndrome = vec![0u8; self.num_detectors];
@@ -3631,6 +3755,7 @@ impl PySampleBatch {
         dem: &str,
         decoder_type: &str,
     ) -> PyResult<Vec<Py<pyo3::PyAny>>> {
+        self.ensure_detector_events()?;
         let mut decoder = create_observable_decoder(dem, decoder_type)?;
         let mut predictions = Vec::with_capacity(self.num_shots);
         let mut syndrome = vec![0u8; self.num_detectors];
@@ -3666,6 +3791,7 @@ impl PySampleBatch {
     ) -> PyResult<usize> {
         use rayon::prelude::*;
 
+        self.ensure_detector_events()?;
         let n_workers = num_workers.unwrap_or_else(rayon::current_num_threads);
         let pool = rayon::ThreadPoolBuilder::new()
             .num_threads(n_workers)
@@ -3719,6 +3845,7 @@ impl PySampleBatch {
     fn decode_count_batch(&self, dem: &str) -> PyResult<usize> {
         use pecos_decoders::{BatchConfig, PyMatchingDecoder};
 
+        self.ensure_detector_events()?;
         let mut decoder = PyMatchingDecoder::from_dem(dem)
             .map_err(|e| PyErr::new::<pyo3::exceptions::PyRuntimeError, _>(e.to_string()))?;
 
@@ -3781,6 +3908,7 @@ impl PySampleBatch {
     fn decode_stats(&self, dem: &str, decoder_type: &str) -> PyResult<PyDecodeStats> {
         use std::time::Instant;
 
+        self.ensure_detector_events()?;
         let mut decoder = create_observable_decoder(dem, decoder_type)?;
         let mut num_errors = 0usize;
         let mut per_shot_seconds: Vec<f64> = Vec::with_capacity(self.num_shots);
@@ -3826,6 +3954,7 @@ impl PySampleBatch {
     ) -> PyResult<PyDecodeStats> {
         use rayon::prelude::*;
 
+        self.ensure_detector_events()?;
         let n_workers = num_workers.unwrap_or_else(rayon::current_num_threads);
 
         // Validate decoder type early.
@@ -4346,20 +4475,23 @@ impl PyDemSampler {
         self.inner.sample(&mut rng)
     }
 
-    /// Sample multiple shots.
+    /// Sample multiple shots into a `SampleBatch` held in Rust memory.
+    ///
+    /// The batch can be decoded by multiple decoders without re-sampling or
+    /// materialized as shots-major Python lists with
+    /// `SampleBatch.detector_events()` and `SampleBatch.observable_flips()`.
+    /// For a raw-measurement sampler, the first set of columns contains raw
+    /// measurements rather than detector events, so the batch data accessors
+    /// work but its decode methods raise `ValueError`.
     ///
     /// Args:
-    ///     `num_shots`: Number of shots to sample.
+    ///     num_shots: Number of shots to sample.
     ///     seed: Optional random seed for reproducibility.
     ///
     /// Returns:
-    ///     Tuple of (`all_detection_events`, `all_dem_output_flips`).
+    ///     `SampleBatch` object with samples held in Rust memory.
     #[pyo3(signature = (num_shots, seed=None))]
-    fn sample_batch(
-        &self,
-        num_shots: usize,
-        seed: Option<u64>,
-    ) -> (Vec<Vec<bool>>, Vec<Vec<bool>>) {
+    fn sample_batch(&self, num_shots: usize, seed: Option<u64>) -> PySampleBatch {
         use pecos_random::PecosRng;
         use rand::RngExt;
 
@@ -4368,20 +4500,25 @@ impl PyDemSampler {
             None => PecosRng::seed_from_u64(rand::rng().random()),
         };
 
-        self.inner.sample_batch(num_shots, &mut rng)
+        if self.inner.mode() == OutputMode::RawMeasurements {
+            let (detection_events, observable_flips) = self.inner.sample_batch(num_shots, &mut rng);
+            return PySampleBatch::from_bool_rows(detection_events, observable_flips, true);
+        }
+        let (det_columns, obs_columns) = self.inner.sample_batch_geometric(num_shots, &mut rng);
+        PySampleBatch::from_columnar(det_columns, obs_columns, num_shots)
     }
 
     /// Sample multiple shots and XOR a known Pauli-frame mask into the outputs.
     ///
     /// Args:
-    ///     `num_shots`: Number of shots to sample.
+    ///     num_shots: Number of shots to sample.
     ///     lookup: Pauli-frame lookup built from the same circuit metadata.
-    ///     `pauli_masks`: Integer array with shape `(num_shots, num_pauli_sites)`.
+    ///     pauli_masks: Integer array with shape `(num_shots, num_pauli_sites)`.
     ///         Values are 0=I, 1=X, 2=Y, 3=Z.
     ///     seed: Optional random seed for reproducibility.
     ///
     /// Returns:
-    ///     Tuple of (`all_detection_events`, `all_dem_output_flips`).
+    ///     `SampleBatch` containing the sampled and XOR-adjusted outputs.
     #[pyo3(signature = (num_shots, lookup, pauli_masks, seed=None))]
     fn sample_batch_with_pauli_masks(
         &self,
@@ -4389,7 +4526,7 @@ impl PyDemSampler {
         lookup: &PyPauliFrameLookup,
         pauli_masks: &Bound<'_, pyo3::PyAny>,
         seed: Option<u64>,
-    ) -> PyResult<PyDetectorObservableRows> {
+    ) -> PyResult<PySampleBatch> {
         use pecos_random::PecosRng;
         use rand::RngExt;
 
@@ -4425,7 +4562,11 @@ impl PyDemSampler {
                 &mut obs_flips,
             )
             .map_err(|e| pyo3::exceptions::PyValueError::new_err(e.to_string()))?;
-        Ok((det_events, obs_flips))
+        Ok(PySampleBatch::from_bool_rows(
+            det_events,
+            obs_flips,
+            self.inner.mode() == OutputMode::RawMeasurements,
+        ))
     }
 
     /// Sample direct tracked-Pauli flips.
@@ -4472,33 +4613,6 @@ impl PyDemSampler {
             .map_err(|e| PyErr::new::<pyo3::exceptions::PyRuntimeError, _>(e.to_string()))
     }
 
-    /// Generate samples and store them in Rust memory as a `SampleBatch`.
-    ///
-    /// The batch can then be decoded by multiple decoders without re-sampling.
-    /// This is the proper way to compare decoders: same samples, different decoders.
-    ///
-    /// Args:
-    ///     `num_shots`: Number of shots to sample.
-    ///     seed: Optional random seed for reproducibility.
-    ///
-    /// Returns:
-    ///     `SampleBatch` object with samples held in Rust memory.
-    #[pyo3(signature = (num_shots, seed=None))]
-    fn generate_samples(&self, num_shots: usize, seed: Option<u64>) -> PySampleBatch {
-        use pecos_random::PecosRng;
-        use rand::RngExt;
-
-        let mut rng = match seed {
-            Some(s) => PecosRng::seed_from_u64(s),
-            None => PecosRng::seed_from_u64(rand::rng().random()),
-        };
-
-        // Use geometric columnar sampler via DemSampler.
-        let (det_columns, obs_columns) = self.inner.sample_batch_geometric(num_shots, &mut rng);
-
-        PySampleBatch::from_columnar(det_columns, obs_columns, num_shots)
-    }
-
     /// Compute statistics without storing individual shots.
     ///
     /// This is the most efficient method for threshold estimation when you
@@ -5153,20 +5267,16 @@ impl PyParsedDem {
         self.inner.sample(&mut rng)
     }
 
-    /// Sample multiple shots from this DEM.
+    /// Sample multiple shots from this DEM into a `SampleBatch`.
     ///
     /// Args:
-    ///     `num_shots`: Number of shots to sample.
+    ///     num_shots: Number of shots to sample.
     ///     seed: Optional random seed for reproducibility.
     ///
     /// Returns:
-    ///     Tuple of (`all_detector_events`, `all_dem_output_flips`).
+    ///     `SampleBatch` object with samples held in Rust memory.
     #[pyo3(signature = (num_shots, seed=None))]
-    fn sample_batch(
-        &self,
-        num_shots: usize,
-        seed: Option<u64>,
-    ) -> (Vec<Vec<bool>>, Vec<Vec<bool>>) {
+    fn sample_batch(&self, num_shots: usize, seed: Option<u64>) -> PySampleBatch {
         use pecos_random::PecosRng;
         use rand::RngExt;
 
@@ -5175,7 +5285,8 @@ impl PyParsedDem {
             None => PecosRng::seed_from_u64(rand::rng().random()),
         };
 
-        self.inner.sample_batch(num_shots, &mut rng)
+        let (detector_events, observable_flips) = self.inner.sample_batch(num_shots, &mut rng);
+        PySampleBatch::from_bool_rows(detector_events, observable_flips, false)
     }
 
     /// Convert to an optimized `DemSampler` for fast batch sampling.
@@ -5734,11 +5845,12 @@ impl PyLogicalSubgraphDecoder {
     /// This runs entirely in Rust — no Python per-shot overhead.
     ///
     /// Args:
-    ///     batch: A `SampleBatch` from `DemSampler.generate_samples()`.
+    ///     batch: A `SampleBatch` from `DemSampler.sample_batch()`.
     ///
     /// Returns:
     ///     Number of logical errors.
     fn decode_count(&mut self, batch: &PySampleBatch) -> PyResult<usize> {
+        batch.ensure_detector_events()?;
         let detection_events: Vec<Vec<u8>> = (0..batch.num_shots)
             .map(|i| {
                 let mut s = vec![0u8; batch.num_detectors];
@@ -5773,6 +5885,7 @@ impl PyLogicalSubgraphDecoder {
         use pecos_decoder_core::logical_subgraph::{LogicalSubgraphDecoder, QubitStabCoords};
         use rayon::prelude::*;
 
+        batch.ensure_detector_events()?;
         // Parse stab_coords
         let mut sc = Vec::with_capacity(stab_coords.len());
         for dict in &stab_coords {
@@ -6021,6 +6134,7 @@ impl PyWindowedLogicalSubgraphDecoder {
 
     fn decode_count(&mut self, batch: &PySampleBatch) -> PyResult<usize> {
         use pecos_decoder_core::ObservableDecoder;
+        batch.ensure_detector_events()?;
         let mut errors = 0usize;
         let mut syndrome = vec![0u8; batch.num_detectors];
         for i in 0..batch.num_shots {
@@ -6248,6 +6362,7 @@ impl PyLogicalAlgorithmDecoder {
 
     /// Decode a batch of samples and count logical errors (wide observable masks).
     fn decode_count(&mut self, batch: &PySampleBatch) -> PyResult<usize> {
+        batch.ensure_detector_events()?;
         let mut errors = 0usize;
         let mut syndrome = vec![0u8; batch.num_detectors];
         for i in 0..batch.num_shots {
@@ -6665,6 +6780,7 @@ impl PyLogicalCircuitDecoder {
     /// Decode a batch and count errors (wide observable masks).
     fn decode_count(&mut self, batch: &PySampleBatch) -> PyResult<usize> {
         use pecos_decoder_core::ObservableDecoder;
+        batch.ensure_detector_events()?;
         let mut errors = 0usize;
         let mut syndrome = vec![0u8; batch.num_detectors];
         for i in 0..batch.num_shots {
diff --git a/python/quantum-pecos/src/pecos/qec/surface/decode.py b/python/quantum-pecos/src/pecos/qec/surface/decode.py
index 6cac5f9c3..c17455d2f 100644
--- a/python/quantum-pecos/src/pecos/qec/surface/decode.py
+++ b/python/quantum-pecos/src/pecos/qec/surface/decode.py
@@ -3482,7 +3482,7 @@ def surface_code_memory(
         require_hosted_operation_order=require_hosted_operation_order,
         max_hosted_tick_separation=max_hosted_tick_separation,
     )
-    batch = ParsedDem.from_string(dem).to_dem_sampler().generate_samples(shots, seed)
+    batch = ParsedDem.from_string(dem).to_dem_sampler().sample_batch(shots, seed)
     num_raw_errors = sum(1 for shot in range(shots) if batch.get_observable_mask(shot) != 0)
     num_logical_errors = batch.decode_count(dem, decoder_type) if decode else num_raw_errors
 
@@ -3751,19 +3751,22 @@ def sample(
             - observable_flips: shape (num_shots, num_observables)
         """
         if pauli_masks is None:
-            det_events, obs_flips = self.sampler.sample_batch(num_shots, seed)
+            batch = self.sampler.sample_batch(num_shots, seed)
         else:
             if self.pauli_frame_lookup is None:
                 msg = "pauli_masks require build_native_sampler(..., twirl=TwirlConfig())"
                 raise ValueError(msg)
             masks_arr = _pauli_masks_as_int64(pauli_masks)
-            det_events, obs_flips = self.sampler.sample_batch_with_pauli_masks(
+            batch = self.sampler.sample_batch_with_pauli_masks(
                 num_shots,
                 self.pauli_frame_lookup,
                 masks_arr,
                 seed,
             )
-        return np.array(det_events, dtype=bool), np.array(obs_flips, dtype=bool)
+        return (
+            np.array(batch.detector_events(), dtype=bool),
+            np.array(batch.observable_flips(), dtype=bool),
+        )
 
 
 def build_native_sampler(
@@ -4061,7 +4064,7 @@ def decode_native_samples(
         raise ValueError(msg)
     weights = (1 << np.arange(obs_arr.shape[1], dtype=np.uint64)).astype(np.uint64)
     obs_masks = (obs_arr * weights).sum(axis=1).astype(np.uint64).tolist()
-    batch = SampleBatch(det_list, obs_masks)
+    batch = SampleBatch(det_list, obs_masks, num_observables=obs_arr.shape[1])
     return batch.decode_count(dem_str, decoder_type)
 
 
diff --git a/python/quantum-pecos/tests/qec/surface/test_circuit_fuzz.py b/python/quantum-pecos/tests/qec/surface/test_circuit_fuzz.py
index e8fd46fb6..7e19847d0 100644
--- a/python/quantum-pecos/tests/qec/surface/test_circuit_fuzz.py
+++ b/python/quantum-pecos/tests/qec/surface/test_circuit_fuzz.py
@@ -588,7 +588,7 @@ def test_pecos_dem_decode_memory(self, patch):
 
         parsed = ParsedDem.from_string(dem_str)
         rust_sampler = parsed.to_dem_sampler()
-        batch = rust_sampler.generate_samples(5000, seed=42)
+        batch = rust_sampler.sample_batch(5000, seed=42)
         errors = batch.decode_count(dem_str, "pecos_uf:fast")
         ler = errors / 5000
         # At d=3 p=0.001, LER should be very low
@@ -628,7 +628,7 @@ def test_pecos_dem_cx_decode(self, patch, nq):
         dem_str = b.build_dem(p1=0.001, p2=0.001, p_meas=0.001)
         parsed = ParsedDem.from_string(dem_str)
         rust_sampler = parsed.to_dem_sampler()
-        batch = rust_sampler.generate_samples(5000, seed=42)
+        batch = rust_sampler.sample_batch(5000, seed=42)
         errors = batch.decode_count(dem_str, "pecos_uf:fast")
         ler = errors / 5000
         assert ler < 0.1, f"CX LER too high: {ler}"
@@ -660,7 +660,7 @@ def _run_threshold(self, builder, _d, decoder_type="pecos_uf:fast"):
         dem_str = str(dem)
         parsed = ParsedDem.from_string(dem_str)
         sampler = parsed.to_dem_sampler()
-        batch = sampler.generate_samples(20000, seed=42)
+        batch = sampler.sample_batch(20000, seed=42)
         errors = batch.decode_count(dem_str, decoder_type)
         return errors / 20000
 
@@ -726,7 +726,7 @@ def test_noisy_h(self, patch, seed):
         dem = c.detector_error_model(decompose_errors=True)
         dem_str = str(dem)
         parsed = ParsedDem.from_string(dem_str)
-        batch = parsed.to_dem_sampler().generate_samples(5000, seed=seed)
+        batch = parsed.to_dem_sampler().sample_batch(5000, seed=seed)
         errors = batch.decode_count(dem_str, "pecos_uf:fast")
         ler = errors / 5000
         assert ler < 0.1, f"LER too high: {ler}"
@@ -767,7 +767,7 @@ def test_noisy_h_cx_h(self, patch, nq):
         dem = c.detector_error_model(decompose_errors=True, ignore_decomposition_failures=True)
         dem_str = str(dem)
         parsed = ParsedDem.from_string(dem_str)
-        batch = parsed.to_dem_sampler().generate_samples(10000, seed=42)
+        batch = parsed.to_dem_sampler().sample_batch(10000, seed=42)
         errors = batch.decode_count(dem_str, "pecos_uf:fast")
         ler = errors / 10000
         assert ler < 0.1, f"H-CX-H LER too high: {ler}"
@@ -805,7 +805,7 @@ def test_noisy_random_composition(self, patch, nq, seed):
         dem = c.detector_error_model(decompose_errors=True, ignore_decomposition_failures=True)
         dem_str = str(dem)
         parsed = ParsedDem.from_string(dem_str)
-        batch = parsed.to_dem_sampler().generate_samples(5000, seed=42)
+        batch = parsed.to_dem_sampler().sample_batch(5000, seed=42)
         errors = batch.decode_count(dem_str, "pecos_uf:fast")
         ler = errors / 5000
         assert ler < 0.2, f"Random composition LER too high: {ler}"
@@ -842,7 +842,7 @@ def test_logical_subgraph_better_than_naive_on_cx(self, patch, nq):
         # Naive: decomposed MWPM via PECOS UF
         dem_decomp = c.detector_error_model(decompose_errors=True, ignore_decomposition_failures=True)
         parsed = ParsedDem.from_string(str(dem_decomp))
-        batch_naive = parsed.to_dem_sampler().generate_samples(20000, seed=42)
+        batch_naive = parsed.to_dem_sampler().sample_batch(20000, seed=42)
         naive_errors = batch_naive.decode_count(str(dem_decomp), "pecos_uf:fast")
         naive_ler = naive_errors / 20000
 
@@ -895,7 +895,7 @@ def test_pecos_dem_logical_subgraph_cx(self, patch, nq):
 
         # Sample and decode
         parsed = ParsedDem.from_string(dem_str)
-        batch = parsed.to_dem_sampler().generate_samples(5000, seed=42)
+        batch = parsed.to_dem_sampler().sample_batch(5000, seed=42)
         errors = batch.decode_count(dem_str, "pecos_uf:fast")
         ler = errors / 5000
         assert ler < 0.1, f"PECOS DEM + logical-subgraph decoder CX LER too high: {ler}"
diff --git a/python/quantum-pecos/tests/qec/surface/test_logical_subgraph_region_comparison.py b/python/quantum-pecos/tests/qec/surface/test_logical_subgraph_region_comparison.py
index 2fb6487b7..cc06ef5c7 100644
--- a/python/quantum-pecos/tests/qec/surface/test_logical_subgraph_region_comparison.py
+++ b/python/quantum-pecos/tests/qec/surface/test_logical_subgraph_region_comparison.py
@@ -140,7 +140,7 @@ def test_from_membership_reproduces_coordinate_path():
     )
 
     assert rebuilt.subgraph_sizes() == coord.subgraph_sizes()
-    batch = ParsedDem.from_string(dem).to_dem_sampler().generate_samples(2000, seed=3)
+    batch = ParsedDem.from_string(dem).to_dem_sampler().sample_batch(2000, seed=3)
     assert rebuilt.decode_count(batch) == coord.decode_count(batch)
 
 
@@ -153,7 +153,7 @@ def test_coordinate_region_beats_raw_backprop_region():
     sc = b.stab_coords()
 
     n = 20000
-    batch = ParsedDem.from_string(dem).to_dem_sampler().generate_samples(n, seed=11)
+    batch = ParsedDem.from_string(dem).to_dem_sampler().sample_batch(n, seed=11)
 
     coord = LogicalSubgraphDecoder(dem, sc, "pecos_uf:fast")
     coflip_membership = _coflip_membership_from_dem(dem, coord.num_observables())
@@ -203,7 +203,7 @@ def test_coordinate_beats_backprop_seeded_groupfill():
     assert sum(len(bp) for bp in backprop_membership) > sum(len(c) for c in coord_membership)
 
     n = 20000
-    batch = ParsedDem.from_string(dem).to_dem_sampler().generate_samples(n, seed=13)
+    batch = ParsedDem.from_string(dem).to_dem_sampler().sample_batch(n, seed=13)
     coord = LogicalSubgraphDecoder.from_membership(dem, coord_membership, "pecos_uf:fast")
     backprop = LogicalSubgraphDecoder.from_membership(dem, backprop_membership, "pecos_uf:fast")
     coord_ler = coord.decode_count(batch) / n
@@ -218,7 +218,7 @@ def _mem_ler(d, p, n, seed, inner=None):
     b.add_memory("A", d, "Z")
     dem = b.build_dem(p1=p, p2=p, p_meas=p)
     sc = b.stab_coords()
-    batch = ParsedDem.from_string(dem).to_dem_sampler().generate_samples(n, seed=seed)
+    batch = ParsedDem.from_string(dem).to_dem_sampler().sample_batch(n, seed=seed)
     dec = LogicalSubgraphDecoder(dem, sc) if inner is None else LogicalSubgraphDecoder(dem, sc, inner)
     return dec.decode_count(batch) / n
 
@@ -270,7 +270,7 @@ def _mem_dem_batch(d, p, n, seed):
     b.add_memory("A", d, "Z")
     dem = b.build_dem(p1=p, p2=p, p_meas=p)
     sc = b.stab_coords()
-    batch = ParsedDem.from_string(dem).to_dem_sampler().generate_samples(n, seed=seed)
+    batch = ParsedDem.from_string(dem).to_dem_sampler().sample_batch(n, seed=seed)
     return dem, sc, batch
 
 
@@ -307,7 +307,7 @@ def _windowed_mem_ler(d, rounds, p, n, seed, step, buffer):
     b.add_memory("A", rounds, "Z")
     dem = b.build_dem(p1=p, p2=p, p_meas=p)
     sc = b.stab_coords()
-    batch = ParsedDem.from_string(dem).to_dem_sampler().generate_samples(n, seed=seed)
+    batch = ParsedDem.from_string(dem).to_dem_sampler().sample_batch(n, seed=seed)
     dec = WindowedLogicalSubgraphDecoder(dem, sc, step, buffer)
     return dec.decode_count(batch) / n, dec.num_windows()
 
@@ -319,7 +319,7 @@ def _nonwindowed_mem_ler(d, rounds, p, n, seed, inner):
     b.add_memory("A", rounds, "Z")
     dem = b.build_dem(p1=p, p2=p, p_meas=p)
     sc = b.stab_coords()
-    batch = ParsedDem.from_string(dem).to_dem_sampler().generate_samples(n, seed=seed)
+    batch = ParsedDem.from_string(dem).to_dem_sampler().sample_batch(n, seed=seed)
     return LogicalSubgraphDecoder(dem, sc, inner).decode_count(batch) / n
 
 
@@ -401,7 +401,7 @@ def test_decode_each_matches_decode_count():
     b = _cx_circuit()
     dem = b.build_dem(p1=0.001, p2=0.001, p_meas=0.001)
     n = 3000
-    batch = ParsedDem.from_string(dem).to_dem_sampler().generate_samples(n, seed=5)
+    batch = ParsedDem.from_string(dem).to_dem_sampler().sample_batch(n, seed=5)
     preds = batch.decode_each(dem, "pecos_uf:bp")
     assert len(preds) == n
     wrong = sum(1 for i, p in enumerate(preds) if p != batch.get_observable_mask(i))
diff --git a/python/quantum-pecos/tests/qec/surface/test_pauli_twirl_handoff.py b/python/quantum-pecos/tests/qec/surface/test_pauli_twirl_handoff.py
index 810b3907e..a0c44873e 100644
--- a/python/quantum-pecos/tests/qec/surface/test_pauli_twirl_handoff.py
+++ b/python/quantum-pecos/tests/qec/surface/test_pauli_twirl_handoff.py
@@ -28,6 +28,7 @@
     _extract_pauli_masks_from_results,
     generate_circuit_level_dem_from_builder,
 )
+from pecos_rslib.qec import SampleBatch
 
 
 def test_extract_pauli_masks_packs_bits_in_row_major_site_qubit_order() -> None:
@@ -228,6 +229,58 @@ def test_native_sampler_accepts_harvested_uint8_pauli_masks() -> None:
     assert decode_native_samples(sampler, 4, seed=123, pauli_masks=masks) == 0
 
 
+def test_sample_batch_with_pauli_masks_returns_sample_batch() -> None:
+    patch = SurfacePatch.create(distance=3)
+    sampler = build_native_sampler(
+        patch,
+        num_rounds=2,
+        noise=NoiseModel(),
+        basis="Z",
+        twirl=TwirlConfig(),
+    )
+    assert sampler.pauli_frame_lookup is not None
+    zero_masks = np.zeros((4, sampler.num_pauli_sites), dtype=np.int64)
+
+    masks = zero_masks.copy()
+    expected_det_xor: list[list[bool]] | None = None
+    expected_obs_xor: list[list[bool]] | None = None
+    for site in range(sampler.num_pauli_sites):
+        for pauli in (1, 2, 3):
+            candidate = zero_masks.copy()
+            candidate[0, site] = pauli
+            det_xor, obs_xor = sampler.pauli_frame_lookup.compute_mask_xor(candidate)
+            if np.asarray(det_xor).any() or np.asarray(obs_xor).any():
+                masks = candidate
+                expected_det_xor = det_xor
+                expected_obs_xor = obs_xor
+                break
+        if expected_det_xor is not None:
+            break
+
+    assert expected_det_xor is not None
+    assert expected_obs_xor is not None
+
+    zero_batch = sampler.sampler.sample_batch_with_pauli_masks(
+        4,
+        sampler.pauli_frame_lookup,
+        zero_masks,
+        seed=123,
+    )
+    masked_batch = sampler.sampler.sample_batch_with_pauli_masks(
+        4,
+        sampler.pauli_frame_lookup,
+        masks,
+        seed=123,
+    )
+
+    actual_det_xor = np.asarray(zero_batch.detector_events()) ^ np.asarray(masked_batch.detector_events())
+    actual_obs_xor = np.asarray(zero_batch.observable_flips()) ^ np.asarray(masked_batch.observable_flips())
+
+    assert type(masked_batch) is SampleBatch
+    np.testing.assert_array_equal(actual_det_xor, expected_det_xor)
+    np.testing.assert_array_equal(actual_obs_xor, expected_obs_xor)
+
+
 def test_canonical_frame_output_reuses_raw_abstract_sampler_topology() -> None:
     patch = SurfacePatch.create(distance=3)
     raw = build_native_sampler(
diff --git a/python/quantum-pecos/tests/qec/test_dem_equivalence.py b/python/quantum-pecos/tests/qec/test_dem_equivalence.py
index 073470eff..eb937a48c 100644
--- a/python/quantum-pecos/tests/qec/test_dem_equivalence.py
+++ b/python/quantum-pecos/tests/qec/test_dem_equivalence.py
@@ -94,7 +94,9 @@ def test_sample_simple_dem(self) -> None:
         dem_str = "error(0.5) D0"
         dem = ParsedDem.from_string(dem_str)
 
-        det_events, obs_flips = dem.sample_batch(10000, seed=42)
+        batch = dem.sample_batch(10000, seed=42)
+        det_events = batch.detector_events()
+        obs_flips = batch.observable_flips()
 
         assert len(det_events) == 10000
         assert len(obs_flips) == 10000
@@ -109,7 +111,7 @@ def test_sample_decomposed_dem(self) -> None:
         dem_str = "error(0.1) D0 ^ D1"
         dem = ParsedDem.from_string(dem_str)
 
-        det_events, _obs_flips = dem.sample_batch(50000, seed=42)
+        det_events = dem.sample_batch(50000, seed=42).detector_events()
 
         det_array = np.array(det_events)
         # Each sub-mechanism fires independently at p=0.1
@@ -124,11 +126,11 @@ def test_sample_deterministic(self) -> None:
         dem_str = "error(0.5) D0 D1"
         dem = ParsedDem.from_string(dem_str)
 
-        det1, obs1 = dem.sample_batch(1000, seed=123)
-        det2, obs2 = dem.sample_batch(1000, seed=123)
+        batch1 = dem.sample_batch(1000, seed=123)
+        batch2 = dem.sample_batch(1000, seed=123)
 
-        assert det1 == det2
-        assert obs1 == obs2
+        assert batch1.detector_events() == batch2.detector_events()
+        assert batch1.observable_flips() == batch2.observable_flips()
 
 
 class TestAggregation:
@@ -432,8 +434,8 @@ def test_raw_decomposed_syndrome_rates_match(
         num_shots = 100_000
         seed = 42
 
-        raw_dets, _raw_obs = raw_dem.sample_batch(num_shots, seed=seed)
-        decomp_dets, _decomp_obs = decomposed_dem.sample_batch(num_shots, seed=seed)
+        raw_dets = raw_dem.sample_batch(num_shots, seed=seed).detector_events()
+        decomp_dets = decomposed_dem.sample_batch(num_shots, seed=seed).detector_events()
 
         raw_array = np.array(raw_dets)
         decomp_array = np.array(decomp_dets)
@@ -464,8 +466,8 @@ def test_raw_decomposed_per_detector_rates_match(
         num_shots = 100_000
         seed = 42
 
-        raw_dets, _ = raw_dem.sample_batch(num_shots, seed=seed)
-        decomp_dets, _ = decomposed_dem.sample_batch(num_shots, seed=seed)
+        raw_dets = raw_dem.sample_batch(num_shots, seed=seed).detector_events()
+        decomp_dets = decomposed_dem.sample_batch(num_shots, seed=seed).detector_events()
 
         raw_array = np.array(raw_dets)
         decomp_array = np.array(decomp_dets)
@@ -494,8 +496,8 @@ def test_raw_decomposed_logical_rates_match(
         num_shots = 100_000
         seed = 42
 
-        _, raw_obs = raw_dem.sample_batch(num_shots, seed=seed)
-        _, decomp_obs = decomposed_dem.sample_batch(num_shots, seed=seed)
+        raw_obs = raw_dem.sample_batch(num_shots, seed=seed).observable_flips()
+        decomp_obs = decomposed_dem.sample_batch(num_shots, seed=seed).observable_flips()
 
         raw_obs_array = np.array(raw_obs)
         decomp_obs_array = np.array(decomp_obs)
@@ -547,8 +549,8 @@ def test_decomposition_equivalence_various_sizes(
         num_shots = 50_000
         seed = 123
 
-        raw_dets, _ = raw_dem.sample_batch(num_shots, seed=seed)
-        decomp_dets, _ = decomposed_dem.sample_batch(num_shots, seed=seed)
+        raw_dets = raw_dem.sample_batch(num_shots, seed=seed).detector_events()
+        decomp_dets = decomposed_dem.sample_batch(num_shots, seed=seed).detector_events()
 
         raw_array = np.array(raw_dets)
         decomp_array = np.array(decomp_dets)
diff --git a/python/quantum-pecos/tests/qec/test_dem_sampler.py b/python/quantum-pecos/tests/qec/test_dem_sampler.py
index 5728e7e4e..7c0acf1d3 100644
--- a/python/quantum-pecos/tests/qec/test_dem_sampler.py
+++ b/python/quantum-pecos/tests/qec/test_dem_sampler.py
@@ -72,7 +72,9 @@ def test_dem_sampler_sampling() -> None:
     assert len(obs_flips) == 1
 
     # Batch sample
-    det_batch, obs_batch = sampler.sample_batch(100, seed=42)
+    batch = sampler.sample_batch(100, seed=42)
+    det_batch = batch.detector_events()
+    obs_batch = batch.observable_flips()
     assert len(det_batch) == 100
     assert len(obs_batch) == 100
 
@@ -97,11 +99,11 @@ def test_dem_sampler_determinism() -> None:
     sampler = builder.build()
 
     # Same seed should produce same results
-    det1, obs1 = sampler.sample_batch(50, seed=12345)
-    det2, obs2 = sampler.sample_batch(50, seed=12345)
+    batch1 = sampler.sample_batch(50, seed=12345)
+    batch2 = sampler.sample_batch(50, seed=12345)
 
-    assert det1 == det2
-    assert obs1 == obs2
+    assert batch1.detector_events() == batch2.detector_events()
+    assert batch1.observable_flips() == batch2.observable_flips()
 
 
 def test_dem_sampler_statistics() -> None:
diff --git a/python/quantum-pecos/tests/qec/test_dem_sampler_modes.py b/python/quantum-pecos/tests/qec/test_dem_sampler_modes.py
index 24a113e61..89ddaf604 100644
--- a/python/quantum-pecos/tests/qec/test_dem_sampler_modes.py
+++ b/python/quantum-pecos/tests/qec/test_dem_sampler_modes.py
@@ -19,6 +19,7 @@
 import pytest
 from pecos.qec import DagFaultAnalyzer, DemSampler, DemSamplerBuilder
 from pecos_rslib import DagCircuit
+from pecos_rslib.qec import LogicalSubgraphDecoder
 
 if TYPE_CHECKING:
     from pecos.qec import DagFaultInfluenceMap
@@ -73,12 +74,28 @@ def test_raw_sample_returns_correct_shape(self) -> None:
         assert len(outputs) > 0
 
     def test_raw_sample_batch(self) -> None:
-        """Test that sample_batch returns the requested number of shots."""
+        """Raw batches expose measurements but reject decoder operations."""
         dag = _build_repetition_code_circuit(2)
         im = _build_influence_map(dag)
         sampler = DemSampler.raw_uniform(im, 0.01)
-        all_outputs, _all_obs = sampler.sample_batch(100, seed=42)
-        assert len(all_outputs) == 100
+        batch = sampler.sample_batch(100, seed=42)
+        assert len(batch.detector_events()) == 100
+        assert len(batch.get_syndrome(0)) == sampler.num_outputs
+        with pytest.raises(ValueError, match=r"raw-measurement.*measurements, not detector events"):
+            batch.decode_count("error(0.01) D0")
+
+    def test_raw_sample_batch_rejected_by_decoder_objects(self) -> None:
+        """Decoder objects taking a SampleBatch reject raw-measurement batches too."""
+        dag = _build_repetition_code_circuit(2)
+        im = _build_influence_map(dag)
+        sampler = DemSampler.raw_uniform(im, 0.01)
+        batch = sampler.sample_batch(50, seed=42)
+        decoder = LogicalSubgraphDecoder.from_membership(
+            "error(0.1) D0 L0\nerror(0.1) D1 L0",
+            [[0, 1]],
+        )
+        with pytest.raises(ValueError, match=r"raw-measurement.*measurements, not detector events"):
+            decoder.decode_count(batch)
 
     def test_raw_zero_noise_statistics(self) -> None:
         """Test that zero noise produces no syndromes or logical errors."""
diff --git a/python/quantum-pecos/tests/qec/test_dem_sampler_vs_stim.py b/python/quantum-pecos/tests/qec/test_dem_sampler_vs_stim.py
index f938dd96c..4cf43bf10 100644
--- a/python/quantum-pecos/tests/qec/test_dem_sampler_vs_stim.py
+++ b/python/quantum-pecos/tests/qec/test_dem_sampler_vs_stim.py
@@ -329,7 +329,7 @@ def test_detector_firing_rates_correlate(
         seed = 123
 
         # PECOS: get per-detector counts
-        pecos_det_batch, _ = pecos_sampler.sample_batch(num_shots, seed=seed)
+        pecos_det_batch = pecos_sampler.sample_batch(num_shots, seed=seed).detector_events()
         pecos_det_array = np.array(pecos_det_batch)
         pecos_det_rates = pecos_det_array.mean(axis=0)
 
diff --git a/python/quantum-pecos/tests/qec/test_from_guppy_dem.py b/python/quantum-pecos/tests/qec/test_from_guppy_dem.py
index ed4e1b5f9..eda9434dd 100644
--- a/python/quantum-pecos/tests/qec/test_from_guppy_dem.py
+++ b/python/quantum-pecos/tests/qec/test_from_guppy_dem.py
@@ -1149,7 +1149,7 @@ def test_constrained_from_guppy_dem_is_consumable_by_pecos_native_decoder() -> N
     assert sampler.num_observables == dem.num_observables
     assert dem.num_observables == 1  # one logical observable for a single patch
 
-    batch = sampler.generate_samples(16, 0)
+    batch = sampler.sample_batch(16, 0)
     assert batch.num_shots == 16
     # Each shot's syndrome covers exactly the DEM's detectors.
     assert len(batch.get_syndrome(0)) == dem.num_detectors
diff --git a/python/quantum-pecos/tests/qec/test_parsed_dem_sampler.py b/python/quantum-pecos/tests/qec/test_parsed_dem_sampler.py
index 11537c084..f67194946 100644
--- a/python/quantum-pecos/tests/qec/test_parsed_dem_sampler.py
+++ b/python/quantum-pecos/tests/qec/test_parsed_dem_sampler.py
@@ -83,7 +83,7 @@ def test_decomposed_all_fire_together(self) -> None:
         parsed = ParsedDem.from_string(dem_str)
 
         # Sample with naive sampler
-        dets, _ = parsed.sample_batch(10000, seed=42)
+        dets = parsed.sample_batch(10000, seed=42).detector_events()
         dets = np.array(dets)
 
         # D0 and D1 should always fire together
@@ -101,7 +101,7 @@ def test_xor_cancellation(self) -> None:
         parsed = ParsedDem.from_string(dem_str)
 
         # Sample - D0 should never fire due to XOR cancellation
-        dets, _ = parsed.sample_batch(10000, seed=42)
+        dets = parsed.sample_batch(10000, seed=42).detector_events()
         d0_fires = sum(1 for d in dets if d and d[0])
 
         assert d0_fires == 0, "D0 should never fire due to XOR cancellation"
@@ -120,7 +120,7 @@ def test_decomposed_matches_stim_semantics(self) -> None:
 
         # PECOS
         parsed = ParsedDem.from_string(dem_str)
-        pecos_det, _ = parsed.sample_batch(10000, seed=42)
+        pecos_det = parsed.sample_batch(10000, seed=42).detector_events()
         pecos_det = np.array(pecos_det)
 
         # Compare statistics
@@ -186,7 +186,7 @@ def test_optimized_matches_naive_sampler(self) -> None:
         sampler = parsed.to_dem_sampler()
 
         # Naive sampling
-        dets_naive, _ = parsed.sample_batch(50000, seed=42)
+        dets_naive = parsed.sample_batch(50000, seed=42).detector_events()
         naive_rate = sum(1 for d in dets_naive if any(d)) / len(dets_naive)
 
         # Optimized sampling
diff --git a/python/quantum-pecos/tests/qec/test_sample_batch.py b/python/quantum-pecos/tests/qec/test_sample_batch.py
index 4170da0ac..7edc66f0f 100644
--- a/python/quantum-pecos/tests/qec/test_sample_batch.py
+++ b/python/quantum-pecos/tests/qec/test_sample_batch.py
@@ -4,7 +4,7 @@
 """Tests for SampleBatch columnar storage and validation."""
 
 import pytest
-from pecos_rslib.qec import DemSampler, SampleBatch
+from pecos_rslib.qec import DemSampler, ParsedDem, SampleBatch
 
 
 class TestSampleBatchConstruction:
@@ -38,6 +38,21 @@ def test_empty_batch(self):
         batch = SampleBatch([], [])
         assert batch.num_shots == 0
 
+    def test_observable_flips_preserves_explicit_width_for_empty_masks(self):
+        batch = SampleBatch([[1, 0], [0, 1]], [0, 0], num_observables=3)
+        assert batch.observable_flips() == [[False, False, False], [False, False, False]]
+
+    def test_observable_flips_infers_width_from_highest_set_bit(self):
+        batch = SampleBatch([[1, 0], [0, 1]], [0, 1 << 3])
+        assert batch.observable_flips() == [
+            [False, False, False, False],
+            [False, False, False, True],
+        ]
+
+    def test_explicit_observable_width_rejects_out_of_range_mask_bit(self):
+        with pytest.raises(ValueError, match=r"mask bit 3.*num_observables=3"):
+            SampleBatch([[1, 0]], [1 << 3], num_observables=3)
+
 
 class TestGeneratedSampleBatch:
     @pytest.fixture
@@ -63,18 +78,45 @@ def d3_setup(self):
 
     def test_num_shots(self, d3_setup):
         sampler, _ = d3_setup
-        batch = sampler.generate_samples(100, seed=42)
+        batch = sampler.sample_batch(100, seed=42)
+        assert type(batch) is SampleBatch
         assert batch.num_shots == 100
 
+    def test_bulk_accessors_match_per_shot_accessors(self, d3_setup):
+        sampler, _ = d3_setup
+        batch = sampler.sample_batch(73, seed=1234)
+
+        detector_events = batch.detector_events()
+        observable_flips = batch.observable_flips()
+        assert len(detector_events) == batch.num_shots
+        assert len(observable_flips) == batch.num_shots
+        for shot in range(batch.num_shots):
+            assert detector_events[shot] == [bool(value) for value in batch.get_syndrome(shot)]
+            mask = batch.get_observable_mask(shot)
+            assert observable_flips[shot] == [
+                bool(mask & (1 << observable)) for observable in range(sampler.num_observables)
+            ]
+
+    def test_parsed_dem_sample_batch_returns_sample_batch(self):
+        parsed = ParsedDem.from_string("error(0.25) D0 L0")
+        assert type(parsed.sample_batch(5, seed=7)) is SampleBatch
+
+    def test_parsed_dem_row_conversion_preserves_exact_column_positions(self):
+        parsed = ParsedDem.from_string("error(1.0) D0 D2 L0")
+        batch = parsed.sample_batch(4, seed=7)
+
+        assert batch.detector_events() == [[True, False, True]] * 4
+        assert batch.observable_flips() == [[True]] * 4
+
     def test_get_syndrome_shape(self, d3_setup):
         sampler, _ = d3_setup
-        batch = sampler.generate_samples(10, seed=42)
+        batch = sampler.sample_batch(10, seed=42)
         syn = batch.get_syndrome(0)
         assert len(syn) == sampler.num_detectors
 
     def test_get_observable_mask_type(self, d3_setup):
         sampler, _ = d3_setup
-        batch = sampler.generate_samples(10, seed=42)
+        batch = sampler.sample_batch(10, seed=42)
         mask = batch.get_observable_mask(0)
         assert isinstance(mask, int)
 
@@ -89,7 +131,7 @@ def test_decode_count(self, d3_setup):
             stim.Circuit(stim_str).detector_error_model(decompose_errors=True),
         )
 
-        batch = sampler.generate_samples(1000, seed=42)
+        batch = sampler.sample_batch(1000, seed=42)
         errors = batch.decode_count(dem_str, "pymatching")
         assert isinstance(errors, int)
         assert 0 <= errors <= 1000
diff --git a/python/quantum-pecos/tests/qec/test_traced_qis_slow_integration.py b/python/quantum-pecos/tests/qec/test_traced_qis_slow_integration.py
index 3b6b55f55..bdc4c2f0c 100644
--- a/python/quantum-pecos/tests/qec/test_traced_qis_slow_integration.py
+++ b/python/quantum-pecos/tests/qec/test_traced_qis_slow_integration.py
@@ -94,7 +94,7 @@ def _decode_raw_measurements(result, circuit, matching, shots):
 
 def _decode_native_dem_samples(circuit, noise_args, matching, shots, seed):
     sampler = DemSampler.from_circuit(circuit, **noise_args)
-    batch = sampler.generate_samples(shots, seed=seed)
+    batch = sampler.sample_batch(shots, seed=seed)
     syndrome = np.zeros(sampler.num_detectors, dtype=np.uint8)
 
     errors = 0
diff --git a/python/quantum-pecos/tests/qec/test_wide_observables.py b/python/quantum-pecos/tests/qec/test_wide_observables.py
index 80cfaf70e..5cd4bb791 100644
--- a/python/quantum-pecos/tests/qec/test_wide_observables.py
+++ b/python/quantum-pecos/tests/qec/test_wide_observables.py
@@ -84,11 +84,23 @@ def test_decode_count_above_64_observables() -> None:
     n = 65
     dem, membership = _wide_dem(n)
     dec = LogicalSubgraphDecoder.from_membership(dem, membership, "pecos_uf:fast")
-    batch = ParsedDem.from_string(dem).to_dem_sampler().generate_samples(2000, seed=1)
+    batch = ParsedDem.from_string(dem).to_dem_sampler().sample_batch(2000, seed=1)
     count = dec.decode_count(batch)
     assert 0 <= count <= 2000
 
 
+def test_observable_flips_matches_wide_per_shot_masks() -> None:
+    n = 70
+    dem, _ = _wide_dem(n)
+    batch = ParsedDem.from_string(dem).to_dem_sampler().sample_batch(73, seed=17)
+
+    observable_flips = batch.observable_flips()
+    assert all(len(row) == n for row in observable_flips)
+    for shot, row in enumerate(observable_flips):
+        mask = batch.get_observable_mask_wide(shot)
+        assert row == [bool(mask & (1 << observable)) for observable in range(n)]
+
+
 def test_u64_observable_getter_rejects_wide_batch() -> None:
     # get_observable_mask returns a u64 and cannot represent observable >= 64, so
     # it rejects a wide batch; get_observable_mask_wide returns the full Python
@@ -137,7 +149,7 @@ def test_decode_each_returns_python_ints() -> None:
     # the value is not truncated; for the <=64 case it equals the historical u64.
     n = 5
     dem, _ = _wide_dem(n)
-    batch = ParsedDem.from_string(dem).to_dem_sampler().generate_samples(8, seed=1)
+    batch = ParsedDem.from_string(dem).to_dem_sampler().sample_batch(8, seed=1)
     preds = batch.decode_each(dem, "pymatching")
     assert len(preds) == 8
     assert all(isinstance(p, int) for p in preds)
diff --git a/scripts/bench_raw_meas_sampling.py b/scripts/bench_raw_meas_sampling.py
index 7622c1cfd..7692c27d6 100644
--- a/scripts/bench_raw_meas_sampling.py
+++ b/scripts/bench_raw_meas_sampling.py
@@ -3,7 +3,7 @@
 # Licensed under the Apache License, Version 2.0
 """Benchmark: raw measurement sampling / detector DEM vs stabilizer simulation.
 
-Compares detector DEM (generate_samples), raw meas_sampling, and stabilizer.
+Compares detector DEM (sample_batch), raw meas_sampling, and stabilizer.
 Quick smoke test by default (~10s). Use --full for stable headline numbers.
 
 Usage:
@@ -54,7 +54,7 @@ def main():
 
         for shots in shot_list:
             t0 = time.perf_counter()
-            _ = sampler.generate_samples(shots, seed=42)
+            _ = sampler.sample_batch(shots, seed=42)
             t_det = time.perf_counter() - t0
 
             t0 = time.perf_counter()
@@ -94,7 +94,7 @@ def main():
 
         for shots in shot_list:
             t0 = time.perf_counter()
-            batch = sampler.generate_samples(shots, seed=42)
+            batch = sampler.sample_batch(shots, seed=42)
             t_gen = time.perf_counter() - t0
 
             t0 = time.perf_counter()
@@ -112,7 +112,7 @@ def main():
         print()
 
     print("Notes:")
-    print("  generate_samples is 3-6x faster after columnar SampleBatch.")
+    print("  sample_batch is 3-6x faster with columnar SampleBatch storage.")
     print("  End-to-end generate+decode is decode-dominated (<1% generation).")
     if not FULL:
         print("  Use --full for larger shot counts and stable headline numbers.")
diff --git a/scripts/compare_meas_sampling_pipeline.py b/scripts/compare_meas_sampling_pipeline.py
index 5fc987a36..7cac295b3 100644
--- a/scripts/compare_meas_sampling_pipeline.py
+++ b/scripts/compare_meas_sampling_pipeline.py
@@ -119,7 +119,7 @@ def run_native_sampler(tc, noise_args, matching, shots, seed):
     num_dets = sampler.num_detectors
 
     t0 = time.perf_counter()
-    batch = sampler.generate_samples(shots, seed=seed)
+    batch = sampler.sample_batch(shots, seed=seed)
     t_sample = time.perf_counter() - t0
 
     t0 = time.perf_counter()
@@ -215,7 +215,7 @@ def main():
     print("  Circuit: Guppy surface code -> traced QIS -> lower_clifford_rotations()")
     print("  Decoder: PyMatching (stim DEM, decompose_errors=True)")
     print("  meas_sampling: geometric raw measurement DEM sampler + Python extraction")
-    print("  native_sampler: DemSampler.generate_samples (detector events directly)")
+    print("  native_sampler: DemSampler.sample_batch (detector events directly)")
     print("  LER differences from different RNG streams, not systematic bias.")
 
 

From fa35982839e4330589278af36e09223a2a785b84 Mon Sep 17 00:00:00 2001
From: Ciaran Ryan-Anderson <ciaranra@gmail.com>
Date: Sat, 8 Aug 2026 15:58:15 -0600
Subject: [PATCH 2/2] Correct qec stub annotations that contradicted the
 runtime binding signatures

---
 python/pecos-rslib/pecos_rslib.pyi | 31 +++++++++++++++++++-----------
 1 file changed, 20 insertions(+), 11 deletions(-)

diff --git a/python/pecos-rslib/pecos_rslib.pyi b/python/pecos-rslib/pecos_rslib.pyi
index 619a48ae4..23843ae48 100644
--- a/python/pecos-rslib/pecos_rslib.pyi
+++ b/python/pecos-rslib/pecos_rslib.pyi
@@ -2146,7 +2146,7 @@ class qec:
         def to_string_source_graphlike_decomposed(self) -> str: ...
         def to_string_terminal_graphlike_decomposed(self) -> str: ...
         def to_string_graphlike_search_decomposed(self) -> str: ...
-        def all_contribution_effects(self) -> list[Any]: ...
+        def all_contribution_effects(self) -> str: ...
         def contribution_effect_summaries(self) -> list[Any]: ...
         def contribution_render_records(self) -> list[Any]: ...
         def contribution_render_records_with_two_detector_direct_policy(self, policy: str) -> list[Any]: ...
@@ -2154,7 +2154,7 @@ class qec:
         def contribution_render_summaries_with_two_detector_direct_policy(self, policy: str) -> list[Any]: ...
         def contribution_source_graphlike_render_records(self) -> list[Any]: ...
         def contributions_for_effect(self, detectors: Sequence[int], dem_outputs: Sequence[int]) -> list[Any]: ...
-        def contributions_for_mechanism(self, detectors: Sequence[int]) -> list[Any]: ...
+        def contributions_for_mechanism(self, detectors: Sequence[int]) -> str: ...
         def to_sampler(self) -> qec.DemSampler: ...
 
     class DemBuilder:
@@ -2419,12 +2419,12 @@ class qec:
     class LogicalAlgorithmDecoder:
         def __init__(self, descriptor: Mapping[str, Any], inner_decoder: str = ...) -> None: ...
         def feed_dense(self, syndrome: Sequence[int]) -> None: ...
-        def feed_sparse(self, detectors: Sequence[int]) -> None: ...
+        def feed_sparse(self, detectors: Sequence[tuple[int, int]]) -> None: ...
         def flush(self) -> int: ...
         def decode(self, syndrome: Sequence[int]) -> int: ...
         def decode_count(self, batch: qec.SampleBatch) -> int: ...
         def reset(self) -> None: ...
-        def accumulated_obs(self) -> list[bool]: ...
+        def accumulated_obs(self) -> int: ...
         def accumulated_obs_mask(self) -> int: ...
         def rounds_fed(self) -> int: ...
         def num_segments(self) -> int: ...
@@ -2441,9 +2441,12 @@ class qec:
         def decode_count(self, batch: qec.SampleBatch) -> int: ...
         def reset(self) -> None: ...
         def num_segments(self) -> int: ...
+        @property
         def can_window(self) -> bool: ...
-        def actual_num_windows(self) -> int: ...
-        def effective_windowing(self) -> bool: ...
+        @property
+        def actual_num_windows(self) -> list[int]: ...
+        @property
+        def effective_windowing(self) -> str: ...
         def has_decision_points(self) -> bool: ...
         def num_decision_points(self) -> int: ...
         def total_detectors(self) -> int: ...
@@ -2481,7 +2484,7 @@ class qec:
         seed: int = ...,
     ) -> None: ...
     @staticmethod
-    def detector_flip_matrix(fired_per_shot: Sequence[Sequence[int]], num_detectors: int) -> list[list[float]]: ...
+    def detector_flip_matrix(fired_per_shot: Sequence[Sequence[int]], num_detectors: int) -> list[float]: ...
     @staticmethod
     def detector_flip_matrices_by_round(
         fired_per_shot: Sequence[Sequence[int]],
@@ -2493,7 +2496,7 @@ class qec:
         fired_per_shot: Sequence[Sequence[int]],
         num_detectors: int,
         max_order: int = ...,
-    ) -> dict[Any, float]: ...
+    ) -> list[tuple[list[int], float]]: ...
     @staticmethod
     def detector_k_body_rates_by_round(
         fired_per_shot: Sequence[Sequence[int]],
@@ -2502,16 +2505,22 @@ class qec:
         max_order: int = ...,
     ) -> list[Any]: ...
     @staticmethod
-    def compare_flip_matrices_rs(sim: Any, dem: Any, num_detectors: int, min_rate: float = ...) -> dict[str, Any]: ...
+    def compare_flip_matrices_rs(
+        sim: Sequence[float], dem: Sequence[float], num_detectors: int, min_rate: float = ...
+    ) -> tuple[float, float, int, int]: ...
     @staticmethod
-    def compare_k_body_rates_rs(sim: Any, dem: Any, min_rate: float = ...) -> dict[str, Any]: ...
+    def compare_k_body_rates_rs(
+        sim: Sequence[tuple[Sequence[int], float]],
+        dem: Sequence[tuple[Sequence[int], float]],
+        min_rate: float = ...,
+    ) -> list[tuple[int, float, float, list[int]]]: ...
     @staticmethod
     def fit_dem_to_marginals(
         mechanisms: Sequence[tuple[float, Sequence[int], Sequence[int]]],
         target_marginals: Sequence[float],
         max_iterations: int = ...,
         tolerance: float = ...,
-    ) -> tuple[list[float], list[float]]: ...
+    ) -> tuple[list[tuple[float, list[int], list[int]]], list[float]]: ...
     @staticmethod
     def mechanisms_to_dem_string(mechanisms: Sequence[tuple[float, Sequence[int], Sequence[int]]]) -> str: ...
     @staticmethod