diff --git a/docs/tutorials/_toc.json b/docs/tutorials/_toc.json
index 117c050b8d4..b91a3d0999f 100644
--- a/docs/tutorials/_toc.json
+++ b/docs/tutorials/_toc.json
@@ -145,7 +145,11 @@
{
"title": "Transverse-Field Ising Model Simulation",
"url": "/docs/tutorials/transverse-field-ising-model"
- }
+ },
+ {
+ "title": "Simulate 2D tilted-field Ising with QESEM Qiskit Function",
+ "url": "/docs/tutorials/qedma-2d-ising-with-qesem"
+ }
]
},
{
diff --git a/docs/tutorials/index.mdx b/docs/tutorials/index.mdx
index 165cadedbcd..c37cf8e5ebb 100644
--- a/docs/tutorials/index.mdx
+++ b/docs/tutorials/index.mdx
@@ -104,6 +104,8 @@ Qiskit Functions are a collection of pre-packaged error management and applicati
* [Transverse-Field Ising Model Simulation](/docs/tutorials/transverse-field-ising-model)
+ * [Simulate 2D tilted-field Ising with QESEM Qiskit Function](/docs/tutorials/qedma-2d-ising-with-qesem)
+
- Experiment with domain-specific problems with **Application functions** -- with familiar inputs and outputs to classical solvers.
* [Quantum Portfolio Optimizer - A Qiskit Function by Global Data Quantum](/docs/tutorials/global-data-quantum-optimizer)
diff --git a/docs/tutorials/qedma-2d-ising-with-qesem.ipynb b/docs/tutorials/qedma-2d-ising-with-qesem.ipynb
new file mode 100644
index 00000000000..f17f401d0bb
--- /dev/null
+++ b/docs/tutorials/qedma-2d-ising-with-qesem.ipynb
@@ -0,0 +1,1200 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "aff344db",
+ "metadata": {},
+ "source": [
+ "# Simulate 2D tilted-field Ising with QESEM Qiskit Function"
+ ]
+ },
+ {
+ "metadata": {},
+ "cell_type": "markdown",
+ "source": "*Usage estimate: 20 minutes on a Heron r2 processor. (NOTE: This is an estimate only. Your runtime may vary.)*",
+ "id": "edb68a2a9b5ee911"
+ },
+ {
+ "cell_type": "markdown",
+ "id": "88c617fe",
+ "metadata": {},
+ "source": [
+ "## Background\n",
+ "\n",
+ "This tutorial shows how to simulate dynamics of the 2D tilted-field Ising model:\n",
+ "\n",
+ "$$\n",
+ "H = J \\sum_{\\langle i,j \\rangle} Z_i Z_j + g_x \\sum_i X_i + g_z \\sum_i Z_i\n",
+ "$$\n",
+ "\n",
+ "with non clifford angles using [QESEM Qedma's qiskit function](https://quantum.cloud.ibm.com/docs/en/guides/qedma-qesem).\n",
+ "\n",
+ "We first use a time estimation feature to estimate the expected QPU runtime for full error mitigation run. Then, we demonstrate the use of [operator backpropagation (OBP)](https://quantum.cloud.ibm.com/docs/en/guides/qiskit-addons-obp) to reduce circuit depth, performing EM for all multiple observables simultanously.\n",
+ "\n",
+ "For more information on QESEM and this model, you can refer to [Reliable high-accuracy error mitigation for utility-scale quantum circuits](https://arxiv.org/abs/2508.10997)."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "4a2ede52",
+ "metadata": {},
+ "source": [
+ "## Requirements\n",
+ "\n",
+ "Install the following Python packages before running the notebook:\n",
+ "\n",
+ "- Qiskit SDK v2.0.0 or later (`pip install qiskit`)\n",
+ "- Qiskit Runtime v0.40.0 or later (`pip install qiskit-ibm-runtime`)\n",
+ "- Qiskit Functions Catalog v0.8.0 or later ( `pip install qiskit-ibm-catalog` )\n",
+ "- Qiskit Operator Back Propagation add on v0.3.0 or later ( `pip install qiskit-addon-obp` )\n",
+ "- Qiskit Utils add on v0.1.1 or later ( `pip install qiskit-addon-utils` )\n",
+ "- Qiskit Aer simulator v0.17.1 or later ( `pip install qiskit-aer` )\n",
+ "- Matplotlib v3.10.3 or later ( `pip install matplotlib` )"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "a675b3b1",
+ "metadata": {},
+ "source": [
+ "## Setup\n",
+ "Let's import the relevant libraries:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "acea2e46",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "%matplotlib inline\n",
+ "\n",
+ "from typing import Sequence\n",
+ "\n",
+ "import matplotlib.pyplot as plt\n",
+ "import numpy as np\n",
+ "\n",
+ "import qiskit\n",
+ "from qiskit_ibm_runtime import EstimatorV2 as Estimator\n",
+ "from qiskit_ibm_catalog import QiskitFunctionsCatalog\n",
+ "from qiskit_aer import AerSimulator\n",
+ "from qiskit_addon_utils.slicing import combine_slices, slice_by_gate_types\n",
+ "from qiskit_addon_obp import backpropagate\n",
+ "from qiskit_addon_obp.utils.simplify import OperatorBudget"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "467f1569",
+ "metadata": {},
+ "source": [
+ "Let's set your [IBM Quantum Platform](https://quantum.cloud.ibm.com/) credentials and load the QESEM function:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "41a53d27",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Paste here your instance and token strings\n",
+ "\n",
+ "instance = \"YOUR_INSTANCE\"\n",
+ "token = \"YOUR_TOKEN\"\n",
+ "channel = \"ibm_quantum_platform\"\n",
+ "\n",
+ "catalog = QiskitFunctionsCatalog(\n",
+ " channel=channel, token=token, instance=instance\n",
+ ")\n",
+ "qesem_function = catalog.load(\"qedma/qesem\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "69eafdcc",
+ "metadata": {},
+ "source": [
+ "## Step 1: Define circuit and observables"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ac9dcff0",
+ "metadata": {},
+ "source": [
+ "We'll start by defining a function that creates the trotter circuit:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "3842021c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def trotter_circuit_from_layers(\n",
+ " steps: int,\n",
+ " theta_x: float,\n",
+ " theta_z: float,\n",
+ " theta_zz: float,\n",
+ " layers: Sequence[Sequence[tuple[int, int]]],\n",
+ " init_state: str | None = None,\n",
+ ") -> qiskit.QuantumCircuit:\n",
+ " \"\"\"\n",
+ " Generates an ising trotter circuit\n",
+ " :param steps: trotter steps\n",
+ " :param theta_x: RX angle\n",
+ " :param theta_z: RZ angle\n",
+ " :param theta_zz: RZZ angle\n",
+ " :param layers: list of layers (can be list of layers in device)\n",
+ " :param init_state: Initial state to prepare. If None, will not prepare any state. If \"+\", will\n",
+ " add Hadamard gates to all qubits.\n",
+ " :return: QuantumCircuit\n",
+ " \"\"\"\n",
+ " qubits = sorted({i for layer in layers for edge in layer for i in edge})\n",
+ " circ = qiskit.QuantumCircuit(max(qubits) + 1)\n",
+ "\n",
+ " if init_state == \"+\":\n",
+ " print(\"init_state = +\")\n",
+ " for q in qubits:\n",
+ " circ.h(q)\n",
+ "\n",
+ " for _ in range(steps):\n",
+ " for q in qubits:\n",
+ " circ.rx(theta_x, q)\n",
+ " circ.rz(theta_z, q)\n",
+ "\n",
+ " for layer in layers:\n",
+ " for edge in layer:\n",
+ " circ.rzz(theta_zz, *edge)\n",
+ "\n",
+ " return circ"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "7d342cd7",
+ "metadata": {},
+ "source": [
+ "We use a hardware-based $R_{ZZ}$ layer mapping taken from the Heron device, from which we cut out the layers according to the number of qubits:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 255,
+ "id": "27402210",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[(16, 23), (24, 25), (17, 27)], [(23, 24), (25, 26), (27, 28)], [(22, 23), (26, 27), (25, 37)]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "LAYERS_HERON_R2 = [\n",
+ " [\n",
+ " (2, 3),\n",
+ " (6, 7),\n",
+ " (10, 11),\n",
+ " (14, 15),\n",
+ " (20, 21),\n",
+ " (16, 23),\n",
+ " (24, 25),\n",
+ " (17, 27),\n",
+ " (28, 29),\n",
+ " (18, 31),\n",
+ " (32, 33),\n",
+ " (19, 35),\n",
+ " (36, 41),\n",
+ " (42, 43),\n",
+ " (37, 45),\n",
+ " (46, 47),\n",
+ " (38, 49),\n",
+ " (50, 51),\n",
+ " (39, 53),\n",
+ " (60, 61),\n",
+ " (56, 63),\n",
+ " (64, 65),\n",
+ " (57, 67),\n",
+ " (68, 69),\n",
+ " (58, 71),\n",
+ " (72, 73),\n",
+ " (59, 75),\n",
+ " (76, 81),\n",
+ " (82, 83),\n",
+ " (77, 85),\n",
+ " (86, 87),\n",
+ " (78, 89),\n",
+ " (90, 91),\n",
+ " (79, 93),\n",
+ " (94, 95),\n",
+ " (100, 101),\n",
+ " (96, 103),\n",
+ " (104, 105),\n",
+ " (97, 107),\n",
+ " (108, 109),\n",
+ " (98, 111),\n",
+ " (112, 113),\n",
+ " (99, 115),\n",
+ " (116, 121),\n",
+ " (122, 123),\n",
+ " (117, 125),\n",
+ " (126, 127),\n",
+ " (118, 129),\n",
+ " (130, 131),\n",
+ " (119, 133),\n",
+ " (134, 135),\n",
+ " (140, 141),\n",
+ " (136, 143),\n",
+ " (144, 145),\n",
+ " (137, 147),\n",
+ " (148, 149),\n",
+ " (138, 151),\n",
+ " (152, 153),\n",
+ " (139, 155),\n",
+ " ],\n",
+ " [\n",
+ " (1, 2),\n",
+ " (3, 4),\n",
+ " (5, 6),\n",
+ " (7, 8),\n",
+ " (9, 10),\n",
+ " (11, 12),\n",
+ " (13, 14),\n",
+ " (21, 22),\n",
+ " (23, 24),\n",
+ " (25, 26),\n",
+ " (27, 28),\n",
+ " (29, 30),\n",
+ " (31, 32),\n",
+ " (33, 34),\n",
+ " (40, 41),\n",
+ " (43, 44),\n",
+ " (45, 46),\n",
+ " (47, 48),\n",
+ " (49, 50),\n",
+ " (51, 52),\n",
+ " (53, 54),\n",
+ " (55, 59),\n",
+ " (61, 62),\n",
+ " (63, 64),\n",
+ " (65, 66),\n",
+ " (67, 68),\n",
+ " (69, 70),\n",
+ " (71, 72),\n",
+ " (73, 74),\n",
+ " (80, 81),\n",
+ " (83, 84),\n",
+ " (85, 86),\n",
+ " (87, 88),\n",
+ " (89, 90),\n",
+ " (91, 92),\n",
+ " (93, 94),\n",
+ " (95, 99),\n",
+ " (101, 102),\n",
+ " (103, 104),\n",
+ " (105, 106),\n",
+ " (107, 108),\n",
+ " (109, 110),\n",
+ " (111, 112),\n",
+ " (113, 114),\n",
+ " (120, 121),\n",
+ " (123, 124),\n",
+ " (125, 126),\n",
+ " (127, 128),\n",
+ " (129, 130),\n",
+ " (131, 132),\n",
+ " (133, 134),\n",
+ " (135, 139),\n",
+ " (141, 142),\n",
+ " (143, 144),\n",
+ " (145, 146),\n",
+ " (147, 148),\n",
+ " (149, 150),\n",
+ " (151, 152),\n",
+ " (153, 154),\n",
+ " ],\n",
+ " [\n",
+ " (3, 16),\n",
+ " (7, 17),\n",
+ " (11, 18),\n",
+ " (22, 23),\n",
+ " (26, 27),\n",
+ " (30, 31),\n",
+ " (34, 35),\n",
+ " (21, 36),\n",
+ " (25, 37),\n",
+ " (29, 38),\n",
+ " (33, 39),\n",
+ " (41, 42),\n",
+ " (44, 45),\n",
+ " (48, 49),\n",
+ " (52, 53),\n",
+ " (43, 56),\n",
+ " (47, 57),\n",
+ " (51, 58),\n",
+ " (62, 63),\n",
+ " (66, 67),\n",
+ " (70, 71),\n",
+ " (74, 75),\n",
+ " (61, 76),\n",
+ " (65, 77),\n",
+ " (69, 78),\n",
+ " (73, 79),\n",
+ " (81, 82),\n",
+ " (84, 85),\n",
+ " (88, 89),\n",
+ " (92, 93),\n",
+ " (83, 96),\n",
+ " (87, 97),\n",
+ " (91, 98),\n",
+ " (102, 103),\n",
+ " (106, 107),\n",
+ " (110, 111),\n",
+ " (114, 115),\n",
+ " (101, 116),\n",
+ " (105, 117),\n",
+ " (109, 118),\n",
+ " (113, 119),\n",
+ " (121, 122),\n",
+ " (124, 125),\n",
+ " (128, 129),\n",
+ " (132, 133),\n",
+ " (123, 136),\n",
+ " (127, 137),\n",
+ " (131, 138),\n",
+ " (142, 143),\n",
+ " (146, 147),\n",
+ " (150, 151),\n",
+ " (154, 155),\n",
+ " (0, 1),\n",
+ " (4, 5),\n",
+ " (8, 9),\n",
+ " (12, 13),\n",
+ " (54, 55),\n",
+ " (15, 19),\n",
+ " ],\n",
+ "]\n",
+ "\n",
+ "subgraphs = {\n",
+ " 10: list(range(22, 29)) + [16, 17, 37],\n",
+ " 21: list(range(3, 12)) + list(range(23, 32)) + [16, 17, 18],\n",
+ " 28: list(range(3, 12))\n",
+ " + list(range(23, 32))\n",
+ " + list(range(45, 50))\n",
+ " + [16, 17, 18, 37, 38],\n",
+ "}\n",
+ "\n",
+ "\n",
+ "n_qubits = 10\n",
+ "\n",
+ "layers = [\n",
+ " [\n",
+ " edge\n",
+ " for edge in layer\n",
+ " if edge[0] in subgraphs[n_qubits] and edge[1] in subgraphs[n_qubits]\n",
+ " ]\n",
+ " for layer in LAYERS_HERON_R2\n",
+ "]\n",
+ "\n",
+ "\n",
+ "print(layers)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "7b6ecaa9",
+ "metadata": {},
+ "source": [
+ "Notice that the connectivity of of the chosen qubit layout is not necessarily linear, and can cover large regions of the Heron device depending on the selected number of qubits."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1d34dd0d",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Circuit 2q layers: 15\n",
+ "\n",
+ "Circuit structure:\n"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "SparsePauliOp(['IIIIIIIIIIIIIIIZIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIZIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIZIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIZIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIZIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIZIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIZIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIZIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIZIIIIIIIIIIIIIIIII', 'ZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII'],\n",
+ " coeffs=[0.1+0.j, 0.1+0.j, 0.1+0.j, 0.1+0.j, 0.1+0.j, 0.1+0.j, 0.1+0.j, 0.1+0.j,\n",
+ " 0.1+0.j, 0.1+0.j])\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Chosen parameters for the Hamiltonian terms:\n",
+ "theta_x = 0.53\n",
+ "theta_z = 0.1\n",
+ "theta_zz = 1.0\n",
+ "steps = 5\n",
+ "\n",
+ "circ = trotter_circuit_from_layers(steps, theta_x, theta_z, theta_zz, layers)\n",
+ "print(\n",
+ " f\"Circuit 2q layers: {circ.depth(filter_function=lambda instr: len(instr.qubits) == 2)}\"\n",
+ ")\n",
+ "print(\"\\nCircuit structure:\")\n",
+ "\n",
+ "circ.draw(\"mpl\", scale=0.8, fold=-1, idle_wires=False)\n",
+ "plt.show()\n",
+ "\n",
+ "observable = qiskit.quantum_info.SparsePauliOp.from_sparse_list(\n",
+ " [(\"Z\", [q], 1 / n_qubits) for q in subgraphs[n_qubits]],\n",
+ " np.max(subgraphs[n_qubits]) + 1,\n",
+ ") # Average magnetization observable\n",
+ "\n",
+ "print(observable)\n",
+ "obs_list = [observable]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "76923674",
+ "metadata": {},
+ "source": [
+ "## Step 2: QPU time estimation with and without OBP\n",
+ "Users would typically want to know how much QPU time is required for their experiment.\n",
+ "However, this is considered a hard problem for classical computers.\n",
+ "\n",
+ "QESEM offers two modes of time estimation to inform users about the feasibility of their experiments:\n",
+ "1. Analytical time estimation - gives a very rough estimation and requires no QPU time. This can be used to test if a transpilation pass would potentially reduce the QPU time.\n",
+ "2. Empirical time estimation (demonstrated here) - gives a pretty good estimation and uses a few minutes of QPU time.\n",
+ "\n",
+ "In both cases, QESEM outputs the time estimation for reaching the required precision for all observables."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 287,
+ "id": "478e18ff",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "precision = 0.02\n",
+ "backend_name = \"fake_fez\"\n",
+ "\n",
+ "# Start a job for empirical time estimation\n",
+ "estimation_job_wo_obp = qesem_function.run(\n",
+ " pubs=[(circ, obs_list)],\n",
+ " instance=instance,\n",
+ " backend_name=backend_name, # E.g. \"ibm_brisbane\"\n",
+ " options={\n",
+ " \"estimate_time_only\": \"empirical\", # \"empirical\" - gets actual time estimates without running full mitigation\n",
+ " \"max_execution_time\": 120, # Limits the QPU time, specified in seconds.\n",
+ " \"default_precision\": precision,\n",
+ " },\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "6357b2b5",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "RUNNING\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Get the result object (blocking method). Use job.status() in a loop for non-blocking.\n",
+ "# This takes a 1-3 minutes\n",
+ "result = estimation_job_wo_obp.result()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 289,
+ "id": "1e3eab49",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Empirical time estimation (sec): 600\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(\n",
+ " f\"Empirical time estimation (sec): {result[0].metadata['time_estimation_sec']}\"\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "75dbab74",
+ "metadata": {},
+ "source": [
+ "Now we will use operator backpropagation (OBP), see [OBP](https://quantum.cloud.ibm.com/docs/en/guides/qiskit-addons-obp) for more details on the add-on. Let's generate the circuit slices for backpropagation:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 290,
+ "id": "cbb1d983",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Separated the circuit into 25 slices.\n"
+ ]
+ }
+ ],
+ "source": [
+ "slices = slice_by_gate_types(circ)\n",
+ "print(f\"Separated the circuit into {len(slices)} slices.\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4f85bd72",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Sets a maximal number of measurement groups for OBP\n",
+ "op_budget = OperatorBudget(max_qwc_groups=8)\n",
+ "\n",
+ "# Backpropagate without the truncation error budget\n",
+ "bp_observable, remaining_slices, metadata = backpropagate(\n",
+ " observable,\n",
+ " slices,\n",
+ " operator_budget=op_budget,\n",
+ ")\n",
+ "\n",
+ "# Recombine the slices remaining after backpropagation\n",
+ "bp_circuit = combine_slices(remaining_slices, include_barriers=True)\n",
+ "\n",
+ "print(f\"Backpropagated {metadata.num_backpropagated_slices} slices.\")\n",
+ "print(\n",
+ " f\"New observable has {len(bp_observable.paulis)} terms, which can be combined into \"\n",
+ " f\"{len(bp_observable.group_commuting(qubit_wise=True))} groups.\\n\"\n",
+ " f\"After truncation, the error in our observable is bounded by {metadata.accumulated_error(0):.3e}\"\n",
+ ")\n",
+ "print(\n",
+ " f\"Note that backpropagating one more slice would result in {metadata.backpropagation_history[-1].num_paulis[0]} terms \"\n",
+ " f\"across {metadata.backpropagation_history[-1].num_qwc_groups} groups.\"\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 306,
+ "id": "cedb7fa1",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "The remaining circuit after backpropagation looks as follows:\n"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "print(\"The remaining circuit after backpropagation looks as follows:\")\n",
+ "bp_circuit.draw(\"mpl\", scale=0.8, fold=-1, idle_wires=False)\n",
+ "None"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "fc1e532a",
+ "metadata": {},
+ "source": [
+ "Now that we have our reduced circuit and expanded observables. Let's do time estimation to the backpropagated circuit:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "88160fbc",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Start a job for empirical time estimation\n",
+ "estimation_job_obp = qesem_function.run(\n",
+ " pubs=[(bp_circuit, [bp_observable])],\n",
+ " instance=instance,\n",
+ " backend_name=backend_name,\n",
+ " options={\n",
+ " \"estimate_time_only\": \"empirical\",\n",
+ " \"max_execution_time\": 120,\n",
+ " \"default_precision\": precision,\n",
+ " },\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "19cd4cc2",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "RUNNING\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Get the result object (blocking method). Use job.status() in a loop for non-blocking.\n",
+ "# This takes a 1-3 minutes\n",
+ "result_obp = estimation_job_obp.result()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 312,
+ "id": "feca3059",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Empirical time estimation (sec): 300\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(\n",
+ " f\"Empirical time estimation (sec): {result_obp[0].metadata['time_estimation_sec']}\"\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "504669f5",
+ "metadata": {},
+ "source": [
+ "We see that OBP reduces the time cost for mitigation of the circuit."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "4d1d5092",
+ "metadata": {},
+ "source": [
+ "## Step 3: Run the QESEM function\n",
+ "\n",
+ "### Run with fake backend\n",
+ "With the improved circuit and measurement strategy, we can launch a full QESEM mitigation job:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "e8e10c9a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Start a job for empirical time estimation\n",
+ "full_job = qesem_function.run(\n",
+ " pubs=[(bp_circuit, [bp_observable])],\n",
+ " instance=instance,\n",
+ " backend_name=backend_name,\n",
+ " options={\n",
+ " \"max_execution_time\": 900,\n",
+ " \"default_precision\": 0.05,\n",
+ " },\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "90820fd4",
+ "metadata": {},
+ "source": [
+ "Let's read the results and compare the ideal, noisy, and mitigated estimates."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 148,
+ "id": "4876b6f3",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "------------------------------\n",
+ "Observable: Average Magnetization\n",
+ "Ideal: 0.8559570312500001\n",
+ "Noisy: 0.7899730159551378 ± 0.004321764424558104\n",
+ "QESEM: 0.8457341051157278 ± 0.01266352434070931\n",
+ "------------------------------\n",
+ "Gate fidelities found: {'ID1Q': 0.9990004566092757, 'RZZ': 0.9935311791048811}\n"
+ ]
+ }
+ ],
+ "source": [
+ "result = full_job.result() # Blocking - takes 3-5 minutes\n",
+ "noisy_results = result[0].metadata[\"noisy_results\"]\n",
+ "\n",
+ "\n",
+ "def calculate_ideal_evs(circ, obs):\n",
+ " simulator = AerSimulator()\n",
+ "\n",
+ " # Use Estimator primitive to get expectation value\n",
+ " estimator = Estimator(simulator)\n",
+ " sim_result = estimator.run([(circ, [obs])]).result()\n",
+ "\n",
+ " # Extracting the result\n",
+ " ideal_values = sim_result[0].data.evs[0]\n",
+ " return ideal_values\n",
+ "\n",
+ "\n",
+ "for en, obs in enumerate(obs_list):\n",
+ " print(\"-\" * 30)\n",
+ " print(\"Observable: \" + [\"Average Magnetization\", \"ZZZZ\"][en])\n",
+ " # print (f\"Ideal: {Statevector(circ).expectation_value(obs).real}\")\n",
+ " print(f\"Ideal: {calculate_ideal_evs(circ, obs)}\")\n",
+ " print(f\"Noisy: {noisy_results.evs[en]} \\u00b1 {noisy_results.stds[en]}\")\n",
+ " print(f\"QESEM: {result[0].data.evs[en]} \\u00b1 {result[0].data.stds[en]}\")\n",
+ "\n",
+ "\n",
+ "print(\"-\" * 30)\n",
+ "print(\n",
+ " f\"Gate fidelities found: {result[0].metadata['gate_fidelities']}\"\n",
+ ") # Some of the data gathered during a QESEM run."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "3da535e9",
+ "metadata": {},
+ "source": [
+ "### Run with real backend\n",
+ "Let's move to larger circuits with 21 qubits and repeat the experiments on real quantum hardware. The number of qubits and required precision can be modified according to the available QPU resources.\n",
+ "\n",
+ "We examine 4 different circuits with precision of 0.05, and compare their ideal, noisy and mitigated expectation values:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7cfb4dbc",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "n_qubits = 21 # can be modified to 10 or 28 qubits\n",
+ "\n",
+ "layers = [\n",
+ " [\n",
+ " edge\n",
+ " for edge in layer\n",
+ " if edge[0] in subgraphs[n_qubits] and edge[1] in subgraphs[n_qubits]\n",
+ " ]\n",
+ " for layer in LAYERS_HERON_R2\n",
+ "]\n",
+ "\n",
+ "\n",
+ "observable = qiskit.quantum_info.SparsePauliOp.from_sparse_list(\n",
+ " [(\"Z\", [q], 1 / n_qubits) for q in subgraphs[n_qubits]],\n",
+ " np.max(subgraphs[n_qubits]) + 1,\n",
+ ") # Average magnetization observable\n",
+ "\n",
+ "\n",
+ "steps_vec = [3, 5, 7, 9]\n",
+ "\n",
+ "\n",
+ "circ_vec = []\n",
+ "for steps in steps_vec:\n",
+ " circ = trotter_circuit_from_layers(\n",
+ " steps, theta_x, theta_z, theta_zz, layers\n",
+ " )\n",
+ " circ_vec.append(circ)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d63c1de0",
+ "metadata": {},
+ "source": [
+ "Again, performing OBP on each circuit to reduce runtime:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 314,
+ "id": "267e030f",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "n.o. steps: 3\n",
+ "Backpropagated 9 slices.\n",
+ "New observable has 205 terms, which can be combined into 6 groups.\n",
+ "After truncation, the error in our observable is bounded by 0.000e+00\n",
+ "-----------------\n",
+ "n.o. steps: 5\n",
+ "Backpropagated 9 slices.\n",
+ "New observable has 205 terms, which can be combined into 6 groups.\n",
+ "After truncation, the error in our observable is bounded by 0.000e+00\n",
+ "-----------------\n",
+ "n.o. steps: 7\n",
+ "Backpropagated 9 slices.\n",
+ "New observable has 205 terms, which can be combined into 6 groups.\n",
+ "After truncation, the error in our observable is bounded by 0.000e+00\n",
+ "-----------------\n",
+ "n.o. steps: 9\n",
+ "Backpropagated 9 slices.\n",
+ "New observable has 205 terms, which can be combined into 6 groups.\n",
+ "After truncation, the error in our observable is bounded by 0.000e+00\n",
+ "-----------------\n"
+ ]
+ }
+ ],
+ "source": [
+ "bp_circuit_vec = []\n",
+ "bp_observable_vec = []\n",
+ "\n",
+ "for i, circ in enumerate(circ_vec):\n",
+ " slices = slice_by_gate_types(circ)\n",
+ " bp_observable, remaining_slices, metadata = backpropagate(\n",
+ " observable,\n",
+ " slices,\n",
+ " operator_budget=op_budget,\n",
+ " )\n",
+ " slices\n",
+ " # Recombine the slices remaining after backpropagation\n",
+ " bp_circuit = combine_slices(remaining_slices, include_barriers=True)\n",
+ " bp_circuit_vec.append(bp_circuit)\n",
+ " bp_observable_vec.append(bp_observable)\n",
+ " print(f\"n.o. steps: {steps_vec[i]}\")\n",
+ " print(f\"Backpropagated {metadata.num_backpropagated_slices} slices.\")\n",
+ " print(\n",
+ " f\"New observable has {len(bp_observable.paulis)} terms, which can be combined into \"\n",
+ " f\"{len(bp_observable.group_commuting(qubit_wise=True))} groups.\\n\"\n",
+ " f\"After truncation, the error in our observable is bounded by {metadata.accumulated_error(0):.3e}\"\n",
+ " )\n",
+ " print(\"-----------------\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "2322b7da",
+ "metadata": {},
+ "source": [
+ "We run time estimation on the deepest circuit to gauge execution costs before dispatching the full jobs."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b6a2ef22",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "run_on_real_hardware = True\n",
+ "\n",
+ "precision = 0.05\n",
+ "if run_on_real_hardware:\n",
+ " backend_name = \"ibm_fez\"\n",
+ "else:\n",
+ " backend_name = \"fake_fez\"\n",
+ "\n",
+ "pubs = [(bp_circuit_vec[-1], bp_observable_vec[-1])]\n",
+ "# Start a job for empirical time estimation\n",
+ "estimation_job_obp = qesem_function.run(\n",
+ " pubs=pubs,\n",
+ " instance=instance,\n",
+ " backend_name=backend_name,\n",
+ " options={\n",
+ " \"estimate_time_only\": \"empirical\",\n",
+ " \"max_execution_time\": 120,\n",
+ " \"default_precision\": precision,\n",
+ " },\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 332,
+ "id": "5554d8f6",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "DONE\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(estimation_job_obp.status())\n",
+ "# print(estimation_job_obp.logs())"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 333,
+ "id": "61e0287f",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Empirical time estimation (sec): 600\n"
+ ]
+ }
+ ],
+ "source": [
+ "result_obp = estimation_job_obp.result()\n",
+ "print(\n",
+ " f\"Empirical time estimation (sec): {result_obp[0].metadata['time_estimation_sec']}\"\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "68bb0915",
+ "metadata": {},
+ "source": "Now we run a batch of full QESEM jobs. We limit the maximal QPU runtime for each of the points for better control on the QPU budget."
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "e60a2fc8",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Running full jobs for:\n",
+ "pubs_list = [\n",
+ " [(bp_circuit_vec[i], bp_observable_vec[i])]\n",
+ " for i in range(len(bp_observable_vec))\n",
+ "]\n",
+ "\n",
+ "# Initiating multiple jobs for different lengths\n",
+ "job_list = []\n",
+ "for pubs in pubs_list:\n",
+ " job_obp = qesem_function.run(\n",
+ " pubs=pubs,\n",
+ " instance=instance,\n",
+ " backend_name=backend_name, # E.g. \"ibm_brisbane\"\n",
+ " options={\n",
+ " \"max_execution_time\": 300, # Limits the QPU time, specified in seconds.\n",
+ " \"default_precision\": 0.05,\n",
+ " },\n",
+ " )\n",
+ " job_list.append(job_obp)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "05c75ada",
+ "metadata": {},
+ "source": [
+ "Checking the status of each job:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 331,
+ "id": "b869fd4f",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "DONE\n",
+ "DONE\n",
+ "DONE\n",
+ "DONE\n"
+ ]
+ }
+ ],
+ "source": [
+ "for job in job_list:\n",
+ " print(job.status())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "426ba0f9",
+ "metadata": {},
+ "source": [
+ "When all jobs are finished running, we can compare their noisy and mitigated expectation value"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 335,
+ "id": "4e9721e5",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "---------------------------------\n",
+ "Ideal: 0.8546084449404764\n",
+ "Noisy: 0.7733082441338024\n",
+ "QESEM: 0.8603361050419767 ± 0.0057385951110160115\n",
+ "---------------------------------\n",
+ "Ideal: 0.7799479166666666\n",
+ "Noisy: 0.6743286862085316\n",
+ "QESEM: 0.7970063156559638 ± 0.013084183430351701\n",
+ "---------------------------------\n",
+ "Ideal: 0.742978050595238\n",
+ "Noisy: 0.6207282012644596\n",
+ "QESEM: 0.7478205911968454 ± 0.025089485167607922\n",
+ "---------------------------------\n",
+ "Ideal: 0.7480236235119049\n",
+ "Noisy: 0.5775631714071018\n",
+ "QESEM: 0.7551863824678515 ± 0.052145546823300824\n"
+ ]
+ }
+ ],
+ "source": [
+ "ideal_values = []\n",
+ "noisy_values = []\n",
+ "error_mitigated_values = []\n",
+ "error_mitigated_stds = []\n",
+ "\n",
+ "for i in range(len(job_list)):\n",
+ " job = job_list[i]\n",
+ " result = job.result() # Blocking - takes 3-5 minutes\n",
+ " noisy_results = result[0].metadata[\"noisy_results\"]\n",
+ "\n",
+ " ideal_val = calculate_ideal_evs(circ_vec[i], observable)\n",
+ " print(\"---------------------------------\")\n",
+ " print(f\"Ideal: {ideal_val}\")\n",
+ " print(f\"Noisy: {noisy_results.evs}\")\n",
+ " print(f\"QESEM: {result[0].data.evs} \\u00b1 {result[0].data.stds}\")\n",
+ "\n",
+ " ideal_values.append(ideal_val)\n",
+ " noisy_values.append(noisy_results.evs)\n",
+ " error_mitigated_values.append(result[0].data.evs)\n",
+ " error_mitigated_stds.append(result[0].data.stds)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "68cabcc4",
+ "metadata": {},
+ "source": [
+ "## Step 4: Visualize results\n",
+ "\n",
+ "Lastly, we can plot the magnetization versus the number of steps. This summarizes the benefit of using QESEM Qiskit function for bias-free error mitigation on noisy quantum devices."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 337,
+ "id": "0f1a44d0",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "Text(0, 0.5, 'Magnetization')"
+ ]
+ },
+ "execution_count": 337,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plt.plot(steps_vec, ideal_values, \"--\", label=\"ideal\")\n",
+ "plt.scatter(steps_vec, noisy_values, label=\"noisy\")\n",
+ "plt.errorbar(\n",
+ " steps_vec,\n",
+ " error_mitigated_values,\n",
+ " yerr=error_mitigated_stds,\n",
+ " fmt=\"o\",\n",
+ " capsize=5,\n",
+ " label=\"QESEM mitigation\",\n",
+ ")\n",
+ "plt.legend()\n",
+ "plt.xlabel(\"n.o. steps\")\n",
+ "plt.ylabel(\"Magnetization\")"
+ ]
+ }
+ ],
+ "metadata": {
+ "description": "This tutorial shows a simulation of 2d transverse-field Ising model with QESEM error mitigation combined with the Qiskit operator backpropagation module.",
+ "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"
+ },
+ "title": "Simulate 2D tilted-field Ising with QESEM Qiskit Function"
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/public/docs/images/tutorials/qedma-2d-ising-with-qesem/extracted-outputs/0f1a44d0-1.avif b/public/docs/images/tutorials/qedma-2d-ising-with-qesem/extracted-outputs/0f1a44d0-1.avif
new file mode 100644
index 00000000000..f5c4f31fd2e
Binary files /dev/null and b/public/docs/images/tutorials/qedma-2d-ising-with-qesem/extracted-outputs/0f1a44d0-1.avif differ
diff --git a/public/docs/images/tutorials/qedma-2d-ising-with-qesem/extracted-outputs/1d34dd0d-1.avif b/public/docs/images/tutorials/qedma-2d-ising-with-qesem/extracted-outputs/1d34dd0d-1.avif
new file mode 100644
index 00000000000..4b0d6c1b477
Binary files /dev/null and b/public/docs/images/tutorials/qedma-2d-ising-with-qesem/extracted-outputs/1d34dd0d-1.avif differ
diff --git a/public/docs/images/tutorials/qedma-2d-ising-with-qesem/extracted-outputs/cedb7fa1-1.avif b/public/docs/images/tutorials/qedma-2d-ising-with-qesem/extracted-outputs/cedb7fa1-1.avif
new file mode 100644
index 00000000000..2e377403290
Binary files /dev/null and b/public/docs/images/tutorials/qedma-2d-ising-with-qesem/extracted-outputs/cedb7fa1-1.avif differ
diff --git a/qiskit_bot.yaml b/qiskit_bot.yaml
index e8fe9a2ba00..35f449c30b7 100644
--- a/qiskit_bot.yaml
+++ b/qiskit_bot.yaml
@@ -658,6 +658,10 @@ notifications:
- "@alexshih"
- "@johannesgreiner"
- "@annaliese-estes"
+ "docs/tutorials/qedma-2d-ising-with-qesem":
+ - "@miamico"
+ - "@oria-qedma"
+ - "@assafb"
"docs/tutorials/global-data-quantum-optimizer":
- "@abbycross"
- "@pandasa123"
diff --git a/scripts/config/notebook-testing.toml b/scripts/config/notebook-testing.toml
index 4af30f0fab3..63271994416 100644
--- a/scripts/config/notebook-testing.toml
+++ b/scripts/config/notebook-testing.toml
@@ -162,6 +162,7 @@ notebooks = [
"docs/guides/qiskit-transpiler-service.ipynb",
# We never run tutorials notebooks
+ "docs/tutorials/qedma-2d-ising-with-qesem.ipynb",
"docs/tutorials/transverse-field-ising-model.ipynb",
"docs/tutorials/fractional-gates.ipynb",
"docs/tutorials/shors-algorithm.ipynb",