Skip to content

QEC (stabilizer decoding: erasure-aware, MWPM, and blind brute-force)

A stabilizer code protects information by encoding it across several physical qubits such that certain "stabilizer" measurements always return the same value on an error-free state -- a syndrome. When an error flips one of those measurements, the syndrome pattern reveals (up to the code's own limits) which error happened, without ever measuring the encoded data itself. This module is generic and code-agnostic: syndrome computation and three decoders work for any stabilizer code you hand them, not just one built in.

Step 1. Do two Pauli operators commute?

from dense_evolution.physics.qec import pauli_commutes

pauli_commutes('XII', 'ZZI')

False

pauli_commutes(p1, p2) checks whether two full-length Pauli strings commute -- False here because X and Z overlap (anti-commute) on qubit 0. This is the building block a stabilizer measurement's outcome is built from: an error anticommuting with a stabilizer flips that stabilizer's measured value, which is exactly the syndrome bit compute_syndrome (Step 2) reads off.

Step 2. The syndrome of a real error

from dense_evolution.physics.qec import compute_syndrome

x_stabilizers = ['IIIXXXX', 'IXXIIXX', 'XIXIXIX']
z_stabilizers = ['IIIZZZZ', 'IZZIIZZ', 'ZIZIZIZ']
stabilizers = x_stabilizers + z_stabilizers

error = 'IIIZIII'
compute_syndrome(error, stabilizers)

(1, 0, 0, 0, 0, 0)

stabilizers here are the Steane [[7,1,3]] code's real 6 generators (3 checking X errors, 3 checking Z errors) -- the running example for the rest of this page. error is a Z flip on qubit 3 alone. compute_syndrome returns one bit per stabilizer: 1 where that stabilizer anticommutes with the error. Only the first Z-type stabilizer (IIIZZZZ, which includes qubit 3) flips -- the other five commute with a lone Z error.

Step 3. Recovering the error, blind

from dense_evolution.physics.qec import blind_minimum_weight_decode

syndrome = compute_syndrome(error, stabilizers)
decoded = blind_minimum_weight_decode(syndrome, n_qubits=7, stabilizers=stabilizers)
decoded == error

True

blind_minimum_weight_decode never sees error -- only syndrome -- and searches every possible Pauli error in increasing weight order, returning the minimum-weight one that reproduces it (or None if more than one error at that weight matches, never a guess). Brute force, O(3**w) at weight w, so it's for small codes and low weights -- but it works on any stabilizer code, including ones (like Steane) the matching-graph decoder below structurally cannot handle.

Step 4. Recovering the error, with a known location

from dense_evolution.physics.qec import erasure_aware_decode

erasure_aware_decode(syndrome, heralded_qubits=[3], n_qubits=7, stabilizers=stabilizers)

'IIIZIII'

If something upstream already knows where an error might have happened (a heralded lost photon in a dual-rail photonic qubit, say) -- not just the syndrome -- erasure_aware_decode uses that location directly instead of searching blind. A distance-3 code like Steane can resolve up to 2 heralded erasures, versus only 1 unlocated error blind -- knowing the location is strictly more powerful, the real result behind quantum erasure-channel codes (Grassl, Beth & Pellizzari, Phys. Rev. A 56, 33 (1997)).

Step 5. pymatching's real limitation

from dense_evolution.physics.qec import pymatching_decode

error_x = 'IIIXIII'
syndrome_z = compute_syndrome(error_x, z_stabilizers)
pymatching_decode(syndrome_z, z_stabilizers, n_qubits=7, error_type='X')

ValueError: pymatching's matching-graph decoder needs every qubit checked by AT MOST 2 stabilizers (a graph edge connects at most 2 detector nodes) -- qubit(s) [6] are each checked by [3] stabilizers here. ...

pymatching_decode (an optional dependency, pip install dense-evolution[pymatching]) is much faster than blind brute-force decoding -- but only for "graph-like" codes, where every qubit sits on at most 2 stabilizers (true for the surface code). Steane's weight-4 stabilizers check some qubits 3 times, so pymatching_decode raises this ValueError rather than silently returning a wrong answer -- blind_minimum_weight_decode (Step 3) or erasure_aware_decode (Step 4) are the ones that work here instead.

Step 6. Is a noise process bursty, or Poissonian?

import numpy as np
from dense_evolution.physics.qec import counts_in_intervals_dimension

rng = np.random.default_rng(0)
window_sizes = [1, 2, 4, 8, 16, 32]

poisson_times = np.sort(rng.uniform(0, 100, 200))
dim, r2, _ = counts_in_intervals_dimension(poisson_times, window_sizes)
dim, r2

(1.011649262707353, 0.9998031727847672)

A separate question from decoding itself: is a stream of error/erasure timestamps temporally clustered (bursty, e.g. a cosmic-ray impact) or Poissonian (uniformly random)? For a homogeneous Poisson process, the mean count of nearby events within radius r scales as r^1 exactly -- dim above lands almost exactly on 1.0, with r2 close to 1 confirming the fit itself is trustworthy. A genuinely bursty stream (150 events packed into [0, 10], 50 more into [50, 52], otherwise identical in count) gives a real, depressed exponent instead:

burst_times = np.sort(np.concatenate([rng.uniform(0, 10, 150), rng.uniform(50, 52, 50)]))
counts_in_intervals_dimension(burst_times, window_sizes)[:2]

(0.7347741100131827, 0.9721804667576291)

0.73 instead of 1.0 -- exactly the signature real burst-like error sources (cosmic rays, correlated hardware glitches) produce. Always check r2 before trusting the dimension itself (rule of thumb: below ~0.98 means don't) -- a narrow or poorly-covered range of window sizes can fit a meaningless slope.

Details

decode_with_erasure_fallback composes Steps 3-4 into the real-world decoding policy, not another raw decoder: use erasure_aware_decode when there are heralded qubits and it resolves the syndrome uniquely, otherwise fall back to blind_minimum_weight_decode. Never worse than always calling blind decoding directly -- promoted from Dense-Evolution-Discovery's cosmic-ray-burst-as-erasure experiment, where this exact fallback logic was first written inline in a Monte Carlo loop testing whether knowing which qubits a real cosmic-ray burst hit (arXiv:2104.05219) lets a Steane code recover better than blind decoding alone.

The naive shortcut that doesn't work: calling erasure_aware_decode with every qubit passed as heralded, instead of using blind_minimum_weight_decode, fails in practice -- with no qubit assumed error-free, many stabilizer-equivalent errors share the same syndrome, so erasure_aware_decode's "exactly one match" criterion is essentially always violated. Minimum-weight selection (Step 3) is what makes blind decoding well-posed at all.

Real validation, not just unit tests: a Steane-specific version of this decoder was checked against STIM's native HERALDED_ERASE noise channel -- 0 decoding failures across every double-erasure shot tested (>60,000 shots total, 40,000 trials x 10 physical error rates), versus a real ~25% failure rate for a standard syndrome-only decoder blind to the erasure locations. See Dense-Evolution-Discovery's Steane [[7,1,3]] investigation and Gu, Vaknin, Retzker & Kubica, "Optimizing quantum error correction protocols with erasure qubits," PRX Quantum 6, 040354 (2025), arXiv:2408.00829.

Never guesses: every decoder here returns None on an ambiguous or unresolvable syndrome rather than a best-effort guess -- more heralded qubits than the code can actually resolve, or a blind search with more than one minimum-weight match, both yield None.

qec

Stabilizer-code quantum error correction utilities: generic Pauli-string commutation/syndrome primitives, an erasure-aware decoder, and a minimum-weight-perfect-matching (MWPM) decoder.

Promoted from Dense-Evolution-Discovery's Steane [[7,1,3]] code investigation (scripts/steane_code_block6_erasure_conversion.py), where a Steane-specific version of erasure_aware_decode was first built and verified: on shots with exactly 2 simultaneous heralded erasures, it achieved exactly 0 decoding failures across 94-12,469 such shots at every tested physical error rate (0 failures out of >60,000 double-erasure shots total, 40,000 trials x 10 p-values), versus ~25% failure for a standard syndrome-only decoder blind to the erasure locations -- a clean confirmation of the real erasure-correction bound below. This version is code-agnostic (works from any stabilizer generator list, not a hand-built Steane-specific table), so it moved here instead of staying Discovery-repo-specific research code.

Erasure-aware decoding exploits a real, foundational fact: Grassl, Beth & Pellizzari, "Codes for the quantum erasure channel", Phys. Rev. A 56, 33 (1997) -- a distance-d stabilizer code can correct up to (d-1) ERASURES (known-location errors, e.g. a heralded lost photon in a dual-rail photonic qubit), versus only floor((d-1)/2) arbitrary (unlocated) errors. Erasure location information is worth roughly twice as much as an ordinary syndrome bit, because knowing WHERE the error is removes exactly the ambiguity a blind syndrome-only decoder has to guess at.

pymatching_decode (prog.txt Sezione 4.3) fills the complementary, far more common case: a real decoder for when NO erasure locations are known at all (the standard setting for e.g. a surface code under generic physical noise) -- erasure_aware_decode's brute force (4**k over heralded qubits) has no answer when k=0 heralded qubits, by design (returns None). Backed by pymatching (Apache-2.0, oscarhiggott/PyMatching, the standard minimum-weight-perfect-matching decoder for stabilizer codes), an optional dependency (pip install dense-evolution[pymatching]), not a required one -- most callers of this module (erasure-aware decoding, raw syndrome computation) never need it. Only works for "graph-like" codes (see pymatching_decode's own docstring for the real >2-checks-per-qubit constraint) -- NOT Steane and similar small non-topological codes.

blind_minimum_weight_decode (prog.txt Sezione 4.4) covers what neither of the above two can: blind (no erasure locations) decoding for codes pymatching_decode structurally can't handle, like Steane. Pure Python, no new dependency -- brute-forces every possible error in increasing WEIGHT order (0, 1, 2, ...), same itertools.product machinery erasure_aware_decode uses, just searching over which qubits are non-identity too instead of only over Pauli letters on a fixed known set. Tried first: reusing erasure_aware_decode itself with every qubit passed as "heralded" (heralded_qubits=range(n_qubits)) -- this does NOT work (verified directly on Steane: returns None even for a plain single-qubit error), because with no qubits assumed error-free, many stabilizer-equivalent full-length errors share the same syndrome, so erasure_aware_decode's "exactly one match total" criterion is essentially always violated. Minimum-weight selection (return the lowest-weight match, not require a totally unique one across every possible weight) is what makes blind decoding well-posed at all -- the same principle pymatching_decode's MWPM implements via a matching graph instead of brute force.

counts_in_intervals_dimension answers a question upstream of all of the above: IS a given stream of error/erasure timestamps actually bursty, or Poissonian? It generalizes the "counts-in-spheres" fractal dimension estimator used to measure large-scale cosmic homogeneity (Scrimgeour et al. 2012, MNRAS 425, 116, arXiv:1205.6812) from 3-D space to 1-D time: D~=1 means homogeneous arrivals, D<1 means temporal clustering (e.g. cosmic-ray-correlated bursts, arXiv:2104.05219 -- the same physical noise source cosmic_ray_burst_profile in dense_evolution.mitigation.zne models). Returns the fit's R^2 alongside D deliberately, never just the number -- a box-counting fit through too few points spanning too narrow a range of scales can report a fractal dimension as absurd as 142 in 3 dimensions purely from fit noise, not from any real structure in the data.

pauli_commutes 

pauli_commutes(p1: str, p2: str) -> bool

Whether two equal-length Pauli strings (each character in IXYZ, no global phase) commute -- the standard symplectic rule: they commute iff the number of qubit positions where the local single-qubit Paulis anticommute (X/Z, X/Y, or Y/Z, in either order; I commutes with everything) is EVEN.

pauli_commutes('XX', 'ZZ') # X,Z anticommute at both qubits -> 2 (even) -> commute True pauli_commutes('XI', 'ZI') # X,Z anticommute at 1 qubit -> 1 (odd) -> anticommute False

Source code in dense_evolution/physics/qec.py 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
def pauli_commutes(p1: str, p2: str) -> bool:
    """Whether two equal-length Pauli strings (each character in IXYZ, no
    global phase) commute -- the standard symplectic rule: they commute
    iff the number of qubit positions where the local single-qubit
    Paulis anticommute (X/Z, X/Y, or Y/Z, in either order; I commutes
    with everything) is EVEN.

    >>> pauli_commutes('XX', 'ZZ')   # X,Z anticommute at both qubits -> 2 (even) -> commute
    True
    >>> pauli_commutes('XI', 'ZI')   # X,Z anticommute at 1 qubit -> 1 (odd) -> anticommute
    False
    """
    if len(p1) != len(p2):
        raise ValueError(f"Pauli strings must be equal length: {len(p1)} != {len(p2)}")
    n_anticommuting = sum(
        1 for a, b in zip(p1, p2)
        if a != 'I' and b != 'I' and a != b and frozenset((a, b)) in _ANTICOMMUTING_PAIRS
    )
    return n_anticommuting % 2 == 0

compute_syndrome 

compute_syndrome(
    pauli_error: str, stabilizers: Sequence[str]
) -> tuple

The syndrome (one bit per stabilizer generator, 1 = anticommutes / detected, 0 = commutes / undetected) a given Pauli error string would produce against stabilizers (a list of equal-length Pauli strings, the code's stabilizer generators -- X-type, Z-type, or mixed; this function doesn't assume a CSS structure).

Source code in dense_evolution/physics/qec.py 
120
121
122
123
124
125
126
def compute_syndrome(pauli_error: str, stabilizers: Sequence[str]) -> tuple:
    """The syndrome (one bit per stabilizer generator, 1 = anticommutes /
    detected, 0 = commutes / undetected) a given Pauli error string would
    produce against `stabilizers` (a list of equal-length Pauli strings,
    the code's stabilizer generators -- X-type, Z-type, or mixed; this
    function doesn't assume a CSS structure)."""
    return tuple(0 if pauli_commutes(pauli_error, g) else 1 for g in stabilizers)

erasure_aware_decode 

erasure_aware_decode(
    observed_syndrome: tuple,
    heralded_qubits: Sequence[int],
    n_qubits: int,
    stabilizers: Sequence[str],
) -> Optional[str]

Erasure-aware decoder for any stabilizer code. Given the observed syndrome and a list of qubits KNOWN to have been erased (e.g. a heralded photon-loss event on a dual-rail-encoded qubit), brute-forces every Pauli assignment (I/X/Y/Z, 4**len(heralded_qubits) combinations) on just the heralded qubits and returns the unique full-length Pauli string reproducing observed_syndrome exactly.

Returns None -- not a guess -- when there are zero heralded qubits, when the observed syndrome is not explained by any assignment on the heralded qubits alone, or when more than one assignment explains it (ambiguous). Both None cases mean: fall back to a standard syndrome-only decoder, or treat as a detected-but-uncorrectable event -- this function will not silently return a wrong-but-plausible correction. The number of heralded qubits this can actually resolve unambiguously is bounded by the code's real distance (Grassl, Beth & Pellizzari 1997: up to d-1 erasures) -- that bound emerges naturally from the brute-force search itself (more heralded qubits than the code can resolve typically yields zero or multiple matches), it is not hard-coded here.

Cost is 4**len(heralded_qubits) syndrome computations, each O(n_qubits * len(stabilizers)) -- fine for the small numbers of simultaneous erasures a real per-shot noise rate produces (verified up to 2 in the original Steane investigation; tractable up to 4-5 for most small codes before it's worth switching to a smarter search).

Source code in dense_evolution/physics/qec.py 
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
def erasure_aware_decode(
    observed_syndrome: tuple,
    heralded_qubits: Sequence[int],
    n_qubits: int,
    stabilizers: Sequence[str],
) -> Optional[str]:
    """Erasure-aware decoder for any stabilizer code. Given the observed
    syndrome and a list of qubits KNOWN to have been erased (e.g. a
    heralded photon-loss event on a dual-rail-encoded qubit), brute-forces
    every Pauli assignment (I/X/Y/Z, 4**len(heralded_qubits) combinations)
    on just the heralded qubits and returns the unique full-length Pauli
    string reproducing `observed_syndrome` exactly.

    Returns `None` -- not a guess -- when there are zero heralded qubits,
    when the observed syndrome is not explained by any assignment on the
    heralded qubits alone, or when more than one assignment explains it
    (ambiguous). Both `None` cases mean: fall back to a standard
    syndrome-only decoder, or treat as a detected-but-uncorrectable
    event -- this function will not silently return a wrong-but-plausible
    correction. The number of heralded qubits this can actually resolve
    unambiguously is bounded by the code's real distance (Grassl, Beth &
    Pellizzari 1997: up to d-1 erasures) -- that bound emerges naturally
    from the brute-force search itself (more heralded qubits than the
    code can resolve typically yields zero or multiple matches), it is
    not hard-coded here.

    Cost is 4**len(heralded_qubits) syndrome computations, each O(n_qubits
    * len(stabilizers)) -- fine for the small numbers of simultaneous
    erasures a real per-shot noise rate produces (verified up to 2 in the
    original Steane investigation; tractable up to 4-5 for most small
    codes before it's worth switching to a smarter search).
    """
    if not heralded_qubits:
        return None

    matches = []
    for assignment_paulis in itertools.product(PAULIS, repeat=len(heralded_qubits)):
        assignment = dict(zip(heralded_qubits, assignment_paulis))
        candidate = _pauli_string(n_qubits, assignment)
        if compute_syndrome(candidate, stabilizers) == tuple(observed_syndrome):
            matches.append(candidate)

    if len(matches) != 1:
        return None
    return matches[0]

pymatching_decode 

pymatching_decode(
    observed_syndrome: Sequence[int],
    stabilizers: Sequence[str],
    n_qubits: int,
    error_type: str = "X",
    weights: Optional[Sequence[float]] = None,
) -> str

MWPM syndrome decoder (via pymatching) for a single physical error type, no erasure/heralding information needed -- the standard decoding setting erasure_aware_decode deliberately doesn't cover (it always returns None with zero heralded qubits).

Restricted to same-purpose stabilizers detecting ONE error type at a time (the standard CSS setup: e.g. Z-type stabilizers decoding X errors, or vice versa) -- NOT the fully general mixed-Pauli case compute_syndrome/erasure_aware_decode accept. A code with both X- and Z-type errors to correct needs two separate calls, one per type (X-type stabilizers -> decode Z errors, Z-type stabilizers -> decode X errors), each contributing its own half of the full correction -- the same two-pass structure any real CSS decoder uses, not a limitation specific to this wrapper.

FURTHER REAL RESTRICTION (verified directly against pymatching, not assumed from its docs): a matching-graph decoder needs every qubit checked by AT MOST 2 stabilizers -- pymatching represents each potential error as a graph EDGE between (at most) 2 detector nodes, a structural fact about topological/surface codes (each qubit sits on an edge between exactly 2 checks), not a general property of stabilizer codes. Steane [[7,1,3]]'s weight-4 stabilizers are a real counterexample -- qubit 6 is checked by all 3 X-stabilizers at once (3 > 2), so pymatching_decode cannot be used for it at all (confirmed: raises ValueError below, not just slow or approximate) -- use erasure_aware_decode for codes like that instead. Repetition codes and the surface/toric code family satisfy the <=2 constraint by construction; small non-topological codes generally do not.

The check matrix pymatching needs is built from stabilizers via this module's own pauli_commutes (column q of row i is 1 iff stabilizer i anticommutes with a lone error_type error on qubit q) rather than a naive "non-identity entry" heuristic -- those disagree whenever a stabilizer has the SAME letter as error_type at some qubit (e.g. an 'X' entry in a stabilizer being used to decode X errors: it commutes with an X error there and must NOT count as detecting it, but is very much non-identity).

Parameters:

Name Type Description Default
observed_syndrome sequence of int

One bit per stabilizer, same convention as compute_syndrome's return value (1 = that stabilizer detected/anticommuted).

 required
stabilizers sequence of str

The stabilizer generators used for this error type (e.g. every Z-type generator, to decode X errors), each a length-n_qubits Pauli string over IXYZ. Mixing generator types that detect different error types in the same call silently produces a wrong check matrix -- this function has no way to tell that apart from a correctly-scoped single-type list, so it isn't validated here.

 required
n_qubits int

Number of physical qubits (columns of the check matrix / length of the returned Pauli string).

 required
error_type str

Which single-qubit Pauli error this decodes for -- one of 'X', 'Y', 'Z'. Defaults to 'X'.

 'X'
weights sequence of float

Per-qubit edge weight for MWPM (e.g. -log(p_q) for a known per-qubit physical error rate p_q) -- forwarded to pymatching.Matching.from_check_matrix. Defaults to pymatching's own default (uniform weight 1.0 for every qubit, i.e. no prior assumption about which qubits are more error-prone).

 None

Returns:

Type Description
str

Length-n_qubits Pauli string (only 'I' and error_type) giving the minimum-weight correction pymatching found.

Raises:

Type Description
ImportError

If pymatching isn't installed (pip install dense-evolution[pymatching]).
ValueError

If error_type isn't one of 'X'/'Y'/'Z', if observed_syndrome doesn't have one entry per stabilizer, if a stabilizer's length isn't n_qubits, or if every stabilizer commutes with every possible error_type error (the check matrix would be all-zero -- almost always means stabilizers is the wrong generator type for the requested error_type, not a real all-zero code).

Examples:

>>> # 3-qubit repetition code, Z-type stabilizers, decoding X errors
>>> stabilizers = ['ZZI', 'IZZ']
>>> syndrome = compute_syndrome('IXI', stabilizers)  # X error on qubit 1
>>> pymatching_decode(syndrome, stabilizers, n_qubits=3, error_type='X')
'IXI'
Source code in dense_evolution/physics/qec.py 
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
def pymatching_decode(
    observed_syndrome: Sequence[int],
    stabilizers: Sequence[str],
    n_qubits: int,
    error_type: str = 'X',
    weights: Optional[Sequence[float]] = None,
) -> str:
    """MWPM syndrome decoder (via `pymatching`) for a single physical error
    type, no erasure/heralding information needed -- the standard decoding
    setting `erasure_aware_decode` deliberately doesn't cover (it always
    returns None with zero heralded qubits).

    Restricted to same-purpose stabilizers detecting ONE error type at a
    time (the standard CSS setup: e.g. Z-type stabilizers decoding X
    errors, or vice versa) -- NOT the fully general mixed-Pauli case
    `compute_syndrome`/`erasure_aware_decode` accept. A code with both X-
    and Z-type errors to correct needs two separate calls, one per type
    (X-type stabilizers -> decode Z errors, Z-type stabilizers -> decode
    X errors), each contributing its own half of the full correction --
    the same two-pass structure any real CSS decoder uses, not a
    limitation specific to this wrapper.

    FURTHER REAL RESTRICTION (verified directly against pymatching, not
    assumed from its docs): a matching-graph decoder needs every qubit
    checked by AT MOST 2 stabilizers -- pymatching represents each
    potential error as a graph EDGE between (at most) 2 detector nodes, a
    structural fact about topological/surface codes (each qubit sits on
    an edge between exactly 2 checks), not a general property of
    stabilizer codes. Steane [[7,1,3]]'s weight-4 stabilizers are a real
    counterexample -- qubit 6 is checked by all 3 X-stabilizers at once
    (3 > 2), so `pymatching_decode` cannot be used for it at all
    (confirmed: raises ValueError below, not just slow or approximate) --
    use `erasure_aware_decode` for codes like that instead. Repetition
    codes and the surface/toric code family satisfy the <=2 constraint by
    construction; small non-topological codes generally do not.

    The check matrix pymatching needs is built from `stabilizers` via this
    module's own `pauli_commutes` (column q of row i is 1 iff stabilizer i
    anticommutes with a lone `error_type` error on qubit q) rather than a
    naive "non-identity entry" heuristic -- those disagree whenever a
    stabilizer has the SAME letter as `error_type` at some qubit (e.g. an
    'X' entry in a stabilizer being used to decode X errors: it commutes
    with an X error there and must NOT count as detecting it, but is very
    much non-identity).

    Parameters
    ----------
    observed_syndrome : sequence of int
        One bit per stabilizer, same convention as `compute_syndrome`'s
        return value (1 = that stabilizer detected/anticommuted).
    stabilizers : sequence of str
        The stabilizer generators used for this error type (e.g. every
        Z-type generator, to decode X errors), each a length-`n_qubits`
        Pauli string over IXYZ. Mixing generator types that detect
        different error types in the same call silently produces a wrong
        check matrix -- this function has no way to tell that apart from
        a correctly-scoped single-type list, so it isn't validated here.
    n_qubits : int
        Number of physical qubits (columns of the check matrix / length
        of the returned Pauli string).
    error_type : str, optional
        Which single-qubit Pauli error this decodes for -- one of 'X',
        'Y', 'Z'. Defaults to 'X'.
    weights : sequence of float, optional
        Per-qubit edge weight for MWPM (e.g. -log(p_q) for a known
        per-qubit physical error rate p_q) -- forwarded to
        `pymatching.Matching.from_check_matrix`. Defaults to pymatching's
        own default (uniform weight 1.0 for every qubit, i.e. no prior
        assumption about which qubits are more error-prone).

    Returns
    -------
    str
        Length-`n_qubits` Pauli string (only 'I' and `error_type`) giving
        the minimum-weight correction pymatching found.

    Raises
    ------
    ImportError
        If `pymatching` isn't installed (`pip install dense-evolution[pymatching]`).
    ValueError
        If `error_type` isn't one of 'X'/'Y'/'Z', if `observed_syndrome`
        doesn't have one entry per stabilizer, if a stabilizer's length
        isn't `n_qubits`, or if every stabilizer commutes with every
        possible `error_type` error (the check matrix would be all-zero --
        almost always means `stabilizers` is the wrong generator type for
        the requested `error_type`, not a real all-zero code).

    Examples
    --------
    >>> # 3-qubit repetition code, Z-type stabilizers, decoding X errors
    >>> stabilizers = ['ZZI', 'IZZ']
    >>> syndrome = compute_syndrome('IXI', stabilizers)  # X error on qubit 1
    >>> pymatching_decode(syndrome, stabilizers, n_qubits=3, error_type='X')
    'IXI'
    """
    _require_pymatching()

    if error_type not in ('X', 'Y', 'Z'):
        raise ValueError(f"error_type must be one of 'X', 'Y', 'Z', got {error_type!r}")
    if len(observed_syndrome) != len(stabilizers):
        raise ValueError(
            f"observed_syndrome has {len(observed_syndrome)} entries but there are "
            f"{len(stabilizers)} stabilizers -- these must match one-to-one"
        )
    for i, s in enumerate(stabilizers):
        if len(s) != n_qubits:
            raise ValueError(f"stabilizers[{i}] has length {len(s)}, expected n_qubits={n_qubits}")

    check_matrix = np.array(
        [[0 if pauli_commutes(s[q], error_type) else 1 for q in range(n_qubits)] for s in stabilizers],
        dtype=np.uint8,
    )
    if not check_matrix.any():
        raise ValueError(
            f"every stabilizer commutes with every possible {error_type} error -- "
            f"stabilizers is very likely the wrong generator type to decode {error_type} "
            f"errors (e.g. passing X-type stabilizers to decode X errors, which they "
            f"cannot detect by construction)"
        )
    checks_per_qubit = check_matrix.sum(axis=0)
    if (checks_per_qubit > 2).any():
        bad_qubits = np.where(checks_per_qubit > 2)[0].tolist()
        raise ValueError(
            f"pymatching's matching-graph decoder needs every qubit checked by AT MOST 2 "
            f"stabilizers (a graph edge connects at most 2 detector nodes) -- qubit(s) "
            f"{bad_qubits} are each checked by {checks_per_qubit[bad_qubits].tolist()} "
            f"stabilizers here. This is a real structural requirement of matching-graph "
            f"decoding (true for the surface code and other topological codes, where each "
            f"qubit sits on an edge between exactly 2 checks), NOT satisfied by every "
            f"stabilizer code -- e.g. Steane [[7,1,3]]'s weight-4 stabilizers check some "
            f"qubits 3 times, so pymatching_decode cannot be used for it; use "
            f"erasure_aware_decode instead for codes like that."
        )

    matching = pymatching.Matching.from_check_matrix(check_matrix, weights=weights)
    correction = matching.decode(np.asarray(observed_syndrome, dtype=np.uint8))

    return ''.join(error_type if bit else 'I' for bit in correction)

blind_minimum_weight_decode 

blind_minimum_weight_decode(
    observed_syndrome: Sequence[int],
    n_qubits: int,
    stabilizers: Sequence[str],
    max_weight: Optional[int] = None,
) -> Optional[str]

Blind (no known erasure locations) minimum-weight decoder for any stabilizer code -- including ones pymatching_decode structurally cannot handle (Steane and other small non-"graph-like" codes; see this module's docstring for why the naive erasure_aware_decode( ..., heralded_qubits=range(n_qubits)) trick does NOT work here).

Searches every possible Pauli error in increasing WEIGHT order (0 non-identity qubits, then 1, then 2, ...), stopping at the first weight with at least one match. Multiple matches at that weight are NOT automatically ambiguous: in a degenerate code (e.g. Shor [[9,1,3]]), two lowest-weight corrections that differ by a stabilizer-group element are the SAME physical correction (applying either one restores the code space identically), not two competing guesses. This checks that directly -- via the corrections' symplectic (GF(2)) representation, not by re-deriving it from scratch each call -- and returns the (deterministic, first-found) representative when every match is in the same stabilizer coset. Returns None only when two matches at the minimum weight differ by an actual LOGICAL operator (not in the stabilizer group), which is genuine ambiguity, or when nothing matches by max_weight. This IS the same principle pymatching_decode's MWPM implements, via brute force instead of a matching graph -- deliberately much slower, in exchange for working on ANY stabilizer code, graph-like or not.

Cost: for a given weight w, C(n_qubits, w) * 3**w syndrome checks, each O(n_qubits * len(stabilizers)) -- explodes quickly for w beyond a handful (e.g. n_qubits=7: 21 at w=1, 189 at w=2, 945 at w=3), so this is for the small-code, low-weight-error regime erasure_aware_decode already targets, not a general substitute for pymatching_decode at surface-code sizes.

Parameters:

Name Type Description Default
observed_syndrome sequence of int

One bit per stabilizer, same convention as compute_syndrome.

 required
n_qubits int

Number of physical qubits.

 required
stabilizers sequence of str

The code's stabilizer generators (any mix of Pauli types -- unlike pymatching_decode, not restricted to one error type per call: this searches full IXYZ Pauli strings directly, same as erasure_aware_decode).

 required
max_weight int

Largest error weight to search before giving up. Defaults to n_qubits (search everything) -- pass a small explicit bound (matching the code's known error-correcting capability, e.g. 1 for a distance-3 code) to fail fast instead of paying the full combinatorial cost on a syndrome nothing low-weight can explain.

 None

Returns:

Type Description
str or None

Length-n_qubits Pauli string (IXYZ), or None if ambiguous at the minimum matching weight or unexplained up to max_weight.

Examples:

>>> # Steane [[7,1,3]], full X+Z stabilizer set -- exactly the code
>>> # pymatching_decode CANNOT handle (weight-4 stabilizers check some
>>> # qubits 3 times), blind single-qubit Z error, no erasure info:
>>> steane_x = ['IIIXXXX', 'IXXIIXX', 'XIXIXIX']
>>> steane_z = ['IIIZZZZ', 'IZZIIZZ', 'ZIZIZIZ']
>>> stabilizers = steane_x + steane_z
>>> syndrome = compute_syndrome('IIIZIII', stabilizers)
>>> blind_minimum_weight_decode(syndrome, n_qubits=7, stabilizers=stabilizers)
'IIIZIII'
Source code in dense_evolution/physics/qec.py 
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
def blind_minimum_weight_decode(
    observed_syndrome: Sequence[int],
    n_qubits: int,
    stabilizers: Sequence[str],
    max_weight: Optional[int] = None,
) -> Optional[str]:
    """Blind (no known erasure locations) minimum-weight decoder for any
    stabilizer code -- including ones `pymatching_decode` structurally
    cannot handle (Steane and other small non-"graph-like" codes; see
    this module's docstring for why the naive `erasure_aware_decode(
    ..., heralded_qubits=range(n_qubits))` trick does NOT work here).

    Searches every possible Pauli error in increasing WEIGHT order (0
    non-identity qubits, then 1, then 2, ...), stopping at the first
    weight with at least one match. Multiple matches at that weight are
    NOT automatically ambiguous: in a degenerate code (e.g. Shor
    [[9,1,3]]), two lowest-weight corrections that differ by a
    stabilizer-group element are the SAME physical correction (applying
    either one restores the code space identically), not two competing
    guesses. This checks that directly -- via the corrections' symplectic
    (GF(2)) representation, not by re-deriving it from scratch each call
    -- and returns the (deterministic, first-found) representative when
    every match is in the same stabilizer coset. Returns `None` only when
    two matches at the minimum weight differ by an actual LOGICAL
    operator (not in the stabilizer group), which is genuine ambiguity, or
    when nothing matches by `max_weight`. This IS the same principle
    `pymatching_decode`'s MWPM implements, via brute force instead of a
    matching graph -- deliberately much slower, in exchange for working
    on ANY stabilizer code, graph-like or not.

    Cost: for a given weight w, `C(n_qubits, w) * 3**w` syndrome checks,
    each O(n_qubits * len(stabilizers)) -- explodes quickly for w beyond
    a handful (e.g. n_qubits=7: 21 at w=1, 189 at w=2, 945 at w=3), so
    this is for the small-code, low-weight-error regime
    `erasure_aware_decode` already targets, not a general substitute for
    `pymatching_decode` at surface-code sizes.

    Parameters
    ----------
    observed_syndrome : sequence of int
        One bit per stabilizer, same convention as `compute_syndrome`.
    n_qubits : int
        Number of physical qubits.
    stabilizers : sequence of str
        The code's stabilizer generators (any mix of Pauli types --
        unlike `pymatching_decode`, not restricted to one error type per
        call: this searches full IXYZ Pauli strings directly, same as
        `erasure_aware_decode`).
    max_weight : int, optional
        Largest error weight to search before giving up. Defaults to
        `n_qubits` (search everything) -- pass a small explicit bound
        (matching the code's known error-correcting capability, e.g. 1
        for a distance-3 code) to fail fast instead of paying the full
        combinatorial cost on a syndrome nothing low-weight can explain.

    Returns
    -------
    str or None
        Length-`n_qubits` Pauli string (IXYZ), or `None` if ambiguous at
        the minimum matching weight or unexplained up to `max_weight`.

    Examples
    --------
    >>> # Steane [[7,1,3]], full X+Z stabilizer set -- exactly the code
    >>> # pymatching_decode CANNOT handle (weight-4 stabilizers check some
    >>> # qubits 3 times), blind single-qubit Z error, no erasure info:
    >>> steane_x = ['IIIXXXX', 'IXXIIXX', 'XIXIXIX']
    >>> steane_z = ['IIIZZZZ', 'IZZIIZZ', 'ZIZIZIZ']
    >>> stabilizers = steane_x + steane_z
    >>> syndrome = compute_syndrome('IIIZIII', stabilizers)
    >>> blind_minimum_weight_decode(syndrome, n_qubits=7, stabilizers=stabilizers)
    'IIIZIII'
    """
    if max_weight is None:
        max_weight = n_qubits
    if not 0 <= max_weight <= n_qubits:
        raise ValueError(f"max_weight must be between 0 and n_qubits={n_qubits}, got {max_weight}")
    if len(observed_syndrome) != len(stabilizers):
        raise ValueError(
            f"observed_syndrome has {len(observed_syndrome)} entries but there are "
            f"{len(stabilizers)} stabilizers -- these must match one-to-one"
        )
    for i, s in enumerate(stabilizers):
        if len(s) != n_qubits:
            raise ValueError(f"stabilizers[{i}] has length {len(s)}, expected n_qubits={n_qubits}")

    target = tuple(observed_syndrome)
    stab_rref, stab_pivots = _gf2_rref(
        np.array([_pauli_to_symplectic(s) for s in stabilizers], dtype=np.uint8)
    )

    for weight in range(max_weight + 1):
        matches = []
        for qubits in itertools.combinations(range(n_qubits), weight):
            for paulis in itertools.product(PAULIS[1:], repeat=weight):
                assignment = dict(zip(qubits, paulis))
                candidate = _pauli_string(n_qubits, assignment)
                if compute_syndrome(candidate, stabilizers) == target:
                    matches.append(candidate)
        if not matches:
            continue
        if len(matches) == 1:
            return matches[0]

        reference = matches[0]
        reference_symplectic = _pauli_to_symplectic(reference)
        same_coset = all(
            _in_gf2_span(reference_symplectic ^ _pauli_to_symplectic(m), stab_rref, stab_pivots)
            for m in matches[1:]
        )
        return reference if same_coset else None

    return None

counts_in_intervals_dimension 

counts_in_intervals_dimension(
    event_times: Sequence[float],
    window_sizes: Sequence[float],
    min_reference_points: int = 5,
) -> tuple

Correlation-dimension estimator for a 1-D point process (event times), generalizing the "counts-in-spheres" statistic used to measure the fractal dimension of galaxy distributions -- Scrimgeour et al., "The WiggleZ Dark Energy Survey: the transition to large-scale cosmic homogeneity," MNRAS 425, 116 (2012), arXiv:1205.6812 -- from 3-D space down to 1-D time.

For a homogeneous (Poisson) point process, the expected number of OTHER events within radius r of a given event scales as N(<=r) ~ r^1 exactly. Real burst-like noise (e.g. cosmic-ray-correlated error events, arXiv:2104.05219, already modelled by cosmic_ray_burst_profile in dense_evolution.mitigation.zne) clusters in time: events arrive in tight groups separated by long, comparatively empty gaps, which depresses this exponent below 1. This function measures that exponent directly from a sequence of observed/simulated event times (e.g. heralded-erasure timestamps, or detected-syndrome timestamps), instead of assuming Poissonian noise or a particular burst model up front -- the measured dimension is then evidence for whether decode_with_erasure_fallback style single-shot decoding or a burst-aware strategy is the right model for a given noise source.

For each candidate radius r in window_sizes, every event at least r away from both ends of the observed time range is used as a reference point (events too close to either edge are skipped for that r, so a partially-empty window near the boundary never silently deflates the count -- the same edge correction real counts-in-spheres analyses use). The mean count of other events within r of each valid reference point is computed, and the dimension D is the slope of log(mean count) vs log(r) over a linear least-squares fit -- with the fit's R^2 returned alongside it, not hidden, because a narrow or poorly-covered range of window_sizes can make the fitted slope meaningless (see the warning below).

D ~= 1 : homogeneous/Poissonian arrivals, no clustering. D < 1 : temporally clustered/bursty (e.g. cosmic-ray-correlated error bursts) -- the smaller D, the tighter the clustering. D > 1 : more regularly spaced than random (suppressed fluctuations, "hyperuniform" arrivals) -- unusual for physical error processes but not excluded by the statistic itself.

A LOW R^2 (rule of thumb: below ~0.98) means window_sizes does not span enough dynamic range or lacks enough valid reference points for the fit to be trustworthy -- widen the range (ideally 2-3 orders of magnitude) and/or supply more events rather than trusting the reported D. This mirrors a real, previously-made mistake: a spatial box-counting fit through only 4 points spanning a narrow range of scales produced a fractal dimension of 142 (physically impossible in 3 dimensions) from noise alone, not from any real structure in the data -- always inspect R^2 before quoting D.

Cost is O(len(window_sizes) * n_events * log(n_events)) -- event_times is sorted once up front, and each radius's per-reference-point count is a pair of np.searchsorted calls (vectorized across all reference points at once) rather than an O(n_events) brute-force scan per point.

Parameters:

Name Type Description Default
event_times sequence of float

Timestamps of observed/simulated events (need not be sorted).

 required
window_sizes sequence of float

Radii r to evaluate counts at, ideally spanning several orders of magnitude and containing at least ~5-8 values.

 required
min_reference_points int

Minimum number of edge-safe reference events required at a given r for that r to be included in the fit. Radii with fewer valid reference points (or a zero mean count) are silently dropped, not zero-padded, so the returned mean_counts may be shorter than window_sizes. Defaults to 5.

 5

Returns:

Name Type Description
dimension float

Estimated scaling exponent D (the slope of the log-log fit).
r_squared float

Coefficient of determination of the log-log linear fit -- inspect this before trusting dimension (see above).
mean_counts dict[float, float]

Mean count of other events within each usable radius r, keyed by the r values from window_sizes that had enough valid reference points and a nonzero count.

Raises:

Type Description
ValueError

If event_times has fewer than 2 events, if any window_sizes entry is not positive, or if fewer than 2 radii end up with enough valid reference points and a nonzero count to fit a slope at all.

Examples:

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> poisson_events = np.sort(rng.uniform(0, 1000, 2000))
>>> radii = np.logspace(0, 2, 10)  # 1 to 100
>>> D, r2, _ = counts_in_intervals_dimension(poisson_events, radii)
>>> 0.9 < D < 1.1 and r2 > 0.99
True
Source code in dense_evolution/physics/qec.py 
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
def counts_in_intervals_dimension(
    event_times: Sequence[float],
    window_sizes: Sequence[float],
    min_reference_points: int = 5,
) -> tuple:
    """Correlation-dimension estimator for a 1-D point process (event
    times), generalizing the "counts-in-spheres" statistic used to
    measure the fractal dimension of galaxy distributions -- Scrimgeour
    et al., "The WiggleZ Dark Energy Survey: the transition to
    large-scale cosmic homogeneity," MNRAS 425, 116 (2012),
    arXiv:1205.6812 -- from 3-D space down to 1-D time.

    For a homogeneous (Poisson) point process, the expected number of
    OTHER events within radius r of a given event scales as N(<=r) ~ r^1
    exactly. Real burst-like noise (e.g. cosmic-ray-correlated error
    events, arXiv:2104.05219, already modelled by `cosmic_ray_burst_profile`
    in `dense_evolution.mitigation.zne`) clusters in time: events arrive
    in tight groups separated by long, comparatively empty gaps, which
    depresses this exponent below 1. This function measures that exponent
    directly from a sequence of observed/simulated event times (e.g.
    heralded-erasure timestamps, or detected-syndrome timestamps), instead
    of assuming Poissonian noise or a particular burst model up front --
    the measured dimension is then evidence for whether `decode_with_erasure_fallback`
    style single-shot decoding or a burst-aware strategy is the right model
    for a given noise source.

    For each candidate radius r in `window_sizes`, every event at least r
    away from both ends of the observed time range is used as a reference
    point (events too close to either edge are skipped for that r, so a
    partially-empty window near the boundary never silently deflates the
    count -- the same edge correction real counts-in-spheres analyses
    use). The mean count of other events within r of each valid reference
    point is computed, and the dimension D is the slope of log(mean
    count) vs log(r) over a linear least-squares fit -- with the fit's
    R^2 returned alongside it, not hidden, because a narrow or
    poorly-covered range of `window_sizes` can make the fitted slope
    meaningless (see the warning below).

    D ~= 1 : homogeneous/Poissonian arrivals, no clustering.
    D < 1 : temporally clustered/bursty (e.g. cosmic-ray-correlated
            error bursts) -- the smaller D, the tighter the clustering.
    D > 1 : more regularly spaced than random (suppressed fluctuations,
            "hyperuniform" arrivals) -- unusual for physical error
            processes but not excluded by the statistic itself.

    A LOW R^2 (rule of thumb: below ~0.98) means `window_sizes` does not
    span enough dynamic range or lacks enough valid reference points for
    the fit to be trustworthy -- widen the range (ideally 2-3 orders of
    magnitude) and/or supply more events rather than trusting the
    reported D. This mirrors a real, previously-made mistake: a spatial
    box-counting fit through only 4 points spanning a narrow range of
    scales produced a fractal dimension of 142 (physically impossible in
    3 dimensions) from noise alone, not from any real structure in the
    data -- always inspect R^2 before quoting D.

    Cost is O(len(window_sizes) * n_events * log(n_events)) -- `event_times`
    is sorted once up front, and each radius's per-reference-point count
    is a pair of `np.searchsorted` calls (vectorized across all reference
    points at once) rather than an O(n_events) brute-force scan per point.

    Parameters
    ----------
    event_times : sequence of float
        Timestamps of observed/simulated events (need not be sorted).
    window_sizes : sequence of float
        Radii r to evaluate counts at, ideally spanning several orders of
        magnitude and containing at least ~5-8 values.
    min_reference_points : int, optional
        Minimum number of edge-safe reference events required at a given
        r for that r to be included in the fit. Radii with fewer valid
        reference points (or a zero mean count) are silently dropped, not
        zero-padded, so the returned `mean_counts` may be shorter than
        `window_sizes`. Defaults to 5.

    Returns
    -------
    dimension : float
        Estimated scaling exponent D (the slope of the log-log fit).
    r_squared : float
        Coefficient of determination of the log-log linear fit --
        inspect this before trusting `dimension` (see above).
    mean_counts : dict[float, float]
        Mean count of other events within each usable radius r, keyed by
        the r values from `window_sizes` that had enough valid reference
        points and a nonzero count.

    Raises
    ------
    ValueError
        If `event_times` has fewer than 2 events, if any `window_sizes`
        entry is not positive, or if fewer than 2 radii end up with
        enough valid reference points and a nonzero count to fit a slope
        at all.

    Examples
    --------
    >>> import numpy as np
    >>> rng = np.random.default_rng(0)
    >>> poisson_events = np.sort(rng.uniform(0, 1000, 2000))
    >>> radii = np.logspace(0, 2, 10)  # 1 to 100
    >>> D, r2, _ = counts_in_intervals_dimension(poisson_events, radii)
    >>> 0.9 < D < 1.1 and r2 > 0.99
    True
    """
    event_times = np.asarray(event_times, dtype=float)
    if event_times.size < 2:
        raise ValueError(f"need at least 2 events, got {event_times.size}")
    event_times = np.sort(event_times)
    t_min, t_max = event_times[0], event_times[-1]

    mean_counts = {}
    for r in window_sizes:
        if r <= 0:
            raise ValueError(f"window_sizes must be positive, got {r}")
        valid_refs = event_times[(event_times - r >= t_min) & (event_times + r <= t_max)]
        if valid_refs.size < min_reference_points:
            continue
        # event_times is sorted (see above), so "count within r of t" is a
        # pair of binary searches instead of an O(n) scan per reference
        # point (prog.txt point 5c) -- both bounds inclusive, matching the
        # brute-force `abs(event_times - t) <= r` this replaces exactly
        # (verified: side='left'/'right' at t-r/t+r reproduces it for
        # every t, including ties exactly on the r boundary).
        lo = np.searchsorted(event_times, valid_refs - r, side='left')
        hi = np.searchsorted(event_times, valid_refs + r, side='right')
        counts = (hi - lo) - 1
        mean_counts[float(r)] = float(np.mean(counts))

    valid_items = sorted((r, c) for r, c in mean_counts.items() if c > 0)
    if len(valid_items) < 2:
        raise ValueError(
            f"only {len(valid_items)} of {len(window_sizes)} window_sizes had at least "
            f"{min_reference_points} edge-safe reference points with a nonzero count -- "
            f"widen window_sizes, supply more events, or lower min_reference_points"
        )

    r_vals = np.array([r for r, _ in valid_items])
    n_vals = np.array([c for _, c in valid_items])
    log_r, log_n = np.log(r_vals), np.log(n_vals)
    design = np.vstack([log_r, np.ones_like(log_r)]).T
    slope, intercept = np.linalg.lstsq(design, log_n, rcond=None)[0]
    predicted = slope * log_r + intercept
    ss_res = np.sum((log_n - predicted) ** 2)
    ss_tot = np.sum((log_n - np.mean(log_n)) ** 2)
    r_squared = 1.0 - ss_res / ss_tot if ss_tot > 0 else 1.0

    return float(slope), float(r_squared), dict(valid_items)

decode_with_erasure_fallback 

decode_with_erasure_fallback(
    observed_syndrome: Sequence[int],
    heralded_qubits: Sequence[int],
    n_qubits: int,
    stabilizers: Sequence[str],
    max_weight: Optional[int] = None,
) -> Optional[str]

The real-world decoding POLICY, not just a raw decoder call: use erasure_aware_decode when there are any heralded qubits and it resolves the syndrome uniquely; fall back to blind_minimum_weight_decode otherwise (zero heralded qubits, or erasure-aware decoding can't resolve it -- ambiguous, or more heralds than the code's real erasure-correcting capacity).

Promoted from Dense-Evolution-Discovery's cosmic-ray-burst-as-erasure experiment (scripts/cosmic_ray_erasure_decoding.py), where this exact fallback logic was first written inline in a Monte Carlo loop. Never worse than always calling blind_minimum_weight_decode directly -- it only uses herald information when doing so actually helps, exactly the policy a real erasure-aware QEC system would run, since a real-time detector (e.g. a cosmic-ray/particle-impact monitor, or a photon-loss herald) only ever adds information, it doesn't obligate a decoder to use it past the point where it stops being useful.

Parameters:

Name Type Description Default
observed_syndrome sequence of int

One bit per stabilizer, same convention as compute_syndrome.

 required
heralded_qubits sequence of int

Qubits KNOWN to have been erased/disturbed this shot -- may be empty (falls straight through to blind decoding).

 required
n_qubits int

Number of physical qubits.

 required
stabilizers sequence of str

The code's stabilizer generators (any mix of Pauli types).

 required
max_weight int

Forwarded to blind_minimum_weight_decode's fallback path only (see its own docstring) -- erasure_aware_decode has no equivalent parameter, its search is always exactly over the heralded qubits.

 None

Returns:

Type Description
str or None

Length-n_qubits Pauli string (IXYZ), or None if neither decoder can resolve the syndrome.
Source code in dense_evolution/physics/qec.py 
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
def decode_with_erasure_fallback(
    observed_syndrome: Sequence[int],
    heralded_qubits: Sequence[int],
    n_qubits: int,
    stabilizers: Sequence[str],
    max_weight: Optional[int] = None,
) -> Optional[str]:
    """The real-world decoding POLICY, not just a raw decoder call: use
    `erasure_aware_decode` when there are any heralded qubits and it
    resolves the syndrome uniquely; fall back to `blind_minimum_weight_decode`
    otherwise (zero heralded qubits, or erasure-aware decoding can't
    resolve it -- ambiguous, or more heralds than the code's real
    erasure-correcting capacity).

    Promoted from Dense-Evolution-Discovery's cosmic-ray-burst-as-erasure
    experiment (scripts/cosmic_ray_erasure_decoding.py), where this exact
    fallback logic was first written inline in a Monte Carlo loop. Never
    worse than always calling `blind_minimum_weight_decode` directly --
    it only uses herald information when doing so actually helps, exactly
    the policy a real erasure-aware QEC system would run, since a
    real-time detector (e.g. a cosmic-ray/particle-impact monitor, or a
    photon-loss herald) only ever adds information, it doesn't obligate a
    decoder to use it past the point where it stops being useful.

    Parameters
    ----------
    observed_syndrome : sequence of int
        One bit per stabilizer, same convention as `compute_syndrome`.
    heralded_qubits : sequence of int
        Qubits KNOWN to have been erased/disturbed this shot -- may be
        empty (falls straight through to blind decoding).
    n_qubits : int
        Number of physical qubits.
    stabilizers : sequence of str
        The code's stabilizer generators (any mix of Pauli types).
    max_weight : int, optional
        Forwarded to `blind_minimum_weight_decode`'s fallback path only
        (see its own docstring) -- `erasure_aware_decode` has no
        equivalent parameter, its search is always exactly over the
        heralded qubits.

    Returns
    -------
    str or None
        Length-`n_qubits` Pauli string (IXYZ), or `None` if neither
        decoder can resolve the syndrome.
    """
    if heralded_qubits:
        result = erasure_aware_decode(observed_syndrome, heralded_qubits, n_qubits, stabilizers)
        if result is not None:
            return result
    return blind_minimum_weight_decode(observed_syndrome, n_qubits, stabilizers, max_weight=max_weight)

nearest_coset_decode 

nearest_coset_decode(
    measured_bits: str,
    coset_a: Sequence[str],
    coset_b: Sequence[str],
) -> int

Nearest-coset (minimum Hamming distance) binary decoding: given a measured bit string, decide which of two disjoint cosets of bit strings it is closer to. Returns 0 if measured_bits is closer to coset_a, 1 if closer to coset_b (ties broken toward coset_a).

This decodes a LOGICAL READOUT, not a syndrome -- distinct from every other decoder in this module (pymatching_decode, blind_minimum_weight_decode, erasure_aware_decode, decode_with_erasure_fallback), which all turn a syndrome into an error correction. For a CSS code whose logical |0>_L is proportional to a superposition over one coset C of a classical linear code, and |1>_L over the complementary coset C+1...1, a physical measurement in the logical basis gives a bit string that exactly matches one codeword in the noiseless case, and the NEAREST one under noise -- exactly the decoding rule Huang, Zhu, Ippoliti, Monroe & Gullans (arXiv:2608.20676, "Continuous-angle logical rotations in the Steane code") use for their real Steane-code logical Ramsey experiment.

Promoted from Dense-Evolution-Discovery's real reproduction of that protocol (scripts/steane_continuous_logical_rotation.py) -- there, coset_a/coset_b were the real 8-codeword-each cosets read directly off this library's own Steane |0>_L / |1>_L statevectors (not assumed from the paper's text), and this decoding rule matched the paper's own theoretical logical-rotation model to within Monte Carlo statistical noise across 6000 real circuit trials (7 data qubits + 3 syndrome-extraction ancillas, real stochastic dephasing, real projective measurement collapse).

Parameters:

Name Type Description Default
measured_bits str

The measured bit string, e.g. '0101101'.

 required
coset_a sequence of str

Two disjoint sets of equal-length bit strings.

 required
coset_b sequence of str

Two disjoint sets of equal-length bit strings.

 required

Returns:

Type Description
int

0 or 1, indicating which coset measured_bits is nearest to.

Examples:

>>> from dense_evolution.qec import nearest_coset_decode
>>> coset_a = ['0000000', '1111000']
>>> coset_b = ['1111111', '0000111']
>>> nearest_coset_decode('0000000', coset_a, coset_b)
0
>>> nearest_coset_decode('1111110', coset_a, coset_b)
1
Source code in dense_evolution/physics/qec.py 
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
def nearest_coset_decode(measured_bits: str, coset_a: Sequence[str], coset_b: Sequence[str]) -> int:
    """Nearest-coset (minimum Hamming distance) binary decoding: given a
    measured bit string, decide which of two disjoint cosets of bit
    strings it is closer to. Returns 0 if `measured_bits` is closer to
    `coset_a`, 1 if closer to `coset_b` (ties broken toward `coset_a`).

    This decodes a LOGICAL READOUT, not a syndrome -- distinct from
    every other decoder in this module (`pymatching_decode`,
    `blind_minimum_weight_decode`, `erasure_aware_decode`,
    `decode_with_erasure_fallback`), which all turn a syndrome into an
    error correction. For a CSS code whose logical |0>_L is proportional
    to a superposition over one coset C of a classical linear code, and
    |1>_L over the complementary coset C+1...1, a physical measurement in
    the logical basis gives a bit string that exactly matches one
    codeword in the noiseless case, and the NEAREST one under noise --
    exactly the decoding rule Huang, Zhu, Ippoliti, Monroe & Gullans
    (arXiv:2608.20676, "Continuous-angle logical rotations in the Steane
    code") use for their real Steane-code logical Ramsey experiment.

    Promoted from Dense-Evolution-Discovery's real reproduction of that
    protocol (scripts/steane_continuous_logical_rotation.py) -- there,
    `coset_a`/`coset_b` were the real 8-codeword-each cosets read
    directly off this library's own Steane |0>_L / |1>_L statevectors
    (not assumed from the paper's text), and this decoding rule matched
    the paper's own theoretical logical-rotation model to within Monte
    Carlo statistical noise across 6000 real circuit trials (7 data
    qubits + 3 syndrome-extraction ancillas, real stochastic dephasing,
    real projective measurement collapse).

    Parameters
    ----------
    measured_bits : str
        The measured bit string, e.g. '0101101'.
    coset_a, coset_b : sequence of str
        Two disjoint sets of equal-length bit strings.

    Returns
    -------
    int
        0 or 1, indicating which coset `measured_bits` is nearest to.

    Examples
    --------
    >>> from dense_evolution.qec import nearest_coset_decode
    >>> coset_a = ['0000000', '1111000']
    >>> coset_b = ['1111111', '0000111']
    >>> nearest_coset_decode('0000000', coset_a, coset_b)
    0
    >>> nearest_coset_decode('1111110', coset_a, coset_b)
    1
    """
    def _min_hamming_distance(bits, coset):
        x = int(bits, 2)
        return min(bin(x ^ int(c, 2)).count('1') for c in coset)

    d_a = _min_hamming_distance(measured_bits, coset_a)
    d_b = _min_hamming_distance(measured_bits, coset_b)
    return 0 if d_a <= d_b else 1

estimate_edge_probabilities_from_detection_events 

estimate_edge_probabilities_from_detection_events(
    check_matrix, events
) -> np.ndarray

Error probability of every qubit, from detection events (Spitz et al.).

Implements the exact inversion of S. T. Spitz, B. Tarasinski, C. W. J. Beenakker and T. E. O'Brien, "Adaptive weight estimator for quantum error correction in a time-dependent environment", arXiv:1712.02360, Eqs. (13) and (16). For a code in which every qubit is checked by at most two checks (repetition and surface codes), each qubit is an edge between two checks, or between one check and the boundary. A qubit shared by checks i and j has probability

p = 1/2 - sqrt(1/4 - (<v_i v_j> - <v_i><v_j>) / (1 - 2 <v_i xor v_j>))

where v are the detection events and <.> the average over cycles. A qubit on the boundary of check i has

p = 1/2 + (<v_i> - 1/2) / prod(1 - 2 p_ij) over the other qubits of i.

The result can be passed as weights (-log(p / (1 - p))) to pymatching_decode.

Parameters:

Name Type Description Default
check_matrix array_like of shape (n_checks, n_qubits)

0/1 matrix, entry [c, q] = 1 when check c detects an error on qubit q. Every column has one or two ones; no two qubits may join the same pair of checks.

 required
events array_like of shape (n_cycles, n_checks)

0/1 detection events: 1 when a check changed value since the previous cycle (the syndrome of one cycle's new errors).

 required

Returns:

Type Description
numpy.ndarray of shape (n_qubits,)

Estimated error probability of each qubit, clipped to [0, 0.5].

Raises:

Type Description
ValueError

If the shapes do not match, events is not 0/1, a column of check_matrix does not have one or two ones, or two qubits join the same pair of checks.
Notes

Valid for independent errors and one error type at a time. Needs about 1 / p cycles per qubit for a stable estimate (the paper's Eq. 18). A pair of checks whose correlation is below the statistical noise gives a probability near zero.

Examples:

>>> import numpy as np
>>> from dense_evolution.physics.qec import estimate_edge_probabilities_from_detection_events
>>> checks = np.array([[1, 0], [1, 1]])
>>> events = np.array([[1, 1]] * 20 + [[0, 0]] * 80)
>>> p = estimate_edge_probabilities_from_detection_events(checks, events)
>>> bool(p[0] < 0.5 and p[1] < 0.5)
True
Source code in dense_evolution/physics/qec.py 
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
def estimate_edge_probabilities_from_detection_events(check_matrix, events) -> np.ndarray:
    """Error probability of every qubit, from detection events (Spitz et al.).

    Implements the exact inversion of S. T. Spitz, B. Tarasinski,
    C. W. J. Beenakker and T. E. O'Brien, "Adaptive weight estimator for
    quantum error correction in a time-dependent environment",
    arXiv:1712.02360, Eqs. (13) and (16). For a code in which every qubit is
    checked by at most two checks (repetition and surface codes), each qubit is
    an edge between two checks, or between one check and the boundary. A qubit
    shared by checks ``i`` and ``j`` has probability

    ``p = 1/2 - sqrt(1/4 - (<v_i v_j> - <v_i><v_j>) / (1 - 2 <v_i xor v_j>))``

    where ``v`` are the detection events and ``<.>`` the average over cycles. A
    qubit on the boundary of check ``i`` has

    ``p = 1/2 + (<v_i> - 1/2) / prod(1 - 2 p_ij)`` over the other qubits of ``i``.

    The result can be passed as ``weights`` (``-log(p / (1 - p))``) to
    `pymatching_decode`.

    Parameters
    ----------
    check_matrix : array_like of shape (n_checks, n_qubits)
        0/1 matrix, entry ``[c, q] = 1`` when check ``c`` detects an error on
        qubit ``q``. Every column has one or two ones; no two qubits may join
        the same pair of checks.
    events : array_like of shape (n_cycles, n_checks)
        0/1 detection events: ``1`` when a check changed value since the
        previous cycle (the syndrome of one cycle's new errors).

    Returns
    -------
    numpy.ndarray of shape (n_qubits,)
        Estimated error probability of each qubit, clipped to [0, 0.5].

    Raises
    ------
    ValueError
        If the shapes do not match, ``events`` is not 0/1, a column of
        ``check_matrix`` does not have one or two ones, or two qubits join the
        same pair of checks.

    Notes
    -----
    Valid for independent errors and one error type at a time. Needs about
    ``1 / p`` cycles per qubit for a stable estimate (the paper's Eq. 18). A
    pair of checks whose correlation is below the statistical noise gives a
    probability near zero.

    Examples
    --------
    >>> import numpy as np
    >>> from dense_evolution.physics.qec import estimate_edge_probabilities_from_detection_events
    >>> checks = np.array([[1, 0], [1, 1]])
    >>> events = np.array([[1, 1]] * 20 + [[0, 0]] * 80)
    >>> p = estimate_edge_probabilities_from_detection_events(checks, events)
    >>> bool(p[0] < 0.5 and p[1] < 0.5)
    True
    """
    h = np.asarray(check_matrix, dtype=int)
    v = np.asarray(events, dtype=float)
    if h.ndim != 2:
        raise ValueError("check_matrix must be 2-D (n_checks, n_qubits)")
    if v.ndim != 2 or v.shape[1] != h.shape[0]:
        raise ValueError(
            f"events must have shape (n_cycles, {h.shape[0]}), got {v.shape}"
        )
    if v.shape[0] == 0:
        raise ValueError("events needs at least one cycle")
    if not np.isin(v, (0.0, 1.0)).all():
        raise ValueError("events must contain only 0 and 1")

    n_q = h.shape[1]
    ends = []
    for q in range(n_q):
        rows = np.flatnonzero(h[:, q])
        if len(rows) not in (1, 2):
            raise ValueError(
                f"qubit {q} is checked by {len(rows)} checks; need one or two"
            )
        ends.append(tuple(int(r) for r in rows))
    pairs = [e for e in ends if len(e) == 2]
    if len(set(pairs)) != len(pairs):
        raise ValueError("two qubits join the same pair of checks")

    mean = v.mean(axis=0)
    p = np.zeros(n_q)
    for q, e in enumerate(ends):
        if len(e) == 2:
            i, j = e
            cov = (v[:, i] * v[:, j]).mean() - mean[i] * mean[j]
            xor = np.abs(v[:, i] - v[:, j]).mean()
            denom = 1.0 - 2.0 * xor
            inside = 0.25 - cov / denom if denom > 0 else 0.25
            p[q] = 0.5 - np.sqrt(max(inside, 0.0))
    for q, e in enumerate(ends):
        if len(e) == 1:
            i = e[0]
            prod = 1.0
            for r, other in enumerate(ends):
                if len(other) == 2 and i in other:
                    prod *= 1.0 - 2.0 * p[r]
            if prod <= 0:
                p[q] = 0.5
            else:
                p[q] = 0.5 + (mean[i] - 0.5) / prod
    return np.clip(p, 0.0, 0.5)

erasure_ml_decode 

erasure_ml_decode(
    observed_syndrome: tuple,
    heralded_qubits: Sequence[int],
    n_qubits: int,
    stabilizers: Sequence[str],
) -> Optional[str]

Maximum-likelihood decoder for erasures at known locations, for any stabilizer code.

With the erased qubits known, the error is supported on them, and the syndrome becomes a linear system over GF(2) in the X and Z bits of those qubits (Delfosse & Zemor, arXiv:1703.01517; Kuo & Ouyang, arXiv:2411.13509). Solving it by Gaussian elimination costs O(n^3), where erasure_aware_decode tries 4**m assignments for m erased qubits. Any solution is a valid correction when every zero-syndrome operator on the erased qubits is a stabilizer element, so degenerate errors (several errors that differ by a stabilizer) are decoded too, not rejected.

Parameters:

Name Type Description Default
observed_syndrome sequence of int

One bit per stabilizer, same convention as compute_syndrome.

 required
heralded_qubits sequence of int

Indices of the erased qubits.

 required
n_qubits int

Number of physical qubits.

 required
stabilizers sequence of str

Stabilizer generators as Pauli strings of length n_qubits.

 required

Returns:

Type Description
str or None

A Pauli string supported on the erased qubits that reproduces the syndrome and is equivalent, up to a stabilizer, to every other solution. None when no error on the erased qubits explains the syndrome, or when the erased qubits contain a logical operator so the correction is ambiguous. With no erased qubits it returns the identity for a zero syndrome and None otherwise.

Raises:

Type Description
ValueError

If the syndrome length does not match the stabilizers, a stabilizer has the wrong length, or an erased index is out of range.

Examples:

>>> stabs = ['IIIXXXX', 'IXXIIXX', 'XIXIXIX', 'IIIZZZZ', 'IZZIIZZ', 'ZIZIZIZ']
>>> syndrome = compute_syndrome('XIIIIIX', stabs)
>>> erasure_ml_decode(syndrome, [0, 6], 7, stabs)
'XIIIIIX'
Source code in dense_evolution/physics/qec.py 
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
def erasure_ml_decode(
    observed_syndrome: tuple,
    heralded_qubits: Sequence[int],
    n_qubits: int,
    stabilizers: Sequence[str],
) -> Optional[str]:
    """Maximum-likelihood decoder for erasures at known locations, for any
    stabilizer code.

    With the erased qubits known, the error is supported on them, and the
    syndrome becomes a linear system over GF(2) in the X and Z bits of those
    qubits (Delfosse & Zemor, arXiv:1703.01517; Kuo & Ouyang,
    arXiv:2411.13509). Solving it by Gaussian elimination costs O(n^3), where
    `erasure_aware_decode` tries 4**m assignments for m erased qubits. Any
    solution is a valid correction when every zero-syndrome operator on the
    erased qubits is a stabilizer element, so degenerate errors (several
    errors that differ by a stabilizer) are decoded too, not rejected.

    Parameters
    ----------
    observed_syndrome : sequence of int
        One bit per stabilizer, same convention as `compute_syndrome`.
    heralded_qubits : sequence of int
        Indices of the erased qubits.
    n_qubits : int
        Number of physical qubits.
    stabilizers : sequence of str
        Stabilizer generators as Pauli strings of length `n_qubits`.

    Returns
    -------
    str or None
        A Pauli string supported on the erased qubits that reproduces the
        syndrome and is equivalent, up to a stabilizer, to every other
        solution. `None` when no error on the erased qubits explains the
        syndrome, or when the erased qubits contain a logical operator so the
        correction is ambiguous. With no erased qubits it returns the identity
        for a zero syndrome and `None` otherwise.

    Raises
    ------
    ValueError
        If the syndrome length does not match the stabilizers, a stabilizer
        has the wrong length, or an erased index is out of range.

    Examples
    --------
    >>> stabs = ['IIIXXXX', 'IXXIIXX', 'XIXIXIX', 'IIIZZZZ', 'IZZIIZZ', 'ZIZIZIZ']
    >>> syndrome = compute_syndrome('XIIIIIX', stabs)
    >>> erasure_ml_decode(syndrome, [0, 6], 7, stabs)
    'XIIIIIX'
    """
    stabs = list(stabilizers)
    if len(observed_syndrome) != len(stabs):
        raise ValueError(
            f"observed_syndrome has {len(observed_syndrome)} entries but there are "
            f"{len(stabs)} stabilizers"
        )
    for i, s in enumerate(stabs):
        if len(s) != n_qubits:
            raise ValueError(f"stabilizers[{i}] has length {len(s)}, expected n_qubits={n_qubits}")
    erased = sorted({int(q) for q in heralded_qubits})
    if any(q < 0 or q >= n_qubits for q in erased):
        raise ValueError(f"heralded_qubits must be in range(0, {n_qubits})")

    syn = np.array(observed_syndrome, dtype=np.uint8) % 2
    k = len(erased)
    if k == 0:
        return 'I' * n_qubits if not syn.any() else None

    sym = np.array([_pauli_to_symplectic(s) for s in stabs], dtype=np.uint8)
    sx, sz = sym[:, :n_qubits], sym[:, n_qubits:]
    a = np.concatenate([sz[:, erased], sx[:, erased]], axis=1)
    rref, pivots = _gf2_rref(np.concatenate([a, syn[:, None]], axis=1))
    if 2 * k in pivots:
        return None

    sol = np.zeros(2 * k, dtype=np.uint8)
    for row, col in enumerate(pivots):
        sol[col] = rref[row, 2 * k]
    free = [c for c in range(2 * k) if c not in pivots]
    stab_rref, stab_pivots = _gf2_rref(sym)

    def embed(v):
        full = np.zeros(2 * n_qubits, dtype=np.uint8)
        for j, q in enumerate(erased):
            full[q] = v[j]
            full[n_qubits + q] = v[k + j]
        return full

    for f in free:
        vec = np.zeros(2 * k, dtype=np.uint8)
        vec[f] = 1
        for row, col in enumerate(pivots):
            if rref[row, f]:
                vec[col] = 1
        if not _in_gf2_span(embed(vec), stab_rref, stab_pivots):
            return None

    full = embed(sol)
    letters = {(0, 0): 'I', (1, 0): 'X', (0, 1): 'Z', (1, 1): 'Y'}
    return ''.join(letters[(int(full[q]), int(full[n_qubits + q]))] for q in range(n_qubits))

peeling_decode 

peeling_decode(
    stabilizers,
    observed_syndrome,
    heralded_qubits,
    n_qubits,
) -> Optional[str]

Peeling decoder for a CSS code whose X-type and Z-type checks each join every qubit to at most two checks (repetition and surface codes).

Source code in dense_evolution/physics/qec.py 
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
def peeling_decode(stabilizers, observed_syndrome, heralded_qubits, n_qubits) -> Optional[str]:
    """Peeling decoder for a CSS code whose X-type and Z-type checks each join
    every qubit to at most two checks (repetition and surface codes)."""
    stabs = list(stabilizers)
    syn = np.array(observed_syndrome, dtype=np.uint8)
    erased = sorted({int(q) for q in heralded_qubits})
    zi = [i for i, s in enumerate(stabs) if 'Z' in s and 'X' not in s]
    xi = [i for i, s in enumerate(stabs) if 'X' in s and 'Z' not in s]
    hz = np.array([[c != 'I' for c in stabs[i]] for i in zi], dtype=np.uint8)
    hx = np.array([[c != 'I' for c in stabs[i]] for i in xi], dtype=np.uint8)
    x_part = _peel(hz, erased, syn[zi])
    z_part = _peel(hx, erased, syn[xi])
    if x_part is None or z_part is None:
        return None
    return ''.join('IXZY'[(q in x_part) + 2 * (q in z_part)] for q in range(n_qubits))

union_find_decode 

union_find_decode(
    stabilizers,
    observed_syndrome,
    heralded_qubits,
    n_qubits,
) -> Optional[str]

Union-Find decoder with erasures and Pauli errors (Delfosse and Nickerson, arXiv:1709.06218): clusters start from the erased qubits and grow by half-edges until every cluster has even syndrome parity or touches the boundary, then each grown cluster is peeled.

Source code in dense_evolution/physics/qec.py 
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
def union_find_decode(stabilizers, observed_syndrome, heralded_qubits, n_qubits) -> Optional[str]:
    """Union-Find decoder with erasures and Pauli errors (Delfosse and
    Nickerson, arXiv:1709.06218): clusters start from the erased qubits and
    grow by half-edges until every cluster has even syndrome parity or touches
    the boundary, then each grown cluster is peeled."""
    hz, sz, hx, sx = _css_split(stabilizers, observed_syndrome)
    erased = sorted({int(q) for q in heralded_qubits})
    parts = []
    for h, s in ((hz, sz), (hx, sx)):
        grown = _uf_grow(h, erased, s)
        if grown is None:
            return None
        part = _peel(h, grown, s)
        if part is None:
            return None
        parts.append(part)
    return ''.join('IXZY'[(q in parts[0]) + 2 * (q in parts[1])] for q in range(n_qubits))

matching_erasure_decode 

matching_erasure_decode(
    stabilizers,
    observed_syndrome,
    heralded_qubits,
    n_qubits,
) -> str

Minimum-weight perfect matching with weight 0 on the erased qubits (Stace, Barrett and Doherty, arXiv:0904.3556), via pymatching.

Source code in dense_evolution/physics/qec.py 
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
def matching_erasure_decode(stabilizers, observed_syndrome, heralded_qubits, n_qubits) -> str:
    """Minimum-weight perfect matching with weight 0 on the erased qubits
    (Stace, Barrett and Doherty, arXiv:0904.3556), via pymatching."""
    import pymatching

    hz, sz, hx, sx = _css_split(stabilizers, observed_syndrome)
    w = np.ones(n_qubits)
    w[list(heralded_qubits)] = 0.0
    xp = pymatching.Matching.from_check_matrix(hz, weights=w).decode(sz)
    zp = pymatching.Matching.from_check_matrix(hx, weights=w).decode(sx)
    return ''.join('IXZY'[int(xp[q]) + 2 * int(zp[q])] for q in range(n_qubits))