{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "qpe-flow-00",
   "metadata": {},
   "source": [
    "# Phase estimation: one routine, different operators\n",
    "\n",
    "**Download Notebook** - {nb-download}`phase_estimation_demo.ipynb`\n",
    "\n",
    "Quantum phase estimation (QPE) extracts an eigenphase:\n",
    "\n",
    "$$\n",
    "U|u\\rangle=e^{i\\pi\\phi}|u\\rangle,\n",
    "\\qquad 0\\leq\\phi<2.\n",
    "$$\n",
    "\n",
    "This notebook uses the same `qpe` routine with three controlled-power oracles:\n",
    "\n",
    "- a single-qubit rotation,\n",
    "- Trotterized Hamiltonian evolution, and\n",
    "- a qubitized walk.\n",
    "\n",
    "## Generic power oracles\n",
    "\n",
    "`qpe` is generic over the type of register acted on by the oracle. Its power oracle has the form\n",
    "\n",
    "```python\n",
    "power_oracle(control: qubit, unitary_regs: UnitaryRegs, power: int) -> None\n",
    "```\n",
    "\n",
    "- `UnitaryRegs` can be a qubit array or a struct containing several registers.\n",
    "- `qpe` requests the required powers; the oracle decides how to implement them.\n",
    "- The examples use the same `qpe` routine for a powered rotation, repeated Trotter steps and repeated qubitized walks.\n",
    "\n",
    "Each example also supplies an input state and a phase register. Six phase qubits give a grid spacing of $1/32$. The phase-to-energy conversion depends on the chosen operator.\n",
    "\n",
    "The examples use [guppy](https://docs.quantinuum.com/guppy/language_guide/language_guide_index.html). The repository default is little endian.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "qpe-flow-01",
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import zixy.qubit.pauli as zqp\n",
    "from matplotlib import pyplot as plt\n",
    "from guppylang import guppy\n",
    "from guppylang.std.builtins import array, output, dagger\n",
    "from guppylang.std.angles import angle, pi\n",
    "from guppylang.std.quantum import (\n",
    "    qubit, h, ry, crz, measure_array, collect_measurements, discard_array,\n",
    ")\n",
    "\n",
    "from guppyalgos.algorithms.phase_estimation import qpe, QubitizationRegs, qubitized_power_oracle\n",
    "from guppyalgos.utils import (\n",
    "    qarray, transversal, binary_fraction,\n",
    "    phase_to_energy_qpe, phase_to_energy_qubitized_qpe,\n",
    ")\n",
    "\n",
    "n_phase = 6\n",
    "n_shots = 128\n",
    "\n",
    "\n",
    "from guppyalgos.utils.python.analysis import measurement_dataframe\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "qpe-flow-02",
   "metadata": {},
   "source": [
    "## 1. Read the phase of one qubit\n",
    "\n",
    "- Choose a rotation whose phase on $|0\\rangle$ is known:\n",
    "\n",
    "$$\n",
    "U=R_z(-2\\pi\\phi_0),\\qquad U|0\\rangle=e^{i\\pi\\phi_0}|0\\rangle,\n",
    "\\qquad \\phi_0=0.75.\n",
    "$$\n",
    "\n",
    "- The **power oracle** must apply controlled $U^k$. For a rotation, multiply its angle by $k$.\n",
    "\n",
    "- Prepare the phase register with Hadamards before calling `qpe`; it then applies controlled powers and the inverse QFT."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "qpe-flow-03",
   "metadata": {},
   "outputs": [],
   "source": [
    "target_phase = 0.75\n",
    "\n",
    "\n",
    "@guppy\n",
    "def rotation_power(control: qubit, qreg: array[qubit, 1], power: int) -> None:\n",
    "    crz(control, qreg[0], -2 * pi * target_phase * power)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "qpe-flow-04",
   "metadata": {},
   "source": [
    "### Build the main circuit\n",
    "\n",
    "- Allocate the registers, prepare the phase register, and call `qpe` with the rotation oracle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "qpe-flow-05",
   "metadata": {},
   "outputs": [],
   "source": [
    "@guppy\n",
    "def rotation_qpe() -> None:\n",
    "    qreg = qarray(1)\n",
    "    phase_qreg = qarray(n_phase)\n",
    "    transversal(h, phase_qreg)\n",
    "    qpe(phase_qreg, qreg, rotation_power)\n",
    "    output(\"phase\", collect_measurements(measure_array(phase_qreg)))\n",
    "    discard_array(qreg)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "qpe-flow-06",
   "metadata": {},
   "source": [
    "### Run and read the phase\n",
    "\n",
    "- Sample the circuit and decode the measured phase with `measurement_dataframe`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "qpe-flow-07",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  bitstring  counts  phase  empirical_probability\n",
      "0    000110     128   0.75                    1.0\n"
     ]
    }
   ],
   "source": [
    "rotation_counts = (\n",
    "    rotation_qpe.emulator(n_phase + 1).with_seed(5).with_shots(n_shots)\n",
    "    .run().register_counts()[\"phase\"]\n",
    ")\n",
    "rotation_results = measurement_dataframe(rotation_counts)\n",
    "rotation_phases = rotation_results.set_index(\"phase\")[\"empirical_probability\"].to_dict()\n",
    "assert rotation_phases == {target_phase: 1.0} \n",
    "print(rotation_results)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "qpe-flow-08",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAEjCAYAAAB+eAaMAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAASOZJREFUeJzt3Qd4FGX3NvCThCT0XkLvRZoIUqSrNAFFOooFlKIiFgQEVJCigKI0RToIFhAEAamCqIB0kCpI772FloQk8133ef+z3+5ms9kkm8wmuX/XtZnd2dnZZ2bbyVPO42cYhiFEREREZAl/a56WiIiIiIDBGBEREZGFGIwRERERWYjBGBEREZGFGIwRERERWYjBGBEREZGFGIwRERERWYjBGBEREZGFGIwRpWJbtmyRdu3aycmTJ1PVcyXWunXrtKwXL15MVceVGJ999pn06dPH6mIQpUkMxohSoOjoaPn555+lX79+8uKLL8r7778vCxculIiICIftzp49q9vdvHkzycuUnM+VWCdOnNCy3rlzJ0HHtWHDBg3Qzpw5I8nlk08+0dc7qfbx999/y5o1axK1fyJKGAZjRCnMpUuX5JFHHpGePXtKUFCQNGnSRLJmzSrvvvuuFC5cWIMH02OPPSYLFiyQ4sWLW1rmlMzVOTx16pQGaLdu3Uq2ciAA/O2335JsHwjox44dm6j9E1HCpEvg44jIIkOGDJG9e/fKP//8Iw8//LBtfe/evbWWzL62p2DBglqDQwmXVs4hgk4isgaDMaIUBkFY7ty5HQIxQO0YamvsmyrR32nMmDF6KVasmK2/1DfffCNfffWVXLhwQWbOnKlNcPXq1ZNXXnlF0qVz/Fq4f/++TJs2TbZv3y45c+bUbTJmzCgDBw7U2pTq1avHWebff/9dli1bprV6+fLlkw4dOnj84x8WFibTp0+XrVu3un3+n376SS8//vijBAYG2h5/6NAh+fDDD3X7atWqxdg/ztGcOXPk7t27Ur9+fXn55ZcdzoHzOcQ5njBhgt6HPlY479C3b1+pVauWy2OwP+fnzp3T58O5xzJ9+vR6XmbNmiUHDx7U565Zs6YG1jhOePvtt/V1x7kwA0PUiv7www96ferUqbYmRn9/f8mRI4ee3+eee06Cg4M92gf6jKEP3ZdffulQ9rjKZv6D8ODBA20GxTGtX79e7+/cubPUqVPHo9eZKC1jMyVRClOiRAm5fv26/jg6w4+l/Y+kq/5OZn8pXAYPHiyVK1fWJjjUrL3++usO+0MtG35Mhw8fLuXLl9cL+hyZj0dgEVffNvxwN2/eXPz8/LRJFRD0fPHFF3EeKwIkBIkff/yxPPTQQ/r8CIAWLVoU4/lxPrAuKirKYR9Xr17V9Qh+nP3yyy8ybNgwbfbFOXjnnXfkqaeeksjIyFjPYYUKFWyBJLbt1KmTXtBEHBvznCPw+eijj/Q4DMPQ50GQiWNDnz8Ec+XKlZMRI0ZI1apVbWVu0aKF5M+fX4Ms8/nat29v2z8CUnN969attSwIQJ988knb+YhrH676jHlSNvjzzz814H7jjTdk37590qBBA31t8DpjPRHFwSCiFOXgwYNGjhw5jEyZMhmvvvqq8e233xp79+51ue2CBQsMfMx3795tWzdt2jRd16ZNGyM6Otq2fujQoUZAQIBx+vRp27pBgwYZ/v7+Dvt/8OCBUa9ePd3H4sWL3T7X+PHjdd2KFSscyjVx4kTd77///uv2WD/88EPDz8/PYZ8RERFGrVq1Yjz/kCFDdN39+/cd9rFhwwZdv2zZshjnoGXLlkZUVJRt/Zo1a3T9hAkT3B7X3Llzdd2+ffvclt/5+Vq0aGF7PpxHHEupUqWMsmXLOpQbr0GGDBmMdu3a2dY1bdrUePjhhw1PHTp0SJ/zxx9/9GgfrVq1MipUqGC7jXJ6WrYGDRroOpwXE95beOzjjz/ucZmJ0irWjBGlMKipOHz4sNawnD59WpvqULtVoEABbQZDjYsnunbtqrVVpkaNGmktCmo2TOi4jpqpSpUqOdS+obbLEzNmzNDHogbJ3quvvqrlXLx4sdbyodnM/oKO5oBmRzx/lSpVbI9FE2SXLl3EG9DkiWY9U+PGjbX2Z/78+ZIUUG7z+XAe0Wx49OhRee2117S50oSaLdRaLVmyRJuJPYGmUNRaPv/883oOP/jgAwkICND+hQmxc+fOeJUN2+C5TXhvPfHEE3qMROQe+4wRpUB58uTRIAwXs4kO19HUiKYvNLd50tzpvE+wz72F3Fp169aN8diSJUt6VE4EjdgvmsQAAZh5QTCCYPLevXvahGfv2Wef1SAMzXuu+hyVKlVKvMHVcWDfu3bt8sr+43o+M3dZ6dKlY2xbpkwZ7YeFZlJX99tDwDR79mzp1q2bBtVZsmTRYGjp0qUSGhqaoLLGt2zoT2cf2AJe+xs3bmg/RvRPIyLXGIwRpQLog4RaLPz4oQO1J8GY84+jWUuGfl7227iqmUEA5QnUYqEmxdVoxI4dO2o/rVy5cmnZ7dWoUcP2/K6ey9U6s6M6fvjta3Lc5T2Lbd/mvrwNQZKr18DVOUZ/OYirLKi9mjJliowcOVIGDBjgcNwImBIqvmVzFWy5ek8RUUwMxohSGHSIRvOPMwQguHiSyNRTaB5ELRFqsuybND2tOUKn72PHjmmncjSZxSa21BHoWO/q+TGy05nZgR6JWNHJ3rRjx45Ynxf32Y+ADA8P12Za1Mq5Y462TGyQYTa/oqO88znAOoyaNY8Lz+nq+cxBDPZNyeCq43xs+0hs2YgocdhnjCiFwQhI9Nmyb05EsDJ+/HgdOYgmPm9Bs+d///0nX3/9tUPT46pVqzwuK4Kj9957zyHlBgIC9Bf7999/43z+I0eO6LGZ0CS7evXqGNsiQEVNDVJImPbv3+82qzzKgKZS09ChQ/Uc4nndMYMQJH9NjKJFi2qgitQU9gHuvHnzNJhCDacZhOI5EXjZj/Q0a0VRK4U+XCacc5wH51q12PaR2LIRUeKwZowohUH6iXHjxkmRIkU05xOavpBLCz+ySC2ANBTegqZEBEz44UWuL+T5QtNX//7946ztAtQwIUh488039UcctTcIyhBg1a5d25avKzbIR4bnR8d05LpCkybyZA0aNChG0IkBDGiqw7ZI04DzgmAEHdlbtmzpcv/IPYb94Lhw/tBHDWVq2LCh23IhtQXKj5xk6NOG53GXZ8wd5HnDfrA/NM+imRSd7jGjgn2zI/qFzZ07V/PLYZABnhOpMtA0PXHiRHnrrbc0JxpuIxj77rvvYtSgxraPxJaNiBLHD0MqE7kPIrIA8jwhqEEHaTQZYUSlc58kBBibN2/WUYLZsmWzdcxG8xxGOGbKlMmhH9DKlSvl0UcftSWINZ0/f15rRxC0IOBAfq62bdvqqEezg7+r57KvCUPzH2qhkOcKI0IRWHnK+fm3bdumARFqtpyDMgQie/bskbx582oAce3aNU1CinKGhITEOAeoVULwhuPH9jiXcZ1D85iwD9yPUagoD7L1uxLbObeHQNBMrIo8XuaACnvoB4b9YBomdJZHQGyfnBXHjQAL5whLvE54Le1Ho8a2DxwjmrhxnPEtG/KM4Xw8/vjjDuvxGFzatGkTo3M/Ef1/DMaIKN5QK4KkrVeuXJHs2bMn+xlEDVBswRgRUUrDf1WIyC2MdLTPao9aKfRHQhOmFYEYEVFqw2CMiNxCkxY6c2MqI0y7gz5SaPJDglkiIko8NlMSUZzQ7wojE9GnCP29nBPGJjdk7ceIPnQsR8d9IqKUjMEYERERkYXYTElERERkIQZjRERERBZKk0lfkQ8HeYvMyXSJiIiInCEV6+3bt7VvalLmykuTwRgCMc6pRkRERJ5AMulChQpJUkmTwZiZpRwnN2vWrFYXh4iIiHxQaGioVt44z27ibWkyGDObJhGIMRgjcg8TZ2NaHWS6d54qiIgoLfBL4i5N7MBPRG5hPsnu3bvrkoiIvI/BGBEREZGFGIwRERERWShN9hkjIqL/P3Q/MjLSYTJ4orQiICBA0qVLZ3maKwZjRORW5syZpUGDBrqk1CUiIkIuXLgg9+7ds7ooRJbJmDGj5M+fX4KCgiwrQ5qcmxJDVbNlyya3bt3iaEoikrSa/PrIkSNaM5AnTx79IbK6doAoOSH8wT8kV65c0Zrh0qVLx0jsmlzxAmvGiCjOH+0HDx5IYGBgkmagpuSFHyG8tsihhJoBorQoQ4YM+t126tQp/UykT5/eknLwm5WI3Prnn3/0CwpLSn0YYFNa5+8D/2RaXwIiIiKiNIzBGBERkY/CJNUVK1aUAwcOWF0USu3B2NmzZ+WPP/6QmzdvevyYEydOyJYtW+TGjRtJWjYiIvItNWrUkB9//NFh3bfffiuVK1eW1atXS2qCjuUIxO7fvy+psQP95MmTpWHDhlKzZk0ZMGCA3LlzJ9btjx49qoGpq8uSJUts27Vr1y7G/UOGDBFfZmkH/m3btsmnn36qQdWlS5dk/fr1+qK4gzdkp06dZN26dVK8eHF9cbCPd999N9nKTURE1jl48KBcu3bNdvvLL7+UgQMHysyZM6Vp06Z8aVKI4cOH62v39ddf67y3/fv3lx07dsjatWtdbo/BJvPmzXNYN3fuXPniiy/k0Ucfta1DXNCsWTPp0qWLbV3OnDnFl1laM3b48GF5+eWXZevWrR4/5uOPP5bdu3fLsWPHZN++ffLTTz9Jnz59NKAjIu/Df5VnzpzRJZGvGTRokHz00UeyePFi6dy5s0PT3po1a6Rnz57y2GOPSevWrWXXrl0Oj0W6AtTGoFamTp06MmzYMAkLC4uxb/vaN+x3586dtnVt2rSRH374wePntGc+ZuXKlfLqq69KrVq1dH/4bXN27tw5t/utV6+e7gu1gy1atJDZs2fH2MeCBQukefPm+jw9evTQEYQmjJgeM2aMVojgXLz11lua8iGp3L17V0aPHi2ffPKJvm4IohFYoaIFLWWuBAcHx6jxwvY4poIFCzpsGxIS4rBdgQIFxKcZPuDMmTPIdWasX78+zm3z5s1rDB061GFd5cqVjR49enj8fLdu3dLnw5KIKC26f/++cfDgQV2mNJkyZTLGjx+v3/s5cuQwNm3a5HD/jRs39Ds+X758xowZM4wtW7YY3bp1M3Lnzu3wvV+/fn2jatWqxurVq42lS5capUqVMjp06GC7/5tvvjHKlClju926dWsjQ4YMxqeffmp7Hn9/f2Pnzp0eP6ercmbJksWYMGGCsXHjRqNLly56TJcvX47Xsfz777/Gvn37jH/++cf47rvvjPz58xsTJ0603b927Vo9b7Nnzza2bdtmzJw502jSpInt/latWhn16tUzVq5cqecT5ShRooRx9+7dWF+Hp59+2qhQoUKsl7p168b62HXr1ulxHT9+3GF9oUKFjMGDBxue2LNnj+4Dr529hx9+2ChXrpy+tijj9OnTjejo6AR9FpIrXkhRwdjZs2d1u19//dVhfdeuXY3q1avH+riwsDA9kebFfD4GY0RxO3bsmNGuXTtdUurh7gfo/PnzGmDYX8wfTWzvfB8upkOHDsW479q1a3ofAgzn+/777794lx1BRZ48eTRo2bt3b4z7zQBm0qRJDr8D6dKlM9asWaO3EYAFBAQYJ0+etG2zdetWfRx+5M0AB7fPnTunP+YIgN59912jcePGev+SJUuM7NmzG1FRUR49Z2zlHDBggG0dnqds2bLGhx9+6PGxuIIApGLFirbbI0eO1GDLXkREhC0wwjkNDQ213YdjQjCG4C02eO0QAMZ2wfmLDfaL43J+/1WvXl1/0z3Ru3dvo0CBAkZkZKTDegTNU6dONTZv3qznLVeuXEb37t19OhhLUUlfzc76zm2/aGt215F/5MiRMnTo0CQvH1FqUmzAcl2GXzwqFxculI2Z6kpwSCldd3JUC4tLR0lpypQpMb4z0ZT03Xff6YCratWqxXiMOZkL+uk4dxtB89MLL7yg3UrefPNNh/uaNGmSoE73RYsW1aa6v//+WypVquRym0ceecShiStXrlxy+fJlvY3HIuM69mM/MABZ1tEVBs195cqV0+YuNJuhqQvH2Lt3b72OZj2sR/OgfZ4qd88Zm8cff9x2HbMgoKkQZfD0WGD79u0yceJEnVUBzZ/oCH/16lXb/U8++aR28+nWrZs8++yzOsVZlixZ9L4NGzboQAEci/k6Yom+3P/991+s5cb5SyjMhwpIuGoPx4ZzG5fw8HD5/vvv5Y033tBZJOzhfYb5JgFNsnnz5tVO/eiTVqrU/77DfE2KCsbMeaOcR5VgXjV3c0qhYyf6ldlPb4COgEREFBP6Jj3zzDMO63LkyKHLQoUKOfSZcoa+SugPZK9YsWK67NChg/Z5smcGBPGF/sYI/PBjjEACS2fmD7I9M9hAGZF93RlmI7AvP4IWDC5DYIMgCQPHMH0UBqAhGEOQ6elzxsa5HM5liGu/hw4d0nK+8847Gmxlz55d/vrrL71tql69uiZuRkCNCor27dtLr169tJ8YfkPxm4j7nKGyIzZ4jxw/fjzW+/GeQaDnCoJJuH79up5PEwZm2HfGjw36CKISBn3tnDmfK5wb2L9/P4Mxb8CXAP4DQUdGe7ht/9+NM0TauBARUdwwaTIurmA2hqpVq8b62LJly8Z6H3507X94EwvBBH54EYihpgWdzj2FGhIMBEOHfXMKHHRYR22Qfe0JAjAELAgSnnjiCds6BAN79uyR6dOne2V0KGqlTEhlUbJkSY8fj9GHKDMyC5gwKMAZavpGjBih11F7icD4+eeflzJlysjp06f1NTeDJE9gFCNqqGLjKoA0mbWrGMDXsmVLvY5zjJo9DJyIy4wZM6Rx48a2QN8d1OZCUs4tmSryjMU14hLVr+Z/CxjlsXTpUtv9qIrFGxEvChERpb1avGnTpml6o7Fjx3r8OIxaRBBm5p9C7Vrfvn216c0MuszAC0EbRkmazYlYh/xY+HGvUqVKoo8Bwd758+dtQRR+01DD5SkEuAimMOrZDO6QMsLerFmztCbPrE0zmwIzZ86sNZaoAcNzouUIEGR98803bqdBw7mKLe8XLgj+YoOauKeeekpHU969e1fLhdcC3ZDQjGrCyFHnHGEnT57UUZTdu3ePsV+UFylOzGZQBNioIUSNZt26dcVXWRqMXbx4Ud8cmzdvtp1E3MaJNn3++efStWtX2228cEju1q9fP1m0aJFWk6JN39WLQkSJly5zLsle/yVdEvmiV155RYMN/C7gN8MTaB5F3yJcEIigRgj/+M+fP9+hRsfsN4bKgPLly+s6BGUIIJz7iyUU0jogeEGAgkAENVzo6+Qp9IdCnzDUcCHoQLmQZ8u5Jgp5vdB0iFo3/HYiYMNjcC4QACJwQWBnNsWihg7XkwpeM3QxyvN/NabLly+XX375RQNEE2rKnFvD8Di8Zq1atYqxT5QXzeh4PUuUKKEtanjtEEy7685kNT/04rfqyVetWiWjRo2KsR79AMxkbfhgIcrHyTehenXSpElanYyOm+iUhw56nkLkny1bNs0x48vVlkS+0IHfFXbgT/nQPIeZTPDjZTbTpRT//vuv5MuXL8ZgLhwP+j+hqRT9qNB0Z39saGlBYIXvfxN+AlEBgE7gRYoUcfl82C/YByb4XcLzY38QHR2t6zx5ThNmnUFwhCAQTb/I+4Ugwv53KT77RTCF3zU03eH1xf4qVKgQ4/cP/bSQl8u587zZVIhyYR/OHeOTyoULF7QvOM4vBjDYQwJXBFP2ecLwemE7d92TUPOHmkI8Lq73t7vPQnLFC5YGY1ZhMEbkeTAWHXZHws4ckPSFK4h/+v/9x8pgLOVLycFYamEfjHnSaZ2Shi8EYz7fZ4yIrPXg5kW5smi4LomIKI2ntiAiIkotUNOCqY98NfcVJR8GY0RERBZA53/O+Ur6XuBpICIiIrIOgzEicssvXZAE5iqiS0p90uAYLiKf+wywmZKI3ArKXUQKdJvEs5TKmGkNkArC1bRARGnFvXv3dOkq1UdyYTBGRJQGIYcU5jA0J5tGLifnHE9Eqb1G7N69e/oZwGchufKqucJgjIjcirh0XC7+8L6EPD9agvKV4NlKRcyEpWZARpQWZc+e3fZZsAqDMSJyyzCixYi4r0tKXVAThsmhMYOJOVchUVoSGBhoaY2YicEYEVEahx8jX/hBIkqrOJqSiIiIKKUFY5jpHROSEhEREZEFwdi4ceN0xvc2bdrIr7/+KlFRUYksBhH5qsBchSTk5XG6JCIiHwnGTp48KcuWLdOOb+3atZMiRYrIwIED5ciRI94vIRFZyj8wvQSHlNIlERH5SDCG+bSaNm0q8+fPl/Pnz8uAAQNk1apVUqZMGalXr57Mnj1bwsPDvV9aIkp2kaGX5dqab3RJREQ+2IE/Z86c8sgjj+glKChIzp07J++9954UK1ZMfvvtN++UkogsE3UvVO7sXq5LIiLyoWAMNWKjRo2SsmXLSqNGjTSL7fLly+XYsWMakPXu3VteffVV75aWiIiIKJVJUDDWsmVL7Sf23XffyWuvvabB17x58zQoQxLB9OnTa9PlmTNnvF9iIiIiolQkQUlfka35r7/+ktq1a7vtV3bixInElI2IiIgo1UtQzdjSpUtjDcRy585tu45+Y0SUsgVkzC5ZHm2lSyIi8pGasWvXrrlcj7nNQkPZyZcoNUmXNbfkfLK71cUgIkq14hWMLVy40OV1iI6Oli1btkipUqW8Vzoislx0xH15cOWkBOYpJv5BGawuDhFR2g7GunXr5vI6IAEsmiUnTpzovdIRkeUeXD8nF7/rp1n4kfyViIgsDMZu3rypS9R+HT161MtFISIiIkp7EtSBn4EYERERUTLXjCGPGHTq1Ml2PTbYhoiIiIi8GIy9+eabtkDLvO6NYAyTi0+fPl0uXboklSpV0iSymTJlcvuYtWvXyooVK+TGjRuafPbll1+WEiVKePycROQ5P/8A8c+QVZdERGRhM+XVq1f1Yn89toun/vnnH53TElMrVa9eXb7//nupX7++RERExPqYkSNHSqtWrSRLliy67eHDh6VixYqya9cuj5+XiDwXlLe4FH7rB10SEZH3+RmGYYhFmjVrptMnrVy5Um8jkENN19ixY6Vnz54uH1O6dGlp06aNjB492rYO82M+88wz8vnnn3v0vMiFli1bNrl165ZkzZrVS0dDlLoUG7A81vtOjmqRrGUhIrJCcsUL8e4z5glPminDw8Nl3bp1MmXKFIfs/U888YROOB5bMIbmyFOnTtlu4wRdv36d+c2IkkjElVNyZdEIydPmQwnKU5TnmYjI6j5j3grGTp8+LZGRkVoTZq9o0aLy559/xvq4b7/9VvuVoXkT26Kp891335Xu3bu7DfxwMXGWACLPGVEPJPLmBV0SEZGFwVh8+oJ5IiwsTJeZM2d2WI/b5n2xdd7ftGmT9OjRQ0qWLCnp06eXmTNnSocOHWKtHUM/s6FDh3q1/ERERESW5RnzBrTBApoYnee9zJ7d9YTE9+/f11qxDz74QD755BN55ZVXtPm0UKFCMnDgwFifC/ehOdO8nDlzxstHQ0RERJTC8owVLlxYg679+/dL8+bNbev37dunKS5cQeB29+5deeihhxzWlytXTnbu3BnrcwUHB+uFiIiIKMWOpkTnerO50rye2CbNXr16abPj9u3bdZTCxo0bpV69erruySef1G2mTp2qucgwUhJFLVCggLRo0UJzkwFqutB/rEmTJjJ58mSPnpejKYk8H00ZHX5Pws/9K8EFHxL/4Iy6jqMpiSgtCPW10ZT2AZa3+o+hqRE1WuXLl9fL33//Lf3797cFYrBt2zbZsmWLBmNIgzFnzhzp3Lmz1p4VL15ctm7dqh35R4wY4ZUyEZEjBGAZSlTjaSEiSo15xiA6OlqDLTMDv3MnfNSaoR8ZcpKZ0FS5e/duXY9ArEqVKvF6TtaMEXleMxZ557rc+WelZK7ylKTLnFPXsWaMiNKCUF+rGXO2efNmGTdunPz77796GzVbSDFRs2bNeO3H399fateuHev9yMzvDNMl1a1bNwGlJqL4irpzXW5t+lEylKppC8aIiMji0ZSzZ8/Wvl3I3dWxY0e9IB1FnTp1tBmRiIiIiJKwZmzw4MHasR6pJezNmDFD73vppZcSslsiIiKiNCdBNWM3b96Udu3axVjfvn37GHnDiIiIiMjLwRg6zLuasuiPP/7QNBNElHr4p88smco31CUREVnYTPnrr7/ariOn13PPPadTEqGDPQZk7tixQ5su3WXCJ6KUJzB7iOR+uq/VxSAiSrU8Tm2BOSA95W5uSV/A1BZEnqe2MCIjJPL2VUmXJbf4pQvSdUxtQURpQWgypbbwj0+A5emFiFKPiKun5fzUHrokIqJUNFE4ERERESUi6av9yMrIyEiHdXHNXUlEREREiQjG0Hbap08fWbBggdy+fTvG/RbPsERERESUupsp+/XrJ0ePHpWlS5fapkYaP3685MqVS0aNGuXtMhIRERGlWgmaKLxgwYKaU6x06dLi5+cnUVFROsfkqlWr5IMPPpCdO3eKL+NoSiLPR1O6wtGURJQWhPraaEp758+fl1KlSul1FM7Muo/5Kg8cOODdEhIRERGlYgkeTYkaMShXrpwsWrRIr69evZqd94lSmQfXzsqFue/pkoiIfCQYq1Chgu06miV79+4tISEhOjfl+++/783yEZHFoh+EScT5w7okIiIfGU25f/9+2/VnnnlGmyYxHRJqyTBvJRERERElU54xQP8xsw8ZERERESVDnzGks+jYsaNUrlxZL506dZKtW7cmdHdEREREaVKCgrHZs2fryMnw8HANyHDBnJR16tSROXPmeL+URGSZdNnySa6W7+mSiIh8pJly8ODBMnXqVHnllVcc1s+YMUPve+mll7xVPiKyWECGLJK5wuNWF4OIKNXyT+h8lO3atYuxHqMpzZxjRJQ6RN27Jbd3/apLIiLykWAMIyb//PPPGOuRlf+RRx7xRrmIyEdEhl6R679N1iUREVnYTPnrr7/arjdp0kSee+456dGjh1SvXl0nBkdqCzRdDhw4MAmKSURERJTGgzFXzZKTJk2KsW748OGaCJaIiIiIvBiMYbQkEREREflInjFve/DgQbwfg+bR6OjoJCkPEf2Pf1AGSV/sEV0SEZEPBWORkZGycOFCbZYcNmyYXse6+MLjc+XKJenTp5eHHnpI1q5dG+dj9u7dq/3W8Ji8efNK3759WXNHlEQCcxaUfB2H65KIiHwkGDt27JhOFv7iiy/K/PnzZcGCBXod63Cfp77++mv5/PPP5eeff5Y7d+7I888/L08//bTbffz333+acBaB25UrV+TcuXNSuHBh2bNnT0IOhYjiYERHSXT4PV0SEZH3+Rlo64unli1bSrp06WT69OmSO3duXXf16lXp1q2b1o7Zj7x0p3Tp0rqvsWPH2tYVLVpUOnTooEGaK8hldvz4cR296efnJwkRGhoq2bJlk1u3bknWrFkTtA+i1K7YgOW6DL94VC5++46EvDxOgkP+NwftyVEtLC4dEVHSS654IUEZ+NevXy9HjhyxBWKA6xhdiQDLE9euXZOjR49K/fr1HdY3aNBAtmzZ4vIxUVFRsmLFCs3yj0AMgR+CQiIiIqI01UwZEBDgssN9RESEx8HRpUuXdJknTx6H9bh9+fJll49Bs+S9e/c0CEPi2QwZMki+fPmkf//++tyxwRyaiG7tL0REREQpNhhr1qyZzkt54sQJ2zo0HXbt2lWaNm0ar305j4bE7diaH80W1TFjxmgtHAIwNInOnDlThgwZEutzjBw5UqsZzQv6mBERERGl2GBswoQJOnqxRIkSOpoRl5IlS2pwhPs8UaBAAV0614Lhdv78+V0+Bk2hgYGB8sILL0jt2rU1aMMMAAgClyxZEutzYVYAtPealzNnzsTreImIiIiSSoI6XIWEhMimTZtk48aNcuDAAQ2KypcvL3Xr1vV4H9mzZ9fRl+vWrbNl90et2O+//64DAeybGLEeTZIIxB577LEYTZLYJigoKNbnCg4O1gsRxV9QnmJSqPf34h+ciaePiMhXgrFKlSrJvn37NPiKTwDmbMCAARp4NWzYUIOszz77TO7fvy+vv/66bZtevXpph/79+/frbUy11LZtW20qrVOnjmzdulVmzZqluc6IyPv8AtJJQMZsPLVERL4UjKGZD8196H+VGGhuRIf8oUOHaod+BHlI+mo2YQISu2bMmNF2G8lev/32WxkxYoScOnVKihQpIqNHj3YI4IjIex7cuCA3fp8mOZ7oLoE5XHchICKiZM4zhs77SLrar18/SYmYZ4wobswzRkRpXagv5xlDtnykk0D2ffQVc+6vhWSwRERERJREwZi/v7907NhRr6MzvbscX0RERETk5WBs3rx5CXkYERERETlJ1FxCyIR/+vRpvY6O9JyaiCj1SZcll+R4/FVdEhGRjyR9RV4vpKVArjAke8UF1wcNGsQmS6JUJiBTDslao7UuiYjIR2rG3n77bZ2w+6uvvtIM+LB9+3adwPvmzZs6VRERpQ5RYXck7OQ/kr5YFQlIn9nq4hARpToJCsZ++OEHzZxvBmKAbPoYWdm4cWMGY0SpSOTNi3J1ySgJeXmcBISUsro4RESpToKaKZGEtVSpmF/KpUuXdkjQSkRERERJEIxhKiJMXWSfLxbXsQ73EREREVESNlNiCqNRo0bJggULpFq1ahqI7dq1S44dOybt27d3mOibCWCJiIiIkjDpK/j5+cmjjz6qFzNDPxGlDv7pgiUoX0ldEhGR9zHpKxG5FZi7sOTvMp5niYjIl/qMEREREZF3MBgjIrciLh2TU2Oe1SUREXkfgzEicktHTUdFOoyeJiIi72EwRkRERGQhBmNEREREKS0Yi46O1nkpH3nkEcmWLZtt/fvvvy9nz571ZvmIiIiIUrUEBWPjxo2TMWPGSPfu3SU0NNS2vly5cjJ8+HBvlo+ILBaYq7Dkf+VrXRIRkY8EY5MnT5affvpJ3njjDYf1mCR88eLF3iobEfkA/8BgCcpTVJdEROQjwdipU6ekUqVKtuz7JkwSbl9TRkQpX+Sty3Jt5QRdEhGRjwRjxYoVk507d8YIxn7++WdtqiSi1CPqfqjc2btGl0RE5CPTIfXp00deeukl+fTTT/X2n3/+KatWrdK+ZFOmTPF2GYmIiIhSrQQFYz179pSIiAh55513dGRlw4YNJXfu3PL5559rkEZERERESRiMQe/eveXNN9/UVBYIyAoXLiz+/kxbRkRERJQswZjZXwxBGBGlXgGZskvWWu10SUREFgZjnTp18nin8+bN83jb+/fva3+zS5cu6QjNOnXqePzYI0eOyPLlyzX5bIMGDTx+HBF5Ll2W3JKjQReeMiIiq4OxzJkze/3JL168qEFUUFCQVKlSRT788ENp2bKlzJ49O87HhoWFSdu2beXYsWOafJbBGFHSiA6/JxGXjkpQvlLiH5yRp5mIyKpgbPr06d5+bhkwYIBkyJBBtmzZIunTp5e9e/dqLVfr1q2lVatWbh+LwQMIwNhPjShpPbhxXi79OEhCXh4nwSGleLqJiLwsUT3uIyMj5fjx43rB9fhAp3/kJevSpYsGYlC5cmVtpkR2f3fwuA0bNshnn32WmOITERERpcxgLDw8XGu1smfPLiVLltQLrg8aNEhTXnji9OnTcufOnRhJYnH74MGDbrP/Yxqm77//XmvVPC0vZgawvxARERGl2NGUb7/9tqxYsUK++uorqV69uq7bvn27DB48WG7evCmTJk2Kcx+3b9/WJYI4ezly5LDd5wy1b88995z069dP+5h5auTIkTJ06FCPtyciIiLy6ZqxH374wdbEWKFCBb3g+sKFC7XGyhOYxxKcAy/UWpn3OcO+UWsWEBCg2f5xuXLliuzevVuvG4bh8nEDBw6UW7du2S5nzpyJ9zETpVV+AekkIHMuXRIRkfcl6NsVwVKpUjE78pYuXTrWQMpZkSJFJDg4WEdDNm7c2LYet7EfV9AciqAPTZUmNIsigDt58qQGY/ZzZZrwPLgQUfwF5SkmhXp9y1NHRORLwVizZs208zzmpjSDHwRCWIf7PBEYGCjNmzfXWrYePXroqEj0I/vjjz9kxowZtu3Wrl2rKTBeeOEFqVu3rl7sYXuMqkTNGBEREVGaCMbu3bsno0aNkgULFki1atU0ENu1a5fWarVv3166devmUUqM0aNHS+3ateWpp56SmjVramBWv359ef755x0SyCL1BYIxIkp+EVdOyuWfhkjeDkO1loyIiHwgGEMtVseOHW23UTv26KOP6gUwStITaI7cv3+/9gVDBv6PP/5YM/2jT5gJTZixNVtC586dXTaZEpF3GFGREnXnmi6JiMhHgrH4THcUl3z58kmfPn1ivd8+6HMFIyuJiIiI0mTSVyIiIiJKnASPVd+4caNs2rRJbty4EeM+9CcjIiIioiQKxoYPH65JVJF41TlpKxGlLoE5Cki+5z7VJRER+UgwNnHiRFm1apU0atTI+yUiIp/iH5xR0hepbHUxiIhSrQT1GQsLC5NatWp5vzRE5HMib1+VG3/O1iUREflIMNaiRQv56aefvF8aIvI5UXdvSuiWhbokIiIfaab88ssvpWLFihqQYYoi5ymIMIE4ERERESVRMPbRRx9pYlc0V547dy4huyAiIiKihAZj8+fPl3Xr1sWYJ5KIiIiIkqHPWKZMmeThhx9OyEOJKIUJyJBVMlduoksiIvKRYKxhw4YyZ84c75eGiHxOumx5JddTb+mSiIh8pJkyIiJC3nzzTe3Aj0m6nTvwT58+3VvlIyKLRT8Il8ibFyVd9hDxDwy2ujhERKlOgmrGgoKCdALv/Pnzy927d7Uzv/2FiFKPB9fOyIWZvXRJREQ+UjM2b94875eEiIiIKA1KUM0YEREREVlYMwb37t2TLVu2yOnTpyUyMtLhvm7dunmjbERERESpXoKCsT179kjLli21f9jNmzclX758cunSJb2vSJEiDMaIUhEdoBOQLsZAHSIisrCZsk+fPvLiiy/KjRs39PbFixfl5MmTmgSWtWJEqUtQvpJStO8vuiQiIh8Jxnbu3Cl9+/bV6/hvGakuihYtKjNmzJBp06Z5u4xEREREqVaCgrFbt25Jzpw59XqePHls81Mi1cXly5e9W0IistSDq2fkwuy3dUlERD44mhJNk5g4HJ35+/fvLxUqVPBOyYjIJ0RHhkvEpWO6JCIiH+nA/8EHH9iujx49Wtq3by+PPfaYNlUyBxkRERFREgdjI0aMsF3HdEi7d++WsLAwSZ8+fUJ2R0RERJRmeSXpK9JbsK8YERERURIHY1evXpVhw4Y5rPvkk08kV65c2kRZtWpVTXNBRKkHJgjP3WqALomIyOJg7IsvvpCMGTPabu/bt08777/77ruyePFiCQgIkOHDhydBMYnIKgHpM0umcnV1SUREFvcZW7RokaxYscJ2e8mSJTp6csyYMXq7cOHC2pk/PlatWiUTJ07UDP6VKlWSIUOGSLFixWLdHmk0vvrqK/n7778lXbp0OpoTSWizZcsWr+clIs9E3b0hdw/8IZkqNJSATDl42oiIrKwZwzyUBQsWtN3evHmzNGrUyHYbgdmFCxfiFYg9/fTTUq9ePa11Q98zBFdYuhIVFSX169eX7Nmzaw1cv379ZNmyZdK4cWNNPEtE3hd5+5rcWD9Dl0REZHHNWIECBWTXrl1Su3ZtHT25adMm6dq1q+1+9BcLCfG8XwlqwTp16iQDBgzQ27Vq1dLHT5482bbOHppBDx48KMHBwbZ1mAsTQeC2bds0kCMiIiJKtTVjaILEnJRffvmlBlGYCqlJkyYONWWo5fIEJhnfvn27NG/e3LYOQRZq2tavXx/r4+wDMUBTJURHR8fnUIiIiIhSXs3Y4MGDtW8XmggxHdL3338vWbNmtd0/depUW/+xuKDvl2EYOoWSPdw+cOCAx2X6+OOPta9ajRo1Yt0mPDxcL6bQ0FCP909ERETkM8EYRlLOmjVLL664q9FyFhkZqcugoKAYNV8PHjzwaB8jR47UUZxr1651m3AW2w0dOtTjshHR/+cfnEkylKqhSyIi8tGkrwmB3GRw7Zpjp2Dczp07d5yPHzt2rOY8++WXX6ROnTputx04cKBObm5ezpzhhMdEngrMkV/yth2sSyIiSkXBGDrqY0DA1q1bHdaj31m1atXcPnb8+PEyaNAgrRVr2rRpnM+F2jY0p9pfiMgzRlSkRN27pUsiIkpFwRj07NlTpk+fLsePH9fbc+bMkf/++0+6devmMOKydevWttvISYaRlgjEmjVrZkm5idKSiCsn5ezEzrokIiIfmSjcW1C7dfLkSSlXrpw2Td67d09mzpwpVapUcejof+TIEb1+48YNefvttyVLlizSt29fvZjQZNmmTRtLjoOIiIgoRQZjSEuB4AsJXzHvJXKGOaeuQJCFIA3QvLh3716X+7JPRktERESUUlgajJly5MihF1fQr8w+6WvFihWTsWREREREqbjPGBEREVFa5xM1Y0Tku4LyFpfC7/wkfoGOXQiIiMg7GIwRkVt+/gHiF5yRZ4mIKImwmZKI3Hpw/Zxcmv+RLomIyPsYjBGRW9ER9yXs5G5dEhGR9zEYIyIiIrIQgzEiIiIiCzEYIyIiIrIQgzEicitd1jySs/FruiQiIu9jagsicisgYzbJUrUlzxIRURJhzRgRuRV1/7bcObBel0RE5H0MxojIrchbl+Tar1/okoiIvI/BGBEREZGFGIwRERERWYjBGBEREZGFGIwRkfsvicD0ElSgrC6JiMj7mNqCiNwKzFVI8r/4Bc8SEVESYc0YERERkYUYjBGRW+EXj8qp0S11SURE3sdgjIiIiMhCDMaIiIiILMRgjIiIiMhCDMaIiIiILMTUFkTkVlDuIlKgx1RJlyU3zxQRURJgMEZEbvmlC5LAHAV4loiIUnMz5ZkzZ2THjh0SGhqapI8hovh7cPOiXF02RpdERJTKgrGwsDBp27atlC1bVl588UUJCQmRiRMnev0xRJRw0WF35O7BP3RJRESprJly6NChsm3bNjl27Jjkz59ffvnlF2ndurXUqFFDatas6bXHEBEREfkqS2vGZs2aJd26ddOgCp599lmpWLGirvfmY4iIiIh8lWU1Y+fPn5dLly5JtWrVHNajhmv37t1eewyEh4frxXTr1i1dsr8ZUeyiw+/9bxkRZlua6/jZIaK0IPT/+qUbhpE6g7Hr16/rMleuXA7rcdu8zxuPgZEjR2rzprPChQsnqOxEadHlHwfYrmcbZ2lRiIiS1bVr1yRbtmypLxgLDAy0dci3d//+fQkKCvLaY2DgwIHSp08f2+3o6GgN3hDE+fn5eRQZI3DDCM6sWbNKWsPj5+vP9z8///z+4/d/Wvz9u3XrlhQpUkRy5syZpM9jWTCGL3d/f385d+6cw3rcxoF76zEQHBysF3vZs2ePd5nxRkyLb0YTj5+vP9///PynVfz+S9vff/7+/qmzA3/GjBmldu3asnTpUtu6u3fvytq1a6Vx48a2dUePHrX1B/P0MUREREQphaWjKUeMGKGpKdCMiAALIyPz5s0rPXr0sG0zatQozScWn8cQERERpRSWBmMNGjSQ9evXy6lTp2T8+PFSoUIF2bhxo2TOnNm2TenSpaVq1arxeoy3oYlzyJAhMZo60woeP19/vv/5+ef3H7//06LgZPr+9zOSerwmEREREfn23JREREREaRWDMSIiIiILMRgjIiIiSqsThVvp7NmzcvHiRR0g4GlWXU8ek5D9WuH06dNy+fJlKVOmjEe5Y6KiouTw4cOaeLd48eKSLp3jWwcpSHDc9jCookqVKuKLTp48KVevXpVy5crFOfjjv//+03NlD+escuXKidqvlU6cOKGJj1HOTJkyxbodupRu2rQp1rx/RYsW1euHDh3S47aH93+lSpXEF+F1wmf10UcflfTp03v0mGPHjsnNmzfloYce0jQ7Cd3GFxw/flynl8NUcu4SZpsiIiL0858lSxbN6eicc2n//v163PaQVBvnwRfhdbpw4YLUrFnTlkw8Nvv27bNNoWfKnTu3fnacHTlyRG7fvi3ly5f3+H1lBfP7+rHHHpOAgIBYt0OC9R07dri8r2TJkrY5ovfu3RtjijRkOcDvi68JCwvT9zKSuBYqVMijxO/mdxwSzGPQYGyfGU+2iZWRxoSFhRnt2rUzMmTIYDz00ENG+vTpjbFjxyb6MQnZrxXu3btntGrVysiYMaNRrlw5Le/XX38d6/bR0dHGsGHDjLx58+pxFS1a1ChcuLCxYsUKh+1effVVI3fu3EadOnVsly5duhi+5s6dO8ZTTz1lZMqUyShbtqyeh+nTp7t9TOfOnfX47Y+tR48eid6vFW7dumU0atTIyJw5s1GmTBldzpkzJ9btw8PDHY4blypVqmDQj/HFF1/Ytmvbtq0REhLisN2bb75p+Jp169YZTZs2NXLlyqXHcOTIkTgfc/36daNBgwZGlixZjNKlSxtZs2Y15s2bF+9tfMHq1av19c+ZM6ce/5kzZ9xuf/fuXaNPnz66faVKlfQ1xvfA1q1bHbZ78sknjUKFCjm8/gMHDjR8zcqVK40nnnjCdvwXLlyI8zF4XfGdZ39sH374ocM2ly5dMmrVqmVkz57dKFWqlJEjRw5jyZIlhq9ZtmyZHg/Kh+O/ceOG2+3Pnj0b4/Nfvnx5fez8+fNt29WsWdMoUqSIw3b43fAl169fN1577TV9jSpXrqy/V1WrVjUOHDjg9nGnT582Hn74Yf3OKF68uJEnTx5j7dq18d4mLmkuGMOHqECBAvomM9+ceGNt2rQpUY9JyH6t0LdvX/3QXLx4UW8vWLDA8PPzM3bs2BHrj/FHH32kb2QzOMOxIujAF5B9MNaxY0fD1yFAKFGihHHlyhW9PXfuXMPf39/Yt2+f22Ds5Zdf9vp+rdCtWzcNFs3Xc9q0aUa6dOmMw4cPe7yPCRMm6GPM95AZjPXs2dPwdePHj9d/JPC59DQYe+GFF4yKFStqIAsTJ040goKCjBMnTsRrG1+AABoB2e+//+5RMIYfGTwGQRk8ePBAPwv58+c3IiMjHYKx999/3/B1n3/+ufHbb78Za9asiVcw9sEHH7jd5tlnnzUeffRR23kaOXKk/kPmyf6T06hRo/S1N3+f4grGXOnfv78Gc/fv33cIxoYPH274sgMHDhjffPON/qaZFSjPPPOM/nPhTsOGDfViPg7vcxy//bnzZJu4pLlgDAGT8381+E8fwURiHpOQ/SY3BFKI3EeMGOGwHm/GXr16ebwfBJz4IONL3YTjRI0bgrpjx44ZUVFRhq/BDwlqLMaMGeOwHv/JvPfee26DsQ4dOuixHT9+PMaxJXS/yQ1fPviB+OqrrxzeE3jvxvVjYw//AbZp08ZhHYKxF1980di+fbtx8uRJ3a8v27x5s0fB2O3btzWosq/lRBCC/6rNHx9PtvE169ev9ygYc2Xjxo36WPsAHsHYG2+8oa8/Ajhff/0RkMUnGOvdu7exbds2l+cL/4DhHy/7mlAEKqgl9cXWEUhoMIbvOtSOvvXWWw7rEYzhH32cI7NCIiX45Zdf9Dxcu3bN5f34vsf9q1atsq3DOQsMDDRmzZrl8TaeSFMd+NHvB/0kqlWr5rAe/SbMKZcS8piE7NcKmOgXM887l7N69erxKuf27dttfQbsrVixQrp27Sq1atXSfmWrV68WX+sng34NCTn+xYsXyyuvvKLblipVSn7//Xev7Dc5oZ/EvXv3HMqJ/hLoN+VpOXfu3Cl79uyR7t27x7hv/vz58uqrr2qSZvQV+euvvySlO3jwoPaXsj9n6GODYzTPmSfbpCb4/KM/DPoM2ps5c6Z069ZN+1Kir6D5PZEaTJs2Td/zFStW1OPbtWuX7T70l4qOjnZ4/dFfDOcgtb3+y5cv175meJ2dTZo0Sc8R+svhvY/z4uu2b9+ufcdy5Mjh8n7z9bN/bTGvNfqEm/d5so0n0lQwhg7LZsdSe7ht3peQxyRkv1bwRjkReL799tvy/PPPOwRjzZs31wnb8QFEYNquXTtp27atdpRO6cePKbfwBYQgBJ1+caytW7fWDuCJ2W9y80Y5Z8yYoR24mzRp4rC+Q4cOcunSJds5atiwoZ4jrEvJUtPn31sB/ccffyx9+/aVDBky2NbjH5UrV67IP//8o59/BCL43Dh3fE+JEGBgcIp5bGXLlpVWrVppR/209vrj849/tp0H5vTq1ct2jvA7gEAdrz/++fNVO3bskC+++EI+/PDDWDvxm68fAra4Pv/utvFEmgrGzFEzGE1hD6MfYhv54MljErJfKyS2nDdu3JBmzZrpCJSpU6c63NemTRvJkyePXsdIy9GjR+sbfNmyZZLSjx+BpflBwz7wAcY+Vq5cmaj9JrfElhPb/fDDD1r75TyaDsEY/hsE7GvcuHH6Q7xmzRpJyVLT598bI7CbNm0qjRo1kmHDhjnch3/OzNHDCNLw+iNwSQ21o507d7aNOMYI2bFjx+o/YpiGLy29/vgnC995rmrFMX+0GZzjfYDvSIzY3rp1q/iif//9V1q0aCEvvPCCvPPOO7FuZ7624eHhcX7+3W3jiTQVjBUsWFCbDxC528Nt/Lef0MckZL9WQFkQICWknBi2jtoQVL+vWrXKbToEMyBDAOP8XFYy0zAk9nXCHGUIPMz9eGu/SS2x5fz555/lzp07WgsSF7w/8KXsS69/QnhyzlLK65/YQAy1nY888oj8+OOPbtMhmLUCrr4TUwOkbLD/Hk0Lrz98++23GnB17Ngxzm3z5cunS198/Q8dOiRPPPGEBmOoVHCX2iK21xb/aMT1+bffxhNpKhhDIFGvXj1ZunSpQ/T622+/SePGjR1y0Jh9Ajx5jKf7tRpyBCGvjn05UdWO/k/25USuHFQ3m1DDgUAMARYCMezHHvpL4Hid3/Doo4Y+Fr4CuYGQ98z++FHbh//e7Y8fTTFmf4fIyMgY//GgKQ7Nteaxebpfq6FGE7mR7MuJ49i8ebNDOfFfI3IrOZs+fbqtZtTegwcPtM+UPfxHjPeNL73+njpw4IBeAP0D0f/R/pyhVgR958xz5sk2KQlee7wH7I/l8ccfl4cfflh++umnGHm58PlAHkJ769at03Up8fXH8eP7y6ztcj42fK9j8Jt5bGiyQ/Bh//ojNyHOYUp8/fH9hvI7Q59A1IA6/yOO7378Btgza8R97fU/fPiwvpfxPYbvM1eBGPp5IQ8b4PcSv3f2ry36mSHQMl9bT7bxiJHGYCQQRjn069dP88A0adLEKFasmG1IOmCIPob/x+cxnmzjK3mWkJZg0KBBOpIEOXeQFwl5skwYuo4Rc4ChusifgxE0y5cvNzZs2GC7mKktkLsMuWe+/PJLHVEyZcoUTZ+BETbmUF9fgbQGAQEBxpAhQ/T469Wrp2W3H6aNFB0oO9y8eVNTFowbN06PDUOjCxYsqI/DyKL47NcXLF68WMuJUX64/thjj+lrHRERYdsGo2IxgsweRh0iBQqOzRneB8hBhZQXGGGL0Zp4v2CEna+NqsVIP7x38R7F1x9GwOG2fZoO5CHDxYRt8JlBuoJFixZpCoPq1as7pHbwZBtfcOrUKT1evFY4fpQVty9fvmzbBq893gNw9epVzZuFnIRIiWD/+cdnA44ePar5miZNmqSvP0YQYtR269atDV+Dkb4oO8qI48d3NW6bKWkAObIwOhgOHTpkVKtWTT/3ODak+UCOsvbt2zvsd8aMGTqiFt+BCxcu1M8DzqOvjSrFyD8c7+jRo/X4kXcNt+1HE+J4MYLc3h9//KHbu0qBtGfPHqNGjRrG5MmT9Rx99tlnmssLo6t97bNfoEAB/Wz++eefDu9lMyUJVKhQwSELAl5zjEJHPk7kVsPvJVJi2PNkm7j44Y+kMagJ+Oqrr7RTNiL3AQMG2DIJw5gxYzQ6/v777z1+jKfb+IINGzboyBfUiuC/XZQTVe+mTz/9VGvHZs2apR0Qn3nmGZf7+eCDD+Spp56yVdHi2FGjhpEpDRo00L5Fzpn6fcH69etl8uTJOrIUzS7vv/++1m6Z0EEZ/SOmTJmit1HDN3HiRK0tQ/MLqri7dOkSo6kmrv36Cvxnj9FheG0xkhLltB9NNGjQIK0xxTGb5s6dq5+HX3/91eVrioEaeP2RiR3HjH5FL730Uoy+ZVabM2dOjP6OgM9Ay5Yt9To6p5vfAyb0lUHnZTTX4z/h/v37x5hhw5NtrIbXHc1NzgYPHmwblNG7d2/9Tx/fAxgp2qNHD5f7wvsD73NATQq+U1AbFBISot8LaM7yNLt5csHn87vvvouxfujQofLkk0/q9TfeeEM/58OHD9fbqCX75ptvdIljwwAeV011v/zyi76/8NmpU6eO9OvXL87uHMltwoQJWrvpbOTIkdq6Axgpiaa3jz76yHY/+oChwzuaqF3B5x7nFu+DAgUK6G8G+hH7kr///ls/k67g+w2124DvLWTQx/eiCcc9b948rQVEzdq7774bY4YFT7ZxJ00GY0RERES+wrf+bSUiIiJKYxiMEREREVmIwRgRERGRhRiMEREREVmIwRgRERGRhRiMEREREVmIwRgRERGRhRiMEVGC5ytcvHixz549JOncsmWLV8uKKU4WLFgQYz0SYmI9lr4OU70g8TMR+Q4GY0Tk0tWrVzWjNC7z58/X+QYxa4N9Ruvu3bv75NnDnKLIAI5s+N4sK+asffHFFx3W9ezZU58LwRjm9XOGeTtxDs2yWA1Z4Vu3bi2nTp2yuihE9H98b64aIvKZmqXnnntOWrRoIZkzZ9ZpoTD5NabJ6dOnj/iy2bNnS4YMGXRC4KSESaQxgTImRsYUKK6EhobqecQUa5hQ3mqYoq1Dhw4ybNgwnb6JiKzHYIyI3Bo3bpyUKlXKNkcd5tx79tlnbfdHR0fr3HRoCsScrMWKFbPdd+PGDVm9erVeDw4OltKlS+s2zm7duiXbt28XzM6G+TLt58o0mwfRBIi5HqtWrapzJ7qDeTIxN6ozb5TVhDloMVcnauG2bdsmly5d0nlL7ed5BbN5FPtGgIu5O3GMq1at0hq1oKAgW9kwbyDm9cQ2mONuyZIlGgyjvIcPH9a5IDHfJ+aAffrpp2Xfvn0aJFeqVEnnE7SHmjiz2RTPlz17dtt9qN1DWTH/pvO5JqLkx2CMiDyGIAwTaSMIgPDwcA0esMyYMaP89ddfOhk1Jts1AxxMoAwILtCHC4EB1gUGBur6P/74Q5vNMDlv1qxZNWBB0Id15sTt48ePl1q1asmdO3d0MmJMyotgwpWzZ89qc6Hz/d4oqz0EXytWrLBNEo9gp3LlyjGCMUwgbr8NgjwEW6gtu3Llim0yeTRnYh36c9WtW1cnnDdrJo8cOaL7xgTWCL569eqlQemDBw90MmI8BjV0nTt31n2hSblt27YapKFWE4EcgmpM4Aw1atTQCd+xXbt27fgJILIaJgonInK2YcMGA18RR44csa1buXKlrtu6davx448/6vXJkyfb7h81apRRoECBWE9maGioUaZMGYfHNGvWzOjTp4/t9u3bt401a9bo9UWLFhkhISHGmTNnbPfjsXiO8PBwl8+xZMkSw8/Pz4iMjLSt81ZZly1bZgQHB9tuX7lyRfe7e/fuWPfjahtcxzrcZ7p//76uw3kHHDNut2nTxuFY5s6dq+unT59uWzdixAijSJEittuNGjUy+vfv73Asv/32m0O5qlev7rANEVmHHfiJyK3ly5drB3TUVr3yyivaNwo1RoAao27dutm2bdiwoTYpogbLhOY3NEGiuQ77KlKkiDbrmdC36/jx47bHoCancePGen3WrFnaVIhaKnSQRzMemvXwHKjtiW3gAWrYAgICHNZ7o6xWQC2Y87GgNqxr164Ox4KmTNTomef02LFjcvfuXb2NZl3UCtpD8yTOFRFZj82UROTWb7/9pj/mOXPmlE8++UReeOEF7bcEzkEP+lpBWFiYBlUnTpyQJk2aaL+q8uXL634QNCCYMI0ePVqDJDTvoSkSzXKvv/66NiWePHlS+5EtXLjQoUwdO3YUPz8/l+XF85pBiD1vlDUuGzdu1GZSU6dOncQbHe6doe+c+RrYHwuaYBGIoS8Y+syZ57Rly5Y66hPn1IRzFFffOyJKHgzGiMjjDvzxNXz4cClTpozWMpkQzNmneUAfqj///FNraX7//XcN+NDZHSMUEUDh8egP5Slsj4AKQVGhQoW8Wta4bN26VWvW4grGzEAKNXEmBIWuxBZ0uoPjQD8y9EnDOR0xYoQG1WYfN0Cga/YxIyJrsZmSiJIMRhyWLVvWIc0DggJ7586d0yU6siPlwoABA2zJWpGaYtGiRQ75zewf4wo6uufJk0c2bdrk9bLG5b333rPlZsMFUOvmHGwVLFjQloDVhA7+3mKeH5wH1CL279/fdk7NQAzbODddEpE1WDNGREk6+hIBCoIC1HJNmTJFRw3aQ3MamkDr1aunebsmTpwo7du31/uQzww1VdWrV5c33nhDRyMi19nmzZttIzpd1TqhbxtGXCIQ8WZZEwLNnBgp+tlnn+kIx3z58mkQ1Lx5c01E+/bbb+vIzLlz54q3vPzyyxISEiJ16tTRc4rRqOY5BSTxrV+/vtZKEpH1GIwRkUtmrUps/YqQ1wp5spw7heMxZh+mHj166Do0OSJIGjJkiHaYR/BhQrCFzvlINYEmOaSyMAMHMwUFAiv0x0Kfr5o1a2rA5s4777yjARA6sZcsWdJrZUWNFmrvTHgs9hFXrq6lS5dqcIc0F8WLF9dgDMeMfGho2kQZkeIDwSfOu3ns2DcCQ3vIjWaf5w0QzGJbM2cZmnkRcKGpEseCJljznCIdBsqCCxH5Bj8MqbS6EERE3mbWNDlPX5TWIfhDsloEaETkGxiMEREREVmIHfiJiIiILMRgjIiIiMhCDMaIiIiILMRgjIiIiMhCDMaIiIiILMRgjIiIiMhCDMaIiIiILMRgjIiIiMhCDMaIiIiILMRgjIiIiMhCDMaIiIiIxDr/D5CETcc7uRRbAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 600x280 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axis = plt.subplots(figsize=(6, 2.8), layout=\"constrained\")\n",
    "axis.bar(list(rotation_phases), list(rotation_phases.values()), width=0.8 * (2 / 2**n_phase))\n",
    "axis.axvline(target_phase, color=\"black\", linestyle=\"--\", linewidth=1,\n",
    "             label=\"Known phase = 0.75\")\n",
    "axis.set(title=\"Single-qubit rotation\", xlabel=\"Phase (half-turns)\",\n",
    "         ylabel=\"Sample probability\", xlim=(-0.05, 2), ylim=(0, 1.1))\n",
    "axis.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "qpe-flow-09",
   "metadata": {},
   "source": [
    "- The phase lies exactly on the grid, so the ideal simulator returns it every time. With a superposition of eigenstates, repeated runs instead sample their phases according to the input's overlaps.\n",
    "\n",
    "- `measurement_dataframe` decodes the raw register counts into phases and empirical probabilities; the following examples use the same helper.\n",
    "\n",
    "## 2. Estimate an energy with Trotterized evolution\n",
    "\n",
    "Use the noncommuting Hamiltonian\n",
    "\n",
    "$$\n",
    "H=\\tfrac12X+\\tfrac12Z,\\qquad E_0=-1/\\sqrt2.\n",
    "$$\n",
    "\n",
    "- Its ground-state energy is $E_0=-1/\\sqrt2$.\n",
    "- Its known ground state is prepared with $R_y(-3\\pi/4)|0\\rangle$.\n",
    "- For larger problems, finding a state with useful eigenstate overlap is a separate task; QPE does not prepare the ground state.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "qpe-flow-10",
   "metadata": {},
   "outputs": [],
   "source": [
    "from guppyalgos.algorithms.time_evolution.trotter import (\n",
    "    ham_sim_trotter, trotter_first_order, cntrl_trotter_first_order,\n",
    ")\n",
    "\n",
    "hamiltonian = zqp.RealTermSum.from_str(\"(0.5, X0), (0.5, Z0)\", 1)\n",
    "controlled_step = cntrl_trotter_first_order(hamiltonian, 1)\n",
    "total_time = 1.0\n",
    "n_steps = 16\n",
    "exact_energy = -1 / np.sqrt(2)\n",
    "\n",
    "# Full, uncontrolled evolution over total_time.\n",
    "simulate = ham_sim_trotter(\n",
    "    trotter_first_order(hamiltonian, 1),\n",
    "    n_steps=n_steps, time_step=total_time / n_steps, n_state_qubits=1,\n",
    ")\n",
    "simulate.compile_function()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "qpe-flow-11",
   "metadata": {},
   "source": [
    "### Prepare the state and controlled evolution\n",
    "\n",
    "Split $t=1$ into $r=16$ first-order Trotter steps:\n",
    "\n",
    "$$\n",
    "U(t)\\approx\n",
    "\\left[\\prod_j e^{-i\\pi(t/r)h_j/2}\\right]^r.\n",
    "$$\n",
    "\n",
    "- `ham_sim_trotter` builds the complete uncontrolled evolution.\n",
    "- QPE needs controlled evolution, so `trotter_power` repeats `cntrl_trotter_first_order`.\n",
    "- A requested power repeats the complete 16-step simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "qpe-flow-12",
   "metadata": {},
   "outputs": [],
   "source": [
    "@guppy\n",
    "def prepare_ground(qreg: array[qubit, 1]) -> None:\n",
    "    ry(qreg[0], angle(-0.75))  # Angles are in half-turns.\n",
    "\n",
    "\n",
    "@guppy\n",
    "def trotter_power(control: qubit, qreg: array[qubit, 1], power: int) -> None:\n",
    "    for _ in range(power * n_steps):\n",
    "        controlled_step(control, qreg, total_time / n_steps)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "qpe-flow-13",
   "metadata": {},
   "source": [
    "### Build the main circuit\n",
    "\n",
    "- Pass the prepared system register and Trotter oracle to the same `qpe` routine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "qpe-flow-14",
   "metadata": {},
   "outputs": [],
   "source": [
    "@guppy\n",
    "def trotter_qpe() -> None:\n",
    "    qreg = qarray(1)\n",
    "    prepare_ground(qreg)\n",
    "    phase_qreg = qarray(n_phase)\n",
    "    transversal(h, phase_qreg)\n",
    "    qpe(phase_qreg, qreg, trotter_power)\n",
    "    output(\"phase\", collect_measurements(measure_array(phase_qreg)))\n",
    "    discard_array(qreg)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "qpe-flow-15",
   "metadata": {},
   "source": [
    "### Convert the measured phase to energy\n",
    "\n",
    "For an energy eigenstate, the simulated evolution contributes\n",
    "\n",
    "$$\n",
    "U(t)|E\\rangle\\approx e^{-i\\pi tE/2}|E\\rangle.\n",
    "$$\n",
    "\n",
    "QPE reports $e^{i\\pi\\phi}$, so\n",
    "\n",
    "$$\n",
    "\\phi=-\\frac{tE}{2}\\pmod 2.\n",
    "$$\n",
    "\n",
    "For the negative-energy branch at $t=1$, use $E=-2\\phi$. The result is\n",
    "\n",
    "| Quantity | Exact | QPE result |\n",
    "| --- | ---: | ---: |\n",
    "| Phase resolution | -- | $1/32$ |\n",
    "| Phase $\\phi$ | $1/(2\\sqrt2)=0.353553\\ldots$ | $11/32=0.34375$ |\n",
    "| Energy $E=-2\\phi$ | $-1/\\sqrt2=-0.707107\\ldots$ | $-0.6875$ |\n",
    "\n",
    "Choose the evolution time and energy range so the energy branch can be recovered without ambiguity from phase wrapping."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "qpe-flow-16",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bitstring  counts  phase  empirical_probability  energy\n",
      "   110100      81 0.3438                 0.6328 -0.6875\n",
      "   001100      26 0.3750                 0.2031 -0.7500\n",
      "   010100       6 0.3125                 0.0469 -0.6250\n",
      "   101100       5 0.4062                 0.0391 -0.8125\n",
      "   010010       2 0.5625                 0.0156 -1.1250\n"
     ]
    }
   ],
   "source": [
    "trotter_counts = (\n",
    "    trotter_qpe.emulator(n_phase + 1).with_seed(5).with_shots(n_shots)\n",
    "    .run().register_counts()[\"phase\"]\n",
    ")\n",
    "trotter_results = measurement_dataframe(trotter_counts)\n",
    "trotter_phases = trotter_results.set_index(\"phase\")[\"empirical_probability\"].to_dict()\n",
    "trotter_energy = lambda phase: phase_to_energy_qpe(phase, total_time)\n",
    "dominant_phase = max(trotter_phases, key=trotter_phases.get)\n",
    "trotter_results[\"energy\"] = trotter_results[\"phase\"].map(trotter_energy)\n",
    "print(trotter_results.head(5).to_string(index=False, float_format=lambda value: f\"{value:.4f}\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "qpe-flow-17",
   "metadata": {},
   "source": [
    "- There are two sources of error: the finite phase grid and the Trotter approximation.\n",
    "\n",
    "- More phase qubits refine the grid; more Trotter steps improve the simulated evolution.\n",
    "\n",
    "- QPE measures the implemented unitary, so refining the grid alone does not remove Trotter error.\n",
    "\n",
    "## 3. Estimate the same energy with qubitization\n",
    "\n",
    "Instead of approximating time evolution, encode $H/\\lambda$ with LCU. For this Hamiltonian,\n",
    "\n",
    "$$\n",
    "B=P^\\dagger SP,\\qquad \\lambda=\\sum_j|a_j|=1.\n",
    "$$\n",
    "\n",
    "`PREPARE` loads the LCU weights, `SELECT` applies the indexed Pauli term, and `UNPREPARE` completes the block encoding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "qpe-flow-18",
   "metadata": {},
   "outputs": [],
   "source": [
    "from guppyalgos.algorithms.block_encoding.lcu import (\n",
    "    LCUCntrl, LCUData, build_cntrl_single_cntrl_select,\n",
    ")\n",
    "from guppyalgos.algorithms.block_encoding.qubitization import QubitizationCntrl\n",
    "from guppyalgos.algorithms.state_preparation import multiplexor_prep\n",
    "from guppyalgos.primitives.subroutines.reflection import ReflectionCntrl\n",
    "from guppyalgos.primitives.gate_decompositions.cnx.cnx import cnx\n",
    "\n",
    "data = LCUData.from_hamiltonian(hamiltonian)\n",
    "prepare = multiplexor_prep(data.amplitudes)\n",
    "controlled_select = build_cntrl_single_cntrl_select(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "qpe-flow-19",
   "metadata": {},
   "source": [
    "### Define the controlled walk\n",
    "\n",
    "Build the walk by following the LCU block with a reflection about the zero preparation state:\n",
    "\n",
    "$$\n",
    "W=(I-2|0\\rangle\\langle0|)B.\n",
    "$$\n",
    "\n",
    "- The two reflections act as a rotation in each invariant subspace.\n",
    "- Each energy therefore produces conjugate walk eigenvalues $e^{i\\theta}$ and $e^{-i\\theta}$.\n",
    "- Their QPE phases are $\\theta/\\pi$ and $2-\\theta/\\pi$.\n",
    "- `walk_power` exposes repeated controlled walks through the QPE oracle interface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "qpe-flow-20",
   "metadata": {},
   "outputs": [],
   "source": [
    "@guppy\n",
    "def unprepare(prep_qreg: array[qubit, 1]) -> None:\n",
    "    with dagger:\n",
    "        prepare(prep_qreg)\n",
    "\n",
    "\n",
    "@guppy\n",
    "def controlled_walk(\n",
    "    control: qubit, prep_qreg: array[qubit, 1], qreg: array[qubit, 1],\n",
    ") -> None:\n",
    "    QubitizationCntrl(\n",
    "        LCUCntrl(prepare, controlled_select, unprepare), ReflectionCntrl[1](cnx)\n",
    "    ).compose(control, prep_qreg, qreg)\n",
    "\n",
    "\n",
    "@guppy\n",
    "def walk_power(\n",
    "    control: qubit, regs: QubitizationRegs[1, array[qubit, 1]], power: int,\n",
    ") -> None:\n",
    "    qubitized_power_oracle(control, regs, power, controlled_walk)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "qpe-flow-21",
   "metadata": {},
   "source": [
    "### Build the main circuit\n",
    "\n",
    "- `QubitizationRegs` groups the LCU preparation qubit with the system register.\n",
    "- The QPE interface is unchanged: pass the phase register, grouped target registers and power oracle to `qpe`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "qpe-flow-22",
   "metadata": {},
   "outputs": [],
   "source": [
    "@guppy\n",
    "def walk_qpe() -> None:\n",
    "    qreg = qarray(1)\n",
    "    prepare_ground(qreg)\n",
    "    regs = QubitizationRegs(qarray(1), qreg)\n",
    "    phase_qreg = qarray(n_phase)\n",
    "    transversal(h, phase_qreg)\n",
    "    qpe(phase_qreg, regs, walk_power)\n",
    "    output(\"phase\", collect_measurements(measure_array(phase_qreg)))\n",
    "    discard_array(regs.prep)\n",
    "    discard_array(regs.target)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "qpe-flow-23",
   "metadata": {},
   "source": [
    "### Convert walk phases to energy\n",
    "\n",
    "QPE now measures a **walk phase**. Convert either member of the conjugate pair with\n",
    "\n",
    "$$\n",
    "E=-\\lambda\\cos(\\pi\\phi).\n",
    "$$\n",
    "\n",
    "- Here $\\lambda=1$.\n",
    "- Decode both peaks and convert each phase to its Hamiltonian energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "qpe-flow-24",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bitstring  counts  phase  empirical_probability  energy\n",
      "   000100      75 0.2500                 0.5859 -0.7071\n",
      "   000111      53 1.7500                 0.4141 -0.7071\n"
     ]
    }
   ],
   "source": [
    "walk_counts = (\n",
    "    walk_qpe.emulator(n_phase + 2).with_seed(5).with_shots(n_shots)\n",
    "    .run().register_counts()[\"phase\"]\n",
    ")\n",
    "walk_results = measurement_dataframe(walk_counts)\n",
    "walk_phases = walk_results.set_index(\"phase\")[\"empirical_probability\"].to_dict()\n",
    "walk_energy = lambda phase: phase_to_energy_qubitized_qpe(phase, data.l1_norm)\n",
    "assert set(walk_phases) <= {0.25, 1.75}\n",
    "for phase in walk_phases:\n",
    "    np.testing.assert_allclose(walk_energy(phase), exact_energy, atol=1e-10)\n",
    "walk_results[\"energy\"] = walk_results[\"phase\"].map(walk_energy)\n",
    "print(walk_results.to_string(index=False, float_format=lambda value: f\"{value:.4f}\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "qpe-flow-25",
   "metadata": {},
   "source": [
    "## From observed peaks to energy\n",
    "\n",
    "- The horizontal axes below show **phase**, not energy. Convert each peak using the operator that QPE measured:\n",
    "\n",
    "| Operator | Observed phase (half-turns) | Conversion | Estimated energy |\n",
    "| :-- | :-- | :-- | :-- |\n",
    "| Trotterized evolution, $t=1$ | $0.34375$ (dominant bin) | $-2\\times0.34375$ | $-0.687500$ |\n",
    "| Qubitized walk, $\\lambda=1$ | $0.25$ | $-\\cos(\\pi/4)$ | $-0.707107$ |\n",
    "| Same walk, second peak | $1.75$ | $-\\cos(7\\pi/4)$ | $-0.707107$ |\n",
    "\n",
    "- The exact ground energy is $-1/\\sqrt2\\approx-0.707107$. Both walk peaks lie exactly on the phase grid and recover it.\n",
    "- The input overlaps both walk eigenstates; finite sampling makes their peak heights slightly different.\n",
    "- The ideal evolution phase, $0.353553\\ldots$, lies between bins. Its dominant measured bin gives an energy about $0.0196$ above the exact value. With 16 Trotter steps, this gap is mainly phase resolution.\n",
    "- The check below allows for both grid resolution and Trotter error. Increasing the phase-register size reduces the grid spacing; increasing the step count reduces Trotter error.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "qpe-flow-26",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trotter: -0.687500\n",
      "Qubitized: -0.707107 (both peaks)\n",
      "Exact: -0.707107\n"
     ]
    }
   ],
   "source": [
    "from guppyalgos.tests.helpers import get_unitary\n",
    "\n",
    "phase_resolution = 2 / 2**n_phase\n",
    "exact_phase = np.mod(-total_time * exact_energy / 2, 2)\n",
    "\n",
    "# Find the Trotter eigenphase closest to the exact ground-state phase.\n",
    "implemented_unitary = get_unitary(simulate, 1)\n",
    "implemented_phases = np.mod(\n",
    "    np.angle(np.linalg.eigvals(implemented_unitary)) / np.pi, 2,\n",
    ")\n",
    "trotter_phase = min(\n",
    "    implemented_phases, key=lambda phase: abs(phase - exact_phase),\n",
    ")\n",
    "nearest_phase_bin = round(trotter_phase / phase_resolution) * phase_resolution\n",
    "np.testing.assert_allclose(dominant_phase, nearest_phase_bin)\n",
    "\n",
    "# Allow for Trotter error and one phase-bin width.\n",
    "trotter_error = abs(trotter_energy(trotter_phase) - exact_energy)\n",
    "np.testing.assert_allclose(\n",
    "    trotter_energy(dominant_phase), exact_energy,\n",
    "    atol=phase_resolution / total_time + trotter_error, rtol=0,\n",
    ")\n",
    "\n",
    "# Both qubitized peaks map directly to the exact energy.\n",
    "qubitized_energies = [walk_energy(phase) for phase in walk_phases]\n",
    "np.testing.assert_allclose(qubitized_energies, exact_energy, atol=1e-10)\n",
    "print(f\"Trotter: {trotter_energy(dominant_phase):.6f}\")\n",
    "print(f\"Qubitized: {walk_energy(next(iter(walk_phases))):.6f} (both peaks)\")\n",
    "print(f\"Exact: {exact_energy:.6f}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "qpe-flow-27",
   "metadata": {},
   "source": [
    "## What changed?\n",
    "\n",
    "| Example | Controlled operation | Readout |\n",
    "| :-- | :-- | :-- |\n",
    "| One-qubit rotation | Multiply a rotation angle by the requested power. | Phase $\\phi$. |\n",
    "| Trotterized simulation | Repeat the controlled Trotter step. | Energy from $-2\\phi/t$, on the chosen phase branch. |\n",
    "| Qubitization | Repeat a controlled LCU-and-reflection walk. | Energy from $-\\lambda\\cos(\\pi\\phi)$. |\n",
    "\n",
    "- These are different oracles for canonical QPE. Each run uses $m$ phase qubits and asks for powers $1,2,\\ldots,2^{m-1}$. The two histograms below use the same Hamiltonian and ground-state input. Trotterized evolution has one dominant energy phase; the reflected walk has two conjugate phases for that same energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "qpe-flow-28",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA48AAAFLCAYAAAB2j7mEAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAigNJREFUeJzt3QeUU1X3NvANQ++9V+kd6SgdQaQrICpVikhTpAhIl66IoKhIb9KRKiBdehGQ3nvvvZf7rWf/v5s3k0kmmZnMpD2/tYaQ5CY5uUnuufuUfaIZhmEIERERERERUSiih3YnEREREREREYNHIiIiIiIicgl7HomIiIiIiMgpBo9ERERERETkFINHIiIiIiIicorBIxERERERETnF4JGIiIiIiIicYvBIRERERERETjF4JPKQcePGSYsWLaL8dRcsWCD169eXp0+fSqDtF0/tcyIif7NhwwatSy5fvux02127dum2p0+fdum5w7q9u+D18Lrbt28XT7L3/sOyv4kiE4NH8kp//fWXHiRd+YtIEPTLL79I69atw3yfO+zZs0eWLFkiUe3IkSMaQL58+VK8UUT3S2ifm6f2ORGRJ61bt066d+8ujRs3llatWslPP/0kV69ejdBznj17VuuS+/fvO9320qVLuu3t27ctt23ZskXrcDyPK9tHBbweXvfixYviSfbef1j2N1FkihGpz04UTrly5ZKPPvrIch0BYpMmTaRs2bLyxRdfBNs2ZsyYEWrdW7lyZZjvc4c2bdpIjRo1Iu35A1Vonxv3OREFkhs3bkjDhg31uPjZZ59JtWrVNPiYN2+e9OzZU3788Ue9PbKVKFFCXzNbtmyW2y5cuKDBUO/evV3anoi8A4NH8krZs2fXP9PDhw/1MlOmTNpS6Q/efPNN/SPucyIid3v+/Lk2UKLH6t9//9VGWVO7du2kR48e2qCWOHFiDTAjU7p06cJUd4d1eyKKOgweyWft3r1bhg4dKoMGDRLDMGTSpEly/vx5GT16tKRJk0bu3r0rU6ZMkf/++0/vL1KkiDRv3lwrSujXr5/OIcB21pXUtGnTZPjw4Q7vixcvnqVinjNnjmzatEkePXokOXLkkJYtW0rGjBkt248aNUqOHz8uv/76q8yfP1/+/vtvSZkypQwZMkTn32FeBcoNBw4ckAEDBth9r0FBQfpaJlde2xz6MnbsWDlz5owG4zhhCAtnr4Oy4z1NnjzZsl9Ma9as0dceNmyYpSHA2WfiSJ8+ffRy4MCBwW6fO3eu/s2aNUt7oEP7TFE+231ucqVc1p8lWssxtDpGjBhSr149effdd8O0X4mIIhuOaehxnDhxYrDA0TR48GBZtGiRdO7cWd5//32JFSuWHttwPJ0xY4bEiRPHsu2JEye0pxJDX4sXLx7iuXbu3ClTp06VBw8e6AihTz/9VI+PJpQD9ep3330nb7zxhr4ujqnQrVs3y7H2q6++krfffjvE9oAA99WrV3bfK479BQoUsFxfv369LF26VIfmpkqVSj788EN56623gj3m2bNnMmHCBK0TkiZNqmV2xa1btzToxl+VKlX0NtQbGB2FOhP7DvsSUHfinAQ9vKg3zf0I0aJF032Mz6ZRo0aSNWtWCY8ffvhBduzYoc/LBmmKCpzzSD7rypUrWtEtXrxYOnTooAfehAkTai/l0aNHJW/evPL7779L4cKF9YD622+/Sf78+fXgDZUqVdLHxI0bV4fImn8IQkK7D65duybFihXTlltsV7lyZdm7d6/ky5dPNm/ebCkjKqXly5dLly5dNJjCUBxzDoPt/LvUqVMHey38Va9eXVasWKF/JldfG8FowYIFZdmyZVK6dGmJHz++fPDBB1rxucKV10GljuANw4tsjRgxQk8AzIrflc/EEZwI4M/W4cOH9TtgnlA4+9zszXl0tVzmZ4lAdu3atXqC9PjxYx0GZh3YExF5AxyXEMBZTwGxbZT85JNPNAGLeUx3NCce9QZuR4OkLdQxCN4KFSqkw0xR31WtWlVevHjhcA5fnjx5tF4CNL6Zx+rMmTPb3d4MHm3ryGPHjul2CFrNIA5BII7Lr1+/1nLgfVaoUEGDUdOTJ0+kXLlyejzPnTu31msIYjE31JnkyZNrsDZ9+nTLbWh4/OOPP7QuxFxO0+zZs7W+QE+q+Viz7Hg/CGjR4InXx3OEBQLVpk2b6ntAcMzAkaKMQeQDHjx4YODr2qhRI8ttS5cu1dtKlixpPHv2zHL706dPjbfeestInTq1cefOHcvtt2/fNlKmTGmUKVPGcluzZs10O3tCu69WrVpGsmTJjKtXrwa7vW7dusYbb7xhvHz5Uq83bNjQiB07tjF48GDLNk+ePNHLNm3aGMmTJ3f4nl+8eGFUqVLFiB49urF48eIwv3apUqWMLFmyGI8ePbJsc+bMGSNx4sS637BPQ+Pq62TPnt0oV65csG0uXLig5e7Xr5/lNlc/E3v75e2339Y/W3h+vBdznzr73Ow9t6vlMj/LH374Idjjy5Yta+TJk8fu6xEReQqOazlz5gx1mwULFugxdPTo0Xp94MCBduuHbdu26e0LFy603DZ58mS97b333jNevXpluX3dunV6+8iRIy234XG4bdeuXZbbZs2apbft3bs3RLnsbW9r6NChuk2HDh0st/3yyy96G84PrP32229aJx08eFCv9+/f34gWLZrx77//WrZ5/vy51jN4/Lx580Ldb59++qmRJk0ay/Xhw4cbmTJlMgoUKGD07NnTcnuOHDmMDz74INTnev36tdZD7777bqjv39zfR44cMW7evKl1Dz7jHTt2hPr8RO7GnkfyecgeZw4RgevXr8vWrVt12GGSJEkst5vDUtDCGpFManh+tLQ2a9ZMewutYUgnUmujh8t6aEzbtm0t162HAoXm888/l9WrV+vQntq1a4fptZGIAD1lWJbCejhplixZLMNs3PUe8RobN24M1kuHYaxmCzCgPJH5mYRXWMuFll58Ltbeeecd7b3E50xE5C0wMgKjMEJj3o9twwvHyujR/3c6WbFiRR25EZkjMjDi5ZtvvpGaNWtahr8ChuiiVxO325YRw0T//PNPy+PR81m0aFHLNhih4urQVRz3MSR2//79eh11NepW/K1atUpvO3funNaL2NYaptegFxSv1aBBA/2zfi5n0NtasmRJuXPnjvaAYkQTUVTinEfyebbZ2My035ifZytnzpx6iTmAGTJkCNfrYd4bAiMETBh6gv8DLs0hNqgczHkhyZIl02AkLDAXBZXgl19+KR07dgzza2NIDtjLVGediMgd7xGBV9++fTVgxFxObIO5NqgwzSFIkf2ZhFdYy4W5tLZzOzGHFe8Zw3yR0ImIyBukSJFCbt68Geo25v3msMrwcFTPoGEuMuB50bCJuekYFophqdaBFRoCzaG6ODabfxjCi3rLPK7bS8jjSv0IqN8QjCJQRF2BhkYsEYW5mwhmMcwXASVYN9hi+gWSGJUqVUrq1q2rjbMoP4bdWg93DQ2GqKIeQmCKz5goqjF4JJ+HeY7WzF5IM4CyhqQvEDt27HC/njl/DvMBbVsUAT1TaBV0VD5nkKwAcxhQsYwcOTJcr232ltlrTXalhTks7zFt2rQ6NxPJEpDQ5p9//tGeSSQzctdngvvMOS22iW4iIqzlsu7hNuEEAjC/hojIW6BnbebMmRowOWrYMgM8BGLWxzuMsnD1WOuonolIPevIyZMnpU6dOpoEBwlxMJfftu5Knz693cAQcwwx+sY8loe3fgS8PnIKIHjEXE/sL+QFQFCHMiDHAYJHvJ51QIq6HfPycZ910IskO65CEiH0uiLJDnpSbfcBUWRj8Eh+BxPPUTFgOAcS6VjDbbgPCVIALZGOTvod3YeMbqggkAzA3anEkZkNQ1kQtGHyvfVQoLC8Nlo/8T6Rnh0LQltDEhtnwvoeMZQViWiQ2AfBL5ICoIIPz2diD7LUmUOBrOH92QrtM7UV0XIREXkrTJdA8IhsnMj4aQtDJXE/sptimCmYmbQxpB+jZkI71lrfV6ZMGct1BFIYgulsOKWZjdXV4zV689BQiXoJycvQcGkLPXpIpIbGV+tsr7aQXAYJ4PDa1vWsK/WjCT2KY8aM0WypeD7Ue4D9iSzkSJSDJHXWkAgI9bt14IjGS0wzcRUSDOE5MDQXZUDm77CObiKKCM55JL+TIEECDWYw38I6QydaAjHPAQsiYxuzosS8ASyabMvRfQiqkIEUPW22WUaR6dV6/kVYYAgKKjxUiGhRtR0eGZbXRm8nhvUgGxxSqJsQkGI+ozNhfY8YhoNyIx35woULdR6qdatzWD4Te9577z3Nrov07tbLb9gbkhXaZ2orouUiIvJWCOi6du0qP//8sx6bzekHZnCIufQIYqyXLkJWUsyDRMZpE4KxlStXOnwdZDw3pwAAls9CYPrFF1+EWj4zUMXcQGcwpxz1I14Hy16h4c8e9OyhrsCSH9bz0BEkov44dOiQXsd0EIyQsR7dg7nr1pnNnUHg9vTpU13uA1ldTfg/GlER7NrmGECQiSGuZtZzlAtZXhMlSiRhgQAVWVrxHvCZYX8TRRm3p+AhiuJsq5s2bQqx/ePHj40mTZoYMWPGNEqXLq2ZR/H/5s2bB8vMee7cOc0omi1bNuP999836tWrZ8lOGtp9yI42ZMgQI2HChJrNrnr16kaxYsWMdOnSGZ07dw6WoTNz5swuZf40M8cVLVpUX8v678MPP7Rs5+prP3z40KhZs6YRK1YszcqG50Um0gEDBriUbdXV1zH16NFDnxd/+/btC/dnYi8jKjL5Icsrtq9YsaJRpEgRo1OnTnazrYb2udl7blfL5eizRBY/lAGZbImIvA0ykKZKlcrImDGj1gnIjo16AZeHDh0KsT0yrwYFBRmFCxc2ypcvb1SuXNlYvny5w2yrf//9t9YvlSpVMvLmzatZqUeNGhXsOe1lD0UdgzIkSZJEy4Vj9ebNm+1uj4ysuJ4+ffoQ9SP+9u/fb3nelStXakZwvOd33nlH30PatGk16ykygZuQWTZGjBiaIRX1CjK3m+cVzrKtmnUH3iu2R4ZZ0+7du/U2ZHdFVlRrJ0+e1AysqIeqVaum9VTfvn2N9u3bG/Hjx3c526rpxIkTmlUdWc9ZB1FUiYZ/oi5UJQofrOOHHi0kYDET0aClDS14yOxmDhexhSEi+/bt03lpmJdgLykA5rahdw69VWgFxHBLc85faPeZw02QdRRJZDDPAlnerLPbYegj7kPPmS0MmcHcxFq1alkm+mNtRnswrMZ2+Iuz1zYdPHhQkwNg3gW2QesqbsOi0NZDZxxx9XXQo4kEO+hxNN9TeD4T2/1iDY/D/B2UAe8HLeL4w76xHnrk6HML7bmdlcvRZ2lmnsVwKnu9xURE3lCH4viG4xyGlPbu3VsTvIwbN87u9jhOYt1BJATD8FMcS7EGItYlNI+N6DHEME8Mo0SdgPmTGJmCOhqPs4a1JHE/5gVaD7HE8Xn37t1aLqwtiWGnSFBmu/29e/csCWjswZqNmIdo/byo51BGPB5rOdpLLoNeShy/kWQHr40RKxhuivmiqO+cQZlQNvTimnPicVqNuYjYJ6gXbGFYL94z9inqGrwO9jV6VdG76mh/We9v61wKeA9ItoPPBZ8PUWRj8EhEREQUQJB0pXv37tKlSxcZMWKEp4tDRD6EwSMRERFRgEFSF2Sxxhy9sM65I6LAxeCRiIiIiIiInGK2VSIiIiIiInKKwSMRERERERE5xeCRiIiIiIiInGLwSERERERERE7FkACE9X+whg7WycGabkREFBiwBhsyTGJNNOu1Qf0F6zciosBkRFH9FpDBIwLHjBkzeroYRETkIRcuXNDFyP0N6zciosB2IZLrt4AMHtHjaO5crm0Uuv/++0/Kly8v//zzjxQuXDhKPh8ioshy//59bTw06wF/4+/1G+skIuIxxrP1W0AGj+ZQVVSs/li5ulOCBAksl9xXROQv/HXKgr/Xb6yTiIjHGM/Wb/434YPcKkWKFNKuXTu9JCIi8iTWSUTEY4xnRTMwuzLAoFs3ceLEcu/ePb9smSWiqPHq1St58eIFd7cXihkzpgQFBQXc8d/f3x8REXn2+B+Qw1bJdY8fP5ajR49K7ty5JV68eNx1RP/fw4cP5eLFi5rdjLxz2A4SBpjDHMk/sE4iIh5jPMsrgsfjx4/rX+nSpSV58uQuPebAgQNy7do1yZs3r6akpciBwLFo0aKye/duKVKkCHcz0f/vcUTgiAaVlClT+u38OV+FgP7GjRv6GeXIkcNuD2RUefTokWzatEnSpEnjctIxlH3//v2SLFkyfQy/X//DOomIIhOPMV4ePKJCHTBggBw7dkwr+fXr10uFChWctvbXqVNH9u3bpycFyLzWq1cv6d27d5SVm4gCG4aqIkBB4Bg3blxPF4fswGdz9uxZ/aw8ETzevXtX+vTpIwsWLJCnT59K9erVZcaMGU4f98svv0i3bt2kQIECcu7cOcmUKZOsWLHC5YZVIiIiv02Yg/WounfvLlu3bnX5MaiMz5w5owHntm3bZMmSJXrbxo0bI7WsRES22CPkvTz92WDOSc6cOeXIkSNSokQJl0fUfPHFFzJ16lTZsWOHnDx5UodpfvXVV5FeXiIiIq8PHhs2bChVqlRxuZJHS//06dOlVatWllZYPB7DKadNmxbJpSUi8k23b9+WDRs2RDgYWrNmTbheG6NKAk3mzJmlY8eOmrzAVeiZxOMaNGig1zFfs3379jJv3jx58uRJJJaWiIjID5fquHTpkty6dUsKFSoU7HbMCcEwVkeePXumGYis/8g10aNH18VGcUlE3gtz6xYtWmT3Pswpj+jQfoz46Nq1a5gfh9fu2bNnhF47UKAes1e/Ydgr9qM9gVa/sU4iIh5jPMsrEuaEZQ4JIImANfRCmvfZM3ToUJ1bSWGHExd/Pxkh8gdXrlyRxo0b67xw8k2ox9KnTx/sNnOUjaM6LtDqN9ZJRMRjjGf5VPAYK1YsvcQcEGu4bt5nD1q9O3fubLmOYChjxoyRWFIiIs9Cj9T27dslSZIkdu9HNk+M5kBiFixpYR5Lly9frlMJELSULFnSpYRAGJqKXjNkZt61a5cOvcyePXuwbR48eCA7d+7U7Nh58uRx6fVQPiRFw5BNHLNxf2jl93Wox+zVb+Z99rB+IyKiqORTwSOyziFrHjKzWrtw4YJkzZrV4eNix46tfxR2hw8f1vk3mHODZVGIyPthmGOZMmXk+fPnOufOel45lhmpV6+eJixLnTq17NmzR3uvmjZtqkNfZ8+erdvdvHlTzp8/r4lbkLk0NBhS+emnn0rSpEn1OREkfv/999KyZUu9H8sqvfPOOxog4vl++uknadSoUaivt3btWvnoo4+kVKlSevzG+0HwGFr5fR3qsdOnT4eo38z77Am0+o11EhHxGONZXh88onUZw3XKlSsnceLEkfLly8uff/4pzZo10/txH04yhg0b5umi+u1JKCprXBKR86Gj+LOGgAon/uZvyZa5fioySCOYspYlS5YQw/RdgWydCMDMXr0ePXrI5s2b9b4//vhD7ty5I19//bVer1y5svTv31+DLzxmypQp2nuI3sRJkybJ3LlzNWmLM1evXtWs12jkw3G7YsWK0rx5c70Pz4WAEs+/dOlSGT16tAaPob3e6tWrNfMosmlbC638vgb1Fz4XBMboIa5WrZq+D3yH0qZNq9vMnz9f50FinUhinUREkYvnvV4ePKIHEUOSkAQHsPQG5usgvTn+AC3UGHp18OBBvY4WZgSSbdu2ldKlS8vYsWP1BMts4SYi8pTff/89xPwzBEnIoonjHYZ12ssiDQi0cKyzhuzSmMcYVocOHZKqVataehzxfzN43Lt3rwZfZo8foHcPsKwEgj6soYvADj2BCApdkStXLg0coWDBgtobht5B8z6z9xKBtHnMD+31sN8++OAD3QeVKlWSLl266Hahld/bIHh//fq13LhxQ4cRL1u2TIflIuA1F6OuVauW1n14Dx9++KH8/PPPuibkl19+qUt3zJw5U5+HiIhIAj14REs7gj+oUaOGbNmyRf8++eQTS/CIFlfrITlYLwst2Hjc4sWL9aSoU6dOXKibiDyuTZs2Urt27RA9j4B5ebt373b4WPTA2et5DA8EYlhg3mT9/1SpUuk8QfTg2UJA1qJFCxkyZIhex7JIZnDrDHrLXrx4ITFjxtT5jehVwzBVzEu0zdZsPmdor4cynjhxQntrFy5cqHUEhseGVn5vbEzAMFuzFxH1VooUKSzBI74beF/mdwTTMtDjigASgSZ6nRH0W8/1JCIiCtjgERWoWYk6gnWybKFV+9dff43EkhERhR2CBDNQsIVh9+YQVXvQO+cu6K3EayGIxHBIjOAwe/4wNxGNchgaiuGSCPYwd7Bs2bIarI4YMULLgkBtwYIFLg1ZBSS1QW8hes3QW4Zew3jx4oX6mNBe799//5WzZ89qzx2GtGLJIGfl9zZo4AwN3jeCRGvx48fXYcZERETeyOvnPJJnvfHGG3oChEsi8l5YUP7999+3DA1ds2aNjB8/XnsBJ06cKH///bfeh7lzyIz622+/6fxDJNVBgxyCL8wlRzbq9evXa+8etjEhCK1SpYrD18+XL5/O18NakxUqVNBhpoDeRww7NSGBj/k8ob0eEufg9hgxYujxZ8mSJU7LT/6PdRIR8RjjWdEMV8ck+RGcrOAE5t69e5IoUSJPF4eIfHBC/ZkzZzRIQ49ioMNcTUwfsJ2z6Y2fkb8f//39/RERkWeP/8EnohDZQPIKJClyNWkGEQUe295FosjCOomIIhOPMc4xeKRQIVviN998Y8maSERkC1lQzaQ3RJGJdRIR8RjjWQweiYiIiIiIyCkGj0REREREboBlipDIi8hfMXgkIgqnAMw35jP42RCRJ2AJoT///JM7n/wWl+qgUCE9f/369fWSiP4P1haMFi2a3LhxQ9dPxP/JuwJHfDb4XPBZkf9gnRT5Hj9+7HSNViJ/xWOMcwweyemaWvPmzeNeIrISFBQkGTJkkIsXL+pC9uR9EDjiM8JnRf6DdVLkQdKrn376SdP9I3j84IMPpH///pIuXbpIfFUi78JjjHMMHilUGLd//fp1SZUqlcSKFYt7i+j/S5AggWYZffHiBfeJF0KPIwNH/8M6KXLMnDlTvvvuO1m+fLm89dZbGkDOnz9fr7dq1UpevXold+7c0W0RWNr2TOI4+ODBA0mWLJn2/GOd1bhx4wbrzcR1R6M0sL2ra+bi+fF6OCfB875+/VqPx9ZlAKx1hzXvQnuNhw8fSvTo0YO9H8xZxHV7z49t7R1XcB9ux/2uevLkSbB9BK6+niuPf/nypb5n/B+ePXsmsWPH1v/fvHlT9xPux2ebMGHCYM9z+/ZtfVwgnvfxGOMCIwDdu3cPE5X0kkK3e/du3Ve4JCLydf5+/Pf398c6KXJ07drVqFy5ssP7T5w4YSRPnlz/4saNa+TOndvYuHGj5f7Vq1cbKVKkMNq1a2ekTJnSiB49ulGrVi3jv//+M4oXL66PSZo0qbF48eJgz7t06VJ9rnjx4hmJEyc2evToYbx69cphOUaMGGEkSZLEiB8/vlGmTBmjQYMGRrNmzSxlwGt//fXXelmgQAHj8ePHRvv27Y1EiRJpGXDbpk2bLM+Hx3bv3j3Ya+TKlcuYN2+e/r9Lly5G9erVjZo1a2r5YseOrbeZHj16ZDRp0sSIEyeOvr/PPvvMyJEjhzFr1iy75X/58qUxcOBAI1WqVFqeTJkyGQsWLLDc7+z1XH08/pIlS2Z06tTJOHPmjPH222/rc+Fx3377rR4jcPv8+fON9OnT6/OaTp8+bQQFBRmnTp0yApEvH2PuRdHxnwlziIiIiAJYlSpVZNOmTTJ06FDZs2eP9lJZy549u/ZW4e/Ro0fy9ddfS8OGDbX3y4T7UqdOLVeuXJELFy7Itm3bpHr16jJhwgTtERs8eLB89tlnlu337t0rH330kfzwww/aA3jw4EFZunSp/P7773bLuHbtWh1Gu2LFCt2+T58+snDhwmDbYK4zetwuXbok+/fvlwEDBujj8H/0SjZu3Fhq1aql27lqzZo18sUXX2iP5K5du2TMmDF6CSjPgQMH5PTp03Lt2jVJmjSpnDhxwuFz4b2il3fz5s26T6ZNmybNmzeX48ePu/R6rjz+77//ls8//1w/jx9//FHfM4Zi4vlOnjwp27dvt2xbp04d7cldtmyZ5Tbs/4oVK+pjiOxh8EhEREQUwKpWrSqrVq2SQ4cOSYMGDTQIQmBx6tSpYNsh0MBwUARgGE6JoMyEIZS9e/fW2zFPsnTp0ppwr2DBgnp/vXr1NMBCUAMTJ06UGjVqSIkSJeTWrVs6pLRly5ayYMECu2WcPXu2BpulSpWylBmPtxYjRgydu2kmykKg1bNnT8mcObOWq1u3bpIoUSL566+/XN431apV0+AaChQooH/YTzBr1izp1auXpE2bVl/z22+/tQwTtWfcuHHy5Zdf6v7FfsiXL5+UK1dOlixZ4tLrufL4SpUq6eeDIcIIords2aKNAti/GJ6KIN56f7Vu3VrGjh1rGbI5efJkvY3IEc55JCIiIgpw5cuX1z+4evWqBhC1a9fWwAWBCuY+olcLAQfmzqEn6/Lly5bHIyiznvOH+XIIcqyvA+bZAXon0cuGnkFr2bJls1s+lKlYsWLBbkNSLPRCmpInT26Z12c+JlOmTJbrCKgyZsyovaOusn4PgOc33wOCYTyf9XtE76sjeM/du3fXgNOaGWA7ez1XHm+d4AjvH58JglvrfWYNvcEILs+cOSM7duzQOZN169Z1+B6IGDwSERERBTD0KFons0mTJo18+umn2nOIIazozUMAg2AESWiwPQISJFsJLwSJ6AmbMWNGiLLYg97DY8eOBbvt6NGjIYIh28ccOXJEKlSooNeRhAbDSrNkyaLX48ePr8NwTXg/CAhdZT5/yZIl9TqGxqK3L7T3/M033+hQ0vCsSxvWxyNwRjCI4ao5c+a07DNrCDbRU4leTQw1btasWUAmyiHXMXikUBUuXFgrDK6VRkREnsY6KXJg/iDm0GF5DgRWWIJo+PDhOvQUPV8ITtDDh+ARGdhHjRoVpiDLng4dOsibb76p8/jQ04VzDWR3RabXgQMHhtgeQ1pRnqlTp2oPKYZqrlu3Tpo2berwNdq2bSuDBg3SOZuYw4f3hMAIPapQvHhx7cVDMIZey++//16H5boKvXaY94jnRjDdt2/fEPNFrfXo0UN7DjG0tWjRorrcE4aJYr9juKozYX081iHG8GPsByzDgrJ99dVXIbZr166dPgf2vTmENVDxGOMcg0cKFYY7WA8BISIi8hTWSZEDvVnjx4/XSwxfxDIO6K3DdcBcRgQYmIuH3romTZrI22+/bTk/QECG4MsahrFaL4GBzw7bmEtPIOBCTxeS2iBJC+bkYQ4j5ijag0Bzzpw5OmcPASECSAR9oZWhU6dO+np4Tiw/UaRIEdmwYYO+B8DjkfCmUaNG2qOKABW9iOZzYjvbYBCLyJtLZHTp0kWD3jZt2mj5EchiaKmj8ybcj/2CQA7zSdFritd89913XXq98DwewTYCRvQuIsDF59iiRYtgZcQ8SfQ2FypUSHLnzi2BjMcY56Ih5aoEGLSs4CCB1iX8CMkxZPBCyxqGM5hDHoiIfJW/H//9/f2xTiJrCGDRs4bsr+Sa1atXay8jjhHmHFUkykEgOnLkyBBDYgONLx9j7kfR8Z89jxQqDFP5559/gk1IJyIi8gTWSYENvWzoNUMv2fTp03VZkT/++MPTxfJq8+bN0+HG6NVFryiytSJANANHZLpFrzMSIX344YcS6HiMcY7BIxERERF5PQzbxPxMrKuI4ZUbN260JL8h+zBcFcOP33vvPR1ai15HzM00vfXWWzq0F2tGMlEOuYLBIxERERF5PcxzRJIcch0CRgxHxZ89thlsiZz534I8RERERERERA4weKRQYY0gjIW3XmSXiIjIE1gnERGPMZ7FbKt+mI2OiIgCMxupv78/ImQGnTRpkly5ckWX+SCiqD3+s+eRQnXz5k2ZMGGCXhIREXkS66TA9fLlS5k8ebLkypVL1yq8e/eup4tEfojHGOcYPFKozp8/L61bt9ZLIiIiT2KdFHhevXolM2fOlLx58+oyHUWLFpX9+/fL6NGj3f5a27dvl8KFC8uuXbvc/tzkG3iMcY7BIxERERF5ldevX8uCBQukYMGC0qhRI+1xxLqO8+fPl/z580fKa+I1kJ20XLlyXD+SyAEGj0RERETkFQzDkGXLlmkPY/369SV9+vSybds2Wbp0qbz55puR+tpJkyaVDRs2SMOGDaVx48bSvXt37fkkov9h8EhEREREHg8aV61aJaVLl9aF7ZHw459//tHbSpUqFWXlQM8j5lZiXcQRI0ZI7dq1NQEJEf0fBo8UqgQJEuiivLgkIiLyJNZJ/glBIs413n33Xb2+evVq7QHE8FFPiBYtmnz11VeyYsUK2bp1q5QsWVKOHz/ukbJQ1OIxxjku1cFU5kREAcPfl7Lw9/dH/gUJavr06SNr1qzRIakDBw6U6tWra/DmLU6cOKG9j1gaZPbs2VKtWjVPF4nILi7VQV4zYf3Zs2d6SURE5Emsk/wDEt/UrFlTh6giKENinN27d0uNGjW8KnCEHDlyaJBbpkwZLd8PP/ygQ2zJP/EY4xyHrVKo/vvvPx3/j0siIiJPYp3k2w4ePCj16tXTZDjo0cMSHPv27ZMPPvjA64JGa+jNX7x4sXz99dfStWtXadasmTx9+tTTxaJIwGOMcwweiYiIiCjSHDt2TD7++GNddmPv3r2akObQoUN6W1BQkE/seZRz6NChGvDOmzdP52hevnzZ08UiinIMHomIiIjI7U6fPi3NmzeXvHnzyubNm+W3336To0eP6m0xYsTwyT2OgBfvBYFjsWLFZMeOHZ4uElHgBY9HjhzRrFrXrl1zafuXL1/K/v37NTvXmTNnIr18RERE4XH79m2tqw4cOODyYy5evCgbN27U4VOYc07kay5cuCBt2rSRXLlyycqVK+XHH3/UYaq4LVasWOLrMOx2165dkjVrVu2BnDZtmqeLRBQYweOjR480LfNbb70l3bp1kyxZssiwYcNCfQx+rNmzZ5c6depI3759pXDhwvocDx8+jLJyExEROfP7779LhgwZpHPnzlKpUiWt6xBMOvL8+XNp0KCB5MmTR3r37q3/x8npunXruLPJJyD5zRdffKHnaUiCg2Ge6H3Ebcif4E/SpEmjv81GjRrpHMguXbpo5waR3zM8qHPnzkaWLFmMGzdu6PWVK1cifZWxadMmh48pUaKEUbt2bePVq1d6/erVq0bSpEmNoUOHuvy69+7d09fBJYXu2bNnxoULF/SSiMjXRdXx/+DBg0b06NGNWbNm6fX79+8befPmNZo3b+7wMePGjTPixIljnDlzxnJb06ZNjWzZsrn8uv5ev7FO8k7Xr183unbtasSNG9dIkiSJMWjQIP3OB4LXr18bo0ePNoKCgox3333XuH37tqeLRAF6jLkXRcf/cPU8jhw5Um7cuBHhwBXd/C1btpQUKVLodfQgYp2fqVOnOnwM1q7Knz+/RI/+f0VPnTq1tv7gdnI/DC9By7k/DDMhIooq06dPl0yZMslHH32k1xMmTCjt27fXdeIcZWlEPYasjpkzZ7bcVqhQIdZvVlgneZc7d+5oL/kbb7whY8eO1UykmE7Uq1cv/c4HAmSJRc8qhufu3LlTSpYsqfM6yTfxGONcuILHUaNGSfr06TW18rJly+TVq1dhfg7M6bh586YGi9ZwHWmbHRk+fLgGnaNHj5ZFixZJu3bt9IfbsWNHh4/BnBEsnGn9R67BcBMMncIlERG5BvMVMa3Ctn5D4Hj8+HG7j/n000+1sQ7JRBYuXCi//vqr1nWYL+ZIoNVvrJO8A75nAwcO1GHV6FDAuRiCxm+//VaSJEkigeidd97RqVUxY8bUAHL58uWeLhKFA48xkRQ8nj17VpYuXao/kPr162vras+ePXUytKvu3r2rl8mSJQt2e/LkybUlK7RJyqiQcbDCH9IlowzogXQEY+7Rmmv+ZcyY0eVyBjp8TvPnz7d8XkRE5Nqx0179Bo7qOGzfsGFDWbFihS5EjjoOc8dKlSrl8HUCrX5jneRZyFXx3XffadA4ePBgbejAyTYa9s1RZIEsW7Zssm3bNk2iU7NmTd0vhoGRhOQreIyJpOARQ0YxxHTOnDmaqrhHjx7aXZ8zZ04pW7asTJkyxWmGOHMY5JMnT4Ld/vjxY4dDJPEDrFatmg6FwMEK2eiwTtCkSZOkf//+Dl8LgS2GA5l/yAJGREQUWVCP2avfzPvsQQ8jgkEMfcNSAGiQRfIc9Gg4qlNZv1FUQI85Rp1heCqGqX744Ydy8uRJvQ1Th+h/EiVKpCPjvvnmGz0/bty4cYhjAVFAZ1tFSymG4uAPFeKlS5c04xQyp65evdrh49A6igVXbQM5DGfFY+3BcyNYbNq0qWVR2VSpUkmNGjU0eHUkduzY+mO2/iMiIoosqMfs1W/mffasWrVKs7Ka92NKBnp2zp0753AOFes3ikzIAIy1GdEDjvmMON86duyY3oYh1uS4k2XQoEHayYIh6OhYMX//RAEbPKLHEctqYA0ftIqiRfWvv/6SU6dOaZCHOYhIhuNI3LhxpVy5cto6Y0Kv4Jo1a7R30YRgccuWLfr/lClTatBoO/8Or8mWLyIi8haox7Zv3y5Xr1613IalCwoWLChp06a11Hlo+DQTvqEes127GPWbeR9RVMGSExjVhXM8JHrCMMzDhw/rbRiySq5BDy3OYa9fvy7FihXTIa1EPi88KVpr1KihKYnz5ctnjBw50rh582aIbbCUhrOn37ZtmxE7dmyjQ4cOxh9//GGUKVPGyJ07t/Ho0SPLNi1bttTXMXXq1EnTQH///ffGggULjDZt2mg69HXr1rlcfn9PZe5OV65cMYYMGaKXRES+LqqO/y9evNClpYoUKWJMnz7d+Prrr7XeXLFiRbA6EGXBJezevduIFSuW0ahRI2P+/Pma/j9VqlRG48aNXX5df6/fWCdFrpcvXxozZswwcuTIod+j+vXr67IzFDFYVg7nuPh9T5o0ibvTi/nyMeZeFB3/o+GfsAacLVq0kFatWumCx84S6zganmPau3evDn+4du2aFChQQBdTtk4ygPH0GK6DFND/P9jVBC5///233Lp1S1Oao4cTjw1LljAkFkBrL4ewEhEFjqg8/j98+FDrsH///VfrNdt6E3Vbp06ddJvcuXPrbZjnOG7cOO1xRDkrVqyoi5CbUzWcYf1G4fH69Wv5888/pV+/ftrDWKtWLc2capsxmCI2BLhDhw4yfvx4+fLLL2XEiBESI0YM7lJym6g6/ocreERGLSyzEdb7vAUr17BlnUJiIgwxDtT020TkP/z9+O/v7491knvhFBBLrvXp00eXSatataoGjVhqgiJnf2MJHgSPaBjCnEjbrMzkWb58jLkfRcf/cM15RI+fPS9evPD7NaYCDeaX1qlTh+s8EhGRx7FOcl8QgwRNWAamdu3aepKME2aM6mLgGHmQBAtzSJFQEiPvSpQoobk9yHvwGONcmPrLMVzU3v/NIQ9IDoCMXERERETkff755x9dbgPLwSB4RKJCZPlFYENRA72Ou3bt0sAdn8HMmTN1qDCR3wWPmK9h7/8QM2ZMnd/4888/u690RERERBRhyPSJ4alr166VIkWKaIb89957j0GjhyBr7datW3X5OYzwwtIeWLeVQTz5VfCIccCA3kUsDktERERE3mv37t3St29fWb58ueTPn18T49StW5dBihdImDChLuEzYMAA6dWrl+zfv1+XQ4kXL56ni0bk3jmPDBwDR5w4cSRv3rx6SURE5Emsk1x34MAB+eCDD3R9QWTvnTVrlibFef/99xk4epHo0aNr8IjpYEuXLpUyZcrI+fPnPV2sgMVjjBuzrc6ePVsvP/roI8v/HcE23szfs9EREVFgHv/9/f2Rc8eOHZP+/ftrJk9MJ8LyG1juhctCeD8E9xjC+uTJE+2RRCBJ5LNLdWAJDsAyHOb/HeFSHURE5I38Pbjy9/dHoWeJxDIb06dPl3Tp0ulQ1ebNm2tOCvIdN27ckAYNGuh8SCzrYZtjhMhnlupAQGgGheb/Hf2R//jvv//0C4hLIiIiT2KdFBKGOH722WeSK1cuXWpj1KhRcuLECWndujUDRx+UMmVKXcoDQSM+w44dO+pSeBQ1eIxxc8IcCjxYguXBgwd6SURE5Emsk/7nypUrMmTIEBk3bpw28g4bNkzatm3LZCt+AL3F6HUsWLCgBo9YC3LevHmSPHlyTxfN7/EY48bg0dk8R1+a80hERETkq8Mahw8fLr/88osm98CcRgQYyNxJ/uXzzz+XPHnySP369aV48eKyePFiKVCggKeLRQHO5eCxQ4cOLj8pg0ciIiIi97lz546MGDFCRo8erRk6v/76a/nqq68kSZIk3M1+rHz58rJr1y5NpFO6dGmZMWOGLrVC5PXBI+cyEhEREUV9EgzMYxw5cqTOfUMvY7du3TiEMYAgay4S6DRr1kyXWkFipN69e3PJFfIIznmkUOXOnVsXGMYlERGRJwVSnfTo0SMZM2aMfPfdd/p/zGfs0aOHpE6d2tNFIw+IHz++zJ07VwYPHqyZdPfv3y9TpkzR28l9AukYE15c55GpzImIAoa/L2Xh7+8vEDx9+lTGjh0rQ4cO1aGqyLrZq1cvSZ8+vaeLRl5i4cKF0qRJE8mePbvOg8ycObOni0RegOs8+sHO9ZcU4JiY3717d8mUKZOni0NEFCH+fvz39/fnz3XS8+fPZeLEiTJo0CC5du2artGIoYkYskhk68CBAzoPEhnxFyxYIOXKleNOCvBjzH2u80jeAHNdkS6ac16JiMjT/LFOevnypUyaNEly5swp7du3l0qVKsmRI0dkwoQJDBzJIWRd3blzp15WrlxZe6sp4vzxGONu0d3+jEREREQUqlevXmnmTCzF0LJlSylRooQcPHhQpk+fLjly5ODeI6dSpEghf//9ty7pgTmx7dq106RKRF4ZPG7btk0aNmyoC5jiD8tz7Nixw72lIyIiIvKzRcjnz5+vPUaYt5Y3b17Zu3evJkPB/4nCImbMmPLzzz/LuHHjtLe6SpUquhYokVcFj8juVLZsWXn27JkGkPjDBO+3335bpk2b5v5SEhEREfkwwzBkyZIlUqRIEWnQoIHOp0KjOxKeFC5c2NPFIx/XunVrWbdunRw+fFiKFy8u+/bt83SRyE+Fa6kOpAhGC0eLFi2C3Y6J3rivadOm7iofeViqVKl0EWJcEhEReZIv1kkIGletWiV9+vTRxd6x6PumTZukTJkyni4a+Rl8p/7991+pW7euvPXWW9qhU69ePU8Xy6f44jHGa5fqsIYMbhcvXgyRyQ1ZfjJkyKCX3szfs9EREVFgHv/9/f35mg0bNmjG1C1btkjp0qVl4MCBmhAnWrRoni4a+bHHjx/Lp59+qkOh+/Xrpx070aMzzYm/u+9t2VatYXjFP//8Y/cg+eabb7qjXOQlHj58qPNbcUlERORJvlInbd26VTNgVqxYUaf1LF++XANI3MbAkSJbvHjxZPbs2TJ48GD59ttvpX79+l7/m/EWvnKM8SSXg8dly5ZZ/qpWrSoff/yxdO7cWWbNmiUzZ87U/3/yySd6H/mP48eP69AHXBIREXmSt9dJu3fvlurVq2sOCCQtwWLuGKr63nvvMWikKIVGim+++Ubn1K5evVp/N2fOnOGn4OPHGJ+a84hWC1tYB8UWhmT06tUr4iUjIiIi8pEF2zE0cNGiRZIrVy7t9UFSHA4VJE+rVauWbN++XerUqaOJdObNm6c94kSR3vOIYReu/hERERH5u6NHj+pSZYUKFZL9+/fL1KlTda1GZKFn4EjeIl++fLJz506dWoalPH755RdN5EQUHpw9S0RERBQGp06dkmbNmulJOeY3/v777xpIItt8jBjhSmRPFKmSJUsmK1askA4dOuhfmzZt5Pnz59zrFGYRPsLdvXtXXr58Gey2FClSRPRpyUugEsTnycqQiIgCvU46f/68DBo0SCZPniwpU6aU0aNH6/p6sWPH9kh5iMICv5tRo0ZJwYIF5fPPP5cjR47IggULuCyFzT7ieW8kLNWBFLBIkINx0w8ePAhxv7d3hTOVORFRYPL347+/vz9PuXLligwZMkTXuMZ+7dmzp7Rt21bixo3r6aIRhQt6zD/44ANt+MBcXa6W4Pvue/NSHd26dZOTJ0/KkiVL9DpS2qL1LXny5DJs2DB3l5GIiIgoyiFjateuXeWNN96QP/74Q/r3768ZK9GAzsCRfBkyiiITMHrQkR0Ya0ISRVrw+Ndff8mECROkQoUKer1EiRLyxRdfyIwZM/jl8zOHDh2S7Nmz6yUREVEg1Em3b9/WzPFZs2bV3sbu3btr0IgexwQJEkTqaxNFlYwZM8qmTZukbt26muSpd+/e8vr164D+AHje61y4Jg1cvnxZD96AblEcZDE+uGzZsgwy/MyzZ880MQAuiYiI/LlOwrCvH3/8UUaOHKn5HNAwjp5HjKwi8kfoQUevOjIGo3EE2YKnT58uCRMmlEDE895IzLaKxUchd+7c8ueff+r///77bybLISIiIp/y6NEjGT58uPY0Dh06VFq2bCmnT5/W/zNwJH+Hc3r0ri9dulTWr18vpUuX1kYaIrcFj0hNbcKwjo4dO0qaNGl0QVx8+YiIiIi83ZMnT7SnEXMa+/Tpo2s24qQZPY+pU6f2dPGIolSNGjVk+/bt2vtWvHhxWbt2LT8Bcs+wVXRpm2rXrq1DVf/991/thSxcuHB4npK8SJYef1n+/+zqSY+WhYiIyN1wcjxx4kQZPHiwXLt2TZo3b67BY+bMmbmzKaDlyZNHdu7cqQ0p7777rjakoJPIHHFI5JaFkjD/0ZwDSf4lZtJ0kqrBAH6+RETkcTjXWLlyZbjrpBcvXsi0adNk4MCBumZj48aNpW/fvqzjiKwkTZpUk2NiNOGXX34p+/fvl19++SUg1jON6DEmEIR7ziOW50BmJiw0ij+0UOzYsSNcz3XixAnZvHmz3Lx50+XHPH/+XHs7Dx8+7PXrSvqy6LHjSdw3inK9MCKicLh7965s2bJFF+MOCwydxPCxx48fc79bQZI+9IaEdQ2zV69eaUb4vHnzSqtWraRkyZI6agqBJE8SiUKKESOG/PDDDzJ16lRNoFOpUiW5evWq3++q8B5jAkm4gscpU6ZoZlUM+0AAib+nT5/qOjE4ELsKlWL16tWlWLFi2iWOlMEjRoxw+rhZs2ZJunTppFmzZvqHJUOuX78enrdCTrx8eFvubv5DF0gmIiLXYUmr9OnTS7t27bR+RL2JYDI0Fy5ckDJlymhwgzWVCxQoYFlTmUTrIqy1aK9OQtKbNWvWBGtQxrID8+bN0/3YpEkTzdnw33//yZw5c3R4HhGFrmnTpvLPP/9oAinMg9y9e3fAHmPo/zPCIWPGjMbEiRND3D5hwgQjc+bMLj9P165djUyZMhnXrl3T63/99ReO+MaWLVscPmbt2rVG9OjRjT/++MNy29atW43Dhw+7/Lr37t3T18ElhZS5+zLLX5pmo3Rf7d69m7uKiHxeVB3/Dx06ZAQFBRkzZszQ63fv3jVy585tfPrppw4f8/TpUyNv3rxGjRo1jMePH1vKO3/+fJdf19/rN9RF9uqk169fGw0bNjSSJk2q/8ff4sWLjUKFCun21apVM3bu3OmxchP5uosXLxrFixc34sSJY8ycOdPwV46OMb7gXhQd/8PV84iW0/r164e4HdlWseajq9AVjuEjqVKl0uvohUTCHfRsOvLtt99KtWrV5JNPPrHchpTCbEEkIiJvgWFeGTJkkEaNGun1xIkTS4cOHWT27Nk6Usce9JAdPXpUxo4dq2uvAYZO1atXL0rL7osmTZqkvYm//vqrrFq1Sntu69SpI8mSJdNpMStWrNBeEyIKH4yiQA8kzv9xDo41ITEcnAJPuIJHBHj4AtnasGGDvPnmmy49x6VLl+TGjRtSpEiRYLfjOoaU2INhspg7giATQSr+f+bMGaevhcdh4V/rPyIiosiCesy2PkT9hqUhjh8/bvcx69at0xwCmJaBoWF4DmwfGtZvovNJMfUFywyMGTNGG5hjxoypywxgn2LIMBFFHBq1MD3t+++/l++++04baHhOHXhczra6bNkyy/+rVq0qH3/8sXz22Wfakof5BUheM27cOG2JcMWdO3f0Eq2C1rAYr3mfLSTUefnypb4W0mtnyZJFW2kxlwEttmYPpi0s8jtgwABX3yoREVGEoB6zXhMZzMXmHdVxmGMTJ04cKVGihNarmMN369YtGT9+vNStW9fuYwK9fkNwXatWLV1GANkhixYtqr2MSHjBpQWI3A+/q65du0r+/Pk1WWapUqVk8eLFkiNHDu7uAOFy8GhvmCqGh9hC+utevXo5fb5YsWLppW2rKpLomPfZQksibNq0SdeaROB57949HbaKL7KjZD0IaDt37my5jlYSJOch56LHSSDx81bQtM1EROQa1GP26jfzPkd1HDKsLliwQD744AO9DctIINHLuXPnQjS2BmL9hroIQ4HNOgn/R2ZanNCiVwS9kAi0sY83btzItaeJIgl6+LHKAnof0eCFYePoXPK3YwxFIHh0NEcjvFC5RY8eXS5evBjsdlxHj6I9KVOmlAQJEmjFYFaimEeCShYZWB3BujSBsDZNZIiZJI2kqNVVsmbN6umiEBH5DNRjZ8+eDXabWd85quNwnDXrNFOLFi20UXbfvn1SsWJFCfT6DfsIS26Y2rZtKwkTJtReXiwtgPOKoKAgiRcvnuTMmdOjZSXyd7ly5dIGL4xGfO+993TFhE6dOvl0r7/tMYYiEDy6G1oIkbZ80aJFmgYYHjx4oHMUBg0aZNkOrYjoXUS3OL6MGIqCFlhrqKBTp04d5e8hEBgvn8vLBze18QDDqYiIyDm0wH/66ae6jJQ5pWLhwoU61Ctt2rSWXsKdO3dqqz0S4+DkCwtxY6iqOcTVDEBZx/0f1EUIwpGMCHVSlSpV9I+IPCNJkiQ6tc0cBYGGLiT98tVzRttjDIUUDSlXJRww9xCBH4I7PAUW3kWPIFr+XLV161ZtScUaWBh6ikrz2rVrmiggfvz4ug2ysaJVA8NUAa+Hbdu0aaOT4NFljkm7KAsmy7sCFTZadxGUchHQkLL0+Mvy/2dXT8rVqZ30M7FNbkRE5Gui6viPOhJ1FRo9cUJ14MABGT58uK7ZiKRvgLoN22zbtk0bSAEBJMqI1nvMeUSvI5LUYSirN70/T9VNZp2UptkoiZ0mu5wd5lq9T0SRDz12OG/HMQuNZWZDmS/Zs2ePzp32xfPe+1F0/A9XtlXML8AQEczDwBhnJKsxF9/Ffa566623NGMq3iSW7UAlipTaZuAICEpxuwlLciBgfPjwoSboQeIBVMCuBo5ERESRDQ2pGEmDQBEnVBcuXNDMn2bgCKjkK1eurJcmNIRiaQ48BktOdO/eXetZIiJv17hxY81LguNdsWLFZNeuXZ4uEnlLz2PNmjW1YpwwYYKkSJHCkgkVrQ1obbXOzOqN/Lll1h3Y80hE/srfj//+/P7Y80jkG5A5GnO39+7dKxMnTrSsd+sL2PMYSXMe169fLydOnLAEjoD/I/sqU/USEREREQUmDFdFrICEVuiNxDxILCuEZFbk+8IVPOLDf/HiRYjbnz9/HqY5j0RERERE5F+QbGbSpElSsGBBXU4PuUtmzpypCXbIt0UP79ouSB9+5swZy22nT5/WzHLIhkr+AwkJMndf5nOThomIyH/rJFwSkXdDwrCvvvpKVqxYYUkMdvz4cfFmON/FjD6e97o5ePzpp580le0bb7yhKcjxly1bNu15xH1ERERERERYugjLEmEdVixNtHLlSu6UQAse06RJo1lSkVEJacSxLiP+j9twH/mPF7cuypXpXeTYsWOeLgoREQU4s07CJRH5DuREweoIZcqU0RUSRowYoT183gbnu1jlgee9joVrgmKBAgV0zSp8AfBH/uv1i6fy/PIxXW+MiIjIG+okXBKRb0EG6MWLF0ufPn2kW7dumkhn/PjxOj/SW+B8F0Euz3vd3POI9VuQBpyIiIiIiMjVpJtDhgzR5Dnz58+XcuXKyaVLl7jz/D14xNot48aNc39piIiIiIjIr3388ceyefNmXROyePHismPHDk8XiSJz2OrDhw/l66+/ljlz5kjevHklVqxYwe6fMGFCeJ6WiIiIiIgCQNGiRWXXrl1Sr149KV++vHZMNW3a1NPFosgIHpEtqWHDhvp/ZFjFH/mnGIlTS/KaXSRLliyeLgoREQU4s07CJRH5PiTaXLdunbRv316aNWum8yCHDx/usXXjcb47ffp0nveGIlyfzOzZs8PzMPJBQXETSoJ8FSVZsmSeLgoREQU4s04iIv8RO3ZsTZxTqFAhXRfy4MGDGmskTZo0ysuC893GjRtH+ev6/ZxH08uXL+X06dP6h/+T/3n1+J482LNMbty44emiEBFRgDPrJFwSkf+IFi2adOzYUf7++28dylqyZEk5cuRIlJcD57u//PILz3vdHTw+e/ZMevToIUmSJJFs2bLpH/7/zTffcAirn3l5/4bcXj1WM+wSERF5Q52ESyLyP5UrV9bgEflUSpUqJX/99VeUvj7Odzt06MDzXncHj19++aWm2B0zZox2LeMP/58xY4Z06tQpPE9JREREREQBDp1S27ZtkwoVKkitWrV0DqRhGJ4uFkVkziMCx7Vr12pqXVO+fPk082qVKlXk119/Dc/TEhERERFRgEuYMKEsXLhQ+vXrp6Md9+/fr6s5xI0b19NFC3jh6nmMFy+eZM+ePcTtOXLk0PuIiIiIiIjCC6s7DBw4UJcGRCBZtmxZuXjxIneoLwaP1apVk++++y5YFzL+j9twH/mP6LHiSpwsb2oLEBERkTfUSbgkosDw4YcfypYtW+T69etSrFgx2bp1a6S9Fs53q1atyvNedw9bffz4sQwbNkzmzZunC3wicNyzZ4+cOnVKGjRoIK1atbJsiy5m8l0xk6WX1A0Haq8yERGRN9RJRBRY3nzzTfn333+lXr16UrFiRfntt9+kRYsWbn8dnO8i4yu5OXhEN3LDhg2DpddFSwD+4OHDh+F5WvJCxutXYrx4Jq9evZKgoCBPF4eIiAKYWSdFixlbokVnnUQUSFKlSqU5V7CkR8uWLXUe5IgRIyRGjHCFM3bhfPfRo0cSP358nvc6EK69jYU7KTA8v35Grk7tJPua7JYiRYp4ujhERBTAzDopTbNREjtNyNwLROTfsITH2LFjpVChQvLFF1/IoUOHdE5ksmTJ3PL8+/bt01GVu3fzvNetcx6JiIiIiIiiGkY8tmvXTlavXi179+6VEiVKaBBJUYPBIxERERER+RTMfdy1a5eu9FCqVClZunSpp4sUEBg8EhERERGRz8maNatmX8U683Xq1JEhQ4YEWw2C3I/BIxERERER+aQECRLI/PnzpW/fvtKrVy/5+OOPdWUIihwMHilUsVJmkQwd/5ACBQpwTxERkVfUSbgkIrJeCaJ///4aRGL4apkyZeT8+fNh3kE438V6kjzvdXPw+Pr1axkzZoyuuZI4cWLL7d27d5eLFy+G5ynJS0ULiiFB8RJLzJgxPV0UIiIKcGadhEsiIltYBxLDWG/fvi3FixeXzZs3h2kn4Xw3ZcqUPO91d/A4atQoXVeldevWcv/+fcvtuXPnloEDuXivP3lx54pcX/CtnDp1ytNFISKiAGfWSbgkIrIHy3ggkU6ePHmkUqVKMn78eJd3FM53a9euzfNedwePWF9l7ty5mibXGiarLly4MDxPSV7q9bNH8uTkTrl3756ni0JERAHOrJNwSUTkCHoPsZRHq1at5LPPPpOOHTvKixcvnO4wnO9i2CvPe90cPJ47d84yFhhrrZiQKte6J5KIiIiIiCiqYQjqr7/+qp1e+Hv33Xfl1q1bllgGw1qvXr3KDyYqgscsWbLI7t27QwSPCxYs0KGrREREREREntamTRtZu3atHDhwQANGXCJny9GjRzWHC0VB8Ni5c2dp2rSpzJ49W6//888/0rNnT/niiy/0PiIiIiIiIm9Qrlw5nQeZKFEiKV26tGzYsEFatmwpv/32mzx6xGHwkR48IoL/6quvpFOnTpp5tUKFCjJhwgT5/vvvNagk/xEjYXJJWrGlpE+f3tNFISKiAGfWSbgkIgrryMktW7bIe++9J++//76Onrx7965MnTrVsg3Od3/44Qee94YimmEYhoQTHoqlORBAZsyYUddY8QWYl4nuakyGRQsEBZelx18hdsnZYTW4m4jI5/n78d+f3x/rJiKKiHnz5snevXu15xFB5PDhwzVIjBMnjhw7dkyCgoJ8egffj6Ljf4SiPUTsCBozZ87sM4Ejhc2rpw/l0dHNcufOHe46IiLyijoJl0REYYEOr8mTJ+tSHAgc06VLJ1euXNFlOSZNmqTb4HwXQSbPex1zeZXdjz76yNVNLXMhyfe9vHtVbi4eJmfONJCkSZN6ujhERBTAzDopTbNREpQmu6eLQ0Q+xJxyd/bsWdm2bZts3bpV1q1bJ0eOHJFp06bp+vVnzpyRDz/8UBOD8rw3gsFjggQJJLLgg7p27ZrkypUrTB8Uxin/999/kiFDBsmenZUIERF5lwcPHsjhw4clWbJkkiNHjjA9dseOHfL8+XMpW7ZspJWPiCiQYNRk1qxZ9e+TTz7R2x4/fqzLepCbg0ckxHG3J0+eaI8m0ufiQzx58qQMGTJEWwZcmW/58ccfy6pVq3Thz1GjRrm9fEREROE1ZcoU6dChg07twHCpIkWKyKJFi3ROijPz58+Xhg0bSty4ceXhQw7RJCKKLFinnlwXoYmKL1++lNOnT+sf/h9W/fv314mrGGuMNVfmzp2rS31s377d6WNHjBih8ywLFCgQztITERFFDqwf1qpVK00Df+jQIV2Q+vLly9KlSxenj8WQKgytQuBJRERRl5Srxk+b9P+4tJeki8IZPD579kx69OghSZIkkWzZsukf/v/NN9/oEJuwtMqick2dOrVer1WrlhQsWFAns4YG67SMHj3aMrmVIk/0GLElVups2vpNRESuwfwZJGNo0qSJXkcdiWBw1qxZWoc6goZYjKpB42pYh7kGUp2ESyIiHmO8eNiqtS+//FKWL18uY8aMkeLFi1sCur59++o8xF9//dXpc1y6dEmuX78uRYsWDXY7rqM3MrQ0tKhYx44dawk6nUFFbV1Z4znINTFTZJS0zUdLnjx5uMuIiFyEesxe/Ya5NUgJj4ZSe3r37i1p0qTRhlXUsc4EWv1m1klERDzG+FDwOHPmTJ2naAaOkC9fPsmbN69UqVLFpeDRTIGLJALWUqRIEWp63DZt2kjVqlWlZs2aLpd36NChMmDAAJe3JyIiigjUY6gXbes38z57Vq9eLTNmzNBEcK5i/UZERF4/bBUTS+1lN8UQG1cnncaKFcuSNMcaWmXN+2wh0cCKFSt0eOuGDRv0D4kEkIgA/0cSHXt69uypC2aafxcuXHCpjCTy/NopOTeibqi9wUREFBzqMXv1m3mfrVevXukQ1+bNm8vBgwe1Tjtx4oTejv9jviTrt//VSbgkInI3HmMiqeexWrVq8t1332lmVKS8BQRuuA33uQLLayDhDYavWsN1ZKazB69RuHBhXdjTdPXqVU2wc/PmTe0NDQoKCvG42LFj6x+FnQbkr146DMyJiCgk1GPnz58PdptZ39mr4xAk5s6dWzZv3qx/5vbII4D5j0gmh4WtA71+Y51ERDzG+GDwiNbTYcOGybx583QOBw7me/bs0aypDRo00Lkazpb4QA/l22+/LUuWLLEkFEAv4po1a+Tbb7+1bIe5IZjDgSGy77//vv5ZQzBZoUIFLtVBREReA1M4UBfeuHFDUqZMaRk9g+kdSKQDqNtQd2IJj0SJEmkPozXMeURyOtvbiYiIfCp4RI8h1p8yofexWLFi+geurkk1ePBgqVy5snTr1k1Kly6tFSUSBbRu3dqyzffff689ixjGQ0RE5AuQ2A3rD2OaRdeuXXU5KmQSX7hwoWWbw4cPS8WKFWXbtm1SqlQpj5aXiIgo0oLH2bNnizuULVtWNm7cqAl2fv/9d219xXMnSJDAsk2uXLl0OI8jCFjtzb8kIiLylJgxY8r69eu1ARQjcJAcbtWqVVKpUiXLNuhtLF++vF7akz59eq0niYiIvEU0IwAns2GoUOLEiTV5jqNKO5BZL4r6+sUzeXn3qpwa8ynXeiQin+fvx39/fn+om8w6KUaSNBI9Zmw5O6yGp4tFRH7C148x96Po+B+unkfAhP4tW7bYTTmO+ZDkH/DDiZUyMwNHIiLymjqJiIjHGM8IV/A4cOBAXTcRyWqSJEni/lKR13h577rc2zpbzrXN7zALLhERUVTWSYnf+khiJE7FnU5EPMb4QvD4888/y8qVK+Wdd95xf4nIq7x6cl8e7l8lt27dYvBIREReUScleLM6g0ci4jHGA6KH50FPnz5lZjgiIiIiIqIAEq7gsUaNGjJ37lz3l4aIiIiIiIj8Z9jqyJEjJX/+/BpAZsuWTdd5tIb1GomIiIiIiCjAg8c+ffrIw4cPdfjqpUuX3F8q8hpB8ZNIolL1JXXq1J4uChERBTizTsIlERGPMT4SPM6ZM0fWrl0rZcqUcX+JyKvESJhCkpZvrotVExEReUOdRETEY4wPzXmMHz++FCpUyP2lIa/z+tljeXp+vzx48MDTRSEiogBn1km4JCLiMcZHgscKFSrItGnT3F8a8jov7lyWa7O+kRMnTni6KEREFODMOgmXREQ8xvjIsNXnz59Lhw4dNGFO9uzZQyTMmTBhgrvKR0RERERERL4aPMaKFUsaNmyo/3/06JG7y0RERERERET+EDzOnj3b/SUhIiIiIiIi/5rzSIEjWlAMCUqQXGLGjOnpohARUYAz6yRcEhHxGBP1wn30ffz4sWzfvl3Onz8vL1++DHZfq1at3FE28gKxUmaRDO2nSoECBTxdFCIiCnBmnURExGOMDwWP+/btk5o1a8rDhw/l7t27uoD8tWvX9L5MmTIxeCQiIiIiIvIz4Rq22rlzZ2nSpIncuXNHr1+9elXOnj0rZcqUYeDoZ57fOCsXf2kmBw4c8HRRiIgowJl1Ei6JiHiM8ZHgcffu3dK1a1f9P5bpwNIdmTNnlokTJ8r48ePdXUbyIOPVS3n18Ja8ePGCnwMREXlFnYRLIiIeY3wkeLx3754kS5ZM/58yZUq5dOmS/j9t2rRy/fp195aQiIiIiIiIfD/bKoaq9unTR5PnfP3115IvXz73lIyIiIiIiIh8O2FOr169LP8fPny4NGjQQEqXLq1DV7kGJBERERERkf8JV/A4aNAgy/+zZ88ue/fuladPn0qcOHHcWTbyAjGTppPUHw+RHDlyeLooREQU4Mw6CZdERDzG+OCwVcByHZzr6J+ix44ncTIVlIQJE3q6KEREFODMOgmXREQ8xnh58Hjz5k359ttvg902ePBgSZ48uQ5ZLVKkiC7bQf7j5YObcuefKZakSERERJ6uk3BJRMRjjJcHjz/88IPEi/e/1j6s/YdkOV999ZUsXLhQgoKCZODAgZFRTvKQV4/uyv3t8+XatWv8DIiIyCvqJFwSEfEY4+VzHv/8809Zvny55frixYs1u+qIESP0esaMGTV5DhEREREREQVwz+P58+clffr0luvbtm2Td955x3IdgeSVK1fcW0IiIiIiIiLyreAxXbp0smfPHv0/sqtu2bJF3n77bcv9mO+YJk0a95eSiIiIiIiIfGfYKoakNmnSRNq3by8bN26UaNGiSdWqVYP1RJYtWzYyykkeEhQ3kSQoWFWTIhEREXlDnYRLIiIeY7w8eOzbt68mTkFSnGTJkskff/whiRL97wA+btw4y/xH8g8xEqeS5O99odl0iYiIvKFOIiLiMcYHgkdkWp08ebL+2bN+/Xp3lYu8xOsXz+Tl3avy5MkTiRs3rqeLQ0REAcysk2IkSSPRY8b2dHGIyM/wGOPmOY8UeF7cuiBXJrWXI0eOeLooREQU4Mw6CZdERDzGRD0Gj0REREREROQUg0ciIiIiIiJyisEjERERERER+UbweObMGdm+fbvcuXPHpe1fvHghBw4c0Hl4z58/j/TyBTIsxyJBMf7vkoiIwuTBgweyY8cOOXHihMuPuXTpkvz7779y+/Zt7m3WSUQUhXje6+XBIzJ41qlTRwoUKCCtW7eWdOnSyY8//uhwe8MwZMCAAZIhQwZp1KiR1KpVS7JkySKLFy+O0nIHklips0nmrovkzTff9HRRiIh8ypQpUyRt2rTSokULKVasmFSsWFHu3bvncPsNGzZI8eLFpUSJEvL5559LxowZpU2bNvLq1asoLbcv1Em4JCLiMSbAgsf+/fvL3r175dSpU9qTOHfuXOncubP2QjoKHvF37Ngx2b9/v5w8eVLatm0rH3/8sVy5ciXKy09ERGTP0aNHpVWrVvLbb7/JoUOH5Ny5c3L58mXp0qWLwx2GbcaOHWvpedy9e7fMmTMn1EZVIiKigAke0SqLyjV16tR6HT2JBQsWdLiOZPTo0TXgTJIkieW2zz77THswEYSS+724eUGuTPmSS3UQEYXBtGnTdDRNkyZN9DrqrQ4dOsisWbPk2bNndh/TrFkzKVq0qOV67ty5pWzZsrJ161bue5s6CZdERO7GY4wXB49oWb1+/XqwihJwPSyB4M6dO/Uye/bsDrdBRX3//v1gf+Sa1y+fyfNrpzRAJyIi16Aes1e/PX78WEfPuAJ11759+yRHjhyhbhNI9ZtZJ+GSiIjHmKgXQzzETI6TLFmyYLenSJHC5cQ5CD7Rkothqzlz5nS43dChQ3WuJEVMlh5/hbjt7LAa3K1ERDZQj+XLly9E/Wbe54quXbvKw4cPpWPHjg63Yf1GREQB0fMYK1YsvbTt0UKrrHlfaFD5VqtWTZPnjB8/PtRte/bsqUkKzL8LFzjchYiIIg/qMXv1m3mfM4MHD5ZJkybJwoULtZ5zhPUbEREFRM8jKkPMYcTwVWu4njlz5lAfe/fuXalatarEiRNHVq5cKfHjxw91+9ixY+sfERFRVEA9dv78+WC3mfWdszpu2LBhGjwuW7ZMypcvH+q2rN+IiCggeh7jxYsnb7/9tixZssRyG4bnrFmzRqpUqWK5DXNDdu3aZbmOnkMEjjFixNDAMWHChFFe9kASI0kaSVGnh2TNmtXTRSEi8hmox7Zt2yY3btyw3LZo0SLJmzevJtIBzE/E8hzW8xS/++47+fbbb2Xp0qVSqVIlj5TdF+okXBIR8RgTQD2PgJbVypUrS7du3aR06dIyZswYSZMmja75aPr+++916Y6DBw/KixcvdKgq0pljqOqePXss2+XKlUvX0yL3CoqTQOLnLiNJkyblriUichHm4o8aNUqziGPuIpajQiZxDEM1HT58WNd+RJBZqlQpXdaje/fu0qdPHwkKCtLAEhIlSiRFihThvreqk4iIIgOPMV4ePCIF+caNG+XXX3+V33//XSvH2bNnS4IECYIFheYCyZgvgiE6efLkkZEjRwZ7LgSgNWoweYu7vXp0Rx4d2iDXrhVz+3MTEfmrmDFjyvr167UBdMKECZocbtWqVcF6ExEUYlgqLs1hrbiOehF/JvRWop6k/9VJ8fNVkKD4bNQkIvfiMcbLg0dAayv+HEFQaEqcOLGlJZaixssHt+TO+oly6VI77nIiojBAnTVo0CCH9yMotK7TQtuWgtdJsTMVYPBIRG7HY4wXz3kkIiIiIiIi38HgkYiIiIiIiJxi8EhEREREREROMXik0L8gseNL3OwldO4OERGRN9RJuCQi4jEmABPmkHeLmTStpKrXV7JlyyYiRz1dHCIiCmBmnURExGOMZ7DnkUJlvHoprx7f0zU2iYiIvKFOwiUREY8xUY/BI4Xq+Y2zcvHnRrrANRERkTfUSbgkIuIxJuoxeCQiIiIiIiKnGDwSERERERGRUwweiYiIiIiIyCkGj0REREREROQUg0cKVaxUWSVjp7lSqFAh7ikiIvKKOgmXREQ8xkQ9rvNIoYoWPUiixY4nQUFB3FNEROQVdRIREY8xnsHgkUL14vYlub16rJxomZN7ioiIvKJOSlblc4mZLD0/DTd5/fq1PH/+nPvTy8SMGZON91GMxxjnGDxSqF4/fyJPz+6VBw8ecE8REZFX1Em4JPdA0HjmzBkNIMn7JEmSRNKkSSPRokXzdFECAo8xzjF4JCIiIgpAhmHIlStXtHcrY8aMEj06U2F402fz+PFjuX79ul5Pmzatp4tEpBg8EhEREQWgly9faoCSLl06iRePc0m9Tdy4cfUSAWSqVKk4hJW8ApuYiIiIiALQq1ev9DJWrFieLgo5YAb1L1684D4ir8DgkUIVI1FKTUyA4SxERETeUCfhktzHV+bTnThxIsJB1MmTJ+XZs2fhem1PJBXylc/GX/AY4xyDRwpVULzEkrBITUmZkhU1ERF5R52ES/JPp0+fdpik77333pNLly5F6Pnr1q2rgWBY4bXPnz8fodcm78djjHMMHilUr548kIeH1svt27e5p4iIyCvqJFySf/rkk09k7dq1ni4GBSgeY5xj8Eihennvmtxa9oOcPXuWe4qIiLyiTsIlBYYnT55oRlh7sLzIuXPnQvRUYumRgwcPyqlTp8I0zNUcmnr//n25efOm3W3wWrblcfZ62P7QoUO6jSvlJ8/hMcY5Bo9ERERE5HVmzZqlWUZLly4t1apVk6dPn1ruW7dunWTNmlXeeecdyZQpk3z99deW+/r16ycfffSRPgb3bdy40eWhqS1btpR8+fLpc7dr1y7Y/YMHD5b8+fNL7ty5pWPHjk5fD9lscVvhwoWlYcOGuo2ZpCi08hN5My7VQURERESWXjLbnrWkSZNqoIPg7fDhwyH2VJEiRfTy2LFj8ujRo2D3ZcmSRZIlSxbmvfvw4UPp0KGDBlnFixeX5cuXS40aNfQ+lKNVq1Yybtw4SZ8+vfZOIjjDkFcEatOmTdPePMyPRCDXv39/fR5XIMfDhQsX5M6dO1KiRAnZsGGDVKhQwfJe0FOIqTxvvPGGDBgwQN+bo9fD/sBzXb16NVjiG2flJ/JmDB6JiIiISP3+++8aFFlr1KiRzJgxQy5evChFixa1u6A9NG/eXLZv3x7svunTp0vjxo3DvHcReCHTOwJHqF69uqROnVr/f/ToUQ1wu3TpEmxJC3OoaadOnWTixImSNm1aiR49epiypDZr1swSMNepU0d27txpCR4R4AECxsyZM2tQiP87er1s2bJJzJgxpWTJkvoc2I+FChVyWn4ib8bgkUIVPWYciZUul8SPH597ioiIvKJOwiVFjjZt2kjt2rWD3YZACjJkyCC7d+92+NgpU6bY7XkMD5x33Lt3z3IdcwnN58Z9+Nu7d68GZ9aQo2Hu3Lly+fJlSZgwoZa3Xr16Lr/u3bt3Lf9H7yN6XE0xYvzvtBk9iZizGNrrxYkTR/bs2aMBNXowMUQVPZOhlZ88i8cY5xg8UqhiJs8gaZv8ILly5cLqSNxbRETk8TqJIg96z/BnD4Ihc4iqPf93ruAeOXPm1GCsR48eurzGpEmTdHgnZM+eXfLmzStNmjTRuYeJE//f0i2Yixg3blx5/PixrF+/Xm/v1atXmF73m2++kWHDhmlQ+Oeff0rfvn1D3T6010Owi2Q6iRIlkrfeekt7JxGclipVymH5rQNUino8xjjHbygREREReQUM9USwheGfS5Ys0UQy3bt3l6ZNm+qw2VixYmmvH+4bMmSI3m/2UP7zzz86tPWnn37SADB58uSa9Gbq1KmW58+RI4cGwY60bdtWH/vs2TPtUcTwVDOYjR07tmU7BLB4ntBe78CBAzq3EeXF8FaUFcl/wFH58RxE3iyaYQ5UDyBIwYxWHvxYcYCi4LL0+Mvy/2dXT8rVqZ10GMYHc0Omyj477P8mrxMR+QJ/P/778/tD3WTWSWmajZLYabKzDoogJG5BzxiGZoYWUAUKBIQrV67US2/Bzyjq+Pox5n4UHf/Z8xjArINEky/9SIiIiIjcxbZ3kYhCYvBIRERERAEPy4EQUegYPJLbezDZe0lERERE5H+ie7oARERERERE5P3Y80ihipUik6T7bJymlBYJmTCHiIgoquukGAlTcKe7UQDmTvQZWEuSog6PMc4xeKRQRYsRS2ImTccsbERE5DV1ErkHFqjHMhI3btyQlClT6v/JewL658+f62eDZUuwRAlFPh5jfCR4vHDhgly7dk2zXLmaWjY8j6Gwe3H3qtzbNEPOtEHPIxERhQUWCT969Kiu8YblECLrMYFWJyUu21hiJknj6eL4vKCgIMmQIYOun3j27FlPF4fsiBcvnmTKlEkDSIp8PMZ4efCItWsaNWokK1as0EVYz507J8OHD5eOHTu69TEUfq+fPpRHhzfInTt3uBuJiMJg+vTpumB4unTp5PLly1KqVClZsGBBqA2e4XlMINZJCYvX9XRR/EaCBAkkR44c8uLFC08XhewE9zFixGCPcBTiMcbLg8cBAwbIzp075dSpU5I2bVpZtGiRvP/++1KiRAkpWbKk2x5D/rfmJNeoJCJvduzYMWnRooWMGzdOPv30U22AQx3VtWtXvc1djyFyV5CCPyIiZzzaBz558mRp1aqVBoFQt25dyZ8/v97uzscQERFFpWnTpkmaNGk0CISkSZNK+/bt5Y8//pBnz5657TFEREQB0fOI4TiYs1i0aNFgt6MHce/evW57DKDSta547927p5f3798Xf5K/398hbjs44F2Ht79+9jjE7dgn1re/fv5ULx8+fOjS9uZt9soTnrI4EtbtKXLZ+6x9gaPvI/kv8zgR2dklUSfZq6seP34sx48flwIFCrjlMYFSv5nHfbNOwiWu++P7JCLP8PVjzP0oqt88Fjzevn1bL5MnTx7sdlw373PHY2Do0KE63NVWxowZxd8lHuWe28uXL+/y9u56TUe3OxLW7Sny+PJn4ctlJ9fdunVLEidOHGm7DHVSvnz5gt1m1l2h1XFhfUyg1m/XZ/XQS/5eiYjHmKit32J4Mj20mQDH2pMnTxymIw7PY6Bnz57SuXPnYGvmoCJGpewoLTWid1S+yOrqa4kKWHbud35fvB9/p56BnjlkLkQm08iE+speXQWh1XFhfQzrN9/B3zz3e6B8Z3y13L5e9ntRVL95LHjEB4O0w5cuXQp2O67jjbvrMRA7dmz9s5YkSRKXyokvjq99eUwsO/c7vy/ej79Tz4jstPfIBo7lD6yZdZej+io8j2H95nv4m+d+D5TvjK+W29fLHj2S67fonly35q233pIlS5YEW9tqzZo1UqVKFcttJ0+etMxndPUxREREnoQ6aevWrXLz5k3LbYsXL5Y8efJI+vTp9fqDBw9k8+bNeunqY4iIiAI22+qgQYN0qQ0Mu0FAiMypqVKlks8++8yyzbBhw6RJkyZhegwREZEnffzxx5I3b16pU6eO1lkDBw6UiRMn6hxF06FDh6Rs2bJ66epjiIiIAjZ4RBKW9evXy7lz52T06NGaKACtsFiw1oSFa4sUKRKmx7gDhgL169cvxHBXX8Cyc7/z++L9+Dv17/2OOYqoqxAcjhkzRg4cOCArVqzQwNCEIVFvv/22ZWiUK49xB373PIP7nfs9UL4zvlpuYNmdi2ZEdj5XIiIiIiIi8nke7XkkIiIiIiIi38DgkYiIiIiIiJxi8EhERERERETeu85jVHv58qVmtIsRI4Zms4sWLZpbHhOe5w0rLBJ9+PBhSZw4sWTPnt2lx9y4cUPXC8uaNWuINS3xfLt37w7xmAIFCuhruHvB0hMnTkiaNGkkQ4YMoW57584dS9ZBa8WLFw8x6Toszxte169f18RMWbJkkZQpU4a67X///ScPHz4McTv2J/aruV7bmTNnQqzFg+Vn3O3y5cty+vRpKViwoMvrFJ0/f17fc86cOR0+xpVtIurUqVNy5coVKVmypC6a7szz58/l6NGjWh6sk2f7G0TSEXxfrKVIkUJy587t9rIfPHhQvwelSpVyuu3OnTu17Nawlp/ten6vX7/W3/+rV68kf/78EhQU5PZyY+r7nj179DhWqFChULe9evWqLqFkT9GiRSVu3LhaZiw5YQvfG2THdqdr167p9wXHOlePX1jiCd8ZLKSMx4V3G28Snt+mt/zmjx8/rr8bJMBzJcHGs2fP5NixY/p54/fiym8ex/BcuXK5tdz4nqPOwu8HZXf229yxY4e8ePHCpd98WJ43PHDswWvg9+rsWIjfF47L9hQrVkzixImjx6dt27aFuB/73Fn9GVbYh//++68ex5FU0RX4fuH3jMegTg/vNhGF19i3b59+5li73BX4DeL7nC1bNl2yzvZcD78FW6VLl3b79wZLCGH/YOmg5MmTOy0z/qwhIViJEiXsfr9wvor3F1mL25vlQQJO231oDb+5LVu22L0PnxfOMQD7wXpJJdvzPXd5+vSpfr7YLzjXdTXGQPlwro/jB/Z7eLdxyAgAW7duNdKnT29kzJjRSJUqlZE7d27j2LFjEX5MeJ43rObNm2ckTpzYyJEjh5EoUSKjTJkyxq1btxxuv337dqNcuXJGypQpjcKFCxtx48Y1WrRoYTx//tyyzZEjR5AkyShevLjx9ttvW/527drl1rL/+OOPRpw4cYw8efJoOerXr288ffrU4fZLly41okWLFqxM+Lty5UqEnjesXr9+bXTs2NGIHTu2kTdvXr3s0qVLqI9p3rx5iHJHjx7daNCggWWb77//3ogfP36wbSpVqmS4044dO4z3339fP398xps2bXL6mMePHxt16tQx4sWLp99h7NNffvklzNtE1IoVK3R/JEuWTMtu+7nbevDggdGpUycjadKkRoECBYzUqVPr5/Xvv/8G2658+fL6G7Xe771793Zr2adOnWoUKVJEy4LfqytQ3pw5cwYr17hx44Jtc+jQISN79uxGmjRp9FiTKVMmt/5O8V0fMWKEHl+SJEliFC1a1OljFi1aFOK7ni5dOiMoKMi4fv265bPBZ1iwYMFg2y1fvtxtZd+4caNRqlQp3Y+FChXS72S7du2Mly9fhvq4adOmGQkSJNB9j8t33nnHuHfvXpi38Rbh+W16y28ev/ESJUrod++NN97Q3/6yZcscbn///n3jiy++0N8Zvlv47PPnz2/s2bMn2Hb4rtn+5vv16+fWsu/fv9/Ili2bkTZtWv1tZs6cOUQ5bCVPntzIlStXsHJNnDgxws8bVmvWrNE6ImvWrFom/H7Onz/vcPsFCxaE+M2jfDFixLCcj9y5c0d/83gu6+3+/vtvt5Ubv8FvvvnGyJAhg/4umzVr5tLjJk+erHUv9j0u3333XT1GhXWbiLh06ZLRvn173W+xYsUyBg4c6NL5H357OO7je473PGzYsGDbTJ8+3YgZM2aIz8edZT948KDRqFEjrYfwGaNczvTq1UvPWa3LVLt27WDb4FiNcyecz5nnWgMGDDDcCfVEjRo19HuOsh84cCDU7Z89exZiX+JcGo/94YcfLNvVq1dP94f1dh06dHBbuW/fvm18/vnnemzEsS5FihR6joFzgtDgd4zfIN4vft/4neP3HtZtnPH74BEVIE5scFJhflmrV6+uH0JEHhOe5w2rc+fO6Y/pp59+slScOIB8/PHHDh+DAwl+LKYTJ07oF8T6B2kGjxcuXDAiy+bNmzUQ/Ouvv/Q6Xgs/tD59+oQaPOL9uvt5w2rChAl6kEYlDjhZR7n++OMPl59j586duo/NcprBI36wkQknIvPnz9fP3dXgsWvXrlo5Xb16Va+jYsA+tg7CXNkmorB/Vq9ebaxatcql4PHs2bPakIDfIqCBpHHjxnrC9erVq2DBIyqyyIQTGuyLn3/+OUzBI36vjuA95MuXTxtHEOQBTpZwMokKzh3wPJ07d9ZGry+//NKl4NEeHPdq1apluW4Gj9u2bTMi87uOxjLrExycrFhX8LaOHj2qJ7yTJk2yVNAInFu3bh2mbbxJeH6b3vKbr1mzplGyZEnLbxgn1Dj2Xrt2ze72p06dMkaPHh3sN//RRx9pOc3fCERGsGj728RJfcOGDS2vi2MPAuAXL144fBzq4lmzZrn9ecMCQR6CbxyzzGNA2bJljYoVK4bpeXBCi4ZK6+fFb97djdDWjh8/bgwaNEjrhsqVK7sUPOJkGw1baBCCmzdvanDetm3bMG0TUTgvw7nc3bt39RjuSvCIhj3rDomVK1dqo7T1eQXqENQlkWnu3Lm6b1D2sASPqHudvT80GJ08eVKvb9iwQT8H6/cXUWjwwrnl7t27XQoe7cHnhjrBPBaawWObNm2MyILv5G+//Wap69FBguAbnSahqVChgv6Zj+vevbv+3vH7DMs2RqAHj3/++adWeJcvXw4WgOBLtHfv3nA/JjzPG1ZDhgzRysa6JR2BDVqZwtIKjsq1SpUqIYJHHMzQoomg1N3Q21msWLFgt/Xo0UNP7J0Fj4cPH9bA7cmTJ2553rB66623tMK2VrduXa2sXIWDClq+rYMYBEcI/vft26efgXVvsLudOXPGpeARJyj4jqFCtoYDFFpJXd3GnRBAuhI82oPKB4/FSaYJFRh6khHQR2aDCYQ1eBw5cqSebNk7WcbIBryX//77z3IbKlnchl5adwtv8IjjHcq0ZMmSEMEjTjIQbISlUopoMGJ9QmvL7LWwNmrUKO1dM0cvuLKNtwjPb9NbfvM4EUMdisYuE4JC9Prgd+QqtJjju4bGVuvg8auvvoq03zzqTrwmGixMOMnHbTh+OYJ9iuAXv3mzl94dzxsWaBRBz5f1OQR6e/EaqDdcgfJje+uRBGbwiHOjqPjNuxo8fv3110aWLFlCBC1opDDrYFe2cSdXg0d7MBKlZ8+ewYJHjHzDdwZ/kXmMwjlZWIJHnEuhfsB32F7jB3obbXvrENQgMHM3s54KT/CIRv8PPvgg2G316tUzmjRpor8FNGRbN15FFoz6wXtwNPrw9OnTej8aGUz4HSJmQM+6q9u4wu8T5uzdu1fSpUsnadOmtdxmjrnGfeF9THieNzxlL1y4cLBx63gNjPe3NzfQHsxDwPPYmyv58ccfS+PGjXXs+ueff67zSNwFr4n5T9ZQdsz7wxh9R1CGmjVrygcffCBJkyaV/v37u+V53VF2Vz/Xx48fy6xZs6RFixY6p9EaPrePPvpIqlSponO/Jk2aJJ504cIFuXXrVoj3i3mm5vt1ZRtvsWvXLp1/kz59+mC3jx8/Xlq3bq1zBjEPFPP7vAG+361atdI5dZUqVdI5tibsW8xBRHlN5pwQb9rvEydO1GNh9erVQ9zXtm1bad68uaROnVoaNmwod+/ejbRyYA7X/v37Q50X7ui3jd8s5t25uo23CM9v01t+85j3hQZs69fA/DvMvwnLa+A3jzlM1nUxjB07Vn/zeD7M48W8dHdB+TA3E89tMueEOit737599TePuVOoB7Cv3fG8YSk75glaz18N67kLfvOY//Xuu++GuK9Nmzb6m0f99sknn4SYexrVHP2eMffQnLvtyjbeAPO78X2xPcZhTvL7778vtWrV0vphxIgR4g22b9+u55gVKlTQOmLOnDnB5vIdOXIkQudaUQH5QXCswrHE1pw5c6Rly5Y6jxK/040bN0ZqWXCsw+eLc2N7zP1mvU+R8wS/d+vYxdk2rvD74PH27dshJvYiCUfChAn1vvA+JjzP646ym9ddfY1+/frpROSvvvrKcluCBAlk2bJlejuCGXwh586dq5Wau4Sn7DjhRxIRTMxHMpzFixfLkCFD9MQ/Is8bFjigYQKxvddAQh+c7Dgzf/58rXQQPFpDQwDeF5KfoALAe8NJRGQfcEJj7jN779f6u+5sG2+A/frtt9/K119/HSzpBg76mNiOk0ckEkIChzp16siDBw88Wl58/jhBR7nOnj2rwQkaFszvGPYtKgrbCfLetN/R2PPHH3/od926kQtB77Rp07RBB8lLcJKA40z79u0jrSz43O/fvy8dO3Z0uI0rx4/IPsa4U3h+m97ym3fHayBB1eDBg6VHjx7BkmuhMdT6N49Gl7p169pNahbesttLGOKs7MOHDw/2m0dghUbciD5vRMtuJilx5TVQP6JxFCfN1o2j2P84FiCQwW8ex2MkHvniiy/Ek/zlN48kSjjOItEO6gkTvttoNEPDFpLk4biLY+G8efM8Wl4kAkRjKH6jaNzv3LmzBpJmgIKGRNR13n5egYYS7POqVasGu/3DDz/UYB6BJRL+IEBGAI/bIgMSRP3www/Su3dvh0lzzP1mm3TI3rE9tG1c4ffBIw5oCAhs4TZH2YVceUx4ntcdZceBG1x5jTFjxsj333+vrSPWGcmQsalGjRqW62iVbdeuncyePdst5Q5v2d98801t2Tbhx4oeSOtyRXSfuFJusPcaOCF2JdMVDjZokbXNoPfOO+/ogd66Vwa9StatcVEttPdr/V13to2n4USsWrVq8t5774VoBGnUqJHEjx9f/48eih9//FEbTjZv3iyehBMBfKcA2QgHDRqkLbV4L6EdY7xpv//55596EoATSWvo/W3SpInl+htvvKEnNGhYQYZqd8NJ+bhx42TBggWhZjB05fgR2ccYdwrPb9NbfvMRfQ00MuI3j96WXr16BbsPJ6lmRkX89keOHKknsvaygYa37OH5beJ3YjayoGcOjV0IsMzex6j4zdt7DfO6K6+B3zAa3mwbR7Gf0dNoQu9Yt27dNIhB4OMp/vCbR5CFRlD0gi1ZsiRYtlBkVbXO8FmvXj09d3Ln+Vx4YCSKmQkf501o4MF1fH985bwCZZk5c2aIhhIzeDRXMkB5R40apY1Bq1atEndD4yvO2XFc69Spk8PtzH1qO4rQ3rE9tG1c4ffBI4aGILW89cELLWMY+ml7ch+Wx4TnecNTdrTYWDOvO3uN3377Tbp06aIHbutA0REMK7N9rcgoOypODF9wlW253PW8juB50ANq7zXMFM2hQc8iehLtDXGIiv0eVmaae3vv1/yOubKNJyHYQqsfhrug5dtZenKctNl7P56G7wKY5cL3DT1p1j2kGJqJ3jxv2O9mQwmG3rmS1h7vD+W3TW8eURiiNWDAAB2pgKG/ET2mRuS4G9XC89v0lt+8eTwNz2ugh6VixYpSpkwZmT59eogTO2e/LXeUHSNRMFrAhJMx9CqGZf/Y+82743mdlT0i32/85hG0u7LMBN4fTko92ZPk6795BI4YCozRYuvWrXNpiSlPn1eEVvea5ULPF0bqeet5BaAx0t4oMnvQeIJRfe7e71hOA/UazuPRQBpaB4ajYypGX1h/151t4wq/Dx5xYoMTsA0bNlhuw0kGIuxy5cpZbkMvBLqeXX2Mq88b0bKjix+9JNavgQ8Y46sBlQzKbj2v4Pfff9fWCQxFrV27tt31y2ytXr1a54O5C8qO57RuVULZsW/MIYWoUFB2swXEtlyYr7l+/fpg5XLled1R9qVLl1qGD6KBANdxuwmfib117FCxYu1JtIbbsn1/OEnAUAR37ndXIMA15//g4I31FNGaaUKwgkrKfL+ubBNVMBwKB1MT1m3CSSTWGkNLq9mTZ8L3BN8ja/j+4LON6v2ONd7Mda/s/QbRYonymycHeF+4ju+e6e+//9YADL3YUQX7D79TNI5Zw5ql+A7Yayhx9P7Qw+rOdR7Ro9SnTx9ZtGiR3e8ivqcouxmAYxv8bq0DWBw/sG6ZOU/WlW28hau/TW/8zWMYP9bTs34NtLCjrNavgSF51uvYobEIvw30uMyYMSNEYxGCFdueLrM3wF2/eZzMIWDFCb1p+fLl2qteuXJly20YSeDsN4+eAHP9SVefNyKwb3HiaD3vG99vzIHEZw54PfxubPMIYP6fo8ZRR+8P9aGzNQHdCedlKLs5RBnvF9dR31q/X/TWmcG7K9tEBRxjUQ6zzkI9hSHYKAt+e1hL3Nl+R523adOmKK/f0LOPaUeOyoXzawxhNcuFIAjfaevfPzpf8H2P6vMKHHdwbmFrwoQJ2lBiu5b4ixcvQqzRjPod5+Hu3O847uFYhzKgLPYCR8QI5rxc/H5x7Lbep5gugsDQ3KeubOMSIwAgIxJSeWOphfHjx2s2xL59+1ruRxYo7ArrDG/OHuPqNhGBLKtYxwwZELHO0vDhwzVdsPWSEcgcZZ2JDZm3kMGuW7dummnT/LNeJwppeVu1aqXpl5EhEZlFkX3NnRkckdIZ2cuqVatmLF68WNO+I5sTMtKaFi5cGCzDG5YgwXqKyCiFDHxVq1bVfWoumeHq87ojHThS/mP9IeyfTz75RNMYI6OWaejQoSGWFcH3CGs4IfurPcgA+O2332oa6hkzZujaQXgvjtLShzeDIT5vZEPDvv3111/1uvUaXshQZ71kyNq1a/V7hSyT2PdYaxFLEzx8+DBM20QU9i/KiuU3UHZ8vrh+48aNYPvQzMSGbIVIYY+MbevXrw/2fTczCWLZBfx+kPIa641hGQekBrdef9MdkCEYr4sMj8jQZ5bDev8g06K5ZAiyG2LtQGQ+RNYzLDWDta5sl5zB7xiPw3b4beP79dlnn7m17Dg2oKzYJ1jjzCy7mSnYzKJom4kNa2VifSh72QiRURLLDeBYhe87MnTi+xOWbG7OjB07VsuFY671Z49sxiYsFWK9ZAhSk+O7jyyAOP7g94jU8PhOm1zZxpu48tv01t/877//rsdRZLPFMQtL09iufYulPPBdMo9vOGYiazUyK1t/7mbWcKS4xxrG+H7gN4+smTh+h7bEVXjgt4611/CdxjIGyJ5sLt1lQv1lLhmCfYg6DdvjN49jgb117Vx53ojCsjpYw3TOnDm6lAHW8ETmZxOOufjd2C4lhLoNS2PZy5yJYysyu5u/eSxzge+PufyFu5ifN7Kuv/fee/p/6+VBcN16yRBkH8X3BWtk4/fcv39//T0ju7vJlW3ckanULDv2IZb+wf+ts3/iM0fZzUy1+C5gH+K7bP1dxzmK9WeJ7Ks4V8E5HZZdwXHZXP7CHbBcEV4XxwSUD5licd06qznOK62XDMH5Dc6TkJF3ypQpmqkZy9DgHM46Ayq+e6gfUH5ktce+cec50cWLF7Ws5r7F9xHXrVdKwHq2tsuKYLkznEvbO+5fu3ZN15bGEh44xowZM0bLjQzA1hn2IwLnbFgOEN/zf/75J9jn/+jRI8t2OGa2bNky2O8QmcHxu8bvG8ds2/U1XdnGmWj4R/wcWgkw/w+t9mjJx5hwZAMzo3i08pQvX14n9GKOnSuPcXUbd7SiYd4i5mokTpxYPv30U81Gat36j7lFGG+N3peBAwdqeWxhzhEmUgM+ckx4R48Gnh+tnkhkYT0fzx3Q0jRs2DBtbULrY4cOHbS12ITWsZ49e+oYeNyPlhy0rqxZs0Y/E7TgYLK9bcufs+d1B0z2x+Rk7F/sF8zdMHt7AcMjkSl17dq1ltvQkovyTp061e6+RKvUL7/8or0ayCqIzwtlN+fjuQNarLFvbH322WfStGlTS6IWtO5Pnjw52Gfx66+/assn5sBifoJtD5Er20QEsiOiJ8EWhiSare6Ym4tWbHzP0VKIeaOOhm2b80DQU4nruMT3BXMxkPnTnb755hu7iY+mTJliyYyH3y0S9Zit9vge4DNAzwRGE2C+EFoZraEHBb8JtDzj/yg73rNtD2tE4HuBYYC20EOL7yla8NHyiXllmFNqQqIPfIcxPN4eHF8wbB69Fzj+4DuI701k73NkqsToC/N3jNfFcB+z1R69C5gjiV5/DJ3C52Hb4urKNt7E2W/TW3/z5rxZDD3F96xs2bL6fbI+JmLIHjKpIjMxWtkdJUSy/ozxueM3jyQieCx+e/Xr13drufF7xGuiBR91KoaV4bdp3ROK30uDBg0sw97Qq4RjAuY4YvgY5mPj3COszxtR6J3C3G+M7MFvHAlYrBP3oK7C62L+uHWiEMzzQiIUR/OucJxCfY5ee9SB+O1ZZ4t2BwxVtoXP2EwQY9YLGAVk9uhilBPqRcwZRP2B75RtT64r20QEenvt1TvI84DPAlasWKEJoFauXKlDIPH9sM7AbUI9gOOf2dOOuvOff/7Rc0/8TnEeYpsQJSJQV2HOui2c85rJGFEGlB/fAcBQ659//ll7tfBeSpUqpZ8L5sNbw28a568YzYXPq3v37i5NEXIVRiThPN2W9fk+9iVGVqC8JhyTcJ6Hcyp79e3Zs2f1eXEeihEUGA2EutTZEPqI7nOzbMjQDnhN1HnYbyac3+N947uBcwp8Rrb73ZVtQhMQwSMRERERERFFjN/PeSQiIiIiIqKIY/BIRERERERETjF4JCIiIiIiIqcYPBIREREREZFTDB6JiIiIiIjIKQaPRERERERE5BSDRyIiIiIiInKKwSORE/v375cNGzZ47X7atm2bHD161K1lxSLba9asCXbby5cvdSHiOXPmyKlTp8TbYYF3LOBLRET2PXz4UBcLf/TokVfuouvXr8vy5cvdWtanT5/q89y7dy/Y7cePH5eFCxfKqlWrxNvdvn1blixZ4uliUICK4ekCEHnaiRMnZPfu3fr/6NGjS9q0aaVw4cKSMGFCvW3mzJkaiFSoUEG8zdWrV6V27dpaPneW9c8//5Rly5bJO++8o9dfv34t5cuX18o2f/78kjx5csmWLVuwx1y7dk3Wr18vH330kXgDnGB8+OGHGgjHiRPH08UhIopyOHbPnTvXcj1RokSSJ08eyZo1q6UO+fjjj+XMmTMSP358r/uEvvzyS8mXL59Ur17dbWW9e/euPs+BAwckceLEetu4ceOke/fuUqlSJcmdO7dUrVo1xOOWLl2qZXnjjTfE01BulDdGjBi6b4iiEoNHCnh///23fPXVV1KvXj2taI8cOSKXL1+WKVOmSK1atbx6/wwaNEiDx8yZM0fq66CS3bp1q7Z2Jk2a1OE2qJC9JXhEsJsmTRoZP368dOzY0dPFISKKcs+fP9fjcunSpSVTpkxy584dHUHyySefyMSJE736E/nvv/+0EROBXWRDPfHNN99It27dHG6DeqR3795eETwGBQVpWXv27MngkaIcg0ciEYkdO7YOYwHDMKRZs2bSvHlzuXHjhmX/PHv2TCszBFDFihWTlClTWu47ffq07Ny5U/+fIEECbZ00W3atXblyRZ8DvZpFixaVuHHjBrv/2LFjGryi9/PNN9+UWLFiOfx8MIRn6tSpliE91txRVhPKM2/ePEugDQi0Y8aMGawl1xwua+7H7Nmzaw8lgkoEuKZbt27J6tWrpX79+tpqikB9y5Yteh3lunDhgvZ4nj17Vns6S5UqZXkvxYsXlxQpUgQrHx6P+9GibrtPmzRpIqNHj2bwSEQB7YsvvrA07CF4xOiUGjVqSKFChSzbYDoCjvcIMgsWLBjs8eZxHUELGitRP1nXAfDixQvZvn273L9/X0fvpE+fPtj9OJ7jfsD9qVOnDrXMY8aMkQ8++MAyCshaRMtqHVxjpM25c+fk5MmT+lhsnytXrmDbYSgrRrPs2rVL602MUsLIFvTqlitXThsqTStXrtQAM2fOnHp98eLF+n5RHoxyypgxo+TIkUP++usvqVOnjly6dEmnnmTJkkVH9tiWb8eOHbpPUa506dJZ7mvQoIG0b99eNm/eLGXKlAl1XxK5E4NHIhvRokXTA/r06dM1kAFclihRQis7VIA40KOCQGsuoOJZtGiR/h8H+U2bNkmbNm1kxIgRluf9/fffpWvXrvqYV69eadCDYaaoEFBBIFhdu3atvg6eDxUx5jSgkrEHwRqeB8GVNXeU1XZYL4I7MB+HYNC6Msbr2G5TpUoVrSzRUmsdPOL50BJes2ZNrYT37NkjjRo1kgkTJmiAiOGwCBJnzJihlSsCTATTCFAxJwWVON4f/Prrrzp056233tL9hWFNs2bNspwQYQhSu3bt9KQAwSwRUaDDqIwkSZJoo5t5rERwiYZFBDCoWz777DMZOXKk5THmcR3H2X379mkjHeoVM0A8f/68Pi8CPTzHoUOHpEWLFtKrVy+9H0EW6hm8HqYRIIgcNmyYfP755w7LieM/trEV0bLaNrRiewSGe/fu1boM9ZJt8IiA+/Hjx5ZtULcheEQPLnpHq1WrZtkW9XyrVq0swSPeNwJc1F+o71EfxosXT+tB/B/vBUEupn0gGPzuu+/0cWhANT8r3I85/Hgu1HmAfY3GYbw+g0eKUgZRgPv555+N+PHjB7tt6NChRvTo0Y3Hjx8b3bt31/9v2bLFcv9HH31kVK9e3eFznjp1ykiQIIGxfft2y21p0qQxpk2bZrl+8eJFY8eOHfr/vn37GiVKlDAePXpkuf/zzz83Kleu7PA1BgwYYBQsWDDYbe4q68CBA42SJUtarq9evdpwdriwt83kyZON9OnTB7tt27Ztut2DBw/0+tKlS/V6v379gm3XpUsXIygoKFi56tevb9SuXdtyPUWKFMbMmTMt18+fP2/s2rUr2PPEihUr2DZERIHiyZMnenydNWuW5barV69qPfHLL78YJ06c0Ptbt25tvH79Wu9fs2aNES1aNK2j7Hn16pVRt25do1WrVpbbevfubZQrVy7YNosWLdL/4zVQx27atMlyP47rceLEMY4fP273NS5duqTl2rlzp+U2d5X1ypUr+jwHDhyw3JY5c2Zj/Pjxoe5Le9ugjlqxYkWw2/Lly2f8+OOPluupU6c28uTJY9y7d89y25EjR7QMbdu2tbyXlStX6ueC8kGPHj2MSpUqWR7z8uVLY/HixcFeC4+33oYoKrDnkej/ZxLFcBUMWcVQGLRionXTHAKJ1lL0bpkw5Of7778Ptu+ePHmiQ1LQ+4XnwzAWDMMsWbKk3o/nQvIW3IfeNLSCmi2hkydPlnfffVeHoKIM+EuVKpVMmjTJsr2tmzdv2p1/6I6yegL2t60iRYoEKxPey88//2y5jn2KFm5zH2E4EP6sodUW+4qIKFAhKzdgBAdGwaDuQc8XphEAegAx6gYwDNMcJWLdW4djLXrJ0EuH+skcgmoei3GcxYgaDK3EsE6M4AGMsMH2qG8wBQL1m5n0BUMu7Y2uMY/Z9uq4iJbVEzCyCFMrQnsvqN+QdwEjZVAnY58i2yymu2D0DXo7rUfxmPuH9RtFNQaPRCI6/BNDV1DhYbgnKjvrg3SyZMlCzJFEum8ThptgCAsO+JjrgIM+hrbgwG9CgIghlGPHjpWyZcvqvEHMyUNFiqGmGF764MGDYK/z/vvva6Bnb84HhtbYS1nujrI6g6xz5muj8kLgGxGY22lbblfeC5IadejQQX755Rc9icC8SQyBxedoQjnt7T8iokCBuXrIiI0ABvUOhpSiYc0MHq2PtZiSgGOoeazFcE1MM8BSUBgmiaAPw1St6wwch3E/pgdgHj2mLSDBDIIeDL/E8ND58+cHKxOCJcyLtwf1Gzir48JTVmfQyIvHm1C/OUoU5yrsB3us3wvqNzDfCxpUMVQVUzkwFxIZYLGfredXsn4jT2DwSGSTMCc8unTpovMuBg8ebLkNB3uzhRUwdwGtoZhziMQzPXr00OuY34C045j8bq/3zRHMp0CGuMgoqzOYd2gmE8K8E0fBIyp1tKRasw7+TGbLa1hhTiMqepycYJ8i+xx6jocMGaL3o4yoXG3nrxARBWrCnLBCIx2Wx0AQhrl6MGrUqGDzERGUog5FIjfMf//pp580eEMvGu5DkBiWOhYjSNCwide1TurjjrI6g8R15rxJwBx7R8FjZNZxCO4xVxT7FD20eB/IB4B9agaaeK+s3yiqMXgkcgMMx7E+gCMoRE+iCZULWj7RYoiJ7wjekC3OHEqEyfZIGIOeSeshqsjCZm+Svxk4odUYFV1YKg9nZXWF9dBR25ZiVJzmuoooOwI49Gya62mh59NdvcUYroOeYgSwSCSAhATWw5NwEoPXRYVLRERhhzoDwZwZjKGhccGCBcG2Mesq1ANoTDSzjWJUDeo3ZL3G8fjtt9+2PAYJ2xB8mXWHNfQoYoQOHlO3bl23ltUZjPjBny2U0zYwxHtGMGdCQyb+3MF6n2If4twhb968ejv2L94bziGwb4miEoNHIjdA5Yasohh2imEzaCG0XsQYgQ4qQgw7wTw+BH0IFvv376/3//DDD3o/MqJiSBFaKTdu3KiX1gs8W0OKcjwfsouaz+OOsoYXFlZGJYceVcxTxDwWZIBDRY7hpGj1xhAcR+8nrJBJD8+PExVksEMgiaHB1j2qc+bM0WVX7M0ZJSIi57DeMUZzdO7cWYMXBGM4lpu9X+ayGv/++68GOegxQ7ZyHJcxjBVBJIbJvvfeezqUFYEPGizRu4cM4/aCR2jdurWOlBk+fHiwqQgRLWt4oScVdQx6UtFAiukfTZs2lYEDB2oDMeprjAayXYIrvFA3Y/gs6ji85rRp07QM5tJa69at09e1F+gSRSbXfo1EfgwVG4IbRzBkpmLFisFuw8HbTAYAaPlD0ITK8+LFixq0oNIz159CKypShmN4KHrGkFQAQ3jMxevRooj1EBE4In06eiWRAhzPE5o+ffpoZYV5ke4qK5hzVkzoMW3YsGGoZcEJw5o1a7Q1FHMi8X5RYW/dulWDSbxvzN3A8FI8l7nUB1pWUQnbwokH5sRYw+NxcgCovLGvcIKA1lczGUPbtm31frw3vBbeGxFRIEKSFRxvUcfYg/nguN+2ARG3mfP0MHIDjZmYt4jjOfIBYP6idY/g0KFDtV5B3YZjPepULG9hBn0TJ07UugaNlngOPDe2QyOoIwiKEDQtXLjQrWVFcIfHoM4yoV5B/RKaH3/8UdedRMCLtRsBDbeYeoIRPBhlgzoISXCsR/fgtTE6xhreF8pg9pBavxdzTiMS3WEqBuoyJLTDfRi5Yw6BxQgg1G/mSB+iqBINKVej7NWIyO0GDRqkyWLMrHP0f9DDiWFRWG+LiIh8D4atolGyX79+ni6KV0GgivUef/vtN7f0qhKFBYNHIiIiIiIicorDVomIiIiIiMgpBo9ERERERETkFINHIiIiIiIicorBIxERERERETnF4JGIiIiIiIicYvBIRERERERETjF4JCIiIiIiIqcYPBIREREREZFTDB6JiIiIiIjIKQaPRERERERE5BSDRyIiIiIiIhJn/h8hSpZxoPlhEAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 900x320 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(9, 3.2), layout=\"constrained\")\n",
    "cases = [\n",
    "    (\"Trotterized evolution\", trotter_phases, [-exact_energy * total_time / 2]),\n",
    "    (\"Qubitized walk\", walk_phases, [0.25, 1.75]),\n",
    "]\n",
    "for axis, (title, probabilities, targets) in zip(axes, cases):\n",
    "    axis.bar(list(probabilities), list(probabilities.values()), width=0.8 * (2 / 2**n_phase))\n",
    "    for i, phase in enumerate(targets):\n",
    "        axis.axvline(phase, color=\"black\", linestyle=\"--\", linewidth=1,\n",
    "                     label=\"Ideal phase\" if i == 0 else None)\n",
    "    axis.set(title=title, xlabel=\"Phase (half-turns)\", xlim=(-0.05, 2), ylim=(0, 1.05))\n",
    "    axis.legend(fontsize=8)\n",
    "axes[0].set_ylabel(\"Sample probability\")\n",
    "axes[1].annotate(\"Same ground energy\", xy=(0.25, walk_phases[0.25]),\n",
    "                 xytext=(0.7, 0.9), arrowprops={\"arrowstyle\": \"->\"}, fontsize=9)\n",
    "axes[1].annotate(\"\", xy=(1.75, walk_phases[1.75]),\n",
    "                 xytext=(1.05, 0.87), arrowprops={\"arrowstyle\": \"->\"})\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "qpe-flow-29",
   "metadata": {},
   "source": [
    "## Go further: THC phase estimation\n",
    "\n",
    "- The {doc}`advanced THC example </examples/phase_estimation/thc_qpe>` builds a chemistry-oriented controlled LCU using alias-sampling PREPARE, QROM and orbital rotations.\n",
    "\n",
    "- It retains the alias workspace across walk steps and reflects its all-zero state; the SELECT record is uncomputed before each reflection.\n",
    "\n",
    "- It is an integration/type-checking example; the full circuit is too large for practical statevector simulation.\n",
    "\n",
    "- Continue with {doc}`canonical QPE </examples/phase_estimation/canonical_phase_estimation>`, {doc}`Trotterized QPE </examples/phase_estimation/zixy_phase_estimation>`, or {doc}`qubitized QPE </examples/phase_estimation/qubitized_phase_estimation>` for larger examples. The library also provides `iqpe` to reverse canonical QPE using inverse controlled powers, and `hadamard_test` as a one-control-qubit phase-sensitive primitive."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "guppyalgos (3.14.x)",
   "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.14.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}