{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "446c0650-d248-4ab1-8dc0-41921ea624da",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"論理ノイズモデルを用いた確率的誤差相殺\"\n",
        "description: \"49キュービットの六角形アイジング回路において、メディエーター・キュービットによるエラー検出とPECを組み合わせることで、PEC単独の場合に比べ、はるかに少ないサンプリングコストでノイズを低減する。\"\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",
        "# 論理ノイズモデルを用いた確率的誤差相殺\n",
        "\n",
        "推定*実行時間：Heron r3 プロセッサで28分（注：これはあくまで推定値です。 （実行時間は状況によって異なる場合があります。）*\n",
        "\n",
        "![ヘビー・ヘキサゴナル・クビット配置に埋め込まれた六角形アイジング格子。 degree-3 のクビット上にアイジングサイトが配置され、それらの間の辺にはメディエーター・クビットが配置されている。](https://quantum.cloud.ibm.com/docs/images/tutorials/probabilistic-error-cancellation-with-logical-noise-models/hex-ising.avif)\n",
        "このチュートリアルでは、 ***ヘビー・ヘキサゴナル*** ・クビットトポロジーを持つHeron QPUを用いて、 ***六角形***&#x683C;子上の22サイト・アイジングモデルの観測量を評価する。 Heron QPUを用いた多くのデモンストレーションでは、システムの接続性に合わせるため、重六角格子上で定義されたシステムに焦点が当てられています。しかし、六角形モデルは自然界に広く見られ、接続密度が高いため古典的なシミュレーションが困難であることから、研究対象としてより興味深い傾向があります。 QPUの量子ビットトポロジーは六角格子の接続性を直接サポートできないため、六角格子モデルを量子ビットトポロジーにどのように効率的に組み込むかを決定しなければならない。 この例では、アイジングモデルの各サイトを、ヘビーヘックスQPU格子の頂点に位置する量子ビットで表現します。 辺にある量子ビット（ ***メディエーター量子ビット*** ）を用いて、頂点にあるイジング量子ビット間の量子もつれを促進します。 さらに、メディエーター量子ビットは、常に基底状態 $|0\\rangle$ に戻ると予想されるような方法で量子もつれを実現している。メディエーター量子ビットの測定結果が $|1\\rangle$ となった場合、それは回路の実行中に論理エラーが発生したことを示している。メディエーター量子ビットでエラーが検出されなかったサンプルのみをポストセレクトすることで、ノイズを含む生の分布よりも高い忠実度を持つ可能性のある、より狭&#x3044;***論理サンプルの***&#x5206;布が得られる。 メディエーター量子ビットを測定することでエラーの発生を検出できるだけでなく、その量子ビットが回路全体を通じて検出可能な論理エラーを正確に特定することも可能です。 ノイズモデルから、メディエーター量子ビット&#x306E;***対称性チェック***&#x306B;よって検出可能なノイズ発生源を除去すると、簡略化されたノイズモデルが残り、これは確率的誤差相殺（PEC）などの手法を用いて低減することができる。\n",
        "\n",
        "このノートブックの例では、前述のエラー検出手法とPECを組み合わせることで、いずれかの手法を単独で用いる場合よりも効果的に量子ノイズに対処します。 49量子ビットを用いて22量子ビットのアイジングモデルを埋め込み、余りの27個のメディエーター量子ビットをエラー検出に用いる。 対称性チェックをすり抜けるノイズを低減するため、事後選択されたノイズチャネルに対してPECを実行します。 エラー検出とPECを組み合わせることに加え、TREXによる読み出しエラーの低減や非マルコフ型エラーチェックなどの手法を用いて、量子ノイズの影響をさらに抑制していきます。\n",
        "\n",
        "ワークフローは、以下のとおりです。\n",
        "\n",
        "1. 22クビットのヘックス・アイジングモデルについて、観測可能な期待値を古典的にシミュレートする\n",
        "2. 49キュービットのエラー検出型ヘックス・アイジングモデルを量子回路で実装し、QPUバックエンドへトランスパイルする\n",
        "3. 回路内の絡み合い層を で指定します `samplomatic`。 これらの絡み合った層に影響を及ぼすノイズについて理解し、その影響を軽減していきます。\n",
        "4. 回路に非マルコフ型のエラーチェックを追加する\n",
        "5. 回路に影響を与えるゲートノイズと読み出しノイズについて学ぶ\n",
        "6. 対称性チェックによって検出された項のノイズモデルを剪定する\n",
        "7. エラー検出回路の動作を確認してください。 すべての対称性および非マルコフ性エラーチェックに合格したサンプルのみをポストセレクトし、すべての期待値計算においてTREX読み出しの緩和策を適用する\n",
        "   * エラー検出回路のサンプル\n",
        "   * PEC を使用してエラー検出回路のサンプリングを行うが、対称性チェックに基づくポストセレクションは行わず、学習されたノイズチャネル全体を軽減する。\n",
        "   * PEC を使用して、エラー検出回路の動作を確認してください。 対称性チェックに基づいてポストセレクトを行い、そのチェックでは検出できないノイズのみを低減する。\n",
        "8. 期待値を計算し、戦略を比較する。 QED+PECは、いずれかの手法を単独で用いる場合よりも、実験に影響を与えるノイズをより効果的に打ち消し、PEC単独の場合に必要なショット数よりもはるかに少ないショット数で収束した期待値を導き出すことに留意されたい。\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1d09a34f",
      "metadata": {},
      "source": [
        "<span id=\"requirements\" />\n",
        "\n",
        "## 要件\n",
        "\n",
        "このチュートリアルを始める前に、以下のコマンドを実行して、必要なパッケージをすべてインストールしてください：\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": [
        "以下の折りたたまれたセルには、このノートブック全体で使用される図のヘルパーが定義されています。\n",
        "\n",
        "<Accordion>\n",
        "  <AccordionItem title=\"図の説明を表示するには、ここをクリックしてください\">\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",
        "## 22量子ビットのヘックス・アイジングモデルにおいて、各サイトごとに $\\langle X \\rangle_i$ を古典的にシミュレートする\n",
        "\n",
        "まず、Qiskit の `Statevector` クラスを使用して、実験のための正確な目標値をシミュレートします。 私たちが解こうとしている22量子ビットの問題は、古典的には解くことが可能ですが、Qiskitの状態ベクトルシミュレータでは、これよりはるかに大規模なデモには対応できません。50量子ビット程度を超える規模に拡張するには、 [パウリ伝播法](/docs/addons/pauli-prop)などの近似的なシミュレーション手法を用いる必要があります。 この49キュービットの実験により、理想的な期待値へのアクセスを維持しつつ、より大規模なシステムにおいてエラー検出とエラー軽減の組み合わせについて調査することが可能となる。\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",
        "## 49キュービットのエラー検出型ヘックス・アイジングモデルを量子回路で実装し、QPUバックエンドへトランスパイルする\n",
        "\n",
        "この回路は、六角格子上の22キュービット横磁場アイジング・ハミルトニアンの4つのトロッターステップをシミュレートする。 22個のデータ量子ビットは、Heron r3 QPUの49量子ビットからなるヘビーヘックス部分グラフに組み込まれている `ibm_boston`。 27個のアンシラ量子ビットは、イジング量子ビット間の量子もつれを仲介することと、回路実行中の論理エラーを検出することという2つの目的で使用される。 この回路の絡み合い層は、メディエーター量子ビットが終端測定の前に基底状態 $|0\\rangle$ に戻るように実装されている。 1つ以上のメディエーター量子ビットが $|1\\rangle$ を測定した場合、サンプルが論理エラーによって破損したことを示している。 これらのサンプルを排除することで、サンプリングされた分布の忠実度を高めることができます。\n",
        "\n",
        "下のグラフでは、緑色の量子ビットは22個のアイジング量子ビットを表し、オレンジ色の量子ビットは27個のメディエーター量子ビットを表しています。\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",
        "## 回路内の絡み合い層を で指定します `samplomatic`。\n",
        "\n",
        "「Samplomatic」パス `generate_boxing_pass_manager` では、回路の複雑に絡み合った層を、パウリ・トゥイールおよびPECノイズ注入用の注釈付きボックスにグループ化します。 こうして得られたテンプレート回路とサンプレックス（テンプレート回路上のパラメトリック分布であり、そのツイールやノイズ注入を規定するもの）を用いて、ノイズ学習およびサンプリングに関するすべてのQPU実験を定義・実行する。 測定ボックスには注釈 `ChangeBasis` も付いているため、X基底でデータ量子ビットを測定するための回転は、回路に組み込まれるのではなく、サンプリング時にサンプレックスへの入力として供給される。\n",
        "\n",
        "以下では、エラー検出用アイジング回路の極小モデルに同じボクシング処理を適用し、ボックス化された回路の構造を可視化しています。\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",
        "## 回路に非マルコフ型のエラーチェックを追加する\n",
        "\n",
        "次に、学習プロトコルでモデル化されていないノイズから回路を保護するために、 [非マルコフ型の誤り検出](/docs/addons/qiskit-mitigation/guides/postselection-with-non-markovian-error-checks)機能を回路に追加します。 これらのチェックは、長パルスのビット反転を実施し、量子ビットがある古典状態から別の古典状態へ正しく移行したことを確認することで機能します。 エッジ上の2つの量子ビットの両方がチェックに合格しなかった場合、そのサンプルは破棄される。 非マルコフ型エラーチェックは、回路の最初または最後（あるいはその両方）で使用できますが、ここでは最後のみに使用します。 また、回路に隣接する未使用の「スペクテーター」量子ビットに対してもチェックを行い、これらのチェックが検出するように設計されたノイズに対する検出範囲をさらに広げています。\n",
        "\n",
        "対称性チェックと同様に、これもポストセレクションの一種であることを忘れないでください。したがって、これらの手法を組み合わせる際には、ポストセレクション率の低下を確実に考慮に入れる必要があります。 具体的には、 $P(\\text{non-Markovian error}) = \\alpha$ かつ $P(\\text{symmetry error}) = \\beta$ である場合、1つの論理サンプルを復元するには、およそ $\\frac{1}{(1-\\alpha)(1-\\beta)}$ 個のノイズの混じったサンプルが必要となります。\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",
        "## 49キュービットのチェックド・アイジング回路のゲートノイズと読み出しノイズについて学ぶ\n",
        "\n",
        "ノイズの影響を軽減するには、ノイズがエンタングルメントゲートにどのような影響を与えているかをモデル化する必要があります。 ここでは、 [qiskit-noise-learning](https://github.com/Qiskit/qiskit-noise-learning) を使用して、回路を構成する 3 つの固有のエンタングルメント層それぞれについて、パウリ・リンドブラッドノイズモデルを学習させます。 その後、対称性チェックによって検出可能な誤差発生源をこのノイズモデルから除去し、残ったノイズチャネルのみをPECを用いて低減します。 下の図は、QPU上の3つの学習済みレイヤーを示しています。各クビットは、 weight-1 のX/Y/Zレートからなるホイールであり、各カプラーは、そのペアにおける weight-2 のパウリレートからなる3×3のグリッドです。 $10$ の $\\gamma$ 値は、 $\\gamma^2=10^2=100$ のサンプリングオーバーヘッドを意味します。\n",
        "\n",
        "読み出し（TREX）の補正値は、ノイズ学習フィットのSPAMパスから直接得られます。 QPUノイズマップの下にあるステムプロットには、観測変数ごとのTREXリスケール係数が示されています。\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",
        "## 対称性チェックによって検出可能な項のノイズモデルを剪定する\n",
        "\n",
        "各対称性チェックは、ノイズモデルに含まれるエラー発生源の一部を検出することができるため、それらのエラーを軽減する必要はありません。 ここでは、の `create_postselected_noise_mask` 関数を使用して、検出可能なノイズ発生源の上にマスクを作成します `qiskit_mitigation`。 このマスクは、PECを実行する前に、検出可能な項のノイズモデルを剪定するために使用されます。 上記のノイズモデルマップに示されているように、剪定されたノイズモデルを用いてPECを実行すると、完全なノイズモデルでPECを実行する場合に必要なサンプリングコストのほんの一部で、収束した期待値を得ることができます。\n",
        "\n",
        "以下の地図は、上記の学習済みモデルと同じ色スケールで描かれており、チェックでは検出できないエラー発生源のみを残しています。 モデルから多数のエラー発生源が除去され、サンプリングのオーバーヘッドが大幅に減少したことがわかります。 剪定されたノイズモデル上でPECを実行するためのサンプリングオーバーヘッドは $\\gamma^2=2.5^2\\approx6.3$ であり、これは完全なノイズチャネルを低減するために必要なオーバーヘッドに比べて、およそ 16x の削減となる。\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",
        "## エラー検出回路のサンプル\n",
        "\n",
        "ここでは、49キュービットのエラー検出型イジング回路をシミュレーションします。 パウリ・トゥイリングおよび非マルコフ型エラーチェックを有効にしていますが、PECサンプリングは実行しません。 これらのサンプルを用いて、ベースラインの期待値およびエラー検出専用の期待値を算出します。\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",
        "## PEC を使用した誤り検出回路のサンプル\n",
        "\n",
        "ここで、再度サンプリングを行い、対称性チェックでは検出できないノイズを軽減するために、PECランダム化を適用します。 これらのサンプルは、PECと誤り検出を組み合わせて期待値を計算するために使用されます。 以下では、PECランダム化ごとの合格ショット数を、平均合格率における二項分布と対比してプロットする。 QED+PECは、すべてのメディエーター量子ビットの測定と可換な、検出不可能なパウリ状態のみを注入するため、受容性はサンプリングされた回路インスタンスとは相関がないことがわかる。\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": [
        "サンプルを収集し、PECランダム化がメディエーター量子ビットの測定と可換であり、ポストセレクションの統計に影響を与えないことを確認する。\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",
        "## 期待値を計算し、戦略を比較する\n",
        "\n",
        "最後に、から提供されているヘルパー `executor_expectation_values` 関数を用いて `qiskit-mitigation`、すべての期待値を計算します。この関数は、測定による反転、ポストセレクションマスク、TREXの再スケーリング係数、および準確率の符号を自動的に適用してくれます。\n",
        "\n",
        "**上図：** どちらのPECバリアントも収束することが確認できるが、QED+PECの方が、より少ないランダム化回数で、より狭い誤差範囲をもって $\\pm 0.025$ バンドに到達している。 QED+PEC計算における残留バイアスは 0.01 を下回っており、これはノイズモデルの不一致、経時的な量子ビットのドリフト、およびモデル化されていないノイズ源による影響の組み合わせに起因すると考えられる。\n",
        "\n",
        "**下段のグラフ** ：両方のPECバリアントについて、無作為化ごとに同一のショット数を使用し、無作為化が蓄積されるにつれて算出される、サイト平均の $\\langle X \\rangle$ の推定値。 PECとエラー検出を組み合わせた場合は、数千回のランダム化処理のうちにバンド内に収まるようになる。一方、PECのみの場合、サンプリングのオーバーヘッドをすべて負担することになるため、収束するにははるかに多くのランダム化処理が必要となる。\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
}