{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "b6d1e3ec",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Q-CTRLの性能管理を備えた横磁場イジングモデル\"\n",
        "description: \"Fire OpalとQiskitを用いて横磁場イージングモデルを構築し解く\"\n",
        "---\n",
        "\n",
        "{/* cspell:ignore fontsize sharey rmse RMSE boxstyle */}\n",
        "\n",
        "<span id=\"transverse-field-ising-model-with-q-ctrls-performance-management\" />\n",
        "\n",
        "# Q-CTRLの性能管理を備えた横磁場イジングモデル\n",
        "\n"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "id": "a6f69b77",
      "metadata": {},
      "source": [
        "*使用時間の目安：Heron r2 プロセッサーで2分。 (注：これはあくまでも目安です。 実行時間は異なるかもしれない)。*\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8bf80006",
      "metadata": {},
      "source": [
        "<span id=\"background\" />\n",
        "\n",
        "## 背景\n",
        "\n",
        "横磁場イジング模型（TFIM）は、量子磁性や相転移の研究に重要である。 量子揺らぎを引き起こす外部磁場の影響を受けながら、各スピンが隣のスピンと相互作用する。\n",
        "\n",
        "このモデルをシミュレートする一般的なアプローチは、トロッター分解を使って時間発展演算子を近似し、1量子ビットの回転と2量子ビットの相互作用を交互に行う回路を構成することである。 しかし、実際のハードウェア上でのこのシミュレーションは、ノイズやデコヒーレンスにより、真のダイナミクスから乖離してしまうため困難である。 これを克服するために、Q-CTRLのFire Opalエラー抑制およびパフォーマンス管理ツールを使用しています。これはQiskit機能として提供されています（ [Fire Opalのドキュメントを](/docs/guides/q-ctrl-performance-management)参照）。 Fire Opalは、動的デカップリング、高度なレイアウト、配線、その他のエラー抑制技術を適用することにより、回路実行を自動的に最適化します。 これらの改善により、ハードウェアの結果はノイズレスシミュレーションにより近くなり、TFIMの磁化ダイナミクスをより忠実に研究できるようになった。\n",
        "\n",
        "このチュートリアルでは\n",
        "\n",
        "* スピン三角形のグラフ上にTFIMハミルトニアンを構築する\n",
        "* 異なる深さのトロッター化回路で時間発展をシミュレートする\n",
        "* $\\langle Z_i \\rangle$ 一量子ビットの磁化を経時的に計算し、可視化する\n",
        "* ベースラインシミュレーションと、Q-CTRLのFire Opalパフォーマンス管理を使用したハードウェア実行結果の比較\n",
        "\n",
        "<span id=\"overview\" />\n",
        "\n",
        "## 概要\n",
        "\n",
        "横場イジング模型（TFIM）は、量子相転移の本質的な特徴を捉えた量子スピン模型である。 ハミルトニアンは次のように定義される：\n",
        "\n",
        "$$\n",
        "H = -J \\sum_{i} Z_i Z_{i+1} - h \\sum_{i} X_i\n",
        "$$\n",
        "\n",
        "ここで、 $Z_i$ と $X_i$ は量子ビット $i$ に作用するパウリ作用素であり、 $J$ は隣り合うスピン間の結合強度、 $h$ は横磁場の強度である。 第一項は古典的な強磁性相互作用を表し、第二項は横場を通して量子揺らぎを導入する。 TFIMのダイナミクスをシミュレートするには、ユニタリー進化演算子のトロッター分解（ $e^{-iHt}$ ）を使用します。この演算子は、接続されたスピン三角形のカスタムグラフに基づくRXゲートとRZZゲートのレイヤーを通して実装されます。 このシミュレーションでは、トロッター段数の増加に伴って磁化 $\\langle Z \\rangle$ がどのように変化するかを調べている。\n",
        "\n",
        "ノイズのないシミュレーションとノイズのあるバックエンドを比較することで、提案するTFIM実装の性能を評価する。 Fire Opalの強化された実行機能とエラー抑制機能は、実際のハードウェアにおけるノイズの影響を軽減するために使用され、 $\\langle Z_i \\rangle$ や相関器 $\\langle Z_i Z_j \\rangle$ のようなスピン観測値のより信頼性の高い推定値をもたらします。\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "55b94021",
      "metadata": {},
      "source": [
        "<span id=\"requirements\" />\n",
        "\n",
        "## 要件\n",
        "\n",
        "このチュートリアルを始める前に、以下のものがインストールされていることを確認してください：\n",
        "\n",
        "* Qiskit SDK v1.4 またはそれ以降、 [可視化](/docs/api/qiskit/visualization)サポート付き\n",
        "* Qiskit Runtime v0.40 またはそれ以降 (`pip install qiskit-ibm-runtime`)\n",
        "* Qiskit Functions Catalog v0.9.0 (`pip install qiskit-ibm-catalog`)\n",
        "* Fire Opal SDK v9.0.2 またはそれ以降 (`pip install fire-opal`)\n",
        "* Q-CTRL Visualizer v8.0.2 またはそれ以降 (`pip install qctrl-visualizer`)\n",
        "\n"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "id": "7db2e559",
      "metadata": {},
      "source": [
        "<span id=\"setup\" />\n",
        "\n",
        "## セットアップ\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ecaa518f",
      "metadata": {},
      "source": [
        "まず、 [IBM Quantum](http://quantum.cloud.ibm.com/) のAPIキーを使用して認証を行ってください。 次に、次のようにQiskit関数を選択します。 （このコードは、 [アカウントが](/docs/guides/functions-get-started#install-qiskit-functions-catalog-client)すでにローカル環境に保存されていることを前提としています。）\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "bc380c46",
      "metadata": {},
      "outputs": [],
      "source": [
        "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
        "from qiskit import QuantumCircuit\n",
        "from qiskit_ibm_catalog import QiskitFunctionsCatalog\n",
        "from qiskit_ibm_runtime import QiskitRuntimeService\n",
        "from qiskit_ibm_runtime import SamplerV2 as Sampler\n",
        "from qiskit.quantum_info import SparsePauliOp\n",
        "from qiskit_aer import AerSimulator\n",
        "\n",
        "import numpy as np\n",
        "import networkx as nx\n",
        "import matplotlib.pyplot as plt\n",
        "import qctrlvisualizer as qv"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "e9d916b6",
      "metadata": {},
      "outputs": [],
      "source": [
        "catalog = QiskitFunctionsCatalog(channel=\"ibm_quantum_platform\")\n",
        "\n",
        "# Access Function\n",
        "perf_mgmt = catalog.load(\"q-ctrl/performance-management\")"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "id": "988ee237",
      "metadata": {},
      "source": [
        "<span id=\"step-1-map-classical-inputs-to-a-quantum-problem\" />\n",
        "\n",
        "## ステップ1：古典的な入力を量子問題にマッピングする\n",
        "\n",
        "<span id=\"generate-tfim-graph\" />\n",
        "\n",
        "### TFIMグラフを生成する\n",
        "\n",
        "まず、スピンの格子とスピン間の結合を定義する。 このチュートリアルでは、格子は直線状に連なった三角形で構成されている。 各三角形は3つのノードが閉じたループで結ばれており、各三角形の1つのノードを前の三角形に連結することで鎖が形成される。\n",
        "\n",
        "ヘルパー関数 `connected_triangles_adj_matrix` は、この構造の隣接行列を構築する。 $n$ 三角形の連鎖の場合、結果のグラフは $2n+1$ ノードを含む。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "829f1b6d",
      "metadata": {},
      "outputs": [],
      "source": [
        "def connected_triangles_adj_matrix(n):\n",
        "    \"\"\"\n",
        "    Generate the adjacency matrix for 'n' connected triangles in a chain.\n",
        "    \"\"\"\n",
        "    num_nodes = 2 * n + 1\n",
        "    adj_matrix = np.zeros((num_nodes, num_nodes), dtype=int)\n",
        "\n",
        "    for i in range(n):\n",
        "        a, b, c = i * 2, i * 2 + 1, i * 2 + 2  # Nodes of the current triangle\n",
        "\n",
        "        # Connect the three nodes in a triangle\n",
        "        adj_matrix[a, b] = adj_matrix[b, a] = 1\n",
        "        adj_matrix[b, c] = adj_matrix[c, b] = 1\n",
        "        adj_matrix[a, c] = adj_matrix[c, a] = 1\n",
        "\n",
        "        # If not the first triangle, connect to the previous triangle\n",
        "        if i > 0:\n",
        "            adj_matrix[a, a - 1] = adj_matrix[a - 1, a] = 1\n",
        "\n",
        "    return adj_matrix"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "53869b13",
      "metadata": {},
      "source": [
        "先ほど定義した格子を視覚化するには、三角形の連なりをプロットし、各ノードにラベルをつければよい。 下の関数は、選んだ数の三角形のグラフを作り、それを表示する。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "2fc293aa",
      "metadata": {},
      "outputs": [],
      "source": [
        "def plot_triangle_chain(n, side=1.0):\n",
        "    \"\"\"\n",
        "    Plot a horizontal chain of n equilateral triangles.\n",
        "    Baseline: even nodes (0,2,4,...,2n) on y=0\n",
        "    Apexes: odd nodes (1,3,5,...,2n-1) above the midpoint.\n",
        "    \"\"\"\n",
        "    # Build graph\n",
        "    A = connected_triangles_adj_matrix(n)\n",
        "    G = nx.from_numpy_array(A)\n",
        "\n",
        "    h = np.sqrt(3) / 2 * side\n",
        "    pos = {}\n",
        "\n",
        "    # Place baseline nodes\n",
        "    for k in range(n + 1):\n",
        "        pos[2 * k] = (k * side, 0.0)\n",
        "\n",
        "    # Place apex nodes\n",
        "    for k in range(n):\n",
        "        x_left = pos[2 * k][0]\n",
        "        x_right = pos[2 * k + 2][0]\n",
        "        pos[2 * k + 1] = ((x_left + x_right) / 2, h)\n",
        "\n",
        "    # Draw\n",
        "    fig, ax = plt.subplots(figsize=(1.5 * n, 2.5))\n",
        "    nx.draw(\n",
        "        G,\n",
        "        pos,\n",
        "        ax=ax,\n",
        "        with_labels=True,\n",
        "        font_size=10,\n",
        "        font_color=\"white\",\n",
        "        node_size=600,\n",
        "        node_color=qv.QCTRL_STYLE_COLORS[0],\n",
        "        edge_color=\"black\",\n",
        "        width=2,\n",
        "    )\n",
        "    ax.set_aspect(\"equal\")\n",
        "    ax.margins(0.2)\n",
        "    plt.show()\n",
        "\n",
        "    return G, pos"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8abc0622",
      "metadata": {},
      "source": [
        "このチュートリアルでは、20個の三角形の連鎖を使う。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "861ab6e3",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/transverse-field-ising-model/extracted-outputs/861ab6e3-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "n_triangles = 20\n",
        "n_qubits = 2 * n_triangles + 1\n",
        "plot_triangle_chain(n_triangles, side=1.0)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ac6f36e3",
      "metadata": {},
      "source": [
        "<span id=\"coloring-graph-edges\" />\n",
        "\n",
        "### グラフの辺の着色\n",
        "\n",
        "スピン-スピン結合を実装するには、重ならないエッジをグループ化することが有効である。 これにより、2量子ビットのゲートを並列に適用することができる。 これは単純なエッジ・カラーリング手順 [\\[1\\]](#references) で行うことができ、各エッジに色を割り当てることで、同じノードで出会うエッジが異なるグループに配置されるようにする。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "c893b88f",
      "metadata": {},
      "outputs": [],
      "source": [
        "def edge_coloring(graph):\n",
        "    \"\"\"\n",
        "    Takes a NetworkX graph and returns a list of lists\n",
        "    where each inner list contains\n",
        "    the edges assigned the same color.\n",
        "    \"\"\"\n",
        "    line_graph = nx.line_graph(graph)\n",
        "    edge_colors = nx.coloring.greedy_color(line_graph)\n",
        "\n",
        "    color_groups = {}\n",
        "    for edge, color in edge_colors.items():\n",
        "        if color not in color_groups:\n",
        "            color_groups[color] = []\n",
        "        color_groups[color].append(edge)\n",
        "\n",
        "    return list(color_groups.values())"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "id": "b4d480b3",
      "metadata": {},
      "source": [
        "<span id=\"step-2-optimize-problem-for-quantum-hardware-execution\" />\n",
        "\n",
        "## ステップ2：量子ハードウェア実行に向けた問題の最適化\n",
        "\n",
        "<span id=\"generate-trotterized-circuits-on-spin-graphs\" />\n",
        "\n",
        "### スピングラフ上でトロッター化回路を生成する\n",
        "\n",
        "TFIMのダイナミクスをシミュレートするために、時間発展演算子を近似する回路を構築する。\n",
        "\n",
        "$$\n",
        "U(t) = e^{-i H t}, \\quad \\text{where} \\quad H = -J \\sum_{\\langle i,j \\rangle} Z_i Z_j - h \\sum_i X_i .\n",
        "$$\n",
        "\n",
        "2次のトロッター分解を用いる：\n",
        "\n",
        "$$\n",
        "e^{-i H \\Delta t} \\approx e^{-i H_X \\Delta t / 2}\\, e^{-i H_Z \\Delta t}\\, e^{-i H_X \\Delta t / 2},\n",
        "$$\n",
        "\n",
        "ここで、 $H_X = -h \\sum_i X_i$ と $H_Z = -J \\sum_{\\langle i,j \\rangle} Z_i Z_j$。\n",
        "\n",
        "* $H_X$ 項は、 `RX` 回転のレイヤーで実装されている。\n",
        "* $H_Z$ 項は、相互作用グラフのエッジに沿って `RZZ` ゲートのレイヤーで実装されている。\n",
        "\n",
        "これらのゲートの角度は、横磁場 $h$、結合定数 $J$、時間ステップ $\\Delta t$ によって決定される。複数のトロッター・ステップを積み重ねることで、システムのダイナミクスを近似する深さの増す回路を生成する。 関数 `generate_tfim_circ_custom_graph` と `trotter_circuits` は、任意のスピン相互作用グラフからトロッター化量子回路を構成する。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "ebe5c712",
      "metadata": {},
      "outputs": [],
      "source": [
        "def generate_tfim_circ_custom_graph(\n",
        "    steps, h, J, dt, psi0, graph: nx.graph.Graph, meas_basis=\"Z\", mirror=False\n",
        "):\n",
        "    \"\"\"\n",
        "    Generate a second order trotter of the form e^(a+b) ~ e^(b/2) e^a e^(b/2)\n",
        "    for simulating a transverse field ising model:\n",
        "    e^{-i H t} where the Hamiltonian H = -J \\\\sum_i Z_i Z_{i+1} + h \\\\sum_i X_i.\n",
        "\n",
        "    steps: Number of trotter steps\n",
        "    theta_x: Angle for layer of X rotations\n",
        "    theta_zz: Angle for layer of ZZ rotations\n",
        "    theta_x: Angle for second layer of X rotations\n",
        "    J: Coupling between nearest neighbor spins\n",
        "    h: The transverse magnetic field strength\n",
        "    dt: t/total_steps\n",
        "    psi0: initial state (assumed to be prepared in the computational basis).\n",
        "    meas_basis: basis to measure all correlators in\n",
        "\n",
        "    This is a second order trotter of the form e^(a+b) ~ e^(b/2) e^a e^(b/2)\n",
        "    \"\"\"\n",
        "    theta_x = h * dt\n",
        "    theta_zz = -2 * J * dt\n",
        "    nq = graph.number_of_nodes()\n",
        "    color_edges = edge_coloring(graph)\n",
        "    circ = QuantumCircuit(nq, nq)\n",
        "    # Initial state, for typical cases in the computational basis\n",
        "    for i, b in enumerate(psi0):\n",
        "        if b == \"1\":\n",
        "            circ.x(i)\n",
        "    # Trotter steps\n",
        "    for step in range(steps):\n",
        "        for i in range(nq):\n",
        "            circ.rx(theta_x, i)\n",
        "        if mirror:\n",
        "            color_edges = [sublist[::-1] for sublist in color_edges[::-1]]\n",
        "        for edge_list in color_edges:\n",
        "            for edge in edge_list:\n",
        "                circ.rzz(theta_zz, edge[0], edge[1])\n",
        "        for i in range(nq):\n",
        "            circ.rx(theta_x, i)\n",
        "\n",
        "    # some typically used basis rotations\n",
        "    if meas_basis == \"X\":\n",
        "        for b in range(nq):\n",
        "            circ.h(b)\n",
        "    elif meas_basis == \"Y\":\n",
        "        for b in range(nq):\n",
        "            circ.sdg(b)\n",
        "            circ.h(b)\n",
        "\n",
        "    for i in range(nq):\n",
        "        circ.measure(i, i)\n",
        "\n",
        "    return circ\n",
        "\n",
        "\n",
        "def trotter_circuits(G, d_ind_tot, J, h, dt, meas_basis, mirror=True):\n",
        "    \"\"\"\n",
        "    Generates a sequence of Trotterized circuits, each with increasing depth.\n",
        "    Given a spin interaction graph and Hamiltonian parameters, it constructs\n",
        "    a list of circuits with 1 to d_ind_tot Trotter steps\n",
        "\n",
        "    G: Graph defining spin interactions (edges = ZZ couplings)\n",
        "    d_ind_tot: Number of Trotter steps (maximum depth)\n",
        "    J: Coupling between nearest neighboring spins\n",
        "    h: Transverse magnetic field strength\n",
        "    dt: (t / total_steps\n",
        "    meas_basis: Basis to measure all correlators in\n",
        "    mirror: If True, mirror the Trotter layers\n",
        "    \"\"\"\n",
        "    qubit_count = len(G)\n",
        "    circuits = []\n",
        "    psi0 = \"0\" * qubit_count\n",
        "\n",
        "    for steps in range(1, d_ind_tot + 1):\n",
        "        circuits.append(\n",
        "            generate_tfim_circ_custom_graph(\n",
        "                steps, h, J, dt, psi0, G, meas_basis, mirror\n",
        "            )\n",
        "        )\n",
        "    return circuits"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "50b94af2",
      "metadata": {},
      "source": [
        "<span id=\"estimate-single-qubit-magnetizations-$langle-z_i-rangle$\" />\n",
        "\n",
        "### 単一量子ビットの磁化を推定する $\\langle Z_i \\rangle$\n",
        "\n",
        "このモデルのダイナミクスを研究するために、期待値 $\\langle Z_i \\rangle = \\langle \\psi | Z_i | \\psi \\rangle$ で定義されるそれぞれの量子ビットの磁化を測定したい。\n",
        "\n",
        "シミュレーションでは、測定結果から直接計算することができる。 関数 `z_expectation` は、ビット列カウントを処理し、選択された量子ビットインデックスに対する $\\langle Z_i \\rangle$ の値を返す。 実際のハードウェアでは、関数 `generate_z_observables` を使ってパウリ演算子を指定することで同じ量を評価し、バックエンドが期待値を計算する。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "3fc929e0",
      "metadata": {},
      "outputs": [],
      "source": [
        "def z_expectation(counts, index):\n",
        "    \"\"\"\n",
        "    counts: Dict of mitigated bitstrings.\n",
        "    index: Index i in the single operator expectation value < II...Z_i...I >\n",
        "        to be calculated.\n",
        "    return:  < Z_i >\n",
        "    \"\"\"\n",
        "    z_exp = 0\n",
        "    tot = 0\n",
        "    for bitstring, value in counts.items():\n",
        "        bit = int(bitstring[index])\n",
        "        sign = 1\n",
        "        if bit % 2 == 1:\n",
        "            sign = -1\n",
        "        z_exp += sign * value\n",
        "        tot += value\n",
        "\n",
        "    return z_exp / tot"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "id": "c9e79716",
      "metadata": {},
      "outputs": [],
      "source": [
        "def generate_z_observables(nq):\n",
        "    observables = []\n",
        "    for i in range(nq):\n",
        "        pauli_string = \"\".join([\"Z\" if j == i else \"I\" for j in range(nq)])\n",
        "        observables.append(SparsePauliOp(pauli_string))\n",
        "    return observables"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "id": "598d11cf",
      "metadata": {},
      "outputs": [],
      "source": [
        "observables = generate_z_observables(n_qubits)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "397d9b05",
      "metadata": {},
      "source": [
        "ここで、トロッター化回路を生成するためのパラメータを定義する。 このチュートリアルでは、格子は20個の三角形が連なったもので、41qubitのシステムに相当する。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "id": "2853d21b",
      "metadata": {},
      "outputs": [],
      "source": [
        "all_circs_mirror = []\n",
        "for num_triangles in [n_triangles]:\n",
        "    for meas_basis in [\"Z\"]:\n",
        "        A = connected_triangles_adj_matrix(num_triangles)\n",
        "        G = nx.from_numpy_array(A)\n",
        "        nq = len(G)\n",
        "        d_ind_tot = 22\n",
        "        dt = 2 * np.pi * 1 / 30 * 0.25\n",
        "        J = 1\n",
        "        h = -7\n",
        "        all_circs_mirror.extend(\n",
        "            trotter_circuits(G, d_ind_tot, J, h, dt, meas_basis, True)\n",
        "        )\n",
        "circs = all_circs_mirror"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a4b0476d",
      "metadata": {},
      "source": [
        "<span id=\"step-3-execute-using-qiskit-primitives\" />\n",
        "\n",
        "## ステップ3: `Qiskit primitives`を使用して実行する\n",
        "\n",
        "<span id=\"run-mps-simulation\" />\n",
        "\n",
        "### MPSシミュレーションを実行する\n",
        "\n",
        "トロッター化された回路のリストは、 `matrix_product_state` シミュレーターを使い、 $4096$ ショットを任意に選択して実行される。 MPS法は回路ダイナミクスの効率的な近似を提供し、その精度は選択された結合次元によって決定される。 ここで考慮した系サイズでは、デフォルトの結合次元で十分であり、磁化ダイナミクスを忠実に捉えることができる。 生のカウントは正規化され、そこから各トロッター・ステップにおける1量子ビット期待値（ $\\langle Z_i \\rangle$ ）を計算する。 最後に、すべての量子ビットの平均を計算し、磁化が時間とともにどのように変化するかを示す1つの曲線を得る。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "id": "3709531f",
      "metadata": {},
      "outputs": [],
      "source": [
        "backend_sim = AerSimulator(method=\"matrix_product_state\")\n",
        "\n",
        "\n",
        "def normalize_counts(counts_list, shots):\n",
        "    new_counts_list = []\n",
        "    for counts in counts_list:\n",
        "        a = {k: v / shots for k, v in counts.items()}\n",
        "        new_counts_list.append(a)\n",
        "    return new_counts_list\n",
        "\n",
        "\n",
        "def run_sim(circ_list):\n",
        "    shots = 4096\n",
        "    res = backend_sim.run(circ_list, shots=shots)\n",
        "    normed = normalize_counts(res.result().get_counts(), shots)\n",
        "    return normed\n",
        "\n",
        "\n",
        "sim_counts = run_sim(circs)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5895aa46",
      "metadata": {},
      "source": [
        "<span id=\"run-on-hardware\" />\n",
        "\n",
        "### ハードウェア上で実行する\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "id": "4b5a2f87",
      "metadata": {},
      "outputs": [],
      "source": [
        "service = QiskitRuntimeService()\n",
        "backend = service.backend(\"ibm_marrakesh\")\n",
        "\n",
        "\n",
        "def run_qiskit(circ_list):\n",
        "    shots = 4096\n",
        "    pm = generate_preset_pass_manager(backend=backend)\n",
        "    isa_circuits = [pm.run(qc) for qc in circ_list]\n",
        "    sampler = Sampler(mode=backend)\n",
        "    res = sampler.run(isa_circuits, shots=shots)\n",
        "    res = [r.data.c.get_counts() for r in res.result()]\n",
        "    normed = normalize_counts(res, shots)\n",
        "    return normed\n",
        "\n",
        "\n",
        "qiskit_counts = run_qiskit(circs)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d30269ee",
      "metadata": {},
      "source": [
        "<span id=\"run-on-hardware-with-fire-opal\" />\n",
        "\n",
        "### Fire Opal搭載ハードウェアで動作する\n",
        "\n",
        "我々は、実際の量子ハードウェア上で磁化のダイナミクスを評価する。 Fire Opalは、標準の IBM Quantum 推定器プリミティブを拡張し、自動エラー抑制とパフォーマンス管理機能を備えたQiskit関数を提供します。 Trotter化された回路を IBM® のバックエンドに直接送信し、Fire Opalがノイズを考慮した実行を処理します。\n",
        "\n",
        "`pubs`各項目には回路と対応するPauli-Z観測変数が含まれている。 これらはFire Opalのestimator関数に渡され、各トロッターステップにおける各クォビットの期待値 $\\langle Z_i \\rangle$ を返す。 その結果を各クビットで平均することで、ハードウェアから磁化曲線を得ることができる。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "ca8348c9",
      "metadata": {},
      "outputs": [],
      "source": [
        "backend_name = \"ibm_marrakesh\"\n",
        "estimator_pubs = [(qc, observables) for qc in all_circs_mirror[:]]\n",
        "\n",
        "# Run the circuit using the estimator\n",
        "qctrl_estimator_job = perf_mgmt.run(\n",
        "    primitive=\"estimator\",\n",
        "    pubs=estimator_pubs,\n",
        "    backend_name=backend_name,\n",
        "    options={\"default_shots\": 4096},\n",
        ")\n",
        "\n",
        "result_qctrl = qctrl_estimator_job.result()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "aa081308",
      "metadata": {},
      "source": [
        "<span id=\"step-4-post-process-and-return-result-in-desired-classical-format\" />\n",
        "\n",
        "## ステップ4：後処理を行い、結果を希望の古典形式で返す\n",
        "\n",
        "最後に、シミュレーターで得られた磁化曲線を実際のハードウェアで得られた結果と比較する。 両者を並べてプロットすると、Fire Opalを使ったハードウェアの実行が、トロッターステップ全体にわたってノイズのないベースラインといかに密接に一致しているかがわかる。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 102,
      "id": "91dd23a5",
      "metadata": {},
      "outputs": [],
      "source": [
        "def make_correlators(test_counts, nq, d_ind_tot):\n",
        "    mz = np.empty((nq, d_ind_tot))\n",
        "    for d_ind in range(d_ind_tot):\n",
        "        counts = test_counts[d_ind]\n",
        "        for i in range(nq):\n",
        "            mz[i, d_ind] = z_expectation(counts, i)\n",
        "    average_z = np.mean(mz, axis=0)\n",
        "    return np.concatenate((np.array([1]), average_z), axis=0)\n",
        "\n",
        "\n",
        "sim_exp = make_correlators(sim_counts[0:22], nq=nq, d_ind_tot=22)\n",
        "qiskit_exp = make_correlators(qiskit_counts[0:22], nq=nq, d_ind_tot=22)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 103,
      "id": "1581f9f0",
      "metadata": {},
      "outputs": [],
      "source": [
        "qctrl_exp = [ev.data.evs for ev in result_qctrl[:]]\n",
        "qctrl_exp_mean = np.concatenate(\n",
        "    (np.array([1]), np.mean(qctrl_exp, axis=1)), axis=0\n",
        ")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "id": "d1f90811",
      "metadata": {},
      "outputs": [],
      "source": [
        "def make_expectations_plot(\n",
        "    sim_z,\n",
        "    depths,\n",
        "    exp_qctrl=None,\n",
        "    exp_qctrl_error=None,\n",
        "    exp_qiskit=None,\n",
        "    exp_qiskit_error=None,\n",
        "    plot_from=0,\n",
        "    plot_upto=23,\n",
        "):\n",
        "    import numpy as np\n",
        "    import matplotlib.pyplot as plt\n",
        "\n",
        "    depth_ticks = [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22]\n",
        "\n",
        "    d = np.asarray(depths)[plot_from:plot_upto]\n",
        "    sim = np.asarray(sim_z)[plot_from:plot_upto]\n",
        "\n",
        "    qk = (\n",
        "        None\n",
        "        if exp_qiskit is None\n",
        "        else np.asarray(exp_qiskit)[plot_from:plot_upto]\n",
        "    )\n",
        "    qc = (\n",
        "        None\n",
        "        if exp_qctrl is None\n",
        "        else np.asarray(exp_qctrl)[plot_from:plot_upto]\n",
        "    )\n",
        "\n",
        "    qk_err = (\n",
        "        None\n",
        "        if exp_qiskit_error is None\n",
        "        else np.asarray(exp_qiskit_error)[plot_from:plot_upto]\n",
        "    )\n",
        "    qc_err = (\n",
        "        None\n",
        "        if exp_qctrl_error is None\n",
        "        else np.asarray(exp_qctrl_error)[plot_from:plot_upto]\n",
        "    )\n",
        "\n",
        "    # ---- helper(s) ----\n",
        "    def rmse(a, b):\n",
        "        if a is None or b is None:\n",
        "            return None\n",
        "        a = np.asarray(a, dtype=float)\n",
        "        b = np.asarray(b, dtype=float)\n",
        "        mask = np.isfinite(a) & np.isfinite(b)\n",
        "        if not np.any(mask):\n",
        "            return None\n",
        "        diff = a[mask] - b[mask]\n",
        "        return float(np.sqrt(np.mean(diff**2)))\n",
        "\n",
        "    def plot_panel(ax, method_y, method_err, color, label, band_color=None):\n",
        "        # Noiseless reference\n",
        "        ax.plot(d, sim, color=\"grey\", label=\"Noiseless simulation\")\n",
        "\n",
        "        # Method line + band\n",
        "        if method_y is not None:\n",
        "            ax.plot(d, method_y, color=color, label=label)\n",
        "            if method_err is not None:\n",
        "                lo = np.clip(method_y - method_err, -1.05, 1.05)\n",
        "                hi = np.clip(method_y + method_err, -1.05, 1.05)\n",
        "                ax.fill_between(\n",
        "                    d,\n",
        "                    lo,\n",
        "                    hi,\n",
        "                    alpha=0.18,\n",
        "                    color=band_color if band_color else color,\n",
        "                    label=f\"{label} ± error\",\n",
        "                )\n",
        "        else:\n",
        "            ax.text(\n",
        "                0.5,\n",
        "                0.5,\n",
        "                \"No data\",\n",
        "                transform=ax.transAxes,\n",
        "                ha=\"center\",\n",
        "                va=\"center\",\n",
        "                fontsize=10,\n",
        "                color=\"0.4\",\n",
        "            )\n",
        "\n",
        "        # RMSE box (vs sim)\n",
        "        r = rmse(method_y, sim)\n",
        "        if r is not None:\n",
        "            ax.text(\n",
        "                0.98,\n",
        "                0.02,\n",
        "                f\"RMSE: {r:.4f}\",\n",
        "                transform=ax.transAxes,\n",
        "                va=\"bottom\",\n",
        "                ha=\"right\",\n",
        "                fontsize=8,\n",
        "                bbox=dict(\n",
        "                    boxstyle=\"round,pad=0.35\", fc=\"white\", ec=\"0.7\", alpha=0.9\n",
        "                ),\n",
        "            )\n",
        "        # Axes\n",
        "        ax.set_xticks(depth_ticks)\n",
        "        ax.set_ylim(-1.05, 1.05)\n",
        "        ax.grid(True, which=\"both\", linewidth=0.4, alpha=0.4)\n",
        "        ax.set_axisbelow(True)\n",
        "        ax.legend(prop={\"size\": 8}, loc=\"best\")\n",
        "\n",
        "    fig, axes = plt.subplots(1, 2, figsize=(10, 4), dpi=300, sharey=True)\n",
        "\n",
        "    axes[0].set_title(\"Fire Opal (Q-CTRL)\", fontsize=10)\n",
        "    plot_panel(\n",
        "        axes[0],\n",
        "        qc,\n",
        "        qc_err,\n",
        "        color=\"#680CE9\",\n",
        "        label=\"Fire Opal\",\n",
        "        band_color=\"#680CE9\",\n",
        "    )\n",
        "    axes[0].set_xlabel(\"Trotter step\")\n",
        "    axes[0].set_ylabel(r\"$\\langle Z \\rangle$\")\n",
        "    axes[1].set_title(\"Qiskit\", fontsize=10)\n",
        "    plot_panel(\n",
        "        axes[1], qk, qk_err, color=\"blue\", label=\"Qiskit\", band_color=\"blue\"\n",
        "    )\n",
        "    axes[1].set_xlabel(\"Trotter step\")\n",
        "\n",
        "    plt.tight_layout()\n",
        "    plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "id": "415760ad",
      "metadata": {},
      "outputs": [],
      "source": [
        "depths = list(range(d_ind_tot + 1))\n",
        "errors = np.abs(np.array(qctrl_exp_mean) - np.array(sim_exp))\n",
        "\n",
        "errors_qiskit = np.abs(np.array(qiskit_exp) - np.array(sim_exp))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 28,
      "id": "d4902d14",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/transverse-field-ising-model/extracted-outputs/d4902d14-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "make_expectations_plot(\n",
        "    sim_exp,\n",
        "    depths,\n",
        "    exp_qctrl=qctrl_exp_mean,\n",
        "    exp_qctrl_error=errors,\n",
        "    exp_qiskit=qiskit_exp,\n",
        "    exp_qiskit_error=errors_qiskit,\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ee41a301",
      "metadata": {},
      "source": [
        "<span id=\"references\" />\n",
        "\n",
        "## 参照\n",
        "\n",
        "\\[1] グラフの色付け。 ウィキペディア 2025 年 9 月 15 日、 [https://en.wikipedia.org/wiki/Graph](https://en.wikipedia.org/wiki/Graph_coloring) より取得\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2cb5785c",
      "metadata": {},
      "source": [
        "<span id=\"tutorial-survey\" />\n",
        "\n",
        "## チュートリアル調査\n",
        "\n",
        "このチュートリアルに関するフィードバックをお寄せください。 あなたの洞察は、私たちのコンテンツの提供とユーザーエクスペリエンスを向上させるのに役立ちます。\n",
        "\n",
        "[アンケートへのリンク](https://your.feedback.ibm.com/jfe/form/SV_3BLFkNVEuh0QBWm)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "id": "a1b8767d",
      "source": "© IBM Corp., 2017-2026"
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3"
    },
    "hours": 1,
    "qpuSeconds": 120
  },
  "nbformat": 4,
  "nbformat_minor": 5
}