A surface code sweep

A distance-3 rotated surface code (17 qubits, 3 syndrome extraction rounds) built with Deltakit under SI1000Noise, benchmarked across three physical noise strengths with direct sampling, multi-worker sinter.collect, and coherent MPS simulation against its Pauli twirl.

At a glance

  • Guides: Sampling and Sinter · Coherent noise

  • Key APIs: MaestroBackend, to_stim(), compile_sampler(), as_sinter_task(), sinter_decoders(), twirled_circuit()

  • Simulators used: Maestro Stabilizer (Pauli frames) and Maestro MatrixProductState (χ=32)

  • Run locally: python examples/example_surface_code.py (~12 s)

Setup: Deltakit surface code with SI1000 noise

We construct a RotatedPlanarCode(3, 3) Z-basis memory circuit and attach superconducting SI1000Noise(p) using Deltakit’s QPU compiler. MaestroBackend accepts the resulting Deltakit Circuit directly (converting it internally via to_stim()).

def build_noisy_surface_code(p: float):
    """Build a Deltakit distance-d rotated surface code memory circuit with SI1000(p) noise.

    SI1000 models a superconducting architecture with 1-qubit depolarizing,
    2-qubit depolarizing, measurement bit-flip error, and data qubit idle noise.
    """
    code = RotatedPlanarCode(width=DISTANCE, height=DISTANCE)
    circuit = css_code_memory_circuit(code, num_rounds=NUM_ROUNDS, logical_basis=PauliBasis.Z)
    qpu = QPU(circuit.qubits, noise_model=SI1000Noise(p=p))
    return qpu.compile_and_add_noise_to_circuit(circuit)


Phase 1: Direct sampling checked against Stim

Because SI1000Noise contains only Pauli channels, MaestroBackend routes the run to Maestro’s Stabilizer simulator and samples it on Pauli frames. For each p in [0.002, 0.005, 0.01], we draw 20,000 shots with compile_sampler(), decode with PyMatching against sampler.task.detector_error_model, and compare against Stim’s own sampler on the exact same circuit.

for p, noisy_circuit in circuits.items():
    start_time = time.time()
    sampler = backend.compile_sampler(noisy_circuit, noise_type="pauli", seed=SEED)
    maestro_dets, maestro_obs = sampler.sample_detection_events(SHOTS)
    elapsed = time.time() - start_time

    dem = sampler.task.detector_error_model
    maestro_lep = decode_with_pymatching(dem, maestro_dets, maestro_obs)

    # Independent reference: Stim's sampler on the exact same circuit
    stim_circuit = to_stim(noisy_circuit)
    stim_dets, stim_obs = stim_circuit.compile_detector_sampler(seed=SEED).sample(
        SHOTS, separate_observables=True
    )
    stim_lep = decode_with_pymatching(dem, stim_dets, stim_obs)

    print(
        f"  {p:>6.3f} │ "
        f"{format_error_rate(maestro_lep):>15} │ "
        f"{format_error_rate(stim_lep):>17} │ "
        f"{elapsed:>11.1f}s"
    )


Output:

       p │         Maestro │ Stim SI1000 (ref) │ Maestro time
  ────────────────────────────────────────────────────────────
   0.002 │ 0.0144 ± 0.0008 │   0.0140 ± 0.0008 │         0.0s
   0.005 │ 0.0692 ± 0.0018 │   0.0708 ± 0.0018 │         0.0s
   0.010 │ 0.1952 ± 0.0028 │   0.1998 ± 0.0028 │         0.0s

Maestro and Stim agree within statistical error at every physical noise rate.

Phase 2: Parallel sweep with sinter.collect

To scale across CPU cores, as_sinter_task() packages each Deltakit circuit with its noise metadata and sinter_decoders(decoder="pymatching") dispatches sampling and decoding across 4 worker processes.

tasks = [
    backend.as_sinter_task(noisy_circuit, json_metadata={"p": p})
    for p, noisy_circuit in circuits.items()
]

start_time = time.time()
collected_stats = sinter.collect(
    num_workers=4,
    tasks=tasks,
    custom_decoders=backend.sinter_decoders(decoder="pymatching"),
    max_shots=SHOTS,
    max_batch_size=500,
)
total_elapsed = time.time() - start_time

Output:

       p │  shots │ errors │             LER
  ──────────────────────────────────────────
   0.002 │  20000 │    271 │ 0.0135 ± 0.0008
   0.005 │  20000 │   1419 │ 0.0709 ± 0.0018
   0.010 │  20000 │   3958 │ 0.1979 ± 0.0028
  (3.4s on 4 workers)

