Add detector error model generation (Issue #4)

- Add dem.py module for DEM extraction and manipulation
- Extract DEM from circuits with decomposition support
- Save/load DEMs in stim native format
- Convert DEM to JSON for analysis
- Build parity check matrix H for BP decoding
- Integrate DEM generation into CLI with --generate-dem flag
- Add comprehensive test suite for DEM operations
- Add make generate-dem target

Features:
- Extract detector error models from circuits
- Save in .dem format (stim native)
- Export to JSON with structured error information
- Build parity check matrix for BP decoder
- CLI integration for automated workflow

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

Author: ChanceSiyuan

Makefile

@@ -1,10 +1,11 @@
-.PHONY: help install setup test test-cov generate-dataset generate-syndromes clean
+.PHONY: help install setup test test-cov generate-dataset generate-dem generate-syndromes clean
 
 help:
 	@echo "Available targets:"
 	@echo "  install            - Install uv package manager"
 	@echo "  setup              - Set up development environment with uv"
 	@echo "  generate-dataset   - Generate noisy circuit dataset"
+	@echo "  generate-dem       - Generate detector error models"
 	@echo "  generate-syndromes - Generate syndrome database (1000 shots)"
 	@echo "  test               - Run tests"
 	@echo "  test-cov           - Run tests with coverage report"
@@ -22,6 +23,9 @@ setup: install
 generate-dataset:
 	uv run generate-noisy-circuits --distance 3 --p 0.01 --rounds 3 5 7 --task z --output datasets/noisy_circuits
 
+generate-dem:
+	uv run generate-noisy-circuits --distance 3 --p 0.01 --rounds 3 5 7 --task z --output datasets/noisy_circuits --generate-dem
+
 generate-syndromes:
 	uv run generate-noisy-circuits --distance 3 --p 0.01 --rounds 3 5 7 --task z --output datasets/noisy_circuits --generate-syndromes 1000
 

src/bpdecoderplus/cli.py

@@ -15,6 +15,7 @@
     run_smoke_test,
     write_circuit,
 )
+from bpdecoderplus.dem import generate_dem_from_circuit
 from bpdecoderplus.syndrome import generate_syndrome_database_from_circuit
 
 
@@ -70,6 +71,11 @@ def create_parser() -> argparse.ArgumentParser:
         metavar="NUM_SHOTS",
         help="Generate syndrome database with specified number of shots",
     )
+    parser.add_argument(
+        "--generate-dem",
+        action="store_true",
+        help="Generate detector error model (.dem file)",
+    )
     return parser
 
 
@@ -117,6 +123,11 @@ def main(argv: list[str] | None = None) -> int:
         write_circuit(circuit, output_path)
         print(f"Wrote {output_path}")
 
