Skip to content

qsimlab.qec

qsimlab.qec

Quantum error correction: memory circuits, detector sampling, error models, decoders.

Provisional (phase 2, see python/API.md §1 and §8).

The pipeline is build → sample → decode → count:

>>> import qsimlab.qec as qec
>>> c = qec.surface_code_memory(3, rounds=3, p=1e-3)
>>> len(c.detectors), len(c.observables)
(24, 1)
>>> dets, obs = qec.sample_detectors(c, 2048, seed=1)
>>> dets.shape, dets.dtype, obs.shape
((2048, 24), dtype('bool'), (2048, 1))
>>> dem = qec.detector_error_model(c)
>>> dem.num_detectors, dem.num_observables
(24, 1)
>>> qec.circuit_distance(c).distance
3
>>> r = qec.logical_error_rate(c, 20_000, seed=2)
>>> r.errors < 100 and r.shots == 20_000
True

Conventions (in addition to python/API.md §5):

  • Detection events are reported relative to the noiseless reference (Stim's convention): a noiseless shot reads all zeros.
  • Bit packing (packed=True) is Stim's bit_packed layout: uint8[shots, ceil(n / 8)], entry k is bit k % 8 of byte k // 8.
  • Seeds: shots are drawn in chunks of 1024 from independent streams derived from (seed, chunk), so a seed gives the same samples for any thread count. seed=None draws a fresh seed from the OS.
  • Noise models (noise=, strength p): "cnot" (DEPOLARIZE2(p) after every CNOT, nothing else), "uniform" (DEPOLARIZE2(p) after CNOTs, DEPOLARIZE1(p) on idle qubits in every moment, readout flip p, reset flip p) and "si1000" (Gidney et al.'s superconducting-inspired model: CNOT p, idle p/10 during gate layers and 2p during measure/reset layers, readout flip 5p, reset flip 2p). Every noise location is an explicit op of the returned circuit (the readout flip is circuit.readout_error), so circuit.to_stim() is exactly the sampled circuit.

CodeLayout dataclass

CodeLayout(
    code: str,
    distance: int,
    rounds: int,
    basis: str,
    noise: str,
    p: float,
    data_qubits: Tuple[int, ...],
    ancilla_qubits: Tuple[int, ...],
    flag_qubits: Tuple[int, ...],
    qubit_coords: Tuple[Tuple[float, float], ...],
    detector_coords: ndarray,
    detector_basis: Tuple[str, ...],
    flag_detectors: ndarray,
    schedule: Optional[Tuple[Tuple[int, ...], ...]] = None,
)

Where everything sits in a generated memory circuit.

  • data_qubits, ancilla_qubits, flag_qubits: qubit indices;
  • qubit_coords: (x, y) per qubit;
  • detector_coords: float64[D, 3] with (x, y, round) per detector;
  • detector_basis: "X" or "Z" per detector (the stabilizer type it checks);
  • flag_detectors: bool[D], true for flag-qubit detectors (colour code);
  • memory_detectors: indices of the non-flag detectors of the memory basis, the sector the logical observable lives in (pass it to circuit_distance as detectors= for a much faster search).

DetectorSampler

DetectorSampler(circuit: Circuit, engine: str = 'auto')

A compiled detection-event sampler for a noisy Clifford circuit.

engine: "fast" (the Poisson-hit FastSampler of research/qec/fast-sampler.md), "symphase" (the plain SymPhase sampler: same distribution, slower) or "auto" (FastSampler, falling back to SymPhase with a note if it rejects the circuit). Compile once, sample many times:

>>> s = DetectorSampler(repetition_code_memory(3, 2, p=0.05))
>>> s.engine, s.num_detectors, s.num_observables
('fast', 6, 1)
>>> a, _ = s.sample(100, seed=5); b, _ = s.sample(100, seed=5, threads=1)
>>> bool((a == b).all())
True

engine property

engine: str

The sampler that runs: "fast" or "symphase".

note property

note: str

Why the engine differs from the request (empty if it does not).

compile_time property

compile_time: float

Seconds spent compiling (SymPhase frame + hit tables).

sample

sample(
    shots: int,
    *,
    seed: Optional[int] = None,
    packed: bool = False,
    transposed: bool = False,
    threads: Optional[int] = None,
) -> Tuple[ndarray, ndarray]

(detectors, observables): bool[shots, D] and bool[shots, O], or Stim-style bit-packed uint8[shots, ceil(D/8)] arrays with packed=True.

transposed=True returns detector-major bit-packed arrays uint8[D, ceil(shots/8)] (bit s % 8 of byte s // 8 of row r is detector r in shot s; padding bits are 0). It skips the bit transposition, so it is the fastest layout (the Rust bench's throughput), and it is what per-detector statistics want. The same seed gives the same shots in every layout.

DemError

Bases: NamedTuple

One error mechanism: independent event with probability that flips the detectors and observables (sorted index tuples).

DetectorErrorModel

DetectorErrorModel(
    num_detectors: int,
    num_observables: int,
    errors: Iterable[
        Union[
            DemError,
            Tuple[float, Sequence[int], Sequence[int]],
        ]
    ],
)

A detector error model: independent mechanisms over detectors/observables.

Built from a circuit by detector_error_model (exact: each depolarizing channel is converted into independent components the way Stim does, mechanisms with equal signatures merged), or parsed from Stim's DEM text format with from_stim_dem.

>>> dem = DetectorErrorModel.from_stim_dem('''
...     error(0.1) D0 L0
...     error(0.2) D0 D1
...     detector D2
... ''')
>>> dem.num_detectors, dem.num_observables, len(dem)
(3, 1, 2)
>>> print(dem.to_stim_dem())
error(0.1) D0 L0
error(0.2) D0 D1
detector D2

merged

merged() -> Dict[
    Tuple[Tuple[int, ...], Tuple[int, ...]], float
]

{(detectors, observables): probability} with equal signatures merged (as independent events).

approx_equal

approx_equal(
    other: "DetectorErrorModel",
    *,
    rtol: float = 1e-09,
    atol: float = 0.0,
) -> bool

Same detectors/observables and the same merged mechanisms with probabilities equal to rtol/atol.

matrices

matrices() -> Tuple[ndarray, ndarray, ndarray]

Dense (H, L, p): uint8[D, E] check matrix, uint8[O, E] observable matrix and float64[E] probabilities.

to_stim_dem

to_stim_dem() -> str

Stim's DEM text (error(p) D.. L.. lines; detector / logical_observable declarations keep the counts).

from_stim_dem classmethod

from_stim_dem(text: str) -> 'DetectorErrorModel'

Parses Stim's DEM text: error, detector, logical_observable, shift_detectors, repeat blocks and ^ decomposition separators (components are XOR-ed back together). detector_separator and coordinates are ignored.

to_stim

to_stim() -> Any

A stim.DetectorErrorModel (needs stim).

from_stim classmethod

from_stim(dem: Any) -> 'DetectorErrorModel'

From a stim.DetectorErrorModel (via its text, loops flattened).

graphlike

graphlike() -> Tuple['DetectorErrorModel', int]

Decomposes every mechanism with more than two detectors into mechanisms of the model with one or two detectors whose XOR (detectors and observables) reproduces it, as Stim's decompose_errors does; returns (graphlike model, number of mechanisms that could not be decomposed and were dropped). Used by the matching decoders.

DistanceResult dataclass

DistanceResult(
    distance: Optional[int],
    count: int,
    count_capped: bool,
    example: List[DemError],
    complete: bool,
    certified: bool,
    lower_bound: int,
    searched_weight: int,
    timed_out: bool,
    node_limit_hit: bool,
    nodes: int,
    wall_time: float,
)

Result of circuit_distance.

  • distance: minimum number of error mechanisms that flip the observable without firing a detector (None: none found up to searched_weight, or the search stopped, see complete);
  • count: number of distinct minimum-weight logicals (count_capped: the count reached the cap);
  • example: one minimum-weight logical as DemError mechanisms;
  • complete: the answer is exact (no timeout, node limit or sector caveat);
  • certified: with detectors= (a sector): the example lifts to the full model, so distance is the full circuit distance (a sector search alone gives a lower bound);
  • lower_bound: a proven lower bound on the distance;
  • nodes, wall_time: search effort.

Decoder

Base class: decode(dets) -> bool[shots, num_observables] predicted flips.

BpOsdDecoder

BpOsdDecoder(
    dem: DetectorErrorModel,
    *,
    max_iter: int = 50,
    ms_scale: float = 0.625,
    osd_order: int = 10,
)

Bases: Decoder

Belief propagation + ordered-statistics decoding (qsim_lab::qec::bposd: normalised min-sum BP, OSD-CS of order osd_order) on the full hypergraph DEM; parallel over shots with the GIL released. At most 64 observables.

PyMatchingDecoder

PyMatchingDecoder(dem: DetectorErrorModel, **options: Any)

Bases: Decoder

Minimum-weight perfect matching via PyMatching (optional dependency).

Mechanisms with more than two detectors are decomposed into graphlike ones (DetectorErrorModel.graphlike); dropped counts those that could not be.

TesseractDecoder

TesseractDecoder(dem: DetectorErrorModel, **options: Any)

Bases: Decoder

Tesseract (A* search, near-optimal; optional tesseract-decoder + stim).

LogicalErrorRate dataclass

LogicalErrorRate(
    engine: str,
    components: List[Tuple[int, int, str]],
    seed: int,
    precision: str,
    wall_time: float,
    explanation: Optional[Explanation] = None,
    errors: int = 0,
    shots: int = 0,
    rate: float = 0.0,
    ci: Tuple[float, float] = (0.0, 1.0),
    decoder: str = "",
    rounds: Optional[int] = None,
    per_round: Optional[float] = None,
    per_round_ci: Optional[Tuple[float, float]] = None,
    stats: Dict[str, Any] = dict(),
)

Bases: Result

Result of logical_error_rate. A shot fails if any observable is mispredicted.

  • errors, shots, rate; ci: Wilson 95% interval;
  • rounds: if known, per_round / per_round_ci convert with p_L = (1 - (1 - 2 eps)^rounds) / 2;
  • decoder; stats: decoder statistics (BP convergence, OSD calls).

surface_code_memory

surface_code_memory(
    d: int,
    rounds: Optional[int] = None,
    basis: str = "Z",
    *,
    p: float = 0.0,
    noise: str = "uniform",
    return_layout: bool = False,
) -> Union[Circuit, Tuple[Circuit, CodeLayout]]

Rotated surface code memory experiment.

d*d data qubits at (2x+1, 2y+1) and d*d - 1 auxiliaries (Stim's surface_code:rotated_memory_* layout and hook-safe CNOT order). Data start in the memory basis, every round resets the auxiliaries (R / RX), runs four CNOT layers and measures them (M / MX); the data are measured in the memory basis at the end.

Detectors, in order: per round, first the memory-basis stabilizers (round 0: the raw outcome; later: XOR with the previous round), then (from round 1) the other basis; finally the memory-basis stabilizers from the data. One observable: Z_L on the y = 1 row (Z memory) or X_L on the x = 1 column (X memory).

>>> c = surface_code_memory(3, 2, "X", p=1e-3, noise="si1000")
>>> c.num_qubits, len(c.detectors), c.readout_error
(17, 16, 0.005)

repetition_code_memory

repetition_code_memory(
    d: int,
    rounds: Optional[int] = None,
    *,
    p: float = 0.0,
    noise: str = "uniform",
    return_layout: bool = False,
) -> Union[Circuit, Tuple[Circuit, CodeLayout]]

Bit-flip repetition code memory: data 0..d-1, auxiliary d+j measures Z_j Z_{j+1} every round; detectors as in surface_code_memory ((d-1) * (rounds+1) of them), observable: the last data qubit.

>>> c = repetition_code_memory(5, 3, p=0.01)
>>> c.num_qubits, len(c.detectors)
(9, 16)

color_code_schedule

color_code_schedule(
    d: int, schedule: Any = "kf"
) -> Tuple[List[List[int]], List[bool]]

Resolves a colour-code schedule spec to (rows, flag_mask).

schedule: "kf" (Kishony–Fowler's colour-dependent schedule), "tri" (Lee et al.'s uniform tri-optimal), "global" (the exact global-search schedules of research/qec/colour-global.md: d = 9, circuit distance 8, and d = 11 with 7+7 layers, circuit distance 10), a path to a .sched file (one line per plaquette: steps for positions a..f, 0 for absent ones, optional trailing F = flag it), or a (num_plaquettes, 6) array of steps.

color_code_memory

color_code_memory(
    d: int,
    rounds: Optional[int] = None,
    basis: str = "Z",
    *,
    schedule: Any = "kf",
    flags: Union[
        bool, str, Sequence[int], Sequence[bool]
    ] = False,
    p: float = 0.0,
    noise: str = "cnot",
    return_layout: bool = False,
) -> Union[Circuit, Tuple[Circuit, CodeLayout]]

Triangular 6.6.6 colour code memory experiment (one auxiliary per plaquette).

The round structure is Kishony–Fowler's (arXiv:2603.28852): CX data→aux over the schedule's layers, M, RX, CX aux→data over the same layers, MX, R; built by the engine's qsim_lab::qec::color (the circuits of research/qec/qec-r4.md and research/qec/colour-global.md). schedule: see color_code_schedule. flags: False, True / "boundary" (a flag qubit on every boundary-touching plaquette, the hook-free boundary of research/qec/colour-flags.md), or plaquette indices / a boolean mask. Default noise is "cnot" (K–F's headline model).

>>> c, lay = color_code_memory(5, 2, p=1e-3, return_layout=True)
>>> c.num_qubits, len(c.detectors), len(lay.memory_detectors)
(28, 36, 27)

sample_detectors

sample_detectors(
    circuit: Circuit,
    shots: int,
    *,
    seed: Optional[int] = None,
    engine: str = "auto",
    packed: bool = False,
    transposed: bool = False,
    threads: Optional[int] = None,
) -> Tuple[ndarray, ndarray]

Samples detection events and observable flips: (dets, obs).

One-shot form of DetectorSampler (which compiles once). Arrays are bool[shots, D], bool[shots, O] (packed=True: Stim's bit-packed uint8; transposed=True: detector-major bit-packed, the fastest layout, see DetectorSampler.sample). Raises UnsupportedOperationError for non-Clifford gates or classically controlled operations.

>>> c = repetition_code_memory(3, 1)            # noiseless: all zeros
>>> dets, obs = sample_detectors(c, 10, seed=0)
>>> int(dets.sum()), int(obs.sum())
(0, 0)

detector_error_model

detector_error_model(
    circuit: Circuit,
) -> DetectorErrorModel

The circuit-derived detector error model (no decomposition).

Every Pauli noise channel of the circuit becomes independent mechanisms (a flip: one; DEPOLARIZE1/DEPOLARIZE2: 3/15 components with Stim's exact disjoint→independent conversion), propagated to the detectors and observables; mechanisms with equal signatures are merged. Equals Stim's circuit.detector_error_model() of circuit.to_stim() up to floating-point rounding. Raises UnsupportedOperationError if a detector or observable is not deterministic.

>>> dem = detector_error_model(repetition_code_memory(3, 1, p=0.1, noise="cnot"))
>>> dem
DetectorErrorModel(detectors=4, observables=1, errors=9)

circuit_distance

circuit_distance(
    circuit: Union[Circuit, DetectorErrorModel],
    *,
    max_weight: Optional[int] = None,
    observable: int = 0,
    detectors: Optional[Sequence[int]] = None,
    count_cap: int = 1000000,
    node_limit: Optional[int] = None,
    timeout: Optional[float] = None,
) -> DistanceResult

Exact circuit-level distance by branch and bound (qsim_lab::qec::distance).

Searches the detector error model for the smallest set of mechanisms that flips observable and no detector, and counts all such sets. detectors: restrict to a sector (e.g. layout.memory_detectors): much faster on CSS memories, a lower bound in general, exact when certified. timeout (seconds) switches to iterative deepening over the weight; each weight gets the node budget that the measured search rate allows in the remaining time, so the limit is approximate (typically within a few tens of percent). node_limit bounds the search nodes. On a stop, lower_bound still reports what was proven.

>>> r = circuit_distance(surface_code_memory(5, 2, p=1e-3))
>>> r.distance, r.complete
(5, True)

make_decoder

make_decoder(
    dem: DetectorErrorModel,
    decoder: str = "bposd",
    **options: Any,
) -> Decoder

A decoder for dem: "bposd" (built in), "pymatching" or "tesseract" (optional packages; MissingDependencyError if absent). options go to the decoder's constructor.

decode

decode(
    dem: DetectorErrorModel,
    dets: ndarray,
    decoder: Union[str, Decoder] = "bposd",
    *,
    packed: Optional[bool] = None,
    threads: Optional[int] = None,
    **options: Any,
) -> ndarray

Predicted observable flips bool[shots, O] for detection events dets (bool[shots, D] or bit-packed uint8).

>>> c = repetition_code_memory(5, 3, p=0.02)
>>> dets, obs = sample_detectors(c, 1000, seed=3)
>>> pred = decode(detector_error_model(c), dets)
>>> float((pred != obs).any(axis=1).mean()) < 0.05
True

wilson_interval

wilson_interval(
    errors: int, shots: int, z: float = 1.959963984540054
) -> Tuple[float, float]

Wilson score interval (default 95%) for errors / shots.

logical_error_rate

logical_error_rate(
    circuit: Circuit,
    shots: int,
    *,
    decoder: Union[str, Decoder] = "bposd",
    seed: Optional[int] = None,
    max_errors: Optional[int] = None,
    rounds: Optional[int] = None,
    engine: str = "auto",
    dem: Optional[DetectorErrorModel] = None,
    threads: Optional[int] = None,
    **decoder_options: Any,
) -> LogicalErrorRate

Monte-Carlo logical error rate of a memory circuit with a decoder.

Samples shots (stopping early once max_errors failures are seen; checked every 65,536 shots so results stay seed-reproducible), decodes with the circuit's own detector error model (or dem) and counts shots where any observable is mispredicted. decoder="bposd" runs sampling and decoding fused in native code; other decoders run in Python chunks.

>>> c = surface_code_memory(3, 3, p=2e-3)
>>> r = logical_error_rate(c, 10_000, seed=7, rounds=3)
>>> 0 <= r.rate < 0.05, r.ci[0] <= r.rate <= r.ci[1]
(True, True)