{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <style>\n",
       "        .output_wrapper, .output {\n",
       "            height:auto !important;\n",
       "            max-height:100000px;\n",
       "        }\n",
       "        .output_scroll {\n",
       "            box-shadow:none !important;\n",
       "            webkit-box-shadow:none !important;\n",
       "        }\n",
       "        </style>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import sys\n",
    "sys.path.insert(0,'..') # allow us to format the book\n",
    "\n",
    "# use same formatting as rest of book so that the plots are\n",
    "# consistant with that look and feel.\n",
    "import book_format\n",
    "book_format.set_style()"
   ]
  },
  {
   "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/gCblAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAA9hAAAPYQGoP6dpAABVMElEQVR4nO3de3xT9f0/8FfuaZumaVNCUwMUirNVUpVUNOhG8AIi46ugm0o70X1tCwMVU1e8TYqron6p29dVZbLvcGoRfhvqRLci2wDnHBMKKCreEbEFUW7lWmj7/v0B55CTJr2R0tvr+X30MXNycs4n8evnfT6f9+eiExEBERHRCfquLgAREXUvDAxERKTBwEBERBoMDEREpMHAQEREGgwMRESkwcBAREQaDAxERKTBwEBERBoMDKT6z3/+g4kTJ2LgwIGwWCzo378//H4/iouLNec99dRTePbZZ7umkCc899xzuOGGG3DWWWdBr9cjIyMj4nkbN27E+PHjMXDgQMTFxSElJQV+vx8vvPDCKd0/EAhAp9NF/ItWlmj+8pe/oLS0NOJ7GRkZuPnmm0+prB21aNEi/PrXv+6Se1PX0nFJDAKA119/Hf/1X/+FQCCAgoICuN1ubN++HevWrcPixYvx9ddfq+cOGzYMqampWLVqVZeV94orrsCOHTtw3nnnYc2aNTh27Bi+/PLLZuetWrUKixcvxiWXXIIzzjgDBw8eRGVlJRYvXoxf/vKXuP/++zt0/0AggG3btqGysrLZexaLBeeff36brzVjxgw8+eSTiPSf4oYNG2C325GZmdmhcp6KH/7wh3j//fcj/q7UuzEwEABg1KhRqKmpwUcffQSj0ah5r6mpCXr9ycZldwgMoWXqSAV20UUXoba2Fl999VWH7h8IBPDdd9/h/fff79DnQ7UUGLoSA0Pfxa4kAgDs2rULqampzYICAE1QyMjIwAcffIDVq1dH7Dqpq6vDXXfdhcGDB8NsNuOMM87AzJkzcfDgQc01dTodZsyYgd/+9rf43ve+B4vFgrPPPhuLFy9uU3lDy9QR0b5rrB06dEj9PaxWK1JSUpCbm4sXX3wRAHDzzTfjySefBABNd5RSGYd3Ja1atQo6nQ6LFi3CrFmz4Ha7YbPZMGHCBHzzzTfYv38/CgsLkZqaitTUVNxyyy04cOCApkxPPvkkfvCDH8DlciEhIQFerxePPfYYjh07pp4TCATw+uuvY+vWrZpyKY4ePYqysjJkZWXBYrGgX79+uOWWW/Dtt9920i9Jp1Pn/5dBPYLf78fvfvc73H777cjLy8Pw4cNhMpmanffyyy/juuuuQ1JSEp566ikAx7tOgOOV4KhRo/D111/j3nvvRU5ODj744AM88MAD2LRpE/72t79pKpdXX30VK1euxIMPPoiEhAQ89dRTuPHGG2E0GnHdddfF9Ps1NTWhqakJe/bswR//+EcsX74cFRUVmnOeffZZ3HLLLVi4cGGb+/UbGhqaHdPr9WrgCgaDeP7551FWVobzzz8fBw8exPvvv49du3YBAH7xi1/g4MGD+NOf/oR///vf6jXcbneL97333nsxevRoPPvss/jyyy9x1113qb/dueeeixdffBEbNmzAvffei8TERDzxxBPqZz///HNMnjxZDd7vvvsuHnroIXz00Uf4/e9/D+B4HqmwsBCff/45Xn755Wa/5dVXX41//vOfKCkpwciRI7F161bMnj0bgUAA69atQ1xcXJt+P+qmhEhEvvvuO7nkkksEgAAQk8kkI0eOlLlz58r+/fs1555zzjkyatSoZteYO3eu6PV6Wbt2reb4n/70JwEgf/nLX9RjACQuLk527NihHmtoaJCsrCwZOnRou8o+fvx4GTRoUIvnFBUVqd/NbDbLU0891eycP/zhD2IwGOQPf/hDq/ccNWqUer3wv//+7/9Wzxs2bJhcc801LV5r+vTpEu0/xUGDBsmUKVPU1ytXrhQAMmHCBM15M2fOFABy++23a45fc801kpKSEvXejY2NcuzYMXnuuefEYDDI7t271fei/a4vvviiAJClS5dqjq9du1YARPxtqWdhVxIBAJxOJ/75z39i7dq1eOSRR3D11Vfjk08+wT333AOv14vvvvuu1Wu89tprGDZsGM477zw0NDSof2PHjoVOp2uWk7jsssvQv39/9bXBYMD111+Pzz77TJPsjoV7770Xa9euxeuvv46f/vSnmDFjBubNm6c556abbkJDQwNuuummNl0zMzMTa9eubfb3i1/8Qj1nxIgR+Otf/4q7774bq1atwuHDh2PyfX74wx9qXmdnZwMAxo8f3+z47t27Nd1JGzZswH/913/B6XTCYDDAZDLhpptuQmNjIz755JNW7/3aa6/B4XBgwoQJmn/P5513HtLS0ro090Sxwa4k0sjNzUVubi4A4NixY5g1axZ+9atf4bHHHsNjjz3W4me/+eYbfPbZZxG7oAA0Cy5paWnNzlGO7dq1Cx6PpyNfIaKBAwdi4MCBAICrrroKAHDPPfdgypQp6NevX4euabVa1d8qmieeeAIejwdLlizBo48+CqvVirFjx+J//ud/cOaZZ3bovgCQkpKieW02m1s8fuTIEdhsNnz11Vf4/ve/j7POOgv/+7//i4yMDFitVrzzzjuYPn16mwLXN998g71796rXDteWhwjq3hgYKCqTyYTZs2fjV7/6VZtG36SmpiIuLk7tp470fqgdO3Y0O0c55nQ6O1DithsxYgTmz5+PL774osOBoS0SEhIwZ84czJkzB998843aepgwYQI++uijTrtvNK+88goOHjyIl156CYMGDVKPb9y4sc3XSE1NhdPpRFVVVcT3ExMTT7WY1MUYGAgAsH379ogJz82bNwMA0tPT1WMWiyXik+UPf/hDPPzww3A6nRg8eHCr9/z73/+Ob775Ru1OamxsxJIlS5CZmRnT1kIkK1euhF6vx5AhQzr1PqH69++Pm2++Ge+++y5+/etf49ChQ4iPj1eT94cPH+70pK2S/FfuCQAiggULFjQ7t6V/z4sXL0ZjYyMuvPDCzissdRkGBgIAjB07Fh6PBxMmTEBWVhaampqwceNGlJeXw2az4Y477lDP9Xq9WLx4MZYsWYIhQ4bAarXC6/Vi5syZWLp0KX7wgx/gzjvvRE5ODpqamvDVV1/hjTfeQHFxsaYiSU1NxaWXXopf/OIX6qikjz76qE1DVj/88EN8+OGHAI63Mg4dOoQ//elPAICzzz4bZ599NgCgsLAQdrsdI0aMQP/+/fHdd9/hj3/8I5YsWYKf//znmtbCc889h5/+9Kf4/e9/36Y8w+HDh7FmzZqI71100UUAgAsvvBA//OEPkZOTg+TkZGzevBnPP/88/H4/4uPj1d8TAB599FGMGzcOBoMBOTk5UbtqTsUVV1wBs9mMG2+8ESUlJThy5Aiefvpp7Nmzp9m5Xq8XL730Ep5++mn4fD7o9Xrk5ubihhtuQGVlJa666irccccdGDFiBEwmE77++musXLkSV199NSZOnBjzstNp1NXZb+oelixZIpMnT5YzzzxTbDabmEwmGThwoPzkJz+RDz/8UHPul19+KWPGjJHExEQBoBm5cuDAAbn//vvlrLPOErPZLElJSeL1euXOO+/UjEACINOnT5ennnpKMjMzxWQySVZWllRWVrapvLNnz446Kmj27Nnqeb///e/l+9//vqSmporRaBSHwyGjRo2S559/vtk1Fy5cKABk4cKFrd6/pVFJAOTYsWMiInL33XdLbm6uJCcni8VikSFDhsidd94p3333nXqt+vp6ufXWW6Vfv36i0+kEgGzZskVEoo9K+uMf/xix7OEjwpTf6dtvv1WPLVu2TM4991yxWq1yxhlnyM9//nP561//KgBk5cqV6nm7d++W6667ThwOh1ouxbFjx2TevHnqdWw2m2RlZUlRUZF8+umnrf5+1L1x5jN1CZ1Oh+nTpzebS0BEXY/DVYmISIOBgYiINJh8pi7BHkyi7qvdLYY333wTEyZMQHp6OnQ6HV555RXN+yKC0tJSpKenIy4uDoFAAB988EGsyktERJ2s3YHh4MGDOPfcc6MmDR977DE8/vjjqKiowNq1a5GWloYrrrgC+/fvP+XCEhFR5zulUUk6nQ4vv/wyrrnmGgDHWwvp6emYOXMmZs2aBQCor69H//798eijj6KoqCgmhSYios4T0xzDli1bsGPHDowZM0Y9ZrFYMGrUKLz99tsRA0N9fT3q6+vV101NTdi9ezecTqdmiWYiIopMRLB//36kp6ef8l4lQIwDg7LOTeiKmcrrrVu3RvzM3LlzMWfOnFgWg4ioT9q2bVtMlpPplFFJ4U/6IhL16f+ee+5BMBhUX+/btw8DBw7Etm3bYLfbO6N4RES9Sl1dHQYMGBCzBQxjGhiUJZN37NihWZBt586dzVoRCovFolnQS2G32xkYiIjaIVbd7zGd4DZ48GCkpaVhxYoV6rGjR49i9erVGDlyZCxvRUREnaTdLYYDBw7gs88+U19v2bIFGzduREpKCgYOHIiZM2fi4YcfxplnnokzzzwTDz/8MOLj4zF58uSYFpyIiDpHuwPDunXrMHr0aPW1kh+YMmUKnn32WZSUlODw4cP42c9+hj179uDCCy/EG2+8wc07iIh6iG63umpdXR2SkpKwb9++FnMMjY2NOHbs2GksWd9hMplgMBi6uhhE1EZtrTfbqsetlSQi2LFjB/bu3dvVRenVHA4H0tLSOJeEqA/qcYFBCQoulwvx8fGsuGJMRHDo0CHs3LkTACJu90lEvVuPCgyNjY1qUOjszeL7MmXf4Z07d8LlcrFbiaiP6VH7MSg5BWWvXOo8ym/MPE7Xqq2txZw5c1BbW9vVRaE+pEcFBgW7jzoff+PuYcGCBVi2bBkWLFjQ1UWhPqRHdSXFkohg1+FdOHD0AGxmG5xxXLSPup+CggLN/xKdDn0uMOw9shd/2PgH/Oad3+DzPZ+rxzOTM3HbiNsw5bwpcFgdXVdAohDp6emYPXt2VxeD+pge2ZXUUcs/Ww7P4x7cufxOfLHnC817X+z5AncuvxOexz1Y/tnymN/75ptvhk6ng06ng8lkwpAhQ3DXXXfh4MGD+PLLL9X3Qv/y8/NjXg4iotb0mRbD8s+WY/yi8RARCJrP6VOOHT52GOMXjcfrk1/H2KFjY1qGK6+8EgsXLsSxY8fwz3/+E7feeisOHjyobmr0t7/9Deecc456vjI6iIjodOoTLYa9R/bi2v93LUQETWhq8dwmNEFEcO3/uxZ7j+yNaTksFgvS0tIwYMAATJ48GXl5eZo9s51OJ9LS0tS/pKSkmN6fiKgt+kRg+MPGP+DQsUOtBgVFE5pw6NghPPfuc51arri4OA4HJaJup9cHBhHBb975TYc++8R/nkBnLSX1zjvvYNGiRbjsssvUYyNHjoTNZlP/NmzY0Cn3JiJqSa/PMew6vEsz+qitBILP93yO3Yd3wxkfm1nWr732Gmw2GxoaGnDs2DFcffXV+M1vfoNDhw4BAJYsWYLs7Gz1/AEDBsTkvkRE7dHrA8OBowdO6fP7j+6PWWAYPXo0nn76aZhMJqSnp8NkMgEAvvzySwDHA8HQoUNjci+i7qy2thYLFixAQUEB0tPTu7o4FKbXBwab2XZKn080x24fiYSEBFb8RDg5oxsA52l0Q70+MDjjnMhMzsQXe76IOEw1Gh10GJI8BClxKZ1YOqK+iTO6u7den3zW6XS4bcRtHfrs7RfezmUyiDqBMqOb3UjdU69vMQDAlPOm4L5/3IfDxw63aciqXqdHnDEON517U8zK8Oyzz0Z9LyMjo9NGPxERtVevbzEAgMPqwNIfL4VOp4O+la+shx466PDS9S9xzSQi6pP6RGAAgLFDx+L1ya8jzhQH3Yn/C6UcizPF4S95f8GYzDFdVFIioq7VZwIDcDw4fB38Gr++8tcYkjxE896Q5CH49ZW/Rk2whkGBiPq0PpFjCOWwOnD7hbfjthG3Yffh3dh/dD8SzYlIiUthopmICH0wMCh0Oh2c8c6YTV4jIuoten9X0s6dXft5IupU3Bc79np3YFi3DjjrLKC8vGOfLy8//vl162JbLqJO1lpl2ZsqU+6LHXu9NzCsWwdccQWwdy9w113tDw7l5cc/t3fv8eswOFAP0lplGe39nhgwCgoKMGHCBM6ijqHeGRh27jwZFBTtCQ5KUFAowYHdStRNrV+/HqNHj8b69esBtF5ZRns/NGD0lCDBWdSdQLqZffv2CQDZt29fs/cOHz4sH374oRw+fLj1C82bJwI0/5s3r3M+18m2bNkiAGTDhg2n5X7t+q2pywUCAYmPj5dAIHBK16mpqZHS0lL1f30+n5SWlsaolNRZWqo3O6J3thgAoLgYmDev+fGWWg7hLQXFvHnHr3cKbr75Zuh0Ouh0OhiNRgwcOBDTpk3Dnj17Tum67RUIBDBz5szTek/qfOXl5RgxYgTKO5pPOyH06ZtdNH1X7w0MQPuCQycGBcWVV16J7du348svv8Tvfvc7LFu2DD/72c9icm3q24YPH46VK1di+PDhMesCYhdN39W7AwPQtuBwGoICAFgsFqSlpcHj8WDMmDG4/vrr8cYbbwAAmpqa8OCDD8Lj8cBiseC8885DVVVVs2t89NFHGDlyJKxWK8455xysWrVK8/7q1asxYsQIWCwWuN1u3H333WhoaABwvNWyevVq/O///q/aelE2CaLeg6N06JTFpEMqhmKWYwgXLXeQnn5acgpTpkyRq6++Wn39+eefy9lnny39+/cXEZHHH39c7Ha7vPjii/LRRx9JSUmJmEwm+eSTT0TkZI7B4/HIn/70J/nwww/l1ltvlcTERPnuu+9EROTrr7+W+Ph4+dnPfiabN2+Wl19+WVJTU2X27NkiIrJ3717x+/1SUFAg27dvl+3bt0tDQ0PE8jLH0HOF5glYjr4h1jmGvhMYRKIHh9OQaJ4yZYoYDAZJSEgQq9UqAASAPP744yIikp6eLg899JDmMxdccIH87Gc/E5GTgeGRRx5R3z927Jh4PB559NFHRUTk3nvvlbPOOkuamprUc5588kmx2WzS2NgoIiKjRo2SO+64o9XyMjDQqWLy+vRh8vlUROtWChXj7qNQo0ePxsaNG/Gf//wHt912G8aOHYvbbrsNdXV1qK2txcUXX6w5/+KLL8bmzZs1x/x+v/rPRqMRubm56jmbN2+G3+/XrPl08cUX48CBA/j666875TtR39He3AWT1z1X3woMwPFKP1oyLT2904ICcHLP55ycHDzxxBOor6/HnDlz1PfDF/ETkTYt7KecE+l8ObEBEBcI7BlONXEcy7kH4ddqb+6ipeR1T5kj0Vf1vcBQXg5E+3/G2tqOL5/RAbNnz8a8efNw4MABpKen46233tK8//bbbyM7O1tzbM2aNeo/NzQ0oLq6GllZWQCAs88+G2+//bZmN7i3334biYmJOOOMMwAAZrMZjY2NnfWV6BSdauI4lonn8GvFsgXABHk3F5MOqRjqzTmG0OSzwufzyfTp0+VXv/qV2O12Wbx4sXz00Ucya9asiMnngQMHyksvvSSbN2+WwsJCsdls8u2334rIyeTz9OnTZfPmzfLKK69oks8iIgUFBXLBBRfIli1b5Ntvv1VzD+GYY+gaHU3YVldXSyAQkKqqqpglfDszeczEdGwx+dxLRiUpKisrxWw2y5dffilz5syRM844Q0wmk5x77rny17/+VT1PCQyLFi2SCy+8UMxms2RnZ8vf//53zfVWrVolF1xwgZjNZklLS5NZs2bJsWPH1Pc//vhjueiiiyQuLk4AyJYtWyKWl4GhZ4nVzGcFK+6ehYGhI5VVa8tcdNNlMLoSA0P3Fl5xKy2G6urqmFw/FiOKGFxOH45Kaq+2TF7ryPIZRF0ovI8+dObzqVCSwhMmTGhzPiFaIrkz8ghMWp8evXsHt/bMaFZeh5+vvO7E0UpE7aVU2LFIBNfW1mLBggUoKChQK3Pg+OCItoj2mViWsbV7UWz13sDQkWUuGByoh1CGgsZCaGXbkcq8MwJAd7hXnxaTDqkYikmO4ZtvRByOjucMIuUcHI7j1+0jmGPoO0JzAeG5ilPJXXTlzOe+lt9g8rmtldXatdrg0N5EcmhwcDiOX68PYWDo3jqr4gsf3dSR0U5K2aqrq7tsuGtfW46DgeHDD+XgwYNtu5gSHDo6umjevD4ZFEREDh48yMDQRcIrvUiVYCwqvkitgaqqKsnIyJCqqqqo57R2jWAwKG63W4LBYIfL1prWvj9bDKemR+UYzGYz9Ho9amtr0a9fP5jN5paXehg2DHj3XcDlAo4caf8Np08HfvSjjn++BxIRHD16FN9++y30ej3MZnNXF6nPCU+whr9ev369+nrChAkdvk9xcTHeeecdFBcXY+XKlQCAFStWoL6+HitWrMDYsWPV0U7tuYby36Tyv+vXr0dxcTHKy8tPedSUorVcQyxzMH1RjwoMer0egwcPxvbt29s3XG3LllO78al+vgeKj4/HwIEDodf3/hHN3U14paf874QJEzBnzhwsX74c7733HiwWC5YtW9bhyra8vFytsIHjo5PefvttNDY2apZVac81ACAYDCIxMVEtd6TgcapY8XcunbT1/wNOk7q6OiQlJWHfvn2w2+0RzxERNDQ0cM2fTmIwGGA0GrnwXjczZ84cLFu2DLm5uXjvvfeQk5ODBx54oNkidaHDT9uz+9qcOXOwdOlSOJ1OVFZWxmznto60GDr6HfqqttSb7dGjWgwKnU4Hk8kEk8nU1UUh6lShFWRoC6KlyjLSWP+2VM7K9S+66CLk5eWhvLwcaWlpp1xBR+uOaqny53yFLhaTTEUMxTqJQtSTdSTJHCnx2tbRRTU1NZKRkSFxcXESCAROOckdWpbwcpWWlorX65VAINAsSdzXksenqk+PSiLqa6JVkJGOt1SZho82qqmpkaKiIvH7/ZoRRaWlpZKVlSUZGRlSXV19yhV0aWmpZGdnS0ZGhhQVFWmCTE1NjQQCAfF6vX1mWGln4VpJRH1I+GY3ylpBjz/+eLN1iCKtTaScv2LFCjidTnU/jwULFqCyshLr1q1D8YkZ/bW1tdi/fz+uuuoq/Otf/8Lw4cPV+2/atAmDBw/G8uXLsX79eowePRrr169vsezr16/H8uXLsXfvXuzYsQPvvfeeZv2l9PR0VFZW4tprr+VM5m6GgYGoB1Aq+PLycixbtgwioqlklUo4Oztb3SoWOBksws8vKChAXl4ecnNzUV5ejtraWuTl5eGNN96A3W5XA5ESBG655RZs3boVU6dO1YwyCi1b+EjB4uJivPvuu0hISAAA5OXlqesxKedG2+WNi+V1rR6ZfCbqaxYsWIClS5fCZrNh1KhRKC4u1lSmSiW8fft2OJ1OAIDdblfnOYQmeCMloufMmYNdu3bB6XRqgs1ll12GI0eO4Mwzz4TBYMAll1yC/Px8PPLII+oQ1UjzLIqLi9XrfPzxx6ivr8e8efPw3nvvYdGiRaitrcVvf/vbFr8vk89dJ+YthoaGBtx///0YPHgw4uLiMGTIEDz44INoamqK9a2I+oyCggI4nU4cOHBA80SvKC8vx4gRIzB//nxMmDABOp0Oy5Ytw7Jly9QncuUpfMaMGZonfuX61157rWaYanFxMQ4fPgyr1Ypnn30Wt956KzZv3owVK1YgEAggLS1N/Wxoa0RpUSxYsAArV67EwoULMWjQIMyfPx+bNm3CkSNHsGnTpla/74QJE3DRRRe1qduKYiwmmYoQZWVl4nQ65bXXXpMtW7bIH//4R7HZbPLrX/+6TZ9n8plIq7W1h5RlKSorK9XlKVpaRqOoqKjZ+cp9lIR0VVWVJjldU1MjwWBQiouLJRgMis/nk+LiYiktLZWysjKxWq1SUVGhKU+kZTSivRctyR3rnel6q24/Kmn8+PHy05/+VHNs0qRJkp+f36bPMzAQaUUaMhpakSqVp8PhaFaJhgaV/Px88Xg86sik8Eq3tLRU4uPjBYDYbDbNaKHQMigjnPLy8sTn84nRaBQAYrVam5W9rauzRhsWG+ud6Xqrbh8Y5s6dK4MGDZKPP/5YREQ2btwoLpdLFi1a1KbPMzBQX9HWoaCtLaIXqcUgcrxSzcjIkOzsbAkEAmKz2cRkMqmBoLq6Wrxer7jdbqmqqpKamhpxu90CQAwGg/j9fnX+gRJUKioqxOFwiNlsFrfbLX6/X0pKStQWg1LWyspKycjIEK/XK/Hx8eL3+1v8rpy3cGq6fWBoamqSu+++W3Q6nRiNRtHpdPLwww9HPf/IkSOyb98+9W/btm0MDNQnnMrksZqaGrVbJ7Qyra6uFr/fL5MmTRK73S5ms1kyMjKkoqJCbDabDBo0SHw+n1RXV0t1dbVYrVYBIBkZGernPR6PpKamqqujlpaWSkJCguj1erHb7WK1WsVisUhcXJy4XC5N+ZXvZLfbBYC4XC4JBAJSWFgY8bu2dz4GRdbtA8OLL74oHo9HXnzxRXnvvffkueeek5SUFHn22Wcjnj979mwB0OyPgYF6u9B++45UguGBRZm1rNPpBIDodDpxOBxqy0Gn04nFYhGj0Sh+v18CgYCYzWaxWq1SVVUllZWV4nA4pKKiQgKBgOj1esnNzRWfzyfZ2dkSFxcnQ4cOFY/HIxMnThSXy6W2KkK/U2lpqUyaNEnMZrPahawcr6qq0rRqlO+g5CuU8/rSXgqx0O0Dg9LcDPXLX/5SzjrrrIjns8VAPVUsnmyLi4vF6XSKx+Npdz96pCUmhg4dqj5c6XQ69ZpKXmDcuHFis9mksLCwWf+9w+EQAKLX60Wv16vXMZlMaiBxuVzidrvV1kp4UFOuGZ68VoTnNZTvEAwGJTMzUxwOh1RWVrZrox+2MHpAYEhJSZGnnnpKc+zhhx+WM888s02fZ46BeorwUT4dSZAGg0Exm81iNBqjjrxpS8WndCH5fD5xu91iMBhk3LhxzSrXlpbYmDRpktraSEpKEoPBILm5uZqRScqyFhMnTpT+/ftLXFycTJo0Sc1xOBwOsVqt4vf7JSMjQ6xWq3i9XsnIyJDKykrNqKfwcjkcDtHpdOpyHEp+pLWWA1sYPSAwTJkyRc444wx1uOpLL70kqampUlJS0qbPMzBQT6FUsn6/v8NDKmtqaqSwsLDZk3UopeILBoMR++ODwaB4PB6Ji4tTk7zKkFJlLaJIC9WF38Pr9cqAAQNEp9NJSUmJVFdXS3Z2tthsNqmsrFTPS0pK0nT7KnmK+Ph4MRgMYrfbpaioSF0jyel0CgCxWCyi0+lk3LhxEYNq6HpOgUBArFarZGRksMXQBt0+MNTV1ckdd9whAwcOFKvVKkOGDJH77rtP6uvr2/R5BgbqaU5lSGVbKjXlnOLiYvXJOPSY2+0Wp9OpPmkrn1Ge7n0+X6sL1dXU1IjP51Mre7vdLoFAQH2dkJAgbrdbsrKyxGQyqccNBoPExcWpAULppqqoqBCr1SolJSWSkJCgdm2Fdk0pI52UlVxDf8vwEVbUsm4fGE4VAwP1BaF9623tBgltXSijfILBoKavPzRIeb1e0el04vF4xG63y6RJkyJ2ISmfz8rKUiv8cePGSVFRkWRkZKijC5XKXUlkK8NQlaGpZWVl6hO/kq8wGo1itVrF4XBISUmJmM1msVgsYjAYxOPxqHmF0LkVymgnZc7FqQzr7SsYGIh6kGiVVaTROG2hJG+jzQsITe6GzklQKvz4+Hi1Eq6urpb+/fuLXq+XlJQUzfmjR48Wm80mOTk5alBQcg9ms1lMJpNaoYdPtAvNN5SVlWme/AOBgBgMBjEYDJKfn99s2W0lqChl9ng8bV6au7WJgL0ZAwNRNxdaGUVLjEarsMK7UiJNWmtpSYmqqiq1+6isrEx9Ug+t3HU6nXrd0K6iqqoqsVgsmm4iJZENQFJTU8Xj8ajvOxwOCQQCai6hsrJSvXe0LiAlSV5UVNSshSNyMteRk5MjHo9H8vPz25Qjifab9pXENAMDUTcXWhlFqqxaykmEL28RusxFS5PBQruklMStw+FQn7QrKys1CWObzSY+n0/ND+j1eqmsrJT8/HzR6XSi1+vVoaNOp1P0er1MnDhRqqqqmg2HzcjIEIvFIjabTWw2m1gsllZHWCmjpQoLC8XtdquT6cK/Y+jr8M2Gwl+3dL9IQaU3tSYYGIi6uZYqnOrqanVIZ3jlGbqIXaQWgxJwJk6cqBnv7/P5JD8/X60klcp66NCh4vP51BFPZrNZrdSVkUJKUFBaAEoZvF6vOmpIWT8JgAwYMEDGjx8vAOSmm24SEZGKigrNvAdlUp3yfSNNaAsEApKVlSV2u12SkpLaNI9DyW0os7TDX7dXb2pNMDAQ9WCBQEAsFoum8lQEg0H16Tk8uFRVValdK8pyE0ofv7ImUVxcnKSmpordbpeSkhLNWkWBQEAqKirEZDJJcnKypiK3Wq2SkJAgqampzRbYU+ZGhLY27Ha7mEwm0ev1UlFRoeYmjEajeDweMZvN6iTX7OxsNbehzIWYNGmS2Gw2cTqdYrVaxWq1akYmRdORFkNL2GKIjoGBqJOEVjzRloQIpQw9VRLSoS0Bl8ulPh0rQ0GVpSvi4+PF5/NplsMwGo3q8dD7KeckJiaqFX1ycrI6pFQZCeR2u8Xr9aqVud/vV4eaTpw4UQ0sVqtVcnJy1PeUoazKSqs2m029j1LxK4EtISGhXUNTQ2dVn8pSIr0R93wm6iFC92BW/nnNmjVYuXKlunNaqGAwiKKiIkyePBl1dXUIBAJYtWoVtm7dCgDqZjdbtmyBzWbDI488goKCAowYMQLPPPMMli5dirPOOgs2mw2lpaUYMWIEgsGg5h7z58/HoEGDcPXVV8Pr9SIhIQF79uzBwYMHAQBerxdTp07Fjh07sH//frWcY8eORXV1Ne677z5UVFTgiSeegNVqxbx587Bw4UJYrVY0NDTAZDIBAH784x+jtrYWY8aMQXx8PIYOHYrvfe97qK2txeWXXw6DwYCxY8eisrISn376KSorKyP+JgA0GwytWbMG11xzDZ577jksWrRIs791W/eipjaISXiJIbYYqLcI7V9vzyqikfY+UHIHgUBAs3SF0WhUu1KUzxUWFqr3DR1CGno/ZThp+Cxmk8kkKSkpAkCuu+468Xq9YrPZJDMzU/3sxIkTxW63q7OhRY5367hcLk2eIfR7RErIK+ULH64a3j0U+r2UJb89Hk+zFkNf3tSHXUlEPURryc1oo5dCk9CRFqALHTIKQNxut4icDDRKv77SFRSavM7KypKMjAw1L6B0IbX2ZzAYZNKkSeL3+9WgZLfbNcNrleS2Mhw2Ulda+FDS8KGooQnl8BFMoUNxQ2d4K9fty5v6MDAQ9RDhFWP4mkjh8x1CK8nQ3IEy7l+ZYVxSUqKptJWF5pQhnUr/v9Pp1LRYiouL1SUolOASmhcwGo3qbmyR/nQ6ncTFxanJ60svvVQcDoc6VFU5LycnR0ROzlkoLCzUJNGVOQ+R9pMIfT+8NRGpRdDSPJFIwbW3YmAg6oFKS0vFZrOJ0WhstjCcUmn7/X61IlSefpWtM202m/j9fnG73WK329W5Bjk5OZpNcDIyMtQKOjExUbOonTL81OfzaUYlKUGhrKxMSktL5corr2yx5eB2u+XSSy8Vg8Eger1ejEajpKena2ZNezwe8Xq9atmVilvpbjIYDGIymaJu+RupNRHazRQaQCK1SpTfO3THut6MyWeiHqigoACTJ09GWloarFarJmm6YMECvPHGGwCAK664AkajEX6/H9u2bUMwGMQFF1yAvLw85OTkAAAuv/xy9OvXD/369cNjjz2GgQMH4t1334XRaMRdd92lXnf//v1oamrCwYMHsW7dOhw7dgwvvPACqqur0dTUpClfQ0MDHnjgAXzxxRd4//334ff7odPp1Pdzc3Ph9Xqh1+uxc+dOrFq1Co2NjWhqaoJer8c333wDABARrFq1Cl9//TU++eQTNDY2wuv1oqCgAADUazY2NuLYsWN49dVXUVtbi+XLl2Pw4MFYvnw5amtrUVdXhzFjxqCyshKrVq1CcnIynnnmGTidTqxZswZTp07F1q1bcd9992H27NlIT0/XJPsLCgqQl5eH3NxclJeXd8K/0V4uJuElhthioNOlK8axR0tCh3abKMtYI2SIZ/hnlb54o9GoGSIausxFXFyc5Ofnq11RoauitvantCCUYaslJSVqK0On02muFbq2UUJCgmb5DeUvJydHqqqqJC8vT9xut4wcOVJMJpO6d7Tb7VZzC+G5F+X3sNvtUWdBd9W/z+6CXUlEUbS3YjgdM1/bWqbQiq6iokLMZrNkZmZq+serq6tl6NChYjKZJDc3V62MExMTRa/Xq3so5OTkiNvt1lSa4XmJ0K4kg8HQbHQSAE0lHboQX2pqqmRnZ4vRaBSXyyVVVVXqEtplZWWSnJysGd0Uej2leygYDKqJb71er3ZFpaenaxbmU7rHzGazZGVl9Yl8QUcwMBBF0d6KvjOfMFtbVjv83n6/X80HKMnivLw8cblc4vV61cARWtFOmjRJ3eVM2eMgkurq6qj5AmWiXEVFRbP3b7rpJnVZbiVwDBo0SIqLi9UhqykpKaLX6yU5OVkto7LXdFVVlWRlZUlcXJxkZ2dLWVmZWK1WGTRokIwcOVJ0Op30799ffD6fZrkOJa+Qn58vbrdbbfVEGooaq5FIPb21wcBAFEV3+o+7tWW1w4OYsnidwWCQQCAgpaWlmspSCQDKiKDExES58sorxeFwqEtkX3XVVeLxeGTSpElqZRq+smroNREypDV82Qvlb/5z8+Wam68RQ6JBjCajmuxVEumh19br9VJWViZ2u12Sk5MlGAxKYWGhWK1WiY+Pl/j4eDVoKC0QvV4vgUBASkpKxGg0isVikUmTJqlbnirfPVoAaO/chWjX6enrJjEwEPUArQWp8PeLi4vF6XSqT8hFRUUyceJEcblc6pIY+fn5mtE3ytISyhyD0CWzjUajpvsHgLq/QuhxnU6njn7S6/UyaNAg0cfrBRdCcDsEpSF/t0OSxibJM889I263W/r16xcxmLjdbnG5XOqOc6FlUloSo0ePFuD4/g5KxR4IBMRkMonNZlMX8gvvEhPRVu5tbTEov1m0lkd3eqjoCAYGol4otOtJScYq+yGHJmZFRF0r6aabbhKj0Sh6vV4dNqoMY7Varc2SwH6/Xy644AIBIOeee6643W7xeDxis9mkrKxMMjMzBZkQ3AvBbAgeCAsMD5w4fi+On9dC8lqn00lZWZmmxZCXl9dsOe3QtaMqKirEYrHIuHHjZNKkSWI2m9VgGHquz+drth90awEi0qzw3iTW9aZORCTCYKUuU1dXh6SkJOzbtw92u72ri0N0WtXW1uLxxx+HiODf//43NmzYAL1ej6SkJCxcuBCfffYZZsyYAQCwWCy4+eab8fvf/x7Hjh2DXq/HmDFj8P7772PHjh1oaGhQr+vxePDNN9/g2LFj0W+eCSDvxD+3NJBdGelaCeDz6KdZLBbMmDED69evx9tvv43y8nJMnz4d69evR3FxMS6//HLMmzcPTz75JCZPnoz09HRs374dNpsNiYmJ2L59OzweD7xeL6qqquBwOFBfX4/s7Gzs2rULcXFxuOqqq5CYmIjly5fj3XffxYgRI7By5cpmv2l5eTl0Oh2CwSDS09Nb+HI9U6zrTc5jIOoGlIXiAGDevHkoLy9HRUWFujidwWDAG2+8geLiYvUzOp0Of//735GcnKwe+8c//oFRo0bhhhtugNvtht/vh91ux9SpU9HY2Bi9AFYA15/459ZqBeX96098DoDRaITBYMAZZ5wBg8EA4PiCfH/+85+xcuVK1NfX48477wQAFBcX46233sL999+PvXv34ic/+QmWL1+Ofv36ATg+p+Kyyy5DRkYGfve73+Ef//gHRAQHDhxQFwz817/+heuvvx4igmXLliEnJwcjRoxAQUFBs4X0FixYgNWrVyMxMbFXBoVOEZN2RwyxK4n6omjJT2XTHafTKampqZKYmKjpIlL2cPb7/dK/f38xGAzqMhB+v19N4Fqt1uNdRdG6fy480U1U2o6/2Sc+B4jL5VL3b3A4HM3mOuDEsFplXaXwmdfKUN3Q5b6VUU7Tpk1TR0+FC88NRNoTu6fnD9qCXUlEMVZbW6vOlu2KJ8r169djxowZyMnJwQMPPNCsDEr30pIlS1BTU4PQ/2Tj4+Nx8OBB1NbW4rbbbsPf/vY35OXlYcGCBWhqaoLJZFJbHHa7HQkJCeoy3hq3A0gGoGv+VlRNAPYCeOJ4t9Hw4cNRUVGBb7/9Fnl5edi1a5d6ql6vxznnnIPPP/8c2dnZ2LlzJ7Zt26a+r9PpYDAYYLPZEBcXh6KiIjz77LPYunUrBg0ahC1btrSpSEo31VlnnYV169ZhwoQJmD17dju+VM/EriSiGAtdSuF0ULqNamtrAQAzZszAunXr8N5772HHjh3NukLS09Mxb948/PnPf0ZCQoJ6XKfT4aqrrkJycjLy8/PxyiuvoK6uDs888wwaGhogIhg5ciQMBgOOHj2K7777rllQsNlsGHzOYCAF7QsKwPHaIwVA3PHun7S0NHz/+99Hfn4+jhw5ojl14MCBGDduHA4dOoTq6mp8/fXXmvdFBA0NDairq8Nll12Guro6PPTQQ+oeFKFLZrT0Ww4fPhwrV67EAw88gAkTJqCgoEA9Rwkad911l3o+RcYWA/V5p7vFMGfOHCxbtkx9mp06dSoqKyuRl5eHjz/+GO+88w7OPfdcjB07FgUFBdixYweKi4vVNX8uu+wy7Nu3DwDU1oNer1fXPzIajQCAhIQEiAjq6uqanQMABoMBBoMB/b7XDzXX1XT4+xh+Y4Bunw5Go1ENCE6nE7t371bL53K5sHfvXhw9ehTA8aA2fPhwbN68GYcOHdJcLzMzEyaTCUeOHMHSpUsxfPhwDB48OGLrIfy3VCr/8vJydeMf5ZzExER8/PHHAICioqJe1ZKIeb0Zkw6pGGKOgXq78D7v0NdKP7vL5ZKsrCwJBoOSkZGhDkW1WCwybdo0dYIZTkwSC93WEzi5hLbS3w8cX8fI5XKJy+WS/v37n+zjj29nbiH8Lx5qHiMuLk7i4uLUIbXKPbKzs2Xo0KHqayVfoGwYVFZWppZTr9ery4Mrs6Dz8vLUbUdb2uch0oS30H0deuuWoFxdlaiHS09PV1cEDX+9Zs0a7N+/H7t378aRI0eg0+lgtVrR2NgIEUF9fb2aPxg6dCgcDgeef/55LF26FKNGjVLv0dDQgIaGBuzduxdZWVmIj48/PrwzNxcGgwHDhw8/uXrqIQC7cXIYagv6HQh50XT8c65EF/R6Pa666irMmjULLpcLM2fOVFsHALB582YYjUakp6fDYDBg7969AIDJkydjz549uO+++3DllVcCgPp9zz77bJSXl2PBggX46KOPcOutt2Ls2LHNuv62b9+O6667DuvXr8d1112HpqYmXHfddc1+7+HDh6O8vBzz5s2L2GVHIWISXmKILQbqydo74znS+8qsX5/Ppy4oV1ZWJiaTSSwWi5SVlUkgEFD3VgjdFlMZtRT6tD5t2jTJyMiQXED2AHJ/XJxYLBbp37+/OvmsLaOSgmMgu60QX0HzUUnKaKOcnBzNJDfln3GiJaAcMxgMmj0ilKf/YDAoHo9HLBaLOoEtfEtSu90ukyZNUkctxcfHi9FolEAgoNkBrqXfu7dtA8qZz0TdWHu282xJWyqu8EpPqRSdTqdmOKhOp5ORZrPs0elEABFAgji+XHZNTc3xzXOsENwH0ZfqowYF5bO7rRDfrSdmQFuP38Pn82kW+VNmPg8dOlT0er26wY/S/eXxeNS1kxISEjS/hzJENysrq9nv5HA4BCfKrmwzajKZNENcPR6PZsZ0tGHAvWkGNLuSiE6j8FEvrSkoKFBHw3TkfeV+d999N0aMGNGuTWbmz5+PQYMG4ciRI2hqalK7ilJF8NrRo3CEjDMpB/DvH/0ImzZtwrFjx2BqNGHEFyOOJ6fDupSCbwPlb5x8nXwEWPEC0K8SwBHAbDYjEAjg4osvhl5/vEpJSEjALbfcgi+//BJNTU1obGzEt99+q87GrqmpQXZ2NqxWK7KysjS/x/Dhw/Gvf/0LN9xwg3pc+V3KysrgcDjw5JNPwuv1QqfToaGhATabDWlpaVizZg2uv/56bN68WR1QEOn3VkYvKQlqChOT8BJDbDFQdxKrVTejbdATfqw991PWVQoGg5rjod0tXq9XdDqdBHHyiT/0766wiWbhayWFthRC/4qdkNGjR6tP3cpWnFlZWWK329VVWD0ej+j1eklJSZFp06apXUnTpk1Tn+5D11CKtAFPpN9FWZY7Pj5e4uLipLCwUD0nGAz2+glt4diVRHQaxWrWbKQKP9Ix5X6VlZWSkZEhFRUVUbs8iouLxe12S3FxsYiITJs2TV3dVNkpTdl32Wg0yv+76KLIlXz4jmtWSP+r+0vwkshBIWiAOkKqsrJSXfwuMzNTsrKy1G4lk8kko0ePliuvvFLsdrtkZWVpVnZVupXMZrOao1COJScni9VqlWnTponH45Fx48apM7qV3y58D+3Qf1e9rauoNQwMRD1QW1sMCiVfoCSHlUoutLILrwhDK3ePx6Mu4W2320Wv14vf75f7rNbIlT1OLs2t0+mkLDm5xfNC/5TK3GKxSCAQUHdmU95TgoHNZmu2HwRO5COysrLU/EHodqVKC0PZqyI/P1/97kVFRZpgEaq3JZdbw8BA1MO1pRWidLNceuml6vj9SJWdEizCN9opKytTg8KAAQPUzW7KysqidispLYdo7wdxckvQ0FFPBoNBzGazDB06VKqqqqSoqEjdjnPatGnidDrFarWK0+mUSy+9tFliXK/Xi8vlkoqKCvF4PHLxxReLxWIRo9Gothji4uLUgNeW35UthlPDwEB0moVvdh+tMistLVX3ZlAmaAUCAamoqFD74ZVg4fP5xO/3q5v1KBv8KJVvRUWFVFdXi8PhEL1eH7Xy/7qFoBD6BO9yuaSsrEzdPyHaBjjKBkChLYX09HRNYAEgJpNJSktLpbS0VA0cDodDRI4Hv9BNe5QFApVhrpF+10h682J6DAxEPYgys7eyslI9FlpBhVZm4S2C8OSsEhhCN+4JDyYpKSkCQFJSUqSqqkpTyQYCAbFarWrO4S69PmIQCP+bZTRKUlKSJCQkqN1CycnJ4vF4JDMzU+Lj4yUrK0vcbrcAx1dRVXIjmZmZotPpZPDgwZpAoLQALBaL2uWlfMfs7Gyx2Wzqbxb6u9TU1KgzwW02W8T8TLSKv6dv39kSBgaibqS1yih03H1rnw+v5MMDhfLa6/WqLYbw3c2UiliZAzBt2jTR6/VSUlKiaXEoAePnrQQHpaVQUVGhjgRyuVzNls0ObwEouZHk5GTR6XRy5ZVXyvjx4wWADB06VE2ul5WVSUZGhjqKqaioSJKTk8Vut6utA5/PJ16vV4qKiiQYDEpWVpZ4PB4pKirS7AB3qv+uejIGBqJupLWn0EgthrYKDxShw1CVyk25vzIDWOlSUo4pFbgyE1j5zKBBg8RkMklycrLsiouLGBRqwrp6QruvwhPISv+/8s85OTlisVjUwJiSkqIpi/LdXC6XABC73a6WXWmVZGRkqK0ch8MhWVlZ4vf7NWsd9bUkczQMDETdSLSn0PYkP9v6JBsIBMRkMmm6UMKHt06bNk2MRqOkpaWpT/FGo1Ht2qmoqBC32y0ZGRliMBii5hqUv/vj4sRkMqldPwCkpKREvF6vpKamSmZmplqRh45OUs5NSEhQgxZOJK/LysrEarWKxWIRl8slZrNZJk6cqOZRQhfMU2ZBDx06VDIyMsTr9WqCcLR5D30NAwNRJ4tFl4PypBu6YX000eYzBINBKSoqUrtQlJ3awpOuyjWysrI0axEp/f5VVVVqElipvFsafRT+98SgQc0SxaFzLJQuotAhqqErpWZkZEhJSYnaclKS4larVV0Lqi0jtCZNmiTBYLDDEwJ7MwYGok4Wi8pGedLNzs5u9TrRZkC73W6x2WySkJAgNptNncimCO2mUkYcGY1G0ev1MnToULWLpaqqSsxmszopraWgsMNojHh8jt2uBhwlb2C1WiPOS3A6nTJu3DjNMWWUVFVVVbMtPKP9fkrZlWSz0qUU2jrozXmD9mBgIOpkp7pCamvnteXzoS0Gn88nLper2dIXSv+9zWaTjIwMycjIkP79+4vb7dbMmFae0I1Go/Tr1y9qUNg2c6YEAgF5pF+/iO+vuPJKsVgskpubKy6XS4YOHaqZyawMrVWS4qEJ6fBcRyAQkLi4OHG5XGq3USgld6B8L5PJpLZCQvMldBwDA1EXa6lF0ZZKv70tktCNZkKvrbQYlKRvcnKy2iIIrTyrqqrE4XDI0KFDpcRgiFjpl5yYVazc69mcnKijlJShssoktISEhIjJdaV8yhBag8GgaTG4XC61ayu8sldaCkrXk06nU1spreUT+trkNhEGBqIu11LlH75+UXuXwgi/T+iOY9ECitPpVFsEFotFrTzDdzr7Y5S1ku468SSuPJUrI4qitSyePjHDWWkRRHqCV5asUFoVoctZKBW9sjx2eIshtBsuPz9fTZYrrZHWhM976AtdTQwMRN1Y+IqnbWkdRKu8Qmc+B4PBiNtSVlRUqF0syoxiZdtM5d5FRUXyZGZmxEr+yczMZjkCpetHr9fLn/z+iJ/b+4tfqNt3KvcLL7syEklJNLe2gqoiNHGvDGc1m83NRiRFE9pi6CvJaQYGom6spf2co4lWeYW2GIqLi9VKPrSbRHlq1+v1mhnRyueLiopkiM0muyNU7psLCqSyslIsFovodDq54IIL1JZH3IlhqtGWzzhktco1I0dKfHy8+P1+zXdURhEpXUg6nU4GDBggLpdL7a6K9HspFXplZaX6HbOystTv5/P52r1fM1sMHcPAQNTFouUQFKG7kinrITkcDqmurtY8tYc+iSv9+9nZ2cdXPbXZ5FhioibRLHIyCex2u8Xj8airmrrdbnWUk8FgkDl2u/rZ3YCMODGbWtliNCsrS80JKCOVnE6nusCf0qoxm81Rlx6PNFmtqKhIXW3V5/Op5/bFPEJLGBiIepDWWhDh6yYpFWzo0toOh0OzNpLdbheDwSB+vz/qfUNHLCkV6M6//EUOWa2y9xe/UM8rKytTF7kzGAyaUUKhW2w6HA4pPhEUfCHdTdXV1er6RXFxcWpZAUj//v3VlsDgwYNFr9fL6NGjI7YYKisrxe12S05OjlRXV2uGq/p8PklISJD8/Hz1t+KMZy0GBqIeJLybqKXXoRWsUuGVlpZKZmamZvx+YWGh2Gw2KSwsjHjPmpoamTRpkrqTWmhZLs/J0QQi+4n5CTqdLuK8AqUVoqxpNGf6dM2Map/Pp+n+qaysFJ/PJ263W12vqbq6WpOoDl0dVWkleTweTTdYaMUfaac6thi0GBiIepC2tBiKi4vVGb0tbcYT7ZrhWspZhHbDeL1emTRpktpVo1TaoX3+ysY9Pp9Pc63QfSDCV4QNXRSvpqZG7Uoym83q7G29Xq9+zuv1itfr1YxOCv0dOpK36WsYGIh6mdZGzrT36bgtgcPr9arDOZUuo8zMTM16RKGb/4QHBiUAhG89qgxFDR3CqizHkZWVJYFAQCZNmiQJCQni8/k0QSpaIIs0WqsvjDRqDwYGol6mtYrc7/eL0WhsMacQLXhEGhoa6Vyl316pvEOf9JWlKOLj49WuqdAAEHq9lu6ntC7CRzEp319ZXC98uG1bA0ZfxsBA1EPEqgIrKioSm80mRUVFUc+JloyN9AQfXuGGd/8oCd+srCw1X6C0HJR9JUIDgN/vV0cNhXZD+f1+KSwsVL9/eEAK7y5SZkFbrdaIvx8DQnQMDEQ9REe7PCJViIWFhVE3vhdpX4sh/PrhwSMQCKgT3JSRQampqWK1WjV7QSiUwOXz+SQ7O1syMjLE5/NF3GUtVGgwUxLQyhLhkbALKToGBqIeoqNPvOEjlYqKisRut4vVau2U4ZnhwaO6uloyMzPFbDZrKulIFbPSMsjPz5fCwkLxeDxitVolKyur2aZCyvnK/s2h+Ym27KvQ3hZDX2ph9IjA8PXXX0teXp6kpKRIXFycnHvuubJu3bo2fZaBgXqrtj7xhs9tiI+PF51OJ3a7PWbDMyNVmmVlZWIwGGT06NHNho+KRG59hK6CGjozO3y0Uuj5yjwHt9utBobOmJfQl1oY3T4w7N69WwYNGiQ333yz/Oc//5EtW7bI3/72N/nss8/a9HkGBuqtOvIEqySBwyd4napIeYvQTXzy8/ObDTtVtt20WCzNhpWWlZWJw+GQiooKNU+hdEOFD8dVWgxer1cNBtG6wk7lqZ8tho6LeWCYNWuWXHLJJR3+PAMDkZZSwQWDwZg9AYeOdFKuX1JSorYYlMo09Km7sLBQfdr3eDya6yl5CofDIW63W93Cs7i4WF13KTx53pZhuH3pqf9UdPvAkJ2dLTNnzpTrrrtO+vXrJ+edd54888wzbf48AwNRZLF8Ag590rdareoqrsrcgvDzlJFDXq9XTCaTuoS2co7SCpg2bZra6lAmxilbfbY03PZ0fOfOuF530e0Dg8ViEYvFIvfcc4+sX79e5s+fL1arVf7whz9EPP/IkSOyb98+9W/btm0MDETt0FLitqWKsKamRrM6q06nU7cFtVgsMm7cOHVVV+UaymqvygxpZXKc2WwWq9WqWWo7ISFBXC6X+Hw+8Xq9zYJOa+XrDL21BRLrwKBHjDU1NWH48OF4+OGHcf7556OoqAgFBQV4+umnI54/d+5cJCUlqX8DBgyIdZGIeqza2lrMmTMHtbW1Uc+ZOnUqtm7diqlTpzZ7b8GCBVi8eDEuvvhiPPnkkxg8eDCWL1+uvte/f39YrVY8+OCDyMzMhNFoxJYtW1BfX4+qqiqICAKBAOrq6lBeXo7Vq1cjMTERy5Ytw7Jly1BcXAyr1YqmpibU19fDZDLBaDSipKQERUVF+MlPfoJXX30VTqcTmzdvRmFhIUaPHo3169erZVi2bBkWLFjQOT9gmIKCAkyYMAEFBQWn5X49VkzCS4iBAwfKf//3f2uOPfXUU5Kenh7xfLYYiKJryxNuay0GZWE+q9XabL+GSPMZUlJS1BaDMhzV5XJJUVGRZkc4ZVXU0EX2lGt4vV5NOaqrq8Xn84nNZju+DPiJ0Ufha0VRx3T7rqQbb7yxWfJ55syZbe5fZI6B6KT2dLVEO1fJE1RUVLQ4VyBSgCktLRWXyyU2m02zuqkifJipMkva7XZHPFfpclLuoSzF0dbd2Siybh8Y3nnnHTEajfLQQw/Jp59+KpWVlRIfHy8vvPBCmz7PwEB9Taz62VtqXYTub9DeSWItPdGHjyxqqfUSur9D6LLjoQv6Ucd0+8AgIrJs2TIZNmyYWCwWycrK4qgkoha0NyEaLZC0FGACgYBYLBYxGo2SmZnZpnud6p4HkcrTlmXFqf16RGA4FQwM1Ne0t3LsyMia0M1yHA5H1IXtQrVnNnLoXg+hs7aVfZpDRzW53W7Jy8trNeicjqDRWwJTrOtNY9ekvIlIkZ6ejtmzZ7f5fGVETXtG1gwfPhyvvPIKpk6divnz5yM9PR0AUFxcjHfeeQfFxcVYuXIlamtrsWDBAhQUFKC8vBzFxcUoLy9v9frK6KJVq1Zh//79qKurw8GDB2E2m3HgwAGsWrUKAKDT6QAAq1evxu7du9X7tnRNAO36fdrjdNyjJ2JgIOpBQitupXJvq7Fjx2LLli2aY+GVf3hFGa3SDqcEqQkTJmDZsmXYv38/Xn31VQCA3+/XDBFNTEzERRddhEceeQR333035syZE/H7RAqAbfn+7fmNOhJk+4SYtDtiiF1JRNG1pRvpVEYyhb5uy4qnLV23LcNQ29st1pbze+sktpawK4moD2vLE257ukfCzw3t1rr44ovViXPhLY3W7NixA9XV1SgvL2/xyX7//v0YNWpUm5/Y2/L92QqIgZiElxhii4Ho1MRi7oNIy0NPW9OW7UhbG6paU1Mj+fn54vF4OlSGviTW9aZORKSrg1Oouro6JCUlYd++fbDb7V1dHKJe41TyE+29x6effoo///nPmDx5Mn77299GPTcvLw+7du3Ctdde26x1M2fOHDz88MM4evQoMjIy2t1q6UtiXW/GfK0kIuqeOrIuUW1tLYqLi3HXXXe1uF5T+D3S0tJw1113tdiVlZ6ejsrKSlx77bURu30KCgrw4x//GB6PB/Pnz29zmenUMcdA1EdE63tXnvKVEUWhLYoFCxbgxRdfBHB8NFFrOYvQe7SlVdLSUN309HQ8//zzrV6DYo9dSUR93Jw5c7Bs2TIkJiZi//79mDBhglpZ19bWory8HDqdDsFgsNO6oOjUxLreZIuBqI8Ln4MQ3qKw2+0dzkucjrwGxR5bDEQU1Zw5c7B06VI4nU5UVla2u3JXWiOhrRCKPSafieiUtWUDIOB4a8LpdGLXrl0d2kyHG+P0TAwMRH1QW0cotTZyiHon5hiI+qD2zA5u7yJ/obhIXc/EwEDUB51KZd8eXJ6iZ2JXElEv1NYcQmdTAhBHJPUsDAxEvVBHZjm3Zv369Rg9ejTWr18fs2tS98SuJKJeqDO6cMI39aHeiy0Gol7oVLtwInVFlZeXY8SIEW3a0Y16NrYYiKiZSKOJhg8fzpZCH8HAQETNtKUristd9F7sSiKiZtrSFdUZCW7qHthiIKIOaa1VwRZFz8UWAxF1SGutCrYoei62GIioU3DWc8/FZbeJiHo4LrtNRESdioGBiIg0GBiIiEiDgYGIiDQYGIiISIOBgYiINBgYiIhIg4GBiIg0GBiIiEiDgYGIiDQYGIiISIOBgYiINBgYiIhIg4GBiIg0GBiIiEiDgYGIiDQYGIiISIOBgYiINBgYiIhIg4GBiIg0GBiIiEiDgYGIiDQYGIiISIOBgYiINBgYiIhIg4GBiIg0GBiIiEiDgYGIiDQ6PTDMnTsXOp0OM2fO7OxbERFRDHRqYFi7di2eeeYZ5OTkdOZtiIgohjotMBw4cAB5eXlYsGABkpOTO+s2REQUY50WGKZPn47x48fj8ssv76xbEBFRJzB2xkUXL16M6upqrFu3rtVz6+vrUV9fr76uq6vrjCIREVEbxbzFsG3bNtxxxx2orKyE1Wpt9fy5c+ciKSlJ/RswYECsi0RERO2gExGJ5QVfeeUVTJw4EQaDQT3W2NgInU4HvV6P+vp6zXuRWgwDBgzAvn37YLfbY1k0IqJeqa6uDklJSTGrN2PelXTZZZdh06ZNmmO33HILsrKyMGvWLE1QAACLxQKLxRLrYhARUQfFPDAkJiZi2LBhmmMJCQlwOp3NjhMRUffDmc9ERKTRKaOSwq1atep03IaIiGKALQYiItJgYCAiIg0GBiIi0mBgICIiDQYGIiLSYGAgIiINBgYiItJgYCAiIg0GBiIi0mBgICIiDQYGIiLSYGAgIiINBgYiItJgYCAiIg0GBiIi0mBgICIiDQYGIiLSYGAgIiINBgYiItJgYCAiIg0GBiIi0mBgICIiDQYGIiLSYGAgIiINBgYiItJgYCAiIg0GBiIi0mBgICIiDQYGIiLSYGAgIiINBgYiItJgYCAiIg0GBiIi0mBgICIiDQYGIiLSYGAgIiINBgYiItJgYCAiIg0GBiIi0mBgICIiDQYGIiLSYGAgIiINBgYiItJgYCAiIg0GBiIi0mBgICIiDQYGIiLSYGAgIiINBgYiItJgYCAiIg0GBiIi0mBgICIiDQYGIiLSYGAgIiINBgYiItKIeWCYO3cuLrjgAiQmJsLlcuGaa67Bxx9/HOvbEBFRJ4l5YFi9ejWmT5+ONWvWYMWKFWhoaMCYMWNw8ODBWN+KiIg6gU5EpDNv8O2338LlcmH16tX4wQ9+0Or5dXV1SEpKwr59+2C32zuzaEREvUKs601jDMrUon379gEAUlJSIr5fX1+P+vp69XVdXV1nF4mIiFrQqclnEUEwGMQll1yCYcOGRTxn7ty5SEpKUv8GDBjQmUUiIqJWdGpX0vTp0/H666/jrbfegsfjiXhOpBbDgAED2JVERNRGPaYr6bbbbsOrr76KN998M2pQAACLxQKLxdJZxSAionaKeWAQEdx22214+eWXsWrVKgwePDjWtyAiok4U88Awffp0LFq0CH/+85+RmJiIHTt2AACSkpIQFxcX69sREVGMxTzHoNPpIh5fuHAhbr755lY/z+GqRETt0+1zDJ08LYKIiDoZ10oiIiINBgYiItJgYCAiIg0GBiIi0mBgICIiDQYGIiLSYGAgIiINBgYiItJgYCAiIg0GBiIi0mBgICIiDQYGIiLSYGAgIiINBgYiItJgYCAiIg0GBiIi0mBgICIiDQYGIiLSYGAgIiINBgYiItJgYCAiIg0GBiIi0mBgICIiDQYGIiLSYGAgIiINBgYiItJgYCAiIg0GBiIi0mBgICIiDQYGIiLSYGAgIiINBgYiItJgYCAiIg0GBiIi0mBgICIiDQYGIiLSYGAgIiINBgYiItJgYCAiIg0GBiIi0mBgICIiDQYGIiLSYGAgIiINBgYiItJgYCAiIg0GBiIi0mBgICIiDQYGIiLSYGAgIiINBgYiItJgYCAiIg0GBiIi0mBgICIiDQYGIiLS6LTA8NRTT2Hw4MGwWq3w+Xz45z//2Vm3IiKiGOqUwLBkyRLMnDkT9913HzZs2IDvf//7GDduHL766qvOuB0REcWQTkQk1he98MILMXz4cDz99NPqsezsbFxzzTWYO3dui5+tq6tDUlIS9u3bB7vdHuuiERH1OrGuN40xKJPG0aNHUV1djbvvvltzfMyYMXj77bebnV9fX4/6+nr19b59+wAc/6JERNQ6pb6M1XN+zAPDd999h8bGRvTv319zvH///tixY0ez8+fOnYs5c+Y0Oz5gwIBYF42IqFfbtWsXkpKSTvk6MQ8MCp1Op3ktIs2OAcA999yDYDCovt67dy8GDRqEr776KiZfsCerq6vDgAEDsG3btj7frcbf4iT+Fsfxdzhp3759GDhwIFJSUmJyvZgHhtTUVBgMhmatg507dzZrRQCAxWKBxWJpdjwpKanP/8tW2O12/hYn8Lc4ib/FcfwdTtLrYzOeKOajksxmM3w+H1asWKE5vmLFCowcOTLWtyMiohjrlK6kYDCIn/zkJ8jNzYXf78czzzyDr776ClOnTu2M2xERUQx1SmC4/vrrsWvXLjz44IPYvn07hg0bhr/85S8YNGhQq5+1WCyYPXt2xO6lvoa/xUn8LU7ib3Ecf4eTYv1bdMo8BiIi6rm4VhIREWkwMBARkQYDAxERaTAwEBGRRrcLDFyu+/gyIRdccAESExPhcrlwzTXX4OOPP+7qYnW5uXPnQqfTYebMmV1dlC5RU1OD/Px8OJ1OxMfH47zzzkN1dXVXF+u0a2howP3334/BgwcjLi4OQ4YMwYMPPoimpqauLlqne/PNNzFhwgSkp6dDp9PhlVde0bwvIigtLUV6ejri4uIQCATwwQcftPs+3SowcLnu41avXo3p06djzZo1WLFiBRoaGjBmzBgcPHiwq4vWZdauXYtnnnkGOTk5XV2ULrFnzx5cfPHFMJlM+Otf/4oPP/wQ5eXlcDgcXV200+7RRx/F/PnzUVFRgc2bN+Oxxx7D//zP/+A3v/lNVxet0x08eBDnnnsuKioqIr7/2GOP4fHHH0dFRQXWrl2LtLQ0XHHFFdi/f3/7biTdyIgRI2Tq1KmaY1lZWXL33Xd3UYm6h507dwoAWb16dVcXpUvs379fzjzzTFmxYoWMGjVK7rjjjq4u0mk3a9YsueSSS7q6GN3C+PHj5ac//anm2KRJkyQ/P7+LStQ1AMjLL7+svm5qapK0tDR55JFH1GNHjhyRpKQkmT9/fruu3W1aDMpy3WPGjNEcj7Zcd1+iLEUeqwWyeprp06dj/PjxuPzyy7u6KF3m1VdfRW5uLn70ox/B5XLh/PPPx4IFC7q6WF3ikksuwd///nd88sknAIB3330Xb731Fq666qouLlnX2rJlC3bs2KGpQy0WC0aNGtXuOrTTVldtr/Yu191XiAiCwSAuueQSDBs2rKuLc9otXrwY1dXVWLduXVcXpUt98cUXePrppxEMBnHvvffinXfewe233w6LxYKbbrqpq4t3Ws2aNQv79u1DVlYWDAYDGhsb8dBDD+HGG2/s6qJ1KaWejFSHbt26tV3X6jaBQdHW5br7ihkzZuC9997DW2+91dVFOe22bduGO+64A2+88QasVmtXF6dLNTU1ITc3Fw8//DAA4Pzzz8cHH3yAp59+us8FhiVLluCFF17AokWLcM4552Djxo2YOXMm0tPTMWXKlK4uXpeLRR3abQJDe5fr7gtuu+02vPrqq3jzzTfh8Xi6ujinXXV1NXbu3Amfz6cea2xsxJtvvomKigrU19fDYDB0YQlPH7fbjbPPPltzLDs7G0uXLu2iEnWdn//857j77rtxww03AAC8Xi+2bt2KuXPn9unAkJaWBuB4y8HtdqvHO1KHdpscA5frPklEMGPGDLz00kv4xz/+gcGDB3d1kbrEZZddhk2bNmHjxo3qX25uLvLy8rBx48Y+ExQA4OKLL242ZPmTTz5p08KUvc2hQ4ea7TtgMBj6xHDVlgwePBhpaWmaOvTo0aNYvXp1++vQWGTHY2Xx4sViMpnk//7v/+TDDz+UmTNnSkJCgnz55ZddXbTTatq0aZKUlCSrVq2S7du3q3+HDh3q6qJ1ub46Kumdd94Ro9EoDz30kHz66adSWVkp8fHx8sILL3R10U67KVOmyBlnnCGvvfaabNmyRV566SVJTU2VkpKSri5ap9u/f79s2LBBNmzYIADk8ccflw0bNsjWrVtFROSRRx6RpKQkeemll2TTpk1y4403itvtlrq6unbdp1sFBhGRJ598UgYNGiRms1mGDx/eJ4doAoj4t3Dhwq4uWpfrq4FBRGTZsmUybNgwsVgskpWVJc8880xXF6lL1NXVyR133CEDBw4Uq9UqQ4YMkfvuu0/q6+u7umidbuXKlRHrhilTpojI8SGrs2fPlrS0NLFYLPKDH/xANm3a1O77cNltIiLS6DY5BiIi6h4YGIiISIOBgYiINBgYiIhIg4GBiIg0GBiIiEiDgYGIiDQYGIiISIOBgYiINBgYiIhIg4GBiIg0GBiIiEjj/wP/v5JP1LVZLgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from kf_book.pf_internal import plot_pf\n",
    "from kf_book.gif_animate import animate\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",
    "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}