Phase 3: Coherent over-rotations (MPS) vs. their Pauli twirl

Passing noise_type="coherent" replaces gate depolarizing noise with unitary over-rotations while keeping SI1000’s idle, reset, and measurement noise intact. Maestro automatically switches to its Matrix Product State (χ=32) simulator. To isolate the effect of coherence, we compare against backend.twirled_circuit() on Maestro’s Stabilizer and decode both with the twirl’s DEM.

def phase_3_coherent_vs_twirl(backend: MaestroBackend, circuits: dict[float, Any]):
    """Simulate coherent unitary over-rotations on the surface code lattice using MPS."""
    print(f"\n{'─' * 70}")
    print("  3. Coherent noise (Maestro MPS) vs its Pauli twirl (Maestro stabilizer), same DEM")
    print(f"{'─' * 70}")
    print("  Coherent mode replaces gate depolarization with unitary over-rotations.")
    print(f"  Simulated via Matrix Product States (MPS) with bond dimension chi={MPS_BOND_DIM}.")
    print("  Both coherent runs and twirled runs are decoded using the twirl's DEM.\n")

    print(f"  {'p':>6} │ {'Maestro coherent':>16} │ {'Pauli twirl':>16} │ {'coherent time':>13}")
    print(f"  {'─' * 60}")

    for p, noisy_circuit in circuits.items():
        start_time = time.time()
        coherent_sampler = backend.compile_sampler(noisy_circuit, noise_type="coherent", seed=SEED)
        coherent_dets, coherent_obs = coherent_sampler.sample_detection_events(COHERENT_SHOTS)
        elapsed = time.time() - start_time

        dem = coherent_sampler.task.detector_error_model
        coherent_lep = decode_with_pymatching(dem, coherent_dets, coherent_obs)

        # Build the twirled circuit and sample it on Maestro's stabilizer
        twirled_circuit = backend.twirled_circuit(noisy_circuit, noise_type="coherent")
        twirl_sampler = backend.compile_sampler(twirled_circuit, seed=SEED)
        twirl_dets, twirl_obs = twirl_sampler.sample_detection_events(SHOTS)
        twirl_lep = decode_with_pymatching(dem, twirl_dets, twirl_obs)

        print(
            f"  {p:>6.3f} │ "
            f"{format_error_rate(coherent_lep):>16} │ "
            f"{format_error_rate(twirl_lep):>16} │ "
            f"{elapsed:>12.1f}s"
        )


Output:

       p │ Maestro coherent │      Pauli twirl │ coherent time
  ────────────────────────────────────────────────────────────
   0.002 │  0.0367 ± 0.0109 │  0.0181 ± 0.0009 │         12.6s
   0.005 │  0.1300 ± 0.0194 │  0.0872 ± 0.0020 │         12.8s
   0.010 │  0.2833 ± 0.0260 │  0.2323 ± 0.0030 │         12.9s

At every noise strength, coherent over-rotations yield a higher logical error probability than their matched Pauli twirl.

Key takeaways

  • Deltakit Circuit objects pass directly into compile_sampler() and as_sinter_task() with no manual translation step.

  • Always compare noise_type="coherent" against twirled_circuit(..., noise_type="coherent") rather than the raw SI1000 circuit, because the default coherent strength calibration differs from SI1000’s 1q/2q gate ratio.

  • Next: Every Deltakit code family runs the same pipeline across five code families, including a non-graphlike qLDPC code.


Full script and complete output

View full script (examples/example_surface_code.py)
"""
Surface Code QEC Benchmark — Riverlane Deltakit & Qoro MaestroBackend
=====================================================================

This benchmark combines Riverlane's Deltakit framework for designing topological
error correction codes with Qoro's MaestroBackend simulation engine.

We simulate a distance-3 rotated surface code memory experiment (17 physical qubits,
3 syndrome extraction rounds) under realistic superconducting QPU noise.

Demonstrates:
    1. Circuit Generation with SI1000 Noise:
       Using Deltakit's `RotatedPlanarCode` and `QPU` to build a surface code
       memory experiment compiled with superconducting hardware noise (SI1000).
    2. Direct Detection-Event Sampling & Decoding:
       Sampling detection events directly using `MaestroBackend.compile_sampler()`
       and decoding them with PyMatching (Minimum Weight Perfect Matching).
       SI1000 noise is all Pauli, so Maestro's stabilizer samples it on Pauli
       frames; results are checked against Stim's own sampler as a reference.
    3. Multi-Worker Sinter Sweeps:
       Distributing Monte Carlo sampling across multiple worker processes
       using `sinter.collect()` with backend tasks.
    4. Coherent Noise vs Pauli Twirl Comparison:
       Simulating genuine non-Clifford coherent over-rotations on the 17-qubit
       surface code lattice using Maestro's Matrix Product State (MPS) tensor
       network engine (bond dimension chi=32), compared to its Pauli twirl, which
       Maestro's stabilizer samples.

Usage:
    pip install -e ".[examples]"    # from project root (deltakit-maestro/)
    python examples/example_surface_code.py
"""

