{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "frontmatter",
      "metadata": {},
      "source": [
        "---\n",
        "title: Low-overhead error detection with spacetime codes\n",
        "description: Automatically insert spacetime Pauli checks with qiskit-paulice and boost sampled Clifford-circuit fidelity by postselecting on check syndromes.\n",
        "---\n",
        "\n",
        "{/* cspell:ignore postsel, frameon, virt */}\n",
        "\n",
        "# Low-overhead error detection with spacetime codes\n",
        "\n",
        "*Usage estimate: 4 minutes on a Heron processor (`ibm_kingston` or equivalent) (NOTE: This is an estimate only. Your runtime may vary.)*\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "learning-outcomes",
      "metadata": {},
      "source": [
        "## Learning outcomes\n",
        "\n",
        "After completing this tutorial, you can expect to understand the following information:\n",
        "\n",
        "* How spacetime Pauli checks detect logical errors in Clifford circuits, and how postselecting on their syndromes boosts the fidelity of a sampled distribution.\n",
        "* How to use the `qiskit-paulice` package to find and insert hardware-efficient checks automatically with `get_check_qubits`, `NoiseModel`, and `add_pauli_checks`.\n",
        "* How to estimate the fidelity of a stabilizer state by sampling its stabilizers and postselecting on check syndromes.\n",
        "* How to run the full error-detection workflow on IBM Quantum® hardware and compare noisy and postselected fidelities.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "prerequisites",
      "metadata": {},
      "source": [
        "## Prerequisites\n",
        "\n",
        "It is recommended that you familiarize yourself with these topics:\n",
        "\n",
        "* [Hardware fundamentals](/learning/courses/utility-scale-quantum-computing/hardware) for utility-scale quantum computing.\n",
        "* The Clifford and stabilizer formalism, including how a stabilizer group describes a pure stabilizer state.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "background",
      "metadata": {},
      "source": [
        "## Background\n",
        "\n",
        "[Low-overhead error detection with spacetime codes](https://arxiv.org/abs/2504.15725) [\\[1\\]](#references) by Simon Martiel and Ali Javadi-Abhari introduces a method for detecting logical errors in Clifford-dominated circuits that sits between full error correction and lighter-weight error mitigation. The idea builds on coherent Pauli checks (CPC) from [Single-shot error mitigation by coherent Pauli checks](https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.5.033193) [\\[2\\]](#references) by van den Berg and others. In both approaches, a Clifford \"payload\" circuit is entangled with ancilla qubits to check certain invariants. Measuring the ancillas produces a syndrome that reports whether an error was detected during execution. Keeping only the samples with no detected error improves the fidelity of the sampled distribution, at the cost of a reduced postselection rate.\n",
        "\n",
        "The key difference between coherent Pauli checks and spacetime checks is the operators they measure. Coherent Pauli checks measure time-localized, high-weight operators. On qubit topologies with limited connectivity, such as heavy hex, those checks need many SWAP gates and often make the circuit too deep to run in practice. Implementing the checks as spacetime codes instead distributes each check across the payload circuit in space and time. This yields a hardware-efficient encoding that stays effective at detecting logical errors while keeping the qubit and depth overhead low.\n",
        "\n",
        "### What the `qiskit-paulice` package does\n",
        "\n",
        "The `qiskit-paulice` package automates the construction of these checks so that you do not have to build them by hand. Its primary role is to find and insert valid spacetime Pauli checks at the locations in a circuit that maximize error detection while minimizing qubit overhead. A check is *valid* when its operators leave the logical action of the payload circuit unchanged, *low-weight* when it uses few entangling gates, and *effective* when it detects a large portion of the errors, relative to the noise the check itself introduces. The package scores candidate checks against a noise model and commits the best ones to the circuit. This tutorial uses three API methods:\n",
        "\n",
        "* `get_check_qubits` inspects a backend coupling map and returns target and ancilla qubit pairs. A check on `target_qubits[i]` uses `ancilla_qubits[i]`.\n",
        "* `NoiseModel.from_backend` builds a rough noise model from backend benchmark data. The model scores candidate checks, so an exact, learned noise model is not required. For a learned Pauli-Lindblad model, see [`NoiseModel.from_pauli_lindblad_maps`](/docs/api/qiskit-addon-paulice/qiskit-paulice-noise-models#from_pauli_lindblad_maps).\n",
        "* `add_pauli_checks` finds and inserts checks into a circuit. It returns a sequence of `CheckedCircuit` objects with an increasing number of checks, and each object provides a `get_postselection_method` that maps a measured bitstring to a syndrome vector. The `cost` argument selects the function that scores a check (`gamma`, the sampling overhead of the postselected inverse noise channel, or `LER`, the logical error rate). The `method` argument selects the search strategy (`windowed`, `genetic`, or `windowed_genetic`). This tutorial uses `cost=\"gamma\"` and `method=\"windowed\"`, which together give deterministic, reproducible check selection.\n",
        "\n",
        "### Estimate fidelity from stabilizer sampling\n",
        "\n",
        "To measure how well the error detection works, you can estimate the fidelity of the stabilizer state $|\\psi\\rangle = U|0\\rangle^{\\otimes n}$ that the circuit ideally prepares against the noisy state $\\rho$ that the hardware actually outputs. The projector onto a pure stabilizer state $|\\psi\\rangle$ equals the uniform average over the $2^n$ elements of its stabilizer group $\\mathcal{S}$:\n",
        "\n",
        "$|\\psi\\rangle\\langle\\psi| = \\frac{1}{2^n}\\sum_{G \\in \\mathcal{S}} G.$\n",
        "\n",
        "Substituting this into the fidelity gives the fidelity of $\\rho$ as the average expectation value of every stabilizer $G \\in \\mathcal{S}$ with respect to $\\rho$:\n",
        "\n",
        "$F = \\mathrm{Tr}(\\rho|\\psi\\rangle\\langle\\psi|) = \\frac{1}{2^n} \\sum_{G \\in \\mathcal{S}} \\mathrm{Tr}(\\rho G) = \\frac{1}{2^n} \\sum_{G \\in \\mathcal{S}} \\langle G \\rangle_\\rho.$\n",
        "\n",
        "For larger problems, enumerating all $2^n$ stabilizers is infeasible, so you can estimate the fidelity from a random sample. Drawing $M$ stabilizers $G_1, \\ldots, G_M$ uniformly at random from $\\mathcal{S}$ gives an unbiased estimation:\n",
        "\n",
        "$\\hat F_M = \\frac{1}{M} \\sum_{i=1}^M \\langle G_i \\rangle_\\rho.$\n",
        "\n",
        "Because a Clifford circuit prepares a stabilizer state, you can estimate its fidelity directly from sampled expectation values of its stabilizers. This tutorial first walks through the workflow on a simulator with a small circuit, then runs the same workflow on hardware with a larger, deeper circuit. As circuits include more non-Clifford operations, the number of valid checks shrinks quickly, so the method works best for Clifford-dominated circuits.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "requirements",
      "metadata": {},
      "source": [
        "## Requirements\n",
        "\n",
        "Before starting this tutorial, be sure you have the following installed:\n",
        "\n",
        "* Qiskit SDK v2.0 or later, with [visualization](/docs/api/qiskit/visualization) support\n",
        "* Qiskit Runtime v0.40 or later (`pip install qiskit-ibm-runtime`)\n",
        "* Qiskit Aer v0.17 or later (`pip install qiskit-aer`)\n",
        "* Qiskit Paulice (`pip install qiskit-paulice`)\n",
        "* tqdm (`pip install tqdm`)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "setup",
      "metadata": {},
      "source": [
        "## Setup\n",
        "\n",
        "Import the required libraries and define the helper functions that are not available as imports. The `random_clifford_circuit` function builds a brickwork random Clifford payload, `find_check_layout` searches a backend coupling map for a low-error qubit path with many available ancillas, `learned_noise_model` turns `NoiseLearner` output into a `qiskit-paulice` noise model, `append_basis_rotation` rotates a circuit so that a stabilizer is measured in the computational basis, `expectation` computes a stabilizer expectation value from sampled counts, and `cum_mean_sem` tracks the running fidelity estimate.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "setup-imports",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Standard library imports\n",
        "import random\n",
        "import time\n",
        "\n",
        "# External libraries\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "from tqdm import tqdm\n",
        "\n",
        "# Qiskit\n",
        "from qiskit import QuantumCircuit\n",
        "from qiskit.quantum_info import Clifford, Pauli, PauliLindbladMap, PauliList\n",
        "from qiskit.result import sampled_expectation_value\n",
        "from qiskit.transpiler import generate_preset_pass_manager\n",
        "from qiskit.visualization import plot_coupling_map\n",
        "\n",
        "# Qiskit Aer\n",
        "from qiskit_aer import AerSimulator\n",
        "from qiskit_aer.noise import NoiseModel as AerNoiseModel\n",
        "from qiskit_aer.noise import ReadoutError, depolarizing_error\n",
        "\n",
        "# Qiskit IBM Runtime\n",
        "from qiskit_ibm_runtime import NoiseLearner, QiskitRuntimeService\n",
        "from qiskit_ibm_runtime import SamplerV2 as Sampler\n",
        "\n",
        "# Qiskit Paulice\n",
        "from qiskit_paulice import add_pauli_checks\n",
        "from qiskit_paulice.layout import get_check_qubits\n",
        "from qiskit_paulice.noise_models import NoiseModel"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "setup-helpers",
      "metadata": {},
      "outputs": [],
      "source": [
        "def random_clifford_circuit(\n",
        "    num_qubits: int, depth: int, rng: np.random.Generator\n",
        ") -> QuantumCircuit:\n",
        "    \"\"\"Brickwork random Clifford on `num_qubits`, with `depth` CZ layers.\"\"\"\n",
        "    qc = QuantumCircuit(num_qubits)\n",
        "    qc.h(range(num_qubits))\n",
        "    for d in range(depth):\n",
        "        for i in range(d % 2, num_qubits - 1, 2):\n",
        "            qc.cz(i, i + 1)\n",
        "        for q in range(num_qubits):\n",
        "            if rng.integers(0, 2):\n",
        "                qc.sx(q)\n",
        "            if rng.integers(0, 2):\n",
        "                qc.s(q)\n",
        "            if rng.integers(0, 2):\n",
        "                qc.sx(q)\n",
        "    return qc\n",
        "\n",
        "\n",
        "def find_check_layout(\n",
        "    backend,\n",
        "    num_qubits: int,\n",
        "    rng: np.random.Generator,\n",
        "    num_trials: int = 200,\n",
        "    max_gate_error: float = 0.03,\n",
        "    max_readout_error: float = 0.2,\n",
        ") -> list[int]:\n",
        "    \"\"\"Find a low-error path of `num_qubits` qubits with many available ancillas.\n",
        "\n",
        "    Builds random self-avoiding walks on the coupling map, excluding the qubits\n",
        "    and two-qubit gates whose reported errors exceed the thresholds, and keeps\n",
        "    the path that offers the most target and ancilla pairs. Ties are broken by\n",
        "    the lower average two-qubit gate error along the path.\n",
        "    \"\"\"\n",
        "    target = backend.target\n",
        "    gate_2q = next(\n",
        "        name for name in (\"cz\", \"ecr\", \"cx\") if name in target.operation_names\n",
        "    )\n",
        "\n",
        "    # Collect per-edge gate errors and per-qubit readout errors\n",
        "    edge_error = {}\n",
        "    for qubits, props in target[gate_2q].items():\n",
        "        edge = tuple(sorted(qubits))\n",
        "        if props is not None and props.error is not None:\n",
        "            edge_error[edge] = min(edge_error.get(edge, 1.0), props.error)\n",
        "    readout_error = {\n",
        "        qubit: target[\"measure\"][(qubit,)].error\n",
        "        for (qubit,) in target[\"measure\"]\n",
        "    }\n",
        "\n",
        "    # Keep only the edges whose gate and readout errors are acceptable\n",
        "    adjacency = {}\n",
        "    for (q1, q2), error in edge_error.items():\n",
        "        if (\n",
        "            error <= max_gate_error\n",
        "            and readout_error.get(q1, 1.0) <= max_readout_error\n",
        "            and readout_error.get(q2, 1.0) <= max_readout_error\n",
        "        ):\n",
        "            adjacency.setdefault(q1, set()).add(q2)\n",
        "            adjacency.setdefault(q2, set()).add(q1)\n",
        "\n",
        "    # Random self-avoiding walks; keep the path with the most check pairs\n",
        "    starts = sorted(adjacency)\n",
        "    best_path = None\n",
        "    best_score = (-1, float(\"inf\"))\n",
        "    for _ in range(num_trials):\n",
        "        path = [starts[rng.integers(len(starts))]]\n",
        "        while len(path) < num_qubits:\n",
        "            options = sorted(adjacency[path[-1]] - set(path))\n",
        "            if not options:\n",
        "                break\n",
        "            path.append(options[rng.integers(len(options))])\n",
        "        if len(path) < num_qubits:\n",
        "            continue\n",
        "        num_pairs = len(get_check_qubits(backend.coupling_map, path)[0])\n",
        "        mean_error = float(\n",
        "            np.mean(\n",
        "                [edge_error[tuple(sorted(e))] for e in zip(path, path[1:])]\n",
        "            )\n",
        "        )\n",
        "        if num_pairs > best_score[0] or (\n",
        "            num_pairs == best_score[0] and mean_error < best_score[1]\n",
        "        ):\n",
        "            best_path, best_score = path, (num_pairs, mean_error)\n",
        "\n",
        "    if best_path is None:\n",
        "        raise RuntimeError(\n",
        "            \"No connected low-error path found. Relax the error thresholds.\"\n",
        "        )\n",
        "    return best_path\n",
        "\n",
        "\n",
        "def learned_noise_model(layer_errors, layout: list[int]) -> NoiseModel:\n",
        "    \"\"\"Build a `NoiseModel` from `NoiseLearner` results.\n",
        "\n",
        "    `NoiseLearner` reports one `PauliLindbladError` per entangling layer, whose\n",
        "    generators are indexed against that layer's own physical qubits, while\n",
        "    `NoiseModel.from_pauli_lindblad_maps` expects `PauliLindbladMap`s indexed the\n",
        "    way `NoiseModel.from_backend` indexes them: by position in `layout`. This\n",
        "    translates between the two and drops generators that fall outside `layout`.\n",
        "    \"\"\"\n",
        "    phys_to_virt = {phys: virt for virt, phys in enumerate(layout)}\n",
        "    maps = []\n",
        "    for layer in layer_errors:\n",
        "        if layer.error is None:\n",
        "            continue\n",
        "        terms = []\n",
        "        for pauli, rate in zip(\n",
        "            layer.error.generators, layer.error.rates, strict=True\n",
        "        ):\n",
        "            label, indices = [], []\n",
        "            for local, phys in enumerate(layer.qubits):\n",
        "                x, z = bool(pauli.x[local]), bool(pauli.z[local])\n",
        "                if not (x or z):\n",
        "                    continue\n",
        "                if phys not in phys_to_virt:\n",
        "                    break  # generator reaches outside the layout, so skip it\n",
        "                label.append(\"Y\" if x and z else \"X\" if x else \"Z\")\n",
        "                indices.append(phys_to_virt[phys])\n",
        "            else:\n",
        "                if label:\n",
        "                    terms.append(\n",
        "                        (\"\".join(label), tuple(indices), float(rate))\n",
        "                    )\n",
        "        # Each map needs a 2-qubit generator to define an entangling layer\n",
        "        if any(len(t[1]) == 2 for t in terms):\n",
        "            maps.append(\n",
        "                PauliLindbladMap.from_sparse_list(\n",
        "                    terms, num_qubits=len(layout)\n",
        "                )\n",
        "            )\n",
        "    if not maps:\n",
        "        raise RuntimeError(\n",
        "            \"No usable layer errors. Check that the learner ran on this layout.\"\n",
        "        )\n",
        "    return NoiseModel.from_pauli_lindblad_maps(maps)\n",
        "\n",
        "\n",
        "def append_basis_rotation(\n",
        "    circuit: QuantumCircuit, pauli: Pauli\n",
        ") -> QuantumCircuit:\n",
        "    \"\"\"Strip measurements, append basis rotations for `pauli`, and re-measure.\"\"\"\n",
        "    out = circuit.remove_final_measurements(inplace=False)\n",
        "    for q in range(pauli.num_qubits):\n",
        "        if pauli.x[q]:\n",
        "            if pauli.z[q]:\n",
        "                out.sdg(q)\n",
        "            out.h(q)\n",
        "    out.measure_all()\n",
        "    return out\n",
        "\n",
        "\n",
        "def expectation(counts: dict, pauli: Pauli) -> float:\n",
        "    \"\"\"Expectation value of `pauli` from counts measured in the Z basis.\n",
        "\n",
        "    Pads with identity on any qubits beyond the support of `pauli`, such as the\n",
        "    check ancillas that appear in the postselected counts.\n",
        "    \"\"\"\n",
        "    if not counts:\n",
        "        return float(\"nan\")\n",
        "    n = pauli.num_qubits\n",
        "    sign = -1 if int(pauli.phase) % 4 == 2 else 1\n",
        "    total = len(next(iter(counts)))\n",
        "    label = \"\".join(\n",
        "        \"Z\" if q < n and (pauli.x[q] or pauli.z[q]) else \"I\"\n",
        "        for q in range(total - 1, -1, -1)\n",
        "    )\n",
        "    return sign * sampled_expectation_value(counts, label)\n",
        "\n",
        "\n",
        "def cum_mean_sem(values: np.ndarray):\n",
        "    \"\"\"Cumulative mean and standard error of the mean, ignoring NaNs.\"\"\"\n",
        "    valid = ~np.isnan(values)\n",
        "    total = np.cumsum(np.where(valid, values, 0.0))\n",
        "    total_sq = np.cumsum(np.where(valid, values**2, 0.0))\n",
        "    count = np.maximum(np.cumsum(valid).astype(float), 1)\n",
        "    mean = total / count\n",
        "    sem = np.sqrt(np.maximum(total_sq / count - mean**2, 0) / count)\n",
        "    return np.where(np.cumsum(valid) > 0, mean, np.nan), sem"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "small-scale-header",
      "metadata": {},
      "source": [
        "## Small-scale simulator example\n",
        "\n",
        "This section walks through the full workflow on a noisy simulator. It uses backend benchmark data to choose a qubit layout and a noise model, finds checks automatically, and uses postselection on the sampled distribution to show the fidelity improvement.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "step1-header",
      "metadata": {},
      "source": [
        "### Step 1: Map classical inputs to a quantum problem\n",
        "\n",
        "The payload circuit is a shallow one-dimensional brickwork random Clifford circuit. Because the circuit is Clifford, it prepares a stabilizer state whose fidelity you can estimate directly from sampled stabilizer expectation values. Start with a shallow circuit so that the checks are easy to visualize in the next step.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "step1-code",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/spacetime-codes/extracted-outputs/step1-code-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 3,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "num_qubits = 12\n",
        "depth = 4\n",
        "seed = 1764\n",
        "rng = np.random.default_rng(seed)\n",
        "np.random.seed(seed)\n",
        "\n",
        "circuit = random_clifford_circuit(num_qubits, depth, rng)\n",
        "circuit.measure_all()\n",
        "circuit.draw(\"mpl\", fold=-1, scale=0.6)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "step2-header",
      "metadata": {},
      "source": [
        "### Step 2: Optimize for quantum hardware execution\n",
        "\n",
        "Mapping the circuit to hardware sets the physical qubit layout, the noise model that scores candidate checks, and the checks themselves.\n",
        "\n",
        "First, select a backend and search its coupling map for a one-dimensional qubit layout with the `find_check_layout` helper defined in the Setup section. The helper builds random self-avoiding walks that avoid the highest-error gates and readouts, and it keeps the path that offers the most target and ancilla pairs. Because the search reads the connectivity and error data from the backend itself, the same code runs on any IBM Quantum QPU. The `get_check_qubits` function then returns the target and ancilla pairs, where a check on `target_qubits[i]` uses `ancilla_qubits[i]`.\n",
        "\n",
        "In the coupling graph that follows, the green qubits are payload qubits and the orange qubits are the ancillas that implement the checks. Qubits with an adjacent ancilla are used as target qubits for checks.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "step2-layout",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Backend: ibm_boston\n",
            "Target qubits:  [105, 107, 108, 123, 125, 141, 143]\n",
            "Ancilla qubits: [104, 97, 109, 122, 126, 140, 144]\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/spacetime-codes/extracted-outputs/step2-layout-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 4,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "service = QiskitRuntimeService()\n",
        "backend = service.least_busy(operational=True, simulator=False)\n",
        "\n",
        "print(f\"Backend: {backend.name}\")\n",
        "\n",
        "# Search for a low-error path, then pair each target qubit with a neighboring ancilla\n",
        "layout = find_check_layout(backend, num_qubits, rng)\n",
        "target_qubits, ancilla_qubits = get_check_qubits(backend, layout)\n",
        "num_checks = len(target_qubits)\n",
        "\n",
        "print(f\"Target qubits:  {target_qubits}\")\n",
        "print(f\"Ancilla qubits: {ancilla_qubits}\")\n",
        "plot_coupling_map(\n",
        "    num_qubits=backend.num_qubits,\n",
        "    qubit_coordinates=getattr(\n",
        "        backend.configuration(), \"qubit_coordinates\", None\n",
        "    ),\n",
        "    coupling_map=backend.configuration().coupling_map,\n",
        "    figsize=(12, 12),\n",
        "    qubit_color=[\n",
        "        \"#4CAF50\"\n",
        "        if i in set(layout)\n",
        "        else \"#FF9800\"\n",
        "        if i in set(ancilla_qubits)\n",
        "        else \"#DDDDDD\"\n",
        "        for i in backend.coupling_map.graph.node_indices()\n",
        "    ],\n",
        "    qubit_size=220,\n",
        "    line_width=2,\n",
        "    font_size=90,\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "step2-transpile-md",
      "metadata": {},
      "source": [
        "With the backend and layout chosen, transpile the payload into an instruction set architecture (ISA) circuit. It is only necessary to set the layout and translate the gates into the native gate set of the backend.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "step2-transpile",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/spacetime-codes/extracted-outputs/step2-transpile-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 5,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "pm = generate_preset_pass_manager(\n",
        "    optimization_level=0, backend=backend, initial_layout=layout\n",
        ")\n",
        "circuit_isa = pm.run(circuit)\n",
        "circuit_isa.draw(\"mpl\", fold=-1, scale=0.6)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "step2-noise-md",
      "metadata": {},
      "source": [
        "Next, model how the gate and readout noise on the backend affects execution. The noise model determines where in the circuit a check captures the most error. A more accurate model improves the detection, but it is usually not necessary to learn one by sampling the QPU. The model that follows infers a uniform depolarizing channel for gate and readout noise from `qiskit-ibm-runtime` benchmark data.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "step2-noise",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "NoiseModel(gate_noise=0.001079865281450939, readout_noise=0.006001790364583333, idling_noise=None)\n"
          ]
        }
      ],
      "source": [
        "noise_model = NoiseModel.from_backend(\n",
        "    backend, layout, uniform_gate_noise=True\n",
        ")\n",
        "print(noise_model)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "step2-checks-md",
      "metadata": {},
      "source": [
        "Now add checks to the circuit. The `add_pauli_checks` function takes the Clifford payload, the list of target qubits, and the noise model. The `ancilla_qubits` argument tells the function which physical ancilla to pair with each target. Checks are added in the order the target qubits appear, so the final layout of the checked circuit is `layout + ancilla_qubits`. To run an output circuit with fewer (`i`) checks, the final layout is `layout + ancilla_qubits[:i]`.\n",
        "\n",
        "The output of `add_pauli_checks` is a sequence of circuits with an increasing number of checks, from no checks up to one check on every target qubit. The visualization confirms that the checks use the specified target and ancilla pairs. For details on finding good checks, see Sections II to IV of the supplementary information in reference [\\[1\\]](#references).\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "step2-checks",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Physical layout of payload and ancillas: [108, 107, 106, 105, 117, 125, 124, 123, 136, 143, 142, 141, 104, 97, 109, 122, 126, 140, 144]\n",
            "Checked circuit:\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/spacetime-codes/extracted-outputs/step2-checks-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 7,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "checked = add_pauli_checks(\n",
        "    circuit_isa,\n",
        "    target_qubits,\n",
        "    noise_model,\n",
        "    ancilla_qubits=ancilla_qubits,\n",
        "    cost=\"gamma\",\n",
        "    method=\"windowed\",\n",
        "    seed=seed,\n",
        ")\n",
        "\n",
        "print(f\"Physical layout of payload and ancillas: {layout + ancilla_qubits}\")\n",
        "print(\"Checked circuit:\")\n",
        "checked[-1].circuit.draw(\"mpl\", fold=-1, idle_wires=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "step3-header",
      "metadata": {},
      "source": [
        "### Step 3: Execute using Qiskit primitives\n",
        "\n",
        "To make the effect of gate noise visible, increase the depth of the payload and sample a subset of its stabilizers. Each stabilizer is generally not qubit-wise commuting with the others, so a single set of checks is not valid for two different stabilizers. Rather than group stabilizers into commuting sets, find a good set of checks for each stabilizer independently. Sampling stabilizers uniformly at random gives an unbiased fidelity estimate.\n",
        "\n",
        "Build the deeper circuit and draw a random sample of its stabilizers.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "step3-stabilizers",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Sampled 20 stabilizers of a 12-qubit circuit with two-qubit depth 24: {ZXIIXZYYXIZZ, XXXYIIZYXIII, ...}\n"
          ]
        }
      ],
      "source": [
        "depth = 24\n",
        "num_stabilizers = 20\n",
        "num_shots = 1_000\n",
        "\n",
        "circuit = random_clifford_circuit(num_qubits, depth, rng)\n",
        "\n",
        "# Build the full stabilizer group, then sample from it uniformly at random\n",
        "circ_no_meas = circuit.remove_final_measurements(inplace=False)\n",
        "stabilizer_group = PauliList([Pauli(\"I\" * num_qubits)])\n",
        "for generator in (\n",
        "    Pauli(label) for label in Clifford(circ_no_meas).to_labels(mode=\"S\")\n",
        "):\n",
        "    stabilizer_group = stabilizer_group + stabilizer_group.compose(generator)\n",
        "\n",
        "keep = np.where(\n",
        "    stabilizer_group.x.any(axis=1) | stabilizer_group.z.any(axis=1)\n",
        ")[0]\n",
        "chosen = np.random.default_rng(seed).choice(\n",
        "    keep, size=min(num_stabilizers, len(keep)), replace=False\n",
        ")\n",
        "stabilizers = [stabilizer_group[int(i)] for i in chosen]\n",
        "\n",
        "two_qubit_depth = circuit.depth(lambda x: x.operation.num_qubits == 2)\n",
        "print(\n",
        "    f\"Sampled {len(stabilizers)} stabilizers of a {circuit.num_qubits}-qubit \"\n",
        "    f\"circuit with two-qubit depth {two_qubit_depth}: \"\n",
        "    f\"{{{stabilizers[0]}, {stabilizers[1]}, ...}}\"\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "step3-findchecks-md",
      "metadata": {},
      "source": [
        "For each sampled stabilizer, rotate the circuit so that the stabilizer is measured in the computational basis, transpile it onto the backend, and find a good set of checks. The target and ancilla pairs are shuffled together for each stabilizer so that each target keeps its ancilla. Remember that checks are committed sequentially in the order the target qubits are given, and a committed check is not changed as more checks are added.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "step3-findchecks",
      "metadata": {},
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "100%|██████████| 20/20 [00:15<00:00,  1.29it/s]"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Added 7 checks to 20 circuits in 15s.\n",
            "On average, two-qubit depth increased from 24 to 33 when adding 7 checks.\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "\n"
          ]
        }
      ],
      "source": [
        "noisy_circuits = []\n",
        "checked_circuits = []\n",
        "depths_2q = []\n",
        "t0 = time.time()\n",
        "for i, pauli in enumerate(tqdm(stabilizers)):\n",
        "    noisy_circuits.append(pm.run(append_basis_rotation(circuit, pauli)))\n",
        "    # Shuffle target and ancilla pairs together so each target keeps its ancilla\n",
        "    targets, ancillas = zip(\n",
        "        *random.sample(\n",
        "            list(zip(target_qubits, ancilla_qubits, strict=True)),\n",
        "            k=len(target_qubits),\n",
        "        ),\n",
        "        strict=True,\n",
        "    )\n",
        "    checked_circuits.append(\n",
        "        add_pauli_checks(\n",
        "            noisy_circuits[-1],\n",
        "            list(targets),\n",
        "            noise_model,\n",
        "            ancilla_qubits=list(ancillas),\n",
        "            cost=\"gamma\",\n",
        "            method=\"windowed\",\n",
        "            seed=seed + 1 + i,\n",
        "        )\n",
        "    )\n",
        "    depths_2q.append(\n",
        "        checked_circuits[-1][-1].circuit.depth(lambda x: len(x.qubits) == 2)\n",
        "    )\n",
        "\n",
        "print(\n",
        "    f\"Added {num_checks} checks to {len(stabilizers)} circuits \"\n",
        "    f\"in {(time.time() - t0):.0f}s.\"\n",
        ")\n",
        "print(\n",
        "    f\"On average, two-qubit depth increased from \"\n",
        "    f\"{circuit.depth(lambda x: len(x.qubits) == 2)} to {int(np.mean(depths_2q))} \"\n",
        "    f\"when adding {num_checks} checks.\"\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "step3-sample-md",
      "metadata": {},
      "source": [
        "Sample the bare payload and the checked circuits with Qiskit Aer. The simulator uses the same depolarizing model that scored the checks, so the noise that the checks target is the noise that the simulator applies.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "id": "step3-sample",
      "metadata": {},
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "100%|██████████| 20/20 [00:17<00:00,  1.13it/s]\n"
          ]
        }
      ],
      "source": [
        "aer_nm = AerNoiseModel()\n",
        "aer_nm.add_all_qubit_quantum_error(\n",
        "    depolarizing_error(noise_model.gate_noise, 2), [\"cz\"]\n",
        ")\n",
        "p = noise_model.readout_noise\n",
        "aer_nm.add_all_qubit_readout_error(ReadoutError([[1 - p, p], [p, 1 - p]]))\n",
        "noisy_sim = AerSimulator(method=\"stabilizer\", noise_model=aer_nm)\n",
        "\n",
        "counts = []\n",
        "for i, checked_circ_result in enumerate(tqdm(checked_circuits)):\n",
        "    noisy_counts = (\n",
        "        noisy_sim.run(\n",
        "            noisy_circuits[i], shots=num_shots, seed_simulator=seed * i + 1\n",
        "        )\n",
        "        .result()\n",
        "        .get_counts()\n",
        "    )\n",
        "    checked_counts_per_variant = []\n",
        "    for k, ck in enumerate(checked_circ_result):\n",
        "        variant_counts = (\n",
        "            noisy_sim.run(\n",
        "                ck.circuit, shots=num_shots, seed_simulator=seed * i + 2 + k\n",
        "            )\n",
        "            .result()\n",
        "            .get_counts()\n",
        "        )\n",
        "        checked_counts_per_variant.append(variant_counts)\n",
        "    counts.append((noisy_counts, checked_counts_per_variant))"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "step4-header",
      "metadata": {},
      "source": [
        "### Step 4: Post-process and return result in desired classical format\n",
        "\n",
        "Each check uses entangling gates between one ancilla and one target. The ancilla starts in $|0\\rangle$, so $Z_\\text{anc}$ stabilizes its input. Propagating $Z_\\text{anc}$ forward through the checked circuit yields a Pauli operator on the output whose non-identity terms define the check's support. A check passes when the bits in its support have even parity. A sample is kept only when every check passes.\n",
        "\n",
        "The `get_postselection_method` of each `CheckedCircuit` returns a function that maps a measured bitstring to a syndrome vector. Keep the samples whose syndrome is zero for every check, and discard the rest. The chart that follows shows that adding more checks lowers the postselection rate. A lower postselection rate needs more shots to reach a target accuracy, so there is a tradeoff between detection capability and sampling cost. The rate appears to converge, which indicates that additional checks contribute less detection capability.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "id": "step4-psr",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/spacetime-codes/extracted-outputs/step4-psr-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "rate_per_variant = []\n",
        "kept_per_stab = []\n",
        "for i, (_, checked_counts_per_variant) in enumerate(counts):\n",
        "    rates = []\n",
        "    kept_at_num_checks = None\n",
        "    for k, variant_counts in enumerate(checked_counts_per_variant):\n",
        "        ps_fn = checked_circuits[i][k].get_postselection_method()\n",
        "        kept = {\n",
        "            bs: n for bs, n in variant_counts.items() if not ps_fn(bs).any()\n",
        "        }\n",
        "        rates.append(sum(kept.values()) / num_shots)\n",
        "        if k == num_checks:\n",
        "            kept_at_num_checks = kept\n",
        "    rate_per_variant.append(rates)\n",
        "    kept_per_stab.append(kept_at_num_checks)\n",
        "\n",
        "max_len = max(len(s) for s in rate_per_variant)\n",
        "rates_arr = np.full((len(rate_per_variant), max_len), np.nan)\n",
        "for i, s in enumerate(rate_per_variant):\n",
        "    rates_arr[i, : len(s)] = s\n",
        "ks = np.arange(max_len)\n",
        "\n",
        "fig, ax = plt.subplots(figsize=(8, 4))\n",
        "ax.plot(ks, rates_arr.T, color=\"#ff8c00\", alpha=0.15, linewidth=1)\n",
        "ax.plot(\n",
        "    ks,\n",
        "    np.nanmedian(rates_arr, axis=0),\n",
        "    color=\"black\",\n",
        "    linewidth=1,\n",
        "    linestyle=\"--\",\n",
        "    label=\"median\",\n",
        ")\n",
        "ax.set_xlabel(\"Checks committed\")\n",
        "ax.set_ylabel(\"Postselection rate\")\n",
        "ax.set_ylim((0, 1.05))\n",
        "ax.set_title(\n",
        "    f\"Per-stabilizer postselection rate ({len(rates_arr)} stabilizers)\"\n",
        ")\n",
        "ax.legend()\n",
        "ax.grid(True, alpha=0.3)\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "step4-fidelity-md",
      "metadata": {},
      "source": [
        "Now compare the fidelity of the bare noisy state with the postselected state. Postselecting only the samples with no detected error raises the expectation value of every stabilizer, and therefore the estimated fidelity. The postselected values use fewer samples than the raw values, yet the expectation values are more accurate and the sampled variance is lower. Notice also that the mean postselection rate is close to the noisy fidelity. This is what you expect when the checks detect nearly all erroneous samples: the fraction of samples that pass every check approaches the fraction of error-free samples, which is the fidelity of the noisy state.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "id": "step4-fidelity",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "ideal fidelity:          1.0\n",
            "noisy fidelity:          0.7899\n",
            "postselected fidelity:   0.9679\n",
            "mean postselection rate: 0.780\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/spacetime-codes/extracted-outputs/step4-fidelity-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/spacetime-codes/extracted-outputs/step4-fidelity-2.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "results = []\n",
        "for i, ((noisy_counts, _), kept) in enumerate(\n",
        "    zip(counts, kept_per_stab, strict=True)\n",
        "):\n",
        "    results.append(\n",
        "        (\n",
        "            expectation(noisy_counts, stabilizers[i]),\n",
        "            expectation(kept, stabilizers[i]),\n",
        "            sum(kept.values()) / num_shots,\n",
        "        )\n",
        "    )\n",
        "\n",
        "fidelity_noisy = float(np.nanmean([r[0] for r in results]))\n",
        "fidelity_postsel = float(np.nanmean([r[1] for r in results]))\n",
        "psr = float(np.mean([r[2] for r in results]))\n",
        "print(\n",
        "    f\"ideal fidelity:          1.0\\n\"\n",
        "    f\"noisy fidelity:          {fidelity_noisy:.4f}\\n\"\n",
        "    f\"postselected fidelity:   {fidelity_postsel:.4f}\\n\"\n",
        "    f\"mean postselection rate: {psr:.3f}\"\n",
        ")\n",
        "\n",
        "evs_ideal = np.ones(len(results))\n",
        "evs_noisy = np.array([r[0] for r in results])\n",
        "evs_post = np.array([r[1] for r in results])\n",
        "idx = np.arange(len(results))\n",
        "\n",
        "\n",
        "def strip(ax, ys, color, label):\n",
        "    m, s = np.nanmean(ys), np.nanstd(ys)\n",
        "    ax.axhspan(\n",
        "        m - s, m + s, color=color, alpha=0.15, label=f\"{label} mean and std\"\n",
        "    )\n",
        "    ax.axhline(\n",
        "        m, color=color, linewidth=1, linestyle=\"--\", label=f\"{label} fidelity\"\n",
        "    )\n",
        "\n",
        "\n",
        "fig, ax = plt.subplots(figsize=(8, 4))\n",
        "ax.axhline(np.nanmean(evs_ideal), color=\"black\", linewidth=1.5, label=\"ideal\")\n",
        "strip(ax, evs_noisy, \"red\", \"noisy\")\n",
        "strip(ax, evs_post, \"green\", \"postselected\")\n",
        "ax.scatter(idx, evs_noisy, color=\"red\", s=22, alpha=0.7, label=\"noisy EVs\")\n",
        "ax.scatter(\n",
        "    idx,\n",
        "    evs_post,\n",
        "    color=\"green\",\n",
        "    s=22,\n",
        "    alpha=0.7,\n",
        "    label=\"postselected EVs\",\n",
        ")\n",
        "ax.set_xlabel(\"stabilizer index\")\n",
        "ax.set_ylabel(r\"$\\langle G \\rangle$\")\n",
        "ax.set_ylim((-0.1, 1.1))\n",
        "ax.set_title(\"Per-stabilizer expectation values\")\n",
        "ax.legend(loc=\"lower left\")\n",
        "ax.grid(True, alpha=0.3)\n",
        "plt.tight_layout()\n",
        "plt.show()\n",
        "\n",
        "M = np.arange(1, len(results) + 1)\n",
        "fig, ax = plt.subplots(figsize=(8, 4))\n",
        "for ys, color, label in [\n",
        "    (evs_ideal, \"black\", \"ideal\"),\n",
        "    (evs_noisy, \"red\", \"noisy\"),\n",
        "    (evs_post, \"green\", \"postselected\"),\n",
        "]:\n",
        "    cm, sem = cum_mean_sem(ys)\n",
        "    ax.plot(M, cm, color=color, linewidth=1.5, label=label)\n",
        "    ax.fill_between(M, cm - sem, cm + sem, color=color, alpha=0.15)\n",
        "ax.set_xlabel(\"number of stabilizers averaged\")\n",
        "ax.set_ylabel(\"running fidelity estimate\")\n",
        "ax.set_title(\"Fidelity convergence versus number of stabilizers\")\n",
        "ax.legend()\n",
        "ax.grid(True, alpha=0.3)\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "step4-gamma-md",
      "metadata": {},
      "source": [
        "The gamma score reports how much of the modeled noise channel remains undetected by the checks. Plotting the gamma score against the number of committed checks shows how the detection capability improves as each check is added. A value of `1.0` means that the checks capture all of the modeled noise. The curves fall toward `1.0` as more checks are committed, which shows that each additional check captures part of the remaining undetected error.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "id": "step4-gamma",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/spacetime-codes/extracted-outputs/step4-gamma-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "stab_scores = [\n",
        "    [variant.cost for variant in checked_circ_result]\n",
        "    for checked_circ_result in checked_circuits\n",
        "]\n",
        "max_len = max(len(s) for s in stab_scores)\n",
        "scores = np.full((len(stab_scores), max_len), np.nan)\n",
        "for i, s in enumerate(stab_scores):\n",
        "    scores[i, : len(s)] = s\n",
        "ks = np.arange(max_len)\n",
        "\n",
        "fig, ax = plt.subplots(figsize=(8, 4))\n",
        "ax.plot(ks, scores.T, color=\"#4682b4\", alpha=0.15, linewidth=1)\n",
        "ax.plot(\n",
        "    ks,\n",
        "    np.nanmedian(scores, axis=0),\n",
        "    color=\"black\",\n",
        "    linewidth=1,\n",
        "    linestyle=\"--\",\n",
        "    label=\"median\",\n",
        ")\n",
        "ax.set_xlabel(\"Checks committed\")\n",
        "ax.set_ylabel(\"Gamma\")\n",
        "ax.set_yscale(\"log\")\n",
        "ax.set_title(f\"Per-stabilizer gamma curves ({len(scores)} stabilizers)\")\n",
        "ax.legend()\n",
        "ax.grid(True, alpha=0.3, which=\"both\")\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "hardware-header",
      "metadata": {},
      "source": [
        "## Large-scale hardware example\n",
        "\n",
        "The same workflow runs on hardware with a larger, deeper payload. This section reuses the backend from the simulator example but builds a new 20-qubit layout with its own target and ancilla pairs and pass manager, then submits the circuits to the QPU in a single job. At this size, most shots trigger at least one check, so the postselection rate is low, and each circuit needs a large shot budget for enough samples to survive. The example therefore concentrates its budget on a few sampled stabilizers; this is still an unbiased fidelity estimate, but is coarser than the simulator example's average over many stabilizers.\n",
        "\n",
        "One thing changes relative to the simulator example: instead of inferring a uniform depolarizing channel from calibration data, this section **learns** the noise model with `NoiseLearner` and builds the `qiskit-paulice` model from the result with `NoiseModel.from_pauli_lindblad_maps`. A learned Pauli-Lindblad model captures the spatial structure of the noise on this specific layout, rather than assuming every edge is equally noisy, so the check placement is scored against noise more similar to that affecting the QPU. Learning the noise requires sampling the QPU and should be accounted for in any overall QPU sampling budget.\n",
        "\n",
        "The parameters that follow set the number of qubits, the depth, the number of stabilizers, and the number of shots. Scale `hw_num_shots` with the inverse of the postselection rate: at a 3% rate, 40,000 shots leave roughly 1,200 postselected samples per circuit. Increase `hw_num_stabilizers` for a more precise fidelity estimate at the cost of more circuits per job, each of which needs the same shot budget.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "hardware-steps-header",
      "metadata": {},
      "source": [
        "### Steps 1-4 (compressed into a single code block)\n",
        "\n",
        "The following cell runs the same four steps as the simulator example. It builds the larger payload and samples a few stabilizers (step 1); chooses a layout, learns the noise model on it, and finds the fully checked circuit for each stabilizer (step 2); submits one Sampler job that contains both the bare and the checked circuits (step 3); and postselects the checked counts to compare the noisy and postselected fidelity estimates, per stabilizer and on average (step 4). At this size, enumerating the full stabilizer group like the simulator example is infeasible, so the cell subsamples a random subset of stabilizers to calculate an estimate of the fidelity.\n",
        "\n",
        "Note that step 2 does more here than in the simulator example: learning the noise model submits its own `NoiseLearner` job before the Sampler job, so the cell runs **two** jobs in total. They carry the tags `TUT_ASPC_LEARN` and `TUT_ASPC` so that you can find them later. See [Organize and search by job tags](/docs/guides/add-job-tags) for more on tagging jobs.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "id": "hardware-submit",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Layout with 11 check pairs: [153, 152, 151, 138, 131, 130, 129, 118, 109, 110, 111, 98, 91, 90, 89, 78, 69, 70, 71, 58]\n",
            "Submitted noise-learner job d9f4mncjeosc73fjfmkg\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Learned 4 layers; readout noise 0.00470\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "100%|██████████| 10/10 [01:01<00:00,  6.18s/it]\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Submitted job d9f4s04jeosc73fjftkg with 20 circuits\n",
            "noisy fidelity estimate:        0.3685\n",
            "postselected fidelity estimate: 0.6869\n",
            "mean postselection rate:        0.2851 (~11404 kept shots per circuit)\n",
            "\n",
            "per-stabilizer results:\n",
            "idx     noisy   postsel      psr  kept shots\n",
            "  0    0.3769    0.6918   0.3247       12987\n",
            "  1    0.3745    0.6760   0.2999       11995\n",
            "  2    0.3659    0.6389   0.3549       14196\n",
            "  3    0.3821    0.7060   0.2660       10641\n",
            "  4    0.3653    0.7475   0.2531       10124\n",
            "  5    0.3752    0.7022   0.2698       10791\n",
            "  6    0.3508    0.7144   0.2711       10842\n",
            "  7    0.3485    0.7087   0.2381        9523\n",
            "  8    0.3825    0.6289   0.2928       11711\n",
            "  9    0.3630    0.6549   0.2808       11232\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/spacetime-codes/extracted-outputs/hardware-submit-5.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# -------------------------Step 1: build a larger payload and sample stabilizers-------------------------\n",
        "hw_num_qubits = 20\n",
        "hw_depth = 36\n",
        "hw_num_stabilizers = 10\n",
        "hw_num_shots = 40_000\n",
        "\n",
        "hw_circuit = random_clifford_circuit(hw_num_qubits, hw_depth, rng)\n",
        "hw_no_meas = hw_circuit.remove_final_measurements(inplace=False)\n",
        "\n",
        "# Enumerating all 2^n stabilizers is infeasible at this size, so draw each\n",
        "# stabilizer by composing a random subset of the group generators\n",
        "hw_generators = [\n",
        "    Pauli(label) for label in Clifford(hw_no_meas).to_labels(mode=\"S\")\n",
        "]\n",
        "sample_rng = np.random.default_rng(seed)\n",
        "hw_stabilizers = []\n",
        "while len(hw_stabilizers) < hw_num_stabilizers:\n",
        "    mask = sample_rng.integers(0, 2, hw_num_qubits).astype(bool)\n",
        "    if not mask.any():\n",
        "        continue  # skip the identity\n",
        "    stabilizer = Pauli(\"I\" * hw_num_qubits)\n",
        "    for generator, chosen in zip(hw_generators, mask, strict=True):\n",
        "        if chosen:\n",
        "            stabilizer = stabilizer.compose(generator)\n",
        "    hw_stabilizers.append(stabilizer)\n",
        "\n",
        "# -------------------------Step 2: find a 20-qubit layout, learn its noise, and add checks-------------------------\n",
        "# A single bad coupler or bad-readout qubit on the path drags every\n",
        "# stabilizer down, so search harder and with tighter error thresholds\n",
        "hw_layout = find_check_layout(\n",
        "    backend,\n",
        "    hw_num_qubits,\n",
        "    rng,\n",
        "    num_trials=500,\n",
        "    max_gate_error=0.015,\n",
        "    max_readout_error=0.05,\n",
        ")\n",
        "hw_target_qubits, hw_ancilla_qubits = get_check_qubits(backend, hw_layout)\n",
        "hw_pm = generate_preset_pass_manager(\n",
        "    optimization_level=0, backend=backend, initial_layout=hw_layout\n",
        ")\n",
        "print(f\"Layout with {len(hw_target_qubits)} check pairs: {hw_layout}\")\n",
        "\n",
        "# ----- learn a Pauli-Lindblad noise model on this layout -----\n",
        "# The simulator example scored checks against a uniform depolarizing channel\n",
        "# inferred from calibration data. Here, learn the noise instead: NoiseLearner\n",
        "# runs its own job on the QPU and returns a Pauli-Lindblad channel per unique\n",
        "# entangling layer, so the checks are placed against the noise this layout\n",
        "# actually has, including its spatial structure. All the sampled stabilizers\n",
        "# share the same entangling layers and differ only in their final basis\n",
        "# rotation, so learning on the bare payload covers all of them.\n",
        "learner = NoiseLearner(\n",
        "    mode=backend,\n",
        "    options={\n",
        "        \"max_layers_to_learn\": 4,\n",
        "        \"num_randomizations\": 32,\n",
        "        \"shots_per_randomization\": 128,\n",
        "        \"environment\": {\"job_tags\": [\"TUT_ASPC_LEARN\"]},\n",
        "    },\n",
        ")\n",
        "learner_job = learner.run([hw_pm.run(hw_circuit)])\n",
        "print(f\"Submitted noise-learner job {learner_job.job_id()}\")\n",
        "hw_layer_errors = learner_job.result().data\n",
        "\n",
        "# To see how much the learned model helps, swap the next line for the\n",
        "# simulator example's uniform model - a one-line change:\n",
        "# hw_noise_model = NoiseModel.from_backend(backend, hw_layout, uniform_gate_noise=True)\n",
        "hw_noise_model = learned_noise_model(hw_layer_errors, hw_layout)\n",
        "# NoiseLearner characterizes gate noise only, so keep the readout estimate\n",
        "# from calibration data rather than leaving it unset\n",
        "hw_noise_model.readout_noise = NoiseModel.from_backend(\n",
        "    backend, hw_layout, uniform_gate_noise=True\n",
        ").readout_noise\n",
        "print(\n",
        "    f\"Learned {len(hw_layer_errors)} layers; \"\n",
        "    f\"readout noise {hw_noise_model.readout_noise:.5f}\"\n",
        ")\n",
        "\n",
        "# ----- add the fully checked circuit per stabilizer -----\n",
        "hw_noisy_circuits = []\n",
        "hw_checked_circuits = []\n",
        "for i, pauli in enumerate(tqdm(hw_stabilizers)):\n",
        "    bare = hw_pm.run(append_basis_rotation(hw_circuit, pauli))\n",
        "    hw_noisy_circuits.append(bare)\n",
        "    variants = add_pauli_checks(\n",
        "        bare,\n",
        "        hw_target_qubits,\n",
        "        hw_noise_model,\n",
        "        ancilla_qubits=hw_ancilla_qubits,\n",
        "        cost=\"gamma\",\n",
        "        method=\"windowed\",\n",
        "        seed=seed + 1 + i,\n",
        "    )\n",
        "    hw_checked_circuits.append(variants[-1])  # keep the fully checked circuit\n",
        "\n",
        "# -------------------------Step 3: submit one Sampler job with the bare and checked circuits-------------------------\n",
        "sampler = Sampler(mode=backend)\n",
        "sampler.options.default_shots = hw_num_shots\n",
        "sampler.options.environment.job_tags = [\"TUT_ASPC\"]\n",
        "\n",
        "pubs = hw_noisy_circuits + [cc.circuit for cc in hw_checked_circuits]\n",
        "job = sampler.run(pubs)\n",
        "print(f\"Submitted job {job.job_id()} with {len(pubs)} circuits\")\n",
        "\n",
        "# -------------------------Step 4: postselect and compare fidelity-------------------------\n",
        "result = job.result()\n",
        "n_stab = len(hw_stabilizers)\n",
        "\n",
        "hw_results = []\n",
        "for i in range(n_stab):\n",
        "    noisy_counts = result[i].join_data().get_counts()\n",
        "    checked_counts = result[n_stab + i].join_data().get_counts()\n",
        "    ps_fn = hw_checked_circuits[i].get_postselection_method()\n",
        "    kept = {bs: c for bs, c in checked_counts.items() if not ps_fn(bs).any()}\n",
        "    hw_results.append(\n",
        "        (\n",
        "            expectation(noisy_counts, hw_stabilizers[i]),\n",
        "            expectation(kept, hw_stabilizers[i]),\n",
        "            sum(kept.values()) / sum(checked_counts.values()),\n",
        "        )\n",
        "    )\n",
        "\n",
        "hw_fidelity_noisy = float(np.nanmean([r[0] for r in hw_results]))\n",
        "hw_fidelity_postsel = float(np.nanmean([r[1] for r in hw_results]))\n",
        "hw_psr = float(np.mean([r[2] for r in hw_results]))\n",
        "print(\n",
        "    f\"noisy fidelity estimate:        {hw_fidelity_noisy:.4f}\\n\"\n",
        "    f\"postselected fidelity estimate: {hw_fidelity_postsel:.4f}\\n\"\n",
        "    f\"mean postselection rate:        {hw_psr:.4f} \"\n",
        "    f\"(~{int(round(hw_psr * hw_num_shots))} kept shots per circuit)\"\n",
        ")\n",
        "\n",
        "# Per-stabilizer breakdown. The postselection rate varies from stabilizer to\n",
        "# stabilizer, so a stabilizer whose postselected value barely moves is usually\n",
        "# one whose checks rejected little; the kept-shot count says how much of the\n",
        "# gap is statistics rather than signal.\n",
        "print(\"\\nper-stabilizer results:\")\n",
        "print(\n",
        "    f\"{'idx':>3}  {'noisy':>8}  {'postsel':>8}  {'psr':>7}  {'kept shots':>10}\"\n",
        ")\n",
        "for i, (noisy, post, psr_i) in enumerate(hw_results):\n",
        "    print(\n",
        "        f\"{i:>3}  {noisy:>8.4f}  {post:>8.4f}  {psr_i:>7.4f}  \"\n",
        "        f\"{int(round(psr_i * hw_num_shots)):>10}\"\n",
        "    )\n",
        "\n",
        "hw_noisy = np.array([r[0] for r in hw_results])\n",
        "hw_post = np.array([r[1] for r in hw_results])\n",
        "idx = np.arange(n_stab)\n",
        "\n",
        "fig, ax = plt.subplots(figsize=(9, 4))\n",
        "ax.axhline(1.0, color=\"black\", linewidth=1.5, label=\"ideal\")\n",
        "strip(ax, hw_noisy, \"red\", \"noisy\")\n",
        "strip(ax, hw_post, \"green\", \"postselected\")\n",
        "ax.scatter(idx, hw_noisy, color=\"red\", s=22, alpha=0.7, label=\"noisy EVs\")\n",
        "ax.scatter(\n",
        "    idx,\n",
        "    hw_post,\n",
        "    color=\"green\",\n",
        "    s=22,\n",
        "    alpha=0.7,\n",
        "    label=\"postselected EVs\",\n",
        ")\n",
        "ax.set_xlabel(\"stabilizer index\")\n",
        "ax.set_ylabel(r\"$\\langle G \\rangle$\")\n",
        "ax.set_ylim((-0.1, 1.1))\n",
        "ax.set_xticks(idx)\n",
        "ax.set_title(\"Per-stabilizer expectation values on hardware\")\n",
        "# Outside the axes so it cannot hide a data point\n",
        "ax.legend(loc=\"center left\", bbox_to_anchor=(1.02, 0.5), frameon=False)\n",
        "ax.grid(True, alpha=0.3)\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "hardware-results-md",
      "metadata": {},
      "source": [
        "For circuits of this size, most samples contain at least one detected error, so the postselection rate is small and postselection discards most shots. The samples that pass every check give a much better expectation value than the bare circuit, and the per-stabilizer values separate cleanly from the noisy baseline. To tighten the fidelity estimate, sample more stabilizers with the same per-circuit shot budget. To raise the postselection rate, reduce the circuit depth or commit fewer checks; to push to larger payloads, scale the shot budget with the inverse of the postselection rate.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "next-steps",
      "metadata": {},
      "source": [
        "## Next steps\n",
        "\n",
        "<Admonition type=\"tip\" title=\"Recommendations\">\n",
        "  If you found this work interesting, you might be interested in the following material:\n",
        "\n",
        "  * The tutorial on [repetition codes](/docs/tutorials/repetition-codes) for an introduction to quantum error correction.\n",
        "  * The [`qiskit-paulice` documentation](/docs/addons/qiskit-addon-paulice) for the full check-finding API, and the package's [GitHub repository](https://github.com/Qiskit/qiskit-paulice) for the source code.\n",
        "  * The paper [Low-overhead error detection with spacetime codes](https://arxiv.org/abs/2504.15725) for the theory behind the checks.\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "references",
      "metadata": {},
      "source": [
        "## References\n",
        "\n",
        "* \\[1] Martiel, S., & Javadi-Abhari, A. (2025). Low-overhead error detection with spacetime codes. *arXiv preprint* [arXiv:2504.15725](https://arxiv.org/abs/2504.15725).\n",
        "* \\[2] van den Berg, E., Bravyi, S., Gambetta, J. M., Jurcevic, P., Maslov, D., & Temme, K. (2023). Single-shot error mitigation by coherent Pauli checks. *Physical Review Research*, 5(3), 033193.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "id": "a1b8767d",
      "source": "© IBM Corp., 2017-2026"
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3"
    },
    "hours": 2.5,
    "qpuSeconds": 10
  },
  "nbformat": 4,
  "nbformat_minor": 5
}