H1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import norm"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Homework 1\n",
    "## Problem 1\n",
    "1. Using neural activity to predict monkey movements\n",
    "\n",
    "a. Imagine that you are recording a neuron inside the monkey motor cortex. You suspect\n",
    "that activity of this neuron reflects whether the monkeys hand is open or closed. Every\n",
    "20 ms, the neuron either fires an action potential (called a spike) or it does not.\n",
    "\n",
    "Let us assume the following probabilities:\n",
    "\n",
    "$p\\left( {{\\text{hand open}}} \\right) = 0.5$\n",
    "<br />\n",
    "$p\\left( {{\\text{spike}}|{\\text{hand open}}} \\right) = 0.6$\n",
    "<br />\n",
    "$p\\left( {{\\text{spike}}|{\\text{hand close}}} \\right) = 0.35$\n",
    "<br />\n",
    "\n",
    "You observe a series of 10 spikes spaced evenly over the course of 200 ms.\n",
    "What is the probability that the hand is open, given the data?\n",
    "\n",
    "*Answer*\n",
    "\n",
    "The probability that the hand is open is\n",
    "\n",
    "$\n",
    "p(\\textrm{open} | \\textrm{10 spikes}) = \\frac{\n",
    "p(\\textrm{10 spikes} | \\textrm{open} )p(\\textrm{open})\n",
    "}{\n",
    "p(\\textrm{10 spikes})\n",
    "}\n",
    "$\n",
    "\n",
    "and also\n",
    "\n",
    "$\n",
    "p(\\textrm{10 spikes}) =\n",
    "p(\\textrm{10 spikes} | \\textrm{open} )p(\\textrm{open}) +\n",
    "p(\\textrm{10 spikes} | \\textrm{closed} )p(\\textrm{closed})\n",
    "$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The probability of the hand being open, given 10 spikes is 0.9955\n"
     ]
    }
   ],
   "source": [
    "pOpen = 0.5\n",
    "pClosed = 1 - pOpen\n",
    "\n",
    "p1Open = 0.6\n",
    "p1Closed = 0.35\n",
    "\n",
    "p10Open = p1Open ** 10\n",
    "p10Closed = p1Closed ** 10\n",
    "\n",
    "p10 = p10Open * pOpen + p10Closed * pClosed\n",
    "\n",
    "pOpen10 = p10Open*pOpen/p10\n",
    "\n",
    "print(f'The probability of the hand being open, given 10 spikes is {pOpen10:.4f}')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "b. What is the probability for (a) if $p(\\textrm{open}) = 0.05$?\n",
    "\n",
    "*Answer*"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The probability of the hand being open, given 10 spikes\n",
      "and a prior proabability of open of 0.05\n",
      "is 0.9202\n"
     ]
    }
   ],
   "source": [
    "pOpen = 0.05\n",
    "pClosed = 1 - pOpen\n",
    "p10 = p10Open * pOpen + p10Closed * pClosed\n",
    "pOpen10 = p10Open*pOpen/p10\n",
    "print(f'The probability of the hand being open, given 10 spikes')\n",
    "print(f'and a prior proabability of open of 0.05')\n",
    "print(f'is {pOpen10:.4f}')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "c. What would be the probability if you had observed only 5 spikes\n",
    "in the same 200 ms?\n",
    "\n",
    "*Answer*"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The probability of the hand being open, given 5 spikes\n",
      "is 0.5665\n"
     ]
    }
   ],
   "source": [
    "pOpen = 0.5\n",
    "pClosed = 1 - pOpen\n",
    "p5Open = p1Open**5*(1-p1Open)**5\n",
    "p5Closed = p1Closed**5*(1-p1Closed)**5\n",
    "p5 = p5Open*pOpen + p5Closed*pClosed\n",
    "pOpen5 = p5Open*pOpen / p5\n",
    "\n",
    "print(f'The probability of the hand being open, given 5 spikes')\n",
    "print(f'is {pOpen5:.4f}')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Problem 2\n",
    "\n",
    "In the lecture we saw equations that allow us to calculate the\n",
    "posterior distribution of the mean if the prior distribution is\n",
    "normal and the likelihood is normal as well. It is also possible\n",
    "to use a method called grid estimation. To do grid estimation, we\n",
    "select points  $\\mu_k$ for $k=1,\\ldots,K$ . For each point $\\mu_k$ ,\n",
    "we calculate a value proportional to the posterior distribution\n",
    "through simple multiplication using the formula:\n",
    "\n",
    "\\begin{equation}\n",
    "p( \\mu_k |  x_i, \\sigma_x ) \\propto p( x_i| \\mu_k, \\sigma_x)p( \\mu_k)\n",
    "\\end{equation}\n",
    "\n",
    "Let’s call this multiplication  $q\\left( {{\\mu }_{k}} \\right)=p\\left( \\left\\{ {{x}_{i}} \\right\\}\\left| {{\\mu }_{k}},{{\\sigma }_{x}} \\right. \\right)p\\left( {{\\mu }_{k}} \\right)$.\n",
    "We can normalize to approximate the posterior by dividing by the total:\n",
    "\n",
    "\\begin{equation}\n",
    "p\\left( {{\\mu }_{k}}\\left| \\left\\{ {{x}_{i}} \\right\\},{{\\sigma }_{x}} \\right. \\right)=\\frac{q\\left( {{\\mu }_{k}} \\right)}{\\sum\\limits_{k}{q\\left( {{\\mu }_{k}} \\right)}}\n",
    "\\end{equation}\n",
    "\n",
    "Usually, this approach yields an accurate picture of the posterior distribution with only 100 points and certainly provides accuracy with 500.\n",
    "\n",
    "a.\tUse this approach to solve again the exercise at the end of the lecture, but this time using a grid approximation. All numbers are in Newtons:\n",
    "Your prior about the box is  $\\mu_0=3.5$ with confidence $\\sigma_0 = 1$.\n",
    "<br >\n",
    "Your confidence in your perception is $\\sigma_x = 1$.\n",
    "<br >\n",
    "Your lift the box multiple times and perceive the weight to\n",
    "be 6, 6, 7, 7, 4 and 5.\n",
    "<br >\n",
    "Plot the posterior of distribution of the weight after each lift on a grid of 500 points between 0 and 10.\n",
    "\n",
    "*Answer*\n",
    "\n",
    "Set up the constants and the underlying grid"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean of inputs is 5.83\n"
     ]
    }
   ],
   "source": [
    "mu_0 = 3.5\n",
    "sigma_0 = 1\n",
    "sigma_x = 1\n",
    "x = [6, 6, 7, 7, 4, 5]\n",
    "mu_k = np.arange(start=0, stop=10, step=10/500)\n",
    "\n",
    "print(f'Mean of inputs is {np.mean(x):.2f}')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now we can generate the prior"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [
    {
     "data": {
      "text/plain": "Text(0.5, 1.0, 'Prior')"
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEYCAYAAAB2qXBEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAj10lEQVR4nO3dfXAc9X3H8fehQ8I1IFtCBuck21LXPevOBhufJDJJGCMeRGu6GYwinwOuiARiYrUTnKGjzrQR+YOplMyQQLD+QFNNoykTLTOalnNjfK4jSBOSmOvFD4x9xT6VE0gHJbICKn6SrOP6R8INwpYtfNKeTvq8Zhi8t7+9/f7uYT+7v9XtOpLJZBIREZErdFWmCxARkeymIBERkbQoSEREJC0KEhERSYuCRERE0qIgERGRtChIROaAa6+9lrfeeivTZYhcEQWJyCxZtWoVixYt4tprr+XGG2/k4Ycf5tSpUxdte+rUKcrKymyuUGRmKEhEZtG///u/c+rUKQ4ePEg4HOapp56aNH9iYiKt5093eZGZoCARsYHL5eLP//zPOXr0KA6Hg46ODlavXs3q1asBcDgc9Pf3AzA6Ospf/dVfUVRUxMqVK3nqqaf4+OOPAfjxj3/Ml770JXbu3ElhYSHf/e53M9UlkRRnpgsQWQgGBwd5+eWX2bJlC3v27OGll17i9ddfZ9GiRRe0/Zu/+RtGR0d56623GBkZ4Z577mH58uU0NjYC8Prrr+P3+3n//fc5f/683V0RuYBD19oSmR2rVq3i5MmTOJ1O8vPz2bx5M08//TR/8id/Ql9fH9XV1am2DoeDaDRKaWkpixYt4vDhw3g8HgCef/55enp6+PnPf86Pf/xjWltbeeeddzLVLZEL6IhEZBa99NJL3HXXXRc8XlJSctH2J0+e5Pz586xcuTL12MqVK4nH45ddViRTdI5EJAMcDsdFH7/hhhu4+uqrefvtt1OPvfPOO7hcrssuK5IpChKROSQnJ4e6ujr+/u//no8++oi3336bH/zgBzz00EOZLk1kSgoSkTnmueeeY/HixZSVlfHlL3+Zr3/96zQ0NGS6LJEp6WS7iIikRUckIiKSFgWJiIikRUEiIiJpUZCIiEhaFuQPEm+44QZWrVqV6TJERLLGwMAAJ0+evOi8BRkkq1atIhwOZ7oMEZGs4fP5ppynoS0REUmLgkRERNKiIBERkbQoSEREJC0KEhERSYutQRIMBnG73RiGQXt7+wXzx8bG2Lp1K4ZhUFVVxcDAQGpeW1sbhmHgdrvZt29f6vEf/vCHeL1e1q5dy7Zt2zh37pwdXRERkT+yLUgSiQTNzc3s3buXSCRCT08PkUhkUpuuri6WLl1Kf38/O3fupKWlBYBIJIJlWRw7doxgMMiOHTtIJBLE43F+9KMfEQ6HOXr0KIlEAsuy7OqSiIhgY5CEQiEMw6CsrIzc3Fz8fj+BQGBSm0AgQH19PQC1tbX09fWRTCYJBAL4/X7y8vIoLS3FMAxCoRAAExMTnD17lomJCc6cOcMXvvAFu7okIiLYGCTxeHzSLUKLi4sn3T70s20+uc/1yMjIlMu6XC6eeOIJVqxYwfLly8nPz+eee+656Po7Ozvx+Xz4fD6Gh4dnoYciIgtTVv+y/YMPPiAQCBCLxViyZAlf+9rXeOGFFy56N7mmpiaampqAS/9CU+a2VX+3x7Z1DbRvtm1dItnMtiMSl8vF4OBganpoaGjSfag/22ZiYoLR0VEKCwunXPZnP/sZpaWlFBUVcfXVV7NlyxZ+/etf29MhsZ2dIfLJ+uxep0g2si1IKioqiEajxGIxxsfHsSwL0zQntTFNk+7ubgB6e3uprq7G4XBgmiaWZTE2NkYsFiMajVJZWcmKFSs4cOAAZ86cIZlM0tfXR3l5uV1dEhtlcoOuMBG5NNuGtpxOJ7t27aKmpoZEIkFDQwNer5fW1lZ8Ph+madLY2Mj27dsxDIOCgoLUX2B5vV7q6urweDw4nU46OjrIycmhqqqK2tpabr31VpxOJxs2bEgNX4mIiD0W5D3bfT6frv6bRebKEYHOmchCdqntpn7ZLnPaXAkREZmagkRkmhRqIhenIJE5ay5uuOdiTSKZpiAREZG0KEhkTprLe/5zuTaRTFCQiIhIWhQkMudkwx5/NtQoYhcFiYiIpEVBInOK9vRFso+CROQKKfRE/kBBIiIiaVGQyJyhPXyR7KQgEUmDwk9EQSIiImlSkMickM179tlcu8hMUJCIiEhaFCQiIpIWW4MkGAzidrsxDIP29vYL5o+NjbF161YMw6CqqoqBgYHUvLa2NgzDwO12s2/fPgCOHz/O+vXrU/9df/31PPPMMzb1RmbKfBgamg99ELlStt2zPZFI0NzczP79+ykuLqaiogLTNPF4PKk2XV1dLF26lP7+fizLoqWlhRdffJFIJIJlWRw7dox3332Xu+66ixMnTuB2uzl8+HDq+V0uF/fff79dXRIREWw8IgmFQhiGQVlZGbm5ufj9fgKBwKQ2gUCA+vp6AGpra+nr6yOZTBIIBPD7/eTl5VFaWophGIRCoUnL9vX18ad/+qesXLnSri7JDNCevEj2sy1I4vE4JSUlqeni4mLi8fiUbZxOJ/n5+YyMjExrWcuy2LZt25Tr7+zsxOfz4fP5GB4enokuiYgI8+Rk+/j4OLt37+ZrX/valG2ampoIh8OEw2GKiopsrE4WCh1dyUJlW5C4XC4GBwdT00NDQ7hcrinbTExMMDo6SmFh4WWX3bt3L7feeis33njjLPdCREQ+y7YgqaioIBqNEovFGB8fx7IsTNOc1MY0Tbq7uwHo7e2luroah8OBaZpYlsXY2BixWIxoNEplZWVquZ6enksOa8ncpD14kfnBtr/acjqd7Nq1i5qaGhKJBA0NDXi9XlpbW/H5fJimSWNjI9u3b8cwDAoKCrAsCwCv10tdXR0ejwen00lHRwc5OTkAnD59mv379/P888/b1RUREfkURzKZTGa6CLv5fD7C4XCmy1jw5usRyUD75kyXIDLjLrXdnBcn20VEJHMUJJIR8/VoBOZ330QuRkEiIiJpUZCIiEhaFCQiIpIWBYnYbiGcQ1gIfRT5hIJERETSoiAREZG0KEhERCQtChKx1UI6d7CQ+ioLm4JERETSoiAR22gPXWR+UpCIiEhaFCQiIpIWBYnILNJwniwEChIREUmLgkRssZD3zBdy32VhsDVIgsEgbrcbwzBob2+/YP7Y2Bhbt27FMAyqqqoYGBhIzWtra8MwDNxuN/v27Us9/uGHH1JbW8uaNWsoLy/nN7/5jR1dERGRP7ItSBKJBM3Nzezdu5dIJEJPTw+RSGRSm66uLpYuXUp/fz87d+6kpaUFgEgkgmVZHDt2jGAwyI4dO0gkEgB861vf4t577+XNN9/kyJEjlJeX29UlERHBxiAJhUIYhkFZWRm5ubn4/X4CgcCkNoFAgPr6egBqa2vp6+sjmUwSCATw+/3k5eVRWlqKYRiEQiFGR0f5xS9+QWNjIwC5ubksWbLEri6JiAg2Bkk8HqekpCQ1XVxcTDwen7KN0+kkPz+fkZGRKZeNxWIUFRXxjW98gw0bNvDII49w+vTpi66/s7MTn8+Hz+djeHh4FnooIrIwZfXJ9omJCQ4ePMg3v/lNDh06xOLFiy967gWgqamJcDhMOBymqKjI5koXNp1s1msg85ttQeJyuRgcHExNDw0N4XK5pmwzMTHB6OgohYWFUy5bXFxMcXExVVVVwB+Gww4ePGhDb0RE5BO2BUlFRQXRaJRYLMb4+DiWZWGa5qQ2pmnS3d0NQG9vL9XV1TgcDkzTxLIsxsbGiMViRKNRKisruemmmygpKeH48eMA9PX14fF47OqSiIgATttW5HSya9cuampqSCQSNDQ04PV6aW1txefzYZomjY2NbN++HcMwKCgowLIsALxeL3V1dXg8HpxOJx0dHeTk5ADw3HPP8eCDDzI+Pk5ZWRn//M//bFeXREQEcCSTyWSmi7Cbz+cjHA5nuowFQecGJhto35zpEkSuyKW2m1l9sl1ERDJPQSIiImlRkMis0bCWyMKgIBERkbQoSERspKM0mY8UJCIikhYFiYiIpEVBIrNCQzgiC4eCRERE0qIgEbGZjtZkvlGQiIhIWhQkIiKSFgWJzDgN3VyeXiOZTxQkIiKSFgWJiIikRUEiIiJpUZCIiEhabA2SYDCI2+3GMAza29svmD82NsbWrVsxDIOqqioGBgZS89ra2jAMA7fbzb59+1KPr1q1inXr1rF+/Xp8Pp8d3ZBL0Enk6dNrJfOFbfdsTyQSNDc3s3//foqLi6moqMA0TTweT6pNV1cXS5cupb+/H8uyaGlp4cUXXyQSiWBZFseOHePdd9/lrrvu4sSJE6n7tr/66qvccMMNdnVFREQ+xbYjklAohGEYlJWVkZubi9/vJxAITGoTCASor68HoLa2lr6+PpLJJIFAAL/fT15eHqWlpRiGQSgUsqt0ERG5BNuCJB6PU1JSkpouLi4mHo9P2cbpdJKfn8/IyMgll3U4HNxzzz1s3LiRzs7OKdff2dmJz+fD5/MxPDw8k10TEVnQbBvami2vvfYaLpeL3/3ud9x9992sWbOG22+//YJ2TU1NNDU1AehciojIDLLtiMTlcjE4OJiaHhoawuVyTdlmYmKC0dFRCgsLL7nsJ/9ftmwZ999/v4a8Mkgnjz8/vWYyH9gWJBUVFUSjUWKxGOPj41iWhWmak9qYpkl3dzcAvb29VFdX43A4ME0Ty7IYGxsjFosRjUaprKzk9OnTfPTRRwCcPn2a//iP/2Dt2rV2dUlERLBxaMvpdLJr1y5qampIJBI0NDTg9XppbW3F5/NhmiaNjY1s374dwzAoKCjAsiwAvF4vdXV1eDwenE4nHR0d5OTk8P7773P//fcDfziC+frXv869995rV5fkU7RnLbJwOZLJZDLTRdjN5/MRDoczXca8oiC5cgPtmzNdgshlXWq7qV+2i4hIWhQkIhmmoznJdgoSERFJi4JE0qY96vTpNZRspiAREZG0fO4gOX36NIlEYjZqERGRLHTZIPn444/5yU9+wubNm1m2bBlr1qxh+fLleDwe/vZv/5b+/n476hQRkTnqskFyxx138D//8z+0tbXxv//7vwwODvK73/2O1157jdtuu42WlhZeeOEFO2oVEZE56LI/SDx//jxXX331pMfGx8fJzc29ZJu5TD9InDk6STyz9ONEmavS+kHixQLiySefTP37V7/6VVaFiIiIzKwrutZWTU0NnZ2dLF68mP/+7//mS1/60kzXJSIiWeJzB0ljYyP5+fkcOnSI2267jba2ttmoS0REssTnDpKuri7Onj3LwYMH+a//+i8ee+wxnn/++dmoTUREssC0f0eybt06HnzwQb73ve/x85//nJUrV3L69GmFyAKmE+0zT6+pZKNpB8l//ud/8uijj7Jo0SIsy2Lt2rW8/PLLs1mbiIhkgWkPbRUUFLBp0yY2bdoEQDQa5amnnpqtukREJEtM+4jkxIkTk6ZXr17NG2+8MeMFSXbQEIyIfGLaQfLYY4+xYsUKvvjFL/LYY49RX1/P2rVrOXPmzLRXFgwGcbvdGIZBe3v7BfPHxsbYunUrhmFQVVXFwMBAal5bWxuGYeB2u9m3b9+k5RKJBBs2bOC+++6bdi0iIjIzph0kr776Ku+88w4vvvgi9913H4ZhcPbsWdavX8+aNWsuu3wikaC5uZm9e/cSiUTo6ekhEolMatPV1cXSpUvp7+9n586dtLS0ABCJRLAsi2PHjhEMBtmxY8ekC0c+++yzlJeXT7crInOajvYk21w2SD57BZUVK1bwl3/5l3znO9+ht7eXEydOEAqFLruiUCiEYRiUlZWRm5uL3+8nEAhMahMIBKivrwegtraWvr4+kskkgUAAv99PXl4epaWlGIaRWufQ0BB79uzhkUcemXanRURk5kzroo3PPfcc77zzzqTHx8fHeeWVV6ivr+df//VfL7uieDxOSUlJarq4uJh4PD5lG6fTSX5+PiMjI5dc9vHHH+f73/8+V1116a50dnbi8/nw+XwMDw9ftl4REZmeywZJMBgkJyeHbdu2pS4fX1payurVq+np6eHxxx/n4YcftqHUC/30pz9l2bJlbNy48bJtm5qaCIfDhMNhioqKbKhu/tLQy+zTayzZ5LJ//nvNNdewY8cOduzYwfnz5zl58iSLFi1iyZIln2tFLpeLwcHB1PTQ0BAul+uibYqLi5mYmGB0dJTCwsIpl929eze7d+/m5Zdf5ty5c/zf//0fDz30kC5rLyJio2mfbN+7dy9f+cpX2LRpE01NTRw4cOBzraiiooJoNEosFmN8fBzLsjBNc1Ib0zTp7u4GoLe3l+rqahwOB6ZpYlkWY2NjxGIxotEolZWVtLW1MTQ0xMDAAJZlUV1drRAREbHZtH+QuGPHDl544QU8Hg+//e1veeKJJ2hubmbbtm3TW5HTya5du6ipqSGRSNDQ0IDX66W1tRWfz4dpmjQ2NrJ9+3YMw6CgoADLsgDwer3U1dXh8XhwOp10dHSQk5NzZT0WEZEZddkbW32iqqqK119/PTV9+vRpqqqqOHr06KwVN1t0Y6srp7F7e+lGVzJXpHVjq0+UlpbyD//wD4yPjwN/uOGV03lFtzMREZF5ZNpBctVVV/Fv//ZvlJSU8OUvfxnDMNi0aRPRaHQ26xMRkTlu2ocUP/nJT4A/XMbk6NGjHDlyhCNHjvDoo4/y1ltvXfA7ExERWRg+99hUXl4eGzdunNZvN0REZP6b9tCWiE6020+vuWQDBYmIiKRFQSIiImlRkMi0aIhFRKaiIBERkbQoSETmOB0NylynIBERkbQoSOSytEcsIpeiIBERkbQoSESygI4KZS5TkIiISFoUJCIikhYFiVyShlTmDr0XMlfZGiTBYBC3241hGLS3t18wf2xsjK1bt2IYBlVVVQwMDKTmtbW1YRgGbrebffv2AXDu3DkqKyu55ZZb8Hq9PPnkk3Z1RURE/si2IEkkEjQ3N7N3714ikQg9PT1EIpFJbbq6uli6dCn9/f3s3LmTlpYWACKRCJZlcezYMYLBIDt27CCRSJCXl8crr7zCkSNHOHz4MMFgkAMHDtjVJRERwcYgCYVCGIZBWVkZubm5+P1+AoHApDaBQID6+noAamtr6evrI5lMEggE8Pv95OXlUVpaimEYhEIhHA4H1157LQDnz5/n/PnzOBwOu7okIiLYGCTxeJySkpLUdHFxMfF4fMo2TqeT/Px8RkZGLrlsIpFg/fr1LFu2jLvvvpuqqqqLrr+zsxOfz4fP52N4eHimuzcvaUx+7tF7InNR1p9sz8nJ4fDhwwwNDREKhTh69OhF2zU1NREOhwmHwxQVFdlcpYjI/GVbkLhcLgYHB1PTQ0NDuFyuKdtMTEwwOjpKYWHhtJZdsmQJd9xxB8FgcBZ7ISIin2VbkFRUVBCNRonFYoyPj2NZFqZpTmpjmibd3d0A9Pb2Ul1djcPhwDRNLMtibGyMWCxGNBqlsrKS4eFhPvzwQwDOnj3L/v37WbNmjV1dEhERwGnbipxOdu3aRU1NDYlEgoaGBrxeL62trfh8PkzTpLGxke3bt2MYBgUFBViWBYDX66Wurg6Px4PT6aSjo4OcnBzee+896uvrSSQSfPzxx9TV1XHffffZ1aV5TWPxc9eqv9vDQPvmTJchkuJIJpPJTBdhN5/PRzgcznQZc5qCZG5TkIjdLrXdzPqT7TLzFCIi8nkoSEREJC0KEhERSYuCRCQLafhR5hIFiYiIpEVBIpNoTzd76L2SuUJBIiIiaVGQiIhIWhQkIiKSFgWJpGjMPfvoPZO5QEEiIiJpUZCIiEhaFCQCaIgkm+m9k0xTkIiISFoUJKI9WhFJi4JERETSoiARmQd0VCmZZGuQBINB3G43hmHQ3t5+wfyxsTG2bt2KYRhUVVUxMDCQmtfW1oZhGLjdbvbt2wfA4OAgd9xxBx6PB6/Xy7PPPmtXV0RE5I9sC5JEIkFzczN79+4lEonQ09NDJBKZ1Karq4ulS5fS39/Pzp07aWlpASASiWBZFseOHSMYDLJjxw4SiQROp5Onn36aSCTCgQMH6OjouOA55dK0Jzt/6L2UTLEtSEKhEIZhUFZWRm5uLn6/n0AgMKlNIBCgvr4egNraWvr6+kgmkwQCAfx+P3l5eZSWlmIYBqFQiOXLl3PrrbcCcN1111FeXk48HrerSyIigo1BEo/HKSkpSU0XFxdfsNH/dBun00l+fj4jIyPTWnZgYIBDhw5RVVV10fV3dnbi8/nw+XwMDw/PVLdERBa8eXGy/dSpUzzwwAM888wzXH/99Rdt09TURDgcJhwOU1RUZHOFc5OGQuYfvaeSCbYFicvlYnBwMDU9NDSEy+Wass3ExASjo6MUFhZectnz58/zwAMP8OCDD7JlyxYbeiIiIp9mW5BUVFQQjUaJxWKMj49jWRamaU5qY5om3d3dAPT29lJdXY3D4cA0TSzLYmxsjFgsRjQapbKykmQySWNjI+Xl5Xz729+2qysiIvIpTttW5HSya9cuampqSCQSNDQ04PV6aW1txefzYZomjY2NbN++HcMwKCgowLIsALxeL3V1dXg8HpxOJx0dHeTk5PDaa6/xL//yL6xbt47169cD8I//+I/8xV/8hV3dEhFZ8BzJZDKZ6SLs5vP5CIfDmS4jozSWPr8NtG/OdAkyz1xquzkvTrbL56MQEZGZpCAREZG0KEhE5iEddYqdFCQLjDYwIjLTFCQiIpIWBYnIPKWjT7GLgmQB0YZl4dF7LnZQkIiISFoUJCIikhYFyQKhIY6FS++9zDYFiYiIpEVBIiIiaVGQLAAa2hB9BmQ2KUjmOW1ARGS2KUhEFgjtVMhsUZDMY9pwyGfpMyGzQUEiIiJpsTVIgsEgbrcbwzBob2+/YP7Y2Bhbt27FMAyqqqoYGBhIzWtra8MwDNxuN/v27Us93tDQwLJly1i7dq0dXRARkc+wLUgSiQTNzc3s3buXSCRCT08PkUhkUpuuri6WLl1Kf38/O3fupKWlBYBIJIJlWRw7doxgMMiOHTtIJBIAPPzwwwSDQbu6kTU0hCFT0WdDZpptQRIKhTAMg7KyMnJzc/H7/QQCgUltAoEA9fX1ANTW1tLX10cymSQQCOD3+8nLy6O0tBTDMAiFQgDcfvvtFBQU2NUNERH5DNuCJB6PU1JSkpouLi4mHo9P2cbpdJKfn8/IyMi0lr2czs5OfD4fPp+P4eHhNHoy92mPUy5HnxGZSQvmZHtTUxPhcJhwOExRUVGmy5k12kCIiN1sCxKXy8Xg4GBqemhoCJfLNWWbiYkJRkdHKSwsnNayIvL5aKdDZoptQVJRUUE0GiUWizE+Po5lWZimOamNaZp0d3cD0NvbS3V1NQ6HA9M0sSyLsbExYrEY0WiUyspKu0rPGtowiEgm2BYkTqeTXbt2UVNTQ3l5OXV1dXi9XlpbW9m9ezcAjY2NjIyMYBgGP/jBD1J/Iuz1eqmrq8Pj8XDvvffS0dFBTk4OANu2beOLX/wix48fp7i4mK6uLru6JJL1tPMhM8GRTCaTmS7Cbj6fj3A4nOkyZpQ2CJKOgfbNmS5B5rhLbTcXzMl2ERGZHQqSLLfq7/boaETSps+QpENBIiKAwkSunIIki+mLLyJzgYIkSylEZDZoqFSuhIJERETSoiDJQtpjlNmmz5h8HgqSLKMvuNhFnzWZLgVJFtEXW+ymz5xMh4IkS+gLLZmiz55cjoIkC+iLLJmmz6BcioJkDtOfYspcos+iTEVBMkfpSytzkXZu5GKcmS5AJtOXVLLBJ59TXTVYQEckc4b29CQb6TMroCOSjNMXUbLdpz/DOkJZmBQkGaIAkflIQ14Lk61DW8FgELfbjWEYqdvoftrY2Bhbt27FMAyqqqoYGBhIzWtra8MwDNxuN/v27Zv2c84VnwxdaQhLFgJ93hcW245IEokEzc3N7N+/n+LiYioqKjBNE4/Hk2rT1dXF0qVL6e/vx7IsWlpaePHFF4lEIliWxbFjx3j33Xe56667OHHiBMBlnzMT9MURmexi3wkdtcwftgVJKBTCMAzKysoA8Pv9BAKBSRv9QCDAd7/7XQBqa2v567/+a5LJJIFAAL/fT15eHqWlpRiGQSgUArjsc840hYTIzPi83yUFz9xlW5DE43FKSkpS08XFxbz++utTtnE6neTn5zMyMkI8Hue2226btGw8Hge47HN+orOzk87OTgDefPNNfD7fFfXjhs/Rdnh4mKKioitaT7aY731U/+YOn+/JK1oum/p4Jezq36dPNXzWgjnZ3tTURFNTk63r9Pl8hMNhW9dpt/neR/Uv+833Ps6F/tl2st3lcjE4OJiaHhoawuVyTdlmYmKC0dFRCgsLp1x2Os8pIiKzy7YgqaioIBqNEovFGB8fx7IsTNOc1MY0Tbq7uwHo7e2luroah8OBaZpYlsXY2BixWIxoNEplZeW0nlNERGaXbUNbTqeTXbt2UVNTQyKRoKGhAa/XS2trKz6fD9M0aWxsZPv27RiGQUFBAZZlAeD1eqmrq8Pj8eB0Ouno6CAnJwfgos85V9g9lJYJ872P6l/2m+99nAv9cySTyWSmixARkeyla22JiEhaFCQiIpIWBcksyZZLt1yJwcFB7rjjDjweD16vl2effTbTJc2KRCLBhg0buO+++zJdyqz48MMPqa2tZc2aNZSXl/Ob3/wm0yXNqB/+8Id4vV7Wrl3Ltm3bOHfuXKZLSltDQwPLli1j7dq1qcd+//vfc/fdd7N69WruvvtuPvjgA9vrUpDMgk8uB7N3714ikQg9PT1EIpFMlzVjnE4nTz/9NJFIhAMHDtDR0TGv+veJZ599lvLy8kyXMWu+9a1vce+99/Lmm29y5MiRedXXeDzOj370I8LhMEePHiWRSKT+eCebPfzwwwSDwUmPtbe3c+eddxKNRrnzzjszsuOqIJkFn74cTG5uburSLfPF8uXLufXWWwG47rrrKC8vT11pYL4YGhpiz549PPLII5kuZVaMjo7yi1/8gsbGRgByc3NZsmRJZouaYRMTE5w9e5aJiQnOnDnDF77whUyXlLbbb7+dgoKCSY8FAgHq6+sBqK+v56WXXrK9LgXJLLjY5WDm24b2EwMDAxw6dIiqqqpMlzKjHn/8cb7//e9z1VXz8ysSi8UoKiriG9/4Bhs2bOCRRx7h9OnTmS5rxrhcLp544glWrFjB8uXLyc/P55577sl0WbPi/fffZ/ny5QDcdNNNvP/++7bXMD+/JWKLU6dO8cADD/DMM89w/fXXZ7qcGfPTn/6UZcuWsXHjxkyXMmsmJiY4ePAg3/zmNzl06BCLFy+eV+fyPvjgAwKBALFYjHfffZfTp0/zwgsvZLqsWedwOHA4HLavV0EyCxbCpVvOnz/PAw88wIMPPsiWLVsyXc6M+tWvfsXu3btZtWoVfr+fV155hYceeijTZc2o4uJiiouLU0eStbW1HDx4MMNVzZyf/exnlJaWUlRUxNVXX82WLVv49a9/nemyZsWNN97Ie++9B8B7773HsmXLbK9BQTIL5vulW5LJJI2NjZSXl/Ptb3870+XMuLa2NoaGhhgYGMCyLKqrq+fd3uxNN91ESUkJx48fB6Cvry/j9/GZSStWrODAgQOcOXOGZDJJX1/fvPpjgk/79KWluru7+epXv2p/EUmZFXv27EmuXr06WVZWlnzqqacyXc6M+uUvf5kEkuvWrUvecsstyVtuuSW5Z8+eTJc1K1599dXk5s2bM13GrDh06FBy48aNyXXr1iW/+tWvJn//+99nuqQZ1dramnS73Umv15t86KGHkufOnct0SWnz+/3Jm266Kel0OpMulyv5T//0T8mTJ08mq6urk4ZhJO+8887kyMiI7XXpEikiIpIWDW2JiEhaFCQiIpIWBYmIiKRFQSIiImlRkIiISFoUJCIikhYFiYiIpEVBIpJBmzZt4s033wRgZGRk0n0mRLKFgkQkg/r7+/mzP/szAN544w3WrVuX4YpEPj8FiUiGvP3227hcrtSl6t944w1uvvlm4A93whPJFgoSkQw5cuRIKjgAfvvb33LzzTdz5swZ8vPzefXVV2lpaZkXt4iV+U1BIpIhhw8fToVENBolEAiwbt06Dh48yKFDhzh+/Djf+973uOaaazJcqcilOTNdgMhCdeTIEa655hpuueUWbr75ZjweD93d3Vx77bVUVVWxePHiTJcoMi0KEpEMeeONNzh48CDXXXfdpMcfffRRnn/+eb7zne/wy1/+kq985SsZqlBkenQZeZEM+Oijj9i4cSMnTpzIdCkiaVOQiIhIWnSyXURE0qIgERGRtChIREQkLQoSERFJi4JERETSoiAREZG0KEhERCQt/w/EGDHILQeODAAAAABJRU5ErkJggg==\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "p_mu = norm.pdf(mu_k, loc=mu_0, scale = sigma_0)\n",
    "p_mu = p_mu/sum(p_mu)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "fig.patch.set_facecolor('white')\n",
    "\n",
    "ax.bar(mu_k, p_mu)\n",
    "ax.set_xlabel('$\\mu_k$')\n",
    "ax.set_ylabel('$p(\\mu_k)$')\n",
    "ax.set_title(f'Prior')\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "And now loop for each posterior"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "fig, ax = plt.subplots(nrows=len(x), figsize=(9,20),\n",
    "                       sharex='all', sharey='all')\n",
    "fig.patch.set_facecolor('white')\n",
    "\n",
    "this_p_mu = p_mu # Set initial prior\n",
    "\n",
    "for i,x_i in enumerate(x):\n",
    "    # Likelihood\n",
    "    p_x_mu = norm.pdf(x_i, loc=mu_k, scale=sigma_x)\n",
    "\n",
    "    # Unnormalized posterior\n",
    "    p_mu_x = p_x_mu*this_p_mu\n",
    "\n",
    "    # Normalized posterior\n",
    "    p_mu_x = p_mu_x / sum(p_mu_x)\n",
    "\n",
    "    ax[i].bar(mu_k, p_mu_x)\n",
    "    ax[i].set_xlabel('$\\mu_k$')\n",
    "    ax[i].set_ylabel('$p(\\mu_k | x_0)$')\n",
    "    ax[i].set_title(f'x[{i}] = {x_i}')\n",
    "\n",
    "    # Now set the prior to the previous posterior\n",
    "    this_p_mu = p_mu_x"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "execution_count": 54,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 648x1440 with 6 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAR+CAYAAADZUZ+AAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAACS3UlEQVR4nOzdf3SU9Z33/9fASIIEQoKwi5NoCFcbTQiNMjH0gPc3cLhNOeh4s7Ak/oAoMVEBUXq6d9xzt5Qtegh7Wrd7N7I1NO4GC4x74pYBJQGMYLE9GAPKr9QywUSTkRaIFQQhIXG+f3iTNYUkE5if1zwf5+ToNfO5ruv9uWbq59XP9WMsXq/XKwAAABMaEuoCAAAAAoWgAwAATIugAwAATIugAwAATIugAwAATIugAwAATIugAwAATIugAyAspaSkaPjw4Vq4cKFP7SsrKxUXFyeLxaKmpqYAVwcgUhB0AIStbdu26ZVXXulZbmlp0YwZM3TjjTfqtttu05tvvtnzXlFRkc6dOxewWv71X/9VEyZM0IgRI3T77bfr2LFjAdsXAP8h6ACIGA888IDuuOMOtbe36/nnn9f8+fN16tSpgO/3V7/6lSorK/XGG2/o3Llzev3113XTTTcFfL8Arh9BB0BIHD9+XImJiTpw4IAk6dNPP9XYsWO1Z8+eq7Y/duyYDhw4oH/6p3/S8OHDNW/ePGVmZuq1114LaJ1fffWV/umf/kn/8i//ovT0dFksFk2cOFGJiYkB3S8A/yDoAAiJiRMnau3atXr44Yf15Zdf6tFHH1VhYaFyc3Ov2v7o0aNKTU3VyJEje177zne+o6NHj/q8z3vvvVejR4++6t+999571XXa2trU1tamI0eOKDk5WRMmTNCPf/xjffXVV4PqL4DQsIa6AADRq7i4WNu2bVNOTo4sFou2bt3aZ9tz584pPj6+12vx8fHyeDw+7+/1118fdI1tbW2SpJ07d+rw4cP6/PPPdc899ygpKUnFxcWD3h6A4GJGB0BIFRcX68iRI3rqqacUExPTZ7u4uDidPXu212tnz57tNcMTCMOHD5ck/e///b81evRopaSk6PHHH9f27dsDul8A/kHQARAy586d0zPPPKOioiKtWrVKn332WZ9tMzIy9NFHH+mLL77oee3gwYPKyMjweX+zZ89WXFzcVf9mz5591XXS0tI0bNgwWSyWnte++e8AwhtBB0DIPP3007Lb7frVr36lOXPm6Iknnuiz7be//W1lZWXpn/7pn3Tx4kX95je/0aFDhzRv3jyf91dTU6Nz585d9a+mpuaq69x4443Kz8/XP//zP+uLL75QW1ubKioq+rymB0B4IegACAmXy6Xa2lr927/9myTphRde0IEDB7Rx48Y+13E6nWpoaFBCQoKeffZZVVdXa+zYsQGvtby8XHFxcbr55pv13e9+Vw8++KAWL14c8P0CuH4Wr9frDXURAPDX0tLSdOLECc2dO1dVVVUDtv/3f/93rVixQhcvXlRjY6NSU1ODUCWAcEfQAQAApsWpKwAAYFoEHQAAYFoEHQAAYFpR+WTkm266SSkpKaEuAwAA+EFLS4tOnz591feiMuikpKSooaEh1GUAAAA/sNvtfb7HqSsAAGBaBB0AAGBaBB0AAGBaBB0AAGBaBB0AAGBaBB0AAGBaBB0AAGBaBB0AAGBaBB0AAGBaQQ06tbW1SktLk2EYKisru+L9jo4O5efnyzAM5eTkqKWlRZK0a9cuTZkyRZmZmZoyZYreeuutnnX279+vzMxMGYah5cuXy+v1Bqs7AAAgzAUt6HR3d2vp0qWqqalRY2OjNm/erMbGxl5tKisrlZCQoKamJq1YsUKlpaWSvv5tqm3btunw4cOqqqrSwoULe9Z58skntX79erndbrndbtXW1garSwAAIMwFLejU19fLMAylpqZq2LBhKigokMvl6tXG5XKpsLBQkjR//nzV1dXJ6/Xqjjvu0M033yxJysjI0IULF9TR0aETJ07o7Nmzmjp1qiwWixYtWqQtW7YEq0sAACDMBS3oeDweJScn9ywnJSXJ4/H02cZqtSo+Pl7t7e292rz22mu68847FRMTI4/Ho6SkpH63CQAAoldE/Xr50aNHVVpaqp07dw563YqKClVUVEiSTp065e/SAABAGArajI7NZlNra2vPcltbm2w2W59turq6dObMGY0ZM6an/dy5c7VhwwZNnDixp31bW1u/27yspKREDQ0Namho0NixY/3aNwAAEJ6CFnSys7PldrvV3Nyszs5OOZ1OORyOXm0cDoeqqqokSdXV1Zo5c6YsFos+//xzzZkzR2VlZZo2bVpP+/Hjx2vUqFHat2+fvF6vNmzYoPvvvz9YXQIAAGEuaEHHarWqvLxceXl5uv3227VgwQJlZGRo5cqV2rp1qySpqKhI7e3tMgxDL7zwQs8t6OXl5WpqatJPfvITZWVlKSsrSydPnpQkrVu3To899pgMw9DEiRM1e/bsYHUJAACEOYs3Ch88Y7fb1dDQEOoyAACAH/Q3rvNkZAAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFpBDTq1tbVKS0uTYRgqKyu74v2Ojg7l5+fLMAzl5OSopaVFktTe3q4ZM2YoLi5Oy5Yt67VObm6u0tLSlJWVpaysLJ08eTIYXQEAABHAGqwddXd3a+nSpdq1a5eSkpKUnZ0th8Oh9PT0njaVlZVKSEhQU1OTnE6nSktL9eqrryo2NlarV6/WkSNHdOTIkSu2vXHjRtnt9mB1BQAARIigzejU19fLMAylpqZq2LBhKigokMvl6tXG5XKpsLBQkjR//nzV1dXJ6/VqxIgRmj59umJjY4NVLgAAMIGgBR2Px6Pk5OSe5aSkJHk8nj7bWK1WxcfHq729fcBtP/roo8rKytLq1avl9Xr9WzgAAIhYEX8x8saNG3X48GHt3btXe/fu1SuvvHLVdhUVFbLb7bLb7Tp16lSQqwQAAKEQtKBjs9nU2tras9zW1iabzdZnm66uLp05c0ZjxowZcLuSNHLkSD344IOqr6+/aruSkhI1NDSooaFBY8eOvZ6uAACACBG0oJOdnS23263m5mZ1dnbK6XTK4XD0auNwOFRVVSVJqq6u1syZM2WxWPrcZldXl06fPi1JunTpkl5//XVNmjQpcJ0AAAARJWh3XVmtVpWXlysvL0/d3d1avHixMjIytHLlStntdjkcDhUVFWnhwoUyDEOJiYlyOp0966ekpOjs2bPq7OzUli1btHPnTt16663Ky8vTpUuX1N3drVmzZqm4uDhYXQIAAGHO4o3Cq3ftdrsaGhpCXQYAAPCD/sb1iL8YGQAAoC8EHQAAYFoEHQAAYFoEHQAAYFqDDjrnz59Xd3d3IGoBAADwqwGDzldffaVNmzZpzpw5GjdunG677TaNHz9e6enp+od/+Ac1NTUFo04AAIBBGzDozJgxQ8ePH9eaNWv0pz/9Sa2trTp58qTeeecdTZ06VaWlpfr1r38djFoBAAAGZcDn6Fy6dEk33HBDvxvxpU044Tk6AACYx3U9R8eXABNJIQcAAESPQf8ExKZNm7R161YNHTpUXq9X9913nx544IFA1AYAAHBdBh103n777V6/QbV06VKCDgAACEuDDjodHR164403lJycrLa2Nl24cCEQdQEAAFy3QT9HZ926dfrLX/6i7du367PPPlN5eXkg6gIAALhuPs/oPP300/r5z3+uG2+8UQ8//HAgawIAAPALn2d0Ro4cKYfDofPnz0uSduzYoWnTpgWsMAAAgOvl84zOc889p02bNik3N1fDhg1TXFycysrKAlkbAADAdfE56NTV1Wn9+vUaMWKETpw4oZdffllpaWmBrA0AAOC6+Hzq6vnnn9fq1au1Z88eVVdXKz8/X2+99VYgawMAALguPs/ofDPUZGZmqqamRvPmzdPvf//7gBQGAABwvQZ9e/ll48ePV11dnT9rAQAA8KsBZ3Q++eQTnzY0evRojRo16roLAgAA8JcBg05hYaEsFov6+5Fzi8WiRx55RIsWLfJrcQAiX8qzb4S6hAG1lM0JdQkAAmTAoLN79+4rXuvo6FBMTExACgJgDpEQcC67XCuBBzCfa7pGZ9WqVT3//rvf/c5ftQAwgZRn34iokPNNkVw7gKu7pqCTl5eniooKbdy4UTU1Nf6uCUCEMktIMEs/AFzDr5cXFRUpPj5e77//vqZOnao1a9YEoi4AEcZs4SDl2Tc4lQWYwKBndCorK/X888/rueee09/8zd/o8ccf93nd2tpapaWlyTCMq/58REdHh/Lz82UYhnJyctTS0iJJam9v14wZMxQXF6dly5b1Wmf//v3KzMyUYRhavnx5vxdNA/A/M5/uMXPfgGjhc9DJzMzUQw89pLVr12rPnj269dZbdf78eb300ks+rd/d3a2lS5eqpqZGjY2N2rx5sxobG3u1qaysVEJCgpqamrRixQqVlpZKkmJjY7V69Wr99Kc/vWK7Tz75pNavXy+32y23263a2lpfuwTgOkVLCIiWfgJm5HPQefvtt1VcXKzhw4fL6XRq0qRJ2r59u887qq+vl2EYSk1N1bBhw1RQUCCXy9WrjcvlUmFhoSRp/vz5qqurk9fr1YgRIzR9+nTFxsb2an/ixAmdPXtWU6dOlcVi0aJFi7RlyxafawIAXxF2gMjkc9BJTExUbm6uli9frqqqKr333nsyDMPnHXk8HiUnJ/csJyUlyePx9NnGarUqPj5e7e3t/W4zKSmp320CCIxoHPijsc9ApPM56Bw7dqzX8re+9S0dOnTI7wUFSkVFhex2u+x2u06dOhXqcoCIxoAPIFL4fNfV448/ruPHj8tms2ny5Mm6ePGiJk2apC+//FI33njjgOvbbDa1trb2LLe1tclms121TVJSkrq6unTmzBmNGTOm3222tbX1u83LSkpKVFJSIkmy2+0D1gvg6qI95HA3FhBZfJ7R2b17tz755BO9+uqruvfee2UYhi5cuKCsrCzddtttA66fnZ0tt9ut5uZmdXZ2yul0yuFw9GrjcDhUVVUlSaqurtbMmTNlsVj63Ob48eM1atQo7du3T16vVxs2bND999/va5cADFK0h5zLOA5A5LB4B7gf2+v19hs2JOns2bM+/aDn9u3b9cwzz6i7u1uLFy/W//k//0crV66U3W6Xw+HQxYsXtXDhQr3//vtKTEyU0+lUamqqJCklJUVnz55VZ2enRo8erZ07dyo9PV0NDQ165JFHdOHCBc2ePVu/+MUvBqzXbreroaFhwHoB/DcG9ysxswOEh/7G9QGDTm5urubNm6f7779ft9xyS8/rnZ2deuedd1RVVaUZM2bokUce8WvRgUTQAQaHkNM3wg4Qev2N6wOeuqqtrdXQoUP1wAMPaPz48UpPT9eECRP0rW99S5s3b9YzzzwTUSEHAPyJEAiEtwEvRo6NjdWSJUu0ZMkSXbp0SadPn9bw4cM1evToIJQHAABw7Xy+GLmmpkZ33323cnNzVVJSon379gWyLgBhghkLAJHM56CzZMkS/exnP9O+fftUUlKiH/zgB9q8eXMgawMQYoQc33CcgPDlc9AZN26cpk2bpoSEBM2aNUs7duzQ888/H8jaAIQQg/fgcLyA8ORz0JkwYYJ++MMfqrOzU5J0ww03yGr1+XmDACIIg/a14bgB4cfnoDNkyBD95je/UXJysqZPny7DMJSbmyu32x3I+gAEGYP19eH4AeHF5ymZTZs2SZI6Ojp05MgRHTx4UAcPHlRxcbE++ugjffLJJwErEgAA4FoM+txTTEyMpkyZoilTpgSiHgCIePweFhA+fD51BcD8OO3iPxxLIDwQdABIYmAGYE4EHQAAYFoEHQDM5gQIxxUIPYIOEOUYjAOL4wuEFkEHiGIMwsHBcQZCh6ADAABMi6ADRClmGYKL4w2EBkEHAIKEsAMEH0EHAACYFkEHiELMLACIFgQdIMoQckKL4w8EF0EHiCIMsuGBzwEIHoIOAAAwLYIOECWYRQgvfB5AcBB0gCjAoBqe+FyAwCPoAEAIEXaAwApq0KmtrVVaWpoMw1BZWdkV73d0dCg/P1+GYSgnJ0ctLS09761Zs0aGYSgtLU07duzoeT0lJUWZmZnKysqS3W4PRjcAAECECFrQ6e7u1tKlS1VTU6PGxkZt3rxZjY2NvdpUVlYqISFBTU1NWrFihUpLSyVJjY2NcjqdOnr0qGpra7VkyRJ1d3f3rLd792598MEHamhoCFZ3gIjBjEH44zMCAidoQae+vl6GYSg1NVXDhg1TQUGBXC5XrzYul0uFhYWSpPnz56uurk5er1cul0sFBQWKiYnRhAkTZBiG6uvrg1U6ELEYQAFEu6AFHY/Ho+Tk5J7lpKQkeTyePttYrVbFx8ervb2933UtFovuueceTZkyRRUVFUHoCRAZCDmRhc8LCIyIvxj5nXfe0YEDB1RTU6MXX3xRv/3tb6/arqKiQna7XXa7XadOnQpylUBwMWhGJj43wP+CFnRsNptaW1t7ltva2mSz2fps09XVpTNnzmjMmDH9rnv5n+PGjdPcuXP7PKVVUlKihoYGNTQ0aOzYsX7tGwAACE9BCzrZ2dlyu91qbm5WZ2ennE6nHA5HrzYOh0NVVVWSpOrqas2cOVMWi0UOh0NOp1MdHR1qbm6W2+3WXXfdpfPnz+uLL76QJJ0/f147d+7UpEmTgtUlICwxKxDZ+PwA/7IGbUdWq8rLy5WXl6fu7m4tXrxYGRkZWrlypex2uxwOh4qKirRw4UIZhqHExEQ5nU5JUkZGhhYsWKD09HRZrVa9+OKLGjp0qP785z9r7ty5kr6eAXrwwQf1ve99L1hdAsIOg6Q5pDz7hlrK5oS6DMAULF6v1xvqIoLNbrdzKzpMiaBjLoQdwDf9jesRfzEygK8RcsyHzxS4fgQdwAQYEAHg6gg6AADAtAg6QIRjNsfc+HyB60PQASIYg2B04HMGrh1BB4hQDH7Rhc8buDYEHQCIEIQdYPAIOkAEYsCLXnz2wOAQdIAIw0AHAL4j6AARhJADie8BMBgEHSBCMLjhm/g+AL4h6AARgEENV8P3AhgYQQcIcwxm6E/Ks2/wHQH6QdABwhgDGHzFdwW4OoIOEKYYuADg+hF0gDBEyMG14HsDXMka6gIA/DcGKlyvy9+hlrI5Ia4ECA/M6ACACRGaga8RdIAwwJ0zCAS+UwBBBwg5BiMEEt8vRDuu0QFChAEIwcJ1O4hmzOgAIUDIQShwihTRiKADBBEDDcIB30FEE05dAUHAwIJww+ksRAuCDhBABByEu29+Rwk9MCOCDuBnhBtEKmZ5YEYEHcAPCDcwk7/+PhN8EMmCejFybW2t0tLSZBiGysrKrni/o6ND+fn5MgxDOTk5amlp6XlvzZo1MgxDaWlp2rFjh8/bBALh8kXFXFyMaMB3HZEsaDM63d3dWrp0qXbt2qWkpCRlZ2fL4XAoPT29p01lZaUSEhLU1NQkp9Op0tJSvfrqq2psbJTT6dTRo0f16aefatasWTp27JgkDbhN4HrwH3agt77+N8GsD8JV0IJOfX29DMNQamqqJKmgoEAul6tXKHG5XFq1apUkaf78+Vq2bJm8Xq9cLpcKCgoUExOjCRMmyDAM1dfXS9KA2wS+ieACBMa1/G+LcIRgCFrQ8Xg8Sk5O7llOSkrSu+++22cbq9Wq+Ph4tbe3y+PxaOrUqb3W9Xg8kjTgNi+rqKhQRUWFJOnDDz+U3W73T8f87NSpUxo7dmyoyzCVbx7Tm0Jci1nwPfW/aDymdvuPA7btaDyegRbOx/Sbl7r8tai5GLmkpEQlJSWhLmNAdrtdDQ0NoS7DVDim/scx9T+OqX9xPP0vUo9p0C5Gttlsam1t7Vlua2uTzWbrs01XV5fOnDmjMWPG9LmuL9sEAADRK2hBJzs7W263W83Nzers7JTT6ZTD4ejVxuFwqKqqSpJUXV2tmTNnymKxyOFwyOl0qqOjQ83NzXK73brrrrt82iYAAIheQTt1ZbVaVV5erry8PHV3d2vx4sXKyMjQypUrZbfb5XA4VFRUpIULF8owDCUmJsrpdEqSMjIytGDBAqWnp8tqterFF1/U0KFDJemq24xkkXB6LdJwTP2PY+p/HFP/4nj6X6QeU4vX6/WGuggAAIBA4NfLAQCAaRF0AJhCbm6uPvzwQ0lSe3u7Jk2aFOKKAIQDgg4AU2hqatK3v/1tSdKhQ4eUmZkZ4ooAhAOCDoCI9/HHH8tms2nIkK//k3bo0CFNnjxZkrR48eJQlgYgxAg6ACLewYMHe4KNJO3fv1+TJ0/Wl19+qfj4eO3evVulpaW6ePFiCKsEEAoEHQAR74MPPugJMW63Wy6XS5mZmTpw4IDef/99/fGPf9TatWsVGxsb4koBBFvU/AQEAPM6ePCgYmNj9Z3vfEeTJ09Wenq6qqqqFBcXp5ycHI0YMSLUJQIIEYIOgIh36NAhHThwQCNHjuz1enFxsV566SX96Ec/0t69e3X33XeHqEIAocIDAwFEtC+++EJTpkzRsWPHQl0KgDBE0AEAAKbFxcgAAMC0CDoAwlJKSoqGDx+uhQsX+tT+zTffVFxcnIYMGaI333wzwNUBiBQEHQBha9u2bXrllVd6ln/0ox8pMzNTVqtVq1at6tV21qxZOnfunG655ZaA1OJ0OnX77bdrxIgRmjhxovbu3RuQ/QDwL+66AhAxDMPQP//zP+uXv/xlUPe7a9culZaW6tVXX9Vdd92lEydOBHX/AK4dMzoAQuL48eNKTEzUgQMHJEmffvqpxo4dqz179vS5TmFhoWbPnn3FbeSB9uMf/1grV67U1KlTNWTIENlsNtlstqDWAODaEHQAhMTEiRO1du1aPfzww/ryyy/16KOPqrCwULm5uQHb55IlSzR69Oir/n3zJyS+qbu7Ww0NDTp16pQMw1BSUpKWLVumCxcuBKxOAP5D0AEQMsXFxTIMQzk5OTpx4oSef/75gO5v3bp1+vzzz6/6d+jQoauu8+c//1mXLl1SdXW19u7dqw8++EDvv/++nnvuuYDWCsA/CDoAQqq4uFhHjhzRU089pZiYmFCXc4Xhw4dLkp566imNHz9eN910k77//e9r+/btIa4MgC8IOgBC5ty5c3rmmWdUVFSkVatW6bPPPgvo/p544gnFxcVd9S8jI+Oq6yQkJCgpKUkWi6XntW/+O4DwRtABEDJPP/207Ha7fvWrX2nOnDl64okn+m1/6dIlXbx4UV999ZW6urp08eJFdXd3+7y/X/7ylzp37txV/44ePdrneo8++qh+8Ytf6OTJk/rLX/6if/mXf9G9997r834BhA5BB0BIuFwu1dbW6t/+7d8kSS+88IIOHDigjRs39rlOcXGxhg8frs2bN+v555/X8OHDez1nJ1B+9KMfKTs7W9/+9rd1++2364477tD/+T//J+D7BXD9+K0rAGEpLS1NJ06c0Ny5c1VVVTVg+7q6Os2bN08dHR3avn27ZsyYEYQqAYQ7gg4AADAtTl0BAADTIugAAADTisrfurrpppuUkpIS6jIAAIAftLS06PTp01d9LyqDTkpKihoaGkJdBgAA8AO73d7ne5y6AgAApkXQAQAApkXQAQAApkXQAQAApkXQAQAApkXQAQAApkXQAQAApkXQAQAApkXQAQAApkXQAQAAphXUoFNbW6u0tDQZhqGysrIr3u/o6FB+fr4Mw1BOTo5aWlokSbt27dKUKVOUmZmpKVOm6K233upZZ//+/crMzJRhGFq+fLm8Xm+wugMAAMJc0IJOd3e3li5dqpqaGjU2Nmrz5s1qbGzs1aayslIJCQlqamrSihUrVFpaKunrH+Hctm2bDh8+rKqqKi1cuLBnnSeffFLr16+X2+2W2+1WbW1tsLoEAADCXNCCTn19vQzDUGpqqoYNG6aCggK5XK5ebVwulwoLCyVJ8+fPV11dnbxer+644w7dfPPNkqSMjAxduHBBHR0dOnHihM6ePaupU6fKYrFo0aJF2rJlS7C6BAAAwlzQgo7H41FycnLPclJSkjweT59trFar4uPj1d7e3qvNa6+9pjvvvFMxMTHyeDxKSkrqd5uXVVRUyG63y26369SpU/7qFgAACGPWUBcwGEePHlVpaal27tw56HVLSkpUUlIiqf+fcwcAAOYRtBkdm82m1tbWnuW2tjbZbLY+23R1denMmTMaM2ZMT/u5c+dqw4YNmjhxYk/7tra2frcJAACiV9CCTnZ2ttxut5qbm9XZ2Smn0ymHw9GrjcPhUFVVlSSpurpaM2fOlMVi0eeff645c+aorKxM06ZN62k/fvx4jRo1Svv27ZPX69WGDRt0//33B6tLAAAgzAUt6FitVpWXlysvL0+33367FixYoIyMDK1cuVJbt26VJBUVFam9vV2GYeiFF17ouQW9vLxcTU1N+slPfqKsrCxlZWXp5MmTkqR169bpsccek2EYmjhxombPnh2sLgEAgDBn8Ubhg2fsdrsaGhpCXQYAAPCD/sZ1nowMAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMi6ADAABMK6hBp7a2VmlpaTIMQ2VlZVe839HRofz8fBmGoZycHLW0tEiS2tvbNWPGDMXFxWnZsmW91snNzVVaWpqysrKUlZWlkydPBqMrAAAgAliDtaPu7m4tXbpUu3btUlJSkrKzs+VwOJSent7TprKyUgkJCWpqapLT6VRpaaleffVVxcbGavXq1Tpy5IiOHDlyxbY3btwou90erK4AAIAIEbQZnfr6ehmGodTUVA0bNkwFBQVyuVy92rhcLhUWFkqS5s+fr7q6Onm9Xo0YMULTp09XbGxssMoFAAAmELSg4/F4lJyc3LOclJQkj8fTZxur1ar4+Hi1t7cPuO1HH31UWVlZWr16tbxe71XbVFRUyG63y26369SpU9fREwAAECki/mLkjRs36vDhw9q7d6/27t2rV1555artSkpK1NDQoIaGBo0dOzbIVQIAgFAIWtCx2WxqbW3tWW5ra5PNZuuzTVdXl86cOaMxY8YMuF1JGjlypB588EHV19f7uXIAABCpghZ0srOz5Xa71dzcrM7OTjmdTjkcjl5tHA6HqqqqJEnV1dWaOXOmLBZLn9vs6urS6dOnJUmXLl3S66+/rkmTJgWuEwAAIKIE7a4rq9Wq8vJy5eXlqbu7W4sXL1ZGRoZWrlwpu90uh8OhoqIiLVy4UIZhKDExUU6ns2f9lJQUnT17Vp2dndqyZYt27typW2+9VXl5ebp06ZK6u7s1a9YsFRcXB6tLAAAgzFm8fV29a2J2u10NDQ2hLgMAAPhBf+N6xF+MDAAA0BeCDgAAMK2gXaMDAP6S8uwboS6hR0vZnFCXAKAfgw4658+fV2xsrIYOHRqIegCgT+EUcC67XBOBBwhPAwadr776Sk6nUxs3btR7772nmJgYdXR06KabbtKcOXP0+OOPyzCMYNQKIEqFY8D5awQeIDwNeI3OjBkzdPz4ca1Zs0Z/+tOf1NraqpMnT+qdd97R1KlTVVpaql//+tfBqBVAFIqEkPNNkVYvYHYD3l5+6dIl3XDDDf1uxJc24YTby4HwF+mBgZkdIHiu6/ZyXwJMJIUcAOEv0kOOZI4+AGYw6IuRN23apK1bt2ro0KHyer2677779MADDwSiNgAAgOsy6OfovP322z0XJ2/atEnvvPNOIOoCgIjHrA4QeoMOOh0dHXrjjTd06NAhbd++XRcuXAhEXQCilNnCgdn6A0SaQQeddevW6S9/+Yu2b9+uzz77TOXl5YGoC0AUMmsoMGu/gEjg8zU6Tz/9tH7+85/rxhtv1MMPPxzImgBEIcIAgEDweUZn5MiRcjgcOn/+vCRpx44dmjZtWsAKAwAzIcgBoeHzjM5zzz2nTZs2KTc3V8OGDVNcXJzKysoCWRuAKBEtISDl2Td4vg4QZD4Hnbq6Oq1fv14jRozQiRMn9PLLLystLS2QtQGIAtEScgCEhs+nrp5//nmtXr1ae/bsUXV1tfLz8/XWW28FsjYAAIDr4nPQeeuttzR9+nRJUmZmpmpqavTDH/4wYIUBgBkxgwUE16BvL79s/Pjxqqur82ctAKIMgz6AQBvwGp1PPvnEpw2NHj1ao0aNuu6CAMDsuCgZCJ4Bg05hYaEsFov6+5Fzi8WiRx55RIsWLfJrcQDMi9kcAMEwYNDZvXv3Fa91dHQoJiYmIAUBQDRgVgcIjmu6RmfVqlU9//673/3OX7UAiBLM5gAIlmsKOnl5eaqoqNDGjRtVU1Pj75oAICoQ+IDA8/mBgZcVFRUpPj5e77//vqZOnao1a9YEoi4AJsXg3hunsIDAGnTQqays1IULF3TgwAG99957evzxx/XSSy8FojYAAIDr4vOpq8zMTD300ENau3at9uzZo1tvvVXnz58fVMipra1VWlqaDMO46u9kdXR0KD8/X4ZhKCcnRy0tLZKk9vZ2zZgxQ3FxcVq2bFmvdfbv36/MzEwZhqHly5f3e3cYgNBiNufqOC5A4PgcdN5++20VFxdr+PDhcjqdmjRpkrZv3+7zjrq7u7V06VLV1NSosbFRmzdvVmNjY682lZWVSkhIUFNTk1asWKHS0lJJUmxsrFavXq2f/vSnV2z3ySef1Pr16+V2u+V2u1VbW+tzTQAAwNx8DjqJiYnKzc3V8uXLVVVVpffee0+GYfi8o/r6ehmGodTUVA0bNkwFBQVyuVy92rhcLhUWFkqS5s+fr7q6Onm9Xo0YMULTp09XbGxsr/YnTpzQ2bNnNXXqVFksFi1atEhbtmzxuSYAwcOsBYBQ8DnoHDt2rNfyt771LR06dMjnHXk8HiUnJ/csJyUlyePx9NnGarUqPj5e7e3t/W4zKSmp321eVlFRIbvdLrvdrlOnTvlcNwAAiFw+X4z8+OOP6/jx47LZbJo8ebIuXryoSZMm6csvv9SNN94YyBr9oqSkRCUlJZIku90e4moAoDfuvgICw+egc/kJyZ988okOHjyoDz74QAcPHlRWVpaGDBmiDz/8sN/1bTabWltbe5bb2tpks9mu2iYpKUldXV06c+aMxowZ0+8229ra+t0mgNDjtBWAUBkw6Hi9Xlkslp7lW265Rbfccovuu+++ntfOnj074I6ys7PldrvV3Nwsm80mp9OpTZs29WrjcDhUVVWl7373u6qurtbMmTN77fuvjR8/XqNGjdK+ffuUk5OjDRs26KmnnhqwFgAIR8zqAP434DU6M2bM0C9+8YsrfsW8s7NTb731lgoLC/Vf//VfA+7IarWqvLxceXl5uv3227VgwQJlZGRo5cqV2rp1q6SvH0bY3t4uwzD0wgsv9LoFPSUlRd///vf1H//xH0pKSuq5Y2vdunV67LHHZBiGJk6cqNmzZw/qAAAILGZzAISSxTvAg2cuXryol19+WRs3btRHH32khIQEXbhwQV999ZXuueceLVmyRHfccUew6vULu92uhoaGUJcBRAWCzuAxqwMMTn/j+oCnrmJjY7VkyRItWbJEly5d0unTpzV8+HCNHj3a33UCMBlCDoBQ8/n28pqaGt19993Kzc1VSUmJ9u3bF8i6ACBqERAB//E56CxZskQ/+9nPtG/fPpWUlOgHP/iBNm/eHMjaAEQwBmsA4cDnoDNu3DhNmzZNCQkJmjVrlnbs2KHnn38+kLUBQNQiKAL+4XPQmTBhgn74wx+qs7NTknTDDTfIah30j58DiAIM0gDChc9BZ8iQIfrNb36j5ORkTZ8+XYZhKDc3V263O5D1AUDUIjAC18/nKZnLD/fr6OjQkSNHdPDgQR08eFDFxcX66KOPrnjODgAAQKgN+txTTEyMpkyZoilTpgSiHgARjlkIAOHE51NXAAAAkYagAwBhjBky4PoQdAAAgGkRdAD4DbMPgcFxBa4dQQcAAJgWQQeAXzDrEFgcX+DaEHQAAIBpEXQAXDdmG4KD4wwMHkEHAACYFkEHwHVhliG4ON7A4BB0AACAaRF0AFwzZhcAhDuCDgAAMC2CDgBEGGbSAN8RdABcEwZbAJGAoAMAEYigCfiGoANg0BhkAUQKgg4ARCgCJzCwoAad2tpapaWlyTAMlZWVXfF+R0eH8vPzZRiGcnJy1NLS0vPemjVrZBiG0tLStGPHjp7XU1JSlJmZqaysLNnt9mB0A4hqDK4AIknQgk53d7eWLl2qmpoaNTY2avPmzWpsbOzVprKyUgkJCWpqatKKFStUWloqSWpsbJTT6dTRo0dVW1urJUuWqLu7u2e93bt364MPPlBDQ0OwugMAYYHgCfQvaEGnvr5ehmEoNTVVw4YNU0FBgVwuV682LpdLhYWFkqT58+errq5OXq9XLpdLBQUFiomJ0YQJE2QYhurr64NVOoD/h0E1PPG5AH0LWtDxeDxKTk7uWU5KSpLH4+mzjdVqVXx8vNrb2/td12Kx6J577tGUKVNUUVHR5/4rKipkt9tlt9t16tQpf3YNAACEqYi/GPmdd97RgQMHVFNToxdffFG//e1vr9qupKREDQ0Namho0NixY4NcJRD5mDUIb3w+wNUFLejYbDa1trb2LLe1tclms/XZpqurS2fOnNGYMWP6XffyP8eNG6e5c+dySgsIAAZRAJEqaEEnOztbbrdbzc3N6uzslNPplMPh6NXG4XCoqqpKklRdXa2ZM2fKYrHI4XDI6XSqo6NDzc3Ncrvduuuuu3T+/Hl98cUXkqTz589r586dmjRpUrC6BAAAwpw1aDuyWlVeXq68vDx1d3dr8eLFysjI0MqVK2W32+VwOFRUVKSFCxfKMAwlJibK6XRKkjIyMrRgwQKlp6fLarXqxRdf1NChQ/XnP/9Zc+fOlfT1DNCDDz6o733ve8HqEgCElZRn31BL2ZxQlwGEFYvX6/WGuohgs9vt3IoODAKnriIHQQfRqL9xPeIvRgYQWIScyMLnBfRG0AHQJwZNAJGOoAMAJkNABf4bQQfAVTFYRjY+P+BrBB0AV2CQBGAWBB0AvRByzIPPEiDoAAAAEyPoAICJMauDaEfQAdCDQRGA2RB0AEgi5JgZny2iGUEHAKIAYQfRiqADgEEwSvA5IxoRdIAox+AXXfi8EW0IOkAUY9CLTnzuiCYEHSBKMdgBiAYEHQCIQgRdRAuCDhCFGOQg8T1AdCDoAFGGwQ3fxPcBZkfQAaJEyrNvMKjhqvhewMwIOkAUYCDDQPiOwKwIOgAASYQdmJM11AUACBwGLgzW5e9MS9mcEFcC+AczOoAJcT0OrhffH5gFMzqAiTA4wZ+Y3YEZMKMDmAQhB4HCDCEiGTM6QARj8EEwMcODSBTUGZ3a2lqlpaXJMAyVlZVd8X5HR4fy8/NlGIZycnLU0tLS896aNWtkGIbS0tK0Y8cOn7cJmBH/DxuhdPn7x3cQkSBoMzrd3d1aunSpdu3apaSkJGVnZ8vhcCg9Pb2nTWVlpRISEtTU1CSn06nS0lK9+uqramxslNPp1NGjR/Xpp59q1qxZOnbsmCQNuE0g0jGYIJz99feT2R6Em6AFnfr6ehmGodTUVElSQUGBXC5Xr1Dicrm0atUqSdL8+fO1bNkyeb1euVwuFRQUKCYmRhMmTJBhGKqvr5ekAbcJRAoCDcygr+8xAQihErSg4/F4lJyc3LOclJSkd999t882VqtV8fHxam9vl8fj0dSpU3ut6/F4JGnAbV5WUVGhiooKSdKHH34ou93un4752alTpzR27NhQl2EqkXJMbwp1AYMQKcc0kpj9mNrtPw7q/sx+PEMhnI/pNy91+WtRczFySUmJSkpKQl3GgOx2uxoaGkJdhqlwTP2PY+p/HFP/4nj6X6Qe06BdjGyz2dTa2tqz3NbWJpvN1mebrq4unTlzRmPGjOlzXV+2CQAAolfQgk52drbcbream5vV2dkpp9Mph8PRq43D4VBVVZUkqbq6WjNnzpTFYpHD4ZDT6VRHR4eam5vldrt11113+bRNAAAQvYJ26spqtaq8vFx5eXnq7u7W4sWLlZGRoZUrV8put8vhcKioqEgLFy6UYRhKTEyU0+mUJGVkZGjBggVKT0+X1WrViy++qKFDh0rSVbcZySLh9Fqk4Zj6H8fU/zim/sXx9L9IPaYWr9frDXURAAAAgcBPQAAAANMi6AAwhdzcXH344YeSpPb2dk2aNCnEFQEIBwQdAKbQ1NSkb3/725KkQ4cOKTMzM8QVAQgHBB0AEe/jjz+WzWbTkCFf/yft0KFDmjx5siRp8eLFoSwNQIgRdABEvIMHD/YEG0nav3+/Jk+erC+//FLx8fHavXu3SktLdfHixRBWCSAUCDoAIt4HH3zQE2LcbrdcLpcyMzN14MABvf/++/rjH/+otWvXKjY2NsSVAgi2qPkJCADmdfDgQcXGxuo73/mOJk+erPT0dFVVVSkuLk45OTkaMWJEqEsEECIEHQAR79ChQzpw4IBGjhzZ6/Xi4mK99NJL+tGPfqS9e/fq7rvvDlGFAEKFBwYCiGhffPGFpkyZomPHjoW6FABhiKADAABMi4uRAQCAaRF0AISllJQUDR8+XAsXLvSp/Ztvvqm4uDgNGTJEb775ZoCrAxApCDoAwta2bdv0yiuvSJJOnjypBx54QDfffLPi4+M1bdo0vfvuuz1tZ82apXPnzumWW27xex1xcXG9/oYOHaqnnnrK7/sB4H8EHQAR4dy5c8rOztb+/fv12WefqbCwUHPmzNG5c+eCsu/Lf3/60580fPhw/f3f/33A9wvg+hF0AITE8ePHlZiYqAMHDkiSPv30U40dO1Z79uy5avvU1FR9//vf1/jx4zV06FCVlJSos7NTf/zjH4NYtfTaa69p3Lhx3KoORAiCDoCQmDhxotauXauHH35YX375pR599FEVFhYqNzfXp/U/+OADdXZ2yjAMn/e5ZMkSjR49+qp/3/wJif5UVVVp0aJFslgsPu8XQOjwwEAAIVNcXKxt27YpJydHFotFW7du9Wm9s2fPauHChfrxj3+s+Ph4n/e3bt06rVu37lrL1ccff6y3335blZWV17wNAMHFjA6AkCouLtaRI0f01FNPKSYmZsD2Fy5c0H333aepU6fqH//xH4NQ4X975ZVXNH36dE2YMCGo+wVw7Qg6AELm3LlzeuaZZ1RUVKRVq1bps88+67d9R0eH/tf/+l9KSkrSSy+9NOj9PfHEE1fcQXX5LyMjY8D1N2zYoMLCwkHvF0DoEHQAhMzTTz8tu92uX/3qV5ozZ46eeOKJPtteunRJ8+fP1/Dhw1VVVaUhQwb/n69f/vKXve6g+ubf0aNH+13397//vTweD3dbARGGoAMgJFwul2pra/Vv//ZvkqQXXnhBBw4c0MaNG6/a/ve//71ef/117dy5U6NHj+6Zidm7d29Q6q2qqtLf/d3fXfHDoQDCG791BSAspaWl6cSJE5o7d66qqqoGbF9XV6d58+apo6ND27dv14wZM4JQJYBwR9ABAACmxakrAABgWgQdAABgWlH5wMCbbrpJKSkpoS4DAAD4QUtLi06fPn3V96Iy6KSkpKihoSHUZQAAAD+w2+19vsepKwAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFpBDTq1tbVKS0uTYRgqKyu74v2Ojg7l5+fLMAzl5OSopaVFkrRr1y5NmTJFmZmZmjJlit56662edfbv36/MzEwZhqHly5fL6/UGqzsAACDMBS3odHd3a+nSpaqpqVFjY6M2b96sxsbGXm0qKyuVkJCgpqYmrVixQqWlpZK+/hHObdu26fDhw6qqqtLChQt71nnyySe1fv16ud1uud1u1dbWBqtLAAAgzAUt6NTX18swDKWmpmrYsGEqKCiQy+Xq1cblcqmwsFCSNH/+fNXV1cnr9eqOO+7QzTffLEnKyMjQhQsX1NHRoRMnTujs2bOaOnWqLBaLFi1apC1btgSrSwAAIMwFLeh4PB4lJyf3LCclJcnj8fTZxmq1Kj4+Xu3t7b3avPbaa7rzzjsVExMjj8ejpKSkfrcJAACilzXUBQzG0aNHVVpaqp07dw563YqKClVUVEiSTp065e/SAABAGArajI7NZlNra2vPcltbm2w2W59turq6dObMGY0ZM6an/dy5c7VhwwZNnDixp31bW1u/27yspKREDQ0Namho0NixY/3aNwAAEJ6CFnSys7PldrvV3Nyszs5OOZ1OORyOXm0cDoeqqqokSdXV1Zo5c6YsFos+//xzzZkzR2VlZZo2bVpP+/Hjx2vUqFHat2+fvF6vNmzYoPvvvz9YXQIAAGEuaEHHarWqvLxceXl5uv3227VgwQJlZGRo5cqV2rp1qySpqKhI7e3tMgxDL7zwQs8t6OXl5WpqatJPfvITZWVlKSsrSydPnpQkrVu3To899pgMw9DEiRM1e/bsYHUJAACEOYs3Ch88Y7fb1dDQEOoyAACAH/Q3rvNkZAAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFoEHQAAYFrWUBcAAKGS8uwboS7hCi1lc0JdAmAqBB0AUSccA85ll2sj8AD+wakrAFElnEPON0VKnUC4C2rQqa2tVVpamgzDUFlZ2RXvd3R0KD8/X4ZhKCcnRy0tLZKk9vZ2zZgxQ3FxcVq2bFmvdXJzc5WWlqasrCxlZWXp5MmTwegKgAgUaeEh0uoFwlHQTl11d3dr6dKl2rVrl5KSkpSdnS2Hw6H09PSeNpWVlUpISFBTU5OcTqdKS0v16quvKjY2VqtXr9aRI0d05MiRK7a9ceNG2e32YHUFAABEiKDN6NTX18swDKWmpmrYsGEqKCiQy+Xq1cblcqmwsFCSNH/+fNXV1cnr9WrEiBGaPn26YmNjg1UuAIQFZnWA6xO0oOPxeJScnNyznJSUJI/H02cbq9Wq+Ph4tbe3D7jtRx99VFlZWVq9erW8Xq9/CwdgCpEcGCK5diDUIv5i5I0bN+rw4cPau3ev9u7dq1deeeWq7SoqKmS322W323Xq1KkgVwkglAgKQPQKWtCx2WxqbW3tWW5ra5PNZuuzTVdXl86cOaMxY8YMuF1JGjlypB588EHV19dftV1JSYkaGhrU0NCgsWPHXk9XACDoCGvAtQla0MnOzpbb7VZzc7M6OzvldDrlcDh6tXE4HKqqqpIkVVdXa+bMmbJYLH1us6urS6dPn5YkXbp0Sa+//romTZoUuE4AiDgEBCC6Be2uK6vVqvLycuXl5am7u1uLFy9WRkaGVq5cKbvdLofDoaKiIi1cuFCGYSgxMVFOp7Nn/ZSUFJ09e1adnZ3asmWLdu7cqVtvvVV5eXm6dOmSuru7NWvWLBUXFwerSwAQVCnPvsGDBIFBsnij8Opdu92uhoaGUJcBIMDMOJtD0AGu1N+4HvEXIwNANDFjeAMCiaADwJQIBAAkgg4AADAxgg4AADCtQQed8+fPq7u7OxC1AAB8wGk5wHcDBp2vvvpKmzZt0pw5czRu3DjddtttGj9+vNLT0/UP//APampqCkadAOAzggCAywYMOjNmzNDx48e1Zs0a/elPf1Jra6tOnjypd955R1OnTlVpaal+/etfB6NWAMD/Q5gDfDPgAwPffPNN3XDDDVe8npiYqHnz5mnevHm6dOlSQIoDgMEiAAD4pgFndK4Wcq6lDQDAvwh1wMAG/RMQmzZt0tatWzV06FB5vV7dd999euCBBwJRGwAMCgM/gL826Luu3n77bTmdTm3cuFGbNm3SO++8E4i6AAA+INwB/Rv0jE5HR4feeOMNJScnq62tTRcuXAhEXQAAANdt0DM669at01/+8hdt375dn332mcrLywNRFwAMCjMbAK7G5xmdp59+Wj//+c9144036uGHHw5kTQCAQUh59g1+1Rzog88zOiNHjpTD4dD58+clSTt27NC0adMCVhgA+IrZHAB98XlG57nnntOmTZuUm5urYcOGKS4uTmVlZYGsDQDgI2Z1gKvzOejU1dVp/fr1GjFihE6cOKGXX35ZaWlpgawNAADguvh86ur555/X6tWrtWfPHlVXVys/P19vvfVWIGsDgAFx2gpAf3ye0flmqMnMzFRNTY3mzZun3//+9wEpDAAA4HoN+vbyy8aPH6+6ujp/1gIAAOBXA87ofPLJJz5taPTo0Ro1atR1FwQAuDZckAxcacCgU1hYKIvFIq/X22cbi8WiRx55RIsWLfJrcQAAANdjwKCze/fuK17r6OhQTExMQAoCAF9xITKAgVzTNTqrVq3q+fff/e53/qoFAHCdCH9Ab4P+UU9JysvLU0VFhUaMGKE//OEPPCEZQNAxoAPwxaCDTlFRkeLj4/X+++9r6tSpWrNmTSDqAgAAuG6DPnVVWVmp559/Xs8995z+5m/+Ro8//rjP69bW1iotLU2GYVz15yM6OjqUn58vwzCUk5OjlpYWSVJ7e7tmzJihuLg4LVu2rNc6+/fvV2ZmpgzD0PLly/u9aBoAogGzXcB/8znoZGZm6qGHHtLatWu1Z88e3XrrrTp//rxeeukln9bv7u7W0qVLVVNTo8bGRm3evFmNjY292lRWViohIUFNTU1asWKFSktLJUmxsbFavXq1fvrTn16x3SeffFLr16+X2+2W2+1WbW2tr10CEKEYyAH4yueg8/bbb6u4uFjDhw+X0+nUpEmTtH37dp93VF9fL8MwlJqaqmHDhqmgoEAul6tXG5fLpcLCQknS/PnzVVdXJ6/XqxEjRmj69OmKjY3t1f7EiRM6e/aspk6dKovFokWLFmnLli0+1wQAZkUYBL7m8zU6iYmJys3NVW5uriTJ7Xbrueee83lHHo9HycnJPctJSUl69913+2xjtVoVHx+v9vZ23XTTTX1uMykpqdc2PR6PzzUBAABz83lG59ixY72Wv/Wtb+nQoUN+LyhQKioqZLfbZbfbderUqVCXA+AaMVMBYDB8ntF5/PHHdfz4cdlsNk2ePFkXL17UpEmT9OWXX+rGG28ccH2bzabW1tae5ba2Ntlstqu2SUpKUldXl86cOaMxY8b0u822trZ+t3lZSUmJSkpKJEl2u33AegEg0vGTEMAgZnR2796tTz75RK+++qruvfdeGYahCxcuKCsrS7fddtuA62dnZ8vtdqu5uVmdnZ1yOp1yOBy92jgcDlVVVUmSqqurNXPmTFkslj63OX78eI0aNUr79u2T1+vVhg0bdP/99/vaJQAAYHIDzuh4vd5eYeOWW27RLbfcovvuu6/ntbNnzw68I6tV5eXlysvLU3d3txYvXqyMjAytXLlSdrtdDodDRUVFWrhwoQzDUGJiopxOZ8/6KSkpOnv2rDo7O7Vlyxbt3LlT6enpWrdunR555BFduHBBs2fP1uzZswd7DABECE5bARgsi3eAB8/k5uZq3rx5uv/++3XLLbf0vN7Z2al33nlHVVVVmjFjhh555JFA1+o3drtdDQ0NoS4DwCARdAaPU1eIBv2N6wPO6NTW1urll1/WAw88oI8++kgJCQm6cOGCvvrqK91zzz165plndMcdd/i9aAAAgOs1YNCJjY3VkiVLtGTJEl26dEmnT5/W8OHDNXr06CCUBwAAcO18vhi5pqZGd999t3Jzc1VSUqJ9+/YFsi4AgB9wug/Rzuegs2TJEv3sZz/Tvn37VFJSoh/84AfavHlzIGsDgB4M2ACuhc9BZ9y4cZo2bZoSEhI0a9Ys7dixQ88//3wgawMA+AEhEdHM56AzYcIE/fCHP1RnZ6ck6YYbbpDV6vPzBgEAAILO56AzZMgQ/eY3v1FycrKmT58uwzCUm5srt9sdyPoAgBkJANfM5ymZTZs2SZI6Ojp05MgRHTx4UAcPHlRxcbE++ugjffLJJwErEgBwffg5CESrQZ97iomJ0ZQpUzRlypRA1AMAvTCbA+B6+HzqCgAAINIQdAAgSjA7hmhE0AEQthiYAVwvgg4ARBHCI6INQQdAWGJABuAPBB0AAGBaBB0AiDLMliGaEHQAhB0GYgD+QtABAACmRdABAACmRdABAACmRdABgCjEdVCIFgQdAGGFARiAPxF0ACBKESoRDQg6AMIGAy8AfyPoAEAUI1zC7Ag6AADAtIIadGpra5WWlibDMFRWVnbF+x0dHcrPz5dhGMrJyVFLS0vPe2vWrJFhGEpLS9OOHTt6Xk9JSVFmZqaysrJkt9uD0Q0AAcDMAoBACFrQ6e7u1tKlS1VTU6PGxkZt3rxZjY2NvdpUVlYqISFBTU1NWrFihUpLSyVJjY2NcjqdOnr0qGpra7VkyRJ1d3f3rLd792598MEHamhoCFZ3AMA0CJkws6AFnfr6ehmGodTUVA0bNkwFBQVyuVy92rhcLhUWFkqS5s+fr7q6Onm9XrlcLhUUFCgmJkYTJkyQYRiqr68PVukAAoyBFkCgBC3oeDweJScn9ywnJSXJ4/H02cZqtSo+Pl7t7e39rmuxWHTPPfdoypQpqqioCEJPAMB8CJswK2uoC7he77zzjmw2m06ePKn/+T//p2677Tb9j//xP65oV1FR0ROETp06FewyAfSBARZAIAVtRsdms6m1tbVnua2tTTabrc82XV1dOnPmjMaMGdPvupf/OW7cOM2dO7fPU1olJSVqaGhQQ0ODxo4d69e+AQCA8BS0oJOdnS23263m5mZ1dnbK6XTK4XD0auNwOFRVVSVJqq6u1syZM2WxWORwOOR0OtXR0aHm5ma53W7dddddOn/+vL744gtJ0vnz57Vz505NmjQpWF0CAABhLminrqxWq8rLy5WXl6fu7m4tXrxYGRkZWrlypex2uxwOh4qKirRw4UIZhqHExEQ5nU5JUkZGhhYsWKD09HRZrVa9+OKLGjp0qP785z9r7ty5kr6eAXrwwQf1ve99L1hdAgBTSXn2DbWUzQl1GYBfWbxerzfURQSb3W7nVnQgDHB9Tvgh6CAS9Teu82RkAEAPwifMhqADICQYUAEEA0EHANALIRRmQtABEHQMpACChaADALgCYRRmQdABEFQMoJGDzwpmQNABEDQMnACCjaADAOgT4RSRjqADICgYMAGEAkEHAACYFkEHANAvZuMQyQg6AAKOgRJAqBB0AAQUIccc+BwRqQg6AACfEHYQiQg6AAKGgdF8+EwRaQg6AAKCAdG8+GwRSQg6APyOgRBAuCDoAPArQg6AcELQAQAMGoEWkYKgA8BvGPyiC583IgFBB4BfMOhFp5Rn3+CzR1gj6AC4bgx04DuAcEXQAXDN+H/z+Ca+CwhHBB0AgN8QdhBurKEuAEDkYTBDfy5/P1rK5oS4EoCgA2AQCDgYDAIPwkFQT13V1tYqLS1NhmGorKzsivc7OjqUn58vwzCUk5OjlpaWnvfWrFkjwzCUlpamHTt2+LxNANePa3FwPfjuIJSCNqPT3d2tpUuXateuXUpKSlJ2drYcDofS09N72lRWViohIUFNTU1yOp0qLS3Vq6++qsbGRjmdTh09elSffvqpZs2apWPHjknSgNsEcG0YnOBPf/19YpYHwRK0oFNfXy/DMJSamipJKigokMvl6hVKXC6XVq1aJUmaP3++li1bJq/XK5fLpYKCAsXExGjChAkyDEP19fWSNOA2AQyMUINgI/ggWIIWdDwej5KTk3uWk5KS9O677/bZxmq1Kj4+Xu3t7fJ4PJo6dWqvdT0ejyQNuE0gGhFcEGmu5TtLOIIvouZi5IqKClVUVEiSPvzwQ9nt9hBXdHWnTp3S2LFjQ12GqUTjMb0pwNuPxmMaaBzTwbPbf9znexxP/wvnY/rNa3r/WtCCjs1mU2tra89yW1ubbDbbVdskJSWpq6tLZ86c0ZgxY/pdd6BtXlZSUqKSkhJ/dikg7Ha7GhoaQl2GqXBM/Y9j6n8cU//iePpfpB7ToN11lZ2dLbfbrebmZnV2dsrpdMrhcPRq43A4VFVVJUmqrq7WzJkzZbFY5HA45HQ61dHRoebmZrndbt11110+bRMAAESvoM3oWK1WlZeXKy8vT93d3Vq8eLEyMjK0cuVK2e12ORwOFRUVaeHChTIMQ4mJiXI6nZKkjIwMLViwQOnp6bJarXrxxRc1dOhQSbrqNgEAACTJ4vV6vaEuAv+toqIiIk6xRRKOqf9xTP2PY+pfHE//i9RjStABAACmxY96AjCF3Nxcffjhh5Kk9vZ2TZo0KcQVAQgHBB0AptDU1KRvf/vbkqRDhw4pMzMzxBUBCAcEHQAR7+OPP5bNZtOQIV//J+3QoUOaPHmyJGnx4sWhLA1AiBF0AES8gwcP9gQbSdq/f78mT56sL7/8UvHx8dq9e7dKS0t18eLFEFYJIBQIOgAi3gcffNATYtxut1wulzIzM3XgwAG9//77+uMf/6i1a9cqNjY2xJUCCLao+QkIAOZ18OBBxcbG6jvf+Y4mT56s9PR0VVVVKS4uTjk5ORoxYkSoSwQQIgQdABHv0KFDOnDggEaOHNnr9eLiYr300kv60Y9+pL179+ruu+8OUYUAQoXn6ACIaF988YWmTJmiY8eOhboUAGGIoAMAAEyLi5EBAIBpEXQAhKWUlBQNHz5cCxcu9Kn9m2++qbi4OA0ZMkRvvvlmgKsDECkIOgDC1rZt2/TKK6/0LM+YMUNjx47VqFGj9J3vfEcul6vnvVmzZuncuXO65ZZb/F5HXFxcr7+hQ4fqqaee8vt+APgfd10BiBj/+q//qvT0dFmtVr377ruaNWuWjh07pvHjxwd0v+fOnev173/7t3+rv//7vw/oPgH4BzM6AELi+PHjSkxM1IEDByRJn376qcaOHas9e/b0uc7kyZNltX79/88sFosuXbqk1tbWYJTb47XXXtO4ceO4VR2IEAQdACExceJErV27Vg8//LC+/PJLPfrooyosLFRubm6/6917772KjY1VTk6OcnNzZbfbfd7nkiVLNHr06Kv+ffMnJPpTVVWlRYsWyWKx+LxfAKHDqSsAIVNcXKxt27YpJydHFotFW7duHXCd119/XZcuXdKbb76pP/zhDz0/5OmLdevWad26dddc78cff6y3335blZWV17wNAMHFjA6AkCouLtaRI0f01FNPKSYmxqd1brjhBs2ePVs7d+70KRz5yyuvvKLp06drwoQJQdsngOtD0AEQMufOndMzzzyjoqIirVq1Sp999tmg1u/q6tLx48d9bv/EE09ccQfV5b+MjIwB19+wYYMKCwsHVSOA0CLoAAiZp59+Wna7Xb/61a80Z84cPfHEE322/fDDD1VTU6MLFy7o0qVL+vWvf63f/va3+v/+v//P5/398pe/1Llz5676d/To0X7X/f3vfy+Px8PdVkCE4RodACHhcrlUW1urw4cPS5JeeOEFZWVlaePGjXrooYeuaO/1erVq1So1NjZq6NCh+ta3vqVXX31Vd955Z1Dqraqq0t/93d9d8cOhAMIbv3UFICylpaXpxIkTmjt3rqqqqgZsX1dXp3nz5qmjo0Pbt2/XjBkzglAlgHBH0AEAAKbFNToAAMC0CDoAAMC0CDoAAMC0ovKuq5tuukkpKSmhLgMAAPhBS0uLTp8+fdX3ojLopKSkqKGhIdRlAAAAP+jvN+84dQUAAEyLoAMAAEyLoAMAAEyLoAMAAEyLoAMAAEyLoAMAAEyLoAMAAEyLoAMAAEyLoAMAAEyLoAMAAEwrqEGntrZWaWlpMgxDZWVlV7zf0dGh/Px8GYahnJwctbS0SJJ27dqlKVOmKDMzU1OmTNFbb73Vs87+/fuVmZkpwzC0fPlyeb3eYHUHAACEuaAFne7ubi1dulQ1NTVqbGzU5s2b1djY2KtNZWWlEhIS1NTUpBUrVqi0tFTS1z/CuW3bNh0+fFhVVVVauHBhzzpPPvmk1q9fL7fbLbfbrdra2mB1CQAAhLmgBZ36+noZhqHU1FQNGzZMBQUFcrlcvdq4XC4VFhZKkubPn6+6ujp5vV7dcccduvnmmyVJGRkZunDhgjo6OnTixAmdPXtWU6dOlcVi0aJFi7Rly5ZgdQkAAIS5oAUdj8ej5OTknuWkpCR5PJ4+21itVsXHx6u9vb1Xm9dee0133nmnYmJi5PF4lJSU1O82L6uoqJDdbpfdbtepU6f81S0AABDGrKEuYDCOHj2q0tJS7dy5c9DrlpSUqKSkRFL/P+cOAADMI2gzOjabTa2trT3LbW1tstlsfbbp6urSmTNnNGbMmJ72c+fO1YYNGzRx4sSe9m1tbf1uEwAARK+gBZ3s7Gy53W41Nzers7NTTqdTDoejVxuHw6GqqipJUnV1tWbOnCmLxaLPP/9cc+bMUVlZmaZNm9bTfvz48Ro1apT27dsnr9erDRs26P777w9WlwAAQJgLWtCxWq0qLy9XXl6ebr/9di1YsEAZGRlauXKltm7dKkkqKipSe3u7DMPQCy+80HMLenl5uZqamvSTn/xEWVlZysrK0smTJyVJ69at02OPPSbDMDRx4kTNnj07WF0CAABhzuKNwgfP2O12NTQ0hLoMAADgB/2N6zwZGQAAmBZBBwAAmBZBBwAAmBZBBwAAmBZBBwAAmBZBBwAAmBZBBwAAmBZBBwAAmBZBBwAAmFZE/Xo5AIRayrNvhLqEHi1lc0JdAhD2CDoA4INwCjiXXa6JwAP0jVNXADCAcAw53xTu9QGhRNABgH5ESoiIlDqBYCPoAEAfCA9A5CPoAAAA0yLoAMBVROJsTiTWDAQaQQcAAJgWQQcATIRZHaA3gg4AADAtgg4A/JVInxWJ9PoBfyLoAAAA0yLoAMA3mGU2xCz9AK4XQQcAAJgWQQcA/h+zzYKYrT/AtSDoAAAA0yLoAAAA0wpq0KmtrVVaWpoMw1BZWdkV73d0dCg/P1+GYSgnJ0ctLS2SpPb2ds2YMUNxcXFatmxZr3Vyc3OVlpamrKwsZWVl6eTJk8HoCgCT4TQPYE7WYO2ou7tbS5cu1a5du5SUlKTs7Gw5HA6lp6f3tKmsrFRCQoKamprkdDpVWlqqV199VbGxsVq9erWOHDmiI0eOXLHtjRs3ym63B6srAAAgQgRtRqe+vl6GYSg1NVXDhg1TQUGBXC5XrzYul0uFhYWSpPnz56uurk5er1cjRozQ9OnTFRsbG6xyAUQRM8/mmLlvgC+CFnQ8Ho+Sk5N7lpOSkuTxePpsY7VaFR8fr/b29gG3/eijjyorK0urV6+W1+u9apuKigrZ7XbZ7XadOnXqOnoCAAAiRcRfjLxx40YdPnxYe/fu1d69e/XKK69ctV1JSYkaGhrU0NCgsWPHBrlKAAAQCkELOjabTa2trT3LbW1tstlsfbbp6urSmTNnNGbMmAG3K0kjR47Ugw8+qPr6ej9XDgCRjdNXiGZBCzrZ2dlyu91qbm5WZ2ennE6nHA5HrzYOh0NVVVWSpOrqas2cOVMWi6XPbXZ1den06dOSpEuXLun111/XpEmTAtcJAAAQUYJ215XValV5ebny8vLU3d2txYsXKyMjQytXrpTdbpfD4VBRUZEWLlwowzCUmJgop9PZs35KSorOnj2rzs5ObdmyRTt37tStt96qvLw8Xbp0Sd3d3Zo1a5aKi4uD1SUAABDmLN6+rt41MbvdroaGhlCXASAMRNNpnZayOaEuAQiI/sb1iL8YGQAAoC8EHQAAYFoEHQBRK5pOW0nR119Auoagc/78eXV3dweiFgAAAL8aMOh89dVX2rRpk+bMmaNx48bptttu0/jx45Wenq5/+Id/UFNTUzDqBAAAGLQBg86MGTN0/PhxrVmzRn/605/U2tqqkydP6p133tHUqVNVWlqqX//618GoFQD8JlpP40RrvxG9BnyOzptvvqkbbrjhitcTExM1b948zZs3T5cuXQpIcQAAANdjwBmdq4Wca2kDAAAQbIN+MvKmTZu0detWDR06VF6vV/fdd58eeOCBQNQGAABwXQZ919Xbb78tp9OpjRs3atOmTXrnnXcCURcABEy0X6cS7f1HdBn0jE5HR4feeOMNJScnq62tTRcuXAhEXQAAANdt0DM669at01/+8hdt375dn332mcrLywNRFwAAwHXzeUbn6aef1s9//nPdeOONevjhhwNZEwAEDKdtgOji84zOyJEj5XA4dP78eUnSjh07NG3atIAVBgAAcL18ntF57rnntGnTJuXm5mrYsGGKi4tTWVlZIGsDAAC4Lj7P6NTV1Wn9+vUaMWKETp8+rf/7f/+v7r777kDWBgB+xWmr/8axQLTwOeg8//zzWr16tfbs2aPq6mrl5+frrbfeCmRtAAAA18XnU1ffDDWZmZmqqanRvHnz9Pvf/z4ghQEAAFyvQd9eftn48eNVV1fnz1oAAAD8asAZnU8++cSnDY0ePVqjRo267oIAAMGR8uwbaimbE+oygIAaMOgUFhbKYrHI6/X22cZiseiRRx7RokWL/FocAADA9Rgw6OzevfuK1zo6OhQTExOQggAAAPzlmq7RWbVqVc+//+53v/NXLQAQMNxOfXUcF5jdoH/UU5Ly8vJUUVGhESNG6A9/+ANPSAYAAGFp0EGnqKhI8fHxev/99zV16lStWbMmEHUBAABcN4u3v6uM+3DhwgUdOHBA7733nv7whz/opZdeCkRtAWO329XQ0BDqMgAECadnBsbdV4hk/Y3rPl+jk5mZqYceekhr167Vnj17dOutt+r8+fODCjm1tbVKS0uTYRhX/Z2sjo4O5efnyzAM5eTkqKWlRZLU3t6uGTNmKC4uTsuWLeu1zv79+5WZmSnDMLR8+fJ+7w4DAADRxeeg8/bbb6u4uFjDhw+X0+nUpEmTtH37dp931N3draVLl6qmpkaNjY3avHmzGhsbe7WprKxUQkKCmpqatGLFCpWWlkqSYmNjtXr1av30pz+9YrtPPvmk1q9fL7fbLbfbrdraWp9rAgAA5uZz0ElMTFRubq6WL1+uqqoqvffeezIMw+cd1dfXyzAMpaamatiwYSooKJDL5erVxuVyqbCwUJI0f/581dXVyev1asSIEZo+fbpiY2N7tT9x4oTOnj2rqVOnymKxaNGiRdqyZYvPNQEAAHPzOegcO3as1/K3vvUtHTp0yOcdeTweJScn9ywnJSXJ4/H02cZqtSo+Pl7t7e39bjMpKanfbV5WUVEhu90uu92uU6dO+Vw3gMjG9Tm+4TjBrHy+6+rxxx/X8ePHZbPZNHnyZF28eFGTJk3Sl19+qRtvvDGQNfpFSUmJSkpKJH190RIAADA/n4PO5Sckf/LJJzp48KA++OADHTx4UFlZWRoyZIg+/PDDfte32WxqbW3tWW5ra5PNZrtqm6SkJHV1denMmTMaM2ZMv9tsa2vrd5sAACB6DRh0vF6vLBZLz/Itt9yiW265Rffdd1/Pa2fPnh1wR9nZ2XK73WpubpbNZpPT6dSmTZt6tXE4HKqqqtJ3v/tdVVdXa+bMmb32/dfGjx+vUaNGad++fcrJydGGDRv01FNPDVgLgOjA6RgAA16jM2PGDP3iF7+44lfMOzs79dZbb6mwsFD/9V//NeCOrFarysvLlZeXp9tvv10LFixQRkaGVq5cqa1bt0r6+mGE7e3tMgxDL7zwQq9b0FNSUvT9739f//Ef/6GkpKSeO7bWrVunxx57TIZhaOLEiZo9e/agDgAAADCvAR8YePHiRb388svauHGjPvroIyUkJOjChQv66quvdM8992jJkiW64447glWvX/DAQCA6MKMzODw0EJGqv3F9UE9GvnTpkk6fPq3hw4dr9OjR/qov6Ag6gPkRcq4NYQeRyC9PRq6pqdHdd9+t3NxclZSUaN++fX4rEAAAIBB8DjpLlizRz372M+3bt08lJSX6wQ9+oM2bNweyNgAAgOvi8+3l48aN07Rp0yRJs2bN0ne/+13l5OTogQceCFhxAAAA18PnGZ0JEybohz/8oTo7OyVJN9xwg6xWn3MSACACcG0TzMbnoDNkyBD95je/UXJysqZPny7DMJSbmyu32x3I+gAAAK6Zz1Mylx/u19HRoSNHjujgwYM6ePCgiouL9dFHH13xnB0AAIBQG/S5p5iYGE2ZMkVTpkwJRD0AcN04/XJ9Up59g9vMYRo+n7oCAACINAQdAABgWgQdAKbCaSv/4DjCLAg6AADAtAg6AADAtAg6AEyD0y3+xfGEGRB0AACAaRF0AACAaRF0AJgCp1kCg+OKSEfQAQAApkXQAQAApkXQARDxOL0CoC8EHQAAYFoEHQAAYFoEHQARjdNWgccxRiQj6AAAANMi6AAAANMi6AAABsTpK0SqoAad2tpapaWlyTAMlZWVXfF+R0eH8vPzZRiGcnJy1NLS0vPemjVrZBiG0tLStGPHjp7XU1JSlJmZqaysLNnt9mB0AwAARIigBZ3u7m4tXbpUNTU1amxs1ObNm9XY2NirTWVlpRISEtTU1KQVK1aotLRUktTY2Cin06mjR4+qtrZWS5YsUXd3d896u3fv1gcffKCGhoZgdQdAGGCWIbg43ohEQQs69fX1MgxDqampGjZsmAoKCuRyuXq1cblcKiwslCTNnz9fdXV18nq9crlcKigoUExMjCZMmCDDMFRfXx+s0gEAQIQKWtDxeDxKTk7uWU5KSpLH4+mzjdVqVXx8vNrb2/td12Kx6J577tGUKVNUUVHR5/4rKipkt9tlt9t16tQpf3YNAACEqYi/GPmdd97RgQMHVFNToxdffFG//e1vr9qupKREDQ0Namho0NixY4NcJQB/4zRKaHDcEWmCFnRsNptaW1t7ltva2mSz2fps09XVpTNnzmjMmDH9rnv5n+PGjdPcuXM5pQUAAHoELehkZ2fL7XarublZnZ2dcjqdcjgcvdo4HA5VVVVJkqqrqzVz5kxZLBY5HA45nU51dHSoublZbrdbd911l86fP68vvvhCknT+/Hnt3LlTkyZNClaXAIQIswqhxfFHJLEGbUdWq8rLy5WXl6fu7m4tXrxYGRkZWrlypex2uxwOh4qKirRw4UIZhqHExEQ5nU5JUkZGhhYsWKD09HRZrVa9+OKLGjp0qP785z9r7ty5kr6eAXrwwQf1ve99L1hdAgAAYc7i9Xq9oS4i2Ox2O7eiAxGMGYXQaymbE+oSgB79jesRfzEygOhCyAEwGAQdAABgWgQdABGD2ZzwwWeBSEHQAQAApkXQAQBcE2Z1EAkIOgAAwLQIOgAiArMH4YnPBeGOoAMAAEyLoAMg7DFrEN74fBDOCDoAwhqDaGTgc0K4IugAAADTIugACFvMEkQWPi+EI4IOAAAwLYIOgLDE7AAAfyDoAAg7hJzIxWeHcEPQAQD4FWEH4YSgAwAATIugAyCsMBtgDnyOCBcEHQBhg8HRXPg8EQ4IOgDCAoOiOfG5ItQIOgBCjsEQQKAQdACEFCHH/PiMEUoEHQAhwwAYPVKefYPPGyFB0AEQEgx60YnPHcFmDXUBAKILAx0ufwdayuaEuBJEA2Z0AAQFpy7w1/g+IBiCGnRqa2uVlpYmwzBUVlZ2xfsdHR3Kz8+XYRjKyclRS0tLz3tr1qyRYRhKS0vTjh07fN4mgNAi4KA/l78ffEcQKEE7ddXd3a2lS5dq165dSkpKUnZ2thwOh9LT03vaVFZWKiEhQU1NTXI6nSotLdWrr76qxsZGOZ1OHT16VJ9++qlmzZqlY8eOSdKA2wQQXAxYuFbf/O5wWgv+ErSgU19fL8MwlJqaKkkqKCiQy+XqFUpcLpdWrVolSZo/f76WLVsmr9crl8ulgoICxcTEaMKECTIMQ/X19ZI04DYBBA6hBoHS13eLAITBClrQ8Xg8Sk5O7llOSkrSu+++22cbq9Wq+Ph4tbe3y+PxaOrUqb3W9Xg8kjTgNi+rqKhQRUWFJOnDDz+U3W73T8f87NSpUxo7dmyoyzAVjqn/XT6mN4W6EBPhe+obu/3HPrXjePpfOB/Tb17q8tei5q6rkpISlZSUhLqMAdntdjU0NIS6DFPhmPofx9T/OKb+xfH0v0g9pkG7GNlms6m1tbVnua2tTTabrc82XV1dOnPmjMaMGdPnur5sEwAARK+gBZ3s7Gy53W41Nzers7NTTqdTDoejVxuHw6GqqipJUnV1tWbOnCmLxSKHwyGn06mOjg41NzfL7Xbrrrvu8mmbAAAgegXt1JXValV5ebny8vLU3d2txYsXKyMjQytXrpTdbpfD4VBRUZEWLlwowzCUmJgop9MpScrIyNCCBQuUnp4uq9WqF198UUOHDpWkq24zkkXC6bVIwzH1P46p/3FM/Yvj6X+RekwtXq/XG+oiAAAAAoEnIwMAANMi6AAwhdzcXH344YeSpPb2dk2aNCnEFQEIBwQdAKbQ1NSkb3/725KkQ4cOKTMzM8QVAQgHBB0AEe/jjz+WzWbTkCFf/yft0KFDmjx5siRp8eLFoSwNQIgRdABEvIMHD/YEG0nav3+/Jk+erC+//FLx8fHavXu3SktLdfHixRBWCSAUCDoAIt4HH3zQE2LcbrdcLpcyMzN14MABvf/++/rjH/+otWvXKjY2NsSVAgi2qPkJCADmdfDgQcXGxuo73/mOJk+erPT0dFVVVSkuLk45OTkaMWJEqEsEECIEHQAR79ChQzpw4IBGjhzZ6/Xi4mK99NJL+tGPfqS9e/fq7rvvDlGFAEKFBwYCiGhffPGFpkyZomPHjoW6FABhiKADAABMi4uRAQCAaRF0AISllJQUDR8+XAsXLvSpfWVlpeLi4mSxWNTU1BTg6gBECoIOgLC1bds2vfLKK1e8/vbbb8tiseiHP/xhz2tFRUU6d+5cQOu52n4BhDeCDoCIcunSJT399NPKycmJiv0CuD4EHQAhcfz4cSUmJurAgQOSpE8//VRjx47Vnj17+l3vZz/7me655x7ddtttQagy9PsFcH0IOgBCYuLEiVq7dq0efvhhffnll3r00UdVWFio3NzcPtf5+OOP9fLLL2vlypXXtM97771Xo0ePvurfvffeG7D9AggdHhgIIGSKi4u1bds25eTkyGKxaOvWrf22X758uVavXq24uLhr2t/rr79+Tetd734BhA4zOgBCqri4WEeOHNFTTz2lmJiYPttt27ZNX3zxhfLz84NYXej2C8A/mNEBEDLnzp3TM888o6KiIq1atUrz5s1TYmLiVdvW1dWpoaFBf/u3fytJOnPmjIYOHarDhw/L5XL5tL/Zs2dr7969V33v7rvvVk1NTUD2CyB0eDIygJC5fEv4q6++qpKSEn3++ef6z//8T0lfP0fnV7/6lWbNmiXp6596OH/+fM+6Tz/9tG6++Wb96Ec/6hWOLBaL3G63DMPwS42+7hdAeGJGB0BIuFwu1dbW6vDhw5KkF154QVlZWdq4caMeeuihK9qPHDmy1492Dh8+XCNGjAh42AjVfgH4BzM6AMJSWlqaTpw4oblz56qqqmrA9v/+7/+uFStW6OLFi2psbFRqamoQqgQQ7gg6AADAtLjrCgAAmBZBBwAAmFZUXox80003KSUlJdRlAAAAP2hpadHp06ev+l5UBp2UlBQ1NDSEugwAAOAHdru9z/c4dQUAAEyLoAMAAEyLoAMAAEyLoAMAAEyLoAMAAEyLoAMAAEyLoAMAAEyLoAMAAEyLoAMAAEyLoAMAAEwrqEGntrZWaWlpMgxDZWVlV7zf0dGh/Px8GYahnJwctbS0SJJ27dqlKVOmKDMzU1OmTNFbb73Vs87+/fuVmZkpwzC0fPlyeb3eYHUHAACEuaAFne7ubi1dulQ1NTVqbGzU5s2b1djY2KtNZWWlEhIS1NTUpBUrVqi0tFTS1z/CuW3bNh0+fFhVVVVauHBhzzpPPvmk1q9fL7fbLbfbrdra2mB1CQAAhLmgBZ36+noZhqHU1FQNGzZMBQUFcrlcvdq4XC4VFhZKkubPn6+6ujp5vV7dcccduvnmmyVJGRkZunDhgjo6OnTixAmdPXtWU6dOlcVi0aJFi7Rly5ZgdQkAAIS5oAUdj8ej5OTknuWkpCR5PJ4+21itVsXHx6u9vb1Xm9dee0133nmnYmJi5PF4lJSU1O82L6uoqJDdbpfdbtepU6f81S0AABDGrKEuYDCOHj2q0tJS7dy5c9DrlpSUqKSkRFL/P+cOAADMI2gzOjabTa2trT3LbW1tstlsfbbp6urSmTNnNGbMmJ72c+fO1YYNGzRx4sSe9m1tbf1uEwAARK+gzehkZ2fL7XarublZNptNTqdTmzZt6tXG4XCoqqpK3/3ud1VdXa2ZM2fKYrHo888/15w5c1RWVqZp06b1tB8/frxGjRqlffv2KScnRxs2bNBTTz0VrC4BiHApz74R6hIkSS1lc0JdAmBaQZvRsVqtKi8vV15enm6//XYtWLBAGRkZWrlypbZu3SpJKioqUnt7uwzD0AsvvNBzC3p5ebmampr0k5/8RFlZWcrKytLJkyclSevWrdNjjz0mwzA0ceJEzZ49O1hdAhChUp59I2xCjhQ+gQswI4s3Ch88Y7fb1dDQEOoyAIRAOIcKZnaAa9PfuM6TkQFEjXAOOVL41wdEIoIOAIQRwg7gXwQdAABgWgQdAAgzzOoA/kPQARAVCA9AdCLoAEAYIpgB/kHQAWB6hAYgehF0ACBMEdCA60fQAWBqhAUguhF0AACAaRF0AACAaRF0AJgWp60AEHQAIIwR1oDrQ9ABYEoEBAASQQcAAJgYQQcAwhyzU8C1I+gAAADTIugAAADTIugAAADTIugAQATgOh3g2hB0AJgOoQDAZQQdAABgWgQdAIgQzFQBg0fQAWAqhAEA30TQAQAApkXQAQAAphXUoFNbW6u0tDQZhqGysrIr3u/o6FB+fr4Mw1BOTo5aWlokSe3t7ZoxY4bi4uK0bNmyXuvk5uYqLS1NWVlZysrK0smTJ4PRFQAICU7NAYNjDdaOuru7tXTpUu3atUtJSUnKzs6Ww+FQenp6T5vKykolJCSoqalJTqdTpaWlevXVVxUbG6vVq1fryJEjOnLkyBXb3rhxo+x2e7C6AiBMEQIA/LWgzejU19fLMAylpqZq2LBhKigokMvl6tXG5XKpsLBQkjR//nzV1dXJ6/VqxIgRmj59umJjY4NVLgAAMIGgBR2Px6Pk5OSe5aSkJHk8nj7bWK1WxcfHq729fcBtP/roo8rKytLq1avl9Xqv2qaiokJ2u112u12nTp26jp4AAIBIEfEXI2/cuFGHDx/W3r17tXfvXr3yyitXbVdSUqKGhgY1NDRo7NixQa4SAACEQtCCjs1mU2tra89yW1ubbDZbn226urp05swZjRkzZsDtStLIkSP14IMPqr6+3s+VA4gEXJ8D4GqCFnSys7PldrvV3Nyszs5OOZ1OORyOXm0cDoeqqqokSdXV1Zo5c6YsFkuf2+zq6tLp06clSZcuXdLrr7+uSZMmBa4TAAAgogTtriur1ary8nLl5eWpu7tbixcvVkZGhlauXCm73S6Hw6GioiItXLhQhmEoMTFRTqezZ/2UlBSdPXtWnZ2d2rJli3bu3Klbb71VeXl5unTpkrq7uzVr1iwVFxcHq0sAEBIpz76hlrI5oS4DiAgWb19X75qY3W5XQ0NDqMsA4CfReNqKoAP8t/7G9Yi/GBkAAKAvBB0AAGBaBB0AiEDReLoOuBaDDjrnz59Xd3d3IGoBgEFjwAfQnwGDzldffaVNmzZpzpw5GjdunG677TaNHz9e6enp+od/+Ac1NTUFo04AAIBBGzDozJgxQ8ePH9eaNWv0pz/9Sa2trTp58qTeeecdTZ06VaWlpfr1r38djFoBAAAGZcDbyy9duqQbbrih34340iaccHs5YB7RfuqK28yB67y93JcAE0khBwAARI9BPxl506ZN2rp1q4YOHSqv16v77rtPDzzwQCBqAwAAuC6DDjpvv/12r59mWLp0KUEHAACEpUHfXt7R0aE33nhDhw4d0vbt23XhwoVA1AUA8EG0X6MEDGTQQWfdunX6y1/+ou3bt+uzzz5TeXl5IOoCgAExyAMYiM+nrp5++mn9/Oc/14033qiHH344kDUBAAD4hc8zOiNHjpTD4dD58+clSTt27NC0adMCVhgAAMD18nlG57nnntOmTZuUm5urYcOGKS4uTmVlZYGsDQDgg5Rn3+B5OkAffA46dXV1Wr9+vUaMGKETJ07o5ZdfVlpaWiBrA4A+cX0OAF/4fOrq+eef1+rVq7Vnzx5VV1crPz9fb731ViBrAwAAuC4+z+h8M9RkZmaqpqZG8+bN0+9///uAFAYAAHC9Bn17+WXjx49XXV2dP2sBAFwjTuUBVzfgjM4nn3zi04ZGjx6tUaNGXXdBADAQBnUAvhow6BQWFspisai/Hzm3WCx65JFHtGjRIr8WBwAAcD0GDDq7d+++4rWOjg7FxMQEpCAAAAB/uaZrdFatWtXz77/73e/8VQsAAIBfDfrXyyUpLy9PFRUVGjFihP7whz/whGQAQcP1OQAGY9BBp6ioSPHx8Xr//fc1depUrVmzJhB1AQAAXLdBB53KykpduHBBBw4c0HvvvafHH39cL730UiBqAwAAuC4+X6OTmZmphx56SGvXrtWePXt066236vz584MKObW1tUpLS5NhGFf9nayOjg7l5+fLMAzl5OSopaVFktTe3q4ZM2YoLi5Oy5Yt67XO/v37lZmZKcMwtHz58n7vDgMAM+O0HnAln4PO22+/reLiYg0fPlxOp1OTJk3S9u3bfd5Rd3e3li5dqpqaGjU2Nmrz5s1qbGzs1aayslIJCQlqamrSihUrVFpaKkmKjY3V6tWr9dOf/vSK7T755JNav3693G633G63amtrfa4JAACYm89BJzExUbm5uVq+fLmqqqr03nvvyTAMn3dUX18vwzCUmpqqYcOGqaCgQC6Xq1cbl8ulwsJCSdL8+fNVV1cnr9erESNGaPr06YqNje3V/sSJEzp79qymTp0qi8WiRYsWacuWLT7XBCCyMGMBYLB8DjrHjh3rtfytb31Lhw4d8nlHHo9HycnJPctJSUnyeDx9trFarYqPj1d7e3u/20xKSup3m5dVVFTIbrfLbrfr1KlTPtcNAAAil88XIz/++OM6fvy4bDabJk+erIsXL2rSpEn68ssvdeONNwayRr8oKSlRSUmJJMlut4e4GgAIjJRn31BL2ZxQlwGEDZ+DzuUnJH/yySc6ePCgPvjgAx08eFBZWVkaMmSIPvzww37Xt9lsam1t7Vlua2uTzWa7apukpCR1dXXpzJkzGjNmTL/bbGtr63ebAAAgeg0YdLxerywWS8/yLbfcoltuuUX33Xdfz2tnz54dcEfZ2dlyu91qbm6WzWaT0+nUpk2berVxOByqqqrSd7/7XVVXV2vmzJm99v3Xxo8fr1GjRmnfvn3KycnRhg0b9NRTTw1YCwAAiA4DXqMzY8YM/eIXv7jiV8w7Ozv11ltvqbCwUP/1X/814I6sVqvKy8uVl5en22+/XQsWLFBGRoZWrlyprVu3Svr6YYTt7e0yDEMvvPBCr1vQU1JS9P3vf1//8R//oaSkpJ47ttatW6fHHntMhmFo4sSJmj179qAOAAAAMC+Ld4AHz1y8eFEvv/yyNm7cqI8++kgJCQm6cOGCvvrqK91zzz1asmSJ7rjjjmDV6xd2u10NDQ2hLgPAIHHXle+4TgfRpL9xfcBTV7GxsVqyZImWLFmiS5cu6fTp0xo+fLhGjx7t7zoBoE+EHADXwufby2tqanT33XcrNzdXJSUl2rdvXyDrAgAAuG4+B50lS5boZz/7mfbt26eSkhL94Ac/0ObNmwNZGwAAwHXxOeiMGzdO06ZNU0JCgmbNmqUdO3bo+eefD2RtAIBrxKk+4Gs+B50JEybohz/8oTo7OyVJN9xwg6zWQf/4OQAMGoM2gGvlc9AZMmSIfvOb3yg5OVnTp0+XYRjKzc2V2+0OZH0AAADXzOcpmcsP9+vo6NCRI0d08OBBHTx4UMXFxfroo4+ueM4OAABAqA363FNMTIymTJmiKVOmBKIeAICf8LtXwCBOXQFAKHB9DoDrQdABAACmRdABAACmRdABAACmRdABELa4PgfA9SLoAAAA0yLoAICJMSuGaEfQARCWGKAB+ANBBwAAmBZBBwBMjtkxRDOCDoCww8AMwF8IOgAAwLQIOgAAwLQIOgAQBTgdiGhF0AEAAKZF0AEAAKZF0AGAKMHpK0Qjgg6AsMJgDMCfghp0amtrlZaWJsMwVFZWdsX7HR0dys/Pl2EYysnJUUtLS897a9askWEYSktL044dO3peT0lJUWZmprKysmS324PRDQAAECGCFnS6u7u1dOlS1dTUqLGxUZs3b1ZjY2OvNpWVlUpISFBTU5NWrFih0tJSSVJjY6OcTqeOHj2q2tpaLVmyRN3d3T3r7d69Wx988IEaGhqC1R0AiEjMmCHaBC3o1NfXyzAMpaamatiwYSooKJDL5erVxuVyqbCwUJI0f/581dXVyev1yuVyqaCgQDExMZowYYIMw1B9fX2wSgcQJAzCAPwtaEHH4/EoOTm5ZzkpKUkej6fPNlarVfHx8Wpvb+93XYvFonvuuUdTpkxRRUVFn/uvqKiQ3W6X3W7XqVOn/Nk1AAAQpiL+YuR33nlHBw4cUE1NjV588UX99re/vWq7kpISNTQ0qKGhQWPHjg1ylQAQPpg5QzQJWtCx2WxqbW3tWW5ra5PNZuuzTVdXl86cOaMxY8b0u+7lf44bN05z587llBYQoRh8AQRC0IJOdna23G63mpub1dnZKafTKYfD0auNw+FQVVWVJKm6ulozZ86UxWKRw+GQ0+lUR0eHmpub5Xa7ddddd+n8+fP64osvJEnnz5/Xzp07NWnSpGB1CQAAhDlr0HZktaq8vFx5eXnq7u7W4sWLlZGRoZUrV8put8vhcKioqEgLFy6UYRhKTEyU0+mUJGVkZGjBggVKT0+X1WrViy++qKFDh+rPf/6z5s6dK+nrGaAHH3xQ3/ve94LVJQB+wmwOgECxeL1eb6iLCDa73c6t6EAYIegEX0vZnFCXAPhNf+N6xF+MDAAYPMIlogVBB0BIMeACCCSCDgAAMC2CDgBEKWbTEA0IOgAAwLQIOgAQxZjVgdkRdACEDIMsgEAj6ABAlCNwwswIOgBCgsEVQDAQdAAABE+YFkEHQNAxqAIIFoIOAAAwLYIOgKBiNgdAMBF0AAQNISe88fnAjAg6AIAehB2YDUEHQFAwgAIIBYIOAKAXQinMhKADALgCYQdmQdABEHAMmpGJzw1mQNABEFAMlgBCiaADIGAIOZGPzxCRjqADICAYIAGEA4IOAL8j5JhLyrNv8JkiYhF0APgVA6J58dkiEhF0APgNA6H58Rkj0hB0AFw3Tm1EFz5rRJKgBp3a2lqlpaXJMAyVlZVd8X5HR4fy8/NlGIZycnLU0tLS896aNWtkGIbS0tK0Y8cOn7cJIHAIONGLzx6RwhqsHXV3d2vp0qXatWuXkpKSlJ2dLYfDofT09J42lZWVSkhIUFNTk5xOp0pLS/Xqq6+qsbFRTqdTR48e1aeffqpZs2bp2LFjkjTgNgH4HwMcLvvmd6GlbE4IKwGuLmhBp76+XoZhKDU1VZJUUFAgl8vVK5S4XC6tWrVKkjR//nwtW7ZMXq9XLpdLBQUFiomJ0YQJE2QYhurr6yVpwG0CuH4EG/jir78nBB+Eg6AFHY/Ho+Tk5J7lpKQkvfvuu322sVqtio+PV3t7uzwej6ZOndprXY/HI0kDbvOyiooKVVRUSJI+/PBD2e12/3TMz06dOqWxY8eGugxT4Zhev5v+aplj6n9mPKZ2+49Dtm8zHs9QC+dj+s1LXf5a0IJOqJWUlKikpCTUZQzIbreroaEh1GWYCsfU/zim/scx9S+Op/9F6jEN2sXINptNra2tPcttbW2y2Wx9tunq6tKZM2c0ZsyYPtf1ZZsAACB6BS3oZGdny+12q7m5WZ2dnXI6nXI4HL3aOBwOVVVVSZKqq6s1c+ZMWSwWORwOOZ1OdXR0qLm5WW63W3fddZdP2wQAANEraKeurFarysvLlZeXp+7ubi1evFgZGRlauXKl7Ha7HA6HioqKtHDhQhmGocTERDmdTklSRkaGFixYoPT0dFmtVr344osaOnSoJF11m5EsEk6vRRqOqf9xTP2PY+pfHE//i9RjavF6vd5QFwEAABAIPBkZAACYFkEHgCnk5ubqww8/lCS1t7dr0qRJIa4IQDgg6AAwhaamJn3729+WJB06dEiZmZkhrghAOCDoAIh4H3/8sWw2m4YM+fo/aYcOHdLkyZMlSYsXLw5laQBCjKADIOIdPHiwJ9hI0v79+zV58mR9+eWXio+P1+7du1VaWqqLFy+GsEoAoUDQARDxPvjgg54Q43a75XK5lJmZqQMHDuj999/XH//4R61du1axsbEhrhRAsEXNT0AAMK+DBw8qNjZW3/nOdzR58mSlp6erqqpKcXFxysnJ0YgRI0JdIoAQIegAiHiHDh3SgQMHNHLkyF6vFxcX66WXXtKPfvQj7d27V3fffXeIKgQQKjwwEEBE++KLLzRlyhQdO3Ys1KUACEMEHQAAYFpcjAwAAEyLoAMgLKWkpGj48OFauHChT+0rKysVFxcni8WipqamAFcHIFIQdACErW3btumVV17pWb4cfuLi4hQXF6d77rmn572ioiKdO3cuIHX0t18A4Y27rgBElG3btmnWrFlRs18A14cZHQAhcfz4cSUmJurAgQOSpE8//VRjx47Vnj17QlsYAFMh6AAIiYkTJ2rt2rV6+OGH9eWXX+rRRx9VYWGhcnNz+13voYce0tixY3XPPffo4MGDg9rnvffeq9GjR1/179577w3YfgGEDqeuAIRMcXGxtm3bppycHFksFm3durXf9hs3btSdd94pr9erf/3Xf1VeXp4+/PBDjR492qf9vf7669dU5/XuF0DoMKMDIKSKi4t15MgRPfXUU4qJiem37bRp0zR8+HDdeOON+sd//EeNHj1ae/fuDXiNodovgOtH0AEQMufOndMzzzyjoqIirVq1Sp999tmg1rdYLBrMM09nz57dc+fUX//Nnj07YPsFEDqcugIQMk8//bTsdrt+9atfqaSkRE888YT+8z//86ptP/nkE7W2tio7O1tfffWVfvGLX+j06dOaNm2az/urqakZdI3+2C+A0CHoAAgJl8ul2tpaHT58WJL0wgsvKCsrSxs3btRDDz10RfsvvvhCTz75pI4fP67Y2FhlZWWppqZGY8aMCWidodovAP/gt64AhKW0tDSdOHFCc+fOVVVV1YDt//3f/10rVqzQxYsX1djYqNTU1CBUCSDcEXQAAIBpcTEyAAAwLYIOAAAwrai8GPmmm25SSkpKqMsAAAB+0NLSotOnT1/1vagMOikpKWpoaAh1GQAAwA/sdnuf73HqCgAAmBZBBwAAmBZBBwAAmBZBBwAAmBZBBwAAmBZBBwAAmFZU3l4OIDqlPPtGqEvo0VI2J9QlAFGBGR0AUSGcQo4UfvUAZkXQAWB64RoqwrUuwEwIOgBMLdzDRLjXB0S6oAad2tpapaWlyTAMlZWVXfF+R0eH8vPzZRiGcnJy1NLSIknatWuXpkyZoszMTE2ZMkVvvfVWzzr79+9XZmamDMPQ8uXL5fV6g9UdAAAQ5oIWdLq7u7V06VLV1NSosbFRmzdvVmNjY682lZWVSkhIUFNTk1asWKHS0lJJX/8I57Zt23T48GFVVVVp4cKFPes8+eSTWr9+vdxut9xut2pra4PVJQBhLlJmSyKlTiASBS3o1NfXyzAMpaamatiwYSooKJDL5erVxuVyqbCwUJI0f/581dXVyev16o477tDNN98sScrIyNCFCxfU0dGhEydO6OzZs5o6daosFosWLVqkLVu2BKtLAAAgzAUt6Hg8HiUnJ/csJyUlyePx9NnGarUqPj5e7e3tvdq89tpruvPOOxUTEyOPx6OkpKR+twkgOjFLAkCKsOfoHD16VKWlpdq5c+eg162oqFBFRYUk6dSpU/4uDQAAhKGgzejYbDa1trb2LLe1tclms/XZpqurS2fOnNGYMWN62s+dO1cbNmzQxIkTe9q3tbX1u83LSkpK1NDQoIaGBo0dO9avfQMAAOEpaEEnOztbbrdbzc3N6uzslNPplMPh6NXG4XCoqqpKklRdXa2ZM2fKYrHo888/15w5c1RWVqZp06b1tB8/frxGjRqlffv2yev1asOGDbr//vuD1SUA8BtOtQGBEbSgY7VaVV5erry8PN1+++1asGCBMjIytHLlSm3dulWSVFRUpPb2dhmGoRdeeKHnFvTy8nI1NTXpJz/5ibKyspSVlaWTJ09KktatW6fHHntMhmFo4sSJmj17drC6BAAAwpzFG4UPnrHb7WpoaAh1GQACKBJnSPj9K+Da9Deu82RkAKYTiSFHity6gXBG0AEAAKZF0AEAAKZF0AFgKpF++ifS6wfCDUEHAACYFkEHAACYFkEHgGmY5bSPWfoBhAOCDgAAMC2CDgAAMC2CDgAAMC2CDgBTMNt1LWbrDxAqBB0AAGBaBB0AAGBaBB0AAGBaBB0AEc+s17OYtV9AMBF0AACAaRF0AACAaRF0AACAaRF0AEQ0rmMB0B+CDgAAMC2CDgAAMC2CDgAAMC2CDgCEMa5BAq4PQQcAAJgWQQcAAJgWQQcAAJhWUINObW2t0tLSZBiGysrKrni/o6ND+fn/f3v3HhxVebhx/FlZQjRAIEg03Q23OTaQCw2yMXREJyAQrbqMJgNBxGAicQRGpWMn/tFSO2iJHZ1qG5wxTuykKqwdZso6IoEawQKWppFbNUU2mpQkUi6Rgg2QkHV/f3RMzY9LNmR3T/bk+/kHdvOek+e8y2QfznmzZ6EMw1B2draampokSW1tbZo1a5aGDx+ulStX9tgmJydHKSkpyszMVGZmpo4fPx6JQwEAAFEgYkXH7/drxYoV2rJli+rr67VhwwbV19f3GFNZWanRo0eroaFBq1atUmlpqSQpNjZWa9as0QsvvHDJfb/11lvav3+/9u/fr8TExLAfC4CBYbAs1B0sxwmEQ8SKTm1trQzD0KRJkxQTE6OCggJ5vd4eY7xerwoLCyVJ+fn5qqmpUSAQUFxcnGbOnKnY2NhIxQUAABYQsaLT2tqq5OTk7sdOp1Otra2XHWO32xUfH6+2trZe9/3www8rMzNTa9asUSAQCG1wAAAQtaJ+MfJbb72lv//979q5c6d27typN95445LjKioq5HK55HK5dOLEiQinBAAAZohY0XE4HGpubu5+3NLSIofDcdkxXV1dOn36tMaMGdPrfiVpxIgReuCBB1RbW3vJcSUlJaqrq1NdXZ3Gjh3bn0MBMAAMtnUrg+14gVCJWNHJysqSz+dTY2OjOjs75fF45Ha7e4xxu92qqqqSJG3cuFGzZ8+WzWa77D67urp08uRJSdKFCxf07rvvKj09PXwHAQAAooo9Yt/Ibld5eblyc3Pl9/tVVFSktLQ0rV69Wi6XS263W8XFxVqyZIkMw1BCQoI8Hk/39hMmTNCZM2fU2dmpTZs2adu2bRo/frxyc3N14cIF+f1+zZkzR8uWLYvUIQEAgAHOFhiEq3ddLpfq6urMjgGgHwbjpZymsrvNjgAMSFd6X4/6xcgAAACXQ9EBEHUG49kcafAeN9AfFB0AAGBZFB0AAGBZfS467e3t8vv94cgCAAAQUr0WnW+++Ubr16/X3XffrcTERE2ePFlJSUlKTU3VT37yEzU0NEQiJwAAQJ/1WnRmzZqlzz//XGvXrtW//vUvNTc36/jx49q1a5dmzJih0tJSvfnmm5HICgCDfkHuYD9+oK96/cDA999/X0OHDr3o+YSEBOXl5SkvL08XLlwISzgAAID+6PWMzqVKztWMAQAAiLQ+3wJi/fr1eueddzRkyBAFAgHde++9WrRoUTiyAQAA9Eufi86HH37Y4x5UK1asoOgAAIABqc9Fp6OjQ5s3b1ZycrJaWlp07ty5cOQCgIuwEBdAX/X5c3ReeeUVnTp1Su+9956++uorlZeXhyMXAABAvwV9RueJJ57QSy+9pOuuu04PPvhgODMBAACERNBndEaMGCG326329nZJ0tatW3XrrbeGLRgAAEB/BX1G59lnn9X69euVk5OjmJgYDR8+XGVlZeHMBgC4hAlPb1ZT2d1mxwCiQtBFp6amRq+99pri4uJ09OhRvf7660pJSQlnNgAAgH4J+tLVc889pzVr1mjHjh3auHGjFi5cqA8++CCc2QAAAPol6DM63y01GRkZ2rJli/Ly8vTRRx+FJRgAAEB/9fnXy7+VlJSkmpqaUGYBAAAIqV7P6Bw5ciSoHY0aNUojR47sdyAAuBQ+LLAnFiQDwem16BQWFspmsykQCFx2jM1m09KlS/XQQw+FNBwAAEB/9Fp0tm/fftFzHR0dGjZsWFgCAQAAhMpVrdF55plnuv++e/fuUGUBAAAIqT7f1FOScnNzVVFRobi4OP3jH//gE5IBAMCA1OeiU1xcrPj4eO3bt08zZszQ2rVrw5ELALqxEPnSWJAM9K7Pl64qKyv13HPP6dlnn9UNN9ygRx99NOhtq6urlZKSIsMwLnn7iI6ODi1cuFCGYSg7O1tNTU2SpLa2Ns2aNUvDhw/XypUre2zz8ccfKyMjQ4Zh6PHHH7/iomkAADC4BF10MjIytHjxYj3//PPasWOHxo8fr/b2dr366qtBbe/3+7VixQpt2bJF9fX12rBhg+rr63uMqays1OjRo9XQ0KBVq1aptLRUkhQbG6s1a9bohRdeuGi/jz32mF577TX5fD75fD5VV1cHe0gAAMDigi46H374oZYtW6Zrr71WHo9H6enpeu+994L+RrW1tTIMQ5MmTVJMTIwKCgrk9Xp7jPF6vSosLJQk5efnq6amRoFAQHFxcZo5c6ZiY2N7jD969KjOnDmjGTNmyGaz6aGHHtKmTZuCzgQAAKwt6KKTkJCgnJwcPf7446qqqtLf/vY3GYYR9DdqbW1VcnJy92On06nW1tbLjrHb7YqPj1dbW9sV9+l0Oq+4TwDRjfU5V8b8AFcWdNE5fPhwj8c33XSTDh48GPJA4VJRUSGXyyWXy6UTJ06YHQcAAERA0L919eijj+rzzz+Xw+HQ1KlTdf78eaWnp+vs2bO67rrret3e4XCoubm5+3FLS4scDsclxzidTnV1den06dMaM2bMFffZ0tJyxX1+q6SkRCUlJZIkl8vVa14AABD9gj6js337dh05ckRvv/227rnnHhmGoXPnzikzM1OTJ0/udfusrCz5fD41Njaqs7NTHo9Hbre7xxi3262qqipJ0saNGzV79mzZbLbL7jMpKUkjR47Unj17FAgE9Pvf/17z588P9pAAAIDF9XpGJxAI9Cgb48aN07hx43Tvvfd2P3fmzJnev5HdrvLycuXm5srv96uoqEhpaWlavXq1XC6X3G63iouLtWTJEhmGoYSEBHk8nu7tJ0yYoDNnzqizs1ObNm3Stm3blJqaqldeeUVLly7VuXPndNddd+muu+7q6xwAAACLsgV6+eCZnJwc5eXlaf78+Ro3blz3852dndq1a5eqqqo0a9YsLV26NNxZQ8blcqmurs7sGAB6wULb4PHBgRjMrvS+3usZnerqar3++utatGiRvvjiC40ePVrnzp3TN998o3nz5unJJ5/UtGnTQh4aAACgv3otOrGxsVq+fLmWL1+uCxcu6OTJk7r22ms1atSoCMQDAAC4ekEvRt6yZYtuu+025eTkqKSkRHv27AlnLgAAgH4LuugsX75cL774ovbs2aOSkhI99dRT2rBhQzizARjEWJ/TN8wXcGlBf45OYmKibr31VknSnDlz9MMf/lDZ2dlatGhR2MIBAAD0R9BndCZOnKif/vSn6uzslCQNHTpUdnvQPQkAACDigi4611xzjf74xz8qOTlZM2fOlGEYysnJkc/nC2c+AACAqxb0KZn169dLkjo6OvTJJ5/owIEDOnDggJYtW6YvvvhCR44cCVtIAIML600AhEqfrz0NGzZM06dP1/Tp08ORBwAAIGSCvnQFAAAQbSg6AADAsig6AGARrG0CLkbRAQAAlkXRAQAAlkXRAQAAlkXRATCgsM6kf5g/oCeKDgAAsCyKDgAAsCyKDgAAsCyKDoABg/UlocE8Av9D0QEAAJZF0QEAAJZF0QEAAJZF0QEwILCuJLSYT+C/KDoAAMCyKDoAAMCyIlp0qqurlZKSIsMwVFZWdtHXOzo6tHDhQhmGoezsbDU1NXV/be3atTIMQykpKdq6dWv38xMmTFBGRoYyMzPlcrkicRgAQozLLOHBvAKSPVLfyO/3a8WKFfrTn/4kp9OprKwsud1upaamdo+prKzU6NGj1dDQII/Ho9LSUr399tuqr6+Xx+PRp59+qi+//FJz5szR4cOHNWTIEEnS9u3bdf3110fqUAAAQJSI2Bmd2tpaGYahSZMmKSYmRgUFBfJ6vT3GeL1eFRYWSpLy8/NVU1OjQCAgr9ergoICDRs2TBMnTpRhGKqtrY1UdAAAEKUiVnRaW1uVnJzc/djpdKq1tfWyY+x2u+Lj49XW1nbFbW02m+bNm6fp06eroqIiAkcCAACiRcQuXYXLrl275HA4dPz4cc2dO1eTJ0/W7bffftG4ioqK7iJ04sSJSMcEcBmsIwEQThE7o+NwONTc3Nz9uKWlRQ6H47Jjurq6dPr0aY0ZM+aK2377Z2Jiou67777LXtIqKSlRXV2d6urqNHbs2JAeGwAAGJgiVnSysrLk8/nU2Niozs5OeTweud3uHmPcbreqqqokSRs3btTs2bNls9nkdrvl8XjU0dGhxsZG+Xw+3XLLLWpvb9fXX38tSWpvb9e2bduUnp4eqUMCAAADXMQuXdntdpWXlys3N1d+v19FRUVKS0vT6tWr5XK55Ha7VVxcrCVLlsgwDCUkJMjj8UiS0tLStGDBAqWmpsput2vdunUaMmSIjh07pvvuu0/Sf88APfDAA7rzzjsjdUgAMOBNeHqzmsruNjsGYBpbIBAImB0i0lwul+rq6syOAUCs0YkEig6s7krv63wyMgDTUHIig3nGYEbRAQAAlkXRAQAAlkXRAWAKLqdEFvONwYqiAwAALIuiAwAALIuiAyDiuIxiDuYdgxFFBwAAWBZFB0BEcVbBXMw/BhuKDgAAsCyKDoCI4WzCwMDrgMGEogMAACyLogMgIjiLAMAMFB0AAGBZFB0AGIQ4w4bBgqIDAAAsi6IDIOw4ezAw8bpgMKDoAAgr3kwHNl4fWB1FBwAAWBZFB0DYcLYgOvA6wcooOgDCgjfP6MLrBaui6AAIOd40AQwUFB0AIUXJiV68drAiig4AoBtlB1ZD0QEQMrxJWgOvI6zEbnYAANGPN0br+fY1bSq72+QkQP9E9IxOdXW1UlJSZBiGysrKLvp6R0eHFi5cKMMwlJ2draampu6vrV27VoZhKCUlRVu3bg16nwDCZ8LTmyk5Fsfri2gXsaLj9/u1YsUKbdmyRfX19dqwYYPq6+t7jKmsrNTo0aPV0NCgVatWqbS0VJJUX18vj8ejTz/9VNXV1Vq+fLn8fn9Q+wQQHrwBDh4UWkSziF26qq2tlWEYmjRpkiSpoKBAXq9Xqamp3WO8Xq+eeeYZSVJ+fr5WrlypQCAgr9ergoICDRs2TBMnTpRhGKqtrZWkXvcJIDR4o8N3/w1wSQvRImJFp7W1VcnJyd2PnU6n/vrXv152jN1uV3x8vNra2tTa2qoZM2b02La1tVWSet0ngL6h0CAYV/p3QgnCQDJoFiNXVFSooqJCknTo0CG5XC6TE13aiRMnNHbsWLNjWApz2jfXBzGGOQ09K82py/VzsyNYaj4HioE8p99d0/v/RazoOBwONTc3dz9uaWmRw+G45Bin06muri6dPn1aY8aMueK2ve3zWyUlJSopKQnlIYWFy+VSXV2d2TEshTkNPeY09JjT0GI+Qy9a5zRii5GzsrLk8/nU2Niozs5OeTweud3uHmPcbreqqqokSRs3btTs2bNls9nkdrvl8XjU0dGhxsZG+Xw+3XLLLUHtEwAADF4RO6Njt9tVXl6u3Nxc+f1+FRUVKS0tTatXr5bL5ZLb7VZxcbGWLFkiwzCUkJAgj8cjSUpLS9OCBQuUmpoqu92udevWaciQIZJ0yX0CAABIki0QCATMDoH/qaioiIpLbNGEOQ095jT0mNPQYj5DL1rnlKIDAAAsi3tdAQAAy6LoDCDcziK0mpubNWvWLKWmpiotLU0vv/yy2ZEswe/3a9q0abrnnnvMjmIJ//73v5Wfn6/JkydrypQp+stf/mJ2pKj361//WmlpaUpPT9eiRYt0/vx5syNFnaKiIiUmJio9Pb37ua+++kpz587VTTfdpLlz5+rUqVMmJgweRWeA4HYWoWe32/Xiiy+qvr5ee/bs0bp165jTEHj55Zc1ZcoUs2NYxhNPPKE777xThw4d0oEDB5jbfmptbdVvfvMb1dXV6ZNPPpHf7+/+xRYEb+nSpaquru7xXFlZme644w75fD7dcccdUfMfcorOAPHdW2TExMR0384CVy8pKUk333yzJGnEiBGaMmVK9ydq4+q0tLRo8+bNeuSRR8yOYgmnT5/Wn//8ZxUXF0uSYmJiNGrUKHNDWUBXV5fOnTunrq4unT17Vt/73vfMjhR1br/9diUkJPR4zuv1qrCwUJJUWFioTZs2mZCs7yg6A8SlbpHBm3LoNDU1ad++fcrOzjY7SlR78skn9atf/UrXXMOPjlBobGzU2LFj9fDDD2vatGl65JFH1N7ebnasqOZwOPTUU09p3LhxSkpKUnx8vObNm2d2LEs4duyYkpKSJEk33nijjh07ZnKi4PDTCpb3n//8R3l5eXrppZc0cuRIs+NErXfffVeJiYmaPn262VEso6urS3v37tVjjz2mffv2KS4uLmouBwxUp06dktfrVWNjo7788ku1t7frzTffNDuW5dhsNtlsNrNjBIWiM0AEc4sM9N2FCxeUl5enxYsX6/777zc7TlTbvXu33nnnHU2YMEEFBQX64IMP9OCDD5odK6o5nU45nc7uM435+fnau3evyami2/vvv6+JEydq7NixGjp0qO6//3599NFHZseyhBtuuEFHjx6VJB09elSJiYkmJwoORWeA4HYWoRcIBFRcXKwpU6boxz/+sdlxot7atWvV0tKipqYmeTwezZ49m/8p99ONN96o5ORkffbZZ5Kkmpoapaammpwquo0bN0579uzR2bNnFQgEVFNTwwLvEPnubZqqqqo0f/58kxMFZ9DcvXygu9wtMnD1du/erTfeeEMZGRnKzMyUJP3yl7/Uj370I3ODAd/x29/+VosXL1ZnZ6cmTZqk3/3ud2ZHimrZ2dnKz8/XzTffLLvdrmnTpkXlp/mabdGiRdqxY4dOnjwpp9OpX/ziF3r66ae1YMECVVZWavz48frDH/5gdsyg8MnIAADAsrh0BQAALIuiAwAALIuiAwAALIuiAwAALIuiAwAALIuiAwAALIuiAwAALIuiA8AScnJydOjQIUlSW1ub0tPTTU4EYCCg6ACwhIaGBn3/+9+XJB08eFAZGRkmJwIwEFB0AES9f/7zn3I4HLrmmv/+SDt48KCmTp0qSSoqKjIzGgCTUXQARL0DBw50FxtJ+vjjjzV16lSdPXtW8fHx2r59u0pLS3X+/HkTUwIwA0UHQNTbv39/d4nx+Xzyer3KyMjQ3r17tW/fPn322Wd6/vnnFRsba3JSAJHG3csBRL0DBw4oNjZWP/jBDzR16lSlpqaqqqpKw4cPV3Z2tuLi4syOCMAkFB0AUe/gwYPau3evRowY0eP5ZcuW6dVXX9XPfvYz7dy5U7fddptJCQGYxRYIBAJmhwCAq/X1119r+vTpOnz4sNlRAAxAFB0AAGBZLEYGAACWRdEBAACWRdEBAACWRdEBAACWRdEBAACWRdEBAACWRdEBAACWRdEBAACW9X9FqFd+8LBFkAAAAABJRU5ErkJggg==\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "b.\tIn these experiments, subjects need to report whether the box is\n",
    "heavier or lighter than the reference box of 6 N. We can model this\n",
    "forced-choice reporting by picking an estimated weight at random from\n",
    "the posterior distribution and saying the box we are lifting is\n",
    "lighter than the reference if our number is less than 6 and heavier\n",
    "if it is greater than 6.\n",
    "What is the probability of claiming that the small box (prior\n",
    "for $\\mu=3.5$) is actually heavier before lifting any box and\n",
    "then after each lift of the box in part (a).\n",
    "\n",
    "*Answer*\n",
    "\n",
    "We can repeat the loop in the last exercise, but this time each time\n",
    "through the loop we will check the probability of picking a random number\n",
    "from the posterior that is greater than 6.\n",
    "\n",
    "This is easy to do, all we need to do is check how much of the posterior\n",
    "is greater than 6."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The probability based only on the prior is 0.01\n",
      "The probability after the 1-th lift of 6 is 0.04\n",
      "The probability after the 2-th lift of 6 is 0.07\n",
      "The probability after the 3-th lift of 7 is 0.22\n",
      "The probability after the 4-th lift of 7 is 0.40\n",
      "The probability after the 5-th lift of 4 is 0.15\n",
      "The probability after the 6-th lift of 5 is 0.09\n"
     ]
    }
   ],
   "source": [
    "p_6 = sum(p_mu[mu_k>6])\n",
    "print(f'The probability based only on the prior is {p_6:.2f}')\n",
    "\n",
    "this_p_mu = p_mu # Set initial prior\n",
    "for i,x_i in enumerate(x):\n",
    "    p_x_mu = norm.pdf(x_i, loc=mu_k, scale=sigma_x)\n",
    "    p_mu_x = p_x_mu*this_p_mu\n",
    "    p_mu_x = p_mu_x / sum(p_mu_x)\n",
    "\n",
    "    p_6 = sum( p_mu_x[mu_k > 6])\n",
    "    print(f'The probability after the {i+1}-th lift of {x_i} is {p_6:.2f}')\n",
    "\n",
    "    this_p_mu = p_mu_x"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "c. Determine whether the results in parts (a) and (b) are consistent with the graphs on\n",
    "slide 14 (copied below). How are the consistent or inconsistent?\n",
    "\n",
    "*Answer*\n",
    "\n",
    "They are consistent because the actual esimate of how heavy the box is\n",
    "and the probability of saying it is heavior goes up on subsequent lifts.\n",
    "\n",
    "## Problem 3\n",
    "\n",
    "You are running an experiment where a participant is facing several\n",
    " speakers and several lights that are placed in a semicircle (±90◦) around\n",
    " them. As the experimenter, you can manipulate the average location and noise\n",
    " of auditory and visual feedback by controlling which lights and speakers\n",
    " will turn on. On a particular trial, you turn on a few lights to light up\n",
    " so that it is hard to pinpoint exactly where\n",
    " (${{\\mu }_{v}}=-{{10.0}^{\\circ }},{{\\sigma }_{v}}={{5}^{\\circ }}$).\n",
    " Simultaneously, you set a few speakers to produce a tone\n",
    " (${{\\mu }_{a}}={{25.0}^{\\circ }},{{\\sigma }_{a}}={{15}^{\\circ }}$). The\n",
    " participant has to point towards a single location that best captures the\n",
    " source of the light and sound. Assuming they have no prior inclination to\n",
    " point in any specific direction, where do you expect them to point? If\n",
    " a particular subject points towards towards  $-{{2}^{\\circ }}$,\n",
    " is their behavior still consistent with Bayesian estimation? Where\n",
    " would they need to point to convince you that their estimation\n",
    " process is not Bayesian?\n",
    "\n",
    "*Answer*\n",
    "\n",
    "We will treat this as a two step process of including information into\n",
    "prior knowledge. The prior knowledge is flat:"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "outputs": [],
   "source": [
    "mu_0 = 0\n",
    "sigma_0 = 10000"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now let's add in the visual estimate:"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Including the visual input we have an estimate of -10.0 with an uncertainty of 5.00\n"
     ]
    }
   ],
   "source": [
    "mu_v = -10\n",
    "sigma_v = 5\n",
    "\n",
    "# This is just to make the equations easier to type\n",
    "tau_v = 1/sigma_v**2\n",
    "tau_0 = 1/sigma_0**2\n",
    "\n",
    "mu_v1 = (tau_v*mu_v + tau_0*mu_0) / (tau_v+tau_0)\n",
    "sigma_v1 = np.sqrt(1 / (tau_v + tau_0))\n",
    "\n",
    "print(f'Including the visual input we have an estimate of {mu_v1:.1f} with an uncertainty of {sigma_v1:.2f}')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can see that if you have no prior information, then the estimate with the new data is just what\n",
    "is indicated by that data. We could have actually skipped this step and just started with incoporating\n",
    "the visual into the auditory information.\n",
    "\n",
    "Now we add in the auditory information."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Including both inputs, we have an estimate of -6.5 with an uncertainty of 4.74\n"
     ]
    }
   ],
   "source": [
    "mu_a = 25\n",
    "sigma_a = 15\n",
    "\n",
    "tau_v1 = 1/sigma_v1**2\n",
    "tau_a = 1/sigma_a**2\n",
    "\n",
    "mu_av = (tau_v1*mu_v1 + tau_a*mu_a)/(tau_v1+tau_a)\n",
    "sigma_av = np.sqrt(1/(tau_v1+tau_a))\n",
    "\n",
    "print(f'Including both inputs, we have an estimate of {mu_av:.1f} with an uncertainty of {sigma_av:.2f}')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "A subject who points at $-2^{\\circ}$ is within the uncertainty and is thus reasonably consistent with\n",
    "the Bayesian estimate.\n",
    "\n",
    "If they pointed at a point more than twice the uncertainty, that would certainly be inconsistent with\n",
    "Bayesian estimate. That is, anywhere to the left of $\\mu_{av}-2\\sigma_{av}$ or to the right of\n",
    "$\\mu_{av}+2\\sigma_{av}$"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Esimates to the right of -16.0 or the left of 3.0 are not Bayesian estimates\n",
      "The probability to point outside one $\\sigma$ of $\\mu$ is 0.317\n",
      "The probability to point outside two $\\sigma$ of $\\mu$ is 0.046\n",
      "The probability to point outside 3 $\\sigma$ of $\\mu$ is 0.003\n"
     ]
    }
   ],
   "source": [
    "print(f'Esimates to the right of {mu_av-2*sigma_av:.1f} or the left of {mu_av+2*sigma_av:.1f} are not Bayesian estimates')\n",
    "\n",
    "p_1_sigma = norm.cdf(mu_av+sigma_av, loc=mu_av, scale=sigma_av) - norm.cdf(mu_av-sigma_av, loc=mu_av, scale=sigma_av)\n",
    "p_outside_1_sigma = 1 - p_1_sigma\n",
    "\n",
    "print(f'The probability to point outside one $\\sigma$ of $\\mu$ is {p_outside_1_sigma:.3f}')\n",
    "\n",
    "p_2_sigma = norm.cdf(mu_av+2*sigma_av, loc=mu_av, scale=sigma_av) - norm.cdf(mu_av-2*sigma_av, loc=mu_av, scale=sigma_av)\n",
    "p_outside_2_sigma = 1 - p_2_sigma\n",
    "\n",
    "print(f'The probability to point outside two $\\sigma$ of $\\mu$ is {p_outside_2_sigma:.3f}')\n",
    "\n",
    "p_3_sigma = norm.cdf(mu_av+3*sigma_av, loc=mu_av, scale=sigma_av) - norm.cdf(mu_av-3*sigma_av, loc=mu_av, scale=sigma_av)\n",
    "p_outside_3_sigma = 1 - p_3_sigma\n",
    "\n",
    "print(f'The probability to point outside 3 $\\sigma$ of $\\mu$ is {p_outside_3_sigma:.3f}')\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}