01_CloudBlog.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Statistical models and PPL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Bayesian Statistics in Climate and Weather\n",
    "\n",
    "In Bayesian statistics, probability is viewed as a degree of belief of an event occurring (in contrast to frequentist statistics where probability is the frequency of an event occurring). The Bayesian view of probability is valuable in climate and weather forecasting. It allows us to quantify how likely we expect events to be, such as the probability of rain on a certain day.\n",
    "\n",
    "In Bayesian statistics, we can incorporate a prior probability: our degree of belief of event occurring before observing any evidence. On an average day, we might expect rain to occur with some probability. We write the prior probability of rain as $P({\\sf event})$. Then we observe some evidence and update our degree of belief about the event occurring. Let's say we look out the window and observe some data: we see there are many dark clouds in the sky.  So we can update our degree of belief of rain occurring this afternoon. This is called conditioning on observations. \n",
    "\n",
    "We write this updated probability as $P({\\sf event}|{\\sf obs})$. This is a conditional probability, or a probability of one event (rain) assuming that another event has occurred (cloud). In particular, we call this the posterior probability of rain, because we have observed some data (clouds) and updated our prior probability. We can use Bayes' theorem to calculate it:\n",
    "\n",
    "$$P({\\sf event}|{\\sf obs}) \\propto P({\\sf obs}|{\\sf event}) \\times P({\\sf event}) $$\n",
    "\n",
    "There are two additional terms here. $P({\\sf obs})$. is the evidence or the overall probability of seeing the cloud, regardless of whether it rains or not. $P({\\sf obs}|{\\sf event})$. is the likelihood, the probability of observing the data if you know the outcome, i.e. the probability of seeing the cloud if it ends up raining later. \n",
    "\n",
    "We are usually concerned with learning the posterior probability $P({\\sf event}|{\\sf obs})$, which is why Bayes' theorem is useful. But it can also be used more generally when we want to learn anything given some data, e.g. we might want to learn parameters of a model. This is called Bayesian inference.\n",
    "\n",
    "This happens to be really useful and has many applications in climate and weather science. It becomes even more powerful when we think about probability distributions, which describe the probabilities associated with all the possible outcomes. One important use of Bayesian statistics is in data assimilation: incorporating observed data from satellites about the current state of the atmosphere into the weather model for prediction of the future state of the atmosphere."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is a probabilistic programming language (PPL)?\n",
    "\n",
    "The traditional approach to constructing statistical models can be time-consuming and requires expertise. Typically it involves statisticians writing down the problem on paper and carrying out a mathematical procedure to design a bespoke statistical model. This would then be translated into code to construct a model ready for use. Since the model is specifically designed for the task at hand, this method must be repeated for each new problem.\n",
    "\n",
    "Probabilistic programming languages are designed to make this process simpler, quicker and less technical. The idea is to remove the need for the expert statistician to hand-design the statistical model, by embedding this step within the programming language itself. Statistical tasks such as encoding observations (conditioning) and learning latent variables (inference) are automated and some of the complicated mathematical and statistical steps are hidden from the user. This makes statistical modelling more accessible to non-experts and tasks that traditionally require a lot of thought and hand-engineering can be automated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What does Probabalistic Programming involve? \n",
    "\n",
    "PPLs are just ordinary programming languages, equipped with tools that make statistical modelling easier for the user, reducing the need for hand-designing programs. The main components of a PPL are:\n",
    "\n",
    "1) Sampling: Our model is probabilistic- it requires drawing values at random from probability distributions.\n",
    "\n",
    "2) Conditioning: We have some observed data that can be used to update probabilities\n",
    "\n",
    "3) Inference: We can learn something from the known data and model. This could be a \"latent variable\", some underlying factor that you can't directly observe but might influence the data. This is usually the hard part.\n",
    "\n",
    "Probabilistic programming is a suitable choice when we have a probabalistic model, that relies on sampling from distributions in order to make predictions. PPL makes it easier to do conditioning on observed data and to learn something about the model. Note that probabilistic programming languages are not new. TThey have been around since Simula in 1966. However, with the rise of deep learning, they are growing in popularity, with many new PPLs designed with inference for machine learning in mind (e.g. Pyro, Edward, Infer.Net, webppl).\n",
    "\n",
    "Let's look at a simple example to show how probabilistic programming works. We'll be using Pyro, a probabilistic programming language built on top of pytorch in python. Most of this example comes from the Introduction to Inference in Pyro tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A simple model for a moving cloud\n",
    "\n",
    "Let's say we have a cloud moving through the atmosphere with some speed in 1 dimension. We observe it after 10 seconds but with some measurement error. We want to use a probabilistic approach to find the distribution of the wind speed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1) Statistical Model\n",
    "The wind speed is something we can't directly observe- its a hidden variable. But we have some prior knowledge for what kind of wind speed to expect. We use this to sample from a prior distribution for the particle speed with a guess wind speed and variance. Then, the particle is moved a distance of speed x time and its final position is observed with some measurement error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "import torch\n",
    "import pyro\n",
    "import pyro.distributions as dist\n",
    "np.random.seed(123)\n",
    "def propagate_cloud(guess, guess_sd, time):\n",
    "    \"\"\" Returns speed and observed position of cloud.\n",
    "    guess = our initial guess for the speed of particle\n",
    "    guess_sd = our guess for the standard deviation of the particle speed\n",
    "    time = length of time before observing \n",
    "    \"\"\"\n",
    "    # sample from prior distribution\n",
    "    speed = pyro.sample(\"speed\", dist.Normal(guess, guess_sd))  \n",
    "    position_final = speed * time\n",
    "    # final position is observed with a bit of noise from our measurement tools\n",
    "    position_observed = pyro.sample(\"position\",dist.Normal(position_final, 10.))\n",
    "    return(speed, position_observed)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The  `pyro.sample` statements are the first component of PPL, allowing us to sample from distributions. We provide pyro with unique names for each variable, so they can be tracked. We can run this model many times to find the distribution of particle speeds and observed positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "def plot_hist(hist1, hist2=None, title1=\"\", title2=\"\", true_dist=None):\n",
    "    fig = plt.figure(figsize=(10,4))\n",
    "    ax1 = plt.subplot2grid((1, 2), (0, 0));\n",
    "    ax1.hist(hist1, 25, density=True, label=\"Guide\", alpha=0.7);\n",
    "    ax1.axvline(np.mean(hist1), c=\"k\", ls=\"--\");\n",
    "    plt.title(title1);\n",
    "    if hist2 is not None:\n",
    "        ax2 = plt.subplot2grid((1, 2), (0, 1));\n",
    "        ax2.hist(hist2, 25, density=True, label=\"Guide\", alpha=0.7);\n",
    "        ax2.axvline(np.mean(hist2), c=\"k\", ls=\"--\");\n",
    "        plt.title(title2);\n",
    "    if true_dist is not None:\n",
    "        # Exact posterior distribution of speed: provide tuple for mu and sigma\n",
    "        exact_mu,exact_sigma =  true_dist\n",
    "        x = np.linspace(-5., 15., 100);\n",
    "        f = 1./(exact_sigma*np.sqrt(2*np.pi))*np.exp(-(x-exact_mu)**2/(2*exact_sigma**2))\n",
    "        for ax in fig.axes:\n",
    "            ax.plot(x, f, lw=2, label=\"Exact\", color=\"red\");\n",
    "            ax.axvline(exact_mu, ls=\"--\", color=\"red\");\n",
    "            ax.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl0AAAEICAYAAABlKUHyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3X28VWWd9/HPNwicUlERHeVBMLFJHaMJsfvOHibU0ClxJrw9VgYNDc19y11mve5wKitmLG1KrdEeMBjJSiysmVODkU89UBNyMDJByRMqIIySED4UGvS7/9gXzHazD6wFZ6+19jnf9+t1Xmfva13rWr+1OfvHb19r7bUUEZiZmZlZa72g7ADMzMzM+gMXXWZmZmYFcNFlZmZmVgAXXWZmZmYFcNFlZmZmVgAXXWZmZmYFcNFlfYakH0h6V9lxmFn/IGmUpKclDdhDn6clHVtkXFZdLrqsV0g6TdJPJW2VtFnSTySdUnZcZmb1JD0s6fepGHpM0r9KOnBfxoqItRFxYETsSGPv9sEvLV/TG7Fb+3PRZftN0sHAd4F/AQ4DhgMfB54tMy4zsx68OSIOBP4COAX4cMnxWD/host6w/EAEXFTROyIiN9HxPcj4l5J09Ks17+kWbAHJE3cuaKkIZLmStoo6VFJ/1Q/VS/pbyXdL2mLpMWSjqlbdkYab6ukawEVutdm1tYi4lHgVuAkSUdL6kwz9d2S/m5nP0kTJHVJejLNjl2V2kdLCkkDJV0OvAa4Ns2iXZv6hKTj0uMhkr4iaZOkRyR9WNIL0rJpkpZI+nTKdw9JOqvo18Ray0WX9YZfATskzZd0lqRDG5afCqwBDgc+CnxL0mFp2XxgO3Ac8ArgTOBdAJLOBf4B+BtgGPBj4Ka07HDgFmqfUA8Hfg28ulU7aGZ9j6SRwNnAz6nllvXA0cAU4BN1HxA/C3w2Ig4GXgJ8o3GsiPgQtRw1Mx1SnNlkk/8CDAGOBV4HvAN4Z93yU4HV1HLap4C5kvxhsg9x0WX7LSKeBE4DArge2JQ+MR6ZujwOXBMRf4iIm6kllb9Ky88CLo6IZyLiceBqoCOt927gkxFxf0RsBz4BjEuzXWcDqyJiYUT8AbgG+K9i9tjM2ty/SfotsAT4ITCHWg77YERsi4gVwJeBC1P/PwDHSTo8Ip6OiJ/l3WCawT8fuDQinoqIh4HP1G0D4JGIuD6dIzYfOAo4crfBrG256LJekQqjaRExAjiJ2qfFa9LiR+P5d1Z/JC0/BnghsFHSb1MS/BJwROp3DPDZumWbqR1CHJ7WX1e3/ah/bma2B+dGxCERcUxE/B9q+WRzRDxV1+cRarkGYDq10ygekLRM0pv2YZuHA4PSuM22AXUfHCPid+nhPp3kb9Xkost6XUQ8ANxArfgCGN4wRT4K2ECtSHoWODwlwEMi4uCIODH1Wwe8u27ZIRHxJxHxU2AjMHLngGn8kZiZ5bcBOEzSQXVto4BHASLiwYi4gNoHwiuBhZJe3GScaNK202+ozZgdU9e2axvWP7josv0m6c8kvV/SiPR8JHABsHMK/gjgPZJeKOk84GXAoojYCHwf+IykgyW9QNJLJL0urfdF4FJJJ6Zxh6T1Af4DOFHS30gaCLwH+NMi9tfM+paIWAf8FPikpAMknUxtdutrAJLeLmlYRPwR+G1abUeToR6jdr5Ws23soHYu2OWSDkqnSVwCfLV398aqzEWX9YanqJ0AulTSM9SKrfuA96flS4Gx1D7pXQ5MiYgn0rJ3UJtyXwVsARZSO4+BiPg2tU+VCyQ9mcY8Ky37DXAecAXwRBr/Jy3dSzPryy4ARlOb9fo28NGIuC0tmwSslPQ0tZPqOyJiW5MxPgtMSd8+/FyT5f8XeIbaF4uWAF8H5vXqXlil6fmn2pj1LknTgHdFxGllx2JmZlYmz3SZmZmZFcBFl5mZmVkBfHjRzMzMrACe6TIzMzMrwMCyA2h0+OGHx+jRo8sOwwq2evVqAF760peWHImVYfny5b+JiGFlx9EbnMP6F+cuy5O/Kld0jR49mq6urrLDsIK9/vWvB+AHP/hBqXFYOSQ9svde7cE5rH9x7rI8+cuHF83MzMwKULmZLuufPvzhD5cdgplZbs5dloeLLquE008/vewQzMxyc+6yPHx40SphxYoVrFixouwwzMxyce6yPDzTZZVw8cUXAz4Z1czai3OX5eGZLjMzM7MCuOgyMzMzK4CLLjNrK5ImSVotqVvSrCbLB0u6OS1fKml0ah8q6S5JT0u6tq7/QZJW1P38RtI1adk0SZvqlr2rqP00s77H53SZWduQNAC4DjgDWA8sk9QZEavquk0HtkTEcZI6gCuB84FtwEeAk9IPABHxFDCubhvLgW/VjXdzRMxs0S6ZWT/iostKM/2GZbseH/jqC3drA5g77ZRCY7LKmwB0R8QaAEkLgMlAfdE1GfhYerwQuFaSIuIZYImk43oaXNJY4Ajgxy2I3fqQnbmqp9wFzl+2OxddVglHjj257BCsPQwH1tU9Xw+c2lOfiNguaSswFPhNhvEvoDazFXVtb5H0WuBXwPsiYl3jSpJmADMARo0alXFXrC9w7rI8XHRZJTz24L2AE5jtlZq0xT706UkHcGHd8+8AN0XEs5L+HpgPvGG3wSPmAHMAxo8fn3VbVoJmM1KN8sxQOXdZHj6R3iph+S2fZ/ktny87DKu+9cDIuucjgA099ZE0EBgCbN7bwJJeDgyMiOU72yLiiYh4Nj29HnjlvodufZFzl+XhosvM2skyYKykMZIGUZuZ6mzo0wlMTY+nAHc2HC7syQXATfUNko6qe3oOcP8+RW1mhg8vmlkbSedozQQWAwOAeRGxUtJsoCsiOoG5wI2SuqnNcHXsXF/Sw8DBwCBJ5wJn1n3z8X8BZzds8j2SzgG2p7GmtWznzKzPc9FlZm0lIhYBixraLqt7vA04r4d1R+9h3GObtF0KXLqvsVp7ynLel9m+8OFFMzMzswJ4pssq4dS3XlJ2CGZmuTl3WR4uuqwSho46vuwQzMxyc+6yPHx40Srh0ZV38+jKu8sOw8wsF+cuy8MzXVYJv/jOPACGnzih5EjMzLJz7rI8PNNlZmZmVgAXXWZmZmYFcNFlZmZmVgAXXWZmZmYF8In0VgmvnuqLfptZ+3HusjxcdFklDDnqmLJDMDPLbU+5qzdvJzR32im9NpaVx4cXrRLWrvgxa1f8uOwwzMxyce6yPDzTZZVw3/e+BsCoca8pORIzs+ycuywPz3SZmZmZFcAzXVZpWc6J8LkOZmbWDjzTZWZmZlYAF11mZmZmBfDhRauE1/7dx8sOwcwsN+cuy8NFl1XCgUOPLDsEM7PcnLssDx9etEpYs/Q21iy9rewwzMxyce6yPDzTZZXwwF23AHDsqWeUHImZWXbOXZZHppkuSZMkrZbULWlWk+WXSFol6V5Jd0g6pm7ZVEkPpp+pvRm8mfU/GfLRYEk3p+VLJY1O7UMl3SXpaUnXNqzzgzTmivRzxJ7GMjPbF3stuiQNAK4DzgJOAC6QdEJDt58D4yPiZGAh8Km07mHAR4FTgQnARyUd2nvhm1l/kjEfTQe2RMRxwNXAlal9G/AR4AM9DP+2iBiXfh7fy1hmZrllmemaAHRHxJqIeA5YAEyu7xARd0XE79LTnwEj0uM3ArdFxOaI2ALcBkzqndDNrB/aaz5Kz+enxwuBiZIUEc9ExBJqxVdWTcfa9/DNrD/LUnQNB9bVPV+f2noyHbg1z7qSZkjqktS1adOmDCGZWT+VJafs6hMR24GtwNAMY/9rOrT4kbrCKtNYzmFmlkWWE+mbfaqLph2ltwPjgdflWTci5gBzAMaPH990bOvb3nDRFWWHYO0hS07JnLPqvC0iHpV0EHALcCHwlaxjOYf1X85dlkeWma71wMi65yOADY2dJJ0OfAg4JyKezbOu2QEHHcIBBx1SdhhWfVlyyq4+kgYCQ4DNexo0Ih5Nv58Cvk7tMOY+jWX9i3OX5ZGl6FoGjJU0RtIgoAPorO8g6RXAl6gVXI/XLVoMnCnp0HQC/Zmpzex5HlzyXR5c8t2yw7Dq22s+Ss93flN6CnBnRPQ4+yRpoKTD0+MXAm8C7tuXsaz/ce6yPPZ6eDEitkuaSa1YGgDMi4iVkmYDXRHRCfwzcCDwzXQqxNqIOCciNkv6R2qJEmB2RPhTou1mZ9Iae9qbSo7EqixjPpoL3Cipm9qsVMfO9SU9DBwMDJJ0LrUPgo8Ai1PBNQC4Hbg+rdLjWGbg3GX5ZLo4akQsAhY1tF1W9/j0Paw7D5i3rwGamdXLkI+2Aef1sO7oHoZ9ZQ/9exzLzCwv3wbIzMzMrAAuuszMzMwK4KLLzMzMrAC+4bVVwpnvu6bsEMzMcnPusjxcdFklDBx8QNkhmJnl5txlefjwolXC/Xcu5P47F5YdhplZLs5dloeLLquEh+6+nYfuvr3sMMzMcnHusjxcdJmZmZkVwOd0mZlZ25t+w7K9dzIrmWe6zMzMzArgosvMzMysAD68aJVw9qwvlh2CmVluzl2Wh2e6zMzMzArgossq4Ze3fpVf3vrVssMwM8vFucvycNFllbDuF0tY94slZYdhZpaLc5fl4aLLzMzMrAAuuszMzMwK4KLLzMzMrAC+ZIRVwoBBg8sOwcwsN+cuy8NFl1XCGy/5bNkhmJnl5txlefjwopm1FUmTJK2W1C1pVpPlgyXdnJYvlTQ6tQ+VdJekpyVdW9f/RZL+Q9IDklZKuqJu2TRJmyStSD/vKmIfzaxv8kyX5ZLlprJzp52Se9wVnXMBGHfO9NzrWv8haQBwHXAGsB5YJqkzIlbVdZsObImI4yR1AFcC5wPbgI8AJ6Wfep+OiLskDQLukHRWRNyalt0cETNbuFvWxpy7LA/PdFklbFi1jA2r9l7QWb83AeiOiDUR8RywAJjc0GcyMD89XghMlKSIeCYillArvnaJiN9FxF3p8XPAPcCIVu6E9R3OXZaHiy4zayfDgXV1z9entqZ9ImI7sBUYmmVwSYcAbwbuqGt+i6R7JS2UNLKH9WZI6pLUtWnTpmx7Ymb9josuM2snatIW+9Bn94GlgcBNwOciYk1q/g4wOiJOBm7nv2fQnj94xJyIGB8R44cNG7a3TZlZP+Wiy8zayXqgfrZpBLChpz6pkBoCbM4w9hzgwYi4ZmdDRDwREc+mp9cDr9zHuM3MXHRZNQw+cAiDDxxSdhhWfcuAsZLGpJPeO4DOhj6dwNT0eApwZ0TscaZL0j9RK84ubmg/qu7pOcD9+xG79UHOXZaHv71ovS7LNxwbTZx5ZQsisb4mIrZLmgksBgYA8yJipaTZQFdEdAJzgRsldVOb4erYub6kh4GDgUGSzgXOBJ4EPgQ8ANwjCeDaiPgy8B5J5wDb01jTCtlRaxvOXZaHiy4zaysRsQhY1NB2Wd3jbcB5Paw7uodhm50HRkRcCly6T4GamTXw4UWrhK5vXkfXN68rOwwzs1ycuywPz3RZJTz+61+WHYKZWW7OXZaHiy4zM7OKa9XdQKxYPrxoZmZmVgDPdJmZmfUBng2rPhddVgkvPvSIskMwM8vNucvyyFR0SZoEfJbadXG+HBFXNCx/LXANcDLQEREL65btAHaeabg2Is7pjcCtb3ndu2eXHYKZWW7OXZbHXosuSQOA64AzqN1eY5mkzohYVddtLbWLBn6gyRC/j4hxvRCrmZmZWdvKMtM1AejeeQNYSQuAycCuoisiHk7L/tiCGK0f+NnXrwLgVW+9pORIzMyyc+6yPLIUXcOBdXXP1wOn5tjGAZK6qN1G44qI+LfGDpJmADMARo0alWNo6ys2r/1V2SGYmeXm3GV5ZLlkRLPbY+zx5rENRkXEeOCtwDWSXrLbYBFzImJ8RIwfNmxYjqHNzMzM2kOWoms9MLLu+QhgQ9YNRMSG9HsN8APgFTniMzMzM+sTshRdy4CxksZIGgR0AJ1ZBpd0qKTB6fHhwKupOxfMzMzMrL/Y6zldEbFd0kxgMbVLRsyLiJWSZgNdEdEp6RTg28ChwJslfTwiTgReBnwpnWD/AmrndLnost0c/Kf7fi6fLwhoZmXZn9xl/U+m63RFxCJgUUPbZXWPl1E77Ni43k+BP9/PGK0fOG3aP5QdgplZbs5dloevSG9mZpWWZTbbrB34htdWCUtu+ARLbvhE2WGYmeXi3GV5eKbLKuHJ/1pbdghmZrk5d1kenukyMzMzK4CLLjMzM7MCuOgys7YiaZKk1ZK6Jc1qsnywpJvT8qWSRqf2oZLukvS0pGsb1nmlpF+mdT4nSan9MEm3SXow/T60iH00s77JRZdVwmGjjuewUceXHYZVnKQBwHXAWcAJwAWSTmjoNh3YEhHHAVcDV6b2bcBHgA80GfoL1O7/Ojb9TErts4A7ImIscEd6braLc5fl4RPprRJe9dZLyg7B2sMEoDvdVgxJC4DJPP9OF5OBj6XHC4FrJSkingGWSDqufkBJRwEHR8R/pudfAc4Fbk1jvT51nU/tVmYf7PW9srbl3GV5eKbLzNrJcGBd3fP1qa1pn4jYDmwFhu5lzPU9jHlkRGxMY20Ejmg2gKQZkrokdW3atCnjrphZf+Oiyyrhh1+6jB9+6bK9d7T+Tk3aYh/67E//3TtHzImI8RExftiwYXlWtTbn3GV5+PCiVcIzWx4vOwRrD+uBkXXPRwAbeuizXtJAYAiweS9j1t/GrH7MxyQdFREb02FI/6Ha8zh3WR6e6TKzdrIMGCtpjKRBQAfQ2dCnE5iaHk8B7oyIHmeu0mHDpyS9Kn1r8R3AvzcZa2pdu5lZbp7pMrO2ERHbJc0EFgMDgHkRsVLSbKArIjqBucCNkrqpzXB17Fxf0sPAwcAgSecCZ0bEKuB/AzcAf0LtBPpb0ypXAN+QNB1YC5zX+r00s77KRZf1C1lvmDt32iktjsT2V0QsAhY1tF1W93gbPRRHETG6h/Yu4KQm7U8AE/cjXDOzXVx0WSUc8ZI/LzsEM7PcnLssDxddtkvW2aBWGH/eRaVt28xsXzl3WR4+kd7MzMysAC66rBLuuPaD3HGtL/RtZu3Fucvy8OFFq4Rnn95adghmZrk5d1kenukyMzMzK4CLLjMzM7MCuOgyMzMzK4DP6eonyrwcRBZHn+CLkppZ+3HusjxcdFkljDtnetkhmJnl5txlefjwopmZmVkBXHRZJSy+6r0svuq9ZYdhZpaLc5fl4cOLVgk7nnu27BDMzHJz7rI8PNNlZmZmVgDPdJmZWWmq/s1qs97kmS4zMzOzAnimyyph5MtPKzsEM7Pc2i13ZZlZnDvN1x5rFRddVgl/ftbbyw7BzCw35y7Lw4cXzczMzArgmS6rhEVX/D0AZ8/6YqlxeOrdzPKoSu6y9uCZLjMzM7MCZCq6JE2StFpSt6RZTZa/VtI9krZLmtKwbKqkB9PP1N4K3Mz6pwz5aLCkm9PypZJG1y27NLWvlvTG1PZSSSvqfp6UdHFa9jFJj9YtO7uo/TSzvmevhxclDQCuA84A1gPLJHVGxKq6bmuBacAHGtY9DPgoMB4IYHlad0vvhG9m/UnGfDQd2BIRx0nqAK4Ezpd0AtABnAgcDdwu6fiIWA2Mqxv/UeDbdeNdHRGfbvW+mVnfl2WmawLQHRFrIuI5YAEwub5DRDwcEfcCf2xY943AbRGxORVatwGTeiFuM+uf9pqP0vP56fFCYKIkpfYFEfFsRDwEdKfx6k0Efh0Rj7RsD8ys38pyIv1wYF3d8/XAqRnHb7bu8MZOkmYAMwBGjRqVcWjrS8ZMOL3sEKw9ZMlHu/pExHZJW4Ghqf1nDes25qMO4KaGtpmS3gF0Ae9vNlPvHNZ/OXdZHllmutSkLTKOn2ndiJgTEeMjYvywYcMyDm19ycveMIWXvWHK3jtaf5clp/TUZ4/rShoEnAN8s275F4CXUDv8uBH4TLOgnMP6L+cuyyNL0bUeGFn3fASwIeP4+7Ou9SPbn93G9me3lR2GVV+WnLKrj6SBwBBgc4Z1zwLuiYjHdjZExGMRsSMi/ghcz+6HI62fc+6yPLIUXcuAsZLGpE+CHUBnxvEXA2dKOlTSocCZqc3seb5/9cV8/+qLyw7Dqi9LPuoEdn5TegpwZ0REau9I324cA4wF7q5b7wIaDi1KOqru6V8D9/Xanlif4Nxleez1nK50TsRMasXSAGBeRKyUNBvoiohOSadQ+7bPocCbJX08Ik6MiM2S/pFaogSYHRGbW7QvZtbHZclHwFzgRknd1Ga4OtK6KyV9A1gFbAcuiogdAJJeRO0bke9u2OSnJI2jdhjy4SbLzcwyy3RF+ohYBCxqaLus7vEyalP1zdadB8zbjxjNzHbJkI+2Aef1sO7lwOVN2n9H7WT7xvYL9zdeM7OdfEV6MzMzswK46DIzMzMrgG94bZUw9rQ3lR2CmVluzl2Wh4suqwQnLjNrR85dlocPL1olbHvqt2x76rdlh2Fmlotzl+XhmS6rhDuvmwXA2bO+WHIkZmbZ9cXcNf2GZXvvBMyddkqLI+l7PNNlZmZmVgAXXWZmZmYFcNFlZmZmVgAXXWZmZmYF8In0Vgl/9pdvKTsEM7PcnLssDxddVgnHnnpG2SGYWS/L+i24dubcZXn48KJVwtNPPMbTTzxWdhhmZrk4d1keLrqsEn50/Uf50fUfLTsMM7NcnLssDxddZmZmZgVw0WVmZmZWAJ9Ib2ZmufSHE+TNWsEzXWZmZmYF8EyXVcJJk95WdghmZrk5d1keLrqsEkaNe03ZIZiZ5ebcZXn48KJVwtaNj7B14yNlh2FtQNIkSasldUua1WT5YEk3p+VLJY2uW3Zpal8t6Y117Q9L+qWkFZK66toPk3SbpAfT70NbvX/WXpy7LA8XXVYJP5n/SX4y/5Nlh2EVJ2kAcB1wFnACcIGkExq6TQe2RMRxwNXAlWndE4AO4ERgEvD5NN5OfxkR4yJifF3bLOCOiBgL3JGem+3i3GV5uOgys3YyAeiOiDUR8RywAJjc0GcyMD89XghMlKTUviAino2Ih4DuNN6e1I81Hzi3F/bBzPopF11m1k6GA+vqnq9PbU37RMR2YCswdC/rBvB9Scslzajrc2REbExjbQSOaBaUpBmSuiR1bdq0aZ92zMz6Pp9Ib5ZTlmsUzZ12SgGR9Etq0hYZ++xp3VdHxAZJRwC3SXogIn6UNaiImAPMARg/fnxjPGZmgIuuPsEXKrR+ZD0wsu75CGBDD33WSxoIDAE272ndiNj5+3FJ36Z22PFHwGOSjoqIjZKOAh7v/V0ys/7CRZdVwsvf/Ldlh2DtYRkwVtIY4FFqJ8a/taFPJzAV+E9gCnBnRISkTuDrkq4CjgbGAndLejHwgoh4Kj0+E5jdMNYV6fe/t3TvrO04d1keLrqsEoafuLfzmc1q52hJmgksBgYA8yJipaTZQFdEdAJzgRsldVOb4epI666U9A1gFbAduCgidkg6Evh27Vx7BgJfj4jvpU1eAXxD0nRgLXBeYTtrbcG5y/Jw0WWV8MTaXwEwdNTxJUdiVRcRi4BFDW2X1T3eRg/FUURcDlze0LYGeHkP/Z8AJu5nyNaHOXdZHv72olXC0q9fxdKvX1V2GGZmuTh3WR4uuszMzMwK4KLLzMzMrAAuuszMzMwK4KLLzMzMrAD+9qJVwivf8n/KDsHMLDfnLssj00yXpEmSVkvqljSryfLBkm5Oy5dKGp3aR0v6vaQV6eeLvRu+9RVHjj2ZI8eeXHYYZma5OHdZHnud6ZI0ALgOOIPabTSWSeqMiFV13aYDWyLiOEkdwJXA+WnZryNiXC/HbX3MYw/eC+DkZWZtxbnL8sgy0zUB6I6INRHxHLAAmNzQZzIwPz1eCExUuryzWRbLb/k8y2/5fNlhmJnl4txleWQpuoYD6+qer09tTftExHZgKzA0LRsj6eeSfijpNc02IGmGpC5JXZs2bcq1A2ZmZmbtIEvR1WzGKjL22QiMiohXAJdQu9nswbt1jJgTEeMjYvywYcMyhGRmZmbWXrIUXeuBkXXPRwAbeuojaSAwBNgcEc+me5cREcuBXwO+QZWZmZn1O1mKrmXAWEljJA0COoDOhj6dwNT0eApwZ0SEpGHpRHwkHQuMBdb0TuhmZmZm7WOv316MiO2SZgKLgQHAvIhYKWk20BURncBc4EZJ3cBmaoUZwGuB2ZK2AzuAv4+Iza3YEWtvp771krJDMDPLzbnL8sh0cdSIWAQsami7rO7xNuC8JuvdAtyynzFaPzB0lI86m1n7ce6yPHxF+gqbfsOyskMozKMr7wZg+IkTSo7EzCw75y7Lw0WXVcIvvjMPcOIys/bi3GV5+IbXZmZmZgXwTJdZC2Q5NDx32ikFRGJmZlXhmS4zMzOzArjoMrO2ImmSpNWSuiXNarJ8sKSb0/KlkkbXLbs0ta+W9MbUNlLSXZLul7RS0nvr+n9M0qOSVqSfs4vYRzPrm3x40Srh1VMvLTsEawPpYsvXAWdQuxPGMkmdEbGqrtt0YEtEHCepA7gSOF/SCdSuIXgicDRwu6Tjge3A+yPiHkkHAcsl3VY35tUR8eli9tDajXOX5eGZLquEIUcdw5Cjjik7DKu+CUB3RKyJiOeABcDkhj6Tgfnp8UJgoiSl9gXp9mQPAd3AhIjYGBH3AETEU8D9wPAC9sX6AOcuy8NFl1XC2hU/Zu2KH5cdhlXfcGBd3fP17F4g7eoTEduBrcDQLOumQ5GvAJbWNc+UdK+keZIObRaUpBmSuiR1bdq0Ke8+WRtz7rI8XHRZJdz3va9x3/e+VnYYVn1q0hYZ++xxXUkHUruDxsUR8WRq/gLwEmAcsBH4TLOgImJORIyPiPHDhg3b8x5Yn+LcZXm46DKzdrIeGFn3fASwoac+kgYCQ6jdE7bHdSW9kFrB9bWI+NbODhHxWETsiIg/AtdTO7xpZrZPXHSZWTtZBoyVNEbSIGonxnc29OkEpqbHU4A7IyJSe0f6duMYYCxwdzrfay5wf0RcVT+QpKPqnv41cF+v75GZ9Rv+9qKZtY2I2C5pJrAYGADMi4iVkmYDXRHRSa2AulFSN7UZro607kpJ3wBWUfvG4kURsUPSacCFwC8lrUib+oeIWAR8StI4aochHwYhK7MnAAAJFklEQVTeXdjOmlVcb90fuD9dKNpFl1lJfNX6fZOKoUUNbZfVPd4GnNfDupcDlze0LaH5+V5ExIX7G6+Z2U4uuqwSXvt3Hy87BDOz3Jy79l9/+gDqoqskvTUt21ccOPTIskMwM8vNucvycNFllbBm6W0AHHvqGSVHYtae+tNsQZU4d1keLrqsEh646xbAicvM2otzl+XhS0aYmZmZFcAzXWZm/YTPJTUrl2e6zMzMzArgma4W8KdJMzMza+SiyyrhDRddUXYIZma5OXdZHi66rBIOOOiQskMwM8vNucvy8DldVgkPLvkuDy75btlhmJnl4txlebjoskpw4jKzduTcZXm46DIzMzMrgM/pMquwrN+E9e1dzMyqzzNdZmZmZgXwTJeZWcX52n9me9eb75NWHT1w0WWVcOb7rik7BDOz3Jy7LA8XXVYJAwcfUHYIZoXzDFb7c+6yPFx0WSXcf+dCAF72hiklR9Kesvzn7ZPtzXqfc5fl4aIr8SfOcj109+2AE5eZtRfnLsvDRZeZmZlVWl+ZGMlUdEmaBHwWGAB8OSKuaFg+GPgK8ErgCeD8iHg4LbsUmA7sAN4TEYt7LfqM+so/llmrtcO3f1qRj3oaU9IYYAFwGHAPcGFEPNeSHTOzPm+vRZekAcB1wBnAemCZpM6IWFXXbTqwJSKOk9QBXAmcL+kEoAM4ETgauF3S8RGxo7d2wAWVWTZ94b3SinyU1ulpzCuBqyNigaQvprG/0Po9NbO+KMvFUScA3RGxJn3CWwBMbugzGZifHi8EJkpSal8QEc9GxENAdxrPzGxftCIfNR0zrfOGNAZpzHNbuG9m1sdlObw4HFhX93w9cGpPfSJiu6StwNDU/rOGdYc3bkDSDGBGevq0pNWZot/d4cBv9nHd3uDt7+f2571zv2rytt//Nt120+3Pe2eu9Y/J2K9V+ajZmEOB30bE9ib9nydjDiv732gnx9GCGNo8d1UlBqhGHIfPe2euGLLmr0xFl5q0RcY+WdYlIuYAczLEskeSuiJi/P6O4+17++22/X60763IR81m/DPnL8iWw8r+N3Ic1YyhKnFUIYaqxNHKGLIcXlwPjKx7PgLY0FMfSQOBIcDmjOuamWXVinzUU/tvgEPSGD1ty8wssyxF1zJgrKQxkgZROxG1s6FPJzA1PZ4C3BkRkdo7JA1O3wIaC9zdO6GbWT/UinzUdMy0zl1pDNKY/97CfTOzPm6vhxfTOREzgcXUvk49LyJWSpoNdEVEJzAXuFFSN7VPlB1p3ZWSvgGsArYDF/XmNxeb2O9DlN6+t9+m2+8X+96qfNRszLTJDwILJP0T8PM09r4q+99oJ8fx36oQA1QjjirEANWIo2UxqPZhzszMzMxaKcvhRTMzMzPbTy66zMzMzArQ54ouSf8s6QFJ90r6tqRDCtruJEmrJXVLmlXENtN2R0q6S9L9klZKem9R226IY4Ckn0v6bgnbPkTSwvTvfr+k/1Hw9t+XXvv7JN0k6YAWb2+epMcl3VfXdpik2yQ9mH4fWvD2S3nftYM9vTaSLk05Y7WkN7Y4jkrlqCL/ZutieV6eSl+eWJpiuDl9kaLVMeyWr0p6LXbLW61+PfLkLtV8Lv293ivpL1ocRyHv0z5XdAG3ASdFxMnAr4BLW71B/fetSc4CTgAuUO2WI0XYDrw/Il4GvAq4qMBt13svcH8J24XaPfO+FxF/Bry8yDgkDQfeA4yPiJOonYjd0eLN3gBMamibBdwREWOBO9LzIrdf+PuujTR9bfT82xJNAj6fckmvq2iOKvJvdqfGPLXzNk9jgS3UbvPUas3yVaGvxR7yVqtfjxvInrvOovYN47HULjzcm7ffahZHIe/TPld0RcT3664g/TNq19ZptSy3JmmJiNgYEfekx09RewM3vWp2q0gaAfwV8OUit5u2fTDwWtK3yiLiuYj4bcFhDAT+RLXrOb2IFl/LKSJ+RO1befXqb33T0tvVNNt+Se+7trCH16bI26RVMUcV9jcLu+cpqfjbPO0hXxX6WiSNeWsjLX49cuauycBXouZn1K6Zd1Sr4ijqfdrniq4GfwvcWsB2mt2apNDCB0DSaOAVwNKCN30N8P+APxa8XYBjgU3Av6bDBl+W9OKiNh4RjwKfBtZSS1pbI+L7RW2/zpERsTHFtBE4ooQYdirqfdeO6l+bIvNGFXNU0X+zjXkq822eelFP+arQ16JZ3gKWU/zrAT3ve5l/sy17n7Zl0SXp9nQcuvFncl2fD1Gb1v5aESE1aSv0WhySDgRuAS6OiCcL3O6bgMcjYnlR22wwEPgL4AsR8QrgGYo5TAFAOv9gMjAGOBp4saS3F7X9qin4fVcZ+5iTiswb/TZHpW03y1NlvCal5qudmuUtaofzGpV5TalS/mZb/T7Ncu/FyomI0/e0XNJU4E3AxCjmQmSl3u5I0gupJbOvRcS3itpu8mrgHElnAwcAB0v6akQUVXisB9ZHxM7ZvYUUm8ROBx6KiE0Akr4F/E/gqwXGAPCYpKMiYmOagn+84O2X8b6rjH3MSUXmjSrmqCL/ZnfLU9Rmvg6RNDDN7hTxmvSUr4p+//aUt4p+PaDnfS/8b7aI92lbznTtiaRJ1K4ifU5E/K6gzWa5NUlLpPMS5gL3R8RVRWyzXkRcGhEjImI0tf2+s8CCi4j4L2CdpJemponUrjhelLXAqyS9KP1bTKScLxTU3/qm8NvVlPS+awt7eG2KvE1aFXNUYX+zPeSpt1HwbZ72kK+Kfv82y1urKOe2Vz3teyfwjvQtxldRO3VjY6uCKOx9GhF96ofaSW7rgBXp54sFbfdsat94+DXwoQL39zRqU5331u3z2SW99q8HvlvCdscBXek1+Dfg0IK3/3HgAeA+4EZgcIu3dxO18zD+QO1T2HRq56fcATyYfh9W8PZLed+1w8+eXhvgQylnrAbOanEclcpRRf7NNsSzK09RO8fq7vRv9M1Wv3fTNnfLV2W8Fs3yVqtfjzy5i9phvevS3+svqX3TspVxFPI+9W2AzMzMzArQ5w4vmpmZmVWRiy4zMzOzArjoMjMzMyuAiy4zMzOzArjoMjMzMyuAiy4zMzOzArjoMjMzMyvA/wd6q2T81Sd2DgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "uncond = [propagate_cloud(guess=5., guess_sd=2., time=10.) for i in range(1000)]\n",
    "speed_uncond, pos_uncond = zip(*uncond)\n",
    "plot_hist(speed_uncond, pos_uncond, \"Speed\", \"Position\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2) Conditioning\n",
    "\n",
    "What if we observe the particle at a certain point, x=30 m? We want to use this information in our model to learn more about the latent variable, speed. This is when we condition on our observations.\n",
    "In pyro, this is done with the pyro.condition statement. This statement returns a new function almost identical to the original, but overwrites the   `pyro.sample(\"position\",...)` by fixing it at the observed value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "position_observed = 30\n",
    "conditioned_propagate = pyro.condition(propagate_cloud,\n",
    "                                       data={\"position\": position_observed})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAEICAYAAAATE/N5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3X28VWWZ//HPN1BswmfRUR5Ekyy1ogTsl2aOT6Fj4kw6YaZSNFQTv9Gp5heaaTJm2pTapJNSEFo+lWYxhYPPNWYihyQV0UBUOECCYT5VGnr9/lj3mbbbfTzrwDnr3vuc7/v12q+z9r3utfa1D+zrXPte91pLEYGZmZmZ5fG63AGYmZmZ9WcuxszMzMwycjFmZmZmlpGLMTMzM7OMXIyZmZmZZeRizMzMzCwjF2PWL0i6Q9LHcsdhZv2DpBGSnpM04DX6PCdp9yrjsubkYsx6naQDJN0l6WlJ6yX9QtLY3HGZmdWS9JikP6Yi6QlJ35E0eGP2FRErImJwRLyU9v2qL4Rp/fKeiN1am4sx61WStgJ+AnwD2A4YCpwNvJAzLjOzTrw/IgYD7wTGAmdkjsf6ARdj1tveBBARV0fESxHxx4i4KSLukzQpjZJ9I42aPSTpkI4NJW0taaakNZJWSTqndshf0kclLZH0lKR5knatWXdY2t/Tki4GVOm7NrOWFhGrgBuBfSTtImlOGtlfJukfO/pJGiepTdIzaTTtgtQ+UlJIGijpS8B7gIvTqNvFqU9I2iMtby3pCknrJD0u6QxJr0vrJkm6U9JXU757VNIRVf9OrPe4GLPe9hvgJUmXSzpC0rZ16/cDlgM7AGcBP5S0XVp3ObAB2AN4B3A48DEASccApwN/DwwB/ge4Oq3bAbie4hvtDsAjwP699QbNrO+RNBw4EriXIre0A7sAxwLn1nxx/Drw9YjYCngj8P36fUXE5yly1NR0aHJqg5f8BrA1sDvwXuAk4CM16/cDHqbIaV8BZkryl8w+wsWY9aqIeAY4AAjgW8C69A1zp9RlLXBRRPw5Iq6lSDZ/m9YfAZwaEc9HxFrgQmBi2u7jwJcjYklEbADOBUan0bEjgQcj4rqI+DNwEfDbat6xmbW4H0n6PXAn8DNgBkUO+1xE/CkiFgHfBk5M/f8M7CFph4h4LiLu7u4LphH/DwKnRcSzEfEY8LWa1wB4PCK+leagXQ7sDOz0qp1ZS3IxZr0uFUyTImIYsA/Ft8uL0upV8cq71T+e1u8KbAaskfT7lBwvA3ZM/XYFvl6zbj3FocihafuVNa8ftc/NzF7DMRGxTUTsGhH/RJFP1kfEszV9HqfINQCTKaZjPCRpgaSjNuI1dwA2T/tt9BpQ84UyIv6QFjfq5AJrPi7GrFIR8RAwm6IoAxhaN9Q+AlhNUTy9AOyQEuM2EbFVROyd+q0EPl6zbpuIeH1E3AWsAYZ37DDtfzhmZt23GthO0pY1bSOAVQARsTQijqf4ong+cJ2kNzTYTzRo6/AkxQjbrjVt//sa1ve5GLNeJenNkj4jaVh6Phw4HugYyt8R+GdJm0k6DngLMDci1gA3AV+TtJWk10l6o6T3pu0uBU6TtHfa79Zpe4CfAntL+ntJA4F/Bv66ivdrZn1LRKwE7gK+LGkLSW+jGA27EkDShyUNiYiXgd+nzV5qsKsnKOaDNXqNlyjmmn1J0pZpusWnge/17LuxZuVizHrbsxQTT+dLep6iCHsA+ExaPx8YRfHN8EvAsRHxu7TuJIqh+weBp4DrKOZJEBE3UHwLvUbSM2mfR6R1TwLHAecBv0v7/0Wvvksz68uOB0ZSjJLdAJwVETendeOBxZKeo5jMPzEi/tRgH18Hjk1nQ/5Hg/X/F3ie4oSmO4GrgFk9+i6saemV03XMqiNpEvCxiDggdyxmZma5eGTMzMzMLCMXY2ZmZmYZ+TClmZmZWUYeGTMzMzPLaGDuALpjhx12iJEjR+YOwzbBww8/DMCee+6ZORJrBQsXLnwyIobkjqMnOH/1T855/VvZHNZSxdjIkSNpa2vLHYZtgoMOOgiAO+64I2sc1hokPd51r9bg/NU/Oef1b2VzmA9TmpmZmWXUUiNj1vrOOOOM3CGYmVXGOc/KcDFmlTr00ENzh2BmVhnnPCvDhymtUosWLWLRokW5wzAzq4RznpXhkTGr1Kmnngp4MquZ9Q/OeVaGR8bMzMzMMnIxZmbWCUnDJd0uaYmkxZJOSe1flLRK0qL0ODJ3rGbWunyY0syscxuAz0TEryRtCSyUdHNad2FEfDVjbGbWR7gYMzPrRESsAdak5WclLQGG5o3KzPoaF2PWpcmzF5TqN3PS2C77nHvuuZsajlkWkkYC7wDmA/sDUyWdBLRRjJ491WCbKcAUgBEjRlQWq+VTny8H739iw/Yy+dL6D88Zs0q9+93v5t3vfnfuMMy6RdJg4Hrg1Ih4Bvgm8EZgNMXI2dcabRcRMyJiTESMGTKkT9xi07ppp1FvY6dRb8sdhjU5F2NWqbvuuou77rordxhmpUnajKIQuzIifggQEU9ExEsR8TLwLWBczhiteT2x9D6eWHpf7jCsyfkwpVXq9NNPB3zNHWsNkgTMBJZExAU17Tun+WQAfwc8kCM+a34Lr/9PAI6cdmnmSKyZuRgzM+vc/sCJwP2SOi6jfjpwvKTRQACPAR/PE56Z9QUuxszMOhERdwJqsGpu1bGYWd/lOWNmZmZmGbkYMzMzM8uo1GFKSeOBrwMDgG9HxHl16z8NfIziatXrgI9GxONp3cnAGanrORFxeWrfF5gNvJ5iyP+UiIhNfUPW3C666KLcIZiZVWa/D306dwjWArosxiQNAC4BDgPagQWS5kTEgzXd7gXGRMQfJH0S+ArwQUnbAWcBYygmui5M2z5FcZ2eKcDdFMXYeODGnntr1oxGjx6dOwQzs8psP+JNuUOwFlDmMOU4YFlELI+IF4FrgAm1HSLi9oj4Q3p6NzAsLb8PuDki1qcC7GZgvKSdga0i4pdpNOwK4JgeeD/W5G655RZuueWW3GGYmVVi1eJ7WLX4ntxhWJMrc5hyKLCy5nk7sN9r9J/MX0a4Gm07ND3aG7S/im8n0recc845ABx66KGZIzEz632//q9ZAAzd29cFts6VGRlrdFp3w7ldkj5McUjy37vYtvQ+fTsRMzMz68vKFGPtwPCa58OA1fWdJB0KfB44OiJe6GLbdv5yKLPTfZqZmZn1dWWKsQXAKEm7SdocmAjMqe0g6R3AZRSF2NqaVfOAwyVtK2lb4HBgXrqNyLOS3pVuN3IS8OMeeD9mZmZmLaXLOWMRsUHSVIrCagAwKyIWS5oOtEXEHIrDkoOBHxS1FSsi4uiIWC/p3ygKOoDpEbE+LX+Sv1za4kZ8JqWZmZn1Q6WuMxYRc6m7/UdEnFmz3Ols7IiYBcxq0N4G7FM6UusTLrvsstwhmJlVZv+TT8sdgrUA35vSKrXnnnvmDsHMrDJb77xr7hCsBbgYs0odduoFAIwY/Z5O+8ycNLaqcMzMetWKRf8DvHbOM3MxZpV64L+vBJyYzKx/cM6zMnyjcDMzM7OMXIyZmZmZZeRizMzMzCwjF2NmZmZmGXkCv1XqwH88O3cIZmaVcc6zMlyMWaUGb79T7hDMzCrjnGdl+DClVWr5/JtZPv/m3GGYmVXCOc/K8MiYVeqh268HYPf9DssciZlZ73POszI8MmZmZmaWkYsxMzMzs4xcjJmZmZll5GLMzMzMLKNSE/gljQe+DgwAvh0R59WtPxC4CHgbMDEirkvtfwNcWNP1zWn9jyTNBt4LPJ3WTYqIRZvwXqwFHPyp87ruZGbWRzjnWRldFmOSBgCXAIcB7cACSXMi4sGabiuAScBna7eNiNuB0Wk/2wHLgJtquvxrR+Fm/cMWW26TOwQzs8o451kZZUbGxgHLImI5gKRrgAnA/xZjEfFYWvfya+znWODGiPjDRkdrLW/pnT8BYNQBR2WOxMys9znnWRllirGhwMqa5+3AfhvxWhOBC+raviTpTOBWYFpEvFC/kaQpwBSAESNGbMTLWlUmz17QZR8nJjPrT5zzrIwyE/jVoC268yKSdgbeCsyraT6NYg7ZWGA74HONto2IGRExJiLGDBkypDsva2ZmZtb0yoyMtQPDa54PA1Z383X+AbghIv7c0RARa9LiC5K+Q918M+u/yoywzZw0toJIzMzMel+ZkbEFwChJu0nanOJw45xuvs7xwNW1DWm0DEkCjgEe6OY+zcx6laThkm6XtETSYkmnpPbtJN0saWn6uW3uWM2sdXU5MhYRGyRNpTjEOACYFRGLJU0H2iJijqSxwA3AtsD7JZ0dEXsDSBpJMbL2s7pdXylpCMVh0EXAJ3roPVk/4NEzq8gG4DMR8StJWwILJd1Mcfb4rRFxnqRpwDQ6mWphZtaVUtcZi4i5wNy6tjNrlhdQHL5stO1jFCcB1Lcf3J1ArW84/F8uyh2CWWlpOsWatPyspCUU+WwCcFDqdjlwBy7GrAHnPCujVDFm1lMGDtoidwhmGyWN8r8DmA/s1DHvNSLWSNqxk218Nng/55xnZfh2SFapJbddx5LbfJ1fay2SBgPXA6dGxDNlt/PZ4OacZ2W4GLNKPXrPLTx6zy25wzArTdJmFIXYlRHxw9T8RM1JSDsDa3PFZ83NOc/KcDFmZtaJdLb3TGBJRNRetHoOcHJaPhn4cdWxmVnf4TljZmad2x84Ebhf0qLUdjpwHvB9SZMp7s17XKb4zKwPcDFmZtaJiLiTxnchATikyljMrO/yYUozMzOzjDwyZpU6ctqluUMwM6uMc56V4ZExMzMzs4xcjFml7r/xe9x/4/dyh2FmVgnnPCvDxZhVauWv72Tlr+/MHYaZWSWc86wMF2NmZmZmGbkYMzMzM8vIxZiZmZlZRr60hVVqwOaDcodgZlYZ5zwrw8WYVep9n/567hDMzCrjnGdllDpMKWm8pIclLZM0rcH6AyX9StIGScfWrXtJ0qL0mFPTvpuk+ZKWSrpW0uab/nbMzMzMWkuXxZikAcAlwBHAXsDxkvaq67YCmARc1WAXf4yI0elxdE37+cCFETEKeAqYvBHxW4tZNGcmi+bMzB2GmVklnPOsjDIjY+OAZRGxPCJeBK4BJtR2iIjHIuI+4OUyLypJwMHAdanpcuCY0lFby1r94AJWP7ggdxhmZpVwzrMyyhRjQ4GVNc/bU1tZW0hqk3S3pI6Ca3vg9xGxoat9SpqStm9bt25dN17WzMzMrPmVmcCvBm3RjdcYERGrJe0O3CbpfuCZsvuMiBnADIAxY8Z053X7vcmz/W3MzMys2ZUZGWsHhtc8HwasLvsCEbE6/VwO3AG8A3gS2EZSRzHYrX2amZmZ9RVlirEFwKh09uPmwERgThfbACBpW0mD0vIOwP7AgxERwO1Ax5mXJwM/7m7w1noGDd6aQYO3zh2GmVklnPOsjC4PU0bEBklTgXnAAGBWRCyWNB1oi4g5ksYCNwDbAu+XdHZE7A28BbhM0ssUhd95EfFg2vXngGsknQPcC/h0k37gkKnn5w7BzKwyznlWRqmLvkbEXGBuXduZNcsLKA411m93F/DWTva5nOJMTTMzM7N+y/emtEq1/eAS2n5wSe4wzMwq4ZxnZfh2SFaptY/cnzsEM7PKOOdZGR4ZMzMzM8vIxZiZmZlZRi7GzMzMzDLynDGr1Bu23TF3CGZmlXHOszJcjFml3vvx6blDMDOrjHOeleHDlGZmZmYZuRizSt191QXcfdUFucMwM6uEc56V4cOUVqn1K36TOwQzs8o451kZHhkzMzMzy8jFmJmZmVlGLsbMzMzMMnIxZpXa6q9HsNVfj8gdhllpkmZJWivpgZq2L0paJWlRehyZM0ZrXs55VoYn8FulDph0eu4QzLprNnAxcEVd+4UR8dXqw7FW4pxnZZQaGZM0XtLDkpZJmtZg/YGSfiVpg6Rja9pHS/qlpMWS7pP0wZp1syU9WvPNcnTPvCUzs54TET8H1ueOw8z6ri6LMUkDgEuAI4C9gOMl7VXXbQUwCbiqrv0PwEkRsTcwHrhI0jY16/81Ikanx6KNfA/WQu6cfS53zj43dxhmPWFq+pI5S9K2uYOx5uScZ2WUGRkbByyLiOUR8SJwDTChtkNEPBYR9wEv17X/JiKWpuXVwFpgSI9Ebi3pmd+u4Jnfrsgdhtmm+ibwRmA0sAb4WqNOkqZIapPUtm7duirjsybhnGdllCnGhgIra563p7ZukTQO2Bx4pKb5S+mb5YWSBnWynZOZmTWViHgiIl6KiJeBb1F8aW3Ub0ZEjImIMUOG+HuomTVWphhTg7bozotI2hn4LvCRlLwATgPeDIwFtgM+12hbJzMzazYpp3X4O+CBzvqamXWlzNmU7cDwmufDgNVlX0DSVsBPgTMi4u6O9ohYkxZfkPQd4LNl92lmVhVJVwMHATtIagfOAg5KJx0F8Bjw8WwBmlnLK1OMLQBGSdoNWAVMBD5UZueSNgduAK6IiB/Urds5ItZIEnAM/mbZL2w34k25QzDrlog4vkHzzMoDsZbknGdldFmMRcQGSVOBecAAYFZELJY0HWiLiDmSxlIUXdsC75d0djqD8h+AA4HtJU1Ku5yUzpy8UtIQisOgi4BP9PSbs+bzrg99OncIZmaVcc6zMkpd9DUi5gJz69rOrFleQHH4sn677wHf62SfB3crUjMzM7M+yLdDskr97LIz+dllZ3bd0cysD3DOszJ8OySr1PNPrc0dgplZZZzzrAyPjJmZmZll5JEx67Mmz15Qqt/MSWN7ORIzM7POeWTMzMzMLCOPjFmldnzjW3OHYGZWGec8K8PFmFVqzHGfyh2CmVllnPOsDB+mNDMzM8vIxZhV6taLP8etFze8J7yZWZ/jnGdl+DClVeqF557OHYKZWWWc86wMj4yZmZmZZeRizMzMzCwjF2NmZmZmGXnOmFVql718tXsz6z+c86wMF2NWqdFHT84dwquUuW2Sb5lkZhujGXOeNR8fpjQzMzPLqFQxJmm8pIclLZM0rcH6AyX9StIGScfWrTtZ0tL0OLmmfV9J96d9/ockbfrbsWY374JTmHfBKbnDMDOrhHOeldFlMSZpAHAJcASwF3C8pL3quq0AJgFX1W27HXAWsB8wDjhL0rZp9TeBKcCo9Bi/0e/CWsZLL77ASy++kDsMM7NKOOdZGWVGxsYByyJieUS8CFwDTKjtEBGPRcR9wMt1274PuDki1kfEU8DNwHhJOwNbRcQvIyKAK4BjNvXNmJmZmbWaMsXYUGBlzfP21FZGZ9sOTctd7lPSFEltktrWrVtX8mXNzMzMWkOZYqzRXK4ouf/Oti29z4iYERFjImLMkCFDSr6smZmZWWsoc2mLdmB4zfNhwOqS+28HDqrb9o7UPmwj92ktbPjbD8gdgplZZZzzrIwyxdgCYJSk3YBVwETgQyX3Pw84t2bS/uHAaRGxXtKzkt4FzAdOAr7RvdCtFb31iA/nDsHMrDLOeVZGl4cpI2IDMJWisFoCfD8iFkuaLuloAEljJbUDxwGXSVqctl0P/BtFQbcAmJ7aAD4JfBtYBjwC3Nij78zMzMysBZS6An9EzAXm1rWdWbO8gFcedqztNwuY1aC9DdinO8Fa65t73icAOHLapZkjMTPrfc55VoavwG9mZmaWkYsxMzMzs4xcjJmZmZll5GLMzMzMLKNSE/jNespu4w7NHYKZWWWc86wMF2NWqbccfGzuEMzMKuOcZ2W4GLNKbXjhTwAMHLRF5ki6Z/LsBV32mTlpbAWRWNUkzQKOAtZGxD6pbTvgWmAk8BjwDxHxVK4YrXm1as6zannOmFXqpgtP5aYLT80dhll3zAbG17VNA26NiFHArem52as451kZHhlrUWVGasxs00XEzyWNrGuewF/uu3s5xT13P1dZUGbWp3hkzMys+3aKiDUA6eeOjTpJmiKpTVLbunXrKg3QzFqHizEzs14SETMiYkxEjBkyZEjucMysSbkYMzPrvick7QyQfq7NHI+ZtTDPGbNKjTrgqNwhmPWEOcDJwHnp54/zhmPNyjnPynAx1mT6+sR8JyZrNZKuppisv4OkduAsiiLs+5ImAyuA4/JFaM3MOc/KcDFmlfrTs78HYIstt8kciVk5EXF8J6sOqTQQa0nOeVZGqTljksZLeljSMkmvup6OpEGSrk3r53ecBi7pBEmLah4vSxqd1t2R9tmxruHZSNa33HbJNG67xJdkMrP+wTnPyuiyGJM0ALgEOALYCzhe0l513SYDT0XEHsCFwPkAEXFlRIyOiNHAicBjEbGoZrsTOtZHhCfAmpmZWb9TZmRsHLAsIpZHxIvANRQXPKw1geLChwDXAYdIUl2f44GrNyVYMzMzs76mTDE2FFhZ87w9tTXsExEbgKeB7ev6fJBXF2PfSYcov9CgeAN80UQzMzPr28oUY42KpOhOH0n7AX+IiAdq1p8QEW8F3pMeJzZ6cV800czMzPqyMmdTtgPDa54PA1Z30qdd0kBga2B9zfqJ1I2KRcSq9PNZSVdRHA69olvRW8t58998IHcIZmaVcc6zMsoUYwuAUZJ2A1ZRFFYfquvTcQHEXwLHArdFRABIeh3FNXgO7OicCrZtIuJJSZsBRwG3bOJ7sRaw+36H5Q7BzKwyznlWRpfFWERskDQVmAcMAGZFxGJJ04G2iJgDzAS+K2kZxYjYxJpdHAi0R8TymrZBwLxUiA2gKMS+1SPvyJrac797AoDB2++UORIzs97nnGdllLroa0TMBebWtZ1Zs/wnOrkCdUTcAbyrru15YN9uxmp9wM+/dRYAR067NHMkZma9zznPyvCNws3MzMwycjFmZmZmlpGLMTMzM7OMXIyZmZmZZVRqAr9ZT9ln/Am5QzAzq4xznpXhYswqNWL0e3KHYGZWGec8K8OHKa1ST695nKfXPJ47DDOzSjjnWRkuxqxSv7j8y/zi8i/nDsPMrBLOeVaGizEzMzOzjFyMmZmZmWXkYszMzMwsIxdjZmZmZhn50hZWqbe//6O5QzAzq4xznpXhYswqNXTvcblDMDOrjHOeleHDlFap3634Db9b8ZvcYZiZVcI5z8pwMWaVmn/VBcy/6oLcYZiZVcI5z8ooVYxJGi/pYUnLJE1rsH6QpGvT+vmSRqb2kZL+KGlRelxas82+ku5P2/yHJPXUmzIzMzNrFV0WY5IGAJcARwB7AcdL2quu22TgqYjYA7gQOL9m3SMRMTo9PlHT/k1gCjAqPcZv/NswMzMza01lRsbGAcsiYnlEvAhcA0yo6zMBuDwtXwcc8lojXZJ2BraKiF9GRABXAMd0O3ozMzOzFlemGBsKrKx53p7aGvaJiA3A08D2ad1uku6V9DNJ76np397FPgGQNEVSm6S2devWlQjXzMzMrHWUubRFoxGuKNlnDTAiIn4naV/gR5L2LrnPojFiBjADYMyYMQ37WOvY9wP/lDsEM7PKOOdZGWWKsXZgeM3zYcDqTvq0SxoIbA2sT4cgXwCIiIWSHgHelPoP62Kf1gftNOptuUPoNZNnL+iyz8xJYyuIxMyaRV/OedZzyhymXACMkrSbpM2BicCcuj5zgJPT8rHAbRERkoakEwCQtDvFRP3lEbEGeFbSu9LcspOAH/fA+7Em98TS+3hi6X25wzAzq4RznpXR5chYRGyQNBWYBwwAZkXEYknTgbaImAPMBL4raRmwnqJgAzgQmC5pA/AS8ImIWJ/WfRKYDbweuDE9rI9beP1/AnDktEu76GnW/CQ9BjxLkd82RMSYvBFZs3HOszJK3Q4pIuYCc+vazqxZ/hNwXIPtrgeu72SfbcA+3QnWzKwJ/U1EPJk7CDNrXb4Cv5mZmVlGLsbMzDZeADdJWihpSv1KX5rHzMpwMWZmtvH2j4h3Utyh5FOSDqxdGREzImJMRIwZMmRIngjNrOmVmjNm1lP2+9Cnc4dg1mMiYnX6uVbSDRR3LPl53qismTjnWRkuxqxS2494U+4QzHqEpDcAr4uIZ9Py4cD0zGFZk3HOszJcjFmlVi2+B4Che4/LHInZJtsJuCHdhncgcFVE/HfekKzZOOdZGS7GrFK//q9ZgBOTtb6IWA68PXcc1tyc86wMT+A3MzMzy8jFmJmZmVlGLsbMzMzMMnIxZmZmZpaRJ/BbpfY/+bTcIZiZVcY5z8pwMWaV2nrnXXOHYGZWGec8K8PFWIUmz16QO4TsViz6HwBGjH5P5kjMzHqfc56V4WLMKvXAf18JODGZWf/gnGdllJrAL2m8pIclLZM0rcH6QZKuTevnSxqZ2g+TtFDS/ennwTXb3JH2uSg9duypN2VmZmbWKrocGZM0ALgEOAxoBxZImhMRD9Z0mww8FRF7SJoInA98EHgSeH9ErJa0DzAPGFqz3QkR0dZD78XMzMys5ZQ5TDkOWJZu/YGka4AJQG0xNgH4Ylq+DrhYkiLi3po+i4EtJA2KiBc2OXKzFlR23uDMSWN7ORIzM2sWZQ5TDgVW1jxv55WjW6/oExEbgKeB7ev6fAC4t64Q+046RPkFpbvt1pM0RVKbpLZ169aVCNfMzMysdZQZGWtUJEV3+kjam+LQ5eE160+IiFWStgSuB04ErnjVTiJmADMAxowZU/+6TcNnSpZz4D+enTuEllDm/5NHz8yan3OelVFmZKwdGF7zfBiwurM+kgYCWwPr0/NhwA3ASRHxSMcGEbEq/XwWuIricKj1cYO334nB2++UOwwzs0o451kZZYqxBcAoSbtJ2hyYCMyp6zMHODktHwvcFhEhaRvgp8BpEfGLjs6SBkraIS1vBhwFPLBpb8VawfL5N7N8/s25wzAzq4RznpXR5WHKiNggaSrFmZADgFkRsVjSdKAtIuYAM4HvSlpGMSI2MW0+FdgD+IKkL6S2w4HngXmpEBsA3AJ8qwffV4/yIcie89Dt1wOw+36HZY7EzKz3OedZGaUu+hoRc4G5dW1n1iz/CTiuwXbnAOd0stt9y4dpZmZm1jeVuuirmZmZmfUOF2NmZmZmGbkYMzMzM8vINwq3Sh38qfNyh2BmVhnnPCvDxZhVaostt8kdgplZZZzzrAwfprRKLb3zJyy98ye5wzAzq4RznpXhYswq5cRkZv2egjRpAAAIHklEQVSJc56V4WLMzMzMLCPPGTNrUb6ZuJlZ3+CRMTMzM7OMXIyZmZmZZdTvD1P6JuDVOvxfLsodgplZZZzzrIx+X4xZtQYO2iJ3CGZmlXHOszJcjFmlltx2HQBvOfjYzJH0D57kb5aXc56V4TljVqlH77mFR++5JXcYZmaVcM6zMjwyZtbPefTMzCyvUiNjksZLeljSMknTGqwfJOnatH6+pJE1605L7Q9Lel/ZfZqZNTvnMTPrCV2OjEkaAFwCHAa0AwskzYmIB2u6TQaeiog9JE0Ezgc+KGkvYCKwN7ALcIukN6VtutrnJvFZkmY9pyc/T31llK1kbjQz61KZkbFxwLKIWB4RLwLXABPq+kwALk/L1wGHSFJqvyYiXoiIR4FlaX9l9mlm1sycx8ysR5SZMzYUWFnzvB3Yr7M+EbFB0tPA9qn97rpth6blrvYJgKQpwJT09DlJD5eIuSfsADxZ0Wt1pc/FMusj43ogFKAP/m56SNPGMusj3dp2154Opgd1mRsz5q+uNNP/j+5o2bhnfWTcK+Lu5ucgl5b9fdM8cZfKYWWKMTVoi5J9OmtvNCJXv8+iMWIGMOO1AuwNktoiYkzVr9uIY+lcM8XjWBprplh6WJe5MVf+6kqr/ps47mo57uqUOUzZDgyveT4MWN1ZH0kDga2B9a+xbZl9mpk1M+cxM+sRZYqxBcAoSbtJ2pxiQv6cuj5zgJPT8rHAbRERqX1iOttyN2AUcE/JfZqZNTPnMTPrEV0epkxzwKYC84ABwKyIWCxpOtAWEXOAmcB3JS2jGBGbmLZdLOn7wIPABuBTEfESQKN99vzb2yTNdGjBsXSumeJxLI01Uyw9prPcmDmsslr138RxV8txV0TFAJaZmZmZ5eDbIZmZmZll5GLMzMzMLCMXYw00yy1OJA2XdLukJZIWSzolVyw1MQ2QdK+kn2SOYxtJ10l6KP1+/k/GWP4l/fs8IOlqSVtU/PqzJK2V9EBN23aSbpa0NP3cNmMs/57+ne6TdIOkbaqIpb/pLF9IerukX0q6X9J/Sdqqk+2zfaZ6IPYsn0FJW0i6R9Kv0+ufndp3U3FrwKUqbhW4eSfbN7xdYDPHLekwSQvTv8lCSQe3Qtw1+xgh6TlJn60q7lIiwo+aB8VE3EeA3YHNgV8De2WKZWfgnWl5S+A3uWKpienTwFXATzLHcTnwsbS8ObBNpjiGAo8Cr0/Pvw9MqjiGA4F3Ag/UtH0FmJaWpwHnZ4zlcGBgWj6/qlj626OzfEFx1ud7U/tHgX/rZPtsn6lNiT3nZ5DiWnOD0/JmwHzgXSmGian9UuCTDbbdK/19GQTslv7uDGiBuN8B7JKW9wFWVfj/ZKPjrtnH9cAPgM9WFXeZh0fGXq1pbnESEWsi4ldp+VlgCX+5g0HlJA0D/hb4dq4YUhxbUfzRnwkQES9GxO8zhjQQeL2Ka+z9FRVfayoifk5xFnOt2luUXQ4ckyuWiLgpIjakp3dTXI/Lethr5Is9gZ+nbjcDH6jfNvdnalNiT7J8BqPwXHq6WXoEcDDFrQGh889fZ7cL7HWbEndE3BsRHb/fxcAWkgb1csgdr70pv28kHQMsp4i7qbgYe7VGtzjJVgB1kDSS4hvJ/IxhXAT8P+DljDFAMWq5DvhOOmT6bUlvyBFIRKwCvgqsANYAT0fETTliqbNTRKyB4g8dsGPmeDp8FLgxdxB9XV2+eAA4Oq06jldeqLZD03ymuht77s9gmrqxCFhLUTA+Avy+5gtIZ39Dsv6t2YS4a30AuDciXui9SF9pY+NO/58/B5xdVazd4WLs1crc/qlSkgZTDK2eGhHPZIrhKGBtRCzM8fp1BlIcCvtmRLwDeJ7iUFzl0lysCRSHGXYB3iDpwzliaXaSPk9xvcErc8fSlzXIFx8FPiVpIcUhwBcbbNYUn6mNiT33ZzAiXoqI0RQjvuOAtzTq1qAt69+aTYgbAEl7U0w7+HjvRNjYJsR9NnBhzchaU3Ex9mpNdYsTSZtRJKcrI+KHueIA9geOlvQYxaHbgyV9L1Ms7UB7RHSMEl5H8Yckh0OBRyNiXUT8Gfgh8O5MsdR6QtLOAOnn2pzBSDoZOAo4IdLEDet5jfJFRDwUEYdHxL7A1RQjCfWyf6Y2Ifam+Aymw7p3UMxh2iYdMoXO/4Y0xd+ajYi7Y8rKDcBJEdHo36TXbUTc+wFfSX/DTgVOV3HR5qbgYuzVmuYWJ5JEMYdjSURckCOGDhFxWkQMi4iRFL+T2yIiywhQRPwWWClpz9R0CMVdHnJYAbxL0l+lf69DKOa75FZ7i7KTgR/nCkTSeIrDA0dHxB9yxdHXdZYvJO2Yfr4OOINigvMr5P5MbUrsZPwMShqidHawpNdTFIZLgNspbg0InX/+OrtdYK/blLjTdj8FTouIX1QRb81rb3TcEfGeiBiZ/oZdBJwbERdXEngZOc4aaPYHcCTF2TyPAJ/PGMcBFMOt9wGL0uPIJvj9HET+sylHA23pd/MjYNuMsZwNPEQxv+W7wKCKX/9qirkyf6b4tj0Z2B64FViafm6XMZZlFHNjOv4PX5rz/05ffXSWL4BTUj77DXAef7nzyi7A3Jrts32meiD2LJ9B4G3AvSnuB4AzU/vuFIXVMooz9wal9qOB6TXbfz79nXkYOKLC3/dGx01RFD9f8++0CNix2eOu288XabKzKX07JDMzM7OMfJjSzMzMLCMXY2ZmZmYZuRgzMzMzy8jFmJmZmVlGLsbMzMzMMnIxZmZmZpaRizEzMzOzjP4/P6Z+BAZopHgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cond = [conditioned_propagate(guess=5., guess_sd=2., time=10.) for i in range(1000)]\n",
    "speed_cond, pos_cond = zip(*cond)\n",
    "plot_hist(speed_cond, pos_cond, \"Speed\", \"Position\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3) Inference\n",
    "The next step is inference- the difficult part. The aim is to find the posterior distribution of the speed, given our model and observation above.\n",
    "$$P({\\sf speed}|{\\sf position}=30) \\sim \\, ?$$\n",
    "Being able to calculate this exactly is rare, so we will estimate the distribution, we can use the stochastic variational inference algorithm. This requires a distribution called a guide (or variational distribution) that is an approximation the posterior distribution we want to know: $P({\\sf speed}|{\\sf position}=30)$. The guide relies on additional parameters that will be tuned to \"guide\" it towards the posterior distribution.\n",
    "\n",
    "For this example, we will sample the wind speed from a normal distribution with mean and variance given by the new parameters, a and b. These are first defined with the statement pyro.param(...) and are stored in the \"Pyro Parameter Store\". The values we provide here for a and b aren't too important, as we will optimise these in the next step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this example, we will sample the wind speed from a normal distribution with mean and variance given by the new parameters, a and b. These are first defined with the statement `pyro.param(...)` and are stored in the \"Pyro Parameter Store\". The values we provide here for a and b aren't too important, as we will optimise these in the next step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from torch.distributions import constraints\n",
    "def guide(guess, guess_sd, time):\n",
    "    \"\"\" The guide (or variational distribution) approximates \n",
    "    the posterior distribution, using extra variational parameters,\n",
    "    a and b, to sample from the speed distribution. \"\"\"\n",
    "    a = pyro.param(\"a\", torch.tensor(guess))\n",
    "    b = pyro.param(\"b\", torch.tensor(guess_sd), constraint=constraints.positive)\n",
    "    speed = pyro.sample(\"speed\", dist.Normal(a, b))\n",
    "    return speed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check the distribution of the guide before optimising a and b. In our simple example we can actually calculate the exact posterior distribution using rules of normal distributions, to find that posterior distribution of the cloud speed has a mean of 3.4 m/s with a standard deviation of 0.89 m/s, shown in red above. Our initial guess is far from the ground truth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact Mean = 3.40, Exact Standard Deviation = 0.89\n"
     ]
    }
   ],
   "source": [
    "speed_obs_mean = 30./10.\n",
    "speed_obs_sd = 10./10.\n",
    "guess = 5.\n",
    "guess_sd = 2.\n",
    "exact_mean = (speed_obs_sd**2 * guess + guess_sd**2 * speed_obs_mean)/(speed_obs_sd**2 + guess_sd**2) \n",
    "exact_sd = np.sqrt(guess_sd**2/((guess_sd**2 + speed_obs_sd**2)))\n",
    "\n",
    "print(\"Exact Mean = {:.2f}, Exact Standard Deviation = {:.2f}\".format(\n",
    "    exact_mean, exact_sd))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASYAAAEICAYAAADyYlmcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl4VOXd//H3l4QQwEAgYZM1KrLZB4SAVgTRogX3tWAf++hl+9BabC9t+3t+Wqu2Lq21Wun1PNqCS2n7UxH1sUWKS8FSl5a1Rssm+xKWJIQQdsjy/f1xZsIYJjNnZs5klvN9XddcJzNz5jt3BD7e933uc46oKsYYk07apLoBxhjTnAWTMSbtWDAZY9KOBZMxJu1YMBlj0o4FkzEm7VgwmYwjIotF5BupbodJHgsmExMRuVBE/i4itSKyT0Q+EpHRqW6XyS65qW6AyRwi0gmYD9wBzAXygHHA8VS2y2Qf6zGZWJwNoKovq2qDqh5V1XdV9VMRuS3Qe/rvQG9qnYh8KfhBEeksIs+LyG4R2Skij4hITsj7t4vIWhGpEZF3RKR/yHuXBurVisj/ANKqv7VpdRZMJhbrgQYR+Z2ITBaRLs3ePw/YDBQDDwL/KyJdA+/9DqgHzgLOBS4DvgEgItcCPwSuB7oBHwAvB94rBl4HfhSouwkYm6xf0KQHCybjmqoeAC4EFHgWqBKReSLSI7BLJTBDVetU9RXgM+CKwPuTgbtU9bCqVgJPAVMDn/sm8DNVXauq9cBPgRGBXtPlwBpVfU1V64AZwJ7W+Y1NqlgwmZgEwuM2Ve0DnAOcjhMWADv182eFbwu83x9oC+wWkf0ish+YCXQP7Ncf+FXIe/twhmu9A5/fEfL9GvrcZCcLJhM3VV0HzMYJKIDeIhI6/9MP2IUTJMeBYlUtDDw6qeqwwH47gG+GvFeoqu1V9e/AbqBvsGCgfl9MVrNgMq6JyGAR+b6I9Ak87wvcDCwJ7NId+K6ItBWRm4AhwAJV3Q28CzwpIp1EpI2InCkiFwU+9xvgXhEZFqjbOfB5gD8Dw0TkehHJBb4L9GyN39ekjgWTicVBnAnupSJyGCeQVgHfD7y/FBgI7AUeBW5U1erAe/+Bs7xgDVADvAb0AlDVN4CfA3NE5ECg5uTAe3uBm4DHgOpA/Y+S+lualBO7UJzxgojcBnxDVS9MdVtM5rMekzEm7VgwGWPSjg3ljDFpx3pMxpi0k7KTeIuLi3XAgAGp+noTi88+c7aDBnlQ6rNAqcRrmcyzcuXKvaraLdp+KQumAQMGsGLFilR9vYnFhAnOdvFiD0pNCJRKvJbJPCKyzc1+NpQzxqQdux6Tie5HP/KwlHe1TPayYDLRTZzoYSnvapnsZcFkoisrc7YjRnhQqixQKvFa6aquro7y8nKOHTuW6qakTH5+Pn369KFt27Zxfd6CyUR3113O1oMJ67sCtbJ58ru8vJyCggIGDBjA5y+24A+qSnV1NeXl5ZSUlMRVwya/jfHYsWPHKCoq8mUoAYgIRUVFCfUYLZiMSQK/hlJQor+/BZNxp6EBPv4Y7BQm0wosmExkK1bA2rXw0UcwciS88kqqW2RcyMnJYcSIEU2Pxx57zLPaZWVlLFiwwLN64djkt2mZKlxzDVRWnnztzTdh6tSWPxPFT3/6Uw8aZqJp37590xFQr5WVlbFixQouv/zypNQH6zGZSLZtg127oGtXmD/fee1vf0toOHfBBRdwwQUXeNRAE4va2loGDRrUdL7izTffzLPPPgvAHXfcQWlpKcOGDePBBx9s+szy5cu54IILGD58OGPGjKG2tpYHHniAV155hREjRvBKsnrQqpqSx6hRo9SkublzVUH1/PNVP/hAtWtX5/nmzXGX/Oijj/Sjjz7ysJHpZ82aNSefODHu/SOKNm3a6PDhw5sec+bMUVXVd999V88//3x9+eWX9ctf/nLT/tXV1aqqWl9frxdddJF+8sknevz4cS0pKdFly5apqmptba3W1dXpb3/7W50+fXps/x2a/nOwQl3kgw3lTMuWL3e2u3c7p6WMGwd/+pPTa4pzfcoPf/hDILvXMaWDloZyl156Ka+++irTp0/nk08+aXp97ty5zJo1i/r6enbv3s2aNWsQEXr16sXo0aMB6NSpU6u134ZypmXLljnbggJne1Hgpibvv5+a9mSiZPWZ4tTY2MjatWtp3749+/btA2DLli088cQTLFq0iE8//ZQrrriCY8eOoaopW/ZgwWTCa2iAlSudn4PBNH68s7VgylhPPfUUQ4YM4eWXX+b222+nrq6OAwcO0LFjRzp37kxFRQVvvfUWAIMHD2bXrl0sD/ScDx48SH19PQUFBRw8eDCp7bRgMuGtWweHDkG/fpCX57w2fLgTUps2wc6dqW2fiejo0aOfWy5wzz33sH79ep577jmefPJJxo0bx/jx43nkkUcYPnw45557LsOGDeP2229n7NixAOTl5fHKK6/wne98h+HDh3PppZdy7NgxLr74YtasWZPUyW+bYzLhBeeXxoyBqirn59xcGDsW3n4bPvggoWUDJrkaGhrCvr527dqmn3/5y182/Tx79uyw+48ePZolS5ac8nqwF5UsFkwmvOBfvNGj4bLLTr4+frwTTO+/H1cwzZgxw6MGmmxmwWTCC058jxnz+cudBNcgBeefYpTNlzsx3rE5JnOq48fhk09ABEaNgoULnQfAWWc5261b4yq9cOFCFgZrGdMC6zGZU61aBXV1MGSIM9n9yCPO6xMnQq9e0Latc5rKkSPQoUNMpR8J1LIrWZpIrMdkTrV5s7MdPPjU99q0gf79nZ+3ubrhhTExs2Ayp9q+3dn26xf+/eD9AOMczhkTjQ3lzKksmDz19dneHlp//rbRrvarqKjg7rvvZsmSJXTp0oW8vDz+67/+i+uuu67Fz1x++eW89NJLFBYWfu71H//4x5x22mn84Ac/SKjtblmPyZwqOESzYMpYqsq1117L+PHj2bx5MytXrmTOnDmUl5dH/NyCBQtOCaVUsGAyp2reY5o503kEJRBMM2fOZGZoLZMU7733Hnl5eXzrW99qeq1///585zvfYfbs2dx5551Nr1955ZVNJ1UPGDCAvXv3AvDoo48yaNAgJk6c2HSpFIBNmzYxadIkRo0axbhx41i3bp3n7behnDlV82AaNOjz7wevLBBHMA1qXsskxerVqxk5cmTcnw/2sD7++GPq6+sZOXIko0aNAmDatGn85je/YeDAgSxdupRvf/vbvPfee141HXAZTCIyCfgVkAM8p6phr9MpIjcCrwKjVXWFZ600refwYaiuds6P697dee3NN53tVVc52wR6TG8Gal0VrGVaxfTp0/nwww/Jy8tj+vTpUff/4IMPuO666+gQWA5y9dVXA3Do0CH+/ve/c9NNNzXte/z4cc/bGzWYRCQHeBq4FCgHlovIPFVd02y/AuC7wFLPW2laz44dzrZvX2dpAMCTTzrbYJj07OkEV2WlE2QdO7ou/2SglgVTcg0bNozXX3+96fnTTz/N3r17KS0tJTc3l8bGxqb3WrrNUrhLnjQ2NlJYWJi0y/YGuZljGgNsVNXNqnoCmANcE2a/h4HHAf/efjQbRDsiB7aWKQNccsklHDt2jF//+tdNrx05cgRw5pHKyspobGxkx44dLAuefhRi/PjxvPHGGxw9epSDBw829XQ7depESUkJr776KuBMsodecM4rboZyvYEdIc/LgfNCdxCRc4G+qjpfRFo8nigi04BpAP0i/cU3qeMmmMAZzm3Y4Aznhg5NdqsymtvD+14SEf74xz9y99138/jjj9OtWzc6duzIz3/+c8aOHUtJSQlf+MIXOOecc8LORY0cOZIpU6YwYsQI+vfvz7hx45ree/HFF7njjjt45JFHqKurY+rUqQwfPtzT9rsJpnCXsGu6hJ6ItAGeAm6LVkhVZwGzAEpLS+0GZekolmACWzKQxnr16sWcOXPCvvfiiy+GfX1ryJ/nfffdx3333XfKPiUlJbz99tuetLElboZy5UDfkOd9gF0hzwuAc4DFIrIVOB+YJyKlXjXStCILJpMG3PSYlgMDRaQE2AlMBb4afFNVa4Hi4HMRWQz8wI7KZahwwfSHP5y6X5zB9IdwtYxpJmowqWq9iNwJvIOzXOAFVV0tIg/h3IplXrIbaVpRuGDq2/fU/eIMpr7hamWhVF7IPx1ogreSd7WOSVUXAAuavfZAC/tOSKhFJnUaG08uFwgNpuB1nadMOflanMEUvEb0lNBaWSY/P5/q6mqKiop8GU6qSnV1Nfn5+XHXsJXf5qSKCjhxAoqLP3+dpeAh59AwCa5lqqqKaS1T8PB1NgdTnz59KC8vpyp4rXQfys/Pp0+fPnF/3oLJnOR24huctUy9e8OWLc4NMYNXtjS0bduWkjhvCGocdhKvOSmWYAKn1wSwZ09y2mN8y4LJnGTBZNKEBZM5KfQ8OTcsmEyS2ByTOSkYML16ff71114Lv38cwfRaS7WMCWHBZE6qqHC2PXp8/vXi4lP3Dd0vhmAqbqmWMSFsKGdOCgZM82CaPdt5NBdHj2n27Nkt3o7amCALJnNSSz2maMEU/JwLFkzGDQsm4zhxAmpqnPVJRUXuPmOT3yZJLJiMo7LS2XbvDjk57j4T7FlVVDinsxjjEQsm42hpGBdJfj4UFjq3E6+pSU67jC9ZMBlHPMEENpwzSWHLBYyjpSNyAAsWnPpaUM+esG6d8/lhw6J+zYJItYwJsGAyjkg9ptArDTQX41qmDpFqGRNgQznjCAZTcGgW6plnnEc4MQ7lnnnmGZ5pqZYxARZMxhGpxzR3rvMIJ8a1THPnzmVuS7WMCbBgMg6b/DZpxILJOCyYTBqxYDIOCyaTRiyYjLNAsrraOR0l1rP/LZhMEthyAXPydJRu3cKfjrJ4ccufLS4GEdi71wm4tm0jftXiSLWMCbAek4l/GAeQm+sEmqpzxxRjPGDBZKIH0xNPOI+WxDCce+KJJ3giUi1jsGAyED2Y5s93Hi2JYS3T/PnzmR+pljFYMBlIbCgHJ4Np925v2mN8z4LJnByChTsdxY3u3Z1tcBLdmARZMJnEe0zdujlbm/w2HrHlAiZ6MLVvH/nzMQRT+2i1jMGCyUD0YHrrrcifDw7lXATTW9FqGYMN5QycnBuyoZxJExZMftfQ4KzahpbvjvLww86jJTEE08MPP8zDkWoZgwWT2bfPWbXdpUvLp5MsWuQ8WhIMpspKp1YEixYtYlGkWsZgwWSCvZxguMSjY0fnjinHjsHhw960y/iaBZPfeRFMIjFNgBsTjQWT33kRTKGft2AyHrDlAn7nJpjc3DLcZTAVub39uPE1V8EkIpOAXwE5wHOq+liz978FTAcagEPANFVd43FbTTK4CabXX49eJ3QCPILX3dQyvhd1KCciOcDTwGRgKHCziAxttttLqvoFVR0BPA780vOWmuSwoZxJQ27mmMYAG1V1s6qeAOYA14TuoKoHQp52BCIfMzbpw00w3Xuv84jEZTDde++93ButlvE9N0O53sCOkOflwHnNdxKR6cD3gDzgEk9aZ5LPTTD94x/R67g8KvcPN7WM77npMUmY107pEanq06p6JvB/gR+FLSQyTURWiMiKKuvypwcbypk05CaYyoG+Ic/7ALsi7D8HuDbcG6o6S1VLVbW0W6L/EIw3LJhMGnITTMuBgSJSIiJ5wFRgXugOIjIw5OkVwAbvmmiSRvXkeXJeBZNdLM54IOock6rWi8idwDs4ywVeUNXVIvIQsEJV5wF3ishEoA6oAW5NZqONR/bvh/p6KCiAdu1a3q9Pn+i1XPaY+ripZXxPNMpJl8lSWlqqK1asSMl3m4D162HQIDjjDNi0KbFaqs4F5Y4fd86X69DBmzaarCIiK1W1NNp+dkqKn4Xe6DJRIjbPZDxjweRnbie+77rLeUTjIpjuuusu7nJTy/ianSvnZ26DqazMXT0XwVTmtpbxNesx+VkwQIKLIxNlR+aMRyyY/MyrNUxBNsdkPGLB5GdeB5NdLM54xOaY/MxtMJ19trt6LnpMZ7utZXzNgsnP3AbTrFnu6rkIplluaxlfs6Gcn9kck0lTFkx+peo+mKZNcx7RuDgqN23aNKa5qWV8zYZyfnXwIJw44Zw6Eu30kfXr3dV00WNa77aW8TXrMfmV18M4gM6dnZtmHjrk3GPOmDhZMPlVMoLJzpczHrFg8qtkBFNoPQsmkwCbY/Kr4AS1m9NRRoxwXzdKMI2IpZbxLQsmv4qlxzRjhvu6UY7MzYillvEtG8r5lZfXYgplQznjAQsmv4rlygK33OI83Ihyvtwtt9zCLW5rGd+yoZxfxTKUKy93XzdKj6k8llrGt6zH5FexTH7HwoZyxgMWTH5lywVMGrNg8qNYzpOLlV3F0njA5pj86OBB5zZLHTpAx47R9//iF93XjtJj+mIstYxv2X3l/GjTJjjrLOjfH7Zu9bZ2Y6Nz88z6eud8uUg30jS+Y/eVMy3z+iYEodq0geJi5+fg7ceNiZEFkx/Furjyhhuch1sRhnM33HADN8RSy/iSzTH5UawT39XVsdWPEEzVsdYyvmTB5EfJHMqBqyNzX5+9/JTXnr9tdHLaYzKODeX8KFnnyQXZWiaTIAsmP0rWGqYgCyaTIBvK+VGsp6N86Uux1Y9wIu+XArW2xlbR+IwFkx/F2mO6//7Y6kfoMd0fqBVujsmYIBvK+VFrTX7bUM7EyYLJb1Rjn/yePNl5uBXhqNzkyZOZHEst40s2lPObWO4nF3T0aGzfEaHHdDTWWsaXrMfkN8kexgF07eqcmrJ/P9TVJe97TNayYPKbZK9hgs+fL2fzTCYOFkx+0xo9JoAePZxtRUVyv8dkJVdzTCIyCfgVkAM8p6qPNXv/e8A3gHqgCrhdVbd53FbjhXgWV155Zezf07Mn/OtfpwTTlYFaa2OvaHwkajCJSA7wNHApUA4sF5F5qromZLePgVJVPSIidwCPA1OS0WCToHiGcj/4QezfE+wx7dnTrJRTy9YxmUjcDOXGABtVdbOqngDmANeE7qCqf1XVI4GnS4A+3jbTeCbZp6ME2VDOJMBNMPUGdoQ8Lw+81pKvA2+Fe0NEponIChFZUWWToqkR7MH07On+MxMmOI9YtBBMEyZMYEKstYzvuAkmCfNa2OvxisgtQCnwi3Dvq+osVS1V1dJuyf4/tgkvnmCKR7C+9ZhMHNxMfpcDfUOe9wF2Nd9JRCYC9wEXqepxb5pnPBcMiiQFU3DuaNiq/XwPWPvxeoYk5ZtMNnPTY1oODBSREhHJA6YC80J3EJFzgZnA1apq9+1JZ63UY6rt1BWATgf2JfV7THaKGkyqWg/cCbyDc5R3rqquFpGHROTqwG6/AE4DXhWRMhGZ10I5k0rHj8O+fZCTA0VFSf2qA50DwVRrwWRi52odk6ouABY0e+2BkJ8netwukwyh12FqE8Pa2q98JeavOljQhUZpQ8fDtc5pKW3bBko5tVbGXNH4iZ3E6yfxzi99+9sxf5W2yeFgQSGdD+xzliicfnqglFPL1jGZSCyY/CTe+aUjgSVqbq9GEHCgU1c6H9jHT559j+39BwFQf/wYALnt8mNrg/EVCyY/CQZTcI2RW5df7mwXL47pYwfCTIC/+9RdTsl7fhNbG4yv2Em8ftJaa5gCDnR2Jtg71dq95ExsLJj8pLWDyZYMmDhZMPlJkhdXNlcbWDLQ2XpMJkYWTH4S7xxTnKzHZOJlk99+Eu9Q7rbbwr4c7ZB/0xxTSDANvDCOazsZ37Fg8hOPgymaptNSai2YTGwsmPzi8GE4dAjatYPOnWP77N69zjZ4HW+XDoSZYzp2cD8A+QWFp+wfrgf2/G2jY/pOkx0smPwiOPHdowdIuCvZRHDjjc42xnVMh04rpFGEjodradNQT2NOLu89fQ9g65hMZBZMfuFyGBe21xLnVzbm5HKwoAudD+yj4OB+agtj63EZ/7Kjcn7RymuYguzInImHBZNftPIapqBgMNlaJhMLG8r5RSuvYQpKxmkpLS1TsIny7GHB5BeJDOXuuCPur91f6ART4X7nyN7gi2+Iu5bxDwsmv0gkmKbEf4vAmi5OD61LjXORujPOuzTuWsY/bI7JLxIJph07nEcc9ndx7oYTDKZD1RUcqrY7p5jIrMfkF4nMMX3ta842xnVMAPu6dAdOBtP7zz4I2DomE5n1mPygoQF2Be641TvSvUq9V9M1EEz77OY5xj0LJj+oqID6eueUkvzWvaTtgU5daZQ2FBysIae+rlW/22QuCyY/KC93tn37Rt4vCRpzctlfWEwb1aYjc8ZEY8HkB8GJ6z59UvL1zSfAjYnGJr/9INEe0/e/n9DX13TpDqymsKaScyb9e0yftds8+ZMFkx8Ee0zxBtNVVyX09TWBI3Ndayrp9+XYgsn4kw3l/CDRodxnnzmPONUEh3L7KqndvY3a3dvirmX8wXpMfpDoUO6b33S2caxjgpM9psKaKj763c8AW8dkIrMekx8kOpRLUE2zRZbGRGM9pmznweLKdXsOAvCLOCeig4ssu9ZUQvHpcdUw/mI9pmy3Z48TTt27O9f7ToHgHFPn/VWIakraYDKL9ZiyXXB+KczEd2sdiq9v246Dp3Wm4FAtxQ31VOW2bZXvNZnLginbeTC/NP+q2xNuRk2X7hQcquWiL05iQ8/+CdcLx+6ykj0smLKdB8G0dtiYhJuxv0t3+u3YwNldunPEg3omu9kcU7aLMJRzq+/29fTdvj6hZgQvf5K3aRXVCdYy2c+CKdt50GOa+tIvmfrSLxNqRvDI3Iklb7M0wVom+1kwZTsPekxe2F/oHJnraZc+MS5YMGW7FC+uDAr2mHrXnUhpO0xmcBVMIjJJRD4TkY0ick+Y98eLyD9FpF5EbvS+mSYu9fUpu3Jlc5XdnR5b/7rjKW2HyQxRg0lEcoCngcnAUOBmERnabLftwG3AS1430CRgzx5obHSu852Xl9Km7Ovak4Y2OfSqr6NdY2NK22LSn5vlAmOAjaq6GUBE5gDXAGuCO6jq1sB79jcunWzf7mwTHMb97w3fTrgpDbm5VBf1pHvVTiZdfD0202QicTOU6w2E3runPPBazERkmoisEJEVVVVV8ZQwsdiwwdmedVZCZTYN/Dc2Dfy3hJsTHM4Nbd8x4Vomu7kJJgnzWlwnPKnqLFUtVdXSbt26xVPCxCIYTAMHJlTmzA2fcuaGTxNuTlUgmPJXL0u4lsluboZy5UDoWKAPsCs5zTGe8iiYrn/9GQB+keA1lII9ptwV78G/J3a5XrfsNJXM5KbHtBwYKCIlIpIHTAXmJbdZxhMeBZNXgsE0wI7MmSiiBpOq1gN3Au8Aa4G5qrpaRB4SkasBRGS0iJQDNwEzRWR1MhttXFBNu2AKDuX6n7BgMpG5OolXVRcAC5q99kDIz8txhngmXVRUwKFD0KULFBWlujUAVHVzjpn0qztOm4Z6GnPsHHITnq38zlZp1lsCONEun925bWkLdN1XkermmDRm/8vKVh4G05yvfi/hGkH7ep9Br22f0a1yJ3u7pXY1uklf1mPKVh4G045+Z7Oj39kJ1wGo6eu0p3vljih7Gj+zYMpWHgbTkNXLGOLR2qONgW33yp2e1DPZyYZy2crDYLryzRcAb65kuWTTKm4DulWWJ1zLZC/rMWUjVdgY6Juk0eQ3wNa2zp1aulswmQisx5SNdu2CI0eguBgKC4HWuyNKNFvznGDqVrXTCVAJd8aT8TvrMWWjNFwqEFSbk8u+nBzyjx+1JQOmRRZM2SiNgwlgdbsOAPTdsSHFLTHpyoIpG61b52w9Cqbf33ovv7/1Xk9qjb31XmpGXgSQ8J1XTPayYMpG//ynsx0xwpNyFb36U9HLm5tUdu7Vn4ohztn9/bZ95klNk30smLJNY+PJYBo1ypOSw8s+YHjZB57U2l72AcuOHASg3w7rMZnw7Khcttm0CQ4cgF69nIcHLnv7RQA+GTEu4Vqr3n6RdarU5ebRrWoX7Y8c5GiHgoTrxsKu0ZT+LJiyzcqVAJR1P5P/TpMlAs3Vi7CzzxkM2LqOvjs2sH7QyFQ3yaQZG8plm0AwbRswOMUNiWx7v0EA9NtmwzlzKgumbNMUTENS3JDItgdOCu5r80wmDAumbKLaNPG9rX9695iCVyvoZ0sGTBg2x5RNNm2C2lr2Fxazv4t3d6F57j9/4lmt8YFa5R1Oo1GE03duJqe+jobctp59h8l81mPKJoFh3FaPe0s1RT2oKerhSa3TinpwWlEPjrXvSGX3PuQ21HP6ri2e1DbZw4IpmyRp4nv00r8weulfPKm1eelf2ByoFRzOlWxZE+kjxocsmLJJkia+J/z1dSb89XVPaq376+usC9Raf/a5AAxZk57LGkzq2BxTtjh2DJYsAWBLSXofkQtadc75AAxdvRRpbEDb5KSsLbboMr1YjylbLF7sXIPp3HM50Lk41a1xpbJHX6qKe3Ha4QP032rnzZmTLJiyxZ//7GyvuCK17YiFCKsDvaZzVv0jxY0x6cSCKRuoZmYwAauHOcE0bPXSFLfEpBObY8oG69bBli3OpXRHj4Z1//S0/K+nP+ZZrUua1Vo7tJSGNjmcself5B89xLH2p3n2XSZzWY8pG8yf72wnT4Yc7yeQDxUUcqig0JNa+QWF5IfUOtqhgM1nDCO3oYHBa1d68h0m81kwZYPgMO7KK5NSfuyH8xn74XxPam34cD4bmtVaM+w8wOaZzEk2lMt0NTXw4YdOT+myy5LyFRcEguSjCxMPvmAoDQypVTZiHNf86VnGLFvI3Cl3caJdfsLf4wVbQpA6FkwZJOw/lL1/g4YGmDix6VZNmWb7gMFsPmMYZ2xezflL3ub9i65NdZNMitlQLoPl1NfBjBnOk7vvTm1jErRw4hQALlk41znKaHzNekwZ7Lwl78DOnZT3PpMHK4ohTa9Y6caK0V9iypwZ9C3fyNmffcz6wel5VUsb3rUO6zFlKGlsZNLb/w+AtyffkvF3tG3IbcvfJlwPwJcWzU1xa0yqWY8pQw0v+4DeOzezr0t3lp2XnEnvoF/dPcOzWpdFqLX44uu5/M+/ZeTKxZy14RM2Dhzu2feazGKR6/ciAAAFc0lEQVQ9pgzUqXYv//G7nwHw7pe/mvSLrJ1ol+/ZkbLcdvnktlCrtrCYd758C220kW88+yD5Rw958p0m81iPKcO0aajnm7+5n84H9rFu8CgWTfxK0r/z4vdeA+Cvl9yYcK21gVpDWqj1p+umcc6qJfTf/hk3v/QUv/36/Ql/Z7KFm3cCm3tKhPWYMog0NjBlzgwGr1vJ/s5FzPzWwzTmJP//LaXLFlK6bKEntbYsW8iWCLUactvy7LSfcKJtOy788E2unPe8HaXzIesxZYpdu/j+L+5kyLqVNLTJYea3Hs2Yy5vEanfvM3jxlv/DrbMf5bo3ZtKnfCO/vf1+jud3SHXTYmJH8OLnKphEZBLwKyAHeE5VH2v2fjvg98AooBqYoqpbvW2qf4T+he65eyvj//ZHLvxwPkMOH6C2U1eenfZQ2h5O98qH46/mQKcuTJv5AKOXL2Lw2pUsvvh63r/oWvYV9Ux180ySRQ0mEckBngYuBcqB5SIyT1VDL9T8daBGVc8SkanAz4EpnrSwqgq2b/ekVKsJN/RQPfVRXw91dXDihHORt8OHoaqKm/7yMT32bKf/tnV0ralsKrFq2Hk8/58PZm1PqblPR4zj0ftf4PbnH+KMzau56s0XuOrNF6gu6smWkqFUdTud/YXdOdyxgGP5HTiRl09DTi4Nubk0Sg6NbdqACIqzlEKbr6hIxRKLlVk4ezJkCHTwtjfrpsc0BtioqpsBRGQOcA0QGkzXAD8O/Pwa8D8iIqoeTA788Y8wbVrCZTLJpJCfj7Q/jeVjJvK3i65lW8nQlLUpVXafXsKjP3qBszZ8wqV/mcPQ1Uspqt5DUfWeVDctPt7dCSt9lJXBcG+Xdki07BCRG4FJqvqNwPOvAeep6p0h+6wK7FMeeL4psM/eZrWmAcGUGQTEcj3VYmBv1L0ymx9+R7DfM5vE+jv2V9WoNz1002MK199tnmZu9kFVZwGzXHznqY0QWaGqpfF8NlP44XcE+z2zSbJ+RzcD3nKgb8jzPsCulvYRkVygM7DPiwYaY/zHTTAtBwaKSImI5AFTgXnN9pkH3Br4+UbgPU/ml4wxvhR1KKeq9SJyJ/AOznKBF1R1tYg8BKxQ1XnA88AfRGQjTk9pahLaGtcQMMP44XcE+z2zSVJ+x6iT38YY09qycFGFMSbTWTAZY9JOxgSTiPxYRHaKSFngcXmq2+QlEZkkIp+JyEYRuSfV7UkWEdkqIv8K/BmuSHV7vCAiL4hIZWA9X/C1riLyFxHZENh2SWUbvdDC75mUf5cZE0wBT6nqiMBjQaob45WQ034mA0OBm0Ukm5d5Xxz4M8yWNT6z+fyCfYB7gEWqOhBYFHie6WZz6u8JSfh3mWnBlK2aTvtR1RNA8LQfkwFU9X1OXbd3DfC7wM+/AzL+1i8t/J5JkWnBdKeIfBroUmZ81zhEb2BHyPPywGvZSIF3RWRl4BSlbNVDVXcDBLbdU9yeZPL832VaBZOILBSRVWEe1wC/Bs4ERgC7gSdT2lhvuTqlJ0uMVdWROMPW6SIyPtUNMglJyr/LtLpQnKpOdLOfiDwLeHPP6vTg5rSfrKCquwLbShF5A2cY+35qW5UUFSLSS1V3i0gvoDLqJzKQqlYEf/by32Va9ZgiCfzhBl0HrGpp3wzk5rSfjCciHUWkIPgzcBnZ9ecYKvQ0rVuBP6WwLUmTrH+XadVjiuJxERmBM8TZCnwztc3xTkun/aS4WcnQA3hDnAu05QIvqerbqW1S4kTkZWACUCwi5cCDwGPAXBH5OrAduCl1LfRGC7/nhGT8u7RTUowxaSdjhnLGGP+wYDLGpB0LJmNM2rFgMsakHQsmY0zasWAyxqQdCyZjTNr5/0fvZ6M+r+kMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Distribution of guide before any optimisation\n",
    "speed_guide = [guide(guess=5., guess_sd=2., time=10.).item() for i in range(1000)]\n",
    "plot_hist(speed_guide, title1=\"Speed\", true_dist=(3.4, 0.89))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stochastic Variational Inference\n",
    "Next we want to optimise the guide function to get a the best approximation to the posterior distribution, according to some loss function- we'll use the Evidence Lower Bound Operator or \"ELBO\" (link to why). We use stochastic variational inference , `pyro.infer.SVI` :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 31.7 s, sys: 628 ms, total: 32.3 s\n",
      "Wall time: 34.2 s\n"
     ]
    }
   ],
   "source": [
    "%%time \n",
    "pyro.clear_param_store()\n",
    "import pyro.infer\n",
    "import pyro.optim \n",
    "adam_params = {\"lr\": 0.001, \"betas\": (0.90, 0.999)}\n",
    "optimizer = pyro.optim.Adam(adam_params)\n",
    "svi = pyro.infer.SVI(model=conditioned_propagate, \n",
    "                     guide=guide,\n",
    "                     optim=optimizer,\n",
    "                     loss=pyro.infer.Trace_ELBO())\n",
    "\n",
    "losses, a, b  = [], [], []\n",
    "num_steps = 20000\n",
    "for t in range(num_steps):\n",
    "    losses.append(svi.step(guess=5., guess_sd=2., time=10.))\n",
    "    a.append(pyro.param(\"a\").item())\n",
    "    b.append(pyro.param(\"b\").item())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAGoCAYAAABbtxOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xl8XHW9//HXJ3uTrmnS0j20lLXIYtlBNkU2wV3xegXFy8WrV73qdfde3PW6e7mKoP5ARQREFhWQtZStQAqF7nSh+5K0TZO02TOf3x/npJymk2SSTuZMMu/n4zGPzJzzPed8zjczZz7zPd/zPebuiIiIiIhIIC/uAEREREREsokSZBERERGRCCXIIiIiIiIRSpBFRERERCKUIIuIiIiIRChBFhERERGJUIIsIiIiIhKhBFlEJMuY2TozazazPZHH9WZ2lZk91cMy88ysJSxbb2bzzezYbmWONrP7wvmNZva4mZ2emb0SERk6lCCLiGSnt7n7yMjjEyks8wl3HwmMB+YBv++aYWazgKeBxcChwGTgbuAhMzst7dGLiAxhSpBFRIYZd+8A/gQcHZl8HfCsu3/F3Xe5e6O7/5wgif5+DGGKiGQtJcgiIsOMmRUB/wQsiEx+C3BnkuJ3AGeYWWkmYhMRGQoK4g5ARESSusfMOiKv/xNo72OZn5vZD4FSoBl4Z2ReBbA1yTJbCRpLxgFNAw9XRGT4UAuyiEh2eru7j408bkphmU+6+1igBLgU+LOZvSGctwOYlGSZSUACqEtL1CIiw4ASZBGRYcbdE+7+JLAauCCc/AjwniTF30vQN1mtxyIiIXWxEBEZWszMSqIT3L0lSaHTCC7SWxpO+jrwgpl9G/gRQXeNq4AP8XoSLSIiqAVZRCRb/bXbOMh3h9NPJ+hfvO9hZl2NHdd3lScYneKr7v4AgLuvAs4EjgPWEfQ9fhfwVnd/OmN7JSIyBJi7xx2DiIiIiEjWUAuyiIiIiEiEEmQRERERkQglyCIiIiIiEUqQRUREREQihtwwbxUVFV5VVRV3GCIiIiIyhCxcuHCHu1emUnbIJchVVVVUV1fHHYaIiIiIDCFmtj7VsupiISIiIiISoQRZRERERCRCCbKIiIiISIQSZBERERGRCCXIadLQ0s7C9bviDkNEREREDlJGEmQzKzGz583sZTNbamZfD6cfambPmdkqM7vdzIoyEc9g+Ogt1bzrl8/S3NYZdygiIiIichAy1YLcCpzn7scBxwMXmtmpwPeBn7j7bKAOuDpD8aTd0s31AHS6xxyJiIiIiByMjCTIHtgTviwMHw6cB/w5nH4L8PZMxCMiIiIi0pOM9UE2s3wzWwTUAA8Da4Dd7t4RFtkETMlUPCIiIiIiyWQsQXb3Tnc/HpgKnAwclaxYsmXN7Bozqzaz6tra2sEMU0RERERyXMZHsXD33cA84FRgrJl13e56KrClh2VudPe57j63sjKlW2iLiIiIiAxIpkaxqDSzseHzEcCbgeXA48C7w2JXAvdmIh4RERERkZ4U9F0kLSYBt5hZPkFSfoe7/83MlgF/MrNvAS8Bv8lQPCIiIiIiSWUkQXb3V4ATkkxfS9AfWUREREQkK+hOeiIiIiIiEUqQRUREREQilCCLiIiIiEQoQRYRERERiVCCnGartjfGHYKIiIiIHAQlyGmyt60TgPfc8GzMkYiIiIjIwVCCnGYdiaR3yxYRERGRIUIJsoiIiIhIhBJkEREREZEIJcgiIiIiIhFKkEVEREREIpQgi4iIiIhEKEEWEREREYlQgiwiIiIiEqEEWUREREQkQgmyiIiIiEiEEmQRERERkYiMJMhmNs3MHjez5Wa21Mw+FU6/zsw2m9mi8HFxJuIREREREelJQYa20wF81t1fNLNRwEIzezic9xN3/2GG4hARERER6VVGEmR33wpsDZ83mtlyYEomti0iIiIi0h8Z74NsZlXACcBz4aRPmNkrZvZbMxvXwzLXmFm1mVXX1tZmKFIRERERyUUZTZDNbCRwF/Bpd28AfgnMAo4naGH+UbLl3P1Gd5/r7nMrKyszFq+IiIiI5J6MJchmVkiQHN/q7n8BcPft7t7p7gngJuDkTMUjIiIiIpJMpkaxMOA3wHJ3/3Fk+qRIsXcASzIRj4iIiIhITzI1isUZwD8Di81sUTjty8AVZnY84MA64F8zFI+IiIiISFKZGsXiKcCSzLo/E9sXEREREUmV7qSXgvf96lk+fuuLcYchIiIiIhmgBDkFTW2dNLd3xh2GiIiIiGSAEuQUufu+vzfOX8OuvW0xRyQiIiIig0EJcgos0nv6xQ27+c79K/j8n1+OLyARERERGTRKkFPk4d/2zgQADS0d8QUjIiIiIoNGCXIKkg2/ISIiIiLDkxLkFLn3XUZEREREhj4lyKkwtSGLiIiI5AolyCIiIiIiEZm61fSQ1tLWyRMbd9PQ0q6uFiIiIiLDnFqQU7ByeyMA//fY6n3T1OlCREREZHhSgtwP2xta+PPCTXGHISIiIiKDSF0s+uGeRVviDkFEREREBplakAdIXZFFREREhiclyBnS0Zngnb94midX1cYdioiIiIj0Qglyhuza28aLG3bzmTtejjsUEREREelFRhJkM5tmZo+b2XIzW2pmnwqnl5vZw2a2Kvw7LhPxiIiIiIj0JFMtyB3AZ939KOBU4ONmdjTwReBRd58NPBq+HhL6O8yb+iyLiIiIDA0ZSZDdfau7vxg+bwSWA1OAy4FbwmK3AG/PRDxx0vjJIiIiItkt432QzawKOAF4Dpjo7lshSKKBCT0sc42ZVZtZdW3t0LnIrb6pnV8/uRZ31x34RERERIaIjCbIZjYSuAv4tLs3pLqcu9/o7nPdfW5lZeXgBdhPrR2dvc7/8t2L+dbfl/P8a7v2TTM1IYuIiIhktYwlyGZWSJAc3+rufwknbzezSeH8SUBNpuI5WEs213PEVx/kkWXbeyxT39wOQFtnAlcvZBEREZEhIVOjWBjwG2C5u/84Mus+4Mrw+ZXAvZmIJx32tgWtx0+82r8uH6ZeyCIiIiJZLVO3mj4D+GdgsZktCqd9GfgecIeZXQ1sAN6ToXgyTn2QRURERIaGjCTI7v4UPQ/gcH4mYsgW6oMsIiIikt10J71BpH7HIiIiIkNPprpY5JS1tXs470dP7HutfsciIiIiQ4dakA9Ssi4T81amduFefVM7J3zjIRaur0tzVCIiIiIyUP1OkM3sU2Y22gK/MbMXzeyCwQhuuOhIJJJ2tli4YRd1Te1c/9iqjMckIiIiIskNpAX5I+FNPi4AKoEPE4xGIT343gMr9j1XZwsRERGR7DaQBLkrx7sY+H/u/jLK+3q1YlsjzW0H3nVPQ7+JiIiIZJ+BJMgLzewhggT5H2Y2CkikN6zh5xePrwbAknRaTjZNREREROIxkFEsrgaOB9a6e5OZlRN0s5BetHQc2IIsIiIiItlnIC3IpwEr3X23mX0Q+CpQn96whrZkPSe6ulNs3t18wDQRERERyR4DSZB/CTSZ2XHA54H1wO/SGtUw1FsyrA4WIiIiItljIAlyh7s7cDnwM3f/GTAqvWENbckSXt1VT0RERGRoGEgf5EYz+xLwz8BZZpYPFKY3LBERERGReAykBfl9QCvBeMjbgCnAD9Ia1RDy4JJtKZVL1sVCbcoiIiIi2affCXKYFN8KjDGzS4EWd8/ZPsg1ja0HTEt6kd7ghyIiIiIiaTCQW02/F3geeA/wXuA5M3t3ugMbbpK1IOviPBEREZHsM5A+yF8BTnL3GgAzqwQeAf6czsCGkr2tHQNaritnXlO7J33BiIiIiMhBGUgf5Lyu5Di0s6/1mNlvzazGzJZEpl1nZpvNbFH4uHgAsWSFT972Ugqleu5ksW5nE/NW1vQ4X0REREQyZyAJ8oNm9g8zu8rMrgL+DtzfxzI3Axcmmf4Tdz8+fPS1jqz16Ir9k9tkXSeS9VWOWl2jVmQRERGRbNDvLhbu/p9m9i7gDIJc8EZ3v7uPZeabWdWAIhwm1iRJgF230hMRERHJOgPpg4y73wXclYbtf8LMPgRUA59197o0rDMr7W3rjDsEEREREUlByl0szKzRzBqSPBrNrGEA2/4lMAs4HtgK/KiXbV9jZtVmVl1bWzuATWXWQNqF1ZgsIiIikh1SbkF297TeTtrdt3c9N7ObgL/1UvZG4EaAuXPnDptU0vd7Pmx2S0RERGRIG8hFemlhZpMiL98BLOmprPRsbe2eAQ8zJyIiIiIHGlAf5P4ys9uAc4AKM9sE/DdwjpkdT9CQug7410zEMtyc96MnOKlqHHdee3rcoYiIiIgMCxlJkN39iiSTf5OJbcdhIHfI+/HDr/LuN06jvKyo38u+sG7YXtsoIiIiknGxdbEYzgbSm7ilPcFX7l4MQHtngpb2vke9uPnp1wawJRERERHpjRLkLNKVFF92/dMc+bUH+yz/88dWD3ZIIiIiIjlHCXKMehrabfnWgYyaJyIiIiLpoAQ5i5gZt7+wIe4wRERERHKaEuQs84W7FscdgoiIiEhOU4Icq/37WCQb/aK+qZ3n1u6kprGFdTv29lleRERERA5ORoZ5k4G76ubneWnD7n2v133vkn3Pde89ERERkfRTC3KWW7YlPRfsPbhkKy9t0HjJIiIiIn1RC3KMehrFostX71lMa0eix/n96WJx7R9eBPZvgRYRERGRA6kFeRDs2tualvX8YYFGtBARERHJNCXIWcR01Z2IiIhI7JQgD4K+uk6IiIiISPZSgpxF1u9sijsEERERkZynBHkQ/GLemgEtt6pmz4C32dLeOeBlRUREROR1SpBjULe3jXf/8hnuX7Itbeu85Zl1aVuXiIiISC7TMG8x+OydL1O9vg7WH9y4xDv3tu173pFQx2cRERGRdFALcgweW1Ez4GVrGlpYsrn+gOkJJcgiIiIiaZGxBNnMfmtmNWa2JDKt3MweNrNV4d9xmYpnqDr7B/O49H+f4tHl2/ebrvRYREREJD0y2YJ8M3Bht2lfBB5199nAo+Fr6UVzeDHe1bdUxxyJiIiIyPCUsQTZ3ecDu7pNvhy4JXx+C/D2TMUz3GjsZREREZH0iLsP8kR33woQ/p2QrJCZXWNm1WZWXVtbm9EAhwpXJwsRERGRtIg7QU6Ju9/o7nPdfW5lZWXc4WSlgbQgN7S08+CSrekPRkRERGQIi3uYt+1mNsndt5rZJGDgwzvkuN7y41c27d7v9cZdTTy5agePrdjOI8trmPe5c6iqKBvcAEVERESGiLgT5PuAK4HvhX/vjTec4aOxpZ1jr3uIb79jDl+5e8l+8664aQGb6pqZOm4E8PqFfyIiIiKS2WHebgOeBY4ws01mdjVBYvwWM1sFvCV8LQMR6WOxqa6JXz2xFuCA5BiCO/l1W6TbqpxfzFtNTWNL+uMUERERyXIZa0F29yt6mHV+pmIYzqK57uXXP73fXfai2joSmFmP6/naPUv4/YL1ADy1agd//JdT0xmmiIiISNYbEhfpSd+ircE9JccAr25v7HU9XckxQF1TO2tr9/Dgkm0HHZ+IiIjIUBF3H2RJk4MZ5q2nrhbLtzZw3o+eAGDd9y6hozOBA4X5+l0lIiIiw5cynWHi2TU7AVi/c2+v5Ras3bnv+ebdzfue1ze3s2Rzfa/LXvq/TzH7Kw8cRJQiIiIi2U8tyMPEixt28+yandzz0uZey33r78uTTv/gr59jcR8J8optvXfPEBERERkOlCAPI1fctGBAyzneZ3IsIiIikivUxUL6bdmWhrStq6ahhUt+/iTb6jWknIiIiGQHJcjCJT9/ql/lL/75k2nb9p9e2MjSLQ3c+tz6vguLiIiIZIASZBkUO/e0sr2h91bh5rZOfvzwqxmKSERERCQ16oMsKfGexoLrwRu/9QgQDA/Xk1/NX3NQMYmIiIgMBrUgS0p+kqSlt6MzwQ//sZKGlna+ds+SpGV609KeSFd4IiIiImmjFmRJyV0vHjh83N9e2cr1j6+mrqmNW5/bAMB/vOXwA8rdOH8N5x05gcMmjNpvevTmJj3f/Hp/jyzbTllxAafNGp968CIiIiL9oBZkSUn0piJd1tbuAdiXHCezt7WD79y/gjf/eD4L1+/iyVW1NLS08+2/L6OptbPfcXz0d9UDHs5OREREJBVqQZYB+/ljqw+Y1plw8vNebw++7PrXR8h41y+f7XFdDy3bzj2LtvDoZ8/edyvrj/1hIdXr63j+y+djlmobs4iIiMjBUQuyDMjVN7+QdPqCtTtpbGnf93pNbe+3vu6yYlsjG3Y1sWtv275pDyzZRm1jK+2dHpZJ3/jLIiIiIj1RC7IMyKMrapJO/6dfP5f2bTnO+p17ufCn6Rt/WURERKQnakGWrHfk1x5MepFgKtydG+evobaxlWfW7KC1o//9nkVERCS3ZEWCbGbrzGyxmS0ys+q445H4XPO7alra909i3eHnj646oOzjK2s49TuPsmJbAw8s3sodL2w8oMzyrY185/4VvPWn8/nATc/xrb8t58lVtexuauOelzZT9cW/U9MY3NCkua2TN/3P4zy7Zufg7JyIiIgMCdnUxeJcd98RdxASr5c31XPk1x7k8uMn91ruzuqN/OefXwHgO/evYP6rtQA0tnZw9ZmH7ivXHCbbXX2bH1m+nd8vWM9JVeOwcHC5v7+ylatOr+KJV2vYsKuJ7z6wnPs+cSYAiUQwGF30wkMREREZ3rIpQRbZ595FW3qd35UcA/uSY4Bv/m0Z7s6vn3yN6y47mmv/8OJ+y22tD1qLX1hXt2/a1/+6jBfW7eL+xdsAeGVT/b55h3/1AToS3usdAX/80EoAPnPBEX3tloiISL89vGw7dXvbeO9J0+IOJWdYf28hPChBmL0G1AEO/Mrdb+w2/xrgGoDp06e/cf369RmNr+qLf8/o9iR+5xxRSfW6Ova0dgBw1elV/Nu5s3hq1Q6mjB3BKTPH82+3LtyXVAM89+XzKcrPY/SIwv1anLfsbmb9ziZOPrR82LdEt3cm2NvawdjSopTK723toK0jwbiy5OUfWbadk6rKGVNamM4wRbLSw8u2c+Qho5hWXhp3KJJluvKQtd+5mAeWbOPCOYeQn2e4u4ZB7QczW+juc1MqmyUJ8mR332JmE4CHgX939/nJys6dO9erqzPbTVkJsgyGL110JLdXb+RdJ07l2rNn9Zg8N7V1MKIw/6APgu7O/Yu3MbOyjDEjCpk8dkTKy7a0d3LPS5vJzzPedtxkNu9uprggj/w844V1dTyxspbyskJuf2EjDS0dnHNEJd+8fA6F+Xmc+t1Hef9J0zh15ngOmzCS5vZOHlm2nV/NX3vAdj53weFcdtwUXtu5l2njRnDej57gDVPH8OdrTyfhzpFfe5BDRpfwzBfPIy9NPzZ2N7WxdEsDh00YycTRJX2Wr29uB4ey4nyWbGng+Glj0xIHQGNLO9Xr6qgcVcycKWMGtA53pzPhdLpTXJC/b/oP/rGCqvFlvGeuWqAGQ0dngudf28Xph1UAwY+/b/19Ge8/aTrjRxZx18LNvOOEKbR0dHL4xFGs27GXcaVF5OcbhfnGqu17uPR/g3Hjf3/1ybjDmw6vPOi43J3WjgQlhflJ57d3JjDALOh0lq7P1cFauqWePDOml5dSVjy8T3av37mX/7p3KZ+94HAK8vK4f/FW/v38w/b7/HbPQz5/4RHcWb2J13bs5cnPn8u08lIWrq+juCCP9s4E9y7aQmfC+dxbj6CxpZ0p4fH+P25fxLiyIj5yxqEUFeQxcXQJbR0JigrSe0laNLdMePZ0UxxyCXKUmV0H7HH3HyabrwRZRLLNe944lRFF+ayp3cOUsSO4o3pTWtZ71elVdCQS/GHBBo48ZBRtHQk6Ek7d3jYaWzuYMb6U9Tub0rKt4ey8Iycws6KMdTubeGT59rjDOSjvPGEKf3lpYKP6xOHE6WM5atJo4PW7rs6sLGNtimPkS+bNrChj7Y7B///c9bHTeOOM8kHfTlR/EuTYf5aZWRmQ5+6N4fMLgG/EHJaISMruXLiJooI82joSaV3vLc+uIz88c7BiW+MB85Ucp+axFTU8U5hHS3t6/z9xGErJMcCLG3azfmcTrZHPhpLj7JaJ5BigozO7Gmi7iz1BBiYCd4enjwuAP7r7g/GGJBKvS46dxN8Xb006b86U0XQmYMrYEfu1hp01u4InV+2gYmQxsyrLWLRxNw4U5edRUpjPhXMmcsK0ccyZMoZtDS0cP3UsRQV5dCQS+27vDex73plwzKCtI0F+ntHc1kl7IsHI4gK2N7RS39xOZyLBprpmSosKqGtqY29rB+5w2ISRlBUXUNPQQktHJ+4wojCf0SMKWbmtEQdGFuezqa6Z5VsbOOXQ8YwtLaRiZPG+0341jS00NHdwzOTRFBfkkQiPpSWFeYwsLiDhYAatHQk6E86yLQ2s2NbA+UdO5Nm1O3jbcZMpLsinub0Td2dkcQHu+59C7kx4cIrZgnrqOqHW1Zulq1tLZ8LpWszD7XZfV5S7s7upneVbG5g+vhR3KC3Kp7ysaN86u87euQcXXzS1dezbVvc+3ImE77etrmW71pVIOG2dCYoLgn3oSDgJd3bubWNzXTOjSgoYV1rExromxpcVUV5WxK69beTnGWXFBZSXFtHU3klNQwu3v7CRX81fS57B+UdN5F0nTqWqopRXt+/hlY27eeOMcazf1cT3HlhxwH739r7tcvqs8Tz/2i46Eql9Ob75qAk8sryG6eWlbNiV/AfBZcdN5r6Xgwt733/SNP70wkZ+9v7jqRxVzNGTRiftE+/ubG9o5dXtjYwfWcSsypGYQXFB/n7/m827m6lpbGXG+FJe2bSbrfUt/P2VrUwcXUJ5WREdnQnW72ri8uMnM3dGOeNHFrGtvoX2TqetI0FJYR61ja0kHA4ZU8LUcSP2ve8K8/MozA8+g4bR0t5JaXE+BXl55Nnr/9+2jgQ79rQyaUzQBWhrfQst7Z2MKy2ivrmdnXtbqRxZwiFjSijMj75Pgr9ravdw+wsb2VjXxKiSQj5wynQqRxYzsriA/Hyjpa2TvDwj34yRJQV0Jpw/LFjPjfPXUtPYCsC1Z8+itCif+xdv5YzDKigrLuDcIyo58pDRFBXkUd/czs3PrGPauBGs3NbIR848lPbOBNPGle5773bv7tH9DPauvW0s3lwfHh8KGDOikImjSnht517mrayhpDCfWZUjeXnjbs4/agKtHQlmVpbR0ek8/9ouNu9u5qUNdTyyPPmNrDLp91efzHNrd7G7uY0rT6vimTU7uejYQygpzGdUcQFmtu9z29TWSUdngtEjCmlobqe5vZMxIwpxD45v9c3ttHcm2LW3jZElBbSH123UN7dTObKY8pFFdHY6O/e2MnpEISMK89lW34IDVePLgvdbIkFTayc79rRSMbKYsaWFJML6Lwm7cjS1d1KUn0dhvtHSnuChZdtYU7OHkqJ85kwew+7mdi6ecwhbdrfQ2tHJH5/fQGtHgvfOnUbCnVkVIykpyqMwL49Od8743mNMGF1Mvhkvb6rnu+88lusfW01nwrn3E2ek1KUtTlnXxaIv6mKRm95y9ESOnzaWGeNLyTdj4pgSGls6mFlRxuSxIwa9f1NHZ9D6kWeWNX30RLJBV6LXk9U1jXzijy+xYlsjj332bGZWjkxa7ht/XcZvn36NVd++iHkra1myuZ5Pv3l2j33vG1vaeXDJNt79xqm6SEmSqm9up7axlZkVZXzuzy9z1elVXHb90/uVee7L51PX1MaM8jIc5+qbq1m4oY5X/vsC9rZ2UFSQR2lRAW/6n8fZvLuZr116NO8+ceoBFw63tHfu+x5q70xQWpQN7Y/S3ZDug9wXJcjZ6bAJI/nvtx1NR8I5rHIk63buZUZ5GeNHFlFadPAXmInI0OXuNLZ2MLpEo5FIvL52zxIeWraNz7/1SOZMGcMRh4yKOyTJoCHVB1ky71/PnsnZsytZU7uH0SMKKSsq4NwjJ+xrCUpHa6yGKRKRLmam5FiywjffPodvvn1O3GHIEKAEOUuVlxVx4vRxPLZiO8dOHcsVJ01jbGkRRx4yCifopzizouygTvd3DUfUJT8v+TBAIiIiIrlECXIKvvfOY/niXxb3OL+sKJ+9bcEtjWdPGMmokgLOP2oic2eM45Hl2znl0PGcObsiuFCpSEmoiIiISDZTH2QRERERGfb60wc5vbdOEREREREZ4pQgi4iIiIhEKEEWEREREYlQgiwiIiIiEqEEWUREREQkYsiNYmFmtcD6GDZdAeyIYbtDkeoqNaqn1KmuUqe6So3qKXWqq9SprlIXR13NcPfKVAoOuQQ5LmZWnerQILlOdZUa1VPqVFepU12lRvWUOtVV6lRXqcv2ulIXCxERERGRCCXIIiIiIiIRSpBTd2PcAQwhqqvUqJ5Sp7pKneoqNaqn1KmuUqe6Sl1W15X6IIuIiIiIRKgFWUREREQkQgmyiIiIiEiEEuQ+mNmFZrbSzFab2RfjjicOZjbNzB43s+VmttTMPhVOv87MNpvZovBxcWSZL4V1ttLM3hqZPuzr08zWmdnisE6qw2nlZvawma0K/44Lp5uZ/Tysj1fM7MTIeq4My68ysyvj2p/BYGZHRN43i8yswcw+rfdUwMx+a2Y1ZrYkMi1t7yEze2P4Hl0dLmuZ3cP06aGufmBmK8L6uNvMxobTq8ysOfL+uiGyTNI66aneh6Ie6iptnzkzO9TMngvr6nYzK8rc3qVPD/V0e6SO1pnZonB6rr+nesoPhv7xyt316OEB5ANrgJlAEfAycHTcccVQD5OAE8Pno4BXgaOB64DPJSl/dFhXxcChYR3m50p9AuuAim7T/gf4Yvj8i8D3w+cXAw8ABpwKPBdOLwfWhn/Hhc/Hxb1vg1Rf+cA2YIbeU/v2903AicCSwXgPAc8Dp4XLPABcFPc+p7muLgAKwuffj9RVVbRct/UkrZOe6n0oPnqoq7R95oA7gPeHz28APhb3PqernrrN/xHwX3pP9ZofDPnjlVqQe3cysNrd17p7G/An4PKYY8o4d9/q7i+GzxuB5cCUXha5HPiTu7e6+2vAaoK6zOX6vBy4JXx+C/D2yPTfeWABMNbMJgHW+Y2OAAAgAElEQVRvBR52913uXgc8DFyY6aAz5Hxgjbv3dofMnHpPuft8YFe3yWl5D4XzRrv7sx58+/wusq4hJ1lduftD7t4RvlwATO1tHX3USU/1PuT08L7qSb8+c2Gr3nnAn8Plh2xd9VZP4X6+F7itt3Xk0Huqp/xgyB+vlCD3bgqwMfJ6E70nhsOemVUBJwDPhZM+EZ4m+W3kNFFP9ZYr9enAQ2a20MyuCadNdPetEBxQgAnh9FyvK4D3s/+Xjd5TyaXrPTQlfN59+nD1EYJWpy6HmtlLZvaEmZ0VTuutTnqq9+EkHZ+58cDuyA+T4fq+OgvY7u6rItP0nuKA/GDIH6+UIPcuWT+XnB0Xz8xGAncBn3b3BuCXwCzgeGArwWkn6LnecqU+z3D3E4GLgI+b2Zt6KZvTdRX2UbwMuDOcpPdU//W3bnKmzszsK0AHcGs4aSsw3d1PAD4D/NHMRpNDdZJEuj5zuVKHV7D/D3q9p0iaH/RYNMm0rHxfKUHu3SZgWuT1VGBLTLHEyswKCd78t7r7XwDcfbu7d7p7AriJ4NQb9FxvOVGf7r4l/FsD3E1QL9vDU0Vdp95qwuI5XVcEPyJedPftoPdUH9L1HtrE/l0OhmWdhRf5XAr8U3hqlrC7wM7w+UKCvrSH03ud9FTvw0IaP3M7CE6XF3SbPmyE+/ZO4PauaXpPJc8PGAbHKyXIvXsBmB1emVtEcCr4vphjyriwz9VvgOXu/uPI9EmRYu8Auq74vQ94v5kVm9mhwGyCTvbDvj7NrMzMRnU9J7hYaAnBfnZdlXslcG/4/D7gQ+GVvacC9eHpqH8AF5jZuPCU5wXhtOFmv9YYvad6lZb3UDiv0cxODT/bH4qsa1gwswuBLwCXuXtTZHqlmeWHz2cSvI/W9lEnPdX7sJCuz1z4I+Rx4N3h8sOuroA3Ayvcfd8p/1x/T/WUHzAcjlcHe5XfcH8QXHH5KsGvwq/EHU9MdXAmwSmNV4BF4eNi4PfA4nD6fcCkyDJfCetsJZErTod7fRJc2f1y+FjatY8E/fMeBVaFf8vD6Qb8X1gfi4G5kXV9hODCmNXAh+Pet0Goq1JgJzAmMk3vqWCfbiM4ddtO0IJydTrfQ8BcgkRoDXA94V1Vh+Kjh7paTdCfset4dUNY9l3h5/Jl4EXgbX3VSU/1PhQfPdRV2j5z4fHv+bD+7wSK497ndNVTOP1m4NpuZXP9PdVTfjDkj1e61bSIiIiISIS6WIiIiIiIRChBFhERERGJUIIsIiIiIhKhBFlEREREJEIJsoiIiIhIhBJkEZEhysw+bWalccchIjLcaJg3EZEhyszWEYwjuiPuWEREhpOCvouIiEjcwjsz3kFwq9V8ghsxTAYeN7Md7n6umV0AfB0oJhhU/8PuvidMpG8Hzg1X9wF3X53pfRARGSrUxUJEZGi4ENji7se5+xzgp8AW4NwwOa4Avgq82d1PBKqBz0SWb3D3kwnuRPXTDMcuIjKkKEEWERkaFgNvNrPvm9lZ7l7fbf6pwNHA02a2CLgSmBGZf1vk72mDHq2IyBCmLhYiIkOAu79qZm8ELga+a2YPdStiwMPufkVPq+jhuYiIdKMWZBGRIcDMJgNN7v4H4IfAiUAjMCossgA4w8wOC8uXmtnhkVW8L/L32cxELSIyNKkFWURkaDgW+IGZJYB24GMEXSUeMLOtYT/kq4DbzKw4XOarwKvh82Ize46gYaSnVmYREUHDvImIDHsaDk5EpH/UxUJEREREJEItyCIiIiIiEWpBFhERERGJUIIsIiIiIhKhBFlEREREJEIJsoiIiIhIhBJkEREREZEIJcgiIiIiIhFKkEVEREREIpQgi4iIiIhEKEEWEREREYlQgiwiIiIiEqEEWUREJIuZ2TlmtimN66syMzezgnStsx/bnm5me8wsf4DLf9nMfp3uuES6U4IsaWFm68yszcwquk1fFB6IqzIczzlmlggPxI1mttLMPpzJGAYqrK/D4o5DRF5nZmea2TNmVm9mu8zsaTM7KZx3lZk9FXeM6RIez5vDY+fucL+vNbODzhncfYO7j3T3zhTiOOCHgbt/x90/erBxJNnWVWbWGX5nNITfXZemezvpFuePneFOCbKk02vAFV0vzOxYYER84bDF3UcCo4EvADeZ2dH9XclAWzrioIOkSPqZ2Wjgb8D/AuXAFODrQGuccaXiII4Jb3P3UcAM4HsEx9DfpC2w7PRs+J0xlmBf7zCz8v6swAJDJrfSd0bPhsw/UYaE3wMfiry+EvhdtICZFZvZD81sg5ltN7MbzGxEOG+cmf3NzGrNrC58PjWy7Dwz+2bYctNoZg91b7FOxgP3AHXA0eG67jSzbWFr0HwzOyaynZvN7Jdmdr+Z7QXONbNLzOylsGVho5ldFynf9Qv+w+G8urC15SQzeyVsgbm+Wz18xMyWh2X/YWYzwunzwyIvhy0Z7wunXxq2aHS15rwhsq51ZvYFM3sF2KsDnkjaHQ7g7re5e6e7N7v7Q+7+ipkdBdwAnBZ+ZncDpHjMuDI8Fu4ws69E5o8Ij0N1ZrYMOCkajJl90czWhMfBZWb2jsi8q8Jj5E/MbBdwnZnlh8fdHWa2Frgk1R1393p3vw94H3Clmc0Jt9PbsXx5tPXVzArCbZ/YvcUzPG4uD/dlrZn9azi9DHgAmBzW6x4zm2xm15nZHyLrvszMlobHxnnh/6Nr3joz+1x4HK43s9vNrCSFfU4AvyVo4JmZ4nfTt83saaApXCbpfoXlzzGzTWb2eTOrMbOtZvZ2M7vYzF614AzFlyPl8yL/851mFk3cu74zdod1dFq4TNLvmHCem9nHzWwVsKqv+shZ7q6HHgf9ANYBbwZWAkcB+cBGgtYHB6rCcj8F7iNohRkF/BX4bjhvPPAuoDScdydwT2Qb84A1BF9WI8LX3+shnnOATeHzPOAdQDtwRDjtI+E2isOYFkWWvRmoB84Ily0J13ds+PoNwHbg7WH5qnAfbwjLXgC0APcAEwham2qAs8PybwdWh/VUAHwVeCayfQcOi7w+MVz+lLBerwzruzhS94uAacCIuN8Leugx3B4EZ6F2ArcAFwHjus2/Cniq27RUjhk3hcey4whao48K538PeDI8Tk4DlnQdz8L57wEmh+t+H7AXmBSJpQP49/D4MgK4FlgRrqsceDzcfkEP+7sOeHOS6RuAj4XPezuW/xdwa2S5S4AV3fa9IDJvFmDA2QQJ5omROtzULYbrgD+Ezw8P9/0tQCHw+fDYWhTZj+fDuioHlgPX9rDP+/6HYb19CmgExpDad9MG4Jhw2cIU9qsjrKdC4F+AWuCP4fqPIfgOmRmW/zSwAJhK8J31K+C2ZPUZTkvlO+bhsE70ndHT5z7uAPQYHg9eT5C/CnwXuDD8ABaEH8aq8ECxF5gVWe404LUe1nk8UBd5PQ/4auT1vwEP9rDsOUAC2A3sIkgg399D2bFhjGPC1zcDv+tjf38K/CR83nWAmhKZvxN4X+T1XcCnw+cPAFdH5uWFB88Z4evuCfIvgW922/5KXk+41wEfifs9oIcew/kRJhs3A5vC5OY+YGI47yq6JchJlk92zJgamf981zEKWAtcGJl3Dd0SxW7rXgRcHollQ7f5jxFJDAl+xA8kQV4AfKWvYzlwGEFyWRq+vhX4r2773tO27wE+FT4/p/t+s3+C/DXgjsi8PGAzcE5kPz4Ymf8/wA09bPeq8P+6G9gR7usBdRCWTfbd9I0+/v/d96sZyA9fjwrr5JRI+YW8/oNqOXB+ZN4kggafgmT1SWrfMefF/ZnK9odOxUq6/Z7glM+hdOteAVQS/AJfaGZd04ygVRQzKwV+QpBcjwvnjzKzfH/9go5tkfU1ASN7iWWLu0/tPtGCPsXfJmiFqSRIpAEqCFqOIWj9ji5zCkGrzhygiOBX/J3dVr098rw5yeuuWGcAPzOzH0U3QdDSvD7JfswgOLX575FpRQStIl02IiKDxt2XEyRRmNmRwB8Ikt4rkpVP8ZjR0/FsMvt/pvc7LpjZh4DPECRHhMtFu5t1Px70ur5+mELQ4NDrsdzdV5vZcuBtZvZX4DLghGQrNLOLgP8maA3OC9e7OMV4JhPZF3dPmNnGMM4u3es4etzsboG7n5kkxlS+m7p/Z/S1XzsjyzaHf3v7zrjbzBKR+Z3AxB72I5XvGH1n9EF9kCWt3H09wcV6FwN/6TZ7B8GH/hh3Hxs+xnhwUQTAZ4EjCH5FjwbeFE430usDwOUELd5jeP1LJrod77bMHwlajKa5+xiC7hQDjWsj8K+ROhjr7iPc/Zleyn+7W/lSd7+tl3hFZJC4+wqC1uQ5XZOSFDuYY8ZWgu4QXaZ3PQn7kt4EfAIY7+5jCbpg9Hb86nF9qbJgxI4pwFP0fSwHuI3gx8PlwDJ3X51kncUEZ9d+SNAaPxa4P7IvfR3XthAkg13rM4L93Nzf/etDKt9N+2JNYb/6ayNwUbfvgBJ330zyOkrlO0bfGX1QgiyD4WqC0zd7oxM9uPDhJuAnZjYBwMymmNlbwyKjCA66u8MLEP57kOIbRdDfbyfBr/rvpLjMLndvMbOTCZLsgboB+JKFFwaa2Rgze09k/nZgZuT1TcC1ZnaKBcosuABo1EHEICIpMrMjzeyzXRdmmdk0guRvQVhkOzDVzIoiix3MMeMOgmPEuHCb0bNHZQTJTW0Yy4d5PVHvbX2fNLOpZjYO+GKqgZjZ6PCCuz8RdG1YnMKxnLD8BcDHCH4sJNPVsl4LdIStrhdE5m8HxpvZmF726xIzO9/MCgkS2Vagp8aGgervd1Nf+9VfNwDfttcv5q40s8vDebUEZ0Fndivf23eMpEAJsqSdu69x9+oeZn+B4OKBBWbWADxC8MscgtOVI3i9/9eDgxTi7whOM20GlvH6l1xv/g34hpk1ElxYccdAN+7udwPfB/4U1sESggt/ulwH3BJelf3esC7/BbieYCSO1YSnekUkIxoJLpJ9zoKRbRYQfG4/G85/DFgKbDOzHeG0gzlmfJ3gGPUa8BBB1zUA3H0Z8CPgWYIE8ljg6T7WdxPwD+Bl4EUOPLuXzF/D2DcS9Dv+MRAdS763YznuvjWM8XTg9mQbcPdG4JMEdVNH8CPivsj8FQQt0WvD4+HkbsuvBD5IMPzeDuBtBMPTtaWwf/3Rr++mvvZrAH4WLv9Q+D9ZQPB+xN2bCLoMPh3W0akpfMdICsxdrewiIiIiIl3UgiwiIiIiEqEEWUREREQkQgmyiIiIiEiEEmQRERERkYghd6OQiooKr6qqijsMEZF+Wbhw4Q53r4w7jnTRsVhEhpr+HIeHXIJcVVVFdXVPI4iJiGQnMxvo3cuyko7FIjLU9Oc4rC4WIiIiIiIRSpBFRERERCIGNUE2s3VmttjMFpnZAefiwtvm/tzMVpvZK2Z24mDGIyKSa8xsmpk9bmbLzWypmX0qSRkdi0VEIjLRB/lcd9/Rw7yLgNnh4xTgl+FfERFJjw7gs+7+opmNAhaa2cPhLYu76FgsIhIRdxeLy4HfeWABMNbMJqVzAx2dCT5w0wJ+MW91OlcrIjIkuPtWd38xfN4ILAemdCs26MfihevruOhnT7JsS0M6VysiMigGO0F24CEzW2hm1ySZPwXYGHm9iQMP3JjZNWZWbWbVtbW1/QqgID+P+uZ25r/av+VERIYbM6sCTgCe6zZr0I/FRfl5LN/awIZdTf1aTkQkDoOdIJ/h7icSnL77uJm9qdt8S7KMHzDB/UZ3n+vucysr+z+M6JmHVbBwfR27m9r6vayIyHBgZiOBu4BPu3v3ZtxBPxZPGlsCwJbdzf1aTkQkDoOaILv7lvBvDXA3cHK3IpuAaZHXU4Et6Y7jomMn0d7pzFupVmQRyT1mVkiQHN/q7n9JUmTQj8Xjy4ooKcxjsxJkERkCBi1BNrOy8IIQzKwMuABY0q3YfcCHwiuoTwXq3X1rumM5dsoYxpYW8uSqnq4VFBEZnszMgN8Ay939xz0UG/RjsZkxvbyU9TvVxUJEst9gjmIxEbg7ODZTAPzR3R80s2sB3P0G4H7gYmA10AR8eDACyc8zTp81nqdX78DdCWMSEckFZwD/DCw2s0XhtC8D0yGzx+JDK8pYXbNnMFYtIpJWg5Ygu/ta4Lgk02+IPHfg44MVQ9SZh1Vy/+JtrKndw2ETRmVikyIisXP3p0jexzhaJiPH4lmVI3l0eQ1tHQmKCuIeRElEpGc5c4Q6a3YFgLpZiIjE5OjJo+lIOK9ub4w7FBGRXuVMgjytvJTp5aU8vVoJsohIHOZMHgPAks31MUciItK7nEmQAc6cXcGCtbto70zEHYqISM6ZMb6UUSUFLFaCLCJZLqcS5LMOq2BPaweLNu6OOxQRkZxjZsyZPEYtyCKS9XIqQT5t1njM4Cn1QxYRicWcKaNZvq1RZ/JEJKvlVII8trSIN0wZw1PqhywiEovDJ46irSPB5jrdMEREsldOJcgQ9ENetHE3jS3tcYciIpJzZowvA2D9Lt0wRESyV84lyGfMqqAz4Ty3dlfcoYiI5JyZlUGCvEpDvYlIFsu5BPnEGeMYX1bEnQs3xh2KiEjOqRhZzOQxJby0QRdLi0j2yrkEuaQwn4uPncSTq3bQ1qGLREREMu2UmeN57jWdxROR7JVzCTIE/ZCb2jp5cUNd3KGIiOScoyaNYseeVuqbdC2IiGSnnEyQT5s1noI84/GVNXGHIiKSc2ZWjARgzY49MUciIpJcTibIo0sKOeOwCv7+ylYSCY87HBGRnNJ1od7a2r0xRyIiklxOJsgA7zhhCpvqmlmobhYiIhk1rbyUgjzjNbUgi0iWytkE+S1HT6SkMI97F22OOxQRkUFjZr81sxozW9LD/DFm9lcze9nMlprZhwc7psL8PKaMG8GGXbpZiIhkp5xNkMuKCzj78ErmrayNOxQRkcF0M3BhL/M/Dixz9+OAc4AfmVnRYAc1vbyUDbpZiIhkqZxNkAFOn1XBprpmNuogLSLDlLvPB3obU82BUWZmwMiwbMdgxzWtvJQNO9UHWUSyU04nyCdVlQOwcL36IYtIzroeOArYAiwGPuXugz5I/IzyUuqa2jXUm4hkpZxOkI84ZBSjigt4fp0GrBeRnPVWYBEwGTgeuN7MRicraGbXmFm1mVXX1h5c97RZlRrqTUSyV04nyPl5xokzxvGC7ugkIrnrw8BfPLAaeA04MllBd7/R3ee6+9zKysqD2uisCWGCXKMEWUSyT04nyAAnH1rOqpo91Da2xh2KiEgcNgDnA5jZROAIYO1gb3TK2BGYwcY6jWQhItkn5xPkN80OWkGeXKXRLERk+DGz24BngSPMbJOZXW1m15rZtWGRbwKnm9li4FHgC+6+Y7DjKirIY9LoEjbpImkRyUIFcQcQt2Mmj2Z8WRHzX63lnSdOjTscEZG0cvcr+pi/BbggQ+HsZ2p5KRvrlCCLSPbJ+RbkvDzjzNkVPLV6J+667bSISKZMHTeCjbpZiIhkoZxPkAHOmFXBjj2tvLpdF4uIiGTKtHGlbG9sobWjM+5QRET2owQZOOvwCgDmrayJORIRkdwxrbwUd1ixtTHuUERE9qMEGZg0ZgRHTxrNo8uVIIuIZMoZh42nKD+Pv768Je5QRET2owQ59JajJ/LC+l0a7k1EJEMmjRnBCdPH8uzanXGHIiKyHyXIoQuOmYi7ulmIiGTSWbMrWLqlgR171DghItlDCXLo6EmjOWR0CY+tUIIsIpIp5xwxAYB5KzUWvYhkDyXIITPj3CMn8OSqHbR3JuIOR0QkJxwzeTTlZUW88NquuEMREdlHCXLEWbMr2NPawZLN9XGHIiKSE8yMaeNGsKVe4yGLSPZQghxx8qHlADyzRheMiIhkyiFjStha3xJ3GCIi+wx6gmxm+Wb2kpn9Lcm8q8ys1swWhY+PDnY8vakYWcwxk0fzhPrCiYhkzKQxI9i6u1l3MxWRrJGJFuRPAct7mX+7ux8fPn6dgXh6dc4RlSzcUEd9c3vcoYiI5ISp40awt62T3U067opIdhjUBNnMpgKXALEnvqk6+/AJdCacp1fviDsUEZGcML28FID1u5pijkREJDDYLcg/BT4P9DYsxLvM7BUz+7OZTUtWwMyuMbNqM6uurR3c7g8nTh/L2NJC/rF026BuR0REAlUVZQC8tmNPzJGIiAQGLUE2s0uBGndf2EuxvwJV7v4G4BHglmSF3P1Gd5/r7nMrKysHIdrXFeTncdGcSTy8bDst7Z2Dui0REYGZFWWUFObxyiaNICQi2WEwW5DPAC4zs3XAn4DzzOwP0QLuvtPdu26fdBPwxkGMJ2WXHDuJprZOnnhVF+uJyNBmZr81sxozW9JLmXPCC6WXmtkTmYwPgoaJY6eM4eWNuzO9aRGRpAYtQXb3L7n7VHevAt4PPObuH4yWMbNJkZeX0fvFfBlzysxyxoxQNwsRGRZuBi7saaaZjQV+AVzm7scA78lQXPs5Yfo4lmxuoLVDZ+5EJH4ZHwfZzL5hZpeFLz8Ztli8DHwSuCrT8SRTmJ/H+UdN4JFl23VXPREZ0tx9PtDbbeo+APzF3TeE5WsyElg3x04ZQ1tngtU16ocsIvHLSILs7vPc/dLw+X+5+33h8y+5+zHufpy7n+vuKzIRTyreeswhNLR0sGCtbhoiIsPa4cA4M5tnZgvN7ENxBHHUpFEArNjaGMfmRUT2ozvp9eBNsysZUZjPg0vUzUJEhrUCgus/LgHeCnzNzA5PVnAwRxSqGl9GUUEeK7crQRaR+ClB7sGIonzOPbKSfyzdTmdCd3cSkWFrE/Cgu+919x3AfOC4ZAUHc0Shgvw8jpg4iqVbNJKFiMRPCXIvLpoziR17Wqle11v3PRGRIe1e4CwzKzCzUuAUYrpg+pjJo1m6pUG3nBaR2ClB7sV5R06gpDCP+17eEncoIiIDYma3Ac8CR5jZJjO72syuNbNrAdx9OfAg8ArwPPBrd+9xSLjBdMzk0exuamdLfUscmxcR2acg7gCyWVlxAW895hD++vIWvnbp0ZQU5scdkohIv7j7FSmU+QHwgwyE06tjpowB4KUNdUwZOyLmaEQkl6kFuQ/vnTuNhpYOHlm+Pe5QRESGtTdMGcO40kIeWqrjrYjESwlyH06dOZ6Jo4u556XNcYciIjKsFeTnccHRh/DYihpa2nXDEBGJjxLkPuTnGe86cSqPrahhdY2GHxIRGUwXHnsIe1o7eGrVjrhDEZEcpgQ5BR89ayZFBXnc8sz6uEMRERnWzphVwaiSAh5epm4WIhIfJcgpKC8r4vyjJnL/4q106NbTIiKDpqggjzdMHaMbhohIrJQgp+htb5jMzr1tPLlap/1ERAZT1fgy1tbuIaGbNIlITJQgp+jcIyupHFXMb596Le5QRESGtZOqymlo6WCJ7qonIjFRgpyi4oJ8rjq9iidX7WD51oa4wxERGbbOml2BGcxbWRt3KCKSo5Qg98MHT5nBiMJ8/t/TakUWERks40cWc+yUMTzxqhJkEYmHEuR+GFNayDtOnMI9i7awc09r3OGIiAxbZ82uYNHG3TS2tMcdiojkICXI/XTV6VW0dSS4d9GWuEMRERm2zppdSWfCeWxFTdyhiEgOUoLcT4dPHMWcKaO5W3fWExEZNCdVlTN13Agda0UkFkqQB+AdJ0xl8eZ63VlPRGSQ5OcZZ82u4MX1dRruTUQyTgnyALztuEnkGfzlRbVsiIgMlhOmj6OhpYN1O/fGHYqI5BglyAMwYVQJ5x4xgdtf2EhLe2fc4YiIDEtHTxoNwNItGlpTRDIr5QTZzMaZ2clm9qaux2AGlu2uPutQdu5t40/Pb4g7FBHJEWZWYmafMbO/mNldZvYfZlbSxzK/NbMaM1vSR7mTzKzTzN6d3qgH7vCJoxhZXMAza3bGHYqI5JiUEmQz+ygwH/gH8PXw73WDF1b2O23meN44Yxy3PqcEWUQy5nfAMcD/AtcDRwG/72OZm4ELeytgZvnA9wmO7VmjqCCP02eN58lVGg9ZRDIr1RbkTwEnAevd/VzgBCCnj1hmxtuPn8yqmj2s3KaL9UQkI45w96vd/fHwcQ1weG8LuPt8YFcf6/134C4g68ZUO/nQcjbVNVPT0BJ3KCKSQ1JNkFvcvQXAzIrdfQVwxOCFNTRcOCe4WE/DEIlIhrxkZqd2vTCzU4CnD2aFZjYFeAdww0HGNihOmD4OgOfX9ZXji4ikT6oJ8iYzGwvcAzxsZvcCOX+njMpRxbz5qIncWb2R9s5E3OGIyDBlZovN7BXgFOAZM1tnZq8BzwIHez3IT4EvuHufVxyb2TVmVm1m1bX/v737Do+qTPs4/r3Tey+EEAglVOmhSRNBBHUt2DuKa10X1rbuq+u6rr6KZXV9rajYKYpdFgQLIJ3QO4QeQkkjIb3M8/4xQ5iQAoGZzCS5P9eVK5MzZ8755czMM/ec85znZDTMQcSerUKJCvJh9oZDDbI+pZQC8DqTmYwxV9luPi0ivwGhwFynpWpEbuifwLwtR/h121Eu7tbC1XGUUk3TZU5cdjIwQ0QAooBLRKTcGPPtqTMaY6YAUwCSk5MbZHBiL08PLu0ex4xVBzheXEawn3dDrFYp1czVe5g3Y8xCY8z3xphSZwRqbIYlRRMb4ssny/a6OopSqokyxuyr6+ccl93WGJNojEkEZgH311Qcu9LlvVpSUm5h/pYjro6ilGomdBzkc+Tl6cGEIW1ZkprFpoO5ro6jlFJViMh0rF0xOolImohMEJF7ReReV2c7U70TwokP82fW6jSM0avqKaWcTwtkB7i+X2v8vT2Zsmi3q6MopVQVxpgbjTFxxhhvY0wrY8wHxph3jDHVTsozxow3xsxyRc66eHgIdwxOZOmuLN5asMvVcZRSzYAWyA4Q6u/NLQNbM3vjIQ5kF7o6jlJKNTkThrRl7Hkt+M/PO8ktLHN1HKVUE6cFsoPcOaQtAlgKXu8AACAASURBVExdssfVUZRSqskREe4Z3p7SCgs/bTns6jhKqSZOC2QHiQv154pe8UxfuZ/sAj1/USmlHK1nq1Diw/yZu0kLZKWUc2mB7ED3Dm9HcZmF6Sv18tNKKeVoIsIl3Vvw+84M8oq1m4VSynmcXiCLiKeIrBWRH2u4z1dEZopIqoisEJFEZ+dxpqTYYAa1i2Taiv1UWPRMa6WUcrQx58VRVmH4davbXRVbKdWENMQe5InA1lrumwDkGGM6AK8Ckxsgj1PdMrANB48VsXCHNt5KKeVovRPCiAn25fVfdmp3NqWU0zi1QBaRVsClwPu1zHIF8LHt9ixgpNgu59RYje4WS3SwLx8vPaex+5VSStXAw0N49fpeHDxWxD2fpmDRo3VKKSdw9h7k14DHAEst98cDBwCMMeVALhB56kwicreIpIhISkZGhrOyOoS3pwd3DE5k4Y4Mlu7KdHUcpZRqcgZ3iOKZK7qxam8OU37X8eeVUo7ntAJZRC4DjhpjVtc1Ww3Tqu0OMMZMMcYkG2OSo6OjHZbRWe4c3JaYYF/eXagNt1JKOcN1yQlc0CmatxfsorC03NVxlFJNjDP3IA8GLheRvcAM4EIR+eyUedKABAAR8QJCgWwnZmoQft6e3DKwDQt3ZJB6NN/VcZRSqskRER68MIncojLeW6TjzyulHMtpBbIx5m+2y5omAjcAvxpjbjlltu+B2223r7HN0yQ6lN00oDU+nh58tFQbbqWUcoa+bcIZmhTFqz/vYH+WXsVUKeU4DT4Osog8IyKX2/78AIgUkVTgIeDxhs7jLFFBvvyhZ0u+XnOQ0vLaumArpZQ6F09f3g0/bw+e/G6Tq6MopZqQBimQjTELjDGX2W4/ZYz53na72BhzrTGmgzGmvzGmSXXavbhbLIWlFaTsbfS9RpRSyi21jw7i1oFtWLQjgy3pea6Oo5RqIvRKek40JCmKYF8vPtcr6ymllNP8aUQS4QHe/O2bjeSX6Al7SqlzpwWyEwX4eHHH4ERmbzjEDC2SlVLKKUIDvHn2yu6sP3CMh79Yp2MjK6XOmRbITjZxVEcGtotg8txt5BaWuTqOUqqZEZGpInJURGrspCsiN4vIBtvPUhHp2dAZHeHSHnH8dUxnftp8hKW7slwdRynVyGmB7GSeHsLfL+vKsaIy3vhtp6vjKKWan4+AMXXcvwcYbozpAfwLmNIQoZzh1kFt8BBYsUcLZKXUudECuQF0axnKtX1b8dHSvaTl6FBESqmGY4xZRB3jyxtjlhpjcmx/LgdaNUgwJwjy9aJvm3A+XrqXI3nFro6jlGrEtEBuIJNGdcRi4JV5O1wdRSmlajMBmOPqEOfi2Su7U1RWwWs/a1urlDp7WiA3kJZh/kwY0pZv1h7k27UHXR1HKaWqEJERWAvkv9Yxz90ikiIiKRkZGQ0Xrh46tQjm1oGJzFx1gL2ZBa6Oo5RqpLRAbkAPXdSRNpEBTJq5jl+2HnF1HKWUAkBEegDvA1cYY2rtwGuMmWKMSTbGJEdHRzdcwHq6Z3g7AKYu2UMTuTirUqqBaYHcgPy8PfnhwSHEBPvy+i87teFWSrmciLQGvgZuNcY0iX4JsSF+XN2nFZ8s28es1WmujqOUaoS0QG5gIX7ePDgyifVpuWw8mOvqOEqpJk5EpgPLgE4ikiYiE0TkXhG51zbLU0Ak8JaIrBORFJeFdaDnx3WnbVQgX6QccHUUpVQjpAWyC1zRqyUBPp7852fdi6yUci5jzI3GmDhjjLcxppUx5gNjzDvGmHds999ljAk3xvSy/SS7OrMjeHl6cHWfeFbtzeHgsSJXx1FKNTJaILtAiJ83k0Yl8cu2oyxOzXR1HKWUapIu69ESgB/Xp7s4iVKqsdEC2UVuPz+RsABvZqzSw39KKeUMiVGB9EwI441fU0nZW+tQ0EopVY0WyC7i6+XJH3q0ZPaGQyzRvchKKeUUr17Xk2A/L16cu93VUZRSjYgWyC70yOhOADz65XoKS8tdnEYppZqedtFBXN+vNSv3ZvM/32zkQLZezVQpdXpaILtQaIA3X947iPTcYt7/fY+r4yilVJN019C2jOsTz6yUNG6YspzcwjJXR1JKuTktkF2sX2IEIzvHMHXJHo4Xa6OtlFKOFujrxb+v68W0Pw7g6PFi7vksRUcQUkrVSQtkNzBpVEeOFZYxdfFeV0dRSqkmKzkxgodHd2L57mxGvrKQDWnHXB1JKeWmtEB2A91bhXJR11jeX7ybY4Wlro6jlFJN1l1D2nL7oDbszixg4ox1HD1e7OpISik3pAWym3hkdCcKSsp5eZ6eaa2UUs7i5enBP684jw9uT2Z/diF//Fi7WyilqtMC2U10ahHMbYMSmbZiP9sPH3d1HKWUatJGdonluSvPY31aLvd+tlpHElJKVaEFshuZODKJUH9v7v40hdSj+a6Oo5RSTdr1/RJ48MIOzNtyhKe/3+zqOEopN6IFshsJD/Thg/H9yCsq497PVlNeYXF1JKWUarJEhIdHd+LmAa35es1B8kt0L7JSykoLZDfTp3U4z4/rQerRfG6YshyLRfvGKaWUM13TN4Fyi+FP09Zof2SlFKAFslsa3TWWVuH+pOzL4Zkft2iDrZRSTtQrIYy7h7VjwfYMft561NVxlFJuQAtkN+ThIcydNIzhHaP5aOle3vwt1dWRlFKqSXt4dEeignyYtmKfq6MopdyAFshuKsjXiw/H9+PSHnG8Mn8Hq/dluzqSUqoREpGpInJURDbVcr+IyOsikioiG0SkT0NndAe+Xp6M69OK37ZnsHx3lqvjKKVcTAtkN+bhIbwwrjuRgT6M/3AVaTmFro6klGp8PgLG1HH/WCDJ9nM38HYDZHJLD13UkYhAH56bvZXfth3Vk/aUasa0QHZzwX7eTPvjQCoshiGTf6PrU3Pp+tRcJs/dxsa0XFfHU0q5OWPMIqCuQ1BXAJ8Yq+VAmIjENUw69+Ln7cnEkUlsSs/ljo9Wcf7zv/Di3G16srRSzZAWyI1Ax9hgPp3Qn+uSW9E6IoDC0greXrCLP7yxmK/XpOlJfEqpcxEPHLD7O802rVm6/fxEUp4YxZOXdqF9TBBvLdjFqr3axU2p5kYL5Eaib5sIXrymJ3MnDWPm3QN58MIOADz0xXru+2wNmw7q3mSl1FmRGqbV+K1bRO4WkRQRScnIyHByLNeJDPLlrqHt+HB8Pzw9hMlzt7EvqwCAA9mFpB7Vq50q1dR5uTqAqr8B7SIZ0C6SiSOTmDhzHbM3HGJxaibzHxpGXKi/q+MppRqXNCDB7u9WQHpNMxpjpgBTAJKTk5v8oauwAB8euKA9r/+aygUvLyA+zJ+0nCIAbhvUhsfGdCbIVz9GlWqKdA9yI+bl6cGbN/VhzsSh5JeUM+j5X0l8fDa/bjvi6mhKqcbje+A222gWA4FcY8whV4dyF5NGdeS925Lp3CKE8grDqC6xBPl68cmyfQyd/Ctv/pZKblGZq2MqpRxMnNV/VUT8gEWAL9Y91bOMMf84ZZ7xwEvAQdukN4wx79e13OTkZJOSkuL4wI3c5yv28cQ3J0dxeuXannh4QNuoIHolhLkwmVIKQERWG2OSXbDe6cAFQBRwBPgH4A1gjHlHRAR4A+tIF4XAHcaY0zayza0tPnGinoeHtUfK419tYMaqk123tz4zBn8fT5dkU0qdmfq0w84skAUINMbki4g3sBiYaDtL+sQ844FkY8yfznS5za1Rrg+LxXDwWBGXv7GYnMKTezQ+v2sAfVqHa+OtlAu5qkB2lubeFpeWW1i9L4ev1qQxa3Uaj4zuyJ8uTHJ1LKVUHerTDjut85SxVt75tj+9bT9Nvs+aK3l4CAkRASx4ZASv/ryD/JJylu/O4ub3VwDw1X3n07dNuItTKqVU4+fj5cGg9pEMah/JscIy3l6wi0HtI+nTOhzr/iGlVGPm1D7IIuIpIuuAo8B8Y8yKGma72nb1plkiklDD/c3mzGlHCQ3w5unLu/HytT35cHw/usaFAHD120v5IuUAFTqmp1JKOcyDF3bAQ4Sr317GE99uorzC4upIdcovKae4rMLVMVQTVtAEXmNO62JRZSUiYcA3wIPGmE120yOBfGNMiYjcC1xnjLmwrmU198N6Z2tpaiY32fYktwz1Iz23mFev78lVvVu5OJlSzYN2sWja8kvKmTh9Lb9sOwrALw8Pp310kEuyGGN4dvZWPli8h7AAb44VlvG3sZ3p2jKEzPwS/jJzPWA9V+Xqvq1YuCODAB9PuseHAtYLpihVH8cKSwny9cLL0wOLxdD+if8SHeTLF/cMolW4P0VlFfh4eXAgu5C2UUF4erjmKItb9EGutiKRfwAFxpiXa7nfE8g2xoTWtRxtlM/e/qxC5mw6xILtGSzbnQXAzLsH0rt1OD5eOqCJUs6kBXLTV1hazpVvLmHHkXxahPix4NELKLcYyisshAX4OH39JeUVPDhtLfO2nP1IRt1ahjD7z0MdmKqqwtJy3v99DyM6xdC9VZ0f943O8t1ZJEYG0iLUz9VRGlROQSnDX/qNvOJyRnWJJb+kjOW7a7+4Tt824YT6e3M4t5gv7x1EgI9ng3VLcosCWUSigTJjzDER8QfmAZONMT/azRN3YjghEbkK+KsxZmBdy9VG2THyS8rp88x8SissdIwN4tkruzN/y2H8vT2ZNKpj5ZnaJxhjKC6znNWJfrmFZfh6e1TulcguKMXHy0PHD1XNihbIzceyXVnc+N5yRnWJJWVfNscKyxiaFMWUW5OdcrL0poO5+Hh5cO07yyqHnLt9UBv+dkkX/Lw9eWDaGmZvsI7c5+vlwby/DGNzeh73f74GsBYsqUfzKx+7+Z8XE1hL+2yM4Z2FuxnSIaqywC2rsLA5PY9X5+/gfy7pQqcWwTU+Nq+4jBunLGdzeh4eAlueGVNtb7Uxhi9XpzGkQxQtw/yxWEy1zyN39MavO3l53g6C/byYde/5vLtoF1+vsQ7Q9eLVPbg2uVW9i8Ci0gr2ZRfQuUWIMyI7zOS523h7wa4q0y7tHsfsjWc2WuQtA1vz7JXdnRGtGncpkHsAHwOeWPs6f2GMeUZEngFSjDHfi8jzwOVAOZAN3GeM2VbXcrVRdpx5mw/z4k/bST2aX2X638Z2pk+bcLak51FuMWxMO8bm9Dx2Hs1n8tXdub5f6zNex+KdmdzywQq8PYX/3NCb6Sv38/vOTNpFBfLfiUP1UN5Z2JKex+LUDAa1i6JzXDB7Mgvw9vQgIsCH0ABvh6zDGMOPGw7RItSPfokR9X7sCdkFpUQG+QLWxj4tp5DFqZl4eQiX9mhJRKDz96qdq7ScQhbtyOSCTtG0DDv7C/Fogdy8PD9nKx8u3gtAqa1PcqCPJ3Fh/iSE+xMZ5EtSTBCX9oijVXjAGS8343gJV721hKN5JQzrGMWKPdkcLy6vvH9oUhQf3N6vylFBi8VQUm7hWFFplYtJ7c8qpLTCQoeYIIwxfLcunUkz1/HghR14eHSnauvem1nABS8vOG3GPc9fUqUYLC23MHnuNqYu2YMxkBDhz4HsIqaOT+bCzrGV8+04cpxdR/O5z1a43zygNbM3HuKXh4ZXtiMN5WheMYdyi+lpGyb1xw3pWAxc3rNllfn2ZRXw9ZqD/OeXnXUub2TnGF64ugfRwVX/j5V7srn70xR8vTz47ZELyC4o5aWfthMT7Mt7v+8BYNUTo6o97lTFZRV8tHQvX61OY392IV/cM4gZq/bTtWUotw5sU99//4ztySxg5CsLGNs9jqt6xVNusWAxcEn3uMp5Vu/L5lBuMUM7RLM3q4Ar3lxC17gQ2kUH8qPti9viv46o9X1w4ktZ26hAxpzX4pzyukWB7CzaKDvejxvSefO3XTxxSRcmz93GxhouW50Q4c+xgjKOl1gb4st7tuS2QW1IToyo8Rv+tsN5XP3WUgpKq3fS9/P2oLjM+oFx/wXteWxMZ4yxDlFXnw+K5mjupkM8MG1trSdaXtojjtFdY1mamsXY7i24oFNMtXlKyy14egieHsLOI8d54ptN7M4sIDbEl8O5xUQH+9KnTTjTVuwH4K4hbXni0i6n3fthjGHy3O28s3AXwb5edIsPYfnubO4e1o47Bicy7q2lHMotrpzf21N49srz2J9dyIC2kfRvG8HCHRkkxQTRzgV9N2t6HVsshsvfXMymg3mV0766bxB929TvSwNogdwclZRX4OPpQUm5hbX7j/HazztYsaf6oedpdw3g/A5RgPU1J0K195sxhoU7Mpg8dztbD518PUYG+jC2ewtC/LwJ9vPm3uHtzvpwtcVi6Pz3uZRWWPjvn4fSJS4YEaGswsLFry5id6b1ctsnzmOx1yshDAOsP3AMABGYOr4fS3Zm8v36dI4eLwHgoYs6cs/wdvR5Zj5FZRW8MK4HK/Zks2JPVuVVCk8VFeTDS9f25MW527nj/ESu65fApoO5zFqdxuiusZXbzhHKKyxMW7mfp77bDMCiR0ewdFcmj3+9EYD5fxmGiPUaA1MW7WbyXOs+PR8vD76693z+8MZiAP73qu6M6xPPB4v38NJP2wHrF6RHL+7EkKRoAnw8eeaHLczdfLhy3c+P686PG9JZkppVJdO/rjyvSpFbXmFtw0Wkyl752jx6cSd6JYTRo1Uoi3ZkcvCYtc3teY7XSCivsHDHR6tYsSebJX+98LRFfE0OZBcy9MXfuO+C9gztEMWg9pGICMcKSwn190ZE+HbtQSbNXAfAtw8MPqdrO2iBrM5aWYWFf3y/mTX7cnji0i5sP3ycK3rFEx3sy8FjRfzP1xtZuKPqSCKtwv2ZODIJL0+ha1woHyzezRcpaQC0iw7klWt7ciCniJ82H2ZYUhTX92vN2wt2VTYsYC2WyioME0cm0SLUj8t6xBHs583cTYfJKSzlqt7x1fY2n+hTfVXveI4VldE+uvaO/6XllkbTz9piMcxancay3Vm0iwokIsiHpbuyWLQjg+PF5bSJDODFq3sw4eMU8kvKCfHzoqC0osai+cELO3B1n1bkFpXh6SHMXHWAn7cewctTOJB98sMo2M+ryl6oEzrFBrP9yHFahvrx9f2DaRHqx/wtRwj196Z/2whSj+YT4udFWIAPr8zbzruLduPj6VG5x8yeCAxLiqZHq1BahPpVubCNvZahftzYvzXvLtrNqC4xvHJdLzw9hNkbDvHWglQARnWJZUTnGIdcBOdQbhHXvL0MizG8fG1P8kvKmbPxEN+uS8fXy1rcjOwcU3nyVafYYOZOGlrvIkQLZGWMYdvh47QI8eOz5fvw8BBe+mk7of7eLHz0Av7v11S+WXuQbi1D8PIQFqdm0r9tBHlF5ZU7LsICvHn5mp4kRgUQGehLqL+3Q7sgrN6Xw9VvLwWgS1wIH93RjwH/+0vl/R+O78eIzjHsysgnzN+borIKooN98fXypKi0giveXMyOI/nVljv2vBZMGtWxsvvFK/O283+/ptaYYebdAzmQU8Sny/cRF+JXpYgE+OTO/jw6az1H8kqIDfFl0WMj8PU686OR2QWlPP7VBoZ3iubmAVX3rj7zwxamLtlT7TGRgT5kFZTWuLyeCWG8fE0PkmKDOZpXjK+XZ5WjeYWl5Vz6+mL22L5g2IsI9OGLewby2KwNrNlv/XJx26A2hAX4MK53POPeXkp2QSmf3NmfQF8vvDyE695dRr/ECP55RTdGvrKwyvLW/P0iDmQXMnHGWsafn8jTP2ypdTvU1ZXmdDKOl/DUd5uYs+kwz111XrXtWB/j3lpS+b8nRgYQH+7PktQsOrcIZsbdAxn2orV/s73Hx3bm3uHt670uLZCVU6XszabCYnjoi/UcPFbzN36Acb3j+ff1vWq9f83+HCbOWMvBnCJOre2GJkXRNiqQT5btAyA8wJufJg0jJsQPYwwzVh3gudlbyS85+abp3CKYG/olUFBawf0XtK8sYL5YdYAnv91UWbRd1DWW5buzKgvCl67pwbXJNY4wWC/GGJbuyqJvm/B6dx0xxrBiTzaHc4tZvS+HT5fvq3G+3q3DePvmvrQI9WPnkeOIQIcY6wfO83O2grGOhx3k61W51+J0buyfwNOXdyOnoIz03CKCfL2YungP9wxvT0SgD49/tYE5mw4TFeRDl7gQft+ZCZz8UmOvU2wwcyYOZX3aMX7eeoTRXVvw/uI9/LTpMO/e1pcRdnu0Z284ZP3S1DGa/dmF/L4zgyO5xdX2THl5COW17DEf3CGS3gnh/OnCDrVu84zjJbyzcBe3D0qkdeTJIxTbDudx54erqq3vVNf0bcWLV/cgp7CUNfuP0blFMAkR9T/SoQWyqsmXKQd4dNaGM5q3d+swpv9xoNO7pi3YfpTxH66qMi3Qx5M5E4dVeQ/VxGIx7MrIZ8uhPCbOWMf1yQmM6hrLRV1jq8xnjOGDxXvIKijlnmHt+GF9OsM7xlRbvjGG6SsPsC+rgP5tI5jw8cnX3ONjO/PCnG2c3z6Sz+8aQOrRfLILSumZEFZlG+3OyCctp4hB7SMxBiZ8vKqyHXt+XHcGtYvk0VnrWbv/GOUWQ3SwL5/fNYA/T1/LtsPHaRMZwPy/DGfN/hwe+HxNZaF8UddYHru4Ex1igk77hbmkvIJR/17IgewiIgJ9yC4o5e+XdeXOwYmICEfzinni202c1zKUB0a0x8vTukPnyW838tny/XUu+5aBrenbJpzwAJ9qRw13HjnOQ1+sJzkxnO/WpXNh5xjKKix8ty6de4a142+XdKlz2TU5klfMmNcWkVNYVrntz+Ukux/Wp/Pg9LV1zvPCuO6s2mu9MM8Jn07oz9Ck6HqtSwtk1SAOZBeyfHcWo7u1YF9WATNWHWDaiv0MTYrixWt6EBnoe0Z7bQtKyskuKCW3qIx/fL+ZQF8vFtn2UkcE+nBDvwTeWrCLG/sncPv5iczbfIR/z99By1A/usSFsHRXFkWnjLf42JhOXN6zJePeWlp5aO9Ev7cTPD2kcq/rVb3jSQj3p01kIFf3tQ59tyHtGNsOHad7q1B+35nB2PPiKgujQ7lF3P/5GkZ1ieXmAa0JC/CpPEnjjsGJ/OMP3c54O1ZYDB8v3cszP578ph8W4M2aJy9iyu+7CfTx5IJOMQT6etWrz26FxTBr9QGKyywcySvm6zUHeeTiTozrHc+HS/cSE+zLsKToM+q3vHBHBpNmrK28QuON/ROYvtJ6mV0PAYuBuFA/Pp3Qv7JgP8EYQ15xOaH+p1+PMYZ/z9/B6n05TB3fj3/+sIXpK60fDj1ahfL8uO6E+nuzOT2PT5btrTwUGR/mzxs39Wbaiv30aBXKrYMS2ZKex9+/28S+rAIy80sr+8GP7hrLv37cwsfLTn4JefiijtzQvzXztxxhy6Fc+rQO57z4UHy9PGgTGXhG2/t0tEBWNTHGcM07y1i9L4ehSVG8e2tfKiyG33dmcmHnGDak5bIrI5+0nEL+PDKpXntKz0VOQSm9/zUfsLan91/QoUHWezqfLd/H37/bxJ8vTOIvF3Xk+TlbeXfh7irz9EoI45v7zwfg6zUHefjL9dWWM6xjdOXnjL2oIF8W/3UEft6eGGMoKqsgwKf6XtZdGfm0iwqsV2F4rLAUDw8hxO/MzxWxWAw3vLeclXZdc967LZnN6bks3ZXFX8d0qneXL2MMj3+1kZkpB5g7aWi9TgL85w+b+XDJXgD6t43gjZt6ExPsmFE7yissfLUmjc4tQugQE8SVby5h59F8LusRx+s39MbDQ7BYDJvSc3n0yw3cM7wd4/rUb6haLZBVo1ZWYeH7dekE+XlxQadofL08q7wpwXq47vUbe+Pt6UF+STleHsJPmw8T6u/N//2ayup9OZXzjusTzyOjO9EyzJ/isgosxpBfUk5MsB+FpeX868etlUUYwAMj2tOtZWjlGd72bh3Yhqv6xHPTe8sr+1EDtAjx43DeyT2R/RMjCPbz4sGRSbV2A6iwGKYu3sNz/90KWPcO9W8bwS9bj/LPy7sx2IH96hzBYjGk5xYRFeRbZe9MhcWwOT2X7vGhDh+qxxjDodxiSsstRAb5EGz3wVJabmFJaiY/rE/nm3UHsW/KzosPobC0gt0ZBcQE+3KjrfjdejivynxPXNKF6GBfLusRV7nHxlm0QFa1KSwtJ7eorMoJdO5g08FcIoN83C5XxvESIgN98PAQyiss9PnX/GqH4Pu3jaBby5DKz40AH08KbefEXJ+cwPPjuvPxsr3868cthAf48ORlXfDz8iQ5MeKs+tI6U3ZBKUfyiukQE8TezAKSYmseJaQ+Mo6X0O+5n+kUG8zsPw+p1v4dyi0iOsi3cnpmfgn/98vOyh0L0/44gPPbu9dn1JnQAlk1OTkFpdz43nK2HT7On0Z04MGRHWrdk3K8uIxbPlhJXlEZk0YlcUWv+NMuf9muLI7kFfPB4j1VTlKMCfalbVQgvRLCeHdR1b0UA9pGkJwYzpu/WYe3iQryZc7Eofxl5joWp2ZWzvfCuO7c0L812w8fZ+XebAa3jyTYz5vxH66sPLFiXO94/ndcdx3V4yyt3pfDHR+u5KGLOvLd+nTW2vqzvXVzH8ae1wIRoai0gie/3cRXa9Lo0zqMd27t67A9H2dCC2SlnCe7oJQAH0+8PITbP1xZeXRpWMdopt6eTIUxFJRUVDsKV9tJkc3B+A9XsmB7Bj6eHnx9//mcFx9KbmEZL8/bXtnNb9KoJK7vl8DFry6q/BJS373O7kQLZKXOUkl5Bc/+uBV/H09Gdo5hQLvIyvsqLIYNaceYveEQmfklPHdVdwJ9vSguqyDjeEll9wtjDBsP5pKZX8I7C3azcm925dWsThUV5Mt//zyEmJDmNbC8MxljOFZYhr+PZ41fODLzS4hq4CGjQAtkpRpKcVkFf/wkhdYRATw+tnOVI0/qpLIKC13+PpdyiyE8wJshSdH8sD691vmHdIji+XHdz+ocDHehBbJSbiIzv4TkZ38GrP1kuRi6GQAACqJJREFUx5+fyPYjx9l++DiPj+3sdt0olPNogayUcjdH8or5YX06z87eWjntxFjaezILuPi1RYT6e/PGTb3rfUKcO6pPO6yXMlPKiaKCfPl0Qn9KyiyMOuUsbqWUUsqVYkP8mDCkLRZjSIoJZkTnk6NgdGoRzN4XLnVhOtfSAlkpJ2sK37pV4yYiY4D/YL2y6fvGmBdOub811iufhtnmedwY898GD6qUanAiwt3D6j+mcFPXOK6coJRS6qyIiCfwJjAW6ArcKCJdT5ntSeALY0xv4AbgrYZNqZRS7kULZKWUatr6A6nGmN3GmFJgBnDFKfMY4MRp6aFA7WfqKKVUM6BdLJRSqmmLBw7Y/Z0GDDhlnqeBeSLyIBAIjGqYaEop5Z50D7JSSjVtNQ3weurwRTcCHxljWgGXAJ+KSLXPBxG5W0RSRCQlI6P6VciUUqqp0AJZKaWatjQgwe7vVlTvQjEB+ALAGLMM8AOqjUFojJlijEk2xiRHR+vJp0qppqvRjYMsIhnAvrN4aBSQedq5nM9dcoD7ZNEc1blLFnfJAe6T5WxztDHGNHhVKSJewA5gJHAQWAXcZIzZbDfPHGCmMeYjEekC/ALEmzo+IM6yLXaX5xDcJ4u75AD3yaI5qnOXLO6SA84uyxm3w42uQD5bIpLiDoP0u0sOcJ8smqM6d8niLjnAfbK4S476EJFLgNewDuE21RjznIg8A6QYY763jWrxHhCEtfvFY8aYeU7I4Tbbzl2yuEsOcJ8smqM6d8niLjnA+Vn0JD2llGribGMa//eUaU/Z3d4CDG7oXEop5a60D7JSSimllFJ2mlOBPMXVAWzcJQe4TxbNUZ27ZHGXHOA+WdwlR2PkTtvOXbK4Sw5wnyyaozp3yeIuOcDJWZpNH2SllFJKKaXORHPag6yUUkoppdRpaYGslFJKKaWUnSZfIIvIGBHZLiKpIvK4E5afICK/ichWEdksIhNt058WkYMiss72c4ndY/5my7NdRC52ZFYR2SsiG23rTLFNixCR+SKy0/Y73DZdROR12/o2iEgfu+Xcbpt/p4jcXs8Mnez+73Uikicikxpqm4jIVBE5KiKb7KY5bBuISF/bNk61PbamK5XVluMlEdlmW9c3IhJmm54oIkV22+ad062vtv+pHlkc9nyISFsRWWHLMlNEfOqRY6Zdhr0iss7Z20Rqf982+OukuTib93I9l69tcfUM2hbXnqPB2+JacjR4O1xHFm2L7RljmuwP1jE/dwHtAB9gPdDVweuIA/rYbgdjHZC/K/A08EgN83e15fAF2tryeToqK7AXiDpl2ovA47bbjwOTbbcvAeZgvRTtQGCFbXoEsNv2O9x2O/wcnoPDQJuG2ibAMKAPsMkZ2wBYCQyyPWYOMLYeOUYDXrbbk+1yJNrPd8pyalxfbf9TPbI47PnAehW2G2y33wHuO9Mcp9z/CvCUs7cJtb9vG/x10hx+6nrtOHAd2haf/jnQtvjktAZvi2vJ4bDngjNsh2vLcsr9zb4tbup7kPsDqcaY3caYUmAGcIUjV2CMOWSMWWO7fRzYCsTX8ZArgBnGmBJjzB4g1ZbTmVmvAD623f4YuNJu+ifGajkQJiJxwMXAfGNMtjEmB5gPjDnLdY8Edhlj6rrilkO3iTFmEZBdwzrOeRvY7gsxxiwz1nfeJ3bLOm0OY8w8Y0y57c/lWC/7W6vTrK+2/+mMstShXs+H7dv4hcCs02WpK4dtOdcB0+sK54htUsf7tsFfJ82EtsUn16ltcTNti92lHT5dFm2LrZp6gRwPHLD7O426G8xzIiKJQG9ghW3Sn2yHAKbaHV6oLZOjshpgnoisFpG7bdNijTGHwPpiBGIaKAvADVR9k7lim4DjtkG87bYjMt2J9dvsCW1FZK2ILBSRoXb5altfbf9TfTji+YgEjtl92JztNhkKHDHG7LSb5vRtcsr71h1fJ02BtsVW2ha753vM1W2xO7XDoG0x0PQL5Jr6mThlXDsRCQK+AiYZY/KAt4H2QC/gENbDFXVlclTWwcaYPsBY4AERGVZXbGdmsfV/uhz40jbJVdukzpj1XLejts0TQDnwuW3SIaC1MaY38BAwTURCHLW+Wjjq+XBUxhup+gHu9G1Sw/u21llrWacrX7uNibbF2hafNmY9191U2mJ3a4dB22Kg6RfIaUCC3d+tgHRHr0REvLE+sZ8bY74GMMYcMcZUGGMswHtYD4vUlckhWY0x6bbfR4FvbOs9YjvMcOKQyNGGyIL1g2GNMeaILZNLtomNo7ZBGlUPxdU7k+3kgcuAm22HfLAdRsuy3V6NtY9Zx9Osr7b/6Yw48PnIxHqYy+uU6WfM9thxwEy7fE7dJjW9b+t4fIO/TpoYbYu1LT7Bbd5j7tAWu1M7bMuvbfEJ5iw6+zeWH8ALa0fttpzszN7NwesQrH1aXjtlepzd7b9g7UsE0I2qHe93Y+10f85ZgUAg2O72Uqz91V6iamf3F223L6VqZ/eV5mRn9z1YO7qH225HnMW2mQHc4YptwiknFThyGwCrbPOe6PB/ST1yjAG2ANGnzBcNeNputwMOnm59tf1P9cjisOcD654p+5ND7j/THHbbZWFDbRNqf9+65HXS1H/qeu04cB3aFteeR9tiN2mLa8jhkna4pix220XbYmOadoFs2zCXYD0rchfwhBOWPwTr7voNwDrbzyXAp8BG2/TvT3kTPGHLsx27synPNavthbve9rP5xDKw9k36Bdhp+33iRSPAm7b1bQSS7ZZ1J9aTAlKxa1jrkSUAyAJC7aY1yDbBemjoEFCG9dvjBEduAyAZ2GR7zBtgvSLlGeZIxdpP6sRr5R3bvFfbnrP1wBrgD6dbX23/Uz2yOOz5sL32Vtr+vy8B3zPNYZv+EXDvKfM6bZtQ+/u2wV8nzeWntteOA5evbXHNWbQtdpO2uJYcDd4O15bFNv0jtC3GGKOXmlZKKaWUUspeU++DrJRSSimlVL1ogayUUkoppZQdLZCVUkoppZSyowWyUkoppZRSdrRAVkoppZRSyo4WyKpZEpFJIhLg6hxKKdWcaVus3JUO86aaJRHZi3X8xExXZ1FKqeZK22LlrrxOP4tSjZuIBAJfYL3EpCfWAdRbAr+JSKYxZoSIjAb+ifWqRbuwDjKeb2u8ZwIjbIu7yRiT2tD/g1JKNXbaFqvGRLtYqOZgDJBujOlpjDkPeA3rtdhH2BrkKOBJYJQxpg+QAjxk9/g8Y0x/rFfgea2BsyulVFOhbbFqNLRAVs3BRmCUiEwWkaHGmNxT7h8IdAWWiMg64Hagjd390+1+D3J6WqWUapq0LVaNhnaxUE2eMWaHiPTFen3350Vk3imzCDDfGHNjbYuo5bZSSqkzpG2xakx0D7Jq8kSkJVBojPkMeBnoAxwHgm2zLAcGi0gH2/wBItLRbhHX2/1e1jCplVKqadG2WDUmugdZNQfdgZdExAKUAfdhPTw3R0QO2fq+jQemi4iv7TFPAjtst31FZAXWL5S17dlQSilVN22LVaOhw7wpVQcdgkgppVxP22LV0LSLhVJKKaWUUnZ0D7JSSimllFJ2dA+yUkoppZRSdrRAVkoppZRSyo4WyEoppZRSStnRAlkppZRSSik7WiArpZRSSill5/8Bu0CeD/gaF+IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(10,6));\n",
    "ax1 = plt.subplot2grid((2, 2), (0, 0), colspan=2);\n",
    "ax1.plot(losses);\n",
    "ax1.set(xlabel=\"step\", ylabel=\"loss\", title=\"ELBO\");\n",
    "ax2 = plt.subplot2grid((2, 2), (1, 0));\n",
    "ax2.plot(a);\n",
    "ax2.set(xlabel=\"step\", ylabel=\"a\", title=\"Mean Parameter\");\n",
    "ax3 = plt.subplot2grid((2, 2), (1, 1));\n",
    "ax3.plot(b);\n",
    "ax3.set(xlabel=\"step\", ylabel=\"b\", title=\"Standard Deviation Parameter\");\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now this has converged, we can use the optimised values of `a` and `b` to sample from the guide. We can also compare this to the analytic speed, which can be calculated for this simple example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean = 3.40, Standard Deviation = 0.89\n"
     ]
    }
   ],
   "source": [
    "print(\"Mean = {:.2f}, Standard Deviation = {:.2f}\".format(np.mean(a[10000:20000]),np.mean(b[10000:20000])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAEICAYAAACdyboFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XmcFNW5//HPwwwwLMMOguwqKhAjUdAEo0FFL+5LYiCJiUQT4pYbveaXqzFq3BITTTQ3atyDiRqXmAUJcRf3BYiIIqCsMizDAMMM2wAzc35/VNXQDN0z1d3V09Pd3/fr1a/urq566vRoH54659Q55pxDRERERNLXJtsFEBEREckXSqxEREREIqLESkRERCQiSqxEREREIqLESkRERCQiSqxEREREIqLEKo+Y2UVmVm5mW8ysZ7bLkwwzm29m41I89t9mdl7ERcpYXBHxmFkHM3vGzKrM7KlslyeMdOoFM7vHzK7JQJkyEldSY5rHqvUws+XAPkAdsAt4C7jQObcyxLFtgWrgi865DzJZznSZ2VSgzDn3s2yXJWBmPwcOcM6dm+2yiOQbM5sJHAr0dc7tiNn+beCHwFjnXK2ZTQa+55z7clYK2khrrBda299I9qYWq9bnNOdcZ6AfUA78PuRx+wAlwPxkT2ge/b8gIpEzsyHA0YADTm/08WDgE+dcbUTnKo4ijkhanHN6tJIHsBwYH/P+ZLxKJ3jfHrgN+Awv6boH6AAcCGzFq7i2AC/7+48FZgFV/vPYmFgzgZuBN4HtwAFAV+BBYA2wCrgJKEpQ1vbAHcBq/3EH0N7/bBxQBvwUWO9/r2/5n03Ba43b6Zf1mcbfHfg58BTwCLAZ+ND/jlcB64CVwImNvsv3/NcHAK/633k98ETMfr/zj60G5gBH+9sn+OXZ5Zfpgzhx2wA/A1b4ZfgT0NX/bIj/tz/P/2+zHrg62/8/6aFHa3gA1/r1zG+B6THbr2/0u7sEqMFrsd8CbPL3i1vv+Z8Fdc3/AmuBP8c5f5jf7hS/HlsDXOF/FqZemOx/t9uBTcBSvHp3sl/XrAPOiynLVOAm/3UvYLp/3EbgdaCN/9mVwBK8+u9j4Cx/+/AEf6OGuP777wOL/bjTgH1jPnPAhcCnQCVwF37vlR7RPNRK0UqZWUdgIvBOzOZf4SUYo/ASiP7Atc65T4CR/j7dnHPHmVkP4F/A/wE98Sq1fzUae/VtvAqlFK/SeRio9WN/ATgR+F6CIl4NfNEvy6HAEXiVV6AvXsXRHy/huM/MDnLO3Qc8CvzaOdfZOXdagvinAX8GugPvA8/hVZD9gRuAexMcdyPwvH/cAPZs8Zvll7cH8BjwlJmVOOeeBX6Bl4R1ds4dGifuZP9xLLAf0Bm4s9E+XwYOAo4HrjWz4QnKKFJIvoP3m38U+C8z2wfAOXcde/7u7sL7B/9t/303//i49V5M/L54v+nBePVZY5Np/rd7LDAMr8670szGh6wXAI4E5uHVs48BjwNj/LKeC9xpZp3jHHcFXlLYG6/H4ad4SQ94SdXReBe71wOPmFk/59yCBH+jBmZ2HPBL4Ot4PR8r/DLFOtUv46H+fv+V4LtJCpRYtT7/MLNNeK0qJwC3gtddh3cVcrlzbqNzbjPej35SgjinAJ865/7snKt1zv0FWIiXsASmOufmO68ZvgdwEnCZc26rc24d3lVYovjfAm5wzq1zzlXg/fi/3Wifa5xzO5xzr+IleV9P4u/wunPuOb9sT+FVPrc453bhVRJDzGyvSgXv6nIw3hVajXPujeAD59wjzrkN/t/jN3hXwgeFLM+3gN8655Y657bgtZ5NatT1cL1zbrvzxrh9gFdpiRQsM/sy3u/xSefcHLyE4ZtJHB+m3qsHrvPrmu1xwoT97W51zn0I/BH4RhJfc5lz7o/OuTrgCWAgXt24wzn3PF6r1wFxjtuFl/gMds7tcs697pzfpOTcU8651c65eufcE3itS0eELM+3gIecc/9x3ni2q4Av+V2ygVucc5ucc58Br+AlrRIRJVatz5n+VUh74FLgVTPri5dYdATmmNkmP/l61t8ez754VyqxVuBd7QViB8UPBtoCa2Li3wv0CRl/hb8tUOmc29rE580pj3m9HVjvV1zBe/CuPBv7CWDAe/6dhucHH5jZFWa2wL8DaRPe1WCvkOWJ932L8a40A2tjXm9LUD6RQnIe8Lxzbr3//jF/W1hh6r0K51xNEzHC/HZXNvo8nboK51zjbfHqglvxuuueN7OlZnZl8IGZfcfM5sZ858+RYl3lJ5Mb2LPuV12VQRro10r5ScTfzOxevC6mv+H9QEc651aFCLEaL1mKNQivUmo4TczrlcAOoJcLN5A0iB8Mlh/kbwt0N7NOMcnVIOCjOOeNlHNuLd4VbnC1/KKZvYZ3Zfi/eN10851z9WZWiZeEhSlT47/nILxu03K8LkcRiWFmHfBaqYvMLPiHvD3QzcwOdfHvXm78O1xP8/VeFL/dgXgt+sHnQV2WybpqM1534BVmNhJ4xcxm4SVb9+PVVW875+rMbC4p1lVm1gmvmzLMvxsSAbVYtVL+nXpn4I0VWuCcq8f7sd1uZn38ffqbWaK+8RnAgWb2TTMrNrOJwAi8wZJ7cc6twRub9Bsz62JmbcxsfzP7SoL4fwF+Zma9zawX3piHRxrtc72ZtTOzo/H69IN5asrxxjpEzszOMbOgsqzEq4Tq8MaR1QIVQLGZXQt0iTm0HK97MdFv4i/A5WY21B8vEYy9iORuJpE8dCbeb28EXlfTKLzB16/jjbuKpxwYYGbtAFKo9+IJ89u9xsw6+gnOd/G69ILyNFUvpMzMTjWzA/zuzmq8v1Ud0Amv3qrw9/suXotVYI+/URyPAd81s1Fm1h7v+77rnFse9XeQ+JRYtT7PmNkWvB/azXh3lAStQv+LdzXzjplVAy+SYIyQc24DXjJzBV4z8E+AU2Oa5OP5DtAO7y6USuCveC098dwEzMYbtPkh8B9/W2CtH2M13qDVC51zwRXhg8AIv5n7H02UJxVjgHf9v+E04EfOuWV4g9//DXyC10xew57N/0HSt8HM/hMn7kN4g+lfA5b5x/8w4rKL5JPzgD865z5zzq0NHngDx7+VYGqEl/FawdeaWVBXha73Egjz233VP8dLwG3+2Chovl5IxzC877IFeBu42zk30zn3MfAbf1s5cAjenYeBeH+jBs65l4BrgKfx7nLcn8RjZSUDNEGoRM68GdQfcc6pi0xEWi1/QPcyoK1anyUqarESERERiYgSKxEREZGIqCtQREREJCJqsRIRERGJSNbmserVq5cbMmRItk4vyVi0yHs+KJkbcRKFWuSHSj+W5J45c+asd84lmtQ2Z6j+yi3b5s4FoOMoTTAuqQtbf2UtsRoyZAizZ8/O1uklGePGec8zZ0YQapwfKv1YknvMrPFqADlJ9VdumdvNW/1qlP6bSRrC1l/qChQRERGJiJa0keb97GcRhooulohIGPU//Wm2iyAFRImVNG/8+AhDRRdLJAwzmwD8DigCHnDO3dLo88l4C+IGa6nd6Zx7oEULKRl12E9+ku0iSAFRYiXN8wd+EsHAz7l+rFF5PIh0165dlJWVUVNTk+2iZE1JSQkDBgygbdu2WS2HmRUBdwEnAGXALDOb5i8bEusJ59ylLV5AaRGLnvCW/jto4sQslyQ3FHodlm79pcRKmnfZZd5zBAPOL/Nj5fPg9bKyMkpLSxkyZAje+qqFxTnHhg0bKCsrY+jQodkuzhHAYufcUgAzexw4A289TCkQ23/wA++FEqtQCrkOi6L+0uB1kYjV1NTQs2fPgquQAmZGz549W8vVbn/2XGy7zN/W2FfNbJ6Z/dXMBsYLZGZTzGy2mc2uqKjIRFlFWoVCrsOiqL+UWIlkQCFWSLFa0fePV5DGy008Awxxzn0eeBF4OF4g59x9zrnRzrnRvXvn/FRcIk1qRb/hFpfud1diJeHU1cH774OWQJLcUgbEtkANAFbH7uCc2+Cc2+G/vR84vIXKJiJ5SImVNG32bFiwAN58Ew47DPxBoNK6FRUVMWrUqIbHLbfc0vxBIc2dO5cZM2ZEFi/DZgHDzGyombUDJgHTYncws34xb08HFrRg+UQkjlyuwzR4XRJzDs44A9at273tmWdg0qSUQ/7iF7+IoGDSnA4dOjTcgRm1uXPnMnv2bE4++eSMxI+Sc67WzC4FnsObbuEh59x8M7sBmO2cmwb8t5mdDtQCG4HJWSuwZMQjZ1wEwO+nzuLByWOyXBoJI5frMLVYSWIrVsDq1dCjB0yf7m179dW0ugPHjh3L2LFjIyqgJKOqqoqDDjqoYb3Gb3zjG9x///0AXHTRRYwePZqRI0dy3XXXNRwza9Ysxo4dy6GHHsoRRxxBVVUV1157LU888QSjRo3iiRxowXTOzXDOHeic2985d7O/7Vo/qcI5d5VzbqRz7lDn3LHOuYXZLbFErfLYs6k89uxsF0PSlCt1mFqsJLFZs7znAw+Erl29BGvVKli+HFK8DfWtt94CKJzkKlMDQJtJbrdv377HXGFXXXUVEydO5M4772Ty5Mn86Ec/orKyku9///sA3HzzzfTo0YO6ujqOP/545s2bx8EHH8zEiRN54oknGDNmDNXV1XTs2JEbbriB2bNnc+edd2bmu4lErPsrfwNQcpUK1WFJU2IliQWJ1Zo13rI2Rx8N//yn12qVYmL1U39piXyex6o1SNSMfsIJJ/DUU09xySWX8MEHHzRsf/LJJ7nvvvuora1lzZo1fPzxx5gZ/fr1Y8wYr+ukS5cuLVZ+kSid/ZT3D+iDSqxyRi7XYeoKlMTee897Li31nr/yFe/5tdeyU55c5FxmHimqr69nwYIFdOjQgY0bNwKwbNkybrvtNl566SXmzZvHKaecQk1NDc65gr7lWkRQHZYCJVYSX10dzJnjvQ4Sq2OO8Z6VWOWs22+/neHDh/OXv/yF888/n127dlFdXU2nTp3o2rUr5eXl/Pvf/wbg4IMPZvXq1czyWy43b95MbW0tpaWlbN68OZtfQ0QKVC7UYeoKlPgWLoQtW2DQIGjXztt26KFekrVkiTfWqn+8CaylNWg8PmHChAmcf/75PPDAA7z33nuUlpZyzDHHcNNNN3H99dfzhS98gZEjR7Lffvtx1FFHAdCuXTueeOIJfvjDH7J9+3Y6dOjAiy++yLHHHsstt9zCqFGjGsY9iIhEKZfrMCVWEl8wvuqIIyBYvqO4GI46Cp59Fl5/Pa1pFySz6urq4m5fsGD3FE2//e1vG15PnTo17v5jxozhnXfe2Wt7cAUoIpIJuVyHKbGS+IL/6caMgRNP3L39mGO8xOq111JKrO64446ICigiEs5T512V7SJIAVFiJfEFA9ePOAJimmMJpkkIxl8lKbZpV0SkJVQfeUK2iyAFRIPXZW87dsAHH3jzlxx+OLz4ovcAOOAA73n58pRCv/jii7wYxBIRaQF9ZvyJPjP+lO1iSIFQi5Xs7aOPYNcuGD7cG6x+003e9vHjoV8/aNvWW+Zm2zbo2DGp0Df5scaPHx91qUVE4jpt+lQAHjz5O9ktiBQEtVjJ3pYu9Z4PPnjvz9q0gcGDvdcrVrRcmURERHKAEivZ22efec+DBsX/fMgQ7znF7kAREZF8pa5A2ZsSq0hdMDXa23ofnDwm1H7l5eVcfvnlvPPOO3Tv3p127drxk5/8hLPOOivhMSeffDKPPfYY3bp122P7z3/+czp37syPf/zjtMouIrknG3VYLtdfarGSvQVdfEqscpZzjjPPPJNjjjmGpUuXMmfOHB5//HHKysqaPG7GjBl7VUoiIi0p1+svJVayt8YtVvfe6z0CaSRW9957L/fGxpKMePnll2nXrh0XXnhhw7bBgwfzwx/+kKlTp3LppZc2bD/11FMbFsUeMmQI69evB7zV4g866CDGjx/PokWLGvZfsmQJEyZM4PDDD+foo49m4cKFLfOlRFL0+IU38fiFN2W7GBJSrtdf6gqUvTVOrA46aM/Phw71nlNIrA5qHEsyYv78+Rx22GEpHx9cIb7//vvU1tZy2GGHcfjhhwMwZcoU7rnnHoYNG8a7777LxRdfzMsvvxxV0UUit/XzY7NdBElCrtdfoRIrM5sA/A4oAh5wzt2SYL+vAU8BY5xzsyMrpbScrVthwwZvfcA+fbxtzzzjPZ92mvecRovVM36s04JY0iIuueQS3njjDdq1a8cll1zS7P6vv/46Z511Fh396TROP/10ALZs2cJbb73FOeec07Dvjh07MlNokXTV1cE//8mPrv4xB6xbxevHnwOTH8l2qSRJuVZ/NZtYmVkRcBdwAlAGzDKzac65jxvtVwr8N/Bu5KWUlrNypfc8cKA3tQLAb37jPQfJUN++XuK1bp2XiHXqFDr8b/xYSqwya+TIkTz99NMN7++66y7Wr1/P6NGjKS4upr6+vuGzmpqauDHMbK9t9fX1dOvWjblz50ZfaJGo3XorXHUVn/ffHvPqP7w5+tq2zWqxpGm5Xn+FGWN1BLDYObfUObcTeBw4I85+NwK/BuJ/S8kNzd0RCJrLKgccd9xx1NTU8Ic//KFh27Zt2wBvHMLcuXOpr69n5cqVvBcsXxTjmGOO4e9//zvbt29n8+bNDS2NXbp0YejQoTz11FOAN8j0gw8+aIFvJJKC558HYHVxW7ab0aFm6+7luqTVyvX6K0xXYH9gZcz7MuDI2B3M7AvAQOfcdDNLeD+jmU0BpgAMauofbsmeMIkVeN2Bn37qdQeOGJHpUuW0sNMjRMnM+Mc//sHll1/Or3/9a3r37k2nTp341a9+xVFHHcXQoUM55JBD+NznPhd3LMNhhx3GxIkTGTVqFIMHD+boo49u+OzRRx/loosu4qabbmLXrl1MmjSJQw89tCW/nkjzdu2Cd70OlFXF7TAHHep2wQsvwFFHZblwuaWl67Bcr7/CJFZ7t6eBa/jQrA1wOzC5uUDOufuA+wBGjx7tmtldsiGZxAo05UIr1q9fPx5//PG4nz366KNxty+P+e959dVXc/XVV++1z9ChQ3n22WcjKaNIxnzwgbfs1rBh7Cpbw6aiIvrV7eLTx/7BLUNOycoFj4SXy/VXmK7AMmBgzPsBwOqY96XA54CZZrYc+CIwzcxGR1VIaUFKrEQkH7z5pvfst05Vtymirk0R+y35iJLtW7JYMMl3YVqsZgHDzGwosAqYBHwz+NA5VwX0Ct6b2Uzgx7orMEfFS6z+/Oe990sxsfpzvFgiIlELEquxY3lkwDEA9Hj8DoYtnsdBi94Hjs1e2SSvNZtYOedqzexS4Dm86RYecs7NN7MbgNnOuWmZLqS0oHiJ1cCBe++XYmI1MF6sPOSci3tXSqFwTj39kkXO7dFitaPtVgA+HnkEwxbPY8R83bzenEKuw9Ktv0LNY+WcmwHMaLTt2gT7jkurRJI99fW7p1uITayeeMJ7njhx97YUE6sn/FgTY2PlmZKSEjZs2EDPnj0LsmJyzrFhwwZKSkqyXRQpVCtWwOrV0L07HHwwg3/kdbJ8POZ4zvjnA4z4ONq17/JNIddhUdRfmnlddisvh507oVcv8CdWAyC45TU2GQrmsqqoSGouq+D22XxOrAYMGEBZWRkVFRXZLkrWlJSUMGDAgGwXQwrVW295z2PHQps2jH/zXwBMnfgjako6su/qZV59t88+WSxk61XodVi69ZcSK9kt7MB18Oay6t8fli2DNWvggAMyW7Yc0rZtW4YGy/6ISMtrNHA9UFdczKr++7H/ko9g0SIlVgmoDkuPFmGW3ZJJrMBrtQJYuzYz5RERSUXQYhVnvqp1vf2WiCVLWrBAUkiUWMluSqxEJNc5B5984r2OM/FjRR8/sVq8uAULJYVEiZXsFrtOYBhKrESktVm3zpsYtHt36Np174/79PdeqMVKMkRjrGS3IEHq12/P7X/9a/z9U0is/poolohIBG6+819cDRAzRujha/7Y8HpdH//CUYmVZIgSK9mtvNx7bjygs1evvfeN3S+JxKpXolgiIhHovW6V92K//Rq21fYb3PC6ImixUlegZIi6AmW3IEFqnFhNneo9GkuhxWrq1KlMjRdLRCQCvdb7K67FtFjt/9BN7P/QTQBUd+lBTfsOsGkTbNyYjSJKnlNiJbslarFqLrEKjgtBiZW0NDObYGaLzGyxmV3ZxH5fMzOndU5zW6+KvROrcbNfZtzsl703ZqzrozsDJXOUWIln506orPTmp+rZM9wxGrwurZyZFQF3AScBI4BvmNmIOPuVAv8NaK2THNc7TotVY7ozUDJJiZV41q3znvv0gaKicMcELVvl5d5yOCKtzxHAYufcUufcTuBx4Iw4+90I/BqoacnCSfTidQU2pjsDJZOUWIknUTdgU0pKoFs32LXLa+0SaX36Aytj3pf52xqY2ReAgc656U0FMrMpZjbbzGYX6lIfrV5tLT02+HXZ4MEJd9OdgZJJSqzEk0piBeoOlNYu3gqyDUvXm1kb4HbgiuYCOefuc86Nds6N7t27d4RFlMisXElRfR2V3Xp7F34J6M5AySRNtyCeRHcEAsyYkfi4vn1h4ULv+JEjmz3NjKZiiUSvDIid8XYAsDrmfSnwOWCmmQH0BaaZ2enOudktVkqJxrJlAFT03pfuMZv/+Ku/7bGbBq9LJimxEk9TLVYdOyY+Lsm5rDo2FUskerOAYWY2FFgFTAK+GXzonKsCGiZXM7OZwI+VVOUoP7Fa33tfDozZXF/abY/dNvbYB9q29RaQ37oVOnVqwUJKvlNXoHiCxCro2ot1993eI54kuwLvvvtu7k4USyRizrla4FLgOWAB8KRzbr6Z3WBmp2e3dBK5ILHqtccwOob/4WqG/+HqhveuTdHuwe1Ll7ZY8aQwqMVKPE21WD35pPd88cV7f5bkXFZP+rEujhdLJAOcczOAGY22XZtg33EtUSbJkJgWq1hjP3wb8DLrBvvv7y3WvGQJHHJICxVQCoFarMSjwesikuuCMVa99m1mR7zECjTOSiKnxEo8SqxEJNf53XqNW6ziGuAPYF+9uun9RJKkxEo8SqxEJJdt2wbl5dQWFVPZPcR0GP39cVhKrCRiSqzEm+BzwwZvOZtevZrfP5YSKxFpDZYvB2BDz37e4PTm7Ou3aimxkohp8LrsXs6md+/4y9nMnJn42F69wAzWr/cStLZtmzzVzKZiiYikatUqADb26LPXRw/e/fLe+weJlX+cSFSUWEnq3YAAxcVeQrZuHVRU7K6sRERakt/ytKmb1w14wdRZTe8f2xXonHeBKBIBdQVK84nVbbd5j0SS6A687bbbuK2pWCIiqVizBoCqbnsPZ/j87f/D52//nz03lpZC586wfTtUVbVECaVAKLGS5hOr6dO9RyJJzGU1ffp0pjcVS0QkFQ0tVnsnVmM+ncuYT+fufYy6AyUDlFhJel2BsDux8q8YRURaXBMtVvFcMHUWC9qUem80gF0ipMRKdnfhxVvOJow+/mDRYBC8iEhLa6LFKpFgPJYSK4mSEitJv8Wqt185VVREUx4RkWQFiVXXJBKr7r32OFYkCrorUJpPrDp0aPr4JBKrDs3FEhFJlnNNdgXuKIr/T92mbn5ru8ZYSYSUWEnzidW//9308UFXYIjE6t/NxRIRSVZlJezYAV26sKOk414fP/L75+MephYryQR1BcrusVHqChSRXBTcONOvX1KHaYyVZIISq0JXV+fNmg7Qs2f8fW680XskkkRideONN3JjU7FERJIVJEYJJige/euLGf3ri/fa3rCmoBIriZASq0K3caM3PqF798TL0bz0kvdIJEis1q3zYjXhpZde4qWmYomIJCtIjBK0WB26fCGHLl+41/aqYKD7mjVQX5+p0kmBUWJV6IJWpt4hVoNPpFMnKCmBmhrYujWacomIhBV0BSa5pFZt23Zs7twVams1lEEio8Sq0EWRWJklNYBdRCRSzbRYNWWTugMlYkqsCl0UiVXs8UqsRKSlpdhiBTED2DXlgkRE0y0UujCJVaJB7bFCJlY9w8QSEUlG7OD1pXt/XNU+8fx5DTO1q8VKIhIqsTKzCcDvgCLgAefcLY0+vxC4BKgDtgBTnHMfR1xWyYQwidXTTzcfJ3YAexOeDhNLRCQZsdMtLN2018d/vf1fCQ+t7O4PY1BiJRFptivQzIqAu4CTgBHAN8xsRKPdHnPOHeKcGwX8Gvht5CWVzFBXoIjkMufSG2Plt1i9OvMDLpg6K8qSSYEKM8bqCGCxc26pc24n8DhwRuwOzrnqmLedgKbvuZfWI0xiddVV3qMpIROrq666iquaiyUiElYw63ppKXTuHHeXsTeez9gbz4/7WTDGqnulLgolGmG6AvsDK2PelwFHNt7JzC4B/gdoBxwXSekk88IkVm+/3XyckHcFvh0mlohIWCEGrg9fsxyAt+J8FrRYdd20PtpyScEK02Jlcbbt1SLlnLvLObc/8L/Az+IGMptiZrPNbHaFuoxaB3UFikguS6MbEKC6i3dDTZfqjVGVSApcmMSqDBgY834A0NQov8eBM+N94Jy7zzk32jk3une6/5BLNJRYiUguS2OqBYDNXboDXmJlmn1dIhAmsZoFDDOzoWbWDpgETIvdwcyGxbw9Bfg0uiJKxji3e53AqBKrZu4KFBGJVDPrBDantm07tnbqQlF9HZ22Vjd/gEgzmh1j5ZyrNbNLgefwplt4yDk338xuAGY756YBl5rZeGAXUAmcl8lCS0Q2bfKWcigthfbtE+83YEDzsUK2WA0IE0tEJCy/xeqJlbt4PsFdfRWdSpsMUdWlB522VtOlekPkxZPCE2oeK+fcDGBGo23Xxrz+UcTlkpYQthvwkUeajxUkZ9u2eY+OHROEChFLRCSs8nIAqromnnx42q3/bDJEdZce7LtmOV2rlFhJ+rSkTSELuu2iGO9mpnFWItLy/MSqukuPlENU+0lZlyoNYJf0KbEqZGFbrC67zHs0J0Riddlll3FZmFgiImGESKzGXXcu4647N+HnwbFd1RUoEdBagYUsbGI1d264eCESq7lhY4mIhBEkVl0TJ1b7V3gD3Gcm+LxKLVYSIbVYFbIgAQom90yX7gwUkZa0cyds3EiARfhAAAAgAElEQVS9tWFL564phwlarDSXlURBiVUhi2oOq4DGWEkrZGYTzGyRmS02syvjfH6hmX1oZnPN7I04a6FKa+VfxFV36Y5rU5RymKC1q4sGr0sElFgVsqgTq5DL2oi0FC0in+ciGLgOu7sCu6rFSiKgMVaFLGxideCB4eKFaLE6MGwskWg0LCIPYGbBIvIfBztoEfkcFjKxWtkt8VQMscerK1CioMSqkIVNrO67L1y8EInVfWFjiURDi8jnsxAD1wGe/8VTTX4eJFal1ZVQXw9t1JkjqdP/PYVMY6wk/2kR+XzW0GLVdItUc+qK27LFX9aGDRpnJelRYlWonAufWE2Z4j2aE+KuwClTpjAlTCyRaGgR+Xy2di3QfFfgiT89hxN/ek6T+wSThAYxRVKlrsBCtXmzd6tyx44Jl59p8Mkn4WKGaLH6JGwskWg0LCIPrMJbRP6bsTuY2TDnXLBwvBaRzyUhlrMBGLip+Vao6i492Hf1Mi/mIYdEUjwpTEqsClXU3YAAXbtC27awZQvU1EBJSXSxRVKgReTzXMgxVmE0JGd+TJFUKbEqVJlIrIL1Alev9uIPHNj8MSIZpkXk81hE0y3sEUNdgZImjbEqVJlIrGLjaXCviGRaJhIrtVhJmtRiVaiCAeZhlrMZNSp83GYSq1HJxBIRSWTXLu8OvjZt2Fzarcldl/Tet9lwVRq8LhFRYlWokmmxuuOO8HGbuTPwjmRiiYgkEtQxvXs3u5zNzOsfaTacWqwkKuoKLFQxlVKk1BUoIi0hSID22SeScA0D4JVYSZqUWBWqIPEJ0xV47rneI4xm1gs899xzOTdsLBGRRJJIrE7/f2dw+v87o8l91BUoUVFXYKFKpiuwrCx83GZarMqSiSUikkgSiVXvrZub3Wdzqd9iVVEBdXVQ1HT3okgiarEqVMkMXk+GugJFpCVE3BVYV1zM5s5dvbUC16+PJKYUJiVWhUrTLYhILos4sQINYJdoKLEqRMmsE5isEOsFioikLRgL1bdvZCGVWEkUNMaqEG3eDDt2eGsEdurU/P5f+lL42M20WH0pmVgiIonEtlg1taw2sKDfkFAhtayNREGJVSFKtrXql78MH7tbNyguhupqL3lr375RqCRiiYgkskditbPJXd+65qFQIRumXNCdgZIGdQUWomSmWkhWmzbQq5f3WgNARSRTgsQq0q5AtVhJ+pRYFaJkJwf96le9R1hNdAd+9atf5avJxBIRaWzXLu/CLfZCrglfu/wUvnb5Kc3uV62uQImAugILUbJdgRs2JBe/icRqQ7KxREQai63DQsw31XXH9lBhq7qoK1DSp8SqEGWyKxBC3Rl4wdRZe217cPKYzJRHRPJLkPhEONUC6K5AiYa6AgtRptYJDGguKxHJpAyMrwJ1BUo0lFgVokzNYRVQYiUimZShFqvNpd29F8GyNiIpUFdgIUp2OZvjj08ufhMLMR/vx1qeXEQRkd2SnBz0gyEHh9qvrrgYevb0xpWuXx954iaFQYlVIUq2xeqaa5KL30SL1TV+rHhjrEREQkmyK3D2T+4OH7tvXy+xWrtWiZWkRF2BhailBq+rK1BEMsFvsbpv0dboL9KCZErjrCRFSqwKjXPJD14/6STvEVYTdwWedNJJnJRMLBGRxvykp2Gm9Gac+8MTOfeHJ4aLrcRK0qTEqtBs3gw7d3rrBHbsGO6Y7du9R1hNtFht376d7cnEEhFpzG+xapgpvRnt62ppX1cbat/ng6n2NJeVpEiJVaHJdDcgQI8e3ozImzZ5MySLiETJb01qWDQ5QprLStKlxKrQZHoOK9hzmQmNsxKRKO3YAZWV1BYVsbVTl8jDay4rSZcSq0LTEi1WoHEKIpIZ/sXh5tIeuDbR/xPW0AqmrkBJUajpFsxsAvA7oAh4wDl3S6PP/wf4HlALVADnO+dWRFxWiUIqk4Oeemry5+nbFz78cK/E6lQ/1oLkI4qINCQ8yXQDzho2KvS+6gqUdDWbWJlZEXAXcAJQBswys2nOuY9jdnsfGO2c22ZmFwG/BiZmosCSplS6An/84+TPE7RYNbrq+7EfS/NYiUhKGgauh7sjEGDe5b8Nva8SK0lXmHbUI4DFzrmlzrmdwOPAGbE7OOdecc5t89++AwyItpgSmUwvZxNQV6CIZELDVAvRD1wHb1mbejOvrqwNdyehSKwwiVV/YGXM+zJ/WyIXAP+O94GZTTGz2WY2u0KDmrMjyaUgABg3znskI0FiNW7cOMYlG0tEJBB0BSbRYnXBxcdxwcXHhdq3rriYrZ26enP+rV+fUhGlsIVJrCzONhd3R7NzgdHArfE+d87d55wb7Zwb3TvTLSYSXyqJVSqC+GqxEpEoZbjFyout7kBJXZjEqgwYGPN+ALC68U5mNh64GjjdObcjmuJJ5JJcYytZF0ydxQVTZ/HbDzYBsOD9TzJyHhEpUCkMXk+W7gyUdIRJrGYBw8xsqJm1AyYB02J3MLMvAPfiJVV7r2MirUcLtVgFzfRdqjdm9DwiUmCSXM4mFRrALuloNrFyztUClwLP4d0l/6Rzbr6Z3WBmp/u73Qp0Bp4ys7lmNi1BOMmmHTtg40YoKoKembvag92VXpcqJVYiEiG1WEkrF2oeK+fcDGBGo23XxrweH3G5JBOCqRb69PFmRw/r619P+lTenTVt6LS1ylvWpm1bP5QXa07SEUVSo3n48kzQYpXE4PW3DvlSUqeo6uqvHLFmTVLHiUDIxEryRKrjqy6+OOlTuTZFbC7tRtfqjd5ty/vu64fyYmkeK2kJmocvz2zfDtXV0K4d2zqWhj5swUU3J3WaTd38xGr1XsOJRZqlxKqQpDq+aps/RVnHjkkdVt2lB12rN3L9/S/z2eCDAKjdUQNAcfuS5MogkpqGefgAzCyYh68hsXLOvRKz/zvAuS1aQgkvuDjcZx+weDesx9dms3czTX1pt1D7b+ru3bX+yfuL+NXUWTw4eUxy5ZSCpsSqkASJVTDHVFgnn+w9z5yZ1GHVcQawP3/7ZV7IK+9JrgwiqYk3D9+RTeyfcB4+aQWCrrkkLw6/+79nA/Dg3S+H2n+T3xXYtWpDUucRASVWhaWl5rDyBfPMdFHlJNmTyjx8X0nw+RRgCsCgQYOiKp8kI+ia69/UHNXpq/K7ArttqvAmChVJQvRLg0vr1dKJlaZckOyLbB4+TXDcCgSJlT9mM1NqOnSipn0H2u/cQYftWzJ6Lsk/SqwKSYYnB22syp9yQc3pkkWahy+frFrlPWc4sYLdA9i7bdKyNpIcJVaFJNUxVilSi5Vkm+bhyzMt1GIFsKmb1yrZbZPWtZXkaIxVIUm1K3Dy5Libm5syoWGMVUxiNezLpyZ3bpE0aR6+PBI7xiqJmRBmjg63AHOsYJxV101qcZfkKLEqJBEnVs1pWNamSomViEQgtsVq9fbQhy05/2dJn2pT7AB2kSQosSoUW7fCli3Qvj107Zrcsev9MQa9eiV1WHWcMVY1/nwyJXHmk4nXAqb5Y0SkwR5dgUtCH1a8xptIv7bf4NDH7O4K1BgrSY4Sq0KR4sR6AHzta95zkvNYbencjXozOm2tok1dLfVFxbx815WA5rESkSRt3QpVVd7FYffuSR163o3fBcLPYwW71wvsqsRKkqTEqlCE7AaM22qU4inri4rZXNqdrtUbKd28qWHMgohI0mJbq5K9OExBMPu6ugIlWborsFC08BxWAd0ZKCKRaME7AkHTLUjqlFgVihaewyoQJFaay0pE0tJCs64HqhqWtVmv2dclKeoKLBQtPIdVIBPL2iSa5kED3UXyWAtODgr+7OslHSmp2eaN7eoWbgFnESVWhSKdrsCLLkr5tJu6eYlV0Jx+8LFfTTmWiBSwNLoCXzzqlJROualrT/rWbPMWf1ZiJSEpsSoU6SRWEyemfNrK7l4LWfdKb6WQ/Y48IeVYIlLA/MTqvsU1vNvM5MSNrfjWFSmdclO33vQtX+mde/jwlGJI4dEYq0KRTmK1cqX3SEFwZ02QWG3ZUM6WDeUpxRKRAuYnVptSuLu4/ZKPaL/ko6SPazjX6iSmeZeCpxarQpHOGKtvf9t7TnIeK4CN3fsAuxOr1+6/DtA8ViKSpCCx8i/WknHub/4bSG4eK9i9rA1r1iR9TilcarEqBHV1LX5HTaCyh59YbVzXoucVkTziXMPg9eBuvZYQzL6uFitJhhKrQlBeDrW13pI0JSUteurqLj2otzaUbq6kqHZXi55bRPLEpk1QUwOlpdR06NRyp1VXoKRAiVUhKCvzngcObPFT1xcVs6lbL9o4p4n2RCQ1LTw5aKChxSqY6kEkBCVWhSAYeD5gQFZO33gAu4hIUrKUWG3s4Y9JTfHmHSlMGrxeCNJtsboitVuVA5Xd+wDz6Va5js9N+FZSxyaaDFRECseDT7zBBcDb29undPy/x389peMqe/Sh3ow2q1Z5wymK9U+mNE//lxSC4Gor1cTqtNPSOn2lf2dgj8p1DPqv5BIrEZFgIeRU7ggEWH32hSkdV1fclqquvei+qcLrDhw8OKU4UljUFVgI0u0KXLTIe6SoMugK3LiOqjUrqFqzIuVYIlJ4em7wposJLtKS1WneW3Sa91ZKx27s6XcHfvZZSsdL4VFiVQjS7Qr8wQ+8R4qCyrBbZQVvPvxL3nz4lynHEpHCEyRW63unNsZq0j0/Y9I9P0vp2A09+3kvVuiCUMJRYlUI0u0KTFNlo0lCRUSS0XO9N0Hnhp4prByRpg1qsZIkaYxVvotgctCFazcDcGuKA8mDSUJ7VK6DXi17V4+I5Djn6LkhSKz6tfjp1WIlyVKLVb5bu9ZLrvr0gfap3VGTrmCMVddNFZhzWSmDiOSoigra79zB1o6lbO/YucVPvzFoJVOLlYSkFqt8F4yvijNwvaWmMqht257NnbtSuqWKXnW1VBS3bZHzikge8FuKNvRq+dYqiOl+VIuVhKTEKt9FML5q+mnnp12Myu59KN1SxVe+NIFP+2bmluV4ieKDk8dk5Fwi0kKWLwdgfRqJ1TOnTk752A2xLVbOgVnKsaQwKLHKdxEkVgtGHpF2MTZ178OglZ9yYPc+bIsgnogUCD+xSmd81bqTv5Pysds7lkKXLlBdDZWV0KNHyrGkMGiMVb5roiswrIGffcLAzz5Jqxgb/TsD2y35iA1pxhKRAuJ3waXTYtXl3Rfo8u4LqZdh0KA9yiLSFCVW+S6CFqtJj/2WSY/9Nq1iBHcG7nznWd5NM5aIFJCGFqvUp1o45+Ffck468+cFM65rALuEoMQq30XQYhWFYJX4vrW7sloOEckxDYPXszhVi1qsJAkaY5Xvsjw5aCBoseq/a2dWyyEiueGCqbPAOe5avJQSsjM5aIMgsVKLlYQQqsXKzCaY2SIzW2xmV8b5/Bgz+4+Z1ZrZ16IvpqSktjbtyUGjsq6P12I2eNeOrJZDRHJHp63VlNRsY3tJJ7Z26pK9ggRdgWqxkhCaTazMrAi4CzgJGAF8w8xGNNrtM2Ay8FjUBZQ0rF0L9fWwzz7Qrl1Wi7KxR1/q2hTRr3YX7evrs1oWEckNwVI263v1y+40B2qxkiSE6Qo8AljsnFsKYGaPA2cAHwc7OOeW+5/pX8zWJKgE0uwG/NtXL067KHXFxWzo2Zc+FauYcOzZaKSViDSnV0RrBP7tnEvTK4harCQJYboC+wMrY96X+duSZmZTzGy2mc2uqKhIJYQk49NPvecDDkgrzJJhn2fJsM+nXZygO3BEh05pxxKR/NewRmCas65XHns2lceenfLx33u+jNqiIigvh5qatMoi+S9MYhWv/TWlBd+cc/c550Y750b37t07lRCSjCCxGjYsrTD7fzqP/T+dl3ZxKvzEqmT+e2nHEpH8t7vFKr3Eqvsrf6P7K39L+XjXpoiNPfxWs2XL0iqL5L8wXYFlQGxf0gBgdWaKI5GKKLE6++m7Abj1ynvSihO0WBXPfhm+dUVascLSMjdiZhOA3wFFwAPOuVsafX4McAfweWCSc+6vLV9KiSdosUpnclCAs5+6E4AH02i1Wtt3MH0qVsGiRTB8eFrlkfwWpsVqFjDMzIaaWTtgEjAts8WSSESUWEUlSKyG6M5AaSG6+Sa39YyoxSoKa/v5A9gXLcpuQaTVazaxcs7VApcCzwELgCedc/PN7AYzOx3AzMaYWRlwDnCvmc3PZKElBOdaXWIVdAUO3qnESlpMw803zrmdQHDzTQPn3HLn3DxAN9+0IlZfzz7l3g046/bJ7gTHAOX7KLGScELNY+Wcm+GcO9A5t79z7mZ/27XOuWn+61nOuQHOuU7OuZ7OuZGZLLSEUF4OW7ZA9+7Qs2e2SwNARW/vnodBu3bQpq42y6WRAqGbb3JUj43ltN+5g6ouPdiWzTmsfGv7encGfvr6nLhDDEQCWtImX7Wy1iqAne1LWFPclrZ4laZIC9DNNzmq35rlAKzZd2hWyxFY289LrPqu1ZQL0jQtaZOvIkysHv/m/6QdI7Cx/370W7GI3utWsb53dmeDl4Kgm29yVL/V3t13a/oNSTvWU+ddlXaMTd16U1PSkdItVXTasinteJK/1GKVryJMrFYOOpCVgw5MOw5A5UCvPH3WrWxmT5FI6OabHBW0WK3ed0jasaqPPIHqI09IL4gZa/1xVn3XaAZ2SUyJVb6KMLEaPv89hkc099Ri/7nPulWRxBNpim6+yV1Bi9XaCFqs+sz4E31m/CntOGv7+omVugOlCeoKzFcRJlanPvMQAAtGHpF2rHeWfMRkoPe6srRjiYThnJsBzGi07dqY17PwugilFem3xkteougKPG36VAAePPk7acUp76txVtI8tVjlI+dgsd821IoGrwMsb9segD5KrEQkkfXrKd2yie0lnajs3ifbpWmgAewShlqs8tHq1bBtG/TqBd26AfFnIM+G5e28xKp3xSovAczmivUi0jotWADAmn6DW1UdsbahxUpjrCQxtVjlo1Y41UKgqqiYjUVFlOzYrikXRCQ+P7EKWohai2CMVZ/ylVCrufgkPiVW+agVJ1YA89t3BGDgyk+zXBIRaZUWLgRgTb/WMYdVYGf7Dmzs3ofiulpYvjzbxZFWSl2B+civlKJKrP4UwRwwgaPOu4rKGQ/DG9MZ+NknfDDq6Mhii0ieCLoCI5hqAeDxC2+KJA54rWg9Ktd5S9sccEBkcSV/qMUqH/3nP97zqFGRhCvvN5jyiJrku/YbTPnwMQAMWqE1t0QkDj+xWh1Ri9XWz49l6+fHRhJrTd8h3ov5mpVD4lNilW/q63cnVocfHknIQ+e+zqFzX48k1mdzX+e9bZsBGLTyk0hiikge2boVVqygtqg4stUZ9v3bPez7t3siibViyMHei1mt44YgaX3UFZhvliyB6mro1897RODEZx8FiKTb7qNnH2Whc+wqbkfvitV02LaZ7R1L046bjHh3SD44eUyLlkFEEvj4YwDK9xlIXXE0/0Sd9OKTADx49oVpx1o+dIT3QomVJKDEKt/MmQPA3D778/tWMsVCY7VmrBqwH0OWL2Tgyk/55KDDsl0kEWkt3nkHiElgWpnV+w6hpn0HSlasgIoK0ILc0oi6AvONn1g1NFe3Up8NOgiAQSvUHSgiMfzEasn+n8tyQeJzbYrUHShNUmKVbxoSq+FZLkjTPvMXdR6ocVYiEqshsTokywVJbPkQrzXtn/f9o9VMviythxKrfOJcw8D1FYNbd4vVSj+xGvSZEisR8a1bB0uXQqdOrO6/X7ZLk9Ayv5ty6LKPs1wSaY00xiqfLFkCVVVs6taLTd2j6/d/4PvXRxbrGD9WWcfO1Jux76qlFNXuoq64bWTnEJEc5bdWMWYM9UXR/fP0yBX/F1ks2J1YDVn2sXdBKxJDiVU+8bsBl0fcWlXZc5/IYnX2Y9UA6/oMoG/5SvZdvayhBUtECliQWH3pS5GG3RHxeK31vfdlc+eudNlcSc8NayONLblPXYH5JEMD18e8+wJj3n0hklhL332BpX6sIJlSc7qIAPD22wD839aekYYd/OhvGPzob6ILaNZw1+IQ1V/SiBKrfJKhgevjXnmaca88HUmsha88zUI/1icHfgGA4R9r8KdIwautbbjLbtl+IyMNPf7NfzH+zX9FGnP3OCvNwC57UldgvqipaWhGXza0dd8RGPjoc18EYMT8d7H6OlyboqyVRZOGimTZ/PmwdSsVvfelumu0LVaZECRWByz+MMslkdZGLVb5YuZM2LYNvvAFqrv2ynZpQlm3z0AqevWj89ZqBi/XuoEiBc3vBmzN0yzE+uSgL1BbVMz+iz+EDRuyXRxpRZRY5Yt/+c3cp5yS3XIkw4z5fqvV5z56O8uFEZGs8uuwYIhAa1fToTOLDj6cNq4eZszIdnGkFVFilQ+cy83ECpg/0kusRs5/N8slEZFsueQPr7Dr2eeoN+P9w47JdnFCmxusnzptWnYLIq2Kxljlg4ULYdky6NULxoyBhf+JNPwfLrklsljHNYq1YMRo6toUsd+SDynZvoWaDp0jO5eI5IbPf/AmbWt38cmBozIylOHha/4YeUzwFqb/1qO3wbPPwo4d0L59Rs4juUUtVvlg+nTv+aSToCj6AeBbSruxpbRbJLFKSrtREhNre8dSlu43kuK6Og5eMCeSc4hIbjl8zisAzDn8uIzEr+03mNp+gyOPu6FXP1YOPAC2bIFXX408vuQmJVb5IOgGPPXUjIQ/6o3pHPXG9EhiffrGdD5tFOvjkUcCGmclUpC2beOQeW8BMOfwYzNyiv0fuon9H7opI7HnjvK7LtUdKD4lVrmushLeeMNrqTrxxIycYuwb0xmbwcQqGKdwxHsv0m5HTSTnicIFU2ft9RCRiD37LO131rB0v5GRrvIQa9zslxk3++WMxP4gdpyVlrcRNMYqp8Sda2n9q1BXB+PHQ7douuta2mdDDmbpfiPZb+l8vvjOs7z2lTOzXSQRaSlPexMGZ6q1KtOWDxkO/fvDypXw3HMwYUK2iyRZpharHFZUuwvuuMN7c/nl2S1Mml4cPxGA4158Uld9IoVi5Ur461+pN2PO6MyMr8o016YNT409C4CFl12d5dJIa6AWqxx25DvPwapVlPXfn+vKe0EOd1XNHnM8Ex+/g4Flizlw0ft8cvBh2S5SXJqhXSRCv/gF7NzJrCNPpKLPgGyXJmUzjz2bU6ZP5eBF//EmOo14EWnJLWqxylFWX8+EZx8B4NmTzgWzLJcoPXXFbXl13NkAHP/Sk1kujYhk3PLl8OCD0KYN0874XrZLk5aaDp15+fiveW9u8aaU0djMwqUWqxx16NzX6b9qKRu79+G9IzMzaD3wu8vviCzWiU3Emnns2Zz8rz9y2JyZHPDpBywedmhk5xWRVubmm2HXLjj3XNb2G5LRU/3xV3/LaHyAF0+YxInP/YV206bBvHkZP5+0XmqxykFdqtbznYd/CcDz//VN6orbZvR8O9uXsLN9SSSxituXUJwgVlW3Xjz3X+fSxtXzvfuvo2T7lkjOKSKtzNtvw9Sp3t3M116b8dPVl3ajPqK5+BLZ3KUHr33lDO/NN79Jux3bM3o+ab3UYpVj2tTV8oN7rqFr9UYWHnw4L43/esbPeezLfwXgleO+lnasBX6s4Qli/fOsKXzuo3cY/NkivvHY7fzxgmvSPmemJWrq19grEU/sb+TBcb3hzDOhthZ+9CMYNgzezGx32fA/eIPKF1x0c0bP8/ezL2J82TyYP5/z/ngz9//gxpwfpiHJU4tVDrH6OiY+fgcHL5zDpq49uffCG6kvynxuPPq9Fxn93ouRxFr23ossayJWXXFb7p9yPTvbtufLbzzDqdMe1F2CInmiw7YtcNppsG4dnHAC3Hpri5x37IdvM/bDzE9AXNOhkzd9RKdOfPHd5zn1mYdUfxUgtVjlitWrueLWSxm+cA51bYq498KbM7KmVmuwpv9+PHru/+O8qTdz1t/vZUDZYv54/jXsKOmY7aIlRXcQiuw2/OP3+O5DN8GGtXDwwVx65pVsf3RutosVvREjvEH5kyZx1t/vZfCKhfz31p+xtXNXQHVAIQiVWJnZBOB3QBHwgHPulkaftwf+BBwObAAmOueWR1vUwhH7D3LfNcs55tV/8OU3pjN8azVVXXpw/5QbWu10BFF545jTqe7SnSn3XsuYWS9x8II5zDz2bF77ypls7Nk328WTHKL6K4t27YIXXuDCu3/LmFkvedsOPxyeeortr67PbtkyaeJE7n59OZMfuonD/vMqB3w6j9ePPp3Xv3IGoMQq3zWbWJlZEXAXcAJQBswys2nOuY9jdrsAqHTOHWBmk4BfARMjKWFFBXz2WSShWky8pl/n9n7U1noVz86dsG0bbN0KFRWc88L77LP2MwavWEiPynUNIT4aeSQPfv+6vG2pamzeqKO5+ZqHOP/BG9hv6XxOe+YhTnvmITb07MuyoSOo6L0vm7r1YWunUmpKOrKzXQl1RcXUFRdTb0XUt2kDZji8MQ6u8VCHbIx9mJOHve/Dh0PH1tmamPX6q74e3n8/klBZEdRlsc/Oed+rvt6rw2proabGq8OqqqC8HFatgg8/9O6Oq65mDFBbVMQzp3+Ps574PbRtC3maWDVcGI8+jhWDD+J7913HsMXzOGXGw5wy42H4/QAYMwb23x/69YNevaBzZ+831K6d97cpKgK//trj0ZjGb6Vvn31gQLRzqIVpsToCWOycWwpgZo8DZwCxFdMZwM/9138F7jQzcy6CzuV//AOmTEk7TC6JXRBhW4fOzDpiPK9+5UxWDB2RtTJly5p9h3Lzzx7igE8/4IQXHmfE/HfpuWEtPTeszXbRUnN9tguQAXPnwqGtdmqM7NZfO3fC6NFph8lpw4fz9IhxvPOlk9jYsy9ntc3sXcytyfre/bnlp/dzwOJ5HDPz7xz2n1fpUFYGZWXZLpoErrgCbrst0pDWXN1hZl8DJp6b9voAAAO3SURBVDjnvue//zZwpHPu0ph9PvL3KfPfL/H3Wd8o1hQgyJIOAhYlUdZeQH5e4uxWCN8R9D3zSbLfcbBzrnemCtOY6q8WVwjfsxC+I+h7xhOq/grTYhWvrbFxNhZmH5xz9wH3hTjn3oUwm+2cy+tLv0L4jqDvmU9y4Duq/mpBhfA9C+E7gr5nOsIM+CgDBsa8HwCsTrSPmRUDXYGNURRQRCQNqr9EpEWFSaxmAcPMbKiZtQMmAdMa7TMNOM9//TXg5UjGJ4iIpEf1l4i0qGa7Ap1ztWZ2KfAc3u3KDznn5pvZDcBs59w04EHgz2a2GO9Kb1IGyppSE3yOKYTvCPqe+aRVf0fVXy2uEL5nIXxH0PdMWbOD10VEREQknDycVEdEREQkO5RYiYiIiEQkZxIrM/u5ma0ys7n+4+RslylKZjbBzBaZ2WIzuzLb5ckUM1tuZh/6/w1nZ7s8UTCzh8xsnT8fUrCth5m9YGaf+s/ds1nGKCT4nnn9u4xSPv+tVH/ltkKow1qy/sqZxMp3u3NulP+Yke3CRCVm2Y2TgBHAN8wsn6dZP9b/b5gvc6RMZc8J8wGuBF5yzg0DXvLf57qp7P09IU9/lxmSd38r1V95YSr5X4dNpYXqr1xLrPJVw7IbzrmdQLDshuQA59xr7D3v0RnAw/7rh4EzW7RQGZDge4qo/spxhVCHtWT9lWuJ1aVmNs9v0svpZslG+gMrY96X+dvykQOeN7M5/hIh+Wof59waAP+5T5bLk0n5+rvMhHz8W6n+yk+FUodF/ptsVYmVmb1oZh/FeZwB/AHYHxgFrAF+k9XCRivUkhp54ijn3GF43QaXmNkx2S6QpCWff5dJK9A6TPWX5KqM/CbDrBXYYpxz48PsZ2b3A9MzXJyWFGbZjbzgnFvtP68zs7/jdSO8lt1SZUS5mfVzzq0xs37AumwXKBOcc+XB6zz8XSatQOsw1V/5Ke/rsEzVX62qxaop/n/YwFnAR4n2zUFhlt3IeWbWycxKg9fAieTXf8dYscuknAf8M4tlyZg8/11GKo//Vqq/8lPe12GZ+k22qharZvzazEbhNTEvB36Q3eJEJ9GyG1kuVibsA/zdzMD7f+8x59yz2S1S+szsL8A4oJeZlQHXAbcAT5rZBcBnwDnZK2E0EnzPcfn6u8yAvKzDVH/lvkKow1qy/tKSNiIiIiIRyZmuQBEREZHWTomViIiISESUWImIiIhERImViIiISESUWImIiIhERImViIiISESUWImIiIhE5P8DPdmfT9A3Af4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Distribution of guide after optimisation\n",
    "new_speed_guide = [guide(guess=5., guess_sd=2., time=10.).item() for i in range(1000)]\n",
    "plot_hist(speed_guide, new_speed_guide, \"Before optimisation\", \"After optimisation\",\n",
    "          true_dist=(3.4,0.89))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is that after optimisation, the distribution of the guide looks almost identical to the exact posterior. So we have a good approximation to the posterior distribution which we can easily draw samples from.\n",
    "\n",
    "This example used stochastic variational inference, which gives us a way to approximate a posterior distribution, but we could have used Markov Chain Monte Carlo to learn a posterior distribution exactly at a higher computational cost."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:infolab] *",
   "language": "python",
   "name": "conda-env-infolab-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}