experiment3_ala2_MC_RNVP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alanine dipeptide Boltzmann Generator\n",
    "------\n",
    "A simple example of how to train a Boltzmann generator on a small system by computing forces and energies in OpenMM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import rcParams\n",
    "import mdtraj as md\n",
    "import matplotlib as mpl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import SVG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bgtorch.nn.flow.base import Flow\n",
    "from bgtorch.nn import DenseNet\n",
    "from bgtorch.nn.flow import (\n",
    "    SequentialFlow, \n",
    "    CouplingFlow, \n",
    "    AffineFlow, \n",
    "    SplitFlow,\n",
    "    MergeFlow,\n",
    "    InverseFlow, \n",
    "    SwapFlow\n",
    ")\n",
    "from bgtorch.nn.flow.transformer import AffineTransformer\n",
    "from bgtorch.distribution.energy import Energy\n",
    "from bgtorch.distribution.sampling import Sampler\n",
    "from bgtorch.distribution import NormalDistribution\n",
    "from bgtorch import BoltzmannGenerator\n",
    "from bgtorch.nn.training import KLTrainer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bgtorch.nn.flow.stochastic import MetropolisMCFlow, LangevinFlow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "rcParams.update({'font.size': 16})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Switch AUTORELOAD ON. Disable this when in production mode!\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ala2\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_2ala_ramachandran(traj, ax=None, weights=None):\n",
    "    if ax==None:\n",
    "        ax = plt.gca()\n",
    "    \n",
    "\n",
    "    if isinstance(weights, np.ndarray):\n",
    "        ax.hist2d(md.compute_phi(traj)[1].reshape(-1), \n",
    "                  md.compute_psi(traj)[1].reshape(-1), \n",
    "                  bins=[np.linspace(-np.pi, np.pi, 64), np.linspace(-np.pi, np.pi, 64)],\n",
    "                  norm=mpl.colors.LogNorm(),\n",
    "                  weights=weights)\n",
    "    else:\n",
    "        ax.hist2d(md.compute_phi(traj)[1].reshape(-1), \n",
    "                  md.compute_psi(traj)[1].reshape(-1), \n",
    "                  bins=[np.linspace(-np.pi, np.pi, 64), np.linspace(-np.pi, np.pi, 64)],\n",
    "                  norm=mpl.colors.LogNorm())\n",
    "        \n",
    "    ax.set_xlim(-np.pi,np.pi)\n",
    "    ax.set_ylim(-np.pi,np.pi)\n",
    "    ax.set_xlabel(r'$\\phi$')\n",
    "    ax.set_ylabel(r'$\\psi$')    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from openmmtools.testsystems import AlanineDipeptideVacuum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "ala2 = AlanineDipeptideVacuum(constraints=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "temperature = 1000\n",
    "n_atoms = ala2.mdtraj_topology.n_atoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from simtk import openmm, unit\n",
    "INTEGRATOR_ARGS = (temperature*unit.kelvin, \n",
    "                   1.0/unit.picoseconds, \n",
    "                   1.0*unit.femtoseconds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bgtorch.distribution.energy.openmm import OpenMMEnergy, OpenMMEnergyBridge\n",
    "ala2_omm_energy_bridge = OpenMMEnergyBridge(ala2.system, unit.nanometers, \n",
    "                                            openmm_integrator=openmm.LangevinIntegrator, \n",
    "                                            openmm_integrator_args=INTEGRATOR_ARGS)\n",
    "\n",
    "ala2_omm_energy = OpenMMEnergy(3 * n_atoms, ala2_omm_energy_bridge)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bgtorch.utils.train import linlogcut\n",
    "\n",
    "class RegularizedEnergy(Energy):\n",
    "    def __init__(self, energy_model, high_energy, max_energy):\n",
    "        super().__init__(energy_model.dim)\n",
    "        self.energy_model = energy_model\n",
    "        self.high_energy = high_energy\n",
    "        self.max_energy = max_energy            \n",
    "\n",
    "    def _energy(self, x, temperature=None):\n",
    "        U = self.energy_model.energy(x)\n",
    "        U_reg = linlogcut(U, high_val=self.high_energy, max_val=self.max_energy)\n",
    "        return U_reg\n",
    "    \n",
    "class AugmentedEnergy(Energy):\n",
    "    def __init__(self, base_model, dim_augment=None):\n",
    "        if dim_augment is None:\n",
    "            dim_augment = base_model.dim\n",
    "        super().__init__(base_model.dim + dim_augment)\n",
    "        self.base_model = base_model\n",
    "    \n",
    "    def _energy(self, x, temperature=None):\n",
    "        U = self.base_model.energy(x[:, :self.base_model.dim])\n",
    "        K = 0.5 * (x[:, self.base_model.dim:] ** 2).sum(axis=1, keepdim=True)            \n",
    "        return U + K"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load training data\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "ala2_pdb = md.load('./data/alanine-dipeptide.pdb').topology\n",
    "training_data_traj = md.load('./data/ala2_1000K_train.xtc', top=ala2_pdb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi = md.compute_phi(training_data_traj)[1].reshape(-1)\n",
    "psi = md.compute_psi(training_data_traj)[1].reshape(-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVsAAAFOCAYAAAAo1L7FAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJztvXuUZXd13/nd912vrnc/qiV1t5AECBss1EYSLJATRGxmRmZ5BTssJraF4wiSkCEB/IjH42Hs2CyPDY7Dij1WvMD2LAMzhjhAYuQBGWSWQSINBAk9kFpSq1v0q7qq63Wr7vs3f1Sh7r33T31OX3WdW7fq+1lLq/W79Tvn/M7j/urUd+/fd0sIAYQQQjaXXK8HQAghOwFOtoQQkgGcbAkhJAM42RJCSAZwsiWEkAzgZEsIIRnAyZYQQjKAky0hhGQAJ1tCCMmAQq8HkBWlsYEwsHfX8+2RQs33kZZql00bAHLQK+46ENcnb/rYNXqNkHfbtMxn7cjvwabpI27PQL1T1PsJenw58ds02n48lmDOs93247OLEaWu+0jH77ewpjfK1dt+v7W6at9w87WXGmpf8OQj31Pt61+x/4r0sddKKuXEsdhtYtvFjn25Y4ntN9YHwwN6m5Z/cDpl830p6ecz9qzZr534Rw1tc7lC0X9fJK8/qz998lwIYdrvTdM3k62I/CiAXwJwI4BxALMAvgrgAyGER5O2H9i7C6+75x8937596gnX50DpnGofLJ5zfUakqdqrwV/CkZzu0wx60jnRGnPbnG2NqPZyZ8D1+V59XLWLOf+0PFXV93ylqZ+eSkGPDQCeW9bjSTMhn18edH06ZgLOP6n75CPfq6mH9S+0oaPnXZ/2Y0+q9heO/IXfUZ/xY6/4FdW+98hvXZE+9lrlr70+cSx2m9h2sWNf7lhi+4316Ry+SbWLZ1dcn+p1+ruwMqOfz9iz1hjRE3Jp2T/nS4d0n/r+hutTGtGfHf2pX3vWH83TN5MtgAkA3wDwB1ifaK8B8MsAHhCRHwwhpDphQgjpBX0z2YYQPgHgExd/JiJfB/A4gLcC+FAvxkUIIWnom8n2BZjb+Nf/bWyo5Jp46ciZS/YpGRGnBC/8lIwYNJLzf2bUzV8nebPfmcKi22Ykt6baJ5qTrs9KoaLa7eB106my/pNr/8CCap+ta7kCAF49dUK1H1vc6/oMl/TfZZ3gteq5+WE9voNaFy+c8vrh3I36Eax8zv9Z+YVO/8sGFvvn85tyP5m4TayPvTb2T/k0x75S19fut3bna1yf0pL5qu6+yfXJ3f8t1V54262uz9i351R7ZWZKtUef8TGZ+Zfp78+5myKOh3n9/S4O++/3VZP6O3XU7yVK32UjiEheREoicj2APwJwGsAnezwsQgi5JH032QJ4EEAdwBMAXgng74cQzsY6isjdInJERI6snY8o5oQQkhH9ONn+NIBbAbwdwBKAL4jIwVjHEMI9IYTDIYTDA+PJKTCEELJZ9N1kG0J4LITw4EbA7I0AhrGelUAIIVuWvg6QhRAWROQogOuS+nYgWOuUnm8XIxnNRbOIYTTn4242ZFYSHyiazpdUe7GjRfZ2J5JNbQJt04Ul12Ukr4NoC+2hyPj0eE7XR1V7IO/PyX728tHTrs98Qx+r3vaPTmlan9fZhWHTw/910TLpup3bfcDk1rfrRBMbHAGAex/ReaBpAk69DLx1c+xY8Mue55U6J5fj+0hynq09duwe2Pt75mafT74Hus/IJx9wfaqR4NvFnLrV77e5yyw2qvjvodTMYommX/BTbZTcZ2nouzfbixGRPQBeBuCpXo+FEEIuRd+82YrIXwL4JoCHsK7V3gDgXwNogTm2hJAtTt9MtgAeAPBTAN4HoATgBIAvA/hgCOFY74ZFCCHJ9M1kG0L4bQC/faX2N11Ydp9VjO9BLbJo4EBBX7J6iJnV6O3KYjwDENGLjRo8mVt1feY6WuDcW1hwfaodrYtOmfM81/KLGlbbyRrUWFHrxa2BZAVqqKh16CeWvWY7dFQb5zx3u9faRo5Hks8TSKNdbpbeeaWw44vp2bnHdDuN1prVecf2a4991dlk74Zzd93mPpt6UPuWDJmVBdY7AfALaKRVdH1CQT9rzTEfkxkt+wUTaehrzZYQQvoFTraEEJIBnGwJISQD+kazfbG0OnmcrV3I+3ym7L1+h3JmSW/e67orQeu6kzmvMeaNRjsIrYmWReufAFAzxi5Wn13vY43BL/935WjeH9vu5yUVv/rZar3TJX9tjq9NqPZqQZ93ftDr26tX62OX5nxeozUVae62+bueNLrkVtNok4j5ut67SeeQxign6fp1a5xj9daYqYztU/nc11V77Rav81pj8PaAjwW0JvT3u1Txz+y5Vf/dTAPfbAkhJAM42RJCSAZwsiWEkAzgZEsIIRmwYwJkQFCFDGNGNDmzsMBWZQCAmtHUY4saitBBnsWOFvhj5jUVM56DkWoO32vrwNB82weKbAXeND8/UPaFLS3jhapqn295E5x9FT1ma4ozMuyDc9W8Pu8GfMDRmorsetbflzf+PZ28b98iugnw9JrNGl/+5ckLCTZrUUgag6CSqd5gKzcAQMV9ohmYi5g9me9lbcJ/D1cGdZ9Gyx9pbNIvJkoD32wJISQDONkSQkgGcLIlhJAM2DGabU4CKhHj7EtxrDXqPnu1SeZvRkxlYHTdimgdaLbtdd6ykY++0/DVddtGA11o++RqW6W3FvTCgsmiT4yvREzSLbWOXlCxmvPmNcvtSytpV495reuRpX2qHTN0Li/od4JWxWtti4f0sUeNAfV9X0o23o6x1XXdJLqt2puGbq5NGoNxS0xjXniV/n5YQ/m1yUvHLl6IXE0/W+0J/zzOV7mogRBCtiycbAkhJAM42RJCSAbsGM22E3JYbV3QGWP5ps82tDnNteUzrs+jTa0NzuS9ybfN4W0aDXck5zXHp1taA7Um4AAwlq+6zyzNoG9pRbSBdy6SO9ww18LuAwDyJgfZGq0DwFRR69kLTa1t1Qt+v9OTepszJ7zp8/yr9PWc/rq/d+N/8jX3WRL9rscCybmu3Z5jmhzabopCWtLsN0axqp/Hc7dMqbY1FweAk3fo7/fwSf9dsO+fy6P+u9rqdPeOyjdbQgjJAE62hBCSAZxsCSEkAzjZEkJIBuyYAFk7CJYbF4Jbjy7vc30Ojz6r2t9rTrg+MAU5x2x1BwC5oCNiQzktxC92kqvFTheW3GerJmgW62OrLoyZAN5sa5fbpih6kYWrWAGgaBZvxIx8qh0d5LvWVOBda8+4bQrm2ozt9RUgFk7rKhEF6wYEX3k2Zl6yHcmqMm43Cx+6HZsNtMWObc1qSks6YHvmDTpgBgCDs/pZiy2O6diCu5EY2sp5LmoghJAtCydbQgjJAE62hBCSATtGs+0EwUrzgqY4UvIVO63ZylREE7UJ/8daPgl/f16baDdNErQ1Cl/fr07UjxmDW5OZXKy6rpehFNbMBnAyNIbMQgjAV/aN6bqWUvG8ah8reHOdiQGtKS+uJdlCA2cP+88mvqO3G79f/7wfzcMtac6hG0Pv2CICW103xmZdPzuemBFNPlJp+GK+FTEeev1bfke1h075hTntsn6O6ot+ihy4gebhhBCyZeFkSwghGcDJlhBCMmDHaLbNVh5n5i6YgU8OeAOZ0w2dgxoz1c6LzvGsiP99VcsZMxjoPNbZji9qaJnOe704Zk6T1GehrQszTua91jVo9NdSxBDdGs804M1grK47Z3TnvSV/TmttnZt7Ijfm+hRGzH2I6Gg299bqfDENshsjlSy13zT6a1L+a5rxdmMgs5nY8dhinoDPo7bnFNtmKEHnBYC6KS7aLvtE25VnfVGBNPTNm62IvFVEPi0iz4rImoh8V0Q+KCIjyVsTQkhv6ZvJFsD7AbQB/AqAHwPwhwD+GYAviEReLwkhZAvRTzLCnSGE2Yva94vIPIA/BfAjAP6mJ6MihJAU9M0boZlov89/2/h3f5ZjIYSQy6Wf3mxj3L7x72NJHUWAfOFC4Gep7pPnzxa0/JsTb3gynNeLIZaD34+tqFA1FW6tWQzgg0uxagkHi9p9/snGnsT9jOT0eKMVIHI6WFgRX/3XXot2xy+OiJ3XxQzm/UKIqqlQUS76YzfKOkBWz/lrvnitMeD5tv75mX/5WrfNno989QXHuhVIE3hL6hMLkG2WqcxmBRxjlZFjAbCLiRkRNY1ZUXWf/y7s/239TJx7522uT20yYeXQC42pq622ACKyH8CvA/hiCOHIC/S5W0SOiMiR9lJySRlCCNks+nKyFZFhAJ8B0ALwjhfqF0K4J4RwOIRwOL9r6IW6EULIptN3MoKIVAB8FsC1AG4PITzX4yERQkgifTXZikgRwKcBvAbAHSGEh9NuGzqCZu3C6dbbPin/xLJOqJ8qeelhsa2TnjsRnfIpo6XahQSxyr52YUE+UgX3ey2f8G+ZMMeyGm4eXode6GgzZGt4AwBD0LppTFO221kz8ZWIvj1S1Jpyven3u7qor7lY5xwARkrHsz+hK6kOnUw2bE+jZdbufM1l76eXhjdpjh0775j5y2bQzUINALgvwYAnNv4vGO03ZsBTNfe3OuP12VBIfpZi9M1ku5FL++cA3gjgfwwhPNDjIRFCSGr6ZrIF8B8A/CSA3wRQFZFbL/rZc5QTCCFbmX4KkL1549//FcDXzH8/36tBEUJIGvrmzTaEcPBF7aAlyM1e0BCroyXXpZjXBizHV70xeM5oqcOR3NFD5bOqbbVLa2YDeG3VGr8AwImmNt+Oab/TBV0w0eqoMS3YMpLz5uFNo01bjRkA5sz4SsYk/VTDG3gsNrUeOzHkDYKaRl9fjeTZNkZ1zuSeIzpf1xYEBHyRyGLEqMQa2FQ+93XXp5cabdKx0+ifMX0zTY5sN3m1abgS1y+N+Xn1Ov/9HjqqDe8LN067PoNnuhtTP73ZEkJI38LJlhBCMoCTLSGEZEDfaLYvFmkDxaULv1tWj3r9cOLlc6o9vzbo+gwXtFa5VvDa7xP5fapdNibk9Y5PFM2Z/NeYcbklptk+15hQ7QMl7adg83ABX7wxloubN/rr6UhByorReheaOi94f1nrYQAw39Ar+5by3huhOmvygL/rr9/aXj3m89frR7tY9ddq4nGdnJtG5+tW39wsrkTBR1sIEUinx16J8+7WjD0przbNvbT6LAAsvErHHcoL/rvQTvbwj8I3W0IIyQBOtoQQkgGcbAkhJAM42RJCSAbsmABZKAD1yQsJ/bH408KyTrAPHf+7qN7Sl+yH9xx3fRZbej+7S/pgNmAGAIMmuLTS9on7q2ZxxFTBB7vONHWFYFvh9vqSz8huQ5tt1CKBt7l2skXlaRMQsws1jtWm3Dan1vR4nzi12/XZtceY65zzhjz2ktoYZGTtCWZfqe/TzNnNCX5lWZG3m+DSVyLb2ABZmsUR3ZxTt0Y5FhsQi+03ZjyThA2iAsCzb/bfzTTwzZYQQjKAky0hhGQAJ1tCCMmAHaPZSgcorF3QJnN1bwrcOaZ1ydwhbx7eauvfTw+evsb1eeX0Kb3foI8VW4xQNIsGCjlvGJOH/iy2n0NlXYTYGs8ca3rd1BaotEUiAb/wYba9y/UZyWvTm2dWtYnH8ao3/lhp6Azxctnr2W2jnUdOGx1j6FxuGh16wt/vYjXZBNoWFowVEkwiKyPutKQxq8nKTCfNsbsZSxpD9NjCh6Ixp2ns8gtoDnxefz+eSjkmvtkSQkgGcLIlhJAM4GRLCCEZsGM0W2DdjOb77HrG63WLRlrL571uao2sh4sR45SWzoett/VlHip4c+7TTZ27d/WAN8mw5ua1mKGN6P2sGteMVw0+67ZpR4pWJhE7ts2rtZpyTIdeWNPjLRXavs/pEdXORwruDZ7U5zB8Uh9r5JO+ZJ3V8E7e4Y2ipx/SOnQ3xi4xw5M0BQptjm832mq3mmgaQ5srYZqeZb5xN+O1mj3g87Px5XRj4pstIYRkACdbQgjJAE62hBCSAZxsCSEkA3ZMgExaQHnuQmL73KsiVWandOK+dHwifNssahja5YNdK00dlBo0AbFa3VeAqJgKBc+uTrg+uyvakOVowwd09lcWdLus283gb/lCRy/mWO4MuD7ThSXVtsEwAFho6/Oy57Da8lUt6nW9n/YJf20wrO9Ved6/I1gjmpUZ0+dtt7ptxr6tK3PYYBgQT2q3JBm7pAlSxfrYz2w14DRj6TYAlWa7LKsIXwm6Mr2JXPPB2eQK1TH4ZksIIRnAyZYQQjKAky0hhGTAjtFs802d6L621+uxrVV9ORoRzTZnEuqX677UZlW0NrmU04n714z4JPeVVnLJzmMrWgOdKK+6PnbhQGzxgWU6r/XYJrzTy2xLG8/Mt3x1XWtCbrXrZ+a8Dt2s6Wse++0/8oS5L74wMg78pTbgae7W44sZyNjlE43rXuP6lJa0GBwzoF6O6MEXc6WMaIpnvVm8HU8vTW+ulF68lbjvS8mG4/KJ96faF99sCSEkAzjZEkJIBnCyJYSQDJAQkg2UtwMjY1eFH7r9Pc+3z97k5er6bq3iDezzGtnqoikSODPv+iyu6j5Tw9qcu972muj0oO4zt+bzTUdKOg94rOzzQvdWtP46bCodThS8IfqUyaGtRApSnjHFHBfbPhf30eV9qj1b07ppTN8+e1YLsLlzXmMuLul3gpz3/sHUw/pDq7XGSGMEniYf1ua/2v2mMZmJkUYDtUYpMV23m2OTddIUm/xi+NQ3QgiHk/r1zZutiFwlIh8Rka+JyKqIBBE52OtxEUJIGvpmsgVwHYCfAnAewFd6PBZCCLks+mmy/dsQwp4Qwv8AoP9zSgghO4q+mWxDCN0tSCaEkC3AjlnU0CkI1iYvBKZGjvvAYHtQ/+6p17xxiilw64JhADBcMZVol7XRizWzAYCXjZ9V7YL4igXn1nTAqZSP9KnrPmsFHXAaLfiFENZ45lzLV849XtcLEmzFYABodHTgz1a1mFvwCyGCWThig2EAUDQxn0j8zlVDiFVOTSIWgLKLBtJUALABsW6CYbH9xBZU3JdQzaEr85WU23XDVl/4kKaCRtVU4MVnP5Vq333zZksIIf3Mtp5sReRuETkiIkdaNZ/yRAghWbGtJ9sQwj0hhMMhhMOFylDyBoQQsknsGM02P1fF+J987UI7osXM/8CkaheOej22vV8bga8u+0R9+5k1rxkZ9osRnljQRuB2AQMA5ETvZz6y8GFPZVm1i6JF5tN17+JizWumisuuz77SomqfqHlTmfPGFH21mWyCM3g0oosbKvP6vAfmvFZtNVp7f2Mabhq9sBvt15JGj01jIBMby5XQaNPo0DGuRGXfXtLt+A7/3Ie7Ot62frMlhJCtQl+92YrIWzf+9+aNf98sIrMAZkMI9/doWIQQkkhfTbbwixn+YOPf+wH8SLZDIYSQ9PTVZBtCJLkzJe3JIZz/8dueb9cm/K4GT+p2zHe7Pa8/bI1HXFGaet+VvToTYvG8D9a1x/Q2sRza5brWN1sRQ5sTVW0YU2/rW/yD4+Yk4U2+qxEj85zRfh9b3OP6VBt6fPMmr7a97C+oVWwj6cWOyue+ntwpBWm0yzRaqt3O5sPGCjXmjWFMLBc3jTF4kqacRnuNcSW01a2kz8boNge59t7XdnU8araEEJIBnGwJISQDONkSQkgGcLIlhJAM6KsA2YvBLmo4f9dtro8NmpUXvFlN9VodEJOiNyML5ndY9awOiA1MeTOYel0HjxbEL6iolLQDS63hA07PLegA2Z5deoHCdyOBrXJBn1Or438Hl3I6cmVNZmIMDtVUe/WUD7zldReMP+kDjhcbCAHxgJPbb4qKBUmBLaB7E5mLKSLZzMRWXACAnAl+xc47tu+L6TZI1c3Ch60eELN0u6hh34e/qtqPpDwe32wJISQDONkSQkgGcLIlhJAM2DGarVTKyF97Qd8afabm+rTLFdXOey8Y5JeNVmnbANpjRnc0XRrHIibae/TBmkW/3+aa/szqvABQKGhtdb6qzWGWV7wWfPVuXSG4nPe66bJd+BDRixeX9bHaS3rJQrHlF5Ks7dW6+HgkR98az8Sq4qbRXy1p+ljSmMrYhQbObBq++m91n9ezR0y7sctf8/u+dPnm4fa805j0bHVTmc2iW5OeGHyzJYSQDOBkSwghGcDJlhBCMkBC8Lmk25HRgX3htmvfcck+VltbPOglbasxNia9c0p+ROtxuRNaC25GzGukYvYT8dwJq0bHzft7J03z+3NYH2t4zOf45o0pebkYMdexxxF/7PNGs83ldJ/a97xWveuoHu/grM9bLlb1Z1bvBLyeac1qYiYuaXJo0+ibljQFH+1+F1416fqMfPIB1a7d+RrXxxa6THNOJD1pNPq/fvSD3wghHE7aF99sCSEkAzjZEkJIBnCyJYSQDOBkSwghGbBjFjVYYonmNiDWtFnlEYoLfvFBZ0X/DmsN6wBP5Xs+Ob12jQ4mFc/5Pq19euFDWPO3L7+qj902ManlWR+kGp7SlSSKkeDcqqkSEQuQdTp6u/qSTtSXSECvtGzOu+oDZNaIpjnk3xFsMOlKVZntxojGbhMzmcFufR+GTvkVNDYQEwsM2oCdPdZ9X0peuHElKunuJNzc8Wi67fhmSwghGcDJlhBCMoCTLSGEZMCO1Wxj+tfgrP7dsxr5XVSZ0+3yYiS5/wa9XX5W/7w+6bcZelJroo3x5MUm+RU/vtaUPq9CSS+WaNX9Nmtr+tjttu/TiJjeWHK2InBD70eaXguuzujP2mWvgdtrPPbtOdfHLi3pxmwlRhqzmiTTltgbjdVjYzGE3P1mzLu9ebg1FLeLO2Ljt9eCeuzl0W11Z77ZEkJIBnCyJYSQDOBkSwghGbBjNdtipCBg65AxjInk2eaNXGj1WQAYPKM1RqtLFla8dlkw/jAh7/u0y1pbLVT9sduD+rNQMR2KXgu2Jt9r4vXZkT36esWMy9sntBGNGP148uuxIpF6PNYoHAC+8plfiGx3ecS0y27yS7sxHI8VarTZxDaXGACGjK57/pC9md4E3xrTNHf7vOrGddrQJo3ZSjemPdvBFCeNni3iv6sx+GZLCCEZwMmWEEIygJMtIYRkACdbQgjJgB0TILv+Fftx75ELgn3UHMQwdNIHk5av0WJ4+bzrgvqY7mP3YwNmAFCb0u3SYmQ8J/TvxqaPfaCwogMt7Y42gylMeMOTVlsHu3IVX6mhVtN9WjUfIMs5oxw9lth5V+b0tYktNrH3KlaJ9oGPv0+105jMWNIEzGKLI5KCZrHxrk4nV6g4d4t+KNp+N5h/mQ6aDZtFDTHTHrso5N4uFzV0W2U2aR/bdZFFX73ZisjVIvIpEVkUkSUR+U8ick2vx0UIIUkkvtmKyH8B8A0A3wLwzRDC8U0fVXwcgwD+BkAdwM9iPWfo3wL4koi8MoRQvdT2hBDSS9LICAMA3g1gHEAQkXkAXwXw0RDCZzZzcIZ/CuBaAC8NIRwFABF5CMCTAN4J4MMZjoUQQi6L1NV1ReQQgFcD+CEAbwLwwwD+MITw7s0bnjr+fQAqIYTXmc/vB4AQwu2X2v7w4cPhyJEjz7djmq018YjpXdbcOtbH6nEdI2/aCr0AkDNSpbSSE6XzXn7F6lVab7XVdiWi2eaMqXeImId3zH5CzSfhizG5CWV9rQpLEZOZObPgwxf/xUO//69VO6bz2YUDdtFKrHqt1XnTkEanjFXytVjjmbkb/XuPXRzTGPH3Jd+49Pc39ozYhSMxY5Ukcx3AX3NrVJ7GBCdGNwsqeomIXJnquht/viOE8EwI4dMhhP8thHArgDsB/LyI+PrKm8MrAHwn8vkjAG7MaAyEENIVaWSEZRF5AsA3L/rvGyGEvxKRPwPwTwB05zl2eUwAiMT+MY91icMhIncDuBsArrmGcTRCSO9IM9n+PICbsC4h/DiAIaxrt8cArAEYEpE3AfhWCOHcZg10g9jfTC/493YI4R4A9wDrMsJmDYoQQpJInGxDCB8D8DEAkHXHhZdhfeJ9NYB/AOAAgL/G+gT8XAjhwCaN9TzW324t44i/8V42Nsfz/PUDrk+nqJWXmMZodbR8Q/+8EikQV5vQvzNiuq41sInlAbeGtC4ajEza6kQMbha1EQ1GfK5rMBpybsj3KU5ovbh+Tl+/1qjP382vaUHbjhcAbvpnOvaZv+s218caslhsQUgAwMd1s9uczySzmpjmePjn9DlNPuqvzcqMvhj1yN9v+bq+Ly3tBRTN1556UH9dqnd6JdCeQ8xMx2q0t779Q6o95g/tiF3fK5G/uxVJk/r181iXDr4TQmgAeGzjvz8Xkf8IYAXA/wzgZqy/AW8Wj2Bdt7XciNT1LQkhpDekkRH+aOPflog8ivWJ9xiAGQB3AfjlEMLTAJ4GsJlLPz4L4HdF5NqN40FEDgJ4HYBf3sTjEkLIiybNZDuGC5rtqwHcAuCnN372KQD/fnOG5viPWM/3/YyI/CrW9dvfAHACF34hEELIliR1nq3aaF27lRCCX9C9iWwszf09rOf5CoD7APyrEMKxpG1HB/aF2659x/PtmKmyxa47jzF80ptdn3qd1toqZ7WuFvM0KBov85hmG3L6s1wkF1fMcJrDJod20I/XeiF0lr3vQX5MC8+Vitds16pa+80XzbGODbltWgNG3677c5r5it5PzD/B5tVanTSNDhjTJXP3fytxuyTNNlbM0Z7DmZt9fKBtHr/mLv9MFJcunY8d81Owz8ihv5h1fez3I3Yd7PVKc63SbJMmz7aXRuX2Wfpi+FSqPNuujGjC+gydeXR/Y6nwP8z6uIQQ8mLpKyMaQgjpVzjZEkJIBnCyJYSQDNgx5uGhVlcmGLbCaIx2KdkMJmpWs3TpBQrWdAYAGqaabqyPtPSx2gPJAZPWoAmQid+mU9OPQXlqzfVp1nWfgZIfYMcsmFg7ZzLsp/w2YgxtJLLo4uxN+tgjx/3Kh7YJZr7+Lb+j2rZSLeCDKjFzojSGLJZYQMxiA0OVQ36hRqGm71V91F+bulk5sLZXx6xL52MVmPV+n37btOtjDWzyN7/W9bnq8zqw1kwR/LILIWJ0Y+KephhAGqOcNIHzon2WUmYQIzjxAAAfTElEQVT5882WEEIygJMtIYRkACdbQgjJgB2j2XYmhrD8o7c+37Ym4ACwNqm1wJgxs9VxF6/1v68aYyZRfy1Z+7X66/Q3/bFXp83iiIiZtD12KJr9NCK/X415eH3eJ9gPTuuqQ0tVv+BDjB6cN4Y2nfPG8AZAKJjxRobXGk5O6bb3xWrpMS3OaoNOi4M3jDmSxjjFGLvEzLlPvVdroJV5f45LB8w5RBY1NHbra5xbiTj5GOyihvq0N8EZfE5PDbGFONaQ3X6nSpFFIvZa1SImOLHtLEkGQda8CPAabeyZsAUEYgtonCZPzZYQQrYOnGwJISQDONkSQkgGdGVE04+MjF0Vfuj29zzfPveDXq6uzOlrsfQSv5+hE1obtMUcAa9v5YwkFjMHsSbP1oQkhh0v4Mds83WtMQ0AYNRWm/R98kWtx5UjRjSry5ETu4jSgN+mcVbn4g6e8JqjNcTec8RrjFZvt2biVosDgKGj2kQ7TTHCbojpks4YfMzr7/beNUb9vq3Wb7eJPZ9iwhWdgr/f5XmT0x25tda8fs/f6kItaa7nlTIPj5kIWaxZUWx8afYz+0od03j437/3yhR8JIQQ8uLhZEsIIRnAyZYQQjKAky0hhGTAjlnU0KoIzl9/4XRjwaXlQ7o9cNoHLaypTCwA0Ro2WePmV9rQM/6y2/2WFpKrMMSCKgOnddsGVcKgDy6hbioGR6rrtpsmwT7vA1mVIV3NYW1JR/kaq5FFDWUdrakeiFUV1se6+D5+n2JVb2eDIY1d3hzGVSN4zHVxxAIo9liW2KKG5ttuVe2Y6dHgrL42pyN1q4um4nJj0jwkpeRiKrLsr6ddQHHg8756sa1kYpP9h+AXiVhiBjL3maCZNRWK8ZXP/ELiNrn7dUAsdi+r+3QkMGb+0y18syWEkAzgZEsIIRnAyZYQQjJgxyxqGNh3dTh013ufb4eIX0drSF8LiVSvbe7SGlhuj9ey2rNmRYIxeinN+YNbs2Y7lth4WkNejxt9Uvepzui2PQ4ArF5jdNzIogYx2mohsqhBzOVqnE9emZFb07/vrS4NxKsIW6xp+szfaQP02KIGa1YzdMpfnCQ9FvDarz1WzGDesjITMTQyent9j9fb7fXLNc0zEjFsd2ZEEVk3v6r77L/fd5q7UWu9k4/q8dmFJoBfbJKGNPfAEjOZsXps7H5bHXrqj77m+px7pzZ6/+//1/u4qIEQQrYKnGwJISQDONkSQkgG7Jg8W2kD5YULWuTqnhTFHHd5nSqMG0PsM16XzO3WWlB7VV/m+nSyTh4rfJir6c9iurPNFQb0sQqrfr+FJb2j1i4vnIZV3adT9NdmcEjrccEUI2zVvG7aMbspLiWbX+cj52CNUqxZyPRDvogloMcT03WLZxOHk6jRxvRYmxccM+fulM29i1yb9m6d29wy+na+7O9lJ28ueohcTxNnOHuTvzZW/7fnPfWgNqYBvOH4yCcf8Ps1Ju62MCfg83Ntccn5N/gCmhOP6+czTfHJ15/0+brNoe5yb/lmSwghGcDJlhBCMoCTLSGEZAAnW0IIyYAdEyDrlICVay7+xAepWoPGZGbEJ5HLSiGxTzGvgxIdazISMWRx7vkDMQMRu0DBC/WVs/qz2m5TvTYSf7Lu/jZRHgDEBv3qfkfLqybK09ZjKS5EqjAM6/O04wf8OdjKFwAwdFz3sQHQWGK8NaexlRuA5AULse1sEChazcMEWWwwDABapgpD7H4HExATE7jsRBaE2KobrYXI82iMh0aOJwd1bXXdc7dMuT6Fmt7PsjHkAYAHPv4+1Y6ZygyZ+7lg9hNbjGCJVYSw46nP+Gd2/MnIA5iCvnmzFZH3isjnROSUiAQR+UCvx0QIIWnpm8kWwD8FsBvAf+71QAgh5HLpJxnhFSGEjogUALyr14MhhJDLoW8m2xBCsgvyJXegjVxsVVIA6AxqrTW/ELk8+7R2WSh6/aZZ19uFNd0u7ln12xgdV1YjCewjRgtu+j9MauaPlbbRAmXIbRLVcS2Fo3qRQMeaVMObl1gd2o4F8As1Vmf8ba7M6v02xv1+Gst6P9YU5bk3T7tthk8a4/LrvMG4xVbFBYDFg3rf+YZdSOL3Y3XcWNVjq6W3JiKmMuYRGBjWz+dg2W8zfz7yEBiC+S7M/4A/78GTxnS+ptsDcxFXIUNpyY/PLljwinLcaOZi7MKIGCfv8M+EpbTs78viwe6mzX6SEQghpG/Z1pOtiNwtIkdE5Ei7Wu31cAghO5ieTLYicsdGRkHSf19+MccJIdwTQjgcQjicH0r+04kQQjaLXmm2XwXw8hT9ImpXd4Sc1mmtngj43ML2hNdjByva+KNR90mUYkw8pKK1q+aCN68pn9G3Imby3Slo3aw9mJybaXU/VxAQwOBxozGPuS7OlGfwhNfwqgf09aqcNvuN/GoPuRS5pOZQ1igcAJojup1GV7Pm0dZcOkZ50V/zFWPQ3inqdsxkxt7f2PMY9uhO0vYXcPfuRdWeW9AHq5S8Jpov6mdARl0XtE0+eXHJH9vqzjEzbovNU47lP7cf04UZa3e+xvWxuc2tQ/o7ZXOdY8T0WEtj5MoVfOzJZBtCWAXweC+OTQghvWBba7aEELJV6JvULxE5DOAgLvyCuFFE3rrx/3+18bZMCCFbkr6ZbAG8G8DPXtT+yY3/AOAQgGNZD4gQQtLSN5NtCOEuAHd1vYMc0B68EIUoLkbc8+1n097df21FB1EGR5KDAtVqxInEUJ9ONrfI1fX4YkELe0NHntHtlWt8YMtWtC3PxaoK27YPLtiAmDWM2fWU2wTWXCcWgFqdtgEovxcbCNzzka+qdszwxAZrWpXkChCx8VliATHL6o2mymykWgJMtY7xyeXE/Y6P6hTHZiuyOKapP5MzPjCYM5vFzH9sQHR2VS98sfcAAOZ+6bWq/ZXP+CoMt779Q6qdJvB25KPvVW27D8Ab5cRUVFtVw1bUAIC2NZZKCTVbQgjJAE62hBCSAZxsCSEkA/pGs33RBCgz65guaU1Q2kteyyoMaXGwujDg+sAYchQntD7XOhfZxhhtWw0S8MYfMWq7dduZaC/5beyiAavhAkChKol9rEFMx/RpL/trbqvexhYWzHxxVrVjJiQuOd4YkcSquJ57p67AGtNjaxP6mscWmyRptLEqzTKv7VU6o14ULQ/rBTSlgr/o80t6ZeTUqL4OI2U/4IU5PeDCTM31CaZqtNWuAWDf3+pr06okG4Pv/22t477+UX/eQ8acxlbOBYBzd+l7F9No3X7NQojoc2TiK92azsTgmy0hhGQAJ1tCCMkATraEEJIBO0azzdeA8UcvaIZLL4nkz41qTWxi76Lrs7Q8qD9oRH5fmY+a88Z4phxxHTHmNZ2IefjKAb3drqP+2LZgotV+rWELsH5tLqbtfXIw9bDW1tYm/fjKC7o9fFJfzzQFFX0upDf1jhlO3/uIz9dMwppUx4o5XvX5OdWO6XyNEa3Bt8wjkp/196n6Uq3HFipeu8wZk55aw49vYpfOqy3m9PU7s+Rv+MAufcPXlvwNl6IxU4po9EsH9HmVF/Q21Rl/3kO33+R3ZDh1q76e++C3sRry1IP6PsWKTVpiJjiAvb9+irTfw7TwzZYQQjKAky0hhGQAJ1tCCMkATraEEJIBOyZA1ikBK9dc1C5Gqpmu6KDP/OmIhb0JZJUnvFlNq2kMWfI6uhBb1DD8tD52rMqsXXxgAzFAxDjFxKQ6kTteMWYbxUi5NhsQq01ETFtMnGXi8cjKDMN9X/oV1f6xV/yK62ODX7E+9rM0AbM0x7aO/zGzGnsf7D1ojEcCKiaw2mr74Fcur7crD0YqhxT1Na61kr/SIib4tey3sYHVwTPJhixVU7Fi4rFIsDNFNQxbGTmGXVxi71NsEYYNbsYWSzSu01UhBmf9OXSK3b2j8s2WEEIygJMtIYRkACdbQgjJgB2j2UpLm89IK1LFtaB1qcaI16nEJJrbxPMY1qw5N95wfapXa5EplkReMcnxMQPvUV2Y1GmrnbLfxi4+WJnxCxasoY3VjwFg5u+8fn0xsQUBFltZFQDelPvJSE9NxyTLp9kmb8xqYjrv4Z/7sGoXahGtv3nphSS1vclJ8Pkhr28XjPFMu+PfjWZXtBGNfarzOX/spZN6oUOkizNojxlmWz1z9Jhuxxa+WC3VPntpsWbhr3/L76j2yCe/7raxVXqHzP0HYgtmvJberlCzJYSQLQsnW0IIyQBOtoQQkgE7RrMtVjvY840LmuJzt/tc12CvRqQI3/CYTkJdPuN1yPyI1n0kZzS9E974o2P0YmuKAwDNiM5ssZpYfVLv15qAA15bixlk2zzL0We84bTF5lQ+8PH3uT5WW7U6KhDXcS02Z/ILnb9Q7TT5uzGdd8DofDGGT+r2/Mv1O4wt1AkAnQGtb1qzbgCom8TlRsTkW4zg2qoZjVEiMYWSEWkj4yst6M9iRTYXr9V9KsaQP/YcpYkPTDyuzzNmGGPvZ8nEA2LPUeVzWseNqcU2z7Y5FDF7mkuO08Tgmy0hhGQAJ1tCCMkATraEEJIBnGwJISQDdkyArFPKqYBNLuJ1UZsxknkkuLA8q4X48rgPWjRWdeXU0NGBg/ZYRJo3xyqcj5iD1JMDEG19aEw+rIMhq9ORqsJmwYKtgAr4gIM1/gCAoVN6QGPf1u75aRYaxIJhNtgV24/tk2a/abDVJc68wVcAsEGfTjm5EnFhSQeGWuPJ5itY9c9E0VTgDU39/hR9jkzefmyBT2zMlqGT+pkdmNMbxYJLlthiCRtYLUYqaNjFBzaIFltAs3iLrsg79eA518eO2T7TALB0IFIdOwV8syWEkAzgZEsIIRnAyZYQQjJgx2i2QeLGzxdTXNA6WnN3xKHDaFmtlk/KFmP6DDHHrfvLXlgzRi+RX4PtAa2R5ev+fGpOUtQ7yje8Dj04qz87f9dtro/V40Y++YAfoKFqFgRUHvN9rBacRltN0meBdPqwra57X4r9vvI9v5fY57o/1Vrg4+/y+nYo6utZHPHaYHNZa5eFea9dWtuUoaO6T206ZkJvzIkiRvpDJ5IX0FjN0967tchz1Kro78vMF2ddnzSGRXYRy8LbblXtWJVm+wzbqs2AX6gBeBfyNObmMfrizVZEbhCR3xeRh0RkRUROichnReRVvR4bIYSkoS8mWwD/AMDfA/CnAO4E8M8BTAN4UERu7uXACCEkDf0iI3wSwH8IITz/946I/A2AYwDeA+BnejQuQghJRV9MtiEElxAXQlgUkScA7E+zDwna+Lky5/usmSKL+fnI5dmntbWhQa+1WZPnpjEPb5T9fsUUmywuec3M5jVWZ/zwbO6tLUZYWvbb2HzYmJZlTTysEXOsjzdi9lidL40em4Y0ubmxgn9J7PmIz0FeNnrhuVu0cL7raCSX9GqTi3t+yPXBHn392hHjdynqZ7Z6wOiJFa9dyjmt60okNNHU/uKYetjrlM+9eVofypz3+J98zW1jnxt7rQCvrUYx+7EabezZs/c7FpuwBviLh7xBUKwgahr6RUZwiMgEgB8AEAm7EELI1qJvJ1sAH8F6FZB/90IdRORuETkiIkeatS5/HRFCyBWgJ5OtiNwhIiHFf19+ge3/DYC3A3h3COHoCx0nhHBPCOFwCOFwsRL5M40QQjKiV5rtVwG8PEW/VfuBiLwLwG8B+NUQwkev9MAIIWQz6MlkG0JYBfD45W4nIj8N4A8AfCiE8JuXs227CKzMXHiRb4z6PsVF/aLfmPRCfbmoAwW5iFnN0qyOSuWXzcKHCR9sKFT1rWiMRYIhbWNEk1wsAcXq5VdYsIGuGF/5zC+4z6x7/he+pNtpFhpsddIE8GylV8AvfKmvJC8ayC/oZ6I97CNZuXnjPGT+Vi3N+q+4XRxTOevH0jTrCqwhzzo6kGpNXGLVEux+SkvJCxhilRpiQdykbe41984uagGARsT0xhKrLpGGvtFsReQnAHwMwB+HEN7f6/EQQsjl0BepXyLyBgCfAPAQgD8RkYtzbeohhMvP4SGEkAzpi8kWwN/H+iLlmwD8nfnZswAOZj0gQgi5HOSiRVnbmoF9V4dDd733+XbMoKM9ojVaWyUXAILRTUM7osQYHbf8lE6Mrr/E66bW9Hnk0ZLrYyvlWjNxwFfBtcbWMaz+6jXHdDpuN3Ruv0m17/uSr4K7WXRT2TeNZmu165i5jj3v+Zf55Pmll+j2rqf8sWyfsqlwa43CgXTP0TW/rhdv2PECXhe1OmqaZyZ2za0RTUxH7co83CxQiC26OPdOvdChOeSvzVWf1+Y5f/3oB78RQjjsOhr6RrMlhJB+hpMtIYRkACdbQgjJgB2j2VZmrg4H7r6g2VptE/CFD2tXe8ErP6Q/KxS8JprL6X2vLWmtKF/227QXjTlI0/8eDGWtMw8e9/FNa1ZzsfkO4E1nAK+1xUw8stRSdwJWFz9/vb+XlXl979rexxrlRd1n6YB+bmKa7fBJWwTUP2uDs8l9ph/Spi2x3FaLLRQaK6hoNdpY4UhrXm8NbmJ5wfbYse/C02/T5jr2+wT4a/7gJ95PzZYQQrYKnGwJISQDONkSQkgG9MuihhdNvgGMHL+gtSz69D40dmuBy+qzADAxqq0ah0oN18eah58wmq1E/BQG9+j9rs56l7LCkl6THdPw7Gere/RYilW/pjy+7v3SxNaVW3PmK2UE3m/Y/N3YdbA5qOMRM/a1yeQ1+FajLS+YHO9F/6zVR3VsorTcXR+r49vzjuXmWq01lmebu1/nJQ9F+thiojbOcPIOrb0Cvrik1XABoLSo27E89ZOvN/flE65LFL7ZEkJIBnCyJYSQDOBkSwghGcDJlhBCMmDHBMjaFR0Ua+zxwa/CoDEGz3uzmmZLi+O1nL+Eq3VtIjM8rgtOrJw3JW8BtGp6P/ldkcBbXQfaSgveJGPlGt0upCi9ZgMFNjgCALe+/UOqPZbCnHmnkiYwmKb6b8kEmGKVXi22unIMG/SJLRpojOjPJh5fc32SAmKxRQ423BQLUhVTLLLxY9YLIfINH9A78wZdyXficW8IVR/V13juRv/9LiSv3YjCN1tCCMkATraEEJIBnGwJISQDdoxmG4oBjX0XtB+J6LGS05+VSr4w4+iA1nlaHf/7au/IsmoX81qpejpiON5o6FthTcpjxBY1FJf0dtXrtN411/ZGzPse0Hrc6vSA6zP27XNmv5cuuAekS+7fCaQpdBlL7s9bzTOi2U4+qp9RawQfW3xitdR7H/F9rEYfI7Zo4WJipunWMMYucojtN6b92iU/1rymXfLfn9Fj+lqdfJ1/zqce1n3Ki/67eupH/NyRBr7ZEkJIBnCyJYSQDOBkSwghGbBjzMPLB68Ke3/tXz7fHpv2OtCBMW3I0uh4I5Dhok5kHC+tuj6doH+HnVrbpdoLNa8VWap1X/DR6sPVJa/h2eKSOSM7t4b8/bY6b0wLtvs58JezvlMCMQ3PslN03W503Ji2mrTf2PW0fayOCnijnJimnGTGHbvfdj/nbplyfaYe1PGBWPFGi9VsY7m5Z27W37vxJ31Mxpr/2IICALA6ozXbZ/8XmocTQsiWgZMtIYRkACdbQgjJAE62hBCSATtmUYPkA0ojF8xdBsve6KWQ04sPrhpccH2mS3rBwmjBB8hW2jpItae8pNpPV31QwAbjvru22/WpFLWg3xzwQYD6S3Q7rOlbnFvzv187Zp1DzOCmqE87aiBiAyTWLX80shDCJuH3I0lBqVgwLE0fGxiKLVBIqnqc5tg/9gq/D1sJIXaf7L59TYNkRp/xZjD2vGOLGmwfG9iKBchspeGVGR8Abw7pZ39tr1/AUFhLXnAUg2+2hBCSAZxsCSEkAzjZEkJIBuwYzbZSbOKle88+3260vV5zYHBetceLXo8dzGmtNw+/SGA4r3Wo0bw2elkq+8UI9Y6+FftHF12fSl5rtqdzI67P+aD1pMExfQ71pr/ljboWbdfKflVDc5fVsiK6bqRybxJWL0yTuN9vpFlYEOuTxgzG6rjFyOKDpG0aXZgKAX6Bgl3EEDOqsQpoTHN2x4qck13EUJvQz+PAnDdcshWC7TYAsHyDWegQ8ZwJue7eUfvizVZERkTk/xWRoyJSFZEFEXlQRP5xr8dGCCFp6Jc32xKAFoAPAjgGoAzgHwH4v0VkOoTwez0cGyGEJNIXk20IYQ7A283HfyUiNwD4OQCcbAkhW5q+mGwvwRzW33JTUZALmYDXjJx3P58s6uqIMyXfpyI6We/q4pzrc7o1dslx5MQLQasdbTwzXfK5hQ8t7lftiQGvKeeNAfpKTV+e8WG/zZmGNsrJT0dyH01ByuJZr4nZXMeWqWuZJlezH7kS5jmx6zCWQn9NIo2JizWdAbzeGjOisRqtvQ6vf8vvuG3sMxA7b2uMEytIWazq59yayiwe9FOb1Wxt3i3g89A7gz57OD/mvx9p6KvJVkQEQB7AKIB/COBHAfyTng6KEEJS0FeTLYB/AeAjG//fBPCeEMKf9XA8hBCSip5kI4jIHSISUvz3ZbPp/wPghwG8GcAfA/iIiLzzEse5W0SOiMiRxoKve08IIVnRqzfbrwJ4eYp+SmAMIcwC+L5r9b0iMgjgd0XkoyEEp8CEEO4BcA8AjL1s985wSSeEbEl6MtmGEFYBPH4FdnUEwM8C2APguUsfU1C7qLLsZCQAtaeoFxLEgl8z+WX3maOgDWxqoXjJNgCM5XXgyprZAMD1w2dVuxn8woxzJR0QKY/pwMHxFZ/APjWmr8Vaw49vyQTImpPe5X72Nt0untfb3PAbPmlk6J23uc8saYJo/V7hIc3ChxhJ1RxiJjM2aLb8tltdH1f1Nlb913xmx+uf4HTnNHRUB6ZjpkdJfaYf8n/JLprqxDZgBgCdUfNcR6pcDw7W3Wdp6ItFDZfgdgArAM4mdSSEkF7SFwGyDV32VgBfxPob7CSAnwLwVgC/HELwfomEELKF6IvJFsDDAN4C4HcBTAA4B+AxAP9TCOG/9nJghBCShh1TXVdEZgE82+txXCGmsP4Lh2wOvL6bz3a6xgdCCNNJnXbMZLudEJEjaUonk+7g9d18duI17vcAGSGE9AWcbAkhJAM42fYn9/R6ANscXt/NZ8ddY2q2hBCSAXyzJYSQDOBk28eIyA0i8vsi8pCIrIjIKRH5rIi8qtdj6zdE5GoR+ZSILIrIkoj8JxG5ptfj2i6IyFtF5NMi8qyIrInId0XkgyLiC+ltUygj9DEi8m4AdwP4UwDfBDAG4BcB3ATgdSGEb/RweH3DhqHRtwHUAfwqgADg3wIYBPDKEEL1EpuTFIjIAwCOA/gM1leB3gTgA1j3SHltCCFSWnF7wcm2jxGRKQBz4aKbKCKjWK/T9rkQws/0amz9hIi8B8CHAbw0hHB047NDAJ4E8IshhA/3cnzbgY1agbPms5/B+ovCG0MIf9ObkWUHZYQ+JoRwLpjfliGERQBPANgf34pE+HEAD3x/ogWAEMIzAP4O68vEyYvETrQb/LeNf3fEs8rJdpshIhMAfgDr3hEkHa8A8J3I548AuDHjsewkbt/4d0c8q5xstx8fASAA/l2vB9JHTADw1T2BeQDeAJi8aERkP4BfB/DFEMKRXo8nCzjZbiFeRLmg72//b7Be8v3dF/9JTFIRC15452jyohGRYawHyloA3tHj4WRGv1gs7hS6KhcEACLyLgC/BeBXQwgfvdID2+acx/rbrWUc8Tde0iUiUgHwWQDXArg9hHDJCivbCU62W4huywWJyE8D+AMAHwoh/OYVH9j25xGs67aWGwE8mvFYti0iUgTwaQCvAXBHCOHhHg8pUygj9Dki8hMAPgbgj0MI7+/1ePqUzwK4VUSu/f4HInIQwOs2fkZeJCKSA/DnAN4I4C0hhAcSNtl2MM+2jxGRNwD4/7D+9vVuABcnhtdDCN/qycD6DBEZwvqihjVcWNTwGwBGsL6owVcHJZeFiPwhgHcB+E0A/8X8+LmdICdwsu1jROQDAP73F/jxsyGEg9mNpr/ZWJr7ewDehPXA2H0A/lUI4Vgvx7VdEJFjAA68wI//jxDCB7IbTW/gZEsIIRlAzZYQQjKAky0hhGQAJ1tCCMkATraEEJIBnGwJISQDONkSQkgGcLIlhJAM4GRLCCEZwMmWEIOIjIlIR0Tu6vVYyPaBky0hnsNYX7L7YK8HQrYPnGwJ8fwwgCV0YXdJyAvByZYQACKSF5FfEpGnsF7GfBeA0yLyQkY/hFwWnGwJWefPAPwCgP8TwCKAjwP4HIAPULslVwJOtmTHIyLvAPA2AHdivZLAONYn2ndivSzOT/RudGS7wMmWEOCfA/irEMLXANy88dm3QghtAMcBDPVsZGTbwMmW7GhEZAbr2Qf3bnx0E4AVAE9utGcAPNWDoZFtBidbstP5ft2xkxv/vhrAt0MIHRE5DGAawF/2ZGRkW8HJlux0qhv/Htr499UAvrnx/78O4Amsl8gh5EXBUuZkp/MQgKMAflFEzmD9TTcvIn8B4A0AbgshNHs5QLI94Jst2dFsBMF+HOsT7sewvnLsZwFUANwSQni4h8Mj2wgWfCRkAxH5JQC/BmDXxiRMyBWDb7aEXOAwgP/OiZZsBpxsCbnAzQCO9HoQZHtCGYEQQjKAb7aEEJIBnGwJISQDONkSQkgGcLIlhJAM4GRLCCEZwMmWEEIygJMtIYRkACdbQgjJgP8fbQJglgoJawEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x18292bf2b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(5, 5))\n",
    "plot_2ala_ramachandran(training_data_traj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "training_data = training_data_traj.xyz.reshape((training_data_traj.xyz.shape[0], \n",
    "                                                3*training_data_traj.xyz.shape[1]))\n",
    "training_data = torch.tensor(training_data)\n",
    "training_data.requires_grad = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Internal Coordinates\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ACE1-H1\n",
      "ACE1-CH3\n",
      "ACE1-H2\n",
      "ACE1-H3\n",
      "ACE1-C\n",
      "ACE1-O\n",
      "ALA2-N\n",
      "ALA2-H\n",
      "ALA2-CA\n",
      "ALA2-HA\n",
      "ALA2-CB\n",
      "ALA2-HB1\n",
      "ALA2-HB2\n",
      "ALA2-HB3\n",
      "ALA2-C\n",
      "ALA2-O\n",
      "NME3-N\n",
      "NME3-H\n",
      "NME3-C\n",
      "NME3-H1\n",
      "NME3-H2\n",
      "NME3-H3\n"
     ]
    }
   ],
   "source": [
    "for atom in ala2_pdb.atoms:\n",
    "    print(atom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "rigid_block = np.array([6, 8, 9, 10, 14])\n",
    "Z = np.array([[0, 1, 4, 6],\n",
    "              [1, 4, 6, 8],\n",
    "              [2, 1, 4, 0],\n",
    "              [3, 1, 4, 0],\n",
    "              [4, 6, 8, 14],\n",
    "              [5, 4, 6, 8],\n",
    "              [7, 6, 8, 4],\n",
    "              [11, 10, 8, 6],\n",
    "              [12, 10, 8, 11],\n",
    "              [13, 10, 8, 11],\n",
    "              [15, 14, 8, 16],\n",
    "              [16, 14, 8, 6],\n",
    "              [17, 16, 14, 15],\n",
    "              [18, 16, 14, 8],\n",
    "              [19, 18, 16, 14],\n",
    "              [20, 18, 16, 19],\n",
    "              [21, 18, 16, 19]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bgtorch.nn.flow.crd_transform.ic import MixedCoordinateTransform\n",
    "mct = MixedCoordinateTransform(rigid_block, Z, X0=training_data, remove_dof=6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Stochastic Flows\n",
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ICEnergy(object):\n",
    "    def __init__(self, crd_transform, molecular_energy, use_jacobian=False):\n",
    "        \"\"\" Single flow layer that does a BD step \"\"\"\n",
    "        self.crd_transform = crd_transform\n",
    "        self.molecular_energy = molecular_energy\n",
    "        self.use_jacobian = use_jacobian\n",
    "\n",
    "    def energy(self, ics):\n",
    "        xyz, dlogp = self.crd_transform(ics)  # evaluate in flow direction\n",
    "        E = self.molecular_energy.energy(xyz) + self.crd_transform.periodic_angle_loss\n",
    "        if self.use_jacobian:\n",
    "            E = E - dlogp\n",
    "        return E\n",
    "\n",
    "    def force(self, ics):\n",
    "        xyz, dlogp = self.crd_transform(ics)  # evaluate in flow direction\n",
    "        E = self.molecular_energy.energy(xyz) + self.crd_transform.periodic_angle_loss\n",
    "        if self.use_jacobian:\n",
    "            E = E - dlogp\n",
    "        F = -torch.autograd.grad(E.sum(), ics)[0]\n",
    "        return F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "class InterpolatedEnergy(object):\n",
    "    def __init__(self, energy_model_0, energy_model_1, lambda_):\n",
    "        \"\"\" Single flow layer that does a BD step \"\"\"\n",
    "        self.energy_model_0 = energy_model_0\n",
    "        self.energy_model_1 = energy_model_1\n",
    "        self.lambda_ = lambda_\n",
    "\n",
    "    def energy(self, X):\n",
    "        return (1.0-self.lambda_) * self.energy_model_0.energy(X) + self.lambda_ * self.energy_model_1.energy(X)\n",
    "\n",
    "    def force(self, X):\n",
    "        return (1.0-self.lambda_) * self.energy_model_0.force(X) + self.lambda_ * self.energy_model_1.force(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Boltzmann Generator\n",
    "-------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RNVP(shift_nhidden, scale_nhidden,\n",
    "         shift_activation=torch.nn.ReLU(), scale_activation=torch.nn.Tanh(), \n",
    "         shift_weight_scale=1.0, scale_weight_scale=1.0, \n",
    "         shift_bias_scale=0.0, scale_bias_scale=0.0, init_downscale=1.0):\n",
    "    \n",
    "    shift_transformation = None\n",
    "    if shift_nhidden is not None:\n",
    "        shift_transformation=DenseNet(shift_nhidden, activation=shift_activation, \n",
    "                                      weight_scale=shift_weight_scale, bias_scale=shift_bias_scale)            \n",
    "    scale_transformation = None\n",
    "    if scale_nhidden is not None:\n",
    "        scale_transformation=DenseNet(scale_nhidden, activation=scale_activation, \n",
    "                                      weight_scale=scale_weight_scale, bias_scale=scale_bias_scale)\n",
    "    transformer = AffineTransformer(shift_transformation=shift_transformation,\n",
    "                                    scale_transformation=scale_transformation,\n",
    "                                    init_downscale=init_downscale) \n",
    "    return CouplingFlow(transformer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "class WrapLeft(Flow):\n",
    "    def __init__(self, leftflow):\n",
    "        super().__init__()\n",
    "        self.leftflow = leftflow\n",
    "\n",
    "    def _forward(self, x1, x2, **kwargs):\n",
    "        y1, dlogp = self.leftflow(x1)\n",
    "        return y1, x2, dlogp\n",
    "    \n",
    "    def _inverse(self, x1, x2, **kwargs):\n",
    "        y1, dlogp = self.leftflow(x1, inverse=True)\n",
    "        return y1, x2, dlogp\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_XV_RNVP_Metropolis(hamiltonian_prior, hamiltonian_target, molecular_energy, crd_transform, \n",
    "                            energy_model_0, energy_model_1,\n",
    "                            n_transform, shift_nhidden, scale_nhidden, dim_x, dim_ic,\n",
    "                            n_relax_z=10, stepsize_z=0.1, \n",
    "                            n_relax_x=50, stepsize_x=0.002, mass=1, gamma=10.0, kT=1.0, \n",
    "                            diffuse_at_0=False, diffuse_intermediate=True, diffuse_at_1=True):\n",
    "    # here we aggregate all layers of the flow\n",
    "    layers = []\n",
    "\n",
    "    # Initial Langevin flow\n",
    "    if diffuse_at_0:\n",
    "        layers.append(MetropolisMCFlow(hamiltonian_prior, nsteps=n_relax_z, stepsize=stepsize_z))\n",
    "    \n",
    "    # split\n",
    "    layers.append(SplitFlow(dim_ic))\n",
    "    \n",
    "    # RealNVP block\n",
    "    for i in range(n_transform):\n",
    "        # ic(x) -> v\n",
    "        layers.append(RNVP(shift_nhidden, scale_nhidden, \n",
    "                           shift_activation=torch.nn.ReLU(), scale_activation=torch.nn.ReLU()))\n",
    "        layers.append(SwapFlow())\n",
    "        # v -> ic(x)\n",
    "        layers.append(RNVP(shift_nhidden[::-1], scale_nhidden[::-1], \n",
    "                           shift_activation=torch.nn.ReLU(), scale_activation=torch.nn.ReLU()))\n",
    "        layers.append(SwapFlow())\n",
    "\n",
    "        # Langevin flow\n",
    "        if diffuse_intermediate:\n",
    "            # Merge, MC, Split\n",
    "            layers.append(MergeFlow(dim_ic))\n",
    "            \n",
    "            lambda_ = i * 1.0/(n_transform-1)\n",
    "            energy_model = InterpolatedEnergy(energy_model_0, energy_model_1, lambda_)\n",
    "            layers.append(MetropolisMCFlow(energy_model, nsteps=n_relax_z, stepsize=stepsize_z))\n",
    "            \n",
    "            layers.append(SplitFlow(dim_ic))\n",
    "        \n",
    "    # add crd transform\n",
    "    layers.append(WrapLeft(crd_transform))\n",
    "\n",
    "    # Final Langevin flow before merge\n",
    "    if diffuse_at_1:\n",
    "        layers.append(LangevinFlow(molecular_energy, nsteps=n_relax_x, stepsize=stepsize_x, \n",
    "                                   mass=mass, gamma=gamma, kT=kT))\n",
    "    \n",
    "    # merge\n",
    "    layers.append(MergeFlow(dim_x))\n",
    "    \n",
    "    # now define the flow as a sequence of all operations stored in layers\n",
    "    flexflow = SequentialFlow(layers)\n",
    "    bg = BoltzmannGenerator(hamiltonian_prior, flexflow, hamiltonian_target)\n",
    "    \n",
    "    return bg    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "dim_x = 66\n",
    "dim_v = 66\n",
    "dim_ic = 60"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reguarlize molecular energy\n",
    "ala2_energy_reg = RegularizedEnergy(ala2_omm_energy, high_energy=500, max_energy=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup full prior and target distribution.\n",
    "hamiltonian_prior = NormalDistribution(dim_ic + dim_x)\n",
    "hamiltonian_target = AugmentedEnergy(ala2_energy_reg, dim_augment=dim_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup Stochastic Flow components\n",
    "energy_model_0 = hamiltonian_prior  # prior model is a Normal on the ICs only\n",
    "energy_model_1 = ICEnergy(mct, ala2_energy_reg, use_jacobian=False)  # posterior model is IC-transformed OpenMM energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "bg = build_XV_RNVP_Metropolis(hamiltonian_prior, hamiltonian_target, ala2_energy_reg, mct, \n",
    "                              energy_model_0, energy_model_1,\n",
    "                              5, [dim_ic, 128, 128, 128, dim_v], [dim_ic, 128, 128, 128, dim_v], \n",
    "                              dim_x=dim_x, dim_ic=dim_ic,\n",
    "                              n_relax_z=20, stepsize_z=0.1, \n",
    "                              n_relax_x=20, stepsize_x=0.002, mass=1, gamma=10.0, kT=1.0, \n",
    "                              diffuse_at_0=False, diffuse_intermediate=True, diffuse_at_1=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "trainer = KLTrainer(bg, train_likelihood=True, train_energy=False,\n",
    "                    optim = torch.optim.Adam(bg.parameters(), lr=1e-3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "batchsize = 128"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add velocities, only 60 dimensions.\n",
    "xv_aligned = torch.cat([training_data, torch.zeros(training_data.shape[0], dim_x).normal_()], dim=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\tNLL: -48.4679\t\n",
      "21\tNLL: -51.5105\t\n",
      "41\tNLL: -52.0141\t\n",
      "61\tNLL: -53.1142\t\n",
      "81\tNLL: -54.4816\t\n",
      "101\tNLL: -54.9355\t\n",
      "121\tNLL: -54.5678\t\n",
      "141\tNLL: -54.5307\t\n",
      "161\tNLL: -54.0902\t\n",
      "181\tNLL: -53.5582\t\n",
      "201\tNLL: -53.9053\t\n",
      "221\tNLL: -53.6568\t\n",
      "241\tNLL: -53.6560\t\n",
      "261\tNLL: -54.2466\t\n",
      "281\tNLL: -53.6887\t\n",
      "301\tNLL: -54.0896\t\n",
      "321\tNLL: -54.2746\t\n",
      "341\tNLL: -53.7516\t\n",
      "361\tNLL: -57.5751\t\n",
      "381\tNLL: -54.1763\t\n",
      "401\tNLL: -52.9791\t\n",
      "421\tNLL: -54.1316\t\n",
      "441\tNLL: -56.0911\t\n",
      "461\tNLL: -53.3274\t\n",
      "481\tNLL: -55.8705\t\n",
      "501\tNLL: -53.0071\t\n",
      "521\tNLL: -55.3959\t\n",
      "541\tNLL: -56.9035\t\n",
      "561\tNLL: -55.4860\t\n",
      "581\tNLL: -53.9594\t\n",
      "601\tNLL: -55.3819\t\n",
      "621\tNLL: -54.9170\t\n",
      "641\tNLL: -55.2599\t\n",
      "661\tNLL: -55.5864\t\n",
      "681\tNLL: -56.4897\t\n",
      "701\tNLL: -54.8998\t\n",
      "721\tNLL: -54.2062\t\n",
      "741\tNLL: -56.4618\t\n",
      "761\tNLL: -54.6318\t\n",
      "781\tNLL: -56.6995\t\n",
      "801\tNLL: -55.7290\t\n",
      "821\tNLL: -55.9380\t\n",
      "841\tNLL: -54.1841\t\n",
      "861\tNLL: -55.2155\t\n",
      "881\tNLL: -54.7161\t\n",
      "901\tNLL: -55.6237\t\n",
      "921\tNLL: -56.3180\t\n",
      "941\tNLL: -55.9921\t\n",
      "961\tNLL: -56.2976\t\n",
      "981\tNLL: -55.3375\t\n"
     ]
    }
   ],
   "source": [
    "trainer.train(1000, data=xv_aligned, batchsize=batchsize, w_likelihood=1, w_energy=0, n_print=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evaluate\n",
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,0,'Iterations')"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEWCAYAAAC9qEq5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJztnXe8FNX1wL/nFR5ViqAIgohiATtoLCiCYu8t9h5rTIyJ0diixqjRJPYakxg1P5OosRcEFbCAihpAFAQpgtI7PN7jlfv7Y2b2zc7O7s7uzu7O7jvfz+d93nuzd++cO+Wee84991wxxqAoiqIohaCi2AIoiqIorQdVOoqiKErBUKWjKIqiFAxVOoqiKErBUKWjKIqiFAxVOoqiKErBiLzSEZFOIvIfEZklIutFZJWIfCwiZ/qUnSsixufnuGLIriiKosRTVWwBAtAGaATuAOYCNcCPgadFpIcx5h5P+VHAzZ5jM/Iso6IoihIAKdXFoSIyAehojNnZdWwu8IExJsEKUhRFUYpP5N1rKVgONBRbCEVRFCU4JaN0xKJKRDYVkYuAQ4F7fYoeLSK1IlIvIhN1PkdRFCU6lMKcjsPlwAP23w3Az40xT3nKvAp8CswBNgd+CrwoImcZY55Jd4Lu3bubfv36hSexoihKK+Czzz5bZozpEaRswed0RORgYHSAouOMMQe6vtcD2AroDhwDXAxcZox5LMW5KoGJQE9jTJ8kZS4CLgLo27fv4Hnz5gVsiaIoigIgIp8ZY4YEKlsEpdMe6BugaK0x5rsU9TwJnAh0M8YkndsRkV8DfwB6GWMWpjrhkCFDzKRJkwKIpiiKojhkonQK7l4zxtQC00OoahJwDpYbbUGKcuKcOoRzKoqiKDlQMoEEPgwD1gFLkhUQkSrgZOA7Y8yiQgmmKIqi+BP5QAIRuRjYGxiDZdFsCpwCnARca4zZaJc7DTgWeAOYj2UBXQ4MBk4rvOSKoiiKl8grHWAqljL5I9ANWAZ8DRxljHndVW4OsBlwt12uFiuS7TBjzKiCSqwoiqL4EnmlY4z5CDgiQLmJwIj8S6QoiqJkSynP6SiKoiglhiodJWuamg3/+XQ+Tc0aGKgoSjBU6ShZ8/SEufz6hSk8M1EX1CqKEgxVOkrWrKi11uSurN1YZEkURSkVVOkoiqIoBUOVjqIoilIwVOkoiqIoBUOVjqIoilIwVOkoiqIoBUOVjqIoilIwVOko2VPgvZgURSl9VOkoiqIoBUOVjpI9IunLKIqiuFCloyiKohQMVTqKAtQ3NrHvHe/wzteLiy2KopQ1qnQUBVi0uo4fVtdxy6tfFVsURSlrVOkoiqIoBUOVjqIoilIwVOkoiqIoBUOVjqIoilIwVOkoiqIoBUOVjpI9mgZHUZQMUaWjKIqiFAxVOoqiKErBUKWjKIrSyjDGsKauoSjnVqWjKIrSynh03Gx2ufltFq2uK/i5VekoiqK0Mt6atgiAhas3FPzcqnSU7NGtDRRFyRBVOoqiKErBUKWjKIqiFAxVOoqiKErBUKWjKIqiFAxVOkr2aBocRVEypOSUjoicJiJGRBYk+fwnIjJdROpFZIaIXFJoGRVFsVhf38j6+sZiixEZ5ixbz38mzS+2GEWlqtgCZIKIdAHuARYl+fwnwGPAHcAY4CDgYRERY8wjBRNUURQABv12FABz7zyyyJJEg2Me+IC19Y2cMqRPsUUpGqVm6dwFTAZGeT8QkSrg98DTxpjrjTHvGWNuAJ4Efici1QWVVClJDOoyLEduffUrjnrg/WKLwVq1+kpH6YjIfsCZwOVJiuwD9ACe8Rx/GtgUGJo/6RRFiTJ/+3AOX36/pthiKJSI0rGtlMeBu40xs5IUG2T//tJzfJr9e2A+ZFMURVGCUxJKB7gGqMGaq0lGN/v3Ss/xFZ7P4xCRi0RkkohMWrp0aW5SKoqiKCkpuNIRkYPt6LN0P2Pt8tsC1wM/NcakSonqJALLyClvjHncGDPEGDOkR48eWbVJURRFCUYxotc+AnYMUK7W/n0/8C4w0Y5eA2gDiP1/vTFmA/EWzUJXPY6FswIlXDThp6JkhTEGyfL9WbKmjvEzl3HS4C1DlqowFFzpGGNqgekZfGUgsBWJbjPsY/cBV9IydzOIeKXjzOV8lZmkSikxf0Utvbu0o6JCFaFS3pzz90/5euEaRuywGd06tCm2OBlTCnM6pwLDPT+jgGX23w/a5SbYx87wfP9MLCvnw0IIqxSe2UvXsf9d7/Hge8liTBQlWuSSzGPpWmuWobG5OSRpCkvkF4caYyZ6j4nIuVhutbGucg0iciPWYtDvsRaHjgDOB64wxmwsjMStiIikwVlo7344cfZyfnbQgCJLoyhKKkrB0gmMMeZR4FLgFCxr6DSsAISHiipYBGloaubtaYswEVEcitKaaM1vXUkqHWPMucYY31k0Y8xjxpjtjDE1xpgBxpiHCy1fKfDAu7O46OnPGDtDw8QVRSkcJal0lNxZsNIKDly+Xr2OilJoWrOHQZWOoihKgWm9KkeVjqIoilJAVOkoACxeU8eNL31JY1NphmEqSimRm3ettNeiqdJRALj+xak8PXEe42e2zsCCVuxiV5SCokpHAaCpufR7XVUcSqnQmvdtUqWjKLTuiV2l8BR9gFREAVTptFbCeOYilvAzYuIoiuKDKh0le4o+XAuP1rxuQlEKiSodRVEUpWCo0mnlqEfKQu0cpZCEYliX6EOrSqfAjJ2xhBF/HMvGxmishynR51ZRlBJFlU6BueGlL5m9bD2L16Taebv10NRseOCdmayrbyyqHGFN6XwyZwUrNZ9dqKypayi7a6oh00qrpdjutdem/MCfRn/DXW9lsplsNGluNpzy2ATO+tvHxRalrNjj1tHs/rvRxRYjVFpz3IoqnVZOsZ/9etvNWLuxqciS5H4lnBq++mFNznUpLTSWwcJlpQVVOoqiKAUmFzVa6uvRVOmUIMYYGkJKzBmV5zcMd0MudYRz/vIakdc3NoX2nCmKgyqdEuTZT+Yz4Po3+WHVhpzrMp7fSvaU2zXc/oa3OOSe8cUWoywptwFKJqjSKUFe/t/3AMxdvr7IkoRHGC6DXOpovV1AauYsK59nLEq05udNlU4JIyE4x8Tzu1hylMPArxzaoCSnsamZ2o3FDe13U6qPmyqdEqYcYv3Lak7J53488f5s5q+ozb1ypeic9+SnDLxpVCh1teYBiiqdEqTUo1daC8vX1XPb619z1l913U458P7MZcUWIY5SVVyqdEqQqD1suVhcUWlKGFaj9744y0uKnW1BiSB5sqwDU8SRqyqdEiaMuRQl/0RtkKAoxUSVjlJUwlSbxV+nk3sdAC9+sYCZi9eGU5kSSfJhWZcKVcUWQMmcEn3WWi2ZejJ+8e/JAMy988g8SKNEgVJVGGGglo4SCcIY+eW0TqfYPnZFyZBSfdpU6ZQgOpPjT7FHj8nOX2y5yo0Xv1jAfybNL7YYOdGaHwl1r5Ug5fjAFjsoQq2U0sFxP54ypE+RJSkupZpKRy2dEkbX68RT7OuRrAsotlxK9MhFYTiPU4nqHFU6rZUwn9dQ0uAUOZonn1mmS7VzUPJHa34kVOmUInl4Yov1EkgZmwFl3DRFyZqSUzoicpqIGBFZ4PPZWPsz78+VxZA13+SWpDMalKpf2g9vS8qoaUrItOZno6QCCUSkC3APsChFsSnAxZ5jc/MlU7aEE6Ib3neLrYSKHUgQBq25I1EKT6k+byWldIC7gMnAQuDgJGXWGmMmFk6k7Mgtb1J4ckSFclin40Xda0oyQpnDDKWOwlMy7jUR2Q84E7i82LKEQU6dXMRGOLk8/GHO6RQ9lU3E7osSYVrxs1ISSkdEqoHHgbuNMbPSFN9dRFaLSIOITBGRCwogYsaE8cyF0V3raNwiKiPPsLjwH5N44v3ZxRZDySNFH2RlSUkoHeAaoAa4I0258cCVwDHAScBM4AkRuSG/4pUupeoX9iNqCrSY13bM14u57fWviyeAkpIyeu0ypuBzOiJyMDA6QNFxxpgDRWRb4HrgeGNMXaovGGNu8hx6WUReBK4XkXuNMet85LkIuAigb9++gdoQBrlEbUVpRB0lir9Ox/O/3iclj4TzdBX+GS1GIMFHwI4Byjl7/N4PvAtMtKPXANoAYv9fb4zZkKKeZ4HjgJ2BCd4PjTGPY7nuGDJkSMHuQCjutRCG9lGzDkoZVTFKUMrJw5ApBVc6xphaYHoGXxkIbAWs9PlsJXAflkstGbGsERmcM+/kcwV8oeUoB8K4DAn3Q6+tkkdK9f0vhZDpU4G2nmPXAoOBk4GERaIeTgc2AFPDFy0XyqdHCmWNTRlGkJfPHVbCJreIT6eO0iTySsdvzY2InIvlVhvrOrY/ljL6L9Zi0M7AOVhBBdcaY9YXQNyCUg7utTBOnxcrJUQ5SrVzUPJHVDwMxRAjNKVjBwj82RizS1h1ZshCrGi8W4HuQANWdoLTjTHPFkmmpETloSs25XQZErxr5dQ4JXKU6vMVpqXTGRgUYn1JMcac63NsFnB4Ic4fBjmtQSzRhy3K5OOSOi6UTKy5cspFpyQnKndZ1+m0IsK42WG6xoqWZTrMuoq+n078VczmHqvOUYJTmg+LKp0ioWs4wqfY63SS3VK904qXqFi0xZBDlY4CRDcCrJTJ5nWORlek5JtSTWETBmnndESkf8C6euYoS6siag9Mbnkuo9GY3Nxr4UevOaNIndNRwsRZopBH4zyvBAkkmEUw2SRgOYXoKZ1yoNjXNNn5MxFLHwulIBTxZQmidM7LuxStkKhYBw6l7F6LqqtCAwnKE2NMJLZZj+pzn460SscY849CCNLayGnSOzwxIkM5tClqAwklPxiTe6Rkax5chBZIICLHisiasOpT0lP8sZZFLmlwwhgwRmUHxXC2IG/FvZGSEaW6B1SY0WtVQIcQ61PSUA7dU6m6CIKg7rXyJJwJ/NJM1hkGGjJdJEr1gckXxbbawrFSvP/rTY4a4WRmLiOFoRkJWg9hdEjF7qjDJLeQ7Wig4c5KISnVx02VToFxHpR8jKwzkyMaT2woczrOepgc6srH6FXda9Gj2O9dmHWEQSTX6YjIiIB17ZyjLEpAoqIwokbULkt2GQki1gglgag8Z6X6rARZpzMG6/1JNY50Pi/Nq1BAwtyAKSoPf7GJ+sgzk0GC3tP8EpXLG5WBYyTX6QDD8y5FKyQqE5oRWONWNiS61zK/P9HoipRUlFOIfjEIsjh0XCEEaW1E5Xkp1Qc3joiGXTtVRmH1umJhDQRyux/FfmfC3dKk8I0JMqdzUyYVGmNuzV6c8ifMBzYq+qLYvuVin98hmRyZudei0RYlv7Tm2xzEvXZzgDLuS6hKJwBRWRQZ5hxTLhS7s83HYr2sotdylkJJRTld36j0IZkSJGS6Os3PnsDbWDbrrPyIWT60mMbRSGERRkeZSxqcMIjuqNEJ5Q5+faLblvIgOh11673RaZWOMabJ7wfoDzwDfAwMBC6yfysBiFrCT6dbLFanl8u8Rz53/cylipY1Wa23g1H8iUqevkiu0/EiIn2A3wJnAyuBXwEPG2M2hixbWROVwZLOcYdHKMpF9VNeKdUkmeVEYKUjIpsB12NZNHVYczf3GGPW50m2siRqgQQJ7rUQ6syGXDrsiBg6odSpHVr0icrYIhxXYTSj1zoD1wBXYHlh7gP+YIxZmWfZypqoelwK/RCWk6WVzL2mczrRoZyub6k2JYilMwfojBUscBuwEOgqIl39ChtjZocnXmnw70+/Y319E+cP3Tpt2Vi0WDllqs2BqIzW8jnBrHM65UVUrJQwiOqcThf796HAIQHKV2YvTmny9rTFLFlbH0jpOITjEspD9FrONRaeqMqczf2JaluUFqIyiIiKHJkSROmcl3cpSpzqygoampoz+k5O0Wt5fNYK/RxHxb0WzlxKfB3ZZZkuzY5EyYzIzN1FMfeaMeYfhRCklKmuqmBjQKUTaiBBRCY0i01U+mm9H9EnKvcoKnIUA91PJwSqK4XZS9dz22tfBR6pRiXGPlHeUnyUIzKnE0KdUVGg5Uo5bBMd5pq6YlhcqnRCoE2ldRmf+GBOWosn5k6K2OLQWN3a6YVGZFwoSuRoze+ZKp0QqK5suYxNzUEtndwJJWqrSN+N1RGRtCSRaYsqqrwSlewV4VCaVpsqnRBwK52GpoBKJ4ebHYKx5BIk5b95p5xGfF6FkVXbyuh6lCua1SA3VOmEQHVVSwhWYxr3WqidbJGzG4QZeFZ8iysanUA0pChfonJ9I2Phq6Xjj4jMFRHj83OcT9mfiMh0EakXkRkickm+5WvjsnQaA7vXipv2JZkchX4Io+LWCoOki0MLK4aSZ6LyvJUqGSf8LCKjSNzbZ4b7HxH5CfAYcAcwBjgIeFhExBjzSL4Ec7vXNjYGCySIiv8/lqolxDozO3803uC85LErUvRaVK5pFCn2PGiYlOr0VCkpnWXGmInJPhSRKuD3wNPGmOvtw++JSC/gdyLyhDGmIR+CuZXO6K8WB8pMUOwRflQCpZ3z5uKqi6p/3JErk7ZFtS1RwBijW3+7KNWxRUm41wKyD9ADa48fN08DmwJD83Xi6sqWF+HW174K9J2ojEbDWFsSBjnN6URkyJd0u+pM6ihRP30hiE6UYmlGjflRjH6olJTO0SJSa8/VTPSZzxlk//7Sc3ya/TtvG8y1qQp+GaOSkcC7wCzUiLhMiMjL5xBmeqJibVcdsUsaGsX2DsTqyL2KUIjKwDVTSkXpvIq1tcKhwBlY+/m8KCJnusp0s397t1xY4fk8dNzutaAUu3OJyuMamcwMebgiRXNZlmhnlI5yaldU3KitYk5HRA4GRgcoOs4YcyCAMeYKTx0vAhOxAgYcd1pWg3URuQhrYzr69u2byVdjZKJ0wshIEKa5lBi9Fo2XIRNKUeZkRKUtUZEjdCLiviwniytTihFI8BGwY4Bytck+MMY0ichzwB9EZAtjzELiLZqFruKOhbMCH4wxjwOPAwwZMiSr++ie0wH46oc19O/RgbbVyXd5CCfyrFQfuxai0oR8dCTZ3J+odEZRuS9uoiJSVKyUMCjGfS640jHG1ALTQ6jKazM4czeDiFc6zlxOsBn+LGjjsXSOuP99jtxlCx46fY+446s3NLBg5QagPDqXYrsIo07x3GtFOnGeicoygzAo9rtbTEplTicOOzz6ZOA7Y8wi+/AEYBnWnI+bM7GsnA/zJY+fe23S3ETD6qaXvTEO2ZHPpAYFXxwakRcnHx2aN0ijlIjIbYkjKgqj2GI4YeP52AOqEER+nY6InAYcC7wBzAc2By4HBgOnOeWMMQ0iciPWYtDvsRaHjgDOB64wxmzMl4xBo9eWr2sRISp+4YQ6i/RGRXVTu0wIIyNBVEbz5eC69SM6HobyvL5BiLzSAeYAmwF3Y83P1AKfAocZY0a5CxpjHhURA/wSuBr4DvipMebhfAoYNJDAve1BLo9cy8g5opork9NHJHlimGNGySFapBz2e8kX5dSuUJpSoq9/5JWOnYVgRAblH8NKhVMw2lT5O1BmLl7LRU9/xguX7ku3Dm3ikoFGZYFZLu61MNxGMRdUKfqgkuDc21LexC0iYoROa55LiQolOacTNZJZOiPvGc+cZesZO2MJEL/tQdTW2BRt0tv5XWT3WjiDgGhYXGEQxY41Mvc5ItZoac7oqNIJhXTutZoqK3S6Ic22B5mSj1FbFDubdERFZq97TaPXypOoXN+oyJEpqnRCIJ3ScQINGuLca7mfN8yRX/HS4ETjzQk1IjAn91pERuKRsblaiMrcXThEQxLdT6dE8a7T8VITUzruO5z93c5P1FpxKFUXgR9RcZmEQUTGAnFERaaIiBHJgUEQVOmEQLVPIIH7BamqsD4P3dLJw8gvq1X0JfrwuwmnQ/Ou08niWkbEAo4iUQkCiEoQUBgU491VpRMCm3VqyylDtow71uTaQXT9xiaufm4yC1fXxY7lFDId5kZwngWMhX4Eo9IJ5INSth6jSFTvczZER4HmXkemqNIJgcoK4a6Tdo075l6T88bUhTz32YLQzheVdy8qi+TCuRz5G71m1sSIRNFF5BkLm6hEnrVmVOmEyJ79usb+drvSPvp2WULZqPj/E+rI89qSJWvq4tcrZX66yJLLmqdcvpMPougyjZ5E2RPZ978AqNIJkX226R77u66hpWNdvKY+oWxUUpV468hnZ7O6toG9bn+H217/OlGOXCqOqJshtl11BgtfS7UjKQRR2SE2KpSqu1GVTohk8hAELWqM4YvvVsbVXehR6OraBh58d2bcPJWboNKsqWsAYMzXi1u+G5GAijAIw70WFT99JPuzqFgHZTQPqdtVlziZ3L/J81cFKvfqlIUc//BHvPy/H7KUylIan83zbqjqj18b7nxrOn98+xtGf7U4bdlM644tqMysqtDJy/xUNPqVrChh0fNOVNYLleo9UqUTIpk8jE98MMc3Q8GStXVxndfspeus38vWJ54v4OnO+fsnnPjIR76WSib95NJ1iW7CjASx8XM35fICRWTQmEA2YoVitUX0euSKBgGUB6p0QiSJ9ykp39sbujl8t7yWvX7/DsP/OJYvvlvJho1NcXNDDkJm+2lMWWBZVY3N/ml4bnl1Gu/NWJr0+13bVwPwzaK1sTxy2eAnr6NgV9Vu5NFx3xbN7RDmxG6zgVcnZ2eZRsXdGBX3j5uoiBQZ92VErkemqNIJkeYMn6T5K2t9/5+7vJbjH/6IHW96i0fHfRtXZvL8VUz9fnVG53FygflaOhj+/uHclv99mrBJO0vpPD1xHuf+/VPG2G42p3PL9NkXH2fa+zOXceeb05kwe3mGtUXn3XNfuyue/aJo0WtRiYwKm3D66Wi4xqKCrtMpcTK9gWs2NMb9n2yi3s2xD7VsgBr0fE4X3+iqf2Njc+A6vCriwqcmBTtxALznj08VlF0dWcmRh3xlUQlwyIaoWBVuomh9ZUu55scLgiqdELlw6NYZlV9f71E6rpfKm8/Nb6I903fwqPs/AOCzeSt488tFvmX8HuSmJCdyDgePxEtfJl1AgTGGfte+zhPvzw52Up/vv/jFAuoamrL6fmbnyuI7OhLPK5GJPMuhilAzkmganNJms03aZlR+nUvpzF9RS7PLEqkI8c44D+l3Kyz33VjX/I3XuPKNMAvpuZyxeG2cPJD40Dufra1roN+1r/Pm1IVxnzvWmnutTyYvzoezlvOLf0/mjjfi1wpFxcdeDp1ivoiKSFGRo1RRpRMyj581OHBZx9IZ/dVi9r/rvYSQ5HQEffi9cyhrNjTE/k63x8/MxWu5e9SMnM7vcPHTnyXWkaSSOXa03oPvzYo73uhyv/1v/ir+/uGcjDpqZ63QkrVJIvFS8OrkH7j8/z5P+rlXjFLunMIORohKap6o3JPorBfKvY5MUaUTMocM6klne+I9Hes2WkrnSzswYFLAtTQO7hf5wXdn8vcP5/DalB/od+3rLFpdR0NTc2zS383auhYLy+tm8j6DzwfIGbdwdR0LPEERDlMXrOazeSvijqVyoTkK0rHAKiviS7tz2h330Ifc8upXsf+XrK1LcFl6cYI9vGHbQd69K579gtenxFtetRuTny9fcxB1DU2c8tgEpi7wDyiJSscaNlGJyiujqaWioEonD1QEXOnodJDO1gepOkwR4nKWuVmxfiN/fPsbbnn1q5iS+Grhau4bM5MLn5oU11EDrHEpnQ3euQ3PG+X9rl/RFz5fwNA/vOdb5ugHP+DERyawurbB9/NE68A60qIc4i+m3zVw6vh26XpOenRCUnndMq/Z0Mjdo6YHCt5IrMMwcfZy7n9nJgNvGhVTqvkaedY3NjHRFdX35fer+WTOCm5+dVrgOgA+mbOCkX8eF2g+K+y2lFdHHQ3FFRVrKVNU6eSBioDJturtNTiVlVb5dS5l4PcwbHv9m3H/1zU0sWh1HT//1xexY04fWt/QzLf2wtKE8zY2uf5OrlQuempSXDh1Lvzq+cm+x73tbGw2GGOYac//eBV4uui2rxeuSfm58+0PZi3jofe+5R07JU+yEfD8FbV8MDM+YWuzgVMfn8ifR38DwBffrbLrjq9j2g+pZfGXL1GO21//mlMfnxhrm3MNqiv9n7NkFsGtr01j5pJ1zFzs/1zE1xEuodQXsJKvF67h8fHf+n5WTsqvVKP5VOnkgQpXT3nC7r154dJ9+O3RA2PbVjs899kCZi1ZG7N01rosnSBrfm58eRp73/FObP4DYPw3VpDAxqZm3/kaY0zci+cd9a6pa+QnT01iydo63k4zx5SJu2NpkjkUbx2NTYanJ87jmhemAlDpUeB+bVqWwfyM90VNZ+jsf9d7nPnXj1PK4HUBOiSbC0stX+IxJwBj5fqNceevrqxgw8am2DxV4HMUyAFXLMV19AMfcPsb00M+u0uOEo0a86MYUqjSyQPuLujOE3dh8FbdOG+/rRl15QEJZY9/+CNfy8g7ovdbUOmw0cdaqW9o9rVimppNnEvpv59/H/f5c5PmM/qrxTz8nv9IMQy8LrN4+Zr5nysvnXeeq9FHS9z3zszA5/Z2GFVBfaEuvJ28U0fQzuj7VRsSrKeYfD7HvPfeUTptKisY/sex7HLz2+kr8aknFZmMor9euIbfvfZVyu8Uci7FeUaas3CdBpIjjDpCqMSJRi01VOnkgYN23Cz2t9sF4te/ra1rTDlvEgS/ke6GhiZfZbS2rjHpuhtoGbWni2p7dNy3vm/f8nX1PDVhLkvW1sVvz52kniCLQ93uwHRygaU4k+EdYSazUlKx1+/f8dRRYdcdjMPuGZ9gPQXBqd9t6SxaU5e0XNJ6AgiaSZ94+l8m8tcP5rDCtsRazhNup5+pdeA3QCkn/LYIcWOM4f53ZrJodeIz4i5TaKoKfsZWwK3H7sTAXp35fN7KuFF9srkev/xqXlK52/y+X9fQ5KvMTn/iYzq0qUxaV7W9KLUxzdzJnW9O5+ID+scdM8Yw+LYxALw2ZSG9u7Rzfxj7U7BGoY3NJsEF1djcnDAib2wy1NhPahClc/XzUzhxjy3j3Jw+YsSEeWrCXGqqUo+//vBWcndNptbS2hQBI36dgPex2ejM6SSROdnE3rIiAAAgAElEQVSjEltUmF7EjEbizjOe6ithdG2Z6pAgCW6zISpBAOmYsXgtfx79DeO+WcoLl+5bgDMGQy2dPFBdWcFZe2/FPT/eLe64n86pEKgPEE2UaZTVvWNmxia43Xy9cE1KS2f6Imv+4Ks0E/KQ+OK4rZTVtQ285cp6MHd5iytg9rL19L/uDd9J/42NzbzweXyYtlvRBE2T88e3/edTvJfxna8Xc9PL0/jj29+krO+RsSncjbEV4v6yJTve3GzY+/Z3OMtl9aTsuO0PGxodSyczZZev7SOcevO9s0Oy6M2k5X0S3D6aJMAgE6IyH5MO5364A5QSyhRIFjeqdAqIn6UjItRubPIcS/xuQ5IM0Q6bb1IT939CKLSLIL7uTJOKes9ZXSVx/6/ekOgC9Fug+c+Pv0s45lY0QTueUdMshdfcbLhn9Dcssd1QXgUwxV7r4gQ6fL9qQ8YuB0cpZvoCj/tmKYvW1PG+a37H79QtFoqJO583VRJYSnRZsi0oYuewUgnd+NKXyctk0BqvfH58MmdFzq6cTN1lfgO1//N5vopBIdxajus41SCzGKjSKSBupdO9Y419DFbWtvjC+3Rrx+hfDEv4bjp3V4c2wT2lYT2EazyKxN2Jffl9ektp3vLEPYKm+Cx4bGhqjnUgmc5/fbt0Hfe9M5O9bn+Hfte+njD56ne+sTOWYoyJm0tKxYaNTUnD0yF54MR5T34aqH7H3Thp7kq+Xboubk7HzcbGZi74xyTOeCLJfJFHjqcnzkt+0oweEX+/nfsxO+OJj/lPirm2IHiVSLqOO19zOtm8Pt8uXRc392qAe8d8w/wsggGCbn/ueH1TRsLqOp3yxnkIKqTlhRERVrkWTs5fsYEONYlzLtMXpe7Ee3SqSfm5mxzjFmL869OWTqRTTRWvZLiHTNC0P/eNmck2173BgpW1/OzZL9J/gZZ36VVPBoH3k0SNuTnvyU959pP5gc/1+PjZHPSncbE1PwmyuF769Aszk/cC970zk4P+NC4WlehVOhtsi9nPqkx1hh9WbeC6F6fGRQ0m48RHPkqwFiTWuXnPE3/A7WLNBu98XrrOP4hLesX6jYGsji9dln82SuegP41j79tbAlDmLlvPvWNmcpFPaqhMSCV7bK4tWoaOKp1C4jwElRUSszY2NjYzc8nauHJtqxKVzoezUu8zc8TOWwSWY8X6zPOOpSObUeXHc1akLwT82x4hf7t0PcvWbUxT2sbAG1MXcr8nnDropP8X361k1LRgStFxEz4zMb3r5sp//S/l56ncaw6O4qquiv8glUvVjbcz3vfOd/m/j7/jxEc+apEjyXc/m7eS616cGi+fU2+a3i3XOaVkSuSFzxb4DsrSPZOzlqxjj9+N5pkALrejHvggmJA+OIrB7UZ32hI023lTs/HNWJLqkjvXe86y9bw7Pf5Znmxb+Jplusxx+jsRiXOHLV7TogQGb9WVttXJo8uSUVUpHLlLMMXjPl9YBIkqy5Wm5ub4iLgUGOCV/yVaXlUBJ9+fC5BzLihr6hrpd+3r9Lv2dd6a5r+lhEOQLsCJVvRG+aVTOk5pv1B6sBbiNjcbmpsNlzwTfAQes3QS3F/+5bKlsdnE5bozWJbKL5+bzPEPfZRQvimJS/o9e/dbZyv4cT674dY3NiV1r17zwhR+9VxLho3/TJrPq5N/YNaStb7LF8Jw893w0lQG/XaUz+Lm5HW7T/vq5BaLP1/rl4KiSqeAOHM6lSI8frZ/NuqHz9gjbfiuH9UVFb4Ty4WiEGsi/vfdKvr36BCo7Jxl6307+GzW5RSSIK4QZ3Ts7nA+nbvCN/no9EVr+HrhGibNXRFznyVTOlWVQv/r3uCMJz6Oi3ycbH8vmSvHUX7p3Fkr1jfwzMR5WU+iNzaZ+C0tXPNufgq3tqGRkX8eF8vS4eC1fj//bhXbXPdGXDDJ9je8xeDfjfGVY+r3q3n+swVs2NjEyvUb+fXzU7ji2S84+M/jOe3xiQnlU6WaSsforxbz8NhZ/Nt2ZXvfM7d1OeJPY3l6wtzY/+7r7H7q3XUUw/Wm63QKSIXLvTaoV2f26teNT+bGu5i6tm/ju74kHVWVUlSlUwjuf3dW+kJpSBMEGEm8gQh1dkfr7uRPfnQCL1y6T8J3D7v3/YRjyUbwjkL2bhl+7EMf8uR5e7LvNt2TyGf9Tudee/YTy431o627MWDzTinLOtw3pkVBNDY3J8wnpQqw+X7lBmYuWZfgDnQum3NdnUWtn81bwWE7bcF+d74LWPtdvTl1IVOSRHLueNNbCcem/bCG92YsYZ/+m1K7sYluHdr4LolIFxDT3GxYv9FKSQVWeHxzk0lor/uSz166nhtfnsZZ+/QDkt8Pv1DyQlISSkdE5gJb+Xx0vDHmJVe5sUBi6Bf8whhzb36kC47jP00VVZLp2guHqsoK2vsEIITBj4f0ic2rBKG6UrLadroQBJ33KDRTFqxiXV0jpyeLPHPhuNe8lsVrnqCJZCQbeaea7/p26Xr22rqb72fOt4K6bVLdgzV1DaxYt5HO7aoZ+od3We8zD+JgCGZhe9ubzNJq7xMBeuk/k++flIzz/v4pW3Ruy8LVdbx42b6+mzs693DOsvXMW76erTaNt+Bve/1r/vbhnNj/1oDVJCgrpw/xu/bJdEucpROoReFSSkPjUcA+np9xPuWm+JT7V4FkTInz0jgjSueGP3LGHrEy3lGt3+jVj+oKYbc+XXIX0oedtuwcuOyoKw8ITeH87thBodTjphDbVGfDMQ9+mKBwDr/vfS7/v88TJuD93GtA4Izg7k74SleG8soU29XWNzYlva/OM+sdWSczfD6evYJ+174eCxeu3djIcntt0TEPfMCBfxzLpHkr4xQO+EevzV2WGHbfIrNV3ptsNtk8SPs2lbHs5rmy0E49M+2HNb6Wjvs5HHb32Njfqzc08M3itXEKB1r6DO86NUd/+K3jcy/FcD9Eyea6CkUpKZ1lxpiJnh+/Xc/W+pRLPXtbIDrYuVx+vGdfoGXEtWnHxHDn0b84gJcu34/BW3Xj3H37pa27urKCzTPcLtvN1Ydun/Szmgzcdt5FqrmwScDN8DLBuxA3U5654EchSZKerxeu4fUpCxNG+M5+SNnsBQTx66tecgVbpLKy6xua4zq8C578lGk/WG4nZ5yUbi2Zw912tohP7OjFI+57P5Y+yQmr9p2j8dy7f0+an3KtU7IBRjIL4KRHJ3BxBgEUQWg2xteyfHz8bN/yRz/wAYfcMz7huOOa91p2jgL1DgiMMXHrtdyJfd0KqrHZ8NB7s3h3+uK0GyCGRSkpnZKnbXUl0393GL+2O3jnMfHzagzYvFPMcrn5mEG8dPl+KeuuqhTap8iplo6RAzdn7p1H+n7mnnxPF3LsKFaHTm2DeXBfu2Io227WMe7YkH7+7pxcyDUzr98aqnzzwaz4tUXOVg5BO3kvyRJFLkyRGLK+sZkbXIt/35m+JBb+LUncxcnCcZ1ABmdtmd/6nbHTEyPKauvjlch7PmXczEhitTQbw6dzV8RbAjazlya3nLLhtSkL49JBJWN9fSNL19YnfT6d125jYzPzV2yIHTe2/vDu3Os3Z+QMct2Dlac+msvdo2Zw/pOTUt7/MCklpXO0iNSKSL2ITBSR45KU211EVotIg4hMEZELCiplGtpWV8YCBa4auR2d2laxfc/0k6q7pnFxVVdWJCidwwb15JkLfsSB2/fgibOHsFPvTZJ+P9XGc+6Jxz7d2vuWee2KoVw1cruEBYt/PWfPlHI77NS7c0Ii0qDh0YWkXQ6KPSyW2+us1hVoZApWZ/emp/Nctq6exqaWTvCYBz8MtMDU4bnPFsTNr8xY1KIk/vvF9wnlve19N43SeWycvzXR0NTMyY9O4NfPTwksa7Z8MmdFoK03Bv12lJW5PQmOC9ObOuqZj+exurYhLmDkifdn+1rBjsXlHqy4UyYVKhCpVJTOq8AVwKHAGUAd8KKInOkpNx64EjgGOAmYCTwhIjckq1hELhKRSSIyaenSpcmK5YX9B/Rg6s2H0qltNVt3Tx0KLCJs2TV5J1xVIbTzTIReOXIAQwd058nz9uLggZvz2hX7071jG9/vJwsl/uXI7Th0UM/Y/yfs3jvu81uPHcS7vxzGTr0787ODBiR8v0+34IpjY5F9zUcGWGBb6AjBI3bumXDMWSCb6eZtueCdYwDL3XXvmPgO9c2pC1lVu5G1dQ1xSsSPVyf/ELdf0qH3JrqV3ITl/km2oWCx+esHidfYwckyscpjnd09aga73hq/n9Kj42b7zr+ttd2ybhfdSlc2FO9i43xR8Og1ETkYGB2g6DhjzIEAxpgrPHW8CEwE7gCecY4bY27y1PGyXfZ6EbnXGJOQIMsY8zjwOMCQIUOK1uu9+fP904ZRpor3r6gQ2nsWlbavDn57k3nNrvAokp+O2JZtNuvIZf/8nN5d2nG2HZ6ZDL/sCsnINIuwQ8eaqrSj/rd/cQCH3Ts+ZXr8w3fuyetTEyPAqiqEG48ayG9fmeYbiZRPBm6xCW9MjbcwnFFskJQ++aS+sZkH34sPY58wezmPJZmv8CNVBmQvC332DnJzyD3j2KRt+nnANRmcM2rc+tpXacs0NDX7vksvffE9vbu2Y7vNO/p8q7wtnY+AHQP8nJ2sAmNME/AcsKWIpBuePgu0BXbOWfI80ra6Mu0Lkyryan19Y4Lrx7s9toW/dknlXov7tghH7LwFT52/F6N+kbgTKhC3uLWmuiJwVF22WQ3c7r9dkrghqyokbX66ZBPzB++4Oefs24+5dx5Jx5oq7jhh5wSr6LS9+mQodTAKreRyxS+BaioyCYZIlyH6m8XrEnaaLTfmBchf19DU7BtK/vs3vuayf37Ounr/fiTZ/kxhU3ClY4ypNcZMD/CTLiFSbBuPkMpFnvoUm72JSEImg6CT+NDiXrt8+DaA1dEO3CL5HNAB2/WgY41//X87t2Uep01lBf84fy/27t+Nh12h4X4EDbX2BjM0NRvuO3U3du/bJWl+rwoRjtutt+9nx9suQ3cH+MbP9ucie5M6byTxaXv15cHTd+fb24+IHdupd/Cw8kzo1Tl681phkmzhpZI9jU2Gc/72SdLPkwU2lLOlkzMiUgWcDHwXIBz6dGADMDVNucjjTgEzdNuW1eG3HbcTBwzojojw6JmDGX/1cD6/cWRCJJmXru1bLCtH6Vx96A7MvfNInjhnCG/8fP/Y57ceO4hrDtshkJz7uWSrqqygc7tq/nXRPmmTkh62U+L8hR+OrIcO2hywfNTH7tabFy/bL2WCr2T1x9ZANJvYdR3YaxN2thWJ39YEIhI3D5Zqa4lXfzo0VXNS0ra6JF/RwHhT0pQi2S7ozhcbm5pjmzH64Z0XcvAGAeWLyD/RInKaiPxLRM4WkeEicirwHjAYuMZVbn8ReV1ELhCRg0TkBBF5GSuo4BZjTLixkEXgmQt/xNn7JCZmOHPvrWId42E79aTvpu3p1sE/YMDh/tN2523Xvj3pvGtn79OPSw/cJrCsJw/eMnDZ3xxuKbPrjtiRz28cmba8o9QuGGpZIskWIXbv2IatNrWi7QzJI+8cy6mp2fDEOUP49PqDgZYQ4CCux1TXO9uItyk3H0JVhh1BuijHsNgkAyu63KnJYN4yCiRbq1aovISRVzrAHGAz4G7gbeAxoB44zBjjzjSwEKs9twJvAE8BPYDTjTF/KKjEeaJ7x5o4CycbnP5zty27xM1xVOaaAtjDXSftwmyX+ykVw7bvYclQIWmVJVhK6p1fDktY1+PFa6F071gTUyhuKlyWTtvqyth1cZRZkCuTSpZsL22HNlWBtmIYc9UwThmyJR9dO4Jbjt0pdtwbaZgr7pD7PbbqGmrd6eiUxmovJv5zp9mz/4Dc3vF05LpAOlcir3TsjAIjjDGbG2OqjTGdjTEHG2NGecrNMsYcbozpbYypMcZ0NMbsa4x5tliy5wNn1J2rjqisjK8n7FGOiPgmLv30+oN5/pL41D6ZZrqtrqxgmx4dfbNxu8+4aYc2PH7WEH48pA9b2VZOj041/PPCH/HomS3zS706W5P1nT0ZEFosnfQyOcrST6Y2lRVs3b0Dw7brkb4iF5UVktTlcf0RO8b+3qRdFXedtCu9urSLV5AhD1z7JrEUC8GF+/cv2rnTkU1W+GQ8f8k+ec/87JeNvJBEXuko8aytt+Lq+20aLMW/F6cfciwb5/9sMltnQ49ONTlnGnAUZLqXfcuu7dm+Zyf+cNIuce3bb9vu7NCzZdR+ybBtuO/U3Tjasx+RE1cQxL3WprKCGbcd5hss0aNTDe/96kB+c0SwOTE37v1/dnAtIj7OZcW491/KV0ADwPer6vjw2hG8/rOhBU+JH+bc1l0n7sKIHTYLrb5Uz6Gf1ZJqTd6Qft1Sby8dAh99m3pDyHyjSqfEGDmwJz87aAC/OWIHBmzWkVOGBJ878cNxQQUNmc4HvX0WvaZyVzkKJNV8x659unDHCcmj5N3traqs4Njdeie445yX3y+QwE+mmqpK3w34nGPJooP27Nc1tgj098fvFPdZtSt07tEzW/Zgcrvd3GuhKisklqvPu8kbwK8PS55jLxWd21Xz84O2pXeXdgzq1Tlpx9ilfXj58tyyuq9rOreqH389Zwg97fDzXfp05sy9++YkmzuyM5V7zS8k/IHTdo/9fcIeLYMHbzqhA7fvwYOn7046BmRxPYqJKp0So2NNFVeN3I72baoYfdUw7jpp15zq81o+hWbunUcmrE/64JrhvHjZvgAM6pUYtu2V1S+44qajBqZclxOouTFLJ0BZG28H1M7VWTqfed14HWqqYhmevXM4bkuna/uW+a5K13Fv9FS1x3Xq5tJh23D3SbtwXQqr6+AdW6wAR4lM/u0hjNhh86TfAWse7+/nBkt7lI7KCuGyA7eN/e+2Jtp5FPsNR+6Il36bxrsCt+/ZiU3aWfNCgmS0aNmNE0Bx8zGDYso9VSCBs/bsV4dsFzvmdpn+8pAWxeoMhBylevWh27N73/RzZ+nWOh2za6+0dRQSVTqtjKG2ue+NqEqR1b7gbNm1PZ1sRfSfi/dhwm9GxH3unn+ae+eR3OqaPL/xqIHs0LNTyjVGEEzpjBy4OTv07JRR1J67c7zzhJ35+neHxf6vsi+yd7Fss2kZ3VZ5boRb6cT97boGXkvM6dT8rBER4eQhfehYk9wiuc41XzT+18OZevMhCWX8DJ1ThvRh1y27cMmwbehvu5DO329r33N8eO2IuMzme23djZEDN+feH+8GkJCHz+0e9S4i9rM0OnkGMt65sZoAW8KftlefuLVYQGxtWkNTcyyDSCpLx0nttFmnlkW+7ud3U1fgjHP4hiMH8vQFezGoV+e4z5Pht62Bm5MCRJIeNijYcoUwiFBXoxSCO07YmXd/OSw22nY6wGK611LRoaaKLTwLJFMFPQzeqitvXXlA2jDlIO3t2qENb115AP17BHdfbLd5JwZs1pH7T9udH+8Zn6WgZ+e2PH7W4ASXiTEmlviyymO1uF1y7s9SXQOnE0y1qZrXOnLPMzjX+yf7b80mbasTOnBIvidNRYVw7eE7xNaUuTtkkZaIui02aRs3L3nh0K35y9lDYm4078JjR9rjd++dMLL3u5fe61hZIS3RiBJ88r+yQmJuOSC21myLzm1j4ek7pEjY22CnrnKvmXPfO7fb0HGHtq2uZP8BPWJ/XzwsdRBFU5Nh6s2HsHtf/8wf3mvhx6NnDU5bJixU6bQyaqoq4zrRv567Jy9cum/BFoaFQRiuwHwp2bbVlYy+ahjH7NrLdy7okEE9EzpxY1osh6qKCn5z+A6xOR73vJV7fqc6hWnq3MtUCR6czu7A7Xsw544jYspt801qaNemkmm3HMpvDk90W7llToXzuXuOZ/zVw/nDSbvwxY0jqaiIX1zrKCdHGXb0rAMSV+CLN8WLnwL2Xp+qComlJBESLZ+5dx4Zs84cFq+xEoO6s8BfcuA2fHLdQfTv0ZFThvTh/V8PZ0i/eBeYe88lxypzb4nhdaHedNRAu40JzQDg2jSLshubDZ3aVidVpF7r2ctlGVjyYVA6PY2SFzq3q2Zwgddc5EoYrsACBesFotkY166ycPGwbXj4DGvk6e6gKpL87cXpuL0uKjeOJVHf0ByXXeGJs605mQ41VSnPkS7CyknAummHNrxw6T6888th9OnWnurKCrraLiN325y5C+c6JLN0AK4YsW3cZ84gpHO7av5y9hCrbs/oXpC4bRT8FNXSdfHZp52s7vef2mKZ1lRVxPLhiQh9urVn323iI9TcSshROn4BJg5OMEEypZMukMVRwn6BI+Df1rtO2iX2968DZhoJC1U6SlH4yf5bp9zfJxVhrCkKEpEWhDFXDeOFS/fN+Hvu1DjGtIRnJ5ufyQTnOzVVFdx/mn/0k2NJODt0OqP5oNkTFqfJ+Ox07+3bVDF4q25s4+OidIIhBmzWMebudTr0c5LslmuAE/bYkjl3tMy1OM/DQTtshnO5vM9IRUWLIrMiDROv68/tjOqzbz+CN3++f2xuq3P76pgl6Hc/Nt+kbdwGiO66nQ3r3Lv6evW1o9yztb7TZWf3W2DsDcYoJNFd5quUNdcfOTCj8h9dO4J973wXCMc1Fpalk034LsDOrnQ1ltXgn3InGwXrdIwbmwxdkmz57XQ6Tuby24/fmWN26xW4Pdv37OS746fD5cO35ZM5n7DjFsnnOxwLxd3G3fp0YfJNh9DZE3rtXBbjE8bufL/JmITjFwzdmp17d6ZT22oeOXMwL37xPf27d/AddFy4f//YItQdPYEoT1+wF//6dH6guSAR4YDtejD+m6WxQIJ21ZX07dae71bUJliJjlIIsvauV+e2/ODZ4dOxDpO9Fn7P0PJ1xdtTSC0dpSTo5dpFNAxLJ6zMDrkw/urhgJWx2+k4/Jp201EDeetKK/nqLccM4uU0W5dvYWdY2KJz25jF0aNTDU/YridosWicPZratalk+PbBF0z+6ZTdUn4+bDtrrihVEIbT2Xrvp1fhgEvp+NTjlO/avk2c0v729iO47ogdOdoOGe7VpR2XD982pnAc6+q/l6W3VH/Uf1Pu+fFugS3kJ8/dk1m/PzzmXqupqogFIXjXMrVvU8VjZw3mqQv2Slvvh9eOSDi219bxi61TRaKNv3o4AzbryEE7pg5/zydq6SglR5iBBMWM2uu7aXs+ue4gunes4bw5K5LKc/7QlrDjZG4nNyN22Iy/nD2E4dv34INZ1kZvO/TsxMEDWzoax9KpT7FHUyo61lTx0+HbJmzi5iZdB10ZwK107eE78OGsls3q/KaSDtyuB3eesDPH7tabT+euSFouUT7rd7bZPVJRUSFUIPx4zz48Pn427WsqufrQ7bn4gP50ad+Gsb86MC4H2qEBQ5a917RX57Y8ZGfBcD46YY/evDUtMfl+j0419N20PaOvshL93nXiLrQtwvbrqnSUkiOUlD12FcUOKHDmMGKRVWE0TYSRtoJx5mpOHuIJ396kLWfu3ZfT9sp+Zf6vDt0+pdJJR0zppLgJlwzbhkuGbcNLX3zv+/nOvTsjIpxqt8NRYEFSyRwwoAevTP4htNxpvz5se776YU3csWsP24GrRm4XW0DqBFH0S7M9fVB269uF9p5tNdxh6pcduA2Dem3CMxf8KMF1esqe+dl4MB2qdJRWSUVM6UQjjM0Y/zmdXNmic7u4SW6HigrhtuNy30z3H+fvlXLDsCBkux3NZzccnNDhtsz9pP/+3SfvwlUjt0u771RQ3BkUHCoqhLYV+bMm3O13otfcz5ATmTY0z5mrM0GVjlIyPHPBj3jjy4Wh1BWFnHNuMtm7x+HwnXom+PMLTaaZs920hImnb3P3jlZKoy1defo27ZiY5qhl7ie91qmpqgzN4igW7rD4VPNeUUKVjlIyDB3QPfQRW7Hdaw5OJpNM5HnkzMKtIs8HTRko2qEDuvOXs4dw4PaplZzjKku1k2up4hf6vJsrC8He/Tfl/ZnL6N3FctmeuEduyYDzRfndGUUJgPP+btm1eHvEuMkko3W54CjaoNGIIwemj7jao29Xrj50+4QURKXOZzccnJBVffQvDmDA5i0h6ZcO24Zjdu1Fn27tfV2qUUGVjtIqad+miofP2CMhhUmxcOYgomJ5ZcK5+/bLajLesXTC3EBQRLh8eOLcSqnj50p0Kxyw5o+SbckeJVTpKK0WZ91EFLj9hJ24e9SMOHdJqXDzMYOy+p6T6+yoXaJzH5T8IyZImEcrYciQIWbSpEnFFkNRWg11DU0p85IpiYyatogKV1h8FBCRz4wxQ9KXVEtHUZQiogonc4IuJI0qmgZHURRFKRiqdBRFUZSCoUpHURRFKRiqdBRFUZSCoUpHURRFKRiqdBRFUZSCoUpHURRFKRiqdBRFUZSCoRkJXIjIUmBell/vDixLW6q80Da3Dlpbm1tbeyH3Nm9ljAm0z4UqnZAQkUlB00CUC9rm1kFra3Nray8Uts3qXlMURVEKhiodRVEUpWCo0gmPx4stQBHQNrcOWlubW1t7oYBt1jkdRVEUpWCopaMoiqIUDFU6OSAifUTkeRFZLSJrROS/ItK32HJlioicJCIviMg8EdkgIjNE5A4R6eQp11VEnhCRZSKyXkTGiMjOPvW1FZG7RWShXd8EETmgcC3KHBF5S0SMiNzmOV5WbRaRI0RkvIiss5/ZSSIywvV5ubV3PxF5W0SW2O39XETO95QJ1BYRqRCR34jIXBGpE5HJInJi4VqTiIhsKSIP2DLX2s9wP59yobdRRH4iItNFpN7uMy4JJLQxRn+y+AHaAzOBL4HjgGOBqcC3QIdiy5dhWyYC/wHOAIYBVwKr7OMVdhkB3gcWAKcBhwHjsGL7t/TU90/7+z8BDgL+C2wAdit2W5O0/zRgIWCA21zHy6rNwMVAA3APMBI4FLgGOKpM27uLLcuhEJMAAAnFSURBVNN79vs5EnjMvs+XZtoW4PdAPfArYLhdVzNwRBHbeCCwGHgDGGW3rZ9PuVDbaNfTbJcfDtxm/39pWpmL/WCU6g/wc6AJ2NZ1bGugEbiq2PJl2JYePsfOth/gEfb/x9r/D3eV6QysAO53HdvVLnee61gVMAN4pdht9WlnF2CR3cl6lU7ZtBnoZ3cyV6YoUzbttWW6HdgIdPQcnwhMyKQtwGZ2Z3yLp653gClFbGOF6+8L/ZRO2G20v7sE+Ien3N+wBijVqWRW91r2HANMNMbMcg4YY+YAH2K9vCWDMWapz+FP7d+97d/HAD8YY95zfW818Crx7T0GazT9b1e5RuBfwKEiUhOi6GFwFzDNGPOsz2fl1ObzsUaij6YoU07tBWiDJecGz/FVtEwtBG3LoXZ9z3jqegbYWUS2Dlf0YBhjmgMUC7uN+wA9fMo9DWwKDE0ljCqd7BmE5VrzMg0YWGBZ8sEw+/fX9u9U7e0rIh1d5eYYY2p9yrUBtg1b0GwRkaFYFt1lSYqUU5uHAtOBU0XkWxFpFJFZInK5q0w5tRfgSfv3/SLSS0S6iIjjXrrH/ixoWwZhWQGzfMpBtN/5sNs4yP7tfVYCXQtVOtnTDVjpc3wF0LXAsoSKiPQGbgXGGGMm2YdTtRda2pyuXLew5MwFEanG8lf/0RgzI0mxcmpzL2AAcDdwJ3AIMBp4UER+bpcpp/ZijPkSa87jWOB7LJkfAi4xxvzLLha0Ld2AVcb2I6UoF0XCbqPz21tnoGtRlVJUJR1+i5yk4FKEiD2afRlrbuo890cEa2/QcsXmGqAd1kRoMsqpzRVAJ+BcY8x/7WPv2pFOvxGR+ymv9iIiA4AXsEbgl2C52Y4FHhWROmPMPymzNich7DY6/2e1yFOVTvasxF+jd8V/VBF5RKQt8ArQHxhmjFng+ngFydsLLW1eAfiFjXd1fV5U7LD267EmXms8cxA1ItIFWEsZtRlYjmXpjPYcfxsrSm0Lyqu9YAUSNGBF5zXYx94RkU2B+0TkWYK3ZQXQVUTEYwlErc1+hN1Gt0Wz0FWum+dzX9S9lj3TaPFtuhkIfFVgWXLGdje9AOyFFR451VMkVXu/M8asc5XbWkTa+5TbSKK/uBj0B9piTYSudP2AFSq6EtiZ8mrztCTHnVFrM+XVXrDu4WSXwnH4BGvCezOCt2UaUANs41MOov3Oh91G51nyPiuBroUqnex5BdhbRPo7B2xXxX72ZyWDiFRgxfEfBBxrjJnoU+wVoLeIDHN9bxPgaOLb+wpQDZzsKlcF/Bh42xhTH34LMuZ/WGsLvD9gKaLhWC9iObX5Rfv3oZ7jhwILjDGLKK/2ghUKv5uItPEc/xFQhzUiD9qWt7A66DM8dZ0JfGlHrkaVsNs4ASs02q/cCqwI3uQUK7681H+ADlgd01QsP/ExwGRgNp51AVH/AR7BXqMC7O352dIuUwF8BMwHTsXqrMbaD1kfT33/wrIWLsRSZM9jveR7FLutaa6Dd51O2bQZy6J5F8vNdglWIMHjdpvPLbf22jKeZLdvlP2OHgI8aB/7c6ZtwQrAqAOuwgpQeATLQjw6Au08yfUeX2r/PyxfbbSfoWa7zzgQK/CoGbg8rbzFfjBK+QfLT/oCsAZrDuAlfFYDR/0HmGs/rH4/N7vKdcNaALYCqMVaNLarT33tgD9jjTTrgI+BA4vdzgDXIU7plFubgU2worcWY41opwCnl2t7bTkPx1KcS+139H9YIfKVmbYFqARuwNpduN6+fidFoI3J3t2x+WwjVoaLb+xyM4HLgsirWaYVRVGUgqFzOoqiKErBUKWjKIqiFAxVOoqiKErBUKWjKIqiFAxVOoqiKErBUKWjKIqiFAxVOorig4ica2/9u639/5UickIR5ekiIjeLyB4+n40VkbFFEEtRMkYTfipKMK4EPsDa5rcYdAF+i7WV9Oeez5LtB6QokUOVjqIUCRGpMSHkKTPGRDnZpKLEoe41RUmDiMwFtgLOsF1uRkSedH2+q4i8IiIrRWSDiHwoIvt76nhSRBaIyD4i8pGIbMDaKhsROVVE3hWRpSKyTkS+EJFzXN/tBzjJFv/ikuFc+/ME95qIbC8iL4rIKlumiSJymKfMzXY9A0Tkdfvc80TkJjsJrFOuo4g8ICLfiUi9iCwWkTEiskOOl1ZphajSUZT0HI+Vs2oU1v7w+wC/A7DnWD7Cyln2E+BErKSaY0RksKeezliJF5/Fygn2f/bx/lgJGM8AjgNeBZ4QkUvszxcCznzSHS4ZXvcTVkR6YbkCdwV+CpwCrAJeF5HDfb7yIlYy0OOw8gfeApzj+vweu45bgJFYyR7/h+XyU5SMUPeaoqTBGPOFiNQDy0zitg93A98BI4wxGwFEZBTW/vE3YnXkDh2BM40xL3vqv93527YwxmJtqnYp8Kgxpl5EvrCLzPaRwctVWBtv7WOMmWXX+wbWPie/B970lP+TMebv9t9jRGQEcBrgHNsH+Kcx5q+u77yIomSBWjqKkiUi0g4YBjwHNItIlb1PiQBjgAM8X2kEXvOpZ4CIPCsi32PtdNmAlYJ++yxFOwCY6CgcAGNME5aFtZu9R44br8X0JfE7TX4KnCsi14nIEBGpzFIuRVGloyg50A0rFfyNtCgL5+enWFv/ut+xJXbnH0NEOmJtIb0rcC2wP7An1vYC7m20M5Vroc/xRVgKsavnuHd74XqsnVUdrgAeA87HUkBLROQen50oFSUt6l5TlOxZhbVx1UPAU34FjDHN7n99iuyDFaSwvzHmA+egbTFlywqgp8/xnrYMKfew92Ksbap/A/xGRLbC2iDsTqw9ea7JQU6lFaJKR1GCUY+1EVYMY8x6EXkfy0r53KNgguJYCw3OARHpirXTpff8eGVIwjjgShHpZ4yZa9dZibU98RfGmLVZyAmAMWYe8CcROQPYKdt6lNaLKh1FCcZXwP4ichSWm2qZ3aFfBYwHRonIX7HcWt2BPbB2p7w2Tb0fYe08+5CI/BZrG/QbsPag7+wqtxgrKu5UEZkCrAfmGGOW+9R5D3AuMNqucw3WAtLtgCMzbDciMgF4BWtr9nVY81i7Av/ItC5F0TkdRQnGb4AZwH+w5jVuBjDGfI41B7McuB94G7gP2BlLGaXEGLMUKyS7Eits+g7gCeAZT7lmrOCCrlhBCp8CRyep8wdgKDANa4/757HmeY40xrwVuMUtjMcKmf4nVtDBScAvjDH3ZVGX0srR7aoVRVGUgqGWjqIoilIwVOkoiqIoBUOVjqIoilIwVOkoiqIoBUOVjqIoilIwVOkoiqIoBUOVjqIoilIwVOkoiqIoBUOVjqIoilIw/h+OxPsfb5e+7AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x6192fce80>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "loss_names, x, ys = trainer.losses()\n",
    "plt.plot(x, ys[0])\n",
    "plt.ylabel(loss_names[0])\n",
    "plt.xlabel('Iterations')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sample\n",
    "n_samples = 10000\n",
    "z = hamiltonian_prior.sample(n_samples)\n",
    "xgen, dW = bg.flow(z, temperature=1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mdtraj\n",
    "\n",
    "torsions_train = mdtraj.compute_dihedrals(training_data_traj, mct.Z_indices)\n",
    "\n",
    "gen_data_traj = mdtraj.Trajectory(xgen[:, :66].detach().reshape(-1, 22, 3).numpy(), ala2_pdb)\n",
    "torsions_gen = mdtraj.compute_dihedrals(gen_data_traj, mct.Z_indices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "def periodic_convolution(x, kernel):\n",
    "    x_padded = np.concatenate([x, x, x])\n",
    "    y_padded = np.convolve(x_padded, kernel, mode='same')\n",
    "    return y_padded[x.size:-x.size]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "torsion_hists_train = []\n",
    "torsion_hists_gen = []\n",
    "xticks = None\n",
    "\n",
    "for i in range(torsions_train.shape[1]):\n",
    "    htrain, e = np.histogram(torsions_train[:, i], 50, range=(-np.pi, np.pi), density=True);\n",
    "    xticks = 0.5 * (e[1:] + e[:-1])\n",
    "    hgen, _ = np.histogram(torsions_gen[:, i], 50, range=(-np.pi, np.pi), density=True);\n",
    "    \n",
    "    htrain = periodic_convolution(htrain, [0.25, 0.5, 1.0, 0.5, 0.25])\n",
    "    hgen = periodic_convolution(hgen, [0.25, 0.5, 1.0, 0.5, 0.25])\n",
    "    \n",
    "    torsion_hists_train.append(htrain)\n",
    "    torsion_hists_gen.append(hgen)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "torsions_simple = [4, 5, 6, 7, 8]\n",
    "torsions_complex = [0, 1, 2, 10, 12]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(-2.84159,0.7,'RNVP + MCMC')"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2kAAADKCAYAAAA/6Ah7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xl4VOXZ+PHvmUkm+wYkISGEkAQISyDsi6hsghsqUuq+tGq1rXZ9W9/6s+6t1qqtb62tdaValyoqgiIgYRFlh7AkQICQEEhC9oXsyZzfH2dmciYzWQaSmYm5P9c1l2ebcx5wmDn3ee7nfhRVVRFCCCGEEEII4R0Mnm6AEEIIIYQQQog2EqQJIYQQQgghhBeRIE0IIYQQQgghvIgEaUIIIYQQQgjhRSRIE0IIIYQQQggvIkGaEEIIIYQQQngRCdKEEEIIIYQQwotIkCaEEEIIIYQQXkSCNCGEEEIIIYTwIj7uutCgQYPUhIQEd11O9BN79uwpVVU1sqfPK59X0Rvk8yr6it76rIJ8XkXPk8+r6Eu6+3l1W5CWkJDA7t273XU50U8oipLXG+eVz6voDfJ5FX1Fb31WQT6voufJ51X0Jd39vEq6oxBCCCGEEEJ4EQnShBBCCCGEEMKLSJAmhBBCCCGEEF5EgrQ+pLKyEoPBwFtvveXppgghhBBCCCF6iQRpfcju3btRVZXp06d7uilCCCGEEEKIXiJBWh+ya9cuQkNDSUlJ8XRThBBCCCGEEL1EgjQv19rayp/+9CeSkpJ4+OGHqa6uZvDgwTz++OOebpoQdpYuXUp0dDSff/65w74bbriBlJQUmpqaPNAyIRxVV1djMBj417/+BcCbb76JwWDg3LlzAMyfP58rr7zSk00Uwka+X0VfIZ/VniNBmpe7/fbb+fOf/8xvf/tbwsLCuPnmm1m8eDGPPfaYjE0TXuXBBx8kIiKC559/3m77hg0b+O9//8tLL72EyWTyUOuEsLdv3z5UVWXixIkA7N27lxEjRhAcHGzbP2nSJE82UQgb+X4VfYV8VnuOBGle7M033+T9999n1apVLF26lIqKChYvXswrr7xCREQEn3zyiaebKITNtGnTuOWWWzh06JBtW3NzM/fffz/Lli1jwYIFHmydEPb27t2Lj48PqamptnVrUHby5EkqKiokSBNeQ75fRV8hn9WeI0GaF3v55Ze58sormTlzJnv27AFg4sSJGI1G4uPjqa2t9XALhbA3ZswYSkpKKCsrA+CFF17g9OnTvPDCCwD88Y9/ZNSoURgMBj799FNPNlX0c3v37mXMmDH4+/tjNpvZv3+/LSjbu3cvgARpwqt09v3a0NDAddddx+jRo0lLS2PRokXk5OR4uMWiv+rqXmD+/PlMmDCBtLQ0Lr74YjIyMjzZXK8lQZqXKigoYPfu3Vx++eWAlnoTHBzMiBEjbPuTkpI82UQhHFiL2hw+fJjTp0/z5JNP8sgjjxAXFwdoX8xffPEFl1xyiSebKYRdz9nRo0epra21re/bt48BAwaQkJDgwRYKYa+r79cf//jHHD58mIyMDBYvXszdd9/tyeaKfqyrz+rHH3/M/v37ycjI4Fe/+hV33nmnB1vrvSRI81LWJ2CxsbGAdkMxYcIEDAYDu3fvpqSkhCVLlniyiUI4GDFiBD4+Phw+fJhf/vKXxMfH84tf/MK2f/r06fJwQXic2WwmOzvbFoS17zlbv349s2bN8lTzhHCqs+9Xf39/Fi1aZDt2xowZ0pMmPKare4GwsDDbcnV1tSea2Cf4eLoBwrmgoCBAGxsB2k2EtdLYI488wsiRI5k/f77H2ieEMyaTicTERP71r3+xe/du0tPT8fX19XSzhLCjqiqgZSSA9v2akJBAREQEX3zxBTt37mTFihWebKIQDlz5fv3b3/7Gtdde6+YWCqHpzmf1lltuYfPmzRgMBr744gsPtdS7SZDmpcaPH09ycjLPPvss0dHR5OTk0NrayrJly9iyZQvbtm2Tm1/hlUaPHs3KlSu58cYbmTt3rqebI4QDo9HI0qVLee2112htbWX79u0EBARwzz338O9//5v77ruP66+/3tPNFMJBd75fn376abKzs9mwYYObWydEm64+q//5z38AeO2113jwwQedluzv7/pGuuORI/DrX8O6dZ5uidsYjUY+++wzkpOT+cEPfoCqqixfvpyGhgZ27Nhhq0gmhLdJTk7G39/fofyuEN5k+fLlPPPMMxw8eJDMzExbIaZNmzbxj3/8w8OtE8K5rr5fn3vuOVasWMGaNWsIDAx0c+uEaNPde4G77rqL9evX24qMiDbe35O2ZQssXgzV1fC3v0FWFiQne7pVbjF69Gi2bt3Kn/70J5544gmqqqowGo2ebpYQncrPz2fixIm28ZRCeCM/Pz9+85vfcP3115OcnMwrr7xiK9QkhLfq7Pv1hRde4L333uOrr74iPDzcA60Tok1Hn9WKigoaGhqIiYkBYMWKFURFRTFgwABPNNOreXdP2urVsGiRFqABNDfDs896tk0esHv3btLS0iRAE33Cnj17mDx5stN9Tz31FHFxcWzbto27776buLg4ioqK3NxCIdpIuX3Rl3T0/Xr69Gl+/etfU1lZydy5c0lLS2PKlCkeaKEQmo4+qxUVFVxzzTWkpqYyYcIEXnnlFVavXo2iKB5opXfz3p60VatgyRJobbXfvnw5PPYY9KOn9Hv27GHx4sWeboYQXaqqqiInJ6fDG96HH36Yhx9+2M2tEqJje/fuZciQIURFRXm6KUJ0qrPv17i4OFtBHCE8rbPPamJiIrt27fJAq/oe7w3SfvObtgBt+HAIDYX9+6GpCV54AZ57zrPtcyMpoyv6irCwMMxms6ebIUS3Pf300zz99NOeboYQXZLvV9FXyGe1Z3hnuuPZs3D0qLbs5wdbt8KTT7bt/+c/obzcM20TQgghhBBCiF7knUHat9+2LU+dihoTQ9mMGahjx2rbamvhpZc80zYhhBBCCCGE6EVeH6Q1Tp7MSy+9xEsvv8zK0aPbjnnxRWho8EDjhBBCCCGEEKL3eGeQ9s03tsUtZjPlltTGA6NHU2ktK1tebt/jJoQQQgghhBDfAd4XpDU0wJ49ttUMf3/bsmo0ckw/R9rGje5smRBCCCGEEEL0Ou8L0vbs0So4ApXR0dQFBdntzk1IaFvZtMl97RJCCCGEEEIIN/C+IE2X6njSMhu5nj5IU3fsgLo6d7RKCCGEEEIIIdzC+4I03Tiz/KFDHXbXBQdTHBkJgNLcLOPShBBCCCGEEN8p3hWkqap9kBYf7/QwfW9a49q1vd0qIYQQQgghhHAb7wrSjh2DkhIA6v39KR040LYrJiYGf0sREX2Q1rxunVubKIQQQgghhBC9ybuCtPapjoa25l122WVMmTIFgDxdkBaYmalNbi2EEEIIIYQQ3wE+nm6AHV3REP14tICAABISEjCZTGzdupW6oCDORkURXVyMobWVpo0bMV19tSdaLIQQQggh+qHHH3+8y2PCwsL4xS9+4YbWuMZsNrNlyxaGDx/OsGHDPN0c4YR3BWk7dtgW9UHa4MGDURSF2NhYQkJCqKmpITchgejiYgCqVq4kUoI0IYQQQgjhJnfddZfd+gcffEB0dDRz5syxbTMajW5uVfeYzWY2b96MwWCQIM1LeU+QZjbD8eO21bODB9uWYyyl+BVFYdSoUezevZvchASm79wJgGHLFve2VQghhBBC9GtxcXF260ajkcDAQIftF6qlpQUfH++5ZRfu4T3/xwsLob4egLqAABoCAmy7BusCtpSUFHbv3m03Li3i+HFaKivxCQ93W3OFEEIIIYTojvz8fL799lvOnDlDfX09YWFhjB07losvvtguAHvttdcwmUxMnTqVzZs3U1payqJFi5g6dSrnzp1jzZo1HDt2DKPRyJgxY0hKSuLDDz/krrvusgsODx06xLZt2yguLsbHx4ekpCQWLlxIaGgoLS0t/OEPfwBg48aNbNy4EYD58+cze/Zs9/7FiA55T5Cm60UrHzDAbleMblLrhIQE/Pz8qIe2cWlmM4WffsqQO+90U2OFEEIIIYTonsrKSmJjY5k4cSImk4mzZ8+yZcsWqqqquO666+yOLS4uZv369VxyySWEhYURFBSEqqq89957lJeXs3DhQsLDw8nMzGSdkyrn27dvZ+3atUyaNIk5c+bQ0NDApk2bWL58Offddx8+Pj7ceeedvPXWW0yePJm0tDRAGz8nvIfXB2m+vr4M0K0bjUZGjBjBoUOHODVsmG1cWu2XX4IEaUIIIYQQwsukpqballVVJT4+Hl9fX1avXs0VV1yBn5+fbX9dXR133HEHkZGRtm1Hjx6loKCAG2+8kVGjRgGQnJzM22+/TVVVle24hoYG0tPTmTJlCldddZVte2xsLH//+9/Zv38/U6ZMYciQIQCEhob2eHqm6BleH6QNHjwYg8F+poCUlBQOHTpEXnw8U3ftAsB/925UVUVRFPe0VwghhBBCiG6or69ny5YtHDlyhOrqasxms21fRUWF3dCegQMH2gVoAKdPn8ZoNDJy5Ei77WPGjCEnJ8e2npeXR3NzM6mpqXbXiIiIICIigry8PNuUVsK79Ykgrb3k5GSMRiOndNVoYvLyKDp1ihipUCOEEEIIIbzIJ598Qn5+PnPmzGHw4MH4+vqSl5fHunXraGlpsTs2JCTE4f3nzp0jKCjIoTMiKCjIbr3WMnfwm2++6bQdERERF/LHEG7k9UGafjyalZ+fHwkJCZxobaUiPJyIykp8W1oo/vJLYu691y3NFUIIIYQQoiuNjY0cO3aMyy67jOnTp9u2nzlzptvnCA4Opra21iFrzBqUWQUGBgJw/fXXM3DgQIfz6NMqhXczdH2IG6iqSz1pAElJSQCcio+3bWtOT++lBgohhBBCCOG65uZmALvhO6qqsn///m6fIy4ujtbWVrKzs+22Z2Vl2a0PGzYMHx8fKioqiI2NdXhZAzfr/G3Wtgnv4x09acXFcO4cAA1+ftRbngIYDAaioqKcvsUWpA0bxoQDBwAIOXCApqYmTCaTGxothBBCCCFE54KDgxk8eDDffPMNgYGB+Pv7s3fvXurq6rp9jpEjRxIbG8vKlSuZN2+erbpjaWkpgK13LSAggAULFrBu3TpqampITk7GZDJRU1PDyZMnGTFiBGPGjEFRFAYNGsTRo0cZPnw4/v7+hIaGEhwc3Ct/B8J13tGT1r4XzfJBi4qK6nCm9sjISEJCQsjT9aQNzcsjVzd4UgghhBBCCE9btmwZUVFRrF69ms8++4zw8HAuu+yybr9fURRuuukmhg8fzrp161ixYgUAl1xyCWCfxjh9+nSWLVtGSUkJH3/8Me+++y6bN29GURS7zo+rrroKg8HAu+++y6uvvkpGRkYP/WlFT/COnjQXUx1B+7AmJSWRUV1NbWAgQXV1BNbXc3bTJkampPRqc4UQQgghhND7xS9+0eG+AQMGcNtttzlsf/TRR+3W77777g7PERwczLJly+y2rVy5Ej8/P7vpqkCrhJ7Sxf1wQkIC9913X6fHCM/xup60ii6KhuglJSWBotiNS2vZtKnHmyeEEEIIIYQn7d27l507d5KTk0N2djaff/45GRkZzJgxw2G6KtH39dmeNIDExERAKx4y+sgRAAZmZVFZWUl4eHgvNFQIIYQQQgj38/X1ZceOHVRWVtLa2kpERAQLFy5kxowZnm6a6AVeEaSpx49jLSaqD9LaT+TXXmBgILGxsXbzpcWfOsWJEyeYPHlybzRVCCGEEEIIt0tNTSU1NdXTzRBu4vm+UVWFY8dsq9Ygzc/PD39//y7fnpiYSOHgwTT5+gIQXlVF4fbt3b682Wxm69atvPvuu2zcuJH6+noX/wBCCCGEEEII0XM8H6SVl6NUVQHQ5OvLOUvpz/DwcIdZ1Z1JTk5GNRrtxqWxeTOtra1dvldVVT788EM2bNjAsWPH2LJlC6+++ipFRUXn92cRQgghhBBCiAvk+SDtxAnbor78fkRERLfeHhcXh8lkIi8hwbZtyLFjHNP1znVkw4YNHLGMZbOqqKjg9ddfd2mCQSGEEEIIIYToKZ4P0jooGhIWFtattxuNRkaMGEGuLkhLyM1l3759nb4vIyODb775xum+lpYWPv30U7a7kDYphBBCCCGEED3Bq4I0ffn97vakAaSlpVEQG2sblxZRWUnxrl3U1NQ4Pf7UqVOsWrWqy/Nu2rSJlpaWbrdDCCGEEEIIIS6U56s7dtCT5koJ/cTERIIjIjgVH0+yJX0y/uRJ9u/fz+zZs+2ObW1t5bPPPsNsNnd53sbGRk6fPk2CrpdOCCGE8FZVVVVkZGRQV1eHwWDAYDAQFhZGamoqAQEBnm6eEEKIbvLaIM2VnjSDwUBaWhp5CQm2IC0hN5et+/Zx0UUX2RUg2bt3L2VlZQ7nGDduHDU1NeTl5dltP3nypARpQgghvJqqquzfv581a9bQ1NTksP+bb77h9ttvZ+DAgR5onRBCCFd1O91RUZSRvdEA9eRJ23KFLjBzdTLqtLQ0h3Fp5eXlnDp1yratsbGRTZs2Obx3yJAhXHvttU7nnsjJyXGpHUIIIYQ71dfXs2LFClauXOk0QAOorq5m+fLllJeXu7l1QgghzocrY9KOKIqyQVGUZYqi9EwPXEMDiqXcvVlRqA4JAbRJqk0mk0unioiIwDRrlt24tLDKSrsCIt988w11dXUO773yyivx8fEhMTHRYd+ZM2dobGx0qS1CCCGEOxw9epR//OMfZGZmdnlsTU0Ny5cvp6Kiwg0tE0IIcSFcCdJ+CAQAHwCnFUX5o6Iowy/o6vn5tsXq0FBUoxFwvRfNasKUKXbzpQ3LzeXAgQPk5ORQXV3Ntm3bHN6TmppKbGwsoAV67a+tqiq5ubnn1R4hhBCiN9TV1fHxxx/z/vvvd1gky5nq6mrefvttefgohBBerttBmqqqb6mqOgtIA1YAPwGOKYrypaIo1yqK4nqlSF3wU6kLjs43SBs9ejSnk5Nt6wm5ubYJq99++22HSo1Go5F58+bZbRs+3DHulJRHIYQQ3qKmpobXXnuNgwcPOuwLqa7m+++/z70rV3Kj2UyYv7/DMRUVFezcudMdTRVCCHGeXA6sVFU9oKrqT4FY4F4gGvgYOKUoymOKokR3+2S6Ih1VunnRzjdI8/X1JXTxYtt6giUIbGhooLS01OH4adOmOVzLWcrjSd24OSGEEMJTmpub+eCDD5ymLPo1NHDrf/7D6CNHGLxvH6OeeIKfPfcci77+Gr+GBrtj9+zZ060qx0IIITzjQuZJSwDGW/7bBBwCfgUcVxRlSbfOoA/SdMGSK5Ud25tw1100+/lp56msZFBxsdPj/P39ufjiix22O+tJKykpcSmdRAghhOhpqqqyatUqzpw547DP0NLCnZ99RtTZs/bbS0uZsWED9//tb4zPyABLYFZVVUV2drZb2i2EEMJ1LgVpiqKYFEW5RVGULcBBYDHwDDBUVdXLgWHAl8AL3TqhLkir7IGeNACjvz/KFVfY1iccOOBwjI+PD0uXLnU6Z0xQUBDR0Y6dgdKbJoQQwpO2bt3qNMXR38+Pn+7bx+CsrLaNP/whDB1qWw2urWXJp59yx/Ll+NfXA7Br165eb7MQQojz40oJ/ueBAmA5UANcAySpqvonVVVLAVRVrQBeRAvWutZBT9qFBGkAPj/4gW059cABFF1KR1RUFPfccw/JurFr7TnrTZMgTQghhKfk5uaSnp7usN1oNHKPqjLg88/bNj71FLz+OuTkwDvvgKU4FkBCXh4L1q8HtPHWJSUlvd52IYQQrnOlJ+024HUgWVXVq1RV/VxVVdXJcUeAHzjZ7kDthZ40AK64AgYNAiCsuto2Nm3q1Kncc889REVFdfp2Z+PScnJycP7HFd9J9fXak+iHHvJ0S4QQ/VxzczOfffaZ033XzpvHgL/+tW3DPfe0fW/5+MAtt8DRo6i//rXtkIn79tnSIqU3TQghvJMrQVqcqqoPqqqa29lBqqqWqqq6vMuztbTA6dO2VWvhkJCQEHx8LnAaNl9fuPlm2+p1VVX8/Oc/t82H1pVhw4ZhMNj/1VRXV1NbW3th7RJ9R1YWvPkmvPQSVFV5ujWir2pp0Xoyyso83RLRh6WnpzstFHLRRReR+uWXYB17PXQovPgiKIr9gcHBKM89R+X06QAYVJWFa9eCqrJ//34pxy+EEF7IlSCtXlGUac52KIoyWVGUVpeufOYMSqv2lnNBQbRaJqG+4F40q9tvty2Grl9PuAuBn8lkIjIy0mF7cQdFSMR3kLVsdU0NvPaaZ9si+q6tW+G22yA62u7BkRDddfr0abZv3+6wfejQocxLSYHnn2/b+NRT4GSstVXA3/+O2RLAJeXkkHz8OE1NTWTpx7IJIYTwCq4EaUon+4yAa7mA+lTHHhyPZjNpEowdqy3X1cHHHzse09wM585BeTlUVtrtcpYSKUFaP6IvHvPii9pnRQhXrVyp/be1FXrqu030Gy0tLay0foZ0fHx8uPbaazE88YT2+wYwYYKW2tgJv8mTyV+0yLa+cO1alNZWjh071qPtFkIIceG6DNIURTEoimK0Hm9Z17+CgCsAx4nIOtPbQZqi2PWm8b//C1ddBUuWwPTpMHgwmEwQEgIDB0JEBNx0E1jGnUmQ1s8NHAjW3tT8fPjoI8+2R/Q9qtoWpAFce63n2iL6pIyMDKdzfM6ZM4eBZ89qxUGs/vxnMBodjm3P8OSTNJpMAESWljL+4EFOnDhBa6tryTBCCCF6V6dBmqIojwLNaPOgqcA3lnX9qxp4BPjQpSv38ETWTt1yC1jHlhUWwhdfwKefws6d0G4uGQDef992My5BWj+nKDTcdVfb+vPP2wJ4Ibrl0CGwVoUNCYE5czzaHNG3mM1mtm3b5rA9NjaWmTNnwgsv2OY8Y+FCuOyybp03duJEdlx6qW09bd8+mpqayNP9JgshhPC8rnrSNgFPAE+ipTu+YVnXv/4fsAT4uUtX7oWJrB0MGQL6G+32DAYICoLAwLZtv/0tNDQ4DdJKSkqkwmM/UVVVxctAs3Us4549sGWLR9sk+hh9L9oVV4Cfn+faIvqc7OxsysvLHbZfffXVGGpq4L332jY+8ki3z2s0Gqlatsw2Ni0hL4/wigpJeRRCCC/TaTUNVVU3A5sBFEVRgddUVT3TI1furfL77b3yCvz851r1q9paaGjQxhsNG6bNHePjo41HS07WKrDl5sJf/0rYgw9iMploamqynaqpqYmqqqqeb6PwOrW1tdT4+3NgwgQm79mjbXz+edA9gRaiU/qS6ZLqKFzkrBctKSmJmJgY+Pvf28aipabCrFkunTt+2jROJCUx4vhxAMbv309mcjKLdOPVhBBCeFa3C4eoqvp4jwVo2M+Rpu9JCw0N7alLaBRFKyAydy5cfTV873tw8cUQH68FaKAN6H/iibb3/PGPKGfPSsqjYNuMGW0rq1drRWaE6EpBAVjmn1J9fODKKz3cINGXnD59mlOnTjlsnzlzppZ2/c9/tm28917HkvtdSE5OZn9amm19wv79lJWWOu25E0II4RldjUl7Q1GU4brlzl6vd3YuO2az0560kJAQjN0Y+NwrfvQjGDNGW66pgd//XsrwC8oiIzkTG6utqCrs3u3ZBom+QdeLdjI+nlc++IAvv/zSgw0SfYmzXrSoqCgSExNh2zZtvCNoqfq33ury+YOCgqiZN48GSwrugIoKhp46JSmPQgjhRbrqSZsLWLu25lnWO3t1T3EximXyzHp/f5osc1KF6dIe3c7HRxuIbfXmm8RaKmDpSZDW/xRYgzSw9Y4I0SndeLSjo0ZRVFTEjh07PNgg0VdUVFRw+PBhh+0zZ85EURQthd/qppvgPH83h48eTea4cbb1tP37OW5JfxRCCOF5nQZpqqoOV1V1v2U5wbLe0Sux21ftoPy+R4M0gEWLtPL8AK2tDLM+rdSRIK1/MOkC9IIhQ2zLqgRpois1NZCebls9OmqUBxsj+prdu3c7FKgKCQkhNTVVS7f+4IO2Hffdd97XGTlyJBkTJtjWx2Rmkp+dbTcOWwghhOe4Mpl1z+mg/H6Pj0c7H0uW2BYjnFTzKy0tlflk+oGQkBDb8hldT1qzkzQkIeysXQuWG93CwYPtxtwK0ZmWlhb27dvnsH3atGnaUIAXXwRLFgqTJsGUKed9rZiYGCpSUigbMAAA/8ZGkrKyOGmdNkIIIYRHdTtIUxRllqIoV+vWByqK8p6iKAcVRXlON+F11zooGuLxnjSwq8Lms2EDob6+drtbW1tlcHU/4Ofnx7BhwwAojYykyfI5MBUXoxYUeLJpXVJVlaysLNLT0/n22285dOgQZ86ckekj3OXbb22Lx0eM8GBDRF+TmZlJfX293Taj0cikSZPg6FF45pm2HQ88cEHXUhSFESNHcnD8eNu2kdnZ5ObmXtB5heh33npLK95jfZlMkJQEDz2kVRTX27RJO8bHB7KzHc8VFwd33qkt/+Uv2rE7d3Z87YsvhuHDtTHzubmO7Rg5En75S6io6Jk/a3fMmaNdv6Oqs3feqe2Pi3PcV1sLTz+tPYQKCQF/fxg1Cu6/H/Tp2O64hhdwpSftGWCybv3PwJVANvBj4KFun6mD8vteEaSlpGgfaoC6OsaXlDgcIimP/cPcudowS9VgoDAmxra9dtMmD7Woe9auXcuHH37I119/zfr161mxYgWvvfYaL730EhXu/KLur6xTNgBndKmy/paxt0J0ZLeTwkRjx44lMCAAfvITWw8t06fD7bdf8PWGDx9Otu5BQtKJE5ySIE2I8/Phh1phn88/14bPPP00/OY3zo9tbe16fsNbbtGCubffdr7/5En45hvtu0Bf4fV3v9PasX69Fqy88oqWJebOB7UhIVob2gc9dXXw8cfa/vYKC2HaNHj2WbjqKvjoI1izBn72M+1cy5a5/xoe5kqQNhrYDaAoii/wPeCXqqouRZvQ+uZun8mbe9LArjdtpJMB3BKk9Q/Dhg0j1pLqqC8eUrd5s6ea1KWKiooOC1SUl5ezevVqN7eonzGbQZeuVqAL7uPj4z3RItFHFBYWcvr0aYftU6dOhf/k2PxgAAAgAElEQVT8p22co9GoleA3XPhohWHDhlEYE0NtYCAAwbW1qBkZMi5NiPORlgYzZsBll8HLL8OCBfD669rvQnsLF8J//wv793d8vqgouPxyeP99aG523P/vf2uBV/sHNomJWjsuvVTrzXvwQdi82e63qVsSEuCxx1x7j9X48dr8w++8Y7/944+1Njubk/G227QgaudOePJJ7Zi5c+GnP9Uqaz/6qPuv4WGufMsHA9WW5WlAEGC949sLdPsORPXmnjSwC9IG797t8A+sxEnvmvhust5Y64uHKLqeEm+zp4u25eTkUF1d3ekx4gIcO6YVDgHOBQVRoxtna02fFcIZZ71ogwcPZkhgIPzqV20bf/5z7WawB4SFhREWEcGJpCTbtqRjx5wGi0IIF02aBPX1UFrquO/++yEmBh5+uPNz3HGH9v41axz3vfMOXHSRllrZmalTtf+6O5XvttscA6h//xuuvx6Cguy379wJGzZoQaWzYQKKAtdd55lreJArQdoZwFoK6grgkKqq1i6lCKCuW2ex5s1aWKs7+vr6EhAQ4EJzetGMGWCZI823rIy4M/ZzeJ89e9YTrRIeEGfJZ9YXDwk5csS9aQPd1FHRgfaysrLc0Jp+ShckF8bE2KWgSJAmOtLQ0MDBgwcdtk+ZMgXl0UfB+mAwLg4ef7xHrx0fH8/x5GTbetLx404n0hZCuCg3V5siY+BAx30BAVqAtno1bN/e8TkWL4aICMeUx2+/1YKuO+7ouh3WYkDuLmJ1222Qk9M2TrugQAuSnKVqf/WV9t9rrvG+a3iQK0Hae8AfFUX5CPgVoA9dJwHdmwWzogLF8qS5ydeXekuaRVhYmDYHjDcwGrV/GBajjhyx211eXk5LS4u7WyU8wBqkVQwYQL1lTJF/bS3NXjjpa1ZWFnV1XT8ryczMdENr+il9kKYL7H19fYnRpT4KoXfgwAGa26Uz+fn5Md5ohH/8o23jiy9CcHCPXjs+Pt6uJ21ofj6F7X7zhBDd0NoKLS1akY433oAVK+Cpp7R7SmfuvltLTXyok5IOfn5w442wahVUVrZt//e/tYIX3/++43vMZq0ddXXauLSnntJ67S6+uOPrqKr2Hv1Lfy7ry5Xq5sOHw+zZWltB6/GKjdXSC9vLz9f+6+rDTHdcw4NcCdIeA/4E+KEVEdHN/MwE4MNunaV9L5olMPOaVEcrXcrjaCc35FKAoX8IDQ3VyvErit24tCrrExkvssvJHG6DBw922Hb69Gmqqqrc0aT+RxektR+PZuiBMUTiuykjI8Nh24Tx4/H97W/bbormzbObIqanDBs2jLrgYNvn1Wg24/v11z1+HSG+81JSwNcXBgyAu+6Ce+/V0ho74uurjfnauLGtl8eZO+7Qpt7473+1devytdc6n8z+3nu1cwcFaWPfkpPhyy+13ruObN6svUf/ysvTxm3pt82f362/Cpvbb9fa2tio9QbeemuPjKd1+zU8pNt/ClVVW1VV/YOqqotVVX1CVdVW3b7rVFX9S7dO5CTVEbxkjjS9yy6zfaAHFhcT3G4cT1lZmSdaJdxMURRbb5o+SGv65htPNcmpjooOXHnllURaUnf1JOWxF5jNsHevbVXfkyZFQ0RHioqKKCwsdNg+s6xMewoO2g3HX/9qX8GthwwaNIiAgABO6FIehx892uPXEeI775NPYNcu+OILrWjIyy+39fB05JZbYOxY+H//r+Njpk/XAkDruVat0nrrOkp1fPhhrR0HDmi9b19/rRXZ6Mzkydp79K+YGLjnHvttr7zS+Xna+/73tWkInngCDh3quCrt0KHaf3U1K7zqGh7i/lCzgyDN63rSAgLaBluCw7g0mSut/7AFabriIT5Onnx7UkdFB+Li4hgzZozDPgnSesHx47aiIbWBgVRL0RDRDc7GkcZFRRH+xBNtG+69F1JTe+X6iqI4jEtL9sJ0biG83rhx2gTzV1yhjTUbOVIrwV9b2/F7DAatt2rnTli5suPjbr9dK7efk6P1Fg0erPWSOTNsmNaO1FTnPW3OhIRo79G/TCYtdVC/bdSo7p3PKjRUGwP2zDPa+0ePdn7cggXaf1etcu387rqGh7gymbVJUZRHFUU5oihKnaIore1e3Ruk1VeCNNAKiFjEteulkJ60/sNZ8ZDwEydcy83uRS0tLR0XHVAUxo4d67BPUh57gT7VMTbW1uthNBoZogvwhbBqaWnhwIEDDtsvO368rRJbeLj2hLgXxcfHczoujgY/PwDCpAKsEBfGzw/+/GcoLtZ61DqzZInWKfD73zsv1w9agQxrj/qaNVoPXEdj3bzN/fdrdR5++9uOj5k2TUul/OMfO65C2VkQ645reICPC8f+GfgpsAb4GGg8ryv2pSBt+nTb4pB2QZr0pPUfMTExGAwGakJDORcURHBtLabGRs7t30/wpEmebh4FBQVOiw6kWp68R0ZGEhkZ6TB1RFZWFjNnznRbO7/z2ld2tIiLi8PHx5WvWtFfHDlyhIaGBrttAa2tDNVXcnv0URg0qFfbER8fj9loJCcxkTFO5gYVQpyHa67Rgq/nntOCiM7GhP3hDx33jIFW2XXePHjpJa3IR3eqOnqL2bO1V1feflvr7Zo6FR54QHuPyQRHjmiFWJqb7epFuP0aHuBKuuP3gEctY9IeUlX18favbp1FH6RFRNiWvTJI0/WkxRYUoOh6TqQnrf/w9fUlOjoaFIUiXSGOCi+Z1DrfWrFIJzExEZPJZFt31pt2WG7GelYHlR2HWvPghWjHWcGQq06cQLE+UImPhx//uNfbERMTg6+vr924NCFED3jqKa037Z//7Py4yy6DOXM6P+aOO7QALS2t19KfPSomBnbs0FJEV63SehgXLYK//AVmztSqZfaFa/QgVx7vBgPbLuhqqoqam4t16LNXFw4BLRd36FDIz8fU3ExUcTFnLU/Ia2pqaG5uxtfX18ONFO4QFxdHYWEhZ6OjST5xAoBmJ9UUPcHZnEbtC1WMGTOGTZs22W0rKCigpaVFenl6QruiIfrKjlJ6XzhTVVXFCct3iZVffT0pq1e3bXjkES1tqpcZjUbi4uI4OnIkLUuWaCX5n3uu168rxHfCnXdqL2cWLrSfV3XOnI7nWd24sfPr3Hqr9upIQkLPzuGq61RxSbt7Dafeesv59uBgbUqCzqYlcNc1vIArPWmrgEsu6Gq6OdIaTSbqLV2/wcHB3nuj2Mm4NEl57D+s49LO6nrSfLyg+Iaqqk570tr33gwaNIjgdvMrtba2UlRU1Kvt6zdOnADLOJ66gACqdZkBzqZBEMJZwZB5+/ZhtI4HGzHCrSlN8fHx1IaEcGDCBGp7eC42IYQQrnMlSPsbcJOiKI8oijJFUZTE9q8uz9BX5kjTk+IhgrYgrSg62rYtLDeXVg8XDyktLaW+vt5um6+vr0NgoJ9KQM9Z2X5vkp+fzxtvvMGLL77Ihg0bHMbeeQ19L5quaIjJZCJCl9YtBIDZbHZIdQysrWWyfn6yxx8HNz68lGkihBDCu7gSpG0DRqBNar0DOObk1bm+VDTESl88RMrw91sREREEBARQNmgQLZaKSmFVVVS0S1dyN2e9aEOGDMHopOqTswqDZ9p9pr1Jbm4uy5cvJz8/n8rKSrZu3cobb7xBZWWlp5vmSHfDXaRLbxw8eDBKL8xtJfq2nJwch+qqs779FqP1gUtqKtxwg1vbFBcXR2RkJBMnTuRaLxo4L4QQ/ZUrj+l+CFxYsqsuSKvy9vFoVpMmaU8zW1qILC3Fv76eBkuapvSk9R+KojB48GBOnjxJSWQkMZY0weqtWxk0cqTH2tWd8WhWznrSvDVIKyws5L333nPoqSwqKuJf//oXS5cuJSkpyUOtc0I3BcJZXW9rtG5ZCKv2qY7+dXVM0xWe4fHHtXLbbmQymfjJT37i1msKIYToWLd/BVRVfUtV1eWdvbo8SV/sSQsI0CrpWOh706QnrX+x3nDrKzw2O5lE2p2cBWkdVROMjY116NWpqKigtrOJNj2gtLSUd955h6amJqf76+vreffddyksLHRzyzqhm+tKH6TJeDTRXm1tLUeOHLHbNn3HDnytpfjHjfOqEtBCCCE8w+VHdYqiGBRFGacoyqWKogS59OaTJ22LfSZIgw7HpUlPWv9iDdL0N+GeLB5y7tw5Kioq7LYpimIfpFVXa708e/ZgMhiIjIx0OI839aa1trby3//+l7q6OgB8GxtZ+uGH3PfyywzTPeQxm8189dVXHmplO5WVYEk7bTUYKBs40LZLgjTR3oEDBzDrJqw1NTQwY+fOtgMeesjtvWhCCCG8j0u/BIqi/BQoAg4A6cAoy/ZPFUX5WZcn6Is9adDhpNa1tbU0Np7fnN6i73HWkxbswTFpznrRoqKi8DOZtLlZBg2CsDAYPx6mTIG773Y6Ls2bioccO3bMNum2saWFGz74gHGZmUQXF3PD++8TpgtKnY3r8QhdqmNJZCRmy3hARVGcBsWi/1JVlb26IjMAU3fvxt86Fm3ECPj+9z3QMiGEEN6m20Gaoij3AC8CnwLfB/R5U18DS7s6h9pBkOb11c/a96Tp5qGQlMf+IzIyEkVR7MrwDywspN4yrYS7OSsaEh8fD888A7//PbTv6V2+nCQn1Si9qScty9IzqZjNXL9iBUk5ObZ9AQ0NLPvwQwwtLbZtB3Rphh7TwXi0QYMGuW8eRbNZm5hTP8eW8DqnT5+mtLTUtu7T1MTMb79tO+B3vwMnRX+EEEL0P670pP0KeF5V1R8Bn7TbdwRLr1qHWlpQzp0D7OdI8/Pzw9/f34VmeEBSElhSmAIaGhigC8wk5bH/8PHxYdCgQTQEBFBlKXbj09pK+bYLm+P9fDnrSRt79Kj9BI3+/qArzDP8yy8d3nPmzBnUnpwA8zw1Nzdz9OhRAK5avZoxhw87HDOkoICF69bZ1jMyMjzfdl2gWOyp8WgZGXDNNVrBCeG12hcMmbJ7N0GW1F7i4zufqFYIIUS/4kqQNhxY28G+WiC8g30aXREA/Rxp4eHh3l+iWlFg2jTbqhQP6b+cpTzWb9/u9nY0Nzc7TEQ9uKCA+IcfbtswZw5UVcF//2vbFPDBB7QfSNrY2OgVDxuOHz9OU1MTQ/PymKxLCTP/7Gfw/PO29ek7d5JiCeDKy8s9n66p70mLirItuzVIswaqu3dzxsPFbIRzTU1NZGZm2tZNjY3M3rq17YAHHwR39bwKIYTweq4EaaVAQgf7RgGd50zpxm71qVRHKwnSBNqYL7BPa1P373d7O4qKiuyKDxhaW7lhxQoU61P55GT46CMwmeCyy7R1QKmqYpYu7djK44EObamO43VBT8GMGRj+8hf45S9hyRLb9rnp6bbApP2kwG6lqh2mO7ozSGvy87Mtl/773267rui+rKwsu4qlM7Zvb+tFS0iAu+/2TMOEEEJ4JVeCtFXAI4qiJOq2qYqiDAJ+iTZWrWO6Hyf9HGnh4Z13wHkNXZAWqwvSvKEHQriPrcKj7gbc35Ki507ty8/HnT5NuPWzGBqqjU+yVhk0GEA3/1Hqli124yrB80GaNdVRaW1ltK5iZusDD2jtVxR4/XVaLWnSUSUlDLWke2ZmZtLc3OyRdpOXB5YxiXUBAZwLCbHtcuccaQ26IC10yxa3XVd0n/5hQkBdnf1YtMce0x6oCCGEEBauBGkPA43AIeArtImt/w84DLQCT3T67r7ekzZ1qm0xprDQVrxAgrT+xdo7UqS7AY/Iy7Pr1XKH9kFa0vHjbSs33AApKfZvuPNObc4/ICQnh6Htio54unjI8ePHaW5uZvjJk7behZrQUGKXLWs7KCICbrrJtjrFktbX2NhoG8vmdu3nR7OkboeEhBAU5NoMJRdCH6QNycqi0RuqXgqb8vJy8vLybOsXffMN/tbfxJQUGYsmhBDCgSuTWZcBU4CnAV/gBOADvATMVFW187uC9mPSLPpMkBYZiTp8OKAVi4guLga0iXWtczqJ776QkBD8/f2pGDCAJsv4keBz56h0UuSiNxUUFNitJ+uDtEWLHN8QEQG33GJbnaqflwkoLi6mRVc10d2sqY7jDh2ybSu+9FKM7cboGH/6U9vymKwsAiwTcR/Svc+tdKmOxZ4ajwa0+PhQNmAAAKbmZspWrHDr9UXn9L1owTU1TNuxo23nk09KRUchRO8pKoKdO+HQIW2+Yrln7TNcmidNVdUaVVWfVFV1tqqqI1VVnamq6uOqqlZ3+eYOetL6TLojoHQwLk160/oPRVGIjo5GNRjsiofUbNrktjY0Nzfb5hIDCDx3jlhrz5rRCPPnO3+jLsAZfeQIvrp/k2azmbNnz/ZKe7tiTXU0tLSQcuSIbbvv7bc7HjxpEg0TJgDaw5I0y83vyZMnaXUyvUCva9+TZuHOVEer7JEjbcvqqlVuv75wzmw2s183bnXGtm34Wh+ITJwI11/voZYJIb7zVqzQxrxOnw6pqZCYqD20XbJEKyomAZtX6zRIUxQl3YXXhs6vZLD1PPTVIK2jcWn6eW/Ed5/1BrwgNta2rXXXLrddv6ioyK7svH4uMWbMgI7+TaWlaV/SgE9LCyOPHbPb3b53zl1yc3Npbm4m6cQJAhoaAKiMiGDIddc5Pd70wAO25cl79oDZTFNTk2fG1el70rwoSAvbutVh3KHwjJycHKqrteeYpoYG7TNr9eij2phLIYToae++qw1/0D2QBbTMtk8/1fYlJWlTuAiv1NWvgwFt0mrrKwWYg1blMcDy3zlo1R07raPfMmIETz/0EH/+n/+xzZEWHBzsvslee0IHPWkSpPUv1hvwQl2Q5qe7We9t7YMpu/Fol1/e+Zu/9z3bor5Ah7PzusuJEycAJ6mOPj5OjzfcdBNNgYEADCwvZ/jJk4A2rs2tGhpQs7Ntq8WRkbZld6c7ApyKj7eNTQsuLcXsxs+k6NgeXVA2ad++trFoI0fC4sUeapUQfVhtLRQUwLlz8jCqI2++qY11tWaYREXB6NEwZIj9cUVFcNtt4KniW6JTnQZpqqrOUVV1rqqqc4EXgWZghqqqiZZUx0RgpmX7i52dq7W1FRSFuuBg2+D6PjMezWriRFTL2IHIkhJMlqf+ku7YvzjrSQs/ftxtPxZ2RUPMZpIsQQ7gfDyani5IG3HsGL66saKeCtJycnLwaW5mlK74h6IrEOIgMJDKa66xrU6yzKnm9iDt8GEUyw9geUQEzZYAycfHhwGW8WHuYjAYMPv4cCIpybatVjc/nvCM6upqW1EbQ2sr0/VzKv761/27F83NxZbEd0BTE/zv/8KAAVqwERKizS3o5u9br/fii/DDH7bdk4wbp6XmZ2XB6dNw9KhWUdbysJNDh+Avf/FYc0XHXPmFeBL4vaqqdhUHVFXdATwGPNXZm52NF+lTqY4AQUE0jxoFaN2G1nFA0pPWv0RaekzKBg6k0VI2O6i6mnp9sNSL9MHU4LNnCbYUz2DQIJg8ufM3jxmjPU1DKzChLzhSUlJiN4+TO1RXV1NSUkJiTg5+lmuXDRhA7JVXdvq+gPvvty2POHYMo2Vy73PnzvVqe+10Mh7N4Oabb2tGgj7lUZFxaR63d+9eW2rymKwswi1VN9XISO3pdX9VWQmzZoGHxsGKPig7W/vM/OlPdoXoaG2FigrPtcubqCo89BD84hdt29LSYONG0KfgjxyppVo//njbtsceAydzqArPcuVOYgRQ0sG+YiC5szc7qxzX53rSAMOMGbZla8pjeXm5Z4oWCI8wmUwMGjTIoXhIVXp6r1+7qanJ7qGAXVXHhQu792ReV9Z+TGambVlVVYqKinqknd2VYxlPl6wbH5c/eTJBwcGdvi9k1iwqLdUU/ZqaSLSc54SbAmUAdBX69J8DT4xHM1keFhwfMQKzJVMh6OBBsFShFe5nNpvZa+nlRVXt5kVT7r/fNiVGv5STo/37mTVLu/kWojPp6TBpEujHc4aHg7+/59rkbWpr4e674emn27bNmqX93Q0a5Pw9P/+5bZw69fXwwAOSPuplXAnSTgL3drDvXiC3szd/J3rSAJ+ZM23L1uIhqqpSXl7uqSYJD4i1pDrqUx6b9JPT9pKzZ8/aFw3pqvS+M7qUx5HZ2fjoctHdnfKYk5MDqmoXbDYvWND1GxWF8jlzbKujLVMguDVI27bNtnh66FDbsieDtLqgIPItbVFUFT77zO1tEZrs7GxqLBOdJ5w8acu8UP397SaX75es32E5OdqNpO7fkrcym82UlpaSmZnJxo0b+eqrrzh48CBVMidh7you1ubHtGaM+PrC889DWZkWWDQ29u+HUaoK772nzbf4xhtt26++GtavpyEggEOHDrF161bS09NZu3YtX3/9tVbN2dcXXnml7T2rV2tj2YTXcD4y37nHgf8oinII+Ag4C0QD30MrKHJLJ+/9zvSk6YuHxOmqyZWWltrS4MR3X0xMDAcOHLAL0nx1ZbZ7iz6IMjU02E9KvXBh904ybpyW7pCdjam5maTjxzlqSYFsP0l2b1JVlZycHAaWlRFRWQlAo8nEgG4WUzDdcINWQhgYdfQoSmsrJ06cQFVVFKXTOkYXrrbWlu6oAmd0g7E9UTTEGqQBHElJYdipUwC0fPQRPnff7fb2CNhtmWwdYPbWrbZl5c47O36y3V8kJ8OZM9pNdlkZzJunPWRKStJuNq+7Drzk97S+vp5t27axZ8+eDudEDQ8PZ8yYMUydOrVPPnz2WqoK99zTFoRFR8OaNdrUFVYmE42hoZ5pnyepKqxdC0884fiQ4447OPPYY+xev57MzEyanRQFSU9PZ9iwYUydOpUx99yD8uqr2o4f/UgrMnL11W74Q3SfqqrU1dVRU1ODv78/YWFhvf877wW6HaSpqvq+oiilaMHa79AmtG4GdgGLVFXttAS/s560PhmkjRlDq8mEsamJ0JoaAs+doy44WMal9TPWnjR9hccwa/GQXvzi0AdRiSdPYrQOvk9Lg+4GB4qi9ab98Y+ANlbGGqS5syetuLiY2tpaxul60XITE0lM7jRz2mbw1VdTHRpKaHU1gfX1DMvLIzcxkcLCQtv/n16za5et8EFJZCSNurQbT/SkKYpCVFQUxcXFHBk9mkXr1gFgSE+H6mrojzcxHlReXm7r1Y0pKLBNk6EaDCj/8z+ebJpXMIeE0PDll/h973soJSXQ0AArV7Yd8MADcPPNWjqWZV5Ed2tsbGTr1q3s3Lmzy7G6lZWVfPvtt2zbto1Ro0Yxe/ZshrSvoidc9/rr9tkAb71lC9CamprIzs4mMzOTY+2mk/nOW7lSC86s6dRWUVGce/hhVkZEcHz58i5Pk5eXR15eHkPHjOH2sWPxyczUxvgtWwbr18Ps2b30B+haXV0dOTk55OTkcOrUKSorK+3iiNDQUBITE0lKSiIlJQWfDqpB93Uu/alUVf0K+EpRFAMwCChVVbVbJZrM7So5GQwGQkJCXLm8d/DxoW7ECEIsY3liioo4kZwsFR77mZiYGBRFoWzAABpNJvyamrSAPTubQEtxmd6gD6L047i44grXTqQL0kZmZ2NobcVsNFJaWkpjYyN+lkqFvcl6E6v/c1TMmNHtaTl8TCYKZswg1BKQjD58mNzERI4fP977QZquSp8+1TE8PNwtf3fODB06lOLiYiojIiiKjmbw2bMYmpu1J8833OCRNvVX+rL7s7/+um2HdV6ifq6oqIg/bdxIxE03ceN77xFV0m64e2Ojlnb15ptawPbCC+DGm7CzZ8/ywQcfUOFiQQpVVTly5AhHjhxh6tSpzJ8/32PfB33epk32BTDuv5/GuXM5dugQhw8fJjs722mG1ndabS389KfQPgDz9cX8wAN8O38+m/bupdXF4Tf5VVX8/eqrua+yEr8zZ7SHJosXw9dfa5k3blRaWsrGjRs5fPiw3dCO9qqrq8nIyCAjI4OgoCCmTZvGlClTCLRWrPyOOK8SZKqqmlVVLe5ugOZMeHi42yug9RTz+PG25cFS4bFf8vX11dJbDQa7lMeqDZ3P6X4h7IqGtBvH5XKQlpYGluDCv7GRoZb0OHBfymNOTg4+TU0k6CpK+biYYmG+9lrbcsqRI2A2u2dcmi69JD8uzrbsiVRHq6G6YPGIpWcU0CYtFW7T0tLCvn37ABhYWmobLwmg/O//eqpZXqliwAD++eMf86977uGjpUtJnzuXkoQE+4P+9jdYskSbE8sNDh06xOuvv+5ygNberl27ePnll/tfL8+F2roVFiyAuXNt49CaEhN5f+JEnn32WVasWEFWVpZbA7SmpiYa208I7W5ZWdpwG32AFhAAP/sZ5zIyeG3UKDbs2nXeRewqAwJ45frrabCm61ZWwlVXafOouUF1dTWrVq3i5ZdfJisrq9MArb3a2lo2btzIX//6V9avX099fX0vtvTCqKqqjQfsJo/1D/blvG3fadO0gZrAYMsHuLS01D1jYYTXiI2Npbi4mMKYGIZbAo2mbdt6rSjAmTNnbF9ckcXFhFVXazvCwkBX0KZbFEXLOf/HPwCtNy1v+HBA661LaH+j1MNaWlrIy8sjIS8PH8uPSsmgQQyZNcul80Rdfz21v/kNQXV1hNbUMOTMGU77+PRub6Cq2hcN0QVpnkh1tNIHaYdTUpizaRMA6uefozQ2gjzRd4usrCzbTcJF33yD9RehddEijLoHfEKjGgwUDhlCoSU98OtLLiEuP59Lvv6aEdYAZ/VquPRS7YGD7nPek1paWvjqq6/Yoava2p7BYGDo0KFER0ejKAr5+fkUFhZ2eENZXV3Nu+++y4wZM1iwYAFGyzyrogP/939aiqtOs78/by5YQJF+/LWblZaW8swzzxAWFkZ0dDRRUVGMGjWKIUOGuOeeLz0drrmmrXgKaFN4PPccZ1WVd999l2rr/YATfn5+jB07lpCQEBRFITMzk5L2vddAxcCBvHXDDfzwrbcwNTbCqVPadTdtaptTrYcVFRWxbds2DvHicp0AACAASURBVB065JBx56rm5ma+/fZb9u7dy0UXXcSMGTO8Ig3SWjk7KyuLrKwslwoNSpB2HgJ0N5Ixll6HxsZGzp071zdTOMV5iY2NJSMjw64nzaSbO6un5et+pEa0L71/Pl9EV11lC9JGHDvGekt1SHeMS8vPz6elpcUu1TFv9Ggmu9gTNTA6msxx4xi3U5u+cczhw5wZOpTc3FxG9Vba6cmTYPmBq/f3p1RXBMKTQVpERARBQUHU1tZSHB1NeUQEAyoqUGpqtB95V3tbxXmxFgwJqapivK6YkPHhhz3VpL5FUTgdH8+7N93EvPR0LrYWXdm7FzUlBeX3v4df/rJHHzqUlpayYsWKDqcgMRqNTJo0iYsuuoiwsDC7fbW1tezdu5ddu3bZqnm2t337dk6dOsXSpUvdPtF9n7Fxo/b/1UI1GjmYlsbGiy6i0kv+zqqqqqiqqiI7O5utW7cSGRnJxIkTmTBhQu+l2a1bB9deq6UggtZ79ve/ww9+wPHjx/noo4867OULDQ1lzpw5jBs3zm4YwSWXXEJubi6ff/65w1CdszExfPi973HTu+9iUFVt/PWtt8KHH0IPPmSora3l888/57Au06Arfn5+BAUFUVlZ2WlA19DQwIYNG8jMzOTmm2/22H15TU0Ne/bs4eDBg+ddAd5jQVqfLBpioYwfj9lgwGA2M7C8HFNjI01+fpSVlUmQ1o84K8Mf3ovFQ/RB2gWNR7OaNw+zvz+GhgYiS0uJKC+nYsAAzlimluhNtvnRdMFm/Zw5Lj+VVBSFc4sWgSVISzl8mPWXXcaJEyd6L0jT9aKdGTLEbm46T6Y7KorCsGHDyMrKAkXhaEoKM61t/eQTCdLcoKioyPbvdPbWrbbCPo1Tp+LnwUH43sZoNOLv709raystLS3Oe6IMBtIXLKAyPJyrPv8cg6qi1NXB736H+tprKHPmwPDhWqXaRYvOqziOqqrs3buXtWvXOq2ABxAZGckNN9zAwIEDne4PCgri4osvZtasWezcuZONGzc6PVdBQQGvvvoqt956q0eKilRVVVFcXExZWRllZWW0tLRgNBrx9fUlODiYhIQEYmNjPZMNlJ+vjde0/Hs5M2QIK5YupaKL4CwgIIDRo0czbtw4HnvsMTc01F5JSQnr1q0jPT2dcePGMW3aNGJiYnruAmvWaKm+1iBsyBBYs4aW0aPZsHYt23Vjo/UMBgMzZszg0ksvtav8a6UoCsOHD+dHP/oRa9asISMjw27/8REjWHPFFVz1xRfahk8+0f6N/ec/9pNin6cjR46watWqDqulWgUHBzNmzBgSExOJj48nwDK3ZGNjI3l5eRw8eJDMzMwOe7KLiop4/fXXufnmm4myzKva21RVJT8/n507d3L48OEL7h2UIO18+PtTExdHmGUcT3RREfnDhlFaWtrraWLCe0RHR2MwGKgYMIB6f38CGhoIqKujdscOgnSTnvcEVVU5bZnywdTYSLxuDBmXX35+Jw0I0EpfW76IR2Rns3PGDCorK6mpqenVBw4nT54koqyMgZanS02+vgSfZxARumQJDc8+i39jIwMqKoguKrIFgb1CXzREl+poMpk8niFgC9LQUh5tQdqnn8LLL7u1+EJ/pO9Fm6SrvGZ6/HFPNckrRUdH8+CDDwJaUbHMzEy2bt1KsZP5rvZOmUJJVBRXrV5NtGW/cuIE6MeeBgZqFeluu00r4R8d3eVnvaqqilWrVnU4hjW8vJwrs7NJamjAsG0bJCTA2LHaTbOTXjyj0cjMmTNJSUlh9erVTr+DGhoaePvtt7n99tt7vbiRqqqcOnWK7OxssrOzuzVuPjAwkKSkJEaNGkVycrJ7ip40NmJeuhSDJTvhXFAQH9xwAzUdBN3BwcGkpKSQkpJCQkKCV6SQtrS02IpYJCQkMG/ePLv0c5c0NGjf16++qmVAWA0dChs3UhwSwopXX3X6bwW0hwY33ngjcbrfpo6YTCauvfZahg8fzqpVq+zG+e2eNo0B5eXMtP7ebdigjWV/7z3QzVHqiqamJr788kvbmN2O+Pn5MXv2bKZPn+60kJifnx8jR45k5MiRzJ8/n+3bt7Nnzx6n4xSrqqp44403WLZsGUm9WLSptbWVw4cPs23bNgoKCgioq2PUyZMk5eQQXVRESWQkW2fPptzF6Vck3fE81Y8aZQvSYnRBmug/fHx8iI6OprCwkJzERMZabo5rP/64x4O00tJSGizpDsNzcmxP6NW0NJQLeHJnuOYaW5A20hKkAZw6dYqxY8deYKuda2hooKCggFmWvy+Ak8OHk5CScl7nSxg1iuyRIxl/8CCgpTxujImhsrKyd75nOpnE2tNjUvUPiU4PHUpNcDAh585p6ZlffXX+Ab3oUmNjIwcs6c4Xb91qG2t5LjWVYPl775DBYCA1NZVx48aRk5Njq9ynTx3Mj4/nlXvvZeru3cxNT8e/fXpXXZ1WUMFaVEFRtClJhg3TgquRI7UALjmZiooK9uzZw65du5yW1g+vqOCSr79mwv79GJwVYRg1SksTnzvX6Z8nIiKCW2+9le3bt/PVV185PElvbGzk7bff5rbbbuuVQK2roLczdXV1HDx4kIMHD2I0GklKSmLcuHGkpKR0XHXXbIbcXMjM1Oa+Ky3VXkajNh9gZKT29z9zpjZ5sk7DgQO03nwzQZZq2WZF4aNly5wGaIGBgSxatIjU1FS3f892Na4prLKSSXv2EFVcTKvRSPkLL3Bu6FAGPfAAkZZhBN2yYQPcfju0H3IwbBhqejq7y8pY9/77HRZNGThwILfccovLnSDjx48nOjraoaLpuoULaTaZuHjLFm1sbVGR9mD3Jz+BP/xBGw/fTWfPnuWjjz7q9D45MDCQqVOnMm3aNAL9/LTPUXGx9iorg4oKqKrS/g1efTUYjYSHh3P55Zcza9YsNm/ezL59+xx61hobG3nnnXeYMmUKCxYs6NGHD9aKrhtXr2bEunVcnJ9P9NmztrlfreLOnGHC/v1kpKWx+dJLu31+6Uk7T2pamjaPBFLhsT+LiYmhsLCQE0lJtiDNR//0q4d0NB5NudAUtiuvtC0m5ObaUnfz8/N7LUjLy8tDVVXGWn6YAXKnTmXUeX4nBAYGcvaii+D/t3feYXGdV/7/3GEYeu9V9A4SIEAF1GXZsixb1bJl2bIdl2TTdjdld5PNetMc/zabTZziOHZiW3JX7FjdstWsDqgBkuhNCNE7AkS7vz/uMBoQSCAGGOD9PM88MHfuvfO+M2dm7nnPOd+jddLCr1zhyJIlFBUVERcXZ5Ax62hrA706o2t6aUsTWY/Wh4uLCxYWFrS3tyOrVFyKjr4VTXvvPeGkjSEXL16kq6sL26YmYvWjaL/85Zj2TpwqSJJEYGAggYGBusyBEydOkJeXByg1SmlJSWTMnInP1as4NDbi0NBAcHExzgNryWQZKiqUmzYSIP/sZ1ydPZuvwsPpMDPDpbeXTo2GGq1KL729JKWmsuzwYdRDpD4CkJurXKg+9hhs2qT0khqQlidVVDB31y5mXrlC5s2blDg4UOHuTrOdHUiSLqK2adMmZsyYYZDXr6amhuzsbDIyMu65/kWfnp4eXRTO1NSUoKAgrK2tsb55E4/CQtzy87HOzESVldVf0GIobGyQly3jRmgotV1dNF69SsRnn2Gu91p/ed99lA6SjRQTE8OKFSsmTF7d1cqK/9i6lRqViorqajIzM6m/fBnva9eIycggNDdXqd/SJzMT9u6lOiCAzmefxeO738VkqPF3d8NLLyltcfTPo1LBqlW0vvwye9PTycnJGXKM/v7+bNiwQZcSOFLc3Nx47rnn+PTTTynou8ZQqTiyZAmlvr6s/fRTrNralPH98Y/w97/Dz36mRNUCA/ul/esjyzLnzp3jwIEDQzqXlpaWLF68mJlBQZgePAjPPgt798IQNZ6A4vj/+MfK51CtxtbWloceeoi4uDjef//9QVMpz549S35+PqtXryYgIGCkL9FtlJWVcfDgQcqLinjq7bfx0WY7DYVKlom7cIHYy5f5l2E+x7g5aRqNBktLS9ra2tBoNPdsSMaCmV6kRF/hUTDGdHfD3/6mrMpFREz0aPDy8uL8+fMU6jVgtsvMVC7mDfiDonPSZNkw9Wh9+PjQHhqKRW4uJr29BBQVkRMe3s8pNDRFRUU41tXhof3cdJuY0PPgg6M6p3rVKjr/+lc0XV241NbiXF1NYWGh4Z20c+cUG0RRo7yp9z3m6+tr2Oe6B/rq0vp+zDNjYm45aZ9+qkQArK0ncIRTE1mWSdPWRSYfP66LojWEheEwStuejkiShI+PD4899hgVFRUcOXJEJ2V/09ycgpAQ3b5fyDKe5eXMvnSJgMpKLJqa0Awiny/JMjPS03kyPb3f9iZbW65EROBeWalT6dWxeDG8+KLSED4vD/7851sXjh98oFN5JioK5s9XbpmZ8Ic/QEcHlsAc7Q2gxdqaci8vSv38OBcXx/bt21m7di0R9/hb1tbWxoULF7h48aJy/SHL2DU2Ym5mRsddfn/sGhuZdeECUZcuoenspNrVlWo3N8o9PSkKDKRD77utq7OT9r17iTx7lrCcHF0mx4hoaUH6xz+wBgZ+A/WoVBxaupQzetdV1tbWzJw5k1mzZuE8whQxg5OTg6mfH54mJnh6eBDf23t7tGsIXIuK4Ec/4sbLL1OxZg323/kOzl5e0NWlOPxHjsDu3bpFRgBcXOh88UXKly/nbGUlOZ98MmRtk0qlYunSpcydO3fUEUYLCwsee+wxDh06xKlTp3Tbi4KC+POLL/Lwzp0E9aUHV1XB888r/1tZKU3G586FefNg9mzw8qKtvZ1du3aRm5s75HOGhITw0KpVWH/4ISxdqnzWhkNenhJ1/Ld/U3q6rVoFixfj5eXFs88+y3vvvUdTVRVu1dW0m5sr9Y2SRFNTE9u3bycxMZFly5YNuy9rHz09PVy5coXU1FRd/f6q/ftvc9B6VCoqPDwoCwvDJjmZsKNHUX/1FQDSIBH8oRg3J83Z2Znvf//7dHZ20tLSMuFpQaPFJiVF979rdTWq7m6ampro6OjA3Nx8Akc2xcnMhBdeUP4fo0jPSOhLV2m2s6PG2RmX2lpMurqQjx5F0otSjZa+ejT3ykqd9H6PjQ0mI5XeHwTpoYeUHwuUlMec8HAqKiro7OwctOh4tBQXFxN56ZLufkFQEL7R0aM6p39kJAVBQURolaIisrNJ8/Wlt7fXsP0Y+5Tm6J/qCBhNPaq+k1bp7q6zS9raYOdO2Lx5gkc49SgoKKC+vh7bxsZ+tWjST38qomijxMPDg8cff5zi4mK+/PLL2/s4ShLXvb3ZpVeDo+ruxqalBfumJhwaGoi4fLm/Iq4eds3Nt+pu+pg1C377W0X2X5/vfhf+5V/go4/6b790Sbm9/vod52LT2kpYbi5hubkknTnDvgcfZMeOHcyfP5+4uLhhKT92dXVRVFREbloaTUeO4FBTQ2xdHW5VVXhev45FRwcyygLNkcWLadJmKKg7O4m5eZOg69fxyM3F7sIFJL2ojW1Li+4CvFeSKPPxodHeHrumJhzr67G9Q1TjhqUl1a6uNDg40GZlRZulJVJvL5ZtbVi3tjKjtBT7pqZBj612ceHTtWup0qbtm5qasnTpUhISEoyvl25PDwwRLbkaFsa58HC6TUww7e7Gv6iIyMuXdQs2Vq2tBG3fDtu33/EpKiIi+GzdOqpNTPrXpA2Cs7Mz69atM6hglUqlYvny5bi5ubFr1y5d37VWW1vee+IJIq5cYcXnn/e3hxs3lN9Gvd/HHktLWhwc8PPx4eqCBbQPWDRQq9Xcf//9xPn6Im3Zovw2DcTBQUlbdnUFJyflvloNH36opD2C4iy//rpyU6kgNBTHWbP4RlkZnDmDiXZRtcHenoLgYC7OnMl1b2/S0tIoKChg7ty5hIeHY2VlNejrIcsynZ2d1NXVkZmZSVZWVr8oXey5c8SfO6e7f2L+fLKiozGNjiYpJYXEiIhbdZOHD8OPfgT29vD558N6P8Y93VGj0QypkjSZMHVxocnJCbu6Okx6e3GtqaHSw4PKykqjuVibkuit7hAZqeTBTyAuLi6o1Wq6u7spDApSLoaB9p07sTSQk9be3q6L0uo3xpUffNAgQhDma9fCr38NKOIhUm8vskpFeXk5/treaYaitbWVmpoa1um9b5ejorh/lM/j4+PD7uhonZMWfuUKxxYupKKiwrBKakeO6P4t1UtTcnR0NBpl137fP5JEZkwMS/t+7N97TzhpY0Bfb62FX32lizJUBwfjun79RA5rSuHv789zzz3H2bNn+fLLL4dUYwToVatpcnCgycGBUj8/LsbG4lpVRdKZM3hUVCBLEr0qFU719VjoN75VqeDf/x1+8hMYbIHK01O5QPzmNxVxh+PHlej6YLVr8fHw9NOKwMmFC8jnzintMLTYNzXx+Pvvkx0WxqXLl/mrnx/2wcH4+vri4eGBra0tbW1ttLW10dLSQktLC42NjZRdvUpMWhr3f/45miFeAwmYmZlJ5OXLVIeE4NTcjOb69X5O2Z1QyTIzrl5lhr5AlR7XvLwonTGDMl9frnl5ccPa+s6LEbKMU10dAYWF2La0YNbRgaazkyp3d9ITEujWRjMCAgJYtWqV8ZXDWFgoGQj6vcUsLZX3eO5c2LoV79BQ6jMzOXr0KE1NTWTMmsUXK1YQe+ECCWlpt/qaDkGPSsVXixZxIjkZeRjOaXx8PCtWrBhxJGi4xMTE4OzszI4dO2jsq6+SJK5ERlIQGEhSWhq+paV4VldjOYgDb9LWhltbG27l5UReusSuhx+mIDgYyxs3iGpqYoG1NVb//d9w8iTop+cGBCiS/+vWQXT04Hb1q18pTe5/9zulVq2P3l7IzobsbAbKyTg0NpKQns7s9HQ+WbeOy9HR1NfXs3fvXvbt24e7uzt2dnZYWVnR0dFBQ0MDjY2NSunAYJ8bWSY0J4eVfQqYQGZ0NEXPPcfKRYvw9fW9PRi1ZIlyDdvcrDhqw0BIfY2ClqAg7LQG4l5RQaW2Pkk4aWOInmgD8+bBxx9P3FhQFL28vb0pKSmhIDCQOdoVWenAAYM9xzW9lbswPSdNvWGDYZ4gMZEOOzvMm5qwvnEDr2vXuObrS1lZmcGdtOLiYpxranQqbV1qNfXz5w+5ijVcTExM6Fy+nO5PPkHd04N7VRUO9fUUFhYazkm7ebPfSmGx3ufcUHUlhsDV1RVzc3Od0ExWdPQtJ+2LL5Qi7HGSI54O1NbWUlhYiFNNDbP0pKxv/uhHIopmYCRJIiEhgaCgIHbu3Elpaemwj612c2P3ww8DinBZfHw8DlFRkJoKn3yiXID/8IcwHNGn5GTlBtDaqrQAOXlSuQAzMYFnnlFUIPXef6mnh7bz58n4zW+YuXMnllrnMDwnh3Bt5LvaxYUSPz9y/fyQJYng/HyCCwtR9faSHxxMS3g4j1y8qKt/HowOMzOdsIq6pwfPofpQSRIsXw5f+xrExCjpdhcuwMGDyOnptzl0NzUaMmbO5GxCAjUj/f6QJOqcnakbJG3RxsaGiIgIIiMj8fb2Ns4sq4gIOHtWUV68fh06OyEoqN8iqQqYNWsWUVFRZGZmkpqaSnV1NSeTkzk9dy7hV66QkJ6uExfpValos7SkzNeXYj8/Svz9b4s2DYa5uTmrV68mPDx8DCes4OnpyQsvvMCuXbv69TPrNDfn+IIFuvvWLS14XbuGT1kZPmVluNTUYNHX1w0lgrz5vfdosLe/TVCjH9/8Jvy//6c4xXfC3h7+8z+VBZXTp5Vm9/v3K4v2A9NCAwORq6qQWlsBZQFjzT/+QbulJUVatUdZlqmoqLg9Sj8YsozXtWssO3gQP73vn1ovLyzffZctdxO2kaQRCa4IJ20U9ERHK1/wKAqPF2F4b7Lg3tF30gyQ6mcIAgICKCkpodTPjy61GtPubixKS+HqVTBAnVJffZhTTQ2u2pW8bo0GtaFEIExMaF60CHNtukFYTo7OSTM0A1Md80JC8DHQj82M6GgKAwMJ1QoNhF+5Ql50NAv0fkxGRVoaaC+s6h0caNZbCTOmhRmVSoWvr69OcKHJwYFSX19lVbynR0nV+ta3JniUU4e+KNqSw4d14gGlYWH4PvnkRA5rSuPg4MBTTz1FXl4eWVlZ5OXl3TGyplKp8PHxISAggICAALy8vG5dSK1YodzuFWtrZYV8yZI772digmVCAvHvvMPut94i8LXXmKUnQgTgqv2OTxxQNwcQd+ECcQOky+sdHCjz8aHeyYlaZ2fk+Hi8k5OJbmzE5he/6P97qU0FIzkZUlKUejt9mfbQUFi/Hn7xC6TqaiU1q62NHk9PrqvV1Nva0tvbS1hbGy51dVRWVt6TQImJiQk+Pj74+/sTGBg4cb3Z7gVzcyXScwfUajVxcXHExsZSWlpKeno6ubm5XI6O5vIo0vqtra2Ji4sjMTFx1IuaI8Hc3JwNGzbcMYLdamNDbng4uX2/5bKM5Y0b+JeUcP/+/VhrhWWGdNC8vZV0xZFmH6nVii2npMArryhpl5mZyoKDvb2y3cMDqbMT+eRJ2p95BsuSEkx6e3n0ww95e+tWKoaxiOtUW0vElSt4XbuGV3m5bj59dDk54XDkCM7BwSMb/3CmaPAzTiM0c+bAm28CtxQehZM2hlRWQnGx8r+5uVI3YAQEBARw+PBhuk1NKZ0x41Ze//79qPrq50ZBn7Okn+rYmpyMvQGFSUzXr9flhIfl5HBw+XLKysoMXtNVXFTE43qpjlciI5lpAJUlgKCgII5HROictMjLlzlVXk5zczO299Do9jb0Uh1LBjhlxuSkgRLZ63PSALJiYm6lLm3fLpw0A9HR0UFGRgae167pUm0Bmn/4w8lz4TlJkSSJ0NBQQkND6ezspKioiNraWrq6uuju7kaSJJydnXF3d8fFxcUo+mmBUvKx5rnn+NLPjzc//ZSQnBz8iovxKi8fkSDH2dmzObBiBaa2tsTGxrI0Pr5/TdvKlYqyZUWFooQXFKT8bg4HV1dFuRIwAXy0t4F0d3fT3d2NqakpKpWKmzdv0tTURJO2XsjGxgZra2s0Gg29vb3IsoyZmZnRvBdjiSRJ+Pn54efnp2ttkJmZyfVhCo6oVCqcnJxwc3PT9YWbqNetL4IdEhLC3r17dSI+dziANmtrLkdFUezvz8YjR5ih7SGJqSkkJCiL7LGxynVcaKhhenhaWSnnHbiAr9EgLV6M5bFj9MyZg8n162i6utj69ttciYggY9YsSmbM6KdQad/QQGhuLtGZmXgN8Z7JajU8/zym//VfY5adIpy0UWCv1yfFraoKqbeX2traMRNcmPborwrGxw9eMzABeHh46NLLCgMDdU5ax86dWI7SSevu7tYpCOk7aSaGSnXUYrd2LZ3PPoumsxOn+nqca2updXGhpqbGYNLy9fX1mObl6er2Ok1NyQ8JYbWBUgWdnJyonT+f7t27lTSfigoc6+rIyckhMTFx9E+g76TppYE6ODgYxgk0IAOdxssRETywfz8mPT2Qng45OUrTX8GoOHv2LF2dnSw7eFC37Up0NMHaC1zB+KDRaAibRPasUqlYsWIFjUlJXLp0iX2XLlFfVoZ3WRl+paX4lpai6u2lxN+ffO3qfOSlS4Tl5CDJMkcffhjWrWNtSAjBwcGD9/GSpDHPNlGr1f2e29zcHHNzc6NoR2JMWFpakpSURFJSEs3NzeTl5VFQUEBzczOyLCPLMmq1Gnd3d7y8vPDy8sLJycnonFk7Ozsee+wxrly5wokTJ6gc2PpiAD4+PqSkpOD7q18pWWddXYqDNlHq7j4+mBw8iJycjFRfj6ari1kZGczKyOCmRkOtszMNDg64VVXprlMGQ7a1RVq1Cumll2AMomf6CCdtFFgEBNBqZ4d1UxNmnZ24VldT5e5OZWWlUchxTzmMMNURlB9cf39/srOzFSn+L74AQHPsmJK7PgpnsrS0lK6uLuwaG/HURml7TEywfvRRg4y9D5WlJRUzZzJDm2YTlpPDCRcXrl69arAf3JycnH6OZl5ICG7+/gZtLBkQF0dBUBBhWrXKqKwswzhpHR397E+/Hs3YomgA7u7u/erSOiwtyQ0JuRXt2bZN6ckjuGe6u7tJTU0lsLBQJ93eK0lUf+tbRAiFX8EwsLe3Jzk5meTkZOrq6rh+/ToVFRUcr6qis7MTS0tLnC0tsbKyQrKxocLWFnt7ex5yczM+5UPBsLC1tWX27NnMnj17oodyT0iSRGRkJJGRkdTU1JCVlUVhYSEdHR2YmZlhZmaGg4MDMTExzJgx41ZGgbFcs4WHI33xBWzdqiiyajHr7MTr+vXBo2amporE/yOPQGIiUkjIkH3hDI1w0kZJQ3g41lqxCJ+rV6lyd6eiokI4aWPBQNEQI6LPSatxcaHRzg77pibUN24oxazaYvV7oa+/iL5zUxMTg/sYqF+1LV+uRFmA0JwcTqSkUFZWRkJCgkHOn52dzYN688gODydEr9+RIQgPD+dkdLTOSYvOyuJ4cTFtbW2ja4R6+rQiHALUOjnRqhc5MybRkD5UKhVBQUFc0vsRypw585aTtn270ojUyFZqJxOZmZm0NTWxQk9K+UJcHDHr1k3gqASTFScnJ5ycnIgeZTsSgWC8cHFxYcmSJSy5Wz2msREfr9SunT2rLFju2KH0fdPH3Fyp23z4Ydiw4baG9eOFWIoZJZ16qyE+2tqhu4WABfdAZ6fygerDWFZltOi610sSWTExuu09b799z+eUZVnnpOmrOvaMwum7Exbr19OrXfXyLi/HurmZwsLCIZtojoSWlhZaMzNx134RdpuYUBAcbHCFKnd3dyoTEujUyhI719XhVlHRrz7rnphE9Wh9DHSA84OCaOsrOL92DY4eHf9BTRFkWebUqVMkpKfr0mJuajRcf/HFYfW6EggEAsEEIklK6uXvw5kWFwAAIABJREFUf6/UbVZVKb+Jb74J+/Yp0v779il9eSfwO104aaNEo1eX5qstzB9uYahgBFy8qKScAfj5KQ0OjQhHR0fstLKqGXpOmrR3LzQ03NM5q6qqaG5uxrG2VmdbvZKE09NPj37Ag+AZGUmpXq1VaG4ubW1tlGhTuUZDdna2TmYaoCggAFsvL5wHkWQeDZIkERQTQ45efUpfyuOoGKIezd7eXve+GxtBQUH9xCt61Wqy9BvAv/POBIxqapCbm0v71ass0rOLrxYuJN6ADewFAoFAMA5IkiL8sXAhPPssPPCA0gfPCBBO2ihxXLJEt2pv39SETXMzNTU1d5QDFtwDRlqP1ockSbpoWp2LC+WengCourruuZdbn2Mx7/Rp+i61r8fEYD5GqbQajYZ6Pbn6SK0KY1ZW1qjPnZOT0y8amB0ePmZ9XsLDw7mklzIUdekShfn5dHZ23tsJ29p0rTagfyTNWKNoABYWFrelXWfoK6J+8gkM0oRUcGdkWebkyZMsPnJE14+qztGR6k2b8NR+7gUCgUAgGC3CSRslVvb2VOldCPlcvYosy1QNzG8VjA49J63NSKT3BxKgJyWfOXPmrQe2bbun8+Xl5WHV0sJMvQa5zS++eM/jGw7WmzfT18LUv6SEGcXF5OTk0N3dfc/nbGtroyYrS5cO3CtJ5IWGjpmT5uPjQ0VUFG1aBSm75mY8iospKCi4txMeO6aoUqE0m71hba17yNA1dYYmeIDyVIWHBw19jkRbG3z66QSManJTUlJCV3o6cefO6bYdWLGCuQsXTuCoBAKBQDDVEE6aAWjWW7X31V6Iin5pBkbPSXu3sJBXX32VPXv2TOCAbsdfLw3uUlQUPX3qP6dOgVaWf7g0NzdTUVFBUmoq6p4eAK55eeGxcaPBxjsYAQsXkhUbq7u/7OBBOtrbKRzh+PXJzc0lNCdHFw286uuLxssL9zFKWVWpVARHRpIdEaHbNjMjg8zMzHs74R//qPu3WM8RNzMzu80JMjZCQ0P7b5AkzumnPP7tb+M7oEmOLMscOXiQh3bv1jWuLggMpDklpd8ijUAgEAgEo0U4aQagVy/9zls4aYanvBy0NVldajVVbm40NDRwTm8l2xiwsrLC29sbgDYrKwqCgm49uH37iM6Vm5uL5uZNErRqiwCXVq7EYYwLWE1NTSl/4QW6tap/3uXlRFy5wmW9BtQjJTs7+7ZUx7CwsDFt9hseHk6mXm3gzIwMys+fp7GxcWQnysoC7WKAjNJAVv85Bu1PZEQ4OTnhMEAJNCM6GrlvAeHYMfjyywkY2eQkLy8Pz08+0ck0d5uYsP+BB0hOSRHNqwUCgUBgUISTZgAs9ORHPSoqMO3sFE6aIdG7iCz38qLXiGXDZ+mlYuqnPMrbtsEIUgbz8vKIO3euX82LxsANrIciYNEi0pKSdPeXHDpE3uXL91Rn2dHRQfnly/gXF+u25YSFjVmqYx+BgYE0RkXpagPVPT0knjnDWX2F0OHwyiu6f3PCwqh1cdHdj4qKMshYxxJJkm5LyWy1taV40aJbG77/fdBGa8cDWZapra0lMzOTY8eOceTIEQ4ePMiRI0fIzc2999rBMUaWZdJ27GDJ4cO6bccWLsRU2zNIIBAIBAJDYtzLwJME97Awql1ccK2pQSXLeJWXU2pmRkdHB+aiqeno0ZOxzzfy9LKoqCgOHDhAV1cXuSEhtJubY9HRgVRcrNSmPfPMXc/R0NDA1exsHtJL8Tw1bx5xeul7Y0lQUBCfL1tG7PnzWHR04FRfT/SZM+Tl5Y34YvT8+fPEnDuHiVbG/7qHB73e3vj4+IzF0HWoVCpmJyRwMjmZjVrhloS0NP586hQLFy7EVCv2c0eKi+HDD3V3TyQn6/63srLql95qzISEhJCqJ3wCsCs+nu+cOYPU1gYZGfDuu/DUU2M2BlmWycvL49y5c5SWltLd1oZLTQ1O9fWYt7dj3tEB3d0Um5uTa2mJVVAQjmvWEBkXh2YUzeANyaWsLJK2bUOjXayocnXl5Lx5PLpkiYiiCQQCgcDgCCfNAFhbW1McFIRrTQ2giIeU+PtTVFRExDhdWE9ZCgvhq68A6FWp+gtyGCFmZmZERESQkZFBj6kpp+fNu7Xy/l//BY89BlpBi6E4duwYC44cwVarvNdqZUVRcjKrxkk5zsTEBP/4eE6kpLBcG8VcfOQI+0+cGJGT1tvbS/qZMzyp5yCcj48nIiJiXC5q4+Li+G1kJHWOjoozcPMmESdOcOnSJWL16u6G5H//VxdhKvL357o2lRUgIiIClWpyJCLMmDEDc3NzOvpaWABNVlZcXb+eGX2iNj/+MWzceFfbHCldXV2cP3+e1NRUWisrmZOayqLsbFyrq3W1lkPR9vvfc2H2bG5u2ULk2rU4OTkZdGwjoaenh6b//E+S8/MBJfV1z0MP4ennZ/R1iQKBQCCYnEyOq4xJQPecObr/+1Ts8rU/6IJRoB9FCwqi1cZGd99Y64H0HYAzSUm06jcQ/tOf7nhsQ0MDFV98wVy9KNrB5cuJnTt3XFfro6KiSE1MpMHeHgDL9nYCXn99RGm8OTk5uKel4aCtA2uzsCAjJobExMQxGfNArKysiIyJ4eT8+bptc0+f5uzJk8iyfIcjUZpb/vWvursn9aJoANF6YkHGjomJCXFxcbdt/yw4GNnVVblz7ZrilBqInp4e0tPTefV3v+PQrl34f/kl3/r971ly+DCeFRV3ddBAsbmk48dZ8PzzVKekcOy117hx44bBxjgScv/930netUt3Pz0hgWs+PixdulRE0QQCgUAwJggnzUDY6zUx9SkrQ9XTQ35+/t0vBgVD09PTz0m7OEB6f8aMGeM8oOHh6+urW/XvMjPjmL409y9/CU1NQx57/KuveHDXLp1yXMmMGeTNmUOSXo3YeODn54e1szP7H3hAty32wgWu/OUvwz7HmTNnSDpzRnf/XHw8AVFR4xoRSUxMJHPmTFq0svk2ra347N1LmXYhZVC6upSIpzbydN3DgyI95T47OzudQMxkITEx8TZnorGnh7Kvfe3Whp/8ZMQCN4NRfO4cOfPnE7FoEf/8gx/wH7/8JQ/t2YNNa2u//Rrs7ckJDeV8bCyn5s7lWEoKaYmJZEZH0zigQXh4djYp3/gGxcnJnPnHP8a1D2XtW28RpufAFvv58cV99xEYGGjUffIEAoFAMLkRTpqB8EpJoUl7YWF+8ya+paXcuHFDCIiMhkOHlBV+oMPGhrwBAgjGKnktSVI/AZFzcXG6iBT19fDyy4Me19DQgOqtt/DRzrlHpWLPqlXMmTsXMzOzMR+3PiqViuTkZPJDQ8kOC9Ntj/rTn6i8k4Ojpby8nK7UVPxKSwFlLumJiczRiziPB56enrjPmMEZPQXW+w4c4MrLLw/d++1739Ol2MqSxMFly0DPwYmKipp00RM7O7tBU1X3eXgg9ylWyrJSl3aPjlp3dzdpr7yC0+LFRKamYtXWplts6KPZxoadq1fzfz/5Cbt++1tKfvtbmn/zG3peeYWGf/1Xjm/cyD/WrePV73yH9zZvJlfvMy8BURcvErVlC59+97ucP3+eXm2t41jRfegQdi+8gEr7PJVubny0aROymRnLli0b0+cWCAQCwfTGOPPFJiEaMzOKExKwO3gQgLDcXEoCAsjPz8dznGqJphxvvaX7NyM6mt4B6Y3G6qSBovJ4+PBhZFmmV63myJIlrO1rHPzKK9DcDL/5DWiFZeTubsq//nVW7t6tO8fJ5GRu+PiMexStj9jYWI4fP87nDzxAYGEhmq4u3KqrufS97+H+0Ud3PDY1NZU5elG0y5GRWAYHT0jkITExkT3FxURnZuJeVYVKlln2xhtc9Pdn9g9+0H/nbdvg1Vd1dw8vXkxxYKDuvkajmbD3Y7TMnTuXS5cu9dtWVVdHyZ/+hP/XvgaZmbcctddeAwcHsLNT6tTMzMDSEpycwNkZ7O2ht1eJdjc303LxInUnT5Jw6RID3dcutZpWa2vOx8VR/fjjJCxcyOrAwEEdXVmWKS8v5/z582RpNBQEB+N+/TqLjxwhRJs+bn3jBhtee40vCwt57cEHWbp8OaGhoYZ3nE+ehFWrMNVG7Rrs7XnviSe4aW7OogULxqzPn0AgEAgEIJw0g9K7ejVonbTQnBw+v/9+8vPzWaif7iYYHnV18I9/6O6e1+t5BUq9kZub23iPathYW1sTExNDRkYGAFlRUSSkpemiZLz2mtLkevNmZFNTGv76V6L0LqDrHB05npJCygRE0fowMTEhOTmZffv2cXTRIu7TiohEffwxzQkJ2H7ve4MeV1hYSNO+ff3mk5qURFJS0oREoCIjIznh5cW7TzzB02+9hVN9PeqeHmJ+/GOaCgqw27JFaY/w2mvQ50gDV8LDOZGS0u9cKSkp2OjVRU4mPD098fX15aq252AfX164wDMHDqBeseKWo6ZXEzkcbLS3PtosLNi9ejV5wcH0qtVERESwcOFCXPtq4IZAkiS8vb3x9vbmvvvuIy0tjZMnT/LB5s0EFBay9pNPdBG6FQcOMOf0aS59+CE7ly8n7pln8PX1HdG4hyQtjZ4VK1BrU15brK3ZvmULrTY2eHp6kjygRlEgEAgEAkMj0h0NiNvGjXRoL6jtm5pwq6ykvLx8word9enu7qa0tJTjx4/zwQcfsG3bNvbt20dRURE949gjadj88Ieg7RHWHBxM9YBV64CAAKNPOVuxYgW2trbKHZWK9554gsv6ap8ZGfCDHyD98z/jqOfQlPr6su2pp1Db2Ex41CY2NhZbW1tS58zhqp5svu33v0/7n/982/4NDQ2c/O1veXz7dp3sfqmvL02hoRMmtmFiYsKaNWtot7Nj+5YtNGudLE1XF3ZvvAELFsCSJbBjh07NsdbNjZ2PPNIvzdHBwWHc0zUNzWDjr6ioYG9qKvLBgzBv3qifoyAwkNe+8Q1ywsPxDQri+eefZ8OGDXd10AZibm7OggUL+Na3vsXs2bMpDgriLy+8wDUvL90+ds3NzD91ikf++7+5sWIFu/74R2pra0c3gQ8/pGfpUky039utVla889RTNDg5YWJiwiOPPIKJEfdqFAgEAsHUQETSDIijmxu50dGEahvmhuXkUOXhQUFBATMnSDpelmUuXLjAoUOHaGtr6/dYcXEx6enpmJubExsby+LFi4fXP2qsOXCgn7LemcWLb9vFmFMd+7CwsGDNmjW88847ANw0N+fvGzZQkp7OigMHblO465Ukji5axImUFGSVitX33TdhUbQ+1Go18+fPZ//+/by/eTNbtm3D6/p1AMy/8Q1uVlVh9sMfgrk5XV1dHPr1r9nw5puYaRsSt1pZsWv1aubMmTOhapweHh4sXryYQ4cOsf3JJ9n87rvYDyHgUuzvz87Vq+kc8NqvWLHCaBVFh0toaChOTk7U1dX1237x4kW8vb2JP3ECcnOhpgYaGxWRm44OZcHkxg2oq+NGaSnV+fm0dXTQq1LRrVbT6OBAnaMjNa6uVLu6YqJWc9/SpcyZM2fUiynW1tY8+OCDJCUlcejQId62smLRkSPEnT+PZXu7br/wnBwC//mfOfnpp9jOmUNkfDzmrq4QEgIuLorD3dYG2dlQXQ0mJqBWg0YDNjZKSudPfwoffECfC9ZmYcG2J5+kTtvEfMmSJbjoNTQXCAQCgWCsmNxXHEZI+333gdZJC83N5avFi8nPz58QJ62jo4M9e/Zw+fLlu+53+vRpiouL2bRpE3YDlNXGlaYm0FOc616zhtOD1PRNBicNFJXElJQUjh8/rmyQJM4mJlLs70/U5ctobt7EpKeHHrWaS1FRVGjnumjRouH18hoH4uLiOH36NI2Njbz7xBM8uW0bHpWVSLKM2Usv0f7nP1O5ahUmX33F2oICnVhE3wWuXUIC8wwQoRkt8+bNIz8/n6vAH775TfyLiwnNyyM4Lw+Tnh6yIyJInz2bmkHSaAMDAwkZIFwzGVGpVDzyyCO8/fbbt0XQ9+3bh0qlYtasWUh6YjF9VFVVcezYMa5YWkJ4+JDP4ebmxiOPPGLwmi1nZ2ceffRRyubN42BAAEcWLyawqIiZFy8SeeUKoERHFx8+DH29CftwdARbWygtVdI5h0GDvT0fbdqks4eYmBjm6gnQCAQCgUAwlggnzcDYbdxI9yuvoO7pwaOyEruGBvLy8mhvb8fCwI1i70RlZSUfffQRjdoeVcM95o033mDjxo2Gq+0YKd//vk7REWdnzj/zDJw7128XZ2fnW2mEk4CFCxdSXFzMtb55AXUuLny1aNGg+8+bN48FCxaM0+jujlqtZuPGjWzbto0OYPuWLWzZvh2PykoALCor8X/zzX7HtJubs/3JJ+kMCWH9+vVG0fi5z0F5/fXXuQkUhIRQMAzHy9LSkpUrVxp9eu1w8fb25oEHHmDPnj39tvf29rJr1y4uXLhASkoKbm5uqFQqCgsLyc7OJjc3947n1Wg0LFq0iMTExDFNB/Tx8WHr1q3k5+dz8OBB/h4SQlppKQ/u2YNrTc3gB9XXK7dhciE2ls9XrKBTK+wTGhrK6tWrp4wNCAQCgcD4EU6agfGJjKQ0MJDAvDxAiaalOTiQnp4+bhfeV69e5f333+emtqZrJNy4cYN33nmH9evXE36H1fIxYds2eOMN3d2eV1/l+CAXhkFBQeM5qlFjYmLC5s2b2blzJzk5OXfcNzExkWXLlhndxaCHhwdbtmxh27ZttANvPvcc8WfPsujo0X4pZzJQ4ufHgRUrqPPx4ZlHH8XS0nLCxj0QBwcHnnrqKXbs2EFDQ8Nd9/fw8GDdunU4OjqOw+jGj/j4eMrLy7lw4cJtj5WVlfH++++P6HyRkZHcd99947Z4IkkSISEhBAUFkZGRwZEjR3jd25u48+fxLS1F3d2Nursbqxs3cK6tRaNVaJRNTOj280MdEACyTE9nJ2319dysq0Pd3k6bpSXHU1LI1fvu8/PzY/369aIOTSAQCATjinDSDIxarVZSHrVOWlhODmlz5pCamsrcuXPHvOarqKiIDz/8cMhmr8HBwQQGBmJubk5+fj45OTm3pT319vby97//fXwdtc8/h2efvXV/3TouBgfTqn0d+5AkiYSEhPEZkwExNzfn0Ucfpby8nKNHj1JQUKB7zMzMjICAAOLi4ozaAfX09OSJJ55g+/btdHZ2kp6URGZMDPNOncKjooKigAAuR0XRYmurRK1Wr8bDw2Oih30bHh4evPjiixw4cIDz588Pud+cOXNYunTppK9DG4qVK1dSXV1NeXn5PZ/D1dWVBx54YMKaOqtUKmJjY4mKiuLMmTOcsLDg7MDvB1nGtrkZzc2bNDg60qNWI0kS8jDSHv38/Ni0adOUtQGBQCAQGC/il2cM8Pmnf4I//AEA/5ISfEtLuTpjBhcvXhxTByM3N5cdO3YMqtao0WhYtWrVLYW9lhZmShL1LS2cOXSIXFdXmvVq0cbVUUtPh/XrFRl0gJgYet94g1PvvnvbrhEREZM6quHl5cXmzZupqqqioqICR0dHvL29jSIdcDh4e3vzzDPPsH//fkpLS7lpYcGRpUv77RMSEkJKSgre3t4TNMq7o9FoeOihhwgLC+PEiRO0tLRgbm6Oubk5Dg4OzJo1Cx89NcupiFqt5oknnmDv3r239U+7GzY2NqSkpBAfH28UtmtqakpKSgqxsbEcOnSIixcv3npQkvp9twHDctASEhJYsWKFiKAJBAKBYEIQTtoYYBcWxrW5c/HW9hpauXcvf3nhBU6fPj1mFzV5eXl8/PHH9Gplz/VxdnZmc1IS9l98AT//OaSlQUkJAI7ASu2t3NOT7PBwLs6axQ0bG52jtmHDBsIGERIwCCdOwJo1inIcgK8v7N9PTkUF9YPUkMyfP39sxjHOuLm5GXWftzvh5ubG1q1baW5upqCggMLCQlpbW/H19SU+Ph57e/uJHuKwCQ4OJjg4eKKHMWGYm5uzbt064uLi2Ldv313l621sbEhOTiYuLs4oo0vW1tY8/PDDzJ49my+//JLS0tIRn0OlUrFy5Uri4+PHYIQCgUAgEAwP4/uVnSKY/eEPdM6Zg6arC7fqahLS0kidO5crV64QFRVl0OcqKCgY0kELUalYf+oUpt/+NgzyuD5e16/jdf06i44eJTMmhlPz5lHn4sKOHTt49NFHDa9u9/rr8K1vQV9qpqMjHDhAr7s7JwaIGoCisGeM6XPTFVtbW+Li4oiLi5vooQhGib+/Py+++CJZWVkUFxdTU1NDbW0tXV1dODk5ERwcTEhICL6+vpMisuTl5cXWrVspKSnh6NGjw3LWJEkiPDycBQsWTNoFFIFAIBBMHYSTNka4xMVxYe1aYj/6CIDFR45wOTKSQ4cOERgYaDClx4KCAj788MN+KY6mnZ0E5+Uxu7gYv4sXkQZrVm1qCv7+yq27G776SpduqO7pIe7CBeIuXCArKopDS5fy8ccfs2nTJsPUTNXVKc2q9Xqh4eICu3cjh4ayd88eKioqbjtsqkTRBAJjxMTEhFmzZjFr1ixASQmUZdko0hnvFT8/P7Zu3UplZSVFRUWUlJRw9epVnaiSRqPB2tqa4OBgkpKScHBwmOARCwQCgUCgIJy0McThpz+l5tAhXGprMevsZOW+fXy8cSOfffYZmzZtGrWC3/nz59mzZ4+uvkLV08OygweZnZ6OaV99lz7Ll8PatZCYCFFRShPXPhoaYPduev/4R1RpabrN0ZcuEZ6dTVpiIl/U19P+9NO36tpGSkcHvPoq/PKXSj+0PmJj4bPPwNeXw4cODSrm4OnpOWHiBALBdESSJKNTGb1X3N3dcXd3Z968eciyTE9PDyYmJlNmfgKBQCCYekzeJdJJwIzgYM49/bTufnhODiv37SMvN5cTJ07c83llWebQoUPs3r1b56CZtbez+d13mXv69O0O2pIlSu3XF1/Aiy9CXFx/Bw3AwQGefBLVmTN0HT5MmV7zbXVPD/NOn+Yb//d/uC5bRv7jj9P91Ve30hT76O6+dbtxA8rKICNDkdZ/9FHw8FAiaPoO2qZNcOIE7S4ufP7550O+LkuWLBEXVAKBYNRIkoRaq/AoEAgEAoGxIiJpY4gkSST+279x7vhx4s+cASDh7Fm6NBoOShKWlpbExcWN6GKhtraWPXv29KuxcKqp4dGPPsJFv+g/IgI2boQNG5T/hz9oTBcvxjU1lb0vvUTMu+/io9eE2a26GrcPPoAPPqDX0hIpJgapvh6uX4fW1uE/T3Aw/OpXNC1ZQsbZs5w6dWrIvm73338/gYGBwz+3QCAQCAQCgUAwiRFO2hjj6OhI9ZtvkvnYY8RkZQEw79QprFpbOVldTVZWFitXrsTV1fWO52lqauKs1pnp7e3FrrGRmMxMwrKz8RxYv/Wzn8GPfgSjWCk2MzNj2UsvsT0gAJuDB5mZkUFQQQFqvfo2VVsbaJ3P4dLj7U3Vk09yJTmZgtJSqn73uzvuv2DBApKSku5pDgKBQCAQCAQCwWREOGnjQFhkJIf+53/I/u53Cc/JAWBmZiYzMzMpDAjg1M6dNM6di+/s2Tg5OaHRaDAxMaG1tZXm5mZKSkp0kTPbxkYWHD/OrAsXMBmg1iibmSG99RY89phBxm1mZsaTTz3FbkdHPoqIwPTmTYILCggqKMC/qAh7/bRFQJYkZK1jKKtUdFhZ0WFlxQ0bGwoDAsjx86Pa1VVxHvXq3oYiPj6eRYsWGWQuAoFAIBAIBALBZEE4aePE4uXLef9734Nf/1rnqAEEFhURWFSE/NlnlHt5UeXmRrWjI012dqh6ezHt6sKro4Po+noc6+vxvXr1NuesV61GWr4c6aWXFFEQA6LRaFi7di0+Pj4cOHCAK5GRXImMBFnGoaEBu8ZGblhZ0WJjQ4eFxaiid32o1Wrmz5/PggULRN2IQCAQCAQCgWDaIZy0cUKlUrFxyxb2OTpyZtcuklJTCcvJQaUV/pAA7/JyvMvLh33O8sBA1N/8Jm5PPw12dmM0cm1tXWIinp6e7N27l8rKSpAkGhwdaXB0NOjzxMbGsnDhQmxtbQ12XoFAIBAIBAKBYDIhnLRxRKPR8MiaNRTNnMmePXv4oqiIqEuXCM7Px6esTOew3Y0yHx9qv/1tor79bUwHqjSOId7e3jz//PPk5eVx7Ngxrl+/bpDzuru7ExISwsyZM3E0oNMnEAgEAoFAIBBMRoSTNgEEBATw9a9/nbS0NC5HRHCyogKLtja8r13DQZvWaNPSQrdaTbdaTadGQ6ODA11+frguWEDo/ffjM0FNVyVJIjQ0lJCQEMrKyigqKqKsrIyysjK6BkryD4KlpSWurq64ubnh7u6Ov78/dmMYBRQIBAKBQCAQCCYbwkmbIExNTZk/fz7z58+nqamJvLw8qqqqaO/spOjmTXp6erCyssLGxgYbGxtm+vri7u5uNDVakiTh6+uLr68vAL29vdy4cYOOjg46Ojro6upCrVbrblZWVlhYWKBSidZ8AoFAIBAIBALBnZDkYabYjfqJJKkGKL3rjoPjDNTeda+pw3Sa72jnOkOWZRdDDaYPYa8jYjrNV9jr5Gc6zXc0cx0TWwVhryNkOs1X2OvkZjrNFcbBXsfNSRsNkiSdlWV59kSPY7yYTvOdinOdinO6E9NpvlNxrlNxTndiOs13Ks51Ks7pTkyn+U7FuU7FOQ3FdJorjM98Re6ZQCAQCAQCgUAgEBgRwkkTCAQCgUAgEAgEAiNisjhpf5noAYwz02m+U3GuU3FOd2I6zXcqznUqzulOTKf5TsW5TsU53YnpNN+pONepOKehmE5zhXGY76SoSRMIBAKBQCAQCASC6cJkiaQJBAKBQCAQCAQCwbRAOGkCgUAgEAgEAoFAYEQYnZMmSVKIJEm/kyQpU5KkVkmSKiRJ2iVJ0sxhHv+2JEnyILffjvXY7zIuH0mS/i5JUpMkSc2SJH0qSZLvMI9D9n32AAAGNUlEQVQ1lyTpf7SvRbskSaclSVow1mO+VyRJWi9J0ieSJJVqx5srSdLLkiTZDOPYwd47WZKkWeMx9pEi7HXQY4W9CnsdV4S9Tj17naq2qh3btLDX6WKrIOz1DscKex0F6tEcPEbcBywG3gHOA/bAD4BUSZLmy7J8bhjnqAFWD9hWYdBRjgBJkiyBw8BN4ClABn4OHJEkKUaW5Rt3OcVfgQeB7wNFwD8BByRJmivL8sWxG/k98z3gKvAfwDUgFngJWCxJ0jxZlnvvcvzbwOsDtuUZeIyGQtjr7Qh7FfY6bgh7nbL2OuVsFaadvU4XWwVhr0Mh7HU0yLJsVDeUDt7SgG12QAOwbRjHvw1cm+h5DBjTd4AeIEhvmz/QDfzLXY6difKheFpvmxrIBXZN9NyGGLPLINue1M5jyV2OlYGfT/QcRjBXYa/9jxX2asQ3Ya+3HSvs1UhvU9FWteOaNvY6XWxVO15hr7cfK+x1lDejS3eUZblW1s5Yb1sTijfqNTGjGjWrgTOyLBf0bZBluRg4CTw8jGO7gI/0ju0GPgRWSJJkZvjhjg5ZlmsG2Zyu/TtZ38NBEfY66LHCXo0UYa+DHivs1QiZorYK08hep4utgrDXOxwr7HUUGJ2TNhiSJDkCUUD2MA9xlSSpVpKkbkmS8iRJ+qEkSSZjOMS7EQlcGmT7ZSBiGMcWy7LcNsixGiBo9MMbFxZq/w7nPfy6JEk3JUlqkyTpsCRJKWM5MEMj7FXY61gOzNAIexX2OpYDMyRTwFZB2Ou0sFUQ9oqw11HbqzHWpA3G7wEJGE4B5UXgHIoRmANrgJeBYOBrYzXAu+CIEvIeSD3gMIpj+x43aiRJ8gJ+ChyUZfnsXXZ/F9gDXAdmoOQxH5Ykabksy0fHdKCGQ9jr4Mf2PW7UCHu9I8JejYxpZq+T3VZhGtvrNLNVEPYq7HW09joOOZ7LUHI173Y7OsTx/659/JlRjOH/tOcIHuv5DvH8ncDLg2z/BdB9l2O/BE4Psn25dk4pEzGnEczdGjirNVzvezjeBigFTozTeIW9CnsV9jqyMQh7naDbZLJXYau6MUxLe51Mtqp9PmGvsrDXibbX8YiknQLCh7HfwHAokiS9CPwS+LEsy38bxRg+AL4LzAbyR3Gee6WBwVcMHBh8lUGfemAwqVMHvceNEkmSzIFdQACwUJblayM9hyzLLZIk7QWeNfT4hkDYq7BXYa8jQ9jrBDAJ7VXYqsK0s9dJaKsg7LUPYa8TaK9j7qTJSi5qzkiPkyRpC/An4H9lWf7FKIch9Q1nlOe5Vy6j5OYOJAK4Moxj10iSZCn3z+uNQFnhKBj8sIlFkiRT4BMgEVgmy3LWaE7HOL13wl4BYa/CXkd4ur7hjPI894qw10lgr8JWdUwre52MtgrCXvUQ9jqK0zHa926iQ4pDhAnXoMh7/sVA5/sd0AsETtB8vqudT4DeNj8U1Zt/vcuxs7Rv8lN629QoRYy7J/q9GmLMKuBjoANYOspz2aL0rfhqoud1hzEKe721n7BXYa/jPR9hr/d+LqO216lmq9oxTBt7nU62qh2jsNf+xwp7HaW9TviLMsjEFmhfoPPAPGCO3i12wL6HgAK9+zOAY8A3UBoLPgT8TWvkr03gnKxQVgyyUCRLVwMZKI39rAeMvxv4yYDjP0QJK38NWAr8XfsaxU30+zXEfF/TfjB/PuD9m4M2t3ewuaI0EnwDeBxYhNI4MQtlxcUoc5eFvQp7FfYq7FXYq7BVYa/T11aFvQp7HSt7nfAXZZAX6SWGLtAsGbDvUf1tKHmzn6EU63UA7doPzDcB1QTPyxclhNoMtGjH6TdgHz/tPF8asN0C+A1QqZ1XKrBoot+rO8y15A7v4UtDzVX7xXQSqEVZpalDyQtOnOg5CXsV9irsVdirMdymi71OVVudTvY6XWxV2Kuw17GyV0n7BAKBQCAQCAQCgUAgMAImRTNrgUAgEAgEAoFAIJguCCdNIBAIBAKBQCAQCIwI4aQJBAKBQCAQCAQCgREhnDSBQCAQCAQCgUAgMCKEkyYQCAQCgUAgEAgERoRw0gQCgUAgEAgEAoHAiBBOmkAgEAgEAoFAIBAYEcJJEwgEAoFAIBAIBAIj4v8D9Itw6CFaIlQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x115d1e0f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#fig, axes = plt.subplots(nrows=2, ncols=5, sharex=True, figsize=(15, 5))\n",
    "fig, axes = plt.subplots(nrows=1, ncols=5, sharey=True, sharex=True, figsize=(15, 3))\n",
    "axes = axes.reshape((1, 5))\n",
    "fig.subplots_adjust(hspace=0.05, wspace=0.15)\n",
    "#for row, torsion_index in zip([0, 1], [torsions_simple, torsions_complex]):\n",
    "for row, torsion_index in zip([0], [torsions_complex]):\n",
    "    for i, ax in enumerate(axes[row]):\n",
    "        ax.plot(xticks, torsion_hists_train[torsion_index[i]], color='grey', linewidth=5)\n",
    "        ax.plot(xticks, torsion_hists_gen[torsion_index[i]], color='red', linewidth=3)\n",
    "        ax.set_yticks([])\n",
    "        ax.set_xlim(-np.pi, np.pi)\n",
    "        if row == 0:\n",
    "            ax.set_ylim(0, 5)\n",
    "        if row == 1:\n",
    "            ax.set_ylim(0, 1.5)\n",
    "        if row == 1:\n",
    "            ax.set_xticks((-np.pi,0, np.pi))\n",
    "            ax.set_xticklabels(('$-\\pi$', 0,'$\\pi$'))\n",
    "    axes[-1,-1].set_yticks([])\n",
    "axes[0, 0].set_ylim(0, 1.5)\n",
    "axes[0, 0].text(0, 1.35, '$\\phi$')\n",
    "axes[0, 1].text(0, 1.35, '$\\gamma_1$')\n",
    "axes[0, 2].text(0, 1.35, '$\\psi$')\n",
    "axes[0, 3].text(0, 1.35, '$\\gamma_2$')\n",
    "axes[0, 4].text(0, 1.35, '$\\gamma_3$')\n",
    "axes[0, 0].set_ylabel('density')\n",
    "#axes[1, 0].set_ylabel('density')\n",
    "axes[0, -1].text(-np.pi+0.3, 1.1, 'Target', color='grey')\n",
    "axes[0, -1].text(-np.pi+0.3, 0.7, 'RNVP + MCMC', color='red')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.41\n",
      "0.20\n",
      "0.16\n",
      "0.08\n",
      "0.10\n"
     ]
    }
   ],
   "source": [
    "eps = 1e-5\n",
    "for i in torsions_complex:\n",
    "    KL_SNF = np.sum(torsion_hists_train[i] * np.log((torsion_hists_train[i]+eps) / (torsion_hists_gen[i]+eps)))\n",
    "    print(\"{:1.2f}\".format(KL_SNF))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}