Distance finding and fault-tolerance tooling: Rust engine, typed inputs, DEM fault distance, docs

[ciaranra]

Distance-finding and fault-tolerance work, consolidated into one PR (previously split as #439 and #440).

1. Expose the Rust distance search and verification workflow

The incremental-weight distance search in pecos-qec and the StabilizerCodeSpec verification machinery were fully implemented in Rust but had no Python bindings and no callers. This exposes them as a supported replacement for the legacy pecos.analysis.VerifyStabilizers development loop.

• StabilizerCodeSpec with a builder: check, logical_z, logical_x, and three build modes — build(), build_verified() (errors name the exact anticommuting generator pair), and build_with_discovered_logicals() (derives paired logicals and destabilizers by stabilizer simulation).
• distance(max_weight=None, css=False, verbose=False) returning DistanceResult; min_weight_logicals(); shortest_logicals(delta) for logicals within delta of the minimum. The search grows error weight from 1, so cost scales with the distance rather than the qubit count, reaching codes the coset enumeration in StabilizerCode.distance() (capped at k + rank <= 30) cannot.
• Xs, Ys, Zs multi-qubit Pauli helpers, so checks read as Zs([0, 1]) * Y(2).
• Typed matrix input: ParityCheckMatrix (pecos-qec, role-neutral) and SymplecticMatrix (pecos-quantum, [X block | Z block]), with CSS orthogonality validated (Hx * Hz^T = 0, errors naming the offending row pair) and width-bearing zeros constructors for single-stabilizer-type codes. Phase-dropping conversions are named honestly (to_positive_paulis, from_pauli_sequence_ignoring_phase) because symplectic form carries no sign.
• Invariant fix: all three StabilizerCodeSpec constructors now reject linearly dependent stabilizers (DependentStabilizers { rank, count }). Previously num_logical_qubits() returned n - stabilizers.len() while documenting "independent generators", so redundant generators silently corrupted k — which matrix input makes easy to hit.
• Removed pecos/tools/fault_tolerance_checks.py and pecos/tools/stabilizer_verification.py, byte-identical dead copies unreachable through the public path (pecos.tools is a deprecation shim re-exporting pecos.analysis).

pecos.analysis.VerifyStabilizers itself is untouched; retiring it is a follow-up now that every capability has a Rust-backed home.

2. Consolidate the duplicate searches and parallelize

Two independent implementations of the same weight-increasing search existed. Their predicates were verified equivalent — both test "commutes with every stabilizer generator AND anticommutes with at least one configured logical" — so StabilizerFlipChecker::{has_undetectable_logical, compute_distance} now delegate to the shared engine via a new has_logical_error_at_weight. The checker's existing tests are unchanged and act as the regression guard. The combinations/pauli_product/build_pauli_string helpers are deliberately retained: analyze_weight needs configurable X/Y/Z subsets the shared iterator cannot express.

The per-weight candidate scan now runs on rayon, partitioned over support combinations. Output is bit-identical to serial rather than merely equivalent — reduction is on enumeration index, so the same operator and the same vector order come back, and tests cover both the serial and parallel branch.

PARALLEL_CANDIDATE_THRESHOLD = 65_536 candidates at one weight, derived from a measured sweep, not intuition. Below roughly 22k candidates parallelism loses (forcing the toric [[18, 2, 3]] weight-3 tier, 22,032 candidates, parallel made that search 4.6x slower); above roughly 193k it stops engaging where the time is spent (the color [[17, 1, 5]] search is dominated by its weight-4 tier, 192,780 candidates, and a higher gate erased the speedup entirely).

Benchmark	Serial	Parallel	Effect
five-qubit [[5,1,3]]	7.70 us	8.46 us	9.9% slower
Steane [[7,1,3]]	15.04 us	16.83 us	11.9% slower
color [[17,1,5]]	41.7 ms	8.6 ms	4.8x faster
shortest_logicals delta=1, color [[17,1,5]]	2.77 s	215 ms	12.9x faster

This is a trade, not a free win: microsecond-scale searches pay about 10%, while searches long enough to wait on improve 5-13x. Small-code figures were confirmed by an A/B/A run after an initial measurement proved to be machine drift. Adds benches/modules/code_distance.rs; no distance benchmark existed before.

3. Detector-error-model fault distance

Code distance is not circuit distance. A distance-5 code whose syndrome extraction spreads one fault across multiple data qubits can have fault distance 3, so code distance alone can overstate real protection. Nothing in PECOS computed the circuit-level number — check_undetectable_logical_errors enumerates failing configurations but never reports a minimum.

This adds the DEM level: the minimum number of fault mechanisms whose XOR flips no detector but flips at least one observable. With H the detector-by-mechanism matrix and L the observable-by-mechanism matrix, that is minimum |e| with H*e = 0 and L*e != 0 — structurally the same problem as code distance, which is why it lands beside the existing fault-tolerance checkers rather than in a separate silo.

• graphlike_fault_distance is exact when every mechanism flips at most two detectors, searching the parity-doubled graph with every detector AND the boundary as a BFS root. Rooting only at the boundary is not exact: a DEM whose minimum cycle avoids the boundary (D0 D1 L0 / D1 D2 / D0 D2, distance 3) leaves the boundary isolated and finds nothing. That is now a regression test.
• exhaustive_fault_distance(max_weight) is correct for any DEM including hyperedges; max_weight is required because the cost is combinatorial in the mechanism count.
• Hyperedges make the graphlike method fail fast with a count rather than being ignored. This is deliberate: DemMatchingGraph silently skips hyperedges, so building on it would have returned quietly wrong distances. The implementation reads to_mechanisms() directly and reuses FaultMechanism::{xor, is_graphlike, is_hyperedge}.
• Both methods return the witnessing mechanism set, mirroring DistanceResult::min_weight_operator — knowing which faults conspire is the point.

Guarded by a seeded property test over 512 random small graphlike DEMs asserting both methods agree on distance and on solution existence. Fixture tests alone shared a blind spot with the original design (every case happened to have a boundary edge); cross-validating an exact special case against a general reference catches the class rather than one instance. Verified by mutation: restricting the roots to boundary-only fails both the boundary-free regression and the property test.

4. Documentation

New docs/user-guide/stabilizer-code-verification.md, replacing the legacy stab_code_verification.rst narrative on the current API: the builder workflow, the ten-qubit design storyline (anticommuting pair diagnosed, then [[10, 3]] at distance 2, then [[10, 1]] at distance 3), logical-operator exploration, matrix input with its error diagnostics, and a note on choosing between the two distance methods. Ten executable doc tests, no skip markers, wired into the mkdocs nav.

Verification

Run on the combined branch:

• cargo test -p pecos-qec -p pecos-quantum: no failures
• uv run --frozen pytest across the stabilizer-code binding suites, the fault-distance suite, and the generated doc tests: 40 passed
• just build-debug and just lint (new files staged first, since pre-commit only inspects tracked files): clean

[ciaranra]

Two fault-tolerance correctness fixes have been folded in (previously #441 and a follow-on), alongside the distance tooling.

Multi-fault propagation

propagate_faults XORed every fault into the initial PauliProp and propagated from min_tick, so a fault at tick 5 was injected as though it existed earlier and intervening gates acted on a Pauli that should not yet have existed. Any weight >= 2 result was untrustworthy. It also ignored before entirely, so it disagreed with propagate_fault even at weight 1 for a before=false fault whose own tick contains a gate on its qubit.

Faults are now injected at their own tick and before/after position. The duplicated Pauli-injection mapping that let the two functions drift is now a single shared helper.

Single-leg fault enumeration

PauliFaultIterator assigned a non-identity Pauli to every qubit of a location, so at a two-qubit gate it generated only 9 of the 15 non-identity two-qubit Paulis — IX, XI, IY, YI, IZ, ZI were structurally unreachable. Each leg now chooses from identity plus the enabled Paulis, with identity-only locations rejected. The weight convention is unchanged: weight counts locations, since a two-qubit gate failing is one fault however many qubits it corrupts. pauli_types() keeps its existing public meaning; identity is handled inside the iterator.

These two bugs were masking each other. Three tests described injecting a single data-qubit X, the iterator actually produced XX, and the buggy propagation pushed XX through the CX a second time, cancelling one leg and accidentally reproducing the intended effect. Those tests now construct the single-leg fault directly and keep their original assertions, including that naive three-qubit syndrome extraction is not 1-fault tolerant.

The DAG path was already correct (possible_faults offers the identity option per qubit), so the DEM builder and the fault-distance work in this PR were unaffected.

What the enumeration fix surfaced

The omission erred toward false confidence: faults that are never enumerated cannot be found to break a circuit.

test_is_fault_tolerant_method previously reported its circuit (CX(0,1) then MZ(1)) as 1-fault tolerant for X errors. It is not. The newly reachable fault is XI after the tick-0 CX: an X on data qubit 0 alone, which the ancilla measurement never sees, leaving an undetected logical error. The test computes the verdict and only prints it, so it passed either way — its own comments enumerate XX and X-on-qubit-1 but never X-on-qubit-0-alone, matching the enumerator's blind spot.

Conversely test_repeated_syndrome_measurement_concept documents in prose exactly this fault class ("X error on data qubit AFTER its CX gate -> no syndrome in this round"); its undetectable count moves 0 -> 3, so the fix makes that description true.

No test containing an actual fault-tolerance assertion fails. Reported counts move widely, as expected when more faults are tested — for example test_fault_checker_three_qubit_code 54 -> 90 configurations, test_steane_code_fault_enumeration 16 -> 64, and the gadget-checker suites roughly triple. The full old/new delta list is available on request.

Two follow-ups worth separate attention, deliberately not changed here:

• Several diagnostics compute a fault-tolerance verdict without asserting it, so they cannot fail. test_is_fault_tolerant_method is the clearest case.
• Two tests now conflict with their own prose: test_syndrome_detection_three_qubit_code says "should be 0" and reports 3, and test_analyze_with_follow_up_resolves_ambiguity shows ambiguity rising from 3 to 15 with follow-up.

Verification

• cargo test -p pecos-qec: 829 passed, no failures
• uv run --frozen pytest python/quantum-pecos/tests/qec -q: 1140 passed, 1 skipped, 1 xfailed
• just build-debug, just lint: clean
• Both fixes mutation-verified: reverting the before/after injection order fails the after-tick equivalence and property tests; removing identity from the per-leg choices fails all four single-leg enumeration tests.

[ciaranra]

Hook-error diagnosis added (crates/pecos-qec/src/fault_tolerance/hook_errors.rs).

What it reports

PauliPropChecker::diagnose_hook_errors(data_qubits, z_ancillas, x_ancillas, logicals, min_data_weight) returns, for each amplifying fault: the responsible gate (SpacetimeLocation — tick, gate type, qubits, gate index), the injected single-qubit fault, the resulting error support restricted to the data block, whether it is detected, and whether it causes a logical error.

A hook error is defined as a fault whose OWN Pauli weight is exactly 1 but whose propagated support on the data qubits has weight at least min_data_weight. Both halves matter: a weight-2 fault landing as a weight-2 data error is not amplification and is deliberately not reported, which is what distinguishes this from a plain output-weight filter. min_data_weight is explicit with no default; 2 is the standard threshold.

detected and causes_logical_error are carried because an amplified error that still trips a syndrome does not reduce distance — only an undetected one does. Without that distinction the report would be a list of alarming faults with no way to tell which ones cost you anything.

This is why circuit fault distance falls below code distance, so the point of the diagnostic is attribution: not "your fault distance is 3 rather than 5", but which gate makes it so.

Design notes

It is a readout over existing machinery, reusing analyze_all_faults, has_syndrome, and anticommutes_with_logical; no new propagator or enumerator. data_qubits is caller-supplied rather than guessed. Output is sorted by tick, gate index, qubits, then Paulis so results are reproducible.

FaultChecker::check_output_weight_expansion is left untouched. It flags configurations exceeding an output weight but returns only the offending configurations — no resulting support, no amplification test, no gate attribution.

This diagnosis is only meaningful because of the single-leg enumeration fix in this PR: an ancilla-only fault on a CX was previously unreachable, so the analysis would have found nothing and looked correct doing it.

Rust-only for now. PauliPropChecker is not exposed to Python, so bindings are a follow-up rather than something bolted on here.

Verification

• cargo test -p pecos-qec: 835 passed, 0 failed (669 unit, 108 integration, 58 doctests)
• just build-debug, just lint: clean, with the new file staged so pre-commit actually inspects it
• Mutation-verified: relaxing the own-weight-equals-one condition fails weight_two_fault_with_weight_two_data_support_is_not_a_hook; replacing the data-qubit restriction with all propagated qubits fails the amplification test on its expected support.

[ciaranra]

Circuit fault distance and DEM search pruning added.

Circuit fault distance as a number, per logical

FaultChecker::check_undetectable_logical_errors enumerates failing configurations but never minimises — create_fault_iterator builds an iterator at exactly config.max_weight and does not scan 1..=max_weight, so the API could answer "does it fail at weight w" but never "what is the smallest w". A single aggregate also hides which logical is weakest.

Added circuit_fault_distance(...) returning CircuitDistanceResult { distance, witness, logical_index }, and per_logical_circuit_fault_distances(...) returning one entry per supplied logical. Both share one increasing-weight loop differing only in the stopping rule, and both take an explicit max_weight because the enumeration is combinatorial. The overall distance equals the minimum over the per-logical values, asserted on a case where they genuinely differ ([1, 2]).

The existing single-weight methods are untouched — they are shipped API with dependents.

Mutation-verified: collapsing the weight loop to a single weight turns the discrimination result from [Some(1), Some(2)] into [Some(2), Some(2)]; stopping the per-logical search on first hit gives [Some(1), None]; inverting the per-logical bit gives [Some(1), Some(1)].

DEM cross-validation was investigated and deliberately NOT added: the two fault models are not directly comparable. FaultChecker creates one location per TickCircuit gate batch spanning all qubits in that batch, whereas the DAG path splits locations per qubit and reconstructs two-qubit noise mechanisms separately; DEM prep/measurement faults are noise-channel-specific while this API enumerates a configured Pauli set; and DEM logical outputs come from measurement metadata rather than anticommutation with supplied final Pauli logicals. A synthesised DEM would have papered over those differences and produced a test that resembled validation without being it.

Connected-cluster pruning for the DEM search

exhaustive_fault_distance enumerates blind k-subsets, which is correct but unusable at real scale — a distance-3 surface memory DEM already has ~1300 contributions at 3 rounds and ~2050 at 5 rounds, and the literature reports ~43,000 mechanisms for an 11-round surface circuit.

Added connected_cluster_fault_distance(dem, max_weight), exact for any DEM including hyperedges, using two provable prunes:

• Connectivity. A minimum-weight undetectable observable-flipping set is connected in the shared-detector graph. If it split into components, each would be individually detector-free (no detector spans components and each appears an even number of times), observable parity XORs across components, so some component alone would be a strictly smaller solution. So clusters grow outward from a seed rather than enumerating arbitrary subsets. This is the published Connected Cluster approach (arXiv:2603.22532), credited as such in the module docs.
• Unique-detector peeling. A detector appearing in exactly one mechanism means that mechanism can never belong to an undetectable set, since the detector would flip an odd number of times. Removal can make further detectors unique, so it iterates to a fixpoint.

exhaustive_fault_distance is retained deliberately as the simple reference implementation used to validate the pruned one.

Measured on a 594-mechanism DEM with non-peelable cycle padding, so the numbers isolate the connectivity gain rather than peeling collapsing the input:

Search	Weight 3	Result
blind exhaustive_fault_distance	25-31 ms	3
connected_cluster_fault_distance	0.8 ms	3

About 38x, same answer. At weight 4 the pruned search completes in 0.8 ms; the blind search would have to consider 5.13e9 candidate subsets.

A bug this work exposed in the existing code

The extended property test caught an inconsistency introduced earlier in this PR. graphlike_fault_distance handled the weight-1 detector-free case BEFORE checking for hyperedges, so a DEM containing both a hyperedge and a detector-free observable-flipping mechanism returned Ok(Some(1)) instead of the documented hyperedge error — the distance was right, but the function sometimes refused hyperedge DEMs and sometimes answered them, depending on whether such a mechanism happened to exist. The hyperedge check is now unconditional and first. Callers wanting an answer regardless have the two methods that handle hyperedges.

Only the randomised generator surfaced this: it needs a DEM with a hyperedge AND a weight-1 detector-free mechanism together, which was case 27 of 512 and which no hand-written fixture had produced.

Tests

The seeded property test now generates hyperedge DEMs as well as graphlike ones, and asserts the blind and pruned searches agree on distance and solution existence for every case, with the graphlike method additionally agreeing where applicable. Plus fixtures for peeling reaching a fixpoint without changing the distance, and peeling preserving a witness whose detectors are all shared.

Mutation-verified: making peeling over-prune (deactivating a mechanism whose detectors are all shared) fails four tests including the property test and both distance-3 witness fixtures.

Verification

• cargo test -p pecos-qec: no failures
• cargo clippy --locked -p pecos-qec --all-targets -- -D warnings: clean
• just build-debug, just lint (files staged so pre-commit inspects them): clean

[ciaranra]

Flag fault-tolerance verification added (crates/pecos-qec/src/fault_tolerance/flag_verification.rs).

What it checks

PauliPropChecker::verify_flag_fault_tolerance(data_qubits, flag_qubits, measured_stabilizer, t) verifies the propagated-fault half of the Chao-Reichardt t-flag condition (arXiv:1708.02246): for every fault configuration of weight v in 1..=t, if no flag qubit is raised then min(wt(E), wt(E * P)) <= v, where E is the propagated data error and P the stabilizer being measured. Violations are returned with the offending configuration, v, and the computed weight.

The min encodes stabilizer equivalence: E and E * P are the same error modulo the stabilizer being measured, so checking wt(E) alone reports violations that are not violations. Worked example from the tests, on the unflagged weight-4 XXXX measurement with an X on the measurement ancilla after CX(a, 0):

E     = X1 X2 X3        wt(E)     = 3
P     = X0 X1 X2 X3
E * P = X0              wt(E * P) = 1

At v = 1, min(3, 1) = 1, so that configuration is correctly not a violation.

Scope limitation, stated rather than implied

The paper's definition also requires that a fault-free run does not flag. That is not checked here, and cannot be: PauliProp is a Pauli-frame model tracking deviations from the ideal execution, so a fault-free run has an empty frame by construction and any flag-outcome field would be permanently false while appearing to verify something. Establishing that half needs stabilizer simulation of the ideal circuit.

The verdict field is therefore named fault_condition_satisfied, not is_t_flag, and both the function and module docs say which half is covered and which must be established separately.

Tests

Six tests, built around a discriminating pair rather than a single circuit: the standard single-flag weight-4 stabilizer measurement satisfies the condition at t = 1, and the same measurement without flag interleaving fails it with a weight-1 fault producing a weight-2 data error. If both circuits returned the same verdict the check would be measuring nothing.

Also: the stabilizer-equivalence case above, restriction of weights to the caller-supplied data qubits, determinism, and a negative t = 2 case. The weight-4 flagged circuit turned out to satisfy the propagated condition at t = 2, so the negative test uses a weight-6 single-flag circuit instead of asserting something false about the weight-4 one.

Mutation-verified:

• replacing min(wt(E), wt(E*P)) with wt(E) fails the stabilizer-equivalence test (computed weight 3 instead of 1) and two others
• inverting flag detection swaps the verdicts of the good and unflagged circuits
• widening the weight computation past the caller's data qubits fails the restriction test

Verification

• cargo test -p pecos-qec: 848 passed, 0 failed (682 unit, 108 integration, 58 doctests)
• cargo clippy --locked -p pecos-qec --all-targets -- -D warnings, cargo fmt --all -- --check: clean
• just lint with files staged: clean

Python exposure is the outstanding follow-up; PauliPropChecker is not currently bound.

[ciaranra]

Exact distance certification for qLDPC-scale codes added, closing the algorithm-import item of the distance roadmap.

Method choice, from evidence rather than reputation

The large-scale empirical study arXiv:2606.12445 was read before choosing: Brouwer-Zimmermann no longer holds its traditional advantage on qLDPC codes; branch-and-bound MaxSAT wins; scalability is governed by cardinality-constraint handling (sequential counter / totalizer), XOR-aware reasoning does not systematically help; exactness comes from an incremental loop where UNSAT at weight w proves d > w and the first SAT weight is d.

DistanceProblem (crates/pecos-qec/src/distance_problem.rs)

One (H, L) GF(2) encoder serves both existing problem shapes: CSS code distance (from ParityCheckMatrix pairs or a CSS StabilizerCodeSpec; non-CSS errors clearly) and DEM fault distance (from to_mechanisms()). Tseitin XOR chains for parity, a Sinz sequential counter for the weight bound, DIMACS and new-format WCNF export with commented variable-range roles.

The trust model is the point of the design: verify_witness checks H e = 0 and L e != 0 natively, so the SAT half of any answer needs no solver trust at all; UNSAT answers (and therefore exactness) rest on the solver, and the docs say so rather than implying both halves are certified.

Tested without any solver: an exhaustive evaluator over the emitted CNF (aux variables are functionally determined, so no search) proves, for every assignment of small instances, satisfiability at bound w iff the native predicate holds — cross-checked against the existing distance searches on Steane and the repetition-triad DEM, with both sides of the sequential-counter boundary pinned and lying-solver mocks rejected. Mutation-verified at every layer, including a Tseitin polarity flip (kills five tests) and counter off-by-one.

In-process backend: certified_distance via batsat

Per the backend decision, a pure-Rust solver: batsat 0.5 (MiniSat 2.2 reimplementation, MIT, sole transitive dependency bit-vec). Fed from the internal clause representation; a fresh deterministic solver per weight. batsat receives no more trust than an external solver — its models pass through verify_witness before being believed, and that guard is mutation-proven: an off-by-one in model extraction is caught by the certification layer as InvalidWitness (OddCheck), not accepted.

Measured capability, not aspiration

Bivariate bicycle codes built in-test from their polynomial definitions, sanity-checked (n, CSS orthogonality, k via F2Matrix rank) before timing:

Code	Result	Total time
BB [[72,12,6]]	d = 6 certified, witness verified	~0.4 s
BB [[144,12,12]] (the gross code)	d = 12 certified, witness verified	~15 min

Per-weight profile on the gross code: UNSAT proofs escalate (w=10: 163 s, w=11: 674 s), then SAT at w=12 in 1.3 s — proving d > 11 is where the time goes. For comparison the study's MaxCDCL does this instance in ~48 s, so the pure-Rust backend is roughly 19x off the state of the art but completes the flagship qLDPC benchmark entirely in-process. The WCNF export is the documented path to external MaxSAT solvers for anyone needing the faster regime; both were measured, neither is guessed.

These are regimes the existing weight search cannot touch (C(72,6) * 3^6 alone is ~1e11 candidates).

Verification

• cargo test -p pecos-qec: 697 unit + 108 integration + 58 doctests, 0 failures (BB probes are #[ignore]d timing tests, run separately)
• cargo clippy --locked -p pecos-qec --all-targets -- -D warnings, cargo fmt --all -- --check: clean
• just lint including the dependency-integrity gate over the new lockfile entries: clean

Python bindings for DistanceProblem/certified_distance are the noted follow-up.