# Deltakit must be imported before pymatching: with pymatching 2.4 loaded
# first, Deltakit's circuit builders fail with "Caught an unknown exception!".
from deltakit.circuit.gates import PauliBasis
from deltakit.explorer.codes import RotatedPlanarCode, css_code_memory_circuit
from deltakit.explorer.qpu import QPU, SI1000Noise

# isort: split
import time
from typing import Any

import numpy as np
import pymatching
import sinter

from deltakit_maestro import MaestroBackend, to_stim

# =============================================================================
# EXPERIMENT CONFIGURATION
# =============================================================================

DISTANCE = 3  # Code distance (d=3 rotated surface code: 17 physical qubits)
NUM_ROUNDS = 3  # Number of syndrome extraction rounds
MPS_BOND_DIM = 32  # Maximum MPS bond dimension (chi) for coherent simulation
NOISE_STRENGTHS = [0.002, 0.005, 0.01]  # Physical noise parameter (p)
SHOTS = 20_000  # Shots per data point for Pauli runs (Pauli frames, fast)
COHERENT_SHOTS = 300  # Shots for MPS simulation (~60 ms/shot per core on 17 qubits)
SEED = 1234


# =============================================================================
# CIRCUIT BUILDER & DECODER UTILITIES
# =============================================================================


def build_noisy_surface_code(p: float):
    """Build a Deltakit distance-d rotated surface code memory circuit with SI1000(p) noise.

    SI1000 models a superconducting architecture with 1-qubit depolarizing,
    2-qubit depolarizing, measurement bit-flip error, and data qubit idle noise.
    """
    code = RotatedPlanarCode(width=DISTANCE, height=DISTANCE)
    circuit = css_code_memory_circuit(code, num_rounds=NUM_ROUNDS, logical_basis=PauliBasis.Z)
    qpu = QPU(circuit.qubits, noise_model=SI1000Noise(p=p))
    return qpu.compile_and_add_noise_to_circuit(circuit)


def decode_with_pymatching(
    dem: Any,
    detection_events: np.ndarray,
    observable_flips: np.ndarray,
) -> tuple[float, float]:
    """Decode detection events using Minimum Weight Perfect Matching (MWPM).

    Returns (logical_error_rate, standard_error).
    """
    matcher = pymatching.Matching.from_detector_error_model(dem)
    predicted_flips = matcher.decode_batch(detection_events)
    failures = np.any(predicted_flips != observable_flips, axis=1)
    logical_error_rate = float(failures.mean())
    std_error = float(np.sqrt(logical_error_rate * (1 - logical_error_rate) / len(failures)))
    return logical_error_rate, std_error


def format_error_rate(rate_err: tuple[float, float]) -> str:
    """Format an error rate and its standard error as '0.XXXX ± 0.XXXX'."""
    return f"{rate_err[0]:.4f} ± {rate_err[1]:.4f}"


# =============================================================================
# BENCHMARK PHASES
# =============================================================================


def phase_1_direct_sampling(backend: MaestroBackend, circuits: dict[float, Any]):
    """Demonstrate direct detection-event sampling and PyMatching decoding, checked against Stim."""
    print(f"\n{'─' * 70}")
    print("  1. Direct sampling + PyMatching: Logical Error Probability (LEP)")
    print(f"{'─' * 70}")
    print("  SI1000 noise is all Pauli, so Maestro's stabilizer samples it on Pauli frames.")
    print("  It is checked against Stim's own sampler on the same circuit, decoded with the")
    print("  same DEM; the two agree to within sampling error.\n")

    print(f"  {'p':>6} │ {'Maestro':>15} │ {'Stim SI1000 (ref)':>17} │ {'Maestro time':>12}")
    print(f"  {'─' * 60}")

    for p, noisy_circuit in circuits.items():
        start_time = time.time()
        sampler = backend.compile_sampler(noisy_circuit, noise_type="pauli", seed=SEED)
        maestro_dets, maestro_obs = sampler.sample_detection_events(SHOTS)
        elapsed = time.time() - start_time

        dem = sampler.task.detector_error_model
        maestro_lep = decode_with_pymatching(dem, maestro_dets, maestro_obs)

        # Independent reference: Stim's sampler on the exact same circuit
        stim_circuit = to_stim(noisy_circuit)
        stim_dets, stim_obs = stim_circuit.compile_detector_sampler(seed=SEED).sample(
            SHOTS, separate_observables=True
        )
        stim_lep = decode_with_pymatching(dem, stim_dets, stim_obs)

        print(
            f"  {p:>6.3f} │ "
            f"{format_error_rate(maestro_lep):>15} │ "
            f"{format_error_rate(stim_lep):>17} │ "
            f"{elapsed:>11.1f}s"
        )


