diff --git a/mwpm/README.md b/mwpm/README.md new file mode 100644 index 0000000..1c07ee5 --- /dev/null +++ b/mwpm/README.md @@ -0,0 +1,48 @@ +# MWPM Decoder + +Minimum-weight perfect matching on the surface-code detector graph. +PyMatching does the work; this package is a thin pipeline shim so +`mwpm` and `diffqec` can be swapped through the same +`decode_hardware_results(decoder=...)` entry point. + +## Quick start + +```bash +python -m mwpm.smoke_test +``` + +## Decode + +```python +from mwpm import build_matching_decoder, decode_hardware_results_mwpm + +matcher = build_matching_decoder(noisy_circuit) # build once, amortize +ler, err = decode_hardware_results_mwpm(syndromes, matcher=matcher) +``` + +Or through the pipeline dispatcher: + +```python +from run_on_hardware import decode_hardware_results + +ler, err = decode_hardware_results(syndromes, decoder="mwpm", + stim_circuit_noisy=circuit) +``` + +## On the diffqec mismatch + +`diffqec/` has a class-based API, training data helpers, and a soft +confidence return. MWPM has none of those — the matching graph is +closed-form and PyMatching already exposes the state. Mirroring the +`diffqec` file layout would add empty scaffolding, so `mwpm/` is just +two functions in one module: + +| diffqec | mwpm | why | +|---|---|---| +| `DiffQECDecoder` class | `decode_hardware_results_mwpm()` function | no stateful model | +| `confidence` per shot | not returned | MWPM is hard-decision | +| `data.py` (training data) | — | no training | +| `model.py`, `training.py` | — | closed-form | + +`build_matching_decoder` is the only expensive call; reuse the +returned `Matching` across batches. diff --git a/mwpm/__init__.py b/mwpm/__init__.py new file mode 100644 index 0000000..f9455d0 --- /dev/null +++ b/mwpm/__init__.py @@ -0,0 +1,4 @@ +"""MWPM decoder package.""" +from .integrate import build_matching_decoder, decode_hardware_results_mwpm + +__all__ = ["build_matching_decoder", "decode_hardware_results_mwpm"] diff --git a/mwpm/integrate.py b/mwpm/integrate.py new file mode 100644 index 0000000..93078c0 --- /dev/null +++ b/mwpm/integrate.py @@ -0,0 +1,76 @@ +"""MWPM decoder — minimum-weight perfect matching via PyMatching. + +Thin pipeline glue. The decoder is just two operations: + + 1. build_matching_decoder(circuit) # graph construction + 2. matcher.decode_batch(det_events) # batched matching + +Anything beyond that (data helpers, model classes, confidence scores) +does not exist for MWPM. DiffQEC needs a class because the model is +stateful; PyMatching's ``Matching`` is already the state. +""" +from typing import List, Optional, Tuple +import numpy as np +import pymatching +import stim + +__all__ = ["build_matching_decoder", "decode_hardware_results_mwpm"] + + +def build_matching_decoder( + stim_circuit_noisy: stim.Circuit, + decompose_errors: bool = True, +) -> pymatching.Matching: + """Build a PyMatching MWPM decoder from a noisy Stim circuit's DEM. + + Edge weights are ``w = -ln(p / (1-p))`` from the DEM's per-error + probabilities. Build once per (distance, rounds, noise) and reuse + across batches — this is the expensive step. + """ + dem = stim_circuit_noisy.detector_error_model(decompose_errors=decompose_errors) + return pymatching.Matching.from_detector_error_model(dem) + + +def decode_hardware_results_mwpm( + syndromes: dict, + stim_circuit_noisy: Optional[stim.Circuit] = None, + matcher: Optional[pymatching.Matching] = None, +) -> Tuple[List[float], List[float]]: + """Decode syndromes with MWPM and return per-observable LER + 1-sigma. + + Parameters + ---------- + syndromes : dict from ``extract_syndromes()``. + stim_circuit_noisy : noisy stim.Circuit, used to build the matcher + on the fly when ``matcher`` is None. + matcher : pre-built ``pymatching.Matching``. Preferred for + sweeps so the graph is amortized. + + Returns + ------- + ler : list[float] per-observable logical error rate. + err : list[float] 1-sigma binomial standard error per observable. + """ + if matcher is None: + if stim_circuit_noisy is None: + raise ValueError( + "decode_hardware_results_mwpm requires either `matcher` or " + "`stim_circuit_noisy`." + ) + matcher = build_matching_decoder(stim_circuit_noisy) + + det_events = syndromes["det_events"].astype(np.uint8) + obs_flips = syndromes["obs_flips"].astype(bool) + + predicted = matcher.decode_batch(det_events).astype(bool) + shots = len(obs_flips) + num_obs = obs_flips.shape[1] + + ler: List[float] = [] + err: List[float] = [] + for i in range(num_obs): + wrong = int((predicted[:, i] != obs_flips[:, i]).sum()) + p = wrong / max(shots, 1) + ler.append(float(p)) + err.append(float(np.sqrt(p * (1.0 - p) / max(shots, 1)))) + return ler, err diff --git a/mwpm/smoke_test.py b/mwpm/smoke_test.py new file mode 100644 index 0000000..bf2d7a8 --- /dev/null +++ b/mwpm/smoke_test.py @@ -0,0 +1,74 @@ +"""One-shot smoke test for the MWPM decoder. + +Usage: + python -m mwpm.smoke_test + +Builds a d=3 noisy surface-code circuit, samples detection events, +runs MWPM, prints LER and a sanity-check that we beat random guessing. +""" +import time +import numpy as np +import stim + +from mwpm.integrate import build_matching_decoder, decode_hardware_results_mwpm + + +def main(): + print("=" * 60) + print("MWPM Smoke Test") + print("=" * 60) + + d = 3 + rounds = 3 + shots = 5000 + noise = dict( + after_clifford_depolarization=0.001, + after_reset_flip_probability=0.001, + before_measure_flip_probability=0.001, + before_round_data_depolarization=0.001, + ) + + # 1. Build circuit + sample. + print(f"d={d}, rounds={rounds}, shots={shots}, p=0.001") + t0 = time.time() + circuit = stim.Circuit.generated( + "surface_code:rotated_memory_z", + distance=d, rounds=rounds, **noise, + ) + det_events, obs_flips = circuit.compile_detector_sampler(seed=42).sample( + shots=shots, separate_observables=True, bit_packed=False + ) + det_events = det_events.astype(bool) + obs_flips = obs_flips.astype(bool) + print(f" det={det_events.shape}, obs={obs_flips.shape}, " + f"raw flip rate={obs_flips.mean():.4f} ({time.time()-t0:.2f}s)") + + # 2. Build matcher. + t0 = time.time() + matcher = build_matching_decoder(circuit) + print(f" matcher: {matcher.num_detectors} dets, " + f"{matcher.num_edges} edges ({time.time()-t0:.2f}s)") + + # 3. Decode. + t0 = time.time() + syndromes = {"det_events": det_events, "obs_flips": obs_flips} + ler, err = decode_hardware_results_mwpm(syndromes, matcher=matcher) + print(f" decode: {time.time()-t0:.2f}s") + print(f" LER = {ler[0]:.4f} ± {err[0]:.4f}") + + # 4. Random-guess baseline. + rng = np.random.default_rng(0) + random_wrong = (rng.random(obs_flips.shape) > 0.5) != obs_flips + random_ler = random_wrong.mean() + print(f" random baseline = {random_ler:.4f}") + + # 5. Acceptance: clearly below 0.4 (random is ~0.5). + if ler[0] < 0.4: + print("SMOKE TEST PASSED") + else: + print("SMOKE TEST FAILED") + raise SystemExit(1) + + +if __name__ == "__main__": + main() diff --git a/run_on_hardware.py b/run_on_hardware.py index d64b2e7..6a1b8f4 100644 --- a/run_on_hardware.py +++ b/run_on_hardware.py @@ -89,17 +89,24 @@ def decode_hardware_results( model_path: str | None = None, syndrome_shape: tuple | None = None, L: int | None = None, + stim_circuit_noisy=None, + matcher=None, ) -> tuple[list, list]: """ Decodes the output of extract_syndromes() via your desired decoder. Parameters ---------- - syndromes : dict output of extract_syndromes() - decoder : "mwpm" (default) or "diffqec" - model_path : path to DiffQEC checkpoint (required if decoder="diffqec") - syndrome_shape : (rounds, D) expected by DiffQEC model - L : number of logical observables for DiffQEC + syndromes : dict output of extract_syndromes() + decoder : "mwpm" (default) or "diffqec" + model_path : path to DiffQEC checkpoint (required if decoder="diffqec") + syndrome_shape : (rounds, D) expected by DiffQEC model + L : number of logical observables for DiffQEC + stim_circuit_noisy : stim.Circuit with noise, used to build the MWPM + matching graph on the fly (only used if decoder="mwpm" + and `matcher` is None). + matcher : pre-built pymatching.Matching. Recommended for sweeps + so the graph isn't rebuilt per batch (mwpm only). Returns ------- @@ -115,7 +122,15 @@ def decode_hardware_results( L=L, ) - # Default / placeholder: MWPM or empty - ler, err = [], [] - return ler, err + if decoder == "mwpm": + from mwpm.integrate import decode_hardware_results_mwpm + return decode_hardware_results_mwpm( + syndromes, + stim_circuit_noisy=stim_circuit_noisy, + matcher=matcher, + ) + + raise ValueError( + f"Unknown decoder {decoder!r}. Expected 'mwpm' or 'diffqec'." + )