nepec/mitiq-qiskit-pec-nepec/mitiq-qiskit-zne-pec-nepec.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5e0dfbd2",
   "metadata": {},
   "source": [
    "# Using Mitiq with Qiskit quantum programs\n",
    "\n",
    "This notebook shows how to use Mitiq to mitigate errors in Qiskit quantum programs.\n",
    "\n",
    "The main sections are:\n",
    "* **Section 1: Computing quantum expectation values without error mitigation**\n",
    "* **Section 2: Applying zero-noise extrapolation (ZNE)**\n",
    "* **Section 3: Applying probabilistic error cancellation (PEC)**\n",
    "\n",
    "Useful links:\n",
    "* Qiskit repository: https://github.com/Qiskit/qiskit\n",
    "* Qiskit documentation: https://qiskit.org/documentation/\n",
    "* IBMQ portal: https://www.ibm.com/quantum-computing/\n",
    "\n",
    "\n",
    "* Mitiq repository: https://github.com/unitaryfund/mitiq\n",
    "* Mitiq documentation: https://mitiq.readthedocs.io/en/stable/\n",
    "* Mitiq white paper: https://mitiq.readthedocs.io/en/stable/\n",
    "* Unitary Fund:  https://unitary.fund"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6575702b",
   "metadata": {},
   "source": [
    "This notebook requires the Python packages: `qiskit` and `mitiq`. If necessary, uncomment the next code cell to install them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1893c653",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install qiskit --quiet\n",
    "!pip install mitiq --quiet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c43ff597",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "warnings.filterwarnings(action='ignore') # Optional warning filter"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fc8ad7a",
   "metadata": {},
   "source": [
    "## Section 1: Computing quantum expectation values without error mitigation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a96f2bc",
   "metadata": {},
   "source": [
    "### Define the circuit of interest"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b89b21a",
   "metadata": {},
   "source": [
    "We first define an arbitrary circuit of interest. For simplicity, in this example we use\n",
    "a two-qubit randomized-benchmarking circuit whose ideal execution is equivalent to the\n",
    "identity operation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0c6c7565",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Randomized benchmarking circuit with 211 gates generated.\n"
     ]
    }
   ],
   "source": [
    "import qiskit\n",
    "import qiskit.ignis.verification.randomized_benchmarking as rb\n",
    "from qiskit.transpiler.passes import RemoveBarriers, RemoveFinalMeasurements\n",
    "from mitiq.interface import convert_to_mitiq, convert_from_mitiq\n",
    "\n",
    "n_qubits = 2\n",
    "depth = 20\n",
    "\n",
    "qiskit_circuit = rb.randomized_benchmarking_seq(\n",
    "    length_vector=[depth], rb_pattern=[range(n_qubits)], group_gates=\"0\", rand_seed=1,\n",
    ")[0][0][0]\n",
    "\n",
    "# Remove Barriers\n",
    "qiskit_circuit=RemoveBarriers()(qiskit_circuit)\n",
    "# Remove Measurements\n",
    "qiskit_circuit=RemoveFinalMeasurements()(qiskit_circuit)\n",
    "# Double conversion trick to rename qubit register\n",
    "circuit = convert_from_mitiq(convert_to_mitiq(qiskit_circuit)[0], \"qiskit\")\n",
    "\n",
    "print(f\"Randomized benchmarking circuit with {len(circuit)} gates generated.\")\n",
    "# print(circuit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c72a789",
   "metadata": {},
   "source": [
    "### Execute the circuit with an noiseless backend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "7a42c100",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Set the number of shots \n",
    "shots = 10 ** 5\n",
    "\n",
    "# Initialize ideal backend (classical noiseless simulator)\n",
    "ideal_backend = qiskit.Aer.get_backend('aer_simulator')\n",
    "\n",
    "# Append measurements\n",
    "circuit_to_run = circuit.copy()\n",
    "circuit_to_run.measure_all()\n",
    "\n",
    "# Run and get count\n",
    "job = ideal_backend.run(circuit_to_run, shots=shots)\n",
    "counts = job.result().get_counts()\n",
    "\n",
    "# Expectation value of A=|00><00|\n",
    "ideal_value = counts[\"00\"] / shots\n",
    "ideal_value"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e57a055",
   "metadata": {},
   "source": [
    "### Execute the circuit with an noisy backend (without error mitigation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8aa79887",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Executing circuit with 211 gates using 100000 shots.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.75779"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from qiskit.test.mock import FakeLima\n",
    "from qiskit.providers.aer.noise import NoiseModel\n",
    "\n",
    "# Set the number of shots \n",
    "shots = 10 ** 5\n",
    "\n",
    "# Select a noisy backend\n",
    "# noisy_backend = qiskit.IBMQ.load_account().get_backend(\"ibmq_lima\")\n",
    "noisy_backend = FakeLima() # Simulator with noise model similar to \"ibmq_lima\"\n",
    "\n",
    "# Append measurements\n",
    "circuit_to_run = circuit.copy()\n",
    "circuit_to_run.measure_all()\n",
    "\n",
    "# Run and get counts\n",
    "print(f\"Executing circuit with {len(circuit)} gates using {shots} shots.\")\n",
    "job = noisy_backend.run(circuit_to_run, shots=shots)\n",
    "counts = job.result().get_counts()\n",
    "\n",
    "# Compute expectation value of the observable A=|0><0|\n",
    "noisy_value = counts[\"00\"] / shots\n",
    "noisy_value"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8dbc201e",
   "metadata": {},
   "source": [
    "## Section 2: Applying zero-noise extrapolation (ZNE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d23b58e",
   "metadata": {},
   "source": [
    "### Define an executor function"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7ce55e6",
   "metadata": {},
   "source": [
    "We rewrite the last code cell of the previous section as a function that we call `executor`.\n",
    "This function takes as input a Qiskit circuit and returns the final noisy expectation value.\n",
    "\n",
    "Hint: It is usually better, but not necessary, to set `optimization_level=0` to avoid automatic circuit optimizations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "72273e28",
   "metadata": {},
   "outputs": [],
   "source": [
    "def executor(circuit, shots=10 ** 5):\n",
    "    \"\"\"Executes the input circuit and returns the noisy expectation value <A>, where A=|00>00|.\n",
    "    \"\"\"\n",
    "    # Select a noisy backend\n",
    "    # noisy_backend = qiskit.IBMQ.load_account().get_backend(\"ibmq_lima\")\n",
    "    noisy_backend = FakeLima() # Simulator with noise model similar to \"ibmq_lima\"\n",
    "\n",
    "    # Append measurements\n",
    "    circuit_to_run = circuit.copy()\n",
    "    circuit_to_run.measure_all()\n",
    "\n",
    "    # Run and get counts\n",
    "    print(f\"Executing circuit with {len(circuit)} gates using {shots} shots.\")\n",
    "    job = noisy_backend.run(circuit_to_run, shots=shots)\n",
    "    counts = job.result().get_counts()\n",
    "\n",
    "    # Compute expectation value of the observable A=|0><0|\n",
    "    noisy_value = counts[\"00\"] / shots\n",
    "    return noisy_value  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "80a0d5c7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Executing circuit with 211 gates using 100000 shots.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.75789"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "noisy_value = executor(circuit)\n",
    "noisy_value"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4f78900",
   "metadata": {},
   "source": [
    "### Run zero-noise extrapolation with Mitiq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "5042f5d9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Executing circuit with 211 gates using 100000 shots.\n",
      "Executing circuit with 421 gates using 100000 shots.\n",
      "Executing circuit with 633 gates using 100000 shots.\n",
      "Estimation error without Mitiq   : 0.24211000000000005\n",
      "Estimation error with Mitiq (ZNE): 0.05457999999999974\n"
     ]
    }
   ],
   "source": [
    "from mitiq import zne\n",
    "\n",
    "zne_value = zne.execute_with_zne(circuit, executor)\n",
    "\n",
    "unmitigated_error = abs(ideal_value - noisy_value)\n",
    "zne_error = abs(ideal_value - zne_value)\n",
    "print(f\"Estimation error without Mitiq   : {unmitigated_error}\")\n",
    "print(f\"Estimation error with Mitiq (ZNE): {zne_error}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "412229f4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ZNE reduced the estimation error by 77.5%.\n"
     ]
    }
   ],
   "source": [
    "print(f\"ZNE reduced the estimation error by {(unmitigated_error - zne_error) / unmitigated_error :.1%}.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d2ea58c",
   "metadata": {},
   "source": [
    "Below we show how to select different options for zero-noise extrapolation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23c14d22",
   "metadata": {},
   "source": [
    "### Select a noise scaling method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "fea50ad8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choise a unitary folding function\n",
    "noise_scaling_function = zne.scaling.fold_global"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae8f2f67",
   "metadata": {},
   "source": [
    "Let's check the effect of the noise scaling function on the input circuit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "c6f22aa5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The input circuit has 211 gates\n",
      "The scaled circuit has 423 gates\n"
     ]
    }
   ],
   "source": [
    "noise_scaled_circuit = noise_scaling_function(circuit, scale_factor=2)\n",
    "\n",
    "print(f\"The input circuit has {len(circuit)} gates\")\n",
    "print(f\"The scaled circuit has {len(noise_scaled_circuit)} gates\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9807fe86",
   "metadata": {},
   "source": [
    "### Select an extrapolation model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "2ff3b526",
   "metadata": {},
   "outputs": [],
   "source": [
    "factory = zne.inference.AdaExpFactory(steps=10, asymptote=0.25)\n",
    "\n",
    "# Examples:\n",
    "\n",
    "# Richardson extrapolation\n",
    "# factory = zne.inference.RichardsonFactory(scale_factors = [1, 2, 3])\n",
    "\n",
    "# Exponential extrapolation\n",
    "# factory = zne.inference.ExpFactory(scale_factors = [1, 2, 3], asymptote=0.25)\n",
    "\n",
    "# Exponential adaptive extrapolation\n",
    "# factory = zne.inference.AdaExpFactory(steps=10, asymptote=0.25)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e64880a",
   "metadata": {},
   "source": [
    "### Apply ZNE with non-default options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "a35eb598",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Executing circuit with 211 gates using 100000 shots.\n",
      "Executing circuit with 423 gates using 100000 shots.\n",
      "Executing circuit with 943 gates using 100000 shots.\n",
      "Executing circuit with 981 gates using 100000 shots.\n",
      "Executing circuit with 983 gates using 100000 shots.\n",
      "Executing circuit with 987 gates using 100000 shots.\n",
      "Executing circuit with 989 gates using 100000 shots.\n",
      "Executing circuit with 987 gates using 100000 shots.\n",
      "Executing circuit with 989 gates using 100000 shots.\n",
      "Executing circuit with 989 gates using 100000 shots.\n",
      "Estimation error without Mitiq   : 0.24211000000000005\n",
      "Estimation error with Mitiq (ZNE): 0.03939425167038613\n"
     ]
    }
   ],
   "source": [
    "zne_value = zne.execute_with_zne(circuit, executor, scale_noise=noise_scaling_function, factory=factory)\n",
    "\n",
    "unmitigated_error = abs(ideal_value - noisy_value)\n",
    "zne_error = abs(ideal_value - zne_value)\n",
    "print(f\"Estimation error without Mitiq   : {unmitigated_error}\")\n",
    "print(f\"Estimation error with Mitiq (ZNE): {zne_error}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "f539d7a1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ZNE reduced the estimation error by 83.7%.\n"
     ]
    }
   ],
   "source": [
    "print(f\"ZNE reduced the estimation error by {(unmitigated_error - zne_error) / unmitigated_error :.1%}.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7212d423",
   "metadata": {},
   "source": [
    "### Analize and visualize ZNE data "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "0d375386",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbkAAAE9CAYAAACWdRzmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA200lEQVR4nO3dd3xUZb7H8c+PhFBDCEUpQYooiFJVQIqgq9JRKRqa5YrIghe5UgTXgqhXFkTFAiu4KgsiRZAmglwwgDTpILhgBEuyi4DSAlKSPPePDNmAASaEyclMvu/X6yznnDkz851nefnjnPOc5zHnHCIiIqEon9cBREREAkVFTkREQpaKnIiIhCwVORERCVkqciIiErJU5EREJGSFex0gq4oXL+6qVq3qdYygdezYMYoUKeJ1jKCl9ssetV/2qQ0zt2HDhgPOudLn7g+6InfllVeyfv16r2MErbi4OJo3b+51jKCl9ssetV/2qQ0zZ2Y/ZrZflytFRCRkqciJiEjIUpETEZGQpSInIiIhS0VORERCloqciIiELBU5EREJWSpyIiISslTkREQkZAVdkUtJSUGzmYuIiD+Crsj9+OOPjBo1yusYIiISBIKuyCUnJzN06FCWLFnidRQREcnlgq7IlShRgtTUVO6//35+/DHT8ThFRESAICxypUqVokWLFvz666907NiREydOeB1JRERyqaArcgBTpkyhUqVKbNiwgb59+6ojioiIZCooi1yJEiWYNWsWBQsW5P3332fjxo1eRxIRkVwo6CZNPaNu3bq8//77REZGcuONN3odR0REcqGgLXIAXbp08TqCiIjkYkF5uTIzK1asoGvXriQnJ3sdRUREcomgPpM74+TJk3Tt2pWEhARiYmIYOXKk15FERCQXCIkzuQIFCjBlyhTCw8MZNWoUM2bM8DqSiIjkAiFR5ACaNm3Kq6++CsDDDz/Mjh07PE4kIiJeC5kiB9CvXz+6du3KsWPHuPfeezl8+LDXkURExEMhVeTMjPHjx1OzZk127drFf/3Xf3kdSUREPBRSRQ6gSJEizJo1iypVqtCrVy+v44iIiIdConfluapWrcrOnTsJDw/JnyciIn4KuTO5MzIWuLi4OH744QfvwoiIiCdCtsid8emnn3LHHXfQoUMHfv/9d6/jiIhIDgr5Ite8eXMqVqzIpk2b6NOnj2YsEBHJQ0K+yEVHRzNr1iwKFSrEhx9+yLvvvut1JBERySEhX+QAateuzYQJE4C0Z+lWr17tcSIREckJeaLIAXTr1o1+/fpx+vRpOnXqxC+//OJ1JBERCbA8U+QAXn31VZo0acLJkyfV21JEJA8I6INkZtYSGAOEAe8550ac83pF4H2gNPAb0N05lxCoPPnz52fGjBmcPHmSihUrBuprREQklwjYmZyZhQHvAK2AGkAXM6txzmGvAv9wztUChgOvBCrPGWXKlDmrwB04cCDQXykiIh4J5OXK+kC8c263c+4UMBW4+5xjagBLfetfZvJ6wDjnGDFiBJUrV+brr7/Oqa8VEZEcFMjLleWBnzNsJwANzjlmC9CBtEua9wKRZlbSOfdrxoPMrBfQC6B06dLExcVlO5xzjhUrVpCUlETLli0ZO3YsZcqUyfbn5nZJSUmXpf3yKrVf9qj9sk9tmDUWqIejzawT0NI519O33QNo4Jx7PMMx5YC3gcrAcqAjcINz7tD5PrdatWpu586dlyXj6dOnadWqFUuWLKFGjRqsWrWKqKioy/LZuVVcXBzNmzf3OkbQUvtlj9ov+9SGmTOzDc65m87dH8jLlYlAhQzbMb596Zxz/3LOdXDO1QX+4tt3KICZzpI/f34++eQTatSowY4dO+jcuTOnT5/Oqa8XEZEAC2SRWwdcY2aVzSwCiAXmZjzAzEqZ2ZkMQ0nraZmjihcvzvz587niiitYvHgxffv21dBfIiIhImBFzjmXDDwOLAK+BaY757ab2XAza+87rDmw08x2AVcCLwcqz4VUrlyZuXPnUrBgQebPn68HxUVEQkRAn5Nzzi0AFpyz77kM658AnwQyg78aNGjAp59+yvXXX58nOqCIiOQFmlU0g5YtW561ffDgQaKjoz1KIyIi2ZWnhvXyl3OOF154geuuu07Df4mIBDEVuUykpKTw1Vdf8csvv9CmTRsOHz7sdSQREbkEKnKZCA8PZ8aMGXq0QEQkyKnInceZRwtKly6tRwtERIKUitwFZHy0YMKECbz66qteRxIRkSxQkbuIhg0bMmnSJADGjh3LsWPHPE4kIiL+UpHzQ6dOnZg4cSKrV6+mSJEiXscRERE/6Tk5Pz3wwANnbR87dkwFT0Qkl9OZXBalpqYyZMgQGjRooEcLRERyORW5LDp27Bhz585l+/bterRARCSXU5HLosjISD777DM9WiAiEgRU5C6BHi0QEQkOKnKXKOOjBYMHD2bmzJkeJxIRkXOpyGVDp06dGDFiBAAvvPACKSkpHicSEZGM9AhBNg0ePJjw8HAefPBBwsLCvI4jIiIZqMhlk5kxYMCA9G3nHElJSURGRnqYSkREQJcrL6vU1FT69evHrbfeyqFDh7yOIyKS56nIXUaHDh1i0aJFbN68mXbt2nH8+HGvI4mI5GkqcpdRiRIlWLx4MTExMXz11Vd06tSJU6dOeR1LRCTPUpG7zCpWrMjixYspVaoUn3/+OQ888IB6XYqIeERFLgCqV6/OwoULiYyMZNq0afTp00ejooiIeEBFLkBuvPFG5s+fT8GCBVm9ejVHjhzxOpKISJ6jRwgC6NZbb2XhwoXUqlWLqKgor+OIiOQ5OpMLsGbNmhEdHQ2kPUO3Zs0ajxOJiOQdKnI5xDlHr169aNSoEVOnTvU6johInqAil0PMjCpVquCco0ePHixYsMDrSCIiIU9FLgcNGTKEQYMGkZycTMeOHVmxYoXXkUREQpqKXA4yM/7617/y6KOPcuLECdq2bcvGjRu9jiUiErJU5HKYmTFu3Djuu+8+jhw5QosWLdi5c6fXsUREQpIeIfBAWFgYkyZN4siRI+zfv5+SJUt6HUlEJCSpyHkkIiKCmTNnkpycTLFixbyOIyISknS50kOFCxdOL3CpqamMHDlSU/SIiFxGKnK5xDPPPMNTTz2lKXpERC6jgBY5M2tpZjvNLN7MhmTy+lVm9qWZbTKzrWbWOpB5crPHHnssfYqejh07aooeEZHLIGBFzszCgHeAVkANoIuZ1TjnsGeA6c65ukAsMDZQeXK7jFP0LFy4kB49emiKHhGRbArkmVx9IN45t9s5dwqYCtx9zjEOONPrIgr4VwDz5HrVq1dn0aJFFCtWjOnTp/PnP/9ZU/SIiGRDIItceeDnDNsJvn0ZDQO6m1kCsAD47wDmCQr16tVj3rx5FCxYkAkTJjBmzBivI4mIBC2vHyHoAnzonBttZrcAk8zsBudcasaDzKwX0AugdOnSxMXF5XzSHDZs2DA++ugjqlatell/b1JSUp5ov0BR+2WP2i/71IZZY4G6HOYrWsOccy1820MBnHOvZDhmO9DSOfezb3s30NA5t+98n1utWjWXV0YIcc5hZgCkpKRgZuTLl72T77i4OJo3b34Z0uVNar/sUftln9owc2a2wTl307n7A3m5ch1wjZlVNrMI0jqWzD3nmJ+AP/kCXgcUBPYHMFNQOVPgUlNTeeSRR+jZs6c6o4iIZEHALlc655LN7HFgERAGvO+c225mw4H1zrm5wABggpn9D2mdUB5y6mnxB99++y0zZszg+PHjJCcn88EHHxAWFuZ1LBGRXC+g9+SccwtI61CScd9zGdZ3AI0DmSEUXH/99SxYsIA2bdowadIkUlJSmDhxIuHhXt9SFRHJ3TTiSZBo1qwZCxcupGjRokyZMoVu3bpx+vRpr2OJiORqKnJBpEmTJnzxxRfpz9HFxsZqZBQRkQtQkQsyt9xyC4sXLyYqKorExEROnjzpdSQRkVxLN3WCUP369Vm+fDlXXXUVkZGRXscREcm1dCYXpGrVqkXx4sWBtGfoXnvtNU6cOOFtKBGRXMbvImdmhQMZRC7dwIEDGTBgAO3bt9c0PSIiGVy0yJlZIzPbAfzTt13bzPLsbAG5Uc+ePbniiitYvHgx7dq149ixY15HEhHJFfw5k3sdaAH8CuCc2wLcGshQkjXXX389cXFxlClThqVLl9KmTRuSkpK8jiUi4jm/LleeGVsyA40tlctcd911xMXFUa5cOZYtW0arVq04evSo17FERDzlT5H72cwaAc7M8pvZQODbAOeSS1CtWjWWLVuWPsP4008/DUBCQgIjRo7ixkbN6PZIb25s1IwRI0eRkJDgcWIRkcDyp8j1BvqSNhdcIlDHty25UNWqVVm2bBmdO3fm5ZdfZtWqVbTqEMu0jXuJbvkEVR9+leiWTzBt415adYhl1apVXkcWEQmYiz4n55w7AHTLgSxymVSpUoXp06eTkJDAY/0HE938YQqXKk/+QkVJPnqUQlFXENP4Xo5Uqslj/Qfz+aypxMTEeB1bROSy86d35Qdm9v65S06Ek+yZPOVjwirdzO4vp7P6zX6cOnb4rNeLla9KeJX6TJ7ysUcJRUQCy5/LlfOBz3zLEqAYoK57QWDG7PkUq1yLQz/u4Ejid6x+qx8nj/x21jElazRhxpzPPEooIhJYFy1yzrmZGZaPgPuAP8y+KrnP4cOHKVauCrf0e4uiV1bk6L92s3n8QJL2/aezbIFiJTh86JB3IUVEAuhShvW6BrjicgeRyy8qKooTh3+lYFQpbun3FsUrXseJg3tZ9fqfOfRjWgfZk0d+I8o3PJiISKjx557cUTM7cuZPYB7wVOCjSXZ1vqctv+5YCUCByGgaPv4G0dfcyKljh1n91hMc/HEHv+74is53t/E4qYhIYPhzuTLSOVcsw5/XOudm5kQ4yZ7uXbuQvHstRxLjAQgvUJgbejxP+ZtbUKzc1eAcybu/pnvXLh4nFREJjPM+QmBm9S70RufcxssfRy6nmJgY3n1jJI/1H8yRKvUpWaMJZhFUa92T/dviOLxiEu++MZIyZcp4HVVEJCAu9Jzc6Au85oDbL3MWCYBGjRrx+aypTJ7yMTPmvMkve/dyZZkydL27Dd3HDKFs2bLcd999VK5cmZEjR5Ivn2ZfEpHQcd4i55y7LSeDSODExMQwZPAghgweRFxcHM2bN09/be3atcybN4/Tp0+zd+9e3n//fSIiIrwLKyJyGfn1z3Yzu8HM7jOzB84sgQ4mOaNBgwZ89tlnFC1alI8++oh27dppBgMRCRn+9K58HnjLt9wGjATaBziX5KA777yTuLg4SpcuzRdffMFtt93Gvn37vI4lIpJt/pzJdQL+BOx1zj0M1AaiAppKctyNN97IqlWrqFKlCuvXr6dx48aapUBEgp4/Re5351wqkGxmxYB9QIXAxhIvVK1alZUrV1KnTh1iYmIoVaqU15FERLLlorMQAOvNrDgwAdhA2riVqwMZSrxTpkwZli1bhnOOggULAuCcw8w8TiYiknX+PAzexzl3yDn3N+BO4EHfZUsJUcWKFSMqKu2K9OnTp+nYsSPTp0/3OJWISNb50/Fkrpl1NbMizrkfnHNbcyKY5A4zZ87k008/JTY2lrfeesvrOCIiWeLPPbnRQBNgh5l9YmadzKxggHNJLnH//fczYsQInHP069ePp59+Guec17FERPziz+XKZc65PkAV4F3SptpR//I8wsx46qmn+OCDDwgLC+OVV17hkUceITk52etoIiIX5e/D4IWAjkBv4GZgYiBDSe7z0EMPMWfOHAoVKsQHH3zAPffcw/Hjx72OJSJyQf7ck5sOfEvaWJVvA1c75/470MEk92nTpg1Lly6lRIkSpKSkkD9/fq8jiYhckD+PEPwd6OKcSwl0GMn9GjZsyJo1ayhbtmx6kdMjBiKSW/lzT26RCpxkdM0111C0aFEATp06RatWrfjkk088TiUi8kcBnVfFzFqa2U4zizezIZm8/rqZbfYtu8zsUCDzyOX38ccfs2jRIjp37szw4cPV81JEcpWAFTkzCwPeAVoBNYAuZlYj4zHOuf9xztVxztUhbQDoWYHKI4HxwAMPMHr0aPLly8fzzz9PbGysOqSISK7hb+/K8mbWyMxuPbP48bb6QLxzbrdz7hQwFbj7Asd3AT72J4/kHmbGk08+ybx584iMjGT69Ok0bdpUgzuLSK5w0Y4nZvZX4H5gB3Dm3pwDll/kreWBnzNsJwANzvMdFYHKwNKL5ZHcqXXr1qxZs4Z27dqxceNGbr75ZtasWUPFihW9jiYieZg/vSvvAao5504GMEcs8Mn5OriYWS+gF0Dp0qWJi4sLYJTQlpSUFND2e+211xg2bBiRkZF8//337NmzJ2Df5YVAt1+oU/tln9owa/wpcruB/EBWi1wiZ0/JE+Pbl5lYoO/5Psg5Nx4YD1CtWjXXvHnzLEaRM+Li4gh0+7Vu3ZrTp09TuHBhAA4ePEhUVBT58gW0n1OOyIn2C2Vqv+xTG2aNP0XuOLDZzJaQodA55/pd5H3rgGvMrDJpxS0W6HruQWZWHYhG0/eEjPz586c/Q/f777/TokULypUrx6RJk4iMjPQ4nYjkJf4Uubm+JUucc8lm9jiwCAgD3nfObTez4cB659yZz4wFpjr1PQ9JO3fuJD4+nnXr1tG4cWPmzp1LpUqVvI4lInnERYucc26imUUA1/p27XTOnfbnw51zC4AF5+x77pztYf5FlWBUp04d1q5dS7t27di2bRs333wzn376KU2aNPE6mojkAf6MXdkc+I60Z97GArv8fIRABEgbIWXNmjW0aNGCAwcOcPvtt/P+++97HUtE8gB/55O7yznXzDl3K9ACeD2wsSTUFC9enPnz59O/f39Onz7NI488wuLFi72OJSIhzp97cvmdczvPbDjndpmZhp+XLAsPD+f111/n+uuvZ8WKFdxxxx1eRxKREOdPkVtvZu8Bk33b3YD1gYskoa5nz5707Nkzffunn37i1KlTVK1a1cNUIhKK/Llc+WfSRjvp51t2+PaJZNuxY8do37499evX1+VLEbns/Jlq56Rz7jXnXAff8nqARz+RPCQ1NZUKFSpw8OBBWrRowQsvvEBKimZ2EpHL47xFzjcjOGa2zcy2nrvkXEQJZZGRkcyePZvnn38egGHDhtG6dWv279/vcTIRCQUXOpN7wvdnW6BdJovIZREWFsawYcNYuHAhpUqV4osvvqBu3bqsXLnS62giEuTOW+Scc//2rfZxzv2YcQH65Ew8yUvuuusuNm/eTOPGjUlMTGTz5s1eRxKRIOdPx5M7M9nX6nIHEQEoX748X375JZMnT6ZPn//8W0qjvonIpbjQPbk/m9k2oNo59+P2ALonJwGTP39+unXrhpkBEB8fT7169diwYYPHyUQk2FzoTG4Kaffe5nL2vbgbnXPdcyCbCAAvv/wymzdvplGjRowbN05ndSLitwvdkzvsnPvBOdfFdx/ud9JmBC9qZlflWELJ8/72t7/Rp08fTp06RZ8+fejevTtJSUlexxKRIODPAM3tzOw7YA+wDPgB+DzAuUTSFShQgHfeeYcpU6ZQpEgRpkyZws0338z27du9jiYiuZw/HU9eAhoCu5xzlYE/AWsCmkokE126dGH9+vVcf/31/POf/6Rp06YcOXLE61gikov5M3blaefcr2aWz8zyOee+NLM3Ah1MJDPVq1dn7dq19OnTh4YNG1KsWDGvI4lILuZPkTtkZkWB5cBHZrYPOBbYWCLnV6RIET788MOz9i1atIirr75agzyLyFn8uVx5N3Ac+B9gIfA9aaOgiHjGzNIfMfjuu+/o3Lkz9erVY+bMmR4nE5HcxJ8i95xzLtU5l+ycm+icexN4KtDBRPx1xRVXcNddd3H06FE6depE//79OXXqlNexRCQX0IgnEvSioqKYMWMGb7zxBuHh4YwZM4Zbb72VPXv2eB1NRDzmz4gn1TMZ8WRbzkUUuTgz44knnmDFihVUqFCBtWvXUrt2bebMmeN1NBHxkD8jnszhjyOedMuBbCJZ1rBhQzZt2kTHjh05ceIElSpV8jqSiHjooiOeAGOA3zLMQJBsZg1yKqBIVpUsWZIZM2awadMmateunb5/y5YtHqYSES/4c09uHJBxDKUk3z6RXMvMuP7669O3p0+fTp06dejbty/Hjx/3MJmI5CR/ipy5DCPiOudS8e/5OpFcY9++feTPn5+xY8dSt25d1q1b53UkEckB/hS53WbWz8zy+5YngN2BDiZyOT3++OOsXbuWGjVqsGvXLho1asSLL75IcnKy19FEJID8KXK9gUZAIpAANAB6BTKUSCDUrVuXDRs20L9/f5KTk3nuuedo2rQpiYmJXkcTkQC5aJFzzu1zzsU6565wzl3pnOvqnNuXE+FELreCBQvy+uuvs3jxYsqXL8+BAwcoXry417FEJED8mWrnWjNbYmbf+LZrmdkzgY8mEjh33HEH27ZtY+7cuRQpUgSApKQk9u3Tv99EQok/lysnAEOB0wDOua1AbCBDieSE6OhorrvuuvTtAQMGcMMNNzBv3jwPU4nI5eRPkSvsnPv6nH26Wy8h5eTJk8THx7N//37at29Pr169NPu4SAjwp8gdMLOrAQdgZp2Afwc0lUgOK1CgAIsXL2b06NFEREQwYcIE6taty5o1mh9YJJj5U+T6Au+SNoZlItCftB6XIiElX758PPnkk6xfv56aNWsSHx9PkyZNeP7550lNTfU6nohcAn96V+52zt0BlAaqO+ea+Ib3EglJNWvWZN26dQwcOJDU1FS++eYbRr46mhsbNaPbI725sVEzRowcRUJCgtdRReQi/OldWdLM3gRWAHFmNsbMSvrz4WbW0sx2mlm8mQ05zzH3mdkOM9tuZlOyFl8kMAoUKMCoUaN4++23+XZPItM27iW65ROUv3sgxf70Z6Zt3EurDrGsWrXK66gicgH+DM81FVgOdPRtdwOmAXdc6E1mFga8Q9p8dAnAOjOb65zbkeGYa0jrudnYOXfQzK7I+k8QCYyEhATGfTiF0nf1plj5qqScPsX2t/tjZtTqMpj8lWryWP/BfD5rKjExMV7HFZFM+HNPrqxz7kXn3B7f8hJwpR/vqw/E+y53niKtWN59zjGPAu845w5C2oPnWQkvEkiTp3xMeJUGFCtfFYCTRw4AcGzfT6we8zg/fjUbi6nN5CkfexlTRC7AnyL3hZnFmlk+33IfsMiP95UHfs6wneDbl9G1wLVmttLM1phZS/9iiwTejNnzKVmjcfp24ZLluPHxt7mm5UNYWDg/rpzDriVTmfDBRA9TisiF+HO58lHSelRO8m2HAcfM7DHAOeeKZfP7rwGaAzHAcjOr6Zw7lPEgM+uFb7zM0qVLExcXl42vzNuSkpLUfn7a+8svFM0XQfLRo+n7XL58lGt6H1HX1mfnp2M4+vNOdh/5jXbt2jFgwAAP0wYH/f3LPrVh1ly0yDnnIi/xsxOBChm2Y3z7MkoA1jrnTgN7zGwXaUXvrHlQnHPjgfEA1apVc82bN7/ESBIXF4fazz9lrryS/KmnKBT1n1vFR48eJTIyksjImlw54F2+W/QP4hd9SKdOndSuftDfv+xTG2aNP70rHzlnO8zMnvfjs9cB15hZZTOLIG0osLnnHDObtLM4zKwUaZcvNY2P5Aqd72nLrztWnvd1yxdG4WJRPDX0aR544IH0/dOmTSM+Pj4nIorIRfhzT+5PZrbAzMqa2Q3AGuCiZ3fOuWTgcdLu330LTHfObTez4WbW3nfYIuBXM9sBfAkMcs79ekm/ROQy6961C8m713IkMfOCdSQxnuTdX9On92OYGQC7du3iwQcfpGbNmowaNUrz1Yl4zJ+HwbsCE4FtwGdAf+fcQH8+3Dm3wDl3rXPuaufcy759zznn5vrWnXPuSedcDedcTefc1Ev/KSKXV0xMDO++MZLflr5HwspZ/H5wHy4lhd8P7iNh5Sx+W/oe774x8qzHB0qUKEHnzp05ceIEgwcPpn79+mzcuNHDXyGSt/lzufIa4AlgJvAj0MPMCgc6mEhu0KhRIz6fNZX765Xl4KI3+f7DgRxc9Cb31yvL57Om0qhRo7OOL1WqFJMmTWLhwoVUrFiRTZs2Ub9+fZ566imOHz/u0a8Qybv8uVw5D3jWOfcY0Az4jnM6hoiEspiYGIYMHsSGlXFM/vvf2LAyjiGDB13wAfAWLVrwzTff0L9/f1JTUxk5ciS9e2vIV5Gc5k+Rq++cWwLplxdHA/cGNpZI8CtatCivv/46q1ev5uabb+aZZzTXsEhOO2+RM7PBAM65I2bW+ZyXHwpkKJFQ0qBBA9auXcu1114LgHOO++67j/fee0+zG4gE2IXO5DLO/j30nNc0MolIFpzpfQmwePFiZsyYwaOPPkqDBg00Z51IAF2oyNl51jPbFhE/3XnnnXz88ceUL1+e9evXc8stt/DQQw+xd+9er6OJhJwLFTl3nvXMtkXET2ZGbGws//znP3n66aeJiIhg4sSJXHvttYwdO9breCIh5UJFrraZHTGzo0At3/qZ7Zo5lE8kZBUtWpSXX36Z7du307ZtW44ePcpvv/3mdSyRkHLesSudc2E5GUQkr6patSrz5s1j0aJFNGvWLH3/4sWLqVq1KpUrV/YwnUhw8+cRAhHJAS1atKBgwYIAHDx4kK5du1KjRg2ef/55PUgucolU5ERyoeTkZFq0aMGJEycYPnw41113HZ988gnO6Xa4SFaoyInkQqVLl2by5MmsWLGCOnXq8NNPP9G5c2fuuOMOtm/f7nU8kaChIieSizVp0oT169czbtw4SpQowdKlS+nevbvO6ET8pCInksuFhYXRu3dvdu3aRZ8+fXjttdfSHy5PSkrSqCkiF6AiJxIkSpYsyTvvvMNtt92Wvu/xxx+nQYMGrFq1ysNkIrmXipxIkDp8+DBLlixh/fr1NG7cmHvvvZdvv/3W61giuYqKnEiQioqK4ttvv+WZZ56hcOHCzJ49mxtuuIFHH32UxMREr+OJ5AoqciJBrGjRorz44ovEx8fTu3dvzIz33nuP6tWr8+uvv3odT8RzKnIiIaBs2bKMGzeO7du307FjRzp37kzJkiWBtKl9Tp486XFCEW+oyImEkGrVqvHJJ58wfvz49H0LFy7k2muv5R//+AcpKSkephPJeSpyIiEoPPw/w9JOnDiRn376iQcffJB69erx+eef6zk7yTNU5ERC3EcffcSHH35IhQoV2Lp1K61bt+ZPf/oT69at8zqaSMCpyImEuLCwMB588EF27drFq6++SnR0NF9++SX169dn6tSpXscTCSgVOZE8omDBggwYMIDdu3fz1FNPcdVVV9GmTZv015OTkz1MJxIYKnIieUzx4sUZMWIE3333HZGRkQAcO3aM6667jueff56jR496nFDk8lGRE8mjIiIi0tfnz59PfHw8w4cP5+qrr2b06NEcO3bMw3Qil4eKnIhw//33s2LFCho1asT+/fsZOHAglStXZtSoUSQlJXkdT+SSqciJCJA2rc9XX33F/PnzqV+/Pvv372fw4MHcfffdXkcTuWQqciKSzsxo06YNa9as4fPPP6dhw4b06dMn/fUDBw5w5MgRDxOKZI2KnIj8gZnRsmVLVq1aRYcOHdL3P/3001SqVImXXnqJw4cPe5hQxD8qciJyXmaWPkFrSkoK33//PQcPHuTZZ5+lUqVKDB8+nEOHDnkbUuQCVORExC9hYWH83//9H0uWLOHWW2/l0KFDPP/881SqVIlhw4ap2EmupCInIn4zM26//XaWLVtGXFwct912G4cPH+aFF17gp59+8jqeyB+oyInIJWnWrBlLly5l+fLlvPjii9SqVSv9tXfffVfz2UmuENAiZ2YtzWynmcWb2ZBMXn/IzPab2Wbf0jOQeUTk8mvatCnPPPNM+nZcXBy9e/emUqVKDB06lAMHDniYTvK6gBU5MwsD3gFaATWALmZWI5NDpznn6viW9wKVR0RyRnR0NC1btiQpKYkRI0ZQsWJF+vXrx549e7yOJnlQIM/k6gPxzrndzrlTwFRAT5WKhLjatWvz+eefs2bNGlq3bs3x48d56623qFq1Kq+99prX8SSPCWSRKw/8nGE7wbfvXB3NbKuZfWJmFQKYR0RyUIMGDfjss8/YsmULPXr0IF++fOkDQgOcPn2a1NRUDxNKXmCBmiHYzDoBLZ1zPX3bPYAGzrnHMxxTEkhyzp00s8eA+51zt2fyWb2AXgClS5e+cfr06QHJnBckJSVRtGhRr2MELbXfpdu/fz+nTp2ifPm0f+t++umnzJkzh/vuu4877rjjrAGj5fz0dzBzt9122wbn3E3n7g9kkbsFGOaca+HbHgrgnHvlPMeHAb8556Iu9LnVqlVzO3fuvNxx84y4uDiaN2/udYygpfbLnozt16xZM5YvXw5AmTJl6NevH7179yY6OtrDhLmf/g5mzswyLXKBvFy5DrjGzCqbWQQQC8w9J1TZDJvtgW8DmEdEcpHFixfzj3/8g1q1arF3716efvppKlSoQP/+/fnxxx+9jichImBFzjmXDDwOLCKteE13zm03s+Fm1t53WD8z225mW4B+wEOByiMiuUtERAQ9evRg8+bNLFq0iDvvvJNjx44xZswYVq5c6XU8CRHhgfxw59wCYME5+57LsD4UGBrIDCKSu5kZd911F3fddRdbtmzh73//O507d05//UzPzJYtW6aPoyniL414IiK5Ru3atXnzzTfJnz8/APv27WPw4MG0bt2amjVr8uGHH3Ly5EmPU0owUZETkVyrUKFCDB8+nHLlyrF9+3YefvhhKlasyHPPPUdiYqLX8SQIqMiJSK4VGRnJoEGD2LNnDxMnTqRmzZr88ssvvPjii1SqVIl///vfXkeUXE5FTkRyvYiICB544AG2bNlCXFwcnTt35s4776Rs2bQO2s45PvroI81aLn8Q0I4nIiKXk5nRrFkzmjVrRkpKSvr+NWvW0L17d4oWLUqPHj3o27cv119/vYdJJbfQmZyIBKWwsLD09TPFLykpiXHjxnHDDTfQvHlzZsyYwenTpz1MKV5TkRORoNewYUPi4uLYtm0bf/7znylSpAjLli3jvvvuo169egRqZCfJ/VTkRCRk3HDDDYwdO5Z//etfvPXWW1SvXp277ror/fm6pKQkli9frqKXh6jIiUjIKVasGI8//jg7duzgpZdeSt//0Ucf0axZM2rXrs3f/vY3kpKSPEwpOUFFTkRClplRqFChs7avvPLK9Mua5cuXp2/fvmzYsEFndyFKRU5E8oxevXrx008/8fHHH9OkSROOHDnC2LFjuemmm3j00Ue9jicBoCInInlKREQEsbGxrFixgq1bt9K/f39KlixJ06ZN04/ZuXMnCxcuPOsxBQlOKnIikmfVrFmT119/ncTERGJjY9P3jxkzhlatWlGxYkWeeeYZvv/+ew9TSnaoyIlInlegQAEKFCiQvl29enWuvvpqEhMTefnll6latSq33XYbkyZN4vjx4x4mlaxSkRMROUe/fv347rvviIuL44EHHqBQoULp63/5y1+8jidZoCInIpKJM6OoTJw4kb179zJ+/HgaNmxIjx490o+ZM2cOo0ePZt++fR4mlQtRkRMRuYhixYrx6KOPsnr1aurVq5e+f9SoUQwcOJDy5cvToUMH5s+fr2HEchkVORGRSzRw4EDatWuHc45PP/2Udu3aUa5cOfr27cvWrVu9jieoyImIXLJ77rmHuXPn8vPPP/PXv/6V6tWrc+DAAcaOHcv69evTjzt16pSHKfM2FTkRkWwqW7YsgwcPZseOHWzatIlBgwbRoUOH9NeffPJJatWqxYgRI/jhhx+8C5oHqciJiFwmZkadOnUYOXIkxYsXB9ImdF2yZAnbtm1j6NChVK5cmSZNmjB27Fj279/vbeA8QEVORCSAzIwtW7Ywd+5cYmNjKVSoECtXrqRv376ULVuWcePGeR0xpKnIiYgEWEREBO3atePjjz9m3759TJ48mdatWwOc1Vtz6dKlzJs3T/fwLiMVORGRHFS0aFG6devGZ599xt69e6lfv376ay+++CLt27enTJkyPPbYY8TFxWn8zGxSkRMR8UipUqXSJ3R1ztGqVStq1qzJwYMHGT9+PLfddhvlypXjscceY8OGDR6nDU4qciIiuYCZMXjwYLZu3crWrVsZMmQIVapUYd++fYwfP57NmzenH3vgwAFOnDjhXdggoiInIpLL1KxZk1deeYX4+Hg2b97Ms88+y913353++gsvvEDp0qWJjY1l+vTpmuH8AsK9DiAiIpkzM2rXrk3t2rXP2v/jjz+SlJTEtGnTmDZtGgULFqRFixZ06NCBdu3aER0d7VHi3EdnciIiQWbu3Lns3r2b0aNH06hRI06cOMGcOXN48MEHGT58uNfxchUVORGRIFS5cmWefPJJVq5cSWJiIm+//Ta33347nTp1Sj/mvffeo3nz5rz11lskJiZ6mNY7KnIiIkHuzKDQS5YsoXHjxun7Z86cybJly+jXrx8xMTE0aNCAl156ic2bN+Oc8zBxzlGRExEJUVOnTmXy5Mnce++9FCpUiK+//ppnn32WunXrnjUvXihTkRMRCVFRUVF069aNWbNmsX//fubMmcOjjz5K2bJladSoUfpxX331FW3atGHcuHH89NNPHia+/NS7UkQkDyhSpAjt27enffv2OOfOmtx19uzZLFiwgAULFgBQq1Yt2rZtS9u2balfvz5hYWF/+LyEhAQmT/mYGbPnc/jwYaKiouh8T1u6d+1CTExMjv2uiwnomZyZtTSznWYWb2ZDLnBcRzNzZnZTIPOIiEjaowkRERHp24MGDeL999+nQ4cOFC1alK1bt/K///u/NGrUiAYNGvzh/atWraJVh1imbdxLdMsnqN5rDNEtn2Daxr206hDLqlWrcvLnXFDAipyZhQHvAK2AGkAXM6uRyXGRwBPA2kBlERGR87vyyit5+OGHmTlzJgcOHOCLL76gX79+VK5cmYYNG6Yf98svv9C4cWM6xnajYJ22xDS+l0LRV2D5wigUfQUxje+lxO09eaz/YBISEjz8Rf8RyDO5+kC8c263c+4UMBW4O5PjXgT+CmiMGhERjxUoUIA777yTMWPG8P333/Pqq6+mv/bZZ5+xatUq9v78A+snDGXp8Fi2TXuVf29ZzunjRwEoVr4q4VXqM3nKxxf8noSEBIY+/RdiKl9LoRJlKFTiSmKqXsfQvzx7WQtkIItceeDnDNsJvn3pzKweUME591kAc4iIyCUwMwoWLJi+fe+991LpmuqUqXUr+QsX4/iBRH5cOYcNf/8Li5+9h5RTJwEoWaMJ02fPP+/nzp49mzoNGjN6wmSSyt9MVMv+RN/zHMk39+BvS/9JkzvbXLZLnp51PDGzfMBrwEN+HNsL6AVQunRp4uLiApotlCUlJan9skHtlz1qv+zzug1PpTiu7TwIgKOJ8RyM38jB+I1YvnCOnzwFJ0+R6sLYum41TZs25aabbuKmm26ifPm0c5wVK1bw8mvv4AoUofitXSFfOIfWfELybwmQchrC8nO8SHHadLif8W+9RunSpbOV1wL1QKCZ3QIMc8618G0PBXDOveLbjgK+B86MLFoG+A1o75xbf77PrVatmtu5c2dAMucFcXFxNG/e3OsYQUvtlz1qv+zzug1vbNSM6JZPUCj6irP2u9RULF/axcF9337N1+MGnPV65cqVady4MfMWLeZUkbLkL3MN+QpFkrTpc4rUvIPCNW4lPLI0yUf3c3zHcpK2fsGN11Vh7Wr/zujMbINz7g+dFwN5uXIdcI2ZVTazCCAWmHvmRefcYedcKedcJedcJWANFylwIiLirc73tOXXHSv/sP9MgQM49VsCg4cMZcKECXTu3JkSJUqwZ88eJk+ezOH9v3D60F7CS8SQtOlziv+pJ1ENOhJRvCz5wsKJiCpD8Vvuo1S7QWzYEc/8+ee/7OmPgBU551wy8DiwCPgWmO6c225mw82sfaC+V0REAqd71y4k717LkcT4TF8/khhP8u6v+e++fejZsyfTp09n3759fP3115StUImwAoVJTT7NsR3LKFLzDo6uncnPb3bht/97N+0D0uaQpUC5ahSteScDnhqarbwBfU7OObfAOXetc+5q59zLvn3POefmZnJsc53FiYjkbjExMbz7xkh+W/oeCStn8fvBfaSmJPP7wX0krJzFb0vf4903Rp71QHhYWBg333wzhSOjiLzqOixfPpJ/S6BQ9caknjiGO/U7uFTf0UbaTTRH4epN2f3zv7OVVyOeiIhIljRq1IjPZ01NG/FkzpscPnSIqOLFuf/uNnR/dep5RzyJioqCyMoc/eVnXMpp8keVodwj75B8ZD+k9w/x/ZmaSnix0rjU1Ew/y18qciIikmUxMTEMGTyIIYMH+f2ezve0ZdJX3xEenp/TqckkH9lHRHQ5wov9pwelc2n/41wqKcd+O+te36XQAM0iIpIjunftQr69O6jS7F7cqd85vn0ZLiUZ51JxvsJGagouNRnLF8bxHcuoUqFstr5TRU5ERHLEmft5YQmbKFXlBpK2fcHJf+9KK3TJp3ApyQBYvnBO7f2OpG2LGf3XV7L1nSpyIiKSY87cz3vy4c5EFwrnwLxRHFk7k+SjvwJG8pF9HFk9jQNzR/I/vR6kbdu22fo+3ZMTEZEclfF+3vz58xnw1FB2T12Y/kB5lQplmTT579kucKAiJyIiHjozb12g6HKliIiELBU5EREJWSpyIiISslTkREQkZKnIiYhIyFKRExGRkBWwSVMDxcyOApo19dKVAg54HSKIqf2yR+2XfWrDzFV0zv1hGvFgfE5uZ2azv4p/zGy92u/Sqf2yR+2XfWrDrNHlShERCVkqciIiErKCsciN9zpAkFP7ZY/aL3vUftmnNsyCoOt4IiIi4q9gPJMTERHxS1AVOTNraWY7zSzezIZ4nSeYmNn7ZrbPzL7xOkswMrMKZvalme0ws+1m9oTXmYKJmRU0s6/NbIuv/V7wOlMwMrMwM9tkZvO9zhIsgqbImVkY8A7QCqgBdDGzGt6mCiofAi29DhHEkoEBzrkaQEOgr/7+ZclJ4HbnXG2gDtDSzBp6GykoPQF863WIYBI0RQ6oD8Q753Y7504BU4G7Pc4UNJxzy4HfvM4RrJxz/3bObfStHyXtPzTlvU0VPFyaJN9mft+iDgFZYGYxQBvgPa+zBJNgKnLlgZ8zbCeg/8iIB8ysElAXWOtxlKDiu9S2GdgHLHbOqf2y5g1gMJDqcY6gEkxFTsRzZlYUmAn0d84d8TpPMHHOpTjn6gAxQH0zu8HjSEHDzNoC+5xzG7zOEmyCqcglAhUybMf49onkCDPLT1qB+8g5N8vrPMHKOXcI+BLdI86KxkB7M/uBtFs1t5vZZG8jBYdgKnLrgGvMrLKZRQCxwFyPM0keYWYG/B341jn3mtd5go2ZlTaz4r71QsCdwD89DRVEnHNDnXMxzrlKpP23b6lzrrvHsYJC0BQ551wy8DiwiLSb/tOdc9u9TRU8zOxjYDVQzcwSzOwRrzMFmcZAD9L+Bb3Zt7T2OlQQKQt8aWZbSfsH62LnnLrBS8BpxBMREQlZQXMmJyIiklUqciIiErJU5EREJGSpyImISMhSkRMRkZClIieSgZk5MxudYXugmQ27yHt6m9kDAQ93EWbWPKuj05vZKN+sAKMu4fuezup7RHKaipzI2U4CHcyslL9vcM79zTn3jwBmCqReQC3n3KBLeG+Wipyl0X9zJEfpL5zI2ZKB8cD/nPuCmVUys6VmttXMlpjZVb79w8xsoG+9n2/Oua1mNtW3r4hvPr+vfXOB/WH2DDMra2bLfQ+Zf2NmTX37W5rZRt88bEt8++qb2WrfZ60ys2qZfJ4/3zkXKApsMLP7zaydma31Hf9/Znal77iiZvaBmW3z/a6OZjYCKOTL+5HvuCd92b8xs/4Z2mynmf0D+Iazh+YTCTznnBYtWnwLkAQUA34AooCBwDDfa/OAB33r/wXM9q0PAwb61v8FFPCtF/f9+b9A9zP7gF1AkXO+dwDwF996GBAJlCZt5o3Kvv0lfH8WA8J963cAM33rzYH5/n7nmd+bYT2a/wwQ0RMY7Vv/K/BGxuMyee+NwDagCGmFcztpMzVUIm3U/IZe/3+rJW8u4VkpiCJ5gXPuiO/Mox/we4aXbgE6+NYnASMzeftW4CMzmw3M9u27i7TBdQf6tgsCV3H25JfrgPd9g0DPds5tNrPmwHLn3B5frjPzAUYBE83sGtLmZMufSQ5/vvNcMcA0MysLRAB7fPvvIG28RHw5Dmby3ibAp865YwBmNgtoStr4sj8659Zc4HtFAkaXK0Uy9wbwCGlnJlnRhrQZ7OsB68wsHDCgo3Oujm+5yjl3VrFxaZPa3krazBofXqQjy4vAl865G4B2pBWwc130OzPxFvC2c64m8Nh5PvdSHLtMnyOSZSpyIpnwnTVNJ63QnbGK/5zRdANWZHyPr1NFBefcl8BTpJ1xFSVtUPH/9s1kgJnVPff7zKwi8ItzbgJpMz/XA9YAt5pZZd8xJXyHR/GfaaYeOs9PuOh3ZiLj5z6YYf9ioG+GrNG+1dO+M09Ia4t7zKywmRUB7uWc9hHxgoqcyPmNBjL2svxv4GHfSPo9gCfOOT4MmGxm24BNwJsube60F0m7pLjVzLb7ts/VHNhiZpuA+4Exzrn9pPV+nGVmW4BpvmNHAq/4jj3fLQd/vvNcw4AZZrYBOJBh/0tAtK9DyRbgNt/+8b7P/8g5txH4EPiatBnT33PObfLjO0UCSrMQiIhIyNKZnIiIhCwVORERCVkqciIiErJU5EREJGSpyImISMhSkRMRkZClIiciIiFLRU5ERELW/wP589sI44M4qAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = factory.plot_fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "b47ba6ea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.75585, 0.59994, 0.40038, 0.39116, 0.39232, 0.39106, 0.38722,\n",
       "       0.39403, 0.38875, 0.38865])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "factory.get_expectation_values()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "4b9ff3ed",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.        , 2.        , 4.46956499, 4.65002695, 4.66240918,\n",
       "       4.67994465, 4.68665802, 4.67488215, 4.69063002, 4.68714552])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "factory.get_scale_factors()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "09b98717",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9606057483296139"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "factory.get_zero_noise_limit()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fd32caa",
   "metadata": {},
   "source": [
    "## Section 3:  Applying probabilistic error cancellation (PEC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "39516289",
   "metadata": {},
   "outputs": [],
   "source": [
    "from mitiq import pec"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fe2687b",
   "metadata": {},
   "source": [
    "### Build quasi-probability representation for CNOT gate"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "086c9a8b",
   "metadata": {},
   "source": [
    "For simplicity, we assume that:\n",
    "    \n",
    "* Only CNOT gates are affected by noise;\n",
    "* The noise is depolarizing.\n",
    "\n",
    "Even within these assumptions, we still need to know how strong the noise is. A good approximation is the error probability of the CNOT gate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "efde7be4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.008339674869236618"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Print error probability of CNOT gate\n",
    "cnot_error_prob = noisy_backend.properties().gate_error(\"cx\", qubits=(0, 1))\n",
    "cnot_error_prob\n",
    "\n",
    "# print(qiskit.transpile(circuit, backend=noisy_backend))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46f4eef1",
   "metadata": {},
   "source": [
    "We now build a quasi-probability representation based on the previous assumptions. I.e., we represent the ideal CNOT as a linear combination of noisy operations:\n",
    "\n",
    "$$ \\mathcal G_{\\rm CNOT} = \\sum_\\alpha \\eta_\\alpha \\mathcal O_\\alpha, \\quad \\eta_\\alpha \\in \\mathbb R, \\quad \\{ \\mathcal O_\\alpha \\} = \\text{Noisy implementable operations}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "5c93da56",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Quasi-probability representation with one-norm 1.01683 generated.\n"
     ]
    }
   ],
   "source": [
    "from mitiq.pec.representations.depolarizing import represent_operation_with_global_depolarizing_noise\n",
    "\n",
    "# Define ideal operation that we aim to represent\n",
    "cnot_circuit = qiskit.QuantumCircuit(qiskit.QuantumRegister(2, name=\"q\")) \n",
    "cnot_circuit.cx(0, 1)\n",
    "\n",
    "# Build representation\n",
    "cnot_rep = represent_operation_with_global_depolarizing_noise(cnot_circuit, noise_level=cnot_error_prob)\n",
    "\n",
    "print(f\"Quasi-probability representation with one-norm {cnot_rep.norm :g} generated.\")\n",
    "# print(cnot_rep)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48c9a537",
   "metadata": {},
   "source": [
    "### Run PEC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "29d77951",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 214 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 214 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 213 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 213 gates using 1000 shots.\n",
      "Executing circuit with 213 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 212 gates using 1000 shots.\n",
      "Executing circuit with 213 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 215 gates using 1000 shots.\n",
      "Executing circuit with 212 gates using 1000 shots.\n",
      "Executing circuit with 214 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 213 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 213 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 214 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 213 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 212 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 213 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 213 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 213 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Estimation error without Mitiq   : 0.24211000000000005\n",
      "Estimation error with Mitiq (PEC): 0.0006080661539451526\n"
     ]
    }
   ],
   "source": [
    "from functools import partial\n",
    "\n",
    "# Set number of probabilistically generated circuits\n",
    "num_samples = 100\n",
    "\n",
    "# Distribute budget of shots among all circuits\n",
    "shots_per_circuit = shots // num_samples\n",
    "pec_executor = partial(executor, shots=shots_per_circuit)\n",
    "\n",
    "# Run PEC\n",
    "pec_value = pec.execute_with_pec(\n",
    "    circuit, pec_executor, representations = [cnot_rep], num_samples=num_samples, random_state=0,\n",
    ")\n",
    "\n",
    "unmitigated_error = abs(ideal_value - noisy_value)\n",
    "pec_error = abs(ideal_value - pec_value)\n",
    "print(f\"Estimation error without Mitiq   : {unmitigated_error}\")\n",
    "print(f\"Estimation error with Mitiq (PEC): {pec_error}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "6641c820",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PEC reduced the estimation error by 99.7%.\n"
     ]
    }
   ],
   "source": [
    "print(f\"PEC reduced the estimation error by {(unmitigated_error - pec_error) / unmitigated_error :.1%}.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e562e94",
   "metadata": {},
   "source": [
    "### Build a noise-agnostic quasi-probability representation for CNOT gate\n",
    "\n",
    "Instead of a assuming a depolarizing nose model with some given error rate, we now try to use a noise-agnostic representation based on the NEPEC framework introduced in  https://arxiv.org/abs/2108.02237 (Section V). \n",
    "More details can be found in the notebook at [this link](https://github.com/unitaryfund/research/blob/main/nepec/noise-agnostic-pec/noise_agnostic_pec.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "70852158",
   "metadata": {},
   "outputs": [],
   "source": [
    "from mitiq.pec import NoisyOperation, OperationRepresentation\n",
    "\n",
    "\n",
    "def get_richardson_coefficients(scale_factors):\n",
    "    \"\"\"Returns the array of Richardson extrapolation coefficients associated\n",
    "    to the input array of scale factors.\"\"\"\n",
    "    \n",
    "    # This is based on the Lagrange interpolation formula\n",
    "    richardson_coeffs = []\n",
    "    for l in scale_factors:\n",
    "        coeff = 1.0\n",
    "        for l_prime in scale_factors:\n",
    "            if l_prime == l:\n",
    "                continue\n",
    "            coeff *= l_prime / (l_prime - l)\n",
    "        richardson_coeffs.append(coeff)\n",
    "\n",
    "    return richardson_coeffs\n",
    "\n",
    "def represent_operation_with_nepec(ideal_operation, scale_factors):\n",
    "    \"\"\"Returns the OperationRepresentation contructed by noise scaling and extrapolation.\"\"\"\n",
    "    \n",
    "    coeffs = get_richardson_coefficients(scale_factors)\n",
    "    \n",
    "    basis_expansion = {}\n",
    "    for scale_factor, coeff in zip(scale_factors, coeffs):\n",
    "        scaled_op = NoisyOperation(noise_scaling_function(ideal_operation, scale_factor))\n",
    "        basis_expansion[scaled_op] = coeff\n",
    "    \n",
    "    return OperationRepresentation(ideal_operation, basis_expansion)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "02e1caf7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NEPEC quasi-probability representation with one-norm 1.02 generated.\n"
     ]
    }
   ],
   "source": [
    "scale_factors = [1, 101]  # Large scale factors reduce the PEC sampling cost.\n",
    "# scale_factors = [1, 3, 9]\n",
    "\n",
    "agn_rep = represent_operation_with_nepec(cnot_circuit, scale_factors)\n",
    "print(f\"NEPEC quasi-probability representation with one-norm {agn_rep.norm :g} generated.\")\n",
    "# print(agn_rep)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a8ce8cc",
   "metadata": {},
   "source": [
    "### Run NEPEC (noise-extended PEC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "d2bfc8ca",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 411 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 311 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 411 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 311 gates using 1000 shots.\n",
      "Executing circuit with 511 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 311 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 311 gates using 1000 shots.\n",
      "Executing circuit with 311 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 311 gates using 1000 shots.\n",
      "Executing circuit with 311 gates using 1000 shots.\n",
      "Executing circuit with 311 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 411 gates using 1000 shots.\n",
      "Executing circuit with 311 gates using 1000 shots.\n",
      "Executing circuit with 411 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 311 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 311 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 411 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 311 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 311 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 311 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 311 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Executing circuit with 311 gates using 1000 shots.\n",
      "Executing circuit with 211 gates using 1000 shots.\n",
      "Estimation error without Mitiq     : 0.24211000000000005\n",
      "Estimation error with Mitiq (NEPEC): 0.011558459442957192\n"
     ]
    }
   ],
   "source": [
    "nepec_value = pec.execute_with_pec(\n",
    "    circuit, pec_executor, representations = [agn_rep], num_samples=num_samples, random_state=0,\n",
    ")\n",
    "\n",
    "unmitigated_error = abs(ideal_value - noisy_value)\n",
    "nepec_error = abs(ideal_value - nepec_value)\n",
    "print(f\"Estimation error without Mitiq     : {unmitigated_error}\")\n",
    "print(f\"Estimation error with Mitiq (NEPEC): {nepec_error}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "5b2cfdfd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NEPEC reduced the estimation error by 95.2%.\n"
     ]
    }
   ],
   "source": [
    "print(f\"NEPEC reduced the estimation error by {(unmitigated_error - nepec_error) / unmitigated_error :.1%}.\")"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "text_representation": {
    "extension": ".myst",
    "format_name": "myst",
    "format_version": 0.13,
    "jupytext_version": "1.11.1"
   }
  },
  "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}