def phase_2_sinter_sweep(backend: MaestroBackend, circuits: dict[float, Any]):
    """Demonstrate parallel Monte Carlo sampling over a noise sweep using Sinter."""
    print(f"\n{'─' * 70}")
    print("  2. sinter.collect (Pauli noise on Maestro's Pauli frames, PyMatching)")
    print(f"{'─' * 70}")
    print("  Distributing tasks across 4 worker processes via Sinter.")
    print("  Each worker compiles and samples tasks independently.\n")

    tasks = [
        backend.as_sinter_task(noisy_circuit, json_metadata={"p": p})
        for p, noisy_circuit in circuits.items()
    ]

    start_time = time.time()
    collected_stats = sinter.collect(
        num_workers=4,
        tasks=tasks,
        custom_decoders=backend.sinter_decoders(decoder="pymatching"),
        max_shots=SHOTS,
        max_batch_size=500,
    )
    total_elapsed = time.time() - start_time

    print(f"  {'p':>6} │ {'shots':>6} │ {'errors':>6} │ {'LER':>15}")
    print(f"  {'─' * 42}")
    for stat in sorted(collected_stats, key=lambda s: s.json_metadata["p"]):
        ler = stat.errors / stat.shots
        err = np.sqrt(ler * (1 - ler) / stat.shots)
        print(
            f"  {stat.json_metadata['p']:>6.3f} │ {stat.shots:>6} │ {stat.errors:>6} │ {ler:.4f} ± {err:.4f}"
        )
    print(f"  ({total_elapsed:.1f}s on 4 workers)")


def phase_3_coherent_vs_twirl(backend: MaestroBackend, circuits: dict[float, Any]):
    """Simulate coherent unitary over-rotations on the surface code lattice using MPS."""
    print(f"\n{'─' * 70}")
    print("  3. Coherent noise (Maestro MPS) vs its Pauli twirl (Maestro stabilizer), same DEM")
    print(f"{'─' * 70}")
    print("  Coherent mode replaces gate depolarization with unitary over-rotations.")
    print(f"  Simulated via Matrix Product States (MPS) with bond dimension chi={MPS_BOND_DIM}.")
    print("  Both coherent runs and twirled runs are decoded using the twirl's DEM.\n")

    print(f"  {'p':>6} │ {'Maestro coherent':>16} │ {'Pauli twirl':>16} │ {'coherent time':>13}")
    print(f"  {'─' * 60}")

    for p, noisy_circuit in circuits.items():
        start_time = time.time()
        coherent_sampler = backend.compile_sampler(noisy_circuit, noise_type="coherent", seed=SEED)
        coherent_dets, coherent_obs = coherent_sampler.sample_detection_events(COHERENT_SHOTS)
        elapsed = time.time() - start_time

        dem = coherent_sampler.task.detector_error_model
        coherent_lep = decode_with_pymatching(dem, coherent_dets, coherent_obs)

        # Build the twirled circuit and sample it on Maestro's stabilizer
        twirled_circuit = backend.twirled_circuit(noisy_circuit, noise_type="coherent")
        twirl_sampler = backend.compile_sampler(twirled_circuit, seed=SEED)
        twirl_dets, twirl_obs = twirl_sampler.sample_detection_events(SHOTS)
        twirl_lep = decode_with_pymatching(dem, twirl_dets, twirl_obs)

        print(
            f"  {p:>6.3f} │ "
            f"{format_error_rate(coherent_lep):>16} │ "
            f"{format_error_rate(twirl_lep):>16} │ "
            f"{elapsed:>12.1f}s"
        )


# =============================================================================
# MAIN ENTRYPOINT
# =============================================================================


