Consolidate DEM batch sampling onto SampleBatch with bulk accessors

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)`,

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,

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] = []

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)

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")
 

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)

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)

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

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(

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:

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)

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 = []

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)