S42_logistics_actions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/gtbook/robotics/blob/main/S42_logistics_actions.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": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "w-CIikP0ZZnm",
    "outputId": "cef8a493-4d65-46ff-915e-fc59e6be77ef",
    "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "uef6Kglzcrha",
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.collections import PatchCollection\n",
    "from matplotlib.patches import Circle, Rectangle\n",
    "from matplotlib.transforms import Affine2D\n",
    "\n",
    "import gtsam\n",
    "from gtbook.display import show\n",
    "\n",
    "N = 5\n",
    "indices = range(1, N+1)\n",
    "u = {k:gtsam.symbol('u',k) for k in indices[:-1]} # controls u_k\n",
    "x = {k:gtsam.symbol('x',k) for k in indices}      # states x_k"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Moving in 2D\n",
    "\n",
    "> Omniwheels are a great way to move omni-directionally in 2D.\n",
    "\n",
    "<img src=\"Figures4/S42-Warehouse_robots-05.jpg\" alt=\"Splash image with robot balancing\" width=\"40%\" align=center style=\"vertical-align:middle;margin:10px 0px\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Omni Wheels\n",
    "> Omni wheels make omnidirectional movement possible.\n",
    "\n",
    "An [omni wheel](https://en.wikipedia.org/wiki/Omni_wheel) is a wheel that, instead of a tire, has small rollers along its rim. The axes for these rollers are perpendicular to the wheel's own axis. Because of this, unlike a conventional wheel, movement perpendicular to the wheel is *not* resisted, i.e., the wheel can slide effortlessly from left to right. The image below (from Wikipedia) shows an example:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a title=\"Rotacaster, CC BY 1.0 &lt;https://creativecommons.org/licenses/by/1.0&gt;, via Wikimedia Commons\" href=\"https://commons.wikimedia.org/wiki/File:Triple_Rotacaster_commercial_industrial_omni_wheel.jpg\"><img width=\"256\" alt=\"Triple Rotacaster commercial industrial omni wheel\" src=\"https://upload.wikimedia.org/wikipedia/commons/e/e3/Triple_Rotacaster_commercial_industrial_omni_wheel.jpg\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apart from being potentially amazing for parallel parking, in robotics omni wheels are very popular because they enable omni-directional movement of robot platforms. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example\n",
    "\n",
    "> Three wheels are enough!\n",
    "\n",
    "As an example, let us imagine that our warehouse robot is circular and has three omni wheels mounted on its rim. The coordinates of each wheel, in body coordinates, are\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "p^i_x &= R \\cos(\\theta^i) \\\\\n",
    "p^i_y &= R \\sin(\\theta^i) \\\\\n",
    "\\theta^i &= i \\frac{2\\pi}{3}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "where $R$ is the radius of the robot. The angle $\\theta^i$ is angle between the robot's x-axis and the wheel's axis.\n",
    "\n",
    "In code, we create a small function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = 0.50 # robot is circular with 50 cm radius\n",
    "r, w = 0.20, 0.06 # wheel radius and width\n",
    "num_wheels = 3\n",
    "\n",
    "def wheel_pose(i:int):\n",
    "    \"\"\"Return x,y, and theta for wheel with given index.\"\"\"\n",
    "    theta = float(i)*np.pi*2/num_wheels\n",
    "    px, py = R*np.cos(theta), R*np.sin(theta)\n",
    "    return px, py, theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This yields the configuration below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAssAAAHSCAYAAAAaOYYVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA/eUlEQVR4nO3dd3Rc533u++edhkHvHQTBAlKESKpBvVGFskTJYbLi4yjF9sl1ro6y7BzHjpKrc32WbyTZjmxHiuVELnJJ3HKURFYiiSKt3ptFNUoEWMAONhAA0evMvPcPQDJEAYMBZg/2lO9nrVnAzOyZ/ePmYPDgnd9+X2OtFQAAAICP8rhdAAAAAJCsCMsAAADADAjLAAAAwAwIywAAAMAMCMsAAADADAjLAAAAwAx8bhcQTVlZmW1oaHC7DAAAAKSxN954o9NaWz7dfY6EZWPMtZLukeSV9CNr7Z3TbLNO0rcl+SV1Wmsvn+15GxoatHXrVidKBAAAAKZljDkw031xh2VjjFfSvZLWS2qX9Lox5mFrbcuUbYokfVfStdbag8aYinj3CwAAACSaEz3L50lqs9butdaOSbpf0sZTtvkjSQ9aaw9KkrW2w4H9AgAAAAnlRFiulXRoyvX2ydumWiGp2BjzrDHmDWPMpx3YLwAAAJBQTvQsm2lus9Ps5xxJV0nKlvSKMeZVa+2ujzyZMTdJukmS6uvrHSgPAAAAmB8nRpbbJS2acr1O0pFptvm1tXbQWtsp6XlJZ0z3ZNba+6y1zdba5vLyaU9KBAAAABaEE2H5dUmNxpglxpiApBslPXzKNg9JutQY4zPG5Eg6X1KrA/sGAAAAEibuNgxrbcgY83lJj2li6rifWGu3G2Nunrz/+9baVmPMryVtkxTRxPRy78W7bwAAACCRjLWnthcnj+bmZss8ywAAAEgkY8wb1trm6e5juWsAAABgBoRlAAAAYAaEZQAAAGAGhGUAAABgBoRlAAAAYAaEZQAAAGAGhGUAAABgBoRlAAAAYAaEZQAAAGAGhGUAmHS0d9jtEgAASYawDCDjdfSN6Ev/9rYu++Yz2nW83+1yAABJhLAMIGONhsK695k2XfH3z+rBtw5rPGx1x6YWt8sCACQRn9sFAIAbtrx7VF/f0qpD3R9uvXhhd6eeaDmu9U2VLlUGAEgmhGUAGaX1aJ9uf6RFr+ztmnGbrz3aostXlCvg48M3AMh0hGUAGaF7cEx3Pb5T979+SOGIjbrt/q4h/eSlfbr58mULVB0AIFkRlgGktVA4op++ckD3PLlLfSOhmB/3T0+36ffPrlN5flYCqwMAJDs+YwSQtp7bdULX3vOC7tjUMqegLEkDoyF967EdCaoMAJAqGFkGkHb2dQ7qjk0tenpHR1zP88Ab7frUBQ1aU1foUGUAgFTDyDKAtNE3Mq6vPdqia/7hubiDsiRFrHTbI9sdqAwAkKoYWQaQ8iIRq3/bekh3Pb5TnQNjjj731gMn9dDbh7XxzFpHnxcAkBoIywBS2m/2deu2R7Zr+5G+hO3jG1t26JqmKmUHvAnbBwAgOdGGASAlHe4Z1ud++aY++YNXEhqUJelI74i+/9yehO4DAJCcGFkGkFKGx8L63rNtuu+FvRoZjyzYfn/w/B79wbmLVFOUvWD7BAC4j5FlACnjv946rCvvelbfebptQYOyJI2MR/T1za0Luk8AgPsYWQaQ9La19+i2R1r0xoGTrtaxadtRfeaibp3bUOJqHQCAhcPIMoCk1dE/olv+4x1tvPcl14Py+257ZLsisyyXDQBIH4wsA0g6o6GwfvziPn33mT0aGJ3bynuJ9t7hPv3HG4f0B+fWu10KAGABEJYBJJXHth/T1ze36kDXkNulzOhbj+3ShjXVyg/63S4FAJBghGUASWHX8X7d/kiLXmzrdLuUqLJ8Ht147iL5PHSxAUAmICwDcFXP0JjufmKXfvnaQYWTvBf4utVV+n83rNKikhy3SwEALBDCMgBXhMIR/eLVA/r2U7vVMzTudjlRraou0FduaNKFy0rn9DhrrcbDVqFIZOJrOKJQxGo8HFFo8vaIlayVrKzs5N8KxkhGRpLkMZLXY+T3euTzTnz1eya+93mN/B6PPB7j9D8ZADCJsAxgwb24u1O3b9quXccH3C4lqpLcgP7qmhW68dx6eacE0kjEanAspOHxsEbGIhoeD2vo/evjYQ1P3jYWWpi5oH0eo2DAq2z/5CXgUdDvVU7A98FtuVle+by0jgDAXBGWASyY/Z2D+uqjrXqy9bjbpUTl9xp96oIG/Y/LlkpG2tc5oL6RkPqGx9U/EtLgaEjJ1DESilgNjIQ0MBJ95pCcgFf5QZ8Ksv3KD/qUH5z4mhfwMToNADMgLANIuIHRkP7xqd3655f2ayy8sCvvzdWZi4r0e2fVKjvg1VM7Otwux1FDY2ENjYV1vG/0Q7d7jJQX9KkkJ6Di3IBKcye++hmJBgDCMoDEiUSsHnijXd98bKc6B0Znf4CLyvKydP2aKq2sKpAkhcJJNHScYBEr9Q2H1Dcc0v4pU/blB30fBOeSyRBNKweATENYBpAQW/d367ZHWvTu4V63S4kq6PfoypUVunBZ2Yf6kiH1j4TUP/LbAO0xE33cFQVBVRZkqTwvi/AMIO0RlgE46kjPsP5uyw498s4Rt0uJykhqbijR+qZK5WXxVhiLiJU6B8bUOTCmliO/Dc+VBUFVFgRVlsfIM4D0w28IAI4YGQ/r+8/t0Q+e26vh8bDb5US1pCxX16+pVk1RttulpLSp4Xn7kT55PVJlQVB1xdmqLcpRdsDrdokAEDfCMoC4PfzOEX1jyw4d7hl2u5SoinL8um51tdbUFrpdSloKR6QjPSM60jMi6aRKcgOTwTlbxbkBt8sDgHkhLAOYt/cO9+q2R7br9f0n3S4lqoDXo8tWlOvSxjJmeFhA3YNj6h4c07b2XuVmeVVblK360hxV5AfdLg0AYkZYBjBnJ/pH9feP7dR/vHEoqeYbPpWRdMaiIn3s9CoVZvvdLiejDY6Gtev4gHYdH1Be0KeG0hwtLs3l/wVA0iMsA4jZWCiif35pn/7p6Tb1j0ZfAMNtdcXZumFNtepLc90uBacYGAnpvcN9eu9wn0py/Wooy1VDaa6CfnqcASQfwjKAmDzZclxf29yqfZ2DbpcSVX7Qp481Vems+iIZw1Rwya57cFzdgz1662CPqgqCWl6Rp9qibFYUBJA0CMsAomrr6Nftm1r1/K4TbpcSlc9jdPHyMq1bWa4sHyOUqcZa6WjviI72jign4NWy8jwtr8hjRg0AriMsA5hW79C4/uHJXfrFqwcUSubGZElN1QXasKZaJcy4kBaGxsJ693Cvth/pVV1xjhor81RZwEmBANxBWAbwIeGI1b++dkB3P7FLJ4fG3S4nqsqCLF2/pkbLK/LcLgUJELHSwe4hHeweUmG2X42VeVpalsvCJwAWFGEZwAdebuvU7ZtatONYv9ulRJUT8OrqVZU6b0mJPPQlZ4Te4XFt3X9S77b3amVVvhor82i3AbAgCMsAdLBrSF99tEWPtxx3u5SoPEY6f0mprl5VSS9rhhoNRbStvVctR/u0vCJPq6oKeC0ASCjCMpDBBkdD+qdn2vTjF/dpLBRxu5yoGivytGFNNb2rkCSFwlY7jvZr17F+LSnLVVNNgfKDzNkMwHmEZSADWWv1wBvt+tZjO9XRP+p2OVGV5ga0YU21VlUXuF0KklDESntODGpv56DqS3K0pq5QBYRmAA4iLAMZ5s2DJ3Xbw9v1Tnuv26VEleXz6IqVFbpoeal8Hk7oQnTWSge6hnSoe0hLy/O0praQ9gwAjiAsAxniWO+I7tzSqofeOSKbxDPBGUnnLC7W+qZKPlbHnEWs1NYxoP2dg1pRla+m6gIFfPyxBWD+CMtAmhsZD+u+5/fq+8/t0dBY2O1yolpcmqMb1taotijb7VKQ4kIRq5YjfWrrGNCq6nytrMxnyjkA80JYBtLYo9uO6u+2tKr95LDbpURVmO3XtaurdEZdkdulIM2MhSJ651Cvdh8f0Nq6Qi0tZ05uAHNDWAbS0PYjvbrtkRb9Zl+326VE5fcaXdZYrstWlMvPqB8SaGgsrFf3dqutY0DNDSWs9gggZoRlII10DYzq7x/fqX97/ZCSfIVqra0r1LWnV6koh9CChdM5MKbHth9TY0We1tYV0c8MYFaEZSANjIcj+unL+3XPU7vVPxJyu5yoaouydcPaai0uzXW7FGQoa6Vdxwd0oGtIZ9YXaRmtGQCiICwDKe6ZHR2649EW7T0x6HYpUeVl+XRNU6XOXlzMEtVICqOhiF6bbM04l9YMADMgLAMpqq1jQF99tEXP7jzhdilReT1GFy0r1RUrKxT0M+8tkk/XZGvGquoCraktlNfDH3MAfouwDKSY3uFx3fPkbv381f0aDyd3Y/KqqnxtWFOt0rwst0sBorJWajnSp8Mnh3XB0hJeswA+QFgGUkQkYvWvvzmou5/Ype7BMbfLiaoiP0vXr61WY0W+26UAc9I7PK7HW44zygzgA4RlIAW8sqdLt29qUevRPrdLiSrb79VVqyp0/pJSQgZSFqPMAKZyJCwbY66VdI8kr6QfWWvvnGG7cyW9KukPrLUPOLFvIJ0d6h7S1ze3ast7x9wuJSqPkc5bUqKrV1UqJ8Df4EgPU0eZ19YWysMfgEBGivu3mjHGK+leSesltUt63RjzsLW2ZZrtviHpsXj3CaS7obGQvvvMHv3whb0aDUXcLieqZeW5un5tjaoKgm6XAjju/VHm430junh5mfKy+GMQyDRO/NSfJ6nNWrtXkowx90vaKKnllO3+QtKvJJ3rwD6BtGSt1X++dVjf+PUOHe8bdbucqEpyA9qwukpNNYVulwIkXNfAmLa8e1TnLylVfWmO2+UAWEBOhOVaSYemXG+XdP7UDYwxtZJ+T9KVIiwD03r7UI9ue2S73jrY43YpUWX5PFq3skIXLyuVjyWqkUHGw1YvtnWqsT9PZ9cX05cPZAgnwvJ07xanzmf1bUn/j7U2bGZZjMAYc5OkmySpvr7egfKA5NbRN6I7t+zQf759WDaJZ4Izks6qL9bHTq9UftDvdjmAa3YfH9CJ/lFdvLxMhdn8LADpzomw3C5p0ZTrdZKOnLJNs6T7J4NymaQNxpiQtfa/Tn0ya+19ku6TpObm5iSODkB8RkNh/eiFffruM20aHAu7XU5U9SU5umFtteqK+fgZkKSeoXE99t4xnbukREvKWLodSGdOhOXXJTUaY5ZIOizpRkl/NHUDa+2S9783xvyLpE3TBWUgU2x596i+vqVVh7qH3S4lqsJsvz52epXOXFTkdilA0glFrF7Z06WugVGdXV/MbBlAmoo7LFtrQ8aYz2tilguvpJ9Ya7cbY26evP/78e4DSBc7jvXptodb9MreLrdLicrvNbpkebkuX1GugI++ZCCaXccH1DcyrouXlynLx5LuQLoxNombJJubm+3WrVvdLgOIW/fgmO56fKfuf/2QwpHk/ZmTpDW1hbpudZWKcgJulwKklLygT5c3lqswhz5mINUYY96w1jZPdx8TRgIJFApH9LNXDuiep3ard3jc7XKiqikM6vq1NfRfAvM0MBLSYy3HdNGyUvr7gTRCWAYS5LldJ3THpha1dQy4XUpUuVk+XbOqUuc0FMszy2w1AKILha1e2N2pNbWFWl3LHORAOiAsAw7b1zmoOza16OkdHW6XEpXXGF24rFRXnlahoJ8+S8Ap1krb2nvVNzyuC5aWcuIfkOIIy4BD+kbG9Y9P7da/vLxf4+Hk7kteWZmv69dUqyw/y+1SgLS1v2tII6GwLm0sl58FfICURVgG4hSJWP3b1kO66/Gd6hwYc7ucqMrzsnT92mqtqMx3uxQgIxzrHdVTrce1biWf4ACpirAMxOE3+7p12yPbtf1In9ulRBX0e3TVaZW6YGkpS/QCC6x7cFyPtxzXFSvLWf0SSEGEZWAeDvcM6+ubW/XotqNulxKVx0jNDSVav6pSuVn8uANuGRgJ6YmWiRHmklymZQRSCb89gTkYHgvre8+26b4X9mpkPOJ2OVEtLcvV9WurVV2Y7XYpACSNjEf0ZOtxXb6iXJUFQbfLARAjwjIQo4fePqw7t+zQ0d4Rt0uJqjjHr+tWVzNtFZCEQmGrZ3Z06OLlZVpUwlzMQCogLAMx+NuHt+tfXt7vdhlRBXwerVtRrkuWl8nHmfdA0opY6aW2TgIzkCL4jQrE4BPn1ClZz4szks5aVKQvXb1C61ZWEJSBFPB+YD7YNeR2KQBmwcgyEIPVtYX6xDl1+vet7W6X8iGLirN1w9oaRqeAFBSx0kt7OmVVqsWlLDMPJCvCMhCjv/7Yadry7jH1j4bcLkUFQZ8+dnqVzlxUJMMS1UDKslZ6ZU+XjIzqS/mjF0hGfF4LxKg8P0ufv3K5qzX4PEbrVpbrS+tX6qz6YoIykAYiVnp5T6cOddOSASQjwjIwB3968RItKXPn49LTawr0xatX6JqmKgV8/OgC6eT9HubDPcNulwLgFPzGBeYg4PPoyxtWLeg+qwuD+rNLluiPz1+sYhYzANJWxEov7e5U58Co26UAmIKwDMzR1U2VurSxLOH7yQl4tfHMGn3uiuVaWp6X8P0BcF8oYvXszhPqHR53uxQAkwjLwDx85YYm+RI0l5zXGF28rFR/tX6lzl9SKg99yUBGGQtF9OzODg2NuX8yMQDCMjAvjZX5+pMLFjv+vCsq8/QXVy3X9WtrlB3wOv78AFLD4GhYz+w4odFQ2O1SgIxHWAbm6YtXr1Bxjt+R5yrLC+jTFy7Wf79oiSryg448J4DU1js8rud3dSoUjrhdCpDRCMvAPBXm+PWl9Svieo6g36MNq6v0hatW6LSqAocqA5AuTvSP6qU9XbLWul0KkLEIy0Ac/uj8xTqtKn/Oj/MY6dyGYn1p/Upd0lgub7KupQ3AdYdPDuvNgyfdLgPIWIRlIA5ej9FXbmia02OaFxfrf17VqN87q055WSyiCWB2O48NqK1jwO0ygIzEb2ogThctL9M1TZV6vOV41O1qi7L119eulNcY9Y9wljuAudm6v1sF2T7OawAWGCPLgAP+9/VNM66ql+336kvrV+ipL12mwqCfoAxgXiJWenF3pwZHeQ8BFhJhGXBAfWmOPnvJkg/dZoz0u2fW6OlbLtf/vKpR24/262jviEsVAkgHI+MRvbD7BDNkAAuIsAw45PNXLFdFfpYk6Yy6Qj1w80X69o1nqbowW3tODGjnsX6XKwSQDroHx/Xq3m63ywAyBj3LgENys3z62985XYOjIX3inDqZyZX3OgdG9fo+frEBcM7B7iG9d7hXq2sL3S4FSHuEZcBBG9ZUf+j6aCisl9o6FWGKVAAOe/dwr8rzs1RZwAl/QCLRhgEk0Gt7uzU4ynK1AJxnrfTynk6NjPMeAyQSYRlIkJ3H+tV+ctjtMgCkseGxiF7Z2+V2GUBaIywDCdA9OKa3WHELwAI42jOiliN9bpcBpC3CMuCwsVBEL9KnDGABbWvvUefAqNtlAGmJsAw47PX93Rpg4REACyhipZfaOjUWYv5lwGmEZcBBbR0DOtA15HYZADLQ4GhYr+2jfxlwGmEZcMjAaEhv0qcMwEWHuoe1r3PQ7TKAtEJYBhzy2t4uhcI0KgNw1xsHTmp4jOnkAKcQlgEH7D7er+N9nFwDwH1joYh+s59VQwGnEJaBOA2MhvTWoR63ywCADxw+STsG4BTCMhAn2i8AJCPaMQBnEJaBONB+ASBZ0Y4BOIOwDMzTIO0XAJIc7RhA/AjLwDy9vr+b9gsASe/NAyc1GqIdA5gvwjIwD4d7hnWkZ8TtMgBgVqOhiLa197pdBpCyCMvAHIUjVm8cYPERAKmjrWNA3YNjbpcBpCTCMjBHLUf6NDAScrsMAIiZtdJWTvYD5oWwDMzBwGhILUf5OBNA6ukcGNOeEwNulwGkHMIyMAdvHDipcMTtKgBgft4+2KOxEG9iwFwQloEYHe4Z1uGTw26XAQDzNnGyX4/bZQAphbAMxCASsXqTk/oApIHdHQPqGeJkPyBWhGUgBntODKifk/oApAFrpXeYSg6IGWEZmEUoHNF7R/jFAiB9HD45rI5+5ooHYkFYBmax41i/hsc4IQZAennnEIMAQCwIy0AUo6GwWo/2uV0GADjuRP+oDvdw0jIwG8IyEMX2I30aD1u3ywCAhHjnUI+s5T0OiIawDMxgcDSk3cf73S4DABKmZ2hc+7uG3C4DSGqEZWAG7x7uZQESAGlvW3uPIhFGl4GZEJaBaQyMhrS/c9DtMgAg4QZHw9rL+x0wI8IyMI3Wo31ioAVApmg52kfvMjADwjJwipHxsPaeGHC7DABYMAMjIR3spncZmA5hGTjFjmP99CoDyDgtR5gmE5gOYRmYYiwUYQYMABnp5NC4jjDvMvARjoRlY8y1xpidxpg2Y8yt09z/x8aYbZOXl40xZzixX8Bpuzv6mVcZQMbazugy8BFxh2VjjFfSvZKuk9Qk6Q+NMU2nbLZP0uXW2rWS7pB0X7z7BZwWCke08xijygAy14n+UXX0j7hdBpBUnBhZPk9Sm7V2r7V2TNL9kjZO3cBa+7K19uTk1Vcl1TmwX8BR+zoHNTJOszKAzEbvMvBhToTlWkmHplxvn7xtJp+VtMWB/QKO2kmvMgDoSM+I+kfG3S4DSBo+B57DTHPbtE2fxpgrNBGWL5nxyYy5SdJNklRfX+9AecDsjvWOqG845HYZAJAUdncM6Oz64g+uN9z6aEL2s//O6xPyvICTnBhZbpe0aMr1OklHTt3IGLNW0o8kbbTWds30ZNba+6y1zdba5vLycgfKA2a3u4NRZQB4394TgwqzMhMgyZmw/LqkRmPMEmNMQNKNkh6euoExpl7Sg5I+Za3d5cA+AccMjYV0+CTTJQHA+8ZCER3oYglsQHIgLFtrQ5I+L+kxSa2S/t1au90Yc7Mx5ubJzb4iqVTSd40xbxtjtsa7X8ApbR0DLG0NAKfYdZyVTAHJmZ5lWWs3S9p8ym3fn/L9n0n6Myf2BTgpErHaw9LWAPAR3YNj6hoYVWleltulAK5iBT9ktPaTwxoeY7o4AJjO7g4GEwDCMjIaJ/YBwMwOdg1pLMSAAjIbYRkZa3A0pON9o26XAQBJKxSxOtjNiX7IbIRlZKz9nOkNALPa1znkdgmAqwjLyFj7+QUAALM60c8ncMhshGVkpJODY+odZjlXAAAQHWEZGWkfLRgAACAGhGVkHGstK1MBAICYEJaRcY73jTK3MgAAiAlhGRlnXyejygAAIDaEZWSUSMTqcM+w22UAAIAUQVhGRukcGGU1KgAAEDPCMjJKO6PKAABgDgjLyCjtJwnLAAAgdoRlZIze4XENjITcLgMAAKQQwjIyRvtJlrcGAABzQ1hGxjhMCwYAAJgjwjIywsh4WF2DY26XAQAAUgxhGRnhWO+IrHW7CgAAkGoIy8gIx/tG3C4BAACkIMIyMsLx/lG3SwAAACmIsIy0NzQWYso4AAAwL4RlpL3jfYwqAwCA+SEsI+110K8MAADmibCMtEe/MgAAmC/CMtLa4Cj9ygAAYP4Iy0hrJxhVBgAAcSAsI611DRKWAQDA/BGWkda6B8fdLgEAAKQwwjLSlrVWJwfH3C4DAACkMMIy0lbfSEihiHW7DAAAkMIIy0hb3YwqAwCAOBGWkbYIywAAIF6EZaQtwjIAAIgXYRlpyVqrk0OEZQAAEB/CMtLSwGhIoTAn9wEAgPgQlpGW+lniGgAAOICwjLREWAYAAE4gLCMt9Y2wch8AAIgfYRlpqZ+wDAAAHEBYRlqiDQMAADiBsIy0EwpHNDgadrsMAACQBgjLSDsDo4wqAwAAZxCWkXZowQAAAE4hLCPtDI3RggEAAJxBWEbaGR4nLAMAAGcQlpF2hhlZBgAADiEsI+2MMLIMAAAcQlhG2qFnGQAAOIWwjLRDzzIAAHAKYRlpJRyxGgtF3C4DAACkCcIy0gqjygAAwEmEZaSVUcIyAABwEGEZaSUUsW6XAAAA0ghhGWmFfmUAAOAkwjLSCiPLAADASYRlpJVQmJFlAADgHMIy0sp4mJFlAADgHMIy0koowsgyAABwDmEZaYWRZQAA4CTCMtIKPcsAAMBJhGWklbBlZBkAADjHkbBsjLnWGLPTGNNmjLl1mvuNMeY7k/dvM8ac7cR+gY8gKwMAAAfFHZaNMV5J90q6TlKTpD80xjSdstl1khonLzdJ+l68+wWmQ1YGAABOcmJk+TxJbdbavdbaMUn3S9p4yjYbJf3MTnhVUpExptqBfQMfQhcGAABwkhNhuVbSoSnX2ydvm+s2kiRjzE3GmK3GmK0nTpxwoDwAAABgfpwIy2aa204d34tlm4kbrb3PWttsrW0uLy+PuzgAAABgvpwIy+2SFk25XifpyDy2AeLmme7PMgAAgHlyIiy/LqnRGLPEGBOQdKOkh0/Z5mFJn56cFeMCSb3W2qMO7BsAAABIGF+8T2CtDRljPi/pMUleST+x1m43xtw8ef/3JW2WtEFSm6QhSX8a736B6RjD0DIAAHBO3GFZkqy1mzURiKfe9v0p31tJn3NiX0A0XvowAACAg1jBD2nF5yUsAwAA5xCWkVYCXl7SAADAOSQLpBXaMAAAgJMIy0grfkaWAQCAg0gWSCt+epYBAICDCMtIKz5GlgEAgINIFkgrfnqWAQCAgwjLSCv0LAMAACeRLJBWgn6v2yUAAIA0QlhGWgn6PWLFawAA4BTCMtKKMUZBPy9rAADgDFIF0k42rRgAAMAhhGWkHfqWAQCAUwjLSDuMLAMAAKcQlpF2cgI+t0sAAABpgrCMtJMd4GUNAACcQapA2mFkGQAAOIWwjLSTHyQsAwAAZxCWkXbysnzysDAJAABwAGEZaccYozxGlwEAgAMIy0hLBUG/2yUAAIA0QFhGWqJvGQAAOIGwjLSUz8gyAABwAGEZaamAkWUAAOAAwjLSUkE2I8sAACB+hGWkpaDfy0p+AAAgbqQJpK3inIDbJQAAgBRHWEbaKsklLAMAgPgQlpG2CMsAACBehGWkLcIyAACIF2EZaSsn4FPQz0scAADMH0kCaY3RZQAAEA/CMtIaYRkAAMSDsIy0VpqX5XYJAAAghRGWkdbK87JkjNtVAACAVEVYRloL+DwsTgIAAOaNsIy0V1lAKwYAAJgfwjLSXkVB0O0SAABAiiIsI+1V5NO3DAAA5oewjLTn93qYQg4AAMwLYRkZoSKfvmUAADB3hGVkhKpC+pYBAMDcEZaRESryg/J5aFwGAABzQ1hGRvB6DKPLAABgzgjLyBi1xdlulwAAAFIMYRkZo7YomynkAADAnBCWkTGCfq9KmUIOAADMAWEZGYVWDAAAMBeEZWSUuqIct0sAAAAphLCMjFKY41de0Od2GQAAIEUQlpFx6ksYXQYAALEhLCPjLCnLdbsEAACQIgjLyDiF2X6V5PrdLgMAAKQAwjIy0uJSRpcBAMDsCMvISA2luSxQAgAAZkVYRkbKDnhVWZDldhkAkPR8HkYWkNkIy8hYDbRiAMCsWMwJmY6wjIy1qCSHERMAmEUDMwghwxGWkbH8Xo8WMecyAMwoO+BRdUHQ7TIAVxGWkdEaK/PcLgEAktby8nx5+AQOGY6wjIxWlpfFnMsAMA2PkZZXMKAAxBWWjTElxpgnjDG7J78WT7PNImPMM8aYVmPMdmPMF+LZJ+C05RX5bpcAAEmntjhb2QGv22UArot3ZPlWSU9ZaxslPTV5/VQhSX9lrV0l6QJJnzPGNMW5X8AxDaU58nv5mBEAplpRyUACIMUfljdK+unk9z+V9LunbmCtPWqtfXPy+35JrZJq49wv4Bif16Ol5XzUCADvK8j2qZIT+wBJ8YflSmvtUWkiFEuqiLaxMaZB0lmSXotzv4CjONEPAH6rkfY04AO+2TYwxjwpqWqau748lx0ZY/Ik/UrSX1pr+6Jsd5OkmySpvr5+LrsA5q0g6Fd1YVBHe0fcLgUAXOXzGC05ZW7l/Xde71I1gPtmDcvW2qtnus8Yc9wYU22tPWqMqZbUMcN2fk0E5V9aax+cZX/3SbpPkpqbm+1s9QFOWVmVT1gGkPGWVeQq4GOyLOB98f40PCzpM5Pff0bSQ6duYIwxkn4sqdVae3ec+wMSpqYoW8U5TCMHIHN5jHRaVYHbZQBJJd6wfKek9caY3ZLWT16XMabGGLN5cpuLJX1K0pXGmLcnLxvi3C+QEE01/JIAkLkWl+YqN2vWD52BjBLXT4S1tkvSVdPcfkTShsnvX5TEvFxICfUlOXon2KuBkZDbpQDAgmPAAPgompKAKYwxaqrmlwWAzFNXnK3CbFrRgFMRloFTLC3LVXaAHw0AmeV0RpWBaZEIgFN4PIYTXABklKrCLJXmZbldBpCUCMvANJZX5DF1EoCM0VRd6HYJQNIiDQDT8Hs9WlXNClYA0l9lQZaqClnaGpgJYRmYwcrKfOUEvG6XAQAJdcaiIrdLAJIaYRmYgc/r0epaPpoEkL7qirNVRq8yEBVhGYhiWXmuCrKZoB9A+jGGUWUgFoRlIApjjNbWFrldBgA4bklZLvMqAzEgLAOzqC/NUUluwO0yAMAxXo+0to42MyAWhGUgBmfVF7ldAgA4prEyXzkBWsyAWBCWgRhUFgRVXcTUSgBSX8DnYbU+YA4Iy0CMzllcLI9xuwoAiM8ZdYXK8jEtJhArwjIQo4KgX6dVMxoDIHWV5Pq1vCLP7TKAlEJYBuZgdU2BcrMYkQGQmpobSmQMH5EBc0FYBubA5/XorEXFbpcBAHO2tDyXBUiAeSAsA3NUX5qj6kJO9gOQOgI+j85kARJgXgjLwDyczcl+AFLIGXWFCvppIQPmg7AMzENhNif7AUgNnNQHxIewDMzT6poC5QeZ1B9A8vIY6bwlpZzUB8SBsAzMk8/r0flLS8TvIADJqqmmQCW5AbfLAFIaYRmIQ0V+UCsq890uAwA+oijHr9U1hW6XAaQ8wjIQpzPqCmnHAJBUPEa6YGmpPJyJDMSNsAzEiXYMAMmG9gvAOYRlwAG0YwBIFrRfAM4iLAMOoR0DgNtovwCcR1gGHOLzenTRslIWKwHgmtW1hbRfAA4jLAMOKs3L0hksKQvABdWFQa2upf0CcBphGXDYquoC1RQF3S4DQAbJDnh04bJSt8sA0hJhGUiAC5aWKifgdbsMABnAGOnCpWUK+nnPARKBsAwkQNDv1UXLS5lODkDCNVUXqKqQT7OARCEsAwlSkR/UGvoHASRQeX6W1tbxPgMkEmEZcJC19kPXT68pUFVhlkvVAEhnWT6PLl5eKsNHWEBCEZYBB33j1zv1uX99U4d7hiVJxhhdtKxMuVn0EgJwjjHSRctLlRNgbncg0QjLgEP2dw7qJy/u06Pbjuqqu57V3U/s0vBYWEG/V5c1lsvHBMwAHHJWfZGqC7PdLgPICIRlwCFffbRVY+GIJGlkPKLvPLVbV931rB56+7CKcwO6YCnTOgGI39LyXJ1WVeB2GUDGICwDDnhh9wk92Xr8I7cf6R3RF+5/W5/43svqHR7X6lp+wQGYv7K8gM5rKHG7DCCjEJaBOIUjVndsaom6zdYDJ/U7976on79yQPlBegwBzF1ulleXrSiXh5YuYEERloE4/eLVA9p1fGDW7ayV/uONdt364Db9Zl+XQpHIAlQHIB34PEaXNpaz8AjgAsIyEIeeoTH9w5O75vSYwdGw/uvtI/rOU7vVcqQvQZUBSCcXLC1VSW7A7TKAjERYBuJw9xO71DM0Pq/Hdg6M6RevHdBPXtqn430jDlcGIF2cVV+k+tIct8sAMhZhGZinncf69cvXDsb9PG0dA/rHp3fr4XeOaGgs5EBlANLFadX5WlXNicGAmzjTCJinOza1KByxs28Yg4iVXt3bpXcO9ejqVRU6f2mpPKzKBWS0hrIcnV1f7HYZQMZjZBmYh8e3H9OLbZ2OP+/weFiPbDuqf3x6t9o6Zj9pEEB6qi4M6oIlzM0OJAPCMjBHY6GIvra5NaH7ON43qp+8tE8/f2W/ugZGE7ovAMmlNC+gSxvLmCIOSBK0YQBz9OMX9+lA19CC7Kv1WL92dQzo4mWlumJlhbKYNgpIa/lBny5fUS6fl7EsIFnw0wjMQUf/iO59pm1B9xmOWD2/u1N3P7FLbxzolrXO9EkDSC65WV5deVoFcykDSYawDMzBt369UwOj7sxY0T8a0q/ePKzvPrtHB7oGXakBQGK8H5Rzs/jAF0g2hGUgRtvae/TAm+1ul6HDPcP6wfN7df/rB9U7PL85ngEkj5zARFDOD/rdLgXANPgTFojRbY+0KJk6ILa196r1aJ8uW1GuyxrL5afHEUg52QGPrlxFUAaSGb9dgRg89PZhvXHgpNtlfMR42Oqp1g79wxO7tK29x+1yAMxBTsCrq1dVqoCgDCQ1wjIQg77hcWUn8Uk3PcPjuv/1Q7rv+T060jPsdjkAZpGb5dVVjCgDKYGwDMTgUxc26OlbLtfGM2vcLiWq/V1DuveZNj34ZrtrJyICiC4v6NNVqyoJykCKICwDMaouzNY9N56lX/35hVpbV+h2OTOykrYeOKm7Ht+pF3afcGxJbgDxK8n165qmSuUx6wWQMgjLwByds7hED33uYn3zE2tVnp/ldjkzGg1FtOW9Y7rnqV3acbTP7XKAjFdVmKWrVlUyjzKQYgjLwDwYY/TJ5kV65pZ1uvnyZQr4kvdHqXNgTD979YD++aV96ugbcbscICMtLs3RuhUVzFoDpCB+aoE45GX5dOt1p+mJL16m9U2VbpcT1e6OAX3n6d16ZNsRDY+F3S4HyBgrq/J08fIyeTzG7VIAzANhGXDA4tJc/fDTzfrln52vlZX5bpczo4iVXtnTpbue2KlX93YpkkwTRwNpaG1doc5ZXOJ2GQDiQFgGHHTx8jJt/sKlun3j6SrKSd4z3YfGwnr4nSP6p6fbtOfEgNvlAGnHY6QLlpZodW3yngwMIDaEZcBhXo/Rpy9s0LO3rNNnLlwsXxJ/9Hqsb0Q/fnGffvHqAXUPjrldDpAWsgMeXbWqUkvL89wuBYADjE3ij2Gbm5vt1q1b3S4DiMuu4/26Y1OLXtjd6XYpUfk8RhcvL9O6leXK8nG2PjAfJbl+XbaiXDkBpoYDUokx5g1rbfN09zGyDCTYisp8/fyz5+uHn25WQ2mO2+XMKBSxem7XCd39xC69eeCkkvkPaSAZLS7N0dWrKgnKQJqJKywbY0qMMU8YY3ZPfi2Osq3XGPOWMWZTPPsEUtX6pko9/sXLdet1pyX1ggT9IyE98Ga7vvfcHh3sHnK7HCDpGSOdsahQFy8vk4+p4YC0E+9P9a2SnrLWNkp6avL6TL4gqTXO/QEpLeDz6ObLl+mZW9bpk811SuJ2ZrWfHNYPntujf996SH3D426XAyQln9fo0sYynV7DiXxAuoo3LG+U9NPJ738q6Xen28gYUyfpekk/inN/QFooz8/SNz9xhh763CVqXjzjBzKus5LePtSju5/Ypad3dGg8HHG7JCBplOT6dd3qKtUVJ297FYD4xRuWK621RyVp8mvFDNt9W9LfSOI3LTDFmrpCPfDnF+k7f3iWagqDbpczo7FwRE+2Hte3n9yldw/3ul0O4LqVVXm6pqlK+cHknSISgDNmbZw0xjwpqWqau74cyw6MMTdI6rDWvmGMWRfD9jdJukmS6uvrY9kFkPJ+54warV9Vqe8/t0c/eH6PRsaT8+/Kk0Pj+j+/OahXy3J1w9pqVRdmu10SsKACPo/OX1KiRSWMJgOZIq6p44wxOyWts9YeNcZUS3rWWrvylG3+TtKnJIUkBSUVSHrQWvsnsz0/U8chEx3uGdbfbW7Vpm1H3S4lKiOpuaFE1zRVKjeJT1gEnFKWF9DFy8t4vQNpKNrUcfGG5W9J6rLW3mmMuVVSibX2b6Jsv07SLdbaG2J5fsIyMtnr+7t12yPb9d7hPrdLiSro9+jK0yp14dJSeZP5jEUgDk01BVpbWygPr3EgLSVynuU7Ja03xuyWtH7yuowxNcaYzXE+N5DRzm0o0cOfu0Tf+P01KsvLcrucGY2MR7T53aO656nd2nms3+1yAEflZnl11aoKnbmoiKAMZChW8ANSQP/IuP7x6Tb9y0v7NZbkM1KsqMzT9WtqVJ6fvAEfiMWKyjyduaiIuZOBDJCwNoxEIywDH7a/c1BffbRFT7Z2uF1KVF5jdMHSEl21qlJBP0tnI7XkZnl1wdJSVRYk7ww1AJxFWAbSzPO7TuiOTS3a3THgdilR5Qa8Wt9UpeaGYnkMH2Ej+TGaDGQmwjKQhkLhiH7+6gF9+8nd6k3yFfaqC4O6fm21lpbluV0KMK28oE/nLylhNBnIUIRlII2dHBzT3U/s0r/+5qDCkeT9eZak1TUFum5NtYpzAm6XAkiSvB6pqbpQTTUFzOYCZDDCMpABdh7r122PbNfLe7rcLiUqn8foksYyrVtRoYCPj7rhnpqioM5ZXMwqfAAIy0Am+fV7x/T1za062D3kdilRFQR9+tjpVTpzUZEM/cxYQHlBn85ZXKzaIlagBDCBsAxkmNFQWD9+cZ/ufbpNg2Nht8uJqr4kRzesrVZdMcsHI7F8HqOmmgKtqqblAsCHEZaBDNXRN6Jv/HqnHnyrXUn8oy4j6cxFRfrY6ioV8JE4EqChNEdrFxUpj6WqAUyDsAxkuHcO9ei2R7brzYM9bpcSVcDn0boV5bpkeRlTd8ER1UVBnVlXpOJcTioFMDPCMgBZa/XQ20d055YdOtY34nY5UZXkBnTd6iqdXlPodilIUaV5AZ21qEgVTAUHIAaEZQAfGBoL6XvP7tF9z+/VaCi5l85eWp6rG9bUqKqQwIPYFGT7dEZdkRaV0AMPIHaEZQAf0X5ySF/f3KrN7x5zu5SoPEY6t6FE61dVKod+U8ygINunpuoCLSnLZXYVAHNGWAYwo9f2dum2R1rUcrTP7VKiyvZ7ddWqCp2/pJSZDPCBktyATq8pUF1xNiEZwLwRlgFEFYlY3f/6Id31+E51DY65XU5UFflZun5NtRor890uBS6qLMjS6TWFtOgAcARhGUBM+kbGdc+Tu/WzV/ZrPJy87w2SdFpVvjasqVZZXpbbpWAB1RVnq6mmgP93AI4iLAOYkz0nBvTVTS16ZucJt0uJyusxumhpqa44rUJBv9ftcpAgAZ9HS8py1ViZxzzcABKCsAxgXp7Z2aGvbmrRnhODbpcSVW6WT9c0VeqcxcXy0LeaNkpy/WqszNfikhzm3QaQUIRlAPM2Ho7oZ68c0D1P7lLfSMjtcqKqKQrqhjU1aijLdbsUzJPXI9WXTIwi02oBYKEQlgHErXtwTH//+E7d/5uDiiTv24YkaU1toa5bXaWiHFZtSxUluQEtKctVQ1mOsny01ABYWIRlAI5pOdKn2zdt16t7u90uJSq/1+jSxnJd1liugI+P8JNRbpZXS8pytbg0V4XZ9CIDcA9hGYDjtrx7VF/b3Kr2k8NulxJVYbZf155epTMWFbldCiRl+TyqL81RQ2muyvNpswCQHAjLABJiZDysH72wV999do+GxsJulxPV4pIc3bC2RrXF2W6XknGyAx7VFuWotjhb1QVBeVhUBkCSISwDSKjjfSP6xpYd+s+3DyuJ31JkJJ1dX6xrTq9UPlOQJVRxjl+1xdmqLcpWKSfqAUhyhGUAC+Ktgyf1t4+06J1DPW6XElWWz6MrVlboouWl8nnoZ3aCz2NUnp/1QUDOzfK5XRIAxIywDGDBWGv14JuH9c3Hduh436jb5URVmhvQdaur1VRT4HYpKcfnMSrLD6giP6iKgiyV5WbRXgEgZRGWASy4wdGQ7n2mTT96cZ/GQhG3y4lqeUWerl9TrcqCoNulJC2fx6g0L6DKAsIxgPRDWAbgmkPdQ/rao6369fZjbpcSlcdI5y0p1fpVlcoOZPY8vz6PUXFuQCW5fpXkZqkkJ6CCbJ8MqyMCSFOEZQCue3lPp25/pEU7jvW7XUpUBUGfbjyvXpcsK9PAWEgDIyGFkn0VlnkyZmKp8PygTwVBv0pyAwRjABmJsAwgKYQjVv/6m4P6hyd2qXtwzO1yolpZma+vfLxJFy8v0+BoSP0jIfWPjKtvZFx9IxMhengsnPRB2mOk7IBXuYGJUJwf9E+E42y/8rN8tFIAgAjLAJJM7/C4vv3kLv38lQNJHzbXN1Xqf1+/SotLc6e9fzwc0dBYWCPjYQ2PhTU8PnEZGQtrNBTReDiiUMROfA1PfJ3vP9mYiRYJv9cjn9fI5/HI7zUK+DzK9nsV9HuVHfAq2+9VTmDietCf2S0lABALwjKApNTW0a/bN7Xq+V0n3C4lqoDXoz+9pEF/cWWj8hyYEi0SsRqPTITniLWy0sT81O+/HU8O9hojeYyRzzN58TLNHQAkAmEZQFJ7esdxfXVTq/Z2DrpdSlTl+Vn664+t1H87p46eXgBII9HCMsMUAFx35WmVeuyLl+nLG1YpP5i8i1mc6B/V3zywTRvvfUkn+pN7DmkAgDMIywCSgt/r0f992VI9c8s63XjuIiX7eWdleQG3SwAALADCMoCkUpaXpTt/f60e/vwlOq+hxO1ypvX/fbyJNgwAyBCEZQBJaXVtof795gv1T390lmqLst0u5wO/c0aNzlmcnCEeAOA8wjKApHbD2ho99VeX64tXr1C2y9OgZfu9+l8bTnO1BgDAwiIsA0h6Qb9XX7i6UU/fcrk2nlnjWh03X75M1YXJM8oNAEg8wjKAlFFdmK17bjxLv/rzC7W2rnBB911blK3/cfnSBd0nAMB9hGUAKeecxSV66HMX65ufWKvy/KwF2eet153GangAkIEIywBSkjFGn2xepGduWaebL1+mgC9xb2fnNZTo42e41/4BAHAPYRlASsvL8unW607TE1+8TOubKh1/fo+RvvLxJsefFwCQGgjLANLC4tJc/fDTzfrFZ8/Xiso8x573k82LtLp2YfujAQDJg7AMIK1c0limLV+4TLdvPF1FOf64nis/y6dbPrbSocoAAKmIsAwg7Xg9Rp++sEHP3rJOn7lwsXzzXDv7L65arrK8hTmBEACQnAjLANJWUU5At21crc1fuFSXNpbN6bFLy3L1pxcvSVBlAIBUQVgGkPZWVObr5589Xz/8dLMaSnNiesyXr18lv5e3SADIdPwmAJAx1jdV6vEvXq5brztNeVm+Gbe7bEW5rlrl/MwaAIDUQ1gGkFECPo9uvnyZnrllnT7ZXKdT25l9HqOv3LDKneIAAEmHsAwgI5XnZ+mbnzhDD33uEjUvLv7g9j+5YLGWV+S7WBkAIJkQlgFktDV1hXrgzy/SPTeeqdNrCvTF9SvcLgkAkERmbtoDgAyy8cxabTyz1u0yAABJhpFlAAAAYAaEZQAAAGAGhGUAAABgBoRlAAAAYAaEZQAAAGAGhGUAAABgBoRlAAAAYAaEZQAAAGAGhGUAAABgBoRlAAAAYAZxhWVjTIkx5gljzO7Jr8UzbFdkjHnAGLPDGNNqjLkwnv0CAAAACyHekeVbJT1lrW2U9NTk9encI+nX1trTJJ0hqTXO/QIAAAAJF29Y3ijpp5Pf/1TS7566gTGmQNJlkn4sSdbaMWttT5z7BQAAABIu3rBcaa09KkmTXyum2WappBOS/tkY85Yx5kfGmNw49wsAAAAk3Kxh2RjzpDHmvWkuG2Pch0/S2ZK+Z609S9KgZm7XkDHmJmPMVmPM1hMnTsS4CwAAAMB5vtk2sNZePdN9xpjjxphqa+1RY0y1pI5pNmuX1G6tfW3y+gOKEpattfdJuk+Smpub7Wz1AQAAAIkSbxvGw5I+M/n9ZyQ9dOoG1tpjkg4ZY1ZO3nSVpJY49wsAAAAkXLxh+U5J640xuyWtn7wuY0yNMWbzlO3+QtIvjTHbJJ0p6etx7hcAAABIuFnbMKKx1nZpYqT41NuPSNow5frbkprj2RcAAACw0Iy1ydsWbIw5IemA23U4rExSp9tFpDGOb+JxjBOL45t4HOPE4vgmFsc3MRZba8unuyOpw3I6MsZstdYyyp4gHN/E4xgnFsc38TjGicXxTSyO78KLt2cZAAAASFuEZQAAAGAGhOWFd5/bBaQ5jm/icYwTi+ObeBzjxOL4JhbHd4HRswwAAADMgJFlAAAAYAaE5QQwxpQYY54wxuye/Fo8zTYrjTFvT7n0GWP+cvK+vzXGHJ5y34aP7CSDxXJ8J7fbb4x5d/IYbp3r4zNVjK/fRcaYZ4wxrcaY7caYL0y5j9fvDIwx1xpjdhpj2owxt05zvzHGfGfy/m3GmLNjfSxiOr5/PHlctxljXjbGnDHlvmnfL/BbMRzfdcaY3ik/+1+J9bGYEMMx/uspx/c9Y0zYGFMyeR+v4USx1nJx+CLpm5Junfz+VknfmGV7r6RjmpjjT5L+VtItbv87kvUS6/GVtF9SWbz/P5l2ieX4SKqWdPbk9/mSdklqmrzO63f64+qVtEfSUkkBSe+8f8ymbLNB0hZJRtIFkl6L9bGZfonx+F4kqXjy++veP76T16d9v+Ayp+O7TtKm+TyWy9yPk6SPS3p6ynVewwm6MLKcGBsl/XTy+59K+t1Ztr9K0h5rbbotwJIocz2+Tj8+3c16fKy1R621b05+3y+pVVLtQhWYos6T1Gat3WutHZN0vyaO9VQbJf3MTnhVUpExpjrGx2a6WY+RtfZla+3JyauvSqpb4BpTWTyvQV6/sZnrcfpDSf9nQSrLcITlxKi01h6VJkKFpIpZtr9RH33Bf37yo8Kf0CbwEbEeXyvpcWPMG8aYm+bx+Ew1p+NjjGmQdJak16bczOv3o2olHZpyvV0f/QNjpm1ieWymm+sx+qwmRvHfN9P7BSbEenwvNMa8Y4zZYow5fY6PzXQxHydjTI6kayX9asrNvIYTxOd2AanKGPOkpKpp7vryHJ8nIOl3JP2vKTd/T9Idmnjh3yHpLkn/1/wqTU0OHd+LrbVHjDEVkp4wxuyw1j7vTIWpzcHXb54m3qz/0lrbN3lzxr9+Z2Cmue3U6Yhm2iaWx2a6mI+RMeYKTYTlS6bczPtFdLEc3zc10U44MHmuwn9JaozxsZjbcfq4pJestd1TbuM1nCCE5Xmy1l49033GmOPGmGpr7dHJj1A7ojzVdZLetNYen/LcH3xvjPmhpE1O1JxKnDi+1tojk187jDH/qYmPuJ6XNJf/n7TkxPE1xvg1EZR/aa19cMpzZ/zrdwbtkhZNuV4n6UiM2wRieGymi+X4yhizVtKPJF1nre16//Yo7xeYMOvxnfIHs6y1m40x3zXGlMXyWEia23H6yCfSvIYThzaMxHhY0mcmv/+MpIeibPuRnqPJgPK+35P0nqPVpb5Zj68xJtcYk//+95Ku0W+P41z+fzJRLMfXSPqxpFZr7d2n3Mfrd3qvS2o0xiyZ/ETpRk0c66kelvTpyVkxLpDUO9kKE8tjM92sx8gYUy/pQUmfstbumnJ7tPcLTIjl+FZNvjfIGHOeJjJGVyyPhaQYj5MxplDS5Zry3sxrOMHcPsMwHS+SSiU9JWn35NeSydtrJG2esl2OJt5ICk95/M8lvStpmyZ+UKrd/jcl0yWW46uJs4nfmbxsl/Tl2R7PZU7H9xJNfDy4TdLbk5cNk/fx+p352G7QxMwhe95/TUq6WdLNk98bSfdO3v+upOZoj+Uy5+P7I0knp7xmt07ePuP7BZc5Hd/PTx6/dzRxAuVF0R7LZe7HePL6f5d0/ymP4zWcwAsr+AEAAAAzoA0DAAAAmAFhGQAAAJgBYRkAAACYAWEZAAAAmAFhGQAAAJgBYRkAAACYAWEZAAAAmAFhGQAAAJjB/w9qpI6SKs6A2wAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "\n",
    "circle = Circle((0, 0), R, alpha=0.4)\n",
    "ax.add_artist(circle)\n",
    "\n",
    "wheels = []\n",
    "for i in range(3):\n",
    "    px, py, theta = wheel_pose(i)\n",
    "    wheel = Rectangle((-w/2, -r), w, 2*r)\n",
    "    wheel.set_transform(Affine2D().rotate(theta).translate(px, py))\n",
    "    wheels.append(wheel)\n",
    "\n",
    "ax.add_collection(PatchCollection(wheels))\n",
    "ax.axis('equal'); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Kinematics of Omni Wheels\n",
    "\n",
    "> From robot velocity to wheel velocity.\n",
    "\n",
    "Kinematics is the study of motion, without regard for the forces required to cause the motion.\n",
    "In particular, for our omni-directional robot kinematics addresses the question \"What should the individual wheel velocities be given that we want the robot to move with a certain velocity?\"\n",
    "\n",
    "The key property of omni-wheels is that they roll without slipping in the steering direction,\n",
    "but roll freely in the perpendicular direction (sometimes called the slipping direction).\n",
    "In the figure below, these directions correspond to the unit vectors $u_\\|$ (steering direction)\n",
    "and $u_\\perp$ (slipping direction).\n",
    "These vectors define an orthogonal decomposition of the robot's motion.\n",
    "The wheel actuators directly produce motion in the direction $u_\\|$,\n",
    "and the robot moves freely in the direction $u_\\perp$.\n",
    "\n",
    "<figure>\n",
    "<img src=\"https://github.com/gtbook/robotics/blob/main/Figures4/omni-wheel-frame.png?raw=1\" id=\"fig:omni-wheel-frame\" style=\"width:14cm\" alt=\"\">\n",
    "<figcaption>Orthogonal coordinate frame for a single omni-wheel.</figcaption>\n",
    "</figure>\n",
    "\n",
    "The vectors $u_\\|$ and $u_\\perp$\n",
    "are defined in the local coordinate frame of the robot, but since our robot\n",
    "does not rotate, we can define them directly in the world coordinate frame.\n",
    "Let $\\theta$ denote the angle from the world $x$-axis to the axle axis of the wheel.\n",
    "The vectors $u_\\|$ and $u_\\perp$ are given by\n",
    "\n",
    "$$ u_\\| =\\left[ \\begin{array}{c} -\\sin \\theta \\\\ \\cos \\theta \\end{array}\\right]\n",
    "\\;\\;\\;\\;\\;\\;\\; \\rm{and} \\;\\;\\;\\;\\;\\;\\;\n",
    "u_\\perp =\\left[ \\begin{array}{c} \\cos \\theta \\\\ \\sin \\theta \\end{array}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When an object moves with a purely translational motion (i.e., there is no rotation), every point\n",
    "on the object moves with the same velocity.\n",
    "Suppose that the robot, and therefore this specific wheel,\n",
    "moves with  the **translational velocity** $v=(v_x, v_y)$, as shown in the figure below. \n",
    "We can decompose the velocity $v$ into components that are parallel and perpendicular to the\n",
    "steering direction by projecting $v$ onto the vectors $u_\\|$ and $u_\\perp$.\n",
    "\n",
    "<figure>\n",
    "<img src=\"https://github.com/gtbook/robotics/blob/main/Figures4/wheel-velocity.png?raw=1\" id=\"fig:wheel-velocity\" style=\"width:14cm\" alt=\"\">\n",
    "<figcaption>Wheel velocity v.</figcaption>\n",
    "</figure>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In particular, we can write\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "v &= (v \\cdot u_\\| )u_\\| + (v \\cdot u_\\perp) u_\\perp \\\\\n",
    "&=\n",
    "(- v_x \\sin \\theta + v_y \\cos \\theta) u_\\|\n",
    "+\n",
    "(v_x \\cos \\theta + v_y \\sin \\theta) u_\\perp \\\\\n",
    "&= v_{\\parallel}u_\\| +  v_{\\perp}u_\\perp\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "which gives\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "v^i_{\\parallel} &= - v_x \\sin \\theta^i + v_y \\cos \\theta^i \\\\\n",
    "v^i_{\\perp} &= v_x \\cos \\theta^i + v_y \\sin \\theta^i\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Note that in the above equation there is some possible confusion in the notation:\n",
    "$u_\\|$ and $u_\\perp$ are unit vectors that denote directions of motion,\n",
    "while $v_{\\parallel}$ and $v_{\\perp}$ are scalar quantities (the speed in the directions\n",
    "of the unit vectors $u_\\|$ and $u_\\perp$, respectively).\n",
    "\n",
    "In code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wheel_velocity(vx:float, vy:float, i:int):\n",
    "    \"\"\"Calulate parallel and perpendicular velocities at wheel i\"\"\"\n",
    "    _, _, theta = wheel_pose(i)\n",
    "    para = - vx * np.sin(theta) + vy * np.cos(theta)\n",
    "    perp =   vx * np.cos(theta) + vy * np.sin(theta)\n",
    "    return para, perp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we use this to show how the desired velocity (in blue) translates graphically to parallel (purple) and perpendicular (green) velocity components at the wheel sites:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAssAAAHSCAYAAAAaOYYVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABb0klEQVR4nO3dd3hc1Z3/8feZot6rJbnJvWGawfRuiiFxIKSRAmmE7LIpu+xuNtn9JaQSNsmGhISEkBACyZIsIQGD6ZhebcBgy03usmxLlmx1adr5/SEZZFtlpDszd6T5vJ5HjzUzd+79chmNPjrzvecYay0iIiIiInI0j9sFiIiIiIgkK4VlEREREZFBKCyLiIiIiAxCYVlEREREZBAKyyIiIiIig1BYFhEREREZhM/tAoZSUlJip06d6nYZIiIiIjKOrV69er+1tnSgx5I6LE+dOpVVq1a5XYaIiIiIjGPGmB2DPaY2DBERERGRQSgsi4iIiIgMQmFZRERERGQQCssiIiIiIoNQWBYRERERGYTCsoiIiIjIIBSWRUREREQGobAsIiIiIjIIhWURERERkUEoLIuIiIiIDEJhWURERERkEArLIiIiIiKDUFgWERERERmEwrKIiIiIyCAUlkVEREREBqGwLCIiIiIyCIVlEREREZFB+NwuQEQkWVhrObjtIPWr6tn10i72rdnHRx/4KOl56W6XJiIiLlFYFpGUZK2lbXcbu1/fTd0rdex4dgcN7zQAYLyGQHuA0rmlGK9xuVIREXGTwrKIpIT2fe3Ur6pn96u72b5yO/ve3kc4EMab5qWnrQdsv40N5EzI4eqVV5OWneZazSIi4j6FZREZd7qau6hfXc/u13qD8d439xLoCOBL9xHoCGDD7yXjUHfoqOen56Xz6ec+TXZZdiLLFhGRJKSwLCJjWk9bD3ve2MPu13az45kd1K+up6u5C3+mn2BnkEgo8u624Z7wsPvzZ/n5xGOfoGhGUTzLFhGRMSImYdkYczFwC+AF7rDW3jTANucAPwX8wH5r7dmxOLaIpKblX1jO5oc309HQ0RuMu4NEAu8F455gz4j36cv08cF7P8jExRNjWaqIiIxhjsOyMcYL/AJYAtQBrxtjHrTW1vTbpgD4JXCxtXanMabM6XFFJHXZiKX20VradrcBowvGR/Jn+VnyoyXMft9sx/sSEZHxIxbzLJ8M1Fprt1prA8C9wLIjtrkKuN9auxPAWtsQg+OKSIoyHsP7bn8f/ix/TPbnz/Jz0vUncdIXT4rJ/kREZPyIRViuAnb1u13Xd19/s4BCY8wzxpjVxphPxeC4IpLCZlw0g6rFVRiPs6ndfJk+Zr1vFhfcdEGMKhMRkfEkFmF5oN9U9ojbPuBE4FLgIuC/jDGzBtyZMdcaY1YZY1Y1NjbGoDwRGa8u+9VleNO8o36+N91L5YmVXH735Rij+ZRFRORosQjLdcCkfrcnAvUDbPOotbbDWrsfeA44dqCdWWtvt9YustYuKi0tjUF5IjJeFc8q5oTPn4A3Y+SB2ePzUDC1gKtWXIXXP/rALSIi41sswvLrwExjTLUxJg34KPDgEds8AJxpjPEZY7KAxcD6GBxbRFLcud85F1/6CK9VNpBZnMk1z1xDeq6WshYRkcE5DsvW2hBwPfAYvQH4L9badcaY64wx1/Vtsx54FHgbeI3e6eXWOj22iEhGfganfOUUPP7o387Sc9O55tlryJmQE8fKRERkPDDWHtlenDwWLVpkV61a5XYZIpKkOvd3cveSu9n71t6on+PL9PGppz7FpFMnDb+xiIikBGPMamvtooEe0wp+IjLmhINhnv7607z0o5cAWPKjJVQuquRPS/9EsDM46PN8mT6u+OMVCsoiIhI1hWURGVM2PriRe5fdC8CMS2Zw5b1Xkp7X23dcfX41tY/UHrbE9SH+LD/n/+B85l4+N6H1iojI2KawLCJjQtOmJm6dfSvQuyjJ9Zuup2h60WHbLL11KbfOufWosOzP8nPidSey+EuLE1aviIiMDwrLIpLUulu6+csH/8K2p7YBcNXDVzFz6cwBt82fnM8pXzmFV376CqGuENDbejHjkhlc+KMLE1aziIiMHwrLIpKUIuEIL978Ik9//WkAzvzGmZzzrXPw+Iae9eLMb5zJ6l+vJtQVwpvuZcJxE/jgnz6oRUdERGRUFJZFJOlsW7mNP5z3BwCqTq7i4498nMyizKiem5adxiU/v4S/X/138ifl84lHP+FolT8REUltCssikjRadrZw2zG30dPaA8AX3/kiZQvKRryfBR9bQOf+TuZdOe/di/9ERERGQ2FZRFwX7AzywKcfYN1f1gFwxZ+uYMFHF4y6dcIYo4v5REQkJhSWRWTEmjsCHOwMMK3U2Qp41lpW/3o1D3/xYQAWfXERF/3PRSNfvjoKWw9ARyDmuz1KZS4UZ8X/OCIikhgKyyIStVA4wl0v7+CWJzdRXZrD3//htFGP/u5+bTd3LL4DgMJphXzmxc/EbfnpjgBccDdk+iFel/mFItAVgmtPgG+cGaeDiIhIwiksi0hUntnYwHcfXk9tQzsAa3Yd5K9v7ObKEyeOaD/te9v53em/48DWAwB87rXPUXVSVczr7S87DSblwfaW+B0j0wdXLYAbTo3fMUREJPEUlkVkSFsb2/nOQzWs3Nh41GM3P7qBSxZMIDuKtolQT4jHvvoYq25bBcClv7qUE689MWFTun14Pvz0VQiEY7tfj4EMH/xkCVwy8PTPIiIyhiksi8iAWruD/OzJzdz18naCYTvgNg1tPdy6spZ/v3jOoPux1rL23rXcf9X9AMz/8HyW3bkMf5Y/LnUP5n2z4JZXY7vPDB9MzIU7l8Hk/NjuW0REkoPCsogcJhKx/HnVLn78+Eb2tw9/RdxvX9jGx06azOQBrmprWNfAbQtuAyA9P50vvv1F8l1KlZPzoSIndq0YmT64fA5862yIw/WIIiKSJPQWLyLvenVrEzcur6FmT2vUzwmEInxvRQ2//uSid+/rau7inovvof71egA+9fSnqD63Oub1jlQsWjHUdiEikloUlkWEugOd/GDFBh5+Z8+onv/Yun28VLufU6YWsfKbK3nh+y8AcP4Pzue0fz0Nj3foJaoTxWkrhtouRERSj8KySArrCoS57Zlabn9+K93BiKN9/ezm53niN28DUH1BNR/560eSbvU8j4GeUY4qq+1CRCQ16S1fJEX9/c3d/PDRDexp6Xa0n9wD3XzwN2sxfdcAXr/peopnFsegwtioa4U73oQ733rvPr8Hov3bQG0XIiKpTWFZJMW8XXeQG5fXsHrHAUf78fWEOfeBLUzc1tvf/OQVM2g/tpSvVOXFokxHBgrIpVnwg/OhugCW/im6/ajtQkREFJZFUkRDWzc3P7qRv75Rhx14JrjoWMv81/dx8so6AN4+uZw3zqrCej3QGeR/ntzEt94/PzZFj8BQAfnsKZDmfe/+ylzYdnDo/antQkREQGFZZNzrCYX57Qvb+OXKLbT3hBztq3xXG0v/tBGAprJMHvvILHqOmC/5nld28IlTJjOjLNfRsaIxkoDc35Vz4ZbXBp4VQ20XIiLSn8KyyDj22Lq9fH/FenY0dTraT1ZrgGW/ryGjqzds//3T8zhQdvS8ygChiOXG5TXc/dnFjo45mNEG5P7ePxt+9trR96vtQkREjqSwLDIObdrXxreX1/BC7X5H+/EGI5z26HZm1DQD8OylU9k6vxiGWaL6+c37ebJmHxfMK3d0/ENiEZD7m5x/dCuG2i5ERGQg+pUgMo4c7Azwkyc28cdXdxKOOGhMtpaZb+/njEd3ALBxYQmvLplM2Bf9fMnfW7Ges2aVkjaC5/QX64B8pEOtGKGI2i5ERGRwCssi40AoHOGeV3bw06c2c7Az6GhfxXs7eP9d6wFoz/Xz0Cfn0pWbNuL9bNvfwe9f2sa1Z02P+jnxDsj9vX82/PfLMKNQbRciIjI4hWWRMe75zY1856EaNu1rd7SfjI4gl/xpIwXNvfMuP/SJOTRW5Tja5+9f3M5nTq/GN8QKfokMyP1NzoeHPgazitR2ISIig9OvCJExavv+Dr77cA1Prm9wtB9POMKilXXMX927n5cvmMSGE8qG7Usecp8GPnLSJG64cPaAQdmtgHykY8oScxwRERm7FJZFxpj2nhA/f2ozd764nUDYwRLV1jJl4wHOe2ArANtnFvD8ZdWEHCbVk6uL+Ob75jG/8vC+BqcB2Vp72PzQxoBxEOhFRESiobAsMkZEIpb7Vtdx82Mb2d/e42hf+fu7uOK36wAI+j38/bPzac9Pd7TPqoJMvr50LpcurHj3voECcnGm5V8XBzi2NEAoEqY7GGHVtgiBcIRwxBIMRwiGLaFI77/hiB1yERVjwOcx+LwGv9eDz+PB7zX4vB78nt5/M/weMv1eMvxestK8ZKZ5yfB58XgUtkVEZGgKyyJjwKrtzdy4vIZ3drc42k9ad4gL7ttM+e4OAB79yCz2THW2PHWm38sXz5nOtWdNA2DNnm7ufMvD3za9d1Fgrj/MB6qbmVXQhc8Dthve2uXosO+yFoJhSzBs6SL6kXZjIN3nISvNS3a6j9wMP3kZvf/mZvjI8CeoF0RERJKawrJIEqs/2MUPHtnA8jX1jvZjIpZjX6rn+Bf3ALD6zEreOaUC62Bk1QDnzC7lihOqaA34+ezfOnlhTy7QGzJz/GEu7xeQk4210B2M0B2M0NwRBLoOezzN5yE3w0dehp+CLD/F2WkUZqfhH+JiRRERGX8UlkWSUHcwzK+e3cKvn91KV3CANZlHoHJbCxf9ZTMAeybl8PQVMwhkOPvRn1iYydlzJrG1u5zrn3lvZDrZA/JIBEIRmtoDNLUHDrs/L9NHUVZvcFaAFhEZ/xSWRZLMg2vq+eEjG9h9sGv4jYeQc7CHK+5Yizfc2/B7/+fm01Kc6Wif2Zk5TKicR21nKfds7zvOOArI0WjtCtHaFWJ73xLixkBhlp+yvAzKctMpy80Y9UIsIiKSfBSWRZLE2t0t3Lh8Ha9vP+BoP75AmLOXb2VybW9/81OXT2fnrMJR7896MzHZ1djsatqB2s7UC8hDsRaaO4I0dwTZsKftsPBcnpdBaU66wrOIyBimsCzissa2Hn702Eb+b/UunKxQjbXMXd3AKU/1Xjm3dlEZq86ZiB1Fi4D1ZkJWNeRU995GATlaR4Znj4HS3HQmFmZRVZhJjlZAEREZU/SuLeKSQCjCnS9u49ana2nrCTnaV1ldO5f+cQMAB0oyeORjs+nJ8o9oH0cGZIBMb4grpx9QQHYgYmFfaw/7WntYveMA+Zl+qgozmViYSUmOs+n6REQk/hSWRVzwZM0+vrdiPdv2dzjaT2ZbgPfftZ6sjiAAD1wzj+byrKifP1BAJtLDyYV1XDrLR5pX8xDHWktXkJauIDX1rWT4PUwqymJKcRZluRlulyYiIgNQWBZJoNqGNr790Hqe29ToaD/eUIRTH9/BzHeaAHh+6VRqFxRHtUT1gAE53I2n9R1Orohw4bwyMtNGNioto9MdjLB5Xzub97WTne5lanE2U0uyyc/U+RcRSRYKyyIJ0NIZ5H+e3MQ9r+wg5KQx2VpmrG3izBXbAdg8v5iXL5pC2D90j8RgAZmWd6CnkVll2Sw9vYLyPI1uuqWjJ8y6+lbW1bdSlO1nSnE2U4uzyXS4/LiIiDijsCwSR+GI5U+v7uAnT2ziQGfQ0b6K9nWy7Pc1AHRm+1l+9Vw6c9MG3X64gGywFGensfT4ycytcLaKn8RW7wWCB3lr10GqCjKZWZ5DRb6zaf9ERGR0FJZF4uSl2v18+6EaNuxtc7Sf9M4gF9+7iaLG3nmXH/74bBom5g64bTQBGXqXeT53djmnzSjG59GVe8nKWqg70EXdgS5yM3zMKMthWmk26T6NNouIJIrCskiM7Wzq5Hsranhs3T5H+zHhCCc+W8cxrzcA8Op5E6lZVH5UX3K0ARl6l6g+YUohF84rJzdDfbFjSVt3iDd3HuSduhYmF2cxqzyXouzBP1kQEZHYUFgWiZGOnhC3rqzlty9sIxCKONrX5E0HOP9vWwDYOT2fZ98/jVC/3tWRBORDphRlcdnCSqoK9XH+WBaKWLY2drC1sYPS3HTmVeZRVaD/pyIi8aKwLOKQtZa/vrGbmx/dQENbj6N95TV188E71gIQ9hru/9wC2gt65+IdTUAGyM/0c/GCCRw7scBRbZJ8Gtt6eHZjIwVZfuZV5DGlOAsTxYwoIiISPYVlEQfe2HmAG5fXsGbXQUf78feEOP/+Wip2tgPw2IdmUj8tf9QBGcDvNZw5s5SzZ5XiH8UqfjJ2HOwM8tKWJt7e3cK8ilyqS3LwehSaRURiQWFZZBT2tnRz0yPreWBNPdbBTHAmYjnmlT2c+Hw9AG+eVsFbZ03H5kwbVUA+ZOHEfC6eP4GCLPW0ppL27hCvbTvAO7tbmDMhj1nluQrNIiIOKSyLjEB3MMxvntvKbc9uoTMQdrSviu2tXPznTQDsq8rliWuWEiid894GIwzIAJUFGVx2TCVTS7Id1SZjW1cgwps7D7JpXxsLqvKZVpKt9gwRkVFSWBaJ0sNv7+EHj6yn7kCXo/1kt/Twgd+tIy3QexHgX7/2D7RMKO19cBQBGSAn3ceF88o5YUohHoUi6dPRE+bVrc2s39PKsRMLmFQU/VLoIiLSS2FZZBg19a3cuHwdr25rdrQfXyDMGSt2Ur2xd4nqp6++ku3HzYNIDzS/PuKADOD1GE6bXsy5s8vI8GvuXRlYa1eI5zfvpzgnjeMmFWilRhGREVBYFhlEU3sPP3p8E39+fSdOVqgGKO9qZenPelsu1p15Mq9fdhaRjhrY+8iIA/IhcybkcukxFRTnpDsrTlJGU3uAp9Y3UFmQwQlTCsnTXNsiIsNSWBY5QjAc4a6XtnPLU5tp6w7FZJ/FhY1kXbGdh4+7gG0dFtP8DKNtlijLTefSYyqYWT7wKn4iw6k/2M3elj3MrchjfmUePs2WIiIyKIVlkX5WbmjgOw/XsLWxI6b7remexuXfuIuP8zbf+8u3gJGHk0y/l/PnlrG4ulgzHIhjEQvr6lvZ3tTBCZML1c8sIjIIDSeIAFsa27nmztf49O9fj3lQ7mX47RPXAbBoxqsjeqbHwOLqIv5lySxOm16ioCwx1dET5vnN+1m5oYGWrqDb5YiIJB2NLEtKa+0O8tMnNnP3K9sJhh02Jg9j74FK9reWcNEJj7Bm+wkEQ8P3Gk8vzebShZVM0AVZEmd7Wrp55J3e1owFVfn6o0xEpI9GliUlRSKWP766g3P/+xl+9+K2uAflQ+5e+RkA3n/yX4fcrig7jU8snsxnz5imoCwJc6g149G1e2lqd7Z0u4jIeKGRZUk5r2xt4sblNazf05rwY3f25PDapsWcPOtV8rMO0NJZeNjjaT4P584q5fQZJbroSlzT0hXk8Zp9zK3IY2FVPh6NMotICtNvY0kZu5o7+eI9q/no7a+4EpQPeertiwD4wsU/f/c+A5wwuYB/XjKLs2eXKSiL66ztnWP8EY0yi0iKi8lvZGPMxcaYjcaYWmPM14bY7iRjTNgYc2UsjisSjc5AiB89tpELfvIsj6zd63Y5RCI+/vrSh/H7QlSX1zK5KIsvnjOdK0+cpHlvJekcGmV+a9dBIk4nHBcRGYMct2EYY7zAL4AlQB3wujHmQWttzQDb/RB4zOkxRaJhreXvb+3mh49sZG9rt9vlHGZD3XwArjr7D2QF78Og1fckeR0aZd7b0sXpM0rI1R91IpJCYjGyfDJQa63daq0NAPcCywbY7p+AvwINMTimyJDe2nWQK257ia/+eU3SBWW/13Du7HJ83bcAEPTc53JFItFp7gjyyNq9bN8fj+kVRUSSUywu8KsCdvW7XQcs7r+BMaYKuBw4DzgpBscUGVBDazc3PbqBv725G5uEnxgvqMrnkgUTKMxKAyASmUPQey/+yGUYsl2uTmR4obDlpS1N7Gvt5sQpheqvF5FxLxZheaDLpI+MKT8F/t1aGzZm6KuqjTHXAtcCTJ48OQblSSroCYW54/lt/HJlLR2BsNvlHKUiP4PLFlZSXXJ4IM4If4NOzyfp9n6fzPD3XKpOZOS2NHawvz3AGTNKyM9SW4aIjF+xCMt1wKR+tycC9Udsswi4ty8olwBLjTEha+3fj9yZtfZ24HaARYsWJeHYoCSbR9fu4fsrNrCzudPtUo6SneblwnkTOHFqIZ4B/lA05OAPf4Sg989EwnV4mOhClSKj09IV5LF1ezlxaiHTS3PcLkdEJC5iEZZfB2YaY6qB3cBHgav6b2CtrT70vTHm98BDAwVlkZHYsLeVby+v4aUtTW6XchSvMZw6vZjz5pSR4R/64j1/5EMEvX+my/9PZAf/Nuh2liAhswqfPQHD8Kv/iSRCKGJ5dWszzR0BTpxcqDmZRWTccRyWrbUhY8z19M5y4QV+Z61dZ4y5ru/xXzk9hkh/BzoC/PiJjfzva7sIJ+FUVrPLc7n0mApKcqMLtAYv6aFv0uO7kZB5GZ899bDHLZaQeZ2g91dYWjDhr+CzZ8ajdJFR27yvndauIKfPKBn2D0QRkbHE2GS8CqrPokWL7KpVq9wuQ5JEKBzhDy/v4JanNtPSFXS7nKOU5qRz6cIKZpXnjur5Hb4PgwmSFfwLht4e0Ai76PH+kojZBqYHLPgil5Ie+VwsSxeJmex0L2fPKqWg7yJWEZGxwBiz2lq7aKDHtNy1jAnPbWrk2w/VUNvQ7nYpR8nwezh/TjmnTCvG6+Aj6MzQz+nyX0fA8wfSIh8h4LmbkOcZIAim749aA2HPOojEpHSRmOvoCfN4zT5OnVbMpKIst8sREXFMYVmS2rb9HXz3oRqe2pB803N7DCyaUsSSeeVkpzv/UfJQjjdyEiHvQ4Q8T9EbkkNHbWfZjSWshUwkaYXCluc372fhxHwWVOW7XY6IiCMKy5KU2rqD/Oypzdz10g4C4eQbRq0uyeayhRVU5GfGbJ9h8zYRswOsAdM1xJZeLPWYwyahEUk+b9e10NET4uTqIoabNlREJFkpLEtSiUQsf1m1ix89vpH97QG3yzlKYZafSxZUxHS0LMIeery/JmI29PYlRyFsavFYhWVJflsaO+gORTh9erEWMBGRMUlhWZLGa9uauXH5OtbVt7pdylHSvB7Onl3KGTNK8MfoF76li4DnXkKeR4EQmChH0E03EbMe7LkxqUMk3nYf6GLlxkbOmlVCuk/tQyIytigsi+t2H+zi+yvW8/Dbe9wu5SgGOG5SARfNn0BeZuxWKQua5wh4b6e3L3nkI+hhz3pd5CdjSmNbD0/WNHDunFKy0vSrR0TGDr1jiWu6AmFue3YLtz+3he5g8iW/SYWZXLawMi5X9Ac9DwBd0Y8mH8GyB0sIox9hGUNauoI8UbOPc2aXkR/DPz5FROJJv2nFFQ+8tZsfPrKB+pZut0s5Sl6GjwvnT+D4SQVxuygpM/xNurzfxFI/qpFl8BNhF16qh99UJIl09IR5omYf584upThHK1GKSPJTWJaE++5DNdzxwja3yziKz2M4Y0YJ58wuI80X3wuRDHlkhm+i23sTEdZHfWHfeywRU4vXKizL2BMIRXh6QwPnzimjRIFZRJKcLk2WhLty0URHi3fEw/zKPL5ywSwunD8h7kH5EEM6GeH/xBc5E+wIA4PpIWJq4lOYSAIEw5aVGxpobBvpH4oiIomlsCwJN2dCHh87OTmmPZuQl8Hnzqjm44unUJSd+OV5DV7SIv+AP/JBsCM7ftizMU5ViSRGMGxZuVGBWUSSm8KyuOJflsx29QKfrDQvy46r5PrzZjCtNMe1OgAMhrTIh0gLXzeiwGxpwBKMY2Ui8RcKW57Z2MD+dgVmEUlOCsviisLsNL58/syEH9dj4LTpxfzLktksri7Gk0SrivntuaSHvw42I8pnpBEx2+NZkkhCHGrJaFJgFpEkpLAsrvnUqVOYUZa4Ud1Z5Tl86byZXLawksy05FwYwWePJTP0A7A5YIf78QwTNpsTUpdIvPW2ZDTS0qlPS0QkuSgsi2t8Xg//77J5cT9OSU4anzp1CtecVk1ZXrSjtu7xMJXM0E8wFIMdYsIaE9BFfjKuBEIRVm5soKMn5HYpIiLvUlgWV501q5Tz55TFZd8Zfg+XLJjAl8+fxZwJeXE5Rrx4KCUz9BM8dsqQfcwRsymBVYnEX2cgzDMbG+kJhd0uRUQEUFiWJPCfl80jzRu7l6IBTppayD8vmc2ZM0uTbpq6aBlyyAh/H689ZtCp5SzNWJJvYRcRJ1q6gjy7sZFQOPlW9hSR1KOwLK6rLsnmmtOnxmRfU4uz+cdzZ3D58RPJSR/7a+4Y0kgPfx1f5NxBAnMaEZN8C7yIOLW/PcALtfuJRKzbpYhIilNYlqTwT+fNoCRn9PMcF2T5+ehJk7j2rGlUFmTGsDL3GTykR76AP/KRAVoygoRNrSt1icRb/cFuXtve7HYZIpLiFJYlKeRm+Lnhwtkjfl66z8MFc8v46gWzWDixIPaFJZG0yOWkha8/PDCbEBGz1r2iROJsa2MH79S1uF2GiKQwhWVJGh9eNIkFVdFfiHf+nDK+csEszptTjj+GPc/JzG/PJCP8X71zMfd9Oh3RyLKMc+/sbmFXc6fbZYhIikqNhCFjgsdj+Ob75g+73cKJ+dxx9SIunF/u6iqAbvHaBWSGfgjkgvVgacHS4XZZInH18pYmDnYG3C5DRFKQwrIklZOmFnHZwooBHyvNTee/r1zI/35+MW3dQVL5QnkPk/vmYi4FEyZitrpdkkhchSKWZzc10h3UlHIiklgKy5J0vr50Lhn+916aaT4P1509nZU3nMPlx1fx/OYmugIpnJT7eCghM/RjvJGTgIGnlhMZTzp6wrywWTNkiEhiKSxL0qksyOQLZ00HYMm8cp746ll87ZI55KT7eG1bM80d+ij2EEM2GeGv47Wz3C5FJCEa2npYteOA22WISAoZ+xPRyrh03dnTWVxdxGkzSt69r6a+le1NushHJNXVNrRTlO1nRlmu26WISArQyLIkpcw072FBubGth7frDrpXkIgkldU7DuiCPxFJCIVlSXo9oTAvbdmP2hRF5JBwBF6o3a8lsUUk7hSWJem9urWZjh5dAS8ih2vtCvH6dvUvi0h8KSxLUtu0r426A11ulyEiSWrb/g62Nra7XYaIjGMKy5K0mjsCvLlTo0YiMrRV2w/Q0hV0uwwRGacUliUpBcMRXqzdn9ILj4hIdEIRy0u1+wnrwgYRiQOFZUlKb+48SFt3yO0yRGSMONAZZI1mzBGROFBYlqSzp6WL2gb1IIrIyGzc20ZjW4/bZYjIOKOwLEklGI7w2rZmt8sQkTHIWnh1W5PaMUQkphSWJam8ufOgpokTkVFr7QqpHUNEYkphWZKG2i9EJBbUjiEisaSwLEkhEFL7hYjEhtoxRCSWFJYlKby1S+0XIhI7rV0h3lY7hojEgMKyuG5/e4/aL0Qk5jbubaOlU4uViIgzCsviKmstq7ar/UJEYi9iYdUOvb+IiDMKy+Kq2oZ2mjs08iMi8bGvtYft+zvcLkNExjCFZXFNdzDMmroWt8sQkXHuzV0HCIYjbpchImOUwrK4Zs2ugwRC+gUmIvHVFYjwzm79YS4io6OwLK5oau9hqz4aFZEE2bS3jZYutXyJyMgpLIsrVu04gNUUqCKSIBELq3Wxn4iMgsKyJNzOpk6a2gNulyEiKWZvSw97WrrcLkNExhiFZUmoSMSyRgsFiIhL1uw66HYJIjLGKCxLQm3d305bd8jtMkQkRTV3BNnRpOslRCR6CsuSMKGwrkgXEfetqWshEtFFEyISHYVlSZiN+9roCmiqOBFxV3t3iC2N7W6XISJjhMKyJERPKExNfavbZYiIALC2voWQFioRkSgoLEtCrN/TRjCsjz1FJDl0BSJs2NvmdhkiMgYoLEvc9YTCbNqnX0oiklw27m3T6LKIDEthWeJu8752QhpVFpEk0xOKsKVRM2OIyNAUliWuQuEIG/VRp4gkqQ17WzUzhogMSWFZ4qq2sZ2ekD7mFJHk1NETZpvmXRaRIcQkLBtjLjbGbDTG1BpjvjbA4x83xrzd9/WSMebYWBxXklskYtmwR6PKIpLc1u9pxVqNLovIwByHZWOMF/gFcAkwD/iYMWbeEZttA8621i4EvgPc7vS4kvy2NXXQGQi7XYaIyJBau0LUHehyuwwRSVKxGFk+Gai11m611gaAe4Fl/Tew1r5krT3Qd/MVYGIMjitJzFrL+j2aV1lExoZ19VpdVEQGFouwXAXs6ne7ru++wXwWeCQGx5UkVt/STWtXyO0yRESi0twRpKGt2+0yRCQJ+WKwDzPAfQM2fxljzqU3LJ8x6M6MuRa4FmDy5MkxKE/csFnzKovIGFO7r52y3Ix3b0/92sMx3f/2my6N6f5EJDFiMbJcB0zqd3siUH/kRsaYhcAdwDJrbdNgO7PW3m6tXWStXVRaWhqD8iTR2ntC7GnRCI2IjC07mzvpDuo6CxE5XCzC8uvATGNMtTEmDfgo8GD/DYwxk4H7gU9aazfF4JiSxDbva0MXlovIWBOxUNvQ7nYZIpJkHIdla20IuB54DFgP/MVau84Yc50x5rq+zf4fUAz80hjzljFmldPjSnIKRyxbtSKWiIxRWxrbNY2ciBwmFj3LWGtXACuOuO9X/b7/HPC5WBxLktuOpg4tQiIiY1ZHT5i6A11MKspyuxQRSRJawU9iarM+whSRMU6tGCLSn8KyxMzBzgBN7QG3yxARcWRPSzftPZr6UkR6KSxLzGzdr15lERkftkf5fpbeGSSrtSfO1YiImxSWJSastexs6nS7DBGRmNjeFF1YPuWJnVz0581xrkZE3KSwLDHR0NZDZ0Dzk4rI+BDNCqSecIRJW1rIPdhDTotGl0XGK4VliYltasEQkRRTsb0Na8AamLq+2e1yRCROFJbFsXDEsqtZLRgiklpmrN2PPxDBF7bMfGfQhWlFZIxTWBbHdh/oIhjWJP4ikjoOtWCYvts5LWrFEBmvFJbFsWgvhBERGS8OtWD0p1YMkfFJYVkcCYUj7GnpcrsMEZGEOtSCcYhaMUTGL4VlcWRPSzdhrW4tIinkyBaMQ9SKITI+KSyLI3UHNKosIqlloBaMQ9SKITL+KCzLqFlrqT+osCwiqeXIFoxD1IohMj4pLMuoNbb30BNSD4aIpI7BWjAOUSuGyPijsCyjtlstGCKSYoZqwThErRgi44vCsozabrVgiEiKGawF4xC1YoiMPwrLMipt3UFau0JulyEikjDDtWAcolYMkfFFYVlGZV9rt9sliIgkVDQtGIeoFUNk/FBYllFpaNWoiYikluFaMA5RK4bI+OJzuwAZm/a1aWRZJBW9vDeHzS0ZcT9OSUaQpVNa4n6caEXbgnHIoVaM9vz0uNaV6rpD8O9PQkcw/se6eiGcOSX+x5Hko7AsI9bSFaQritEVERl/ag5kUtuSGeejWE4oSa73mJG0YBwydX0za0+piE9BAkBHAB7cCPF+taR7e8OypCa1YciINWpUWSRlnVjaQZonXtHE4vdE+ND0Jj40I7l6fqNtwThErRiJUZwFx5THb/8ZPphRCE9+UqPKqUxhWUZsn/qVRVLW3MIuIiMdYo2Cz0QoTA9x/TF7OaG0M+b7d2KkLRiHaFaMxPjIfMiMw+fkGT64Yg6suAom58d+/zJ2KCzLiDVoZFkkZaV7LVPzYvse4PdEWFDcyVcW7qUsM/mmpBxNC8YhmhUj/i6eDuEYftjhMZDlh59eCD84H9LVsJryFJZlRFq71a8skupi14rR23bxgepmPjKjmTSvjcE+Y2+kLRiHqBUjMYqzYG5pbPaV4YNpBfDYx+GSmbHZp4x9CssyIs3tAbdLEBGXxaIVI5nbLvobbQvGIWrFSIxYtGKo7UIGo7AsI9LcqbAskuqctmIke9tFf05aMA5RK0b8OWnFUNuFDEcvCRkRjSyLpLYD3V6e35PXN32chRGNuVr8HssHqpuTejS5v7Ld7XhDlkDawGNLvmAEj2XQxz1hy4RdbZpCLk7CEXhlN3zzGRhNh2CGDybmwp3LNJosg1NYlhE5oJFlkZRzKCC/vC/33fvy/CHaQ14iUbYZ+0yE3LQw18xpTPrR5P7ePrWCbXOLBn384ns3ktkZ4uFPzB10m85cfzxKS1n9A/LmfoP2M4ugrhW6onx5HWq7+NbZGk2WoenlIVFr7Q4SDCfnBTgiEluDBeQPTGtmZn43Pg/8dn1pVAuU+D0R5hd1cnn1gaS9iG8wYb+Hg6WD/zf2ZPjI7AwNuY04N1hAvnw2XH8yzCiCpk445XfD78tjeoPyT5boIj6JjsKyRE0tGCLjWzQBub8TSzvY2ZZOIDLY5S9jr+1Ckkc0Abm/4iyYWwJr9g2+T7VdyGgoLEvUdHGfyPgz0oDc39zCLv46yNVvY7XtQtw10oB8pI/Mh01NA7diqO1CRksvF4naQYVlkXHBSUDu79CsGEe2YozltgtJPKcBub+Lp8O3nj38PrVdiFMKyxK11mivmhCRpBOrgHykw1sx1HYh0YllQO7vyFYMtV1ILCgsS1RC4QidgbDbZYjICMQrIPd3qBVDbRcynHgF5CN9ZD6s3987oqy2C4kFvXwkKm3d+uUnMhYkIiD3l+61LCprpyds+IDaLuQIiQrI/V08He55B750ktouJDYUliUqCssiySvRAflIy6oPxPcAMqa4EZD7K86CR66K7zEktSgsS1Rau4NulyAi/bgdkEX6czsgi8STwrJERWFZxH0KyJJMFJAlVSgsS1TUhiHiDgVkSSYKyJKKFJYlKu0KyyIJo4AsyUQBWVKdwrIMKxyx9IQibpchMq4pIEsyUUAWeY/CsgyrO6j5lUXiQQFZkokCssjAFJZlWFqMRCR2FJAlmSggiwxPYVmGpZFlEWcUkCWZKCCLjIzCsgyrS2FZZMQUkCWZKCCLjJ7CsgyrS20YIlFRQJZkooAsEhsKyzIsjSyLDE4BWZKJArJI7Cksy7DUsyxyOAVkSSYKyCLxpbAswwqGrdsliLhOAVmSiQKySOIoLMuwQmEtSCKpSQFZkokCsog7FJZlWMGIRpYldSggSzJRQI5exIaxNoLX43e7FBlnFJZlWBpZlvFOAVmSiQLy6Ly48wdsbl7BlPyzmFX8PipzT1ZwlphQWJZhhdSzLOOQArIkEwVk5zqDTYRtD1sPPsGu1pewhJmcd6aCszimsCxDstYSUhuGjBMKyJJMFJDjJxjpAFBwlphQWJYhaSYMGesUkCWZKCAnnoKzOKWwLEMKa1RZxiAFZEkmFiCtmAvuVkB2m4KzjEZMwrIx5mLgFsAL3GGtvemIx03f40uBTuAaa+0bsTi2xFfEKizL2KCALMnkUEAmbz74e1+Tm5sVkJOJgrNEy3FYNsZ4gV8AS4A64HVjzIPW2pp+m10CzOz7Wgzc1vevJDlFZUlmCsiSTAYKyAB01kF7LTu+d45LlclwFJxlKLEYWT4ZqLXWbgUwxtwLLAP6h+VlwB+stRZ4xRhTYIypsNbuicHxJY6sRpYlySggSzLKmHAs+HP6boWhqwY6tmHCnb13eaE7dNCt8lJCxAZjsh8FZzlSLMJyFbCr3+06jh41HmibKkBhWUSiFrFw81tVgAKyJIfuGR68Ozv5l7OWDbvt3W/fNOw2klz6B+etB59gYu5pXDLz5y5XJYkWi7BsBrjvyOHIaLbp3dCYa4FrASZPnuysMhEZVzwGvnbCbrJ9YQVkSQqPnjOLyaV5sPKaYbf932tPiX9BKez13b+gofOdmO7TY9IwGLLTyphV9D6mF10U0/3L2BCLsFwHTOp3eyJQP4ptALDW3g7cDrBo0SL1ALis99pMkeSRnxZ2uwSRd1njY8f+iVFtW5l7UpyrSW0ZvsKY7OfogHwheemThn+ijFuxCMuvAzONMdXAbuCjwFVHbPMgcH1fP/NioEX9ymODorKIiIx3CsgyFMdh2VobMsZcDzxG79Rxv7PWrjPGXNf3+K+AFfROG1dL79Rxn3Z6XEkMDSyLiMh4pIAs0YrJPMvW2hX0BuL+9/2q3/cW+MdYHEsSy+tRWhYRkfFBAVlGQyv4yZD8Hl1FJSIiY5cCsjilsCxD8ngMXg+EI25XIiIiEh0FZIklhWUZls/jIRxRWhYRkeTVu1iIIS99ogKyxJTCsgzL5zX0hNyuQkREZHBnTfkvTq76JwVkiTmFZRmW3+sBNLetiIgkrzRvLmneXLfLkHFIV2/JsHyaEUNERERSlMKyDKt3ZFlEREQk9SgFybDS/XqZiIiISGpSCpJhZfq9bpcgIiIi4gqFZRlWZprCsoiIiKQmhWUZlkaWRUREJFUpLMuwNLIsIiIiqUphWYalkWURERFJVQrLMiyFZREREUlVCssyLJ/Xg9+rhUlEREQk9SgsS1Ry0rUyuoiIiKQehWWJSm6G3+0SRERERBJOYVmikpuhkWURERFJPQrLEpW8TI0si4iISOpRWJaoaGRZREREUpHCskRFYVlERERSkcKyRCXd5yXdp5eLiIiIpBalH4ma+pZFREQk1SgsS9QKshSWRUREJLUoLEvUCrPS3C5BREREJKEUliVqxdkKyyIiIpJaFJYlavmZfrx6xYiIiEgKUfSRqHk8hgK1YoiIiEgKUViWESlSK4aIiIikEIVlGRFd5CciIiKpRGFZRkQX+YmIiEgqUViWESnI8uP3GrfLEBEREUkIhWUZEWMMpbnpbpchIiIikhAKyzJi5XkZbpcgIiIikhAKyzJiCssiIiKSKhSWZcQK1bcsIiIiKUJhWUbMGEOZRpdFREQkBSgsy6iU6SI/ERERSQEKyzIqEzSyLCIiIilAYVlGpTA7jaw0r9tliIiIiMSVwrKMWlVhptsliIiIiMSVwrKMWlWBwrKIiIiMbwrLMmrleRn4NIWciIiIjGMKyzJqXo/RhX4iIiIyriksiyMT1bcsIiIi45jCsjhSWZCJUSeGiIiIjFMKy+JIht9LaY4WKBEREZHxSWFZHJtaku12CZJglojbJYiIiCSEwrI4NrkoC69eSSnBYunx3EW391/cLkVERCQhFHHEsTSfh4p8Xeg33llC9Hh/TMjzAOB3uxwREZGEUFiWmKhWK8a4Zumi2/tfhM3rAHjsApcrEhERSQyFZYmJyoJM/FqgZFyK0EyX7wYiphZMAMjAa2e5XZZIXGWm6dejiPTSu4HEhNdjmFyU5XYZEmMR6ujy/TOWfWBCffeG8NgZrtYlEm9TivVpmYj0UliWmFErxvgSNuvp8v0b0Aom3O8RH4Zit8oSSYipCssi0kdhWWKmLC+D3Ayf22VIDATNy3R7vwWmC4w97DGPrcaglhsZv4qy0yjKTnO7DBFJEgrLElMzy3PcLkEcCniWE/D+tK8/+QjWg1cX98k4p/cxEelPYVlialpJDj6PRh3HIkuEHs/vCHr+OHBQBiADj52Z0LpEEinN52GKrr8QkX4chWVjTJEx5gljzOa+fwsH2GaSMWalMWa9MWadMebLTo4pyS3N52FysX7RjDWWID3e/ybkeRxMzxBbBvDq4j4Zx6pLsvFplSUR6cfpO8LXgKestTOBp/puHykE/Iu1di5wCvCPxph5Do8rSWxWea7bJcgIWDro9v4nYfPGMEEZIANDQSLKEnGFWjBE5EhOw/Iy4K6+7+8CPnDkBtbaPdbaN/q+bwPWA1UOjytJrCg7jeIcXRwzFkRo6ptDeesQrRfv8dhpCahKxB0V+RnkZWh1ShE5nNOpC8qttXugNxQbY8qG2tgYMxU4Hnh1iG2uBa4FmDx5ssPyxC0zy3Joam92uwwZQoSddPn+C2gHExn+Cdari/tkXJtRdvio8vabLnWpEhFJJsOOLBtjnjTGrB3ga9lIDmSMyQH+CnzFWts62HbW2tuttYustYtKS0tHcghJIlOKs8nwq+8vWYXNWrp8/07vHMpRBGUA0vHq4j4Zp7LTvUwszHS7DBFJQsOOLFtrLxjsMWPMPmNMRd+ocgXQMMh2fnqD8h+ttfePuloZM7wew+wJuazZ1eJ2KXKEoHmegPfWqNouDhfQyn0ybs2ryMMYzeQjIkdzOvT3IHB13/dXAw8cuYHpfff5LbDeWvsTh8eTMWRmWS5+r375JJOA52+jDMoA2Rh08ZOMP5lpHqaV6rUtIgNzGpZvApYYYzYDS/puY4ypNMas6NvmdOCTwHnGmLf6vpY6PK6MAWk+j2bGSBK9cyj/mqDnz6MMymjKOBm3Zpfn4dX88CIyCEcX+Flrm4DzB7i/Hlja9/0LoLVxU9XsCbls3NtGKGKH31jiJmxWE/I8dtTS1VGzPjy6uE/GIb/XaLo4ERmSrsCSuMrwe5lWmu12GSnPa4/DH/kE2HSwo/kbOU0jyzIuzSrPxa9FSERkCHqHkLibW5GHPuF0l8FPWuQKskK/wmtPB5sGIxpk7sFjp8erPBFX+PouRBYRGYrCssRddrqPKcUaXU4GhgIywl8hI/QDjJ3WO9Ic5fMMmlZLxpfpZdlk+L1ulyEiSU5hWRLimIn5Gl1OIl6mkRn+EWnh63vvsEP/z/FofmUZZ3xew/zKfLfLEJExQGFZEiIn3aeLaJKMweC11QB4bHVfa8YAbwnWj8fOT3B1IvE1Z0KuRpVFJCoKy5Iw8yvz8Wne5aTS5e8dWc4I30xm6Bd47aLe0HwYny7uk3El3edhbkWe22WIyBihsCwJk+H3MneCfkEli5B5FYD00H9i8OKhhIzwf5AR/hbGVoHN6Nsy0DvyLDJOzK/K0wwYIhI1vVtIQs2pyCXdp5ed2ywhenw3gfXgsyce9pjXziUz9DPSwp8Fm4WhDEN0FwKKJLvsdC8zyzQDhohEz9GiJCIj5fd6WFCVz+odB9wuJaUFPPcAkBm6dcDHDR789gJ8odMY4RxzIkntmKp8rdYnIiOiIT5JuJllOWSn68Iat1haCHkfwBs5AQ8VQ25ryMKgaf9kfCjI8lNdoteziIyMwrIknMdjOGFyodtlpKwu3zcBSA/f4HIlIol14pRCjNGosoiMjMKyuGJSURYVBRnDbygxFWYb1uzAH75Gi4xISplSnEV5nt5zRGTkFJbFNYumFKIL0hPHYun2/zMA/sj7XK5GJHF8XsPxkwvcLkNExihFFXFNboafOZpKLmFCZiUAGaHvYfSjLylkQWU+WWm6nl1ERke/McVV8yvzdLFfAlgCBHw/B5uP185zuxyRhMnL9DFngqaKE5HRU1gWV/m8Hl3slwA93l8DkBn6kcuViCTWSVOL8GiqOBFxQGFZXKeL/eIrQjNhz9N4I2fjocTtckQSRhf1iUgsKCxLUjh5ahF+r0Z/4qHb968ApIe/6HIlIomT7vNw4hR9aiUiziksS1LITvdxvNoxYi5sNmBNM2mhf9SS1ZJSTppaRIZf10OIiHMKy5I0ZpTlUJGvj0xjxRKh2/cfAPjs+S5XI5I4k4uymFyc5XYZIjJOKCxLUjm5Wu0YsRL0rAAgI/gjDDqnkhrSfR4WTdWnVCISOwrLklTUjhEbli6C3t9i7ES8THe7HJGEUfuFiMSawrIkHbVjONfj/SkAmaHvuFuISAKp/UJE4kFhWZLS4mlFpPv08hyNCPsIe17DF74UQ4Hb5YgkRGaa2i9EJD6URiQpZaX5WDytyO0yxqQu3z8BkBa52uVKRBLDGDhteonaL0QkLhSWJWlNLMxitpapHZGQeQtMkPTQv2Hwu12OSEIsqMzX4iMiEjcKy5LUjp9UQFF2mttljAmWMD2+GwHw2VNdrkYkMcrz0llQled2GSIyjiksS1LzeAynzyjWdHJRCHr+D4DM4M9crkQkMdJ9Hk6bXoIxen8QkfhRWJakl5vhZ3F1sdtlJDVLO0Hvn/FE5uJhktvliCTEqdOLyUxTn7KIxJfP7QJEBhKJWPa2dlNZkAnA5OIsZrTmUNvQ7nJlyanb+30AMsJfd7mS5NUdMoRs/Ecg070R/BqGiLu5Fbnvvj+IiMSTwrIkpXtf38W3H1rHtWdO44vnzCAzzcuiKYUc7Aywvz3gdnlJJcIuIp71+MMfwZDjdjlJqSds+N7qiYCN2zEshrA1nDGhlUunHozbcQQq8jM4blKB22WISIrQ+IckndbuID9+fCPdwQg/e7qW8378DH9/czcej+GsWaVkp+tj1/66/F8CwB/5kMuVJK80jyXLHyZkPXH78hjLjLwuzpvY4vZ/7riWm+HjtBnF6lMWkYRRWJak87MnN9PU8d7o8Z6Wbr7y57f44G0vsWlfG2fOLMXn0S9KgJB5GYCM0Dcx6I+IwRgDx5d04DXxGVn2eyJcMLGFz8xtJNMXv9HrVOf39v7BnO7Ta11EEkdhWZLK1sZ27np5+4CPrd5xgGW/eJHvr1jPjLLsxBaWhCxBenw3g/Xjtce5XU7SO66kE0+Mw7LXWLJ9YT4/r4GzKtvQYGf8GAOnzSghP1Pzh4tIYqlnWZLKdx6qIRgePNBYC/etruPRtXu5/PgqZpbl4POm5t98Ac9dAGSGfu5yJWNDeWaQTF+EYCA2rxe/J8KUnB6umrVfo8kJcOzEAqp0QZ+IuEBhWZLGyo0NrNzYGNW27T0h7n5lB2W56Vw4r5x5lflxri65WA4S8j6MN3ISHsrdLmdMMAZOKOng+T15hB3OiuH3RFgysYUzKjSanAhTi7OYV6mFR0TEHak5JCdJJxiO8N2Hakb8vIa2Hu55dSe/fWEre1u741BZcury/ScA6eGvulzJ2HKsw1aM/m0XZ6rtIiHKctNZPE3zrIuIexSWJSnc9dJ2tjR2jPr5Wxo7uPXpzTzw1m46A6EYVpZ8wmzBmt34w5/FoI+lR6IoPcRorw31eyJU53bzL8fVMylH0xcmQkGWn7NmleLVBb0i4iK1YYjrmjsC/OypzY73E7Hw6rZm3q5r4fy5ZSyuLh53v2Qtlm7/DQD4I0tdrmZsCIQNa5sz+euWYiL0vh482He/j4baLhIvO93LubPLSPNpTEdE3KWwLK770eMbae2O3WhwVzDMQ2/v4bVtzVy6sIKZZbkx27fbQuYpADJCP8Dog6FBDRSQAS6YeJCpuT3ctbGUSGT41Os1lgxvhKvnNGo0OYHSfR7OnVOmpaxFJCkoLIuraupbufe1nXHZd0NbD3e+uJ05E3K59JgKinPS43KcRLH0EPD9AmOL8No5bpeTdIYKyIvL28nxR4DeGVWyfBFahpkVw++JMCW3h6tmaraLRPJ5DOfMLiUvQ1PEiUhyUFgWV337oXVE4pxDNuxtY3NDO6dNL+a82WWk+8fmaFWP9zYAMkL/7XIlySPagNzfoQVKnqvPG7QVQ20X7vAYOHNWyZj/w1ZExheFZXHNinf28MrW5oQcKxyxPL95P2/uPMiF88o5YUohnjGUgiLsJ+x5Fm/kfDwUuV2Oq0YTkI90bEknL+7NPaoVQ20X7vEYOG16CRX5umhVRJKLwrK4ojsY5vsr1if8uO09Ie5/czevbmvmsoUVTCkeGysBdvl6L+pLD1/rciXuiEVA7m+gBUrUduEeY+DU6cVMLs5yuxQRkaMoLIsr7nh+K3UHulw7/u6DXfz6ua0snJjPxfMnUJCV5lotwwmbGjAtpIW+hCF564y1WAfk/g4tUHKoFUNtF+4xBk6dVjxm/nAVkdSjsCwJt6+1m18+s8XtMgB4u66F9XtaOXNmKWfPKsWfZEtnWyJ0+74BgM+e424xCRDPgHykY0s6ebY+j2xfWG0XLjkUlKeWKCiLSPJSWJaE+8XKWjoDYbfLeFcwbHl6QwNv7DjARQsmcOzEArdLelfQ8yAAGcGfYEYwL/BYksiA3F95ZpCPz9zPtPxutV24wBg4bbpGlEUk+SksS8L928VzyE738dsXthEIxScIjcbBriB/fn0Xr2xt4rKFlVQVuHuhkaWLoPcujJ2Cl2pXa4k1twJyf8bA/GL3WoFS2aGL+dSjLCJjgbE2eUdUFi1aZFetWuV2GRInO5s6+e7DNTxes8/tUo5igBOnFLJkXjm5Ls332u39DmHPG2QF78KQ50oNsZQMAVnc5/MYzphZQqXLf4yKiPRnjFltrV000GMaWRbXTC7O4vZPLeLF2v18e3kNG/e1uV3SuyywascB3tndwrmzyzhtRjE+T+L6mSPsIex5A1942ZgOygrI0l+6z8PZs0sp0TzKIjKGaGRZkkI4YvnTqzv4yRObONAZdLucoxRnp7H0mArmViQmuHb4PggmQlbw/zBj7G9aBWQZSHa6l3Nml5GfqZX5RCT5aGRZkp7XY/jkqVN5/7FV/M+Tm7jnlR2E4r203wg0dQS4+5UdzCzLYekxFZTnZcTtWCGzGkyE9NDXxkxQVkCWoRRk+TlndilZaWPj9Swi0p9GliUpbd7XxrcfquH5zfvdLuUoHgOLq4u5YG45mWmxXTrbEqbTfyUA2cG/xXTfsaaALNEozU3n7FmlpPmSa1pGEZH+NLIsY87M8lzu/uxinqzZx3cfrmF7U6fbJb0rYuHlrU2sqTvIBXPLObm6KGZLZwc89wKQGbw1JvuLNQVkGYnJRVmcOr0Yr2d8TnsoIqlBYVmS2gXzyjlrVil3vriNW5+upa0n5HZJ7+oMhHlwTT2vbWvm0oUVTC/NcbQ/Sxsh7314IsfgoSpGVTqngCyjsXBiPguq8t0uQ0TEMbVhyJjR2NbDfz+2gftW15FE7czvmleRx9JjKijKHt2S1F3efyPi2UxW8B4M7i7UoIAso+XzGk6dVsykIs2hLCJjh9owZFwozU3n5iuP5VOnTuXG5et4ffsBt0s6TM2eVjbta+P0GSWcM7uUdF/0/cwRdhLxbMYf/rhrQVkBWZzKyfBx1swSCrJG9wejiEgycjSybIwpAv4MTAW2Ax+21g6YYIwxXmAVsNtae1k0+9fIsgzlwTX13LRiPfUt3W6XcpTcDB8XzZvA8ZMLMMP0M1ssnf4rAMgK3ochthcNDkUBWWKlPC+d02eUkOFP3OtXRCRW4jmy/DXgKWvtTcaYr/Xd/vdBtv0ysB7G8AoLklTef2wlF84r57ZntvDr57bQHUyeYNfWHeK+N+p4ZVvv0tmTh/hIOmReACAjdGNCgrICssTa7Ak5HD+pEI8u5BORccjpyPJG4Bxr7R5jTAXwjLV29gDbTQTuAr4H/LNGliXW6g928YNHNrB8Tb3bpRzFAMdOKuCi+ROOWpDBEqTT/2GwWWSH/hi3GhSQJR7SfB4WVxepP1lExryhRpadhuWD1tqCfrcPWGsLB9juPuAHQC5ww1Bh2RhzLXAtwOTJk0/csWPHqOuT1PP69mZuXL6Otbtb3S7lKGleD2fNKuHMmaX4vb1zzvZ4fk3I+yiZwdvxUBrT4ykgSzyV5KRx+owSstN16YuIjH2O2jCMMU8CEwZ46BtRHvwyoMFau9oYc85w21trbwduh96R5WiOIXLISVOLePAfz+C+1XXc/NhG9rf3uF3SuwLhCE+ub2D1jgNcvKCCBVWWkPdRvJFTYxaUFZAlEeZV5rGwKl9tFyKSEoYNy9baCwZ7zBizzxhT0a8No2GAzU4H3m+MWQpkAHnGmHustZ8YddUiQ/B4DB8+aRKXHDOBW5+u5c4XtxMIJ09IPNAZ5H9f28mXLruFXD+kh7/kaH8KyJIoGX4Pp04vpiI/0+1SREQSxunnZw8CVwM39f37wJEbWGv/A/gPgL6R5RsUlCURcjP8/MfSuXzs5Ml89+H1PLl+n9slvauicDe5WU08svoyvOH9LJk3gZwRfJytgCyJVlmQweLq4pgv8S4ikuychuWbgL8YYz4L7AQ+BGCMqQTusNYudbh/EcemlmRzx9WLeH5zI99eXsPmhnaXK7J8ZsmvAXhjy0nAAd7Z3cJ5s8s4dXrJoEsDKyCLG/xew4lTCpnmcIVKEZGxSiv4SUoJhSPc88oO/ufJzbR0BV2p4fhpr7N00XLufPJa6psnHvZYSU46lx4zgdkTemdYVEAWN2k0WURShVbwE+nj83q45vRqlh1XxU+e2MSfXttJOIFrZ/u8AZYuWs6B9oKjgjLA/vYefv9KHRPKprDPMwergCwu0GiyiMh7FJYlJRVmp/GdDyzgE6dM4dsPrePF2qaEHPeyk/4OwF1Pf/6w+63xQsYEyF8IxsPevvsn+3bxobkRSrI164AkhkaTRUQOp7AsKW32hFz++LlTeGzdXr6/Yj07mjrjdqy8rIPMn7yW1bWL6OjOPSogv6ttI3TuxEQC7AJ+vc/LBfPKOWlqEZ5hls4WGa3sdC8nTilkYqEWGBER6U9hWQS4aP4Ezpldym9f2MYvnq6lIxCO+TE+f+EvAHh84+ewE04cNCAfqSMQ5oG36nltWzOXLqxgWok+GpfY8XpgXkU+8yrzBr24VEQklekCP5EjNLR288NHN3L/m3XE6sdjSnYdH/Q/xGOV17H2wKW9dw4RkIcyvzKPpQsqKMxOi01xkrIqCzI4cUohuRn+4TcWERnH4rbcdbwpLIub1uw6yI3L1/HGzoOO9zWtYz9n37odgCeunEPddGejwz6P4YyZJZwzq4w0n2f4J4j0k5Ph48QphVQVaHERERFQWBZx5O9v7uamRzawt7Xb0X5yD3Rz5e1r3739188voLUow9E+8zJ8XDR/AsdNKsCon1mGkeH3ML8yn5llOVqqWkSkH4VlEYc6AyF+9cwWbn9+K91BZ1O3TdxykCX31QKwe0ouKy+fQTDd2cwDkwozuWxhJZOKdHGWHM3nNcydkMecilz8Xn0SISJyJIVlkRipO9DJD1Zs4OF39jjaj4lYjn9hN8e+3DtJ3Kqzqli7eALWwWifAY6bVMBFCyaQpx5UATwGZpbnML8ynwy/poITERmMwrJIjL26tYkbl9dQs6fV0X7Su0Jc+JdNlOztnbLukY/NYu/kPEf7TPN5OGdWKWfMKMGnUcSU5DEwpTibYybmk5OuSY9ERIajsCwSB5GI5c+rdvGjxzbS1DGyGS2OVNjYyQd+VwNAT7qXBz4zj468dGf7zPJzyYIKFlTlO9qPjB1eD0wrzWFuRZ5CsojICCgsi8RRa3eQnz25mbte3k4w7ODnyVqqa5o556FtAGydU8gLS6cSdvjx+bTSbC47ppIJ+c4uJpTk5fcaZpbnMmdCrtotRERGQWFZJAG2Nrbz3YfX8/SGBkf78YQiLH5yJ3PW7AfgpQsns/G4UnAw24XHwKKpRVw4t5wsjTiOGxl+D7Mn5DKzLFdTCIqIOKCwLJJAz2xs4DsP1bClscPRfjLbgyy9Zz15Lb0tHss/NZf9FdnO9un3ct6cMk6ZVqzV2sawouw0ZpbnMLU4W/8fRURiQGFZJMFC4Qh3vbyDW57cRGt3yNG+Surbed/dGwBoLUhnxcfn0JXjbLaL0tx0Lj2mglnluY72I4nj8xgmFWUxqzyH4hxn/ewiInI4hWURlzR3BPjx4xu59/VdhCPO+plnv9XIaY/vBGD9caW8dv4kIg4/ep9dnsulCysoUfhKWjkZPmaW5TCtNJt0n/qRRUTiQWFZxGXr97Ry4/J1vLK12dF+fIEwZ6zYTvXGAwA8875qts0tctTP7PUYTp1WzHlzynRxWJLweQ2TCrOoLsnWhZkiIgmgsCySJB55Zw/ff2Q9u5q7HO0nu6WHD/xuHWmB3tUE//aZ+RwszXS2z3QfF84r58QphXi0dHbCeQxUFGQytTiLqoJMzZEtIpJACssiSaQ7GOaO57fyy2e20BkIO9pXxY5WLr53EwCNFdk8/qGZBDKdzXZRmZ/BpQsrqS5xdjGhRKckJ43qkmwmFWVpZF9ExCUKyyJJaF9rNz98ZAN/e2s3Tn4MTcSy4NW9LHpuNwBvnVrBW2dUOlo6G+CYqnwuWTCBgqw0R/uRw/k8hrK8dCYW9o4gZ6YpIIuIuE1hWSSJvbnzADcur+GtXQcd7cffE+K8+7dQubMNgCeunEHd9AJn+/QazpxZylkzSzWPrwMZfg+VBZlUFWRSkZ+hFgsRkSSjsCyS5Ky13P/Gbm5+bAP7Wnsc7SuvuZsP/mZt734N/PXzC2grdHaRWH6mn4vnT+DYSQWO9pMqfB5DaW46ZXnplOdlUJydhlEfuIhI0lJYFhkjOgMhfrGyljue30ZPKOJoX5M2H+SC+2sBqKvOY+Wy6YTSnX3kP6Uoi8sWVlJV6OxiwvHG5zGU5KZRlpvxbjj2aLEQEZExQ2FZZIzZ1dzJ9x5ez6Pr9jrajwlHOPG53Rzz2j4AXjt3IutOKnc01ZwBTphcyIXzy8nNcLY4ylhkTO9Ie2FWGsU5aRRlp1GUpXAsIjKWKSyLjFEvb2nixuXr2LC3zdF+0juDXPTnTRQ39E5Zt+Kq2eyb5Gz1vqw0L8uOq+TUacV0ByOEnCy6kqR8HkNepo/8zN5gXJiVRmGWXz3HIiLjjMKyyBgWjlj+97Wd/OSJTTR3BBztq7Chkw/cWQNAV6aPB6+ZR2ees9kuphZn8fWlczl9Rglt3SHauoO0dgdp7Q7R1h2ioyfkaLaPePMYyEr3kZfhIzfDT16Gj7xMP7kZPrLSnE3DJyIiY4PCssg40NIV5JYnN3P3K9sJhp0tnT19XRNnPbwdgNp5Rbx08VTCfmejpWfOLOG/LpvHrPLDR6yttXQHI3QFw71fgTDdR3wfiliC4QihcO+/TgapvR7weTz4vAa/14PPY8jwe8lM85Lp95Lh95LV931mmpd0n0cX34mIpDiFZZFxpLahne88VMOzmxod7ccbinDKEzuZ9fZ+AF64eAqbF5Y46mf2eQwfXzyZry6Z5Wh+5kjEEoz0hudQ2GKxWAv9360MvaUaY/B5TG849njUOywiIiOmsCwyDj29YR/ffWg9W/d3ONpPZluAy+5eT05bEIAHr55L0wRnq/cVZPn56gWz+MQpU/AqvIqISJJTWBYZp4LhCL9/cTs/e3ozbd0hR/sq3d3OZfdsAOBgUQaPXDWb7mxns13Mr8zjgX88XRfEiYhIUhsqLOs3mMgY5vd6+PxZ03jmhnP42MmTcDKI21iVw53/diIvXzCJguZuPnbrGk5+ciee8Ojne15cXaygLCIiY5pGlkXGkXX1Ldy4vIbXtjU72o8vEOash7YxZfNBAJ5eNo0dswtH1M9clJ3GyhvOIT8z9eZiFhGRsUUjyyIpYn5lPn/5wqncetXxVBWMfpW9UJqXp6+Ywf9ddwwhn+G8B7by6ZtXk7+/K+p9/POSWe8G5TfueIP7PnIfEYerEoqIiCSaRpZFxqnuYJjbn9vKbc9soSsYdrSvym0tXPSXzQDsq8rmyStnEsgYfA7iuRV5PPxPZ+DxGDoaOrhl2i3YsGXulXO5/A+Xa6o2ERFJKhpZFklBGX4vXzp/Jk/fcDbLjqt0tK/66nx+/68n8sYZlZTv7uDjt7zFcS/sxgwyIfL/u2zeu1O4PX7D40SCEULdITbcv4GV/7XSUS0iIiKJpLAsMs5V5Gdyy0eP569fPJWFE/NHvR/rMaw5vZI/fvk49k7M4fgX93DNf6+mclvLYdtdsmACp04vBmDvW3upua+GcKB3ZDvYGeSV/3mFN37zxuj/g0RERBJIbRgiKcRay32r67j5sY00tvU42ld+UxdX3LEOgLDXcP/nFhAsyeTJfz6bSUVZWGv5zaLfsOeNPUc915fp40P3fYhZS2c5qkFERCQW1IYhIkDvancfWjSJlTecw3VnTyfNN/q3gJbiTO7890U8/YHpeMOWD/36HT7xZB3lab29zDX/V8P+jfsHfG6oK8R9H7qP+lX1oz6+iIhIImhkWSSF7Wjq4LsPr+eJmn2O9uMJRzjzpT1Me6l3FHnJj5bw4g9fpLOxc8jnpRek84XVX6BwWqGj44uIiDihkWURGdCU4mx+86lF/PFzi5ldnjvq/US8Hpb9fCk3NNxA6bxSnrjhCXpah2/zCLQGuPOsO+ncP3SoFhERcYvCsohw+owSVnz5TL69bD4FWSNfROT4yQVcfnwV2aXZfOKxT+DL8hHuGX66OhuxdDR0cNd5dxHsDI6mdBERkbhSWBYRALwew6dOncozN5zDNadNxRfl2tnGwDffN//duZMf+dIjRALRLz4SCUZo2tzEvR+4l4iDpbVFRETiQWFZRA5TkJXGt94/n0e+fCZnziwZdvvLj6/iuEkFAOx6eRe1j9aOeKW+cHeYXS/uYvnnl5PM11GIiEjqUVgWkQHNLM/l7s8u5jefWsTU4qwBt8lO8/K1i+cAvS0VD372QUJdoVEdL9gZZN2f1/Hcd54bdc0iIiKxprAsIkNaMq+cx796Nv9xyRxy0w9f4vofzp1BWV4GAGvuWkPLzpaBdhG1YGeQF256gTd//6aj/YiIiMSKwrKIDCvN5+ELZ0/n6RvO4cOLJuIxMKkok8+dWQ1AT1sPj/3zYwQ7nF+kF+oKseIfVlD7WK3jfYmIiDjlG34TEZFepbnp3HzlsXzylKm094RI93kBeOOONwi0B2J2nFBXiL988C98+vlPU3F8Rcz2KyIiMlIKyyIyYsdMzD/s9oKPLiASirDjmR3Ur6qn62AX/gw/wc7giC/2OyTYEeQP5/+BL7zxBQqmFsSgahERkZHTCn4iEnNdzV3Ur6pn92u72f7Mdva8sYdgZxBfuo9AewAbie59x3gMuRNzue7N68gsyoxz1SIikqqGWsFPYVlEEqJ9bzv1q+qpe7WOHc/sYO+avURCEbx+Lz1tPTDIW5HH76Fsfhmfeekz+DNHvmCKiIjIcBSWRSTpWGtprWul/vV66l6uY/tz22lc2wiA8RoCbe/1QPsyfUw9eyofe+hjeLy6LllERGJLYVlExgRrLQe2HqB+VT27XtzFjud3sH/9fozHEOoKcdq/ncaSHy5xu0wRERlnhgrLusBPRJKGMYai6UUUTS9iwUcWABAJR2ja1ET9qnqKZha5XKGIiKQahWURSWoer4fSuaWUzi11uxQREUlBav4TERERERmEwrKIiIiIyCAchWVjTJEx5gljzOa+fwsH2a7AGHOfMWaDMWa9MeZUJ8cVEREREUkEpyPLXwOestbOBJ7quz2QW4BHrbVzgGOB9Q6PKyIiIiISd07D8jLgrr7v7wI+cOQGxpg84CzgtwDW2oC19qDD44qIiIiIxJ3TsFxurd0D0Pdv2QDbTAMagTuNMW8aY+4wxmQ7PK6IiIiISNwNG5aNMU8aY9YO8LUsymP4gBOA26y1xwMdDN6ugTHmWmPMKmPMqsbGxigPISIiIiISe8POs2ytvWCwx4wx+4wxFdbaPcaYCqBhgM3qgDpr7at9t+9jiLBsrb0duB16V/Abrj4RERERkXhx2obxIHB13/dXAw8cuYG1di+wyxgzu++u84Eah8cVEREREYk7p2H5JmCJMWYzsKTvNsaYSmPMin7b/RPwR2PM28BxwPcdHldEREREJO4cLXdtrW2id6T4yPvrgaX9br8FLHJyLBERERGRRNMKfiIiIiIig1BYFhEREREZhMKyiIiIiMggFJZFRERERAahsCwiIiIiMgiFZRERERGRQSgsi4iIiIgMwlibvCtKG2MagR0OdlEC7I9ROalM5zE2dB5jQ+cxdnQuY0PnMTZ0HmNH53LkplhrSwd6IKnDslPGmFXWWi2G4pDOY2zoPMaGzmPs6FzGhs5jbOg8xo7OZWypDUNEREREZBAKyyIiIiIigxjvYfl2twsYJ3QeY0PnMTZ0HmNH5zI2dB5jQ+cxdnQuY2hc9yyLiIiIiDgx3keWRURERERGbUyHZWNMkTHmCWPM5r5/CwfYZrYx5q1+X63GmK/0PfYtY8zufo8tTfh/RJKI5lz2bbfdGPNO3/laNdLnj3dRviYnGWNWGmPWG2PWGWO+3O+xlH5NGmMuNsZsNMbUGmO+NsDjxhjzs77H3zbGnBDtc1NJFOfx433n721jzEvGmGP7PTbgz3iqiuJcnmOMaen3M/v/on1uKoniPP5rv3O41hgTNsYU9T2m12QfY8zvjDENxpi1gzyu98h4sNaO2S/gZuBrfd9/DfjhMNt7gb30zqUH8C3gBrf/O5LhK9pzCWwHSpz+vxivX9GcB6ACOKHv+1xgEzCv73bKvib7fj63ANOANGDNofPSb5ulwCOAAU4BXo32uanyFeV5PA0o7Pv+kkPnse/2gD/jqfgV5bk8B3hoNM9Nla+RngvgfcDT/W7rNfneuTgLOAFYO8jjeo+Mw9eYHlkGlgF39X1/F/CBYbY/H9hirXWy0Ml4NdJzGevnjxfDngdr7R5r7Rt937cB64GqRBWYxE4Gaq21W621AeBees9nf8uAP9herwAFxpiKKJ+bKoY9F9bal6y1B/puvgJMTHCNY4WT15Vek+8Z6bn4GPC/CalsjLHWPgc0D7GJ3iPjYKyH5XJr7R7oDSBA2TDbf5SjfwCv7/uo4nep2jrQJ9pzaYHHjTGrjTHXjuL5492IzoMxZipwPPBqv7tT9TVZBezqd7uOo/+IGGybaJ6bKkZ6Lj5L70jUIYP9jKeiaM/lqcaYNcaYR4wx80f43FQQ9bkwxmQBFwN/7Xe3XpPR03tkHPjcLmA4xpgngQkDPPSNEe4nDXg/8B/97r4N+A69P4jfAX4MfGZ0lSa/GJ3L06219caYMuAJY8yGvr90U0YMX5M59P5C+Iq1trXv7pR6TR7BDHDfkdP1DLZNNM9NFVGfC2PMufSG5TP63Z3yP+P9RHMu36C3ta+97xqDvwMzo3xuqhjJuXgf8KK1tv/oqV6T0dN7ZBwkfVi21l4w2GPGmH3GmApr7Z6+jxkahtjVJcAb1tp9/fb97vfGmN8AD8Wi5mQVi3Npra3v+7fBGPM3ej/aeQ4Yyf+LMS0W59EY46c3KP/RWnt/v32n1GvyCHXApH63JwL1UW6TFsVzU0U05xFjzELgDuASa23TofuH+BlPRcOey35/6GKtXWGM+aUxpiSa56aQkZyLoz4B1mtyRPQeGQdjvQ3jQeDqvu+vBh4YYtujeqD6wswhlwMDXl2aIoY9l8aYbGNM7qHvgQt575yN5P/FeBbNeTTAb4H11tqfHPFYKr8mXwdmGmOq+z4J+ii957O/B4FP9V3xfQrQ0tfuEs1zU8Ww58IYMxm4H/iktXZTv/uH+hlPRdGcywl9P9MYY06m9/dqUzTPTSFRnQtjTD5wNv3eN/WaHDG9R8aD21cYOvkCioGngM19/xb13V8JrOi3XRa9b175Rzz/buAd4G16XzQVbv83JfO5pPcq2jV9X+uAbwz3/FT7ivI8nkHvx19vA2/1fS3teyylX5P0Xsm9id6rtr/Rd991wHV93xvgF32PvwMsGuq5qfoVxXm8AzjQ7/W3qu/+QX/GU/UrinN5fd+5WkPvxZKnDfXcVP0a7jz23b4GuPeI5+k1efj5+F9gDxCkdxT5s3qPjP+XVvATERERERnEWG/DEBERERGJG4VlEREREZFBKCyLiIiIiAxCYVlEREREZBAKyyIiIiIig1BYFhEREREZhMKyiIiIiMggFJZFRERERAbx/wHLJ+ADgSI3XgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "vx, vy = 0.3, 0.1 # change these !\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "\n",
    "circle = Circle((0, 0), R, alpha=0.4)\n",
    "ax.add_artist(circle)\n",
    "\n",
    "wheels = []\n",
    "for i in range(3):\n",
    "    px, py, theta = wheel_pose(i)\n",
    "    wheel = Rectangle((-w/2, -r), w, 2*r)\n",
    "    wheel.set_transform(Affine2D().rotate(theta).translate(px, py))\n",
    "    wheels.append(wheel)\n",
    "ax.add_collection(PatchCollection(wheels))\n",
    "\n",
    "\n",
    "plt.arrow(0, 0, vx, vy, head_width=0.05, color=\"dodgerblue\")\n",
    "for i in range(3):\n",
    "    px, py, theta = wheel_pose(i)\n",
    "    para, perp = wheel_velocity(vx, vy, i)\n",
    "    plt.arrow(px, py, vx, vy, head_width=0.05, color=\"dodgerblue\")\n",
    "    plt.arrow(px, py, np.cos(theta)*perp, np.sin(theta)*perp, head_width=0.05, color=\"yellowgreen\")\n",
    "    plt.arrow(px, py, -np.sin(theta)*para, np.cos(theta)*para, head_width=0.05, color=\"purple\")\n",
    "\n",
    "ax.axis('equal'); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Jacobian Matrix\n",
    "\n",
    "As mentioned above, the motion *perpendicular* to the wheel is not resisted by the rollers: they will roll freely\n",
    "at exactly the right speed to accommodate any movement in this direction. Because this allows the wheel to \"slide\", we also call this the **sliding direction**.\n",
    "\n",
    "On the other hand, the contact of each wheel must move exactly opposite to the parallel or **driving direction** of the velocity, because it cannot slide in that direction. As a consequence, if the radius of each wheel is $r$, the wheel's angular velocity $\\omega^i$ (i.e., the angular rate at which the wheel spins about its axle)\n",
    "must satisfy\n",
    "\n",
    "$$\\omega^i = \\frac{v^i_\\parallel}{r} = \n",
    "\\frac{1}{r} (-\\sin \\theta^i, \\cos \\theta^i) \\begin{pmatrix}v_x \\\\ v_y \\end{pmatrix}\n",
    "$$\n",
    "\n",
    "Because we have three wheels, we can stack three instances of the above equation\n",
    "to obtain a matrix equation\n",
    "that maps a desired robot velocity $v=(v_x,v_y)$ to the appropriate commanded wheel velocities:\n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\\omega^1 \\\\ \\omega^2 \\\\ \\omega^3 \\end{pmatrix} \n",
    "= \\frac{1}{r}\n",
    "\\begin{pmatrix}- \\sin\\theta^1 & \\cos\\theta^1 \\\\ - \\sin\\theta^2 & \\cos\\theta^2 \\\\ - \\sin\\theta^3 & \\cos\\theta^3 \n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}v_x \\\\ v_y \\end{pmatrix}\n",
    "$$\n",
    "\n",
    "in which the $3\\times2$ matrix is called the **Jacobian matrix**.\n",
    "\n",
    "For the three regularly arranged omni wheels with $\\theta^1=0$, $\\theta^2=2\\pi/3$, and $\\theta^3=4\\pi/3$ this becomes\n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\\omega^1 \\\\ \\omega^2 \\\\ \\omega^3 \\end{pmatrix} \n",
    "= \\frac{1}{r}\n",
    "\\begin{pmatrix}0 & 1 \\\\ -0.866 & -0.5 \\\\ 0.866 & -0.5 \n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}v_x \\\\ v_y \\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def jacobian_matrix():\n",
    "    \"\"\"Compute matrix converting from velocity to wheel velocities.\"\"\"\n",
    "    rows = []\n",
    "    for i in range(num_wheels):\n",
    "        _, _, theta = wheel_pose(i)\n",
    "        rows.append([-np.sin(theta), np.cos(theta)])\n",
    "    return np.array(rows)/r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can verify that doing the matrix-velocity multiplication yields the same result as using the `wheel_velocity` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "wheel 0: 0.5\n",
      "wheel 1: -1.549038105676658\n",
      "wheel 2: 1.0490381056766573\n",
      "\n",
      "Using matrix method:\n",
      "[ 0.5        -1.54903811  1.04903811]\n"
     ]
    }
   ],
   "source": [
    "for i in range(3):\n",
    "    para, perp = wheel_velocity(vx, vy, i)\n",
    "    print(f\"wheel {i}: {para/r}\")\n",
    "\n",
    "print(\"\\nUsing matrix method:\")\n",
    "print(jacobian_matrix() @ gtsam.Point2(vx, vy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Omni-directional Movement in Practice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Additional transforms might be required.\n",
    "\n",
    "To summarize, given a velocity, $v$, in the robot's body coordinates,\n",
    "we compute the corresponding vector angular wheel velocities using the Jacobian mapping, $\\omega = Jv$.\n",
    "The above exposition can be directly generalized to an arbitrary number of wheels, mounted in arbitrary locations. Some care needs to be taken that the chosen configuration is not singular, which would result in the robot being able to move only in one particular direction (or not at all!).\n",
    "\n",
    "In general, we would like to be able to specify the robot's direction of motion relative\n",
    "to arbitrary coordinate frames, e.g., to a fixed warehouse coordinate frame,\n",
    "or a *world coordinate frame*.\n",
    "In this chapter we have assumed that the robot orientation is aligned with the warehouse orientation, so conveniently velocity vectors in either frame have the same coordinates.\n",
    "However, if we were to include the robot's orientation in the state representation, we would need to rotate the desired world velocity into the body coordinate frame. We will see how to do this later in this book.\n",
    "\n",
    "If we *do* care about orientation, omni wheels can actually be used to *rotate* the robot as  well as to induce pure translational motion. Think about applying the same constant angular velocity to all three wheels: the robot will now simply rotate in place. To represent this mathematically, we could add a *third* column to the Jacobian matrix that maps a desired angular velocity into wheel velocities. This would allow us to command any combination of linear and angular velocity to the robot.\n",
    "\n",
    "Finally, there is a more general type of wheel, the \"mecanum wheel\", for which the rollers are *not* aligned with the wheel rim, but are rotated by some (constant) amount. That complicates the math a little bit, but in the end we will again obtain a Jacobian matrix that maps robot velocity to wheel angular velocities. We will not discuss this further."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Gaussian Dynamic Bayes Net\n",
    "\n",
    "> Modeling stochastic action in continuous spaces can again be done using DBNs.\n",
    "\n",
    "Modeling with continuous densities can be done with *continuous* dynamic Bayes nets (DBNs). \n",
    "Omni wheels are notoriously flaky in the execution, especially on uneven ground or carpet.\n",
    "Hence, it is important to model the uncertainty associated with commanded movements.\n",
    "Because we are modeling our robot as a discrete-time system,\n",
    "we let $u_k = v \\Delta T$ denote the displacement that occurs when the robot\n",
    "moves with velocity $v$ for time $\\Delta T$.\n",
    "Since the velocity can be in any direction, $u_k\\in\\mathbb{R}^2$.\n",
    "In the absence of uncertainty, our motion model would simply\n",
    "be $x_{k+1} = x_k + u_k$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Given the uncertainty inherent to the use of omni-wheels, we will instead develop a probabilistic\n",
    "motion model.\n",
    "Recall that we created a DBN with discrete actions, states, and measurements in the previous chapter. Here we do the same, except with *continuous* densities, specifying a DBN that encodes the joint density $p(X|U)$ on the continuous states $X\\doteq\\{x_k\\}_{k=1}^N$,\n",
    "\n",
    "$$\n",
    "p(X|U) = p(x_1) \\prod_{k=2}^{N} p(x_{k}|x_{k-1}, u_{k-1})\n",
    "$$\n",
    "\n",
    "where we assumed the **control variables** $U\\doteq\\{u_k\\}_{k=1}^{N-1}$ as *given* to us."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us assume we *know* that the robot starts out in the warehouse at coordinates $\\mu_1 =(20,10)$, with some noise/error in position of magnitude on the order of plus or minus half a meter. We can encode this knowledge into a continuous Gaussian prior as\n",
    "\n",
    "$$p(x_1) = \\mathcal{N}(x_1;\\mu=\\mu_1,\\Sigma=P) \\propto \\exp\\{ - \\frac{1}{2} (x_1-\\mu_1)^TP^{-1}(x_1-\\mu_1) \\}$$\n",
    "\n",
    "with the covariance matrix $P=\\text{diag}(0.5^2,0.5^2)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A Gaussian *motion model* is a bit trickier. \n",
    "Here we need a **conditional Gaussian density** that depends on the current state\n",
    "and the commanded motion.\n",
    "The usual way to incorporate uncertainty into a motion model is to assume that\n",
    "a random disturbance is added to the nominal motion of the robot.\n",
    "We can model this as\n",
    "\n",
    "$$x_{k+1} = x_k + u_k + w_k$$\n",
    "\n",
    "in which $w_k$ is a random vector whose probability distribution is typically based on empirical\n",
    "observations of the system. In our case, we will assume that this disturbance\n",
    "is Gaussian, $w_k \\sim N(\\mu,\\Sigma)$.\n",
    "In this case, it can be shown that the conditional distribution of $x_{k+1}$\n",
    "given $x_k$ and $u_k$ is given by a **conditional Gaussian density** of the form\n",
    "\n",
    "$$\n",
    "p(x_{k+1}|x_{k}, u_k) = \\mathcal{N}(x_{k_1}; x_{k} +  u_k, \\Sigma)\n",
    "$$\n",
    "\n",
    "i.e., the conditional distribution for the next state is a Gaussian, with mean equal to the\n",
    "nominal (i.e., noise-free) next state $x_k + u_k$, and variance equal to the noise variance.\n",
    "\n",
    "It is easy to generalize this motion model to linear systems of the form\n",
    "\n",
    "$$\n",
    "x_{k+1} = A x_k + B u_k + w_k\n",
    "$$\n",
    "\n",
    "which is the general form for a linear time-invariant (LTI)) system with additive noise.\n",
    "The matrix $A \\in \\mathbb{R}^{n \\times n}$ is called the system matrix,\n",
    "and the matrix $B$ maps control input $u_k$ to a vector in the state space.\n",
    "The motion model for this general LTI system is given by\n",
    "\n",
    "$$\n",
    "p(x_{k+1}|x_{k}, u_k) = \\mathcal{N}(x_{k+1};\\mu=A x_{k} + B u_k, \\Sigma=Q)\n",
    "$$\n",
    "\n",
    "where $Q$ is the traditional symbol used for motion model covariance. If we assume the robot executes actions $u_k$ with a standard deviation of 0.2m, then we will have $Q=\\text{diag}(0.2^2,0.2^2)$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gaussian DBNs in Code\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In GTSAM, a Gaussian density is specified using the class `GaussianDensity`, which corresponds to a negative log-probability given by \n",
    "\n",
    "$$\\frac{1}{2 \\sigma^2} \\|x - \\mu\\|^2.$$\n",
    "\n",
    "with a named constructor `FromMeanAndStddev` that takes $\\mu$ and $\\sigma$. For their part, Gaussian *conditional* densities are specified using the class `GaussianConditional`, which correspond to a negative log-probability given by \n",
    "\n",
    "$$\\frac{1}{2 \\sigma^2} \\|x - (A_1 p_1 + A_2 p_2 + b))\\|^2.$$\n",
    "\n",
    "i.e., the mean on $x$ is a linear function of the conditioning variables $p_1$ and $p_2$. Again, a named constructor `FromMeanAndStddev` is there to assist us. Details in the GTSAM 101 section at the end of this section, as always."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of the above is used in the following piece of code, which builds the dynamic Bayes net:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "gaussianBayesNet = gtsam.GaussianBayesNet()\n",
    "A, B = np.eye(2), np.eye(2)\n",
    "motion_model_sigma = 0.2\n",
    "for k in reversed(indices[:-1]):\n",
    "    gaussianBayesNet.push_back(gtsam.GaussianConditional.FromMeanAndStddev(\n",
    "        x[k+1], A, x[k], B, u[k], [0, 0], motion_model_sigma))\n",
    "p_x1 = gtsam.GaussianDensity.FromMeanAndStddev(x[1], [20,10], 0.5)\n",
    "gaussianBayesNet.push_back(p_x1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n<!-- Generated by graphviz version 2.50.0 (0)\n -->\n<!-- Pages: 1 -->\n<svg width=\"350pt\" height=\"116pt\"\n viewBox=\"0.00 0.00 350.00 116.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 112)\">\n<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-112 346,-112 346,4 -4,4\"/>\n<!-- varu1 -->\n<g id=\"node1\" class=\"node\">\n<title>varu1</title>\n<polygon fill=\"none\" stroke=\"black\" points=\"54,-108 0,-108 0,-72 54,-72 54,-108\"/>\n<text text-anchor=\"middle\" x=\"27\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">u1</text>\n</g>\n<!-- varx2 -->\n<g id=\"node6\" class=\"node\">\n<title>varx2</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"99\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x2</text>\n</g>\n<!-- varu1&#45;&gt;varx2 -->\n<g id=\"edge2\" class=\"edge\">\n<title>varu1&#45;&gt;varx2</title>\n<path fill=\"none\" stroke=\"black\" d=\"M45.17,-71.83C54.73,-62.27 66.53,-50.47 76.65,-40.35\"/>\n<polygon fill=\"black\" stroke=\"black\" points=\"79.4,-42.55 83.99,-33.01 74.45,-37.6 79.4,-42.55\"/>\n</g>\n<!-- varu2 -->\n<g id=\"node2\" class=\"node\">\n<title>varu2</title>\n<polygon fill=\"none\" stroke=\"black\" points=\"126,-108 72,-108 72,-72 126,-72 126,-108\"/>\n<text text-anchor=\"middle\" x=\"99\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">u2</text>\n</g>\n<!-- varx3 -->\n<g id=\"node7\" class=\"node\">\n<title>varx3</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"171\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"171\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x3</text>\n</g>\n<!-- varu2&#45;&gt;varx3 -->\n<g id=\"edge4\" class=\"edge\">\n<title>varu2&#45;&gt;varx3</title>\n<path fill=\"none\" stroke=\"black\" d=\"M117.17,-71.83C126.73,-62.27 138.53,-50.47 148.65,-40.35\"/>\n<polygon fill=\"black\" stroke=\"black\" points=\"151.4,-42.55 155.99,-33.01 146.45,-37.6 151.4,-42.55\"/>\n</g>\n<!-- varu3 -->\n<g id=\"node3\" class=\"node\">\n<title>varu3</title>\n<polygon fill=\"none\" stroke=\"black\" points=\"198,-108 144,-108 144,-72 198,-72 198,-108\"/>\n<text text-anchor=\"middle\" x=\"171\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">u3</text>\n</g>\n<!-- varx4 -->\n<g id=\"node8\" class=\"node\">\n<title>varx4</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"243\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"243\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x4</text>\n</g>\n<!-- varu3&#45;&gt;varx4 -->\n<g id=\"edge6\" class=\"edge\">\n<title>varu3&#45;&gt;varx4</title>\n<path fill=\"none\" stroke=\"black\" d=\"M189.17,-71.83C198.73,-62.27 210.53,-50.47 220.65,-40.35\"/>\n<polygon fill=\"black\" stroke=\"black\" points=\"223.4,-42.55 227.99,-33.01 218.45,-37.6 223.4,-42.55\"/>\n</g>\n<!-- varu4 -->\n<g id=\"node4\" class=\"node\">\n<title>varu4</title>\n<polygon fill=\"none\" stroke=\"black\" points=\"270,-108 216,-108 216,-72 270,-72 270,-108\"/>\n<text text-anchor=\"middle\" x=\"243\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">u4</text>\n</g>\n<!-- varx5 -->\n<g id=\"node9\" class=\"node\">\n<title>varx5</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"315\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"315\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x5</text>\n</g>\n<!-- varu4&#45;&gt;varx5 -->\n<g id=\"edge8\" class=\"edge\">\n<title>varu4&#45;&gt;varx5</title>\n<path fill=\"none\" stroke=\"black\" d=\"M261.17,-71.83C270.73,-62.27 282.53,-50.47 292.65,-40.35\"/>\n<polygon fill=\"black\" stroke=\"black\" points=\"295.4,-42.55 299.99,-33.01 290.45,-37.6 295.4,-42.55\"/>\n</g>\n<!-- varx1 -->\n<g id=\"node5\" class=\"node\">\n<title>varx1</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"27\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x1</text>\n</g>\n<!-- varx1&#45;&gt;varx2 -->\n<g id=\"edge1\" class=\"edge\">\n<title>varx1&#45;&gt;varx2</title>\n<path fill=\"none\" stroke=\"black\" d=\"M54.22,-18C56.64,-18 59.11,-18 61.6,-18\"/>\n<polygon fill=\"black\" stroke=\"black\" points=\"61.74,-21.5 71.74,-18 61.74,-14.5 61.74,-21.5\"/>\n</g>\n<!-- varx2&#45;&gt;varx3 -->\n<g id=\"edge3\" class=\"edge\">\n<title>varx2&#45;&gt;varx3</title>\n<path fill=\"none\" stroke=\"black\" d=\"M126.22,-18C128.64,-18 131.11,-18 133.6,-18\"/>\n<polygon fill=\"black\" stroke=\"black\" points=\"133.74,-21.5 143.74,-18 133.74,-14.5 133.74,-21.5\"/>\n</g>\n<!-- varx3&#45;&gt;varx4 -->\n<g id=\"edge5\" class=\"edge\">\n<title>varx3&#45;&gt;varx4</title>\n<path fill=\"none\" stroke=\"black\" d=\"M198.22,-18C200.64,-18 203.11,-18 205.6,-18\"/>\n<polygon fill=\"black\" stroke=\"black\" points=\"205.74,-21.5 215.74,-18 205.74,-14.5 205.74,-21.5\"/>\n</g>\n<!-- varx4&#45;&gt;varx5 -->\n<g id=\"edge7\" class=\"edge\">\n<title>varx4&#45;&gt;varx5</title>\n<path fill=\"none\" stroke=\"black\" d=\"M270.22,-18C272.64,-18 275.11,-18 277.6,-18\"/>\n<polygon fill=\"black\" stroke=\"black\" points=\"277.74,-21.5 287.74,-18 277.74,-14.5 277.74,-21.5\"/>\n</g>\n</g>\n</svg>\n",
      "text/plain": [
       "<gtbook.display.show at 0x10579b190>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "position_hints = {'u': 2, 'x': 1}\n",
    "show(gaussianBayesNet, hints=position_hints, boxes=set(list(u.values())))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Note that we used a constant motion model noise here, which is perhaps unrealistic: if a robot is commanded to go further, it could be expected to do a worse job at going exactly the distance it was supposed to go. How could you modify the above code to account for that?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulating a Continuous Trajectory\n",
    "\n",
    "> Bayes nets are great for simulation.\n",
    "\n",
    "Let us now create a simulated trajectory in our example warehouse, by specifying a \"control tape\" and using the dynamic Bayes net created above to do the heavy lifting for us:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table class='VectorValues'>\n",
       "  <thead>\n",
       "    <tr><th>Variable</th><th>value</th></tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr><th>u1</th><td>2 0</td></tr>\n",
       "    <tr><th>u2</th><td>2 0</td></tr>\n",
       "    <tr><th>u3</th><td>0 2</td></tr>\n",
       "    <tr><th>u4</th><td>0 2</td></tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "VectorValues: 4 elements\n",
       "  u1: 2 0\n",
       "  u2: 2 0\n",
       "  u3: 0 2\n",
       "  u4: 0 2"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "control_tape = gtsam.VectorValues()\n",
    "for k, (ux,uy) in zip(indices[:-1], [(2,0), (2,0), (0,2), (0,2)]):\n",
    "    control_tape.insert(u[k], gtsam.Point2(ux,uy))\n",
    "control_tape\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table class='VectorValues'>\n",
       "  <thead>\n",
       "    <tr><th>Variable</th><th>value</th></tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr><th>u1</th><td>2 0</td></tr>\n",
       "    <tr><th>u2</th><td>2 0</td></tr>\n",
       "    <tr><th>u3</th><td>0 2</td></tr>\n",
       "    <tr><th>u4</th><td>0 2</td></tr>\n",
       "    <tr><th>x1</th><td>20.6469  10.199</td></tr>\n",
       "    <tr><th>x2</th><td>22.8706 9.90054</td></tr>\n",
       "    <tr><th>x3</th><td>24.8089 10.1832</td></tr>\n",
       "    <tr><th>x4</th><td>25.1275 12.2733</td></tr>\n",
       "    <tr><th>x5</th><td>25.0941 14.5563</td></tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "VectorValues: 9 elements\n",
       "  u1: 2 0\n",
       "  u2: 2 0\n",
       "  u3: 0 2\n",
       "  u4: 0 2\n",
       "  x1: 20.6469  10.199\n",
       "  x2: 22.8706 9.90054\n",
       "  x3: 24.8089 10.1832\n",
       "  x4: 25.1275 12.2733\n",
       "  x5: 25.0941 14.5563"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gaussianBayesNet.sample(control_tape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GTSAM 101\n",
    "\n",
    "> The GTSAM concepts used in this section, explained."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A `gtsam.GaussianDensity` class can be constructed via the following named constructor:\n",
    "\n",
    "```python\n",
    "FromMeanAndStddev(key: gtsam.Key, mean: np.array, sigma: float) -> gtsam.GaussianDensity\n",
    "```\n",
    "\n",
    "and a two similar named constructors exists for `gtsam.GaussianConditional`:\n",
    "\n",
    "```python\n",
    "- FromMeanAndStddev(key: gtsam.Key, A: np.array, parent: gtsam.Key, b: numpy.ndarray[numpy.float64[m, 1]], sigma: float) -> gtsam.GaussianConditional\n",
    "      \n",
    "- FromMeanAndStddev(key: gtsam.Key, A1: np.array, parent1: gtsam.Key, A2: np.array, parent2: gtsam.Key, b: np.array, sigma: float) -> gtsam.GaussianConditional\n",
    "```\n",
    "\n",
    "Both classes support some other methods that can come in handy, as demonstrated here on the prior created above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GaussianDensity: density on [x1] \n",
      "mean: [20; 10];\n",
      "covariance: [\n",
      "\t0.25, 0;\n",
      "\t0, 0.25\n",
      "]\n",
      "isotropic dim=2 sigma=0.5\n",
      "\n",
      "unweighted error vector = [-1.  2.]\n",
      "weighted error vector = [-2.  4.]\n",
      "error 0.5 w'*w = 10.0\n"
     ]
    }
   ],
   "source": [
    "print(p_x1) # printing\n",
    "\n",
    "values = gtsam.VectorValues()\n",
    "values.insert(x[1], [19,12])\n",
    "\n",
    "e = p_x1.unweighted_error(values)\n",
    "print(f\"unweighted error vector = {e}\")\n",
    "\n",
    "w = p_x1.error_vector(values)\n",
    "print(f\"weighted error vector = {w}\")\n",
    "\n",
    "E  = p_x1.error(values)\n",
    "print(f\"error 0.5 w'*w = {E}\")\n",
    "assert E == 0.5 * w.dot(w) # check 10 == (4+16)/2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above, all error functions take an instance of `gtsam.VectorValues`, which is simply a map from GTSAM keys to values as vectors. This is the equivalent of `gtsam.DiscreteValues` from the previous sections."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "include_colab_link": true,
   "name": "S42_logistics_actions.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3.8.12 ('gtbook')",
   "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.12"
  },
  "latex_metadata": {
   "affiliation": "Georgia Institute of Technology",
   "author": "Frank Dellaert and Seth Hutchinson",
   "title": "Introduction to Robotics"
  },
  "vscode": {
   "interpreter": {
    "hash": "9f7376ced4243bb13dfcffa8a3ba834e0602aa8334cd3a1d8ba8d285f4628083"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}