Restricted_Boltzmann_Machine.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 🧑‍🎓 Students:\n",
    "\n",
    "- Diyun LU\n",
    "- Marwan MASHRA\n",
    "- Yihan ZHONG"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  Deep Probabilistic Generative Models - Restricted Boltzmann Machine\n",
    "\n",
    "\n",
    "In this lab exercise, you will code a Restricted Boltzmann Machine with Gaussian observed random variables and Bernoulli latent variables. The explanation of the model and the derivation of all formulas are given in the PDF available from the course website, lab exercise section. Please, read it carefully to understand what is going on!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats\n",
    "import matplotlib.pyplot as plt\n",
    "import itertools\n",
    "import random\n",
    "import math\n",
    "import time\n",
    "import sklearn.datasets\n",
    "\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "import torch.autograd as autograd\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# samples1-2 shape must be (n samples, 2)\n",
    "def plot_samples(samples1, samples2=None):\n",
    "    fig, ax = plt.subplots()\n",
    "\n",
    "    ax.scatter(samples1[:,0], samples1[:,1], marker=\"x\", color=\"blue\")\n",
    "    if samples2 is not None:\n",
    "        ax.scatter(samples2[:,0], samples2[:,1], marker=\"x\", color=\"red\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data generation\n",
    "\n",
    "The data we will use in the lab exercise looks like a ring.\n",
    "You can see the it would be difficult fit this dataset with a Gaussian Mixture Model: we would need to use many Gaussians and each one would imperfectly fit a subset of the data space.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/FklEQVR4nO19XYxl2VXed7qHmcEIuttTnm53wzCeBIUY4aq4q6tmnEyXlUTCduQ2JEaCBwKSkbHirqqu4sWWqaoeeEAQ6L8BJUEjBHkhDrxknB+PsLHHPNg9rpZ6/AMamLFiYeQZPME9KFKMu+ruPOy7fNZZd+2fc+4593d90lHVuefcffY59+xvr/2ttdcunHMwGAwGw+zjyLgrYDAYDIbRwAjfYDAY5gRG+AaDwTAnMMI3GAyGOYERvsFgMMwJ7hl3BUJYWFhwDz/88LirYTAYDFOFW7duveqce4N2bGIJ/+GHH8b+/v64q2EwGAxThaIovho6ZpKOwWAwzAmM8A0Gg2FOYIRvMBgMcwIjfIPBYJgTGOEbDC1BpqWyNFWGSYMRvsHQAi5fBra2SpJ3zu9fvjzOWhkMVRjhGwxDwjngzh3g+vWS9Le2/P6dO2bpGyYHExuHbzBMC4oCuHrV/3/9ut8AYHPTf14U46ubwcBRTGo+/OXlZWcTrwzTBOeAI2zM3OsZ2RtGj6IobjnnlrVjJukYDC2AZBwOrukbDJMAI3yDYUhwzX5z01v2m5tVTd9gmASYhm8wDImiAI4fr2r2pOkfP26yjmFyYBq+Ya7gXJWA5f6klj2O6ximE6bhGwyYjVj5WbgHw/hghG+YC3QdKz8KIrZ4f8OwMA3fMNFoKl9o5x07BiwtVWPll5b851qZudfmRAx4/Z6IeGOj+r1h5BeL9zcMC7PwDWNFLP9MU6tZ+96lS8AzzwC3b1fPvX0beO21wXrUuTYRMUXmHDni/66uDt7bsFY/J32Ckb0hF0b4hrEhRqpN5Qv5vV7P/71xA1hZARYXq+cvLQFXrgxa8nWvrRHx6qq/bpvyi8X7G4aCc24it7NnzzrD7KLXc25z0znA/9X2+We00bHcsmnb2HBufb36WazMutcOXXNjo379h3lmBgOAfRfg1bETe2gzwp995JBqr1c9nktqh4fV7x0eOreyMkj2S0vO7e6G6yevrXUMMSKWhD8sKe/tVZ8RXWtvb7hyDbMDI3zDSKGRYuzcECG2aeFLyztlHWtlrK6W3+Xn7O3pRLyx4b/TloXP6xbbN8w3jPANI0MdCzRG6E3lC3metPRDhB0rg8g71WHIzoq+M075xTqH+UOM8C0s09AaXCQ8cXNzMDyRH+Pn0nebpCvgaQ6uXAG2t8P1pTJlWVqqhGvX/LGbN70j9sYNvx8KiSwK4MSJ4dItOFeWpe2ncPmy/z2uXPGRQ/TMjx0DnngirwzDjCHUE4x7Mwt/+lDX0ZozGsixULVzDg8HLf26FnZIr6+jyUurX+6HsLdXlZC4RJSj1/PfYWmpev8xv4Vh+gGz8A1dg6zJq1f9RpY6EI4Tp/BLOqZZ3JrlHLpuUZRW7PHj7Sc0c25wxLC1FY+Dp895PaksqqeMy3cO+Nu/9aOJmzfLz5580v+/ulp9bqHrXrkCPPusn2tw9Kj/fGnJ76+tpcswzCBCPcG4N7PwpwcyKkVGpnDdXPtubF87Rn81qz2mq6fK55Cjj8NDbxk30eTr+iP29sJhpLFnqUH6MHgZpu/PJtC10xbA7wL4GwBfChwvANwA8CKALwB4a6pMI/zpgnRsSvIPERvJLVTG5qYuN9C5u7tViWZ3tyTitiJhYgRN8gg/r67Ekhv/L58nhZfSefJ72jVllBCVu7bmr3NwUL2uST3Tj1EQ/nkAb40Q/rsA/K8+8T8K4GaqTCP86UMobj0VCZPSmOW52t9cXT3Xqg0RNBFu6vuhMnPqqXWenKyp00v5PkJlAM4tLvq/CwtleW3o+zZqGD86J3x/DTwcIfz/BOCn2f4LAN4YK2/eCH+SG0quNRkLsdRwcDBI2FI2iZUvySsle9SdtFTXQas9G3mtnJFIrzco5/B9ut+Njbgzend38NnQdvGiJ/ucZ58LmxQ2GZgEwv/vAP4Z2/8kgGXlvPcD2Aew/9BDD3X8WCYHk9xQpDVJBCOtybox83TPBweDZBQjwpDFSiRPVu3qarWzqVvHOgQt70lek1vQudeWUszqqnMPPujcyZO+vNyUDdq5sa0p2Td5BwzdYGoIn2/zYuFPUkORFjXXx8maDFmBTSdcaZLMwUEeAacIX+uscvT+Jr9J6js5Mowsh555zDdCm5SaeHl1CD9UTg6adJKG9jEJhG+STgTjaCiybE5Ie3teQgjp5KS5S7mGPgtdg0MjX17++npJhjENX5MsNKmD7+fINE1GXanfMVe2i0UIpTo8joOD6jO4e9e5U6eq39NknZSfIuUjyXm+hu4wCYT/r4TT9rlUefNE+M5131CkNc4JghOKtB41ojk8DE8MWlnJnzi1s1Mt9+LF6vV4HWNROqEEZRoB14noaeJXaet3lN/TLH0aDUkZyznnzp8vHbBE9pzcT5+u7vNnz59JquOTslmb2UENzTCKKJ0/APB1AHcBfA3A+wB8AMAH+scLAL8N4CUAX0zJOW7OCL9rC583Wk3njskrIWuSOxFlJ7G+Xl6Lk0OoHtrGyYs/J/k3RTJaHLocJbT1rLW6yNw92v00LZueNWUBXVnx5x0clGS+tFR1ji8sOPfBD1bLoO9J/0yOTEWdPu+QyM9gGv54MBILv+1tXgi/aw1fK08jDrKaU2TPRwKh/PJacjEuqWgdhCwr15qOPbuQjCOPt+Ec547WHH9CjkTEy87R4eV1X/e6wd+O/COhZy3/DxkjcsQhM4NyWW0Sgg/mCUb4E46uo3S0Rqs5/WJO0fX18jskFUji0MrlmrCmRYcck0Qq8j60Z8etaHI0a7Hl8n6pzGEsb/6dGOEPkzmTl7246Ose6mz5dvHi4D3njiblaEzrIEIdUWxUY+geRvhTgDZIJ1V+ipg165scqEQgNPGn19NJR/tMklyoo5Ejg5xZrZzwNas65HjmDmpObPw6dX+DmKSjzQiuM8mJ7pNIW15HOq8vXhz8jP+WsTUB+OepVBlap28kP14Y4c84Up1FiGxDBEvkwsn9zBkf4UGSQEzD1zqS0L4WsqiFf0oLWcoNMbmKlxWSlsgKp336v+4oS7OGQ88/18LnZUsyDt0zafinTnnnON+njps6NdmRhnwjMlVGyJFcN9+PoV0Y4c8wcqIoJClK6YHOkekMNAKmTuF7vqcklrt3q45C6QxOWfyafKNJTNIZTJ2RJB05yojp5znEVoeQQ3KJJmdpIZA5kKOSmLzDpbiFBR+9wzs1KWmlpBz+7KWGry0Sk/PMYvuG+jDCn1HkOnw12YJbeFRWDKFY8HPnqo19cbGMoU+RXMpXIQlHy7nDCSw2ioj5AzRZYhiyl79FbASk1S33evSXOnB53w8+WEpAcuauNoLJlXJoC4Xm5ubsn+QZ5tMMI/wpRI5ME5ILNKJK6dU5kFE80rLkkTualRzSi0P3H3Mix5KnaZ1DiPRj16lrbcYILKTh7+wMT3Jc35cjm5g8w38LPpqT566u6vl7yGCQnWeuZd9ldNo8wwh/ypCyfKRlJRu6bIC8jKYNLGSth8h2d7e5BafVLZQbPiRnSPmH0gHzuvBQwpRzUtYvtB/qmPmoRPs7LMkdHoat+FA+Hbo/OalN5jc6d859Z/S2ulpeZ9g65xorhnowwp8ipIhZWmA8BUKI7DQ9P0RqWmPTNPzUJC1pOcbK15AzSYsczFpmSd6pSEufP4fV1arOTSQZIrWmnRh9T/NNtEGcMT2dW+fyt5cSGTCYbmFhYTDEcxjHbExOM7IfHkb4U4aU5RMiQB6GF3I8arJMDnmtrVVz6BDJnj49GP5HseJtPAfNOuYjibt39dh+raMJPdNcuWvYURJ/zm2THI+s4mVTeGroWErWkuTPO9u6kM85NP9i2OfR1MiYFRjhTyFSpKA5GrX88rwBaToyUE6tT5EXb+RkWceckXXlkBC0eHoiIp5CQLPk6z7TVN2GlSFC32+zg5T3yA0EGb4aisjh75T2ed2FUmSHtLGRZ6DUhTmCjfCnDm1Y+JLUcmWZuuRVZ0LRMI1RJk2T10tF/rSpFze10OVoRcopwy4vGLvH2LOPWfjynZJZVFOzoeUIIhX9Mww5DzsCmxUY4U8R6mr4skHKIbiWdkCSZRPyyqmr1GrrnCuvI//yOueQTltEMGzHwX8D+j3bcNzm3GOKjOldovo88ID/e/KkPmt3Z6daXtNOJfbbD/Mchu3YpxVG+B2iC72wTpSO1OQ5kfAIFKm9c4ttmMZRx2rPaYx1JpLVqXMbQ/22Oo4uHLfODedQpugbntvozJky8+Xjj1frS76TtbXqtVIdjkb4bRNyFz6SaYIRfkfoUi/M0ZdD8sbOTjVkcWXF71MaXbnJFAs5kTt16iqPxay71OimDadp7n1paOs374qUmtwjlwhluoaNjcFc+nwkqclpGplrMmRbur28l5QcWvf5TBuM8DtAG9Ze3ZdPkmNoYlMqhTE1ai3FAll6XXViKcs2dc4kOOWGJY1JlB1CdeKdrEb6BweD5ciOjJddN2dRXWMi1iZTKTZyrzPpMMLvCMM03LrEpZ0fImeeGjhE9nQ+b2yxTmRYQqrTQaas32lunG3JQm3VhSM0gY/ePRmxI53MsfYQC3sN3XOTzj30ndhC8nwSY+51JhlG+B2iydC8bqOPnR+SX0J6qXyx5Xe7tD5zGvAkWr9tYxJGKbIOWtQWP66F/IbknJiGzxH6LFReKGNqkzJjstK4OuC2YITfEYYhp7rfDWmg2vk5emlqeFu3E8tFzDqfJOu3a4xzlCKfqxayy49zspchvXwyXs7kNY4mDvph5ceU5DQLRoYRfgfItUBSZeQSK01ckS+/fLl5Yw2lDYgNb7U46VG+/JNg/c4DNJKLzWfgM62dK0l/ba2ZXJPbucs2Ip3K/LOUcRAj9i6NnFHDCL8j8BedzyTkcccxZ1SuVaFZ7JqVT/Xhk5Q4wfMGOMnD23Fav/MESXKp+QzyOFn2TUdloTYQi/oJva8agcu6ylHNpBg5bcMIv0PIIWHK2qh7Pn0nRfjaMF3ua3WXVo1Z2POBNmWMVFkpGY9/T5tRTbOQZRsIWebyndUMIdrXEuVNu4xohD8C5DSgYUYElLuGly+zQtbV+WPDW3muYXbQha8kJInUnX3LZx3zWclaiHHoMz5K5TPTQ4bQrBk5RvgjQkwHHFbzj+Wsoe/m6vzz5Bw16GiT5EKGhjZRTvqTtHdQy34amhhGhk8sPw+34kPGWFMjZxKNIyP8ESDHwm86jJZREjs7+qLUqZm0odEGr9u0WjWG+miDrHq9wQl8RL6rq/pscDqPjBiK56d3cHd30HBKzUMJOXe5dKpJP8M8h0ltQ0b4HaOOxVw3GoDrj9LyWVgoVyMC9GyZ0tqSEpC8D4OhDlKEn0PGvI0cHNSTGns9vVOR15PGEI0atJXQckh7kkfJRvgjQBeTiuSQ+O7d8Esd2siCavtltM7CQIj5jmIBB/J9bJJJNPYd2vhqXTIiR5Zfp53ktudRtxUj/BEhFY2QWoSCgzqQlPVSZ2uL7Cd1KGsYH6QVTyGb0trXpBX6PtfwtSgd7ZryO1KrB7zsSQkE+echyalOO0mN2MfRVozwJwB7e9WcHdyZlHKqyhwm2rJ+0qrn0Qmhl7EJJnkoaxgPQtE2nIB52uUQuYYsZi2kOHZtwLm3vGVwRrBsR/y6TdpJysIfV1sxwh8BUtY9dx7x/fX16gtNnUHoRebD0KWlqja5ulpdsLqrySRNnc+G2YMkMU1SkRE7MQ2/CfnK74QWa0lFEtV5l3PJfBxtxQi/A8iIl1S2Pe2H54tOcKt/ZcV/V2Yv5A4tslhkPD9ZUdKR1tSyCHVkTa0iw+xBS8SmvRtNcuek3lftO5LYDw7CjuVTp9KknXvfIblm1G3FCL9lyNwhUqfMjdCR+UBSscVc0+fLy1GkAt9vQzsMlaGlXzYLf76Ra9Fqo1n6v25aD+0czUHMo3Qk4Z85M1xq5BwH7cxZ+ADeAeAFAC8C+JBy/OcAfAPA7f7286kyJ5Xwc18yjew1SyQUwSDXF6W4+81N537pl/T4Zc2yiO3XvU8+uhi1LmmYXIQ0ay2Cht7TVFw9Ly9GvrnGV2hWLu905P2k7jd1/kxq+ACOAngJwCMA7gXwPIA3i3N+DsBv1Sl3UgnfubxhZIjs5Q8fskikBkl/eSOS0QwyjcKwL1TIOomlVjbMFyRx88lUMsqGv6cha77JO5wrr7YhrdQdOc9clA6AxwA8w/Y/DODD4pypIPw6PX1MnuEvcki/lFo736QlojmWpGOsq1V7Qo1kmNGDYTaQsmBjDtEupQ7t3Wzjek0t9lG3la4J/70AnmL7PyPJvU/4XwfwBQB/BOAHAmW9H8A+gP2HHnqo26ciUKcn1l4eoIz1lflCuFOKriXlnMVFnfx545GdgNbhtDl0HIf+aJgupN6RmFWtHasjl+ScJ+s4bPuYhjYxCYT/AID7+v//AoA/SZU7Sgu/zguhHSOi5pn9tLSr8vukWZJev75e/i81SC1xmiT8NsMwx6U/GqYPsVFgHQs/d5RaVybp9eqvxtXkficFY5d0xPlHAbyWKnfUkk6dnlt74fhCI6kyYteiHN1c+glNHT95skr0Ugoa9kW0GbWGFELvspYpM6bhy4i1kIFR1xCJEX2T9mEWPnAPgK8AeBNz2v6IOOeN7P+fAPC5VLnj0vBzCTM0pMwtIzXU5ZBO0t3dqgREEQgUU9zmi2havSGEFPnGnPshoymVxlheNyQjpepHk654mal3e1pGvaMIy3wXgL/oR+t8pP/ZLwO40P//VwF8ud8ZfArAD6fKnBQLn3KC8PNiJD6shR+rHwfVS3PmataVwdAFciZTccQMG2pbTQwmIm/t+pp0RP42umYozUnd+w3VNbbfNmziVQKxnptkk9SLkdv7t20l5L74BkNXaIvQhjGYQoaO1olw6VP6vkIrxDW935kLy+xqG7WFr4VN8pwccubr+npcL3RueKdT6sXKHdoaDJOOpgaTlruHztei4QBv3Ut/Vy7Zd3FPbcMIPxOabKJNjCLyD4VsxvZzz6uTd2SS9USDIRdNDSYtdw+NxkOWvEb4XbSZcTh5jfCHgHyZ6OXp8ofLJXOLojHMGuoaTDFC7fV8rhzZdhcXnTt9Wm/XXbXnrjsWDiP8huj1dEsg9HKEnFCxc2LXztUzm5RvMEw7UoaRNkLX8unkaPhdt+M2ESP8e2BQ4Rzw6KPAc8/5/fV14Mknq8c5Ll8G7twBrl4FisIff+wxf+yzn/Wf9XrA9jZw/Lg/3zn/uYai8GVdv15+RmXL82L7BsOsoih8W9rcLNvG1av+2PHjelt48klgZcX//+ijwLVr5bGbN4ETJwa/p7Xtra2yHYdA512/XtaR9gG9PXeNuSR8SbQh4qXPFhcHj1NHQN+/c6f6Q1665F8gwP9//Djw9NPA7dv+x5fkr9Vxa6v62dbWeF4Sg2FSIQ0nTvpbW8CNG8DGhv9LePRR4MoV4MiR8nvUrmTb0to2J/GU0Va3Q+ocIdN/3FtXkk7dnDlyCMijdWjYyMuJDRdj4WPyuuaQNRiGQyhKh+fApy02+zbUtvm5chnG1HyDLgHT8D2aEKnmcKHvra3lRQtoPoAUeZtD1mAYDjxSh9oSN+Bo8hV9Rv/zNsbbX4jweUpoOnecbdUIn6GOEyV2rrTSQ2t5hhY4yfHWm0PWYBgOoXz9585V2+LKyuCqWNT+tRxZdJ5cKH0SZrnPPeFr8fUp4tWsA7mvdQj0o0vPf2yCSBf3aJ2DweChhXDGDDGZvE2u8KZ9t8u2XRdzTfjaJI3Qj8PDKHd3fQoFmlF7eOj/52kVNMmGa/qU/pj3+nL1nzZeDJN/DPOKuoZOiLA52Utu0JLAyfNyjMhRIUb4Mx2l41zVw37lCnD2rI+UWVoCbt3ykTLXrwOf+5w//9FHqxE1X/2q3/+d3wFeecWHZ+7t+XN5JA5he7v0yH/2s/4z7q2n6ADn/LFhvfXyHutEERgM0wwZLlkn7DkEapcEHooZK+vs2er+xEbUhXqCcW9tWfgx2YWO81XtteGZNtzLzd8t6xLbb/MeLZrHMMuQWvnOTnjBdA4aeUs+WF+vckCoHcl2vrPj3MJCWc7BwaAENGpgniSd0GxXTXaR58SGetpwr6t1ZJtg1NO3DYZxgSTM0FrPBwdVUj448N+TGj45cGmfonXkanOSuLXrLyzoC7ePA3ND+MMsqkDn1yH8YVIntAmz8A3zAi19Qqy9kuW/u+u/u7vrFwpaWSl9c8QRa2vVLLjEH7EkidqIQQvvHCVihH9kzIpSa3BMy97aKjX2Gze8zr6x4TW+zc3qOfz7ly6Fy19crM7WAwZnwgKj1+ycq2r2sXs0GKYdRQEcO+Z9cNevA0ePhs91DvjYx7wv7umnfbv4+MeBl1/2M+UvXSrb682bwLe+VX5PlqPVg7C6Onj8ySc9H/V6+fc2EoR6gnFvTSx8zdLNkV2knCM1fJ4XPzbca2pRDztKsCgdw7xAa+O0kZYuN95+AecuXhwcrfM5MzK0mkYIUgqm+mhSMIVij6MdYl6idLSEYzxShp/De+ii8EmTVlf9dvw4cP68750//3ngx34MuP/+sjxuZRw7Vi+hkkTTxEyyDOfi92gwzAKKwke6Pfust9wJCwvAq6/6NvnudwO/8ivlseefr5Zx5IhvXxzXrvk29JnP+HJp5LC05PfPnwfe+lb/vU9/2h+jtqrl6wHKYxMVLRfqCca9tWXh19GyuSYv4/Ll/7FJWbG8GqH6Ws4cgyGNmIW/tOTc3bvh4/JcjSd2dvQRwsmTVeufEMrXM05fGubBaVuHPNuQULQfeXVVn6QRG9aZw9VgqAceXcMJfGenGoUTWsuCyy4yckfKP/IamjEnjT3+HU0G6hoxwp8Zp61MRQr4v5ub1eHb5ctVZ6Zz3nnD5RM6pu07B3zzm4PDN8Dn2ZaOY3Ko3rkTdv5QfQkmxxgMOpwDXnutmmZ8c9Pv/93feYl1Y8OfS7nvFxerZayseOn22rWy/W1seMft8897eUjDrVteDuIgOZgmVW5vV49vb+cFTsQ4p1WEeoJxb03j8A8Pq07Mw0M9CRIdp5AqLWGSc9WYW+f0iVqaoybXYjcL32Coh1D4NU95Qm2dUqJo7ZSDvhObdLm4GLbYtWSKufJs20EXmAdJx7lBqYX+J1IPTdbQhndyAofMgaNlxgwN60zDNxjahZwIpUXicYknl4R3dgZlnZiGz9GEuLvggLkgfKnFaWFX/IfQyFo6dGJpj6WFT9fUUqny5GwSFlJpMDRHijBzfWrSwj9zxhM9kf/6uj+2tpauT2w/dQ9tjPLngvCdi3vw5Q8eI/zQKECeL0cFfDUsiv+Xo43YrL3QvsFgCCNFmLH21esNKgOUW4dy5sfi8Nu8hxxVIAdzQ/jOhdMjhIZ0WjhVasp2bDJXSlaSS6MZuRsMw6MJYcq2SlZ8Su/vou5m4Te08GMz6OQQj58vZ9BqYV+8w4gRN3f+5nYUBoOhGZoQZkoCbsPSrlt30/BrgD+olJwSc/hwzZ46iJy0qxrkKIGiBcxBazC0g2EIMyYBD2Np15VoRxmlU/jjk4fl5WW3v79f6zs8TQGB0hTs7YVj252rHtvb87G+wyyssLdXLqJCWFoC7r3XJ24i0LwBi7s3GJphmPQkzg3G1gNlLD+lRqCFi/j3tDbbtC6yvByOCaEoilvOueXAhcZvzWtbm/nwuygntE+x/3JhY74/yiGjwTAPaBodo0Xo0UiclIG1tfxIn9RoYxQBGpiHmbYE2SvGrPrYfqwcbbbu1hbw9rf7v0UBXLhQJl46erRcVvHee6vlWgpjg2F45LZ7ArXZGzf8rNuNDb/dvu2/u7HhZ9BevVqmYk7NnqdZu5Se/MiRMm351avAE0/ovJE7y78VhHqCOhuAdwB4AcCLAD6kHL8PwEf7x28CeDhVZltLHGoYRjOL9eJyaTM5w880fINhcsB5gAda7O0NBmHkOIVD83xkEEeo/bel5aNLpy2AowBeAvAIgHsBPA/gzeKcfwfgP/b//ykAH02V2xXhD+sVD0XgxGbxWpSOwTCZkCGXoRDMVNhnbNlF3qnkHBvWIOya8B8D8Azb/zCAD4tzngHwWP//ewC8CniHcWjr0sJvGvcqLQLtBdBm8coQMIvDNxiaQ7aZ3HTkErkWdc7ELjnKl39TvJFznVzECL8NDf8MgL9i+1/rf6ae45w7APAagAdkQUVRvL8oiv2iKPa/8Y1vtFA1HU0yVDpXZsK8dGlwOcRLl3w0j7bsIZW/uel1Qbn4isFgyAP3nznno+HOnvV/6TOpi2vg7TmmzfPPQ0uIcu2eovIom+etW9WMvZIfuKY/ksy5oZ4gdwPwXgBPsf2fAfBb4pwvAfh+tv8SgIVYueOw8FOWQmpil+zN5SIpZs0bDM2hzbWhXDeLi4OJEnPk2RyLem2tmjSNou5kXp3YqD83gqdrC3/uJJ2U05WnQQ4N70IOmdzwLYPB0Aya0SWzW0rZNFVeTJvP1dZzyDomIU2Thn8PgK8AeBNKp+2PiHM+iKrT9r+myh1llI6cXRuyzlNOF+cGLXmz7A2GdiFJWm51yD53pJ+r4afIOhaHPxVROr58vAvAX/Slmo/0P/tlABf6/98P4A/hwzKfA/BIqswuCd+5wQfPc99LS4E6ADlc1ByxodzcZuUbDMOj1yvbYWirK+fIPFmUHVOeH+tU2iLrWIeQi84Jv4uta8Ln4CFV/EddWdF1Q/6jrqyUidd6vTLWXiZjs1h7g2E4SDlHS3Ao57rEcP68cwsL5fkHB9X9OnINnRfbHxWM8CPQQqpCVgMtl0jQEqtpZRjZGwztYG9PT2N86pSrBFGkLOtYu9d8eW1o66OCEX4CGlHTkma0ySyZqRemiaZoMBjS6PW8sRVqn8No+FQWx7StShcj/JnLpZMD56r7lPuG49Sp6v7Zs9UYXS32lvD449V9y5djMLSL116r7m9v+4yWly/nx60XBfCbvzn4+fZ22V6d8zH+Mh7+6tV0rP9EItQTjHvrysIPrXivOYJkuJe2gHFoZSxaKWfSh38GwzShTYlld7fU7Gm7//6yrMNDzwvkj+PXn1Tr3jmz8L8D5wZn121u+ox5N2/6DHmHh8CZ/jzhtbXq97/v+6o5sXs9b/lznDxZ/r+1BRw7Vs60s1m1BsNwKArflvg6EjTSrtPGej2/XsWrrwILC8Ddu/7vt74F3H8/8L3f69vvzZt+u3TJcwOfjdvrdXefnSHUE4x7G8XEKx6NQxb5+rpzp087Vcc/dUp35vDoABoZrKxULQWDwdAe2oiI4UEXtEmLX1v+kCJ6JtXSxzxb+E7RzmW+ikcf9Z+97W3Ak0/6nh4AXnmlet7LL1fzZ5D1fuuW7/VplZznn/erWpEVoq2oYzAYmoNWk9L2ZZvXOADw+elv3ap+9vLLg+cdO1bd/43fAJaX9Zz4k457xl2BLqEtN3bpkh+icdy44TfC2pof1vGlCNfXfRmU/IyWQaSlz7ThnS1daDB0g9BSgrTQUM4Sg855By3HslgYkPMC4bu+y/89dcq3f7rOVLT1kOk/7m3YJQ41Zw4fmtEEKi0ZWswJ61yz+Htz2BoM7UCb5U77FGiRk+KAH9OWJQ1xAW3r69WUK5Mi72Be4vBjOXJo47NlNcKXq1TFIm80kr94Ue8QJumFMBimHVqKBdLb5eeh6B3JF9Rud3d1btCi+CYxEm8uCD8WrsU3Inr+QxGRaz/wxsZguoRQ2RcvVl84enlyX4hJmZptMEw6YqNz+bmWwDA08qbZ9NRm19edu3t30Jkr9yeF7J2LE37hj08elpeX3f7+fq3vOFeGTRFoIXHC5qbX3ba3S73v2LHS+fLkk8Dp08Ab3gDcd5/X8Tc2/LHXv77UAQ8PgXuEB4SuRc5brv/RdbkD1zHdL6RJatqjwTDvcH1/nNTYtba3uuq3a9fKz1Jt6/Jl4OMfB1ZWfHvUtHyOXm9yNPyiKG4555bVg6GeYNzbMBq+tLxjy47xLJmrq75H55YC5eXgvXdIswd86Ob6+uCErPPnR5ML22CYdUgNX8os1G5k2CW17dQiKVLupfDsBx/0f9/ylsF2P0ntFPMg6TinSy1LS3rq49xseNpQkGv4BwfV7507V+1U5BAwROi5GfkMBkOZhlwS/pkzuq+N5sVw+Sel64d0fOpUKGGbXEtj3JgLwk9ZydLi1og8dpyDljy7e1fPoa+lZNA6AS1yILcOBsM8g5OxjMKLBWXItiXbX6hM6SfgbXjSgjJihD8zcfihKdeA1+hlvO3W1qBeHjrO4ZzX6q9fB86dK/0At2/7uNwbN7zux/Hud/vY/mPHqv4EiuGlcuXC6KE6GAzzDpoTw9v7tWslDxRF3PdFaU9ee22wja2u+vbNfYEhcK6ZinYa6gnGvQ0bh0/Qwqak5V9XPw+tjkXRPDIlQ2g0QPG+3JKwhVMMhnzEom9C1r3my+NhnvJ7tBCSVt4ktk/Mg4VPkL3skSNhy//48fL4xkb1uHN6Mqa3v91bBfv7VQvg9m3gU58CPvMZn1qBQNY/jQYoWufsWb9/9qwvi2b/rq76a/M6ToXlYDCMAbJt5LSVxx8Hzp/3bWttbdCaf/bZ6vnf/nZ1Jv3qKvDZz1YjAqfFwp85wtdw+bInUfpBmg7Bej1P9iTfcPzZnwG7u4O5OB5/3Id33n9/tdO5daskfQrvpE5ne9u/jNPyEhkMkwaSfDY2qiGVGxv+8729Us7lZL+4qBtsv/iLpWFI0tFUGmUh03/c2ziWOMyRdPi6l7RRDm35OW18CjZdk8ri503irD2DYVoRcrrGovPkLF1tpTt5jUkD5iFKZxjUDYmUy6ulNv6S8QXTtdm6RvYGQxoyZFrbr+u7C2n005bePEb4lrgX1eEZISSn0MIJIdBMP75PGTadKxdgoSUTl5YGy+DROwaDoYrLl4HHHvNRbUTNly75zygyJ7VQivTtcayu+na+uenbKF/ycNoxFxp+Cs7lhWXSClcUirm/77X8V18tz5EOH8DrhUC1YyHdkId1Era3Tb83GJyrtgEi3W9+s1yJikA6PQU9AGnfnTxOmv9Ua/QphEz/cW+TquHTpCsemrmw4FMnyPSqsdSsoaGjTPtKkMNKk3wMswxt7WmeikSTX4ZNV5ySiaYFMEknjLprZH760z7C5ujR8nsvv+wt+wsXvLV+4YIfMmrlaKOJpaUysdrVq95K4dbL3p4fWdBIgcqwpGqGWQSXPmmFOb6WLFBNhEbgkyjrrkZ1+XK1XVI5Tzwx1K1MHkI9wbi3UTptnWvufdcmeoWOx6Zua/k7QgszWCSPYZQYR2RKKr9VKld9nbYxa8kLYVE63aBuQ9jb0xdgWV2tDj+1CJ5UHh6DoQvEpJWuoeWW0shey1xLhlduG52l5IUxwp97SacpaAhIQ0YaAoZkFudKZxOB8nmvrpZyzeXL3ml75Ur1++9+d3U/x6lLdQvtGwwxuIS00uX7RNfioP0TJ3ybWV/3TlYe8EDY3vYTIXkb7fXCbbROpN5UI9QTjHubZAu/yRAwlM6VcufIckPplXOtj3FaZobZwTgs35z2RZvMgstXrqM2Q+mSY2mM58XCHzuxh7ZJJnzn6r0g/NzUzD9tcRV6cXM1/FnTJA3jxTjSducYLKGZsuvrZYpyLZGh5mebpfZihN8R6jSEnFAyWl1LzuTliytTWSlrfZYsFsP4MKr3SLO4ZSoS/lcaUZSpljZtVSqu7UvM0oi4M8IH8HoAfwzgL/t/TwTOOwRwu789nVP2pBN+3YYQiyzY2ChX5OFLLPJyDw4Gywtdh/8/asvMMDsYleUbI9vcYzmRO1RvSknO71OLuZ/W9tIl4f86gA/1//8QgF8LnPd/65Y9yYTfpCH0eoOJmbTogtj6uzlRQLEGYBa+oS5ooiFZxiQ5rq21U36sLXGfV0zLp3KkgRNqW7TudCxSbprRJeG/AOCN/f/fCOCFwHkzRfjOpYeAofh70uw1a2SY5dNCDSXUcAyGFEZl4adi7lMj6VBAhPb+8+AHrZ3MQtvokvDvsP8Lvi/OOwCwD+BzAH48Ut77++ftP/TQQ10/l6ERGgJqncHqahmRs7fnLQxp4cuFlqV2mVMf2Th4zH9uxxHbN8wXmkiXsf3Q57E1p7VIHF4OJ23Zpsh5S+S+sxM2tmblXR+K8AF8AsCXlO09kuABfDNQxpn+30cA/G8A/yB13Um38ENIDVF52JhmZQwrwcR0+zqSEL+XWRnqGpoh1xeU+/7I87TINDpOAQtSmqEABrqOlEvX1/22u1vNwaPdz6z5t8Yu6Yjv/B6A96bOm1bCdy5uaTvnX3BpzdNLPcyatsNEVMxaaJqhHeS+U7nvj/w8ljpkY6O00OmYPJdfXxK4NkKOyamhEfW0vftdEv6/F07bX1fOOQHgvv7/C/2Injenyp5mwndu8AWkFyqkN54/X44A6Pt1rOs2CHuYDsMwe6j7TjXpHLjBo737uRZ+7nVjo2uSfKZ9lNsl4T8A4JN9Ev8EgNf3P18G8FT//7cB+CKA5/t/35dT9jQTvvYC0jAzJN8sLTl3927Y8ZuDNiSZ3OG7YT5Q953KfX+kLh8LO05p+HU6JZ7Pane3JP3VVU/2MipuGke5NvFqhNBeklR8MEUOnDo1/MslJ5bUWZ7NLHyDhjqO2Jz3J8dqr1Nmk05JSkd8f9oTFRrhjxjaC5gzKURqiW1cN/Xi8/9TMc8GQwi5lramw8d0+VzrvYnuHupMYhFD04AY4dsShx1ALp1WB3I929xynCszGwI+0x9l49zc9MeB6vJud+5UMwLevOmzEPKFYIAZW+LN0AlCCwkB1ffnyBG/QBDgs1wePer/5wsHceSUSdeX9cmp89WrZZsBfPvb3q6epy13OrUI9QTj3qbZwueQYZkpSz/kvKp7LTli4GWlQkel5TTtUQuG0SH3XZFWNMXHy1ExD6dMldmkrprzeNpHuTBJZ7zgUsvurnfeyuROy8vVoS1fM5eGmTkIRQdpU9Jz9FaLzTe0De3dI/+VnKvC04d3VQdJ7sMYXZMAI/wJAA/JlJn9aPvRHy1ffKkp5s6SjY0gNN0zplVabP58I8eyrmt9y3dod7eMtW/bp5VCyJiRzuNpe8+N8CcA/EWXObr5C//BD1aPyWgCGa7Gy6fZhnSe7Fgk2ccsfNkIpjlqYZoxLjktZ2TXdPRHhg9/t9/ylsHJiOvr8TWi28AsypUxwrclDkcEchBtbgKvvFI99vzz/u/KSunEIhw96p1KZ84Azz5bOl97PeDsWeDtb9ev5xzwp39a/YyWe3OuXKpuc9OXtblZLmW3t1ddGk4utzgzDqwJR91lNNuCc+mlDXPOCZVNKAr/bm1sAF/4AvDcc9Vzn3uu+/tv4uydaoR6gnFvs2bhE6SMIoevodAw0vxjoWyapEPHZchlbDhLVhdNUJFDbrkso7w/w/AYt5yWM7KrO/qTs8zpf2nZh2QdkxPzAJN0JgNaA5FEKoe6/KWXxBtark1KQtQRyJzfIbLWOo7FxcFp6NQ5mEN3OMR+h3HKaSkfT+45dJ4k+hjBkxwp5U9NduTXMBjhTwRi4Zl8n6df4A5bPgLgZC6voa2ZK3Px59Y31Bg52Y/LAp0V5KyrECLULgmvDQufn5vKWb+yUn33ifTlcp9UpkWPhWGEP0LEGiF/SXkD4LHxoWXdNjbiCzJLp21ssfScewiNRHhc/rgt0GlHSrbRRnopSa4NwsuRk7RoGx7DziPL+LkyjTEf3dJolL4Xun/+uRkbgzDCHxFyGqG0kEIWm4zGSaWJlYQfSzglIa8bG3prE7NyhvQGHaFOM0ZqOUv/NakHx/nz1XeLL20orWyZh0Yb/WmjT6A6H4VPEox1OiHZ08jewwh/BOAvNL14bS4rmLO2aKhRxeqgdVKrq36ILYmfOpS6k7ckupIiplXT1fwuzpVWM8WFcwOiybMPPZ/YgiTkV5LvgszomiPvhN5N8g/J0Uno3aSgATM2dBjhjwAhjbLNmYI5mTA1yz5k/aWsKO6U5QSfSs+QitToQoqYVk1XI0q5EEhIF69DeDlRWfJ3lEECqbWRU/6GkJyzuxuu+zCjz3mFEX7HkETZ1ktYx2KNWdyxmbqpoXGuT4KOxUi2aSeRQlfldg1Ns89NzRv7rTXZsImvQBoPsTrFLHwZeCAnBK6v5/9G2kgh1RHNG4zwR4DQkLUp4dchUy2+eX29KsFoowGuwUqLKxd1ZZSmMlBOPaZR09XklJTVzu9VSi5LS/63l4EAOZJLaGQoN22kGaoP/b+25utFMfcUkUP7dUbC2khiWkZ0o4AR/ggQsjyaEE8dizXW+KXeGrqGthgFWYqy3m1Mde9Ke51WTTdXB+eQnfz6evV31Kze0PPRyJ3K0lZp05KLra1V/TwrK/6806f9/unTXiI6c6Y8zo2SXMKPPaNp+b27hhH+CKBplDHCzSkvV2rhM2nrNAJNQqB9sr5WV8vl52JOxGHvKzcbaN1yc0kktj8KNJGlQjJe6D2Q521s+N+W/+b8vAcfrB67eNG5172uSvq8jgcHYcfs/ff7v7S6W5NR8LRKd6OGEX7HkBq+tt+0XM0i0+L5pS6aK8vIiS2nT/tGubxcfrawUIbphRp73QYrU0AvLfmc6PL8nP1hiGCSnL1N6yLfk5DcwTt0ysrKf1MuBXFrXabnkKSdG40T2uqMGifp95pUGOGPAG2/iCkHrBy2y01LuyDLT1mHFy8ONu6m633KUENO+nyfOir5PKlzoOe5s1M9TuVSp0H3Fuv42rYY25C/2vCJaL8P99do0iMZDHxEEHoHDw50Q4TqU5fw6+afn4QR2STDCH9EaOtFlNZzbF/biKRjztomSawAvbHH7rPXq84h6PWc+/a3yyG+7Ei0+z1/vnpPd+9W96UFy+cqxMgjJnVIEktBdlBtyF8pyJFkKF2HlPdCI4JQ9I20xmMaeh2y536HpqNGwyCM8KcIRBwUH02NgPb5hJeU1URWMCcZGVHBnXLSete27/7uQXIM+SmoY6F4bppgI1NEaISjEfHCQnVEEBp95FjrnKDls+SEn2txyg5ZzozuKmxQk/e0dB1aXeWmkX3Imas948PDqgzEf6OjR/1f8gHQ85HOZiP74WGEPyXQGpDcl+eFLPy7dwcJjw/rQw354sXBBkiNVetUYlYxtz7lJJ4HHhgsT5NAYh2aHG3kTrePWcbyfrROQytP+01iaaW13z62H4N87qHRSeq+tY5OI/eYBEOdDhkoOzv+OTz+eDni4cf39vJCUQ35MMKfIoR0U60hxjR8+V2tM+Hn7ewMDrHPnSststgWIrJUxxQKB835Lnf48nvNJY+Ulp2qFye42EghRKgco3REauGc3EHL6xiqF/eLaB0M35yr5uOhc/i5OZ20IR9G+FMGSRwxOYI7QbUtZOFpHcPh4aDezvXynGuk7oU2GjWQ1RfT8A8OButAIwSZ5ydHHohNOAvp0CEJg3e+obWKY4TftuM4BzTa4+sv0H3nrpngXLhDWFvL68DGce/zACP8KUKuxSOtI2q4KcKTJKc5yjQrO9aphBpniDzvu69K+kTqmuXMNfsHHqjOdXjgAX3iUYw8YvdEMpGss5ZThpenWfI8CyS39lOkPyorl9ebrlXHz5CSfHJJnEYc8nlamGVzGOFPCZpYPJouH0tRG5JBQhICt76lFRvLYcKvdeaMrtnLumlT9p0rRx0HBzrB8tFOjmUZm3AmyUpL+CUJSXaiBweDjmme7kJ7zqlRHT8vtl8HdYwLbT/1PuV0YEb47cMIf4rQRM+VmqumN8fkCP4dDZxI9/ZK3TcUCaLdi3SwcmIjYpVr5fJyef1i5JhLiHLCGT0PLkeECB+opugNyWQkQ4U6E7qvXOKt+27kPIs68qF2zdD3czowk3S6gRH+lGEYK66J5ppjTUlSjZGs/F6MFLX48VjDj5Fj7nOLWab8ezy8Um6Li1XHMUlq/BweKcVHEnISGs1/kP6IkBTVllSS6mhS1wxZ8nUs/HHIWbMOI3zDd9CmJJB7PantaiSqRRzFyF4SEDkgU51ZXeKksEJJ9nxGsExax0cBshPgpMjTHKyv+5BY/iy0VBO5EoyMuqG6aSOTXJ9HiNSH0fDpGqnRgCEfRviGsYI0eJoIxh2xUtpJNXxtlKJZxSnrN2ekI0MYpZzj3GB4pzxfdhYh3wFtlBuenhFf0YzqmiOVhMJ1c8IuNZ8HL4NGKMNG6fBjZuG3h84IH8BPAvgygB6A5ch57wDwAoAXAXwop2wj/NmAJgOEJjrJsMaYRchRRwcPlaHtx+Y4kISzuemtfS1TKpG97IQkecoOgFv/qbkJsWeUMw8gtt/r6WG/MmWE9v1c/4Fp+O2jS8L/xwD+EYBPhwgfwFEALwF4BMC9AJ4H8OZU2Ub4s4MQ2Z08OZhv/eRJPcIohNjEJ67rc4tc1i2Ew8NwGgi+5i91ZPLaMmooFh2kbadO6WSfQ465hB97rnIi3vp6NctmzNGfi1FOOpsXdC7pJAj/MQDPsP0PA/hwqkwj/NmCJGMep76+XiVWbj2nopO401QS3LlzZVkbG2WK5zoJzWQ0jzYKiTkpQ4TNyZ+ieWTsvlz6L5cccyWdnOdKC5mELPw2MGq/0qxj3IT/XgBPsf2fAfBbgXPfD2AfwP5DDz3U8WMxjAqahb+4WDopQ4SU0/A1ciMrlF8LiC/eoV0zZCVLSYZ3OjFHssygSVlAKa2FfBaaFZ0rlaSctjnPNRRZ1YZlb+gOQxE+gE8A+JKyvYed0wrh880s/NlAzLLV5IwmZCJHD7EUB9KJysleykOSJEO+B0oIxucmpBzJvKPQMkZqGn4dtDGhST5X7ZkZJg/jtvBN0plzSClC07Kb6MzO6ZZoLK+NprXLcqiu3Ermzmau30s5iXcYKSerJgXxTkuL0qmDYaSS0OimqzTPhvYwbsK/B8BXALyJOW1/JFWmEf5sQUolUs6RKRvqyg6hWHBNJolJSDFHJ7eatXNDzlOtg+G5gvhxCl2lZ9I2qdaNnpFzDPhIxjCZiBH+PRgCRVH8BIAnAbwBwP8oiuK2c+7HiqI43Zdx3uWcOyiK4iKAZ+Ajdn7XOfflYa5rmD4URXX/858fPL6+Djz3HHDixOD5oTKPHwc2N4GrV/3+lSvAs88Ct28Pnr+w4D9fWgJu3QLe9jbgxg1/7No1/RrXrpV1uXzZ0zLtX7tWfh8o60BwDtjaqpa3teXreOcOcP26ryvHxz4G7O/7/48fz3sOubh82V+X6kn1O368em/0XDc2qs+FPt/ba7dehhEi1BOMezMLfzYhI2uk07Ophs//p1QFJL3Q6OHMmXKWbK5GzyWa0H1oFn4qjFKbfCZTRLeJVH20xHNaKgbD5ANdWfgGQ11oVvnVq/7Y8ePAkSP1yuMWN+Gd7wQefXTQOj1xAtjdLa9Bx2/cqF53ZcV/H/BWuFcmff2eeKK0jK9fL++D9oHyvmL3efQocOFCdSRy6xawvd3sOaTAr3/9elnXzU0/4tjertZ/a8s/l83NquVvmHKEeoJxb2bhzzbaiL2OxaXnlq9F+Ei9mi+uIhOfpWL6YzNRQzlqukTIp5AasRimB7BcOoZZQxvT8kMRPlr0TO7iKqHrcMQSj3VJsjkylNYZGKYLRviGmcQwVmmqw4itDFbnWm0kGGsDOfdrFv5swAjfMLMYxioNkTHP5SPlljrXyiFZeX6XyLlfS2I2/TDCN8wk2tCd68gtOQukd1HHNhHyKVgSs9mBEb5h5tCWhq/tx2YGN7lWnVFIG87sphjntQ3tIUb4FpZpmEqkwh5pYpGcCMUnUeVMQgJ8iOSFC8DaWvhaIVC5HFtbg5O0curUNWR9LAxzBhHqCca9mYVvyEETiaLp6KCuBVznOm2MWAwG58zCN8wwNKvUuTJ1AVCdGLW5WX4GDE5C0izv2LVSdUuNQvi5TepkMNRB4TuEycPy8rLbp6QiBkNNkBxCxAkMkqdz1RmtvV43xBqTlrRzR1Enw+yiKIpbzrll7VjLE7gNhskAt5gJkuw1bb0L+yd3ZDDKOhnmE0b4hplEjDy59b+56a3ozU2/Py6CncQ6GWYPpuEbZg6SPLXkZrna+qhQR+83GJrCNHzDTCInxLGOtj4qTGKdDNOFmIZvhG+YWRh5GuYR5rQ1zCVsIpHBUIURvsFgMMwJjPANBoNhTmCEbzAYDHMCI3yDwWCYE0xslE5RFN8A8NVx1yOABQCvjrsSY4Dd9/xhXu99mu/7B51zb9AOTCzhTzKKotgPhT3NMuy+5w/zeu+zet8m6RgMBsOcwAjfYDAY5gRG+M3wO+OuwJhg9z1/mNd7n8n7Ng3fYDAY5gRm4RsMBsOcwAjfYDAY5gRG+BkoiuIni6L4clEUvaIogqFaRVG8oyiKF4qieLEoig+Nso5doCiK1xdF8cdFUfxl/++JwHmHRVHc7m9Pj7qebSH1+xVFcV9RFB/tH79ZFMXDY6hm68i4758riuIb7Df++XHUs20URfG7RVH8TVEUXwocL4qiuNF/Ll8oiuKto65j2zDCz8OXAPxrAJ8JnVAUxVEAvw3gnQDeDOCni6J482iq1xk+BOCTzrkfAvDJ/r6G/+ecW+pvF0ZXvfaQ+fu9D8A3nXP/EMBVAL822lq2jxrv7UfZb/zUSCvZHX4PwDsix98J4If62/sB/IcR1KlTGOFnwDn35865FxKnrQB40Tn3FefctwH8FwDv6b52neI9AH6////vA/jx8VWlc+T8fvx5/BGAf1EUU590eRbf2yw45z4D4G8jp7wHwH92Hp8DcLwoijeOpnbdwAi/PZwB8Fds/2v9z6YZJ51zX+///zKAk4Hz7i+KYr8ois8VRfHjo6la68j5/b5zjnPuAMBrAB4YSe26Q+57+2/6ssYfFUXxA6Op2tgxc23a1rTtoyiKTwA4pRz6iHPuv426PqNC7L75jnPOFUURiuH9QefcXxdF8QiAPymK4ovOuZfarqthbPgYgD9wzv19URS/AD/K+edjrpOhAYzw+3DO/cshi/hrANzy+f7+ZxON2H0XRfFKURRvdM59vT+U/ZtAGX/d//uVoig+DeCfAJg2ws/5/eicrxVFcQ+AYwD+z2iq1xmS9+2c4/f4FIBfH0G9JgFT2aZjMEmnPXwewA8VRfGmoijuBfBTAKY2YqWPpwH8bP//nwUwMNIpiuJEURT39f9fAPBPAfzZyGrYHnJ+P/483gvgT9z0z1xM3rfQrS8A+PMR1m+ceBrAv+1H6zwK4DUmcU4nnHO2JTYAPwGv3/09gFcAPNP//DSA/8nOexeAv4C3bj8y7nq3cN8PwEfn/CWATwB4ff/zZQBP9f9/G4AvAni+//d94673EPc78PsB+GUAF/r/3w/gDwG8COA5AI+Mu84juu9fBfDl/m/8KQA/PO46t3TffwDg6wDu9tv3+wB8AMAH+scL+Aiml/rv9vK46zzsZqkVDAaDYU5gko7BYDDMCYzwDQaDYU5ghG8wGAxzAiN8g8FgmBMY4RsMBsOcwAjfYDAY5gRG+AaDwTAn+P9E+1ZQ79NoLQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "target_samples, y = sklearn.datasets.make_circles(n_samples=1000, noise=0.1, factor=0.2)\n",
    "target_samples = target_samples[y == 0]\n",
    "target_samples = torch.from_numpy(target_samples).float()\n",
    "\n",
    "plot_samples(target_samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([500, 2])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "target_samples.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Restricted Boltzmann Machine code\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The class below represents a restricted Boltzmann Machine.\n",
    "There are 3 functions you need to code:\n",
    "\n",
    "- forward: it should return the log-partition of the input, i.e. c(x)\n",
    "- p_x_given_z: it should return the mean and variance for each observed variable\n",
    "- p_z_given_x: it should return the bernoulli parameter for each latent variable\n",
    "\n",
    "Note that due to batching operations (i.e. each function takes as input a batch of data, not a single data point), computations that contains the matrix W will be tricky to implement!\n",
    "Each time you have a sum with W appearing inside (look at the PDF!) it should be implemented as a matrix multplication, i.e. with the symbol @.\n",
    "As an example, in the function p_x_given_z bellow you will have the following operation: \"z @ self.W.T\".\n",
    "You should take time to to understand how operation broadcasting work and how to implement batched operation in this RBM. You will have to play with the transpose operator...\n",
    "You should probably read this page on operation broadcasting: https://pytorch.org/docs/stable/notes/broadcasting.html\n",
    "\n",
    "**Advice:**\n",
    "\n",
    "- for each tensor, think about what each dimension contains\n",
    "- how does it compares to the way we wrote it mathematically in the PDF?\n",
    "\n",
    "**Sampling functions in Pytorch:**\n",
    "You will need the following sampling function from Pytorch:\n",
    "\n",
    "- torch.normal\n",
    "- torch.bernoulli\n",
    "\n",
    "**WARNING:** Read the documentation! torch.normal takes as argument the **standard deviation** and not the **variance**, so you should pass sigma_squared.sqrt()! Always be careful when you use functions from a library: take time to read the documentation.\n",
    "\n",
    "**Stable operation**:\n",
    "\n",
    "To compute the value log(1 + exp(...)), you should use the softplus functino which is more stable: https://pytorch.org/docs/stable/nn.functional.html?highlight=softplus#torch.nn.functional.softplus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import axis\n",
    "from regex import R\n",
    "\n",
    "\n",
    "class RBM(nn.Module):\n",
    "    def __init__(self, n_visible = 2, n_hidden=10):\n",
    "        super().__init__()\n",
    "        self.n_visible = n_visible\n",
    "        self.n_hidden = n_hidden\n",
    "        \n",
    "        # Create parameters:\n",
    "        # See the PDF file for a description of each parameter\n",
    "        \n",
    "        self.b = nn.Parameter(torch.empty(n_visible))\n",
    "        # because sigma squared should be greater than zero (variance),\n",
    "        # we instead take parameter log_sigma_squared which is unconstrained.\n",
    "        # as usual, to get sigma_squared you can do:\n",
    "        # sigma_squared = self.log_sigma_squared.exp()\n",
    "        self.log_sigma_squared = nn.Parameter(torch.empty(n_visible))\n",
    "        \n",
    "        self.W = nn.Parameter(torch.empty(n_visible, n_hidden))\n",
    "        self.d = nn.Parameter(torch.empty(n_hidden))\n",
    "        \n",
    "        self.init_params()\n",
    "    \n",
    "    # Initialize parameters:\n",
    "    # unfortunately, the RBM and training loop we are going to code\n",
    "    # is really sensitve to initialization...\n",
    "    # this is not perfect, but at least it works! :)\n",
    "    def init_params(self):\n",
    "        with torch.no_grad():\n",
    "            self.b.fill_(0.)\n",
    "            self.log_sigma_squared.fill_(math.log(0.1**2))\n",
    "            self.W.normal_(0., 0.01)\n",
    "            self.d.normal_(0., 0.1)\n",
    "\n",
    "    # returns the log-partition function with fixed x, i.e. c(x)\n",
    "    # - shape of input x: (batch size, n_visible)\n",
    "    # - shape of return tensor: (batch size,)\n",
    "    def forward(self, x):\n",
    "        part_1 = ( (x - self.b)**2 / (2*self.log_sigma_squared.exp()) ).sum(axis=1)\n",
    "        part_2 = self.d + ((x/self.log_sigma_squared.exp()) @ self.W)\n",
    "        log_partition = part_1 - ( F.softplus(part_2) ).sum(axis=1)\n",
    "        \n",
    "        return log_partition\n",
    "    # compute the parameters of the conditional distribution p(x | z)\n",
    "    # - shape of input z: (batch size, n hidden)\n",
    "    # - shapes of the two returned tensors:\n",
    "    #   1. mean of the Gaussians mu: (batch size, n_visible)\n",
    "    #   2. variance of the Gaussians sigma: (batch size, n_visible)\n",
    "    def p_x_given_z(self, z):\n",
    "        sigma_squared = self.log_sigma_squared.exp()\n",
    "        mu = self.b + (z @ self.W.T)\n",
    "        \n",
    "        # note the variance is independent of the hidden variable,\n",
    "        # it is fixed, so we just \"extend\" the dimension of the vector\n",
    "        # to fit the output shape\n",
    "        return mu, sigma_squared.unsqueeze(0).repeat(z.shape[0], 1)\n",
    "    \n",
    "    # compute the parameters of the Bernoulli conditional distribution p(z | x)\n",
    "    # - shape of input x: (batch size, n visible)\n",
    "    # - output shape: (batch size, n hidden)\n",
    "    def p_z_given_x(self, x):\n",
    "        mu = self.d + ((x/self.log_sigma_squared.exp()) @ self.W)\n",
    "        return torch.sigmoid(mu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check shapes of your RBM!\n",
    "batch_size = 15\n",
    "\n",
    "machine = RBM()\n",
    "\n",
    "x = torch.ones((batch_size, machine.n_visible))\n",
    "f = machine(x)\n",
    "\n",
    "if f.shape != (batch_size,):\n",
    "    print(\"Expected shape: \", (batch_size,))\n",
    "    print(\"Returned shape: \", f.shape)\n",
    "    raise RuntimeWarning(\"The shape of the log-partition tensor is wrong!\")\n",
    "\n",
    "z = torch.ones((batch_size, machine.n_hidden))\n",
    "mu, sigma_squared = machine.p_x_given_z(z)\n",
    "\n",
    "if mu.shape != (batch_size, machine.n_visible):\n",
    "    print(\"Expected shape: \", (batch_size, machine.n_visible))\n",
    "    print(\"Returned shape: \", mu.shape)\n",
    "    raise RuntimeWarning(\"The shape of the Gaussian mean tensor is wrong!\")\n",
    "\n",
    "if sigma_squared.shape != (batch_size, machine.n_visible):\n",
    "    print(\"Expected shape: \", (batch_size, machine.n_visible))\n",
    "    print(\"Returned shape: \", sigma_squared.shape)\n",
    "    raise RuntimeWarning(\"The shape of the Gaussian variance tensor is wrong!\")\n",
    "\n",
    "\n",
    "mu = machine.p_z_given_x(x)\n",
    "if mu.shape != (batch_size, machine.n_hidden):\n",
    "    print(\"Expected shape: \", (batch_size, machine.n_hidden))\n",
    "    print(\"Returned shape: \", mu.shape)\n",
    "    raise RuntimeWarning(\"The shape of the Bernoulli parameter tensor is wrong!\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training loop\n",
    "\n",
    "During training, we must approximate the loss via Markov Chain Monte Carlo.\n",
    "It usual to rely on contrastive divergence in the case of RBM training:\n",
    "\n",
    "- we use a single sample per batch data point in the Monte Carlo approximation of the partition function contribution to the loss\n",
    "- we start one Markov Chain per data point in the batch\n",
    "- we only do a few steps, e.g. in general one! this is parameterized by k below\n",
    "- we take the last sample from each chain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkRUlEQVR4nO3deXhU5f3+8fdnZrInEEIChCUS9kVZY0RBcUFEa0WLWqhfSy0Vv27VWluXbmrbn7Yu1bqCK/ZrcautWqVILQIqAkH2JUDYBIGEsCQsCVme3x8ZLFJSErKcOcn9uq65MnNm4tyHA7fPPHMWc84hIiL+E/A6gIiInBgVuIiIT6nARUR8SgUuIuJTKnAREZ8KNeabpaamus6dOzfmW4qI+N7ChQt3OufSjl7eqAXeuXNncnJyGvMtRUR8z8w2HWu5plBERHxKBS4i4lMqcBERn1KBi4j4lApcRMSnVOAiIj6lAhcR8SlfFPi89YU8N2c9lZU69a2IyGG+KPC/L93Gb95bxcQ/LUTnLxcRqeKLAr9vdF9uO78H/1y1g8/W7/I6johIRPBFgZsZE8/qQkpCNJNn52kULiKCTwocIDYqyIRhmczMLeCmqYs0Jy4izV6jnsyqrq4f3pX8ohKmzv+C95ZuY1PhAX5+cW9iQkGvo4mINDrfjMABAgHj3tEnk/ubUUwYlsmfPtvEhY/NYfveEq+jiYg0Ol8V+GFmxi8u7sOL15xKflEp331hHsUlZV7HEhFpVL4s8MPO6dmGSVcPJq9gP7e9voSDhyq8jiQi0mh8XeAAQ7ul8vNv9GbGyh2MeGQWM1bu8DqSiEij8H2BA1wzNJPXJg4hMSbEtS/n8OEqlbiINH1NosABTuvSmrdvGkrv9Bb85M2lbNl9wOtIIiINqskUOFTtK/7EdwZSVlHJD6bkkF+svVNEpOk6boGbWayZzTezJWa2wszuDS9/ycw2mNni8G1Ag6etga5piTx91WA2Fu7nm49/zMJNu72OJCLSIGoyAi8FznXO9QcGAKPMbEj4uZ845waEb4sbKGOtDeueyl9vGEpMKMjYyXOZs7bA60giIvXuuAXuquwLP4wK3yL+GPbe6S1496ZhdE1LZOLLC7l56iIKiku9jiUiUm9qNAduZkEzWwzkAzOcc/PCT/3WzJaa2R/MLKaa351oZjlmllNQ0Lgj4ZbxUbx0TTYXnZLOjJXb+c6zn/HQ9FztLy4iTUKNCtw5V+GcGwB0BLLN7GTgLqAXcCqQAtxRze9Ods5lOeey0tLS6id1LbRrGcvDV/Zn8tVZ7C8t54mZ63hz4ReNnkNEpL7Vai8U59weYCYwyjm3LTy9Ugq8CGQ3QL56c1aPND696zz6tm/B1Plf6JS0IuJ7NdkLJc3MksP344DzgdVmlh5eZsClwPKGi1l/xmZnsHJbEZNn63S0IuJvNTmdbDowxcyCVBX+6865v5vZv8wsDTBgMfC/DRez/lwxuCOzcgu4f9pqZq0p4MnvDKJVQrTXsUREas0acyohKyvL5eTkNNr7Vcc5x6sLvuBX76zglA4teeUHpxEbpXOKi0hkMrOFzrmso5c3qSMxa8rMGJedwR+uHMDnm3fzo9cWazpFRHynWRb4Yd/ol87PLurNtOXb+X/vr/I6johIrfjqkmoNYcKwTLbsPshzH29gwabdPPrtAWSmJngdS0TkuJr1CBz+fXWfX32zD+sL9nHvuyu8jiQiUiPNvsABggHjmqGZ3HRONz7KLeDluRu1n7iIRDwV+BHGn9GZYd1S+eXbK3hg2mqv44iI/Fcq8CPERgV5+fvZXHVaBpNmr+dPczd6HUlEpFoq8KMEAsY9l/RlRO82/OLtFVzyxMcs27LX61giIv9BBX4MUcEAT101mB+N6MGW3Qd54B/axVBEIo8KvBrRoQC3jOjOhGGZfLKukFtfXcSCjbu8jiUi8hUV+HGMy84gKTbE20u+5Lo/LWRdfrHXkUREABX4caUkRDP3rvP44NazOFReyYhHZvPkzHVexxIRUYHXRGJMiO5tk/jwx8O56JR2PPRBLp/m7fQ6log0cyrwWmjbIpYHL+9PZmoCt766mMJ9usamiHhHBV5LCTEhHh83kD0Hy7j1tcWUV1R6HUlEmikV+Ano274lvx7dlzlrd/KE5sNFxCMq8BP07VMzuLhfOs/MymNHUYnXcUSkGarJNTFjzWy+mS0xsxVmdm94eaaZzTOzdWb2mpk1u+uS3TGqFxWVjisnzWVuXqHXcUSkmanJCLwUONc51x8YAIwysyHA74A/OOe6AbuBCQ2WMkJ1Sonn2e9mUVHp+PnflumqPiLSqI5b4K7KvvDDqPDNAecCb4aXT6HqyvTNztk923D7yJ7kFezngkdnc9dby6hQkYtII6jRHLiZBc1sMZAPzADygD3OufLwS7YAHar53YlmlmNmOQUFBfUQOfJcdEo6HZLj2F5UwtT5m5kwZQFf7jnodSwRaeJqVODOuQrn3ACgI5AN9KrpGzjnJjvnspxzWWlpaSeWMsJFhwL8/eZhzLv7PO4b3Zf5G3YxYUqORuIi0qBqtReKc24PMBM4HUg2s8PX1OwIbK3faP7SKiGa+OgQ3z29M78b049V24p46/MtXscSkSasJnuhpJlZcvh+HHA+sIqqIr88/LLxwNsNlNF3Lu6XzoBOyTz0QS77S8s1EheRBlGTEXg6MNPMlgILgBnOub8DdwC3mdk6oDXwfMPF9Bcz42ff6M2OolIG/XoG17y0QNfYFJF6FzreC5xzS4GBx1i+nqr5cDmGUzun8L0zOvPZ+kJmryngk3WFDOue6nUsEWlCdCRmA7rnkr68fdNQOiTH8Zv3VlJaXuF1JBFpQlTgDSwmFOTeS/qyensxN77yOR+v1WloRaR+qMAbwYg+bfnhud2Ym1fI/zw/j5//bZnmxEWkzlTgjeS2kT1Z9MuRfH9oJv/32WZembfZ60gi4nPH/RJT6k90KMDPv9GbdQX7+OXby4mLCjJmcEevY4mIT2kE3sgCAeOZ/xnE6V1bc+dbS1m9vcjrSCLiUypwD8RHh3h83CBaxkXxw6mL2HPgkNeRRMSHVOAeSUmI5rGxA9m484DOmyIiJ0QF7qGh3VJ5YMwpLNy0m1cX6EtNEakdFbjHLhvYgdMyU3hwei6792sqRURqTgXuMTPjvtEnU1xSzoQpC3jkg1ztIy4iNaLdCCNAz3ZJXD+8K89/vIHPN+8hFAxwzdDOJMVGeR1NRCKYRuAR4vYLerLyvgsY2actj8xYwzkPfURBcanXsUQkgqnAI4iZ8eRVg5h89WD2Hizj/vdXeR1JRCKYCjzCRAUDjOzbjuvO6spbi7YyN6/Q60giEqFU4BHqxnO60bFVHLe9vpgnZ66jrKLS60giEmFU4BEqLjrIY2MH0KZFLA9Oz+X7Ly3Q+cRF5Gtqck3MTmY208xWmtkKM7slvPweM9tqZovDt4saPm7zMvikFN6+cSi/G3MKc9bu5O63lnsdSUQiSE12IywHfuyc+9zMkoCFZjYj/NwfnHMPNVw8Afj2qRls3nWAJ2fmMfGsLvRsl+R1JBGJAMcdgTvntjnnPg/fL6bqivQdGjqYfN34MzpjBtOWb/M6iohEiFrNgZtZZ6oucDwvvOgmM1tqZi+YWatqfmeimeWYWU5BQUHd0jZjbZJiyTqpFdOWbdeRmiIC1KLAzSwR+Atwq3OuCHga6AoMALYBDx/r95xzk51zWc65rLS0tLonbsZGD+hA7o5izvz9TEY/+YnOnSLSzNWowM0siqryfsU59xaAc26Hc67COVcJPAtkN1xMAfhOdgYPX9GfXu1asGzLHp6eled1JBHxUE32QjHgeWCVc+6RI5anH/GyywDtItHAAgFjzOCOPDc+i8sGduSlTzfy9Ed5VOpc4iLNUk32QhkKXA0sM7PF4WV3A+PMbADggI3AdQ2QT6px10W92LW/lN/9YzVm8L/Du3odSUQa2XEL3Dn3MWDHeOr9+o8jNZWaGMML3zuVG175nIem59I7vQXDe+g7BpHmREdi+piZ8cCYfvRom8TEl3N4c+EWryOJSCNSgftcy7go/u8HpzEwI5nb31jCc3PWex1JRBqJCrwJSEmI5pUfDOHCk9vxm/dWMWuN9rcXaQ5U4E1EMGD84dsD6N4mkTveXMraHcVeRxKRBqYCb0Jio4I8fGV/9peWc8Gjs3lmVp6O2hRpwlTgTUy/jsnM+uk5XHhyOg9MW83rOV94HUlEGogKvAlKSYjm8XEDye6cwv3TVrOpcL/XkUSkAajAm6hAwPjtZSdTUeEY9egcXZpNpAlSgTdh3dsm8cFtZ9E+OZabpy4iv6jE60giUo9U4E1cess4nrxqEPtLyxn/4gL2HizzOpKI1BMVeDPQq10Lnrl6MOvyi7l2Sg4lZbq2pkhToAJvJob3SOPhKwewYNMu7vv7Sq/jiEg9qMnZCKWJuKR/e1Zs3cuk2etpnRDNzed2Jzqk/4eL+JUKvJm5bWQPtheV8Pi/1hEMGLeO6OF1JBE5QRp+NTMxoSCPjR3IiN5tmPLpRg4cKvc6koicIBV4M3Xd8K7sPlDG7W8sobhEe6aI+JEKvJk6tXMKd4zqxT+Wb+ebj3/M+oJ9XkcSkVqqyTUxO5nZTDNbaWYrzOyW8PIUM5thZmvDP1s1fFypT9ef3ZVXJ57OnoNl3PGXpTrxlYjP1GQEXg782DnXBxgC3GhmfYA7gQ+dc92BD8OPxWeyM6tG4gs27ubHry/R0ZoiPlKTa2JuA7aF7xeb2SqgAzAaODv8sinAR8AdDZJSGtSVWZ1YvnUvbyzcwoFDFTxz9WCvI4lIDdRqDtzMOgMDgXlA23C5A2wH2lbzOxPNLMfMcgoKdKWYSBQMGL+97BSuO6sL01duZ9W2Iq8jiUgN1LjAzSwR+Atwq3Pua//CXdXk6TEnUJ1zk51zWc65rLQ0XTU9kn3vjM4kRoe45ImP+eHURazL1xebIpGsRgVuZlFUlfcrzrm3wot3mFl6+Pl0IL9hIkpjaZ0Yw99/OIyrh3RmZm4+lzzxMVPnb6ayUl9uikSimuyFYsDzwCrn3CNHPPUOMD58fzzwdv3Hk8Z2UusEfvnNPvzztuH069iSu95axu1vLtEeKiIRqCYj8KHA1cC5ZrY4fLsIeAA438zWAiPCj6WJaNsilqnXDuHmc7vx1udbmTpfl2YTiTQ12QvlY8Cqefq8+o0jkcTM+NGIHszfsItHZuRy6cD2xEfr9DkikUJHYsp/FQgYP7mgJzv3HWLy7PVexxGRI6jA5biyOqdwcb90Hv3nWibNyqO0XBeEEIkEKnCpkYeu6M95vdpw/7TVXPvyQu2ZIhIBVOBSI7FRQZ4bn8WvvtmH2WsKmDJ3o9eRRJo9FbjUmJnxvTM6fzUS1xGbIt5SgUutmBkPjOlHUkyIi/44h3vfXeF1JJFmSwUutZaWFMPfbhzKmEEdefGTjbz1+RavI4k0S9qpV05Ip5R4HvjWKWwuPMBP3lxKpYPLB3f0OpZIs6IRuJywUDDAi9ecypAuKdz91jLW7Cj2OpJIs6IClzpJiAnx2NiBJMaGuHLSXJ6bo4N9RBqLClzqLDUxhj9NyKZfx2R+894qXluw2etIIs2CClzqRd/2LXlhfBbDuqVy37srKSgu9TqSSJOnApd6EwoGuG90X0rKK/nVO8tV4iINTAUu9apLWiI3nt2V95dt58LHZpNXoKv6iDQUFbjUu9tG9mTaLWcC8O1Jc/lk3U6PE4k0TSpwaRC901vw6sTTSY6P5n+en8dTH63zOpJIk6MClwbTrU0i79w0lIv7tef3/8jl/WXbvI4k0qTU5JqYL5hZvpktP2LZPWa29ahLrIn8h/joEA9f0Z+BGcn8+PUlfLa+0OtIIk1GTUbgLwGjjrH8D865AeHb+/UbS5qS6FCAyVdn0aFVHD+YksPGnfu9jiTSJBy3wJ1zs4FdjZBFmrC0pBimfD+bYMC4/pXPKdynXQxF6qouc+A3mdnS8BRLq+peZGYTzSzHzHIKCgrq8Hbidx2S4/jjuIGsL9jHpU99wqw1+vsgUhcnWuBPA12BAcA24OHqXuicm+ycy3LOZaWlpZ3g20lTMbxHGlMnDiEqEGD8C/P54dRFlFVUeh1LxJdOqMCdczuccxXOuUrgWSC7fmNJUzYooxXTbj2TW87rzjtLvuSed1boGpsiJ+CEzgduZunOucP7hF0GLP9vrxc5WkwoyI/O70FJeQWTZq2noLiUx8YOJC466HU0Ed+oyW6EU4G5QE8z22JmE4Dfm9kyM1sKnAP8qIFzShN156he/OLiPsxYtYNrXpqv6RSRWjjuCNw5N+4Yi59vgCzSDJkZE4Zl0jIuitvfWMJzczZw/dldvY4l4gs6ElMiwuWDOzKyT1senL6aa16cz77Scq8jiUQ8FbhEjAcv7891w7vy0ZoCfv3uSq/jiEQ8XdRYIkbL+CjuGNULA576KI9e6UlcMzTT61giEUsFLhHntvN7kFewj3vfXUlpeSXXndUFM/M6lkjE0RSKRJxQMMAfxw3kG/3SeWDaaq54Zi5Lt+zxOpZIxFGBS0SKCQV5fOxAHvjWKWws3M+lT37CzNx8r2OJRBQVuESsQMAYm53Bv24/m57tWnDL1EVMmpVHufYVFwFU4OIDLWKjmHz1YHqlt+D+aauZOn+z15FEIoIKXHyhU0o8r00cQnZmCo99uI69B8q8jiTiORW4+IaZcdeFvdhz4BBjnvmU95dtwzmdBEuaLxW4+MrAjFa8PCGbg4cquOGVz/m/eZpOkeZLBS6+c0bXVOb89BzO6Nqah6bnsi6/2OtIIp5QgYsvBQLGfaP74pxj5B9mc+7DH7F8616vY4k0KhW4+Fa3Nkn86/azueW8HhQdLOeOvyylQheGkGZEBS6+lpoYwy0junPPJX1Y8WURv/vHaq8jiTQaFbg0Cd84JZ3vnn4Sk2evZ8zTn7Jmh+bFpelTgUuTYGbc882+/PLiPmwqPMBVz81j5ZdFXscSaVA1uaTaC2aWb2bLj1iWYmYzzGxt+Gerho0pcnyBgPH9YZn8+drTcA4ueeJj3l681etYIg2mJiPwl4BRRy27E/jQOdcd+DD8WCQi9GibxIwfncWgjFb85M2lzM0r9DqSSIM4boE752YDu45aPBqYEr4/Bbi0fmOJ1E2rhGieuXownVrFMf6F+fxtkUbi0vSc6Bx4W+fctvD97UDb6l5oZhPNLMfMcgoKCk7w7URqLyUhmr9cfwYDM5K59bXFXPbUJyzcdPRYRMS/6vwlpqs6GUW1O9865yY757Kcc1lpaWl1fTuRWkmOj+ZPE07j1hHd2b63hJv+vIiiEp0IS5qGEy3wHWaWDhD+qTPtS8SKDgW4dUQPnrpqEDuKSrj0iU+Yv0EjcfG/Ey3wd4Dx4fvjgbfrJ45IwxmY0YoXr8mmwjm+8+xnTPl0o85mKL5Wk90IpwJzgZ5mtsXMJgAPAOeb2VpgRPixSMQb3iONd28exvAeafzqnRU8OD3X60giJ+y4V6V3zo2r5qnz6jmLSKNoERvFs9/N4u6/LuOpj/LompbImMEdvY4lUms6ElOapaqzGZ7M6V1a8+M3lnD/tFWUlFV4HUukVlTg0mxFhwK89P1TGZfdiUmz1jN28mfaQ0V8RQUuzVpMKMj93+rH01cNYvnWvVz46ByenLmOjTv3ex1N5LhU4CLAhaek8+drh9AqIYoHp+dy6VOf6GRYEvFU4CJh2Zkp/P3mM5l5+9nERQW54ZWFmheXiKYCFzlKZmoCD13Rn42FB7j7r8t4b+k2lm3R5dok8hx3N0KR5mhot1RuOLsrT8/K463Pt5IUG+K9m88ko3W819FEvqIRuEg1fjqqFx/cehbPj8/CgBv//Dml5ZpSkcihAhf5L7q3TeK83m15+MoBLNu6l9++t8rrSCJf0RSKSA2c36ct156ZybNzNhAXFWRk33YMPkkXohJvaQQuUkM/HdWL07u0ZtLs9YydPJd563WlH/GWClykhqKCAf587Wl8fMc5dEiO49uTP+Pq5+fpoB/xjApcpBbMjI6t4nnz+jO4fWQPFn+xhwsenc2v3l7OR7k6Lb40LhW4yAlITYzhpnO788/bhnNe7zZMnf8FE6bksOJL7S8ujUcFLlIHbVvE8tRVg5n/s/NoFR/FtVNyeGh6LvtLy72OJs2AClykHiTHR/PHcQPplBLPkx+t4xt/nMMXuw54HUuaOBW4SD05o2sqr113Oq9eO4Rd+w9x+TOfMmtNgdexpAmrU4Gb2UYzW2Zmi80sp75CifjZaV1a89p1p5MUG8X4F+Zz+xtL2Fyo0bjUP6vLRV3NbCOQ5ZzbWZPXZ2VluZwc9bw0DyVlFTz6z7W88PEGzODOC3sRHx3kisGdCATM63jiI2a20DmXdfRyHYkp0kBio4LceWEvvndGZ37w8gLufXclAMu27uVHI3rQOjHG44Tid3UdgW8AdgMOmOScm3yM10wEJgJkZGQM3rRp0wm/n4hfHThUzsovi3hv2TZe/GQj0cEAN57TjZvP7abRuBxXdSPwuhZ4B+fcVjNrA8wAbnbOza7u9ZpCkebOOcfKbUU8M2s97y75khvO7spPR/XyOpZEuAaZQnHObQ3/zDezvwLZQLUFLtLcmRl927fkj2MHkBgT5KmP8theVMKQzNa0SojmvF5tNCKXGjvhAjezBCDgnCsO3x8J3FdvyUSaMDPj16NPJjk+mknhi0YA3H1RLyae1dXjdOIXdRmBtwX+amaH/zt/ds79o15SiTQDoWCAO0b14ofndmfnvlJ++94qHpyeS/vkOIb3SCMpNsrriBLh6jQHXluaAxep3t4DZXzvpfks2ryHgMFt5/fgxnO6ER4kSTOm3QhFIlzL+CimXjuE95Zu41+5+Tz0wRq27inh3F5tOKtHKjGhoNcRJcKowEUiSGxUkDGDO/KtQR3okBzH5NnrmTp/M4NPasVDV/QnMzXB64gSQTSFIhKhnHOs2bGP5Vv38ou3l3OovJIJZ2Zy63k9iIvWaLw50RSKiM+YGT3bJdGzXRJn9kjloem5TJq1nmnLtjOyT1u+c1oGXdISvY4pHtIIXMRHPs3byUPTc1n+ZRHOObqmJXLb+T0Y2bed19GkATXIkZi1pQIXqR8FxaU8N2c9M3PzySvYT7+OLbnh7G6c36et19GkAajARZqgfaXl3P/+KuauL2Rz4QG+c1oGaYkx9Epvwdk904gK6pT/TYHmwEWaoMSYEL+97BSKS8q4+6/LeW3BF5SWVwIwondbJl09mKAOzW+yNAIXaULKKiqpqHRM+XQj909bTbc2iXRuncClA9tzcb/2XseTE6QRuEgzEBUMEBWE64Z3pUVcFP9Yvp3cHUXc9OcdTFu+nayTWnFm91S6tUnyOqrUA43ARZq4sopKHv9wLc9/vIH9hyqIjw7y0wt6kpYUS7+OLemUEu91RDkOfYkp0syVlFWwbW8Jt7y6iKVb9gLQIjbELy7uw5AurVXkEUwFLiJA1RGeGwsPsHNfKbe/sYRNhQcIBowxgzpwSf8OdElLoH1ynNcx5QgqcBH5DyVlFWzYuZ/Xc77glXmbORTeg+W0zBQuOiWdM7un6mjPCKACF5H/qqC4lLX5xSzcuJu3Fm1lw879AGR3TuHSgR3Izkyha1qCTm/rARW4iNSYc44tuw/y3rJtvDp/MxsLDwCQEB2kQ6s4+qS3YFx2Bid3aElCjHZma2gqcBE5Ic458gr2s3DTLlZvL2br7oPMzSukuLScmFCAb/Zvz/l92jKwUzIVzpGWGENIR4DWqwbZD9zMRgGPAUHgOefcA3X574lI5DEzurVJpFubf8+FF5WU8VleITNz83l3yTbeXLjlq+eSYkKMGdyR3ulJZKQkkNE6nnYtYnVEaAM44RG4mQWBNcD5wBZgATDOObeyut/RCFyk6TlUXknOxl2s2VFMVCjAp3mFTF++nfLKf3dLdDBAx1ZxZKYm0LdDS3bvP0R2Zgrtk+Po1S5J0zDHUe9TKGZ2OnCPc+6C8OO7AJxz91f3OypwkeahvKKSL/eUsHnXATbt2s/mXQfYXHiA3B3FrC/YT2xUgJKyqj1eQgFjUEYrdh04RGJMiAOHyokOBRgzqCPpLeMoq6ikoLiU9Jax9G3fkkAA2iTFUlJeQVxUkA0797OvtJwWsVG0iAsRHx1i5ZdFtIyLome7qiNOS8oqWJe/j9LyStokxeBc1aeIUNBISYimTVIsRSVl5BeV0rZFzFcXlK6odBQUl7JzXyl90ltUvaa4lGDAyGydwN6DZSz6Yjct46Lok96SYMDYvreE9smxlFU4dh84RElZBZ1bJxCowyeQhphC6QB8ccTjLcBpx3jjicBEgIyMjDq8nYj4RSgYIKN1PBmt4xlG6teeKymrIBQwVnxZREFxKQs27uKTvJ10bp1ASVkFbZJi2FFcyr3vVvthnphQgNLySqJDga92fTyWhOggsVFB9hwso6Ky+sFqUmyIg4cqvvrUkBgTIhQ09hwo++o10cEAhyr+/V5xUUEOhc89c1jAoNKBGRw5Nk5vGcvDV/TnjG5f/7Ooqwb/3OKcmwxMhqoReEO/n4hEttioqsvB9e+UDMCIY5zD3DnHki17CQWM6FCAlIRotuw+yKptRTgHa/OLaZ0Qza79ZXROjadTSjzFJeXsPVjG/tJyTkqJp2BfKZsKD1BSVkFyfBS901uQEBMiv6gEM6NlXBSVlY4v95awuXA/CTEhurdNJL+olO1FJRwqryQtKYbEmBAt46LI3V5Mu5axtG0Ry8FDFazeXkxiTJCh3VIpLiln1bYiSssraZ8cx5d7DpIQEyI5vmok/8+VOxrkSNe6FPhWoNMRjzuGl4mI1ImZMSBc8IelJsb8x7JIcqz/ER02LrthZh/qsq/PAqC7mWWaWTQwFninfmKJiMjxnPAI3DlXbmY3AdOp2o3wBefcinpLJiIi/1Wd5sCdc+8D79dTFhERqQUdLiUi4lMqcBERn1KBi4j4lApcRMSnVOAiIj7VqKeTNbMCYNMJ/noqsLMe43hJ6xKZtC6RSesCJznn0o5e2KgFXhdmlnOsk7n4kdYlMmldIpPWpXqaQhER8SkVuIiIT/mpwCd7HaAeaV0ik9YlMmldquGbOXAREfk6P43ARUTkCCpwERGf8kWBm9koM8s1s3VmdqfXeWrLzDaa2TIzW2xmOeFlKWY2w8zWhn+28jrnsZjZC2aWb2bLj1h2zOxW5Y/h7bTUzAZ5l/zrqlmPe8xsa3i7LDazi4547q7weuSa2QXepD42M+tkZjPNbKWZrTCzW8LL/bhdqlsX320bM4s1s/lmtiS8LveGl2ea2bxw5tfC10/AzGLCj9eFn+9c6zd1zkX0japzjecBXYBoYAnQx+tctVyHjUDqUct+D9wZvn8n8Duvc1aT/SxgELD8eNmBi4BpgAFDgHle5z/OetwD3H6M1/YJ/z2LATLDf/+CXq/DEfnSgUHh+0nAmnBmP26X6tbFd9sm/OebGL4fBcwL/3m/DowNL38GuD58/wbgmfD9scBrtX1PP4zAs4F1zrn1zrlDwKvAaI8z1YfRwJTw/SnApd5FqZ5zbjaw66jF1WUfDbzsqnwGJJtZeqMEPY5q1qM6o4FXnXOlzrkNwDqq/h5GBOfcNufc5+H7xcAqqi4y7sftUt26VCdit034z3df+GFU+OaAc4E3w8uP3i6Ht9ebwHlmVqtL1/uhwDsAXxzxeAv/fQNHIgd8YGYLzWxieFlb59y28P3tQPUX1Is81WX347a6KTyt8MIR01i+WY/wx+6BVI32fL1djloX8OG2MbOgmS0G8oEZVH1C2OOcKw+/5Mi8X61L+Pm9QOvavJ8fCrwpGOacGwRcCNxoZmcd+aSr+gzly/05/ZwdeBroCgwAtgEPe5qmlswsEfgLcKtzrujI5/y2XY6xLr7cNs65CufcAKou8p4N9GrI9/NDgW8FOh3xuGN4mW8457aGf+YDf6Vqw+44/DE2/DPfu4S1Vl12X20r59yO8D+4SuBZ/v1RPOLXw8yiqCq8V5xzb4UX+3K7HGtd/LxtAJxze4CZwOlUTVkdvnzlkXm/Wpfw8y2Bwtq8jx8KfAHQPfxNbjRVk/3veJypxswswcySDt8HRgLLqVqH8eGXjQfe9ibhCaku+zvAd8N7PQwB9h7xkT7iHDUPfBlV2wWq1mNseC+BTKA7ML+x81UnPE/6PLDKOffIEU/5brtUty5+3DZmlmZmyeH7ccD5VM3pzwQuD7/s6O1yeHtdDvwr/Mmp5rz+5raG3+5eRNW303nAz7zOU8vsXaj61nwJsOJwfqrmuj4E1gL/BFK8zlpN/qlUfYQto2r+bkJ12an6Fv7J8HZaBmR5nf846/GncM6l4X9M6Ue8/mfh9cgFLvQ6/1HrMoyq6ZGlwOLw7SKfbpfq1sV32wboBywKZ14O/DK8vAtV/5NZB7wBxISXx4Yfrws/36W276lD6UVEfMoPUygiInIMKnAREZ9SgYuI+JQKXETEp1TgIiI+pQIXEfEpFbiIiE/9f1JqQ7ymtwU2AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "batch_size = target_samples.shape[0] # we just feed all data at once, you can try with smaller batch size\n",
    "n_epoch = 300 # number of epochs\n",
    "k = 1 # number of step to take in the Gibbs sampler for contrastive loss\n",
    "\n",
    "n_visible = target_samples.shape[1]\n",
    "n_hidden = 100\n",
    "\n",
    "\n",
    "machine = RBM(n_visible, n_hidden)\n",
    "# we fix the variance of the observed variables,\n",
    "# its easier to train this way!\n",
    "machine.log_sigma_squared.requires_grad_(False)\n",
    "# optimizer\n",
    "optimizer = torch.optim.SGD(machine.parameters(), lr=1e-3)\n",
    "\n",
    "\n",
    "losses = list()\n",
    "for _ in range(n_epoch):\n",
    "    for i in range(0, target_samples.shape[0], batch_size):\n",
    "        optimizer.zero_grad()\n",
    "        data_samples = target_samples[i:i+batch_size]\n",
    "        \n",
    "        with torch.no_grad():\n",
    "            # generate sample from your model!\n",
    "            # and store the in a variable called model_samples\n",
    "            mc_samples = []\n",
    "            mc_samples.append(data_samples)\n",
    "            for _ in range(k):\n",
    "                probs = machine.p_z_given_x(mc_samples[-1])\n",
    "                z = torch.bernoulli(probs)\n",
    "                mu, sigma_squared = machine.p_x_given_z(z)\n",
    "                x = torch.normal(mu, torch.sqrt(sigma_squared))\n",
    "                # print(x.shape)\n",
    "                mc_samples.append(x)\n",
    "\n",
    "        # the loss function is quite simple\n",
    "        loss = (machine(data_samples) - machine(mc_samples[-1])).mean()\n",
    "        losses.append(loss.item())\n",
    "\n",
    "        loss.backward()\n",
    "        # we use gradient clipping to stabilize training\n",
    "        torch.nn.utils.clip_grad_norm_(machine.parameters(), 1)\n",
    "        optimizer.step()\n",
    "    \n",
    "# plot the training loss\n",
    "plt.plot(np.arange(len(losses)), losses)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sampling from the trained RBM\n",
    "\n",
    "We now turn to sampling from our trained model to check if it models the data correctly.\n",
    "We will run many markov chains in parallel and do several visualization.\n",
    "In the first one, we will just plot the samples at the timestep of each markov chain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAAsTAAALEwEAmpwYAABlZklEQVR4nO29bZAkyXke9tbdidjdC/IWt43twS55AmDTlkHvdJm47QIoYkdhK4IQHQZJC4ygQqLFCMggzZvKwsL+AQrbVUf+UZCMQN8CkBfEwQzJWGlNgH8EOnYDIfJ6SEfsAPRheQeAOnxyKS651IEgeYc428JNV77+kfVmvflR1dUz09M9M/lEVMxUd3XWV+aTb76fESJCQEBAQMDRxwPLvoCAgICAgINBIPyAgICAY4JA+AEBAQHHBIHwAwICAo4JAuEHBAQEHBM8tOwLaEKv18M3vOENy76MgICAgEOFz3/+899CxNf5vltZwn/DG94Azz777LIvIyAgIOBQIYqif9/0XVDpBAQEBBwTBMIPCAgIOCYIhB8QEBBwTBAIPyAgIOCYIBB+QMABwE5ZFVJYBSwDgfADAhaMJ58EuHy5JnlEtf/kk8u8qoDjiED4AQELBCLASy8BXL1ak/7ly2r/pZeCpB9wsFhZP/yAgKOAKAIYj9X/V6+qDQAgy9TnUbS8aws4fohWNR/+448/jiHwKmApQDSZ2N7fZZMPsPW0lIHsAxaDKIo+j4iP+74LKp2AAI4FKNypCQ5+ioCAg0Ig/IAAwgIU7ryJLFOSfZaZpzg0CK5Ghx5Bhx8QQFiAwj2KAE6fNpugU5w+fYjUOk8+qSY9ugmayU6fDu5GhwhBhx8QwIAIEMH+K9wXYBY4uHPYy5Tx2N0/NDPX0UebDj9I+AEBFZ58EuCv/wrhqegycPr67Nsuw1u3V5vUFiqAB1ejI4Ogww84tuCLW0SAWzcR3vjhyxB96CqgyCBLJTwFGbz1c1cB37t7hbtjB5YIly8DFIV5IbtdbB+Irz8nfUIg+0OHIOEHHHrsRpXhk4hf3YkghufgDyCGH/zQBwEgAoAPwrt6vwvnn3vOaLTrOTkZAwCMH3kSfu/TL8HV58YQxxEUOcID//NlwEdOw+WXn9yVRB5FAOMPIgBETABHyLJo/zi5ydUokP7hAiKu5PaWt7wFAwIIUvr3iwIxy+p9KdV+UbS3lWWIAOpvWdK+xOs99cUYMgSQOInZgdVJ5j1nfT6JY1DtXe+p9req9rdigQCybte+4TZUFyRLiYqZ1Xlk3vIQ5oH9wHz7ASsDAHgWG3h16cTetAXCDyA0EWye756H+LG0CYG4+URNyrRJUTe2W+6TEg0y5u1vQ4JjEJgJWZP9rFnLcyNbcWa0vxWrSWBfsJuZNWApCIQfcGgxi2Br6bzeugqdZWn+bnOzJmX+hUil0Z5vsphF9ubxZvtjEObEMqf0LKelXinQRuS/rwJ40zIrYKUQCD/gUGMWwUqTPzvxkK9NAGmodWi7seZKyr5z+vjQ4e9SOuQ8BqFJf+5Zq5K8i5E5e8k8DwL4MUUg/ICVxDwCYxOpzytt27+hVUIBhUm8QqAUAu/2E6dR3zmTRKmEfBoPrQ0p6x/6dPh7mrXi2Px9HKOclrPbmAfzvrCuxwbsKwLhB6wc5lEJN5E6V+fMq8On85clYiZqvfc2JLg1UGSv1CxCMXl1Yb5zCMbVRPr2MfpaqhPnI1lPAkKgHCadZy2jvbJE6SH7tjYaubiNpOd5YUHfv1QEwg9YKUipSNYgamvfODZrJnUy3HbmFtY4kT2A+itFZhCn1qlbROjjM5FKTAzObtGfs5UC5jniYFCfryxbCbsozNWElIg3L+Ym4dNMmOfOqflEx5/XM5c8N0UPch5LdfDoWToC4QesFGReoBSZJlvyLLmVFF4+mCUwdtEeSGk2RN4w2wk77yxjANuXEmu9fqHup5zO6RZpq2T4kiWOHcL2rSZE6nr9YJYhTqfOzNd8Oma3aCLpeXRnu9GzBewbAuEHrAyKvNZbS5EhdyPk7o82HKNpg7shfa4nh1JiUShiJAmeJpvthJ13lruPNevQfUjmG0p6+XncIvnzoG0rroywHpQlJ/36XPfODGomB9CrBpto+QJCGalpX+Jk4BqOjetvmxB9k+O8NomAfcHCCR8Afh0AvgkAX2r4PgKADwHA1wHgCwDwg7PaDIR/9FALfrIOaNKeKpnj/qh/lxcG+ZDx05agJxsFbsW1frycqvN87PW5MqwOG9Q2s4wB1vfc02YrzrDcKR3Pni5ukfx58N82BWAVhVI9cXfSAgrchgSlEChTk7DlMDGXRRVqI3WmSV+kEm/0rd+Xsl4k2MsLAMQ0NaV/vuQKEv7ScBCEfwkAfrCF8H8UAG5VxP9WAPjcrDYD4R9N1LxhkpxI/Tp8LgFvxUptwveLvJ4E6HOStImEXzgRI0BZncpyvm+Q4LuQGJG6r12DLOn39oOwrpu3K0vpqKDIuHytn7PDpfYssiX0rYFQKwWasIoCZSlRbJZ6ZaAjim3pnt2fSGVtVE5Tk/iHw3q/YXJ09gMWigNR6QDAG1oI/9cA4B+w/a8AwOvb2guEPxuH1fNNlv5o1kxIw9jKJwdbgiZSNzjaQ553IFZkGGcIUOp9r+QprRWGZ9+45mk9iTjtCmHo+A2yk7RCqdVBYxBoROFyYiXSr1Ynd6CewOKBIvxtsLx8APA+9PHu2SHKQVw9Y4FbgxS3IcFrZ3NX9w+A//qseR2TOEOxWeK984m+L0sv5D5Hig3ImSdStSILzjqLxyoQ/v8JAD/M9n8HAB5vay8QfjtW0fPNdmjxOLgYpEwkp3X6qXBUGOQV2ar2sNo3jysdYrsDMYrNstbpV+fVHizT2g5gqDQsCf9eL0aAKd6HNb0vp0zJniSGxKsld7b/zNtzvNtP9LMAkMotNEnM3BFE4msxcsn8X/fV918+sW7do7mPcawmkWp/DALTJ8xVyTPrAgfrpuH55rBQfC6kdlX1buT2w96DoSWz9w+JcHIYcWgIHwDeAwDPAsCzjz322IIfy+HFsj3ffCsLn7ugEOozWytSQIF3+0rvrDx1GMl5/N19nihjyDQx6+vwRrFmaKtbxKaSjpOhxD85n+DdtQTLKRGSWk08c6mon6lguv8sQzkt8V7PknABcAypSYyJ0qE3qW0ApJbQka1yjFQLVv6HLHUnsG0Y4kun1izCHzSS8+2hQLFZeuwoNOnQR1IL9V4vJntC8XgVBVX+wWMVCD+odPYZBzmYDJ164UaUpikaPui2mpdse/r3qWSBTRlmaYnbiaUHrtrmZE9SLc84mY+kQ6rXewK5Gui5yCTnSUW2ybDWXU8q+4CdLVM/U9tLx05lwKRnkyWr5+SsPCRmQiohXrjG5O0k83oOjUE9L5Pc1f29CD0EmLrqJWuT05rsJ3GGoyvmM37i56V7K6Wnw/V6iNNpHfhFvp5W3+E/CWS/eKwC4f+3ltH292e1Fwh/NhY1mJoInjtqVMKrJnqb5Gkjmx4P9BECcXSlYcaq/A5lXuhjC8g1Afd6iAASP3oiwwIKvNYvtDun8tIReHso8FZSYLlTS+K2y+SkUq/IUuKNNfM6aFVgPFOun/LNtvbN27otMyoLx1B51RSFX3r2eA5pN1JLf74NF/FF6HnJ/YUTA5fwU4GTt49wEme4fkHiww+jfjYFFNUzrjfu3bSdZPi1178d5YkTxkGvnOohXrpk9B/uhWQIJfthuD2sBqwDwEF46dwAgD8HgB0A+FMAeDcA/BwA/Fz1fQQA/xwAvgEAX5ylv8dA+DPhlfCFx8g4J7gg6yP4JjUukbrvcy780TX3eoijD5g/KEalJsatgcByqqRwRUS51bZSh9yuSFAbNlOmDpFSG0bNoCihXTqlRBSb1nWwCcaQbpuiTn0PxTdLgpsoTaap0vv7Jj87jJinUaiMwJNBVrW76VzDFx4c1Ofp9fD5k8N6fxDj//q6+pkqt/2amCvhXV/6raTQqzGfsZeeuXaJJccgKPRKjMhfxS8Uc/dNbyell7hsA9YKIQReHTH4OOdWovyqd5VTvaXdJi5r2+fb5qbr0EFSuu15Y0vASr+dVfsZMyiiFoIxz72Jw7g+mRsMKUkaPSefPt707GH6d3oo9qw4rIh0OPTOkjIv8EbfzYh5D87j0yfT2iWy8qW/u8YStlkqEhyNNNmr05SNXjoFjHAbhnjvjDLgXu8JHEOqjcxjyLB/VrVTT4hqo8fH3UtpRZSlrleSjGPMr5RGVojRFbUSQ3Bdah3CtjviPJ00WIMNBMI/gjCEHGkuufcyCLpoK3zxN/Y+P2bTET5tXXm7PpyCsuzzZEI6wUaGhFy5VXLj6+2hqNvcrMn+y6di5B49jlGWPUc9odJOYrksWl46Uiq1SFMe/G1IcCuuA6fGINS+LQXzRDhSYj6S+LGTNdnbqp47EGP8A6/iHRjgPTjvfAdQYlmiJuICCvN2WcAbBX3JtDKue565NkTH6h3b79r3LPdcsszX5jFHIPwjClt9Yxv+djsIbHVym8ai6budHfV3OHTU1wiAWI4KTfa2imDLCgIi1Qo/b60uECgpjUC1ffvhviEZaztErnT8tqR9H9YwXieferWqePpcbhxju2rmI0tH7zPaWgFPtlfMZCDw9kXzWsaQYTyQZvv8xDSR1P6qePd1F1FWivdvn1B//x9QOnbS7fsMuXQuUguRXYOrXhC495Dav9sfuhL+YFDn4K9WUWQHKXcagt32Iq0Ha3AjAuEfF3QcBG0raJ+hDQB12oMZwqze7/UQX3213rdsfJVapyb7F07EFTENPMQkMBlKw/FFedQ065JkWgc+2bl07KjS62dS5B45PBK1PkZo180CCrx+RtSk7CN8brBg5Gx7Gtl58Om589Wb8XI8Ewt5PKm4gB2vimoM5u/u981jnhmoFU+R16uiIndndG6DmAwEilTWK4xUYLlT6u/vwxoCTFtXS7uS1oOE34pA+KuI/fYy6DgI2lbQXHKmICFF8Eo6NTxniro9bdhEZegjHX0cI7797YgPPliT/He+U5N/HJPUnSHpovn187z0vlw75VTWBUpsaXOUG3ZPgObgodsXU8zS0qseK6cSnz+ZaLIDKPH6mUqy7ic6n33TxKNVS/o+a1XHVuzmwSepuywbXlbTubIMiytTdfi0tNpUZH81Mt1V+SY2S9yGRN8TPQvtLqsnpFJL8DIVKEe5nrTkMKnsLqnTPq08vKmfWwQVZ1gchA7/kHsABcJfNey3l0HHQdDlMDubJVdFaM8XiYbKhHS9dJ7pjiIVrhN+6CHEJ54wVwAXL6pbLqe13pcP/HIqNQHfSgr38UhXh09kzPXJly7VGSG3YWgcT+Qk+/3aUlkURi58miS2IUGeu8bZfK5KlQ6fVk36uYu6HVvqJ197PTm1WdBpq2Y2+eqOo8u/D2u4DUMcg8Azj0oc/dOpVgHRxldW24kw7ts8jmwdou4j09KZuOSmS/qm+07h9l36Letn3mHRNBHaeTl2gyPgARQIf5WwWwllltSR52bBjrLUft40eEhya1wI0MRQ1j7fehD2ephtTlW/p1wqGxut98OzMdq8SIbc2sDbkkRMSu9jkaWsPVo0gSuVQwEF8ngg9XWpPVRo4+oGI7UDPTOLnGtSEu5NtaRY9vHIraTAT55T3kh0fWMQeK1f4DOXCpRp6pY+NB9avdFNVsunV04pQid/edLl3x6Kxonj3pmBdzKbDNTKxtbbk+G7nErH7oKWXYUmP0Q0/1bPy3b55PUSvMOCf+CLBtwNSR/E6uEAEAh/1TCvDnKW1FEUuJ1UUmFeVCXzhni3n2AxKnE7UUZIUsvYQqiUaEi1UqITzUnSXZaWbmSl534o0KlJICWOFAJx45LprmekImh4LmaKYmH5pNcpArJMqZm4i6fPgPmp16dGHAMJjWSziAfm6oNWEl4S5pIsuwdnzs4L3Bqoa705LHAyUG6TBeTazmCvSHAwqPX49ADpIi2pnd7Ti9DDCVxy7BI3zppSeAE5itS8z62BwOmOPyPotbM53ljLMLnoquN814NxXBdm2dhQ75p0bqXKbXTzYq4DwDrx7H6T9BGwDwTCX0W06C2d49o6dFknARuDcCQtkmC1b/kgYxJlxRupqbf21Ul1IjntItnW/ZA6pdbp+wlf35alFjKCnRpA+e9p5UJpGLgqifgFgOeAdyczp+Ri9Zvh0JwsbDLkKQfuna/cK6vnQpPS5FLe+l65Sya9R58eXD80HpCllyGZ8tG3jh9DhvGFqXfSskl6DMIxaI9BYO+Mcin1JaHjfykrp0PynPytZddWnBmBZ5TVdNawaBwj+0HSXcfmiiIQ/qph3g4663hmaGzbahVLTfik2/aRmoxN33Tju2lZ83GDhD/dqQOeej3Ec+f8/KVvexfGMltNZU9mcYz4gQ/wSccvjcpUsJKL7BGXdTSvXgFVEjXpug07A+R6ElLXVDZHlvreq4dQvTOkHZDlUyeBSrZWlmi4VfJtsp6ira6hpHa1TaH2zPnoicwg/9a8PWtr5ozLHmy5UzrqqsnATt42B297SHpmd/IdECT8QPj7iq5LULuDefUwVrszCB9AGpWOiOxJKs6E2Ubtmz51JHzlAlg6mSTt+ylyVUCD6+z5oLZyjHlvvctY4wG3tlCphcv10tC9b62nxr5MTcLR5+VG3OqCpMdIaKZ+Vs94lmqq7b056hxS4dirHmPGMwlYVsssvgrkE94dGOhV4csn+3gf+nhrmFfdre4b1G92XjXTJxeQm+dLU6MwOwqhSJ8dU+RSpWqwDO5ke+GxFp341kPS24myAfCSl1mGumBOo+E3ScyTBh1+IPw9o4NO3viejKRNUkdnCb8e7K5niHAMkaQOIPUOeWiQVKfzv7fdT0WWRS6NotukdqEcPV0fTRuYOtgR0nSStTxHmSS4tZ4qCT3PK5tHgreSolmwsyfihonZV9yllexbJHz1DlKzNKNvhqSb5y6PrM9IkeGtYW6s8BwjqxA6wvl2khkF5uvDpNUN3dWf5DaG9XXlvdPvG8dsraeNKqvJOlOJ2STtg4eUtXoSEq12Iy8puw6xQ+w8x7dxEUXzNawYAuGvIprEWI+OHm0jaYMOfxsu1v7O1fadB00vDVvCsz8jSZbUGGPIcAKXVJi+lvgV6ctLG+33w+5FisxIdUz7fMzNa39rOiX9z4nKWIWUpsujnJaNOnzjvDNmo/p3ptTuzQ7JT2LVo52spyir/DxjEJhcLDvZNLzJ1ipVEwBqo7TvGvW1au8Yaa2YpF4pra0hbj5hupXefHyEr5w8ow7o9xHX1822z5xBefGi3r/XG3jPf7c/VO+myZbje+me9yKF0DEa3tVWi0fVrpaYK4RA+CuGTrpFuzPyXOO21LGxgff7MW4NKqnJspK+cqqHT69dwcmAJGuXkJ4+XxjBTeVUsuMVyZtBnqWpSqYB6pOIrXtp87PuqkLtEkBmtCE8qyDWwMyVxYzZyJeSQEu05F5qg5+0UCqOu/1h5e+vbA03L1butnRjs+DpXM69eQqy6E1KbQy/9PbK8L6jjv/1x1Tq5Le/HfHp8wXeg/N4BwYIsIOTddU5duBBf7sAVVoLN6un/r4qv8j7iPHM7QnNfunWfbettqREV512yIi9CYHwVwhFgY77n13LVX/BO6NtpPMQqmQh/MY2nWJZIqabfm+L2wllopQ6RYKdLpc2W79K0bKc9IrcI5U1DKwuam3fBNnEvfwROLws2hve1UTMJ418fvdS4zO28iBvIsc+sksYE7HvYdHG3EspeI0XjQeQGK+bhu9vRmdbK2wh1MVZvHEFJAikoiqG43nGHjfXNttX22pLj0FLALidZDP75mFAIPwVgZRuGmPyhriVFP5ByTq9LSW2Hm9J1Ep9Uhf7oMFLUbS+9MOc4HkOHSm5psklhM5L54Zn1OXwtuO80rpwg8l2RaIzZqPduJca9yT8Bd73lX0oUIlnpOMGS89EQGqRW0mBtxMVcDW50OA6Oscm49gJuGpyUPBW3fKsfKlYum+1RUXZ7QAv4/wNq8/DgkD4qwJZSxW20dSIkm2QZLhqwKty8Awo8sO/sabOd3NY6JQJJI1zPa+P7BHrxGhcaK+5vGHp3FUqY7cwrw6/iXuNY2W7J1FnMp1nNmrbn9G0vRJxCtvsB+yZUUq9QpF54TxcMvbqwu9J4pZ17LjdHgq8Ncy1M0Bd5jJDOcq9z7jIq7QUpaefkzDElnb5lbI2TFurLdeIXcdakOC1m+6xKgiEv0qQ7lLSdvmzRVTqrAUUfr7yJPCSaao7/B2IVRQlETiXmKrBTm54PpdJuhyfCaHmhQb2ndPtpuvhXbl37oabMO9stAs0qRq2LVXDvsKjAvFdg67xK/w5dpq2OxDjyyctL51YZR410kdLVSe4SSVGCfZsid12UsA4RnllhJlQBeq5Gxhfbal+W/d1ilXo3J9WGIHwVwi0bDckuNTK1OgxQvGqTUSwmuzpC3IpqyQvKSivTK47tUZR6OPLKQVeqWIaBRQ6gKle3UvLy0O5WNJ3ra6Ic0q8XQTmXXHvLiRvA3udNGbBsxLZLx1+x9NbAVeK5CeWyk48YfrVN233+wO8e1Z55sg4NjJrGiUgpaxXnFDoVBlcLaN+y6/LrVRG270zA2b8FeZzk2Yf5t3Vrvq1L7V3l4BA+CuCJh2tzlFjV63y/B6gThGgOyT5YOe5HkC3kwyv9QtDitGV/6S5IrCLgtztq8yK8UCqJqskZZQtksh+Emc6spQTwkxD5T5g0dzbiL1OGrOwtBurT8eD83hNX5kXfltIwybTFGWV1C9LSyXHlFIJGiR5M3sCtU8BWMYKFhHlpQ0dCwIgUfxPr6J8sNkraBbZcwGGau+a6auzvdXeXRIC4a8KZD1YaBD4MhD6dLZ8yV0Xz8hMCb8iirKsUwyQTUunUKj2UUrcvujqYLcGqS4DqL0WrCLcXJdL0lerl87iHmfr/qHFkm7MlHzrFZ5WA/IVBwvU8m3bkNT6cMFiHLj60TYOs3gEmnC0GpGpd7YhwTGk7WkdADTBI6KeSG0vHZkXVTGdzBhTXIg5bJJ+IPxVAgvTlxKdrJQ+sndiraYeP2pLmrZre5MEQ1WaVKd2dbFNOlxfaUAy/qoGJf9zhNj3eKCzmowk8hkG2y0rL46xaPF521hCjzamxjwdj1vlzEnqZ/XPTMhaAGEEzgWWIlclJZ34iYEwo3wPSZ8OhL9qYMv1VgMu1it8O21APjLtAHZnLAo+JutVAamLmgxvpD/1eWkciPdIwNLQWZtUlrULZ5qa5D8cOuUceffU5/B527B+ZQeYj0b1KrVVqgfAF15T5/WnFXM5lU5NYVJBakcl65oMlVCLam3VVpqB8FcRjOybXDRtKYu7QXol/Gr9Sy6Udh4UX7ZEdX5Lek9TZ0KwU/jyySPg6KAzeZGkz/XwaapWsCXZdwrdXcz6MC2RvlCrErkqkvqwTdq6L8IQJ+sC750ZGH365rDQtmE71z+RPak4vfmoZmRyW7LJxYtA+CsIKbsFYbkeBR6jKB0Qx1iMSm+SMvqtPUg+cUZVM6LSfc+fHOpyddseHT4NgIP0HglYUVgGURJSqA4ycaVR6H7KM4qaggZXG5K6UOZkQJ5qgeXbD6/hyyfOGn353pkBPvPDI5TTkvX7+hp8VcvGUBWmFxI/dT4zhS1LuPH19c5qsANGIPwVRdc0C5Z2xYnm5OtfWqZSRC3PkGhL+J84I7CAkT5GFboo8fmTyiNHXxuNYp7asqMos2rL3YDFoiZ89BK+lMhInNJ112UpyX402Si0JxhCHUvy5VNx3YfPnlWNV9k47/ViFJul7u9c0CE3zes9oYUbBMDnTynPs2t9dR2m8GWRvqfzugKZS/YHPQYC4a8wZnWGxswEtucAy5pZSzAu2dvqI/J0sD/X1Zz4hc3Zc1dxuRuweNiLQk7+rpoGkRM06dLr37v+9i+fMgO5EFQAFsWb1C7DhUP45G5MtQa2YYjxQKoEcQNh/MZxamgQ2x2BjB2yjDEQCP+QoiXFudP3aLXAex4NoukUtQTD9aEUlGVL/oaP/y6xqsvdgIOBTYKkw+eG2Om0ucwDXyn4SlL6vIQofQj1Z6rORf39YyddNc0YhHKAqDonBXXpMWHr8C3f/jYJf1ljIBD+otEm+bZ853xVWh1JmL70TFWvAqiQdSzhqmx4ZasCcrx+RuDoCuU5L/UyOt30G7P22iG7LHcDjh58790OAOTlL6k/22SoJo0GZ4NNK82yNXNQ3x8MzFWFTfZ8n9ShAMqN+W4/qT11aNmSJJ76B82EvowxEAh/kWhbs3ly4tB33MkBUek1ySeY2vElNYtjxJ0d81T5qNnjh6JhXzgRs05d4hceUvuTgetfb7ir8ULlM+Cb29qWuwFHDz7SuzkskLtollN10GSj4M5lxtCxhRivv/1w2DS7IC/YQhMOBVfRZhtmuaDDdfjaVdrOu43dVDYHPQYC4S8KbVO8lZid50rn1Z6SBP2RfVl9rJ2X3pb4ez3EZy4V1rFVJCHLQeKrdyo31dL4xlqG8tIlXTxlK86wuFKtuTc2Zj4KX8fnGXeDhH98wH3tyUWTJGrby6ycVsIQ6zecW29dHKGs+uQdiDHbnKI8dcoYCHLq6vm5tM4FoE+dUzp8nxeOnY321jDH+2tmu8byGtG4dt/+kZTwAeAdAPAVAPg6ALzf8/3PAMBfAMBz1fZPZrV5KAgfsbsSz+pUbtyT65tcl5urP7aKWdUFumMqcFJ/l24qj4PRFVe6QQAsRqWO/JU7TKHa66GsI13MnOMzHgHduseTM+jwjwnIE4eXS5xYNXS3YkX2WabkCVtYMJzCRiOdStkhXtbZdNBgta88clRZTlLRTAbqMzLIvnxqzcgHxcduY53ozLVx+frzkdThA8CDAPANAHgTAHwXADwPAG+2jvkZAPjIPO0eGsJHbF+zWd/Z6YeNn1mRfk11Vu3NNxb0OVJZVa9yDV9aoqHr5VY1eznR4RHY13cE6kEHzAk7xfPtJMN4IHV+JtrykTRsVADV71hWWJ24DdEtikJ9krK+DodKaq+qht3tD/EenFd6+FIaKpq7/QQLyLWH2lacKdWl3TkbxkMxKjt73hw5Lx0AeBsAfIbt/wIA/IJ1zNElfI8k4C1mUm1cl2kSvivhbyd1UJaU/v4HoLwd7LEgBM96ONU6fNr0ErjqjXqOahpY3R6FM++1LXcDjhYMCVa0SMhACdCETqOQCeU6SRlZtesxzytiSdmOJwP7e78/YFKPkvx1zWdrPNr5oJyb8QhKhr4/M4aS97m07e83Fk347wKAj7P9n7bJvSL8PweALwDAbwLA9zW09R4AeBYAnn3sscd2d7cH8XTZdM2zX3Ly305MHT7X03PSHw7RCDCZxEq/z9stcmmYBGzSt9U8WaYMY7avPRm+aF9Wy2ItdYz2V8IPqpvjB7MfWBJARbzcSUCmagnIdepjEEaCM+qTuvwhZ1hP/IluWy1tWzdjcrFhZyCsrkXm+Ur39VUg/DMA8Jrq/58FgGdmtbsrCf8g1k/sHGTYuQOxKtnGPAtuJYVT3Jvn2CZjJvWnm8MCJ3GGyVAaElJdn1NNDrz/bW6aOnxX4Gly1SyrPCM5+w0zfJk5lXelww/6+kOOPQhOaqXnqc2bCq3Tp2I7Zt/0ZGQloYSVQdQSf+Xt5hSnt5e6vpWtIa2Xbu77GZ1allaxlBXq40tX6VjHPwgAL89qd27CPwjW8bRpu0NqspbsN2hlCcTa1mSuVq284bI+1jbekm6cbKvkOGAHZ5VTczCIVGK6aRaZ0Oe8tGGSOzW2Sy+doK8/pNjDy+RCD42LTEgjz73W1VtETZ87knpF9nahc/tczu9S4RrKAHSCNSJ/W7056zlwIew4SvgPAcAfAcAbmdH2B6xjXs/+/wkA+Oysdncl4R+EXsFzDtt1sul0TUKTo/u2c+VUEwsnaN6WLXzTJOKrdTsGoT0mDN/osqGxOXX4bfsBhwB7EJz4obeSwvAyy4SqwqYTA/psXx6Ctouc2PWfqQ06TqsrB+bvjDbWB3i3P/TyhGN3su6XG5VXdTV7EG6ZPwoAX628dT5QffZLAPDO6v9/BgB/WE0GEwD4W7Pa3LXRts1jZr9gncNX/Hg6NX/SxJvu/OFmw/SuIrL2WzOKTLBKQpz4uZ+y0d4q9NqA5WEPgpMhFEtWQrCo93n73uyUwtLht43nyq2YUoQAlDr+hGfUlELg1npq9H9jEshz7S6tg2kb7nc3C6CDFIaOT+DVkiR8I8Vx9VWvV5N+k2akSZgqIHd0jYadIOt4a7xnWpOUTlZVRfg6nhFBF3O8sQfBqRO5+ViTOd9Lid4SoHaBICmxypCJuk+LzbLyUsvx7tmLeL9vtoEAeB/MBGwvn+xXPvp1zWYpssZx4Pjht+SeOmh15/Eg/CXp8G0vnenUNKTyBFG0z+F0hrK5OIRdGMWXRtnZJ7L3uoeWRoQvLwW3MuvTgIPHQQhOdB57n/o3082Tasinw6eaEkTyZCwuQKVX9ua1t7Y7MGDjQtRk33TfczD4QdCSjeNB+IgH7qVD57AlAV8WQF6bcyZRT133SGl5ynTug4zESZLnJM9Jf6GDO+BwYBkM5UFR1HYAPglsJ1kdmluWtcA1GKDMc2YsrqJuy9KJ0rWle/Jas72GvPe7i+dzUPMn4fgQPuLBKMs6nMMMhpKsELPlxWP/3GPMmtVLZvbByghc5KbahrwNivwA7B4Bhwcr4nKlV6h83+eH3zcJnNw7bw3z1uAve7XrxA00jYNdMPhBmBYJx4vwVwDcNbKOdjX962nw+MbWzY5Jmzhm9sHqH2e/nL/zBhwDHKSVcV74Oruh+hSayG9d9OThAcDnTw6R572fDIST78dnLzB2OjJ4kPCPMOGXJeLaWsXRA4ml5WFTMt2gz8UrGZq6Sl4DlH7ThLmliBVZvgcEzA27s3tIP0tLncNnDJlT6epjpwRefEuJk4FAWUVC3u0nmAxr2xyNA8f7yF457Hb1fcA6/IcgoDMQAaKoeR9A7T/4oPr/h98ewWUcwxsB4L1wFd4LVwE+BIAiA/jgGKIHIhiP1bFXr6oNIIJ3wGnYijOYPDKGTz8ewXPPjWErBrj0yGm4/L4ITp8GePJJ99ouXzY/u3wZYDx2r9G42NOnAbKsPpAu6PTplh8GBCwRvs4+GAA8/7zejQBg/NpfhCj6HHx2KOCP3zqGz38Q4EMPvVUf8T/+vx+C7y8juPTsGOBv/xAAAPzxWgI/8iMAyZNjgPcBwOnTgBDBSy9V4xMRxnAZog9dhacggz8WY7WvBq8z4FZuiDXNBMveVk3Cn0etKaVdkMeURpKhNLJI2vnKyNuAa3MolaxPMtizFLHKy/eAAA67c+c5oi9nznBYB3EJYayUVQxK6aQMJ7dMylulixFZpyU1rUjNuAIjNQO6Y7Tpu/0GBJXO3uCoXjyqGOc3muPd0G+qvZllisjtos424XfR/a2InS0gYPGw40uI2CmBPiN6O3KXUoKr8WsKYpzs7bEmS2lpkUpDaMtH0ik9uqzxFwh/D/AlQCM/3zay58fqjjbKtQGJXCLHkOEENvB6T+nrbV9im/DbJIOuUkQQ6AMOPYjsEVHXC22qlch1+5UNLL/iSvjcuSIeyLq5vMDJQOjEhoavv/DnGVymKSwQ/i7ByxJilhkph4082miqZwwJIS9wMlDkrkscVilW7fTFVPqQG23jgZu2YS8dKKwEAo4E7HrRVgETnXzQ4x4jRyMdyX69p8YmjUFF5LmR5I28d7YhMYKy+ARBZH9QnjhtCIS/C3Ap3Q5MohdNL5SSleW5WlVevKjUiqMrEodDxP5Z9bcYubm77YhXLYmwxFO0VGyVGjqI7cEpJ+BIwOq4tmBGmThtTxv9m4qZXzgRI3fNvAMxFjAyhDpK5WwvtXldatpsW9yyxlMg/F2Ckz5/k9yAmiS1zYjP8AUUlfSgOsVgXUkGN4e52XEa9IncFcy3UjWwy1DvZUsiAQG7hqcj2/WiKVLXqGhFkpkni2z6hNLLc9fomuCFM25twrdX47om7wEjEH5H+IRkVaEncxiSZ9Zzl3Om6iceuNWnjLamzekNWgV30mNSx+9YLfwgo/4CAhYGqyPz7K9JogQlEqiKohKcpKwybLqEf3NYYFlWOvmB2faWFZR1dy3ReXfStB7z1/pFtW8OyYNEIPwO8AnJIpV4Yy3TxG2nLibLvSloUKcrnYmC6+qNCFdesaSrnqXJU2GG2B4k/IAjAU/wU13NLcNffbgwPG6oWlw+kg7Zc9IXqTTsaDbp2zUmZKo+u9HPDJ4QAjFNZ9jGFuQ9EQh/BnwcS++0AFV60FemkEC6O/LPrXNzm9XFJ3DJSKWszfksGZRxQU29pe2CW8T2oMMPOBJgHZdSkNgFULhBdTBgateBxMm6n/AnA4HJxRI/eT7TghlKqaNwtxNhCFhyaJdorM9JzhvkKOQMLl2laP+9JwLhd4BP8k2SaknGyhJmWRWQUb2o2iWrdNQ2f/mgWV1cxjGKJ6Zo6BenrJwgf+Ez2NebA6eD2B68dAKOBKrCJ2LTXUnb3jNaxx5zwczSyadKpQOA+KsP10JeWSKmmxJv9AVONoraG5T+8aqU6nMJoepVbyeZ6cq3m1V9RwTC7wifbtu76qpYk0v913sZPnMpx4+dyrwETJLIJFadlDwIdpOH3qiPy86xPeyuw2/bDwg4FJCKkPtnzXFw5lF/0CJXu9qETxXg6qBdZQeoE3HWapokMdOO08adNIj0Vf1odZxOxmarcvdZtxoIvwO8um1hGU1tI6nHJazcsXyzSIKYljiJ62ANWUqctOShbwuaanIX3YqF6X8cxPaAIwwp0VvEfAwZrl+QjpazJvvMIGjuYFFOlQu1/bteD9nnZT12swyLUcnSNpikr3nEl2xtQX6cgfBnwCcUUyUdTfqSpTVusnzyEDtr204yFE9MrU5ovnBHdVQ0XG8pnQ62FQvzBubpPEHkDziMkKowOh8HRN6UvkQIN9UOT1muqtOZkbPcy8beCsiVumdU6eCrQhdylKvVelN0vK0+8HFFkPAPDk3pT0n3ZpcytNUpWLKgKh5nDaDTKdyBmNXfdOt1mnk+2t+/zM08+7uuRRuU+gGHFD7BjCT+p8/7vXRqKV2aPMvGnS8XG5G9HrNE9laqFE/6HrNkItf38EHuGfS75f5A+B3BOa/J7YvCrcldU0vnIsNbwxxvJxlOLqnZXk7VJCCv1KHc5DpGHedeL8Zyx10StpI9U+vYnXbuGw5uOwGHGEVhql7JnZoW4pRmR0rE6Y6t5pE69w2i+uuSvZk/h8asTd7FlalB9vS/Xd/C0eFXJyfVMHn/7UXuCoTfEdRBtFbEkuJVqgNT52fn1hhcqIspC4EoNssqOMtN1nSvp0K7swx1nh5aEs4m+33i6OCYH3DIwbsqJ3j6TgjEZy4VuqYz2dBIBcO7e56rVcBgYKp+Tp1CHL5l6gZOAtRqnbxwwmNoFWJ46dCYq1JrcgFuP8Z0IPwOoI6iZ+jUNapKoV4+z4fNZ/9r/cKY3esOoyaB4eOmzl5sqs/rPDlMYvcZjNm17qsWxueeFBBwyNAsDNU2L1Kb0oHaS65CWaKOnKVV/Cd6KrbGkewb1DNek9gMO9l+yl2B8GeAP2wh0Ii080nxipjd5aFphzHTKYgnXOlgDBmmm3U7ZPflekkdDm4x+r7ZWYOEH3CE0NidG2o3U7Q8olKrJEldO6XcKXEyMH8jT5xoJn3SDe3h2vdD7gqE3wF2R+HLubKsjUG3kkIvE+137npauYZZqdwCtCF3O1ESBq3wuAvX7URNNjzr39x6+q43HXT4AUcEjcTp+UKvlnMVbEWr7IuPS5wMBBaQ+wne3kajPV9zkPAPGHZ/oJdvkLFEg/ApI16ThV9sevzyq5lBu3miVTTB47dL+Xf23YEmeOkEHCHMI+HXKchrlQ/55m+DSpvg1dn7tj0ISPstdwXC7wCf1G6XOyunph5d17REM1qaUjI0JWEycuY0XoM5+3R119z1zbftBwQcAnTR4duMyknfO1ahSo3sc86PY22wtQfmvENqP+WuQPgzwDuKJmvLvWqyUThvhEvoiGY+JFnWUjqlVZil7yPDsa/zdXHXDAg47mgiTt/4JUatV/amkLUNF13yj2PES5fqGBuROXmw2si7bSLYL7krEH4H+Nyp9Eviy8EZay7jJeV1DVvMlF8+t+jzyNpaujcNxvb+vurwAwKOIBqJ0/MFXwXYQpbtlXO/H+OtYa6aqVSytxLliskl+yaq0MkYF6w9DYTfEa0zbINycCYBV2I/GWap9qbMlf8vz4zMUzKbOTlmF07vdhMBAQEEn8rnY6cyBCgbXTDHIIyoXp+K1UcVtnq4SWbcj+G7cMIHgHcAwFcA4OsA8H7P968BgN+ovv8cALxhVpurUOLQRpH79eqzZmheEpFsATzYjhuGebY+OlWSYDcdfjDABgTMBe6lg5lKniZSidtDf0KdrYEwxmYTWfs8hWZ54uzX8F0o4QPAgwDwDQB4EwB8FwA8DwBvto75eQD4aPX/TwHAb8xqd9UInycs0y+/cuNqk7p5XVrfy/bnWzNTrJKxuLUDBBfLgIBdgciYq2VkXriGWiFQ5oVD5IjdC9A1uYzu5/BdNOG/DQA+w/Z/AQB+wTrmMwDwtur/hwDgWwAQtbW7UoTPnj6RPOn7dORe+888Pvqmsw59NqvC1kx1TgiiCjhGsLu27QvRpet7JWshVYUry1uH1DhNQpudPM3nANI0PPdr+C6a8N8FAB9n+z8NAB+xjvkSAHwv2/8GAPQ8bb0HAJ4FgGcfe+yx+e5y0WBFT7henZc61LCMuHYnod+bncV0DSO/e8edbBaaRIiAgCMGW6rO86pube46XjTBK1kLN+EZT8NMEfBcHetbqZOBliT+qlLiTB3+XofvoSF8vq2UhF+B69GJtO2q9KQLtN03eZ5smdeJnFSnqaV4r+/+PGQfJPyAYwCf9wvZvwYDP8F2aYs2KkGqPXlEnS+LkztNMIjIhMFqnLOVv5SzdfSHRcI/+iod9M/I1KH0TO6RyH05eHgip41L5m9kKTETu5jmgw4/4JjBpyt3It1F9+HjGFktwm5S5yCaQhwJg1wVy9vx7R8mHf5DAPBHAPBGZrT9AeuYJyyj7SdntbtqhI/YbJgxvWhc4y5fBtIxvkyc1LC3HNqulZFZ8NIJOLLIc3Oo2FtXsifhzda986HTaHAt3WycTnbODjgUXjqqffhRAPhqpar5QPXZLwHAO6v/TwDApyq3zN8HgDfNanMVCR/R7ECkM7RLpuVXTOtsMjQLNAyHqjACPyYeSMxHpnF43uLmzgX69gMCjgi4Hr1p66LOaTKy2nr4VpfK3O/Fx1Mvd8F+DN8QeLUA1B3Ajowt8YUTZi+cDISe5fMcMb4w9aZOiAfSWBpmGe6+dGFAwBEGJ2Af6aepG+jka6QozNxXvA1KjNhV3WIXTFpWVHwg/AWhljBcQ+sdUNWs7q7Vin5ZSrx1cYQvQs+YICgjH9fzG500SOkBAQ4o95TtLk91a31qGf7jOnFaXXeWO1dwF8+Z6pYVcpgIhL8AuO/XnN0BSgSQKFMlZtw7n6j8G5U48iL0EGCKY1DfyzjGAnLURckDxwcEzERZtuvfjXHksZAS6dM4bBS60BXA9L6R2zxbusNEIPx9gv3eRqNmCZ86TjyQOAZTBJHWGnQMAuP1UrfTKWdOQMAxB+dV0rXb+xq2iO4xAIxBpU3YedWsd2tkM2ftSMkugpJirYDDRBvhPwQBnfDkkwAvvQQwHgNEEUBZAvzWbwE89xzCVnwZNp67Cv+ql8HXvvUI/P0HPw3vLa8CAMDl5z8I/0svUo6oFX5Qfh7+AB7U+ydfA/DcFyK43rsM//BbV+F34wz+zlUEgEifLyAgwEQUAZw+DZBl9bgcj9V3p0+zcYOoBu9VNSZhPAZ43/sAnnvO2+7/9j2X4Wf/41WAOINnvgfh8uUInn8eIB4gjOEliD50FRABLsMYfupzl+Gtn7uqLuKDH6xPyi+GA3G5A7ppJlj2tkoSvr1CI+8cMhiVI2Vofex7S/zoCXXgHVAqGloq2vp9x8rELPv5aAGVrQICjig6ebZ4dOz2SptvHz2R4c6rUquL6FCRSrfOLblU21iSmzQElc78sN+fL3TaLHgjcW1N/b3esw4EULr8stTFUGQcY5aWzoTAUy7s+aKDTijgmMPR4Vvj8naSYbzuH4e218/mJu2b7WRCuhy+xEDINsJ/YHlri9XFk08CXL6s3iaA+vu+9wE88oh53HPPATz0kFopDgYR/ORPAgBE8I++ZS7lngIB74WnAKMH4Cv/2TvhT3sxwH/3Thg/FUE8MNvcii/D+IMID/A3QxfStN900Zcvq88DAo44fEPEGBKIgO+9bB4Ux/C2H3kE7mxcBgCzgW14m/5MCIAPfxjgIx8BAEAYg9nOGz50GV76awQp2Yek0skyRRAPPKD+cv3TMtA0Eyx7W5aE3zYxzwryEAJxeNE14E4GAotcsjz4ZVWonDwDhGn45UahLsvCkFYh4BiD3DP5EElTlkvHGGuZsS/7a3qcjiE1JH0al7zOdFM1OvL28Ur6nCQOYCxCUOnMhu1+ZefosDPj2eRve+tQh6AUCvfOJ7g1EIz0JW5DgtuQIM+/Manyb1Duju0qS9/MFHsr5AccEHBQ4GkRiPRp7A6HlXsmFLgNSV2tqlQzghwmePfs0Cu9vfCaAcpR7vDABDZ0jI0at8rV+mvnN5Q6VtRR9V4iOYAxGQh/BnxCtP2eqEA5L0UYx8o1k/OsW56wxMlA4N2+6pVbA4H5lZJNDAI3n5BsApFmOTThya3Di6DzAspLkCYCApYJ31jl46icch97JXDpyKw0xfxKiVsD0236+hk1Rnnb588jbv58PW7VBJIabd/tD1EOWXiukWxrNXT4Syf2pu2gCN/3HnwdSCibq5YokqSaBITE6dQuTag6282hkixuvG4T5Sh36mSSJw/3+LGFgbJEr7EJs8wsii6EDvIKEn7AcYI9ZqmIUDmVeCspcAyp1zNOCoFZWuL9fux8TuOG1EXlqMDrvQwBXAPvGIQZa0NkwZceRB7BS2f5Er5PG1IJAI2JldJUETrlwAFAXOvXObMBJE4G7Md23lbdWTKkpGs7O+bXly75JfxXTvXMA+MYZRVfvp1kQYcfcKxgEn6tVp0MhNctGgFQ9nrez/lgp1w4Vz4gtbv1Rx7KcHCh9P5mK/YIXUT2fBW+YATC7wBbiLar1TiCc+nmujcMOGmJxahs9fW9AwOsdYGuhN8747bplVSmZSXJZLUOcYmRfgEBBwX/itxTRKjLJoT2vXz5u887Kc0/9GB7mwASs9RTx/SAha9A+DPgk/Btqz//jhtlbOlbigwxz+tgDLuQrSPhKynkl08WttCOAIjX+qZNINucetpQ3xsGI7rwgICjCJbPRmlOpGW09ahB9XjZdD6TvZ4Kqqlmj7t95UxhOmmUVjum7n9rkOL20L+iaCT7BcTOBMJvwSyPRl/AFS+CYlenytKyngR8qfwsCf/Zh4YGadPXZAxWBR5qA7C9YuCZNpeVjjUg4EBheVkUee3dhqhW33aBId+YadpkKlCk0krKJr3qoTGkSuc/ZN4+pMaxJf0Z97Ffq/JA+DPQ9NzJM8c3GVD+bJ8BR5F+w0zfIOX3Hi1NCX8gcWeHTzZMvdProdyZ6k79IvRwApeCuj7g6GMeCc2qZrINCd6BASIAfvtEsw4/GUr9cxp725DosQpQ6v1tSNR4L5k3R5uU2PU+9jCQA+F3QNPKqnESzqVKeVxJBCjdVMj8hX/rb6xpcrY72LWzI0PCVxk2yfhbTz7XzuZaOpnEGQKU+OVTat8omhJIP+AooynmhKRoGrRlqfLepwInA6FJmsYgec3ZThBUp4J73tFvJwMKknRVP7KUrs6+TYe/oNiZQPh7hHcykIzwK+UhSfX3zieOJen6o2mjx8AXHoxZR5M6OdOnzmeW4Qhx7WyJN/qZ8ftJSLgWcNxgG9ZoOU6DtRo0Mi+qIieyWpFnxu9IRSp7PSxgpL9/4QTVp1CHKgcKc/xuD4WuZW2Mu3lUNQuInQmEvwCQ/t6n0tlaT43ZWm6m1XemPh+gxGssu6ZB4gOlk+f9t55DSqedXSdcCwg4bPBJxjwU3lKPFDkPZjQJVo5GmvSbK9BVSRKFn5y9HN3FGBsk/MND+Ih+o+2tYY6SimQyv85754Y62tZeOqqiJ7nxnUhrsjdWqJA7k8O9XoxylC/zUQQEHAw8um7tJOGLXJTS8Nyxx+BkIHByYVOrbMyxWTtReCPeZ5FzG+kHHf7hI3yvWyaLfOUv0ZYieM4d8YRb1JwCqHg/iNdr6YPyeWjy5+kWAgKOMpjKZGND2bzIHbpJPVLkErdYIKQUwti30ysAyDojAlvJa3frWRXSu6h1gpfOISJ82RJ4ZUkaUmR4rZ/rqFzlqSm1ZP+JM2Y72qWs6gw8WVujhJ8HCT/g6MIRlktzXMTrpTcmhrsqy7xwPHcQAHE4dLzqJnH9W5kXeLef4HZirtoxSZr18l2l9+CHfzgIX0p0oluVJJDhraEraZAkQm75aar2H35Y5f643suMTJpbzK8Y0RZeXB2+EGh07rJEQ+LX0bcBAYcMbYJwWdZebSQ0DS7U5L8VZ1iM2MqXCJuPT0ta1xNHpUc14mps3+ymMbUg/XwXBMLfT/BlYoHKB5d3xLRUkoDdoSpJnXfecip1ucR8JA0VTpEzyaSUrO+4huIbaxk+fa6O7s1zxGv9XFXXynM9GW0nWXDjCThU6CIsK9tWoVfMkzjDcqesCpHnpvDkI2JmcysKpq8vCnPszEveC/DA6YJA+PsFj6hBHUPvNyXnJoUgVn9YWzojZiW2GCu+vKh87qVRuGEbEsOQRGoeKTJD11+XUsz090HSD9grFqCJaD1Xk7Bcloj9Pn1u2sYAFOnrfs/8m/U4YNJ6LVjVaUr05NLgobOri14wAuHvB7qIGlI2V2OoCF+rVpo6njWhkGRBOkU5VO1/+eRALVVLFdV7tz90PA1sXX8g+4D9wDJqc/uE5bI0k9AKgbh+wfW1305qtSsWhVmdSi+xlQ1MlrJaGRQ1T8/robNAD5wuCIS/X+gya3OSt1Q6xkDxePjYOkFZShSbpdOBKTxcDtVSdPtis5+/M0paOttBSm0BhxPL4LK2YbexoUjfHHLm7JBuSuM3IpUukVs3U7tlNh/TesPLmBUrBMLfT3TRy3mO8Q4Ue5lopei8dz7Bu/0ExaZJ3E+fy528274iD/ehb+zLVFQ13wrnkpfYPwMOGQ5SW9FVh19nMnEFpOu9TJM+bVrD6rkZSlNCH91KCnN17BscPmlpSRJUIPz9QlcJv+EY8yt/9j0ife4mdu/MwDjmXi9GueOmSUYAfAVOVeS/rj97+dQabsPQOYfvkp1BJaTZT4PYH4ALsEe2kOMsYUTvM5flG2sZpk+U+MKJar9vkrgxBKyb4TmpOo2DFZOWFkb4APAoAPxbAPha9fe1DceVAPBctX26S9srR/hddfh2QIYdFSgV2ZO+XZdT86RRlkI4ZM+NsT7Cp+35k0N85cSj7nfMeNx0e3NJNgHHDvsu4VuEKUuzn9HQsuUq2sjQSvErFJRIkj6VEvUuphskfO2DP6vLL1lf78MiCf9XAOD91f/vB4BfbjjulXnbXjnCR3Rncsvgo422nFRZj+Gdkwh/ayCUtO5JpyxLiTcv5sZn8YWpEW07WTcnilfgZOtEYHdAexDVh+5SdxlwpOGTaWYFnc5skPWrIpfaswazTHvObGw0C9F6WOa5t16Ej+wBENNNs4+XU8/NSDn7ng5Sx9UBiyT8rwDA66v/Xw8AX2k47mgQPqLb43yEaAdk+HT4pRnqjQAoB2ZnlULgjb55zGSQYQFXDH39ZD3VhtyZG+vIfP7iA1cfKub0Tgg48igKw21d95umoNNOaJGy6WOdG0247pJSslWBp8/beXFEKvHcOcS1NdTRt+VU4rV+gZ86z+xc86xo913HtXsskvBfYv9HfN86bgoAzwLAZwHgx1vae0913LOPPfbYop/L3jHnzO6o+jx6eCquzIssTwbKR5jSJqsOPHXcMCfrAuWZM26n7/dVrU5idCFQbqZ4Kyn0/NRUrH1u/+Om59S2H3Ao4JNpmmScudsuXT06H1JlWUe263KfwoqDKUsVbGhJ+KputLrGm0PVRu9MXXtCCoHXzuYsVkX4BbouD2YFBKM9ET4A/DYAfMmz/ZhN8ADw1w1tnK/+vgkA/hgA/pNZ511ZCd/GnDO7/tpXFYdIP44x21SE/vzJpI7IrULInz5fKB/kdfPclIa5cSPSv3hRHT9MdK5wTvZaNyoaIgyrlA2G0Qsbbn3FDFqHHkuePOflti6Xa6hxPFI56dpvV32RvqPV53aSmSkQPBK+IniJk2pVfXuYVeNJ6GN8+Xicvtv2QPZ7Ftwllq7SsX7zLwDgXbOOW3nCNy1G883spPun43d2jDYApqq/bJZYTnl6V0X6Zan0j7b72d1+gvLsWfN6LlxAtD9j7C6nlJ65MOcr2aLDj2MsRqVWTWGWqUITmcXjK2jQOtRYkcmzq4zT5XIp0AlBqXHKqcTrvcwg9ixT+eh5zhxO5tqDJs+9Ev61szmORmqFsJ0Ix6V5G1TVKh3xzm+O8lG1PeuF6Ll2j0US/q9aRttf8RzzWgB4TfV/r/LoefOstlea8IvC1HvQ/xRlO0vk4YrJnR3EnlliTS1Dp7owc1FUydZ0X647PkXgbsV1J5br63X7AE77CIDY76Oc1qHnfHDpS7dHLEtPSL7KfLB6Syyu2HL30GJFJs+ur3Oey5V5Yfm+K9InIYJ35XhgEjIJQPwEUmTGGKG0CkYSNGPSUGULdZplfnNV2nFH+LFvdl8t2XvDIgn/DAD8TkXivw0Aj1afPw4AH6/+/yEA+CIAPF/9fXeXtleW8HnPpVndJv5Zs7pn1MgTJxChrrdJfyeDTEcG8o0ybKabyvgqS4l315K6nq49Kn2kX222gczoo3Zn9bTtnSzs++0iEga0Y8mT51xzjnTjTtr4UpZmHy+n0jjnxobfieB2pc5BxDptQloagtGtpNBqIa83XCowXq/dOG8PMxSbpePxQzUqZj6cJQs2IfBqP2B53bguLVm9nOvaHifddbfmLS01G3gaB+vmIJIlO79Nsi0b+RxTemZHLWNft9W2oWv13eeKDIQjgSVPnp20SuwgdblSqw3ty7UXy7TZns1cxVjXn1A/kiSFl1InCiSSVzaq2jPOPhGpdyYDVbPik+drASu5aEa4O8FXNlZEsAmEv1f4erndQ+d5uT4SBHDcNCcDgXyZax6ujFBrayztvT0KG2eKgUPApN555lIx+74Tfzm4VnVO0OHvHSsyebYaYi3VCtWIoH7CSVNKNFavQrj7/Fij/kRe4NZA1HmlqrQh984n+tx8lVEUVSNW36XKVwUUmKZqUqByiE7B87bkgyvybhAD4e8NPpLyLAs7v9yu7THCp1zfnPy3IcFtSPDi42YGTpnnZvtm/li9ffs1yn3zXi82/pJ6x/BrztRnPBpYXdt0tg5/RQyNhx6HafKUruqFyN+43EoFkwzrfj2GDG8OC0czKiXWuepZ39waCNweeqR29jz4+DBmEj0xCX1tdB0Ty7Y1U4e/Qu8mEP5e0SQxz+uvS/AZQy2pm6T9ybrA6z1hdD6u+rl+RgWN6AEQC5y8fVSTPUk0Z84gjkZ6X/b7eLc/NM6pUzATJ+ee+x4OUaYCk8en+OVTMcpLl7SXjm+Qmv807Ad0w0FNnvvxvjzqDUfabjCy3k6UbUrfZiV1UyyUVis2rWJ5bWd+7S3Pz75c8tdv9Pm3sUKCTSD8/YDdI+z15rwv17YJWEvN7UQtNWVe4OhK7apWq1EE/toJU7K5u5bgNgzV76alK92nKeK5c3pfVv74tInUjQgup+Z9y4GKETBy+kynWhe7Av39aGPRk2cbcXU9d4t6w1b/2CuBMWS41pf1EGMFgJIEzTw7TXaqwUC5Pvs6oeceGuU5ezx00eF3eT4LRiD8vcLXIzzpE/bctr3UTFMscqUm+cA/NTs3gMT1C+ZnRj58IZzCzPZ2rzdwfm94U7SEq2uyL8tVW9EG7BZtLzJJ6gjU6livxDtHZ5ASnUhuM42xqf8XqTRdK+3+zb1qeJL8lk7ouzy7htFhE14C4e8FB8FmPqkqSVAmdSTshx8wO/eHHxSOV88YBEpP1k18/PFW4pdxnYSNfJJ9Ony7wAqPuF0hm9WxwMKEyQbhhuo0a2Moj3K19eWWl04jY0o3knsSK3XOUGsbXeOpQfZJ4johzNkJm3wT2hYGq4xA+HvFQejn7NWCoAEm8PqZmtjHIAzivXFWGIOC9P2OJN5C+OKJaRW5K3AbEiygUNJUXmiy52oc3q6OuJW8CIUxHwTsMxbeHS1ViSyl4W1T90XT68a4Lp/O3j4H0+HziNtnBnbsicceUAVrFbmn41nHdr3ltv3DhED4+4GD7hEeSYukb/LS+eLDyk+fdO3Xe8IsdNJxexF6CDDVK4a7/QTHoDwXyqkiewoEkwNVFJ3I/15PJacaDs0VNa2wKXP0Mh7hUcTCF5wNSzUiff657WI573VNNgq8sZZpXXk5lfiJnllPtsk9kla+3oyux3yZGQj/sMKStEyffInJUHkvUMAVV+eolUC7Dh+HQ3zlpJldU26mOsEUJaWiZbfs9RCnUzWYU5WZUF7acGq38DQQlZp/lZwYDg2aJsiFqc/aWNtjE7L90ue5rnoRWxto+T71cx5s1T9r+fSnZbNO346AJ71j2wM+IgiEvyqYp8M1eDBQGtd4UA2UypAl80Ilh4LEyA0yBoH34PV4H9bMkTgYYDnKMb5gpmi+3lP+9eSXbwzu6dS8dKbD980nJOEHo+786FLWj5Oj1pXjHnmtQaEtK0tm7RqcNZJ+k0bFvo6ybAxB0RtVXROpyhhrR+1ONoq6QlVRoBwO8V4vxsmlvJ5VkqS9gsoRQyD8VcAMH2AObjDdTtTylVLDkj7+dvU5Sd868ElIFi2IqEq9qR3ysyeR/HpPuIZYAF0H1NjK0ssk5GnBpTIiIp3Uim5XmNHCM93cjilmTZBtueEnG8XueM1mZk7ipTSjXKt3aeepqa+z7gO0ArX1+TQcRiPVRw1ZZF3iaGS6Q06nfPXgrgrIq4ySCOpJQKeZjY+NxBEIf9loGcFE6PYAnWwUtacDqWxEhpjnTnpXzDLtvkljVXG6Gkh1pO4UhVC60q1Y6AIqZdpB/VMliqO8O3Q9Xz23oV3mVLCKaquAQiWgIte9okApzFQRrYEsxwEtorhPPSJElbtOmBI25Y+hSR2gnI/XuDAiWV9jvveXLtWpuaWsEqcOJG5sIMq80Hr1AgqcDAROqnQFFLV6c1gYEjtN/r98IrfKENbRtpSnzw4etzRN1irBza2vZ8lj4kYWCH8V4BnBtFTlfc9wHa6WqfZgbMrjowsvlxI//90beAdiRbqlRPGE8rL5/Hdv6NQJtZRYouMJwSV7dr6tqvoW3cuXT8WKaAbmpPHFh4d1cFZa638pIVyTSoA/rrb93b6C/W5z1+hg1LDVI6SKLgqlyptYxLadZDonEjd8tvGa4X6bJChT5oKZpmqSz4tKQJZuglh2HRQsSCfmfcLWuUuR4eBCadRn5oXHedS373HxctKOGqm0P2DPuEnfdIQQCH9VYHW4IpftmQJ9IlpTHh+SYDY2ENMU7/crsh0MENMUZV/p8L/9cF+J/1JW3myldsf0En4lVkoh8G7fPGarkiZt4uHb/bVYkdDAvW7brY+wCAPvntvcz9mig1HDJ+Fr0tcqO48fbPUjp2KUB0ZuGqtfbcPQyK8EUOJ1y4OGq1Qcqdp6zwDS7dPVZrv73lir3C0bXsEsCd/pj0HCD4S/ELSRgqejk97d5nAjiLdp5FOwCXeH4X/TFKWVNE2uD/Tv6kyG9cntRGq4tmaQfrljEgzXnTYNduN+LdLPhNxLoOZcr2VPbS5iBvK9Vw/Z2/N8rZ4rG0nWLCZSq4L4qfkK79YwbyyPSZK3XSjHueyWPmDnwnfr15r9iufRaXt8zrhhqw2tw6cHGXT4gfD3FbNykDQwzsQanA7hU1v2AR6JZStWS3rvBFFNAgahV9vdfoI3+gILyHFrIBSxC4E63DFJsBiVzu/GkGF+peV81jbLrY/f7n4LY7tucxEzEG+bXxBry+ssk9aqMC0Ve4LqAErtDktpAnh6GbvU3x1Y7/T+7ImENp68z7dN4lp1KVKJk4F5rC3hj0Fl0Gx9tFJaj0/tGwFZ/P0FL51A+PuGLqTgGcFS1MYszuUGlzRZ70ifb0lSs6StrUHqSFQAJYpUokil0t/SJOURo144YepayaNnayDwNktTKzdTf4Rvxwyji1C37rrNg56ByCbDVnn1q/CkGrCe8/WeImava3rK7SdiduF7q39Van2DaHViP2HWi50MBN7oq/0bfVUk/OmH1bEfeUjtf/Ehde1feDDG5GKJN/rq+xtrWe0gYKMoDMOvXdnK+R0ned/nRwyB8A8CXUjB6mCylHaSTLNKYttE4lH+k67V52qpiXgw0HlzaLvXU3p2b85vZ1Kpy8ddO5HhBDY0wVzrF1imArcvplrff38tRjnKa/1/h4xUTY+S3Dyd1c8+vZ623+7rDNT2Xu1i2JVgcCsp6sN9E7qlp26SxgFwdnQqqER633zQVAn+2kkW/yEQNzcRT52qy22WU1ITCdyKBZajQk8wpPu3aztcO+t66VAt56Znx4Ox+L002YSOGwLhHxTmIAWbu337iDi76pRFGDrlAZxHORgooy27pv/4Nx7WgwOgdJbTTo5jp+AzDdba1xpAapNC/6zE4cW6aIvYLGvDXpXu2XgIFrimKo4Rr1ypBdhr/QJvD7Na9y9N104mHBv7bc96Fulrw6Zt1EzE3kjf53YyGJipGtl71u67vgympL+hm84UYTbxuU+Pbm8vvGag//+Tc0MdfU1Bf9x5S5VLricXfZ1Yq6LM5u190w+fS+jeRyzdpGvaLbTY/Ss5KgiEfxDYhQjZ2RboW4paP6bkUxPYqI16AyXJyzjGW8NcS/5kiFOXaUbaGtJlpfOX586rUnAsoyb52vOfubwh9YDWkiliY1Y18rcejWqSp9rrvTNuqTwdjCaUzpbUDaSRGg5rj6eiUN+lqWlWyUfNy3wpFXmRB5MUopJY1c3KJGlWO8x4hTwojUutXret6p14q5B5Zi5ekL5ZwndfmE+XP4ZU51TiKw3exadWF7In3VnRtPbWKS2xR7g67pI9IRD+otFFh9/y07b9medFLhXXidSyTC2fyZCrxgctr9XyuYDclfABFEvmufbSIWPtHYjx7tmhVs/4PDbsbQIqHkBOVTyAitiJlRGN3Swn+SxDfPXVmuz5BOLL2ChSnk63VjfYBDK8KI19Ocr9Bj6+Yiglu18VOEbSrhwmnaRKe2IfXZFGYrlyKo2oaWfjEyTNXPYqjF1zk9bP0eFXrrrkzUX1YflmqEikbSht93ach+x5+nr72p0htBf93DFAIPyDwJKyg9FpSTq2oxOL3E5/QH/rdMf/7kSM+ZXSUf8Q2QNMdYAVJ9r1C+1kz1VGL59aw8mFTby/VrUTx4qVi0IH0QCg1hFzkrf3bUJquwZql+uOhUCU0xLvnVH364ThZ1bkcmm6ryKoYDKnRqsH9tz/zKUCP3pCXUccI053pPZvt7NRGjdBpE+qH5otLMLnfUJKNGI9igK1O+Z2kqlEeFbUthyY75l7Uvl4lt6XlnPSUu+XJep6tfaKjTa7IJs9QbSS/TFwsdwNAuEfFPYkru/udLyvc2mL9qVEIw8KhcSPIcMJXHIkfJ+U59vSTal18+0EXTq5eSSRPVRSckUS1/pusXZu8GtKldu0yiDvVV4HmNxgKY00kb7ekgTltPRIySYZG7liLBWG/Y7q91RfP5E+r2PgeMykzNMpjhF3dsx9nmDGIjspUTO/jlhl+n1DQrAenC9BGlcT0amoX5FOH/Ncqw+LXGpd+81hgflIanVdv4948aJ6dvlIqeKSpJ6Q+OV4h9CShKvDgkD4RxgzV7fMwEWpHLgu3Nbh+yJifdskVqohU6ViEzRt0+a2Kle+p8/l1nWZ1xmvl44On3tn2IY/7kKopHt/viA5HDrXg2mqi7LTpOVTJenVAjNQct6xs2I0TVqqUll1MjI8VKutrXXXvVWvkhrFYGyXhLm+xGr7ReiheGKqhIcqR8/NYaFXC9xdmJ7/7SqNNrUl+2uGuC5ToVM0IKp5Sko0CuzQSqpzWogDFq4OEwLhH3F47FfOAbYb3mSgdPt2EIwm4akbaGVHXZKahyRAinI0pXSPJxAn14oUnhlkjcTaO6OyJ95KCp0lNM9RJxH76rkNnAz4OSUzKkt9jzfOCut+Bu41nT2LNPGVU4lc9XU7yfDmxdruQRMNuRE2rbbIWMyvj58ToFQ+5ESS1SR0r6fcFbcubBrHZ5vWBNpgBPdKA+T2OZ3WkdQ0kVTpOO72E5R5rgSEzVJJ3zlbHbT0K7vP3OjXvvhE7qRqtA3Q1K8oj47Xay1gJgLhH2F0tl9ZswL3uyeC1t4om7V0dgdiQ/VDvv5bcYZylHtz9tcum2ZyLJu49LWIjOnuzevsnal1xNOpug/it6JQUmip70Pg4EKJX3y4NrL2X1ei2PRPOo46x7iPTYwvTNn1D7B3eqdWc8QxfsyzKonj2vOnznkj9SojfaLEj54wn9kLJ2IsYFSrvuoZots2y4hgT7L011rdUAGdbRjqNMOkAtNpHHyxJdbk5btGm8S115bVeXkFLDJPkLonoBsC4R9RNK/apTkuPa4UUmSVq2Y9wCi3zq2k0CHqtaqkIvlKkjU8W4wBL3VOkwIKnMAG/sVDa3jvtRcaCYt0+GRbsCcQsj1kWRUhzzxHSPXwh9/tJn+7AwPchqE/MVwLocqHH3bauf5oituQ4LW+b5Kr0z4rzyiBcqR05OUOJadT5PryqTWkCRAB8C8fNK2YOgjOMqaqc1iSPdPh+wqReKUBn9unRfpUPtCJ6K06FVcX2u3bEj7vF975yeo/rdljg5TfCYHwjzAc+5WdS4TpVm1d7naSqZS13PZF6pKsnjgct0/SIHj87iaDOhXyJFYqge3Ha8l+eygQL5jkvz0UKDZr9ZLPaEgeMcqwK6xCHCnePWtJq5bHkZ1ITuubh8OZvoM8cvn2UDiGTl5r+KMnqrY8Lq12u0p6brFvGOfYdG0Q1cvYTqoEaNZ71G6e1nv3+eDzCaE2nLq6Qt3fmHS+FWdYjOp+dgdiJ9qbT4ptE4aZMiGQ/W4QCP+IwzuAbEOdJxoTi6LR9tXqCGG6nvj96Sri0JKmnUOCHK+rz8kIPBlkLDLTNAJnQqIccp/4upqXQ+ie7fZQmDnZBwMdTIaW/rxx890rI7S2QjKOa+fUtVnw7ddes9noRSUr/0aZCsw2p0YuGXotVGyEv8TbQzfNtXFNVWxBe/Hwqm/oVWC1ohzleK8X69XMGAROBsKJ23DUOc6ElCFfEQSynw+B8I8TmpT6tmGvwyjyTgaGk3dRqwhUwVF0Rin/y78nn9HqeicbBQ6H6ORq4QO/nLorCmcSafjuDsQYX5jW6YArlVYBhTJSCmG6PXq2LSsH0SfPCfxw5KqGPvl6v7rIJnyS/slbylbvcPsH/5yMq5SKoZ4EMyQVnJ6s2XvPRxKvn1HXsDUQRtF7vk+5cBDqnDykw3clcHsVWBp1aGnl4zNsk5cO7yNN0byB9LtjYYQPAD8JAH8IABIAHm857h0A8BUA+DoAvL9L24Hw94CZbjt7aNeWyhrc+7xRYE2jWEqzqYGqYVo3yQy3r5Y60pW2u/2ho8IxJNOqIXJnpNUBAJppBqrjXjl5xtvO9tA8r5qcSsfbx1sTmEi+MhTrnEe9GMtXpzrj5Ivgkr6t9pmsp3jjrDmpyFRoN1mn2Eie1ys0KPD6GbUaIUPy9TOVy6RUaiCy4fi8amReOK77hhGdXmv1TpOkSm/R4Lpq2x3sYLYQUzU/Fkn4/wUA/OcAsNVE+ADwIAB8AwDeBADfBQDPA8CbZ7UdCH+X6Oy2s4/tNxVhsf/6osOsAW8X5z51UuIz68oXvIACP3rCX3jd2IbNE4Dsr+E2DLXkqSNsbZUTbYNBo5rndpJh79HSq8YhlYuWxjlhDwQO3zLFLz6cKBuClChHeRUrYHoyFTBqrkbGNuUzL6oJrYpTSNlEVpaGYZ4Hs5VTabxfQ0VY/8EiV/EW9ty+tmYG/trdpUlt6EOIqdo7Fq7SmUH4bwOAz7D9XwCAX5jVZiD8XaBNh7/fpM9HvJ3XwZ4MGvM+FE67tquornE6EHj9zAyiP3UKcX1dnWdnx5eMR2/bMMQxpHVeHDvKlYy5m5sor4xqdVGVjY2u8z70q+s1Sf+vT51TxuW0dAu/VIXo6XMplK6bdN9846uHOzBotBFsQ1KpZaxzkW2DErBZz7dr15DSdGyKY3e/KRxgN92rbT+gHcsm/HcBwMfZ/k8DwEcajn0PADwLAM8+9thjC34sRwT2aCAyXZSINMtGYE8GTZK9hx1qCdRsvyanEu/Dmtu+vbW4HirSNKXoTzyaOmqY2xdTnf7hbj/BydtH2hBNF/sn54b4SvSwQ8w0EWSbUyyn0ole1gFrjPT5duNs6hB7Y1UqoeoPGPduPwvLBZWTve3ANevV+xY7PBFcwPKxJ8IHgN8GgC95th9jx+wL4fMtSPgd0LT+tUffInX4TeGlnA1s8vdcj9GUMCcNbrgtwDIOp6naZnnXsO36o/6Vgm00rQlSuYFSPhgp0ciiaf/uFTiJCIBlWhs/r/cyzEdSRyM3RR/LalJ5+lxuXYN5zfd6gzothBAon9j0q6Xo+bDPKPcNd+DqSti++dyYu4N4vnS0Ef4DMAOI+HcR8b/0bP9m1m8r/BkAfB/b/97qs4C9ABHgpZcArl4FuHxZ7V++rPZfflntE6Jof84ZRQCnTwNkGcB4rPbHY7X/yCMA73ufOn+WAUip/j73nNkG/a6paYEwhsvmT+AyACAAILwDPuO/rne8A0AI97s4dj76h3/1Ye/tnSn/AzwFZhtPgYA/Fk/BB8cRvPztSD9uAIDPQgIAAI+W3zJ+8zS8G/4ABnDn158H+T2vhX91RsA/+tYYXv52BJeeHcNWnMHnz73Tew0RAGAp4R2vftr4/L1gXvP3fut5+N3vvBUwFQDPPw/RAxHAcOg2+OEPq43hzsZleCBC/ere+U6AX/xF7+UYQAR473vdz6n7wZNPsh2o++STT85uPOBg0DQTzLNBu4T/EAD8EQC8EWqj7Q/MajNI+B2waANt23l9+/aKgwd9dbk+pmN2E72ZxtqJ5VboJMFn21ZMPvK2asMNDvJ9pj1MdKoE2kqnXCRCbbAlN0kdvUwLMU8kLSUYU9K7emZt3j4IldcQz5i51qLusnT4PGq2i1DOvaioOWM/ddsOLjbLASzQS+cnAOBPAeA7APAiVMZZADgHADfZcT8KAF8F5a3zgS5tB8LviEW5YO7leujvbgzIbNLIc8S1fh18RS6FOhkXbUliqi0sNhpD6jWe0ncq9UGzJ8xkIHTeIDMoLG0M+JK9HgJMjWhhYlcixm1IjCyeslIBTWADt+IM8yv+vDR6giC7jG9itQ3WDWWkZnnF8NdUFLUtm5oyct0sSwAJMLAwwl/kFgi/A1Z9gO3Wx66SPOnWdK6fTBGtSD2TXFG4BcCFwHvnhjrQShGvwHi91AQ/GZBPfonffGgNSaofrEvtwbMNiVHTdVZErTupZDoVBHc9pdqvZKi+lRS6Wpkvc6iXwAl24BtArZy3E9Iwd8u2Odln/6cYO/aqzO62agLIMUQg/KOIg3LB3Cv2YMTzzhfCUwqQExlvv9ovS7Ox0agOuJpsFHpB8MylHG8PBcYDtbqQpSJ9SrNsE/DLJ5V0P1lPcbLePAHYwURFgUYCOLovLbBPTVUWX5lssVxFxn37jOXTKTtB5p1o25yu5u5eqy6AHBMEwj+qOAZRKrb0OK+e2Bf9qR2ZmKTL3S15gXGzApMpvd4c5jhZT/HpU4qAJwOBcnMTv/2wqUundAT8tXjnQXZ9FAULIKtUBQK3BmqCMt5z14m/hXSbhPK5+PuwCCDHAIHwjzKOmxvcHJNcZw5qaJPnfmkqoA5Q6nq5YrO2E9yB2FIlZd53w6+B2whoBUJ5ZTLhRsTu5pm0PSMfqc+loTkGAshhQCD8gKOFOSa5mVJqy6ygM0EKM0LVLa8oNVET+eu00y2k50s1fGMtw2RoeivZq4O9PhPfs+kaVjFTWD9uAsgKIhB+wLHGTCmVrJGWUdSo4VpJ4HZ5RcpJY/5cmudokOwNsmUl/mx2XSRntsXuBQ3N4UQg/IBji056aErzbHvB2PUCpBmhSoZXXxaHrqkKzGtbjodLk1AeNDSHE4HwA44lOunwmxgbKl93xoY+IuQ/r/KqzSUJ16uP5pKC/Nz2tbR+uA8TRtDQHD60Ef5DCw7kDQhYGpoyQQCoz6MIAKVKZ+DD7/1fAJewOg5VloDTp1WmAMQ6Q8RrX6uyOjz1lP8cTaA2AVQ6iffCVfjdOINLnx9D9L4qTQYAPPnIGF56OdL3YF+LOuhJlWqj9aD5YV//fmXpCFgSmmaCZW9Bwg/YL8xUWeSFI+Vvxcr/vqv+usnNsukgp027XGDpegq1rlKCwj2gAgSVTkCACYMThfTmttlTMe0OCnDnEKouRYcwdc5MO8RcTvMBRxmB8AMCPJASjfz7tsvlrotpzyFxd9WRd/KHn8tpPuCooo3wZ6ZHDgg4qogigPFTEbwEp+EpyOAyjGH8VATRU2NAkcFnPncauIafZ/6d3fBYGQ+uXgV44IE6bbSVHrqLjhyxTsnceC2dDgo49miaCZa9BQk/YNEwtSB1vpuyxP0ppr0PEndnT6Ogww+oAMFLJyDABAnEteAd6X0AgEdOR83ePYBg+PYguqJ5k8TdUACmCV08jQA6HRQQABGu6JLv8ccfx2effXbZlxFwhDHLk9HmcUSA6Bdn/IgONGcTd39OEvZei91Ep4MCjjqiKPo8Ij7u+y5I+AHHFjapk2DM9zkiwLqsJIBL5NRYN7F8LnTyhw9O8wEzECT8gIB5wKV3QpPUHiTugCWgTcIPhB8QMC8QlecNQcpA5AErgzbCD26ZAQHzILg/BhxiBMIPCOgK2xgrZe1rH0g/4BAgGG0DArpiAcbYgICDRNDhBwTMi2CMDVhhBB1+QMB+Irg/BhxSBMIPCAgIOCYIhB8QEBBwTBAIPyAgIOCYIBB+QEBAwDHBynrpRFH0FwDw75d9HfuIHgB8a9kXcYA4Tvcb7vXo4jDe799ExNf5vlhZwj9qiKLo2SZXqaOI43S/4V6PLo7a/QaVTkBAQMAxQSD8gICAgGOCQPgHh48t+wIOGMfpfsO9Hl0cqfsNOvyAgICAY4Ig4QcEBAQcEwTCDwgICDgmCIS/IERR9JNRFP1hFEUyiqJGt64oit4RRdFXoij6ehRF7z/Ia9xPRFH0aBRF/zaKoq9Vf1/bcFwZRdFz1fbpg77OvWDWu4qi6DVRFP1G9f3noih6wxIuc1/Q4V5/Joqiv2Dv8p8s4zr3A1EU/XoURd+MouhLDd9HURR9qHoWX4ii6AcP+hr3C4HwF4cvAcB/DwC/13RAFEUPAsA/B4C/BwBvBoB/EEXRmw/m8vYd7weA30HE7weA36n2ffj/EDGutnce3OXtDR3f1bsB4K8R8T8FgDEA/PLBXuX+YI5++RvsXX78QC9yf/EvAOAdLd//PQD4/mp7DwBcO4BrWggC4S8IiPgCIn5lxmFDAPg6Iv4RIr4KAP8HAPzY4q9uIfgxAPiX1f//EgB+fHmXshB0eVf8GfwmAPw3UXQocycfpX45E4j4ewDwVy2H/BgA/O+o8FkAOB1F0esP5ur2F4Hwl4vzAHCP7f9p9dlhRB8R/7z6/z8AQL/huBNRFD0bRdFnoyj68YO5tH1Bl3elj0HEKQC8DABnDuTq9hdd++Xfr1QcvxlF0fcdzKUtBUdmnIYSh3tAFEW/DQBrnq8+gIj/5qCvZ9Fou1++g4gYRVGTv+/fRMQ/i6LoTQDwTBRFX0TEb+z3tQYsHL8FADcQ8TtRFP0sqJXNf73kawqYgUD4ewAi/t09NvFnAMAlo++tPltJtN1vFEUvRlH0ekT882q5+82GNv6s+vtHURRtAcB/BQCHgfC7vCs65k+jKHoIAB4BgL88mMvbV8y8V0Tk9/VxAPiVA7iuZeFQjdM2BJXOcvF/A8D3R1H0xiiKvgsAfgoADpXnCsOnAeAfV///YwBwVjhRFL02iqLXVP/3AOBvA8C/O7Ar3Bu6vCv+DN4FAM/g4YxsnHmvlg77nQDwwgFe30Hj0wDwP1TeOm8FgJeZ+vJwARHDtoANAH4ClK7vOwDwIgB8pvr8HADcZMf9KAB8FZSU+4FlX/ce7vcMKO+crwHAbwPAo9XnjwPAx6v/fwgAvggAz1d/373s657zHp13BQC/BADvrP4/AQCfAoCvA8DvA8Cbln3NC7zXfwYAf1i9ywkA/K1lX/Me7vUGAPw5AOxUY/bdAPBzAPBz1fcRKK+lb1T99vFlX/Nut5BaISAgIOCYIKh0AgICAo4JAuEHBAQEHBMEwg8ICAg4JgiEHxAQEHBMEAg/ICAg4JggEH5AQEDAMUEg/ICAgIBjgv8f6kDuCQT3zvYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Sample from the distribution!\n",
    "# just\n",
    "\n",
    "n_chains = 500\n",
    "n_steps = 1000\n",
    "with torch.no_grad():\n",
    "    # the vector x contains sample of the observed random variables at each time step for each markov chain\n",
    "    # at each time step you must do an inplace operation, i.e. x[timestep] = ...\n",
    "    x = torch.empty(n_steps, n_chains, machine.n_visible)\n",
    "    \n",
    "    # we don't plot the hidden random variables so we won't keep track of them,\n",
    "    # just use a temporary variable z\n",
    "    # we start with a random vector of binary variables\n",
    "    z = torch.bernoulli(torch.empty(n_chains, machine.n_hidden).fill_(0.5))\n",
    "    for timestep in range(n_steps):\n",
    "        # do one step of Gibbs sampling here and fill x!\n",
    "        \n",
    "        mu, sigma_squared = machine.p_x_given_z(z)\n",
    "        x[timestep] = torch.normal(mu, torch.sqrt(sigma_squared))\n",
    "        probs = machine.p_z_given_x(x[timestep])\n",
    "        z = torch.bernoulli(probs)\n",
    "\n",
    "# plot the samples at the last timestep of each markov chaine\n",
    "plot_samples(target_samples, x[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Single Markov Chain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now try to use only a single markov chain, and get several values of $x$ by fixing a $b$ and $s$ as seen in the class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAAsTAAALEwEAmpwYAABSeElEQVR4nO29bawk13km9tSQMCkSMEecFu8lR0tRyhq70WKn2+ZMN+lYc5PsAtEqAGXHMuJN4NiAFlrDvF01l/sjEsTbdSn/WNgL6HKutKHWEvYrE4wVGQjCJKSElXivFECjkYf0UB/xSqbM7FKiJJsrjrgbJNJ0nTc/Tr1V55w69dHf3bffByjcW91Vp05V13nOe57zvu8JiAgCgUAgOP44segKCAQCgWA+EMIXCASCNYEQvkAgEKwJhPAFAoFgTSCELxAIBGuCWxddgTK0Wi164IEHFl0NgUAgWCk8//zzrxHRW3zfLS3hP/DAA7h27dqiqyEQCAQrhSAI/k3ZdyLpCAQCwZpACF8gEAjWBEL4AoFAsCYQwhcIBII1gRC+QDBluOmpJF2VYFkghC8QTBF7e8DOTk7yRHp/b2+RtRIINITwBYIpgQi4cQO4eDEn/Z0dvX/jhlj6gsVjaf3wBYJVQxAA+/v6/4sX9QYAUaQ/D4LF1U0gAIBgWfPhnz17liTwSrCKIAJOGGNnpYTsBfNDEATPE9FZ33ci6QgEUwTLOCZMTV8gWCSE8AWCKcHU7KNIW/ZRZGv6AsEiIRq+QDAlBAFw8qSt2bOmf/KkyDqCxUM0fMFagsgmYHd/WctexHUEqwXR8AUCA3t7wIULtq/8hQvT8ZWfFwmLv79gHAjhC4413AGsUsBnPwscHOSkf+GC3v/sZyfT2ctIOI6r6zQqxN9fMC6E8AUrgXHSFZQRMOPgQLtPHhzo/V5v/GtXkfDTT+uOxqzDJJY4zw3whPCJE/lEsfj7CypBREu5PfjggyRYLyjl349joijK95XS+3FcXVYUEQH6b5Lk+2GoN02/+Wfu9Ue9tnlN3jqdvA5unXzXGwVK2deatDzB8QCAa1TCqwsn9rJNCH+9UEasg0GRJJuSpo+Ao4hoOCRqt+sJ33etumu7JGx2NGYdpkH2syhXsPoQwhcsNeqIdRLSTBL7PB/ZNyX9umuXHevWYZpkP4uRg2C1IYQvWHrUEes48oWvTB/Bs7zT6zWz2pUqHmd2Si4Js6wzTUt8HJlLsB4QwhcsDGW6fNmxPlIfR77wjRLKrHkm/aa6fK9XPD+KiLa27Holyew1/Kp9wXpCCF+wEEw64enKOaOSJl/fJwm58k2dRMPXMid7uQy3YzHBcxDjWuLuaMI3umhSRtW+4HhBCF8wd4yiM9cdOwlpuh2GrwOpgq/TCkNt5TcdcYxLuHFsjybMazftMMxOj6/Nz1RwPFFF+JJLRzB1EI2WGz4A4eTJIP8OhP19fcDJk9pnncsE8rJdf3PzGN4/cWKy/DZ7e3Z6Y/PaZgrkKv/3qshbd9/8/PXXgatX9cYwYwbKzjXL4NiAL34ReP554LHH9H6no4954ony8wXHEGU9waI3sfBXE65FXOuhYpygFFnme50lbH4fx35dnS3hSazsSS38qrKqRipKEfX79nXKvInKYM4j8Mb7YViUoETuWX2gwsKXSFvB1GBalDs72jJ+8EH7GCtNsHECXbDDU+n1G8aBxeuYUbRKaWv44AB4+GGrmCzVgGsJN4lGde+HKE/DcPUqEIbNUyD7yqpLh/DEE/56jpJ588QJbdmbuH5d1/3qVeAXfzG/NtdJ8vEcY5T1BKNsAP4pgL8A8I2S7wMABwBeAvA1AL9QV6ZY+KuJsmjTUu1cKbrSs09QYURhX3m1arbkefKU/+92i5bwrAKcyrx06nT1UX363Whg93k2GbWUlWN+tr1t122UEUTZfdbVSzA7YNaTtgDOA/iFCsJ/D4BnU+J/CMDVujKF8FcXvmhT/twlRqWIwr59QthXXvnCJSWXyFz5o85ts2q/6n58njJN95vGE5TJObx1u7kkw/XxPVt+Xr5gszAkOneu+FmZi2oTSHzA4lFF+FORdIjoSwB+VHHIewH8y7Q+XwFwMgiCe6dx7eMIour9RYFpoWyfP3OX+HvsMXsi15QMAhCeDOwT3v6xHQCEMASefLI4YRpFWlbhCcwyuCmQGb70yGVShu9+3H2um1l+WdbMpssfBgFw991At2t/vrEB3HOP/v/gQEtmDz2kpSxXHgoC4K679ATtiy8C/b5dllLAT39avPbBwXhZN2kM2UowZ5T1BKNuAB5AuYX/vwP4JWP/CwDOeo77AIBrAK7df//9M+sBlxnLYCFZVmc6mZq5CPZVqYvgKK6YfALLOSqMCFC0D72/j4hUohpFvrpWvilj9Hq2W6f73STuok19992I21HKKpN1treLk7FlUsxg4B8RlW2TSDqjyFaC2QDz8MOfBuGb2zpKOuMQzDSuacLyeY9jUmFE/W2VeqVoQn6mG1vkwSRLVO8t41772V5M+4iofYZlHH2NGLFXXvCRYL9PtLGRk2qS2KTO/7M/uo/4yp7vpJk6zfJHiSfw/fZ1hD0cFsthycec83A9p7a3p0f4XPem0ppg+lgGwv8nAP6usf8tAPdWlbeOhE80HwvJJGczKMdKBRAqUqGuyD406R928n0gIVNLZmucySWObWIrq0OSkEH2ejt1t/ISkPlsmMiZyPr98gyYvmfqkue09H7zGF/5o5RV5hJaRvidjh1QtbWVd4DcEZ85Q3T6dPE8c9+XU2iUOQqx8BeLZSD8/xL2pO1X68pbV8Inmq2F5Es3wERRTPaVSyy8aet7kJK+SglI0eVN2/Lv923r3LwHi8iUot1d8xKqknzMc93Jyiqff/c7d0J0UqvWhI+Yzc7HPbauLPN/30Tuo4/mvxkTvNl582dnztj1cS17/s1ckq4b5fC7U/VOCenPDzMnfACXAXwfwE0A3wXwfgC/DeC30+8DAP8YwHcAfL1OzqE1JvxZWkiuTOALyikmGrN7HyAxdPZQW/7tfB9Q1O3aejUTwWBg1+GZrpaMOu2inFPV2flIs+q5+b4zic9HcpM8Yw7KcrNxnj49njsnH1sn6bRaRLu7+fH9PtHmpn3Mxobd+bgdpFunutQUu7v2b727q+vhHi9eOvPDXCz8aW/rSPjz0PCryA9wE435Lfz+owm9vNFzPg9pH6FF1q7/vek/HvbtCdpOW1FiSEhs6ZuSk3kPTZ+b67PvWvpcH19+mXGedxXhM/mO+9sOBnkZYUh086a/AzOv6cYnuG6cZR0kUXXyOZ9Fb37Pv/W4z1EwPoTwVwjz8NJxLTu3oQJFDd8k/8N2SMm2rS3sIyyQtbfcyCQkf4fSaauCxFTly8/PzbSeB4NcVjCDtcpSIriePOZ1RiWsMknHtIbN5zJKIjP2uPGRcJnV715vOKw3LHyjQbfjMH8Pn+EgWAyE8FcM40wSNj3X10BLLfFBTIcdTeC9HlEy1BO3MQZ0qVXUFy5vhF6yB/S59nVzGcc8MAqVfxI5sonUZ426WSWBfA7ATWtctsatOxk8biCS26lWEfSoozefzOKTe5jsNzdtqaXTIXr8cbuzcTvSOinMrLPbGYxzT4LpQQh/TVA3OnCtNraC3SG6aXHGA2VZzsnNhP7gDl3IYdtmmX2EFPYVnTlDdOed+VcxYjrqRJQMba3+me4gG0XwpsLIYgofUbIbqCnFNHFhrEqJUGaVj6Pvl8klw2G19DEK3N96OCxa8/wbmJOxrRbR+fN2Z+d6EJll+yx781mX3dMoz2wSA0dQhBD+qqNBi2iq/7tE4U6q1U2SKkX07fu26AV0aB+2rPODW0/T9qPKaviDXUVHqTvnpVZEpozzSquTkXwUGjl1HKZwrWXfiGR3t1mqhapHWSZ1jUv27u/gEuOk1rBJ1GUTuu12rtu7clZZnINv5OMbHflkN9dLpw7LEGh43CCEv8pIg5/MFqHCyGoRru5cRyauRef+rYRSpPo5C/Bkrbnf6+pRQb+vG30yVHSpZVeMpaGjTpRH1Hpae90kM5OoGWDlfuc+C9/zqLrOqGRcRmLnz1dr+ONatiaZu51ct5tb6j7y9r0jvg6k283jHMIwnwfgZR1dTb8J2Tc1UgSjQQh/VaHs1AMZ2QP6c6XylAehfwheRuaTWFYqUZVeOsNhThhMaINd23xWidb040G5ye0jgLL88HWZJbkcn6xjTuaWWbVN5kaq9l0ppExOG9eyjWNbuuGt18vJ3pdR1CfpcAfoGy2127pMfl5NvKiq0NRIETSHEP6KQintLePzZOHJTTfFgS+/ik+2Gdeyyq0/mxH0JKw9Yevz+DFbtUrqW7UbaOUjdrODcTdT7zeJ3b3vXs8mTNerp0oWa9pZNgl6m0TeMTtZ31/fCIAJn+tWlndnc3O6wWpVUpqQ/WQQwl9hMOmbLSIKlUU2VdatRbxRte7rs/Z89el1/Z3QubNFDx2L7Mcct7uSi08WefzxYmqFfr+oT1dZlCbxuTLZNDpL87nOwrKt8pt3OzR3BMDHu8/QjM51t1Enm90O083zM00Lv27EdZwhhL/KUEXr2PVkGQzKyb6sEbEbo0v47LFSarE6MlM8UFmOncubUcHKTxLymsS+eYg6mDII/2XPlLJ8775y6yzKJvvTIOsy3/ZJYHYobtmuZu8mlSubwyjLyz9KDEHZ6MaUmaYV9bzuE8FC+KsKVR78xKTvS4/g2+om5UZqdM5Ecu6fH5e7HRoFxTEVRilNG6Qvb4uP5Mus7mmR9aQyhNl5jUOgdXXz3WNdtk73nsoI35QOfW6ldfMXvnfWTXQ3LjlPOgI7DhDCX1EolacPZoJkTf/ZXlwacOPb6uScquNLK2dgsKsKk5BlUtKkDbLMiq3LUzMtMpi00zCfDRNmFYGOgrp7LCNn3z2Zm5kfyeygOHcPo8y6Nucpmhgkk2AWctkqQQh/TpiFbui1hkOVERl76fiG3W6UqalN8+RckwbYFK4FWeZ5MkqDLHumvjJcCcdX3qTD/Wl1GrO08Ee9R/MezGhcgOiOO/Jn+0u/ZNf35k193NZWsZyyzmZsA2NETDoCW2UI4c8Bs9QN6zoSn6Vv5kFn64yzVbKLnhuI42uAo3ZiTY9v0iCbWIuTTJw2vaem9Rr1t56Fhs8Y9R7d+REzgrbTIfrwh4uRvGbnUJVbp0yCq1zcfsJ7rzIoZmGYLROE8GeMacoUZftlx1ddm4OfTAu4Kn2Au1+VimASNLHw657pKKtHzQKTksYyyg4+Uo4iO9maj/TdlbbcznzSGIRRnvW47800MqUuC4Tw54BJGvCoFqN7vJkZ0j2/TK/PVqpSdjZJ85wyn/VJiGmUzvG4WmrTkoWmWR/zf9/Ii9+54dD+3keUTSaMm6T0MK87Ssc+zsiwSUbWVYEQ/pwwjm44auNvenxVI+a6mcewT7S7PwsrdJRGfFy12GVxHawLbDPrWJf8rYmGb6Ku4/aV585HTVKmeQ/L0vlOA0L4c8Ak5Djqub6GWeVz7rPwzeOryGdWhNvEOl9G2WOaWPQIxZX/yqQ9tuzrPIt4tGh2IKOOVt3jfe8Av7tVI9u6+y6TnI7DeyaEP2M0sUSaaPJNibUqaMqEL/VCVV543xzALKIgm2DZZI/jiipCdQnUXBSdKCf9ra1qoq8yXJqOVt36+TqcJu9HmRHhTqCv8vslhD8HmC+8TxOvm5RqamGUWexlVvtgYIfVsztmr1dvRbmd1rwJd1lkj+OOKmPD/Z1dCcWV/kZ9V6o6HP6+yn24aUZUrqsrMZkavlj4QvgjwbRszBe3alJylOP5nDrC92mnbhi9r1y30Tci3BnrEouWPY47piWb1ZVT9Tv6LHh3TsG3chnLMe5763tPTQPIJHt3EfbjMJIUwp8zmjSiSUYEW1vF3DHtdh4Aw3UYRecvq28l4YoJvtKYtmxWNlKomyNy3z3zXTXTVvviTdx24BuZVln2ZidwXF5jIfwFoG6YPK7mX5VOgRuF2Ym437ua6kiN3j1x0ZqPYGJMq88uMzB82rqbgtr3XVU5robfbpcbN3w/vshy0wDySVVN77tqfxEQwp8zmlj44wylzZGAj/BNrd632AWQf+8b7lY2et9BprP+cRA/1xSTkpZribtSjM8LpsrTxn2/y0YKfO7ubrksxNf1pR/hdjBup7esA1wh/DliFIu5ahRQVW7Zcn4myfvSBfNxZXVxr9fopspap2BtUEf4LKtUEbI5sq1aaMV9vcoWbHHPL0vxbHqujTJQnbYcNk0I4c8ZTXr+US18dzJ2km2sF9J38RK/zSYrWQmOF6rmjEz5pep9bNJxVF2zjPirVkRz5aJR2kjTNjxv2UcI38Scnn7VZXzyd1UOerMD8SXcapoieWIj3DduBrKFyPntP+pExbVqBccevtfD5y9fNigclfDdc3hz8/10OnoB+arYlVFG21X3XNV2+fhZyz5VhH8C64S9PWBnR/82gP67s6M/nyKIgCAo33/iCeDqVSAMgf39/PNeDzh5snjujRvAxYvAhQu6uiYuXKiuS7utr+Oew4+gMfhZmR995Sq+2Anxn17fx85jAeij+/hiJ8Lh9ZO48eNg9GsIVhae1wMHB8AttwDXrwOdDvD88/a7zeAmGQTAlSv6fT04AE6c0H/DUH/uO5fPMfHaa/q9TxJ93evXgTfeAJQqnh/H/rqbNDHKPbv0wm2XP9/Z0fs3bozRBqeBsp5g0dvULfw5iW6jhIq7Oma/b3sH8DDXN1zu96szYfZ69nqk/L1vUe5aVDw71Q8pCs1lDdXCNUzBfOF7PXz+8lXvvttmmlrbPlnFZ+Hz/JZvCmqWGv44zhmTAiLppJjR0/f9uKO+BL2enY6YG023W57ThoezPNx1/fmTxP7enAcYeUgZx/Zauipfl3bc4bDg+KBpEjbXIHKzZo6ip1e1N7ezaZoWwix3Wl46824fQvgmpvz0fS9wkzDtEjnca627OXDMMpPEHhW4UbTcAN1rj3WfYXHlLd/SdWLhrx9cg8Knw5vvLJG9sDmRf1nMOmu7SftrYsWP0kaq2pevnGNn4QN4N4BvAXgJwAc93/8WgL8EcD3d/l5dmatg4Te1MJq8BFWuZexiyS8yD1lNv/smEbrjouo+j1tYumB0lPnHDwblUeTTzE3fZIQ9rUSAo4wG5qQiFzBTwgdwC4DvAHgHgJ8B8CKAdzrH/BaAj49S7qI1/KY9vo+8qyz8upfSJXufv72rS/qsGDdScNKXq6yvPG5h6YLR0JRkfZp9mXwzybtbl8ZhksH9OAS+bF460yD8hwF8ztj/EIAPOccsnvCJGj/9UX8k90XyvdDmi+wbhvqCVlnOcfVI3yhiHiRc1mCmIRkJVhdVA+e6QXXZ3JRbftl1ffu+z6c1uB+nnHm3j1kT/vsAfMrY/w2X3FPC/z6ArwH4IwB/paSsDwC4BuDa/fffP5unUfP0R+3Fyyx8Nysf75veOkT+dAntdvkC41xmGeG7OUemNYxchBYpWB1UWc9VhoLPecGNrPUZK4uUViYdKcway0D4pwDclv7/9wE8V1fuIiNtm5JbU23bl0DKPZ+9aViv7/fz/9nLxi3bt5VlFZwW2YtWL3AxjoVflVStzlgZxyhzvYjGHfWuguGzcEnHOf4WAD+uK3fRqRWa9uJllkZT75WqF4jzhJjSD5O9OyowOwFf7pBJX8hlTRQlWCyaavg+YvZJj03z8dV1Mowqoh91rmBVDJ8qwg/09+MjCIJbAXwbwN8C8D0AfwzgvyGibxrH3EtE30///xUA/z0RPVRV7tmzZ+natWsT1W1cEOURcYwo0lGxRDoCkKGUjvbzRdbqYwlAkB8LKoQMlpVplsWIY+DHP9Z1eeIJ4PXX9ecnT+rovYMD/z1x/X3Rik3h1sXdF6wn9vb0u8fvF7efkyf192Xf7e353ymgvD2Y8LUbIC8/jvN2zO8/7/d6envySV22UsBjj+X1GudepxywPzaCIHieiM56vyzrCUbZALwHmvS/A+DD6WcfAfBI+v8/BPBNaA+eQwB/va7MRVn4Vb346dPlwRu+cp7txbSPiAAdiRqFebCS73rjTAKZFo07Kpmmhi8QVKFqamyUSctR5FRfbpyqxU/MeQKWS1leKssiO+q9TuP4SQEJvBoNPvmi3yfa3NRPrN0uLsDgpkSIQkX7iIgAUmFU2CelpjpE9DUUdxEIkV8Ey4ym7aFMNnLfezOwy/ye3ZnduTDTmJsWjp1b5qy2ZdDw3f0kIdrYsF8SXm3H60XAFr1xwpVeZBVe+kIM7Aqoof0mmimIfQ2j6QpaAsEyoSlButq8O19V5eTQ6RBtb9ufAdMn+0Vp/kL4U0LZMLJq2Nlk9rdw3iC2Cjw8P6BXWh1S6bJAKlF01InocCvOzpFJVcFxQVMJpMoLiA00H9n7OoFZkfAivHqE8KeEJPGvJGV60pgwc8T7fm3vj+6YAWqY0Cst/Xa+0uqQGiZ01NHfZ3nojVPdogSC44omFvTurt38hkM/2ddp+JO0raYef9OCEP4UkCS2hm/+gBsbxeFgPFAZMVMU0WBX0aVWvp8MlWWBWy+Bxyxg0ufNJXuBYB1Rl0rBtbdcOdac1DXX1m16jTqIhb+ihK+U9tLxEf6dd3ombSOiGDEddTS5awtCk/5z52NvYJb1AjlmgRraM09C9gKBhs/6dq19d+J2e9tO0BaGfrKfRIcXDX8JCX+UoVou6SjLSuBFP1zS1/p+fiyTvusVUHgJxMIXCCaG6aXjWtmue7W73KGJJlZ6GY+Il84SEf44P8Zz523f+mSopZt/9rbYGzXoegGUDS+9ZD+ihi8QCGy4BpWp4Xc6WuOvHW1TuScQUTGvv8sjy+SHv15r2hogGn29SVKEE2/cwAVcxD52ABD+z7M72Lp+EW9/8w1cvEhZWRcuFKNeH3wwK8ktOY+CDQIdtpeGBga3nMBL73wE3211cPp3HkFwywmcf16vHUt3nURwQkJdBYIynDhhNSfccgvwyCNAq6XXuv3d39V/Wy3grrt0tK3LAUoV1659+GHdxpXKeeTBB+3oXi7DjRJeaHR6WU+w6G2WFr4vgKNuQiU/1piMNeQVnoT1efCY7mFPbcSUhBF12izt6ICsZ3txceLW3K3ww/fdW9m+QLCOqGrzPo+dqlW43CRvTVe5mxcgkk4OXy56n8uUG6w0HBrJnhJ7fBcPlDcAxHTXHAyIOm1FX+5FRADtIyrsm0sITuPexBdfsA4Y1cjxEXSZxGp68LiTvGUpzE0eWQSE8FP4Zu59PfNgYOflPn9eLyt47hzRYDcnad5UGFHYV4UXwPXPTxJdiSue86NQTUTMi/IIEAgWiapsmNUj9fLNxwm+zJqugScW/pIRPpH/B3e9ZcwffHs7X0MWUPTcGX3QPiJqnVKUhPk+oJolLPPMAE3j5ViEz69AsCiYbdZc35mNLV/KE6J8zV2znbRaegLX1xk07Tg6HXtSeFFtb20Jv2yo5/KtO7vOeeh9PXeMOCN7U4P/5Om40Uo9s2bleUf1CQSLgOkZ40t3wutB8DE+7zlu6yZB+yJxfS6YrnHIRqHpAeTz658H1pLwx12YpKxTsDdFw6G972r+ZlnWBzPUXcTCF6wD6oKqeOv3c+ONA6s4gHJzMzf0mOTPn2/ePM0Oxyfl9vuLmztbO8Kv4lV3uOX7QX0+9O6WyzwjEuuMZlZFwxesE3yyjE+LN1d+63T0PJxJzG472d0trjJXJg2xkVeVUHHaGTiboIrwb12gR+jMEATa5xbQ/rC8clUUaV/bra18xRo+7uTJfAUb04d+exv4wz8EXntN7z/6KPDpT+v9Vgv4wQ+Af/AP8muYK+F4/W3dZX64Eum+e15pOZ57Nv2NffcmEBwHEOlV365frz7u+nX7GP6/39ft4eAgb+fMDZ/7nF4JC9A+9hcuAF/9KvDud+vrAnbT5c9MHBwAYQh89KPNVtGaK8p6gkVv09LwfXp2nfQSx+VeOnGsrYtWS39OlFsbrNmNa7RPw/gXP3zBOqDMtXJ7W7db05J3N5/u737W69lrX7BM0+vZEbQ8Omiy+ta8gHWTdIgm17N9fvjuPh/n/rDeIV16claGG1glkoxA0AhVrpVhqNtmlSTrS3HObZVXw3LlW/NzU/KJY/+qW4ucQ1s7wh+FPCe1iMt+8F7PCNRKFzRRSZoSeeA33WXSVSBoBp+Gz7lxfO7VZZ2Da6Wbmn/VOQyfgeeOHOaNKsI/lrl0fHr2Rz+q902tfm8vz6MD5Pq9qbfxd759IuD114s5cwCg2+UcG4QvPX0DuHgRX3pwBxcvEv6zp/1Je0zdnWHI+wKBALaGH0Vaa48ivf/GG1qL73T0fq+nNftOxy6j1wOefDJvc8wNdeBzGEEAvPnN+nzW7E089phf5/fdU9X+1FDWEyx6m4aGz70ra+OuDMPuVGxF875v8W+l/FF9YegfBnKODc6/wwuYV5nuYuELBM3gm+8yvWnYHVMp3Q5dy90XBV/nmec7zzzfdRVtKslO23EP6ybpEPlzYfPEqjmE8/3I7lCvLIqPOwh36Ogf4pXMIKcQDV8gGA2uX7wb+Mi+92Z7r/Ot73aLGv8dd5Rr+CbGIe5ZtPsqwj9+bplEIARZylIQ4aP7Ab74RdNNixBFQSafHBwQAD1OC0PbXSsM+Rj9P3/PCEPgK18pViMM83SrAKXplI1qXthB8GSu14hbpUAwGrhNEGl19OAgbzc7O3qf3S2btCulgJ/8BHjxRS35EAEvvQT86EfAQw/p7epVLeH42mONx3XpPZS5kM9Ezi3rCRa9jWXhG12sUkRRqKWUGHHWO3NqBJX4j3EtflPayT/PV63iVa3CvrKO62+r7LvLm/rDw05EgKLDNL3ylV7ktfSr9gUCQRF1cmhZu+K/5iQw58RhS5+VAV80/TTrXyEAjASshaTjGQspJ7GZqaUfdnSGS97njJeuvGNq9maHwaTf6yp6eaNHqtezOporvYgOt3R+nWe6+pywrzuBXldf95luXOkxJBAImmNU0vTJvm4EPZP/rOs9zbm79SB8Iu+TM8m508mXJDSPUaG2+NkiT7k7m3DN3bUSu8PYTmgfxpAgFfe4o2FXTO4EzB+011UFzXEKGRYEgrXEqKTp2odlgVyzdquct4a/cGIv28aetHW6+ShUhYx4ydBvChxuxXRorBOrEi2/fPJ0TDFiurypSf7Q6TCunOsXtSDn13L9c81sfjJBKxCMj3FJsyqAaxJLe1RZVrx0xiV8zy+owohIqXxx4kGJKeDzo0r3r/RCOmqHWXlq6GFvO30m0e6uVx80h4quO6eQvUAwHsYlTVcGctuoK+u6545aj7q5hKqym2I9CL9Cw8/2kxpTwCR9h4XV0POduTni33+4o0XPdgdWfm1eZMGcHJrWRI1AsO4YlTR9Fn6rpW0315W7SYdSN9LIIu9nLOOuB+ETNftVSo5Rgzjfd1mYyxkM7O+ShK6c3S6Q/xu3a/J/AR2K+gnt7vong8TCFwgWAx8ZW7Jvkn8+imRUNpdQISA08iYaBTMnfADvBvAtAC8B+KDn+9sAfDr9/iqAB+rKnEjDr9r3fBYPdI4blajqtHeOSa7CiKJ+QldQTM33SqtDz5wblA4IRMMXCBYL1/bLZN9Y7zchcp+XXZm3UF0Z09LyZ0r4AG4B8B0A7wDwMwBeBPBO55jfAfCJ9P9fB/DpunJntaatC9Mbh713jjqh7X0D5E65aXfNC5GrMKSX7ykSvtreJjWISwnfTL9s1kO8dASC+cH1wqnyyqly+/S5ePpIvUlnMKkROGvCfxjA54z9DwH4kHPM5wA8nP5/K4DXAARV5c6L8IkMsjV87KPQ+WUMAY7dLPcR0ssbnkQ66XbU7lMepKUKi6X7Mu8JBIJ6uG3FJeombWkUi7qKyE1ydufnfDJRVWcwDX/8WRP++wB8ytj/DQAfd475BoC3GvvfAdDylPUBANcAXLv//vtHu8sJkfe8OkArm/B1hTjj+ChUdIjz9EM4An267SPUbqG7ii61dDSv+cOLNS8QjA43iaG5AJFJwNPKYePT930ePFXavZnTp07DLxtFNMXKEL65LcLCZ8I3o2+rfplkmB/7AjoFwn95o2t5Bh11Iu0WSmLNCwTjwGyOLIuy2tpu+7PgNimrSpc3rXdfMkbz2DK5hqiZy+YqWPgrLen4+PyTp/NUCErpydyjjk6VQEREcawnbENO1+Do/QC9es+Zwi/HAV0CgWB8KI9vhZvhsiyjpa8sl6R9fh9l1nvdMXU+JKuo4d8K4M8BvN2YtP0bzjGPOpO2/3NdufO08N0hovagyVMs8D578pg5eqJ+Qsp521S/2AHwLyaWvUAwOVwP6ZLmVgml8tGA2VGYefXNY8vKnxZZr4SXji4f7wHw7VSq+XD62UcAPJL+fzuAz6RumV8F8I66MudJ+ET2sMu31Fm3a2h1oaLPnI5K3zZljjMdCz/qJ6LdCwQToCzvzShSiDtKcNfG8HnQVZU/LbJeCT/8WWxTIfxRnp6y0xu7mTNd/c714vlyN6RO25Z3vof76AV06NLduvc4OtOnF9ChPzu9JZa+QDAGfB4x5mYuVFRH+ltbugzXuDt9ejypZRpkPQ1UEf7xWwCFsbenV0XgVQSI9KoIJ0/ai9Yaxwb7+9jfD0CK8PaP7SDGSTwBfez168Att+jDo7C4oMlPfwJcfxFA+0k89OJVAITb8FP8PK4DPwKeRB//9dc+g3vxA9BbOghIAcGxXFJYIJgZeKGgMNRNWi9opNHtAh/7WL5QUdXCQUR6nduLF4Hz5+3vfvVXi9drsoCKe62lXLSorCdY9DaRhT9qt2x8xxkyWaMHFJ2z4qpyDT9bwCQ1KY7SBU7yLSHlRufyUEEgEIyNJKnW35tq+FVLnLrHVu0vE7CWkk7DaXOl/MeaefTdPDhPbWgvnihU1vm+yNrOGTuzZtRPlvplEQiWHZb0Gvr3m5bjyjnb2zZVuK6V5rnLivUkfKLaKAZrosU5NuyrQqZLs6j2Ge2Hz0iGijptW/dv3Z0U/PM5odoyvzACwbJjGpOkgwHR5qbdrlstonPn8lFCGGqt3+wAxp2QnRfWk/BrLHxLyekn9OWuPbZL+jpK9vx5oo2NYmbLT56OM0mHPXf2EdFTG7FOgXwmJ/tXNzTJv7rJ++2iE69AIBgJk1jdpqcPL2PIRh2nSGbJiI9j7z0zsGsZ1dn1I/yGGn6SED3bHdCr0N38PkLqnElItfUvfAU96p5TBQ9LXzQuJ1NLsoAsoj++c4u+l5I9dyxqY8PvBrCs5oJAsIJo0hn4FiZyR/Ls+bPtZEHf3l7eFCnrR/hEhTFfPNDRspz3nq1y07LfR99yq3zxTV3vwuZM+oftyPqQl0e0J4OSvJ8x0y9PGk4nEAi8xD5qUjSzXbsL1/kmddnqP306b7rLZOmvJ+ETeeQbVeBZ13feJH9ewNz1BgAU9fvpconGF1x+khRflOzFbBqDLRAIKsFBkiaxu5JLUwc9k8ibEL5p/S+bpb+2hO/+uO6Pl2W8c4KomLx9gR2cQjncTgoZNbVbZuJ9MSyLI6meTF4plwCBYAEw0yIw6XP77naLbb3GG7ug6fuMNt9mLoe4LHZbFeEf28ifvT0dZ0VUfswjj6RBVHSh8N1h+wKuX3dPJpzEDVzARfzWxx9EcHAR//72FgDgu60Otq5fxAt4EDFiAMD2NjAY5AEeOzuASghfetAO2rIq6lacSO+7wWICwZqj19N/Dw6AEyf0XwB46KE8OIrx0Y8WA6HMgKoTJzQfdDr6r3ns6dO6LQNFPrj1Vt22zcCspUZZT7DobRIL3zeM8wZY9JWV6OzyRmjt/4+nQsqDqFLLoV+UgF5pdQgYZl457MPfbufePZ0O0WA3X1WL9X6uqAqj5gtfCgQCb7sOQ/9CI2a8o6npu03K1OLjWEs2w2Eee8NcoKVgvcYFsDoa/rG08Dn8OYp072v2/v2+DrsGgIOPBfjs1Tfjxp2n8Sfo4Mqv7YP2n8TRmT6+j0289d+9iD+4dw+XWjsACF/9KkCK0Gnb13vra9dBuBU/j+t4EhH6N/fR6QR48UXgJz/R17x+HfjI7wY4vH4Sf7gZ4X95135WUQojfO7qSex95ASO/uQufLfVySt+8SK+2+rg6E/uWgHzQSBYLHhAzFb3cKit9uvXgQcfBJIk//7GjeL5J07YqgARcPZBwv/3Qz2y/8TtO7jnLYTffHEHF3ARJ3EDAOGxx6rVhKVBWU+w6G1aGr7b+3PgRG4Z5BY7u1iq9Msr6GZW+2FH6/avbrSLQwVH++fdzU1t4bueAN1zyjbiUzfOsJ+ndXC3Ix4RCASC0lG7qeHzKlNmAGWZps/wrU17551EG/foVevMQj5+q7b43dWvFg2s46StbwbendHPvysuaaj6+ULmvtWsKAxJDRP613fY35nDvu1HVTactIaXbeWEdOfeQ2qY0Cst+4RXWh1SwyUaMwoEC0bVpG2vZ0+mhmHR3dKUd8wyTa/pYhrmcq888dJZIOF7465CZf9AYa7NA0VPnXw1q8hL9pTkXjovoEPfvu98RtSa9BO61NKRt+ap/+jOfDWtfISRxgf4einXJBEPHoGAiPxumebiJWan4G69Xj4CcM/3n1Pkgn3ko+5laoZVhH9sNXwrpekTe9j/ysM46lzAXT9LeGyH8MDBBby8+TDUYA9RSHjgwPaceeBgB1EIRMOPll4kePNJUBjhn/3ONfwfr/483vradfwJOvhr5+7CpdZj+G9fu4gT//4GkM3uE87+Va0FPvwZPS+wD60FfvHpG/ot+tlUwzfw3VYH9LN3AU88IR48AkGKvT3gySfzqa0g0Ptmc2BPHgbP3129Cjz9dO49x03p4KB4jtlO/6dWhAAKTyLCBVzElx7cASlanem1sp5g0dvU/PBdsS8M6UrX3mfPHM6N8+U0TYLq9z16TMeyuOOBypZAdC2AS63IGkVkOr0ToXt5Ux/X3841/EutiAa7ig7bum5HHcODSDx4BIJalHrnhbbG71r+biqVVovoC++Ks/bc6ehkiUcd7aWzbM0P6ybpFFD1y6ex2DxhG8da3lH9fGz3SqtDg8eLy9UrlU/A+jS+U3cXV80CFPW69nEs75w+TfQH9+Uv1jPdmA7bIV3eCPUC6r613ZbtbRMIlgCmPVTmusnHmd8x2Xc6RLu7eW6dTodo93HFTV+fm6il0u4ZQvhExV8WKMzYWEZzqEgN4ix6NoqIFGdb2t0lIt3LX2qxL65f4zMt/BgxXToVFjJzql7PydmjLI+dw3Zo+eyX3gPfZ9W+QLAmYI3fJfxez05zXGfhmyTv+tsvY/MSwq+z8J1D7ZdA5YnPSqK5LrVyjx4meSb/nPQTO2DLkZZUyEFe9rWPXDfNKgt/WispCwQrAJd8fWTs5irkfbMjcBVSlyqWKaiqCdab8D0afuny9MYpXiPaYxKoMKT+tqIYu/QCOhRjNyP4w05EMQYUY5Dl33nltO0CsI+QDtshPdONC/0R4Mm7U6bhS5SuYI3AC5AzGbPaubWVHpC+75kNZHjTsA3ks498Xjqr1nzWm/CJ9C/b6xUjr3q9gvVbm8zS6Q0O2yHFGOSLnWCT9rFN+wjpf3jLgPYR0qvY0Md2Iu1Pb+r324ll2bt+v4VArLIY8UaVFwhWCz6F0pfozNxXg9he7Cix24nre2/+fxyylwvhExV92D0+7V7/fXO/TEdPtx/CXj3BDNh65VTbe44Ko4KUE4bam+fypr7WkS/vjmGWWLdROjwRCFYLVQqlz3+BvWfGGumm7cgcEZg6/yqpokL4I6D0JRuUa/i1WxiSulkM9ePALSspU0e/zEpRNmk82M0rc6UX0TPd2C/Ti4UvOCaoNb6Uf/ESpTRZF+a+qtpBrJcrjUKVna/CiJ7txd4Ea8sOIfwRUeroYobmKVVc8r5s86TvU6HW9PehXyzfcJJIZ9g0L8nePOxpYHkVcXqIVR6PCpYSi3D+qrJfkkSvNW1+d+qUqcE3HOkqZRleZjvK9t0cVkveloTwpwEzq5LJzml2NOVbLYU39vNylsO6vJkupFIyrcB+vtmLPkycQUVi87kbK+4m+VjyF1WwnFik85dPoUwS23UyDDXZV859mZq+2wxUMZdWRvbOfMAqaDxC+JPCZy3zVH6/T1d6moWVswLyq5tGJ9BuZ52F6of08kaPYsSFxE+Z01BK3skwDe7IvIAGOl9IOlEcQweAZXVkh2Gzg6LqiSuBoAxNpJV5XNvl7a0t3aSKqqoRDxNF1hoUVnuKPRczCgKU3xV7BUbNQvjTgO/tC3VA1LNd7Xb57NnHSfV6dPQ3H6V9hPTsuV1S7ZT0u12rjMN2WFggnWUa4gCvdMI2+ekwmxB+AR3qnc0XW3kBHYr6iS3nOC6avIA77y9rhKBgObGIqaGmGr47AgB0gONhOvfF0e2HnYjUIPanMq6y8H2OGktM9kRC+PVoKlD6xpdxTCrU5K5XwHmcXkBHk30YET3+eCF0T0fvKm/It/mCuWmSXS8gV0ZyvXfcF5W9fVbASBEsGRbh/FUnJfn8JjodonPnNMlzs9Nav+36nAVTNdXw533zE0AIvwpNBcoSC5/fOLWxaZGyMpNwOCaISlTBuk/VIW25DD1eBgABtqeP69Pv7bicYeqKGCmCJcK8LHzf68ub7zOuU6+rExhyU2udUgbRF7dC5KzjpcOdwLO92PbOW5HGMzPCB3A3gH8F4M/Sv28uOS4BcD3dnm5S9lwIv6lAaX5uRngAxcQbvs1581Q/TAO24sKh3S4XnRTKKSzEUpVmwXD75M10/1zi91WwRJiXhl9md21tldtjcUzawy3UI9ckKa41q7dixlqT9DNXTPNelGj4PiL/fQAfTP//IIDfKznuP4xa9tws/Kbmi2/tM5ZK2oWZI//mpHXYR54/x+buhP70dpvM/x/cTgTQG7enIwdzBGFo/hRFRMNhNmGswojUMMm8FjTpDwsapkBQhto0BhOiqlNxNXc3ENKcuzL32biJYS8+vnGP7hCe2ogthztP0L3GCuanmiXhfwvAven/9wL4Vslxy0v4RM0FSvNHN45nL52yTZ1pF8TGK9AZMos6ZO5l8AI6BAwL5H/UCUnt7urRxblzuccQp202VuJ65tyADtOMn/uIspW5YgwyV09x3RSUYV4Wfpnd5bGvisaKc0BO8HlbOuxE1N+2kxqGfVWXViu/RtX+kmGWhH/D+D8w953jhgCuAfgKgF+uKO8D6XHX7r///lk/F41RBUrf8UDujdOyNfw37tRyjtq2g7Q+ed+A53zJXWclt0q0n31y05Z3rnRD2xnZkwNWhRFF20P7BX90mE30vtLq0OG7dm0TyvHbX7H3fDQc65ubLsZpIlX7ZZ+7kbMl9lVR13cmVTttO9W4u/g4LzjkDr6PyyswEeED+DyAb3i297oED+D1kjJOp3/fAeD/BvAf1V13qTT8quPZuu529Vuzu5tb2qm/PQH04zs2rZfulVaHkptJwcLIg3fzNW+POhUjiE4njyl3WoX+yLMurzcJSdGVc8VGss2xgsP0RaPpILjpo3WPK1vbJ0mK2SvZvilLJ8IWfr+v29Ng1668SooOE8eF7IkmJPyqramk45zzzwG8r+64pfLSMd8GjvgwBc12W39uvr3peWqYFNwrs/00EKTb9fP47uO5t461xKG5bW/rt9ppFSqM0oXYNelb57imlNPKfK6bKzBX1Qzz0iiOEZpa+OP4QDCpu9kvXQ3fN5AN+0UN3xzRqsTvR3+YukWXWfirPvibJeH/I2fS9vc9x7wZwG3p/63Uo+eddWUvjR++2SGYgrubzMYdExr/q/PnrReOpZWX37ZVuiLPcKiLfu68TqAWD5T+sIyojdbBss0+IjsoizuDqjQQab3n5Yq3EBzrm5suJhkEN+0c3IGm+f3WVvWqVdw+cyMl99Jx80upRFnOCz4NfzAokv+qDf5mSfinAHwhJfHPA7g7/fwsgE+l//8igK8DeDH9+/4mZS9FpK3v7fZlyDRz7ZvnxbHXwmf3ysO2b5UrojNnjERqodIpX31j3hILnZJEZ9U8l0/gcrK2LN0DS0EVrp1Nh/EriWN9c9PFqApY00fry3bplsN/WZ60yuT31B2RGvvstsnHxgOddpwzzrKDW7ebk73Pnlsle0ACryaB6fLoWNKWaeAxgVQYZi6brob/Ajq0j77XF5+3Z7rpy2qOcVstPU9QY6HHAx3cxcmfONjr8PwgX5fXHDdz4jXDGjq2RrBY+CNjlInYJo+2rFlxKii3zGd7tntlZr2nvU5lp+Rc3M11z024zJ5btVdDCN/EKAJdiUdOYfONOdO35HAr1vq746VzBVq4dxc6zzcn3bHZQqos/TAkNUzSBhJmoeFRqHR+n54edRBRPhwecgJ+fb9qEGf6/7GTuUXDnxmaPtq6FatMWUcpyt5dHqm6++aI1K3PKHV27bdVHfwJ4TPG8c4ok1I4DwK/LW5ufCOUTw0Tr+lgBl55ubtfTOqU1SUMyZrtNV18NjZInT5tXSNrIL2eld87Hihr2Mudw2dO53n6mz6qlcG6e+nMcFay6aP1WfitVm7hm+dlRpPjtPDKfd3SebNR4EpQJbbbykAIn2g8y67KwjfXx+12iTZtySYbnzoSj3lMUcPPc4BoF3tNvla5LL24/mrsFsr7Z8546622+xSFKlvJJ+onmWdDbj3l++ZzcXOQrNrCEBZmSHpLjSpGntIzaVqMz1mMm1UmtZhGjy9TWtn6ziPU1dfE3XqsEukL4TPG0W7rxEZ3PGpOhPJ+HHtln6/f2aPtR3VWPzMEnLP7XToV0r89XXS3zCSZxD9yYFnn8Iz93cv3dLNOJor0yEOFUZbe2b2O+VwKPDEY5N5DXJd1spJXEVVGT4XjwayrwptrM2X+9oMScX1CZvadxjaUac+t2msthG9iFO+MpqOCsk6BJ0bN6f9I+97zZO6lUyEBZkRsSEBCT92Wv+D7CLXHjZEjP95N8kAv14pPFPW6uYzDW+Yd1IlI7eZ11gupD/3PRSn7tkNFaneQeR4dddIOyJtoXDAtTG1AUiZalzoeRNbFygbCo9TNbUaDgT8HISuWUaiKMSi+5UXHePd8A54sqKvh/SwbhPAZ41j4Td+Imzftcm/eLL7Vhj+/ShS9vNGjV3A6I/lLp+yX+gp62XccsfvG7S2KMaBDTth2xx1Fwt/epssbfovIzbipOh2KtofFz8O8zuzVEPZ1J3IFui6uu6mdaFwwLUx9ysFn9HjahpUyuOS649RNqXzQaw5S7WUKc15PhsofbuuS/owlqFWBED7RZN4Z5nf8pppvuE/D93Uo/KanGrwyXlhN7vYLrHPpcBIoj0VjbJdO9enbd+YErO65J/t/H2F2PnsH+TbV6diavjGiGOwqI2glLIweCMilJs9jE4yHSV7b2gLdd9TpCOo8tbx1c85xK8ivPw96k4Sov62cDOL5vFavm793vMJctlZEVXrwNYYQPmNSU8n3hptWBmv25ktoRpM4Q+Won9ALaJcSMK9QpS0gP+mrjQ26dKpPgKLOmYSyJG4G2es0sQNS7Ta9fE854VOSZI32Sk/LPu5CLLkbaXHGjes7zqM99pjAjBxnYFpbkMviZQn4QlV5XbMInocyFxIxXwLzWObqpzb0ORv35GlAzHz2Gxu6XHu1tsRKT7Kys6szghC+iUnHb77Wl69aUtw8LaTgaulsQJJl+OMIWdO10j52qMm+k6opjutDPkLIz728EeqMmyV1NS0zNxMhk70rAVma/ryWUVyVsXhDQ6PqdqqmnkZ6DGUSJUsmHg2/EOHqQNfN493leQnsdMd27htfOuO8c7BXa3vuXQPbwEgt/8OtuFjBNYMQ/rTha30ey94KnFK5Flo4v2Apa93+sBPRt+/bon2EuWbvbDpnfqLJ3tMZ7UM3nHA7ybR3a2u3C1G2BWJyRgzcebzS6ujJ4/QYzrNv9R+u6yaXW7XfBKviS99Qk6nzliyz8Md6DL7n7ylIhToFgXldX8ooH4HzZqU1iDM10zikeA4nN+NBx2Bg33veaSjnkdoZXtcVQvjThK/1eYbDV3pRnrgsjvOGOPCcD71qlqnpH3VCih8f0pWePlY5yyR+7UTbIOG+lnMMf+Wt83kecJ78jeF49HS7WpIyWYWZwrhPXnSdGyaXF24n1tBdveu8vYxi4mGfaRD11IXtCVHXgdVoMr7q8+tU4UBT+d3Yko/xb6+rsuuYk6u9HmWrRbl+CZ22bcxYMR/pscXpLvscHiHGsX8BFK7PVGSuYwgh/GmhqmUab6K16n0/sSyQTBPnFpISueobrarX0zlw0pkxVwK61GKvnj6pe++jF+/oEaDoqY2Y1HafVKhTKAAJHbZzrxpeYD3beJlE03vItN4GcWZt8bD6Ukvrq7zIRHbqMMnujUk/W/qRD3Jy7k/EUBUkOlelp2kHVqXJkP922Jpmr19fFOqo+r732RgdD/+jVO4WyTF9bI/wfhw7K3+GRWvd1PR9QeutU/5RAXuGua+H6SfvBm4J2WsI4U8Tvgbe6xn5WikjaU26TkMcOOczyZuWtTvEdshi67wmYPX4btaCLrUi2vqlPOhLdXtpkFVCPs092wxXyjimghte71xCvZ4R4DtUWV9VJB1VmOQ9amvPiuwgUz6qIepaeEh04gHEKL1F05FGQ1Z2b6cqpsiVVZoQn+/ZPNuL6UrPWPAmHZWpQUxhqFfQNMvu94ujCu403PVkfRHbtjyTH3N5Q6cr5hEtF+42A3OU4cZiHadVqyaBEP604SMF4zMmTtcCyTIBlpGKjyB8idnC0CJ7X9AXy0Pc+ErdMdNKceM3LTJ20eRUsqYHjnsbWcN0Jnmv9MKiZe8xzcYiZodETa+SsQYQ4zqVV5F5w07BV0wTq72phe+tRmhPkpod9pGho1fWx3kJVN9MF2JntPTVlT172KL3PXO3c/MRvdsJrTOE8OcI86UuW0nQC5dsfOPfXi8LScwWMWk5Mg1A6uaQVKLo8kZUziCuWeRYaDEG2ajAStFcQX4qUcUAGZclEs+CLB6ibvyQHRJt4kpYKKumzNpK1ZnYNR2J71JNMja6doFvv+wWs2cTFiVDNhJ8tgZvSWLfFxsER53QSvnNxgMv6MN1s2OmVJaL0KrwCI97ZIPhGEMIf87wcXVl5oEqsjGJMQ0vV50OPXtuN1/MxNkOO3p47E6GvQodjKU4jp3/9np6SO9p/K9udmySdqxS6x7KF+jVm9Hq3URt2ciiaUK2ChJtKnE06mSbkH1TE7tiv0y6qKtGHBtpcBJlLegRD4ojJ++zcT5kC903cZy9z+38XVGhft8ub+r9T94ZUaetsp+/39ev2tvepv/nuQGWh3jfVEXHedzrbtkzhPAXgDJXslILxPcmO+RTuq6tsbE/vC8aNluEhVtwao2rXi+zil0vC3O7vBFlhOzej1JUmiTOaplGignNM3mgjRrEdu6WOrPNQ6JN+dfXyXrjI5qS/aST0GRbqvwz1RXJnUOM2Ojo9XM97ORuldnEqnOL/e1yC7/XszV77sPZYcwkfd6ea0fZhD5nfeXsHxsbOolrmrk7cxLjQWEd4U/5cR9bCOHPGWWkU5tmxjW/0hPNcHJf8JXbQbyATqbZZwTKrda3WpfSyyh22kWPCdfy63aJBrvKGrHwRG4cp/fgavQeDV9Z6o/KCIqgXVrNxaezezBHFp79kQnB80O56/1e6flONDCDeIBRi7Tfj4jMyVAeOZk/Qd4n58ddaunzPnG73mcNv9vNvXL4p+33TddK+51Vib5Wngwtl9e0+qh8KqQt6dBYAztBCiH8OWJsK6TMwk97CZXYict4e+P2FgE3swb/tVt08BNPhnXaepifle8OPVQqo4R5KmZNGkXPHtb2mVQ6HaLhTU02MWLLy8Pb26XeH6ZPN6A9Qdjt0+1wsgReqXfT8KZKPUs0az13PraWu2NJY3fXfqy+pfOyAxz2qVsXwFtG1f4Y54xapH4/osLvBeTeT6YXFqtX5u/JhHzUiTIvHTddMI8EtVzjif5OX/IksVN+M/mbaRM8dkd2ryMO7AQGhPDnjJGtkIpegn2SzWXd9qFz2JvJzYAkI+RWK03Lnw6tNzeNtAtugFg3pKOO/uzlzTw7J5PHC+hkpMD7PksyS5PseuMYZuXg8STTmwe7WuN112npnrMJmOciuEP7xG22VGXumxO25gR5qZzm6WQ5chkwFt8wJlfr4COjwmeua+6EZmp+G/azy9Z/5UvFtmSWDFX2zlh8bXhjmY+ICTzsqyxzKkEH4pnzPyqMqHMmsUYZvjQKZaQvMs1kEMJfAEa2Qjy9xJWe1mCZo5/pap2WIyA5klcN4sLw2SLRrm7E1spBDvHrYC6VdRK/96Y4s/TzfDxs4dvn7kPn+CcqIbNej/7t6V5WNqDosB3RUxuxU+eixZiTeLXcZJKIO0HuJRDjC87J4vqR80SoeXgVJ/vmgM0QDaVs+WVauj+TfZZVMt1MTd+akC/IhPmz88mO5lwLlwuoLGX3ES+qEypKQv5didpninVyyd43ty9kPxmE8FcFxhtuDp/ZKmbdNVsJKCYrUZlv297OyYaDbEgpneLBYw2ahNk6pZxheTERG59rhrtnq2ClN8ITzSwxuSka3HKvoFfQ8PmYshs169807F4NYsvfXFuttuRQsDpDOzDNnEcwedtdE4cnP7VEUgxQYxIueR0qEQ/y8q70NMkz0XKHmc2vOA+GtXvztzdJ3z7F0+mmKYvNOZytLa3h8/vq/k48ics6/xRT2wtICH9lUTb5m1mLDsmUuYNac6ahouRmUpH2mOjxx3NNnfPx8PefuM0me5M0Sl1Plc8TJJdiDts6wRvPT+wj1HJKYkYsN7fw+Vk0IZBk6AYXqVKr86mNXBLhDvlKL8rMfl9H0+kUy8pGDsaHWaTrmCqPObLikYTppmq9VA4BR5G9Mic/P6/SGCrvg3Wfb5Jo2c59zy619NoKnBSNO8dVXzh8mSCEv8Jw5xTdRmDmMinzhrQakWF6svbv6vJRqGh3Nyf9gjWKPM++2yGUNlQP0XgzJbZDK2lXnYb/lLGfy12qkTt99uyGthQ22FWWY5Nr3fomdZmw3aSpvs03sX1krjVsjO5GIr46HdFjQfDiNlGkJ+CZgLmzsWQqT8ddNjIxRx0Upct6Gvuue6/VqUTlv5mgHkL4K4oqC989LpXKCchd3ArDZJNoOh1KbuYpk19Ah9TuICOyZ3txtiSva43mMkxOhCyDlJF94vH1vv22okyTDO1goYyUd+NsJPF7b4rp8oYm+kNs0VO3hXTY1hGe/e08wVsVgWTXQOzIGvr8ZHdQIPz2mfJEX3wNn8uhvXk0/J6thfMzfbYXT4/wKjQn7vjZO8d9Xpl05brJOnMgfBmWrooLl+QeQCY4fMP8bdxVRAXNIYS/gvAOp519F2bDKVqoBnl7RP999K2JPdb6fccfpgnR8vLtxSnK5Jx9RHTHm0xvjWJwGJOJqyMrRVnSuOFQ30e4nVhlmfVvsvpWMiyOUHj/T2/X3jr8zC5vcqdWlGKIdJ1csr/9dns/m3tJO5nBrv1DqV6vMJKYqolrDgeNjt9c2axypObMSvPEr69zNeecCu+fp1gh/OlBCH9FMW4uL7fhFbT+oU1abvSuFehknBgP7MlA9gBhL5TSDimOM396LXnYMQVZ8qyUrD9zOipdPNvtCLza/ggToIPd4vnu6l15h1XsoA6NjmV3t9CP0pveRNZ6rWHIk5U5qaphQsrpna2lAj2/cWG/Qs4pdMCOxW5eM1vLoKw4R7P3OQ2YIwNrhFlyLyLpTBdC+CuMOlm2DKWdxcBj4Vdl6vJYdcy+PsustENSypnUtD2ATC8ZdutrIslwWVaH1ZDs8zIcWelmcfb10qmc7HO/8jDr/FSijLQG+Xbzpn9+ha3dw62YDtt6GUvzgE/eN/CSnu93Nb2v3B/CyGRBRE5Mgndepfis6wwPH7E3lSOLv6WQ/aQQwl9TFDqLxGZPjrAtsJTb2ppaj55983NDRSAgt3StvqfvdgzVZO9L+HbYQM4xyb6Qw5+DyBwy5AW3ubO61Arpy92Q1CCmXtfurHjb2Mg9XixSHMSkwjBbuvLrt7StA15pdSjcTqwIYa81HJavJfvlXmQtVGN6cfmem2894rpJVd9vNc5KXE1GA4JmEMIX5DDNNVPor9VlJsfWVp72gD1aXG8ak9SrCCBLE8DRnf3QytjIC69U3Y7rSWIFRTkVyyd2be8iFYakhkmWinof2lo3+1EzARmPSHiu4LAd0gto22R/qp0ReBSqUtLPOkMPebtRx+bxVtKzVKpzn4Nrwfty0Fd1BlkGzwajP7Hwp4uZET6AXwPwTQAKwNmK494N4FsAXgLwwSZlC+HPEK5GMEoikwku6VqGZYOL4bCZX7ZSZOUHToZ6kvkoXdaRPYeycz0FuJHBaphYDukqyectfHlhzEp+8s6I8hQXcebPz95TtpVdHldw6VSfVLeXZbo0lzfk+y50hp4PfUQKFHPZE9lSXfZsU/j85N0lDn2vT5PRn2j408csCf8/BvDXAByVET6AWwB8B8A7APwMgBcBvLOubCH8OWLciYIxLuMjoI0N2wpmb5ds0BE6Pumu3JSemE+y2lq7UpQRe+b/PbSZybplR/hOhor+yR3FpF/tv2nLPrbXUJSt4uRLuaBPSYoPI9320SdAZVkpeQLXlKKqLHy+oBuAVtcBuhgMzMyXRdLnOtUUUwnJgDldzFzSqSH8hwF8ztj/EIAP1ZUphH884RqiZvI0XijD3Occ+YX0AAYbWFJMRphRRornzuZW+mEnovPvyiUVn4SRVZT8BFtm4bvXNpUzk8zy1MSeIU5WRh7r8Ew3dvz9c3/2Kg1fhZGRv4iy0YLZgbj36/5WVcF801xDdk42x1pg0YT/PgCfMvZ/A8DHS479AIBrAK7df//9M34sgnnDZ+G32zrfj0sm/b49yeySmTved4PDTOuYCTojeZZOWlGWy98t0pU0bB3c9d+3zWg1TLxatkn2PNfw9Tt7tA87Qm4ffWM0Ua+3+7x08pQUxfQaT23E9OUSrx7fb+ZbtVL85JcXExE+gM8D+IZne69xzFQI39zEwj9e8Gn4ZWuwA4a/fZMZPc8x+9D5eUyS9qcPLqYwsOa1jQCzXi/fz11Kk4Klzwu4mPnkuTzeZ6s9uZnQpVOhU3czktm/spSlv5fIM2acVfaYQqXJ3iyjRjT3yULTtO4F08WiLXyRdARE1GwJWS+hVLnsGITFmS/NNL79bfb3D7NVwDJiRpf2oVMyuMYuoP3bVZino+b8+DyZunmPrdmbOeBZgvH1V9m9JflawYcdu/PYR0TbjypLgil9BhXwSiUjuMWUJeSrOEWwYCya8G8F8OcA3m5M2v6NujKF8I8nXK525Rwz/0+WVbKGnMw0xxzQxCH/Oh2zbUEXV/JKs3OaOnvfXrAjWwwFeQriJNGdgpWpkydQzTmGEq6O47xTMdei3UfeqfT7ui5XetXPYOTfwZHAfCkPzH7BlYVczyHB8mCWXjq/AuC7AH4C4IdsyQO4D8AzxnHvAfDt1Fvnw03KFsI//vDpw6YroxUVXCM/uOmFk6GZMVMVlobMiK7Xs7Jzulq7L22D6dWjyTyXX1x5pcyYthYD6+u1Y035ib2W4oGzcI1Z4Ji6ijUPkG5HHT0/4BbHDkuWLBQJ2S8zJPBKsHRwectZDTHX8Bv67LnugUlCGYl2u5pUvYTvavQG2TP5llnCdcqIr3+qi0Uwo1QzPjdiDaweqdcbeebU9GhiN1bev7wZFUY6PnIXGWe5IYQvWEo09r92191z9svK2drKZaHCpKfB0Hbmz3xCtns28S4byPnyXUPbN/iYJA9N9hxM/0guzN0fAaYExp3aYSeXkSQAarUhhC9YWtT6X9cwZq0VPTQ+cN2CnIhaU0Mn5Fkz+e+/vkP/ZXdODlLOlpss6bBcV09z0rpshFB4DiNMtDZBfr9pR1MywSxkv3oQwhesJurY3HBLrCQq7jRMQdrYNxf9MC1es0B2wXQXV3EN7TJyrPJQamxRV3krjflYfR3NFC4hWCCE8AWri4ZmZy1RlZjZvgRrHM1qFsjyhxmo1dQKruq3zEXDSyWtEZ5DHarqYs4fiIW/uhDCF6w2ath8Ui4sqEYVqRp4gDCqFVzlreMeV3nyFAR2n0rmCxITDX81IYQvWF3UsPm0uDA7ziiAJzZNr5amC6SXXWNsuWTKGcZ8cyeSxOx4QAhfsJpoyOaTElWB/NLMmplvfxr8pQbx6Lq751bqOorSiezaGe7JMYdLCGaMKsK/FQLBsiIIgJMngSgC9vf1/v6+/u7kSb0PYG9PU2i6mx3G++Z37v7eHnDjRn48EbDz4z2cvIuw90SQHqsLDIIAjwDY2qqsTgFEwM4OcPFifiu8D9h19dZnR5e/t+dcoOyCE8AtcgaXECwQQviC5UYdmyP/2LdfRaBxrL8ziTcn5sDuKNJ/nniiUXUKdWnQb4Goqj7FjksgGBllpv+iN5F0BJOiiSI0T9/zJnKJ+MILJgUqJJ1Af798OHv2LF27dm3R1RCsOEw5hWFa2nzMiRP590ot1pJetvoIVgtBEDxPRGd9353wfSgQHBeY8gnDJfudHfv7nR39+dxB5NSHFlsfwbGDEL7gWKOK0N3JVKX034sXF0Cye3ugCzvYuUC6PiFBhTt4tre3mPoIjiVk0lZwbNHEO6bJZOpcKnrjBoKDi/j1HoBwH/vYQXBwEf9FGCHqEU6eDETWEUwM0fAFxxrVbo76mCq3zbmhYrKBIGQvaI4qDV8IX3DssRSE3gQyWyuYAmTSVrDWWIlgoqWaPRYcVwjhCwSLxlLNHguOM2TSViBYNJqG4goEE0I0fIFgWbAykw2CZYZo+ALBKmAlJhsEqwwhfIFAIFgTCOELBALBmkAIXyAQCNYEQvgCgUCwJlhaL50gCP4SwL9ZdD0aoAXgtUVXYgFYx/tex3sG1vO+V/me30ZEb/F9sbSEvyoIguBamQvUccY63vc63jOwnvd9XO9ZJB2BQCBYEwjhCwQCwZpACH9y/MGiK7AgrON9r+M9A+t538fynkXDFwgEgjWBWPgCgUCwJhDCFwgEgjWBEP6ICILg14Ig+GYQBCoIglK3rSAI3h0EwbeCIHgpCIIPzrOOs0AQBHcHQfCvgiD4s/Tvm0uOS4IguJ5uT8+7ntNA3W8XBMFtQRB8Ov3+ahAEDyygmlNFg3v+rSAI/tL4bf/eIuo5bQRB8E+DIPiLIAi+UfJ9EATBQfpcvhYEwS/Mu47ThBD+6PgGgP8KwJfKDgiC4BYA/xjA3wHwTgB/NwiCd86nejPDBwF8gYh+DsAX0n0f/l8i6qTbI/Or3nTQ8Ld7P4DXieivAtgH8HvzreV0McL7+mnjt/3UXCs5O/xzAO+u+P7vAPi5dPsAgKfmUKeZQQh/RBDRnxLRt2oO6wJ4iYj+nIh+CuAPAbx39rWbKd4L4F+k//8LAL+8uKrMFE1+O/NZ/BGAvxUEK53L+Di+r41ARF8C8KOKQ94L4F+SxlcAnAyC4N751G76EMKfDU4DeMXY/2762Spjg4i+n/7/AwAbJcfdHgTBtSAIvhIEwS/Pp2pTRZPfLjuGiIYAfgzg1FxqNxs0fV9/NZU1/igIgr8yn6otHMeqLcsShx4EQfB5AJuerz5MRP/rvOszL1Tdt7lDRBQEQZk/79uI6HtBELwDwHNBEHydiL4z7boK5o7/DcBlIvpJEAR/H3qE858vuE6CESGE7wER/e0Ji/geANMCemv62VKj6r6DIPhhEAT3EtH30yHtX5SU8b30758HQXAE4OcBrBLhN/nt+JjvBkFwK4C7APy7+VRvJqi9ZyIy7+9TAH5/DvVaBqxkWy6DSDqzwR8D+LkgCN4eBMHPAPh1ACvpsWLgaQC/mf7/mwAKI50gCN4cBMFt6f8tAP8JgP9rbjWcDpr8duazeB+A52i1Ixhr79nRrR8B8KdzrN8i8TSA/y711nkIwI8NaXP1QESyjbAB+BVoHe8nAH4I4HPp5/cBeMY47j0Avg1t3X540fWewn2fgvbO+TMAnwdwd/r5WQCfSv//RQBfB/Bi+vf9i673mPda+O0AfATAI+n/twP4DICXAHwVwDsWXec53PM/BPDN9Lc9BPDXF13nKd33ZQDfB3AzbdfvB/DbAH47/T6A9mD6TvpOn110nSfZJLWCQCAQrAlE0hEIBII1gRC+QCAQrAmE8AUCgWBNIIQvEAgEawIhfIFAIFgTCOELBALBmkAIXyAQCNYE/z/7OFoR/mkeUQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "b = 50\n",
    "s = 5\n",
    "plot_samples(target_samples, x[b::s,-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualization of the Gibbs sampling Markov Chain Monter Carlo sampling process\n",
    "\n",
    "We can visualization the sampling process in two different ways:\n",
    "\n",
    "- first, we can plot at the sample of each markov chain at each time step\n",
    "- second, we focus on a single markov chain and observe how it evolves\n",
    "\n",
    "Note that we use the values x of before for the animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.lines as mlines\n",
    "from matplotlib import animation, rc\n",
    "from IPython.display import HTML\n",
    "\n",
    "\n",
    "# this is just code for building the animation\n",
    "# http://louistiao.me/posts/notebooks/embedding-matplotlib-animations-in-jupyter-as-interactive-javascript-widgets/\n",
    "def animate_mc(animation_data):\n",
    "    # First set up the figure, the axis, and the plot element we want to animate\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.set_xlim(( -1.5, 1.5))\n",
    "    ax.set_ylim((-1.5, 1.5))\n",
    "\n",
    "    sct = ax.scatter([], [], lw=2)\n",
    "\n",
    "    # initialization function: plot the background of each frame\n",
    "    def init():\n",
    "        sct.set_offsets([])\n",
    "        return (sct,)\n",
    "\n",
    "    # animation function.  This is called sequentially\n",
    "    def animate(i):\n",
    "        #fig = plt.figure()\n",
    "        #ax = plt.axes(xlim=(-1, 1), ylim=(-1, 1))\n",
    "\n",
    "        data = animation_data[i]\n",
    "        sct.set_offsets(np.empty_like(data).normal_())\n",
    "        #sct.set_array(dev_labels)\n",
    "        \n",
    "        return (sct,)\n",
    "    \n",
    "    # call the animator.  blit=True means only re-draw the parts that have changed.\n",
    "    anim = animation.FuncAnimation(fig, animate, init_func=init,\n",
    "                                   frames=len(animation_data), interval=20, blit=True)\n",
    "\n",
    "    return anim.to_jshtml()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# show the samples of each markov chain at each time step\n",
    "# you can clic on - to slow down the animation\n",
    "HTML(animate_mc(x[:20].numpy()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is just code for building the animation\n",
    "# http://louistiao.me/posts/notebooks/embedding-matplotlib-animations-in-jupyter-as-interactive-javascript-widgets/\n",
    "def animate_mc_persitent(animation_data):\n",
    "    # First set up the figure, the axis, and the plot element we want to animate\n",
    "    fig, ax = plt.subplots()\n",
    "\n",
    "    ax.set_xlim(( -1.5, 1.5))\n",
    "    ax.set_ylim((-1.5, 1.5))\n",
    "\n",
    "    sct = ax.scatter([], [], lw=2)\n",
    "\n",
    "    # initialization function: plot the background of each frame\n",
    "    def init():\n",
    "        sct.set_offsets([])\n",
    "        return (sct,)\n",
    "\n",
    "    # animation function.  This is called sequentially\n",
    "    def animate(i):\n",
    "        #fig = plt.figure()\n",
    "        #ax = plt.axes(xlim=(-1, 1), ylim=(-1, 1))\n",
    "\n",
    "        data = animation_data[:i]\n",
    "        \n",
    "        sct.set_offsets(data)\n",
    "        #sct.set_array(dev_labels)\n",
    "        \n",
    "        return sct,\n",
    "    \n",
    "    # call the animator.  blit=True means only re-draw the parts that have changed.\n",
    "    anim = animation.FuncAnimation(fig, animate, init_func=init,\n",
    "                                   frames=len(animation_data), interval=20, blit=True)\n",
    "\n",
    "    return anim.to_jshtml()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# animation of the evolution of a single chain\n",
    "# we show the 300 first steps of the first chain\n",
    "HTML(animate_mc_persitent(x[:300,0].numpy()))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.4 64-bit",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.4"
  },
  "vscode": {
   "interpreter": {
    "hash": "33184612fe23ea789693bd8266f7c5cf5860109e6e433b1efdda378c5f430ce8"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}