{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "dac81d40",
   "metadata": {},
   "source": [
    "(examples/pyquil_demo)=\n",
    "# Zero Noise Extrapolation with pyQuil and parametric compilation\n",
    "In this example we demonstrate how to use Zero Noise Extrapolation (ZNE) in Mitiq with a parametrically compiled program in pyQuil. Parametric compilation saves time by compiling the program ansatz once, instead of compiling each time the parameters are updated. In pyQuil, a `memory_map` argument is used to substitute in values for previously-declared Quil variables in the pre-compiled executable. \n",
    "\n",
    "When adding the noise mitigation function `zne.execute_with_zne` the `memory_map` is still passed to the executor function and the circuit parameters are updated without the need for repeated compilation. \n",
    "\n",
    "The example shown here is adapted from the [pyQuil parametric compilation tutorial](https://pyquil-docs.rigetti.com/en/v2.28.1/migration3-declare.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1e508002",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from pyquil import get_qc, Program\n",
    "from pyquil.gates import RX, RY, H, CNOT, MEASURE\n",
    "import mitiq\n",
    "from mitiq import zne\n",
    "from scipy import optimize"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f96ca56a",
   "metadata": {},
   "source": [
    "Use the get_qc command to initialize the simulated noisy device where the PyQuil program will run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5e802eb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize the quantum device\n",
    "qc = get_qc(\"2q-noisy-qvm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e04607e",
   "metadata": {},
   "source": [
    "Set up the quantum circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1f780754",
   "metadata": {},
   "outputs": [],
   "source": [
    "program = Program()\n",
    "theta = program.declare(\"theta\", memory_type=\"REAL\")\n",
    "ro = program.declare(\"ro\", memory_type=\"BIT\", memory_size=1)\n",
    "program += RX(theta, 0)\n",
    "program += RY(theta, 0) \n",
    "program += RX(theta, 1)\n",
    "program += RY(theta, 1)\n",
    "program += H(0)\n",
    "program += CNOT(1,0)\n",
    "program += H(0)\n",
    "program += MEASURE(0, ro[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "380f094e",
   "metadata": {},
   "source": [
    "Mitiq's `zne.execute_with_zne` function takes a quantum program (e.g. PyQuil program) \n",
    "and an executor function as arguments. The executor function must output the expectation value\n",
    "resulting from running the quantum program."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3434165b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def expectation(theta, executable: Program) -> float:\n",
    "    bitstrings = qc.run(\n",
    "        executable.write_memory(region_name=\"theta\", value=theta)\n",
    "    ).readout_data.get(\"ro\")\n",
    "    result = np.mean(bitstrings[:, 0])\n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5a8606f",
   "metadata": {},
   "source": [
    "Run the parameter scan without noise mitigation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "34f2332d",
   "metadata": {},
   "outputs": [],
   "source": [
    "program.wrap_in_numshots_loop(1000)\n",
    "quil_prog = qc.compiler.quil_to_native_quil(program, protoquil=True)\n",
    "\n",
    "thetas = np.linspace(0, 2 * np.pi, 21)\n",
    "results = []\n",
    "for theta in thetas:\n",
    "    results.append(expectation(theta, quil_prog))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85a8f3d5",
   "metadata": {},
   "source": [
    "Plot the energy landscape without noise mitigation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f3de2176",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEeCAYAAABPMvhnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABF80lEQVR4nO3deXxU9dX48c/JJCEbWUkCBELCTgRkx60qm9TaKlqttbWt7ePPtlbtSovda58+2mL3xS52sbWtVYvUuqEIuLKD7IlC2BKWBEISIHtyfn/MjY4hgUwymTvLeb9e82Jy7525Z8LNnPvdRVUxxhhjeiPG7QCMMcaEP0smxhhjes2SiTHGmF6zZGKMMabXLJkYY4zpNUsmxhhjes2SiQk7IvKsiHzC7TiihYgUiIiKSKzbsZjQZcnEBJ2I7BORChFJ9tl2q4is6s7rVfVKVX0owDGpiJwWkVM+j68G8hyBZl/yJpRYMjFu8QCfdzuIDs5X1RSfx48CfQL74jeRypKJccti4Csikt7ZThG5SETWi0iN8+9FPvtWicitzvORIvKSc9wxEfmXs/3XIvLjDu/5pIh80d9AReS7IvKoiPxVRE6KyA4Rmeazf7CI/FtEKkVkr4jc1eG1j4vIwyJSC9wiIoUi8rLzXsudWB92jn9aRO7scP6tInKtnzHPEJHVIlItIodF5FciEu+zX0XkMyLylnPMr0VEnH0eEbnf+X2WAld1eO9bRKTUiX+viHzUZ9//E5Fdzr6dIjLF2b5IRPb4bL+2w/u95sRYIyLFIjLHZ3+aiPzR+RzlIvK/IuLx5/dhgkBV7WGPoD6AfcBcYAnwv862W4FVzvNM4ATwMSAWuMn5OcvZvwq41Xn+T+AbeG+MEoBLnO0zgENAjPPzAKAOyO0iJgVGdrHvu0AD8D68Jap7gTXOvhhgI/BtIB4YDpQC831e2wwscI5NBFYD9zvHXwLUAg87x38IWOtz7vOB40B8J3EVOHHHdrJvKnCB8/srAHYBX+jweZ8C0oF8oBJ4r7PvM0AxMNT5v1jZfh4g2Yl3jHPsIOA85/kNQDkwHRBgJDDMZ99g53dwI3AaGOTsuwVoAb4IxDn7a4BMZ/8TwO+cc+cA64BPu30d26PDNed2APaIvgfvJJPxzpdGNu9OJh8D1nV4zWrgFuf5Kt5JJn8Ffg8M6eQ8u4B5zvM7gGfOEpM6X5LVPg/fhLDc59gioN55PhM40OG97gb+7PPal3325TtfnEk+2x7mnWSSgDdxjnJ+vh/4TRcxF9BFMunk2C8AT3T4vJf4/PwosMh5vgL4jM++Kzokk2rgg0Bih3MsAz7fzWvgDeAa5/kteBO/+Oxf51wHuUCj77nw3lysdPs6tse7H1bNZVyjqtvx3h0v6rBrMLC/w7b9QF4nb/NVvHfB65zqp0/57HsIuNl5fjPwt3OENEVV030ey3z2HfF5XgckOO0fw4DBTlVRtYhUA1/H+yXY7mCHz1alqnWd7VfVBuBfwM0iEoP3i/NccZ9BREaLyFMicsSpXvs/vKUzXx0/U4pPjL4xv/1/oaqn8ZYcPgMcdqrlxjq7hwJ7uojn4yLyhs/vaHyHeMrVyRQ+5xyM9/cb55yr/bW/w1tCMSHEkolx23eA/8e7E8UhvF8ivvLxVqG8i6oeUdX/p6qDgU8DvxGRkc7uh4FrROR8YBywNMCxg/dLd2+HJNRfVd/nG6bP88NApogk+Wwb2uE9HwI+CswB6lR1dQ/iegBvVdUoVU3Fm+Ckm6893CGmfN+dqrpMVefhreIqBv7g7DoIjOj4ZiIyzDnmDrxVlenA9g7x5LW32fic85Dzno3AAJ/fb6qqntfNz2KCxJKJcZWq7sZ7J36Xz+ZngNEi8hERiRWRG/FWLT3V8fUicoOIDHF+PIH3i7vNee8yYD3eO/t/q2p9H3yEdcBJEfmaiCQ6jdfjRWR6Zwer6n5gA/BdEYkXkQuBD3Q4ZrXzGX5M90ol/UQkwecRA/THW213yik5fNaPz/QocJeIDBGRDHxKjiKSKyLXiLdbdyNwyokV4EG8nSqmitdIJ5Ek4/1/qXTe45N4Sya+cpxzxonIDXiT/zOqehh4HvixiKSKSIyIjBCRy/z4PCYILJmYUHAP3i8cAFT1OPB+4Mt4G5+/CrxfVY918trpwFoROQU8ibfOvtRn/0PABLr3pbxF3j3O5GfneoGqtjqxTgL2AsfwfqmmneVlHwUuxPvZ/hdvMm3scMxfnbgf7kbcp4B6n8ds4CvAR4CTeEsF/+rG+7T7A972jy3AJrwdJdrFAF/CW2qoAi7DSVSq+hjwA+AfznmX4m1E34k3Ma4Gjjqf67UO51wLjML7+/sBcL1zHQB8HG9nhZ14bxgex1sqMiFE3l1NaUxkEZFL8X4hD9MQvdjF2525WFW/47Pt48BtqnqJe5EFh4jcgrdDRcR/1khmJRMTsUQkDu/AyAdDKZGIyHSnqiZGRN4LXINPe47TnnI73l5qxoQFSyYmIonIOLxdWAcBP3M1mDMNxNu9+RTwC+CzqroZQETm421bOIq3usiYsGDVXMYYY3rNSibGGGN6LWonnRswYIAWFBS4HYYxxoSVjRs3HlPV7I7bozaZFBQUsGHDBrfDMMaYsCIiHWenAKyayxhjTABYMjHGGNNrlkyMMcb0miUTY4wxvWbJxBhjTK9FbW8uY0x4W7q5nMXLSjhUXc/g9EQWzh/DgsmdLXljgsGSiTEm7CzdXM7dS7ZR39wKQHl1PXcv2QZgCcUlVs1ljAk7i5eVvJ1I2tU3t7J4WYlLERkrmZg+Y9UQpq8cqu58nbOutpu+ZyUT0yfaqyHKq+tR3qmGWLr5jJV3jfHb4PREv7abvmfJxPQJq4YwfenG6UPO2BYj8JV5o12IxoAlE9NHrBrC9JVTjS0s2VROemIsg9ISECA9KY42hdrGFrfDi1rWZmL6xKD0BA5VN5yxXQR+s2o3H50xjLSkOBciM+Hu2//ZzoGqOh657UJmFGYCoKp84s/r+eFzxcwem8PQzCSXo4w+VjIxAaeqDM04s+463hPDiOxkfvRcCRfe9yLffXIHB47XuRChCVdPbC5jyaZy7poz6u1EAiAi3HvdBAS4e8k2bNG/4LNkYgLuFy/uZu3eE1xRlENeeiIC5KUn8qPrJ/LCly7n6bsu4b3nDeThNfu5/P6VfPbhjWzcf8LtsE2I23/8NN98YjvTCzK4Y9bIM/bnpSey6H3jeHX3MR7dcNCFCKNb1C7bO23aNLX1TALvsQ0HWfj4Vq6bksePbzgfEeny2CM1DTy0eh9/X7Of2oYWpuSn8//eM5wrzhvIf7ccsm7F5m1NLW3c8NvX2XvsNM9+4VLyuui11dam3PSHNew8VMsLX7qMgWkJQY408onIRlWddsZ2SyYmUF5+s5JP/WU9FwzP4k+3TCc+tnsF39ONLTy24SB/fG0vB6vqyUyO42RDC82t71ybiXEe7r1ugiWUKHXvs7v43UulPPDRKVw5YdBZj9137DTv/fnLXDxiAA9+YtpZb2iM/7pKJlbNZQJi56Fabv/7JkbmpPDAzVO6nUgAkvvFcsvFhaz6yiwe+OiUMxIJWLfiaPbym5X87qVSPjIz/5yJBKBgQDJfuWIMLxZX8OSWQ0GI0IAlExMAh6rr+eRf1tE/IZa/fHIG/RN61kvLEyNcOWEQLa2dl5atW3H0OXaqkS89uoVROSl866qibr/ukxcXMmloOt99cgfHTjX2YYSmnSUT0ys19c3c8ud11DW28udPTg9IHbWNbjbgbf/48qNbqG1o5pcfmUxivKfbr/XECIuvn8jpxla+8+SOPozStLNkYnqssaWVz/xtI3uPneZ3H5vK2IGpAXnfhfPHkBj37i+OxLgYFs4fE5D3N+HhT6/t5aU3K/nWVeN6dG2Nyu3PXXNG8vTWwzy3/XAfRGh8WTIxPdLWpnzt8a2sLj3Oj66fyEUjBwTsvRdMzuPe6ya8q8fOHbNHWeN7FNleXsMPnytmXlEuN18wrMfv8+nLRlA0KJVvLt1BdV1TACM0HVkyMT2y+PkSlr5xiIXzx3Dt5DPnSeqtBZPzeG3RbNZ9Yw7gTV4mOpxubOGuf24mK7kfP/rgxF71xorzxLD4holU1zVxz1M7Axil6ciSifHbw2v288CqPdw0I5/bLx/Rp+fK6Z/ApKHpLC+u6NPzmNDx3Sd3sPf4aX724UlkJMf3+v3OG5zGZy4bwZJN5awsseuor1gyMX5ZvvMo3/7PdmaPzeH715wXlD78c8flsOVgNRUnz5zry0SW/7xRzmMby7hz1kguGJ4VsPe9c85IRuak8PUl2zjZ0Byw9zXvsGRium3LwWru/Odmzhucxi9vmkysJziXz5xxuQCstNJJRDtwvI5vPrGdqcMyuGvOqIC+d79YD4uvn8jR2gbufbY4oO9tvGzWYHNWvqslikBaYhx/vGUayf2Cd+mMHdifvPRElu+q4Mbp+UE7r+l7vtdXrEeIEfj5hyf1yY3K5PwMPnVxIQ++upf3TxzERSMC12nEWMnEnEXH1RLbFOqaWnl99/GgxiEizBmXwytvVdLQYcEtE746Xl/NrUqbwoZ9fTfp55evGMOwrCQW/XsbdU229kkgWTIxXepstcTGljZXpjWZMy6XhuY2Xt9zLOjnNn2js+uruVX79PpKjPfwww9O5EBVHZ/920Yuvm8FhYue5uL7VtiS0r1kycR0KZRWS7xgeCbJ8R6W77J2k0jh1vV1wfAsLh6RyUtvHXu7VFReXc/dS7ZZQumFsEgmIvJeESkRkd0isugsx31QRFREzpjR0vgvlKY16Rfr4T2jslmxq8IWPooQbl5fpcfOXJTNJhPtnZBPJiLiAX4NXAkUATeJyBkzvolIf+DzwNrgRhi5vnLFaDp2/E2M87g2rcmccTkcqW1gx6FaV85vAusrV4ymY8/yYF1fR2o672Zuk4n2XMgnE2AGsFtVS1W1CXgEuKaT474P/BCwwQgBkp+VjALpSXFvr5bo5pois8bmIALLdx115fwmsLL7J6Dq7SEY7OsrlErdkSIcugbnAb5rcJYBM30PEJEpwFBVfVpEFnb1RiJyG3AbQH6+dTE9l8c3lpEY5+HVr80mJYhdgbsyIKUfk4em8+KuCr4wd7Tb4Zhe+uWKt8hN7cdLC2eRENf9GYEDYeH8Mdy9ZNu7OgC4WeqOBOFQMjkrEYkBfgJ8+VzHqurvVXWaqk7Lzs7u++DCWENzK09tOcSVEwaGRCJpN2dcLtvKazhaawXQcLZubxVr91Zx26Ujgp5I4J3JRAc5SyYkxdtKnr0VDsmkHBjq8/MQZ1u7/sB4YJWI7AMuAJ60RvjeWbbjCCcbW7h+auAnceyNuc5o+BetV1dY+9XK3WQlx3PTjKHnPriPLJicx+q753CDc43PLcp1LZZIEA7JZD0wSkQKRSQe+DDwZPtOVa1R1QGqWqCqBcAa4GpVtQXee+HxjWXkpSdyQWHg5kcKhNG5KQzJSORFazcJW1vLqnn5zUr+5z2FJMW7X+r98Iyh1DV5S+Km50I+mahqC3AHsAzYBTyqqjtE5B4Rudrd6CLToep6Xt19jA9OHUJMTN9P5OgPEWHuuFxe3X2M+iYbDR+OfrViN6kJsXysF+uUBNKU/AxG5aTwz/UHz32w6VLIJxMAVX1GVUer6ghV/YGz7duq+mQnx15upZLeeWJzOapw/ZTQquJqN2dcDo0tbby220bDh5viI7U8v/Mon7y4kP4JcW6HA3hvUG6cPpQtB6vZddi6nfdUWCQTEzyqyuMby5hZmEl+VpLb4XRqZmEWKf1iebHYqrrCza9X7iE53sMnLy5wO5R3uW7KEOI9MfzLSic9ZsnEvMumAyfYe+x0yDW8+4qPjeHS0QN4cVeFrcAYRvZUnuKprYf42IUFpCf1ftGrQMpMjmf++IEs2VRmk4n2UK+TiYjEiMjlIjLHGa1uwtjjG8tIivfwvgmD3A7lrOaOy6XiZCPbD9W4HYrppgdW7aFfbAy3vqfQ7VA69eHpQ6ltaOG57UfcDiUs9SiZiEiSiFwnIn8FKoAVwAtApYg8LCLXi0hKIAM1fa++qZWnthzmyvGDgrpeSU/MGpNDjHhXfjSh72BVHU9sLuemGfkMSOnndjidunB4FvmZSfxz3QG3QwlL3U4mIpIrIreKyFPAMeBxvGM6/gxcClwMPAhMAx7Fm1ieEZHbRGRg4EM3gRaqY0s6k5Ecz9RhGTaLcJj47Ut78Ihw26XD3Q6lSzEx3ob4tXurKK085XY4Ycefkskh4HdANt55sMY7PawWquqrqrpaVb+qqmPxTsj4PSAdeAAoE5HQ6LphuvT4xjKGZCQyszDT7VC6Zc64XHYerrXJ+ULckZoGHttQxvXThjAoLbTnvrph6hA8McK/NlhDvL/8SSa3A0NUdaaq3quqO7s6UFWLVfU+Vb0I79xanwWsVSuElVfX89qeY3xwSuiNLenK3HE5ALxoa8OHtN+/XEqrKp+9bITboZxTTmoCs8fm8O+NZTS1tLkdTljpdjJR1d+p6mF/T6CqR1T1D6pq/zMh7IlNZajCB0N0bElnRmSnMCwryUbDh7Bjpxr5x7r9LJiUx9DM0Oxq3tFNM4Zy7FQTK6zruV+sa7AJi7ElnRER5ozN5fU9x2097xD1x1f30tjSxu2zQr9U0u7SUdkMTE3gn+usqssfvUomIjJaRK4VkU87De3XisioQAVngmPj/hPsO17HDdPcm3Svp+aOy6GppY1X3rLR8KGmuq6Jv76+j6smDGJEdvh07oz1xPChaUN4+a1Kyq09rtv8TiYiMk5Efi4i5XjnynocbyP7b53nxSJySER+JiLjAhuu6QvtY0uuHB9+ne6mF2bSPyHWqrpC0F9e38fpplY+N2uk26H4rf3G6lEbEd9t/nQNHiEijwPbgf8BtuLtsfVx4H3AVc7ze4AtwK3AdhF5TERCtz9glKtvauWprYd534TQH1vSmThPDJeNzmZFcaWNhg8hpxpb+PNr+5hXlMu4Qaluh+O3oZlJXDJyAI9tOEirXVfd4s+3x05gG3ALsERVT5/tYBFJBq7Huy77TiChhzGaPrRsxxFOhcnYkq7MHZfLU1sPs6Wsmsn5GW6HY4CH1+ynpr6ZO8KwVNLuphn53P73Tbz8ViWzxuS4HU7I86ea6wZnlcK/nSuRAKjqaVV9SFWnADf2PETTlx7beJChmYnMKAiPsSWduXxMNp4YsQWzQkR9UysPvlLKpaOzOX9outvh9NjccblkJcfziI2I7xZ/ugafMd27H6/9T09fa/pO2Yk6Xt9zPKzGlnQmPal9NLy1m4SCR9Yf4NipJu6cHb6lEvBOKPrBqUN4cVcFFSdtmehz6UkDfIqIXCAi80XkIpsqJXw9sak87MaWdGXuuByKj5yk7ESd26FEtcaWVn73UikzCzOZHsal3XY3Th9KS5vy743l5z44yvnTAN9PRB4AKoHXgGeAV4ByESkXkb+KyFUiEr63uFFEVXl8UxkXDM8Mm8FkZzPHWRt+hY2Gd9W/N5ZzpLaBO2dHxgiBEdkpzCjI5F/rD6BqDfFn40/J5H7g08Aq4OvAV5xtgrdx/Wa8a7NvF5E5gQ3TBNqG/SfYf7yO66eG39iSzozITqFwQDIv2CzCrmlubeM3q3YzaWg6F4/McjucgPnwjKHsO17HmtIqt0MJaf4kkxuBP6nqlar6Q1X9KfAjZ98NwHBgEZAILBORLwY2VBNIj28I37ElXZkzNoe1pVWcarTR8G548o1DlJ2o587ZI4mkCor3TRhE/4RYHllvDfFn408ySQRWd7VTVfep6mJgDPBT4H4RmdfL+EwfqGtq4elt4Tu2pCtzxuXS1NrGK29Wuh1KVFm6uZyL7nuRLz+2hdgYoba+2e2QAiohzsO1k/N4dvsRquua3A4nZPmTTDYAs851kKo2q+pC4D/AN3oamOk77WNLbgjjsSWdmVaQQVpinK1xEkRLN5dz95JtHKr29nZqaVO+/sR2lm6OrAbrD0/Pp6mljSci7HMFkj/J5D7gI35UXz0DTPU/JNPXHt9YRn5mUkT0tvEV54nh8jHZrCypsFHLQbJ4WQn1HdZMr29uZfGyEpci6htFg1M5f0gaj6w7aA3xXfBnnMkyYCHe6qt1InIz0P8sL7kMsH6aISZSxpZ0Zc64XKpON/HGwRNuhxIVulqYLBIXLLtxej4lR0/yxsFqt0MJSX6NM1HVHwPzgWTgr8CbgAK3i8jXROR2EVkkIiuBjwB/D3TApneWOGNLrpuS53YofeKy0dnExohVdQXJ4PTOV07sans4u3rSYJLiPTxiU9N3yu9Bi6q6HBiPd2LHfwAHgeuAe4FfAf8HTAd+CXwtYJGaXmtft+TC4VkRMbakM2mJcUwvyLRZhINk4fwxxHneXcJNjPOwcP4YlyLqOyn9YvnAxMH8d+sh6zHYiR6tZ6Jez6rqLapagHdd+JnAHLyJJFtVv6CqkdWtI8yt33eCA1V1YT2pY3fMGZfDm0dPcbDKaln72oLJeUwakk6MeAec5aUncu91E1gwOTJLvjfOGEpdUyv/3XLI7VBCTkBWWlTV46q6XlVXqupGVY28CtMI8PjGgyTHe7hyQuSMLenMXGc0vM3VFRwn6puZNSaHvfddxWuLZkdsIgGYPDSdMbn9bfLHTvgznUqPR7WLyNyevtb03tLN5Vx474s8uqEMBZ7fEdlfsgUDksnpH899zxZTuOhpLr5vRcR1VQ0VJxua2VN5iolD0t0OJShEhA/PGMqWshp2Hqp1O5yQ4k/J5DkRWSEi7xcRz7kOFpE4Zxnfl/B2EzYuaB8HcLjGOw6grqmVu5dsi+gv16Wby6k63UxjSxsKlFfXR/xndsu28hpU4fyhaW6HEjTXTs7DI3DDb1+3mxUf/iSTyUAL3vm3DonI30Xk805yuUhELhaRD4jIl0TkUeAI3mV864BJAY/cdEu0jAPwtXhZCS0dxplE+md2y9ayGoCoKZkArCqpBITTTa12s+Kj23NpqOp24AoRuRC4HbgGuAlv12BfAtQCS4AHVHV9gGI1PRBN4wDaReNndsuWg9XkZyaRmRzvdihBs3hZCa3a+c1KJLcXnYvfEzOp6mpgtVPVNRUowtubS/FOT78d2KyqbYEM1PTM4PREyjv5Eo3EcQDtovEzu2VrWQ2T89PdDiOo7Galcz3uzaWqraq6TlX/oqqLVfV+Z5nejZZIQsfC+WPwxETHOIB2C+ePITHu3c16kf6Z3VB5spHy6nomhfHSvD0RTQM1/eFXMhGRfJ/HEBE523QqJgQsmJzHoNR+xHtiomIcAHg/873XTWBAirfqJTMpLuI/sxu2llUD0dVeAnaz0hV/q7n20aGNREQqgWXAT1R1S4DiMgFSU99MeU0Dd80exRfnjXY7nKBZMDmP900YxMTvLePqSXmWSPrAlrIaYgTG56W6HUpQtV9LP3hmF5UnG8lIiuM7Hzgv6q8xf6u5XvZ5vAJsxTvZ48eA9SLy+cCGZ3prw74qVGHm8MiaIbg74mNjmDosg7V7bYW8vrDlYDWjc/uTFB85a+J014LJeaxeNJukeA9Xnz846hMJ+D/R4+WqOst5XK6qk4E04ENAGfATZzbhgBKR94pIiYjsFpFFnez/jIhsE5E3RORVESkKdAzhau3eKuI9MUzJz3A7FFfMLMyi+EgtNXU2s08gqSpby6qZOCR6xpd0FOuJYXJ+OhsP2AzVEIDpVFS1RVUfxzsn1x7gpyIS1+vIHE6vsV8DV+LtOXZTJ8niH6o6QVUn4V1K+CeBOn+4W1N6nElD00mIO+c404g0szATVVi3z0ongXSwqp4Tdc2cH2WN7x1Nzc9g1+GTnLaJHwMzNxd45+cCvgNkAtMC9b7ADGC3qpaqahPwCN4xLr7n9p3XIJkzx75EpZMNzWwvr4nKKq525w9NJz42hrWlx90OJaJscRrfz4+yxveOpgzLoLVN2WJrnPjXAC8ibXTvi/pVEQHvBMO9rVDNwzvNfbsyvDMUd4ztc8CXgHhgdmdvJCK3AbcB5Ofn9zKs0Ldh/wnaFC4YnuV2KK5JiPMwaWi6tZsE2JaD1cTHxjBmYHR36JzsVB9v3H+Ci0YOcDkad/n7Rf8MZ08mWXi/6NfjHcAYNKr6a+DXIvIR4JvAJzo55vfA7wGmTZsW8aWXNaXHifNI1LaXtLugMJNfrdxNbUMzqQkBq4GNalvLajhvcCpxnoBVboSltMQ4RuemWLsJfiYTVX3/2fY768PPBG5T1a29CcxHOTDU5+chzrauPAI8EKBzh7W1pVVMHJJOYnx0tpe0mzk8i1+s2M3GfSeYNTbH7XDCXktrG9vKa7hx+tBzHxwFpg7L4Omth2lr04hcCru7AnZbISLnAd8CSgKYSMBbyhklIoUiEg98GO9kk77nHuXz41XAWwE8f1g63djCtvIaLoji9pJ2U/IziPMIa/Zau0kg7K48RX1za1TNFHw2U4dlUtvQwu7KU26H4qpeJxNnJPzdwGogEadNIlBUtQW4A+/AyF3Ao6q6Q0TuEZGrncPuEJEdIvIG3naTM6q4os2G/SdobVNmFkZve0m7xHgPE4eks87aTQJi60HvTMHR3vjebuqwd9pNopm/DfClvj8CqUC687wGuElVXw1YdA5VfYYOa6Ko6rd9nttgyQ7Wlh7HEyNvX+jRbmZhJr9/uZS6ppaoHGQXSG+UVdM/IZaCrGS3QwkJBVneWZM37j/BTTMiv2NPV/wtmRT4PPKd128A7gHGqurTAYzN9MLavVVMHJJGcj/74gRvu0lLm0b93WMgtA9WjOb2AV8i3k4um6L82vJ3BHyMz8OjqhmqOlNVv6uqkb0WbBipa2phy8Fqq+LyMXVYBp4YYW2pVXX1RkNzK8WHT1oVVwdTh2VQeuw0Vaeb3A7FNdHdry9CbdpfTUubWuO7j5R+sYwfnMpaa4TvlZ2Ha2lp06ibKfhcphV4q5OjuXRiySQCrXHaS6YVWDLxNXN4FlsO1tDQYRlj033tI72jbQ2Tc5mQl0acR9gQxcmk2xXqIrIiAOdTVZ0TgPcxZ7F273HGD04lxdpL3qW9EX7TgRNcNCK6Ryv31NayGnL692NgWoLboYSUhDgP5w1Os5KJH8dKLx9WEupj9U2tbDlYE9VTqHRlWkEmIli7SS9sOVgd9ZM7dmXqsAy2lFXT1BKdC812+9ZVVS/vwzhMgGw+cIKm1raontyxK2mJcRQNSrXxJj1UU99M6bHTXDfF1u7ozNRhGfzx1b3sPFwbldWAVlKIMGv2VhEjWHtJF2YWZrHpwAkaW6zdxF/bypzBilH4Rdkd0T54MeDJRETsr9RFa0uPc97gNJvQsAszh2fS2NLGVueL0XRf+7TzE/PSXY0jVOWmJjAkI5GN+6Oz5NsXJRMbyeSShuZWNh+sZmahlUq6MsMpsdn6Jv7bWlZN4YBk0pLsRqUrU4dlsHH/CVQjflLyM5wzmYjI6yIyyY/3jL7fYoh446C38c8a37uWkRzP2IH9bX2THthysCaql+ntjqnDMjha20h5db3boQRdd0omFwDLu0ooInJNZ9tN8K0trUIEplvJ5KxmFmaycf8Jmlujs9dNTxytbeBIbYMNVjyHKfnR227SnWRyGu/6IMtF5PxO9v8tsCGZnlpTepxxA1NJS7RqiLOZUZhFXVMr28qt3aS73hmsaCWTsxk7sD9J8Z6oHG/SnWSiqvot4LfAi50kFGsjCQGNLa1sOnDCqri6YUZhe7uJVXV119ayGjwxQtEgSyZnE+uJYXJ+elSOhO92A7yqfhP4Hd4SykTfXQGPyvhta1kNjS02vqQ7svv3Y0R2ss3T5YctZdWMye0f9at2dsfU/Ax2Ha7ldGOL26EElb+zBn8D7xrqL4rIhL4JyfTEmj3eL8YZNr6kW2YOz2LDPu8CYubsVNUZ+W6lku6YMiyDNn2najBadCeZvKsay0kofwBWdCihGBet3VvF2IH9yUiOdzuUsDCzMJNTjS3sPFTrdighb9/xOmobWmza+W6aHKWN8N1JJos6blDVr+NNKC8C/QIdlPFPU0sbG/dbe4k/2n9XVtV1blvbBytaMumWtMQ4RuemsPGAJZN3UdVfd7H968CDgHUdctm28mrqm1tt/RI/5KYmUJCVxBprhD+nNw5WkxAXw+jcFLdDCRtTh2Wyaf8J2qKoGrVXI+BV9W6gs+7CJojavxBn2MqKfplZmMX6fVVR9QffE1vLahg/OI1Yj03l111Th2VQ29DC7spTbocSNL2+OlR1WyACMT23dm8Vo3NTyLT2Er/MHJ5JTX0zxUdOuh1KyGpubWPHoRqb3NFP0Tjpo91qhLnm1jY27Kuy9pIeeHu8ibWbdOnNoydpaG6zaVT8VJCVRGZyvCWTsxGRS0XEOpuHiO3lNdQ1tTLTqrj8NiQjibz0RBu8eBbtsytbTy7/iAhT8jOiaiR8T0omK4Fu3aaIyPUi8jMR+ZSIxHbY93QPzm06aJ+wcIbNx9UjM4dnsm5fVVTO8todWw5Wk5YYx7CsJLdDCTtTh2VQeuw0Vaeb3A4lKHqSTLo1fYqI3AH8CkgCFgKviYjvN957enBu08Ga0uOMzEkhu7/10O6JCwqzqDrdxO6K6Gko9ceWMu9MwSI2a5K/phVEV7tJX7aZ3AHMV9XbgAnAG3gHOrYnFLs6e6mltY0N+07Y+iW90D79zBqbkv4M9U2tvHn0ZFQuQRsIE/LSiPOIJZMAGKyqWwBUtUVVPw2sAFaKSBY2p1ev7Txcy6nGFmt874X8zCQGpibYYlmd2HGohtY2tcGKPZQQ5+G8wWlR027Sl8mkUkQKfTeo6pfwtrmsBGI7fZXptjXOF6BN7thzIsLM4Zms3WvtJh294cwtdb715OqxqcMy2FLmXbQu0vUkmXT3L24FcMsZL1b9ArAKSOjBuY2PtaVVDB+QTE5/+1X2xszCLCpPNrL32Gm3QwkpW8tqGJSWQE6qXV89NXVYBo0tbew8HPlzwPVZAzzwOeCHne1Q1buAgh6c2zha25R1e6uYaVVcvdZesrOlfN9tS1m1dQnupWgavNiTZPIvoPFcB6lqk6rWnWX/gR6c2zh2Ha7lZGOLzccVAMMHJDMgpZ+1m/iormti//E6Jtq0872Sm5rAkIxENu6P/BsVv5OJqt6kqlYf4LK320tssGKviQgzC63dxJcNVgycqcMy2Lj/RMRfWzadSphaU1pFQVYSA9OsPjsQZg7P5HBNA2Un6t0OJSS0L+w0wRrfe23qsAyO1jZSXh3Z15YlkzDU1qas31dlpZIAav9drrGqLsA7WHF4djKpCbbCRG9NiZLFsiyZhKHiIyepqW+2LsEBNConhYykOGuEx1mmt6yaSVbFFRBjB/YnKd4T8eNNepVMRCReREoDFYzpnnfGl1jJJFBiYoQZhZk2gzBwpLaBypONNlNwgMR6Ypicn84GSyZnJVgX36Bbu/c4QzMTyUtPdDuUiDKzMIuDVfUcivC67XNpby+xNUwCZ2p+BrsO13K6scXtUPrMOZOJiJR29QBKCMK0KCLyXhEpEZHdInLGmvQi8iUR2SkiW0XkRREZ1tcxuaXNGV9ygbWXBNw7402iu3SypayG2Bhh3KBUt0OJGFOGZdCm7yTqSNSdKU2yga8DnY0LiQceCWhEHThrp/wamAeUAetF5ElV3elz2GZgmqrWichngR8BN/ZlXG55s+IkJ+qarYqrD4wdmEpqQixrS6u4dvIQt8NxzZaD1YwblEpCnC1bFCiTfRrhLxo5wOVo+kZ3kskbwFFV/U/HHSLSj76f/XcGsFtVS51zPgJcA7ydTFR1pc/xa4Cb+zgm17Qv5GQzBQee5+12k+hthG9rU7aV1XD1pMFuhxJR0hLjGJ2bwsYDkdtu0p02k18AXZX7m4FPBi6cTuUBB31+LnO2deV/gGc72yEit4nIBhHZUFlZGcAQg2ft3uPkpScyNNMWK+oLMwoz2XvsNBW1DW6H4orSY6c52dhi7SV9YOqwTDbtP0FbW2QOXjxnMlHVx1T1xS72tanqQx02u7ZOiYjcDEwDFne2X1V/r6rTVHVadnZ2cIMLAFVlbWmVdQnuQ+3jTaK1dLK1rBqwke99YeqwDGobWthdGZkLsQV8nImqBvo9y4GhPj8Pcba9i4jMBb4BXK2q55w7LBztrjjF8dNN1vjeh84bnEpKv9iobYTfWlZDUryHkTkpbocScSJ90sdwGLS4HhglIoUiEg98GHjS9wARmQz8Dm8iqXAhxqBoXw3QFsPqO7GeGKYOy3i7bSravHGwmvF5aXhibCHUQCvISiIzOT5ik0m3F6gSkRUBOJ+q6hw/X9DirCe/DPAAf1LVHSJyD7BBVZ/EW62VAjzmrFV9QFWvDkC8IWPp5nL+72lvn4Ob/rCahfPHsmDy2ZqOTE/1T/DwVsUpChc9zeD0RBbOHxMVv+smZ92NT1wYsT3rXSUiTMnPiNiR8P6sdhhD78eU9Oh2R1WfAZ7psO3bPs/n9jKukLZ0czl3L9lKfbN3tbby6gbuXrINICq+5IJp6eZyXtjpLdwqUF5dHzW/65IjJ2lqabPG9z40dVgGy3cdpep0E5nJ8W6HE1DdTiaqenkfxmHOYvGykrcTSbv65lYWLyuJ+C+4YFu8rITGluj7XS/dXM53ntwBwPef2klLq0b053XLtIJ32k3mFeW6HE1ghUObSdTranqPaJ/2oy9E4+/aW/LdRk19MwBHaxu5e8k2lm4+o5+L6aUJeWnEeSQi200smYSBwV3MwdXVdtNz0fi79pZ8W9+1rb00ZgIrIc7DeYPTIrLdpLezBo8WkWtF5NPOgMBrRWRUoIIzXnfOHnHGtsQ4Dwvnj3Ehmsi2cP4YEjtMIxLpv+toLI25aeqwDLaUVdPUoTo13PnTAA+AiIwDPgNcDwxs3+z8q84xR4FHgd+p6q4AxBnV0hK9DXUDUuI5fqopqnoYBVv773TxshLKq+uJ8wj3Xjchon/Xg9MTO10FMJJLY25qbWujsaWNMd98NqL+lv3pGjwC+CFwLVAPvIJ3bMcevNOtCJAJjAQuAG4F7hSRJcDX2ufWMv5bUVxBWmIca+6eQ6zHaib72oLJeSyYnMcvX3yLnyx/k4tGRva4noXzx/DlR7fQ6rNGeaSXxtyydHM5j6zzzg4Vab0F/flm2ol37ZJbgFxVvVJV71HVv6vqc6r6rPP8e6p6JZALfAoYgc+kjMY/bW3KypJKLh2dbYkkyOYW5aIKK4sjdhwsAB84fzD9YoXEOA8C5KUnRnxpzC2Ll5XQ0EVvwXDnTzXXDc4AwW5R1dPAQ8BDInKN35EZALYfquHYqUZmjw2/ucTC3diB/clLT+SFnRXcOD3f7XD6zBsHq6lrbuMXN03m6vNttuC+FMntU92+1fUnkXTy2jOmrzfds6K4AhG4bHSO26FEHRFh7rgcXt1dSUOH3k6RZGVxBZ4Y4bJRdsPS1yK5t6Df9SYikiIiF4jIfBG5SEQGnvtVpqdWFlcweWh6xI2WDRdzi3JpaG7jtd3H3A6lz6wormBqfgZpSXFuhxLxIrm3YLeTiYj0E5EHgErgNbzTm7wClItIuYj8VUSuEmdyLNN7lScb2VJWw+yxVipxy8zCLFL6xbJ811G3Q+kTR2oa2Hm4lll2jQXFgsl53HvdBPKckognRvi/a8dHRPuUPyWT+4FPA6vwLuP7FWebAAl4Vzd8EtguIn5N5mg6t6rE2/Brf+juiY+N4bIx2SzfVRGRixqtdK4xu2EJngWT83ht0WzuvW4CrW3K2EGpbocUEP4kkxvxzth7par+UFV/inetdYAbgOHAIiARWCYiXwxsqNFnZUkFA1MTKIqQiy1czR2XQ+XJRraW17gdSsCtKK4gLz2R0bm2fkmwzRmXgwi8sDMySr3+JJNEYHVXO1V1n6ouBsYAPwXuF5F5vYwvajW3tvHKm8eYNTYbqzl016wxOXhihBcjrKqrsaWV13bbNeaWnP4JTBqaHpXJZAMw61wHqWqzqi4E/oN35UPTA+v3VXGysYVZY6z6wW3pSfFMG5YRMX/07daWVlHX1GpVXC6aV5TLtvIaDtdEUddg4D7gI35UXz0DTPU/JAPeXlzxnhguHjnA7VAM3j/64iMnOVhV53YoAbOiuIJ+sTFcONyuMbdcUeTtDLs8Am5U/BlnsgxYiLf6ap2I3Az0P8tLLgMi5y8vyFYUVzBzeCbJ/fyePs30gTnjvGtPREpVl6qysqSCi0ZkkRjvOfcLTJ8YmZPC8AHJPB9NyQRAVX8MzAeSgb8Cb+KdYuZ2EfmaiNwuIotEZCXwEeDvgQ44Ghw4XseeytNW/RBCCgckMyI7meW7ImNqldJjp9l/vM6usRAwryiXNaXHqW1odjuUXvF70KKqLgfGA1cB/wAOAtcB9wK/Av4PmA78EvhawCKNIiuKvXcp9oceWuYW5bJ2b/j/0cM7841Zt3P3zSvKpblVWVVS6XYovdKjmQPV61lVvUVVC4BsYCYwB28iyVbVL6hq+P/VuWBFSSXDs5MZlpXsdijGx7xx3j/6l98M7z968Fajjs5NYUhGktuhRL3J+RlkJceHfQePgExDq6rHVXW9qq5U1Y2qGv5dE1xyurGFNXuOM9t6cYWcyfkZZCbHh31j6cmGZtbtrbJSSYjwxAhzx+WyqrgirBfM8mc6lR6PaheRuT19bbR5bfcxmlrbrIorBHlihFljclhZUklLa/j+0b/61jFa2tRuWELIvKJcTja2sKb0uNuh9Jg/JZPnRGSFiLxfRM7Z/UNE4pxlfF/C203YdMPKkgpS+sUyrSDT7VBMJ+YV5VBT38yGMF7De0VxBakJsUwdluF2KMZxyagBJMZ5wrqqy59kMhlowTv/1iER+buIfN5JLheJyMUi8gER+ZKIPAocAR7H2z14UsAjj0CqysriSt4zagDxsbYQVih6z6hs4j0xYVvVZYuthaaEOA/vGTWA5buOohqec8B1exCDqm4HrhCRC4HbgWuAm3DWffchQC2wBHhAVdcHKNaIt/NwLUdqG6wuO4Ql94vlwhFZvLDrKN+4alzYTUPyzmJrdo2FmnlFuTy/8yjby2uZMCTN7XD85veIOFVdDax2qrqmAkV4e3Mp3unptwObVTV8K5Vd0t5d8/IxtkhRKJtblMu3lm5nT+UpRuacbdxu6HlnsTW7xkLNnHG5xAg8v/NIWCaTHpdzVbVVVdep6l9UdbGq3q+qDzm9uSyR9MCK4gomDkkjp3+C26GYs5g7zntX/8LO8BvAuLK4gklD08lK6ed2KKaDzOR4phVkhm27iV/JRETyfR5DRCS8bstCWNXpJjYfrLaJHcPAoLRExuelht3UKm8vtmbXWMi6IozngPO3ZLIP2Os89gPVInJERB4SkfMDHVw0eenNClRt1Hu4mDM2l40HTnD8VKPboXTbSltsLeTNK/LOAReOc3X5m0xe9nm8AmzFO9njx4D1IvL5wIYXPVYUVzIgpR8T8sKvrjQazSvKRdVbNRkuVhZXkJvaj/MG22JroWpYVjKjc1N4YecRt0Pxm78TPV6uqrOcx+WqOhlIAz4ElAE/cWYTNn5oaW3jpZIKLh+TTUxMePUOilbnDU5lYGpC2KwN39TSxitvHWPWmJyw64EWba4oGsi6vVWcON3kdih+6XVHc1VtUdXH8c7JtQf4qYjE9TqyKLLpQDW1DS1WxRVGRIS5RTm88tYxGppb3Q7nnDbsq+JUY4tVcYWBeUW5tIVZqRcCNDcXeOfnAr4DZALTAvW+0WBFcQWxMcIlo2yRonAyZ1wudU2trA6DKTBWOIutXWKLrYW8CXlp5Kb2C7teXf725moTkdauHsDDzqGvOttaAh9y5FlZXMGMwkxSE6xAF04uHJ5FUrwnLEbDryixxdbCRYwz8ePLb1WGRam3nb8lk2fO8VjrHLfe+fnZwIQZucqr6yk5etKquMJQQpyHS0dl8+KuipCeAmP/8dOUVp62budhZF6Rt9T7+p5jbofSbX7dpqjq+8+231kffiZwm6pu7U1g0WKFLVIU1uYW5fLcjiPsOFTL+BDtidd+jdkNS/i4cEQWKf1ieX7HUWaPzXU7nG4JWJuJiJwHfAsoCXQiEZH3ikiJiOwWkUWd7L9URDaJSIuIXB/Ic/e1lcUVDMtKYvgAWwgrHM0ak02MENL12yuKKxg+IJkCu8bCRr9YD5eNyWb5rgra2kK31Our18nEGQl/N7AaSARu63VU735/D/Br4Eq884DdJCJFHQ47ANyCdxnhsNHQ7C3GWnfN8JWV0o8p+Rkh20X4dGMLa0urrFQShq4oyuXYqUY2H6x2O5Ru8bcBvtTnsVdEjuMdCf8DoBW4XlVfDXCMM4Ddqlqqqk3AI3hnLH6bqu5zSkNhNSfY6j3HaWi2hbDC3dyiXHYcquVQdegtMGqLrYWvy8fkEBsjIV3q9eVvyaTA55HvvH4DcA8wVlWfDmBs7fKAgz4/lznb/CYit4nIBhHZUFnp/jreK4orSIr3MHO4LYQVzuaO89Zph+JcXbbYWvhKS4xj5vDMsBkN7+8I+Bifh0dVM1R1pqp+V1VD7y+pA1X9vapOU9Vp2dnuTsGtqqworuDikQPoF3vOhStNCBuRnUzhgGSW7wqtQWa22Fr4u6JoIHsqT7On8pTboZxTOFxh5cBQn5+HONvC2lsVpyivrrfqhwggIswZm8PqPcc51Rg6Q6tssbXwN9eZ+DEcqrq6nUxE5HYRGeLvCURkkIh8WkR6mrjWA6NEpFBE4oEP4106OKy93SXY+v5HhLlFuTS1tvHKm+5Xn7azxdbCX156IucNTo2sZAL8HNjvtDl8U0QmdnWgiIwXkW+IyFq8bRy/AnpUl6OqLcAdwDJgF/Coqu4QkXtE5GrnfNNFpAy4AfidiOzoybmCaUVxBUWDUhmYZgthRYJpwzJIS4zjhRBqN7HF1iLDvKJcNh04QeXJ0F7uwJ9kkgt8Cm833EXAZqdH189EZLaIXCYiPxGRPcAW4OvAYeBWYJCqNvc0SFV9RlVHq+oIVf2Bs+3bqvqk83y9qg5R1WRVzVLV83p6rmCoqWtm4/4TVsUVQWI9Mcwem8PK4gpaQ2BcgC22FjneWe4gdG5UOtPtZKKqVc6yvNcBA4AFwAq81U7Lnec3Ay8B1wEDVHWBqv5ZVcNnToAgeOmtSlrb1OqyI8yccTmcqGtm04ETbodii61FkKJBqeSlJ/L8jtBOJj2a9U1VG4D/Av8V72i7C/Emptdt/fdzW1lcQWZyPJOGprsdigmgS0dnE+cRlu88ynSXu+LaYmuRQ0SYV5TLP9cdoK6phaT40JysMxDrmaiqvq6qr1oiObfWNmVVSQWXjc7GYwthRZTUhDguGJ7leruJLbYWea4oyqWxpY2X3wzdSp5w6BocUd44WM2Jumar4opQ2f37UVp5msJFT3PxfStYujn4vdhtsbXIM70wk9SE2JDu1WXJJMhWFlfgiREuG2XdNSPN0s3lPL31MACKd3mBu5dsC3pCscXWIk+c08FjRfFRWlpDswLIkkmQLN1czsX3reBXK3fjEWFlSWiNlja9t3hZCY0t7/5Dr29uZfGykqDGsbK4gukFtthapLnivIGcqGtmw373O3h0xpJJECzdXM7dS7ZR7kwE2NTa5sodq+lbXU30GMwJIG2xtch16ehs4j0xIVvVZckkCBYvK6G+w/Kbbtyxmr41OD3Rr+2BtnRzOe/7+SsA/OGVUrtZiTAp/WK5aGQWL+w8GpIre1oyCYJQuGM1fW/h/DEkxr17ogdPjLBw/pg+P3d76bem3js2uOJko5V+I1B2/34cqKpj+N3PuNbBoyuWTIJgUBdTpgTrjtUEx4LJedx73QTy0hMRICneQ2ubMjq3f5+f20q/kW/p5nL++8YhwN0OHl2xZBIEhQOSztiWGOcJyh2rCa4Fk/N4bdFs9t53FasXzSErOZ5v/Wd7ny+9aqXfyLd4WQkNIdDBoyuWTPrY8p1HeW1PFbPGZL99x5qXnsi9101gweQerfFlwkRaUhyLrhzLxv0neHxjWZ+eKyslvtPtVvqNHKF+wxCa4/IjRMXJBr76760UDUrltx+baotgRaEPThnCoxsOcu+zu5hXlEtGcudf+r1RdbqJppY2BG/1Rzsr/UaWwemJb/cI7bg9FFjJpI+0tSlfeWwrpxtb+MVNkyyRRKmYGOH7C8ZT29DCj/qgOkJV+cpjW2hobuNLV4y20m8E66yDB8CHpvu9zFSfsJJJH3lo9T5efrOS7y8Yz8icvm+ANaFr7MBUPnlRAX98bS8fmjaEyfkZAXvvP766lxXFFXzv6vP4xEUF3Dl7VMDe24SW9huDxctKOFRdT25qAk2trfz19f0smJTHsKxkV+OTUOyvHAzTpk3TDRs29Ml7Fx+p5epfvcZ7Rg7gwU9MwzuxsolmpxpbmPPjVWT378d/PndJQCb5fONgNdc/8DpzxuXw25un2nUWhfZUnuKDD7xORlI8//7sRWT2QTVqRyKyUVWnddxu1VwB1tDcyuf/+QapCXH88PqJ9gduAO+As2+9v4jt5bX8fe3+Xr9fTX0zd/xjE7mpCfzog+fbdRalRmSn8ODHp1FeXc+tD62noUP38GCyZBJgP3yumJKjJ7n/hokMSOnndjgmhFw1YRCXjBzA4mUlvVqCVVX52uNbOVLTwC8/Mpm0JJuDK5pNK8jk5zdOYvPBaj7/yGbXVvq0ZBJAq0oq+PNr+7jlogIut+VSTQciwj3XnEdjcxv3PrOrx+/ztzX7eW7HEb763jFMCWD7iwlfV04YxDevKmLZjqN8/6mdrky3YskkQI6fauQrj21lTG5/Fl051u1wTIganp3CbZcOZ8nmctaUHvf79dvLa/jfp3Yxa0w2t14yvA8iNOHqfy4p5FMXF/KX1/fxx1f3Bv38lkwCQFX52r+3UtvQzM9vmkRCJ933jGn3uVkjGZKRyLeWbqfZj7UpTjW2cMc/NpGZHM+PPzTJVlE0Z/jmVeO4cvxA/vfpXW+vrRMslkwC4O9rD7B8VwWL3juWsQNT3Q7HhLjEeA/f/cB5vFVxij918w5SVfn6km0cqKrjFzdNDkqvHRN+YmKEn944ianDMvjio2+wfl9V8M4dtDNFqN0VJ/nfp3dy6ehsbrmowO1wTJiYW5TL3HG5/Gz5W92aDuPRDQd5csshvjRvNDMKM4MQoQlXCXEeHvz4NIakJ3LrQxvYXXEqKOe1ZNILjS2t3PXPN0iKj+X+6ydatYPxy3c+UISifP+pnWc9ruTISb7z5A4uGTmAz14+MkjRmXCWkRzPXz45gziPcMuf11FxsqHPz2nJpBd+8vyb7Dxcyw8/OJGc1M6nmTemK0Mzk7hz9iie3X6EVV0s41zX5G0nSekXx09uPD8ggx1NdMjPSuJPt0zn+Kkm/ucvGzjd2NKn57Nk0kOv7z7G718p5aMz85lXlOt2OCZM3fqeQoYPSOY7T+7odMDZd5/cwe7KU/zsxknk9LcbFuOfiUPS+dVHJrPjUA13/GMTLX50+PCXzc3lh6Wby9+eF0cEBqTE882ritwOy4SxfrEe7rlmPDf/cS2/fWkPX5g7+u19T2wu49ENZdw5eySXjBrgYpQmnM0Zl8v3F4znG09sZ+L3nqe+qZXB6YksnD8moBOBWsmkm9qXRS2vrkeBNoXa+haW7TjidmgmzF0yagAfOH8wv1m1h/3HTwNQWnmKbzyxnRkFmXx+jk3eaHonOT6W2Bihrqm1z1ZptGTSTZ0ti9rY0hYyq5yZ8PbNq8Yhqlzx05cpXPQ0V/z0ZVDl5zdNItZjf6amdxYvK6GlwzQrgV6l0a7Sbgr1Vc5MeFu95zhteG9QFGhpU1raYG1p8MYJmMgVjO8vSybd1NVqZqGyypkJb4uXldDc+u47x6ZWK/mawAjG95clk27qbJUzWxbVBIqVfE1fCsb3l/Xm6qaOq5z1RW8IE71CfX1vE96C8f1lKy0aEwLaewv6dvJIjPPYOu4m5HS10qKVTIwJAVbyNeHOkokxIWLB5DxLHiZsWQO8McaYXguLZCIi7xWREhHZLSKLOtnfT0T+5exfKyIFLoRpjDFRK+STiYh4gF8DVwJFwE0i0nFCrP8BTqjqSOCnwA+DG6UxxkS3kE8mwAxgt6qWqmoT8AhwTYdjrgEecp4/DswREZur2xhjgiQckkkecNDn5zJnW6fHqGoLUANkdXwjEblNRDaIyIbKyso+CtcYY6JPVPXmUtXfA78HEJFKEdnfw7caABwLWGDuCPfPYPG7L9w/Q7jHD+58hmGdbQyHZFIODPX5eYizrbNjykQkFkgDjp/tTVU1u6cBiciGzgbthJNw/wwWv/vC/TOEe/wQWp8hHKq51gOjRKRQROKBDwNPdjjmSeATzvPrgRUarUP7jTHGBSFfMlHVFhG5A1gGeIA/qeoOEbkH2KCqTwJ/BP4mIruBKrwJxxhjTJCEfDIBUNVngGc6bPu2z/MG4IYghvT7IJ6rr4T7Z7D43RfunyHc44cQ+gxRO9GjMcaYwAmHNhNjjDEhzpKJMcaYXrNk4qdzzRMWykTkTyJSISLb3Y6lp0RkqIisFJGdIrJDRD7vdkz+EJEEEVknIluc+L/ndkw9ISIeEdksIk+5HUtPiMg+EdkmIm+ISNgtbCQi6SLyuIgUi8guEbnQ9ZiszaT7nHnC3gTm4R2Jvx64SVV3uhpYN4nIpcAp4K+qOt7teHpCRAYBg1R1k4j0BzYCC8Lo/0CAZFU9JSJxwKvA51V1jcuh+UVEvgRMA1JV9f1ux+MvEdkHTFPVsBy0KCIPAa+o6oPOkIkkVa12MyYrmfinO/OEhSxVfRlv1+mwpaqHVXWT8/wksIszp9cJWep1yvkxznmE1R2diAwBrgIedDuWaCQiacCleIdEoKpNbicSsGTir+7ME2aCxFlqYDKw1uVQ/OJUEb0BVAAvqGpYxQ/8DPgq0OZyHL2hwPMislFEbnM7GD8VApXAn52qxgdFJNntoCyZmLAkIinAv4EvqGqt2/H4Q1VbVXUS3qmBZohI2FQ5isj7gQpV3eh2LL10iapOwbu0xeecKuBwEQtMAR5Q1cnAacD19ltLJv7pzjxhpo85bQ3/Bv6uqkvcjqennKqJlcB7XQ7FHxcDVzttDo8As0XkYXdD8p+qljv/VgBP4K3CDhdlQJlPifZxvMnFVZZM/NOdecJMH3IasP8I7FLVn7gdj79EJFtE0p3niXg7cxS7GpQfVPVuVR2iqgV4r/8Vqnqzy2H5RUSSnc4bONVDVwBh08NRVY8AB0VkjLNpDuB6B5SwmE4lVHQ1T5jLYXWbiPwTuBwYICJlwHdU9Y/uRuW3i4GPAducdgeArztT7oSDQcBDTs/AGOBRVQ3L7rVhLBd4wlk/Lxb4h6o+525IfrsT+LtzU1sKfNLleKxrsDHGmN6zai5jjDG9ZsnEGGNMr1kyMcYY02uWTIwxxvSaJRNjjDG9ZsnEmBAiIvNFZJWInBKRShH5lYgkuB2XMediycSYECEiXwaeAw4DXwT+C3wO+LmbcRnTHTbOxJgQICJzgeeBr6rq/T7bnwNmAdnhNgeZiS5WMjHGZSISg7f0sRn4cYfdq4B4IGwmgzTRyaZTMcZ984Ei4BY9s6qgyfk3LbghGeMfSybGuO9GoBV4RUQGdNiX6/x7MrghGeMfazMxxmUish/IP8dheap6KBjxGNMTlkyMcZFTEqnEu6bGbzo55FGgUVUHBTUwY/xk1VzGuGu48+96VV3uu0NECoEM4B9Bj8oYP1lvLmPcleL821mbyPXOv/8KUizG9JglE2Pc1T52JNV3o7Po0WeBEuDpYAdljL8smRjjrp1AHd7uwb5+ABQAd6lqa7CDMsZf1mZijItUtU5EHgTuEpGHgZeAK4FrgYWq+ryrARrTTdabyxiXOVVaPwI+CiQBG4H/C8N1yU0Us2RijDGm16zNxBhjTK9ZMjHGGNNrlkyMMcb0miUTY4wxvWbJxBhjTK9ZMjHGGNNrlkyMMcb0miUTY4wxvWbJxBhjTK/9f7SHFRCeSkIHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "_ = plt.figure(1)\n",
    "_ = plt.plot(thetas, results, 'o-')\n",
    "_ = plt.xlabel(r'$\\theta$', fontsize=18)\n",
    "_ = plt.ylabel(r'$\\langle \\Psi(\\theta) | \\frac{1 - Z}{2} | \\Psi(\\theta) \\rangle$', fontsize=18)    \n",
    "_ = plt.title('Noisy Energy Landscape')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1326fbbd",
   "metadata": {},
   "source": [
    "Run optimization without noise mitigation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "97f1a06e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     fun: 0.021\n",
      "   maxcv: 0.0\n",
      " message: 'Optimization terminated successfully.'\n",
      "    nfev: 20\n",
      "  status: 1\n",
      " success: True\n",
      "       x: array([-0.08955312])\n"
     ]
    }
   ],
   "source": [
    "quil_prog = qc.compiler.quil_to_native_quil(program, protoquil=True)\n",
    "init_angle = [0.1]\n",
    "res = optimize.minimize(expectation, init_angle, args=(quil_prog), method='COBYLA')\n",
    "print(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69e97df8",
   "metadata": {},
   "source": [
    "Now run parameter scan with ZNE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6479ce2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "quil_prog = qc.compiler.quil_to_native_quil(program, protoquil=True)\n",
    "init_angle = [0.1]\n",
    "results_zne = []\n",
    "# here we use a linear fit for the zero noise extrapolation\n",
    "fac = mitiq.zne.inference.LinearFactory(scale_factors=[1.0, 3.0])\n",
    "for theta in thetas:\n",
    "    results_zne.append(\n",
    "        zne.execute_with_zne(quil_prog, lambda p: expectation(theta, p), factory=fac)\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2132113f",
   "metadata": {},
   "source": [
    "Plot the energy landscape with ZNE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "5aacde8e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEeCAYAAABPMvhnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABHOklEQVR4nO3deXzcdZ348dc799E0aZP0SpqeSUpp6d1CC6UVpCAqBZFDAVl10Z+g7KpVwFVXXVfWul67rouiLqeACAUELVcppS090vSgR3qnbdI2d5M0d/L+/THfYEiTJpOZzHeO9/PxmEeT70zm+55pMu/v53p/RFUxxhhjfBHldgDGGGNCnyUTY4wxPrNkYowxxmeWTIwxxvjMkokxxhifWTIxxhjjM0smJmBE5H9F5Nvnuf8BEXk4kDH1RkT+T0T+ze04QpGI/KuIPO52HCawLJkYn4nIURFpEZGMbscLRURFZDyAqn5RVX/g3LdERE50fbyq/ruqft5PMamITPbHc/Xw3HeKSLuI1He7jRmM8/mLfcibwWTJxPjLEeDWzm9EZDqQ5F44g26jqg7pdiv15wlEJMafz2fMYLJkYvzlMeCOLt9/Bni06wM6u45EJBn4KzCm61V99ytnEblDRIpFpFJEvu20gK507psvIhtFpEZETorIf4tInHPf285T7HCe+2bn+EdFZLvzMxtE5KIu55olIttEpE5EngYSBvpGOHF+XUR2isgZEXlaRBK63H++OI6KyDdFZCdwVkRiensfRGSUiDSISHqXn58tIuUiEutlzPeJyCHn9e8Rkeu73HeniLwjIj8RkWoROSIi13S5f4KIrHV+9jUgo8t9CSLyuBN7jYhsEZGRzn3DReQPIlLqPO8q5/gwEfmL8zqqna+zuzznWyLyIxHZLCK1IvKCiAzvcv/FzvtaIyI7RGSJN++FGRhLJsZf3gWGisgFIhIN3AL02KWiqmeBa4DS3q7qRWQq8D/Ap4HRQCqQ1eUh7cA/4/ngugS4AviS8/yLncfMcJ77aRGZBfwe+AKQDjwEvCgi8U4SWoUnIQ4H/gR8wpc3A7gJuBqYAFwE3Om8rl7j6PKztwLXAmlAXm/vg6qeAt5yztXpduApVW31Mt5DwGXO838PeFxERne5fwFQhOf9/jHwOxER574ngQLnvh/guZDo9BnnOcc6r/eLQKNz32N4Wq8XAiOAnznHo4A/AOOAHOfx/90t3juAz+J5T9qAXwKISBbwMvBveP4vvw78WUQyvXw/jJcsmRh/6mydfBjYC5T48Fw3Ai+p6juq2gJ8B3i/kJyqFqjqu6rapqpH8XwoX36e57sLeEhVN6lqu6o+AjQDFzu3WODnqtqqqs8CW/qI72Lnyrfzdqjb/b9U1VJVrQJeAmb2I46uP3tcVRv7eh+AR4DbAJwkfiue/wevqOqfnHg7VPVp4AAwv8tDilX1t6ra7pxzNDBSRHKAecC3VbVZVd92Xm+nVjxJZLLzegtUtdZJVNcAX1TVaud9X+vEUqmqf1bVBlWtA37Iuf+3j6nqe86FybeBm5zXfxvwiqq+4ryW14CtwEe8fU+Md6xP1vjTY8DbeK7GH+3jsX0ZAxzv/EZVG0SksvN7EckDfgrMxXN1G4Pn6rg344DPiMiXuxyLc86jQIl+sOppcR/xvauql57n/lNdvm5wztNXHJ2Od/n6vO8D8ALwvyIyAcgHzqjq5j5iP4eI3AF8FRjvHBpCl+6qrq/HiaHrY6qdD/VOxXhaIuD5nRgLPCUiaXhaq99yjlWpanUPsSThaaVcDQxzDqeISLSTzOCD71ExnouBDDzv7ydF5GNd7o8F1pz/HTC+spaJ8RtVLcYzEP8R4Lm+Ht7H/SeBrv3kiXiucDv9GtgH5KrqUOABQOjdceCHqprW5Zakqn90zpXVpdsGPN0rg+F8cXTq+t6c931Q1SbgGTxX5LczgFaJiIwDfgvcA6SrahrwHud/P7vGN0w842Cd3n/vnBbH91R1KrAQ+Cie1utxYLiTYLr7Gp7EuMD5v+3stuwaz9guX+fgaQFVOM/7WLf3N1lVH+zHazE+sGRi/O1zwIe6Xan25DSQLiKpvdz/LPAxEVnojGn8Kx/8MEkBaoF6EZkC/L8enn9il+9/C3xRRBaIR7KIXCsiKcBGPP3uXxGRWBG5gQ928fjT+eLoSV/vA3hagXcCH6fvZBLlDIp33uKBZDwJrBxARP4BmNafF+NcQGwFvicicSJyKfB+q0BElorIdKcLqhbPh36Hqp7EMwnjf5wB91gR6UwaKXjGSWqcgfXv9nDq20RkqtOK+T7wrNNqeRzP+7VMRKKd17ik6wC+GRyWTIxfqeohVd3aj8ftA/4IHHbGHMZ0u3838GXgKTxXv/VAGZ7xBfAMrH4KqMPzAf10t1P8K/CI89w3OTH9I56B3GrgIM6guDMWcYPzfRVwM323rC6Rc9eZzOvH6+41jl4e39f7gKquBzqAbc6H+/nciueDuvN2SFX3AP+JJ6meBqYD6/t6LV18Cs8AfRWeD/6uXZyj8CTEWjzjaGv5e8K7HU9y2ee8pn9yjv8cSMTT0ngX+FsP53wM+D883W8JwFcAVPU4cB2elmo5npbKCuyzbtCJbY5lQoGIDAFq8HRrHXE5HNf09j6IyJvAk6oaFBUEBpOIvAU8HgmvNZRYtjZBS0Q+JiJJTn/8T4BdwFF3owq8vt4Hp0U0m3NbZ8YEjCUTE8yuA0qdWy5wi0ZmU7rX90FEHgFeB/7JmUZrjCusm8sYY4zPrGVijDHGZxG7aDEjI0PHjx/vdhjGGBNSCgoKKlT1nPI0EZtMxo8fz9atfc5gNcYY04WI9Dj93Lq5jDHG+MySiTHGGJ9ZMjHGGOOzkEgmInK1iBSJyEERua+H++8Uz0Y6252bX7Z+NcYY0z9BPwDvFIj7FZ49Mk4AW0TkRaeeUFdPq+o9AQ/QGGNM8CcTPNVbD6rqYQAReQrPiuDuycQEmVWFJaxcXURpTSNj0hJZsSyf5bOy+v5BY0zICYVuriw+uBHOCT64fWunT4hnz+1nRWRsD/ebAFpVWML9z+2ipKbRs/NUTSP3P7eLVYW+bL5ojAlWoZBM+uMlYLyqXgS8hmdb0XOIyF0islVEtpaXlwc0wEizcnURja3tHzjW2NrOytVFLkVkjBlMoZBMSvjgrmrZdNtb3NkzunN/h4eBOT09kar+RlXnqurczMxzFnAaPyqtafTquDEmtIVCMtkC5IrIBGenuVuAF7s+QERGd/n243g24TEuGjE0vsfjY9ISAxyJMSYQgn4AXlXbROQeYDUQDfxeVXeLyPeBrar6Ip7tVj+OZ+vVKs6zc50ZfJX1zXT0UI06ITaKFcvyXYjIGDPYgj6ZAKjqK8Ar3Y59p8vX9wP3Bzouc67aplY+84fN1Da2cc/SSTxfWEqJ07V189yxNpvLmDAVEsnEhIam1nY+/8hW9p2s47d3zGXplBF8fdkUOjqUxSvXsP90vdshGmMGSSiMmZgQ0NrewZee2MaWo1X8500zWDplxPv3RUUJt87PYePhSg6VW0IxJhxZMjE+6+hQvv6nHby5r4wfXDeN62ae25V109yxxEQJf9x0zIUIjTGDzZKJ8Ymq8t0Xd/PC9lJWLMvntovH9fi4zJR4lk0bxbPbTtDUbf2JMSb0WTIxPvnPV/fz2LvF3LV4Il9aMum8j/30ghxqGlp5ZdfJAEVnwtmqwhIWPfgmE+57mUUPvmnVFVxmycQM2G/fPsx/rznILfPGcv81UxCR8z7+konpTMxI5gnr6jI+snI9wceSiRmQZ7Yc54ev7OXa6aP54fXT+0wkACLCpxbkUFBczb5TtQGI0oQrK9cTfCyZGK/9dddJ7ntuJ5flZvDTm2cQHdV3Iun0idnZxMVE8aS1TowPrFxP8LFkYryy7kA59z61nZlj03jo9jnEx0R79fPDkuP46PTRPLethLPNbYMUpQl3Vq4n+FgyMf1WUFzNXY8WMDEzmT/cOZ+kuIGtef30xTnUN7fx0o5SP0doIkFBcRX1TT1fiHzh8okBjsZ0smRizqvrjJkbf72B5LgoHv3cfFKTYgf8nLNzhpE/MoUnN1tXl/HOC9tLuPW3mxgxNIEHPjKFrLREBM/U85goeHF7KS1tHW6HGZGsnIrpVeeMma4DnXXN7Ww4WOlTjS0R4dMX5/CdF3az80QNF2Wn+SFaE85UlV+8cYCfv36A+ROG89BtcxiWHMddi/8+Hf2lHaV8+Y+FfO+l3fzw+ukuRhuZrGVietXTjJnmtg6/zJhZPiuLxNhoG4g3fWpqbeefn97Oz18/wCdmZ/P45xYwLDnunMd9bMYYvnj5JJ7YdMx+r1xgycT0ajBnzAxNiOW6mWN4YXsptU2tPj+fCU+V9c3c9vAmVjkVFn7yyYuIi+n9Y2vFsnwW52Xy3RffY+vRqgBGaiyZmF71NjPGXzNmPr1gHI2t7bxgC81MDw6W1XP9/2xgV8kZfvWp2dy9dHKf65mio4T/umUWWWmJfPHxbZw60xSgaI0lE9OrFcvyz1lDkhgb7bcNrqZnpzI9K5UnNh1De9hMy0Sudw5UcP3/rKehpZ2n7rqYay8a3fcPOVKTYvnNHXNpbGnjC48XWC24ALFkYnp13cwxJMdFkRgbhQBZaYn86Ibpft3g6tMLcth3qo5tx6r99pwmtP1x8zE+84fNZKUlsuruhczKGeb1c+SNTOGnN89kx/Ea/mXVe3axEgA2m8v0quh0HbVN7fz4Exdx07yxg3KOj80Yww9f3ssT7x5jzrjhg3IOE7xWFZawcnURpTWNjE5LIH9kCmuKylmSn8l/3TqLlISBT0FfduEo7r0il1+8cYBpY4Zy56IJfozcdGctE9OrtUXlACzOyxy0cyTHx3D97Cz+susk1WdbBu08Jvh0L9ZYWtPEmqJyLpuczsN3zPUpkXS694pcrrxgJD94eS8bDlX4HrTplSUT06u1+8uZMiqFUakJg3qeTy3IoaWtgz9vOzGo5zHBpaep5wCHKxqIifbPR1NUlPCzm2cwISOZu5/YxvGqBr88rzmXJRPTo7PNbWw5WsXlg9gq6TRl1FDmjBvGkzYQH1ECVawxJSGW39w+h7YO5QuPFdDYYgPyg8GSienRxkOVtLZrQJIJeAbiD1ecZePhyoCcz7hvsKeedzUxcwi/vGUWe0/V8o0/77SLlkFgycT0aO3+cpLiopkz3vuZNAPxkemjSUuKtY2zIsiKZfnEdevO8ufU8+6WThnBimX5vLSjlN+8fXhQzhHJbDaXOYeq8tb+MhZOSve6xPxAJcRGc+PsbB7ZeJTyumYyU3ouMW7Cx/JZWfx52wnWHahA8LRIVizL9+vU8+7+3+WT2F1Sy3/8bR9VDS38ZcdJSmsaA3LucGctE3OOo5UNHK9qDFgXV6dbF+TQ2q78qeB4QM9r3KGqHCqr58oLRnLkwWtZf9+HBv3DXERY+cmLGJkSz0NrD9u2v35kycScY21RGQCX540I6HknZQ7hkonpPLnpGB0d1qcd7naX1lJ6pomrLhwZ0PMmxcXQU5F62/bXN5ZMzDnW7i9nQkYyOelJAT/3py/O4UR1I28fKA/4uU1gvbrnNFECV0wJ7EULQFltc4/HbdvfgbNkYj6gqbWdjYcrA97F1emqqaPIGBJnA/ER4LU9p5k7bjjpQwI/PhbImWSRwpKJ+YAtR6toau1wLZnExURx09yxvLmvjJNn7CoxXB2vamDvyVo+PDWwXVydVizLJzH2g5NLBnMmWSTwOZmISJSILBGRK0QkMFN/zKBZW1ROXEwUCya6Vyfr1vk5dKjy9BYbiA9Xr+05DeBaMlk+K4sf3TCdMU51h8TYaL8XMY00A0omIpIkIjeIyKNAGfAm8BpQLiKPi8iNIjLEn4GawFi7v5wFE4aTFOferPGxw5NYnJvJU5uP09Zu+3mHo1f3nCJv5BDGZyS7FsPyWVlsuP8Kbr94HB2qLM0P/NhNOOl3MhGRkSLyeRH5C1ABPAtcDPwBWAwsAh4G5gLP4Eksr4jIXSIyyv+hG38rqWnkQFm9a11cXU3OTOZUbROTv/VXFj34pk3ZDCPVZ1vYcrSaq6YGx8fCzfPG0tzWwQs77HfMF960TEqBh4BM4AfANFXNU9UVqvqOqm5U1W+o6hRgKvA9IA34NXBCRHwvAWoG1dv7PTOoluS7m0xWFZbwxOa/D8DbGoDw8ua+Mto71LUuru6mZaVy4ZihPLXZulV94U0y+RKQraoLVPVHqrqntweq6j5VfVBVFwJZwP8DrLpakFtbVE5WWiKTMt3toVy5uoim1g92b9kagPDx2p7TjBqawPSsVLdDed8t88ay52Qt75WccTuUkNXvZKKqD6nqSW9PoKqnVPW3qjrgzm8RuVpEikTkoIjcd57HfUJEVETmDvRckaq1vYP1BytYnJfZ5z7bgy1Q1WRN4DW1trN2fzlXTh1BVJS7v2ddfXxmFvExUTy1xaakD1TQTw12Zoj9CrgGT/fZrSIytYfHpQD3ApsCG2F42FZcTV1zW1CMl9gagPD1zoEKGlvbg2a8pFNqYizXTh/NC4WlVqJ+gHxKJiKSJyLXi8gXnIH260Uk11/BOeYDB1X1sKq2AE8B1/XwuB8A/wE0+fn8EWHt/nJiooSFk9PdDqWXNQBRtgYgDLy25zQp8TFcPNH937Pubpo3lrrmNl7Z5XUHjGEAVYNF5ALgi8CNQOflRWd7VZ3HnMYzo+shVd3rY4xZQNeRsRPAgm4xzQbGqurLIrLiPLHfBdwFkJOT42NY4WXt/nJmjxvGUD9sleqrzrn+K1cXUeJ0bX3lilxbAxDi2juU1/eeZsmUEcTFBF+nyIIJwxmfnsTTW47ziTnZbocTcryZGjxJRJ4F3gM+B+zEM2PrDuAjwLXO198HdgCfB94TkT+JyER/B94lrijgp8DX+nqsqv5GVeeq6tzMTPe7c4JFWV0Tu0trg6KLq9PyWVmsv+9DvPPNpQABK4VvBk/hsWoqz7YEzSyu7kSEm+flsPloFYfK690OJ+R4c3mwBxgP3AmMVNVrVPX7qvqEqv5NVf/qfP09Vb0GGAl8Fpjk/OxAlQBju3yf7RzrlAJMA94SkaN41r68aIPw/bdufwVAUCWTTtnDkhiXnsSGQ7YDY6h7dc9pYqPF9ann5/OJOVlERwnPbLVpwt7yJpl80rmqf0xVz/b1YFU9q6qPqOps4OaBh8gWIFdEJohIHHAL8GKX85xR1QxVHa+q44F3gY+r6lYfzhlR1u4vJ2NIPFNHD3U7lB4tnJTBpsOVtho+hKkqr+4+xcUT04OiK7U3I1ISuGLKCP5ccIJW+33zijdTg1/s+1G9/uwLPvxsG3APsBrYCzyjqrtF5Psi8vGBPq/xaO9Q1h0oZ3FeRlBN1exq0eR06prb2GVrAELWwbJ6jlY2cNWFwTWLqye3zB9LRX0Lb+wtczuUkDKQAfgheLqVUoE64LCqnvJ3YF2p6ivAK92OfaeXxy4ZzFjCza6SM1Q3tAZlF1enS5yZPxsOVTIrJzB70hv/erWzsOMFwTle0tXi3ExGDo3n6S3HuHpa8Ce/YOHNAHy8iPwaKAfW4/lwXweUiEiJiDwqIteK2yvejFfWFpUjApflBm8ySR8Sz5RRKaw/WOF2KGaAXt1zmhnZqYxyqvQGs5joKD45Zyxr95fbQlkveDNm8hPgC8BbwAPA151jAiQAt+EZy3hPRK7wb5hmsKzdX8ZF2WkMT45zO5TzWjQ5g63F1TS12oKyUHO6tokdx2uCdhZXT26aO5YOhWcLTrgdSsjwJpncDPzemcX1H6r6M+DHzn2fBCYC9wGJwGoR+Wf/hmr8raahhe3Ha4K6i6vTosnptLR1sK242u1QjJc69y4JhfGSTjnpSSyanM4zW4/T0aFuhxMSvEkmicDG3u5U1aOquhLIB34G/EREPuxjfGYQvXOwgg4NzinB3c2fkE50lLD+kHV1hZrX9pxmXHoSuSNCa4ujm+flcKK60aal95M3yWQrsLSvB6lqq6quAF4AvjXQwMzgW1tUTmpiLDOyg6d6a2+GxMcwIzuV9QftDzuU1DW1suFQBVdNHel6AVFvXTV1JGlJsVb8sZ+8SSYPAp/yovvqFWCO9yGZQFBV1u4v59LcDGKig6+0RU8WTc5g54kaapta3Q7F9NPa/eW0tisfDrLCjv2REBvN8plZvLr7NFVnW9wOJ+h5s85kNbACT/fVZhG5Dc/q895cDjT4GJ8ZJPtO1VFW1xwSXVydFk7KoENh8+Eqt0Mx/fTq7tMMT45jzrjQnNJ987yxtLR38LxtzNYnry5JVfU/gWVAMvAosB9Pcccvicg3ReRLInKfiKwBPgU84e+AjX+sdXZVDKVkMntcGvExUTZuEiJa2jpYU1TGFVNGEB2kC2L7csHoocwYm8bTW46hagPx5+P1okVVfV1EpgFX45nhtQS4wbl1agD+C/imH2I0g2BtUTlTRqUwcmjwz/vvFB8Tzbzxw9lg4yYhYdORSuqa2kJqFldPbpk3lvuf28X24zW2aPY8BtRZrh5/VdU7nXpYmXjKwl8BzAMyVfWfVNU6t4NQfXMbW4uruDyIC+71ZuHkdIpO11Fe1+x2KKYPr+05TUJsFJdOznA7FJ98bMYYkuKieXqLFX88H7+MvKpqpapuUdU1qlqgqrZsNIhtPFRJa7uGVBdXp0WTPB9MGw9b6ySYeQo7nmZxbiaJcaG9fcCQ+BiunT6aF3eUUt/c5nY4QcubcioDXtUuIlcO9GeN/63dX0ZSXDRzxw13OxSvTctKZWhCDBustEpQ21VyhlO1TSG16v18bpk/loaWdl7eWep2KEHLm5bJ30TkTRH5qLMv+3mJSKyzje9auhVpNO5RVd4qKmfhpIyg3O2uL9FRwsUT020QPsi9tuc0UQJXhEBhx/6YnTOMySOG8JR1dfXKm0+TWUAbnvpbpSLyhIjc6ySXhSKySEQ+JiJfFZFngFPAs3gG42f6PXIzIEcqznKiujEkx0s6LZyUzvGqRo5X2czzYPXq7tPMHT886Gu+9ZeIcMu8sRQeq2H/6Tq3wwlK3qwzeU9VrwIWAa8CH8NTNuUFPNWD3wZW4Sn+eJVz/GKnlpcvOy0aP3p/SnAQVwnuyyJnQNeqCAen4sqzFJ2u46ow6eLqdP2sLGKjxQbiezGQqcEbgY1OV9ccYCqe2VyKpzz9e0Chqto2ZUFo7f5yJmYkk5Oe5HYoAzZ5xBBGpMSz4VAlt8zPcTsc0837hR1DcNX7+aQPieeqqaN4btsJvnF1PvExoT2xwN+8TiadVLUd2OzcTAhoam3n3cOV3DIvtD+ARYSFk9J552AlqhpyNZ/C3at7TjNlVEpIX7D05qZ5Y3l510le23Oaj140xu1wgopXI7AiktPlli0i5yunYoLM5iNVNLV2sCSEx0s6LZyUQUV9M/tP17sdiumi6mwLW49Whc0sru4unZxBVlqidXX1wNvpPEeBI86tGKgRkVMi8oiIzPB3cMa/1u4vJz4mioudbXBD2cLJntdg4ybB5Y29p+nQ8Ovi6hQdJXxybjbrDlTYBJBuvE0mb3e5rQN24in2eDuwRUTu9W94xp/W7i9nwcR0EmJDv683e1gS49KT2GBThIPKq3tOMzo1gWlZQ90OZdB8cu5YAD7yi3VMuO9lFj34JqusEKR3YyaquqT7MRGJAZbj2XXxpyJSqaqP+yU64zcnqhs4WFbPrWE0YL1wUgZ/2VFKW3tHyJTRD2eNLe2sO1DOTXPHhvU41pYjVUQJ1Dmr4UtqGrn/uV0ALJ+V5WZorvL5L1BV21T1WTw1uQ4BPxORWJ8jM36zqrCEj/7yHQAeWnsobK6iFk1Op665jV0lZ9wOxQDrDpTT1NoRtuMlnVauLqL7Tr6Nre2sXF3kTkBBwm+Xc6paCXwXGA7M9dfzGt+sKizh/ud2UdPoqblZVtfM/c/tCouEcokz9mPbqgaH1/acJiUhhgUTQn9M7nxKa3ouPdjb8Ujh7WyuDhFp7+0GdHZvveMcs6poLlu5uojG1vYPHAuXq6j0IfFMGZVig/AuW1VYwsIH3+BPBSdoa1de2XXS7ZAG1Zi0RK+ORwpv15m8gmdxYm/S8ZSi34JnAaNxWbhfRS2anMFj7xbT1NoeFhMLQk1ny7fzgqWxtT3sxw9WLMv/wGsGSIyNZsWyfBejcp+3A/AfPd/9zv7wC4C7VHWnL4EZ/xiTlkhJD4kjXK6iFk1O53fvHGFbcTULQ3zfjFB0vpZvuCaTztf14F/3caq2iaEJMXz/umlh+3r7y29jJiJyIfBtoMgSSfBYsSyf6G4za8LpKmre+OFER4lVEXZJuLd8e7N8VhbvPnAF+SNTmJ6dGvGJBPyQTJyV8PcDG4FE4C6fozJ+s3xWFikJMSTGRiFAVloiP7phetj88qckxDIjO5X1tpWvKyJ9/GBxXgZbjlTT0GLDw94OwB/ucjsiIpV4VsL/EGgHblTVdwYjUDMwp840UdPYyteXTeHIg9ey/r4PhU0i6bRocgY7T9RQ22S7RAfaimX5xEWHb8u3L5flZtLS3sGmI1Vuh+I6b1sm47vccpyf3wp8H5iiqi/7MTbjB9uOVQMwZ9wwlyMZPAsnZdChsPmw/UEH2vJZWSxwpmiHY8u3L/MnDCc+Jop1+62b1dsBeFtmHGK2Hq0mPiaKqaPDt7zF7HFpxMdEsf5QBVeG+YK5YFRZ38LCSek8+Y8Xux1KwCXERjN/wnDePmCTVy05hLmCY9XMGJsWklv09ld8TDTzxg9ng42bBNyZhlb2nqpl/oThbofimsW5mRwsqw/7SQd9Cd9PGENTazu7S86EdRdXp4WT0yk6XUd5XbPboUSUrcVVqBL2q97PZ3GeZ0uHdw5EdldXv7u5RORNP5xPVfUKb39IRK4GfgFEAw+r6oPd7v8icDeeSQD1eNa5RPxWwTtPnKGtQ5mTE/7JZNGkDKCIDYcquG5mZPTXB4PNR6qIjRZm5aS5HYpr8kZ6dv5ce6Ccm+aNdTsc13jTMonCM8bmy83rlpCzPfCvgGvwbBF8q4hM7fawJ1V1uqrOxKle7O15wtHWYs+A9OwIaJlMy0olJSGGjVanK6A2HaliRnZaRFcfEBEuy81k/cEK2rtXgIwg/W6Z9FR+PkDmAwdV9TCAiDwFXAe83/JQ1douj0/m/CVfIsa24momZiYzPDnO7VAGXXSUcPHEdFu8GEBnm9t4r+QMdy2e6HYorlucl8Gft53gvZIzzBib5nY4rgiFMZMsoOsemSecYx8gIneLyCE8LZOv9PREInKXiGwVka3l5eE9+0JVKSiujogurk6LJqVzvKrRdsALkMJjNbR1aEQPvne6dHIGIvD2/vD+XDkfvycTp3pwwKnqr1R1EvBN4F96ecxvVHWuqs7NzAz9fdDP53DFWaobWpk7PoKSiVOby6oIB8bmI5VESXivYeqv9CHxTBuTyroIHoQfjJaJv7dYKwG6jmplO8d68xSenR8jWkFx+C9W7G7yCM9A6HobNwmITUequHBMKikJthcewGW5GWw7Vk1dhFZi6DOZiMgGEZnpxXP6e7xiC5ArIhNEJA64BXix6wNEJLfLt9cCB/wcQ8jZVlxNamIsEzOGuB1KwIgICyels/FQBao2bDaYmtvaKTxeY11cXVyWm0lbh0bsJJD+tEwuBl7vLaGIyHV+jagbVW0D7gFWA3uBZ1R1t4h8X0Q+7jzsHhHZLSLbga8CnxnMmELB1uJq5owbRlRU+O7F3ZOFkzKoqG9h/+l6t0MJaztPnKGlrcOSSRdzxg0jKS46Yru6+jOb6yzwazwJ5QpV3dHt/seAQa3Voaqv4NmYq+ux73T5+t7BPH+oqWlo4WBZPddHSH2krhZO9iyeW3+wgvxRKS5HE742O4UN5423ZNIpLiaKSyamR2xplf60TFRVvw38L/CGiMzodn9kXfqGgMJjNQDMjqCZXJ2yhyUxLj2JDTZFeFBtOlJF3sghETHt3BuX5WZQXNlAceVZt0MJuH4PwKvqvwAP4WmhXNT1Lr9HZXyytbiK6ChhZoTOd184KYNNh6toa+9wO5Sw1NbeQcHRKuvi6kFnaZVI7OryajaXqn4L+A2eFsr0wQnJ+KqguJoLxwwlMS4yVyUvmpxOXXMbu0rOuB1KWNpzspazLe3Mj+B6XL2ZkJFMVlpiRK436U8y+UA3lpNQfgu82a2FYoJAa3sHO46ficgurk6XOPtrbIjQWTWDrXO8ZL6Nl5xDRFicl8nGQ5W0RljLuD/J5L7uB1T1ATwJ5Q0g3t9BmYHbe7KWxtb2iFqs2F36kHhGD43nl28cYMJ9L7PowTdZVXi+pUnGG5uOVDEuPYlRqQluhxKUFudmUNfcxvbjNW6HElB9JhNV/VUvxx8AHgZsxVIQicTFit2tKiyhrL6F5rYOFCipaeT+53ZZQvGDjg5ly9Eqa5Wcx8JJGUQJrIuwri6fVsCr6v1A99ldxkUFxdWMSU1gdGqi26G4ZuXqonOqtza2trNydZFLEYWPA2X11DS0vr9VrzlXalIsM8em8XaEDcL7XE5FVXf5IxDjHwXF1RFRcv58etvxLtJ3wvOHTUc841ALbCbXeV2Wm8nOEzXUNLS4HUrAhELVYNNPpTWNnDzTxNwITyZj0npulfV23PTfpiNVjE5NIHuYvZfnszgvkw6F9RG0lfRANqta7GxYZYLM38dLIvuqccWyfBK7bdaUGBvFimX5LkUUHlSVzUc860tEbK3y+czI9mzWti6CVsMPpGWyBkjtzwNF5EYR+bmIfFZEYrrd9/IAzm3Oo6C4msTYaKaMjuwyIstnZfGjG6aT1aUlcvsl41gegeVl/OloZQPldc22WLEfYqKjWDQpg7f3l0dM0dGBJJN+XZKIyD3AfwNJwApgvYh0/S28bADnNudRUFzNzLFpxEZb7+XyWVmsv+9D7P+3axg7PJH1ByvpiOAtVf1hs42XeGVxXialZ5o4VB4ZpVUG81PnHmCZqt4FTAe241no2PmbaO1kP2poaWPPydqInhLck7iYKP75yjx2l9byynsn3Q4npG06UsXw5DgmZUbOtga+uCzXs1lbpHR1DWYyGdNZYVhV21T1C8CbwBoRScdqevnV9uM1tHeoJZMeXDczi7yRQ/jpq/utXpcPNh/xrC+x8ZL+GTs8iQkZyRFTWmUwk0m5iEzoekBVv4pnzGUN/St/b/ppmzP4HsllVHoTHSV8/ap8Dlec5dmCE26HE5JKaho5Ud1o4yVeWpybwbuHq2huc2U384AaSDLpb4viTeDOc35Y9Z+AtwCrxeBHBcXV5I4YQmqSFSToyYenjmRWThq/eOMATa3h/4ftb1s663FZMvHKZbmZNLa2U3C02u1QBt2gDcADdwP/0dMdqvoVYPwAzm160NGhFDg7K5qeiQgrluVz8kwTj79b7HY4IWfTkSpS4mO4YPSg7oMXdi6ZlE5stETEaviBJJOngea+HqSqLaracJ77jw3g3KYHh8rrqW1qs2TSh4WTMrgsN4NfrTlIXVOr2+GElM1HKpk7fhjREbYNtK+S42OYnTMsIgbhvU4mqnqrqkbGXLcQYcUd+2/FsnyqG1p5eN0Rt0MJGRX1zRwqP2v7lwzQ4rxMdpfWUl7X5zV4SLMFCWGgoLia4clxTMhIdjuUoHdRdhrXTBvFw+sOU1kf3n/c/mLjJb5ZnOvZfXH9wfDu6rJkEgYKiquZnTPMpmz209euyqOxtZ1frTnkdighYdORKhJio5ie1a/CF6abC8cMZXhyHG+HeVeXJZMQV3W2hcMVZ62LywuTR6Rw45xsHn+3mBKrJNynzUeqmJ0zjLgY+7gYiKgo4dLJGaw7UBHWpVV8+u0QkTgROeyvYIz3ttl4yYDce2UeAL94fb/LkQS3M42t7D1Va11cProsN4Pyumb2napzO5RB4+ulhmBTfF21tbia2GjhomzrgvBGVloit108jmcLTnCwrN7tcIJWQXEVqrDABt99cpkzbhLOq+H7TCYicri3G1CElUVx1bbiai4ck0pCrO0K4K27l04iMTaan75mOzD2ZtPhKmKjhVk5aW6HEtJGpSaQPzKFdWG83qQ/JU0ygQeAntaFxAFP+TUi028tbR3sOFHDbRePczuUkJQ+JJ7PXTaRX75xgF0nzjDdWnfn2HSkihnZaXax4geX5Wbw6LvFNLa0kxgXfu9nf7q5tgOnVfWF7jfgRaz6r2t2l56hua3Dxkt88I+XTWBYUiw/Xr3P7VCCztnmNt4rOWPjJX5yWV4mLW0d7299HG76k0x+CfT26luBf/BfOMYbtljRdykJsXxpyWTWHahg46Hw/CMfqMJjNbR1qCUTP1kwYThxMVFh29XVZzJR1T+p6hu93Nehqo90O2wtlQDZdqya7GGJjBxqNTN9cfsl4xg1NIEfr94X1lM3vbX5SCVRYhcr/pIQG82CCcPDtrSK3yeOq6pNRg8AVWXrUSvu6A8JsdHce2UuhcdqeH1vmdvhBI1NR6q4cEwqKQlWidpfFudmsv90PSfPhN/6JvvgD1Enqhspq2tmriUTv/jknGwmZCTzk9VFtNv2vjS3tVN4vMa6uPzssrzO3RfDr6ur3xtUicibfjifquoVfnieiLftmLMZliUTv4iJjuKrH87jy38s5MUdJVw/K9vtkFy188QZWto6LJn4Wf7IFEakxLPuQAU3zR3rdjh+5U3LJArPeIgvN2sJ+cnWo9Ukx0WTPzLF7VDCxrXTR3PhmKH89LX9tLRF9va+m53ijvPGWzLxJxHhstxM3jlQHnYt4H63TFR1ySDGYbxUUFzNrJxhxERbfvaXqCjPBlp3/mEL8374OrWNrYxJS2TFsnyWz8pyO7yA2nSkiryRQxieHOd2KGEnKS6K6oZWJj/wSlj9foXEJ5GIXC0iRSJyUETu6+H+r4rIHhHZKSJviEhYr+Krb25j36la6+IaBNVnW4gST00qxbP3+f3P7WJVYYnboQVMW3sHBUerrItrEKwqLOFPBScAwu73K+iTiYhEA78CrgGmAreKyNRuDysE5qrqRcCzwI8DG2VgbT9WQ4falM3B8JNX99O996GxtZ2VqyOn5Mqek7WcbWm3zbAGwcrVRTS1frALNVx+v/rdzdUTEckDLgRG4Em05cB7qnrAD7F1mg8cVNXDzjmfAq4D9nQ+QFXXdHn8u8Btfjx/0CkorkYEq5c0CEp7KUnf2/Fw1DleMt/GS/wunH+/vE4mInIB8EXgRmBU52HnX3Uecxp4BnhIVff6GGMWcLzL9yeABed5/OeAv/Z0h4jcBdwFkJOT42NY7ik4Vk3+yBSG2vx/vxuTltjjHidj0hJdiMYdm45UMS49iVGpthjW38L596vf3VwiMklEngXew/OBvRP4HnAH8BHgWufr7wM7gM8D74nIn0Rkor8D7yXG24C5wMqe7lfV36jqXFWdm5mZGYiQ/K69QyksrrbxkkGyYlk+id2KGibGRrNiWb5LEQVWR4ey5WiVtUoGSTj/fnnTMtkD7ALuBJ5T1bPne7CIJONpvdzr/OxAL3NKgK4TsrOdY93PdyXwLeByVQ3bzb0PlNVR19xmixUHSeesmh+v3kdpTROJsdH86IbpYTHbpj8OlNVT09DKgok2XjIYOn+PVq4uoqSmkSiBf18+LSx+v7wZgP+kc1X/WF+JBEBVz6rqI6o6G7h54CGyBcgVkQkiEgfcgqda8ftEZBbwEPBxVQ3rehhW3HHwLZ+VxYb7ruBjM8YwJCGG62aOcTukgNnsVLRdYDO5Bs3yWVmsv+9D/OzmGXQoTA6TtWL9Tiaq+mLfj+r1Z1/w4WfbgHuA1cBe4BlV3S0i3xeRjzsPWwkMAf4kIttFZMCxBruCo9VkDIkjZ3iS26GEvcXOVqt7T4bvVqvdbTpSxejUBLKHhX4ffrBbnJuJCKwpCo/r34EMwA8BpgGpQB1wWFVP+TuwrlT1FeCVbse+0+XrKwfz/MGk4Fg1s3OGIWLFmQfb4jxnq9UD5UwdM9TlaAafqrLpSBULJ6Xb71cApA+JZ0Z2GmuKyvjKFbluh+Mzbwbg40Xk13im/67H8+G+DigRkRIReVRErhX7LRw05XXNFFc2MHe8dXEFwsihCUwZlRLW+3Z3dbSygfK6ZlusGEBL80ew/XgNVWdb3A7FZ96MmfwE+ALwFp5tfL/uHBM8g+u34RnLeE9ErJjjIOgs7mjjJYFzeV4mW45Wcba5ze1QBtWqwhKu/9V6AH7x+oGwWJEdCpZOyUSVsLhg8SaZ3Az8XlWvUdX/UNWf8feV5p8EJgL3AYnAahH5Z/+GGtlWFZbwtWd2APDlJwvtjz1AFudl0tquvHs4fHdhXFVYwv3P7aKmsRWAsrrmsCnxEeymjUklY0hcWIybeJNMEoGNvd2pqkdVdSWQD/wM+ImIfNjH+Ax//2Ovd66OS8802R97gMwdP4zE2OiwuHLszcrVRTS2tn/gWLiU+Ah2UVHC5XkjWLs/9KsIe5NMtgJL+3qQqraq6grgBTzrPoyP7I/dPfEx0Vw8cThvh+FmRp3CucRHKFg6JZOahla2H692OxSfeJNMHgQ+5UX31SvAHO9DMt3ZH7u7FudlcqTiLMcqG9wOZVD0VsojHEp8hILLJmcSHSWs2RfarV9v1pmsBlbg6b7a7JQuOd9qm8uB8PzrCzD7Y3dX5xThtQdC+4+9NyuW5RPdbRJmuJT4CAWpSbHMyRkW8uMmXpWgV9X/BJYBycCjwH48xR2/JCLfFJEvich9IrIG+BTwhL8DjkT39jAH3f7YA2diRjLZwxLDdtzkupljSIyLIjE2GgGy0hIjqoRMMFgyJZPdpbWcrm1yO5QB83rRoqq+LiLTgKvxzPBaAtzg3Do1AP8FfNMPMUa8xDhPYbiMIXFU1reE1e5soUBEWJyXyYvbS2lt7yA2zHa33H+6nvrmdv7jE9O5eV7oVtMOZUvzR/DjvxWxtqicm+aF5t7wA9rPRFUVT5n3vwKISDqeqcFDgFpgj6pah76frCosYdTQBNbf9yGio2xNqBsW52by5KZjbCuuDrsiiOsPeiYXLJyU4XIkkWvKqBRGpyawpqgsZJOJXy6xVLVSVbeo6hpVLbBE4j+V9c2s3V/OdTPHWCJx0cLJ6URHCW+H4bjJhkMV5AxPYqzVe3ONiLAkfwTrDlTQ2t7R9w8EIW/KqQx4VbtTHt4MwF92nqStQ7l+tnVpuWloQiyzc9JYG2bjJm3tHWw6XMWiyeHV2gpFS/MzqW9uY+vR0Jwi7E3L5G8i8qaIfNTZl/28RCRWRK4XkbV0K9Jo+u/5whKmjEphyqjwLzQY7C7Py+S9kloq6sNnu5ydJWeoa25j0WTr4nLboskZxEYLb4XorC5vksksoA1P/a1SEXlCRO51kstCEVkkIh8Tka+KyDPAKeBZPIPxM/0eeQQ4UnGW7cdruN4G2oNC5xThd8JoAeMGZ7zkkjAbBwpFyfExLJiQHrJThPs9AK+q7wFXicglwJeA64BbcfZ970LwDMI/B/xaVbf4KdaI83xhCSJw3UxLJsFg2phUhifH8fb+8rCZSbf+YCUXjB5K+pB4t0MxwJL8TP7t5b2cqG4ge1hojWENZGrwRmCj09U1B5gKZOJJKuV49ogvVNXQHEUKEqrKqsISFk5KZ1TqQHc8Nv4UFSVcOjmDtw9U0NGhRIX4hIim1nYKjlVzx8Xj3A7FOJZOGcG/vbyXt4rKuS3E/l8GNDUYQFXbgc3OzfjZtmM1HKtq4Msfmux2KKaLxXmZvLijlD0na5mWlep2OD7ZerSalrYOGy8JIhMzkskZnsSafWUhl0y8mhosIjldbtkiEh6bFwehVYUlxMdEcfW0UW6HYrpYnOv54A2HKcLrD1UQEyW2GVYQERGW5mey/lAFTd2KuwY7b9eZHAWOOLdioEZETonIIyIyw9/BRaqWtg5e2lnKVReOIiUh1u1wTBcjhiZwweihYVFaZcPBCmaOTSM5fsAdFGYQLJkygqbWDjYdqXI7FK94m0ze7nJbB+zEU+zxdmCLiNzr3/Ai09r95dQ0tHL9rDFuh2J6sDgvg4Li6pDeffFMYyu7Ss6w0Lq4gs4lE9OJj4lizb7QmtXlbaHHJaq61LktUdVZQCpwE3AC+KlTTdj4YFVhCenJcVyWm+l2KKYHl+d6dl/ceCh0d19893AlHQqLJtmU4GCTEBvNwknpIbfexOdyKqrapqrPAvOAQ8DPRMT6ZgaotqmV1/ae5mMzxoRdQcFwMcfZfTGUV8NvOFhBYmw0s3KGuR2K6cHSKSM4WtnAkYqzbofSb377tFLVSuC7wHBgrr+eN9L8dddJWto6wmYdQziKj4nmkknpIT0Iv/5QJfMmDCcuxi5YgtGSvBEAIdXV5e1srg4Rae/tBjzuPPQd51jodiq75PnCEiZkJDMjO7SnnYa7y/MyKa5soLgydK4cO52ubeJgWb11cQWxnPQkJmUmh9RqeG+ncbzCuSveu0oHFgBb8CxgNF4oqWnk3cNVfPXDeYiE9oK4cNdZWuXt/eXcfkmyy9F4Z8MhTwkVW18S3Jbmj+DRjcU0tLSRFBf8M+68ilBVP3q++5394RcAd6nqTl8Ci0QvbC8BYLmVTwl649OTGDs8kbX7K7j9kvFuh+OV9QcrSUuKZepoKx4azJZOGcHD7xxhw8FKrpw60u1w+uS3DlMRuRD4NlBkicR7qsrz20qYM24YOemhVZMnEokIi3Mz2Xiogpa20KkcpKpsOFjBJRPTQ74cTLibN344yXHRIdPV5XMycVbC3w9sBBKBu3yOKgLtOVnLgbJ6qxAcQhbnZXK2pZ2C4tDZf+JIxVlKzzTZ+pIQEBcTxaW5GbxVVI5nc9vg5u0A/OEutyMiUolnJfwPgXbgRlV9ZzACDXerCkuIjRaunT7a7VBMPy2clE5MiO2+uN5ZG2OD76Fhaf4ISmoaOVBW73YoffK2ZTK+yy3H+fmtwPeBKar6sh9jixjtHcoL20tZkj+CYclxbodj+iklIZbZ44aFVGmVDQcrGJ2awISM0Jo0EKmW5HumCL8ZAlOEvV0BH9XlFq2qw1R1gar+q6qeHqwgw92GQxWU1TVbF1cIujwvk92ltZTXBf/uix0dysbDlSyclGGzBUPEqFRPLbhQWG9iK5aCwPOFJaQkxPChKSPcDsV4abFT8mZdCHR17TlZS01Dq+33HmKW5meytbia2qZWt0M5r34nExH5kohke3sCERktIl8QEUtcPWhoaeNv753i2umjSYiNdjsc46ULxwwl3dl9MditP2jrS0LR0ikjaO/QoN8u2psP+F8AxSKyVUT+RUQu6u2BIjJNRL4lIpvwFID8b2DAn5QicrWIFInIQRG5r4f7F4vINhFpE5EbB3oeN7y25zQNLe1WPiVERUUJl+ZmsM7ZfTGYrT9UyeQRQxg51HbuDCWzxqYxNCEm6Lu6vEkmI4HPAseA+4BCZ0bXz0XkQyJyuYj8VEQOATuAB4CTwOeB0ao6oDaasz3wr4Br8GwRfKuITO32sGPAncCTAzmHm54vLCErLZH5422DolB1eV4mlWdb2HOy1u1QetXS1sGWI1U2iysExURHsTgvk7f2lwf1BUu/k4mqVqnqI6p6A5ABLAfeBG4BXne+vg1YC9wAZKjqclX9g6r60j6bDxxU1cOq2gI8BVzXLbajzkLJ0Fk9BpTXNbPuQAXXzRxjC8hCWOdWAcFcRbjwWDWNre22viRELc0fQXldc1BfsAxoHENVm1T1JVX9HDAauBS4HBilqp9V1RdUtdFPMWYBx7t8f8I55jURucvppttaXu7+H/5LO0pp71CbxRXiMlPimTp6aFAnk/WHKokSuHiitUxC0eX5nguWYO7q8sd+JqqqG1T1HVUN6paBqv5GVeeq6tzMTPc3nlq1vYRpWUPJHZnidijGR4vzMtlWXE1dkM642XCwgulZqaQm2lZDoShjSDwzslODurRKKMywKgHGdvk+2zkW0g6W1bPzxBkr6hgmFudl0NYRnLsvnm1uY/vxGuviCnFL8kdQeLyGqrMtbofSo1BIJluAXBGZICJxeMZoXnQ5Jp+tKiwhSuDjM2yf93Awd9xwkuKig7K0yuYjVbR1KIsmWTIJZUunjEA1eNc0BX0yUdU24B5gNbAXeEZVd4vI90Xk4wAiMk9ETgCfBB4Skd3uRdy3jg5l1fYSLs3NZIRN0wwLcTFRLJyUztr9wVeUb/3BCuJiopg73rboDWUXZaWSnhwXtOMmwb/jCqCqr+DZmKvrse90+XoLnu6vkFBwrJoT1Y187ao8t0MxfrQ4L5PX95ZxtLIhqGpfrT9UyZycYbYoNsRFRQmX52WypqiM9g4lOshmgAZ9yyQcPV9YQmJsNFdNHeV2KMaPOkurBNNq+Mr6ZvaerLUSKmFiSEIM1Q2tTH7gFRY9+CarCoNn+NiSSYA1t7Xz8s6TLLtwJMnxIdEwNP00PiOZnOFJQZVMNh72TAiwwffQt6qwhGe2elZJKJ5tvu9/blfQJBRLJgG2Zl85ZxpbuX52yPTKGS9cnpfJxsOVQbP74vqDlaTEx3BRVqrboRgfrVxdRFPrB3+vGlvbWbm6yKWIPsgujQNkVWEJK1cXUVLTSJRAZV2T2yGZQRAXE0VDSzt5//JXstISWbEs39W6axsOVbBg4nBiou26MdSV1vS8Dry344Fmv2EBsKqwhPuf20WJ85/eofCtVbuDpnlq/GNVYQlPbCp+/3u3uyFOVDdQXNnAQpsSHBbGpCV6dTzQLJkEwMrVRTS2tn/gWDA1T41/BFs3xIaDzha9Nl4SFlYsyyexhxl5lwRJ8U5LJgEQ7M1T4x/B9v+8/lAFGUPiyRs5xJXzG/9aPiuLH90wnay0RAQYk5rA1NEpPFtw4v2BeTfZmEkApCXFUt1wbs2mYGmeGv8Yk5b4fldm9+OBpqpsOFTJwknptkVvGFk+K+sDY3BNre3846Nb+eafdxIfE8V1LpZnspbJIHt9z2lqGlrpvr4oMTaaFcvy3QnKDIpeuyEmBn6vmgNl9ZTXNdv6kjCXEBvNb26fy4IJw/nqMzv4666TrsViyWQQbT5Sxd1PbuOi7FR+dP3fm6dZaYn86IbptrtimDmnGyLN0w3x/PbSgG+52rlFrw2+h7/EuGh+95l5zBybxpf/WMjre067EocEWx2hQJk7d65u3bp10J5/T2ktNz+0kRFD4/nTFxcyPDlu0M5lgld9cxuf+J8NnDzTyPN3L2JSZmDGLz7/yFYOlNWxdsXSgJzPuK+2qZXbH97E3pN1/PYzc7k8b3C22RCRAlWd2/24tUwGQXHlWe74/WZSEmJ47HMLLJFEsCHxMTz8mbnERkfx+Ue2UtMw+OXD29o72HS40lolEWZoQiyPfnYBk0cM4a5Ht7LhUGBbw5ZM/KystonbfreJ9o4OHv3cAhtkN4wdnsT/3j6HkupG7n5yG63tg7s6flfJGeqa22y8JAKlJsXy+OcXMC49ic/931a2HK0K2LktmfjRmYZW7vj9ZirrW/i/f5jP5BE2JdN4zBs/nH+/YTrrD1byvZcGd4eEDc4GXZfYFr0RaXhyHE98/mJGpyXwD3/YQuGx6oCc15KJnzS2tPO5R7ZwuPwsv7l9LjPGprkdkgkyN87J5guXT+Txd4/x6Majg3aedw5UcMHooaQPiR+0c5jglpkSz5Ofv5jhyXHc8fvNvFdyZtDPacnED1rbO7j7yW0UHKvm57fM5NJc66s2PfvGsilcecEIvvfSnkGpLtzU2k7BsWoWBcmqaOOeUakJPPmPCxiaEMvtv9vEvlO1g3o+SyY+6uhQvvHsTt7cV8a/LZ/GR6aPdjskE8Sio4Sf3zKL3BFDuPvJbRwsq/fr8289Wk1LW4eVUDEAZA9L4sl/XEB8TDQ3/noD83/4OhPue3lQ9kKxZOIDVeUHL+/h+cISvn5VHp9eMM7tkEwI6JzhFRcdxecf2eK3GV6rCkv44uMFADzwfPDsc2HcNS49mX+4dDz1ze2U1TUP2l4olkx88Ks1B/nD+qN8dtEE7l462e1wTAjJHpbEQ7fPobSmiS894fsMr87K1PXNbQCcPNMUVBsnGXc9uqH4nGP+LkJqtbm80LknSWlNI6mJsdQ0tnL9rCz+5doLrP6R8drc8cP50Q3T+dqfdvDdF3fzw+XTBvx79OO/7eu1MrVVWjCBKEJqyaSfOq/8Ov9gaxo99bYunZxOVPfCW8b00yfmZHOgrJ7/XXuIvBFDuHPRhH79XF1TK9uO1bD1aBWbj1RReqbnzdasMrWBwBQhtWTSTz3tSdKh8NPXDvCJOWNdisqEg28sy+dgWT3/+tIe/uvNg1SdbWFMt10ay+qa2Hq0ms1HqthaXMWe0lo61DOgf+GYoSTHR3O2uf2c57ZFswY8RUi7XgyD/4vNWjLpp2Dbq8KEj6go4coLRvDG3tNUnvUMxpfUNLLi2R08uamYsrpmjlY2AJAQG8WsscO450O5zBs/jFk5wxgSH3NOyxmsMrX5u86Lks5u+u4XK/5gyaSfgmmvChN+/uvNg3QvudrarmwtruaKC0byqQU5zBs/nAvHpBIXc+68mUB8WJjQ1n0vFH+zZNJPgWgmmsjVWwtXFX57xzkFWns02B8WxpyPTQ3up+57VdieJMafemvhWsvXhAprmXjBrvzMYLGWrwl1lkyMCQI25mFCnSUTY4KEtXxNKLMxE2OMMT6zZGKMMcZnlkyMMcb4zJKJMcYYn1kyMcYY4zNR7V7EITKISDlwbpH//skAKvwYjhtC/TVY/O4L9dcQ6vGDO69hnKpmdj8YscnEFyKyVVX7V+MiSIX6a7D43RfqryHU44fgeg3WzWWMMcZnlkyMMcb4zJLJwPzG7QD8INRfg8XvvlB/DaEePwTRa7AxE2OMMT6zlokxxhifWTIxxhjjM0smXhKRq0WkSEQOish9bsfjDRH5vYiUich7bscyUCIyVkTWiMgeEdktIve6HZM3RCRBRDaLyA4n/u+5HdNAiEi0iBSKyF/cjmUgROSoiOwSke0istXteLwlImki8qyI7BORvSJyiesx2ZhJ/4lINLAf+DBwAtgC3Kqqe1wNrJ9EZDFQDzyqqtPcjmcgRGQ0MFpVt4lIClAALA+h/wMBklW1XkRigXeAe1X1XZdD84qIfBWYCwxV1Y+6HY+3ROQoMFdVQ3LRoog8AqxT1YdFJA5IUtUaN2Oylol35gMHVfWwqrYATwHXuRxTv6nq20CV23H4QlVPquo25+s6YC8QMpuAqEe9822scwupKzoRyQauBR52O5ZIJCKpwGLgdwCq2uJ2IgFLJt7KAo53+f4EIfRBFm5EZDwwC9jkcihecbqItgNlwGuqGlLxAz8HvgF0uByHLxR4VUQKROQut4Px0gSgHPiD09X4sIgkux2UJRMTkkRkCPBn4J9UtdbteLyhqu2qOhPIBuaLSMh0OYrIR4EyVS1wOxYfXaqqs4FrgLudLuBQEQPMBn6tqrOAs4Dr47eWTLxTAozt8n22c8wEkDPW8GfgCVV9zu14BsrpmlgDXO1yKN5YBHzcGXN4CviQiDzubkjeU9US598y4Hk8Xdih4gRwokuL9lk8ycVVlky8swXIFZEJzqDXLcCLLscUUZwB7N8Be1X1p27H4y0RyRSRNOfrRDyTOfa5GpQXVPV+Vc1W1fF4fv/fVNXbXA7LKyKS7EzewOkeugoImRmOqnoKOC4i+c6hKwDXJ6DEuB1AKFHVNhG5B1gNRAO/V9XdLofVbyLyR2AJkCEiJ4Dvqurv3I3Ka4uA24FdzrgDwAOq+op7IXllNPCIMzMwCnhGVUNyem0IGwk877kuIQZ4UlX/5m5IXvsy8IRzUXsY+AeX47GpwcYYY3xn3VzGGGN8ZsnEGGOMzyyZGGOM8ZklE2OMMT6zZGKMMcZnlkyMCSIiskxE3hKRehEpF5H/FpEEt+Mypi+WTIwJEiLyNeBvwEngn4GXgLuBX7gZlzH9YetMjAkCInIl8CrwDVX9SZfjfwOWApmhVoPMRBZrmRjjMhGJwtP6KAT+s9vdbwFxQMgUgzSRycqpGOO+ZcBU4E49t6ugxfk3NbAhGeMdSybGuO9moB1YJyIZ3e4b6fxbF9iQjPGOjZkY4zIRKQZy+nhYlqqWBiIeYwbCkokxLnJaIuV49tT4nx4e8gzQrKqjAxqYMV6ybi5j3DXR+XeLqr7e9Q4RmQAMA54MeFTGeMlmcxnjriHOvz2Nidzo/Pt0gGIxZsAsmRjjrs61I0O7HnQ2Pfp/QBHwcqCDMsZblkyMcdceoAHP9OCufgiMB76iqu2BDsoYb9mYiTEuUtUGEXkY+IqIPA6sBa4BrgdWqOqrrgZoTD/ZbC5jXOZ0af0Y+DSQBBQA/x6C+5KbCGbJxBhjjM9szMQYY4zPLJkYY4zxmSUTY4wxPrNkYowxxmeWTIwxxvjMkokxxhifWTIxxhjjM0smxhhjfGbJxBhjjM/+P/HpXTn4dn4/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = plt.figure(2)\n",
    "_ = plt.plot(thetas, results_zne, 'o-')\n",
    "_ = plt.xlabel(r'$\\theta$', fontsize=18)\n",
    "_ = plt.ylabel(r'$\\langle \\Psi(\\theta) | \\frac{1 - Z}{2} | \\Psi(\\theta) \\rangle$', fontsize=18)    \n",
    "_ = plt.title('Mitigated Energy Landscape')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ec6804a",
   "metadata": {},
   "source": [
    "Run optimization with ZNE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "5f2b37ea",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     fun: 0.016000000000000007\n",
      "   maxcv: 0.0\n",
      " message: 'Optimization terminated successfully.'\n",
      "    nfev: 20\n",
      "  status: 1\n",
      " success: True\n",
      "       x: array([-0.04267812])\n"
     ]
    }
   ],
   "source": [
    "def mitigated_expectation(thetas, executable: Program, factory) -> float:\n",
    "    mitigated_exp =  zne.execute_with_zne(quil_prog, lambda p: expectation(thetas, p), factory=factory)\n",
    "    return mitigated_exp\n",
    "\n",
    "quil_prog = qc.compiler.quil_to_native_quil(program, protoquil=True)\n",
    "res_zne = optimize.minimize(mitigated_expectation, init_angle, args=(quil_prog, fac), method='COBYLA')\n",
    "print(res_zne)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e85d131",
   "metadata": {},
   "source": [
    "## References\n",
    "[1] [Parametric compilation tutorial in pyQuil.](https://pyquil-docs.rigetti.com/en/v2.28.1/migration3-declare.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a13da8d",
   "metadata": {},
   "source": [
    "display information about Mitiq, packages, and Python version/platform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "4054ce1b",
   "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.11.4"
   }
  },
  "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
}