{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "446c0650-d248-4ab1-8dc0-41921ea624da",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Cancelación probabilística de errores con modelos de ruido lógico\"\n",
        "description: \"Combinar la detección de errores mediante qubits mediadores con la PEC en un circuito de Ising hexagonal de 49 qubits para mitigar el ruido a una fracción del coste de muestreo que supone la PEC por sí sola.\"\n",
        "---\n",
        "\n",
        "{/* cspell:ignore xslow edgecolor zorder fontsize facecolor markerfacecolor ncols unconverged frameon DSATUR NNLS fracs unmodeled */}\n",
        "\n",
        "<span id=\"probabilistic-error-cancellation-with-logical-noise-models\" />\n",
        "\n",
        "# Cancelación probabilística de errores con modelos de ruido lógico\n",
        "\n",
        "Estimación *de tiempo de ejecución: 28 minutos en un procesador Heron r3 (NOTA: Se trata únicamente de una estimación. (El tiempo de ejecución puede variar.)*\n",
        "\n",
        "![Red de Ising hexagonal integrada en una disposición de qubits «heavy-hex», con sitios de Ising en los qubits « degree-3 » y qubits mediadores en los bordes que los unen](https://quantum.cloud.ibm.com/docs/images/tutorials/probabilistic-error-cancellation-with-logical-noise-models/hex-ising.avif)\n",
        "En este tutorial evaluaremos los observables de un modelo de Ising de 22 sitios en una red ***hexagonal*** utilizando una QPU Heron, que cuenta con una topología de qubits ***«heavy-hexagonal*** ». Muchas demostraciones realizadas con las QPU Heron se centran en sistemas definidos en redes hexagonales densas para adaptarse a la conectividad del sistema; sin embargo, los modelos hexagonales suelen resultar más interesantes de estudiar, ya que son más habituales en la naturaleza y más difíciles de simular mediante métodos clásicos debido a su mayor densidad de conectividad. Dado que la topología de qubits de la QPU no admite directamente la conectividad de una red hexagonal, debemos decidir cómo integrar de manera eficiente el modelo hexagonal en la topología de qubits. En este ejemplo, representamos cada sitio del modelo de Ising con un qubit situado en un vértice de la red «heavy-hex» de la QPU. Utilizamos los qubits de los bordes ( ***qubits mediadores*** ) para facilitar el entrelazamiento entre los qubits de Ising de los vértices. Además, los qubits mediadores implementan el entrelazamiento de tal forma que siempre se espera que vuelvan al estado fundamental $|0\\rangle$. Si un qubit mediador da como resultado $|1\\rangle$, esto indica que se ha producido un error lógico durante la ejecución del circuito; al realizar una postselección que incluya únicamente las muestras en las que no se hayan detectado errores en los qubits mediadores, se obtiene una distribución más reducida de ***muestras lógicas*** que podría tener una mayor fidelidad que la distribución bruta, afectada por el ruido. No solo podemos detectar que se ha producido un error midiendo un qubit mediador, sino que también podemos determinar exactamente qué errores lógicos es capaz de detectar ese qubit a lo largo de todo el circuito. Al eliminar de un modelo de ruido los generadores de ruido que ***las comprobaciones de simetría*** de los qubits mediadores son capaces de detectar, se obtiene un modelo de ruido reducido, que puede mitigarse mediante técnicas como la cancelación probabilística de errores (PEC).\n",
        "\n",
        "En este ejemplo de cuaderno, combinaremos la técnica de detección de errores descrita anteriormente con el PEC para contrarrestar el ruido cuántico de forma más eficaz de lo que lo haría cualquiera de las dos técnicas por separado. Utilizaremos 49 qubits para integrar el modelo de Ising de 22 qubits y emplearemos los 27 qubits mediadores adicionales para la detección de errores. Aplicaremos el método PEC al canal de ruido preseleccionado para mitigar el ruido que escapa a los controles de simetría. Además de combinar la detección de errores con el PEC, utilizaremos técnicas como la mitigación de errores de lectura TREX y las comprobaciones de errores no markovianas para contrarrestar aún más el efecto del ruido cuántico.\n",
        "\n",
        "El flujo de trabajo es el siguiente:\n",
        "\n",
        "1. Simulación clásica de los valores esperados de las observables para el modelo hex-Ising de 22 qubits\n",
        "2. Implementar el modelo hex-Ising de detección de errores de 49 qubits mediante un circuito cuántico y transpilarlo al backend de la QPU\n",
        "3. Especifica las capas de entrelazamiento del circuito con `samplomatic`. Analizaremos y mitigaremos el ruido que afecta a estas capas entrelazadas.\n",
        "4. Añadir comprobaciones de errores no markovianos al circuito\n",
        "5. Conoce el ruido de puerta y de lectura que afecta al circuito\n",
        "6. Recortar el modelo de ruido de los términos detectados por las comprobaciones de simetría\n",
        "7. Prueba el circuito de detección de errores. Selecciona a posteriori únicamente las muestras que superen todas las comprobaciones de simetría y de errores no markovianos, y utiliza la mitigación de lectura TREX para todos los cálculos de valores esperados\n",
        "   * Prueba el circuito de detección de errores\n",
        "   * Realiza un muestreo del circuito de detección de errores con PEC, pero no apliques la postselección basada en comprobaciones de simetría y mitiga todo el canal de ruido aprendido.\n",
        "   * Prueba el circuito de detección de errores con PEC. Realiza la postselección basándote en comprobaciones de simetría y mitiga únicamente el ruido que no pueda detectarse mediante dichas comprobaciones.\n",
        "8. Calcula los valores esperados y compara estrategias. Obsérvese que el método QED+PEC elimina con mayor eficacia el ruido que afecta al experimento que cualquiera de los dos métodos por separado y produce valores esperados convergentes utilizando un número de disparos mucho menor que el que requiere el método PEC por sí solo.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1d09a34f",
      "metadata": {},
      "source": [
        "<span id=\"requirements\" />\n",
        "\n",
        "## Requisitos\n",
        "\n",
        "Antes de empezar este tutorial, ejecuta el siguiente comando para instalar todos los paquetes necesarios:\n",
        "\n",
        "`%pip install networkx numpy \"qiskit[visualization]\" \"qiskit-ibm-runtime[visualization]\" samplomatic qiskit-mitigation matplotlib \"qiskit-noise-learning @ git+https://github.com/Qiskit/qiskit-noise-learning.git@ieee-demo-2026\"`\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "60f6ee26",
      "metadata": {},
      "source": [
        "La celda colapsada que aparece a continuación define los elementos auxiliares de las figuras que se utilizan a lo largo de este cuaderno.\n",
        "\n",
        "<Accordion>\n",
        "  <AccordionItem title=\"Haz clic para ver las ayudas de las figuras\">\n",
        "    <CodeCellPlaceholder tag=\"id-figure-helpers\" />\n",
        "  </AccordionItem>\n",
        "</Accordion>\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "4055a12e",
      "metadata": {
        "execution": {
          "iopub.execute_input": "2026-09-04T05:06:09.329390Z",
          "iopub.status.busy": "2026-09-04T05:06:09.329197Z",
          "iopub.status.idle": "2026-09-04T05:06:09.589989Z",
          "shell.execute_reply": "2026-09-04T05:06:09.589585Z"
        },
        "jupyter": {
          "source_hidden": true
        },
        "tags": [
          "id-figure-helpers"
        ]
      },
      "outputs": [],
      "source": [
        "# Figure helpers for this notebook (collapsed): all plotting and styling lives here.\n",
        "# The chip maps follow the style of Fig. 30 of arXiv:2607.25998.\n",
        "\n",
        "from math import comb\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "from matplotlib import cm\n",
        "from matplotlib.colors import BoundaryNorm\n",
        "from matplotlib.patches import Circle, Rectangle, Wedge\n",
        "\n",
        "# one typography scheme for every figure in the notebook\n",
        "plt.rcParams.update(\n",
        "    {\n",
        "        \"font.size\": 13,\n",
        "        \"axes.titlesize\": 15,\n",
        "        \"axes.labelsize\": 13,\n",
        "        \"xtick.labelsize\": 11,\n",
        "        \"ytick.labelsize\": 11,\n",
        "        \"legend.fontsize\": 11,\n",
        "    }\n",
        ")\n",
        "\n",
        "\n",
        "def x_labels(layout, n_data):\n",
        "    \"\"\"Per-observable tick labels carrying the physical qubit indices (site i = qubit i).\"\"\"\n",
        "    return [f\"$X_{{{layout[i]}}}$\" for i in range(n_data)]\n",
        "\n",
        "\n",
        "def plot_layout(backend, layout, n_data):\n",
        "    \"\"\"Chip cartoon of the embedding: data qubits green, check qubits orange.\"\"\"\n",
        "    from qiskit.visualization import plot_coupling_map\n",
        "\n",
        "    return plot_coupling_map(\n",
        "        num_qubits=backend.num_qubits,\n",
        "        qubit_coordinates=None,\n",
        "        coupling_map=list(backend.coupling_map.get_edges()),\n",
        "        figsize=(9, 9),\n",
        "        qubit_color=[\n",
        "            \"#4CAF50\"\n",
        "            if q in layout[:n_data]\n",
        "            else \"#FF9800\"\n",
        "            if q in layout[n_data:]\n",
        "            else \"#DDDDDD\"\n",
        "            for q in range(backend.num_qubits)\n",
        "        ],\n",
        "        qubit_size=220,\n",
        "        line_width=2,\n",
        "        font_size=90,\n",
        "    )\n",
        "\n",
        "\n",
        "def plot_exact(obs_exact, tick_labels, title):\n",
        "    plt.figure(figsize=(12, 4))\n",
        "    plt.plot(obs_exact, \"o-\")\n",
        "    plt.title(title)\n",
        "    plt.xticks(np.arange(len(obs_exact)), tick_labels)\n",
        "    plt.xlabel(\"Observable\")\n",
        "    plt.ylabel(r\"$\\langle X \\rangle$\")\n",
        "    plt.grid()\n",
        "\n",
        "\n",
        "def draw_toy_circuit(\n",
        "    generate_ed_ising, zz_coeff, x_coeff, include_checks=True\n",
        "):\n",
        "    \"\"\"The boxing pipeline on a 3-plaquette-ring miniature (1 Trotter step) so the box\n",
        "    structure is legible; every box carries the Twirl / InjectNoise annotations, and with\n",
        "    ``include_checks`` the terminal xslow non-Markovian error check pattern is appended as in the\n",
        "    production pipeline.\"\"\"\n",
        "    import networkx as nx\n",
        "    from qiskit.circuit import ClassicalRegister\n",
        "    from qiskit_mitigation.postselection import XSlowGate\n",
        "    from samplomatic.transpiler import generate_boxing_pass_manager\n",
        "\n",
        "    toy, _, _ = generate_ed_ising(nx.cycle_graph(3), 1, zz_coeff, x_coeff)\n",
        "    toy.add_register(\n",
        "        ClassicalRegister(3, \"data\"), ClassicalRegister(3, \"check\")\n",
        "    )\n",
        "    toy.barrier()\n",
        "    toy.measure(range(3), range(3))\n",
        "    toy.measure(range(3, 6), range(3, 6))\n",
        "    toy_boxed = generate_boxing_pass_manager(\n",
        "        enable_gates=True,\n",
        "        enable_measures=True,\n",
        "        inject_noise_targets=\"gates\",\n",
        "        inject_noise_strategy=\"individual_modification\",\n",
        "        inject_noise_site=\"after\",\n",
        "        twirling_strategy=\"active_circuit\",\n",
        "        measure_annotations=\"all\",\n",
        "    ).run(toy)\n",
        "    if include_checks:\n",
        "        toy_boxed.add_register(\n",
        "            ClassicalRegister(3, \"data_ps\"), ClassicalRegister(3, \"check_ps\")\n",
        "        )\n",
        "        toy_boxed.barrier()\n",
        "        for qb in range(6):\n",
        "            toy_boxed.append(XSlowGate(), [qb])\n",
        "        toy_boxed.measure(range(3), toy_boxed.cregs[2])\n",
        "        toy_boxed.measure(range(3, 6), toy_boxed.cregs[3])\n",
        "    return toy_boxed.draw(\"mpl\", fold=-1, scale=0.6)\n",
        "\n",
        "\n",
        "# --- chip-level noise maps ---------------------------------------------------------\n",
        "\n",
        "BANDS = [1e-5, 2e-5, 3e-5, 4e-5, 6e-5, 1e-4, 2e-4, 3e-4, 4e-4, 6e-4, 1e-3]\n",
        "CMAP = plt.get_cmap(\"YlOrBr\")\n",
        "NORM = BoundaryNorm(BANDS, CMAP.N, extend=\"both\")\n",
        "\n",
        "\n",
        "def _color(rate):\n",
        "    return (\n",
        "        \"white\"\n",
        "        if rate < BANDS[0]\n",
        "        else CMAP(NORM(min(rate, BANDS[-1] * 0.999)))\n",
        "    )\n",
        "\n",
        "\n",
        "def _layer_sparse(mit, weights):\n",
        "    \"\"\"Per-layer labelled terms from the saved run, box-local -> physical qubits.\"\"\"\n",
        "    phys = np.sort(mit[\"layout\"])\n",
        "    return [\n",
        "        [\n",
        "            (p, tuple(int(phys[q]) for q in qs if q >= 0), r)\n",
        "            for p, qs, r in zip(\n",
        "                mit[f\"label_paulis_{i}\"],\n",
        "                mit[f\"label_qubits_{i}\"],\n",
        "                w,\n",
        "                strict=True,\n",
        "            )\n",
        "        ]\n",
        "        for i, w in enumerate(weights)\n",
        "    ]\n",
        "\n",
        "\n",
        "def _aggregate(layers):\n",
        "    \"\"\"w1[qubit][P] and w2[(a,b)][PaPb]: rates summed over the 3 layers.\"\"\"\n",
        "    w1, w2 = {}, {}\n",
        "    for terms in layers:\n",
        "        for p, qs, r in terms:\n",
        "            if len(qs) == 1:\n",
        "                w1.setdefault(qs[0], dict.fromkeys(\"XYZ\", 0.0))[p] += r\n",
        "            else:\n",
        "                (a, pa), (b, pb) = sorted(zip(qs, p, strict=True))\n",
        "                w2.setdefault(\n",
        "                    (a, b), {x + y: 0.0 for x in \"XYZ\" for y in \"XYZ\"}\n",
        "                )[pa + pb] += r\n",
        "    return w1, w2\n",
        "\n",
        "\n",
        "def draw_noise_map(mit, backend, reduced=False):\n",
        "    \"\"\"Chip map of the learned model: X/Y/Z wheel per qubit, 3x3 two-qubit Pauli grid per\n",
        "    coupler, log-banded colors, dashed outlines for unused hardware.\n",
        "\n",
        "    With ``reduced=True``, keeps only the error terms the checks cannot see: detectability\n",
        "    depends on circuit position, so each layer's 0/1 site scales are averaged over its uses.\n",
        "    \"\"\"\n",
        "    from qiskit_ibm_runtime.visualization.embeddings import Embedding\n",
        "\n",
        "    if reduced:\n",
        "        scales = [\n",
        "            mit[\"site_scales\"][mit[\"site_layer\"] == i].mean(axis=0)\n",
        "            for i in range(3)\n",
        "        ]\n",
        "        weights = [mit[f\"rates_{i}\"] * scales[i] for i in range(3)]\n",
        "        title = f\"Reduced noise model ($\\\\gamma$ = {mit['gammas'][1]:.1f})\"\n",
        "    else:\n",
        "        weights = [mit[f\"rates_{i}\"] for i in range(3)]\n",
        "        title = f\"Full noise model ($\\\\gamma$ = {mit['gammas'][0]:.0f})\"\n",
        "    w1, w2 = _aggregate(_layer_sparse(mit, weights))\n",
        "\n",
        "    xy = np.array(\n",
        "        [(c, -r) for r, c in Embedding.from_backend(backend).coordinates]\n",
        "    )\n",
        "    fig, ax = plt.subplots(figsize=(13, 7))\n",
        "    for a, b in {tuple(sorted(e)) for e in backend.coupling_map.get_edges()}:\n",
        "        if (\n",
        "            (\n",
        "                a,\n",
        "                b,\n",
        "            )\n",
        "            in w2\n",
        "        ):  # 3x3 Pauli grid laid along the bond (columns: qubit a, rows: qubit b)\n",
        "            d = xy[b] - xy[a]\n",
        "            u = d / np.hypot(*d)\n",
        "            v = np.array([-u[1], u[0]])\n",
        "            cl, cw = (np.hypot(*d) - 0.6) / 3, 0.17\n",
        "            for i, Pa in enumerate(\"XYZ\"):\n",
        "                for j, Pb in enumerate(\"XYZ\"):\n",
        "                    ax.add_patch(\n",
        "                        Rectangle(\n",
        "                            xy[a] + u * (0.3 + i * cl) + v * ((j - 1.5) * cw),\n",
        "                            cl,\n",
        "                            cw,\n",
        "                            angle=np.degrees(np.arctan2(u[1], u[0])),\n",
        "                            facecolor=_color(w2[a, b][Pa + Pb]),\n",
        "                            edgecolor=\"black\",\n",
        "                            lw=0.4,\n",
        "                            zorder=2,\n",
        "                        )\n",
        "                    )\n",
        "        else:\n",
        "            ax.plot(\n",
        "                *zip(xy[a], xy[b], strict=True),\n",
        "                ls=\"--\",\n",
        "                lw=0.7,\n",
        "                color=\"black\",\n",
        "                alpha=0.5,\n",
        "                zorder=1,\n",
        "            )\n",
        "    for q, (xq, yq) in enumerate(xy):\n",
        "        if q in w1:  # three-sector wheel: X top, Y lower left, Z lower right\n",
        "            for P, t0 in ((\"X\", 30), (\"Y\", 150), (\"Z\", 270)):\n",
        "                ax.add_patch(\n",
        "                    Wedge(\n",
        "                        (xq, yq),\n",
        "                        0.3,\n",
        "                        t0,\n",
        "                        t0 + 120,\n",
        "                        facecolor=_color(w1[q][P]),\n",
        "                        edgecolor=\"black\",\n",
        "                        lw=0.6,\n",
        "                        zorder=3,\n",
        "                    )\n",
        "                )\n",
        "            ax.annotate(\n",
        "                str(q),\n",
        "                (xq + 0.39, yq - 0.39),\n",
        "                fontsize=9,\n",
        "                color=\"gray\",\n",
        "                ha=\"left\",\n",
        "                va=\"top\",\n",
        "                zorder=4,\n",
        "            )  # southeast, clear of the bond grids\n",
        "        else:\n",
        "            ax.add_patch(\n",
        "                Circle(\n",
        "                    (xq, yq),\n",
        "                    0.24,\n",
        "                    facecolor=\"none\",\n",
        "                    edgecolor=\"black\",\n",
        "                    ls=\"--\",\n",
        "                    lw=0.7,\n",
        "                    alpha=0.5,\n",
        "                    zorder=3,\n",
        "                )\n",
        "            )\n",
        "    ux = xy[sorted(w1)]\n",
        "    ax.set_xlim(ux[:, 0].min() - 2.2, ux[:, 0].max() + 2.2)\n",
        "    ax.set_ylim(ux[:, 1].min() - 1.6, ux[:, 1].max() + 1.6)\n",
        "    ax.set_aspect(\"equal\")\n",
        "    ax.axis(\"off\")\n",
        "    ax.set_title(title, fontsize=15)\n",
        "    cb = fig.colorbar(\n",
        "        cm.ScalarMappable(norm=NORM, cmap=CMAP),\n",
        "        ax=ax,\n",
        "        fraction=0.035,\n",
        "        pad=0.02,\n",
        "        extend=\"both\",\n",
        "        ticks=BANDS,\n",
        "    )\n",
        "    cb.set_label(\"coefficient (log bands; white < 1e-05)\", fontsize=11)\n",
        "    cb.ax.set_yticklabels(\n",
        "        [f\"{b:.0e}\".replace(\"e-0\", \"e-\") for b in BANDS], fontsize=10\n",
        "    )\n",
        "    _noise_map_legends(fig, ax)\n",
        "\n",
        "\n",
        "def _noise_map_legends(fig, ax):\n",
        "    \"\"\"Weight-1 wheel and weight-2 grid keys, in a reserved band left of the lattice.\"\"\"\n",
        "    fig.subplots_adjust(left=0.17)\n",
        "    axl = ax.inset_axes([-0.185, 0.70, 0.13, 0.24])\n",
        "    for P, t0 in ((\"X\", 30), (\"Y\", 150), (\"Z\", 270)):\n",
        "        axl.add_patch(\n",
        "            Wedge(\n",
        "                (0.5, 0.45),\n",
        "                0.38,\n",
        "                t0,\n",
        "                t0 + 120,\n",
        "                facecolor=\"white\",\n",
        "                edgecolor=\"black\",\n",
        "                lw=0.8,\n",
        "            )\n",
        "        )\n",
        "        axl.annotate(\n",
        "            P,\n",
        "            (\n",
        "                0.5 + 0.2 * np.cos(np.radians(t0 + 60)),\n",
        "                0.45 + 0.2 * np.sin(np.radians(t0 + 60)),\n",
        "            ),\n",
        "            ha=\"center\",\n",
        "            va=\"center\",\n",
        "            fontsize=9,\n",
        "        )\n",
        "    axl.set_title(\"weight-1\", fontsize=10)\n",
        "    axl.set_xlim(0, 1)\n",
        "    axl.set_ylim(0, 1)\n",
        "    axl.set_aspect(\"equal\")\n",
        "    axl.axis(\"off\")\n",
        "    axm = ax.inset_axes([-0.185, 0.32, 0.14, 0.30])\n",
        "    for i, Pa in enumerate(\"XYZ\"):\n",
        "        for j, Pb in enumerate(\"XYZ\"):\n",
        "            axm.add_patch(\n",
        "                Rectangle(\n",
        "                    (i / 3, 1 - (j + 1) / 3),\n",
        "                    1 / 3,\n",
        "                    1 / 3,\n",
        "                    facecolor=\"white\",\n",
        "                    edgecolor=\"black\",\n",
        "                    lw=0.6,\n",
        "                )\n",
        "            )\n",
        "            axm.annotate(\n",
        "                Pa + Pb,\n",
        "                ((i + 0.5) / 3, 1 - (j + 0.5) / 3),\n",
        "                ha=\"center\",\n",
        "                va=\"center\",\n",
        "                fontsize=7.5,\n",
        "            )\n",
        "        axm.annotate(Pa, ((i + 0.5) / 3, 1.08), ha=\"center\", fontsize=8.5)\n",
        "        axm.annotate(\n",
        "            \"XYZ\"[i],\n",
        "            (-0.13, 1 - (i + 0.5) / 3),\n",
        "            ha=\"center\",\n",
        "            va=\"center\",\n",
        "            fontsize=8.5,\n",
        "        )\n",
        "    axm.annotate(\"qubit a\", (0.5, 1.27), ha=\"center\", fontsize=9)\n",
        "    axm.annotate(\n",
        "        \"qubit b\", (-0.33, 0.5), rotation=90, va=\"center\", fontsize=9\n",
        "    )\n",
        "    axm.set_xlim(-0.35, 1.05)\n",
        "    axm.set_ylim(-0.05, 1.35)\n",
        "    axm.set_aspect(\"equal\")\n",
        "    axm.axis(\"off\")\n",
        "\n",
        "\n",
        "# --- run diagnostics ---------------------------------------------------------------\n",
        "\n",
        "\n",
        "def plot_trex(trex_rescale, tick_labels):\n",
        "    _fig, axt = plt.subplots(figsize=(12, 3))\n",
        "    axt.stem(np.arange(len(trex_rescale)), (trex_rescale - 1) * 100)\n",
        "    axt.set_xticks(np.arange(len(trex_rescale)), tick_labels)\n",
        "    axt.set_xlabel(\"Observable\")\n",
        "    axt.set_ylabel(\"Readout correction (%)\")\n",
        "    axt.set_title(\"TREX rescale factors\")\n",
        "    plt.tight_layout()\n",
        "\n",
        "\n",
        "def plot_postselection(mit):\n",
        "    \"\"\"Accepted shots per PEC randomization against a single binomial at the mean\n",
        "    acceptance rate: agreement means acceptance is independent of the sampled circuit\n",
        "    instance, the condition under which pooling accepted shots across randomizations\n",
        "    is a consistent estimator.\"\"\"\n",
        "    _fig, ax = plt.subplots(figsize=(7.5, 3.8))\n",
        "    counts = mit[\"acc_counts_post\"]\n",
        "    K, p = 64, counts.mean() / 64\n",
        "    ks = np.arange(K + 1)\n",
        "    pmf = np.array([comb(K, k) * p**k * (1 - p) ** (K - k) for k in ks])\n",
        "    ax.hist(\n",
        "        counts,\n",
        "        bins=np.arange(-0.5, K + 1.5),\n",
        "        density=True,\n",
        "        alpha=0.6,\n",
        "        color=\"#da1e28\",\n",
        "        label=\"measured\",\n",
        "    )\n",
        "    ax.plot(ks, pmf, \"k-\", lw=1.5, label=f\"Binomial(64, {p:.3f})\")\n",
        "    ax.set_xlim(-0.5, max(int(counts.max()) + 3, 20))\n",
        "    ax.set_xlabel(\"Accepted shots per randomization\")\n",
        "    ax.set_ylabel(\"Probability\")\n",
        "    ax.legend()\n",
        "    ax.grid(alpha=0.4)\n",
        "    plt.tight_layout()\n",
        "\n",
        "\n",
        "def plot_convergence(mit, obs_exact, n_data):\n",
        "    \"\"\"Running site-averaged estimate vs. randomizations for both PEC arms (the S5-consistent\n",
        "    signed-ratio estimator, evaluated on growing prefixes of the sweep).\"\"\"\n",
        "\n",
        "    def running(prefix):\n",
        "        bits = np.squeeze(\n",
        "            np.unpackbits(mit[f\"data_{prefix}\"], axis=-1)[..., :n_data]\n",
        "            ^ mit[f\"flips_{prefix}\"]\n",
        "        )\n",
        "        mask = np.squeeze(mit[f\"mask_{prefix}\"])\n",
        "        signs = 1 - 2 * (np.squeeze(mit[f\"signs_{prefix}\"]).sum(axis=-1) % 2)\n",
        "        qv = (1 - 2 * bits.astype(int)) * mit[\"trex_rescale\"]\n",
        "        u = (\n",
        "            (signs[:, None, None] * mask[..., None] * qv)\n",
        "            .sum(axis=1)\n",
        "            .mean(axis=1)\n",
        "        )\n",
        "        v = signs * mask.sum(axis=1)\n",
        "        Rs = np.arange(500, len(u) + 1, 500)\n",
        "        est, err = [], []\n",
        "        for R in Rs:\n",
        "            e = u[:R].sum() / v[:R].sum()\n",
        "            est.append(e)\n",
        "            err.append(\n",
        "                np.sqrt(((u[:R] - e * v[:R]) ** 2).sum()) / abs(v[:R].sum())\n",
        "            )\n",
        "        return Rs, np.array(est), np.array(err)\n",
        "\n",
        "    _fig, ax = plt.subplots(figsize=(12, 5))\n",
        "    ideal_avg = float(np.mean(obs_exact))\n",
        "    ax.axhline(ideal_avg, color=\"black\", label=\"ideal\")\n",
        "    ax.fill_between(\n",
        "        [-400, 24000],\n",
        "        ideal_avg - 0.025,\n",
        "        ideal_avg + 0.025,\n",
        "        color=\"grey\",\n",
        "        alpha=0.22,\n",
        "        label=r\"$\\pm 0.025$\",\n",
        "    )\n",
        "    for prefix, label, color in (\n",
        "        (\"van\", \"vanilla PEC\", \"#8a3ffc\"),\n",
        "        (\"post\", \"PEC + error detection\", \"#da1e28\"),\n",
        "    ):\n",
        "        Rs, est, err = running(prefix)\n",
        "        ax.errorbar(\n",
        "            Rs,\n",
        "            est,\n",
        "            yerr=err,\n",
        "            marker=\"o\",\n",
        "            linestyle=\"\",\n",
        "            markerfacecolor=\"none\",\n",
        "            color=color,\n",
        "            alpha=0.85,\n",
        "            capsize=3,\n",
        "            label=label,\n",
        "        )\n",
        "    ax.set_xlim(-400, 24000)\n",
        "    ax.set_ylim(ideal_avg - 0.08, ideal_avg + 0.08)\n",
        "    ax.set_xlabel(\"# randomizations\")\n",
        "    ax.set_ylabel(r\"Site-averaged $\\langle X \\rangle$\")\n",
        "    ax.legend(ncols=2)\n",
        "\n",
        "\n",
        "def plot_final(\n",
        "    obs_exact, baseline, ed, pec, post, gammas, tick_labels, title\n",
        "):\n",
        "    \"\"\"Per-site <X> for every method, with an rms-deviation inset.\n",
        "    baseline/ed/pec/post are (values, errors) pairs; gammas is (gamma, gamma_post).\"\"\"\n",
        "    x = np.arange(len(obs_exact))\n",
        "    _fig, ax = plt.subplots(figsize=(13, 5))\n",
        "    h_ideal = ax.errorbar(\n",
        "        x, obs_exact, fmt=\"-\", capsize=4, label=\"ideal\", color=\"black\"\n",
        "    )\n",
        "    h_base = ax.errorbar(\n",
        "        x,\n",
        "        baseline[0],\n",
        "        yerr=baseline[1],\n",
        "        fmt=\".--\",\n",
        "        capsize=4,\n",
        "        label=\"baseline\",\n",
        "        color=\"#0f62fe\",\n",
        "    )\n",
        "    h_ed = ax.errorbar(\n",
        "        x,\n",
        "        ed[0],\n",
        "        yerr=ed[1],\n",
        "        fmt=\"^\",\n",
        "        capsize=4,\n",
        "        label=\"error detection\",\n",
        "        color=\"#009d9a\",\n",
        "        alpha=0.8,\n",
        "    )\n",
        "    h_pec = ax.errorbar(\n",
        "        x,\n",
        "        pec[0],\n",
        "        yerr=pec[1],\n",
        "        fmt=\"x\",\n",
        "        capsize=4,\n",
        "        label=f\"PEC ($\\\\gamma$={gammas[0]:.0f})\",\n",
        "        color=\"#8a3ffc\",\n",
        "        alpha=0.7,\n",
        "    )\n",
        "    h_post = ax.errorbar(\n",
        "        x,\n",
        "        post[0],\n",
        "        yerr=post[1],\n",
        "        fmt=\"d\",\n",
        "        capsize=4,\n",
        "        label=f\"PEC + error detection ($\\\\gamma$={gammas[1]:.0f})\",\n",
        "        color=\"#da1e28\",\n",
        "    )\n",
        "    ax.set_xticks(x, tick_labels)\n",
        "    ax.set_xlabel(\"Observable\")\n",
        "    ax.set_ylabel(\"Expectation value\")\n",
        "    # Set the y-axis range using the ideal, baseline, error detection, and PEC + error\n",
        "    # detection curves, including error bars. Exclude PEC-only values when setting the\n",
        "    # range so large fluctuations do not make the other curves difficult to distinguish.\n",
        "    # PEC-only values may fall outside the visible range. Leave room above for the inset.\n",
        "    series = [np.asarray(obs_exact)] + [\n",
        "        np.asarray(v) + s * np.asarray(e)\n",
        "        for v, e in (baseline, ed, post)\n",
        "        for s in (-1, 1)\n",
        "    ]\n",
        "    lo, hi = min(a.min() for a in series), max(a.max() for a in series)\n",
        "    span = max(hi - lo, 0.1)\n",
        "    ax.set_ylim(lo - 0.1 * span, hi + 1.05 * span)\n",
        "    # legend ordered to match the curves' vertical positions in the chart\n",
        "    ax.legend(\n",
        "        handles=[h_ideal, h_pec, h_post, h_ed, h_base],\n",
        "        ncols=5,\n",
        "        loc=\"lower center\",\n",
        "        bbox_to_anchor=(0.5, 1.02),\n",
        "        frameon=False,\n",
        "    )\n",
        "    ax.set_title(title, pad=44)\n",
        "\n",
        "    # inset: rows bottom-to-top so it reads top-to-bottom: PEC, PEC+ED, ED, baseline\n",
        "    axi = ax.inset_axes([0.36, 0.68, 0.28, 0.26])\n",
        "    methods = [\n",
        "        (\"baseline\", baseline[0], \"#0f62fe\"),\n",
        "        (\"QED\", ed[0], \"#009d9a\"),\n",
        "        (\"PEC+QED\", post[0], \"#da1e28\"),\n",
        "        (\"PEC\", pec[0], \"#8a3ffc\"),\n",
        "    ]\n",
        "    for k, (_nm, vals, color) in enumerate(methods):\n",
        "        axi.barh(\n",
        "            k,\n",
        "            np.sqrt(np.mean((vals - np.array(obs_exact)) ** 2)),\n",
        "            color=color,\n",
        "            alpha=0.9,\n",
        "        )\n",
        "    axi.set_yticks(range(4), [m[0] for m in methods], fontsize=10)\n",
        "    axi.set_title(\"RMS deviation from ideal\", fontsize=11)\n",
        "    axi.tick_params(labelsize=10)\n",
        "    axi.patch.set_alpha(1.0)\n",
        "    axi.set_zorder(5)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4a7c1057-0ffd-4ce3-804c-bfb00d11dde0",
      "metadata": {},
      "source": [
        "<span id=\"classically-simulate-$langle-x-rangle_i$-for-each-site-in-the-22-qubit-hex-ising-model\" />\n",
        "\n",
        "## Simulación clásica de « $\\langle X \\rangle_i$ » para cada sitio en el modelo hex-Ising de 22 qubits\n",
        "\n",
        "En primer lugar, utilizamos la clase `Statevector` de Qiskit para simular los valores objetivo exactos de nuestro experimento. Aunque el problema de 22 qubits que estamos resolviendo es manejable desde el punto de vista clásico, el simulador de vectores de estado de Qiskit no se adapta a demostraciones mucho más grandes que esta, y para ampliar la escala más allá de unos 50 qubits, tendríamos que utilizar técnicas de simulación inexactas, como [la propagación de](/docs/addons/pauli-prop) Pauli. Este experimento de 49 qubits nos permite estudiar la combinación de la detección de errores con la mitigación de errores a escalas más grandes, al tiempo que se mantiene el acceso a los valores esperados ideales.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "0c82c211",
      "metadata": {
        "execution": {
          "iopub.execute_input": "2026-09-04T05:06:09.591394Z",
          "iopub.status.busy": "2026-09-04T05:06:09.591303Z",
          "iopub.status.idle": "2026-09-04T05:06:14.417462Z",
          "shell.execute_reply": "2026-09-04T05:06:14.417054Z"
        }
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/probabilistic-error-cancellation-with-logical-noise-models/extracted-outputs/0c82c211-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "import warnings\n",
        "\n",
        "import networkx as nx\n",
        "import numpy as np\n",
        "from qiskit.circuit import QuantumCircuit\n",
        "from qiskit.quantum_info import Pauli, Statevector\n",
        "\n",
        "# Silence two harmless upstream warnings (a Samplomatic default-change notice and a\n",
        "# numpy datetime timezone notice from the noise-learning circuit generator)\n",
        "warnings.filterwarnings(\n",
        "    \"ignore\", message=\"The default of the 'inject_noise_site'\"\n",
        ")\n",
        "warnings.filterwarnings(\n",
        "    \"ignore\", message=\"no explicit representation of timezones\"\n",
        ")\n",
        "\n",
        "# Hexagonal lattice\n",
        "data_graph = nx.convert_node_labels_to_integers(\n",
        "    nx.hexagonal_lattice_graph(2, 3)\n",
        ")\n",
        "\n",
        "# Number of Trotter steps and rotation angles\n",
        "depth = 4\n",
        "zz_coeff = -np.pi / 4\n",
        "x_coeff = 3 * np.pi / 8\n",
        "n_data = data_graph.order()\n",
        "n_checks = data_graph.size()\n",
        "n_qubits = n_data + n_checks\n",
        "\n",
        "qc_data = QuantumCircuit(n_data)\n",
        "qc_data.h(range(n_data))\n",
        "for _ in range(depth):\n",
        "    for edge in data_graph.edges:\n",
        "        qc_data.rzz(zz_coeff, *edge)\n",
        "    qc_data.rx(x_coeff, range(n_data))\n",
        "psi_exact = Statevector(qc_data)\n",
        "\n",
        "# X on site i = qubit i\n",
        "observables = [\"I\" * (n_data - 1 - i) + \"X\" + \"I\" * i for i in range(n_data)]\n",
        "obs_exact = [psi_exact.expectation_value(Pauli(o)).real for o in observables]\n",
        "\n",
        "plot_exact(\n",
        "    obs_exact,\n",
        "    x_labels(range(n_data), n_data),\n",
        "    f\"Statevector simulation, {n_data} qubit hex-Ising, {depth} Trotter steps\",\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "77e23625-91dc-4e0f-bb9a-72ad5066b084",
      "metadata": {},
      "source": [
        "<span id=\"implement-the-49-qubit-error-detecting-hex-ising-model-with-a-quantum-circuit-and-transpile-to-qpu-backend\" />\n",
        "\n",
        "## Implementar el modelo hex-Ising de detección de errores de 49 qubits mediante un circuito cuántico y transpilarlo al backend de la QPU\n",
        "\n",
        "Este circuito simula cuatro pasos de Trotter de un hamiltoniano de Ising de campo transversal de 22 qubits en una red hexagonal. Los 22 qubits de datos están integrados en un subgrafo «heavy-hex» de 49 qubits de la QPU Heron r3`ibm_boston`. Los 27 qubits auxiliares se utilizan con dos fines: mediar en el entrelazamiento entre los qubits de Ising y detectar errores lógicos durante la ejecución del circuito. Las capas de entrelazamiento del circuito se han implementado de tal forma que los qubits mediadores vuelvan al estado fundamental $|0\\rangle$ antes de que se realicen las mediciones finales. Si uno o más qubits mediadores registran un valor de « $|1\\rangle$ », esto indica que la muestra se ha visto afectada por un error lógico. Descartar estas muestras puede mejorar la fidelidad de la distribución muestreada.\n",
        "\n",
        "En el gráfico que aparece a continuación, los qubits verdes representan los 22 qubits de Ising, y los qubits naranjas representan los 27 qubits mediadores.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "3675348f",
      "metadata": {
        "execution": {
          "iopub.execute_input": "2026-09-04T05:06:14.418804Z",
          "iopub.status.busy": "2026-09-04T05:06:14.418667Z",
          "iopub.status.idle": "2026-09-04T05:06:29.475545Z",
          "shell.execute_reply": "2026-09-04T05:06:29.475067Z"
        }
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/probabilistic-error-cancellation-with-logical-noise-models/extracted-outputs/3675348f-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 3,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "from qiskit.circuit import ClassicalRegister\n",
        "from qiskit.transpiler import generate_preset_pass_manager\n",
        "from qiskit_ibm_runtime import QiskitRuntimeService\n",
        "\n",
        "backend_name = \"ibm_boston\"\n",
        "\n",
        "\n",
        "def generate_ed_ising(graph, depth, zz_coeff, x_coeff):\n",
        "    \"\"\"Build the mediated error-detecting Ising circuit for data-qubit `graph`; returns (circuit, CZ layers, hardware graph).\"\"\"\n",
        "    hw_graph = nx.Graph()\n",
        "    for i, (a, b) in enumerate(graph.edges()):\n",
        "        hw_graph.add_edges_from(\n",
        "            [(a, i + graph.order()), (b, i + graph.order())]\n",
        "        )\n",
        "    coloring = nx.coloring.greedy_color(\n",
        "        nx.line_graph(hw_graph), strategy=\"DSATUR\"\n",
        "    )\n",
        "    layers_coupling = [\n",
        "        [e for e, c in coloring.items() if c == i]\n",
        "        for i in set(coloring.values())\n",
        "    ]\n",
        "    circuit = QuantumCircuit(hw_graph.order())\n",
        "    circuit.h(range(hw_graph.order()))\n",
        "    for _ in range(depth):\n",
        "        for angle, qubits in (\n",
        "            (zz_coeff, range(graph.order(), hw_graph.order())),\n",
        "            (x_coeff, range(graph.order())),\n",
        "        ):\n",
        "            circuit.barrier()\n",
        "            for layer in layers_coupling:\n",
        "                for edge in layer:\n",
        "                    circuit.cz(*edge)\n",
        "            circuit.barrier()\n",
        "            circuit.rx(angle, qubits)\n",
        "    circuit.barrier()\n",
        "    circuit.h(range(graph.order(), hw_graph.order()))\n",
        "    return circuit, layers_coupling, hw_graph\n",
        "\n",
        "\n",
        "service = QiskitRuntimeService()\n",
        "backend = service.backend(backend_name)\n",
        "\n",
        "circuit, layers_coupling, hw_graph = generate_ed_ising(\n",
        "    data_graph, depth, zz_coeff, x_coeff\n",
        ")\n",
        "\n",
        "circ_meas = circuit.copy()\n",
        "data_reg, check_reg = (\n",
        "    ClassicalRegister(n_data, \"data\"),\n",
        "    ClassicalRegister(n_checks, \"check\"),\n",
        ")\n",
        "circ_meas.add_register(data_reg, check_reg)\n",
        "circ_meas.barrier()\n",
        "circ_meas.measure(range(n_data), data_reg)\n",
        "circ_meas.measure(range(n_data, n_qubits), check_reg)\n",
        "\n",
        "# Physical qubits hosting the 22 Ising sites on ibm_boston, one heavy-hex row per lattice chain;\n",
        "# the mediator for each edge is the physical qubit sitting between its two data qubits\n",
        "data_layout = [\n",
        "    *[95, 93, 91, 89, 87],\n",
        "    *[115, 113, 111, 109, 107, 105],\n",
        "    *[135, 133, 131, 129, 127, 125],\n",
        "    *[155, 153, 151, 149, 147],\n",
        "]\n",
        "coupling_graph = nx.Graph(list(backend.coupling_map.get_edges()))\n",
        "layout = data_layout + [\n",
        "    next(\n",
        "        iter(\n",
        "            nx.common_neighbors(\n",
        "                coupling_graph, data_layout[a], data_layout[b]\n",
        "            )\n",
        "        )\n",
        "    )\n",
        "    for a, b in data_graph.edges()\n",
        "]\n",
        "\n",
        "circ_trans = generate_preset_pass_manager(\n",
        "    backend=backend, optimization_level=1, initial_layout=layout\n",
        ").run(circ_meas)\n",
        "\n",
        "plot_layout(backend, layout, n_data)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d4bc5be6",
      "metadata": {},
      "source": [
        "<span id=\"specify-the-entangling-layers-in-the-circuit-with-samplomatic\" />\n",
        "\n",
        "## Especifica las capas de entrelazamiento del circuito con `samplomatic`.\n",
        "\n",
        "El paso `generate_boxing_pass_manager` «Samplomatic» agrupa las capas entrelazadas del circuito en recuadros anotados para el giro de Pauli y la inyección de ruido PEC. El circuito modelo resultante y el «samplex» (una distribución paramétrica sobre el circuito modelo que especifica sus giros y la inyección de ruido) se utilizan para definir y ejecutar todos los experimentos de aprendizaje con ruido y muestreo en la QPU. Las cajas de medición también llevan una anotación `ChangeBasis` , por lo que las rotaciones que miden los qubits de datos en la base X se suministran al samplex como entrada en el momento del muestreo, en lugar de integrarse en el circuito.\n",
        "\n",
        "A continuación, se aplica el mismo «boxing pass» a una versión reducida del circuito de Ising con detección de errores, para que podamos visualizar la estructura de un circuito «encajonado».\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "28dd30dc",
      "metadata": {
        "execution": {
          "iopub.execute_input": "2026-09-04T05:06:29.477553Z",
          "iopub.status.busy": "2026-09-04T05:06:29.477172Z",
          "iopub.status.idle": "2026-09-04T05:06:29.750482Z",
          "shell.execute_reply": "2026-09-04T05:06:29.749989Z"
        }
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/probabilistic-error-cancellation-with-logical-noise-models/extracted-outputs/28dd30dc-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 4,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "from samplomatic.builders import build\n",
        "from samplomatic.transpiler import generate_boxing_pass_manager\n",
        "from samplomatic.utils import find_unique_box_instructions\n",
        "\n",
        "boxed = generate_boxing_pass_manager(\n",
        "    enable_gates=True,\n",
        "    enable_measures=True,\n",
        "    inject_noise_targets=\"gates\",\n",
        "    inject_noise_strategy=\"individual_modification\",\n",
        "    inject_noise_site=\"after\",\n",
        "    twirling_strategy=\"active_circuit\",\n",
        "    measure_annotations=\"all\",\n",
        ").run(circ_trans)\n",
        "template, samplex = build(boxed)\n",
        "unique_instructions = find_unique_box_instructions(\n",
        "    boxed, normalize_annotations=None, undress_boxes=True\n",
        ")\n",
        "\n",
        "# Measure X on the data qubits and Z on the mediators\n",
        "basis_key = next(\n",
        "    s.name\n",
        "    for s in samplex.inputs().get_specs()\n",
        "    if s.name.startswith(\"basis_changes.\")\n",
        ")\n",
        "meas_basis = np.array(\n",
        "    [2 if q in set(layout[:n_data]) else 1 for q in sorted(layout)],\n",
        "    dtype=np.uint8,\n",
        ")\n",
        "\n",
        "draw_toy_circuit(generate_ed_ising, zz_coeff, x_coeff, include_checks=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8ed4dfbb",
      "metadata": {},
      "source": [
        "<span id=\"add-non-markovian-error-checks-to-the-circuit\" />\n",
        "\n",
        "## Añadir comprobaciones de errores no markovianos al circuito\n",
        "\n",
        "A continuación, añadimos [comprobaciones de errores no markovianos](/docs/addons/qiskit-mitigation/guides/postselection-with-non-markovian-error-checks) al circuito con el fin de protegernos contra el ruido que no se ha modelado en nuestro protocolo de aprendizaje. Estas comprobaciones funcionan mediante la aplicación de un «bit-flip» con pulso largo y la verificación de que el qubit haya pasado correctamente de un estado clásico a otro. Si los dos qubits de un borde no superan la comprobación, la muestra se descarta. Las comprobaciones de errores no markovianas pueden utilizarse al principio o al final del circuito (o en ambos); en este caso, las utilizamos únicamente al final. Además, colocamos qubits de «espectador» sin utilizar junto al circuito, lo que proporciona una mayor protección frente al ruido que estos qubits están diseñados para detectar.\n",
        "\n",
        "Recuerda que, al igual que en las comprobaciones de simetría, esto también es una forma de postselección, por lo que debemos asegurarnos de tener en cuenta la disminución de la tasa de postselección al combinar estas técnicas. En concreto, si $P(\\text{non-Markovian error}) = \\alpha$ y $P(\\text{symmetry error}) = \\beta$, entonces se necesitan aproximadamente $\\frac{1}{(1-\\alpha)(1-\\beta)}$ muestras con ruido para recuperar una muestra lógica.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "686235cc",
      "metadata": {
        "execution": {
          "iopub.execute_input": "2026-09-04T05:06:29.752265Z",
          "iopub.status.busy": "2026-09-04T05:06:29.752154Z",
          "iopub.status.idle": "2026-09-04T05:06:29.986205Z",
          "shell.execute_reply": "2026-09-04T05:06:29.985860Z"
        }
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/probabilistic-error-cancellation-with-logical-noise-models/extracted-outputs/686235cc-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 5,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "from qiskit.transpiler import PassManager\n",
        "from qiskit_mitigation.postselection import PostSelector\n",
        "from qiskit_mitigation.postselection.passes import (\n",
        "    AddPostCircuitNonMarkovianErrorChecks,\n",
        "    AddSpectatorPostCircuitNonMarkovianErrorChecks,\n",
        ")\n",
        "\n",
        "add_checks = PassManager(\n",
        "    [\n",
        "        AddPostCircuitNonMarkovianErrorChecks(x_pulse_type=\"xslow\"),\n",
        "        AddSpectatorPostCircuitNonMarkovianErrorChecks(\n",
        "            backend.coupling_map, x_pulse_type=\"xslow\"\n",
        "        ),\n",
        "    ]\n",
        ")\n",
        "template_checked = add_checks.run(template)\n",
        "selector = PostSelector.from_circuit(template_checked, backend.coupling_map)\n",
        "\n",
        "draw_toy_circuit(generate_ed_ising, zz_coeff, x_coeff, include_checks=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "421379f5-c277-45e8-86db-9922d675701b",
      "metadata": {},
      "source": [
        "<span id=\"learn-the-gate-and-readout-noise-for-the-49-qubit-checked-ising-circuit\" />\n",
        "\n",
        "## Averigua el ruido de puerta y de lectura del circuito Ising comprobado de 49 qubits\n",
        "\n",
        "Para mitigar el ruido, debemos modelar cómo afecta este a nuestras puertas de entrelazamiento. Aquí aprendemos un modelo de ruido de Pauli-Lindblad para cada una de las tres capas de entrelazamiento únicas del circuito utilizando [qiskit-noise-learning](https://github.com/Qiskit/qiskit-noise-learning). Más adelante, eliminaremos de este modelo de ruido los generadores de error que sean detectables mediante las comprobaciones de simetría y mitigaremos únicamente el canal de ruido restante mediante PEC. El siguiente mapa muestra las tres capas aprendidas en la QPU; cada qubit es una rueda con sus tasas X/Y/Z weight-1, y cada acoplador es una cuadrícula de 3×3 con las tasas de Pauli weight-2 de ese par. El valor « $\\gamma$ » de « $10$ » implica una sobrecarga de muestreo de « $\\gamma^2=10^2=100$ ».\n",
        "\n",
        "Las correcciones de lectura (TREX) proceden directamente de las trayectorias SPAM del ajuste basado en el aprendizaje del ruido. El gráfico de tallos situado debajo del mapa de ruido de la QPU muestra los factores de reescalado de TREX por observable.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "d6b90573",
      "metadata": {
        "execution": {
          "iopub.execute_input": "2026-09-04T05:06:29.987523Z",
          "iopub.status.busy": "2026-09-04T05:06:29.987451Z",
          "iopub.status.idle": "2026-09-04T05:52:16.828241Z",
          "shell.execute_reply": "2026-09-04T05:52:16.827803Z"
        }
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "learning shots removed by checks: 0.2393\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/probabilistic-error-cancellation-with-logical-noise-models/extracted-outputs/d6b90573-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/probabilistic-error-cancellation-with-logical-noise-models/extracted-outputs/d6b90573-2.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "from qiskit.quantum_info import PauliLindbladMap, QubitSparsePauli\n",
        "from qiskit_ibm_runtime import Executor, Session\n",
        "from qiskit_ibm_runtime.quantum_program import QuantumProgram\n",
        "from qiskit_mitigation.trex import TREX\n",
        "from qiskit_noise_learning.analysis import (\n",
        "    ComputeObservables,\n",
        "    CurveFitObservables,\n",
        "    FlipPostSelect,\n",
        "    LeastSquaresSolve,\n",
        ")\n",
        "from qiskit_noise_learning.circuit_generator import ExecutorCircuitGenerator\n",
        "from qiskit_noise_learning.experiment_builder import (\n",
        "    BindFragmentDepths,\n",
        "    CompleteSequences,\n",
        "    EvenDepthVanillaPaths,\n",
        "    Experiment,\n",
        "    GenerateInstructionSequences,\n",
        "    IdentifyRelations,\n",
        "    MergeInstructionSequences,\n",
        "    SPAMPaths,\n",
        "    VanillaInstructionSequences,\n",
        ")\n",
        "from qiskit_noise_learning.gate_sets import QiskitGateSet\n",
        "from qiskit_noise_learning.models import PauliLindbladModel\n",
        "from qiskit_noise_learning.models.utils import split_pauli_lindblad_model\n",
        "from samplomatic.annotations import InjectNoise\n",
        "from samplomatic.utils import get_annotation\n",
        "\n",
        "gate_set = QiskitGateSet(\n",
        "    target=backend.target,\n",
        "    qubit_subset=sorted(\n",
        "        {\n",
        "            boxed.find_bit(q).index\n",
        "            for i in unique_instructions\n",
        "            for q in i.qubits\n",
        "        }\n",
        "    ),\n",
        ")\n",
        "ref_to_qubits = {}\n",
        "for inst in unique_instructions:\n",
        "    if ann := get_annotation(inst.operation, InjectNoise):\n",
        "        gate_set.add_box_as_gate(inst, name=ann.ref)\n",
        "        ref_to_qubits[ann.ref] = sorted(\n",
        "            boxed.find_bit(q).index for q in inst.qubits\n",
        "        )\n",
        "\n",
        "fidelity_model = PauliLindbladModel.k_local(\n",
        "    gate_set, gate_k={**{r: 2 for r in ref_to_qubits}, \"M\": 1, \"P\": 1}\n",
        ")\n",
        "experiment = (\n",
        "    EvenDepthVanillaPaths()\n",
        "    + VanillaInstructionSequences()\n",
        "    + IdentifyRelations()\n",
        "    + SPAMPaths()\n",
        "    + GenerateInstructionSequences()\n",
        "    + MergeInstructionSequences()\n",
        "    + CompleteSequences()\n",
        "    + BindFragmentDepths([2, 4, 8, 12])\n",
        ").run(Experiment(fidelity_model=fidelity_model, shots=384, randomizations=64))\n",
        "\n",
        "circuit_generator = ExecutorCircuitGenerator(\n",
        "    gate_set, pass_manager=add_checks\n",
        ")\n",
        "program_learn, data_mapper = circuit_generator.generate(experiment)\n",
        "\n",
        "session = Session(backend)\n",
        "fit = circuit_generator.collect(\n",
        "    Executor(session).run(program_learn).result(), data_mapper\n",
        ")\n",
        "fit = (\n",
        "    FlipPostSelect()\n",
        "    + ComputeObservables()\n",
        "    + CurveFitObservables()\n",
        "    + LeastSquaresSolve()\n",
        ").run(fit)\n",
        "print(\n",
        "    \"learning shots removed by checks:\",\n",
        "    f\"{float(fit.raw_data.datatree['0']['data_mask'].mean()):.4f}\",\n",
        ")\n",
        "\n",
        "# TREX factors from the fit's SPAM paths: the fit only identifies the product of\n",
        "# state-prep and measurement error, so hand TREX the composition of the two maps\n",
        "spam = fidelity_model.to_pauli_lindblad_maps(\n",
        "    fit.model_data, include_spam=True\n",
        ")\n",
        "spam_map = spam[\"P\"].compose(spam[\"M\"])\n",
        "z_terms = [\n",
        "    QubitSparsePauli((\"Z\", [layout[i]]), num_qubits=backend.num_qubits)\n",
        "    for i in range(n_data)\n",
        "]\n",
        "trex_rescale = np.array(\n",
        "    [TREX.calculate_trex_factor(spam_map, z) for z in z_terms]\n",
        ")\n",
        "\n",
        "# learned maps come back in backend qubit indexing; samplex wants box-local order\n",
        "plm = split_pauli_lindblad_model(fit.model).model\n",
        "noise_maps = {}\n",
        "for ref, m in plm.to_pauli_lindblad_maps(fit.model_data).items():\n",
        "    box = ref_to_qubits[ref]\n",
        "    noise_maps[ref] = PauliLindbladMap.from_sparse_list(\n",
        "        [\n",
        "            (p, tuple(box.index(q) for q in qs), r)\n",
        "            for p, qs, r in m.to_sparse_list()\n",
        "        ],\n",
        "        num_qubits=len(box),\n",
        "    )\n",
        "# Full-PEC cost: each learned layer acts twice per Trotter step (compute + uncompute)\n",
        "gamma = float(\n",
        "    np.exp(2 * depth * sum(2 * sum(m.rates) for m in noise_maps.values()))\n",
        ")\n",
        "\n",
        "# Collect the learned model in the form the figure helpers expect\n",
        "mit = dict(\n",
        "    **{\n",
        "        f\"rates_{i}\": np.asarray(m.rates)\n",
        "        for i, m in enumerate(noise_maps.values())\n",
        "    },\n",
        "    **{\n",
        "        f\"label_paulis_{i}\": np.array([p for p, _, _ in m.to_sparse_list()])\n",
        "        for i, m in enumerate(noise_maps.values())\n",
        "    },\n",
        "    **{\n",
        "        f\"label_qubits_{i}\": np.array(\n",
        "            [\n",
        "                list(qs) + [-1] * (2 - len(qs))\n",
        "                for _, qs, _ in m.to_sparse_list()\n",
        "            ]\n",
        "        )\n",
        "        for i, m in enumerate(noise_maps.values())\n",
        "    },\n",
        "    layout=np.array(layout),\n",
        "    trex_rescale=trex_rescale,\n",
        "    gammas=np.array([gamma, np.nan]),\n",
        ")\n",
        "draw_noise_map(mit, backend)\n",
        "plot_trex(mit[\"trex_rescale\"], x_labels(layout, n_data))"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "211b9e87-7f14-4108-bef9-18ba95128fed",
      "metadata": {},
      "source": [
        "<span id=\"prune-the-noise-model-of-terms-detectable-by-the-symmetry-checks\" />\n",
        "\n",
        "## Eliminar del modelo de ruido los términos detectables mediante las comprobaciones de simetría\n",
        "\n",
        "Cada comprobación de simetría permite detectar un subconjunto de generadores de errores en el modelo de ruido, lo que significa que no es necesario corregir esos errores. Aquí creamos una máscara sobre los generadores de ruido detectables utilizando la función `create_postselected_noise_mask` de `qiskit_mitigation`. Esta máscara se utilizará para depurar el modelo de ruido de los términos detectables antes de realizar el PEC. La aplicación del PEC al modelo de ruido podado puede proporcionar valores esperados convergentes por una fracción del coste de muestreo que se requeriría para aplicar el PEC al modelo de ruido completo, tal y como se ilustra en los mapas del modelo de ruido anteriores.\n",
        "\n",
        "El mapa que se muestra a continuación se ha elaborado con la misma escala de colores que el modelo entrenado anterior, conservando únicamente los generadores de error que las comprobaciones no pueden detectar. Observamos que se han eliminado del modelo un gran número de generadores de errores y que la sobrecarga de muestreo se ha reducido considerablemente. La sobrecarga de muestreo para ejecutar el PEC en el modelo de ruido podado es de $\\gamma^2=2.5^2\\approx6.3$, lo que supone una reducción de aproximadamente un 16x respecto a la sobrecarga necesaria para mitigar el canal de ruido completo.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "758392e9",
      "metadata": {
        "execution": {
          "iopub.execute_input": "2026-09-04T05:52:16.829410Z",
          "iopub.status.busy": "2026-09-04T05:52:16.829349Z",
          "iopub.status.idle": "2026-09-04T05:52:17.402729Z",
          "shell.execute_reply": "2026-09-04T05:52:17.402235Z"
        }
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "gamma PEC 10.08 | gamma QED+PEC 2.52\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/probabilistic-error-cancellation-with-logical-noise-models/extracted-outputs/758392e9-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "from qiskit.quantum_info import Pauli\n",
        "from qiskit_mitigation.noise import create_postselected_noise_mask\n",
        "\n",
        "SHOTS_PER_RAND = 64  # shots per PEC randomization\n",
        "N_RAND = 23_054  # PEC randomizations per experiment\n",
        "MAX_PAIR_RATE = 0.02  # Max value of any coupler's summed two-qubit error rate\n",
        "\n",
        "# The detectors are virtual Pauli-Z's representing the measurements on the mediator qubits\n",
        "# We can find the set of detectable noise generators in the model for each detector by conjugating\n",
        "# it backward through the circuit and calculating what generators it anti-commutes with.\n",
        "detectors = [\n",
        "    Pauli(\"I\" * (boxed.num_qubits - 1 - q) + \"Z\" + \"I\" * q)\n",
        "    for q in layout[n_data:]\n",
        "]\n",
        "# Detectable generators are handled by postselection; prune them from the model for PEC\n",
        "local_scales, gamma2_post = create_postselected_noise_mask(\n",
        "    boxed, noise_maps, detectors\n",
        ")\n",
        "gamma_post = gamma2_post**0.5\n",
        "print(f\"gamma PEC {gamma:.2f} | gamma QED+PEC {gamma_post:.2f}\")\n",
        "\n",
        "# Kill the job if any coupler's summed rate exceeds MAX_PAIR_RATE\n",
        "pair_lam = {}\n",
        "for m in noise_maps.values():\n",
        "    for _, qs, r in m.to_sparse_list():\n",
        "        if len(qs) == 2:\n",
        "            pair_lam[tuple(sorted(qs))] = (\n",
        "                pair_lam.get(tuple(sorted(qs)), 0.0) + r\n",
        "            )\n",
        "assert (\n",
        "    max(pair_lam.values()) < MAX_PAIR_RATE\n",
        "), f\"kill: degraded coupler {max(pair_lam, key=pair_lam.get)}\"\n",
        "\n",
        "# Detectability depends on circuit position, so record each site's 0/1 scales by layer\n",
        "site_ref = {\n",
        "    a.modifier_ref: a.ref\n",
        "    for inst in boxed.data\n",
        "    if inst.operation.name == \"box\"\n",
        "    and (a := get_annotation(inst.operation, InjectNoise))\n",
        "    and a.ref\n",
        "}\n",
        "sites = sorted(local_scales, key=lambda s: int(s[1:]))\n",
        "mit[\"site_scales\"] = np.stack([local_scales[s] for s in sites])\n",
        "mit[\"site_layer\"] = np.array(\n",
        "    [list(noise_maps).index(site_ref[s]) for s in sites]\n",
        ")\n",
        "mit[\"gammas\"] = np.array([gamma, gamma_post])\n",
        "draw_noise_map(mit, backend, reduced=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "e969d2ef-2ba4-4e9d-8462-fee5f91c0b00",
      "metadata": {},
      "source": [
        "<span id=\"sample-the-error-detecting-circuit\" />\n",
        "\n",
        "## Prueba el circuito de detección de errores\n",
        "\n",
        "Aquí analizamos el circuito de Ising de 49 qubits con detección de errores. Habilitamos el «Pauli twirling» y las comprobaciones de errores no markovianas, pero no realizamos muestreo PEC. Utilizaremos estas muestras para calcular los valores esperados de referencia, así como los valores esperados correspondientes únicamente a la detección de errores.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "aa3fba3c",
      "metadata": {
        "execution": {
          "iopub.execute_input": "2026-09-04T05:52:17.404375Z",
          "iopub.status.busy": "2026-09-04T05:52:17.404294Z",
          "iopub.status.idle": "2026-09-04T05:52:18.394296Z",
          "shell.execute_reply": "2026-09-04T05:52:18.393058Z"
        }
      },
      "outputs": [],
      "source": [
        "# Baseline = the ED+PEC template itself at noise scale 0, i.e. twirling only\n",
        "program_tw = QuantumProgram(shots=100)\n",
        "program_tw.append_samplex_item(\n",
        "    template_checked,\n",
        "    samplex=samplex,\n",
        "    shape=(1000, 1),\n",
        "    samplex_arguments={\n",
        "        \"pauli_lindblad_maps\": dict(noise_maps),\n",
        "        basis_key: meas_basis,\n",
        "        **{f\"noise_scales.{k}\": 0.0 for k in local_scales},\n",
        "    },\n",
        ")\n",
        "job_tw = Executor(session).run(program_tw)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "869cabc2-5a6e-497b-906c-5e2fa634f0a8",
      "metadata": {},
      "source": [
        "<span id=\"sample-the-error-detecting-circuit-with-pec\" />\n",
        "\n",
        "## Prueba el circuito de detección de errores con PEC\n",
        "\n",
        "Ahora volvemos a realizar un muestreo e incluimos aleatorizaciones PEC para mitigar el ruido que las comprobaciones de simetría no pueden detectar. Estas muestras se utilizarán para calcular los valores esperados mediante la detección de errores en combinación con el PEC. A continuación representamos gráficamente los tiros aceptados por aleatorización PEC frente a una distribución binomial con la tasa de aceptación media. Observamos que la aceptación no está correlacionada con la instancia del circuito muestreada, ya que QED+PEC solo inyecta estados de Pauli indetectables, que conmutan con las mediciones de cada qubit mediador.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "15e09e01",
      "metadata": {
        "execution": {
          "iopub.execute_input": "2026-09-04T05:52:18.398262Z",
          "iopub.status.busy": "2026-09-04T05:52:18.398004Z",
          "iopub.status.idle": "2026-09-04T06:07:43.151360Z",
          "shell.execute_reply": "2026-09-04T06:07:43.149954Z"
        }
      },
      "outputs": [],
      "source": [
        "def launch(session, template, samplex, samplex_args, n_rand, shots_per_rand):\n",
        "    \"\"\"Sample `n_rand` randomizations of `template` on `session` in <=100k-randomization jobs; returns the jobs.\"\"\"\n",
        "    jobs = []\n",
        "    for start in range(0, n_rand, 100_000):\n",
        "        program = QuantumProgram(shots=shots_per_rand)\n",
        "        program.append_samplex_item(\n",
        "            template,\n",
        "            samplex=samplex,\n",
        "            samplex_arguments=samplex_args,\n",
        "            shape=(min(100_000, n_rand - start), 1),\n",
        "        )\n",
        "        jobs.append(Executor(session).run(program))\n",
        "    return jobs\n",
        "\n",
        "\n",
        "# Reduced PEC\n",
        "args_post = {\n",
        "    \"pauli_lindblad_maps\": dict(noise_maps),\n",
        "    basis_key: meas_basis,\n",
        "    **{f\"noise_scales.{k}\": -1.0 for k in local_scales},\n",
        "    **{f\"local_scales.{k}\": v for k, v in local_scales.items()},\n",
        "}\n",
        "# PEC on the full noise model\n",
        "args_van = {\n",
        "    \"pauli_lindblad_maps\": dict(noise_maps),\n",
        "    basis_key: meas_basis,\n",
        "    **{f\"noise_scales.{k}\": -1.0 for k in local_scales},\n",
        "}\n",
        "\n",
        "# Run sampling jobs\n",
        "jobs = launch(\n",
        "    session, template_checked, samplex, args_post, N_RAND, SHOTS_PER_RAND\n",
        ")\n",
        "jobs_van = launch(\n",
        "    session, template_checked, samplex, args_van, N_RAND, SHOTS_PER_RAND\n",
        ")\n",
        "outs = [j.result()[0] for j in jobs]\n",
        "outs_van = [j.result()[0] for j in jobs_van]\n",
        "(out_tw,) = job_tw.result()\n",
        "session.close()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "16fe6335-2173-4fb0-aab7-dfd2f3ad6b46",
      "metadata": {},
      "source": [
        "Recoge las muestras y comprueba que las aleatorizaciones de la PEC sean conmutativas con las mediciones en los qubits mediadores y no afecten a las estadísticas de la poselección.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "id": "5317dcaf",
      "metadata": {
        "execution": {
          "iopub.execute_input": "2026-09-04T06:07:43.155309Z",
          "iopub.status.busy": "2026-09-04T06:07:43.155072Z",
          "iopub.status.idle": "2026-09-04T06:07:43.973520Z",
          "shell.execute_reply": "2026-09-04T06:07:43.973108Z"
        }
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/probabilistic-error-cancellation-with-logical-noise-models/extracted-outputs/5317dcaf-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "def bitflip_mask(out, selector):\n",
        "    \"\"\"Per-shot True/False for job result `out`: shot passes `selector`'s non-Markovian error checks.\"\"\"\n",
        "    regs = {\n",
        "        k: np.asarray(out[k])\n",
        "        for k in out\n",
        "        if not k.startswith((\"measurement_flips\", \"pauli_signs\"))\n",
        "    }\n",
        "    return selector.compute_mask(regs, \"edge\", mode=\"post\")\n",
        "\n",
        "\n",
        "def symmetry_mask(out):\n",
        "    \"\"\"Per-shot True/False for job result `out`: every twirl-corrected symmetry check reads 0.\"\"\"\n",
        "    return ~(\n",
        "        np.asarray(out[\"check\"]) ^ np.asarray(out[\"measurement_flips.check\"])\n",
        "    ).any(axis=-1)\n",
        "\n",
        "\n",
        "def keep_mask(out, selector):\n",
        "    \"\"\"Per-shot True/False for job result `out`: shot passes both check types.\"\"\"\n",
        "    return symmetry_mask(out) & bitflip_mask(out, selector)\n",
        "\n",
        "\n",
        "def fracs(out_list, selector):\n",
        "    \"\"\"Fractions of shots across the job results in `out_list` passing [no, non-Markovian error, symmetry, both] checks.\"\"\"\n",
        "    bf = np.concatenate([bitflip_mask(o, selector) for o in out_list])\n",
        "    sy = np.concatenate([symmetry_mask(o) for o in out_list])\n",
        "    return [1.0, bf.mean(), sy.mean(), (bf & sy).mean()]\n",
        "\n",
        "\n",
        "mask_post = np.concatenate([keep_mask(o, selector) for o in outs])\n",
        "\n",
        "# Collect the sampled data alongside the learned model, in the form the figure helpers expect\n",
        "mit.update(\n",
        "    ps_fracs=np.array(\n",
        "        [\n",
        "            fracs([out_tw], selector),\n",
        "            fracs(outs_van, selector),\n",
        "            fracs(outs, selector),\n",
        "        ]\n",
        "    ),\n",
        "    acc_counts_post=np.squeeze(mask_post).sum(axis=-1),\n",
        "    data_tw=np.asarray(out_tw[\"data\"]),\n",
        "    flips_tw=np.asarray(out_tw[\"measurement_flips.data\"]),\n",
        "    mask_tw=keep_mask(out_tw, selector),\n",
        "    data_post=np.packbits(\n",
        "        np.concatenate([np.asarray(o[\"data\"]) for o in outs]), axis=-1\n",
        "    ),\n",
        "    flips_post=np.concatenate(\n",
        "        [np.asarray(o[\"measurement_flips.data\"]) for o in outs]\n",
        "    ),\n",
        "    signs_post=np.concatenate([np.asarray(o[\"pauli_signs\"]) for o in outs]),\n",
        "    mask_post=mask_post,\n",
        "    data_van=np.packbits(\n",
        "        np.concatenate([np.asarray(o[\"data\"]) for o in outs_van]), axis=-1\n",
        "    ),\n",
        "    flips_van=np.concatenate(\n",
        "        [np.asarray(o[\"measurement_flips.data\"]) for o in outs_van]\n",
        "    ),\n",
        "    signs_van=np.concatenate(\n",
        "        [np.asarray(o[\"pauli_signs\"]) for o in outs_van]\n",
        "    ),\n",
        "    mask_van=np.concatenate([bitflip_mask(o, selector) for o in outs_van]),\n",
        ")\n",
        "\n",
        "\n",
        "# Plot the PEC bias check\n",
        "plot_postselection(mit)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0bdbdc3e-f75f-4667-951d-27e936a415c3",
      "metadata": {},
      "source": [
        "<span id=\"calculate-expectation-values-and-compare-strategies\" />\n",
        "\n",
        "## Calcular los valores esperados y comparar estrategias\n",
        "\n",
        "Por último, calculamos todos los valores esperados con la función auxiliar `executor_expectation_values` de `qiskit-mitigation`, que se encarga de aplicar por nosotros los cambios de signo debidos a la medición, las máscaras de postselección, los factores de reescalado TREX y los signos de cuasiprobabilidad.\n",
        "\n",
        "**Gráfico superior:** Observamos que ambas variantes del PEC convergen, pero el modelo QED+PEC alcanza la banda « $\\pm 0.025$ » con menos aleatorizaciones y barras de error más estrechas. El sesgo residual en el cálculo QED+PEC es inferior a 0.01 y puede atribuirse a una combinación de discrepancias en el modelo de ruido, la deriva de los qubits con el paso del tiempo y los efectos de fuentes de ruido no modeladas.\n",
        "\n",
        "**Gráfico inferior** : Estimación actualizada de la media por centro del « $\\langle X \\rangle$ » a medida que se acumulan las aleatorizaciones, utilizando el mismo número de disparos por aleatorización para ambas variantes de PEC. La detección de errores PEC+ alcanza la banda en unas pocas miles de aleatorizaciones; el método PEC solo, que soporta toda la sobrecarga de muestreo, requiere muchas más aleatorizaciones para converger.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "id": "11bbb9c9",
      "metadata": {
        "execution": {
          "iopub.execute_input": "2026-09-04T06:07:43.974738Z",
          "iopub.status.busy": "2026-09-04T06:07:43.974659Z",
          "iopub.status.idle": "2026-09-04T06:07:44.524608Z",
          "shell.execute_reply": "2026-09-04T06:07:44.524148Z"
        }
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "gamma PEC 10.08 | gamma QED+PEC 2.52 (model) / 2.51 (data) | QED survival 0.304 | QED+PEC survival 0.305 | mean QED+PEC err 0.0037\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/probabilistic-error-cancellation-with-logical-noise-models/extracted-outputs/11bbb9c9-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "from qiskit.quantum_info import SparsePauliOp\n",
        "from qiskit_mitigation.utils import executor_expectation_values\n",
        "\n",
        "gamma, gamma_post = mit[\"gammas\"]\n",
        "basis_map = {Pauli(\"X\" * n_data): [SparsePauliOp(o) for o in observables]}\n",
        "rescale = dict(zip(observables, mit[\"trex_rescale\"], strict=True))\n",
        "\n",
        "\n",
        "def evs(bits, basis_map, **kwargs):\n",
        "    \"\"\"Per-site (means, standard errors) from boolean shot data `bits` of shape (rands, 1, shots, n_data).\"\"\"\n",
        "    out = executor_expectation_values(\n",
        "        bits, basis_map, None, avg_axis=(0, 1), **kwargs\n",
        "    )\n",
        "    return np.array([m for m, _ in out]).ravel(), np.sqrt(\n",
        "        [v for _, v in out]\n",
        "    ).ravel()\n",
        "\n",
        "\n",
        "def unpack(mit, prefix, n_data):\n",
        "    \"\"\"Restore `mit[f\"data_{prefix}\"]` from packed bytes to booleans of shape (rands, 1, shots, n_data).\"\"\"\n",
        "    return np.unpackbits(mit[f\"data_{prefix}\"], axis=-1)[..., :n_data].astype(\n",
        "        bool\n",
        "    )\n",
        "\n",
        "\n",
        "unmit_tw, unmit_tw_err = evs(\n",
        "    mit[\"data_tw\"], basis_map, measurement_flips=mit[\"flips_tw\"]\n",
        ")\n",
        "ed_tw, ed_tw_err = evs(\n",
        "    mit[\"data_tw\"],\n",
        "    basis_map,\n",
        "    measurement_flips=mit[\"flips_tw\"],\n",
        "    postselect_mask=mit[\"mask_tw\"],\n",
        "    rescale_factors=rescale,\n",
        ")\n",
        "post, post_err = evs(\n",
        "    unpack(mit, \"post\", n_data),\n",
        "    basis_map,\n",
        "    measurement_flips=mit[\"flips_post\"],\n",
        "    pauli_signs=mit[\"signs_post\"],\n",
        "    postselect_mask=mit[\"mask_post\"],\n",
        "    rescale_factors=rescale,\n",
        ")\n",
        "pec, pec_err = evs(\n",
        "    unpack(mit, \"van\", n_data),\n",
        "    basis_map,\n",
        "    measurement_flips=mit[\"flips_van\"],\n",
        "    pauli_signs=mit[\"signs_van\"],\n",
        "    postselect_mask=mit[\"mask_van\"],\n",
        "    rescale_factors=rescale,\n",
        ")\n",
        "\n",
        "# effective ED+PEC overhead measured from the data: accepted shots per signed accepted shot\n",
        "signs_post = 1 - 2 * (np.squeeze(mit[\"signs_post\"]).sum(axis=-1) % 2)\n",
        "gamma_eff = (\n",
        "    mit[\"mask_post\"].sum()\n",
        "    / (signs_post * np.squeeze(mit[\"mask_post\"]).sum(axis=-1)).sum()\n",
        ")\n",
        "\n",
        "print(\n",
        "    f\"gamma PEC {gamma:.2f} | gamma QED+PEC {gamma_post:.2f} (model) / {gamma_eff:.2f} (data) | \"\n",
        "    f\"QED survival {mit['mask_tw'].mean():.3f} | QED+PEC survival {mit['mask_post'].mean():.3f} | \"\n",
        "    f\"mean QED+PEC err {post_err.mean():.4f}\"\n",
        ")\n",
        "plot_final(\n",
        "    obs_exact,\n",
        "    (unmit_tw, unmit_tw_err),\n",
        "    (ed_tw, ed_tw_err),\n",
        "    (pec, pec_err),\n",
        "    (post, post_err),\n",
        "    (gamma, gamma_post),\n",
        "    x_labels(layout, n_data),\n",
        "    \"\",\n",
        ")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "id": "3ef3ea12-3b91-407a-90ad-9302581d5343",
      "metadata": {
        "execution": {
          "iopub.execute_input": "2026-09-04T06:07:44.525762Z",
          "iopub.status.busy": "2026-09-04T06:07:44.525676Z",
          "iopub.status.idle": "2026-09-04T06:07:44.857803Z",
          "shell.execute_reply": "2026-09-04T06:07:44.857414Z"
        }
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/probabilistic-error-cancellation-with-logical-noise-models/extracted-outputs/3ef3ea12-3b91-407a-90ad-9302581d5343-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "plot_convergence(mit, obs_exact, n_data)"
      ]
    },
    {
      "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
}