mitiq Intro.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install mitiq --quiet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "warnings.filterwarnings(action='ignore') # Optional warning filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit import IBMQ\n",
    "IBMQ.save_account('0238b0afc0dc515fe7987b02706791d1719cb89b68befedc125eded0607e6e9e9f26d3eed482f66fdc45fdfceca3aab2edb9519d96b39e9c78040194b86e7858', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Section 1: Computing quantum expectation values without error mitigation\n",
    "## Define the circuit of interest\n",
    "We first define an arbitrary circuit of interest. For simplicity, in this example we use a two-qubit randomized-benchmarking circuit whose ideal execution is equivalent to the identity operation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Randomized benchmarking circuit with 211 gates generated.\n",
      "     ┌───┐┌─────┐┌───┐          ┌───┐┌───┐ ┌───┐ ┌───┐┌───┐     ┌─────┐┌───┐»\n",
      "q_0: ┤ Y ├┤ Sdg ├┤ H ├───────■──┤ H ├┤ S ├─┤ X ├─┤ H ├┤ S ├──■──┤ Sdg ├┤ H ├»\n",
      "     ├───┤└┬───┬┘├───┤┌───┐┌─┴─┐├───┤├───┤┌┴───┴┐├───┤└───┘┌─┴─┐└┬───┬┘├───┤»\n",
      "q_1: ┤ Z ├─┤ S ├─┤ H ├┤ S ├┤ X ├┤ H ├┤ S ├┤ Sdg ├┤ H ├─────┤ X ├─┤ H ├─┤ S ├»\n",
      "     └───┘ └───┘ └───┘└───┘└───┘└───┘└───┘└─────┘└───┘     └───┘ └───┘ └───┘»\n",
      "«          ┌───┐┌───┐┌───┐      ┌───┐ ┌───┐     ┌───┐┌───┐┌───┐          ┌───┐»\n",
      "«q_0: ──■──┤ X ├┤ H ├┤ S ├──■───┤ H ├─┤ S ├──■──┤ X ├┤ H ├┤ S ├───────■──┤ S ├»\n",
      "«     ┌─┴─┐├───┤├───┤└───┘┌─┴─┐┌┴───┴┐├───┤┌─┴─┐├───┤├───┤├───┤┌───┐┌─┴─┐├───┤»\n",
      "«q_1: ┤ X ├┤ X ├┤ H ├─────┤ X ├┤ Sdg ├┤ H ├┤ X ├┤ H ├┤ S ├┤ H ├┤ S ├┤ X ├┤ Y ├»\n",
      "«     └───┘└───┘└───┘     └───┘└─────┘└───┘└───┘└───┘└───┘└───┘└───┘└───┘└───┘»\n",
      "«     ┌───┐┌───┐          ┌─────┐┌───┐┌───┐                 ┌─────┐┌───┐┌───┐»\n",
      "«q_0: ┤ H ├┤ S ├───────■──┤ Sdg ├┤ H ├┤ X ├──────────────■──┤ Sdg ├┤ H ├┤ H ├»\n",
      "«     ├───┤├───┤┌───┐┌─┴─┐└┬───┬┘├───┤├───┤┌─────┐┌───┐┌─┴─┐├─────┤├───┤├───┤»\n",
      "«q_1: ┤ S ├┤ H ├┤ S ├┤ X ├─┤ H ├─┤ S ├┤ Z ├┤ Sdg ├┤ H ├┤ X ├┤ Sdg ├┤ H ├┤ Z ├»\n",
      "«     └───┘└───┘└───┘└───┘ └───┘ └───┘└───┘└─────┘└───┘└───┘└─────┘└───┘└───┘»\n",
      "«               ┌─────┐┌───┐     ┌───┐┌─────┐┌───┐          ┌─────┐┌───┐     »\n",
      "«q_0: ───────■──┤ Sdg ├┤ H ├──■──┤ Y ├┤ Sdg ├┤ H ├───────■──┤ Sdg ├┤ H ├─────»\n",
      "«     ┌───┐┌─┴─┐├─────┤├───┤┌─┴─┐├───┤└┬───┬┘├───┤┌───┐┌─┴─┐└┬───┬┘├───┤┌───┐»\n",
      "«q_1: ┤ S ├┤ X ├┤ Sdg ├┤ H ├┤ X ├┤ H ├─┤ S ├─┤ Y ├┤ H ├┤ X ├─┤ H ├─┤ S ├┤ X ├»\n",
      "«     └───┘└───┘└─────┘└───┘└───┘└───┘ └───┘ └───┘└───┘└───┘ └───┘ └───┘└───┘»\n",
      "«           ┌───┐ ┌───┐     ┌───┐┌─────┐┌───┐     ┌───┐┌───┐┌───┐┌───┐┌───┐»\n",
      "«q_0: ──■───┤ H ├─┤ S ├──■──┤ X ├┤ Sdg ├┤ H ├──■──┤ H ├┤ S ├┤ X ├┤ S ├┤ H ├»\n",
      "«     ┌─┴─┐┌┴───┴┐├───┤┌─┴─┐├───┤└┬───┬┘└───┘┌─┴─┐├───┤├───┤├───┤└───┘└───┘»\n",
      "«q_1: ┤ X ├┤ Sdg ├┤ H ├┤ X ├┤ Z ├─┤ H ├──────┤ X ├┤ X ├┤ H ├┤ S ├──────────»\n",
      "«     └───┘└─────┘└───┘└───┘└───┘ └───┘      └───┘└───┘└───┘└───┘          »\n",
      "«     ┌───┐     ┌───┐┌───┐┌───┐           ┌───┐ ┌───┐     ┌───┐┌─────┐┌───┐»\n",
      "«q_0: ┤ S ├──■──┤ Y ├┤ H ├┤ S ├───────■───┤ H ├─┤ S ├──■──┤ Z ├┤ Sdg ├┤ H ├»\n",
      "«     └───┘┌─┴─┐├───┤├───┤├───┤┌───┐┌─┴─┐┌┴───┴┐├───┤┌─┴─┐├───┤└┬───┬┘├───┤»\n",
      "«q_1: ─────┤ X ├┤ H ├┤ S ├┤ Y ├┤ H ├┤ X ├┤ Sdg ├┤ H ├┤ X ├┤ H ├─┤ S ├─┤ Y ├»\n",
      "«          └───┘└───┘└───┘└───┘└───┘└───┘└─────┘└───┘└───┘└───┘ └───┘ └───┘»\n",
      "«                      ┌───┐┌───┐     ┌───┐               ┌─────┐┌───┐     »\n",
      "«q_0: ──────────────■──┤ H ├┤ S ├──■──┤ Y ├────────────■──┤ Sdg ├┤ H ├──■──»\n",
      "«     ┌─────┐┌───┐┌─┴─┐├───┤├───┤┌─┴─┐├───┤┌───┐┌───┐┌─┴─┐└┬───┬┘├───┤┌─┴─┐»\n",
      "«q_1: ┤ Sdg ├┤ H ├┤ X ├┤ H ├┤ S ├┤ X ├┤ X ├┤ H ├┤ S ├┤ X ├─┤ H ├─┤ S ├┤ X ├»\n",
      "«     └─────┘└───┘└───┘└───┘└───┘└───┘└───┘└───┘└───┘└───┘ └───┘ └───┘└───┘»\n",
      "«     ┌───┐┌───┐┌─────┐┌───┐     ┌───┐┌───┐┌───┐┌───┐     ┌─────┐┌───┐     »\n",
      "«q_0: ┤ H ├┤ S ├┤ Sdg ├┤ H ├──■──┤ X ├┤ S ├┤ H ├┤ S ├──■──┤ Sdg ├┤ H ├──■──»\n",
      "«     ├───┤├───┤└─────┘└───┘┌─┴─┐├───┤├───┤├───┤└───┘┌─┴─┐└┬───┬┘├───┤┌─┴─┐»\n",
      "«q_1: ┤ X ├┤ S ├────────────┤ X ├┤ S ├┤ H ├┤ S ├─────┤ X ├─┤ H ├─┤ S ├┤ X ├»\n",
      "«     └───┘└───┘            └───┘└───┘└───┘└───┘     └───┘ └───┘ └───┘└───┘»\n",
      "«     ┌───┐┌───┐┌───┐┌─────┐┌───┐     ┌───┐┌───┐┌───┐┌─────┐┌───┐     ┌─────┐»\n",
      "«q_0: ┤ H ├┤ S ├┤ Y ├┤ Sdg ├┤ H ├──■──┤ H ├┤ S ├┤ X ├┤ Sdg ├┤ H ├──■──┤ Sdg ├»\n",
      "«     ├───┤├───┤└───┘└─────┘└───┘┌─┴─┐├───┤├───┤├───┤└┬───┬┘└───┘┌─┴─┐└┬───┬┘»\n",
      "«q_1: ┤ Y ├┤ S ├─────────────────┤ X ├┤ H ├┤ S ├┤ Y ├─┤ H ├──────┤ X ├─┤ H ├─»\n",
      "«     └───┘└───┘                 └───┘└───┘└───┘└───┘ └───┘      └───┘ └───┘ »\n",
      "«     ┌───┐     ┌─────┐┌───┐     ┌─────┐┌───┐     ┌───┐ ┌───┐ ┌───┐┌───┐     »\n",
      "«q_0: ┤ H ├──■──┤ Sdg ├┤ H ├──■──┤ Sdg ├┤ H ├──■──┤ H ├─┤ S ├─┤ Z ├┤ S ├──■──»\n",
      "«     ├───┤┌─┴─┐└─────┘└───┘┌─┴─┐└┬───┬┘├───┤┌─┴─┐├───┤┌┴───┴┐├───┤└───┘┌─┴─┐»\n",
      "«q_1: ┤ S ├┤ X ├────────────┤ X ├─┤ H ├─┤ S ├┤ X ├┤ Z ├┤ Sdg ├┤ H ├─────┤ X ├»\n",
      "«     └───┘└───┘            └───┘ └───┘ └───┘└───┘└───┘└─────┘└───┘     └───┘»\n",
      "«      ┌───┐ ┌───┐     ┌───┐┌───┐\n",
      "«q_0: ─┤ H ├─┤ S ├──■──┤ H ├┤ Y ├\n",
      "«     ┌┴───┴┐├───┤┌─┴─┐├───┤└───┘\n",
      "«q_1: ┤ Sdg ├┤ H ├┤ X ├┤ X ├─────\n",
      "«     └─────┘└───┘└───┘└───┘     \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",
    "\n",
    "# print(qiskit_circuit)\n",
    "\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",
   "metadata": {},
   "source": [
    "# Execute circuit with noiseless backend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 13,
     "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",
   "metadata": {},
   "source": [
    "## Executing circuit on a noisy simulated backend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Executing circuit with 211 gates using 20000 shots.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.75465"
      ]
     },
     "execution_count": 28,
     "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 = 20000\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",
   "metadata": {},
   "source": [
    "## Section 2: Applying zero-noise extrapolation (ZNE)\n",
    "\n",
    "### Define an executor function\n",
    "We rewrite the last code cell of the previous section as a function that we call executor. This function takes as input a Qiskit circuit and returns the final noisy expectation value. The executor is viewed as a blackbox by Mitiq.\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": 36,
   "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": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Executing circuit with 211 gates using 100000 shots.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.75832"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "noisy_value = executor(circuit)\n",
    "noisy_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "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 633 gates using 100000 shots.\n",
      "Estimation error without Mitiq   : 0.24397000000000002\n",
      "Estimation error with Mitiq (ZNE): 0.028410000000000046\n",
      "ZNE reduced the estimation error by 88.4%.\n"
     ]
    }
   ],
   "source": [
    "from mitiq import zne\n",
    "\n",
    "# By default, noise is scaled by locally folding gates at random, and the default extrapolation is Richardson. Scale factor = 3.\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}\")\n",
    "print(f\"ZNE reduced the estimation error by {(unmitigated_error - zne_error) / unmitigated_error :.1%}.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Decorator Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mitiq\n",
    "\n",
    "@mitiq.zne.zne_decorator()\n",
    "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": 35,
   "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 633 gates using 100000 shots.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.0453900000000003"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "noisy_value = executor(circuit)\n",
    "noisy_value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select a noise scaling method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose a unitary folding function\n",
    "noise_scaling_function = zne.scaling.fold_global"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The input circuit has 211 gates\n",
      "The scaled circuit has 739 gates\n"
     ]
    }
   ],
   "source": [
    "noise_scaled_circuit = noise_scaling_function(circuit, scale_factor=3.5)\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",
   "metadata": {},
   "source": [
    "### Select Extrapolation Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Examples:\n",
    "\n",
    "# Richardson extrapolation\n",
    "# factory = zne.inference.RichardsonFactory(scale_factors = [1, 1.2, 2, 3, 5.5])\n",
    "\n",
    "# Exponential extrapolation\n",
    "factory = zne.inference.ExpFactory(scale_factors = [1, 1.5, 2, 3, 3.5], asymptote=0.25)\n",
    "\n",
    "# Exponential adaptive extrapolation\n",
    "# factory = zne.inference.AdaExpFactory(steps=10, asymptote=0.25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Executing circuit with 211 gates using 100000 shots.\n",
      "Executing circuit with 317 gates using 100000 shots.\n",
      "Executing circuit with 423 gates using 100000 shots.\n",
      "Executing circuit with 633 gates using 100000 shots.\n",
      "Executing circuit with 739 gates using 100000 shots.\n",
      "Estimation error without Mitiq   : 0.24168\n",
      "Estimation error with Mitiq (ZNE): 0.03625839672298803\n",
      "ZNE reduced the estimation error by 85.0%.\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}\")\n",
    "print(f\"ZNE reduced the estimation error by {(unmitigated_error - zne_error) / unmitigated_error :.1%}.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbkAAAE9CAYAAACWdRzmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA3sklEQVR4nO3deZzNdf//8cfLjLUxlkIYSdbIxUVZJpcsF1FSGUJJRI1KuMSkfhVX5VuRErlSKVSylUoiV1nKlq3FWEuyjMKVfZkwvH9/nGMa0+BgznzmnHneb7fPzWc75zzfTp2Xz/Z+m3MOERGRcJTL6wAiIiLBoiInIiJhS0VORETCloqciIiELRU5EREJWypyIiIStiK9DnC+Chcu7CpUqOB1jEx3+PBhLrnkEq9jZDq1K7SEY7vCsU2gdqW3cuXK351zxdKvD7kiV6JECVasWOF1jEw3f/58GjVq5HWMTKd2hZZwbFc4tgnUrvTMbEtG63W6UkREwpaKnIiIhC0VORERCVsqciIiErZU5EREJGypyImISNhSkRMRkbClIiciImFLRU5ERMJWyBW5lJQUTp486XUMEREJASFX5LZs2UK/fv1wznkdRUREsrmQK3InTpzg5Zdf5sUXX/Q6ioiIZHMhV+RKliwJQEJCAuPHj/c4jYiIZGchV+QKFizIK6+8AkC3bt2YOXOmx4lERCS7CrkiB9CrVy8GDBjAiRMnaNeuHd98843XkUREJBsKySIH8H//93907dqV5ORkVq1a5XUcERHJhkJu0NRTzIw33niDzp07h+XAgSIicvFC9kgOIDIy8rQC9+OPP7J3717vAomISLYS0kUurZUrVxIbG0vr1q1JTk72Oo6IiGQDYVPkihcvTr58+Vi4cCEdOnQgJSXF60giIuKxsClyZcqUYfbs2RQuXJjp06fzwAMPqFcUEZEcLmyKHEC1atWYMWMG+fLlY8yYMQwcONDrSCIi4qGwKnIA119/PZMnTyZXrlw888wz/Oc///E6koiIeCTsihxA69ateeONN4iIiCBPnjxexxEREY+E7HNy59KtWzcaNmxIxYoVvY4iIiIeCcsjuVPSFrjExER++OEHD9OIiEhWC9sjubQSExP5xz/+Qf78+Vm8eDHlypXzOpKIiGSBsD6SO6VSpUrUrl2bHTt20Lx5c3bt2uV1JBERyQI5osjlzZuXjz76iL///e9s3LiRm2++mUOHDnkdS0REgixHFDmA6OhoZs6cSbly5VixYgVxcXEcO3bM61giIhJEOabIAVx++eX897//pVixYvz3v//l3nvvVa8oIiJhLEcVOYAKFSowa9YsihYtSsuWLTEzryOJiEiQ5Ii7K9OrXbs2mzZtolChQl5HERGRIMpxR3KnpC1wy5YtY9iwYR6mERGRYMiRR3Jp7dmzhxtvvJF9+/ZhZvTt29frSCIikkly7JHcKUWLFuWll14C4JFHHlGHziIiYSTHFzmArl27MmrUKAAeeugh3n77bY8TiYhIZlCR83vwwQdTr8t1796d999/3+NEIiJysYJa5MyshZltMLONZjYgg+1lzWyOma0ys/lmFhPMPOfSt29fBg8ejHOOe++9l+3bt3sZR0RELlLQbjwxswhgFNAMSAKWm9l059zaNLu9CLzjnBtvZk2A54C7g5UpEI8//jjHjx/nmmuuoXTp0l5GERGRixTMuyvrABudc5sAzGwScCuQtshVBU7dzjgP+DiIeQI2cODA05aTk5PJnz+/R2lERORCWbC6tTKztkAL51x3//LdQF3nXM80+7wPLHXOvWJmbYAPgcucc7vTvdf9wP0AxYoVqz1lypSgZM5IYmIigwYN4sknn6RmzZpB+5xDhw4RFRUVtPf3itoVWsKxXeHYJlC70mvcuPFK59y1f9ngnAvKBLQFxqRZvht4Nd0+pYBpwHfAK/hOaxY+2/tWqlTJZaWHH37YAS4qKsotXrw4aJ8zb968oL23l9Su0BKO7QrHNjmndqUHrHAZ1Ixg3niyHSiTZjnGvy5tgf3VOdfGOfd34P/51+0LYqbz9vLLL3PXXXdx6NAhWrZsybfffut1JBERCVAwi9xyoKKZlTOzPEAHYHraHczsMjM7leExINs9oBYREcG4ceNo27Yt+/fvp1mzZiQmJnodS0REAhC0IuecSwF6ArOBdcAU59waM3vazFr7d2sEbDCzH4ESwOBg5bkYkZGRTJgwgVtuuYU9e/bQtGlT1q9f73UsERE5h6D2XemcmwnMTLfuqTTzHwAfBDNDZsmTJw9TpkyhdevWzJs3j/Xr11OlShWvY4mIyFnk+A6az0e+fPn4+OOPWb58OTfccIPXcURE5BzUrdd5KlCgwGkFbtmyZeoZRUQkm1KRuwhLly6ladOmNG3alJ07d3odR0RE0lGRuwgVKlSgXLlybNiwgWbNmrF79+5zv0hERLKMitxFuPTSS/nyyy+pUqUKiYmJNG/enH379nkdS0RE/FTkLlLx4sWZM2cO5cuX59tvv+XGG29k7969XscSERFU5DJFqVKlmDt3LuXKlWPZsmW0bNmSkydPeh1LRCTHU5HLJFdccQVfffUVVatWZcCAAeTKpb9aERGv6Tm5TFSmTBl++OEHIiP//Gs9efKkCp6IiEf065vJ0ha4BQsWUKtWLbZt2+ZhIhGRnEtFLkicczz11FP88MMPNGzYkF9++cXrSCIiOY6KXJCYGR999BF169Zl8+bNNGzYkJ9++snrWCIiOYqKXBAVLlyY//73vzRo0ICkpCQaNmzI2rVrvY4lIpJjqMgFWXR0NJ9//jlNmjRhx44dNGrUiFWrVnkdS0QkR1CRywKXXHIJM2bMoEWLFuzbt4+kpCSvI4mI5Ah6hCCL5M+fP3WYngYNGngdR0QkR9CRXBbKmzfvaQVu3rx5fP311x4mEhEJbzqS88jatWtp1aoVANOnTyciIsLjRCIi4UdHch6pXLkyd9xxB0eOHKFVq1YsW7bM60giImFHRc4jERERvPXWW8THx/PHH3/wxBNP8Omnn3odS0QkrKjIeShXrly89tpr9OrVi+PHj9OmTRs+/PBDr2OJiIQNFTmPmRnDhw+nffv2pKSk0LlzZ3bt2uV1LBGRsKAbT7IBMyM+Pp6rr76aWrVqUbx4ca8jiYiEBRW5bMLMGDhw4Gnrfv31V0qVKuVRIhGR0KfTldnU/PnzqVixIi+99JLXUUREQpaKXDa1fv16jhw5wiOPPMKjjz6Kc87rSCIiIUdFLpvq0aMH7777LpGRkQwZMoR7772XlJQUr2OJiIQUFblsrFOnTkyfPp0CBQowbtw4br/9do4cOeJ1LBGRkKEil821bNmSuXPnUrRoUWbMmMHtt9+uU5ciIgFSkQsBdevWZeHChZQvX56+fftiZl5HEhEJCXqEIERcffXVrF27ljx58qSuS05OJn/+/B6mEhHJ3nQkF0LSFrgvvviC8uXLs3TpUg8TiYhkbypyIertt9/mt99+o0mTJnz++edexxERyZZU5ELUO++8wz333MORI0e45ZZbmDBhgteRRESyHRW5EJU7d27Gjh1LQkICKSkpdOrUiZdfftnrWCIi2YqKXAgzM1544QVefPFFAPr27cvTTz/tcSoRkexDRS4MPPLII7z77rtERUXRsGFDr+OIiGQbeoQgTHTq1IkWLVpw2WWXpa5zzumZOhHJ0XQkF0bSFrgZM2bQuHFj9uzZ42EiERFvBbXImVkLM9tgZhvNbEAG268ws3lm9p2ZrTKzm4KZJ6dISUmhX79+fPXVVzRs2JCkpCSvI4mIeCJoRc7MIoBRQEugKtDRzKqm2+0JYIpz7u9AB+A/wcqTk0RGRvLFF19QtWpV1qxZw/XXX8/69eu9jiUikuWCeSRXB9jonNvknDsGTAJuTbePA6L984WAX4OYJ0cpU6YMCxYsoH79+mzdupUGDRqodxQRyXGCWeRKA9vSLCf516U1COhkZknATODhIObJcYoWLcqXX37JTTfdxO7du9U7iojkOBbosC1mVsA5F/BgZmbWFmjhnOvuX74bqOuc65lmn77+DMPMrD7wFnCNc+5kuve6H7gfoFixYrWnTJkSaIyQcejQIaKiooLy3ikpKQwdOpSvvvqK4cOHU6VKlaB8TkaC2S4vqV2hIxzbBGpXeo0bN17pnLv2Lxucc2edgFhgLbDVv1wD+E8Ar6sPzE6z/BjwWLp91gBl0ixvAoqf7X0rVarkwtG8efOC+v4nT550q1ev/su6YAt2u7yidoWOcGyTc2pXesAKl0HNCOR05cvAjcBuf1H8AQjkiePlQEUzK2dmefDdWDI93T5bgaYAZnY1kA/4XwDvLefJzKhWrVrq8uTJk2nVqhUHDhzwMJWISHAFdE3OObct3aoTAbwmBegJzAbW4buLco2ZPW1mrf27PQLcZ2Y/ABOBLv6KLEH0xx9/0K9fP2bOnEmDBg3YunWr15FERIIikB5PtplZLODMLDfQG1/ROifn3Ex8N5SkXfdUmvm1wPWBx5XMkC9fPubPn8/NN99MYmIiderU4dNPP+W6667zOpqISKYK5EiuB/AQvjsjtwM1/csSwsqXL8+SJUto0qQJO3fu5IYbbmDatGlexxIRyVTnLHLOud+dc3c550o454o75zo553ZnRTgJriJFivD555/TrVs3kpOTiYuLY/To0V7HEhHJNOc8XWlmY/E9tH0a59y9QUkkWSp37ty8+eabVKpUicGDB1O/fn2vI4mIZJpArsnNSDOfD7gd9UwSVsyMhIQEunbtSrFixVLXHz16lLx583qYTETk4gRyuvLDNNME4A7grw/cSchLW+DGjx9PzZo1+fnnnz1MJCJycS6kW6+KQPHMDiLZx4kTJ3j11VdZv349devWZeHChV5HEhG5IOcscmZ20MwOnPoT+BR4NPjRxCsRERHMmTOHli1bsnv3bpo2bcqECRO8jiUict4COV1Z0DkXnebPSs65D7MinHgnOjqa6dOn07NnT44dO0anTp0YNGgQelZfRELJGW88MbNaZ3uhc+7bzI8j2UlkZCQjR46kUqVK9OnTh3//+98cOXKEIUOGeB1NRCQgZ7u7cthZtjmgSSZnkWzq4Ycfpnz58vTo0YOuXbt6HUdEJGBnLHLOucZZGUSyt5tuuomffvrptEcKdu3aRfHif70HKSkpiffen8jUj2ewY+dOLi9Rgna3taLTnR2JiYnJytgiksMFdHelmV1jZneYWedTU7CDSfaTtsC99tprVKpU6S+DsC5evJiWbTow+dsdFGnRmwpdX6RIi95M/nYHLdt0YPHixVkdW0RysEDurhwIjPRPjYEhQOuzvkjCmnOO+fPns3//fm666Saee+45nHMkJSUR3yeBok26E3P97eQvUhzLFUH+IsWJuf52ijbpTnyfBJKSkrxugojkEIEcybXFN+bbDudcV3yDphYKairJ1syMiRMn8u9//xvnHI8//jjt2rXjrbFjibyqLtGlK2T4uujSFYi8qg7vvT8xixOLSE4VSJFLds6dBFLMLBrYBZQJbizJ7nLlysVTTz3Fp59+SnR0NB9++CHPP/c8+UuUO+vrLq3agKmffJZFKUUkpwukyK0ws8LAm8BK4FtgSTBDSeho1aoVy5Yto0qVKvyRfIT1M94467N0eaOLsn/fvqwLKCI52jk7aHbOPeifHW1mnwPRzrlVwY0loaRy5cosXbqUchUqUblVPGZ2xn2PHthDocKFsy6ciORogdx4Mt3M7jSzS5xzm1XgJCPR0dH07/cIyTs3Ab6bU375+kNSjh45bb/daxfS7tabvYgoIjlQIKcrhwENgLVm9oGZtTWzfEHOJSGo050dSdm0lAPbN7J5wTTWfDCc7177F4d2bQXgwPaNpGxaRqc7O3qcVERyikBOV34FfGVmEfh6ObkPeBuIDnI2CTExMTG8PnwI8X0SOFGsIpcUi+Hw/7axYOh9XFGnGZf88TuvDx+iB8JFJMsE+jB4fiAO6AFcB4wPZigJXbGxscyaNol7bqhKhavKkq9AFCeOHmHzwunc2LAe9erV8zqiiOQggVyTmwKsw3cU9ypQ3jn3cLCDSeiKiYlhQEJ/vv9mITNnTGfw4MEADBs2jDZt2nDw4EGPE4pIThHIkdxb+ApbD+fcPP8zcyIBMTMef/xxPvvsMwoXLsyWLVuIiIjwOpaI5BCBXJObnRVBJLy1bNmS5cuXExkZSYECBQA4efIkuXJdyOD0IiKB0S+MZJkKFSpw5ZVXAr5HDO666y4GDRrEyZM6OSAiwXHOIzmRYFi2bBmTJ0/GOcd3333HO++8Q6FC6hJVRDJXoHdXljazWDNreGoKdjAJb3Xr1mXWrFkULlyY6dOnU7duXdavX+91LBEJM4HcXfkCsAh4Aujvn/oFOZfkADfeeCMrVqygevXqbNiwgTp16jB58mSvY4lIGAnkdOVtQGXn3NEgZ5EcqHz58ixZsoRu3boxefJkOnToAED79u09TiYi4SCQIrcJyA2oyElQXHLJJUycOJEbbriBd999l1tvvdXrSCISJgK5JncE+N7MXjezEaemYAeTnMXMeOCBB1i4cCH58vm6Rt23bx8ffvihx8lEJJQFUuSmA88Ai/GNJ3dqEsl0p56bc87RrVs32rZty/33309ycrLHyUQkFJ2zyDnnxgMT+bO4ve9fJxJUzZs3J2/evLz55pu6+1JELkggd1c2An4CRgH/AX7UIwQSbGZGfHw833zzDRUrViQxMZFrr72WCRMmeB1NREJIoOPJNXfO3eCcawjcCLwc3FgiPjVr1mTlypV07NiRw4cP06lTJ7p3765eUkQkIIEUudzOuQ2nFpxzP+K721IkSxQsWJAJEybwxhtvkDdvXvLly6c+L0UkIIE8QrDCzMYA7/mX7wJWBC+SyF+ZGffddx+xsbGUL18+df2ePXsoWrSoh8lEJDsL5J/DDwBrgV7+aa1/nUiWq1at2mmPGFx77bXce++9HDlyxONkIpIdBXJ35VHn3EvOuTb+6WX1fiLZwYoVK/jtt98YO3YsderUYe3atV5HEpFs5oxFzj8iOGaWaGar0k9ZF1EkY//85z9ZtmwZVapUYc2aNVx33XWMGzfO61giko2c7Uiut//PVsAtGUznZGYtzGyDmW00swEZbH/ZzL73Tz+a2b7ziy85XfXq1Vm+fDl33303R44coWvXrnTp0oXDhw97HU1EsoEzFjnn3G/+2Qedc1vSTsCD53pjM4vA92xdS6Aq0NHMqqb7jH8552o652oCI4FpF9gOycGioqIYP348b7/9Nvnz52f8+PHMmzfP61gikg0EcuNJswzWtQzgdXWAjc65Tc65Y8Ak4Gw973bE17OKyHkzM7p27cry5ct5/vnnadWqldeRRCQbONs1uQfMLBGonO563C9AINfkSgPb0iwn+ddl9FllgXLA3MCji/xVtWrVePTRR1OXly1bRqtWrfj11189TCUiXjHnXMYbzAoBRYDngLTX0w465/ac843N2gItnHPd/ct3A3Wdcz0z2PdRIMY59/AZ3ut+4H6AYsWK1Z4yZcq5Pj7kHDp0iKioKK9jZDov2+Wco2fPnqxdu5bo6GgeeeQRGjbMnB7p9H2FjnBsE6hd6TVu3Hilc+7av2xwzgU0AcWBK05NAexfH5idZvkx4LEz7PsdEBtIjkqVKrlwNG/ePK8jBIXX7dq+fbtr3ry5AxzgunTp4vbv33/R7+t1u4IlHNsVjm1yTu1KD1jhMqgZgXTQfIuZ/QT8AnwFbAZmBVBYlwMVzaycmeUBOuAbtif9+1fBd8S4JID3FDkvpUqVYtasWYwYMYJ8+fIxbtw4atasyaJFi7yOJiJZIJAbT54F6gE/OufKAU2Bb871IudcCtATmA2sA6Y459aY2dNm1jrNrh2ASf5KLJLpcuXKxcMPP8zKlSupWbMmv/zyC7feeiuHDh3yOpqIBFkgfVced87tNrNcZpbLOTfPzIYH8ubOuZnAzHTrnkq3PCjQsCIXo2rVqixdupSBAwdSo0aNsLyeISKnC6TI7TOzKOBrYIKZ7QL0pK2EpDx58vDcc8+dtm7kyJFERETwwAMPYGYeJRORYAjkdOWtwBHgX8DnwM/4ekERCXnbtm2jX79+PPTQQ9x8883s2LHD60gikokCKXJPOedOOudSnHPjnXMjgEfP+SqREFCmTBkmTJhAkSJFmDVrFtWrV+fjjz/2OpaIZJJg9ngiEhLatm1LYmIizZo14/fff+f222+ne/fuHDx40OtoInKRAunxpEoGPZ4kZl1EkeArXbo0n3/+Oa+88gp58+blrbfeIj4+3utYInKRznbjyfv4noe7oB5PREJNrly56NWrF02bNuWBBx7g2Wef9TqSiFyks41CsN85txl4Bdjj/hyBIMXM6mZVQJGsVq1aNb7++muuuuoqwNcrUP/+/dmwYYPHyUTkfAVyTe41IO1Ts4f860RyhLFjx/Liiy9Ss2ZNXnzxRU6cOOF1JBEJUCBFztL2RuKcO0lgz9eJhIW4uDg6d+7MH3/8Qf/+/YkuXIQ77u5K7dgbeH7IUJKSkryOKCJnEEiR22Rmvcwst3/qDWwKdjCR7KJQoULEx8dTtuLV5M4fxZFDB/n9tyT25S/NpBXbadmmA4sXL/Y6pohkIJAi1wOIBbbjGxOuLv5hb0RygqSkJOL7JFC2TQJNBk3litjWuBMpbJo7kRNEULRJd+L7JOiITiQbOudpR+fcLnydKIvkSO+9P5HIq+oSXboCAH/r0J/CVeqxc/lMrqjfilyRuTlwVR3ee38iAxL6e5xWRNIKZKidSmY2x8xW+5f/ZmZPBD+aSPYw9eMZXFr1+tPWFSlfk+vue45ckbkBKFi2OoMHD2bZsmVeRBSRMwjkdOWb+AY8PQ7gnFuFjuwkB9m/fz/5Cl161n22LZ3JoQP7qV+/PgkJCSQnJ2dROhE5m0CKXAHnXPp/nqYEI4xIdlSoUCH+2L/7rPtc+Y82FC8VA8DQoUOpUaMGCxcuzIp4InIWgRS5382sPOAAzKwt8FtQU4lkI+1ua8XutWcfSXzfT8v5V+9eLFmyhGrVqvHTTz/RsGFDevXqpcFZRTwUSJF7CHgdXx+W24E++O64FMkROt3ZkZRNSzmwfWOG2w9s30jKpmV0urMjderUYeXKlTz55JNEREQwcuRIvv/++6wNLCKpArm7chPwTzO7BMjlnFPX7JKjxMTE8PrwIcT3SeDAVXW4tGoDnOUhee8udq9dSMqmZbw+fAgxMb7TlXnz5uXpp5+mTZs2zJ07lwYNGqS+19GjR8mbN69XTRHJcQK5u/JSMxsBLADmm9krZnb2q/AiYSY2NpZZ0ybRvlZJ9s4ewc/j+rF39gja1yrJrGmTiI2N/ctratasSd++fVOX58yZQ4UKFfjss8+yMrpIjhbI6cpJwP+AOKCtf35yMEOJZEcxMTEMSOjPykXzee+t0axcNJ8BCf1Tj+DO5c033yQpKYlWrVrRuXNndu8++80sInLxAilyJZ1zzzjnfvFPzwIlgh1MJNxMmDCBYcOGkS9fPt59912qVKnC+PHjSdM1rIhkskCK3H/NrIOZ5fJPdwCzgx1MJNxERETQt29fVq1aRaNGjfj999/p0qULjRs3ZuvWrV7HEwlLgRS5+/ANoHrUP00C4s3soJkdCGY4kXBUsWJF5s6dyzvvvEOxYsVYv3490dHRXscSCUvnLHLOuYLOuVzOudz+KZd/XUHnnP7PFLkAZsbdd9/N+vXr+fjjjylcuDAAycnJzJ6tEyUimSWQuyu7pVuOMLOBwYskknMULVqUevXqpS4/99xztGjRgjZt2rBt2zYPk4mEh0BOVzY1s5lmVtLMrgG+AQoGOZdIjlS8eHGioqL46KOPuPrqqxk2bBjHjx/3OpZIyArkdOWdwHggEfgM6OOc6xfsYCI5Uc+ePVm3bh1xcXEcPnyYfv36ce2117JkyRKvo4mEpEBOV1YEegMfAluAu82sQLCDieRUMTExfPDBB3z22WeUK1eOVatWERsby4oVK7yOJhJyztmtF/Ap8JBzbo6ZGdAXWA5UC2oykRzupptuYvXq1QwePJjExERq167tdSSRkBPINbk6zrk5AM5nGHB7cGOJCECBAgUYPHgwn3zyCb5/Y8LatWtp1qwZa9eu9TidSPZ3xiJnZgkAzrkDZtYu3eYuwQwlIqc7VeAAnnjiCb788ktq1KjBY489xpEjRzxMJpK9ne1ILu3o34+l29YiCFlEJABjxowhPj6elJQUnn/+eapVq8aMGTO8jiWSLZ2tyNkZ5jNaFpEsUrRoUUaPHs2SJUuoUaMGmzdv5pZbbqF169Zs2bLF63gi2crZipw7w3xGyyKSxerVq8eKFSt46aWXiIqK4osvvjjttKaInL3I1TCzA2Z2EPibf/7UcvUsyiciZxEZGcm//vUvfvrpJyZPnswVV1wBwMmTJ5k6dSonTpzwOKGIt85Y5JxzEc65aH8flZH++VPLubMypIic3eWXX07r1q1Tl8eNG8cdd9zBddddx4IFCzxMJuKtQB4hEJEQc+mllxITE8N3331Hw4YNad++va7XSY6kIicShm699VY2bNjAwIEDyZ8/P1OmTKFKlSo89dRTHD582Ot4IllGRU4kTBUoUIBBgwaxfv16OnTowB9//MEzzzzDO++843U0kSyjIicS5q644gomTpzIggUL6NChA927d0/dtnv3bg+TiQRfIH1XXjAzawG8AkQAY5xzz2ewzx3AIHyPJfzgH/VARDJZgwYNaNCgQeryzp07qVy5Ms2bN6dCpcrMnvs1O3bu5PISJWh3Wys63dmRmJgYDxOLXLygFTkziwBGAc2AJGC5mU13zq1Ns09FfL2pXO+c22tmxYOVR0ROt2jRIo4cOcLUqVPJFZmbco3u4KpOvckbAZO/XcSEDzrw+vAhxMbGeh1V5IIF83RlHWCjc26Tc+4YMAm4Nd0+9wGjnHN7AZxzu4KYR0TSqFOnDldWqc6lFf/OyZTj/PzlBFa8+hD7tq6ndOxtFG3Snfg+CSQlJXkdVeSCBbPIlQa2pVlO8q9LqxJQycwWmdk3/tObIpIF3nt/IpdUbUT9h0dQ96GXKViyHH/s2cHKt/4faz8aSXTpCkReVYf33p/odVSRCxbUa3IBfn5FoBEQA3xtZtWdc/vS7mRm9wP3AxQrVoz58+dnbcoscOjQIbUrhIRDu8a+N4miLftw8OBB8pWqzN8fGMH2ZZ+xZe77FK5S37e+7N95+72XqVfnOq/jXrBw+K4yonYFJphFbjtQJs1yjH9dWknAUufcceAXM/sRX9FbnnYn59wbwBsAlStXdo0aNQpWZs/Mnz8ftSt0hEO7Tpw4waWlrsByRaSus/qtqdQojog8+QC4pEB+Erds4qOPPuKJJ56gWLFiXsW9YOHwXWVE7QpMME9XLgcqmlk5M8uDb+ie6en2+RjfURxmdhm+05ebgphJRPwKFSrEH/v/+gjBqQIHsPeX1Rw+sJ8RI0ZQvnx5nn32WT1MLiElaEXOOZcC9ARmA+uAKc65NWb2tJmd6mRvNrDbzNYC84D+zjk9uCOSBdrd1ordaxeddZ/knZvo1edftGzZkoMHD/Lkk09SoUIFXn/9dY4fP55FSUUuXFAfBnfOzXTOVXLOlXfODfave8o5N90/75xzfZ1zVZ1z1Z1zk4KZR0T+1OnOjqRsWsqB7Rsz3H5g+0ZSNi2j/yN9mTlzJnPnzuW6665jx44d9OjRg9q1a6vQSbanHk9EcqiYmBheHz6EPXPHkLRoGsl7d+FOnCB57y6SFk1jz9wxvD58SOoD4Y0bN2bp0qVMnTqVihUr0qRJE3Ln9g1I4pyGmJTsSUVOJAeLjY1l1rRJtK9Vkr2zR/DzuH7snT2C9rVKMmvapL88CG5mtG3bljVr1vDss8+mrp88eTKtWrUiMTExq5sgclYqciI5XExMDAMS+rNy0Xzee2s0KxfNZ0BC/7N26ZU7d26ioqIA31HcCy+8wGeffUaNGjXo2rUrW7duzar4ImelIiciF8XMmD17Nj179iQiIoJx48ZRqVIlEhIS2LNnj9fxJIdTkRORi1a8eHFGjhyZOqzP0aNHGTp0KOXLl2f58uXnfgORIFGRE5FMU758eSZOnMiKFSto2rQpBQsWpHr16qnbdYOKZDUVORHJdLVr1+bLL79k+fLl5Mvne7h8z5491KhRgzfffFOPHkiWUZETkaApUaJE6vy4ceNITEzk/vvvp1KlSrz99tsqdhJ0KnIikiX69OnDpEmTqFKlCps3b6Zbt25cffXVjB8/npSUFK/jSZhSkRORLJErVy7at2/P6tWree+996hUqRI///wzXbp04Z577vE6noQpFTkRyVIRERHcddddrFmzhvHjx1O+fHnuvffe1O379u3jxIkTHiaUcKIiJyKeiIyMpHPnzqxfv54mTZqkro+Pj+dvf/sbU6dO5eTJkx4mlHCgIicinoqMjMTMAN+AmcuWLWPt2rXccccd1KxZk2nTpqnYyQVTkRORbCMqKooNGzYwevRoYmJiSExMJC4ujlq1avHJJ5/oOTs5bypyIpKt5MmTh/j4eDZu3MioUaMoVaoUP/zwA3FxceoTU86bipyIZEt58+blwQcf5Oeff2bEiBH079+fsmXLAnDy5ElmzZql05hyTpFeBxAROZt8+fLx8MMPn7bu008/5bbbbuPqq69mwIABdOzYMXVsO5G0dCQnIiEnJSWFMmXKsG7dOu655x4qVqzIqFGjSE5O9jqaZDMqciIScuLi4ti4cSNjx46lcuXKbNmyhZ49e3LllVcyevRor+NJNqIiJyIhKU+ePHTp0oU1a9bwwQcfULt2bXbt2sX+/fu9jibZiK7JiUhIi4iIIC4ujjZt2vDFF19Qr1691G0vv/wyCxYsoFy5cqk3rUjOoiM5EQkLZkbz5s2Jjo4G4Pjx4zz//PN89NFHVKhQgS5durBu3TqPU0pWU5ETkbCUO3duvvzyS5o2bcrJkycZP3481apVIy4ujhUrVngdT7KIipyIhK3q1avzxBNP8OOPPxIfH0/u3LmZNm0a1113HYsXL/Y6nmQBFTkRCXvly5dn9OjR/PLLL/Tr14/rr7+e+vXrp25fsWKFRj4IUypyIpJjlCpViqFDh7JgwYLUTqE3bdpEvXr1qFKlCqNGjeLw4cMep5TMpCInIjnOqQIHsHnzZmJiYti4cSM9e/akTJkyDBgwgKSkJA8TSmZRkRORHK1JkyZs3LiRqVOnUr9+ffbu3csLL7xAuXLl6NKli0Y+CHEqciKS40VGRtK2bVsWL17MkiVLuOOOOzh58iRHjx5NPeo7efKkrtuFIBU5EZE06tWrx+TJk/n5558ZPHhw6vrPP/9c1+1CkIqciEgGrrzySq666qrU5ffff1/X7UKQipyISADGjRvH1KlTiY2NPe26XadOnfj++++9jidnoCInIhKAU9ftFi1adNp1uwkTJjBz5kyv48kZqINmEZHzdOq63ZYtW3j11VeJj49P3fb2229z6NAh7rnnHgoVKuRhSgEdyYmIXLCyZcsydOhQLr30UgCOHTvGk08+Se/evSldujQ9evRg1apVHqfM2VTkREQySUREBCNGjKBx48YcPnyY119/nRo1avCPf/yDSZMmcezYMa8j5jgqciIimeTU2HZz585lzZo19OzZk4IFC7Jw4UI6duzIvHnzvI6Y46jIiYgEQdWqVRk5ciTbt2/ntdde45ZbbqFZs2ap24cNG8acOXPUo0qQ6cYTEZEgKliwID169KBHjx6p67Zv386jjz7KiRMnqFy5Mg8++GBAN6okJSXx3vsTmfrxDHbs3MnlJUrQ7rZWdLqzIzExMcFuSkjSkZyISBbLnz8/AwcOpHTp0mzYsCGgG1UWL15MyzYdmPztDoq06E2Fri9SpEVvJn+7g5ZtOmh8vDMIapEzsxZmtsHMNprZgAy2dzGz/5nZ9/6pezDziIhkB0WLFuXJJ5/kl19+4YMPPqBJkyapN6rUqlWLnTt3nrZ/UlIS8X0SKNqkOzHX307+IsWxXBHkL1KcmOtvp2iT7sT3SVAPLBkI2ulKM4sARgHNgCRguZlNd86tTbfrZOdcz2DlEBHJrnLnzk1cXBxxcXGsW7eO1157jQMHDlCiRAnA1yn0sGHD2L1nL5FX1SW6dIUM3ye6dAUOXFWH996fyICE/lnZhGwvmNfk6gAbnXObAMxsEnArkL7IiYjkeFdffTUjRow4bd2cOXNISEgAoGiFmuQqUJjLqzcgV2Tuv7z+0qoNmPrJCBW5dIJ5urI0sC3NcpJ/XXpxZrbKzD4wszJBzCMiElIuv/xyOnbsCGbs2fg93459ii+fvJ21H4/i8K6tp+2bN7oo+/ft8yZoNmbBun3VzNoCLZxz3f3LdwN1056aNLNLgUPOuaNmFg+0d841yeC97gfuByhWrFjtKVOmBCWzlw4dOkRUVJTXMTKd2hVawrFd4dCm7g/2JqXU3/h99QIO79wCQJ7oy6jXfxyWy3escnT/7+z5fDivvzLMy6gX7UK/r8aNG690zl2bfn0wT1duB9IemcX416Vyzu1OszgGGJLRGznn3gDeAKhcubJr1KhRpgbNDubPn4/aFTrUrtARDm3q3qUTk7/dwTWPv8u+revYtmQGEQUvI9r/yMHRg3tZ/eFQ2rdsxA033JA60GsoyuzvK5inK5cDFc2snJnlAToA09PuYGYl0yy2BtYFMY+ISEjqdGdHUjYt5eCvP1OkbFX+1iGBKxq2S93+85fvs/vnRP7z6khq1qzJyJEj2bNnj4eJs4+gFTnnXArQE5iNr3hNcc6tMbOnzay1f7deZrbGzH4AegFdgpVHRCRUxcTE8PrwIeyZO4akRdNI3rsLd+IEyXt3kbRoGhH/20CHDh247LLLWLVqFb169aJUqVLceeedOb4rsaD2eOKcmwnMTLfuqTTzjwGPBTODiEg4iI2NZda0Sb4eTz4Zwc4dOyhx+eW0v/VmOr04nZiYGI4ePcr06dMZM2YMX3zxBRMnTmTv3r00btzY6/ieUbdeIiIhIiYmhgEJ/RmQ0D/Da1d58+alXbt2tGvXji1btjB27Fjq1auXun3hwoW88MILdOnShVatWpE3b94sbkHWU5ETEQlDZcuWZdCgQaetGzNmDDNmzGDGjBkUKVKE9u3b07lzZ+rVqxfSN6ucjfquFBHJIYYMGcKwYcOoWbMme/fuZfTo0cTGxlKpUiVGjRrldbygUJETEckhihcvTt++ffnuu+9YtWoV/fv3p2TJkmzcuJFt2/7su+PgwYPs3bvXw6SZR0VORCQHql69OkOGDGHbtm3Mnj2b+Pj41G3jxo2jZMmStGvXjk8//ZTjx497mPTiqMiJiORgERERNG/enHLlyqWu++mnnzh27BgffPABrVu3plSpUvTq1Yvly5eH3CCvKnIiInKaESNGsHXrVp5//nmqVq3K77//zsiRI6lTpw7333+/1/HOi4qciIj8RUxMDI8++iirV69m5cqV9O7dm+LFi5/22ML333/P2LFj2b9/v3dBz0FFTkREzsjMqFWrFsOHDycpKYk77rgjdduoUaO49957KV68OLfddhuTJk3i8OHDHqb9KxU5EREJSO7cucmd+8+x7Bo3bkyjRo04fvw4n3zyCR07dqR48eJ06NCBL774wsOkf1KRExGRC3Kqb8ykpCReeeUV6tevz5EjR5g8efJpRS45OZljx46d8X2SkpJ4fshQasfewF3delA79gaeHzKUpKSki86oIiciIhfl1N2XixcvZvPmzQwZMoR77rkndftbb71FyZIlue+++5gzZw4nTpxI3bZ48WJatunA5G93UKRFbyp0fZEiLXoz+dsdtGzTgcWLF19UNnXrJSIimaZs2bL079//tHXLly9nz549jBkzhjFjxlCiRAnatWtH48aNeeq5YVza9D6iS1cAIOXgQfIXKk7M9bdz4MrqxPdJYNa0ScTExFxQHh3JiYhIUI0fP57ExESeeOIJKlSowM6dO3n11VeJi4vj1//tSS1w6UWXrkDkVXV47/2JF/zZKnIiIhJ011xzDc888ww//vgjK1eupH///uTJk5eStf6Zus/ezWvY89PK0153adUGTP3kswv+XBU5ERHJMqceSRgyZAgx5StTtsFtqduSls4iT1SR0/bPG12U/fv2XfDnqciJiIgnChcuzLFDfz5IXqbezVxyebnT9jl6YA+FChe+4M9QkRMREU+0u60Vu9cuSl0uXPbqv4xrt3vtQtrdevMFf4aKnIiIeKLTnR1J2bSUA9s3Zrj9wPaNpGxaRqc7O17wZ6jIiYiIJ2JiYnh9+BD2zB1D0qJpJO/dhTtxguS9u0haNI09c8fw+vAhF/z4AOg5ORER8VBsbCyzpk3ivfcnMvWTEezcsYMSl19O+1tvptOLF/583CkqciIi4qmYmBgGJPRnQEJ/5s+ff9pIBxdLpytFRCRsqciJiEjYUpETEZGwpSInIiJhS0VORETCloqciIiELXPOeZ3hvJjZQWCD1zmC4DLgd69DBIHaFVrCsV3h2CZQu9Ir65wrln5lKD4nt8E5d63XITKbma1Qu0KH2hU6wrFNoHYFSqcrRUQkbKnIiYhI2ArFIveG1wGCRO0KLWpX6AjHNoHaFZCQu/FEREQkUKF4JCciIhKQbFvkzKyFmW0ws41mNiCD7XnNbLJ/+1Izu9KDmOctgHZ1MbP/mdn3/qm7FznPh5m9bWa7zGz1GbabmY3wt3mVmdXK6owXIoB2NTKz/Wm+q6eyOuP5MrMyZjbPzNaa2Roz653BPiH3fQXYrlD8vvKZ2TIz+8Hfrn9nsE/I/RYG2K7M+S10zmW7CYgAfgauAvIAPwBV0+3zIDDaP98BmOx17kxqVxfgVa+znme7GgK1gNVn2H4TMAswoB6w1OvMmdSuRsAMr3OeZ5tKArX88wWBHzP4bzDkvq8A2xWK35cBUf753MBSoF66fULxtzCQdmXKb2F2PZKrA2x0zm1yzh0DJgG3ptvnVmC8f/4DoKmZWRZmvBCBtCvkOOe+BvacZZdbgXeczzdAYTMrmTXpLlwA7Qo5zrnfnHPf+ucPAuuA0ul2C7nvK8B2hRz/d3DIv5jbP6W/kSLkfgsDbFemyK5FrjSwLc1yEn/9DzZ1H+dcCrAfuDRL0l24QNoFEOc/TfSBmZXJmmhBFWi7Q1F9/ymXWWZWzesw58N/Wuvv+P4VnVZIf19naReE4PdlZhFm9j2wC/jCOXfG7yuEfgsDaRdkwm9hdi1yOdmnwJXOub8BX/Dnv9Ak+/kWX1dCNYCRwMfexgmcmUUBHwJ9nHMHvM6TWc7RrpD8vpxzJ5xzNYEYoI6ZXeNxpEwRQLsy5bcwuxa57UDaqh3jX5fhPmYWCRQCdmdJugt3znY553Y75476F8cAtbMoWzAF8n2GHOfcgVOnXJxzM4HcZnaZx7HOycxy4ysEE5xz0zLYJSS/r3O1K1S/r1Occ/uAeUCLdJtC8bcw1ZnalVm/hdm1yC0HKppZOTPLg+9i6vR0+0wH7vHPtwXmOv/VymzsnO1Kd+2jNb5rC6FuOtDZf9dePWC/c+43r0NdLDO7/NS1DzOrg+//p2z94+LP+xawzjn30hl2C7nvK5B2hej3VczMCvvn8wPNgPXpdgu538JA2pVZv4XZsoNm51yKmfUEZuO7I/Ft59waM3saWOGcm47vP+h3zWwjvpsDOniXODABtquXmbUGUvC1q4tngQNkZhPx3bl2mZklAQPxXUjGOTcamInvjr2NwBGgqzdJz08A7WoLPGBmKUAy0CG7/7gA1wN3A4n+6yEAjwNXQEh/X4G0KxS/r5LAeDOLwFeUpzjnZoT6byGBtStTfgvV44mIiISt7Hq6UkRE5KKpyImISNhSkRMRkbClIiciImFLRU5ERMKWipxIGmbmzGxYmuV+ZjboHK/pYWadgx7uHPy97M84z9cM9fcCP/QCPu/x832NSFZTkRM53VGgzfn0hOGcG+2ceyeImYLpfuBvzrn+F/Da8ypy/ofL9ZsjWUr/wYmcLgV4A/hX+g1mdqWZzfV3GDvHzK7wrx9kZv38873MN6bZKjOb5F93ifnGpltmZt+Z2V9GnjCzkmb2tfnGzVptZv/wr29hZt/6OxWe419Xx8yW+N9rsZlVzuD9AvnM6UAUsNLM2pvZLeYbj+w7M/vSzEr494sys7FmluhvV5yZPQ/k9+ed4N+vrz/7ajPrk+bvbIOZvQOs5vTuwkSC72LH6tGkKZwm4BAQDWzG1wdgP2CQf9unwD3++XuBj/3zg4B+/vlfgbz++cL+P/8P6HRqHb6xzi5J97mPAP/PPx+Bb0y0Yvh6ly/nX1/U/2c0EOmf/yfwoX++Ef7x0gL5zFPtTTNfhD87iOgODPPPvwAMT7tfBq+tDSQCl+ArnGvwjQRwJXCSdGOFadKUVVO27NZLxEvOuQP+I49e+Lp/OqU+0MY//y4wJIOXrwImmNnH/NnLfXOg9amjPSAfvu6m0vbFtxx423ydDH/snPvezBoBXzvnfvHnOjW2XSF8XSJVxDcGV+4McgTymenFAJP9fQbmAX7xr/8nabqKcs7tzeC1DYCPnHOHAcxsGvAPfP0qbnG+celEspxOV4pkbDjQDd+Ryfm4GRiFb0Tx5ebrFd6AOOdcTf90hXPutGLjfAO0NsTXo/y4c9zI8gwwzzl3DXALvgKW3jk/MwMj8Y3EXB2IP8P7XojDmfQ+IudNRU4kA/6jpin4Ct0pi/nziOYuYEHa1/hvqijjnJsHPIrviCsKX4fcD6fpAf/v6T/PzMoCO51zb+IbVqQW8A3Q0MzK+fcp6t+9EH8OfdPlDE0452dmIO373pNm/RfAQ2myFvHPHvcfeYLv7+I2MytgZpcAt5Pu70fECypyImc2DEh7l+XDQFczW4Wvx/ve6faPAN4zs0TgO2CE842V9Qy+U4qrzGyNfzm9RsAPZvYd0B54xTn3P3x3P04zsx+Ayf59hwDP+fc90yWHQD4zvUHAVDNbCfyeZv2zQBH/DSU/AI3969/wv/8E59y3wDhgGb4Rucc4574L4DNFgkqjEIiISNjSkZyIiIQtFTkREQlbKnIiIhK2VORERCRsqciJiEjYUpETEZGwpSInIiJhS0VORETC1v8HbrULpPuaGisAAAAASUVORK5CYII=",
      "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": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.75907, 0.66326, 0.59669, 0.5017 , 0.45726])"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "factory.get_expectation_values()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1. , 1.5, 2. , 3. , 3.5])"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "factory.get_scale_factors()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.963741603277012"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "factory.get_zero_noise_limit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 3: Applying probabilistic error cancellation (PEC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mitiq import pec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Build quasi-probability representation for CNOT gate\n",
    "\n",
    "For simplicity, we assume that:\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": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.008339674869236618"
      ]
     },
     "execution_count": 64,
     "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 # a very rough approximation of the error acting on a CNOT gate\n",
    "\n",
    "# print(qiskit.transpile(circuit, backend=noisy_backend))"
   ]
  },
  {
   "cell_type": "markdown",
   "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",
    "$$\\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}.$$\n",
    "\n",
    "Key idea: Represent a unitary channel by a linear combination of noisy channels which a computer can actually implement.\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Quasi-probability representation with one-norm 1.01683 generated.\n",
      "q_0: ───@───\n",
      "        │\n",
      "q_1: ───X─── = \n",
      "\n",
      "1.008\n",
      "q_0: ───@───\n",
      "        │\n",
      "q_1: ───X───\n",
      "\n",
      "-0.001\n",
      "q_0: ───@───────\n",
      "        │\n",
      "q_1: ───X───X───\n",
      "\n",
      "-0.001\n",
      "q_0: ───@───────\n",
      "        │\n",
      "q_1: ───X───Y───\n",
      "\n",
      "-0.001\n",
      "q_0: ───@───────\n",
      "        │\n",
      "q_1: ───X───Z───\n",
      "\n",
      "-0.001\n",
      "q_0: ───@───X───\n",
      "        │\n",
      "q_1: ───X───────\n",
      "\n",
      "-0.001\n",
      "q_0: ───@───Y───\n",
      "        │\n",
      "q_1: ───X───────\n",
      "\n",
      "-0.001\n",
      "q_0: ───@───Z───\n",
      "        │\n",
      "q_1: ───X───────\n",
      "\n",
      "-0.001\n",
      "q_0: ───@───X───\n",
      "        │\n",
      "q_1: ───X───X───\n",
      "\n",
      "-0.001\n",
      "q_0: ───@───Y───\n",
      "        │\n",
      "q_1: ───X───X───\n",
      "\n",
      "-0.001\n",
      "q_0: ───@───Z───\n",
      "        │\n",
      "q_1: ───X───X───\n",
      "\n",
      "-0.001\n",
      "q_0: ───@───X───\n",
      "        │\n",
      "q_1: ───X───Y───\n",
      "\n",
      "-0.001\n",
      "q_0: ───@───Y───\n",
      "        │\n",
      "q_1: ───X───Y───\n",
      "\n",
      "-0.001\n",
      "q_0: ───@───Z───\n",
      "        │\n",
      "q_1: ───X───Y───\n",
      "\n",
      "-0.001\n",
      "q_0: ───@───X───\n",
      "        │\n",
      "q_1: ───X───Z───\n",
      "\n",
      "-0.001\n",
      "q_0: ───@───Y───\n",
      "        │\n",
      "q_1: ───X───Z───\n",
      "\n",
      "-0.001\n",
      "q_0: ───@───Z───\n",
      "        │\n",
      "q_1: ───X───Z───\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) # full quasi-probability representation of the CNOT gate\n",
    "\n",
    "import mitiq.pec.representations.optimal.find_optimal_representation \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run PEC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 214 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 214 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 213 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 213 gates using 200 shots.\n",
      "Executing circuit with 213 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 212 gates using 200 shots.\n",
      "Executing circuit with 213 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 215 gates using 200 shots.\n",
      "Executing circuit with 212 gates using 200 shots.\n",
      "Executing circuit with 214 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 213 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 213 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 214 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 213 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 212 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 213 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 213 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 213 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Estimation error without Mitiq   : 0.24168\n",
      "Estimation error with Mitiq (PEC): 0.008878595141676104\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": "markdown",
   "metadata": {},
   "source": [
    "From above, we see that most of the times, we execute the original circuit, but sometimes, we apply one of those noisy operations in the previous section. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PEC reduced the estimation error by 96.3%.\n"
     ]
    }
   ],
   "source": [
    "print(f\"PEC reduced the estimation error by {(unmitigated_error - pec_error) / unmitigated_error :.1%}.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even under very strong assumptions on the noise model, we are able to get a good improvement."
   ]
  },
  {
   "cell_type": "markdown",
   "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). 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).\n",
    "\n",
    "We try to do something like PEC without full knowledge of the noise model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "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": "markdown",
   "metadata": {},
   "source": [
    "Similar to ZNE, we take a circuit, scale it and then take a linear combination of the results (ideal gate as a linear combination of channels). But we are doing this locally, at the level of an individual gate.\n",
    "\n",
    "This will work for different levels and types of noise, based on the extrapolation idea. This method is independent from the type of noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NEPEC quasi-probability representation with one-norm 1.0202 generated.\n",
      "q_0: ───@───\n",
      "        │\n",
      "q_1: ───X─── = \n",
      "\n",
      "1.010\n",
      "q_0: ───@───\n",
      "        │\n",
      "q_1: ───X───\n",
      "\n",
      "-0.010\n",
      "q_0: ───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───@───\n",
      "        │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │   │\n",
      "q_1: ───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───X───\n"
     ]
    }
   ],
   "source": [
    "# scale_factors = [1, 3, 5]  \n",
    "scale_factors = [1, 100] # Large scale factors reduce the PEC sampling cost.\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",
   "metadata": {},
   "source": [
    "### Run NEPEC (noise-extended PEC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 407 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 309 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 407 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 309 gates using 200 shots.\n",
      "Executing circuit with 505 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 309 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 309 gates using 200 shots.\n",
      "Executing circuit with 309 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 309 gates using 200 shots.\n",
      "Executing circuit with 309 gates using 200 shots.\n",
      "Executing circuit with 309 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 407 gates using 200 shots.\n",
      "Executing circuit with 309 gates using 200 shots.\n",
      "Executing circuit with 407 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 309 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 309 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 407 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 309 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 309 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 309 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 309 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Executing circuit with 309 gates using 200 shots.\n",
      "Executing circuit with 211 gates using 200 shots.\n",
      "Estimation error without Mitiq     : 0.24168\n",
      "Estimation error with Mitiq (NEPEC): 0.0004337331703063718\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": 99,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NEPEC reduced the estimation error by 99.8%.\n"
     ]
    }
   ],
   "source": [
    "print(f\"NEPEC reduced the estimation error by {(unmitigated_error - nepec_error) / unmitigated_error :.1%}.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.12 ('qiskit-env')",
   "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.8.12"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "c5e597ca477812040ed5a1629c6451aa891ac46335823a5f6165f6c73dfc5551"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}