ScaLERQEC

Figure 1: Our logo.

ScaLERQEC is a scalable framework for estimating logical error rates (LER) of quantum error-correcting (QEC) circuits at scale. It combines an optimized C++ backend (QEPG with SIMD acceleration and OpenMP parallelism) with high-level Python interfaces for QEC experimentation, benchmarking, symbolic analysis, and Monte Carlo fault injection.

ScaLERQEC is compatible with Stim circuits, but uses a completely different approach -- stratified fault sampling with S-curve fitting -- to estimate logical error rates orders of magnitude faster than brute-force Monte Carlo.

Citation

If you use ScaLERQEC in your research, please cite our paper:

@misc{ye2026scalabletestingquantumerror,
      title={Scalable testing of quantum error correction},
      author={John Zhuoyang Ye and Jens Palsberg},
      year={2026},
      eprint={2602.04921},
      archivePrefix={arXiv},
      primaryClass={quant-ph},
      url={https://arxiv.org/abs/2602.04921},
}

Documentation

We use Sphinx to automatically generate the documents:

python -m sphinx.cmd.build -b html docs/source docs

You may visit the current documentation through the following link:

Documentation website: https://yezhuoyang.github.io/ScaLERQEC/

Installation

Option 1 -- Install via pip (recommended)

pip install scalerqec

This installs:

• The Python package scalerqec
• The compiled C++ QEPG backend (scalerqec.qepg) with SIMD and OpenMP support
• All Python modules for LER calculation, sampling, symbolic analysis, etc.

You can then immediately import all modules in Python:

import scalerqec
import scalerqec.qepg

Option 2 -- Install from source

1. Clone the repository:

git clone https://github.com/yezhuoyang/ScaLERQEC.git
cd ScaLERQEC

2. Build and install:

pip install -e .

This compiles the C++ backend using pybind11 and installs the package in development mode.

Prerequisites

All platforms:

• Python >= 3.10
• A C++20-compatible compiler

No external C++ libraries required. The C++ backend is self-contained -- Boost has been removed. All dependencies (pybind11, stim, pymatching, etc.) are handled automatically by pip.

Platform-specific notes:

Project Structure

scalerqec/
├── qepg                # C++ QEPG backend with SIMD acceleration (via pybind11)
├── Clifford/           # Clifford circuit representation, STIM parser, Python QEPG
├── Monte/              # Monte Carlo LER estimation (standard + LDPC codes)
├── QEC/                # High-level QEC circuit construction from stabilizers
├── Stratified/         # ScaLER: stratified sampling with S-curve fitting
├── Symbolic/           # Exact symbolic LER computation (ground truth)
└── util/               # Binomial utilities, Pauli helpers, output formatting

Bundled Stim circuits (under stimprograms/):

Quick Start

Sample once, plot LER versus p

New in 1.1.0: reuse a fixed circuit-and-decoder profile across physical error probabilities. The expensive sampling and fitting run only once:

import numpy as np
from scalerqec.Stratified import Scaler, LERProfile

profile = Scaler(time_budget=60).profile_from_file(
    "stimprograms/surface/surface3", codedistance=3,
    decoder_reference_p=0.001,
)
profile.save("surface3-profile.json")

p_values = np.geomspace(1e-5, 0.02, 100)
logical_error_rates = profile.evaluate(p_values)
ax = profile.plot(p_values)
ax.figure.savefig("ler-vs-p.pdf", bbox_inches="tight")

# Reuse in a later session without sampling or decoding.
profile = LERProfile.load("surface3-profile.json")
print(profile.evaluate(0.0005))

This interface requires a noiseless Stim circuit, one logical observable, and the paper's uniform independent single-qubit depolarizing (SID) model. The decoder is fixed at decoder_reference_p; it is not rebuilt for each p. codedistance must be the circuit-level distance for that circuit and decoder. The S-curve still has systematic extrapolation error. In particular, evaluating a broad p range does not establish accuracy throughout that range.

Use profile.curve(p_values) to inspect contributions from fitted weights, and profile.sample_counts / failure_counts to inspect the measured evidence. These diagnostics are not confidence intervals. See the profile guide for the precise assumptions and API. After an existing calculate_LER_from_file() run, use scaler.get_profile().

1. Define a QEC code with StabIR

StabIR is our stabilizer-level intermediate representation. You define the code structure (stabilizers, logical operators, measurement scheme, rounds) independently of any noise model. The IR is then compiled into a noiseless Stim circuit.

A detailed tutorial is available in Tutorial.ipynb. Below is a smaller example using the [[3, 1, 3]] Z-repetition code.

from scalerqec.QEC.qeccircuit import StabCode

# Step 1: Define the code structure
qeccirc = StabCode(n=3, k=1, d=3)
qeccirc.add_stab("ZZI")          # stabilizer generators
qeccirc.add_stab("IZZ")
qeccirc.set_logical_Z(0, "ZZZ")  # logical Z operator

# Step 2: Configure the measurement scheme
qeccirc.scheme = "Standard"       # also supports: Shor, Knill, Flag
qeccirc.rounds = 2

# Step 3: Compile to a noiseless Stim circuit
qeccirc.construct_circuit()

# The compiled circuit is available as a stim.Circuit object
print(qeccirc.stimcirc)

You can inspect the intermediate representation:

qeccirc.show_IR()

Output:

c0 = Prop[r=0, s=0] ZZI
c1 = Prop[r=0, s=1] IZZ
c2 = Prop[r=1, s=0] ZZI
d0 = Parity c0 c2
c3 = Prop[r=1, s=1] IZZ
d1 = Parity c1 c3
c4 = Prop ZZZ
o0 = Parity c4

Once compiled, you can combine the code with any noise model to estimate the logical error rate. The noise model is applied separately -- it is not part of the code definition:

from scalerqec.QEC.noisemodel import NoiseModel, SD6NoiseModel, SI1000NoiseModel

# Option A: Simple depolarizing noise
noise_model = NoiseModel(0.001)

# Option B: Standard depolarizing (6 noise locations per round)
noise_model = SD6NoiseModel(0.001)

# Option C: Superconducting-inspired noise
noise_model = SI1000NoiseModel(0.001)

# Estimate LER with Monte Carlo
from scalerqec.Monte import MonteLERcalc

mc = MonteLERcalc(time_budget=30, samplebudget=500000, MIN_NUM_LE_EVENT=50)
mc.calculate_LER_from_StabCode(qeccirc, noise_model)
print(f"LER = {mc._estimated_LER:.2e}")

All LER calculators (ScaLER, Monte Carlo, Symbolic) accept a StabCode + NoiseModel pair. You can also pass a custom decoder via decoder=my_decoder -- any object with a decode_batch() method works.

2. ScaLER -- Time-budgeted S-curve LER estimation (main method)

ScaLERQEC estimates the LER by stratified fault sampling and S-curve fitting. The Scaler class runs a three-phase adaptive algorithm within a wall-clock time budget:

1. Phase 1 (Initialization): Binary search for error-onset and saturation weights, fit an initial S-curve.
2. Phase 2 (Sweet-spot exploration): Uniformly sample weights near the theoretical sweet spot.
3. Phase 3 (Adaptive refinement): Iteratively add samples at weights needing more logical error events.

Figure 2: Diagram for the main method in ScaLERQEC.

From a Stim circuit file (recommended):

from scalerqec.Stratified import Scaler, ModelType

scaler = Scaler(
    error_rate=0.001,
    time_budget=120,               # wall-clock seconds
    model_type=ModelType.OUR_MODEL,
    gamma=1,
    num_subspaces_phase2=12,
)

ler = scaler.calculate_LER_from_file(
    filepath="stimprograms/surface/surface7",
    pvalue=0.001,
    codedistance=7,
    figname="Figures/Surface7_ScaLER",
    titlename="Surface Code d=7",
)

print(f"LER: {ler:.6e}")
print(f"R²: {scaler._R_square_score:.4f}")
print(f"Sweet spot: {scaler._sweet_spot}")
print(f"Total samples: {sum(scaler._subspace_sample_used.values()):,}")

Output:

Figure 3: Subspace error rate in log-logit space	Figure 4: Fitted S-curve in probability space

3. ScaLER for LDPC codes

For LDPC codes (e.g., BB codes), use ScalerLDPC which integrates a belief-propagation + OSD decoder:

from scalerqec.Stratified.ScalerLDPC import ScalerLDPC

calculator = ScalerLDPC(
    error_rate=0.001,
    time_budget=60,          # seconds
    max_bp_iters=20,
    osd_order=0,
)
calculator.calculate_LER_from_file(
    filepath="stimprograms/ldpc/bbcode_72_12_6_rounds18",
    pvalue=0.001,
    codedistance=6,
    figname="BBCode",
    titlename="BB Code [[72,12,6]]",
)

4. Monte Carlo LER estimation

Standard Monte Carlo fault injection with adaptive batching. All calculators accept an optional decoder parameter -- any object with a decode_batch() method (pymatching, BPOSD, or your own). If omitted, pymatching is used by default. See Section 1 above for how to use calculate_LER_from_StabCode with a StabCode + NoiseModel pair.

From a Stim circuit file (uniform noise):

mc = MonteLERcalc(time_budget=30, samplebudget=500000)
ler = mc.calculate_LER_from_file(
    samplebudget=500000,
    filepath="stimprograms/surface/surface7",
    pvalue=0.001,
)
print(f"LER = {ler:.2e}")

From a Stim circuit string (non-uniform noise):

ScaLERQEC supports circuits with mixed noise types (DEPOLARIZE1 at varying rates, X_ERROR, Y_ERROR, Z_ERROR, DEPOLARIZE2):

mc = MonteLERcalc(time_budget=30, samplebudget=500000)
stim_str = open("my_noisy_circuit.stim").read()
ler = mc.calculate_LER_from_stim_circuit(stim_str)
print(f"LER = {ler:.2e}")

Monte Carlo for LDPC codes:

from scalerqec.Monte import MonteLDPC

mc_ldpc = MonteLDPC(time_budget=60, samplebudget=100000, max_bp_iters=20)
ler = mc_ldpc.calculate_LER_from_file(
    samplebudget=100000,
    filepath="stimprograms/ldpc/bbcode_72_12_6_rounds18",
    pvalue=0.001,
)

5. Symbolic LER analysis (exact ground truth)

ScaLERQEC can compute the exact symbolic polynomial representation of the logical error rate:

from scalerqec.Symbolic import SymbolicLERcalc

sym = SymbolicLERcalc()
exact_ler = sym.calculate_LER_from_StabCode(qeccirc, noise_model)
print(f"Exact LER = {exact_ler:.6e}")

# Or from a Stim file
exact_ler = sym.calculate_LER_from_file(
    filepath="stimprograms/surface/surface3",
    pvalue=0.001,
)

This is useful for validating Monte Carlo and ScaLER estimates on small circuits.

6. Hotspot analysis -- identify dominant error sources

ScaLERQEC includes a decoder-agnostic hotspot analysis module that reveals which noise categories contribute most to logical failures. It uses a three-phase C++-accelerated pipeline:

1. C++ labeled sampling: sample noise via QEPG with per-shot category bitmask tracking (near-zero overhead).
2. User's decoder: decode detector outcomes with any decoder (pymatching, BPOSD, custom).
3. C++ aggregation: compute P(category fired | logical error), lift, and multi-error configuration breakdown.

The noise label system classifies each DEPOLARIZE1 source into one of four QStab IR error types using positional classification (relative to the CX schedule):

import pymatching
from scalerqec.Analysis.hotspot import HotspotAnalyzer
from scalerqec import qepg as qepg_cpp

# Build circuit and decoder (any decoder with decode_batch works)
stim_circuit = code.stimcirc
dem = stim_circuit.detector_error_model(decompose_errors=True)
matcher = pymatching.Matching.from_detector_error_model(dem)

# Compile QEPG and auto-label noise sources
prog_str = code.circuit.stim_str
graph = qepg_cpp.compile_QEPG(rewrite_stim_code(prog_str, keep_noise=True))
cc_cpp = qepg_cpp.CliffordCircuit()
cc_cpp.compile_from_rewrited_stim_string(rewrite_stim_code(prog_str))
label_map = qepg_cpp.auto_label(cc_cpp, num_data_qubits)

# Run hotspot analysis with your chosen decoder
analyzer = HotspotAnalyzer(graph, label_map, decoder=matcher)
result = analyzer.analyze(noise_probs, num_shots=1_000_000)
analyzer.print_report(result)

Output:

======================================================================
  Hotspot Analysis  (613 logical errors, LER = 0.000613, 1000000 shots)
======================================================================
  Category                  Count   P(fire)  P(fire|err)  P(err|fire)    Lift
  ------------------------- -----  --------  -----------  -----------  ------
  data_qubit_error              4    0.0200       0.3622       0.0111   18.10
  ghost_error                  11    0.0534       0.9494       0.0109   17.77
  measurement_error             2    0.0101       0.0326       0.0020    3.23
  hook_error                   16    0.0772       0.1974       0.0016    2.56
======================================================================

7. Using the C++ QEPG backend directly

The QEPG (Quantum Error Propagation Graph) is a binary model of how errors propagate to flip detector outcomes.

Figure 5: Illustration of how we compile a QEPG graph in ScaLERQEC.

import scalerqec.qepg as qepg

# Compile a Stim circuit into a reusable QEPG graph
stim_str = open("stimprograms/surface/surface7").read()
graph = qepg.compile_QEPG(stim_str)

# Sample at fixed error weight (stratified sampling)
det_outcomes, obs_outcomes = qepg.return_samples_many_weights_separate_obs_with_QEPG(
    graph,
    weight=[3, 5, 7],
    shots=[10000, 10000, 10000],
)

# Monte Carlo sampling at a given error rate
det, obs = qepg.return_samples_Monte_separate_obs_with_QEPG(
    graph, error_rate=0.001, shot=100000
)

# Non-uniform noise sampling (per-source probabilities)
import numpy as np
from scalerqec.Monte.noise_model_parser import extract_noise_model

noise_model = extract_noise_model(stim_str)
det, obs = qepg.return_samples_nonuniform_to_numpy(
    graph,
    noise_model.noise_probs,
    np.array([p.source_a for p in noise_model.correlated_pairs], dtype=np.int64),
    np.array([p.source_b for p in noise_model.correlated_pairs], dtype=np.int64),
    np.array([p.prob for p in noise_model.correlated_pairs], dtype=np.float64),
    shot=100000,
)

LogiQ -- A high-level, fault-tolerant quantum programming language

We introduce LogiQ -- which supports users to define their own logical QEC block and implement logical Clifford+T operations.

# 1) Define a family of surface codes (sugar -> CSSCode core)
code surface(d: Int) as CellComplex over Z2 {

  cells {
    faces     F[x,y]  in 0..(d-2), 0..(d-2);
    edges_x   Ex[x,y] in 0..(d-2), 0..(d-1);
    edges_y   Ey[x,y] in 0..(d-1), 0..(d-2);
    vertices  V[x,y]  in 0..(d-1), 0..(d-1);
  }

  boundary {
    d2(F[x,y]) =
      Ex[x,y]   +
      Ey[x+1,y] +
      Ex[x,y+1] +
      Ey[x,y];

    d1(Ex[x,y]) = V[x,y]   + V[x+1,y];
    d1(Ey[x,y]) = V[x,y]   + V[x,y+1];
  }

  css {
    hx = matrix(d2);
    hz = transpose(matrix(d1));
  }
}

code five_qubit as StabilizerCode {

  # Number of physical qubits (optional if implied by generator length)
  n = 5;

  generators {
    S0 = "XZZXI";
    S1 = "IXZZX";
    S2 = "XIXZZ";
    S3 = "ZXIXZ";
  }

  logical_z {
    LZ0 = "ZZZZZ";
  }
}


surface q1 [n=40, k=1, d=5]   # First surface-code block (distance-5)
surface q2 [n=40, k=1, d=5]   # Second surface-code block (distance-5)
surface t0 [n=84, k=1, d=7]   # Magic-T ancilla block (distance-7)

q1[0] = LogicH q1[0]

t0 = Distill15to1_T[d=25]     # returns a magic_T handle (see MagicQ below)
InjectT q1[0], t0

q2[1] = LogicCNOT q1[0], q2[1]

c1 = LogicMeasure q1[0]
c2 = LogicMeasure q2[1]

MagicQ -- A high level fault-tolerant quantum programming for dynamic protocol with Post-selection

We introduce MagicQ -- which allows the user to construct a magic state factory. MagicQ also has the full power to express all code-switching protocols.

protocol Distill15to1_T(surface f, int d):
  Repeat:

      # ---- X-type stabilizer checks ----
      c_x1 = LogicProp IIIIIIIXXXXXXXX
      c_x2 = LogicProp IIIXXXXIIIIXXXX
      c_x3 = LogicProp IXXIIXXIIXXIIXX
      c_x4 = LogicProp XIXIXIXIXIXIXIX

      # ---- Z-type stabilizer checks ----
      c_z1  = LogicProp IIIIIIIIZZZZZZZZ
      c_z2  = LogicProp IIIZZZZIIIIZZZZ
      c_z3  = LogicProp IZZIIZZIIZZIIZZ
      c_z4  = LogicProp ZIZIZIZIZIZIZIZ
      c_z12 = LogicProp IIIIIIIIIIZZZZ
      c_z13 = LogicProp IIIIIIIIZZIIIZZ
      c_z14 = LogicProp IIIIIIIIZIZIZIZ
      c_z23 = LogicProp IIIIIZZIIIIIIZZ
      c_z24 = LogicProp IIIIZIZIIIIIZIZ
      c_z34 = LogicProp IIZIIIZIIIZIIIZ

      Success = c_x1 == 0 && c_x2 == 0 && c_x3 == 0 && c_x4 == 0 &&
                c_z1 == 0 && c_z2 == 0 && c_z3 == 0 && c_z4 == 0 &&
                c_z12 == 0 && c_z13 == 0 && c_z14 == 0 &&
                c_z23 == 0 && c_z24 == 0 && c_z34 == 0
      Until Success

      return

Performance: QEPG Sampling vs Stim

ScaLERQEC's C++ QEPG (Quantum Error Propagation Graph) backend provides significantly faster detector sampling than Stim's built-in detector sampler. The QEPG compiles the circuit into a sparse error-propagation graph once, then generates detector and observable outcomes by traversing only the active graph edges — avoiding full tableau simulation on every shot.

Two QEPG sampling modes are available:

• QEPG Monte Carlo: Sparse Poisson error injection (same semantics as Stim's detector sampler). Uses O(k) sparse sampling where k ~ Np is the expected number of faults, instead of iterating over all N noise locations.
• QEPG Fixed-Weight: Injects exactly w faults per shot (used by ScaLER stratified estimation). Cost is O(weight) per sample via Floyd's algorithm.

Both modes use a fused sampling pipeline that writes directly into NumPy buffers with zero intermediate allocation, SIMD-accelerated XOR accumulation, and OpenMP parallelism.

Single-threaded comparison (OMP_NUM_THREADS=1, 100K shots, p=0.001)

Stim's detector sampler is single-threaded with excellent SIMD optimization. For a fair apples-to-apples comparison, we benchmark with OpenMP disabled:

The single-threaded advantage comes from the algorithmic difference: QEPG Monte Carlo costs O(k) per sample (k ~ Np expected faults), and QEPG Fixed-Weight costs O(w), while Stim simulates the full circuit at O(circuit_size).

Multi-threaded comparison (32 cores, 100K shots, p=0.001)

With OpenMP enabled, QEPG parallelizes across shots:

Decoding (pymatching) excluded to isolate sampling performance. Reproduce with python benchmark_sampling_speed.py. Note: Stim's detector sampler is single-threaded with SIMD; QEPG uses OpenMP for multi-threaded results.

How ScaLER works

Main method

ScaLERQEC estimates the LER by stratified fault-sampling and curve fitting:

We propose a novel method which tests the logical error rate by stratified sampling and curve fitting. See the tutorial for a detailed explanation. With a fixed QEC circuit and the noise model, we provide a simple interface for this method.

Roadmap

Completed

• Support installation via pip install
• Higher-level, easier interface to generate QEC programs
• Add cross-platform installation support (Windows, macOS, Linux)
• Python interface to construct QEC circuits from stabilizers
• Write full documentation (Sphinx)
• Support LDPC codes and LDPC code decoders (ScalerLDPC with BP+OSD)
• Remove Boost dependency -- use custom DynamicBitset with SIMD acceleration
• SIMD support (AVX2/SSE2/NEON) with cache-line aligned FlatBitTable
• Non-uniform noise model support (DEPOLARIZE1/2, X/Y/Z_ERROR, PAULI_CHANNEL_1/2)
• OpenMP parallel sampling across threads
• CI/CD pipeline with GitHub Actions (lint, build, test on Linux/macOS/Windows)
• Support toric codes, color codes, and BB LDPC codes
• Hotspot analysis -- C++-accelerated noise attribution with pluggable decoders
• Pluggable decoders -- all LER calculators accept user-provided decoders (pymatching, BPOSD, custom)

In Progress

• CUDA backend support

  • Port FlatBitTable::xor_row_into to CUDA kernel for GPU-accelerated GF(2) matmul
  • Implement GPU-side Poisson+CDF sparse sampler
  • Benchmark against Stim's SIMD sampler on large circuits (d >= 21)

• Enum-based gate dispatch for C++ backend

  • Replace Gate::name (std::string) with enum class GateKind in clifford.hpp
  • Convert string comparisons in backward_graph_construction() to switch statement
  • Expected 5-10x speedup on the backward traversal hot loop

• Refactor adaptive batching in monteLER.py

  • Extract the 5x copy-pasted adaptive loop into _adaptive_monte_carlo(sample_fn, decode_fn)
  • Reduces ~400 lines of duplication across calculate_LER_from_* methods

Planned

• Magic state distillation / cultivation

  • Implement 15-to-1 and 20-to-4 distillation protocols as Stim circuits
  • ScaLER estimation of magic state factory output error rate
  • Support post-selection in LER estimation

• Qiskit compatibility

  • Convert Qiskit QuantumCircuit to ScaLER StabCode or Stim circuit
  • Import noise models from Qiskit NoiseModel objects

• Advanced noise models

  • Decoherence noise (T1/T2 relaxation as Pauli channels)
  • Spatially correlated noise (crosstalk between neighboring qubits)
  • Leakage errors and leakage reduction units

• Lattice surgery and code switching

  • Support split/merge operations between surface code patches
  • LDPC code switching protocols
  • Estimate LER of multi-patch logical operations

• Visualization

  • Interactive QEPG graph visualization (networkx or D3.js)
  • Stim circuit diagram rendering

• Dynamic circuits

  • Support mid-circuit measurement and classical feedback
  • Compatible with IBM dynamic circuit model

• Static analysis pass

  • Detect symmetries in the circuit structure
  • Exploit symmetry to reduce sampling cost

• Pauli-measurement-based fault tolerance

  • Support circuits using Pauli measurements instead of CNOT+measure
  • Compile Pauli-based schemes to QEPG

Development Notes (for contributors)

Building from source

git clone https://github.com/yezhuoyang/ScaLERQEC.git
cd ScaLERQEC
pip install -e .

This compiles the C++ QEPG backend via pybind11. No external C++ libraries are needed -- the build system handles everything automatically.

Running tests

# Full test suite
pytest tests/ -v

# Quick smoke test
python -c "import scalerqec; import scalerqec.qepg; print('OK')"

Code formatting

We use ruff for Python formatting:

pip install ruff
ruff format src/scalerqec/
ruff check src/scalerqec/

Rebuilding the C++ backend

After modifying C++ source files under QEPG/src/:

pip install -e . --no-build-isolation