Motion model/ex4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 4\n",
    "\n",
    "**_DISCLAIMER:_** The notation used in this exercise follows the one of the Probabilistic Robotics book (refer to Chapter 5.4 in case you have doubts).\n",
    "\n",
    "## 4.1 Inverse motion model\n",
    "\n",
    "The odometry model uses the _relative motion information_. The odometry readings are $u_t = [{\\overline{x}}_{t-1} , {\\overline{x}}_{t}]$, where $\\overline{x}_{t-1}$ and  $\\overline{x}_t$ are poses in a robot-internal coordinate frame (different from the map).\n",
    "\n",
    "The function `inverse_motion_model` takes as input an odometry reading $u_t$ that consist in:\n",
    "\n",
    "- the initial pose of the robot in the odometry coordinate frame $\\overline{x}_{t-1} = [\\overline{x},\\overline{y},\\overline{\\theta}]$\n",
    "- the estimated pose of the robot in the odometry coordinate frame $\\overline{x}_t = [\\overline{x}',\\overline{y}',\\overline{\\theta}']$\n",
    "\n",
    "The output is the relative motion $\\delta_{rot1}, \\delta_{trans}, \\delta_{rot2}$.\n",
    "\n",
    "Implement the function `inverse_motion_model` and verify that it is correct for some test input. **[2.0]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.7853981633974483, 1.4142135623730951, -0.7853981633974483)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from ex4 import *\n",
    "import math\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "inverse_motion_model([[0,0,0],[1,0,0]]) #Pure translation\n",
    "inverse_motion_model([[0,0,0],[1,1,0]]) #translation + rotation\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.2 Odometry-based motion model\n",
    "\n",
    "The function `motion_model_odometry` computes the posterior $p(x_t | u_t, x_{t-1})$ from odometry readings.\n",
    "\n",
    "This function takes as input:\n",
    "\n",
    "- the initial pose of the robot $x_{t-1} = [x,y,\\theta]$ _(**map** coordinate frame)_\n",
    "- the hypothesized (or query) final pose $x_{t} = [x', y', \\theta']$ _(**map** coordinate frame)_\n",
    "- the odometry readings $u_t = [\\overline{x}_{t-1} \\overline{x}_t]$ _(**odometry** coordinate frame)_\n",
    "- the noise parameters $\\mathbf{\\alpha} = [\\alpha_1, \\alpha_2, \\alpha_3, \\alpha_4]$\n",
    "\n",
    "The output is the probability $p(x_t | u_t, x_{t-1})$\n",
    "\n",
    "Assume that a robot starts at pose $x_0 = [2.0, 3.0, 0.0]$ in the map frame and moves according to a motion model with $\\mathbf{\\alpha} = [1.0, 1.0, 0.01, 0.01]$.\n",
    "\n",
    "The robot excecutes one motion command and the odometry readings are:\n",
    "\n",
    "1. $\\overline{x}_0 = [0.0 , 0.0 , 0.0   ]$\n",
    "2. $\\overline{x}_1 = [0.5 , 0.0 , \\pi/2 ]$\n",
    "\n",
    "Implement the `motion_model_odometry` function and verify that it is correct for some test input. **[1.0]**\n",
    "\n",
    "---\n",
    "\n",
    "Consider a 150x150 grid map the world with a resolution of 0.01, centered in the original position of the robot.\n",
    "\n",
    "Plot the posterior $p(x_t | u_t, x_{t-1})$ for all possible $[x, y]$ values from the grid. **[2.0]**\n",
    "\n",
    "**Note that** the query input is a position, not a pose. Therefore, to plot the posterior belief over the gridmap, you can assume the term $\\hat{\\delta}_\\mathrm{rot2}$ to be zero and, for each position, integrate over all possible orientations. This can be implemented by considering $p_3 = 1.0$ in the equations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.0, 6.804138174397719, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAABECAYAAACCuY6+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAGgklEQVR4nO3dXYicVx3H8e/PRKtNoY02YEyKiSiNiyIxi1YDRZoULEp6oUILSiuW9cLaWgRfwQuvoogvFyKERClaaiEWjRJ8I+2VELrbRrSJsTFKkxjtNn3RCrZGf17ME3cYJ9nNPs/O2c75fWDZ5+XkOX8Omd/MPjNzjmwTERHj7yWlC4iIiNFI4EdEVCKBHxFRiQR+REQlEvgREZVI4EdEVKJV4Et6paRfSHqs+b36PO3+LelQ87OvTZ8REbE4avM5fElfBp6yvVPSZ4DVtj89pN1zti9rUWdERLTUNvCPAu+yfVrSWuBB21cPaZfAj4gorG3gP2P7imZbwNPn9gfanQUOAWeBnbZ/eJ7rTQFTAKtWrdqyadOmRdfWhZmZmaL9n7Nly5bSJUTEi8TMzMyTttcMOzdv4Ev6JfDqIac+D9zdH/CSnrb9f/fxJa2zfUrS64ADwDbbf7hQv5OTk56enr5gbUut9xxWXqa/iIiFkjRje3LYuZXz/WPb2y9w4b9KWtt3S+eJ81zjVPP7uKQHgc3ABQM/IiK61fZjmfuAW5rtW4AfDTaQtFrSJc32lcBW4HDLfiMi4iK1DfydwPWSHgO2N/tImpS0u2nzRmBa0q+BB+jdw0/gR0SM2Ly3dC7E9hlg25Dj08BtzfavgDe36SciItrLN20jIiqRwI+IqEQCPyKiEgn8iIhKJPAjIiqRwI+IqEQCPyKiEgn8iIhKdBL4kt4t6aikY828+IPnL5F0X3P+oKQNXfQbEREL1zrwJa0AvgncAEwAN0uaGGj2EXpTJ78e+Brwpbb9RkTExeniFf7bgGO2j9t+Afg+cONAmxuBu5vtvcA2LZe5hyMiKtFF4K8DTvTtn2yODW1j+yzwLPCqwQtJmpI0LWl6dna2g9IiIuKcZfWmre1dtidtT65ZM3TBloiIWKQuAv8UcFXf/vrm2NA2klYClwNnOug7IiIWqIvAfwh4g6SNkl4G3ERvYZR+/QulvB844KzbFxExUq3mw4fePXlJtwM/A1YA37b9qKQvAtO29wF7gO9KOgY8Re9JISIiRqh14APY3g/sHzj2hb7tfwIf6KKviIhYnGX1pm1ERCydBH5ERCUS+BERlUjgR0RUIoEfEVGJBH5ERCUS+BERlUjgR0RUYlQLoNwqaVbSoebnti76jYiIhWv9Tdu+BVCupzc18kOS9tk+PND0Ptu3t+0vIiIWZ1QLoERERGFdzKUzbAGUtw9p9z5J1wK/B+6yfWKwgaQpYKrZfU7S0Za1XQk82fIaxXW0ONhYjEVHMhZzMhZzxmUsXnu+E51MnrYAPwbutf28pI/SW+7wusFGtncBu7rqVNK07cmurvdilrGYk7GYk7GYU8NYjGQBFNtnbD/f7O4GtnTQb0REXISRLIAiaW3f7g7gSAf9RkTERRjVAih3SNoBnKW3AMqtbftdoM5uD42BjMWcjMWcjMWcsR8LZaXBiIg65Ju2ERGVSOBHRFRibAN/vukeaiHpKkkPSDos6VFJd5auqSRJKyQ9IuknpWspTdIVkvZK+p2kI5LeUbqmUiTd1Tw+fivpXkkvL13TUhjLwO+b7uEGYAK4WdJE2aqKOQt80vYEcA3wsYrHAuBO8imxc74B/NT2JuAtVDouktYBdwCTtt9E78MnN5WtammMZeCT6R7+x/Zp2w8323+n96BeV7aqMiStB95D77sgVZN0OXAtsAfA9gu2nylaVFkrgVdIWglcCvy5cD1LYlwDf9h0D1WGXD9JG4DNwMHCpZTydeBTwH8K17EcbARmge80t7h2S1pVuqgSbJ8CvgI8DpwGnrX987JVLY1xDfwYIOky4AfAJ2z/rXQ9oybpvcATtmdK17JMrATeCnzL9mbgH0CV73VJWk3vDsBG4DXAKkkfLFvV0hjXwJ93uoeaSHopvbC/x/b9pespZCuwQ9Kf6N3iu07S98qWVNRJ4KTtc3/t7aX3BFCj7cAfbc/a/hdwP/DOwjUtiXEN/Hmne6iFelNt7gGO2P5q6XpKsf1Z2+ttb6D3/+GA7bF8FbcQtv8CnJB0dXNoGzC4hkUtHgeukXRp83jZxpi+gT2q2TJH6nzTPRQuq5StwIeA30g61Bz7nO395UqKZeLjwD3Ni6LjwIcL11OE7YOS9gIP0/tU2yOM6TQLmVohIqIS43pLJyIiBiTwIyIqkcCPiKhEAj8iohIJ/IiISiTwIyIqkcCPiKjEfwHiBhP17omTNgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "u = [[0.0,0.0,0.0],[1.0,0.0,0.0]] # u = (x,x')T\n",
    "x_0 = [0.0,0.0,0.0]\n",
    "alpha = [1.0,1.0,0.01,0.01]\n",
    "p = []\n",
    "x= []\n",
    "gridmap =[[0 for i in range(10)]]\n",
    "# for i in range(100):\n",
    "#     if i <=50:\n",
    "#         query = [1.0-(0.01*i),0.0,0.0]\n",
    "#         x.append(query[0])\n",
    "#         probs = motion_model_odometry(x_0, query, u, alpha, )\n",
    "#         print(x,probs)\n",
    "#         p.append(probs)\n",
    "#     else:\n",
    "#         query = [1+(0.01)*(i-50),0.0,0.0]\n",
    "#         x.append(query[0])\n",
    "#         probs = motion_model_odometry(x_0, query, u, alpha)\n",
    "#         p.append(probs)\n",
    "for i in range(len(gridmap[0])):\n",
    "    gridmap[0][i] = motion_model_odometry(x_0,[i,0,0], u, alpha)\n",
    "print(gridmap)\n",
    "plot_map(gridmap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQEAAAD8CAYAAAB3lxGOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbU0lEQVR4nO2da4wk13Xff2d6Hj09s8slRYah+fBSNq2AMZCIWDAM5CiBmRiUImsdWBAoGLEsEyCC0LEUObBI64P90YwTKTYSyFAsRlJA6xFZhglDjkXTsgUDJqMlRYkvW1zRokRiyeVjZ+c9PdNz8qHrFO/Udu/sTk+/5v5/QKOrqqun71TV/d9zzj33XnN3hBD5MjHsAgghhotEQIjMkQgIkTkSASEyRyIgROZIBITInL6JgJndZmZ/a2Ynzezufv2OEKI3rB95AmZWA74D/CvgBeAbwPvc/el9/zEhRE/0yxK4GTjp7s+5exP4PHC8T78lhOiByT793auBHyT7LwD/pNvJl19+uR89erRPRRFCADz66KOvuvsV1eP9EoFdMbM7gTsBrrvuOk6cODGsogiRBWb2fKfj/XIHXgSuTfavKY6VuPsn3f2Yux+74opzxEkIMSD6JQLfAG4ws+vNbBq4HXigT78lhOiBvrgD7r5lZr8E/ClQA+5z96f68VtCiN7oW0zA3b8CfKVff18IsT8oY1CIzJEICJE5EgEhMkciIETmSASEyByJgBCZIxEQInMkAkJkjkRAiMyRCAiRORIBITJHIiBE5kgEhMgciYAQmSMRECJzJAJCZI5EQIjMkQgIkTkSASEyRyIgRObsWQTM7Foz+5qZPW1mT5nZB4vjl5nZg2b2bPF+6f4VVwix3/RiCWwBv+LuNwK3AHeZ2Y3A3cBD7n4D8FCxL4QYUfYsAu5+yt0fK7aXgGdor0F4HPhMcdpngJ/psYxCiD6yLzEBMzsKvBV4BLjS3U8VH70EXLkfvyGE6A89i4CZzQN/AHzI3RfTz9zdAe/yvTvN7ISZnXjllVd6LYYQYo/0JAJmNkVbAO539y8Xh182s6uKz68CTnf6rhYkFWI06KV3wIBPAc+4+8eSjx4A3l9svx/4o70XTwjRb3pZi/BtwL8FnjCzx4tjvwb8JvBFM7sDeB54b08lFEL0lT2LgLv/FWBdPr51r39XCDFYlDEoROZIBITIHImAEJkjERAicyQCQmSORECIzJEICJE5EgEhMkciIETmSASEyByJgBCZIxEQInMkAkJkjkRAiMyRCAiRORIBITJHIiBE5kgEhMgciYAQmSMRECJz9mPxkZqZfdPM/rjYv97MHjGzk2b2BTOb7r2YQoh+sR+WwAdpr0MY3At83N1/FDgD3LEPvyGE6BO9rkB0DfCvgd8r9g34SeBLxSlakFSIEadXS+C/Ab8KbBf7bwIW3H2r2H+B9krFQogRpZdlyN4FnHb3R/f4fS1IKsQI0Isl8Dbg3Wb2PeDztN2A3waOmFmsbHQN8GKnL2tBUiFGgz2LgLvf4+7XuPtR4Hbgz93954CvAe8pTtOCpEKMOP3IE/gI8GEzO0k7RvCpPvyGEGKf6GVV4hJ3/wvgL4rt54Cb9+PvCiH6jzIGhcgciYAQmSMRECJzJAJCZI5EQIjMkQgIkTkSASEyRyIgROZIBITIHImAEJkjERAicyQCQmSORECIzJEICJE5EgEhMkciIETmSASEyByJgBCZIxEQInMkAkJkTq/LkB0xsy+Z2d+Y2TNm9k/N7DIze9DMni3eL92vwgoh9p9eLYHfBv6vu/8D4B/RXpj0buAhd78BeKjYF0KMKL0sQ3YJ8HaKdQXcvenuC8Bx2guRghYkFWLk6cUSuB54BfhfZvZNM/s9M5sDrnT3U8U5LwFX9lpIIUT/6EUEJoGbgE+4+1uBFSqmv7s74J2+rAVJhRgNehGBF4AX3P2RYv9LtEXhZTO7CqB4P93py1qQVIjRoJcFSV8CfmBmbykO3Qo8DTxAeyFS0IKkQow8va5F+B+A+81sGngO+ABtYfmimd0BPA+8t8ffEEL0kZ5EwN0fB451+OjWXv6uEGJwKGNQiMyRCAiRORIBITJHIiBE5kgEhMgciYAQmSMRECJzJAJCZI5EQIjMkQgIkTkSASEyRyIgROZIBITIHImAEJkjERAicyQCQmSORECIzJEICJE5EgEhMkciIETm9Log6X80s6fM7Ekz+5yZ1c3sejN7xMxOmtkXipmIhRAjSi9rEV4N/DJwzN1/HKgBtwP3Ah939x8FzgB37EdBhRD9oVd3YBKYNbNJoAGcAn6S9mpEoAVJhRh5elmB6EXgvwDfp135zwKPAgvuvlWc9gJwda+FFEL0j17cgUtpL0N+PfBDwBxw20V8XwuSCjEC9OIO/Evg79z9FXffBL4MvA04UrgHANcAL3b6shYkFWI06EUEvg/cYmYNMzPeWJD0a8B7inO0IKkQI04vMYFHaAcAHwOeKP7WJ4GPAB82s5PAm4BP7UM5hRB9otcFSX8d+PXK4eeAm3v5u0KIwaGMQSEyRyIgROZIBITIHImAEJkjERAicyQCQmSORECIzJEICJE5EgEhMkciIETmSASEyByJwB5xd9x92MUQomd6GkCUG50qvbvTHkktcqH6HIz7/ZcI7ELc8E4tf3rzx/1BELvTzfIbd1GQCJyHtOJfiPk/bjdf7I1uz0Hc/3GzDiUCHehU8be3t3fcfDPbcbOr++Jgkj4bnRhHIZAI7EIIgbuzvb1dHk8r/8TEBNvb25jZ2Nx4cXGkz0Hsp0QjENtxzjg8DxKB85De6Kp1MDEx0fE8cfDoJACd7vm4Vf5AXYRdOJ8AhFXQaVuCcLC52DjROCBLoAth3lVv+Pb29g7TP9yB+E7sj1NLILrTyR3crfKnruI4sKslYGb3mdlpM3syOXaZmT1oZs8W75cWx83MfqdYjPTbZnZTPwvfL3a7eakYdHodlBYid1IB6Gb1HYT7fCHuwKc5d2Whu4GH3P0G4KFiH+AdwA3F607gE/tTzOFQDfSlN73Vap1T6Q/awyE6C8GFWgLjwq4i4O5fB16vHD5Oe7FR2Lno6HHgs97mYdqrEV21T2UdKLvdyFarxebmJltbW2xtbdFqtUphiHcxvnSy9jpZA1XGTQBg7zGBK939VLH9EnBlsX018IPkvFiQ9BRjSOr3p68gKnwwMTGxY1/dhuNJt5Y/zRWp3tPqfR6ne95z74C3r8pF277jsCBpp4qfEg9GWAKd4gNi/DifC3AQM0f3KgIvh5lfvJ8ujr8IXJucdyAWJO1mEbh76RY0m002NzdLtyB1DxQfGB/2EvQdVwsg2KsIPEB7sVHYuejoA8DPF70EtwBnE7dhLInuvvO5BiEAzWaTZrO5I0YQ22L02U0AOqWOp92B4+r67RoTMLPPAf8CuNzMXqC99uBvAl80szuA54H3Fqd/BXgncBJYBT7QhzIPnGpqcBxLb3qr1SozxdJxBPHa3t7ekWUoRovzdQGmxzpV8nHLC6iyqwi4+/u6fHRrh3MduKvXQo0i1eSg9AWwtbW14/xOQcGpqamxfVAOOlUrIPZD3FMRr1oA6fs4oozBC6CaGRhjB2q1Gq1Wq3zf2toqH5zJycmy9YhtM6NWq8kiGDHSbt1qxY9td6dWq3V0CVOXcRyRCFwEExMTO8z8EIL05m9ubpbb1ZTitPJLCIZP+PfV1j+2I64DO+9X1SVMj40jEoGLICpyWAJhItZqtdIa2NjYOKc/OboK0wdGIjB8ooWvWgHpdqvV2jU4PM4CABKBiyKtvGlFr9VqTE5Osr29zeTkJJubm2U2YavVYmpqqhSCNH8gjSmIwZJW+Gq2Z/oKSy+EIH0fdzcgkAjsgdQVCF8xXlNTU2V8YH19HXjD7AwrIrZrtVq5LQZD1QXoJAapS5e6fN0sgnFHInCRpPGAcAcmJyfLYGAE/tydZrNZugMRS0itgPh7EoHBUm39I+tze3u7TPiC3XNEJAKZE604UEb/gdIScHeWl5dZWVkpXYTt7e3y8+3tbaanp+UaDJBO5n5U/HDdtra2dlh6IeqdXgdBAEAisGfSByVak8nJyXMq+fr6Os1ms6zsMzMz50xUEl2H8RL7S1zntAs3zeoMEUjvYycBqMYCDgoSgR6IyhsBpMgJSN0DoHzYwh1IE0+ipyE4aA/YsEl7AFLfPx0CHi5AmuDVqdU/KIHAKhKBHogHIip76v8DOwKEi4uLNJtNZmZmSstgZmaGmZmZcntqagp3Z2pqSq7BPpGa/qnJn84Fsbm5WYrxzMzMDksgrIHUKpAIiHOoxgempqYAdkSat7a2WF5eZnNzs8wliIFHEZgKVyF6HCYnJw/cAzcoOvn/MbgrFYGwAMKCS1+d4gEH8X5IBPaBND4Q+QLRokfFnp6exszKBzB6ECImENlq8ZDFOIODaoL2i3TQT7X1bzab52yHa5YGArsFAw/qPZAI7BPxkKQVOx7E2A5rYGlpiY2NDRqNBrOzs+VDOTMzU1oP4SpMT08rYHiBVCd5SVv8GOodxzc2NoC2gIcrNjU1VVoBnXoHDioSgX3EzEpXIBWECBpubW2VD9PCwgKLi4ulIMzPzzMzM1M+sBsbG8zOzjI7O1s+mGFNiHNJzf4I9MV17DTxC7Strenp6VIA4pVaBNWxIQcRicA+Ez0GQNlN6O5ly+7uNBoNVlZWyklIzp49y9bWFo1GY0drFiIS34u/rQSjNumY/2rrH6+IA6SxmHDVouUPNy6XGEAViUAfiIpanYgiBqNAezai6DE4c+YM6+vrpTCEOxCmbLgKs7Oz5YOruQneSP1Ng3zp7E5R8dPtqPCpCxDWQCoGOblgEoE+MTExwfT0NMCOVOJ4uNLWKLILFxYWWFlZ4dChQ8zOznL48GHm5uZKAdjY2GB6erq0KtIodi6EpZTO5xjCGYIZ+1URmJ2dLa9dxFsieNvJIsgFiUAfiRyCsApmZmbK7eibBsoHdW1tjdXVVVqtFmtra+WDPDs7Wz70ETcIVyNGKaYJLQeNNOJfbf0jyBdBv9hOA4GRxFVt+asxgJxcgBSJQJ9JA0tppmAEp8yszCWo1Wq89tprpdm6sbFR9iLMz8+XbkKIQrRm9Xq9fJgPYnditfWPil8VgdiOz2IAV7T89Xq9vGbxSq2pnFyAFInAAJiYmCh9+HjooiehXq+XmWqNRoNarcbZs2dZXV3lzJkzzM/PU6/XmZub45JLLqHRaDA3N8fc3Fz5UDcajfKhDl83hjWPK936+VP/PjX5o+Kvr6+ztrZGrVYrr1sIQFynuO5R+VMBzZELmW34PuBdwGl3//Hi2G8BPw00ge8CH3D3heKze4A7gBbwy+7+p/0p+niR9hpE3kCY7zFJaYxDAMqkozBt19bWaDabzM7O0mg0SkGINOR6vV5aBSEE0SWZJsOMKtXpvaKHJO3nrwb5wuTf2NhgfX29FIzoTm00GuX1iOsTIlDNCsyZC7EEPg38d+CzybEHgXvcfcvM7gXuAT5iZjcCtwP/EPgh4M/M7MfcXRPvww6/PRWBSCiqikCr1WJhYYH19fWyr7ter1Ov11lfX6fRaNBoNHaIQAhDtHZh8qYClLZ4w2r90vn70wE+aevfKdIfohCjM0MAIv16cnKSubm5Mgg4Ozt7jpWUpgYf9ESgC+FCphz/upkdrRz7arL7MPCeYvs48Hl33wD+zsxOAjcDf70/xR1/orJHElCYomkX4MzMDIcPHy5bspWVFRYXFzlz5kwZ5Ao3oV6v7zj38OHDO6yCiBdUXYW0FYwBUIOgOolHbKc5/dXWPzX7o+Kvra2Vn7VarbJH5ZJLLtnRo5JaADMzM+X/ri7WN9iPu/+LwBeK7atpi0IQC5KKChEfiFYoHs7YjnThiYkJlpeXy2zBqAhnzpwpH+bl5eXygV9YWCjjBfEeFSD84BCFaou42xRa3WbWrc6PkG7HPArdovuxnXbxVbv70oqfHotWfX5+nkOHDpXxkrAC0kBgtSdAAvAGPYmAmX0U2ALu38N37wTuBLjuuut6KcbYEqZodQ7CdHgywPT0dOlCrKyssLS0xPLyMs1mk1qtRrPZLCt6VPwIhFVbxU5WQafMuU6TaXbreahW/E5detWJPNLx/Gk6b9rPHxU/9f3TacDD2qm2/uEChBDE/yXzvzN7FgEz+wXaAcNb/Q0H76IWJAU+CXDs2LGsV+yMyri9vV26BpEYFC36oUOHmJubY2lpiaWlJRqNBsvLy6yurvLqq6+W4+HTShAVPypDah5Ha1kdONNtGO2F5iB0CvBVff00vz/171utFuvr6+e4AGEVxP/XaDQ4cuQI8/PzZfdpCF8aGA03IFweVf7O7EkEzOw24FeBf+7uq8lHDwC/b2Yfox0YvAH4fz2XMgOiokXLlSazpK344cOHWVxcpF6vs7S0xMrKCjMzM2XliQlMgLJSRKsYMYfYT7vLQgjiWFScEINOrkD1WCdLoNtcfp3M/nSgTxwPFyasmni/9NJLywof1k81GJp2x0oAurPXBUnvAWaAB4sH4GF3/3fu/pSZfRF4mrabcJd6Bi6ctBsxNcPTfuy0dUvN3bAKtre3S7M5fOhozaOCRPAxKn6aNpsOXU4nNknjBcA5C66mVBf1SEUg9tOMv/RVjSuk/f1HjhwpRSBGXYa4xf8T72nrL////FjaVTMsjh075idOnBh2MUaSqNDNZpPV1dWyxV9cXGR9fZ3V1VVef/11FhcXWV1dLYcnr6+vs7S0xOLiIsvLyywuLpbdh+EOREtZdQc6uQLV1jTmVazGDGBnjKC6pFe1dyAVgImJidLtqdfrzM/Plz0fkSwVFX9ubq4UxBC31IJSxT8XM3vU3Y9VjytjcMSJacmjBZ+bm2NjY4N6vV6a//Pz86yurrKyssLCwgLLy8usra2VE5jEmITwvdfW1lhYWCh973ARqpWoGkmPcQ+wc1HWTkIQVN2DeE/779NWfX5+vozwhwiEO5RmRsb3w4VJg34SgItDIjDiRGWMKcrSwS7paMIwmWdnZ1lcXGRtba1MJlpbW2NlZYXl5WXW19dLF6JWq5WrJDWbTQDW1taAnSv1VCtxtOpVK6IahU97GtJyR3nDjI/uvRCB6POPGEBq9odgpWm/6d8WF49EYEyI2EAMPZ6cnCzN6Egd3tjYKM3mjY2NsichLIHFxcXSelhZWSldi9XV1TL5JoQi8vHTgToRvEtFIB2bH5UzKm6Y8WmXZTWYF+KVtv7R4odQpH+7U7em6A2JwBgSI+Ni1qJ6vV5W1KjAzWazjBnEdsQKVldXyzyDEIFwDcJl6JSxl4pAjIiMYFwamAvzPM1YTAc6pYlMUfHTVOfz5TNU3Q/ROxKBMaVaCcLMjooUFTqGHcc8hiEIy8vLpVWwtrZWDsAJEUin6EqTe9L+/3SG5XTGo4hfpK13GsyLoF/aRZlaE+HihOWTDoCSz7//SAQOAOk4+HS+vRCDNOMuxtzHoJsYoRj7qQiERZCKQJoIFJUzHQORpuhGAC9EIB0vkeYpVCf5SH19Jfn0H4nAASNiB+GrAzu64tKx+emsvOmgnWrFT7+fBgcj6Jf2y3cKEqazJcfxdHRjOsV3mi6tFn8wSAQOKJ1chRibEF104dtXK32a5FPdhp0LqaZjDKpz9Vcn7kxN+6jsncYniMEiEciAtIsuJc3tT9N80/789LO0zz/9253EoNO4g3RgkqL6o4NEIGOqrXKVTtmk3TJMO7XgatXHA4mA6Eq3il0VAlX28UYiIC4aVfqDhfpehMgciYAQmSMRECJzJAJCZI5EQIjMkQgIkTkSASEyRyIgRObsKgJmdp+ZnTazJzt89itm5mZ2ebFvZvY7ZnbSzL5tZjf1o9BCiP3jQiyBTwO3VQ+a2bXATwHfTw6/g/ZaAzfQXl3oE70XUQjRT3YVAXf/OvB6h48+TnsBkjSR/DjwWW/zMHDEzK7al5IKIfrCnmICZnYceNHdv1X56GrgB8m+FiQVYsS56AFEZtYAfo22K7BntCCpEKPBXiyBHwGuB75lZt+jvejoY2b297nIBUnd/Zi7H7viiiv2UAwhxH5w0SLg7k+4+99z96PufpS2yX+Tu79Ee0HSny96CW4Bzrr7qf0tshBiP7mQLsLPAX8NvMXMXjCzO85z+leA54CTwP8E/v2+lFII0Td2jQm4+/t2+fxosu3AXb0XSwgxKJQxKETmSASEyByJgBCZIxEQInMkAkJkjkRAiMyRCAiRORIBITJHIiBE5li3BSYHWgizV4AV4NVhlyXhclSe3Ri1Mqk85+eH3f2c0XojIQIAZnbC3Y8NuxyByrM7o1YmlWdvyB0QInMkAkJkziiJwCeHXYAKKs/ujFqZVJ49MDIxASHEcBglS0AIMQSGLgJmdpuZ/W2xYMndQyrDtWb2NTN72syeMrMPFsd/w8xeNLPHi9c7B1im75nZE8XvniiOXWZmD5rZs8X7pQMqy1uSa/C4mS2a2YcGfX06LYTT7ZoMYiGcLuX5LTP7m+I3/9DMjhTHj5rZWnKtfne/y7Nn3H1oL6AGfBd4MzANfAu4cQjluIr2PIkAh4DvADcCvwH8pyFdm+8Bl1eO/Wfg7mL7buDeId2zl4AfHvT1Ad4O3AQ8uds1Ad4J/AlgwC3AIwMqz08Bk8X2vUl5jqbnjdJr2JbAzcBJd3/O3ZvA52kvYDJQ3P2Uuz9WbC8BzzCa6yUcBz5TbH8G+JkhlOFW4Lvu/vygf9g7L4TT7Zr0fSGcTuVx96+6+1ax+zDtGbdHmmGLwMgtVmJmR4G3Ao8Uh36pMO3uG5T5XeDAV83s0WKNBoAr/Y3Zm18CrhxgeYLbgc8l+8O6PkG3azIKz9Yv0rZGguvN7Jtm9pdm9s8GXJauDFsERgozmwf+APiQuy/SXkvxR4B/DJwC/usAi/MT7n4T7fUd7zKzt6cfetvGHGjXjplNA+8G/k9xaJjX5xyGcU26YWYfBbaA+4tDp4Dr3P2twIeB3zezw8MqX8qwReCCFyvpN2Y2RVsA7nf3LwO4+8vu3nL3bdpTqN88qPK4+4vF+2ngD4vffjlM2uL99KDKU/AO4DF3f7ko29CuT0K3azK0Z8vMfgF4F/BzhTDh7hvu/lqx/SjtWNiPDaI8uzFsEfgGcIOZXV+0MrfTXsBkoJiZAZ8CnnH3jyXHUx/y3wDnLM/ep/LMmdmh2KYdbHqS9rV5f3Ha+4E/GkR5Et5H4goM6/pU6HZNhrIQjpndRnuh3ne7+2py/AozqxXbb6a9cvdz/S7PBTHsyCTtKO53aCvjR4dUhp+gbUZ+G3i8eL0T+N/AE8XxB4CrBlSeN9PuKfkW8FRcF+BNwEPAs8CfAZcN8BrNAa8BlyTHBnp9aAvQKWCTto9/R7drQrtX4H8Uz9UTwLEBleck7VhEPEe/W5z7s8W9fBx4DPjpQT/n3V7KGBQic4btDgghhoxEQIjMkQgIkTkSASEyRyIgROZIBITIHImAEJkjERAic/4/JU5C3LYfYSkAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "gridmap = np.zeros((150,150))\n",
    "u_t      = [[0.0,0.0,0.0],[0.5,0.0,math.pi/2]]\n",
    "x_init   = [2.0,3.0,0.0]\n",
    "alpha    = [1.0,1.0,0.01,0.01]\n",
    "\n",
    "for i in range(gridmap.shape[0]):\n",
    "    for j in range(gridmap.shape[1]):\n",
    "        gridmap[i][j] = motion_model_odometry(x_init, map2world(gridmap,i, j, x_init,), u_t, alpha)\n",
    "\n",
    "plot_map(gridmap)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.3 Sample odometry motion model\n",
    "\n",
    "The `motion_model_odometry` requires high computation complexity and does not scale well to large real-world environments. \n",
    "\n",
    "One effective approach to approximate $p(x_t | u_t, x_{t-1})$ is to use **sampling**.\n",
    "\n",
    "The `sample_motion_model_odometry` function defines the sampling-based odometry motion model. \n",
    "\n",
    "This function takes as input:\n",
    "\n",
    "- the initial pose of the robot $x_{t-1} = [x,y,\\theta]$ _(**map** coordinate frame)_\n",
    "- the odometry readings $u_t = [\\overline{x}_{t-1} \\overline{x}_t]$ _(**odometry** coordinate frame)_\n",
    "- the noise parameters $\\mathbf{\\alpha} = [\\alpha_1, \\alpha_2, \\alpha_3, \\alpha_4]$\n",
    "\n",
    "The output is a new (sampled) pose predicted by the motion model.\n",
    "\n",
    "Implement the `sample_motion_model_odometry` function and verify that it is correct for some test input. **[2.0]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_0 = np.array([0.0,0.0,0.0])\n",
    "alpha = np.array([0.1,0.1,0.01,0.01])\n",
    "u = np.array([[0.0,0.0,0.0],[2.0,0.0,0.0]])\n",
    "\n",
    "x =[]\n",
    "y =[]\n",
    "for i in range(10000):\n",
    "    pose = sample_motion_model_odometry(x_0,u,alpha)\n",
    "    x.append(pose[0])\n",
    "    y.append(pose[1])\n",
    "x = np.array(x)\n",
    "y = np.array(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x1b9221281f0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnzklEQVR4nO3df5RcdZnn8ffTlQpUg0t1IIOhSAi6HNBsJK09gJs9s4BKGJHQg8iPEUUHJsc9w+ygTh/DypqAeGgnR8HZmR1lHFdcWAggtmHQiWhgZ0/GMDTTHUKAKL8SUuCQMenMaFqodD/7R91qqqvrVld3/bi3uj6vc3K66t5bVd/bnbrP/f56vubuiIhI++qIugAiIhItBQIRkTanQCAi0uYUCERE2pwCgYhIm5sXdQFm47jjjvOlS5dGXQwRkZbyxBNP/Iu7Lyzd3pKBYOnSpQwODkZdDBGRlmJmu8ttr0vTkJl9y8xeM7OnQvabmf25mT1nZk+a2buL9l1lZj8P/l1Vj/KIiEj16tVH8G3g/Ar7fxc4Jfi3BvgrADNbAKwDzgTOANaZWVedyiQiIlWoSyBw978H9lc45CLgO563DUib2SJgFfCwu+939wPAw1QOKCIiUmfNGjWUAV4uer432Ba2fQozW2Nmg2Y2uG/fvoYVVESk3bTM8FF3v93de9y9Z+HCKZ3eIiIyS80aNZQFFhc9PzHYlgXOLtn+aJPK1NYGhrJs2LyLV0ZGOSGdom/VqfR2l62Micgc16wawSbg48HoobOAg+7+KrAZOM/MuoJO4vOCbdJAA0NZrn9gB9mRURzIjoxy/QM7GBjKRl00EYlAXWoEZnY3+Tv748xsL/mRQEkAd/868APgg8BzwCHgk8G+/Wb2ReDx4K1ucvdKnc5106w74jjeeW/YvIvR3NikbaO5MdZv2hm7sopI49UlELj7FdPsd+CPQvZ9C/hWPcpRrcIdceFiWLgjBmZ14Stc7LMjoyTMGHMnk05xzmkL+e4T2bp9TvFn1XKxfmVktOz2kdEcI6O5upVVRFqDteLCND09PV7LzOKV/VvIlrkYZtIptq49d0bvVRpUihlQ7rc7088pDjSl75lKJrjl4uVVXayL36das/mdiEg8mdkT7t5Tur0lU0zUKuyOOGx7JeWaWQrCQmx2ZJQbBnZwc+/yae/wSwNN6XuO5sa48cHpm3QqBaxKZhI0RKQ1tWUgOCGdKnuBOyGdqvo9ZnN3XezObXt4cd+v+Kc9Byc1HX32vu18euMwDiTMOGKeMZobr/heBw7lOHDozSadvvu2A5ObdCoFrEJzVjlmMz0zEWk1LTOPoJ76Vp1KKpmYtM3IX0RX9m+ZdvTMwFCWvvu213y3vPX5/VMuzmPjPnHXP+bOoWmCQDm5cWf9pp0TZQ1rCpv4zArNg6W7Cu938tqHqvpdiUj8tWUfAdTW7r7ixh9NdKrGWVdnkoOHcsw8lEyWCZqbgCnNSzPpoxCRaIX1EbRtICiYTcfx0rUP1eWzW0kqmeDIZMdEE1Sxrs4k6y5cNqmf4pzTFvLIs/s0FFUkRhQIQpy89qGynboGvNh/QdnXtGMgqFUyYRw1fx4HR3Mck0piBiOHcgoSIk0UFgjaso+gWFgHcaWO467OZE2fmTBj5dsXkGijntjcmDMymsPJz1c4cCinWc0iMdH2gaBcx3EqmZhoEy9n3YXLSCZmfhE34LbLVvCVS0/nn/YcrNhJW9BRxcdcedYSMjMY8RQ3o7kxPnvv9ll3QKsDW6Q2bd80BLObrVv8muKmjnRnkt/kxqYM+TTgo2ct4ebe5dOO4il+za2Xrag4TPWU3zqKhz9zNgNDWa7bOFzlGcdbssPY8JHTJ/4Glf4+082P6Ex2cEQyoWYoEdRH0HSVLl5h/RKlSjusP/rXP2Xr82+mYlr59gXc9YfvZWAoy2fuHWa89f6UFXV1JrngXYsmpemAySOVqg2q5V4r0m40s7jJerszoRebsAltxZIJm9I8ddcfvnfKcQNDWfru3z7nggDkJ8rduW3PlO2juTE2bN5Fb3dmxrPBi18bx4SAIlFo+z6CKJTrlyg1NlbdlX3D5l3kqjx2Limk6ZjNmb8yMqpU3CJFFAgi0Nud4ZaLl5NJpzDKp3EYh4nZwZXMJj/SXFGutlCNE9Kp0FTcGzbvmniuTmhpFwoEEentzrB17bm82H/BlDQOBdXMXp5JfiTJ9xGcc9rC0Ka5wnbVGKSdKBC0uL5Vp85qKGs7MuDEriMr1iSM/Gzz6zYOT1tjEJkrFAhiIGyCWjUT13q7M2y45PSq5hsI/Py1X1fcX7j7D9POTXEyd9Vrqcrzga8BCeCb7t5fsv9W4JzgaSfwW+6eDvaNATuCfXvcfXU9ytRK1l24jL77t0/q9E0mjHUXLqvq9YWRLrNZb6Cd1KNLvdAUNzCU5cYHd07kXkqnknzo9EXKryQtqeZ5BGaWAH4GfADYS3794Svc/emQ4/8Y6Hb3Pwie/8rdj57JZ7bCPIKZqsdQxoGhLJ+9d3vZGcuZKoasyvS6OpO8c9FbJs3nCDPTOQsaziqN1rAJZWb2XmC9u68Knl8P4O63hBz/D8A6d384eK5AUEflZtoWLkhhQUIaJ2HGuPu0F/ZKfzcFA6mXRk4oywAvFz3fC5wZUoiTgJOBLUWbjzSzQeAw0O/uAyGvXQOsAViyZEntpZ6jCheNcneWg7v3h3aUVlqlTGav8DvNjozy6Y3D3De4h5d+OTrlb1NpOKsCgTRaPWoElwDnu/s1wfOPAWe6+7Vljv0ccKK7/3HRtoy7Z83sbeQDxPvc/flKn6kaweyVpqmA/J2n+haiUbooUtgxaiqSemhkjSALLC56fmKwrZzLgT8q3uDu2eDnC2b2KNANVAwEMnuF3ESlNYZa1l+W2avmNqx4HsPg7v088uw+siOjE7W4jIKE1KgegeBx4BQzO5l8ALgc+P3Sg8zsNKAL+GnRti7gkLu/bmbHASuBP6tDmaSCsDxIpSOXJF5Gc2PctW3PpDWt4c0gASgYyKzUPI/A3Q8D1wKbgWeAe919p5ndZGbFQ0EvB+7xyW1R7wAGzWw78Aj5PoKyo42ksQrzEYrnLmhqQvyEhWlNdpNaKA21hCpkNi2d33DZby+eGC/fev975jb1J0glWo9AZmW6se0zXQ9AmiPZYRx95DwtyCOTKBBIQ4SNf//wezJsfPxl9TnEhOYkCGhhGmmQSvMWek5aoNFIMVHch1CaGmP96mUKEG1ONQJpODUfxVvpGtEyd4XVCJR9VBqu3IpsyQ5T+uyYyI27Rhy1OTUNScOFNR8Vtqm2ED2l125vCgTSFGGT2ArbbhjYwV2P7ZlYra0z2cGh3Hgzi9jWTkinlP20jamPQGJLfQvRMuCjZy3h5t7locfcMLCDux97mTF3EmZccebiisdLtDRqSFpO36pTKy62o2R5jeXAndv28LfbX2X96vwiScU1hqXHpiYlMBxzn8huq2DQWlQjkFgrbq5IdyZxh4OjOSXLa7Jkh4FR1bwQA17sv6DxhZIZU41AWlJY30IxLdHZeLnx6m8YHVi69iFlR20hCgTS0opHJKlmEC/Kjto61DQkc8oNAzsmpWqW+EinkgyvOw94s8lP6yo0l3INSdsYGMpqfeaYuvKsJfSctKBic15XZ5J1FyrtRSOoj0DaRuECUnqxKSwLmU4lGRnNRVO4NnfXtj088MReRivMETlwKKempCZTigmZk3q7M9xy8XIy6RQGZNIpbr1sBS/1X8DwuvO48qwlURexLTlUNVFQC+00l2oEMmdVGnF0c+9yXtz3q0nj4CVeSjv/NfO5cerSR2Bm5wNfAxLAN929v2T/J4ANvLmo/V+4+zeDfVcBNwTbb3b3O6b7PPURSL0Uz4wtVcjhDxqVFJVKzXhaY2HmGtZZbGYJ4GfAB4C95Bezv6J47eEgEPS4+7Ulr10ADAI95GuNTwDvcfcDlT5TgUAaoZo7zqVrH4qodFJOJp1i69pzoy5Gy2hkZ/EZwHPu/kLwQfcAFwHVLEK/CnjY3fcHr30YOB+4uw7lEpmRaiavZdKpsjWDsO3SWPqd10c9OoszwMtFz/cG20p92MyeNLP7zWzxDF+Lma0xs0EzG9y3b18dii0yc+XWVkglE/StOpVMOhVRqdrb0rUPsbJ/CwND2ekPlrKaNWroQWCpu78LeBiYth+glLvf7u497t6zcOHCuhdQpBrlRiMV2qnLBQlpjuzIKJ/eOMwNAzuiLkpLqkfTUBZYXPT8RN7sFAbA3X9Z9PSbwJ8Vvfbsktc+WocyiTTMdGsrqGM5Gk5+nkLPSQvUgTxD9agRPA6cYmYnm9l84HJgU/EBZrao6Olq4Jng8WbgPDPrMrMu4Lxgm0hL6u3OsHXtuWgRzmg4aP7BLNRcI3D3w2Z2LfkLeAL4lrvvNLObgEF33wT8VzNbDRwG9gOfCF6738y+SD6YANxU6DgWaWUnqPM4MtmR0dBhwUpfUZ5yDYk0wMBQVumxYyqZMDZccnpbBgMlnRNpstJFdQ4cUn6juEinkhx1xLy2m6WspHMiTVbaqfyO//7DisnWpHlGRnMTM5azI6P03bcdaN8kd0o6J9Ikt1z8rvySjxI7uXGn777hqIsRGdUIRJqkdHhpYUGWdCrJv71+mLEZLAcp9Zcbz09Oy6RTnHPaQh55dl/bNB2pj0AkBopX7OowUEyIl7mS4E59BCIxVtqfMDCUZf2mnVpAJyYK6yO0eiAIoz4CkRjq7c5MrO8r8ZAdGZ2z+YwUCERiTIns4qXvvu1zMhgoEIjEmBLZxUtu3Fm/aWfUxag79RGIxFjYSKPCT2m+udhvo0AgEnPlsp0qhYXUk5qGRFpQYV0EicZcW/dA8whEWpjWUI6OkU97XZiA9rfbX51oNoprltOweQSqEYi0MI0qik7hFjo7Msqd2/ZM6js4cChH3/2tM8JIgUCkhWlUUXzlxrxlFslRZ7FIC9PymPH2Sov8TRQIRFpcYVSRRhLFT4fZRPNQYW2KOCaxq0sgMLPzga+RX6rym+7eX7L/M8A15Jeq3Af8gbvvDvaNAYUu+D3uvroeZRJpN6odxM+YO333bwfPT0aDfJ/C9Q/kL3lxCQY1jxoyswTwM+ADwF7y6w9f4e5PFx1zDvCYux8ys/8CnO3ulwX7fuXuR8/kMzVqSGR6pYnrjpjXweuHtTBOXCTM+MqlzV0ys5HZR88AnnP3F4IPuge4CJgIBO7+SNHx24Ar6/C5IlJB2ES0z967XbOSY2DMnU9vHGZw935u7o12Tkg9Rg1lgJeLnu8NtoW5Gvhh0fMjzWzQzLaZWW/Yi8xsTXDc4L59+2oqsEi76u3O8JVLT9dIo5hw4M5teyKfoNbU4aNmdiXQA2wo2nxSUFX5feA2M3t7ude6++3u3uPuPQsXLmxCaUXmpsKs5K7OZNRFkcBd2/ZEOuegHoEgCywuen5isG0SM3s/8Hlgtbu/Xtju7tng5wvAo0B3HcokIhX0dmcY+sJ53HbZChKmdZSj5sB1G4dZ2b8lkoBQj0DwOHCKmZ1sZvOBy4FNxQeYWTfwDfJB4LWi7V1mdkTw+DhgJUV9CyLSWL3dGcbVXxAbhRFFzQ4GNQcCdz8MXAtsBp4B7nX3nWZ2k5kVhoJuAI4G7jOzYTMrBIp3AINmth14BOgvHm0kIo13gtJUxEphWcxmqss8Anf/AfCDkm1fKHr8/pDX/QOgFIoiEepbdaomosVMs2cka2axSJsrnYhWyKop0Wl2LU2BQEQmzTm4YWAHd27bE3GJ2lvfqlOb+nnKPioikzzyrObptBsFAhGZpFUyZs5l120cZsWNP2ra6CEFAhGZRKOI4mFkNNe0oaQKBCIyiRa7iY9mDSVVZ7GITKJRRPHSjKY6BQIRmaJ4FNHAUJYbH9zJgUO5aV4ljdCMpjo1DYlIRb3dGTrn654xCgacc1rjk2wqEIjItDSSKBqFNNXdNzV2BJECgYhMSyOJonXgUGNHECkQiMi0+ladSqJD6aqj1MgRRAoEIjKt3u4MbzlC/QRRy46MNmTNAgUCEanKwVGNGoqDRqxZoEAgIlVRP0F81LuZSIFARKqiGcfxUs+RXGr0E5GqlM44lmjVs4amGoGIVK23O8PWtedqwfuIpZKJuq5ZUJdAYGbnm9kuM3vOzNaW2X+EmW0M9j9mZkuL9l0fbN9lZqvqUR4RaawrzlwcdRHaVsKMWy5ePlFDq4eaA4GZJYC/BH4XeCdwhZm9s+Swq4ED7v7vgVuBLwevfSdwObAMOB/4n8H7iUiM3dy7nJVvXxB1MdrSuHtdgwDUp0ZwBvCcu7/g7m8A9wAXlRxzEXBH8Ph+4H1mZsH2e9z9dXd/EXgueD8RiamBoSwr+7ew9fn9URelLTVi9FY9AkEGeLno+d5gW9lj3P0wcBA4tsrXAmBma8xs0MwG9+3TUnoiURgYynL9AzvUWRyhRqxn3DKdxe5+u7v3uHvPwoWNz8YnIlNt2LyL0dxY1MVoW12dybo3C0F9AkEWKO45OjHYVvYYM5sHHAP8ssrXikhMqCYQHQMueNeihrx3PQLB48ApZnaymc0n3/m7qeSYTcBVweNLgC3u7sH2y4NRRScDpwD/WIcyiUidDQxl0aDR6Djw3SeyDclAWnMgCNr8rwU2A88A97r7TjO7ycxWB4f9DXCsmT0HfAZYG7x2J3Av8DTwd8AfubvqnSIxtGHzLi1ZGbFGZSCty8xid/8B8IOSbV8oevwb4CMhr/0S8KV6lENEGkeL08RDI/4OLdNZLCLRUtK5eEh3Juv+ngoEIlIVJZ2Lh4OHclqPQESi0dud4ZaLl6M0Q9Eah7r3EygQiEjVerszpFP1b5qQmal3P4HSUItIqIGhLBs27+KVkVFOSKfoW3UqI4e0UlnU6t1fo0AgImUV0kkUZhIXlkhMdyY5oGAQmWSH1T3NhAKBiJRVLp1E/rlmE0QllezglovfVfc0EwoEIlJWWDv0aG68ySWRhBlfufT0huQZAnUWi0gIzRuIj0YGAVAgEJEQ55y2ULmFYqKRQQAUCESkjIGhLN99IqvegDahQCAiU2jdgfhoxrwNBQIRmUIJ5uJj/eplDf8MBQIRmUIdxfHQqBXJSmn4qIhMuGFgB3c/9jJjrt6BqCUTxroLG18bAAUCEQncMLCDO7ftiboYEjhq/rym1AZATUMiErj7sZejLoIUOTjavDQeNQUCM1tgZg+b2c+Dn11ljllhZj81s51m9qSZXVa079tm9qKZDQf/VtRSHhGZPTUHxUsz+2lqrRGsBX7i7qcAPwmelzoEfNzdlwHnA7eZWbpof5+7rwj+DddYHhGZpYQWGohM6YI/qWSi7onlKqk1EFwE3BE8vgPoLT3A3X/m7j8PHr8CvAYsrPFzRaTOznrblAq9NEE6leSWi5eTSacwIJNOccvFy5vWPwC1dxYf7+6vBo9/ARxf6WAzOwOYDzxftPlLZvYFghqFu78e8to1wBqAJUuW1FhskfZTbm2BwsVmYCjL1uf3R1zC9pPsMNavXkZvd6apF/5S09YIzOzHZvZUmX8XFR/n7k6F/LRmtgj438An3b2QvvB64DTgt4EFwOfCXu/ut7t7j7v3LFyoCoXITBTWFsiOjOK8ubZAYe3b9Zt2RlvANlJogMukU2z4SGOTyVVr2hqBu78/bJ+Z/bOZLXL3V4ML/Wshx/074CHg8+6+rei9C7WJ183sfwF/OqPSi8i0BoayfPbe7VM6g0dzY1y3cZjrNg5HU7A2detlK2Jx8S9Wax/BJuCq4PFVwPdLDzCz+cD3gO+4+/0l+xYFP418/8JTNZZHRIoUagIaERQPmXQqdkEAau8j6AfuNbOrgd3ApQBm1gN8yt2vCbb9DnCsmX0ieN0nghFCd5nZQvK1pWHgUzWWR6TtaXZwfDVzJNBM1BQI3P2XwPvKbB8Ergke3wncGfL6c2v5fBGZ7KN//VN1+sZUs/IGzYZSTIjMAQNDWdZv2slIE2ejysw0K2/QbCgQiLSwgaEsNz64kwOHFADiLJ2Kb20AFAhEWlahI1gLyMRbYa5AnCkQiLSYwsSwrBaPib10KjkxYSzOFAhEWohqAfGWSiaanh6iHpSGWqSFrN+0U0EgRoz8aKCocgTVi2oEIi1iYCirUUExNPSF86IuQs1UIxBpETc+qHxAcTNX1nZWjUAkhoo7hDsMxjVJODIJs9BZ2nGdKTxTqhGIxMwNAzv49MbhiVFBCgLRyaRTXHHm4rL7Vr59QUv2B5SjGoFIjGgB+fgorBJWuNgX8jclzLjizMXc3Ls84hLWj3kLJqbq6enxwcHBqIshMmPFTT6FJofCz3Qqqc7gCHV1JumcP6/swj1zhZk94e49pdtVIxBpktI5AIV258JPBYForbsw/hO/GkV9BCJNsmHzLs0BiLF2DQKgGoFIQxWvE9x6jbDto6szGXURIqVAINIgSgfRGpIJi3WK6GZQIBBpEKWDiL/MHO0UnqmaAoGZLQA2AkuBl4BL3f1AmePGgB3B0z3uvjrYfjJwD3As8ATwMXd/o5YyiTRbcfNPYbQJqPM3askOY8NHTgeYUjNr1eRwjVJrjWAt8BN37zeztcHzz5U5btTdV5TZ/mXgVne/x8y+DlwN/FWNZRJpmtLmn+zIKNc/sANTj0CkzGDDR06fdKEvDdYKAm+qaR6Bme0Cznb3V81sEfCou0+Zc21mv3L3o0u2GbAPeKu7Hzaz9wLr3X3VdJ+reQQSFyv7t2hdgJjpMPjqpSt0oS8jbB5BrcNHj3f3V4PHvwCODznuSDMbNLNtZtYbbDsWGHH3w8HzvUDoX87M1gTvMbhv374aiy1SHwoC0bKS50fNTygIzMK0TUNm9mPgrWV2fb74ibu7mYVVL05y96yZvQ3YYmY7gIMzKai73w7cDvkawUxeK1Kr4n6AdGcSdzioPoBYeKn/gqiL0PKmDQTu/v6wfWb2z2a2qKhp6LWQ98gGP18ws0eBbuC7QNrM5gW1ghOB7CzOQaShSvsBtFB88xmU7XWZK2mgo1Zr09Am4Krg8VXA90sPMLMuMzsieHwcsBJ42vOdE48Al1R6vUjUNCM4WgbcetkKUsnEpO2FpHBSu1pHDfUD95rZ1cBu4FIAM+sBPuXu1wDvAL5hZuPkA0+/uz8dvP5zwD1mdjMwBPxNjeURmZFySeBSyQ5+kxvXuJ+YOCGdmmjz18ifxlD2UWkLYWP9NfM33jTev76UfVTaVmmO/+zIKH33b+eo+fMUBCJmwEfPWkLPSQum1Mw067d5FAhkThsYypZd6CU35pr5G7F0Ksn61W+mftYFPzoKBDInFbf9SzwNrzsv6iJIQIFAWl5x+/8xqSSH3jjMG2Ot1/fVTjIa9hkrCgTS0krH+Ku5J1rJhJGbJggbaNhnzCgQSKyEje4JGzaoVM/xsuGS0/nsvdsnlt8sVegcVn9AvCgQSGyUy+TZd992MCbuMgvZPQEGd+9XDSBG0qnkxAW+3LDc0s5hiQ8FAomNcjN4c+NT7yxHc2Nct3G4SaWSo+Yn+PUblWtdyQ5j/er8Kl+a/NV6FAgkchrhE29f+r3l3PjgztAcS+XG+/d2Z3ThbyEKBBIpresbH4kOY6yoBlbanq9VvuYuBQJpiHKdvuUuGEroFg9dnUnWXbgs9G+m5p65TbmGpO7K3eUX3z0ODGVZv2mnOnpjQnf27UO5hqQuymXrLM0NU+4ufzQ3xobNuwDou2972U5gabxMOsU5py3kkWf36c5eJqhGIFWrR3t+Jp1Sp3CDJTuMo4+cx8ihnC70MolqBFKzerTnKwg0joEu/DIrCgRStVd0EY+tTDrF1rXnRl0MaVG1LlUpbUTrw8aTlmyUWtUUCMxsgZk9bGY/D352lTnmHDMbLvr3GzPrDfZ928xeLNq3opbySGP1rTpVdw4xkU4lMfI1AY34kVrV2jS0FviJu/eb2drg+eeKD3D3R4AVkA8cwHPAj4oO6XP3+2sshzTJeNQFaHOFSV439y6Puigyh9QaCC4Czg4e3wE8SkkgKHEJ8EN3P1Tj50oTKQVEPGjpRmmUWgPB8e7+avD4F8Dx0xx/OfDVkm1fMrMvAD8B1rr76+VeaGZrgDUAS5YsmX2JZUaUAqKxOgyKp1QkEwY+OdmeJnxJo007j8DMfgy8tcyuzwN3uHu66NgD7j6lnyDYtwh4EjjB3XNF234BzAduB55395umK7TmEVSvdPUuMxg5lCPdmcQdDo7mJm0/oWTCEUDrzTSJp85kB6OHx3GHhBlXnLl4YtH2atdfEKlF2DyCmiaUmdku4Gx3fzW4qD/q7mWHL5jZnwDL3H1NyP6zgT919w9N97lzPRBUm6enmvfR3Xy01JwjcdKoCWWbgKuA/uDn9yscewVwfUmhFgVBxIBe4Kkay9Pyyi3OUliIZaYXEyV0a7ywNBu68EsrqTUQ9AP3mtnVwG7gUgAz6wE+5e7XBM+XAouB/1vy+rvMbCH5wRDDwKdqLE/Lq5Snp/TiUqnmMDCUVedugxnwlUtP10VfWl5NgcDdfwm8r8z2QeCaoucvAVO+Le7e1lMhy13Iw2bvlm4PqzkM7t7PQ0++GrqIiNSH1t6VuUQpJiISdiFPdybLXsRLZ/WG1Rzu3LancYVuQ4Umn66iznV14Mpco0DQINN1+IZdyI+Y10EqmZiSy/+c0xaysn/LxPup2aexNGRT2okyBjRA4W4/OzKK8+bd/sBQduKYsCagg6M5brl4OZmgBpAwYzQ3xl3b9kx6P6mPdCo56XcNStsg7Uc1ggaopsM37K7+hHSq7BqxGss/e+lUkg+dvoi7tu2Z9HtMJROsX71MF3xpewoEDVCpw7dSuoZkwiZNKNLQz+okzBh355hUkn/9TW7yTN0Om7jYl5u8pSAgokDQEMekkmXX4z0mlaw8wcthcPf+ts7rU7yIernfgcGUu/riZpxKfTO93Rld+EXKUCCos4GhLL9+4/CU7ckOIzc2XvEuPzfubT/qp3P+vLJNY5C/6H/4PZmK6+3qYi8ycwoEM1BN6ocNm3eRG5vaoj9/Xge/fqN1m3pKk6PNlAG3XraC9Zt2lq0tFRSa1Qq/VzXliDRe2wSCai7i083UDUv9AG9esMKula0cBCAfBG67bAWDu/dz12N7mGmKqkIneOH3ubJ/S2hneYHu7kWaoy0CQTX5e6Y7Jmwk0PpNO3n9cOUmn7liw+ZdbF177pRFUZaufaji68otpdi36tSyTT9aclGk+dpiHkGl4ZzVHhM2EmhkNDdtEEglE6RTybL7LPhXGMMeZ2G/gyvPCl8fImxMfm93ZmK+hJZcFIlWW9QIqsnfM90xs53Na8CH35MfuljuDrhw8RsYynLdxuEZv38zhS1eX6gh3P3YyxOZOK84c/G0yymq6UckHtqiRhB2ASvePt0xfatOJZVMTNqXSibo6ix/p1/gwCPP7pv2Dri3O1PxvcJqDAmzactQD9M129zcu5znb/kgL/VfwPO3fFBr6oq0kLYIBGEX8eIL23THhF3I1124jOkadYpHwmxdey4v9l/A1rXnTrkbXnfhsrJluO2yFYyH9M6Ou4e+7sqzlkyUd6ZSyQ4124i0ibZoGqpmKGK1x5S7GA7u3j8lfUGxsNrGTMoZNsGqOCVFpbKHjdIJ85vcOFvXtnWWcJG2UdNSlVGJ41KVA0NZbnxw55QU0vXKYllu2cmZvHfY64+Y11F2XH8mnVIgEJljGrVUpQQKtYV6rTdc7v1h9hOswl4P5WfwahinSPuodfH6jwDrgXcAZwQrk5U77nzga0AC+Ka79wfbTwbuAY4FngA+5u5vTPe5cawRtLJGBS8RiZdG1QieAi4GvlHhgxPAXwIfAPYCj5vZJnd/GvgycKu732NmXweuBv6qxjLJDGkYp0h7q2nUkLs/4+67pjnsDOA5d38huNu/B7jIzAw4F7g/OO4OoLeW8oiIyMw1Y/hoBni56PneYNuxwIi7Hy7ZLiIiTTRt05CZ/Rh4a5ldn3f379e/SKHlWAOsAViyJDylgYiIzMy0gcDd31/jZ2SBxUXPTwy2/RJIm9m8oFZQ2B5WjtuB2yHfWVxjmUREJNCMpqHHgVPM7GQzmw9cDmzy/HClR4BLguOuAppWwxARkbxah4/+HvA/gIXACDDs7qvM7ATyw0Q/GBz3QeA28sNHv+XuXwq2v4185/ECYAi40t1fr+Jz9wG7Z13w6hwH/EuDP6MZdB7xMRfOAebGecyFc4CZn8dJ7r6wdGNLzixuBjMbLDfettXoPOJjLpwDzI3zmAvnAPU7j7ZIOiciIuEUCERE2pwCQbjboy5Aneg84mMunAPMjfOYC+cAdToP9RGIiLQ51QhERNqcAoGISJtru0BgZt8ys9fM7KmQ/X1mNhz8e8rMxsxsQbDvJTPbEeyLNA92FedxjJk9aGbbzWynmX2yaN9VZvbz4N9VzSt12XLWch5jRX+rTc0r9ZQyTncOXWb2PTN70sz+0cz+Q9G+881sl5k9Z2Zrm1fqsuWs5Txi8d0ws8Vm9oiZPR38f/mTMseYmf158Dt/0szeXbQvFt+NOpzHzL4b7t5W/4DfAd4NPFXFsRcCW4qevwQcF/U5VHMewH8Dvhw8XgjsB+aTn7z3QvCzK3jc1WrnETz/VdR/hyrPYQOwLnh8GvCT4HECeB54W/C32Q68s9XOI3gei+8GsAh4d/D4LcDPSn+nwAeBHwIGnAU8FmyPzXejlvMI9s3ou9F2NQJ3/3vyF5NqXAHc3cDizFoV5+HAW4J030cHxx4GVgEPu/t+dz8APAyc3+jyhhZy9ucRG1WcwzuBLcGxzwJLzex4QlK0N7q8YWo4j9hw91fd/Z+Cx/8GPMPUrMYXAd/xvG3kc54tIkbfjRrPY8baLhBUy8w6yf8n+G7RZgd+ZGZPBNlQ4+wvyK8c9wqwA/gTdx8nPC14XIWdB8CRZjZoZtvMrDeqAlZhO/kFnDCzM4CTyCdZbLW/Rdh5QAy/G2a2FOgGHivZFfZ7j+XfYxbnATP8bmjN4nAXAlvdvfgO6T+5e9bMfgt42MyeDe6i4mgVMEx+8Z+3ky/v/4u0RLNT9jzc/V/J503JWj5n1RYz2+Huz0dY1jD9wNfMbJh8MBsCxiq+Ip4qnUesvhtmdjT5m7jrgv8rLamG85jRd0M1gnCXU9Is5O7Z4OdrwPfIV+3j6pPAA0G18TngRfLtumFpweMq7DyK/x4vAI+Sv2uKHXf/V3f/pLuvAD5Ovq/jBVrsb1HhPGL13TCzJPmL513u/kCZQ8J+77H6e9RwHjP+bigQlGFmxwD/maK02GZ2lJm9pfAYOI/8ms1xtQd4H0DQjnsq+S/tZuC8YARIF/nz2BxZKadX9jyC8h8RbD8OWAk8HVkpKzCztOVTsANcA/x9cHdXNkV7VOWcTth5xOm7EfQl/Q3wjLt/NeSwTcDHg1E3ZwEH3f1VYvTdqOU8ZvPdaLumITO7GzgbOM7M9gLrgCSAu389OOz3gB+5+6+LXno88L3834d5wP9x979rVrlLVXEeXwS+bWY7yI8q+Jy7/0vw2i+SvwgB3FTS/NVUsz0PM/uPwDfMbJz8DU2/u0cSCKo4h3cAd5iZAzuBq4N9h83sWvIXm0KK9p3NP4O82Z4H8fpurAQ+BuwImrAgP/JsCUycxw/Ij7h5DjhEvtaJu++P0Xdj1udB/u80o++GUkyIiLQ5NQ2JiLQ5BQIRkTanQCAi0uYUCERE2pwCgYhIm1MgEBFpcwoEIiJt7v8DlQjhWrZiSaAAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.4 Evaluate sample odometry motion model\n",
    "\n",
    "Assume that a robot starts at pose $x_0 = [2.0, 3.0, 0.0]$ in the map frame and moves according to a motion model with $\\mathbf{\\alpha} = [0.1, 0.1, 0.01, 0.01]$.\n",
    "\n",
    "The robot obtains the following odometry readings:\n",
    "\n",
    "1. $\\overline{x}_0 = [0.0 , 0.0 , 0.0   ]$\n",
    "2. $\\overline{x}_1 = [0.5 , 0.0 , \\pi/2 ]$\n",
    "3. $\\overline{x}_2 = [0.5 , 0.5 , 0.0   ]$\n",
    "4. $\\overline{x}_3 = [1.0 , 0.5 , 0.0   ]$\n",
    "5. $\\overline{x}_4 = [1.0 , 1.5 , \\pi/2 ]$\n",
    "6. $\\overline{x}_5 = [1.0 , 2.5 , \\pi/2 ]$\n",
    "\n",
    "Evaluate the `sample_motion_model_odometry` by considering 1000 samples and plot the resulting positions for each sample in one unique plot. **[3.0]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "commands = [[0.0,0.0,0.0],[0.5,0.0,math.pi/2],\n",
    "            [0.5,0.5,0.0],[1.0,0.5,0.0],\n",
    "            [1.0,1.5,math.pi/2],[1.0,2.5,math.pi/2]]\n",
    "alpha = np.array([0.1,0.1,0.01,0.01])\n",
    "\n",
    "trajectory = np.array([0.0, 0.0, 0.0])\n",
    "\n",
    "x_0 = np.array([2.0,3.0,0.0])\n",
    "\n",
    "samples =np.array([x_0]*1000)\n",
    "\n",
    "for i in range(len(commands)-1):\n",
    "    init = x_0 + np.array(commands[i])\n",
    "    trajectory = np.vstack((trajectory,init))\n",
    "    \n",
    "    u = [commands[i],commands[i+1]]\n",
    "    itr = (i)*1000\n",
    "    \n",
    "    for j in range(1000):\n",
    "        pose = sample_motion_model_odometry(samples[itr+j],u,alpha)\n",
    "        samples = np.vstack((samples,np.array(pose)))\n",
    "trajectory = np.vstack((trajectory, x_0 + np.array(commands[-1])))  #append the last commmand to trajectory for visualization\n",
    "# print(trajectory)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1b922178fd0>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxJklEQVR4nO2de5Ac9XXvv2dGLWlWtjWLEb5iQEhgR9hE1m60F4iVopBMkIMUecNLJCY3VOVasZ04loB1rXwVLBEM65JtlFynnODkVtkRNoslmCBkW9iWKGLZgqzYXWSBlBjrAQP3IqFdHtqRNDt77h8zPert6edM9/TrfKqmdne6p/s3sz3fPr/zOw9iZgiCIAjRJxX0AARBEARvEEEXBEGICSLogiAIMUEEXRAEISaIoAuCIMSEKUGd+Pzzz+e5c+cGdXpBEIRIsm/fvhPMPMtoW2CCPnfuXAwMDAR1ekEQhEhCREfNtonLRRAEISaIoAuCIMQEEXRBEISYIIIuCIIQE0TQBUEQYkJgUS6CoJIfLGDTzkN4bbSIC7MZ9Cybj+7OXNDDEoTIIYIuBEp+sIB1j+1HsVQGABRGi1j32H4A8FTU5aYhJAERdKElmAnqpp2HamKuUiyVsWnnIc8E181Nw0vhl5uI0GooqHroXV1dLIlFyUAvqCrtbQpGxkqGryEAh/uWmx5PFcqZGQVEwOhYyVQ0F/ftQmG0WHecXDaDPb1LLceZUdJ44MYFroXYy2MJghYi2sfMXUbbZFFU8B0jKxyAqZgDwIXZjOHzqlAWRotgAKPFEkbGSmCcs7zzg4VJrzESc+3z+cECFvftwpr+IdPZglusZh6C4BfichF85zUTQbVi7Ow48oOFOmvW7OagonfX6MVdS5rIdPagRT9+J64Us/dcGC1icd8uccMIviAWuuA7Zta2FSNjJUNr28nN4TWN5a36yo0oM9veIIDJ488PFnBn/1BthlAYLeLO/qG6cZq9Z6q+xmpGIQiNIha60DCqpVoYLSJNhDIzcgZW55LLZ2HL3mOuj18slbFx+4FJ1nDWwu+ukm1TAAAbtx+wFWszd4xKRkmjZ9n82t/rHnsBE7p9JqrPa99zz7L5dZY/AdCvWHm9ACwkGxF0oSHW5/fj4b3HagJVri6uG0WQ7D54vOHzjIyVagJeGC1CSRGUNKFUNl/MHxkr4SN/8yOMlfTS6440Ud0iZtHkmPrn1deoN6PpSsr0tdpZR36wgA1PHMBosfKe29sUfPkPrxDBFxwhgh5zvAid0x4j26bgdKlsKk5AvdVpZwW7oTTByCgpS0EH0LSYA8AEM7o7c5NmIlaovvFsmwJm4K1iJfLmY5edhz0vnzR9neqeyQ8W0PODYZQmzr23kbES1vQPYeDoSdzXvaDp9yTEGwlbjDFehM45WTQ0I5fN4LWqvziqtCkplCbY9gbSKEqKsOmWhejuzJmGVwIVd82DqzrEUhcswxbFQo8Remt87Oy4Zeic3nI3es7JoqEZXlrmQeGFpW9FaYJr/w+rBV8GxNcu2CIWekxoxpIGKpYiCJMs0YySbvh4gjuMFkyN9jFLthKSg1joCaAZSxrAJL+tSrFUrkWvCP7i5BNuJPxTSBYShx4TGknecUKZGRkl7cuxBXcsudywL7Ag1BALPSZcmM1Y+qydTOnNmDYlhelKyjb+W/AXs/BPp/kAQvwRCz0m9CybDzLZ1t6mNDVdHy2WcLo0gfZqwo4QDEazMG1tG6A+H0CyUJOFCHpM6O7M4VNXzzEU9ZGxkqOIk2xGQTZjLNrFUhmnZYE0ULK6G2p+sIC7Hh02XTuRYmDJw5HLhYiOAHgHQBnAuH6FlYiuBfBvAA5Xn3qMme/1bJSCI+7rXoCuS86rTb/duFkyShobVl6BTTsP1bIU9VglE1mRThHKBouugjuYJ7tXnPx//VpbEcKJGx/6EmY+YbH935l5RbMDErzD7suuJv5oM0jX9A95Pg4Rc28YLZYmhaZKZIygRxZFY4RaCdCNHa3PPpQwxXDjJjRVX1hMiD+OEouI6DCAEVSMgn9i5od0268FsA3AqwBeA3A3Mx8wOM5qAKsBYM6cOYuOHj3a5PAFLR/+mx+5douo03YR8niRzSjYsFKKesURLxKLfo+ZC0R0AYCfENFBZn5Gs/15AJcw87tEdAOAPIAP6Q9SvRE8BFQyRd28iThgVSir0W1aGvFx66slCuZML50GAJxWpgc8EnveOT0e9BCEAHCd+k9EGwC8y8xfs9jnCIAuK5970lL/rQplAXC1DTAuqzq3d4efbyHxPPK9XgDAbX/SF/BInJECMF1J1erRSCneeNCUhU5EMwCkmPmd6u/XA7hXt89/A/D/mJmJ6EpUrqU3mx96fLDrMelmG3Cuow9Q37neCHGpJI8JTC4uNjJWQs/WYQDOrhkhejiJQ/8AgJ8T0TCA5wDsYOYfE9FniOgz1X1uBvCr6j5/D+A2DqrqV0ix6jFpFiP+2mjRMuzMTZyxiLkAVIqvSWx6fLG10Jn5NwAWGjz/j5rfvwngm94OLV7YpeabvQawLkNbGC2i896nIHotOEXfIanZBihCeJCwxRZh1GPSCm3Imd3rpMaK4IYUEeb17kC2TcG7p8drlTYLo0Ws7R/Cmv4hqQUTUUTQW4S+x6SVQW30ZdL2mRQap5kiZVsXXOflUAJDdb8ZGQLqZ2PUG1YIP9LgIiA6Nj5lKNDZjIKhL18PoH46PPf9Gez9zYj4w4WWkctmsKd3adDDEDRIg4sQQialEdXn9WGOVounQmtoH3sLADDSNjPgkbSOwmhxkvFhFvoovvhwIIIeEKMmfu+RsRLm9e5ASsIMQ8e38g8AiE4culdoZ5JGoY9Gxoe4a4JByucGhFXRJIaEGXqBkgKUtFmVeKFR9KGPdjkWQusQQfeZ/GABi/t2YV7vDizu21VrONCzbH5sWru1KY1fRhklhc2rOnCkwebHZnKdzSj4r/uXY9PNdRG3ggdoQx/NciWkdG/rEZeLjziZijqJegk77TOm4f5l8xsqvXvejGm137MZxXUkz8yMgjPjE5MsRNUqn9e7AxdmM2hvUyS002O0M0yzHAsp3dt6xEL3EbupaHdnDnt6l+Jw33LkInzxvzZaRHdnrqH3oG2VtmHlFa5f/1axhAduXIBcNgNCZdEOXPH7cvX470qhKk9R0lTLkcgPFjB2tv7zldK9wSCC7iNupqJGLpioeH9VS6xRN5J6k2tkAe3CbGbSjbFt6pRaooyK/u9G2dJ5A7Z03uDJsaLK1DRhxtQpWNs/hI6NT6Fn63Dd7CebUfDAjQtkQTQARNB9xGzKafR8d2dukqWZy2Yi4YbRWmLdnTnctCiHdDX2MkWVhUknqDc5Myt/xtR03c3CyAp067e1u2lqtz/54Wvw5IevcXV8N5j1cw0TZ8tcm/2MFksoleuv0tFiCZt2HpIG1QEggu4jRharXoS0i6abdh5Cz7L5ONy3HHt6l7bUDdPeZt4g2opiqYyN2w8gP1hAfrCAbfsKtQidCQampNO4/eo5tRtV2iQA38rKJwCnzpYxXUkhm1FqNzwjK9DsJprNKIb/i09Vx2YG49xN5sK3j2P228dN922WM+MTkRB1J2hdaULrkExRn9E29VVL2Kqp/YB5HfTuzhzW5/djy95jLRlnRknjpkU5bNtXcNXmTEWpTsWNFjW12YZWdeH1DT2MGiHr99VjV3feLPllcd8uw4U97dhPdP0ufv3GKdz2Jw84/2BcErcF3PY2BW1Tp0jCkYdIpmiAqBdvz9bh2vS0MFo0jQjRLppu29c666ZYKmP3weN44MYFuOvRYddx8KXqVNwIrRtEH91j9CXv7syhuzNnKLJ2/na745u9zqh4mn42df57ptWeb+Sm54TRsRLaNE0pos7IWKl2g5KEI/8RQW8BG7cfMPQ1mvHaaNEwQsZv1GgVwLjCIxEaKtOrd4Oogu1kPG6ed3t8/WsA6xsNUBH1Rm96TlDP2/ODYc8Wc8NEMwvggj0i6C3A7RQ626YEUrdFFV4rcTNzTZhBQMPha62Ob3Z6IzC76anuofY2BcyVkMqMC2tbnRGox7/z0SGYaXo2o+DU2XFXhkJYkIQj/xBB9xAvChQpaQokblrvXjATN7fWI6Px6bUTN0hQOLXoAfM1gTYlhWlKGqNjJcPXp4kwoZsFKGnCppsXYsMT7mZ9YWJmTBZ+w4gIukeYZYUOHD3p+Bi5bAanzowHUvfcadywuo/T+uzNROq4Ec2WcNddk/5s1qJnkGnT5k07DxneNGdMrXxlo1wb/+3TJeQHC+J28QGJcvEIM1eEG79zkI2cG62lYuWCsYtISRJOomi0zOvdYZiHQGisnWHYUCO9QnOzjhBWUS4Sh17FrIiWU8z8gm702amYq3HYt189x3FmZi6bMbWWm7GizbJD29timC146FDl0QBuF3itktKsfNCpiKQXqzPYQrWOkcSte4MIOs65S5q5uFpZiEgd42P7Xp2UmWmG6nd2m+jk5MZmlOG6eVUHBu+5Pl5iDgB/8ReVRwO4yRoGrJPSzF7T3qbgG7d2RKZkhJTc9R7xocO6iJZTUXLbBNoLxkoT2LL3GNqUFEoTmLRIpi6+GfUnNZvmNtqooJEwwaThdoHXbv3A6FiqP37g6MmWJaR5jUTANIf40GHtrzzswresj3Jp5QKnkiK8Z/oU04gJJ7j18yaOa6+t/Hz66YZe7mWbNrtjze3d0dBxg0aN7BGfujmSKWqDV/HOekvVKA1dj5Mu9CmCaTyySmmC0TZ1Cgbvud75gHVIowJ/8XImY3esXEQXTtV1JMkqbYzE+dCNfMROfMtujwlM9i8D9ZX9MkoaH7vsPFsfuNOEwWaF162fVwgvPcvmR779nvjU3ZMoC93MR/zAjQvwwI0LXE+H84MFbNx+YFImqN6y0FpS+mnykstnTapO2CzNCm+YE3lCwfr1QY/AMeo1p78+o4bMDt3hyIdOREcAvAOgDGBc778hIgLwdwBuADAG4A5mft7qmEH40L30Edu5U5wc000afa56A3j42WOGoZBexXx76ecVwoORMfHk8OuhT1CS9Zt6vPKhL2HmEybb/gDAh6qPqwB8q/ozVHjpI7YrnuXkmE7PmyLUhPW+7koZWL+EVyJWLBgaqvzs6AhyFA1h9H/dffB4qAVdZofu8crl8kkA3+WKub+XiLJENJuZX/fo+J7gZbEnO8s6RVRrUqwPDVSFOOUwMzRtkC0iwhsAa9ZUfjYY5RI2wuzOMAq3FexxKugM4CkiYgD/xMwP6bbnALyi+fvV6nOTBJ2IVgNYDQBz5sxpaMDN4MZHbGUB5wcLttEpRqv1wOT4Yae+81KZsXH7AQAhqmsiRB4zA0cNHcxqqkZmW9R4Y8bUNA7c+wnfzxNXnAr67zFzgYguAPATIjrIzM+4PVn1RvAQUPGhu319szgt9mSXYLNp5yFTMTeq3aJdrTdy0zip4TIyVqprkiFhXUIzmBk4ancn7baRsRLSKULZ5xrtY2db2wMgbjgSdGYuVH++QUSPA7gSgFbQCwAu1vx9UfW50OHEVWGXOWo5VTW53q1eM1FtS2fnxtGXS5VmAUIz2NW9138HyhOMGVPTGDtb9q2BebZNSus2g20cOhHNIKL3qr8DuB7Ar3S7PQHgf1CFqwG8FTb/uRvsFk/NfO4E8wvywmzGMs7brMhVo2MVBCd0d+awp3dprTG5KvJm19XY2TIO9y33rZn1yFgJnfc+5biOUrNF9eKGk8SiDwD4ORENA3gOwA5m/jERfYaIPlPd54cAfgPg1wC+DeBzvoy2Rdgl2PQsm29YAIlRcbeYJSkZibaSJpw6M461/UOYNiWFdpcWiiT9tJD77688EoDdd2DDyit8O/fIWAlr+4ewPr/fcj8viurFDVtBZ+bfMPPC6uMKZv5K9fl/ZOZ/rP7OzPyXzHwZMy9g5nAUaWkQu8zR7s6c6ZTzrWKprvqgGh+ur0zY3qYAXGlWwKj8PF2aMCyLq6SoLvNPwrpazMc+VnkkACffAb+sdKBiHD2895ilOFu5RpNK7DNFG4nXNsqymzZl8r3PzOd9YTZj6afXblvct6sucqBYKmP3weOGmauARLkEyi9+UfmZAFF3EkCwYeUVWNM/5Oq4bpq4MCrfQaMx5AcLpmtOSXZDxlrQGy0Hq3Ja09x3tFia9Fov0uStfPVmNwUR8AD50pcqP2MSh26HXQCB21K9BODrty6s9Vd1wshYqWb0aNs6bttnbrkn2Q0ZO0G3S9xxGhliF+niRb9Lq0QnScEXosB93QvQdcl5NZG2ssCnpIC1/UNNRcgUS2V8/9lXTM+RdDdkrOqhOylXq6L2ZjQTSq9qpFthNN6MksZNi3LYtq9gGB8soh4gTdZDTwpuahR5zeZVHbH/jiSmp6hdfRUtdqvirSgla9S+7YEbF2D3weOy2CNEFi982FYlpc225arrV0kmVi6XRi4kMxdMq0rJGvkp15osNCV5scct4rIKDjNXolOyGQVvWRQN++OrLjacwSbZ1aISKwu9UevZSCjNrOdWiII0mmgO3+KTN2+uPARLzPI0nKKG8Zqx44XXUSyVa5Z6K7+bYSdWgm6VbZlR0qZxs2ZCaZZF5zfNdlBKOr7FJ3d0RLJ0bquxytPwAjXqpcxc+16ImFeIlctF/afe9eiwYXTLdCWFjJIO/VTNiwiaOGHkPhk4erIW7ZAmwh9fdXGtVrxvvVF/+tPKz+uua+44CaBVPU2lntFkYhXlomIVofLgqg4RyghhFAmUAjBhsn8um8GpM+OGjRusut848rlLlItj3ESceYFd1Fqc8KpjUWSwiu+WxhDhRi+sp86M14mCmZgDFX+5kqY60VdSZDoTayQBTRZdrdHPMtumpnHKx9K42rUS7fmTRiwFXZodR5P1+f2Tsg4bnbLrywwDgHaVzslNw2gqf+LdMzh2sogbe3dManAiQmKM3nhan98/yU3mVXN0LUl3wcRS0MUHHQzNWK35wYLjFPJGKJW5tiiqbxRiRmG0iPxgoVY75MLjpzBRFSG9FCVdSJxwX/eC2joH4F8CUpLDe2Mp6MBk60AVmrX9Q5ZdiuQG4A7tZzYzo+DU2fGGOyq1ImmqMFp0XUxqbf8QBo6exO6Dx/F1G4tSFRK5lpxhNpOeNiXVVPNqbdRa0v4XsRV0FSf+0WaLeMUZsy+E/jMz+gK6sVq9tKrs+r26QS3jygC+tOyvLPdVa/DIteQMs5m0WWKdE7Su1ST+L2IZ5aLFbFqnjXhwsk8Ssao1Y1UgSc8RTe0bsxuEV9Pv9jYFH5n9Xux5+WTTx9LixOe7uRpBJddSc3RsfKohC10NX9UWCzOivU3B4D3XNzvMwEhMLRcjnMQk+xa3HHHMEnS27D3mWMwJqGVo5gcL6Nk6PCmDs2frMPKDBSy5fJbtsdoU+8v1dKmMwWOjjsbmhjIzPv7rZ/HxXz9ruL29TUF3Z85URAqjRWmR5hCLMi6WlJmxZe8xrOkfsjQORsZKsf0/xF7QnaTRS6q9MV7c0Bjn/OMbtx+oi0AplRn/6/H96H/uFdtjjTvoOF8sTWCsZBXY2Diffu5xfPq5xw23jYyVcNm6H1q+vjBaRM8Phifd4KQfZj2jY86sc6sCXnbEtdBd7AXdSRq9pNob49UNTY0W0XdnUjl1toySA7E+axSOGCKczFpKE4wNTxyQfpgW2F13GSWNzas6mgp7jOvsO/aC7qTIVpCFuMKMVW0cI6wMJnUxSqgsIEs/THOMrjv10splM7hpUQ4btx9o6hxxnX3HPsoFsG+l5XSfpKGNQnCyYGllMLUqBTwqyLqNOVZ5JF6UFIjz7DsRgi40jvrlamVdjrgzbUoK579nmml5CsHcwHLTxEaLGsqaJpo0E4qbERcpQU9akkBYaPRLFDfWrrjLk+OcGZ9A29T6yp8AMHZ2vJadKtTT6AyGURF11e8e15j0yPjQZRGptWgjMILqDxk2Xn/fLLz+PvvwSif81xuncNOiXF2N/pGxklzXFjQzgzEr1xAnIiPosojUOvQ3T6HCipeewYqXnvHseE8Ov44Z0+onycVSGWv6h2qhjBLeeA63C/V2xG3NIjIuF1lEah1OXCxeptdHhdsHK3HmT374Gk+ON1osWWZEqnHrIDRcIydu6BdMZ2YUEFVi162uR7PrNW5rFo4FnYjSAAYAFJh5hW7bHQA2AVBNh28y8z97NUjAusa54C1WN0m1kYC4YVqDUXx+0is7mi2YmpWPUEsCJKGxtBuXyxcAvGSxvZ+ZO6oPT8UckOSfRmlkum52k8xlM3hwVYfHIxQaQWam9ZhpxNdvXYj7uhckItfEkYVORBcBWA7gKwDu9HVEJkiNc/c0Wm3OrKzpkstnSfhiSJCZaT12GpGEXBOnLpfNAL4I4L0W+9xERNcA+E8Aa5m5rjgHEa0GsBoA5syZ426kaPwfktRwR6uFZKv3b/bFkPDFcKCkzdvpJR0nGhFnPbAVdCJaAeANZt5HRNea7LYdwPeZ+QwR/QWA7wCoqxXKzA8BeAiolM9tdNBuSGJNZBWnC8lmF7i+QUjS/eaf7V4X9BAqJG012kPirgdOfOiLAawkoiMAHgGwlIi2aHdg5jeZ+Uz1z38GsMjTUTZBksMdnVSRtIvv125POiNtMzHSNjPoYdQKfAnu2bj9gKEeaMNEo4ythc7M6wCsA4CqhX43M9+u3YeIZjPz69U/V8J68bSlJC3cUd8WTknTpJK1+oVkO7eMuFnOcfP+nwIAti64LuCRVEIeJaPUnEZmlWqY6MbtBzA6VoqkO6bhOHQiuhfAADM/AeCviWglgHEAJwHc4c3wmidJ4Y5GbeGUFKG9TTG9QM1ubGrJW6svBJF1Qa64ESZBB2C7FhJnX7EVzRTwKk1wrcxzFN0xrjJFmflpNQadme+pijmYeR0zX8HMC5l5CTMf9GOwjZCUcMf8YAF3PTpcdxGXJhhtU6fUwg3X6qaWVje2nq3Dptty2QwevLUDSrrxJgNCc1jNMpNcKsPLWWXU3LORSf1vlCTUOle/vGYF/9Uvs9GXu2fZfFNR1ncX0rLk8lno7sxh1X+/2Iu3IDTATF0dGC1JXjvy2p0aJfds7AUdqIj6nt6lONy3HHt6l8ZKzAF7i0QtGapF6yefMdW95+3J4cqSye6Dx12/VvCGt0+b98ZM2tqRFq/dqVFyzyZC0OOO1Zc0o6RNLXf1dW810GFdrUEi0S/BMcEwjc5Icp/cRt2pGSVVN1uNmntWBD0GmH1J00Q1d5PV67Jt5lN3K/KDBSTJg37HLRtwxy0bgh5GHYXRInq2Dk8S9SWX15f5TUpCUndnfVliFQJw+9Vz6oRbSRMeuPGj2HTzwki7Z4kDClPo6urigYGBQM4dN4xW9TNKunYxWm0HgJ4fDDtq0qxlxtQ0TpcmmmrUK/hDe5uCt4ol6P+lKQDfWNURKYFqFKtr3iycMZfNYE9vXT5k6CCifczcZbRNLPQYYLfwa7V9085DrsUcAM6OJ0/Mb39+B25/fkfQw7BlZKxezAFgAkjEoihgfc1bheoCjRW0CwuRqYcumJMfLGDDEwdqfu2xs+N1+5jVuLCMM4dxljnBuKyrH6QIhuIUBCsO/jsAYMvvLA94JI2TpDUPs2veLDeFAKzP759UZjdqsehioUec/GABPT8YntQoYWSsVOdTNXutmQ88TYRPXT3HMIbfSl+VNEFJeeNZzyhpfOPWDhzpW47Nqzpq1lZ7gz5/oSJaUbI4/aBn2XzD654BfP/ZVyId7imCHnHMXCalMttehJt2HjK1wK1qSJsJKgGYMXUKShOMNFW+Mtlq+QH9fnYYuY3U0NO2BsIshQqM5LhdzOjuzJkaJXYRYWFHvhkRQ5/ObTWFtrsIzbYzMElI9VNNs8JQjHPhjGVmEIAVC2ej65LzHI9Zu5jrdtxWqO6jJLbO01MYLWJx367ElAIwImdyHaaJDEU9KuGeYqFHCKN0bitr1+4iNMs0NAtzVHEat84AHt57DAAmJXaZhZSpYZZ2ImP1vtSZQXubgmxGqc0sHqy6bJoV8ynpFNqU6H9tklQKwIieZfPrXINKqtKqLsqlQsRCjxBGGaFmAmUXc5wfLOCUweKpkrKPVXbTU1Sd4mtrq5udd9MtCx1ZjGYdlexuBmv7h0y3OVl8vetzf4c9vUvxIoBPffuX2PPySduxhpmk9yats4YI6LrkvLoZZZRmMiLoEcLK1ZDNKDV3R3ubgi//4RWWF+GmnYcMa7W8Z/oU24vXSFCtXBnacTdzXpVG2xGa3Yi08ced9z5Vq7anR3ujO/JmNHyqdkTFN+w1Rtehuu4U5fIgIugRwokgAef87Gv7h0zFzuyLPGoiZlqMBHXJ5bPw8N5jhqKudZE0c179GNx+6cwse61QW42j+2ffB34G4O67YyOEUfENe01ca92IoEcIJ4LktMVWs3XizQRVL+r68QVZn96JZW9108STT1b+uPtuV26nsBIl37DXxLVPQvRXdxKEk1LATsum+lEn/r7uBbXFR7PxBV2f3q7yptPxGe0XJdrblEn/myhnRzZC0NehX0gtl5gxr3eHaWz54b7JGY5BdbQJeycd0/Fde21lh6efNt1v4OhJU9dTmEgTYYK55i7TZkcCzhaZo47+/7fk8lnYffB4aK9LFataLiLoMWNx365IFx4KNTpBN0Pf15UItRaAb7xdRGnC95G6wmxBO0nXjF2BuzBhJejiQ48ZTvzsQoNkmltfyA8Wqm39wmW/O4lOijt2zdKjggh6zDCLQLGLehEc8KMfNfVys5DNsBL1BUI3xCXqRQQ9hmgtRKdRL4L/REkckjari0vUi0S5xJwkNwv2nL/928qjQazEgULU+okA3LTIfZx/lIlL1IsIesyJy1QyFPzsZ5VHg5iJxuZVHXjw1g7Pyg43CyMZzb+1oZqbdh7CTYtykW4/B4jLJfbEZSoZB5wkNqnbgGCXTuN+wzdyRW7bV4ikiGsRQY85EvUSLqxKFmi3zesNttVd3G/4cYlq0SOCHnMaLWQlBEuQpQWScMO3c0Xq2zo6KXgXBhwLOhGlAQwAKDDzCt22aQC+C2ARgDcBrGLmIx6OU2iCRgpZCQa8//0tO5XRzKpVJGFB1MoVqbZ11HYCU9s6AuGODnOzKPoFAC+ZbPtzACPM/EEADwL4arMDE4TQsW1b5dECjOr2LL7svJacOwkLolZRLVZtHTduN+7WFRYcWehEdBGA5QC+AuBOg10+CWBD9fetAL5JRMRB1RUQhBhgNLNan9+PLdUuUH4R9wVRwNoVadUIZWSshPxgIbRWulOXy2YAXwTwXpPtOQCvAAAzjxPRWwDeD+CEdiciWg1gNQDMmTOngeEKQoCsW1f5+cADgQ3hvu4F2H3wuK/+dXVBNOxF1JrFzBVpt34R5oVTW5cLEa0A8AYz72v2ZMz8EDN3MXPXrFmzmj2cILSWX/6y8ggYM3eBFy4Z1e1g1L82KT1IjfqNagnzDMaJhb4YwEoiugHAdADvI6ItzHy7Zp8CgIsBvEpEUwDMRGVxVBAEjzGr17NtX3Nim9PU/TGyUOMQ1ucE9f2t7R+y7cAVNmwFnZnXAVgHAER0LYC7dWIOAE8A+DMAvwRwM4Bd4j8XBP/QuwsW9+1qOiJm7vszdXXR9YTZOvUS9bONWg5Hw3HoRHQvgAFmfgLAvwD4VyL6NYCTAG7zaHyCIDjAC6Hd8/JJ233CbJ16TRRzOFwJOjM/DeDp6u/3aJ4/DeAWLwcmCE5o6cLdRRf5c1wPaEUiEgGhtk79IGo5HFKcS4gsLV+427Kl8gghZgulXsIId1KNIIIuRJT8YAF3PTospYGrmDUQ95JcgtwtUUVquQiRQ7XMyybr7r65HtasqfzcvNmf4zeJkXvALGLFLRkljSWXz8Livl2R8ScnERF0IXJs3H7ANqJjXu8OZNsUMANvFUveCNDQUOOvDQgvasLkNGGR0vlqMmFLvhKXixAp8oMFjIyVbPdjVNK0R4ulxCXGaNG6YhqBAOzpXYrdB4+Le0tHGJOvxEIXIkUzAlIslWvFlYysqvxgARu3H6jdMLIZBRtWhr9kqh2qK0bf1MEJapiidL6qJ4w11UXQhUjRrICoZVBL5Yr/XbWqBo6eRP9/vFJ7HgBGiyX0/CD8JVOdoo2rLowWkSYyXYcAACVFtTBF6XxVTxhvciLoQqTwIt5aK9pAxaoyq2BYmuBzFtdv/VZT5w0D+oVTS6udgIGjJ2s3AMLktnhhz5r0mzDe5MSHLkSKnmXz0epWyoXRIub17sDiS29D/rNfbvHZ/UX1saep/lMtlRkP7z1WEy0Gap99VJsoe4lVTfWgEEEXIkV3Zy6Q5slhWfTyg+7OHCZMXC/6ZxkVMd/TuzTRYg6Yx/4H+bmIy0WIHLmA+m3e/+P/DQDYNOOLsRMzN64sdcYShjC9oAlbaQCx0IXIYTTVbQWXnizg0pMFFEaLuGzdD7E+v7/lY/CLJZe7608Q5xlLlBFBFyJHd2cONy3K1fy+RMDUdGs962VmbNl7DHN7d2Bx367Ii1qjfUSTHoseNkTQhciRHyxg275CLeSOuT5ypZXEwVK1CrUjwHIhOsmx6GFDBF2IHEYJHUF3UymWyrjr0eHIirpVqB0R8OCqDtNs0yTHorslP1jA4r5dlagpH2Z2sigqRI6gLMIXL7jUcnuZOZL1TfKDBZw6M266fYIrnXt+Z85MvFZNc1cJOkwvSuhj/v2ohyMWuhA5grII771uNe69brXlPlHzKasiM1q0ro9TLJWx5+WTk8ScANy0KFxRHmHGqlSAV4igC5EjqCgXp0TJp2wkMk5hnFtM9duVEAdaUSpABF2IHPoKgq2Kb3lw+9fw4Pav2e6XIoqMoDUrJq+NFkNZdTCMmM0svZxxiqALkaS7M4c9vUtxpG95bcGOAMMUdq+Y/c4JzH7nhO1+qi89CoLWrJhcmM20xJUQB1pRKkAEXYg8qrgf7luOr9+6sOW1XoyIiqCZiUybYi8NqhiFsepgGGlFqQCJchECxW3HF7v9uztzWNM/1IKR2xMFQdOW1NV+pgAsa6eniWpiZNbmTsIZ6/G7VIAIuhAYRmFca/qHJgmykgKmTknj1Nl6YdHun9OIe1C1XvRERdCsRGbDEwfqImAySnqSZWnU5k7CGYNBBF1oKevz+/G9Z49hwmEmUGkCKBmIuR6jm4HXPJ+73NX+Y2fHkR8sRDasT9vpyG5WBBh3gRJaC7FFxxI/6erq4oGBgUDOLbQWVRDCYDW3GiVF2HTLQhE3wTOIaB8zdxltEwtd8JX8YCE0Pu0gKE0w1j32QuwEPWzd7oUKtkvZRDSdiJ4jomEiOkBEGw32uYOIjhPRUPXxP/0ZrhA17oyRmH/r8fvxrcfvd/26YmkiVqV2Je48vDix0M8AWMrM7xKRAuDnRPQjZt6r26+fmf/K+yEKUSLO7pX24tsNv/bhvcfQdcl5sbBiw9jtXqhga6FzhXerfyrVR9DF7YQQorXchMkwEIm4dCdI3Hl4ceRDJ6I0gH0APgjgH5j5WYPdbiKiawD8J4C1zPyKd8MUokAzdUHsWPHSM7h98Id1z3+2ex1G2mbi5v0/xc37f1q3/Y5bNuC0Mh23P78DKw7+e9322/6kDwDw6Wcfw8dffm7SttNTpuGOWysexs/v+T4+8sZhvHjBvIbfQ6GaJh91KzaM3e6FCo4yRZm5zMwdAC4CcCUR/bZul+0A5jLzRwH8BMB3jI5DRKuJaICIBo4fb6xDihBe/LDQZr99HLPfDse18uIF8/BvH7m2qWP0bI1uzXSVMHa7Fyq4DlskonsAjDGzYZWiqjV/kplnWh1HwhbjR8fGp2zLsLrlke/1AjhnSceB9jYFg/dcH/QwmkKiXIKjqbBFIpoFoMTMo0SUAfD7AL6q22c2M79e/XMlgJeaHLMQIfKDBcOMQsGYkbHof05h63YvVHDiQ58N4DtVyzsF4FFmfpKI7gUwwMxPAPhrIloJYBzASQB3+DVgIVzo0/cFQQgOW0Fn5hcAdBo8f4/m93UA1nk7NCEK+LkQGleyGSXoIQgxRcrnCk0hoWruWbFwdtBDEGKKCLrQFH6Hqn37yj/Ct6/8I1/P0Wp2vPC6/U6C0ABSy0VoCqPSqV7ysw9e5ctxgyQOi6JukaiY1iAWutAURl1YNq/qwOZVHQ13DvrQBTNwpG85Nq/qQMfY/8Wlb77q5ZCFFiO1X1qHlM8VfMMoAkbbHMGR1XbttTjx7hl88uav4LXRImZmFJw6O45Suf66nTYlhTPjE36/LU/YvKojMRbq4r5dhpmluWwGe3qXBjCiaCPlc4VAsGt84DSW+fz3TJv0xbe6EazP78eWvcd8eDfekqRCVlL7pXWIoAu+4kcCitUx7+tegK5LzsNdjw6jbDD7VN1AQVeXS5KYSe2X1iE+dCF2dHfm8PVbFxrWG3lwVQcO9y1He1uwseBJEjOp/dI6RNCFWGK0WKttbPzlP7wisLElTczs/heCd8iiqBBuflotiXvddZ4fen1+Px7ee6zl7pckLYgK3iOLokJ08UHIVVR/u7rACrTGt75x+wEAEFEXPEcEXQg3Q0OVnx0dvhxeu8DaqkJjI2Ml9Gwdrp1fELxCfOhCuFmzpvJoAaqvtxWUyhyblnRCeBBBFwQNrbSYpfeq4DUi6IKgI9eikMI0NVocQRCMEUEXBB09y+Y3XIfGDUaJT4LQDCLogqCjuzOHT109x/fztGomICQHiXIRws399wdy2vu6K4ujfsWpJy25SGgNIuhCuPnYxwI7tRqnvqZ/yNPjpokkU1LwBXG5COHmF7+oPAKiuzPnuWvkj6+6WMRc8AURdCHcfOlLlUeAeO0a2X3wuKfHEwQVEXRBsKG7M+dp1IvEnwt+IYIuCA7wMupF4s8Fv5BFUSE0qJ2ItBbsI795EwBwW++O2nMzpqbxlT9q7aKiGvXyvWePYUIX9kJwV9RL4s8FvxALXQgF2kbCdpw6W8aa/qGWNxm+r3sBZs+sXyBlwJVLRuLPBb8QC10IBZt2HjKscnjvx1ebvmbdYy+0PFrErHUcA2hvUzAyVrJ8vcSfC34iFroQCsyE8sUPXIoXP3Cp4bZiacLPIRli1joul81g8J7rLV8rnXoEv7G10IloOoBnAEyr7r+Vmb+s22cagO8CWATgTQCrmPmI56MVfGeuxletcqRvue/nNWskvPjIEABgz9wO38fghJ5l8+tqpmut7pzJ+8hlM9jTu7Rl4xSSiRML/QyApcy8EEAHgE8Q0dW6ff4cwAgzfxDAgwC+6ukohZZgJOZWz3uJmRvi8794BJ//xSO+n98pdv0xpSGyECS2FjpXmo6+W/1TqT70y/SfBLCh+vtWAN8kIuKgGpYKTfPI93onP7F3E3DrrcDnPgeMjQE33FD/ojvuqDxOnABuvrl++2c/C6xaBbzyCvCnfzppUzeA7edfg5998Cpc+uaruH/nNwEAH3njMF68YJ4Xb8kztF2OjLYBqLW1uzCbQc+y+eJmEVqCo0VRIkoD2AfggwD+gZmf1e2SA/AKADDzOBG9BeD9AE7ojrMawGoAmDPH/2p2QvR58YJ5+LePXBv0MFxhJfiC4CfkxogmoiyAxwF8npl/pXn+VwA+wcyvVv9+GcBVzHzC8EAAurq6eGBgoNFxCz5g5VpphR+9EddOK8YlCGGCiPYxc5fRNldRLsw8CmA3gE/oNhUAXFw92RQAM1FZHBUEQRBahK2gE9GsqmUOIsoA+H0AB3W7PQHgz6q/3wxgl/jPo4eZtStWsCBEAyc+9NkAvlP1o6cAPMrMTxLRvQAGmPkJAP8C4F+J6NcATgK4zbcRC74i4i0I0cVJlMsLADoNnr9H8/tpALd4OzQhabQpKYy5SBZafNl5Po5GEKKHZIoKoeH+Gz/quCbKhy6YgYc//bu+jkcQoobUchFCgxrqt3H7AdOaKNOmpPDVmz4qYYGCYIAIuhAqJIZbEBpHXC6CIAgxQQRdEAQhJoigC4IgxAQRdEEQhJgggi4IghATXBXn8vTERMcBHA3k5PWcD11lSEE+ExPkc6lHPhNj/PpcLmHmWUYbAhP0MEFEA2bVy5KKfCbGyOdSj3wmxgTxuYjLRRAEISaIoAuCIMQEEfQKDwU9gBAin4kx8rnUI5+JMS3/XMSHLgiCEBPEQhcEQYgJIuiCIAgxIdGCTkT/h4jeqDa5FgAQ0cVEtJuIXiSiA0T0haDHFDRENJ2IniOi4epnsjHoMYUJIkoT0SARPRn0WMIAER0hov1ENEREAy09d5J96ER0DYB3AXyXmX876PGEASKaDWA2Mz9PRO8FsA9ANzO/GPDQAoOICMAMZn6XiBQAPwfwBWbeG/DQQgER3QmgC8D7mHlF0OMJGiI6AqCLmVuebJVoC52Zn0GlB6pQhZlfZ+bnq7+/A+AlAIkuUM4V3q3+qVQfybWENBDRRQCWA/jnoMciJFzQBWuIaC4q/WSfDXgogVN1KwwBeAPAT5g58Z9Jlc0AvgjAeTPY+MMAniKifUS0upUnFkEXDCGi9wDYBmANM78d9HiChpnLzNwB4CIAVxJR4l10RLQCwBvMvC/osYSM32Pm3wHwBwD+surabQki6EIdVT/xNgAPM/NjQY8nTDDzKIDdAD4R8FDCwGIAK6s+40cALCWiLcEOKXiYuVD9+QaAxwFc2apzi6ALk6guAP4LgJeY+RtBjycMENEsIspWf88A+H0ABwMdVAhg5nXMfBEzzwVwG4BdzHx7wMMKFCKaUQ0mABHNAHA9gJZF0SVa0Ino+wB+CWA+Eb1KRH8e9JhCwGIAf4qKtTVUfdwQ9KACZjaA3UT0AoD/QMWHLiF6ghEfAPBzIhoG8ByAHcz841adPNFhi4IgCHEi0Ra6IAhCnBBBFwRBiAki6IIgCDFBBF0QBCEmiKALgiDEBBF0QRCEmCCCLgiCEBP+Pxj1JTIITWJdAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#red line : actual motion without any errors\n",
    "#scatter plot : the sample of the position probilities\n",
    "plt.scatter(samples[:,0],samples[:,1])\n",
    "plt.plot(trajectory[1:,0],trajectory[1:,1],'--r')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.10 64-bit",
   "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.8.10"
  },
  "vscode": {
   "interpreter": {
    "hash": "767d51c1340bd893661ea55ea3124f6de3c7a262a8b4abca0554b478b1e2ff90"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}