S46_logistics_learning.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github",
    "tags": [
     "no-tex"
    ]
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/gtbook/robotics/blob/main/S46_logistics_learning.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "JoW4C_OkOMhe",
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "%pip install -q -U gtbook\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "10-snNDwOSuC",
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.random import default_rng\n",
    "rng = default_rng()\n",
    "\n",
    "import plotly.express as px\n",
    "try:\n",
    "    import google.colab\n",
    "except:\n",
    "    import plotly.io as pio\n",
    "    pio.renderers.default = \"png\"\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "nAvx4-UCNzt2"
   },
   "source": [
    "# Some System Identification\n",
    "\n",
    "> We can \"learn\" measurement and motion models.\n",
    "\n",
    "<img src=\"Figures4/S46-Warehouse_robots-09.jpg\" alt=\"Splash image with various evolutions of warehouse robots\" width=\"40%\" align=center style=\"vertical-align:middle;margin:10px 0px\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section, we will limit ourselves to learning linear conditional Gaussians from data.\n",
    "In Section 4.1, we discussed three different ways to represent densities over a two dimensional state\n",
    "$x\\in\\mathbb{R}^2$: Gaussian densities, finite element approximations, and samples.\n",
    "In Section 4.4 we presented three different Bayes filters stemming from each of these choices.\n",
    "However, in each case the motion model was the same: a conditional Gaussian where the mean was a *linear* function of the previous state. In this section we will investigate whether we can learn both the linear mapping as well as the noise parameters, from data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maximum Likelihood Parameter Estimation\n",
    "\n",
    "Let us first consider a one dimensional system (i.e., $x_k \\in \\mathbb{R}$ is a scalar),\n",
    "and a sensor model that is linear with additive Gaussian noise:\n",
    "\n",
    "$$\n",
    "z_k = h(x_k) + n_k = C x_k + n_k\n",
    "$$\n",
    "\n",
    "Here, $n_k$ denotes zero-mean, Gaussian noise with covariance $R$. \n",
    "\n",
    "In the one dimensional case, both $C$ and $R=\\sigma^2$ are scalars, and estimating them will be similar to the Gaussian parameter estimation from Section 2.6. Remember that there we estimated a sample mean $\\hat{\\mu}$ from a set of $N$ data points, by simply computing the average of the data values:\n",
    "\n",
    "$$\\hat{\\mu} = \\frac{1}{N} \\sum_i z_i.$$\n",
    "\n",
    "The underlying reason is that this minimizes the negative log likelihood (nll) of $\\mu$ given the data $X$ by\n",
    "\n",
    "$$\\hat{\\mu} = \\arg \\min_\\mu \\sum_i \\text{nll}(\\mu; z_i) = \\arg \\min_\\mu \\sum_i \\frac{1}{2\\sigma^2} (\\mu - z_i)^2,$$\n",
    "\n",
    "hence this is a *maximum likelihood* (ML) estimator of $\\mu$. \n",
    "\n",
    "Suppose we are given a set of ground truth values $\\{x_k\\}$ for the state along with corresponding\n",
    "measurement values $\\{z_k\\}$.\n",
    "If we generalize the scheme from Chapter 2, the negative log likelihood of $C$ given the data for our measurement model becomes\n",
    "\n",
    "$$\\hat{C} = \\arg \\min_C \\sum_k \\text{nll}(C; x_k, z_k) = \\arg \\min_C \\sum_k \\frac{1}{2\\sigma^2} (C x_k - z_k)^2.$$\n",
    "\n",
    "Taking the derivative with respect to $C$ and setting to zero we have\n",
    "\n",
    "$$\\sum_k x_k(\\hat{C} x_k - z_k) = 0,$$\n",
    "\n",
    "yielding\n",
    "\n",
    "$$\\hat{C} = \\frac{\\sum_k x_k z_k}{\\sum_k x_k^2}.$$\n",
    "\n",
    "Fitting the variance then be done using\n",
    "\n",
    "$$\\widehat{\\sigma^2} = \\frac{1}{N-1} \\sum_k (\\hat{C} x_k - z_k)^2.$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ML Code Example\n",
    "Consider a one dimensional version of the warehouse problem in which the robot\n",
    "moves only along the $x$ axis.\n",
    "Let us simulate some states randomly from the $x$ direction in the warehouse, i.e., between 0 and 100 meters, and then the measurement $z$ from those sampled $x$ coordinates, and then plot them against each other:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArwAAAH0CAYAAADfWf7fAAAgAElEQVR4XuydC5xVZbn/f2vtPcwMigygMmgCExlgJFASghWIBZKFR88Jg4qbqZilBy0vlVrgJS3F9OCtkiMlipbmpRDzmqX4LxXLEvSYYBqMFxjuc9t7/T/v2rNn9uzZw957PWvtZ609v/U5J2X2+zzvWt/n3fLl5V3vazmO44AXCZAACZAACZAACZAACZQpAYvCW6aV5WORAAmQAAmQAAmQAAm4BCi8HAgkQAIkQAIkQAIkQAJlTYDCW9bl5cORAAmQAAmQAAmQAAlQeDkGSIAESIAESIAESIAEypoAhbesy8uHIwESIAESIAESIAESoPByDJAACZAACZAACZAACZQ1AQpvWZeXD0cCJEACJEACJEACJEDh5RggARIgARIgARIgARIoawIU3rIuLx+OBEiABEiABEiABEiAwssxQAIkQAIkQAIkQAIkUNYEKLxlXV4+HAmQAAmQAAmQAAmQAIWXY4AESIAESIAESIAESKCsCVB4y7q8fDgSIAESIAESIAESIAEKL8cACZAACZAACZAACZBAWROg8JZ1eflwJEACJEACJEACJEACFF6OARIgARIgARIgARIggbImQOEt6/Ly4UiABEiABEiABEiABCi8HAMkQAIkQAIkQAIkQAJlTYDCW9bl5cORAAmQAAmQAAmQAAlQeDkGSIAESIAESIAESIAEypoAhbesy8uHIwESIAESIAESIAESoPByDJAACZAACZAACZAACZQ1AQpvWZeXD0cCJEACJEACJEACJEDh5RggARIgARIgARIgARIoawIU3rIuLx+OBEiABEiABEiABEiAwssxQAIkQAIkQAIkQAIkUNYEKLxlXV4+HAmQAAmQAAmQAAmQAIWXY4AESIAESIAESIAESKCsCVB4y7q8fDgSIAESIAESIAESIAEKL8cACZAACZAACZAACZBAWROg8JZ1eflwJEACJEACJEACJEACFF6OARIgARIgARIgARIggbImQOEt6/Ly4UiABEiABEiABEiABCi8HAMkQAIkQAIkQAIkQAJlTYDCW9bl5cORAAmQAAmQAAmQAAlQeDkGSIAESIAESIAESIAEypoAhbesy8uHIwESIAESIAESIAESoPByDJAACZAACZAACZAACZQ1AQpvWZeXD0cCJEACJEACJEACJEDh5RggARIgARIgARIgARIoawIU3rIuLx+OBEiABEiABEiABEiAwssxQAIkQAIkQAIkQAIkUNYEKLxlXV4+HAmQAAmQAAmQAAmQAIWXY4AESIAESIAESIAESKCsCVB4y7q8fDgSIAESIAESIAESIAEKL8cACZAACZAACZAACZBAWROg8JZ1eflwJEACJEACJEACJEACFF6OARIgARIgARIgARIggbImQOEt6/Ly4UiABEiABEiABEiABCi8HAMkQAIkQAIkQAIkQAJlTYDCW9bl5cORAAmQAAmQAAmQAAlQeDkGSIAESIAESIAESIAEypoAhbesy8uHIwESIAESIAESIAESoPByDJAACZAACZAACZAACZQ1AQpvWZeXD0cCJEACJEACJEACJEDh5RggARIgARIgARIgARIoawIU3rIuLx+OBEiABEiABEiABEiAwssxQAIkQAIkQAIkQAIkUNYEKLxlXV4+HAmQAAmQAAmQAAmQAIVXOAb+/f5eYYbiw/v36YWqXjFs3dGExpZk8QkYISJg2PeujGHrzmZRHgZ7I2DY96qIoWEX+XsjKIvqUx0HLAs797TIEjHaE4G++1WgJeFgT2Orp3gGyQiY33/3NCXQ2JyQJfIQfciAag9RDEkToPAKxwKFVwgwguEUXt2iUXh1+VN4dflTeHX5U3h1+Ut6p/BK6AGg8AoBRjCcwqtbNAqvLn8Kry5/Cq8ufwqvLn9J7xReCT0Kr5BeNMMpvLp1o/Dq8qfw6vKn8Oryp/Dq8pf0TuGV0KPwCulFM5zCq1s3Cq8ufwqvLn8Kry5/Cq8uf0nvFF4JPQqvkF40wym8unWj8Oryp/Dq8qfw6vKn8Oryl/RO4ZXQo/AK6UUznMKrWzcKry5/Cq8ufwqvLn8Kry5/Se8UXgk9Cq+QXjTDKby6daPw6vKn8Oryp/Dq8qfw6vKX9E7hldCj8ArpRTOcwqtbNwqvLn8Kry5/Cq8ufwqvLn9J7xReCT0Kr5BeNMMpvLp1o/Dq8qfw6vKn8Oryp/Dq8pf0TuGV0KPwCulFM5zCq1s3Cq8ufwqvLn8Kry5/Cq8uf0nvFF4JPQqvkF40wym8unWj8Oryp/Dq8qfw6vKn8Oryl/RO4ZXQo/AK6UUznMKrWzcKry5/Cq8ufwqvLn8Kry5/Se8UXgk9Cq+QXjTDKby6daPw6vKn8Oryp/Dq8qfw6vKX9E7hldCj8ArpRTOcwqtbNwqvLn8Kry5/Cq8ufwqvLn9J7xReCT0Kr5BeNMMpvLp1o/Dq8qfw6vKn8Oryp/Dq8pf0TuGV0KPwCulFM5zCq1s3Cq8ufwqvLn8Kry5/Cq8uf0nvPU54WxMJ2JYN27a6cNu5aw/M5/369un02aNPP4/RRwzDQQNqusT8+/29Ev6eYs0XzkjX1h1NaGxJesrBIO8EKLze2fkRSeH1g6L3HBRe7+z8iKTw+kHRew4Kr3d22pE9Snj3NjbjlDO+j9O/8gV8/rMT2tnv2duICy67BY//6UX3Z0ceMQw3XHY2Duzf1/31uOkLcd3ib+CYcaMovNojNgT9U3h1i0Dh1eVP4dXlT+HV47+twUKyOY7+A5KwYomS38ghA6pL3mc5ddhjhPfHN6/C8rtWu7W76rtndBLen638Le558En84obvorqqF868cCnqBg/CkvMXUHjLabT79CwUXp9AekxD4fUIzqcwCq9PID2mofB6BCcI29sI3LUqhjc2pf5muKoamP7ZBMaOcQRZiw+l8BbPLDOixwhvw/ZdaGxuxuyvL8G5p8/sJLz/ddqlmDZ5HE778uddNmue/H849/s34uUnlsOyrE4zvO9v24GLrrgVE8eNwryZx4NLGmQDMIrRFF7dqlF4dflTeHX5U3hLz3/1GhvPPmd36riqCvjO+a0lvRkKrwx3jxHeNKZps76Nby44uZPwmiULl11wqiu95vrHqxvxxdO/j2ceXIa+ffZrF95RI+ow75wr3dnfqy9eiHgsRuGVjb9IRlN4dctG4dXlT+HV5U/hLS3/F9dZWPN7G3v2WHDaXv1JvwE0f04CdUNLN8tL4ZXVvscLr+M4GHXsfNx45SJMmjDapfn6xrcxY9538eiqazBo4ABXeC+/8GtYcc8a9O/XB9dcehYq4jG3bdIp3WBPl9q2Ul83jb5lw608otP/sSt95WX8ur6mKcvHaBIggfImELX/xvldjQcfTsL8f6fLQGn7j+kVl8RxYH+/e+0+X/r3/tL1WF499XjhNeVMC+3USUd1O8NrPjAvt62+4yoMPnRg+yjYrLRLQ2XbLg1N3KWh5N9Iw97MMm7b2VzyviUdlstvXpzhlYwCeSxneOUMJRlKOcPb0/6Q/Mp6s3Qh9dTbdwDvb7XQNr/UpWRDhzg4dW5pX1wbxJfWJF8dUHgBmDW8xx/7CXxt9gkuzFxreL8wdSI217+PTW9twcplF6Om7/5uW67hFY2/SAZzSYNu2Si8uvwpvLr8Sym8uk9a2t7/8oKFBx6MIT2Ba5YvWA7cZQzZ4n/MRAef/mQC1VWlvUcuaZDx7jHCa/bXdZIOPj/nIiycMwOf/8wEVFTEXXo/veMh/Oqhp9xdGnpXV2LhBdfm3KVh7KjDcep5V7sxP7/mfLcthVc2AKMYTeHVrRqFV5c/hVeXP4U3GP6X/zCOpqaO5QrpXtJ/M5aWXiO5P/ieg8bm0s7umvuh8Mpq32OE1+y6YGZuM6+HVqReQNu9pxHfWnwT/rD2JffjUcPrcMPl5+DgA1MHTZglD9cv+SYmHPURmN0eZp+1BIMPPRjLrliE+gbzDSntxYMnSss7uzcKry5/Cq8ufwqvLn8Kr//839hoYfmKjtndLj20T/sCXznFwpFHJim8/pch8Iw9RngLIbl95260tLS2HzhRSAxneAuhVF5tKLy69aTw6vKn8Oryp/D6y98cJtHQAFd4u7vMu+lmhnfwYQ7OPyeGPU0JCq+/ZShJNgqvEDOFVwgwguEUXt2iUXh1+VN4dflTeP3h/+QfbDz1RxuJVsDswxC3gEQSOV9S27+Pg+EfdjD1uCQOPagXhdefEpQ8C4VXiJzCKwQYwXAKr27RKLy6/Cm8uvwpvHL+t90ew8a2U9NMNnedrgOY3UaN9Kavqkpg+rTOJ6qZJYWc4ZXXQCMDhVdIncIrBBjBcAqvbtEovLr8Kby6/Cm83vmbQyRefdXG39d33XDNSG/dEAcL8mw1RuH1zl87ksIrrACFVwgwguEUXt2iUXh1+VN4dflTeAvjb+S2YXtKbMeMdvDiSxaefMpOTefm2GCYwlsY1yi3ovAKq0fhFQKMYDiFV7doFF5d/hReXf4U3vz8770/hnUvdXNsxj6E9+QZnZcv5OqJM7z5+Ye1BYVXWBkKrxBgBMMpvLpFo/Dq8qfw6vKn8O6b/95G4MqrU3vsd3e5uy5k+LBtAcdOTmLSp7KOEc6RgMKrO/4lvVN4JfR40pqQXjTDKby6daPw6vKn8Oryp/Dum396T93sVuaFtNbMsyIcoG+NgymTkhgxwin41DQKr+74l/RO4ZXQo/AK6UUznMKrWzcKry5/Cq8ufwpvbv5GdLfUW+jXz8HKu7ruqVs3BNi8BWhsOyvK7MAw65QE6oamz1IrrK4U3sI4hbEVhVdYFS5pEAKMYDiFV7doFF5d/hReXf4U3q78V66ysX6D3b69WGWlg+bmjjULRm7PPCM1vWsOmTBXbW3hs7qZPVJ4dce/pHcKr4QeZ3iF9KIZTuHVrRuFV5c/hVeXP4W3M38za3vTrV3X7B54kIOjxjroVwMMHZoseMlCvupSePMRCu/nFF5hbTjDKwQYwXAKr27RKLy6/Cm8uvwpvMDjT9mpLcYAVPQCmptyn5A2a2YSI0fkfxGtmIpSeIuhFa62FF5hPSi8QoARDKfw6haNwqvLn8Kry7+nC293L6Xl2m1s8qSk+1KanxeF10+apc1F4RXypvAKAUYwnMKrWzQKry5/Cq8u/54uvL9dE8PatW3rc62OMySytxozVTp+ahITj6bw6o7Y8PRO4RXWgsIrBBjBcAqvbtEovLr8Kby6/Huy8D6z1sbDj6SWMpgrvb+C0V/bBpIZbpt+Ua1fTXG7MOSrLmd48xEK7+cUXmFtKLxCgBEMp/DqFo3Cq8ufwqvLvycL7yWLuz9QYsG8JF540XJ3YaipASaMT2BQrf+1ovD6z7RUGSm8QtIUXiHACIZTeHWLRuHV5U/h1eXfU4W3u7W7lZXAbA/76XqtIoXXKzn9OAqvsAYUXiHACIZTeHWLRuHV5U/h1eXfU4W3uyODRwxPYvYp/q7T3VeFKby641/SO4VXQo/78ArpRTOcwqtbNwqvLn8Kry7/chZeM4u7cZOFxkZgUK2DMaM7r79NHzCRWYH5c4o/LU1SQQqvhJ5uLIVXyJ8zvEKAEQyn8OoWjcKry5/Cq8u/XIX3lfU27ry744U0Q9kI78knpk5IS1/mxbX08cB1Q5yijwaWVo/CKyWoF0/hFbKn8AoBRjCcwqtbNAqvLn8Kry7/chXe226PubO72dfiS1p1gWf1TuENVTmKuhkKb1G4ujam8AoBRjCcwqtbNAqvLn8Kry7/chJes4QhfWLa5i1W+8xtJuGLzm/17VhgPypH4fWDok4OCq+QO4VXCDCC4RRe3aJReHX5U3h1+ZeL8G5rsLD0+lg7zFwnpZm9dL9zAWd405AOGVCtO/gi3juFV1hACq8QYATDKby6RaPw6vKn8Oryj7LwZs7obt9hYeu2Dpau8GZZbxAnpUmrxxleKUG9eAqvkD2FVwgwguEUXt2iUXh1+VN4dflHVXgLmdE1vjv2SAdjxyTdwyP8PiXNj8pReP2gqJODwivkTuEVAoxgOIVXt2gUXl3+FF5d/lESXrN37rqXbOzcBazfYOO9d9vYWaljgbu+ogaceXprICek+VU1Cq9fJEufh8IrZE7hFQKMYDiFV7doFF5d/hReXf5REt5lt8RRXw84DmBl2a35mTHezB9PnpTElEmlO0TCSyUpvF6ohSOGwiusA4VXCDCC4RRe3aJReHX5U3h1+YddeF9cZ7mzutt3Wtj6fseMbja1tASHfUY3+74pvLrjX9I7hVdCjyetCelFM5zCq1s3Cq8ufwqvLv8wC695KW35io6dFwyp9Flp2csXjPAeOzn8M7oUXt3x7mfvFF4hTc7wCgFGMJzCq1s0Cq8ufwqvLv8wC+/jT9nt++p2omRWKXQ+RC3nKWq6ZAvrnTO8hXEKYysKr7AqFF4hwAiGU3h1i0bh1eVP4dXlH2bhXb3GxrPPZZmtwWWE10zxtk3zDhwIzD4lEcpdGPJVl8Kbj1B4P6fwCmtD4RUCjGA4hVe3aBReXf4UXl3+YRZes373vgc6L2kwtMaPd3DE8CSqqhxUVVmRFN101Sm8uuNf0juFV0KPa3iF9KIZTuHVrRuFV5c/hVeXv7bwvrLexpZ6oKoKGDHc6SKvmbO85qS06dMSGDsmvZJXl50fvVN4/aCok4PCK+TOGV4hwAiGU3h1i0bh1eVP4dXlrym82UsWjPTOnxPufXP9rhaF12+ipctH4RWypvAKAUYwnMKrWzQKry5/Cq8ufy3h3Vxv4aZbui5XiMWAGSeU1yzuvipM4dUd/5LeKbwSelzSIKQXzXAKr27dKLy6/Cm8uvw1hNcsY7hzlZ3zaLT0frqLzo7mS2jFVpPCWyyx8LSn8AprwRleIcAIhlN4dYtG4dXlT+HV5a8hvObENLNuN9dRwGkaJ83oGbO8FF7d8S/pncIroccZXiG9aIZTeHXrRuHV5U/h1eWvIbyXLI7nfmgHcNqOB6bwBj8uDhlQHXwnZdwDhVdYXM7wCgFGMJzCq1s0Cq8ufwqvLv9SCK9ZwrD2OQt7myxUVzpo2G5hW0Pquc0SBneqt+2f6VlfLmkIflxQeGWMKbwyfqDwCgFGMJzCq1s0Cq8ufwqvLv8ghdeI7kMPW9i5o/PihXgF0NKSsaTBAfrWANu3A0OHOJgw3sHIEeZ0ifK/uKQhujWm8AprR+EVAoxgOIVXt2gUXl3+FF5d/kEJ77YGCzfeEkNjU+61up+cmIQRX3PVDXFQN7R89tYtpqIU3mJohasthVdYDwqvEGAEwym8ukWj8Oryp/Dq8g9CeB9/ysaTT6Z2YUivWMh+yp6yRjdfdSm8+QiF93MKr7A2FF4hwAiGU3h1i0bh1eVP4dXl77fwbt4C3HRrvMN0uzHeM0/vWQdMdFdlCq/u+Jf0TuGV0OMuDUJ60Qyn8OrWjcKry5/Cq8vfb+F1Z3efslMzu22ym+m8FRXAcccmMfHonrFGN191Kbz5CIX3cwqvsDac4RUCjGA4hVe3aBReXf4UXl3+fgvvi+ss3PdA2wlqZpsxI74W0Gd/4OMfT2LKJIpuZsUpvLrjX9I7hVdCjzO8QnrRDKfw6taNwqvLn8Kry18ivObFtGefs1C/xcK27cCggQ4+9jEHv7439bJa5tVTthkrtpoU3mKJhac9hVdYC87wCgFGMJzCq1s0Cq8ufwqvLv98wru3Ediwweyba6FfjYMxox2Ydbp33h3Dtm2WO3ubedXUODhpRtIV4cbG1IfHTkr22F0Y8lWXwpuPUHg/p/AKa0PhFQKMYDiFV7doFF5d/hReXf77El4ju0uvj6OxseMeDzwIqIglsXmLjaQD2DnOB+YLaYXXlMJbOKuwtaTwCitC4RUCjGA4hVe3aBReXf4UXl3++xLeZ9baePgRu+MGc+y4YE5Ky57lpfAWXlMKb+GswtaSwiusCIVXCDCC4RRe3aJReHX5U3h1+XcnvGZ21yxb2PiGBSc9i5tDbrMduKoS+M4FrboPFaHeKbwRKlbWrVJ4hbWj8AoBRjCcwqtbNAqvLn8Kry7/7oT32p/E0LA9a71Cd6dItD1CTV8HJ53I9brFVJTCWwytcLWl8ArrQeEVAoxgOIVXt2gUXl3+FF5d/rmE9xe/jOG1f3ZdnJtr+YK5+8mTuN2Y1ypSeL2S04+j8AprQOEVAoxgOIVXt2gUXl3+FF5d/tnCe/evY3j5Zcs9FrjLZbbQzVjSS9mV147CK2eolYHCm0F+z95GNDe3oqbv/p3q8ejTz2P0EcNw0ICaLnWi8GoNXb1+Kbx67E3PFF5d/hReXf7ZwvuDy+NoTbT5rjk4IsN9zQyvEWHjwh8a5uBTx3D5grR6FF4pQb14Ci+A+ne34bLrVmDtC6+4lRg+7DB895yvYOThQ9xfj5u+ENct/gaOGTeKwqs3VkPTM4VXtxQUXl3+FF49/mY/3TfeqMA/Nzro3Tvpiu7fXrY7Te6ml+3WDgTmz029jFZdpXfP5dYzhTe6FaXwAjh/yc1o2LELy674b1i2hR9cczvefX8bbr7qPApvdMd2YHdO4Q0MbUGJKbwFYQqsEYU3MLTdJjbbjf3pWQs7d1rIXpdrBNdcmSsaTJsll3LnhSAqReENgmppclJ4AXzlG5djyAcG4vILv+ZSv2/107jhtnvx+D1Luwjv+9t24KIrbsXEcaMwb+bx4JKG0gzUMPVC4dWtBoVXlz+Ft7T839hoYfmKWEenufbWzRLeoz7uYMYJidLeaA/pjcIb3UJTeAE8/scX8M3vXY/jPvUxnDT9U/jRjXdhwZc+h//6/KROwjtqRB3mnXMl6gYPwtUXL0Q8FqPwRnfse75zCq9ndL4EUnh9weg5CYXXMzpPgb9bY2Ptc6k3z7rbZSxz1tdsNXbuOZRdT7ALCKLwFgAppE0ovADe3vIeTvvWj/DhDx6GP/35ZVRVVmD50gvxobpD24XXzP6uuGcN+vfrg2suPQsV8dSfuFsT5jXY0l4x23ZPykkkk+5fb/EqLQHLstzjORPmnE5eJSdg+JvxnyT/krM3Hdptx3Ql+R+fQPk/9pSDhx5JYs/e3BswZHaeFt4PHGJhziwLgw/NtWVDoLfbY5LHbMs9otlRGP/xWNaWGz2Guj8PSuEFcMoZP8CkiWPw9bknYueuPbj0x/+Lp5/7K559aJk7i2teWjOX2cVh9R1XYfChA9vpv9PQ5E8lishSs18FelXYaNjdguaW0gt3Ebdalk0rK2yYWd7tu1vK8vnC/lDVvWz3D5w79pC/Rq32q4q5Z9Pu3ss1okHx39pg4cdLrdSMbp7DI8w9GPcyL6Z9+9wkX1ALqihtec0uGY3NCTQp/N57cE1lwE9X3ul7vPDu3tOIT3xuIW647GxM+eTH3Gr/fcNGzDzj+/jN8stweN0HXOH9wtSJ2Fz/Pja9tQUrl13cvnUZ1/CW9xck19NxSYNuzbmkQZc/lzQEz/+OO2PY8FrHLG16MtFMrrv/bv7fTPaZo4MBVFUB06clMHYM/9Yp6OpwSUPQhIPL3+OF16CdNuvbqBtci6u+txC9qypx3U9/hSeeeREP3H5F+wyv2ZZs7KjDcep5V7vV+Pk156N3dSXX8AY3NkObmcKrWxoKry5/Cm+w/M3WYzfdGu/SSa5T046fmsTIEQ761VB0g61KR3YKb6lI+98PhRfAK69twk0r7sdjT7+A3tVVOGr0cHd5w0dHftAlbmZ4r1/yTUw46iNo2L4Ls89agsGHHoxlVyxCvcKSBvOFM9K1dUcTGhX+WsX/YRitjBRe3XpReHX5U3iD5f/4UzaefKrrWk2jtJ+bmoSTtJFwgA8cmkDdUIpusNXomp3CW2ri/vVH4c1gaZY3tLYm0PeA/QomzCUNBaMqm4YUXt1SUnh1+VN45fzNVmPmqqpyMKi2cz6z5+7Dj3QVXrNc4aQZCWSftCa/G2YohgCFtxha4WpL4RXWg8IrBBjBcAqvbtEovLr8Kbze+e9tBJaviGPLlo4ckyclMWVSx8vH2xos3HRLDI1Z70MvOjvhLl2g8Hrn70ckhdcPijo5KLxC7hReIcAIhlN4dYtG4dXlT+H1zr+72dszT2/tNNNrZoCffc5CY6PlzgJPGO+0L1+g8Hrn70ckhdcPijo5KLxC7hReIcAIhlN4dYtG4dXlT+H1zv/e+2NY91LXPXLnzyl8PS6F1zt/PyIpvH5Q1MlB4RVyp/AKAUYwnMKrWzQKry5/Cq83/mY5w33321i/oev6XAqvN6YaURReDer+9EnhFXKk8AoBRjCcwqtbNAqvLn8Kb/H8zVZjZu2ukd7s+d2BA4Gzzij8EA/O8BbP388ICq+fNEubi8Ir5E3hFQKMYDiFV7doFF5d/hTewvm7L6ndHsOW+pTmmk3EzOtp6Tne6VOTGDumuNPRKLyF8w+iJYU3CKqlyUnhFXKm8AoBRjCcwqtbNAqvLn8Kb+H877grhg2vdp7TzTwpuJilDOleKbyF8w+iJYU3CKqlyUnhFXKm8AoBRjCcwqtbNAqvLn8Kb+H8r/lJHNu352jfZr0Xnd+K6qrC85mWFN7iePndmsLrN9HS5aPwCllTeIUAIxhO4dUtGoVXlz+Ftyv/Da9ZePNNCwcd6GDIELQf9XvVNXHs3p1beCccncT0aR377xZaVQpvoaSCaUfhDYZrKbJSeIWUKbxCgBEMp/DqFo3Cq8ufwtvB36zRXXZzHDt2dK5JeqnC7XfE8PrrWUsaHOD4aUkcc3TxsssZXt2xb3qn8OrXwOsdUHi9kmuLo/AKAUYwnMKrWzQKry5/Cm8H/9VrbDy71u6y9cLQIQ4WzE1gWwNw2+0xNDSkpNeygPHjHZwwLeG5iJzh9YzOl0AKry8YVZJQeIXYKbxCgBEMp/DqFo3Cq8ufwpvi787u3hTHjp1d69GrF/C9Czu2GjPbkplT02pqOpY7eK0ihdcrOeztL/0AACAASURBVH/iKLz+cNTIQuEVUqfwCgFGMJzCq1s0Cq8ufwpvSnZvvjWGrdus1L66WZvrfuBQB6ef6n0Wd18VpvDqjn8Kry5/Se8UXgk9ABReIcAIhlN4dYtG4dXl3xOF1wju6jUxbNyUYm9bFrZuS/175jZj6crMmpnEyBHe1ujmqy6FNx+hYD+n8AbLN8jsFF4hXQqvEGAEwym8ukWj8Ory74nC+0uzn+4Gy53JdSdzsy3XARwLMEsZZp6cwPAPmwbBXBTeYLgWmpXCWyip8LWj8AprQuEVAoxgOIVXt2gUXl3+PUl4N9db+PV9Nt55p2PNQq4Z3XRFRgxPYvYpwczspvug8OqOfwqvLn9J7xReCT0uaRDSi2Y4hVe3bhReXf7lLrzmBbPHnrDxzzdstHa8d9YJei7prawEvn5Gon0P3qCqROENimxheSm8hXEKYysKr7AqnOEVAoxgOIVXt2gUXl3+5S68S66Io6VNdLubzXUcoKoSaGpO1eLgg4D/PKkVg2qDrw2FN3jG++qBwqvLX9I7hVdCjzO8QnrRDKfw6taNwqvLv5yF95HHYvjjn/IvX6g92MH8eYmijwX2o3IUXj8oes9B4fXOTjuSwiusAGd4hQAjGE7h1S0ahVeXf7kK7zNrbfzuERt2Bt5cM7w1fR2ce04wW44VUlkKbyGUgmtD4Q2ObdCZKbxCwhReIcAIhlN4dYtG4dXlX07C+9JfbfxujY29e1IbL5gre1/d9v0WHOCwDzj46pd1ZnbTVafw6o5/Cq8uf0nvFF4JPS5pENKLZjiFV7duFF5d/uUivBtes/DLO2MdktuGNekgNcvbtrLBPQ74E0nUDUFge+sWU1EKbzG0/G9L4fWfaakyUniFpDnDKwQYwXAKr27RKLy6/KMuvOlT0rY1tBlt2x66Vtv6BfNCmpFc888++wOnnRr8zgvFVJTCWwwt/9tSeP1nWqqMFF4haQqvEGAEwym8ukWj8Oryj7LwvrLedvfVbW7pzDC9VtddvtAmvOZfgzwxzWsVKbxeyfkTR+H1h6NGFgqvkDqFVwgwguEUXt2iUXh1+UdJeM0s7qZNQGWV5a7V3b6tba1u+sS0bJQOcECNg4+NcVA3xEHd0OBOTPNaRQqvV3L+xFF4/eGokYXCK6RO4RUCjGA4hVe3aBReXf5REd4X11m474HUGt32K71soe0HHRuQpX5g9tY9swSHR0gqSOGV0JPHUnjlDLUyUHiF5Cm8QoARDKfw6haNwqvLPyrCe/lVcTQ1ZfluxnKFbA8eOzqJk08M9lhgPypH4fWDovccFF7v7LQjKbzCClB4hQAjGE7h1S0ahVeXfxSE17yYduXV8X2Ccl9OAxCvAGbObMWID+lyLbR3Cm+hpIJpR+ENhmspslJ4hZQpvEKAEQyn8OoWjcKryz8KwmteTrtzld2+tViaWHoHBvPrWAz40hcTGDzYUTkxzWsVKbxeyfkTR+H1h6NGFgqvkDqFVwgwguEUXt2iUXh1+YddeDdvAX7/WAyvvm65M7id1ulmHJ325S8lMPzD4XspLV91Kbz5CAX7OYU3WL5BZqfwCulSeIUAIxhO4dUtGoVXl38YhdfsxtDQADz3Zxv/eKVDcY3OtjtuEhg5MonaWmDsaAf9aqInu6byFF7d8U/h1eUv6Z3CK6HHk9aE9KIZTuHVrRuFV5d/2ITXLF+47wEbjY25uaSFd+BA4KwzWnXh+dA7hdcHiIIUFF4BPOVQCq+wAJzhFQKMYDiFV7doFF5d/mET3mW3xFFfnzmV25lP7UAH06clQ7mnrpdKUni9UPMvhsLrH8tSZ6LwColTeIUAIxhO4dUtGoVXl3/YhPeSxfvejWHE8CRmnxL+7cYKrSqFt1BSwbSj8AbDtRRZKbxCyhReIcAIhlN4dYtG4dXlHzbhbZ/h7QbL/DmJspndNY9I4dUd/xReXf6S3im8EnpcwyukF81wCq9u3Si8uvzDJryPP2Xjyafsdihm67GhQ1PHAkf55bTuqkzh1R3/FF5d/pLeKbwSehReIb1ohlN4detG4dXlHzbhNTTMi2ubzTpeAIMGAiNHlM8ShuxqU3h1xz+FV5e/pHcKr4QehVdIL5rhFF7dulF4dfkHKbzmhLQtW1LbitXUILJbhwVZIQpvkHTz56bw5mcU1hYUXmFluIZXCDCC4RRe3aJReHX5ByW85sCI5SvinbYXmzUzWdaztV4qSeH1Qs2/GAqvfyxLnYnCKyRO4RUCjGA4hVe3aBReXf5BCe/KVTbWb+hYi2ue0uyhe9ihDr765USkjv8NskIU3iDp5s9N4c3PKKwtKLzCylB4hQAjGE7h1S0ahVeXf1DCe+1PYmjY3ukg4PYHnTA+6e6ly4u7NGiPAQqvdgW890/h9c7OjaTwCgFGMJzCq1s0Cq8u/yCE17x0tnqNjYbtuZ/NHB7x9TMSug8ekt45w6tbCAqvLn9J7xReCT0Kr5BeNMMpvLp1o/Dq8vdbeNP76JrtxKzcE7yg8HbUnMKrO/4pvLr8Jb1TeCX0KLxCetEMp/Dq1o3Cq8vfT+E1M7t33p2xbtcBHCO9ZvEuOgT4uClJTPoklzQYJhRe3fFP4dXlL+mdwiuhR+EV0otmOIVXt24UXl3+UuF9Y6OFjZtSU7lmG7K1z3V+Uc383PVdB7BtYPiHy+toYGn1KLxSgrJ4Cq+Mn2Y0hVdIn2t4hQAjGE7h1S0ahVeXvxfhNaehbdpoof5dC7t3m6lb9/9SZptjGcNJMxKorXUwqFb3WcPYO4VXtyoUXl3+kt4pvBJ6nOEV0otmOIVXt24UXl3+xQpv9tG/7t1niG5FBdDS0vFMAwcCC+a2chuybspM4dUd/xReXf6S3im8EnoUXiG9aIZTeHXrRuHV5V+o8G5rsNDQAHf3hS31Xadx0847dIiDCeMd92jgfn0djBjhUHb3UWIKr+74p/Dq8pf0TuGV0KPwCulFM5zCq1s3Cq8u/30Jr3kJ7YmnLFdw3V0X0pO5OZYtpIV3xHCu0S2mohTeYmj535bC6z/TUmWk8ApJcw2vEGAEwym8ukWj8Ory70541/8fsHJlvH21Qlpocy7TbfthVSUw65QE6oa2bcug+2iR6J3Cq1smCq8uf0nvFF4JPc7wCulFM5zCq1s3Cq8u/+6E97IfxtHc3HFvaYVNz/Sm99jts7+DkSOAQw9JYuhQoF8NZbeYilJ4i6Hlf1sKr/9MS5WRwiskzRleIcAIhlN4dYtG4dXln0t4//KChfsfjHXsvtB2i0Z2zWVk16zVXTCXp6VJq0fhlRKUxVN4Zfw0oym8WfRbWlrxzvsNOKh/X/TqVeF++ujTz2P0EcNw0ICaLrWi8GoOX52+Kbw63NO9Unh1+WcKr3kh7dm2fXTdVQrm4IiMAyPS++ka4Z01M4mRI3h4hLR6FF4pQVk8hVfGTzOawttG/403N+OSHy3HC3971f3JxYvm4EsnTnH/fdz0hbhu8TdwzLhRFF7N0RqSvim8uoWg8OryTwvvH55pxX0PxDrdTObiBPOemmUDR49LYsRwh+t0fSobhdcnkB7TUHg9ggtBGIUXQP272zDli4swfcp4zD7pOIw8fCgam5rQr28fCm8IBmnYboHCq1sRCq8u/7Tw3rI8gfUbup6SZqZ4jfge0Af42oJW9Ov6F2O6DxDx3im8ugWk8Oryl/RO4QVw9bI78eDvn8ETv74O8VjnGYvsGd73t+3ARVfcionjRmHezOPBJQ2S4RfNWAqvbt0ovHr8zVHAL7wQx+v/BLbvAN59L8e9OMCUY5OY/GkuXwiiUhTeIKgWnpPCWzirsLWk8AKYMfc7qK6qxKCBA7C5/n2MPHwIFs6dgdqD+nea4R01og7zzrkSdYMH4eqLF7pyTOEN25AO/n4ovMEz3lcPFF49/tf+JIaG7alNdd3dF7L313WAqipg0Tk8KS2oKlF4gyJbWF4Kb2GcwtiKwgvgI5PnYfzYkThp+qfQq1ccP73jt9iztxH3L78cFRVxdw3v5Rd+DSvuWYP+/frgmkvPQkU8NRPc1FL6t44r4jZsy0JLaxLJ9GvQYRxdZXpPhr1tW2hNcAZLo8Tkr0EdeO3/gOtu7jzmzdKF/fez0NTkoLISGHMkMO1YGwNScwW8AiAQty2YKiST3M4tALx5U8Zjtste4/feyoqufwOd94bZoJ0AhbdNeK9fcjaO+9THXDDmBbbPz7kI9/58CYYPO8wVXnMZCV59x1UYfOjAdoDvbW8q+XA6YL8K9Irb2LG7Bc2tlK5SF6BXhY2qChs79rSWumv2B6Cqlw3zm86uveQf5IAwyxce+p2Ff29JTeMO6Ofg5Ve6Hpl2xEgHc2ZRvoKsRWbu/arjSCQcNDaXfrKlVM8Y5n4O6B1HY0sSzS2l/733wL6VYUYT+nuj8AL4r9MuxQnHHY35X5ruFuz1jW9jxrzv4q6bL8VHR9S5wvuFqRPd5Q6b3tqClcsuRk3f/d22XNIQ+jHu+w1ySYPvSItKyCUNReHy3HjlKrvTS2k5T0wDcPzUJCYeXfrf/D0/WMQDuaRBt4Bc0qDLX9I7hRfAbXf9DsvvWu0K7v77VWPpLffgsT8+j0fuugbVVb3atyUbO+pwnHre1S7vn19zPnpXV1J4JaMvorEUXt3CUXhLw/+SxfEuHVVWOWhq7JjlNYdJmKOBq6tKc0/sBaDw6o4CCq8uf0nvFF4Azc0t+M4Pf4bVjz/nshx4UD9c94Nv4Mgjhrm/NjO81y/5JiYc9RE0bN+F2WctweBDD8ayKxahvqH0SxrMF85I19YdTe5frfAqLQEKb2l5Z/dG4S0N/1zCW9PXwTfOsLBtm4Wk1YJBtaW5F/bSQYDCqzsaKLy6/CW9U3gz6O3YtQe7d+9F7cH9YXV5/Tg3Zi5pkAy/aMZSeHXrRuGV83/uzxYef9LGnr0W3AMiLGDI4NThEGNGO+hX4+C222PYuKnzmt0Rw5M4Y57tBuzc0yK/EWYomgCFt2hkvgZQeH3FWdJkFF4hbgqvEGAEwym8ukWj8Mr4L7slji31cEU3+0qv050/J4GqKgf33h9HfX2qlVm+cNKJSQweFKPwykogiqbwivCJgym8YoRqCSi8QvQUXiHACIZTeHWLRuH1zv+NjRaWr4jl3kPX7K1rZnvb5HbB3Ny7AKRPWuMMr/c6SCIpvBJ68lgKr5yhVgYKr5A8hVcIMILhFF7dolF4vfN//CkbTz5lU3i9I1SPpPDqloDCq8tf0juFV0KP25IJ6UUznMKrWzcKb3H8zX66616ysXEj3ANr/r4+dXCBnSNNeobXrNWdfUruF2I5w1scf79bU3j9JlpcPgpvcbzC1JrCK6wGZ3iFACMYTuHVLRqFtzj+2S+fpY8ETh/SaF5Yy/z3qkq4W42ZF9hyXRTe4vj73ZrC6zfR4vJReIvjFabWFF5hNSi8QoARDKfw6haNwls4fzO7e8XVqf10zdrc9AzuwQcCR3wkiUEDgaFDk2hoAN7YaKO6ysHQoXB3aejuovAWzj+IlhTeIKgWnpPCWzirsLWk8AorQuEVAoxgOIVXt2gU3sL5/+peG399uWPxQlpj64Y46O6ltHzZKbz5CAX7OYU3WL75slN48xEK7+cUXmFtKLxCgBEMp/DqFo3CWxj/bQ0Wll4fy9nY7LV78om5d2HIl53Cm49QsJ9TeIPlmy87hTcfofB+TuEV1obCKwQYwXAKr27RKLyF8U9vQZbdOhYDzj4rsc9lC/vqgcJbGP+gWlF4gyJbWF4Kb2GcwtiKwiusCoVXCDCC4RRe3aJReLvyf/4FC03NFgbVpk5LM1d3M7wTxicxfZr3I8kpvLrjn8Kry5/Cq8tf0juFV0KP25IJ6UUznMKrW7eeLLzmJbTVa2LYUp86J63vAQ42vGa2WUjVxPxjwADgjFNbUV0F3Ht/DOte6jhTzezAMH9uKwbVeq8hhdc7Oz8iKbx+UPSeg8LrnZ12JIVXWAHO8AoBRjCcwqtbtJ4svNlbjOWqhJHe6VOTmHh0ahb3lfU2NtcDRnZHjnA8L2VI90Xh1R3/FF5d/hReXf6S3im8Enqc4RXSi2Y4hVe3bj1ZeC9ZnNpiLPNKbzWW+TPpsoV9VZjCqzv+Kby6/Cm8uvwlvVN4JfQovEJ60Qyn8OrWjcKbJbwOYA6PaL8cYPLkJKZM8r5Ol8KrO8b31TuFV7c2FF5d/pLeKbwSehReIb1ohlN4devWk4V32S1x1Nfn57/obO+7MOTLzhnefISC/ZzCGyzffNkpvPkIhfdzCq+wNlzDKwQYwXAKr27RerLwbt4C3Pa/cTQ1p2pgjgR2lzS0zfD2rnaw8LSkeJ0uZ3h1xzhneMPLn8Ib3trkuzMKbz5CeT6n8AoBRjCcwqtbtJ4mvGZ7MbPTwsaNFoYOdVBXB/x2td0+01vZC5gyOekeESzZfaHQqnKGt1BSwbTjDG8wXAvNSuEtlFT42lF4hTWh8AoBRjCcwqtbtJ4kvGYbsqXXx9HY2MG8qgpYdHar+4PGRivQ2dxclabw6o5/Cq8ufwqvLn9J7xReCT2u4RXSi2Y4hVe3buUovCtX2Vi/wXbBmpPQRg538KlPJvDGRhsPP5L6eeY1a2YSI0cE81JavupSePMRCvZzCm+wfPNlp/DmIxTezym8wtpwhlcIMILhFF7dopWb8D6ztqvUuutyHeCoo5L4y/NdhXfypOB2YchXXQpvPkLBfk7hDZZvvuwU3nyEwvs5hVdYGwqvEGAEwym8ukUrN+HNnN3NJGteSBs2zME//5m551iqxfw5ifYjhEtdDQpvqYl37o/Cq8ufwqvLX9I7hVdCj0sahPSiGU7h1a1buQivmdltagL+9ncL775noYvWOsDgwxzYMWDjpo5Px4x2cPKJCbUiUHjV0LsdU3h1+VN4dflLeg9ceP931cMYelgtPjn+o4ibxWlldnGGt8wKWsDjUHgLgBRgk6gKr9ltYdOmFJjn/mLj7bezFDfzyLS2f//EUQ4+/7kEzHZk5gW1qiqnJDsx7Kt8FN4AB3cBqSm8BUAKsAmFN0C4AacOXHh/cO3tuPuBJzDwoH6YO/N4/Me0T6LvAfsF/FilS0/hLR3rsPRE4dWtRBSF98V1Fu57IOsP/FlnApslDO7V5sFm5e7CMxIYNDD9gS73dO8UXt06UHh1+VN4dflLeg9ceM3N/e2Vf+Ku+x/Hbx7+o3uvM2cciy+dOAXDhx0mufdQxFJ4Q1GGkt4EhbekuLt0FkXhveKqOBqbOj+KEdzMI4H79gWSCQc7d1no1w84/rN6OzFwhld3jO+rdwqvbm0ovLr8Jb2XRHjTN7i1YSfuf/iP+MWvH0H9u9swbswIfPU/p2LSxNGRXe5A4ZUMv2jGUnh16xZ24TV75/7mfhuvtG0zVlsLbNmSn9mE8UlMn6az1Vj+u+towRneYmj535bC6z/TYjJSeIuhFa62JRXe7Tt244FH/oTlq1a7wtu7ugp79jaif00fLJxzIr588mfCRaeAu6HwFgCpzJpQeHULGnbhvff+mHsyWqcra/mC+Sxzhremr4P5c4M9EtivqlF4/SLpLQ+F1xs3v6IovH6RLH2ekgjvyxvewKr7n8C9v/uD+4RTjhmL2Sd9BuM/dgQ2vP4mfvGrR7D2hX/g8XuWlp6AsEcKrxBgBMMpvLpFC7vw/s9Ncbzzbl7fdZcsDBqUWp9bW+ugukqXa6G9U3gLJRVMOwpvMFwLzUrhLZRU+NoFLrzpl9bMbK6Zwf3iFybj0NoDu5DYvnM3+vaJ3stsFN7wDeqg74jCGzThfecPo/CaXRTuvDuGhm0W3MncrlvnwnzgWMC4jyfx0Y84avvoSqtH4ZUSlMVTeGX8pNEUXilBvfjAhfemFffjA7UH4bOTjkJVZS+9Jw2oZwpvQGBDnJbCq1ucMArvtT+JYdv2DsvN9t3M5QuLzk6gX024dl4opqIU3mJo+d+Wwus/02IyUniLoRWutoELb7ge1/+7ofD6zzTsGSm8uhUKk/CavXXXb7Cweo3tzuymLyO8ruSaf1qpf7ctwKzVPfccvUMj/KgchdcPit5zUHi9s/MjksLrB0WdHBReIXcKrxBgBMMpvLpFC4vwmmUMy1fE0diYwSPHy2nmU7PE4UMfSuIzxybVD46QVo/CKyUoi6fwyvhJoym8UoJ68RReIXsKrxBgBMMpvLpFC4vwZu7GkJ7NNVO6rvO2ia+Z4f3c9ATGj4vuEobsalN4dcc/hVeXP4VXl7+kdwqvhB4ACq8QYATDKby6RQuL8N52ewwbN6VW63YzsYuo7K1bTEUpvMXQ8r8thdd/psVkpPAWQytcbSm8wnpQeIUAIxhO4dUtmqbwmjW7W+otVFU6ePnvFv78vDkAuOMyM71DBjuwbWDEcAcTjw7/QRLFVpPCWywxf9tTeP3lWWw2Cm+xxMLTnsIrrAWFVwgwguEUXt2iaQjvi+ssPPi7GFpbU89uxLaqCnASQHPbz8zPy3FGN7vaFF7d8U/h1eVP4dXlL+mdwiuhxyUNQnrRDKfw6tat1MJrjgpeen3Wy2lt0vvhDycxaqTjbkk2cnj0X0grpLIU3kIoBdeGwhsc20IyU3gLoRTONhReYV04wysEGMFwCq9u0UotvG9stLB8RazLQ5t1u3VDHCyYG+1txoqtJoW3WGL+tqfw+suz2GwU3mKJhac9hVdYCwqvEGAEwym8ukULi/AaCkMpvLqDoQf2TuHVLTqFV5e/pHcKr4QelzQI6UUznMKrW7eghPe+B2Iwa3XNZdnA+I8n8bnpSZglDTfdEkNDxklqbiMHmHVKEiNHlN+LafuqMGd4dcc/hVeXP4VXl7+kdwqvhB6FV0gvmuEUXt26BSG8RnSN8GZeZsnCjBMcjPt4AuaQiSeesvHeexZaWi3U1jqYOD6JuqHls79uoVWl8BZKKph2FN5guBaalcJbKKnwtaPwCmvCJQ1CgBEMp/DqFs1P4TWzt7f8NIat26yce+l+6IMO5nylZ63RzVddCm8+QsF+TuENlm++7BTefITC+zmFV1gbCq8QYATDKby6RZMK7yvrbWzcBOzabeGVV8yMrVnD4P5f6nIAp+3XtQc7+PpCCm9mxSm8uuOfwqvLn8Kry1/SO4VXQo9LGoT0ohlO4dWtm0R4cy1dyDW1mz4q+LgpSUz6VM9ao5uvuhTefISC/ZzCGyzffNkpvPkIhfdzCq+wNpzhFQKMYDiFV7doEuG99ecxvPV2+1xut2cCG+EdNNDB/HkJVFfpPm/Yeqfw6laEwqvLn8Kry1/SO4VXQo8zvEJ60Qyn8OrWzYvwmiOB77w7hs2bASu/7+Lww5P46izO7OaqNIVXd/xTeHX5U3h1+Ut6p/BK6FF4hfSiGU7h1a1bocJrXkjbsiVlt2v/bLvrdc2eChm+637mLl/I+uFJMxIYO6bn7cBQSGUpvIVQCq4NhTc4toVkpvAWQimcbSi8wrpwSYMQYATDKby6RStEeM02YstXdD4O2Iit+3JahvXG48D++ztobbGQSAID+pttyJKU3X2UmMKrO/4pvLr8Kby6/CW9U3gl9DjDK6QXzXAKr27dChHelatsrN9gd7rRtOe6/2ybvF10TgL9ajiTW0xFKbzF0PK/LYXXf6bFZKTwFkMrXG0pvMJ6cIZXCDCC4RRe3aLlE16zlOEn/xPHnj37vs+eeCywH5Wj8PpB0XsOCq93dn5EUnj9oKiTg8Ir5E7hFQKMYDiFV7do+YT3hz+KY/eerutyzV0byXX/OdTBhPFJ7sDgoZQUXg/QfAyh8PoI00MqCq8HaCEJofAKC0HhFQKMYDiFV7doaeF99i+tWPuchTc2maN+gWM/nURrK3D3vbb7Ipqd9SKaEdzp07jzgrR6FF4pQVk8hVfGTxpN4ZUS1Iun8BbA/tGnn8foI4bhoAE1XVpTeAsAWGZNKLy6BTXC++rrNn52u4Nk1iFovSocNLekTNdIrzkxzazkra52cNG3eWKaH5Wj8PpB0XsOCq93dn5EUnj9oKiTg8KbxX3prffgZyt/i2cfuhEH7N/b/XTc9IW4bvE3cMy4URRenXEaql4pvLrlsJIxLLnKQnNL1/vI3nYs/eshgx2cOo/C60flKLx+UPSeg8LrnZ0fkRRePyjq5KDwZnC/b/XT+N5VP3d/QuHVGZBR6JXCq1ulv78cw12/trrsnevO6mbts5vekWHBvFYMHaJ73+XSO4VXt5IUXl3+FF5d/pLeKbxt9P68bj2+ftF1WPzt+fjW4pu6Fd73t+3ARVfcionjRmHezOPBJQ2S4RfNWAqvXt3MiWl/eNrGn1+03KUKXa4s4+1VAZx1Jrce87NiFF4/aRafi8JbPDM/Iyi8ftIsbS4KL4BNb9Xjv0671F22MPDAfjhx/ndzCu+oEXWYd86VqBs8CFdfvBDxWIzCW9rxGoreKLw6ZfjZ8hje/FfHm2juCWnmVsz/OEBNP6ChofO98UU1/2tF4fWfaTEZKbzF0PK/LYXXf6alytjjhXf7jt2Yecb3MXfm8Zh90nH4vzfezim8l1/4Nay4Zw369+uDay49CxXxmFujPU2tpapVez+VFTHEbAtNLQkkktw0v9QFMOzjto2mVq4JLRX7O1Y5eO753BO6RnaP/RRw/GctPPx7B2/9O9XuQ8OAyZ+00Lu6VHfZM/qpiKXm1lvM0XS8Sk6gV9yG+c9+K/mXnL3psDIeQ2syqfJ7b+/KuMozl0unPV541zz5/3Du92/EnC9OcyeKtm7fiQcfeQannDgFX/z8JIw8fIj70port3sbsfqOqzD40IHt9W/YnePNmYBHx35VcVTELOxubEVLgsIbMO4u6c1v+L0qbJc/r2AJvPx34JHHgbf/nXtfXTP6Rw4HvjY32Ptg9g4CVRW2W4zGZv6BT2Nc9O4VQ6tjdiPhHzg0+Jvffw17jT/w1exXofHIZdNnjxfe5x4TdQAAIABJREFU1ze+jcf++EJ7Qd/buh133PsozvjqF3DCcUdj2NBDXeH9wtSJ2Fz/Pja9tQUrl12Mmr77uzFcw1s234WCH4RLGgpGVVTDZ9baaGpKhYwZnfqD3NLrY+6LaOZ/rKx9dd0fO8DIEUnMPoW/+RcFW9CYSxoE8HwI5ZIGHyAKUnBJgwCecmiPF95s/t0taTDre8eOOhynnne1G/Lza85H7+pKCq/yANbonsLrL3VzFPD1y+LYtbstrwNUVQHjP+HgD09bKeE1S3VzbMFg9tmdPjWJiUdTeP2tSvfZKLylIp27HwqvLn8Kry5/Se8U3ix63Qnv9Uu+iQlHfQQN23dh9llLMPjQg7HsikWob2ibkpJUochY84Uz0rV1RxMa+ddaRdKTN6fwyhlmZlhxRwyvvd5xWER6Ird2oIP6d9p+1XaIhPt+Wttsr1HcY3h6mr/FKCAbhbcASAE2ofAGCLeA1BTeAiCFtAmFV1gYLmkQAoxgOIXX36JdfU0MO3dZHTsuZKRPT+pmLmuo7AWMPtLCcVNaUF3l770wW34CFN78jIJsQeENkm7+3BTe/IzC2oLCK6wMhVcIMILhFF5/ivbGRgsbN1l44mkbSHbM3GZnr+nroGG7BfPPYyclccx489JgDA27mv25EWYpigCFtyhcvjem8PqOtKiEFN6icIWqMYVXWA4KrxBgBMMpvPKiGdldviK1tZ/ZYsm2uhfexZd03g2jd2WMwisvgecMFF7P6HwJpPD6gtFzEgqvZ3TqgRReYQkovEKAEQyn8MqLdu/9Max7qfO6XXf5QtZODFWVwHcuoPDKifuXgcLrH0svmSi8Xqj5F0Ph9Y9lqTNReIXEKbxCgBEMp/DKi3bb7TF3OYO52tfnphfsZkjv8Tl2YOAMr5y/JAOFV0JPHkvhlTOUZKDwSujpxlJ4hfwpvEKAEQyn8BZWtFfW21j7XMpe9+8D9O8PvLnJbDnmoKXFwv/9M+OY4LYTgvv3d1BZCRx2qIOPfyyJQbVd+6LwFsY/qFYU3qDIFpaXwlsYp6BaUXiDIht8XgqvkDGFVwgwguEU3vxFe/4FG/c/lDqCNn1lbqNr/r13FWD24E1fkyclMWVS/v10Kbz5+QfZgsIbJN38uSm8+RkF2YLCGyTdYHNTeIV8KbxCgBEMp/Duu2iPPGrjj890lt1M602v1R06JLXrgrlqaoB+NYUdk03h1f3SUHh1+VN4dflTeHX5S3qn8Ero8WhhIb1ohlN4c9dtW4OFTZuA+x6IuTsv5DgJOLVg1+zIAKBuiIMFcxNFDwIKb9HIfA2g8PqKs+hkFN6ikfkaQOH1FWdJk1F4hbg5wysEGMFwCm/Xov1prY2H19iu5BqZNf9vthrLvtoPknCAkSOSmH1K/iUM2TkovLpfGgqvLn8Kry5/Cq8uf0nvFF4JPc7wCulFM5zC27luZmb3uhti7pG/6cuV3rb9dd2fmaOBM7Ydq652sPC0ZMHLGDJ7pPDqfm8ovLr8Kby6/Cm8uvwlvVN4JfQovEJ60Qyn8Hau2zNmdveRrmt2XQFOz/KaJQ5t/96rF3Def7d6PhaYwqv7vaHw6vKn8Oryp/Dq8pf0TuGV0KPwCulFM5zC27lujz9l48mnugqvEdxDD3HQvwZobnXQ2Gi5W5KZF9VybTdW6Gig8BZKKph2FN5guBaalcJbKKlg2lF4g+FaiqwUXiFlruEVAoxgOIW3c9HMfrt33t1VeCdMSGL6Z4tfo5tvSFB48xEK9nMKb7B882Wn8OYjFOznFN5g+QaZncIrpEvhFQKMYHhPFN4X11l44g82Ghos1NQ4OPbTSYwd07FoN/OoYFPSI0Y6+NIXi9+BoZDhQOEthFJwbSi8wbEtJDOFtxBKwbWh8AbHNujMFF4hYQqvEGAEw3ua8G7eAtx0a7xLpRadnej00pl5ea2hAaitdTyvzy1kOFB4C6EUXBsKb3BsC8lM4S2EUnBtKLzBsQ06M4VXSJjCKwQYwfCeJrzdrdE9fmoSE4/2f8lCviFB4c1HKNjPKbzB8s2XncKbj1Cwn1N4g+UbZHYKr5AuhVcIMILhFN5U0U6akei0rKFUpaTwlop07n4ovLr8Kby6/Cm8uvwlvVN4JfS4S4OQXjTDy1149zYCGzZYMEsUzPG/ZmeFQpY0lKqaFN5Skabw6pLO3TuFV7cqFF5d/pLeKbwSehReIb1ohpez8BrJvenWGPbuTW2h61hAzQEOjvoY8PIrFurrgYEDgSmTku5JaRoXhVeDekefnOHV5U/h1eVP4dXlL+mdwiuhR+EV0otmeDkLr7te90m748CIthL1qgTOO8f7YRF+VprC6yfN4nNReItn5mcEhddPmsXnovAWzywsERReYSW4hlcIMILh5SS8ZgeGX90Xx7vvpgoRjwHNCaDrrrrA/DkJ1A3NOD9YqXYUXiXwbd1SeHX5U3h1+VN4dflLeqfwSuhxhldIL5rh5SS8P14ax46dneuQeSJw5iezZuotY8i8Dwqv7veGwqvLn8Kry5/Cq8tf0juFV0KPwiukF83wchBeM7NrrlwvozkOYI4Fzr6y993Vqh6FV4t8ql8Kry5/Cq8ufwqvLn9J7xReCT0Kr5BeNMOjKLxm54V1L9n40zMWduzssNmccps1xVtVCUyfprMFWa4RQuHV/d5QeHX5U3h1+VN4dflLeqfwSuhReIX0ohkeNeE1s7nPPheDOR4418xtd0sYTHUWX9IauiJReHVLQuHV5U/h1eVP4dXlL+mdwiuhR+EV0otmeFSE1wju6kdiaGxMce5ObN1Z3vTnVke7yZOS7vZjYbsovLoVofDq8qfw6vKn8Oryl/RO4ZXQo/AK6UUzPArCa/bTXXp9zJVc8z9mZrc74R1+uIOmZuCgAx0cMshB//5wD5sYVBvO+lB4detC4dXlT+HV5U/h1eUv6Z3CK6FH4RXSi2Z4FIT3jY0Wlq+IdQKcS3jN+txFIdlft9DRQOEtlFQw7Si8wXAtNCuFt1BSwbSj8AbDtRRZKbxCytyHVwgwguFREN4bb45hyzs5pnWd1ExvPA4cdlgS06cmQzuT293QoPDqfmkovLr8Kby6/Cm8uvwlvVN4JfQ4wyukF83wMAqvmdFdv8FCQwNQ3dvC8y9asNJnROTYYiyML6MVOhoovIWSCqYdhTcYroVmpfAWSiqYdhTeYLiWIiuFV0iZM7xCgBEMD5vw5lq+YKZxnTbRzfbdoUMcLJibiCD51C1TeHVLR+HV5U/h1eVP4dXlL+mdwiuhxxleIb1ohodNeFeusrF+Q9fDgDtN8DpARQUwbFgSJ52YRHVVNNlTePXrRuHVrQGFV5c/hVeXv6R3Cq+EHoVXSC+a4WET3ttuj2HjphzrFjLeUoviy2ndjQ7O8Op+byi8uvwpvLr8Kby6/CW9U3gl9Ci8QnrRDA+b8N57fwzrXuoqvGYf3fQ1YXy0Z3UzRwqFV/d7Q+HV5U/h1eVP4dXlL+mdwiuhR+EV0otmuIbwmn111z5nYcsWC61JYEB/oF+NgxHDk6iqsrD8dhsN2zuk96QZ4TkK2O8qU3j9JlpcPgpvcbz8bk3h9ZtocfkovMXxClNrCq+wGnxpTQgwguGlFt4Nr1q470Ebe3ZbMHO2Rmsz53Pnz0mgbqgD8/KauWpqUjJcrheFV7eyFF5d/hReXf4UXl3+kt4pvBJ6nOEV0otmeKmEN/to4DSt7N3GzCzv7FPCdwRwUNWl8AZFtrC8FN7COAXVisIbFNnC8lJ4C+MUxlYUXmFVOMMrBBjB8FIJ7xVXxbG3qfNsbqb0pmd5o77NWLFDgMJbLDF/21N4/eVZbDYKb7HE/G1P4fWXZymzUXiFtCm8QoARDC+F8LbvrZvrPGAAjgNYbcY7ZrSDk0+M7r66xQ4BCm+xxPxtT+H1l2ex2Si8xRLztz2F11+epcxG4RXSpvAKAUYwvBTCa15SW3p9zD0GOMeGY+3UzHZjZ56RKOs1u9lDhMKr+6Wh8Oryp/Dq8qfw6vKX9E7hldDjGl4hvWiGByG8exvh7sCQ+cJZen/dzNlcQ8wcIDHl2CQOqXVQW+tE+hAJLyOAwuuFmn8xFF7/WHrJROH1Qs2/GAqvfyxLnYnCKyTOGV4hwAiG+yW8RnLXPBLDa/+0sGNHx0xueomC+fzJp2xs3mKhsQk4cABw4IEOxo52etSMLmd4w/UlofDq1oPCq8ufwqvLX9I7hVdCjzO8QnrRDPdLeK++No5duzozSC9hmDUziZEjes7OC8WMBM7wFkPL/7YUXv+ZFpORwlsMLf/bUnj9Z1qqjBReIWnO8AoBRjDcD+E1L6WZJQvpF8/SGNLC+8ljkph6HIU31/Cg8Op+aSi8uvwpvLr8Kby6/CW9U3gl9DjDK6QXzXA/hPeZtTYefsTuFsDBBzn4xpk9Z+eFYkYChbcYWv63pfD6z7SYjBTeYmj535bC6z/TUmWk8ApJc4ZXCDCC4X4Ib3czvO62DGbbMQtYcklrBOkEf8sU3uAZ76sHCq8ufwqvLn8Kry5/Se8UXgk9zvAK6UUz3A/hNU/+s9tiePNfVsfbahmya7YiW0zhzTlAKLy63xsKry5/Cq8ufwqvLn9J7xReCT0Kr5BeNMP9El7z9OtftfDoYzbeecdyZ3XTe+6a/XW/cwFneHONEAqv7veGwqvLn8Kry5/Cq8tf0juFV0KPwiukF81wP4XXEDCHTKxcFUN9fYqHkd2TTuQuDd2NDgqv7veGwqvLn8Kry5/Cq8tf0juFV0KPwiukF81wv4U3TcGIb2Ojg0G10eRSqrum8JaKdO5+KLy6/Cm8uvwpvLr8Jb1TeCX0KLxCetEMD0p4o0mj9HdN4S0988weKby6/Cm8uvwpvLr8Jb1TeCX0KLxCetEMzxRes9vCxk2plbf9ahyYU9J4BUuAwhss33zZKbz5CAX7OYU3WL75slN48xEK7+cU3rbatLQm8N77Dejf7wBU9qroVLFHn34eo48YhoMG1HSpJLclC+/gDurO0sL72NMtuO+BWKduJoxPYvo0HhgRFHuTl8IbJN38uSm8+RkF2YLCGyTd/LkpvPkZhbUFhRfAT+94CNf99FftNZo2eRwuPXce+h6wn/uzcdMX4rrF38Ax40ZReMM6kkt4X2nh/fH/JNpndzO753ZiwRaDwhss33zZKbz5CAX7OYU3WL75slN48xEK7+cUXgD3PPQkDjvkYIw+4kP417/fwannXoVTZ52AeaccT+EN79hVu7M//yWOR58A9uxBl6OBzU1ReIMtDYU3WL75slN48xEK9nMKb7B882Wn8OYjFN7PKbw5anPx1bfh7c3v4ralF3QR3ve37cBFV9yKieNGYd7M48ElDeEd3EHc2eYtwE23xt3UZrVuet/cdF/cPzcI6p1zUniDZ7yvHii8uvwpvLr8Kby6/CW9U3iz6Jm1vNNmfQsnHDcB5y2c2Ul4R42ow7xzrkTd4EG4+uKFiMdiFF7J6Itg7IvrrI51u+n30zKs96QZCYwdwxfXgiwthTdIuvlzU3jzMwqyBYU3SLr5c1N48zMKawsKb1ZlLv3xcvzusefw21/8EAcfmHpJzazhvfzCr2HFPWvQv18fXHPpWaiIp15W2rm39KdhVfeKIR6zsLcpgdYk5SroL9eevcBvHgK2bXOwazewpe2AiPZZXgcYMxr4/HQLA/oFfTfMXxGzEIvZaGxOEIYCgcq47f7VRlMLX85UwI+qChsJB2hpJX8N/ub335ZEEq2mCCW+3D9s8vJMgMKbge7G//0Nlv3vb3DXzZfioyPq2j8xwmuuPXsbsfqOqzD40IHtn+3c0+IZvtfA6sp4m/C2qnzpvN53VON+9BPg31tSd+8uY8ixluG8s4EPDIrqE0brviviNmK2ReFVKluvCttdzNPcwj9waJSgslcMyaRD4dWAD6C6MoaWVgetidL/gaNP7847SCkhiGy3FF7A/Y/HNTevwt0PPonbf3Ihjvjw0E4FNcL7hakTsbn+fWx6awtWLrsYNX33d9twDW9kx/4+b9ys1W3YbuPd94DfP2p3Wat7wAFA//6pP+FPGO9g5IjS/8evPMnnfyouacjPKMgWXNIQJN38ubmkIT+jIFtwSUOQdIPNTeEF8L2rfo77Vj+Nm686Dx8c0jFNN/Cgfu463fS2ZGNHHY5Tz7varcjPrzkfvasrKbzBjk+V7CtX2Vi/wcxipa7spbrm14d/0MKcr5R+dl8FSMg6pfDqFoTCq8ufwqvLn8Kry1/SO4UXwLRZ38Zbm9/twvF3v7wKQz4w0BXe65d8ExOO+ggatu/C7LOWYPChB2PZFYtQ39Ak4e8p1nzhzF6wW3c0oZHr6DwxzA7a2wise8nGKxuAjRs7ZLddep2OLciM8P7nDAtjxlB4fYFfZBIKb5HAfG5O4fUZaJHpKLxFAvO5OYXXZ6AlTEfhFcLmkgYhwJCE33Z7LHWIRK69xrK2IKupAb5zro3G1uaQ3H3Pug0Kr269Kby6/Cm8uvwpvLr8Jb1TeCX0uIZXSC8c4WZ298qr295+LUB4Z57s4NMTYti6k8KrUUEKrwb1jj4pvLr8Kby6/Cm8uvwlvVN4JfQovEJ6pQ9/Y6OFVzbY2L4d6NsXGDs69aZ5+jAJ8+9OxvIF90MHcKzUIRMjhiex4KsWjHRReEtfP9MjhVeHe7pXCq8ufwqvLn8Kry5/Se8UXgk9Cq+QXmnDV6+x8eza1B6imdfkSUmsXWujsW05tvuSmgPU1prdFxyMHe2goSH16+oquOunKbylrV1mbxRePfamZwqvLn8Kry5/Cq8uf0nvFF4JPQqvkF5pwy9ZHO86ewugqgqYPjXRcYIagJq+DubPTaJfTdfNxSm8pa1bdm8UXl3+FF5d/hReXf4UXl3+kt4pvBJ6FF4hvdKFb2uwsPT6WLcvpS2+pBWmjZnJNVfd0O5P0aHwlq5uuXqi8Oryp/Dq8qfw6vKn8Oryl/RO4ZXQo/AK6ZU2vNsZ3krgOxcUfkQ0hbe0deMMry7v7N4pvLr1oPDq8qfw6vKX9E7hldCj8ArplTa8uzW8x09NYuLRhZ+URuEtbd0ovLq8Kbzh4k/h1a0HhVeXv6R3Cq+EHoVXSM+f8GfW2mhqe+FszGgn57rbdE/myOAXXophzy4H++1vdmlIYlBtcfdB4S2Ol9+tuaTBb6LF5eMMb3G8/G5N4fWbaHH5KLzF8QpTawqvsBo8eEIIUBjefmBEWx7zAtqZpyf2Kb3CLrlLgxSgMJ7CKwQoDKfwCgEKwym8QoDCcAqvEKBiOIVXCJ/CKwQoCDd76i5fEeuSwWwzNmVS4UsUir0FzvAWS8zf9hRef3kWm43CWywxf9tTeP3lWWw2Cm+xxMLTnsIrrAWFVwhQEL7+VQt33NkmvBl765p9c08+MXWgRBAXhTcIqoXnpPAWziqIlhTeIKgWnpPCWzirIFpSeIOgWpqcFF4hZwqvEKDHcDO7e8eqGJoaASvrIIkxRzo4+T8ovB7Rhj6MwqtbIgqvLn8Kry5/Cq8uf0nvFF4JPb60JqTnPXzlKts9IjjLdd2E4z6exBdO4JIG73TDHUnh1a0PhVeXP4VXlz+FV5e/pHcKr4QehVdIr7hwczDEupcsbNxo4c23LCTM1rk5jHfoEAcL5nKGtzi60WlN4dWtFYVXlz+FV5c/hVeXv6R3Cq+EHoVXSK/w8L+9bOHe+20kWlOG65iD0Oycvgu+tFY41yi2pPDqVo3Cq8ufwqvLn8Kry1/SO4VXQo/CK6RXWPhfX7Zxz71dly+kD//NnOSt6evgzDMSqK4qLLeXVnxpzQs1/2IovP6x9JKJwuuFmn8xFF7/WHrJROH1Qi0cMRReYR340poQYFa4ORji4TUxvLHJQk2NA7Pjwl//ZuP9rV37MbO8Q4c6OOAABzV9gA99yEFtrROo7Jq7oPD6W/Nis1F4iyXmb3sKr788i81G4S2WmL/tKbz+8ixlNgqvkDaFVwgwI3xvI7D0+jgaGzvntC0gmZ7OzfjI/Gh6kccC+3G3FF4/KHrPQeH1zs6PSAqvHxS956DwemfnRySF1w+KOjkovELuFF4hwLbw1WtsPP+ijeZmwIis5f5P6sOKXkBLc9d+elcDF37bvLlW2ovCW1re2b1ReHX5U3h1+VN4dflTeHX5S3qn8ErocQ2vkF4q/MV1Fu57oPOJaZnrcw85BHj77Tb/bZNg8/mCOQnUDc0x9evLXXWfhMIbMOA86Sm8uvwpvLr8Kby6/Cm8uvwlvVN4JfQovEJ6qfCfLY/hzX913V/MrNE1h0qcNCMB27bw3J8t7NgB9K0BPnNsUkV2zf1SeH0pu+ckFF7P6HwJpPD6gtFzEgqvZ3S+BFJ4fcGokoTCK8TOJQ0ygK+st3Hn3XbOJEZ2J306iSmTgjtEwsvdU3i9UPMvhsLrH0svmSi8Xqj5F0Ph9Y+ll0wUXi/UwhFD4RXWgcLrHaB5Se2Wn8bdHRiyjwc2C3lr+jk49+zgDpDweucUXq/k/Imj8PrD0WsWCq9Xcv7EUXj94eg1C4XXKzn9OAqvsAYU3sIBmpPSGhpS7c2WYzfd2rEjg7sSt20JQ3pVbl3AJ6YVfuedW1J4vZLzJ47C6w9Hr1kovF7J+RNH4fWHo9csFF6v5PTjKLzCGlB4CwOY/WJaVZWDvY1WrpOB2xMGfWJaYXfetRWF1ys5f+IovP5w9JqFwuuVnD9xFF5/OHrNQuH1Sk4/jsIrrAGFtzCAV1wVR2NT923TM7ymRf/+wORPJTB2TOl3YCjkaSi8hVAKrg2FNzi2hWSm8BZCKbg2FN7g2BaSmcJbCKVwtqHwCutC4c0P8P7f2nj++a4vpqV3YcjM8MWTk/joqHC9pJb9hBTe/DUPsgWFN0i6+XNTePMzCrIFhTdIuvlzU3jzMwprCwqvsDIU3twA39hoYft24J13LTz9jJ1z6UIsBiQy3kkbGtI1uxRe4ZfE53AKr89Ai0xH4S0SmM/NKbw+Ay0yHYW3SGAhak7hFRaDwtsV4MpVNtZv6JjRzTxEIrP1f/5HEr16AZvrAfOCmsYhEl7KzxleL9T8i6Hw+sfSSyYKrxdq/sVQeP1j6SUThdcLtXDEUHiFdejpwmt2Xli/wUJjI2BmaM21fEXnU9PMz8zyBXOZ7cfMv35uahITjg730oXuhgaFV/ilEYZTeIUAheEUXiFAYTiFVwhQGE7hFQJUDKfwCuH3ZOE1yxZ+eVcMLeZltDaRPeAAYOeOrlCN5KbPUquqBBad04rqKiF8pXAKrxL4tm4pvLr8Kby6/Cm8uvwpvLr8Jb1TeCX0evjRwu7ShfV2h8m2scz1MlpVFdpngadPS2BQrRC8YjiFVxE+AAqvLn8Kry5/Cq8ufwqvLn9J7xReCb0eLLybtwDLfxFH494cs7ltB0ikP4n6jG72E1J4hV8aYTiFVwhQGE7hFQIUhlN4hQCF4RReIUDFcAqvEH5PXNJgjgT+8XVxNDflOBI4Y73usZNTa3QnjE9GdvlCruFB4RV+aYThFF4hQGE4hVcIUBhO4RUCFIZTeIUAFcMpvEL4PVF4X1lvwyxnMItyrczFuW0s00sa5s9JRGbnhWKGAYW3GFr+t6Xw+s+0mIwU3mJo+d+Wwus/02IyUniLoRWuthReYT16mvBurrdw+y9i2LMnBc7M4RrpNbsvuNsvpN9MA3DR+dF9MW1fw4LCK/zSCMMpvEKAwnAKrxCgMJzCKwQoDKfwCgEqhlN4hfB7kvCamd27f213OizC4Et7bubLaiOGJzH7lGhuO5ZvSFB48xEK9nMKb7B882Wn8OYjFOznFN5g+ebLTuHNRyi8n1N4hbXpScJ77U9iMPvuurO5Oa4hbYdH9OvrYOyY9HETQsAhDKfw6haFwqvLn8Kry5/Cq8ufwqvLX9I7hVdCr4ft0nDJ4rh7gEQu4TUvpk2eVF4vp3U3NCi8wi+NMJzCKwQoDKfwCgEKwym8QoDCcAqvEKBiOIVXCL9cZ3jNoRLmqqkB+tWkZmuvuCqOvWZnhixm5oS1BXMTQpLRCafw6taKwqvLn8Kry5/Cq8ufwqvLX9I7hVdCr4xmeNc8ZuOvf7WxZy/gJIHWJGC3sTEzt1MmJbF6jY1nn7Pbjwk2Hw8d6mD2KYmy2nYs35Cg8OYjFOznFN5g+ebLTuHNRyjYzym8wfLNl53Cm49QeD+n8AprUw4zvEZ2//gnu9PMrZnTTb+MZmZ00zsuvLjOwhubUipcNyRZ1mt1uxsaFF7hl0YYTuEVAhSGU3iFAIXhFF4hQGE4hVcIUDGcwiuEH2XhNS+gPbzGwuuv22hu7QqifZcxB5g/tzz31PVSfgqvF2r+xVB4/WPpJROF1ws1/2IovP6x9JKJwuuFWjhiKLzCOkRZeM2uCw3brW5fRMvcVvfM01sxqFYIq0zCKby6haTw6vKn8Oryp/Dq8qfw6vKX9E7hldCL8BrezVuAm26Np54+x2lp7o/bdmQ4oI+Dby3qOS+l5RsSFN58hIL9nMIbLN982Sm8+QgF+zmFN1i++bJTePMRCu/nFF5hbaI6w2t2YVi+Itb+9NnOm7n92AnTExg/rnz31S12CFB4iyXmb3sKr788i81G4S2WmL/tKbz+8iw2G4W3WGLhaU/hFdYiqsJrHttsM9bY1AHAKO0HhwLvvuNg1x4LVZXA0Uendmjg1UGAwqs7Gii8uvwpvLr8Kby6/Cm8uvwlvVN4JfQisKRhbyOwYYPlnpBm9sutG9oxU2uOCr7vfrtdes1xwCed2DMOj5CUncIroSePpfDKGUoyUHgl9OSxFF45Q0kGCq+Enm4shVfIP6wzvGbJwrq/2njprxaSyY5luiOHJzHrFM7YSspO4ZXQk8dSeOUMJRkovBJ68lgKr5yhJAOFV0JPN5bCK+QfRuHNXp/rPmLGIl3uuCArOoVXxk8aTeGVEpTFU3hl/KTRFF7t9sgoAAAUdklEQVQpQVk8hVfGTzOawiukH0bhXbnKxvoN6XPSOh4w7bzz53BPXUnZKbwSevJYCq+coSQDhVdCTx5L4ZUzlGSg8Ero6cZSeIX8wyi8t90ew8ZN5ny0rKvNeDnDKys6hVfGTxpN4ZUSlMVTeGX8pNEUXilBWTyFV8ZPM5rCWyD9nbv2oDWRQL++fTpFhFF4770/hnUv5RBewH1xbcFc7qlbYNlzNqPwSujJYym8coaSDBReCT15LIVXzlCSgcIroacbS+HNw3/P3kZccNktePxPL7otjzxiGG647Gwc2L+v+2tN4f3D2mb869+pXRcyd2AwOzIsv912T1FLX8OHO/jYaAcjR/CFNelXjsIrJSiLp/DK+EmjKbxSgrJ4Cq+MnzSawislqBdP4c3D/mcrf4t7HnwSv7jhu6iu6oUzL1yKusGDsOT8BarCu+Yx4MGHO8vr8VOTmHh0x8/My2vmqqkB+tXw4Ai/vmYUXr9IestD4fXGza8oCq9fJL3lofB64+ZXFIXXL5Klz0PhzcP8v067FNMmj8NpX/6823LNk/8P537/Rrz8xHJYllXyGV6zr+6rGypw34OOu91Y5lVbC3z99NbSj6Ie1iOFV7fgFF5d/hReXf4UXl3+FF5d/pLeKbx56I2bvhCXXXCqK73m+serG/HF07+PZx5chr599iu58C67JY76+u5vevElFF7JF6KQWApvIZSCa0PhDY5tIZkpvIVQCq4NhTc4toVkpvAWQimcbSi8+6iL4zgYdex83HjlIkyaMNpt+frGtzFj3nfx6KprMGjggJJW9c23HSy+uk1oM/bVTd/Eh4dZOP/seEnviZ2RAAmQAAmQAAmQQNgJUHgLmOG9/MKvYeqko3LO8JaywBtec/Cj/8ktvNVVwIIvxzH2yNy7M5TyPtkXCZAACZAACZAACYSJAIU3TzXMGt7jj/0Evjb7BLel5hpes373yqs7ZnDdSV4H+MQ4B8dMTPLFtBJ9s7ikoUSgu+mGSxp0+XNJgy5/LmnQ5c8lDbr8Jb1TePPQ++kdD+FXDz3l7tLQu7oSCy+4VnWXhhfXWbjvgVj7XZvtyGadkoCZ4eVVGgIU3tJw7q4XCq8ufwqvLn8Kry5/Cq8uf0nvFN489HbvacS3Ft+EP6x9yW05angdbrj8HBx8YI37a619eDdtstGcbMaAA7mvruQL4CWWwuuFmn8xFF7/WHrJROH1Qs2/GAqvfyy9ZKLweqEWjhgKb4F12L5zN1paWtsPnEiHaQmvka6tO5rQ2ELhLbCEvjWj8PqG0lMiCq8nbL4FUXh9Q+kpEYXXEzbfgii8vqEseSIKrxA5hVcIMILhFF7dolF4dflTeHX5U3h1+VN4dflLeqfwSugpLmngDK+wcIJwCq8Ang+hFF4fIApSUHgF8HwIpfD6AFGQgsIrgKccSuEVFoAzvEKAEQyn8OoWjcKry5/Cq8ufwqvLn8Kry1/SO4VXQo8zvEJ60Qyn8OrWjcKry5/Cq8ufwqvLn8Kry1/SO4VXQo/CK6QXzXAKr27dKLy6/Cm8uvwpvLr8Kby6/CW9U3gl9Ci8QnrRDKfw6taNwqvLn8Kry5/Cq8ufwqvLX9I7hVdCj8IrpBfNcAqvbt0ovLr8Kby6/Cm8uvwpvLr8Jb1TeCX0KLxCetEMp/Dq1o3Cq8ufwqvLn8Kry5/Cq8tf0juFV0KPwiukF81wCq9u3Si8uvwpvLr8Kby6/Cm8uvwlvVN4JfQovEJ60Qyn8OrWjcKry5/Cq8ufwqvLn8Kry1/SO4VXQo/CK6QXzXAKr27dKLy6/Cm8uvwpvLr8Kby6/CW9U3gl9Ci8QnrRDKfw6taNwqvLn8Kry5/Cq8ufwqvLX9I7hVdCj8IrpBfNcAqvbt0ovLr8Kby6/Cm8uvwpvLr8Jb1TeCX0KLxCetEMp/Dq1o3Cq8ufwqvLn8Kry5/Cq8tf0juFV0KPwiukF81wCq9u3Si8uvwpvLr8Kby6/Cm8uvwlvVN4JfQovEJ60Qyn8OrWjcKry5/Cq8ufwqvLn8Kry1/SO4VXQo/CK6QXzXAKr27dKLy6/Cm8uvwpvLr8Kby6/CW9U3gl9BhLAiRAAiRAAiRAAiQQegIU3tCXiDdIAiRAAiRAAiRAAiQgIUDhldBjLAmQAAmQAAmQAAmQQOgJUHhDX6KuN9jc3IJt23fh4ANrYFlWBJ8gOrecTDrY2rADFRVx9O2zX84b37lrD1oTCfTr2yc6D1ZGd/re1u3oXV2F3tWVZfRU4XoU8z2of28b9u9dhT779+50c+QfbK327G1Ec3Mravrun7Mj8veff0trAhXxWNG8zffknfe34cD+fRGP5Y73/26ZsVACFN5CSYWgneM4uGnFA1i2/D73bvrX9MH/XPHfGH3EsBDcXfndwrN/+TvOvvgGmN9wzDVuzAh868xTMGp4nftr8/MLLrsFj//pRffXRx4xDDdcdrb7Hzte/hJ46tmX8PWLluLGKxdh0oTRbvI3367Hwguuxaa36t1fn/y5T+OSc+d2+xuVv3fUM7KZP8xdfv0v8eAjz7gPPG3yOFz7/bPIvwTlr393Gy67bgXWvvCK29vwYYfhu+d8BSMPH0L+AfJ/8+13MP3L5+P3d/0Yh9Qe2N5Tvv/emP9GfWvxTe2/X1x67lzMnHFsgHfK1MUSoPAWS0yx/Ysvv4avfONy/OKG7+CjIz6I639+L3772LN4dNW1sG3O9PpdmrUv/APvvteAT08YjcbGZixeejvMn+Bv+uEit6ufrfwt7nnwSfzihu+iuqoXzrxwKeoGD8KS8xf4fSs9Ot+G1//ljnvzB4xM4T392z/G/vtV4/ILT8OWd97HzDN+gEsWzcEXpk7s0bz8engz1mee8X3EbBsLZk3Hp8aPxq7de92/WTIX+ftFOnee85fcjIYdu7Dsiv+GZVv4wTW34933t+Hmq84j/4DQz/r6Evz1H6+72bOFd1/jfW9jMz590tn4/+3de3BV1RXH8TUTkeeYYgiPgijWFjvSFh1EoTKACOURKSBICVIDBowgxgo2DARD5GWaaIyIRqAEIhYi4aUpNFpeRSniFIRSbRTKSwYGKEqtigKdztr03iYYyElyz40363v/5WTvvT7r5uZ3ztn38PCogTJ80F2yaet7kjx1jhQvzZRWLWJ9Wi3DVlaAwFtZsRo8/uncV+WDvQdlQdbjbhXHT34q3Qc/KoXz04Nn/TW4vFo/tV7lmjRrnuxa/1t3u2rw6DR3xWv08DhXe/Gm7fLYtBdkz8Y8tpqE6N1w4p+fytCkdHlszL2S/sxiyXriIXeF9/Rnn0vnu8fJkuenyM3tvu9mm5nzshw7fkrmzEwO0ey2h9E7F+On5MjaJRlybatmZTDw9/+9oSd56j5zUqKbbNW6LTJn4UrZsDyb979P/Po3VU+eNfiWDrwVvd8Dd6B2vjFfrryyjltd3/tSXPgdPqinT6tl2MoKEHgrK1aDx+vtksbRjWRK8ojgKm7qllDmqlcNLq/WT61hd+/+I+4EQ1+39kmSGSkPuNCrr/c/PCBDxkyTra/PveR+31qPFMIC9apJQvJs6XLbj92VE/UOBN59B45I/4QpsmnFsxIbc+GK48uFb8ia4reD/QnhUkwOlTF3qRQWbZbe3TvK3gNHJDYmWh4Y1s9tocLf/7fEhrd2yPjU56RHl1tkYJ8ukvnCMhn1i74yOK4r/j7y61aSO4f8qkzgrej9/urrm2RRwTp3chh46cnidde0kAlJ9/q4WoaujACBtzJaNXys3lJp+73WZX6BNARMm5gg/XrcXsOrq93TB67u6tX1Th1uEt1P3a77yDInG4EPxT8WPC0tmsXUbhCfq9Pb6XqCpy8Nubplp3TgDWzvKX1yoX90cvPXuCtgvKov8MjU56Rk72FJGNpbmjVpLH/YuF1+v36bFOXPdrfa9Qok/tV3vtQIR46dlNETM+UH118jb7+7R+rVrSN52ZPkhjYthfe/f+7lBd6KvHV7m/5+BC6G6Or086tRg/ru7zOvb4cAgffb0QdPq9BfIP2i2uRH7gsezxVeT3TVOkj/2OjJxsVfQtAAprcbe3Xt4MbnCm+1mMv8cGC7jl7Nali/nvu3xcuLpVvn9tK/10/lhuu+667wbl6ZE/ySIFd4Q+evI2ngbdk8VlLGDXMDnz//H+l2T7I8dP8Aue3mG/EPLfc3Rhv6YLp07dxext7/c9EvD6ZlLZIt7+yWPxfNlYOHj+Hvk//lrvBe6vOGK7w+NSPEwxJ4Qwzq53C6h7dk3yGZlznRTcMeXj+1L4wd2JerWxf0tmLpl+7h1du9ifH9yhzLHt7q90W/oLZkxZtlBspZsELienaSuLs6uSdiXLyHd3p2vhw/+Ql7eKvP70bQz5uP9h8OfklKA+/tcWNl3MgB7ncB/xBBlzPM51+ckY59k9xTX+684xZ3xN9KDrgvEa7OmyFNmzTG3yf+8gJveXt4S3/eBPbwvvfmAvcIS339bNjj8sshvdjD61OfqjIsgbcqajX0M/9/SsMU+dEPr5ecBYWydv02ntLgUz90P+jk2fNl0sPxwT86OpXuo9bnvs5/pcjtcdSnNOgzYPURWTylwadm/G/PdGAPr86SODFTrmrU0F1l5ykNoXff9f4+iR873Z1gd2x/o6wufkumZS0KfkkW/9Cblx5RA1Ob1s0lIzVJGtSrK8/OL5SNW3fKa4tnuS/N4h96f33+rn6W9I7/tduPq48lCzyP93LeX3z5ldza50F3NySepzSEvjEhGpHAGyLIcAyj+0afz1slufmvuek0dM3LnBD8lno41mBpjiez86VgzYZvlBy42qtXYXSbyZ+27XLH6PN59QkBgcc2WbIKR62l9/DqfPsPHXUnGR8fPeGmH9D7Dpk2ISF4hSUca6rtc+QtWydZuQXBMkvf6cDf3+5/8NFBeTF/jazfssN91nf4SVu3vUEvdvD+98deP2MCz13XGXQL4ZbVczx5B55qElhZ6qMjZNiAHv4slFGrJEDgrRJbzf7Qma++llOf/EuaN43h+bs12wo3u97uOnv2HP/hRA31Qm9B6vN4Gza4sNeXV2gF9PNGHw+nnzfl/e9T+IfW++LR9MT63LnzEn1V+f/TI/7++l88+uW8ddvPsROnpGnMdzjxDm9bPM1G4PXExEEIIIAAAggggAACkSpA4I3UzrFuBBBAAAEEEEAAAU8CBF5PTByEAAIIIIAAAgggEKkCBN5I7RzrRgABBBBAAAEEEPAkQOD1xMRBCCCAAAIIIIAAApEqQOCN1M6xbgQQQAABBBBAAAFPAgReT0wchAACCCCAAAIIIBCpAgTeSO0c60YAAQQQQAABBBDwJEDg9cTEQQgggAACCCCAAAKRKkDgjdTOsW4EEEAAAQQQQAABTwIEXk9MHIQAAggggAACCCAQqQIE3kjtHOtGAAEEEEAAAQQQ8CRA4PXExEEIIIAAAggggAACkSpA4I3UzrFuBBBAAAEEEEAAAU8CBF5PTByEAAIIIIAAAgggEKkCBN5I7RzrRgABBBBAAAEEEPAkQOD1xMRBCCCAAAIIIIAAApEqQOCN1M6xbgQQQAABBBBAAAFPAgReT0wchAACCFRdYM7ClbJ9599l5qREad2yqRuoZN9hmZ6dL0P7d5e7e3Wu+uD8JAIIIIBAhQIE3gqJOAABBBConsDJU6dl4KhUaRZ7tfxubqqcPXdehoxJk5jG0bIwO0XqXBFVvQn4aQQQQACBywoQeHmDIIAAAmEQ2PHXD2XE+FkyfFBPOf3Zv2Xru3tk1cIZ0uTq6DDMzhQIIICAbQECr+3+Uz0CCIRRYPHyYvnN3KVuxoKX0qRd2zZhnJ2pEEAAAbsCBF67vadyBBAIs8CWd3ZLUsozbta1SzLk2lbNwrwCpkMAAQRsChB4bfadqhFAIMwCHx89IQNHTZXe3TvKX3aXyBVRUbIsN00a1K8b5pUwHQIIIGBPgMBrr+dUjAACYRb48szXEj/2SYmKinJfWvvHoaNyT+IT7ukMT00eE+bVMB0CCCBgT4DAa6/nVIwAAmEWSMvKk8KizbLulQxp3fLCNoZlaza4x5KlTxwpg+O6hnlFTIcAAgjYEiDw2uo31SKAAAIIIIAAAuYECLzmWk7BCCCAAAIIIICALQECr61+Uy0CCCCAAAIIIGBOgMBrruUUjAACCCCAAAII2BIg8NrqN9UigAACCCCAAALmBAi85lpOwQgggAACCCCAgC0BAq+tflMtAggggAACCCBgToDAa67lFIwAAggggAACCNgSIPDa6jfVIoAAAggggAAC5gQIvOZaTsEIIIAAAggggIAtAQKvrX5TLQIIIIAAAgggYE6AwGuu5RSMAAIIIIAAAgjYEiDw2uo31SKAAAIIIIAAAuYECLzmWk7BCCCAAAIIIICALQECr61+Uy0CCCCAAAIIIGBOgMBrruUUjAACCCCAAAII2BIg8NrqN9UigAACCCCAAALmBAi85lpOwQgggAACCCCAgC0BAq+tflMtAggggAACCCBgToDAa67lFIwAAggggAACCNgSIPDa6jfVIoAAAggggAAC5gQIvOZaTsEIIIAAAggggIAtAQKvrX5TLQIIIIAAAgggYE6AwGuu5RSMAAIIIIAAAgjYEiDw2uo31SKAAAIIIIAAAuYECLzmWk7BCCCAAAIIIICALQECr61+Uy0CCCCAAAIIIGBOgMBrruUUjAACCCCAAAII2BIg8NrqN9UigAACCCCAAALmBAi85lpOwQgggAACCCCAgC0BAq+tflMtAggggAACCCBgToDAa67lFIwAAggggAACCNgSIPDa6jfVIoAAAggggAAC5gQIvOZaTsEIIIAAAggggIAtAQKvrX5TLQIIIIAAAgggYE6AwGuu5RSMAAIIIIAAAgjYEiDw2uo31SKAAAIIIIAAAuYECLzmWk7BCCCAAAIIIICALQECr61+Uy0CCCCAAAIIIGBOgMBrruUUjAACCCCAAAII2BL4L4I+TjDTEuf3AAAAAElFTkSuQmCC"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#| caption: Plot of location versus measured location. \n",
    "#| label: fig:location_vs_measured_location\n",
    "C = 100 # scaling from meters to centimeters\n",
    "R = 30**2  # 30 centimeter standard deviation\n",
    "N = 200  # number of samples\n",
    "X = rng.uniform(low=0, high=100, size=N) # warehouse length\n",
    "Z = rng.normal(C*X, np.sqrt(R))\n",
    "px.scatter(x=X, y=Z)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is pretty clear that the slope is 100. Can our math corroborate?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "estimated C = 100.03303142771483 (cm/m)\n"
     ]
    }
   ],
   "source": [
    "estimated_C = np.sum(X*Z) / np.sum(X*X)\n",
    "print(f\"estimated C = {estimated_C} (cm/m)\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pretty close! The variance and standard deviation are likewise well estimated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "estimated R=966.2397937984271 (cm^2), stddev=31.084397915971078 (cm)\n"
     ]
    }
   ],
   "source": [
    "estimated_R = np.sum(np.square(estimated_C*X-Z))/(N-1)\n",
    "print(f\"estimated R={estimated_R} (cm^2), stddev={np.sqrt(estimated_R)} (cm)\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extending to multi-dimensional systems\n",
    "\n",
    "> For multi-dimensional systems, if the undertainties for the various dimensions are independent, we can solve the parameter estimation problems independently.\n",
    "\n",
    "In the warehouse example, for the GPS-like sensor, we assumed that the $x$ and $y$ were measured directly and that the noise on them was uncorrelated, and that the noise for each dimension had variance $\\sigma^2$.\n",
    "When this assumption holds, the parameter estimation problem is to estimate the\n",
    "independent 1D scaling factors $C_x$ and $C_y$, and the scalar variance $\\sigma^2$\n",
    "We can use the code above to do so for $x$ and $y$ independently."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multivariate Linear Regression\n",
    "\n",
    "For more general measurement processes or motion models, we typically need heavier machinery. The scalar case above is an instance of *linear regression*, and it can be generalized to the multivariate case. Let us assume a multivariate linear measurement model\n",
    "\n",
    "$$z_k = H x_k + n_k$$\n",
    "\n",
    "where now $z_k\\in\\mathbb{R}^m$, $x_k\\in\\mathbb{R}^n$, and $H$ an $m \\times n$ matrix. Again, $n_k$ is assumed zero mean Gaussian, with covariance matrix $R$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The negative log likelihood of $H$ given the data for our measurement model is\n",
    "\n",
    "$$\\hat{H} = \\arg \\min_H \\sum_k \\text{nll}(H; x_k, z_k) = \\arg \\min_H \\sum_k \\frac{1}{2} (H x_k - z_k)^T R^{-1}(H x_k - z_k).$$\n",
    "\n",
    "Unfortunately, the estimation process for arbitrary $R$ is out of scope for this text. But if we assume *isotropic* noise with covariance matrix $I \\sigma^2$, then the optimization problem simplifies to \n",
    "\n",
    "$$\\hat{H} = \\arg \\min_H \\frac{1}{2\\sigma^2} \\sum_k \\sum_i (H_i x_k - z_{ki})^2,$$\n",
    "\n",
    "where $H_i$ is the $i^th$ row of $H$. Making use of the fact that the scalar $H_i x_k$ can be re-written as $x_k^T H_i^T$, and then setting the derivative with respect to $H_i^T$ to zero we have, for all $i$:\n",
    "\n",
    "$$\\sum_k x_k (x_k^T \\hat{H}_i^T - z_{ki}) = 0,$$\n",
    "\n",
    "yielding the following expression for the $i^{th}$ row of $H$:\n",
    "\n",
    "$$\n",
    "\\hat{H}_i^T = (\\sum_k x_k x_k^T)^{-1} \\sum_k x_k z_{ki}.\n",
    "$$\n",
    "\n",
    "Note the resemblance with the scalar case, but be mindful that the quantity in parentheses is an $n\\times n$ matrix whose inverse will exist only if its rank is equal to the dimensionality of the system, $n$.\n",
    "This rank condition is farily easy to satisfy, and is essentially a requirement that the data set\n",
    "be sufficiently rich.\n",
    "It is interesting to note that the $i^{th}$ row of $\\hat{H}$ is formed as sum of the input vectors $x_k$, weighted by the $i^{th}$ element of the measurement. Finally, we can collect the above for all $i$ and write:\n",
    "\n",
    "$$\n",
    "\\hat{H}^T = (\\sum_k x_k x_k^T)^{-1} \\sum_k x_k z_k^T.\n",
    "$$\n",
    "\n",
    "After this, we can again recover the variance $\\sigma^2$ by\n",
    "\n",
    "$$\\widehat{\\sigma^2} = \\frac{1}{(m N-m n)} \\sum_k \\sum_i |\\hat{H}_i x_k - z_{ki}|^2.$$\n",
    "\n",
    "in which $N$ is the number of data points.\n",
    "Above we would divide by $m N$ as that is how many elements there are in the double sum. However, to get an unbiased estimate we need to subtract $mn$ from the denominator as that is how many entries we already estimated with $H$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multivariate ML Code Example\n",
    "Let us simulate some *2D* states randomly in the warehouse, and multiply with a random $3\\times 2$ $H$ matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = np.array([[10,20],[30,40],[50,60]]) # scaling from meters to centimeters\n",
    "R = 30**2  # 30 centimeter standard deviation\n",
    "N = 200  # number of samples\n",
    "X = rng.uniform(low=[0,0], high=[100,50], size=(N,2)).T # warehouse dimensions\n",
    "Z = rng.normal(H @ X, np.sqrt(R))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The multivariate estimator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H_optimal=\n",
      "[[10.1  19.86]\n",
      " [30.06 39.85]\n",
      " [50.03 59.92]]\n"
     ]
    }
   ],
   "source": [
    "H_transpose = np.linalg.inv(X @ X.T) @ (X @ Z.T)\n",
    "H_optimal = H_transpose.T\n",
    "print(f\"H_optimal=\\n{np.round(H_optimal,2)}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pretty close! And the estimate for the standard deviation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "29.935915630018545"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sqrt(np.sum(np.square(H_optimal @ X - Z))/(3*N-6))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application: Motion Model Estimation\n",
    "\n",
    "> Linear regression can be used to estimate all of the parameters of a general linear motion model.\n",
    "\n",
    "Assuming we have a dataset of triplets $x_k, u_k, x_{k+1}$, i.e., the current state $x_k\\in\\mathbb{R}^n$, the control $u_k\\in\\mathbb{R}^p$, and the next state $x_{k+1}\\in\\mathbb{R}^n$. Note the similarities with the reinforcement learning setup in Section 3.6! If we assume a *linear process model* then we should have:\n",
    "\n",
    "$$x_{k+1} = A x_k + B u_k + w_k$$\n",
    "\n",
    "where typically the covariance on the *process noise* $w_k$ is denoted as $Q$. We only have to note that we can write this in the form needed by the multivariate ML estimator, as\n",
    "\n",
    "$$\n",
    "x_{k+1} = \\begin{bmatrix} A & B \\end{bmatrix} \\begin{bmatrix} x_k \\\\ u_k \\end{bmatrix} + w_k\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we generate some random states, controls, and next states:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 500  # number of samples\n",
    "X = rng.uniform(low=[0,0], high=[100,50], size=(N,2)).T # states\n",
    "U = rng.normal(scale=5, size=(N,2)).T # controls, zero mean, 5m standard deviation\n",
    "A = np.eye(2)\n",
    "B = np.eye(2)\n",
    "R = 4  # 2 meter standard deviation\n",
    "Z = rng.normal(A @ X + B @ U, np.sqrt(R)) # the next state\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The multivariate estimator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4, 500)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "XU = np.vstack([X,U])\n",
    "XU.shape\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A_optimal=\n",
      "[[ 1. -0.]\n",
      " [-0.  1.]]\n",
      "B_optimal=\n",
      "[[ 1.01 -0.02]\n",
      " [-0.02  0.99]]\n"
     ]
    }
   ],
   "source": [
    "AB_transpose = np.linalg.inv(XU @ XU.T) @ (XU @ Z.T)\n",
    "AB_optimal = AB_transpose.T\n",
    "A_optimal, B_optimal = AB_optimal[:,:2], AB_optimal[:,2:]\n",
    "print(f\"A_optimal=\\n{np.round(A_optimal,2)}\")\n",
    "print(f\"B_optimal=\\n{np.round(B_optimal,2)}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pretty close! And the estimate for the standard deviation (note the changed denominator!):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.9365058899020415"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sqrt(np.sum(np.square(AB_optimal @ XU - Z))/(2*N-8))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Expectation Maximization\n",
    "\n",
    "> We can estimate unknown parameters even if we don't have access to ground truth state values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above assumed that we had perfect knowledge of the state $x_k$ at all times. But what if we only had access to the controls $U$ and measurements $Z$? A natural idea is to try and estimate the states with the techniques from Section 4.4, and then somehow use the posterior density $p(X|U,Z)$ to estimate any unknown parameters $\\Theta$. However, the estimation problem itself *uses* the parameters, so we have a chicken and egg problem!\n",
    "\n",
    "The **expectation-maximization** or **EM algorithm** is an algorithm that allows us to estimate measurement and motion models *even when the state is only observed from measurements*. We do need *some* initial estimate $\\Theta^1$ for the parameters $\\Theta$, and the idea is to then alternate between two steps:\n",
    "- E-step: calculate the density $p(X|U,Z;\\Theta^t)$ on hidden states $X$, *given* the guess $\\Theta^t$.\n",
    "- M-Step: maximize the *expected* log-likelihood with respect to $\\Theta^{t+1}$ under this density."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We do not discuss the detailed derivations here, but we do want to mention that a particularly simple variant of EM is to use a particle filter in the E-step. We can then simply \"pretend\" that the samples form a large dataset, and perform the above estimation procedures in the M-step. This is an instance of a so-called **Monte Carlo EM** algorithm, and the amazing fact is that this then actually converges to an ML estimate for the unknown parameters $\\Theta$."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "S46_logistics_learning.ipynb",
   "provenance": []
  },
  "interpreter": {
   "hash": "c6e4e9f98eb68ad3b7c296f83d20e6de614cb42e90992a65aa266555a3137d0d"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.18"
  },
  "latex_metadata": {
   "affiliation": "Georgia Institute of Technology",
   "author": "Frank Dellaert and Seth Hutchinson",
   "title": "Introduction to Robotics"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}