applications/chemistry/qpe_for_molecules/qpe_for_molecules.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "75e2c1fe-8a06-4f0b-8fd6-d1b6d10ce82f",
   "metadata": {},
   "source": [
    "# Quantum Phase Estimation (QPE) for Solving Molecular Energies\n",
    "\n",
    "\n",
    "Quantum Phase Estimation (QPE) is a key algorithm in quantum computing, allowing one to estimate the phase (or eigenvalue) of an eigenvector of a unitary operation. For a given Hamiltonian $H$, and an eigenvalue ${|\\psi\\rangle}$, the output of the algorithm is $\\epsilon$ where\n",
    "\n",
    "$U{|\\psi\\rangle} = e^{2\\pi i\\epsilon}{|\\psi\\rangle} , U = e^{2\\pi iH}$ .\n",
    "\n",
    "\n",
    "Therefore, by measuring the phase accumulated, the QPE algorithm allows calculating the energies relating to the chosen initial state.\n",
    "When using QPE for chemistry problems, it is common to search for the lowest energy of a given molecule. As the molecule can be written in the form of a Hamiltonian (Hermitian matrix representing the energetic forces of the structure), one only needs to insert the ground eigenvector in order to obtain the minimal energy value using QPE. However, obtaining the ground state is not a trivial problem. In order to overcome it, it is sufficient to use a state with big overlap with the ground state.\n",
    "\n",
    "We define a state ${|v\\rangle}$ which will be chosen as the algorithm's initial state. Let us define {$\\psi_i$} to be the set of (unknown) eigenvalues of $H$. Generally, any vector can be rewritten as a superposition of any basis set, thus -\n",
    "\n",
    "${|v\\rangle} = \\sum_i a_i{|\\psi_i\\rangle}$\n",
    "\n",
    "and -\n",
    "\n",
    "$U{|v\\rangle} = \\sum_i a_i e^{2\\pi i\\epsilon_i}{|\\psi_i\\rangle}$.\n",
    "\n",
    "where ${\\epsilon_i}$ are the eigenvalues of $H$, i.e. the span of energies relating to the molecule. Using execution with enough shots, one will obtain this set of $\\epsilon_i$, i.e., a subset of the Hamiltonian's eigenvalues. As we are specifically interested in $\\epsilon_0$, the ground state of $H$, it is important that there is a large overlap between ${\\psi_0}$ and ${|v\\rangle}$ so the probability to measure ${\\epsilon_0}$ is high, i.e.\n",
    "\n",
    "$P(\\epsilon_0) = |\\langle v|\\psi_0\\rangle|^2 > \\zeta $.\n",
    "\n",
    "How large is $\\zeta$? After execution, we will obtain a set of ${E_i}$. If we have done 1000 shots of execution, and $P(\\epsilon_0)>1\\%$, we should sample $\\epsilon_0$ about 10 times.\n",
    "\n",
    "A common choice for ${|v\\rangle}$ (the initial state) is the Hartree-Fock (HF) state, which typically has a large overlap with the ground state. However, other guesses for the initial state are possibly good or even better fit, and choosing the right initial state is a sort of art and an active field of research.\n",
    "\n",
    "For further reading about QPE we recommend [[1](#NC)].\n",
    "\n",
    "\n",
    "**What are the benefits of using QPE algorithm for finding a molecule's ground state?**\n",
    "\n",
    "The two most prominent methods to solve ground energy for molecules are quantum variational algorithm (VQE) and QPE. Those promise better scalability compared to classical counterparts as the molecules becomes more complex, with larger number of electrons, referring to a physical problem with more degrees of freedom.\n",
    "\n",
    "The number of parameters in VQE is closely related to the number of electrons. This may create inherent difficulty achieving high-precision calculations through sampling statistical estimators, and perhaps even not converge for very large systems. On the other hand, the number of parameters in QPE is a flexible value which is directly related to the resolution of the problem, but is not bounded with the number of electrons.\n",
    "\n",
    "Furthermore, it is known that advanced quantum algorithms based on QPE can perform electronic structure calculations in sub-exponential time with accuracy that rivals exact diagonalization methods. This guarantee of simultaneously achieving high accuracy, efficiency, and generality is a feat that is believed to be impossible for classical algorithms. For those reasons, VQE is applicable in the near term (NISQ) era, while QPE is suited for fault-tolerant design.\n",
    "\n",
    "**In this tutorial, we follow the QPE algorithm steps as follows:**\n",
    "\n",
    "1. Define a molecule and convert it into a Hamiltonian.\n",
    "\n",
    "2. Prepare the Hamiltonian for QPE: including normalization and trimming of negligible terms.\n",
    "\n",
    "3. Construct a quantum model: initializing the state for the HF state and leveraging the `qpe_flexible` function.\n",
    "\n",
    "4. Executing the circuit to find the related phases and analyzing the results to find the ground state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "11eaf8f8-c06e-443e-ac74-5aa364f5198c",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T20:32:05.487269Z",
     "start_time": "2025-03-04T20:32:03.276286Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "## Imports\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from classiq import *\n",
    "\n",
    "# for chemistry\n",
    "from classiq.applications.chemistry import Molecule, MoleculeProblem\n",
    "from classiq.applications.combinatorial_helpers.pauli_helpers.pauli_utils import (\n",
    "    pauli_operator_to_hamiltonian,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac537458-568d-4718-9020-408cc7232641",
   "metadata": {},
   "source": [
    "## Defining a molecule with Classiq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00cfd366-f819-440b-90b8-7d03eb27e2a6",
   "metadata": {},
   "source": [
    "In this tutorial we choose to work on the LiH molecule:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2b67b2e1-94c2-4c8f-8ee9-577c79a7234f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T20:32:07.229353Z",
     "start_time": "2025-03-04T20:32:07.224866Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# build your molecule\n",
    "\n",
    "\n",
    "molecule_LiH = Molecule(atoms=[(\"H\", (0.0, 0.0, 0.0)), (\"Li\", (0.0, 0.0, 1.596))])\n",
    "\n",
    "# Example of some other molecules:\n",
    "molecule_H2 = Molecule(atoms=[(\"H\", (0.0, 0.0, 0)), (\"H\", (0.0, 0.0, 0.735))])\n",
    "molecule_O2 = Molecule(atoms=[(\"O\", (0.0, 0.0, 0)), (\"O\", (0.0, 0.0, 1.16))])\n",
    "molecule_H2O = Molecule(\n",
    "    atoms=[(\"O\", (0.0, 0.0, 0.0)), (\"H\", (0, 0.586, 0.757)), (\"H\", (0, 0.586, -0.757))]\n",
    ")\n",
    "molecule_BeH2 = Molecule(\n",
    "    atoms=[(\"Be\", (0.0, 0.0, 0.0)), (\"H\", (0, 0, 1.334)), (\"H\", (0, 0, -1.334))]\n",
    ")\n",
    "\n",
    "molecule = molecule_LiH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5d27e5df-e16e-44b2-b78a-187905577908",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T20:32:16.972965Z",
     "start_time": "2025-03-04T20:32:07.977442Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Hamiltonian is defined on 6 qubits, and contains 231 Pauli strings\n"
     ]
    }
   ],
   "source": [
    "# define your molecule problem\n",
    "\n",
    "gs_problem = MoleculeProblem(\n",
    "    molecule=molecule,\n",
    "    basis=\"sto3g\",\n",
    "    mapping=\"jordan_wigner\",  #'bravyi_kitaev'\n",
    "    z2_symmetries=True,\n",
    "    freeze_core=True,\n",
    ")\n",
    "\n",
    "operator = gs_problem.generate_hamiltonian()\n",
    "gs_problem = gs_problem.update_problem(operator.num_qubits)\n",
    "\n",
    "mol_hamiltonian = pauli_operator_to_hamiltonian(operator.pauli_list)\n",
    "problem_size = len(mol_hamiltonian[0].pauli)\n",
    "\n",
    "print(\n",
    "    f\"The Hamiltonian is defined on {problem_size} qubits, and contains {len(mol_hamiltonian)} Pauli strings\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e1e299b-bd9d-40e8-a1a6-c610e61f7354",
   "metadata": {},
   "source": [
    "Finally, we calculate the ground state energy as a reference solution to our quantum solver (note that the Hamiltonian obtained from the Molecule object excludes the nuclear repulsion energy)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ca1d9c36-2d54-45c4-b3ce-f6ee8abbb28f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T20:32:17.052329Z",
     "start_time": "2025-03-04T20:32:16.972582Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expected energy (without nuclear repulsion energy): -1.0789776863285236 Ha\n"
     ]
    }
   ],
   "source": [
    "mat = hamiltonian_to_matrix(mol_hamiltonian)\n",
    "w, v = np.linalg.eig(mat)\n",
    "classical_sol = np.real(min(w))\n",
    "print(f\"Expected energy (without nuclear repulsion energy): {classical_sol} Ha\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4f40418-500a-4d86-88e4-ca64bdaf24ac",
   "metadata": {},
   "source": [
    "## Preparing the Molecule for QPE"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "034e17f6-6008-4ea7-b5de-af7626461dbb",
   "metadata": {},
   "source": [
    "### Trimming the Hamiltonian\n",
    "\n",
    "As we can see, the Hamiltonian may contain a large number of terms. In many cases we can compress the Hamiltonian by trimming small terms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c7d00124-e791-43f7-b910-bc1031865e37",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T20:44:26.990326Z",
     "start_time": "2025-03-04T20:44:26.711352Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnsAAAHYCAYAAAA1VC5KAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABSI0lEQVR4nO3de3wU5fn///cmgQRyWAinBAgQgyAhSImGgoBFFEVsFESpWAQ8fCwHreIJtfKF1BPSb/WLGqnVfrSKVguKmhYjINg2BQUJh8agAgYMEAUSSMIphM38/uC3W0J2k93NHia7r+fjkYdmdmbnmp2Z7MXc933dFsMwDAEAACAkRQQ7AAAAAPgPyR4AAEAII9kDAAAIYSR7AAAAIYxkDwAAIISR7AEAAIQwkj0AAIAQRrIHAAAQwkj2AAAAQhjJHgAAQAiLCnYAAAD4S1VVlfbv3+/Vtueff74iIyN9HBEQeCR7AICQ9dFHH+mWW27xatuDBw+qY8eOPo4ICDyacQEAIatVq1bBDgEIOpI9AEDImjBhgrKzsx2/P/zwwzIMw60fnuohVJDshbH58+fLYrEEOwyXNm7cqEsuuUSxsbGyWCzasmWLy+Wvv/66LBaLdu/e7dE+vN3O7IJ1XK7OWbhx9vkH+1oL9v6DJSoqSn/96181cuRISdKCBQu0cOHC4AYVIAsXLtQFF1ygurq6YIfiU+F8f/3hD39Qjx49VFNT49F2JHte+s9//qMbbrhBPXv2VExMjLp166bRo0frhRde8Ol+1q1bp/nz5+vIkSM+fV+zq62t1Y033qiKigo999xzevPNN9WzZ0+Xy80sXM5hoM+N/Y/rl19+6fT1kSNHKiMjw2/7b65wuS4k6ZNPPpHFYtHbb7/t9PVrr71WsbGxfktKYmJilJeXp8GDB0uS5syZo5dfftkv+3LGfvz2n8jISPXq1UuzZ8/W0aNH/bLPqqoqPfPMM5ozZ44iIup/1fsqHvs9aP+JiYlRnz59dNddd+nHH3/09SEFjJnvzWnTpunUqVOeX78GPPbvf//baN26tdG7d2/j8ccfN1555RXj//yf/2NceeWVRlpamk/39bvf/c6QZJSUlPj0fQ3DMObNm2eY9RLYvn27Icl45ZVX3Fp++vRp48SJE0ZdXZ1H+/F2O0/48xy6EojjOperc+Mvr732miHJ2Lhxo9PXf/aznxn9+/cPSCzncvb52+O1XweBvi6CcU3YPfPMM4Yk45tvvnH6eteuXY3Bgwf7PY7y8nIjIyPDkGREREQYf/nLX/y+T8P47/E/99xzxptvvmm88sorxtixYw1Jxp133umXfT733HNGQkKCceLECb/FY7+mf/vb3zreZ+rUqUZERISRmppqHDt2zJeH1GC/Z987vry+vbk3A3l/PfTQQ0bPnj092hejcb3w5JNPymq1auPGjWrXrl291w4cOOCTfRw7dkyxsbE+ea+WyP45uvp8z10eGRnpVYkEb7czu2Acl6tz0xwt9T4w43UVzJi2bdumhIQEnX/++Q1e++GHH7R//35dc801fo8jMTFRK1eu1IgRI7Rr1y5NmTJF8fHxft/3tm3bFBsbq1//+teOp2y33nqr0tLS9OGHH/rlKeNrr72ma6+9VjExMX6P5+qrr9bFF18sSbrjjjvUoUMHPfvss/rwww81adKk5h+MG4J1fdv/RgVy/xMnTtTChQu1du1ajRo1yq1taMb1wq5du9S/f3+nX2qdO3eu9/vmzZt19dVXKyEhQXFxcbr88sv1+eef11vH3neuuLhYN998s9q3b6/hw4dr/vz5evDBByVJqampjkflZ/cJ2Ldvn2677TZ16dJF0dHR6t+/v/73f/+3QVwFBQXKyspSTEyM0tLSPL6Z9+3bp9tvv11du3ZVdHS0UlNTNWPGDJ06dcqjY3Un5mnTpulnP/uZJOnGG2+UxWLRyJEjXS6XXPeXaCpuZ9u585naz9nOnTs1bdo0tWvXTlarVbfeequOHz9eb73GzmF1dbXuvfde9erVS9HR0ercubNGjx6twsJCl+fCnW3OPa7du3fXa24598eT8+NMY+dGat594Et79uzRzJkz1bdvX7Vp00YdOnTQjTfe2OC6scfy7bffavLkybJarerUqZPmzp0rwzBUWlqq6667TgkJCUpKStLvf//7ets31X/HnXvbk8+sqevQWUyefhbu7MOVrVu3atCgQU77CG/atEmSNHDgQLfeq7mSk5O1evVqdevWTbW1tbrhhhv0j3/8w61t3333XWVmZqpNmzbq16+fVq9eLcMw1L9/fz355JMut9u6dasuvPDCes2pkZGR6ty5s6qrq5t9TOcqKSnRtm3bdMUVVwQlHnsCUlJSIsn9a23atGnq1atXg/dzp3+5u33mmvr72di92djfKGf79+Te+eyzz3TxxRfX+452ddwXXXSREhMT9eGHHzZ6rGfjyZ4XevbsqfXr16uoqKjRPkFfffWVRowYoYSEBD300ENq1aqVXn75ZY0cOVL/+Mc/9NOf/rTe+jfeeKPOP/98PfXUUzIMQ8OHD9e3336rv/zlL3ruueccI8M6deokSfrxxx81ZMgQWSwW3XXXXerUqZM+/vhj3X777aqqqtK9994r6Uz/wiuvvFKdOnXS/Pnzdfr0ac2bN09dunRx63j379+vwYMH68iRI7rzzjt1wQUXaN++fVq2bJmOHz+u1q1bu32s7sT8q1/9St26ddNTTz2lX//618rKylKXLl0UFxfndHlz4j6Xu5+p3cSJE5Wamqqnn35ahYWFevXVV9W5c2c988wzkqTrr7++0XM4ffp0LVu2THfddZfS09NVXl6ugoICbd++XZmZmU6Py5ttOnXqpDfffLPestraWs2ePbve5+Dp8du5OmdS8++DplRWVurQoUMNltfW1jZYtnHjRq1bt0433XSTunfvrt27d2vx4sUaOXKkiouL1bZt23rr/+IXv1C/fv20YMEC/f3vf9cTTzyhxMREvfzyyxo1apSeeeYZvfXWW3rggQeUlZWlSy+9tMl4paavC08/s6auQ2c8/Sy82YcknTp1St98840mTZqknTt3Nnh9zZo1kqQLL7zQ6fa1tbWqrKxsdB92iYmJDfqnOdOrVy+tWrVKl156qQ4dOqTs7GytWbPG8XTKmccee0xPPvmkpk2bpttvv12LFi3SlClT9Morr2jv3r266667nG5nP/7/+Z//qbf8xx9/1FdffeXynm2OdevWSZLT9w5EPLt27ZIkdejQQZLn15o/NfX3s6l7U/L8b1RT987mzZs1ZswYJScnKycnRzabTb/97W/r7fNcmZmZ+ve//+3+gfunRTm0rVy50oiMjDQiIyONoUOHGg899JDxySefGKdOnaq33rhx44zWrVsbu3btcizbv3+/ER8fb1x66aWOZfa+c5MmTWqwr8b6Dtx+++1GcnKycejQoXrLb7rpJsNqtRrHjx93xBETE2Ps2bPHsU5xcbERGRnpVp+9KVOmGBEREU77Rtn7DLh7rO7GvHbtWkOSsXTp0nrruVrurA+HO3Gfu5278dnP2W233VZvvfHjxxsdOnSot6yxc2i1Wo1Zs2Y1WN4Yd7Zx9nmca+bMmUZkZKSxZs0axzJ3j98ZV+fGF/dBY8fY2M+5ffacxb9+/XpDkvHGG280iOXs/kunT582unfvblgsFmPBggWO5YcPHzbatGljTJ06tUFsZ3/+nvTZ8/Qzc+c6PHf/nn4W7uzDmc2bNzd5niQZR44ccbq9/bpy58fT/o+bNm0yEhISDElGhw4djK+++srpev/85z8NScacOXMcy5YtW2ZIMjIyMuotd3X8v//9742DBw8a+/btM1auXGlcfPHFRmRkpLFq1SqPYnbHY489Zkgyqqur/RqP/ZpavXq1cfDgQaO0tNR45513jA4dOhht2rQx9u7daxiG+9fa1KlTjZ49ezZY99z+5e7cX6648/fT1b3Z2N8oZ/t3997Jzs422rZta+zbt8+xbMeOHUZUVJTL7+g777zTaNOmTaPHcTaacb0wevRorV+/Xtdee622bt2qhQsX6qqrrlK3bt300UcfSZJsNptWrlypcePG6bzzznNsm5ycrJtvvlkFBQWqqqqq977Tp093OwbDMPTee+8pOztbhmHo0KFDjp+rrrpKlZWVKiwslM1m0yeffKJx48apR48eju379eunq666qsn91NXV6YMPPlB2drbTf/VaLBa3j9XdmH3BnbjP5U18556zESNGqLy8vMG5daVdu3b64osvPJrOyZttzvXGG2/opZde0sKFC3XZZZdJ8u74m+Lv+0CScnNztWrVqgY/zp4UtWnTxvH/tbW1Ki8vV+/evdWuXTunx3bHHXc4/j8yMlIXX3yxDMPQ7bff7ljerl079e3bV999951Hcbvii8/MnevQ08/C22t927ZtkqT/9//+n5YuXdrgp0OHDurVq5esVqvT7QcOHOj0/Dr7SUpKajSWc2VmZmrs2LGSpPLycr3//vtO11u0aJHat2+vxx57zLHM3ty4c+dOzZ4927G8pqZGXbp0cXwu9uO///771alTJ3Xr1k1XXnmlqqur9dFHHzmaWs/drjnKy8sVFRWluLi4Bq/5I54rrrhCnTp1UkpKim666SbFxcVp+fLl6tatmyTPrzV/8sXfT0//RjV279hsNq1evVrjxo1T165dHev07t1bV199tcv3bN++vU6cOOF2Vwqacb2UlZWl999/X6dOndLWrVu1fPlyPffcc7rhhhu0ZcsWJSYm6vjx4+rbt2+Dbfv166e6ujqVlpaqf//+juWpqalu7//gwYM6cuSI/vjHP+qPf/yj03UOHDiggwcP6sSJE047Rvft21crVqxocj9VVVWNNlcfPHjQrWPt1KmTWzH7gjtxO9vG0/jOTqClMzegJB0+fFgJCQlN7nPhwoWaOnWqUlJSdNFFF2ns2LGaMmVKvS95X2xzti1btmj69OmaNGmS7rvvPsdyb46/Ke5eG97eB5I0ePBgpwl9+/btGzTvnjhxQk8//bRee+017du3r14TjLOmwnPPr9VqVUxMTINiu1arVeXl5R7F7Yo3n5k312FzPwt3r/WtW7cqMjJS06dPV3R0dIMYDh8+3GjfzPbt27vse9ZcTzzxhN555x1J0tSpU+slc3b25Ds7O9tp8nTrrbfW604SHR1dr+zI1q1bFRUVpRUrVshisah169ZKSUlpcJ2fu52/+COe3Nxc9enTR1FRUerSpYv69u1brznd02vNn5r791Py/G9UY/fOsWPHdOLECfXu3bvBds6W2dk/Q3dr5ZLsNVPr1q2VlZWlrKws9enTR7feequWLl2qX/3qVx6/19n/+mmKvR7V5MmTNXXqVKfrXHjhhaYqpuluzMHiTXyuRl+d/cesMRMnTtSIESO0fPlyrVy5Ur/73e/0zDPP6P3333f5rzpvtrE7fPiwJkyYoD59+ujVV1+t95pZzo8n94Gn7r77br322mu69957NXToUFmtVlksFt10001O7xVn57e559wfvInJF59FU/uQzjxJSktLa5DoSdL27dtVV1fX6HV16tQpVVRUNLoPu06dOrk9IvLVV1/V3LlzJUnXXHNNg/vB7rvvvlN1dXWDvmwHDx6UJM2aNavR/Wzbtk29e/fW6NGj3YrLFzp06KDTp0+rurpa8fHxfo/H1T+47Ny91lwlLjabzWexNufvp52nf6P88Tfj8OHDatu2rduxkOz5kP1iLysrU6dOndS2bVt98803Ddb7+uuvFRERoZSUlCbf09XF36lTJ8XHx8tmszX6r16bzaY2bdpox44dDV5zFpuz/SQkJKioqKjRddw51tjYWLdi9gV34na2jT/ia+pfXsnJyZo5c6ZmzpypAwcOKDMzU08++WSjf3i82aaurk6//OUvdeTIEa1evbpBp2h/HL+v7gNfWbZsmaZOnVpv9OzJkyeDUjy1sXs7EJ9ZoD6Lbdu2adiwYU5f++qrryQ1PhJ33bp1jq4GTSkpKXE6mvNcH374oaNp7ZJLLtHSpUsVFeX869Ce1J37NPfpp592unzRokXatm2b/vSnP0k6c/wjRoxoMqZzt2uOCy64QNKZz+PcRDoY8bh7rbVv397p9bdnz55mx3C2pv5+BnJmqc6dOysmJsbp4CVny+xKSkrUr18/t/dDnz0vrF271mlGbm8S7du3ryIjI3XllVfqww8/rDcU+8cff9Tbb7+t4cOHu9XMZ68xdu4NEBkZqQkTJui9995zmtDY/0BFRkbqqquu0gcffKDvv//e8fr27dv1ySefNLn/iIgIjRs3Tnl5eU5nKjAMw+1jdTdmX3An7nP5Kz5X59BmszVowujcubO6du3qciocb7axy8nJ0SeffKK//OUvTpsh/HH8vroPfCUyMrLBuX/hhRd8+uTAXY3d24H4zALxWfzwww86cOBAvSbns7mT7Pm6z15BQYFuuukm2Ww29e/fX3/7298afTpi70t49j3x9ttv65///KekM0nL2bZt2+ZIsOzHn56e3mRcZ2/XXEOHDpWkBn/7ghWPu9daWlqaKisrHf0KpTMPT5YvX+6TONz9++nq3vSHyMhIXXHFFfrggw/q9SPcuXOnPv74Y5fbFRYW6pJLLnF7PzzZ88Ldd9+t48ePa/z48brgggt06tQprVu3Tu+++6569eqlW2+9VdKZ/iCrVq3S8OHDNXPmTEVFRenll19WTU2N23MzXnTRRZKk3/zmN7rpppvUqlUrZWdnKzY2VgsWLNDatWv105/+VP/zP/+j9PR0VVRUqLCwUKtXr3Y0feTk5Cg/P18jRozQzJkzdfr0ab3wwgvq379/vZvKlaeeekorV67Uz372M915553q16+fysrKtHTpUhUUFKhdu3ZuH6u7MfuCO3Gfyx/xuTqHtbW16t69u2644QYNHDhQcXFxWr16tTZu3NigbptddXW1x9tIZ8rvPP7447r00kt14MABLVmypN7rkydP9tvx++I+8JWf//znevPNN2W1WpWenq7169dr9erVjhIRgdTYvR2IzywQn8XWrVslqdFkLy4uTmlpaS7fw5d99oqKipSdna2TJ08qJSVF+fn5jv5TrvTr10+pqal6/vnn1bZtW0VERGjBggWaOHGi/vrXv2r+/Pm67777NGDAAElnkiT7/dTU8Z/t7O2a67zzzlNGRoZWr16t2267zbE8WPG4e63ddNNNmjNnjsaPH69f//rXOn78uBYvXqw+ffr4ZCCHu38/Xd2b/jJ//nytXLlSw4YN04wZM2Sz2fTiiy8qIyPD6fzimzZtUkVFha677jr3d+L2uF04fPzxx8Ztt91mXHDBBUZcXJxj6rS7777b+PHHH+utW1hYaFx11VVGXFyc0bZtW+Oyyy4z1q1bV28d+/DsgwcPOt3f448/bnTr1s2IiIhoMLT7xx9/NGbNmmWkpKQYrVq1MpKSkozLL7/c+OMf/1jvPf7xj38YF110kdG6dWvjvPPOM/7whz94NF3anj17jClTphidOnUyoqOjjfPOO8+YNWuWUVNT49GxuhuzL0qvuBO3s+3cic/VOXMVh7NzWFNTYzz44IPGwIEDjfj4eCM2NtYYOHCg8dJLL7k8D+5uc24cTZWvOJu719S5XJ0bw/DNfXAub6ZLO3z4sHHrrbcaHTt2NOLi4oyrrrrK+Prrr42ePXvWK53iKpapU6casbGxTe7L3dIQjd3bzfnM3Nl/cz8Ld0pdLFy40JBkbNu2zenrqampxpAhQ1xu70t79uwxunXr5iizsn37dre33bx5szFkyBAjOjraaN++vfGb3/zGqKurM2677TYjKirKeP311w3DMAybzWa0bdvWUbqoqeO3O3c7X3j22WeNuLi4emVPfB1PU/egnbvXmmGcKW2WkZFhtG7d2ujbt6+xZMkSn5Ve8eRvrrN7s7G/UY2VXnHn3vn000+NQYMGGa1btzbS0tKMV1991bj//vuNmJiYBvuaM2eO0aNHD4+mS7MYRhB7FQMA4Gfl5eUaPny4vv76a8XGxurTTz9tUJjaF7799luNHDnS47Ie3m7XmMrKSp133nlauHBhvVJBwYoHnhs3bpy++uqren3ua2pq1KtXLz388MO655573H4v+uwBAELWsWPHdM011+jrr79WVFSUli5d6pdET/K+n5sv+8fZWa1WPfTQQ/rd737ncVUGf8SDxp04caLe7zt27NCKFSvqTTspnZnzuFWrVh7X+iPZAwCErHvuuUdffPGFJOn06dMaO3Zso/NEn/3z+uuve7Sv//znP14lSd5u15Q5c+Y4RnCbIR64dt555+mRRx7RK6+8oscee0xDhgxR69at9dBDD9Vbb/r06fr++++dljJqDM24AICQNWbMGLcqDzjz2muvadq0ab4NCHDi1ltv1dq1a/XDDz8oOjpaQ4cO1VNPPeWzuZNJ9gAAIWvv3r06evSoV9smJye7nMYNaElI9gAAAEIYffYAAABCGEWVTaiurk779+9XfHx8QKdtAQAA3jMMQ9XV1eratavHA2P8iWTPhPbv3x/Q+UIBAIDvlJaWqnv37sEOw4Fkz4Ti4+MlnblYAjlvKAAA8F5VVZVSUlIc3+NmQbJnQvam24SEBJI9AABaGLN1wTJPgzIAAAB8jmQPAAAghJHsAQAAhDCSPQAAgBBGsgcAABDCSPYAAABCGMkeAABACCPZAwAACGEkewAAACGMGTTCiK3O0IaSCh2oPqnO8TEanJqoyAhzVfkGAAC+RbIXJvKLypSTV6yyypOOZcnWGM3LTteYjOQgRgYAAPyJZtwwkF9UphlLCuslepL0Q+VJzVhSqPyisiBFBgAA/I1kL8TZ6gzl5BXLcPKafVlOXrFsdc7WAAAALR3JXojbUFLR4Ine2QxJZZUntaGkInBBAQCAgCHZ87GdO3dq7NixiouLU8eOHTVz5kwdO3YsaPEcqHad6HmzHgAAaFkYoOFDlZWVGjVqlLp27aqlS5eqoqJC9913n3788Ue99957QYmpc3yMT9cDAAAtC8meD7388ss6ePCgvvzyS3Xu3FmS1KZNG02YMEGbNm3SRRddFPCYBqcmKtkaox8qTzrtt2eRlGQ9U4YFAACEHppxfWjFihUaNWqUI9GTpGuvvVZxcXH629/+FpSYIiMsmpedLulMYnc2++/zstOptwcAQIgK+WTvm2++0QsvvKBp06ZpwIABioqKksVi0RNPPOHW9kuXLtXIkSPVvn17xcbGauDAgVq4cKFqa2sbrFtcXKx+/frVWxYVFaU+ffpo+/btPjkeb4zJSNbiyZlKstZvqk2yxmjx5Ezq7AEAEMJCvhl38eLFWrRokVfb3nvvvVq0aJGioqI0atQoxcXFac2aNZozZ47y8vK0cuVKtWnTxrH+4cOH1a5duwbv0759e1VUBHe065iMZI1OT2IGDQAAwkzIP9nLyMjQAw88oLfeekvbt2/XLbfc4tZ2H3zwgRYtWqS4uDh98cUX+uSTT/Tee+9px44dGjBggAoKCjR37lw/R+9bkREWDU3roOt+0k1D0zqQ6AEAEAZC/sneHXfcUe/3iAj38tunnnpKkvTwww8rMzPTsbxjx4566aWXNGLECL344ouaO3eurFarpDNP8I4cOdLgvQ4fPqzzzz/fyyMAAADwXsg/2fPGvn37tHHjRknSzTff3OD14cOHKyUlRTU1NVqxYoVjeb9+/Rr0zbPZbPr2228b9OUDAAAIBJI9JzZv3ixJSkxMVGpqqtN1Lr744nrrStLYsWO1du1aHTx40LEsLy9PR48e1TXXXONyfzU1Naqqqqr3AwAA4Aske06UlJRIknr06OFynZSUlHrrStKvfvUrJSYm6rrrrtPHH3+st99+W9OnT9d1113nSA6defrpp2W1Wh0/9vcGAABoLpI9J6qrqyVJsbGxLteJi4uTpHpP4dq1a6c1a9YoISFBN9xwg+6++26NGzdOS5YsaXR/jzzyiCorKx0/paWlPjgKAACAMBigEWh9+vRRfn6+R9tER0crOjraTxEBAIBwxpM9J+Lj4yVJx44dc7nO0aNHJUkJCQkBiQkAAMAbPNlzolevXpLUaHOq/TX7ui2Jrc7QhpIK/VB5QhXHTikxLlpJCRRZBgAgFJHsOTFo0CBJUnl5uUpKSpyOyP3yyy8lqV4NvpYgv6hMOXnFKqs82eC1ZGuM5mWnM30aAAAhhGZcJ7p3766srCxJ0ttvv93g9YKCApWWlio6Olpjx4712X5zc3OVnp7u2Lev5ReVacaSQqeJniSVVZ7UjCWFyi8q88v+AQBA4JHsufDoo49KkhYsWKDCwkLH8vLycs2cOVOSdNdddzlmz/CFWbNmqbi42FHQ2ZdsdYZy8opluLFuTl6xbHXurAkAAMzOYhhGSH+rFxYWOpIzSdq1a5cOHTqk7t27q1u3bo7ly5cvV3Jy/ebLe+65R88//7xatWqlyy+/XLGxsfr000915MgRDRs2TKtWrVKbNm18HnNVVZWsVqsqKyt9NgBk/a5yTXrlc7fX/8v/DNHQtA4+2TcAAOHAH9/fvhDyffaqqqr0xRdfNFi+d+9e7d271/F7TU1Ng3UWLVqkYcOGKTc3V+vWrVNtba3S0tL08MMPa/bs2WrdurVfY/elA9XOm259tT4AADCnkE/2Ro4cqeY8vJw4caImTpzow4iCo3N8jF/XBwAA5kSfvTAxODVRydamEziLzozKHZya6P+gAACA35HsmYg/R+NGRlg0Lztd7lTRm5edTr09AABCRMgP0GiJ/NnBkzp7AAD4BwM0YApjMpI1Oj2JGTQAAAgTJHthKDLCQlkVAADCBH32AAAAQhjJHgAAQAgj2QMAAAhhJHsm4s/SKwAAIDxResWEzDp0GwAAuGbW72+e7AEAAIQwkj0AAIAQRrIHAAAQwkj2AAAAQhgzaEC2OoPp0wAACFEke2Euv6hMOXnFKqs82eC1ZGuM5mWna0xGchAiAwAAvkAzrokEus5eflGZZiwpdJroSVJZ5UnNWFKo/KKygMQDAAB8jzp7JhSIOj22OkPDn1njMtGzs0hKssaoYM4omnQBAGgEdfZgKhtKKppM9CTJ0JknfBtKKvwfFAAA8DmSvTB1oLrpRK856wMAAHMg2QtTneNj/Lo+AAAwB5K9MDU4NVHJ1qYTOIvOjModnJro/6AAAIDPkeyFqcgIi+Zlp8udIRfzstMZnAEAQAtFshfGxmQka/HkTJdP+JKtMVo8OZM6ewAAtGAUVQ5zYzKSNTo9iRk0AAAIUSR7JpKbm6vc3FzZbLaA7jcywqKhaR0Cuk8AABAYFFU2IbMWZQQAAK6Z9fubPnsAAAAhjGQPAAAghJHsAQAAhDCSPQAAgBBGsgcAABDCSPYAAABCGHX20ICtzqDIMgAAIYJkD/XkF5UpJ69YZZUnG7yWbI3RvOx0pk8DAKAFoRkXDvlFZZqxpNBpoidJZZUnNWNJofKLygIcGQAA8BbJHiSdabrNySuWO9Op5OQVy1bHxCsAALQEJHsmkpubq/T0dGVlZQV83xtKKlw+0TuboTNP+DaUVPg/KAAA0GwkeyYya9YsFRcXa+PGjQHf94HqphO95qwPAACCg2QPkqTO8TF+XR8AAAQHyR4kSYNTE5VsbTqBs+jMqNzBqYn+DwoAADQbyR4kSZERFs3LTpc7VfTmZadTbw8AgBaCZA8OYzKStXhypssnfMnWGC2enEmdPQAAWhCKKqOeMRnJGp2exAwaAACECJI9NBAZYdHQtA7BDgMAAPgAzbgAAAAhjGQPAAAghNGMi0bZ6gz67wEA0IKR7MGl/KIy5eQVO51GLdkao3nZ6YzMBQDA5GjGhVP5RWWasaTQ5Xy5ZZUnNWNJofKLygIcGQAA8ATJHhqw1RnKySuW4ca6OXnFstW5syYAAAgGkj00sKGkwuUTvbMZOvOEb0NJhf+DAgAAXiHZM5Hc3Fylp6crKysrqHEcqG460WvO+gAAIHBI9kxk1qxZKi4u1saNG4MaR+d459OlubL70HE/RQIAAJqLZA8NDE5NdDk/rjP/b/W3DNQAAMCkSPbQQGSERfOy0+VJFT0GagAAYE4ke3BqTEayFk/OdOsJHwM1AAAwL5I9uDQmI1kFc0bprst6u7U+AzUAADAfkj00KjLComG9O7q1rqcDOwAAgP+R7KFJ9gEbrvrwWXRm+rTBqYmBDAsAALiBZA9Nsg/YkNQg4bP/Pi87XZERngzpAAAAgUCyB7fYB2wknTNgI8kao8WTMzUmIzlIkQEAgMZEBTsAtBxjMpI1Oj1JG0oqdKD6pDrHn2m65YkeAADmRbIHj0RGWDQ0rUOwwwAAAG6iGRcAACCEkewBAACEMJI9AACAEEayBwAAEMJI9gAAAEIYyR4AAEAII9kDAAAIYSR7AAAAIYxkDwAAIISR7JlIbm6u0tPTlZWVFexQAABAiLAYhmEEOwjUV1VVJavVqsrKSiUkJAQ7HAAA4Aazfn/zZA8AACCEkewBAACEMJI9AACAEEayBwAAEMJI9gAAAEIYyR4AAEAII9kDAAAIYSR7AAAAIYxkDwAAIISR7AEAAIQwkj0AAIAQRrIHAAAQwkj2AAAAQhjJHgAAQAgj2QMAAAhhUcEOAC2Trc7QhpIK/VB5QhXHTqld29Y6ctx3/02Mi1ZSQowGpyYqMsIS7MMFAKDFItmDx/KLypSTV6yyypN+31eyNUbzstM1JiPZ7/sCACAU0YwLj+QXlWnGksKAJHqSVFZ5UjOWFCq/qCwg+wMAINSQ7MFttjpDOXnFMgK8X0PS/I++kq0u0HsGAKDlI9mD2zaUVATsid65fqiq0YtrdgZl3wAAtGQke3DbgergJHp2z63+luZcAAA8RLIHt3WOjwl2CMrJK6Y5FwAAD5DswW2DUxOVbA1uwldWeVLzPyrS8s37tH5XOYkfAABNINnzsZ07d2r69OnKzMxUq1at1KtXr2CH5DORERbNy05XsKvevfn595r97hZNeuVzDX9mDU27AAA0gmTPx7766iv97W9/U69evZSRkRHscHxuTEayFk/ODPoTPjtKswAA0DiLYRi0g/lQXV2dIiLO5NDTp09Xfn6+du/e7dF7VFVVyWq1qrKyUgkJCX6Isvn8MYNGxbEavbh2lypP1HoUi0VSkjVGBXNGMdsGACBozPr9zQwaPmZP9EJdZIRFQ9M6+Px9UxLbavqSQo+2MXTmCd+Gkgq/xAQAQEvW4jKTb775Ri+88IKmTZumAQMGKCoqShaLRU888YRb2y9dulQjR45U+/btFRsbq4EDB2rhwoWqrfXsaRL8Y0xGsm4f1surbYNdGgYAADNqcU/2Fi9erEWLFnm17b333qtFixYpKipKo0aNUlxcnNasWaM5c+YoLy9PK1euVJs2bXwcMTx1RXqS/vTv3R5vZ4bSMAAAmE2Le7KXkZGhBx54QG+99Za2b9+uW265xa3tPvjgAy1atEhxcXH64osv9Mknn+i9997Tjh07NGDAABUUFGju3Ln1tnn99ddlsVia/Fm2bJk/DjVseVrixSIp2RqjwamJ/gsKAIAWqsU92bvjjjvq/e5uH7mnnnpKkvTwww8rMzPTsbxjx4566aWXNGLECL344ouaO3eurFarJGn8+PEaMmRIk+/drVs3d8OHG+wlXmYsKXR7Ht552ekMzgAAwIkWl+x5Y9++fdq4caMk6eabb27w+vDhw5WSkqLS0lKtWLFCkyZNkiRZrVZH4ofAspd4yckrbnQ+3mRrjOZlp2tMRnIAowMAoOUIi2Rv8+bNkqTExESlpqY6Xefiiy9WaWmpNm/e7Ej2EFxjMpI1Oj3JaYmXxLhoJSWcabrliR4AAK6FRbJXUlIiSerRo4fLdVJSUuqt663jx49rxYoVkqTvvvtOx48fd/Tpy8rKUs+ePRtsU1NTo5qaGsfvVVVVzYohlPirxAsAAOEiLJK96upqSVJsbKzLdeLi4iQ1P9E6cOCAbrzxxnrL7L+/9tprmjZtWoNtnn76aeXk5DRrvwAAAM6ERbIXSL169ZKnk5I88sgjuu+++xy/V1VVOZ40AgAANEdYJHvx8fGSpGPHjrlc5+jRo5IUlOlNoqOjFR0dHfD9AgCA0BcWyV6vXr0kSaWlpS7Xsb9mXxfmZp+b90D1SXWOZ6AGAACuhEWyN2jQIElSeXm5SkpKnI7I/fLLLyWpXg0+mFN+UVmDkiyUYAEAwLkWN4OGN7p3766srCxJ0ttvv93g9YKCApWWlio6Olpjx44NdHgOubm5Sk9Pd8SKhvKLyjRjSWGD2ntllSc1fUmhcj4q0vpd5bLVedZvEgCAUBUWyZ4kPfroo5KkBQsWqLCw0LG8vLxcM2fOlCTdddddQS2iPGvWLBUXFzsKQKM+W52hnLziRmfVeG3dHk165XMNf2aN8ovKAhYbAABmZTE8HToaZIWFhY7kTJJ27dqlQ4cOqXv37vWmLVu+fLmSk+s36d1zzz16/vnn1apVK11++eWKjY3Vp59+qiNHjmjYsGFatWqV2rRpE7BjcaWqqkpWq1WVlZVBGTBiVut3lWvSK5+7vb5F0uLJmTTtAgACwqzf3y2uz15VVZW++OKLBsv37t2rvXv3On4/u0ix3aJFizRs2DDl5uZq3bp1qq2tVVpamh5++GHNnj1brVu39mvsaJ4D1a6nTXMlJ69Yo9OTGLwBAAhbLS7ZGzlypMd17M42ceJETZw40YcRIVA6x8d4tL6hM335NpRUMAsHACBshU2fPbR8g1MTlWz1LOGTvHsiCABAqCDZMxFG4zYuMsKiednpHm/n6RNBAABCSYsboBEOzNrB0yxWbNuvu/6yWU1VV7FISrLGqGDOKPrsAQD8zqzf3zzZQ4sz9sKuenGSe8Wv52Wnk+gBAMIayR5apLEXJusPkzNd9uFLtsZQdgUAALXA0biA3ZiMZI1OT9KGkgr9UHlCFcdOKTEuWkkJzJULAIAdyR5atMgIC2VVAABoBM24AAAAIYxkz0QovQIAAHyN0ismZNah2wAAwDWzfn/zZA8AACCEMUADIcVWZzA6FwCAs5DsIWTkF5UpJ69YZZUN58JNtsZoXnY6dfcAAGHHo2Rv1KhRPg9g2rRpmjJlis/fF+Elv6hMM5YUylUH1LLKk5qxpJBCywCAsONRsvfZZ5/5PICRI0f6/D0RXmx1hnLyil0memfLySvW6PQkmnQBAGHDo2SvpKTE5wG0a9fO5++J8LKhpMJp0+25DJ15wvf6v0s0bVgqCR8AICx4lOz17NnTX3FAZ+rs5ebmymazBTuUFuVAddOJ3tke//t2vVpQQh8+AEBYoM6eCZm1To9Zrd9VrkmvfO7VtrOvOF93jTqfp3wAgGYz6/e3z0bjHjhwQJ07d/bV2wFuG5yaqGRrjFtNued6bvUO/amgRFemd9Gw8ztRpgUAEHJ8VlS5a9eu2rx5syTpt7/9rfLy8rRnzx5fvT3gUmSERfOy0+VtelZ18rSWFe7T7He3aNIrn2v4M2uUX1Tm0xgBAAgWnzXj/v3vf9dFF12kpKQkpaSkaN++fbJYLEpISNCFF16ogQMHOn4yMjIUExPji92GJLM+Bja7xursecoiUaYFAOARs35/+63PXkVFhbZu3Vrvp7i4WKdOnVJkZKRqa2v9sduQYNaLpSWw1Rl6/d8levzv25v1PhZJSdYYFcwZRZMuAMAtZv3+9tsMGomJibrssst02WWXOZadPn1a27dv17Zt2/y1W4S5yAiLpg1L1asFJc16wmcv07KhpEJD0zr4LkAAAAKsWX32SktLPVo/KipKAwYM0C9/+cvm7BZoVHP78J3N07IuAACYTbOSvb59++r+++9XeXm5r+IBfGJMRrIWT85UUkJ0s96nczx9SwEALVuzkr0xY8boueeeU1pamh5//HEdO3bM6XovvviiHn300ebsCvDYmIxk/fvhyzX7ij4eb2uRlGw9U4YFAICWrFnJ3vvvv69NmzZpxIgRmj9/vtLS0vTCCy/UG3xRU1OjJUuWaOHChc0ONtTl5uYqPT1dWVlZwQ4lZERGWHTPFefrD5MzlWz17CndvOx0BmcAAFq8Zo3GtdlsWrRokT788ENt3LhRJ0+elMViUadOndS5c2dFRkZq7969Ki8v109+8hMVFhb6MvaQZdbRPC2drc7QhpIK/VB5Qv/eeUirth9Q5YmGo8KTrTFMpQYA8JhZv7+bNRr3ySefVE5OjgzDUNu2bdWtWzdJ0vHjx1VUVCSLxaKoqChdcskleumll3wSMOCtyAiLY2Tt+Mzu9ZK/imOnlBgXzQwaAICQ06xk780331SbNm20fPlyjR492rH89OnTevbZZ/XEE0/o+PHjuv7663XhhRc2O1jAl85O/gAACFXN6rO3d+9eXXnllfUSPelMiZWHHnpIW7du1YUXXqgHH3xQS5cubVagAAAA8Fyzkr0uXbro+PHjLl9PTU3VqlWrlJiYqN/+9rfN2RUAAAC80Kxk78orr9TatWu1ZcsWl+t06NBBo0aN0nfffdecXQEAAMALzUr2HnroIUVGRuryyy/XG2+84XK9HTt2KCKiWbsCAACAF5qVgfXu3VsvvPCCjhw5oltvvVUDBgxQTk6OPvvsM+3YsUMbN27UtGnTtGXLFvXv399XMQMAAMBNzRqNK0m333672rVrp7vvvltfffWViouL671uL+N33333NXdXAAAA8JBP2lYnTJig3bt36y9/+YsmT56svn37KiYmRm3bttWQIUO0dOlSTZw40Re7AgAAgAeaNYMGfCs3N1e5ubmy2Wz69ttvTVeBGwAAuGbWGTRI9kzIrBcLAABwzazf3wyRBQAACGEeDdAYNWqUzwOYNm2apkyZ4vP3BQAAgIfJ3meffebzAEaOHOnz9wS8YasztKGkQj9UnlDFsVNq17a1jhz/738T46KVlBCjwamJioywBDtcAADc4lGyV1JS4vMA2rVr5/P3BDyVX1SmnLxilVWebHLdZGuM5mWna0xGcgAiAwCgeRigYUJm7eAZqvKLyjRjSaE8uREskhZPziThAwA4mPX7u9lFlYGWzFZnKCev2KNET5IMSQ8t26rSihNKjKWZFwBgXiR7CGsbSircarp1puqkTU+u2N5guTUmSqPTu2hoWkf6/AEAgo5kD2HtQLV3iV5jKk+e1rLCfVpWuM/p6/T5AwAEEnX2ENY6x8cEfJ9llSc1Y0mh8ovKAr5vAED4IdlDWBucmqhka+ATPknKySuWrY7xUQAA/yLZQ1iLjLBoXna6At2DztCZJ3wbSioCvGcAQLgh2UPYG5ORrMWTM4PyhM8ffQYBADgbAzQAnUn4RqcnNZhBo+JYjV5cu0uVJ2r9st9g9BkEAIQXkj0Tyc3NVW5urmw2W7BDCUuRERYNTevQYHlKYluPiy43xSIpyXqmDAsAAP7EDBomZNYK3OHMk+nU3MEMHAAQesz6/c2TPcANrpp57UWT1+86pFXbD7jV3EudPQBAIJHsAW5y1cwrSRMu6i5bneEyGWQGDQBAsJDsAT7SWDIIAECwUHoFAAAghPFkDwgQmnkBAMFAsgcEgCejeRnAAQDwJZpxAT/LLyrTjCWFbpdtKas8qRlLCpVfVObnyAAA4YBkD/AjW52hnLxirwoy5+QVy1ZHGUwAQPOQ7AF+tKGkwqtCzIbOPOHbUFLh+6AAAGGFZA/wowPVzZtxo7nbAwBAsgf4Uef4mKBuDwAAyR7gR4NTE5Vs9Txhs+jMqNzBqYm+DwoAEFZI9gA/ioywaF52urypmjcvO516ewCAZiPZA/xsTEayFk/OdPsJX7I1RosnZ1JnDwDgExRVBgJgTEayRqcnMYMGACDgSPaAAImMsGhoWodghwEACDM04wIAAIQwkj0AAIAQRrIHAAAQwkj2TCQ3N1fp6enKysoKdigAACBEWAzDYKZ1k6mqqpLValVlZaUSEhKCHQ4AAHCDWb+/ebIHAAAQwkj2AAAAQhjJHgAAQAgj2QMAAAhhJHsAAAAhjGQPAAAghJHsAQAAhDCSPQAAgBBGsgcAABDCSPYAAABCGMkeAABACCPZAwAACGEkewAAACGMZA8AACCEkewBAACEMJI9AACAEBYV7AAAnGGrM7ShpEI/VJ5QxbFTate2tY4cr//fimM1OnKiVoYhtW/bWomxZ15LjItWUkKMBqcmKjLCEuxDAQCYCMkeYAL5RWXKyStWWeXJZr1PsjVG87LTNSYj2UeRAQBaOppxgSDLLyrTjCWFzU70JKms8qRmLClUflGZDyIDAIQCnuwBQWSrM5STVyzDh+9pSHpo2VaVVpygmRcAQLIHBNOGkgqfPNE7V9VJm55csb3Bcpp5ASD80IwLBNGBat8neo2hmRcAwg/JHhBEneNjgrLfnLxi2ep82XgMADArkj0giAanJirZGtiEz9CZJ3wbSioCul8AQHCQ7PnYsmXLNH78ePXo0UNt27ZV//799fvf/161tbXBDg0mFBlh0bzsdAVjyMS7G7/Xe5v26k//+k7LN+/T+l3lPO0DgBBkMQyDv+4+NGTIEPXq1Uvjxo1Tly5dtG7dOj3xxBOaOHGi/vznP7v1HlVVVbJaraqsrFRCQoKfI4YZ+KrOXnNZY6I0Or2LhqZ1rFfQmdG8ANA0s35/k+z52MGDB9WpU6d6y5544gnNnTtXP/zwg7p06dLke5j1YoF/eTqDhrVNK7302S5VngjcU2NG8wKAa2b9/qb0io+dm+hJ0kUXXSRJ2r9/v1vJHsJTZIRFQ9M6eLRNzw5tNWNJoU/r9DWmrPKkpi8p1NShPdQjMdZpQursvzwZBIDgaXHJ3jfffKOVK1dq06ZN2rRpk7Zv3y6bzabHH39cjz32WJPbL126VLm5udq6datOnTql3r1765e//KVmz56tVq1a+SXmf/7zn2rdurXS0tL88v4IX2MykrV4cmbAm4D/vP57r7bjySAABF6LS/YWL16sRYsWebXtvffeq0WLFikqKkqjRo1SXFyc1qxZozlz5igvL08rV65UmzZtfBpvcXGxFi1apDvvvNNUj3QROsZkJGt0elKDJuCCHQe1fMv+YIdXj73O3+LJmSR8ABAgLW40bkZGhh544AG99dZb2r59u2655Ra3tvvggw+0aNEixcXF6YsvvtAnn3yi9957Tzt27NCAAQNUUFCguXPn1tvm9ddfl8ViafJn2bJlTvd56NAhjRs3Tr1799aCBQuafeyAK/Ym4PGZ3XX7iPM04aLumpjVI9hhuUSdPwAInBb3ZO+OO+6o93tEhHv56lNPPSVJevjhh5WZmelY3rFjR7300ksaMWKEXnzxRc2dO1dWq1WSNH78eA0ZMqTJ9+7WrVuDZdXV1br66qt16tQpffbZZ4qNjXUrTsBX7DX8gj3C91xn1/nztI8iAMBzLS7Z88a+ffu0ceNGSdLNN9/c4PXhw4crJSVFpaWlWrFihSZNmiRJslqtjsTPEzU1Nbruuuu0e/duFRQUqGvXrs07AMAL9hp+gRzA4YlATxUHAOGqxTXjemPz5s2SpMTERKWmpjpd5+KLL663rrdsNptuuukmbdy4UStWrFDfvn2b3KampkZVVVX1fgBfsA/gCPQsHe4I1lRxABBuwuLJXklJiSSpRw/XfZhSUlLqreutWbNm6YMPPtDjjz8um82mzz//3PFaenq600EaTz/9tHJycpq1X8AVVwM47KVR1u86pFXbDwS0Xl9s6wj9c8cBFe2rVMd4yrIAgD+FRbJXXV0tSY32m4uLi5OkZj9Vy8/PlyTNnTu3wYCPtWvXauTIkQ22eeSRR3Tfffc5fq+qqnIkn4AvNFbDb8JF3V0WdP6+4rj+vH6Pz+M5dqpOiz/7rt4yyrIAgH+ERbIXSLt37/Z4m+joaEVHR/s+GMBNjSWDQ9M6BKSOHwWbAcA/wiLZi4+PlyQdO3bM5TpHjx6VJGrhAedoqhnY1X8rjtXoxbWeT+dGwWYA8K2wSPZ69eolSSotLXW5jv01+7oA/subqdzW7yoPaD9A+5PBWy/pqSv7J/OkDwD+f2ExGnfQoEGSpPLycpcDML788ktJqleDL9Byc3OVnp6urKysoMUA+EqwSqu8tm6PJr3yuYY/s0b5RWVBiQEAzCQskr3u3bs7Eqi33367wesFBQUqLS1VdHS0xo4dG+jwHGbNmqXi4mJHTUCgJQt2aRX71GwkfADCXVg040rSo48+qvHjx2vBggW6+uqrHU/wysvLNXPmTEnSXXfd5VURZQANmWEGD0PSQ8u2qrTihBJjGfABIDxZDMMwY3F9lwoLCx3JmSTt2rVLhw4dUvfu3etNW7Z8+XIlJ9fvqH3PPffo+eefV6tWrXT55ZcrNjZWn376qY4cOaJhw4Zp1apVatOmTcCOxZWqqipZrVZVVlYyYAQtWn5RmWln8GgKAz4AeMqs398tLtn77LPPdNlllzW5XklJidPBFn/961+Vm5urLVu2qLa2VmlpaZo8ebJmz56t1q1b+yFiz5n1YgG8kV9UFpDSLf5gkbR4ciYJHwC3mPX7u8Ule+HArBcL4K2zizYfOlqjIydqtfPAUX3y1Y/BDq1RFklJ1hgVzBlFky6AJpn1+zts+uy1BLm5ucrNzZXNZgt2KIBPuSrdYvanfobODPTYUFLhcekZADALnuyZkFn/ZQD4g6up2pr6r7+mcnNm0U0/0XU/6db0igDCmlm/v3myByCovCnYbBeoqdyCXUYGAJojLOrsAQhNYzKSVTBnlN66/aeytmnl8/e36Myo3MGpiT5/bwAIFJI9AC1aZIRFw87vqGcmDJA/hlDMy05ncAaAFo1mXAAhYUxGshZPzvRZs641Jkqj07voWI1Nf/rXdxRdBtBiMUDDhMzawRNoCbwd8NGubWut33VIq7YfUOWJ2ib3Q9FlAOcy6/c3T/ZMhNIrQPN5O+Ajv6hM7xXuc3u2D/vcuxRdBmB2PNkzIbP+ywAIVbY6Q8OfWeNx8y9FlwGczazf3wzQABD2NpRUeNXP7+yiywBgViR7AMLegermDeho7vYA4E/02QMQ9ppbNHnt9h91oKpGibGNDwJhJC+AYCDZAxD2BqcmKtka43XJlg+2lumDrWVur89IXgCBRDMugLAXGWHRvOx0vxRldsY+kje/yP0EEQC8RbIHAPpvUeZka+Dmwc3JK5atjoIIAPyLZlwToc4eEFxjMpI1Oj2pQVHmLaWH9ebn3/t0X2eP5PWmLiAAuIs6eyZk1jo9QLj6cMs+3fPOFr+897ifdNWI8zsxgAMIAWb9/ubJHgA0obmjdRvzwZb9+mDL/nrLGMABwJfoswcATbCP1g0UBnAA8CWSPQBoQqBH69oxgAOAL9Bnz4TM2uYPhLv8ojLl5BV7XY/PG7cM6aGfpLRvtFgzRZsBczDr9zfJngmZ9WIBINnqDMdo3UNHa3TkRK0MQ2rftrXTGTQKdhzU8nP65PmbNSZKo9O7aGhaxyaTxIpjTR8DiSbgHrN+fzNAAwA8EBlh8ahUStd2bQKe7FWePK1lhfu0rHBfQPcreZZo+vK/JJuAayR7AOBHzZ2KraUJZqIpMZIZcIYBGgDgR8Ea3BGuGMkMNMSTPRNhBg0gNNmnYgv04I5wZUh6cNlWHT15WpUnamkCRthjgIYJmbWDJ4DmOXtwhz+nYoN3aAJGc5n1+5tkz4TMerEA8D1bnaHhz6zhiZ9JWCQtnpxJwgevmPX7mz57ABBE9OkzH4pZI9TwZM+EzPovAwD+E4yCzXDt7GLW9OmDu8z6/U2yZ0JmvVgA+JezPn3u1ppbv+uQVm0/oMoTtcE+jJBFnz40xazf3yR7JmTWiwWAuXmTLPpiBo1wSjTp04fGmPX7m2TPhMx6sQCAK815KunNfxPatNJv/1as6pOnA3qcFklJ1hgVzBlFky4aMOv3N3X2AADN5uk0cr4QHxOlGUsKFcgnFobOFG7eUFIR8OMFvMVoXABAi2QvVp1sjQn4vg9UM5AGLQdP9gAALdaYjGSNTk/yuAm5ucWsO8cHPsEEvEWyZyJMlwYAnvOmCXncoG5avf2AV6VuEmIiVbSvUvuPnKA0C1oEBmiYkFk7eAJAKMkvKvNpnz9Ks8Cs39/02QMAhCVf9/krqzypGUsKlV9U5pP3A3yFZlwAQNhy1eev4liNXly7y+PagYakh5ZtVWnFCY9qFnpb75AmZLiDZA8AENac9flbv6vc6yLRVSdtenLFdl+E5jaakNEYmnEBADhHSyutQhMyGsOTPQAAztESS6s4a0KmmRcSyR4AAA0MTk1UsjXGq9IsweSqCZlm3vBG6RUTMuvQbQAIJ74uzWIGU4f2UI/EWLfnIObJoGfM+v1NsmdCZr1YACDc5BeVKSevuMU94fM1a0yURqd30dC0jiSDjTDr9zfJngmZ9WIBgHBkqzN8VpolVNFMfIZZv7/pswcAQCNcTceWktg25Jp5vWUfDbx4cmbYJ3xmROkVAAC84OsZOEJBTl6xbHWkv2ZDM64JmfUxMACgIWfNvO7OnOHNDBpmb0L+5U9T1K5ta1l05onokPM6hE1/PrN+f5PsmZBZLxYAgDm0pJHC7dq20oLrB4RF865Zv79pxjWR3NxcpaenKysrK9ihAABMrCU1IR85XqvpzO4RVDzZMyGz/ssAAGAurpqQv684rj+v3xPs8OpJtsaoYM6okG7SNev3N6NxAQBooVyNFJakoWkdTFUjsKzypDaUVLiMF/5DsgcAQAgak5Gs0elJXg8eWb/rkFZtP+DTgSAHqs2ReIYbkj0AAEJUY0/+mjLhou5Om4m3lB7Wm59/79V7do43fx/DUESyBwAAnHKWLI4b1E2rtx/wuHk42XpmWjUEHqNxAQCA2yIjLJqXnS5Ph1nMy04P6cEZZkayBwAAPOJJ6Zf2bVvpD0yjFlSUXjEhsw7dBgDgbM769HkyG0hj/02Mi1ZSwpmm35byRNCs39/02QMAAF5xNQAkv6jMZ2Vfkq0xmpedzpPBZqAZFwAA+Ix9Kjdf1fcrqzypGczA0SwkewAAwCdsdYZy8or9MmdvTl6xbHX0PPMGyR4AAPCJDSUVfpmxw9B/Z+CA50j2AACAT/h7hgxm4PAOyR4AAPAJf8+QwQwc3iHZAwAAPjE4NdGt2nuesogZOJqDZA8AAPiEt7NruIMZOLxHnT0AAOAz9tk1fFVnzxoTpdHpXXSsxqY//eu7kC7K7C/MoGFCZq3ADQCAu5zNruHuDBrrdx3Squ0HVHmittlxBLIos1m/v3myBwAAfM7V7BpNyS8q03uF+3xWq89elHlxGM/PS7IHAABMwV9FmQ1JDy3bqtKKE465esOpmZdkz0Ryc3OVm5srm80W7FAAAAg4fxVllqSqkzY9uWJ7g+XhMPcuo3FNZNasWSouLtbGjRuDHQoAAAEXjKLJ4TD3LskeAAAwhWAWTQ7luXdJ9gAAgCn4qyhzU0J97l2SPQAAYAr+LMrsjlCde5dkDwAAmIa9KHMwnvCF6ty7jMYFAACmMiYjWaPTkzwuylxxrEYvrt3lcTFmi6SkEJ57l2QPAACYjrdFmVMS22rGkkKPa/WF8ty7NOMCAICQ4WkzcLI1JuRn1+DJHgAACCnuNAMzgwYAAEAL5m0zcCiiGRcAACCEkewBAACEMJI9AACAEEayBwAAEMJI9gAAAEIYyR4AAEAII9kDAAAIYSR7AAAAIYxkDwAAIIQxg4YJGcaZ6ZurqqqCHAkAAHCX/Xvb/j1uFiR7JlRdXS1JSklJCXIkAADAU9XV1bJarcEOw8FimC39hOrq6rR//37Fx8fLYvHt5MxVVVVKSUlRaWmpEhISfPrecB/nwRw4D+bAeTAHzkPzGYah6upqde3aVRER5ukpx5M9E4qIiFD37t39uo+EhARuZhPgPJgD58EcOA/mwHloHjM90bMzT9oJAAAAnyPZAwAACGEke2EmOjpa8+bNU3R0dLBDCWucB3PgPJgD58EcOA+hiwEaAAAAIYwnewAAACGMZA8AACCEkewBAACEMJK9MLF06VKNHDlS7du3V2xsrAYOHKiFCxeqtrY22KGFjGnTpslisTT6c/LkSafbbtq0STfeeKO6dOmimJgYpaam6u6779aBAwcCfBQtwzfffKMXXnhB06ZN04ABAxQVFSWLxaInnniiyW1Xr16tsWPHqmPHjmrTpo0uuOAC/eY3v9HRo0cb3W7nzp2aNm2aunfvrujoaHXv3l3Tpk3Td99956vDanG8OQ/z589v8j75+uuvXW7PeaivtrZWn376qR588EFlZWWpXbt2atWqlZKSknTttdfq73//e6Pbcz+ECQMh75577jEkGVFRUcaVV15pXH/99Ua7du0MScbw4cON48ePBzvEkDB16lRDkjFs2DBj6tSpTn9OnTrVYLulS5caUVFRhiQjKyvLmDhxonHeeecZkowuXboYO3bsCMLRmJv9mj735/HHH290u2effdaQZFgsFuPSSy81brzxRiMpKcmQZPTt29c4ePCg0+0KCgqMtm3bGpKM/v37G7/4xS+M/v37G5KM2NhYY/369f44TNPz5jzMmzfPkGQMHDjQ5X2yf/9+p9tyHhpatWqV43NPSkoyrrnmGmPixIlGRkaGY/mdd95p1NXVNdiW+yF8kOyFuOXLlxuSjLi4OGPTpk2O5QcPHjQGDBhgSDLuv//+IEYYOuzJ3muvveb2Nvv27XP80Xz55Zcdy0+fPm1MnjzZkQA6+0Mdzl555RXjgQceMN566y1j+/btxi233NJkklFYWGhYLBYjMjLSWLFihWP5sWPHjMsvv9yQZEyYMKHBdseOHTO6du1qSDIeeeSReq898sgjhiQjJSUlLP/R5M15sCd78+bN82hfnAfnPv30U2PChAnGP//5zwavvfPOO0ZkZKQhyfjzn/9c7zXuh/BCshfisrKyDEnGE0880eC1f/3rX4YkIzo62jhy5EgQogst3iR7Dz74oCHJuOKKKxq8Vl1dbVitVkOSkZ+f78NIQ4/9s28sybjxxhsNScYdd9zR4LXdu3cbERERhiRj+/bt9V7Lzc01JBl9+vQxbDZbvddsNpvRp08fQ5Lxhz/8wTcH04K5cx68TfY4D965/fbbDUnG5ZdfXm8590N4oc9eCNu3b582btwoSbr55psbvD58+HClpKSopqZGK1asCHR4kLR8+XJJzs9PXFycrr32WknS+++/H9C4Qs2pU6ccfZecfdY9e/bUsGHDJP33nNjZf7/pppsaTGweERGhX/ziF5I4R/7GefDOoEGDJEmlpaWOZdwP4Scq2AHAfzZv3ixJSkxMVGpqqtN1Lr74YpWWlmrz5s2aNGlSIMMLWWvXrtV//vMfVVdXq0OHDho8eLDGjh3boCp9dXW1du7cKenMeXDm4osv1ptvvuk4l/DOt99+q+PHj0tq/LP+17/+1eCztv/e2HZnrwf3FBYW6uGHH1ZFRYWsVqsGDRqk7OxsxcfHO12f8+CdHTt2SJKSk5Mdy7gfwg/JXggrKSmRJPXo0cPlOikpKfXWRfO98cYbDZYlJyfrf//3fzVmzBjHst27dzv+39U54vz4hv3za9eunctkwtlnXV1drfLycklNn6ODBw/q2LFjio2N9VncoSwvL095eXn1llmtVj3//POaMmVKveWcB+/88MMPev311yVJEyZMcCznfgg/NOOGsOrqaklq9GaLi4uTJFVVVQUkplA2cOBALVq0SEVFRaqqqtKPP/6olStX6pJLLlFZWZmuvfZaffbZZ4717edHcn2OOD++4e294Mk5OndbOJeWlqannnpKmzdvVkVFhSoqKlRQUKCf//znqqys1NSpU/XWW2/V24bz4LnTp09r8uTJqqys1IABA/SrX/3K8Rr3Q/jhyR7gI7Nnz673e3x8vEaPHq0rrrhC48eP14cffqh7771XW7ZsCU6AgAnccsstDZYNGzZMeXl5+vWvf60XXnhBs2fP1o033qjWrVsHIcLQMH36dH366afq0KGDli1bxmcZ5niyF8Lsj+ePHTvmch174cyEhISAxBSOLBaLcnJyJElbt251dJQ+u/nE1Tni/PiGt/eCJ+fo3G3hufnz5ysyMlIHDx7UF1984VjOefDMPffcoz/96U9q3769Vq1apT59+tR7nfsh/JDshbBevXpJqj8K61z21+zrwj/69evn+P+9e/dKOjPize777793uh3nxzfsn9+RI0fqNUWdzdlnHR8fr8TERElNn6OOHTvSP6mZEhMT1blzZ0n/vU8kzoMn7r//fj3//PNq166dVq5c6RiNezbuh/BDshfC7Dd5eXm5yw7+X375pSQpMzMzYHGFI3unZum//zpOSEhQ7969Jf33PJyL8+Mbffv2Vdu2bSV5/lnbf+cc+Z/NZlNlZaUkNRg4wHlo2kMPPaRnn31WVqtVK1eudDlilvsh/JDshbDu3bsrKytLkvT22283eL2goEClpaWKjo7W2LFjAx1eWHnnnXcknUnw+vbt61g+fvx4Sc7Pz9GjRx2jFa+//voARBm6WrdurWuuuUaS8896z549WrdunaT/nhM7++/vvPOO6urq6r1WV1end999VxLnyBc++ugjHT9+XBaLpUGiwnlo3MMPP6zf/e53slqtWrVqleNvvzPcD2Eo2FWd4V+upks7dOgQ06X50ObNm40PP/zQqK2trbfcZrMZr776qhETE2NIMh577LF6r589Xdof//hHx/LTp087pp5iurSmuTNzw6ZNmxzTQ3388ceO5Z5MD/Xoo4/We+3RRx81JBndu3dneiij6fOwZ88e48033zROnDjR4LXly5cbiYmJhiRj8uTJDV7nPLj2m9/8xpBktGvXztiwYYNb23A/hBeLYRhGwDNMBNQ999yj559/Xq1atdLll1+u2NhYffrppzpy5IiGDRumVatWqU2bNsEOs0X74IMPNH78eLVv316ZmZnq0qWLjhw5oqKiIkfflkmTJumNN95QVFT9QfBLly7VpEmTZLPZ9NOf/lS9evXSxo0b9d1336lLly4qKChwNPfijMLCQs2cOdPx+65du3To0CF1795d3bp1cyxfvnx5vWKyzz33nO677z5ZLBb97Gc/U+fOnfWvf/1LZWVl6tu3rwoKCtSxY8cG+/v3v/+tK6+8UsePH1dGRoYyMjJUVFSkoqIixcbGavXq1RoyZIh/D9qEPD0PW7Zs0aBBgxQXF6dBgwapW7duOnHihIqLix3Ffy+77DJ99NFH9Up42HEeGvroo4903XXXSTpT0Lh///5O1+vYsaP+7//9v/WWcT+EkWBnmwiMd99917j00kuNhIQEo02bNkZGRoaxYMECo6amJtihhYTvvvvOuPfee43hw4cb3bp1M2JiYozo6GijR48exg033GD8/e9/b3T7L7/80rj++uuNTp06Ga1btzZ69uxpzJo1y/jhhx8CdAQty9q1aw1JTf6UlJQ02HbVqlXGmDFjjMTERCM6Oto4//zzjUceecSoqqpqdJ87duwwpkyZYnTt2tVo1aqV0bVrV2PKlCnGzp07/XSU5ufpeTh06JAxZ84cY9SoUUaPHj2M2NhYo1WrVkZycrLx85//3Hj77bcbzLd6Ls5Dfa+99ppb56Bnz55Ot+d+CA882QMAAAhhDNAAAAAIYSR7AAAAIYxkDwAAIISR7AEAAIQwkj0AAIAQRrIHAAAQwkj2AAAAQhjJHgAAQAgj2QOAADt69KiioqKUkJAgm80W7HAAhDiSPQAIsC+++EI2m01DhgxRZGRksMMBEOKYLg0AAqyyslJlZWVq3769unTpEuxwAIQ4kj0AAIAQRjMuAASQYRhq166dIiIiVFlZGexwAIQBkj0ACKBdu3apsrJSaWlpslqtwQ4HQBgg2QOAANq0aZMkKTMzM8iRAAgXJHsAEECFhYWSpIsuuijIkQAIFyR7ABBA9id7JHsAAoXRuAAQQImJiTp8+LAqKirUvn37YIcDIAzwZA8AAqSkpESHDx9WamoqiR6AgCHZA4AAoQkXQDCQ7AFAgDASF0AwkOwBQIAwEhdAMDBAAwACpGPHjiovL9ehQ4fUoUOHYIcDIEzwZA8AAmDPnj0qLy9Xjx49SPQABBRP9gAAAEIYT/YAAABCGMkeAABACCPZAwAACGEkewAAACGMZA8AACCEkewBAACEMJI9AACAEEayBwAAEMJI9gAAAEIYyR4AAEAII9kDAAAIYSR7AAAAIYxkDwAAIIT9f419AqouA7QGAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "coeffs = [term.coefficient for term in mol_hamiltonian]\n",
    "plt.semilogy(np.sort(np.abs(coeffs))[::-1], \"o\")\n",
    "plt.ylabel(r\"$|\\alpha_i|$\", fontsize=16)\n",
    "plt.xlabel(r\"$i$\", fontsize=16)\n",
    "plt.tick_params(axis=\"both\", labelsize=16)\n",
    "plt.title(\n",
    "    r\"Sorted coefficients size for Hamiltonian $H = \\sum \\alpha_i P_i$   ($P_i$ Pauli string)\"\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21828537-6d50-4807-bf98-73f0a7434418",
   "metadata": {},
   "source": [
    "We define some threshold and trim the Hamiltonian accordingly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e84941d6-18e3-430b-b3e0-570c2858df50",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T20:44:27.521414Z",
     "start_time": "2025-03-04T20:44:27.515198Z"
    }
   },
   "outputs": [],
   "source": [
    "def trim_hamiltonian(hamiltonian, threshold):\n",
    "    return [\n",
    "        PauliTerm(pauli=term.pauli, coefficient=term.coefficient)\n",
    "        for term in hamiltonian\n",
    "        if np.abs(term.coefficient) > threshold\n",
    "    ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a3e74354-df3d-481e-a00d-e36b69279d5a",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T20:44:27.856262Z",
     "start_time": "2025-03-04T20:44:27.848941Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Length of trimmed Hamiltonian: 49\n"
     ]
    }
   ],
   "source": [
    "THRESHOLD = 0.03\n",
    "\n",
    "trimmed_mol_hamiltonian = trim_hamiltonian(mol_hamiltonian, THRESHOLD)\n",
    "print(f\"Length of trimmed Hamiltonian: {len(trimmed_mol_hamiltonian)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6ef93ba-a682-474d-befd-0b0b5a691de5",
   "metadata": {},
   "source": [
    "### Normalizing the Hamiltonian for QPE\n",
    "\n",
    "Since we are working with QPE, the ground state energy is inferred as a phase, thus, we shall normalize the Hamiltonian such that its eigenvalues are in $\\left[-\\frac{1}{2},\\frac{1}{2}\\right)$. This is done by finding a bound on the maximal absolute value of eigenvalues $\\tilde{\\lambda}_{\\max}$, and normalizing the Hamiltonian by $2\\tilde{\\lambda}_{\\max}$. A simple bound is given by the sum of Pauli coefficients of the Hamiltonian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "0639c8de-a611-4fd3-bb8a-bcd4da64485e",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T20:44:28.271155Z",
     "start_time": "2025-03-04T20:44:28.265440Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The normalization value of the Hamiltonian is 17.61546220154492\n"
     ]
    }
   ],
   "source": [
    "def normalize_hamiltonian(hamiltonian):\n",
    "    approx_lambda_max = sum(np.abs(term.coefficient) for term in hamiltonian)\n",
    "    normalization = 2 * approx_lambda_max\n",
    "    normalized_mol_hamiltonian = [\n",
    "        PauliTerm(pauli=term.pauli, coefficient=term.coefficient / (normalization))\n",
    "        for term in hamiltonian\n",
    "    ]\n",
    "\n",
    "    return normalization, normalized_mol_hamiltonian\n",
    "\n",
    "\n",
    "normalization, normalized_mol_hamiltonian = normalize_hamiltonian(\n",
    "    trimmed_mol_hamiltonian\n",
    ")\n",
    "print(f\"The normalization value of the Hamiltonian is {normalization}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "817d2e3b-ec15-4998-847c-f7942b37dcc1",
   "metadata": {},
   "source": [
    "In our QPE models we will define the `phase` variable as a fixed point quantum number in $[-1/2,1/2)$: `QNum[QPE_SIZE,SIGNED,QPE_SIZE]`. Thus, we need to postprocess the algorithm result by multiplying back the normalization factor (see Measurement and Anlysis section below)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf78c5bb-240f-429e-aa47-d7c9720d90c1",
   "metadata": {},
   "source": [
    "## Designing the Quantum Model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00d22272-d61f-41e6-bcbd-a2801db7ee1f",
   "metadata": {},
   "source": [
    "### Defining a powered Hamiltonian simulation\n",
    "\n",
    "We will now create a quantum model of the QPE algorithm using the Classiq platform. In particular, we will use the open library  [qpe_flexible](https://github.com/Classiq/classiq-library/blob/main/functions/qmod_library_reference/classiq_open_library/qpe/qpe.ipynb) function (see [notebook](https://github.com/Classiq/classiq-library/blob/main/tutorials/advanced_tutorials/high_level_modeling_flexible_qpe/high_level_modeling_flexible_qpe.ipynb) as well).\n",
    "\n",
    "We need to approximate the Hamiltonian simulation $e^{2\\pi i H}$, for this we will use Classiq built-in implementation for [Suzuki Trotter formulas](https://github.com/Classiq/classiq-library/blob/main/functions/qmod_library_reference/qmod_core_library/hamiltonian_evolution/suzuki_trotter/suzuki_trotter.ipynb). For a given Suzuki Trotter order $o$, we can specify a repetitions parameter $r$, which controls the level of approximation. The literature provides some lower bounds for $r$ as a function of the operator error $\\epsilon$ (defined by the dimond norm [[3](#Dimond)]). For example, Eq. (14) in Ref. [[2](#Bounds)]  states that Suzuki Trotter formula of order 2 approximates $e^{i \\sum \\alpha_m P_m t}$ up to an error $\\epsilon$, given $r$ repetitions that satisfies:\n",
    "$$\n",
    "r \\leq \\left(\\frac{2^5\\gamma_2}{3\\epsilon}\\right)^{1/2}  t^{3/2} ,\n",
    "\\tag{1}\n",
    "$$\n",
    "where $\\gamma_2 \\equiv \\sum_{l,m,n} |\\alpha_m\\alpha_n\\alpha_l| \\left |\\left[P_l,\\left[P_m, P_n\\right]\\right]\\right|_\\infty$. **In particular, we can see that the number of repetitions grows as $t^{3/2}$**.\n",
    "\n",
    "\n",
    "In QPE, we apply powered Hamiltonian simulation:\n",
    "$$\n",
    "\\left(e^{2\\pi i H}\\right)^ p = e^{2p \\pi i H},\n",
    "\\tag{2}\n",
    "$$\n",
    "and each power should be approximated with Suzuki Trotter with an appropriate order and repetitions parameters, keeping the same error per QPE iteration. We can thus use the bound above to define a powered Suzuki Trotter `qfunc` for our specific molecule.\n",
    "\n",
    "First, we first define some classical auxiliary functions that helps us to evaluate the right-hand-side in Eq. (1):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "740ccdac-65a1-4b70-bbee-a340950683ab",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T20:44:29.178026Z",
     "start_time": "2025-03-04T20:44:29.176709Z"
    }
   },
   "outputs": [],
   "source": [
    "import itertools\n",
    "\n",
    "\n",
    "def multiply_single_qubit(pa, pb):\n",
    "    \"\"\"\n",
    "    Multiply two single-qubit Pauli operators pa*pb.\n",
    "    Returns (phase, pc) where:\n",
    "       - phase in {1, -1, 1j, -1j}\n",
    "       - pc is the resulting single-qubit Pauli.\n",
    "    \"\"\"\n",
    "    table = {\n",
    "        (Pauli.I, Pauli.I): (1, Pauli.I),\n",
    "        (Pauli.I, Pauli.X): (1, Pauli.X),\n",
    "        (Pauli.I, Pauli.Y): (1, Pauli.Y),\n",
    "        (Pauli.I, Pauli.Z): (1, Pauli.Z),\n",
    "        (Pauli.X, Pauli.I): (1, Pauli.X),\n",
    "        (Pauli.X, Pauli.X): (1, Pauli.I),\n",
    "        (Pauli.X, Pauli.Y): (1j, Pauli.Z),  # X*Y = iZ\n",
    "        (Pauli.X, Pauli.Z): (-1j, Pauli.Y),  # X*Z = -iY\n",
    "        (Pauli.Y, Pauli.I): (1, Pauli.Y),\n",
    "        (Pauli.Y, Pauli.X): (-1j, Pauli.Z),  # Y*X = -iZ\n",
    "        (Pauli.Y, Pauli.Y): (1, Pauli.I),\n",
    "        (Pauli.Y, Pauli.Z): (1j, Pauli.X),  # Y*Z = iX\n",
    "        (Pauli.Z, Pauli.I): (1, Pauli.Z),\n",
    "        (Pauli.Z, Pauli.X): (1j, Pauli.Y),  # Z*X = iY\n",
    "        (Pauli.Z, Pauli.Y): (-1j, Pauli.X),  # Z*Y = -iX\n",
    "        (Pauli.Z, Pauli.Z): (1, Pauli.I),\n",
    "    }\n",
    "\n",
    "    return table[(pa, pb)]\n",
    "\n",
    "\n",
    "def multiply_pauli_lists(pA_list, pB_list):\n",
    "    \"\"\"\n",
    "    Multiply two multi-qubit Pauli lists A and B.\n",
    "    Returns PauliTerm(pauli, coefficient), where:\n",
    "       - coefficient in {1, -1, 1j, -1j}\n",
    "       - pauli is a Pauli list representing AB.\n",
    "    \"\"\"\n",
    "    if len(pA_list) != len(pB_list):\n",
    "        raise ValueError(\"Pauli lists must have the same length.\")\n",
    "\n",
    "    phase_total = 1 + 0j  # start with complex 1\n",
    "    product = []\n",
    "\n",
    "    for pA, pB in zip(pA_list, pB_list):\n",
    "        phase_local, pC = multiply_single_qubit(pA, pB)\n",
    "        phase_total *= phase_local\n",
    "        product.append(pC)\n",
    "\n",
    "    return PauliTerm(pauli=product, coefficient=phase_total)\n",
    "\n",
    "\n",
    "def commutator_infinity_norm(pA_list, pB_list):\n",
    "    \"\"\"\n",
    "    Compute the 'infinity-norm of [A,B]' . The commutation relation between two Pauli lists\n",
    "    can be inferred by counting the number of qubit positions i where A[i] != B[i],\n",
    "    ignoring position with  A[i] or B[i] equal to Pauli.I.\n",
    "    Then:\n",
    "      - If the count is even, A and B commute => norm = 0\n",
    "      - If the count is odd, A and B anticommute => norm = 2\n",
    "\n",
    "    \"\"\"\n",
    "    if len(pA_list) != len(pB_list):\n",
    "        raise ValueError(\"Pauli strings must have the same length.\")\n",
    "\n",
    "    # Count positions where A and B are both in {X,Y,Z}, differ, and none is I\n",
    "    difference_count = 0\n",
    "    for pA, pB in zip(pA_list, pB_list):\n",
    "        if pA != Pauli.I and pB != Pauli.I and pA != pB:\n",
    "            difference_count += 1\n",
    "\n",
    "    return 2 * (difference_count % 2)\n",
    "\n",
    "\n",
    "def calculate_gamma_2(hamiltonian):\n",
    "    \"\"\"\n",
    "    Compute the $\\gamma_2$ value appearing in the bound for Suzuki Trotter of order 2\n",
    "    \"\"\"\n",
    "    gamma_2 = 0\n",
    "    for triplet in itertools.combinations(range(len(hamiltonian)), 3):\n",
    "        terms = [hamiltonian[index] for index in triplet]\n",
    "        factor = np.abs(\n",
    "            (terms[0].coefficient) * (terms[1].coefficient) * (terms[2].coefficient)\n",
    "        )\n",
    "        inner_commutator = commutator_infinity_norm(terms[1].pauli, terms[0].pauli)\n",
    "        if inner_commutator != 0:\n",
    "            outer_commutator = 2 * commutator_infinity_norm(\n",
    "                terms[2].pauli,\n",
    "                multiply_pauli_lists(terms[0].pauli, terms[1].pauli).pauli,\n",
    "            )\n",
    "            gamma_2 += factor * outer_commutator\n",
    "\n",
    "    return gamma_2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13fd6fa0-b39e-491b-bf84-7d4996f122a1",
   "metadata": {},
   "source": [
    "In QPE, the powers of the Hamiltonian simulation grows exponentially with the phase variable size. Let us examine the number of repetitions needed per QPE iteration, according to the bound above for QPE of size 7."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "204977f9-3242-41d7-a6e2-8862aed2b031",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T20:44:29.664418Z",
     "start_time": "2025-03-04T20:44:29.662393Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The powers of the Hamiltonian simulation along a QPE routine of size 7: \n",
      "[ 1  2  4  8 16 32 64]\n"
     ]
    }
   ],
   "source": [
    "QPE_SIZE = 7\n",
    "qpe_powers = 2 ** np.arange(QPE_SIZE)\n",
    "print(\n",
    "    f\"\"\"The powers of the Hamiltonian simulation along a QPE routine of size {QPE_SIZE}: \n",
    "{qpe_powers}\"\"\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96671a1c-286b-4e20-a3ed-8b42a84170c2",
   "metadata": {},
   "source": [
    "These powers enter as an evolution coefficient for the Hamiltonian simulation (see Eq. (2) above). Using the theoretical bound we find"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "210fa2d5-2e1d-450e-9c20-bc3781f07b9d",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T20:44:30.203858Z",
     "start_time": "2025-03-04T20:44:30.104505Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The theoretical bounds for the repetitions for QPE size 7, keeping an error 0.1 per QPE iteration are:\n",
      "[   4.   11.   30.   85.  240.  677. 1915.]\n"
     ]
    }
   ],
   "source": [
    "EPS = 0.1\n",
    "gamma_2_LiH = calculate_gamma_2(normalized_mol_hamiltonian)\n",
    "theoretical_r0 = np.sqrt(2**5 * gamma_2_LiH / (3 * EPS)) * (2 * np.pi) ** (3 / 2)\n",
    "\n",
    "print(\n",
    "    f\"\"\"The theoretical bounds for the repetitions for QPE size {QPE_SIZE}, keeping an error {EPS} per QPE iteration are:\n",
    "{np.ceil(theoretical_r0*qpe_powers**(3/2))}\"\"\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f27cd6b0-ac6f-4f87-8fcd-4ea8495da3f9",
   "metadata": {},
   "source": [
    "We note that applying a naive QPE, i.e., assuming a single unitary approximated with Suzuki Trotter, $e^{iHt} \\approx {\\rm ST}(H, o, r ,t)$, and simply taking its powers, gives $\\left(e^{iHt}\\right)^p \\approx \\left({\\rm ST}(H, o, r ,t)\\right)^{p} = {\\rm ST}(H, o, pr ,pt)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "bb103796-9d0a-4fac-b281-8983d7e9ebd3",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T20:44:30.533308Z",
     "start_time": "2025-03-04T20:44:30.524656Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The repetitions for QPE size 7, taking a naive QPE, per QPE iteration:\n",
      "[  4.   8.  15.  30.  60. 120. 240.]\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    f\"\"\"The repetitions for QPE size {QPE_SIZE}, taking a naive QPE, per QPE iteration:\n",
    "{np.ceil(theoretical_r0*qpe_powers)}\"\"\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd52c197-23cc-4786-b937-9f0420aa9361",
   "metadata": {},
   "source": [
    "While this naive QPE results in a shallower circuit, compared to taking repetitions according to the theoretical bounds (due to smaller values of repetitions), it is unclear whether it keeps the same operator error per phase bit.\n",
    "\n",
    "In practice, the bounds given in the literature are quite loose. In this tutorial we thus take a more experimental approach. We assume that the scaling of the bound with the evolution time $t$ is similar to Eq. (1), but take a smaller prefactor. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "8a0424ec-2f10-4ccb-a84e-c5ccd255b727",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T20:44:31.191893Z",
     "start_time": "2025-03-04T20:44:31.190067Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The experimental repetitions for QPE size 7, per QPE iteration are:\n",
      "[ 1.  1.  1.  2.  4. 10. 26.]\n"
     ]
    }
   ],
   "source": [
    "experimental_r0 = 0.05\n",
    "print(\n",
    "    f\"\"\"The experimental repetitions for QPE size {QPE_SIZE}, per QPE iteration are:\n",
    "{np.ceil(experimental_r0*qpe_powers**(3/2))}\"\"\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6774add1-08df-4b43-b6b5-8b9d5d9d19e1",
   "metadata": {},
   "source": [
    "Now, we use this approach to define our powered Suzuki Trotter function, for the specific Hamiltonian at hand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "31b0c72b-9fcb-4403-a19d-9181cc716d85",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T20:44:31.827222Z",
     "start_time": "2025-03-04T20:44:31.815368Z"
    }
   },
   "outputs": [],
   "source": [
    "from classiq.qmod.symbolic import ceiling as ceiling_qmod, pi\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def powered_st2_for_LiH(p: CInt, state: QArray[QBit]):\n",
    "    suzuki_trotter(\n",
    "        pauli_operator=normalized_mol_hamiltonian,\n",
    "        evolution_coefficient=-2 * pi * p,\n",
    "        order=2,\n",
    "        repetitions=ceiling_qmod(experimental_r0 * p ** (3 / 2)),\n",
    "        qbv=state,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a17ed790-e728-4787-bbda-5ba671d2c124",
   "metadata": {},
   "source": [
    "### Defining and synthesizing the phase estimation model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "346b55fc-c89b-4dab-8f04-558d81bdce5f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T20:45:38.912360Z",
     "start_time": "2025-03-04T20:44:32.740923Z"
    }
   },
   "outputs": [],
   "source": [
    "@qfunc\n",
    "def main(\n",
    "    state: Output[QArray[QBit]], phase: Output[QNum[QPE_SIZE, SIGNED, QPE_SIZE]]\n",
    ") -> None:\n",
    "\n",
    "    allocate(problem_size, state)\n",
    "    molecule_hartree_fock(molecule_problem_to_qmod(gs_problem), state)\n",
    "    allocate(QPE_SIZE, phase)\n",
    "    qpe_flexible(lambda p: powered_st2_for_LiH(p, state), phase)\n",
    "\n",
    "\n",
    "qmod = create_model(\n",
    "    main,\n",
    "    preferences=Preferences(timeout_seconds=600),\n",
    "    out_file=\"qpe_for_molecules\",\n",
    ")\n",
    "qprog = synthesize(qmod)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7e16994-cca8-4b2b-b65d-79868780f749",
   "metadata": {},
   "source": [
    "## Measurment and Analysis\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "669c971b-82e1-4553-aa60-04c281248e5f",
   "metadata": {},
   "source": [
    "First we execute on the default simulator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "97391a32-334b-4c85-922a-767e65eee836",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T20:45:51.048568Z",
     "start_time": "2025-03-04T20:45:38.913523Z"
    }
   },
   "outputs": [],
   "source": [
    "res = execute(qprog).result_value()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eacde622-296a-4c80-8369-8af68648c973",
   "metadata": {},
   "source": [
    "Next, we draw a histogram for the energies, by taking the output of the `phase` variable and multypling back the normalization facor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "7bffbb3e-f3ba-44f2-b8f8-4a0432cb90e7",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T20:45:51.135994Z",
     "start_time": "2025-03-04T20:45:51.049352Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Energy with maximal probability: -1.1009663875965574 Ha\n",
      "Precision: 0.13762079844956968 Ha\n",
      "Classical solution:, -1.0789776863285236 Ha\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkoAAAHTCAYAAADYojePAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABNKklEQVR4nO3deVxUZeP///coOCAI7miigksuH/eFNLAMbdE2NbVQ79K7NDUj7W7RSkut9PZuw7K0XMsls3Ircrk1LVxR0TTL8hYUd0VZNEWR8/vDH/N1ggMzDMuAr+fjMY+HnHOd61znDMy8vc51rmMxDMMQAAAAsilT3A0AAABwVwQlAAAAEwQlAAAAEwQlAAAAEwQlAAAAEwQlAAAAEwQlAAAAEwQlAAAAEwQlAAAAEwQlAMXujTfekMViKe5mlFirVq1Sq1at5OXlJYvFouTk5OJuElBqEJSAQjB37lxZLBbT19atW4u7iYXKYrFoxIgROa7LOjc7duxwaR/Hjx/XG2+8od27d7tUT0mXlJSkvn37ytvbW9OmTdMXX3whHx+f4m5WnjZt2qSePXsqICBAVqtVQUFBGjp0qBITE7OVzQrSWa/y5curadOmeu2115Sammord7P/3aFweBR3A4DSbMKECQoODs62vEGDBsXQGvf12muvafTo0U5tc/z4cY0fP15BQUFq1apV4TSsBIiNjVVaWpomTpyorl27FndzHPLhhx/queeeU7169fTss8+qZs2a+u233zRz5kwtXrxYP/zwgzp06JBtu08++US+vr66cOGC1qxZo7feekvr16/Xpk2b7Hok+btDQSIoAYWoW7duateuXXE3QxcvXnTrXgYPDw95eJSsj6OMjAxlZmaqXLlyxdqO06dPS5IqVqyYZ9m//vpL5cuXL+QW5W7Tpk0aOXKkwsLCtGrVKrv2DBs2TKGhoXrkkUf066+/Zjum3r17q2rVqpKkoUOH6pFHHtG3336rrVu3qmPHjrZy7vJ3h9KBS29AMUpISJDFYtE777yjTz/9VPXr15fValX79u0VGxubrfzvv/+u3r17q3LlyvLy8lK7du20YsUKuzJZlx82btyo4cOHq3r16goMDLStnzZtmurVqydvb2+FhITo559/VufOndW5c2dJ0oULF+Tj46Pnnnsu2/6PHj2qsmXLatKkSQV6HnIao7R27VqFhYWpYsWK8vX1VaNGjfTKK69IkjZs2KD27dtLkgYNGmS7tDJ37lzb9kuWLFHbtm3l7e2tqlWrasCAATp27Fi2fS9ZskRNmzaVl5eXmjVrpqVLl2rgwIEKCgqylbnxffrggw9s79P+/ft15coVjRs3Tm3btpW/v798fHzUqVMn/fjjj3b7ubGOrPegfPnyuueee5SYmCjDMDRx4kQFBgbK29tbDz/8sM6dO5freevcubOeeOIJSVL79u1lsVg0cOBA27pmzZpp586duuOOO1S+fHnb+Tt9+rSefPJJBQQEyMvLSy1bttS8efMKvb2SNHHiRFksFs2bNy9baKtfv76mTJmi48eP69NPP82zrvDwcElSfHx8nmWB/CpZ/4UDSpiUlBSdPXvWbpnFYlGVKlXsli1cuFBpaWl6+umnZbFYNGXKFPXq1UuHDh2Sp6enJOnXX39VaGioatWqpdGjR8vHx0dfffWVevTooW+++UY9e/a0q3P48OGqVq2axo0bp4sXL0q6fulixIgR6tSpk0aNGqWEhAT16NFDlSpVsoUpX19f9ezZU4sXL9Z7772nsmXL2upctGiRDMNQ//798zz2y5cvZzt26XoQy8uvv/6qBx54QC1atNCECRNktVp18OBBbdq0SZLUpEkTTZgwQePGjdOQIUPUqVMnSdLtt98u6XpYHDRokNq3b69Jkybp1KlTioqK0qZNmxQXF2frqfj+++/16KOPqnnz5po0aZLOnz+vJ598UrVq1cqxXXPmzNHly5c1ZMgQWa1WVa5cWampqZo5c6YiIiI0ePBgpaWladasWbr33nu1ffv2bJcFFyxYoCtXrujZZ5/VuXPnNGXKFPXt21fh4eHasGGDXn75ZR08eFAffvihXnjhBc2ePdv0PL366qtq1KiRPv30U9vlpvr169vWJyUlqVu3bnrsscc0YMAABQQE6NKlS+rcubMOHjyoESNGKDg4WEuWLNHAgQOVnJycLSAXZHv/+usvrVu3Tp06dcrx0pgkPfrooxoyZIhWrlypl156ybQuSfrf//4nSdn+nhz9uwMcYgAocHPmzDEk5fiyWq22cvHx8YYko0qVKsa5c+dsy5cvX25IMlauXGlb1qVLF6N58+bG5cuXbcsyMzON22+/3WjYsGG2fYeFhRkZGRm25enp6UaVKlWM9u3bG1evXrUtnzt3riHJuPPOO23LVq9ebUgyfvjhB7vjatGihV05M2bHfuMrNjbWVv711183bvw4ev/99w1JxpkzZ0z3ERsba0gy5syZY7f8ypUrRvXq1Y1mzZoZly5dsi3/7rvvDEnGuHHjbMuaN29uBAYGGmlpabZlGzZsMCQZdevWtS3Lep/8/PyM06dP2+0vIyPDSE9Pt1t2/vx5IyAgwPjnP/+ZrY5q1aoZycnJtuVjxowxJBktW7a0e18iIiKMcuXK2b3fOcl6v288n4ZhGHfeeachyZg+fbrd8g8++MCQZMyfP9/unHXs2NHw9fU1UlNTC629u3fvNiQZzz33XK7H1KJFC6Ny5cq2n7N+Pw4cOGCcOXPGiI+PN2bMmGFYrVYjICDAuHjxot25yOvvDnAGl96AQjRt2jStXbvW7vXDDz9kK/foo4+qUqVKtp+zekgOHTokSTp37pzWr1+vvn37Ki0tTWfPntXZs2eVlJSke++9V3/++We2y0qDBw+26w3asWOHkpKSNHjwYLvxQP3797fbtyR17dpVt9xyixYsWGBbtm/fPv3yyy8aMGCAQ8f+8MMPZzv2tWvX6sUXX8xz26wen+XLlyszM9Oh/WXZsWOHTp8+reHDh8vLy8u2/P7771fjxo31/fffS7o+GHzv3r16/PHH5evrayt35513qnnz5jnW/cgjj6hatWp2y8qWLWsbp5SZmalz584pIyND7dq1065du7LV0adPH/n7+9t+vu222yRJAwYMsHtfbrvtNl25ciXHy4WOslqtGjRokN2y6Oho1ahRQxEREbZlnp6eioyM1IULF7Rx48ZCa29aWpokqUKFCrm2u0KFCrayN2rUqJGqVaum4OBgPf3002rQoIG+//77bJfwHP27AxzBpTegEIWEhDg0qLROnTp2P2cFl/Pnz0uSDh48KMMwNHbsWI0dOzbHOk6fPm13yejvlzYOHz4sKfudPx4eHnbjcSSpTJky6t+/vz755BPbAOAFCxbIy8tLffr0yfN4JCkwMDDHu7COHj2a57aPPvqoZs6cqaeeekqjR49Wly5d1KtXL/Xu3VtlyuT+/7us42zUqFG2dY0bN1ZMTIxduZzuhGrQoEGOIcfsctG8efP07rvv6vfff9fVq1dzLf/39zorhNSuXTvH5Vm/A/lRq1atbIPNDx8+rIYNG2Y7j02aNLGtL6z2ZgWknELQjdLS0lS9evVsy7/55hv5+fnJ09NTgYGBdpcZb+To3x3gCIIS4AZu7Pm5kWEYkmTrVXnhhRd077335lj271/43t7eLrXp8ccf13/+8x8tW7ZMERERWrhwoR544AG73oXC4u3trZ9++kk//vijvv/+e61atUqLFy9WeHi41qxZY3q+iqJdfzd//nwNHDhQPXr00Isvvqjq1avbBrxnjaG5kVnb8/odKKj2Oqsg29uwYUN5eHjol19+MS2Tnp6uAwcOKCQkJNu6O+64w3bXG1BUCEpACVCvXj1J1y+R5HeunLp160q63jt111132ZZnZGQoISFBLVq0sCvfrFkztW7dWgsWLFBgYKCOHDmiDz/8MJ9H4LwyZcqoS5cu6tKli9577z29/fbbevXVV/Xjjz+qa9eupjN5Zx3ngQMHbHdFZTlw4IBt/Y3n4+9yWmbm66+/Vr169fTtt9/aten11193uI6iVLduXf3yyy/KzMy061X6/fffbesLS/ny5dWlSxf997//1eHDh3Pc11dffaX09HSHey6BwsYYJaAEqF69ujp37qwZM2boxIkT2dafOXMmzzratWunKlWq6LPPPlNGRoZt+YIFC0wvl/zjH//QmjVr9MEHH6hKlSrq1q1b/g/CCTndZp5191h6erok2eaF+vvjOtq1a6fq1atr+vTptrKS9MMPP+i3337T/fffL0m65ZZb1KxZM33++ed2d+Jt3LhRe/fudbitWT0rN/akbNu2TVu2bHG4jqLUvXt3nTx5UosXL7Yty8jI0IcffihfX1/deeedhbr/1157TYZhaODAgbp06ZLduvj4eL300kuqXbu2/vGPfxRqOwBH0aMEFKIffvjB9j/1G91+++22XiJHTZs2TWFhYWrevLkGDx6sevXq6dSpU9qyZYuOHj2qPXv25Lp9uXLl9MYbb+jZZ59VeHi4+vbtq4SEBM2dO1f169fPsYemX79+eumll7R06VINGzbMNlVBYZswYYJ++ukn3X///apbt65Onz6tjz/+WIGBgQoLC5N0fc6dihUravr06apQoYJ8fHx02223KTg4WP/+9781aNAg3XnnnYqIiLBNDxAUFKRRo0bZ9vP222/r4YcfVmhoqAYNGqTz58/ro48+UrNmzRyaxkCSHnjgAX377bfq2bOn7r//fsXHx2v69Olq2rSpw3UUpSFDhmjGjBkaOHCgdu7cqaCgIH399dfatGmTPvjggzwHWrsqLCxM77//vkaOHKkWLVpo4MCBqlmzpn7//Xd99tlnKlOmjJYtW+bQBJpmCvLvDiAoAYVo3LhxOS6fM2eO0x/YTZs21Y4dOzR+/HjNnTtXSUlJql69ulq3bm26n78bMWKEDMPQu+++qxdeeEEtW7bUihUrFBkZaXeHWJaAgADdc889io6OLtL/4T/00ENKSEjQ7NmzdfbsWVWtWlV33nmnxo8fbxsj5enpqXnz5mnMmDEaOnSoMjIyNGfOHAUHB2vgwIEqX768Jk+erJdfflk+Pj7q2bOn/v3vf9t9AT/44INatGiR3njjDY0ePVoNGzbU3LlzNW/ePP36668OtXXgwIE6efKkZsyYodWrV6tp06aaP3++lixZog0bNhTC2XGNt7e3NmzYoNGjR2vevHlKTU1Vo0aNNGfOHNtklYUtMjJSbdq0sU3gmZSUJMMwVL16de3Zs0c1atRwqf6C/LsDLIYrIwUBlHiZmZmqVq2aevXqpc8++yzb+p49e2rv3r1Ojdsp6Vq1aqVq1app7dq1xd2Um8bEiRM1btw4vfrqq3rzzTeLuzmADWOUgJvI5cuXs92V9Pnnn+vcuXO2R5jc6MSJE/r+++9L7XiRq1ev2o3Xkq4/HmXPnj05ng8UnrFjx2ro0KF66623HHp8CVBU6FECbiIbNmzQqFGj1KdPH1WpUkW7du3SrFmz1KRJE+3cudM25058fLw2bdqkmTNnKjY2Vv/73/9cvhzijhISEtS1a1cNGDBAt9xyi37//XdNnz5d/v7+2rdvH4+8AMAYJeBmEhQUpNq1a2vq1Kk6d+6cKleurMcff1yTJ0+2m5hw48aNGjRokOrUqaN58+aVypAkXZ/Ys23btpo5c6bOnDkjHx8f3X///Zo8eTIhCYAkepQAAABMMUYJAADABEEJAADABGOUXJSZmanjx4+rQoUKpo9UAAAA7sUwDKWlpemWW27J9WHbBCUXHT9+PNtTtAEAQMmQmJiowMBA0/UEJRdlTfefmJgoPz+/Ym4NAABwRGpqqmrXrp3nY3sISi7Kutzm5+dHUAIAoITJa9gMg7kBAABMEJQAAABMEJQAAABMuH1QWrJkiTp37qxKlSrJx8dHLVu21JQpU3T16lWn67p48aImTZqkdu3ayc/PT56enqpRo4YeeOABrVixohBaDwAASjK3foTJyJEjFRUVJQ8PD4WHh8vX11fr169XcnKywsLCtGbNGnl7eztUV1JSku644w7t379fvr6+uv3221WxYkUdPHhQu3btkiRFRkYqKirKqTampqbK399fKSkpDOYGAKCEcPT72217lJYtW6aoqCj5+vpq27ZtWr16tb755hv9+eefat68uWJiYjR27FiH65swYYL279+vtm3b6vDhw1q9erUWL16snTt36vvvv5eHh4emTp2qrVu3FuJRAQCAksRtg9Lbb78tSRo9erTatGljW161alV9/PHHkqSPPvpIKSkpDtW3fv16SdLLL7+sypUr263r3r277rrrLknSli1bXG47AAAoHdwyKB07dkyxsbGSpH79+mVbHxYWptq1ays9PV3R0dEO1enl5eVQuapVqzreUAAAUKq5ZVCKi4uTJFWuXFnBwcE5lmnXrp1d2bx069ZNkvTvf/9b586ds1sXHR2tH3/8UTVq1NBDDz2U32YDAIBSxi1n5o6Pj5ck1alTx7RM1vPVssrm5eWXX9b27du1evVq1a1bV6GhobbB3Dt37lRoaKhmzZolf3//XOtJT09Xenq67efU1FSH9g+g4FzLNLQ9/pxOp11W9QpeCgmurLJleCg1gILnlkEpLS1NkuTj42NaxtfXV5LjQcXHx0crV67UK6+8onfffVerV6+2ratSpYq6du2qWrVq5VnPpEmTNH78eIf2CaDgrdp3QuNX7teJlMu2ZTX9vfT6g011X7OaxdgyAKWRW156KwwnTpxQaGioPvzwQ7355ps6dOiQLly4oO3bt6tt27YaP368wsLCbCHNzJgxY5SSkmJ7JSYmFtERAFi174SGzd9lF5Ik6WTKZQ2bv0ur9p0oppYBKK3cMihlPcn34sWLpmUuXLggSQ7PXfTEE08oNjZWEydO1CuvvKLg4GD5+Pioffv2+u6779S8eXPt2bNH77zzTq71WK1W2wNweRAuUHSuZRoav3K/cpr4LWvZ+JX7dS3TbaeGA1ACuWVQCgoKkqRce2uy1mWVzc2xY8e0du1aSVJERES29Z6enurdu7ck6b///a+TrQVQFLbHn8vWk3QjQ9KJlMvaHn/OtAwAOMstg1Lr1q0lXZ9N22yw9o4dOyTJbo4lM0eOHLH926wHKGsQ99/viAPgHk6nmYek/JQDAEe4ZVAKDAxU+/btJUkLFy7Mtj4mJkaJiYmyWq3q3r17nvXdOEh727ZtOZbJmpHbbDoCAMWregXH5kJztBwAOMItg5IkvfLKK5KkyZMn257FJl3vZRo+fLgkacSIEXa38y9dulSNGzdWly5d7OqqU6eOLXg999xzSkhIsFs/f/58LV68WFLOE1wCKH4hwZVV099LZpMAWHT97reQ4MomJQDAeW45PYAk9ejRQ5GRkZo6dao6dOigLl26yMfHR+vWrVNycrJCQ0M1ceJEu21SUlJ04MABXb6cvet99uzZuuuuu/Tbb7+pSZMm6tChg6pWrarffvtNv/76qyRpwIAB6t+/f5EcHwDnlC1j0esPNtWw+btkkewGdWeFp9cfbMp8SgAKlNv2KElSVFSUFi9erI4dO2rz5s2Kjo5WYGCgJk+erPXr18vb29vhupo1a6Z9+/bp5Zdf1q233qrY2FgtW7ZMp0+f1r333qvFixfriy++kMXChyzgru5rVlOfDGijGv72l9dq+HvpkwFtmEcJQIGzGIbBvbQuSE1Nlb+/v1JSUpgqACgizMwNwFWOfn+77aU3ADBTtoxFHetXKe5mALgJuPWlNwAAgOJEUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADBBUAIAADDh9kFpyZIl6ty5sypVqiQfHx+1bNlSU6ZM0dWrV/Nd5/Lly/XQQw+pRo0aKleunKpXr67bb79dEyZMKMCWAwCAks5iGIZR3I0wM3LkSEVFRcnDw0Ph4eHy9fXV+vXrlZycrLCwMK1Zs0be3t4O13flyhUNGDBAS5Yskbe3tzp27KiAgACdPHlSv/76q65du6azZ8861cbU1FT5+/srJSVFfn5+zh4iAAAoBo5+f3sUYZucsmzZMkVFRcnX11cbN25UmzZtJElnz55VeHi4YmJiNHbsWL3zzjsO1zl48GAtWbJEPXr00GeffaaqVava1mVmZmr79u0FfhwAAKDkctsepZCQEMXGxurNN9/Uq6++arcuJiZGnTp1ktVq1alTp+Tv759nfevWrVPXrl3VrFkz7dq1S56engXSTnqUAAAoeRz9/nbLMUrHjh1TbGysJKlfv37Z1oeFhal27dpKT09XdHS0Q3V++OGHkq5fziuokAQAAEo3t7z0FhcXJ0mqXLmygoODcyzTrl07JSYmKi4uThEREbnWd+3aNa1bt06SdMcdd+jkyZP68ssvdeDAAVmtVrVu3VqPPPKIfH19C/ZAAABAieaWQSk+Pl6SVKdOHdMytWvXtiubm0OHDunChQuSpK1bt2r48OG2n7O8+OKL+vLLLxUeHp7fZgMAgFLGLS+9paWlSZJ8fHxMy2T1/qSmpuZZX1JSku3fTz75pNq2bavY2FilpaVp9+7d6t69u86cOaOHH35Yf/75Z651paenKzU11e4FAABKJ7cMSgXtxvHqtWrV0urVq9WuXTv5+vqqZcuWWrFihZo1a6YLFy5o8uTJudY1adIk+fv7215ZPVsAAKD0ccugVKFCBUnSxYsXTctkXTpz5E6zrPokaeDAgbJarXbry5Ytq6efflqS9N///jfXusaMGaOUlBTbKzExMc/9AwCAksktxygFBQVJUq4hJGtdVtm86rNYLDIMQ/Xq1cuxTNbyEydO5FqX1WrNFrQAAEDp5JY9Sq1bt5Z0fWyR2WDtHTt2SJJtIsrc+Pr6qlGjRpJkOvN21nLufAMAAFncMigFBgaqffv2kqSFCxdmWx8TE6PExERZrVZ1797doTr79OkjyfzS2tq1ayVdn+gSAABActOgJEmvvPKKJGny5MnatWuXbXlSUpKGDx8uSRoxYoTdrNxLly5V48aN1aVLl2z1RUZGqlKlSoqOjtaMGTPs1n355ZdasGCBrRwAAIDkxkGpR48eioyM1IULF9ShQwd169ZNvXv3VoMGDbR3716FhoZq4sSJdtukpKTowIED+t///petvqpVq2rx4sXy8vLS0KFD1axZM/Xp00dt2rRRRESEDMPQ2LFjHe6hAgAApZ/bBiVJioqK0uLFi9WxY0dt3rxZ0dHRCgwM1OTJk7V+/Xp5e3s7Vd/dd9+tPXv26IknnlBycrKWL1+uI0eOqHv37lq9erUmTJhQSEcCAABKIrd9KG5JwUNxAQAoeUr0Q3EBAADcAUEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADABEEJAADAhIerFVy7dk3btm3Trl27dOrUKZ0/f16VKlVSQECA2rZtq5CQEJUtW7Yg2goAAFCk8h2UYmJiNG3aNH3//fe6ePGibblhGLJYLLaffX19df/99+uZZ55RaGioa60FAAAoQhbDMAxnNvj55581atQoxcXFyTAMlSlTRs2aNdP//d//qUqVKvLz81NKSoqSkpK0b98+7d+/X5mZmbJYLGrTpo3ee+89derUqbCOp8ilpqbK399fKSkp8vPzK+7mAAAABzj6/e1Uj9Jjjz2mJUuWyMPDQw899JAGDhyo8PBwVahQIdeGrFu3TnPnztWqVavUuXNn9e3bV4sWLXJm1wAAAEXOqR4lq9WqIUOG6LXXXlNAQIDTOzt16pQmTJigWbNm6fLly05v747oUQIAoORx9PvbqaCUkJCgoKAglxsXHx+v4OBgl+txBwQlAABKHke/v52aHqAgQpKkUhOSAABA6ebSPErnz58vqHYAAAC4HZeCUmBgoAYPHqzdu3cXUHMAAADch0tB6cqVK5o1a5batm2rTp06afHixcrIyCiotgEAABQrl4LSkSNHNHbsWAUEBGjTpk3q16+f6tSpo/Hjx+vkyZMF1UYAAIBi4fSEkznJyMjQ119/rWnTpmnTpk2yWCzy8PBQr1699MwzzygsLKwg2uqWuOsNAICSp1DuejPj4eGhxx57TD///LP27Nmjp556SuXKldPixYt15513qnXr1qVq7iQAAHBzKJCgdKPmzZtrxowZOnr0qF544QUZhqFffvlFQ4YMUa1atTR27FilpqYW9G4BAAAKXIEHJen68+CGDBmiqKgoSVK5cuUUEhKi5ORkvf3222ratKl+/fXXwtg1AABAgSmwoHTp0iV9+umnatmypTp37qwlS5aoatWqmjBhgo4cOaItW7bo999/V58+fXT8+HH961//KqhdAwAAFAqnHoqbk4MHD2ratGmaN2+eUlJSZBiGQkJCFBkZqb59+8rD4//tomHDhvryyy91+PBhbd261dVdAwAAFCqXglK3bt20du1aZWZmytPTU4899pgiIyN122235bpd06ZNtX37dld2DQAAUOhcCkqrV69WtWrVNGTIEA0fPlw1a9Z0aLsePXqoTp06ruwaAACg0Lk0j9LcuXPVr18/lStXriDbVKIwjxIAACWPo9/fLvUoDRw40JXNAQAA3FqhTA8AAABQGrjUoxQeHu5QuXLlyqlq1apq166dIiIiFBAQ4MpuAQAAioRLY5TKlLneIWWxWCRJOVVlsVhsyy0Wi7y8vPTJJ5/o8ccfz+9u3QpjlAAAKHmKZIzSjz/+qO+++07vvvuu2rdvr379+ikoKEgWi0UJCQlauHChtm/frueff16tWrXS+vXrNW/ePD311FNq3LixQkJCXNk9AABAoXIpKJUrV05RUVF67733NHLkyGzrIyMjFRUVpRdffFEbNmzQgAED1LFjRz399NOKiorSggULXNk9AABAoXLp0lv37t119OhR/fLLL7mWa9GihQIDAxUdHS3DMFS/fn1du3ZNhw8fzu+u3QaX3gAAKHkc/f526a637du3q3nz5nmWa968ubZt2ybp+jilpk2b6vTp067sGgAAoNC5FJQuXbqkEydO5FnuxIkTunz5su1nHx8fu2fAAQAAuCOXglKTJk30888/23qLcrJt2zb9/PPPatq0qW3ZsWPHVLVqVVd2DQAAUOhcCkrDhw/XtWvXdM8992js2LH67bffdOnSJV26dEm///67xo0bp3vvvVeZmZkaNmyYJOmvv/5SXFyc2rZtWyAHAAAAUFhcGswtXQ9L06dPt82l9HeGYejpp5/WJ598Iknav3+//vOf/+ixxx7Tvffe68qu3QKDuQEAKHkc/f52OShJ0ooVKxQVFaUtW7bYxiJZrVZ17NhRkZGR6tGjh6u7cFsEJQAASp4iDUpZrl27prNnz0qSqlSpclMM2CYoAQBQ8hTJ9AD16tXTfffdZ/u5bNmyCggIUEBAwE0RkgAAQOnmUlA6deqUKleuXFBtAQAAcCsuBaW6desqNTW1oNoCAADgVlwKSr1799ZPP/2kM2fOFFR7AAAA3IZLQWnMmDFq0qSJ7rnnHm3evLmg2gQAAOAWXBpxff/996ts2bLas2ePOnXqpOrVqysoKEje3t7ZylosFq1bt87pfSxZskTTpk3Tnj17dOXKFTVo0ED9+/fXqFGj5Onp6UrzFR0drfvvv1+S1KVLF/33v/91qT4AAFC6uDQ9QJkyjndIWSwWXbt2zan6R44cqaioKHl4eCg8PFy+vr5av369kpOTFRYWpjVr1uQYyhxx/vx5NWvWTCdOnJBhGPkOSkwPAABAyePo97dLPUo//vijK5vnatmyZYqKipKvr682btyoNm3aSJLOnj2r8PBwxcTEaOzYsXrnnXfyVf+zzz6rU6dOaejQobZZwwEAAG5UoBNOFqSQkBDFxsbqzTff1Kuvvmq3LiYmRp06dZLVatWpU6fk7+/vVN1Lly5Vr1699OKLL6pp06YaNGgQPUoAANxEimTCycJy7NgxxcbGSpL69euXbX1YWJhq166t9PR0RUdHO1X32bNnNXToUDVq1EgTJkwokPYCAIDSqUCCkmEYio6O1muvvaann35as2fPtq07c+aM/vjjD6fGJ8XFxUmSKleurODg4BzLtGvXzq6so4YNG6azZ89q1qxZ8vLycmpbAABwc3H5OSN79uzRo48+qj///FOGYchisejq1av65z//KUlau3at/vGPf2jZsmV68MEHHaozPj5eklSnTh3TMrVr17Yr64gvv/xSX3/9tZ577jmFhoY6vB0AALg5udSjdPToUXXt2lV//PGHunXrpilTpujvQ5569OghT09PLV++3OF609LSJEk+Pj6mZXx9fSXJ4ZnBT548qWeeeUb169fX22+/7XBb/i49PV2pqal2LwAAUDq5FJTefvttJSUl6YMPPtB3332nF154IVuZ8uXLq2XLlrYxR8VlyJAhOn/+vGbOnKny5cvnu55JkybJ39/f9srq2QIAAKWPS0Fp1apVaty4sSIjI3MtFxQUpBMnTjhcb4UKFSRJFy9eNC1z4cIFSXLoTrN58+Zp5cqVGjp0qDp37uxwO3IyZswYpaSk2F6JiYku1QcAANyXS2OUjh8/rocffjjPchaLxalLVEFBQZKUawjJWpdVNjdLly6VJMXGxmYLSidPnpQk7dy507buyy+/VI0aNXKsy2q1ymq15rlPAABQ8rkUlHx8fBx6IG58fLwqV67scL2tW7eWJCUlJSk+Pj7HO9927NghSbaJKB2RtU1OkpOTtXHjRknS5cuXHa4TAACUXi5demvevLl27typs2fPmpY5fPiw9uzZo7Zt2zpcb2BgoNq3by9JWrhwYbb1MTExSkxMlNVqVffu3fOsb9myZTIMI8fXnDlzJF1/1lvWMkd6qQAAQOnnUlAaMGCA0tLS9NRTT+mvv/7Ktv7KlSsaPny4rl69qgEDBjhV9yuvvCJJmjx5snbt2mVbnpSUpOHDh0uSRowYYTcr99KlS9W4cWN16dIlP4cDAABgx6VLb4MGDdKCBQu0YsUKNW7cWPfdd5+k63MrRUZGasWKFTpy5Ii6du2qRx991Km6e/ToocjISE2dOlUdOnRQly5d5OPjo3Xr1ik5OVmhoaGaOHGi3TYpKSk6cOAAl84AAECBcKlHqWzZslq5cqUiIiJ07NgxzZw5U9L12bI/+ugjHTlyRI888oi+/fbbfNUfFRWlxYsXq2PHjtq8ebOio6MVGBioyZMna/369fL29nal+QAAALkqsIfi/v7774qOjtahQ4eUmZmp2rVrq1u3bmrVqlVBVO+2eCguAAAlj6Pf3y4/wiRL48aN1bhx44KqDgAAoNgVyENxAQAASqMC61E6duyYjh07lutA6jvuuKOgdgcAAFDoXA5Ky5cv1+jRo/XHH3/kWs5isSgjI8PV3QEAABQZl4LSDz/8oEceeUSZmZny9/dXvXr1GNAMAABKDZeC0ltvvaXMzEy98cYbGj16tMqVK1dQ7QIAACh2Lk0P4Ovrq1tvvdVu5uybDdMDAABQ8jj6/e3yhJNMCQAAAEorl4JSixYtdPTo0YJqCwAAgFtxKSiNHDlSmzZt0o4dOwqqPQAAAG7DpaD0yCOPaOzYsbr33nv18ccf68iRIwXVLgAAgGLn0mDusmXLOr6jUjqPEoO5AQAoeYrkWW/OZKwCevYuAABAkXEpKGVmZhZUOwAAANwOD8UFAAAwQVACAAAw4VRQ+vzzz7V58+Yc16Wmpury5cs5rlu0aJGef/5551sHAABQjJwKSgMHDtTMmTNzXFepUiU988wzOa5bs2aNoqKinG8dAABAMSqwS2+GYXBnGwAAKFUYowQAAGCCoAQAAGCCoAQAAGCCoAQAAGCCoAQAAGDC6UeYHDx4UJ9//rlT6w4ePOh8ywAAAIqZxXDinv4yZcrIYrE4vRPDMGSxWHTt2jWnt3V3jj59GAAAuA9Hv7+d6lGqU6dOvoISAABASeRUUEpISCikZgAAALgfBnMDAACYICgBAACYcCooXbp0qUB2WlD1AAAAFCanglL9+vU1ffr0fN+9lpGRoWnTpql+/fr52h4AAKAoORWUbrnlFg0fPlxBQUF67bXX9Oeffzq03YEDBzRmzBgFBQXp2WefVa1atfLVWAAAgKLk1DxKhmFo5syZeu2113TmzBlZLBYFBgaqY8eOatKkiapUqSI/Pz+lpqYqKSlJ+/fv15YtW3Ts2DEZhqFq1arprbfe0pNPPllqphlgHiUAAEoeR7+/nQpKWS5fvqwvvvhCH330kfbu3Xu9ohyCT1bVLVq00IgRI9S/f395e3s7uzu3RlACAKDkKdSgdKOEhAStX79ecXFxOnXqlFJSUlSxYkVVr15dbdq00V133aWgoCBXduHWCEoAAJQ8RRaUbnYEJQAASp5CeYRJloMHD+rbb79VQkKCrFarWrVqpb59+5a6y2oAAODm5nRQ+uCDD/TSSy9lmyJg7Nixio6OVrNmzQqscQAAAMXJqekBYmJi9K9//UsZGRkqX768Wrdurfr168tisejo0aN65JFHlJmZWVhtBQAAKFJOBaWPPvpIhmHoiSee0MmTJ7Vjxw798ccf2rVrl+rXr6+DBw9q1apVhdVWAACAIuVUUNqyZYsCAwM1Y8YM+fj42Ja3aNFCUVFRMgxDW7duLfBGAgAAFAengtKpU6fUrl07lStXLtu6sLAwSdLp06cLpmUAAADFzKmgdOXKFVWsWDHHdVm31l25csXlRgEAALgDp4ISAADAzcTp6QEOHjyozz//PF/rH3/8cWd3BwAAUGycmpm7TJky+X6YrcViUUZGRr62dWfMzA0AQMlTKDNz16lTJ99BCQAAoKRxKiglJCQUUjMAAADcD4O5AQAATBCUAAAATBCUAAAATBCUAAAATBCUAAAATBCUAAAATBCUAAAATBCUAAAATBCUAAAATBCUAAAATBCUAAAATBCUAAAATBCUAAAATBCUAAAATBCUAAAATBCUAAAATBCUAAAATBCUAAAATBCUAAAATLh9UFqyZIk6d+6sSpUqycfHRy1bttSUKVN09epVp+qJi4vTpEmT1KVLFwUEBMjT01OVKlVSp06dNG3aNKfrAwAApZ/FMAyjuBthZuTIkYqKipKHh4fCw8Pl6+ur9evXKzk5WWFhYVqzZo28vb3zrCcjI0Oenp6SJF9fX7Vv314BAQE6evSotmzZomvXrikkJESrV69WxYoVnWpjamqq/P39lZKSIj8/v/wcJgAAKGKOfn+7bY/SsmXLFBUVJV9fX23btk2rV6/WN998oz///FPNmzdXTEyMxo4d63B9bdu21VdffaWzZ89q/fr1WrRokX7++WfFxcWpZs2a2r59u55//vlCPCIAAFDSuG2PUkhIiGJjY/Xmm2/q1VdftVsXExOjTp06yWq16tSpU/L393dpX/Pnz9c//vEPeXt7KyUlxdb75Ah6lAAAKHlKdI/SsWPHFBsbK0nq169ftvVhYWGqXbu20tPTFR0d7fL+WrduLUm6dOmSzp4963J9AACgdHDLoBQXFydJqly5soKDg3Ms065dO7uyrvjzzz8lSeXKlVPlypVdrg8AAJQObhmU4uPjJUl16tQxLVO7dm27svllGIamTJkiSXrggQdktVpdqg8AAJQeHsXdgJykpaVJknx8fEzL+Pr6Srp+jdEV48eP15YtW+Tr66vJkyfnWT49PV3p6em2n13dPwAAcF9u2aNUVD7//HNNmDBBZcqU0ezZs9WwYcM8t5k0aZL8/f1tr6yeLQAAUPq4ZVCqUKGCJOnixYumZS5cuCBJ+b7TbMmSJfrnP/8pSfrss8/Up08fh7YbM2aMUlJSbK/ExMR87R8AALg/t7z0FhQUJEm5hpCsdVllnfHtt9+qX79+yszM1IwZM2yByRFWq5VxTAAA3CTcskcp63b9pKQk08HaO3bskCS1adPGqbqXLVumxx57TNeuXdMnn3yiwYMHu9ZYAABQarllUAoMDFT79u0lSQsXLsy2PiYmRomJibJarerevbvD9a5cuVJ9+/ZVRkaGPvnkEz399NMF1mYAAFD6uGVQkqRXXnlFkjR58mTt2rXLtjwpKUnDhw+XJI0YMcJuVu6lS5eqcePG6tKlS7b6oqOj1bt3b2VkZGj69OmEJAAAkCe3HKMkST169FBkZKSmTp2qDh06qEuXLvLx8dG6deuUnJys0NBQTZw40W6blJQUHThwQJcvX7Zbfvr0afXq1UtXrlxRYGCgNm/erM2bN+e433feeUdVq1YttOMCAAAlh9sGJUmKiopSaGiopk2bps2bN+vq1auqX7++Ro8erVGjRqlcuXIO1fPXX3/Z5j46evSo5s2bZ1r2jTfeICgBAABJbvxQ3JKCh+ICAFDylOiH4gIAALgDghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJj+JuAAC46lqmoe3x53Q67bKqV/BSSHBllS1jKe5mASgFCEoASrRV+05o/Mr9OpFy2baspr+XXn+wqe5rVrMYWwagNODSG4ASa9W+Exo2f5ddSJKkkymXNWz+Lq3ad6KYWgagtCAoASiRrmUaGr9yv4wc1mUtG79yv65l5lQCABxDUAJQIm2PP5etJ+lGhqQTKZe1Pf5c0TUKQKlDUAJQIp1OMw9J+SkHADkhKAEokapX8CrQcgCQE4ISgBIpJLiyavp7yWwSAIuu3/0WEly5KJsFoJQhKAEokcqWsej1B5tKUrawlPXz6w82ZT4lAC4hKAEose5rVlOfDGijGv72l9dq+HvpkwFtmEcJgMuYcBJAiXZfs5q6u2kNZuYGUCgISgBKvLJlLOpYv0pxNwNAKcSlNwAAABMEJQAAABMEJQAAABMEJQAAABMEJQAAABMEJQAAABMEJQAAABMEJQAAABMEJQAAABMEJQAAABMEJQAAABMEJQAAABMEJQAAABMEJQAAABMEJQAAABMEJQAAABMEJQAAABMEJQAAABMexd0AAMjLtUxD2+PP6XTaZVWv4KWQ4MoqW8ZS3M0CcBMgKAFwa6v2ndD4lft1IuWybVlNfy+9/mBT3desZjG2DMDNgEtvANzWqn0nNGz+LruQJEknUy5r2PxdWrXvRDG1DMDNgqAEwC1dyzQ0fuV+GTmsy1o2fuV+Xcu0L3Et09CW/yVp+e5j2vK/pGzrAcAZXHoD4Ja2x5/L1pN0I0PSiZTL2h5/Th3rV5HEZToABY8eJQBu6XSaeUjKqVxRXaajxwq4udCjBMAtVa/g5XC5vC7TWXT9Mt3dTWu4dLccPVbAzYceJQBu51qmoUzDUEVvT9MyFl0PKSHBlZ26TJdfDCwHbk5uH5SWLFmizp07q1KlSvLx8VHLli01ZcoUXb16NV/17dy5U3369FFAQIC8vLwUHBysZ599VqdPny7glgMwk9vlq1X7Tijs3+vVf+Y2JV/K+e88q0/o9QebqmwZi8OX6TYdPJOvS2X5HVgOoOSzGIbhtn/ZI0eOVFRUlDw8PBQeHi5fX1+tX79eycnJCgsL05o1a+Tt7e1wfV9//bUiIiKUkZGh9u3bKzg4WDt27NChQ4cUEBCgmJgYNWjQwKk2pqamyt/fXykpKfLz83P2EIGbTk6Xr2r4WRURUkepl65q1qaEPOv4++WuLf9LUsRnWx3af9a+gqr6ODx5paP1LxrcwTawHIB7c/T7223HKC1btkxRUVHy9fXVxo0b1aZNG0nS2bNnFR4erpiYGI0dO1bvvPOOQ/UdP35cTzzxhDIyMjRjxgwNGTJEknTt2jUNHDhQ8+fPV79+/bRt2zZZLMz4CzjLkdmzsy5f/f1/ZydT0/X+f//Mcx8Vy3tqWkQbdahfxa7ukODKqunvpZMpl3Ps9cltX38PXTkdh7MDywGUHm7boxQSEqLY2Fi9+eabevXVV+3WxcTEqFOnTrJarTp16pT8/f3zrO+ll17Sf/7zH3Xt2lVr1661W3fhwgUFBgYqJSVFq1at0r333utwOwujR4nHNRS94jznf99327qVtPPweYfa4uy2hXWcufUSBVX1UVVfqzIzDT27KM70cpqjzHptskKYpDzDUk6eDA2Sn7enFm0/opOp6bblNfysCmtQTV/vOppnHQuevE1l/v9Lgc6+l3nJ670rqZ8bJbXdhY3zUvjnwNHvb7cMSseOHVNgYKAk6dChQwoODs5Wpk6dOkpMTNTChQsVERGRZ50NGzbUwYMHNXv2bA0aNCjb+scff1xffPGFhgwZohkzZjjc1oIOStxVU/SK85zntO8yFunGoS5mbXF228I6TrNeosIS9VgrPdyqlmlb/n6MRcEiyb+8p7w8yupkqvPvZV7yeu9K6udGSW13YeO8FM05cPT72y0Hc8fFxUmSKleunGNIkqR27drZlc1NWlqaDh48aLedK/UVFu6qKXrFec7N9v338cA5tcXZbSdF7y+U48xtkHNhyW3agPua1VTMy+EacZdzYw1dYdH1Hqzkv67ahSTJsfcyL3n9jhbWe1vY+LzLGefF/c6BWwal+Ph4Sdd7jczUrl3brmxuEhISbP82q9OZ+goDd9UUveI8584EjL+3xdltDUmf/RxfKMeZ1235BenG6QByU7aMRaENqhZJmyQpwM+qiuXNpzG4kbPn25Hf0cJ6bwsTn3c547y45zlwy6CUlpYmSfLx8TEt4+vrK+l615mj9eVWp6P1paenKzU11e5VEIpiHhjYK85z7mzAuLEt+QknuX2muHKcRTV4+e/TAeQla3B3YY/oGHt/E73bt5WS/3J83JUz59uR39HCem8LE593OeO8uOc5cMug5M4mTZokf39/2yurJ8pV3FVT9IrznOe3ztNplwvtdyA/9To6e7aravh76ZMBbRwem1C2jEWvP9hUkgo1LFWtYNXZC+l5F8yBI+e7oN5rd/vc4PMuZ5wX9zwHbhmUKlSoIEm6ePGiaZkLFy5IkkMDqLPqy61OR+sbM2aMUlJSbK/ExMQ89+8IZx7XgIJRnOc8v3VWr+BVaL8D+am3sHpusur7Z2iQFg3uoJiXw50ewHlfs5r6ZEAb1fAvvL8ZV94PR7YrqPfa3T43+LzLGefFPc+BW86jFBQUJEm5hpCsdVllc1O3bl3bv48cOaLmzZvnuz6r1Sqr1ZrnPp2V1zwwFl3/X3Ve4zPguOI8587M+5NTW5zZVrp+95Vh5HzbvCvHmdVzM2z+LtugZmdULO+pJzoGaXFsot1A6BoFdHfLfc1q6u6mNWy3GCec/ev/v/3ftf+NuvJ+OHO+HfkdtVjML7+56+cGn3c547y45zlwyx6l1q1bS5KSkpJMB1fv2LFDkmwTUebGz8/PNuN21nau1FcYcrtU4Oz4DDimOM+5M5eG/t4WZ7e1SBrcKTjH8gVxnPnpuclq1+RezTXq7lu1aXS4Fg3uoKjHWuW7B8lM2TIWdaxfRQ+3qqXnuja07eufoUG2tjjD1ffjxm0daXtev6ODOwXbzqcr+ypKfN7ljPPinufALYNSYGCg2rdvL0lauHBhtvUxMTFKTEyU1WpV9+7dHaqzZ8+epvVduHBBK1eulCT16tUrv812mdkXjrPjM+C44jznZvv++99/Tm1xdtsx3ZsW6nFm3ZafFXZGdb1VNfzMg9Pf93tjmOn4t1m3C1rWvsY9+H+ansM5qenvpVFdbzU9joJ4P5w533n9jhb2e1tY+LzLGefF/c6BW044KV1/hEnPnj2zPcIkKSlJd911l/bu3at//etfdo8wWbp0qcaMGaNatWpp3bp1dvUdP35cDRs21F9//aVPP/1UgwcPlnT9ESaDBg3SF198ofbt2zv9CBNm5i4dmJm7cI+rqq9VMqSzF9Pd7ne6IGe8duW9LMp2upOS2u7CxnlhZm6HPPfcc5o6dao8PT3VpUsX+fj4aN26dUpOTlZoaKjWrl1r91DcuXPnatCgQapbt67d3ElZlixZooiICF27dk233XabgoKCFBsby0NxAQC4yZTombmzREVFafHixerYsaM2b96s6OhoBQYGavLkyVq/fr1dSHJEnz59tG3bNvXq1UuHDh3S0qVLde3aNT3zzDPas2eP0yEJAACUbm7do1QS0KMEAEDJUyp6lAAAAIoTQQkAAMAEQQkAAMAEQQkAAMAEQQkAAMAEQQkAAMAEQQkAAMCER3E3oKTLmoYqNTW1mFsCAAAclfW9ndd0kgQlF6WlpUmSateuXcwtAQAAzkpLS5O/v7/pembmdlFmZqaOHz+uChUqOPUwXbiX1NRU1a5dW4mJicywjjzx+wJn8PvingzDUFpamm655RaVKWM+EokeJReVKVNGgYGBxd0MFBA/Pz8+yOAwfl/gDH5f3E9uPUlZGMwNAABggqAEAABggqAESLJarXr99ddltVqLuykoAfh9gTP4fSnZGMwNAABggh4lAAAAEwQlAAAAEwQlAAAAEwQlIAcWiyXX12OPPVbcTUQxWLJkiTp37qxKlSrJx8dHLVu21JQpU3T16tXibhrcyMCBA/P8DLl8+XJxNxMOYsJJIBdPPPFEjstvu+22Im4JitvIkSMVFRUlDw8PhYeHy9fXV+vXr9fLL7+slStXas2aNfL29i7uZsKNhIaGqkGDBjmuK1u2bBG3BvlFUAJyMXfu3OJuAtzAsmXLFBUVJV9fX23cuFFt2rSRJJ09e1bh4eGKiYnR2LFj9c477xRzS+FOnnrqKQ0cOLC4mwEXcekNAPLw9ttvS5JGjx5tC0mSVLVqVX388ceSpI8++kgpKSnF0j4AhYegBAC5OHbsmGJjYyVJ/fr1y7Y+LCxMtWvXVnp6uqKjo4u6eQAKGZfegFy89957OnjwoCwWi+rUqaO7777brkcBpV9cXJwkqXLlygoODs6xTLt27ZSYmKi4uDhFREQUZfPgxn788Uft3btXaWlpqlKlikJCQtS9e3dm6C5hCEpALv71r3/Z/Tx69Gjdd999mjt3rgICAoqpVShK8fHxkqQ6deqYlqldu7ZdWUCSPv/882zLatasqdmzZ+u+++4rhhYhP7j0BuSgX79+WrZsmRISEnTp0iX98ccf+uijj1SlShWtWrVKd999N7f33iTS0tIkST4+PqZlfH19JUmpqalF0ia4t5YtWyoqKkr79u1TamqqTp06pTVr1uj222/XiRMn9NBDD2nDhg3F3Uw4iB4llCovvfSSVqxY4fR2M2fOVFhYmO3nBQsW2K1v2LChGjZsqO7du6t169bau3evpk+frpEjR7raZAClzKhRo+x+rlChgu6++2517dpVPXv21PLlyzVy5Ejt3r27eBoIpxCUUKocP35cBw4ccHq7CxcuOFQuODhYgwYN0gcffKCVK1cSlG4CFSpUkCRdvHjRtEzW74+fn1+RtAklk8Vi0fjx47V8+XLt2bNHiYmJtsu2cF9cekOpMn/+fBmG4fTLmfECTZo0kSQdPXq0sA4DbiQoKEiSlJiYaFoma11WWcBM1ueHxGdISUFQApyUlJQk6f/1NKB0a926taTr77vZYO0dO3ZIEndEIk9Znx8SnyElBUEJcEJmZqa++uorSVJISEgxtwZFITAwUO3bt5ckLVy4MNv6mJgYJSYmymq1qnv37kXdPJQwX375paTrl2kbNWpUzK2BIwhKwN8sWLAgx3FOp0+fVv/+/bV79255enrq2WefLYbWoTi88sorkqTJkydr165dtuVJSUkaPny4JGnEiBHy9/cvlvbBfezevVsrVqxQRkaG3fLMzEzNmjXL9rsUGRkpT0/P4mginGQxDMMo7kYA7qRHjx5avny5GjZsqKZNm8rHx0dHjhzR7t27deHCBZUvX15z585Vnz59irupKELPPfecpk6dKk9PT3Xp0kU+Pj5at26dkpOTFRoaqrVr1/JQXGjZsmXq2bOnKlWqpDZt2iggIEDJycnat2+fjhw5IkmKiIjQ559/Lg8P7qcqCQhKwN8sXbpU3377rXbv3q2TJ08qOTlZ3t7eql+/vrp06aJnnnnGdIZmlG5fffWVpk2bpt27d+vq1auqX7++BgwYoFGjRqlcuXLF3Ty4gfj4eE2dOlU7duxQfHy8kpKSZBiGAgICFBISokGDBnGJtoQhKAEAAJhgjBIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghIAAIAJghKAXAUFBcliseT5mjt3bnE3tcSJi4tT2bJlsz03cMOGDbbzmpuEhARbuYSEBJfbEx8fr3Llyqlv374u1wWUFjxoBoBDQkND1aBBA9P1ua1Dzp599ll5e3tr7Nixxd0USVJwcLCGDBmiadOmaePGjbrzzjuLu0lAsSMoAXDIU089pYEDBxZ3M0qNr7/+Wps2bdKLL76o6tWrF3dzbF577TV9+umnGjVqlHbt2lXczQGKHZfeAKAYvP/++5KkJ598sphbYq9GjRrq3r274uLi9NNPPxV3c4BiR1ACUChuHGPzzTffKCwsTH5+fvLx8VFoaKiio6NNt83IyNDMmTPVuXNnVa5cWVarVcHBwRo2bJgSExOzlc8a09O5c2f99ddfGjdunJo0aaLy5csrKCjIVs4wDM2ePVvt2rVT+fLlVaVKFXXr1k2bN2+2qyPLnDlzZLFYdO+995q29fjx4/L09JS3t7eSkpIcOjdxcXHavHmzOnTooEaNGjm0jbMOHz6sf//73woPD1edOnVktVpVsWJFhYWFacaMGcrMzDTdNqvncNq0aYXSNqAkISgBKFSvv/66+vTpI0nq3r27GjZsqM2bN+uBBx7Q0qVLs5VPS0vT3XffrcGDB2vnzp1q0aKFHnroIVmtVk2fPl2tW7dWXFxcjvu6fPmyOnfurPfee0/BwcF66KGH1LBhQ9v6Z555Rk8++aTi4uIUEhKie+65R4mJibrjjjv03XffZauvX79+qlatmtauXas//vgjx33OmDFDGRkZioiIUJUqVRw6J8uWLZMkde3a1aHy+fHFF19o9OjRSkhI0K233qpevXqpVatWio2N1dChQ9WnTx8ZhpHjtuHh4SpTpoy+//57Xb16tdDaCJQIBgDkom7duoYkY86cOU5tJ8mQZFSsWNHYunWr3brXX3/dkGTceuut2bbr16+fIcl44IEHjFOnTtmte//99w1JRsOGDY2MjAzb8h9//NG2vxYtWhgnTpzIVu/y5csNSYavr6+xadMmu3Xvvvuubfs777zTbt2rr75qSDIiIyOz1XnlyhWjRo0ahiRj586deZ6TLGFhYYYk4/vvv89x/Y3Hk5v4+Hhbufj4eLt127dvN/bu3Zttm2PHjhktW7Y0JBlfffWVad0tWrQwJBk///xz3gcElGIEJQC5ygpKeb3Onz9vt13W8qlTp2ar8/Lly4a/v78hyThy5Iht+f79+w2LxWLccsstRmpqao7t6d69uyHJWLlypW3ZjcHip59+ynG78PBwQ5IxZsyYHNe3b98+x6B07Ngxw9PT0/D39zcuXLhgt27RokWGJKNjx4451mnGx8fHkGQcOnQox/U3Ho+jr78HpdysXr3akGT06dPHtExERIQhyYiKinLq2IDShrveADgkr+kBypUrl+PyBx98MNsyq9WqevXqKS4uTseOHVPt2rUlSdHR0TIMQ926dVOFChVyrK9z586Kjo62Xb67UfXq1dWpU6ds22RkZGjz5s2SpP79++dYb79+/RQbG5tt+S233KLevXtr0aJF+uKLLzR06FDbuqwxPCNGjMixzpxcvHhRFy9elCSHLtU98cQTpusuXLigb775xnR9enq61qxZo9jYWJ0+fVrp6ekyDENpaWmSpAMHDphum9W2U6dO5dlGoDQjKAFwSH6nB6hTp06Oy/38/CRdH1eU5dChQ5KkWbNmadasWbnWe+bMmWzLbhy4faOzZ8/a9mNWxmy5JEVGRmrRokWaNm2aLSj98ssviomJUUBAgHr37p1rW2+UkpJi+7dZGLxRbhN5JiQkmAalrVu36tFHH9WRI0dMt09NTTVdl/X+nD9/Ps82AqUZQQlAoSpTxvF7RrLuxGrVqpVatmyZa9nbbrst2zJvb2/nGneD3GbB7tChg0JCQrR9+3bbRIxZvUlDhgwx7U3LScWKFW3/TktLswWSgvTXX3+pR48eOnXqlAYNGqRhw4apQYMG8vPzU9myZfXHH3+oUaNGpoO5pf8X6CpVqlTg7QNKEoISALeRdQkuNDRUH330UYHVW6VKFVmtVqWnp+vw4cNq2rRptjJ5PQIkMjJSAwYM0EcffaSWLVtqwYIF8vDwsLsU54jy5cvLx8dHFy9eVFJSUqEEpZ9++kmnTp1SmzZtNHv27Gzr//zzzzzryJrqICAgoMDbB5QkTA8AwG1069ZNkrRixQq7S3Ku8vT0VMeOHSVJCxcuzLHMokWLcq2jb9++qlmzppYtW6a33npLFy9eVM+ePXXLLbc43Z42bdpIkvbv3+/0to44d+6cJPPLnvPnz8+zjn379kmS2rZtW3ANA0ogghIAt9G6dWs98sgjSkxMVK9evXLs5bl48aIWLFjg9CDjyMhISdLUqVO1detWu3VRUVHatm1brtt7enpq2LBhysjI0DvvvCPJuUHcN7rrrrskSVu2bMnX9nlp0qSJJGndunXZwtinn36qxYsX57p9SkqK9u/fL19fX4WEhBRKG4GSgktvABwyc+ZMbdiwwXT9Pffco379+rm8nzlz5ig5OVk//PCDGjVqpJYtWyo4OFiGYSghIUF79uzRlStX9Ntvvzl1Wahnz54aMmSIPv30U4WFhalTp06qWbOm9u7dq99++02jRo3S+++/n+t4o6efflpvvfWW0tPT1aJFC91xxx35OsYePXpowoQJWrt2rd5888181ZGb1q1b6+GHH9by5cvVunVr2wznu3fv1oEDB/TKK6/orbfeMt1+/fr1yszMVPfu3eXp6Vng7QNKEoISAIds2rRJmzZtMl1fsWLFAglKFSpU0Jo1a7R48WLNnz9fO3fu1O7du+Xn56eaNWuqf//+euihh1S/fn2n654+fbrat2+vTz75RFu3bpWXl5dCQkL08ccf23qvqlatarp99erV1apVK23btk3PPPNMfg9RrVu31u23367Nmzfrt99+s/UAFaQlS5YoKipKn3/+uWJiYuTl5aV27dpp6tSpatiwYa5BKetOO1eOESgtLEZutz0AwE3in//8p+bMmaN3331Xzz//fI5l/vjjDzVu3Fj+/v46duyYypcvn+/9ff311+rTp4+ef/55vfvuu/mup6CdPHlSderUUbNmzbRr167ibg5Q7BijBOCm8euvv9ome8ySmZmpzz77THPnzpWXl5ciIiJMtx83bpwMw9CwYcNcCkmS1Lt3b4WGhmrGjBluNanjxIkTdfXqVb333nvF3RTALdCjBOCmMXDgQH311Vdq3bq1atWqpYsXL2r//v1KSEhQ2bJl9dlnn2nQoEF226xYsULLly/Xr7/+qm3btqlGjRr67bff7OZDyq+4uDi1a9dOw4YNK9DpEPLr0KFDaty4sXr06KGvvvqquJsDuAXGKAG4aTz66KNKTU21jXvKyMhQ9erV9eijj2rkyJHq0KFDtm127dql2bNnq0KFCuratavee++9AglJ0vWxSteuXSuQugpCvXr1dOXKleJuBuBW6FECAAAwwRglAAAAEwQlAAAAEwQlAAAAEwQlAAAAEwQlAAAAEwQlAAAAEwQlAAAAEwQlAAAAEwQlAAAAE/8fwHDBu4jI/h8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "phase_counts = res.parsed_counts_of_outputs(\"phase\")\n",
    "num_shots = res.num_shots\n",
    "\n",
    "energy_results = {\n",
    "    sample.state[\"phase\"] * normalization: sample.shots / num_shots\n",
    "    for sample in phase_counts\n",
    "}\n",
    "\n",
    "plt.plot(energy_results.keys(), energy_results.values(), \"o\")\n",
    "\n",
    "max_prob_energy = max(energy_results, key=energy_results.get)\n",
    "\n",
    "print(f\"\\nEnergy with maximal probability: {max_prob_energy} Ha\")\n",
    "print(f\"Precision: {(2**(-QPE_SIZE))* normalization} Ha\")\n",
    "print(f\"Classical solution:, {classical_sol} Ha\")\n",
    "plt.xlabel(\"Energy (Ha)\", fontsize=16)\n",
    "plt.ylabel(\"P(Energy)\", fontsize=16)\n",
    "plt.tick_params(axis=\"both\", labelsize=16)\n",
    "plt.title(\"Energy Histogram from QPE\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d295b803-a203-47e5-82e8-5820b0e31b26",
   "metadata": {},
   "source": [
    "Recall that we are looking for a signal from the smallest eigenvalue, under the assumption that the initial state has some overlap with the ground state. Below we estimate the energy as the first peak of the histogram, such that the corresponding probability is larger than `ASSUMED_OVERLAP`*0.4 (0.4 is the case for ASSUMED_OVERLAP=1).\n",
    "\n",
    "*Note that this is a very rough and simplistic analysis of the QPE algorithm result. One can utilize more complex tools of spectral analysis, such as Gaussian mixtures, etc. Additional assumptions, such as difference between adjacent eigenvalues or number of overlapping eigenstates, can facilitate the analysis further.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "342f2fdd-1f07-4005-aecd-10155222dab4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of maxima: 1\n",
      "maxima over the threshold at [-1.1009663875965574] Ha\n",
      "\n",
      "Lowest eigenvalue: -1.1009663875965574 Ha\n",
      "Precision: 0.13762079844956968 Ha\n",
      "Classical solution:, -1.0789776863285236 Ha\n"
     ]
    }
   ],
   "source": [
    "from scipy.signal import find_peaks\n",
    "\n",
    "ASSUMED_OVERLAP = 0.05\n",
    "\n",
    "\n",
    "def estimate_energy(data_dict, assumed_overlap):\n",
    "    max_prob = assumed_overlap * 0.4\n",
    "    data = tuple(data_dict.items())\n",
    "    data_sorted = sorted(\n",
    "        data, key=lambda x: x[0]\n",
    "    )  # sort the data according to the energy value\n",
    "    probs_sorted = [data[1] for data in data_sorted]\n",
    "    maxima = find_peaks(probs_sorted, height=max_prob)[0]\n",
    "    print(f\"Number of maxima: {maxima.size}\")\n",
    "    if maxima.size == 0 and np.all(np.array(probs_sorted) <= max_prob):\n",
    "        print(\n",
    "            \"\"\"No probabilities above threshold were found, try to increase the assumed_overlap. \n",
    "Returning energy with max probability\"\"\"\n",
    "        )\n",
    "        return max(data_dict, key=data_dict.get)\n",
    "    elif maxima.size == 0:  # strictly increasing or decreasing function\n",
    "        return max(data_dict, key=data_dict.get)\n",
    "    else:\n",
    "        print(\n",
    "            f\"maxima over the threshold at {[data_sorted[maxima[k]][0] for k in range(maxima.size)]} Ha\"\n",
    "        )\n",
    "        return data_sorted[maxima[0]][0]\n",
    "\n",
    "\n",
    "measured_energy = estimate_energy(energy_results, ASSUMED_OVERLAP)\n",
    "print(f\"\\nLowest eigenvalue: {measured_energy} Ha\")\n",
    "print(f\"Precision: {(2**(-QPE_SIZE))* normalization} Ha\")\n",
    "print(f\"Classical solution:, {classical_sol} Ha\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9797b512-d41f-47c1-8c0f-4134f0500b80",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "<a name='NC'>[1]</a>: [Michael A. Nielsen and Isaac L. Chuang. 2011. Quantum Computation and Quantum Information: 10th Anniversary Edition, Cambridge University Press, New York, NY, USA.\n",
    "](https://archive.org/details/QuantumComputationAndQuantumInformation10thAnniversaryEdition)\n",
    "\n",
    "<a name='Bounds'>[2]</a>: [M. Hagan and N. Wiebe. Composite Quantum Simulations, Quantum 7, 1881 (2023).\n",
    "](https://quantum-journal.org/papers/q-2023-11-14-1181/)\n",
    "\n",
    "<a name='Diamond'>[3]</a>: [Diamond Norm (Wikipedia)\n",
    "](https://en.wikipedia.org/wiki/Diamond_norm)\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}