{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "576260b7",
   "metadata": {},
   "source": [
    "# Variational Quantum Eigensolver improved with Zero Noise Extrapolation \n",
    "\n",
    "In this example we investigate how Zero Noise Extrapolation (ZNE) can improve \n",
    "convergence when applied to a variational problem. ZNE works by computing the \n",
    "observable of interest at increased noise levels, i.e. beyond the minimum noise \n",
    "strength in the computer, and then extrapolating back to the zero-noise limit. \n",
    "The two main components of ZNE are noise scaling and extrapolation. You can read more about \n",
    "ZNE in the Mitiq Users Guide.\n",
    "\n",
    "\n",
    "The Variational Quantum Eigensolver (VQE) is a hybrid quantum-classical algorithm used to\n",
    "solve eigenvalue and optimization problems. The VQE algorithm consists of a quantum \n",
    "subroutine run inside of a classical optimization loop. In this example, the goal of the \n",
    "optimization is to find the smallest eigenvalue of a matrix H, which is the Hamiltonian \n",
    "of a simple quantum system. The quantum subroutine prepares the quantum state \n",
    "|Ψ(vec(θ))⟩ and measures the expectation value ⟨Ψ(vec(θ))|H|Ψ(vec(θ))⟩. By the \n",
    "variational principle, ⟨Ψ(vec(θ))|H|Ψ(vec(θ))⟩ is always greater than the smallest \n",
    "eigenvalue of H, which means a classical optimization loop can be used to find this \n",
    "eigenvalue. \n",
    "\n",
    "\n",
    "The VQE example shown here is adapted from the `VQE` function in Grove \n",
    "[[1]](#references) and the pyQuil / Grove VQE tutorial [[2]](#references). \n",
    "\n",
    "\n",
    "## Defining the quantum system using pyQuil"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "826a27ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from pyquil import get_qc, Program\n",
    "from pyquil.gates import RX, RY, S, T, Z, CNOT, MEASURE\n",
    "from pyquil.paulis import PauliTerm, PauliSum, sZ\n",
    "from pyquil.noise import pauli_kraus_map, append_kraus_to_gate\n",
    "from typing import List, Union\n",
    "from collections import Counter\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy import optimize\n",
    "import mitiq\n",
    "from mitiq import zne\n",
    "from mitiq.zne.scaling.folding import fold_gates_from_right"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b376fe37",
   "metadata": {},
   "source": [
    "Use the `get_qc` command to initialize the simulated backend \n",
    "where the pyQuil program will run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b6edbe45",
   "metadata": {},
   "outputs": [],
   "source": [
    "backend = get_qc(\"2q-qvm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "477be8fb",
   "metadata": {},
   "source": [
    "Define example ansatz, consisting of a rotation by angle `theta` and a layer of static gates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "df13534d",
   "metadata": {},
   "outputs": [],
   "source": [
    "program = Program()\n",
    "theta = program.declare(\"theta\", memory_type=\"REAL\")\n",
    "program += RX(theta, 0)\n",
    "program += T(0)\n",
    "program += CNOT(1, 0)\n",
    "program += S(0)\n",
    "program += Z(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "deb0c6bc",
   "metadata": {},
   "source": [
    "Simulate depolarizing noise on the static gates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8bfd9336",
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_noise_to_circuit(quil_prog):\n",
    "    \"\"\"Define pyQuil gates with a custom noise model via Kraus operators:\n",
    "        1. Generate Kraus operators at given survival probability\n",
    "        2. Append Kraus operators to the gate matrices\n",
    "        3. Add custom gates to circuit\n",
    "\n",
    "    Args:\n",
    "      quil_prog: the pyQuil quantum program to which the noise model will be added\n",
    "      \n",
    "    Returns:\n",
    "      A quantum program with depolarizing noise on the static gates.\n",
    "    \"\"\"\n",
    "    prob = 0.8\n",
    "    num_qubits = 1\n",
    "    d = 4 ** num_qubits\n",
    "    d_sq = d ** 2\n",
    "    \n",
    "    kraus_list = [(1 - prob) / d] * d\n",
    "    kraus_list[0] += prob\n",
    "    kraus_ops = pauli_kraus_map(kraus_list)\n",
    "    \n",
    "    k_list = [(1 - prob) / d_sq] * d_sq\n",
    "    k_list[0] += prob\n",
    "    k_ops = pauli_kraus_map(k_list)\n",
    "    \n",
    "    T_gate = np.array([[1, 0], [0, np.exp(1j * np.pi / 4)]])\n",
    "    CNOT_gate = np.block(\n",
    "        [[np.eye(2), np.zeros((2, 2))], [np.zeros((2, 2)), np.flip(np.eye(2), 1)]]\n",
    "    )\n",
    "    S_gate = np.array([[1, 0], [0, 1j]])\n",
    "    Z_gate = np.array([[1, 0], [0, -1]])\n",
    "    \n",
    "    quil_prog.define_noisy_gate(\"T\", [0], append_kraus_to_gate(kraus_ops, T_gate))\n",
    "    quil_prog.define_noisy_gate(\"CNOT\", [1, 0], append_kraus_to_gate(k_ops, CNOT_gate))\n",
    "    quil_prog.define_noisy_gate(\"S\", [0], append_kraus_to_gate(kraus_ops, S_gate))\n",
    "    quil_prog.define_noisy_gate(\"Z\", [0], append_kraus_to_gate(kraus_ops, Z_gate))\n",
    "    \n",
    "    return quil_prog"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "265a65ba",
   "metadata": {},
   "source": [
    "## Set up VQE: define Hamiltonian and energy expectation functions\n",
    "\n",
    "\n",
    " Hamiltonian in this example is just `sigma_z` on the zeroth qubit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d27e21fc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pyquil.quil.Program at 0x7fa796f19f90>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hamiltonian = sZ(0)\n",
    "pauli_sum = PauliSum([hamiltonian])\n",
    "\n",
    "for j, term in enumerate(pauli_sum.terms):\n",
    "    meas_basis_change = Program()\n",
    "    marked_qubits = []\n",
    "    for index, gate in term:\n",
    "        marked_qubits.append(index)\n",
    "        if gate == \"X\":\n",
    "            meas_basis_change.inst(RY(-np.pi / 2, index))\n",
    "        elif gate == \"Y\":\n",
    "            meas_basis_change.inst(RX(np.pi / 2, index))\n",
    "        program += meas_basis_change\n",
    "    readout_qubit = program.declare(\"ro\", \"BIT\", max(marked_qubits) + 1)\n",
    "\n",
    "samples = 3000\n",
    "program.wrap_in_numshots_loop(samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d510176",
   "metadata": {},
   "source": [
    "Compute expectation value of the Hamiltonian over the over the distribution \n",
    "generated from the quantum program. The following function is a modified version \n",
    "of `expectation` from the `VQE` function in Grove [[1]](#references). Here the noisy \n",
    "gates are defined inside the executor function via `add_noise_to_circuit`, since \n",
    "pyQuil custom gates cannot be folded in Mitiq."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "790a342c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def executor(\n",
    "    theta,\n",
    "    backend,\n",
    "    readout_qubit,\n",
    "    samples: int,\n",
    "    pauli_sum: Union[PauliSum, PauliTerm, np.ndarray],\n",
    "    pyquil_prog: Program,\n",
    ") -> float:\n",
    "    \"\"\"\n",
    "    Compute the expectation value of pauli_sum over the distribution generated from \n",
    "    pyquil_prog.\n",
    "    \"\"\"\n",
    "    noisy = pyquil_prog.copy()\n",
    "    noisy += [\n",
    "        MEASURE(qubit, r) for qubit, r in zip(list(range(max(marked_qubits) + 1)), readout_qubit)\n",
    "    ]\n",
    "    noisy = add_noise_to_circuit(noisy)\n",
    "    expectation = 0.0\n",
    "    pauli_sum = PauliSum([pauli_sum])\n",
    "    for j, term in enumerate(pauli_sum.terms):\n",
    "        qubits_to_measure = []\n",
    "        for index, gate in term:\n",
    "            qubits_to_measure.append(index)\n",
    "            meas_outcome = expectation_from_sampling(\n",
    "                theta, noisy, qubits_to_measure, backend, samples\n",
    "            )\n",
    "            expectation += term.coefficient * meas_outcome\n",
    "    return expectation.real"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06b44329",
   "metadata": {},
   "source": [
    "The following function is a modified version of `expectation_from_sampling` \n",
    "from the `VQE` function in Grove [[1]](#references). It is modified to follow pyQuil \n",
    "conventions for defining custom gates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ac3330cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "def expectation_from_sampling(\n",
    "    theta, executable: Program, marked_qubits: List[int], backend, samples: int\n",
    ") -> float:\n",
    "    \"\"\"Calculate the expectation value of the Zi operator where i ranges over all \n",
    "    qubits given in marked_qubits.\n",
    "    \"\"\"\n",
    "    bitstring_samples = backend.run(\n",
    "        executable.write_memory(region_name=\"theta\", value=theta)\n",
    "    ).readout_data.get(\"ro\")\n",
    "    bitstring_tuples = list(map(tuple, bitstring_samples))\n",
    "\n",
    "    freq = Counter(bitstring_tuples)\n",
    "\n",
    "    \n",
    "    exp_val = 0\n",
    "    for bitstring, count in freq.items():\n",
    "        bitstring_int = int(\"\".join([str(x) for x in bitstring[::-1]]), 2)\n",
    "        if parity_even_p(bitstring_int, marked_qubits):\n",
    "            exp_val += float(count) / samples\n",
    "        else:\n",
    "            exp_val -= float(count) / samples\n",
    "    return exp_val"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44cc3833",
   "metadata": {},
   "source": [
    "Calculate the parity of elements at indexes in marked_qubits. The function is a \n",
    "modified version of `parity_even_p` from the `VQE` function in Grove [[1]](#references)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "4aa4dba8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def parity_even_p(state, marked_qubits):\n",
    "    mask = 0\n",
    "    for q in marked_qubits:\n",
    "        mask |= 1 << q\n",
    "    return bin(mask & state).count(\"1\") % 2 == 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09c420e0",
   "metadata": {},
   "source": [
    "## Run VQE first without error mitigation and then with ZNE, and compare results\n",
    "\n",
    "Scan over the parameter `theta` and calculate energy expectation, without mitigation. \n",
    "In a later section we will plot these results and compare them with the results from \n",
    "ZNE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "aa174530",
   "metadata": {},
   "outputs": [],
   "source": [
    "thetas = np.linspace(0, 2 * np.pi, 51)\n",
    "results = []\n",
    "\n",
    "for theta in thetas:\n",
    "    results.append(executor(theta, backend, readout_qubit, samples, hamiltonian, program))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c729c9e",
   "metadata": {},
   "source": [
    "Optimization routine without mitigation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "94f260a3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " final_simplex: (array([[2.9953125 ],\n",
      "       [2.99560547]]), array([-0.432     , -0.42733333]))\n",
      "           fun: -0.432\n",
      "       message: 'Optimization terminated successfully.'\n",
      "          nfev: 31\n",
      "           nit: 11\n",
      "        status: 0\n",
      "       success: True\n",
      "             x: array([2.9953125])\n"
     ]
    }
   ],
   "source": [
    "init_angle = [3.0]\n",
    "\n",
    "res = optimize.minimize(\n",
    "    executor,\n",
    "    init_angle,\n",
    "    args=(backend, readout_qubit, samples, hamiltonian, program),\n",
    "    method=\"Nelder-Mead\",\n",
    "    options={\"xatol\": 1.0e-3, \"fatol\": 1.0e-2},\n",
    ")\n",
    "print(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "606cbbe3",
   "metadata": {},
   "source": [
    "The result on the unmitigated noisy circuit result in loss of accuracy (relative to \n",
    "the ideal expectation value of -1.0) and additional iterations required to reach \n",
    "convergence. \n",
    "\n",
    "\n",
    "### Now we introduce ZNE and compare results.\n",
    "\n",
    "This is done by wrapping the noisy executor into a mitigated executor. \n",
    "We will fold the gates from the right and apply a linear inference (using a Linear \n",
    "Factory object) to implement ZNE. You can read more about noise scaling by unitary folding in the Mitiq user guide."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "cd551693",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mitigated_expectation(\n",
    "    thetas, backend, readout_qubit, samples, pauli_sum, executable: Program, factory\n",
    ") -> float:\n",
    "    \"\"\"\n",
    "    This function is the ZNE-wrapped executor, which outputs the error-mitigated \n",
    "    expectation value. \n",
    "\n",
    "    Args:\n",
    "      thetas: the input parameter for the optimization\n",
    "      backend: the quantum computer that runs the quantum program\n",
    "      readout_qubit: declared memory for the readout\n",
    "      samples: number of times the experiment (or simulation) will be run\n",
    "      pauli_sum: the Hamiltonian expressed as \n",
    "      executable: the pyQuil quantum program\n",
    "      factory: factory object containing the type of inference and scaling parameters\n",
    "\n",
    "    Returns:\n",
    "      The error-mitigated expectation value as a float.\n",
    "    \"\"\"\n",
    "    mitigated_exp = zne.execute_with_zne(\n",
    "        executable,\n",
    "        lambda p: executor(thetas, backend, readout_qubit, samples, pauli_sum, p),\n",
    "        factory=factory,\n",
    "        scale_noise=fold_gates_from_right,\n",
    "    )\n",
    "    return mitigated_exp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0a4f3bc",
   "metadata": {},
   "source": [
    "Here we use a linear inference for the extrapolation. \n",
    "See the section on [Factory Objects](../guide/zne-3-options.md#extrapolation-methods-factory-objects) \n",
    "in the Mitiq user guide for more information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f6fe2457",
   "metadata": {},
   "outputs": [],
   "source": [
    "fac = mitiq.zne.inference.LinearFactory(scale_factors=[1.0, 3.0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8d2628f",
   "metadata": {},
   "source": [
    "Scan over the parameter `theta` and plot the energy expectation with error mitigation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "1b9658c5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEeCAYAAAC+OaPqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABl6ElEQVR4nO2dd3hUZdbAfyeTHiChE5rAiihIqKIiKoqCCouCorCrgq6Luiro5+Laja6uuLoiqGvvFUUEFJUisGJBujRBiigldBJaElLO98e9EyYz9yYzYVImeX/PM09m3vvee8+U3HPPe5qoKgaDwWAwlDdRlS2AwWAwGGoGRuEYDAaDoUIwCsdgMBgMFYJROAaDwWCoEIzCMRgMBkOFYBSOwWAwGCoEo3AMVQoReVFEHihh+70i8mpFyuSGiLwpIo9WthyRiIiki8i7lS2HoWIxCsdQIYjIZhE5KiIN/MaXiYiKSCsAVb1JVf9pb+stIlt956vqv1T1hjDJpCJyYjiO5XDsESJSICKH/B5Ny+N84cIoAkN5YhSOoSL5FRjmfSEiHYHEyhOn3PlBVWv5PbaH8wQiEh3O4xkM5YlROIaK5B3gWp/Xw4G3fSd4l6lEJAn4Emjqax3434GLyLUi8puI7BWRB2xL6gJ7Ww8R+UFEMkUkQ0SeE5FYe9s39iF+so99lT0+QESW2/t8LyJpPufqIiJLReSgiEwE4sv6Qdhy/l1EVohIlohMFJF4n+0lybFZRP4hIiuAwyIS7fY5iEgTETkiIvV99u8qIrtFJCZEme8WkY32+18jIoN8to0QkW9F5CkR2S8iv4rIxT7bW4vI/+x9ZwENfLbFi8i7tuyZIrJIRBrb2+qJyBsist0+7hR7vK6IfG6/j/328+Y+x5wnIo+LyEIROSAiU0Wkns/2M+zPNVNEfhKR3qF8FoayYRSOoSJZANQRkVNExAMMBRyXb1T1MHAxsN3NOhCR9sB/gT8DqUAy0MxnSgFwB9bF7UygD/A3+/jn2HM62ceeKCJdgNeBG4H6wEvANBGJsxXVFCylWQ/4GLj8eD4M4ErgIqA1kAaMsN+Xqxw++w4D+gMpwElun4Oq7gDm2efycg3woarmhSjvRuBs+/gPA++KSKrP9tOBdVif97+B10RE7G3vA0vsbf/EutnwMtw+Zgv7/d4EZNvb3sGygjsAjYBx9ngU8AZwAtDSnv+cn7zXAtdjfSb5wAQAEWkGTAcexfou/w58IiINQ/w8DCFiFI6hovFaORcCPwPbjuNYVwCfqeq3qnoUeBAoKg6oqktUdYGq5qvqZqwL97klHG8k8JKq/qiqBar6FpALnGE/YoBnVDVPVScBi0qR7wz7Dtr72Oi3fYKqblfVfcBnQOcg5PDdd4uqZpf2OQBvAVcD2Ip+GNb3EBKq+rEtb6GqTgTWAz18pvymqq+oaoF9zlSgsYi0BE4DHlDVXFX9xn6/XvKwFM2J9vtdoqoHbGV2MXCTqu63P/f/2bLsVdVPVPWIqh4EHiPwu31HVVfZNy8PAFfa7/9q4AtV/cJ+L7OAxcAloX4mhtAw67+GiuYd4Busu/q3S5lbGk2BLd4XqnpERPZ6X4vIScDTQHesu+RorLtsN04AhovIbT5jsfZ5FNimxavd/laKfAtUtVcJ23f4PD9in6c0Obxs8Xle4ucATAVeFJHWQDsgS1UXliJ7ACJyLfB/QCt7qBY+S2O+78eWwXfOfvvC7+U3LIsGrN9EC+BDEUnBsnrvs8f2qep+B1kSsaydi4C69nBtEfHYCg+Kf0a/Yd0wNMD6fIeIyB99tscAc0v+BAzHi7FwDBWKqv6GFTxwCTC5tOmlbM8AfNftE7DulL28AKwF2qpqHeBeQHBnC/CYqqb4PBJV9QP7XM18lojAWsopD0qSw4vvZ1Pi56CqOcBHWHf211AG60ZETgBeAW4F6qtqCrCKkj9PX/nqiuWX81L02dmWy8Oq2h7oCQzAsoK3APVsJeTPnVjK83T7u/UukfrK08LneUssS2qPfdx3/D7fJFUdG8R7MRwHRuEYKoO/AOf73fE6sROoLyLJLtsnAX8UkZ62jyWd4hec2sAB4JCInAzc7HD8Nj6vXwFuEpHTxSJJRPqLSG3gByw/wCgRiRGRwRRfTgonJcnhRGmfA1jW5AhgIKUrnCjbke99xAFJWEpuN4CIXAecGsybsW8yFgMPi0isiPQCiqwLETlPRDray10HsBRDoapmYAWO/NcOEogREa9iqY3lt8m0gwEecjj11SLS3raGHgEm2dbPu1ifVz8R8djvsbdv0IGhfDAKx1DhqOpGVV0cxLy1wAfAJtsH0tRv+2rgNuBDrLvoQ8AuLH8HWM7gPwEHsS7iE/1OkQ68ZR/7Slumv2I5n/cDG7Ad+bZvZLD9eh9wFaVbaGdKYB7OaUG8b1c5XOaX9jmgqt8BhcBSWwGUxDCsi7n3sVFV1wD/wVK8O4GOwHelvRcf/oQVVLAPSzn4Lqc2wVKaB7D8ev/jmFK8BksBrbXf0+32+DNAApbFsgD4yuGc7wBvYi31xQOjAFR1C3AplsW7G8viGYO5HpY7YhqwGaoLIlILyMRaQvu1ksWpNNw+BxGZA7yvqlWiUkN5IiLzgHdrwnuNJIxGN0Q0IvJHEUm0/QNPASuBzZUrVcVT2udgW1ZdCbTyDIYKwygcQ6RzKbDdfrQFhmrNNNtdPwcReQuYDdxuhxAbDJWCWVIzGAwGQ4VgLByDwWAwVAgm8bMEGjRooK1atapsMQwGgyGiWLJkyR5VDSgVZBROCbRq1YrFi0uN3jUYDAaDDyLiGHpvltQMBoPBUCEYhWMwGAyGCsEoHIPBYDBUCMaHY4gI8vLy2Lp1Kzk5OZUtiqEKEh8fT/PmzYmJCamnnKGCMQrHEBFs3bqV2rVr06pVK4oXbDbUdFSVvXv3snXrVlq3bl3Z4hhKIKKW1ETkdRHZJSKrXLaLiEwQkQ1ite7t6rNtuIistx/DnfYPB1OWbeOssXNoffd0zho7hynLjqe/mMFLTk4O9evXN8rGEICIUL9+fWP9hoHyvn5FmoXzJlYFXbfGXRdjlfVoi1WZ9gXgdJ/y5d2xSqwvEZFpTo2djocpy7Zxz+SVZOdZ/Z+2ZWZzz+SVRdufnLGO7ZnZNE1JYEy/dlzWpZnboQwOGGVjcMP8NkJjyrJtAdcjwPX6Fa5rVUQpHFX9RkRalTDlUuBtu4bUAhFJsdvU9gZm2a18EZFZWJ0CP3A9Uhl4csa6oi/LS3ZeAQ9NW8XR/EKy8wqB8vkiDQaDIRicbozvnryCWE+U4/XryRnraqbCCYJmFG8ru9UecxsPQERGYvWUp2XL0Bo6bs/MdhzPys4PGPN+kWAsn0ihVq1aHDp0KGB8xIgRDBgwgCuuuCLkY6anp1OrVi3+/ve/BzVeFo5HPkNk42TJON0Y5+QVkmPfEPvjdl0rC9VN4Rw3qvoy8DJA9+7dQ6ps2jQlgW0hfDnbMrP5xycryM03lk+4cfpHM5+poSbhZMncNeknjhaEVrC5aUpC2GSKqKCBINhG8T7mze0xt/GwMqZfOxJiPMXGEmI81E10D9X0KhsvvpaPoWx4/9G2ZWajHFPk4XKAqiq33nor7dq144ILLmDXrl1F25YsWcK5555Lt27d6NevHxkZGQC88sornHbaaXTq1InLL7+cI0eOBH0+t31HjBjBqFGj6NmzJ23atGHSpEmlynf33XfTvn170tLSiqynnTt3MmjQIDp16kSnTp34/vvvAbjsssvo1q0bHTp04OWXXy46Rq1atbjjjjvo0KEDffr0Yffu3QBs3LiRiy66iG7dunH22Wezdu3asny8hjDhZMmUpGxSEmIcr19e/044qG4WzjTgVhH5ECtoIEtVM0RkBvAvEalrz+sL3BPuk3vvoEtzxoH1Rfr/GLyE04Stjjz82WrWbD/gun3Z75kcLQhU5HdNWsEHC3933Kd90zo89McOQZ3/008/Zd26daxZs4adO3fSvn17rr/+evLy8rjtttuYOnUqDRs2ZOLEidx33328/vrrDB48mL/+9a8A3H///bz22mvcdtttQZ2vpH0zMjL49ttvWbt2LQMHDuSKK65wlW/v3r18+umnrF27FhEhMzMTgFGjRnHuuefy6aefUlBQULRs+Prrr1OvXj2ys7M57bTTuPzyy6lfvz6HDx+me/fujBs3jkceeYSHH36Y5557jpEjR/Liiy/Stm1bfvzxR/72t78xZ86coN6jIfyUdB3xv/4kxHhIH2j9/stzZSCiFI6IfIAVANBARLZiRZ7FAKjqi8AXwCVYPeCPANfZ2/aJyD+BRfahHvEGEISby7o0c/2CnNZSnZbgasdH88mSLTw9a71ZEioD/sqmtPFQ+eabbxg2bBgej4emTZty/vnnA7Bu3TpWrVrFhRdeCEBBQQGpqakArFq1ivvvv5/MzEwOHTpEv379gj5fSftedtllREVF0b59e3bu3FmifMnJycTHx/OXv/yFAQMGMGDAAADmzJnD229bgZ8ej4fk5GQAJkyYwKeffgrAli1bWL9+PfXr1ycqKoqrrroKgKuvvprBgwdz6NAhvv/+e4YMGVIkW25uboifrKGs+C8hX3Vac6I9Qp6DRdPM5/rjdH0pz+tMRCkcVR1WynYFbnHZ9jrwennIFQxuisjf8okSOJCTz98/XoH3p2J8O8UpzRI5a+wcR0XeLCWBiTeeWV5ioap06NCBH374IWDbiBEjmDJlCp06deLNN99k3rx5QR+3pH3j4uKKnb8koqOjWbhwIV9//TWTJk3iueeec7VA5s2bx+zZs/nhhx9ITEykd+/ernkuIkJhYSEpKSksX7486PdlCA9OvpqnZ60nziPEeqTYMpp3iaykG+PypLr5cCKKy7o04/HBHWmWkoBgXRCfvrIz9ZJi8b90GN9O8Lj50sK1Fn3OOecwceJECgoKyMjIYO7cuQC0a9eO3bt3FymcvLw8Vq9eDcDBgwdJTU0lLy+P9957L6Tzhbqvm3yHDh0iKyuLSy65hHHjxvHTTz8B0KdPH1544QXAssqysrLIysqibt26JCYmsnbtWhYsWFB0/MLCwiJ/0fvvv0+vXr2oU6cOrVu35uOPPwYs5ec9vqF8cfLVANSrFce/r+hU7Pry+OCOlXrTGlEWTnXE6U7jjonLHeca305wuPnSwvWPNmjQIObMmUP79u1p2bIlZ55pWU2xsbFMmjSJUaNGkZWVRX5+PrfffjsdOnTgn//8J6effjoNGzbk9NNP5+DBg0GfL9R93eQ7ePAgl156KTk5OagqTz/9NADjx49n5MiRvPbaa3g8Hl544QUuuugiXnzxRU455RTatWvHGWecUXT8pKQkFi5cyKOPPkqjRo2YOHEiAO+99x4333wzjz76KHl5eQwdOpROnTqF9NkaQsfturAjK6fSLBk3pDQzvCbTvXt3rYwGbG5LQo3rxPHjvRdUuDxVgZ9//plTTjmlssUw4J6PVNnUxN+IqpKWPpODuYG5fs1SEvju7vMrQSoQkSWq2t1/3Fg4VZAx/doF+HYADufmM/7rX/ho0VYTTGAw1EB8gwNSk+NpkhzPwdx8PFFCQWGgr6aqYXw4VRAn3849F59MjCeKcbPWl1t+icEQDFXRuqkJ+OeXbc/KYenvmfTv2ISnrkirUr4aN4yFU0VxWnt947vNQF6xsXDXOjIYDFUTt+CA5VuyeP7P3RjUtXklSBUaxsKJIHYecA5LNcEEBkP1x+3/PJL+/43CiSDcahqFs9aRwWComjRJjnccj6T/f6NwIgin/BKA/h2bVII0BoOhojhyNJ8YT2DPn6oaHOCGUTgRhH8wQWpyPM1T4nl7wW8s+a1cKvUYfBAR7rzzzqLXTz31FOnp6SXu8+KLLxaVjSkrb775JrfeeutxHcNLeno6Tz31VFiOZagYjuYXctO7S9m6P5vrz2oVEcEBbpiggQjDP5hgz6Fchrz4A39+ZQF1EmLZfTDXhEsDrPgIvn4EsrZCcnPo8yCkXXlch4yLi2Py5Mncc889NGjQIKh9brrppuM6p6Fm4hv+HG8X2nzi8o5cdVpLHgyyyGxVxFg4EU6DWnFcfUZLcvKVXQdzTbg0WMrms1GQtQVQ6+9no6zx4yA6OpqRI0cybty4gG2bN2/m/PPPJy0tjT59+vD771ZVal+LYsKECUWtAYYOHUphYSFt27YtKu9fWFjIiSeeWPTaic8++4zTTz+dLl26cMEFFxQV7ExPT+f666+nd+/etGnThgkTJhTt89hjj3HSSSfRq1cv1q07Vh7JXx6wQp6vu+46OnbsSFpaGp988gkAN998M927d6dDhw489NBDRcdo1aoVd911Fx07dqRHjx5s2LABgN27d3P55Zdz2mmncdppp/Hdd9+F/oHXUPzDn7PzCoiOEuKiA5fTIw1j4VQDXv92c8BYtQ6X/vJu2LHSffvWRVDgV6k4Lxum3gpL3nLep0lHuHhsqae+5ZZbSEtL46677io2fttttzF8+HCGDx/O66+/zqhRo5gyZUqxOWPHjuXXX38lLi6OzMxMoqKiuPrqq3nvvfe4/fbbmT17Np06daJhw4au5+/VqxcLFixARHj11Vf597//zX/+8x8A1q5dy9y5czl48CDt2rXj5ptvZsWKFXz44YcsX76c/Px8unbtSrdu3RzlAauMTnJyMitXWp/v/v37AUtp1atXj4KCAvr06cOKFStIS0sDKJr/9ttvc/vtt/P5558zevRo7rjjDnr16sXvv/9Ov379+Pnnn0v9fA3O4c/5hVot/p+NhVMNqA7hkmHFX9mUNh4CderU4dprry1mQQD88MMP/OlPfwLgmmuu4dtvvw3YNy0tjT//+c+8++67REdb93rXX399kY/n9ddf57rrrivx/Fu3bqVfv3507NiRJ598sqg4KED//v2Ji4ujQYMGNGrUiJ07dzJ//nwGDRpEYmIiderUYeDAgSXKM3v2bG655VjB9bp1rRZSH330EV27dqVLly6sXr2aNWvWFM0ZNmxY0V9v4dLZs2dz66230rlzZwYOHMiBAwdMwmiQVOf/Z2PhVAPcWltHUrhkSJRmiYw71V5O8yO5BVw3/bhPf/vtt9O1a9dSlYM/06dP55tvvuGzzz7jscceY+XKlbRo0YLGjRszZ84cFi5cWGo16Ntuu43/+7//Y+DAgcybN69Y0IJvqwKPx0N+fmB9rdLkceLXX3/lqaeeYtGiRdStW5cRI0YUa1UgIgHPCwsLWbBgAfHxzqG8Bnca1Ylj54HAm6Pq8P9sLJxqgFu49B87pVaCNFWAPg9CjN8/Z0yCNR4G6tWrx5VXXslrr71WNNazZ08+/PBDwKqafPbZZxfbp7CwkC1btnDeeefxxBNPkJWVVXTHf8MNN3D11VczZMgQPJ6S1+mzsrJo1sxaVnnrLZflQR/OOeccpkyZQnZ2NgcPHuSzzz4rUZ4LL7yQ559/vmj//fv3c+DAAZKSkkhOTmbnzp18+eWXxc7hrRY9ceLEosrUffv25dlnny2aY/rkBEdufgGxnsDLcqSFP7sRUQpHRC4SkXUiskFE7nbYPk5EltuPX0Qk02dbgc+2aeUl4/RN0+k7qS9pb6XRd1Jfpm86/jvq0nAKl25SJ46PFm+tFmZ4yKRdCX+cYFk0iPX3jxOOO0rNlzvvvJM9e/YUvX722Wd54403SEtL45133mH8+PHF5hcUFHD11VfTsWNHunTpwqhRo0hJSQFg4MCBRc760khPT2fIkCF069YtqEi5rl27ctVVV9GpUycuvvhiTjvttBLluf/++9m/fz+nnnoqnTp1Yu7cuXTq1IkuXbpw8skn86c//Ymzzjqr2Dn2799PWloa48ePLwqomDBhAosXLyYtLY327dvz4osvliqrAR79/Ge27M/murNOqJTw5/K+fkVMewIR8QC/ABcCW7HaRQ9T1TUu828Duqjq9fbrQ6paK5RzhtqeYPqm6aR/n05OwbHlhnhPPOk90+nfpn8opz5uNu4+xKXPfUfbxrWYOPJMYqMj6t4igOpcen7x4sXccccdzJ8/v7JFCZlWrVqxePHioMPEy5NI/41MXb6N0R8u569nt+a+/u0r/PzhvH5Vh/YEPYANqroJQEQ+BC4FHBUOMAx4yGVbuTB+6fhiXxZATkEO45eOL9q+4/AOmiQ1YXTX0fRv05/pm6Y7jh8vf2hYiycuT+OW95fyl7cWsWn3YdPSoAoyduxYXnjhhZC7gBqqB775NgCtGyRy10UnV4osJV2/wnXDHEkWzhXARap6g/36GuB0VQ1IwRaRE4AFQHNVLbDH8oHlQD4wVlWnuJxnJDASoGXLlt1+++23oGVMeysNDWgObRETFUNe4bFKz/GeeC498VKmbpjqeEcBzgoqVK597Ue+Wb+n2FhCjCfiMpQj/e7VUP5E2m/Em2/jGwIdHx3F2MvTKuV/0+36JQgrhq8I6VjVwcIJhaHAJK+ysTlBVbeJSBtgjoisVNWN/juq6svAy2AtqYVy0iZJTcg4nOG4zVfZgHXnMHHdxIB5OQU5/POHf1KgBUWKKONwBunfpwOErHQ27AoMRa3WOToGQ4TglG+Tk19Ypv/NUFdKfOc3TmzMqfVPdZ3bJCl8tRojaWF/G9DC53Vze8yJocAHvgOqus3+uwmYB3QJt4Cju44m3lM8DNT/dTAczj9c4tJcKGRkmZYGBkNVJFz5Nl7fS8bhDBQtukF1c/j7z99xZAezt8ymWVIz4jxxxebGe+IZ3XV0SPKURCQpnEVAWxFpLSKxWEolINpMRE4G6gI/+IzVFZE4+3kD4CzcfT9lpn+b/qT3TCc1KRVBSE1KLXrtRJSE9vHvOLwjZJlMSwODoWrSuE542g2U5Htxijr7z+L/BMwHKKCAh3s+HHD9CmfAU8QsqalqvojcCswAPMDrqrpaRB4BFquqV/kMBT7U4s6pU4CXRKQQS8mOdYtuO176t+nv+AU5RX+4+XDio+PJzM0MOEZSTBJTN0zl+eXPB206j+nXLmCdOC46qlrE9BsMkYqq0qBWLDv8miqWJd/G7UY043AG9393P/mF+UWv751/L4UUuh7H7foVLiLJwkFVv1DVk1T1D6r6mD32oI+yQVXTVfVuv/2+V9WOqtrJ/vua/7HLEzfL5/4z7nccv7vH3QFLcR7xcCjvEA9890DQpjME5uhECTSqHceAtBqaFHociAhXX3110ev8/HwaNmzIgAEDAJg2bRpjx1pVEKZMmVKs/MuDDz7I7Nmzy3Te5cuX88UXX4S8X+/evXEK63cbLwutWrUqlo9kCI6py7ezavsBBnZKDSnfxt9iuXf+va6BSkCRsvFSSCFCYF8dCK+vxo2IsXAiHbc7h5LuKPydgGMXjg2wfIIJW/RtaTB9RQa3vL+Ul77ZxC3nnVj2N1TFKY9w86SkJFatWkV2djYJCQnMmjWrKOsfrAROb62yKVOmMGDAANq3t/IpHnnkkTKfd/ny5SxevJhLLrnkuOQ3VA12Hsjhwamr6HZCXcZd1QVPlLMC8Mc/TybjcAafbfqM1MRU9uXuI9enVmC8J95x2QxA0YDt4fbVuBFRFk5Non+b/sy8YiYrhq9g5hUz6d+mP1m5WY5zQ/Ht9E9LpX9aKs/M/oW1Ow6ES9wqRahO1FC45JJLmD7dOs4HH3xQVLgSjjVK+/7775k2bRpjxoyhc+fObNy4kREjRjBp0iQAvvjiC04++WS6devGqFGjiiykhQsXcuaZZ9KlSxd69uzJunXrOHr0KA8++CATJ06kc+fOTJw4kcOHD3P99dfTo0cPunTpwtSpUwHIzs5m6NChnHLKKQwaNIjs7NKdzyW1HXjooYfo2rUrHTt2ZO3atQDs3buXvn370qFDB2644Qa8K9eHDx+mf//+dOrUiVNPPbWo3M2iRYvo2bMnnTp1okePHhw8eJDNmzdz9tln07VrV7p27cr3338PwLx58zjnnHPo378/7dq146abbqKw0Fr+mTlzJmeeeSZdu3ZlyJAhEVsIVFW5Z/JKjhYU8uQVaUErG3D21QAgOPpe3HzHvtvLy1fjhrFwIgi3sOtQTeF/XnoqP27ay1/eXIQCGZk5EZUQ+sTCJ1i7b63r9hW7V3C08GixsZyCHB787kEm/TLJcZ+T653MP3r8o9RzDx06lEceeYQBAwawYsUKrr/++oAKAT179mTgwIEMGDCAK664orgcOTnceOONfPPNN7Ru3bqYwjr55JOZP38+0dHRzJ49m3vvvZdPPvmERx55hMWLF/Pcc88BcO+993L++efz+uuvk5mZSY8ePbjgggt46aWXSExM5Oeff2bFihV07dq11PdTUtuBBg0asHTpUv773//y1FNP8eqrr/Lwww/Tq1cvHnzwQaZPn15UT+6rr76iadOmRco4KyuLo0ePctVVVzFx4kROO+00Dhw4QEJCAo0aNWLWrFnEx8ezfv16hg0bVrTEt3DhQtasWcMJJ5zARRddxOTJk+nduzePPvoos2fPJikpiSeeeIKnn36aBx8MT228isCb4Oktsjuoc1PaNAyp8InrjWVJvhcn37HX2q/o6idgLJyIwinsGqBPyz4hHadeUix/7NSUbZk5bM/MqXZN2/yVTWnjoZCWlsbmzZv54IMPyrTEtXbtWtq0aUPr1q0BiimcrKwshgwZwqmnnsodd9xRrPWALzNnzmTs2LF07tyZ3r17k5OTw++//84333xT5GNKS0srUhwlUVLbgcGDBwPQrVs3Nm/eDFDsHP379y9qX9CxY0dmzZrFP/7xD+bPn09ycjLr1q0jNTW1qH5bnTp1iI6OJi8vj7/+9a907NiRIUOGFDtnjx49aNOmDR6Ph2HDhvHtt9+yYMEC1qxZw1lnnUXnzp156623CCUhu7Lxbajm5avVO0L+X2uU2Mhx3O2G0813XBmKxouxcCII7w/FN2ErxhPDxHUTOaf5OZzZ9MygjzVz9c6AsUhJCC3NEuk7qa+jJZialMobF71x3OcfOHAgf//735k3bx579+497uN5eeCBBzjvvPP49NNP2bx5M71793acp6p88skntGt3fJGGpbUd8LY7CKbVwUknncTSpUv54osvuP/+++nTpw+DBg1ynDtu3DgaN27MTz/9RGFhYbEWBr6tDryvVZULL7yQDz74wP9QEYFTgmd2XmgJnkcLjjrebJbme6ksS8YNY+FEGL6+nVlDZvFB/w9oldyKW2bfQu+JvYOu8lqdmzy5JeCGyyl6/fXX89BDD9GxY0fXObVr1+bgwYMB4+3atWPTpk1FFoPX1wHFWw+8+eabrsfq168fzz77bJH/ZNmyZYDViuD9998HYNWqVaxYUXI5ktLaDjjhe44vv/yyqCPo9u3bSUxM5Oqrr2bMmDEsXbqUdu3akZGRwaJFiwA4ePAg+fn5ZGVlkZqaSlRUFO+88w4FBccuxgsXLuTXX3+lsLCQiRMn0qtXL8444wy+++67ovbVhw8f5pdffilV1qrC8f6vqSrp36fz28HfGNZuWJWyWELFWDgRTnJcMleedCX/+vFf7M2x7raDKYVTnZu2+VuC4SyKCtC8eXNGjRpV4pyhQ4fy17/+lQkTJhQFCwAkJCTw3//+l4suuoikpKSi5SaAu+66i+HDh/Poo4/Sv/8xWc8777yiJbR77rmHBx54gNtvv520tDQKCwtp3bo1n3/+OTfffDPXXXcdp5xyCqecckpRK2k3fNsOtGjRIqDtgBMPPfQQw4YNo0OHDvTs2ZOWLVsCsHLlSsaMGUNUVBQxMTG88MILxMbGMnHiRG677baiyL7Zs2fzt7/9jcsvv5y333676HPwctppp3HrrbeyYcMGzjvvPAYNGkRUVBRvvvkmw4YNIzfXisR69NFHOemkk0qVtyrQJDneseJHSf9rvlGWtWJqcTDvILd2vpUbO93IvdxbnuKWKxFTvLMyCLU9QWVR0hLSzCtmOu7jWDgwJoqxgyuncGBpRFphxpI4dOgQtWrVQlW55ZZbaNu2LXfccUdli1XpzJs3j6eeeorPP/+8TPtX1d/IiNd/ZN4vwRfQdWoT4BEPj571KAP+MKDc5Q0HbsU7zZJaNaCk6BU3/BNCATo3T6mSyqa68corr9C5c2c6dOhAVlYWN954Y2WLZCgn1u44wPwNezmzTb2gEzydwp8LtIAJyyZUgMTli1lSqwa4hUvXT6hf4n6+CaGPTV/DK/N/Zenv++nasm65yGmwuOOOO4xF40Dv3r1dAyUiEVXlgSmrqBMfzX//3I26SbEBc/wTlG9Mu9G14nxZailWNYyFUw1wcpILwsHcg/y0+6fgjnHBSTSuE8cDU1ZRUFg1l1nN8q/Bjar42/hk6TYWbd7P3Ref7KpsAhKUf0h3PV5FlJ4pb4zCqQY4xdvfddpdNEpqxMiZI3l22bOl9imvFRfNAwPas3r7Ad5dUPVyHOLj49m7d2+VvLAYKhdVZe/evcXCqyubrCN5PP7Fz3RtmcKQbi0c57hVDkiKTirXKMvKxAQNlECkBA24sevILq76/Cr2ZBd3WLr1KVdVrn19IQs37aVuUhw7D1SdCgR5eXls3bq1WJ6IweAlPj6e5s2bExMTU6ly+FcU+Hu/k7j1vLaOc0vqsPn42Y+XW5RlRRD2jp8iEgWcg9UqYJ5fd01DFaBRYiM84gkYdyv4KSKc3bYB89fvKSqb7q1AAFSq0omJiSnKzjcYqiJOkZ/Pz9lI85REx/+dRomN2HkkMAG7SVKTKpewGS5CWlITkUQRGSwibwO7gDnALGC3iLwrIleISGgFggzlyq4juxzHdxze4dic6a3vA5fTvBUIDAaDO84VBdz/d1rWbhkwVl2WztwoVeGISGMRuUFEPgf2AJOAM4A3sCycs4BXge7AR1jK5wsRGSkike/linDcHI1REsVD3z8UUFF5V+H3jvOrQwUCg6E8CaWiwKo9q1i8czG9mvaK6MoBoRLMktp2++9i4J/AVIdumT8Ad9ntnS8DBgIvAP8VkQRVzQuHsCJyETAeaxnvVVUd67d9BPAk4K2K95yqvmpvGw7cb48/qqpvhUOmqs7orqMDksiio6LJL8wvVlIErKW2hMYzyTvQJeA41aECgcFQnqQmx7M9iIoChVrIYwseo35CfZ4890lqxdacRaFgltT+BjRX1dNV9fGSWjOr6lpVHauqPYFmwM1AWHw7IuIBngcuBtoDw0SkvcPUiara2X54lU094CHgdKAH8JCI1IhkE6cItkfPetS1659GZ5IQU9zvEx9jWlIbDKXRqUVKwJhTy+jJ6yezau8q7ux+Z41SNhCEhaOqL5XlwKq6A3ilLPu60APYoKqbAETkQ+BSwFUB+tAPmKWq++x9ZwEXAZFZfjZEnByQ45eOdymH04S/De7IkzPWsT0zGwV6n9So0qPUDIaqzK4DOfzvl92kNavD3sN5bM/MdozwzMzJ5Jmlz9CtcTf6t66+S2duRFKlgWbAFp/XW7EsFn8uF5FzgF+AO1R1i8u+jldQERkJjASKChNWR5yW2o41ZzpWgeDmd5cwb91udh7IoXGdqpPnYDBUJZ6e9Qt5BYVMGNaVVg2SArZ7Kwp4b/J6NesV0IqhJlCmxE8ROUlEBonIjXZwwCARcQ42r1g+A1qpahpW9FzIfhpVfVlVu6tq94YNG4ZdwKqC71Kbl9u63BZgCd198cnkFxbyn5kmSs1gcOLnjANMXLyFa85o5apsvBUFvLz000thaXkeaQStcETkFBEZLyLbgJ+xotVeAF60n68Vke0i8oyIlEfJ1m2Ab8puc44FBwCgqntVNdd++SrQLdh9ayLe3jpzhswh3hPv2Lb5hPpJDD+zFR8v2cqa7QcqQUqDoeqiqvzri5+pEx/DqD4nOs5xqijgzYWraQQTFv0HEZkErAL+AqwAHgauBS4B+tvPHwF+Am4AVonIxyLSJoyyLgLaikhrEYkFhgLT/GRN9Xk5EEsxAswA+opIXTtYoK89ZgAaJjZk2CnD+HzT52zYvyFg+23ntyU5IYZ/ffGzKS1jMPgw75fdzF+/h1F92pKSGFgvDajWxThDJRgfzhpgJTACmKyqh0uaLCJJwBXAaHvfsCz8q2q+iNyKpSg8wOuqulpEHgEWq+o0YJSIDATygX22zKjqPhH5J5bSAnjEG0BgsLi+w/V8vO5jnl/+POPOG1dsW3JiDKPOb8sjn6+h+6Oz2Xf4aJUpeWMwVAZTlm3j3zPWsj0zB0+UkBxvXUr9qz+fXPdk12NUh2KcoVJqLTURGWhfzEM/uMilqjq1TJJVAcpUS23FR/D1I5C1FZKbQ58HIe3K8hEwzLzw0wv8d/l/+aD/B5za4NRi2yYt3sKYSSuKVX4qqYmUwVBdcSphkxDjYeh5u/l8+4SA5bPWdVqTcTgjIECnUpM8y/k6VeYGbGVVNva+EatsysSKj+CzUZC1BVDr72ejrPEVH8G4UyE9xfq74qPKljaAa9tfS924ujy77NmAbeNmrw8oM2hK3hhqIm4lbD759RXH6s85BTkBuXCVrmzcrlPlTNBh0XaNtFOBZOAgsMnOtTF4+foRyPMrY5GXDV+MgYLcY9u8X7DvflXAIkqKSeIvHf/CU4ufovfE3uzL2VdUqXZ7pvM+puSNoabh9psv9Ox3TKfecXhH5RXjdLJk3K5TXz9S7teeUhWOiMQBz2D5Q2L9tu0AvgYmAl9oTfcoZ211Hs/JDBzLy4bP/w8K8yDfvivyVUSVpHTqxNYBYG/OXoCiGmsNmgxh944OAfNNyRtDTaNJcjwZDiVsogpS0OjMwPmV5avxWjK+N7pTbobCfOf5WVvLfaktmLDop4AbgXnAvcDf7THBCgi4GitabJWI9AmbZJFIcvPQ5h89eEzZePHeaVQSL/z0QsBYTkEOcY1mmJI3BgPQrnFgOZqEGA8n1AlstFZh1Z+dluydLBk3ZQOAwqc3lutSWzAK5yqsiLCLVfUJVR0H/NveNgRoA9wNJAAzRKTmNmvv8yDE+N3xxyRAQr3QjuNmKVUAbqGaB/J28/jgjjRLSShaNuh1YgMTMGCoUWzPzOb7Tfvo0apu0f9Cs5QEru1ziM1HVtK7ee+K99U4+WSm3mK/dsH/OhWdADGJoIXFx8N8AxyMDycBqxq0I6q6GXhSRJ4B/gU8JSKrVHVWWCSMJLymp79JCsVNW7C+8OgEyHaIzvbEwsJX4LvxFe7baZLUxDFvoElSEy7rcqzkzagPljFrzU52HcyhUW1T8sZQM5jw9XpQePqqzjSvmwjA3uy9DJ42mJPrnczTvZ8mxlOOXUeD9ckUHHU/RnKLY/v5HmfySOf5YbwBDkbhLAbOA14raZLdgmCMiPwBuA+rtEzNI+1Kd8UQjCLyxEBBHnzx92NjwQQZhGnt1anGGkCjhEZM3TCV55c/z47DO2iQ0JjCxHN5bk5zHrn0VJejGQzVh027D/Hxkq1cc8YJRcpGVXn4h4c5ePQgr/Z9tfyVjb9PZuotJSuXmITAG13vtcH/+vD1I85WUaiughIIRuGMBb6046rHlTobvgCCmVezCEURzXwADvktbZUU7fb7AvjpfecouBCVjtf8901eS2uQxozfZrByz0oKsUzu3dk7SGj6KRN/hr+e3YYW9RJDOo/BEGn8Z9YvxEVHcev5JwYU47yk9SW0rVvO5STDZcm4XRP6POi8EuO9OQ4DpSZ+AojInVh+myXABOBbYBNwgarO8Zv7DtBXVRuHTcpKokyJn+EgPQUCsl7KQHILuGPV8R8HOPvDs8nMzQwY17wU+tZ+lqev6hyW8xgMVZFV27IY8Oy33Hb+iZx84gbHSuth9df4r1icf7/l0HfDyZL544TQVznCtFLilvgZVB6Oqv5HRH7C6rb5NlbpGAX+JiKnYeXl1MHqO3OOPc9QVpKbl+zwC5asLZCTBb/MOO4fUVZuluO4xGTx6fJt3HjuH2jXpPbxy2wwVEGenLGOlMQY/npOG674/DbXYpwhKxynCzwELp19epP7MUK1ZEqipJWYMBB04qeqzhaRU7Eal10F9AYG2w8vR4BngX+EUcaah5tp6xZkIB5Ql8aqT7YDzTsWDlnG5Ta3YILGiY05GiVc+ty35OYXmhprhmrFlGXbeHT6GvYcOkqd+Gjm/LzLNZIz5GKcjj6ZWy0/rv/SGQoxSUBh8D6ZKkhI/XDU4ktVHaGqrYCGWE3Q+gCnAQ1V9XY7gMBQVtKutMzh5BaAWH//OAEufsI57LrbCOfx8+63sqX8Y++9oY4hlNsZ3XU08Z7i0WjxnnjOrn8tBQo5+YUosC0zm3smr2TKshrf/cEQ4Vg101aw55DlJzmQk889k1dSK9q5O33ICZ6OPplcOHrIeX7eEefrQgQoGi/H1fFTVfcCe8Mki8GXUIIM0q6Elmc4j899zPkYWVtg2ijIDy7QwDeYwGvpDGo7iOlzm5FfWPyfxltjzVg5hkjGqplWPC8lO+8o0TkScOUsU4JnqOHGyc0jxpJxI5jSNn1U9euyHFxELlDV2WXZ1+CC2w/Obbwkf1C+391VKfWUvPWgCgoLuOKzK/h++/dszzwVq1tEcUyNNUOk4/Qbjq3/DUTvZVi7YczbOq8oktNqzR6k/0YVfnge18CghHrW/2Y5RotVFsFYOF+JyHzgaeBLVTdngYWIxAADgNuBM/Grv2aoYNz8QQFrxDZB3HV5ojyM6jKKUXNH0SB1BbszugTMMTXWDJFO3aRY9h0+FnYcFZdBbMPZRGd35t4z7uVe7g3+YL7BATGJkHcYmnaFXT8Xv/GLSbCWzqHKFPUNJ8EonC5YymYasFtEZgMLgY1YTc4EqAe0Bc7A8uekADOBzuEUVkQuwoqA8wCvqupYv+3/h9VxNB/YDVyvqr/Z2wqwGskB/K6qA8MpW5XFrfrBcSZ59W7Rm04NO7EpejYJe9PIPnrMyomNNjXWDJGNqhIfHYXgtUPyiW/6ERQk8Pdu94R2MP/ggLzDEBUNp98EIu6KpRooGH9KVTiqugqrPfOZwN+AS4FhBNqDAhwAJgMvqOoiwoiIeIDngQuBrcAiEZmmqmt8pi0DuqvqERG5GSt36Cp7W7aqdg6nTBGD23Kbv+UTHR+02S4i3N71dq6bcR2XnLWJBcvS2J6ZjQg0S4nn0s5NwyS8wVDxzFi9k+1ZOZzT+XeWH/7Aaj0gcGb9Qfz5tPahHcytiOacf1p5ctVQsbgRSlj0D8AP9oW/G9AeK0pNsayJVcAyVf/qb2GjB7BBVTcBiMiHWMqvSOGo6lyf+QuwKlkbnPC3fABqNYEOg9338aN7k+70ataLBbs/5sv/u546sXV454fNPDB1Nd9u2MPZbRuWg+AGQ/lSWKg8M/sXmjZbw8/5H6HROUUFa5dlfsn0TacH768pLHT3oVZikd7KIuQoNduHs9B+VCTNAN9vbitWSLYbfwG+9HkdLyKLsZbbxqrqFKedRGQkMBKgZcuWxyNv1cfX8ln9KXw8Aub/B3oHn0Y1uutohnw2hH6T+nE47zCNE5vQoEkfxs1KodeJDRBxakllMFRdvlq9g7U7DtI87Sty8kJM8PT11dRJhXjnEGogrDXKIoWg8nBEpKXPo7mIVOmUchG5GugOPOkzfIJdauFPwDN2kdEAVPVlVe2uqt0bNqxBd+gdBkHHK+F/T8C2JUHvtjFzI1ESxaG8QyjKjiMZFNT7mBWZc/nfL7vLUWCDIfx4rZs/NEziQJ7z79c1wdO/TcCB7bBrNbTs6ZwnVw2izkIl2MTPzcCv9uM3IFNEdojIWyLSqbyE82Mb4NvhqLk9VgwRuQCrWvVAVc31jqvqNvvvJqxmcoGhVTWdS56E2k3g/aHwdIegEkLHLx1Pod8qar7mkth4JuNm/UJNbwJriCymr8zgl52HuP2Ck2iY6HzD6Zrg6eSrAUsBRXjCZrgIdkntG44FCQiQDJwEXAMME5Exqlre9dMWAW1FpDWWohmKZa0UISJdgJeAi1R1l894XeCIquaKSAPgLI41kTN4SUiBtKvg26ePjZWSEOp2t6fRmfy0NYs5a3fR55SIr+NqqAEU2NbNSY1r0b9jKh9va8auI7uKzSkxwdPNJ5O1NeITNsNFsMU7e/uPiUg0cBnWhftpEdmrqu+GVbriMuSLyK3ADKyw6NdVdbWIPAIsVtVpWEtotYCPbd+BN/z5FOAlESnEsurG+kW3Gbys/DhwzLfrn18Ip1uNtdSkJsQkxXDTu0vIL1BTY81QZZmybBtPzljHNjvRc0TPE1i5dwXLdi+jd/PerNu/rvQEz52rISoKCh3SFGugr8aNoNoTlHgAkfpYHUHrAk2rUx21SmtPUJmU1BohOiEgSW36WX8lfetXxarnCsJlze/ggzmNySs4dqyEGA+PD+5olI6hymDVS1tJdt4xRZEQAyd0fJ2j7OezQZ+RFJMUuKNvcEBiPcg5CDHxkJ9r1UPzUtY2ARGOW3uCkIp3OmHXU3sIK/kz4ASGCKOkuzGHUjj9l31Kes/0oj7udeProiifr/uxmLKBYzXWDIaqglUvrbhVkpe4hG3Zv3B7t9vdlY1vcMCRvaD5cP4DcOlzxldTAkEtqdlLUcGYQt/aS1mqqsdVGNRQSZShFI63xpqXpxc/zRur3yC6dkvyDxaPKTE11gxViYDfY1QucY2+oiC7BQPaDHDeySk4QAvh+wk1LpEzVIJVCl9QssKpj5UTswgrCdQQqYShFM5tXW/j7WVziW/6EVowHYk+gOalkLu7H42jepbzGzAYgqdpSkKR7wYgtv4coqIPkph5A1HisgBUUnCAoUSCDRpwUfUWInIHlsIZqaorwiGYoRIJthQOQIN28NOHMOfRIgUV0+dB+p1wEdO3v0hUzAEAJDaT+NTJ9D3hhAp4AwZDcNx54Un8Y8abxDacgcRkAlB4pBX3nN/PeYc9G0CinBsemuCAUjnuZS8R6QA8AKwzyqYaE2D5NLOUzcbZsGmOtaQARWHUy1qfiH+RAYnK47t972BF0xsMlc+qA/OIS52MRB2LdYpJ2k5M8nKs4iY+bF0C7w+xag5qvhUgULRTzUzkDJUyKxwRaY515bgHiMEuB2OoxjhZPmNPgJzM4mN52ew4mkmAxgEyQm3DazCUE4WFyie/vox4igfW5uvRY+VrfKPRwOpVc9N8qxpHNWwfUN4EGzSwyfclUAerBYEAWcAwVf027NIZqj45WY7DTfILyIgJ/HlFF9ZFVU2NNUOl89XqHRRE7cfpl7jj8I7AtgJgtRbYtsQkcpaRYMOiW/k8Wtr7LQYeAU5W1enlIJshEnBZtx6d6yHeE19sLAoPBzMu5NsNeypCMoPBlcJCZcLX6/EUOhfXbJLUxDkaLT/nWBK0IWSCUjiqGuXz8KhqXVU9XVXTVXVneQtpqML0eTCwMKFE0f/MfxTLz0mMTqSQAurFNmH87PWmxpqhUpm5ZidrdxzkjMbnBWwrKl9j2gqEneNO/DTUcNKuLF6YMLG+FUCwejL9W17AzCtmsmL4CuZcOYemSU1JaDqJxb/v4oeNeytbckMNRdWybk5oCL8c+obGiY1pktQEQUhNSiW9Zzr9o5LdD2Ci0cqMSc40HD/+69lL34Fpt8LrF8HhXZC1jaTk5jx02tXcuPE9ajeZw4g34sgrKDQ11gwVhn/NtJM7fsbO3Eze6/8e7ev7dPHMWAFv9rcaEuZkBZRzMtFoZadUhSMic8JwHlXVPmE4jiES6HoNbPofrPIpBJq1hZ7/e4az25zLN8lziUlaRGz0ITLzUrh35sXAcKN0DOWGf8206Fqr2Zb/Hec1+jPtt62Ct6+0lspqN4HcQxBfB/4yE3773kSjhZFgLJwogitrUxImJKmmsWVB4FheNt22LeebBrFExRwCrIRQbTSJx/4XzWVd7qpgIQ01hSdnrKNj4nvsbbiU3dHW5Sj2aAIpSw7DUp9ItIN25fNzx1gKxkSjhZVSFY5TawKDoVRcHKsTU6IcE0KPJH0GGIVjKB+a6xtsarKUnKhjbmuNzqZT7FtwxKG+38JX4CyXvjeGMhNsi+mnRcQUwTIEj4tjdUe0x3E8yi4rYjCUB3saFlc2AEejhJfqOVSDBhOJVk4EG6V2O/C1iFxcjrIYqhNO4dJAkwKHGlRAcmyj8pbIUIPxLqP543YDZCLRyodQwqLjgMkiUlohz64ictPxieV67ItEZJ2IbBCRux22x4nIRHv7jyLSymfbPfb4OhFxqcxnCBv+4dLJLeCifzM6J5r4wsJiUz2q3NO0p5XZPe5UqwncuFOt1wZDqPj+jp7uAG8MoEm+841OkwINvDEykWjlRihh0U8C1wOTRORKu6WzE6cAzwMvHq9wvoiIxz7uhcBWYJGITPNrFf0XYL+qnigiQ4EngKtEpD0wFOgANAVmi8hJqk4lXw1hw8Hh2v/7CXBkH+PrprAj2kNioXLYE0Xikvcg58gx561dBLToOAZDMPiXozmwFQ5s5cSGLcmIKT41vlAZ/YfBcFpHE4lWQYSicH4CzgPmAB+JyFBVnVIuUjnTA9igqpsARORD4FLAV+FcCqTbzycBz4lVtOtS4ENVzQV+FZEN9vF+qCDZDV4ObKM/Sv/DRwA4Cgxr2oSHU+LptG0/9Xzn5mVbFwLzz28IFodyND/HxvBDktJCGpKfv5sdUdCkEEa3GUT/3v+0JpnfWIUQUuKnqq4SEa/SmSgiw1R1cvmIFkAzwLfWxFasHjyOc1Q1X0SysJrDNQMW+O3rmPQhIiOxK1+3bNkyLIIbfEhuXqxkSCzw+O69DG3WhPQG9Ri/a0/xGHrjvDWEgt/vJVuEuxs2oG5BIa9dPonU2vUrSTADlKG0jaquxrJ09gMfisiQsEtViajqy6raXVW7N2zYsLLFqX44BBO01WhG7stmblIiZ7dsRlqrFvRt3pTpSYnGeWsIjeTmTE9KpG/zpqS1akHvls3YFBvDfVlqlE0VoEy11Gy/SW9gL/CeiFwVTqFc2Aa08Hnd3B5znCMi0UCyLWMw+xoqAodgAhk4gZ31LyZKlSyPBxUhIyaa9Ab1mN5lUGVLbIggprfqRnqDemTERKMiHImKwqPKgU6XV7ZoBo6jeKeqrsVSOruBd0RkWLiEcmER0FZEWotILFYQgH/gwjRguP38CmCOWmWJpwFD7Si21kBbYGE5y2twI+1KuGMVpGdaf9OuZH7iOgr9MkJzoqIYnzG3cmQ0RB7ZmYzfvyQg36ZAhBf2m3/3qsBxFe9U1XUi0huYC7wtIuVWfdr2ydwKzAA8wOuqulpEHgEW21Fzr2Epvw3APiylhD3vI6wAg3zgFhOhVrXYecS5y8WO3Ew4uMOqcWUwlMTM+9nhcgUynWarBsEqnLnAYacNqrreR+m8iaUQygVV/QL4wm/sQZ/nOYCjT0lVHwMeKy/ZDMeH5Keg0fsDxhvnF8DUW+DPkxxbVhsMAGycA8veoXaLEzkQfTRgs+SnVLxMhgCCUjilVXpW1Q0ici6W0rmE4y/2aahhZO/sS1zqZCTqWH95VaHZno6QMwP+3RqyM4vnSfj2mzf5EzUP3+9foiCpMft29sfT9NNi9yZaGEPOzr6VJ6ehiLD1w1HVTT6WTouSZxsMxWkU1ZOdGRDXcAYSk4kWJBAVnQ2xiZAbBdm29eNNCP19Afz0vkkUran4J3hqAeRk0qnOOlYLFOYnIZ7DaF4Kubv70TjKlIKsCoS1AZuq/ioiXQDz7RpCYky/dtwz+SiHN3axRwpJavUy6+suZs9haODrccvLhsWvE2BIm0TRmoNDgudWyWdD8mryD6SRve3PReMJMR7GDG5X0RIaHAi7k19V96vq9HAf11C9uaxLMx4f3JFmKQl24mcUJ0VdTw7wr/r1HPZwWbU1iaLVD/8aewteLJY8DNav4fH69fCokrNzALXiohGgWUoCjw/uaJr7VRFMx09DleGyLs2KLgwPTl3F+z/+zpjmynNJiZzTshmZUVE0yS9g9P5M+h/JtZZR/DGJotUL/6WzrC3w1T8Cps1JTOCbxAT+tj+fZ6PqMf+u86ibFFvBwhpKw3T8NFRJbjr3D3yw8Hd+jz8H0fns91hl5K2E0PpQuyP9180NWFbhjJsrQVpDWHAKAplxX+B3DBCXzPQ4YXydRHZEexCsiMY1uwdyba9WRtlUUUzHT0OVpGlKAld2b8GUvauI8vgnhArj2Uv/P044doGq3cSKYlv5MZx2A0THVY7ghrLhZMlMHonbve706HzSGzQhR62oRgX2eGL4uk4y889uXTEyG0Km3BI1DYbj5ebef0CiMx237Ti8o3jFgjvXwuWvwvZlMDu9IsU0hAOHIABQK9zZgfH16xUpGy8FUkjtprOoX8vcbFRVwhqlZjCEk+Z1E0mMakC27gnY1iTJofLAKQPg9JtgwX+tO+Yje01+TqTgFuyhhVaxV19lFJPADo/zKn124d5yEM4QLspk4YjISSIySERuFJGR9vO24RbOYLi1821oYfHOWbFRsYzuOtp5h9ROgMCRPYAey88x3UOrNskuUWTJLQI7x/5xAk2SUh2nO96IGKoMQVs4InIKcBNWUUzvt+q9zVB7zk7gI+AlVf05jHIaaijXpg3mvR9/YyuTkZhMBKGWpwGXtL7EeYe5/8Lk50QgLc6ArEnFx7ytnh06x46ulcQ98+9Bfb7raIlzvxExVAmCCYv+A1ar5kFANjAfeAnYiFX6X4B6wInAGcANwG0iMhn4h7dDp8FQFqYs28bm304mN/9uAGJSfmRf6qc88c1k7j7XoeS829KMyc+puhzcCb/MgManQk5WUKWK1m0TFKUwPwHxZFsVBfZdTF5W54qV3RASwVg4a4CVwAhgsqo6FvH0IiJJWFbQaHvf+OOU0VCDeXLGOnLzC4te52V2J7b+N3yw/iXuOmcQUf5OZb+OosfGTeJfleXrRyA/B4a8BQ1OLHW6qvL22hcplDoc3jgG9NiS65Mz1pkkzypMMD6cIXYHzHdKUzYAqnpYVd9S1a5ARTRmM1Rjtmf6Ry55yN3dl8KYDL749YvAHRw6igKQ1BgK8stFRsNxsHUJLH8XzvxbUMoG4Lvt31EQ+ytH95xfTNmA0+/FUJUoVeHYfWbKhKpOLeu+BgNY+Tj+5B/oSFReM55f9jx5BcVDY506inLqENi+BD69EX76sHiZFBNMUHkUFsKXd0GtxnDOmKB2UVWeW/Yckl+PvMzuAdudfi+GqkMoQQO1gFOx2jYfBDapqulqZChXrKKeK8nOO1bGxhPl4c9tb+SdzQ8yef1krjrZz5B2cDLTpIOVn7P602MlcUyF6cqhqKKAvfTZ7TqIqx3UrnO3zGX13tVc0OQ2Pl1f/PKVEONhTD9TpLMqU6qFY7dlfgGrlfR3WA3Q5gPbRGSbiLwtIv1Fyq87lojUE5FZIrLe/lvXYU5nEflBRFaLyAoRucpn25si8quILLcfnctLVkN48S/qmRjrAVX+nNaPE2qfwL8W/ou0t9LoO6kv0zeVUDO21x0QlxxYf80bwWaoGLwVBXz9bCs+DMrSLNRCnlv+HCfUOQE92I3oKGhSJ94U6YwgRLXkMmki8ixwC1Ynz3nAUayw6DHAfqAuVhzqWmCUqn4ddiFF/g3sU9WxInI3UFdV/+E35ySsIqHrRaQpsAQ4RVUzReRN4HNVnRRw8BLo3r27Ll68OEzvwhAOtuw7wvn/mUfPTr+xKvdVjhYe6+4Y74knvWc6/dv0d945PQXnUiliVSswlD/jTnUJ6mhhVY1wYPqm6YxfOp6MwxkAXNRiMJ/M6sENZ7fh3ktOKU9pDWVERJaoasCaZzBBA1cBr6vqxar6hKqOA/5tbxsCtAHuBhKAGSJyR7iE9uFS4C37+VvAZf4TVPUXVV1vP98O7AIaloMshkqkRb1EruzegsUH3i+mbAByCnIYv3S8+85ulaSTmwWWwPfecbuNG8pGiGHr0zdNJ/379CJlAzBzy+fEp/zEyHPalIeEhnIkGIWTAPzgtlFVN6vqk0A7YBzwlIhcGCb5vDRWVe8vbgfQuKTJItIDiMXKFfLymL3UNk5EXIst2ZUTFovI4t27dx+34Ibwc8t5J5ZcY80Ntwi2nEMw7Vb7ztunOsHn/+ez/GOqFoSMv7Kefqf73OTmTN80nb6T+hZbIh2/dDw5BTnFphZylNpNZ9HA1EyLOIJZUpsLbFPVq33G6mP5dC5Q1Tl+8ycD9UKtMi0iszlWwcCX+4C3VDXFZ+5+VQ3w49jbUrGW/oar6gKfsR1YSuhlYKOqlrpwb5bUqi493j7PscZaalIqM6+Y6b6jfwn89pfCghece+tIlFXLy58Sln8MNv7Vn70kNoSjB628Gy8xCUw/66+kb/2qmHKJjYoNsGK9CMKK4SvKQ3JDGHBbUgsmSm0s8KV9gHFBzP8Cy9IJCVW9wG2biOwUkVRVzbCVxy6XeXWA6cB9XmVjH9trHeWKyBvA30OVz1C1GN11FI8v+icSVTws+k8n/6nkHZ0i2H543nmuk7IBU7XAH6c+No7Vn4GYOLjoXwHzx//yaoAl46ZswNRMi1SCycOZgRUg8JSILBSRq4GSYhjPBY6EST4v04Dh9vPhQEB+j4jEAp8Cb/sHB9hKCjuS7jLA3J5GOH8+dRCtdDiFR1NQBSlIJkbimLF5BnmFeaUfwBc3345LaXxqNTK+HS/Fos7sZccpNzsHBgBkbSveVuKOVZB2ZYlLoTF+K+CmZlrkElQejqr+R0R+AsYDbwP5WOE+fxOR07DycuoA/YBz7HnhZCzwkYj8BfgNuBJARLoDN6nqDfbYOUB9ERlh7zdCVZcD74lIQ6y6b8uxipAaIpgpy7axYdNJ5Ng11gASU1azSt/hpZ9e4tYutwZ/sD4PBi7/xCRApz/BT+8H3qkf2gmf3lTz8nn8LZned1u5Tf6fT2EJFR1clHuTpCbFAgOKpsc0Yt+WPki9L5GYTFMzLcIp1YdTbLJlIVyEFbnWG2jpN+UI8CowRlVDvM2sehgfTtXlrLFz2OZQxqTeCZMpSFrMWxe9RedGnYM/oNOyUNqVgePn/N1qe3z0UOAxqotvx+mzAGefTEk49LHhjxMclfLnGz/n3m/vLVb9Od4TD3uGsHtHh4D5zVIS+O7u84OXxVChuPlwQlI4DgetjxUWXQs4AKxR1WpTzMgonKpL67unO2fUROXQPO0ZDh49SKEW0iSpCaO7jnbPzSkL1Tmfx8nZH2UvhDhZLiUFVnh9OUFUf16ycwkjvhpBndg6HDx6sOh7u/Vl10+aX8eG8Ts1hJXjCRpwRVX3YrUoMBgqlKYpCY4WToNGGzmSd4QCe7kr43AG6d+nA4RP6bhVpK5VYrR+ZODk7C9picylI6dbHxs3Xl/1OnXj6jLjihkkRB8LXf9Xytdsy8wJmG9qpkUmwZS26VPWg4uIa+SZwXA8jOnXjoQYT7GxKIG4RjNCTwgNFbd8niN74bPbIzeYoLDQ3dnvhktHzlB8Wev3r+ebrd/wp1P+VEzZAPyxU9OA+aZmWuQSjIXzlYjMB54GvlR1Slg4hojEAAOA24EzsXJfDIaw4q2Z9eSMdWzPzKZ2fDQHcvLJynNO1i0xITRUvBdT3+WinqNg0auw5I1j86pyMIG/n6bXHfDLV+7zE+pBfvZxWzJOvLn6TRKiExh28rBi46rKdxv2kpIQQ2Ksh4ysHJqmJDCmXztTMy1CCUbhdMFSNtOA3XaC5kKsLP59HOv42Rar42cfIAWYCXQOu8QGg81lXZoVXXhy8gro/eQ8cgvrkh+1L2CuRzxsObCFFnVahOfkThfZ7x2sKN/ioEH6M8odfz9N1haY/n8gHug0DNZMCVQsFz9hPQ/ze8g4lMEXm75g6MlDSY5LLrbty1U7WLkti6eGdOKKbi6h64aIolSFo6qrgL4icibwN6y6ZsMI9OUJVuDAZOAFVV0UZlkNBlfiYzyM6tOWB7++kDotppBXmFu0LTYqFo94GPbFMIa0HcL0X6ez4/CO8AcUZG1zGd8C0247ll1f2ZaPW1JmUgMY9CL84Xx3xRJmed9e8zYA17a/tth4fkEhT81cR9tGtRhkrJlqQ9BBA6r6A/CDiHiAbkB7rOKYilXmZhWwTNUtPdtgKF+GdG/Oy9+cRcGBeBIazyymVNIapDH8q+G8uurVovlhDyhwCyaA4qVcoOyWT7Dh227jve92l/GQXcDjOJfIgiUrN4tP1n/Cxa0vJrVWarFtk5duY9Puw7x4dTc8UeXW+cRQwYQcpWb7cBbaD4OhyhDjieL/+rZj1AdHeKbnsIB1/iiHygHegIKwKBy3BFK33JWsLTBtlOUb8b72Wj5Qei6Md/7vC4onqJY0PvUWd/ndKi6UEx+s/YDs/GyuO/W6YuM5eQU8M/sXOrVIoV+HahD5ZygiKIUjIr4JnoVAlqoeLB+RDIayM6BjKmO/+Jk7P/6JOyYuL+Zk3nXEsQRf+AIKnIIJinJR3CwfP2WUlw1f/qO4g967JCeeQOWVlw2LXws8rts4QGxt0HznAIAKYPqm6Tyz5Bl2HNlBnCeOX/b/Qtu6bZmybBtPzlhXFO5+WZdmlGNfR0MlEKyFsxk/n42I7MZqyva0qv4UZrkMhjIx7aft7Dl0lIJC6+e6LTObeyavBNzLp4S1EKTbcpS/5ROdEKhsvGQHBj0ELMkdD0cPweCXKyWIwdvfxluoM7cgl/Tv01m8eR8fzm1YrJX4G99t5qTGtU1EWjUimH44AN/4POYDK7AKeF4DLBIRU0nPUCV4csY6jhYUdyNm5xXw5Ix1jO462iqX4kOURDGqyyjKlbQrA3NVBnpfhwHxhDae3NyxgGZF4NTfJqcgh09+faWYsoFj35uh+hBs8c7e/mMiEo1VefnfwNMisldV3w2rdAZDiGx3qD7gHff6acYvHc+OwzuoHVubA0cPsCc7sK9O2AnW8olJsKwfJyvHLRfGqchoSeMVtHTmj6o6WpgAhZ79juNu36chMilzaRtVzQcm2Q3afgDGicjE6lC00xC5uJW88ZZC6d+mf5HiUVXu/N+dPLP0GdIaptG1cdcKldXV5wPOiqikXJiWZ4Q2XgF4O3buOLyDxomNqRdfz3VuVIFjP0VTwqaacVzFO4sOIjIMeBfoZYdPVwtM8c7IY8qybdwzeWWx5RlPlPCfIZ0cfQEHjx5k6OdD2Z+zn8SYRHYd2VU+BT9DxS3MOULw99V46dygM2v3ry02Hu+J56yUm5jybXFfWkKMh8cHdzQ+nAjkuIp3ikghzkVb/fnWjipRVT2uwqAGQ1nwL3mTGOfhcG4B7Zo49wysHVuby068jAnLJnAwzwq8LJeCn6FSQbkw5YWTrwZgZ/ZO0numF1k+XuX+4dwGJMTsIyUxlh2mhE21JSgLR0Q+p2SFUx84HViElQSKqv4xHALa568HTARaYUXMXamqAYu+IlIArLRf/q6qA+3x1sCHtpxLgGtU1b1/rY2xcCKfzCNHOffJeaQ1T+bt63s4htn2ndTX0beQmpTKzCtmVoSY1Y60t9KK9bbxIggrhq8oNjZv3S5GvLGI+/ufwg1nt6koEQ3liJuFE1SUmqoOUNU/uj2Aj+2pI33GwsndwNeq2hb42n7tRLaqdrYfA33GnwDGqeqJwH7gL2GWz1BFSUmMZVSftsxfv4d5v4RW2DOsBT9rGG7+Gv8Q9PyCQv71xc+cUD+Ra89sVQGSGSqTYMOiXRGRDsADwDpVXVHa/DJyKfCW/fwtrOi4oLC7lJ4PTCrL/obI55ozTqBV/UQem/4z+X4h0+CehxPW/JwaxPr96zmUdwihuDUZ74lndNfiGRQTF2/hl52HuOfik4mNPu7LkaGKU+ZvWESai8g9WBFqCcDIsEkVSGNV9a557ADc6l3Ei8hiEVkgIpfZY/WBTDuqDmAr4LowLCIj7WMs3r3b+Y7YEFnERkdx98WnsGHXIT5cFJjx75SfEy3RARdHQ+lsO7SNm2bdRHJsMmNOG0NqUiqCkJqUSnrP9GI+sYM5eYyb9Qs9WtWjXwej3GsCwQYNbPJ9CdTBakEgQBYwTFW/PR5B7LYHTr+6+3xfqKqKiJs/6QRV3SYibYA5IrLSli9oVPVl4GWwfDih7GuouvTr0Jgererx+Bc/8/zcDX6O6eL5OfHR8eTk59AquVXlCh0h+IY/R0kU0VHRfND/A9rWbcs17a8JmO9fwubaM+ubEjY1hGAjyVr5PFesNgSLgS+xWhHsPF5BVNW1O6iI7BSRVFXNEJFUwLEolqpus/9uEpF5WL18PgFSRCTatnKaAy515A3VFRGhV9v6LNy8j8NHrZBp37I3l3U5lp+TlZvF4GmDuW/+fUz840TiPHGVJndVxz/8uUAL8KinqD6aP05h6y/M20TLekkmIq0GEGzQQJTPw6OqdVX1dFVND4eyCYJpwHD7+XBgqv8EEakrInH28wbAWcAatcLw5gJXlLS/ofozcdHWgDGn8inJcck83PNhNmZt5Pnlz1eUeBGJU/jz0cKjri29n5yxzpSwqcFEipduLHChiKwHLrBfIyLdRcTb4OQUYLGI/ISlYMaq6hp72z+A/xORDVg+HZcyuobqTEllb/zp1awXl7e9nDdWvUHvib1JeyuNvpP6Mn3T9PIWM6IINcIvlO/AUP0oVeGIyN9EJORGGSKSKiI3ijg0IQkRVd2rqn1Uta2qXqCq++zxxap6g/38e1XtqKqd7L+v+ey/SVV7qOqJqjpEVXPdzmWovriVSXEb79igI4KwN2cvihYlhBqlcwz/ttBe3CL8UlPiHcdNCZuaQTDKYDzwmx25db+IpLlNFJFTReQ+EfkRKxrsOcClZK3BULGM6deOhJjiP8f4mCjG9GvnOP+lFS8FJC96G7ZN3zSdvpP61mjLZ/We1RzMPUiU32XEKfzZS/cTAmumJcR4XL8DQ/UiGIXTGLge+B0r4XKZiPwqIs+IyPkicq6IPC0iG4GfgHuBDOAGINUU8zRUFS7r0ozHB3ekWUpCUYbIuW0bujqr3ZaFMg5ncN+395FxOKPGWj57s/cyeu5oGic15r4z7isx/NlLRlY2s3/eRfvU2jRLiUeAZikJpl5aDSKk4p0iEg9ciJU42R9ohBW1thf4HMsZP1NVq8WCrCltU7257YNlzFi9g5m3n0OrBkkB291K3gjiWLalppTCySvM468z/8qqPat4++K3aV+/fVD73fLeUmb/vJNZd5xLy/qJ5SyloTI5rtI2XlQ1R1U/U9W/AKlAL+BcoImqXq+qU6uLsjFUf+7vfwqxnigemrYapxsvp4TQeE+8o7KB6l8Kx7uM2PWdrizZuYRLT7w0aGXzzS+7mb4yg1vOO9EomxpMmR36avG9qn6rqoH1QgyGKk7jOvHcceFJ/O+X3cxYHags+rfpT3rP9IDlotSkVMfjVedSON58G1+Lb9qGaUEtI+bmF/DQtNW0qp/IyHNMcc6aTFj64VRXzJJa9Se/oJABz37LjqxsEmKjgyqN79TrJc4Tx8M9H67cHjrlSFkqavtXFLjx3Dbcc/Ep5SqnoWoQliU1g6G6Ee2Jom+HxmRm55ORlYNyrALBlGXOBSn8LR9B+EPyH6qtsjlw9IBra2i3ZURvRQHf7qtvf/+b62dqqBmYJmmGGs8nSwIvgt7sdzcrx7dV9Ys/vcjzy59n5e6VdGzYsVxlrQh8a6PVj69PflHd20DclhFLqihgItJqLsbCMdR4jjf7/Zr211A3ri4Tlk0Ip1iVgq+vRlH25OwhMzeTPi36OAZQuOXbmIoCBieMwjHUeEKtQOBPUkwSN3S8gQUZC/gx48dwilbhuLWGXrNvjWMAhdsyYuM6pqKAIRCjcAw1HqcKBHHR7hUInLjq5KtonNiYCUsnOIZYRwol1Ubr36Y/M6+YyYrhK5h5xUxXZaOq1EuKCRg3FQUMRuEYajz+FQiiBBrWiqV/mnP4sxNxnjhu7nQzK/as4NyJ5waUvImEUjiLdixy3RZKyPcnS7exJuMgl3VuWvSZmooCBjBBAwYDYCkd78Xwq1UZ3PTuUl6Yt5FRfQJ7urgRExWDIOzP3Q9QVPJm6c6lTN04ldyC3GLjXrwO+iZJTRjddXRYo918AwB8j+8/fmbTM5m2YRr1E+pz8OjBIlmhZF+NPxlZ2Tz82Wp6tKrH01d2JirKNFYzHMPk4ZSAycOpuYz+cBnTV2Qw9daz6NDUuSKyP265Km7ESAwqSn7hsSiweE98ib4RN5wUCxCQLxTviefSEy9l6oapAb6atiltefPiN5m/dX6ZlKCqMuKNRSz8dR9fjj7bsVyQoWbglodjFE4JGIVTc9l/+Ch9n/mGGPsOPSOIhNC0t9Jcy96EQknJlMEqlmiJRkTIKwy+dm5Za8H5J3gO7tqMp6/sHPJxDNUHN4VjltQMBgfqJsVyaaemvPrtr0VjxVtSByqdJklNHC2cKImiMITqT26Oe/8KBxmHM3jwuweJkqgAayVf8wlV95WlFpxTy+gvV2ZwTglVuA01l4gIGhCReiIyS0TW238DmmqIyHkistznkSMil9nb3rRbKni3da7o92CIPL5cFXgBLqkdsluxzyEnDXEcT4lLcTxOg4QGjuNu7ZydwphLIsqlJ2JZasE5J3gWmpbRBkciQuFg9eH5WlXbAl/br4uhqnNVtbOqdgbOB44AvusDY7zbVXV5BchsiHBCTV50K/Z5/xn3O47f3ePuAEUEkJOfw28HfgsYD8U/BJAcmxySAgw2MMAXk+BpCIVIWVK7FOhtP38LmAf8o4T5VwBfquqR8hXLUJ1pmpJQrBaY77gbviVvghmH4lFqV7a7krdXv82wz4cRHx3Pnuw9NExsSJNEd+sjOTaZ3ILcgOCAe06/J+D43gCALo26hCU6rl5SLHsPHw0YNwmeBiciImhARDJVNcV+LsB+72uX+XOAp1X1c/v1m8CZQC62haSquS77jgRGArRs2bLbb78F3mkaagZO/onoKOGpIZ3K1T/x6spXGb90fMB4h3od2Ji1MUCxpPdMB8o3vNqJjKxsLvjP/zhytKCYuyghxmNybmo4VT5oQERmA063cff5vlBVFRFXLSkiqUBHYIbP8D3ADiAWeBnLOnrEaX9VfdmeQ/fu3au+NjaUG94L5pMz1rE9M5v4GA/ZeQU0qhNXruf9aN1HjuP7cveR3jPdVbFUZLXqo/mF3PLeUgDuvuRk3v7+N7ZnZpcayWeo2USKhbMO6K2qGbZCmaeqjjUyRGQ00EFVR7ps7w38XVUHlHZeExZt8OVwbj6XPv8dmUeOMn3U2a71wo4Xt/BqQVgxfEW5nDNUHv5sNW98t5nn/9Q1pIoMhppBlbdwSmEaMBwYa/+dWsLcYVgWTREikmorKwEuA1aVk5yGakxSXDQvXt2Vgc99x9CXfyA3v5CMzNLzc0LFLby6sjuK+ufbnNO2gVE2hpCIlCi1scCFIrIeuMB+jYh0F5FXvZNEpBXQAvif3/7vichKYCXQAHi0IoQ2VD9ObFSby7s159c9R9ieGVzDtlBxC68uSxRZuHBqqLZw8z7TUM0QEhFh4ajqXqCPw/hi4Aaf15uBgNtMVT2/POUz1Czm/LwrYCyczcW8vpiKDgIoCad8mxw738b4awzBEhEKx2CoSlRE7klJYdSVgcm3MYSDSFlSMxiqDMfbsC3SyC8oJC7a+VJRXd+zoXwwCsdgCBGnhm0AI3q2qnhhyhlV5b5PV5GTX0iMp3irAdNQzRAqRuEYDCHi37CtUe04kmI9vPXDZnYeCK2uWVXnyRnrmLh4C7edfyJPXtHJNFQzHBcRkYdTWZg8HEOwrNiaybCXF1A73moLsCOIdgZVFf/w5zP/UI/3bzgDK6vAYCgdtzwcY+EYDGEgrXkKw3u2YseBXDKyyidcuiJwCn9e/nsmU5dvr0SpDNUFo3AMhjDhdFEuqZ1BVcS0GzCUJ0bhGAxhItJDh1XVsTo2RM57MFRtjMIxGMKEW4hwg9rlW+wzHBzNL+TOj35y3W7Cnw3hwCR+GgxhYky/dgHtDATIPHyURz9fw5erdlSZisrewIDtmdk0SY4nMdbDxt2HueTUJsxdt4vsvGMtsU34syFcGAvHYAgT/uHSzVISePjS9jRJjufVb39lW2Z2lQgm8A0MUCAjK4eNuw/z5zNa8N+ru/H44DQT/mwoF4yFYzCEkcu6NAu4OL8wb1PAvHDWXgsVp8AAgHlr9wDO78FgCAfGwjEYypkdWc7JoJXliI/04AZD5GIsHIOhnGmakuAY/RUXHcWb3//KK9/8Wm6+HV9fTWpKPD3/0KBEOQ2G8sRYOAZDOeNUey06SsjNLyR92ppy8+34+2q2Z+YwaclWGtaODSjGaQIDDBWBsXAMhnLGa7F4LQ2vJfP4lz+z80Busbm+iaL+80uyfHwtGe/8f89Y6+irifFE8cTlaSEd32AIBxFRS01EhgDpwClAD7vxmtO8i4DxgAd4VVW9nUFbAx8C9YElwDWqerS085paaobypPXd03H774v1RHG0oHho8uODOwKBiggICMeOjhLyC52PLsCvY6tOrx1D9cOtllqkWDirgMHAS24TRMQDPA9cCGwFFonINFVdAzwBjFPVD0XkReAvwAvlL7bB4I6bbwcopmzAsnwemLKSvEIlx86R2ZaZzT8+WUGMJyrAknFTNt7zGgyVQUT4cFT1Z1UtrZhTD2CDqm6yrZcPgUvFKnF7PjDJnvcWcFm5CWswBImTb8epz46Xg7kFRcrGS25+IYdy8133cTq+8dUYKouIUDhB0gzY4vN6qz1WH8hU1Xy/cUdEZKSILBaRxbt37y43YQ0Gp0RR7+tw4Hs8k8RpqApUmSU1EZkNNHHYdJ+qTq0oOVT1ZeBlsHw4FXVeQ83ELcnS3yeTEOMhPiaK/UfyAuamJMSQm18YMN8bCGAUjKGqUGUUjqpecJyH2Aa08Hnd3B7bC6SISLRt5XjHDYYqiVtUGzgrovSBHRznG0VjqGpUGYUTBhYBbe2ItG3AUOBPqqoiMhe4AsuvMxyoMIvJYCgLJVkmborFKBhDVSdSwqIHAc8CDYFMYLmq9hORpljhz5fY8y4BnsEKi35dVR+zx9tgKZt6wDLgalXN9T+PPyYs2mAwGELHLSw6IhROZWEUjsFgMISOm8KpTlFqBoPBYKjCGIVjMBgMhgrBKByDwWAwVAhG4RgMBoOhQjBBAyUgIruB38q4ewNgTxjFqWgiXX6I/PcQ6fJD5L+HSJcfKuc9nKCqDf0HjcIpJ0RksVOURqQQ6fJD5L+HSJcfIv89RLr8ULXeg1lSMxgMBkOFYBSOwWAwGCoEo3DKj5crW4DjJNLlh8h/D5EuP0T+e4h0+aEKvQfjwzEYDAZDhWAsHIPBYDBUCEbhGAwGg6FCMAqnHBCRi0RknYhsEJG7K1ueUBCR10Vkl4isqmxZyoKItBCRuSKyRkRWi8joypYpVEQkXkQWishP9nt4uLJlKgsi4hGRZSLyeWXLUhZEZLOIrBSR5SIScVV8RSRFRCaJyFoR+VlEzqx0mYwPJ7yIiAf4BbgQq531ImCYqq6pVMGCRETOAQ4Bb6vqqZUtT6iISCqQqqpLRaQ2sAS4LFI+fwARESBJVQ+JSAzwLTBaVRdUsmghISL/B3QH6qjqgMqWJ1REZDPQXVUjMvFTRN4C5qvqqyISCySqamZlymQsnPDTA9igqptU9ShWH55LK1mmoFHVb4B9lS1HWVHVDFVdaj8/CPwMRFRnMrU4ZL+MsR8RdWcoIs2B/sCrlS1LTUREkoFzgNcAVPVoZSsbMAqnPGgGbPF5vZUIu+BVF0SkFdAF+LGSRQkZezlqObALmKWqkfYengHuAgorWY7jQYGZIrJEREZWtjAh0hrYDbxhL2u+KiJJlS2UUTiGaomI1AI+AW5X1QOVLU+oqGqBqnYGmgM9RCRiljdFZACwS1WXVLYsx0kvVe0KXAzcYi83RwrRQFfgBVXtAhwGKt2fbBRO+NkGtPB53dweM1QQtt/jE+A9VZ1c2fIcD/YyyFzgokoWJRTOAgbaPpAPgfNF5N3KFSl0VHWb/XcX8CnWcnmksBXY6mMZT8JSQJWKUTjhZxHQVkRa2466ocC0SpapxmA73F8DflbVpytbnrIgIg1FJMV+noAVgLK2UoUKAVW9R1Wbq2orrN//HFW9upLFCgkRSbKDTrCXovoCERO5qao7gC0i0s4e6gNUeuBMdGULUN1Q1XwRuRWYAXiA11V1dSWLFTQi8gHQG2ggIluBh1T1tcqVKiTOAq4BVto+EIB7VfWLyhMpZFKBt+yIxyjgI1WNyNDiCKYx8Kl1/0I08L6qflW5IoXMbcB79o3vJuC6SpbHhEUbDAaDoWIwS2oGg8FgqBCMwjEYDAZDhWAUjsFgMBgqBKNwDAaDwVAhGIVjMBgMhgrBKByDIcIQkX4iMk9EDonIbhF5TkTiK1sug6E0jMIxGCIIEbkT+ArIAO4APgNuAcZXplwGQzCYPByDIUIQkQuAmcBdqvqUz/hXwHlAw0isG2eoORgLx2CIAEQkCsuKWQb8x2/zPCAWiJgCn4aaiSltYzBEBv2A9sAIDVyWOGr/Ta5YkQyG0DAKx2CIDK4CCoD5ItLAb1tj++/BihXJYAgN48MxGCIAEfkNaFnKtGaqur0i5DEYyoJROAZDFce2aHZj9WT5r8OUj4BcVU2tUMEMhhAxS2oGQ9Wnjf13karO9t0gIq2BusD7FS6VwRAiJkrNYKj61LL/OvlorrD/TqwgWQyGMmMUjsFQ9fHm1tTxHbQba90MrAOmV7RQBkOoGIVjMFR91gBHsEKjfXkMaAWMUtWCihbKYAgV48MxGKo4qnpERF4FRonIu8D/gIuBQcAYVZ1ZqQIaDEFiotQMhgjAXj77N/BnIBFYAvxLVb+qVMEMhhAwCsdgMBgMFYLx4RgMBoOhQjAKx2AwGAwVglE4BoPBYKgQjMIxGAwGQ4VgFI7BYDAYKgSjcAwGg8FQIRiFYzAYDIYKwSgcg8FgMFQIRuEYDAaDoUL4fytDWhbjQ33yAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "results_zne = []\n",
    "for theta in thetas:\n",
    "    results_zne.append(\n",
    "        mitigated_expectation(theta, backend, readout_qubit, samples, hamiltonian, program, fac)\n",
    "    )\n",
    "\n",
    "_ = plt.figure()\n",
    "_ = plt.plot(thetas, np.cos(thetas), \"o-\", label=\"Ideal landscape\")\n",
    "_ = plt.plot(thetas, results, \"o-\", label=\"Noisy landscape\")\n",
    "_ = plt.plot(thetas, results_zne, \"o-\", label=\"Mitigated landscape\")\n",
    "_ = plt.xlabel(r\"$\\theta$\", fontsize=18)\n",
    "_ = plt.ylabel(r\"$\\langle \\Psi(\\theta) | Z | \\Psi(\\theta) \\rangle$\", fontsize=18)\n",
    "_ = plt.legend()\n",
    "_ = plt.title(\"Mitigated Energy Landscape\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10c88eff",
   "metadata": {},
   "source": [
    "In the energy landscape plot, we can see that the noise has flattened the unmitigated \n",
    "landscape and with error mitigation it has become peaked again. Therefore, we expect \n",
    "the optimization loop to have better convergence with ZNE applied.\n",
    "\n",
    "Run VQE routine with ZNE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "d87720c2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " final_simplex: (array([[3.08437042],\n",
      "       [3.084375  ]]), array([-0.79566667, -0.79533333]))\n",
      "           fun: -0.7956666666666666\n",
      "       message: 'Optimization terminated successfully.'\n",
      "          nfev: 45\n",
      "           nit: 16\n",
      "        status: 0\n",
      "       success: True\n",
      "             x: array([3.08437042])\n"
     ]
    }
   ],
   "source": [
    "res_zne = optimize.minimize(\n",
    "    mitigated_expectation,\n",
    "    init_angle,\n",
    "    args=(backend, readout_qubit, samples, hamiltonian, program, fac),\n",
    "    method=\"Nelder-Mead\",\n",
    "    options={\"xatol\": 1.0e-3, \"fatol\": 1.0e-2},\n",
    ")\n",
    "print(res_zne)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c331fd0c",
   "metadata": {},
   "source": [
    "We can see that the convergence to the minimum energy is enhanced by applying ZNE.\n",
    "\n",
    "\n",
    "## Conclusion\n",
    "\n",
    "While the VQE algorithm is generally considered to be robust to noise \n",
    "[[2]](#references), at the noise level modeled in this example, the \n",
    "accumulation of errors results in loss of accuracy and additional iterations \n",
    "required to reach convergence. Adding ZNE then improves the convergence of the \n",
    "algorithm to the minimum energy. The result is also demonstrated in the energy \n",
    "landscape plot, where the noisy landscape is noticeably flatter than the landscape \n",
    "generated with ZNE.\n",
    "\n",
    "Note: In this example, a small ansatz was used to keep the runtime within acceptable \n",
    "limits. ZNE generally performs better on longer circuits, but there is a tradeoff with \n",
    "execution time.\n",
    "\n",
    "\n",
    "## References\n",
    "[1] Rigetti Computing (2018) Grove (Version 1.7.0) \n",
    "[[Source code].](https://github.com/rigetti/grove/blob/v1.7.0/grove/pyvqe/vqe.py)\n",
    "\n",
    "[2] [[VQE tutorial in pyQuil / Grove].](https://grove-docs.readthedocs.io/en/latest/vqe.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f18a18b2",
   "metadata": {},
   "source": [
    "This final block displays information about Mitiq, installed packages, and Python version/platform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "4533412f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Mitiq: A Python toolkit for implementing error mitigation on quantum computers\n",
      "==============================================================================\n",
      "Authored by: Mitiq team, 2020 & later (https://github.com/unitaryfund/mitiq)\n",
      "\n",
      "Mitiq Version:\t0.13.0dev\n",
      "\n",
      "Core Dependencies\n",
      "-----------------\n",
      "Cirq Version:\t0.13.1\n",
      "NumPy Version:\t1.20.3\n",
      "SciPy Version:\t1.7.3\n",
      "\n",
      "Optional Dependencies\n",
      "---------------------\n",
      "PyQuil Version:\t3.0.1\n",
      "Qiskit Version:\t0.32.1\n",
      "Braket Version:\t1.14.0\n",
      "\n",
      "Python Version:\t3.7.7\n",
      "Platform Info:\tLinux (x86_64)\n"
     ]
    }
   ],
   "source": [
    "mitiq.about()"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "text_representation": {
    "extension": ".myst",
    "format_name": "myst",
    "format_version": 0.13,
    "jupytext_version": "1.10.3"
   }
  },
  "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
}