{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "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\n", "plt.rcParams.update({'figure.figsize': (10.0, 8.0), 'font.size': 18})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, We just handcraft some GMM, generate data from it and ML train GMM using standard EM algorithm as a contrastive system" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xl4VdW5+PHvm5kxCSQkQCADCVOEhBmkDAIKIoh1BCesFltbr1r1/mqtFa9Wr/be1qH1qrTibMVZrCgoFhAKSIAwJoEQQshEQkamzOv3R845PcSEnCRnyPB+nuc8OWfvtdd+z85J3rPXXnstMcaglFJKeXk6AKWUUu2DJgSllFKAJgSllFIWmhCUUkoBmhCUUkpZaEJQSikFaEJQSilloQlBKaUUoAlBKaWUhY+nA2iJkJAQExUV5ekwlFKqQ9m5c+dJY0xoc+U6VEKIiooiKSnJ02EopVSHIiLHHCmnTUZKKaUATQhKKaUsNCEopZQCNCEopZSy0ISglFIK0ISglFLKQhOCUkopQBOCUsrF9uzZwx/+8Ae2bNni6VBUMzrUjWlKqY7lgw8+YPHixdTV1QHw/PPPc88993g4KtUUPUNQSrlEbm4uP/3pT5k0aRLHjh1j0aJF3H///SQnJ3s6NNUETQhKKZd45JFHqKys5M0332Tw4MG89tprBAUF8eijj3o6NNUETQhKKac7ceIE77zzDnfccQexsbEABAcH8x//8R98/vnnHDhwwMMRqsZoQlBKOd3f/vY3qqqqfnC94Je//CW+vr6sXLnSQ5GpC9GEoJRyuvfee49p06YxbNiw85aHhIQwb948Vq1aZbvQrNoPTQhKKadKSUlh//79XH/99Y2uv/HGG8nJyWHz5s1ujkw1x6GEICLzRCRNRNJF5KFG1k8XkV0iUiMi19otv0REku0eFSJylWXd6yJy1G5dovPellLKUz799FMArrnmmkbXz58/Hx8fH7788kt3hqUc0GxCEBFv4EXgcmAksERERjYolgXcBrxrv9AY809jTKIxJhGYBZwF1tkV+U/remOM9kVTqhNYt24diYmJ9O/fv9H1vXv3ZurUqXz11Vdujkw1x5EzhIlAujEmwxhTBbwHLLIvYIzJNMbsBS7UKHgt8KUx5myro1VKtWtnzpxhy5YtXHrppRcsN3fuXJKTk8nPz3dTZMoRjiSEgcBxu9fZlmUttRj4e4NlT4rIXhF5VkT8G9tIRO4UkSQRSSosLGzFbpVS7vLdd99RXV3NnDlzLljOun7Tpk3uCEs5yJGEII0sMy3ZiYj0B0YBa+0W/wYYDkwA+gC/bmxbY8wKY8x4Y8z40NBm54hWSnnQ119/jb+/P9OmTbtgucTERLp3767jG7UzjiSEbGCQ3esIILeF+7ke+MQYU21dYIzJM/Uqgdeob5pSSnVgmzdvZtKkSXTr1u2C5Xx9fZk4caImhHbGkYSwA4gTkWgR8aO+6Wd1C/ezhAbNRZazBkREgKuA/S2sUynVjlRUVLB7924mT5583vLKykpycnKoqqo6b/nUqVNJTk7m9OnT7gxTXUCzo50aY2pE5G7qm3u8gZXGmAMi8jiQZIxZLSITgE+AYGChiPyXMSYeQESiqD/D2Nig6ndEJJT6Jqlk4OdOek9KKQ/YvXs31dXV5yWE5ORk1qxZQ3V1Nf7+/ixatIgRI0YA9QmhtraW77//nlmzZnkqbGXHoeGvjTFrgDUNlj1q93wH9U1JjW2bSSMXoY0x+glQqhPZvn07AJMmTQLqb1D77LPPiIqKYvz48Wzbto0PP/yQW2+9lcjISKZMmQLAtm3bNCG0E3qnslLKKbZt28bgwYMZMGAAZ8+e5R//+AcDBgzgxhtvJD4+nptvvpnAwEA+//xzampqCAoKYsiQIezevdvToSsLTQhKKafYtm2brblo48aNnDt3jiuvvBJfX18A/P39mT9/PkVFRezcuROAsWPHsmvXLo/FrM6nCUEp1Wb5+fkcO3aMyZMnc/bsWXbt2sXo0aMJCws7r1xsbCyDBw9m69at1NXVMWbMGDIyMigpKfFQ5MqeJgSlVJtZv+WPGzeOHTt2UFNTw8UXX9xo2YsvvpiysjLS0tIYO3YsgM6i1k5oQlBKtdmePXsAGDVqFLt27SI2NpZ+/fo1WjYuLo6ePXuyZ88exowZA6DXEdoJTQhKqTbbs2cPkZGRlJaWUl5eTkJCQpNlvby8GDVqFIcPH6Znz54MHDhQryO0E5oQlFJttmfPHhISEti7dy9+fn4/mBinoVGjRlFXV8ehQ4cYNWqUTqnZTmhCUEq1yblz52z/2A8ePMjIkSNtPYuaEh4eTq9evTh06BDx8fGkpKRQW1vrpohVUzQhKKXaZP/+/dTV1dG/f3+qqqpsdyJfiIgQFxfHkSNHGD58OJWVlRw5csQN0aoL0YSglGqTvXv3AhAQEICPjw/R0dEObTds2DCqqqro06cPgDYbtQOaEJRSbbJnzx569uxJeXk50dHRzTYXWUVHR+Pt7Y23tzegCaE90ISglGqTPXv2MHz4cMrLyxk6dKjD2/n6+jJw4EAKCgqIiorShNAOaEJQSrWaMYa9e/cSEVE/tmVcXFyLto+KiiIvL48RI0ZoQmgHNCEopVotPz+f0tJSAgMDCQkJITAwsEXbR0VFYYwhIiKCtLQ0ampqXBSpcoQmBKVUq6WlpQH1N5tFRka2ePuIiAi8vb0JCgqiqqqK9PR0Z4eoWkATglKq1Q4dOgRAYGAgUVFRLd7e19eXAQMG4Ofnd159yjM0ISilWi0tLQ1/f3969+7dqoQAMHDgv+fP0oTgWQ4lBBGZJyJpIpIuIg81sn66iOwSkRoRubbBuloRSbY8VtstjxaR7SJyWERWWeZrVkp1IGlpaYSFhdGvXz969uzZqjoGDhyIr68vffv21YTgYc0mBBHxBl4ELgdGAktEZGSDYlnAbcC7jVRxzhiTaHlcabf8GeBZY0wcUALc0Yr4lVIedOjQIXr37t2q6wdW1jOEiIgITQge5sgZwkQg3RiTYYypAt4DFtkXMMZkGmP2AnWO7FREBJgFfGhZ9AZwlcNRK6U8rqqqioyMDIKDgx2+O7kxQUFBdO/enZCQEE0IHuZIQhgIHLd7nW1Z5qgAEUkSkW0iYv2n3xcoNcZY+5i1tE6llIdlZGRQW1tLSEgIgwcPbnU9IkJERAQ9evQgLy+PU6dOOTFK1RKOJARpZJlpwT4GG2PGAzcCz4nIkJbUKSJ3WhJKUmFhYQt2q5RyJeu3+cjISHr16tWmugYMGEBAQAAAhw8fbnNsqnUcSQjZwCC71xFArqM7MMbkWn5mABuAMcBJIEhEfJqr0xizwhgz3hgzPjQ01NHdKqVczHoPQmJiYpvrGjhwIH379gW0p5EnOZIQdgBxll5BfsBiYHUz2wAgIsEi4m95HgJMBQ4aYwzwT8DaI2kp8FlLg1dKec7+/fvp0aNHi8Yvakp4eLht1FM9Q/CcZhOCpZ3/bmAtkAK8b4w5ICKPi8iVACIyQUSygeuAV0TEOijJCCBJRPZQnwCeNsYctKz7NXC/iKRTf03hVWe+MaWUax04cIC+ffsyaNCg5gs3o2fPngQFBREaGqpnCB7k03wRMMasAdY0WPao3fMd1Df7NNzuX8CoJurMoL4Hk1KqAzpy5AhDhgwhPDzcKfWFh4frvQgepncqK6VarLS0lNLSUqKjo/Hxceh7ZbPCwsLo3bs3aWlp1LcqK3fThKCUarGUlBQARo5seI9q64WHhxMcHExZWRnFxcVOq1c5ThOCUqrFkpKSABg3bpzT6gwLCyM4OBhA51f2EE0ISqkW2717NyLC5MmTnVZnSEgIISEhQP1Nb8r9NCEopVosLS2Nvn374sx7g7y8vBg2bBigCcFTNCEopVrs2LFjREREUD8smfMMGjSIXr166UQ5HqIJQSnVIlVVVRQUFBAbG+v0usPCwggKCtKE4CGaEJRSLZKcnEx1dTUXXXSR0+sODQ0lODhYLyp7iCYEpVSLfP/994BzexhZWRNCXl4elZWVTq9fXZgmBKVUi+zbtw+AsWPHOr3uXr160a9fP4wxZGZmOr1+dWGaEJRSLXL48GECAgLo37+/0+sWEWJiYgDtaeQJmhCUUg6rqakhKyuLQYMGOb2HkdXw4cMBTQieoAlBKeWwEydOcPLkSacMed2UoUOH4uPjQ2pqqsv2oRqnCUEp5bCjR49SVlbGqFGNDmLsFP369SM4ONg2AY9yH00ISimH7dq1C4DRo0e7bB+hoaH06dNHm4w8QBOCUsph+/fvB/7dzu8KvXv3JiQkhOzsbB0G2800ISilHFJTU2P71h4XF+ey/YgIgwYNorKykhMnTrhsP+qHNCEopRxSUFDAyZMnCQsLo2fPni7dl3VYDG02ci+HEoKIzBORNBFJF5GHGlk/XUR2iUiNiFxrtzxRRLaKyAER2SsiN9ite11EjopIsuWR6Jy3pJRyhby8PE6ePGkbkdSVRowYAfx7Ih7lHs0mBBHxBl4ELgdGAktEpOE0SVnAbcC7DZafBW41xsQD84DnRCTIbv1/GmMSLY/kVr4HpZQb5OTkUFRU5NRZ0ppi7cV08OBBl+9L/Zsjk6FOBNKNMRkAIvIesAiw/aaMMZmWdXX2GxpjDtk9zxWRAiAUKG1z5Eoptzp8+DAVFRVuOUMYMGAAvXr14vDhwy7fl/o3R5qMBgLH7V5nW5a1iIhMBPwA+2EMn7Q0JT0rIv4trVMp5R61tbW2G8XckRCCgoIIDg7m2LFjLt+X+jdHEkJj96e3qC+YiPQH3gJ+YoyxnkX8BhgOTAD6AL9uYts7RSRJRJIKCwtbslullJMUFhZSUFAAuCcheHt7Ex4eTk5Ojsv3pf7NkYSQDQyyex0B5Dq6AxHpDXwBPGKM2WZdbozJM/Uqgdeob5r6AWPMCmPMeGPMeGdO16eUclxeXh5FRUX4+fkRGRnpln1GRERQXFysw2C7kSMJYQcQJyLRIuIHLAZWO1K5pfwnwJvGmA8arOtv+SnAVcD+lgSulHKfvLw8iouLiY2Nxdvb2y37jImJwRjD0aNH3bI/5UBCMMbUAHcDa4EU4H1jzAEReVxErgQQkQkikg1cB7wiIgcsm18PTAdua6R76Tsisg/YB4QAv3fqO1NKOU1+fj6lpaVuaS6ysu7rwIEDzZRUzuJILyOMMWuANQ2WPWr3fAf1TUkNt3sbeLuJOme1KFKllEfU1dWRk5NDYWGhWxOCdYrO/fv3c80117htv12Z3qmslLqg4uJiCgsLqampcemw1w2NHDkSb29v7XrqRpoQlFIXZL2gDO7pYWTVq1cvgoODdSpNN9KEoJS6oLy8PEpKSgD3JgQRISwsTLueupEmBKXUBeXl5XHmzBn69OlD37593brviIgIHfHUjTQhKKWaZIzxSA8jq+joaM6dO4felOoemhCUaoPa2lqysrJISUkhLy+v003oUlpaSkVFBXl5eR5JCNZ9Jifr2Jfu4FC3U6XU+YwxJCcn8+2333L69Gnb8pCQEObNm8eQIUM8GJ3z5OXlUVlZSWFhoVt7GFnFx8cD9fciXHrppW7ff1ejCUGpFqqrq+Pzzz8nOTmZwYMHM3/+fIKCgsjPz2fLli28/fbbzJkzh6lTp3o61Daz3qEM7r2gbJWQkADAoUOHmimpnEETglItYIzhiy++IDk5mWnTpnHJJZdQP/oK9OzZk8jISNavX88333yDiHDxxRd7OOK2yc/Pp7q6GvBMQujXrx/dunXT4SvcRBOCUi2QlJTErl27+NGPfsSsWbOoqanhlVde4S9/+YtteOghQ4YwZcoUqqurCQ0Nden8w65kjLH1MPLy8rJNa+lu/fr1Izs72yP77mo0ISjloMLCQtatW0dsbCyzZs2iqKiIq6++mk2bNjFlyhSeeuopvLy8WLduHW+//TYDBw6kqqqK5cuX06NHD0+H32KnT5/mzJkznDx5kujoaPz9PTNlycCBAzl8+DDGGNvZmHINTQhKOcAYwz/+8Q98fX1ZtGgRpaWlTJ8+nSNHjvDmm29y88032/5Z/frXv+aLL77g1ltv5aWXXmLAgAHcd999Hn4HLZeXl2f7OXz4cI/FER0dzfbt2zl16hS9e/f2WBxdgXY7VcoBBw4cICsri9mzZxMQEMCPf/xj0tPTWbNmDbfccssPvrleccUVbN68GV9fXx599FG+//57D0Xeenl5edTV1ZGRkeHRhBAXF0dtbS0pKSkei6Gr0ISgVDNqamr4+uuvCQ8PZ8yYMTz55JNs3LiRv/3tb8ya1fSgvSNGjGDt2rVUVFRw0003UVVV5cao2y4vLw8vLy+3zaPcFGvX0/37dcoUV9OEoFQzdu/eTXl5OZdeeinJyck88cQT3HLLLdxyyy3NbjthwgSeeOIJ0tPTefDBB90QrfPY9zDy5BmCdRjstLQ0j8XQVWhCUOoCamtr2bx5MxEREURFRXHPPffQt29fXnjhBYfr+H//7/8xduxY/u///o+9e/e6MFrnOXv2LGVlZZw6dQrwbEKIjo5GRMjIyPBYDF2FJgSlLmDPnj2Ul5czffp0Vq1axZYtW3jqqacICgpyuA4R4bnnnsPf35/bb7+9QwxvYR1htKioiD59+hASEuKxWPz9/QkODub48eMei6Gr0ISgVBOMMWzfvp2wsDAiIyN5+OGHGTNmDLfddluL65o6dSoLFy5k586dfPjhh84P1slyc3MBOH78OMOGDfN4d88BAwZw4sQJ6urqPBpHZ+dQQhCReSKSJiLpIvJQI+uni8guEakRkWsbrFsqIoctj6V2y8eJyD5LnS+Ipz9xSjWQlZVFQUEBEydO5J133iEzM5MnnniiVZPMe3l58eCDD9KvXz/uv/9+KioqXBCx8+Tm5hISEsKhQ4c82lxkFR0dTXFxMWVlZZ4OpVNrNiGIiDfwInA5MBJYIiIjGxTLAm4D3m2wbR9gOTAJmAgsF5Fgy+qXgDuBOMtjXqvfhVIu8P333xMQEMCIESN46qmnGDt2LPPnz291fWPGjGHRokVkZ2fz0ksvOTFS5zLGkJOTYxufqT0khNjYWE6dOqWT5biYI2cIE4F0Y0yGMaYKeA9YZF/AGJNpjNkLNDyfmwt8bYwpNsaUAF8D80SkP9DbGLPV1Deovglc1dY3o5SzlJeXk5KSwpgxY/j8889JT0/nkUceaVPTibe3N0uWLCE6OppnnnmGc+fOOTFi5ykvL+fMmTO2brLtISGMHFn/HfTAgQMejqRzcyQhDATsr+ZkW5Y5oqltB1qet6ZOpVxu9+7dGGOYMGECL7zwAkOGDGHRokXNb9iMcePGMXPmTE6cOMGKFSucEKnzWb+FW5tnPHkPgpU1Bu166lqOJITGvhI52k2iqW0drlNE7hSRJBFJ0lmTlDsYY9i7dy9RUVFkZmayZcsWfvnLX+Ll1fY+GL179+byyy8nJiam3Z4l5Obm4uXlRX5+Pj4+PsTExHg6JNv8EkeOHPFwJJ2bI5/wbGCQ3esIINfB+pvaNtvyvNk6jTErjDHjjTHjQ0NDHdytUq13/PhxiouLSUhI4M9//jPdu3fnJz/5idPqnzhxItOmTSMvL4+33nrLafU6S25uLmFhYRw6dIjY2Fh8fX09HRLh4eH4+vqSlZXl6VA6NUcSwg4gTkSiRcQPWAysdrD+tcBlIhJsuZh8GbDWGJMHnBKRyZbeRbcCn7UifqWcbs+ePfj6+hIeHs7f//53brnllhbdd9CcwYMHM27cOCIjI3n22WfbVVdKYwy5ubkMGDCAgwcPtovrB1DfS6t///6cOHHCdve0cr5mE4Ixpga4m/p/7inA+8aYAyLyuIhcCSAiE0QkG7gOeEVEDli2LQaeoD6p7AAetywDuAv4G5AOHAG+dOo7U6oVqqurOXDgACNHjuSjjz6ioqKCn//8507dh4gwZswYxo0bR2pqKl999ZVT62+LoqIiKisrCQkJIT093TZsRHswePBgSkpKbDO4KedzqFHUGLPGGDPUGDPEGPOkZdmjxpjVluc7jDERxpgexpi+xph4u21XGmNiLY/X7JYnGWMustR5t+kIt2+qTi8tLY3KykoSEhJ47bXXSExMJDEx0en7SUhI4KKLLiIkJIQ//vGPTq+/taw3pJ09e5ba2tp2lRBiY2MpLS2lqKjI06F0WnqnslJ2Dh48SK9evTh16hQ7d+5s1V3JjujVqxfDhw9n0qRJfPvttyQnJ7tkPy2VnZ2Nr6+vradRe0oIw4cPp6KiQqfTdCFNCEpZVFdXk56ezrBhw3jjjTfw9fXlpptuctn+EhMTueiii+jevTvPPfecy/bTEtnZ2URERJCSkoKPj0+7mv7TGot1qlLlfJoQlLI4cuQI1dXVxMbG8vbbb7NgwQKXDuo2dOhQ+vbty9SpU3nvvfc83hRSVVVFfn4+ERER7N+/n2HDhuHn5+fRmOxZu7+mp6d7OJLOSxOCUhapqakEBARw6NAhCgoKXNZcZOXt7c3o0aMZOnQolZWVvP766y7dX3NycnIwxjBo0CD279/frpqLoH48I0C7nrqQJgSlgLq6Og4dOsTQoUP56KOPCAoKYt481w+vlZCQQGhoKKNHj+bll1/2aBdU6/DSwcHBHD161DZTWXsRGBhI7969KSgo4OzZs54Op1PShKAUcOzYMc6dO0dMTAyffvopixYtcktzSXh4OGFhYUyYMIH09HTWr1/v8n025fjx44SGhtou2ra3MwSAQYMGUVJS4vHmtc5KE4JSYLuImpmZSVlZGddff73b9p2YmEj//v3p27evx0ZBNcbYLihb5y5ujwlBu566liYE1eUZY0hNTSU2NpZPPvmEoKAg5syZ47b9jxo1Cj8/P2bNmsXq1as9MsTzyZMnqaiosF0/CAgIaBdjGDUUFxdHaWkpOq6Za2hCUF1ebm4up06dcntzkVWPHj2Ii4tjyJAh1NbW8te//tVt+7ayXj8YPHgwycnJXHTRRa2aCMjVYmNjqa2t1Z5GLqIJQXV5KSkpiAjHjx93e3ORVUJCAgEBAUybNo2//vWv1NTUuHX/x48fp1u3bgQHB7Nr1y7Gjh3r1v07ytrT6PDhwx6OpHPShKC6vNTUVKKioli9erXbm4us4uLi6NatG1OmTCE3N5cvvvjCrfvPzMwkMjKSrKwsSktLGTNmjFv37yhrQjh27Bg62o3zaUJQXVphYSFFRUUeay6y8vHxsd213L9/f15++WW37bu0tJTS0lKioqLYvXs3QLtNCJGRkYgIRUVFOr+yC2hCUF2adRiEvLw8jzUXWVkH0VuwYAFr164lMzPTLfu1djONjo5m9+7deHl5MWrUKLfsu6X8/Pzo37+/dj11EU0IqktLTU1l4MCBfPHFFx5rLrLq378/oaGhxMXFISJuu7icmZlJ9+7dCQ0NZdeuXYwYMYLu3bu7Zd+tERMTo11PXUQTguqyysrKyM3N9XhzkZWIkJCQwNmzZ7n00kt59dVXXT4ZjDGGzMxMoqKiEBF2797dbpuLrPReBNfRhKC6LGtzUX5+vsebi6xGjx6NiDB9+nROnDjBZ5+5diLBkpISysvLiYqKoqCggNzc3HafEGJiYigvLycvL8/ToXQ6mhBUl5WamkpISAhr1671eHORVa9evRgyZAj+/v4MHjzY5ReXrdcpOsIFZStrTyO9F8H5NCGoLuns2bMcO3aMmJgYPvvsM483F9lLSEjg9OnTXHPNNaxfv96lfe7T09Pp1asXISEhJCUlAe0/IVjvoM7KynL7/RqdnUMJQUTmiUiaiKSLyEONrPcXkVWW9dtFJMqy/CYRSbZ71IlIomXdBkud1nX9nPnGlLqQQ4cOYYyhsLCw3TQXWQ0fPhx/f3/i4+Px9vZmxYoVLtlPbW0tGRkZtovYW7duZcSIEQQFBblkf85iPUPQ+ZWdr9mEICLewIvA5cBIYImIjGxQ7A6gxBgTCzwLPANgjHnHGJNojEkEbgEyjTH2cwXeZF1vjClwwvtRyiGpqan07t2b9evXt5vmIisfHx/i4+M5ceIECxcu5LXXXqOystLp+zl+/DiVlZXExsZijGHbtm1MmTLF6ftxtvDwcAICArTrqQs4coYwEUg3xmQYY6qA94BFDcosAt6wPP8QmC0i0qDMEuDvbQlWKWeoqqriyJEjDBkypN01F1klJiZSXV3NZZddRlFRER999JHT93H48GG8vLyIiYnh8OHDFBUVMXnyZKfvx9lEhKioKO1p5AKOJISBwHG719mWZY2WMcbUAGVA3wZlbuCHCeE1S3PR7xpJIEq5RHp6OjU1Ne2yucgqIiKCPn364O/vT0xMDK+88orT95Gens7gwYPx9/dn27ZtAB3iDAHqryOUlZVpQnAyRxJCY/+oGw4icsEyIjIJOGuM2W+3/iZjzChgmuVxS6M7F7lTRJJEJEmHvFXOkJqaSrdu3di4cWO7ay6yst6TcPz4cW655RY2bdpESkqK0+ovKyujoKDANnH91q1b6d27NyNHNmwNbp9iYmIoKSnh5MmTng6lU3EkIWQDg+xeRwC5TZURER8gELC/2rOYBmcHxpgcy89TwLvUN039gDFmhTFmvDFmfGhoqAPhKtW02tpaDh06RExMDKtXr26XzUVWCQkJtp++vr5OPUtIS0sDOC8hTJw4ES+vjtHxMCYmhnPnztmG7VbO4chvfwcQJyLRIuJH/T/31Q3KrAaWWp5fC3xrLEMRiogXcB311x6wLPMRkRDLc19gAbAfpVwsMzOTyspKW3PRDTfc4OmQmhQYGEh0dDRZWVlcffXVvPHGG5w7d84pdR88eJDQ0FBCQ0MpKSlh7969TJ061Sl1u8PQoUMByM7O1vmVnajZhGC5JnA3sBZIAd43xhwQkcdF5EpLsVeBviKSDtwP2HdNnQ5kG2My7Jb5A2tFZC+QDOQA7p8VRHU5KSkp+Pr68t1339GnT5922VxkLyEhgdLSUq666ipKS0t5//3321znqVOnOHbsmK15aOPGjRhjmDVrVpvrdhdrQigqKtLZ05zIofNDY8waY8xQY8wQY8yTlmWPGmNWW55XGGOuM8bEGmMm2v/zN8ZsMMZMblDfGWMtmq4HAAAgAElEQVTMOGPMaGNMvDHmXmNMrTPfmFINGWNIS0sjMjKSzz//nKuvvhpfX19Ph3VBI0aMwM/PDz8/P4YNG+aUZiPrkB3WhLB+/Xq6devWIXoYWUVFReHj48PJkyc1IThRx2gwVMoJsrOzOX36NAUFBZw+fbpd9i5qyM/Pj8TERA4cOMDSpUvZunUre/bsaVOdBw4cICQkBOs1uW+//ZZp06a122spjfH19bVdWNaE4DyaEFSXkZKSgpeXF1u2bCE0NJRLLrnE0yE5ZPLkyRhjGDZsGN27d+fZZ59tdV3FxcUcO3aMUaNGISLk5eVx8OBBZs+e7cSI3WPYsGGUlpZqQnAiTQiqSzDGkJqaSv/+/fnyyy+55ppr8PHx8XRYDgkODmb48OEcPnyYpUuX8u6775KTk9Oqunbv3o2I2Cbj+fbbbwE61PUDq6FDh1JYWMiJEyc8HUqnoQlBdQknTpywNS+cPXu2XfcuasyUKVOoqKhgzpw51NbW8sILL7S4jrq6OpKTk4mNjaV3794ArFmzhpCQkHY/oF1jhg4dSlVVFbm5uU7rfdXVaUJQXUJKSoptALfw8HCmTZvm6ZBaZNCgQURFRZGRkcHVV1/Nyy+/THl5eYvqOHToEKdPn2bs2LEAVFdXs2bNGhYsWIC3t7crwnapYcOGAdrTyJk0IaguISUlhZCQENatW8d1113XIf8Bzpo1izNnzjB37lzKy8tbPMXmv/71LwIDA21dNjdt2kRpaSmLFjUcmqxjsL4P7WnkPJoQVKdn/YdhHd3z5ptv9nRIrTJo0CDi4uIoKipi+vTp/PGPf3T4pqysrCyOHz/OlClTbHcjf/bZZwQEBHDppZe6MmyXCQ8Pp2fPnnph2Yk0IahO7+DBg0D9DVjDhw9nwoQJHo6o9WbNmkVlZSVXXHEFeXl5vPTSS81uY4zhn//8J927d7ddK6iurmbVqlVcfvnl9OjRw9Vhu4SIMHToUMrKyjQhOIkmBNXppaSk4Ofnx9atW1m6dCkdeWDd8PBwJkyYwLlz55g2bRpPP/00p06duuA2hw4dIjMzkxkzZtjuNVi3bh0FBQUsXbr0gtu2d8OGDdMmIyfShKA6tZKSEvLz8zl8+DAi0mGbi+zNmjWLXr16MWnSJE6ePHnB+xIqKyv56quv6Nu3L+PGjbMtf+ONNwgJCeHyyy93R8guY+16WlxcTEVFhafD6fA0IahO7eDBg9TV1bF+/Xpmz55NRESEp0NqM39/f6666ipbUnj66afJysr6QTljDF9++SVlZWVceeWVtgvpOTk5fPrpp9x8880d6u7kxgwdOhRjDMXFxToUthNoQlCd2sGDBykvLycrK6vDN4/Yi4mJYdasWUyZMoXa2lp+9atfnbfeGMPGjRvZs2cP06dPZ/DgwbZ1zz//PLW1tdxzzz3uDtvprOMxFRYWUlCgs/C2Vce4VVOpVigqKiI3N5fk5GT69u3Ltdde6+mQnGrq1KlUVlayd+9ePv74Y1asWMGyZcsoLi5mw4YN7N+/n4SEBGbMmGHbpqSkhJdffpnrr7/eNll9RzZs2DC8vLwoLi7WO5adQBOC6rT27dtHeXk53333Hb/61a8ICAjwdEhOJSK26wnXXHMN9957L6mpqQQGBuLt7c3MmTOZPn36eRfRH3vsMc6cOcPDDz/swcidp1u3brbpNDUhtJ02GalOyRjDvn37yMjIoLa2lp/97GeeDsklRISJEyfy1VdfISKsWbOGmTNncs899zBjxozzksGWLVt48cUXufPOOxk1apQHo3au+Ph4CgoKOHHiBJZ5uVQraUJQnVJubi4FBQX861//Yu7cucTGxno6JJeKj49n1apVpKen81//9V/U1NSctz4jI4MbbriByMhInnnmGQ9F6RojR44kNzeX06dPt3g4D3U+TQiqU9q3bx8HDx6koKCAe++919PhuMXChQt544032Lx5M+PGjeP1119n//79/O1vf2PKlCmcO3eOTz75xDawXWcRHx9PbW0txcXF5OfnezqcDk0Tgup0amtr2bt3L9u3b2f06NHMmzfP0yG5zU033cTGjRvp0aMHP/nJTxg1ahTLli1j0KBBfPfdd4wePdrTITpdfHw8gK3ZSLWeQxeVRWQe8DzgDfzNGPN0g/X+wJvAOKAIuMEYkykiUdTPw5xmKbrNGPNzyzbjgNeBbsAa4F6jDYDKCQ4dOsTu3bvJycnhD3/4Q4e+M7k1pkyZwt69e0lKSiIzM5Po6GjGjx/faY+DtafRqVOnNCG0UbMJQUS8gReBS4FsYIeIrDbGHLQrdgdQYoyJFZHFwDOAdcD5I8aYxEaqfgm4E9hGfUKYB3zZ6neilMWuXbv417/+RVRUVIeYJtMVvLy8mDhxIhMnTvR0KC5n7WlUWlqqTUZt5EiT0UQg3RiTYYypAt4DGo6Xuwh4w/L8Q2C2XODriIj0B3obY7ZazgreBK5qcfRKNVBWVsbq1avJysrid7/7XYeZFU21TXx8PPn5+RQXF1NVVeXpcDosRxLCQOC43etsy7JGyxhjaoAyoK9lXbSI7BaRjSIyza58djN1AiAid4pIkogk6QBWqjlJSUl88803DB8+vFPdmawuLD4+npycHGpqavSO5TZwJCE09k2/YVt/U2XygMHGmDHA/cC7ItLbwTrrFxqzwhgz3hgzPjQ01IFwVVdVW1vLc889R3FxMf/7v//bISfBUa2TkJBATU0NhYWF2mzUBo4khGxgkN3rCCC3qTIi4gMEAsXGmEpjTBGAMWYncAQYailvP8pYY3Uq1SJffvklX331FXPmzGH+/PmeDke5UWJi/WXKoqIiTQht4EgD6w4gTkSigRxgMXBjgzKrgaXAVuBa4FtjjBGRUOoTQ62IxABxQIYxplhETonIZGA7cCvwZ+e8JdXeVFRUkJqaSlZWFiJCaGgo8fHx9OrVy2n7qKmp4YEHHsDHx4fXXnut0/aoUY0bMmQIPXr0oKysjNxc/W7ZWs0mBGNMjYjcDaylvtvpSmPMARF5HEgyxqwGXgXeEpF0oJj6pAEwHXhcRGqAWuDnxphiy7q7+He30y/RHkadzrZt23j22Wf54osvOHPmzHnrRITJkyfz4x//mKVLl9KvX7827evee+/l0KFDLF++vFMMca1axtvbm9GjR5Ofn8+JEyeoqanRDgWtIB2p6//48eNNUlKSp8NQzcjPz+cXv/gFn3zyCX379uW6665j1qxZREVFISLk5eWxc+dO/vGPf7Bz5078/Py44YYbuO+++xg7dmyL9/fqq6/y05/+lIkTJ7J582Z8fX1d8K5Ue/fLX/6SN998kwceeIBly5YxcGCj/VS6JBHZaYwZ32xBY0yHeYwbN86o9m39+vUmJCTE+Pv7m9///vfm1KlTFyyfmppq7r77btOzZ08DmJkzZ5rVq1eb2tpah/b3l7/8xYiIGTJkiNmwYYMz3oLqoFasWGEAc88995jvv//e0+G0K9S35jT7P1aHrlBO89ZbbzF37lzCwsLYvXs3v/3tb+nZs+cFtxk2bBh//vOfyc7O5n/+539IT0/nyiuvZOTIkbzyyiucPXu20e1ycnJYvHgxd999NwkJCdxxxx1cfPHFrnhbqoOwXlguKSnR6witpE1Gyinefvttbr31VsaPH8+yZcs4e/Ys1dXV9OzZk4iICOLj421NRhdSXV3NBx98wB//+Ed27dpFt27duOSSSxg7diz9+vWjvLyc7du38+WXX+Ll5cWdd95Jnz59WLBgARMmTHDTu1Xt0blz5+jZsyeLFi1izpw5/OIXv/B0SO2Go01GetVFtdmnn37KbbfdRnR0NJdddhmVlZUMGjQIX19fysrK2LdvHzt37iQkJITZs2czbNiwJhODr68vN954I0uWLGHTpk189NFHrFu3jq+++oq6ujoAYmNjuffee1m2bBlr1qyhW7du500gr7qmbt26MXz4cE6cOMHJkyepqqrq8HNGu5smBNUmO3fuZPHixYSHh/Of//mfLFiw4Ae9fKqrqzl48CCbN29m1apVREdHs3DhQoKDg5usV0SYMWOGbfrHmpoaSkpK6NGjB927dwdg7dq1lJeXc8011+Dlpa2fCsaOHcu6deu47LLLyMvLIzIy0tMhdSj6V6RaLTc3l3nz5uHr68uLL77Iz372s0a7fPr6+pKQkMBdd93F/PnzycnJ4aWXXmLHjh0Oz3Dl4+NDaGioLRkcOXKEbdu2MX78+PMmkFdd28SJEykoKKCsrIzs7OzmN1Dn0TME1SqVlZXMnTuXkpIS3n77bRYtajje4Q95eXkxYcIEhg4dyurVq1mzZg0pKSnNni00VFRUxEcffURoaCiXXXZZW96G6mQmT54MQGlpKcePH2+mtGpIzxBUixljuOuuu9i/fz8PPfQQixcvbn4jO4GBgdx8880sWLCgxWcLJSUlvPPOO4gIixcv1nsO1HkSEhLw9/enpKSErKwsnWO5hTQhqBZ79913efPNN5kyZQpPPPFEq+oQEcaNG8cvfvELBg8ezJo1a3jttdc4duxYk9scOXKElStXcu7cOW688Ub69OnT2regOik/Pz/Gjh1LZmYm586do6ioyNMhdSjaZKRa5MiRIzz44IP06NGDTz75pM1jBgUGBnLTTTeRnJzMt99+y+uvv27rptq/f398fHwoKipi3759pKenExISwvXXX4+OfKuaMmnSJF555RVqa2vJysoiJCTE0yF1GJoQlMOqq6u55557yM/P58MPPyQsLMwp9YoIY8aM4aKLLiIpKYnk5GTWrl17XpkePXpwySWXMGXKFG0mUhc0efJknnvuOcrKyjh+/HirhkPpqjQhKIe98cYbrF27lgULFnDNNdc4vX5fX1+mTJnClClTKC0t5eTJk9TV1REYGEhoaKh2LVUOmTRpEgCnT58mKyvLw9F0LJoQlEMyMzN58sknCQgI4K9//avL9xcUFERQUJDL96M6n8jISMLCwsjNzWXw4MGcPn262SFUVD39yqWaVVdXx/Lly8nMzOSZZ54hPDzc0yEp1SQR4Uc/+hH79+8H4OjRox6OqOPQhKCa9c033/DBBx8wZswY7rrrLk+Ho1SzLrnkErKzszl79iwZGRmeDqfD0ISgLujMmTM8/PDDVFVV8cYbb2g7vuoQZs6cCcCpU6fIyMjQ+xEcpH/d6oKef/55du7cyd13382oUaM8HY5SDhk5ciShoaFkZWVRXl6u9yM4yKGEICLzRCRNRNJF5KFG1vuLyCrL+u0iEmVZfqmI7BSRfZafs+y22WCpM9nyaNscisrpjh49yp/+9Cf69+/PU0895elwlHKYiDBz5kz27NmDMUabjRzUbEIQEW/gReByYCSwRERGNih2B1BijIkFngWesSw/CSw0xowClgJvNdjuJmNMouVR0Ib3oZzMGMN9991HUVERK1assA0qp1RHMXPmTLKzs6mrq9OE4CBHzhAmAunGmAxjTBXwHtBwJLNFwBuW5x8Cs0VEjDG7jTHWqYsOAAEi4u+MwJVrrV69mi+++IJ58+axYMECT4ejVItdeumlABQWFnL06FFqa2s9HFH750hCGAjYDxuYbVnWaBljTA1QBvRtUOYaYLcxptJu2WuW5qLfSVvHQFBOc+7cOR544AH8/PxYuXKlp8NRqlXi4uKIjY0lJSWFqqoqMjMzPR1Su+dIQmjsH3XDS/YXLCMi8dQ3I/3Mbv1NlqakaZbHLY3uXOROEUkSkaTCwkIHwlVt9cgjj3DkyBEee+wx+vfv7+lwlGq1K664gu3bt2OMITU11dPhtHuOJIRsYJDd6wig4QzWtjIi4gMEAsWW1xHAJ8Ctxpgj1g2MMTmWn6eAd6lvmvoBY8wKY8x4Y8x4HdDM9VJSUnjppZcYMWIEDz74oKfDUapN5s+fT0VFBZWVlaSlpWn302Y4khB2AHEiEi0ifsBiYHWDMqupv2gMcC3wrTHGiEgQ8AXwG2PMFmthEfERkRDLc19gAbC/bW9FtZUxhmXLluk9B6rTmD59Ot27dycjI4NTp06Rm9vwu6yy1+xfvOWawN3AWiAFeN8Yc0BEHheRKy3FXgX6ikg6cD9g7Zp6NxAL/K5B91J/YK2I7AWSgRzA9QPkqAtauXIlW7Zs4dZbb2XChAmeDkepNgsICGDu3Lls2rQJYwwpKSmeDqldk450CjV+/HiTlJTk6TA6pbKyMoYMGYKXlxdHjx6lR48eng5JKad47733WLJkCb/97W8JCQnh3nvvbfM8Hh2NiOw0xoxvrpy2CSgA7rjjDoqKinj++ec1GahOZeHChXTv3p3Dhw9TVlZ2wVn5ujpNCIrVq1fz8ccfM2/ePJYsWeLpcJRyqh49erBgwQL++c9/4uXlxd69ez0dUrulCaGLO3XqFD/96U8JDg7m9ddf93Q4SrnE4sWLKSwspKqqioMHD1JdXe3pkNolTQhd3E9+8hMKCwv57//+b6dNialUe3PFFVfQr18/tm/fTmVlpW2uBHU+TQhd2DvvvMNHH33EvHnzWLZsmafDUcpl/Pz8uP322/n222/x9fXl+++/13sSGqEJoYs6cuQId955JwMHDmTlypVdrteF6nqWLVtGbW0tx44dIz8/n+PHjze/URejCaELqqioYOHChdTV1fHCCy/o8BSqS4iJiWHevHl89tlneHl5sXXrVk+H1O5oQuhijDHccccdpKSksGzZMq666ipPh6SU2/zmN7/hxIkTFBQUkJqaSn5+vqdDalc0IXQxTz/9NO+++y6zZ8/miSee0OEpVJcyffp0pk+fzscff4y3tzcbN270dEjtiv436EI++OADHn74YeLj4/nLX/5CYGCgp0NSyu0effRRcnNzyc3NJTU1Va8l2NGE0EWsWbOGG2+8kUGDBvH0008zfPhwT4eklEfMnj2bK664grfffhtjDGvWrKGurs7TYbULmhC6gPXr1/PjH/+Y0NBQHnvsMa644gpPh6SURz333HNUVlaSlJREfn4+O3bs8HRI7YImhE7ugw8+YP78+QQHB/Pggw9y8803axdT1eXFxsby8MMP88UXX1BQUMA333xDQYFO664JoZMyxvDss89yww030L9/fx544AF+9rOf4efn5+nQlGoXfvvb3zJ16lTeeustysvL+eijj6iqqvJ0WB6lCaETKi8vZ/Hixdx///0MHz6c++67j7vuuktHMVXKjo+PD++++y4BAQG8++67ZGRk8OGHH3bp6wk+ng5AOdeGDRtYtmwZGRkZzJkzhxtvvJHFixfTrVs3T4emVLszePBg1qxZw8yZM/n444+pra0lICCARYsW4e3t7enw3E7PEDqJnJwcbr/9di655BJKS0tZunQp9957L7feeqsmA6UuYPz48Xz++ecUFhby1ltvsWbNGt577z0qKio8HZrb6RlCB3f06FGef/55Xn75ZWpqavjRj37E3Llzufrqqxk5cqSnw1OqQ7jkkkvYvHkzV155Ja+99hppaWlkZWWxZMkSoqKiPB2e+xhjmn0A84A0IB14qJH1/sAqy/rtQJTdut9YlqcBcx2ts7HHuHHjjDKmvLzcvP/+++byyy83ImK8vLzM2LFjzX333Wc+//xzc+bMGU+HqFSHVFpaam6//XYDmO7du5sZM2aYP//5zyY/P9/TobUJkGQc+V/fbAHwBo4AMYAfsAcY2aDML4CXLc8XA6ssz0dayvsD0ZZ6vB2ps7FHV00Ip06dMhs2bDBPPvmkufTSS42vr68BTGBgoJkxY4Z58MEHTXx8vCkqKrJts3z5cjNjxozzXhtjbMvs10dGRpoZM2b8YJvIyEjbdg3LWB/GGFsc1jKBgYFm+fLltuWBgYHG29vbeHt727aPjIw0kZGRJjAw0AQGBhp/f//z6p0xY4bx9/f/QQzLly833t7epv67jLHVZa3XWk9gYKCJjIw0/v7+58VnX05EbOv9/f1t9Vof1tfW92U95oAREVv8gG299WGt29vb24jIeQ/rNvb7se7L39/f9tN6jLy9vW3vxRqPNWbr8bM/FtZtrMfbWpe3t7cJDAw0InLeMZoxY4bx9vY2/v7+tjqtcVuXNTw2rnrY70dEbM+tv1frMbD/XVvfn/WYWH/Hy5cvtz23fibtP9/Wz5n182Utb4wxu3btMvPnz7fFMXDgQLNgwQLzpz/9yezbt8/U1dW16W+6OdY4nMXRhCCmmTHBRWQK8JgxZq7l9W+o/2v8b7syay1ltoqID5APhAIP2Ze1lrNsdsE6GzN+/HiTlJR0wXg7gtraWs6dO8fZs2dtj5KSEnJycsjLyyMvL4/c3FwOHz5MRkbGef2jQ0JCiIuLY8SIEVx88cUkJCQQHx9P9+7dzxvf3XqvgXWZiNT/wu1+Wtc3vC/BfpumytiXdcd9DU3FqfdU1NNj4biGn2/rsoZ/MwCZmZmsXLmSTz75hIMHD9p6IAUEBBASEkJYWBjh4eGEhIQQHBxMYGAgQUFBBAUF0a1bN3x9fS/48Pb2tj28vLxszyMjI8nNzT1vXe/evVt9oVtEdhpjxjdXzpFrCAMB+8E+soFJTZUxxtSISBnQ17J8W4NtB1qeN1en01x++eXs3r0bS3w/+Gn/AWjsdWuXN1ZnbW0tNTU1zcbcrVs3+vTpQ3h4OKNGjSI2NpaJEycSExPDwIEDiYyM1HsKlHKxqKgoHn/8cR5//HHKysrYtGkTmzZtYv/+/eTl5ZGfn8/hw4epqKhw+j0MAwYMOO/1tm3bmDTJZf8mAccSQmNfOxqeVjRVpqnljfVuavRURUTuBO6E+i5irTF48GBKS0sREdu3APvn9iN+Wpc3LNtwm+bKNLUvPz8//P398ff3JyAgAH9/f7p3706vXr3o168fYWFhhIWF0atXL9uje/fujX77mzlz5nmjNTZWxn5Zw1gd2aapMo6scyZH4uzK9Fg4ztG/k+XLl/PYY4/ZlgcGBrJw4UIWLlxoW2aM4fTp05SVlVFeXk5ZWRmlpaWUlJRQUVFBZWUl1dXVjT5qa2upra2lrq4OYwxbt27l+++//0FsiYmJJCQkuGWKW0cSQjYwyO51BJDbRJlsS5NRIFDczLbN1QmAMWYFsALqm4wciPcHXnnlldZs1u5t2LDB9tzaFGT/GrTJqKvQY+G4ljQZNUdEbF/cnKnh37O7OHIfwg4gTkSiRcSP+ovGqxuUWQ0stTy/FvjWciFjNbBYRPxFJBqIA753sE6llFJu1OwZguWawN3AWup7B600xhwQkcepv3K9GngVeEtE0qk/M1hs2faAiLwPHARqgF8aY2oBGqvT+W+v65gxY8Z5r5cvX37eGcTy5cvPK2e/PjIykqioKGbOnHneNpGRkdx2220XLAP1p9KJiYlkZmYSFRVFcnIy9913H8899xyJiYkkJydz+vRpAB555BE2bNhAZmYmAKWlpUD9tJ6TJ08+r95t27YRHh5+Xgy33XYbv//976mtrbW9n8zMTG677TY2bNjAtm3bmDx5MsnJyQQFBZGfn09AQIAtPvtyVVVV+Pn5ERAQQEVFBTU1NbZ6Aby9vamtrbW9940bNxIYGEhZWRkiYmuKLCsrY8aMGT9ovvPz86OmpuYHQyH07t2bsrKy8/YD9deX/P39qaysxN/fn/DwcACys7OJiIggPz+fyspKIiMjKS0tpaKigoCAgPM+A9bjmp2dbTveTz/9NOHh4WRnZ9OzZ0/Ky8uZPn267RhFRUWxefNmfHzq/x0EBARw+vRp6urqbNepGh4bV7Eeczj/W7K/vz+TJ09m27ZtVFZWMmPGDNvv2vo5sR6Thx56iA0bNjBz5kyg/kza+pm0/+xa/xas5ax/I+2Bp2JptpdRe9JZehkppZQ7OdrLSIeuUEopBWhCUEopZaEJQSmlFKAJQSmllIUmBKWUUkAH62UkIoXAsVZuHgKcdGI4zqJxtYzG1TIaV8t01rgijTGhzRXqUAmhLUQkyZFuV+6mcbWMxtUyGlfLdPW4tMlIKaUUoAlBKaWURVdKCCs8HUATNK6W0bhaRuNqmS4dV5e5hqCUUurCutIZglJKqQvoVAlBRK4TkQMiUici4xus+42IpItImojMbWL7aBHZLiKHRWSVZWhuZ8e4SkSSLY9MEUluolymiOyzlHP5iH4i8piI5NjFNr+JcvMsxzBdRB5yQ1z/IyKpIrJXRD4RkaAmyrnleDX3/i1Dva+yrN8uIlGuisVun4NE5J8ikmL5/N/bSJmZIlJm9/t91NVxWfZ7wd+L1HvBcrz2ishYN8Q0zO44JItIuYjc16CMW46XiKwUkQIR2W+3rI+IfG35P/S1iAQ3se1SS5nDIrK0sTIt5sjEyx3lAYwAhgEbgPF2y0cCewB/IBo4Ang3sv37wGLL85eBu1wc7x+BR5tYlwmEuPHYPQY82EwZb8uxiwH8LMd0pIvjugzwsTx/BnjGU8fLkfcP/AJ42fJ8MbDKDb+7/sBYy/NewKFG4poJ/MNdnydHfy/AfOBL6mdXnAxsd3N83tTPAR/pieMFTAfGAvvtlv0BeMjy/KHGPvNAHyDD8jPY8jy4rfF0qjMEY0yKMSatkVWLgPeMMZXGmKNAOjDRvoDUT5c0C/jQsugN4CpXxWrZ3/XA3121DxeYCKQbYzKMMVXAe9QfW5cxxqwzxlgnod5G/ex6nuLI+19E/WcH6j9Lsy2/a5cxxuQZY3ZZnp8CUvj33OXt3SLgTVNvGxAkIv3duP/ZwBFjTGtveG0TY8wm6ueQsWf/GWrq/9Bc4GtjTLExpgT4GpjX1ng6VUK4gIHAcbvX2fzwD6YvUGr3z6exMs40DThhjDncxHoDrBORnVI/r7Q73G05bV/ZxGmqI8fRlW6n/ttkY9xxvBx5/7Yyls9SGfWfLbewNFGNAbY3snqKiOwRkS9FJN5NITX3e/H0Z2oxTX8p88TxAggzxuRBfbIH+jVSxiXHzZE5lbEJ92EAAAMbSURBVNsVEfkGCG9k1W+NMZ81tVkjyxp2r3KkjEMcjHEJFz47mGqMyRWRfsDXIpJq+TbRaheKC3gJeIL69/wE9c1ZtzesopFt29xNzZHjJSK/pX7WvXeaqMbpx6uxUBtZ5rLPUUuJSE/gI+A+Y0x5g9W7qG8WOW25PvQp9VPaulpzvxdPHi8/4ErgN42s9tTxcpRLjluHSwjGmDmt2CwbGGT3OgLIbVDmJPWnqz6Wb3aNlXFKjCLiA1wNjLtAHbmWnwUi8gn1zRVt+gfn6LETkb8C/2hklSPH0elxWS6YLQBmG0sDaiN1OP14NcKR928tk235PQfywyYBpxMRX+qTwTvGmI8brrdPEMaYNSLyfyISYoxx6bg9DvxeXPKZctDlwC5jzImGKzx1vCxOiEh/Y0yepfmsoJEy2dRf57CKoP7aaZt0lSaj1cBiSw+QaOoz/ff2BSz/aP4JXGtZtBRo6oyjreYAqcaY7MZWikgPEellfU79hdX9jZV1lgbttj9uYn87gDip743lR/3p9moXxzUP+DVwpTHmbBNl3HW8HHn/q6n/7ED9Z+nbppKYs1iuUbwKpBhj/tREmXDrtQwRmUj9336Ri+Ny5PeyGrjV0ttoMlBmbS5xgybP0j1xvOzYf4aa+j+0lv/fvh2jNBBEARj+t1JsRDtT5gxWYh0wt1CbFLmBnXcQLCwES/tACu9gEizEtRPEIwQbi/eEYQlY6E7j/8GQZLKBty/DPJiZhVHTNHu5vDvKvt/pexe9ZiMmsjdgDXwA8+K7C+KEyDNwUvTPgEG+HxKFogXuga2e4rwFJp2+ATAr4lhkeyKWTvrO3R2wApY5IA+6ceXnMXGK5bVSXC2xVvqY7bobV818bbp/4JIoWADbOXbaHEvDCjk6JpYLlkWexsDke5wB08zNgticP6oQ18b/pRNXA1xlPlcUpwN7jm2HmOB3i77q+SIK0jvwmXPXObHn9AC85Ot+XnsI3BS/Pctx1gKnfxGPTypLkoD/s2QkSfqBBUGSBFgQJEnJgiBJAiwIkqRkQZAkARYESVKyIEiSAPgC3XROxMSHaZwAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#Handcraft some GMM parameter \n", "mus = [-4.0, 0.0, 4.0, 5]\n", "sigmas = [1.0, 1.4, 1.2, 1]\n", "pis = [0.1, 0.4, 0.2, 0.3]\n", "t = np.linspace(-10,10,1000)\n", "true_GMM_pdf = sps.norm.pdf(t[:,np.newaxis], mus, sigmas).dot(pis)\n", "\n", "# Generate N datapoints from this GMM\n", "N = 1000\n", "Nc = sps.multinomial.rvs(N, pis) # Draw observation counts for each component from multinomial distribution\n", "x = sps.norm.rvs(np.repeat(mus, Nc), np.repeat(sigmas, Nc))\n", "np.random.shuffle(x)\n", "\n", "#Choose some initial parameters\n", "C = 6 # number of GMM components \n", "mus = x[:C] # we choose few first observations as the initial means\n", "sigmas = np.repeat(np.std(x), C) # sigma for all components is set to std of the the training data\n", "pis = np.ones(C)/C\n", "\n", "# Run several iteration of 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, sigmas) + np.log(pis)\n", " gammas = np.exp(log_p_xz - logsumexp(log_p_xz, axis=1, keepdims=True))\n", "\n", " #M-step\n", " Nc = gammas.sum(axis=0)\n", " mus = x.dot(gammas) / Nc\n", " sigmas = np.sqrt((x**2).dot(gammas) / Nc - mus**2) # we use std, not variance!\n", " pis = Nc / Nc.sum()\n", " \n", "# Plot the true GMM distribution and ML estimated GMM distribution\n", "EM_GMM_pdf = sps.norm.pdf(t[:,np.newaxis], mus, sigmas).dot(pis)\n", "plt.plot(t, true_GMM_pdf, color='gray');\n", "plt.plot(t, EM_GMM_pdf, 'k');\n", "plt.plot(x, np.zeros_like(x), '+k');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Gibbs Sampling" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# parameters of NormalGamma prior over means and precisiona\n", "m0, kappa0, a0, b0=[0.0, 0.05, 0.05, 0.05] \n", "alpha0=np.ones(C)\n", "\n", "# Lets initialize latent variable z using the responsibilities gammas from EM algorithm\n", "# Random initialization would be fine, but we would need some burn-in Gibbs Sampling itterations\n", "# For convenience, we use one-hot encoding for z\n", "z = np.zeros((len(x), C))\n", "z[np.arange(len(x)), np.argmax(gammas, axis=1)] = 1\n", "\n", "def NormalGamma_rvs(m, kappa, a, b, N):\n", " # Sample from NormalGamma distribution\n", " lmbd = sps.gamma.rvs(a, scale=1.0/b, size=N)\n", " mu = sps.norm.rvs(m, 1.0/np.sqrt(lmbd*kappa), N)\n", " return mu, lmbd\n", "\n", "def GMMParamsPosteriorParams(alpha0, m0, kappa0, a0, b0, x, z):\n", " # Function estimates parameters of posterior distributions of GMM parameters given:\n", " # - priors distributions (Dirichlet, NormalGamma)\n", " # - observations 'x'\n", " # - assignments of frames to gaussians 'z'\n", " # Input:\n", " # alpha0 (C-dim. vector) - parameters of Dirichlet prior for GMM weights\n", " # m0, kappa0, a0, b0 (4 scalars) - parameters of NormalGamma prior for Gaussian component params\n", " # x (N-dim. vector) - vector of observations\n", " # z (NxC matrix) - one-hot encodings assigning observations to Gaussian components\n", " # Output\n", " # alphaN (C-dim. vector) - parameters of Dirichlet posterior for GMM weights\n", " # mN, kappaN, aN, bN (4 C-dim. vectors) - parameters of NormalGamma posteriors for all GMM components\n", " #\n", " # We use an alternative formulas for calculating the NormalGamma posteriors!!!\n", " # (see function NormalGammaPosteriorParams2 from bayesian_inference_for_gaussian.ipynb)\n", " N = z.sum(axis=0)\n", " f = x.dot(z)\n", " s = (x**2).dot(z)\n", " alphaN = alpha0 + N\n", " kappaN = kappa0 + N\n", " mN = (kappa0*m0 + f) / kappaN;\n", " aN = a0 + 0.5*N;\n", " bN = b0 + 0.5 * (s + kappa0 * m0**2 - kappaN * mN**2);\n", " return alphaN, mN, kappaN, aN, bN\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To simulate burn-in, one can re-run the following cell to \"reset\" the expected_pdf" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd83Vd9+P/XuXvpSrq6utp72Zbjbcd2nDg7gZCQhCQkBRoKv6bQHxRKB1D4kZaWx5eWDr6U0UAZhUKzaMjAiZM4duI43rasYe2959Xd+57fH1cJiuLEki1Lxj7Px+M+fO9nnGFb963PmUJKiaIoiqJolrsAiqIoysVBBQRFURQFUAFBURRFmaECgqIoigKogKAoiqLMUAFBURRFAVRAUBRFUWaogKAoiqIAKiAoiqIoM3TLXYCFcDqdsrS0dLmLoSiK8nvl+PHjE1LK7LNd93sVEEpLSzl27NhyF0NRFOX3ihCidz7XqSYjRVEUBVABQVEURZmhAoKiKIoCqICgKIqizFABQVEURQFUQFAURVFmqICgKIqiACogKIqiKDN+ryamKYry+6ej4WV++8oR9PZ8yvNyufbaazGZTMtdLOUM5vWEIIS4VQjRKoToEEJ86QznvyCEOC2EqBdC7BFClMw696AQon3m9eCs4xuFEA0zaX5HCCEWp0qKolwsRl/9Di/95nlWO21UhOqZDAZ4/fXXl7tYyrs4a0AQQmiB7wHvA1YBDwghVs257CSwSUq5BngS+KeZex3Aw8CVwBbgYSFE5sw9PwAeAqpmXreed20URbloyKE6XmocIrtqE9c98Blu2b6Z6b4m9HodY2Njy1085Qzm84SwBeiQUnZJKaPAo8AHZ18gpdwrpQzOfDwEFM68vwV4SUo5JaV0Ay8Btwoh8gC7lPKglFICPwfuXIT6KIpykRg68ksCyRwqalag0WhgzYe5q9hLc2cHDQ0Ny1085QzmExAKgP5Znwdmjr2bTwLPn+Xegpn3Z01TCPGQEOKYEOLY+Pj4PIqrKMqym+zkwKSJHFcek3oX//hCCy+cHiW/eBUh3zhSJonH48tdSmWO+QSEM7XtyzNeKMRHgU3At85y77zTlFL+UEq5SUq5KTv7rKu3KopyEQg3PkVrSM+hZCt1A4e5wbKXxuNtvKK7ls3pU/hDQTo6Opa7mMoc8wkIA0DRrM+FwNDci4QQNwJfAe6QUkbOcu8Av2tWetc0FUX5PSQlz4+2MGTpJeBPZyj4It7YFt63po3/ebGbNdkWugYHGRwcXO6SKnPMJyAcBaqEEGVCCANwP/DM7AuEEOuBR0gFg9m9RbuBm4UQmTOdyTcDu6WUw4BPCLF1ZnTRHwJPL0J9FEVZZoHhOp4LutkodrA+38hdfI7DthdZufYTXLuqmUd7N6ALTyNlklQXonKxOGtAkFLGgc+Q+nJvBh6XUjYJIb4uhLhj5rJvATbgCSFEnRDimZl7p4C/JxVUjgJfnzkG8GngP4EOoJPf9TsoivJ77PGT36csuJ5ev59a3XaOOAbpn87hv377FW5dW8yJgJ4K8zTRSISJiYnlLq4yy7wmpkkpdwG75hz72qz3N77HvT8BfnKG48eA1fMuqaIoF71IIsJQYIoc7QZ6ifLlZIyrtHFirjX8bDLBlvrjuJwVWCIFtI+Nkt3Tg+obvHiopSsURVk0e7tfoCSZRyDoo6V0LeVjU6QNdVHc20mexcmfGSu4Kb+eg1MuPOEQPp9vuYuszKICgqIoi6a5dy/hSAU9WgMTYTMPFU8QL53GbJwgOqLFrPFxJJpgWAisgQmkTC53kZVZVEBQFGVReKNebP4ppoM6TmeXUTYd4DlDCy+1vUH7iTeoOvZrwr0TvCpqMaW5yYwkCIaCBAKB5S66MkMFBEVRFsXBoYOsFBkQDTFtsPN+McDLra/whYab+CN/LtGVO/DZirD2jmO0ncaTLGTC7WZgYODsiStLQgUERVEWRetkM+a4lRadBXs0ySvyNf586AEcESdGzU3cctpFus+POz3MpNdEjzQzGQgyMjKy3EVXZqjlrxVFOW9SSghM0B7NojutgJqIl9UxDXX6AfYVDGC1jvGBuIU72q/mZ5skDmFGL4NoQ0GSSdWPcLFQTwiKopy3Nncb1dEogyELAbOJ7ZZfcHq6kKsLD1JQ0Ysj40Ze1htYk9FMIq7Hm24jHBvFPO0lEYstd/GVGSogKIpy3o6NHmNd0kBvBFaaT3Hs9K1cae5mcNRJZWIA6X6NY13r+I1WgytgJmmKE9LEEHET034/wWDw7JkoF5wKCIqinLfJ0CQyoWXUrMERihKMWLFnd2HIGmRX3Q5Gm/ew1bCHWFxQ6hZ4wj58WWlMatKY9nlVP8JFQvUhKIpy/mSSwUiCpD1Kd/tGbin8JYgIj/yXhTf2fXNmzaJ21u/4IFfdbiCsyyaUHmC0L4oxEGR0dJTy8vLlrsVlTz0hKIpyXob9w+QlBW0JKxGNGYMvTlF+C4ceixLYu5d9UjKl1fI/QP/rTxNqfolI3EHQlk4s5EYEwoRCoeWuhoIKCIqinKeTYydZn4AerQefL59sbQdRqefAi23sA64BMhMJ7ie1ZeKTTx7FFk5gcYcZyzBj8PhArXp6UVABQVGU89Lp6aQ04KUrmYEYN7CufC+/2JXNN4MB0oGpvALabtyA32FhHfA3fg/avY8TjGnxuCwQkvgDATX89CKgAoKiKOdFSslE3IsnZicZiZObMUDxr/awHfDodPzo7rV8e6eZn30oF4DPAuMv72fCkEkgQ08kGGY6GMTr9S5rPRQVEBRFOQ+xRAyD1sCg7GEyXEAaEwz2hXkgFAbg0LZixovGsI2up2ZwJX1r0jADn5yYJNbVjtSbGIkZ8Hi9am+Ei8C8AoIQ4lYhRKsQokMI8aUznL9GCHFCCBEXQtwz6/h1MxvmvPkKCyHunDn3MyFE96xz6xavWoqiLIWO6Q4qDZl0YcEnHGzJfJ2WR3xsATxaLYPXb+XO7rvZ5txKe2kOfVXFADwIhA/uRe+x4TZqCfvCKiBcBM4aEIQQWuB7wPuAVcADQohVcy7rAz4O/Gr2QSnlXinlOinlOuB6IAi8OOuSv3rzvJSy7tyroSjKcmiZaqE6FqMhVoz0x1mddoJVLalddBuqCqlJwHD4OGPuV6ixbqa/Jo/JTBO5wKZjR4lOh5nIsqL3B9TktIvAfJ4QtgAdUsouKWUUeBT44OwLpJQ9Usp64L16he4BnpdSqn91RblEDPoHsXrqaAzWYAh6GWzz8/54HIDkVdcRGcnBnrOBrcNWmrVhnNM7GNpWBMCH/H6mRzsYz8nANDb2XtkoS2Q+AaEA6J/1eWDm2ELdD/zPnGPfEELUCyH+TQhhPIc0FUVZRhKJO9rNZDgPV7yH6Zfd5ANDegOTpeOMD9YTGxhmX9k68rrr0MtM3BtTP+rvB7wd9UTNZpKBIMGg2hdhuc0nIIgzHFvQoGEhRB5wBbB71uEvAyuAzYAD+OK73PuQEOKYEOLY+Pj4QrJVFOUCSiQTaISGQDJOWGNgpbGJorppANpLy7C0TSG1hYxvvJPhDcdprhphMtFNTLOJ8TQzWcCahlOIsAlfQofb51MT1JbZfALCAFA063MhMLTAfO4DnpJSvrWsoZRyWKZEgJ+Sapp6BynlD6WUm6SUm9Rm3Ipy8ej19VJizaU+nk4sEqXM2MLmQASA8e1ORLdErN9BX96/8lkTXJUf5GB6B4aJCsaucAJw7UAfUXeEAWMmvkBAdSwvs/kEhKNAlRCiTAhhINX088wC83mAOc1FM08NCCEEcCfQuMA0FUVZRq1TrZQmp9gf2op5eoxwQy+rgIAQCMcUbSvW06xpof3Xn+Tz//whjjR+isqEYEwzQbA2E4CdUhJqb2Q0005gOsDk5OTyVuoyd9bF7aSUcSHEZ0g192iBn0gpm4QQXweOSSmfEUJsBp4CMoHbhRB/J6WsBRBClJJ6wnh1TtK/FEJkk2qSqgM+tUh1UhRlCfR4e6iOtjMUuAVHpBmxzw9AY3oGfs8aMDv50Y//lKGB1IS0J/fDysrt3H/353EWFpKgns1AvPMYnlUPYJzsxePxLGONlHmtdiql3AXsmnPsa7PeHyXVlHSme3s4Qye0lPL6hRRUUZSLSyKZYDIwRiScpNA8RWlfKiB0lWUifIJ/euHTDA3lku+coOb2OIefsNHcUcoTv/5brvqTb9HvSKN0yse6ztMM2tKxtg6TSCSWuVaXNzVTWVGUcyMlPXGIxULkm4dYF4wC4K8q4dmuP6Z/KI+srCnu/NqrfPqql7j5L57HaA7T2LmO5w9eQ7A01Wy0dXSYaDhJMqjWM1puKiAoirJgkUQEqwgwGskkoNGSduIYRcC0EEiZz9OnrgbgfX/8Iv961w3c84mP8sj7S9n6gdcA+MlL9+Ffm2pK2p5M4u/txSf0eHy+5aqSggoIiqKcg25PN3naKMPjpeiSUcxHUgMPm9LSeK7+k8Tieko29fKdmmbiIxE8z3eTvWEjV28+gSNvErc/k2e9dwOwAYh1tzBidOAJBAiHw8tYs8ubCgiKoixYx3QHiagX/7gJswiQ25/qP+hxOtjdsQOAT971NJnvvx/rphxsV+bh29fPvbfdyoadBwH4xYEHGDIZMAOlHcfpc5USmPYwNTW1XNW67KmAoCjKgvV7+/FEE4zoLBRGB6jxpH6rb5I7iSb05K3p4YvrByB7BQC6LDPaDCM1tnJWrz2NxRKkd6SYwbQ8ANYMdDKR6SIxMo7b7V62el3uVEBQFGXB4skY0bCPKYONgpFGapOSBPDkxJcBuPbWkxiM6SB+t9CBZUMOofoJCovSqFh/GoDjxlsAWO/z4pYC8+S4ekJYRiogKIqyYLrENJEpCGiNZJ1sQg+c1hlo99VgsHv5mxtOQMHGt90jNAJ9npVrSzZQu6EBgMcm7gfgSsAz2IvePU18ZnE8ZempgKAoyoKE4iFkbAJDtyRq0pLVNQJAnS4fgPyrGlhhMEDpVcTj8bcNJbWsyabEl0VuyTjOrCneCG8nBtQA9LUSklLNRVhGKiAoirIgXZ4udDJKeMyMnjC5o6mhoodj2wC4bkczmrhk38HjvP7667zyyiucOnUKAKHXoNfqyHSaKatpIYqRLm0mGqC8uwGv1qqGni4jFRAURVmQzulOhMbCsMmJ0z9MhT81Ie1g4k7M9ikeunaEnr4BVq9ezbXXXsuNN96IVqvl9OlUv4GpxsFaQx6FV3QAcJIrAFg13MlAZinTHg+xWOzMmSsXlAoIiqIsyKB3AH0oTFdmMYVd9VQmJTGgkQ/g3HQSR1iPqexKnE7nW/esXr2akZER/H4/hqI0qnXF5JWPYrP4OZq4EYAat5uuzGLC4xNMT08vU+0ubyogKIqyIJHYFIwHGLK5KGw8jgZowkIUC5Xbm5C9veRvvfsd923bto0jR44gNIIcmxOzTUNlZSd1bAdgbTLJUMhHYnRMBYRlogKCoigLEg6PYugPE9RpyO5PdSgfpxijIcD11a0UZNnBnPmO+8xmM0ajEZ/Ph7ksg3yTldySbk6xFoA1gHewC/3klAoIy0QFBEVR5i0pk2gSE+gn9AgRoWAitVz1Ca4gp7yOLKMFm832rvevX7+ekydPYqpIpzLpwFXZxyRO+knDBuT01aMNBPH7/UtUI2U2FRAURZm3kcAIJnSEo3ZMST8VMzOUT3ANztomKm0VkFX5rvdbLBai0ShSKyhJy8ecFyE7c4I6UvfUDrURPeOuvcpSmFdAEELcKoRoFUJ0CCG+dIbz1wghTggh4kKIe+acSwgh6mZez8w6XiaEOCyEaBdCPDazG5uiKBexHk8PpjCMZhSRPTFEZSxBAqjnNirXN7Aty/KOCWlzlZaW0tPTQ64jF73FQE1VE/VsAKByYowpk4NwJLIEtVHmOmtAEEJoge8B7wNWAQ8IIVbNuawP+DjwqzMkEZJSrpt53THr+D8C/yalrALcwCfPofyKoiyhkxNtWAMRWtOLqD59DD3QhhZ7hqC22Ist0A+OivdMo7y8nK6uLjKqXeTpdBSVd3KaqwCoCoRodZQyNT6+BLVR5prPE8IWoENK2SWljAKPAh+cfYGUskdKWQ/Ma3eLmX2UrweenDn0X6T2VVYU5SI26G3D5JX0WtIpbE/NK2ggk5LqI+SlFwASNO/9taLRaNBqtWjzzJQlbGQVDXKadQDUAp2RJJ6RUTUXYRnMJyAUAP2zPg9whi0x34NJCHFMCHFICPHml34WMC2lfHPRkoWmqSjKMtDExtFOaAhowTk8AUATpaSVNSGMq0E7v5bf4uJi+gb6yTZlgi3BiC2PJFAFTE0OER0dVyONlsF8AsKZenjkAvIollJuAv4A+LYQomIhaQohHpoJKMfG1WOkoiyrjOgUwpeBVobIcQcAaGQt+TVNlGvMkHvFvNIpKSmht7eXQkc+MVsaxWUDdJGBDigaOIVwT6uAsAzmExAGgKJZnwuBoflmIKUcmvmzC9gHrAcmgAwhhO5saUopfyil3CSl3JSdnT3fbBVFWWSReARTMIJGONAnQ5T4Uh2/Y9ZqynNjrJB9kL9hXmnpdDoSiQT5lcVYTEbKSls4TRkA5UOdaEMhFRCWwXwCwlGgamZUkAG4H3jmLPcAIITIFEIYZ947gauA01JKCewF3hyR9CDw9EILryjK0jk21YM1KpmyuciaHKU4GScCWMttOHILcQkPpOXMO72srCxiDkExBgoq22iamaBWNjlOMi7UXIRlcNaAMNPO/xlgN9AMPC6lbBJCfF0IcQeAEGKzEGIAuBd4RAjRNHP7SuCYEOIUqQDwTSnl6ZlzXwS+IIToINWn8OPFrJiiKIvr0GgTZqml0eCgoukEAK1oyanuxWxeg06rXVB6paWl9Az3USDMZBWMcnpmCYvKQIghYzoyOa8xKsoi0p39EpBS7gJ2zTn2tVnvj5Jq9pl73xvAGRsVZ5qQtiyksIqiLJ9h93Eqw0bajVZubGsGoJEsbMUtxLgBdAtbtjozM5Pp6WnSjXbCQSMeRzFMQa2E/xYWrpycvBDVUN6DmqmsKMq8pAX60YdsTBs0OPq9ALRpCzBkuHGEA5BTe07p5ufk4Ra5iApBEqgGBvweRgaGF6/wyryogKAoyrxkhqcwxpyQjJA3lQoIY+kuitLSWCV65z3CaDan04k1NwMM6eSW9tODHT2QNd5KcGJSzUVYYiogKIpyVt54AmMoioiko0+GqYqkmoc0ZWYyXPkUaN1gX/hUorKyMsZjblzCQH5JG80zAxpLx3phyq1GGi0xFRAURTmrY+5xjFLiETbs024KiRECstf6SNjXkGHWg1j4onR2ux1fwE+JxoCztJvW1O7KFE5MoAlHVEBYYiogKIpyVkdGT2FER7PORsHBVIdyMzrM+V4i+ho0moWNMJrLobUi9YIRc2qZtFJPEOJSBYQlpgKCoihnFRg7hAkbp3VmXI19ALRq7PiFwBpOQNZ7L2j3Xmw2G5npmbjjDqI5xQBUJuKMxnVqLsISUwFBUZSzSnM3kmEqYtyoIW8sNfqn35qBRm+kKNoNuWvOOe2ioiJ0aQa85KGvSK1gUwM0Sg1+X2Axiq/MkwoIiqK8p1hSYo94sMQzics4ZYFRAPx5NnLTnVRohsFZdc7p5+TkENTHEEkn9hVD+NDhBHz+UcbU0NMlpQKCoijvqSMYRp9MEpu0QCTBiqQbAPOqBGbnWrKsWjiPPgStVovUQr7GTG5ZG204AXCNduMZUwtaLiUVEBRlgYY9IXY3jeANXx5j5Ju8AXQk8Ab1aFs15BMjAJgqNXgsm0kzzmvBg/ek0WjIFQY0dj+d2tTQ09yhCRKeaTUXYQmpgKAoC/B6+wSPHx0g127iX19so3fy0m/jbh89iQYDnRjJOdwJQDNaIhYTRq0FkZZ33nnk5ORg0RqIaXUM2lId1MXT02iDQbxe73mnr8yPCgiKMk9j3jAHuyb43I1VrC3K4Cu3reTHr3eTSC5ke5DfP9qhVzHYMmnSaMnvbgOgw2AhIASuQCfkrT3vPAoKCtAadfgSGQRdqWWwy0J+YtGkGnq6hFRAUJR5+sWhXj59beVbn/VaDR/aUMhTJweXsVQXlpQSi7eF7IwyhsxaSjxdAIw6zFi1VorjPeCau8X6wtntdrDp8MpczGVWAGpIMhRCBYQlpAKCoszD4HSIdLMe25z28rVFGbSN+kheok8Jo9E4xoQHe8SFz2tiZSI1ByFSasRlLyfHLMFgWZS8jHYzMpmLdWVqa84qoCUWYGR0YlHSV85OBQRFmYffnBzk3o1FZzy3szqb19ovzdEwjf4QeiGJjZpwdzupJfXlLMuskH0dWbb57aE8H3qDHoe0Y68cYAgzJkAzMYx7TC2DvVTmFRCEELcKIVqFEB1CiC+d4fw1QogTQoi4EOKeWcfXCSEOCiGahBD1QogPzzr3MyFEtxCibua1bnGqpCiLK5GU+CNx0i36M57fXpHFoa6pJS7V0jjt9qAVCXxuHboWDS7ieAFNkYmAtQKzLWPR8nK5XDiTWqKGGB1aFwCOwTHcI6OLlofy3s4aEIQQWuB7wPuAVcADQoi5jYZ9wMeBX805HgT+UEpZC9wKfFsIMft/0F9JKdfNvOrOsQ6KckEd7p5ka3nWu54XQmAzavFH4ktYqqWRGD5JWKejPxIjv2VmDSOhIWw24op2n9OS1+8mLy8PMxoiOuixppawyB+fJOK+NIPtxWg+TwhbgA4pZZeUMgo8Cnxw9gVSyh4pZT2QnHO8TUrZPvN+CBgDshel5IqyRA51TbHtPQICwA0rc9jTfOn9JqsfPYDe7qAhmqBk/BQAPWYDIXS4gm3ntWTFXOnp6WhMOiLCxLhzZqRRYBpCARKJxKLlo7y7+QSEAqB/1ueBmWMLIoTYAhiAzlmHvzHTlPRvQgjjQtNUlAtNSkk8kcSge+8flRW5abSOLGwLyYtdIJ7AHO3AkVlMsyebVckGAMZdZlxk4NIGwfregXIhhBAYHGb80kWiNPUVU5kIEIpKNRdhicwnIJxpkfMFDakQQuQBvwD+SEr55lPEl4EVwGbAAXzxXe59SAhxTAhxbHz80uy4Uy5ejYNeVhekn/U6IQRajbik5iQ0B8LoCJIVLWRwqIBaugHwF9jJSt+0qB3Kb7LYrUQTLtJXpJrfqpB0+4Nq6OkSmU9AGABmD68oBIbmm4EQwg78FviqlPLQm8ellMMyJQL8lFTT1DtIKX8opdwkpdyUna1am5SldaBzgqurnPO6tjbXypGOS2cxtkZfEEEU42QW7q5sVpMa7ZMos+HP3onDbl30PHNycrDEXKRXDRAHSoDeqSB9w+qXwaUwn0VIjgJVQogyYBC4H/iD+SQuhDAATwE/l1I+MedcnpRyWAghgDuBxgWVXFGWgD8cJ8105tFFiaTkUNckTf2TeId7qNF4aTnRwIjDQgESy4oq4hYL+dUrcJWWL23BF8HYcB+5esn0iAZTRwIHCdyAtsiC1eTDkHv+E9LmysvLIz1mQJpC9Ig0KqUPw8AkQ32DsH3Rs1PmOOsTgpQyDnwG2A00A49LKZuEEF8XQtwBIITYLIQYAO4FHhFCNM3cfh9wDfDxMwwv/aUQogFoAJzAPyxqzRTlPEXiiXftO2gf9fH1Z5vw+vykezpYlwhzuE9wouY6NK4s2vwx2uua8Q2f5tHXf8y//+vnaBqqX+IanB/z6HECBgOnB3yUTZwAoEUj8JsMuMLtizrC6E2ZmZmYNZKAFrpNqaGn2aOjeMbU5LSlMK9lCqWUu4Bdc459bdb7o6Sakube99/Af79LmtcvqKSKssRO9k2zvvid4+xbRrw8UzfEl2+tZs/LL3GDJonn4EvYtFo8nSU4LDYKrr2FQx2nSXbo+dD2LZhvq+HJn3yLAzsr+H82/gkacXHPCU1IicnbSDA3jTf6MqjlZQB6rQaicROZkWHILF30fDUaDSabiYDQMZhRAKFOirwThMYvvRFcF6OL+3+loiyjYz1TbCpxvO3YhD/CE8cG+Mubazi++wDrOmHqP5/EXXUTWbkfwFBZxPNXXsHzr/6Aq1Y4iGzPpu9AL5E2Dx/55FfIrgvw3ZPfRcqLu/O5KxjBqhsh05xD62gRq0kNOR3NtpEfzyHDogNxpvEm58+YYSEsrEwVlQBQHplC+D0XJC/l7VRAUJR3EY4lMRt+t/GLlJL/2NfJn99UzeC+VsLtffTs+RGv3bqNdK2eIfkChuQEaxpfIVCwladah2jpdfPjKws4+koz0akoG7bcSMWYnWc6n1nGmp1dkz+ERobJpZThwd+NMAoUppOVuZUs24UbJZ6d6yKUcCFWpQaRVEk/k1FBMpk8y53K+VIBQVHO4Ez9B7saRrh+pQvv0R52H34R64Hf0FJZyHWjPk6e3keXv5zOsr1EnVo+6n2NDTETO2Nd3NTcynS2h5eebSazaiVpk5KByR5GAiPLVLuz6xifIClCOII1uLtd1JKaLZwss5LMLiAj/9y3zDybgoICtDEXuspUAKgCGqZi+HyX1jyPi5EKCIpyBs3DPlbl2d/6HIknqOt3UxUM8/IbL1HdVcdodhZXWyRtbiNhZybZqxvw9P8BfZnrsH3+++T596EzFHCl34jO7SE70MCTj52k5qbbWD+cw69a5q70chEZaiVqhLojJnKCo6STZBygMI0Mwyi6/MWboTyXy+XCHMpA5/QQQkMeMD4QoH9EdSxfaCogKMoZnOxzs25Wh/KjR/q5b1UOLz//PLb4JMZIDJfTy4G+u/lxww6eOrqTbz3xOdKGI3gbuvj0P+3iR+V/xPjAbwlWWdgczmRMSqr8bn7UNU0ymqBYn0/zZPMy1vLdmSfr8Jn07D+SYDV7AGjWCrxGPc5QN2SvuGB5m0wmbBoNRqbp0qUBkDY4TFfvpbvvxMXi/DdDVZRL0FQginOmnTwSTzDsCdN8/A2cOQasr3ZxXOTxs198lqP9b18S+zVWsmJDHZ948O8pT2jwV4aJRb+Lu6iUwsgqxkc72HnCTP01m9nYeYpd0V2szFq5HFV8V+PRGOmaAQb0Vhr7s7ieRwHosRkJR+0YRAJ0iz9LeTaLxYBHo6HP4qJkH06GAAAgAElEQVTW6yF/coTR7h5So9iVC0UFBEU5g9mDgJ4+OcQObYDxkI+4PsnhoZU8vPdP8UWtWE0+tuxoxiiHiI2u4UBrES0n1vHN7p/y5x97EodzkuG0Z4mFpxk3TlGW00csXIrem2DQ7aG4uoguTxfl6RfPxLUmfwibzotZZ6dvpJDVpOZPjGenkR8tuKAdym8yZ6Ux5dcxnF0I3nZKgqNERvvPfqNyXlSTkaLMMe6LvLVOTzIpaRvx0Hf8KFvvvYnWx4f48p4/wxe1sr6qje/+7V+yfkc7RsOTRLN3UFtdjjXjc0y5J/jWz+4i5p1ix8jn+XAkzrFAC71jaxmPnGTtiT7qildQM+Hguc7nlrnGb9c45UEjouTIGiaHXdTSC0CwKINCezWZztwLXobCkiJCyQw8lalF7iriUyQ8gQue7+VOBQRFmaOuf5r1xZlAai2j1ZP95GSl8fKv9vJ3h/6aSNzIziuG+fSdf87TL3r5v//wcZ5+/pfs3zvM8aYBAtPfAarxer/KP//0Y0Tz/hl39DY+Fc3kkPFF7Bo9/cFT7LC7ONzUillnZjp88SzeFhjoI2YIMdJ4FSIJq3ADkKy0Yc42kVG+4YKXobi4GBlzEl6TB0A1ATr8F7aZSlEBQVHe4fSQ960RRq+3jREc6mHd7TfwL9/dynQojXWlE9xz4z/wgyfa+c0rj5KIx/lATQFP1eTRkmWn02BhDwm+yPeJTN7NN37wMBbLc1jJ484uPU9rvKATiENPMpqRzZXGDezq3nWWUi0dw2gT07oYh/akUUY7FiT9gMi349QPIS7AkhVzORwOtCEnkYJUq3YVUD+UuOgn9P2+UwFBUeaIzex/MO6LkN/bgctl5//76x5OjdTisEzzJzcc5pEnn+ZkTztOaxrH1m7mmY5h7mwdpmbSS3k0yPXAN4EOWtnc+gDfP3A3/lg/RdkrKPC0E44JfP40KuMGXn/jFCPBi2NOQjCRxBZtx2cw09xkYPXMkhWtOg1TWgOGRABMZ18O/HxpNBqsyTQMdi8etGQC4f4phicuniepS5EKCIoyi5TyrRUZnj45gM49RPHWm/np46k1Gf/kumYeO/Y3NPYPUpxuob5QsvHUUUQiSWtpNftuvI29H3iQl+7/AsdW1mIH/h03Nz//1+wO1mCINXLHxHoOGluQGjPZupcZHp6izF5Ny1TL8lV8RksgRJZlkiBWuofzWM0BAPrsJkIhB1bj0o1DsZoMWPHQbUg9rTlHhmjrvTgC56VKBQRFmWXEGyY33YSUkmhbD2lWAx/+fIRw1MSOslP4tD9g36l6bHotBzM05LX6CRpNPPGx2zn1uSJabrFRX5lB3N+Ib931fPuDN+BD8FG8rP6P7/O6WEXIfpCbx3Lwmibwd25m04ownad07Onbs9zVp37Kg82QwDOZg8+XzhUiNcLInWMnN+IkIytnycqSnpuBSROlPz216mmhZ5D+lpNLlv/lSAUERZnlzf6D+gEPtpFOfFU3cGpPATpNjA++71f810upbT1eLkgjv9eP12zhy//v+/l1jY5n4nZeSR/k4MY9/PtNYV5zTVLt2sAf3/FRPMCHouOYfvkkx8RmSjyjjMaGScRMZKZNMN3aQCwAieTy7h082tuD3ixxt6c2H7iCPgCCxQ5WmK1kVW9dsrKUVJQR1+oZy0+NNCqNjBDo6Vqy/C9HKiAoyiwtIz5W5No5cKQbrTbO3/6zBYng9m27eX7/C/hCYb6ak86VPdNEdFr+8aHbidigMllFEgcZ3uvJGfsrtlirOViWoC42yIP2Ej6S81niwCf6W0mv28X+NBvrPE7izj58dVeyqaqH8cE1nBg7saz11w02EzRE6DxWi54oVdJPEpCVaRRkR7AWr1+yspSVlRGLZeKrTgWEKjnB5Lj6yrqQ1N+uoswSjqUWtTO21XMoZxud+5zoNDHWFe/ilYZ6Vuu0PDyRWmTtqTvuAn0O+d5KNL44VYF01tgaWZX+a3ISPu6zQ+6OKMGaF7nnjpX8Bal+iDteOMpk1Ecw0kAkECPkNpGVlSQen+TQyVPLVvdQIokt3smwjNN9IoNqGtADXUA0LwODVQMGy5KVx2q1oom4cK9IzQavJsSpsaXL/3I0r4AghLhVCNEqhOgQQnzpDOevEUKcEELEhRD3zDn3oBCifeb14KzjG4UQDTNpfmdmK01FWXaHOycIyzB7f1mERMOdW3bz2L5U+/6jNh26RJIjV6xkyLUCrTDRWNhHe/UQ0eopenNuIVbwKUprv0xioIzJESc94etxrvwFJ1c9zHPoSZeS+x59nFfNLkwRH2bLFP7+q1iv28PIYDaxWHxZ6t3sD+G0TDASNtE/mM1qXgKgzaBjVJpYjh9RYywLd3Fq7+ZKoGMouuRluJycNSAIIbTA94D3AauAB4QQczdT7QM+Dvxqzr0O4GHgSmAL8LAQInPm9A+Ah0gNMa4Cbj3nWijKIvCGY9iMOlr2NfBqRhUt+zLQaOLU2PZzeriNTxm11E5H8JvNHNp6I4OZQQ4VNZLlgrxMJ/2mtcQCnRyePs7zB5+iK5qkxVTCkYo8/sP6WbZ+rI5P8VXcwMphPzsbD3BcN4BWBvC128k2TTG2wsWLLx5elvqfnHSTYUrQ115NPK5nveE1AAYyzAQCWZizCpa8TFaDDaMpwLjQYwG0gyP4QiooXCjzeULYAnRIKbuklFHgUeCDsy+QUvZIKeuBuTtY3AK8JKWcklK6gZeAW4UQeYBdSnlQpmaa/By483wroyjno3nIS3WOjdHpASYPVJFI6rixtIX/PvVL0oD/M/O/e9/VO5m0ONhTG8eVv4F0bSGt4ZuJhMtIDJeSdVpH0JDLti2F3B3/X+4aeZKtU/sJJNIxXH8nXyT1xbrjuSZCcoh+7QDWRIx4ZBVbYidoGhwmHn33zuVk8sJMzpro7cGeYaOzvhaAtZwGwOtKpzJiIrPyyguS73tx5jux4aHXlBp6mjveT0tfx5KX43Ixn4BQAMxeVWpg5th8vNu9BTPvz5qmEOIhIcQxIcSx8fHxeWarKAt3ethL2B3imCOL40+nJl+tLnqC3vFBvqLXkhFL0J9XyPHqDTx2tYvVfjtbxnvod76f7U4NX3v23/lKz8/5atnzPGx5lRfdgo68W8jPNVKWrKNy7DD3FOzmx/yYOsAVifGJXb28YG0g1xBlsHsFpYFDNFRpaNx/5qWej/dO8fAzTfz0QPei118z2E7SFKPnVKoBYHU8NeY/WJHFanOSvKoLv2TFXKVV5dhlgEFHave0En8fvacOLXk5LhfzCQhnajic768o73bvvNOUUv5QSrlJSrkpOzt7ntkqysJNBaI0n2wmcjgLT8DOaucYuzp+hgv4s0Tq8eDlq7fw1I4V1PbWcaWpnnjhdrZ2HOQDu5/CtnInjr9+jMw/fIHgNSv5x+0VtIUrecr+Eba6ojhsQ6yOH6KmqpjPkRqtU32omzsHvRyxNuGYBBGyk2lNMjg5SiQYe1v5pJQ8Vz/M1z9YSySeZMwbXrS6BxNJ0hK9jCfCDDWVkskUhckIQSBU5SI/X4fOaF60/OaroqICA5Lx4nwAypODjNaPLXk5LhfzCQgDwOxF3wuBoXmm/273Dsy8P5c0FeWCkMkkB6xJWnenVvO8qvhntAz28XdaDeakpK2snP+87S42D7/GLem9pFtWMtkNFVNR+nRp9FoNHD20n0NPPom3dRWdnkN886YiYoMxflj7XWrWRNHljPKxnB/zGt/nKcCcTLL2hWH2W/ZTpDHS2ZvPDbKOo4XdNO1/+49E05CXjSWZCCG4Z2Mhz5xavB+ZJn8Il3WUo+0uwgELm6yvANAITDvT0Zgu/JLXZ6LX60lGnbgrSgCoZoKWwQu/dMblaj4B4ShQJYQoE0IYgPuB+e4Qvhu4WQiROdOZfDOwW0o5DPiEEFtnRhf9IfD0OZRfURZFNJ5kfCqItdNDR18pmeYIjb6fkwd8IimRwPc/chc5E/u5IdrLtNnF0cFi1vvGMOXms33VWq79whe56r6PsO3ej1C5bSe+nhr27X6Uq51xQqc6ObDqu2jKBWsrD1CU6+Thmc1eahp7eKjdR729AcNYHlneCRqlllg0QTT8uxFHr7SMccOK1Exhp83IZGDxOldPjoyTYY1z8kQlANtMqR/HVquB8ZCFhGv1ouW1UIZYNpMr3xx6GqFJDT29YM4aEKSUceAzpL7cm4HHpZRNQoivCyHuABBCbBZCDAD3Ao8IIZpm7p0C/p5UUDkKfH3mGMCngf8EOoBO4PlFrZmiLEDHmJ/ehJeJA6lBcNcX7uFwZzN/ARik5HhtLf1ZAzw01U13kZXW6R1UZzuwrF5P5cAw/WsreebZH/Gbp7/Db3/7CCfrX0PklaIrfD99bae5TvZw/HAHU7UP0J7v4hO5T9DAt3kMMEoJBySnjbtYocmmr1uwUrRjW6en+cDwW2WMxBOYDdq3PqeZdPjCMRbDdFc7GZl22utS/QcbYscBGM22Y/Sl41y5fDuVmTXZjBSnkwTKgaEh/7KV5VI3r5WqpJS7gF1zjn1t1vujvL0JaPZ1PwF+cobjx4Dl+7VDUWY53D+JK1bPb9/4YwCi+m9jTyb5tBAgJftudXG7183rOQU0hDZwe24Bzsw0ogd2c7hUQ4VsYuc12zCZckgkgvj9rXh9BygszGLMejOalj38kXeC703dyE3VZjJ6T2Fr+VP+LryDe3md7Sfr8V1zLf1pI3gaTVy3eYofT5zmPl8l8WgCXyxBhvnt+wFsKnFwom+andXn37emHWnBWKRnoKkagJXh1FiQQImTyriWitLS887jXOUVFqD3jzGsMVCQjJI21IcvME6aVfUpLjY1U1lRgL0DU8SOgD9soTZnmEN9B/gMYJGShpoKxnLDvBS4lUFnCX+als2Q1cjosRdY4dByz4N/z4YNf0J6+hqMxhwsljJcrluprPgrsrNvptR2hPRNdgyZGdzWepK96fej3ZbkpsIXaeYbPAoYZJL4CfDzW8pMVxIcnWYsMMyKbXm0HBzmcPck2yqy3lbmNYXp1Pef/3LQQ+EoWYZB+kfT8I85sRq9lEf9JIDYylzKMiUG3fJ9VZRVV5Ce9DJgSQ09LZjuob2nYdnKcylTAUG57CWkxJl4g/0v3QxAWdr/IegP8vmZmblv3Bind/iPKS9qxDIS5LhMkjFwhA8EddQ89PBbM3i7gxF+MTTBt7qH+Vb3ML8cmmSMbK5Y+VVGgzlY1zRTbZygqmOczpKNbFjxMnAV/0Rqw5mbX3+NloxaIroY3ldGWc9xBsxRvJNhTg96WTmzac+bTHotscTcqT8Ld3Taj8s6zt4jxQBcmfcqOqAV8JfkkVXkOu88zkdxcTEZ0stothOAsmgX7Udal7VMlyoVEJTL3gujbrIHO2nqWolJn6A3/L98HMiUkt6SLPbbb2dVcD8d2asotmQwmK7lJpFLxi23IXQ6TngCPNw+yBvTfq532PnL0lz+ojSXaxxpvDrl4/+OTHA4eh1FRQ9ivWaUjaNHGBnJQffhJJtcJznFN9gNWJNxtj57iJNZvVji11GcHOGHXc1UrHehGQyi1VyYpSNau7vIyrSxd39qJM827bMANBt19GCG0u0XJN/5MhgM2GNxRopT5atigMYj+mUt06VKBQTlsnfk9B4O7U51mm7If4GmgUE+M3PuuWuzsU/kMrS5kNyxEUwbC1gTsZE5No5+yxb+uXuEjlCEhyvz+Uh+FgUmA0IINEJQZDLwhwVOHq7Ix4Pkp/51WBxXsvbGENuPn+JY3lVsr3kOuI1/EakuuGveOEhSF2Q03UXG6VHC/jocxTbEVOyM20dmWAxMne9oo94OrBk2Wk7UpP4OAqmlMwYdVvCYKKlee37pLwJDLIvxyjeHnk5yetB+ljuUc6ECgnJZa/aHyHe/wvFjOwAwaL/NDcAKYNJuYn/6hyjU95Cw6yFjio74Om5oa0b/4ft5uGOIe3MzuS/XgeY9Fn7TCMFajYGP5jr4eeJDDLoy2V4xQfEbw+g/EcBlGecl+Q8cB5zREJufO8SBkij+rutYHXiBx3sncJSmMdzpeUfaq/LtNA97z7n+4UQSc7CbcbcL32gOOlOUGl+qQ9lXmEV6MJMVecs/7t+qyWN8ZWpyWjUxOseXZ17EpU4FBOWy9mLvIcaPFDA46SLN6OeU9+BbTwcvbK9g48Qwx3euI3+yg6or7qSktwOD1cY/RbV8sSyXEvP8vpgqXTZ8ngh/V1XAbuNDyG1WVgyO0ZFTzjXlzwMP8C/CBsCagye5t/tF9pesp3RgiqcHuli/NY+eUxPvSHdFbtp5BYR6X5D8tEHeOFAKQE71IOXBVEd1YGUB2RbbsnYovyk/t4rxwgxiQDEwPTRNLKb2V15sy/8vrSjLZDoWxzH+Anv3XQvASscPyJgI8AEgphU0Z21jqLqCMk0zQWsRTYFcbm+q55Ht1/PVinwy9PPfX3hlXuo3eaNGw9+vWM0vjbdRWdzLVa/1k/EHvWiFhsflF+gGSgLT6BvjaJNDDE3cRV7/EzizLNidJqZHg29LN8NiwBs697kIh3p6cWUmeXFfKQBXZ72KWUo6AHfNCsxlNeec9mIqra4kqdcwqEsFYOd4FxOTbctcqkuPCgjKZet/+5tI65M0tqbayKfSH+VPSf1Q7K+tJDPNQEdtOmX+OJp0CyX1DewvqeTPVpSRptO+Z9pzlTgs9E6mvsx1GsE31+7ktS3XkxsSBCrMbCp4gwSf5Tsi9SNZ8PoR7h3czetZFayYGObXLW3UbM2j9fAZNpk/j30KQu0tZGRlc+J4KQA7vL8F4KRBS58pg/UbNp5z2ospPz8fe9LPiC0NgJJQKw11ajvNxaYCgnJZSkgJE89x+JXVjHnScBj7mexu4pMz5/sqr6CldBVb5AFC+rVMmyqobW5g+30fIse48BEuGo142+qNJr2dj63MR1oGuK7ZQ9VNJwAnP9F9GA9wxVgvE8OF3DRwiK6MNXQd2o3bM4lWpyEcWJzZyZFkErOnneHxYvzudAxpYWpHUju29TptMAk7apZ3yOmbzGYzjqSHsazU0h3VtFK3f3n3n74UqYCgXJZeGhulaNrPiYbU00FW8Xf5UCRCJtCb68KXnsVwqYbcmJ4BMUTF6x3E77mPVWnnt47O7L0MyvM/iGtdOraEFetVo+Sl9eONfZUfzpzPeeMNVvqP4Y5UkZvewLP7XqV4rf1ty1kAWA1a/JGF77J23BOkwD7MnldS6wQ5a0YonUgtmOevyMHuTceVZjq3il4AmZEEI8WpslYxSENz2jKX6NKjAoJyWeoZ/g2aY/kcbU0t1RAO/JbPzpxrWlXLoSt2cjUvYDe9n0mZoNDn5Y6rt55XniUOC/3u3/UBCKFl47X3YBadbJw0s2X9fmAV37dsIg6s66yjx7OBB1te5rhpGyflOEfrDhLwhUjEfjchrcxppXs8sODyHOnuJccR5eU9ZQBUl7dRHAkQBXxrVpA0Os+rvostlwwGVs2UFTcdIyogLDYVEJTLzmlfkLxwH3vqS3EHzGSkn6JsqIk1gM9soje7nClXDIeI0RUZY93BKd73xb847z2FVxek0zj49hFB2a5r0Dohm2yKbzmFQRuhJ/gNHge0EkpOHCeuO8EVPRrGDDqcmVYmaaf1yO+eEsqzbXRNLHzBt3BHE6a0AhrqSwH4QGIXGqAOmKxYh7Mw99wrewHkO1cyvjI1X6OaOH1jBuJx3zKX6tKiAoJy2dnb9wpVpzM42ZIaQZOZ8chbQ03by0t4betN7BRPk2G6m97hBrJySnHknH9bekW2lY6xd35xV2+7myLzMSo12WysOAjcxHfMeQBc0fAqvb5rubP3FB5NJq927qV27QrqG+vemqhW7LDQNxl8R7rvxROLk+5v5sjxNcSietKKJ6npPAFAk1XPWELH9Zuqzq/Ci6yk8grcLitBBC4gMTyKz9u+3MW6pKiAoFxWpmNxHKHDuNsKOdhSCEis009yF5AQgtfW3shkpiTPMIyWImpPjfCRv/7zRclbp9WQSL5z7aGqDbcyhJb1Rgdrd+4HBEcTf8OrgCUWp6ipk4G0w3ziZAtPOzYx3rSXrAIbR/alFngz6DTEFrjP8r4JD8WZQ+zalWqTz7liiNz+1KidkdwMNBMx3rf54goIubm5CK2WQUNq6GlRsJH2033LXKpLiwoIymXlf/sbqZrQcrDThidiITPzFT7sGUcHdBe5OLj5Kq5KPkuMHZze/xjBq29Db1zcWbFzl6DQaLSkOzdRmPM6tVlQ7OgmGf0k39anOnTX1+3ijdAq1nqmIW6le/okFSvyOd3UjM93bk0mzc2nyc6ycvhoqk1+ZXkL1e7UkNaJ2jKMfgeWcxhNdSHZbDYypJeRtNTM6WrqOPz6wjvTlXc3r4AghLhVCNEqhOgQQnzpDOeNQojHZs4fFkKUzhz/iBCibtYrKYRYN3Nu30yab567OMa3/f/s3Xd8HMXd+PHPXNOd6ql36dQtq9qS5d4bpjeH3glpPIGEHwQSEhzCk5CQACFAAg8QejfFNHfLvcmWbFmyZKtavetUrt/N7w/JxBiDDRhLcvb9eumlu93Z1XdX0n13Z2ZnFGctt5SInk8I3xXKplITAMG+f+e24fVbZ8yl26gl03cfIdrZ+HS28JPrbzqtMUQaDbSYvzwXcuqkObQYIsgOCmXOxFWAgU99buEwENRvZWKFoMt7J3/YsZqHQy/DsvkpJmZPYtUna3G73SDlCcc6OhGPlGhq93HoSCyd7Ua0vlYutRTi73FTD3TmzcTj5XNaj/t0EEIQ7BygLfRo19PD7N87+uIcy06aEIQQauApYAkwHrhKCDH+uGK3AD1SymTgMeDPAFLK16SUuVLKXOA6oE5KWXLMdtccXS+lVGbOVnyv1rS1YnLZcA42saY6DXAxu201YUCzvw8f551PnmUDXSRh/uRTSufPwF+nO9luv5HMqAAONH15TKLwpBScgwnExuxkQmIteq0FR+/veHz4QbW8ko940ldPjB0CBgUfB0TRZ95GgIxn8+bNBPro6LWc2vMJJf0WYn2qWLOmAIDQzBaiK4aGk96uEXSEpWEMDf66XYyYKLcXTQnRAKTQTGWj0tPodDqVO4QCoEpKWSOldABvAhcdV+Yi4KXh1+8C88WXu2RcBbzxXYJVKL6Lupb3Gbe1mQ/2TcHm8iIw+BV+5LADUDYvG7PRwPSgQjTNE6k3dPD/zrvtJHv85lIjfKls/XI1jxACvU8A6gmzSA71Y27KeiCMN8MX0gXENHcz62A8/V7F3Fe6izdFLlH923AKM/4+QfgMtlDXdWpdTzdU12E0dLFl+1CjenLmIQIPVwFQG+GHb2cPsyeYTtMRn14pQak0Z5oASMXMkfZgZUyj0+hUEkI00HDM+8bhZScsMzwHsxk4/hLjCr6cEP49XF302xMkEACEELcJIYqEEEUdHR2nEK5C8WUHByzE2g5h6FXzXsVQY+k08SiTgD6Nhhen3UJq314GNL4EFzewaXI84wNP/xSNXho1jq+Y1CZ50lT6WmJIG1fFhVkrAejreph/Da9fsHc1zwabSR0cwNRn5d7ghSTbPsRcL9G5LVTWNZ3050spcZduo6kzgZrDIaj1DuYF7CSxa2jbtsw0NO0eFheMjjGMjhebkEvXuKOjnrppb1XTZ64Y4ajOHqeSEE70QX18ZeXXlhFCTAYsUsoDx6y/RkqZBcwc/rruRD9cSvmslDJfSpkfGqrMoar4djbUrSV7fwm1rVPZ0ZgCmLmmqwyAsonxdAQHsSToE9Q7wtk0IY6L0uZ9b7EIIXCfoFdQcEws3S1NBC+6geAQB+NDy3E5c3k+KgMHkFTWSFBDGlbNIS6xWumsdVKXOAFL0ycUZE+k7uB+7Hb71/7s4j4LYa4iysqvAiA0p4lJ1naC3C6agapp52JX+RMZODrr5sPCwhkICKBHqPAHggbLqSxtOOl2ilNzKgmhEYg95n0M0PxVZYQQGiAA6D5m/ZUcd3cgpWwa/t4PvM5Q1ZRCcdr1OF0kt7yKV9Nk3qmKwu3WkBH0EJdLiRt4/eJbiTYfQqVyoXbFsifUzE9yF39v8aSF+52w2gjA2z8AL10uMel2rsxcDkD74IO8ztA/65UlO1geXc30qnKm+Wh4sCuUKXG9bHj/fdTR6WzYsAHPCbq2HrWyvIJA324+WT30gFdGSimB+4faDzZqBLboTIRu9AxXcTyj0YhaBfX6oSFEMtjEnu2jqzfUWHYqCWE3kCKESBBC6Bj6cF9xXJkVwA3Dry8H1svhLg9CCBWwlKG2B4aXaYQQIcOvtcD5wAEUiu/BqvLVxLZ1M9gVx5u1GQDcMPgiWmB/RAiHwlI51/A+cr8ftTMXEmMIQKM+9aGtv6l8UyB76rtPuC51ynQO7dxKbcZ1zMnYSIh3J4PmS3gqaOhBtbRtVdh6M3GIRuIDAsncXcHvopYyly2460rJysqisLDwhD2OBpwu7EUfIrwuo/KAHxpvB3P8txJQNnSnVBnpS6jZjM53dN4dAKhUKowOC82BQzXSGezlQInvCEd19jhpQhhuE7gdWAUcBN6WUpYJIR4UQlw4XOx5IFgIUQX8Eji2a+osoFFKeexYtV7AKiHEfoaelG8C/u87H41CcRyn2010/WP4VE1hiwyhrj4KnajgWvvQZDMf/+BiIruPEGoz02CawzrXIX456dLvNaZwfz3t/Seu2vEPCaO/q5MWRyya8ECuHrcCEFTyK9YAepeHS/fvZUXsQeZu3sLiay+krnA/m+bdxdz+zZQUFxMeHs7GjRu/cKfgcrn4c+FWJobVsHrdUJ+Q2LxqEv3tJPe04QYqp04lqMlMSmL493r831W4S0VzwtCwGhlUUXnEFymVkU9Ph1N6DkFK+amUMlVKmSSl/N/hZb+TUq4Yfm2TUi6VUiZLKeSHq0sAACAASURBVAuO/fCXUhZKKacct79BKWWelDJbSpkhpbxDKr9RxffgQOFjqGQohztDeedwFFKquNn7DiKBGi89W3MWs7DjY8z6AERYJhrZzcSYzBGNOSA0DNnXRe7P/8YPCl7DVzdAf/fP+If30JVw6qeH0ajjcHiq6Kxr57q0YN7aUU5R8AKmyVJ6e3ro6+tj+fLlNDU1sX//flavWYPa0kVEVAwvvTzU5Dcndh2p9R3okOwCGmZeicWlYl7e6GxQPio70ERdVhIAGXRQ1xqK1ao8sXw6KE8qK85eHYfosO1Fuz+GJj8je8qTAA+3WtYDsGJ6Pqn1dYRFt3DEMJOtWgPTjXFnJLRwfz3NvdYTrgvMyMOnpRyDfxxBxn4uMm0CNGwy/IwywGi1k7+vmveiGplVuJHx8xcR0tlBoZ8fW1oTmC53sXjRItLS0ti4cSMtLS044hJIG/iIFSuXMGDW4x9pJiXwAOweGv5il58WL79gBrS+FCSN7s4bGSm5tOSPA2A8DpqbfOnsUHoanQ5KQlCcnZw22rY8yYBHz4fmZCrN0NwcxTzVY+RJF50ICi+5lZy6PbhifPG4wul2lHLL9CvOSHizUkLZdOjE3ahLOz1E6FxIKQm/+G4uLliOXmvF3PV7ntAZAEh4q5qYCC12+wGKVm/n0TtuxlpVy8v6YMpEDl6b/pfs9BSuvvpqpk2bRmlVKWFGFctXTQYgP2sPicE6THVDDcpVqZGEme2g1mLQfX/tJ6dDVFQUg4GBdA33NIpw7aR4u/IswumgJATF2WnLo3wYlkZ5RRRJ3t1sLTcBcB8PA7A8KYXsqlrUC6pot2VTE5FBjNtFqN+ZeUI3LviLcyMcq7Ktn8wJOTRXHiRo0qVkRB3gqoT1gBfvGX5BOxDf10vgJj1rxvUwc9cuKsx27r/+HPrrKnj4sGBP3BWw9vdQtY51Fhc57o+weO6nco8eL4OdBXEfk9LsJMTloBrYdv5VxLV2oQsY/Q20Pj4+CCGO6Wm0gX07A0Y4qrODkhAUZ5+KT2g1pmE+UkxMr4MBly/FB5PJYhsLPJ0MAiuvuIm6qAQS/aC5N4YubSfnxc0+o2FqVCrsri82nUkpcbk9JOUVUFO8G1RqApNSmTXxAH76fjrNv+Nfw4Pepb2zCR+9Dwd9NtL41mdMiU8iJjkBjauZdzY28YB5Kp+W1rK/8DEq+oL45UNDP2vqhJ3offvxbNkCwCd6Nb6JBbjdTsKjRneD8lG+DietgYEAZFDMgQPGEY7o7KAkBMXZxdyIp3Ev9xxwk2J24K33pqLJl772aO4RdwLwdkAwYd5usr030OdIoiFiHN49hzhn0pwzGuqctFA2VHxxCK+y5j7GR/mjUqvR6Q3YBgYImnENycGbuTXxM8CLf6nvxgyk9/QQs85GXayaMOsqnl29kfEhqczOTabB1sOihCzeDDJh06rRh9xN9e4kVGo3k1MqyQhMIPHwUHVRTVIkXtJApzCQm3L8IASjU6RLS2PC0HiYGVRxqNGoTJZzGigJQXH2cNrwbPobD/QvYZbmIzxdaloGA1lZFUsc9Vwhd+MC3ltyAd3BkRREV1LWkogmLhqTCEWjVZ/RcHNjjRQ3fLHue0NFO3OGJ7YfN30WFVs3oktbQKxJxayJzUQHN9Niu59HNP4AZKxYQ0/NJA7pOklsep6Kzje5ckIuSxZO5Jmm95nf+CRpxjD++AcrUgqy8iqJ9G7C/9AugtxODgAHl84gwuzApdGRHRt4Rs/BtzUpKo2G3ERgqKdRfWs4ff0HRziqsU9JCIqzg5Q41j3E4/YLEINlxDqcNLtDGLCoaanI5n5+gRZ4W60mKjOQ8b2NCE8wVeGpaDu2ccWEy894yEIIYowGajuHBqWzOd1YnG58vYYadQPCIjB3tCFVaoJMSeg8pTyY8jKg43HXP2kGEvrtTD/0MtUNU2itUbPA3cRfN95OV+dfmK9uYnf5BP74mom2XclovRycl32QccEh+KzZCsDbfgYcqecS1WXGqdWREjb62xAAJqVn0jRxqHvweByYe7VUlVePcFRjn5IQFGeFvs3/4qXOdLJDdMSrC2kwNyP7ovmw0UBUdx838QEu4OXJs2iPjKEgpoiSI8n4p6WjMztITI4ckbiX5sfy4tZaXG4Pz2+p5Yr82C+sj0oZR/OhCvzyLyM0qouYzCBm52xikKtYphoapC/3vSOYAzbT12BCVBoIKpmK8UAOZmsOmkXn0/HeXACWTK9luk83A9pPyWrvxAq0TMlG4wrDy2kmMCwYH6/R3cPoqLCwMGxGIx1ChS8Qz3qKNo2N2EczJSEoxryKDW+wsc5CRnIyRft2kThYzWDdZPpsRuqax3E/96NB8goQMzeGNHMXId5udkRNIMZVztzARXzFYLvfO71WzU3TE3hkdSW5sUZMIV8cNiJp0mSqi3ZA3FTiEwxo3If5Y9Q7+Cf28pznLbYBIS4317yyl4PhxfTVRDK+rYRBszfWVgurfhhGT6+e+Ngufph3iHqfI8S8MTQlyWsqgVgyjgCLhw7hQ0zMyCTFb0OlUiGljlrvoZ5GOayntHh0Pz8xFigJQTFmOVwePlj+BuaKImjz4HjtNcZpt1JU74UPmbxrcBBW5Mf1vIUTeH5cFl1xmeT51FHSEEdMcixtDfVMm5wzosdhCvHhviXpTE8O+dI6lUqNX3AI5s4OQtLG49B14Eyewl1Tn0Hvn8Rt3IQTKKjqJLqtg3W5a2lt8SOkxcyHz15GxeFYAn3t3HdZOfsdOwg3fkp+RQsOYN2E8VgiC4jssOLS6/EznN6pQr9vXjY1DaFD3U1z2U1FjT9u94m78ipOjZIQFGOKu7+fgc2bWfvY8zz7v3/EVPkprdp08vMnsWHBYjr19XgPzmGvKpDBFi1/dP8aNZIXgORZmSTp9xDgbWNt7HQm+vWT5S7AN3B0fxBmzJ5P2cZ1qNKWEB08gNa7g3nODqbdsZ0q7z/yJ4JRATe/UoT1QD4fHong569fw/aDYRi8LCy69l163e+RFbaRvL9XoAZeVgmSFk+h0xNLsKMHjY+BxNDRO6jdicSr/WgYN3RXM4EqqltC6e8vH+Goxjal0k0xqkmPB1tZGYNbt+GxWjmiN/KJVzzTEr1J7xZ4cn7BOcHBFK/4EH3YKpwVeWiZSLnfYaY878WFfEIf8GhoOLk56eQPHKJUHcuE8QYqGw9wlenWkT7Ek9IZvBFChT1wHMnjQ9hR3oLeNI1balZgvs/NU/94ghnt1zLP4+Cel99mESUM4EdsUCf3XvkCk4J20dvbQUJRHfFNFrqBT6ZlcE5oLLXdAr1KRURkCOMixtZ0lAvS0llekAUri8ilm6bWcJrq1mPMzR/p0MYsJSEoRiV7bS39K1fisVjRZmVRMmUJu1sGCdXYmN69E8NAHxOv+AWavj6an3+WJ1O6mV3li0ssZJuqD9+uWu7pXQPA/wKTps9kYsA69I5wNgfmcbMfDNZNJGHGl6tpRqOMOQsoXb+a/KQcVHt2IwL6iQpIYtLAAcJ/Bn/+YD7JxWuZSg3rdBn8/JKr+X+WSlK086nrcxOjbSb+g6FZ0X5t0DJjxlXsjg4iqMJDHVqSAkO/1H4x2k1KT+Ppw5lYgHg8BHgq2V2oIiN3pCMbu5SEoPheuKVkX7+FikEbdo8kQqdhor8P4V7Dk5lICS47qHWg+k/NpWVvMf2rV6NLSKD33Mv45GA7lq52YgfLyPdRE27ex7i0BHS5P8ZRX8/B15/h31kuZu7tQOtOptdipz6un2v+3EYmZdQALxgDOW9aGpE2K2t0M7k6vofyzn7maZei0Z3ZZw++Lf+QUOyWQexJ55EVUcihSImwmFhUu5onp1+A6rqDPKCL4A87WylwNPDZ8kdYWxBNg+UgGb1tJO7pRQX8TQjk9ReTbOtnryWRcDFAr48OlUqgVY+tGmS9Xo9Qe1Gl8yLbYSeHjynenc+NIx3YGKYkBMVpZfd4eLW5ixa7k3x/H2YYfdEh2b53L/9XVoTFZSPK3Ue6qw3h5Y/W258Eo4qQgT56t9fSmzuDTSkTaevqwbhpK9lxISSNTybKV6DZ+SRMuwqichncsYM1u9+iMdtI7IZBjIk91O6bwbZIXy7Z/ALXdg31s/8xsHjqeZwXvpk+TRxalQe3zofsrixSCyJG9mR9Q9kLlrB/60YmJkfTv7ecsJR59PldyOXr3uCjaVfTenMbf/Ir4ZK1LSxweVi67T9TS7qAB1WC2gUZXGaYx/opHkIOG4h39RGeFo1z5A7rO/E4vWkI9CW7zc4ENlN6+BzcbgtqtfdIhzYmnVJCEEKcA/wdUAPPSSkfPm69F/AykAd0AVdIKeuEECaGJtWpHC66Q0r54+Ft8oAXAQPwKXCHPNE0T4oxo9Zi59nGDm6JCSHZW4/damXdG8/TcriC+J56LjYYKDfNoDjIxIqgUHL6jrC4aQ2flQbQ4h3KYOQUZKeaKRRxdawd3/iJkJgDB1dAdQPMfwCPW7D/n3/i46Aj5CUW4CwsI2bcbhoO3EBtcDBpvatYsqEKPXZeALb7+/Pr3AAcGi8KmcMt49sotQvizJEYw8fWh4ZfcAhOuw13+nwyDiynOcKGqgYCUmcwv2w965OupXWpk83jA/jg9VoyOu2EAYeBPRE60jMymDzhZnq9WlF5ArFpNBz0CG4uSKe4eWz2zolweNOQGAJtXeRSzjsNYfT3l2M0Ku0I38ZJE4IQQg08BSxkaO7k3UKIFVLKY5vzbwF6pJTJQogrgT8DR8cRrpZSnqhW75/AbcAOhhLCOcBn3/pIFCNqt3mQdV19PJgcjdbjpPqJe9lT0YqfoZesiC7Wpc3A0OEgsWkbeQe0ZNg0tPv48ETGIkjwI9rfxvi8iZjb2jnQ3EhlwyDZ9dvJW/13gmJTUEXnUfjWn9jed5jwxCmk7I+joXoL7rReetvnYtc6sKnauW3lx6QMNNCMmrtwc1nBUrLHl7HbO4Nz6ms44AphofcFaBPGVjI4Knfx+RSt/Ii8MBuH6+pIO38xG7f6MmhuZn59CSvVl3Mk4wP0D0Wh2+RGWmtxRemY3hZDXfwMrtXo+GN2JMll3vgzyEGDjg6rJCt6bI4WOsOUwN4JqbC9klxaaGqLpKF6FcY8JSF8G6dSaVgAVEkpa6SUDobmRr7ouDIXAS8Nv34XmC++5kkfIUQk4C+l3D58V/AycPE3jl4xKhSZB9na08+vEiJQNZVR+sML2W7uwH98C9YgDVtt85noE0Fknj9r8iazMz+X5ggfGrwCuHz/OuaWrWLQ1Y++/Dkc3m0E+RUT5beFYn81f1JP4u5DKu5Zu5vddb3EdkQQ/uY2ROkmCvPSCNcbsHWFU6z15po1z5JX2YALFVfiRhsSwZIMLZXBEfiUR+E3zUCEbyzd5S5MOWOjMfl43v4BeBuDsJsmkGppobVhD7F6PYmxBQwaormsfA+G+pnonbGUzXXxcU4egUcysCdM4hc+CTydYGF2Ty0dmmC83U4y4kMpb+4jPdJ/pA/tW5mXnUrD5Mm4gXQc6GQP21aPjXah0ehUqoyigYZj3jcCk7+qjJTSJYQwA0cHlk8QQhQDfcD9UsrNw+Ubj9vn2BhmUfEF9VY7qzvN3JcYibNkPZWPP8TerETCK7qwdgTTbUwnxv8IH3WDrLXjcm+jBieGaEmQbOegyoDGrWfW1s8Y0PrhG/8aA043R1y+DGp68MOXjFATvlFJuDZuQ9VwhD2pSez2S+JKz0o8jfls85LctOUlZhcfAOA3aNiMg3umX0PJQohvtOLvcXHAreanQUuoD+hGPcYaUI+VPX8xm16sYZLnID6+04nLDGXDRivplmgaAsIItNZg2qzDxxCJWach9AdZXNwQzesRdnx8OlHtScMYLSlVq/j9kgKe396I/gwP7He6BAYE4DIYqVWrSXa7yeAjivfEjHRYY9apJIQTXekfX9f/VWVagDgpZddwm8EHQoiMU9zn0I6FuI2hqiXi4s7M9IaKUzPodvOvhg5+nxyFY886qv75B/alJ+C2BoCpm74wI37xzbxZZyRicBAfdzQ9tf5cHVGKv18rjkAbhmqJydaJOdSHOkswA9UhoNLgZwpF6qNw96voLWlgcOdhalNSsAbEUqoJ4drIj/BUT6K8p53rDpYyq2gXXk43/0ccf+EIyclTCZ0XTcLAB8j2LMonWrgl8384tKGN3AVj++9IpVaTNnshXWv2YeqoZm97K5eF5POkWnLJ4R6qjEG8khmFVuPFdVFeWKo6+EOSmqD+epbsCeeT4DAKWnrZGKYjwHdsVp0dJYTAZQ+gNsCb5O5+CviM8urfYrO1oNePnaE4RotTSQiNwLEjbsUAzV9RplEIoQECgO7h6iA7gJRyjxCiGkgdLn9sGj/RPhne7lngWYD8/Hyl0XkU+Ud9O3fGhyMaqqh7YhmluSbaPSYmOwrpSA2hzOFFT5GLvF4nNrsL0VdOdHgf3d1u1K1gkG76fNSscUwl1F/DpKXXYUyYyK6173FkyyYCe+oJ1fahnmNij683+9s0OLwDmBy2j95NU1G3bmdOs5mFuzahczhYrU/lJ7ZDCLUfD85bSmDEGzR3TKRer2Xa+LkEE0aLugmdYex3rotMTmNP2QJ8i59h4g9fpnzbJiZajLyfEcO1+21M6NPTo7JwqFlPiN7J/H3raTVOptbth1e0kTJPF9fOzaStz064v36kD+c7CXT5UZ8YCN39TGYP7zeG0ttTRETkeSMd2phzKv8Zu4EUIUQC0ARcCVx9XJkVwA3AduByYL2UUgohQhlKDG4hRCKQAtRIKbuFEP1CiCnATuB64B+n55AUZ8LKDjO5ft6EOmwc+fWtVEwPon3ARJJrF/sDTWyv6iemrY3gPh07zDX0tNYx2DeIRjoJ8lUTGh5ManIEKREhxDc6kLU2djzxL2zSjVOrxennTVNqJLtlJl67+/F3DnKRvhKH04J0B6JS7SL3YBUzS/cCUJxg4tzaWtzA1VNvIC5rLdWWWMJFNdsCp7HYtJi9q+pJn372XDVOuOByyso/JGrtq4RnzmOy6wj7Br1ZO3MCCQcqcSPxUx2mblCLT+QCUg9r2DormHP2HuIDo4FJuemsPNBKTszYnm1salQYVfkZUHSEyTTQ1hlOdclBJSF8CydNCMNtArcDqxjqdvqClLJMCPEgUCSlXAE8D7wihKgCuhlKGgCzgAeFEC7ADfxYStk9vO4n/Kfb6WcoPYzGjG6niz19g/wmKYrm239A7TQ17ZZxVOoPU2GD2MONGFvNLN97iEPVfaQB5wI5DF0R+AFaejFTTY9KxUCYN/4TdHhHx2CNCaHDmEqrOYM+syTA2YmvtwGdxxudn43QjnCCGj5h0raD+FqtuNQqtswIZ+HWXtw4CTRexS/me3EkaJDIgVbe0E7j7oV3YO134HF78AkY3eMWfRMqlZq0O56lcdlsoi+4FdtAP5dvKucfnmyip+YSaoM4lwqVWkXjln1UZneTa0ths1cd100bmlymvKWPBelhI3sg39HCTBPr2mZg/9dnpOMggHq2rotk+pKRjmzsOaV7Zynlpwx1DT122e+OeW0Dlp5gu+XA8q/YZxGQ+U2CVYwO/zrSzk/jwuh96WlafeupdE+gSN3J4kYLta3w7J5m6ssPchVDXc6yvm5nHg+0DgxfDgxdK7hUhfQEBuMM0YMeXE4NclAQ3NWJ70Df55vWJaXxygwbf1rehcs1ACzgqWty6Er+GLtwsC84hfFHriI9NohdH9eSPefsa2zUefsScuED1P9tKTF3vIFQqbltzV6esGaTpNPzwNxoDr65lt3pWlKMSQRX7KDax4cp0yYB4PFINGO4gR0gKiIcpzqEg1oNuU4Xk3iLvftm4PHYUanOnguAM2HsV6YqzqidvQOk+ujxPlLPkdUvsWVBJPVtXlzR0sgrrb28++k+Jg8M8BlDdwMALoOWmsxxtMRMxZIQyy7vIqItXUztdrO/2R9zZTVeTfUkuJwkAzEeD6FdHUOPOB7H7q3DnJLJroJFLOt6mX2vteNyuYC5nDv1DhKy/06J9yCtA6EsCniIklQXve0WdHo1el/tmTtRZ5D/jAtxb/+U9vceoj98FpH5UdxTW87THdH8/LNOklMDuSW0gxpdMKsGJL/92TyEEAzYXXh7jc3eRcdSqVR47MHUhHqT29zHZNbyQc3l9PWVKg+ofUNKQlCcMpdHsqK9lz8kRdJ4901sXaSnoTWKGa1tPFHTxNrPSnnE7eb24fJHNCbeDPsp5bmLuDSunsgAG6AidzCdiuDlPJJQRV6CjqDp4ykMvJa1NU2UbNpOffEBTBLigVBfHdHR0Xj7xxKYYKPSE8v+1n52vv03rOajAy7cip/PH/nV4vv42LeH0O5g7rv6Q956pYK5FyRQWtjItEuSR+aknSEBP/0b4tEfEp/eTLU6h5ZDFcz0kuRmJVNS08GHQTMxfPwut10wnaCooTulvfU95MWNjTmUTyZe60d9WgQ09zGZUv5yJIbm6veUB9S+ISUhKE7Z223dXB0VTO+Lz1EZ38Tu3gTyujw8e7iC4k/3sxbJdMCBloe4n4dd9+Js1kEzvEQOwdF1pM9bxYwJzYxvu5rUARt1xm0ckj1c21qHRR3HjNn3UJ3vYGfVKrbuW4e5qxsqa4Haoe4N7Po8Hv/IKBw9f8dmu5y75rzDiqhyplojOP/ad6jc3kFLkAZtsw1Tdghq7diuFjkZlY8PXpf8hsHCt0lPKyQ9z4s124txFdczL3MG7eufIX7qdNJnzf18m5KGXm6blTiCUZ8+8xLD+WxaPmw4xGTacbq0bF/pzfi8kY5sbFESguKUWNweqi12lnps1Kx8lXfP8WXq4WBe33uYynWlbEGSDLSoInku9WZ6pibxE7GTfRYVFUei6SqOpqvJxJZXfsT+Df1kXX4YU0QPGdYLSLU5We51kKagBpZalnNhbCOTpnjRdlE6A9ub2OqyUXFEj5fFisXHHyLHkxNtpGnNL9jTkkOuqYa4xS8QbMhk0eVPIN1q2jsshBn1dDcPMvnCs+ND72QMmRl4BpbQvnMnwbfcSmKmpKSqjbA9n5F8zgX4z5v3hfI2p3vMPpB2vKnpcbxYM5kOXicMD4lsZ8fuNG5S2hG+ESUhKE7JK82dXB8VTPuv7uSFBRaWVKbx0n5J9botbMJDPNAYEEfLRZeSkj6BquAw3E4nWR2DTIppwTN5D+UVgWzZmkNfYzBbH59IS34hnRfuRBOXiE4fQpA2ite8YulT1xA20IPJ10z7eak4tFGYZAiDTj2R3VWk+zRT+u457DmUg9Fg5We3/ZPQiEmce9kyVCoVuz6updIXUtqc5F2ZOtKn7ozymTIFXXw83S+9iJ/Thbami9A7b8Qr8YtJsdVsG/PPHxzLaDQirCHs89WxYMDBbF6hvPrnmM0lBAYeP7CC4qsoCUFxUj1OF30uN0E7trPas5ukliQ+rDKxY9WzbMdBPNAbaWL/RdfxRs40pH2QgYYK2gP0mCKM2D2+xLnauTrZj2sXH2T52jhWrIyhpmgODZVTmXDdLmbFtNDuNUii2p8MzwTajM3s8GjIbvRlhtWFM+4jYgVEdmt5543L2bTzSoSQ3Pmjx5g8eyFZ0xYB0HGkH41BjaO4lyk3ZqM5S66AvwltZCShP/sZAIdXV3Jh4pfvkNZXtDM3bWx3Nz2WEAI/VRBVicEs2N/CHDbwTs0jdNY+piSEb0BJCIqT+ndTJzeGGWn+3T/YMjkA46YYlr+/go10YQK6YyP55xU/oyg1kd6+NnqkD5cYvQlr92KzoY1fXZJDlNc4Pih9lw7LPiYk2Jk6M483/n0VxRVx7Hp6Jl0xrfx86Qtkp9dhsM8ixnMe1zS0MGBuwW7tQ9VyMR5vDc8UT+O5nUMdWX9yRSF33HcHxrChmb5cTjflW5up7Bxg7uIE/EMMI3fSRoloo4G6zsEvzYZ2pNtCXPDYHrbieAVB3lRMy4X9LcymhgGLLzvXeZMycaQjGzuUhKD4Wo02B14qFeLNN3gtsZZJJeP40YeDvEsZE4DO4CAev/F+dkZ50dYNC5uauNDtz76oJoJT13Gdn52WTZJ6NESIYPJ8ognQh+GKPMjCu+5mbek0/vbqtVQ3RnDHY7+mIMvKwkmtxPi3EJUQTtTkcVQ09VJVrOOt92KobPBHJSS//oWZ3z8yB5VqaFgsKSUbX6/E5fLQHaUjN+vsufr9LhZnRPBWUQM/np30+bLGHgvRgWdfslyUm8DqnoV0/esz4nESzz42b8vhKrcVtfrsO97vg5IQFF/r5aZObvdRs3Pdq0SYErj3wwX8xnELi4F+g55nfnonK8LC8O10cVvQGpIK2qkPbEOqJVvcF9DhcyE/TI9lphakxYK9pgNbeSc+U+PxSooi91ZffvoneOQR+OtfYVepgV2lCWi1JqLDHGi1ktaOWPoHh6p+IiM8PHhfLxmx3RR92oOXQYPD5qSxopeIxACKAyTXToga2ZM2igT66OixOJBScnRE+vf2NnH91PgRjuz0i4uPRmMPZa+PloWDTubwb/ZW3E5Px0ZCIs4Z6fDGBCUhKL7SwQErsQYd3U8/zkcZHtrenUlG5z3cC7gFvPjTm3g/NJlcWwUzEncT6tGxNtjEPuvVeAem8ZvseCYH+wHgsbsZ2NuPOiCSkNtyEar/DHjr7Q0PPAB33gn//je8/joUFQnqmv7TOyQ+XpIyq4NHfudNbnIQEITH7cHS72D/ukYmX5hAiw78W/qICTy7qkK+qxnJIayvaGd+ejjNvVbUKoHRWzfSYZ12KpUKH00wh+ONLCzvYDbreKP6zzRXPK0khFOkJATFV3qztZtf9nfwds1q4p05vFeyjuLhx4c/u3A2q9PTmNBzkKyIUoJ7M3gv+wJKOhz8IiOFm1IjP78itR3uwVbZg9+saNT+X90FMCBgKCncpEfkggAAGk5JREFUeSeYzVBXBw4HGP0cdB9qRKgFa3a10m8NIzPUj6bKHsydNrLnxtDscLK2pIl7FqediVMzpsxIDuHhlRUMOtzsru3m1+emj3RI35vJRh9KC/KgfCVzOYTDqWPz6iSyZv/nDknx1ZSEoDihnb0DTPQzUHHXr6lNCeHVZxN5necIBcpSI3nt8kWoqzxMTC4lqXkpj+Qn09k+wIZzphPiPfSh77G56N/YiDbSh4DzEr7RP2RAAOTkHH2ng3GJ2K0uQg/1sHVfGzstzeRMCCcoP4iX9zWiVgnuXpSm/NOfgBCCXy0eR3lLH4vGh581zx6cyIIJcazsWEDHiysx4SCNIrYXTeTGgSp8/FJOvoP/ckpCUHyJlJKPO3q58pnH+CSki8q1F/PD/t8zC+g2aHn8rh/RWh7ERYkbiWu4jWV5IRi629lz6RJUqqEngq0V3dhrzEN3Bb6np3rCy6AhISeUhJxQ7C43+xvN9NlcXF0QR6DP2VcFcjqpVILMMTpv8jcRnxCLjyuM7X5aLux3soSn+bTiQTprX8Qn+7cjHd6od3Y/z6/4VlZ19jFxzw46922la2Acffve4DfYcQOP/OqnNNf5Mj3sANGtN/C3CSH4DjSy9rJzhwcZc2H+rBbcHoznJpy2ZHA8L42aSaYg5qaFKclA8Tm1Wk2k2kjpuKE5vc5lDYcaYjlUZB7hyMYGJSEovsDlkWzYd4C4V19gXaKBTz/x5lVKUQFvXDKPGk8YQT5Wxg0sYXmmEWlt4NOLL0AIgaOhn75V9fjOiMaQMTYnsVeMfXMSwimeuwQPMIsmfOhn5do5OOwnGD5X8QWnlBCEEOcIISqFEFVCiHtPsN5LCPHW8PqdQgjT8PKFQog9QojS4e/zjtmmcHifJcNfSsfxUeDVbUVM+vhdjqg72b1jKo8PPkM4sDcunPcnnceAyod5jkR2pBhpopWPz1+CGjX9m5twNA0QcEEiaj/lil0xcqZPH4chJp09KoEXMI9X2bo/g47qF0Y6tFHvpAlBCKEGngKWAOOBq4QQ448rdgvQI6VMBh4D/jy8vBO4QEqZxdAUm68ct901Usrc4a/273AcitOg68AB9paXk9pezjpNKgX7n2M+kk6NlsduvwO1w8Z8m4rDyTEUazt5c9Fs9FYVvR9V45VsxHdKpNKoqxhxAQH+BFv92RU11GZyMS+ytzKGpvLKEY5s9DuVO4QCoEpKWSOldABvAhcdV+Yi4KXh1+8C84UQQkpZLKVsHl5eBuiFEMrQg6OQs62Npzfu4ILeSjY5HNR/2s7v6MADPHrXrUQ76glzeGNJzGSzdy9Pzs4jpF1F/5YmApYkoIv0OenPUCjOlLygcDbPmA3AxexBumD12oU4bB0jHNnodioJIRpoOOZ94/CyE5aRUroAMxB8XJnLgGIppf2YZf8eri76rVAuLUeMx2ql7Kl/0hpsQLtyDR+VpfGMbQ1q4KXZM/AOsHDYnkN6UDIbg63cNSmBlGo1zpZBApYkoNKdvd0YFWPT3AWZMHk6ZUAQbubzLoV7JtBR+c+RDm1UO5WEcKIPavlNygghMhiqRvrRMeuvGa5Kmjn8dd0Jf7gQtwkhioQQRR0dSnY/3aTHQ9ujj/HP1CRuev853tWP4xcVLxIF7PQPpfKCHEockzgfL1akeHNuagBzy3zRBBvwnRqlVBEpRqUYUxTRg35sCPMHYClPsmV/HM2Hqkc4stHtVBJCIxB7zPsYoPmrygghNEAAwzOmCyFigPeB66WUn/82pJRNw9/7gdcZqpr6Einls1LKfCllfmho6Kkck+Ib6HrueTb5eUg0N7Lf4k3Qmt2cg51O1Cxf9gP6BqPIcap5oyCVvEA3N1aF4zM1Cn2ScaRDVyi+Vn5gLBtnzgTgEnYhHYIVn/2AgY5dJ9nyv9epJITdQIoQIkEIoQOuBFYcV2YFQ43GAJcD66WUUghhBD4B7pNSbj1aWAihEUKEDL/WAucDB77boSi+KfNHH3Okr4514TFM+GQDm/dqecA11PD22C1XEN1rxeaIZkdeChM8/dwzmIzxvEQ0RqUZSDH6zbogD92UKZQAQbi4iJdYtTOTtvInRjq0UeukCWG4TeB2YBVwEHhbSlkmhHhQCHHhcLHngWAhRBXwS+Bo19TbgWTgt8d1L/UCVgkh9gMlQBPwf6fzwBRfz7J3L62lW1gX5MvCwm2s7Q3jj80fogaeSp5CVkIXRZ45GFN6yLW4uC88m4DFJsRZPjex4uwRHRvBpCbBe9FDPdpv4W/sqYihujICt2NghKMbnU5p6Aop5afAp8ct+90xr23A0hNs9xDw0FfsVpn+eoQ4Gpto/eAVyoxWWq2pGJpKWbLjQyKBTdpgvO4KZVvL+cQFHcKmmsf/FIzHmBI+0mErFN9YRv5U/uo4hO3pl1nIIWI8Dbz56U1k5j5IVMFfRjq8UUe53Psv4x4YoPnpv1Di20eFVwHjN22HDbuZKy20omLdsvNobZyKWeOiMyaPm+amEKEkA8UYNfOy6RQYI1ihUaECbuP3fFQYRWdrDdLjHunwRh0lIfwXkW43TX/5HaX+/QTZU9FUNtGz+wC3D1ThBn6/+HpirD4Uq+IJy+7kB9OyyIyJGemwFYpvTa/3Is8SyvL8oaFzf8LrWM1aPlr3P3SVPjrC0Y0+SkL4L9L0+O8p9+pmvCWR99Vh+Owq5Ke1OwBYFjqf6Qv62SJmkpa1jazsS5gfnXSSPSoUo9/kO68lcG4B24AgHNzM47z2oYmOlvVIl3OkwxtVlITwX6Lllcep6qshzZPBS9pQ0rZ9wqU7C/EGXlAnEvOAP2U95+KJrcKYtoTrE7JGOmSF4rSIjI9gSYOB55KHpg29i0eoro9g/Y576Sx6YISjG12UhPBfoPmDp6g9sIMY1QS2dloJPLCLOSs/IhIPhfhw6FdL6Gucxd5QB9lZau7JPm+kQ1YoTqsJ/+8n6JdkUwbE08ePeYxnX42hs+8AHkvnSIc3aigJ4SzXtPpJjmzYTJg6j4ojR2gwdzP//dcY73ZQiZp/XngXsapg3g8K4YKJJfy/yfeNdMgKxWkXl5PK+S3hPBw91EHitzxIfXUw762/n6Zt/zPC0Y0eSkI4ix0pfJTm97cQosujrKmcfcKLc95+lmynnXoEt457lNkFZl7yG8dluVv48bRHlKEoFGetggd+RfLiZDYCIVhZxt3834sx9MgobBVvj3R4o4KSEM5CHo+Dw5/dS8dr+zGSxub2cso1wVz4+pNMsFppBS4J+itLr6/gRc1kzsvYy62T/weNxm+kQ1covjchmcnMdOfz5LgIXMDPeY7otlr+8PcfUN/0IrK34WS7OOspCeEs43B0UrHil/S/b8PXYeRDezu9eHP560+SZbXSBJxr+F+uuaOeN1RTmZdexU8LlhDgnzHSoSsU37uZjz7Awuzx/FUlUAH/5gdsXJ3Ey0W30bDjdnBaRzrEEaUkhLNIV9cmDr38WzyronB2tfGOyoeAvm5+8PozpNjsVAFztH/nop/38b5+ErPHtXD7xCRCQ+eddN8KxdlAGxTIOYtvpur8REqBVJp5xnk97/9zAp8MTqRz3Y/A/d/bFVVJCGcBt9tOTe0TtP7jM1RFCZQ17WC9MZu0mvVcu/xtol1udqNirvY5Lv1pN6tCs7k4o56fTQwnJubqkQ5foTijYm+6mhu85vCrWG96gUv4jJvqn+C9xy9mucObrtU/BqdtpMMcEUpCGOO6ujZSsetP2P/qpLG8i/UdFXQHZjB/9eMs3bQDA/AiBuZ7reOin7RRmJTBz3IquDLHRGzsDSfdv0JxthFCMOXJP3Jdzlx+6K3GCdzNo8zb8jqbX76BN62Sqo+uhIH/vll9hZTHz3UzeuXn58uioqKRDmNUsFhqaax/A9fKUFr2NlLftZ/ekBzGV62iYPduglxuBoBfkMZ7gW+z5IaN1GRH8mj6HhKT5hAWunikD0GhGFHW0lKe+s2d7PlsA6+4JBrgDZay9ryfMPnqFSTp9pGeeC1RE28e6VC/MyHEHill/knLKQlhbLFaG2g68g7ObUEc3t5Je8NO3L5hRHc1krVrNfGDQ7e6hQhu5pdoUq5n3FVFRKWZ+Z+MZuKSbsHPd9wIH4VCMToMbNnKU3/+LTtXFvKSS+IHHMbE09m/JPWGAPyS3iTS6sJv/FXkZVyFWq0f6ZC/FSUhnEWklPT0bKOrajddm4PZdrCNwCPb0GkCCGuvJrNsJ5E2BwB1wN2MZ7l4k5lLWvC+vI2fROwma1wqcXG3jtk/aIXi+zK4azev/OZ+Xtu1gaf7nBwdtGWlYR7FF88keKIJT9xGvHX7cRjSSU1YRH7MPAyG6DHz3M5pTQhCiHOAvwNq4Dkp5cPHrfcCXmZojoMu4AopZd3wuvuAWwA38HMp5apT2eeJ/DclBCkl/f1ltFds58g2PTvqnfjU7yKxvQ7fASvBDWWM7+76fEKLauBvRPNvHiMmLRPTtQdYmr6ZWeNDiYu/Dm9v0wgejUIxujnb2tj2s7u4v3wdsw+28mvAe3hdhT6endOnM5idhCoqFmf4brQ+e7HqvPB4GQnwjmZcSA7pkfMJ8ktFiNHXNHuqCeGkE+QIIdTAU8BChuZO3i2EWCGlLD+m2C1Aj5QyWQhxJfBn4AohxHiGptzMAKKAtUKI1OFtTrbP/ypSujF3VdNYcpDD+weoarRgOFJOYk0Rvp3tnDvQQ3R/D6Eu1+fbOIEPgZfI532smJI/YPal+7l+2mNsX1PLK6v6uW3z8Gimy5axbNky5syZQ2FhIcuWLaOwsJDCwkJMJhMmk+nzdYWFhQCYTCZuvPFGli1b9qUyRxUWFmI0GsnNzaWurg6TyURJSQl33nknjz/+OLm5uZSUlDAwMDRD1f33309hYSF1dXUA9Pb2AmCz2ZgyZcoXzsmOHTuIiIj4Qgw33ngjDz30EG63Gyklc+bMoa6ujhtvvJHCwkJ27NjBlClTKCkpwWg00trail6v/zy+Y8s5HA50Oh16vR6bzYbL5cLt/s8Y+Wq1GrfbTXx8PCaTiY0bNxIQEIDZbEYIQVxcHL29vZjNZmbPns3GjRs/31YIgU6nw+Vy4fF4vnBc/v7+mM3mL/wcALfbjZeXF3a7HS8vLyIiIgBobGwkJiaG1tZW7HY78fHx9Pb2YrPZ0Ov1n5/Ho+fi6DZHz/fDDz9MREQEjY2N+Pr60tfXx6xZsz4/RyaTiS1btqDRDH0c6PV6BgYG8Hg86HQ6gC+dm+/L0XMOQ+fw6AWrl5cXU6ZMYceOHdjtdmbPnv357/ro38nRc3LvvfdSWFjInDlzgKG/0aN/k8f+fR9dP2fOHJYtW8b6jCSmqV0Yja3MqSrhkg4zPwTG2eoZt64e1oFZ5cXBsFjaE+Owx5jwBAYiAi3UR62iPuQFXDrJoFaPTeeFVOvQanwJNEQR5J9KuDGHcN8Egg3BeGu9+TpH/1/PtJPeIQghpgLLpJSLh9/fByCl/NMxZVYNl9kuhNAArUAow1NpHi17tNzwZl+7zxP5tncIh7aU43G5kW43HrcH6fEgPRIpJXL4PdKD9IDHffT1cBmPB6QcXi6Hlw8v80jwfLnsf8p7wO1G5XCAzYpnYABnn5nBtnYsXd04+gfQWPvxtvbhZ7fg67Th77IS6bYRLF0nPJZOYA2ClYzjE67DEXABqVNs7Fk1ieKdvyUu+VwCAyehUg39cx/9/R795zr2+9H1x9/2HrvNV5U5tuyZuG3+qjjHyi379005F6fu+L/vo8uOvvd4PPRt2MTbf/4rr5VX4d9cy0UeO3OArxsQ3oagV61lQK1lQK3DotZh0WiwaVRItQep8eBRuxFqNx6NxKUFp1aNU6vCoVHj0mhwa3RIrZafrznIU5fOQuPtj04fiK9PKBMv/QGJsyZ/q2M+bXcIQDRw7DPdjcDxUX1eRkrpEkKYgeDh5TuO2zZ6+PXJ9nnaxMzM4Ovz8ehjA+qBWgSlBLKHBPYyk1b/mcSYYpk6uY0XFquZvUCPn9841GrILXhwpMNWKMY8IQQB82bzw3mzuVVKBquqWffMi9z90QEam7qJtjaQ4mllPDbiGKr6iAQCkES4HeB2AIPfKYafA/+/vbONseIq4/jvz8tdVqlbkFYoLwtEXlq/IBJSXyHSUEpMUYMGYyKISVMrif1gIg1Jl7R+qVY/mKgEhahNtbTal01DQ9HWNDEppVJeBWSpJG4hUECgBtgu8PjhnFtmLzN7Z9k79wL7/JKbO3POMzP/eebcee55zrl3Hnj2tR5l6/5+hMlb/9Cv/VYjT0BI+9pR2a3IsskqT0uypXZVJN0H3AcwYcKEbJW9cIzBDOMShrB4oJ6vtPI8Zb2th+VLwDkN4jziPIM4zyC6GEqXhvK+SpxRM6eGDOe90nDODbuFCx8Zy7lR02D0dCZMHM/06a3cPnYI88fD1KnQ3Hz5vObOncuixT3TFCn+u2I5rSxrmyybPHW1JI/OgYz7Ij95PydtbW2sXr2a4VM+zqLHf8Sixy/bd3XB3r3neOutd9m46wiH3znO2eMnuXTiBJdOnaR0/hRNXado7j7LsItnGXLxAiV7n5J1U6KbJrtAiW5KdpEmuilxkRfsJM9zOZVYVvRNmllGM5fGTqVo8gSETmB8Yn0ccDjDpjOmjFqAk1W2rbZPAMxsLbAWQsooh94rmJiRfrneSebyk/nW8jp4ymig4L7IT7WUUZ6JNk1NMGNGMzNmTACu7otqJcsSy5Wf53qRZzh8KzBF0iRJJcIgcXuFTTtQ/tnrYuAVC2fTDiyR1CRpEjAFeCPnPh3HcZw6UrWHEMcEVgCbCFNE15vZHkmPAG+aWTuwDnhCUgehZ7AkbrtH0tPAP4ELwPfM7CJA2j5rf3oDhzlz5vRYb2tr69GDaGtr62GXrC/PoqmcQdTa2sqyZct6tQFoaWkpfJZRUk95llH5fHyW0eU24LOMqs8ySvoLLs82Kn9GrgUapcV/mOY4jnODk3eW0bX3CwrHcRynIXhAcBzHcQAPCI7jOE7EA4LjOI4DeEBwHMdxItfVLCNJ7xL+0eFqGEX4K6BrDdfVN1xX33BdfeNG1dVqZrdUM7quAkJ/kPRmnmlX9cZ19Q3X1TdcV98Y6Lo8ZeQ4juMAHhAcx3GcyEAKCGsbLSAD19U3XFffcF19Y0DrGjBjCI7jOE7vDKQeguM4jtMLN1RAkPQ1SXskXZI0q6LuIUkdkvZLujtj+0mStkg6IGlD/GvuWmvcIGl7fB2StD3D7pCkXdGu8H/0k7Ra0jsJbQsz7BZEH3ZIWlkHXT+RtE/STknPSbo5w64u/qp2/vGv3jfE+i2SJhalJXHM8ZJelbQ3tv/vp9jMlXQ6cX0fLlpXPG6v10WBn0d/7ZQ0sw6apiX8sF3SGUkPVtjUxV+S1ks6Jml3omykpM3xPrRZ0oiMbZdGmwOSlqbZ9Bkzu2FewO3ANOBvwKxE+R3ADqAJmAQcBAanbP80sCQurwG+W7DenwIPZ9QdAkbV0XergR9UsRkcfTcZKEWf3lGwrvnAkLj8GPBYo/yV5/yBB4A1cXkJsKEO124MMDMu3wT8K0XXXODFerWnvNcFWAi8RHhA2J3AljrrG0x4BnxrI/wFfAGYCexOlP0YWBmXV6a1eWAk8HZ8HxGXR/RXzw3VQzCzvWa2P6VqEfCUmXWZ2b+BDmB20kDhcUlfBP4Ui34HfLkorfF4Xwf+WNQxCmA20GFmb5vZ+8BTBN8Whpm9bPbBI+9eJzxdr1HkOf9FhLYDoS3Ni9e6MMzsiJlti8vvAXu5/Ozya51FwO8t8Dpws6QxdTz+POCgmV3tD177hZm9RniGTJJkG8q6D90NbDazk2b2X2AzsKC/em6ogNALY4H/JNY7ufID81HgVOLmk2ZTSz4PHDWzAxn1Brws6R8Kz5WuBytit319Rjc1jx+LZDnh22Qa9fBXnvP/wCa2pdOEtlUXYorqk8CWlOpPS9oh6SVJn6iTpGrXpdFtagnZX8oa4S+Aj5nZEQjBHrg1xaYQv+V5pvI1haS/AKNTqlaZ2QtZm6WUVU6vymOTi5wav0HvvYPPmtlhSbcCmyXti98mrpredAG/Ah4lnPOjhHTW8spdpGzb72lqefwlaRXhqXtPZuym5v5Kk5pSVlg76iuShgN/Bh40szMV1dsIaZH/xfGh5wmPtC2aatelkf4qAfcCD6VUN8pfeSnEb9ddQDCzu65is05gfGJ9HHC4wuY4obs6JH6zS7OpiUZJQ4CvAp/qZR+H4/sxSc8R0hX9usHl9Z2kXwMvplTl8WPNdcUBsy8B8ywmUFP2UXN/pZDn/Ms2nfE6t3BlSqDmSBpKCAZPmtmzlfXJAGFmGyX9UtIoMyv0f3tyXJdC2lRO7gG2mdnRyopG+StyVNIYMzsS02fHUmw6CeMcZcYRxk77xUBJGbUDS+IMkEmESP9G0iDeaF4FFseipUBWj6O/3AXsM7POtEpJH5Z0U3mZMLC6O822VlTkbb+ScbytwBSF2VglQne7vWBdC4AfAvea2dkMm3r5K8/5txPaDoS29EpWEKsVcYxiHbDXzH6WYTO6PJYhaTbhs3+iYF15rks78K042+hO4HQ5XVIHMnvpjfBXgmQbyroPbQLmSxoR07vzY1n/KHoUvZ4vwo2sE+gCjgKbEnWrCDNE9gP3JMo3ArfF5cmEQNEBPAM0FaTzt8D9FWW3ARsTOnbE1x5C6qRo3z0B7AJ2xgY5plJXXF9ImMVysE66Ogi50u3xtaZSVz39lXb+wCOEgAUwLLadjtiWJtfBR58jpAt2Jvy0ELi/3M6AFdE3OwiD85+pg67U61KhS8Avoj93kZgdWLC2DxFu8C2Jsrr7ixCQjgDd8d71HcKY01+BA/F9ZLSdBfwmse3y2M46gG/XQo//UtlxHMcBBk7KyHEcx6mCBwTHcRwH8IDgOI7jRDwgOI7jOIAHBMdxHCfiAcFxHMcBPCA4juM4EQ8IjuM4DgD/B24+Zs8SI2kGAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#Iterations of Gibbs sampling\n", "iters = 1000\n", "expected_pdf = 0\n", "for ii in range(iters):\n", " #Given z, calculate parameters of posterior distributions for all GMM parameters\n", " alphaN, mN, kappaN, aN, bN = GMMParamsPosteriorParams(alpha0, m0, kappa0, a0, b0, x, z)\n", "\n", " #Given z, sample GMM parameters from NormalGamma posteriors \n", " mus, lmbds = NormalGamma_rvs(mN, kappaN, aN, bN, C)\n", " \n", " #Given z, sample GMM weights from Dirichlet posteriors\n", " pis = sps.dirichlet.rvs(alphaN)\n", "\n", " #Given GMM parameters, sample z\n", " log_p_xz = sps.norm.logpdf(x[:,np.newaxis], mus, sigmas) + np.log(pis)\n", " gammas = np.exp(log_p_xz - logsumexp(log_p_xz, axis=1, keepdims=True))\n", " z = np.vstack([sps.multinomial.rvs(1, gamma_n) for gamma_n in gammas])\n", "\n", " #We need to average PDFs of the \"sampled GMMs\" in order to evaluate the empiricl posterior predictive distribution\n", " pdf = sps.norm.pdf(t[:,np.newaxis], mus, 1./np.sqrt(lmbds)).dot(pis.T)\n", " expected_pdf += pdf\n", " if ii % 10 == 0: # plot only every 10-th sampled GMM\n", " plt.plot(t, pdf, lw=0.5)\n", "\n", "#plot the true GMM, ML EM trained GMM and the empiricl posterior predictive distribution\n", "plt.plot(t, true_GMM_pdf, 'b', lw=2);\n", "plt.plot(t, EM_GMM_pdf, 'k', lw=2);\n", "plt.plot(t, expected_pdf/iters, 'r', lw=2);\n", "plt.plot(x, np.zeros_like(x), '+k');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Collapsed Gibbs Sampling" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xl8nFW9+PHPmZlMZpJM9sm+p0nbpPtGy1IWgbayimUTBRQvinKv9ydeQP25Lz9FvSiKCqgICAKCQFmkQGmhK23apE2ardn3fTJZZz+/P54pxNCSSTvJlPa8X6+8ZuZZzjlPoPPNc57vOUdIKVEURVEUXagboCiKopwaVEBQFEVRABUQFEVRFD8VEBRFURRABQRFURTFTwUERVEUBVABQVEURfFTAUFRFEUBVEBQFEVR/AyhbsB0JCYmypycnFA3Q1EU5WNl//79fVJK61THfawCQk5ODiUlJaFuhqIoyseKEKI5kONUl5GiKIoCqICgKIqi+KmAoCiKogAqICiKoih+KiAoiqIogAoIiqIoip8KCIqiKAqgAoKiKIri97EamKYoyseDY3SEnc/8Denzsuqqa4m2JoW6SUoAArpDEEKsF0LUCCHqhBD3HmP/14UQlUKIQ0KILUKI7An7bhFCHPH/3DJh+3IhRLm/zAeEECI4l6QoSqjtf/VFzr7uJi645XYO/OulUDdHCdCUAUEIoQceBDYARcCNQoiiSYeVAiuklIuA54D7/OfGA98DzgJWAd8TQsT5z/kDcDtQ4P9Zf9JXoyhKyI2PDGMIM2KOsmAIC8OanUdfS1Oom6UEIJA7hFVAnZSyQUrpAp4Grpp4gJRyq5RyzP9xD5Dhf78OeFNKOSCltAFvAuuFEKlAtJRyt5RSAo8DVwfhehRFCbEj7+2kcM25738uXHMute/tCmGLlEAFEhDSgdYJn9v8247nNuBfU5yb7n8/ZZlCiNuFECVCiJLe3t4AmqsoSigNdncRl5L2/ucwYzhejzuELVICFUhAOFbfvjzmgUJ8FlgB/GKKcwMuU0r5sJRyhZRyhdU65eytiqKEkNvhINwc8aHtcanpDHS0h6BFynQEEhDagMwJnzOAjskHCSEuBr4NXCmldE5xbhsfdCsdt0xFUT5e2murSJs7/0Pbc5csp6lMTV1/qgskIOwDCoQQuUIII3ADsGniAUKIpcBDaMGgZ8KuzcClQog4/8PkS4HNUspOYFgIsdqfXXQzoFIRFOVjrvNINalz5uLxjNDQ+ABNTb/H53MTGRvH+PBQqJunTGHKgCCl9AB3on25VwHPSikPCyF+KIS40n/YL4Ao4B9CiDIhxCb/uQPAj9CCyj7gh/5tAHcAfwLqgHo+eO6gKMrHlNftwWA00tj0O1LM12A6vJCaN3+Oz+UFQMshUU5VAQ1Mk1K+Brw2adt3J7y/+CPO/Qvwl2NsLwEWBNxSRVFOaT6fF6HTMTJSS4Q5D8eecexzrXicFna8sI2yIR0Zbe1kZ2ZMXZgSEmrqCkVRgsLW2UFcahrdPa8Q3beafd5aPB7Bc+/G0OZ7hYPRqdzz0hb22UZC3VTlOFRAUBQlKHqaGrBmZwGChsoW0ubksH3XCC5nGB5nF9/2RpLvHuLn5S04vL5QN1c5BhUQFEUJCltHG7qILqLDl9I63o29MZyuri2Ed5UQ7sulrP1N5vrGMA64+EfnwNQFKrNOBQRFUYJC+nzYh0pw1CQSm5lM9b49rDqwmeWUElP1Kq74SswdTuLjDLzTYQt1c5VjUAFBUZSTdjR7SOKluqkBoyuRhMpXqLwxnt9d+S1ExhcxhHVglMnMadzHiM1J45hzilKV2aYCgqIoJ23MPogpNhyDPhq310Pb1vcYWTTKU7Gf5e5+N00WH/oePaMWG1F1jehG3LzSo+4STjUqICiKctIGOtowJwyjG8rBFB2Js7uJtxauJMdm5s2aFmodRgbtcRDfzPjoECneHqr7R0PdbGUSFRAURTlpA+1t6CL6aavyoR/W4YtqpM+biqVmnLRYB0uHxtisX4swdWElioSeevoGHYx6vKFuujKBCgiKopy0EVs/YeEGbPYRBg/WsefsdGKa04iO6qau7ve898az5Lf14BofxWCIJ8VVhn7Mw167uks4laglNBVFOWkSqc1hLGB8eIReUwLuMcm6x+5mQ5MdCfxx616eWp3JdWs/gd5cj7G/m/dsVi5MiA518xU/dYegKMrJ0w0hvHHoPBKnYRBXXzw3vXgvlzfZEWh/ed4JXLanlR01Wxgbs5I7fIjOIUeIG65MpAKCoignxe1yIsL7GeqIxNA7Ru1cI5mHO7jm4GGcwL2LLuSRL/0Mt07wX0Dj5m3IrlTy9e/SPeTApya8O2WogKAoykkZ7OzAGDNKb7OOsQ47NRmx3PjCbwH4rSmKguJLsHuTeOXiTwNwr8PN23uOEGboQY67aRhX4xFOFSogKIpyUgZ7ugiPNOJx6RkfHyGqzsXqji5swKE1adSuepOmvKfYM28No2EG1gPtpdsZ8OiJ6G/m4NDYVFUos0QFBEVRToq9u4vwyEiky8WwUbDh1X8C8JRBT/6lS4mJyiYhMouu3Fd4d/UFAFwzOsSuCifLB1+m3K4CwqkioIAghFgvhKgRQtQJIe49xv61QogDQgiPEGLjhO0X+hfMOfrjEEJc7d/3VyFE44R9S4J3WYqizBbH2BBgRNrHaU42sH7/PgB2LIojds4gJpzMj+slPszKv85aCMDngPL3RkkJq6VbPVg+ZUyZdiqE0AMPApegrYW8TwixSUpZOeGwFuBW4BsTz5VSbgWW+MuJR1sd7Y0Jh/yPlPK5k7kARVFCS+oHGB+MQGdzMsYgMR4vh4HkK+fgeVZH/KCXIWMCl+XW87tCSWdcNKm2IeIqmnH5oukfcSGlRFtNVwmlQO4QVgF1UsoGKaULeBq4auIBUsomKeUh4KMmOd8I/EtKqe4PFeU0IvV92LuNuO1Olux6B4BXjTpiuix0F68g/JrPETHfSsNoMTc2t7B95RoANrjG2VM7iBiz0+l0h/ISFL9AAkI60Drhc5t/23TdAPx90rafCCEOCSHuF0KEn0CZiqKEkJQSqe9jrD+CQccIF5aWAFA7L4bY4Rjme8dYkKdjzor5RHltxNmN7FyYD8AVwK6Do+QObKd6VHUbnQoCCQjHuo+bVuKwECIVWAhsnrD5m8A8YCUQD9xznHNvF0KUCCFKent7p1OtoigzbHTQRpjZgG/MR1+4h4yxcboAXVIKrsxIrrj2ehZ7trE6bYBPnOVlWJ9Klq+V8TA98wFnzSiFYzupGhkP9aUoBBYQ2oDMCZ8zgI5p1nMd8IKU8v37Qillp9Q4gUfRuqY+REr5sJRyhZRyhdVqnWa1iqLMJHt3F+HmCDxDo5j7mgB4E8hPtFKcE0NiyzNw6Y/hvLtI+uJjFEdVs2Ssg8Pz8wBY3DGIfryL1mF1h3AqCCQg7AMKhBC5QggjWtfPpmnWcyOTuov8dw0I7UnS1UDFNMtUFCXEBnu6kAYDxiEXcyr2ALA/PhKRlsgnMs2MrbybBx8y8tWvwp8eiyDssp+QmAIVRfMBuAA4UD1A77AanHYqmDLLSErpEULcidbdowf+IqU8LIT4IVAipdwkhFgJvADEAVcIIX4gpSwGEELkoN1hvDOp6CeFEFa0Lqky4MtBuiZFUWbJUF8jwhqJZ8jFWc11AAzkRbI23c2gr4hPXJBEVdUHx69YcTn/u/w+bEmRAFwE/KlmiKyhXnyyAJ3KNAqpgGY7lVK+Brw2adt3J7zfh9aVdKxzmzjGQ2gp5UXTaaiiKKceD114BqPoHxsm3uWmFYjMjWFxjI7137qJqiqYNw9uvRUeeghKSuCuob/yoxtuZjQ8jDynm9GaEayjFXQ4V5JhMob6ks5oaqSyoignzKfvx2GLRD/QCGjdALl5Gfzxza9xqFzPnDmwYwfccw9s3QpJSbCvNo+SQxdRN0d7NDmve5wU+wG1xvIpQAUERVFOnM6Ot99AYv1BAPbHRFCoi+C3L54PwOOPQ0KCdmh2Nvz619r7X27+P9Tl5wCw0idx1B6ifkw9WA41FRAURTkhXo8boRO4bMPMbdGeH7QWWPnHjltwOnVcey2sWfPv59xwA6xYAYPjcewb1XauBppr2qlXmUYhpwKCoignZLi/H324CZdtmLyRIZyAPt/Kc7uvRgj4/vc/fI4QcK9/NrSXD14PwHKgqnkY27CaxCDUVEBQFOWEDPf3InV6dPZBAEoBXfeFuDxhXHYZFBUd+7yrr4a0NKjsW0i7JZ5IwNQwzri9abaarhyHCgiKopyQob5uvB4d+o4mAPbpBDsPadnjt99+/PP0evjsZ7X3peYVAMwddOLsqcSrVk8LKRUQFEU5IUMDDbjGzaTUa2NKD8Yk0DqQR0YGbNjw0efefLP2ust2HqBNiSxr99HmcM1gi5WpqICgKMoJ8fh6Ge8zUNDZDECFaTEAN94IhilGOBUXw8IFXva7VwJaQBg9UkuDSj0NKRUQFEU5IVJvZ7xulCTnOENApe1TAFx7bWDnX3W1njJtuRSWAN0NbdSqSe5CSgUERVFOiE8MIttHATgEDDsuJitLSysNxBVXQA/JdIo4ogHZYqfBPjpj7VWmpgKCoignxCvdGDu6ADgk9EABGzdqqaWBWLEC4mOGKJVaBMnodjLS3zxDrVUCoQKCoijT5nKM4wViWhoAOKTLAnRcc03gZeh0cNGFY5ShBYRFPknv4d3Bb6wSMBUQFEWZtuG+Ptw+QXa3tphiqXcZ0ZFjnHXW9Mq57KrE958jLAV6a8qRKvU0ZFRAUBRl2ob6ehi1DZEzZMcLlHMeF68dnTK7aLKLLjJwiEWAtqTiSFMj/W5v0NurBEYFBEVRpm1ooAl9vRMDkiPAOEvYcIVl2uVkZYE9IRInBnIAZ3MrLQ6VehoqAQUEIcR6IUSNEKJOCHHvMfavFUIcEEJ4hBAbJ+3zCiHK/D+bJmzPFUK8J4Q4IoR4xr8am6IoHwPDQ40YW7T32jynxay73HRCZS1eWEsNuQAktvXTosYihMyUAUEIoQceBDYARcCNQojJs5S0ALcCTx2jiHEp5RL/z5UTtv8cuF9KWQDYgNtOoP2KooSADxvGhn4ADhJJRuoYmZlTnHQcK5Z2Uul/jpBtc3KovTVYzVSmKZA7hFVAnZSyQUrpAp4Grpp4gJSySUp5CPAFUql/HeWLgOf8mx5DW1dZUZSPAa+wYW1pB+AgeVxaVH7CZZ27RnKYhYD2F2dd6Y5gNFE5AYEEhHRgYshu4xhLYn4EkxCiRAixRwhx9Es/ARiUUnpOsExFUUJo2NlPdm8PABUs5fJzT3wtg7lzc6iPzAO0gNBVURKMJionIJCcgGMNM5lOXliWlLJDCJEHvC2EKAeGAi1TCHE7cDtAVlbWNKpVFGUmSCkZ7+wlwelgBGhlJRdcnXzC5cWnFDKa3QuVUAz0N6vBaaESyB1CGzCxdzAD6Ai0Aillh/+1AdiGlm7cB8QKIY4GpOOWKaV8WEq5Qkq5wmq1BlqtoigzxDE6grldm5W0GohJSSEuf84JlxcVn0BC0Sgu9OQA4y3tahrsEAkkIOwDCvxZQUbgBmDTFOcAIISIE0KE+98nAucAlVIbebIVOJqRdAvw0nQbryjK7Bvu6yW8S+s4qARyi4wQdeJ3CHqDgUWLW9/PNEpu6aDL6Q5GU5VpmjIg+Pv57wQ2A1XAs1LKw0KIHwohrgQQQqwUQrQB1wIPCSEO+0+fD5QIIQ6iBYCfSSkr/fvuAb4uhKhDe6bw52BemKIoM8Pe14ahaRiASqI4f0574BMYHce8uW0cRps+O6d/hPqhkZNupzJ9AY0rlFK+Brw2adt3J7zfh9btM/m8XeBPH/jwvga0DCZFUT5GervLiG4eAKCKXL5d0HDSZUYYIymNLoAhmC9h5/59rN2w7qTLVaZHjVRWFGVaWjr3kt6tBYQ6cwFLCqNPukzhS6AvPRHQuhUO7Hn3pMtUpk8FBEVRpsVl7yDNOYoTGMzJIaxg6UmXabHMQZendTvNBTpqT3xcg3LiVEBQFGVa4nvtANQASQvjEda5J11mTEIOmUvt+BDkAfZWNVo5FFRAUBRlWsI7tJTQKqBoCRCXc9JlRiclU7zIRjPJGICo1s6TLlOZPhUQFEUJmPT5cNUcTTmN5KI4O+jDTrrc6EQrRrODBmM+AGldg2pdhBBQAUFRlIC1Nh8kqkmbaKDOlMwyQ3C+QsLCTfikpCVOG+CW73JypLsnKGUrgVMBQVGUgO3Z+yKZtl4AWpNTKUwMD1rZwhePLU3LNJoL7N6v5jSabSogKIoSsO7GSrLdQ3gA24I0InLmBa1s4U2AXO0raS6wd9eWoJWtBEYFBEVRAhbbYkOPpB5IWx6PsBYGrWzhjSdz6SAAhcChisqPPkEJOhUQFEUJmO6I1kVUCWSkp0Bi8AJCRFQyyYWCMcJIBvrruoNWthIYFRAURQmI1+vB16DNclqjM7HAp4fw6a+jfDzR1iSM0RE0GLVZcJLau4JWthIYFRAURQlIW/N+4rq1B8oN0TEsCt7zZEBLPdUZoCVOWywne6gfny+gRRiVIFEBQVGUgJTVvEHemLZsSXtuLEVxwS0/2pqEy+2mP1VbfqVAuunoCHjpFSUIVEBQFCUgZQfszJH9AAwtTyM55cTXQDgWU2QUPqcFb6YJ0DKNDh48FNQ6lI+mAoKiKAEZ3hOBER+NQHROPrqkk5/DaDLhTSCmYBzQAsI7W7cHvQ7l+AIKCEKI9UKIGiFEnRDi3mPsXyuEOCCE8AghNk7YvkQIsVsIcVgIcUgIcf2EfX8VQjQKIcr8P0uCc0mKosyEiFovoM1hlBKXFdQMo6OEN56ofG3KigJg+87SoNehHN+UC+QIIfTAg8AlaOsr7xNCbJqw8hlAC3Ar8I1Jp48BN0spjwgh0oD9QojNUspB//7/kVI+d7IXoSjKzBq127B0tANQazCSpo8ES0rQ6xGEE54USbc+mmTvEM7GtqDXoRxfIHcIq4A6KWWDlNIFPA1cNfEAKWWTlPIQ4Ju0vVZKecT/vgPoAaxBabmiKLOmoXo3aaPaymit8RHMDfee9LKZxxJuNhMRE0WzRUs9TR5U8xnNpkACQjowcXLyNv+2aRFCrAKMQP2EzT/xdyXdL4QIchKbopwaBgdLqK//JfX1v6S//51QN+eEvLH9CPOl9k+3Iy+SRfEzMxNptDUJg8lAT3IWAHnOfjwez4zUpXxYIAHhWH8GTOv/BiFEKvAE8Hkp5dG7iG8C84CVQDxwz3HOvV0IUSKEKOnt7Z1OtYoSch0dzzE21khe3l3k538Dj2eEtrYnQ92sadt7IIZ5aCOHe4pzKIg1zkg90YlJ+BA4MhIAKMRHY2PTjNSlfFggAaENyJzwOQMIODlYCBENvAr8XynlnqPbpZSdUuMEHkXrmvoQKeXDUsoVUsoVVqvqbVI+PgYGdgGQlnYtwt+9kpx8GQaDhb7+bSFs2fSNVkQRgZcOIDxnPmbrnBmpJ9qahHssgrAM7aupENizU81pNFsCCQj7gAIhRK4QwgjcAGwKpHD/8S8Aj0sp/zFpX6r/VQBXAxXTabiinMo8nmH6+reSmvrpD+1LSbkSm20PXu9YCFo2fV6Pm+hGra2VQEpMzoxkGAGYo2NwjZgJT9U6EgqBdzaXzUhdyodNGRCklB7gTmAzWsbZs1LKw0KIHwohrgQQQqwUQrQB1wIPCSEO+0+/DlgL3HqM9NInhRDlQDmQCPw4qFemKCHU2vY42Vm3v39nMFlm5q20tT0xy606MUfKW0kb0f5KbwgzEhMWDfG5M1KXEALhjUfkhuFFkAMcOHBgRupSPmzKtFMAKeVrwGuTtn13wvt9aF1Jk8/7G/C345R50bRaqigfAzU1NTQ0VDHmGOPtugqWLFrIUmsCukmBwRSegsc7hsczisEQGaLWBubNNxooQluspi3WQprJF5RlM49HEIYp1kSbKY5sxwC6joYZq0v5d2qksqIESVlZGS6Xi/6sEcqKNhI3fwEvl5Zzd3k9ewZHPnR8Wuq1dHb+4xglnVp27nUxn1oAOtItzI+Z+QnnIuMstMenAZAx1jnj9SkaFRAUJQg6OjpwOp0cjk9g0KPnVkM7ix1lfGVpAcvL9rKnvI3Xm/v+7RyzOQOns4cPEu9OTeVV6RShTUXdU5zKopiZSTk9ymA0Yo6NZjhDCwhzfP04na4ZrVPRqICgKCfJ7XZrdwdz5lNR+zqF3a/T13+A0tLtbH/rXhIsB4io2MYbZVUc2PfvI2/j4s7CNvheiFo+Na9X4myMIAY3fYBrznwy46c9DGlatGmwTfiSYwEoRLJ7d+OM1qloVEBQlJO0b98+ilau4omK7VwgXyI6fQlvd5TQJAUxC28iJn01+XMPcVm0m9/0NDNS3f/+ufHx5zIwsCOErf9oZfsGyXVoA9IqgfS4AsQMTGo3UXRiEj63BV2i9pyiEHhnW/WM1qloAnqorCjKh0kJL7wwxjPPZHFoyMnChFHePXsJ1hKQ4+vwhY/QNbCP/UNt5KZnIH2PsNyzjvvGfHw/ew06swEh9Oh1plP24fKW1/spQgtYDWHhRBotkFgwo3VGW5PobgvHl/nBWIRf7a5g0ow5ygxQAUFRTkB5Odx6Kxw4EAFEAFDNdfAkzJ3XwN0bX8ATB5UDcTh9qTTW72ZpRhzFWXb+3N1H9eYKiq7WMrATEi+kv38bycmXhe6CjmPXnjE+yV4AOqPjMJl1YJzZwBUZF8e4TYcvW49D6EmTXporVOrpbFBdRooyTS+/DKtXw4EDkJDgYPHn9vCp27Yyd8VhTOHj1FTn8bXffImSTgsRkeOMiSFMrKF10EJr9zZuXzGX/x1sxjOiPSi1RBUzPHJqjsY9VBVLETUA9CUnkhflnfE6dTo9SD0REUZaorQpLCz9R2a8XkUFBEWZli1bYONGGBuDdeu6+M4TP2Pxea9TX3s/upFPkZWxCqv5S4wMN/PUr24k+/mD3DHwOtbYVlpHl9HZOo/a0l8irdHse3YroA3G0utMeL3jIb66fzc6Ck0tH2QY9RUksTDCPWv1R8ZG02PVptjOcnYgZza5SUF1GSlKwFpa4NprweWCr3zFw6rzn+Fnj7xD7Uvb8fkm/uVcATzCqPtmbi//Kdeu/DN3vvsqvoUGaihmWdcACSsNPNk1yFkeLzqDnoSE8+nvf4ekpPWhurwP2bPbRRIDxOPCBgwU5bPYmjBr9UfFx9KVlgINhyign+ZmBzk5plmr/0yk7hAUJQAeD3zmM2CzwWWXwSWXv8F3v/tTql/YhpQ+FsxfzO2X3sC5d3yDZWtXkaqDcB7DYbuYkpc/wf2Xfpccmx3fgKRuvIA5FX+lITmBqld3AmCxLGR4+NSazmvrZjvFHAS0DKMYax5RKfNmpW69Xo851oo3PgbQHiy//bYasTzTVEBQlAD85jewcyekp8NvfzvE1/7zS7TU9GCJieGmm/6br111B1nL4vn7rkfYu2s/HT6JAyijjktrr8TwchnvZOayLMbD8EA+thEL6bpRnuvQppTW5jwSp9Qgtd173BTxLgANxnDizNEzNqndZNHWJAQJkKxNs10I7N6tUk9nmgoIijKF9nb4/ve19w8/LLnhS9fQUt9GRHwi11z/f7nImsqKNx7m2z/7PRkH7eg9XlwRBrw6wWLg1wzxnZd/woLOXRyWlRTEj9Pbt5wFtZvZmmJlbMAOgCV6wSlzlyAllJVbKGIfAH3RCRgi9GCOnZX6Y5JTcQ2bcOZovdqFQPmBU+N3czpTzxAUxc87PMzgP55D6HXEfPrT6KOiALj7bhgZgauvhkMNv2fvm1swmox876rvkSk6ueSJP5DYZ8NtMlF+WS61awoZ8GTjsXnIqT1C4QvvMA8X6b98hr988wYORIeTbVqDyxmHdayfN1/YwlW3XUN83Nm0tT1JdPSiEP8moLUVBuwWitH+Kh9JyiTTPPMZRkfFJqfQ3VQLKRK7wUisx4W94eCs1X+mUncIigL4Rkfpvf9+Yq68Asv6DfQ+8ACOqirKy+GppyA8HO6+u43v33s3AP+17j+JEj1seOa3JPbZ6M9O4/WfruXQnCJGK0eRVS0YbILapBR+uf4qHsNIJD5u//nTzHW2cSj2CJbhBazqLOdPkVpev8Fgwev98CR4obB7l9Z1Nd+fYTSSncui8NnLgoqIiWXcPozZFE5rTDwASUMq9XSmqYCgKMDAE38j4fbbMSQmEpacRPK99zL89la+9GVtkfflN4zy2W/dinN0jHPmL0GfOpcbX/gNsYPDdMzJYv8dq+jdIrEdcdKr8xGb0sd8026W9o2wMiOfX2TfxSNAuM/Hf/3i7+R563k6HTJG0jAO9zPUNwBAmDEBl6vvI1o6O97ZOoaVHqy4sQNDBbksiouetfqPriMRGWPBnpQEQL6vlb7Q/2pOawEFBCHEeiFEjRCiTghx7zH2rxVCHBBCeIQQGyftu0UIccT/c8uE7cuFEOX+Mh8Qx1tJRFFmmKevD4QgLCXl/W1SCH6QuZHdu5IwGb3856Xbadi2hYhwE+ef/wW+9vq3ibMN0ZmRQsMXzqLycAI7r1qL46ocPnPtVyiYfxOlUek0r2shxVzLV+f3cAc3shmIGHVx1x+fItNzmBdjw1jTPcADW7QJ7uLjznl/6c1Q2rPbRxGlgLYqli89jaTM+bPeDktCLO5EbdbTQgYpLf14rDL3cTVlQBBC6IEHgQ1AEXCjEKJo0mEtwK3AU5POjQe+B5yFtmby94QQcf7dfwBuBwr8P6dOArZyRrG/8goxn7r637b9ua2P2ke1roqbL6zle3ffCsBV567hlsP3kdrUi91ioeWOlWxuSOfpNUtpGIui1ngOj4wf5r3OWooc88k+uBZDdi8LljhYl7GBzxJPG5BaO8AtezfRGh2FySdotg0BEBlZwOhoaLtGnE6oqI6giLcBaDIZMcQZILl4llsisCTk47NqXWqFwLvv1s1yG84sgdwhrALqpJQNUkoX8DSTZpmSUjZJKQ8Bk3Pm1gFvSikHpJQ24E1gvX895Wgp5W4ppQQeR1tXWVFmlZQSb19/rUvXAAAgAElEQVQfYf5uCemTvLO1me6XB9mySY/BAIvOfp7a9m4SzCa+mNNL4fY23Ho9FXcu563OFZQl5XFJQz0XHumG/Xb27C/gWecKjsQnIy5fh8v4KVpzXNx6TROD4kfc6q/7gqdLWT28mxfTzMzrtdFpH/J3lUhkCIfllpSA221gmUEbI9EbG4vBEgbmuCnODK7I2DjCTBm4M8MBLSDs3187q2040wQSENKB1gmf2/zbAnG8c9P970+kTEUJGmftEcLnfdAVUvZWKzsNHuoPJeDzCRatreRHf/w1ALdcksM5z2hfSOXrixlyXchQmJ4VyzsR4ZEkx3r5zTlO3vjBDdz1hfP4S1Qh9+3tYk1BIctjv4rnrF2sX7SMLSzjMcDg8fGNPz3BSEQU7cCzf98CgNmcw/h482z/Kt63wz8b90KpzWE0as0iWT/7C9TEJqfgHjbjzNE+zwFqK1Xq6UwKJCAcq28/0D9fjnduwGUKIW4XQpQIIUp6e3sDrFZRAjO6/V2izjsXgKG+cWxeD+Z4PU8/pc1gutr8LN2d/cTEGvlil53wERcNuTmMLFvMoxFptOdWI9+L5XznIQqXprDXUsTbb7/NQpeLvXduIOvcRVxa5qG/aTef8H6bDbc+DPyKuwAbkFbRx39UPcPe1Eg8fb1Ij4/4+LMZsIXuOcI729yAZK7PBoAvYzFF+tmfZyk2OZXh3iFEjI5ukxkToO9UqaczKZCA0AZkTvicAXQEWP7xzm3zv5+yTCnlw1LKFVLKFVarNcBqFSUw3qFh9DHa9Ag173WxL8fIe0+24hg2sXrZCFsPPAbAt5dmM39vJ26DgfJbo2lwFEDMfpa0L+I67zM4briUnuJ0jiQcoSK2gudan+Mf/3qGb+XF8Z3rzuKLuoXsrnydS3XnsXa5iX7W8f/8bbj60bfwxRg4ItyUv1mJyZSG0xHoP7Hg8vlg1y5BBm3ESi99gC5nAYtT0ma9LZZEK0N9PUQYzXTGad1VWa467PZZb8oZI5CAsA8oEELkCiGMwA3ApgDL3wxcKoSI8z9MvhTYLKXsBIaFEKv92UU3Ay+dQPsV5YR5h4fRR1sA7dmBw+XFa3Dz7rPaw9NF5u9T1dlEfISOL5Rrd6c7PrmMWNtGSsNz+eRABrda/kXqfzzEZcu/yvXzrueOJXfwteVf446z78CT6+GR6kcwN5Xwx8+s4dsJa9jX0cinL9wG/JjfomVjJLQPcc/Wh9mTnM7uw6VItxeELiTTWFRWgn3IwGrzZgDKBXhTrWTnL5n1tugNBnxeH5a4aEYStAywQpqpqpr1ppwxpgwIUkoPcCfal3sV8KyU8rAQ4odCiCsBhBArhRBtwLXAQ0KIw/5zB4AfoQWVfcAP/dsA7gD+BNQB9cC/gnplijKFsf37MS9fDkBPyzANMToGdx6isz6exKhB6rvfBOD/pWWS0DdIR3IywwutNBiyMKX8k41Re2jKu464uSs+VHa0MZqbi2/m3ovupdJUyZ69f+KnVy/jf+OXYLGYmFcYgYNP8V3/8Wuf3Ut09DhHxrux7+vAElUUkjUSjj4/OCfsRQCaLOGMJ3jQJc6Z9bYclZhehEjRZlktZJj9+4dD1pbTXUDjEKSUr0kpC6WU+VLKn/i3fVdKucn/fp+UMkNKGSmlTJBSFk849y9Syjn+n0cnbC+RUi7wl3mnDGVahXJGclQcxrxgAQCtVQO0phjY/oQ2vfLqjFfZ2lCBVcDNLe0A/OvWFcQPnkNr2vPcNTDEobRPsujqmz+yDovRwt1n383CBQvZWnofnzyniJe8EZy9sgX4AX8DmoC4/jHu2Pk4O+IyKampJNayCptt98xd/HEcDQiLPFpffX9yNOMJPtCHbpabuPTFuNO1/y6FwJ49KvV0pqiRysoZS3rciDBtIXf7uJvChlco230BADrTc/h8Pv5fnBWTy8OBxQtIFi7M0QMUerpx5f2U+IJijOaIgOraULCBL573Rbrbfo0oXkRc3G6iLVl4uYKf+Y/5xLM7MCZL9tsbcB924XHPfmf5jh3a32U5Lq2LzJGVQbYYmvV2HGU0mzEaM3BlfbC+ckWFSj2dKSogKGck78goOv8cQj6vjzqHk62vG3CMhzPXWsE7jVsoAG619eETgjdvWkLKyPn8Kbqa5Su+QfNIJXPXnDetOotTirnr/LtIkI+ze+7ZLF24H/gmfwW6gKSuEW7c/zwHHG76WnsQwojP5wzylR9fays0NwviTX1ke5x4gbC8xSwwh+5rIjYllZG+YRyp4XiBHKC9QaWezhQVEJQzkuPwB91FA51jYOjmtbe1L/j56U9jtw/zK2M4eil5+8I1RPbp6ZcOIpfNxVkazuLrruREZlspSCngP1d9idSkA+jPtSPEWTg5j1/691/xzJsM5MZT5m0nqn8pdntpsC55Ske7i85PehU9UAOkWM+mOG/hrLVhsvjUDAY62wkzmeiIikQPJAyVMzoasiad1lRAUM5IjooKTP6AUF9vw9R8kJaq+ZgNDqpG/sYK4AqXE5chjNc3ruK8kSweXWljox5MGXFExcWfcN2Lcxdze+5yHAVO5hQcAb7JQ8AwkNFkZ2XLXkrqKjB0JTFo3x+Myw3I29pMFZwt/glArcmAxxxP0pwPPzSfLbEpKdi7OomMiKAnVvudF1KjMo1miAoIyhnJNzqC3qKlnJa1DrJph5bFsjr/bWrqW/iB/6//zRvOwdzjpis8grS8uRgarOR/6pyTrv/SlZdyc/o4uvN6gPU49Av4s3/f9c+9SGVKDD1ho8hZnN1zizZQmsVjWhDqTIhgKNaLCI+avUZMojeE4fV6SEguwJWojUMqpJXDh0PWpNOaCgjKGW+gv4rtuz8BwFjkn1glJZ+UEldYGK9efi5rhpJ4dPkgF3eVkJf1KURYcP7ZfOHiz7JkbQ1R8UN4vN/kd2iTgRUdaCPG18PBoSYMDal4PDOfZtnYqP1EmUdIHtUyw0ey4nBZTo31GZLz1yCytQBeyCj79w+GuEWnJxUQlDOOp78ffaw28rVjZJShmmF6ezOINQ9T3ryd7/mPe3vdSsy9YxgiEjCFp5IZn45lWcbxC56msLAwfn3p5cStaQSupdWQxKtAmEfy6Zc38U53PXpzCgOde4NW5/Fs3aq9LiwoJ3vcAYAsSsUa5ZnxuqcmiLYW4crWMsIKgZIStVjOTFABQTnjjJeXY1qoPT/49dvPU12WDUBy9vMs7O/jk4ArLIznLr+EbE8Cz2d28NWsZCzuYsKsgaWZBiolJYV7Pj+K0Otwef6bB/zbz91SRn+Sng6LnrHSrqDWeSxHu4vOtm4mRko6gYTc+SzNCt2AtKMsCYk4RpyMpmq/+0Kgtlalns4EFRCUM46zuhrT/Pk4nU766/vZUaGtYWzX/Y3v+I/ZeelC9HYvGbpURuKziaaJpPyLZ6Q9X756NamL24AvsVVnpBKIGXZz3q7dvFGyE0YMMzodtpQfPFBeOaSNUC4LM2AJy2f+grNmrN5AxadnMNDehishEqdOkAq4+isYVgOWg04FBOWMI10udCYTD775Z3SNxfQNJxFh6SapYTuXAS69nic/eRlpTsn+mFbu2/AZPLZhzHOSZqQ9er2en90TDsSjN3yW3/m3r3/lberSXDgTJSO1LTNSN0B1NXR1QWyknbSORgCaEyPxyXRM5sgZqzdQ8WkZDHS0E2GMos2fCFBAGZWzP7PHaU8FBOWMIqW2+ExtbS0HRCsNO7QvvJjMP/HfDm3O/4q1+TjGzKQY8+kuKCKss4TYmFUnNO4gUJ+9NoWYFDsu1908AQwBcxoGSOjqZuvBdobKa2as7qPdRUvzqoi2uQEYzo/GHR7c7rETZbZEMz5sJyFuAfZ4bWbaQo6oTKMZoAKCckZxt3egS0qmtK4Uty6Skmpt0FX8+N+4CS3L58kNl5E30oMtUnLP8vPpa9hCyrKZXeFVCPjO/5iAuQjLeTzu3/7JTW/RHNmF09eNd9Q9I3W/9pr2urp4D3n+B8oszsCQqJ+R+qZLC8SC3EUXI9O1FNhC2qlQA5aDTgUE5YziqCjnCJKuhB502xYyNB5NZEIFn2usxgjULk6lWZ9GXEQejTmJ5EgjwhSG3hg+42374m3hGIwehoe/z4P+befsrKLPEsmetnrGSruDXufY2AcZRmebnsYiJc1ASlEaeQsnL50eWjGZ8/DlHQ0IDkpL+0PcotOPCgjKGWX4cCWu5GTK7ZW0vqalkCYn/5ov+/e/cNk6cgfaMJitfKFgHl1lr5FSNLN3B0fFxMDNnwG4kA5LDluAcI+XNVtKaBgfY6y5OugPl7duBYcD8lJaMFfXA1BuCifCYKZ4XkFQ6zoZRrMZt3Oc0ZSjXUZw+LBKPQ02FRCUM0pHexsix4PDl0tJbTEgubHl78QAR9Is7E1cRHJUPOWZcSyLz2DMVE1swuxN3fCVOw2AYNT5zfcfLm94dSe2LA+lHQdw1g181OnT9uqr2uuK3IPQpI05aEqMQTqTiQoP3ZTXkyVmZtPX0sx4ytGHytDbW8ugGp8WVCogKGeMocFBdHo9bzc9ieuNVThdJuJSX+dLI2MAvHPFeaT1txNpTGdjShzDJY2EpycgxOz9M1m+HFYs8+J13cJbJgstQGrfILF1I+wfdWB/+52g1SXlB88PLl2wi2Sb9vxgJD8OvS83aPUEQ2JWNn2tTegiMhgy6IkFkilTD5aDLKD/04UQ64UQNUKIOiHEvcfYHy6EeMa//z0hRI5/+01CiLIJPz4hxBL/vm3+Mo/um5mcPkXxO/zmmyQsnUuv10TZa9r/bjf6vkkm0BwRxuuFF5Eepqcuw8qalDnYo3eTnHrZrLfzjq/qgXA8pjt4yL/tkpd3IBLrONLrxNvZHJR6KiuhuRliIkdYnfYW850OPIDp3HSMiQuCUkewRMUlMGIbICVpGd2x2nOEYkpVQAiyKQOCEEIPPAhsAIqAG4UQk5823QbYpJRzgPuBnwNIKZ+UUi6RUi4BPgc0SSnLJpx309H9UsqeIFyPohyTz+dDVtewPeo93PZLaagvAFx8vucQANvOKSLaPkCcMYmChCi85TY8Sf1ERubPeluvvx5ion04Bv+Hv+oNOIFVpVUM6uMo1dlpfW6z9uf9SXrhBe11UX4VfSUt6IEyvR5rVizFS0+tB8pHU34Lz72M4dSjAaFeBYQgC+QOYRVQJ6VskFK6gKeBqyYdcxXwmP/9c8AnxIeTtm8E/n4yjVWUE1VfX0+Cd5wa8wiHn0nA6wnjyvjvskJKBgS8fdF1ZLrHGUxJ5sqYHHzZdiIjQ9NtEhkJN9+iAxKxWa/nH2j/UFe9cZhB33Z6hlPwVGw66Xqee057XZFTymCVD4Dy2FjChIGsOSc+vfdMCk9IwJWnzUNVTBcVFWrl3WAKJCCkA60TPrf5tx3zGCmlB7ADCZOOuZ4PB4RH/d1F3zlGAAFACHG7EKJECFHS29sbQHMV5cOampoY89Yz4jiLtgOpAPzH6MMAvJKXgkuC1WilL11PRIdk0LKTpKRPhqy9X/anPTn7v83v/dvWvb4Tkayj1eCifEcDOE58ic0jR+DgQYg0u7huzlskdmvPD9oy09ELIzrjqTEGYSJLfCKjtgHcOdpXSzFuysvVd0IwBRIQjvVFPTksf+QxQoizgDEp5cShJDdJKRcC5/l/PnesyqWUD0spV0gpV1it1gCaqyj/zm63E2UUNDr7qNmXSVvzHHI4yAanDTewc8NG8vqa8MXE8GlvHhFnWfH5nBgMoVsHoKgI1q4Fn3s+B9MuZD9gcTiILXdzaHQz0rmAkW2/PeHy3787mFdFiqWRxaPaNNeehXkYRFYQriD4ErOy6W1pZMSaBkAxWqZRvxqOEDSBBIQ2IHPC5wyg43jHCCEMQAwwMT/uBibdHUgp2/2vw8BTaF1TihJ05eXlJAzv4WDSYjylYfi8eu4xfRU98HKUmb7cXJIjcqmaE84CUyp2/T4S4teGutnccYf26hn51vsD1da9tJOYmHYGnePs646F1hObGvuZZ7TX4vxS6g93EQXUCUFCURxZyctPuu0z4WjqqT5qIUN6HXFAKu+p5whBFEhA2AcUCCFyhRBGtC/3yR2Ym4Bb/O83Am9L/wgaoeXsXYv27AH/NoMQItH/Pgy4HFAD0ZWg83q9OBxjtB7Yx8HRuRysWU4Uw3zGsRuAV1asIb+1CoMxnALCiDonncHBvcTGrgxxy2HjRsjKcuEaupgXrcX0AwUtbbgdiex17sYqltK6+5/gnd6UFmVlWndRtMXN+SlbaTuo3eC/l5BAbMwYBXNDP8PpsRjNEbidDvKKL6IjxgRAMbtUQAiiKQOC/5nAncBmoAp4Vkp5WAjxQyHElf7D/gwkCCHqgK8DE1NT1wJtUsqGCdvCgc1CiENAGdAOPHLSV6Mokxw5coT0dBvb7NlkOZ301qVxm7iPaHzsENC9fh2ZEZnU5Zi5Ofd8HJ4WzObsGZ3ILlAGA9x1l7YozLjnG/zFv33V69Xg2stwRz/l+sX43nvo+IUcw6OPaq/nL6tkhamJwjatz6V2Tj7hYZKIzORgXULwSUnu2uXYU6IBKOYw5eUhbtNpJKBxCFLK16SUhVLKfCnlT/zbviul3OR/75BSXiulnCOlXDXxy19KuU1KuXpSeaNSyuVSykVSymIp5deklN5gXpiiALS2tjLkraJtPIF3jpyDTvo4ugzNYxnZ5Pc1Y9BHIZPCSVyQTVfXJlKSJyfRhc5ttwniYsdx2D7Hny0p+IBzdu7HaLVQOV5HYVw2B9rGoTewBWNcLnjySe398oKdDHnHWDE+hgcYXL4ao96IznTqjFCezGAMxyc9OHO1LKhi2jh4MMSNOo2okcrKactmsxEVZePPOyKYG+aipySRK9hEjhyiETi0/jKKdLG0xen43NwMPJ4REAKDIfRrABwVGQlfvmMI0NMR/n94DTD6fKS+56B99E1693YwlrSSgXcfAY9ryvJefBH6+2FOzhCLrDt4q9RDGFAaHk50pIVES8pMX9JJScrJo6epEW92DgDFjFJWZsfnC227ThcqICinrfLycmR4I+6acbbHnUVXUybf4psA/DbcyKJ4E8KYQN/CRBZnnk9X10ukJF85Ramz7+tfTyLSPMZw33/zsFnLwb/g5Z2Yk9z02AfJSuphj1iGd+cDU5QE99+vvV6ysox8RwtJ1dp40JKMTOKjB8iwntq5Hcn5BXQ3HGEsqhDQMo3GxqpoaPjo85TAqICgnJa8Xi9udw9/KYmiwOtjdK+HC9nKKmroBV5dez6r3dAV7uOCFBdC6Bl3tBIRkRPqpn9IYqLgC7ftBYxsC/sWdUDqyAiJnWHUunZQ/+IYi1YtYE93ODTtPG45u3fDnj0QG+NlTcEWBkUKF3R3AlC9cjUJMcNYc0P/MP2jRETHMD5kJzZtEYMGHTFAOjsoK5vyVCUAKiAop6Xq6mq6fe3IdhdHMhJpLF3Kt/ghAL8BshctAHMGzUtyuCS1kJ6eV0lOmv15iwJ1z73xREfZGR76Lx4O02b8nP/8IUYjevHYdXS3P0rU3HOpeW8z2NuOWcavfqW9bjinnCzdTp5t8ZAlJZ1C0L34fCyRHiLST80xCJMtufZyOmPMABSzRz1HCBIVEJTTUu/hFnbvSSUnOpuCxm6Suvq4mHcYAp4oyOc6t4eucA9nZw0TH3cWw8OHiY5eGOpmH1da2kI+e/1LgJHH9d9mDFjd2kG2HKFRHKHtuUSMprfoTjqfjtd/De7xfzu/tBSefx6MRsmF8zbjGoKi3fsB2JaeQWa/jZgoCzr9qf+VYDRHYIiNxPZ+plGFCghBcur/11eUaeo52MKujiEszn5aRAPl1bl8i58C8AcgfO1izJZCjuRncZHFSf/AuyQkXhjaRk9BCMEX/6OF5KROuh1f5ym9dpeQ/FQdDe4yTPoc6p9+HWtmE9XRa+l48fvgcb5//r3+RPAN5x5kbvx2XrWs5dIO7U7ivbUXE6/rIz46c3K1p6Sk3Hx6GusYS9PuZhbRprqMgkQFBOW04nN6KH+3np7xHrrzDVxRs4m+fVl8ihcZA35nNnNjUjL1pjHOnushJeVybLbdxMWunrLsUJu78HN85fOPINDzY+//4gLWN3ZQPNJBi66VUec6Wl56HNJqaYpfy4Env4/PMczzz8Mbb0BUlJeNix6mxm3BdqCEORJ6dIKqCz9FWvII6bmhH50diOS8OXQ31OFKXgzAUkZpbR1mILhrB52RVEBQTiv2nW086eggljCSmuup8hZw9+h9ADwAmM5ewJyw+TRnZnBBdCvDw+VYretOiYFoU4mIyOacs2qZu6qMZm7jLyINHZDyUgMNw3uJkLkkOxdi2rqNRtfj1OXO4dEHfsFtt2kLAF1+9vPkFPSwy1jMuu1vAPBG/lxy2keISxghLvfUzjA6yhQZhXNsFE9mJj60OfmNHFTdRkGgAoJy2vA5PDxW3kzRiI2W6HHWDb5Cz84CLuEtBoH7gCuXLqE6ws3FhWaiLQuw2w8QdwpMUxGolOzlXHP1ZubFV/Nj+RgOYF33IDnVFXTr63lrdDFzvGmc1zmX1JaXeeCpW7DbI1iUU83PLr6P34edR3ufnStHtO6kzTfcRupgDZFRZgxhptBe3DQtu/Y8Wk0GwoAi3lIBIQhUQFBOG7U7WznS14Hb3MXq1hJ6CsxcX/1PAH4BJM9bQH5MISPJVpYlVOH1DpOack1oGz1NBQs+z4r0Copu3oHTvJQH0Lq6lm9rYKC1hBRh4MHRxXi6RvjZz3/EoYP5JFn6+ON1P+buJffTpUviwiceJAooj4ygM38laREerP4ZRD8uYpJTiUzKoz1Oe7C8lF0qIASBCgjKacHp9vKHPXVcYNtDc2QU7swm3H/OYhX76ETHb4Bblp1LRaSDDZkmTOGJjI+3YrHMD3XTp8VojCUxKYLRdBN/uPgr3G96mnYMrPB46P9nOYeqo9n77hKW/uoHvFaxkOgoNz/4n508ctMP6Rlyk1H6Rz7fp2Ug/eOam5jTNkBs6jgZcy8J8ZVNT8a8Ytprq7Elaw+Wl1CpHiwHgQoIymnhN/84xLlNNWydn0zU/2/vzMOjKPI+/qmZyUwSkkxCSEjIHc4QhBC5EYIgCHiLB3ji8Xqs7qvrq666ruiuJ3uhqy6rooKKgoIcgiAI4Q4SINw3BAgkhNz3JDNT7x/TE4eYwABJhpD6PE8/0131q5pvVdf0r7u6pkqmEu1/imszHAvS/xE7/sFhGBND6RQUTkzEYWqsxURFPeBh1RdG154PcnPARtIGd+PfQ17nDf9nAXi5bD9z5+xkyc+DKSnwJ7h9HmP+N5WfE4LI3HmQXjkfEP3ZesKAoyYvVo26i/DyrfiHlBAYcmlOed0Q5vZhFOfmUB3WDYDenGLnTjsWyzkSKs7KpTuLlULhJkt35RC4YTcnYvdTEFDO9fpMYl4y0E7msx4fvqSSF/reyBHvcm6KCkCvz8HbFIvJ1DIXXArtMJCOIe/xWekIJnZaQEnQOJYt68zIwgN8zas8GgrhY/sTkJjG4ah8QqxWwqsz2PyfI/xYYQVg5i33EFYk6Si8CGpnRqczerhU54dzEEBll7awBHphxWY9zvbtMfRtOa+ELjnUE4KiRbPjeCHr/jOHQL9s1nRtw+8q1+L3i52Eo/uxoeNJKmlrbktU7yDah5gJC9uHzVZJhw53elr6RRHT+RqS2hxjecchRFsyOD3iZrK9vekLTMp9leDNb9E/O5SbD4QTt7+Sw28W8PddpxwL4bQLZNnoibSTe/EJMRHROdnTxbkgfM2BhKVEk6vXEQDE8RObNnlaVctGPSEoWizHs04zdfIMrgyoYE3vCp7Jn091vo6+Hx9Dh+QtQthKLo8NvJETbSvpF9YByCQy4iF0upbd9Dt1v58b9t7NX0vHM6Pfy1StDmTtTQ8w6vtPuKm6hsAdq/n9zjUUBQViKipmmt1Of6DE6MU/H3mfCEsR8Tmn8EsuJCz8fzxdnAsiolsihgMWMv39CS0qpjcrSE9vmWW5VHDrCUEIMVoIsU8IcVAI8UI98SYhxCwtfqMQIlYLjxVCVAohMrRtqkuaK4UQO7Q074mWMBBccUkgpWTlj7N57p1PMLXbQnWvPK4+tQrvvd7EzS7Gr6KCXUTyKrm0DQik8xWRlIfbiAk9iL9fIv7+iZ4uwkWj0xnomjiSkKBiVptGInuW4lVcw7y7HyLPvw0pwHYp2VRQyAG7naFAmbcvH0x4k7SebelceoxE/04ERprw8grydHEuiHbRMZTkV3K6bSQASWxTTwgXyTkdghBCD3wAjMHxH5AJQojudcweAgqllJ2AfwHvuMQdklImadtjLuH/AR4BOmvb6AsvhqI1UF5Tzrzd3zH1jxP4JH0vE/vOpU+ffHTVuSQWFVCdYyF2WzYWjNyPpBq4bsRVYK6gY0gBBi8T0dETPV2MRqNj14ncE7qez20pJETs5FhcBwzler566T6WDE2hxqCnPWDV68jqFMWzD7zEquEhjN+/m5ByHfa408TGXtpTdpwNnU6PXdopinJMudGHo+zeLSkv97CwFow7Twj9gINSysNSymocayPXXVLqJmC6tv8dMOJsd/xCiHAgQEq5QVt7eQZw83mrV1z2WGwWlh7+mecXv8Snr95J9dTpbAntwj0xKygxRJO/fyzJW4+w5kR3+n27F4CnuIHNnCCiXVsSE/pzNHYbQyPC6NTxeRz3N5cHOp2BYf3GYQ4rZU7pQ1zT8Xty/IPxOWwi44Fkbvj7XD666zl+uf9vfH7Vffj0rMA3O58Ei5HQ4HjahGUSEtKyhpvWJSg8glOJoQD0pQK7vZwtWzwsqgXjTkdqBHDc5TgLqLsKd62NlNIqhCgGgrW4OCHEVqAEeFlKuUazd52jN0sLU7QCpJQUFxeTl5dHZWUlOp0Oo6+R07rTHKk4wlF+mqsAACAASURBVJHCk+w9dYCakkxCqsvolwU3HNKxclAnCgMrGWSbx6Ly/rS3BxKYv4/PjNfxj//8Fb3Nzrf+o/lv6Y8A3HjLYEK9ctBFxhATeW+LHVV0NsLCrmV8l8lMz45kiL4LfbutZlH2cOL2FzAuYD4/9bmOuYFtifE6TZesSu6s8GFPBzMhHU8QHZPc4t+lxPbsTdD6ueTqdITa7XRiMZs23cGQIZ5W1jJxpzXUd6cv3bTJBqKllPlCiCuBeUKIRDfzdGQsxCM4upaIjm4Zc7UrHBf99YfySTucD9iprDhGWX4m1fmlhHpVEhJgok1AOKdlEQX2Yo4V6iiynSRAHKadvYYxxmjCq0dx+kQ+2V2qOdUjhzblVVTvjWZD7CD8ZBGbCzsTVZXLWx/9Ce+qGvZEXMG9JyRQQc/uCSSGXs2+7ht5tPMEzObenq6SJmN88oMsK17EBxtu5dmgD7ix/RLmFQ0iMMaLQQEb6VawFmumgfncwjXtu3Gomw8dvTYRE/2up6VfNH5tg4kIjmRfmwBCS4vozw9s2nSHp2W1WNxxCFmA67y4kcDJBmyyhBAGwAwUaN1BFgAp5WYhxCGgi2YfeY480dJ9BHwE0KdPn3qdhuLSYueJYr7bnMWgOCPJ1r+zf6cdyn3wrS7C21qFn6093jIcncwnHDtBAkL1ldRIX+wVEeirc7HGHuVo8DFOx4ewo6grHY770ll3gp86dsZ02EZyXgeS8hdzw48LCCitoCAsjgk1o7HwN/R6P+4dP5IgfSYiOoa4DnV7OC8vjMZ23JVwBUuL1/LF4f/j6QHpPLp1Lr8cvoLImDWUlngxjad4o017lif406fqG5L6PIFO5+Vp6Y2CX1A0x0LDoLSI/vzCe+meVtRyccchbAI6CyHigBPAeOCuOjYLgPuBDcBtwAoppRRChOBwDDYhRDyOl8eHpZQFQohSIcQAYCNwH/DvximSwpN8tfEo+ZknGXh4MYe2nORIfglHt+0gPw/KLDZqsOLvL4kMhsQYO11ivfD2DiLXtw9Wfy8C7RZEdTClp3vicySUYEshg22H2REYwjp9KDFb08gNisBWcITxc2djrKmhMDKB53oPZ9vCKQDcfcdI4uxxrOyxlb/3+a+Ha6R5GBOdzM+FksjKBUze0IvuVxlpW7Ofn9PG4B1YzP92WkNpcDyB5fsZ0PdhAgJa/kgrJxEJPVjXMRwO7aU/mRw8CPn5EBx87rSKMzmnQ9DeCTwJLAX0wKdSyl1CiL8A6VLKBcA04AshxEGgAIfTABgK/EUIYQVswGNSSues5Y8DnwM+wI/apmihSCl5c9ZGAvceALmRbw8fYcXyjRzJya+18QECcDSEzcBXgJevNxH9gugxNI9ObdpzqiaWUn0Ap6PbUuxlxssWgZ+tO6F2iZ9VR49dGdzwyUzaZjlWVU9PGsrUvgnMmPYVUENk2DjuHNCbjfrdPDz4D3gbWtYMnhfD77v0YLZvCGN2prJxbRIHjX3pEOzNaLMf+b4lbDFWM6zfU4SZ/T0ttVGJTEikomMEdiAJCyZOsWFDe66/3tPKWh5uvVGSUi4GFtcJe8Vlvwq4vZ50c4A5DeSZDvQ4H7GKS5PS0lL+/O5SIi0nKZUrmDE7gz0HjhIPvKLTM9boTTebFXPNmRPNlAo4XlHFsdTNHE3dzEmTDyUJXdBd3Z0UXxNdS0rwtldhL7Fhziwgfvch2uY77icq/b35buQ4lgX6MPeLr7HaSxDiWt54ugtFthp8kvuSFJrkgdrwHDE+JuL8/TAOuYnXbg7AeqoCu8WGsYMfeywWdAWldL/MnAGA3uBFaLCRQ3oDnW1WevEt69Y9qRzCBdCyhxgoPM7u3buZMXsHAfpcTucv5P3pa0moqOInoWektIHdBlWOgeFWvZ4aXy+EXaKvsuFvs9Idx59bALBUQsY2x9YA1W3asG9oFLO7XM/q7elsmL+WGpsVGMm44Y/RxniEeZFlzBj8cFMX/ZLktrC2fHz8NPNyi7i5fSBCCLYUl7Mor5iX4sM9La/JSO43kAMB8+hcWMAAlrJu3ZOeltQiUQ5BccGkpaWxPDUHqy6Xwh3zmTVnJf+y2x1DwqSNCuHDD8FjONkxkHZXZ2PoaeHY3ms4XdkJX0sJMfn78LOZ8bFU4F+cje3kdizZe2lXXE47HMPOanSCPLORgna+7OgYwEbvWE4etbBv1YeUVVRoSh4nsM1rDBz+PUuiDfxj1HMtYgW0puJ/okJYW1jK20dy0AExPkb+FB+O7jKuk+SRd/Hf8HegsIAhbOG/v4DFAiaTp5W1LJRDUFwQGzdu5Id1uRRWHMey8yc2f/8zW4B4wIKRKTzNZPk8BXnBkAdshNDgbG5L/oYxA6fg096MvmsIpflV5GVXUmYOIrPfE9jwo6S0nO3bVrN5zyLy8nOh0OLYDhQCR2s1RLZrS17pR1RZxnHD3am0HZHM+IAoQv1bz3uDhrgqyJ+rgi6/7qGG8PL2JatLAuw+yFBysFgkW7YIBg70tLKWhXIIivNm8/YdfLI2G1NxFhVb0slbtIg1gB+QYezJUyl/ZZRpJa8YXmKXX1fS06/i4OFu5OaH8+GyP/DlnnvoefMaRnovJTbYQmR3A9iPYazcgMkm0Qf70qeLL3fYRnDwcDl7dheQk3Maa1Upft7+xARGE9zNzoIVH1KV14WUoVYefbYHB3JLGZHQ3tPVo/AQomc02fMgHDvdWMXatcOUQzhPlENQnBdHT5zkteV7SSwtIHtrFjWLvmIeEj3wfcJoPpl4A29veJljA+LJ9Alj+LpUkod/i3+inm3545iZcTcnskJZ+/6tZA/uSdKodfQozCIuoBdRtghs5ZJTORlUVh0mXFRgIwR7tyu4p7OZEK9gTkgL1dKHhXv6cDCzC6GhkrtePMWB3GruHxTr6epReJDBA+LY6utNeEUVKXzJunXDeO45T6tqWSiHoHCb8soqnvwmlWFVNvbsLKXyhw+ZiQ098PGY29ncqxvPbvmajOt6oi/Qcd3qU3R99z/owkJZuXMauct288cO97M2czxzfxnPoXWdOLW1Pf43z8eUvI0s/2xi2tqJLLFiscZz0suI3urFdYcqWWE0s9jLn2E+1VSUpTBnZRxCJxn9xDF6dfGjf7ya+aS1M3zgbfwnfApjD2WRwip+vxbsdtCpVV/cRjj+TNwy6NOnj0xPV39D9BT3ffAtsVYbeVtLODr9TeZzFAPw2dhb2Nu1I33tR8hNhrh9lXS39iD6zTcR+l8nk7Pb7XyQNp/NG47SNyOHbzbcy9pDjj9IhUUdYcSE9wiJ20iV3UTnU5CyKQ/vSlgb3Y+lMf2Z+tTtTHrdl6lTHG8K3/hnFS/9Qb0vUPzKI7ffxkffzeEEBiKpZutWQVLrGn1cL0KIzVLKPueyU08ICrd488cNVIkKdEcsbJy5glTNGcwePpbd3eKI9aqgoIeFnunVJCTfQ/B99/4mD51Ox+8H3cKRxCImr0mla79Uem/aytyF13LieBxfTf4XfcO3MzFxNld334bx/57g7aNGMNp5KOh6rh9jYNMmxx3fhx/Co48qZ6A4E1v/K8j7bg4RWOnIZpYv76McwnmgHILinOzMKWT9zp1MKBL87Ssb82pm4Q+s7pLE0iED6GItwrv7fnqledHlumcIuPbas+YXZw7kg+tuYknWcZbFb2XMwJ85tLw7639MYFN2TzZl90S30k6b2eW0C9CTl+3LJ6WOtJGRMG0ajBrV9OVWtDwGR8EvPkbGVlYzgs9YvrwPzz7raVUtB9VlpDgrdrudW/4xkz5tykn7MoI/bhjHUKrZ5d+ON194nk7l+cQkbuTK9T50nPAn/AYPPq/8pZTsLask9fghTucUkb6gHTtWRHB8tz8226/j5jt1gocfhscfh4CAxi6l4nLh9IkM3ht1M3/dfZR5dOIun/0UFopW/38E1WWkaBRe+uEXvNsXoF/Rln4b/shQqskRBv7+h6cJryqlQ48NDNjgQ+yDb+CbfP5TTAshSPD3JaH7FY6/LA93hJeVQWYmVFdDhw4QFtaoxVJcpoR06MW2/sNg93RGcBhrZQ1paUZSUjytrGWg3r8rGmR/fhmZBzZw204TG+Zs4mV2Ywf+PeFefHUWorusJ2WTgdiJb16QMzgbfn7QowckJytnoDgPhCCwRzy7AH/sDGYpy5d7WlTLQTkERYO8PHMxocF6Zq304sOyf6MHvujam+M94oiL3co128qJGv8WvsnJnpaqUNQyItbOikAfAMYwjSVLWk63uKdRDkFRL19mHKONaT/dlgUxIf15opBsNvqz8P7b6Ra8i1F7jxNz6z9p06+fp6UqFGdw85AJpCd2AmA0q0lPF5ysd/ktRV2UQ1D8hqpqKwt+WsigLD8OzJ3JOPIpRvDuE78jwT+TkSe2Ez/6XdqoeQEUlyDmkK6c6j+KMqAnhURxlIULPa2qZeCWQxBCjBZC7BNCHBRCvFBPvEkIMUuL3yiEiNXCRwohNgshdmifw13SpGp5ZmhbaGMVSnFxPD9vA+1CC9j4fTFvVi0B4I3hNxIcaWVUaSrdBv0NP/WWTnEJ06OTDz8ZHKPUbmUq8+fbPayoZXBOhyCE0AMfAGNwjAOZIIToXsfsIaBQStkJ+BfwjhaeB9wgpbwCxxKbX9RJd7eUMknbci+iHIpGYm9uKWUnUun9QwBP7XwLH2BmUCyFo3swUvxAry5/xv/a6zwtU6E4K3cPSGJZrGOiw9v4hhU/S8rKPCyqBeDOE0I/4KCU8rCUshr4Bqi7avlNwHRt/ztghBBCSCm3SimdvXe7AG8hRCsfEXzpIqXktVkLSSyvggVTuQIL+4UXi5+8g5Ft5jHI/AgB4+7xtEyF4pwk9byZHYOuphK4ikyCq3P46SdPq7r0ccchRADHXY6ztLB6baSUVqAYqLvE9Thgq5TSdR3Fz7Tuoj+L1ryiySXCrK3HaK/fTMGnJ/ifmr1YgDdvuZtropYzvGosgQ8+42mJCoVbCL2e7omx/KhdVW7lc2bNKvWsqBaAOw6hvgt13XFcZ7URQiTi6EZ61CX+bq0raYi2/XbyG0faR4QQ6UKI9NOnT7shV3EhVNXYWPjzHBIWSZ45MgOAt6P7Mui6Q4w6kUi7ZyZ7WKFCcX5MGBDHkqi2ANzDJyyYb1TdRufAHYeQBUS5HEcCdQdx1doIIQyAGSjQjiOB74H7pJSHnAmklCe0z1JgJo6uqd8gpfxIStlHStknJCTEnTIpLoDX56+lV9kGEpZMJxjJj/oAgv8cwIhdPnR4dYan5SkU582wQRPZcdUIioH+ZBJvOcicORXnTNeaccchbAI6CyHihBBGYDywoI7NAhwvjQFuA1ZIKaUQIhBYBLwopVznNBZCGIQQ7bR9L+B6YOfFFUVxoWTmlVFwdC6BH2Yx1F5ANoK0P/Tjqq2lxP99iaflKRQXhDB4kdQzkm+0/osHeJ//fnDMs6Iucc7pELR3Ak8CS4E9wGwp5S4hxF+EEDdqZtOAYCHEQeAZwDk09UmgE/DnOsNLTcBSIcR2IAM4AXzcmAVTuM9r38ym57xDPJyXhh14p89QRpacotd761v1YvWKls/j1w7k+07hANzLl2xOj+Xg3nwPq7p0UbOdtnIWZRwj47M7efC9dMKxMtmnI4kPeDP2X+kIo1pvQNHCkZKrHnqAqZ9Npwcwga/wHR/AtK+v97SyZsXd2U7VP5VbMWUWK98ue4Wh/z1MOFZSMWEcH87Yv69XzkBxeSAEdw4N4wOTY2LnZ3ideUsHcnCLurGsD+UQWjEvffMx1324mSGWXE4Ds4aO5tHJMxA+asEBxeXDo3c8R+rAPpwG+rKHxMJdzPx6M1VqyNFvUA6hlbJmXyaJM97n9syd2IDn2o/llU8ew6ddnKelKRSNitE3mJGDuvOhdvwCr/HVghGkza87cYJCOYRWSGW1lUVTxnHPij0AvKjvyROzxhDeebSHlSkUTcNfnv4js3t0ohQYywpC92ezNxt2r17haWmXFMohtEJefPsOnvhoB22QfE4QUW/dS5+hT3halkLRZASGdKHf1cn8Qzt+i+d599OrKS05xKnDBz2q7VJCOYRWxj9nTOJ3bywhyl7DerxYd9vLPP6HR9XwUsVlz7uT3uDz+AjygKtIo8eeHaz55RS71qygorjI0/IuCZRDaEUsWPsdo574F12qK9mFjkmJ/+bd6bdhMPh7WppC0eQEBHdiwm0pvKwdv8eTTJ02nqCIctZ/O5PqqkqP6rsUUA6hlbBs9Vyir3uIHmWlHAXubzuFz5f2x9c32tPSFIpm4803PuenhE6sA8LJ5YWTbzNlegjdh/Zizczp2Kw1npboUZRDaAXMXfgV4TfcT1JJCVnALb7v8PmqZCIikjwtTaFoVoTBi+n/eJnHDHqqgIeZhm2WnR8WzKDXqJGs/vIzrDWt1ykoh3AZI6Xk4z//maQ7HqJHSRmHgWu83mXKj8Po0WOwp+UpFB5hyJj76XvrMJ7SjqfWPM7X79/GmiUv0fva61j1xbRW232kHMJlSlZBIQvGjeX2N98kvsrCdnSMNH7CvxelMHRovRPLKhStho9nLmVZRAhfAn5UMPPofXz62QRWLXie/jeOY/WXn1Gan+dpmc2OwdMCFI2LxWZjwZR38Z45lZu2HADge7x5xHsRc5Z2ZejQumsbKRStD71ez7pfMriiSzzh5RZGkMOH63/Pa4EvYbM/w6gb32bLkoXEJV1JVGJPT8ttNtTkdpcJNrudFfMXUbxwCsmz1hNfUUUN8Aodmd5+Pkt+iqVnzzaelqlQXFKkpaVx3ZDBzLfauQooJJDJKU/R55499I58hAqTnbKCfK684Ra8jC139V93J7dTDqGFU1NdQ+rXX1ORPpO28zIYknUKgN3AvTyB96BXmDcvFLW2kEJRP6tXr+bGEVfzqdXOrVrY7MibMTxrg4Jkhj/8IFtXLCW6Ry/ik/u2yP/sKIdwmXN01x4OfPM+uSd20GnxbvqdcszxXgK8Rkc+1H3HK6/34PnnDej1ntWqUFzqpKenc/01wxlfXMo7OBZsKRb+bEjpT+XI9vjaO9Nx1HCO795DfHI/oq/o1aIcQ6M6BCHEaOBdQA98IqV8u068CZgBXAnkA3dKKTO1uBeBhwAb8L9SyqXu5FkfrdkhSKuVI5vSOT7/I4pz9iF+OcmQ/ccJtNkAKAY+IJQp/I3e197Fv98z0KWLZzUrFC2J06dP8+Cdd3BkZSr/AK7Vwi3CwP6ErhQlh5Pp1412A3yRVQZCIrvR/arraWMO8qRst3DXISClPOuG44J9CIgHjMA2oHsdm98BU7X98cAsbb+7Zm8C4rR89O7kWd925ZVXytaA3W6XFSdOyv2f/VduuGuYnD+qu1ycECl3+beRNSCly7YF5FMkSj0Jcvhwm1y2zJHHpEmTZEpKSm2ekyZNklLK2jDX+JiYGJmSkvKbNDExMbXp6to4NymlNJvNMiUlpdbGbDbLSZMm1YabzWap1+ulXq+vTR8TEyNjYmKk2WyWZrNZmkymM/JNSUmRJpPpNxomTZok9Xq9dDRdWZuXM19nPmazWcbExEiTyXSGPlc7IURtvMlkqs3XuTmPneUCpNlsloAUQtTqB2rjnZszb71eL4UQZ2zONK7f4/wuk8lU++msI71eX1sWpx6nZmf9udaFM42zvp156fV6aTabpRDijDpKSUmRer1emkym2jydup1hdeumqTbX7xFC1O47z6uzDlzPtbN8zjpxnuNJkybV7jvbpGv7drYzZ/ty2v+wcKHsGhslh4OcD9Lm8nuzgTwWECS3JfSQqYNS5Ndjhsj3Hxgp337ldvnu7Mfkgh3/lJuOLpRHC/fL8upyabfbL+ga4NTUWADp8hzXVynluZ8QhBADgVellNdqxy9qjuQtF5ulms0GIYQByAFC0JbSdNo67bRkZ82zPi70CeHA6p3YaqxIu13bpLbZsdvsCJz7jk+kdLH7NQ12O0iHHUikzV5rj5TYba5p7SBB2u0ImxWdpQpdZTlUlmIrLcZWlE91QS7VBblQWoy+ugbv6moCqi0E1lQTZbVg4rfnpgZYj2AJHfmB28kNf5yJ90YxebLA9Vw6H2edYUI44l0/nfF1H31d0zRk42rbHI/ODelsSY/tTYmqC/ep276dYa7tXUrJypUrefedKexc/hM32S3cCAzCcQdbH5XoKNAZKNQbKNIbqDQYqDLosRj01OgFNQY91To9NQYj0qBHGAzgZUB6GRBeXgijAZ2XEWEy8sjiH/nszjvQ+3hjMPng62PmiltuJX5o/wsqs7tPCO4MO40AjrscZwF1VdXaSCmtQohiIFgLT6uT1jnu8Vx5NhoRKVfg21SZNyGngEMIDuLPZmLYzEBOtL+Bnv2HMXy4H1+kQM+eoNPB5MmeVqtQXD4IIRg+fDjDhw+npqaGTZu2Mv3zVB5buoqQ7H30rjlBN6qIx9HNEQP4YCfCXk2Evdpx53YRPAJMnDX7jLBP1x4jftPXF5fxOXDHIdR321H31rUhm4bC6/tDXL2PKkKIR3DUD9HRFzbvTjZ6fJANPKcKt8LsvwkTDdpTJ8yKoBIdVeiwCB2VwguLMDo2nYkygx8lXgFU+IRg8Y/AFhiFpX03wjt3o1u3aDp00PNgFLzdBXx8fi3XsGHDWLVqlWtd1Vd/v9mvL6yhNA3ZuBPXmLijszWj6sJ93P2dTJo0iVdffRUvLy8GDerHoEH9gOcBqK6G3bsr2br1ND/vzuXkiTyq8gqQeafxKs7FtzIPr+pSDDWVmKxVGG3VeNst+EoLXtKKHisGaUePDYO0YcDGApnP95T8qkP7vAcfJuKDrUPXJqsTJ+44hCwgyuU4EjjZgE2W1mVkBgrOkfZceQIgpfwI+AgcXUZu6P0NHaX1QpJd8qSmptbuO7uCXI9BdRm1FlRduI87XUbnwmiEpCQfkpKigcaZIHKiy37d3zPAiEb5lrPjztQVm4DOQog4IYQRx0vjBXVsFgD3a/u3ASu0FxkLgPFCCJMQIg7oDPziZp4KhUKhaEbO+YSgvRN4EliKY3TQp1LKXUKIv+B4c70AmAZ8IYQ4iOPJYLyWdpcQYjaO/0lZgSeklDaA+vJs/OK1HlJSUs44njRp0hlPEJMmTTrDzjU+JiaG2NhYhg0bdkaamJgYJk6ceFYbALPZTFJSEpmZmcTGxpKRkcHTTz/NlClTSEpKIiMjgzJtQfOXX36Z1NRUMjMzASgqcixMUlVVxYABA87INy0tjbCwsDM0TJw4kddffx2bNtw2JSWFzMxMJk6cSGpqKmlpaQwYMICMjAwCAwPJycnB29u7Vp+rXXV1NUajEW9vb6qqqrBarbX5gmN6A5vNVlv2VatWYTabKS4uRghBdHQ0RUVFFBcXk5KS8pvuO6PRiNVqxW63n1GugIAAiouLz/geAJvNhslkwmKxYDKZCAsLAyArK4vIyEhycnKwWCzExMRQVFREVVUV3t7eZ7QBZ71mZWXV1vfbb79NWFgYWVlZ+Pn5UVJSwtChQ2vrKDY2lrVr12IwOC4H3t7elJWVYbfbMRodr1Dr1k1T4axzOPMu2WQyMWDAANLS0rBYLKSkpNSea2c7cdbJCy+8QGpqKsOGDQMcT9LONunadp2/Baed8zdyKeApLeqPaQqFQnGZ4+4oIzXbqUKhUCgA5RAUCoVCoaEcgkKhUCgA5RAUCoVCoaEcgkKhUCiAFjbKSAhxGjh6gcnbAZfimnhK1/mhdJ0fStf5cbnqipFSnnNVlBblEC4GIUS6O8Oumhul6/xQus4Ppev8aO26VJeRQqFQKADlEBQKhUKh0ZocwkeeFtAAStf5oXSdH0rX+dGqdbWadwgKhUKhODut6QlBoVAoFGfhsnIIQojbhRC7hBB2IUSfOnEvCiEOCiH2CSGubSB9nBBioxDigBBiljY1d2NrnCWEyNC2TCFERgN2mUKIHZpdk8/oJ4R4VQhxwkXb2AbsRmt1eFAI8UIz6PqbEGKvEGK7EOJ7IURgA3bNUl/nKr821fssLX6jECK2qbS4fGeUEGKlEGKP1v6fqsdmmBCi2OX8vtLUurTvPet5EQ7e0+pruxAiuRk0dXWphwwhRIkQ4uk6Ns1SX0KIT4UQuUKInS5hbYUQy7Tr0DIhRFADae/XbA4IIe6vz+a8cWfh5ZayAQlAVyAV6OMS3h3YBpiAOOAQoK8n/WxgvLY/FXi8ifX+A3ilgbhMoF0z1t2rwLPnsNFrdRePY2nZbUD3JtY1CjBo++8A73iqvtwpP/A7YKq2Px6Y1QznLhxI1vb9gf316BoG/NBc7cnd8wKMBX7EsUDYAGBjM+vT41gDPsYT9QUMBZKBnS5hk4EXtP0X6mvzQFvgsPYZpO0HXayey+oJQUq5R0q5r56om4BvpJQWKeUR4CDQz9VAOJZLGg58pwVNB25uKq3a990BNO0iqY1LP+CglPKwlLIa+AZH3TYZUsqfpKxd8i4Nx+p6nsKd8t+Eo+2Aoy2N0M51kyGlzJZSbtH2S4E9/Lp2+aXOTcAM6SANCBRChDfj948ADkkpL/QPrxeFlHI1jjVkXHFtQw1dh64FlkkpC6SUhcAyYPTF6rmsHMJZiACOuxxn8dsfTDBQ5HLxqc+mMRkCnJJSHmggXgI/CSE2C8e60s3Bk9pj+6cNPKa6U49NyYM47ibroznqy53y19pobakYR9tqFrQuqt7AxnqiBwohtgkhfhRCJDaTpHOdF0+3qfE0fFPmifoCaC+lzAaHswdC67FpknpzZ03lSwohxHIgrJ6oP0kp5zeUrJ6wusOr3LFxCzc1TuDsTweDpZQnhRChwDIhxF7tbuKCOZsu4D/AX3GU+a84urMerJtFPWkvepiaO/UlhPgTjlX3vmogm0avr/qk1hPWZO3ofBFC+AFzgKellCV1orfg6BYp094PzcOxpG1Tc67zQyQENgAAAoVJREFU4sn6MgI3Ai/WE+2p+nKXJqm3FucQpJTXXECyLCDK5TgSOFnHJg/H46pBu7Orz6ZRNAohDMCtwJVnyeOk9pkrhPgeR3fFRV3g3K07IcTHwA/1RLlTj42uS3thdj0wQmodqPXk0ej1VQ/ulN9pk6WdZzO/7RJodIQQXjicwVdSyrl1410dhJRysRDiQyFEOyllk87b48Z5aZI25SZjgC1SylN1IzxVXxqnhBDhUspsrfsstx6bLBzvOZxE4nh3elG0li6jBcB4bQRIHA5P/4urgXahWQncpgXdDzT0xHGxXAPslVJm1RcphGgjhPB37uN4sbqzPtvGok6/7S0NfN8moLNwjMYy4njcXtDEukYDfwRulFJWNGDTXPXlTvkX4Gg74GhLKxpyYo2F9o5iGrBHSvnPBmzCnO8yhBD9cPz285tYlzvnZQFwnzbaaABQ7OwuaQYafEr3RH254NqGGroOLQVGCSGCtO7dUVrYxdHUb9Gbc8NxIcsCLMApYKlL3J9wjBDZB4xxCV8MdND243E4ioPAt4CpiXR+DjxWJ6wDsNhFxzZt24Wj66Sp6+4LYAewXWuQ4XV1acdjcYxiOdRMug7i6CvN0LapdXU1Z33VV37gLzgcFoC31nYOam0pvhnq6Coc3QXbXeppLPCYs50BT2p1sw3Hy/lBzaCr3vNSR5cAPtDqcwcuowObWJsvjgu82SWs2esLh0PKBmq0a9dDON45/Qwc0D7barZ9gE9c0j6otbODwAONoUf9U1mhUCgUQOvpMlIoFArFOVAOQaFQKBSAcggKhUKh0FAOQaFQKBSAcggKhUKh0FAOQaFQKBSAcggKhUKh0FAOQaFQKBQA/D+1QP0eNAcz5wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#Iterations of Collapsed Gibbs sampling\n", "plt.figure()\n", "iters = 100\n", "expected_pdf = 0\n", "for ii in range(iters):\n", " # Given z, posterior predictive distribution is mixture of Student's t-distribution\n", " alphaN, mN, kappaN, aN, bN = GMMParamsPosteriorParams(alpha0, m0, kappa0, a0, b0, x, z)\n", " weights = (alphaN)/alphaN.sum()\n", " pdf = sps.t.pdf(t[:,np.newaxis], loc=mN, df=2*aN, scale=np.sqrt(bN*(kappaN+1)/aN/kappaN)).dot(weights)\n", "\n", " # The PDF that we have evaluated here is the predictive distribution only for a given z. \n", " # To approximate \"full\" posterior predictive distribution, we need to average the PDFs\n", " # obtained for different samples of z\n", " expected_pdf += pdf\n", " if ii % 1 == 0:\n", " plt.plot(t, pdf, lw=0.5)\n", " \n", " # In the inner GS iterations, we ...\n", " for n in range(len(x)):\n", " # remove one observation x_i and corresponding z_i\n", " z[n] = 0\n", " \n", " # calculate P(z_i| x_i, x_\\i, z_\\i), given all remaining observation x_\\i and z_\\i\n", " alphaN, mN, kappaN, aN, bN = GMMParamsPosteriorParams(alpha0, m0, kappa0, a0, b0, x, z) #2)\n", " weights = (alphaN)/alphaN.sum()\n", " p_zi = sps.t.pdf(x[n], loc=mN, df=2*aN, scale=np.sqrt(bN*(kappaN+1)/aN/kappaN)) * weights\n", " \n", " # resample value of z_i using this probability\n", " z[n] = sps.multinomial.rvs(1, p_zi/p_zi.sum()) #3)\n", " \n", "#plot the true GMM, ML EM trained GMM and empiricl posterior predictive GMM from GS\n", "plt.plot(t, true_GMM_pdf, 'b', lw=2);\n", "plt.plot(t, EM_GMM_pdf, 'k', lw=2);\n", "plt.plot(t, expected_pdf/iters, 'r', lw=2);\n", "plt.plot(x, np.zeros_like(x), '+k');" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "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 }