{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import division, print_function\n",
    "%matplotlib inline\n",
    "import sys\n",
    "sys.path.insert(0,'..') # allow us to format the book\n",
    "sys.path.insert(0,'../kf_book') \n",
    "# use same formattibng as rest of book so that the plots are\n",
    "# consistant with that look and feel.\n",
    "import book_format\n",
    "#book_format.load_style(directory='..')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from numpy.random import randn, random, uniform, seed\n",
    "import scipy.stats\n",
    "\n",
    "class ParticleFilter(object):\n",
    "\n",
    "    def __init__(self, N, x_dim, y_dim):\n",
    "        self.particles = np.empty((N, 3))  # x, y, heading\n",
    "        self.N = N\n",
    "        self.x_dim = x_dim\n",
    "        self.y_dim = y_dim\n",
    "\n",
    "        # distribute particles randomly with uniform weight\n",
    "        self.weights = np.empty(N)\n",
    "        self.weights.fill(1./N)\n",
    "        self.particles[:, 0] = uniform(0, x_dim, size=N)\n",
    "        self.particles[:, 1] = uniform(0, y_dim, size=N)\n",
    "        self.particles[:, 2] = uniform(0, 2*np.pi, size=N)\n",
    "\n",
    "\n",
    "    def predict(self, u, std):\n",
    "        \"\"\" move according to control input u with noise std\"\"\"\n",
    "\n",
    "        self.particles[:, 2] += u[0] + randn(self.N) * std[0]\n",
    "        self.particles[:, 2] %= 2 * np.pi\n",
    "\n",
    "        d = u[1] + randn(self.N)\n",
    "        self.particles[:, 0] += np.cos(self.particles[:, 2]) * d\n",
    "        self.particles[:, 1] += np.sin(self.particles[:, 2]) * d\n",
    "\n",
    "        self.particles[:, 0:2] += u + randn(self.N, 2) * std\n",
    "\n",
    "\n",
    "    def weight(self, z, var):\n",
    "        dist = np.sqrt((self.particles[:, 0] - z[0])**2 +\n",
    "                       (self.particles[:, 1] - z[1])**2)\n",
    "\n",
    "        # simplification assumes variance is invariant to world projection\n",
    "        n = scipy.stats.norm(0, np.sqrt(var))\n",
    "        prob = n.pdf(dist)\n",
    "\n",
    "        # particles far from a measurement will give us 0.0 for a probability\n",
    "        # due to floating point limits. Once we hit zero we can never recover,\n",
    "        # so add some small nonzero value to all points.\n",
    "        prob += 1.e-12\n",
    "        self.weights += prob\n",
    "        self.weights /= sum(self.weights) # normalize\n",
    "\n",
    "\n",
    "    def neff(self):\n",
    "        return 1. / np.sum(np.square(self.weights))\n",
    "\n",
    "\n",
    "    def resample(self):\n",
    "        p = np.zeros((self.N, 3))\n",
    "        w = np.zeros(self.N)\n",
    "\n",
    "        cumsum = np.cumsum(self.weights)\n",
    "        for i in range(self.N):\n",
    "            index = np.searchsorted(cumsum, random())\n",
    "            p[i] = self.particles[index]\n",
    "            w[i] = self.weights[index]\n",
    "\n",
    "        self.particles = p\n",
    "        self.weights.fill(1.0 / self.N)\n",
    "\n",
    "\n",
    "    def estimate(self):\n",
    "        \"\"\" returns mean and variance \"\"\"\n",
    "        pos = self.particles[:, 0:2]\n",
    "        mu = np.average(pos, weights=self.weights, axis=0)\n",
    "        var = np.average((pos - mu)**2, weights=self.weights, axis=0)\n",
    "\n",
    "        return mu, var"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAGGCAYAAAB/gCblAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde1xUdf4/8NeZ+8AwDAyODo06ChVUg+WYhV0cu3iNSnK3Vlit/QbkamVMS/eElrL6Sdu2ZBb7yC47pltWu3TBalfd3VorMbOLtV28BdrFG4iKIO/fH/o5nTPMwACD3N7P74PHxsyZcw7Tt8/7fD7vz+f9kYiIwBhjjB2n6ekbYIwx1rtwYGCMMabCgYExxpgKBwbGGGMqHBgYY4ypcGBgjDGmwoGBMcaYCgcGxhhjKhwYGGOMqXBgYLL3338f06dPx7Bhw2A0GjF48GBkZmbC7/erjlu8eDGeeeaZnrnJ45577jlcc801OPXUU6HRaOB2u0Met3HjRkybNg3Dhg2D2WxGYmIiMjMz8Ze//KVL1/f5fJAkKeRPuHsJ54033kBxcXHI99xuN6699tou3WtnLVu2DI8++miPXJv1LIlLYjAAeP3113H55ZfD5/MhLy8PTqcTO3fuxPr167F8+XJ899138rFnnHEGkpKSsGbNmh6730svvRS7du3CmWeeiXXr1qGpqQlbt25tddyaNWuwfPlynH/++TjppJPQ0NCAQCCA5cuX4/e//z3uvvvuTl3f5/Nhx44dCAQCrd4zGo0466yzIj7XvHnz8PjjjyPUf4offfQRrFYrUlJSOnWfXXHZZZfh008/Dfm9sv6NAwMDAIwfPx41NTX44osvoNPpVO+1tLRAo/m5c9kbAoPynjrTgJ177rmora3F9u3bO3V9n8+Hn376CZ9++mmnPq/UVmDoSRwYBi4eSmIAgN27dyMpKalVUACgCgputxufffYZ1q5dG3LopK6uDrfeeitGjBgBg8GAk046CfPnz0dDQ4PqnJIkYd68eXjyySdxyimnwGg04rTTTsPy5csjul/lPXVGuL812g4ePCh/HyaTCYmJiRgzZgxeeOEFAMC1116Lxx9/HABUw1GiMQ4eSlqzZg0kScKyZctw2223wel0wmKxICsrC99//z3q6+uRn5+PpKQkJCUl4brrrsOBAwdU9/T444/jwgsvhMPhQGxsLDweDx5++GE0NTXJx/h8Prz++uvYtm2b6r6EI0eOoLS0FGlpaTAajRg0aBCuu+46/Pjjj930TbITqfv/y2B9QmZmJv785z/jpptuQk5ODkaPHg29Xt/quFdeeQUzZsxAfHw8Fi9eDODY0AlwrBEcP348vvvuO9x5553IyMjAZ599hnvvvReffPIJ3nnnHVXj8ve//x2rV6/Gfffdh9jYWCxevBi/+tWvoNPpMGPGjKj+fS0tLWhpacHevXvx4osvYtWqVSgvL1cd88wzz+C6667D0qVLIx7Xb25ubvWaRqORA1dhYSGef/55lJaW4qyzzkJDQwM+/fRT7N69GwBwzz33oKGhAS+99BL++9//yudwOp1tXvfOO+/EhAkT8Mwzz2Dr1q249dZb5e9u1KhReOGFF/DRRx/hzjvvRFxcHB577DH5s9988w1mzpwpB++PP/4Y999/P7744gs8/fTTAI7lkfLz8/HNN9/glVdeafVdXnHFFfj3v/+NoqIijBs3Dtu2bcOCBQvg8/mwfv16mM3miL4/1ksRY0T0008/0fnnn08ACADp9XoaN24cLVy4kOrr61XHnn766TR+/PhW51i4cCFpNBr68MMPVa+/9NJLBIDeeOMN+TUAZDabadeuXfJrzc3NlJaWRqmpqR2692nTptHw4cPbPKagoED+2wwGAy1evLjVMc8++yxptVp69tln273m+PHj5fMF//zf//2ffNwZZ5xBV155ZZvnmjt3LoX7T3H48OE0e/Zs+ffVq1cTAMrKylIdN3/+fAJAN910k+r1K6+8khITE8Ne++jRo9TU1ETPPfccabVa2rNnj/xeuO/1hRdeIAC0cuVK1esffvghAQj53bK+hYeSGADAbrfj3//+Nz788EM8+OCDuOKKK/C///0Pd9xxBzweD3766ad2z/Haa6/hjDPOwJlnnonm5mb5Z9KkSZAkqVVO4uKLL8bgwYPl37VaLa6++mp8/fXXqmR3NNx555348MMP8frrr+M3v/kN5s2bh0WLFqmOmTVrFpqbmzFr1qyIzpmSkoIPP/yw1c8999wjHzN27Fi8+eabuP3227FmzRocOnQoKn/PZZddpvo9PT0dADBt2rRWr+/Zs0c1nPTRRx/h8ssvh91uh1arhV6vx6xZs3D06FH873//a/far732Gmw2G7KyslT/ns8880wMGTKkR3NPLDp4KImpjBkzBmPGjAEANDU14bbbbsMf/vAHPPzww3j44Yfb/Oz333+Pr7/+OuQQFIBWwWXIkCGtjhGv7d69Gy6XqzN/QkjDhg3DsGHDAABTp04FANxxxx2YPXs2Bg0a1Klzmkwm+bsK57HHHoPL5cKKFSvw0EMPwWQyYdKkSfh//+//4eSTT+7UdQEgMTFR9bvBYGjz9cOHD8NisWD79u244IILcOqpp+KPf/wj3G43TCYTPvjgA8ydOzeiwPX9999j37598rmDRfIQwXo3DgwsLL1ejwULFuAPf/hDRLNvkpKSYDab5XHqUO8r7dq1q9Ux4jW73d6JO47c2LFjsWTJEnz77bedDgyRiI2NRUlJCUpKSvD999/LvYesrCx88cUX3XbdcF599VU0NDTg5ZdfxvDhw+XXN27cGPE5kpKSYLfbUVVVFfL9uLi4Lt8n61kcGBgAYOfOnSETnps3bwYAJCcny68ZjcaQT5aXXXYZHnjgAdjtdowYMaLda/7jH//A999/Lw8nHT16FCtWrEBKSkpUewuhrF69GhqNBiNHjuzW6ygNHjwY1157LT7++GM8+uijOHjwIGJiYuTk/aFDh7o9aSuS/+KaAEBEqKioaHVsW/+ely9fjqNHj+Kcc87pvptlPYYDAwMATJo0CS6XC1lZWUhLS0NLSws2btyIsrIyWCwW3HzzzfKxHo8Hy5cvx4oVKzBy5EiYTCZ4PB7Mnz8fK1euxIUXXohbbrkFGRkZaGlpwfbt2/HWW2/B7/erGpKkpCRcdNFFuOeee+RZSV988UVEU1Y///xzfP755wCO9TIOHjyIl156CQBw2mmn4bTTTgMA5Ofnw2q1YuzYsRg8eDB++uknvPjii1ixYgV+97vfqXoLzz33HH7zm9/g6aefjijPcOjQIaxbty7ke+eeey4A4JxzzsFll12GjIwMJCQkYPPmzXj++eeRmZmJmJgY+fsEgIceeghTpkyBVqtFRkZG2KGarrj00kthMBjwq1/9CkVFRTh8+DCeeOIJ7N27t9WxHo8HL7/8Mp544gl4vV5oNBqMGTMG11xzDQKBAKZOnYqbb74ZY8eOhV6vx3fffYfVq1fjiiuuwPTp06N+7+wE6unsN+sdVqxYQTNnzqSTTz6ZLBYL6fV6GjZsGP3617+mzz//XHXs1q1baeLEiRQXF0cAVDNXDhw4QHfffTedeuqpZDAYKD4+njweD91yyy2qGUgAaO7cubR48WJKSUkhvV5PaWlpFAgEIrrfBQsWhJ0VtGDBAvm4p59+mi644AJKSkoinU5HNpuNxo8fT88//3yrcy5dupQA0NKlS9u9fluzkgBQU1MTERHdfvvtNGbMGEpISCCj0UgjR46kW265hX766Sf5XI2NjXT99dfToEGDSJIkAkBbtmwhovCzkl588cWQ9x48I0x8Tz/++KP8WmVlJY0aNYpMJhOddNJJ9Lvf/Y7efPNNAkCrV6+Wj9uzZw/NmDGDbDabfF9CU1MTLVq0SD6PxWKhtLQ0KigooK+++qrd74/1brzymfUISZIwd+7cVmsJGGM9j6erMsYYU+HAwBhjTIWTz6xH8AgmY71Xh3sM//rXv5CVlYXk5GRIkoRXX31V9T4Robi4GMnJyTCbzfD5fPjss8+idsOMMca6V4cDQ0NDA0aNGhU2afjwww/jkUceQXl5OT788EMMGTIEl156Kerr67t8s4wxxrpfl2YlSZKEV155BVdeeSWAY72F5ORkzJ8/H7fddhsAoLGxEYMHD8ZDDz2EgoKC6Nw1Y4yxbhPVHMOWLVuwa9cuTJw4UX7NaDRi/PjxeO+990IGhsbGRjQ2Nsq/t7S0YM+ePbDb7aoSzYwxxkIjItTX1yM5ObnLe5UAUQ4Mos6NsmKm+H3btm0hP7Nw4UKUlJRE8zYYY2xA2rFjR1TKyXTLrKTgJ30iCvv0f8cdd6CwsFD+ff/+/Rg2bBh27NgBq9XaHbfHGGP9Sl1dHYYOHRq1AoZRDQyiZPKuXbtUBdl++OGHVr0IwWg0qgp6CVarlQMDY4x1QLSG36O6wG3EiBEYMmQI3n77bfm1I0eOYO3atRg3blw0L8UYY6ybdLjHcODAAXz99dfy71u2bMHGjRuRmJiIYcOGYf78+XjggQdw8skn4+STT8YDDzyAmJgYzJw5M6o3zhhjrHt0ODCsX78eEyZMkH8X+YHZs2fjmWeeQVFREQ4dOoTf/va32Lt3L8455xy89dZbvHkHY4z1Eb2uumpdXR3i4+Oxf//+NnMMR48eRVNT0wm8s4HFYDBEZdobY6z7RdpuRqrP1UoiIuzatQv79u3r6Vvp1zQaDUaMGNEtm8Uwxnq3PhcYRFBwOByIiYnhRXDdoKWlBbW1tdi5cyeGDRvG3zFjA0yfCgxHjx6Vg0J3bxY/0A0aNAi1tbVobm6GXq/v6dthjJ1AfWoQWeQUxF65rPuIIaSjR4/28J0MbLW1tSgpKUFtbW1P3wobQPpUYBB4aKP78XfcO1RUVKCyshIVFRU9fStsAOlTQ0nRRETYfWg3Dhw5AIvBAruZi/ax3icvL0/1v4ydCAMuMOw7vA/PbnwWf/rgT/hm7zfy6ykJKbhx7I2YfeZs2Ey2HrxDxn6WnJyMBQsW9PRtsAGmTw4lddaqr1fB9YgLt6y6Bd/u/Vb13rd7v8Utq26B6xEXVn29KurXvvbaayFJEiRJgl6vx8iRI3HrrbeioaEBW7duld9T/uTm5kb9PhhjrD0Dpsew6utVmLZsGogIhNZr+sRrh5oOYdqyaXh95uuYlDopqvcwefJkLF26FE1NTfj3v/+N66+/Hg0NDfKmRu+88w5OP/10+Xiz2RzV6zPGWCQGRI9h3+F9uOqvV4GI0IKWNo9tQQuICFf99SrsOxzdRXRGoxFDhgzB0KFDMXPmTOTk5Kj2zLbb7RgyZIj8Ex8fH9XrM8ZYJAZEYHh247M42HSw3aAgtKAFB5sO4rmPn+vW+zKbzVzWgzHW6/T7wEBE+NMHf+rUZx97/zF0VympDz74AMuWLcPFF18svzZu3DhYLBb556OPPuqWazPGWFv6fY5h96HdqtlHkSIQvtn7DfYc2gN7THRWWb/22muwWCxobm5GU1MTrrjiCvzpT3/CwYMHAQArVqxAenq6fPzQoUOjcl3GGOuIfh8YDhw50KXP1x+pj1pgmDBhAp544gno9XokJyfLpSa2bt0K4FggSE1Njcq1GOvNamtrUVFRgby8PCQnJ/f07bAg/T4wWAyWLn0+zhC9fSRiY2O54WcMP6/oBsDrNHqhfh8Y7GY7UhJS8O3eb0NOUw1HgoSRCSORaE7sxrtjbGDiFd29W79PPkuShBvH3tipz950zk1cJoOxbiBWdPMwUu/U73sMADD7zNm465934VDToYimrGokDcw6M2aNmhW1e3jmmWfCvud2u7tt9hNjjHVUv+8xAIDNZMPKX66EJEnQtPMna6CBBAkvX/0y10xijA1IAyIwAMCk1El4febrMOvNkI7/n5J4zaw3442cNzAxZWIP3SljjPWsARMYgGPB4bvC7/Do5EcxMmGk6r2RCSPx6ORHUVNYw0GBMTagDYgcg5LNZMNN59yEG8feiD2H9qD+SD3iDHFINCdyopkxxjAAA4MgSRLsMfaoLV5jjLH+ov8PJf3wQ89+njHWrXhf7Ojr34Fh/Xrg1FOBsrLOfb6s7Njn16+P7n0x1s3aayz7U2PK+2JHX/8NDOvXA5deCuzbB9x6a8eDQ1nZsc/t23fsPBwcWB/SXmMZ7v2+GDDy8vKQlZXFq6ijqH8Ghh9++DkoCB0JDiIoCCI48LAS66U2bNiACRMmYMOGDQDabyzDva8MGH0lSPAq6m5Avcz+/fsJAO3fv7/Ve4cOHaLPP/+cDh061P6JFi0iAlr/LFrUPZ87ARYsWECjRo06Idfq0HfNepzP56OYmBjy+XxdOk9NTQ0VFxfL/+v1eqm4uDhKd8m6S1vtZmf0zx4DAPj9wKJFrV9vq+cQ3FMQFi06dr4uuPbaayFJEiRJgk6nw7BhwzBnzhzs3bu3S+ftqK1bt0KSJGzcuPGEXpd1r7KyMowdOxZlnc2nHad8+uYhmoGr/wYGoGPBoRuDgjB58mTs3LkTW7duxZ///GdUVlbit7/9bVTOzQa20aNHY/Xq1Rg9enTUhoB4iGbg6t+BAYgsOJyAoAAARqMRQ4YMgcvlwsSJE3H11Vfjrbfekt/fvn07rrjiClgsFlitVvzyl7/E999/3+o8Tz75JIYOHYqYmBj84he/wD5FLqWlpQX33XcfXC4XjEYjzjzzTFRVVcnvjxgxAgBw1llnQZIk+Hy+qP19rHfgWTqsq/p/YADaDg4nnXRCgkKwb7/9FlVVVfIubkSEK6+8Env27MHatWvx9ttv45tvvsHVV1+t+tzXX3+Nv/71r6isrERVVRU2btyIuXPnyu//8Y9/RFlZGRYtWoRNmzZh0qRJuPzyy/HVV18BOLbXNAC888472LlzJ15++eVu+xtZz+gtQ0B9JXnNQohKpiKKopZ8DiVcYvkEJJpnz55NWq2WYmNjyWQyEQACQI888ggREb311luk1Wpp+/bt8mc+++wzAkAffPABER1LPmu1WtqxY4d8zJtvvkkajYZ27txJRETJycl0//33q6599tln029/+1siItqyZQsBoI8++qjN++XkM+sqTl6fOJx87opwPQelbuwpTJgwARs3bsT777+PG2+8EZMmTcKNNx7bRGjz5s0YOnQohg4dKh9/2mmnwWazYfPmzfJrw4YNg8vlkn/PzMxES0sLvvzyS9TV1aG2thbnnXee6rrnnXee6hyMdUZHewC9pefCOm5gBQbgWKMfLpmWnNytw0diz+eMjAw89thjaGxsRElJCYBjQ0mhiviFe10Q7ymPCT6+vXOw3qOrwy/RHL4JPldHcxdtJa95mKl3G3iBoawMCPf/jLW1nS+f0QkLFizAokWLUFtbi9NOOw3bt2/Hjh075Pc///xz7N+/H+np6fJr27dvV/3H9N///hcajQannHIKrFYrkpOT8Z///Ed1nffee08+h8FgAAAcPXq0O/801kldTRxHM/EcfK5o9gA4Qd7LRWVAKor6c47hiiuuaPW61+uluXPnUktLC5111ll0wQUXUHV1Nb3//vvk9Xpp/Pjx8rELFiyg2NhYuuSSS2jjxo30r3/9i0455RS65ppr5GP+8Ic/kNVqpeXLl9MXX3xBt912G+n1evrf//5HRERNTU1kNpuptLSUdu3aRfv27Qt5v5xj6BnKBWYdUV1dTT6fj6qqqjr1+WjeS0+feyCKdo5h4ASGcEEhOfmEBIdwgSEQCJDBYKDt27fTtm3b6PLLL6fY2FiKi4ujX/ziF7Rr1y75WLHyefHixZScnEwmk4mys7Npz5498jFHjx6lkpISOumkk0iv19OoUaPozTffVF2zoqKChg4dShqNRhV4lDgw9C3RWvkscMPdt3Bg6Exj1V6Zi15cBqOncGDo3YIbbtFjqK6ujsr5ozGjiIPLicOzkjoqksVrnSmfwVgPCh6jV6587gqRFM7Kyoo4nxAukdwdeQROWp8Y/XsHt46saBa/Bx8vfu/G2UqMdZRosKORCK6trUVFRQXy8vLkxhw4NjkiEuE+E817bO9aLLr6b2DoTJkLDg6sjxBTQaNB2dh2pjHvjgDQG641oEVlQCqKopJj+P57Iput8zmDUDkHm+3YeQcIzjEMHMpcQHCuoiu5i55c+TzQ8hucfI60sfrwQ3Vw6GgiWRkcbLZj5xtAODD0bt3V8AXPburMbCdxb9XV1T023XWglePgwPD559TQ0BDZyURw6OzsokWLBmRQICI6ePAgB4YeEtzohWoEo9HwheoNVFVVkdvtpqqqqrDHtHeOwsJCcjqdVFhY2Ol7a097fz/3GLqmT+UYDAYDNBoNamtrMWjQIBgMhrZLPZxxBvDxx4DDARw+3PELzp0L/OIXnf98H0VE+PHHHyFJklz9lZ04wQnW4N83bNgg/56VldXp6/j9fnzwwQfw+/1YvXo1AODtt99GY2Mj3n77bUyaNEme7dSRcwSXadmwYQP8fj/Kysq6PGtKaC/XEM0czEDUpwKDRqPBiBEjsHPnzo5NV9uypWsX7urn+yBJkuByuaDVanv6Vgac4EZP/G9WVhZKSkqwatUqbNq0CUajEZWVlZ1ubMvKyuQGGzg2O+m9997D0aNHQUSdOgcAFBYWIi4uTr7vUMGjq7jh714SRfr/ASdIXV0d4uPjsX//flit1pDHEBGam5u53k830uv1HBR6mZKSElRWVmLMmDHYtGkTMjIycO+997YqUqecftqR3ddKSkqwcuVK2O12BAKBqO3c1pkeQ2f/hoEqknazI/pUj0EQQxw8zMH6O2UDqexBtNVYhprrH0njLM5/7rnnIicnB2VlZRgyZEiXG+hww1FtNf68XqGHRSVTEUXRTqIw1pd1JskcKvEa6eyimpoacrvdZDabyefzdTnJrbyX4PsqLi4mj8dDPp+vVZJ4oCWPu2pAz0pibKAJ10CGer2txjR4tlFNTQ0VFBRQZmamakZRcXExpaWlkdvtpurq6i430MXFxZSenk5ut5sKCgpUQaampoZ8Ph95PJ4BM620u3CtJMYGkODNbkStoEceeaRVHaJQtYnE8W+//TbsdjvWrVsnHxsIBLB+/Xr4j6/or62tRX19PaZOnYp3330Xo0ePlq//ySefYMSIEVi1ahU2bNiACRMmYMOGDW3e+4YNG7Bq1Srs27cPu3btwqZNm1T1l5KTkxEIBHDVVVfxSuZehgMDY32AaODLyspQWVkJIlI1sqIRTk9Pl7d4BX4OFsHH5+XlIScnB2PGjEFZWRlqa2uRk5ODt956S97wSZx3woQJuO6667Bt2zbccMMNqllGwfennC3o9/vx8ccfIzY2FgCQk5Mj12MSx4Xb5Y2L5fWsPpl8ZmygqaiowMqVK2GxWDB+/Hj4/X5VYyoa4Z07d8JutwMArFarvM5BmeANlYguKSnB7t27YbfbVcHm4osvxuHDh3HyySdDq9Xi/PPPR25uLh588EHVFFVlsjgrKwt+v18+z5dffonGxkYsWrQImzZtwrJly1BbW4snn3yyzb+Xk889J+o9hubmZtx9990YMWIEzGYzRo4cifvuuw8tLS3RvhRjA0ZeXh7sdjsOHDigeqIXysrKMHbsWCxZsgRZWVmQJAmVlZWorKyUn8jFU/i8efNaPfHn5eXhqquuUk1T9fv9OHToEEwmE5555hlcf/312Lx5M95++234fD4MGTJE9XnRIxE9ioqKCqxevRpLly7F8OHDsWTJEnzyySc4fPgwPvnkk3b/3qysLJx77rkRDVuxKItKpkKhtLSU7HY7vfbaa7RlyxZ68cUXyWKx0KOPPhrR5zn5zJhae7WHRFmKQCAgl6doq4xGQUFBq+PFdURCuqqqSpWcrqmpocLCQvL7/VRYWEher5f8fj8VFxdTaWkpmUwmKi8vV91PqDIa4d4Ll+SO9s50/VWvn5U0bdo0+s1vfqN6LTs7m3JzcyP6PAcGxtRCTRlVNqSi8bTZbK0aUWVQyc3NJZfLJc9MCm50i4uLKSYmhgCQxWJRzRZS3oOY4ZSTk0Ner5d0Oh0BIJPJ1OreI63OGm5abLR3puuven1gWLhwIQ0fPpy+/PJLIiLauHEjORwOWrZsWUSf58DABopIp4K2V0QvVI+B6Fij6na7KT09nXw+H1ksFtLr9XIgqK6uJo/HQ06nk6qqqqimpoacTicBIK1WS5mZmfL6AxFUysvLyWazkcFgIKfTSZmZmVRUVCT3GMS9BgIBcrvd5PF4KCYmhjIzM9v8W3ndQtf0+sDQ0tJCt99+O0mSRDqdjiRJogceeCDs8YcPH6b9+/fLPzt27ODAwAaEriweC/ckXV1dTZmZmZSdnU1Wq5UMBgO53W4qLy8ni8VCw4cPJ6/XS9XV1VRdXU0mk4kAkNvtlj/vcrkoKSlJro5aXFxMsbGxpNFoyGq1kslkIqPRSGazmRwOh+r+xd9ktVoJADkcDvL5fJSfnx/yb+3oegwWWq+vrrpixQr85S9/wbJly3D66adj48aNmD9/PpKTkzF79uxWxy9cuBAlJSXRvg3Ger28vDzU1dWhvr4etbW1HSo5UVlZifr6elURvdraWlx11VXYtm0biAiSJCE+Ph4rV67EVVddhYaGBjQ1NaGmpgbz5s2D0WhES0sLTCYTlixZgmXLlmHu3LkoLS3FSy+9hEcffRT/+te/QEQYNmwYtm7dCofDgcOHD+Pss8/Gu+++i5SUFNUaBPHPmzZtwmuvvYaJEyfi+eefl8tfiGSymBElZh/V19fLhfd4RlIvEJXwoiC6m0q///3v6dRTTw15PPcYWF8VjSdbv99PdrudXC5Xh8bRww0vpaamEgACQJIkyecUeYEpU6aQxWKh/Pz8Vr0Om81GAEij0ZBGo5HPo9frKTMzk3w+HzkcDnI6nXIS2u/3q+5BnDM4eS0E5zXE31FYWEgpKSlks9koEAh0aKMf7mH0gaGkxMREWrx4seq1Bx54gE4++eSIPs85BtZXBM/y6UyCtLCwkAwGA+l0urAzbyJp+MQQktfrJafTSVqtlqZMmdKqcW2rxEZ2djZJkkQAKD4+nrRaLY0ZM0Y1M0mUtZg+fToNHjyYzGYzZcTUPWcAACAASURBVGdnyzkOm81GJpOJMjMzye12k8lkIo/HQ263mwKBgGrWU/B92Ww2kiRJLsch8iPtDbUNtN3aQun1gWH27Nl00kknydNVX375ZUpKSqKioqKIPs+BgfUVopHNzMzs9JTKmpoays/Pb/VkrSQavsLCwpDj8YWFheRyuchsNstJXjGlVNQiClWoLvgaHo+Hhg4dSpIkUVFREVVXV1N6ejpZLBYKBALycfHx8XJvAsdnIwGgmJgY0mq1ZLVaqaCgQK6RZLfbCQAZjUaSJImmTJkSMqgq6zn5fD4ymUzkdru5xxCBXh8Y6urq6Oabb6Zhw4aRyWSikSNH0l133UWNjY0RfZ4DA+trujKlMpJGTRzj9/vlJ2Pla06nk+x2u/ykLT4jnu69Xm+7hepqamrI6/XKjb3VaiWfzyf/HhsbS06nk9LS0kiv18uva7VaMpvNcoAQw1Tl5eVkMpmoqKiIYmNj5aEt5dCUmOkkKrkqv8vgGVasbb0+MHQVBwY2ECjH1iMdBlH2LsQsn8LCQtVYvzJIeTwekiSJXC4XWa1Wys7ODjmEJD6flpYmN/hTpkyhgoICcrvd8uxC0bi73W6SJEmehiqmppaWlspP/CJfodPpyGQykc1mo6KiIjIYDGQ0Gkmr1ZLL5ZLzCsq1FWK2k1hz0ZVpvQMFBwbG+pBwjZUYShGrhyNtzETyNty6AGVyV7kmQTT4MTExciNcXV1NgwcPJo1GQ4mJiarjJ0yYQBaLhTIyMuSgIHIPBoOB9Hq93KAHL7RT5htKS0tVT/4+n4+0Wi1ptVrKzc1tVXZbBBVxzy6XK+LS3O0tBOzPODAw1sspG6NwidFwDVbwUEqoRWttlZSoqqqSh49KS0vlJ3Vl4y5Jknxe5VBRVVUVGY1G1TCRSGQDoKSkJHK5XPL7NpuNfD6fnEsIBALytcMNAYkkeUFBQaseDtHPuY6MjAxyuVyUm5sbUY4k3Hc6UBLTHBgY6+WUjVGoxqqtnERweQtlmYu2FoMph6RE4tZms8lP2oFAQJUwtlgs5PV65fyARqOhQCBAubm5JEkSaTQaeeqo3W4njUZD06dPp6qqqlbTYd1uNxmNRrJYLGSxWMhoNLY7w0rMlsrPzyen0ykvpgv+G5W/B282FPx7W9cLFVT6U2+CAwNjvVxbDU51dbU8pTO48VQWsQvVYxABZ/r06ar5/l6vl3Jzc+VGUjTWqamp5PV65RlPBoNBbtTFTCERFEQPQNyDx+ORZw2J+kkAaOjQoTRt2jQCQLNmzSIiovLyctW6B5vNFraHI+7X5/NRWloaWa1Wio+Pj2gdh8htiFXawb93VH/qTXBgYKwP8/l8ZDQaVY2nUFhYKD89BweXqqoqeWhFlJsQY/yiJpHZbKakpCSyWq1UVFSkqlXk8/movLyc9Ho9JSQkqBpyk8lEsbGxlJSU1KrAnlgboextWK1W0uv1pNFoqLy8XM5N6HQ6crlcZDAY5EWu6enpcm5DrIXIzs4mi8VCdrudTCYTmUwm1cykcDrTY2gL9xjC48DAWDdRNjzKHEC4YSQx9VQkpJU9AYfDIT8di6mg5eXlcgPu9XrlJ2jRSIvXldcTx8TFxckNfUJCgjylVMwEcjqd5PF45MY8MzNTnmo6ffp0ObCYTCbKyMiQ3xNTWUWlVYvFIl9HNPwisMXGxnZoaqpyVXWoVdcDGe/5zFgfodyDWfzzunXrsHr1arm+kVJhYSEKCgowc+ZM1NXVwefzYc2aNdi2bRsAyJvdbNmyBRaLBQ8++CDy8vIwduxYPPXUU1i5ciVOPfVUWCwWFBcXY+zYsSgsLFRdY8mSJRg+fDiuuOIKeDwexMbGYu/evWhoaAAAeDwe3HDDDdi1axfq6+vl+5w0aRKqq6tx1113oby8HI899hhMJhMWLVqEpUuXwmQyobm5GXq9HgDwy1/+ErW1tZg4cSJiYmKQmpqKU045BbW1tbjkkkug1WoxadIkBAIBfPXVVwgEAiG/EwCqDYbWrVuHK6+8Es899xyWLVum2t860r2oWQSiEl6iiHsMrL9Qjq93pIpoqL0PRO7A5/OpSlfodDp5KEV8Lj8/X76ucgqp8npiOmnwKma9Xk+JiYkEgGbMmEEej4csFgulpKTIn50+fTpZrVZ5NTTRsWEdh8OhyjMo/45QCXlxf8HTVYOHh5R/lyj57XK5WvUYBvKmPjyUxFgf0V5yM9zsJWUSOlQBOuWUUQDkdDqJ6OdAI8b1xVCQMnmdlpZGbrdbzguIIaT2frRaLWVnZ1NmZqYclKxWq2p6rUhui+mwoYbSgqeSBk9FVSaUg2cwKafiKld4i/MO5E19ODAw1kcEN4zBNZGC1zsoG0ll7kDM+xcrjIuKilSNtig0J6Z0ivF/u92u6rH4/X65BIUILsq8gE6nk3djC/UjSRKZzWY5eX3RRReRzWaTp6qK4zIyMojo5zUL+fn5qiS6WPMgKrQqg4Xy/eDeRKgeQVvrREIF1/6KAwNjfVBxcTFZLBbS6XSqwnDhhlXE72LrTIvFQpmZmeR0OslqtcprDTIyMlSb4LjdbrmBjouLUxW1E9NPvV6valaSCAqlpaVUXFxMkydPbrPn4HQ66aKLLiKtVksajYZ0Oh0lJyerVk27XC7yeDzyvYuGWww3abVa0uv1Ybf8DdWbUA4zKQNIqF6J+L6VO9b1Z9FuN7XFxcXFUUpXREVjYyMefPBB3HHHHTAajT19O4xFxSmnnII9e/Zg586diImJwZEjR+Dz+fDII4+gsrISdXV1SE5OxsiRIxEbG4tZs2YBAO68807U1NTg0ksvxeDBg/Hll19i/Pjx2L17N2JjY7F48WI0NDRgxYoVuOiii3DRRRfhjTfeAAAcOXIEANDc3IyWlhY0NTVh7dq18kY+Si0tLVizZg2GDh2K999/H6effjpqamrk98eMGQO73Y7du3fjwIED2Lp1K1paWkBE0Ol0qK+vlzcH2rp1K+rq6rB//340NzfD6/XioYceQlxcHMrKytDQ0AAiQktLC7Zu3YrZs2fjvffew6WXXor09HTExMRgzZo1SE1NxaJFi7BmzRpceOGF2LNnD44cOYLY2FiUlpZi27ZtqK6uxiuvvIK4uDj5uzxw4ADy8vKwd+9eaLVa/PGPf4TT6TxB/6Z7RtTbzaiElyjiHgM7UXpiHnu4lb3KyqmijDUUUzyDPyvG4nU6nWqKqLLMhdlsptzcXHkoSlkVtb0f0YMQ01aLiorkXoYkSapzKWsbxcbGqspviJ+MjAyqqqqinJwccjqdNG7cONLr9fLe0U6nU84tBOdexPdhtVrDroLuqX+fvQUPJTEWQmcahROx8jXS+1I2dOXl5WQwGCglJUU1Pl5dXU2pqamk1+tpzJgxcmMcFxdHGo1G3kMhIyODnE6nqtEMzksoh5K0Wm2r2UkAVI20shBfUlISpaenk06nI4fDQVVVVXIJ7dLSUkpISFDNblKeTwwPFRYWyolvjUYjD0UlJyerCvOJ4TGDwUBpaWkDIl/QGRwYGAuhM418dz5htldWO/jamZmZcj5AJItzcnLI4XCQx+ORA4eyoc3OzpZ3ORN7HIRSXV0dNl8gFsqVl5e3en/WrFlyWW4ROIYPH05+v1+espqYmEgajYYSEhLke5QkSV6VnZaWRmazmdLT06m0tJRMJhMNHz6cxo0bR5Ik0eDBg8nr9arKdYi8Qm5uLjmdTrnXE2oqarRmIvX13gYHBsZC6G3/YbdXVjs4kInidVqtlnw+HxUXF6saSxEAxIyguLg4mjx5MtlsNrlE9tSpU8nlclF2drbcmAZXVlWeE4oprcFlL8TPkueW0JXXXknaOC3p9Do52SsS6cpzazQaKi0tJavVSgkJCVRYWEj5+flkMpkoJiaGYmJi5KAheiAajYZ8Ph8VFRWRTqcjo9FI2dnZ8pan4m8PFwA6unYh3Hn6et0kDgyM9QHtBarg9/1+P9ntdvkJuaCggKZPn04Oh0MuiZGbm6uafSNKS4g1BsqS2TqdTjX8A0DeX0H5uiRJ8uwnjUZDw4cPJ02MhnAOCDeBUKz4uQkUPymennruKXI6nTRo0KCQwcTpdJLD4ZDzJsp7Ej2JCRMmEHBsfwfRsPt8PtLr9WSxWORCfsFDYkTqxj3SHoP4zsL1PHrbg0VHcWBgrB9SDj2JZKzYD1mZmCUiuVbSrFmzSKfTkUajkaeNimmsJpOpVRI4MzOTzj77bAJAo0aNIqfTSS6XiywWC5WWllJKSgohBYQ7QVgAwr1BgeHe46/fiWPHtZG8liSJSktLVT2GnJycVkl3Ze2o8vJyMhqNNGXKFMrOziaDwSAHQ+WxXq+31X7Q7QWIUKvC+5Not5sSUdC8tR5WV1eH+Ph47N+/H1artadvh7ETqra2Fo888giICP/973/x0UcfQaPRID4+HkuXLsXXX3+NefPmAQCMRiOuvfZaPP3002hqaoJGo8HEiRPx6aefYteuXWhubpbP63K58P3336OpqSn8xVMA5Bz/57aqqLUc/98AgG/CH2Y0GjFv3jxs2bIF//nPfzBkyBAsXboUAOD3+3HJJZdg0aJFePzxxzFz5kwkJydj586dsFgsiIuLw86dO+FyueDxeFBVVQWbzYbGxkakp6dj9+7dMJvNmDp1KuLi4rBq1Sp8/PHHGDt2LFavXt3qOy0rK4MkSSgsLERycnIbf1zfFO12k4voMdYLiEJxALBo0SKUlZWhvLxcLk6n1Wrx1ltvwe/3y5+RJAn/+Mc/kJCQIL/2z3/+E+PHj8c111wDp9OJzMxMWK1W3HDDDTh69Gj4GzABuPr4P7fXKoj3rz7+OQA6nQ5arRYnnXQStFotgGMF+f72t7/h5Zdfxg8//IBNmzbB7/fD7/fjP//5D+6++27s27cPv/71r7Fq1SoMGjQIwLF1FxdffDHcbjf+/Oc/45///CeICAcOHJALBr777ru4+uqrQUSorKxERkYGxo4di7y8vFaF9CoqKrB27VrExcX1y6DQLaLS74giHkpiA1G45KfYdMdut1NSUhLFxcWphojEHs6ZmZk0ePBg0mq1chmIzMxMOYFrMpmODRWFG/455/gwUXEHfhYc/xxADodD3r/BZrORJEmtEt06nY68Xi8FAoFWK6/FVF1luW8xy2nOnDny7KlgwbmBUHti9/X8QSR4KImxKKutrUVFRQXy8vJ65Ilyw4YNmDdvHjIyMnDvvfe2ugcxvLRixQrU1NSoVi3HxMSgoaEBtbW1uPHGG/HOO+8gJycHFRUVaGlpgV6vl3scVqsVsbGxchlvlZsAJACQOnDjLQD2AXjs2LDR6NGjUV5ejh9//BE5OTnYvXu3fKgkSTjjjDPwzTffID09HT/88AN27Nihel+r1cJiscBsNqOgoADPPPMMtm3bhuHDh2PLli0Rf5d+vx+nnnoq1q9fj6ysLCxYsKADf1TfxENJjEWZct+EE0EMG9XW1gIA5s2bh/Xr12PTpk3YtWtXq6GQ5ORkLFq0CH/7298QGxsrvy5JEqZOnYqEhATk5ubi1VdfRV1dHZ566ik0NzeDiDBu3DhotVocOXIEP/30U6ugYLFYMOL0EUAiOhYUgGOtRyIA87HhnyFDhuCCCy5Abm4uDh8+rDr0pJNOwpQpU3Dw4EFUV1fju+++U71PRGhubkZdXR0uvvhi1NXV4f7775f3oFi1ahVGjBiBVatWtfldjh49GqtXr8a9996LrKws5OXlyceIoHHrrbfKx7PQuMfABrwT3WMoKSlBZWWl/DR7ww03IBAIICcnB19++SU++OADjBo1CpMmTUJeXh527doFv9+PsrIyAMDFF1+M/fv3A4Dce9BoNGhpOZYV1ul0AIDY2FgQEerq6lodAwBarRZarRaDThmEmhk/10XqKO2ftJD2S9DpdHJAsNvt2LNnj3x/DocD+/btk+s3SZKE0aNHY/PmzTh48KDqfCkpKdDr9Th8+DBWrlyJ0aNHY8SIESF7D8HfpWj8y8rK5I1/xDFxcXH48ssvAQAFBQX9qicR9XYzKgNSUcQ5BtbfhauXJGoAiXUAaWlpVFhYSG63W56KajQaac6cOfICMxxfJKbc1hP4uYS2GO8HjtUxcjgc5HA4aPDgwT+P8cd0MLcQ/BMDOY9hNpvJbDbLU2rFNdLT0yk1NVX+XeQLxIZBpaWl8n1qNBq5PLhYBZ2TkyNvO9rWPg+hFrwpK9j21y1BeWtPxvq45ORkLFiwQO6dKH9ft24d6uvrsWfPHhw+fBiSJMFkMuHo0aMgIjQ2Nsr5g9TUVNhsNjz//PNYuXIlxo8fL1+jubkZzc3N2LdvH9LS0hATE3NseueYMdBqtRg9ejQk6fjY0UEAe/DzNNQ2DDqg+KXl2OcccQ5oNBpMnToVt912GxwOB+bPny/3DgBg8+bN0Ol0SE5Ohlarxb59+wAAM2fOxN69e3HXXXdh8uTJACD/vaeddhrKyspQUVGBL774Atdffz0mTZrUauhv586dmDFjBjZs2IAZM2agpaUFM2bMaPV9jx49GmVlZVi0aFHIITumEJXwEkXcY2B9WUdXPId6X6z69Xq9ckG50tJS0uv1ZDQaqbS0lHw+n7y3gnJbTDFrSfm0PmfOHHK73TQGoL0A3W02k9FopMGDB8uLzyKZlVQ4EbTHBPLmtZ6VJHZxy8jIUC1yE/+M4z0B8ZpWq1XtESGe/gsLC8nlcpHRaJQXsAVvSWq1Wik7O1uetRQTE0M6nY58Pp9qB7i2vu/+tg0or3xmrBfryHaebYmk4Qpu9ESjaLfbVdNBJUmicQYD7ZUkIoAIoEIcK5ddU1NzbPMcEwh3gTTFmrBBQXx2jwnkvf74CmjTsWt4vV5VkT+x8jk1NZU0Go28wY8Y/nK5XHLtpNjYWNX3IabopqWltfqebDYb4fi9i21G9Xq9aoqry+VSrZgONw24P62A5qEkxk6g4Fkv7cnLy5Nnw3TmfXG922+/HWPHjpUTzpFYsmQJhg8fjsOHD6OlpUUeKkoiwmtHjsCmmGdSBuC/v/gFPvnkEzQ1NUF/VI+x3449lpwOGlIqfA8oe+vn3xMOA2//BRgUAHAYMBgM8Pl8OO+886DRHGtSYmNjcd1118kb+hw9ehQ//vijvBq7pqYG6enpMJlMSEtLU30fo0ePxrvvvotrrrlGfl18L6WlpbDZbHj88cfh8XggSRKam5thsVgwZMgQrFu3DldffTU2b94sTygI9X2L2UsiQc2CRCW8RBH3GFhvEq2qm6GGNEK91pHribpKhYWFqteVwy0ej4ckSaJC/PzEr/y5NWihWXCtJGVPQfnjt4MmTJggP3WLrTjT0tLIarXKVVhdLhdpNBpKTEykOXPmyENJc+bMkZ/ulTWUQm3AE+57CQQCZDAYyGQyUX5+vnxMYWFhv1/QFoyHkhg7gaK1ajZUwxbqNXG9QCBAbrebysvLww55+P1+cjqd5Pf7iYhozpw5cnVTsVOa2HdZp9PRX889N3QjH7zjmgk0+IrBVHh+6KBQqIU8QyoQCMjF71JSUigtLU0eVtLr9TRhwgSaPHkyWa1WSktLU1V2VRb6S0tLo0AgIA81JSQkkMlkojlz5pDL5aIpU6bIK7oF5d4PwbmI/jZU1B4ODIz1QZH2GATR6InksGjklI1dcEOobNxdLhc5HA7yer3kdDrlUhl3mUyhG3v8XJpbkiQqTUho8zjlj2jMjUYj+Xw+eWc28Z4IBhaLpVWZDPEjGnjxGbF1qAgcYq+K3Nxc+TsK17sg6n/J5fZwYGCsj4ukFyKGWS666CJ5/n6oxk4Ei+CNdkpLSykjI4P0ej1NnTpVDiilpaVhh5VEzyHc+4X4eUtQ5awnrVZLBoOBUlNTqaqqigoKCuTtOOfMmUN2u51MJhPZ7Xa66KKLVHtHiz0lrFYrlZeXk8vlovPOO4+MRiPpdDq5x2A2m+WAF8n3yj2GruHAwNgJFrzZfbjGrLi4WN6bQSzQ8vl8VF5eLj8pi2Dh9XopMzNT3qzH4XCQy+WSG+BAIEDV1dVks9lIo9GEbfy/ayMoKJ/gHQ4HlZaWyvsnhNsAR2wApOwpJCcnqwILANLr9VRcXEzFxcVy4LDZbER0rMHPzs6WcxeiQKCY5hrqew2lPxfT48DAWB8iVvYGAgH5NWUDpWzMgnsEwclZERiUG/cEB5PExEQCQImJiVRVVaVqZH0+H5lMJjnncKtGEzIIBP/cptNRfHw8xcbGysNCCQkJ5HK5KCUlhWJiYigtLY2cTicBoLi4ODk3kpKSQpIk0YgRI1SBQPQAjEajvGud+BvT09PJYrHI35kyQCpXglsslpD5mXANf1/fvrMtHBgY60Xaa4yU8+7b+3xwIx8cKMTvHo9H7jEE724mGmKxBmDOnDmk0WioqKhI1eMQAeN37QQH0VMoLy+nmpoays3NJYfD0apsdnAPQORGEhISSJIkmjx5Mk2bNo0AUGpqqpxcLy0tJbfbLc9iKigooISEBLJarVRVVSX3DrKzs+VeQlpaGrlcLiooKFDtANfVf1d9GQcGxnqR9p5CQ/UYIhUcKJTTUEXjJq4vVgCLISXxmmjAxUpg8Znhw4eTXq+nhIQE2m02hwwKNUFDPcrhq1AJZDF0BYAyMjLIaDTKgTExMVF1L+Jvczgcco5B3LvolbjdbrmXY7PZKC0tjTIzM1W1jgZakjkcDgyM9SLhnkI7kvyM9EnW5/ORXq9XDaEET2+dM2cO6XQ6GjJkiPwUr9Pp5KGd8vJycjqd5Ha7SavVhs01iJ+7zWbS6/Xy0A8AKioqIo/HQ0lJSZSSkiI35MrZSeLY2NhYOWjhePK6tLSUTCYTGY1GcjgcZDAYaPr06XIeRVkwT6yCTk1NJbfbTR6PRxWE25qZNJBwYGCsm0VjyEE86So3rA8n3HqGwsJCKigooIKCAiosLJR3agtOuopzpKWlqWoRiXH/qqoqOQksGu+2Zh8F/zw2fHirRLFyjYUYIlJOUVVWSnW73VRUVCT3nETpDJPJJNeCimSGVnZ2NhUWFnZ6QWB/xoGBsW4WjcZGPOmmp6e3e55wK6CdTifFxsaSwWAgh8PR6jximKq8vJz8fj+5XC7S6XSk0WgoNTVVHmKpqqoig8EgL0prKyjs0ulCvl5yfFqpVquV8wYmkynkugS73U5TpkxRveZwOOQn++AtPNv6DjMzM+UaS2JISdk76M95g47gwMBYN+tqhdT2jovk86LH4PV6yW63U2ZmZqvjxfi9wWAgi8VCU6ZMkRvcQCAgD2WJJ3SdTkeDBg0KGxR2zJ9PPp+PHhw0KOT7b0+eTEajkcaMGUMOh4NSU1NVK5nFzCGRFFcmpINzHT6fj8xmszytNtwiNVFRVq/Xy70QZb6EHcOBgbEe1laPIpJGvyM9krY2pRE9BjFNND4+nrRaLWm1WlUytqqqimw2G6WmplKRVhuy0S86vqpYXOOZjIyws5TEVFmxCC02NjZkcl3cn5hCq9VqVT0Gh8MhD20FN/ai/LhYuS1JktxLaS+fMNAWtxFxYGCsx7XV+AfXL+poKYzg6yh3HAsXUOx2u9zwarVaslqtVF1d3SqovBimVtKtx5/ExVO5mFEUrmfxxPEVzqJHEOoJXjTsolehLGchyn2I8tjBPQblMFxubq6cLBe9kfYoZyoNlKEmDgyM9WLBFU8j6R2Ea7yCF3aF2payvLxcHmIRK4rFtpni2gUFBfR4SkrIRv7xlJRWOQIx9KPRaOilzMyQn9t3zz3y9p3iesH3LmYiiURzexVUBWXiXkxnNRgMrWYkhaPsMQyU5DQHBsZ6sbb2cw4nXOOl7DH4/X65kVcOk4indo1Go1oRLT5fUFBAIy0W2hOicd+cl0eBQICMRiNJkkRnn322nIswH5+mGq58xkGTia4cN45iYmIoMzNT9TeKWURiCEmSJBo6dCg5HA55uCrU9yUadGV+JC0tTf77vF5vh/dr5h5D53BgYKyHicZLPOEGN2LKXclEPSSbzUbV1dWqp3blk7gY309PTz9W9dRioaa4OFWimejnJLDT6SSXyyVXNXU6nWSz2eTqqCVWq/zZPQCNPb6aWmwxmpaWJq9gFjOV7Ha7XOBP9GoMBkPY0uOhFqsVFBTI1Va9Xq987EDMI7SFAwNjfUh7PYjgukmigVWW1rbZbKraSFarVS6jHY6YsWSxWOQG9Ic33qCDJhPtu+ce+bjS0lK5yJ1Wq1XNElJusWmz2ch/PCh4FcNNIpfhdrvJbDbL9wqABg8eLPcERowYQRqNhiZMmBCyxxAIBMjpdFJGRgZVV1fLDb9YgxEbG0u5ubnyd8UrntU4MDDWhwQPE7X1u7KBFQ1ecXExpaSkqObv5+fnk8Viofz8/JDXDK5GqryXSzIyVIFIlL2WJCnkugLRCxE1jUrmzlWtqPZ6varhn0AgIM8kEvWaqqurVYlqZXVU0UtyuVyqYTBlwx9qpzruMahxYGCsD4mkx+D3++UVvW1txhPunMHaylkoh2E8Hg9lZ2fLQzWi0VaO+YuNe7xer+pcyn0ggivCKovi1dTUyENJBoNBXr2t0Wjkz3k8HvJ4PKrZScrvoTN5m4GGAwNj/Ux7M2c6+nQcSeDweDzydE4xZJSSkqKqR6Tc/Cc4MIgAELz1qJiKqpzCKspxpKWlkc/no+zsbIqNjSWv16sKUuECWajZWgNhplFHcGBgrJ9pryHPzMwknU7XZk4hXPAINTU01LFi3F403sonfYvFQrNmzZLLdxOpA4DyfG1dT/Qugmcxib9fFNcLnm4bacAYyDgwMNZHRKsBKygoIIvFQgUFBWGPCZeMDfUEH9zgBg//iIRvWlqa/FQvSlpoNBrVZ8TQkJg1pByGyszMpPz8fPnvDw5IwcNFYhW0yWQK+f1xQAiPAwNjfURnhzxCNYj5+fmUmZkZdjipIz2G4PMHBw+fzycvcBOBzrD3FQAAFXRJREFUpqioiCRJonHjxrVqmEXg8nq9lJ6eTm63m7xeb8hd1pSUwUwkoEWJ8FB4CCk8DgyM9RGdfeINnqlUUFBAVquVTCZTt0zPDA4e1dXVlJKSQgaDQdVIh2qYRc8gNzeX8vPzyeVykclkorS0tFabConjPR4POZ1OVX4ikn0VOtpjGEg9jD4RGL777jvKycmhxMREMpvNNGrUKFq/fn1En+XAwPqrSJ94g9c2xMTEkCRJcg2kaAjVaJaWlpJWq6UJEya0mj5KFLr3IZ76RS9BrMwOnq2kPF6sc3A6nXJg6I51CQOph9HrA8OePXto+PDhdO2119L7779PW7ZsoXfeeYe+/vrriD7PgYH1V515ghVJ4OAFXl0VKm+h3MQnNze31bRTse2m0WhsNa20tLRU3htC5CnEArXg6biix+DxeORgEG4orCtP/dxj6LyoB4bbbruNzj///E5/ngMDY2qigRP7M0TjCVg500mcv6ioSO4xiMZU+dSdn58vP+27XC7V+USeQpQBF1t4+v1+ue5ScPI8kmm4A+mpvyt6fWBIT0+n+fPn04wZM2jQoEF05pln0lNPPRXx5zkwMBZaNJ+AlU/6JpNJruIqZiEFHydmDnk8HtLr9XIJbXGM6AXMmTNH7nWIhXFiq8+2ptueiL+5O87XW/T6wGA0GsloNNIdd9xBGzZsoCVLlpDJZKJnn3025PGHDx+m/fv3yz87duzgwMBYB7SVuG2rIaypqVFVZ5UkSd4W1Gg00pQpU+SqruIcotqrWCEtFscZDAYymUyqUtuxsbHkcDjI6/WSx+NpFXTau7/u0F97INEODBpEWUtLC0aPHo0HHngAZ511FgoKCpCXl4cnnngi5PELFy5EfHy8/DN06NBo3xJjfVZtbS1KSkpQW1sb9pgbbrgB27Ztww033NDqvYqKCixfvhznnXceHn/8cYwYMQKrVq2S3xs8eDBMJhPuu+8+pKSkQKfTYcuWLWhsbERVVRWICD6fD3V1dSgrK8PatWsRFxeHyspKVFZWwu/3w2QyoaWlBY2NjdDr9dDpdCgqKkJBQQF+/etf4+9//zvsdjs2b96M/Px8TJgwARs2bJDvobKyEhUVFd3zBQbJy8tDVlYW8vLyTsj1+qyohBeFYcOG0f/93/+pXlu8eDElJyeHPJ57DIyFF8kTbns9BlGYz2QytdqvIdR6hsTERLnHIKajOhwOKigoUO0IJ6qiKovspaamkiRJ5HQ6W01T9Xq9ZLFYjpUBPz77KLhWFOucXj+U9Ktf/apV8nn+/PkRjy9yjoGxn3VkqCXcsSJPUF5e3uZagVABpri4mBwOB1ksFlV1UyF4mqkooxETE9MqmPl8PnnISVxDlOKIdHc2FlqvDwwffPAB6XQ6uv/+++mrr76iQCBAMTEx9Je//CWiz3NgYANNtMbZ2+pdKPc36Ogisbae6EOVuRAltUMFKLG/g7LsuLKgH+ucXh8YiIgqKyvpjDPOIKPRSGlpaTwribE2dDQhGi6QtBVgfD4fGY1G0ul0lJKSEtG1urrnQaj7iaSsOOu4PhEYuoIDAxtoOto4dmZmjXKzHJvNFrawnVJHViMr93pQrtoW+zQrZzU5nU7KyclpN+iciKDRXwJTtNtNXQ/lvBljxyUnJ2PBggURHy9m1HRkZs3o0aPx6quv4oYbbsCSJUuQnJwMAPD7/fjggw/g9/uxevVq1NbWoqKiAnl5eSgrK4Pf70dZWVm75xezi9asWYP6+nrU1dWhoaEBBoMBBw4cwJo1awAAkiQBANauXYs9e/bI123rnAA69P10xIm4Rl/EgYGxPkTZcIvGPVKTJk3Cli1bVK8FN/7BDWW4RjuYCFJZWVmorKxEfX09/v73vwMAMjMzVVNE4+LicO655+LBBx/E7bffjpKSkpB/T6gAGMnf35HvqDNBdkCISr8jingoibHwIhlG6spMJuXvkVQ8beu8kUxD7eiwWCTH99dFbG3hoSTGBrBInnA7MjwSfKxyWOu8886TF84F9zTas2vXLlRXV6OsrKzNJ/v6+nqMHz8+4if2SP5+7gVEQVTCSxRxj4GxronG2geithfOtSeS7Ujbm6paU1NDubm55HK5OnUPA0m0202JiKing5NSXV0d4uPjsX//flit1p6+Hcb6ja7kJzp6ja+++gp/+9vfMHPmTDz55JNhj83JycHu3btx1VVXterdlJSU4IEHHsCRI0fgdrs73GsZSKLdbka9VhJjrHfqTF2i2tpa+P1+3HrrrW3Wawq+xpAhQ3Drrbe2OZSVnJyMQCCAq666KuSwT15eHn75y1/C5XJhyZIlEd8z6zrOMTA2QIQbexdP+WJGkbJHUVFRgRdeeAHAsdlE7eUslNeIpFfS1lTd5ORkPP/88+2eg0UfDyUxNsCVlJSgsrIScXFxqK+vR1ZWltxY19bWoqysDJIkobCwsNuGoFjXRLvd5B4DYwNc8BqE4B6F1WrtdF7iROQ1WPRxj4ExFlZJSQlWrlwJu92OQCDQ4cZd9EaUvRAWfZx8Zox1WSQbAAHHehN2ux27d+/u1GY6vDFO38SBgbEBKNIZSu3NHGL9E+cYGBuAOrI6uKNF/pS4SF3fxIGBsQGoK419R3B5ir6Jh5IY64cizSF0NxGAeEZS38KBgbF+qDOrnNuzYcMGTJgwARs2bIjaOVnvxENJjPVD3TGEE7ypD+u/uMfAWD/U1SGcUENRZWVlGDt2bEQ7urG+jXsMjLFWQs0mGj16NPcUBggODIyxViIZiuJyF/0XDyUxxlqJZCiqOxLcrHfgHgNjrFPa61Vwj6Lv4h4DY6xT2utVcI+i7+IeA2OsW/Cq576Ly24zxlgfx2W3GWOMdSsODIwxxlQ4MDDGGFPhwMAYY0yFAwNjjDEVDgyMMcZUODAwxhhT4cDAGGNMhQMDY4wxFQ4MjDHGVDgwMMYYU+HAwBhjTIUDA2OMMRUODIwxxlQ4MDDGGFPhwMAYY0yFAwNjjDEVDgyMMcZUODAwxhhT4cDAGGNMhQMDY4wxFQ4MjDHGVDgwMMYYU+HAwBhjTIUDA2OMMRUODIwxxlQ4MDDGGFPhwMAYY0yl2wPDwoULIUkS5s+f392XYowxFgXdGhg+/PBDPPXUU8jIyOjOyzDGGIuibgsMBw4cQE5ODioqKpCQkNBdl2GMMRZl3RYY5s6di2nTpuGSSy7prkswxhjrBrruOOny5ctRXV2N9evXt3tsY2MjGhsb5d/r6uq645YYY4xFKOo9hh07duDmm29GIBCAyWRq9/iFCxciPj5e/hk6dGi0b4kxxlgHSERE0Tzhq6++iunTp0Or1cqvHT16FJIkQaPRoLGxUfVeqB7D0KFDsX//flit1mjeGmOM9Ut1dXWIj4+PWrsZ9aGkiy++GJ988onqteuuuw5paWm47bbbVEEBAIxGI4xGY7RvgzHGWCdFPTDExcXhjDPOUL0WGxsLu93e6nXGGGO9D698ZowxptIts5KCrVmz5kRchjHGWBRwj4ExxpgKBwbGGGMqHBgYY4ypcGBgjDGmwoGBMcaYCgcGxhhjKhwYGGOMqXBgYIwxpsKBgTHGmAoHBsYYYyocGBhjjKlwYGCMMabCgYExxpgKBwbGGGMqHBgYY4ypcGBgjDGmwoGBMcaYCgcGxhhjKhwYGGOMqXBgYIwxpsKBgTHGmAoHBsYYYyocGBhjjKlwYGCMMabCgYExxpgKBwbGGGMqHBgYY4ypcGBgjDGmwoGBMcaYCgcGxhhjKhwYGGOMqXBgYIwxpsKBgTHGmAoHBsYYYyocGBhjjKlwYGCMMabCgYExxpgKBwbGGGMqHBgYY4ypcGBgjDGmwoGBMcaYCgcGxhhjKhwYGGOMqXBgYIwxpsKBgTHGmAoHBsYYYyocGBhjjKlwYGCMMabCgYExxpgKBwbGGGMqHBgYY4ypcGBgjDGmwoGBMcaYCgcGxhhjKhwYGGOMqUQ9MCxcuBBnn3024uLi4HA4cOWVV+LLL7+M9mUYY4x1k6gHhrVr12Lu3LlYt24d3n77bTQ3N2PixIloaGiI9qUYY4x1A4mIqDsv8OOPP8LhcGDt2rW48MIL2z2+rq4O8fHx2L9/P6xWa3feGmOM9QvRbjd1UbinNu3fvx8AkJiYGPL9xsZGNDY2yr/X1dX9//buN7TK8o/j+Od01LMZekJFz5ZzHCGYujTdKsqlgjUoCSQIksxFj4y5Ng6UK4NWoMesfOJqMgkfZKIPmrWKoFG5OSS0sdXQ0CDbJBqjqG1oTdy+vwf7beNi689+v/vsOnbeLzgPzrWz+/7sEq7PzrnvXaY6EgDgL6T04rOZKZFIqKSkRIWFhZO+JplMKhqNjj3y8vJSGQkA8DdS+lFSeXm5Pv74Y7W2tmrx4sWTvmaydwx5eXl8lAQA/9AN81FSRUWFGhsb1dLS8qelIEmRSESRSCRVMQAAUxR4MZiZKioqdOLECZ08eVLxeDzoUwAAUijwYigvL9fRo0f1wQcfaM6cOerp6ZEkRaNRZWdnB306AEDAAr/GEAqFJh0/fPiwnnzyyb/9fm5XBYCpSftrDCn+swgAQIqxVxIAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcFAMAwEExAAAcKSuGt956S/F4XFlZWSoqKtKpU6dSdSoAQIBSUgzHjx9XVVWVdu3apfb2dt1333168MEH1d3dnYrTAQACFDIzC/qgd999t9asWaO6urqxsWXLlmnz5s1KJpN/+b39/f2KRqPq6+vT3Llzg44GAP86Qa+bMwLI5Lh27Zra2tpUXV3tjJeWlur06dMTXj84OKjBwcGx5319fZJGflAAwN8bXS+D+j0/8GL4+eefNTQ0pEWLFjnjixYtUk9Pz4TXJ5NJvfzyyxPG8/Lygo4GAP9qv/zyi6LR6P99nMCLYVQoFHKem9mEMUl6/vnnlUgkxp7/9ttvys/PV3d3dyA/4I2sv79feXl5unz5csZ/rMZcjGMuRjAP4/r6+rRkyRLNmzcvkOMFXgwLFixQOBye8O6gt7d3wrsISYpEIopEIhPGo9Foxv9jj5o7dy5z8V/MxTjmYgTzMO6mm4K5nyjwu5JmzZqloqIiNTU1OeNNTU269957gz4dACBgKfkoKZFI6IknnlBxcbHuuece1dfXq7u7W9u3b0/F6QAAAQrX1NTUBH3QwsJCzZ8/X3v27NHrr7+u33//Xe+8845WrVr1z0KFw9qwYYNmzEjZJZAbBnMxjrkYx1yMYB7GBTkXKfk7BgDAjYu9kgAADooBAOCgGAAADooBAOBIu2Jgu+6RbULuvPNOzZkzRwsXLtTmzZt14cIF37G8SyaTCoVCqqqq8h3Fix9//FFbt27V/PnzNXv2bN1xxx1qa2vzHWvaXb9+XS+++KLi8biys7O1dOlSvfLKKxoeHvYdLeVaWlr08MMPKzc3V6FQSO+//77zdTNTTU2NcnNzlZ2drQ0bNujcuXNTPk9aFQPbdY9obm5WeXm5vvzySzU1Nen69esqLS3VlStXfEfz5uzZs6qvr9fKlSt9R/Hi119/1dq1azVz5kx98sknOn/+vN544w3dcsstvqNNu1dffVUHDx5UbW2tvv32W+3bt0+vvfaaDhw44Dtayl25ckWrVq1SbW3tpF/ft2+f9u/fr9raWp09e1axWEwPPPCABgYGpnYiSyN33XWXbd++3RkrKCiw6upqT4nSQ29vr0my5uZm31G8GBgYsNtuu82ampps/fr1VllZ6TvStNu5c6eVlJT4jpEWNm3aZE899ZQz9sgjj9jWrVs9JfJDkp04cWLs+fDwsMViMdu7d+/Y2B9//GHRaNQOHjw4pWOnzTuG0e26S0tLnfE/2647k4xuRR7UBlk3mvLycm3atEn333+/7yjeNDY2qri4WI8++qgWLlyo1atX69ChQ75jeVFSUqLPPvtMFy9elCR9/fXXam1t1UMPPeQ5mV+XLl1ST0+Ps4ZGIhGtX79+ymto2vy54FS3684UZqZEIqGSkhIVFhb6jjPtjh07pra2Nn311Ve+o3j1/fffq66uTolEQi+88ILOnDmjZ555RpFIRNu2bfMdb1rt3LlTfX19KigoUDgc1tDQkHbv3q0tW7b4jubV6Do52Rra1dU1pWOlTTGM+qfbdWeKHTt26JtvvlFra6vvKNPu8uXLqqys1KeffqqsrCzfcbwaHh5WcXGx9uzZI0lavXq1zp07p7q6uowrhuPHj+vIkSM6evSoVqxYoY6ODlVVVSk3N1dlZWW+43kXxBqaNsUw1e26M0FFRYUaGxvV0tKixYsX+44z7dra2tTb26uioqKxsaGhIbW0tKi2tlaDg4MKh8MeE06fnJwcLV++3BlbtmyZ3nvvPU+J/Hn22WdVXV2txx57TJJ0++23q6urS8lkMqOLIRaLSRp555CTkzM2/r+soWlzjYHtuseZmXbs2KGGhgZ9/vnnisfjviN5sXHjRnV2dqqjo2PsUVxcrMcff1wdHR0ZUwqStHbt2gm3LF+8eFH5+fmeEvlz9erVCf/vQDgczojbVf9KPB5XLBZz1tBr166publ56mtoIJfHA3Ls2DGbOXOmvf3223b+/Hmrqqqym2++2X744Qff0abV008/bdFo1E6ePGk//fTT2OPq1au+o3mXqXclnTlzxmbMmGG7d++27777zt59912bPXu2HTlyxHe0aVdWVma33nqrffTRR3bp0iVraGiwBQsW2HPPPec7WsoNDAxYe3u7tbe3myTbv3+/tbe3W1dXl5mZ7d2716LRqDU0NFhnZ6dt2bLFcnJyrL+/f0rnSatiMDN78803LT8/32bNmmVr1qzJyFs0JU36OHz4sO9o3mVqMZiZffjhh1ZYWGiRSMQKCgqsvr7edyQv+vv7rbKy0pYsWWJZWVm2dOlS27Vrlw0ODvqOlnJffPHFpGtDWVmZmY3csvrSSy9ZLBazSCRi69ats87Ozimfh223AQCOtLnGAABIDxQDAMBBMQAAHBQDAMBBMQAAHBQDAMBBMQAAHBQDAMBBMQAAHBQDAMBBMQAAHBQDAMDxH62sJlCI2Op5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from pf_internal import plot_pf\n",
    "\n",
    "seed(1234)\n",
    "N = 3000\n",
    "pf = ParticleFilter(N, 20, 20)\n",
    "xs = np.linspace (1, 10, 20)\n",
    "ys = np.linspace (1, 10, 20)\n",
    "zxs = xs + randn(20)\n",
    "zys = xs + randn(20)\n",
    "\n",
    "def animatepf(i):\n",
    "    if i == 0:\n",
    "        plot_pf(pf, 10, 10, weights=False)\n",
    "        \n",
    "    idx = int((i-1) / 3)\n",
    "    x, y = xs[idx], ys[idx]\n",
    "    z = [x + randn()*0.2, y + randn()*0.2]\n",
    "\n",
    "    step = (i % 3) + 1\n",
    "    if step == 2:\n",
    "        pf.predict((0.5, 0.5), (0.2, 0.2))\n",
    "        pf.weight(z=z, var=.6)\n",
    "        plot_pf(pf, 10, 10, weights=False)\n",
    "        plt.title('Step {}: Predict'.format(idx+1))\n",
    "    elif step == 3:\n",
    "        pf.resample()\n",
    "        plot_pf(pf, 10, 10, weights=False)\n",
    "        plt.title('Step {}: Resample'.format(idx+1))\n",
    "\n",
    "    else:\n",
    "        mu, var = pf.estimate()\n",
    "        plot_pf(pf, 10, 10, weights=False)\n",
    "        plt.scatter(mu[0], mu[1], color='g', s=100, label='PF')\n",
    "        plt.scatter(x, y, marker='x', color='r', s=180, lw=3, label='Robot')\n",
    "        plt.title('Step {}: Estimate'.format(idx+1))\n",
    "        #plt.scatter(mu[0], mu[1], color='g', s=100, label=\"PF\")\n",
    "        #plt.scatter([x+1], [x+1], marker='x', color='r', s=180, label=\"True\", lw=3)\n",
    "        plt.legend(scatterpoints=1, loc=2)\n",
    "    plt.tight_layout()\n",
    "\n",
    "from gif_animate import animate\n",
    "animate('particle_filter_anim.gif', animatepf, \n",
    "        frames=40, interval=800, figsize=(4, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='particle_filter_anim.gif'>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}