{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "59c663c5",
      "metadata": {},
      "source": [
        "---\n",
        "title: Low-overhead error detection with spacetime codes\n",
        "description: Experiment with error detection on a small-scale random Clifford circuit, then walk through the task of GHZ state preparation.\n",
        "---\n",
        "\n",
        "{/* cspell:ignore unflagged Tunables parallelizable uncomputation  */}\n",
        "\n",
        "# Low-overhead error detection with spacetime codes\n",
        "\n",
        "*Usage estimate: 10 seconds on a Heron r3 processor (NOTE: This is an estimate only. Your runtime might vary.)*\n",
        "\n",
        "## Introduction\n",
        "\n",
        "[Low-overhead error detection with spacetime codes](https://arxiv.org/abs/2504.15725) [\\[1\\]](#references) by Simon Martiel and Ali Javadi-Abhari proposes synthesizing low-weight, connectivity-aware spacetime checks for Clifford-dominated circuits, then post-selecting on those checks to catch faults with far less overhead than full error correction and fewer shots than standard error mitigation.\n",
        "\n",
        "This paper proposes a novel method for error detection in quantum circuits (specifically Clifford circuits) that strikes a balance between full error correction and lighter-weight mitigation techniques. The key idea is to use spacetime codes to generate \"checks\" across the circuit that are capable of catching errors, with significantly lower qubit and gate overhead than full fault-tolerant error correction. The authors design efficient algorithms to select checks that are low-weight (involving few qubits), compatible with the physical connectivity of the device, and cover large temporal and spatial regions of the circuit. They demonstrate the approach on circuits with up to 50 logical qubits and \\~2450 CZ gates, achieving physical-to-logical fidelity gains of up to 236x. Also note that as circuits include more non-Clifford operations, the number of valid checks shrinks exponentially, indicating the method works best for Clifford-dominated circuits. Overall, in the near term, error detection via spacetime codes may offer a practical, lower-overhead route to improving reliability in quantum hardware.\n",
        "\n",
        "This error detection technique relies on the notion of coherent Pauli checks and is based on the work [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 et al.\n",
        "\n",
        "More recently, the paper [Big cats: entanglement in 120 qubits and beyond](https://arxiv.org/abs/2510.09520) [\\[3\\]](#references) by Javadi-Abhari et al., reports the creation of a 120-qubit Greenberger-Horne-Zeilinger (GHZ) state, the largest multipartite entangled state achieved to date on a superconducting-qubit platform. Using a hardware-aware compiler, low-overhead error detection, and a \"temporary uncomputation\" technique to reduce noise, the researchers achieved a fidelity of 0.56 ± 0.03 with about 28% post-selection efficiency. The work demonstrates genuine entanglement across all 120 qubits, validating multiple fidelity-certification methods, and marks a major benchmark for scalable quantum hardware.\n",
        "\n",
        "This tutorial builds on these ideas, guiding you through implementing the error detection algorithm first on a small-scale random Clifford circuit and then through the task of GHZ state preparation, to help you experiment with error detection on your own quantum circuits.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "fce0de6c",
      "metadata": {},
      "source": [
        "### Requirements\n",
        "\n",
        "Before starting this tutorial, ensure that 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.2 (`pip install qiskit-aer`)\n",
        "* [Qiskit Device Benchmarking](https://github.com/qiskit-community/qiskit-device-benchmarking) (`pip install \"qiskit-device-benchmarking @ git+https://github.com/qiskit-community/qiskit-device-benchmarking.git\"`)\n",
        "* NumPy v2.3.2 (`pip install numpy`)\n",
        "* Matplotlib v3.10.7 (`pip install matplotlib`)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4286334d",
      "metadata": {},
      "source": [
        "### Setup\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "8fdec072",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Standard library imports\n",
        "from collections import defaultdict, deque\n",
        "from functools import partial\n",
        "\n",
        "# External libraries\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "\n",
        "# Qiskit\n",
        "from qiskit import ClassicalRegister, QuantumCircuit\n",
        "from qiskit.circuit import Delay\n",
        "from qiskit.circuit.library import RZGate, XGate\n",
        "from qiskit.converters import circuit_to_dag, dag_to_circuit\n",
        "from qiskit.quantum_info import Pauli, random_clifford\n",
        "from qiskit.transpiler import AnalysisPass, PassManager\n",
        "from qiskit.transpiler.passes import (\n",
        "    ALAPScheduleAnalysis,\n",
        "    CollectAndCollapse,\n",
        "    PadDelay,\n",
        "    PadDynamicalDecoupling,\n",
        "    RemoveBarriers,\n",
        ")\n",
        "from qiskit.transpiler.passes.optimization.collect_and_collapse import (\n",
        "    collect_using_filter_function,\n",
        "    collapse_to_operation,\n",
        ")\n",
        "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
        "from qiskit.visualization import plot_histogram\n",
        "\n",
        "# Qiskit Aer\n",
        "from qiskit_aer import AerSimulator\n",
        "from qiskit_aer.noise import NoiseModel, ReadoutError, depolarizing_error\n",
        "\n",
        "# Qiskit IBM Runtime\n",
        "from qiskit_ibm_runtime import QiskitRuntimeService\n",
        "from qiskit_ibm_runtime import SamplerV2 as Sampler\n",
        "\n",
        "# Qiskit Device Benchmarking\n",
        "from qiskit_device_benchmarking.utilities.gate_map import plot_gate_map"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "3888d8aa",
      "metadata": {},
      "source": [
        "## First example\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "704eb75f",
      "metadata": {},
      "source": [
        "To demonstrate this method, we start by constructing a simple Clifford circuit. Our goal is to be able to detect when certain types of error occur in this circuit, so that we can discard erroneous measurement results. In error detection terminology, this is also known as our payload circuit.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "f8384a7d",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/error-detection/extracted-outputs/f8384a7d-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 2,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "circ = random_clifford(num_qubits=2, seed=11).to_circuit()\n",
        "circ.draw(\"mpl\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "aa41f799",
      "metadata": {},
      "source": [
        "Our goal is to insert a coherent Pauli check into this payload circuit. But before we do that, we separate this circuit into layers. This will be useful later when inserting Pauli gates in between.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "04bb6ea7",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/error-detection/extracted-outputs/04bb6ea7-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Separate circuit into layers\n",
        "dag = circuit_to_dag(circ)\n",
        "circ_layers = []\n",
        "for layer in dag.layers():\n",
        "    layer_as_circuit = dag_to_circuit(layer[\"graph\"])\n",
        "    circ_layers.append(layer_as_circuit)\n",
        "\n",
        "# Create subplots\n",
        "fig, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(1, 5, figsize=(10, 4))\n",
        "\n",
        "# Draw circuits on respective axes\n",
        "circ_layers[0].draw(output=\"mpl\", ax=ax1)\n",
        "circ_layers[1].draw(output=\"mpl\", ax=ax2)\n",
        "circ_layers[2].draw(output=\"mpl\", ax=ax3)\n",
        "circ_layers[3].draw(output=\"mpl\", ax=ax4)\n",
        "circ_layers[4].draw(output=\"mpl\", ax=ax5)\n",
        "\n",
        "# Adjust layout to prevent overlap\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "dda1b961",
      "metadata": {},
      "source": [
        "Now we are ready to add coherent Pauli checks into the payload circuit. To do this, we need to construct a \"valid check\" and insert it into the circuit. A \"check\" in this case is an operator that is capable of signaling whether an error occurred in the circuit by making a measurement on an ancilla qubit. It is considered a valid check when the additional operators inserted into the quantum circuit don't logically change the original circuit.\n",
        "\n",
        "This check is capable of detecting types of errors that anticommute with it, and the check will trigger a measurement of $\\ket{1}$ state in the ancilla qubit instead of $\\ket{0}$ through phase kickback. Therefore, we will be able to discard measurements where an error was signaled.\n",
        "\n",
        "In general, coherent Pauli checks are controlled-Pauli operators inserted into \"wires\" - spacetime locations between gates. The ancilla qubit responsible for signaling the error is the control qubit.\n",
        "\n",
        "Below we construct a valid check for the Clifford circuit we created above. We can demonstrate that this check doesn't change the circuit operation by showing that when these Pauli checks are propagated to the front of the circuit, they cancel each other out. This is easily shown because a Pauli operator through a Clifford gate is another Pauli operator.\n",
        "\n",
        "In general, one can use a decoding heuristic as outlined in [\\[1\\]](https://arxiv.org/abs/2504.15725) to identify valid checks. For the purposes of our initial example, we can also construct valid checks by using analytical Pauli and Clifford gate multiplication conditions.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "ca5727c1",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Define a valid check\n",
        "pauli_1 = Pauli(\"ZI\")\n",
        "pauli_2 = Pauli(\"XZ\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "e1c1e16d",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/error-detection/extracted-outputs/e1c1e16d-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 5,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "circ_1 = circ_layers[0].compose(circ_layers[1])\n",
        "circ_1.draw(\"mpl\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "bebdd7df",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "Pauli('-ZI')"
            ]
          },
          "execution_count": 6,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "pauli_1_ev = pauli_1.evolve(circ_1, frame=\"h\")\n",
        "pauli_1_ev"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "246a42b8",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/error-detection/extracted-outputs/246a42b8-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 7,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "circ_2 = circ.copy()\n",
        "circ_2.draw(\"mpl\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "415c9471",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "Pauli('-ZI')"
            ]
          },
          "execution_count": 8,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "pauli_2_ev = pauli_2.evolve(circ_2, frame=\"h\")\n",
        "pauli_2_ev"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "f7368e97",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "Pauli('II')"
            ]
          },
          "execution_count": 9,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "pauli_1_ev.dot(pauli_2_ev)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "dba1bc80",
      "metadata": {},
      "source": [
        "As we can see, we have a valid check, since the inserted Pauli operators simply have the same effect as an identity operator on the circuit. We can now insert these checks into the circuit with an ancilla qubit. This ancilla qubit, or the check qubit, starts the in the $\\ket{+}$ state. It includes the controlled versions of the Pauli operations outlined above and is finally measured in the $X$ basis. This check qubit is now capable of capturing errors in the payload circuit without logically altering it. This is because certain types of noise in the payload circuit will modify the state of the check qubit, and it will be measured \"1\" instead of \"0\" in case such an error occurs.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "id": "52add690",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/error-detection/extracted-outputs/52add690-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 10,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# New circuit with 3 qubits (2 payload + 1 ancilla for check)\n",
        "circ_meas = QuantumCircuit(3)\n",
        "circ_meas.h(0)\n",
        "circ_meas.compose(circ_layers[0], [1, 2], inplace=True)\n",
        "circ_meas.compose(circ_layers[1], [1, 2], inplace=True)\n",
        "circ_meas.cz(0, 2)\n",
        "circ_meas.compose(circ_layers[2], [1, 2], inplace=True)\n",
        "circ_meas.compose(circ_layers[3], [1, 2], inplace=True)\n",
        "circ_meas.compose(circ_layers[4], [1, 2], inplace=True)\n",
        "circ_meas.cz(0, 1)\n",
        "circ_meas.cx(0, 2)\n",
        "circ_meas.h(0)\n",
        "\n",
        "# Add measurement to payload qubits\n",
        "c0 = ClassicalRegister(2, name=\"c0\")\n",
        "circ_meas.add_register(c0)\n",
        "circ_meas.measure(1, c0[0])\n",
        "circ_meas.measure(2, c0[1])\n",
        "\n",
        "# Add measurement to check qubit\n",
        "c1 = ClassicalRegister(1, name=\"c1\")\n",
        "circ_meas.add_register(c1)\n",
        "circ_meas.measure(0, c1[0])\n",
        "\n",
        "# Visualize the final circuit with the inserted checks\n",
        "circ_meas.draw(\"mpl\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ece387e6",
      "metadata": {},
      "source": [
        "If the check qubit is measured in \"0\", we keep that measurement. If it's measured in \"1\", then this means that that an error occurred in the payload circuit, and we discard that measurement.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "id": "78679785",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Stabilizer simulation result: {'0 11': 523, '0 01': 477}\n"
          ]
        }
      ],
      "source": [
        "# Noiseless simulation using stabilizer method\n",
        "sim_stab = AerSimulator(method=\"stabilizer\")\n",
        "res = sim_stab.run(circ_meas, shots=1000).result()\n",
        "counts_noiseless = res.get_counts()\n",
        "print(f\"Stabilizer simulation result: {counts_noiseless}\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "id": "4741910b",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/error-detection/extracted-outputs/4741910b-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 12,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Plot the noiseless results\n",
        "# Note that the first bit in the key corresponds to the check qubit\n",
        "plot_histogram(counts_noiseless)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2bbbdbc4",
      "metadata": {},
      "source": [
        "Note that with an ideal simulator, the check qubit doesn't detect any errors. We now introduce a noise model to the simulation and see how the check qubit captures errors.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "407c2a8a",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Noise model simulation result: {'1 01': 5, '0 11': 478, '1 11': 6, '1 00': 2, '1 10': 1, '0 01': 500, '0 00': 5, '0 10': 3}\n"
          ]
        }
      ],
      "source": [
        "# Qiskit Aer noise model\n",
        "noise = NoiseModel()\n",
        "p2 = 0.003  # two-qubit depolarizing per CZ\n",
        "p1 = 0.001  # one-qubit depolarizing per 1q Clifford\n",
        "pr = 0.01  # readout bit-flip probability\n",
        "\n",
        "# 1q depolarizing on common 1q gates\n",
        "e1 = depolarizing_error(p1, 1)\n",
        "for g1 in [\"id\", \"rz\", \"sx\", \"x\", \"h\", \"s\"]:\n",
        "    noise.add_all_qubit_quantum_error(e1, g1)\n",
        "\n",
        "# 2q depolarizing on CZ\n",
        "e2 = depolarizing_error(p2, 2)\n",
        "noise.add_all_qubit_quantum_error(e2, \"cz\")\n",
        "\n",
        "# Readout error on measure\n",
        "ro = ReadoutError([[1 - pr, pr], [pr, 1 - pr]])\n",
        "noise.add_all_qubit_readout_error(ro)\n",
        "\n",
        "# Qiskit Aer simulation with noise model\n",
        "aer = AerSimulator(method=\"automatic\", seed_simulator=43210)\n",
        "job = aer.run(circ_meas, shots=1000, noise_model=noise)\n",
        "result = job.result()\n",
        "counts_noisy = result.get_counts()\n",
        "print(f\"Noise model simulation result: {counts_noisy}\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "id": "51a8f3d9",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/error-detection/extracted-outputs/51a8f3d9-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 14,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "plot_histogram(counts_noisy)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "e7d863b1",
      "metadata": {},
      "source": [
        "As we can see, some measurements captured the error by flagging the check qubit as \"1\", which are visible in the last four columns. These shots are discarded.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "dd4baf08",
      "metadata": {},
      "source": [
        "Note: The ancilla qubit can also introduce new errors to the circuit. To reduce the effect of this, we can insert nested checks with additional ancilla qubits to the quantum circuit.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "7027e13d",
      "metadata": {},
      "source": [
        "## Real-world example: Prepare a GHZ state on real hardware\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c2b7f050",
      "metadata": {},
      "source": [
        "### Step 1: Map classical inputs to a quantum problem\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "96ee4f90",
      "metadata": {},
      "source": [
        "Now we demonstrate a significant task for quantum computing algorithms, which is preparing a GHZ state. We will demonstrate how to do this on a real backend using error detection.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "87317d7d",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Set optional seed for reproducibility\n",
        "SEED = 1\n",
        "\n",
        "if SEED:\n",
        "    np.random.seed(SEED)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0eca54de",
      "metadata": {},
      "source": [
        "The error detection algorithm for GHZ state preparation respects the hardware topology. We begin by selecting the desired hardware.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "78cafb7c",
      "metadata": {},
      "outputs": [],
      "source": [
        "# This is used to run on real hardware\n",
        "service = QiskitRuntimeService()\n",
        "\n",
        "# Choose a backend to build GHZ on\n",
        "backend_name = service.least_busy(\n",
        "    operational=True, simulator=False, min_num_qubits=133\n",
        ")\n",
        "\n",
        "backend = service.backend(backend_name)\n",
        "coupling_map = backend.target.build_coupling_map()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "12188865",
      "metadata": {},
      "source": [
        "A GHZ state on $n$ qubits is defined as\n",
        "$\\lvert \\mathrm{GHZ}_n\\rangle \\;=\\; \\frac{1}{\\sqrt{2}}\\Big(\\lvert 0\\rangle^{\\otimes n} \\,+\\, \\lvert 1\\rangle^{\\otimes n}\\Big).$\n",
        "\n",
        "A very naive approach to prepare the GHZ state would be to choose a root qubit with an initial Hadamard gate, which puts the qubit to an equal superposition state, and then to entangle this qubit with every other qubit. This is not a good approach, since it requires long-range and deep CNOT interactions. In this tutorial, we will use multiple techniques alongside error detection to reliably prepare the GHZ state on real hardware.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9dfcf756",
      "metadata": {},
      "source": [
        "### Step 2: Optimize problem for quantum hardware execution\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "cbb7cd7f",
      "metadata": {},
      "source": [
        "#### Map the GHZ state to hardware\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "181ce021",
      "metadata": {},
      "source": [
        "First, we search for a root to map the GHZ circuit on hardware. We remove edges/nodes whose CZ errors, measurement errors, and the T2 values are worse than the thresholds below. These will not be included in the GHZ circuit.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "id": "5033c96c",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "17 bad edges: \n",
            "[[30, 31], [112, 113], [113, 114], [113, 119], [120, 121], [130, 131], [145, 146], [146, 147], [111, 112], [55, 59], [64, 65], [131, 138], [131, 132], [119, 133], [129, 130], [47, 57], [29, 38]]\n",
            "5 bad nodes: \n",
            "[1, 113, 131, 146, 120]\n"
          ]
        }
      ],
      "source": [
        "def bad_cz(target, threshold=0.01):\n",
        "    \"\"\"Return list of edges whose CZ error is worse than threshold.\"\"\"\n",
        "    undirected_edges = []\n",
        "    for edge in backend.target.build_coupling_map().get_edges():\n",
        "        if (edge[1], edge[0]) not in undirected_edges:\n",
        "            undirected_edges.append(edge)\n",
        "    edges = undirected_edges\n",
        "    cz_errors = {}\n",
        "    for edge in edges:\n",
        "        cz_errors[edge] = target[\"cz\"][edge].error\n",
        "    worst_edges = sorted(cz_errors.items(), key=lambda x: x[1], reverse=True)\n",
        "    return [list(edge) for edge, error in worst_edges if error > threshold]\n",
        "\n",
        "\n",
        "def bad_readout(target, threshold=0.01):\n",
        "    \"\"\"Return list of nodes whose measurement error is worse than threshold.\"\"\"\n",
        "    meas_errors = {}\n",
        "    for node in range(backend.num_qubits):\n",
        "        meas_errors[node] = target[\"measure\"][(node,)].error\n",
        "    worst_nodes = sorted(\n",
        "        meas_errors.items(), key=lambda x: x[1], reverse=True\n",
        "    )\n",
        "    return [node for node, error in worst_nodes if error > threshold]\n",
        "\n",
        "\n",
        "def bad_coherence(target, threshold=60):\n",
        "    \"\"\"Return list of nodes whose T2 value is lower than threshold.\"\"\"\n",
        "    t2s = {}\n",
        "    for node in range(backend.num_qubits):\n",
        "        t2 = target.qubit_properties[node].t2\n",
        "        t2s[node] = t2 * 1e6 if t2 else 0\n",
        "    worst_nodes = sorted(t2s.items(), key=lambda x: x[1])\n",
        "    return [node for node, val in worst_nodes if val < threshold]\n",
        "\n",
        "\n",
        "THRESH_CZ = 0.025  # exclude from BFS those edges whose CZ error is worse than this threshold\n",
        "THRESH_MEAS = 0.15  # exclude from BFS those nodes whose measurement error is worse than this threshold\n",
        "THRESH_T2 = 10  # exclude from BFS those nodes whose T2 value is lower than this threshold\n",
        "\n",
        "bad_edges = bad_cz(backend.target, threshold=THRESH_CZ)\n",
        "bad_nodes_readout = bad_readout(backend.target, threshold=THRESH_MEAS)\n",
        "dead_qubits = bad_readout(backend.target, threshold=0.4)\n",
        "bad_nodes_coherence = bad_coherence(backend.target, threshold=THRESH_T2)\n",
        "bad_nodes = list(set(bad_nodes_readout) | set(bad_nodes_coherence))\n",
        "print(f\"{len(bad_edges)} bad edges: \\n{bad_edges}\")\n",
        "print(f\"{len(bad_nodes)} bad nodes: \\n{bad_nodes}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "637c4584",
      "metadata": {},
      "source": [
        "Using the function below, we construct the GHZ circuit on the chosen hardware starting from the root and using breadth-first search (BFS).\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "id": "5046c8ac",
      "metadata": {},
      "outputs": [],
      "source": [
        "def parallel_ghz(root, num_qubits, backend, bad_edges, skip):\n",
        "    \"\"\"\n",
        "    Build a GHZ state of size `num_qubits` on the given `backend`,\n",
        "    starting from `root`, expanding in BFS order.\n",
        "\n",
        "    At each BFS layer, every active qubit adds at most one new neighbor\n",
        "    (so that two-qubit operations can run in parallel with no qubit conflicts).\n",
        "\n",
        "    It grows the entanglement tree outward layer-by-layer.\n",
        "    \"\"\"\n",
        "\n",
        "    # -------------------------------------------------------------\n",
        "    # (1) Filter usable connections from the backend coupling map\n",
        "    # -------------------------------------------------------------\n",
        "    # The coupling map lists all directed hardware connections as (control, target).\n",
        "    # We remove edges that are:\n",
        "    #   - listed in `bad_edges` (or their reversed form)\n",
        "    #   - involve a qubit in the `skip` list\n",
        "    cmap = backend.configuration().coupling_map\n",
        "    edges = [\n",
        "        e\n",
        "        for e in cmap\n",
        "        if e not in bad_edges\n",
        "        and [e[1], e[0]] not in bad_edges\n",
        "        and e[0] not in skip\n",
        "        and e[1] not in skip\n",
        "    ]\n",
        "\n",
        "    # -------------------------------------------------------------\n",
        "    # (2) Build an undirected adjacency list for traversal\n",
        "    # -------------------------------------------------------------\n",
        "    # Even though coupling_map edges are directed, BFS expansion just needs\n",
        "    # connectivity information (so we treat edges as undirected for search).\n",
        "    adj = defaultdict(list)\n",
        "    for u, v in edges:\n",
        "        adj[u].append(v)\n",
        "        adj[v].append(u)\n",
        "\n",
        "    # -------------------------------------------------------------\n",
        "    # (3) Initialize the quantum circuit and BFS state\n",
        "    # -------------------------------------------------------------\n",
        "    n = backend.configuration().num_qubits\n",
        "    qc = QuantumCircuit(\n",
        "        n\n",
        "    )  # create a circuit with same number of qubits as hardware\n",
        "    visited = [\n",
        "        root\n",
        "    ]  # record the order qubits are added to the GHZ chain/tree\n",
        "    queue = deque([root])  # BFS queue (start from root)\n",
        "    explored = defaultdict(\n",
        "        set\n",
        "    )  # to track which neighbors each node has already explored\n",
        "    layers = []  # list of per-layer (control, target) gate tuples\n",
        "    qc.h(root)  # GHZ states start with a Hadamard on the root qubit\n",
        "\n",
        "    # -------------------------------------------------------------\n",
        "    # (4) BFS expansion: build the GHZ tree one layer at a time\n",
        "    # -------------------------------------------------------------\n",
        "    # Loop until we’ve added the desired number of qubits to the GHZ\n",
        "    while queue and len(visited) < num_qubits:\n",
        "        layer = []  # collect new (control, target) pairs for this layer\n",
        "        current = list(\n",
        "            queue\n",
        "        )  # snapshot current frontier (so queue mutations don’t affect iteration)\n",
        "        busy = (\n",
        "            set()\n",
        "        )  # track qubits already used in this layer (to avoid conflicts)\n",
        "\n",
        "        for node in current:\n",
        "            queue.popleft()\n",
        "\n",
        "            # find one unvisited neighbor of this node not already explored\n",
        "            unvisited_neighbors = [\n",
        "                nb\n",
        "                for nb in adj[node]\n",
        "                if nb not in visited and nb not in explored[node]\n",
        "            ]\n",
        "\n",
        "            if unvisited_neighbors:\n",
        "                nb = unvisited_neighbors[\n",
        "                    0\n",
        "                ]  # pick the first available neighbor\n",
        "                visited.append(nb)  # mark it as part of the GHZ structure\n",
        "                queue.append(\n",
        "                    node\n",
        "                )  # re-enqueue current node (can keep growing)\n",
        "                queue.append(nb)  # enqueue the newly added qubit\n",
        "                explored[node].add(nb)  # mark that edge as explored\n",
        "                layer.append(\n",
        "                    (node, nb)\n",
        "                )  # schedule a CNOT between node and neighbor\n",
        "                busy.update([node, nb])  # reserve both qubits for this layer\n",
        "\n",
        "                # stop early if we've reached the desired number of qubits\n",
        "                if len(visited) == num_qubits:\n",
        "                    break\n",
        "            # else: node has no unused unvisited neighbors left → skip\n",
        "\n",
        "        if layer:\n",
        "            # add all pairs (node, nb) scheduled this round to layers\n",
        "            layers.append(layer)\n",
        "        else:\n",
        "            # nothing new discovered this pass → done\n",
        "            break\n",
        "\n",
        "    # -------------------------------------------------------------\n",
        "    # (5) Emit all layers into the quantum circuit\n",
        "    # -------------------------------------------------------------\n",
        "    # For each layer:\n",
        "    #   - apply a CX gate for every (control, target) pair\n",
        "    #   - insert a barrier so transpiler keeps layer structure\n",
        "    for layer in layers:\n",
        "        for q1, q2 in layer:\n",
        "            qc.cx(q1, q2)\n",
        "        qc.barrier()\n",
        "\n",
        "    # -------------------------------------------------------------\n",
        "    # (6) Return outputs\n",
        "    # -------------------------------------------------------------\n",
        "    # qc: the built quantum circuit\n",
        "    # visited: order of qubits added\n",
        "    # layers: list of parallelizable two-qubit operations per step\n",
        "    return qc, visited, layers"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0167c0f2",
      "metadata": {},
      "source": [
        "We now repeatedly search for the best root, where the GHZ circuit will originate from.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "id": "90787353",
      "metadata": {},
      "outputs": [],
      "source": [
        "ROOT = None  # root for BFS search\n",
        "GHZ_SIZE = 100  # number of (data) qubits in the GHZ state\n",
        "SKIP = []  # nodes to intentionally skip so that we have a better chance for finding checks\n",
        "\n",
        "# Search for the best root (yielding the shallowest GHZ)\n",
        "if ROOT is None:\n",
        "    best_root = -1\n",
        "    base_depth = 100\n",
        "    for root in range(backend.num_qubits):\n",
        "        qc, ghz_qubits, _ = parallel_ghz(\n",
        "            root, GHZ_SIZE, backend, bad_edges, SKIP\n",
        "        )\n",
        "        if len(ghz_qubits) != GHZ_SIZE:\n",
        "            continue\n",
        "        depth = qc.depth(lambda x: x.operation.num_qubits == 2)\n",
        "        if depth < base_depth:\n",
        "            best_root = root\n",
        "            base_depth = depth\n",
        "    ROOT = best_root"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9b8d8394",
      "metadata": {},
      "source": [
        "We now construct the GHZ circuit starting from a specific node - as in, the best root - searching for the shortest depth using breadth-first search.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "id": "a04e921b",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "base depth: 17, base count: 99\n",
            "ROOT: 50\n"
          ]
        }
      ],
      "source": [
        "# Build a GHZ starting at the best root\n",
        "qc, ghz_qubits, _ = parallel_ghz(\n",
        "    ROOT, GHZ_SIZE, backend, bad_edges, SKIP + bad_nodes\n",
        ")\n",
        "base_depth = qc.depth(lambda x: x.operation.num_qubits == 2)\n",
        "base_count = qc.size(lambda x: x.operation.num_qubits == 2)\n",
        "print(f\"base depth: {base_depth}, base count: {base_count}\")\n",
        "print(f\"ROOT: {ROOT}\")\n",
        "if len(ghz_qubits) != GHZ_SIZE:\n",
        "    raise Exception(\"No GHZ found. Relax error thresholds.\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "21bca6bf",
      "metadata": {},
      "source": [
        "We need one final consideration before inserting valid checks. This is related to the concept of \"coverage\", which is a measure of how many of the wires in a quantum circuit a check can cover. With a higher coverage, we can detect errors on a wider part of the circuit. With this measure, we can select among the valid checks with the highest circuit coverage. In other words, we will be using the `weighted_coverage` function to score different checks for the GHZ circuit.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "a31f2dfb",
      "metadata": {},
      "outputs": [],
      "source": [
        "def weighted_coverage(layers, parities, w_idle=0.2, w_gate=0.8):\n",
        "    \"\"\"\n",
        "    Compute weighted fraction (idle + gate) of wires that are\n",
        "    covered by at least one parity to all active wires.\n",
        "    \"\"\"\n",
        "    wires = active_wires(layers)  # defined below\n",
        "    covered_by_any = {n_layer: set() for n_layer in range(len(layers))}\n",
        "    for parity in parities:\n",
        "        trace = z_trace_backward(layers, parity)  # defined below\n",
        "        for n_layer, qs in trace.items():\n",
        "            covered_by_any[n_layer] |= qs\n",
        "    covered_weight = 0\n",
        "    total_weight = 0\n",
        "    for n_layer in range(len(layers)):\n",
        "        idle = wires[n_layer][\"idle\"]\n",
        "        gate = wires[n_layer][\"gate\"]\n",
        "        total_weight += w_idle * len(idle) + w_gate * len(gate)\n",
        "        covered_idle = covered_by_any[n_layer] & idle\n",
        "        covered_gate = covered_by_any[n_layer] & gate\n",
        "        covered_weight += w_idle * len(covered_idle) + w_gate * len(\n",
        "            covered_gate\n",
        "        )\n",
        "    return covered_weight / total_weight if total_weight > 0 else 0\n",
        "\n",
        "\n",
        "def active_wires(layers):\n",
        "    \"\"\"\n",
        "    Returns per-layer dict with two sets:\n",
        "    - 'idle': activated wires that are idle in this layer\n",
        "    - 'gate': activated wires that are control/target of a CNOT at this layer\n",
        "    \"\"\"\n",
        "    first_activation = {}\n",
        "    for n_layer, layer in enumerate(layers):\n",
        "        for c, t in layer:\n",
        "            first_activation.setdefault(c, n_layer)\n",
        "            first_activation.setdefault(t, n_layer)\n",
        "    result = {}\n",
        "    for n_layer in range(len(layers)):\n",
        "        active = {\n",
        "            q\n",
        "            for q, n_layer0 in first_activation.items()\n",
        "            if n_layer >= n_layer0\n",
        "        }\n",
        "        gate = {q for c, t in layers[n_layer] for q in (c, t)}\n",
        "        idle = active - gate\n",
        "        result[n_layer] = {\"idle\": idle, \"gate\": gate}\n",
        "    return result\n",
        "\n",
        "\n",
        "def z_trace_backward(layers, initial_Zs):\n",
        "    \"\"\"\n",
        "    Backward propagate Zs with parity cancellation.\n",
        "    Returns {layer: set of qubits with odd parity Z at that layer}.\n",
        "    \"\"\"\n",
        "    wires = active_wires(layers)\n",
        "    support = set(initial_Zs)\n",
        "    trace = {}\n",
        "    for n_layer in range(len(layers) - 1, -1, -1):\n",
        "        active = wires[n_layer][\"idle\"] | wires[n_layer][\"gate\"]\n",
        "        trace[n_layer] = support & active\n",
        "        # propagate backwards\n",
        "        new_support = set()\n",
        "        for q in support:\n",
        "            hit = False\n",
        "            for c, t in layers[n_layer]:\n",
        "                if q == t:  # Z on target: copy to control\n",
        "                    new_support ^= {t, c}  # toggle both\n",
        "                    hit = True\n",
        "                    break\n",
        "                elif q == c:  # Z on control: passes through\n",
        "                    new_support ^= {c}\n",
        "                    hit = True\n",
        "                    break\n",
        "            if not hit:  # unaffected\n",
        "                new_support ^= {q}\n",
        "        support = new_support\n",
        "    return trace"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "631f4797",
      "metadata": {},
      "source": [
        "We can now insert checks to the GHZ circuit. Finding valid checks is very convenient for the GHZ state, since any two-qubit Pauli $Z$ operator $Z_i Z_j$ acting on any two qubits $i,j$ of the GHZ circuit is a support and therefore a valid check.\n",
        "\n",
        "Also note that the checks in this case are controlled-$Z$ operators neighboring Hadamard gates from the left and right on the ancilla qubit. This is equivalent to a CNOT gate applied to the ancilla qubit. The code below inserts the checks into the circuit.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "b7e798f7",
      "metadata": {},
      "outputs": [],
      "source": [
        "# --- Tunables controlling the search space / scoring ---\n",
        "MAX_SKIPS = 10  # at most how many qubits to skip (in addition to the bad ones and the ones forced to skip above)\n",
        "SHUFFLES = 200  # how many times to try removing nodes for checks\n",
        "MAX_DEPTH_INCREASE = 10  # how far from the base GHZ depth to go to include checks (increase this for more checks at expense of depth)\n",
        "\n",
        "W_IDLE = 0.2  # weight of errors to consider during idle timesteps\n",
        "W_GATE = 0.8  # weight of errors to consider during gates\n",
        "\n",
        "# Remove random nodes from the GHZ and build from the root again to increase checks\n",
        "degree_two_nodes = [\n",
        "    i\n",
        "    for i in ghz_qubits\n",
        "    if all(n in ghz_qubits for n in coupling_map.neighbors(i))\n",
        "    and len(coupling_map.neighbors(i)) >= 2\n",
        "]\n",
        "\n",
        "# --- Best-so-far tracking for the randomized search ---\n",
        "num_checks = 0\n",
        "best_covered_fraction = -1\n",
        "best_qc = qc\n",
        "best_checks = []\n",
        "best_parities = []\n",
        "best_layers = []\n",
        "\n",
        "# Outer loop: vary how many GHZ nodes we try skipping (0..MAX_SKIPS-1)\n",
        "for num_skips in range(MAX_SKIPS):\n",
        "    # Inner loop: try SHUFFLES random choices of 'num_skips' nodes to skip\n",
        "    for _ in range(SHUFFLES):\n",
        "        # Construct the skip set:\n",
        "        #   - pre-existing forced SKIP\n",
        "        #   - plus a random sample of 'degree_two_nodes' of size 'num_skips'\n",
        "        skip = SKIP + list(np.random.choice(degree_two_nodes, num_skips))\n",
        "\n",
        "        # Rebuild the GHZ using the current skip set and bad_nodes\n",
        "        qc, ghz_qubits, layers = parallel_ghz(\n",
        "            ROOT, GHZ_SIZE, backend, bad_edges, skip + bad_nodes\n",
        "        )\n",
        "\n",
        "        # Measure circuit cost as 2-qubit-gate depth only\n",
        "        depth = qc.depth(lambda x: x.operation.num_qubits == 2)\n",
        "\n",
        "        # If we failed to reach the target GHZ size, discard this attempt\n",
        "        if len(ghz_qubits) != GHZ_SIZE:\n",
        "            continue\n",
        "\n",
        "        # --- Build \"checks\" around the GHZ we just constructed ---\n",
        "        # A check qubit is a non-GHZ, non-dead qubit that has ≥2 neighbors inside the GHZ\n",
        "        # and all those incident edges are usable (i.e., not in bad_edges).\n",
        "        checks = []\n",
        "        parities = []\n",
        "        for i in range(backend.num_qubits):\n",
        "            neighbors = [\n",
        "                n for n in coupling_map.neighbors(i) if n in ghz_qubits\n",
        "            ]\n",
        "\n",
        "            if (\n",
        "                i not in ghz_qubits\n",
        "                and i not in dead_qubits\n",
        "                and len(neighbors) >= 2\n",
        "                and not any(\n",
        "                    [\n",
        "                        [neighbor, i] in bad_edges\n",
        "                        or [i, neighbor] in bad_edges\n",
        "                        for neighbor in neighbors\n",
        "                    ]\n",
        "                )\n",
        "            ):\n",
        "                # Record this qubit as a check qubit\n",
        "                checks.append(i)\n",
        "                parities.append((neighbors[0], neighbors[1]))\n",
        "                # Physically couple the check qubit 'i' to the two GHZ neighbors via CNOTs\n",
        "                # (This is the actual \"check\" attachment in the circuit.)\n",
        "                qc.cx(neighbors[0], i)\n",
        "                qc.cx(neighbors[1], i)\n",
        "\n",
        "        # Score this design using the weighted coverage metric over the GHZ build layers\n",
        "        covered_fraction = weighted_coverage(\n",
        "            layers=layers, parities=parities, w_idle=W_IDLE, w_gate=W_GATE\n",
        "        )\n",
        "\n",
        "        # Keep it only if:\n",
        "        #   - coverage improves over the best so far, AND\n",
        "        #   - the 2q depth budget isn't blown by more than MAX_DEPTH_INCREASE\n",
        "        if (\n",
        "            covered_fraction > best_covered_fraction\n",
        "            and depth <= base_depth + MAX_DEPTH_INCREASE\n",
        "        ):\n",
        "            best_covered_fraction = covered_fraction\n",
        "            best_qc = qc\n",
        "            best_ghz_qubits = ghz_qubits\n",
        "            best_checks = checks\n",
        "            best_parities = parities\n",
        "            best_layers = layers"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8da8af69",
      "metadata": {},
      "source": [
        "We can now print out the qubits used in the GHZ circuit and the check qubits.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "id": "e11e2392",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "GHZ qubits: [50, 49, 51, 38, 52, 48, 58, 53, 47, 71, 39, 46, 70, 54, 33, 45, 72, 69, 55, 32, 37, 73, 68, 34, 31, 44, 25, 74, 78, 67, 18, 24, 79, 75, 89, 57, 11, 23, 93, 59, 88, 66, 10, 22, 92, 90, 87, 65, 12, 9, 21, 94, 91, 86, 77, 13, 8, 20, 95, 98, 97, 14, 7, 36, 99, 111, 107, 15, 6, 41, 115, 110, 106, 19, 17, 5, 40, 114, 109, 108, 105, 27, 4, 42, 118, 104, 28, 3, 129, 117, 103, 29, 2, 128, 125, 96, 30, 127, 124, 102] 100\n",
            "Check qubits: [16, 26, 35, 43, 85, 126] 6\n",
            "Covered fraction (no idle):  0.4595959595959596\n"
          ]
        }
      ],
      "source": [
        "# --- After search, report the best design found ---\n",
        "qc = best_qc\n",
        "checks = best_checks\n",
        "parities = best_parities\n",
        "layers = best_layers\n",
        "ghz_qubits = best_ghz_qubits\n",
        "if len(ghz_qubits) != GHZ_SIZE:\n",
        "    raise Exception(\"No GHZ found. Relax error thresholds.\")\n",
        "\n",
        "print(f\"GHZ qubits: {ghz_qubits} {len(ghz_qubits)}\")\n",
        "print(f\"Check qubits: {checks} {len(checks)}\")\n",
        "\n",
        "covered_fraction = weighted_coverage(\n",
        "    layers=layers, parities=parities, w_idle=W_IDLE, w_gate=W_GATE\n",
        ")\n",
        "print(\n",
        "    \"Covered fraction (no idle): \",\n",
        "    weighted_coverage(\n",
        "        layers=layers, parities=parities, w_idle=0.0, w_gate=1.0\n",
        "    ),\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4b0e4e74",
      "metadata": {},
      "source": [
        "We can also print out some error statistics.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "id": "040d4a25",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "cz errors: \n",
            " mean: 0.002, max: 0.012\n",
            "meas errors: \n",
            " mean: 0.014, max: 0.121\n",
            "t1 errors: \n",
            " mean: 267.9, min: 23.6\n",
            "t2 errors: \n",
            " mean: 155.9, min: 13.9\n"
          ]
        }
      ],
      "source": [
        "def circuit_errors(target, circ, error_type=\"cz\"):\n",
        "    \"\"\"\n",
        "    Pull per-resource error numbers from a Qiskit Target\n",
        "    for ONLY the qubits/edges actually used by `circ`.\n",
        "\n",
        "    Args:\n",
        "        target: qiskit.transpiler.Target (e.g., backend.target)\n",
        "        circ:   qiskit.QuantumCircuit\n",
        "        error_type: one of {\"cz\", \"meas\", \"t1\", \"t2\"}:\n",
        "            - \"cz\"   -> 2q CZ gate error on the circuit's used edges\n",
        "            - \"meas\" -> measurement error on the circuit's used qubits\n",
        "            - \"t1\"   -> T1 (converted to microseconds) on used qubits\n",
        "            - \"t2\"   -> T2 (converted to microseconds) on used qubits\n",
        "\n",
        "    Returns:\n",
        "        list[float] of the requested quantity for the active edges/qubits.\n",
        "    \"\"\"\n",
        "\n",
        "    # Get all 2-qubit edges that appear in the circuit (as undirected pairs).\n",
        "    active_edges = active_gates(circ)  # e.g., {(0,1), (2,3), ...}\n",
        "\n",
        "    # Intersect those with the device coupling map (so we only query valid edges).\n",
        "    # Note: target.build_coupling_map().get_edges() yields directed pairs.\n",
        "    edges = [\n",
        "        edge\n",
        "        for edge in target.build_coupling_map().get_edges()\n",
        "        if tuple(sorted(edge)) in active_edges\n",
        "    ]\n",
        "\n",
        "    # Deduplicate direction: keep only one orientation of each edge.\n",
        "    undirected_edges = []\n",
        "    for edge in edges:\n",
        "        if (edge[1], edge[0]) not in undirected_edges:\n",
        "            undirected_edges.append(edge)\n",
        "    edges = undirected_edges  # (not used later—see note below)\n",
        "\n",
        "    # Accumulators for different error/physics quantities\n",
        "    cz_errors, meas_errors, t1_errors, t2_errors = [], [], [], []\n",
        "\n",
        "    # For every active (undirected) edge in the circuit, fetch its CZ error.\n",
        "    # NOTE: Uses active_gates(circ) again (undirected tuples). This assumes\n",
        "    # `target['cz']` accepts undirected indexing; many Targets store both directions.\n",
        "    for edge in active_gates(circ):\n",
        "        cz_errors.append(target[\"cz\"][edge].error)\n",
        "\n",
        "    # For every active qubit, fetch measure error and T1/T2 (converted to µs).\n",
        "    for qubit in active_qubits(circ):\n",
        "        meas_errors.append(target[\"measure\"][(qubit,)].error)\n",
        "        t1_errors.append(\n",
        "            target.qubit_properties[qubit].t1 * 1e6\n",
        "        )  # seconds -> microseconds\n",
        "        t2_errors.append(\n",
        "            target.qubit_properties[qubit].t2 * 1e6\n",
        "        )  # seconds -> microseconds\n",
        "\n",
        "    # Select which set to return.\n",
        "    if error_type == \"cz\":\n",
        "        return cz_errors\n",
        "    elif error_type == \"meas\":\n",
        "        return meas_errors\n",
        "    elif error_type == \"t1\":\n",
        "        return t1_errors\n",
        "    else:\n",
        "        return t2_errors\n",
        "\n",
        "\n",
        "def active_qubits(circ):\n",
        "    \"\"\"\n",
        "    Return a list of qubit indices that participate in at least one\n",
        "    non-delay, non-barrier instruction in `circ`.\n",
        "    \"\"\"\n",
        "    active_qubits = set()\n",
        "    for inst in circ.data:\n",
        "        # Skip scheduling artifacts that don’t act on state\n",
        "        if (\n",
        "            inst.operation.name != \"delay\"\n",
        "            and inst.operation.name != \"barrier\"\n",
        "        ):\n",
        "            for qubit in inst.qubits:\n",
        "                q = circ.find_bit(\n",
        "                    qubit\n",
        "                ).index  # map Qubit object -> integer index\n",
        "                active_qubits.add(q)\n",
        "    return list(active_qubits)\n",
        "\n",
        "\n",
        "def active_gates(circ):\n",
        "    \"\"\"\n",
        "    Return a set of undirected 2-qubit edges (i, j) that appear in `circ`.\n",
        "    \"\"\"\n",
        "    used_2q_gates = set()\n",
        "    for inst in circ:\n",
        "        if inst.operation.num_qubits == 2:\n",
        "            qs = inst.qubits\n",
        "            # map Qubit objects -> indices, then sort to make the edge undirected\n",
        "            qs = sorted([circ.find_bit(q).index for q in qs])\n",
        "            used_2q_gates.add(tuple(sorted(qs)))\n",
        "    return used_2q_gates\n",
        "\n",
        "\n",
        "# ---- Print summary statistics ----\n",
        "cz_errors = circuit_errors(backend.target, qc, error_type=\"cz\")\n",
        "meas_errors = circuit_errors(backend.target, qc, error_type=\"meas\")\n",
        "t1_errors = circuit_errors(backend.target, qc, error_type=\"t1\")\n",
        "t2_errors = circuit_errors(backend.target, qc, error_type=\"t2\")\n",
        "\n",
        "np.set_printoptions(linewidth=np.inf)\n",
        "print(\n",
        "    f\"cz errors: \\n mean: {np.round(np.mean(cz_errors), 3)}, max: {np.round(np.max(cz_errors), 3)}\"\n",
        ")\n",
        "print(\n",
        "    f\"meas errors: \\n mean: {np.round(np.mean(meas_errors), 3)}, max: {np.round(np.max(meas_errors), 3)}\"\n",
        ")\n",
        "print(\n",
        "    f\"t1 errors: \\n mean: {np.round(np.mean(t1_errors), 1)}, min: {np.round(np.min(t1_errors), 1)}\"\n",
        ")\n",
        "print(\n",
        "    f\"t2 errors: \\n mean: {np.round(np.mean(t2_errors), 1)}, min: {np.round(np.min(t2_errors), 1)}\"\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a8f0a54f",
      "metadata": {},
      "source": [
        "As previously, we can simulate the circuit first in the absence of noise to ensure correctness of the GHZ state preparation circuit.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "b905e875",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Stabilizer simulation result:\n",
            "{'000000 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111': 525, '000000 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000': 475}\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/error-detection/extracted-outputs/b905e875-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 23,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# --- Simulate to ensure correctness ---\n",
        "\n",
        "qc_meas = qc.copy()\n",
        "\n",
        "# Add measurements to the GHZ qubits\n",
        "c1 = ClassicalRegister(len(ghz_qubits), \"c1\")\n",
        "qc_meas.add_register(c1)\n",
        "for q, c in zip(ghz_qubits, c1):\n",
        "    qc_meas.measure(q, c)\n",
        "\n",
        "# Add measurements to the check qubits\n",
        "if len(checks) > 0:\n",
        "    c2 = ClassicalRegister(len(checks), \"c2\")\n",
        "    qc_meas.add_register(c2)\n",
        "    for q, c in zip(checks, c2):\n",
        "        qc_meas.measure(q, c)\n",
        "\n",
        "# Simulate the circuit with stabilizer method\n",
        "sim_stab = AerSimulator(method=\"stabilizer\")\n",
        "res = sim_stab.run(qc_meas, shots=1000).result()\n",
        "counts = res.get_counts()\n",
        "print(\"Stabilizer simulation result:\")\n",
        "print(counts)\n",
        "\n",
        "# Rename keys to \"0 0\" and \"0 1\" for easier plotting\n",
        "# First len(checks) bits are check bits, rest are GHZ bits\n",
        "keys = list(counts.keys())\n",
        "for key in keys:\n",
        "    check_bits = key[: len(checks)]\n",
        "    ghz_bits = key[(len(checks) + 1) :]\n",
        "    if set(check_bits) == {\"0\"} and set(ghz_bits) == {\"0\"}:\n",
        "        counts[\"0 0\"] = counts.pop(key)\n",
        "    elif set(check_bits) == {\"0\"} and set(ghz_bits) == {\"1\"}:\n",
        "        counts[\"0 1\"] = counts.pop(key)\n",
        "    else:\n",
        "        continue\n",
        "\n",
        "plot_histogram(counts)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a5db4e3b",
      "metadata": {},
      "source": [
        "As expected, the check qubits are measured as all zeros, and we successfully prepared the GHZ state.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "14770785",
      "metadata": {},
      "source": [
        "### Step 3: Execute using Qiskit primitives\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d7a6e8a8",
      "metadata": {},
      "source": [
        "Now we are ready to run the circuit on real hardware and demonstrate how the error detection protocol can capture errors in GHZ state preparation.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "e2d37555",
      "metadata": {},
      "outputs": [],
      "source": [
        "BAD_QUBITS = []  # specify any additional bad qubits to avoid (this is specific to the chosen backend)\n",
        "SHOTS = 10000  # number of shots"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "eef98bd0",
      "metadata": {},
      "source": [
        "We define a helper function to add measurements to the GHZ circuit.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "5f99ba1f",
      "metadata": {},
      "outputs": [],
      "source": [
        "def add_measurements(qc, ghz_qubits, checks):\n",
        "    # --- Measure each set of qubits into different classical registers to facilitate post-processing ---\n",
        "\n",
        "    # Add measurements to the GHZ qubits\n",
        "    c1 = ClassicalRegister(len(ghz_qubits), \"c1\")\n",
        "    qc.add_register(c1)\n",
        "    for q, c in zip(ghz_qubits, c1):\n",
        "        qc.measure(q, c)\n",
        "\n",
        "    # Add measurements to the check qubits\n",
        "    c2 = ClassicalRegister(len(checks), \"c2\")\n",
        "    qc.add_register(c2)\n",
        "    for q, c in zip(checks, c2):\n",
        "        qc.measure(q, c)\n",
        "\n",
        "    return qc"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "19eae8b4",
      "metadata": {},
      "source": [
        "Before execution, we draw the layout of the GHZ qubits and the check qubits on the selected hardware.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "d0faf114",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/error-detection/extracted-outputs/d0faf114-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Plot the layout of GHZ and check qubits on the device\n",
        "plot_gate_map(\n",
        "    backend,\n",
        "    label_qubits=True,\n",
        "    line_width=20,\n",
        "    line_color=[\n",
        "        \"black\"\n",
        "        if edge[0] in ghz_qubits + checks and edge[1] in ghz_qubits + checks\n",
        "        else \"lightgrey\"\n",
        "        for edge in backend.coupling_map.graph.edge_list()\n",
        "    ],\n",
        "    qubit_color=[\n",
        "        \"blue\"\n",
        "        if i in ghz_qubits\n",
        "        else \"salmon\"\n",
        "        if i in checks\n",
        "        else \"lightgrey\"\n",
        "        for i in range(0, backend.num_qubits)\n",
        "    ],\n",
        ")\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 29,
      "id": "78c8c2b6",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/error-detection/extracted-outputs/78c8c2b6-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 27,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "qc.draw(\"mpl\", idle_wires=False, fold=-1)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "681d4f29",
      "metadata": {},
      "source": [
        "We now add the measurements.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 30,
      "id": "75169ff2",
      "metadata": {},
      "outputs": [],
      "source": [
        "qc = add_measurements(qc, ghz_qubits, checks)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9bf787f5",
      "metadata": {},
      "source": [
        "The scheduling pipeline below locks in timing, removes barriers, simplifies delays, and injects dynamical decoupling, all while preserving the original operation times.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 31,
      "id": "673d1a87",
      "metadata": {},
      "outputs": [],
      "source": [
        "# The scheduling consists of first inserting delays while barriers are still there\n",
        "# Then removing the barriers and consolidating the delays, so that the operations do not move in time\n",
        "# Lastly we replace delays with dynamical decoupling\n",
        "collect_function = partial(\n",
        "    collect_using_filter_function,\n",
        "    filter_function=(lambda node: node.op.name == \"delay\"),\n",
        "    split_blocks=True,\n",
        "    min_block_size=2,\n",
        "    split_layers=False,\n",
        "    collect_from_back=False,\n",
        "    max_block_width=None,\n",
        ")\n",
        "\n",
        "collapse_function = partial(\n",
        "    collapse_to_operation,\n",
        "    collapse_function=(\n",
        "        lambda circ: Delay(sum(inst.operation.duration for inst in circ))\n",
        "    ),\n",
        ")\n",
        "\n",
        "\n",
        "class Unschedule(AnalysisPass):\n",
        "    \"\"\"Removes a property from the passmanager property set so that the circuit looks unscheduled, so we can schedule it again.\"\"\"\n",
        "\n",
        "    def run(self, dag):\n",
        "        del self.property_set[\"node_start_time\"]\n",
        "\n",
        "\n",
        "def build_passmanager(backend, dd_qubits=None):\n",
        "    pm = generate_preset_pass_manager(\n",
        "        target=backend.target,\n",
        "        layout_method=\"trivial\",\n",
        "        optimization_level=2,\n",
        "        routing_method=\"none\",\n",
        "    )\n",
        "\n",
        "    pm.scheduling = PassManager(\n",
        "        [\n",
        "            ALAPScheduleAnalysis(target=backend.target),\n",
        "            PadDelay(target=backend.target),\n",
        "            RemoveBarriers(),\n",
        "            Unschedule(),\n",
        "            CollectAndCollapse(\n",
        "                collect_function=collect_function,\n",
        "                collapse_function=collapse_function,\n",
        "            ),\n",
        "            ALAPScheduleAnalysis(target=backend.target),\n",
        "            PadDynamicalDecoupling(\n",
        "                dd_sequence=[XGate(), RZGate(-np.pi), XGate(), RZGate(np.pi)],\n",
        "                spacing=[1 / 4, 1 / 2, 0, 0, 1 / 4],\n",
        "                target=backend.target,\n",
        "                qubits=dd_qubits,\n",
        "            ),\n",
        "        ]\n",
        "    )\n",
        "\n",
        "    return pm"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6b6c4078",
      "metadata": {},
      "source": [
        "We can now use the custom pass manager to transpile the circuit for the selected backend.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "8b8c537b",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Transpile the circuits for the backend\n",
        "pm = build_passmanager(backend, ghz_qubits)\n",
        "\n",
        "# Instruction set architecture (ISA) level circuit after scheduling and DD insertion\n",
        "isa_circuit = pm.run(qc)\n",
        "\n",
        "# Draw after scheduling and DD insertion\n",
        "# timeline_drawer(isa_circuit, show_idle=False, time_range=(0, 1000), target=backend.target)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 33,
      "id": "0c4b279b",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/error-detection/extracted-outputs/0c4b279b-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 31,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "isa_circuit.draw(\"mpl\", fold=-1, idle_wires=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4e0712a5",
      "metadata": {},
      "source": [
        "We then submit the job using the Qiskit Runtime Sampler primitive.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "803e5790",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "d493f17nmdfs73abf9qg\n"
          ]
        }
      ],
      "source": [
        "# Select the sampler options\n",
        "sampler = Sampler(mode=backend)\n",
        "sampler.options.default_shots = SHOTS\n",
        "sampler.options.dynamical_decoupling.enable = False\n",
        "sampler.options.execution.rep_delay = 0.00025\n",
        "\n",
        "# Submit the job\n",
        "print(\"Submitting sampler job\")\n",
        "ghz_job = sampler.run([isa_circuit])\n",
        "\n",
        "print(ghz_job.job_id())"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c39bfe38",
      "metadata": {},
      "source": [
        "### Step 4: Post-process and return result in desired classical format\n",
        "\n",
        "We can now retrieve and analyze the results from the Sampler job.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 33,
      "id": "23495738",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Retrieve the job results\n",
        "job_result = ghz_job.result()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 34,
      "id": "55e07409",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Get the counts from GHZ and check qubit measurements\n",
        "ghz_counts = job_result[0].data.c1.get_counts()\n",
        "checks_counts = job_result[0].data.c2.get_counts()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 35,
      "id": "d7162ed7",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Post-process to get unflagged GHZ counts (i.e., check bits are all '0')\n",
        "joined_counts = job_result[0].join_data().get_counts()\n",
        "unflagged_counts = {}\n",
        "for key, count in joined_counts.items():\n",
        "    check_bits = key[: len(checks)]\n",
        "    ghz_bits = key[len(checks) :]\n",
        "    if set(check_bits) == {\"0\"}:\n",
        "        unflagged_counts[ghz_bits] = count"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 36,
      "id": "d7902975",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/error-detection/extracted-outputs/d7902975-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 36,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Get top 20 outcomes by frequency from the unflagged counts\n",
        "top_counts = dict(\n",
        "    sorted(unflagged_counts.items(), key=lambda x: x[1], reverse=True)[:20]\n",
        ")\n",
        "\n",
        "# Rename keys for better visualization\n",
        "top_counts_renamed = {}\n",
        "i = 0\n",
        "for key, count in top_counts.items():\n",
        "    if set(key) == {\"0\"}:\n",
        "        top_counts_renamed[\"all 0s\"] = count\n",
        "    elif set(key) == {\"1\"}:\n",
        "        top_counts_renamed[\"all 1s\"] = count\n",
        "    else:\n",
        "        top_counts_renamed[f\"other_{i}\"] = count\n",
        "        i += 1\n",
        "\n",
        "plot_histogram(top_counts_renamed, figsize=(12, 7))"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "90505f1c",
      "metadata": {},
      "source": [
        "In the histogram above, we plotted 20 bitstring measurements from the GHZ qubits that were not flagged by the check qubits. As expected, all-0 and all-1 bitstrings have the highest counts. Note that some erroneous bitstrings with low error weights were not captured by the error detection. The highest counts are still found in the expected bitstrings.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5423c559",
      "metadata": {},
      "source": [
        "## Discussion\n",
        "\n",
        "In this tutorial, we showed how to implement a low-overhead error detection technique using spacetime codes, and demonstrated its real-world application to preparing GHZ states on hardware. Refer to [\\[3\\]](https://arxiv.org/abs/2510.09520) to further explore technical details of GHZ state preparation. In addition to error detection, authors utilize readout error mitigation with M3 and TREX and perform temporary uncomputation techniques to prepare high-fidelity GHZ states.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "7bbdd161",
      "metadata": {},
      "source": [
        "### References\n",
        "\n",
        "* \\[1] Martiel, S., & Javadi-Abhari, A. (2025). Low-overhead error detection with spacetime codes. *arXiv preprint arXiv: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",
        "* \\[3] Javadi-Abhari, A., Martiel, S., Seif, A., Takita, M., & Wei, K. X. (2025). Big cats: entanglement in 120 qubits and beyond. arXiv preprint arXiv:2510.09520.\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"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}