calculate_config_entropy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f9114e15",
   "metadata": {},
   "source": [
    "# $ \\textit{Calculating Configurational entropies from GMM objects}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7edd888",
   "metadata": {},
   "source": [
    "### $\\textit{Import Libraries}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "44289e0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import pyemma\n",
    "from shapeGMM import gmm_shapes\n",
    "from shapeGMMTorch import torch_sgmm\n",
    "from shapeGMMTorch import scripts\n",
    "import MDAnalysis as md\n",
    "import torch\n",
    "import pickle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "ecf0b00b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDAnalysis.lib.distances import calc_dihedrals\n",
    "def calc_phi_psi(trajData):\n",
    "    phis = calc_dihedrals(trajData[:,0,:], trajData[:,1,:], trajData[:,2,:], trajData[:,3,:])*180/3.14\n",
    "    psis = calc_dihedrals(trajData[:,1,:], trajData[:,2,:], trajData[:,3,:], trajData[:,4,:])*180/3.14\n",
    "    return phis,psis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f863a4d6",
   "metadata": {},
   "source": [
    "### $\\textit{Metadynamics simulation (BF=10)}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8d3f553",
   "metadata": {},
   "source": [
    "### $\\textit{Path to files}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "2b891de5",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_path = \"/scratch/projects/hockygroup/ss12902/metad_sims_ala2_using_same_tpr_TREX/bf10/\"\n",
    "gro_file = \"run_ala2_metad_sigma0.3_height1.2_pace500_bf10.gro\"\n",
    "traj_file = \"run_ala2_metad_sigma0.3_height1.2_pace500_bf10_wrapped.trr\"\n",
    "colvar_file = \"run_ala2_metad_sigma0.3_height1.2_pace500_bf10.colvars.dat\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e38df593",
   "metadata": {},
   "source": [
    "### $\\textit{Load Metad fbias from colvar}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "a8d64d5f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x146e11121210>]"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAN0AAADFCAYAAADQQffsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAY2ElEQVR4nO3de1gTZ74H8G8CJIDkwkUSgiAgN0W8l4CCq5LlUtd6262ltFXX6rbFuhYrLVsr9rRPdW2P9axHcXfPKu3uKq37tPbsaulaVGgVaGW961K0WLAY6qUQ8CiivOcPy8BICATCJAO/z/PkeUjmzcwvQ76Z+7wSxhgDIUQwUnsXQMhgQ6EjRGAUOkIERqEjRGAUOkIERqEjRGAUOkIE5mzvAh7U2tqK2tpaKBQKSCQSe5dDBjDGGBobG6HT6SCVCrf8cbjQ1dbWIiAgwN5lkEGkpqYGw4YNE2x6Dhc6hUIB4P6MUCqVdq6GDGQmkwkBAQHcd04oDhe6tlVKpVJJoSOCEHozhnakECIwCh0hAqPQESIwh9ums6T+/+5g3H8c6HK4VukKo+m2xXHInaVovttq07rWzRqF/z50Adea7vR6HC8YwvHOZ193en1CoBr/qq7nnv90lAZJozRY/bdTvZ5WX3nIndHUfBcAkK4PxF/Lqju1GTF0CC5evWn1uGVOUty5Z/n/88spwdhxpIp7nhDmg88rr3XZ/v1lsdCHeFtdS3+RONr1dCaTCSqVCg0NDZ12pAS9vM9OVRGx27ciHlE6Fe81S9+1/kSrl2RQ+HPJt/YugUOhI4OCI63PWRW69evX46GHHoJCoYCvry/mzJmDiooKXpvbt28jIyMD3t7e8PDwwPz581FXV2fTogmxVqsDpc6q0BUVFSEjIwOlpaU4cOAAWlpakJSUhJs32zeYX3jhBfz973/Hnj17UFRUhNraWsybN8/mhRNijU/PGu1dAseqvZcFBQW853l5efD19UV5eTmmTp2KhoYG/OlPf8KuXbswY8YMAMDOnTsxcuRIlJaWIjY21naVE2IF0+279i6B06dtuoaGBgCAl5cXAKC8vBwtLS0wGAxcm8jISAQGBqKkpMTsOJqbm2EymXgPQgayXoeutbUVK1euxJQpUzB69GgAgNFohEwmg1qt5rXVaDQwGs0v3tevXw+VSsU96AoDMtD1OnQZGRk4c+YM8vPz+1RAdnY2GhoauEdNTU2fxkeIo+vVGSnLly/HP/7xDxQXF/OuQ9Jqtbhz5w7q6+t5S7u6ujpotVqz45LL5ZDL5b0pgxBRsmpJxxjD8uXL8dFHH+HgwYMIDg7mDZ84cSJcXFxQWFjIvVZRUYHq6mrExcXZpmJCRM6qJV1GRgZ27dqFjz/+GAqFgttOU6lUcHNzg0qlwpIlS5CZmQkvLy8olUo8//zziIuLoz2XhPzIqtDl5uYCAKZNm8Z7fefOnVi0aBEA4J133oFUKsX8+fPR3NyM5ORkbNu2zSbFEjIQWBW6npwb7erqiq1bt2Lr1q29LoqQgYzOvSREYBQ6QgRGoSNEYBQ6QgRGoSNEYBQ6QgRGoSNEYBQ6QgRGoSNEYBQ6QgRGoSNEYBQ6QgRGoSNEYBQ6QgRGoSNEYBQ6QgRGoSNEYBQ6QgRGoSNEYBQ6QgRGoSNEYBQ6QgRmdeiKi4sxa9Ys6HQ6SCQS7N27lzd80aJFkEgkvEdKSoqt6iVE9KwO3c2bNzF27FiL97VMSUnBlStXuMfu3bv7VCQhA4nVHYikpqYiNTXVYhu5XN5lhyGEDHb9sk13+PBh+Pr6IiIiAs8++yyuX7/eZVvqFJIMNjYPXUpKCt577z0UFhbit7/9LYqKipCamop79+6ZbU+dQpLBplf901ny2GOPcX9HR0djzJgxGDFiBA4fPozExMRO7bOzs5GZmck9N5lMFDwyoPX7IYOQkBD4+PjgwoULZofL5XIolUreg5CBrN9Dd/nyZVy/fh1+fn79PSlCRMHq1cumpibeUquqqgonTpyAl5cXvLy88Nprr2H+/PnQarW4ePEisrKyEBoaiuTkZJsWTohYWR26Y8eOYfr06dzztu2xhQsXIjc3F6dOncK7776L+vp66HQ6JCUl4fXXX6d+xQn5kdWhmzZtmsXOIT/99NM+FUTIQEfnXhIiMAodIQKj0BEiMAodIQKj0BEiMAodIQKj0BEiMAodIQKj0BEiMAodIQKj0BEiMAodIQKj0BEiMAodIQKj0BEiMAodIQKj0BEiMAodIQKj0BEiMAodIQKj0BEiMAodIQITVehcXURVLiFm2bwnVsYY1q5dCz8/P7i5ucFgMKCystImxf5sjM4m4yHEnmzeE+vGjRvxu9/9Dtu3b0dZWRmGDBmC5ORk3L59u8/FxoV493kchNibTXtiZYxh8+bNWLNmDWbPng0AeO+996DRaLB3715eN1ptmpub0dzczD231ClkgJe7teUS4nBsupFUVVUFo9EIg8HAvaZSqaDX61FSUmL2PdZ0CimR2LJaQuzDpqEzGo0AAI1Gw3tdo9Fwwx6UnZ2NhoYG7lFTU2PLkghxODbvidVacrmcevQhg4pNl3RarRYAUFdXx3u9rq6OG9YXtHZJBgKbhi44OBharRaFhYXcayaTCWVlZYiLi7PlpAgRLZv2xBoYGIiVK1fijTfeQFhYGIKDg/Hqq69Cp9Nhzpw5fS42ZKhHn8dBiL3ZtCfWvLw8ZGVl4ebNm1i2bBnq6+sRHx+PgoICuLq69rlYryGyPo+DEHuTMEvdqtqByWSCSqVCQ0MDlEplp+FBL++zQ1VkILi0YSbveXfftf5CJzMSIjAKHSECo9ARIjAKHSECo9ARIjAKHSECo9CRQcGRrlCh0JFBQadys3cJHAodGRQc6RwQ0YVutL9wZw4Q0h9EF7plU0fYuwRC+kR0oZs1xs/eJRARcpyVSxGGTuJIu6EI6QXRhY4QsaPQESIwCh0ZFBzoiAGFjhChUejIoOBIt/qg0JFB4Q9PTbR3CRwKHRkUtMq+3xjLVih0ZFBwpOO7ogzdihmh9i6BkF6zeejWrVsHiUTCe0RGRtp0GhFaOumZiFe/dCASFRWFzz77rH0izradjAOtKRCRcKSvTL+EztnZ2SYdhnTFxUmUa8WEAOinbbrKykrodDqEhIQgPT0d1dXVXbZtbm6GyWTiPbrjLHWk3y0iBo60dmTz0On1euTl5aGgoAC5ubmoqqpCQkICGhsbzba3pidWQgaCfu/LoL6+HsOHD8emTZuwZMmSTsPN9TkeEBBg8f7yh/79PRbnfdVvNZOBp2r9w50OG9irL4N+74lVrVYjPDyc171WR9QTKxHCoDpO19TUhIsXL8LPj674JgToh9C9+OKLKCoqwqVLl3D06FHMnTsXTk5OSEtLs9k0mENdfE+IdWy+enn58mWkpaXh+vXrGDp0KOLj41FaWoqhQ4faelKEiJLNQ5efn2/rURIyoIjyKLPEoc4vIMQ6ogwdIWJGoSNEYKIMHe29JGImytBZ8tw0uu06cWwDLnRZKba9do8QWxtwoXNkU0K97V0CcQAUOgFRj0MEEGnoFK4uFoe/9fMxvRrvihmh+PqN1F69l5CeEmXoJg33xKLJQV0OnxLqw/39x6cm9WicTlIJMpMiIHPmz5Jhnv3Xbe7q5IhOr6ndXRDk7d5v07SFEUOH2G3a+1ck2G3atiLK0EkkEqx7JKrL4Tq1G96YMxr/9dg4i+P5LHMq9/dTccPNtjm4alpvSjTLVyHH4ilB3PPYEC8ceXkGr81flujxRKz5WhyFv6f9fhRG6ZQYM0xlt+nbgihD12bfivguhz0ROxyzx/lbfL/3kPbr+KRdXG/14JKvL0b6KZGubw8UY4C/mr8kHe2vsrgU7624kL7vxPGQ3z9Vd8Gk3l/db27pbq2/PTO5z+OwJ1GHLkqn4r4IAPCCIbxTG4Wr+XO6w3w94Nnh/vYhZlaZfBU9v7h2/oRhVrXvaGyAmvfc2UnaKXguTn073zT7Yf6hFMNIDXw8rKv32BoDPl05FQ9Ht9906s250dixqGer8ADwhL7vS3GZs5T3fxcbUYcOAJ5OCAYApI7W4teGsE7D9cFeWDQ5CG/MGc295quQ46OMKQCA95fF4teJYWZ/vVckdh6fOV++koj/fHQsyn6TiKyUCPzhSevum795wTiMC1Dz3rcqif8DUvYbg8VxWHuvpj8+NZF3s55fTBzW7XtcXZwQoVXwrsJWuDpjRqTG4vtigr24v1XuLvjnC1OxsIvVeXNeTApHQpgP0mLa/0ezxup6/H5HI96fix+tmBGGaRG+GOVn/h4XHbf/1uw9A+D+l6Dtl1If4g39A6teB16YirKqG0iLCeSe//SdYrg4SdByr/MpaL4KV25az02zfPfpjquTbd/dYJ8h2Pvjj0CbB/fQOnVI1bu/jMHCHV/yhj+dEII/FH/T5XSDfdqX5E/HB3e6fYGrixNigrzw5aUbFutv8/rsKJR/+wMeju7+jgBb0sZj3rajmPljf/HhGgWmhPrg3ZJvO7WVOUlx514rAODfr6fAWSqBs5MUyx9oZ+4ujFvSxuP53cd7VL89iT50UqkE4x5YPeurMI0CYRoF7/mlDTMBAEEv7+vTuN1kTh2eWV48SSVAq5nTTCO1ik6vdXULkNPrknDnbisvxB337na05fHxWPe/Z/FvYyOqrt20WNuTcUF4Mi7IYps2KjcXfPHS9B7dp2T9vGis2nMSwP0fAmt0tV3uaEQfOqF9unIqPq+8itPfNeDjE7X2LoejU5k/tNExbF+8NB2VdU2YFmH+Kn6N0hW5T0xE9oenuw1dTxlG+loVHrW75WOwlojlRHjRb9MJLUKrwNMJId3eZfqvT+vx6KRhvO2QBw2RW/dL3lHbNliERoFFk4PwuD6w2/cM83TH9EhfQe+M9T8LH+q2Tcft0YSwoYjUKjBvguU9zx05SyXwkDuL5s7fg3JJZ4svnblVvI6mhPpgSqgP3tx/vtOwdbNGwWhqRmQfOkJ5c1400mOHI9pfxdves5ey3yTim6s38fzu47jW1Nxt+47/g0nD27clZc5SfPLrhG7/R9PCffGX0mrInKU4lZMEiQQo/voaN/zLVxKx+bNK7CqrhtyGh31sYVCGzhaeigvC7ZZ7SAiz/oZLi6YE93n6Lk5Sm2/L9oVG6QqNjTpe7MmPYuJIX+xeGoswjYfZ1VdfhSteeXgkgrzdkRzVf/1q9IZj/QSIiMxZiuUzwjodYxvsOh6a6U8SiQRxI7x5xxp1an7oh8idsWzqCAz3tt9pa+ZQ6BzY4h+XiMlRGt4qZFerk2/OjUZsiBd+Em790nfO+PbjXisSe9/pZspoLfQdjssJKUqnwtu/GItdS/V2mX5P0eqlA3s5NRKJI30xIdATri5OWJoQjLutrMszSR7XB+JxfSA2/bMCRV9ftWpaE4e3B8VP5YYPn5uMeduO9qruYJ8hKKvq2fE+W/t5Dw7y21u/Lem2bt2KoKAguLq6Qq/X48svv+z+TYTHxUmKySN8uG2WV2aOQs6srk/0bmOLHUXjhqmREGb+eF53slIi8ZPwodj6+IQetY8dcf/kBEfYISSEflnSvf/++8jMzMT27duh1+uxefNmJCcno6KiAr6+vv0xSdLB4ilB+Fv5ZfxsTO/7j5BKJfjzEj2OXryGZ//yL7xuxbaa1xAZ3v1lTI/bZ0wfAV+FvFerxWLUL0u6TZs2YenSpVi8eDFGjRqF7du3w93dHTt27OjUtjedQvZW24m6S+L7vvfQkandZfjipenIfnhkn8c1eYQPTqz9KR6x8bmOHTv2lDs74YnY4QjwcuzrCG3F5qG7c+cOysvLYTC0n6ArlUphMBhQUlLSqb2QnUJufXwCTuYkCbqrfUl8MBRyZ6tO8LWFnqxitl3blxhpee2jPw6mx4f5YEKgGuk9OKg/0Ni8U8ja2lr4+/vj6NGjiIuL417PyspCUVERysrKeO170ymk2NxrZQ65vdLaynDycj1G6ZSQO/f+7BixGrCdQnZnMHQK6YiBA+5vt40P9LR3GYOOzVcvfXx84OTkhLq6Ot7rdXV10God68wAQuzB5qGTyWSYOHEiCgsLuddaW1tRWFjIW90kZLDql9XLzMxMLFy4EJMmTUJMTAw2b96MmzdvYvHixf0xOUJEpV9Ct2DBAly9ehVr166F0WjEuHHjUFBQAI3G8mX9ANC2X6c/Dx0QArR/x2y8L7FbNt972VeXL1/u18MGhDyopqYGw4YJd/qYw4WutbUVtbW1UCgUnY4PtR1OqKmpGTCHE2yF5o15luYLYwyNjY3Q6XSQSoU799/uhwweJJVKu/3VUSqV9MXqAs0b87qaLyqV8DeupUt7CBEYhY4QgYkqdHK5HDk5OQP+DJbeoHljniPOF4fbkULIQCeqJR0hAwGFjhCBUegIERiFjhCBUegIEZioQjeQ7jC2bt06SCQS3iMysr3jxtu3byMjIwPe3t7w8PDA/PnzO12jWF1djZkzZ8Ld3R2+vr5YvXo17t69y2tz+PBhTJgwAXK5HKGhocjLy+tUi73na3FxMWbNmgWdTgeJRIK9e/fyhjPGsHbtWvj5+cHNzQ0GgwGVlZW8Njdu3EB6ejqUSiXUajWWLFmCpqYmXptTp04hISEBrq6uCAgIwMaNGzvVsmfPHkRGRsLV1RXR0dHYv3+/1bV0i4lEfn4+k8lkbMeOHezs2bNs6dKlTK1Ws7q6OnuX1is5OTksKiqKXblyhXtcvXqVG/7MM8+wgIAAVlhYyI4dO8ZiY2PZ5MmTueF3795lo0ePZgaDgR0/fpzt37+f+fj4sOzsbK7NN998w9zd3VlmZiY7d+4c27JlC3NycmIFBQVcG0eYr/v372evvPIK+/DDDxkA9tFHH/GGb9iwgalUKrZ371528uRJ9sgjj7Dg4GB269Ytrk1KSgobO3YsKy0tZZ9//jkLDQ1laWlp3PCGhgam0WhYeno6O3PmDNu9ezdzc3Njv//977k2R44cYU5OTmzjxo3s3LlzbM2aNczFxYWdPn3aqlq6I5rQxcTEsIyMDO75vXv3mE6nY+vXr7djVb2Xk5PDxo4da3ZYfX09c3FxYXv27OFeO3/+PAPASkpKGGP3v6hSqZQZjUauTW5uLlMqlay5uZkxxlhWVhaLiorijXvBggUsOTmZe+5o8/XB0LW2tjKtVsveeust7rX6+noml8vZ7t27GWOMnTt3jgFgX331Fdfmk08+YRKJhH333XeMMca2bdvGPD09uXnDGGMvvfQSi4iI4J4/+uijbObMmbx69Ho9+9WvftXjWnpCFKuX1t5hTCwqKyuh0+kQEhKC9PR0VFdXAwDKy8vR0tLC+7yRkZEIDAzkPm9JSQmio6N51ygmJyfDZDLh7NmzXJuO42hr0zYOMczXqqoqGI1GXo0qlQp6vZ43L9RqNSZNau/73GAwQCqVcjfCKikpwdSpUyGTtfcz33Yv1h9++IFrY2l+9aSWnhBF6K5du4Z79+51ughWo9HAaDTaqaq+0ev1yMvLQ0FBAXJzc1FVVYWEhAQ0NjbCaDRCJpNBrVbz3tPx8xqNRrPzo22YpTYmkwm3bt0SxXxtq8NSjUajsdNNjJ2dneHl5WWT+dVxeHe19ITDXdozWKSmpnJ/jxkzBnq9HsOHD8cHH3wANzfzvaqSgUEUS7rBcIcxtVqN8PBwXLhwAVqtFnfu3EF9fT2vTcfPq9Vqzc6PtmGW2iiVSri5uYlivrbVYalGrVaL77//njf87t27uHHjhk3mV8fh3dXSE6II3WC4w1hTUxMuXrwIPz8/TJw4ES4uLrzPW1FRgerqau7zxsXF4fTp07wv24EDB6BUKjFq1CiuTcdxtLVpG4cY5mtwcDC0Wi2vRpPJhLKyMt68qK+vR3l5Odfm4MGDaG1thV6v59oUFxejpaWFa3PgwAFERETA09OTa2NpfvWklh7p8S4XO8vPz2dyuZzl5eWxc+fOsWXLljG1Ws3beycmq1atYocPH2ZVVVXsyJEjzGAwMB8fH/b9998zxu4fMggMDGQHDx5kx44dY3FxcSwuLo57f9shg6SkJHbixAlWUFDAhg4davaQwerVq9n58+fZ1q1bzR4ysPd8bWxsZMePH2fHjx9nANimTZvY8ePH2bfffssYu7+bXq1Ws48//pidOnWKzZ492+whg/Hjx7OysjL2xRdfsLCwMN4hg/r6eqbRaNiTTz7Jzpw5w/Lz85m7u3unQwbOzs7s7bffZufPn2c5OTlmDxl0V0t3RBM6xhjbsmULCwwMZDKZjMXExLDS0lJ7l9RrCxYsYH5+fkwmkzF/f3+2YMECduHCBW74rVu32HPPPcc8PT2Zu7s7mzt3Lrty5QpvHJcuXWKpqanMzc2N+fj4sFWrVrGWlhZem0OHDrFx48YxmUzGQkJC2M6dOzvVYu/5eujQIQag02PhwoWMsfu76l999VWm0WiYXC5niYmJrKKigjeO69evs7S0NObh4cGUSiVbvHgxa2xs5LU5efIki4+PZ3K5nPn7+7MNGzZ0quWDDz5g4eHhTCaTsaioKLZv3z7e8J7U0h26no4QgYlim46QgYRCR4jAKHSECIxCR4jAKHSECIxCR4jAKHSECIxCR4jAKHSECIxCR4jAKHSECOz/AWXePpQQZ7xsAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 200x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "metad_fbias = np.loadtxt(data_path+\"colvar_reweight\", usecols=3)\n",
    "metad_fbias -= metad_fbias.min()\n",
    "\n",
    "plt.figure(figsize=(2,2))\n",
    "plt.plot(metad_fbias)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "27082f42",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x146e10f6afe0>]"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQQAAADFCAYAAACyy4yvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAhd0lEQVR4nO3de1hTV7438C8RSdA2QUslxKZK31LvBYs1hmPr9JgxVp6OzOs5BYep1qHSzkhbB1tFq+B0fA6+WFulOqLttMzF+2lLO6hRClq8xIBcLHIrKgqoARSTAMotWe8fDku3BEgoCuLv8zz7waz923uvbMI3O3vvLF0YYwyEEAJA1NsdIIT0HRQIhBCOAoEQwlEgEEI4CgRCCEeBQAjhKBAIIZxrb3fgQWOz2XD58mU8+uijcHFx6e3ukH6MMYa6ujooFAqIRPfnvZsCwUmXL1+GUqns7W6Qh0hFRQWeeOKJ+7ItCgQnPfroowBu/ZKkUmkv94b0ZxaLBUqlkr/m7gcKBCe1fUyQSqUUCOS+uJ8fTemkIiGEo0AghHDdCoTNmzdj5MiRkEgkUKlUyMzM7LR+7969GD16NCQSCSZMmID9+/cL5jPGEBMTA29vb7i7u0Oj0aC0tFRQU1tbi7CwMEilUnh4eCA8PBz19fV8fmNjI15//XVMmDABrq6uCA4ObtePI0eOwMXFpd1kNBq7sxsI6XecDoTdu3cjKioKsbGxyMnJgZ+fH7RaLaqrq+3WnzhxAnPnzkV4eDhyc3MRHByM4OBgnDlzhtfEx8cjISEBiYmJMBgMGDx4MLRaLRobG3lNWFgYCgoKkJqaipSUFGRkZCAiIoLPt1qtcHd3xzvvvAONRtPpcygpKcGVK1f4NGzYMGd3Q59ytroOScfL0GK1OVhfjxvNrfzxdsNFbDlyrof7VI9F23NQbLQAAC6ZbmJNSiEqam+gudWGY6VX0dhixdHSGpytru9ibe2duWRGenFVj/a5M1bb7VECztXU4/Oj59HYYu2wvsrSiJvN7eczxtDq4O+pVzAnTZ48mS1atIg/tlqtTKFQsLi4OLv1r776KgsKChK0qVQq9uabbzLGGLPZbEwul7N169bx+SaTiYnFYrZz507GGGOFhYUMAMvKyuI1Bw4cYC4uLuzSpUvttjl//nw2e/bsdu2HDx9mANj169cdfr6NjY3MbDbzqaKiggFgZrO5y2ULL5tZXvl11txqZZXXbzi8zTt9fvQ8Szxy1u68zYdLWcL3P7ERy1LYiGUpbOsP9uvulFV2jY1YlsLU//M9b2tbvvxaA2OMsbrGFpZeXMU+OljM1qQUMMYYs1ptTH/uKqtrbHGo31P+53s2YlkKG7PqAGOMMc36I2zEshQ2LT6drfwmn41YlsJ+se4w33abzLJrLL24qsv1ty1XWmXhbbX1TaypxepQ/+7WarW1e242m42VGC2s+IqFjYvRsfUHiwXbbnt8t/JrDWzEshTm/6eDgnW1Pe/xMTpmamjusk9ms9nh11pPceoIobm5GdnZ2YJ3YJFIBI1GA71eb3cZvV7f7h1bq9Xy+rKyMhiNRkGNTCaDSqXiNXq9Hh4eHpg0aRKv0Wg0EIlEMBgMzjwFAIC/vz+8vb3xy1/+EsePH++0Ni4uDjKZjE+O3oNgtTG8vPEoZm8+jlc+PYb/WJuOo6U1+OhgCY6fverQOm42W/HnlELEHShGbUOzYF5qYRXidSVYn/oTbztdYe5ynfvyrwAALptvHX3d+S5349/vaOFJWVjwZRY+TT+Lz46WoaL2Bv6uv4DQbSfx34nC33PZ1QYk/nAO5pstvK3EWIcr/15/2zpL/30UcOHaDfzj5EW+7N3+O1GPBV9mobqusd08AGhqteK705f544vXbqDYaMGkNd9j4p9Tofn4BwDAtfomsE7G/tGduYKUH2+v578ST2B87EEYzbe3uy//CmZ8kgHthgzUN7UiIf2sYB3Z5dftrvto6a3f7/Ubt/fJ6Uozf951Ta3w+/AQ8iu7/n3db04FwtWrV2G1WuHl5SVo9/Ly6vBzuNFo7LS+7WdXNXcf1ru6umLo0KFOff739vZGYmIivvrqK3z11VdQKpX4xS9+gZycnA6XWb58OcxmM58qKiq63M71hmY0t94+LCw21gEAXvtrJjYdPouwzw240dwKw/lrqKi9gX+dviw4JG3Taru9jqZW4eHn0dIau9t+daseI6P3oaL2Rpf9bGq1YmyMrl27oaxW8LjZasPXuZcAAEVXLIJ5L310BGsPFMPvT4f4IXJQwtEut21P9sXb271W32y3ZuP3pXhnZ66gLeLv2bha3wQAKK+9ge9OX0bAmu+x+rsCMMbaBUNjixVv/TMHkTtykZBWioraG8gtNwEA/nrsPP/otTOzvFvPw54bTa3t2l7ZdKzH1t9THqr7EEaNGoVRo0bxx4GBgTh37hw++eQT/OMf/7C7jFgshlgsdngbJcY6aDdkwE/p0Wnd619mIfOOPzzTjWa8ph7p8HY60rbOGZ9koOjPMzutXf51Pu7OoTvfNe9UVtP+3fxuY2J0OBH9n2i1E25dYYxhzpbbRx8vbzwK3eIXMFouvNfjf7Mr2y1ruiEMj7X7iwAAf9NfRNGVOmReqMXJ5dMhl0kAAE13hPXHqT9h6w+3z598drQMWReuI3nRfzj9HDpTabrZo+u7V5wKBE9PTwwYMABVVcKTOVVVVZDL5XaXkcvlnda3/ayqqoK3t7egxt/fn9fcfdKytbUVtbW1HW7XUZMnT8axYz2T1EZzI7QbMgAApytMndZm3vUuvOH7Uly4dgMZP9Vgy2+fw9PDhHen6c9dg2TgAOzKqsBLox5HQ1P7E1aGsmv83zdbrBgZvQ/B/gosnTkav0vKwij5o/g27/Yf/Nc5lwTLJ/5wDt/kCtsA4P/+5QTq7niHM99sgcx9oN3ntWZfYSfPumN5dvbXzA1HsStiClxFLpg0cihuNLeiuq6py3Xd2dfMC7f285S4NHw4exz251/BtGeER5sNd538y6swIbWwCsfPXoOjTpy7inUHS7AyaAxWfJMvmLf58FmsO1ji8Lp6k1MfGdzc3BAQEIC0tDTeZrPZkJaWBrVabXcZtVotqAeA1NRUXu/j4wO5XC6osVgsMBgMvEatVsNkMiE7O5vXpKenw2azQaVSOfMU2snLyxME0c+xZG9et5e91tCMvx4rQ2l1PTQfZ7Q7NI/acxp/2J6DjJ9q8Kd/FeKrnPbvlFftHGYn511G4Np0FBvrBGFgj70wACA4PwAAfn86hKX/e9rux5L9+d27hPvrv5yw2x667ST+K1GP5lYbxsUedGhddY3tD88BIObbApw8X4v/pyvuch0L/37KbvuH/7odeMfPXkOr1YYV3+TjN58ZkFtuEhzltNV3FgadnefoDU5/ZIiKisL8+fMxadIkTJ48GRs2bEBDQwMWLFgAAJg3bx6GDx+OuLg4AMC7776LadOmYf369QgKCsKuXbtw6tQpbNu2DcCt2zIXL16MNWvWwNfXFz4+Pli1ahUUCgW/l2DMmDGYOXMmFi5ciMTERLS0tCAyMhKhoaFQKBS8b4WFhWhubkZtbS3q6uqQl5cHAPxIY8OGDfDx8cG4cePQ2NiIzz//HOnp6Th06FB395/AeQcOqx318sajyF89o8fW19P2nKrEDz/ZP49xL7RYbbD3t1NaXQ9LBwFwL3xxvEzweO2BYuwwdHyu4e76u/3pX4VY/atxPdK3nuB0IISEhKCmpgYxMTEwGo3w9/eHTqfjJwXLy8sFX9UMDAzEjh07sHLlSqxYsQK+vr5ITk7G+PHjec3SpUvR0NCAiIgImEwmTJ06FTqdDhKJhNds374dkZGRmD59OkQiEebMmYOEhARB32bNmoWLFy/yxxMnTgRwO4Wbm5uxZMkSXLp0CYMGDcKzzz6L77//Hi+99JKzu8EueycGf45LffxzZ5Wl68N3AKi83vUJzu5ae6Drd/t76fNjnf/BdyXpxIU+FQgurK8ds/RxFosFMpkMZrNZ8OWm5lYbnll5oBd71r/lrvolJv45tbe7cU9cWBtkt72j19q9RN9l6CH5l/reNeX+pL+GQV9DgUAI4SgQegx98iIPPgqEHtLD5xMJ6RUUCIQQjgKBEMJRIBDSy85W1/V2FzgKBEJ6WeSO3K6L7hMKhB5y95eVCHFU29fj+wIKhB5yP+/rJ+ReoUAghHAUCIQQjgKBEMJRIBBCOAoEQghHgdBD7t9/x0nIvUOBQAjhKBAIIRwFAiGEo0AghHAUCIQQjgKhh9CASaQ/oEAghHAUCIQQjgKBEMJRIBBCOAqEHkK3LpP+gAKBEMJRIBBCOAqEHkL3IZD+gAKBEMJRIBBCOAoEQghHgUAI4SgQCCFctwJh8+bNGDlyJCQSCVQqFTIzMzut37t3L0aPHg2JRIIJEyZg//79gvmMMcTExMDb2xvu7u7QaDQoLS0V1NTW1iIsLAxSqRQeHh4IDw9HfX09n9/Y2IjXX38dEyZMgKurK4KDg+325ciRI3juuecgFovx9NNPIykpqTu7gJB+yelA2L17N6KiohAbG4ucnBz4+flBq9Wiurrabv2JEycwd+5chIeHIzc3F8HBwQgODsaZM2d4TXx8PBISEpCYmAiDwYDBgwdDq9WisbGR14SFhaGgoACpqalISUlBRkYGIiIi+Hyr1Qp3d3e888470Gg0dvtSVlaGoKAgvPTSS8jLy8PixYvxxhtv4ODBg87uBkL6JRfGmFOX0FUqFZ5//nls2rQJAGCz2aBUKvH2228jOjq6XX1ISAgaGhqQkpLC26ZMmQJ/f38kJiaCMQaFQoElS5bgvffeAwCYzWZ4eXkhKSkJoaGhKCoqwtixY5GVlYVJkyYBAHQ6HWbNmoXKykooFArBNl9//XWYTCYkJycL2pctW4Z9+/YJwig0NBQmkwk6nc6h52+xWCCTyWA2myGVSm8/z616GOg/fCXddGFtULu2jl5r95JTRwjNzc3Izs4WvAOLRCJoNBro9Xq7y+j1+nbv2FqtlteXlZXBaDQKamQyGVQqFa/R6/Xw8PDgYQAAGo0GIpEIBoPB4f531Rd7mpqaYLFYBJM9dGMS6Q+cCoSrV6/CarXCy8tL0O7l5QWj0Wh3GaPR2Gl928+uaoYNGyaY7+rqiqFDh3a4XWf6YrFYcPPmTbvLxMXFQSaT8UmpVDq8PUIeNHSVoQvLly+H2WzmU0VFhf1COkQg/YBTgeDp6YkBAwagqqpK0F5VVQW5XG53Gblc3ml928+uau4+adna2ora2toOt+tMX6RSKdzd3e0uIxaLIZVKBRMh/ZVTgeDm5oaAgACkpaXxNpvNhrS0NKjVarvLqNVqQT0ApKam8nofHx/I5XJBjcVigcFg4DVqtRomkwnZ2dm8Jj09HTabDSqVyuH+d9UXQh52rs4uEBUVhfnz52PSpEmYPHkyNmzYgIaGBixYsAAAMG/ePAwfPhxxcXEAgHfffRfTpk3D+vXrERQUhF27duHUqVPYtm0bAMDFxQWLFy/GmjVr4OvrCx8fH6xatQoKhYLfSzBmzBjMnDkTCxcuRGJiIlpaWhAZGYnQ0FDBFYbCwkI0NzejtrYWdXV1yMvLAwD4+/sDAN566y1s2rQJS5cuxe9+9zukp6djz5492LdvX3f33200QgrpB5wOhJCQENTU1CAmJgZGoxH+/v7Q6XT8ZF15eTlEotsHHoGBgdixYwdWrlyJFStWwNfXF8nJyRg/fjyvWbp0KRoaGhAREQGTyYSpU6dCp9NBIpHwmu3btyMyMhLTp0+HSCTCnDlzkJCQIOjbrFmzcPHiRf544sSJAG7d+ATcOhrZt28f/vjHP2Ljxo144okn8Pnnn0Or1Tq7G9qhPCD9gdP3ITzs6D4Eci88kPchEEL6NwoEQghHgdBDXOgkAukHKBAIIRwFQg+hU7OkP6BA6CGUB6Q/oEAghHAUCIQQjgKhh9BFBtIfUCAQQjgKhB5CJxVJf0CBQAjhKBAIIRwFQg+hk4qkP6BA6CF0DoH0BxQIhBCOAoEQwlEgEEI4CgRCCEeB0EPoKgPpDygQCCEcBQIhhKNAIIRwFAiEEI4CgRDCUSD0ELp1mfQHFAiEEI4CgRDCUSD0FPrMQPoBCgRCCEeB0FPo3mXSD1AgEEI4CgRCCEeBQAjhKBB6Cl1lIP1AtwJh8+bNGDlyJCQSCVQqFTIzMzut37t3L0aPHg2JRIIJEyZg//79gvmMMcTExMDb2xvu7u7QaDQoLS0V1NTW1iIsLAxSqRQeHh4IDw9HfX29oObHH3/ECy+8AIlEAqVSifj4eMH8pKQkuLi4CCaJRNKdXUBIv+R0IOzevRtRUVGIjY1FTk4O/Pz8oNVqUV1dbbf+xIkTmDt3LsLDw5Gbm4vg4GAEBwfjzJkzvCY+Ph4JCQlITEyEwWDA4MGDodVq0djYyGvCwsJQUFCA1NRUpKSkICMjAxEREXy+xWLBjBkzMGLECGRnZ2PdunVYvXo1tm3bJuiPVCrFlStX+HTx4kVnd4F9dJWB9AMujDGnDnZVKhWef/55bNq0CQBgs9mgVCrx9ttvIzo6ul19SEgIGhoakJKSwtumTJkCf39/JCYmgjEGhUKBJUuW4L333gMAmM1meHl5ISkpCaGhoSgqKsLYsWORlZWFSZMmAQB0Oh1mzZqFyspKKBQKbNmyBR988AGMRiPc3NwAANHR0UhOTkZxcTGAW0cIixcvhslkcn5P/ZvFYoFMJoPZbIZUKuXtr27VI7OsttvrJQ+3C2uD2rV19Fq7l5w6QmhubkZ2djY0Gs3tFYhE0Gg00Ov1dpfR6/WCegDQarW8vqysDEajUVAjk8mgUql4jV6vh4eHBw8DANBoNBCJRDAYDLzmxRdf5GHQtp2SkhJcv36dt9XX12PEiBFQKpWYPXs2CgoKOn3OTU1NsFgsgomQ/sqpQLh69SqsViu8vLwE7V5eXjAajXaXMRqNnda3/eyqZtiwYYL5rq6uGDp0qKDG3jru3MaoUaPwxRdf4Ntvv8U///lP2Gw2BAYGorKyssPnHBcXB5lMxielUtlhLSEPuofqKoNarca8efPg7++PadOm4euvv8bjjz+OrVu3drjM8uXLYTab+VRRUXEfe0zI/eVUIHh6emLAgAGoqqoStFdVVUEul9tdRi6Xd1rf9rOrmrtPWra2tqK2tlZQY28dd27jbgMHDsTEiRNx9uxZ+08YgFgshlQqFUz20DlF0h84FQhubm4ICAhAWloab7PZbEhLS4Narba7jFqtFtQDQGpqKq/38fGBXC4X1FgsFhgMBl6jVqthMpmQnZ3Na9LT02Gz2aBSqXhNRkYGWlpaBNsZNWoUhgwZYrdvVqsV+fn58Pb2dmY32OVCiUD6Aac/MkRFReGzzz7D3/72NxQVFeH3v/89GhoasGDBAgDAvHnzsHz5cl7/7rvvQqfTYf369SguLsbq1atx6tQpREZGAgBcXFywePFirFmzBt999x3y8/Mxb948KBQKBAcHAwDGjBmDmTNnYuHChcjMzMTx48cRGRmJ0NBQKBQKAMBvfvMbuLm5ITw8HAUFBdi9ezc2btyIqKgo3pcPP/wQhw4dwvnz55GTk4Pf/va3uHjxIt54441u70BC+hNXZxcICQlBTU0NYmJiYDQa4e/vD51Ox0/glZeXQyS6nTOBgYHYsWMHVq5ciRUrVsDX1xfJyckYP348r1m6dCkaGhoQEREBk8mEqVOnQqfTCW4a2r59OyIjIzF9+nSIRCLMmTMHCQkJfL5MJsOhQ4ewaNEiBAQEwNPTEzExMYJ7Fa5fv46FCxfCaDRiyJAhCAgIwIkTJzB27Fhnd0M7zl28JaRvcvo+hIddR9eGQ7bqYaD7EEg3PZD3IZCOSQYO6O0uEPKzUSD0kMk+Q3u7C4T8bBQIhBCOAoEQwlEgEEI4CgRCCEeBQAjhKBAIIRwFAiGEo0AghHAUCIQQjgKBEMJRIBBCOAoEQghHgUAI4SgQCCEcBQIhhKNA6CE0yCrpDygQCCEcBUIPoZEpSX9AgUAI4SgQCCEcBQIhhKNAIIRwFAiEEI4CgRDCUSAQQjgKBEIIR4FACOEoEAghHAUCIYSjQCCEcBQIhBCOAoEQwlEg9BAaIIX0BxQIhBCuW4GwefNmjBw5EhKJBCqVCpmZmZ3W7927F6NHj4ZEIsGECROwf/9+wXzGGGJiYuDt7Q13d3doNBqUlpYKampraxEWFgapVAoPDw+Eh4ejvr5eUPPjjz/ihRdegEQigVKpRHx8vNN96a76xtYeWQ8hvcnpQNi9ezeioqIQGxuLnJwc+Pn5QavVorq62m79iRMnMHfuXISHhyM3NxfBwcEIDg7GmTNneE18fDwSEhKQmJgIg8GAwYMHQ6vVorGxkdeEhYWhoKAAqampSElJQUZGBiIiIvh8i8WCGTNmYMSIEcjOzsa6deuwevVqbNu2zam+dJdk4ICfvQ5CepsLY84N/qVSqfD8889j06ZNAACbzQalUom3334b0dHR7epDQkLQ0NCAlJQU3jZlyhT4+/sjMTERjDEoFAosWbIE7733HgDAbDbDy8sLSUlJCA0NRVFREcaOHYusrCxMmjQJAKDT6TBr1ixUVlZCoVBgy5Yt+OCDD2A0GuHm5gYAiI6ORnJyMoqLix3qiz1NTU1oamrijy0WC5RKJcxmM6RSKW//JPUnbEwrtbcKQrp0YW1QuzaLxQKZTNbutXYvOXWE0NzcjOzsbGg0mtsrEImg0Wig1+vtLqPX6wX1AKDVanl9WVkZjEajoEYmk0GlUvEavV4PDw8PHgYAoNFoIBKJYDAYeM2LL77Iw6BtOyUlJbh+/bpDfbEnLi4OMpmMT0ql0m5d4P95rMN1EPKgcCoQrl69CqvVCi8vL0G7l5cXjEaj3WWMRmOn9W0/u6oZNmyYYL6rqyuGDh0qqLG3jju30VVf7Fm+fDnMZjOfKioq7NapnnoMQwe72Z1HSGeOvPeL3u4C59rbHejrxGIxxGKxQ7U5q355j3tDyL3l1BGCp6cnBgwYgKqqKkF7VVUV5HK53WXkcnmn9W0/u6q5+6Rla2sramtrBTX21nHnNrrqCyEPO6cCwc3NDQEBAUhLS+NtNpsNaWlpUKvVdpdRq9WCegBITU3l9T4+PpDL5YIai8UCg8HAa9RqNUwmE7Kzs3lNeno6bDYbVCoVr8nIyEBLS4tgO6NGjcKQIUMc6gshDz3mpF27djGxWMySkpJYYWEhi4iIYB4eHsxoNDLGGHvttddYdHQ0rz9+/DhzdXVlH330ESsqKmKxsbFs4MCBLD8/n9esXbuWeXh4sG+//Zb9+OOPbPbs2czHx4fdvHmT18ycOZNNnDiRGQwGduzYMebr68vmzp3L55tMJubl5cVee+01dubMGbZr1y42aNAgtnXrVqf60hWz2cwAMLPZ7OyuI8QpvfFac/ocQkhICGpqahATEwOj0Qh/f3/odDp+sq68vBwi0e0Dj8DAQOzYsQMrV67EihUr4Ovri+TkZIwfP57XLF26FA0NDYiIiIDJZMLUqVOh0+kgkUh4zfbt2xEZGYnp06dDJBJhzpw5SEhI4PNlMhkOHTqERYsWISAgAJ6enoiJiRHcq+BIXxwIUAC3jmIIuZfaXmPsPv63YE7fh/Cwq6ys7PDSIyH3QkVFBZ544on7si0KBCfZbDZcvnwZjz76KFzu+kZT201LFRUV9+1GkgcB7ZeOdbZvGGOoq6uDQqEQHHXfS3TZ0UkikajLtJZKpfTCt4P2S8c62jcymey+9oO+7UgI4SgQCCEcBUIPEovFiI2NdfjOxocF7ZeO9bV9QycVCSEcHSEQQjgKBEIIR4FACOEoEAghHAUCIYSjQOghzo5E3ZetXr0aLi4ugmn06NF8fmNjIxYtWoTHHnsMjzzyCObMmdNunIny8nIEBQVh0KBBGDZsGN5//320tgpHpj5y5Aiee+45iMViPP3000hKSmrXl97erxkZGXjllVegUCjg4uKC5ORkwXzWh0YMd6QvXbpv36vsx3bt2sXc3NzYF198wQoKCtjChQuZh4cHq6qq6u2udUtsbCwbN24cu3LlCp9qamr4/LfeeosplUqWlpbGTp06xaZMmcICAwP5/NbWVjZ+/Him0WhYbm4u279/P/P09GTLly/nNefPn2eDBg1iUVFRrLCwkH366adswIABTKfT8Zq+sF/379/PPvjgA/b1118zAOybb74RzF+7di2TyWQsOTmZnT59mv3qV7+y+9V9Pz8/dvLkSXb06FH29NNPC766bzabmZeXFwsLC2NnzpxhO3fuZO7u7u2+uj9gwAAWHx/PCgsL2cqVK+0OI9BVX7pCgdADJk+ezBYtWsQfW61WplAoWFxcXC/2qvtiY2OZn5+f3Xkmk4kNHDiQ7d27l7cVFRUxAEyv1zPGbv0RiUQiPkYGY4xt2bKFSaVS1tTUxBhjbOnSpWzcuHGCdYeEhDCtVssf97X9encg2Gw2JpfL2bp163ibyWRiYrGY7dy5kzHGWGFhIQPAsrKyeM2BAweYi4sLu3TpEmOMsb/85S9syJAhfN8wxtiyZcvYqFGj+ONXX32VBQUFCfqjUqnYm2++6XBfHEEfGX6m7oxE/SAoLS2FQqHAU089hbCwMJSXlwMAsrOz0dLSIni+o0ePxpNPPikYJXvChAmCAW21Wi0sFgsKCgp4TWcjYD8I+7UvjRjuSF8cQYHwM3VnJOq+TqVSISkpCTqdDlu2bEFZWRleeOEF1NXV8f/3wsPDQ7DM3aNkd3cEbIvFgps3bz4Q+7UvjRjuSF8cQV9/Ju28/PLL/N/PPvssVCoVRowYgT179sDd3b0Xe0buNTpC+Jm6MxL1g8bDwwPPPPMMzp49C7lcjubmZphMJkHN3aNkd3cEbKlUCnd39wdiv/alEcMd6YsjKBB+pu6MRP2gqa+vx7lz5+Dt7Y2AgAAMHDhQ8HxLSkpQXl4uGCU7Pz9f8IeQmpoKqVSKsWPH8prORsB+EPZrXxox3JG+OMTh04+kQ12NRP2gWbJkCTty5AgrKytjx48fZxqNhnl6erLq6mrG2K3Ljk8++SRLT09np06dYmq1mqnVar5822XHGTNmsLy8PKbT6djjjz9u97Lj+++/z4qKitjmzZvtXnbs7f1aV1fHcnNzWW5uLgPAPv74Y5abm8suXrzIGOtbI4Y70peuUCD0kE8//ZQ9+eSTzM3NjU2ePJmdPHmyt7vUbSEhIczb25u5ubmx4cOHs5CQEHb27Fk+/+bNm+wPf/gDGzJkCBs0aBD79a9/za5cuSJYx4ULF9jLL7/M3N3dmaenJ1uyZAlraWkR1Bw+fJj5+/szNzc39tRTT7Evv/yyXV96e78ePnyYAWg3zZ8/nzF263LfqlWrmJeXFxOLxWz69OmspKREsI5r166xuXPnskceeYRJpVK2YMECVldXJ6g5ffo0mzp1KhOLxWz48OFs7dq17fqyZ88e9swzzzA3Nzc2btw4tm/fPsF8R/rSFRoPgRDC0TkEQghHgUAI4SgQCCEcBQIhhKNAIIRwFAiEEI4CgRDCUSAQQjgKBEIIR4FACOEoEAgh3P8HmQISOEEJLbwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 200x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "kt = 0.596161  # @ 300K in kcal/mol unit\n",
    "fbias_weights  = np.exp(metad_fbias/kt)\n",
    "fbias_weights /= np.sum(fbias_weights)\n",
    "\n",
    "plt.figure(figsize=(2,2))\n",
    "plt.plot(fbias_weights)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "accb22a7",
   "metadata": {},
   "source": [
    "### $\\textit{Load the GMM objects for different choice of weights}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7dae7751",
   "metadata": {},
   "source": [
    "Loading uniform objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "5cb2fd98",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(5,)\n"
     ]
    }
   ],
   "source": [
    "objs_uniform = []\n",
    "#for i in range(3):\n",
    "for i in range(5):\n",
    "    with open(data_path+\"run_with_v1.6.1/wsgmm_%dstate_wout_weights_nattempts_40.pickle\"%(i+2), \"rb\") as fo:\n",
    "        obj_ = pickle.load(fo)\n",
    "    objs_uniform.append(obj_)\n",
    "    \n",
    "objs_uniform = np.array(objs_uniform)\n",
    "print(objs_uniform.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb6fc13f",
   "metadata": {},
   "source": [
    "Loading objects from other choice of weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "041e5e1e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3, 5)\n"
     ]
    }
   ],
   "source": [
    "wsgmm_list = []\n",
    "\n",
    "# loop over different weights (w1:= bias, w2:= rbias and w3:= fbias)\n",
    "for i in range(3):\n",
    "    objs = []\n",
    "    \n",
    "    # loop over different clusters (2,3,4)\n",
    "    #for j in range(3):\n",
    "    for j in range(5):\n",
    "        \n",
    "        # load pickled object\n",
    "        with open(data_path+\"w%d/run_with_v1.6.1/kabsch_thresh_1e-1/wsgmm_%dstate_nattempts_40.pickle\"%(i+1, j+2), \"rb\") as fo:\n",
    "            obj = pickle.load(fo)\n",
    "        \n",
    "        objs.append(obj)\n",
    "    \n",
    "    wsgmm_list.append(objs)\n",
    "        \n",
    "wsgmm_list = np.array(wsgmm_list)\n",
    "print(wsgmm_list.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "e53f7f69",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4, 5)\n"
     ]
    }
   ],
   "source": [
    "#all_wsgmm = np.concatenate([objs_uniform.reshape(1,3), wsgmm_list], axis=0)\n",
    "all_wsgmm = np.concatenate([objs_uniform.reshape(1,5), wsgmm_list], axis=0)\n",
    "\n",
    "print(all_wsgmm.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3f36f10",
   "metadata": {},
   "source": [
    "### $\\textit{Lets create the modified uniform objects}$\n",
    "\n",
    "Hard code the cluster populations from fbias weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "ee42af4e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 5)\n"
     ]
    }
   ],
   "source": [
    "## making a copy of uniform objects \n",
    "#uniform_modf_objs = all_wsgmm[0,:].copy().reshape(1,3)\n",
    "uniform_modf_objs = all_wsgmm[0,:].copy().reshape(1,5)\n",
    "\n",
    "print(uniform_modf_objs.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "b4547703",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original Weights for each of the uniform objects:\n",
      "\n",
      "2 cluster:  [0.9187048 0.0812952]\n",
      "3 cluster:  [0.35271946 0.3324601  0.31482044]\n",
      "4 cluster:  [0.33975473 0.27061529 0.23671324 0.15291675]\n",
      "5 cluster:  [0.29853009 0.27504356 0.21523758 0.13976213 0.07142664]\n",
      "6 cluster:  [3.14788718e-01 2.59534446e-01 2.39811606e-01 1.85795758e-01\n",
      " 5.94719141e-05 9.99990000e-06]\n"
     ]
    }
   ],
   "source": [
    "print(\"Original Weights for each of the uniform objects:\\n\")\n",
    "#for j in range(3):\n",
    "for j in range(5):\n",
    "    print(\"%d cluster: \"%(j+2), all_wsgmm[0,j].weights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "224d8478",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original Weights for each of the fbias weight objects:\n",
      "\n",
      "2 cluster:  [0.55865813 0.44134187]\n",
      "3 cluster:  [0.5131445 0.4683594 0.0184961]\n",
      "4 cluster:  [0.4322426  0.3621806  0.18708089 0.01849591]\n",
      "5 cluster:  [0.28425596 0.27785684 0.22222823 0.19716292 0.01849604]\n",
      "6 cluster:  [0.41034481 0.30439397 0.16123026 0.10553483 0.01219767 0.00629845]\n"
     ]
    }
   ],
   "source": [
    "print(\"Original Weights for each of the fbias weight objects:\\n\")\n",
    "#for j in range(3):\n",
    "for j in range(5):\n",
    "    print(\"%d cluster: \"%(j+2), all_wsgmm[3,j].weights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "23a56438",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Modified weights for uniform objects:\n",
      "\n",
      "2 cluster:  [0.99434779 0.00565221]\n",
      "3 cluster:  [0.23604636 0.54791274 0.2160409 ]\n",
      "4 cluster:  [0.57484136 0.01834223 0.34939254 0.05742387]\n",
      "5 cluster:  [0.45953646 0.20583064 0.25952673 0.06031051 0.01479567]\n",
      "6 cluster:  [5.48468245e-01 3.43984125e-01 1.84959719e-02 8.90516554e-02\n",
      " 4.37242368e-14 2.78614494e-09]\n"
     ]
    }
   ],
   "source": [
    "print(\"Modified weights for uniform objects:\\n\")\n",
    "# loop over sizes 2,3 and 4 (uniform GMM objects)\n",
    "#for i in range(3):\n",
    "for i in range(5):\n",
    "    obj_i = all_wsgmm[0,i]\n",
    "    wgts = np.empty(obj_i.n_clusters)\n",
    "    \n",
    "    # loop over num clusters in each GMM object\n",
    "    for j in range(obj_i.n_clusters):\n",
    "        wgts[j] = np.sum(fbias_weights[[x[0] for x in np.argwhere(obj_i.cluster_ids == j)]])\n",
    "    \n",
    "    # assign the newly calculated weights !\n",
    "    uniform_modf_objs[0,i].weights = wgts    \n",
    "    print(\"%d cluster: \"%(i+2), wgts)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e794388",
   "metadata": {},
   "source": [
    "### $\\textit{Combine all the GMMs together now}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "c8f4a1e8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(5, 5)\n"
     ]
    }
   ],
   "source": [
    "# uniform -> bias -> rbias -> fbias -> uniform_modf (This is the order along 0 axis)\n",
    "all_wsgmm_combined = np.concatenate([all_wsgmm, uniform_modf_objs], axis=0)  \n",
    "print(all_wsgmm_combined.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ce9c042",
   "metadata": {},
   "source": [
    "### $\\textit{Function for performing 2D numerical integration}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "e107f0e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def numerical_integration_2d_simp(fxy,dx,dy):\n",
    "    #dx and dy could be backwards, but probably identical here\n",
    "    int_x = np.array([ np.trapz(fxy_i,dx=dx) for fxy_i in fxy ])\n",
    "    return np.trapz(int_x,dx=dy)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b7571ec",
   "metadata": {},
   "source": [
    "### $\\textit{load the Metad reference FE for BF=10 simulation}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "9b964740",
   "metadata": {},
   "outputs": [],
   "source": [
    "metad_fe_sumhills = np.loadtxt(os.path.join(data_path,\"2d_fe_phi_psi_sum_hill_bf10_nbins_100.txt\"))[:,2]\n",
    "metad_fe_sumhills = metad_fe_sumhills.reshape((100,100))\n",
    "\n",
    "#  calculate the nomrlaized probability density \n",
    "dx = dy = 2*np.pi/len(metad_fe_sumhills)\n",
    "prob_metad = np.exp(-metad_fe_sumhills/kt)\n",
    "prob_metad = prob_metad/numerical_integration_2d_simp(prob_metad,dx,dy)  # normalized"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "9682793d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAESCAYAAADXMlMiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACbp0lEQVR4nO19a5gU1bX2O8zQPd0z04xch4mKoEa8gShxHJMYjHxcJCYaYxIOHrx7NJBEMRrNMYIaHe+3SCQ5UTCP+hE9R02O+hEJiMaIJKCceAuPEgSMDh5AGGa6p5uZ6e9H9a5etWrvXbu6q2cGut7nqae7qnZV7erLfmutd621K7LZbBYhQoQIESJEgRjQ1x0IESJEiBD7NkIiCREiRIgQRSEkkhAhQoQIURRCIgkRIkSIEEUhJJIQIUKECFEUQiIJESJEiBBFISSSECFChAhRFKr6ugO9gZ6eHnz88ceoq6tDRUVFX3cnRIgQASKbzWLPnj1obGzEgAH+n407OzuRyWSM20ciEVRXV/u+zv6MsiCSjz/+GAcddFBfdyNEiBAlxNatW3HggQf6OqazsxPDYjG0+zimoaEBmzZtCsmEoCyIpK6uznozaCuQTgCdYk8brJU9ufU95H0naSPg5+cGALXkfYK8pz/AOsV72To/VnZuD/DD6XoNea0m7/n2WgDx3OuQ3Gsk193hVpuBw9sQqd6LwfEdiCGFQdiNIdiBCDIYjJ2oxy4MxmcYhF32vkHYjcHYiUHdu5H4373ATgD/C+AzAJ/mXncCSOa2i/dpa70tCXSmra7SbxRwf8rV0dwnNyx3L9Fc3w/IrQ8GMCy3PhzA56xtbcMGYnflIOzEYOzGIOzAEOzGIOxGPZKI4TMMRhoR7MRgZBBBCnEkEUMnYkgjghRiSCEGAEgl4wCATOdAx1cSqd5rv4/Fk9YrUvZrFBlUI4U4UoghiQgyubNa+2JIohqdiCOFCNL2vggy9jExdDr2xXPnH9idQbzDun5VEp7osm4ByZqB2FsZsd7bZ40hgyhSud7uwiB7+y7UoxPV2IV6pBDDbgxCCnHsJJ/f7p2DsLcjDuyoAjoAtALI5F7Tud9EO4A9bcALB+X/5z6QyWTQDuBqWD8BL6QB3NnaikwmExIJQVkQie3OqkhYiwMDAWRhDTm1AHoApGCNJm25VzEcxXKvKc3VYuR9N/JDWBecA75o10O2i2o1dWwdkJNFTLLNg1Q6yfsYrH8GYBFEEtZHkAJQSdpX5bZnYH1cGVjEsRfWbXXn2nfnjqkC9iYTqEi04bOqOgyIb8deDMRu1GAotqMbQDcGoAPVqEYkN+AORA2qkEElejAASOxCoi5jXbcG1tcQy/VzR+7aA3PvOwA0AIlPgbYKINVpfYJtcFK5QB2AWO4nkUjBUgrjuXvNAKiHNap0w/p64gAqrOvFK/eiOtGJSqQxEHtRgb2oQheyyKIa3ciiB5+hGoMAJFGFvYijCnFUIoZKRDEAMdQASCGOmgSQTMYRdX1l1gAVjychfgsDUJkb7GtRiTSqMABVqMBAAANRhQiACLKIYACi6EF1rj8RZBFDN+LoQgRZxHM9snpTiRgGII4BiMP6QCLdFYjnvvuqSt4vN7pyDxlVNRXIVFrnqMSA3FKJNCoxAFVIIYE4BgKIII06VCOKDgxGBBF0I4aBiGM3hqISMaRRb/WoKwEMiAPpKosw4rB+m9Ww1iO5JeeVKsZtHYX8ES2EGcqCSBywOaANbkJoU7yXnkCznw7we5AnhjbkB3pxnhi5VoIcAzifo3kb3pcYayegIRbahU7J/lq2vRbWoE3BHwKJ0ZZptx5XU/G4vS2JuKN5LPfUHEcSScQRQwpJxBCvTKKrJoOqTuQtIo5P2fpwi0xs5PpOPxFP2y0N6x6pJdYBa6TpAKqqgUhNBvHKFFKI2/2OI4UMorn+x+17jcP9WJ9EDDEkrePjSSSTcVebeDx/XIycI2ozf/CIdJvrBECeRCj498u3JSUPP9JjxGfS7jFEef0dDVEH+WMZx0DvJmWJ8iIS2WCJNuQHbsD9y9TtU4GO0OIcMjIRbTkJcEIBOx6sHe8b/UsYDKOcUMQTHyAnE4o9UKO6ChnEkazOAHEgjhR2YAgA5NwxSXtdigRQj3brR/oprEG+BnCMpR5kkuo0dPyJcwt0wPocamB9FjWwCSbekUEmkUQSMaQRRRxJZBCxB8lUzo1FSUSQZH7dSSYqCBKJk2NjOYdZf4X4HNKI2m48Dr5dkEla5mASv8U0ew3Rb1BeROIbOhKRWSx8yKIkwcmEtudWjMz60FkpvG2BpCK6wcnEL8RYUAsAVWhHPTA07+8XT+o7MBRDsB0pxOwne8uPTwZJSiYcHbD0CwmZoEPybeQIMZbzYST4E3UH8laXxBop1CoBYLcRiOcsL0EmMnAS4dZIPKeP6EglniMd0S6eUzCCsGyqOuRWCQDlPflpJyxaJQr9fUoQg5lFEkKOMI/Ehs6tZUIiqu30WP743sbaya7Dz0kDArza0vOanJs0FVYI/bN25tbbYQ2sKfJ+D5xPjHty+zoBtFch0xnBjuRQW4jdjqFEjI5jB4ZgF+qxC/XYjqG2BJxEHMmaSN46iMIa7AfDFvZtN1QtaZPblqjJE0asOk8iSoj77YDbghVuLwbSU0RyH0LcJgH1IC4IwhK/3QttI45RWSMxJBHNXSdmk4eeLGIaAioFZC4sGUxJyIbkO+kPaGlpwRe+8AXU1dVh+PDhOPPMM7FhwwZ7/86dO/H9738fRxxxBGKxGA4++GD84Ac/wO7du7XnPf/881FRUeFYpk2bVurb0SK0SBwwcV2pSITuN7VMZO25W4xeU2ahQHI+AVk/BFTuNNZU5+aSgXomFHqJGByFBZLvkVMvsVEJYPguJMB8+MIi8UKHxPoQUG2n16DuLQDodLu3ANhRUZncOr2HDKyIJkoClltMfBbu52EZieT3JR2fXaGIG/3mewdSN5j4vfEwvH5KHhQvv/wy5syZgy984Qvo6urCT37yE0yZMgXvvvsuampq8PHHH+Pjjz/GXXfdhaOOOgqbN2/GZZddho8//hj/+Z//qT33tGnTsHjxYns9GjWJOSsdQiIpCWRkQuFFJoDb3SXaQdJW5vai7WXHcMKS9EFFJoDcrcD1C/rbZnrJEOxwDLSWi8c9kEjFd0oe2+BNBgKmg4/QS4RO4su9lczZE3mhPUUIU9wzJQVKKBycRLg1EifWS6Gg54x3WP2skn1WYlD3sOrEvUj1DgkEyfYlRES7FwwC2WwsW7bMsb5kyRIMHz4c69atwymnnIJjjjkG//Vf/2XvP/TQQ3HLLbfg3HPPRVdXF6qq1MNzNBpFQ0ODj96UFmVKJLKILQ6V+8jPNVSiuji/CZkAxRMKPYYfp4seg5xMxDoHH3xq+fsqZGqtQSMWTxIXj1MnMNJLBGmNkPTDL2REJCwRHjVGorpUoju3SmQDvRDfU4jZnwMdeKmGwUlEZo34dWv5gpcVSkCtySjSSrE9MPSRZdLW5hwLotGop1UgXFaDBw/WtkkkEloSAYBVq1Zh+PDhOOCAA/DVr34VP/vZzzBkiCZwpcQoUyIR8EMMhZ7fL5kAZtaJ7Pz0vAI6UtFFjxmQiRdo2xypeIUER3PJdXQgHIrtiCKNJFKWFYCM84dLB3qTQY8OPuJYr4fnNJzuLYI4nFaJRYJpZfivACUTqwvuwV9GIvnr6q0RldAugzL0V/Z5imAMOAX3SHcGmUq1dSGsNde1kTHWT0qFBMwsEvG745Uy5s+fjwULFiiP6+npwRVXXIEvfvGLOOaYY6Rttm/fjptvvhmXXnqptg/Tpk3DN7/5TYwePRobN27ET37yE0yfPh2rV69GZaUfmyk4lDmR9AdwMgGCsU7o+QVUri9Z9JikD2QA8SSUOnLpKGkvCQkWg4gICaZPsTS/ZAeGWL6Fml1IdGScg/o2yMNChThPQ0j9uMJI2K+9jbm3kpV5d5WwSgBItRJ6X6bgJMItNWqNmCCKtJaApG4tDvpbgGWdJWvcJEIj2PYnbN26FYlE/v/hZY3MmTMHb7/9Nl599VXp/ra2NsyYMQNHHXWUlpAA4Lvf/a79/thjj8W4ceNw6KGHYtWqVTjttNPMbyJAlHnUlk7HkD3J+yhFooSpuOk3IVIVtUWhi/iSnb/NvVlEblHQSC56qTR7ZVFcyaQY/ixNQURxfZYrm5FBFEnEsQv12IEhSCOKJGLIVEbQNjxilWcRUVojYEVy1SC/vQbOSC/xXrYAZsI7RU50p5YBdTvRCC5RuoTul1kZFMKKoG2cwQnWefkxVsCCd64JbyP0EXp/WuT2exGPPMrM/TuWEqFXlJ3pQ0HASCQSjkVHJHPnzsVzzz2Hl156SVoLbM+ePZg2bRrq6urwzDPPYOBAf2mPY8aMwdChQ/HBBx/4vo+gUKYWCc8uLyW8hHdAbpXojlW5usQx8LimV+a8wjKhm1V6Cf9IqTVSl3tfbbm4ItWZnEvIHcUlXFwC4ik8ijS2AxhaucMpvn8KeT0UFfjgJwYk1TmEJQI4s99huXSE6C5yYrhVIkMml7To5QLjJMJdWlqXlU+3lk0KKhKh1pkGwt1H+8j1I8AiQ3e1A/c21MFpcYpIuoBQDbM8kh4f58xms/j+97+PZ555BqtWrcLo0aNdbdra2jB16lREo1H8/ve/L6h+10cffYQdO3Zg5MiRvo8NCmVukXD0VkqS7A+tSxE3PQdFoRaKOEZxfppjAhArA+78EsBpjWisknTO+rAyxON2gUMxBGYQwXYMzRX5s4bGXYladFHrYzAsIhB5JCrLpJZt5wOjzktBB68O59N4jAzW1CqRkYSwTOhAz/UPuo1bMVY3M65jvKwR4dbiUWIua4SigyyAk4TJbyHekbFIFaLQpLwPskAAWX8jtZLj/ddl7DPMmTMHjz32GJ544gnU1dWhtbUVra2tSKWsz76trQ1TpkxBR0cHHn74YbS1tdlturu77fOMHTsWzzzzDACgvb0dV199NV5//XV8+OGHWLFiBb7xjW/gsMMOw9SpU/vkPoGytUj8Iob8wJpA6UV6ChOLxut4AZ2OIosgEyaIQi8B9FaAGHDok7ywSuAW3qM5gRqA62mWZorbegksjSLRkclbJSKSa5umXyaQEQzVV8iTOY3eipI8GDE4imx3CirOR5CR5pkATlcWHaCjDheXnDgoSckQR9LbGlG5rZhlospyFwEE4j6T7HsVpWTE/Vu/AfZARzUuGQJwb9XBzKD18+T90EMPAQAmTZrk2L548WKcf/75eOONN7BmzRoAwGGHHeZos2nTJhxyyCEAgA0bNtgRX5WVlfjb3/6GRx99FLt27UJjYyOmTJmCm2++uU9zSUpqkbzyyis444wz0NjYiIqKCjz77LOO/SYZmjt37sSsWbOQSCRQX1+Piy66CO3tAdq0NhKS96rHnyC0kiDg1zXXBrWlorKIxDUkegmQt0T8WiWAXZBPFOgTVolYUuT9LtTb76V6ichwFzrJCKgtE6qVUMi2yaB6Ircz1fPRUkDeKhHbhV5iHZO3OGSLaCMjEe6u0lkjIqNeZiU4rBETEuH7JW4wWb6L7DUqIUoA8vpjUfbqx5XZB8hms9Ll/PPPB2ARjKqNIBFxHnFMLBbDH/7wB3z66afIZDL48MMP8atf/QojRgQRB184SmqRdHR0YPz48bjwwgvxzW9+U9rGK0Nz1qxZ+OSTT7B8+XLs3bsXF1xwAS699FI88cQTAfaUhhipQK0SP/BDOiqtRAedXqKDzNLRVSpm16AhwToUYJUA+cFTuLgAp5tIqpfwa8oskxp4D5A+o7qqOvLRWzEk7ZwSbpUIeIUEy9rF2MArc2lR6LQROyxYZo3Qe5O9F5BYZlUdQBxWNF2mMgKemKjSSUR/83XJiL5S2wV0VlmkIfuLxgDskmz3CVOLJJxfVY6SEsn06dMxffp0bRtdhuZ7772HZcuW4a9//SsmTpwIAPj5z3+O008/HXfddRcaGxv9daga1i8hBTgHUfqELghDuLBUJCOO7w03l5d7qxgyATu3rlKxoh86A1H8O9NkfQ95324NEslkHLG4NQiLkFEqvANuF5dAEql8siIV3wE3KZgmsMlcWxQ8+93uX8xOxBODvcgrEaD70og4CjfKrQk5iXCXVpxZQhRG2oiwLGRaCAcX3Vk4ML92/mEgTxjie6biegypfMmZ6oy7cKP4/dR59C9Er6LPxXaRoXnEEUfg8ssvx44dO+x9q1evRn19vU0iADB58mQMGDDA9i3KkE6n0dbW5ljM4GU9+J1Iyq8LrK9cabLCkBwKF5dXiCgt6ijW0+Q9LKsk05mfTZCGAwvhPUmer6n7S7i4RHFHh/guwoJZIUdHoUfVImAYRCMGZO7esrZZAzh1cQkIclAVbtSRiIBwaTn6o7BGxD6XNWJCIh2QWyrExVXV4RbdeX9Fn+n98PcOF5z4DurYaw17DdFn6FMimTZtGn7zm99gxYoVuP322/Hyyy9j+vTpdsRCa2srhg93VuWrqqrC4MGD0draqjxvS0sLBg0aZC/e87Un4BzEY2Q7FPv48bJFhVJEh6VQXDizn4rGCsi0EgrqfaEzGku0EiBvdaTsV5FbErNJxpq29QBbL7ErBQ+HPIeEkgklFAq+TWWdMI1ADM4xOMu/cAuB6iVCM+H5IBQqEuEuLXEtbtVQa4TnjVTpiEG8lxEIb6uwDsTUvrT/Mp2EuywBSyexI7dUfqeAdJI6qP/FdNmHgsZ6FX0atVWqDM3rrrsO8+bNs9fb2tosMqlB3snpcG/RAVgVlUVdXIXqJeLYUsKrX7rrU9eVzMWlcKHpIrhEd8RALM12d+aVWLudT6wxNohyF5cooYJK5PWSYp9UuTWiEeKrRP0t5MuA0KxuOoALfUCURxEuOx2ZyEgkv88tsAuykpVdiXRnzAR2E7GdZ/53WoOK0EpE1j9191GdhH5GnjqJdSLrdxNF/rcUzpHb5+hX4b80Q/O0005DQ0MDPv3UOWtRV1cXdu7cqa18aVJATQ4ZWaiIpRAy0Q3iveXWUpVakcEgUdJUnlGJ7oA1kNUCmc58QUeulVCrBMi7RiK2aygfEqydDMsvTMvMQ4QBO0VmL2GdkonYxqEiEe7SUoX7UmtE69KSWReq2o/ie6RkArg+L9GXFOKO75PXIxPuS0GqLp2kGnl9pAQ6SXUUiBko6XuzCGdolKDPNRIKnqHZ3NyMXbt2Yd26dXablStXoqenB01NTf4vIB0UErBGQln4r87FRfeboBASKSVUrrAigweoe0sGmXsLcMzNTUmDaiVAfqCVubhcegnXPUTCIg3/pX52lVbiw7KhA7XlUlJlnKddBAE4ExGdYb9uEqHHcPcZDfeVCexSlxZ/rxsw02Q/P55oJfn+6N1b/D0NVXYkJqp0khB9ipJaJO3t7Y76L5s2bcL69esxePBgDB48GDfeeCPOPvtsNDQ0YOPGjbjmmmscGZpHHnkkpk2bhksuuQSLFi3C3r17MXfuXHz3u9/1H7ElQN0v0v+4Vygw308JIqXYroOORHojX8XUpCjAvSWp/gsg/1RJ21TnQ4Hjcbd7SGWV8EQ32yVWmbRmSJRlbIsnaBnZySoDWydWg0RuAe4qtzKrROwXT+CCKFRzstDz5IVqHrXlLbC7XFrUChHvKYGonvjFZ5NGfq4Wuj2HSHcGqLTu18u9JeZtiSNfuDMeT9qWqp2YyN1b/TyfpBxQUiJZu3YtTj31VHtd6BbnnXceHnroIaMMzccffxxz587FaaedhgEDBuDss8/GAw88UFiHqqHxRslmDGyD3MWlIhs/FoqXFdKXSY8qrUTVxgC8BLsYCARy7i0gHwoMALx2FQ0VlT3d5oV6Sy9pGw4kPi1yJkHRzxqoBy02b4nX5E5if4bknACKwoWQk4iAH4EdIC6tQkmE7hPzwkTZdsChU/HPg5fbB5zfrQgDTsHSzjI8n4S6tz7T9NMQiTiQMPHP9CB0bUlQUiIRmZsq/OEPf/A8x+DBg4NNPuQDgS2685yKPfAmE9HOD0zcWL1NIoXmofg8vAPO20+z9VxOiXVKd9kUapWIwZcORM7cklyHqPhu0j+V5cGjuAwhE7ut/ol7SttkQrdTuIV0GrXlzAxXZbALa8ShiwByVxbd5pUjJDQScSzXSDr0onsy995yTUbI9rh9H0nEEalNItOesPojCjiK1/2vSv0+h34ltpccNdBkwfIILpkLSya+05FQRip+9A+/BKJyq/Vj7IFzsKbJichHbwkI0R3IV8sF8gMutUbo06yDcGoicE2GVQjoQ4iEcERmdyahntxJIH8fMZtM6HYObonIXFrua+QFdldRRi6ucxIxqULUDjmZEAirRCa6iz5y95a4J2P3VhCogZli3INALKD9DeVFJIA8VNA1BnMLhQvQqqz2QkTzYshDtj0oQpG5rjx0ElNwSwSwy8sDzugtIK8bULIQobIya4QOqsLFZQ1mPshEpe8YWCZCF7D6Lu9XvkhjyiYTa7t7ZFSRCHdp6QR2QOLSArxJROfeEqXcKZlIQK0SXkJG5d4SupGReyvUSPoc5UUkNfDIxFZZJdTFRfcXWiYlKPJQtQ3aOvFZA4w/zcoGGNkA0OncRucqAZzzlaRtzUTklbhLcNDtvsikQ9I3eh9Rtk3j6qIuJiqkRwgRCjKxtstdYTISye+T54xQa0Tq0tKRiC6ZlIf+UjKRwdZJ3CVkAnFv+ZhTPkRpUF5EAqifXhxjLxXeOZmI/TyJsRQo9Km/EDIp0soQEIMKRYfi1DLXBNFJgDwRUPeWte60VgDhJsm7laiLC0BhlgkHFd4LBLWsxPs4I0RnezmJmOSMKF1aspBfXo1AJSqL7TRaC3B+75JijiLHhueU+IreosmJwi36v4p++kEcthWpRbd3k3JEv8ojKTl0UTc2dOXk6b4Ygs9Sj7Gl2HOVsj3gssR0T4YqXpMMVqL2lvPwmMPSSCFfg0sg6WoTs10jYp8oo9Klq7kV1SwUPmty8SKLssmqrGHUudA2smN0JeIB5tLiRCEjEZojAjjJh0d40XbiHGl2nKSEDL0HWupFfFP03sRnYGtnvPZWP3VttbS04Atf+ALq6uowfPhwnHnmmdiwYYOjTWdnJ+bMmYMhQ4agtrYWZ599NrZt26Y9bzabxQ033ICRI0ciFoth8uTJeP/990t5K54oLyIRqEU+U1a8d4yjCbJBRyaAv4GfE0WQxKG6Xj+CSsDdI98nam85tiHusDoAN8kItwhNVKSWiZZMAL27SrNPNrGTgMztBOitC75dNte7rkS81KUF5Ad7FYmINpx4+DnoMfxcdBvyCYriW+ElcGTv87W4JMmJgjyCEtt1Dw+6hwkNXn75ZcyZMwevv/66PQ2GmBFR4Morr8R///d/46mnnsLLL7+Mjz/+WDnlhsAdd9yBBx54AIsWLcKaNWtQU1ODqVOnorOz73x85eXaqgYwMPdeNqi5PEI8Uou6uQC3NtLPBm4AwYvwRYKHAFMQwV2A194S4KI7IBeqncJ7zLHNEc3F8kA8UUTlWeraku2TQUciupwRG5QYZFYEIA//9RLbPcJ/7eNz32ukO+MQ3U3dW3GkkOLurX4uti9btsyxvmTJEgwfPhzr1q3DKaecgt27d+Phhx/GE088ga9+9asArNkTjzzySLz++us46aSTXOfMZrO47777cP311+Mb3/gGAOA3v/kNRowYgWeffdZRv7A3UV4WSS17z60SB0zKpJhU+u0v6Eckp0u9UTxUWVaF+x6Ee4tbJQLcKqEuLpdlwqsBe5WZF6/kt5Os8Q79FaCWhmwOEb5PVlbEKGdEZlmIzH4VidBjOtkCSVtxDnretPMcoiwLtTS83FuUMGNIIlKdsawS6t4y/8gDA5+mIp32zlIU0+UOHjwYALBu3Trs3bsXkydPttuMHTsWBx98MFavXi09x6ZNm9Da2uo4ZtCgQWhqalIe0xsoL4sEcI6nqhBP++Gdh/5SywRwjoh9nYluApUIHwDJ+NXqZSHAAu1VyMCZTyKgCvMVkFklKsSRRKYyYgvwgOYPoSgBIgY06tbKVPojE5pwqWtHX3mUllZgp/dAB3sB7s4SUHlKVNt55QKxTqwSWthSiO68ZEoKMbu0vMgdsh4UYnLRXS8pmCEOs9Gwy3rhU1PMnz8fCxYsUB7W09ODK664Al/84hdxzDHHALCmyYhEIqivr3e0HTFihHKaDLGdT62rO6Y3UF5EEgcc/1VVXSjHeMtDgmXZ7YD/DPcg4ScMmZNJH1oqtKS8BDyfhINbIKryIrK2Ag4y6cigqyb35CyeeMWgqSAQIE8i1BpRXY+XjveCTEOgJGJUBoW6tHQ5JNQCoW1kIGXjUS1pV8O2idD7Grd7i+eU8JkTpTkloiKwCAXuZWzduhWJRP7h0ava+Jw5c/D222/j1VdfLXXX+gTlRSSAmT9V+sRFtQbZwO31a+4NolGVvOfoZ7qJzjqBWidxt1OTIh1opfWvcqHBVltSal17RTeJ+LFGeN901hYlEQHjMigcqugqkON02gjdLwiimm0H2U+2V1WLUOB8pntROSW1yE+rXAxq4csiSSQSDiLRYe7cuXjuuefwyiuv4MADD7S3NzQ0IJPJYNeuXQ6rZNu2bcppMsT2bdu22VXSxfpxxx1n1J9SoPw0EhNxVBrFJRBj201dWnWGS7Hw42LbR3STHIROIp5SvWDrISyKK3+uuC3Ac91EEENXjX4B1LqILCdEh7jdSz2JUF1ENs+IDZk1woNMdCQiC/vl0Vz0WBn5cHFfEgos7kuI7uJeRSAFzeCPI2XNnFidsdxbQL8tJZ/NZjF37lw888wzWLlyJUaPHu3Yf8IJJ2DgwIFYsWKFvW3Dhg3YsmULmpubpeccPXo0GhoaHMe0tbVhzZo1ymN6A+VnkQDeVon4s0ldXLoMd4pC5/UIwrIxtUz6GUhJeRPwMGAOniVOyUe4SjhsVxcgL0AozsXIg1oiMpJTVQHWgeeKAHDpIrx/WoFdQObm0s3ZrrMyaOSWTjshFgyfTVKI7vmJrWL2e9mEVy6rpJ8SyZw5c/DEE0/gd7/7Herq6mwNY9CgQYjFYhg0aBAuuugizJs3D4MHD0YikcD3v/99NDc3OyK2xo4di5aWFpx11lmoqKjAFVdcgZ/97Gc4/PDDMXr0aPz0pz9FY2MjzjzzzD6603Ijkijyg5RXVVPA+mO49Gme4Q7egLSTodgB3lSX2UfJBMgNSFVFh3Vy4Z2WJQGcZBKDew4RSg6R7ozU8uCurHztL+tcKUYqOvcbhY5EaC0t3wK7LHpLRSI6IpJBNaBTN1hOfwpUdN+uuK4fRJBPDdDBhw/noYceAmBVQadYvHgxzj//fADAvffea0+PkU6nMXXqVPziF79wtN+wYYMd8QUA11xzDTo6OnDppZdi165d+NKXvoRly5ahurrv5hwuLyKh9YB0gxR9MvYkE8BfFd6grZe+FPkJSuAlExNdAWrBnSLNLBRnDolTf5CRiQqcVGSQ1fcSJOLHGnGXPHGTiKqtkTVidUhdJkUW0isDr7dFj5GBWS5VACI1AYru/RC6KTQEqqursXDhQixcuND4PBUVFbjppptw0003Fd3HoFBeRAI4I4RkVkk1268lE8A7KdFE0C6GXPoRmZQQ/MleBlmJkHx5dmetLhmZWMUE4w7rJI4UTPQYwKmJyPprao0IqCe5IpNa+bVGODGoJriCpC0k+2i9LRrJRdfFuWlb2yoJoJCjppvGqIGZRbI3iIvtfygvIqmFVXQtTdY5eG6JlkwAJwnIBv9CiMXkvBTlQSZ+QEvPA86sdhmZyMDdXQK0Wq/7uvnzUUuEk4jJBFYCvq0RLwjC0IX5iv+IygVcS9rx2RFlFkotOScLBY5X5mdP1GW6A/l5SkSmOwBkasorZqg/ojy/gSissZcXcaRZ7rVkG5C3VLQPlSZRXIXU1ipV9nyKLPs+5IUc4w6tROzLkHbilUZ0CfcJXfLncH53sjbFkkhR1ogsUkund6jcWTodsZ3sl+WncFKjfRChwCzTnb+nwRIx5GeDpJnuADCwxtvtGaK0KC+LJAIrX4A/vIv/Of3hU1FeWCa0vXbsNbUm/Forhc5/wiG7TkBl5PsA1MKgSCFmDzrUxSWigHRRX1w7EW4vHbgWInNlqSwg0ScOY2vEBFQbkRVyFG0A78mtaki7WtaOu7TotjS0orsf95YQ3dNBPA/HYVZqJRA/2v6H8iKSOgBJmKVriD+SIBFquSjdXDL4cVGZJgrKIrJK4d6SzZLI+0HQy6Inn25XIM2IhQ7SyZxbhJKJmAxL5eoyEeQF/JBHMXDljZiI7BwqTUQ1uRU/lpIJ4HZxAW43Fw8XlmS6+yrkiDgi1XtD6aKPUV5EAuR/7F5jLg0B5hYJ1UwAH54hv5ZKKVxOunPuG1ZJPmTXfS9CLM+vx+yBCEBuWlc1mQDuiaf8iuQCpSAQI5iI7Ly9gAmJ0OOo9SGrACzOQ0mGiO5VuUKZxeSUdARRSz4OsxLx5TdiGqGkGskrr7yCM844A42NjaioqMCzzz7r2G8yQcvOnTsxa9YsJBIJ1NfX46KLLkJ7u855qwH9gdfBrZXE4NZN6FO2antB44wfPaUQyIjKhJiKJK8Sl/TWZYvzgZtrH0DexUUz2kX2u3B10cz5fK553LGdX0fXzmu/qv8csgREX1C5tQDvCC3ZIs4JyDPm+Tk6JO878/dBs9pVugktXmm7/QxCw0OUFiUlko6ODowfP14ZI20yQcusWbPwzjvvYPny5Xa9mksvvbTwTgkC8YIQ3uk6oCeToghFRyqyE/P2XiZWLwrqJSITnUbB9Y40IQWBvNCeJxNKKOI8nFC8iIHChDB422JgrI9wyNxgMmvEJASYt6eiukx0p32AdQ989kRAXzJFgJNrwaj2sYRwoaSG2vTp0zF9+nTpPpMJWt577z0sW7YMf/3rXzFx4kQAwM9//nOcfvrpuOuuu9DY2OivQzVwah8Ce6CekEfm4qKuLuHmEijKK6UT041FGQn8HkddXEIn2QNftcAC/sMlk3HENVWA6VOrcG+lEbEHH+ristxZFplwVxfgDBG22jgz4k2hEvN5uXiZ0E4DBTxRzMR4XomLfsC1RNl1aiTvc78VL/cWIEKyuXtr/4g63JfRZ+G/JhO0rF69GvX19TaJAMDkyZMxYMAArFmzRnnudDrtmnjGhuyHrhofZdYHtUxoaDAfOIvxShUd6msS1dVGlhJBfEYlrIWkcmfJtqWY9QGorRNqLQgrxe+igmw/D0UWMM6KL9QysS7iPg/fJivYSNvJrBKvfBWNe8t6L596l273qggdonfQZ0RiMkFLa2srhg8f7thfVVWFwYMHaydxaWlpwaBBg+yFT0Jju7dkpEK1EroNcJMJfy97Ci/K5RUE+NOajDxk21RPeQUSTwx5UjF1LxKorBEKOjjny5M4XVwpRhbWcVEXochIxa8G4nWcF+FweFYTVuWP8LBfGemoZEeTelt0LNedR6ankPN4VQQGAFlF4EAskjj0M2KKpY/iJ/o79suExOuuuw67d++2l61bt1o7aKKhABfdKajw7odMVITiG5xMTE5CB3oZiZge2wuog/XZBRB0Q5EnDTWZeBGKtS/mIhZ+nUI0Dhmh0P7JggQCRxAuLVkBSC/ILBiYJSfyRETrfejW6g/os2A2kwlaGhoa8Omnzllrurq6sHPnTuXEL4A1W5l2xjJa3kG3P0XWZU9adLsqRJiiZFG9JvkjpiQhyx3xyicJAIJUfAj1VEdQJfQJCL2EtuPHA0If4VWD07k2wQzscUmfeT4LBw9rDgReROAnBJiWSlFZJTSnRJRMkSQnAiIBNF/Whpa4oXlCVnJiAFmCNTDT9SqLv9T+iD6zSEwmaGlubsauXbuwbt06u83KlSvR09ODpqYm/xet6XKuC/eWl1UiwK0Suo2/F+2KdncVOoBTtirU0uijpz3FH9pkIOVP986yJe7wXmqd0H1cRyl2oZC5zmjfaf+tfgdsslH4sT5kGgkg1z1U5xDgEWJkwiu/YcD9GV4pEBUVFdLlzjvvVJ5zwYIFrvZjx44t8Z3oUVIiaW9vx/r167F+/XoAlsC+fv16bNmyxTFBy+9//3u89dZbmD17tmOCliOPPBLTpk3DJZdcgr/85S/485//jLlz5+K73/2u/4gtAT/RRHSwV7m46Db+3uuagXsuTAljj2QxOY9puwBQ2+XdxhAyMgGcg7UgFBWpmITz6iAjFZWFw/slPV9uJscu0ydpHagFIXM9ychBRiiclHgNLtV5We0twH8YcH+FVwrEJ5984lgeeeQRVFRU4Oyzz9ae9+ijj3Yc19dzwZfUtbV27Vqceuqp9vq8efMAAOeddx6WLFliNEHL448/jrlz5+K0006zJ4B54IEHiuuYGOxVHiFe+oGGLFJQ95XOzcXbUhQc1evHVUWhummf4b29gEit2dMmdwfJ3ETUNSTIRLi6BKjLi4L74Qslk7wbSx56rMqyF0hJ+mrDZE6QQmB6PprlzkPnxbr4W6dJW+beAtyuLPFZWJ+XMwzYOm0AlnMEZlpdj7/T6lIgALhc9L/73e9w6qmnYsyYMdrzVlVVad37vY2SEsmkSZO0k7uYTNAyePBgPPHEE4H0Z2BNEnuTg/MbVOMmHWvpQC/LI1GhV8hEBdXJCqnFFVDZFB6xRYV28V6CSLXzqVOXgCbTF1RkAjgJRbTNX8dZZqUYcB1GnD+DqJZM+PVLGupaTGIjL5OiGpDpPCVCJ2EoRCdp64OyPo6UAhjosgbYtm0bnn/+eTz66KOebd9//300Njaiuroazc3NaGlpwcEHH1zU9YvBfhm1pYUflwn9oXPtxMvFJVvn7QOBjhzoj92ERGgbL4vHkPlqYM5DooQ/g0noL4dKb+Dl3gGLUMRCjykkpFcGldvM6qfl7uIZ9rQPvhBEzo5fUim0vUInAfT5JPx9INntJqG/YgFw0EEHOVIMWlpaiu7Co48+irq6Onzzm9/UtmtqasKSJUuwbNkyPPTQQ9i0aRO+/OUvY8+evpuTqDxLkJnM265yZ3lZGny/qWVStFUiBv5i/1QqF1cBkVuFlkopQh+RRULx7YDTTSSzUihMooJMBvw48gUgaSFJ0S+eYc81ACq6a0VmOhNhMWDnaGPrCVk1CJ2LTUZytIhjbr+oBgzwzylluwWjyHiW9S8ltm7dikQi/38o1hoBgEceeQSzZs3ynHudusrGjRuHpqYmjBo1Ck8++SQuuuiiovtRCMqKSCLVe7GX/jdrIX9Q5xXZVYO8SViwqp3KzRUYvKwRekMlcA0UWW+L6yOFhr5yMhFQkYrsWjJy8QPuOlMRCnVzCSQdhMe2VwKoARKy4o20HJApNJ4zTiJim4tMxHlk4yrRQRw6CUEVIRRRVh6AQzeKuci2uO+nECQSCQeRFIs//elP2LBhA37729/6Pra+vh6f//zn8cEHHwTWH78oO9eWS8ClWe6yjGuZe6sQl1WJq+L6A2fFgLPY/bjvNPo+10cKhSx7XOc24jMjFrJQcNcZDz8W25zZ9XHjrPeuEpagAeQk4gKLwFJCVi5FQmC6xESBQPNqRPKx11Kiz/rhhx/GCSecgPHjx/s+tr29HRs3bnTk4/U2yo5IAAQaWqq+hkEbWX2uwMGtERVp0O1exxQZ9ssz2ukUxwxUHxFP6vlZA+X3IiMIVTmSQnUPL6iIhROK1dZd5p72W/THmikwKnfp5Ob1cKCYhxdD15gRyZieW2S3EwuLBhioyCXwRM0AoUuBEGhra8NTTz2Fiy++WHqO0047DQ8++KC9/qMf/Qgvv/wyPvzwQ7z22ms466yzUFlZiZkzZ5b0XnQoK9eWb1CdxMS95eWuMon20kJ0QAzke9i6lz7SCwmGusHLpMZWbVdgbi0ZZFV9VSiWTAQsV5ZTk+GZ9jQ6KWOghQhx3nUPXNsTWeclhO3iohqJ7Jqyp/l2OGdN9Ik4kugIwrVlWkerwt9pvVIgAGDp0qXIZrNKIti4cSO2b99ur3/00UeYOXMmduzYgWHDhuFLX/oSXn/9dQwbNsxf5wJEWRFJLJ7E3q4h1ko1ChvUaehiISHAvYJiMtm9zKIC8018DhLFurVkeRgcfgomFgKVNsPL3MvIRICTGbfCMpURt07CyaMI8b0gi0MFLsbLHio6823EVLrW+3xpf/pZ9Xd4pUAAwKWXXqqdY+nDDz90rC9dujSIrgWKsiISF3QDfSmmQO910Bvwa430fYJiIW4tDq8aXH5hGp0FuImKCsa645zE47xPOoDSfV01QBW1hjl5qKIQSwGddUEFd7qNkYqfyK19IcN9f0fZEUmkOoNMO/kzyyK3onA+0fWZZVFqyEJ6uVXiI+xX5dYqcaCBF1EE5aIyhex69Ik64rBC8kmSNCoJcOsl4jz56zD3lhig6cRR9Hdbg6IjBVPk+BghBKl7yws8csvjWFXkViAwzWzvBXl1X0TZEYkU/MHbryXiRyfpE+jE8l6o7KsCrfhbDWkQBNVH+mISoyAGK1rmg5dBofqJqmSLqj/is3G5t6j1UQunJVKNwCyTVKeTTIxB3FcyqEKAZYghid19mE8SwkJZEkmkNolMex8NniWDSog3QQBk4rcwpWJsiNQmPcuiULdWb7qtCvHL88RDlbVCyUR3LVn+i9jmcG9RnYRrJio3F7fEGVKFPCCJ66i4gJeUJzCpuRUYSNZ6CP8oKyKJIdX/jIWi4VdYLyZ0NyDrxcDVZVIWJSgSUeeTqInDq7S7sJ5k5xCDoltA1lsjvK/iGkmqE1CLI6p4DzgHTRo5FRQ4gZi6sWQaigK05laIvkVZEYkDtV1Ap8Htm4QABwVf5za1OExPWgIXl2yg8KHfy9xaXiK7LuKJww+BqIhDVgo+jpSrPS1AyPurIhMVxD2Ka8SRcrq36G+W6iR+yv4YEgt1bykz3QUKDPFVFW8M0X9QvkSyzyDIOT9Mz+Wz4m8J/9eqgnz8iV32BK9yJanAB3pOBqYDGG9HiUUMiiqxOE8SPOIr73PiorxVVoW5t1Q6iU/BPVFTwhBggXb3NloF2AuBztnuBX0kb9mivImk0FwSgf02mkugiBBgmfuqyLp2YsAwIRG6TzZg8206EuHE4KdYYMzh53dbKqo+8iRMHsHlrNeVO74SiCOTd2/RnA2qk/DoLo5i3VymlgcP+/U4bl/NJelv2LRpE/70pz9h8+bNSCaTGDZsGCZMmIDm5mbPgpEqlB2RxONJZDqL9Kv2Zkx+r6KPIrhyEVsyoR3QR2uZ6CRelgkdkHRWCCeQZFJPKPG43FUlCMXtsjFxaznnNomxpL1kTUTt3hK3RrcLUKtA8XHHquWCu+/oLS8dRJGUWFKYWiQ+J7bqT3j88cdx//33Y+3atRgxYgQaGxsRi8Wwc+dObNy4EdXV1Zg1axZ+/OMfY9SoUb7OXXZE4gv9IikxgZJOadtbUBk2ssmN4knXE7nMGinlfN0qApGRh+zBJFKdcbRVkQonE0CecU/dWPY5JVZJBBmnewuwiMHLvVWETtKXsKLcSjin/X6CCRMmIBKJ4Pzzz8d//dd/4aCDDnLsT6fTWL16NZYuXYqJEyfiF7/4Bc455xzj84dEUmr3lElJ+X0G3FpRWC+mCYge7XQTFslIhGY488HY5IlfWCMmJOJl1ar2C0KhLi9XG8m9cTcOvR9hlQjR3XZvAe75PkzdWwRUJ1FZJUqI66m+azrVbl/yQTXMosX2lrojpcFtt92GqVOnKvdHo1FMmjQJkyZNwi233OIqy+KFkEj6E3xFhfW1ueTDP11k9JaAygKRlciIIGNUS0vnZxckoiIQR4UEDSK1EndqnL6Vaycy8DIqwiqJ5uYyAYBITQZxZJyiu6l7C5I2PuCI3KKaB40cA9RJiT4nuAphBh2JcAwZMgRDhgzxdf6QSFTorTIpRWfC97brqw69qaPI9BHn/BTqOkucTHRWCbdGdCTiIpB2zd+otstuzwmFurtMyIT3n+aj8FBgl+gu3FuAM6eDkkWtZNs+4N4K4Q0+x7wOhUzYFRJJiGBh4h7w+TCpymQ3KdZnapl4QUoinED4A0G1s00GckIR3KCL6gKckWY0u1u8CqtEKrqrkhO9psbNbfft3tIlHHKLRWS3E1Spji8RuuJmE4R1dZe+L6VAfX09Kir0NfCz2SwqKirQ3e3/Jvt8YqsFCxagoqLCsYwdO9be39nZiTlz5mDIkCGora3F2WefjW3btgXbCe6/LSTitc9nQBQuGv40sX+GSJay4iu3Rlwk0l6VJ4hOsnDwfbnjMu1xa+mMINMZQTJpTX7l1Gby761JreJI5ya5AqxwYNFebMtbVda6PTDWkKWWbQNZFzOFivelgsrCUUWLlTCootR45ZVXcMYZZ6CxsREVFRV49tlnHfvPP/981/g3bdo0z/MuXLgQhxxyCKqrq9HU1IS//OUv2vYvvfQSVq5cqV1Em0LQLyySo48+Gn/84x/t9aqqfLeuvPJKPP/883jqqacwaNAgzJ07F9/85jfx5z//uS+6ug9BpqH4dYMFUEbeJ8HKIrYA7wgtvt/59G5mlcjEbymJCFDy8OX2tM6RQdxpnRBjJ6qoJ2WVThdVf9N2pBa1SgAgXpl0FnIsRnSXuLd8i+4qyER2H2VSgkKyZiCqarxnrUp2Z+FHce/o6MD48eNx4YUX4pvf/Ka0zbRp07B48WJ7PRrVs/hvf/tbzJs3D4sWLUJTUxPuu+8+TJ06FRs2bMDw4cOlx3zlK18x7nMh6BdEUlVVhYaGBtf23bt34+GHH8YTTzyBr371qwCAxYsX48gjj8Trr7+Ok046KbhO1CG4meSK0T2kgrsgALFTrActuPeC9qEgFp5DoovY4ihFGLA0R4RaIfY23oat17Lttbnjc26vDCcLjX5Pa0vFiOAs7l9oJcpQ4CJFd57lbkwmssgtThaS7Pb9AdOnT8f06dO1baLRqHT8U+Gee+7BJZdcggsuuAAAsGjRIjz//PN45JFHcO211xqdY9euXXj44Yfx3nvvAbAe5i+88EIMGjTIuB8Ufe7aAoD3338fjY2NGDNmDGbNmmXPZ7xu3Trs3bsXkydPttuOHTsWBx98MFavXq08XzqdRltbm2MxhtcTdBCeIn6Nop++BAGoOke3m5BFgTcZsHuPC+1iwKRuLRWJ+N2ugzI6q528tkNulfB94r1weeXISbi6qJvLmp895tBChIuLb6evwjWXqcxZOuL3xd1XYhv//dH9IO0M4SipQt+bWG3M8uHfV3+ZxIqPL+l04U+hq1atwvDhw3HEEUfg8ssvx44dO5RtM5kM1q1b5xgTBwwYgMmTJ2vHRIq1a9fi0EMPxb333oudO3di586duOeee3DooYfijTfeKOge+pxImpqasGTJEixbtgwPPfQQNm3ahC9/+cvYs2cPWltbEYlEUF9f7zhmxIgRaG1tVZ6zpaUFgwYNsheefKPEfhNRKMhC5ZqSkUkfltWXzEMSBExJgwrcNG/EIYhza6SdvQKWcahb2uEmFgMyASySyBCNxCKVmK2ViH2p3DHJHNEkayKWVkI1EfEaZdtUZEG286KMsWqfWe10vO2UbCOId/QeaeytjCBjsOzNkfNBBx3kGGNaWloKuu60adPwm9/8BitWrMDtt9+Ol19+GdOnT1cK3tu3b0d3dzdGjBjh2O41JlJceeWV+PrXv44PP/wQTz/9NJ5++mls2rQJX/va13DFFVcUdB997tqiZt+4cePQ1NSEUaNG4cknn0QsVtiT8XXXXYd58+bZ621tbeZkUgoEGjpcqHuL+8yoXqIjkb6dt0VXkM+PdVFo9JbDGlGRCP/4ZQNjlLSrg+T3UOUIFQYAxEWJkPz/IJKbXtbaLVxYTq1EmaDIQ4Flme4yMin0YdtP5JZ4LzhdkmfitwhnKbF161ZHmKyXrqHCd7/7Xfv9sccei3HjxuHQQw/FqlWrcNpppxXdTxnWrl2L//iP/3Bo0VVVVbjmmmswceLEgs7Z5xYJR319PT7/+c/jgw8+QENDAzKZDHbt2uVos23bNq1PMRqNIpFIOJZ9ClL+9LoH7t6SWSX8xAmD85YfXPqILE9ERiJpqAdduk9moRDLxNEX4uIC4NsqAYBkTY5AuSVCo7Oq4e2+yu3Xloo3BRf3Ve6wXoJtwXku1mfKx5dCiYRjzJgxGDp0KD744APp/qFDh6KystIVueo1JlIkEglbPqDYunUr6uoKC7Dpd0TS3t6OjRs3YuTIkTjhhBMwcOBArFixwt6/YcMGbNmyBc3Nzb3XqVKIgIHrJDr4+XEIsgkgYisglLKmlgrKEijclSWQZm1UCycU2j5HJioXlxjMMojYWgng1kioVpJEHJnKnHtLkIX4PctCgSnB0FDgYv8DaejJwwOyxNRAysf3M3z00UfYsWMHRo4cKd0fiURwwgknOMbEnp4erFixwnhM/M53voOLLroIv/3tb7F161Zs3boVS5cuxcUXX4yZM2cW1O8+d2396Ec/whlnnIFRo0bh448/xvz581FZWYmZM2di0KBBuOiiizBv3jwMHjwYiUQC3//+99Hc3BxsxFYx6JPaWSr3liq6S3as3+v1D5RKbPWca0QXqZVWbKdEU6doI/Y5trNoLubiEi6emP1qRXBZpUMirrIp2gRF6lYqMBRYBm2pFFm4r2ImxUh3BvHKXqoAXAK0t7c7rItNmzZh/fr1GDx4MAYPHowbb7wRZ599NhoaGrBx40Zcc801OOywwxwlTU477TScddZZmDt3LgBg3rx5OO+88zBx4kSceOKJuO+++9DR0WFHcXnhrrvuQkVFBWbPno2uLkufHDhwIC6//HLcdtttBd1nnxPJRx99hJkzZ2LHjh0YNmwYvvSlL+H111/HsGHDAAD33nsvBgwYgLPPPhvpdBpTp07FL37xiz7udS9AGwZcCKiWIgZNfoH+Qxj9Ae5SKOS9+ChlJCKTrDiptCM/VzndZiNPJvF4vsAjrTmVsl/zeSXWa8ZuCwBJpIBKoKsm464KDDgzzCHZD7JdEQrsC5QsaH9YCLAuuz3oedtTiKHSwEGT8llHfu3atTj11FPtdaHdnnfeeXjooYfwt7/9DY8++ih27dqFxsZGTJkyBTfffLPDVbZx40Zs377dXv/Od76D//3f/8UNN9yA1tZWHHfccVi2bJlLgFchEong/vvvR0tLCzZu3AgAOPTQQxGPF/559jmRLF26VLu/uroaCxcuxMKFCwO7ZqQ6Y1xwT4l+V4PIxCrhwryKOPqPW6vPwXULOthzEvErunNwqyXnekom48FaJYA8QZGXS+GWg7if3G/fiExMBHdxbnp9TVLivjTd7qRJk5DNqqdV/MMf/uB5Dlkl3rlz59oWSqGIx+M49thjizqHQJ8TSQgCX4mMOvcVh45MvI6Tvd93oZ/gSrGPfydeFocuyknsE4Qii+CykxiLs0oKTlA0rL9VEmjKycdyqg+PvgtCK7EskkqDdvtosS2Czs5O/PznP8dLL72ETz/9FD09TiurkFySsiKSGFKBJa/7hkxLKUhf8UMcKj2kr0vQ9w8EMlWrSi8B9JZHmq0rfwdVQG1hVomsbIpRVWCTUGAPq8Qosou61Dym2KU5NXzfnt6up7KP46KLLsKLL76Ib33rWzjxxBM9izmaoKyIhCJSm0SmvZ8/ZRvr4iY5JZyAVG11bi0Pl1c73O4IXvpCJrb2AkxzSJTT55oQBo/GEuDaA80nUWoksAd46oaNxa0INmGViLwSUYNLDLq8xLxtldC5SmjfxPdiIrozq6QovYSCVwHutJISMwnFJGEI520vBM899xxeeOEFfPGLXwzsnGVLJCVH0fOMmEDHNDKrREYmQN7Hwo8X1ygvGIX+8nUZiejCXWuQ/9gpmUjPX2UPsCLRkE/PG0OeYABnKLCvBEWviafovUrEd5cl4mWZCOLiLi2t0J5B0JNbZRBF2sC1ldkPXFuf+9znCs4XUaHf5ZGUGo7KsialOfx83kHUmjKy0nWWlKzDXrW2+lhc100MBaeuEcTcIiZQBmNwa0NlAPIaU7wWl9hPrRmRpEjft+f7IwiOl04B3FoJLTFvlKBI1wHv36FkkA8kURFQakxRpKU5RdF+Un9rX8Hdd9+NH//4x9i8eXNg5yw7InHBj3u1FHOOFH1O1TwkfJufwo0lLpnSZ0KVG9pZCVUEx/vPrRHqLpJZMe2sHScjQSYCJOtdJCkCedIQ5GqaoAiwuUroq6z+Fs2A53OVBJ2sK9GKIt2WBSJDUEmJKVQjhZjBsu/rMRMnTkRnZyfGjBmDuro6O6dFLIWgrFxbUWSQhhUBo3RfOA8IdtArRFwvJH/QoX+oxPlenKK3A55GT6Y97igjL5BCrG8zmFXuSVXyYYdiv7iFmGS/cHMJUI+jQHXeSqIRXIIYRO6IeDqP52wRaSgwnavE71S8vN+6KK8ASKY3c0nKBTNnzsQ///lP3HrrrRgxYkQotgcG8Qfif2aV24IP7iqCKKlOIiMCrzBfWeFGQF28UWXtaHSTTphZWXvc7TKdEcTj5uVQhJis2x8IVOSv+n3Q9pwD6XqMtOVkQiHmMFFEcImQWDFfCS3mCBiEAgNy0V30SwVBOF5tCyUURS6JuN99JZekv+G1117D6tWrMX78+MDOWXauLWEiy55+UYd8jSEO2Z9Bts1rEJXtL7rultcfSuXiovtLHMFmGG0sK6EeNAoiGBWZUIuVP5Gn2D7V/nbJPn5thVZivcpdWUKMpqHA2rlKAPVUvDLItBUOrsOYQvJ5xJCUWqdBTMUr3H5ey75aqoVi7NixSKWCtfLLjkgoIrX7yFzQSp7w0kXqFNtNnuR4Gx+CvKn7zqebr1RuDOMn2z2QuzopSYp74iRC31NS4e3ENajwLsC0EgCOya9EMUfrtHLR3bFNJ7rzqsCyYo78GE4+AeonQU5UVu647bbbcNVVV2HVqlXYsWNH4ZMAEpSVa6saKbQDtnlvDOolUrmx/Li3TJITfbvFdGG+Xu1l+wQCtFTSkIe4+rTAMojYmoCXe6skkH3PlCxUJCI7poa0V1mzCq0kFc8nHIrPQSQoApYQLZurxCjTHXCGAnOIvoqwXUp4nEBMLJcCEGQIsEWuAw3a7Q3ken2JadOmAYBrvpNsNouKigrlpFo6lBWRSFHbBXQW8DGIcVhWMqIUFYGV474XaeiE94KUfHIuDXgiooBEGzHRVfwShl/rxbfLQmad8O9cN80sHYhrFG0oBAEzrSQWF9qHOkFRJrp75pSA3J+M6Gk7ngNCy8AIUNcZnQ/FgFxkSYmCGEP4x8qVKwMR2CnKlkgckVvVUP+Jg4jc6nOrxItMVOBk4dM66fA4PUV7lZlI39/AtR/x8auit8C205kKxWdFCUgTwRWpzjgSFHWiO89015aXF/2jREB/h7w6sFexx30AGURQZWCRZBDsANybeOSRR/D1r38dkyZNCvzcZaeRyCbIAWD9oWvZOnfDeD098eMpZE91fgdOI61E5ZZS6SUmKCJhUWbwKAZXKiQL0NkBASokRxzb6NMpf1LNsPwKFZJJ9/WNoHJftbP3ssREHuElCyuWaSWAPfkVnfgKkIvuxjkl1HLgCYrcghB6CRfoZYtopyIZHxnx3DKtLn0Jif0Cjz32GA488ECcfPLJuP322/Hee+8Fdu6yIxIKbeSWbLuA+OOI8U31oy+ETAqO4PJLDsUkJQYAjZVHI7f8gBOKDroaTZn2uGUh+RmfZFnrsn10Wzt5T8V3sc7dZ3voMVV5rYQI7DrRXRAqvXdX9Bbg/P1TNxQX3mlb3UMR30fPwaf4lZzHDgrYB/HKK6/gjDPOQGNjIyoqKvDss8/a+/bu3Ysf//jHOPbYY1FTU4PGxkbMnj0bH3/8sfacCxYsQEVFhWMZO3asZ19WrlyJTz75BN/73vewbt06NDU14fDDD8dVV12FV155xVUF2A/KmkhKAi8ro1gyMRpbTawSsc9P5JdxB5xQhbfSgdEAfeITFwO+GNTFexmo9dXOXgU4OcnIhINbJe1wWCWAPBSYZ7rT98qSKdQ60f2WVWQiW0Qbryl7g86S9wFBwt7hv/5+/x0dHRg/frx0PqVkMok33ngDP/3pT/HGG2/g6aefxoYNG/D1r3/d87xHH300PvnkE3t59dVXjfpzwAEH4Nxzz8WTTz6J7du34+c//zlSqRRmzZqF4cOHY/bs2fjP//xPdHT4q8JZVhpJHCnshkjkymUEm1QBpjqJ8CULPUMmulOtQ6aFmGgmOsHeSHg30UvocV5QtPHzv9oDubVnELlFM9yFr59Gb/U6RIiurjAj/f46Je+rSTsqXMvAtZJa2LqSSnQH8kUbhXuLRm+lcvK7Hb0lri9Ed0A9uOtKv1MrSia4y8KEuRWku3Y/AA+TjUajjlkNBaZPn47p06dLzzFo0CAsX77cse3BBx/EiSeeiC1btuDggw9WXr+qqgoNDQ0F9DyPSCSCadOmYdq0afjFL36BtWvX4ve//z1uvvlmvPfee/jpT39qfK7QIvFCneK9DCoTXfZk59cy4e2L1ktMdI8SFHNUubRyORK8WGKxiYm6Io+BJTzqqgIDahdZJ9knq8tFIbNKAKl7C4DUvSW2u7ZVks+I19ziWgl3R/FXcaxY+Hllugsk7Rh0GlgQSOXSHb0W8bkddNBBGDRokL20tLQE0o/du3ejoqIC9fX12nbvv/8+GhsbMWbMGMyaNQtbtmwp+toTJ07ETTfdhP/5n//Btdde6+vYsrJIlKjtAlDl/CPT6BUOlVWivYbkXH4tk4JKrniVkfea91UcB826T4hLKsZxW/COu7OWeRhwn1olFPz7F99Zp2SbQC1pU61oQ0GtElI2JVNrfV6xeFKaU2IavWXPU8KjsgDnDIog/fWYlMoBSiaURChpsbZd/dQq2bp1KxKJ/P9AZo34RWdnJ3784x9j5syZjnNzNDU1YcmSJTjiiCPwySef4MYbb8SXv/xlvP3228ry8GKueC9UVFTg7rvvxsCB3hFsFGVLJNKkRPpD5uOraRiwysUlWwfMQ4Nl7QvKLVHNSWKCgGsbic/TMHpN5t4CzMiERy3ltwd4Tx3wDvfl2yiZAO7PgubjUAImDxc8FFjl3hLvqXuLbhPnMwrj5XWwTAlFkAhdp6A6Sq6dsJj6U3mSRCKhHez9Yu/evfj2t7+NbDaLhx56SNuWusrGjRuHpqYmjBo1Ck8++SQuuugi6TFvvvmmUT8KzS8pWyLxBfrgLt7rrJIgyIRCl1/iWy/h+7xg8GcptrI2tbgUSCNqh27LkhM5mfiZt8R4gKKiu652GLdGVOG+XEvj+wHrq+N5JZSAc8mctOAlzSnRJSdabfMTXtkVgVXVfAH35FPit6g7hruzxDahjahcaBLw0vlBIYMoBhicsxR5JIJENm/ejJUrV/omqPr6enz+85/HBx98oGzz0ksvFdtNLfYZjWThwoU45JBDUF1djaamJvzlL3/pnQv7yfWgD7e6kEaeswK4B1GvSC7VdR3gegnXTPyEAPdCpVWik4gcCV0YsCxfRCw6yCJvHDkkIvSXR2xxCAtE59ZUkYhsnbvBOiTnltX1IvW3ePQWkJ+nxNqurr0FMFeS10BPLQdOEDxMWNZGBrKvqyYfUaazHDv34TlCBIm8//77+OMf/4ghQ4b4Pkd7ezs2btyIkSNHlqCHZtgnLJLf/va3mDdvHhYtWoSmpibcd999mDp1KjZs2IDhw4f3XkeEe0tllQDmlolsG7dMAtFLvMqiBJ1/wuD1VOsTfH4Sk9IpKreWDMqZEb3A3VqqMF8O4d4SrzIXF9fuxLqY5524twC4oreEG5fX3sq3zbl5K5EvmULvS4D2ic5ZItqZzqqoEvTp74Kci4Ys0/Ug3ZJJxNBj8MPsNJiOl6K9vd1hKWzatAnr16/H4MGDMXLkSHzrW9/CG2+8geeeew7d3d1obW0FAAwePBiRiEWip512Gs466yzMnTsXAPCjH/0IZ5xxBkaNGoWPP/4Y8+fPR2VlJWbOnGncr7Vr1+LJJ5/Eli1bkMk43cJPP/20r3sE9hGL5J577sEll1yCCy64AEcddRQWLVqEeDyORx55pLQXpr8rLylBZTFwy8RvAqJpJJf2P+VjPhEpCvzDqgZPWt3WB2iUVbGRO6UqUS+1QDjB0GgtVTthDXEipm418UpySmSRWtR1l5F8hnS/K/mPWiV8G+AkB1lYL9+uIxHWNlkTcekjquz9/oq1a9diwoQJmDBhAgBL9J4wYQJuuOEG/POf/8Tvf/97fPTRRzjuuOMwcuRIe3nttdfsc2zcuBHbt2+31z/66CPMnDkTRxxxBL797W9jyJAheP311zFs2DCjPi1duhQnn3wy3nvvPTzzzDPYu3cv3nnnHaxcuRKDBg0q6D77vUWSyWSwbt06XHfddfa2AQMGYPLkyVi9erX0mHQ6jXQ674soqDSyKqBJZZVQ8Id+3sZLNynUMtFGj6lqbHmFm2kqAQfxQEhdRrKs5twkTkB+LhlhldDoJGu/2zLpk4FGFanFrRSaT0ItE7FNgH5F4oFG1MKyj8vXKkvmCjMCYHpIzBGgwGdOBODMKVGBBgBQy0Q1wRUPDRbbVPP+VOddbEni2iyVPlJKTJo0CdlsVrlft0/gww8/dKwvXbq0qD7deuutuPfeezFnzhzU1dXh/vvvx+jRo/Fv//ZvBbvH+r1Fsn37dnR3d2PEiBGO7SNGjLDNQI6WlhZHjPdBBx0EwIcp7EefULXjl+J/LhPLpFqzv2DLRGadcP1Ets0nfFob9jFMJ+HQWRKyTGS6D3A+qRflGvHSR0xzSeg+mYYiq70lyFeSU0LnKVFluXMdhX6mjpwSGSgBCEtCNieJrJSKjES4NZIDt0a4WwvIW1dBRHN1Gs7Zvi/rMQIbN27EjBkzAFhJiR0dHaioqMCVV16JX/3qVwWds98TSSG47rrrsHv3bnvZunVrMCfm7i2xrhLWCyETnatLR2i+yqiodA4deXhoI71QuZcO/Ck2GBYKX4OQKmJLF/YLOEkkxRZZm3a2LkR3eg2Fe8u6hDx5T/Ze5d7qorkeNMJKgBIB4ExClC2irY5ENNaIu/9OKyWEOQ444ADs2WP9cD73uc/h7bffBgDs2rULyWRhc/v0e9fW0KFDUVlZiW3btjm2b9u2TVkiQFWuwAjcjSTLH5FtUwnrJm4uwNzVpSu/YuzmAtzztZu07TukkJ+wiYYCA/7nbudWjXiKV0LU2PKCV1SW7PtIIc/dqlwSig6o3VuQhwFb760wYMB6kpfNUULdW0rQxETx+9OF/gLyBy0NiVBrBHBWgKbTBgfptkwhjm4DayOzHzx7n3LKKVi+fDmOPfZYnHPOOfjhD3+IlStXYvny5a7JrkzR7z+VSCSCE044AStWrLC39fT0YMWKFWhubi74vIElN6lEeC/LpBBXl2wft2K4ZWJkneiKN/YtiYgwYBlk5eU56HZtxd/OiDtiS1egUQURNkzXrYurIcuKt918cIcCy9xbgDIMmGoKsil4rdM43VtKq0RAvDcpI0/b8fwRCTKVEdsaoSK7zK1l9X3f0Uz6Ax588EF897vfBQD8+7//O+bNm4dt27bh7LPPxsMPP1zQOfu9RQJYkQ7nnXceJk6ciBNPPBH33XcfOjo6cMEFF/g+F6166gvi6UuI7F5WCT0GkGvbxYjwfkKDTUq4APBFGuJ+gnQZO57gq5BBHJFaXiIlhjhStlXCs91NofwN0ByS3oSwTGi5FBl4pju1aBVhwIDIbneWS3F3IW4HNCjBRXZqhdM+8m30WEoiEmuE9kdmjVC3lvWdF17+vBwxePBg+/2AAQN819WSYZ8gku985zv43//9X9xwww1obW3Fcccdh2XLlrkEeC90lsqfSiO8vMjBr6uLr3My4W2LJpOA4KWZyBL8PEKsqXvLL3rdl14IEVE3l6rGGyB3b5Esd+tUeVIwKZcCQO/e4u4r8SAlc2upCEQcpyERao0IyKyRoKPx0ogga5TZvu+T1gsvvIDKykpMnTrVsf3FF19Ed3e3slqxDv3etSUwd+5cbN68Gel0GmvWrEFTU1Nfd0mdZ+LltirW1aXLNeHtfLu6DKA6h6ng7tNdxN1bYpAR7hgvkqD7qQun12o3+SVwVYQXd29RkCx3Gu3mlU9C3V5a95YAj8Ti272ithQkYl+TuLRks2NygkkjEgruPnHttdeiu7vbtb2np6dg62SfIZJeQ1AuDRMy8SIUrn/oSELVDnATTzGEQo8T5y1hxBbXLVSDvwjPlG0vKXTz/3hWHpAEOlDS4VoJ3Q4ow4ABp04CwEgnsfZpyFWmkfBy8Rw8aku8sv+HbBZEOtkUn1o56O81lSMv72XfJ633338fRx11lGv72LFjtfW6dCgrIkkjYj/teEbq6ECtD132u4xMCiUU1Tn8iPDFEorf/1CxBNOu9rxyq0SAx/1T5HMRyFNtMb8DMejTAd8IbciTCH1Pzis7F810F688DBhO8k0zNxCfNVEWAZXMDeBaq4S+AvqwX36MeK9wadG8EcCpjWRc3/f+Mbj3JgYNGoR//OMfru0ffPABamoUERAeKCsi4ZBG6hQLGZkUQyiqY1QkoWsHyIVcE0Lh+1WCcFDiu2QgVQ36JSt3ooNpuXiXW8tnlQVZUiJ3bymy57l14RU2azQg67LWVe1lxFMtn2uEurS4NgLkrRGaZBnCH77xjW/giiuuwMaNG+1tH3zwAa666iqjaX5lKGsisaF58jUCH8dkorHsD2dCKDp3V5DWCb22bKGQXceP9VHAxItOv7+7aJ8XmfQq2RTkHvWwSkzcW2Q/10lk5OGlk7ggs0rEe9nvW+YKY9aNSmDP91FtjQDu8OZCUao52/sj7rjjDtTU1GDs2LEYPXo0Ro8ejSOPPBJDhgzBXXfdVdA594moraCQQgxplQ9YVR9JBV6Li4cD8zm2AecfyyspkYcL8+gYVWKiV1QXWFsBU7eMSaY93V6IpUwS7IJAn1gsRhBfrviy2+AZgk3DyUX0lisMWF13i89PIkCTOgVR0/1dNUAVraXVQV5BtssgsUh0Li1Ta6Q/TXS1L2HQoEF47bXXsHz5cvzP//wPYrEYxo0bh1NOOaXgc5YVkQiIuScCd2vJcktkhAJ4h02q8k/ocXTwYPkEAMxyTgS8SEXnsqo1aFMEMu35fJJkMm5nbgP5vBIgTxg0692VxW7nIJAIJklNLy10ArsW1OpIsfcSMpFFe4nPOpV7L3KaAAehZNrjSFZnEMt9VpQo0qSQY5IRCoUoLY8aIN5R4HTGnEQULi0VdNZIkMgYhv/uRVdJ+9Eb+Oijj3DggQdiypQpmDJlimPf66+/jpNOOsn3OUPXFnVrBRGxpfq910Ht8jLRUWQuL9pW9t5LO5G51uhxdFH1zcQqkXkDihwXdE+jaUTtpc+gtfBUtVIMYKLNsGurPoeCSrLz34LKrQXJdolLi/dTFalF+yysEdGmP1smr7zyCs444ww0NjaioqICzz77rGN/NpvFDTfcgJEjRyIWi2Hy5Ml4//33Pc9bzER/U6ZMwc6dO13b//znP2PatGnG56EoKyJJISYXbGVuLWFByBLnvKAbv1SEAvgnFFkJCtkx3B2lKhbp5U7yEvz5Nn6sCQogc9NKBYHOz14SCDLxEOT5/CVpyT77jHLBXd0Dw89IJ7qryqnAbY3wUigyiATEfVGf6OjowPjx47Fw4ULp/jvuuAMPPPAAFi1ahDVr1qCmpgZTp05FZ6f6SURM9Dd//ny88cYbGD9+PKZOnYpPP/3UqE8nnXQSpkyZYhduBCzCO/300zF//nx/N5hDWRGJgLKukl/oRGOv+Po6tlDIrBQZodD2snb0vSzvxItU+MLbCagSJL1cGCaiO7EYfbugfCBwN6cLBcyJA+RFdx69JQsDZnW3dIK7eO9HcPd0SckehIhLy762JGeEg5dDKbSNKYSF4y22+/udTJ8+HT/72c9w1llnufZls1ncd999uP766/GNb3wD48aNw29+8xt8/PHHLsuFotiJ/n7961/j4IMPxhlnnIF0Oo2XXnoJM2bMwE033YQrr7zS1/0JlCWR2Cg2WssEsph6GbxIRYAO6irrBHBbDipCoftNLRKTSDEBFeEVAWFVOsulmw8mdtRXMTkkZhcy3SjZryAeXpJeBvZQxD+bQlxB0sG/iO/Tj8hu94G5tfoSbW1tjoVOpGeKTZs2obW1FZMnT7a3DRo0CE1NTcpJ+8REf/QYr4n+OAYMGIClS5di4MCB+OpXv4qvf/3raGlpwQ9/+EPf9yBQVmJ7KhlHBQ+KCdqtpYOMTLxqTtEaXoC7SF47nIK8Soznx6jEdVMXlMoKkVkj9JycKHWWToDon24tShY+CmbyqDaF4A5YZBkjwQmqGRMB78gtG7yOmwqMaPyI7BS94dZKIoYqA3LqgjWjoZgsT2D+/PlYsGCBr2uKifn8TNqnm+jv73//u/Jaf/vb31zbFixYgJkzZ+Lcc8/FKaecYrcZN26cr/sAyoxIgCLcWqo2qil5TcHJxSvqSxcGXAihAPJoKz5QqCKyZJqJuHYvj900gmvfAw39pVFccFoe9LurIa9iG/k90UrAgHreFpPIrUh3gVFbgG+3Fu1XqaO1CsXWrVuRSOTJv+D5j3oJxx13HCoqKhxT+4r1X/7yl/jVr36FbDaLiooKaR0uL5QdkdhQubUKIQVViG8hoL9HXV6KjihopV8TQqHbBLxCeVWWhKrGGC+nUSLIyKR/WiMyGOaRyL6bNJwVgZkLkpaUt9YjSvIIBBrrw3M6X3gHBeT1kTg6A6jI24kYKg1+J925ayUSCQeRFAIxMd+2bdscc6Vv27YNxx13nPSYQib6Ayw3WilRvkQC6JMQNZEwShSQsa0lHxmpmBKKSf6JaAsU7lJSJR/K/pN1ivdFEIustDwlE18k0huamS8oiEWUi++AMzERcBJKDnyekaTkMysYPClRbGOQRWvJYJKp3tf6SFAYPXo0GhoasGLFCps42trasGbNGlx++eXSY+hEf2eeeSaA/ER/c+fOVV5r1KhRQXffgfIW2/sDePSWKjyYC/a8nSoUmJc54RFhXtFZHLK29Hz0Wr1sjVBYhQf7syVSYBSXDpKHEq9IN1n2eKDQWLb8+zG5dn/OGZGhvb0d69evx/r16wFYlsH69euxZcsWVFRU4IorrsDPfvYz/P73v8dbb72F2bNno7Gx0SYJADjttNPw4IMP2uvz5s3Df/zHf+DRRx/Fe++9h8svv9xzor/XX3/duM/JZBLvvPOOr/vsb49gJUWmcyAq+P+qt2fCM4VMcAfcVgpvp7M8VFaKgGy2OxVkrgs6LnASkZGjSoQvkeAOuAeiUoYUFwemk3AI95Usw13sJ6VSBOgkVybbSw3ZvCMqlLJIYxoRVBo85XRjr6/zrl27Fqeeeqq9Pm/ePADAeeedhyVLluCaa65BR0cHLr30UuzatQtf+tKXsGzZMlRX5xl448aN2L59u71eyER///qv/4oxY8bg4osvxumnny6t8vvuu+/isccew+LFi3H77bfj6KOPNr7PsiKSfRZepCJze5lEeQmoiMUE/FwyV5foV5Stm6C2l0tSGJeCNwDVqoyhcGeJ85iQLIvcyp/CHaUV22eDE/YNTJo0ySFwc1RUVOCmm27CTTfdpGzz4YcfurbNnTtX68riePfdd/HQQw/h+uuvx7/8y7/g85//PBobG1FdXY3PPvsMf//739He3o6zzjoLL774Io499ljjcwNlSCTaxDOZXlEqi8VkQJBdW0YqJoQCqPWRIDxAKq1ERSK1bDs/R4lR8hwSTyRQkHuL11XjkVsK0OKN+yqsYpOl+d5SiGGAwR+hB/4jmvoDBg4ciB/84Af4wQ9+gLVr1+LVV1/F5s2bkUqlMH78eFx55ZU49dRTHfO5+0HZEUmfwu9AKWtPyYWTiqnbC3CTiuz8fvoFuAczTiIqyCyUEhWA3PdgEMlFwUKAAXcuCQC7CrAMqf2AdEKoMXHiREycODHQc5YdkURqk/7KYdA8jUIR5JO2KmSXR3OprBTRphh9RED2FMxJQTaDpMwaAdsnAc2JEAgs+ohC930HEeJdLExCgHPguSQyZBDt0/wbQVomYnsMyf0mamt/QlkRSaR6r14qKza5kKNQAvFKUpSdXwx+XlaKrA1QXPkSVZSZbD8nEU4uFDl9RJSQF4jHS5j/EBRc2d8FurK8IEKAZRBhwgwZRQ6Jarvn9YtADElXFn0QE1X5heXa8iaofdW1VWr0afjvIYccgoqKCsdy2223Odr87W9/w5e//GVUV1fjoIMOwh133NFHvfUJUxLRzXOta8PbycJ3eTix7HhV+LGfxevctH+AmkRoPbASwJd/vb9G86nAB3RF/3vlaZ72RRO8oLKCvKLHSppIGaIg9LlFctNNN+GSSy6x1+vq8o+vbW1tmDJlCiZPnoxFixbhrbfewoUXXoj6+npceumlwXQgCNeV7JxeCCKfQmW5mFgquj6oLCCTPutqaencXBzMGvFyzxSD0lf+pSgolMsNhbXhFzySq1So6nAmJcY9XFR560T+o4sgEygpppJxVFR5ny+bDC0SGfqcSOrq6pSp/Y8//jgymQweeeQRRCIRHH300Vi/fj3uueee4IjEBH7IxuvPbUogJiGy3A0nc2N5kYrsXIWI4xT8M1DNaS+r0+VhjVC3Vkn0EYr+oIcUAolWYm2O9l64r5iKV4JId8aoTAoARJBWEoYgm/YwMqPP0eeZ7bfddhuGDBmCCRMm4M4770RXVz5vYPXq1TjllFMQieR/dFOnTsWGDRvw2WefKc+ZTqddZZ6l0JVAl7mOdDDJCi90fpJCjpG5wHQZ7H5dWbLzcfeVyoUGON1ZnESKsEZKLhqbfDe9GMbcJwgy18YAXhZTdZgLY4TTTz8du3fvttdvu+027Nq1y17fsWMHjjrqqILO3acWyQ9+8AMcf/zxGDx4MF577TVcd911+OSTT3DPPfcAsMosjx492nGMyN5sbW3FAQccID1vS0sLbrzxRv3Fa7vktZV0grvKMimWQIICP5fKhSWzVoqBbnphDpkVAihJhMLLGhEkoqqzRYXdoiByN2Qohau0v8JUaO+E/f3GOzK+KgALxJAqqb5jVb3w7le2c2DJ+lBq/OEPf3DMm3Lrrbfi29/+Nurr6wEAXV1d2LBhQ0HnDpxIrr32Wtx+++3aNu+99x7Gjh1rlwsArBr4kUgE//Zv/4aWlpaiyjJfd911jnO3tbU55g8wDgGOwq0X+B18/QyyhV6DQpVnAnhrI8VAdT9e8404rEI3iQhrxJRE+LaS1dsS98HDpn09rZcokqu/gLm3uE7iF3yulBD+wDPsdRn3fhE4kVx11VU4//zztW3GjBkj3d7U1ISuri58+OGHOOKII9DQ0CAtlwxAWzI5Go1KiSgWT6Kiqs5dX6kUT5GFWCFBWAj8HDpiEfDSAkytJqVwztZdk2Ll3Zkyd5aORErqyjIJB6dzgZQ7NLqIDDqXlS4EuBS5JHs74sAAg3N29HLJnn0EgRPJsGHDMGzYsIKOXb9+PQYMGIDhw4cDAJqbm/Hv//7v2Lt3LwYOtEzK5cuX44gjjlC6tXyDxvsLQqEDiMwq8YJfK0RHIIWQCx3UTLLXiyUKkzYa8gDkVgjQhyRSEsgsEL7NMItdpy/3xRxLdA55XkGhBg73FuAtuIsQYJVFKUgoTE40h0ix4NuCQJ9pJKtXr8aaNWtw6qmnoq6uDqtXr8aVV16Jc8891yaJf/mXf8GNN96Iiy66CD/+8Y/x9ttv4/7778e9997bV93Ww68VYvoE7xde5FHM+b2OVc6k6H6S4zqIikCAwknExL0ldXUGZaVKo339hAAbEouBJSAre2ISxRXpziDekUEV1URU7+m6ok9+dJK83qWO3AptBDNks1mcf/75trems7MTl112mV0JuJB55wX6jEii0SiWLl2KBQsWIJ1OY/To0bjyyisd2sagQYPw4osvYs6cOTjhhBMwdOhQ3HDDDYGE/kZqk8ggnhfc+cDh1yopVHA22Q74S9KT+emDjiTS9UdRsVcmoPOIrKAIZN+CjCw05KcqkKlqYwhhBfgKETYR3ImVUqhOUnILpKMKqDAYDjv6PGOiYMyePdthgZx77rnSNoWgzz6V448/3miylXHjxuFPf/pTYNcVA5VyHgrVkyivXcW3y2BKIrJtxYTGex3rJQibXFtT3l1GGIA8jFdW7sQveYgn7bTiy/AtugdljbjKpBSBWO58OoIwcGn5TT4sKlmRu7mYeyt/jfz36zUvyb4muB9yyCHYvHmza/v3vvc9LFy40LV9yZIlrgmqotEoOjuL/yEtWbKk6HOosO/SawGIIYW0bHRX/eFlYmsh2d1AcARiOj+HbtpYPyQVEGEA6hpZplFYgNw9I9unIhURAhyPJ0tbSl5JRkIT0bm3YqStxzXoK4fie+YlSCJ+RMAAc0gi3Rmg0mllUGvIS3D3auMLnTAbDX3e/1//+ld0d+ez4d9++238n//zf3DOOecoj0kkEo4w3KB0jH/84x8YPXp0YOejKCsikYLmk4g/Px0E/BZyLJREfGgLRghwQigVYQDmVgagzkTXWRyFljOPIq0kEy2CsCJ6I5fE0D2q+u5MqgJ4uhGLKNhoopN4Ce6iTawfqyQ88Oi2227DoYceiq985SvKYyoqKrRRqYXi8MMPxyeffGIHM33nO9/BAw88oJ1Z0RRlSyRiALRFVt0AYkImheZRiGu7jtP/OXSDO4cuZ8b0PH6sDL+EYUoWOt+9zCVCyUTl3opUZ4KbbleVS6I0Pvjcx3Rbwr1JgP5exD76+2O/xXg8qcizSbreG7myVAQiE9mFe6sIncSP66u3wKtlqFIOKDKZDB577DHMmzdPaxW0t7dj1KhR6OnpwfHHH49bb73V17S3KvC8kRdeeAEtLS1FnxcoMyKJIYUBSAJxNkOel1UCqMnET0iviRXiQ6g2RbFkIWBKGoUQRrE1oGJBDDaqagcCXmHSMitE+oDCQ34N+iv0EVMRXdLOz2fsu44Zj+JSEYVCJzGF83sO0DVpal3l2tEEZwCYP38+FixYoD302Wefxa5du7R5dkcccQQeeeQRjBs3Drt378Zdd92Fk08+Ge+88w4OPPBAw072PsqKSGRwhH7yP72MTExQCIkUSSBBVsctVssoJMy0UFFXFsUTQ8omE5lV4qtUCv3O/bqrtC4uXVa7z5Bf8eqqvOz+TfHPOZDKv7JBmJKJRwFHVMr3eUVoid+jdo6hEmHr1q1IJPLfk0kljocffhjTp09HY2Ojsk1zczOam5vt9ZNPPhlHHnkkfvnLX+Lmm28uqs/7ZR5JX8MVvSWzSvh7EwREIn6FbL/wmhzKj3uKE4eMNLwGLNMBjQ8u9Dgu2srIxHEsEdwDzyXhxxqnjjASEdFa/NxcaI+y1xwi1RlJFJw8Kk5sN3I1FqKPMELx0kl4P1Ul5QNBBwCTiiG5LiUSCQeReGHz5s344x//iKefftpXtwYOHIgJEybggw8+8HWcDF55JAJ++wiUGZFEkUFl7k9Dn0ilVkkhZBIAiRRLIIXMHlgMaVjHO9vJSEFHFF4TGdHIHNl5BIHwXANKJvnjDUOBiyERbtkq9TePgYh2U5CHTB+hbSBpB/f3ppqvncORjCi7Bx2haCyRQvJJaHTZvpbRvnjxYgwfPhwzZszwdVx3dzfeeustnH766UX34bzzznOsy/JICkVZEYkMDrFV5yPXDSwB1JgCzBL2KPyShpffOyji4Os6ovCyRJKIK48XBMMJhE6aJMiEWyXUveUS3AuN3BLft+p3UsicVio9gQ7CCqHdNUUxko7Pkn6PsmTEkk14VaROAjj7lkEA7pk9gNEsugV8JD09PVi8eDHOO+88VFU5x5fZs2fjc5/7nC1633TTTTjppJNw2GGHYdeuXbjzzjuxefNmXHzxxf4vzLB48eKiz6FCWRKJPWAS0d3TKgHcA4WfGlk+ScRvLgZFMaGdpSAO1YDk1c8U4spjOcHQ+caTueO8yASAO5/ES3Cn4E/UsqfzQi0blTXCBXeZPkL2yyK26OfuK4cE8O/SkkRtceh0EqBv5nAPEn/84x+xZcsWXHjhha59W7ZswYAB+WmhPvvsM1xyySX2NBknnHACXnvttYLnCektlBWRVCOFvZLBRGuVyAYCr8iZEpCIjkD6E3Hwfaq+mbpWZKAEI8gigozDQuFk4uy/IhRYN72AzD3vFRJOfzvi4cTEKhFd80xOZX1zTGLm/I3pvlNZ6G/JStH41EkAuUXbX8KATTBlyhRlyfZVq1Y51u+9997+W0tQg7IiEo4YCwVWRnD5earsJRIpxE0lUCx56CwO3i8ZYRTjMkkjYl+DE4rol7BOKJnIrBKte0tAfJ9+p92V/X4omVg34AYnEZk1UkPauSK18m9ddcwMI7ZKPoVxDsXMT0J/rzH0FN+ZDLRWkaNdCBfKjkjoD5A+lUqtEr9kUmRhxGJJpD+QR1Rjnciu6QdJxKXEJCMUHZlYfctbJdJyKWIgLyTklx7Hfzf0N6V6qKYkogN1a9UiTyrk4SWGpOO7jyNZPLl3sFcB+hOj1pHMvRWAThKi/6CsiCSOFDISR6227hInE8Dc1VVghBbtF0UQGeNeobmmbisT8jCJ5ioUMlJJEdJQkYnVX7lVYsOPTiIgSIEOrqoHEZ2Yz0lEZY1wtxacx0Vqk47fj9vl6JyWWKDQkjSuw9KSvgko9JJCJ7oKBO0wE9v3x+LTAaCsiEQGPpDY7i2vCC4veDxtebm0TEikVLpHoeThdQ3fwq4BVFaKjExEH2VWCSApm6MDHyRlDxfcKtGRiWPKYfaqAnVrKfQR+hux6lJJIgMV/hrpPCQlQLwjA9TAaKIrjn0tDHh/RdkSiRhsle4tikAK+ZkXliuURIq1QPwQiMzycJ7L2ZdSzyMiCIVqKJxMxDbAbZVQrcwF2YCu00x0hMH3qa5VS9rIrJEazTFwPpioCJ5/R9rQ306oI7ZUzwfUKlFEbwUxj3sGA7wbhigpyopIYkiiyhZqWaKaLMuZWiV+yCSgUF/RZ44g8j1ULixTApENOhGHL750+QgZRF26h1jXkQnti8oqsVHbBfvv4UcnEYOl63xQ6yZiG3/vpSFEkbdKiD4ifl8x9n3rEkV9R2yZWik6F1d/0knaAaMiwgGW0d+fUFZE0isoNtHKI0/ExKddrBurUAJRkYdX5nopICMT8R7I34csFBwwdG9RiME/xbbRaC1OJrrz8Kgtbo3I3Fosf8TalXZ8FzKLsigU6q3UFXYMsc8hJBJ4hIEWapUYwo81UioSkVkhJi4sTiBe59b1y8vXTd1TOgjrBHCTCb0uvZ4sQdUYYkDnT+iqKD+TCgmmUVt1cLu1NPoIh+x78wz9LUQvEVaJB3moEhN7RQdJA0YJ8sHLfPsFyp5ICkpOU8GgLLwuUktnjfQliZgSiOycpiHA8rpY8oRCiiRiiCBtu7sAN5nQfnAXl9XGoBigalCXDSzUvcWtEt25AHnUlswaoV1mYb8qfSTGvk+OgiO2/EKik8ShTkzk312I/oeyJxLfUFklBbi0/GojOvQmiXgRiC4E2GuqXNNZDWUWiopMrHVnuRV5+XlLdM+fLw7fOgnNXOfRWl4WrYxEdKiFy62l0kf4e/EwEKjbkVsrJq4rHzqJyYyJIfoGIZEUghK4uPzAq2y7n8RC6/jiSIQTCO0P7WsQGdNJqOtvccjIREB2jqIGKDFotpN1mRtIDJqmob8qa0RWWytn/VJ9RMDrM+O/IW3VX3Ffack2FWTurQJ1EisfLOCS8u0wm9gkdG1JUbK4uVtuuQUnn3wy4vE46uvrpW22bNmCGTNmIB6PY/jw4bj66qvR1eV0Ba1atQrHH388otEoDjvsMCxZsqRUXS4chvOtFzrLoZ/Q2UKFdb8kEkHGRSIxpBBDClGk7SWGpL3EkfJcaHuRlW2dxzq3iD4SC7WGRP8iOYHZSWiZXB+SwT6F80GdQ2VdVJNFBa9Bll6bubVk+gf/jmUWpGwdQOHzj4QoC5TMIslkMjjnnHPQ3NyMhx9+2LW/u7sbM2bMQENDA1577TV88sknmD17NgYOHIhbb70VALBp0ybMmDEDl112GR5//HGsWLECF198MUaOHImpU6eWqut9gkLmEQH0JUdMhW7VOXUkIt5T4hFPwDF7gPIfAmyii3hB6CbObXprhHXChsO9JYMsn0S4t7hVYjoVAX2vskYk7blbK2YTqzlxFpTvU8icJP0taisJs/DfsNaWFCWzSG688UZceeWVOPbYY6X7X3zxRbz77rt47LHHcNxxx2H69Om4+eabsXDhQmQy1re1aNEijB49GnfffTeOPPJIzJ07F9/61rf2yeqYQaFQQdSPS4u+9yIRaoWIc+YtDGI5dGeMF6fV4bRUdBaKsJBklkkcSdsqKRlqINc1+KDvtS7O5eVpc0RtyUdBbnmqSqMYIbQwfGPBggX2FLdiGTt2rPaYp556CmPHjkV1dTWOPfZYvPDCC73U28LRZxrJ6tWrceyxx2LEiBH2tqlTp+Lyyy/HO++8gwkTJmD16tWYPHmy47ipU6fiiiuu0J47nU4jnc4PuG1tqrmx9y/hzitaSrddwDnZkdPHLiMRALYLy9qXypNQN3Gldfh4lKNPqpX56+ssFZl2YhdlLIQ84pBXOVChDk7rhFsl3BrREQ7PWqd5I4I8qERA3GMyfUQHX6G/FIU8z+iSE31gX/vPHn300fjjH/9or/PJrShee+01zJw5Ey0tLfja176GJ554AmeeeSbeeOMNHHPMMb3R3YLQZ0TS2trqIBEA9npra6u2TVtbG1KpFGIx+Q+qpaUFN954Ywl6vW+iUGvEeu8W8p1+dzWJUPLwU68p0ZGxy2Y4wkJNynxLUKgwaycnerm3AOfAWkwOiSAR2U9bVqBRUhaF6yP8VRax1Wuhv/0V7YBRdHEBrq2qqio0NDQYtb3//vsxbdo0XH311QCAm2++GcuXL8eDDz6IRYsW+b94L8GXa+vaa691mWl8+fvf/16qvhrjuuuuw+7du+1l69at2vau6q/wyGzmUSwBRHD5ToQrAbysER7WK0iEurJs91SORKo6kI/8ocsOtrD94riqDsuaiXdY563v3uVyeQlXl+grd3OpXFzOdkmHi65gUGsiJtlWq1houxhpy7PYAUnUVpenPlLMRGJamDwg7Icusba2NsdCPSAc77//PhobGzFmzBjMmjULW7ZsUbZVeWFWr14dWN9LAV8WyVVXXYXzzz9f22bMmDFG52poaMBf/vIXx7Zt27bZ+8Sr2EbbJBIJpTUCANFoFNGo+yk0hXju6VINX+4ML7RXOXzXmfZ4wZFbAmlECx7s/Lp4ZH501YArSIQSCIA8ydLBRDawcPGVrIsfaVwMhjUo2DrhUM2gmHsDwGc1YOHeEtaHSnjn4MIzd33JRHbAt1vLNEFUiT4mhZJluXfCrIx8LkT4oIMOcmyeP38+FixY4Gre1NSEJUuW4IgjjsAnn3yCG2+8EV/+8pfx9ttvo67OHe6n8sIIL01/hS8iGTZsGIYNGxbIhZubm3HLLbfg008/xfDhwwEAy5cvRyKRsOcnbm5udglNy5cvR3Nzc1HX9jVNJy8lr7I+CihAl+mMKJMSU4iT8h6xwsRRQ+i0FV24bN4akZAIJxAvIhHbNZE8DkIR7RihpBHVDox8IJINrK76W7VJf9FbKtLQbReIKbZzKNxaFFxol+1TwZem5YX+FJ0VELZu3YpEImGvyx5cAWD69On2+3HjxqGpqQmjRo3Ck08+iYsuuqjk/ewtlEwj2bJlC3bu3IktW7agu7sb69evBwAcdthhqK2txZQpU3DUUUfhX//1X3HHHXegtbUV119/PebMmWN/KZdddhkefPBBXHPNNbjwwguxcuVKPPnkk3j++eeL7h8fLIxcS14uLBmZ+LBKksl4wWHAfsD1EQ7ZvCHc/cNJRMAmEU4gdADVGVS60uO51yrorRPx3coEetUASstv0Cd630/AVHTnVonorwrUpUXPR60RjVsLkOeP5E+vn1KAo9TzkPhFfyqRkkgkHERiivr6enz+85/HBx98IN2v8sKYaix9hZKF/95www2YMGEC5s+fj/b2dkyYMAETJkzA2rVrAQCVlZV47rnnUFlZiebmZpx77rmYPXs2brrpJvsco0ePxvPPP4/ly5dj/PjxuPvuu/HrX/86sBwSVfSH7xpbFDKy0cy2Z+pKo32lJOjLuiohHFnQlDzEks4tO9l22bKTtJft7wSqduT1k3zIsFwzUeklFFQrcdxXPIlIdcbbJcnrXwnCEKTg9TVxEpFpIxRi5kTST1ktLVnSoZd1Sx8M+kvZdB2hp4MgGK/fJF2KQHt7OzZu3IiRI0dK9zc3N2PFihWObUF4YUqNklkkS5Ys8cxCHzVqlGeM9KRJk/Dmm28G2DMnqNDuGtQpAfj5Q3m4ubhVQl1c1Cqh7i0T0BBYP6VEAFmZlZT2eO7ScpXSoAQi1gH/c6ALa0RmqYBYJz50kwyxWAToQCU+C9eDBp2jhIPeFw8Fzp/YQkqyDZCH/8qsEaaTULcWtaZkLkkeRMHhO1TaS/uhCLiqSSpnFwM9wZ44QPzoRz/CGWecgVGjRuHjjz/G/PnzUVlZiZkzZwIAZs+ejc997nNoaWkBAPzwhz/EV77yFdx9992YMWMGli5dirVr1+JXv/pVX96GJ8qq1lYSMekNe7q1pFYGeS/LCeBkUqCLS6WVFCO6+yEoWVST9lhuhYht7eQ9Bb8FCVl4wQ+ZWHOTuF1XMlKhortxEEYU+XtS5ZDIrBNOIjprxD5G7tay3rtLo3ihJDqc6jv0oZvIrJHARfcOmI2G5hOdAgA++ugjzJw5Ezt27MCwYcPwpS99Ca+//rqtNW/ZsgUDBuQdQyeffDKeeOIJXH/99fjJT36Cww8/HM8++2y/ziEByoxIAOspRldh1nZrCWuEk4jsiVpFKuJYQSgaMtEJ715IIUaeot2WSAYRRJDxbaWYwuHSUpGIaaE/SiA+Cv2ZkAm/f9k6kB+E6e/ElVNCv3P6cxKWiLBKKJkA6vwSCkEiKmuEubUAdbVf1fa+mGzMBQMykbluAy/YWEIsXbpUu3/VqlWubeeccw7OOeecEvWoNCg7IhEQWQUCRk+cJm4ZUTacglonhmTi1yqhZJK/Rz1xpBEpKr+AJx26ICMRv9FbOhJJk201ejLhlYCTZF016VX+s5U8/dLvmLuxqFUCyOcn4ecQ7eg5DBCpzjgCNLTzrmvgaWXqYBKNFpVsA5Qu4BTixlMKhOh7lC2RCFC3ltYa8ePblw0UvUAm9j0VYHnIiEgHZYVYYY2oSITnl6jAXYMyQpG4FGVkYn12+Xuz5jJxW3ApRijiaZjOU6INxBD94VaJqh2FzKUlyxuxz9Hlco3qXJ06gunrrPauGiBZE0GmMmI7UWXQ7SsaHTDT2ExyTcoQZU8kgMIa0ZGIjlToH58TSonIREBnlQTp3pL60nVPrZREZAmKdF2WlEg/Ny/thFkmyUqnHuScv91plYh2IsxUJrprc0roeCysEurikt03vRfq0qIwcGtxfUQ2ZUCgUFldsu+G308BeSWq0N80IkgZzZEbopQoKyLpRAxdiEqfalzWiEC74r0KMmuEurs4mQA2oRRKJjIXlxdhOM9TJLnIJj4S1oiKRFQuLk4odHAySfr8FMDwPJlkEvn7SpMZFGVWiQgjzRNNXnSXurdU4JYIXVcNopxEdNZIDtytZfWVWxsZtt/5EOArqx3QR2l5EUSN4r0HZMJ6PmIrIHQitEiKQFkRCUUKcXe0Fndp6UhE5rKgbglOKJxMAKl1EgSZCIgBklslQaKKkwD3kqhIRCa+U/BKsZxQdJZJ7txVACI1GYeLC5BbJSnE7EGXaia+3Fuiv+1k3cRrxCO06LrMGpG4taw+cpJwtuGJpiWxVDh0MofHQ4GKLEpWJiVEwShbIhHIdEacA4MXicgIRLZPDAaUUAxdXX7JRECQiZfeQS0QLrjzP2nBxKMjFUoiOitPkK8ih0R77RziHU4Xl8wqyeTKqlACKci9JbNCqIuLtwHcpCHgNV87QVxR+8zPeslRw15zkOkjOqFdFrEVSFJuB8zSs/tvykqfomSZ7f0RaUTsH6rSGgHyg1s7eb8H7kEgzRYK3l5FTA4tJt8HSm5Uw6H9zj9hu7PexZ8rSZ7CAaevmbtr5KGWERe50Ot1yQZ1Tg6y0GBBIl5ZxNty7XbAnRUvziu2S85nVw/O6Qg0612QpHhSF5WBxT4/2eAO6Eigji0C1KXFtylAs9kp+kVoL6DufwEaiUxot/SR0DrpDygrIuFQWiOA2gpREYdqXz8kE/4+LSELjoLdCSp/Oo/ooiVUxMLLpLSzdRWZSHQZUfYjP3eK+1UQjC4nw7MWGh08VeK57hjAXWtL4tZS5RyZ1NhSEQ2vmxYYZNaIYYHTJOIuoT1wfSRE0Sgr15aVjCj5AbZXyV1anEAcx0guQAcA0T4KZzioytWlcXMBOXdKjkwi1Rmpmys/I2BeM5GJ7zSCCzCfGc/SEawbozqDElw/EdYIFeJpO1kkF40GSkveDyftxUC1DcCI/Dmqqp3COw/zlYnuOveWdsIrXrBRwEQv8WGNyGBS8kTl0nJEfYmSNzIEWA5FuLU40oj60kcCKebYARgFf2WLv9T+iLK0SJLJeN4aUbm0ZCTSDqe7i0O23491IiZ2Aqx+eVgnyWTeRUdNfDHgySwTmYvLahtBvuxhHJlcdBt/IiwoSUxGEJREckubbPk0t/9Tchy1PMT2ztx7ao2k8/uqOmAXd9RZJQCCcW8J6MghKtnGrREDUHedc7uZdWH0ICEjENVPQfcTkbi1uD4iQxIxlz6SIr/nEH2LsiMSaV2tTshJhLqpZFFbfKGghELPQ9sG7OoyIRNrn9PFZeJnpseI9pnKXB+EJeXl++biOhHe23KvqU73IkhFmvDIyYRbPcTFJbLwqTjNrZKi3Vu8CrDYJhCFnEBUpMHmHaHg+SN8H4esrTYZMYjKv0W4tQCVRhcPpuIvBZ/BU7eEcKHsXFsAnNaIjkQAKF1dMqiitsRAQKOP9kDu6vLIN/FydYlBXri6uJsLcM7TIdxTOjJxJu3l3TxxJJGsiSDBy6TUwjhBkRIIALTJ2rOqwolcSRQHxPXEZ0ZLcohw4GorHJhGcIn7EsUcuXuL3re29pYAr7tmknskm0qXg+kjOjKTua9kc8z4ziExgaHBqnJrAXkNREYg/IEogwhS/WieknJFeRFJMo69EYlLC9CTiE4r4aBEAbjJgmsnnHC4bgIYhwgLy0QQiimZiHUVVHkVDp3ETzlxhlRnnkC0XE0IxXNKIZrQGIWdIR/vcGologZXJhfRR3NKqJaSRFxfe4tDaCLiO1WVS3FMUqXY7gMy3SqQcvEm8MpgN3BryVyn3K0VRmv1P5QVkbggrBFOIqYEQslC1k4mtHPrxFSIByxC8WmdeJEJYCZWUuKJ5TK940ghUxlBV00GVZ3IkwnNTvdBLuKj4FZJGyzisHk393m4yERihaAdnlaJIETq3uIlU2TJicpCnzrSkFmtAaHXc0MoTHURn9FaTpeqU4AX+5NBKOChy6oolJVGkukc6HZpeZGITN/gmohqOz8W5Jqq/QC02olEjM+0x21SoUI84A4Ppk983H0g9BIqvNN9gnDy4j2J6+cDhEfWOWC5tYQ1wkkkRRaxXbRrQ147cWgmO5AnMvpeopXQvBKZ6M5zSgCLYBwTR4mZE2u7cq4nOK0K8VFzS4PnkNA2BSclyknE7/S6Lng9CMj0HhnI78FvtFYmR/YCXGTvLFUhxyLR0tKCL3zhC6irq8Pw4cNx5plnYsOGDdpjlixZgoqKCsdSXe1DVOojlJVFsrcjDlQwEhFQkQhdF/CqNcQtFZX1UYh1Amj1E2qd2IjTtymp+wCQu0W4+yuvIUTz7q2ajKWTCO2Ci+rUQhD7JWiDc/JAYYkoPfnCMsnV1wLg1Geqc9djAr/MKqGiO810j5F79lV7S2WV+IUBmfjROnotWVGRyU7hFa3FkxC5yJ5/EOqf6eYvv/wy5syZgy984Qvo6urCT37yE0yZMgXvvvsuamrUH0wikXAQTkVF/y9KWVZEIsUeuElERiAmLhrapgZOYpBpIzrtBJI2BvoJd3fJhHgBmhPiBSFE0/bW0zxzbwmIAV1U7/WBNsV7KQSZUFeWgCAw7vKqcWslcQdhxFyiuy/3Fv2uuFYig2piK4Gc0M4hK43i3N+Hri4ZCkxC9LKO09gbaDeDwrJlyxzrS5YswfDhw7Fu3TqccsopyuMqKirQ0NBQ6u4FivIiko4q646pS8uLRCg5mETgiEGBPoFzQtGRRTGEItFPKLhuIkDrTMkQyz2lA26rxBW9xXUSAbFd8RnKCENWlkoZ1fUpkBiOPGmImRX5nCa57SqrhNbfoqK7SOA0Sk4U98itEhmZyFxhCtCILa/cD1kOSa8Ri0Zk57W1gPwkVtytpRPZqVs2kFpbPtHW5vwlRqNRRKP6L3D37t0AgMGDB2vbtbe3Y9SoUejp6cHxxx+PW2+9FUcffXRxHS4xyotIaPa6ikQ4gfA/vsqLIH7LPPxTRigCQREKhYRQXO6uOIzJxBpI3Tko1CpJIoUIMuiqQV50T5NXTixRSC28FPTRW14RXQl6Th6CLFxcBlaJCAX2cm/Rz8EFmVUitpvAhz5SLOI51SHw8igGIrt0SgeNyO62RqIIpra7p+3raHfQQQc5ts6fPx8LFixQHtXT04MrrrgCX/ziF7Xzrx9xxBF45JFHMG7cOOzevRt33XUXTj75ZLzzzjs48MADDfvY+yiZ2H7LLbfg5JNPRjweR319vbQNF5UqKipccxyvWrUKxx9/PKLRKA477DAsWbKkuI4VQiJc+ZVB1obqMHQw1Qn5/DivNu1sXSLICzGeCvEpe/iI2U+DQP5JT/xZecZ7Es4Yf2GVOJITa2ANnmKu8RqykAEy4ZHAmJIs4mMQegoV4B3iezvy9bd40UhS0JHW4IrZMq5bdJfV5AIsK89R90qI7gJ+622pwoF9wLTkjQzaqZMBf3OOMEgLfOYgyx2hIrvKGrG+jd63SLZu3Yrdu3fby3XXXadtP2fOHLz99tuec7g3Nzdj9uzZOO644/CVr3wFTz/9NIYNG4Zf/vKXQXY/cJTMIslkMjjnnHPQ3NyMhx9+WNlu8eLFmDZtmr1OSWfTpk2YMWMGLrvsMjz++ONYsWIFLr74YowcORJTp0713ymZziEjEUogumNlqGHHxeC0IoRuoNNPxHYvC0XWjlspwkKRuLtkugkdJPnTtszFZXUpZ5UI0d0ELDRYK6oT6NokqF4ic3H5sEpUojtgfS7WPecHMLu0PM9PolDV2yqg6owKJrW0KIohHS0UYrtJyXgqsntbI1Z+WPHoBDDQsJ0liCcSntlMAIC5c+fiueeewyuvvOLbqhg4cCAmTJiADz74wNdxvY2SEcmNN94IAJ4WRH19vVJYWrRoEUaPHo27774bAHDkkUfi1Vdfxb333ls4kXQj/8RKSURmhQBuAvESTLk7RZyHEgrd74dQOHy5vXKhwnC6umRkogJ3cdGw2CRiiFcmnTkl3L3FCzACiFUDbTkCiEHtYBARXBTSqT1ULi4qvBtqJV6iu9a9RSc0o1oJJxMZiVBLxEfkZ5+Xj1dZHEWI7ICXNRJHf522MJvN4vvf/z6eeeYZrFq1CqNHj/Z9ju7ubrz11ls4/fTTS9DD4NDneSRz5szB0KFDceKJJ+KRRx5BNptPLlq9ejUmT57saD916lSsXr1ae850Oo22tjbHAsDp0lKRiPChdLDtPFyYQ9aGnoP6ZkQ77u7yyl2hx6ra8Tbilbi7qKtLFH7kbi7q4hL5JdzFlSFuL9spVBMxc2+JbXAShO4ZT+SSiI9T3LZwcyldXLxOVzscBR1FXgkgItGSdoKirJAjAKl7y5FTIiBzT/mut+WO2OLQCemy8ihFQUYYsml0FSK7gInI7mWNWAVYTSyJ3secOXPw2GOP4YknnkBdXR1aW1vR2tqKVCr/25k9e7bDLXbTTTfhxRdfxD/+8Q+88cYbOPfcc7F582ZcfPHFfXELxuhTsf2mm27CV7/6VcTjcbz44ov43ve+h/b2dvzgBz8AALS2tmLEiBGOY0aMGIG2tjakUinEYnLfaEtLi20RucBJhG7nVogqMVAHbg1wwV1noehgaqXApE3eOnHAwzDhhQzpn7xQqyRRYyUXKiOyJNtpfklM0s7l4qqFVdiRg1kl8cp8fgi1SlSiO3dvuUB/C155JaoExRxUc5DECNn5hbZgowD9vri1bQKJNcLzQwSoyE5DfvXWSC4/rGjsgVmNePM/7UMPPQQAmDRpkmP74sWLcf755wMAtmzZggED8s/zn332GS655BK0trbigAMOwAknnIDXXnsNRx11lPF1+wK+iOTaa6/F7bffrm3z3nvvYezYsUbn++lPf2q/nzBhAjo6OnDnnXfaRFIorrvuOsybN89eb2trs6IsZJa/0ERSbF2AZ5arYOdzkG08VNeLUFRRXiq3l9gn01Jk+3l/26tcri4dmVAXGH0i5+XWkQDq0W7NadEBp3uJEkxuW6zDcm8Ja0QQQgxqTUSXrOhwcdUgPz8Jd3HthG0lCa1EWGW0bIoIBda5t8BzSqqR/7141dui8BDYReivVw6J3d4nySRrIt6Cuw4aa0QGLrJza0QkIKqtkUjBNd5KDepdUWHVqlWO9XvvvRf33ntviXpUOvgikquuuspmUhXGjBlTcGeamppw8803I51OIxqNoqGhAdu2bXO02bZtGxKJhNIaATQx3e0AuuB0aalIhLuEvMDbVcM5wAdBKBwm4jzYfpd+UmVsmVAfvHhajORcQdsxFEOx3Sm8D8m1N7RKAKcWIm4tAbm1Imsbo+87WKKiTHjPkQm3SpwVkvOhwDynRFkR2A7BhvPzls3drnJpFVkZoyTVfXXgVgrrv05kl4X8Wq9xh5VivVrH52c4DSJgIHiLpJzgi0iGDRuGYcOGlaovWL9+PQ444ACbBJqbm/HCCy842ixfvhzNzc2FXaAd1ixofknEz2+HljEBnISigqnLy9RKoX0xco9LXF0e3gKhGVgkMgRDsUPt4hL5Vx5WCTrdJcx0VokAdXGJtnWwyMkmEuHikgnvn0JqlVD3lrhXEQ5tVBFYZpUIyMijlr0CRvpIIQicZGQaCfQhv4C73hsP+bW25bU4ao0AgGN20xB9hpJpJFu2bMHOnTuxZcsWdHd3Y/369QCAww47DLW1tfjv//5vbNu2DSeddBKqq6uxfPly3HrrrfjRj35kn+Oyyy7Dgw8+iGuuuQYXXnghVq5ciSeffBLPP/984R3TkYgXgXgRgkwfoT9ylZWiQrFWim9UeZIJd5XkM73jtovLFQ4sBm4hdymsEgCoy31eKTg9QV5WCR8WHVYJkLdANOHAKquEZ7rrSqYISK0S0/lJAkBJsthNqjkL4iDWiCrkl0cBqkJ+na/cGjHoU4iSo2REcsMNN+DRRx+11ydMmAAAeOmllzBp0iQMHDgQCxcuxJVXXolsNovDDjsM99xzDy655BL7mNGjR+P555/HlVdeifvvvx8HHnggfv3rXxcW+gvIn869SET1x6eFEwVU+ohYl1kpxbi9RLtCQ4i5O6wT8CIT+vQYQxI7MASA08VloxLA8F1IdGacwrsosCiiqHKEkoB6gqsU1GQiAxXeE9QCUVklLK9EZpWockqsV4WrlVolgJpMuDZCflti3hlTeE2xa0wyvO+Ac54Xui6BlzUCeIf8UpdXaa2RPTAr/hiylgwlI5IlS5Zoc0imTZvmSERUYdKkSXjzzTeD6ZTQSIQ1IiMRGYGofqyq7dydpa2PpWgjIxQOUytFgOokWpeXmkwEcQB5dw59spS6uIZk8j80MmjbhEIgc3Hlr+dNJvQ40dbIKhHZ7gZWiZfoLqDUSmiOCScQhdgui9iSwSvBsCBLRWaJ6EKAJdYIhZ+QX/GaIrUHHNaI+O5C9Cn6NPy315GEFblFSURARiKcKEzcyjGYubP8EgpQvJXiC2oyUdWXEoOu1sUl9BKN8A7kXVyAO4qLkolX3glgaJWQ/lCrZBfq7fOZiu4uiCd7mZWqA9NHdNPruiLnehuMWLg1Qt1aFCYhv4DGGhHh/EWjE2ZpdaEgI0N5EYkM7dCTiF9NkpdHCep3V6yVUhCheGsmIpIrioxawBUuLmTy7ixNAdTEp3kXF+BMPqRkwsHzSgSMrBJJtjsq83kapqI7kNOLeASXikwoBLFoorVkFocpgRSV+W6ijwAua0QFmcgutstCfq1Xq73DGrEahuhjlBeRtMOqpkBdWpxECiUQGVSk0htWSlDorEKmPQE05J7vFdFcMr1EGsVF3VmaKK4RgDXLIZyuLJNILlnplJSXVUJEd2GVxBPmojuQf2qOx5P5J2fAsixoDS6vcvIeEDkkfpIRS2qtaKwRvyK79Rq3iSaZOyaJWD5vJHBrBLB+ZSYRcv1sjpd+gvIikk5YTy8Fk4ip1Ct5XqaDfVCEAqitlEIivjiIDz/T6iYT8RQeRxI7MBRDsN0xSLjAExVHqJsCar1ERSaqKzsiuARhcatE9IlYJXGfojsPBbatEvEEzcOBDeBXaC8GmcqIvpS8l1ZiaI0AepEd8GGN0CreIfoM5UUkKmhJRBY/pIIYyugxjFSCJhSxrrJSVBCE4sPlJSMTAHYyIvV/a/WSIWywkuklyH1ynwKf83APqlxarnY0251n3IsBSewbYpWY9yu6A/lQYKlVIouEoigyCVEGVeXfguBRmFFmjQgUIrKLNgCpHCA+P1pbLkSforyIhBbtE9aIkkQoGZj+EWXPyApZOAhC4e3EehCuL1lEEdxkIquAqw0/9aOXpHOzHkrIpA35PBMVifBHgASsgo52DS5hlQDuysCdZqHAstkTbdAIrva4OZkADqHdNGJLhqILNpr0lYFbI0GJ7Jn2nN4k/reBWiPC7+2FXq4WsI+g/IikE04SEZCSCP/R+HVt0RE8YELh7XhbsV6C0EhOJkOwwxEWLBBDylwv0UBEcunIhEJVcaQNOdIRVgjNiTBIUKQQFojKFcMRqU0654ORDdAlsEYKRVeNNfGXERTWiAwqkR2Qu7VoyC8AtzVCp8cO0WcoLyKRoRMeJMLJwyRZg7u2AiYUsd6HpCIjk3xtKoMnaKGXiHWaX0LFd9H8U+u1zufTMYen6K5JUFRNemVildgDIRXeS0wcvRIOrLkHnchu7ddnsktDfgF3hGUgVkmhXocQQLkRSQeAvXC6tIxIxE9oCK/GR4kjAEIB/JMKR0CkQp+wY5IchxjRDWSQ5pfI0AFgeJ5M6L36zXgHcu4tleiuSFCUhQLTSa/Ee6tP7nt2CO88iosj59aiQrsuh6TkMHBvCWvES2QH3G4tAZrJbiSy09cQfYo+n9iq/8CERNoMF5Bj97BjxTVSku0EtIlwxwm0w+ma4+udbOFt6CKevsX16DY6b0s7e98OYHuVPUHWjuRQO4FPROTswFB8hnqkcoKqcGuIibB2VdajbXgEGAIry30wrEiuwXBvi1rriRpgRLVF0QlYfBvLvReUnGDr4lMWE18ByBMpLdPSSbaTbdbEVyntpFdiO2ARjtU296xNSMAmB1UxxhIVaSwJFNaIX7cWzR0R69arh8i+xz6gX2PhwoU45JBDUF1djaamJvzlL3/Rtn/qqacwduxYVFdX49hjj3UVru2PKC+LRFggUmsEcA7unEBk7ThUUVs0K5BbKPQYw7BhwHv+E7rNJNO+YOSTFlNx6zWKtCOSy1h8B+Ql5wW2wa7JBcC+L6F9qBIVOVzuLSBPngWEAlvddlol1ra8/1+QSaYzktdLBGnYZVT6MYlwq4SQiMwa8ePWch6jyWQH3CJ7YHkke2C5K7zgz7f629/+FvPmzcOiRYvQ1NSE++67D1OnTsWGDRswfLhbJHzttdcwc+ZMtLS04Gtf+xqeeOIJnHnmmXjjjTdwzDHH+Lp2b6IsiMSeYKa7zarLlgWbeqAN6phC+kv1+hHRRyPxTxOPuKLQlZCCxfs9pL2m8AflrmrJOu+qGIDb2Dbq0uom23eRV/F/yuSWNPKElMq9T+aWDgBdQKamC5+1DUDP4G70oAr16EEEwF5EkEEU9ehEBSrRg0pYtfwj6EEVuhFFNwaiqzuDqh5Yg3cd8qyQzvVhACwrpRvWd9gDxLJAezpvVouSe3y8E7c8IHflCgB724FELHcPu2BZPuJed5PGuRN3VqfRXdmBvYgii4HowkBUoAoDUIEKpFGBAahCBfYihi4AXciiCj2IIoIexJBCDNlkNwYCyHQOxEC0uWf2I9/fwJoksm1ApHovshmgp8salnsAdCOFbmTQhRS6kMJeJFGBDDJIoRKdqEAGNdiFTmTQgwwqkAGwF13Ym+tZF3rQjb3oRjd6cov1hxjYnUVX7gOrkj3pd+Y+VDLmduV+R8lu6xx7K7NIogcp9OR62p3raRd2oQIpVKETGexCHGl0IoMB2I0o9iKJLkTQjmp0ox0diCCb3IPMjkFAR5uVoJpE3nrenfu+PrN+5CYTSalh6iOz2tnTd+egmgNJFKK94IILAACLFi3C888/j0ceeQTXXnutq/3999+PadOm4eqrrwYA3HzzzVi+fDkefPBBLFq0yMf99DKyZYCtW7cK6giXcAmX/XTZunWr77EhlUplGxoafF2ntrbWtW3+/Pmuc6fT6WxlZWX2mWeecWyfPXt29utf/7q0PwcddFD23nvvdWy74YYbsuPGjfN9b72JsrBIGhsbsXXrVtTV1aGiosL7gH0IYhrhrVu3IpEwcezs/wg/Ezf2588km81iz549aGxs9H1sdXU1Nm3ahEzGPFcnm826xhGZNbJ9+3Z0d3djxAhnCYcRI0bg73//u/Tcra2t0vatra3G/esLlAWRDBgwAAceeGBfd6OkSCQS+90AUSzCz8SN/fUzGTRoUMHHVldXo7q6HyXx7IMIo7ZChAgRogQYOnQoKisrsW3bNsf2bdu2oaGhQXpMQ0ODr/b9BSGRhAgRIkQJEIlEcMIJJ2DFihX2tp6eHqxYsQLNzc3SY5qbmx3tAWD58uXK9v0FZeHa2p8RjUYxf/58qY+2XBF+Jm6En0nfYN68eTjvvPMwceJEnHjiibjvvvvQ0dFhR3HNnj0bn/vc59DS0gIA+OEPf4ivfOUruPvuuzFjxgwsXboUa9euxa9+9au+vA1PVGSzRcXMhQgRIkQIDR588EHceeedaG1txXHHHYcHHngATU1NAKypxA855BDHtORPPfUUrr/+enz44Yc4/PDDcccdd+D000/vo96bISSSECFChAhRFEKNJESIECFCFIWQSEKECBEiRFEIiSREiBAhQhSFkEhChAgRIkRRCIlkH8Ett9yCk08+GfF4HPX19dI2W7ZswYwZMxCPxzF8+HBcffXV6OpyVpRdtWoVjj/+eESjURx22GGOaJH9AX5Ldu/LeOWVV3DGGWegsbERFRUVePbZZx37s9ksbrjhBowcORKxWAyTJ0/G+++/72izc+dOzJo1C4lEAvX19bjooovQ3h5Ogh7CH0Ii2UeQyWRwzjnn4PLLL5fu7+7uxowZM5DJZPDaa6/h0UcfxZIlS3DDDTfYbTZt2oQZM2bg1FNPxfr163HFFVfg4osvxh/+8Ifeuo2SQpTsnj9/Pt544w2MHz8eU6dOxaeffup98D6Ijo4OjB8/HgsXLpTuv+OOO/DAAw9g0aJFWLNmDWpqajB16lR0dubrI8+aNQvvvPMOli9fjueeew6vvPIKLr300t66hRD7C/qwYGSIArB48eLsoEGDXNtfeOGF7IABA7Ktra32toceeiibSCSy6XQ6m81ms9dcc0326KOPdhz3ne98Jzt16tSS9rm3cOKJJ2bnzJljr3d3d2cbGxuzLS0tfdir3gEAR5XZnp6ebENDQ/bOO++0t+3atSsbjUaz//f//t9sNpvNvvvuu1kA2b/+9a92m//3//5ftqKiIvvPf/6z1/oeYt9HaJHsJ1i9ejWOPfZYR+XQqVOnoq2tDe+8847dZvLkyY7jpk6ditWrV/dqX0uBTCaDdevWOe5vwIABmDx58n5xf36xadMmtLa2Oj6PQYMGoampyf48Vq9ejfr6ekycONFuM3nyZAwYMABr1qzp9T6H2HcREsl+AlX5abFP16atrQ2plGrWx30DupLd/b0Edykg7ln3ebS2trpm6auqqsLgwYPL8jMLUThCIulDXHvttaioqNAuqnkLQoQIEaK/ICza2Ie46qqrcP7552vbjBkzxuhcDQ0NrgglUY5alKBWlahOJBKIxZxzau9rKKRk9/4Mcc/btm3DyJEj7e3btm3DcccdZ7fhgQhdXV3YuXNnWX5mIQpHaJH0IYYNG4axY8dql0gkYnSu5uZmvPXWW46BYfny5UgkEjjqqKPsNvtiiWoTFFKye3/G6NGj0dDQ4Pg82trasGbNGvvzaG5uxq5du7Bu3Tq7zcqVK9HT02MXFQwRwgh9rfaHMMPmzZuzb775ZvbGG2/M1tbWZt98883sm2++md2zZ082m81mu7q6ssccc0x2ypQp2fXr12eXLVuWHTZsWPa6666zz/GPf/wjG4/Hs1dffXX2vffeyy5cuDBbWVmZXbZsWV/dVqBYunRpNhqNZpcsWZJ99913s5deemm2vr7eEcm2P2HPnj327wBA9p577sm++eab2c2bN2ez2Wz2tttuy9bX12d/97vfZf/2t79lv/GNb2RHjx6dTaVS9jmmTZuWnTBhQnbNmjXZV199NXv44YdnZ86c2Ve3FGIfRUgk+wjOO++8LADX8tJLL9ltPvzww+z06dOzsVgsO3To0OxVV12V3bt3r+M8L730Uva4447LRiKR7JgxY7KLFy/u3RspMX7+859nDz744GwkEsmeeOKJ2ddff72vu1QyvPTSS9LfxHnnnZfNZq0Q4J/+9KfZESNGZKPRaPa0007LbtiwwXGOHTt2ZGfOnJmtra3NJhKJ7AUXXGA/nIQIYYqwjHyIECFChCgKoUYSIkSIECGKQkgkIUKECBGiKIREEiJEiBAhikJIJCFChAgRoiiERBIiRIgQIYpCSCQhQoQIEaIohEQSIkSIECGKQkgkIUKECBGiKIREEiJEiBAhikJIJCFChAgRoiiERBIiRIgQIYrC/wfKBCDEvKmstAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 400x300 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "grids_ = np.linspace(-180, 180, 100, endpoint=True, dtype=float)\n",
    "plt.figure(figsize=(4,3))\n",
    "plt.contourf(grids_, grids_, metad_fe_sumhills, cmap=plt.get_cmap(\"jet\"), levels=100)\n",
    "cbar = plt.colorbar()\n",
    "cbar.set_label(\"FE (kcal/mol)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a540362d",
   "metadata": {},
   "source": [
    "## $\\textit{S$_{config}$ for Ground Truth}$"
   ]
  },
  {
   "cell_type": "raw",
   "id": "1f24e1b2",
   "metadata": {},
   "source": [
    "S_ref = numerical_integration_2d_simp(prob_metad*np.log(prob_metad), dx, dy)\n",
    "print(\"Configurational Entropy for GT= \", S_ref)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "c05a8578",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Conigurational Entropy for GT=  -1.0879662775166619\n"
     ]
    }
   ],
   "source": [
    "locs = np.nonzero(prob_metad.T.flatten())\n",
    "S_ref = dx*dy*np.sum(prob_metad.T.flatten()[locs]*(np.log(prob_metad.T.flatten()[locs])))\n",
    "print(\"Conigurational Entropy for GT= \", S_ref)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcf1dbd9",
   "metadata": {},
   "source": [
    "## $\\textit{calculate S$_{config}$ for all GMM objects by integrating all space}$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "d491d189",
   "metadata": {},
   "outputs": [],
   "source": [
    "from shapeGMMTorch import similarities as sm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "636bbc2a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Configurational Entropies (all): \n",
      " [[43.12544939 55.02468127 55.96272221 55.33515733 55.90894127]\n",
      " [58.71385013 59.7263472  60.41080449 61.01982639 61.44010854]\n",
      " [58.61144214 60.01207617 60.64230191 61.06923379 60.94503501]\n",
      " [58.67150763 60.04222813 60.34279791 60.93067555 61.06024796]\n",
      " [43.12532024 55.02682317 55.95976985 55.33266165 55.90882048]]\n"
     ]
    }
   ],
   "source": [
    "num_samples = 1000000  # 1M samples\n",
    "S_config_all = np.zeros([all_wsgmm_combined.shape[0], all_wsgmm_combined.shape[1]])\n",
    "\n",
    "for i in range(all_wsgmm_combined.shape[0]):\n",
    "    for j in range(all_wsgmm_combined.shape[1]):\n",
    "        all_wsgmm_combined[i,j].device = \"cpu\"\n",
    "        S_config_all[i,j] = sm.configurational_entropy(all_wsgmm_combined[i,j], n_points=num_samples)[0]\n",
    "\n",
    "print(\"Configurational Entropies (all): \\n\", S_config_all)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b83e4a0",
   "metadata": {},
   "source": [
    "## $\\textit{calculate S$_{config}$ for all GMM objects by integrating 2D phi-psi space}$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "9406eede",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Configuartional Entropies (2D): \n",
      " [[-3.51434295 -3.27133606 -2.56835631 -3.32795525 -2.61198669]\n",
      " [-0.74024087 -0.28176202 -0.2614884  -0.28619679 -0.33923085]\n",
      " [-0.77568305 -0.29112475 -0.26107979 -0.2998121  -0.25785576]\n",
      " [-0.74728867 -0.28067347 -0.266921   -0.3944218  -0.2834189 ]\n",
      " [-3.5157787  -3.27315366 -2.56988602 -3.32937728 -2.61190938]]\n"
     ]
    }
   ],
   "source": [
    "num_samples = 1000000  # 1M samples \n",
    "nbins = [100,100]\n",
    "range_ = [-180,180]\n",
    "\n",
    "S_config_2d = np.zeros([all_wsgmm_combined.shape[0], all_wsgmm_combined.shape[1]])\n",
    "\n",
    "for i in range(all_wsgmm_combined.shape[0]):\n",
    "    for j in range(all_wsgmm_combined.shape[1]):\n",
    "        all_wsgmm_combined[i,j].device = \"cpu\"\n",
    "        \n",
    "        # generate trajetory\n",
    "        traj = all_wsgmm_combined[i,j].generate(num_samples)\n",
    "        \n",
    "        # calculate dihedrals\n",
    "        phis, psis = calc_phi_psi(traj)\n",
    "            \n",
    "        # histogram\n",
    "        hist = np.histogram2d(phis, psis, bins=nbins, range=[range_, range_], density=True)[0]\n",
    "        hist = hist/numerical_integration_2d_simp(hist,dx,dy)  # normalized\n",
    "        \n",
    "        # get the S_config_2d\n",
    "        #S_config_2d[i,j] = numerical_integration_2d_simp(hist*np.log(hist), dx, dy)   # giving nan values mostly!\n",
    "        \n",
    "        locs = np.nonzero(hist.T.flatten())\n",
    "        S_config_2d[i,j]  = dx*dy*np.sum(hist.T.flatten()[locs]*(np.log(hist.T.flatten()[locs])))\n",
    "\n",
    "print(\"Configuartional Entropies (2D): \\n\", S_config_2d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "66ae8db0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dS_config_all: \n",
      " [[44.21341567 56.11264755 57.05068849 56.42312361 56.99690755]\n",
      " [59.80181641 60.81431348 61.49877077 62.10779267 62.52807482]\n",
      " [59.69940841 61.10004245 61.73026819 62.15720006 62.03300129]\n",
      " [59.75947391 61.13019441 61.43076419 62.01864182 62.14821424]\n",
      " [44.21328652 56.11478945 57.04773613 56.42062793 56.99678676]] \n",
      "\n",
      "dS_config_2d: \n",
      " [[-2.42637667 -2.18336978 -1.48039003 -2.23998897 -1.52402041]\n",
      " [ 0.34772541  0.80620426  0.82647788  0.80176949  0.74873542]\n",
      " [ 0.31228323  0.79684153  0.82688649  0.78815418  0.83011051]\n",
      " [ 0.34067761  0.80729281  0.82104528  0.69354448  0.80454738]\n",
      " [-2.42781243 -2.18518738 -1.48191974 -2.241411   -1.5239431 ]]\n"
     ]
    }
   ],
   "source": [
    "dS_config_all = S_config_all - S_ref\n",
    "dS_config_2d = S_config_2d - S_ref\n",
    "\n",
    "print(\"dS_config_all: \\n\", dS_config_all, '\\n')\n",
    "print(\"dS_config_2d: \\n\", dS_config_2d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "1fe00a76",
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_tabular(dS):\n",
    "    dS = dS.T\n",
    "    y = dS.shape[0]\n",
    "    x = dS.shape[1]\n",
    "    \n",
    "    x_labels = [\"uniform\", \"bias\", \"rbias\", \"fbias\", \"uniform_modf\"]\n",
    "    #y_labels = [\"Sizes\", \"K=2\", \"K=3\", \"K=4\"]\n",
    "    y_labels = [\"Sizes\", \"K=2\", \"K=3\", \"K=4\", \"K=5\", \"K=6\"]\n",
    "    \n",
    "    for i in range(y+1):\n",
    "        if i == 0:\n",
    "            print(\"size|     \", \"     \".join(x_labels))\n",
    "            print(\"--\"*33)\n",
    "        else:\n",
    "            print(y_labels[i], \"|\", dS[i-1,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "d875d241",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dS_config_2d :\n",
      "\n",
      "size|      uniform     bias     rbias     fbias     uniform_modf\n",
      "------------------------------------------------------------------\n",
      "K=2 | [-2.42637667  0.34772541  0.31228323  0.34067761 -2.42781243]\n",
      "K=3 | [-2.18336978  0.80620426  0.79684153  0.80729281 -2.18518738]\n",
      "K=4 | [-1.48039003  0.82647788  0.82688649  0.82104528 -1.48191974]\n",
      "K=5 | [-2.23998897  0.80176949  0.78815418  0.69354448 -2.241411  ]\n",
      "K=6 | [-1.52402041  0.74873542  0.83011051  0.80454738 -1.5239431 ]\n"
     ]
    }
   ],
   "source": [
    "print(\"dS_config_2d :\\n\")\n",
    "print_tabular(dS_config_2d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "40b5e7e5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dS_config_2d * -1: \n",
      "\n",
      "size|      uniform     bias     rbias     fbias     uniform_modf\n",
      "------------------------------------------------------------------\n",
      "K=2 | [ 2.42637667 -0.34772541 -0.31228323 -0.34067761  2.42781243]\n",
      "K=3 | [ 2.18336978 -0.80620426 -0.79684153 -0.80729281  2.18518738]\n",
      "K=4 | [ 1.48039003 -0.82647788 -0.82688649 -0.82104528  1.48191974]\n",
      "K=5 | [ 2.23998897 -0.80176949 -0.78815418 -0.69354448  2.241411  ]\n",
      "K=6 | [ 1.52402041 -0.74873542 -0.83011051 -0.80454738  1.5239431 ]\n"
     ]
    }
   ],
   "source": [
    "print(\"dS_config_2d * -1: \\n\")\n",
    "print_tabular(dS_config_2d*-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "147afeaa",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dS_config_all :\n",
      "\n",
      "size|      uniform     bias     rbias     fbias     uniform_modf\n",
      "------------------------------------------------------------------\n",
      "K=2 | [44.21341567 59.80181641 59.69940841 59.75947391 44.21328652]\n",
      "K=3 | [56.11264755 60.81431348 61.10004245 61.13019441 56.11478945]\n",
      "K=4 | [57.05068849 61.49877077 61.73026819 61.43076419 57.04773613]\n",
      "K=5 | [56.42312361 62.10779267 62.15720006 62.01864182 56.42062793]\n",
      "K=6 | [56.99690755 62.52807482 62.03300129 62.14821424 56.99678676]\n"
     ]
    }
   ],
   "source": [
    "print(\"dS_config_all :\\n\")\n",
    "print_tabular(dS_config_all)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "9ab59458",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "S_config_2d :\n",
      "\n",
      "size|      uniform     bias     rbias     fbias     uniform_modf\n",
      "------------------------------------------------------------------\n",
      "K=2 | [-3.51434295 -0.74024087 -0.77568305 -0.74728867 -3.5157787 ]\n",
      "K=3 | [-3.27133606 -0.28176202 -0.29112475 -0.28067347 -3.27315366]\n",
      "K=4 | [-2.56835631 -0.2614884  -0.26107979 -0.266921   -2.56988602]\n",
      "K=5 | [-3.32795525 -0.28619679 -0.2998121  -0.3944218  -3.32937728]\n",
      "K=6 | [-2.61198669 -0.33923085 -0.25785576 -0.2834189  -2.61190938]\n"
     ]
    }
   ],
   "source": [
    "print(\"S_config_2d :\\n\")\n",
    "print_tabular(S_config_2d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9b63acc",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65f10259",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8c51911",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "644f20be",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5e9fa1f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ca1b2b1",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}