def main():
    n_qubits = 2 * DISTANCE**2 - 1

    print("=" * 70)
    print("  Surface Code QEC Benchmark — MaestroBackend & Deltakit")
    print("=" * 70)
    print(f"  Code:          d={DISTANCE} rotated surface code ({n_qubits} physical qubits)")
    print(f"  Rounds:        {NUM_ROUNDS} syndrome extraction rounds")
    print(f"  Simulator:     Maestro Stabilizer (Pauli), Maestro MPS chi={MPS_BOND_DIM} (coherent)")
    print(f"  Shots:         {SHOTS} per Pauli point, {COHERENT_SHOTS} per coherent point")

    backend = MaestroBackend(max_bond_dimension=MPS_BOND_DIM)
    circuits = {p: build_noisy_surface_code(p) for p in NOISE_STRENGTHS}

    # 1. Direct sampling & decoding compared to Stim
    phase_1_direct_sampling(backend, circuits)

    # 2. Parallel Sinter sweep
    phase_2_sinter_sweep(backend, circuits)

    # 3. Coherent over-rotations vs Pauli twirling
    phase_3_coherent_vs_twirl(backend, circuits)

    print(f"\n{'=' * 70}")
    print("  Notes:")
    print("    - Pauli runs: SI1000's noise is all Pauli, so Maestro's stabilizer samples")
    print("      it exactly on Pauli frames. Coherent runs need its MPS simulator, since")
    print("      no Clifford simulator can run a coherent over-rotation.")
    print("    - Coherence Gap: Coherent mode replaces gate noise with unitary")
    print("      over-rotations while retaining idle, reset, and readout noise.")
    print("      Any divergence from the twirl column reveals coherent interference.")
    print("=" * 70)


if __name__ == "__main__":
    main()
View complete output
======================================================================
  Surface Code QEC Benchmark — MaestroBackend & Deltakit
======================================================================
  Code:          d=3 rotated surface code (17 physical qubits)
  Rounds:        3 syndrome extraction rounds
  Simulator:     Maestro Stabilizer (Pauli), Maestro MPS chi=32 (coherent)
  Shots:         20000 per Pauli point, 300 per coherent point

──────────────────────────────────────────────────────────────────────
  1. Direct sampling + PyMatching: Logical Error Probability (LEP)
──────────────────────────────────────────────────────────────────────
  SI1000 noise is all Pauli, so Maestro's stabilizer samples it on Pauli frames.
  It is checked against Stim's own sampler on the same circuit, decoded with the
  same DEM; the two agree to within sampling error.

       p │         Maestro │ Stim SI1000 (ref) │ Maestro time
  ────────────────────────────────────────────────────────────
   0.002 │ 0.0144 ± 0.0008 │   0.0140 ± 0.0008 │         0.0s
   0.005 │ 0.0692 ± 0.0018 │   0.0708 ± 0.0018 │         0.0s
   0.010 │ 0.1952 ± 0.0028 │   0.1998 ± 0.0028 │         0.0s

──────────────────────────────────────────────────────────────────────
  2. sinter.collect (Pauli noise on Maestro's Pauli frames, PyMatching)
──────────────────────────────────────────────────────────────────────
  Distributing tasks across 4 worker processes via Sinter.
  Each worker compiles and samples tasks independently.

       p │  shots │ errors │             LER
  ──────────────────────────────────────────
   0.002 │  20000 │    271 │ 0.0135 ± 0.0008
   0.005 │  20000 │   1419 │ 0.0709 ± 0.0018
   0.010 │  20000 │   3958 │ 0.1979 ± 0.0028
  (3.4s on 4 workers)

──────────────────────────────────────────────────────────────────────
  3. Coherent noise (Maestro MPS) vs its Pauli twirl (Maestro stabilizer), same DEM
──────────────────────────────────────────────────────────────────────
  Coherent mode replaces gate depolarization with unitary over-rotations.
  Simulated via Matrix Product States (MPS) with bond dimension chi=32.
  Both coherent runs and twirled runs are decoded using the twirl's DEM.

       p │ Maestro coherent │      Pauli twirl │ coherent time
  ────────────────────────────────────────────────────────────
   0.002 │  0.0367 ± 0.0109 │  0.0181 ± 0.0009 │         12.6s
   0.005 │  0.1300 ± 0.0194 │  0.0872 ± 0.0020 │         12.8s
   0.010 │  0.2833 ± 0.0260 │  0.2323 ± 0.0030 │         12.9s

======================================================================
  Notes:
    - Pauli runs: SI1000's noise is all Pauli, so Maestro's stabilizer samples
      it exactly on Pauli frames. Coherent runs need its MPS simulator, since
      no Clifford simulator can run a coherent over-rotation.
    - Coherence Gap: Coherent mode replaces gate noise with unitary
      over-rotations while retaining idle, reset, and readout noise.
      Any divergence from the twirl column reveals coherent interference.
======================================================================