+        # Generate DEM if requested
+        if args.generate_dem:
+            dem_path = generate_dem_from_circuit(output_path)
+            print(f"Wrote {dem_path}")
+
         # Generate syndrome database if requested
         if args.generate_syndromes:
             syndrome_path = generate_syndrome_database_from_circuit(

src/bpdecoderplus/dem.py

@@ -0,0 +1,179 @@
+"""
+Detector Error Model (DEM) extraction module for noisy circuits.
+
+This module provides functions to extract and save Detector Error Models
+from Stim circuits for use in decoder implementations.
+"""
+
+from __future__ import annotations
+
+import json
+import pathlib
+from typing import Any
+
+import numpy as np
+import stim
+
+
+def extract_dem(
+    circuit: stim.Circuit,
+    decompose_errors: bool = True,
+) -> stim.DetectorErrorModel:
+    """
+    Extract Detector Error Model from a circuit.
+
+    Args:
+        circuit: Stim circuit to extract DEM from.
+        decompose_errors: Whether to decompose errors into components.
+
+    Returns:
+        Detector Error Model describing error mechanisms.
+    """
+    return circuit.detector_error_model(decompose_errors=decompose_errors)
+
+
+def save_dem(
+    dem: stim.DetectorErrorModel,
+    output_path: pathlib.Path,
+) -> None:
+    """
+    Save Detector Error Model to file in stim format.
+
+    Args:
+        dem: Detector Error Model to save.
+        output_path: Path to save the DEM (.dem file).
+    """
+    output_path.write_text(str(dem))
+
+
+def load_dem(input_path: pathlib.Path) -> stim.DetectorErrorModel:
+    """
+    Load Detector Error Model from file.
+
+    Args:
+        input_path: Path to the DEM file (.dem).
+
+    Returns:
+        Loaded Detector Error Model.
+    """
+    return stim.DetectorErrorModel.from_file(str(input_path))
+
+
+def dem_to_dict(dem: stim.DetectorErrorModel) -> dict[str, Any]:
+    """
+    Convert DEM to dictionary with structured information.
+
+    Args:
+        dem: Detector Error Model to convert.
+
+    Returns:
+        Dictionary with DEM statistics and error information.
+    """
+    errors = []
+    for inst in dem.flattened():
+        if inst.type == "error":
+            prob = inst.args_copy()[0]
+            targets = inst.targets_copy()
+            detectors = [t.val for t in targets if t.is_relative_detector_id()]
+            observables = [t.val for t in targets if t.is_logical_observable_id()]
+
+            errors.append({
+                "probability": float(prob),
+                "detectors": detectors,
+                "observables": observables,
+            })
+
+    return {
+        "num_detectors": dem.num_detectors,
+        "num_observables": dem.num_observables,
+        "num_errors": len(errors),
+        "errors": errors,
+    }
+
+
+def save_dem_json(
+    dem: stim.DetectorErrorModel,
+    output_path: pathlib.Path,
+) -> None:
+    """
+    Save DEM as JSON for easier analysis.
+
+    Args:
+        dem: Detector Error Model to save.
+        output_path: Path to save the JSON file.
+    """
+    dem_dict = dem_to_dict(dem)
+    with open(output_path, "w") as f:
+        json.dump(dem_dict, f, indent=2)
+
+
+def build_parity_check_matrix(
+    dem: stim.DetectorErrorModel,
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """
+    Build parity check matrix H from DEM for BP decoding.
+
+    Args:
+        dem: Detector Error Model.
+
+    Returns:
+        Tuple of (H, priors, obs_flip) where:
+        - H: Parity check matrix, shape (num_detectors, num_errors)
+        - priors: Prior error probabilities, shape (num_errors,)
+        - obs_flip: Observable flip indicators, shape (num_errors,)
+    """
+    errors = []
+    for inst in dem.flattened():
+        if inst.type == "error":
+            prob = inst.args_copy()[0]
+            targets = inst.targets_copy()
+            detectors = [t.val for t in targets if t.is_relative_detector_id()]
+            observables = [t.val for t in targets if t.is_logical_observable_id()]
+            errors.append({
+                "prob": prob,
+                "detectors": detectors,
+                "observables": observables,
+            })
+
+    n_detectors = dem.num_detectors
+    n_errors = len(errors)
+
+    H = np.zeros((n_detectors, n_errors), dtype=np.uint8)
+    priors = np.zeros(n_errors, dtype=np.float64)
+    obs_flip = np.zeros(n_errors, dtype=np.uint8)
+
+    for j, e in enumerate(errors):
+        priors[j] = e["prob"]
+        for d in e["detectors"]:
+            H[d, j] = 1
+        if e["observables"]:
+            obs_flip[j] = 1
+
+    return H, priors, obs_flip
+
+
+def generate_dem_from_circuit(
+    circuit_path: pathlib.Path,
+    output_path: pathlib.Path | None = None,
+    decompose_errors: bool = True,
+) -> pathlib.Path:
+    """
+    Generate and save DEM from a circuit file.
+
+    Args:
+        circuit_path: Path to the circuit file (.stim).
+        output_path: Optional output path. If None, uses circuit_path with .dem extension.
+        decompose_errors: Whether to decompose errors into components.
+
+    Returns:
+        Path to the saved DEM file.
+    """
+    circuit = stim.Circuit.from_file(str(circuit_path))
+
+    if output_path is None:
+        output_path = circuit_path.with_suffix(".dem")
+
+    dem = extract_dem(circuit, decompose_errors=decompose_errors)
+    save_dem(dem, output_path)
+
+    return output_path