{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BAYa class Assignment 2019\n",
    "\n",
    "In this assignment, your task will be to implement and analyze inference in Bayesian Gaussian Mixture Model (GMM) as described in the slides on [Approximate Inference in Bayesian Models](http://www.fit.vutbr.cz/study/courses/BAYa/public/prednasky/4-Approximate%20inference.pdf). The preferred and easiest way of accomplishing this task is to complete this Jupyter Notebook, which already comes with training data definition and a code for Maximum Likelihood (ML) training of a contrastive GMM system. If you do not have any experience with Jupyter Notebook, the easiest way to start is to install Anaconda3, run Jupyter Notebook and open this notebook downloaded from [BAYa_Assignment2019.ipynb](http://www.fit.vutbr.cz/study/courses/BAYa/public/notebooks/BAYa_Assignment2019.ipynb). You can also benefit from reusing pieces of code from the Jupyter Notebooks provided for this class, namely: [vb_gmm_training.ipynb](http://www.fit.vutbr.cz/study/courses/BAYa/public/notebooks/vb_gmm_training.ipynb) and [gs_gmm_training.ipynb](http://www.fit.vutbr.cz/study/courses/BAYa/public/notebooks/gs_gmm_training.ipynb).\n",
    "\n",
    "\n",
    "The following cell contains a code with the definition of training data and the contrastive system training. You should not edit this part! The code does the following:\n",
    "\n",
    "1. We \"handcraft\" a GMM (i.e. define GMM parameters) with 4 Gaussian components, which represents the \"true distribution\" of training data.\n",
    "2. We have pre-generated training data from this GMM (see definition of variable x), so that everyone works with exactly the same data.\n",
    "3. GMM with C=6 components is trained on the training data using ML training (the standard EM algorithm). You will use this GMM as a contrastive model. You will compare it to your implementation of VB GMM.\n",
    "4. The following plots are made:\n",
    " * The true training data distribution (in grey).\n",
    " * The training observation generated from this true GMM (black + at the x axis).\n",
    " * The ML GMM estimate obtained using the EM algorithm (in black)\n",
    " $$\n",
    "\\DeclareMathOperator{\\aalpha}{\\boldsymbol{\\alpha}}\n",
    "\\DeclareMathOperator{\\bbeta}{\\boldsymbol{\\beta}}\n",
    "\\DeclareMathOperator{\\NN}{\\mathbf{N}}\n",
    "\\DeclareMathOperator{\\ppi}{\\boldsymbol{\\pi}}\n",
    "\\DeclareMathOperator{\\mmu}{\\boldsymbol{\\mu}}\n",
    "\\DeclareMathOperator{\\llambda}{\\boldsymbol{\\lambda}}\n",
    "\\DeclareMathOperator{\\diff}{\\mathop{}\\mathrm{d}}\n",
    "\\DeclareMathOperator{\\zz}{\\mathbf{z}}\n",
    "\\DeclareMathOperator{\\XX}{\\mathbf{X}}\n",
    "\\DeclareMathOperator{\\xx}{\\mathbf{x}}\n",
    "\\DeclareMathOperator{\\YY}{\\mathbf{Y}}\n",
    "\\DeclareMathOperator{\\NormalGamma}{\\mathcal{NG}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm4AAAHdCAYAAACkHvVVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3XlwXOd95vvvi8ZCLCQA7gI3LAT3neACkdopL5KsyKOrazuRx/bYsuMZO9eem1Sia4/tSTKJx1UapZJxTWLfuuN4q0jeZCmSJYtSSEqiKO4kwE0EBRAgQGKhCJAgtmbjvX80TwsAu4HTjT7oBc+nSkXy9DnnfSGJxMPfuxlrLSIiIiKS/DIS3QERERERcUfBTURERCRFKLiJiIiIpAgFNxEREZEUoeAmIiIikiIU3ERERERShIKbiIiISIpQcBMRERFJEQpuIiIiIikiM9Ed8MrMmTNtaWlporshIiIiMqZDhw51WGtnjXVf2ga30tJSDh48mOhuiIiIiIzJGHPezX0aKhURERFJEQpuIiIiIilCwU1EREQkRSi4iYiIiKQIBTcRERGRFKHgJiIiIpIiFNxEREREUoSCm4iIiEiKUHATERERSREKbiIiIiIpQsFNREREJEUouImIiIikCAU3ERERkRSh4CYiIiKSIhTcRERERFKEgpuIiIhIilBwExGZQNbaRHdBRFKYgpuIyAT5b//tv5GRkYHf7090V0QkRSm4iYhMkL/8y78E4OWXX05wT0QkVSm4iYhMkMrKSgCOHj2a4J6ISKpScBMRmQDWWpqamgBob29PcG9EJFUpuImITICrV69y9epVADo6OhLcGxFJVQpuIiIT4MqVK6GfK7iJSKwU3EREJkBXV1fo5wpuIhIrBTcRkQngBLd58+YpuIlIzBTcREQmgBPcKioqFNxEJGYKbiIiE8AJbvPnz6e3t5dAIJDgHolIKlJwExGZAE5wKykpAaCnpyeR3RGRFKXgJiIyAUYGt+vXryeyOyKSohTcREQmwNWrV8nOzmbGjBkAdHd3J7hHIpKKFNxERCbAtWvXmDp1KgUFBYAqbiISGwU3EZEJ0NPTQ15eHvn5+YAqbiISGwU3EZEJ4AQ3VdxEZDwU3EREJsDIipuCm4jEQsFNRGQCaKhUROJBwU1EZAJoqFRE4kHBTURkAvT29qriJiLjpuAmIjIBnIpbXl4eEAxyIiLRUnATEZkAPT095ObmkpmZic/no6+vL9FdEpEUpOAmIjIBnIobwJQpUxTcRCQmCm4iIhNgaHDLzc1VcBORmLgKbsaYDGPM140xp40xfcaYJmPMU8aYfBfPLjHG/KUxZp8xpt0Yc80Yc9QY841wzxtjvmOMsRH++dNYvkgRkUQKBAL09/cPq7hpjpuIxCLT5X1PA38C/AZ4Clh+89frjTE7rLWDozz7H4D/BDwP/AzwA/cAfw38n8aYrdbacH+CfR3oGHHtkMv+iogkDae6pqFSERmvMYObMWYl8FXg19baR4dcrwf+Hvgk8PNRXvFL4G+ttV1Drv2jMeYs8A3g88D/DPPcc9bahjG/AhGRJNfT0wMEh0hBwU1EYudmqPRTgAH+bsT1HwI9wOOjPWytPTgitDmeufnjqkjPGmOmGWPcVgVFRJJSf38/ADk5OYDmuIlI7NwEt03AILB/6EVrbR9w9ObnsZh/88fWCJ8fB7qAPmPMXmPMR2NsR0QkoUYGN1XcRCRWboJbCdBhre0P81kzMNMYkx1No8YYH/At4Aa3DrN2Aj8gODz7B8CTwCLgRWPMZ8d47xeNMQeNMQfb29uj6ZKIiGfCBTctThCRWLgZhswDwoU2gL4h9wxE0e7fAVuB/8dae2boB9bakUOyGGP+P6AWeNoY80trbdizYqy1PyAY+qiqqrJR9EdExDPhgpv+cikisXBTcesBciJ8NmXIPa4YY/4K+ArwA2vt37p5xlp7GfhHoAi43W1bIiLJQEOlIhIvboJbC8Hh0HDhbR7BYVRX1TZjzHeAbwL/G/hjt528qeHmjzOjfE5EJKG0OEFE4sVNcDtw877NQy8aY6YA64CDbhoyxnwb+DbwY+AL1tpohzIrb/4YaTGDiEhSUsVNROLFTXB7BrDA10Zcf4Lg3LafOReMMRXGmGUjX2CM+RbwHeAnwOcibdhrjMk0xhSGub4A+DJwGdjros8iIklDixNEJF7GXJxgra0xxnwf+Iox5tfAS3xwcsJuhq8KfY3gClDjXDDG/CfgvwKNwE7gD40xQx6h1Vr76s2fFwD1xpjngFPAFWAp8IWbn30qwikLIiJJSxU3EYkXt5vbfo3gHLMvAg8SPIrqH4BvjXHcFXywz9tC4J/DfL4bcIJbL/ArYAvwCMGw1kEw8H3PWrs/zPMiIkltZHDLyclhYGAAay0j/iIrIjIqV8HNWhsgeEbpU2PcVxrm2meBz7psp59gdU1EJG2MDG7Z2dlYawkEAmRm6nAYEXHPzRw3EREZh3DBbeh1ERG3FNxERDwWKbgNDESzb7mIiIKbiIjnws1xAwU3EYmegpuIiMdUcROReFFwExHxWH9/Pz6fD5/PByi4iUjsFNxERDzW398fqraBgpuIxE7BTUTEYwpuIhIvCm4iIh6LFNy0HYiIREvBTUTEY6q4iUi8KLiJiHhsZHDTdiAiEisFNxERj6niJiLxouAmIuIxBTcRiRcFNxERjym4iUi8KLiJiHhMwU1E4kXBTUTEY9oORETiRcFNRMRjqriJSLwouImIeEzBTUTiRcFNRMRj2sdNROJFwU1ExGOquIlIvCi4iYh4bGRwy8rKAhTcRCR6Cm4iIh4bGdwyMjLIzMxUcBORqCm4iYh4bGRwg+BwqbYDEZFoKbiJiHgoEAgQCATCBjdV3EQkWgpuIiIecqpqCm4iEg8KbiIiHooU3HJychTcRCRqCm4iIh5SxU1E4knBTUTEQ044c/Zucyi4iUgsFNxERDzk9/sBBTcRiQ8FNxERDznhzNl016HtQEQkFgpuIiIeUsVNROJJwU1ExEOjVdwU3EQkWgpuIiIeilRx03YgIhILBTcREQ+p4iYi8aTgJiLiIW0HIiLxpOAmIuIhLU4QkXhScBMR8ZC2AxGReFJwExHxkCpuIhJPCm4iIh7S4gQRiScFNxERD2k7EBGJJwU3EREPqeImIvGk4CYi4qFI24FkZWXh9/ux1iaiWyKSohTcREQ85AyVhqu4Ady4cWPC+yQiqUvBTUTEQ6NV3IZ+LiLihoKbiIiHIi1OcIKb83kqe/nll9mxYwcXLlxIdFdE0p6Cm4iIh5yKms/nG3Y9nYLbZz7zGV577TVee+21RHdFJO0puImIeMjv95OdnY0xZth1pwKXDsHNCadnz55NcE9E0p+Cm4iIhwYGBm5ZmADpM8etu7ubzs5OAN59990E90Yk/Sm4iYh4yKm4jZQuQ6X19fVhfy4i3lBwExHx0FgVt1QPbhcvXgSgtLSUy5cvJ7g3IulPwU1ExEMDAwNhK27OtVQfKr1y5QoAlZWVCm4iE0DBTUTEQ+k+VPr+++8DUFFRwdWrV1P+6xFJdgpuIiIeSvehUqfitnjx4mG/FhFvKLiJiHhoMlTccnNzKSkpAdBwqYjHFNxERDwUqeKWTnPciouLmT59OvDB0KmIeEPBTUTEQ5Oh4jZ9+nQFN5EJouAmIuKhdJ/j1tnZSVFREVOnTgXg2rVrCe6RSHpTcBMR8VCk7UDSJbhdu3aNadOmKbiJTBAFNxERD/n9/rSe49bd3U1BQQEFBQWhX4uIdxTcREQ8lO4Vt5HBTRU3EW+5Dm7GmAxjzNeNMaeNMX3GmCZjzFPGmHwXzy4xxvylMWafMabdGHPNGHPUGPONSM8bY5YaY54zxlwxxlw3xrxhjLk3mi9ORCTR0n1xwrVr1ygoKMDn85GXl6fgJuKxaCpuTwP/AzgJfBX4BfAnwAvGmLHe8x+ArwPngL8E/gw4A/w1sNcYkzv0ZmNMBbAXqAa+d/P+AuAVY8yOKPosIpJQY20HksrBzVpLd3d3aH7b1KlTNVQq4rFMNzcZY1YSDGu/ttY+OuR6PfD3wCeBn4/yil8Cf2ut7Rpy7R+NMWeBbwCfB/7nkM/+FigCNlprj95s68fACeD7xphl1lrrpu8iIok0VsUtlee49fX1MTg4GBomLSgoUMVNxGNuK26fAgzwdyOu/xDoAR4f7WFr7cERoc3xzM0fVzkXbg6dPgzsckLbzXd0A/8vsATY5LLfIiIJlc7bgTjVNSe4TZ06VcFNxGNug9smYBDYP/SitbYPOErsQWr+zR9bh1xbA+QAb4e5f9+Q/oiIJL10nuMWLrhpqFTEW26DWwnQYa3tD/NZMzDTGHPrn0yjMMb4gG8BNxg+zFoy5L3h2gKYF+GdXzTGHDTGHGxvb4+mOyIinhir4pbKQ6WquIlMPLfBLQ8IF9oA+obcE42/A7YC37LWnhnRFhHaG7Uta+0PrLVV1tqqWbNmRdkdEZH4i7QdiM/nIyMjIy0qbs7ihNzcXHp6ehLZJZG05za49RAcvgxnypB7XDHG/BXwFeAH1tq/DdMWEdqLui0RkUSKNFQKwapbOgS3vLzg36Vzc3Pp7e1NZJdE0p7b4NZCcDg0XJiaR3AY1VW93xjzHeCbwP8G/jhCW857w7UF4YdRRUSSyuDgIIFAIOxQKaR+cHNCWm5ucEenvLw8BTcRj7kNbgdu3rt56EVjzBRgHXDQzUuMMd8Gvg38GPhChC09aggOk1aH+WzrzR9dtScikkhOKItUccvOzk7pOW4jg5sqbiLecxvcngEs8LUR158gON/sZ84FY0yFMWbZyBcYY74FfAf4CfA5a+1guIZubvvxAnC3MWbtkOcLgC8AZxmxulVEJBk5oWyyVNwU3ES852oDXmttjTHm+8BXjDG/Bl4ClhM8OWE3w1eFvgYsIrjvGwDGmP8E/FegEdgJ/KExZsgjtFprXx3y6yeB+4DfG2OeBq4SDInzgAe1+a6IpIKxKm6pHtz6+oLrxYYGN7/fz40bN8jMdPXtRUSiFM3vrK8BDcAXgQeBDuAfCK4KDVs9G8LZd20h8M9hPt8NhIKbtbbOGLMN+C7wF0A2cBj4iLV2ZxR9FhFJmMlYcXOuOytNRSS+XAc3a20AeOrmP6PdVxrm2meBz0bTMWvtKeAPonlGRCSZOMFtssxxc1aXKriJeCeaQ+ZFRCQKTjUtnStuxphQMB1acRMRbyi4iYh4ZKyKWzoEt9zcXJw5ywpuIt5TcBMR8Yib7UBSPbhNmTIl9GsFNxHvKbiJiHjEzeKEVJ/j5oQ1+CC46dgrEe8ouImIeCTdtwOJFNxUcRPxjoKbiIhH0n07kL6+vmHBbeiqUhHxhoKbiIhHJsMcN1XcRCaWgpuIiEcm6xw3BTcR7yi4iYh4ZLJsB+LQ4gQR7ym4iYh4RIsTRCTeFNxERDwy1lBpOhx5pcUJIhNLwU1ExCOToeI2dAPerKwsMjIyFNxEPKTgJiLikXTfDmRkxc0YQ25uroKbiIcU3EREPDIZKm5Dgxug4CbiMQU3ERGPpPMcN2vtLRvwQjC4aVWpiHcU3EREPJLO24H09/cD3BLc8vLyVHET8ZCCm4iIR5xQNtoct8HBQQYHByeyW3HhhDMNlYpMLAU3ERGPOBW3zMzMsJ87lbhUrLopuIkkhoKbiIhH/H4/2dnZGGPCfu5U4lJxnpuCm0hiKLiJiHhkYGAg4jApfBDc0q3ipsUJIt5RcBMR8YhTcYskHYLb0A14nV87CxdEJP4U3EREPDJWxS0d57hNmTKFvr6+RHRJZFJQcBMR8Yjf73c1VJqKc9yccDYyuOXk5KjiJuIhBTcREY8MDAyQk5MT8fNUHip1wtnIr0/BTcRbCm4iIh5J58UJTsVtZHDTUKmItxTcREQ8MjAwMOriBOezVBwqVcVNJDEU3EREPJLOq0pHC25+vz8lT4MQSQUKbiIiHhmr4pYOwW3kdiBOkEvFKqJIKlBwExHxSDrPcYtUcXOCnOa5iXhDwU1ExCOTdY7b0M9FJL4U3EREPDJZ57gN/VxE4kvBTUTEI+k+xy0zM5OMjOHfRpzgpqFSEW8ouImIeCTdg1u4zYWdOW6quIl4IzPRHRARSVduzypN1Tlu4YJbKg6VBgIBDh48yOHDh7l8+TIFBQWsWrWKO+64Y9STL0QSQcFNRMQj6T7HLR2C2/Xr13nmmWdoampiwYIFbN68mffff5+9e/dSW1vL448/zsyZMxPdTZEQBTcREY+k81BpX1/fqEOlqTDHrb+/n5/+9Kd0dHTw6KOPsnLlSowxADQ1NfHMM8/wox/9iCeeeILCwsIE91YkSHPcREQ84nY7kFQMbqlecbPW8txzz9HW1sYnPvEJVq1aFQptAwMDnDhxgoKCAi5dusSzzz5LIBBIcI9FglRxExHxiNsNeDXHbeIdOXKE06dPc//997N48eLQ9d27d/P4449z4cKF0LV169ZRVlbGjh07EtFVkWFUcRMR8Ug6D5X29/ffctwVpEZw6+3tZefOnSxatIjq6urQ9Z07d3L//feTl5fHCy+8QG1tLX/2Z3/GsWPH+NKXvkRbW1sCey0SpIqbiIgHrLWTcnFCKsxx27VrF319fXz0ox8NDY82Njby6KOPsmzZMvbs2UNRUREA3/ve9ygvL+fLX/4yn/vc53jxxRcT2XURVdxERLxw48YNAFfBTUOlE+f999/nwIEDbNy4kTlz5gDBkP35z3+ewcFBfvvb34ZCm+OP//iPeeSRR3jppZd4/vnnE9FtkRAFNxERDzhhbLQ5bsYYMjMz06riluzB7a233iIjI4M777wzdO35559n586dfPe736WsrCzsc//0T/9EcXExX/3qV7VQQRJKwU1ExANOcBut4gbBYJdOwS2Zh0qvXr3KsWPHWL9+PVOnTgWCm+8++eSTLFu2jC996UsRn509ezZPPPEEjY2N/OhHP5qgHovcSsFNRMQDThibbMEtmStu+/btY3BwkNtvvz107cUXX+TUqVN85zvfITNz9GnfX//615kzZw5//dd/jbXW6+6KhKXgJiLiAbcVt+zs7LSa4+YcPJ9swc3v93PkyBGWL19OcXFx6PrTTz/NggULePTRR8d8x9y5c3nooYdoaGhg586dXnZXJCIFNxERD7iZ4+Z8nk4VNwhW3ZItuJ04cYK+vj6qqqpC144ePcquXbv46le/Oma1zfEf/+N/JC8vj//+3/+7V10VGZWCm4iIB9J9jlukI68gOM8t2ea4HTx4kJkzZ1JaWhq69o//+I/k5eXxhS98wfV71qxZw5YtW3j99ddpaGiIf0dFxqDgJiLigck6xw2Sr+LW1tZGc3MzGzduDO3b1t/fz7PPPsvHP/7xYUOnY8nMzOSTn/wk1lp+8pOfeNVlkYgU3EREPJDOc9ystSkV3I4fP44xhtWrV4eu/e53v+PKlSs8/vjjUb9vx44dLFiwgH/+53/WIgWZcApuIiIeSOc5bk5/wx15BcHglixDpdZaamtrqaioID8/P3T9pz/9KbNnz47p/NHS0lI2b97MuXPnqKmpiWd3Rcak4CYi4oF0nuPmVNNGm+OWLBW3xsZGurq6hlXbrl69yr/+67/yiU98wvWihKEyMjJ49NFHycjI4Mc//nE8uysyJgU3EREPuJ3jlp2dnXbBLZmGSmtqasjKymLZsmWha6+88gr9/f089thjMb93y5YtlJeX84tf/ELDpTKhFNxERDwQTcUt1ea4pUpwCwQCnDx5kqVLlw777/D8888zY8YMqqurY373okWLWLlyJY2NjZw+fToe3RVxRcFNRMQDk32oNBnmuNXX19Pb28uqVatC1/x+Py+++CIPPfRQTMOkDp/PxwMPPACgg+dlQim4iYh4IJ0XJ6RKxe306dNkZWVRUVERuvbWW29x5coVHn744XG/f9u2bcydO5df/epX436XiFsKbiIiHohmjpuGSuPPWsuZM2eorKwcVll74YUXyMnJ4UMf+tC42ygvL2fp0qUcOnSIjo6Ocb9PxA0FNxERD0zmodJkCG7Nzc10d3ezdOnSYddfffVVtm/fTkFBwbjbyM3NZdu2bQwODvLqq6+O+30ibii4iYh4IJ2DmzN/LZnnuJ0+fZqMjAwqKytD11pbW6mpqYlp77ZI7rnnHqZMmcIrr7wSt3eKjEbBTUTEA5rjltiK2+nTpyktLSU3Nzd07fXXXweIa3BbsmQJZWVlvPLKK9oWRCaEq+BmjMkwxnzdGHPaGNNnjGkyxjxljMkf+2kwxjxpjPmFMeY9Y4w1xjSMcu+Pbt4T7p//w+XXJSKSUJrjlrjg1tHRweXLl28ZJt25cydFRUWsX78+bm3NmzePyspKLl26RF1dXdzeKxKJ27XQTwN/AvwGeApYfvPX640xO6y1g2M8/zfA+8BhoMhlm58Oc22/y2dFRBIqnYdKnVAW6cirRA+VOgFq6DCptZadO3dy77334vP54taWz+fjnnvu4fnnn2fnzp3D2hTxwpjBzRizEvgq8Gtr7aNDrtcDfw98Evj5GK+psNa+d/O5WmDMWaHW2p+OdY+ISLKa7EOlAwMDWGsxxkxk14BgcJsxYwbFxcWha+fOnaOxsZE///M/j3t71dXVFBYW8uKLL/LlL3857u8XGcrNUOmnAAP83YjrPwR6gMfHeoET2qJhgqYZYzQPT0RSzsDAABkZGWNWd9I1uAEJGQL2+/2cP3+exYsXD7u+c+dOIL7z2xwVFRWUlZXx1ltvaZ6beM5NKNoEDDJimNJa2wccvfm5F7pu/tNrjHnVGLPFo3ZEROJuYGBgzGFSSN85bkPvm0jnz5/nxo0btwS3f/u3fwvNR4u3GTNmUFlZSWdnJ2fOnIn7+0WGchPcSoAOa22434HNwExjzNh/Orl3ieCcui8DHyc4P64KeMMYE/+/KomIeMDv97sKbk7FLZUqNW6OvAISMs+trq6OzMxMFi1aFLpmreXNN9/kjjvu8GTo1hjDtm3bAHjjjTfi/n6RodwEtzwg0l+b+obcExfW2r+w1v5na+3PrLXPWWv/K7AZ8AP/a7RnjTFfNMYcNMYcbG9vj1eXRESi5rbi5syBCwQCXncpbpK54lZXV0dpaemwuYXnz5+npaWF7du3e9buli1byM/PD205IuIVN8GtBwj/uxOmDLnHM9bas8CzwGJjzJJR7vuBtbbKWls1a9YsL7skIjKqgYGBMRcmwAerTlNpnluyBrcrV65w+fLlYWeTQvB8UiBUFfNCaWkpCxcuVMVNPOcmuLUQHA4N9zt0HsFh1ImYoNFw88eZE9CWiMi4RFtxS6V5bv39/fh8vogLL5zgNtFDpefOnQO4ZX7bW2+9xdSpU1m9erVnbc+ePZuKigqam5tpbm72rB0RN8HtwM37Ng+9aIyZAqwDDnrQr3CcGaWtE9SeiEjMopnj5tyfKvr6+iJW2+CDOW4TXXGrq6ujsLCQGTNmDLv+5ptvUl1dHdf920bKyMhg69atofZEvOImuD0DWOBrI64/QXBu28+cC8aYCmPMslg7Y4zJvxkIR15fDzwGnLLWnov1/SIiEyXailsqBbf+/v5Rg1sihkoHBwdpaGigvLx82AKEzs5OamtrPR0mdWzfvp3s7GzNcxNPjbkBr7W2xhjzfeArxphfAy/xwckJuxm++e5rwCKC+76FGGM+ffM6wCwg2xjzzZu/Pm+t/cnNn1cCvzPGPAecBa4Da4H/AASAL0b9FYqIJEC6z3FLtuB26dIl+vv7KSsrG3Z93759WGs9XZjgqKioYP78+ezZs8fztmTycnvk1dcIzjH7IvAg0AH8A/AtF8ddAXweuGvEtb+6+eNuwAlul4CdwD3AHwG5wEWCVb+/tdaedtlfEZGESvc5bpGOu4LEbAdSX18PcEtwe/PNN/H5fGzZ4v1WoHPnzmXBggW88cYb9PT0kJcXtw0XREJcBTdrbYDgGaVPjXFfaYTrd7ts5xLhzygVEUkp6TzHLRkrbg0NDcyaNYuCguEnKu7bt481a9aQn5/veR8yMzNZs2YNu3fv5vDhwxNS5ZPJR8dJiYh4QHPcJi64BQIBzp8/T2lp6bDrg4ODHDx4kM2bN4d/0AO33347AG+//faEtSmTi4KbiIgHop3jlmpDpcm0qrS5uRm/33/LMGldXR1dXV1s2uTVyYy3WrVqFUVFRdrPTTyj4CYi4gFV3CZujpszv23oMVcABw4cAJjQ4DZ//nzmzZvHwYMTtVOWTDYKbiIiHtAct4mruDU0NDB37txbFgMcOHCA3NxcVqxYMSH9ACgsLKSsrIyLFy/S2qptRyX+FNxERDwwmStuE7mq1O/309TUdMswKQSD24YNG8jMdLuBwvgZY9i4cSMA+/fvn7B2ZfJQcBMR8YDb4JaOc9wmsuJ24cIFAoHALcHN7/dz5MiRCR0mdWzfvh1jjE5QEE8ouImIeMDt4oRUrLj19fWNuo9bZmYmxpgJCW719fUYY1i4cOGw6ydOnKC3t3dCV5Q6Fi9ezJw5cxTcxBMKbiIiHpjMc9yMMUyZMmVChkobGhooKSm5pT+JWJjgKCkpYd68edTU1GCtnfD2Jb0puImIeGAyz3GD4HCp1xU3v99Pc3PzLatJIRjciouLqaio8LQP4WRnZ7N06VKuXbvG+fPnJ7x9SW8KbiIiHpjMc9xgYoJbc3Mzg4ODEYNbVVXVsAPnJ9KGDRsAOHToUELal/Sl4CYi4oF0nuPmJrhNxFCpU81asGDBsOu9vb3U1NQkZJjUsWXLFowxOkFB4k7BTUQkzgKBAIODg2k5VGqtHXNxAkxMxa2xsZE5c+aQm5s77PqxY8cIBAJUVVV52v5oysvLmTVrlrYEkbhTcBMRiTMnhEUzVJoqwc3pZ6KHSgcHB2lqarplNSnA0aNHgQ+GKxNh9uzZlJSUaIGCxJ2Cm4hInDnz1aKpuKXKHDcnjCV6qPTixYv4/f6Iwa2oqCjsZxMlMzOTpUuX0tnZycWLFxPWD0k/Cm4iInHmhLB0nOPmNrh5XXFrbGwEbj2fFILBbd26dQlbmOBYv349AIcPH05oPyS9KLiJiMRZLBU3BbfoNDY2UlxczNSpU4ddDwQpB5OMAAAgAElEQVQCHD9+PBSaEmnr1q0A7N27N8E9kXSi4CYiEmfRzHHz+XwYY1JmqNQZ/hxrcYKXQ6XWWhobG8NW29599116e3tZt26dJ21Ho7KykhkzZvDOO+8kuiuSRhTcRETizKk0uRkqNcaQlZWlilsUOjo66OnpGXVhQjIEN2eBQm1tbaK7ImlEwU1EJM6cwDJWVcqh4BYdZ35bpOCWnZ3N8uXLPWk7GpmZmVRWVtLW1kZHR0eiuyNpQsFNRCTO3IYbRzoGNy+HShsbG8nPz2f69Om3fHbkyBFWrVrlqto5EZy5djpBQeJFwU1EJM6iDW7Z2dkpM8ctGSpu58+fZ9GiRbesGrXWhlaUJovq6moA9u3bl+CeSLpQcBMRibN0rrg5VbREBbeuri66urrCDpNevHiR9vb2pFhR6li+fDnTpk3jwIEDie6KpAkFNxGROHMbbhypFNzczt/zaqjUOZ803IrSI0eOAMmxMMExe/Zs5syZw8mTJxPdFUkTCm4iInGmxQnBzwcGBuJ+3FNjYyM5OTnMnj37ls+cFaVr1qyJa5vjkZmZSVlZGY2NjZ6f3SqTg4KbiEicaY7bB5/H++tqbGxkwYIFZGTc+u3r6NGjLF68mGnTpsW1zfFatWoVgUCA06dPJ7orkgYU3ERE4iyd57hFs6oUiOtwaU9PD+3t7RHPIE22hQmOjRs3Amiem8SFgpuISJylc3CLZnECENfhwdHOJ7169Sp1dXVJGdyqqqrw+XxaWSpxoeAmIhJn0S5OyM7OTpng5nb+nlfBzefzUVJScstnx48fB5JrYYJj/vz5zJ49O9RHkfFQcBMRibNYFiek2xw3L4ZKz58/z/z588nMzLzls5qaGiC5FiY48vLyWLBgAWfPnk10VyQNKLiJiMRZOg+VOl9bdnb2qPfFu+I2MDDAxYsXI85vq62tpbCwkPnz58elvXhbunQpnZ2dtLa2JrorkuIU3ERE4qy/vx9jTNjKUDipFtyysrLCruocKt7B7cKFC1hrIwa3mpoaVq1adctpCsli7dq1ABw+fDjBPZFUp+AmIhJn/f395OTkuA4RqbQdSF9fn6tKYryHShsbGzHGsGDBgls+s9ZSU1PD6tWr49KWF7Zu3QrA22+/neCeSKpTcBMRibO+vj7X89sg9Spubr62eFfcGhsbmTt3btjQ2NzcTGdnZ1IHt+XLlzN16lRV3GTcFNxEROLMqbi5lWrBzc3XFs/gFggEuHDhQthqG3ywMCGZg1txcTG33Xabjr6ScVNwExGJMwW3+A6VXrp0Cb/fH3b/NvgguK1atWrcbXnFGENFRQWNjY0pMywuyUnBTUQkzqINbqk0xy0RFTfnYPnRFibMmzeP4uLicbflJR19JfGg4CYiEmfpXHFzuzghnsGtqamJ6dOnU1BQEPbz2trapB4mdWzatAmAvXv3JrgnksoU3ERE4kyLE+I3VGqtpbGxMWK17caNG5w6dSolgtvmzZvx+XwcPHgw0V2RFKbgJiISZ+lccZvoodLLly/T09MTMbidPXuW/v7+lAhuJSUlzJw5k9ra2kR3RVKYgpuISJxpjlv8gpub+W2Q3AsTHFlZWSxYsIBz584luiuSwhTcRETiLJaKWyAQwFrrYa/iY6JXlTY1NZGfn8/06dPDfl5TU4PP52P58uXjameiLF26lI6ODq5evZrorkiKUnATEYmzWIIbkBLDpW6/tszMTDIyMuJScVu4cGHEUyhqamqorKyMak5hIjlDutqIV2Kl4CYiEmfRLk5wDmxPheAWzdeWk5MzruB29epVOjs7Iw6TQuqsKHVs3rwZ0NFXEjt3JyCLiExi/f391NfX09bWRn9/P9nZ2cydO5dFixaFDTGxVtwGBgbIz8+PW7+9EM3XNmXKlHENlTY2NgKR57ddv36d9957j8985jMxtzHRNm7cSFZWFkePHk10VyRFKbiJiERw9epVdu3aRU1NDTdu3ADA5/MRCASA4HDg6tWrueuuuygsLAw9p6HSoPFW3BobG8nKymLu3LlhPz9x4gTW2pRYmOCYNm0at912G6dOnUp0VyRFKbiJiIxgreXAgQO8+uqrWGtZu3Ytq1evpqSkhOzsbPx+Py0tLRw/fpzjx49TU1PDPffcQ3V1NcYYBbeb4hHcFixYQEZG+Fk9qXBGaThlZWWa4yYxU3ATERnixo0bPPfcc5w4cYLKykoeeOABioqKht2TlZXFokWLWLRoEXfeeSe/+93vePXVV2lqauKRRx6JaTsQSL/gNp6h0r6+PlpbW7n77rsj3lNTU0NeXh7l5eUxtZEoy5YtY/fu3Vy6dCliNVEkEi1OEJFJw1rLq6++yp//+Z/zgx/8gJ6enmGfDwwM8POf/5wTJ05w33338alPfeqW0DZSYWEhn/jEJ/jwhz/MmTNn+MlPfhLTyQlO+8lscHAQv98/IYsTmpqagMjz2yAY3FauXBmxIpesNmzYAMBbb72V4J5IKkqt/9tFRGLU39/PY489xoc+9CGeeuopvvSlL7F161YaGhoACAQCPPvsszQ0NPDII4+wffv2iFtQjGSMYevWrTz22GNcuHABay2Zme4HNFJlqNQJYRMxVNrY2EhGRgbz5s2LeE+qrSh1VFdXA3DgwIEE90RSkYKbiEwKTzzxBL/61a/4m7/5G65fv86LL75IU1MTH/vYx+ju7uaFF17g3LlzfOxjH2Pt2rUxtbF8+XIefvhhAM6cOcPg4KCr59I1uI1nqLSxsZHbbrstNIw8UltbG21tbSm1MMGxfPly8vLyOH78eKK7IilIwU1E0t4vfvELfvKTn/Dtb3+bJ598kpycHB544AGeeeYZTpw4wec//3mOHTvGXXfdxfr168fV1qJFiwC4cuUKO3fudPWME06Sfah0oipuN27coLm5ecxhUki9hQkQXI08f/58zp49m+iuSApScBORtNbX18ef/umfsmHDBr75zW8O++xDH/oQn/70p/nlL39Jbm4ud911V1zaA6ioqODtt9/m9OnTYz6TrhW3WINbS0sLgUAgbYMbwOLFi2lqakqJY84kuSi4iUha+6d/+icaGxv53ve+d8u8s4GBAZYuXUp2djb79+93PadtNE5Q2bBhA7fddhu//e1v6ezsHPWZVAluTij1eqjU2Xh3wYIFEe+pqalh1qxZzJkzJ+r3J4NVq1bR39/vKtiLDKXgJiJpKxAI8PTTT3PHHXdw33333fL566+/jt/v5zOf+QzPPfcc58+fH3ebTnDLy8vjscceY3BwkN/+9rejVlZSJbg5X1s0q0pjDW4zZ84c9RSJmpqalJzf5ti4cSMAe/fuTXBPJNUouIlI2nrhhRc4f/48X/va12757MKFC7zzzjtUVVXxjW98A2MM3//+98fdZm9vLxAMbsXFxXzoQx+ioaGBQ4cORXwmXee45ebmRh3cBgcHaWpqGrXaNjg4yIkTJ1J2mBRg27ZtAKP+fyESjoKbiKStH/7whyxYsCC00tMRCAR4/vnnmTZtGjt27GDBggU89NBD/PjHPw4dbRUrZ2+43NxcIDhkWlZWxquvvkpXV1fYZ1Kt4hZNcHOCrFutra309fVRWloa8Z76+np6enpSOrjNmzePwsJCTp48meiuSIpRcBORtNTe3s4rr7zCH/3RH90yt+3gwYO0t7fz0Y9+NBRCPvvZz9La2sorr7wyrnad4JaXlwcE93j72Mc+hrU24rvTNbhNmTIl6uDmDFePFtychQmpPFQKwa/x3Llzie6GpBgFNxFJS7/85S8JBAL84R/+4bDrvb297Nq1i7KyMpYuXRq6/sADDzBz5kx+9rOfjavdkcENoLi4mO3bt3Pq1Cnee++9W55JleAW7eIEZ6g0mpWTDQ0NFBcXM23atIj31NbWArBy5UrX701GS5cu5dKlS+M6z1UmHwU3EUlLP//5z1m1atUtw2m7du2iv7+fD3/4w8NWkWZlZfHwww/z0ksvjWuu2dA5bkPdfvvtFBcX8/LLLxMIBIZ9ls5z3IY+NxZrLefPnx+12gbB4FZWVsbUqVNdvTdZrVmzhhs3bugEBYmKgpuIpJ2mpibefPNNPvWpTw27/v7773Pw4EHWr18fdhuJRx55hK6uLnbt2hVz2+EqbhDcdPXDH/4w7e3tHDx4cNhnXlXc6urq+O53v8vhw4fj8j6n4uYEsrE497kdLnXmtzmbGEdSW1ub8sOkAFu2bAFg3759Ce6JpBLXwc0Yk2GM+box5rQxps8Y02SMecoYE3m99vDnnzTG/MIY854xxhpjGsa4f4sxZqcx5pox5qox5mVjzDq3/RWRyevFF18E4OMf//iw63v27CEjI4O777477HM7duwgLy+P5557Lua2Ry5OGGrJkiWUlZWxZ8+eYVUoL4JbU1MTW7Zs4cknn+T222+/JSzGwglgXgU359zY0SpuAwMDnDlzJi2CW3V1NcYYjh07luiuSAqJpuL2NPA/gJPAV4FfAH8CvGCMcfOevwHuBc4BV0a70RizFdgNlAHfAr4NVAJvGGNSdxmRiEyIF198kbKyMpYtWxa6dvnyZY4fP05VVVXEIbbc3Fw+8pGP8Nvf/tb1OaMjRaq4QXChwn333UdPTw9vv/126LozVBrP4PZf/st/obe3lz179jB9+vSwW6JEy+vgdv78eYqLiyksLIx4z5kzZ7hx40ZaBLf8/HxmzZrFmTNnEt0VSSGugpsxZiXBsPZra+2/s9b+0Fr7n4H/DNwDfNLFayqstTOstfcDLWPc+/fAAHCntfZpa+3TwJ2ABZ5y02cRmZx6e3t57bXXePDBB4fNYdu9ezeZmZmh/bMiefjhh2lpaYm5CuIEt0ib1M6bN4/ly5fz9ttvc/36deCDilu85rh1dXXxzDPP8NnPfpY77riDv/iLv+Ctt94a955hTgBzuwGvc5+b4ObMb3MzTAqpv6LUUV5eHqo0irjhtuL2KcAAfzfi+g+BHuDxsV5grb11KVUYxpjFwCbgF9ba5iHPNxOs8u0wxsx12W8RmWR27dpFb28vDz74YOhae3s7NTU1bNq0iYKCglGfv//++wFcHxA/Um9vL7m5uWRkRP7j9d5778Xv9/PGG28A8R8q/Zd/+Rf6+vr43Oc+B8C///f/nqysLP7lX/5lXO/1co5bW1sbvb29rhYmZGZmDqumprIVK1bQ0dFBR0dHorsiKcJtcNsEDAL7h1601vYBR29+Hi/Ou94O89k+ggFyYxzbE5E08tJLL5GXlzdsHtsbb7xBdnb2mNU2gJKSElasWBFzcOvp6Rkz2MycOZN169Zx8OBBrl69Gvfg9stf/pJly5ZRVVUFQFFREXfddVdo7l+sent78fl8of6Oxfn34Ob0BKfq5KbitmTJktDwcqrbsGED1lreeuutRHdFUoTb4FYCdFhrw63pbgZmGmPi9buoZMh7w7UFMC/cg8aYLxpjDhpjDra3t8epOyKSSl577TXuvPPO0DBdZ2cntbW1bNy4Mey8s3B27NjBG2+8EdM5mz09Pa7aufPOO0PfsDMyMvD5fHEZKu3v7+fNN9+8ZbuTBx98MOI+cm451US3oqm4NTQ0UFRURFFR0aj3pcuKUsftt98OEJfFIzI5uA1ueUCkjXj6htwTD857wrU3alvW2h9Ya6ustVWzZs2KU3dEJFW0trZy6tQp7rnnntC1ffv2YYwJbb3gxo4dO+jt7R22gMAtt8GtqKiINWvWcPjwYbq7u8nKyopLxW3fvn309fVx7733DrvuDB2//PLLMb+7t7fX9fw2cB/c3O7f1t3dzXvvvZdWwW3VqlX4fL7Q3D2RsbgNbj1ApB0Xpwy5Jx6c94RrL95tiUgacfZfc4ZJe3t7OXz4MKtXrx51peJId911Fz6fL6bh0t7eXteVve3btxMIBNi7dy/Z2dlxqbi9/vrrZGRkcOeddw67vnjxYubOnRtTGHX09fV5UnFz5reNNUzqnOuZymeUjpSVlcW8efM4e/ZsorsiKcJtcGshOBwaLkzNIziMGq8tv50Vp+GGQ51r4YZRRWSS27VrF1OnTmXDhg0AHDhwAL/fT3V1dVTvmTZtGlVVVezZsyfqPriZ4+aYMWMGq1at4uDBg+Tk5MQ0NDvSG2+8wbp1624ZcjTGsHXr1nFt9urVUKmb/dsg/VaUOiorK7lw4QI3btxIdFckBbgNbgdu3rt56EVjzBRgHRDPwXnn7I9wf9JuJbglyPjWtItIWtq1axd33HEHmZmZ3Lhxg/3791NZWRn2lISxbNu2jQMHDkR9jqTboVLHHXfcgd/vxxgz7jMrrbUcPnyYTZvCrxfbunUrdXV1Ma9gjDa4ud0OpL6+3vX8ttzcXMrKylz3IRWsXr2arq4uHTgvrrgNbs8QDEwjd3B8guB8s9CpzMaYCmNMzOu0rbV1BIPgY8YYZ6ECN3/+GPC6tfZSrO8XkfR08eJFTp8+HRomra2t5fr161FX2xzbtm2jv78/6uOiog1us2bNYsWKFQQCgdC+brF677336OrqYuPG8Avvt27dCsD+/fvDfj4WL+a4DQ4O0tDQQHl5+Zjvq62tZcWKFfh8Ptd9SAVO0B7PMLZMHq6Cm7W2Bvg+8O+MMb82xnzBGPMUwZMUdgM/H3L7a8Cpke8wxnzaGPNNY8w3gVlAofNrY8ynR9z+fxGc4/aGMeZrxpivAW/c7O//HeXXKCKTwO7duwG45557sNayf/9+Zs2aNebwWyTO1iHRbtMQzRy3oW35fD5aW1ujem4kZ4NdZ6h4pI0bN2KMiflQ81jnuI02BNzc3Ex/f7+r4FZTU5N2w6RA6C8X490gWSaHaI68+hrwp8BKgiHuk8A/AA9Za92cDfN54K9u/jMbKBry688PvdFauxe4G2gA/vrmPXUET1LQoW4icotdu3Yxbdo01q1bR3NzMxcvXmTTpk3DtsSIxpw5c6ioqIg6uEUzx81RUlJCfn4+7e3tBAKBqJ4d6tChQ2RlZUUMNwUFBVRUVFBTUxPT+6MdKs3KyiIjI2PUipuzPclYw58dHR1cunQprRYmOEpLS5kyZUpo8YXIaDLd3mitDRA8bmrUI6estaURrt8dTcestW8D90XzjIhMXm+++Sbbtm0jMzOT/fv3k5OTw9q1a8f1zu3bt/PSSy9hrXUdAKMdKnUUFxdz5coVTpw4wZo1a6J+HuDo0aOsWrWKnJxImwDAmjVrOH78eEzvjza4GWPIzc0dM7jddtttY/47O3HiBJB+CxMg+O+ptLSU+vr6qP5fk8kpmoqbiEhS6urq4uTJk1RXV9Pd3c2JEydYt27duHfX37ZtG+3t7VFt1RBrcCsqKsIYw969e7HWRv08wKlTp1ixYsWo96xZs4a6urrQmarRiHaOGzBqcBsYGODChQuu57dBegY3gOXLl3Px4kWuXr2a6K5IklNwE5GU984772Ctpbq6mkOHDjE4OBhxZWU0op3nZq2NaY4bBFdg5uTk0NraSn19fdTPd3d309TUxPLly0e9b/Xq1VhrQxWsaEQ7xw2CX1ek4Hb+/HkGBwddB7eioiJKSkrGvDcVrV+/nr6+vpiroTJ5KLiJSMp7++23McawceNGDh06REVFBTNmzBj3e5ctW8b06dNdB7doD2EfKicnh4yMDPLz89m7d2/Uz58+fRpgzODmDMPGEhCiHSqF0Stu586dIzMzk4ULF475HmdhQroOIzoLFN55550E90SSnYKbiKS8t99+m1WrVtHa2sq1a9fiUm0DyMjIYPPmza5XYTrDj7FW3Pr7+9myZQvnzp2LeoXpqVPBxfxjBbfy8nJycnJCQS8asQa3SKtK6+vrWbhwIZmZo0+3ttZSW1ublgsTHM5K4GPHtP5ORqfgJiIpbXBwkH379lFdXc2RI0coKCigsrIybu+vqqrixIkTruaEjSe45eTk0N/fT1VVFVlZWVGfcHDq1CkyMzNZvHjxqPdlZGRQWVnJmTNnonq/tZa+vr64zXHr7u6mra3N1Wa6zc3NdHV1pe38NoDp06dTXFzMu+++m+iuSJJTcBORlHb69Gm6urpYv349Z8+eZc2aNWRkxO+Ptk2bNhEIBDh69OiY94634ubMIVu/fj3Hjx/n2rVrrp8/deoUixcvJisra8x7lyxZEnVAcE51iNdQqbMNSEVFxZjvSPeFCY7FixfT2NiI3+9PdFckiSm4iUhKc3abLywsxFrL+vXr4/r+qqoqAFfDpc7JB+OpuEHwhIPBwcGoTjg4derUmMOkjqVLl3Lu3LmoAoITvuIV3M6dO0dubi5z584d8x2TJbitWrWK9vZ2Ll68mOiuSBJTcBORlLZv3z6mT5/O5cuXmT9/PjNnzozr+0tKSigpKeHgwbGPZHYqZNOmTYu6HWeOGwT3dFu2bBmHDh1iYGBgzGcHBgaoq6tzHdyWLFnCjRs3Qoe7uxFrcAu3qtRaS11dHZWVla4WG9TW1lJSUsL06dOjajvVVFVVEQgEYj7ZQiYHBTcRSWn79+9nzZo1XL58Oe7VNkdVVZWrb6bOHlxTp06Nug2n4ubs4VZdXU1vb6+ryep1dXUEAoGoghsQ1Tw3J3zFY45bc3MzPT09Y87Hc6TrUVcjOWfJ6ugrGY2Cm4ikrJ6eHk6cOMGcOXPIzMxk5cqVnrSzadMmzpw5M+bmqE7FLdbgBoQqbAsWLKCkpCS0R91ozp07B+B6UcbSpUsBoprnFutWJ/n5+bcs7Dh79izGGFfz2wKBACdPnpwUwW3FihUYY2LaY08mDwU3EUlZx44dIxAIkJ2dzcqVK0c96mk8nO1FxqqEjCe4OZUsJyAZY6iuruby5ctjntzgDHm6WaEJMGPGDKZPnx5VcIt1qDQ/Pz80989RV1fH/PnzXc0FfO+99+jr65sUwW3KlCmUlJRQV1cX8+kZkv4U3EQkZTnzzmbNmsW6des8a2fjxo3D2otkPHPcnNDpzHOD4J5s06ZNCy3AiKS+vp68vDxmzZrlur2lS5fGNFQabXArKCigu7s7FES6u7tpaWlxPUw6WRYmOJYuXaqjr2RUCm4ikrIOHjxIUVERixYtYtGiRZ61M3PmTMrKysac5+Z8s83Pz4+6jZEVNwCfz8fmzZtpaGjg0qVLEZ9taGigtLQ0qlMFot0SJNY5bvn5+aE94CBYbQP3w7q1tbUYY8Y8gzVdrFu3jitXroS2SxEZScFNRFLWO++8w+zZs1m7dq3nRyG5WaBw7do1CgoKYtpHLlzFDYLVvrE25K2vr6e0tDSq9iorK2lpabllGDOS8cxxg2ClDYLBraCgwNU2IBBcmFBeXh5TGE5FW7ZsAYhqKxiZXBTcRCQldXd38+6771JSUuLpMKlj06ZNNDQ00NHREfGea9euxTS/DT6oZI0MblOmTGH9+vXU1NRE3JC3vr7e9fw2h3Owu9stQcYzVArBPe4GBwepq6tj8eLFroN2bW2tZ4tOkpFz9JWbDZ9lclJwE5GUdPjwYay1bNy4kcLCQs/bczbiHW2BwniCm1NxC3eu52gb8nZ2dtLV1RVzcHM7JOdUzJwg5tbQitv58+fp7+8PbUcylr6+Ps6cOcOaNWuiajOVOWfJOmfPioyk4CYiKen3v/89AA899NCEtOdUQkZboNDZ2UlRUVFM7480VAqjb8hbX18PEPVQqRP0nOfH4gypRjtkObTi5pyn6mYbEIATJ04wODjI2rVro2ozlWVkZFBWVkZjY6OrzZdl8lFwE5GUtGfPHgoLC9m+ffuEtFdYWMiSJUtGDW5XrlyhuLg4pveHW5wwVKQNeaPdCsQxa9Ys8vPzJ6zidu3aNU6fPs3ixYvJzs529ezx48cBJlVwg+B+bm1tbbS1tSW6K5KEFNxEJOX09fVx6tQpVqxYQWZm5oS1W1VVNepQ6XiC22gVN4i8IW+sFTdjDOXl5VEFN5/P5zp0OZyg19TUxLVr11i2bJnrZ48dO0ZeXl5oWHey2LBhQ6hCKTKSgpuIpJx9+/bR0dHBnXfeOaHtbty4kaamJlpbW8N+7mXFLdKGvA0NDUybNi2mdsvKyqIaKi0oKIh69a5TcXNOS3A7vw2CFbfVq1fj8/miajPVbd68GdDKUglPwU1EUs6//uu/AnD33XdPaLujLVCw1no2x80RbkNeZyuQWLZDcSpubnbp7+7ujmlLDqfi5vTT7apUay3Hjh2bVAsTHM7X7Gw+LDKUgpuIpJS2trbQ3CcnSE2U9evXY4wJO8/t2rVrBAIBzypuEH5D3oaGhqjntznKy8vp6elxNZfKqbhFywl7V65ciWqYtLm5mffff3/SzW8DmDNnDoWFhZw9e1ZHX8ktFNxEJKUcPXqUixcvsmjRImbOnDmhbU+dOpWlS5eGrbh1dnYCeDbHzTF0Q15rbUyb7zqiWVna3d09ruA2MDAQVXCbrAsTHEuWLOHixYuh/69EHApuIpIyAoEAx48fp62tLXTw+0SrqqoKW3F7//33Ae+D29ANeRsaGrh+/fq4Km7gbi+369evxzRUmp2djc/nIycnJ6ozXJ3Vs6tXr466zXSwZs0a2tvbuXjxYqK7IklGwU1EUsbZs2dpb2+nra1twodJHVVVVbS0tNDS0jLsujPcOHv27Jje62ao1OFsyPvSSy8B0a8odTjPuQlusVbc2trayMrKivrZ48ePU1paOiGbKyejTZs24ff7OXz4cKK7IklGwU1EUsbRo0e5cuUKEBwyTIRICxSclaZz5syJ6b1uK27wwYa8e/fuBaLfw82Rl5fH3LlzXQ2Vxlpxq62tJTs7O/T1uTVZFyY43JzUIZOTgpuIpATnbNIbN24AiQtu69atIyMjI+7BzefzkZmZ6ariBsENeZ0qX6wVN8D1Xm6xVNystdTW1pKfnx/VKQC9vb2cOXNm0s5vA1i5ciU+n4+TJ08muiuSZBTcRCQlHD9+HGstHR0dlJeXxzyXbLzy8/NZvnz5LfPcWltbo57HNV+FXpwAACAASURBVFJOTo7r4LZgwQL8fj95eXkxn48K0QW3aCtuzc3NXLlyhaKiotDJC26cPHly0h11NdKUKVNYtGgR9fX1rqqwMnkouIlI0rPWcvToUebPn09tbW3Cqm0OZ4HC0K0aWltbmTNnTkz7qTny8vLo7e11da8xhhs3blBYWMiZM2dibrOsrIwLFy6MWRGLZTuQmpoafD4fM2fODJ116oazMGEyD5VCcGHGpUuXdPSVDKPgJiJJr6Wlhfb29tBO/86B74myceNGWltbaW5uDl27ePEit91227jem5eXR09Pj+v7Ozo6mD17Nnv27Il5v6/y8nIGBwdpbGyMeI/f72dgYCCqituNGzeoqalh6dKlTJ06NaqK2/Hjx8nLy3N9GH262rRpE9euXdNwqQyj4CYiSe/IkSPD5n8luuLmHEk09ASDhoaGcc01g+iCm7WW8+fPs2bNGi5evEhdXV1MbbrZy82plkVTcTtz5gy9vb2sX7+egoKCqCtuq1evJiNjcn+L2rp1K6Cjr2S4yf27QkSSnt/vp7a2lhUrVlBTUwOQ8Irbhg0byM/PZ/fu3UBwf7nz58+PO7jl5+e7DjiXLl2ir6+PTZs2UVhYyO7du2OqujnBbbR5bk61LJrgduTIEaZNm0Z5eTkFBQVcu3bN1XPOUVeTeX6bY926dcAHmxGLgIKbiCS5U6dO0d/fz/r16zl8+DClpaXMmDEjoX3Kyspi27Zt7Nq1CwgOk/r9/pi35XBEU3FraGgAoKKigu3bt9Pc3OxqkcFI8+bNIysry1XFze1QaWdnJ+fOnQutwC0sLKSrq8vVs86CBgU3mDFjBrNmzdLRVzKMgpuIJLWjR49SVFTEokWLOHToUMKHSR133303J06coL29ndOnTwOMe05WNBU3J2iVlpaybt06pk2bFtNcN5/PF1q9GEm0FbejR48CwbNdIbjv3PXr1/H7/WM+q4UJwy1fvpyWlpbQ/oUiCm4ikrQ6Ozupr69n3bp1dHV1ce7cuYQPkzruvfdeAH7/+9+H9nQbb9+iqbgNDW6ZmZls27aNxsbGUCUuGmNtCeIENzcVt0AgwKFDh6ioqKCoqAgg9KObczcn+1FXI23YsIGOjo6Y/rtKelJwE5Gk5VRu1q1bFzr6J1kqbps2baKkpIRnn32Wd955h7KyMqZPnz6ud0ZTcWtoaGD27Nnk5eUBwW/w06ZNY+fOnVFX3ZzVupFEszjh5MmTdHd3s2XLltA1Z889N8Ht0KFDLF68eNIedTVSdXU11lreeeedRHdFkoSCm4gkJWeSenl5OYWFhUkX3DIyMvjMZz7D888/z29+8xseeOCBcb8z2orb0Dl1mZmZ3HPPPbS0tHDq1Kmo2i0rK+Py5ctcvXo17OduK27WWvbt28eMGTNYvHhx6LpTcXMz3Hf48OGk+W+cDJyjr3RmqTgU3EQkKdXX19PZ2RmaJ3Xo0CEWLlzIzJkzE9yzD3z961+ntLSUoqIivvKVr4z7fdFW3EauYl2zZg2zZs3i9ddfJxAIuG63vLwciLwliBPoxjoV4sKFC7S0tLBly5ZhGxE7Fbexgtv7779PQ0ND0gyHJ4PS0lJyc3OjDuOSvhTcRCQpHT58mNzcXJYtWwaQVAsTHLNmzaKmpoaGhoZQP8fDqbiNNdQZCARobGy8ZRVrRkYG9913H5cvX+bIkSOu2x1rLzdnRehYw5f79u0jJyfnlhWhbodKnT4ruH0gIyODyspKGhsbXR+HJulNwU1Ekk5PTw+nT59m9erVZGZm0tXVxdmzZ5PyG3pBQUHc5mM5Q5FjHXvV0tIScfuRJUuWsHDhQnbt2uX6G/1Ywa2zsxNjzKhnora1tXHy5Ek2bdpEdnb2sM/cDpU6izycKqsErV27lkuXLg07qUMmLwU3EUk6NTU1BAKBUFBzFikkW8Ut3pyFBmPNcxu6onQkYwwf+chH6Onp4d/+7d9ctTt9+nSmTZs2asVt2rRpo55ksGfPHrKzs6murr7lM7cVt2TZpy/ZVFdX4/f7OXDgQKK7IklAwU1Ekoq1lsOHD1NSUsKcOXOADyoxkyW4jTXPzdkaItKGv7fddhtVVVUcOHCAixcvjtmuMYaysrKIW4J0dXWFqmbhtLW1ceLECbZs2RL6GobKzc0lJydnzIrb4cOHk7KqmmjOCl0dfSWg4CYiSaalpYW2trZhw2WHDh1i/vz5zJ49O4E9854zVOqm4maMYeHChRHvuffee8nNzeWll15ytT3IaFuCdHZ2jjocvGvXLrKzs0Nna4ZTVFQ0asUtmYfDE23FihX4fD5qa2sT3RVJAgpuIpJUnAPlV61aFbqWjAsTvOC24lZfX09JSQk5OTkR75kyZQr3338/Fy5cCFUsR1NeXk59fX3YkDdaxa2hoYFTp05x++23h622OYqLi0etuDnD4Qput5oyZQqlpaW89957rlcdS/pScBORpDEwMEBtbS0rV65kypQpAFy7do133313UnxDd1txC7cVSDhr166lvLyc3//+97z//vuj3ltWVkZvby+tra23fNbV1RW24jY4OMjLL79MYWEht99++6jvH6vi5uxTNhn+O8di48aNtLS00NLSkuiuSIIpuIlI0jh58mToQHnH0aNHsdZOqoqbm6FSNwfaG2P4gz/4A3w+H7/5zW8YHByMeO9oK0uvXLkStuK2f/9+Wltb2bFjB1lZWaP2ZayK26FDh5g3b15oXqMMd8cdd9D7/7d35+FRVOnix78nK2tYQlgjDEgIyI5hE1nCLjCIjAiiIzBcRkfHbX5zvTOjjzDXq85cx+ssjIjL6IyICihEBQQStrCEkCBbIIQlCYQlkBBCEsh+fn90d6YTutOdpNOVTr+f5+knnao6p97T1V39dtWpU7dvO3X0VDRukrgJIRqMH374geDg4Ep9t7zlwgT49xG36k6HlZSUkJGR4dQRNzANmjtt2jQyMjKqvcq0ukF4s7Oz7xj4+Pr168TExBAWFkbfvn0dxtGmTRuHR9zkaJt9lv6DcXFxBkcijCaJmxCiQbh69Srnz59n8ODBlUbdT0xMpHPnznTs2NHA6NzDmSNuFy5coLy83Kkjbhb9+/dn8ODB7NmzhxMnTthcxpIIVr2ytLCwkPz8/EqJW1lZGRs2bMDX15cZM2ZU2l72tG3bluzsbJvz8vLySE5O9orkvLYsYxoePXrU6FCEwSRxE0I0CAkJCfj6+t4x+GpiYqLXHIlx5oibZSgQZ4+4WUybNo3Q0FA2bNhgc4iQpk2b0rlzZ86cOVNpuiXZsk7ctm3bxoULF5gxY4bD22BZdOjQgZycHIqKiu6YFx8fj9a62qtSvV1gYCA9e/YkLS3N7j1lhXeQxE0IYbji4mKOHDlC3759K12ZmJubS3JyMsOGDTMwOvdx5oib5VRmTY64gekm9I888ghNmzZl1apVXLt27Y5lwsPDOXXqVKVplsTNMihuQkICBw4cYPjw4ZWu/HXEcsT06tWrd8yznP7zlu1cWxEREVy6dImMjAyjQxEGksRNCGG4Y8eOUVxcTERERKXpiYmJaK295gvdmeFA0tLS8PX15a677qpx/S1btuSJJ57Ax8eHf/7zn3fcQsmSuFkPCZKVlQWYjrgdOXKEjRs3EhYWxqRJk2q0bstFB1euXLljXlxcHH369Km4w4KwbfTo0RQVFZGQkGB0KMJAkrgJIQyltebgwYN06NCB0NDQSvMsI8UPHTrUiNDczt/fH39/f4dH3EJDQ/Hz86vVOoKDg1mwYAH+/v588sknFVftgilxy8nJqUjW4N+J27lz59iwYQPdu3dnzpw5+Pr61mi9liNuVRM3rTVxcXFymtQJljsoHDhwwOBIhJEkcRNCGCojI4PMzEyGDh16Ryf3AwcOEBYWRtu2bQ2Kzv2aNWtW7RG31NTUGvdvq6pdu3YsXryYzp07ExUVxerVq7l06RLh4eEAFadLtdYkJycDcPLkSQYNGsT8+fMdDv1hiyVxq9q/7uzZs2RlZUni5oS+ffsSEBDAsWPHnLobhmicaveTTQghXCQhIYGAgAD69+9/x7z4+HgiIyMNiMo4QUFB5OXl2Z1/7tw5HnjggTqvp0WLFixYsIADBw6wa9cuPvjgg4rEecOGDeTk5HDhwgW2bt2Kr68v8+fPZ+DAgU5dQWpLp06d8PX15cKFC5WmW/q3SeLmmJ+fH7179yY9PZ2cnByv+kEj/k2OuAkhDJOfn09SUhIDBw4kICCg0ryLFy9y6dIlr+nfZtGqVStyc3Ntzrt16xZXrlzh7rvvdsm6fHx8GDlyJC+++CKTJ0+ma9eu+Pr6sn//fk6ePEnz5s0JCgoiNDSUQYMG1TppA1PSERoaWnFVrEVsbCxBQUFOjQUnTN0GLl++THp6utGhCINI4iaEMMzBgwcpKyur6LtjzdK/zda8xiwoKMhu4ma5otQyWK6rBAYGMnLkSBYuXMiAAQNo3rw5L730EgsWLKCoqIguXbq4ZD0/+tGP7kjctm/fztixY2vcZ85b3X///ZSUlEg/Ny8miZsQwhAlJSUkJCTQq1eviqEmrB04cAB/f38GDhxoQHTGadWqld1xus6ePQvgsiNutgwZMoRDhw5V9KG6ePEinTt3dkndP/rRjyrdmeH8+fOcOXOG8ePHu6R+byAXKAhJ3IQQhjh27Bi3bt1i5MiRNufHx8czcODAipvNe4vqjrhZ7mrg6iNu1oYMGUJ2djYXLlygrKzM6fuiOiM8PJyLFy9WtM9yCy5J3JzXq1cvmjZtyokTJyguLjY6HGEASdyEEG6ntWb//v107NiRbt263TG/rKyMhIQEr+vfBo6PuAUFBdk8QukqlrtUHDp0iPT0dIqLiyuuNq0rywUox48fB+Cbb76hY8eONRrI19v5+vrSv39/MjIyuHTpktHhCAM4nbgppXyUUi8qpZKVUoVKqQtKqbeVUs1dXV4ptVMppe08ImzVL4TwHJYhIEaOHGmzw3tycjJ5eXle178Nqr844dy5c/To0aNOFwk4MmDAAHx8fDh06BApKSkALkvcBgwYAMAPP/zAzZs32bhxI3PmzMHHR44h1MSYMWO4cuXKHbcnE96hJsOBvAM8B6wH3gb6mP8frJSaqLUud3H5LOBFG/WcszFNCOFB9u/fT4sWLexeSWjpv+MtA+9aCwoKorCwkOLi4juutD179my9H51q1qwZgwcPZvv27RWnqe+55x6X1H3XXXfRrVs3oqOjCQoKoqioiHnz5rmkbm8yZswY/vSnPxEbG8vEiRONDke4mVOJm1KqL/As8LXW+idW01OBvwLzgNUuLl+gtV7lZDuEEB7i0qVLnDt3jvHjx9u9knDfvn20adPGZUd6PEmrVq0A031aQ0JCKqaXl5eTmprKzJkz6z2GBx54gDfeeIPc3FwGDBjgsvHClFJMnz6dd999l507dxIWFibjt9WC5TVLSEhAa12vR2BFw+Ps8elHAQX8ucr0D4BbwOP1Ud58ejVIybtSiEYjNjaWJk2aVHs0be/evdx3331eeQrN0n/t+vXrlaZfvHiR4uLier2i1OKJJ56gvLyc48ePM336dJfW/etf/5oWLVqQn5/Pn//8Z6/cxnUVEhLCXXfdxblz58jJyTE6HOFmzn5ihgLlQLz1RK11IXDYPN/V5bsA+UAukK+U+lop1dvJeIUQDVBmZibJyckMGzbM7tWi2dnZJCcnM2rUKDdH1zBYErfs7OxK091xRalFWFgY7733HrNnz+Y3v/mNS+vu3r07KSkpnD17lmnTprm0bm8yfPhwLly4wPnz540ORbiZs4lbZyBLa11kY95FoJ1SKsDGvNqWTwX+F1gEzAHeBR4ADiil7rwvjhDCI8TGxhIQEFDt6bF9+/YBpoFGvZElcbO+0TvA6dOnAejZs6db4njyySf56quvCAoKcnndnTp1snk1sXDeuHHjKCgoICEhwehQhJs5m7g1A2wlXQCFVsu4pLzWepHW+mWt9Zda63Va6/8EJgMtgP+ztxKl1M+VUglKqYRr165VE44Qwt2uXbtGUlISQ4cOpWnTpnaX27t3L/7+/kREeOcF5O3atQPuPOKWnJxMkyZNJOERwL9/2Fh+6Ajv4WzidgsItDOvidUy9VUerXUssBuIVErZ3Otrrd/XWkdorSOsO/UKIYy3a9cu/P397Q64a7F3717uvffeapO7xszeqdKTJ08SHh4ufcIEAP369aNp06YkJSVx+/Zto8MRbuTsHuASptOZtpKvLphOg1Y3hHNdy1ukAb5AGyeWFUI0EJcuXSIpKYkRI0bQvLn9oR+Lioo4ePCg1/ZvA2jRogX+/v53nCpNTk6md2/p5itMfH19GTx4sPRz80LOJm4HzctWGsZcKdUEGAQ4Osle1/IWYUApcN3RgkKIhkFrTXR0NM2aNXOYkB06dIiioiKvTtyUUoSEhGDd3aOwsJDU1FRJ3EQlY8eOrbjgR3gPZxO3LwENvFBl+hJMfdM+s0xQSt1t4+rPmpRvpZS6Y3AnpdR0YBSwzXw1qhDCA5w9e5bU1FTGjBlDYKC9HhMmsbGxANx3333uCK3B6ty5c6XbGZ0+fRqttSRuopL777+f8vJydu/ebXQowo2cGoBXa31MKfV34JdKqa+BTfz7zge7qDx4bgzQDdO4bbUpHwn8n1LqW0x3SSjFdKTucUx3U6ia/AkhGqjy8nKio6Np06aNUxcb7Nixg3vuuYcOHTq4IbqGq3PnzqSmplb8f+LECQBJ3EQlo0aNqrg9WVFRkcMfRqJxqEkv1xeAXwN9gb9jutvB34AZTtzuqiblTwGJwAzgdUxXkd4PvAcM0lqn1CBmIYSBEhMTyczMZMKECXbvkmBRUlJCbGwskZGRboqu4erSpUulI26HDh0iICDAZbeeEo1Dq1at6NOnD2lpaVy4cMHocISbOH2vUq11GaZ7jL7tYLkf1bH8SUxjtwkhPFh+fj4xMTH06NHDqYTj4MGDFBQUSOKG6YhbdnY2hYWFNGnShMTERPr373/HvUuFmDBhAu+++y6nT5922xh/wlhyXbkQol5s27aN0tJSpk2b5tS9FHfs2AGYOlx7u65duwKQnp6O1prExETuvfdeg6MSDVFkZCSlpaXSz82LOH3ETQghnHXmzBmOHj3K6NGjK8Ylc2THjh0MHDiwYgBab2bpy5acnIyfnx83btyQxE3YNHr0aMB0xLqkpAR/f3+DIxL1TRI3IRqovLw8Ll++TE5ODiUlJQA0bdqU4OBg2rdvT7Nm1d2sxDi3bt0iKiqK9u3bM2bMGKfKFBUVsXfvXp566ql6js4zhIeHA6bE7dYt09jkkrgJW4KDg+nVqxepqalkZGTQvXt3o0MS9UwSNyEakLy8PI4cOcLRo0dxdNu2jh07EhYWxsCBA50+quUOmzZt4tatWzz22GP4+Tm3i4mLi6OwsFD6t5m1atWK0NBQfvjhB5KTk2nTpg2DBg0yOizRQE2YMIGPPvqIlJQUSdy8gCRuQjQA+fn5xMbGkpiYSFlZGV27dmXSpEmEhobStm3bisv8CwoKuH79OhkZGZw7d449e/YQGxtLt27duO+++wgLC3OqP1l9SUxMJCkpifHjx9OxY0eny0VHR+Pj4+P0ETpvEBkZyddff01AQABTpkxxeFWu8F7jx49nxYoV7Ny5kylTphgdjqhnkrgJYSCtNYcPH2bLli0UFxczePBgRo0aRdu2bW0u37p1a1q3bk2PHj0YM2YMeXl5HD16lIMHD/L5558TEhLCuHHj6NOnj9sTuAsXLrBp0yZ69uxZ4zsfbNq0ifvuu4/WrVvXU3Se56GHHuLTTz+loKBATiGLall+8CQkJHD79m2vvc+vt5DETQiD3L59mw0bNpCSkkK3bt2YMWNGjTvmt2zZklGjRjFixAiSkpKIjY1l7dq1hIaGMmnSpIqrE+tbbm4ua9asoVWrVsyePbtGN0K/fPkyhw4d4s0336zHCD3PrFmzeOONN2jXrp1caSuq1b59e3r27ElaWhqpqaky3l8jJ4mbEAbIysri888/58aNG0yZMoXhw4fX6QiZr68vAwYMoF+/fhw+fJidO3fy8ccfEx4ezsSJE+v1Ss28vDz+9a9/UVJSwk9/+tMa/9r//vvvAZg2bVp9hOexlFL89re/NToM4SGmTp3KypUrSU5OlsStkZNx3IRws4yMDD788EMKCwtZsGABI0aMcNlpTR8fH4YMGcKzzz5LZGQkqampvPvuu2zatImCggKXrMNafn4+q1atIi8vj8cee4z27dvXuI6NGzcSGhpK//79XR6fEN5iypQplJSUsH37dqNDEfVMjrgJ4Ubp6emsXr2a5s2b88QTT9Rbny5/f3/GjBnDvffey86dO0lISODIkSOMHj2a4cOHu2Ssp6ysLD777DPy8/OZP38+d911V43rKCkpYevWrTz66KOGXlQhhKcbN24cfn5+HD58mJycHNq0aWN0SKKeyBE3IdwkNTWVVatWERQUxKJFi9zSEb958+ZMnz6dp59+mu7duxMTE8Py5cs5fPgwZWVltapTa82RI0f44IMPKCkpYeHChbUegmDv3r3k5eXJaVIh6qhFixYMGzaMs2fPcvbsWaPDEfVIEjch3ODy5ct88cUXtG3bloULF9KyZUu3rr9du3bMmzePBQsW0Lx5c6KiovjrX//Kvn37uH37ttP1XL16lc8++4wNGzbQsWNHlixZQpcuXWodV1RUFAEBAUyYMKHWdQghTKZNm8aVK1c4dOiQ0aGIeqS01kbHUC8iIiJ0QkKC0WEIwfXr1/nHP/6Bn58fixcvdnvSVpXWmtOnT7N//37S0tLw8fGhZ8+ehIWF0bVrV4KDgyvGDNNac+PGDdLS0jh+/Djnzp0jMDCQcePGMWzYsBpdPVpVeXk53bp1Y8iQIURFRbmqeUJ4rcTERCIiInjkkUdYvXq1jP3nYZRSiVrrCEfLSR83IerR7du3+eyzzygvL+fxxx83PGkD09WKvXr1olevXly5coWjR4+SlJRESkpKxTLNmjXDx8eHwsJCSktLAdNo/pGRkURERLjkdlsHDhwgIyODN954o851CSFg8ODBtG7dmlOnTpGenk6PHj2MDknUA0nchKgn5eXlrFu3jtzcXBYsWNAgb57esWNHOnbsyKRJk8jJyeH8+fPcuHGDvLw8tNYEBgbSrl07QkNDad++vUsvIFi7di0BAQHMnDnTZXUK4c18fHyYPHkymzZt4tSpU5K4NVKSuAlRT6Kjozl37hwzZ86s1RWX7qSUom3btnbv2OBqlqR28uTJtGrVyi3rFMIbTJ06lTVr1rBz506mTp0qV2s3QnJxghD14OjRo+zfv5+hQ4cyePBgl9d/+/Ztrl27RkZGRsXRMU+ya9cuLly4wLx584wORYhGxXKv0oSEBLKysgyORtQHOeImhItdvXqVb7/9lm7durnkhs9lZWXExcWxbds2YmNjOXnyJJcvX660TGBgIH369GHQoEFMnDiRqVOnEhwcXOd115ePP/6YoKAgHnroIaNDEaJR6dy5M0OGDOHUqVOkpKQQEhJidEjCxSRxE8KFSktL+eqrrwgMDOThhx+u01VdFy9eZMWKFfzrX//iwoULKKUYPHgwU6ZMoWfPnrRq1YqAgAByc3PJzMzk+PHjfPfdd3zyySf4+PjwwAMP8POf/5zp06c3qKvLbt68ybp16/jpT3/qkoschBCVzZ49m1deeYX4+HhGjRpldDjCxSRxE8KFtm3bxtWrV5k/fz4tWrSoVR3nz5/nD3/4Ax999BGlpaVMmTKFP/7xj0ydOtXhaOjl5eUkJCSwfv16/vnPf/Lggw8SHh7Oq6++yty5cxtEAvfll19y+/ZtFi1aZHQoQjRKM2fO5JVXXmHHjh08+eST8gOpkZE+bkK4yOnTp4mPj2fYsGGEhYXVuHxRURGvvfYa4eHhfPjhhyxcuJAzZ86wadMmHn30UaduYePj48OwYcN48803OX/+PGvWrCEgIIDHHnuMAQMGEBMTU5umuYzWmr/97W/079+f4cOHGxqLEI1Vv3796Nq1K8nJyZWG+RGNgyRuQrhAfn4+UVFRtG/fnkmTJtW4fGxsLP379+fVV1/lxz/+MWfOnGHlypW1vpUUgJ+fH3PmzOHw4cOsXbuWwsJCJk6cyLx58+7oI+cuMTExHDt2jBdeeEGudhOiniilmDVrFqmpqXIXhUZIEjch6khrzTfffENhYSE/+clP8PNzvgdCSUkJr7zyCuPGjaO8vJwtW7awZs0aunbt6rL4fHx8ePjhh0lKSuL3v/89UVFR9O3bly+++MJl63DWO++8Q/v27Zk/f77b1y2EN5k1axalpaVs3bq1Rre1Ew2fJG5C1FF8fDynT59m8uTJtG/f3ulyGRkZjB49mtdff52FCxdy+PBhJk+eXG9xNmnShFdffZXDhw/Tq1cvHn30UebOnUt2dna9rdPa4cOH2bRpE08//TRNmjRxyzqF8FajR48mODiY48ePk5ycbHQ4woUkcROiDjIzM9m2bRthYWEMHTrU6XJ79+4lIiKCpKQk1qxZw0cffVTrixlqKjw8nD179vD666+zfv16t/V9e/nll2ndujXPP/98va9LCG9n6Spx+vRpEhMTjQ5HuJAkbkLUUklJCV9//TVNmjThwQcfdLrP1gcffEBkZCQtW7YkLi6OOXPm1HOkd/Lz8+N3v/sdBw4cICgoiEmTJvHSSy9RXFxcL+vbs2cPmzZt4r/+679o3bp1vaxDCFHZvHnzKC4uZuvWrdy6dcvocISLSOImRC1FR0dz9epVZs2aRfPmzR0ur7Xm5Zdf5uc//zmRkZHEx8fTt29fN0Rq3+DBg0lMTOTJJ5/krbfeYsSIES4/rVJSUsKzzz5Lp06dePbZZ11atxDCvvvvv58OHTpw7NgxTpw4YXQ4wkUkcROiFlJSUoiPPiG2VQAAHahJREFUj2f48OH07NnT4fKlpaUsXryYN954gyVLlrBx40anhvdwh2bNmrFixQqioqK4cOECQ4YMYeXKlS67jdbbb7/N4cOHWb58uVMJrhDCNXx9fZk7dy5nzpwhLi7O6HCEi0jiJkQNWYb+6NChAxMnTnS4fEFBAbNmzeLjjz9m6dKlrFy5skZXnrrLzJkzOXr0KKNHj+app57ioYceqvO9Dg8dOsSyZcuYPXs2s2fPdlGkQghnzZs3j9LSUmJiYuTepY2EJG5C1IDWmqioKIqLi5k9e7bDBCwvL4+pU6eyefNm3nvvPZYtW9agxy/r1KkTmzdv5p133mHz5s3079+fbdu21aquzMxMHnzwQdq3b8+KFStcHKkQwhkjRoygR48eHD58mCNHjhgdjnABSdyEqIEDBw5w5swZJk2a5HDoj9zcXCZPnkxcXBxffPEFTz75pJuirBsfHx9eeOEF4uPjadu2LZMnT+ZXv/oV+fn5TteRlZXF1KlTyc7OrhiYWAjhfkopFi9eTFpaGjExMZSXlxsdkqgjSdyEcFJmZibR0dH06tXL4dAfOTk5TJw4kcTERNauXWvIlaN1NXDgQBISEnjmmWd455136NmzJytWrHA4mOfBgwcZNmwYycnJfP311wwePNhNEQshbHniiSfw8fFh7969nDt3zuhwRB1J4iaEE0pKSvjqq69o0qQJM2fOrPZ0Z1ZWFuPHj+fo0aOsX7+eWbNmuTFS12ratCnLly9n//79hIWF8fTTTxMaGsqLL75IdHQ0169fR2vNjRs3iImJ4fHHH2f48OEUFxezY8cOpk6danQThPB6oaGhTJ48mSNHjhAfH290OKKOJHETwglbtmzh2rVrDof+uHr1KuPHjyc5OZlvvvmG6dOnuzHK+jNixAh2797N9u3bmTBhAsuXL2fSpEkEBwfj5+dHmzZtmDhxIlFRUfzqV78iKSmJESNGGB22EMLsZz/7Gbm5uWzevJkbN24YHY6og4Z3aZsQDczJkydJTExk5MiR1Q79cfnyZSZMmEBaWhrfffcdEyZMcGOU9U8pRWRkJJGRkdy8eZM9e/Zw6tQpsrKyaNu2LeHh4YwbN85td4AQQjhv5syZhISEEB8fT0JCglNXxIuGSRI3IaqRm5vLN998Q6dOnapNxC5evMj48eO5ePEimzdvZuzYsW6M0v2CgoKYNm0a06ZNMzoUIYQTAgMD+cUvfsF///d/s3XrVsaNG9cghyUSjsmpUiHsKC8vZ/369ZSXl/OTn/wEX19fm8udP3+eMWPGcPnyZbZs2dLokzYhhGf6xS9+gb+/P7t37+bYsWNGhyNqSRI3IeyIjY0lPT2dadOmERwcbHOZtLQ0xo4dS3Z2Ntu2bWPUqFFujlIIIZzTsWNH5s6dy5EjR4iJiXHZ3VGEe0niJoQNZ8+eZdeuXQwYMICBAwfaXWbMmDHk5uYSHR3N8OHD3RylEELUzHPPPUdhYSHbtm1z+X2JhXtI4iZEFbm5uXz11Ve0a9fO7lWhKSkpjB07llu3bhETE0NERISboxRCiJobOnQokZGRxMXFsX37djnq5oEkcRPCSmlpKWvWrKGsrIy5c+cSEBBwxzLJycmMGzeOoqIitm/fLgPMCiE8yquvvsrNmzfZuHGjDMjrgSRxE8LK999/z6VLl5g1a5bNfm1JSUmMGzeOsrIydu7cyYABAwyIUggham/s2LHcf//97Nu3jy1btshRNw8jiZsQZocOHSIxMZFRo0bRp0+fO+bHx8czduxYlFLs3LmTvn37GhClEELUjVKKV199ldzcXDZu3MiJEyeMDknUgCRuQgCpqals3LiRu+++m/Hjx98xPzo6mvHjxxMUFERsbKzNxE4IITzFxIkTGTduHLt37+bbb7+lrKzM6JCEkyRxE14vKyuLNWvWEBwczMMPP4yPT+WPxbp165g+fTo9evRg79691d49QQghPIFSirfffpuCggK+/fZbDh48aHRIwkmSuAmvduvWLVavXo2vry/z58+nSZMmlea///77PPLII0RERLBr1y46depkUKRCCOFaQ4YMYcGCBcTHx7Nu3Tpu3rxpdEjCCZK4Ca9VVFTEqlWryMvLY968ebRu3bpiXnl5OS+99BJPPvkkDzzwANu2baNNmzYGRiuEEK73+uuvExgYyDfffMP3339vdDjCCZK4Ca9UUlLC559/TmZmJnPmzCE0NLRiXkFBAQ8//DBvvfUWTz/9NFFRUTRr1szAaIUQon507tyZN998k9OnT/PFF19w6tQpo0MSDkjiJrxOWVkZa9asIT09nYceeohevXpVzMvIyGDMmDFERUXxl7/8heXLl8uNmIUQjdozzzzDyJEj2bp1K6tXryY/P9/okEQ1JHETXsUywO6ZM2f48Y9/TL9+/SrmRUdHM2TIEFJSUoiKiuK5555DKWVgtEIIUf98fHz4xz/+QWlpKZ9//jnr16+Xsd0aMEnchNcoLi5m9erVpKSkMH36dIYMGQKY+rO99tprTJ48mZCQEA4ePMiMGTMMjlYIIdynd+/evPPOO5w+fZpPPvmE/fv3Gx2SsEPOAQmvUFBQwOeff86lS5d46KGHKu54cOXKFRYtWsT333/P/PnzWblyJS1atDA4WiGEcL+nnnqK3bt38+WXX/L+++8TEhJCWFiY0WGJKuSIm2j0rl69yocffkhmZiaPPPJIRdK2bt06+vXrx86dO3n33XdZtWqVJG1CCK+llOL9998nPDyctWvX8ve//53MzEyjwxJVSOImGrWUlBQ++ugjSktLWbRoEb179+bq1as89thjzJkzh+7du/PDDz/wi1/8QvqzCSG8XsuWLfn+++8JCgrik08+Yfny5WRlZRkdlrAiiZtolMrKytiyZQuff/45wcHBLFmyhA4dOvDee+9V/JpctmwZ+/bto3fv3kaHK4QQDUa3bt3YvHkzZWVlrFixgj//+c9cv37d6LCEmfRxE43O1atXiYqK4tKlSwwdOpTJkyezZ88e/vM//5OEhAQiIyN59913JWETQgg7Bg0axNatW5k6dSrLly+nuLiY5557rtKYl8IYcsRNNBqlpaXs2LGDlStXkpOTwyOPPEKnTp2YMWMGkZGRXLlyhVWrVhETEyNJmxBCODBy5Eh27dqFv78/y5cv5+WXX+b48eNGh+X15Iib8Hhaa06cOEFMTAw5OTn079+fZs2a8etf/5pNmzbRpk0b3nrrLZ555hmaNm1qdLhCCOExBg0aRHx8PA8++CCffvop6enpPPfcc8yYMYOAgACjw/NKkrgJj6W1JiUlhdjYWC5evEhQUBABAQH84Q9/4ODBg4SEhPDaa6/xy1/+stJ9SIUQQjive/fuxMXF8dRTT/Hpp59y6tQp9u3bx5IlSwgPDzc6PK+jGuvoyBERETohIcHoMEQ9KCoq4vjx4xw4cIArV66QnZ1NRkYGMTExFBQU0Lt3b55//nkWLFggR9iEEMKFNmzYwJIlS8jOzqZ///4sXLiQhx9+mLvuusvo0DyeUipRax3hcDlnEzellA/wPPAk8CPgGrAGeFVrXeDq8kqpacArwECgCIgBXtJapzoTryRujUtZWRnp6ekcPXqUQ4cOkZKSQnp6OidPniQ3N5fmzZszb948Fi9ezIgRI2RoDyGEqCc3btzg9ddf569//SulpaX06dOH2bNnM3fuXHr37o2vr6/RIXqk+kjc/gI8B6wHNgN9gGeBWGCi1rrcVeWVUrOBdcAR4AOgFfACUAZEaK0vOYpXEjfPprXm+vXrnD9/nri4OPbt28fZs2c5f/48ly9fpry8nNatWzNjxgxmzZrF1KlTad68udFhCyGE18jIyOCdd97h/fffJz8/n+DgYAYPHsyDDz7I5MmT6dGjB35+0iPLWc4mbmitHT6AvkA58FWV6c8CGpjvqvKAP3ARSAdaWE0fhClxe9+ZmO+9914tPENpaam+ePGi3rRpk/7jH/+oFy1apEePHq27d++umzZtqs3vEe3v76/HjBmjly5dqnfs2KGLi4trtJ6lS5e6LGZX1uWKepcuXVpR1tk6arp8bWKqzTKO4rJuqzPrsp5n/Xzs2LGVnluvtzb12SpreW69rqVLl1asz/rhTFus59t77WzNd1S/ozL21lG1XVXLW7+ujjjzHrZVX9X/rWOyrtfWa22rDbbKa611t27dKv1f9f1jL15bddlrZ23eC7bW7eg97Go3b97UK1as0CNGjNBKKQ3oZs2a6X79+un58+frP/3pTzouLk7n5OTo8vJyt8XlaYAE7UR+49QRN6XU/wAvA2O01rFW05sA2cAurfU0V5RXSk0EtmE6hfpalXpigAignda6pLqY5Yibe2mtKSsro7CwkIKCAnJzc8nKyuL69esVjxs3bpCTk8PVq1fJzMwkOzub69evc/PmTfLy8rB+LwYGBtKjRw/uvfdeRo4cyZAhQxg4cGCd+qwppXDm/e7uulxRr+XUsNba6TosyxnZFlvLOIrLuq3OrMt6XnXPLXXaqt+Z+mzFZqst9k7j26vT3vrtvXbW7XDU7urqdOb1t/caVm2ns+9HW/XaW6Zq7NX9X5W95aur395fW+t0FJu9bWjv/WOPvfeqo/dwfbp8+TKbNm3iu+++Y8+ePZXuutCyZUvatm1LSEgIXbp0oXPnznTo0IGQkJCKvyEhIbRp04agoCCaNGmCj4/3jFrm7BE3Z49hDsV0xCzeeqLWulApddg831XlLc/326gnDhgP9AKSnIy9XiQmJrJo0aJK0ywfkqp/7c13ZhlbHzxn669LWXvzysrKKC0tpbS0tOJ5WVkZZWVlTu8kfHx8aNmyJa1bt6Z9+/b069ePu+66i3vuuYeBAwfSu3dvOnbsKP3UhBDCw3Tq1InFixezePFitNZkZmaSmJjIvn37OHLkCOfPn+f8+fP88MMPlJWVVVuXr68vfn5++Pn54e/vj5+fH76+vvj4+FQ8lFIopWz+b3le9bvE1neLo2UmTZrE22+/XctXxbWcTdw6A1la6yIb8y4C9ymlArTWxS4o39lquq1lAbpgI3FTSv0c+DlA165d7TbGFUpLS8nPz7es1zqGO6Y5M93yvLo3ma311GUZZ2Oy/mv5EAUEBODv709AQEDF88DAwIq/LVq0oHXr1rRu3Zq2bdsSHBxMcHAw7dq1o0OHDm7rvLps2TJ+//vf39GepUuXsmzZMsPqckW9VctZl7VXh711ubMt9paprh1jx45l165dNpextS7A7jqc+YzUtr7q5jv6IWJv3cuWLXO43ezV5Wy7HZWp7vWvyZeio/ejrXq7detGenp6tetwJoa6lHdm32l53qpVK3Jzc2tcl731O3ov2KvDVll3U0rRsWNHpk+fzvTp0yvNKy8v5+rVq1y5coXLly9XjBCQnZ1Nfn4+t2/f5tatWxQWFlJcXExRURFFRUWUlpZSXl5u82E5lWh5bvlrUd1BE2u25hUWFrrgFXENZ0+VngX8tdZ3ZENKqX8BPwXaaK1v1LW8Uuoj4GfA3Vrrc1WW/RnwEfCQ1npDdTHLqVJRlaPTDkbV5Yp6rU+NOFuHvdM1ruJMvbaWcRSXrdNA1a3L3imtqs8tdTo6zWSvDlux2WqLvS9ZOVUqp0rtvX/ssfdedfQeFg2TcvJUqbMnj28BgXbmNbFaxhXlLX9tLe/MuoQQQgghGiVnT5VeAu5RSgXaON3ZBdNpUHunSWta/pLV9JM2lgXbp1GFqJbllFdDq8sV9VqXc7YOy3JGtsXWMo7iqq6Mo+Wtn48dO7bS83Hjxjmsq7r6qvvfel1Lly5l586dFeurrv7q5jt6HZx5Xl2dzrz+VdtVdbnq2uls7Nast5O9Za1jqq4uW8tbYrZVvlu3bpX+r/r+cVS3dVl729DR/1XZiteZcsKzueqq0t1a6wdcUV7JVaVCCCGE8DKuPlX6JaaxtF6oMn0J0Az4zGrFdyulete2PLALuAz8h1KqhVW9A4FxwFpHSZsQQgghRGPk1KlSrfUxpdTfgV8qpb4GNmG688FzmBKt1VaLxwDdAFWb8lrrEqXU85iSvVil1AdAEPAipttkyTFgIYQQQnilmtyL4gUgDdNwG9OBLOBvmE5pVnu7q5qW11qvVUrdxnSv0j/x73uV/pfWWvq3CSGEEMIrOX2vUk8jfdyEEEII4Slc3cdNCCGEEEIYTBI3IYQQQggPIYmbEEIIIYSHkMRNCCGEEMJDSOImhBBCCOEhJHETQgghhPAQkrgJIYQQQngISdyEEEIIITyEJG5CCCGEEB5CEjchhBBCCA8hiZsQQgghhIdotPcqVUpdA9LdsKp2QJYb1tMQSdu9lze335vbDt7dfmm793JH+7tprUMcLdRoEzd3UUolOHNT2MZI2u6dbQfvbr83tx28u/3Sdu9sOzSs9supUiGEEEIIDyGJmxBCCCGEh5DEre7eNzoAA0nbvZc3t9+b2w7e3X5pu/dqMO2XPm5CCCGEEB5CjrgJIYQQQngISdyEEEIIITyEJG5CCCGEEB5CEjcblFJPKqU+U0olK6XKlFLVdgRUSoUrpTYopXKUUgVKqVil1PgarrPOddQXpdRCpZR28OjiRD2fVFP+YXe0pTaUUmnVxN2uBvUMV0pFK6XylFI3lVLfK6UG1WfsdaWU6qKU+q1SapdS6rL5vZmklHpLKRVcg3oa7LZXSvkopV40f94LlVIXlFJvK6Wau6O8kZRSvZRS/62UilNKXTO/Nw8rpV6uQft3VrNtG8S4V/ZUE3d+DeqYppTaZ/5sXFdKrVVKda/PuF1BKbXMwT69xIk6Gvy2N++/1iqlzpnjSnOwfJ330/W9r/dzVUWNzG+BYOAHoDkQam9BpdTdwD6gFPhfIBdYAmxRSj2gtY52tDJX1FHPdgM/tTG9E6Z4D2utL9agPlt1xdcmMDdKBl63MT3PmcJKqRHATuAi8Kp58i+BWKXUfVrrY64Ish78GFgGbATewtTeYcALwFyl1DCt9ZUa1NcQt/07wHPAeuBtoI/5/8FKqYla6/J6Lm+knwHPAN8AnwElQCTwP8AjSqkRWuvbTtSTBbxoY/o5VwVaj2K584pBh0kLgFJqNrAOOAL8J9AK02djr1IqQmt9yZWButjXwBkb0wdgasu3TtbT0Lf9G8B14BDQuroFXbGfdsu+XmstjyoP4EeAj/n5d6aXye6ya4AyYJDVtBaYbrd1CvOVuw7WV+c6DHqdfgto4Bknl/+kuteyoT6ANGBnHeuIB24CXaymdTFP22p0G6uJuy/Q0cb0/zBv+z958rY3t68c+KrK9GfN7Ztfn+WNfgARQCsb0//HHP8vnahjJ5BmdFtq2X4NfFLLsv6YvpzTgRZW0weZ9+fvG92+WrZrpfl1md4Ytj3Qw+r58eridcV+2h37ejlVaoPWOk078SvZfCphJqYv9cNW5fOBD4FewND6rsMISimF6df6bUy/1GtUVikVpJTyqPefUspPKRVUi3I9MW3DtdrqyKT5+VpgolKqo+sidR2tdZK2fUTtS/PffjWprwFu+0cBBfy5yvQPgFvA4/Vc3lBa6wStda6NWTXevuZTxkHmfYNHUUoFKKVa1LDYWKAz8KF5fw2AeT++E9MRaX/XRVn/lFLNgHmYEtLva1CuwW57rbVTR/5csZ92176+oew8PdUAIBDYb2NenPmvo6TLFXUYYSzQE9ORhhs1LJtrftxWSm1TSg13eXSuNxzTF3GuUuqGUuqfSqnOTpa1bD9721gB97ogRneydB/IrGG5hrbth2I6YlbpdK3WuhA4jOPPXl3LN1Q13b5dgHxM2zZfKfW1Uqp3vUTmeg9j+mznKaWuKqX+ppRq5UQ5R5/rIEw/vD3JI5ji/lhrXeZkGU/e9tZcsZ92y75e+rjVjeWL21b/Lss0R532XVGHERab/35YgzJXMPUHSgQKgIGY+oPEKqWmaeP78tmThKmdyZg+M+MwnSqcYO7j5agfi6du4+r83vz3n04u31C3fWcgS2tdZGPeReA+pVSA1rq4nso3OEopX0x9c0qB1U4USQX2AkcxnSIcjqlPzwSl1P264fbfBFPCvRZTX68gYBqm2Mea+yNVd5GCs5/rJBfF6g6LMZ0m/YeTy3vytq/KY77PG23ippRqjemLwVl/1Vpfr+Fqmpn/2tppF1ZZpj7rcIqrXhNzPT/BtLPb7WxlWuvfVJm0QSm1GtORiRVAWA1iq5G6tF1rPb3KvC+UUrsxnSL+PaYLSarjtm1sjys/D0qp/wfMwdSHZ7szlRm57R1ohu3tApW3jb3Eq67lG6I/AyOA32mtTzlaWGu9qMqkdUqpbzCdLvw/YJLLI3QRrXXVI77/UkodxXQh0vPYviDJwvDPtSsppcKB+4EYrXWqM2U8edvb4DHf5402ccN09cjSGiy/CtOVJzVxy/w30Ma8JlWWqc86nOWq12Q+0BT4SJt7XtaW1vq0UmoNsFAp1UtrnVKX+qrh0veD1nq1Uup1oGpSZ4s7t7E9Lmm/Uuo/MF1duhHTL+tac+O2r84toL2dec5sm7qWb1CUUq9h2q7va63frG09WutY84+bSKVUU+3clakNxVuYPivTqT5xawifa1eqzVmUO3jwtveY7/NG28fNfIGBqsHD1mXRjlhOkdk69GmZ5miYDFfU4RQXviaLMZ1G+cQVcWG6ahPA6THRaqqe3g9pTsbstm1sjyvar5T6GaZhE7YCP9FaOzVkggNp5r/1tu0duAS0U0rZ2tF2wXQatLqjZXUt32AopZYBrwAfA0+5oMo0wBdo44K63Mb8vr6E4/ek4Z9rV1FK+QFPYPqxtt4FVabhedveY77PG23i5ibHMB0SHWlj3gjz3wQ31OE25kEEhwAb7VxtWBuW02Q17ehutJ44F/NB819721hj6vvVYCmlFmG6UjIamGWnT1dtGL3tD2LaDw6znqiUaoJpWAdHn726lm8QlFJLMR1l+hfwH3U9km4WhukHXk3PZBjKvO1CcfyedPS5vgkYcRS5Nn4MdAA+ddFn2xO3vSv20+7Z1zsaL8TbHzgex20tpk6ZA62mWcZgS8FqDDZMgzP2BtrVtg6jH8By85vvx9Us087czlZW05oDTWwsOxhT4nrC6LbZaUtbO9OfMb8O7zpqu3n6QUw78s5W0zqbp0Ub3U4Hr8FC8/szGmjqYFmP2vZAf6ofh+1xq2l3A71rW76hPjBdiKAxJW0+1SzXybxtm1lNawX42lh2urnOTUa3r5r2BNuZ/pY59pcctN0f0xGWquO4DTR/Xj40uo01eC2+M7e5f2Pe9jgex83p/bSR+3rDX8iG+MD06+MV8yPZ/Ca0/P/LKsv2xPSrIhP4DfA0pjsulAJTqiy70FzXstrWYfDr0sQc50VbH1ir5ZaZ27nQatog4DKmjui/Ap4E3sXUYfMWcL/R7bPTlhcwHRV9C1Oy9jymUwka08UZIY7abp5+H6Yk5ay5zhfMz/OxStgb2gPTGINlQA6mizAer/KY5enbHvibOeavMV0t/DamkfN3YpXIYDr9o2tbviE++PcPkHRMp8qqbt9JVst+Yl52nNW0WZhGyP+L+bPxDKYrjcuAa0Avo9tYTdvfwTRswxuYTg3/GthubmMcVj9SbLXdPH0OpsT9B0z77d9g2o9fwWoA1ob8wJRUlAIHqlnGY7c9pru1WL6/M837Msv/P62yrNP7aVv7uprWUes2Gf2iNsSH1ZvU1iPNxvJ9gCjgBqYvoj3ARBvLLcRG4laTOgx+Xeab43/dwXJ3vKGBjsCnmBLhm5i+2M6bP+i96zPuOrZ5FKbbAZ3HNNhwIXAS+APQ2pm2W80bCcSYP8B5wBZgiNFtdHJbOvV58MRtj6kvzv/DdJeSIkw/TP4Pq6Mo5uXSsJ24OVW+IT4c7Os0VncMwfaXdx9MZwwsX0yWL6y/08ATF+BB82fwovlzXYDpKuffUeUIsa22W82bgSnRu4UpKVgH3G10+2rwOvzO3LYlTrxPPG7bY/oB5fD9bbW8U/tpW/u6mtZR24cyr0QIIYQQQjRwcnGCEEIIIYSHkMRNCCGEEMJDSOImhBBCCOEhJHETQgghhPAQkrgJIYQQQngISdyEEEIIITyEJG5CCCGEEB5CEjchhBBCCA8hiZsQQgghhIf4/1DGWWVkGI05AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as sps\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.special import logsumexp\n",
    "from scipy.special import digamma\n",
    "%matplotlib inline \n",
    "#%matplotlib qt5 #get an extra window for the plots\n",
    "plt.rcParams.update({'figure.figsize': (10.0, 8.0), 'font.size': 18}) #Make the plots bigger\n",
    "\n",
    "\n",
    "def GMM_pdf(x, mus, vars, pis):\n",
    "  # Function for evaluating univariate GMM at values given by vector 'x'\n",
    "  # 'mus', 'vars', 'pis' are vectors of means, variances and weights of individual Gaussian components\n",
    "  # 'x[:,np.newaxis]' makes 'x' column vector to evaluate every coefficient of 'x' w.r.t. every Gaussian\n",
    "  # component (row vectors 'mus' and 'vars') at once\n",
    "  return sps.norm.pdf(x[:,np.newaxis], mus, np.sqrt(vars)).dot(pis) \n",
    "\n",
    "\n",
    "#1. Handcraft some GMM parameter \n",
    "mus_true = [-4.0, 0.0, 4.0, 5]\n",
    "vars_true = [1.0, 1.96, 1.44, 1]\n",
    "pis_true = [0.1, 0.4, 0.2, 0.3]\n",
    "\n",
    "#2. Pregenerated N datapoints from the GMM distribution above \n",
    "x=np.array([ -3.07371088e+00, -4.14348725e+00, -3.01490821e+00, -3.54303388e+00, -3.48708234e+00, -3.59262207e+00,\n",
    " -5.76100178e+00, -5.02726789e+00, 8.28817153e-04, 4.99898450e-01, 2.83745605e-01, 6.71947042e-01,\n",
    " 7.66679495e-01, 6.96995763e-01, 1.33546855e+00, 6.03847649e-01, -1.05992122e+00, 2.91024229e+00,\n",
    " -2.12682520e+00, 8.33533885e-01, 1.77147857e+00, 7.37629536e-01, -1.25040836e+00, 1.87318623e+00,\n",
    " -4.14582880e-01, 5.05680493e-01, 1.63091140e+00, 6.63219549e-01, -3.30841863e-01, -1.21874646e+00,\n",
    " 2.64384057e+00, -4.32674840e-01, -1.79034947e+00, 3.13567565e-01, -5.43911715e-01, 2.28572951e+00,\n",
    " 9.55653291e-01, -5.43582974e-01, -2.73850574e-01, -1.50383720e+00, 1.15925073e-01, 3.92541838e+00,\n",
    " -1.57267817e+00, 4.43581114e-01, -8.11856886e-01, 2.62442641e+00, -4.36298436e-01, -6.72286580e-01,\n",
    " 1.52223784e+00, 1.25033658e+00, 4.88645989e+00, 2.96110183e+00, 4.74249957e+00, 2.17531545e+00,\n",
    " 3.43072143e+00, 3.49315547e+00, 2.51223591e+00, 2.55369053e+00, 2.93122261e+00, 6.40247818e+00,\n",
    " 5.12748233e+00, 4.08277439e+00, 4.96716209e+00, 1.56304959e+00, 4.31869585e+00, 2.07957592e+00,\n",
    " 4.56265393e+00, 3.74342366e+00, 4.36177483e+00, 5.21824922e+00, 4.94100019e+00, 4.70062989e+00,\n",
    " 6.12111884e+00, 6.69125720e+00, 5.03104495e+00, 5.72199065e+00, 4.29367941e+00, 3.72747772e+00,\n",
    " 4.41461701e+00, 5.48741263e+00, 4.56782193e+00, 6.45701533e+00, 5.49118936e+00, 4.25947605e+00,\n",
    " 3.39235348e+00, 4.10586407e+00, 2.76696554e+00, 6.66059909e+00, 6.00107916e+00, 5.92828295e+00,\n",
    " 4.97460855e+00, 2.77746143e+00, 2.99416076e+00, 5.24459233e+00, 6.44720235e+00, 4.71084807e+00,\n",
    " 5.62475093e+00, 3.76422931e+00, 5.79482964e+00, 5.11432194e+00])\n",
    "\n",
    "\n",
    "\n",
    "#3. Choose some initial GMM parameters for EM training\n",
    "C = 6                             # number of GMM components \n",
    "mus_ml = x[:C]                    # We choose few first observations as the initial means\n",
    "vars_ml = np.repeat(np.var(x), C) # Variances for all components are set to the global variance of the training data\n",
    "pis_ml = np.ones(C)/C             # All component weights are set to the same value 1/C\n",
    "\n",
    "\n",
    "#EM algorithm for ML GMM training\n",
    "for _ in range(1000):\n",
    "  #E-step\n",
    "  log_p_xz = sps.norm.logpdf(x[:,np.newaxis], mus_ml, np.sqrt(vars_ml)) + np.log(pis_ml)\n",
    "  log_p_x  = logsumexp(log_p_xz, axis=1, keepdims=True)\n",
    "  gammas_ml = np.exp(log_p_xz - log_p_x)\n",
    "    \n",
    "  #M-step\n",
    "  Nc = gammas_ml.sum(axis=0)\n",
    "  mus_ml  =  x.dot(gammas_ml) / Nc\n",
    "  vars_ml =  (x**2).dot(gammas_ml) / Nc - mus_ml**2\n",
    "  pis_ml  = Nc / Nc.sum()\n",
    "\n",
    "#4. Plot the true GMM, ML trained GMM and the observations\n",
    "t = np.linspace(-10,10,1000)\n",
    "true_GMM_pdf = GMM_pdf(t, mus_true, vars_true, pis_true)\n",
    "ml_GMM_pdf  = GMM_pdf(t, mus_ml, vars_ml, pis_ml);\n",
    "plt.plot(t, true_GMM_pdf, 'gray')\n",
    "plt.plot(t, ml_GMM_pdf, 'k')\n",
    "plt.plot(x, np.zeros_like(x), '+k');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For for the Bayesian GMM, use the following setting of the prior parameters!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters of NormalGamma prior over means (mu) and precision (lambda)\n",
    "m0, kappa0, a0, b0=[0.0, 0.05, 0.05, 0.05] \n",
    "\n",
    "#Parameters of Dirichlet prior weights\n",
    "alpha0=1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variational Bayes inference\n",
    "\n",
    "Now, include code the Variational Bayes (VB) approximate inference in Bayesian GMM.\n",
    "\n",
    "* Use the responsibilities obtained from the ML training (EM algorithm implemented above; variable 'gammas_ml') to initialize $q(\\zz)$. This implies that we use the same number of Gaussian components $C=6$ as for the ML training.\n",
    "\n",
    "* Run iteratively the updates of $q(\\mmu,\\llambda)$, $q(\\ppi)$ and $q(\\zz)$ until the algorithm converges.\n",
    "\n",
    "Once the VB inference converges:\n",
    "\n",
    "* Plot the estimated approximate distributions $q(\\mu_c,\\lambda_c)$ for all the Gaussian components $c \\in \\{1\\dts C\\}$. You can reuse the code for plotting NormalGamma distribution from [vb_gmm_training.ipynb](http://www.fit.vutbr.cz/study/courses/BAYa/public/notebooks/vb_gmm_training.ipynb) notebook.\n",
    "\n",
    "* Print also the vector of parameters $\\aalpha^*$ of the approximate posterior distribution of the weights $q(\\ppi)$.\n",
    "\n",
    "* Analyze these results. What can you say about the obtained posteriors distributions? What do they represent? How do these posterior distribution compare to the parameter estimates obtained from the EM algorithm?\n",
    "\n",
    "* Generate several (say 50) samples from the approximate posterior distribution $q(\\ppi,\\mmu,\\llambda) = q(\\ppi) \\prod_{c=1}^C q(\\mu_c,\\lambda_c)$. Each such sample is set of parameters of one GMM. Plot the GMM distributions corresponding to all the samples into a single figure. Comment on this plot. What do the individual GMM distributions represent?\n",
    "\n",
    "* Now, average all the GMM distributions from the previous step. Attention! Average the GMM distributions, not their parameters as sampled from $q(\\ppi,\\mmu,\\llambda)$. What can you say about the obtained average distribution? What does it represent?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code goes here :)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Now estimate and plot the VB (approximate) posterior predictive distribution.\n",
    "How does the resulting posterior predictive distribution compare to:\n",
    "- the true training data distribution\n",
    "- the GMM obtained using ML training (i.e. using EM algorithm)\n",
    "- the average of GMM distributions obtained in the previous step by sampling from $q(\\ppi,\\mmu,\\llambda)$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code for posterior predictive distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gibbs sampling inference\n",
    "Now repeat the whole excercise, but this time using the inference based on Gibbs sampling (not the collapsed one). \n",
    "\n",
    "* Start from fixed $\\zz$ initialized as the most likely assignments of observations to GMM components according the the ML training (i.e. assignment will be again derived from variable 'gammas_ml').\n",
    "\n",
    "* In each iteration, plot the GMM distribution corresponding to the sample from $p(\\ppi,\\mmu,\\llambda|\\zz, \\xx)$. What do the individual GMM distributions represent? How do they compare to the samples from the VB approximate posterior $q(\\ppi,\\mmu,\\llambda)$?\n",
    "\n",
    "* Average the distributions from the previous step to obtain an approximation to the posterior predictive distribution. Plot this distribution. How does it compare to the VB posterior predictive distribution? Is it any different? If yes, why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code for Gibbs sampling inference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## More training data\n",
    "Generate larger number of training observations (e.g. 1000 samples) from the true distribution (mus_true, vars_true, pis_true) and repeat the whole experiment once more (for both VB and Gibbs sampling inference) with the larger amount of training data. Regenerate all the plots and comment on how they change from the previous experiments with smaller training data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try to implement as much as you can and answer as many questions as you can. It is fine if you do not understand some questions or if you do not know some of the answers."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}