lectures/Lecture_10.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Bayesian Statistical Inference: III\n",
    "\n",
    "*S. R. Taylor (2021)*\n",
    "\n",
    "Material in this lecture and notebook is based upon the lectures of A. Connolly's & Ž. Ivezić's \"Astrostatistics & Machine Learning\" class at the University of Washington (ASTR 598, https://github.com/dirac-institute/uw-astr598-w18). Also the \"Inference2\" lecture of G. Richards' \"Astrostatistics\" class at Drexel University (PHYS 440/540, https://github.com/gtrichards/PHYS_440_540), Karen Leighly's [Bayesian Stats](http://seminar.ouml.org/lectures/bayesian-statistics/), and J. Bovy's mini-course on \"Statistics & Inference in Astrophysics\" at the University of Toronto (http://astro.utoronto.ca/~bovy/teaching.html). \n",
    "\n",
    "##### Reading:\n",
    "\n",
    "- [Textbook](http://press.princeton.edu/titles/10159.html) Chapter 5.\n",
    "- [MCMC Sampling](https://twiecki.io/blog/2015/11/10/mcmc-sampling) by Thomas Wiecki.\n",
    "- [Sampler, Samplers, Everywhere...](http://mattpitkin.github.io/samplers-demo/pages/samplers-samplers-everywhere/) by Matt Pitkin.\n",
    "- [MCMC Interactive Demo](https://chi-feng.github.io/mcmc-demo/app.html) by Chi Feng.\n",
    "\n",
    "---\n",
    "\n",
    "## Contents\n",
    "* [Monte Carlo Methods & Markov Chains](#one)\n",
    "* [Markov Chain Monte Carlo (MCMC) Techniques](#two)\n",
    "* [The Metropolis-Hasting Algorithm](#three)\n",
    "* [Worked Example](#four)\n",
    "\n",
    "---\n",
    "\n",
    "***Exercises required for class participation are in <font color='red'>red</font>.***\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Please pause for a few minutes and install these two packages `emcee`, `pymc3`, and `PTMCMCSampler` before going through today's notebook. Make sure this notebook is in the correct Python kernel for the class conda environment before executing each of the following cells in turn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "!pip3 install emcee"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "!pip3 install pymc3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip3 install git+https://github.com/dfm/acor@master"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip3 install git+https://github.com/jellis18/PTMCMCSampler@master"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Monte Carlo Methods & Markov Chains <a class=\"anchor\" id=\"one\"></a>\n",
    "\n",
    "### Motivation\n",
    "\n",
    "Consider the problem of estimating location and scale parameters for homoscedastic data drawn from a Gaussian distribution that we looked at in `Lecture_8`. We had a two-dimensional posterior pdf for $\\mu$ and $\\sigma$:\n",
    "\n",
    "![Ivezic, Figure 5.5](http://www.astroml.org/_images/fig_likelihood_gaussian_1.png)\n",
    "\n",
    "It was easy to *numerically* integrate the posterior pdf to find the marginal distributions, and find its maximum, using a brute-force grid search because it was only a $2$D problem. With $100$ grid points per coordinate it was only $10^4$ values. \n",
    "\n",
    "However, **what about high parameter dimensions?**\n",
    "- Even in a case of rather simple $5$D problem (as we'll discuss later today), we'd have $10^{10}$ values! And often we work with models of much higer dimensionality (it can be thousands!). **Brute-force grid methods are not feasible beyond a few dimensions**, and even then can be waste of time. \n",
    "- You could **simply randomly sample the grid at every point**, and try to find the minimum based on that. But that can also be quite time consuming, and you will spend a lot of time in regions of parameter space that yields small likelihood.\n",
    "\n",
    "A better way is to adopt a ***Markov-Chain Monte Carlo (MCMC)***. MCMC gives us a way to make this problem computationally tractable by sampling the full multi-dimensional parameter space, in a way that builds up the most sample density in regions that are closest to the maximum probability. Then, you can post-process the “chain” of points to infer the posterior distribution and uncertainty regions.\n",
    "\n",
    "Ivezic, Figure 5.22 shows a problem similar to the one above, done with a Markov Chain Monte Carlo algorithm.  The dashed lines are the known (analytic) solution.  The solid lines are from the MCMC estimate with 10,000 sample points.\n",
    "\n",
    "![Ivezic, Figure 5.10](http://www.astroml.org/_images/fig_cauchy_mcmc_1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### What does Monte Carlo mean?\n",
    "\n",
    "In case you are not familiar with Monte Carlo methods, it derives from the name of an administrative area of the Principality of Monaco ([Monte Carlo](https://www.wikiwand.com/en/Monte_Carlo)) known for high-stakes casinos and gambling. Gambling and random sampling go hand in hand together.\n",
    "\n",
    "![](https://www.thomascook.com/.imaging/mte/thomascook-theme/og-image/dam/legacy/heros/citybreaks/monte-carlo/montecarlohero.jpg/jcr:content/montecarlohero.jpg)\n",
    "\n",
    "**EXAMPLE:** \n",
    "- You have forgotten the value of $\\pi$, but you know the formula for the area of a square and how to draw a circle. \n",
    "- We can use the information that we *do* know to numerically compute $\\pi$.\n",
    "\n",
    "1. We start by drawing a square and circumscribing a circle in it (actually it suffices to just do a quarter of a circle and scale accordingly).\n",
    "2. We put down random points within the square and note which ones land in the circle.\n",
    "3. The ratio of random points in the circle to the number of random points drawn is related to the area of our circle, allowing us to calculate $\\pi$.\n",
    "4. Using more random points yields more precise estimates of the area."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<font color='red'>Try it below. How many draws does it take to get the error down to 1 part in a thousand?</font>\n",
    "\n",
    "We'll consider one quadrant of a square of sides [-1,1] in which we inscribe a circle (actually just one quarter of a circle)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Execute this cell\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format='retina'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(subplot_kw=dict(aspect='equal'))\n",
    "ax.axis([0, 1, 0, 1]);\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Lay down M random points.  Tally how many are within a unit circle.\n",
    "M = 1000\n",
    "x = stats.uniform(___,___).rvs(___) # M random draws between 0 and 1\n",
    "y = stats.uniform(___,___).rvs(___) # M random draws between 0 and 1\n",
    "\n",
    "r2 = ___ # equation for radius of cirle in x,y\n",
    "\n",
    "fig, ax = plt.subplots(subplot_kw=dict(aspect='equal'))\n",
    "\n",
    "# Plot points in square\n",
    "plt.plot(x, y, '.k', \n",
    "         markersize=3, c='blue') \n",
    "\n",
    "# Plot points also in circle\n",
    "plt.plot(x[r2 < 1], y[r2 < 1], '.k', \n",
    "         markersize=3, c='red') \n",
    "\n",
    "#ax.axis([0, 1, 0, 1], aspect='equal');\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "    \n",
    "# A_square = d^2, Acircle = pi*d^2\n",
    "piEst = 4 * np.sum(r2 < 1) / M # pi = 4A/d^2, where A is d^2 times the ratio of points \"in\" to total points\n",
    "\n",
    "print(\"Estimate of pi is {0} for {1} draws with fractional error {2}.\".format(piEst, M, \n",
    "                                                                              np.abs((np.pi-piEst)/np.pi)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "See also [this Khan Academy simulation](https://www.khanacademy.org/computer-programming/monte-carlo-finding-the-value-of-pi/6530004791197696/embedded?embed=yes&article=yes&editor=no&buttons=no&author=no&autoStart=yes&width=610&height=420)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- In general, ***Monte Carlo methods*** use random sampling to obtain a numerical result (e.g., the value of an integral), where there is no analytic result or it is difficult to obtain.\n",
    "- In the case of the circle above, we have computed the intergral: $\\int\\int_{x^2+y^2\\le 1} dx dy.$\n",
    "\n",
    "Recall ***Monte Carlo integration***. \n",
    "- Our goal was to integrate $\\int g(x)f(x)dx$. \n",
    "- We could do this on a grid of $x$, such that $\\int g(x)f(x)dx \\approx \\sum_{i=1}^N g(x_i)f(x_i)\\Delta x$, where we use the probability distribution at each grid point $f(x_i)$ as weights for the sum over $g(x_i)$.\n",
    "- But we might waste terms in the sum over low probability areas. If we have random draws from $f(x)$ then the *density* of those points in $x$ directly represent the weighting we desire. We will get most weight in regions of high probability. \n",
    "- We can think of our random samples as $f(x)$ being a sum of $\\delta$ functions at the sampled points, all overlapping to give high density at the high probability regions. The integral is then \n",
    "\n",
    "$$ \\int g(x)f(x)dx \\approx \\frac{1}{N}\\int g(x)\\sum_i^N\\delta(x_i)dx = \\frac{1}{N}\\sum_i^N \\int g(x)\\delta(x_i)dx = \\frac{1}{N}\\sum_i^N g(x_i).$$"
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0IAAADYCAYAAADGS1PMAAAgAElEQVR4AezdC7i8VXUecJvUaMWqKQZjC8Yb1cSUxGhoorFqFMTgXVSUO2oARUFBUQmIURIFqsYQiOIFBamIClFrMdFIjRY18ZYYTDXWEmyamNQaGxNLL1+f37bvYfz4z3/uc2bmrP083zPnzHzXtfe313rXetfat+iqlQRKAiWBkkBJoCRQEigJlARKAiWBHSaBW+yw563HLQmUBEoCJYGSQEmgJFASKAmUBEoCXQGhGgQlgZJASaAkUBIoCZQESgIlgZLAjpNAAaEd1+X1wCWBkkBJoCRQEigJlARKAiWBkkABoRoDJYGSQEmgJFASKAmUBEoCJYGSwI6TQAGhHdfl9cAlgZJASaAkUBIoCZQESgIlgZJAAaEaAyWBkkBJoCRQEigJlARKAiWBksCOk0ABoR3X5fXAJYGSQEmgJFASKAmUBEoCJYGSQAGhGgMlgZJASaAkUBIoCZQESgIlgZLAjpNAAaEd1+X1wCWBkkBJoCRQEigJlARKAiWBkkABoRoDJYGSQEmgJFASKAmUBEoCJYGSwI6TQAGhHdfl9cAlgZJASaAkUBIoCZQESgIlgZJAAaEaAyWBkkBJoCRQEigJlARKAiWBksCOk0ABoR3X5fXAJYGSQEmgJFASKAmUBDZfAt/5zne6v/qrv+r+8i//svvbv/3b7n/9r//V/d//+383/8HrCceWQAGhsUVVO5YESgIlgZJASaAkUBIoCayLBICgj370o93v/u7vdl/96le7//E//kcDQ+ty/3Wfi5dAAaHFy3hpV/g//+f/dP/zf/7P5vX47//9v3d/8zd/0914441Lu35dqCRQEigJlARKAiWBksAqSEDk50tf+lL3ute9rjvrrLO6Sy+9tLv22mtbdKgiQ6vQQ6txDwWEVqMf5nIXQNB/+2//rXk9PvvZz3Z//Md/3MDQXE5eJykJlARKAiWBkkBJoCSwBhIAgjiH/8N/+A/dYYcd1v3ET/xE97CHPax7wQte0MBQRYbWoBOXdIsFhJYk6EVexgvPu/Ff/+t/7a655pru7W9/e3fBBRd0l19+eXf99de3yaA4sYvsgTp3SaAkUBIoCZQESgKrIgGO4b/+67/u3vWud3U/+7M/2/3Df/gPu3/8j/9x98AHPrA777zzGkCSN1SRoVXpse27jwJC2yf7uV3Zi8y7wfNxyimndI94xCO6n//5n++e85zndJ/85CcbXY5npFpJoCRQEigJlARKAiWBTZfAN77xje4P/uAPGugRDfoH/+AfNDB0pzvdqXvAAx7QIkNspooMbfpIGP18BYRGy2il9xDpQYfzwv/ar/1a93M/93PdP/kn/6RtwNDb3va27j/9p//U8oZW+kHq5koCJYGSQEmgJFASKAnMQQLYMNgxz3zmM7u73/3u3S1ucYu2JTIEDFVkaA6C3oBTFBBa404MB/YLX/hCd/bZZ7dI0N577939wA/8QNvufe97t6jQFVdc0d1www1r/KR16yWBkkBJoCRQEigJlATGk8Af/uEfdi9+8YsbLW7PPffcAkIVGRpPfjtpr6mBUPJS1GUXkZCfwtiGwr/2ta81buaiarYPXhvH0/Vs7sG9LOq6qzQwyOBb3/pW9x//43/sLrnkku7Rj35090//6T/t/tE/+kdbL/wP/dAPdf/qX/2r7vTTT+9+//d/v/v7v//77n//7/+9So9R91ISKAmUBEoCJYGSQElgLhKQKsAG/PCHP9w9+clPvpldtKvI0Lnnntv2V2QKgPr4xz/ets997nOt6JTCU9lUoWNvsj2V5k51XjZZtfWUwFRASIfLOWGIf+UrX+k+9alPdf/u3/27TuRBecJ3v/vdLV/Fb/ax7zwHSXJi1IRXG/4973lP29yDe1nUdVeliyN/IOj888/vDj300G7fffdtIOj7v//7t4DQrW51qw4YetSjHtX65M///M8bGFqV56j7KAmUBEoCJYGSQEmgJDAvCQBB0gHe+ta3dg960INuZhcFCCUytNdee3X/8l/+y+6pT31q95KXvKR70Yte1D3hCU/oHve4x3XPf/7zu1/+5V/+nk0KAhuX7fmxj32sASTLlczbzp2XPOo8oyUwMRBihCcS8ZGPfKTloJxzzjndaaed1j372c/unvGMZzQ61q/8yq+03+zDYHfMtGDIAEsFkOuuu64Nvg984APdW97ylu6lL31p99znPrdt7sG9yIuZx3VHi2/6PchCdGZw3Z+sfgyw/Nmf/VkDdH/yJ3+y5Ykgx//8n/9zi7rxSgCAhx9+ePfP//k/b9VQ8oL3P3/sx36sycmLy4sxbT9M/7R1ZEmgJFASKAmUBEoCJYHFSuDrX/96s//Yhmyfvj3U/5/z+Da3uU23zz77tBzrn/mZn+lQ6X7wB3+w+6mf+qlWclvZ7WyPf/zjm43r/FISVOh9//vf35zw7DXMpG9/+9sFjBbbzXM9+0RAiAENlCQSceyxx3YPfehDu/3226+Tj3LPe96zu9vd7tY+/8W/+BftN/uIWjhmWsSc9XFU+Hj5y1/eIiAPf/jDu5/+6Z9uA100xOYe3It7msd15yrpgZNFjqhqyjuKbAnJ8i4AK//23/7b7t/8m3/Tvf71r2/ALh6Jf/2v/3WjwV122WXt5TvppJO6n/zJn2wgSAJg/wXP/zweZGJBsT/6oz+auh8GHqH+LAmUBEoCJYGSQEmgJLBSEvjTP/3T7sILL+ye8pSnNHATO2jYp8hQwNAd73jHBoKSZw0MsZ8Gt3/2z/5Zs3GBrB//8R/v9t9//+6ggw5qRRnYaJz0UkTYrWzeaqsvgYmAEJQrUnHllVe2SMSP/uiPNuqVcoSADyRt87fv0LLsI2rhGMc6x6QNBxM/881vfnOjeQFcQI+SiOrDC3/69P+P/MiPtMFvgB555JENVPzFX/zFtlLCBoFPFjxV5e13fud3WlRHZOu1r31t8y7wMljw67jjjmuh2l/4hV/Y8kQ85jGP6Z7+9Ke337zkKsR5cYe94Ple3pD8oSc+8YlNHsCXl7RaSaAkUBIoCZQESgIlgXWXADtL2gTb6uSTT27RHEAmdtC4nwFGt73tbZstea973avZm2xO7BvOfvaUcwNMqULnN7naL3zhCxstDyvpM5/5THN0s2Hdm3ustnoSmAgISRDDjXze857X3fe+921A5y53uUsnOoMKJ4ph87fv/AYM2dcxjnWOSdt/+S//pXvve9/bAAKwY0AeffTRrfShxbLkBvl89atf3UAXhC7MCZS5l3//7/99J1y6XS3UPoBMDhPuKhohkKOYgciWyA7wxsvgxVPukfzufOc7b3kjvHxeQr95PiBIHtCoF5y3AxjCg33Na17TJgo19quVBEoCJYGSQEmgJFASWHcJJHf8gx/8YHOYB6iMso/6v7OXbn3rWzeWkfzrU089tVWfkzskZ4gzmlMabc41kmtksVY2GvuNDWofdi8Hvsq+tV7R6o6wsYEQJItWdcYZZ7TFOhnk6GgQMLoa2hpupC0UNr/Zx77WtHGsc0yKitHqhBwf+9jHNnBwv/vdr4EeyF90IxQz6Puiiy7qjjnmmO6ud71r2/ewww5r30meW1YL8OEFAOJEs4Cxd7zjHZ3qJO4P6PHi8CbkRcJTtQZQwrAJwfI0iHQBleMAn/6Lnf+d44gjjmjykGNUHopljYi6TkmgJFASKAmUBEoCi5IAoIEW98Y3vrG7//3vP9JBHLuo/xkgxO7CZvrVX/3VTjqCPKCrrrqqpSewR4EcYIfDXUoGm/P2t7/9lk0HJAFLmEnSQ1Sx++IXv9ic8mzWSe3gRcmtztt1YwEhHca4l8NyyCGHNGCD/gYtK90M7Ur017k2f/vOb/axLzDkWOdwrkkGAcBjUSzREogbH/N973tfJ6qB4pWiA/73IkDgoixAhU9Us89//vNL6+/kNAFAIlmvfOUrWxnHYS9MQqsiQKI2/aQ8Xggv04Mf/OAGhvov7rj/A1rAEHAoF6k8FEsbEnWhkkBJoCRQEigJlAQWJAFFpoAVERxRmXHtov5+ocYBMmzOo446qgEgtqXcHwWrOOc53tHfhjm40eacg82FFSWdQcU5TnHsoEnt4AWJrU7bjQmEAA0ABwVNPg6AwWDHhZTk3+/QACe/2ce+jnGsczjXOOvZOC9QwWg/+OCDG03MAFcIwbo4/ZbrGmhKRougABfKIl577bUTga/+ucf5/zvf+U4DgQCQhDnVRICYRz7ykd0P//APb0V+bne72zVq233uc58mG9Ey0TNgh7xSHAGNDdWQh+PEE0+cGQjlhfdS8lAAq6JW1UoCJYGSQEmgJFASKAmsowTYfqItnM6xt2LvTPsZB/UDHvCAloaB6aTqbpg0sU8HUx7YaVIe5KxL4WCDSksQZbLYvTxt6RoKYgFT3/zmN29mP6+j/Nf9nseKCAEjKGjycHSwzhWdOPPMM4dS3QxMNDj7JJLhWOcYN1k/kRWFFuTSoIrJPQIUgI1dNdcFkn7xF3+xFWpwryJI11xzzcIHnEjYRz/60Vb44ElPelILiwqvAkE4p3mxFJDgHVDF7e1vf3srmuC4T37yky1y5dlsgIrFvUTWvNyeZRZqXCYEHgqgS+RMxbpqJYGSQEmgJFASKAmUBNZNAnGAWwSVXRV7K/bOtJ9JWVD4CxjCLAKGwqTJdTn2FcESKRIlkqN08cUXtyVl2L5YTMCQDRjCDGJ//fqv/3qz89i5QFW17ZPAWEBIR8v9EZ3QiToWNxJPEqod1pLbY1/HONY5nMs5R7U+ADMgDzzwwIaoof9hTSRKaWmFFe5whztMHIkadt7+914E3gGLaX35y1/urr766sYnhfoN+BQpUNRArlTKLB5//PGtvCNwJoeo/yLkvORkH8BP9GjaF7p/XLir1l8yeViAzHNUKwmUBEoCJYGSQEmgJLAuEujbiX17Z9b/Q3HjUL/iiiuazfZ3f/d3Q8WTBV2xggQCLM4qZ4lNCAzJDecg5yx/wxve0Bz3okp9O3DoBeqHuUtgKiAk30flDJEKKHhY85t97OuYSYFQn5InEU1SGgMe2BnW5BSJCKHRGXSTRqKGnbf/PfDAO/C5z32ue93rXte4pJ4xA96gBwCVuXbPWXhLxArfVCU7L1SfWpjz8j6ccsopW9TCWV/oHJ8X+4ADDuguvfTStgqzl7daSaAkUBIoCZQESgIlgXWRgNxwNt95553XnN+xc+b1GXtJCoP1ieSbo7QNa+w39pScJYye3/qt32pLo8RBHmaQ4goiRiJN0jlElYChasuXwFhAqF+NA5qV+6KDRTSGNb/Zx76OgYrluwABzjmqJfT4e7/3e93jHve4liOE1iX8KUmtn2uUAag6x5Of/OStkOSkAGzUfeU61kVS/OE3fuM3WiEIzwh4ocEJzyqJbRXi008/va0XhO4metQHPv3rZb0m0TOFE1Jdbl4vds6TflTQQV+R96o29za4rep91n2VBEoCJYGSQElgkyQwqHtXzU6Iw/0Zz3hGc7jHvpn1M3acCr/ssOc85znd5Zdf3uhs4+RWkxNbT7QH0JHzzY5ld2HlcJRz7suhx/pBqVPduBg6y39zxgJCcnpEJ5QRVAFOEpia6qhg6FvDmt/sY1/HONY5nMs5x2kGk+gPMKUCnHV1HvKQh7QBCXEPUuwSkrROj8IMBhp62ryBUK4D5In0yF+C7gNYgCA5PS95yUsaAOJBuOGGG1phghtvvHEk4OAZ+MQnPtGdc845Lc+IBwFfddYXu3882qDCCRYf+/SnPz0SoI3TX4vYJxOKCKFtFJBcxD3UOUsCJYGSQEmgJLDTJLDq+hcjB3OGox3A6Ns50/4fO479qny2HG5V44Agdtw4jezYqNg/7pNz270qqw0EsVEVEpP6IG9IDn3A0Djnr33mI4GxgJBOhGiTqyLZ/5d+6Zdakj+0O6z57Xd+53favo7R2c7hXJMscCqCJCQpEmQhUZQ31dWgc5EhVDObSBAQxLB3vQAIQMi+Kn6o7DZtM6CBO2FYL8auBrTrWu/o7LPPbhXyRFomNdwT6hU983Kg9sl3spDqnnvuuQXwpn3Bcxz5WD2Zt8MLCqztjvs6qdwygZo0nNcEoqAEACua9pWvfKX7kz/5kxZqNsmIrhkb+nRwUzVQ337oQx9qm799N7iPvx3rHCk6IUeNt0gf6DfgWyRS+Bmgcn/VSgIlgZJASaAksGkS2HT9mxSC3/7t327Mm1mLSSlE5Rz3vOc9m9Mdq0gVOrk+ANCsOTx9B/oDH/jAxhyKnSqPnGP9ne98Z8u9r8jQ8t7IsYAQw1zEQHEEBvm0QMixzuFczjluA5oYui996UtbXXeRF2BIZAg4Uk7bZuCKBLk/uTmJ0Kj48e53v3tmDmYiXECOctf9EKeIVWhwolj2BwAmNbi9cOQDKCj7zUuA/2otJaDIs/EkBNBM+5mqKGrlA3WjqI7j9lf2AwA9Cz7t1772tRZSBlSUjgS8Xv/617eoF1AL8JGpqoD6dXDDo93VNriPvx07WIbcWJOj5rlEJkUigWqACKidFKDmueqzJFASKAmUBEoCqyyBTde/SdlQ/RaoAGS+7/u+b2q7CAjC7lEC2zk5VqUzsOOkK8xqL/RTKtgn7jt2Kqc0MGS5F2CoIkPLe7vGAkLf+ta3muf+N3/zN7v73e9+DQBYVXeU4ZwcIfsCDY51DlEA5xy3BUkDMwCPCI+iCbibqF2qsRlQihLgWwJc97jHPbq99tqrvRy+twCrF8dgnLRlAIsEAUGpUw/JG8RJegMmeCf6lL1Jr5f949FhuLu2UtsvfvGLtwAfACNCJhwsoW9aQEROEgGV8wbgRM289OM09xj5uE+A50tf+lILAwMeojjGiaIMcqle8YpXtOTA4447rr3wau7rM+s9ZUKY9jnSH1mYVrVChTrQKoW30TJF2QBL61mZ6EQSTXbKiItgAm1AkohRtZJASaAkUBIoCayqBHay/mVnzbKAaoogyDtnM7IXMJ1UhgNCZo0ADRsz7CW2KCe3Qgns0xTYkspRkaFhklvc92MBoSTvv+1tb2ugw8A54ogjuje96U1twAy7PYPJPvZ1DMDiHGhRzjlui6HtOOWes5Iv5G7wKqktGiSs+OpXv7oZ9L5jEN/mNrdpKB9Fz3kmjc64xwAxdDgRh9SpZ7gDeINlEAOC5mVIu99EiIRnVSEJBVCETEQM33QWbixPCjCE0qdwAvraOBRC9wYwAbWiV4APsGr1ZOCDZ0M/AKnAjvsEXkXz0Pz0j5wvHFn9lBDxtEAoES7nck6RM9UK9VFy1PCIRdVEDi26q7oguuZb3vKWBmLlc6UPxx2ftV9JoCRQEigJlASWKYGdrn+ts3jGGWc0XW6JkkntBjYTm+Twww9vdiNQZR0gKQJsvlkjQMPGgn5ji0rVAIbYS494xCO2mD4VGRomucV9PxYQYhTrNHk2ojEMTEYz0DFqHSH72Ncxs+bqGJhAAc/9pz71qXY/IjToaAoLiJj4/j3veU/jjDK0hTsZvCrPTdoCwDyjUCWgBa0n8oCC97SnPW3pC2MFmKELyp2SE+UZGflARoDapBODHK6Xvexlbe2iwRyuTLgiJUASQIpiBpShugFPAK8Vk1VWUSnPBNOn8MUD4/78BhwDKa7LIyNiCCwbJ9NsqIkq9QFbomW4vhbhNQYAo37EKUBWNOyoo45qxS2UQUfbA5xROD2nsU/m0wLpScdd7V8SKAmUBEoCJYFIIMYzPUQHy321YP0m61+5vgoM0MEKSIWpETvQIvQcwZypKryNsneyrqP1KBXuwkZhO2GJAEDLLl/NnhAZCtOIDcdxmyJfbKOkLLC5lml/DI43ckEPBBA547F+sH/mZRPFzpYOwu5k34+z2dcxuY+8K9N8jgWEgA/C4O2Xb8OQhWAZvrtb2NRv9rGvY2bN1Rk0yN2PwXHddde1CQHdzuRgkkB7skYOAxhY4fU32CZtBCyqBQSJbgBB0HoMaCCI0QwooVQtyoPQv+8MHAPBJOHZlF70Qp966qlb1L1RE0P/d4UYREvQ74CcNM+l4IGXweQjIgeAAaBq44v4AGEml4AP3pa80LlOPDDGA3AsUog2KSKDL4s2KWJIptNsojqvfe1rG31RtCygTH6RiJSJRd/lfgJoRcPQG40VVEuhapPkCSec0J5TYQbUuWmplZFjfZYESgIlgZJASWASCcTuCfOCU5eeFA055JBDNk7/siekH4jUAAJsDVV0GcfAUOzRq666qlH6gSAgJ3p92Cd7hANWLjG7VLoExzkbz7ImzsvWWVYL2GDLim7JSwKG4kCOvcQBD/Auy/7ojzcyYlOjDEpxgAOwf7CAjMlZ7d7Y2RzP0hU4ocfZ7OsY/eccs7SxgFCQq+oZBhEvO4PZIIXag9RzI2hhvhP2O+mkk7bKXjvWOSbtUOdihKfKGKNUBbI+QjaQIVU5IO4P+LIasMofuwNsue98BmiMigSdf/75DQTNOhBy3Wk/0z9AkfCu5+UpIQORlnGrzd3ylrfslNQGVIA/uT6AJm/JNddc08DJy1/+8u5Zz3pWAzEKFACbmXgS8TE+ACLACEACgCUhqqGPMmcSEikURZI/lAENzAK3+nqazbiQ42TSAF5C0wOKTjvttAaITTTJJQN6gB8gCBhKxCgAyW/k6Hggy9h1bpOW+5uHJ2LaPq/jSgIlgZJASWDzJBBDVKGlODvlstI/9BAdLPeVTt1E/ZvqbRzPGCJ0MJsGI4nDV2qAgkvYQNgfsT/6n4kAkRFqPDq8FAryY8gDIGzG7W7pb3aFZ2J/JWdIiW3AkKOYPc3Rv0h7070AN2xfjCOAG9uK/fTsZz+7s1YTBzMbzm/2sa9jHDtNM8adBxCUxy1Pa5zNvo5xrHPM0sYCQh6Qsc0IPPTQQ5uRi4IlNwYq7+dUAC6+gyABEaFLNCjHOkcfwIx6AEanql+pMubhDZo+oBImEx3JCsMJLWbB0FHXye9BqONEgmYZALnerJ/pH/Igd6BPLtWk1eZUXAGGUNUAlkRXgFkRHMDKxJN8G2BHhCwTUDwYqHFelhQnAEhMYO5JqNv9eXkgeQU1DGJ9R5Zyx4yfaTb9xrNjskBnS+EGY0WEizIBukyiaJRocNZ6QotDj0vEKLlGt7vd7VoxCjQ7tDtA3hg+88wzm4dkHp6IWfu+ji8JlARKAiWBzZEAQ5eBTjdyQIqIHH/88S1KEh2MQsXABxo2Tf/K8VUESsEr0RF/c6pyqNLTokVsDEyPQSAYOeQzESD7MeQtrQJEce5ymC87ArS7EcqGA8zYr3KGOGuzzhCHPjDkGfy+qPt2D8Ye24yTX2Gyhz70oa0wGZuPvWfcxcntN/vMGhDgwDfGjznmmPbcHNDjbGTkGMc6xyxtLCCUCzAm0aYga170lItm6DIyv/zlL7fN375TFhlgMlgNRgYkozRNhzKAGZTCbwlT+s5vaaIEqFIqfwFfRx55ZBO+Qc2odrwkd+E7q/cSDoPdNd/whjeMvRJwJiDnAyKG5QTN2vF5rkV9ZkB72b04k1abEx0xyEVwUAzl3QA5WdRV5ES/yumRCyTq40VNxMeLrP95XQzQVfO8AEyUjDGjyh8vG3qe8uQmWeMa+JFjJrcoFL9hE6sxXRGiRY3mOm9JoCRQEthcCURfcwLSI6qYctjR2yq5iogwDMNYSMSEUcrO2TT9K9LF7gB8sEuif+Og5HyV59MHgvkdgGCfosOjDq6L4zLMnlSTY1MBQZ7fM8WeBeQwoubdUhTtyiuvbLREY46Mk1OVvG2gNPK3DwqjY9jpkxRBy/0b7+wv61myPbGS2Jz95VH6/8feFJmSFjNLmwgISZTycip9zHvOSCQIdCP5IjzsNn/7jrHs5YUmoVlAxYuelggOihSvvM3fjHe/pYly8OI7r84QXUJrMkEwui+66KJWhtCgN/hNDn4HiFCkxl0JGPhitIteAFu7ywlahUhQ5LOrT5NrgOak1eZMtAakSciANzAHy3N7OQEGQFHUR5+JuuEuo6ah1CWhTpTKfQCZq9DIxYSDdiCvC+8Ypc7LiIcs8ggYocMddthhDRSFszss1L4uE+0qyL/uoSRQEigJlARukkAcsPQ0+hsaGAoXWwY7Q0QEOyFVVRmnwAJ6HJrXpulfeSioYHJ0gaHoX5GegB15QeyUwXWDQmln/7H9OMXZnOviqIxtkmpyWWcokSHgTiCAI17QYd6NzcaBLW+b/WyccQYnpyp522hxvksxMvs6xrHOMWnrAyEgSN+hvO1uC8MICJoGgA3e50RAKMDlggsuaMnk4THmxXzCE57Q2cJdhWTtA7kRYp9KxAiVFIUPCFDZ/O07v6Xlur/+67/eUKPreQmcW7J+DFa5MLwkPCRQ6iQDxiBEqRLF0NFQb1465wS6Vj0SFHkN+wzlDzi0po8iEiI6AO0g0EloOc8fiphJOV4WEzDOrsE4r6S5Yfe9rO/jkQGMyCgAGx1OSF6+FRAuOkZe/QgRzwRPHiCYKNiqAMBlybCuUxIoCZQESgK7l8CgoxJrghMRdQuLhpORs3GY/mVjKZRA3zhu0/QvQCiFgjw40Ifp39gp+UyOMsoUBy1dPGvuyO57cTG/xg4hA85mNhqbQ1BBAELAgMOWvTrOMifj3KXxqAqhccXWYfMKBFguJjlVghg2NrLv/GYf+zrGsc7hXJO0PhAS+QGAltkmAkKJMKCwQatACDCSUK0Igq0PVKA7uSEMcZ2cNi4Qcl00LwPbi9EHWqEw6ZCDDz64RY6U0B43hKjjGKw8BzoY2hUJiYcBTUxe0qxJYXnu7foke31A7goTkJHCB8ANMJQJJZ95fgAVEBSmFxHUD2QlFKpf0q+TvgDbJYdh13X/kRGvTCiXKuXJhVPgAcDOxNSPEIUzK7fKOxIwNOx69X1JoCRQEigJ7CwJxN4IdV21VxVUFQbAogGCRDyG6V95Q9IUUiV30/Qv7z4j35sAStkAACAASURBVPN5zmH6N3ZKPuU3c9pK3cAU6jOL1mWUxQ5hXykmpSgZAJTxgBUF6JELOc3aMh4BL6wqdrRAAAB6ySWXNFvPdZK37W/2n9/sY1/HONY52NLOOW5bOyCUByMInYAGJxoDmAhhhkPob9/5zT7DOoyhyGBUNplAbf6OEZnrpaMMDOE3VCRCdw3oH5fSp3CysLLwMq/CuJQsHQwcGHS8LQnFAgdAAi8NytekHZz7X5VPJbDRBMkGDQwFTF5MgFA8UOHYygFSwUSipoQ0E7DcmnHluirPPet9eF5jVfU8SsvaA8pr9yNE8diIgKq0kgXahItN7pNMDrPecx1fEigJlARKAqsjgVTT5WTjOUf35yTOgvOlf3fdV+PqX45JhRb222+/RtVKxIxDE0V/3fQvgOveL7vsshZ9CUAOE4rDFYiY1R7LuEQjBCJH2b2xx9nEbOPYj451Dva0c47bgCqpCHKgBFFWPiKUBxOOA4bwFBnUcmoAlHAI/e07v9nHvrsK4QWAQP0mBZu/gRK/DTbCZ0wyKglOeA6FTi14AMaxIhUGhvCdfccFLsKnjlW6EK0uyXmoYJK4/GafdXuRBuXnbyBIlMw6QACQCjQqs4QaF49DOLaifmSr+IXwPRnIrRlXrv3rr+v/ntdkoyKdiQlQN+76EaLIT1TUZPyYxzymUT29D0D8TpPbuvZ33XdJoCRQEpi3BNg08p3RrrErMGoYkejWDPjoj9K/3yv5cfUvIAkMceSKoEi+58RFcWe/rJv+DTtFZIttGspkUk7kC3H6z8o8YduITmbZGWAkxc2GUd3Ywn4TlAiAYU86h3M557iNzY/pFUC1NkCo/4AeWmeEQzhux2SAQ+vAkm2S5HoDxf4MfMDHfTjnuC1IWEdA16H6zXugjXs/i9jPgDUBAzEWn7r44ovb5CAEb+I1eQjD8zaYgFVrUS8eP9cLaFCv2wSyCDkOnjPjvR8hSjWVAGkTiqglni+vDg+Kfph0nA5eu/4uCZQESgIlgfWQwKD+ZWeorsX7jUYd5knp38n6clz9CxBxSCp8xZ5ZV/3LrhZpsdSHQgKxL9gb8tkFBdjO0zrq2YeuIZCB1WVccuJyhEsHGdb8Zh/7OsaxzuFc/UDGsHP4np0pX10KBptJvpG0jaQmcMTb2K/AV5zyoYTu7tzj/jZRjtCwkwbQeHjbuIaejnMsQCNiZPP3uIa34+2P8gXUjHtcnsO98tCYnPogyP/zCj3mesv+jHxF2ES1RDAMWlEgNC4gCL8WCOJtUBrd2kkGpoIBilToy2lfsGU/77Kul/HejxClmkoUXHLnvOC8U8LI+mFcR8GynqeuUxIoCZQESgLzlUBf/8bZamkGVVlj0Jb+nUzuk+jfO97xjq1q8TrrXywcKQnAgfWTYl/IY5fPLq8dS2pS+zdSZwcPAiH5PioSygHaXVlqv9nHvo6ZFggp+KHKNBAk358z/txzz90qVmWtIpt8dmBQmoZiCoASMDSPNhcgNI8bWfY5TFI63zoyynIzVhOa1qHzTEbbjmcDEBnqiiKgDTLSA/ZUV9lzzz3b4OVhkBfld14HL9y0L9Syn3NVrhcPVaqppIx7n/Ig9ItqKZyNwjmMMroqz1X3URIoCZQESgKTSYBtsTv9GzsjFOrSv5PJt7/3MP2r3LSlP9AORYfWVf8GUAM7KHIHHnhgK+YFSANFxo8Uhl2llPRltav/sapEWeQ/Sw1BK3z+85/fzskeHNb85rr2dYxjncO5nHPcBtSgw+V5rOOovkCqQVt+R766WgD2U7hLVFX0SRG2edhROxIIZWAJlZ566qltrRh5Mkl2X0R5wnEHxTz2MwkbiAaJ6ilCw8CdJDuDDQhS0OIZz3hG+x1QsrgoYMj7QD7VxpdAPFRohHKwAMpdJcEaY/LOcHt5B4cVERn/yrVnSaAkUBIoCaySBEbp39gZKapT+ne23tuV/k3+LjAkb0hFuXXWv2wyTBJsnfPOO68VakrVWgAEQEou+6TSZLdw4nL+o9sBHi960YtaiW424bDmN2W87esYxzqHcznnuC1AKM+jn0SGUg2a3SrFAAi65z3v2dbVElX1/rCz5mFH7UgglFAjGpgwnImJl4aHBrpW6W4RC1aNOzCm3c/LYhLOglzWChKdgNY9IxBkUElug+KFNYGlompNK/FdH5d1ryRp8mqoLieMTf5e9lR9UeHwQx/6UAs/C/EWAN21POvbkkBJoCSw6hIYV/9ySnJOMuJSIVdUo9p8JED/yil5wxve0D3taU/bGP0bgK3YBrsijJPYE5zeX/ziF5sNOIktIW8ZiMJW4ajFjtrVep793hGBGlwHtF9crL//sP+9A4IPrmsDqIC7VJ8GeNirqdYLBCmTnudmr1u7cZJn7t/LjgRCkCyvvZwYBQJSNABNjKceOEArW7eWF+Xaa6/tXvCCFzQUzSMSoJeBEwMcx9JzmoR5VarNRwLkaTL2cpq0eHCEs8PtDSBV8j2FKfRFURLnI/86S0mgJFASWLYExtW/aOiiQGjrcUKW/p1fb22q/g3QZldwsqbseuwJlLEPf/jDjQ1kLI7b2CoKESh8wGk7LRByrHM4l3OO29g+b3rTmxr4AsBElVDsUn364x//eLPJU63Xwq3AUJ4bKBSFmsV+2lFAyEAiLPkZwEDKEeKQQtcEynBFK5tkII3b4Yvaz72KKEhes9aPZ1PKEC/WYBFqNLhTAAL6h+ZnGTiLepZNOm8Uo3C20LVER2A76wFYLdrirEqZX3HFFWu/YO8m9V09S0mgJFASGEcCpX/HkdLy99lU/QtksCmADrkzceTf97737c4///xWOEEl5XHbt771rQbKreHpfBhEz3ve88bOEbKvYxzrHAC+c47b8jyiSzagRp5RnAQ5D4DrO+szKvqluIKCVBzKHAuzVJHbUUCI4U+YgIDFV2OQAkHQNZQNbXuBZgmzpeOW9WkAQNUQ9DHHHNNCi/KAgCBgSPKZghCq4ynhKRSqUsg6PeOyZDnP68SD4+WV6Kjqi6onAeCKVgBDQsFAuElMScoCqPPshTpXSaAkUBJYnARK/y5OtrOceVP1LxtWDo5CBRa8D+OHHYty+eY3v7lV/R1Xdpaesc6hBWj333//FhRgD4vSsCuHtURyEplyrHM4l3OO2xLB45y3eTbBCN+zhdJiv1u2RORI3hAwpAjG29/+9lboS9rLNG1HASGdc/3117ea8jotKzlDstA1lD1JSG8agc/zmHiiGM/ynaxXs++++249F6Cn0gbEboFbpcIBoGrLlYAJ2UusygpKpnKXJrBwfE1kPCqHHnpoy9sCmoCndQPky5VqXa0kUBIoCWyfBEr/bp/sJ7nypunfPA9AwLbDKuFQDcNELg0aWR9IDJOZZWvklV9++eVb6wg99rGPbcutjFpHSFEK+6L9y31zDudyzkU19zR43WkB2OD97SggJDdI6WLruWRBUUYooxS6hkQNnnVp8UQBQU996lMbCLrtbW+7VQYcCALwcCy/9rWvNRBkvaVqy5eAyYu3whiUw6WQRTwpKacKxD760Y9uNDqgfN0omsuXal2xJFASKAlsjwRK/26P3Ke56qbpX8/Dsfre9763O+WUU5ojNQyTX/iFX+je97733YxaNkxubF6OV4wiBQo40NnEctkUXxjW/GYf+zrGsc7hXIu0o/vXlQYitWCWMto7AggFQfdzg9DGFEuQsM7wXBdK0jBPVAxqBRKEDRVMYHSXQT3sVd6e78OJff3rX9+q+gFAACwqI8+KicxvKdph/FYrCcxDAsZS5g8KCzC/4YYbWqScs4QziIG3yGikBbBx2Cku3kP3IVK9LvPvPPqhzrG+Esj702dilP5djz7dFP0bIP7Wt761e9CDHtTsB1VpFS0QMRExGqfol/HMRhQksECr6snABYD1yU9+8mYOdM508zXb8qSTTmqpF3e+853bsc4xrr2Z89BB7G8pKx/72MfaUi4ocq4x6LiPHZ9ImJSPaa67q1G6I4BQuIX93CCAQU6NgfTVr351bXJm8gL0I0GhWCmZLeJgoDI0FmnU7GpQ1Xe7lwBvicn4j/7oj7p3vOMdW5TG1NFXK99aQ8phSjwsA3H38qxfx5NAFIlE1q985SuNPoGqyZt26aWXNm+eRFW/2WdR4w4I4pSi9BSnQeOg+LwXrlmtJLDKEij9u8q9M/reNkX/BpCrFCfnPVVpUe4PP/zwVnltdzk+kVScY+Zj9Hz2B8YUG+Sqq666WUoFgCLNgt446KCDGhhxTcc6x7j2Zs4DPB133HENSGHEcOD/9m//9s2uGzveukOKJQBB01w3zz34uSOAEIFTtIxOPMZ1zQ3KgFUdTggyOUHxRKkMp24+A1o0YVxkPjgg6u/lSWCUQtWXVk/meWSY6v9qJYFpJGDsGEPGEkUiqfWcc87pTjvttBYRt7jyc57znEZ18Jt95jXuXJsSA4Cszwb8WONDJUUrhKs0pIhLzVfT9GwdsywJlP5dlqSXc51N0b9yis8888y2PqTiAdPmCpmDLSkjuiRIkCJbbE3rMpm7bf72XVJMXNPalO7BvaQFcAJjnF02f3MC+y1ASBEvFY1d03qLuW6Ke6EAiho59hOf+ESLdol6qYZsX8Dp85//fC471eeOAEKoFzrh3HPP7fbbb7+tHJp1yw1KCFOJbJGsUKoSCYrhLIogJDouMp9q5NRBM0sgHp1hFIsA26omN7Ood/QJAkSMM2Pp2GOP7R760Ie2udCK4DyAd7vb3dqnCoZ+s8+8xl08edddd107p7mL4kRLth4ET6Lyp7NwvHd0B9fDL0UCpX+XIualXWRT9C+QILLPqWU+nzZXSCEx1ddEZxRvyrIrj3rUo1rVYQuX2lQg9p2oEdvTNV3bPbiXtFAQVZ+zYKrN3ylKFmocAGb9TkwmYCjXzXIv8vevvvrq7uKLL+6e9axntQJgd7zjHbu73OUurdquas/jRL5yX7v63HggxAhAe0N/o4AJ2qq0jMzjjz++oVuK2n6r2uKJQnNDd7NOEIMlkS3PZBAVlWpVe3D39zWuZ0ppd17zaiWBSSSQ8qg8bCgT5j5ePN4384gouc3fvvObfezrmEnLoebeoujQKHgRLZKXld7x0ClR0WzXFvmkRHkJq5UEVkkCpX9XqTfmfy/rrn9z//1cIdXkLrjggkZDHmddoQAXx8hT3nvvvVveEX0g5/wJT3hC21K2OguaPu5xj2vzNzDiXtKwsKwLpNQ1fWLzt+/8lsYB9tGPfrR71ate1T3kIQ/p9tprr26PPfZotjq79vnPf373ohe9qMNaEAG6+93v3pxorvua17ymHescs7SNBkImMCAHAlUQgQdSgYR99tmnswrvhRde2BZuWmUQpHPjiQKChAElslknKJQ4A1NOUCXXz/IqbN+x43qm9LGFxqqVBCaRgCIIqAxKrVp0j2LjTZMYq+oPEGLzt+/8Zh/7OsaxzjFpC/UB31virXOLyPMkojbwOgJDBYQmlWztv0wJlP5dprSXf61117+5/36ukPkViFBZDr1sVAuVDXtKtWERGWDIoqX0AeeVzd++85t9gBG2JxDkXtLGBUJKbQMymE4iSz//8z/fwA6nHCe/5xB1Cph62MMe1oARneW682ASbDQQ0il48ZKxLGIZD+R97nOf7qyzzmqCt7joKjcgDc1NhzOEgR5gDhrPOkFVHW6Ve3D8e4tnZ1gRDKFlIeJ5vPjj31Xtuc4SMH8oynHGGWc0BYMClzLt1rNSHEEkxuZv30lYtY99KSXHOsekDqMURbDAHyoFEAQM8Swq93/AAQe0CFQBoXUeYZt976V/N7t/B59u3fWvOVrE5ed+7uc61DH2IbsXg0hxmnGZT+wLERo0OAu9szmVxjZP2/ztO7/Zx76O6bekpMgBRX+2+RvQ8lu/xWH3y7/8y42aTQ8BPa7lmegO55CLhIbXj0D1zzfJ/xsNhAxsFZDe8pa3tCiKCAo62TxWop1EyNPuaxI2eOX8oL2FQwkEAUNZJ6iqw00r4dU6Lp6dYTlDJqFf/dVfHTrxrNbT1N1stwQyf6jkY+4AbNDfKJNLLrmkJbZSYCI3Nn9LdvWbfezrGMc6x7iKNM/tfI6jLEXg0eLQ43C+lYfH977HPe5REaEIrD5XSgJ5f0r/rlS3LOxm1l3/AhLvec97topo3frWt25gSEU5FZNFfMzho1oiNAojKE7wwQ9+sLECwhzAEPCd3+xjnt/VAqp0iqiQHCA5oDZ/+85v/WadRc8AtKFS/97v/V4rrONa1kUS0KBP6CiOu34Eqn++Sf7faCBkTQxeTpQPSj1UMvxHYcRxB8YkAp3nvhlIBtCRRx65RSVJtQy0FYulVrWleUp9+88Vz9Rll13WvPM8Oya1hKIlFppQjI/BOvvbf+d1B6skgeTo8J6pBGTe4ARS7eezn/3szYBNDD+/2ce+jnGsc0w63uIRpEBf8YpXtEIJCibIOaLcwh2viNAqjZq6l0ig9G8ksbM+11X/JgIvx2f//fdv9i7H/6wLjrKTzeVhDvjbd6Ma0GU/9imwZPP3uHZ39JH9MbsWaeduNBBKqE1ZWJWRUOPwDZ/5zGe28NykHs5RHT/v39H2IHlAjrEQah+upModjAnPyJNh0FTbDAnEMyW3TYlh4W1gKMmJeLmS2CWhU9bVSgK7kgAFwhn0rne9qylD3O6UOR1GdTOP+G2wHCtF6hzO5ZzjthiS8tp48XgPKWvjdjCJtoDQuBKt/ZYpgdK/y5T26lxrXfWvhaqlUbz//e9v6/vEXpTrqQKoef2b3/zmxIIOoDGf2yYFMuQpYmTz9yR2d8CQ42yLsnM3GgihGEn6wjW02BRj0gJQOJNf/OIXJx4QyzogntyUFUzSWhnCy+qB1bgODwxPPM+5aimZ2AKEs+jYJBPLajxZ3cUyJEBpGUMiMsCGOdBCdOZEc+OwlnnTvo6ZFqhEgfYVZz+JdtrzD7v/+r4kMIsESv/OIr3NOXbd9G9Agxwczv4Uo/GpWM24RRM2pwfHf5KNBkI86us4IBgwvKa8/rsCQUWNGn+Ar/OeOLOqvZjAjjjiiK2JLdRIkxsubd/QXOdnrnufnwT6QEi+z9Of/vSWA2RR5mHNb/KE7OuYaYFKFHMfqBcQGib5+n4VJFD6dxV6YfvvYR31rzlXTts6BgC2s8c3EgjFE2kF84MPPnjLkz5riHAZHWUgWy/IvVt13YKDyW1ikKDJocuterW7Zchq068RQ3JYsq6KXm9729vmWj1l02W6k54PpxotTYGC+9///g1IK6eqWMHuyqn6zT725U10rHM4l3PO2goIzSrBOn5REij9uyjJrt9511X/mr+vuuqq7uSTT27VP9ctJWQ7RspGAiEecgldIipK76Va3KxJY4vuoLx4uJzKe1vh3eJSGcgKJiicMKzqxqLvr86/fAkYE8PKp6utL//tiiuu6G644Ybl31xdcaUlkGIxKg0qFmO8WJhOCXa0j2HNb/axr2Mc6xwKzzjnrK2A0KwSrOMXIYHSv4uQ6nqfcx31L3vhM5/5TPfqV7+6rdcWR/q6FAnbjhGzkUBIhQk8dyVa73e/+3Xf//3f3xLNDzzwwO5DH/rQylKJADiGBrCDEoefLy/IIlZZRVd0oE812Y6BU9dcngQkCfLEq9d//PHHby0MLPkd0JdDZOIritzy+mQdrpRkb+syWDtt2Mre/WfpAxXHOse8ItH9809Lvevfd/2/3hIwf1ndnhMTPW3Zeq7073qPn0Xd/brp31A75YZaf2ddAgGL6r9xzruRQAi1DOBJeVZgAqiwiB+v5rIn2HE6wj6UgLwmK/UqXeu+gTgLEZ5zzjlbC1HxUlTbORLQ3yZj1KSLL764O+aYY9qKy1nxWVK7SjHjlrXcOZLb2U9qPvn0pz/d+OI/8RM/MTUQcizOuXM556ytgNCsEtzM46P/JHsbI8t27OT6pX83c3xN+1Trpn9T7ENlTlGgFFlah9SQafto1uM2Egj1k32TXP6CF7yg+/znPz+rzOZ+vBcNOAufX2L8j/zIj7S1Y1S6sxK7KFEZunMX/VqdMIra4ro8PZngUJfkjgH5avUXUF6rbl3YzYqMiyBbzVtkXL6PtcfGzRGyr2Mc6xzO5ZyztgJCs0pw8443Z1n8XC7aeeed1xZwtGCj1ePNe0DRoua10r+bN54W8UTrpn8VUjr22GO7u9/97t0ee+zRHGEveclLOtVmzcHVbpLARgKhfvnXu9zlLt1hhx3WXXTRRW1ivenxV+MvIMhEbyXdgKDb3OY2rdy3NWQYuQDcsj1kqyGduotIQP+bjHl6jj766GakAkN3utOduoc//OHdy1/+8rZey6pGPPMc9bkcCXz7299ui5cqqGGBPc4V88ub3vSm3c6DjE/7ZC5yrHNYCNU5Z20FhGaV4GYdHyBy7bXXtoXDAW809mc961ktAo4lgTK+qHmt9O9mjadFPc266d8//uM/bpRmRZXkmt/1rndtdgNWyVe/+tVFiWktz7uRQEixAcj3gQ98YHfHO96xu9e97tWJBqEPKUu9ao1xcf3113dvfetbm8FiNWC8zh//8R/vXvziF7e1ZHaX3Lxqz1P3s1gJfOlLX+qUUJc3Jn8s1M8nPelJDSTx2qPSVdvZErCAHZrw5ZdfvrWO0GMf+9iWRDtqHSGJtvbNOkLO4VzOOWsrIDSrBDfr+FB5rJn2oAc9aKtKKv2HBoyWKYoJEDHgEiGalxRK/85LkjvjPOuif70rQA+nKRAEDAFF8j2BpGo3SWAjgRAuO2+SJN/b3e527fNlL3vZ3JJ9bxLffP4Ccj7wgQ90L3zhCxtoS5WPAw44oHvnO9/ZSt2qaV+tJEACKJJ49PLG5I+lGAgjwsSHYvK3f/u3JawdLgEeTGPl3e9+d/eABzygRZgf8YhHtAjz7haU9psotH1Rcx3rHPOi5hYQ2uEDs/f4xqmIz7ve9a5OZdc4Am9/+9s3A06EyFh89rOf3ZyFABEwNK9W+ndektwZ51kX/WueRYMTFFAoBz0OTQ5dDm2u2k0S2EgghFv8tKc9bavq2k/+5E82D/pnP/vZVor4psff3r9CCYDOX/nKV3ZocAwPCoA3zMRv0l8UJWB7n76uPq0EYuDKG5M/Zszc+ta3bqWOTzvttIogTivYDTsu1Y44WVAnRQ8Zmhbi/eQnP9kqc/HGp8Uzj6J00kkndT/90z/d3fnOd27HOofKhfOING4yEDKnk6N3VPUmm4iDKC3jXSU/z29jgFvzQ06rCB09sLvtC1/4QveHf/iHTSegUX/kIx+ZaFP172Mf+1jre1Tr3V3Lb+7JvblH95r7lof4N3/zN+25OOjynJ6Zrpq09de7usUtbtENbnEMJkKkmIHoEXm4L06facZl6d9Je6r2J4F10b/KaLN5sUfYwHmPHv3oR7c0DOO/2nclsJFASJnhgw46aKvqGoX+lre8pYXVV8lTHm4yBfXkJz95y6AVxkQJQJUT3qwBW6/roAQybhgzPPcB0MppP/jBD+7OPPPMZiQMHlN/7zwJmDcYiIzfQw89tLvnPe/ZIs4olBbcQxNmxKb523fWpTJ/AkHyihzrHM41j7loU4EQ2Xg3yVGEg5Fuk1ul0ASmgrlejp/NWk0oX5dcckmjf6Gs7G6zwDZHxy/+4i92j3rUo7qHPOQhE23AMCPIenTYB7u7lt9Q0tybe3SvuW9jAVAaBEmec9qIIVkp9JL1rgZBkL/7ESIVVQ855JA2zwFEctqm0euZR0v/Zgaoz3EkkHGz6vrXO8F+ZPuygfMePexhD+t+93d/txzsA529cUCIMtLJDMJMqKu6TgVvGm/b29/+9lYuOwMVFYASmjcFYKDf688NkACjVd7bqaee2gzc5AoB1R/+8Ien9pRugGjqEQYkIIog1xB1kpOFUjz99NMb3Q1F4stf/nLb/I0CF4ruMGDNIyq6wQBF0bQNVvfKpQEnyrgfCZHD+Y53vKMtBgyc/dRP/VR3wQUXdJ/73OeacZ3IA8MaqGB4zAOA5b6GfQbIuG+5UP1Ijrnac4qUoA96DuAGA4Hzjd7x3jHO0bys42FTaELVPXM6gGFZB5sFa5///Od3T3/60zsl8Bkou9voNCBAJT9FUqLfxv2MR1jxIOfZ3bX85p7cm3t0r7nvV7ziFd1v/MZvNBou3ZXnNHY++MEPbkWpRK3oMONP9KYfYQKARINQeVHALQ5tPIx6nsxz5AEYchiS+3XXXdcA6LiRKdc29lWqI49b3vKWnSUJ7n//+7c8utK/w96U+p4EVl3/JnIVanTeK3rAHGVuHWQE7ORe3SggFEVmUlwHIIQmcc0113RnnXVWy2OKosKH5oWjKAzmaiWBXUmgTykZzBViHEzrKd3Vteq79ZWAQiwM1uOOO64Z0ZYTwBkXVQCIXve617XN375TXIahfe9737sZmpQmj38aYMNIVFnuqKOOapu/+4YjEGQM9iMhQJA57/GPf3yj6zHsGdvG7Pve976tyMOy15OJp9d7hf7Vj+SYkz2nQhIisbj3clHRsEXQ6BwbQwMFkQPOpuoe55b1mLKwLfmTr2e/293u1mjckpl3twGm+i5l82PYjPsZR5uKpM6zu2v5TaEM9+Ye3at7tqGoAa8ANQCR55RLZoHnRKqMJdErYEU0qx9hEgWyZAQqJmp4inOMep7Mc+Thvsib80fVTOccNzKVZTa8F/rlDne4Q+ecBx98cPfe97639G9e+PrcpQRWXf9mPuOc8F7mvTI3cdSUfXlTt24UEArHneI2Oa5yxwNtlMCFF17YPeUpT+n22WefVtiBojn++ONbMtuyPKE3DYf6a50kwHNtMjbRMTp+8Ad/sPuBH/iBoQbsOj1b3ev8JBDgIupigb2999670YYZfRSkPDObv33H424fVQl5+/uAuk9tM2eJFqBO+S2No0cey5vf/OZGe0pEgecfo4601gAAIABJREFUCGJMG7Ny3NA7GaSqe2Y/URTRBGPcWJ+0mWMHozvyWgCcfmRHJOrjH/94iyp4l4Ax9MB+JMe6SkqKM9g5q1QlZUADDGQWfTPs07vpefsAZBjgCPDof1o3TNQiAGTcz2GArH/+/C86I6+sf7+ZZ4Y9Z74H2IAoYAlA7EeYUOFEYwAkYFIeAzCS48f9TIQI9Q9A1X+JUsqT6EeI4jDV78Ybw9B4RwV1DyeeeGID9aV/J33jdtb+66J/zcGcE3mf5r0cwib0+kYBIRMelJvqM+n4VQsFZiJGRQF6KB7V7YAhoAg4ApLsV60kMEwCMfR44k8++eQtw5Ixy/CoXKFhkttZ35sXgSHGIaPziU98YjP80ICMFcauzd++YxTaR1I6Y1FkZxCIjAuE4uiR72gB4L6BHYNaoQ9gSEUj0ajsJzdJQRDgZZqy3TFUEt2R1yK/pR/ZEY0C+rwzwKB7ZRz3gQOjXs4U4OJ+Lc1g3maIi1JE3wz79LzAX5+S1gcIAYLDPgMgQkkb97MP7IadP98HsPbvNwB22HPm+zAcRJ+MrT7gC6ATKUPXA4LQ03L8uJ+JEFlPzTmBffMhkPWZz3zmZhEiAOfGG29s+VoAuX4UJRP9AsgA8NK/O2uOnOZp10X/9oGQ6PQ8F8ieRnardsxGASGcRxQOioGXLBMpBYcutyoeHgpaFSGVmHhCQ3XgXUQZQZfjTa1WEhhHAjz2lP7hhx/eDLV4SBmSchf6HtFxzln7bJ4EAALjARXOAtOMfnSmRBRCbfKbfezrmH5L+VjK1Biz+RvQ8luaRF3RIB5281yuM+6n/BTlX/sUDvO4MQ2g8fhb3+iGG25ouSaKEiTCI0I1GN2R1yK/pR/ZGWbY96lkwA9jHhgCin7sx36s0cREEdDE8lyRYyhi+QS2gK5+kYI+ZSxFCYZ9hlJG102y9al+w86f7+UW/Nqv/drN7jfAMc9Fv9pCCRQBIg9yRaPrUwAnBZDR4+N+BnCaD9EY5VECRMaIMSP3S34HOwHg1c/OXfo3b259TiKBVde/ot3mnjieOB7MOSjLHGTVum7HACGJrKsSYTEZSxBVzcNEzHNmMqZAcPlRN2rdoHo9x5UAY5Cip/RNcvGQWjxNdbBxOfPjXq/2W08JiKoANhLEJfijgTF2E1FIsrvf7GPfXUViOJxEhdDWRGxs/vad39JQ0STJm3uvvPLKrevkeqM+gSD5Tf2kXiDImAa0lId1fsn2r3/969vaWjHU+9GdGOb9yE4MhL6hnYhGigugw6HFoccBU/Jezj777Fae1lye54kc++WtGSRAWr9sdb+IADnubgMM0QXJZZKtX/xhd9fw29e+9rXO4pH9+w3QzPP1i0TQYdYzC/DsF4WYlFLY75dR/4eCCLCaD0XcRLnQHRmt9Ctno5wiQC3nK/2bN7c+J5HAqutf76s5MQ4fgJ8zRrXEcrh/t6cLCE0y4ue0LwNDRZ1MxKGo4OmboMuDPydB75DTxDAEeoCfUHVMfPJCGDKM0molgUEJZNwkqjAuYE5EhjFuLrP5uz9vhZrmvKI2uc6wTwaqqAUnkQpgxi1DXJQH4El1tgA44MMaGSigihY89alPbbSoKPwYuPkMsBGRSGlwoKif/J+ITr/ctPwTzgYFE9DrAABUO/cGlOW5xpXjYF+s898cjMYEUAak6UeFCPQf+fTLhKfIhAIHokj77rtvd9vb3nYLkKS/5vWJlgcUi0yKUAKs+lIkEx299O86j77tv/fMo6uqf/tAyHwnMu7dNHdX20ERoVWixlES+NqoH3jJmajPOOOM5lldFQpfvSDrIYEYpqhMlDv6DjDEyOMJUgGJcVKtJDAogYybRBX6QGZw38G/Y/gCOiJGNn/3561w6J3XPrnOsE90YWAC8BE9QWm67LLLGtgRgUl1tlD60K/Q0lDU7nGPe7Q8EwBnWIQnyfsiEiqDJbLTLwediA6PKWPB4rPWDFEyW/QmUQXgjiHBIyzKn+caV46DMl33v/V1FpLFZhB9SnEKfRpwRH4pO04nqxT43Oc+t4GheQGf/nkCdBR/MGb0v/yvVEcs/bvuo2977z/z6Krq3z41jqNItLYcpDeNm42KCAWZU2S8epkQV6lYAoVBEQA9FLpJGBgCioAjIKlaSWAaCaAnWS8G3XLPPfdsgMgCivjADJBqJYHtkEAAEQP5m9/85lZOj4R01DnFPkTIUeFQ6BjHIi/KeZ9wwglbFeYC8DOvhwKqYIEoj4IPjN1EeBi73gXFD8y1KefsvOecc85WZAfYYRSYexnsOy2isx1jwjVFj+Q7ic4ocpB+XfSnyOAee+zRNn+X/t2uEbBZ111V/dsvlsAZgLpqvuOwqLZhEaEgc9QJii8TKmW4CnXTGQTukeK1MjYQxFtFcaPJMQYqVFmv5bQSyHoxz3zmM1sFrnjA/c/YrFYS2A4JhCIHZIisJKeHVxKtzVozQEoiPQCM3A7V40R6ABxRnlA+M6+nKAjOOy+/CmAqnSXCw7EkLyRrE2WBT44oUadEdsy5oheMAmBtJ0Z0tmNcyEHitBx3IdX0+6yfKYIBBPm79O929P7mXXNV9W8fCFUu3M3H3kZFhPJ4/Y6nWClF/PPtRMAMAvQPVCWKOxO6NSHwznHjUSyqlQSmkQBPtspd5557brfffvu1IhzAkHKyKkGh7tRK0tNIto4ZJoHBaE+oUEmyT1K9sQeMoLmpRpecnic96UkN/KBwGqeZDxPpuf3tb98piZwqbSm3zLElym8RU3lB8oPknahwxrBOhMd8D+RMuw7RsGeu7+cjARFBYNg4ULI9/b/sz9K/8+nPnX6WVdW/fXtYdUt5ciirbIJqGxYRSof2O35V6qanWhzQY/LNhA8UAUdAErBUrSQwjQRMaqo9veMd72ilfOP5RA0y8SkZWxPfNJKtY3YlgUS4zWuiPUmOB0aAkn71Np7I+973vls5PYxfUXEgiHc+82EiPcD8gQceuJXLk/VzRPctmq2wzLXXXtvKwIryKKwAhA1GeJK/5F6rrZYE0CKB4lDE0//L/iz9u1rjYl3vZlX1b98eLiB08xG2kRGhfnLYqtRN73OiU+ZTjXe0jX6y8c27q74pCQyXQKihkpBFgWJgMj554tGS5GhUKwmMI4FEeyh4URXz12C0BwgBRpTPVjI55ZKzEGe/eluAucUrLbLZz+kBlID2RHok0csfSZW2rJ/jPtxTzZfj9OJq7mNsWVAcLbIfEZwWCCWSmJwxhTOUPs9CvT6zeKp9S/+u5thY17taVf1bQGj0iNpIINQvF7gqddP7nGi8d8YCz6kSrOW1HD1ga4/hEjB+TMbocfKCYmDItUAdkowuYlStJDCOBJLbY8yItgAig9EetDSgBXgBYsxl8i0ULEheD2Mzhq2oD3CedXn6OT3AFFCVSI85UaW25PIAQCnTjeJZ8+U4vbh6+2SeUpVPblgcNhkn034mkpicMc4gY/S4447rXvCCF7RPlErjDxgv/bt6Y2Od7yjjetX0bwGh0aNqI4EQ5cmTaAE8q4FTzLyUlDgwsl2NMYEywgBgKNj87Tu/VSsJzEMCKEoWTLOm0F577dXd9a537Y466qi2gK/1TqqVBEggHsx+tEfkUMnoRHtQLd/whje0OXUw2tOv4hagg+6G9mbeldeDmoyO8ZCHPKQZvkceeWT3whe+8GY5Peh1aHYV6dns8QnEiuihNwImfeCTSI25y5bIoTLpgLaS6SpfGVOqw6YqIFCuxHpyxlQGtMguWjCAfeGFF7bFVeWcAUKlf6cbZ5k39KGtiot8rxxXRf8Oe8+8N5dffnmr3ok6XG1Dc4SUClZ+1SrglPWqTHj9SNWqALR6ETZLAsCOxSaBHyCIMQEUAUcm6WolARLIcgOS1gejPUAKsMJbD7xQnMAMUDMY7elXcePZF4WU8yEBHmjikAKigKlEe7IuTz+nR8U2UaiK9Gz2+DTuRPfkeil80QdCidQ87GEP62yJHL70pS/tzj777LaulPmNMce5maqA1nGx6O6nP/3ptkSFyoAKECmaIT/Sfo985CM71Llb3vKWK+MgXbfezryxUxcQHtVfq6J/gVTjXpEaEfu8Z5wH3ptaJuCmntzIiFB/wdJMrCeffHKjDS3bg5GQKWNDqD6eVIaFBGDfUwzLarkfhsfXv/71Vk8eeLRlLQ2TXJKOvVBVbWxZvTP7dVCZ0ODQ4dDiYqDi4+Pll6E5u4xX+Qx5v81zoj0UnveZggaEs26PZQYoxDe+8Y1tHkq0h9JEHwplaVhuT7+KW9bpkQAvX8i5zW0AjzHpXqoYzCqPnOXcGy803WL8WT8vVQAf/OAHdzY5s+jiHDe2VANUdl1FQJRJY9mCtsb2qHHlfTDu5E4CVimbvV36N1JOZCURWe/IOBud/Y1vfKNFT7fDcUDuZCnSZvO376p9VwKron8t9PyZz3ymBQXkySdizxlw9dVXj3xvdlJ/biQQMskqPmBtHqH0hNpxhnmFlo2EM+F96EMfalWQ4klVOY4RwiNrMlxWy/2goqC/XHLJJW3RTQtv+vu3fuu32ouSMrRVbWxZPTOf6yiIgN6kQIJCCZkAGap4+fq/wNB8ZL2KZ8n7zVgyt+CsoyHxojMs++v2mIcGoz0KGaAOxWDM+Onn9vSruBlbgLZqYK6Lhpy8Hs6U7TDaVrF/dvo9GQf0nWgNMJQqgAxqYEexI+wJoN2WyCFjm25n4KFQAlTA/qhxlevR/ah4gD3v+Hbp3/R/IiuJyCozP87mPRP1Ij+VZpc9n7s2+qFIsc3fvqv2XQmsiv5l31111VWdAMC+++675RCVPyxyOuq92Un9uZFAyCQpwdYK5Sa+AA+oGGUOSjaZLqvdeOONrVoXgDG4ftCyQ5QiQIwTKyADQJKTzzrrrJZY/4QnPKEt8mold9QY0YQsTMizy8C54YYbisO/rEEzw3X0s0lQ/w6GxL0LDOKK8M0g3BU4NEDHPGce44FMtMd7iiLEoHrPe97THC1yJU477bRGlUSRTBGNUCVSbUvkHI3YIqbyMYBo42dYbo9oDyMO2DGmlm2QrUBX1C1MIYGMX2AIsFn0+AGYgKh3vvOdLafIeEeNQ+G88sorl+4Yjcg8v/fHuyrq9Uu/9EtjbZwZnJYcGx/4wAeaowMNcFnrwnCsoF1jG9j87btq35XAquhfDB/j5NGPfnRjIcmXl0f3yle+shWhqf66SQIbCYQgXUYCD9OTn/zkLSqaJMnDDz+8GQeA0rJaAMill176PYapiZhHbFlUvVStO+OMMzrAh5Gjwg7D6F73uleLnskFSD6AaJpKUCIJp59+elulndyWdb/L6p9Nu07Gv5wMicQxeFFQeF8ZHvqw2npKIJ5k4AdNCA3y4osvbtGeY489tik+gIejhdfbejz3vve9t/LFQnnLuEi1Le+6fAweXvQ2EUX88mG5PYn2uJ+q4raeY2k77jrUTfNUojqLHD9sAe8KxyhHpGinPCEMETbCdukz7w9nAhBExwZYjPqkszl1zedoTjz8DF50J2Bo0Y0NoKR9gJu/l2lPLfr5Zj3/quhfFOjB6rHsvFNOOaWtWclRWu0mCWwkEMrjZcE2vGOVjAZzhQwSA3YZFCErrgvx480zNkLVw4VGA1h0iwdOJMzkRR7WWNhzzz0bJSYeX+DIxoPAkAKIJNsLq+JWP+95z2seNFSFZVL5Fi2fTTt/DI1+edrtpoJsmpzn/Tz6zZwkuuL9GpXbwyP82te+tkVvjz766PbOWiulD3Qy38hNvNvd7tYMrkR7GIYWLZW7KC+DUcZzjqKEXklhuhf3Va0ksI4SQFWy4C79y/l3+9vfvuXAcYoq3LFdTSRXNIhOBn6AM5FYDgxOUnp5cKO3fQ8AeX85N/bZZ592jMJQqKqec9F2Dcot+yk0Pn/7rtp3JbDd+lf/A/+c7P31BM8///w2Rmo9we8drRsNhHhHUIHQQkwaMQh4XwySZSl44X/5NqhmoiwBZAwPXOhFt3iQ3//+97fyoUCQheVMqMKk8fgCZbZQa1R7YmDZDxhSstSaDLzPPGzVVlcCJmM0qcEFC7c7OXh1pbUad0aBmZMYSJwNo3J7vI9KCTOiUh1wjz322MrtScQn881jHvOYlpyO9iraI0KNrit3kWfaXOS6KLCoRJSlaPaiDavVkH7dxaZKIDnD9K/oKN0HEKGAA/vb1fpACAjC1kBFZ59Y/2VwE73y/RVXXNExaE888cQW8eXgkL93yCGHNHCyaLuGPQH4uH+bv31X7SYJbKf+BYJE6NhpbLfkevrb2PE+iMRWu0kCGw2EMiD6uUKMB55Uih/H3qBdZFNsANiwqBsKGq6mNY6UlrXm0aJbIlIXXHBBi0jd4Q53aKH1Zz/72a2oRH/iTAQJx5TX+fjjj28giNwOPfTQ5jVmMKVlYuxXs8nv/c/sP2wize/988XT8u1vf3urChZAxmvtGR3n3vttXufLeYbdd/+6+X/a43L8tJ9ywTgBKH0J8MaeMWgsGpPVliOBjFu5goCFsUoZoap6j8xDHBA8rBK6s27PuLk9PNx3utOdGhgCirLOSjzLqcKVYiiodAxA1wd2hr03y5FOXaUksFgJbLf+HfZ09MhgREjhAVXYhuX6RC+bO0R+5IAec8wxTZfT6aMi/uYhTg36U7Q58w8bJBVjnZs9kDlhV7bRKH2W33elv11fgQclzeVHMdj1T2yQXM9nzpMiGeZNc+ioluMm1dOjzjvp78vWv+lfsuXgQoNDh0uxm+2mgk4qv2Xuv9FAyEsHDPVzhSQEo4ChgUgyNMHkBVyE8A1ME5xoCuoKD45QtsINJqBFN0Dhve99b6O+eDEkQ8sDsA4DINH3+JIFmZiwVKbJsfIPJGryWJlE0xIqTzUbkxu5D2vZf1hoPb/3zxdFcP3117ckUR4Pm/tDPRxWDXBe58t5ht33tM877LhZv0cNtfZGqKHGnjFoLBqT1ZYjgYxboMN7Y6yKEHuPzEEiw8AK2oucrqzbM25uj/3Q20RvRXs4eTgw4lkGsoCtGDsMhER7GBaLnv+WI+W6Sklg1xLYbv2767vqWjRlEiAUvcyb7/01jwBDdDmdztF16qmndpgfu3J00fMAB3moGpb5h0OWLfK2t72tOUYBFLoOoDA39NsoPZjf+/o712dzsT04ZuUXuV/XHARDriuHCu0uZdM9LzA0quX6k+rpUeed9Pdl69/Il2557nOf25ze2ABZPkO+EJZBzfc378mNBkJ5XC/Qy172su6hD31oW1wyA0N0Q7loA5YR3QcEOX7WT2DCCtd4+ELYijYcccQRbRJYRpIh48ukZ50QCyIOAiGgYtiLkYnXxOl4kxZPQ58T7Bn8hqP85je/uYXzeYOGtew/LNkyvyt/jsLj2vrQZKksugRuv3keLzdD0qSKNuBYIExfps3rfJHjeeed18YT6uB11103tAob+bkPXv+LLrqoHeNYzzMIJHOf8/50XYa2BHjg39gzBo1FY7LabBLw3jAWeFg5DZLTYy0uoCORnlBavDvoaHIV0HSybo+8QQorVLZ8hso7KreH0mPMcAqI9mSdlWGOgdmeuo4uCayXBLZb/w6T1qQRof55GPzKVov0Kp4wzK7IPMUhCpyweV71qle1iAHn5tOe9rRmj2CImEfocs5F53aNfos+HaW/nYuOBHTYWCJP5kK6kM5WbY7zBnizjIh7Q8sVNeeg4aQ1X7IrUPj9rS+H2WkKbsivpJNd1xyLQrhdDIhl6l+2hvmebUbnc6hhCiiEgzaNNsmW+vKXv9zvzvq/67odAYSEVnkVeMfxcBMqlPeiMIDVqg2geCTmPTJ4QNBchCbl54jKKDzAMDI5LbqlOo2JSZ5IqHEMMcmiwzw/7ssLxtBjuJsA3a/J0TFpKafp5ROqv/DCCxtoye/9z+yf6jj98pv5HcVHyWfnNDmbLAFIhRvkO/CCex7yZEweeeSRrSpawFCuO6/zZWI77LDDWkKrioTq9PO+mYD7zYRtTBl7Ocanicq5Ft0C3AKAjT1j0Fg0JqvNJgHvAOXDUUCexpmcxKzLlUhPkpwVJkBTNFblCnJKAKhAENATAJTPcXN7AB9GBsOBcYXuyxmwu/d6tievo0sC6yOB7da/wyQ1KxDifPFsv/mbv9kKGymcQCdxCnIapmWekgsoZ3T//fdvuVJyDNk/2DEHHHBAq6iHXicijbqvUiwKXr9Fn47S387zpCc9qeViAT/03lOe8pR2HXrdXCh3WyTLdS2uC7SwMYAh9N0XvehF3cMf/vBGzaff5UwNs9PoYLqYTiYHjm9z8HYxIJalf+OwBjbZsvpu7733biAIGMI04HxLKki/P+v/HQKE8oLwcqgUAwDd9ra3bQOFt9XA8ZJee+21LUl4mMdh2gFjAHohvewSNXlveCsUJcDJXXRjFAEHQt8mBwaWnBEviEkTCGRETVsO1+Skus1ee+3VQAr6nAjOsJb9Y/A51ndp+T2AleFIfin7zZgXCeLN8jzpT8alSJGo0aBc53U+ygXIM5lTOqPWpSJ3ckVTonQc49hRQDFymPVz0AGQJGFjkCyNyWrfK4FE8FBPgP/k8lDM3h9GB8PA+2KMiUCK7nmvGCPJ6aHQAfVxIz2MAf2jSiPjgUPBu1m5Pd/bP/VfSWAaCWy3/h12z7MAIXMVHafIiYgJh6CKkOYe85GodBpnDfBifjLPpDodPWqxTXRa9DqASDU6eoou97s5rs+wiD4dpb/DvHEegMZ1sBNSpRJYAYiiv/v61Fwr6mQ+ZDfRoXQpneqe+o0DSGVcND/n4iCl69h9y3A49+9nGfrXOACIMRA44LJmkEhQHGkYA+hywwBk/7534v87IiKUkClDVtUMA8PLl4UEoWcvm0RyYGjeAwZH3ySTEt48w0K2y0LogJ2JQ+gZeGB4eUk8t3UIVM8Rdp52gcRMjPMGQlYAB4aE/MlPRM99Cq+b2D2PST/9OUwR5P5mPR+6n3MlsmgMUSQ8ULuaaAcnQpFIm2OdY3fUwXlNRBQlg52i5AXU58YgWRqT1b5XAongkZsoDzBPgYjcikgyJHhJeVVVngRWFCQgWyCGkeHdMg45WMh7nEiPghbOjc6hUIKCCfjtldvzvf1T/5UEppHAduvfYfc8LRBKBADNiQMmlPthjkDsA3mDHCsi0KlOZ56hRwEONDh6zPwGuABDqoyh2/YZFtGno4BQHJnmRQWXULPYX6lSSY8rZBX93denIioWhc7v0Z8YFnRrv9HBnoFOdi6OKI5e8/g4uUX98836/6L1b8YBEMSelOrhuYEgti0gqGAC/UWfzdvBP6t8Vun4HQGEInBgwEttnQyTh0EjMpQXFkDBtxWpAZpCATPgZmnXXHNNi1xYzVp5WwaUXBoTGdC1rOZ6JiPPjnZmsmO0mfBEWBjM6GcW3UTxYQSaJHmXvNR9z1DuOxPjvIFQJtrktoikoACF8pP+JEveIkY+2VoMUpg4Lfc36/lyPcqBzMiP0shKzf1xYqyhBKDzycsaplhyn/P+jIdMtJPBbuwZg6JoxuSmtigIE7/ojn4jC2OYAqVgGQcxACQOp2Q844CyBrCTy4PGqg8PPvjgRuugYHg7M57y2c/pGTfS411jkLgnxtG8HTGb2s/1XCWBcSSwKvq3f699IBTDnYOUzlU9jnGfnEMGvQ0VloOL3kMbp/vQnkVX6CZ5NOa8NLYMQEDHc3xy5AAj/RzCPpAYBqz6+nQYoyPzYoo4iMx4luhv85x7Y3NgKnAeiWiwwRj3wIvnNQ+TDTsAtZxOpVsHmzmfbUAXm6fZNdZCBIz6zzl43CL/XpT+Tc4XfcbOEQkCgqKX2LRsW3YeW7cPZBf5zOt67h0FhBIZ8ZIZIPE0JFKgrDUwxHthQkGBQRcbVkxg3E7PxPF93/d9DXQxSkU2khg47nlm3c+LaZJVacazM4gZ8xZWVVFMeN2ExGA3qUqwQ0djFJp4h71Qeb5FASH3BaSJWImkpD/Sn5IxRfRudatbNTDkvu2blvvLxDzt+XI9dAEAyPOSH7obr2PuK9fFcRZl1N+iWo7ZFdUg+8/7kzJEIzCWUa4AfmOwr7jmfd3tPl8UBUVLWZAB4wGw501EZwMOQwmRLHzQQQc1EG0cGf+iPJQvY4CC0X9ZfwsIIsuMp3yGipD1esaN9DB4KGv3i8ZrnPVB9XbLtK5fElhXCWT+327925dfHwiZZ1C5REk4ZK6++urmzU/OIcq5zT6iOwCQRHjzkvlJPo5IMqqUOSQtgIPDxzkxKYCeAJLsF8M91DJgYldUu8gz815fn/R/Z1MAOxyxwE30pOtzNov60JGh0gFsbC8FE9zTlVde2X5H5/fMdOrg+k9xfNHBdLH7Zs8ABxxc/efM8y76c1H61/PQF2wcqQDAI+ATveRv9h0bl63rPgbHw6Kfex3Pv6OAUDrIwGDUJzJkAmLECClKLmP88LS85jWvacl7jCheFt7kDKpJDJX+xMAo3d16AbnPeX96GUyKch2UrwRygAbGn8RF/GAJkCYS3hkvlMIEJkMv3LAXK8+3KCDk3lSKG1aUINcfd2Ke9XxoBCeccMLWauDAM29XFFAMcZOwyVhkgFHtGMcuqzGqGdnGmjE3TD7Lup9R14lCM05Fcty/zTtHISaiQ4l7f03yvID9yE6qtCW6g1eOGiKRFDUx1dr0233uc5821r37kU8+EynGT4+xYa4YdBh4d0QhFUTo5/RUpGdUj9fvJYHFS6CvH7ZL//aftA+EOGJFbCxxwNiXVwP0JOeQLrbRJeyVeP45LeXbqASnYlqAxuD1zK0KEJg72TKcQz7ZAok0KZwUehw7gA2wq2U++vIcBoSSeqDYwQc+8IGbRbqjJzkHzaGZbxn2nFbRB32AA+y4B7pB2gN9QffSweZ087V5WX4R3bBdbR76lww8Hx3ICSyChj1Atyl+Rbb6if6ySXegh+hLow0iAAAgAElEQVT8YY7r7ZLHKl93RwKhDKxEhvAoTS4pN+iTp4XB/MQnPrGVdxSO5U02gQATzjFu608c2zUR56UyaeCMmiRMOCq1iBKhzfGSZ+LlGQeGACPgcFioNc+3KCBkomdUZuLryz3XjwE7bGLO77Oez2Q0WDSBx0sOSSgJ8dhQKsLzjGdRh2UVSYh85jER51zL+Ixi9H4BPUCcrR/RSc4OzyVl0I/spEpboju8iKigxjCeeaq1UZgShykQSjvjI5/xUFKu6BbocWhyoZCi0FFIFDnwRWGLuJpXQm2pSM8yRk5doyQwXAJ9/bBd+rd/h30gdOtb37otti6agU7GgUZ3xDFJv9qSexjPvwiS+QeoEXGh5wdbAIXoi0iKEvscoZZRkJuYSBM7SFTJfCn6Mmzh9748h+lb8yojnTOQncFmGry33Jd7N2eHmcO4F7XS7NOnvLme+49jFEige+lgujhyEWlCg96uNqv+jXzYa56P3NkQcqCwSzjkLKQd/QUEsVc58OmhOO236/nX6bo7EgilgwwUqJknweDibWEsmQBMSnmRTQy4qbzJvMu8G15s3n1gArffC8d4MxGpOMX74kWwxeMRA2tVJuK8aOTAU5RJMqH4E088sdGohN1vc5vbtIkZ3QcoYaCmZWJcFBDqT7S5bj5z/ci3v/+o33OefI7an2cG7/yss85qEYUkcTrOb+E242UD2Pl9WJJnrjvvz2ETMe+bMZnxOe0nwII6SqHvautHbuJ5xGu26naqr5GbrR/J4dWy9SM6ydnhrVTxcVRkJ55JoAf4oZwZF4wMxsawCI+CCPj0EogZDBQtEJaiIt53Mtgu6sW8x0udrySwiRLoz+eron/7QCjV3NDEzNH02ODGDjFvBQhx7oh0M3pFzIc1+kikCC1OcQV666STTmopANYOoudtylOjy5sXOYGG6fO+PIfpW88j4q4i3e6YEKPOR4+w0UTH3Fu/upx5GHUfa4X9htmieADbzMKz29XG1b/D1qFLLphnp3vYn6JhiQhGrwFDdBgQazxgopBJtfElsKOBUCJDcnWEG9GvhKOFp4GhDDTeEYl6DFqTEY8FL4cXnDdavoG8A+CIsccb7OWNR9tvJrcY6qsyERsmwBA5AG4mDZNzkjNNUKJEPOHAkGIEPDfDihEMmzj7w3HUxDfq9+0+XygGomiihskTSsSnP3EbM7sr+9l/nnn9P2wiNhaNyYzPaT9FxoCCrODd/+xHbuJ5NKED1Km+FmXfj+R4T2zkNxjRoQyNR4DGezoqshOHBqUMNCWHh2L1vg+L8FCuwBqFytvKaaJvU2Yefc+7sysqyrz6sM5TEigJzCaBvj5ZFf3bB0LsC3R1czNHlfse3OR6qvoaQ1g1tnPPPbfZLiLPw5q5C20ODd7c7zhgQVSB/ULH28yD5sSUYB6mz/vyHAaEzMvmW9cwhw5ro87XdyyK+AxWa2WvyGPGZAEUMVlE64EBc/R2tXH176h16MhQpMv4YGsECEevsUcxIzjo2Z/b/dzbJe9ZrrujgVAEF0rO9ddf3waTyiM4t4xcxhMQBAwlBJkBGMNKbXx5B9YGUuHEhMVIjkfb5GYQbxcQyvMx4HkfvKC4taOaUDpDUKiVLLyAXsR+EmImslETXyJQCYVHHsMm0mG/9+871x+2/6jfJz1f5Om8jGqTL+UqUmEdA6H88847r62ZwGgHrHnjJqVU9u9r0v+HTcRRuBmf035KgKU8s8J4/7MfuQnHHeAxnuJ5TL8N+4xDIhEd4AcIAoacw/PwBAJLQJO+oPCTv6MQAmpiHBeqEol6Amr6pSI8k46s2r8kMFoCcbKJmDJIo38Ytoy1RJFRnBiD6D+DOSuJIGMqyF/haESZMu9OstFf5obML6sKhEbdF4qwSnEiN+Y+lUjNaWhibJe+QybyJzc625zJdmHLmH85ojjo4sBSlACLQY4Set6sQGjY8f2RM0o/96nm9C2qMjvNGMEuYGNZFFZuk+UhMDa2O1I/rv4F4sZZhy5VSdlgdF5yw0TC6DBg2DNXm1wCBYT+f1TEJAIk8I57uUy4wIwwMlocAysTSd8wU5tf3gFjjLfGviYaE5vNoGXEbddEnImEN5v3JNzaUcPFccAQKhxvUrwR/bKUmchGhcIDIPpUwXUDQgF0QA9jXx+bnIEisjAp8U4ZC9aXUQVHxZ4+R3qU/Gf9fdhErJ+MyYzPaT95qYTks8J4/7MfuaEYbSKLlHGq3OS9GPbZdzygw6HFJWdHhEkhBLQAip1hwKsqf4eSRxUg/1BZk8OTqG1FeGYdaXV8SeDmEjDfcf4w0BQaiv6hXxnnMcDlqHLG9HNWEkG2zh3jX6SC9zsR5HE/gaA73OEO26Z/by6Z737TjwiZh8mBDWLu7jcAEq0Y7UtExzOhiXHCAop9wz/yFykgtzhyOXDJVgntlOx3L/oFGFKoQY7SMCATfZ/5epj+HnZ8/7lGnS92Q4oPAWn0LWc1Ow3oAX7oH2AIKAKO+sCwf91F/z+u/h13HTqOaDYFehydl9wwtEc6rN//i36+TTp/AaFd9GZePBO3l0yhBJMDr4oJxcAVXjZZeCklqTHuxjXsRk14u7ilqb6KwW6SM2EwEL08yk3zII16cSKHTEAAn2fEX5aImJcPhUg+BRlkzR9goD8R8QryBApbC1+PmkiH/d4XxqiJdNTvk54v+/eLJgA9FsXktQvXmsGOvsDTuew2bCKOXGf9TJUfoEbfB+jksx+56QOlYY4D0ZxBjnw/ooMGoFBCcnbkXjGscKp5TRldKVLAGKhWEigJTC+B6IF+RKcfyelHbgAdzghOCbqHYxFbIutzJYJMt4pCKMaDrpXIcT45FemVcfXrqHltWfp3lMQnBUIBNuY7zlkRIXMsXUPnADWAZ5r53zUs0uyZE1EgY/3SjyCIyuknUSbnNqeL3jO6d7cu36KBUJ4Hvc6YEenHPom+pXP9zRbjdDTO2Dfb3SbVv+kfDndAL+vQcTjSiRytHK8o+Gyaqgo3vx4uILQLWQZAmPglvxtwJnnGFu8KWg3ObSZwSWoMt3GpPsuaiKPAfv/3f7/dKxDnhVIaXOWY/kTYFwWgZJ9UP6OIeKFMpJScFx3FDn1OmJ2ha/I88MADm7epD7TkIJnM4tGKwho2kQ77vX+fo4DOqN8nPV/2Nzacm6IQYVFp8Oijj27yFSEULTFxkxWFtOw26UQceY/7aTwY8wwVCjOGSz77kZsYPvkcRiXtc+T7ER3UQ1Gd5OwAPgC2hGGFPwBuY4/h4F2uVhIoCUwvgeiBfkSnH8npR27oREYrnUPnmSewJXi1OUHiGIlDEdhRpj6OlHxytCjWw/Gisti489Ow/Zalf0dJfFIgZC4zpwE8nLOiZAALXUOnA5kcQWlAkX0DFBJRkF8jJ7qvn0UWRNZFVSTgh0Inn1OEJW2UPs3v+o8uEH0SbRrWsn/6q28P5DgFqThgs/5j9C2d6+9xc6ZyvkV/Tqp/0z8AD0ZJ1qFDQffc5ISFwgHL9qiqcPPrwQJCE8gyHhmeE1SbKAIeLgUE+snfXmjbdoXmA4RMeiY/CsjLltA4MASYmGBEd4TeTZ4pmiCK4eXjgWfUx/DldXHORHxMciY7lDmTXzw1DFYvLAOVEgXI8GHjcRo28Y07MabrRu0/6vecJ5/j7p+qgxSQCBcQyHtj87fvULVw3+277DZsIgZmjcmMz2k/U1XNhA0MhsqSz37kJlSYfHIsGDv94iKjAPqy5VjXKwlskgRiUJsfzPfoz+ZpeoDhTCegkdJxnBCozJwT/YhOHIGJ5IwbuYnnO0BnVOSYvkIBo1ewMQCZSbbtpqYPGzuTAqGcp5+7CzwCmACCEv4xkH3SPSkrHUOb/pUbqUgUJ5J+1+eiQWyY6C/7i0oAXFgkoeyN0o/5fd5AKEUT2FuxZXx6djKgxzzDuNT/yHNRn5Pq32Hr0AE++rwPXBd13zvxvAWEJuj1KBC5RIy1UAMMVKHjfjlgE4Jtu5I13S+wYsJTgtE6ATz48fQcddRRjbZFyaGyMUqt/pwy2jjbymyKJPEQCdniaUs2B5KcX+OFYvA7H89MPDVCuM7LiyGqJLkRvZBSG+Rsm8DIKS0T6TCglP3yOWr/Ub/nPPkcd3/AmLJhKIiS7bnnng0AAUH+ViRBBDELreb8y/ocNhEDQcak55xlS1U13kLj3/gZ3PqRG4p/cBPFody8T6Kv7tdWE/6yRkhdZydKIA497yLgA/SgTtMDKqD21+US2RHV6Ud0+pGccSM3McjHjRxzrNA5Eso50iYt7rLdxYqGjbFpgZD5keOSI9a6MQonAAIiOZdeeukWZUo/Dy40GgDKQcdxa4F0NDvRBlF6dDvOrQBa/amvAA0yBJ6AoVH6Mb/PGwjdeOONzakqj8l4zPN4djLoL7Q6TO7L+n5S/bu7deiiE2NzLesZdsp1CggtoaczMcSwX3ZonrGORoQrrOgBMJRQMi+EZEveexVjhNyzsCpPECUoGsSrhupEKfXXEeJNlEulghiKlPKWJlueQt8BRPJklCp2TpxmHNhwvtcVCGXoMCROOOGERvVQiGBY5Cz7L+tz2ETcl/ey7qeuUxIoCSxOAnF8MZpEejjrGK5ZDsF6M4wtEVk5IvQB2hSQYV6nB1RAHbUuVwzQRHRSLIiO2F31RvOOLZ7vcSPHQJo5lvMNG8MzTbIBeAzn7dK/w3p8WiCUfuZ8UjEsjkp6nbOT0xMDgxMUdf3aa69tepfOt1i8okccdyJDAJFjOO0wOjgz0ezkZepL+4vY+Z0jk3Orb8/09Ul+5xDcVY5RXx7ZP/3TP1/2z3Mbw0CP+7LeI8AGrIlQcuKuSiv9uyo9Mfo+CgiNltHMe/Rf9GUDoUQuhL9FeSSmKr1o4pRgmHwWYWZhcRVp8Ld5/EyGgMvxxx/fcnvQmUyyPPhp/vadyIjJiKcPGML5BqKUugSCgCxeLImdJjLXBxr6E19fXv3fc918jtp/1O85Tz4n3b9fNGFY5CznX9ZnTcTLknRdpySw/RIIFRqNxpyEvsxplQWyrREDhCR3hx5CNxOhRz+jB4AaiegcOqkyFgM1n/2ITpaPQI/dXfVG86otnu9Ej0dFjjnaPJPoRpZ/MLeNu/WrlC5b/w4bGdMCIecDCkTUMS5CXU+/ZBkHIMiYQI+TZwLg0OcALDAkvwh4pfcBRSAIzduyAiKD9geegA2gA5sDqBqlH/N7mCf9HKO+PLJ/xtfu9L3nNl6AMuDMfbEzRLPYFxy+q9JK/65KT4y+jwJCo2U08x79F327JuJwi3G9X/ziF7eIjcIGElpFcNzX4MZDZPKzr2NQ/4TkTa4mpLQoYJOQinDC6MpZOi/vItBjYuXZM/nz7PlfkqNwPs8Uj1+av32Xe+n/nv3yOWr/Ub/nPPmcdH8eV8ZGFnQD8Bgd8czlvMv+rIl42RKv65UEFieBzLMi/PIuzaVoyQCF/Ev5PKI9qE4M33POOaclXJuXJGAnryIGZz4T4UExYlRyjnGCcWJxjNEN5nIAioe/H9Exl4u6rGr1xlXRv/2RIWIHrKL80bP6Sd5uqrH29+//L/LnHMpgy+3RP/S5ohUqaAJC9HT0voqxJ554YnNM2jf9Kjok7xcTRNQwOdByi4AjwMQmz9h5R+nH/C5ShQFiWQMUzGEt+4+r7wE7Dl2ME+MViLfMCVYK0LwqrfTvqvTE6PsoIDRaRjPvsSoTsYnTpKi0pLwO3iSTKJoEANPnXsvroVzt6xjHOscgCCIc/1PSFLR8JPublER/eJcsoMqLY5IyOVDiJlyTPqXtd0mbaf72Xe6n/3v2y+eo/Uf9nvPkc9L9VTBDL2FsiAYxONAMGShyYLar1US8XZKv65YE5i+BGL7AD4PU/MmDj2J27LHHtvxNzivefQBGFS3J7mjI5iURA6AnACifiSQwiBmXGAMi+2jNABXdQEeIrJiL+xEdRq6Ee/PmKlZvXBX92x8R5mdglr6kK/UnPet7IGZUCzDGxuBklDtDn8bhGCAUvc9hJ18XYLVv+lWfAmTyg+nm5EBzbPqe/GzGnN9H6cf87vyKLDgHwDasZf9x9T1HKpBvvALvcpo4avtMlWHXW9b3pX+XJenZr1NAaHYZjjxDfyLGZ02pS+Bhu1omUnxyCqzPuzZ5mUQnbTmvCTWV6NDz0pKs6/wmwf51ovBzP/3fc558jtp/1O85Tz7H3Z+iMdmhCygqoaIRygGvHAA46r5zvUV9kj/lxxPI+6f8rMWAUR0B1WolgZLA6kgg86b5WFQgVdyS2wOIiLwoTMPLrqiBqDrwI4qTnMsAnER6gKDBdUnQ4XjfE+VJhIdxydPOIGWYc5Sh2NEN7sm86B7Xra2L/tXn9AkZ952Nu5N5X5+OOk/027T9muOH6edRv/efZdz97ccZqz9Ftoxh45oD0nd+X6XxWfq339Or+38BoSX0DaOT8Zk1EFTtQieg4HhYtqslkmMitfaKSXhwm3ZiyXkBBZVenH9wYve375w/a74MTmD+9lvuZdR9jNp/1O99+Y+7v/tTOZDnS3IpEISD/YxnPKMZEaPuu3/def+vDDpvo3KjvMR77LFH43ujLPAaVisJlARWRwLmC84T4AMISRW3wdwe7zFHGh0ij1MOh5wP73Z/nZ1EevrrkiiQoFCCKA8qXSI8PO1oR6hR5jUGNQBkjjZfmxcH5/HVkdzu72Rd9K/+TxRn90/0vb/uSp/u7jzRb9P2a44fpp9H/f69d9+1ceXZh50v+4fih8KHqgcEAUPJiVq18Vn6Nz23+p8FhJbQRxSNhNJQE3CvJThSdpRNtfWQAGNABE8IHohlRPCeKgQBADFKVOUzUaMZbHfj8VPiGsWF8WT8SS61ALCxV60kUBJYngRisCouk7XaBnN8QllCR/POpopbP7eHQw3oAX6AIPMOUAQcec+BJTkgifSg7cphVDgBTQmVDSULyDGnbXqjf+WsSt4XJSv9ux49Hn0rMioHDu3PWD7y/y8kqyCTvCqMB8yHVWt9/WvJkNK/q9ZL372fAkJL6BdGM6+FamxeBpVaVNgRzl2l5L4liGKtLwEE4U0rgsBLi6OMXgIEqbbE+BB92e5IX4TMq4viInnWmFMNiqGkXKp7rFYSKAksTwIMO+CDIyVrtQ3m+KiuaQ6R38NYTxW3ONBCeQOC0ODQ4dDi0OPQ5NDl0OZEktDoEumRw8hQRGUScXIPvO/uZx0jPJP2mLnOfB1AWfp3Ugluz/7RtwAQlgWAr8KhancAhbFvXSk6eTtTDIZJZ1f6Fy3VWCz9O0xq2/N9AaElyJ0iOu+889piZZL7KDjePgqLt6PaekiAF1cJcuW/LWQHAKX8KLqZBGP5QqviaWX8DFazkyuAJsOrZkxWKwmUBOYrgcGoD9oz0CMCw/Cx9IAFli+77LKttdp2leMjj0/pagAoa/SkghtjkAFoYWsFEhRKAKYk2ktmF2Ha7tzE+Up09rP1i9mU/p1dpss4Q/QtKqe8W4umqsgK0NK3mBjeq1WjxEU2ff3L9rNgLVuw9G+ktBqfBYSW0A/9F4JyQ5WTJ7S7spJLuLW6xAQSyMR8+eWXt6iKdRZEVyQZo8gBQYo/rIqntW8A4FSjEgBHxmS1kkBJYH4SAIIYZaHPqtCl8qZ5nhdYiWQGnCpXWattVzk+QFCqT2aNHlRqlDnecZRc1TzRhVSmBH5Q3QAv197u3MT5SXQ+Zyr9Ox85Lvss0bfy2RQjOuKIIxqzJvr2E5/4REstWNWoZunfZY+Y6a9XQGh62Y19JK+gKmKUIY/GuCsuj32B2nEpEkioXt6NkLx1FvSrCVmC8apEgiIM0SuUGWs5ocVZG0TUSklTYftqJYGSwOQSCOCRbM5YQzmTEygiY25AR1XO99WvfnVbewWY4cUGboCcYTk+wJHqjrzfyljTF1mjJxXcsoQBwFVtPAmU/h1PTqu217rp2778Sv/2JbK6/xcQWkLfUJRXX311C+Va12HcFZeXcGt1iQkkkORNSZC8jLjJlCwqilKZqxIJyiNloToLzlkp3ufLXvayVjGuctMipfosCUwmASBE1MU8wNhBdzv//PMbcFGgQN6gqI98BgndGABobkAQ2tuwHB/nUIHS2m6AD9pP1uhJBbdEe1bVCz6ZJJezd+nf5ch53ldZN33bf/7Sv32JrO7/BYSW0DeUGMWmEhCut0RXkSGrQaM3lFJbQifsoEvEY23M8SzHE62qlITqLIy3g0RSj1oSmFgC3qMYYyivHB+KHHivgBWgBXhRjES1yCTj94sa3OlOd2oJ3snzURBBYYR+jg9QBVwF7Ex8w3XALiVQ+neXYqkvFySB0r8LEuwCT1tAaIHCzalRpqwNgeOtIlA8ghLnrOeyqsl+uf/6XC8JxGPNWEOziSdaxThVpeQUoB1UKwmUBIZLAAgyd3tflLa+6KKLGrWUc8F7hcaGzoZyuvfee28tjxAgFIfXAQcc0D3rWc9qjjA64Morr2zn6+f4oNmh25U+GN4n0/xS+ncaqdUx00qg9O+0ktu+4woILUH2PHy8Uu9617ua8oyiRKHAJ2eUUrrVSgLzkIDxJpdAYrUqUxlv/lapsKpKzUPKdY5NkYDFJ5WTTpU3Sc5/8Ad/0FmI0/zMeWC5g8MOO6yVn0+EtZ/roxiJSpKcXfLyOLpE/c8+++xOgRWRJA6xev+WO3JK/y5X3jv9aqV/128EFBBaQp9F0VrQ7kEPetCWYYpHLqlWnolF9qqVBOYhgb4HNEAIJcf6IiZqXqtqJYGSQNdAkOIhqfKmvPwzn/nMBmSsYC+Sqkz+Xe5yl5bfmQhrIvtZz0dFxpe//OWtiIq1xkT7RX2uu+66Vkyln+dTsl+OBEr/LkfOdZXvSqD07/qNhAJCS+izcEYpRl5ClAlKFK3ida97XeOdo0VUKwnMQwLhxCu5Ky/BauqStSVyW2SxWklgp0kg+T4cTv3Ij3LUyt+nyps1euRwZp7Ouj6K3ABDP/qjP9rmbhF9iyon10clSVXjvvKVr3Tf+ta3iuK2IoOs9O+KdMQOuY3Sv+vX0QWEltRnJmPeQVSJKFkKVXlji+FZB6JaSWAeErBILwqcRXstHggESeSuFa3nId06x7pJIIYwCrLoez/yA8wANanyhvoWEIT+JgLkO5Eh9DhzNgcW8AREJddHJUnFDkJ1dt1qqyGB0r+r0Q874S5K/65fLxcQWmKfWfzOSuBWBt9rr706i+kdddRRWwnsS7yVutQGS0DJXeuPWL9E6d5awHeDO7se7XskwOBNpbeUuVeN7WMf+9jN1vfpR35CdUuVt/3226/bf//92yKoFkO1BpeCCRxXqseh08ktqrYeEij9ux79tO53Wfp3/XqwgNAS+0z1IYm3wE9WFAeKgCOTdLWSwDwkYP0R0Z+U8xUVEh0SJeKtqlYS2EQJJPKDloaeloWPzzjjjFbeur++Tz/y06/y9qpXvaq75JJL2jpBgBTwo4S26D0qMxAk/6Taekig9O969NO632Xp3/XrwQJCS+wzCpQ3EbUCLS6KV2LuJz7xiZVbkHOJoqlLzUECSQr+yEc+0j3lKU/ZWj+IZ5tRZ4ypWFWtJLDuEthd5Eehgje96U2tcIECBoqE/NAP/dBWkZpEfn74h3+40ZTlaiqFjSJnLk6VN+/Ln/3Zn7US2lXVc91HTNcAbOnf9e/HVX2C0r+r2jOj76uA0GgZzW0PXkReRfxyyjcK+ZGPfGR39dVXl8Kdm6R35ol4qNF1lOpV7SrVrdB7eLYZdd/+9rd3pnDqqTdGAqMiP0pXK2GtlLWS1kDQrW51qy0gFAfUQQcd1J1yyiltUVSLo374wx/uPvWpT21VeeM08L4AQa5Zbb0lUPp3vftv1e++9O+q99Dw+ysgNFw2c/8lL4okWzQNQEgybi10OXdR78gTDip6tDjV4hh9BbR35HDYiIeeJfID/ABBaMip9DYY+XnlK1/ZiR7JIZJPVGXlN2LIDH2I0r9DRVM/zEECqlEqxKL6pMIr5eieg1CXdIoCQksStMskdKrSEBoGI9XLYo0KC/ZZvfwv//Ivl3hHdalNkgDqpXWCTj/99A4dTrnfe9/73o3uU9TLTerpnfEss0Z+gCC0uKOPPrrRkc8///xuMPLzxS9+sVWR40D4zne+U+WuN3xYlf7d8A7e5sdTkfKqq67qTj755G7fffet1Idt7o9JLl9AaBJpzWnfT3/6090JJ5zQvJS3u93t2toUyrKqSCQZt1pJYFIJMBpVq0G7POSQQxolSDI4ityZZ57ZfeELX5j0lLV/SWApEgjgEZGxGGG/2tuonJ9EftDg0OHQ4tDjvAcKJVx88cUtN7MiP0vpzpW/SOnfle+itbtBcxjHiigzBob8Qxv6re/8Vm11JVBAaBv65k//9E+7Cy+8sCW077PPPs1zjx73/9o795gpyiuMazDVAEVaDcYbKpragnfxUkq1VMVL8Yp4QVAE1GoVK4qirVpLEUVRIaBgvQEiXkFEFAW8oSJRC2q9lTRWbag2jZE/sIbYvM3v6FmX4Vt3dnd2Z2b3mWTy7bc7887M887Me573nPMc5FnJIdIiBCpBwA3Jl19+OQwcODBst912oUOHDkawBwwYYDkQq1atqqRJbSsEGobA//73PwtL++yzzwLvxqjaW7mcH/f8IIwwZswYE0qAPKH0BvlBLcyV3uT5aVi3ZvZAGn8z2zW5PDEffylWjkgRBGizzTYzIRZyEHkX4S3Skl0ERIRS6BvC3wiDIxyOsDgv2EetihdeeEGx6in0SZ4PuW7dukB88ty5c0OvXr0KSeG77rpruOSSS8L8+fNNRCHP16hzzz8CbjBARghHwzhgphSBAoQK5syZE+64446Sam/lPD8zZswwEoV0NhLaECyOqUUIFCOg8bcYDX2uFYFS4y+CWITjMhnD+05LdhEQEUqhb6g87jUN8AS1a/lCUXIAACAASURBVNfOyBA1hYgxRa2IMBEtQiAOAmvWrLE6VFOmTLHaQRtttJGRIUKECLfEG0TIkRYhkCYC7vkh9A3jgJlSwkaQrCZnEhL/XWpv5Tw///jHPyysjverlN7S7OlsH1vjb7b7J29nV2r8RZiFnETed0z+aMkuAiJCKfQNgzSG6RNPPGFVy13mGKURFEdQHmGGX4sQiIMARVKZTR8xYoQlaSLA0bFjR4tPfvLJJyXLHgdEbZMYAnhheMd98cUXbXp+MA6YKSVshBj6H/3oR5ZY7ARenp/EukINtYGAxt82QNFXVSMQHX8JiyM87qSTTgrPP/+8InyqRrZxO4oINQ7rwpHcUCCGnXA4ktohQzvssEMYNGiQhYdINKEAlz6UQYDZ9d///vcmyb7lllsaCUK1hpyJZcuWaXa8DH76OVkE3NAk9I2cx6jnh5lSwkYgQBgMrp7pREien2T7Q62tj4DG3/Xx0H+1IRAdf3mnIZiAt/vtt99WiG5t8DZkbxGhhsDc9kFQ8kLR6+CDD7Z6Fz/4wQ8stAn5ReLmVdeibdz07dcIYHAS5kF+BbNPTqjxBkGEUM0i5+LZZ5+t6/rqq68GZsU8JEn90xoIuEFJv5N3weQNxIcJHrzd9913n6kYRj0/1E7Da4mgR5cuXUzhsC21N+X8tMZ9lNZVavxNC/nmOG6p8ddzcx9//HHl5uakq0WEUuyof/3rX4HQpcsuu8zqvVAAEzJEvDwJ7soVSrFzcnBoDFCMz+nTp4ef//zn5lUk38xD41CPI+cCCe16rsOHD7fQPM6Fc9LSGgi454d8R8RfyEdD+RIv96GHHhoOOOCANj0/3J94gbp16xbIiyyl9qacn9a4j9K6So2/aSHfHMctNf4y5jIBiTqhcnPz0dciQin2U6kHSblCKXZKTg7NbPzq1att5t2JtIcWNfovL/4777wzoNYlIpSTGyjmabrXh3wf8hYJd3v//ffDa6+9Fp577rmwYMGCcPfdd5sCJrXQEH9xz6R7fiA9fIeXkndb7969bbJn6NCh4Y9//GOQ5ydmZ2izRBHQ+JsonC3VWFvjr09k9+vXzwqbQ4KYLNKSfQREhFLso1KuVeUKpdgpOTg0L2EUuN588831QisbTYD8eNR5IQ8EyWK9+HNwA8U8Rb/PMBghQIi4oGo5YcIEU3pjwMfTCPmhDEDXrl2tJpqLv7jnh1ygo48+2iquIwYze/bssGjRokBIJTH08vzE7BBtligCGn8ThbNlGvP3YnT8LU5tYKKI+4tttWQfARGhDPSRYpUz0Ak5OgVyx/7zn//YbHz//v0LM/BOTBr119W9yEUiL0R1W3J0ExWdKoM1g7arvBEyhOT6ihUrrC4PXp9Zs2aZoiX5i5AaV3qD7Lj8/w9/+EMjQz/5yU/CvvvuGw466KDwq1/9ykgT5AkSBZnCk/nf//636Az0UQikh4DG3/Swz+ORS42/msDOY29+fc4iQhnoO8UqZ6ATcnQKn332mYUm3XzzzZaH4TPwjSJAfhxX97ryyiutLoxmv3J0ExWdKiSIMA5XeUPoAHnrkSNHhuOOO85C2cj3IayN8DbC3FzpjfA37j++wzNEeBwKhtOmTbMCv+QOMTtKOB3tE14HCfrqq6+KzkAfhUB6CGj8TQ/7PB651PirlIY89ubX5ywilIG+81hlYuX79OlTMDJ23333cO2114aXX37ZPAAZOFWdQooIuEueJEySMQcPHmyS6y6OgHQ2AgnId1LLwAlLvf7utNNOlujOfUt4k5ZsIuD3DTOZEB68iaj8RXN9XOXt0ksvNcXBffbZx8Rb/P7xUDfuLzxCSGBDkBBGQCABoQQEExYvXmxtS+wlm/eDzmp9BDT+ro+H/msbAX+PRsdfrxt0wgknWJSG3ntt45flb0WEMtA7HqtM8S0SiD3sZKuttjIjY8yYMeGdd97JwJnqFNJEgNAzjNmlS5cWSFD79u2tbtDOO+9sM/eEqVHDAGPVDdh6/UXymHuT86F6tpZsIuD3DTOZDOJMrDzyyCMb5Pq4yhvyrxBqYt5JAPb7By8Q7yaKoCKJPWnSJJPIxoNEaCTS2SgHIqVNvhj3KsfWIgSyjIDG3yz3TnbOzd+j0fHX6wYxac07UO+97PRZ3DMREYqLVAO2+9vf/hZuueUWC0fZdtttCyEnGLfIbGPI8JBpaU0EmM3nHrn99tstB8MNVIxWZqOuuOIKC2kaO3ZsIY8DY5aZfEKYfPuk/kotLhv3oc9U8m4o9vhwr6xcudKID16aOXPmmCeRAfuCCy4o3CPFYW7u9SHUbZdddgl4pelnBDEQRzj77LOtUCDiGAz6hLtxTIxJLUIgzwho/M1z79X/3KPjr6tiMiF4+eWXm41GmKWW/CEgIpShPsOlSiHV8ePHhz322KOQhMxMLfkgxNpDhrS0JgLUa7nnnnvCkCFDwo477mjEhpcxiem33XabJbczmHOfEObErD3hTW7oJkWAvB2pxWXjPvSZyqjHZ+LEiZbnc+qpp4a+ffuGXr16GamB3EByork+ToLw+iCIAFmCNCGNDvHBY43S27vvvmsE6PPPPzeBBUgQZEyLEMgzAhp/89x79T/36Pjr78vDDjssPPjgg/ZORHBGS/4QEBHKUJ8xo8vLmIrExxxzTCAZHWUu1EjIByH+HkNXhkeGOq0Bp7Ju3bqwZs0aC0G76KKLwv777x/IB+rQoYOFMGHoEprEdhjFHvM+b948M2YptgpxSpoQkSNCeBTk/b333gsff/yxhUXJQE7upoBgICxQ7O3BC0MIGtLTEJMXXnjB5KijHp/jjz/eiDCEB0EDJ7Be7wKvM4SISRfuKaSw8fqcddZZFjZH+BxhdNSHItSNe0uEJ7m+VUvZQkDjb7b6Iytn4+Mv71kmhzx3snPnziYgM2LEiPD666/r/ZiVDqviPESEqgCtXrv4zC4SsygvIT8LGSIPBDKEIhMqTApFqVcPZLNdSBBGL4Urmanffvvtw+abb27kBlUvarNARNxQ9Zj3Dz/80AjSlClTLPndc8/cIK71L2QMkuWGM7Ni3J8KmUrmPvKQN1TWEDjw/B4ILl4aipGSU8g98ctf/nIDjw9Eh9DIqKqg17uAKDGwX3/99WHmzJmW6EuRVFd5Q1CB40KsNfmSTJ+qlewioPE3u32T5pn5+MtE9CGHHFLInSQv97zzzjNv0EcffaRJojQ7qcZjiwjVCGA9dscAYTaWmQZqcniOB/KMN954oxkqeI40O1sP9LPTJv2LAcqM/L333hvOPfdcy9mgXguz/Mzgkw8E+SBBPbqwL6SZPBHCpPAqInGclKqchwZ4KBWDwlVXXVVIop8/f755Kl566SU7B1TKUJcjjpr7Fy/D2rVrWy651A0urh0MCGn797//bUpueHzfeustm2FctmxZgJiQH/jwww9bfs+4ceMCXkH6EvLTrVu3DTx9UY+P5/kQGsfkCuQZaWzuCd4zr7zySmAg1wRL9AnS/62IgMbfVuz1Da85Ov6SH8n71t+vkCLIEZOUkCUt+UVARCiDfYeBhGFCngcPm4c0oSJHrD/GL0Xg3AOQwUvQKSWAgBOZZ599Npx++umBGSg8QXgJITTnn3++EQ1I0JdffrnBEf1FTu0WPEbkeeBpTEpVzpNFuT8hZpxf9+7dC7LKvXv3Xs/wpqgmngcMew+n4z6HFBGW0iqLh+Bw7fQLYRXk30BKICeIXkAqCXlEoQ2v209/+tNCfg/EEw9xly5dLDzSJ0rcwxf1+ECekFun/aeeeio4MYV0YfSBP+8ceX1a5Q7UdX4XAhp/vwud1vktOv5CgghH9/cr3nTC5SBBhM9pyS8CIkIZ7DsMWEgOld2p6YFYAl4AL1xISAv5AOQKKDkvgx2YwClxD2CgEqaEFxBvoBMPiloOGzYsTJ8+PVb9Hr+fuF8gIbWqynXq1Cn06NHD8kr22muv4HLL3KOEcZYyzMl747wvu+wyEwSZOnVqoAaRG+iQAXKdli9fbuF17kHivPEiQfiQ6YbYEa5FyBgrJJBBq14TA44fx+BYflyePYwmBkLO65NPPrHVc3g4f4QFCBXkmrg2rhEywjVz7WAAQSTMjYGVZxsRCvCN5va0a9fO3gEMxIS9lSOe7vEhz4ewulYjnAk8hmqiBRHw513jbwt2/jeXnOT427oo5ufKRYQy2lc8iITLYDhhJGEYuSGE8cmsMUat5Boz2oE1nJYPxHj9IC14AfEGeigaIVEY0STMQwjiLLSJ4c79UquqHPfiNddcY94djG0vwAlh79q1a8GD6R4KDyXAsKcI649//GNL0EftjiR9D9miMCf5Lni/aNM9SOTELFy40EIAIRMYKCj4cC2skBDCuvC0QIaSXmiTtjkGx/LjQngI9SM0gvNatGiRrZ7DQ+4WqmtcC9fEtXGNhKdxzVw7GEBywRThAggOuVeQzWhuj0+EkKwLYSoXiugeH/J86olP0nirPSGQNgIaf9PugfSOX4/xN72r0ZHjICAiFAellLbxEBoKFg4YMMAUwjCGvIAXxSzJI1B9oZQ6qA6H5SVM3gjeBMLIMJ7dM4DHBbKBYQ0ZqNYDUq2qHAqGhOWRYwIxQSUObwf35+TJk9cz+IuT91Elg/xAgrgWPBqQIydK/teJHmSK63QPEvkso0ePDldffbURQ2ptIRwxe/ZsWx966KFAPhKeliVLlgRCCZNcaZO2OQbH8uPOmjXL5MwRo/jTn/5kYYeEHnoOz7HHHhsOP/xwu5a2CKJftxPFrbfe2oike3qY8CAEkrA4wuMIkyNczvN7JE5RhwdQTQqBbxDQ+Nt6t0Ijxt/WQzX7VywilOE+8ploZpwx/vr3729kaLPNNjMyhJGFEcYMv+oLZbgjY56az0RBgiAWp5xySiAMzj0DeA3wDuIlxFvI9tUshHhBhipVlYME4c3A2IeI4WFCKhsPyapVqwohYFE5Z+piEQ5HWBzkxuVHnQj4Xw/9I7zOBSEgT+TEQKTIP6J4nROEAw88MLBCFDwf6eCDDzYRCYQkklppk+vmGBzLjwtZg6hwPZwXwiasnsMD6WPSolTIoF+3x5wfddRRRnTc0+OED6EEvL8IJzDxgfgF3h7JlVdz92sfIRAPAY2/8XBqlq0aNf42C17NdB0iQjnoTcJxMC6vu+46M8Z8Rt3rC02bNs2UpjBuMXK15BMB99TgdUAqHYMaIYKOHTsaISK8Cs9EUrke3CuETMVVlYOUQGYIyyMkrNTiBoQX+EQYAUMeDxfhbng08CpBMFhdDABSAcmBVEB6PFQMUQDWUiFjTijS+uueLAgPZNE9OyTXkj8VJXDu4fHrdxU3CCPPsnt6PASQ0DbykqrxADK4O/HlvoG0QqAgwUnLY/txXAXPc6bK/YXUs4/eX6WeKH2fJgIaf9NEv3HHbvT427gr05HKISAiVA6hDPxOgjYv48WLF4ff/OY3hRl1ry904oknhvvvv7+inJEMXJZOIYIA5AIBBOryIHnsaoF4hZBSxzNAXlBSuTBuJMdVlYOgEI65dOlSux8jp1/412fWPK/GDXCuD28XxAvlMmS/CTuDJHFthLvhBSEvijA4Fw9AOZG1lIhAWgTIj0s/QVrxEEFu3LNzzjnnhFGjRm0Q0uceHq4dDFzFDQU5pNLd0+OiEOBIUVVwrWTxfiDUknadkEK0kWNHsCHJgqluSLgKnudMlfuLh5N9Ksl5qwQHbSsEakFA428t6OVn30aPv/lBpvnPVEQoR31c7kFlJhlDUzOrOerUEGy2n5l68m4Ii/LQMULitttuO/OekJdSL0PRDeZyqnI9e/a0Qp4Y1dxjtS5+XLwdeD0QQMALgvAABAFDHTEGwgFZXVZ64MCBli/joW8etoYAgYetJf2Xtj38zo/rf/v16xeoMUH4HwTOPTv02YIFC+xaikUeavHwxMUcbD3XjHwpvHgeoshkyvDhw41oIubAb2zDu4N92LeaBc8O7dx1110ByW5CKOOsbMs+7EsbWoRAFhHQ+JvFXqn9nBh/0hx/a78CtVArAiJCtSLYwP19xrVU6BSJ1JChehnMDbzUljqUv4QpRkrYlIc+QoLIC7v55pvNi1JPgovxW05VDo8MtYgwlpMKweS4eDvwenB9eEHwfiKVTegWuTDkyLF6oVEKgBIqiuHMiscBIQOI0wMPPFAQM3BRg1r/0iZte4FYP67/xUNGnt6bb75p9b3cs0OdoNWrV28g+12thyfuQwGmhNF5rtnQoUNDnz59TKmPXCtCDglz5C+eR35jG/LS2KeaEDzODYnu2267LZx55plGGj1nqtxfCCb7sC9taBECWURA428We6X2c8rC+Fv7VaiFWhAQEaoFvQbvGzVwosn0HkIlz1CDO6aKw9GXDKzMMuIJggR5oVMXw0Bq+YYbbjCjH3LQqMUHfEjPhRdeGA477DDzUpGjRLJ+tYZyvc7fc50IwYNAMbAludImbZNPlRQBrBcWtOsFIefOnRsGDRpkAg7kLiHBDvFxbxmf+Y7fICtsyz4QONqodIGs4rkjjJGcrs6dO1vOl3vOSv31HKlyuWeVno+2FwJJIqDxN0k0020ry+Nvusi05tFFhHLW7zzAHvISlVf2pHp5hrLdqT6g4rkjJ4hwODxBKIw5CYIUYVR6Pg5x6o1aPOmdXBXySshpoajrHXfcYbkmnH+WFs6Hc8bTUlzw1Auf1vqXNmmbY2Tt2tvqBzxpeLCQ8d57772N6CDfDbEmFM49ZHzmO36DDLEt+7AvbVS6RIkQwhd4M91zVuqv50gxKVANAav0PLW9EKgWAZ5/jb/VopeN/ehDJvOyOv5mA6XWOgsRoRz2tz/IXnCT5GzCXKLJ9fIMZa9ziwdSQhwRRvCcICdBzJz/4Q9/ME8Rno20FgYLCACqXxAiybSn1RPxj8v9RQjhlVdeGajlRAgcnmLqUSF0gTiCe8v4zHf8xjZsyz7sSxu0VckSJULcx5AfLUKgmRDQ+Jvf3szT+JtflPN35iJC+eszO2MeaMJ1ME6R3KXgalRuWZ6hbHWuD6Ceu4FENuFJnhPkhXIhQdSMwWDFm5HW4ufLOXCvIXEMMdKSTQS8vxBmQEkSYsP9RQgt3mMmTgixdA8Zn/mO39iGbdmHfWmj0hBIEaFs3hc6q+QR0PibPKb1btHfj3kZf+uNh9r/FgERoW+xyN0nQnWK68CgpkWsf9QzdN999xlhwrAmIZ4XgpbGIQDe9BVkotjwdOLq6nCEKREOR85Qmp6gxiGjIyWJAKITkBzkyKnNRF0jJL0vvfRSU+OLEhs3DFDqYxu2ZR/2pQ3aos24C/c2qnlIiBNqJ49QXOS0XR4R0Pibj17T+JuPfkrzLEWE0kS/xmP7A+51YIj9hwy5ge05Q4S+UJsF4waJ5KhBVONpaPcyCPiAifeOfvBQJCesrg6HMAI5QWl7gspcjn7OKAJ465DmJqeLnDPICKQEIY5SoW68Q/iNbZzAsC9t0FYlHkBU85AQd0IlIpTRG0WnlQgCGn8TgbHujWj8rTvEuT+AiFDuuzCYhwdy4y7fqGfIQ66o6YESGPK+kCdeELzMtdQHgXXr1oU1a9aYwAC5EggO9O3bN2yzzTYBL5ATVVSzSCpHErqR6nD1uWq1mhYCeHAg0UyIoAzHfXbMMceECRMm2Luh1Hnx3mAbtmUf9qWNSgk5RJ96SnilIWHkG82ZM8eSkt94443w6quv2krxVMgXUtnUDUKlkHeRFiGQRwQYQzX+Zq/nNP5mr0+yekYiQlntmQrPi5exq9lEPUOehL/ffvuFwYMHW5HKlStX5kYOuEIoMrM5JIi8iXvvvTcgPb3nnnuaXDEkqF27dpagPmLECDM66Q9IUCPV4TIDlE4kEQSiRIh8n2HDhlkOEIpspRZ+I0+IbdmnWiKEN5PispAglOiQXcfL+ec//zmMGjXKahVRrwiVRIrjUjeICQLVPSvVM/o+Lwho/M1eT2n8zV6fZPWMRISy2jNVnJfPTFGEcsqUKeGMM84wo4bQK4xvQrEImyMZmt/xQDArqyT4KsBuYxfHn7ytDz74wPDFCGSWvFu3bmHjjTcOm2yyieVh9OjRw8gRIUgyBNsAU19VjAD3Hc8zMuc9e/a0Z33kyJFh3rx5FhJbqkHCZdmGbXk/sC9t0BZtxl0gNYTD8a7Bs8TEC4IgrHymmOu+++4bKMzLdieffLLlFDFxo4mAuChru6wi4O9/jb/p9JDjr/E3HfzzfFQRoTz3XhvnzsvADSKS7seNGxf69+8fIEMY4ZChHXfcMRx00EHh17/+dbjnnnski9wGjtV8RXgEORUMhOB6wQUXWHFJSFCHDh0K+JNDgSjCU089VSBBCg2qBnHtU4wAOT1IYvPMowAH8Rg9enRZ8Q1C4HhXsC37sC9t0BZtxl2cCOHthAwhvIBniJXPvIN470CCkPvffvvtQ/fu3YNCQ+MirO2yjoDG3/R6SONvetjn/cgiQnnvwe84f0Kt8PrccsstRoaI3ccgwUjZfPPNw2677RbOPPNMC5VbtGiRGfAom2HM81LRUh4BBj7wwqvGDDrhQbfffnsYMmRIoT7Q9773PZPIplYLeRModD3//POmIldJMnr5s9EWrYwA+TbcV5BsPI487+QF8mxTC6rUwm9sw7bsw760QVu0GXeh1hReaNpghVDhXerVq5cRIAgPXqff/va3gdpnkKBOnToZQWKyZtKkSfYO4pnSIgTyjoDG3/r3oMbf+mPcCkcQEWriXibfhJcxak5z584Nv/vd7yxEBRLkZAjvUO/evc1whzBhzDgZamJoErk0XsIQRmbNSRQnnIgcLMJ/wNXrA/GXoqnkYMyYMcOkjDEwRTgT6QY18g0CkHGECBA+IB+tWiLEvrRBW7QZdyHE88477zRCBanCq8Qz8cgjj5j386WXXrIQON4x8+fPt8KtkCHeRXiLCKHDC8UzJTIUF3Vtl1UENP7Wt2c0/tYX31ZqXUSoBXqbJOrVq1eHp59+Olx88cXmlSAEhqTmTTfdtJA7xIzt+PHjw4IFCwIqT+xDmJ3U5da/ScADpSvwgWSCF6pvkKAddtihkAvUuXPnsPPOO1t4HGFy06dPt1A4GXrr46n/kkEAsRQUIadOnWpknHyfiy66KHaOENuyD0SeNmiLNuMukCYmBPAusUJq8JJGJ1aYAOC7xx9/3JTq/D1Ubf2iuOen7YRAGgho/E0WdY2/yeKp1kIQEWqBu8ALLWK4I2GLVwLDnHh9jBDPHdp2223DHnvsEU444YRw7bXXmqHiCdPKYfn2RoEEMfv9xBNPWP0VwnrI+4EEtW/fvoAnM+soZCGYQIgi6lzsq9nub7HUp+QQWLt2bfjoo4/s+d5///3tfoSc46Xhfi21uCfHiTz78o6gLdqMu3iIKKF2rHhKmUiJej6ZCOC7v/zlL+Y98vdQtfWL4p6fthMCaSCg8TdZ1DX+JounWhMRaql7wA0QDJ+HHnrIjBBqhxDLv9VWW1mIConOXnfokksusdAWkvqZ6WU/6g9FDZtmBhHS4vUI/vnPf1r9kyVLlph3h+KRFKH0ukAuU07uFdLBSGM/+OCDJqGNlKcWIVBPBAjF+fTTT8MDDzxQqCN07LHHhptuuqlsHSG2YVuvI0QbtFVPOXfqF6Vx3Hr2gdoWAqUQ0PhbCpnS32v8LY2NfkkOAXmEksMy8y15TC0zKh9//LHNyBKegvfn0EMP3cCg33XXXS3ZmdogJDgT788sbjTUJfMXXuUJOl6QP+oBURyS3IeTTjopEMbj4YXkOBQTyMsvv9wIEDkWzKpDgiBTWoRAPRFggoJnk5wcBAqY0Dj88MPt+X733XdLHprfeAewLfuwL23U+zlP67glgdAPQqCOCPh4ovE3HsiOl8bfeHhpq+oREBGqHrvc7+mGE7H8Y8aMCQMGDGjTwPdk/0GDBtkM7qOPPmqhXuTHfPjhh2YwEQed99wXD2HgxUttFQw1ErwRmqDuEh4ehCXcA+RqcITEEQbnIYVPPvmk7Z93PHJ/g7fYBRC+SigaIZtMbBDqSrgZeYHLly8PPKPc4774/b5s2bJw4YUXmpDK1ltvbfvShucH+val/no7yHDjOUZt7sUXX7Q8Q0Lkosd1A4dJFfKSqDFUzXFLnY++FwJ5QEDj7/q95O8Rjb/r46L/6o+AiFD9Mc7sEdxVj7LcX//6V0twpv5NNOQravAjAX3KKacEPB+zZs0ypTkMnryHzHlSKwYaZO+6666zoo/MkKP6hvz1lltuWQghjBJEF5nAIPziiy+UC5TZO785TwyCARmChPB8UqsHry4THNzP5Ahyj/vi9zthskcccYSREUg9+9IGbdFmucXbgTydc845RqSOPvroMGrUKBNoiR7X3zvUHSI0FxJUzXHLnZd+FwJZRsCfA42/X/eSv0c0/mb5rm3OcxMRas5+reiqfIYWhai///3vAY/GVVddVTIEjFAwvCLMNp911lnhhhtuCPfff38gdwZpXMLIPvjgA0uY/vzzz40UxDWqKjrxCjf264SwMdtNMjd5PxRApbI9hhm5EeQtEAp45JFHWqjQRhttFFg9B6jVQwYrhF2bNxgBPLVMUhC+iYw7Hhek8wl3Qyxl1apVtvKZ76hrxT2NcAo5bzz7TIz4wvOCIhw5gjzfrHzmO35zAwbPKcIhHJOcQz8u33NOeFmZJGDfV155xZ4zPKnUNmNbiBNqlVqEQCsh4OOSxl+Nv61032fpWkWEstQbKZ4LL2PICvHLGCsYQqVEAbxy/BZbbBG6detmSnMHHnigFUmEGFGM8e6777bZYAwbDCDaTTtUzGfgMOBQwyMkEENw4sSJVujx5JNPtqr3GGd4f8iXgPw4EZKIRIo3qA4dGwHCVfHU4p1BDhuiQU0hcv0gRBQuZeUz30GCvv/971vOG95gJkJ4B/jistioz1EwlZXPhMHxm4e0QHaoRXbiiScaFBloUQAAB2ZJREFUGfLjQo74ft68eWHhwoUBrzNqij/72c/Mw9q1a1erIYS6IiRJixBoNQQ0/o606Itf/OIXFmau8bfVnoB0r1dEKF38M310yAuGiXuICJkhLI4wMQwsvEIYUMhvb7zxxoV6RGyDsXTFFVeYwYVRRtgYuQPMQmMwQbTwxOCRQZ2K8ABigzkmM8zVrOxLGxQrJRwHuWqUqd566y0rDgnxWbx4sYkeIPxAgjgy4scff7xdk+f+uJw4xIfr3Hvvvc37hUgC4XIITEhWPNO3bkufnBOXW2+91SYnKFaKFxePD1LV5LKxumy1FzSljtjs2bPtmedZ8oWwV+oCIRQCoSpVqJVnGJn466+/3iYUunTpEjp06GCkCHI0cuTIMHr06DB8+HDzADGJgsIix6UOF/vShhYhIASCjYUafzX+6lmoPwIiQvXHOLdHiHqIMFTINcDAwqghDwCi4GTICQQGEOExGEyQCGrsIDLA7PPZZ59tOUhXX321eWLwyDzzzDOWk7BixQoLqWM2upqVcDzagHBBVmbOnBkmTJhghIwZ6FNPPTX07dvXVLF69uxpsuHkUZBUTr6Pq79xPVwXeRMkmk+ePNnkxmmXsD9IlieSM5OnRQhkCQEPZSOEjfsfjwxkiOLJkCHud1Y+8x2/sQ1khBBRSBDPvi9xiRBS2xCZ5557zp5xJkQgO4TJ8T7gmUJp0cnUIYccYsQI8sVx2beect1+PforBPKAgMZfjb95uE+b4RxFhJqhFxt0DRhYSOpCBh577DEzsgiFO+qoo4zoQHowcjB6IENOkDy0zAkGxIhcBDwxeGQgRWPHjrXwGULqMIyqWdmXEJxrrrnG8g2GDRtmydiE4PTo0aOg9ubn4yIQGIUQIorJUkyS3CCuC+8P14lHCSON6ye8TosQyAMC3LNMXhAGd9ppp5kHCOEPwlhZ+YxXiN/YppRHhmceUjV16lQTUsAzzGe+47fogpeXCQ5CZIcOHWoTJpAejsWECIp2tEEuEt5mZr2LPVDR9vS/EBACwcYfjb8af/UsJI+AiFDymDZti55jQ10c8n4IOyNPgLAZiAvGFOFwhMUxG+zeIice7jEid8BnpiEg3bt3txCZvfbay0Jm3FCr9C8J17QB6SHvYaeddjLyg9Jbp06dCh4fPx9XfXNCRkgPXiTyGEjmRj6b60TwgZnqtHOcmvbG0oXVBQH30CCOwP1MYWQIik8y8Jnv+I1tSnlkCFPFK0RIKyGurHzmO36LLigmQoaYMCEUdunSpeb15Vjz5883eW9U6QiPxfMb9UBF29P/QkAIBBt/mIzT+Ku7QQgki4CIULJ4tmRrGEOEiyGM8PTTT5tQArPBhMERDod6FQpzHo4GSSFkhiRppHMhRu3bty/kGjlRqfSvizhAeiA/5PxAhjge5MiJls+Ek5tAiB9iCRiFGIQUQPWwt5bsTF100yLgHl0PO2V2me/KLT4BwnMBWWLlc1wPqatisT3KWHq+yiGu34VAfAQ0/sbHSlsKgbYQEBFqCxV9VxECrhqF54RZYnJ1fDaYvBrCXx5++OGAQMG4ceMsbA1FK3J2CKsjVA5SFA2lq5QIuaw3pIdwOGqUEB6HLC/hcoTNET7nM+EUSyU3wUUbMAzXrl0bu35KRSBpYyGQMgJOaDCcWCslMuQs4GVi5XMlHlInQ+zHqty6lG8GHb5pEND42zRdqQtJCQERoZSAb4XDuvGD0UXNHpesJrSG8Jxp06aF8ePHWx2T008/3XIJyB9AQrN4JZ+orbV4Gz4jbjBw4ECT5iVEj0RxQt0QToCQuRhD3JnwVugjXaMQEAJCQAg0HwIaf5uvT3VF9UFARKg+uKrVbxDgZcyMFbPPhMRAiAijIwSN4q3vvfeehdQtX77clOMgLBQ2LV5RlaOmEdLXrHzmu+Jt+EyyN+Ftr7/+ugkckMOEhDbHQ1Lb5bnjzoSrE4WAEBACQkAI5BUBjb957TmddyMREBFqJNo6VlUI+MwWhIq1kpCcqg6onYSAEBACQkAICAELY2XM1firm6FZERARataebbLrggwVr012ebocISAEhIAQEAKZRKB47OWzFiHQTAiICDVTb+pahIAQEAJCQAgIASEgBISAEIiFgIhQLJi0kRAQAkJACAgBISAEhIAQEALNhICIUDP1pq5FCAgBISAEhIAQEAJCQAgIgVgIiAjFgkkbCQEhIASEgBAQAkJACAgBIdBMCIgINVNv6lqEgBAQAkJACAgBISAEhIAQiIWAiFAsmLSREBACQkAICAEhIASEgBAQAs2EgIhQM/WmrkUICAEhIASEgBAQAkJACAiBWAiICMWCSRsJASEgBISAEBACQkAICAEh0EwIiAg1U2/qWoSAEBACQkAICAEhIASEgBCIhYCIUCyYtJEQEAJCQAgIASEgBISAEBACzYSAiFAz9aauRQgIASEgBISAEBACQkAICIFYCIgIxYJJGwkBISAEhIAQEAJCQAgIASHQTAiICDVTb+pahIAQEAJCQAgIASEgBISAEIiFwP8BQ4KuTHZkQOoAAAAASUVORK5CYII="
    }
   },
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### What is a Markov Chain?\n",
    "\n",
    "A Markov Chain is defined as \n",
    "\n",
    "> a sequence of random variables where a parameter depends *only* on the preceding value. Such processes are \"memoryless\".  \n",
    " \n",
    "Mathematically, we have\n",
    "\n",
    "$$p(\\theta_{i+1}|\\theta_i, \\theta_{i-1}, \\theta_{i-2}, \\cdots) = p(\\theta_{i+1}|\\theta_i).$$\n",
    "\n",
    "For  equilibrium, or a stationary distribution of positions, it is necessary that the transition probability is symmetric:\n",
    "\n",
    "$$    p(\\theta_{i+1}|\\,\\theta_i) = p(\\theta_i |\\, \\theta_{i+1}). $$\n",
    "\n",
    "This is called the ***principle of detailed balance*** or reversibility condition (i.e. the probability of a jump between two points does not depend on the direction of the jump).\n",
    "\n",
    "**EXAMPLE** \n",
    "\n",
    "If you are an (*electromagnetic*) astronomer and you want to know how likely it is that tomorrow night will be clear given the weather tonight (clear or cloudy).  From past history, you know that:\n",
    "\n",
    "$$p({\\rm clear \\; tomorrow} \\, |\\,  {\\rm cloudy \\; today}) = 0.5,$$\n",
    "\n",
    "which means that\n",
    "\n",
    "$$p({\\rm cloudy \\; tomorrow} \\, |\\, {\\rm cloudy \\; today}) = 0.5.$$\n",
    "\n",
    "We also have\n",
    "\n",
    "$$p({\\rm cloudy \\; tomorrow} \\, |\\, {\\rm clear \\; today}) = 0.1,$$\n",
    "\n",
    "which means that\n",
    "\n",
    "$$p({\\rm clear \\; tomorrow} \\, |\\, {\\rm clear \\; today}) = 0.9.$$\n",
    "\n",
    "***NOTE:*** These numbers could be somewhat related to the mercurial weather we have in Nashville.\n",
    "\n",
    "- We can start with the sky conditions today and make predictions going forward more and more into the future.\n",
    "- This will look like a big decision tree. \n",
    "- After enough days, we'll reach equilibrium probabilities that have to do with the mean weather statistics (ignoring seasons) and we'll arrive at\n",
    "\n",
    "$$p({\\rm clear}) = 0.83,$$\n",
    "\n",
    "and \n",
    "\n",
    "$$p({\\rm cloudy}) = 0.17.$$\n",
    "\n",
    "You get the same answer for day $N$ as day $N+1$ and it doesn't matter whether it was clear or cloudy on the day that you started. The steps that we have taken in this process are a **MARKOV CHAIN**.\n",
    "\n",
    "Here is an illustration of this process from an article in [towarddatascience.com](https://towardsdatascience.com/introduction-to-markov-chains-50da3645a50d).\n",
    "\n",
    "![image.png](attachment:image.png)\n",
    "\n",
    "Sampling for $10,000$ days using this prescription gives a chain that looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's start off on a cloudy day\n",
    "# Cloud = 0\n",
    "# Clear = 1\n",
    "weather = [0]\n",
    "\n",
    "for ii in range(1,20000):\n",
    "    # implementing forecasts based on above probabilities\n",
    "    if weather[ii-1] == 1:\n",
    "        weather.append(np.random.choice(np.array([0, 1]),\n",
    "                                        p=np.array([0.1, 0.9])))\n",
    "    elif weather[ii-1] == 0:\n",
    "        weather.append(np.random.choice(np.array([0, 1]),\n",
    "                                        p=np.array([0.5, 0.5])))\n",
    "weather = np.array(weather)\n",
    "\n",
    "running_clear = np.cumsum(weather)[1:]/np.arange(weather.shape[0])[1:]\n",
    "# plot the running average\n",
    "plt.plot(running_clear)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The above figure is called a **trace-plot**, showing how our estimate of $p(\\mathrm{clear})$ evolves as the chain samples. \n",
    "- A histogram of the above plot reveals the distribution of $p(\\mathrm{clear})$. We can use this to determine the most likely value and an error on our estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(running_clear, bins=30, density=True, alpha=0.3);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- In MCMC the process must be **stationary** which basically means that the chain statistics look the same no matter which chunk you look at, e.g. first half, second half, or every other point, etc.  \n",
    "- Obviously that isn't going to be the case in the early steps of the chain. In our example above, after some time the process was stationary, but not in the first few days.\n",
    "- So, there is a **burn-in** phase that needs to be discarded. How one determines the number of early steps to discard as burn-in is tricky, but ***you should always start with a traceplot of your samples!!!***\n",
    "\n",
    "<font color='red'>Try the above example again but this time start with a clear day. Also experiment with chopping off different numbers of initial points as burn-in from `running_clear`.</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Markov Chain Monte Carlo (MCMC) Techniques <a class=\"anchor\" id=\"two\"></a>\n",
    "\n",
    "> The use of **Markov chains** to perform **Monte Carlo integration** is called ***Markov Chain Monte Carlo (MCMC)***.\n",
    "\n",
    "Given such a Markov chain of length $N$ that corresponds to draws of $p(\\theta)$, integrals can be estimated as\n",
    "\n",
    "$$ \\int g(\\theta) \\, p(\\theta) \\, d\\theta \\approx {1 \\over N} \\sum_{i=1}^Ng(\\theta_i).$$\n",
    "\n",
    "- To estimate the expectation value for $\\theta_1$ (i.e., $g(\\theta)=\\theta_1$), we simply take the mean value of all $\\theta_1$ in the chain.\n",
    "- To visualize the posterior pdf for parameter $\\theta_1$, marginalized over all other parameters, $\\theta_2, \\ldots, \\theta_k$, we can construct a histogram of all $\\theta_1$ values in the chain, and normalize its integral\n",
    "to 1. \n",
    "- To get a MAP (*maximum a posterori*) estimate for $\\theta_1$, we find the maximum of this marginalized pdf.\n",
    "\n",
    "In order for a Markov chain to reach a **stationary distribution** proportional to $p(\\theta)$, the probability of arriving at a point $\\theta_{i+1}$ must be proportional to $p(\\theta_{i+1})$,\n",
    "\n",
    "$$ p(\\theta_{i+1}) = \\int  T(\\theta_{i+1}|\\theta_i)  \\,   p(\\theta_i) \\,    d \\theta_i, $$\n",
    "\n",
    "where the transition probability $T(\\theta_{i+1}|\\theta_i)$ is called the **proposal distribution** (and it is assumed that we know how to compute $p(\\theta_i)$). The proposal distribution is exactly what it sounds like-- it proposes new positions to jump to in parameter space.\n",
    "\n",
    "This requirement will be satisfied when the transition probability satisfies **the principle of detailed balance**\n",
    "\n",
    "$$ T(\\theta_{i+1}|\\theta_i)  \\,  p(\\theta_i) = T(\\theta_i|\\theta_{i+1})  \\,  p(\\theta_{i+1}). $$\n",
    "\n",
    "Most of the difficulty in MCMC comes from either determining the burn-in or deciding how to step from one position to another.  In our circle example above, we drew points in a completely random manner. However, that's sub-optimal. ***We want random sampling that preferentially steps into regions of high probability density.***  \n",
    "\n",
    "The most commonly used algorithm for stepping from one position to another is the [Metropolis-Hastings](https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm) algorithm.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Metropolis-Hastings Algorithm <a class=\"anchor\" id=\"three\"></a>\n",
    "\n",
    "**The Metropolis-Hastings algorithm** adopts the following acceptance probability for newly proposed points to step to:\n",
    "\n",
    "$$ p_{\\rm acc}(\\theta_i,\\theta_{i+1}) = { p(\\theta_{i+1}) \\over p(\\theta_i) }, $$\n",
    "\n",
    "where the proposed point $\\theta_{i+1}$ is drawn from an *arbitrary* symmetric density distribution $T(\\theta_{i+1}\\,|\\,\\theta_i)$. Since it is symmetric, the **ratio of transition probabilities cancels out and detailed balance is ensured**. A Gaussian distribution centered on the current point $\\theta_i$ is often used for $T(\\theta_{i+1}|\\theta_i)$.\n",
    "\n",
    "This algorithm guarantees that the chain will reach an equilibrium, or stationary, distribution, and it will approximate a sample drawn from $p(\\theta)$! \n",
    "\n",
    "**In summary, the Metropolis-Hastings algorithm consists of these steps:**\n",
    "\n",
    "1. Given $\\theta_i$ and $T(\\theta_{i+1}|\\theta_i)$, draw a proposed value for $\\theta_{i+1}.$ \n",
    "2. Compute acceptance probability $p_{\\rm acc}(\\theta_i,\\theta_{i+1})$.\n",
    "3. Draw a random number between 0 and 1 from a uniform distribution; if it is smaller than $p_{\\rm acc}(\\theta_i,\\theta_{i+1})$, then accept $\\theta_{i+1}$. \n",
    "4. ***If $\\theta_{i+1}$ is accepted added it to the chain. If not, add $\\theta_{i}$ to the chain.***\n",
    "5. Use the chain (of $\\theta$ values) for inference; e.g. a histogram of $\\theta$ is an estimator of the posterior pdf for $\\theta$, $p(\\theta)$, and the expectation value for $\\theta$ can be computed from \n",
    "\n",
    "  $$ I = \\int g(\\theta) \\, p(\\theta) \\, d\\theta = \\frac{1}{N} \\sum_{i=1}^N g(\\theta_i).$$\n",
    "\n",
    "where N is the number of elements in the chain (e.g. the expectation value for $\\theta$ is simply the mean value of chain elements). \n",
    "\n",
    "***Step 3*** is very important. If we only accepted points of higher probability then we would only move up the posterior density to find its maximum. But we want to map out the full posterior density, so sometimes we must settle for less!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Corner Plots\n",
    "\n",
    "If you are using MCMC, then you probably have multiple parameters (otherwise, you'd be doing something easier).  So, it helps to display the parameters two at a time, marginalizing over the other parameters.  An example is given in Ivezic, Figure 5.24, which compares the model results for a single Gaussian fit to a double Gaussian fit. These are called `corner` or `triangle` plots, which we've discussed before.\n",
    "\n",
    "The most popular packages to make these are [corner](https://corner.readthedocs.io/en/latest/) and [chainconsumer](https://samreay.github.io/ChainConsumer/).\n",
    "\n",
    "![Ivezic, Figure 5.24](http://www.astroml.org/_images/fig_model_comparison_mcmc_1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Worked Example <a class=\"anchor\" id=\"four\"></a>\n",
    "\n",
    "We'll end by going through the example given at\n",
    "[http://twiecki.github.io/blog/2015/11/10/mcmc-sampling/](http://twiecki.github.io/blog/2015/11/10/mcmc-sampling/).\n",
    "\n",
    "Let's generate some data points and plot them.  We'll try a normal distribution, centered at 0 with 100 data points, using [numpy.random.randn](https://numpy.org/doc/stable/reference/random/generated/numpy.random.randn.html). Our goal is to estimate $\\mu$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import scipy as sp\n",
    "import pandas as pd\n",
    "from scipy.stats import norm\n",
    "\n",
    "np.random.seed(123)\n",
    "\n",
    "data = np.random.randn(100)\n",
    "plt.figure(figsize=(4,4))\n",
    "plt.hist(data, bins=10, density=True, alpha=0.3)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('data density')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Now we have to pick a model.\n",
    "\n",
    "- We'll assume a Gaussian likelihood distribution: $\\mathscr{N}(\\mu,\\sigma=1)$ (i.e., with $\\sigma=1$).  \n",
    "- We'll also assume a Gaussian distribution for the prior on $\\mu$: $\\mathscr{N}(0,1)$.\n",
    "\n",
    "We can use that to write an analytic function for our posterior distribution as follows (*recalling that the product of two Gaussians is also a Gaussian*):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Execute this cell\n",
    "def calc_posterior_analytical(data, xgrid, mu_0, sigma_0):\n",
    "    sigma = 1.\n",
    "    n = len(data)\n",
    "    \n",
    "    # Error-weighted sum of prior and average of data\n",
    "    # Recall the product of 2 Gaussians\n",
    "    \n",
    "    # This is the Gaussian for the prior: mu_0 / sigma_0**2\n",
    "    # This is the Gaussian for the likelihood: data.sum() / sigma**2\n",
    "    # We are deriving the mean from the sum of the data values\n",
    "    \n",
    "    mu_post = (mu_0 / sigma_0**2 + data.sum() / sigma**2) / (1. / sigma_0**2 + n / sigma**2)\n",
    "    var_post = (1. / sigma_0**2 + n / sigma**2)**-1\n",
    "    \n",
    "    return norm.pdf(xgrid, mu_post, np.sqrt(var_post))\n",
    "\n",
    "plt.figure(figsize=(8,4))\n",
    "xgrid = np.linspace(-1, 1, 500)\n",
    "posterior_analytical = calc_posterior_analytical(data, xgrid, 0., 1.)\n",
    "plt.plot(xgrid, posterior_analytical)\n",
    "plt.xlabel('mu')\n",
    "plt.ylabel('post probability density')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's try this instead with MCMC. \n",
    "\n",
    "- We'll start with a trial $\\mu_0 = -1$ and evaluate.\n",
    "- Then we'll jump to a new position using the Metropolis-Hastings algorithm mentioned above.  \n",
    "- Next we evaluate whether that jump was \"good\" or not -- by seeing if the value of the likelihood $\\times$ prior increases.  \n",
    "\n",
    "*We always work with the un-normalized posterior (i.e. likelihood $\\times$ prior) in these MCMC algorithms (or equivalently the un-normalized log-posterior which equals log-likelihood + log-prior).*\n",
    "\n",
    "<font color='red'>Read through the first function below to understand how it is implementing the Metropolis-Hastings algorithm. But don't worry about the second function, which is just setting up plots.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "## Execute this cell\n",
    "## See https://github.com/twiecki/WhileMyMCMCGentlySamples/blob/master/content/downloads/notebooks/MCMC-sampling-for-dummies.ipynb\n",
    "def sampler(data, samples=4, mu_init=.5, proposal_width=0.5, \n",
    "            plot=False, mu_prior_mu=0, mu_prior_sd=1.):\n",
    "    mu_current = mu_init\n",
    "    posterior = [mu_current]\n",
    "    for i in range(samples):\n",
    "        # suggest new position\n",
    "        # changing the width of this distribution changes how big the jump is between trials.\n",
    "        mu_proposal = norm(mu_current, proposal_width).rvs()\n",
    "\n",
    "        # Compute likelihood by multiplying probabilities of each data point\n",
    "        likelihood_current = norm(mu_current, 1).pdf(data).prod()\n",
    "        likelihood_proposal = norm(mu_proposal, 1).pdf(data).prod()\n",
    "        \n",
    "        # Compute prior probability of current and proposed mu        \n",
    "        prior_current = norm(mu_prior_mu, mu_prior_sd).pdf(mu_current)\n",
    "        prior_proposal = norm(mu_prior_mu, mu_prior_sd).pdf(mu_proposal)\n",
    "        \n",
    "        p_current = likelihood_current * prior_current\n",
    "        p_proposal = likelihood_proposal * prior_proposal\n",
    "        \n",
    "        # Accept proposal?\n",
    "        p_accept = p_proposal / p_current\n",
    "        \n",
    "        # Usually would include prior probability, which we neglect here for simplicity\n",
    "        # Accept if p_accept > 1 (since random limited to [0,1]) and if p_accept > rand as well\n",
    "        accept = np.random.rand() < p_accept\n",
    "        \n",
    "        if plot:\n",
    "            plot_proposal(mu_current, mu_proposal, mu_prior_mu, \n",
    "                          mu_prior_sd, data, accept, posterior, i)\n",
    "        \n",
    "        if accept:\n",
    "            # Update position\n",
    "            mu_current = mu_proposal\n",
    "        \n",
    "        posterior.append(mu_current)\n",
    "        \n",
    "    return posterior\n",
    "\n",
    "## Function to display\n",
    "def plot_proposal(mu_current, mu_proposal, mu_prior_mu, \n",
    "                  mu_prior_sd, data, accepted, trace, i):\n",
    "    from copy import copy\n",
    "    trace = copy(trace)\n",
    "    fig, (ax1, ax2, ax3, ax4) = plt.subplots(ncols=4, figsize=(16, 4))\n",
    "    fig.suptitle('Iteration %i' % (i + 1))\n",
    "    xgrid = np.linspace(-3, 3, 5000)\n",
    "    color = 'g' if accepted else 'r'\n",
    "        \n",
    "    # Plot prior\n",
    "    prior_current = norm(mu_prior_mu, mu_prior_sd).pdf(mu_current)\n",
    "    prior_proposal = norm(mu_prior_mu, mu_prior_sd).pdf(mu_proposal)\n",
    "    prior = norm.pdf(xgrid,mu_prior_mu, mu_prior_sd)\n",
    "    ax1.plot(xgrid, prior)\n",
    "    ax1.plot([mu_current] * 2, [0, prior_current], marker='o', color='b')\n",
    "    ax1.plot([mu_proposal] * 2, [0, prior_proposal], marker='o', color=color)\n",
    "    ax1.annotate(\"\", xy=(mu_proposal, 0.2), xytext=(mu_current, 0.2),\n",
    "                 arrowprops=dict(arrowstyle=\"->\", lw=2.))\n",
    "    ax1.set(ylabel='Probability Density', \n",
    "            title='current: prior(mu=%.2f) = %.2f\\nproposal: prior(mu=%.2f) = %.2f' % (mu_current, \n",
    "                                                                                       prior_current, \n",
    "                                                                                       mu_proposal, \n",
    "                                                                                       prior_proposal))\n",
    "    \n",
    "    # Likelihood\n",
    "    likelihood_current = norm(mu_current, 1).pdf(data).prod()\n",
    "    likelihood_proposal = norm(mu_proposal, 1).pdf(data).prod()\n",
    "    y = norm.pdf(xgrid,loc=mu_proposal, scale=1)\n",
    "    \n",
    "    ax2.hist(data,alpha=0.5,density='True')\n",
    "    ax2.plot(xgrid, y, color=color)\n",
    "    ax2.axvline(mu_current, color='b', linestyle='--', label='mu_current')\n",
    "    ax2.axvline(mu_proposal, color=color, linestyle='--', label='mu_proposal')\n",
    "    \n",
    "    ax2.annotate(\"\", xy=(mu_proposal, 0.2), xytext=(mu_current, 0.2),\n",
    "                 arrowprops=dict(arrowstyle=\"->\", lw=2.))\n",
    "    ax2.set(title='likelihood(mu=%.2f) = %.2f\\nlikelihood(mu=%.2f) = %.2f' % (mu_current, \n",
    "                                                                              1e14*likelihood_current, \n",
    "                                                                              mu_proposal, \n",
    "                                                                              1e14*likelihood_proposal))\n",
    "    \n",
    "    # Posterior\n",
    "    posterior_analytical = calc_posterior_analytical(data, xgrid, \n",
    "                                                     mu_prior_mu, \n",
    "                                                     mu_prior_sd)\n",
    "    ax3.plot(xgrid, posterior_analytical)\n",
    "    posterior_current = calc_posterior_analytical(data, mu_current, \n",
    "                                                  mu_prior_mu, mu_prior_sd)\n",
    "    posterior_proposal = calc_posterior_analytical(data, mu_proposal, \n",
    "                                                   mu_prior_mu, mu_prior_sd)\n",
    "    ax3.plot([mu_current] * 2, [0, posterior_current], marker='o', color='b')\n",
    "    ax3.plot([mu_proposal] * 2, [0, posterior_proposal], marker='o', color=color)\n",
    "    ax3.annotate(\"\", xy=(mu_proposal, 0.2), xytext=(mu_current, 0.2),\n",
    "                 arrowprops=dict(arrowstyle=\"->\", lw=2.))\n",
    "    \n",
    "    ax3.set(title='posterior(mu=%.2f) = %.5f\\nposterior(mu=%.2f) = %.5f' % (mu_current, \n",
    "                                                                            posterior_current, \n",
    "                                                                            mu_proposal, \n",
    "                                                                            posterior_proposal)) \n",
    "    if accepted:\n",
    "        trace.append(mu_proposal)\n",
    "    else:\n",
    "        trace.append(mu_current)\n",
    "    ax4.plot(trace)\n",
    "    ax4.set(xlabel='iteration', ylabel='mu', title='trace')\n",
    "    plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Each row below is a single iteration through our Metropolis-Hastings sampler.\n",
    "\n",
    "- **1st column = prior distribution**. The distribution is static and we are only plugging in our $\\mu$ proposals. The vertical lines represent the current $\\mu$ in blue and the proposed $\\mu$ in either red or green (rejected or accepted, respectively).\n",
    "\n",
    "\n",
    "- **2nd column = likelihood**.  The likelihood function changes in response to the proposed $\\mu$. The blue histogram is the data. The solid line in green or red is the likelihood with the currently proposed mu. The more overlap there is between likelihood and data, the better the model explains the data and the higher the resulting probability will be. The dashed line of the same color is the proposed mu and the dashed blue line is the current mu.\n",
    "\n",
    "\n",
    "- **3rd column = (normalized) posterior distribution**.\n",
    "\n",
    "\n",
    "- **4th column = trace** (i.e. the posterior samples of $\\mu$ we are generating) where we store each sample irrespective of whether it was accepted or rejected (in which case the line just stays constant).\n",
    "\n",
    "Note that we always move to relatively more likely $\\mu$ values (in terms of their posterior density), but only sometimes to relatively less likely $\\mu$ values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "np.random.seed(123)\n",
    "sampler(data, samples=8, mu_init=-1., plot=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "What happens when we do this lots of times? You might need to wait for a minute. Consider this a lesson in that running MCMCs can sometimes feel like watching paint dry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "posterior = sampler(data, samples=15000, mu_init=-1.)\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(posterior)\n",
    "_ = ax.set(xlabel='sample', ylabel='mu');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Making a histogram of these results is our estimated posterior probability distribution. Let's remove the first $500$ points as burn-in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "ax = plt.subplot()\n",
    "\n",
    "burn = 500\n",
    "ax.hist(posterior[burn:], bins=30, alpha=0.5, \n",
    "        density='True', label='estimated posterior')\n",
    "\n",
    "xplot = np.linspace(-.5, .5, 500)\n",
    "post = calc_posterior_analytical(data, xplot, 0, 1)\n",
    "ax.plot(xplot, post, 'g', label='analytic posterior')\n",
    "_ = ax.set(xlabel='mu', ylabel='belief');\n",
    "\n",
    "ax.legend(fontsize=10);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute the expected value of $\\mu$, credible regions, and the MAP value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# expected value\n",
    "ex_mu = np.mean(posterior[burn:])\n",
    "\n",
    "# 68% equal-tailed credible region\n",
    "lower, upper = np.percentile(posterior[burn:],q=16), np.percentile(posterior[burn:],q=84)   \n",
    "\n",
    "# MAP value\n",
    "# many samplers output the likelihood x prior, \n",
    "# so it's usually than this\n",
    "hist, bins = np.histogram(posterior[burn:], bins=30)\n",
    "bin_mid = 0.5 * (bins[1:] + bins[:-1])\n",
    "idx = np.argsort(hist)[-1]\n",
    "map_theta = bin_mid[idx]\n",
    "\n",
    "print('Expected value of mu = ' + str(ex_mu))\n",
    "print('68% credible region = ' + str(lower) + ' to ' + str(upper))\n",
    "print('MAP value = ' + str(map_theta))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Our algorithm for deciding where to move next used a normal distribution where the mean was the current value and we assumed a width of 0.5. <font color='red'>Make this much smaller and much larger and see what it does to the sampling efficiency through the behavior of the traceplot.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "posterior = sampler(data, samples=15000, mu_init=-1., proposal_width=___)\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(posterior)\n",
    "_ = ax.set(xlabel='sample', ylabel='mu');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ideally, our tracplot in each parameter would be **mixing well** (moving across parameter space without getting stuck), and carving out the same patch of parameter space on average. **The acceptance rate of new samples should be somewhere between $\\sim20-50\\%$ depending on the type of problem you're trying to solve**.\n",
    "\n",
    "Take a look at these plots below for an arbitrary problem. The 1st column is the trace, the 2nd is the histogram of the chain, and the 3rd column is the acceptance rate of newly proposed samples. \n",
    "\n",
    "- **In the top row, the proposal width was too small**. The chain hardly moves, giving it a high acceptance rate but poor global exploration of the posterior surface.\n",
    "- **In the middle row, the proposal width was too big**. The chain proposes new points that are far away from regions of high probability density, meaning that they're hardly ever accepted. The chain gets stuck.\n",
    "- **Only the bottom row shows reasonable sampling. This is the Goldilocks scenario.**\n",
    "\n",
    "![](./figures/fig_taylor_mcmc.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### More Complex Models & Sophisticated Sampling\n",
    "\n",
    "- The example you worked through above was overkill in that we were only trying to estmate $\\mu$. Note also that we can do this in fewer than 10 lines using the ${\\tt pymc3}$ module.\n",
    "\n",
    "\n",
    "- The process is essentially the same when you add more parameters. Check out this [animation of a 2-D process](http://twiecki.github.io/blog/2014/01/02/visualizing-mcmc/) by the same author whose example we just followed.\n",
    "\n",
    "\n",
    "- Also, check out [this](https://chi-feng.github.io/mcmc-demo/app.html?algorithm=AdaptiveMH&target=multimodal) amazing visualization tool, and have some fun exploring the many different types of sampling algorithms. Metropolis-Hastings is only the start!\n",
    "\n",
    "\n",
    "- In fact, the most basic type of MH algorithm I would ever personally use is **Adaptive Metropolis-Hastings**, where the proposal width updates itself in response to the currently-estimate covariance matrix of the MCMC chain. The chain learns how to optimize itself!! \n",
    "\n",
    "\n",
    "- We'll look at other algorithms and samplers later this week. There are plenty of tools that a professional astrophysicist can have in their toolbox to tackle research problems."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:astr8070] *",
   "language": "python",
   "name": "conda-env-astr8070-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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}