{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Chapter 5: Bayesian Model fitting I"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1 (20 points)\n",
    "\n",
    "We observe $r$ heads in $n$ tosses, and we adopt a beta prior with parameters ($\\alpha_{prior}$ , $\\beta_{prior}$).\n",
    "What is the expectation value (the mean) of $p$ and what is its most likely value (the mode), for\n",
    "\n",
    "(a) the general case of this beta prior, and\n",
    "\n",
    "(b) a uniform prior (a special case of the beta prior, as explained in the script.)\n",
    "\n",
    "What form does the posterior PDF have if you have no data, and what is its mean and mode, for\n",
    "\n",
    "(c) the general case of this beta prior, and\n",
    "\n",
    "(d) the special case of the uniform prior. \n",
    "\n",
    "(e) Discuss briefly whether the mean or the mode is to be preferred as an estimator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2 (20 points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A production chain of electronic components produces items that have a probability $p$ of being defective. The chain manager does not know $p$, but from past experience and quality assesments, the manager expects this probability to be equal to  $4\\%$. In addition, the manager found some uncertainty about $p$ and assumes a standard deviation of  $2\\%$. \n",
    "\n",
    "The manager comes to you, expert in statistics, and after discussion you convince the manager to use a Beta distribution to model the uncertainty about  $p$. \n",
    "\n",
    "1) Explain why you think a Beta prior is reasonable or not.\n",
    "\n",
    "2) What values do you need to set the two parameters of the distribution in order to match the manager's priors about the expected value and the standard deviation of $p$?\n",
    "\n",
    "After choosing the parameters of the Beta distribution so as to represent the priors about the probability of producing a defective item, the manager now wants to get a better estimate of the posterior pdf of $p$ by observing new data. He decides to inspect a production lot of $100$ items, and finds that $3$ of the items in the lot are defective. \n",
    "\n",
    "3) How should the manager change the parameters of the Beta distribution in order to take this new information into account?\n",
    "\n",
    "4) After updating the parameters of the Beta distribution, what are the new expected value and the new standard deviation of the probability of finding a defective item?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3 (20 points)\n",
    "\n",
    "You have a 3D posterior PDF, $p(\\omega, \\mu, v_r)$, over the parallax, $\\omega$, proper motion, $\\mu$, and\n",
    "radial velocity, $v_r$ of a star. Write down an expression of the 3D posterior of this in terms of the\n",
    "distance, $r$, tangential velocity $v_t$, and radial velocity, $v_r$, where $r = 1/\\omega$ and $v_t = a\\,\\mu/\\omega$ where a is $a$ known constant. Check your result by checking the dimensions (probabilty is dimensionless)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercice 4 (40 points)\n",
    "\n",
    "\n",
    "A lighthouse is somewhere off a piece of straight coastline at a position $x$ along the shore and a distance $y$ out at sea. It emits short flashes at random times and hence at random azimuths, $\\theta$, from its position, i.e., $p(\\theta) = constant$. You see these flashes on the coast while walking along the shore. You record your position $D_k$ along the coastline at the instant you see a flash, but you do not record the direction the flash came from.\n",
    "\n",
    "* You will find the data, $\\{D_k\\}$, in the file `lighthouse.dat`. \n",
    "* We suppose all distances to be in $km$.\n",
    "\n",
    "We will suppose to limit the exploration space that the lighthouse is to be somewhere in a rectangle $-2 < x < 2$ , $0 < y < 2$ , with uniform prior where $x=0$ coincides with $D=0$. We take the origin in distance at the coast, and simply center the box of interests along the coastline.\n",
    "\n",
    "Below is a schematic to help you understand the situation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzAAAAFqCAYAAADFgr+jAAAAiXpUWHRSYXcgcHJvZmlsZSB0eXBl\nIGV4aWYAAHjaVY7rDcQwCIP/Z4qOQICYME5VNVI3uPEPmj7uPkVgWZFx2T/HKEtSiYs263CAAnV1\nXkN0mghRZaq5Y06uLTUUv3YRngLejfT9qJd/0wQdw9QMDRs2jnTehYVjZlCmUtbwn2vnS/HvY30K\nzBZfQzcsSfEGfW8AAAoGaVRYdFhNTDpjb20uYWRvYmUueG1wAAAAAAA8P3hwYWNrZXQgYmVnaW49\nIu+7vyIgaWQ9Ilc1TTBNcENlaGlIenJlU3pOVGN6a2M5ZCI/Pgo8eDp4bXBtZXRhIHhtbG5zOng9\nImFkb2JlOm5zOm1ldGEvIiB4OnhtcHRrPSJYTVAgQ29yZSA0LjQuMC1FeGl2MiI+CiA8cmRmOlJE\nRiB4bWxuczpyZGY9Imh0dHA6Ly93d3cudzMub3JnLzE5OTkvMDIvMjItcmRmLXN5bnRheC1ucyMi\nPgogIDxyZGY6RGVzY3JpcHRpb24gcmRmOmFib3V0PSIiCiAgICB4bWxuczpleGlmPSJodHRwOi8v\nbnMuYWRvYmUuY29tL2V4aWYvMS4wLyIKICAgIHhtbG5zOnRpZmY9Imh0dHA6Ly9ucy5hZG9iZS5j\nb20vdGlmZi8xLjAvIgogICBleGlmOlBpeGVsWERpbWVuc2lvbj0iODE2IgogICBleGlmOlBpeGVs\nWURpbWVuc2lvbj0iMzYyIgogICB0aWZmOkltYWdlV2lkdGg9IjgxNiIKICAgdGlmZjpJbWFnZUhl\naWdodD0iMzYyIgogICB0aWZmOk9yaWVudGF0aW9uPSIxIi8+CiA8L3JkZjpSREY+CjwveDp4bXBt\nZXRhPgogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgCiAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAKICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAKICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIAogICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgCiAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAKICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAKICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIAogICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIAog\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgCiAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAKICAgICAgICAgICAgICAgICAgICAgICAgICAg\nCjw/eHBhY2tldCBlbmQ9InciPz4uWXaiAAAABHNCSVQICAgIfAhkiAAAIABJREFUeNrs3WdAFNff\nhuFn6SCCoqDYC/beS+y995bExMRoNJr4TzPFGKMmpve8Eo2JiVGj2HvvvWPHXlARRUF62YV9P5ig\nhkUhIku5r0/DnrM7M7+dXebZOTNjWOt/3SwAAAAAyAZsKAEAAAAAAgwAAAAAEGAAAAAAEGAAAAAA\ngAADAAAAAAQYAAAAAAQYAAAAACDAAAAAAAABBgAAAEDuYEcJkJW8PWyAju3bSiEAAMATU71uY301\nbSGFIMAAj+/Yvq3yXX2WQgAAgCfCzcVeA5uXohDZGEPIAAAAABBgAAAAAIAAAwAAACDX4hwY/Cfx\n8XFau2SOdm5apYtnTio6KlJ58rqplE9FNWrRXp16DZSjkzOFAgAAAAEG1nXmxGF98s7LuhF05YHH\nw8Nu68j+nTqyf6fm/T5ZH3wxVVVrN6BgAAAAyDAMIUO6XDhzQu8M65sivPxb6K0bendYXx3et4Oi\nAQAAIMNwBAZplpSUqC/HjlJsdJQqVauj7k8PVrXajZSvQEFFht+R/97t+mvad7py6ZwkyWQyatI7\nL2v60l3K656PAgIAAOCxcQQGabZ781pdPHtS/V58Vd//uUItO/ZSwULesrOzV/4CnmrVqZcmz12v\nBs3aJj8nIjxMi2f/QvEAAABAgEHm2rxmsbr2e0Evjfog1T6Ojk4a8/kUeRcrmfzYri1rKR4AAAAI\nMMhcgRfPatCIdx/Zz8nZRX0HjUz++/rVSxQPAAAAGYJzYJBmvyzYkua+dRo1T562d3CgeAAAAMgQ\nHIHBE1HQyzt5upB3MQoCAACADMERmEf49L3h2rp2aartPZ8ZquGjJz7w2MwpX2vW1G8eeKx5++4a\n8/mU/7QM7Wt5Z/h6rfW//kTrFhUZnjxds35TNiQAAABkCI7APMKoMV9o5LuTUgyDqtu4pf5a558i\nvEjSwGFvafaaQypRprzsHRw0+uMf9daE73NV3YKuXEyebtOlLxsSAAAAMgRHYB7B1c1d3QYMlr2j\no76f+Hby483bdVMBz8IWn2MwGFSwkLc8C3mrbZe+j70D/6SPljwJ+7ZvTA56pctVYkMCkOsZE+Ll\nv3Otdq2dr7NH92ryqjMP7b/096+1ZflMtej6nLq/+HYmL2uc/Hf8vazH9j1yWXPa/AEQYHKEjj2f\n1da1y+S/d5sk6cTh/WrXfUCq/ePjYnX+9AmN/WparqtVYqJJm9cslo2NrV763wdsPABytSvnT2rX\n2vnat2mZYqMj0vy8zcv+VEJcrLYsn5lpAebK+ZPauWae9m9enq5lzSnzB5A9MIQsHQaNGH3vH8ua\nxYqKCE+1785Nq1WnUQu55Mmb6+q0ceVCBV8LVPenB6tM+SpsOABytXk/f6yb1y7JxsaQrue17D5I\nDk7OatHt+Uxb1jk/faiQoMsyGAxWqZW15w8ge+AITDpUql5X5avU1JkThxUfF6vVi2er76ARqezE\nz1evgcNzXY1ioiP1589fqmiJMnph5HtsNAByvbe+niNJOn14l354f1Can9f9hbfU/YW3MnVZ3/l+\noSQp4NAO/fTBi5leK2vPH0D2wBGYdOr5zNDk6WV+vyspKTFFn7DbIbp84YxqNWiS6+rz6/cfK+x2\niN7//Gc5ObuwwQDA30qUq5qNlrXaE3ndG1cvakTHclabP4CcgSMw6dSsXVdN+26iQm/d0M3rV7Vz\n02o1bdPlgT5b1ixRs7ZdZWNjmyHzzC6XUd61eY1WLZyl/439SuUqVWdjAYD7ODhmnx91nJzzPJHX\n9d+xxqrzB5AzcAQmvYnPzl5d+t4bj7zkr5Qn6W9YOV+tO+euSwdfPn9aX304Sl37v6iOvZ5lQwGA\nf//DtbXN1ctqNpu1b9OSHFcrAFbYH6cE6de5z/Oa89sPMiYk6Lj/Pp0LOCafSncPdwdePCuT0aiy\nFTLu5PWsfhnl0Fs3NPa1garXpJVGjP6YDQQALMhOJ6Y/iWXdtXaegq9cyHG1ApD5OALzH+TzKKgW\n7Xsk/734vqMwG1csUKtOvXJNLSLD72jMiKdVonQ5vfPJTzLYsEkBAB505dwJLZz2OYUAQICxph73\nncy/Ze0Shd0Okdls1pa1S9SqY+9cUYPoqAiNGfm0XPO6a9y302VnZ8+GAQB4wJHd6/XD+4MUFxNF\nMQBkCIaQ/Uc+Fauqau0GOn5or0xGo1bMn6Ga9ZvIy7uYPAsXyRXh5f1XBshgMOjjH2fK0dHpof3N\nZrOCr12Wd7FSbDwA8AjGhHgd2bVeO9fO15kju9N1J/rw0BBtXT5TJw9uV0hQoIwJcXL38FR+T29V\nqddC9Vp0lYfXf/s/ZTab5b9jjXavW6jLZ4/JmBCnQsXKqHmXgWrULuWPdz+OGaRT/rtSPP7vK5G9\n+vFvqly3WYbPP7XXOHlwmw5sWaELAf66cytYBhsbueUrqOI+lVW9YRvVad451R/lUruKmu/qs4/V\n935XLwRo05I/dPboXoWHhsjJOY+K+1RR9YatVb1ha127dFq+44Y+9HXiYqK0a90Cndi/VUGXzyg6\nPEz2jk4qWLi4KtZqrOZdn/vP2wFAgMnGejw9RMcP7ZUkrVzwp4KvBapN5z45fr0jwsP0/iv9ZZBB\nn0+dL+c8ro98zozJX+iY/x5989sSNhwASMWV8ye1e90C7du0TDFR4el+/q51C+Q3eYKKlamovsPH\nqmS5aoqOuKM9GxZpxawfdO74AS39/esUz/PwKqpPZmx56GtH3rmt3798M0UgCTx7XDO/e09Bl8+o\n99D3H2gb9ekMizvzj9qBz6j5/9u1i6c067sxSkpKVOveL6nH4HeUJ6+7wm4F6/i+LVo3b4r8d6zV\nqtk/acCrE1SpdsrbIfzfilOKuHNLa+b+rG0rZj90fj8uP6mwkOtaN+8X7Vzjl6b13LTkD63+6//U\nYcAIdRk4Sm4eXoqJvKOzx/Zpzdwp8vOd8MjX2LN+kRZO+1TlqjVQ1+dfV9HSFRUdeUeHtq/Wylk/\nav2CX7V5yQz1fOldtewxiA8esh2GkD2Gp1p1lFfhopLu3vtl+/rlatq2S45e51s3g/XW4B6ytbXT\n51PmyTWvW4o+ZrNZxoQEhd0O0YFdWzT+9Rc057cf1LZrPzYaAHiIOT+NU9ClM0pMNKb7udtWzNas\n796XZ5GSev2LWSpbuY7s7B3kXsBL7fsP18DXP0vuazDY6APf5fJdfVa+q88+MrwEXzmvL1/vLXcP\nL7351V/6btERfTl3r3oNeTf5imEbF01X4NnjT6QuGTH/M0f26Ou3+suzSEmN/n6B6rfsJncPT9nZ\nO8jTu4Radn9eH05drTKVayvkeqAmfzhEu9ctTLnjZGurfAUKqcfg0Y9cbjs7e3l6l1CPF99O03oe\n37dFC6ZO0ovvfqfWvQbLw6uo7Ozs5ZbfU3Wadda7Py565JGqpb9/rT+/fVctug/Syx9OVqkKNWTv\n4Kh8BQqpVY8X9Pa38+To7CKTyaj5Uz/R+vnT+OAh2+EIzOOkPxtbde3/on774RNJUoPm7eSSJ2+O\nXd+rly/o/Vf66+b1q5Kk3s0rpvm5jo5Oata2KxsNADzEO98vkCSdPLhd/zd2cJqfd+PqRc2fevd/\nUYcBw2Xv4JiiT4PWPbRzjZ/OHT8gszlJW5fP1jOj0nblyCkTXtEzoyapUu2n7n2vO7uoTe8hiomK\n1Jq5vskhauAbn2V4XR53/rdvXNPUj0coT958euZ/n6Q6PMzF1V3DP5qiSSO6KPz2Tc36fowKFCqq\n8jUapuibnnvV5HHLl6Z+a+dNlSSVLF891UD0wttf6aOX2ig2OjJF+54Ni7R23lSVr95AnZ99zeJr\neJfwUcenR2rJ9K/uBp4/vlG1hq1VuHgZPoDIPvvglODxdOo1UI5OzpKk1p1y7sn7AccO6s0XuiaH\nl/Rq3KpTjg53AJCRSqbzTvQbFv6qRJNJkuRTtV7q38Xt7t2j7MzRvWl+/eEf/fxAeLhf7aYdkqcv\nnTn2ROrxuPP38x2v2OhINWzT65HBw9Utv54eeXeYltmcpD+/fVcmkzFT3vcr505Iknasnpv68rl7\nqEmnASkeT4iL1aK/r/TWovvzD51P1XotkqeTkhK1ZekMPnQgwOQmrm7uatOlj9zzeajuUy1z7Hq+\n+3Jfhd8J/c/Pb9OlLxsLAKSRo0v67kQfcGhn8rRbvoKp9itVsWbydFhI2u8xVrh42VTb8nveOxE8\nPPTGE6nH48w/8OxxHd+3RZJUrUHa/k9Xa9hahYqVliSF3gzS/k3LMud9d3aRJC3741st/eMbGRPi\nLPZr2Lpnisf2blqqqIiwuyG2St2Hzqdg4eIP/H3q8G4+dMhWGEKWAUqVrShbW7scfRnhZbsv8EYD\nQCaxtU3fv+eIsJDkaZMpQQ62zhb75StY6IntdEtSfGxM5oe9R8x/78bF94Ud7zS9psFgUN3mXbRy\n9k+SpMO71qX5KmePw6dqXfnvWCuzOUlr/aZoz/qFatVzsJp07C/n+0YxeJcsl+JCCAGHtidPvzOg\nQbrme+dWMB86EGBym02rF2v42xMoBAAgQ6T3TvR53PIp/PZNSXcvo+zpXcJiPwcHp3TvzKcnbCVm\n0lCr9Mz/9JF7Q+Xy5M2f5tctU7lO8nTgueOZsi7dX3hbp/x3KzY6Ivm9XPzbF1r91/+pUfs+atl9\nUIqjJ/+4cv5k8vQPS49bPA8KyCkYQvaYrl+9pMjwUFWsVptiAACsoniZysnTF08dTrVfVMS9ocAV\nazXO8LBlNputGvYszf/W9cDkaaMxPs2vW+i++5ZFh4dlyrp4FS2ld75foFIVajzweFxstDYvmaHx\nL7XV71++qdCbQSmee/8yxkSG86EAAQap27hyoVrl4JP3AQBZX9POTydP71m/MNV+F0763/3nb2Or\nVrnk/h9JSYn3Alx42s/ldMl778ph9o+4WXNGKlSstEZ/N1+D3/1O3iXLpViX/ZuX65PhnXR0z8YH\n2kzGhOTpqxcD+FCAAAPLzElJ2rRqkVrngptXAgCyrmoNWqlu87v3ITvlv0sHtq5IuSOfmKh18+9e\nprfroDfkVbR0rqhN3nwFkqevnDuZ5ufZ2zs+ECoyk8FgUN0WXTT255Ua/tEUi0dkfvlk5ANH2+4P\nXCcPbONDAQIMLNu8Zomq12mkwkVLUAwAgFU99+Znati2lyTpz6/f0dp5UxUeGiKTyahLp4/qxw9e\n0OUzx9Xluf+pfb9huaYupe+78trpI2m/2lZcbHTydIWajVMJGk92N8pgMKh6w9Z65/sFGv7RFBUo\nXOyBQLpsxnfJf99/H5fd6xYqmmFkIMDkbrOmfqPezStqQJvqmvrNR7oRdEUHd2/VollTNeSNcRQI\nAGB19g5OatXjReX39FaVei20e91CjXuxld7qXVu/f/mmPL1L6gPfZer0zKu5qi61m3ZMnj64beUD\nQ60e5lZwYHKIaGDhssWSku8DJyld94pJSky0+PiIjuV089pFi23VG7bWB5OXq0zle+fcXjh5KHm6\nYq2nHghfs757P03nJB3ds1HfvfMsHyBkK1yFLA2WzPlVURF3f8lYNOsXLZr1i0qUKa+Pf5wp17xu\nFAgAYHV7NizSnJ/G6fm3vlCdZp2z3PIZDDYym5MkSYkmk2ztMmcXpOZT7VSoWGnduHpR0RF3tGXZ\nn2rTe8gjn/fP+UK1mnRI9S71bh5eivs7cNy5FZzqFcKCLp1+4O/4uJgHLot8v32blqrLc69bbHNy\ncdXTr07QpBFdJUkOjveGuTVu30er5/gm3zvmyO71mvHV23r29U9TvSLZtYunNPO797Lk9gI8DEdg\n0mD0xB9VonQ52dnbq1ipsho8aowmz1nH0DEAQJoZEx68Atb9J5c/Tl9J2rd5mf785l0ZE+Ll7JJX\niSbTYy3rv48Q/BM8LPn3r/ypHeFwcnFNnr4VfCXT5m9jY6vn3vxCNja2kqRVf03WjasXHz7/pETt\nXOMnV7f86vPymFT7FS977+pvh3eus9jn6O4NmjpxxAOPBZ49luprbloyQ6E3r6XaXrDwvX2P8jUa\n3gtT+T3VYcArKbaLiS+317YVs3U7+KpMJqOiIsJ0/uRBzZ/6ib58va9sbe1y3VE5ZH8cgUmDBs3a\nqkGzthQCAPCf/fvSt+G3b6Z6L5Z/971zK1geXkVTfe0l079Knv6/D1965LLkK1hYdZp1UueBo+Tk\nnCdFe9i/bmwYHRkuVzfL91CJiYr417Jfs3iBgCKlyuv8iQOSpGUzvtUzoz6RKSFeGxdP163rV/Ty\nh5Of2PzLVKql59/+Un9+/Y7iYqLk+9FQvfrxb/IsUtJiIFo47TOF3gzSyI9/Vb6ChVOtY9NOA3Rw\n20pJ0vI/v5M5KUn1WnaVi6u7Lp0+oi3LZuraxVMa+fFvGj/k3n7Er5+OUulKtRQTGa63v533wGvG\nxURp8riheu2T6RbnfWzfJkmSg6Ozuj7/xgNtHQa8ouAr57R/8/Lkx27fuKa5k8enGiqHjftZbvkL\n8gFFtmL73PC3x1MGZBWzpn6jzgNHUQgAOUaiyaSwm0Fa8Msk3bx26d6O9o1rKlG+mhydXZKPDiSa\nTAoLuf5333tHCW7duKpSFWrIwele3/sd2LJC4aE307xMcTFRuhjgryO71qtu8y5y+PsywcaEeN24\nckGLf/vigaMU4bdvqGjpinJ0cpGN7d35JyUl6s6tG1r06+e6dvH0fWHrhkpWqJ5iWV3z5tOBrXd3\n9oMDz2n9/F+0cdF0xcXG6KX3v5ejc54nOv+ipSuoRPmqOn/ioG4HX9We9QtlMiYob74CcnJ2UUJc\njM4e3as5P43T7RtXNWLCLyrhU+WhdSxQuJiiI+/o0umjSkpM1Cn/ndq4+Het8ftZezYsVoFCxTTy\nk+nKX7CwVs7+Kfl5dvaO8qlaVx2fHiF3D6/kx//pExUeql3rFirRZLq7fC6uigwP1b6NSzR/6iTl\nzVdAIyb88sARIOnu+To1G7eTjY2tLgQceuiRu8LFy2jExF9VslzVXPeZdLS31aI/vtdzw9/mCyqb\nMqz1v26mDMgq2tfylu/qsxQCQI4xomO5R/b553svPX3vFxMVrvlTJmnvxsXpXr6mnQbo6dc+zpRl\n9d+xRitn/59Cgi7Jw6uI6rXsptY9B8vR2SXTapUQF6v9W5bryK71Crp8RhFht2RjY6v8noVVulIt\n1WvRTZVqP5Xm+pnNZu3duFibl/6pG1cuyMHRSSXLV1fj9n1Uq0mH5H4jO5VXuWr11bh9X9Vq0sHi\neSkjOpbT65/PVL6ChXTiwHad8t+haxfPKCIsRHZ29ipSqrxqN+2opzr0T65ZakJvBmnH6rk6cWCb\nbgdfVXxctFzd8qu4TxXVbtJRdVt2lZ2dfa78TLq52Gtg81Ja63+dLygCDECAAQBrmf3DB4qOvKOX\nx05OtY8xIU4xkRG6ffOaLpw8pI2LfpMkfTZ7FwUEAQbZBifxAwCQzR3ctlI718xT9YatH9rP3sFJ\n7gW8VKZSLbXp/ZKGj5+q6Ig7FBAAAQYAAGSO6Mjw5JO0S1esla7nFilZPtULCQAAAQYAAGS4dfOm\nJB9Fye9ZOF3PvXL+pKrWb0ERARBgAABA5rj/krmhN9M+pt9sNmvdvKlq22coRQRAgAEAAJnDZDIm\nT0///HXdDr6apuf89eNY1W7W6aH3OQEAAgwAAMhQPQe/I4PBIEm6eiFAE15urz+/fVeHd61TyPVA\nJcTFKikpUXExUbp64ZQ2L5mhz1/rKZ+qdVW/ZTcKCCDbsaMEAABkX43a9ZabR0HN852okOuBMhkT\ntGf9Iu1Zv8hi/1IVaujZ/32i0hVrUjwABBgAAJD5qtRtrnG/rNXx/Zt1eOd6XT5zVHdu31BCfKzy\n5M0ndw9P+VStr+oNW6tircYUDAABBgAAWJetnZ1qNGqrGo3aUgwAORrnwAAAAAAgwAAAAAAAAQYA\nAABArsU5MAAAIMsyGKSCeRNVOH+iCuRNlLODWY72ZtnZmpVgNCg2waDoeBvdDLfV9TBbRcfx2yxA\ngAEAAMhkBfImqkrxBPl4G+Vob07z88KibHTyqoNOX3NQvNFAIQECDAAAwJOT3zVJTSvHqqiH6T8/\n/6mKcWpYPk5HLjnq0HlHGRMJMgABBgAAIAPZ2EgNysWpeql42WRA3rC1kWqXiVf5IkZtO+GsyyHs\n8gA55vuCEgAAAGtycjCra91o1SydMeHlfq5OSepYO1q1ysRTaCCH4OcIAABgNa5OSepeP1puLklP\nbB4Gg9SwfJzyOidp2wlnig5kcxyBAQAAVuFgZ1anOjFPNLzcr0rxBI7EAAQYAACA/6Z19VgVyJuY\nqfNsWD5OZQoZKT5AgAEAAEg7H2+jSnlZJ0g0qxKbrkszAyDAAACAXMzBzqynKsZabf7ODmY1LB/H\nGwEQYAAAAB6tcvEEuTha9whIxWLWXwYABBgAAJDFGQx3T6a3+g6QQapSnBP6AQIMAADAQxT1MGXa\nVccepVIxTuYHCDAAAAAPCzAFTFlmWfI4JSlfniTeFIAAAwAAYJl3/sQstTxFPEy8KQABBgAAwLKC\nblkrwBTMm8ibAhBgAAAAUrKzNcveNmtd+cuZK5EBBBgAAACLYcEh64UFJ25oCRBgAAAALLHNgnsd\ndrYEGIAAAwAAYEG80ZDllikuwcAbAxBgAAAALISFLBhgYhPYFQIIMAAAABaYzVJkbNba9YiK4wgM\nQIABAACw4OKpw9q4fmeWWqbrYXa8MQABBgAA4J6oiDDN+n6Mvn6zv1av3JxllivJLAWH2fIGAQQY\nAAAAKSkpUdtWzNZHg9to19r5MpuTtHz5NplMWePmkVdv28mYyBAygAADAAByvQsB/vpiVC/NnTxe\nsdERyY9fv35LixZtyhLLePyyI28UkA0x8BMAAGSYyDu3tXj6l9q7YbHM5pT3WLGzs9eeYzHq18+6\nyxkeY6PAW+wGAQQYAACQKyUlJmrbytla/ucPDxxxuV+FGo3Ub8Q4eZfw0bnrMfLxNlpteXcEOMvM\nPSwBAgwAAMh9zp84ID/fCbp64ZTF9nwFC6v30PdVp1mn5Md2nnJW8YImOdpnfoo4H2yvwBB2gQAC\nDAAAyFXCQ0O0ZPqX2rdpaarDxVr1fFEdnx4pR2eXB9pi4g3adMxZHWrFyJCJ59GHx9ho+0ln3jyA\nAAMAAHKLpMREbVk+Uytn/ZTqcLFKtZ9Sv1c+UqFipVN9nUs37bU9wFnNKsdmynLHJhi08kAexSZw\n5TGAAAMAAHKFs8f2yc93goIunbHY7uFVRL2HjlGtJu3T9HonAh1kazCrccW4J3okJjLWRqsPuSg8\nhguwAgQYAACQ44WHhmjRr5/pwJYVqQ4Xa9NniDr0f0UOTukbonX0sqPuxNiqbY0YOdhl/DkxwXds\nteYQR14AAgwAAMjxkhITtWnJH1o1+yfFxUZb7FO5bjP1f2WcPIuU/M/zCQyx05ztedW4YqzKZdDV\nyeKNBu0/56QTgQ5K4opjAAEGAADkbGeO7JGf7wRdDzxnsd3Dq6j6DBujmo3bZcj8YuIN2nDERSev\nmFS9ZIJKehll8x8OmsQmGHTqmoOOXHTkqAtAgAEAADndnds3tGja5zqwdYXFdnsHR7XpPUQdBgyX\nvYNThs8/KNROQaF2yuOYpHJFjPLOn6hC+Uxydkj9MEpYlI2C79jpyi07Xbxpr6Qk3keAAAMAAHI0\nk8mozUtmaNVfPyk+NsZin6r1W6jvsLGPNVwsraLjbXT4oqMOX/x7+aKva/onQ+Th4a7RX81QnNGg\nOKONYuINijdypAUgwAAAgFzj9OFd8vOdqOAr5y22FyxcXH2GfaDqDVtbbRmj4mx0+PDdq591u2KS\ngyP3cwEIMAAAIFcJvRmkRb9+rkPbV1tst3dwUru+Q9Wu38tPZLhYejg63bsZZqLJJDny/gEEGAAA\nkCuYTEZtXDRda+b6pjpcrFqDVuo7fKwKFi5OwQAQYAAAgHUEHNohP9+JunntosV2zyIl1XfYWFWt\n34JiASDAAAAA6wi9GaQFv0zS4Z3rLLbbOzipff/hattniOwdst74LFu7e7st8XExcs6TlzcVIMAA\nAICcxmQyasOCaVo952cZE+Is9qnRqK36Dv9AHl5Fs+x63H/SflJiIm8sQIABAAA5zYkDWzXv548V\nEnTZYrtnkZLq/8o4Va7bLFuszz/n49jasgsDEGAAAECOcfvGNS2YOklHdq+32O7g6KwOA15Rmz5D\nZGdnn23Wa+Lvm3hzAQIMAADIKYwJcVq/4Fet9Zua6nCxWk3aq/fQMfLwKkLBABBgAACAdRzbu1kL\nfpmU6nCxQsVKq98r41SpdhOKBYAAAwAArCPkeqAWTJ2kY3stD69ydHZRx6dHqlXPF7PVcDEAIMAA\nAJCDGBPitG7eL1o3f1qqw8XqNOukXkPeU35PbwoGgAADAACs4+juDZr/yyTdDr5qsd27hI/6vfKh\nKtRsTLEAEGAAAIB1hARd1vwpn+j4/i0W2x2dXdT52dfUovsghosBIMAAAADrMCbEac3cKVq/YJpM\nxgSLfeq26KLeQ96XewEvCgaAAAMAAKzj8K51WjB1kkJvBlls9y7ho/4jx6t89QY5vhYrZn6vfZuW\nqX6rbury3OtsHAABBgAAZBU3r13SvJ8n6uTB7RbbnZzzqPPAUWrZfZBsbG1zRU2iwsN0K/iKosLD\n2EAAAgwAAMgKEuJitXqurzYu/E0mkzFFu8FgUL0WXdVzyHty9/DMVbVp12+YGrfvK1d3DzYUgAAD\nAACszX/HGi2c9lmqw8WKlKqgASM/kk/VermyPh5eReThVYQNBSDAAAAAawq+ckHzfp6gU/67LLY7\n53FTl4Gj1LzrwFwzXAwACDAAAGQx8bExWjXn/7R58R+pDher36qHegweneuGiwEAAQYAgCzk4LaV\nWjjtc925FWyxvViZSuo/8iOVrVyHYgEAAQYAAOu4HnhPY98aAAAgAElEQVROfr4TdObIHovtznnc\n1PX519Ws8zMMFwMAAgwAANYRHxujlbN/1OalM5RoMqVoNxgMati2l3q8OFp58xWgYABAgAEAIPOZ\nzWYd3LpSC6d9qvDQEIt9ipetrAGvTlDpijUpGAAQYAAAsI7rl8/eHS52dK/FdhdXd3Ub9Iaadn5a\nBoMNBXuEI7vX62LAYZWuVFM1GrWlIAABBgAAZITY6EitnP2Ttiz7U0mJiSnaDQYbNWrXWz0Gj5ar\nW34KlkYBB3do28q/1KzzMwQYgAADAAAel9ls1v7Ny7T4ty9SHS5Wslw19R/5kUpVqEHB0qlSnSZy\ncnFV6UoMtQMIMAAA4LFcu3hKfr4TdO74AYvtrm751XXQm2rSsR/Dxf6jGo3acuQFIMAAAIDHERsd\noRUzf9DWFbNTHS72VId+6vbCmwwXAwACDAAA1mE2m7V342Itmf6VIsJuWexTqkJ19R85XiXLVaNg\nAECAAQDAOq6cPyk/3wm6cPKQxXZXt/zq/uLbaty+D8PFAIAAAwCAdcRGR2jZjO+0fdWcVIeLNe00\nQN1eeFMuru4UDAAIMAAAZD6z2azd6xdqyW9fKioizGKf0hVrasCrE1S8bGUKBgAEGAAArOPKuROa\nO3m8Lp46bLHd1d1DPQaPVqO2vWUwGCgYABBgAADIfNGR4Vo+41ttXzVXZnNSinYbW1s16/yMuj7/\nupzzuFGwTBB6M0hR4aFydfeQh1cRCgIQYAAAgNmcpF1r52vp79+kOlysbOU66jdiHMPFMtm6eVO1\nbeVfatb5GQ14dQIFAQgwAADkbpfPHJWf7wRdOn3UYrtbfk/1GPy2GrTuyXAxACDAAABgHVERYVr2\nx7fauWZeqsPFmncdqC4DRzFcDAAIMAAAWIfZnKQdq+dp2YxvFB1xx2Ifn6r11H/EOBUtXZGCAQAB\nBgAA67h46rD8fCco8Oxxi+3uHp7q+dK7qteyG8PFAIAAAwCAdURFhGnJ9K+0e93CVIeLtej2vLoM\nHCUnF1cKBgAEGAAAMp/ZnKTtK+do6R/fKjY6wmKf8tUbqP+Ij+RdshwFAwACDAAA1nEhwF9+k8fr\nyvmTFtvdPTzVe+gY1W3RhWIBAAEGAADriLxzW4unf6m9GxbLbDanaLe1s1PL7i+o87OvydHZhYJl\ncb2GvqduL7wpO3sHigEQYAAAyDmSEhO1beVsLf/zh1SHi1Wo0Uj9RoyTdwkfCpZNODg6y8HRmUIA\nBBgAAHKO8ycOyM93oq5eCLDYnq9gYfUe+r7qNOtEsQCAAAMAgHWEh4ZoyfQvtW/TUovDxezs7NWq\n54vq+PRIhosBAAEGAADrSEpM1JblM7Vy1k+pDherVPsp9R0+ToWLl6FgAECAAQDAOs4e2yc/34kK\nunTaYruHVxH1HjpGtZq0p1gAQIABAMA6wkNDtOjXz3Rgy4pUh4u16TNEHfq/IgcnTvoGAAIMAABW\nkJSYqE1L/tCq2T8pLjbaYp/KdZup3/AP5VW0FAUDAAIMAADWceboXvlNHq/rgecstnt4FVWfYWNU\ns3E7ipVD7VwzTwGHdqhS7SZ6qkM/CgIQYAAAyHrCb9/Uwmmf6cDWFRbb7R0c1ab3EHUYMFz2Dk4U\nLAe7cu6EDm1fLVe3/BQDIMAAAJC1mExGbV4yQ6v++knxsTEW+1St10J9h4+VZ5GSFCwXqNeyq4qV\nrcTNRwECDAAAWcvpw7vk5ztRwVfOW2wvULiY+g4bq+oNW1OsXKRslboqW6UuhQAIMAAAZA1hIde1\ncNpnOrR9tcV2ewcntes7VO36vcxwMQAgwAAAYB0mk1EbF03Xmrm+qQ4Xq9aglfoOH6uChYtTMAAg\nwAAAYB0Bh3bIz3eibl67aLHds0hJ9R02VlXrt6BYAECAAQDAOkJvBmnBL5N0eOc6i+32Dk5q33+4\n2vZ5ieFiAAACDADAOkwmozYsmKY1c6coIT7WYp8ajdqq7/AP5OFVlIIBAAgwAADrOHFgq+b9/LFC\ngi5bbPcsUlL9XxmnynWbUSwAAAEGAGAdt29c04Kpk3Rk93qL7Q6OzuowYLja9BkqOzt7CgaLrpw7\noZDrgfL0LqHiPlUoCECAAQAgYxkT4rR+wW9a6zdFxoQ4i31qNWmv3kPHyMOrCAXDQ+1cM0/bVv6l\nZp2f0YBXJ1AQgAADAEDGObZ3sxb8MinV4WKFipVWv1fGqVLtJhQLaVKoeFlVrNVYhYqXpRgAAQYA\ngIwRcj1QC6ZO0rG9myy2Ozq7qMOAEWrdazDDxZAuLbs/r5bdn6cQAAEGAIDHZ0yI07p5v2jd/Gmp\nDher06yTeg15T/k9vSkYAIAAAwCwjqN7Nmr+1E90O/iqxXbvEj7q98qHqlCzMcUCABBgAADWERJ0\nWfOnfKLj+7dYbHd0dlGnZ15Tyx6DGC4GACDAAACsw5gQpzVzp2j9gmkyGRMs9qnbvIt6D31f7gW8\nKBgAgAADALCOw7vWacHUSQq9GWSx3buEj/qPHK/y1RtQLAAAAQYAYB03r13SvCkf6+SBbRbbnZzz\nqPPAUWrZfZBsbG0pGACAAAMAyHwJcbFaPddXGxf+JpPJmKLdYDCoXouu6jnkPbl7eFIwAAABBgBg\nHf471mjhtM9SHS5WpFQF9R8xTuWq1adYeOJCbwZp7KDmkiTf1WcpCECAAQDgruArFzR/ykQFHNpp\nsd05j5u6DByl5l0HMlwMAECAAQBYR3xsjFbN+T9tXvxHqsPF6rfqoR6DRzNcDFaVaDLJ1o5dGYAA\nAwDItQ5uW6mF0z7XnVvBFtuLlamk/iPGqWyVuhQLVuHg5Jw8nRAfK2e7vBQFIMAAAHKb64HnNM93\nok4f2W2x3TmPm7o+/7qadX6G4WKwKltbdl0AAgwAINeKj43Rytk/avPSGUo0mVK0GwwGNWzbSz1e\nHK28+QpQMAAAAQYAkPnMZrMObl2phdM+VXhoiMU+xctWVv+R41WmUi0KBgAgwAAArOP65bPy852g\nM0f3Wmx3cXVXt0FvqGnnp2Uw2FAwAAABBgCQ+eJiorRi1o/asuxPJSUmpmg3GGzUqF1v9Rg8Wq5u\n+SkYsiRHJxeN+nRG8jQAAgwAIIcxm83av3mZFv/2RarDxUqUq6oBI8erVIUaFAxZmo2trSrWakwh\nAAIMACAnunbxlPx8J+jc8QMW2/O45VO3QW+pScd+DBcDABBgAADWERsdoRUzf9DWFbNTHS72VId+\n6vbCmwwXAwAQYAAA1mE2m7V342Itmf61IsIsDxcrVaG6+o8cr5LlqlEwAAABBgBgHVfOn5Sf7wRd\nOHnIYrurW351f/FtNW7fh+FiAAACDADAOmKjI7RsxnfavmpOqsPFmnYaoK6D3lSevO4UDABAgAEA\nZD6z2azd6xdqyfSvFBUearFP6Yo1NeDVCSpetjIFAwAQYAAA1nHl3AnNnTxeF08dttju6u6hHoNH\nq1Hb3jIYDBQMOUrAoR26cv6kipetrEq1m1AQgAADAMiqoiPDtXzGt9q+aq7M5qQU7Ta2tmra6Wl1\nG/SGnPO4UTDkSEd2rde2lX+pWednCDAAAQYAkBWZzUnatXaBlv7+taIiwiz2KVu5jvqNGMdwMeR4\nZavUVVJSkspWqUsxAAIMACCruXzmqPx8J+jS6aMW293ye6rH4LfVoHVPhoshV6jXsqvqtexKIQAC\nDAAgK4mKCNOyP77VzjXzUh0u1rzLs+ry3P8YLgYAIMAAAKzDbE7SjtXztGzGN4qOuGOxj0/Veuo/\nYpyKlq5IwQAABBgAgHVcOn1EcyePV+DZ4xbb3T081fOld1WvZTeGiwEACDAAAOuIigjTkulfafe6\nhakOF2vR7Xl1GThKTi6uFAwAQIABAGQ+szlJ21fO0dI/vlVsdITFPuWrN1D/ER/Ju2Q5CgYAIMBQ\nAgCwjgsB/vKbPF5Xzp+02O7u4aneQ8eoTvPODBcDAIAAAwDWEXnntpb8/pX2rF8ks9mcot3Wzk4t\nu7+gzs++JkdnFwoG/EtURJjiYqLk5OIqV7f8FAQgwAAAnoSkxERtWzlby//8IdXhYhVqNFK/EePk\nXcKHggGpWPHn99q28i816/yMBrw6gYIABBgAQEY7f+KA/Hwn6uqFAIvt+QoWVu+h76lOs84UCwAA\nAgwAWEd4aIiWTP9K+zYtsThczM7OXq16vqiOT49kuBiQRgNencCRF4AAAwDISEmJidqyfKZWzvop\n1eFilWo/pb7Dx6lw8TIUDAAAAgwAWMfZY/vk5ztRQZdOW2z38Cqi3kPfV60mHSgWAAAEGACwjvDQ\nEC3+9XPt37I81eFibfoMUYf+r8jByZmCAQBAgAGAzJeUmKhNS/7Qqtk/KS422mKfynWbqd/wD+VV\ntBQFAwCAAAMA1nHm6F75TR6v64HnLLZ7eBVRn2EfqGbjdhQLAAACDABYR/jtm1o47TMd2LrCYru9\ng6Pa9B6iDgOGy97BiYIBAECAAYDMZzIZtXnJDK366yfFx8ZY7FO1Xgv1HT5WnkVKUjAAAAgwAGAd\npw/v0ryfP051uFiBwsXUd9hYVW/YmmIBT9DaeVPlv2ONajXpoPb9hlEQgAADALhfWMh1LZz2mQ5t\nX22x3d7BSe36DlW7fi8zXAzIjM/kzSAFnj2uUuWrUwyAAAMA+IfJZNTGRdO1Zq5vqsPFqjVopT7D\nPpCndwkKBgAAAQYArCPg0A7N+3mibly9aLHds0hJ9R02VlXrt6BYAAAQYADAOkJvBmnBL5N0eOc6\ni+32Dk5q33+42vZ5ieFiAAAQYADAOkwmozYsmKY1c6coIT7WYp8ajdqqz7APVKBQUQoGAAABBgCs\n4+SBbfL7eaJCgi5bbPcsUlL9XxmnynWbUSwAAAgwAGAdt29c04Kpk3Rk93qL7Q6OzuowYLja9Bkq\nOzt7CgYAAAEGADKfMSFO6xf8prV+U2RMiLPYp+ZT7dTn5Q/k4VWEggFZTI3GbeVRqKiKl61MMQAC\nDADkbMf3bdH8qZ+kOlysULHS6vfKOFWq3YRiAVlUpdpN+IwCBBgAyNlCrgdqwdRJOrZ3k8V2R2cX\ndRgwQq17DWa4GAAABBgAsA5jQpzWzftF6+ZPS3W4WO2mHdV76PvK7+lNwQAAIMAAgHUc3bNR86d+\notvBVy22e5fwUb9XPlSFmo0pFgAABBgAsI6QoMuaP+UTHd+/xWK7o7OLOj3zmlr2GMRwMQAACDAA\nYB3GhDitmTtFGxb+KmNCvMU+dZt3Ue+h78u9gBcFAwCAAAMA1nF41zotmDpJoTeDLLZ7l/BR/5Hj\nVb56A4oFAAABBgCs4+a1S5o35WOdPLDNYruTcx51evY1terxgmxsbSkYkAMEX7mgO7eCla9gYRUu\nXoaCAAQYAMj6EuJitcbvZ21Y8KtMJmOKdoPBoHotuqrnkPfk7uFJwYAcZMvSGdq28i816/yMBrw6\ngYIABBgAyNr8d6zRwmmfpTpcrEipCuo/YpzKVatPsYAcyNU9vwoWLi5X9/wUAyDAAEDWFXzlguZP\nmaiAQzsttjvncVPnga+pRdfnGC4G5GBdnntdXZ57nUIABBgAyJriY2O0es5kbVr8e6rDxeq36qEe\ng0czXAwAAAIMYFlCQrzWLJqtTasX69K5AMlgUNHipdS0bVd17v288rrno0h4bAe3rdTCaZ/rzq1g\ni+3FylRS/xHjVLZKXYoFAAABBrDs3Knj+mT0UDk6OWnQiHdVs/5TMsigi+cCtODPn7V0znSN/vgH\n1W7YnGLhP7keeE7zfCfq9JHdFtud87ip6/P/U7POzzJcDAAAAgyQuqMHd2vsyGdUvkpNTfKdI0dH\np+S2yjXqadw39TTl63H6cNRzGvf1b2rQrC1FQ5rFx8Zo5eyftHnpH0o0mVK0GwwGNWzbSz1eHK28\n+QpQMAAACDBA6oKvBeqj/w2SvYOjPvz61wfCy/2GvjFOh/Zs06R3Xpbv3A0qVqosxcNDmc1mHdy6\nUgunfarw0BCLfYqXraz+I8erTKVaFAwAgFzKhhIgPb4Z/4ZioiPV67lhcs+f+q/ftrZ2enboG4qP\nj9NX40YpKSmR4iFV1y+f1Q/vPafpX7xhMbw453HTgJHj9d5PiwkvAADkchyBQZrt27FRRw/sksFg\nUPtuAx7Zv3HLjnLO46pTxw5p06pFatOlL0XEA+JiorRi1o/asuxPJSWmDLkGg40ateutHoNHy9WN\n+z0AAACOwCAdFs2cKkkqUaa8ChbyfmR/ewcH1WvcUpK0zO93CohkZrNZ+zYt1YSh7bRp8e8Ww0uJ\nclU1+rt5Gvj6p4QXAA9YMfN7jXuxlVbM/J5iALkQR2CQJrduXNfh/TskSRWr1U7z8ypWq61t65fr\n9HF/XTwboNLlKlHMXO7axVPy852oc8f3W2zP45ZP3Qa9pSYd+8lg4DcWAClFhYfpVvAVRYWHUQyA\nAANYdnD3FpnNZklS0RKl0/y8+wPLsYO7CTC5WGx0hFbM/EFbV8xOdbjYUx36qdsLb3LEBcBDtes3\nTI3b95WruwfFAAgw+LdP3xuurWuXptre85mhGj564gOPzZzytWZN/eaBx5q3764xn0/5T8vQvpZ3\nhq/XWv/r6ep//PC+5GnvYqXS/DzPwkWTpwOOHVS3AYPZqHIZs9msvRsXa8n0rxURZvnqYqUqVFf/\nER+pZPnqFAzAI3l4FZGHVxEKARBgYMmoMV+oas36+uW7CTImJCQ/XrdxS705/lsV8Cyc4jkDh72l\njj2f1fsjBuj61Ut6/cOv1bRt12xdh6uXzidPF0lHgMnv4Zk8fencaTaoXObK+ZOa5ztR508etNju\n6pZf3V98W43b92G4GAAAIMBkBFc3d3UbMFj2jo76fuLbyY83b9fNYniR7t5or2Ahb3kW8lbbLn0f\n++pb6T1a8iTcCLqSPO2WL+3Dexyc7t0nJiriDhtULhEbHaFlM77T9lVzUh0u1rTTAHUd9Kby5HWn\nYAAAIM34yTONOvZ8VrUaNEv++8Th/Q/tHx8Xq/OnT6hLv0E5Yv1joiOTp51c8qT5efZ29snTkZHh\nbEg5nNls1q51C/TRS221dfksi+GldMWaeu/HRRrw6gTCCwAAIMA8SYNGjE6e3rxmsaIiUt8h37lp\nteo0aiGXPHlzxLrHx8UlTzs7pz3AGE3Ge68RG8tGlINdOX9SX7/ZT7O+e19R4aEp2l3dPTTwjc/0\n9rfzVNynCgUDAAAEmCetUvW6Kl+l5t879LFavXh2qn03rpyv1p375Jh1t3dwSJ62tbNN8/NMxnvn\nDTk6O7MR5UDRkeGa+38f6fPXeuriqcMpv2RsbdW860BN+G29GrfrI4PBQNEAAAABJrP0fGZo8vQy\nv9+VlJRyiEzY7RBdvnBGtRo0yTHrff+RpPQcSYmPu9fXxcWVDSgHMZuTtHPNPE0Y0lbbVv4lszkp\nRZ8ylWvr3R8Wqf+Ij+Scx42iAQCAx8ZJ/OnUrF1XTftuokJv3dDN61e1c9NqNW3T5YE+W9YsUbO2\nXWVjY5sh88wKl1Eu4OmlsNs3JUlxcTFyzpO2MHIn9HbytJd3MTagHOLy2WPymzxel04ftdjult9T\nPQa/rQate3LEBUCGiwgL0abFf0iSegweTUGAXIYjMOlNfHb26tL3+eS/l/w1LUWfDSvnq3Xnvjlq\nvUv53LsBZciNtIef8LB7AaZYyTJsQNlcVESY/vrxQ335vz4Ww4uNra1adn9eH01bo4ZtehFeADwR\ncTFRWjf/F61fMI1iALlxf5wSpF/nPs9rzm8/yJiQoOP++3Qu4Jh8KlWTJAVePCuT0aiyFTLuJOWs\ncBnlcpWra8OK+ZKk4KuXVb5yjTQ9L+jKpeTpyjXqsfFkU2ZzknasnqdlM75RdCqXw/apWlf9R3yk\noqUrUjAAT3bnxd7x7+8mM8UAciGOwPwH+TwKqkX7Hsl/L77vKMzGFQvUqlOvHLfODZu1S56+cOZk\nmp8XeOFM8nStBk3ZeLKhS6eP6Iv/9dacnz60GF7cPTz1wuiv9caXfxFeAAAAASar6nHfyfxb1i5R\n2O0Qmc1mbVm7RK069s5x61u4aInkK7AdObArzc87feLuVakq16gn72Il2XCykaiIMM36foy+eqOf\nAs8eT/nlYWurVj1f1EfT1ql+q+4MFwMAAASYrMynYlVVrd1AkmQyGrVi/gwd998rL+9i8ixcJEeu\nc/8XX5UknTlx+IFzW1ITEx2p08f9JUm9Br7MRpNNmM1J2rZitsa/1Fa71s63eHWx8tUb6IPJy9Xn\n5TFy4upyADKZvYNj8nTsfTdaBpA7cA7MY+jx9BAdP7RXkrRywZ8KvhaoNjno3i//9lSrTvKpVE3n\nAo5p7dK56vfCyIf237pumRITTSpfuYaatO7MBpMNXAjwl9/k8bpy3vIwQXcPT/UeOkZ1mnfmiAsA\n6+282DtQBCAX4wjMY+3Qd5RX4aKS7t77Zfv65WratkuOXV+DwaA3xn0rO3t7LZ49TTEP+dUrISFe\nftP/T46OTnprwvfs7GZxkXdua+Z37+mbt/pbDC+2dnZq0/sljf91g+q26ML7CQAACDDZsng2tura\n/8Xkvxs0b/fADR9zIp+KVfXGuG8Udvumfvj4HYs38jSZjPp63P8UcuOa3v9iikr5cGJ3VpWUmKgt\ny2Zq/JB22r1uocUr+lSo0UhjJi9XryHvydHZhaIBsDpbWztVrNVY1Ru2lq0dg0mA3IZP/WPq1Gug\nZk39RvFxsWrdqXeuWOc2XfrKzt5e3098W+8N76+BL78pn0rVFBsTrSP7d2nOr9/LaDTqq2kLuXRy\nFnb+xAH5+U7U1QsBFtvzFSys3kPfU51mDP8DkLU4ODlr1KczKARAgMl+YqOj0nxH+CfF1c1dbbr0\n0Y4NK1X3qZa5ZsNp0b6HqtZqqEWzpuiHSe8o+FqgHByc5FOxqnoNHKY2XfrK3iFrjFGOi4nK8iea\nJ8THysHROVPmFR4aoiXTv9K+TUssHnGxs7NXy54vqNPTr3LEJYu8Z1l5PXJKHdhmWFcgu4qLjZFT\nLvt/na2HkG1ctTBLLEepshXVvH132dnZ56qNp6BXYb385nj9tniHVu4L1OIdZ/TVr4vUsdezWSa8\nhARd1tlj+7J8LQ9tX62E+NgnOo+kxERtXjJDE1/uoL0bF1sML5VqP6UxvivUc/A7hJeHiAoP1Yn9\nW7P9ekSEhejkwe2P9Rr7Ny9TUmIiG0Uaa2Xpqn45ze3gqzp/4hBvOJBJdm9dq+ioCAJMtgkwKxem\n6XK+T9qm1YvVOgdffSw7u3jqsC6eOpLll/P8iUMW77WSUc4e26dPX+2u+VM/UWx0yi85D68iGvrB\nT3pt0h8qXLwMG04atqtLp49k+/W4EOCvS4/5+Th9ZE+qwxDxoJMHdyg48HzO/3ycPqyLp/x5w4FM\ncvqYf/JtKwgwWZzJaNS5gKMKOHrQqstx/eolRYaHqmK12nyCsuiO5uUzWX9H8/KZI7p85liGv254\naIj++PItff/uQAVdOp2i3c7OXu37D9e4qWtUq0kHNpg0unT6aLYIxo8SePb4Y293l88c0+UnGL5z\n2vfRxVOHc8d6BhzmDQcyycmjBxRwLHcd9cy2Aeb86eNKSIjXKSu/YRtXLlSrXHLyfrbc0Tx1RJdO\nH83SwzbiY2MUdOmsLgRk3K8nSYmJ2rhouiYMaat9m5dZHC5WuW4zjZ2ySt1feEsOToxVT1+AOazA\ns8ez/dCpiwGHdenMEYtXE0yLqPBQhQRd5tf2NIgIC9GdW8EZ+jnPutuVvy6dOcKbDmSChIR4nT99\nXGdO5K4fDbLtSfynjt8NLtYMMOakJG1atUif+s7hE5QFGRPide3iKZlMRgUHnpd3yXJZc2f47x3I\njPpl9szRvfKbPF7XA89ZbPfwKqI+wz5Qzcbt2Ej+y+fenKRLp48qIT5WVy8EqES5qtlyPZKSEnXp\nzJHkAF2sTPovd/7PNnsxgADzKP8El5x+BMZkMurq+QCZTEaFBF2WZ5GSvPnAE3Qu4KhMRqPVRyRl\neoC5EXhO1y5mv/HLuzeukCQFHD2gg1uWWeXGeocP7lWRIkUUdO6Igs7xa1NGaNq0qa74r8uQ17oa\neEEmk1GSdHDdbNWo2zhLrvPhLWslSXduBevklnnK657vP71OVES4Nq5aqJOpfInZ2tmpYdM2atyi\ng+zslWF1zm1u3biu2L9v4uq/0U+GqKBsuR43r19TfGyMJOnQ+r9krt8k3a9xZMsySdKNa5d0Ztdi\nOTvnYQNJxdGtSyRJwVfO6/ze5XJwcMyR6xl09VLy9+6BtTNVtWZ93nxkOQ52NmratKkObV2e7ddl\n59YNkqTwsNtav/A3FSjolTsCTGToDV08eSDbLfj5M3fvFh4XFyv/XRuUP3/+JzavgwcP6vjx47Kx\nsZGPj4+qVq2q8PBw7du3T126dMmW9cuqKlWqpNuXMiYMnj16NHn6YsAhFSuYNXeuLp+5t5ynD21R\n6dKl0/X8pKQkHTt2TIcOHZLRaLTYp0SJEmrcuLHc3NwUfu0kG9pjOH3q1APbVclCbtlyPc7etx6X\nTvmrhFf6b8IbeO7E3QmzWWcOblGxYsXYQFL7nP99npA5KUmnD2xSkSJFcuR6nj1+73yoCycOyDvf\nkwtqO3bs0MmTJ1W5cmU1adKEjQzp3t/ICftv9x95Obhzo8qVK5cr3j87zzzZ7xezuLg4RUTcu5LS\njRs3nmiAOX78uOLj4yVJx44d07Fjx5Q/f3516NBBDlnkcsFI6ebNmxans/pypifAXLt2Tbt27VJY\nWJjF9rx586px48YqWZJhHBklJCQkW2xX6dnubty4ke7nm83mFNsuASb1Wt2/3dy4cSPHBpj7t6X/\nsl0B+O/f5SEhIbknwGT3L8h/3ryKFSs+sfm1bNlSe/bsUUREhNzc3FS+fHlVq1ZNtra2fHKyyYc6\nNDRURqNR9vZZ6149kZGRiomJSfcOcVRUlPbs2VhUcqQAACAASURBVKMLFy5Y/mDb2alGjRqqUaOG\n7Ozs2Bie0HZ1584dxcfHy9Ex+w0Huv979M6dO0pISEjXDzL/PIed1UcLDQ2VyWTKEcE3vd+7iYmJ\nT+x/ZcOGDVWvXj3+FyPXiomJUWRkZK74bkmxn1O4cHHVqdM0Wy10UFDIv3bmYp7oOtSpI/Xs+Syf\nlGwkIiJcUVG/JP9tNpuVL19hVaxYOUst5/79ex74+/bt26pVq7FsbCz/QzaZTNqwYbVWrVqm+Pg4\ni32qV6+lfv0GytPTiw0hg8XHx+nXX3994DE3N09VrVojW61HbGyM7tyZ9q/PRyFVqpT2CxLs3Ln1\nXzvpYapdu4lVzkfM6rZt25SiVtnt/25aREZGKjLy3vduUlKSChQoqrJly7ERAE+Av/+Bf323hKp6\n9YZZ7sfaJxJgfHyqyMenSrZa6Nmz/3zg7+vXr6lVqx7Kk8eVrRmSpNWrF6d4zMkpr7p2HZillvPA\ngQdPuE9ISFDp0lVVvXqdFH23bl2nsWNH6fz50xZfq1QpH02Y8J3atu3CBvCE7NixSUlJD16S29Ex\n621Xj96hXp/i0tr29nnStR5btz4YYKKjo1SlSn2VLVuBDeVfNm3a+K8d/QjVrNlExYuXylHruX59\nyhOinZ3ds93nA8guTpx48CJcJpNJJUpUUN0setGijJTt7gOTlJSkw4f3PfBYYmKi/P33sSUj2cGD\ne1I8dvjw/iy3nJa2238v57VrgRo6tI+efrq9xfDi5OSs0aMnatOmo4SXJ+zf3z2SdOjQnmy4Hvsf\n+/OREa+RWxw4sDtHbDf/5Xv34MHdbADAE2Lpe8TS55AAkwWcPRugiIjwXPHPAI8TDPZm+X+kRmOC\njlm4j9E/y240Juj/27vzqKautQ3gD/OMiAKCEyKIyBQGFbVWrVO9zkNxav0carW1tVrb2qvWq71W\nq71a517rgFPhqli1arVaRByqBZRRcURUBFERZYYwfH/QxIQkEDBAQp7fWqwVkp3knJN9znnfffbe\nZ926b/Hmm244fvyg3M94++0RiIi4jrlzv4Yxb0bZIAnn1at/yb1RqKYFmjU5hubl5eLmzUQeh5WQ\nnf0S9+7d1oogQ97vz8ZForqhqPFeWxqSNC6BkReY8iBJlXfqeDn3Qnny5DEePkxRm+VMTIyVO44l\nOvoSwsNPolcvd6xcuQgFBfkyZRwdnbF372/YseNQo+uGot6Bv2wS/OLFcyQn39L4QDMz8ylSUu4o\nnciVlpZqRVD+umJjI2W6HSqqS5qsrKxM7nn40aMHePLkMSsCkYrdvJmIvLxcpY7vTGDUIoGJrFFi\nQ9yp1XXHVhTA3L17ExMnDpIbTJqYmGL+/GWIiLiGt94axB+7Hj18mKIwENOkwD0l5Q4yM5++1noo\n2o+uX4+Tm3BrM3ndxwDg2rU4CIXFjWY9b99OUnjcZTcyItVTdKWlqnMVE5gGDfrknzjVrXWd1K+O\nVJUAN4S4uJrdQOsf/xiFc+eS8OmnC2FgwPsP1beqkt+4OM25ZF/VPqBs1wNF5UpKSuR2i2S9kVVc\nXITExFit2D94ZY6ofmMdbWg00JgEJi3tIY4dC5Xb71okKGgjkpLipebbJ+3x8mUWIiJO4dChYIVl\nzpw5gd9/P4KMjLQGW87c3BxcuhSBS5cilCrv5NQBISG/Y9u2g2jZsg1/6HpWUlKC69fjcPToAYVl\nzp37A5GRF5Cfn6fW63HtWmw163EaUVEXFV5FycnJxsWL4YiKuqjwM0JD9+DWretyu01pk+zsl7hw\n4UyVQcb+/Ttx584Njd5WQqEQiYkxOHYsVGGZs2d/x5Url1BYWKDS7w4O3oYZMwIRHLyNByrSGnl5\nubh8+RwuXTqrsMyRI/uQlBQvt6tvY6GTlqZ+o0+fP3+G2NgoxMVFIzY2CrGxkXj6VPmbpJmamsHD\nwwcCQWd4e/tDIOgMR0dn3p+gEcnPz0NiYgxiY6MQExOJuLhopfvvi7Ro0RICQWepetKkSVOVLmdx\ncRGuXYtDbGwU4uMr6vPt20lKByyGhobw8ekqtZyOjs6sAHWkvLwcKSl3pOpVYmKM0t2i9PT00KFD\nJwgEXcR1ys3Nq97n5C8vL0dy8i3xMTQuLhoJCVeVDiD19fXh7OyG1q0dYWJigry8XKSk3MW9e7eV\nrrvm5hbw9PSFj08XeHtX1N/GOl6rsLAA167Fird1TEwkkpNvKT25g4WFJby8/ODt7Q9v787w8emC\nVq3aqt16lpWVITn51t/n5SjExUUpHMunqF65unrA29sfPj5dIBB0hqurR61vthscvA0REafQq9cA\nTJjwPg9g1OgIhUJcuxYrPpbHx0fj1q3rSicmjTkebvAEJjc3BwkJVyUSlkg8eHBP5d9jaWn194nU\nX/xD2tu34t6hITtwUlK81EmzJjtwTTg6OkvVEw8PH5iamin13tLSUty5cwOxsZHi+lwX/dybNm0m\nVY8Fgs6ws3NgRamF9PRU8YlBlLBkZ79Q6XcYGhrB3d377yS0om45O3eErq6uStdDtPyiuqfq9VAF\na+vmEvW2IoC1sbHTqDpTUlKCW7euiY9HsbFRuHEjQeVX/ps3txXv315e/g2yrVJT70vUqYp6lZOT\nrdLvMDY2gbu7QFwvfHy6wMmpAxscSeuUlZX9HUNU7G+xsVG4di0OxcVFjIcbOoERCotx7VocYmIi\nxa3Rd+7caLBLXHZ29vDyetUS5O3tj6ZNm3EvUoMdWLK1ISEhRuU7sLJELeqS9UTUon7/fjLi4l4F\nvgkJVxUOYq37uuwgVY/r4mqSpsvKypS4qltxgsjISG+QZRFdnaj4vSpa3JW9OiFaj5iYSPE+0lDr\noQoODq2lglcvL39YWjZRi2UTXZETbeuaXpFTtVat2koEHZ3h5eWnsm2VmflUvG+Izs816fmg2gCr\nyd9XpF5deVbHK1JEr+PhwxRxL6PY2CgkJFxFbm5Og8XDov1NU+LhektgcnNzcOTI/xAd/WeDJy6K\nEhmBoDOGDx8LJ6cO3LMaKFiIirqI8+fDGvwEWtmrLh6d4ecXgLy8PJw/fxqxsVG4e/em2vRhNzEx\nhYeHD7y9/TFixHj4+nZlxULFYMfDh0P+btGKVXlf/NrS1dWFs3NH+Ph0QZ8+b2PQoFFVdjeLirqI\nX3/dj9jYKCQmxijddUcTiLo6+Pp2xahRE+Hh4dOgy1NYWIBff92Pc+dOIyYmEikpd9Tmfj9mZubi\nBFgV2+ry5XM4evQAYmMj66TFVxUNSP37D8HAgcOhp6fHAxppNKFQiBMnfkF4+Mkadyuva5JXRIcN\nC0Tnzj2YwMhLaCS7jsXE/FUvs4hZWlqJW7BEl9DYlUx9PXr0QKrFvC66MFRmZGQs02fUyamDwi4/\nOTnZiI+/Iu4CGRsbhdTU+3W+bQwMDODq6iFuMREIOqNDB/da9yfXFiUlJbh5M1GqVf3mzcR6mfyj\nTZt2UuNjPD19YW5uUeuT4M2bieL9IyYmErdvX9eISUwMDAzg7i6Q6iLVoUMntQ5Os7NfyHQ3TE9P\nrYdtZQh3d2+pbnfOzh3rbFuJuuxKXuGrqy67lTk6Oksdz2rShZdIUzVUPCw5Jk20z3Xs6KkxMYRa\nDeKXHrwf+dot8KLWaMm+fhzMr9kkB5GK6klNBpEqswO7unq89qDrZ8+eSE1CERcXjWfPntT683R1\ndeHk1EHi8m5neHgIYGRkzEqhAgUF+UhMjJHqYnbv3u3XanGXvLIrql/W1s3rZT0kx2e87pUDXV3d\n12od1NPTg7NzR6mxHO7u3jA0NNL4epORkY74+GipZDgrK1Ml20qU6Lq7ezf4tOn5+XniAEu0njWd\nNKWy+phEhUhTqToe1tHRgZNTB6kGfHd3AYyNTTR2G6nlLGSS0tIeSrQw/oWLF8MVlm3f3hU9evQR\n/0AuLp3YGq0FJFuiY2IiERZ2XOFNnCwsLNG//1DxSdPDwwcmJqb1spypqfelEpq//jqvsLW8WTMb\nBAS8CR+frvD29oeXlx8sLCz5Y9cjyRb32NgoXLoUgRcvnssta2BgiICAN6XGcqjLld3s7BdSCc3l\ny+eqXI9u3XpJBZb29q3E471EY28iIy8oTGpsbVsgIKCXSq4yaaKUlLviAbhxcdFVbis7OwepeuPp\n6QszM3ONWM8XL55LJTSXLkUonDjC2NhEaj058QhR7ePhiivAUfjrrwsKJwmytm4u3ucEgi4qHS/H\nBKaWBAJ7hcHp8eN/wcenC2u5llu/fjm++26h3NcCA/8Pa9fuVIvlHDWqFy5fPif3tdWrt2H8+Gn8\nMdXIli1rsHTpPLmv9e8/FLt2/aoR67F58yosWzZf7muDBo3E9u2/VPsZb7/tj/j4K3Jf+/HHEAwf\nPo4V5m/9+glw/Xqc3Ne2bg3F4MGjG8V6rl69FKtXL5H72ujR72LDhj2sDEQqNG7cAJw7d1rua99+\nuxFTpsxq1Ouvq2kL7OMjf1CygYEhPDwErNEEX98Aha95e3fWiOWs6jVqGH5+AbV6TZP2D2XXo6py\nrLs12VaNZ5KNqhoP2bBIpHr+/t0axTlJixIY+QdCDw9Bg/cTJvWpI4oGuKrTTq1oWSwsLOHi4sYf\nUs14evoq7JKqSUG7t7f/a6+HonK2ti0a7c0qa0sgkN9oYm/fqlFNIOPr21Xh+FJFDY9E9Dr7nPzj\nsJGRMdzcvJjAaMoP5ufXjbWZAFRMMers3FHmeRMTU3Tq5K1WJ3xFAaYqb3JIqmFkZCy3/ujq6mpU\nC7OJiSk6dHCXeV5PTw+enr6vFZSr0xVO9T9nNa4WUisra7m3IKiY1ZG9I4jqK4bw9PR97YmImMDU\nAUWt69pwuYxeL2jw8vJTq0kd7Owc5LbAsguO+pJ3nHFxcdOYgddVrYerq4fS6+Hk1AFWVtY8DivB\nxcVN7uBZgaDxdauSF1BVBFPsHUFUF40Gjo7OWnsc1rgExszMHB06dJJ5ni1/JEle31BFrcbqFkgy\nCNSsxFgTE055y1yTeqejoyM3WGXyLeckq6srN1lpjDeZ1Zb1JFLvGEI7eiRpZD+VygfJZs1s4OjY\nnjWZqgnQ1G+nlrdM7C+uvuR1FatqIKUmJfg1TT4ql9fT0+NgbYXnLOnGE319fXh7+2tFMMU6QVS/\nsY46NtYygVFw8mWrH1Xm7NxRpjuMOl7ZqHygad3aEc2b2/IHVFPt2rnIdJ3SxK5A7dq5yHRrqmmg\nWbm8vH2O5B97OnXyrrf7T9UneevF8zNR/R1bbG1boFWrtkxgNCXjZJcbqqxya7Ci8SYNrfK4HJ7s\n1VvlrlOKurSq/YG/UrcmS8smcie+qC6BkZx1inVX+XNWY7z6AlRcWZKcCIKz0hHVLTc3TxgZGUvE\nw9ozoZVGJjCVB83yEjVVFzSoaz/syjOjMRnXrHpV1ZTd6k7ySrZA0KXGM99VnnVKE7vS1ZfK3Zwb\nc7IneT5md1iiumVgYCjVaKBNMYRGJjCS05bq6upyAD8pCDS7asROLblsbMXWrABNk38vyeCytsmH\n5PsYrFZN8jzVmJM9yRZgNi4S1W8M4eXlrzXrrbE3mxAFDoqmqCSS7CKjzoGm6IBjYGAId3dv/nAa\nkMCIuk5pcoAmOf6qtl2aRHXXzMycN19VMshQdL+UxhhMsUGGqP7iYW2bSEVjExjRQZIHSFJE1P9a\nT09Prfuci1pjPTwEUn1ZST1JBqCafNWholuT82sdR0XHYU3uSld/56xufyd9fgrvWN8Y2Nu3gp2d\nPWelI6onosaomtzLqzHQ19wfrIv4xElUVcuEhUUTtZ7xR3RTQG0afKfp/P27obi4CLa2LTQ8qA4Q\nJzO1IZp1ig1J1fPwEMDQ0Egrxgr5+gbg3r3bnJWOqB60bu0IW9sWWjeGVmMTGBsbO7Ru7cibZFE1\nJ9Kuat/FUDSzFZNxzeHt3RmFhYWNYv8oLy+v/QlEXx9eXn6su0oQdRHVhmTPzy9AZrpxIqrLfa6b\n1tz/ReMTGADo2bMfXF09WHNJIR+fLrC0tFL75RQIOjMI1LDAv7i4SOPXQyDogtLSUhXUXTYkKRtk\naEOjm69vAJo2bcYfnKgeYx1tuxKuN2/ekiWauvCtW7fT+C4cVLeaNbNBkyZWat8aaGZmgY4dmYxr\nTr2yhampmcYff5o3t4G5uQVsbOxq/RlNmzZD+/YdWCmUYG3dXGo65ca8npaWTWBt3Zw/OlE9MDY2\ngZubZ6MeX1eZTlraa/QfICIiUmNFRYU4fvwgQkK2488/z+LRo7Jq37N8+T8RFLQRU6Z8jAULVnAj\nEhExgSEiIqpbiYkxCAnZjoMHf0Z29gvx82lp1Z/ynJ3NkZ+fBzMzc9y+ncONSUSkZnS5CYiIqLFZ\ntGg2kpNvQ1e35qe5adNmw9TUDFOnfsINSUSkhngFhoiIGq0LF8IQGNhP/L8yV2CIiEi98QoMERFp\nnLt3b8LBofoBq15eftxYRERMYIiIiBrW8eMHlSpnYmLGjUVExASGiIio4ZSXlyM0dI9SZfX19bnB\niIiYwBARETWc4OBtuHPnhlJltem+CERETGCIiIjUTELCVXzzzefcEERETGCIiIjU28mThxEY2A85\nOdncGERETGCIiEidCIVC7N+/C1OnjoS/fxs4OhrDxcUC/foJsGzZl3j06IFSn1NeXo4zZ05g9uxJ\n6N7dBe3amcDZ2RzdurXH9OljEBq6B0JhcY2WLT7+Cr744gP06tUJLi4WcHQ0hr9/GwQG9sOePVtQ\nWFhQ7WdcuxaLuXOnoGvXdnB0NIaHhw3Gjx+IoKBNSEt7iLCw36RmGRs7tj+mTh2Jly+zpD7HwUFH\n6i88/ORrbfeiokIcPhyCwMB+aNlSt9pyY8f2lykXHn4SU6eOhI9PSzg7m6N/fx/s2vUjymtw14Kc\nnGxs3boW48cPhI9PS7Rta4SOHa0wYIBvjX5/IqLGiPeBISJSMxcvhuPzz99H+/aumD17ATw8fFBU\nVIjTp4/hq68+RGFhAUxNzbBq1RaMGjVR4eckJcXjs8+mobS0FDNnzkOPHm+haVNrpKWlIizsODZu\n/A4ZGelwdGyPFSs2o1evAdUu28qVi7Bhwwp07OiJ5cs3wtPTFxkZ6QgN3Y3165ejpKQE7u4CHDx4\nFpaWTeR+xtata/HDD//Gp58uxODBo2FnZ4+srOe4fDkC69YtR1JSvLisvPu2SCY2ytzXRZnyiYkx\nCAnZgV9++VkqSapcvqpyT59m4IsvPsCpU7/K/Y6JE6fj++9/qnZ59+3biaVL56Fbt16YNWs+3N29\nkZWViaNHD2D16iXIzn4JAwNDfP31Krz//qfcYYiICQwRETWckJDtmD9/JsaMmYTVq7fJDEJfs+Yb\n/Oc//6o4gOvoICjoMAYMGCY3CZo8eRj69RuC9et3wcDAUKbMy5dZeO+9IYiO/hN6enr4z3+2YuzY\nKQqXbceODVi0aDYA4PLlZLRp007q9e++W4j165cDAGbOnIfFi/8j8xl//HEckyYNQXDwSfTuPVDm\ndaGwGJMnDxdfSamvBGbIkAAYG5sgLi4aeXm5CssPG9YDhoaGMuXOn7+B8eMHwte3K959dwa8vHyh\nq6uH48dDMX/+hyguLgIAhIaGo3v33gqXdfnyf2Ljxu8wb94SzJv3L5nXb926jsGDu4q/e9Gilfjo\noy+54xARExgiIqp/586dxsSJg9CyZRtERFyHkZGxTJm4uGgMGtRZ/L+rqzvCwxOlyjx8mIL+/QWw\nsGiC8PBEmJtbKPzO58+foW9fb2RkpEFXVxf79v2BHj36yC3bpYsjUlPvw9DQCCkphTKvJyUloG9f\nLwBAu3YuuHjxlkyZESN6IjLyAq5fz4SVlbXc78nMfIoePVyQnf2y3hIYkbNnf8eECW9XW/7cudMY\nN+7VFav27V3x7bcb8Oab/WXKrlr1NdauXQYAePvtEdix45Dcz9y/fxfmzJmM7t17IzQ0XOEybtq0\nEt9++xUAQE9PD+HhiXB27sgdiIi0BsfAEBGpAaGwGF988QFKS0sxZcrHcpMXAPD09EX79q7i/1NT\n78uUWbjwY2Rnv8TYsZOrTF4AwNq6Ob77bjMAoKysDHPmTFY4Jub582cAFN9bpW1bJ/HjJ0/S5ZZJ\nSLgKANi7V3FXqmbNbPDuuzMa5Hfw9vavVbmgoMNykxdR0iKSmBgjt0x+fp54drVp02ZX+d19+w4W\nPy4tLcX27eu5AxERExgiIqpfhw6F4OHDFABQGAgDgK6uLvbu/Q12dg6wtGyCzz9fKvV6fPwV/PHH\ncQBAv35DlPruAQOGiZOiR48e4JdffpZbbvr0OTAzM8fMmfPkvm5q+uqu94oG8puZmQOo6G62YsUC\nheXeeWdSg/wO1SV8r9bVXOr/qq6AtGvnIn789GmG3DIHD+4VJ4hdurxR5XdX7rp3/nwYdyAiYgJD\nRET16/Tpo+LHrVs7Vlm2bVsnxMQ8wo0bLzBjxmdSrx04sFv82MGhtVLfraOjgxEjxon//+03+V2c\n5s9fhtu3c2SSpszMp/jppx/Qr59A/Fxpaancz+jatSeAiqs9GzasQLdu7bF58/fIzn4pVc7V1V2p\n7mGqpq9voGQ5faU/U5S0ARCPhaksIuKU+LGnp63M7GqSf87O0slTenoqdyAiYgJDRET1Kz7+ivix\nou5jyvjzz1djJ5o2bab0+zp37iF3WaoSHf0nZs4ch549O+LOnRtYs2Z7te9ZsGAFLC2txP9nZKRj\n2bIv4e/fGosXz8H9+8kN+jtUnjThdcsBFVfNRBRNpSzZtezevQKkpZUr/Xf3bh53ICJiAkNERPXr\n2bMn4scvXjyv9eekpNwVPy4qKlT6fU5OHcSPRV2ZFLl0KQLDh7+Bd98dDC8vX1y+nIxVq7bAy8uv\n2u9p184Fx49fhq9vV6nnc3NzsG3bOrzxRgfMmjVR6+5zIrnNX+f3JyJiAkNERPVCsmX++vW4Wn+O\nZNetzMynSr9PckYwExMTuWUKCvIxd+5UjB7dG6amZoiIuI6PPvpS4f1eFGnf3hVHj17C5s3BcHV1\nl1n+Q4eC0aePh8L7qTRGRUWvupZduxbHHYKIiAkMEZF6k+zuFRb2m9LvS0yMwZAhAeL/mze3lXpN\nWZLd1iRnOXsVYBdi3LgB2LcvCAEBb2L37qOws7Ov9fpWjLsZjzNnErBz5xG5V2Tef380rly5rCW/\n/6sEUnQPHCIiYgJDpDKFhQUICtqEwMC+8PCwQdu2hvD0tMU777yFbdvWoaAgnxuJaqRjRw/x49DQ\nPTKD2hXZs2eLVPLj5/cqmbl48YzS35+XlyN+/MYbfWVe37BhBaKiLgIAvvrqW7k3xqxtIjNgwDAc\nO3YZO3cekZphq6SkBCtXLtKK319yFrN9+4LYjYyIsQIxgSFSndjYKPTq1QkLF36MCxfO4PnzZxAK\nhcjMfIqLF8OxePEcdO/ujL/+Os+NRUqTnDr5xYvnWLlyYbXvSU9Pxf79u/DWW/8QPzdkyDvix0eO\n7FM461VlosHzOjo6cqcwDg3dI37s6elb6/V0cNBBcvItua8NGDAMp0/Hwt+/u/i56Og/te73z83N\nwWefTVM44F/SqVO/YvTo3tyBiBgrMIEhIvmuX49DYGBf8f06FMnISEdgYF9cuHCGG42UUvmmk0FB\nm7BjxwaF5UtLS/HFFx9AV1cXo0dPFD//j3+MEncBy8rKxPbtG5T6flGiMGTIGLn3NJGcqjcjI12p\nz1Q0lfLBg3sVvsfCwhIrV/4o/t/YWHY8juSsXkKhsFH8/uPGTZVa15MnD+OTT96rciKGpKR4zJ07\nVWYcERExVmACQ0TiYOyTTyYhNzcHfn4B2LTpZ1y58hD37xcjLu4xNm7cKxX4CYVCzJgRyK4gpBQr\nK2uZ+6ssWjQbU6aMwIULYXj5Mgv5+Xm4fTsJu3f/FwMG+OLMmROYPXuB1LTEenp6WLs2CHp6egCA\ntWv/jbt3b1Zbt3/+eSusrZtjyZI1cstYW9uIH+/du0Xm9ZycbHzyyXtSzym6SeW2beuQmnpf4fK0\naeMkftyjRx+Z183NLcWPHzyoetrlygmAoqSqpuWFwmKp/8vKyhR+ZuUrKfKuitnatsDs2Quknvvl\nl5/x5ptu2LlzMx48uAehsBjPnz9DVNRFLF48B4MHB0BfXx9z5y7mDkTEWIEJDBHJ+v33I0hKises\nWfNx9OgljBw5Afb2rWBgYAAbGzuMGjURp05dRf/+r+5+npWVia1b13LjkVI++GAuxo6dIlPvAgP7\nwc3NGs7O5ujVqxO++upDJCXFo3fvgfjooy9lPsfPrxvWrdsFPT095ORk4733hiAl5Y7C4Hrp0nl4\n9OgBgoIOw96+ldxykvV6y5Y1+OGHf+Pp0ww8efIYQUEb0bOnK8rKymBoaCQud+lSBJKTb8kkNhXL\nNFjhDRj/+OMYAMDExBRffvlvmdclxwutXLkIL148R0ZGGr755nNMmzZKqmzlRCkjI63K36ByeUXL\nKDldNVD1jG9ZWZlS/6elPZRb7tNPF2LkyAlSzz18mIIFC2YhIMAJbdsawcPDBsOHv4Ft29bBwMAA\nO3Ycho2NHXceIsYKWkUnLa28nJuBSJng8h00b26L5cs3VVkuPz8Pfft6iccUuLl5ISyM06KScsrL\ny7F27TKsXbtMppVffODW0cGkSTOxZMmaKm96GRb2G7766kM8evQA5uYWmDHjMwwdGog2bdqhuLgI\nV65cxo8/fo+CgnysWbMDLi5uCj8rIyMdgwZ1xuPHj2Rea9bMBsuXb8LQoe/grbc8ceNGovi1Jk2a\nYteuX9GlyxsAKsbAiFhaWmHGjM8wZMgYtG3rhKys5zhx4hBWrPgnLCyaYMuW/VKTEkgGCFOmjJB5\n3s3NEyEhp2Br2wJCoRDp6an4+uvZOH36mLjMoEEj8fXX36Nly9ZSExEIhUI8fvwIixbNxunTR8XP\nDxw4HEuWrIGDQysYGBiipKQE9+/fxeLFjbjXwAAABtBJREFUc6RmCxs1aiIWLVoJW1t7cRe30tJS\npKenYvnyf+Lw4RCpZfjXv1bD3r6lzGQIot9/3bpvqxy/5OzcEf/97//QqZM3dxoixgpMYIhIvj59\nPHDo0Dmp+2UosmfPFsyfPxMAYGpqhjt3crkBqUbu30/G7t0/4syZE0hNvY/S0lI4OLRGz559MWnS\nTLi5eSn1Ofn5eTh8OAQnThzCjRuJePo0A3p6enBwaA1//24YOXKC1ADyqqSnp+K77xYiLOw35Obm\noFWrthg0aCRmzfpSvF+cOXMC8+a9j6KiQgwYMBTz5y+Tuqrj4KCDAwfOwN6+JcLDT+LcudNISkrA\nkyePYWhoCFdXDwwd+g4mTpwOMzNzhcty7FgofvjhGyQn30arVm0xcuQEzJjxmfg9komSImlp5VLL\npUz5mnxuTZdB0qNHD7BnzxaEh5/Egwf3kJ+fC2vr5vD09MXgwWMwatQElc0ER0SMFZjAEBEePLiH\ngICKfvxWVta4fj2TG4WIiIgYK6gAx8AQ1QF7+5bix61ateUGISIiIsYKKqLPTaA9qurOIK8bw/Tp\nY3D8+EGly7/O99dWbZajPrx8+UL8uGfPvqx8RETEGIMxBmMFFeEVGC1y40YWNm7cK9MvU9Ggsc2b\ng7Fu3S7x/y1btsH69buRkPCEG7MakjM+jRkziRuEiIgYYzDGYKygIrwCo0UsLa0watREODl1wIgR\nPcUz3EREnJY7INjAwBBDhozBnDmT0b17b2zdGqrUoDRF1PVqSV3444/jAIA+fd6Gm5snKx8RETHG\nYIzBWEFF9ObNW7KEm0G7tGjRElZW1ggL++3vg8sp6OsbICDgTZmyP/zwDZo3t8WOHYeqnBGIXikp\nKcGXX85Abm4OfvrpAGxsWnCjEBERYwzGGIwVVIRdyLTU5MkfYdiwseL/V636GidOHJIq8+uv+3D/\nfjJ+/PF/nK6zBg4e3IsHD+5h6tRPeI8GIiJijMEYg7GCinEaZS1WWFiAUaN6ITY2CgBgbGyCAwfC\n4OfXDWfP/o6QkO3YtCkY+vrsaaisnJxs9OnjASMjY5w+HQNTUzNuFCIiYozBGIOxAhMYUpXKd9e2\nsrLG4sXf4+TJI/jpp/0wNDTiRqqB+fNn4n//C8LRo3/Cy8uPG4SIiBhjMMZgrMAEhlQtIeEqRozo\niYKCfACAkZExoqMfoFkzG5V+j7pMcajsctT0s0+ePIxp00Zh1aotmDhxOisWERExxtCyGIOxQv3g\nGBiCp6cv1q3bBR2dip2/qKgQCxd+jHLmtkq7efMaPv30/zB58iwekIiIiBhjMFaoQ7wCQwCAhw9T\nMHZsf6k5yefNW4J58/7V6NZV1VdgMjLSMXRoN/j6BmDz5mDo6rJdgIiISBtjDMYK9YNbj5CRkY7P\nPpuKAwfCMHbsZPHza9YsxbFjoY1ufdPSypX6U8aLF88xYcJAuLi4YcOGPTwgERERaXGMwVihfvAK\njJbLysrEjBmBWLnyv2jXzgXFxUUYPboPrly5BAAwMTHF4cPn4enpy41VSXb2S4wb1x/GxiYIDj4J\nY2MTbhQiIiLGGIwV6hhTQC2Wk5ONDz8cj6VLf0C7di4AAENDI+zYcQj29q0AAAUF+Zg8eTiePHnM\nDVbpgDR+/ADo6Ohg9+5j1R6QysvLkZJylxuOiIgYY2hJjMFYgQkMqZhQWIzp08fgiy++gZubl9Rr\nNjZ22LnziHhHS09PxZQpw1FcXMQNh4oWpXfeeQulpaUICfkd5uYW1b5n5cpFmDt3CjceERExxtCC\nGIOxAhMYqgPr1n2LTp284OcXIPd1T09fbNmyX3x33JiYSKxYsUDrt9vjx48wcuSb0NfXx759p2Fp\naSVTpry8HMXFRXj6NANnz/6OyZOHY/365QgM/D9WPCIiYozRyGMMxgp1j2NgtFBcXDSGDeuOyMgU\n2Nk5VFn24MG9+OST9yoqi44Otm07iEGDRmrldktOvoVx4wYgNfV+jd9rbGyCuLjHsLCwZAUkIiLG\nGI00xmCsUD94BUbLJCRcxfvvj4ZQKERw8Ha8fJklt1xZWRmePHmM3NwcqdaCWbMmYteuH5GZ+VSr\nttuVK5cxbFiPWh2QAGDQoJE8IBEREWOMRhxjMFaoP7wCo2Xk3QNF3pTBCxd+jKCgTVV+Vl3coVZd\ntW9vJr6LcG0EB59E794DWQGJiIgxRiONMRgrMIEhIiIiIiKSwS5kRERERETEBIaIiIiIiIgJDBER\nERERMYEhIiIiIiJiAkNERERERMQEhoiIiIiImMAQERERERExgSEiIiIiImICQ0REREREWkJflykM\nERERERFpiP8HSKuc9RxV1VwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image('lh_schema.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Schematics**: A lighthouse at a position $(x,y)$ emitted $n$ flashes observed at $D_k$ on coast (indicated by star symbols). Our prior imposes the lighthouse to be within a rectangular area aligned with the coast. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Work out an analytic expression for the posterior probability distribution of the position of the lighthouse, $(x,y)$. State your assumptions at each of the following guided steps:\n",
    "    1. Given the geometry of the problem, write down the probability on the azimuth angle, $\\theta_k$, of emission of a $k$th-flash from the lighthouse that can be visible on the coast, $p(\\theta_k|x, y)$ .\n",
    "\n",
    "    2. Work out the relation of $\\theta_k$ with $D_k$ through simple geometric considerations.\n",
    "\n",
    "    3. Deduce an analytic expression for the probability of seeing a flash at $D_k$ given the (unknown) position of the lighthouse (i.e. the likelihood).\n",
    "\n",
    "    4. Finally, write down the posterior probability distribution of the position $(x,y)$ of the lighthouse given the ensemble of the measurements $\\{D_k\\}$. Give an expression for your normalization constant (but you don't need to do the integration).\n",
    "    \n",
    "2. Inferring the position of the lighthouse from the data involves the estimation of both $x$ and $y$. The full procedure is beyond the scope of this chapter and methods will be learned in later chapters. We will therefore assume that the position along the coast is known, $x=1.25$ and reduce it to a single parameter example. State your assumptions at each of the following guided steps to find the distance between the coast and the lighthouse:\n",
    "\n",
    "    1. Write down the posterior distribution of $p(y | \\{D_k\\}, x)$, and differentiate it with respect to $y$. It leads to a condition that is not easily solvable analytically. However, there is nothing to stop us tackling the problem numerically. \n",
    "\n",
    "    2. The most straightforward method is to use brute force and ignorance: gridding the values. Grid the $y$ 1d-space allowed by your prior with a sensible step size and compute at each point the posterior value $p(y |\\{D_k\\}, x)$. Be careful of the influence from your stepsize, make a sensible choice. The normalization of the posterior will be computed numerically on that grid. Plot the posterior. Locate and report the maximum of this function $y_{max}$.\n",
    "    \n",
    "    3. Numerically compute the expectation value $E[y]$, and the standard deviation $\\sigma_y$ of the posterior. On top of your posterior distribution of $y$, plot a Gaussian of mean $E[y]$ and standard deviation $\\sigma_y$. Does this agrees with $y_{max}$, given the width of the distribution?\n",
    "    \n",
    "    4. Make a quadratic approximation of the posterior distribution $p(y |\\{D_k\\}, x)$, around the peak value $(x, y_{max})$. Show your working and report the parameters of this approximation. Compare your result with the ones from the previous question (2.3). Why might they differ?"
   ]
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}