examples/fault_tolerant/qsp_state_prep.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_KpZQ0nC_lXP"
      },
      "source": [
        "# State preparation using quantum signal processing\n",
        "\n",
        "We are going to apply a delta function to the Hamiltonian ($H$), such that all states with eigenvalues outside a specific energy band ($\\lambda \\pm \\Delta/2$) are supressed. This is performed by applying the operation\n",
        "\n",
        "$f(H,\\Delta,k)=\\sum_{n=0}^{k}\\frac{T_n(-1+2\\frac{H^2-\\Delta^2}{1-\\Delta^2})}{T_n(-1+2\\frac{-\\Delta^2}{1-\\Delta^2})}$\n",
        "\n",
        "using quantum signal processing, where $T_n$ refers to the Chebyshev polynomial of the first kind, defined by $T_n(cos(\\theta)) = cos(n\\theta)$.\n"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Installation of tangelo if not already installed.\n",
        "try:\n",
        "    import tangelo\n",
        "except ModuleNotFoundError:\n",
        "    !pip install git+https://github.com/goodchemistryco/Tangelo.git@develop --quiet"
      ],
      "metadata": {
        "id": "Aps1JmYB_1an"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "As an example, we will use the polynomial with approximation order $k=20$ for $\\Delta = 0.1$ which is depicted below"
      ],
      "metadata": {
        "id": "cLKMJJNsQa4j"
      }
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "id": "wzsbLOu4_lXR",
        "outputId": "83cd80f9-3f5e-4a21-fb45-5e627a20ad2d",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 447
        }
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7f954d58b310>]"
            ]
          },
          "metadata": {},
          "execution_count": 2
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABjX0lEQVR4nO3deXhTVd4H8O9N2qZ7C7R0gbLv+yYFZB0qi8uACyo6sojg64CjMm7MjMuo74uOjs6MMqjIoqOI4iBuiAsKiBSwZV+lpRtLC23pkrRNmuS8fyT3tqFNm7TZ+X6eJ89Dk5ubc0ly87vn/M7vSEIIASIiIiIfpvJ2A4iIiIiaw4CFiIiIfB4DFiIiIvJ5DFiIiIjI5zFgISIiIp/HgIWIiIh8HgMWIiIi8nkMWIiIiMjnBXm7Aa5gNptx/vx5REVFQZIkbzeHiIiIHCCEQGVlJZKTk6FSNd2HEhABy/nz55GSkuLtZhAREVELFBQUoGPHjk1uExABS1RUFADLAUdHR3u5NUREROSIiooKpKSkKL/jTQmIgEUeBoqOjmbAQkRE5GccSedg0i0RERH5PAYsRERE5PMYsBAREZHPcypgWb58Oa655hpERUWhffv2mDlzJk6dOmWzTU1NDRYvXox27dohMjISt956K4qKiprcrxACTz/9NJKSkhAWFoa0tDScPn3a+aMhIiKigORUwLJjxw4sXrwYe/bswXfffYfa2lpMmTIFOp1O2eaRRx7BF198gY0bN2LHjh04f/48brnllib3+7e//Q3/+te/8Oabb2Lv3r2IiIjA1KlTUVNT07KjIiIiooAiCSFES5986dIltG/fHjt27MD48eNRXl6O+Ph4rF+/HrfddhsA4OTJk+jbty/S09MxatSoBvsQQiA5ORl//OMf8eijjwIAysvLkZCQgHXr1uHOO+9sth0VFRWIiYlBeXk5ZwkRERH5CWd+v1uVw1JeXg4AaNu2LQAgMzMTtbW1SEtLU7bp06cPOnXqhPT09Eb3kZOTg8LCQpvnxMTEIDU11e5z9Ho9KioqbG5EREQUuFocsJjNZjz88MO49tprMWDAAABAYWEhQkJCEBsba7NtQkICCgsLG92PfH9CQoLDz1m+fDliYmKUG6vcEhERBbYWByyLFy/G0aNHsWHDBle2xyHLli1DeXm5cisoKPB4G4iIiMhzWhSwLFmyBF9++SV+/PFHm9r/iYmJMBgMKCsrs9m+qKgIiYmJje5Lvv/KmURNPUej0ShVbVndloiIKPA5FbAIIbBkyRJ8+umn+OGHH9C1a1ebx4cPH47g4GBs27ZNue/UqVPIz8/H6NGjG91n165dkZiYaPOciooK7N271+5ziIiI6OriVMCyePFivP/++1i/fj2ioqJQWFiIwsJCVFdXA7Akyy5YsABLly7Fjz/+iMzMTMyfPx+jR4+2mSHUp08ffPrppwAs6wc8/PDDeOGFF/D555/jyJEjmDNnDpKTkzFz5kzXHSkRERH5LacWP1y5ciUAYOLEiTb3r127FvPmzQMAvPbaa1CpVLj11luh1+sxdepU/Pvf/7bZ/tSpU8oMIwB4/PHHodPpsGjRIpSVlWHs2LHYunUrQkNDW3BIRBQIMvMu46vDFyBQV3lBLUmYObQDBnSI8WLLiMgbWlWHxVewDgtR4Jny2g78WqRtcP+Izm3wyQNjvNAiInI1Z36/nephISLylHOXLUPNc0Z3RlRoEIorDfgoowDnyqq93DIi8gYGLETkc6oMRugMJgDAY1N7Iyo0GOfKqvFRRgGKtXoIISBJkpdbSUSexNWaicjnFFcaAACaIBUiNZbrqnYRIQCAWpNAeXWt19pGRN7BgIWIfM4lrR4AEB+lUXpSQoPViA61BC/F1seJ6OrBgIWIfI4ckMRFamzuj4uy/H3J2gNDRFcPBixE5HPsBizWvy+xh4XoqsOAhYh8zqXKuiGh+uS/iysZsBBdbRiwEJHPkXtY4iNDbO6Pt/awMIeF6OrDgIWIfI48Sygu6sohIUsAc4k9LERXHQYsRORzlFlCkXaGhNjDQnTVYcBCRD5HSbpt0MMiByycJUR0tWHAQkQ+R06qtTtLiENCRFcdBixE5FPql+W3N0uoRGcpz09EVw8GLETkU+SE29BgFSJC1DaPtYtkeX6iqxUDFiLyKZfqFY27coFDTVBdeX4OCxFdXRiwEJFPuWQnf0UmDwux2i3R1YUBCxH5lGJt41VuZZwpRHR1YsBCRD7F3jpCsroFENnDQnQ1YcBCRD5FWUfoirL8MpbnJ7o6MWAhIp/S3JAQF0AkujoxYCEinyLnpthNuo1k0i3R1YgBCxH5FGWWkL2k2yjLUBGHhIiuLgxYiMinFNtZ+FCmzBKq5CwhoqsJAxYi8hk6vRFV1rL89npY6q/YbDazPD/R1YIBCxH5DLl3pbGy/LJ2EZaAxWhmeX6iqwkDFiLyGfVnCF1Zll8WEqRCTFiwzfZEFPgYsBCRz7hU2fQMIVk8i8cRXXUYsBCRz2iuyq0szlpUjlObia4eDFiIyGcoVW7tJNzKuJ4Q0dWHAQsR+QxHe1g4JER09XE6YNm5cyduuukmJCcnQ5IkbN682eZxSZIavb388st29/nss8822L5Pnz5OHwwR+be6GiyNryMki+N6QkRXHacDFp1Oh8GDB2PFihWNPn7hwgWb25o1ayBJEm699dYm99u/f3+b5+3atcvZphGRn3N0SIgLIBJdfYKcfcL06dMxffp0u48nJiba/P3ZZ59h0qRJ6NatW9MNCQpq8Fwiuro0t46QjENCRFcft+awFBUV4auvvsKCBQua3fb06dNITk5Gt27dcPfddyM/P9/utnq9HhUVFTY3IvJ/js8SYg8L0dXGrQHLu+++i6ioKNxyyy1Nbpeamop169Zh69atWLlyJXJycjBu3DhUVlY2uv3y5csRExOj3FJSUtzRfCLyoPpl+ZudJWRdALFEa2B5fqKrhFsDljVr1uDuu+9GaGhok9tNnz4ds2bNwqBBgzB16lRs2bIFZWVl+PjjjxvdftmyZSgvL1duBQUF7mg+EXmQ3FsSFqxGhKbp0er65fnLWJ6f6KrgdA6Lo3766SecOnUKH330kdPPjY2NRa9evZCVldXo4xqNBhpN01dgRORflOGgqKZnCAGW8vyx4cEoq6pFsVaPthHNP4eI/JvbelhWr16N4cOHY/DgwU4/V6vVIjs7G0lJSW5oGRH5ImWGUDP5KzIlj4WJt0RXBacDFq1Wi4MHD+LgwYMAgJycHBw8eNAmSbaiogIbN27Efffd1+g+Jk+ejDfeeEP5+9FHH8WOHTuQm5uL3bt34+abb4Zarcbs2bOdbR4R+alLDs4QksmBDcvzE10dnB4SysjIwKRJk5S/ly5dCgCYO3cu1q1bBwDYsGEDhBB2A47s7GwUFxcrf589exazZ89GSUkJ4uPjMXbsWOzZswfx8fHONo+I/JTcUxLXTMKtLI5Tm4muKk4HLBMnToQQTWflL1q0CIsWLbL7eG5urs3fGzZscLYZRBRgLjk4pVkmL4DI9YSIrg5cS4iIfEKxg1VuZSweR3R1YcBCRD7B0XWEZCweR3R1YcBCRD7B2SEhridEdHVhwEJEPqG40pKLwiEhImoMAxYi8jqd3ojqWktZfseTbi3blehYnp/oasCAhYi8Tu4lcaQsv6ydNdfFxPL8RFcFBixE5HVKwq2Dw0EAEKxWoU14MAAOCxFdDRiwEJHXKesIOThDSMaZQkRXDwYsROR1cvG3dg7mr8jaKcXjGLAQBToGLETkdVq9EQAQHRrs1PPk7eXnE1HgYsBCRF6nrbEEHFGhzq0WEmndXn4+EQUuBixE5HVyD0mkgzOEZFHW7dnDQhT4GLAQkddVWntIIlvYw1LJHhaigMeAhYi8rrLGUkfF2R6WSE2w9fkMWIgCHQMWIvI6eUjH2RwWeXutnoXjiAIdAxYi8rrWByzsYSEKdAxYiMjr5Fk+8hCPo+QhJM4SIgp8DFiIyOsqWzhLSN6+kj0sRAGPAQsReR3rsBBRcxiwEJFXGU1mVNeaALSkDgsr3RJdLRiwEJFX6fQm5d8Rzg4JWXtYqgwmmMzCpe0iIt/CgIWIvKrSOiVZE6RCSJBzp6QIjVr5N3tZiAIbAxYi8qqWTmkGAE2QWglyGLAQBTYGLETkVXVTmp0PWIB66wkx8ZYooDFgISKvauk6QrK69YRY7ZYokDFgISKvkmuoRDlZNE4mDyWxFgtRYGPAQkRepW1tDwuHhIiuCgxYiMir5IULo1qYwxLJWixEVwUGLETkVa3tYYlitVuiqwIDFiLyqpauIyTjekJEVwenA5adO3fipptuQnJyMiRJwubNm20enzdvHiRJsrlNmzat2f2uWLECXbp0QWhoKFJTU7Fv3z5nm0ZEfqjVOSzsYSG6KjgdsOh0OgwePBgrVqywu820adNw4cIF5fbhhx82uc+PPvoIS5cuxTPPPIP9+/dj8ODBmDp1Ki5evOhs84jIzyiF41rZwyLnwhBRYHL6DDF9+nRMnz69yW00Gg0SExMd3uerr76KhQsXYv78+QCAN998E1999RXWrFmDJ5980tkmEpEfkQOWVuewcEiIKKC5JYdl+/btaN++PXr37o0HHngAJSUldrc1GAzIzMxEWlpaXaNUKqSlpSE9Pb3R5+j1elRUVNjciMg/KYXjWliHRclh4ZAQUUBzecAybdo0vPfee9i2bRteeukl7NixA9OnT4fJZGp0++LiYphMJiQkJNjcn5CQgMLCwkafs3z5csTExCi3lJQUVx8GEXmI1kVJt+xhIQpsLTtDNOHOO+9U/j1w4EAMGjQI3bt3x/bt2zF58mSXvMayZcuwdOlS5e+KigoGLUR+Si6p35LFDy3PC7buhwELUSBz+7Tmbt26IS4uDllZWY0+HhcXB7VajaKiIpv7i4qK7ObBaDQaREdH29yIyD/Js3taHrBwlhDR1cDtAcvZs2dRUlKCpKSkRh8PCQnB8OHDsW3bNuU+s9mMbdu2YfTo0e5uHhF5kcksoDNYhos5JERETXE6YNFqtTh48CAOHjwIAMjJycHBgweRn58PrVaLxx57DHv27EFubi62bduGGTNmoEePHpg6daqyj8mTJ+ONN95Q/l66dClWrVqFd999FydOnMADDzwAnU6nzBoiosCkM9QFGa2uw6I3wmwWLmkXEfkep88QGRkZmDRpkvK3nEsyd+5crFy5EocPH8a7776LsrIyJCcnY8qUKXj++eeh0WiU52RnZ6O4uFj5+4477sClS5fw9NNPo7CwEEOGDMHWrVsbJOISUWCRh3FC1CpogtQt2kf9nhmdwajktBBRYHE6YJk4cSKEsH8V88033zS7j9zc3Ab3LVmyBEuWLHG2OUTkx1pbgwUANEEqBKsl1JoEtHoGLESBimsJEZHX1NVgaXnAIklSXR4LE2+JAhYDFiLymtbWYJHJPTRcAJEocDFgISKvae3ChzK5Si57WIgCFwMWIvIaecHCli58KIvi1GaigMeAhYi8prKVReNkLB5HFPgYsBCR11S6akjI+vwKa5l/Igo8DFiIyGvqkm5bNxWZ1W6JAh8DFiLymtauIySL5JAQUcBjwEJEXuOqac1MuiUKfAxYiMhrKl1Vh0XDOixEgY4BCxF5jdaaJNv6pFvWYSEKdAxYiMhr5CGc1tZhYdItUeBjwEJEXuOqSresw0IU+BiwEJHXuDqHhT0sRIGLAQsReYUQom6WkIumNVeycBxRwGLAQkReUWUwQQjLv6NDW1c4ThkS0hsh5J0SUUBhwEJEXiGX5Q9SSdAEte5UFGWtlGsWlkCIiAIPAxYi8gp5pebI0CBIktSqfYUGq6BWSdb9Mo+FKBAxYCEir1AWPmxlwi0ASJJUVzyOM4WIAhIDFiLyCleV5ZdxphBRYGPAQkRe4aqFD2WsxUIU2BiwEJFXuKoGi6yuh4VTm4kCEQMWIvKKuiq3rZvSLKurxcIeFqJAxICFiLyCOSxE5AwGLETkFcrCh8xhISIHMGAhIq+Qh25au1KzLMo6tMQeFqLAxICFiLzCVesIyZQ6LAxYiAISAxYi8gqtdaFCl+ewcEiIKCAxYCEir6h0cR0WrthMFNgYsBCRV9TNEnLNtOYozhIiCmgMWIjIK5S1hFzew8KAhSgQOR2w7Ny5EzfddBOSk5MhSRI2b96sPFZbW4snnngCAwcOREREBJKTkzFnzhycP3++yX0+++yzkCTJ5tanTx+nD4aI/AfrsBCRM5wOWHQ6HQYPHowVK1Y0eKyqqgr79+/HU089hf3792PTpk04deoUfvvb3za73/79++PChQvKbdeuXc42jYj8hBDCfXVYGLAQBSSnzxTTp0/H9OnTG30sJiYG3333nc19b7zxBkaOHIn8/Hx06tTJfkOCgpCYmOhsc4jID9XUmmEyCwCu7GGx1mGpMUIIAUmSXLJfIvINbs9hKS8vhyRJiI2NbXK706dPIzk5Gd26dcPdd9+N/Px8u9vq9XpUVFTY3IjIf1RaFyiUJCA8RO2Sfco5LEazgN5odsk+ich3uDVgqampwRNPPIHZs2cjOjra7napqalYt24dtm7dipUrVyInJwfjxo1DZWVlo9svX74cMTExyi0lJcVdh0BEbqAsfKgJcllPSHiwGvKumHhLFHjcFrDU1tbi9ttvhxACK1eubHLb6dOnY9asWRg0aBCmTp2KLVu2oKysDB9//HGj2y9btgzl5eXKraCgwB2HQERuIueZRLtopWYAUKkkJt4SBTDXDB5fQQ5W8vLy8MMPPzTZu9KY2NhY9OrVC1lZWY0+rtFooNFoXNFUIvKC+j0srhSlCUJljZHVbokCkMt7WORg5fTp0/j+++/Rrl07p/eh1WqRnZ2NpKQkVzePiHxApYvXEZIptVj0rHZLFGicDli0Wi0OHjyIgwcPAgBycnJw8OBB5Ofno7a2FrfddhsyMjLwwQcfwGQyobCwEIWFhTAYDMo+Jk+ejDfeeEP5+9FHH8WOHTuQm5uL3bt34+abb4Zarcbs2bNbf4RE5HMq3dTDoiyAyB4WooDj9NkiIyMDkyZNUv5eunQpAGDu3Ll49tln8fnnnwMAhgwZYvO8H3/8ERMnTgQAZGdno7i4WHns7NmzmD17NkpKShAfH4+xY8diz549iI+Pd7Z5ROQHlIUPXd7DUje1mYgCi9Nni4kTJ0IIYffxph6T5ebm2vy9YcMGZ5tBRH5MKRrnhhyW+vsnosDBtYSIyOMqXVyWX8ZZQkSBiwELEXmc1sULH8q4ACJR4GLAQkQe5+qFD2V1PSycJUQUaBiwEJHHyT0srlr4UKYsgMgeFqKAw4CFiDyuUlmp2XWVbi37Yw4LUaBiwEJEHueuSrfyis3MYSEKPAxYiMjjtG6udMseFqLAw4CFiDzOXXVYOK2ZKHAxYCEij5OHhCJYmp+IHMSAhYg8Sm80wWAyA3B9wBKhUQNgDwtRIGLAQkQepdOblH+7Ouk2ypp0azCaUWsNiogoMDBgISKP0ll7P8KC1VCrJJfuW+5hqf86RBQYGLAQkUfJwzWuHg4CgCC1Cpoglc3rEFFgYMBCRB6lU8ryq5vZsmXkYab6Q09E5P8YsBCRR1W6sYel/n65nhBRYGHAQkQepfNYwMIeFqJAwoCFiDxK56aicbIoZUiIOSxEgYQBCxF5lNzz4b4eFtZiIQpEDFiIyKM8NSTEHhaiwMKAhYg8ynOzhBiwEAUSBixE5FHurMNSf79MuiUKLAxYiMijtEoPC6c1E5HjGLAQkUe5O4dFHmpi4TiiwMKAhYg8yt09LJHWBRA5S4gosDBgISKPkns+3DckJPewMGAhCiQMWIjIo9w/JMRZQkSBiAELEXlU3Swh90xrrku6ZcBCFEgYsBCRR+ncnsPC1ZqJAhEDFiLyGLNZQGdwd2l+9rAQBSIGLETkMTpDXRDh9qRbgxFCCLe8BhF5ntMBy86dO3HTTTchOTkZkiRh8+bNNo8LIfD0008jKSkJYWFhSEtLw+nTp5vd74oVK9ClSxeEhoYiNTUV+/btc7ZpROTj5GGaIJUETZB7rpeirNOahQCqDBwWIgoUTp8xdDodBg8ejBUrVjT6+N/+9jf861//wptvvom9e/ciIiICU6dORU1Njd19fvTRR1i6dCmeeeYZ7N+/H4MHD8bUqVNx8eJFZ5tHRD6sfll+SZLc8hqhwSqorLvmTCGiwOF0wDJ9+nS88MILuPnmmxs8JoTAP/7xD/zlL3/BjBkzMGjQILz33ns4f/58g56Y+l599VUsXLgQ8+fPR79+/fDmm28iPDwca9ascbZ5ROTD3J1wCwCSJDGPhSgAubRPNicnB4WFhUhLS1Pui4mJQWpqKtLT0xt9jsFgQGZmps1zVCoV0tLS7D5Hr9ejoqLC5kZEvk/n5inNMs4UIgo8Lg1YCgsLAQAJCQk29yckJCiPXam4uBgmk8mp5yxfvhwxMTHKLSUlxQWtJyJ3q3Rz0TiZvP9KLoBIFDD8cpbQsmXLUF5ertwKCgq83SQicoAnhoSAuoCFPSxEgcOlAUtiYiIAoKioyOb+oqIi5bErxcXFQa1WO/UcjUaD6OhomxsR+T5lSCjEvQFLJNcTIgo4Lg1YunbtisTERGzbtk25r6KiAnv37sXo0aMbfU5ISAiGDx9u8xyz2Yxt27bZfQ4R+SetvPBhqLsDFibdEgUap88aWq0WWVlZyt85OTk4ePAg2rZti06dOuHhhx/GCy+8gJ49e6Jr16546qmnkJycjJkzZyrPmTx5Mm6++WYsWbIEALB06VLMnTsXI0aMwMiRI/GPf/wDOp0O8+fPb/0REpHP8PyQEAMWokDh9FkjIyMDkyZNUv5eunQpAGDu3LlYt24dHn/8ceh0OixatAhlZWUYO3Ystm7ditDQUOU52dnZKC4uVv6+4447cOnSJTz99NMoLCzEkCFDsHXr1gaJuETk39y98KGMKzYTBR6nA5aJEyc2We5akiQ899xzeO655+xuk5ub2+C+JUuWKD0uRBSYdB6eJaRl0i1RwPDLWUJE5J+0HhoSqsth4bRmokDBgIWIPEbroVlCESHyLCH2sBAFCgYsROQxStKtu2cJhVoWQOQsIaLAwYCFiDxG7vFw/5AQ67AQBRoGLETkMVqPJ90yYCEKFAxYiMhjdAY56da905qVOiwGBixEgYIBCxF5hBAC2hrP9LAos4RqGLAQBQoGLETkEXqjGUazpYaTp4aEOEuIKHAwYCEij6ifAOv2xQ+t+zeYzDAYzW59LSLyDAYsROQRcm9HeIgaapXk1teqX/qfM4WIAgMDFiLyCE/NEAKAILUKocEqm9clIv/GgIWIPKJuhpD7A5b6r8OZQkSBgQELEXmEp1ZqlkVwxWaigMKAhYg8QpnS7OaEW5n8OpWc2kwUEBiwEJFH6Dy0UrMsklObiQIKAxYi8ghPJt1aXofrCREFEgYsROQRysKHbl6pWcYVm4kCCwMWIvIIz88SYg8LUSBhwEJEHqEMCXk46VbLac1EAYEBCxF5hI7TmomoFRiwEJFHyNOaPT1LiCs2EwUGBixE5BGenyVkDVg4rZkoIDBgISKPUJJuPTZLiENCRIGEAQsReYQyrdnTs4SYdEsUEBiwEJFHeG2WEHtYiAICAxYi8ghPl+bnLCGiwMKAhYjczmQWqDJYhoQ8Na2Zs4SIAgsDFiJyu/p5JJ6eJaQzmGA2C4+8JhG5DwMWInI7eVgmSCVBE+SZ0079oaeqWk5tJvJ3DFiIyO2U/JXQIEiS5JHXDA1WQa2SbF6fiPyXywOWLl26QJKkBrfFixc3uv26desabBsaGurqZhGRF8nF2zw1QwgAJElCRIja+voMWIj8ncvPHr/88gtMprru16NHj+K6667DrFmz7D4nOjoap06dUv721BUYEXmGp2cIySI1QaioMbKHhSgAuPzsER8fb/P3iy++iO7du2PChAl2nyNJEhITE13dFCLyEVoPL3woqyvPz4CFyN+5NYfFYDDg/fffx7333ttkr4lWq0Xnzp2RkpKCGTNm4NixY03uV6/Xo6KiwuZGRL5LnlrsqRlCsghObSYKGG4NWDZv3oyysjLMmzfP7ja9e/fGmjVr8Nlnn+H999+H2WzGmDFjcPbsWbvPWb58OWJiYpRbSkqKG1pPRK6irCPkhSGh+q9PRP7LrQHL6tWrMX36dCQnJ9vdZvTo0ZgzZw6GDBmCCRMmYNOmTYiPj8dbb71l9znLli1DeXm5cisoKHBH84nIRTy9UrNMHoLiis1E/s9tZ4+8vDx8//332LRpk1PPCw4OxtChQ5GVlWV3G41GA41G09omEpGHeC/pNtjm9YnIf7mth2Xt2rVo3749brjhBqeeZzKZcOTIESQlJbmpZUTkaZ5eqVmmrNjMgIXI77klYDGbzVi7di3mzp2LoCDbE9ScOXOwbNky5e/nnnsO3377Lc6cOYP9+/fjd7/7HfLy8nDfffe5o2lE5AXeGxLiLCGiQOGWs8f333+P/Px83HvvvQ0ey8/Ph0pVFyddvnwZCxcuRGFhIdq0aYPhw4dj9+7d6NevnzuaRkReUDck5J1pzexhIfJ/bglYpkyZAiEaX2xs+/btNn+/9tpreO2119zRDCLyEd7qYYlkDwtRwOBaQkTkdt4fEuIsISJ/x4CFiNxOHpKJ8lYdFvawEPk9BixE5HbyLCFvDQkxYCHyfwxYiMjtvF84jgELkb9jwEJEbiWE8OpqzQB7WIgCAQMWInIrvdEMo9kya5CrNRNRSzFgISK3qh8sRIR4Z5ZQrUlAb+RMISJ/xoCFiNxKHo4JD1FDpZI8+toRIXU9OjpObSbyawxYiMittF7KXwGAILUKYcFcT4goEDBgISK38tbChzLmsRAFBgYsRORWOi9NaZZxxWaiwMCAhYjcqq4Gi2dnCMnYw0IUGBiwEJFbeTOHBWDAQhQoGLAQkVt5f0iIxeOIAgEDFiJyK2+V5ZdFcsVmooDAgIWI3Kqi2hKwRIcGe+X1o0ItAUtlTa1XXp+IXIMBCxG5VYU1UIgO804PS3SYJVAqr2bAQuTPGLAQkVvJgUJMmHd6WGIYsBAFBAYsRORWvhKwVDBgIfJrDFiIyK0qfCZg4SwhIn/GgIWI3EoOWLyVdCu/LoeEiPwbAxYicitfGRJiwELk3xiwEJHb1JrM0Bks9U8YsBBRazBgISK3qaypyxuR66F4mhywVNeaYDCavdIGImo9BixE5DZyr0akJghBau+cbiLrBUoVLB5H5LcYsBCR23g7fwUA1CpJ6d3hsBCR/2LAQkRuIwcI0V4MWADmsRAFAgYsROQ2dVOavZO/IpOnNrN4HJH/YsBCRG7jC0NC9V+fPSxE/osBCxG5ja8FLOxhIfJfDFiIyG28XZZfpgQsNSzPT+SvXB6wPPvss5AkyebWp0+fJp+zceNG9OnTB6GhoRg4cCC2bNni6mYRkRfI04i9nXQbHcZZQkT+zi09LP3798eFCxeU265du+xuu3v3bsyePRsLFizAgQMHMHPmTMycORNHjx51R9OIyIN8bUiovIoBC5G/ckvAEhQUhMTEROUWFxdnd9t//vOfmDZtGh577DH07dsXzz//PIYNG4Y33njDHU0jIg/yuYCFPSxEfsstAcvp06eRnJyMbt264e6770Z+fr7dbdPT05GWlmZz39SpU5Genm73OXq9HhUVFTY3IvI9vhKwRCs5LAxYiPyVywOW1NRUrFu3Dlu3bsXKlSuRk5ODcePGobKystHtCwsLkZCQYHNfQkICCgsL7b7G8uXLERMTo9xSUlJcegxE5BoV1ZYkVzmHxFui2cNC5PdcHrBMnz4ds2bNwqBBgzB16lRs2bIFZWVl+Pjjj132GsuWLUN5eblyKygocNm+ich1fKWHhUNCRP7P7Zc9sbGx6NWrF7Kyshp9PDExEUVFRTb3FRUVITEx0e4+NRoNNBqNS9tJRK5lNgufmSXEgIXI/7m9DotWq0V2djaSkpIafXz06NHYtm2bzX3fffcdRo8e7e6mEZEbVeqNEMLyb7k0vrfIr6/VG2E2C6+2hYhaxuUBy6OPPoodO3YgNzcXu3fvxs033wy1Wo3Zs2cDAObMmYNly5Yp2z/00EPYunUr/v73v+PkyZN49tlnkZGRgSVLlri6aUTkQXLROE2QCqHBaq+2Re5hEQKoZPE4Ir/k8iGhs2fPYvbs2SgpKUF8fDzGjh2LPXv2ID4+HgCQn58PlaouThozZgzWr1+Pv/zlL/jTn/6Enj17YvPmzRgwYICrm0ZEHuQr+SsAEBKkQliwGtW1JpRX1yIm3PttIiLnuDxg2bBhQ5OPb9++vcF9s2bNwqxZs1zdFCLyIl8pyy+LCQtWAhYi8j9cS4iI3EIODLydcCuTp1azFguRf2LAQkRuIQcGvtTDAnCmEJG/YsBCRG7hSzksAAMWIn/HgIWI3MLXAhalPD8DFiK/xICFiNxCyWEJ9W5Zfplci4U9LET+iQELEblF3TpCvtHDwiEhIv/GgIWI3MLXhoQYsBD5NwYsROQWvjet2ZrDwkq3RH6JAQsRuYUvFo4D2MNC5K8YsBCRW/hqHRbOEiLyTwxYiMjlhBDMYSEil2LAQkQuV11rQq1JAPClHBZraf7qWgghvNwaInIWAxYicjm5F0OtkhARovZyayzkHhajWaDKYPJya4jIWQxYiMjl5BosMWHBkCTJy62xCAtWI1htaQuHhYj8DwMWInI5X8tfAQBJkuoSb7liM5HfYcBCRC7na2X5ZUp5/ioGLET+hgELEbmcrxWNk0VzphCR32LAQkQu52tF42Sc2kzkvxiwEJHL+WIOC1CveBzL8xP5HQYsRORyvjskZMmpYQ8Lkf9hwEJELufrQ0Isz0/kfxiwEJHL+do6QjLmsBD5LwYsRORyddOafStgkdvDHhYi/8OAhYhczteTbtnDQuR/GLAQkcsxYCEiV2PAQkQuV38tIV/CwnFE/osBCxG5lMFoRnWtZTVkeRqxr+BaQkT+iwELEblU/d6LKF9LurUGLDW1ZuiNJi+3hoicwYCFiFxKDliiQoOgVklebo2tKE0QJGuTOCxE5F8YsBCRS/lqDRYAUKmkelObWZ6fyJ+4PGBZvnw5rrnmGkRFRaF9+/aYOXMmTp061eRz1q1bB0mSbG6hoaGubhoReYCv1mCRsTw/kX9yecCyY8cOLF68GHv27MF3332H2tpaTJkyBTqdrsnnRUdH48KFC8otLy/P1U0jIg/w1bL8MpbnJ/JPLk/h37p1q83f69atQ/v27ZGZmYnx48fbfZ4kSUhMTHR1c4jIw3y1BouMtViI/JPbc1jKy8sBAG3btm1yO61Wi86dOyMlJQUzZszAsWPH7G6r1+tRUVFhcyMi31ChrNTsW1OaZUoOC6c2E/kVtwYsZrMZDz/8MK699loMGDDA7na9e/fGmjVr8Nlnn+H999+H2WzGmDFjcPbs2Ua3X758OWJiYpRbSkqKuw6BiJzkNz0sVQxYiPyJWwOWxYsX4+jRo9iwYUOT240ePRpz5szBkCFDMGHCBGzatAnx8fF46623Gt1+2bJlKC8vV24FBQXuaD4RtYDfBCwcEiLyK27rs12yZAm+/PJL7Ny5Ex07dnTqucHBwRg6dCiysrIafVyj0UCj0biimUTkYr4esLA8P5F/cnkPixACS5YswaeffooffvgBXbt2dXofJpMJR44cQVJSkqubR0RuJtc3ifbxgIU5LET+xeU9LIsXL8b69evx2WefISoqCoWFhQCAmJgYhIWFAQDmzJmDDh06YPny5QCA5557DqNGjUKPHj1QVlaGl19+GXl5ebjvvvtc3TwicjOlDouPBiwcEiLyTy4PWFauXAkAmDhxos39a9euxbx58wAA+fn5UKnqOncuX76MhQsXorCwEG3atMHw4cOxe/du9OvXz9XNIyI38/UhobqAhZVuifyJywMWIUSz22zfvt3m79deew2vvfaaq5tCRF7gy6X5ARaOI/JXXEuIiFzGZBaorLHmsPhqaf5Qy3UaAxYi/8KAhYhcpn5eiK/2sMSGhwAAKvVGGIxmL7eGiBzFgIWIXOZ8WTUAIC4yBCFBvnl6aRMeDI21bYXlNV5uDRE5yjfPKETkl85etgQsHWLDvNwS+yRJUtp3tqzKy60hIkcxYCEilzln7WHp2Cbcyy1pWoc21oDFGmARke9jwEJELnP2sqXHQg4IfFVHa/vOMWAh8hsMWIjIZc75wZAQUNc+uUeIiHwfAxYichk5APD5gIU9LER+hwELEbmMErD4+JBQh1hLjg17WIj8BwMWInIJnd6IsipLHRZfD1jkHJYL5dUwmZuvzk1E3seAhYhcQu6tiA4N8tkqt7KE6FAEqSTUmgQuVrIWC5E/YMBCRC5RN0PIt6c0A4BaJSExJhQA81iI/AUDFiJyCX+ZISTjTCEi/8KAhYhc4qxSNM5PAhYWjyPyKwxYiMgl/K2HpSN7WIj8CgMWInKJc37WwyIvH8AeFiL/wICFiFxC6WHxk4ClrngcF0Ak8gcMWIio1WpqTbhYqQfgP0NC9ZNuhWAtFiJfx4CFiFrtQrmllklosAptI0K83BrHJMVapjXX1JpRqjN4uTVE1BwGLETUavUTbiVJ8nJrHKMJUqN9lAYAE2+J/AEDFiJqtXNlljyQjn5QNK6+jpzaTOQ3GLA0ocpgxOYD5/DKN6e83ZSrQmF5DY6fr/B2M6gF/C3hViZX5WW1W/905Gw5Lllzp8i9/m/LCXx5+Dxqak1ea0OQ117ZD1RUG/HwRwchScDcMV0Qb+0+9jVCCLd0wxtNZryzKweJ0aH47eBkqFTu6+ov1upx/b9+QqnOgD9d3weLxnd322vVJ4TA+fIaHMi/jEMFZegSF4G7Uzu7/XUra2qRkXsZAgIqSYJKkqBWSeibFO2RHJB3d+fix1MXcf3AJNw4KAnhIa07FchF4/wl4Vbmymq3lTW1+PzQeXxzrAg3DkzC7dektHqfzblUqcepwkqYhYBJCAhh+Tyldm2HsBC1W19bCIE1P+fiYkUNhqTEYminNspyB+4mhMBr3/2Kf/2QhcToUHzzyHjEhLlv/SqjyYxNB86horoW917b1S3nQnedx10ht1iHt3eeQZBKQuZT1yE02L2fLXsYsDQhMSYUgzrG4PDZcmw7UYQ7R3Zy6f7zSnSIj9K06sfCaDLjjrf3oESrx7O/7Y+Jvdu7rH2bDpzDi1+fBACs3pWDp27sh5Fd27ps//X99YvjSuLj/205ibKqWjw2tbfbvsAFpVV48euTyMgrRVGF7RVat7hIjO7ezi2vC1hOTA+8vx+7soobPNapbTi+fWS8W08IJwsr8NcvjsEsgO2nLuH5L45jxtBk3HlNJwzoENOifco9FP5Sg0XW2mq3QggcLCjDhn0F+OLweVQZLFefe86UYFS3dujUzn1DZBU1tbjhXz8ps7Pqmz4gESt/N9xtrw0A3x4vwvNfHre5LzkmFNd0bYs/Xd8XCdHuCV7MZoFnvziG99LzAACFFTV4aetJ/N/NA93yej+dvoTnvzyOX4u0ACxB7vSBSS7b/9ajhXj+y+PoFh+B9+4d2apzXmVNLcqqapHS1rWfu++OFwEAUru1dWtg2BwOCTXjur4JAOreMFc5XVSJSa9sxw3/2tWqLs3MvMvIzLuM3JIqzFv7C5ZtOozKmlqXtPH9PZYTgiQBR86V4/a30rH4g/0oKHVt3YptJ4rwxaHzUEnA3amWoPDf27Px581HYTK7frqpEAIPf3QQXx25gKIKPdQqCQM7xGBISiwA4Lkvj7vldWXbTlzErqxihKhVGNQxBgM6RKNvUjQiNUHIL63C6l05bnttIQSe//I4zAIY0CEanduFo1JvxPt78nHj67uwfMuJFu33rJ9VuZW1ttrtnz49ipv/vRsfZRSgymBCt/gI9EmMgsFoxv+18P/SUf/+MRsXK/WIDg1C36RoDOgQjUEdY6BWSfj6aCH2nClx22vrjSbl+EZ0boO+SdFQScD58hp8dvA8nvjvYbe8bq3JjEc+Poj30vMgScDvRlnOF+v35mOvi483+5IWC9b9gntW78OvRVrIccT7e/Ncsv+yKgMe2nAA//N+Js6VVeOn08U4caGyxfs7V1aNqa/txOS/71AWInUV+fdP/j30FgYszbiuv+UN2pVVjCqD0WX73fHrJZgFkFOsw7y1+1DRwiBD/iAlWbtiP9xXgGn/+Ak/N3L17oxDBWU4fLYcIWoVvnl4PO5K7QSVBHx15ALSXt2Bo+fKW7V/mVZvxF82HwUA3DeuG/735oFYfstASJLlJPSHDQdgMJpd8lqy9OwSZOZdRkiQCh/cl4qjz07FFw+OxZp51yA6NAgnLlTgo18KXPqaslpT3Q/ZfeO64vMlY/Hlg+Pw9UPj8NyM/gCAlduzUax1z7j8d8eL8HNWCUKCVFh593D8+MeJ+OC+VNwwyHLFuHpXDi5W1Di1T6PJjELrc/wvh6XlxeNyinX4cF8+JAm4ZWgHfHz/aGxbOgH/vHMoVBKw9Vghdme37ntoT0FpFdb8bAls/377EHz90Dh8+eA4fL5kLGaPtAxFvfDVcZjdFHiv/TkXeSVVaB+lwbp7R+Lrh8bhyLNTsWbeCKgkS8/d4bNlLn3NaoMJ9/8nE58dPI8glYR/3DEEL8wcqBzvsk1HXJZf8UtuKab9Yye2nbyIIJWE+dd2wRdLxkKSgJ+zSpB9Sduq/W87UYTrXtuJzw5aLtQSrb1R359o2YVxqc6Ae1bvxfnyGhhM5laf/6/cd0ZeKQAgrR8DFp/WOyEKHduEQW8046fTrvsQZOReVv597HwFFr6b4fSXTQiB76wf8Kdv7IcNi0YhpW0YzpVV4+539mLT/rMtbp/cu3L9wET0SojC/908EFseGofBKbHQG81Y+3Nui/dd38tbT+JCeQ06tQ3HI2m9AACzR3bCG7OHIVgt4avDFxp0O7fWv344bXmda1JwbY84Zay/bUQIHra24e/fnmpxENmUD/bk4UyxDnGRIXhgom2ezswhHTCwQwy0eiNe++5Xl7+23mjCC19ZgqWF47oipW04VCoJ1/aIw4q7hmFE5zYwmgU2OBmsFVXqYTILBKsltI/yTA6Dq8g9QhU1Rqff7w+s35GJveLx6h1DMLJrW0iShN6JUfjdKEse1HNfHIfR5NqAGwBe/uYUDEYzRndrh7S+tsPAD6f1QpQmCEfPVeDTA+dc/tqXKvV444csAMDj0/ogUmMZ0o7QBOE3fRIwY0gHAMDr1m1cZdmmw/jh5EWEBquwas4I5XWenN4X7aM0OFOsU9rVWqt2nkGtSWBk17b45pHxeOam/hjQIQa/sQ65f7Anv8X7fi89FwvezcClSj26x0fgvw+MwdLrLOedlgQsWr0R89fuw5lLOuW++r8vrfXDyYswC6BvUrTXZwEyYGmGJEm4rp9rh4WEEMjMt3ygnrmpHyI1QdibU4o/fHjAqZNb9iUt8kqqEKJWYVyveIzq1g5bHxqPW4ZavsgtDSrKqgz4/NB5AMA9o+sSUPskRuPpG/sCAL4+eqHVPU6ZeaV4z3rSX37LQJskwRsGJWHFXcMAAB9nFKC8yjXBwy+5pdhzphTBagn3T2iY2HvP6M7oHh+BEp0Br2877ZLXlJVX1eIf1n0+cl0vRIXajgWrVBL+coPl//fDffk4XdTy7uHGrNmVi/xSy1Xx7yf2aPC4/F6v35vv1OdQzl9JigmD2o2J2e4QoQlCm3DL++DMTKFqgwkbMy0XBHNGd2nw+CNpvRATFoyThZVOB4DNOVhQhs8PnYckAX++oW+DnIe4SA1+P8ny/r78zSmX9gwDwCvfnIJWb8TgjjHKuaa+xZN6QJIs58sTF1wz66+wvAZfHL4AAFgz7xpM6lMXpMWEBeO5GQMAAG/uyG71a5bqDPjx1EUAwPMzBqB7fKTymByIfpJZgGqD8705QgissQ75/m5UJ3z1h3EY2qkNJvVpD0kCDp8tR2G54z2ceqMJ//OfTBw6W4424cH40/V9AED5fXGF744XAoDyO+hNDFgcIL9RP5y82CC3QQiBv2w+gttW7nY4t6OgtBqXKvUIVkuYPbITVs0ZgZAgFb49XoQ/f3rU4TLh3x23fKlGd29nc5Wz7Pq+St7JhXLnx+Y/yTwLvdGMPolRGNapjc1jwzq1Qed24agymPDNsUKn9y3TG0144r9HIAQwa3hHXNsjrsE21/VLQN+kaOiNZvy3Fb1F9clXfbcNT0FyI/kWwWoVnrqxHwBg3e5c5BTrGmzT8tc+jbKqWvRKiMQdIxqfQZLarR2m9EuAWcClORAXK2vwhrVn6YlpfRChaZjoPW1AItpFhKCwogbfn7jo8L7lGiz+lr8iqxsWcvy78sXh8yivrkVK2zCM7xXf4PE2ESF4JK0nAEtvnasCbiEEXrD2ON4ytKPdJOn513ZBh9gwFFbUYNVO1+VEHT1Xjo8zLQHY0zf1a3S2TI/2kbjempT6xo+u6fH46JcCmMyWHo8x3RueK6YNSMS0/okwmgWe/O/hVuWgfXn4PGpNAgM6RKN3YpTNY+N7xSOlbRgqaoz4wnpR54ysi1rkWi8yn5zeV0muj4/SYKg1h87RXhaTWWDpR4ewK6sY4SFqrJ0/ErdbzytnLukcrt6cdbESM97YheVfNzzf1NSasPNXy8jCFAYs/mFkF0tmdKnOgMw828h1+6lLeH9PPjLyLuPOt/c4lLwnjwcO6BCD0GA1Rndvh9dnW8a9P8oowJYjjgUC9iLf+CgNhlsDje+d7BUymwU+2Gvp7rxndOcGV2+SJOGWoR0BAJv2t7y7+f09+ci6qEVcZAj+bO1VuJIkSUoS7gd781q93svBgjLs/PUS1CoJDzTSuyKb2Ls9JvWOR61J4H+/cs1wVF6JDu+m5wIA/nxDPwSp7X/1npzeB0EqCT+euoSfTl9yyeu/vPUUdAYTBqfE4uZGrooBS+XXO6xTcf+zJ9fhfftrDRZZS6Y2y0Omd6d2tturdPeozujZPhKXq2rxTxf11n19tBAZeZcRGqzCY1N7290uNFiNJ6dbrrbf3JGNIifzkhojhMBzXxyHEMBvBydjeGf7MwYf/I2lh2fLkQvIuti6nkKjyYwP91nOSfL5oDF/ndEfUaFBOHS2HP/NbPkFzn+t5zX5PFefWiXhrpGWXpaWJN9+az0fX9uj7iJTJueHOBqwfJJZgK+OXECwWsLb94zAkJRYxIaHoEd7S4/Qlb9VjTlzSYvZq/bi0NlyvLXjDH7JLbV5/OesYlTXmpAcE4r+ydEOtcud3BawrFixAl26dEFoaChSU1Oxb9++JrffuHEj+vTpg9DQUAwcOBBbtmxxV9OcFqRW4TfWLsj6HyajyaxEpcFqCefKqjH77T3N9mpkWD9IIzrX9V5M7Z+Ie6/tCqAuEGnKpUo9DhSUAQAm9204lVkOYr51MmDZnV2CnGIdIjVBmDmk8R82+QdvV1axU92XMiGEcgJ6KK0XYsPt1x2ZObQDIkLUyL6kw96cUrvbOULuYZg5pEOz003/cmM/BKkkfH/iInb+2vqg4cWvT6LWJDC+VzwmNHJFXl+3+EhleOZ/vzrR6hlLR8+VK8MXz9i5KpbdldpJSSzMuuhYYqG/zhCSdYi1Fo9zMGBREtKDVMoVbWPq99a9l56LM61M1NQbTUqZgUXjuzdb8+TGQUkY2ikW1bUm/P3b1he/3HKkEPtySxEarFKCIXv6JEZjSr8ECAGs+DG7Va/7w8mLKKyoQduIEEwbkGh3u4ToUCUv7MNfWpZjknVRi0MFZVCrJPx2SHKj29w+oiNC1CocPluOQ9ZzsKPk8/F1/Roeh9yDsTurBDp988N43x6z7GvJpJ4Y27Ou10n+XZEvjO3JK9HhrlV7ld5+wNKrW//CUP69S+uX4BM1YtwSsHz00UdYunQpnnnmGezfvx+DBw/G1KlTcfFi493Mu3fvxuzZs7FgwQIcOHAAM2fOxMyZM3H06FF3NK9F0upNb5bf0E8yz+LXIi1iw4Px5YPj0KltOPJLq3DXqr1NzrTItCZEXXmFIgcZO08XN5vd/+PJixACGNghBkkxDX8opvS3fCHSs0tQXu14d7R8ZX3LsA6NDhsAQKd24RjZpS2EQIuS+g4UlCHrohahwSrMsHNSkEVqgjDDGiDJPT8tcex8Ob4/cREqCVg8qfmidN3jI5XchNYm8h3Iv4yvjxZCJQF/vr7x3qQrPTS5J6JDg3CysLLViZNv7TwDALhpcHKDIb4rdWwTjsnW4PwDB68g5R96f6vBIqsrz+/YkO5/rL0rNw5MarbI3/he8ZjUOx5Gs8A7rZyuvn5vPvJLqxAfpcH947s1u70kSUrAtDHzLH5tRU6UEAKvWwP+ReO7NzqceqUHf2MZEvvs4DnktmJoVf7ezxrREZqgpusT3Ta8I9QqCQfyy1qUAyZPVJjYKx5xkY0XCm0XqcH1Ay3nV7mnzRFFFTU4VFAGSQLS+jW8yOweH4ku7cJhMJmbvUiqqTVhd7ZlGveVPezDrQFLZhOJtwWlVZj99h4UVtSgV0IkvnhwLMJD1DiQX6b08JvNQhkaTvPydGaZWwKWV199FQsXLsT8+fPRr18/vPnmmwgPD8eaNWsa3f6f//wnpk2bhsceewx9+/bF888/j2HDhuGNN95wR/NaZELveISoVcgp1iH7khZVBiNetc7kePA3PdE7MQrrF6aiQ2wYcop1mL1qT6P1Vcqra/GrtYt0eOcr8kM6t0GUJgilOgMONzNtuC5Sb/yD1DUuAj3aR8JoFth+yrF8hAvl1coHVE4us+eWYZYgYtP+s04P1WzMsIyBXz8gCdGhzRchustasG/r0QstrlkjBx03DkpGt3pJdE25f0I3BKkk7MstxZGzLZ/G/bY1YLh5aMcGY+L2xIaH4AFrYuzbO7NbPBxWWF6Dr49YkhUd+ZED6icWnnUoYdPvh4ScyGG5rDMouQu/G930d0QmV23etP9si3NZjCYz3vnJEvA8NLmn3YuJKw3r1AZT+1t6Ot756UyLXhuwXPicLKxEWLAaC6w9wc0Z2DEGk3rHwyyAf29vWdCfX1KFndZh0bscKNzZPioUk6wzeTY6OSxkNgvl4uCWYQ2Hg+qTe0A/P3Te4fdUnrQxJCW20dl0kiTVXRg3Myz0S24pqmtNaB+lQd8k23PKiC6WC+HD58qhNzZMDD5XVo3Zq/bgfHkNusdH4IP7RqFPYjQWjrOcH/72zUkYjGYcPFuGS5V6RGmCMKqb+wppOsPlAYvBYEBmZibS0tLqXkSlQlpaGtLT0xt9Tnp6us32ADB16lS723tDpCZIqX767fEivPNTDi5W6pHSNkwpXtSxTTg2LBqFpJhQZF/SYdmmhsWT9udfhhBAl3bhDUr9B6tVSvLpjlP2I+xqgwm7siyPNxX5Oju76cN9dYltvRKa/mG9flASQoJUOH1Ri6PnHM/KrzaY8MUhyw/orCa60+sbYC3qVmsS2Jjp/IyLrIuV+Pqo5aphyW8azo6xJyE6FDda65Os/bllV8f5JVVKcvIiBwMG2V2pnRARosavRVrsaOGw1Pt78mA0C4zs0tbhKrbje8ajU9twVNYY8fnBphMLhRB1PSyx/rXwocyZHBY5Ib1/crSSJNmcUd3aok9iFGpqzdjQwqGKr48W4lxZNdpFhOC24U3/mF5J/txtPnC+xQG/XPPltuEdERPueKXTBydbelk27T+H8y0ozvfhL/kQAhjXMw6d20U49JzbR9Tl2NU6Mdttz5kSXCivQVRoUKPD7PUN69QGfRKjoDeaHT4nfdfMRSZQl8fy48mLTc7Uk38fJvSKbzBU06VdONpFhMBgNDdaL+uRDQdx9nI1usZF4MOFo5TfoUXjuyEuUoO8kip8sDdPae+E3vEICfKNdFeXt6K4uBgmkwkJCbZvSkJCAgoLG8/NKCwsdGp7vV6PiooKm5snyB+0T/efw1s7LOOyj0/tY9NNmdI2HO/eOxIAsO3kReSV2HaF2hsOkk3sbclv2P6r/V6Rn7OKUVNrRofYsAbRdX3ymOj2U5cajbTrqzWZ8ZH1ZNpc7woARIcGK/t3ZgbP10cvQKs3olPbcKQ6UeZfTrZbvzff6WJY8pXpdf0Smg3ErnTvWMvV5BeHzztdUA2wnOjNwjI04GjviiwmLBh3XGM5bvkYnFFTa8J6a67QvGu7OPw8lUpSgvD30ptOdi7WGqA3miFJ8Ng6Mq4mDwkVaw1N1kIym4WSaHnPqIYJ6fZIkqXwGGD5/3S2LosQQukd+d2ozk4v2zC8c1sM7RQLg8mM/1gTv52RU6zDtpOW85EznyPA8sM+uls7GM0C63Y799p6owkfW6eEO7O+16Q+7REXGYJirR7bm7jwu5KcbHvjoORm/48lSbIpA9BcD2hlTS3SrUM4UxrJX5GN6NwGseHBuFxV22TS7HbrBUxjS7FIkoRhch7LFcNCR8+VY1+upazDe/eORPt6yydEaILwyHWWAPNf205ji7Vn1hemM8t8I2xy0vLlyxETE6PcUlLcv8gYUPfGnb6otcy46BijXIHX1yshChN7x0OIhmOcciLUiC6N5xJMsAYshwrKcNnOtDQ5Eeq6ZhKhBneMRfsoDbR6I/acaToBa9uJiyiq0CMuMgTT+tv/QtV3q7Xb9PND5x2+kvnYOhw0a3hHpxYQu3FQMqJDg3D2crXSReyIS5V6bLJ28zo6JFLfoI6xGNG5DWpNwqnxasBSd0U+3oXjHOtGv9L8a7tArZKwK6vY6ZWsPz90HqU6AzrEhjk9JXHW8BSEBKlw/EIF9jdR00HulUiICvWZqzBnxYQFI8JaA6ipXpadpy8hr6QKUaFBdhMy7ZkxpAPahAfjXFm108XBfsm9jENny6EJUtnURXLGfWMtn/3/7Mlzun7Iu7tzIQQwqXe8TU0SR8k9POv35jtVnO+bY0Uo0RmQEK1ptsejvmC1ShnSkb9/zakyGPH1UcsP9K3DGp9scKUZQywTAs4U65DezLIAO369BIPJjG7WoXp7gtQqpTidvc/J2ctVyLqohUoCxjZSDgKon3hr+939j3X9pWkDkhpdb+iOESnoYZ3ZlldShSCV5NL16VrL5WeYuLg4qNVqFBXZ/mcXFRUhMbHxH8LExESntl+2bBnKy8uVW0GBe8qoXykhOhSDO9Z1q1vqnTT+ozvHemL56Je6AkO1JjMOFVi66EZ0bjxgSYoJQ++EKJgF8FMj5ZWdSYRSqSSli7G5mUfylfht1h8qR4zrGYe4yBCU6gxNDmHJ8kp02HOmFJIE3Opkt3ZYiFp5jjPJt/9Jz4XBaMbQTrENcoYcJfeyvL8336lqxOv35aPKYEKfxCi7J5bmpLQNx3TrzIh3djmegyCEUAoH3jO6c5PTqBvTJiIEvx1s+VFuqndHTlT11/wVwHJF6kgei/z/MGt4itMLloYGq3GXtZdwjZMFHeUcqFuGdbSbCNqcqf0T0LFNGC5X1WLTAcd7RMur64Ju+XvgrAm94tGzfSS0eiM+2uf4uVquJHzHNZ0Q7OTnd5b1XPHDyYu4WNl8z+g3xwpRZTChU9twh88T9ScErG/mnKQMB/Vv/sIhrd5QfmM9N/Lw8LBObewOz8kXxPvzLiv7KKsy4LNDlou3uXYC3yC1Ck9Oq5sBNqpbO68udngllwcsISEhGD58OLZt26bcZzabsW3bNowePbrR54wePdpmewD47rvv7G6v0WgQHR1tc/MUuSBSWt+EJhORJvRqj05tw1FRY8Tn1g/JiQsVqK41ITo0qMkrFWVYqJFk2YNny1CstSRCObJycv08FntDKQWlVUq9D3ldDkcEqVVKeWxHhoU+sSbBjesZ79AsgyvJw0LbThQ5NKOj2mBSKukuGtetxdPypvRLQIfYMJTqDM3mdMgMRjPW7bb8wC0Y27VVUwLlZLgvDp13uJ7G3pxSnLhQgdBgFe68pmU9kPLrbj1WaHeWh7+u0nwlueS4vVWbj54rx66sYqhVEu4d26VFr/G7UZaaLftySnHsvGNJ3GcuabHtpOXHbkELAwbA8l2Vyyas/inH4WHVjdZFHXu2j2xx0K1SScpnac3POQ71xp4qrMTenFKoJLTo89szIQpDO8XCZBbY7MAsO7mm1C3DOjj1XZUTgb85Vmh3/a9akxk/WIfUHOnpHN/LMsEjt6Sq0TWL5GEu+XeiMQM6xCAkSIUSnQG5JZZz5caMs6ipNaNvUnSTQdnkvu0xqpvlt2X6QMd62z3FLX24S5cuxapVq/Duu+/ixIkTeOCBB6DT6TB//nwAwJw5c7Bs2TJl+4ceeghbt27F3//+d5w8eRLPPvssMjIysGTJEnc0r1XuHdsVK+4ahn/eOaTJ7dT18gDe3W3JA8hQ8lfaNDkcItfp2Plrw+nNn1q/WBP7tHeoJ2RM93aICFGjqEKPI3ZmHn24z/nENpk8LPTt8aImpxGazEIJWOSkOGf1aG/pqTAL4G9bm68r8UlmAcqqatGpbbgyzbslgtQqzB1juSJZ83OOQzN2vjx8HkUVesRHaZwePrjS4JRYjOzSFrUmx/MA1lmv4m8e2rHJOjdN6Z1YN7Rpr3fn1yLLCdVfa7DI5Pbbm/q7yppDcsPApBavp5IUE6b0lq1zsJdl9a4cCAFM7tO+yaEER9x+TQqiQoNwplin/IA2xWiqWzPs3lYG3TOGJiMuUoML5TX4ylpivykvWutbTemX2KKLGwBKjZyPM5qeyXiwoExZJ66xYnFNGdAhBoOtEwI+sTMrae+ZUlTWGBEXqcGQlOZ7b+pP8PjvFcU5DUYzdlt73if0sj9UowlSY5A1yT4jt9Qm/2puIwVB65MkCW/dMwIr7x6G2dc0PzPLk9wSsNxxxx145ZVX8PTTT2PIkCE4ePAgtm7dqiTW5ufn48KFug/tmDFjsH79erz99tsYPHgwPvnkE2zevBkDBgxwR/NaJVitwg2DkhyaVnj7iBRo6uUByElU8rQze0Z0aYuIEDWKtXocr7cuRmbeZeVDZ6+0+5U0QWplDPLbRoaFDEaz0uXbVBVJe/olR+O6fgkwmQWe/+qE3RPDrqxiXCivQWx4cKuSuJ6c3geSZMnPyMi1n5djqlf3YsHYrq1e4+aOEZ0QHqLGycJKJXnOHkuSpOW1543p0mztCEfcZ82B+WBPXrNFpQpKq5T3er6TSZJXkvMPNmacRckVV5CHz5bhU+vwQmNLK/gTuf3r9+Y3qMx69nIVvrT+yDo70+tK8629HJ8dOt/g//NKpTqD8iN437jWvS5g+SGUewRWOTDF+fsTRThXVo024cF2qyM7ShOkxjxr0L/qpzNNBhA/nryIH09dQrBawhPNFKhryo2DkhAarELWRa1SZPNKluq9xwBYeleaKyjZmLtHNj0hQP4upvVt7/B5SO5VeuenMzhZWPcbkJFXCp3BhLjIkGYrzw63Dgtl5l3GDmv+VXRokNIr3pSYsGBMH5jkVJ6hJ7gtS27JkiXIy8uDXq/H3r17kZqaqjy2fft2rFu3zmb7WbNm4dSpU9Dr9Th69Ciuv/56dzXNY2LDQ5TCaO/uzlMSbpsbIw0JUmGM9QQqDwvV1Jrw+CeHIISlV6N+ZcPmNDW9+bvjRSjWGhAfpcHkFhYH+vP1fRGslrDz10vKomFXkoOimUM6tOoHfECHGCVY++sXx+12bX93vBB5JVWICQvGrBb26NQXEx6sTCdd08wU5/TsEhy/UIGwYHWLgsDGpPVNQNe4CFTUGJU6Nva8vycPZmEp/+3srKgrje7WDgM7xEBvNOO99Lqk41qTGY9/chhma5l2fw9YpvZPwKTe8TCYzHjyv0dsPldrduXCZBa4tkc7h6eG2zOsUywGdYyBwVhXbt6e9/fkQW80Y0CHaKWLvrXmXdsFQSoJe3Oary20ZlcuAMv0emdnJjXm7tTOCAtW49j5CrtBv8FoxvPW5TDmX9sVXeOc6/GtLyo0WBnCt/ed+fzQeezPL0N4iBpPTGtZcHTj4CREaYKQX1qFn7Nt8w6FEMryKM5cqE0bkIjr+iWg1iTw2MbDyswyOVdwfK/4ZoMJeXmWjLzLeM/aMztrRIrNIrP+xj/T+v2IXC1VHiIIUkkY3DG22efJw0JygtW/tp1G9iUd4qM0eOpGx6qlyib1bo8glYRfi7QNiket31fXY+NsYpusS1yEMj7+wpcnGoxRf3+8CN9Y66A0VcrcUX+c0htRmiAcOVeOT+zkzsiJiveM6ux0gqQ988Z0AWCZrv7l4cZzWaoNJrxiLYM+a0TLh2OupFJJStLjqp9y7NbT+CW3VFmzaP6Yluc8yCRJUnoV3kvPVRLI3955BicLK9EmPBjP3NSv1a/jbZIk4YWbByIiRI2Mej2Z5VW1Su0UuQBca19H7vV6a8eZRutkAJbZSu9af2QWtiL/6kpJMWHKzEa5QFhjNuzLx77cUgSpJNwzqotLXrtNRIgyHPy2nR4eyxIGOsRFhjhVM8ke+Xzz3/3nlKEUWbWhbqmD30/sjoTolk3LDw8JUgpp1k++FULgH9+fxvnyGoSHqJ0K6iVJwgszByA61HKeW2XtsZV/D5pb3gOouzDOuqhVpkHf40DJCl/GgMXNBnSIwbBOsZAv2Pp3iHEowpUTqvbnl+HnrGKlvPoLMwc4/SMYEx6M+ydYfnRe+OoEnrP2TOQU6/BzVgkkCbjTiWTbxiz+TQ+0iwjBmWKdzZX49lMX8fsP9sNoFrhlaAf0c8ECWvFRGvzBWpDqb1tPofKKqZKZeaXYn1+GELUKc8a47gvaLT4Sc0Z3hhDAwxsONuixqjIYce+6X5QrNnkqqavcNqwjkmNCca6sGne+nd4gAffouXLcu/YX1NSaMbF3vLL+VWtNH5CozDD5JLMA2Ze0ymJ+T9/UD+1aOHPF13SIDVOGIF76+iTOlVXj/b15ykyv8U70ajblxkHJGNm1LSr1RsxZs6/BEFRBaRXueCsdJToDusVHKL0ErvLAxB4ICVLhp9PF+P0HmQ1qNG0+cA7LPj0CwFLt2ZX1de4d2xUqyZI4emW+ULFWr3yuHpva26Eq2M1J7doW1/VLgMFoxoJ3M7Cv3npkb+3MxoXyGnSIDWv1kNtd1jox3x4vwsWKGhhNZizbdEQ5ngd/09PpXqqE6FBlaYXXvv8VP2cV42RhJVSSpbhjc9pFatDN2kMlhCXI6dKKHitfwIDFA+ReFsD+dOYrdWwTjh7tI2EyCyx8LwMms8BNg5MxtYXJo49O6Y1l1pPxmp9z8OCGA8oV3IRe8S1OJJRFhwbjUevqsf/8/leU6gz4OasYi/6TCYPJjOkDEvG32wa16jXqmzumC7rGRaBYq1eWsK+oqcWKH7Ow6L1MAJZFGhsrgd0az9zUHzOHJMNoFlj8wX5lzQ+t3oh5a39B+pkSRISo8d69I1s0Ht6UsBA11i8chWRrJeU73kpXqoeeLqrEPav3olJvxMiubbHy7uEuG38OUqtwX73enWX/PQKD0RIU2Vsg01/9LrUzRnRuA53BhGWbjihJzovGu66XI1itwuq5IzCoYwxKdQbc/c5eFJRaZnLkFOtw+1vpSiXS9xektrjn057eiVF4Z84IaIJU+P7ERdz/n0xluv7Woxfwx42Woee7Uzvh0Sn2V4Ruic7tIpQFDH/3zl6889MZZfmHv3/7KyprjBjQIRq3DXdNbS1JkvDGXUMxoVc8qmtNmL92HzLzLuN8WTXetBb/XHZ9n1YPefVOjMLwzm1gshbIW/heBjb8UgCVZLnIlBdldNZtwztiQq94GIxmLHwvA4AlCb9NM2tYyeqnH8x14cWbt0iipYuU+JCKigrExMSgvLzco1OcHaU3mnDtiz+gWGvAm78b3uSKo/U9/+VxrLYmjraNCMF3j4xv9dXsZwfP4dGNh1BrqnvbV80Z4ZJqhiazwI2v78KJCxUY1zMOv+SWoqbWjLS+CVj5u2EuP/H+cLII967LQLBawpzRXfBxRgEqaywnvy7twvH+famtDsQaYzSZ8eCHB/D10UKEBquw4q5hWLk9Gxl5lxGlCcK7C0Y2u8hgaxSUVmH2qj04e7kaHduE4aVbB2HpxwdRVKHH4I4xeP++VES54Oq0viqDEWNe/AFl1nVTwkPU+PaR8W75//W2rItaXP/Pn2CwDm0mxYRi5+OTXP75vawz4I630/FrkRYpbcPw4i2D8MhHB3GxUo8e7SOx/r5Um0qkrvZzVjEWvGvpkRvXMw53jeyEP2w4gFqTwK3DOuLl2wa5Jeky+5IW97yzF+etK723CQ/GLcM6Yq21KvTG/xmNa5qZmOCsmloTFrz7C37OKkGUJgj9O0Rjz5lSXNOlDT6+f7RLgtH/Zp7FHzceUv4ODVbh9dnDWn1uPVdWjamv7YTWmmz/cFpPPJzWy6HnfpxRgMc/OYyUtmHY/uikVk8+cAdnfr8ZsHjIvpxS7D1Tgt9P6uHwh+an05dwz+p9AIDXZw/FTYNbNz1Wtju7GPe/l4lKvRGJ0aHY9cQkpwuL2ZOeXYLZq/Yof0/sHY+37hnukpkyVxJCYO7aX2xWNu3RPhK/n9gdvx2c7LJjaozBaMb9/8nAj/UK5kWHBuE/C1Ix2ME1ZlrjfFk17lq1R6mxAAC9E6KwYdEoh6++nPX3b0/hdesiks/e1A/zHFwEzx+9vu00/m5d3PTP1/fFwlbODrLnYkUNZr2Vjrx672OfxCi8f19qi4vEOWPPmRLcu+4XVNWrfnvjoCT8886hbv1x0xtN+HT/OazckW1z7DcOSsIbdw1zy2tWGSy9oPKwkCQBny8ei4EdW5dILaupNSH1/7ahvLoWbSNC8M7cES67cPlgbx7+/OlRAMDmxddiiIPnGIPRjDd+zMLkPu09cl5qCQYsAUKeiZEQHYonpvV2WZc0AJwsrMBLX5/E7SNSMN3FY+S//yATW44UYlzPOKyaM8IlMwzsyb6kxV2r9iAhOhS/n9gdU/olemwqXv2rtpiwYHxwX2qrZ5E442JFDe56Zy+yLmrRpV04Pr5/tFuvyIu1etz875/RIz4S78y9xiev1lzFYDTjntV7UazVY/Pia13eY1VfQWkVbn8rHRfKazCgQzT+c2+q24LOxvySW4p5a/ZBZzDhun4J+Pfdru8NtcdoMuOrIxfw5o4zqKypxUf3j3ZrTR+t3oi5ayzDQneMSMFLLhymBiwrcn91+AL+fENfh1eFd4TZLPDXL45BbzTj/24e6HPTjVuDAQt5ld5owp4zpRjVra1belZ8SZXBiE37z+HaHnGtmoLZUqU6A746fB5TByS6PF+HPOd8WTW2nSjCjKEdXJJs6qzTRZX4Jfcybh3eurID/qCm1oS9OaW4tns7t/bCkmMYsBAREZHPc+b3m+ElERER+TwGLEREROTzGLAQERGRz2PAQkRERD6PAQsRERH5PAYsRERE5PMYsBAREZHPY8BCREREPo8BCxEREfk8BixERETk8xiwEBERkc9jwEJEREQ+jwELERER+bwgbzfAFeQFpysqKrzcEiIiInKU/Lst/443JSAClsrKSgBASkqKl1tCREREzqqsrERMTEyT20jCkbDGx5nNZpw/fx5RUVGQJMml+66oqEBKSgoKCgoQHR3t0n37gkA/PiDwj5HH5/8C/RgD/fiAwD9Gdx2fEAKVlZVITk6GStV0lkpA9LCoVCp07NjRra8RHR0dkB9CWaAfHxD4x8jj83+BfoyBfnxA4B+jO46vuZ4VGZNuiYiIyOcxYCEiIiKfx4ClGRqNBs888ww0Go23m+IWgX58QOAfI4/P/wX6MQb68QGBf4y+cHwBkXRLREREgY09LEREROTzGLAQERGRz2PAQkRERD6PAQsRERH5vKs+YPnf//1fjBkzBuHh4YiNjXXoOUIIPP3000hKSkJYWBjS0tJw+vRpm21KS0tx9913Izo6GrGxsViwYAG0Wq0bjqB5zrYlNzcXkiQ1etu4caOyXWOPb9iwwROHZKMl/9cTJ05s0Pb/+Z//sdkmPz8fN9xwA8LDw9G+fXs89thjMBqN7jyURjl7fKWlpXjwwQfRu3dvhIWFoVOnTvjDH/6A8vJym+28+f6tWLECXbp0QWhoKFJTU7Fv374mt9+4cSP69OmD0NBQDBw4EFu2bLF53JHvpCc5c3yrVq3CuHHj0KZNG7Rp0wZpaWkNtp83b16D92ratGnuPowmOXOM69ata9D+0NBQm238+T1s7HwiSRJuuOEGZRtfeg937tyJm266CcnJyZAkCZs3b272Odu3b8ewYcOg0WjQo0cPrFu3rsE2zn6vnSauck8//bR49dVXxdKlS0VMTIxDz3nxxRdFTEyM2Lx5szh06JD47W9/K7p27Sqqq6uVbaZNmyYGDx4s9uzZI3766SfRo0cPMXv2bDcdRdOcbYvRaBQXLlywuf31r38VkZGRorKyUtkOgFi7dq3NdvX/DzylJf/XEyZMEAsXLrRpe3l5ufK40WgUAwYMEGlpaeLAgQNiy5YtIi4uTixbtszdh9OAs8d35MgRccstt4jPP/9cZGVliW3btomePXuKW2+91WY7b71/GzZsECEhIWLNmjXi2LFjYuHChSI2NlYUFRU1uv3PP/8s1Gq1+Nvf/iaOHz8u/vKXv4jg4GBx5MgRZRtHvpOe4uzx3XXXXWLFihXiwIED4sSJE2LevHkiJiZGnD17Vtlm7ty5Ytq0aTbvVWlpqacOqQFnj3Ht2rUiOjrapv2FhYU22/jze1hSUmJzbEePHhVqtVqsXbtW2caX3sMtW7aIP//5z2LTpk0CgPj000+b3P7MmTMiPDxcLF26VBw/fly8/vrrQq1Wi61btyrbOPt/1hJXfcAiW7t2rUMBi9lsFomJieLll19W7isrKxMajUZ8+OGHQgghjh8/LgCIX375Rdnm66+/FpIkiXPnzrm87U1xVVuGDBki7r33Xpv7HPmgu1tLj2/ChAnioYcesvv4li1bhEqlsjmprly5UkRHRwu9Xu+StjvCVe/fxx9/LEJCQkRtba1yn7fev5EjR4rFixcrf5tMJpGcnCyWL1/e6Pa33367uOGGG2zuS01NFffff78QwrHvpCc5e3xXMhqNIioqSrz77rvKfXPnzhUzZsxwdVNbzNljbO78Gmjv4WuvvSaioqKEVqtV7vO191DmyHng8ccfF/3797e574477hBTp05V/m7t/5kjrvohIWfl5OSgsLAQaWlpyn0xMTFITU1Feno6ACA9PR2xsbEYMWKEsk1aWhpUKhX27t3r0fa6oi2ZmZk4ePAgFixY0OCxxYsXIy4uDiNHjsSaNWscWiLclVpzfB988AHi4uIwYMAALFu2DFVVVTb7HThwIBISEpT7pk6dioqKChw7dsz1B2KHqz5L5eXliI6ORlCQ7fJhnn7/DAYDMjMzbb4/KpUKaWlpyvfnSunp6TbbA5b3Qt7eke+kp7Tk+K5UVVWF2tpatG3b1ub+7du3o3379ujduzceeOABlJSUuLTtjmrpMWq1WnTu3BkpKSmYMWOGzfco0N7D1atX484770RERITN/b7yHjqrue+gK/7PHBEQix96UmFhIQDY/JDJf8uPFRYWon379jaPBwUFoW3btso2nuKKtqxevRp9+/bFmDFjbO5/7rnn8Jvf/Abh4eH49ttv8fvf/x5arRZ/+MMfXNb+5rT0+O666y507twZycnJOHz4MJ544gmcOnUKmzZtUvbb2HssP+Yprnj/iouL8fzzz2PRokU293vj/SsuLobJZGr0//bkyZONPsfee1H/+ybfZ28bT2nJ8V3piSeeQHJyss3Jf9q0abjlllvQtWtXZGdn409/+hOmT5+O9PR0qNVqlx5Dc1pyjL1798aaNWswaNAglJeX45VXXsGYMWNw7NgxdOzYMaDew3379uHo0aNYvXq1zf2+9B46y953sKKiAtXV1bh8+XKrP/eOCMiA5cknn8RLL73U5DYnTpxAnz59PNQi13P0GFururoa69evx1NPPdXgsfr3DR06FDqdDi+//LJLfvDcfXz1f7wHDhyIpKQkTJ48GdnZ2ejevXuL9+soT71/FRUVuOGGG9CvXz88++yzNo+58/2jlnnxxRexYcMGbN++3SYp9c4771T+PXDgQAwaNAjdu3fH9u3bMXnyZG801SmjR4/G6NGjlb/HjBmDvn374q233sLzzz/vxZa53urVqzFw4ECMHDnS5n5/fw99QUAGLH/84x8xb968Jrfp1q1bi/admJgIACgqKkJSUpJyf1FREYYMGaJsc/HiRZvnGY1GlJaWKs9vLUePsbVt+eSTT1BVVYU5c+Y0u21qaiqef/556PX6Vq834anjk6WmpgIAsrKy0L17dyQmJjbIcC8qKgIAl7yHnji+yspKTJs2DVFRUfj0008RHBzc5PaufP/siYuLg1qtVv4vZUVFRXaPJzExscntHflOekpLjk/2yiuv4MUXX8T333+PQYMGNbltt27dEBcXh6ysLI//2LXmGGXBwcEYOnQosrKyAATOe6jT6bBhwwY899xzzb6ON99DZ9n7DkZHRyMsLAxqtbrVnwmHuCwbxs85m3T7yiuvKPeVl5c3mnSbkZGhbPPNN994Nem2pW2ZMGFCg9kl9rzwwguiTZs2LW5rS7jq/3rXrl0CgDh06JAQoi7ptn6G+1tvvSWio6NFTU2N6w6gGS09vvLycjFq1CgxYcIEodPpHHotT71/I0eOFEuWLFH+NplMokOHDk0m3d544402940ePbpB0m1T30lPcvb4hBDipZdeEtHR0SI9Pd2h1ygoKBCSJInPPvus1e1tiZYcY31Go1H07t1bPPLII0KIwHgPhbD8jmg0GlFcXNzsa3j7PZTBwaTbAQMG2Nw3e/bsBkm3rflMONRWl+3JT+Xl5YkDBw4o03YPHDggDhw4YDN9t3fv3mLTpk3K3y+++KKIjY0Vn332mTh8+LCYMWNGo9Oahw4dKvbu3St27dolevbs6dVpzU215ezZs6J3795i7969Ns87ffq0kCRJfP311w32+fnnn4tVq1aJI0eOiNOnT4t///vfIjw8XDz99NNuP54rOXt8WVlZ4rnnnhMZGRkiJydHfPbZZ6Jbt25i/PjxynPkac1TpkwRBw8eFFu3bhXx8fFem9bszPGVl5eL1NRUMXDgQJGVlWUzjdJoNAohvPv+bdiwQWg0GrFu3Tpx/PhxsWjRIhEbG6vMyLrnnnvEk08+qWz/888/i6CgIPHKK6+IEydOiGeeeabRac3NfSc9xdnje/HFF0VISIj45JNPbN4r+RxUWVkpHn30UZGeni5ycnLE999/L4YNGyZ69uzp0eC5Ncf417/+VXzzzTciOztbZGZmijvvvFOEhoaKY8eOKdv483soGzt2rLjjjjsa3O9r72FlZaXyWwdAvPrqq+LAgQMiLy9PCCHEk08+Ke655x5le3la82OPPSZOnDghVqxY0ei05qb+z1zhqg9Y5s6dKwA0uP3444/KNrDWq5CZzWbx1FNPiYSEBKHRaMTkyZPFqVOnbPZbUlIiZs+eLSIjI0V0dLSYP3++TRDkSc21JScnp8ExCyHEsmXLREpKijCZTA32+fXXX4shQ4aIyMhIERERIQYPHizefPPNRrd1N2ePLz8/X4wfP160bdtWaDQa0aNHD/HYY4/Z1GERQojc3Fwxffp0ERYWJuLi4sQf//hHm2nBnuLs8f3444+NfqYBiJycHCGE99+/119/XXTq1EmEhISIkSNHij179iiPTZgwQcydO9dm+48//lj06tVLhISEiP79+4uvvvrK5nFHvpOe5Mzxde7cudH36plnnhFCCFFVVSWmTJki4uPjRXBwsOjcubNYuHChS38IWsKZY3z44YeVbRMSEsT1118v9u/fb7M/f34PhRDi5MmTAoD49ttvG+zL195De+cI+Zjmzp0rJkyY0OA5Q4YMESEhIaJbt242v4mypv7PXEESwsPzUImIiIicxDosRERE5PMYsBAREZHPY8BCREREPo8BCxEREfk8BixERETk8xiwEBERkc9jwEJEREQ+jwELERER+TwGLEREROTzGLAQERGRz2PAQkRERD6PAQsRERH5vP8HnF5nFvvvns4AAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ],
      "source": [
        "from scipy.special import eval_chebyt\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "k = 20\n",
        "delta = 0.1\n",
        "\n",
        "def filter_func(x, delta, k):\n",
        "    val = 0\n",
        "    for i in range(k):\n",
        "        arg = -1 + 2*((x**2-delta**2)/(1-delta**2))\n",
        "        topk = eval_chebyt(k, arg)\n",
        "        arg = -1 + 2*(-delta**2/(1-delta**2))\n",
        "        botk = eval_chebyt(k, arg)\n",
        "        val += topk/botk\n",
        "    return val\n",
        "\n",
        "# Plot the polynomial\n",
        "xs = np.linspace(-1, 1, 100)\n",
        "ys = np.zeros(len(xs))\n",
        "\n",
        "for i, x in enumerate(xs):\n",
        "    ys[i] = filter_func(x, delta, k)\n",
        "\n",
        "plt.plot(xs, ys)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UUVDpk-o_lXT"
      },
      "source": [
        "This polynomial needs to be converted to phases to use Quantum Signal Processing. For convenience, here are the phases for $k=20$ and $\\Delta = 0.1$. At the end of the notebook, we show two ways to calculate the phases for a given polynomial."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "id": "ytwUrMTo_lXU"
      },
      "outputs": [],
      "source": [
        "# phases below are calculated using QSPPACK for n=20, delta = 0.1\n",
        "phases = [ 0.76974538,  0.00646461, -0.00773143,  0.00908568, -0.01051591,  0.01200849,\n",
        "          -0.01354772,  0.01511600, -0.01669402,  0.01826109, -0.01979540,  0.02127448,\n",
        "          -0.02267559,  0.02397626, -0.02515479,  0.02619086, -0.02706603,  0.02776435,\n",
        "          -0.02827282,  0.02858185, -0.02868552,  0.02858185, -0.02827282,  0.02776435,\n",
        "          -0.02706603,  0.02619086, -0.02515479,  0.02397626, -0.02267559,  0.02127448,\n",
        "          -0.01979540,  0.01826109, -0.01669402,  0.01511600, -0.01354772,  0.01200849,\n",
        "          -0.01051591,  0.00908568, -0.00773143,  0.00646461,  0.76974538]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mvtr4xEo_lXU"
      },
      "source": [
        "## Hydrogen molecule example\n",
        "\n",
        "We are going to prepare two different states using the phases above, as well as the Quantum Signal Processing circuit `get_qsp_circuit_no_anc` as defined in [arXiv:2002.11649](https://arxiv.org/abs/2002.11649).\n",
        "\n",
        "We are going to apply it to $H_2$ in a STO-3G basis with the \"scbk\" qubit mapping. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "id": "I1BkclWt_lXU"
      },
      "outputs": [],
      "source": [
        "from openfermion import get_sparse_operator\n",
        "import numpy as np\n",
        "\n",
        "from tangelo.molecule_library import mol_H2_sto3g\n",
        "from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping\n",
        "from tangelo.linq.helpers.circuits.statevector import StateVector\n",
        "\n",
        "fe_op = mol_H2_sto3g.fermionic_hamiltonian\n",
        "qu_op = fermion_to_qubit_mapping(fe_op, \"scbk\", n_spinorbitals=4, n_electrons=2)\n",
        "\n",
        "ham = get_sparse_operator(qu_op).toarray()\n",
        "eigs, vecs = np.linalg.eigh(ham)\n",
        "alpha = eigs[3]-eigs[0]\n",
        "\n",
        "init_state = np.random.random(4)\n",
        "init_state /= np.linalg.norm(init_state)\n",
        "\n",
        "sv = StateVector(init_state)\n",
        "init_circ = sv.initializing_circuit()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "D5nQUKJ2_lXV"
      },
      "source": [
        "### First excited state\n",
        "\n",
        "We need to shift the Hamiltonian such that the desired eigenvalue ($\\lambda$) is centered at zero and the spectral range is in [-1, 1]. This is done by defining\n",
        "\n",
        "$\\tilde{H} = \\frac{H-\\lambda}{\\alpha + \\left|\\lambda\\right|}$\n",
        "\n",
        "where $\\alpha$ is the approximation to the spectral range, defined in the previous cell."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "id": "W5tNf2eE_lXV"
      },
      "outputs": [],
      "source": [
        "from tangelo.linq import Circuit, Gate, get_backend\n",
        "from tangelo.toolboxes.circuits.qsp import get_qsp_circuit_no_anc\n",
        "from tangelo.toolboxes.circuits.lcu import get_uprep_uselect, get_lcu_qubit_op_info\n",
        "\n",
        "lamb = eigs[1]\n",
        "# Shift Hamiltonian so that spectrum is in [-1, 1] and desired lambda is at 0\n",
        "qu_op_tilde = (qu_op - lamb)/(alpha+abs(lamb))\n",
        "\n",
        "# Generate LCU circuit that applies the Hamiltonian\n",
        "uprep, uselect, qu_op_qs, m_qs, alpha = get_uprep_uselect(qu_op_tilde)\n",
        "cua = uprep + uselect + uprep.inverse()\n",
        "\n",
        "# Ancilla qubits\n",
        "n_ancilla = len(m_qs) + 1\n",
        "\n",
        "# Total number of qubits needed to apply algorithm\n",
        "n_qubits = len(qu_op_qs) + n_ancilla\n",
        "\n",
        "# Feed in LCU circuit to generate QSP circuit with phases for eigenvalue filtering\n",
        "eig_filt_circ = get_qsp_circuit_no_anc(cua, m_qs, phases)\n",
        "\n",
        "# Add the mesurement gates to the ancilla qubits for the application of the Linear Combination\n",
        "# of Unitaries plus the extra qubit for the QSP\n",
        "# Measure gates for LCU ancilla qubits\n",
        "full_circuit = init_circ + eig_filt_circ + Circuit([Gate(\"MEASURE\", m) for m in m_qs])\n",
        "# Measure gate for the QSP qubit where the phases are applied\n",
        "full_circuit += Circuit([Gate(\"MEASURE\", m_qs[-1]+1)])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "id": "5rKQQ-Vg_lXV",
        "outputId": "d7cbcbb0-6e40-4934-8ebe-f441fed2156e",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "initial overlap = 0.18367656363390905\n",
            "final overlap = 0.9887577317376988\n"
          ]
        }
      ],
      "source": [
        "sim = get_backend()\n",
        "\n",
        "f, sv = sim.simulate(full_circuit, desired_meas_result=\"0\"*n_ancilla, return_statevector=True)\n",
        "# Reorder statevector to \"msq_first\" if statevector in backend is \"lsq_first\"\n",
        "if sim.backend_info()[\"statevector_order\"] == \"lsq_first\":\n",
        "    sv = np.reshape(sv, [2]*n_qubits).T.flatten()\n",
        "# Pick out part of wavefunction that corresponds to \"0000\" on ancilla qubits\n",
        "sv = np.reshape(sv, (2**4, 4))[0, :]\n",
        "# Reorder statevector to correspond to openfermion ordering\n",
        "sv = np.reshape(sv, (2,2)).T.flatten()\n",
        "\n",
        "print(f'initial overlap = {abs(np.dot(vecs[:,1], init_state))}')\n",
        "print(f'final overlap = {abs(np.dot(vecs[:,1], sv))}')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rtG_tBsq_lXW"
      },
      "source": [
        "### Second excited state\n",
        "\n",
        "We can then reapply the same approach to the second excited state."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "id": "YRWXx4o7_lXW"
      },
      "outputs": [],
      "source": [
        "lamb = eigs[2]\n",
        "# Shift Hamiltonian so that spectrum is in [-1, 1] and desired lambda is at 0\n",
        "qu_op_tilde = (qu_op - lamb)/(alpha+abs(lamb))\n",
        "\n",
        "# Generate LCU circuit that applies the Hamiltonian\n",
        "uprep, uselect, qu_op_qs, m_qs, alpha = get_uprep_uselect(qu_op_tilde)\n",
        "cua = uprep + uselect + uprep.inverse()\n",
        "\n",
        "# Ancilla qubits\n",
        "n_ancilla = len(m_qs) + 1\n",
        "\n",
        "# Total number of qubits needed to apply algorithm\n",
        "n_qubits = len(qu_op_qs) + n_ancilla\n",
        "\n",
        "# Feed in LCU circuit to generate QSP circuit with phases for eigenvalue filtering\n",
        "eig_filt_circ = get_qsp_circuit_no_anc(cua, m_qs, phases)\n",
        "\n",
        "full_circuit = init_circ + eig_filt_circ + Circuit([Gate(\"MEASURE\", m) for m in m_qs])\n",
        "full_circuit += Circuit([Gate(\"MEASURE\", m_qs[-1]+1)])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "id": "Ba0LhGfn_lXX"
      },
      "outputs": [],
      "source": [
        "sim = get_backend()\n",
        "\n",
        "f, sv = sim.simulate(full_circuit, desired_meas_result=\"0\"*n_ancilla, return_statevector=True)\n",
        "# Reorder statevector to \"msq_first\" if statevector in backend is \"lsq_first\"\n",
        "if sim.backend_info()[\"statevector_order\"] == \"lsq_first\":\n",
        "    sv = np.reshape(sv, [2]*n_qubits).T.flatten()\n",
        "# Pick out part of wavefunction that corresponds to \"0000\" on ancilla qubits\n",
        "sv = np.reshape(sv, (2**4, 4))[0, :]\n",
        "# Reorder statevector to correspond to openfermion ordering\n",
        "sv = np.reshape(sv, (2,2)).T.flatten()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "id": "xOEVtX68_lXX",
        "outputId": "479076d4-f830-44df-c6e3-c907132bdd25",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "initial overlap = 0.18413991077227115\n",
            "final overlap = 0.9969619094421538\n"
          ]
        }
      ],
      "source": [
        "print(f'initial overlap = {abs(np.dot(vecs[:,2], init_state))}')\n",
        "print(f'final overlap = {abs(np.dot(vecs[:,2], sv))}')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ueW7rS0T_lXX"
      },
      "source": [
        "## Code to generate phases\n",
        "\n",
        "This section shows 2 options for converting a polynomial into phases, with either QSPPACK or pyqsp.\n",
        "\n",
        "The cells are formatted in markdown as additional dependencies are required to run the following phase calculations, but the code below can be run once they are installed."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nXKQtCYa_lXX"
      },
      "source": [
        "#### Using QSPPACK"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "a49aIXnd_lXY"
      },
      "source": [
        "```python\n",
        "from sympy.abc import x\n",
        "from sympy import Rational\n",
        "from sympy.functions.special.polynomials import chebyshevt\n",
        "\n",
        "k = 20\n",
        "# Necessary to work in exact representation for higher order polynomials\n",
        "delta = Rational(delta)\n",
        "\n",
        "def filter_func(x, delta, k):\n",
        "    val = 0\n",
        "    for i in range(k):\n",
        "        arg = -1 + 2*((x**2-delta**2)/(1-delta**2))\n",
        "        topk = chebyshevt(k, arg)\n",
        "        arg = -1 + 2*(-delta**2/(1-delta**2))\n",
        "        botk = chebyshevt(k, arg)\n",
        "        val += topk/botk\n",
        "    return val\n",
        "```"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "z3MmVAXP_lXY"
      },
      "source": [
        "We need to convert the `filter_func` polynomial to coefficients of expansion in Chebyshev polynomials to use with QSPPACK"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aYhwSmKb_lXY"
      },
      "source": [
        "```python\n",
        "# Use chebyshev interpolation ak=2/(n+1)*sum_{j=0}^n p(x_j)T_k(x_j), x_j = cos(\\frac{2*j+1}{2*n+2}\\pi)\n",
        "poly = filter_func(x, delta, k).as_poly()\n",
        "\n",
        "n = 2*k\n",
        "xjs = np.zeros(n+1)\n",
        "aks = np.zeros(n+1)\n",
        "for j in range(n+1):\n",
        "    xjs[j] = np.cos((2*j+1)/(2*n+2)*np.pi)\n",
        "\n",
        "for ki in range(0, n+1, 2):\n",
        "    aks[ki] = 0\n",
        "    for j in range(n+1):\n",
        "        aks[ki] += 2/(n+1)*poly.eval(xjs[j])*eval_chebyt(ki, xjs[j])\n",
        "\n",
        "aks /= k*np.sqrt(2)\n",
        "aks[0] /=2\n",
        "aks = aks[0: n+1: 2]\n",
        "```"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "749he8U0_lXY"
      },
      "source": [
        "We can then obtain the phases using `oct2py` (with octave installed) and the folder location of `QSPPACK` cloned from https://github.com/qsppack/QSPPACK. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WUapwQuC_lXY"
      },
      "source": [
        "```python\n",
        "from oct2py import octave\n",
        "\n",
        "eps = 0.01\n",
        "folder = 'path to QSPPACK'\n",
        "octave.addpath(folder)\n",
        "opts = octave.struct(\"criteria\", eps)\n",
        "phases, _ = octave.QSP_solver(aks, 0, opts, nout=2)\n",
        "```"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "plF-Dd20_lXY"
      },
      "source": [
        "#### Using pyqsp"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2F5C7BhL_lXZ"
      },
      "source": [
        "For $k=20$, pyqsp with the `laurent` method works properly. However, it often fails for higher order polynomials. In that case, you can install tensorflow and use `method=tf`. It however seemed much slower than `QSPPACK`. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iVG8FmVD_lXZ"
      },
      "source": [
        "```python\n",
        "import pyqsp\n",
        "from pyqsp import angle_sequence\n",
        "from pyqsp.angle_sequence import AngleFindingError\n",
        "from sympy.polys.polytools import primitive\n",
        "from pyqsp.completion import CompletionError\n",
        "\n",
        "# Compute phases for real part Cos(Ht) of Exp(iHt)\n",
        "pg = pyqsp.poly.PolyCosineTX()\n",
        "prefac, poly = primitive(filter_func(x, delta, k))\n",
        "poly = poly.as_poly()\n",
        "polydict = poly.as_dict()\n",
        "\n",
        "coefs = np.zeros(2*k + 1, dtype=float)\n",
        "for term, coeff in polydict.items():\n",
        "    coefs[term[0]] = prefac*coeff/k/np.sqrt(2)\n",
        "\n",
        "n_attempts = 100\n",
        "method = 'laurent'\n",
        "eps=0.01\n",
        "\n",
        "for i in range(n_attempts):\n",
        "    try:\n",
        "        phases = angle_sequence.QuantumSignalProcessingPhases(\n",
        "                coefs, eps=eps, suc=1-eps/10, method=method, tolerance=0.01)\n",
        "    except (AngleFindingError, CompletionError):\n",
        "        if i == n_attempts-1:\n",
        "            raise RuntimeError(\"Real phases calculation failed, increase n_attempts or eps\")\n",
        "        else:\n",
        "            print(f\"Attempt {i+2} for the real coefficients\")\n",
        "    else:\n",
        "        break\n",
        "```"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "qsdkmain",
      "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.9.16"
    },
    "orig_nbformat": 4,
    "colab": {
      "provenance": []
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}