From 15983a196c5cbb272ef9ed64815638670bada9d2 Mon Sep 17 00:00:00 2001 From: Adel Date: Mon, 17 Dec 2018 22:29:56 +0100 Subject: [PATCH] finish updated and resampling --- SLAM/FastSLAM1/FastSLAM1.ipynb | 401 +++++++++++++++++++++++++++++++-- 1 file changed, 376 insertions(+), 25 deletions(-) diff --git a/SLAM/FastSLAM1/FastSLAM1.ipynb b/SLAM/FastSLAM1/FastSLAM1.ipynb index ecdc3024..705bba8b 100644 --- a/SLAM/FastSLAM1/FastSLAM1.ipynb +++ b/SLAM/FastSLAM1/FastSLAM1.ipynb @@ -80,7 +80,7 @@ "metadata": {}, "source": [ "### 1- Predict\n", - "The following equations and code snipets we can see how the particles distribution evolves in case we provide only the control $(v,w)$, which are the linear and angular velocity repsectively. \n", + "The following equations and code snippets we can see how the particles distribution evolves in case we provide only the control $(v,w)$, which are the linear and angular velocity repsectively. \n", "\n", "$\\begin{equation*}\n", "F=\n", @@ -133,12 +133,24 @@ "\n" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The following snippets playsback the recorded trajectory of each particle.\n", + "\n", + "To get the insight of the motion model change the value of $R$ and re-run the cells again. As R is the parameters that indicates how much we trust that the robot executed the motion commands.\n", + "\n", + "It is interesting to notice also that only motion will increase the uncertainty in the system as the particles start to spread out more. If observations are included the uncertainty will decrease and particles will converge to the correct estimate." + ] + }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ + "# CODE SNIPPET #\n", "import numpy as np\n", "import math\n", "from copy import deepcopy\n", @@ -153,16 +165,24 @@ "\n", "DT = 0.1 # time tick [s]\n", "SIM_TIME = 50.0 # simulation time [s]\n", + "MAX_RANGE = 20.0 # maximum observation range\n", + "M_DIST_TH = 2.0 # Threshold of Mahalanobis distance for data association.\n", "STATE_SIZE = 3 # State size [x,y,yaw]\n", + "LM_SIZE = 2 # LM srate size [x,y]\n", "N_PARTICLE = 100 # number of particle\n", + "NTH = N_PARTICLE / 1.5 # Number of particle for re-sampling\n", "\n", "class Particle:\n", "\n", " def __init__(self, N_LM):\n", " self.w = 1.0 / N_PARTICLE\n", - " self.x = np.random.uniform(-0.1,0.1,1)[0]\n", - " self.y = np.random.uniform(-0.1,0.1,1)[0]\n", + " self.x = 0.0\n", + " self.y = 0.0\n", " self.yaw = 0.0\n", + " # landmark x-y positions\n", + " self.lm = np.zeros((N_LM, LM_SIZE))\n", + " # landmark position covariance\n", + " self.lmP = np.zeros((N_LM * LM_SIZE, LM_SIZE))\n", "\n", "def motion_model(x, u):\n", " F = np.array([[1.0, 0, 0],\n", @@ -194,6 +214,7 @@ "def pi_2_pi(angle):\n", " return (angle + math.pi) % (2 * math.pi) - math.pi\n", "\n", + "# END OF SNIPPET\n", "\n", "N_LM = 0 \n", "particles = [Particle(N_LM) for i in range(N_PARTICLE)]\n", @@ -208,17 +229,6 @@ " history.append(deepcopy(particles))\n" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The following snipets playsback the recorded trajectory of each particle.\n", - "\n", - "To get the insight of the motion model change the value of $R$ and re-run the cells again. As R is the parameters that indicates how much we trust that the robot executed the motion commands.\n", - "\n", - "It is interesting to notice also that only motion will increase the uncertainty in the system as the particles start to spread out more. If observations are included the uncertainty will decrease and particles will converge to the correct estimate." - ] - }, { "cell_type": "code", "execution_count": 3, @@ -227,7 +237,7 @@ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "b5fd2bf875bd43f6bf6dd92d638d2973", + "model_id": "3f279cf1dcaf4ec886c01942c36e601d", "version_major": 2, "version_minor": 0 }, @@ -257,7 +267,7 @@ " plt.plot(x, y, '.r')\n", " plt.axis([-20,20, -5,25])\n", "\n", - "interact(plot_particles, t=(0,499,1));" + "interact(plot_particles, t=(0,len(history)-1,1));" ] }, { @@ -266,7 +276,353 @@ "source": [ "### 2- Update\n", "\n", - "work in progress" + "For the update step it is useful to observe a single particle and the effect of getting a new measurements on the weight of the particle.\n", + "\n", + "As mentioned earlier, each particle maintains $N$ $2x2$ EKFs to estimate the landmarks, which includes the EKF process described in the EKF notebook. The difference is the change in the weight of the particle according to how likely the measurement is.\n", + "\n", + "The weight is updated according to the following equation: \n", + "\n", + "$\\begin{equation*}\n", + "w_i = |2\\pi Q|^{\\frac{-1}{2}} exp\\{\\frac{-1}{2}(z_t - \\hat z_i)^T Q^{-1}(z_t-\\hat z_i)\\}\n", + "\\end{equation*}$\n", + "\n", + "Where, $w_i$ is the computed weight, $Q$ is the measurement covariance, $z_t$ is the actual measurment and $\\hat z_i$ is the predicted measurement of particle $i$.\n", + "\n", + "To experiment this, a single particle is initialized then passed an initial measurement, which results in a relatively average weight. However, setting the particle coordinate to a wrong value to simulate wrong estimation will result in a very low weight. The lower the weight the less likely that this particle will be drawn during resampling and probably will die out." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "initial weight 1.0\n", + "weight after landmark intialization 1.0\n", + "weight after update 0.023098460073039763\n", + "weight after wrong prediction 7.951154575772496e-07\n" + ] + } + ], + "source": [ + "# CODE SNIPPET #\n", + "def observation(xTrue, xd, u, RFID):\n", + "\n", + " # calc true state\n", + " xTrue = motion_model(xTrue, u)\n", + "\n", + " # add noise to range observation\n", + " z = np.zeros((3, 0))\n", + " for i in range(len(RFID[:, 0])):\n", + "\n", + " dx = RFID[i, 0] - xTrue[0, 0]\n", + " dy = RFID[i, 1] - xTrue[1, 0]\n", + " d = math.sqrt(dx**2 + dy**2)\n", + " angle = pi_2_pi(math.atan2(dy, dx) - xTrue[2, 0])\n", + " if d <= MAX_RANGE:\n", + " dn = d + np.random.randn() * Qsim[0, 0] # add noise\n", + " anglen = angle + np.random.randn() * Qsim[1, 1] # add noise\n", + " zi = np.array([dn, pi_2_pi(anglen), i]).reshape(3, 1)\n", + " z = np.hstack((z, zi))\n", + "\n", + " # add noise to input\n", + " ud1 = u[0, 0] + np.random.randn() * Rsim[0, 0]\n", + " ud2 = u[1, 0] + np.random.randn() * Rsim[1, 1] + OFFSET_YAWRATE_NOISE\n", + " ud = np.array([ud1, ud2]).reshape(2, 1)\n", + "\n", + " xd = motion_model(xd, ud)\n", + "\n", + " return xTrue, z, xd, ud\n", + "\n", + "def update_with_observation(particles, z):\n", + " for iz in range(len(z[0, :])):\n", + "\n", + " lmid = int(z[2, iz])\n", + "\n", + " for ip in range(N_PARTICLE):\n", + " # new landmark\n", + " if abs(particles[ip].lm[lmid, 0]) <= 0.01:\n", + " particles[ip] = add_new_lm(particles[ip], z[:, iz], Q)\n", + " # known landmark\n", + " else:\n", + " w = compute_weight(particles[ip], z[:, iz], Q)\n", + " particles[ip].w *= w\n", + " particles[ip] = update_landmark(particles[ip], z[:, iz], Q)\n", + "\n", + " return particles\n", + "\n", + "def compute_weight(particle, z, Q):\n", + " lm_id = int(z[2])\n", + " xf = np.array(particle.lm[lm_id, :]).reshape(2, 1)\n", + " Pf = np.array(particle.lmP[2 * lm_id:2 * lm_id + 2])\n", + " zp, Hv, Hf, Sf = compute_jacobians(particle, xf, Pf, Q)\n", + " dx = z[0:2].reshape(2, 1) - zp\n", + " dx[1, 0] = pi_2_pi(dx[1, 0])\n", + "\n", + " try:\n", + " invS = np.linalg.inv(Sf)\n", + " except np.linalg.linalg.LinAlgError:\n", + " print(\"singuler\")\n", + " return 1.0\n", + "\n", + " num = math.exp(-0.5 * dx.T @ invS @ dx)\n", + " den = 2.0 * math.pi * math.sqrt(np.linalg.det(Sf))\n", + " w = num / den\n", + "\n", + " return w\n", + "\n", + "def compute_jacobians(particle, xf, Pf, Q):\n", + " dx = xf[0, 0] - particle.x\n", + " dy = xf[1, 0] - particle.y\n", + " d2 = dx**2 + dy**2\n", + " d = math.sqrt(d2)\n", + "\n", + " zp = np.array(\n", + " [d, pi_2_pi(math.atan2(dy, dx) - particle.yaw)]).reshape(2, 1)\n", + "\n", + " Hv = np.array([[-dx / d, -dy / d, 0.0],\n", + " [dy / d2, -dx / d2, -1.0]])\n", + "\n", + " Hf = np.array([[dx / d, dy / d],\n", + " [-dy / d2, dx / d2]])\n", + "\n", + " Sf = Hf @ Pf @ Hf.T + Q\n", + "\n", + " return zp, Hv, Hf, Sf\n", + "\n", + "def add_new_lm(particle, z, Q):\n", + "\n", + " r = z[0]\n", + " b = z[1]\n", + " lm_id = int(z[2])\n", + "\n", + " s = math.sin(pi_2_pi(particle.yaw + b))\n", + " c = math.cos(pi_2_pi(particle.yaw + b))\n", + "\n", + " particle.lm[lm_id, 0] = particle.x + r * c\n", + " particle.lm[lm_id, 1] = particle.y + r * s\n", + "\n", + " # covariance\n", + " Gz = np.array([[c, -r * s],\n", + " [s, r * c]])\n", + "\n", + " particle.lmP[2 * lm_id:2 * lm_id + 2] = Gz @ Q @ Gz.T\n", + "\n", + " return particle\n", + "\n", + "def update_KF_with_cholesky(xf, Pf, v, Q, Hf):\n", + " PHt = Pf @ Hf.T\n", + " S = Hf @ PHt + Q\n", + "\n", + " S = (S + S.T) * 0.5\n", + " SChol = np.linalg.cholesky(S).T\n", + " SCholInv = np.linalg.inv(SChol)\n", + " W1 = PHt @ SCholInv\n", + " W = W1 @ SCholInv.T\n", + "\n", + " x = xf + W @ v\n", + " P = Pf - W1 @ W1.T\n", + "\n", + " return x, P\n", + "\n", + "def update_landmark(particle, z, Q):\n", + "\n", + " lm_id = int(z[2])\n", + " xf = np.array(particle.lm[lm_id, :]).reshape(2, 1)\n", + " Pf = np.array(particle.lmP[2 * lm_id:2 * lm_id + 2, :])\n", + "\n", + " zp, Hv, Hf, Sf = compute_jacobians(particle, xf, Pf, Q)\n", + "\n", + " dz = z[0:2].reshape(2, 1) - zp\n", + " dz[1, 0] = pi_2_pi(dz[1, 0])\n", + "\n", + " xf, Pf = update_KF_with_cholesky(xf, Pf, dz, Q, Hf)\n", + "\n", + " particle.lm[lm_id, :] = xf.T\n", + " particle.lmP[2 * lm_id:2 * lm_id + 2, :] = Pf\n", + "\n", + " return particle\n", + "\n", + "# END OF CODE SNIPPET #\n", + "\n", + "\n", + "\n", + "# Setting up the landmarks\n", + "RFID = np.array([[10.0, -2.0],\n", + " [15.0, 10.0]])\n", + "N_LM = RFID.shape[0]\n", + "\n", + "# Initialize 1 particle\n", + "N_PARTICLE = 1\n", + "particles = [Particle(N_LM) for i in range(N_PARTICLE)]\n", + "\n", + "xTrue = np.zeros((STATE_SIZE, 1))\n", + "xDR = np.zeros((STATE_SIZE, 1))\n", + "\n", + "print(\"initial weight\", particles[0].w)\n", + "\n", + "xTrue, z, _, ud = observation(xTrue, xDR, u, RFID)\n", + "# Initialize landmarks\n", + "particles = update_with_observation(particles, z)\n", + "print(\"weight after landmark intialization\", particles[0].w)\n", + "particles = update_with_observation(particles, z)\n", + "print(\"weight after update \", particles[0].w)\n", + "\n", + "particles[0].x = -10\n", + "particles = update_with_observation(particles, z)\n", + "print(\"weight after wrong prediction\", particles[0].w)\n", + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 3- Resampling\n", + "\n", + "In the reseampling steps a new set of particles are chosen from the old set. This is done according to the weight of each particle. \n", + "\n", + "The figure shows 100 particles distributed uniformly between [-0.5, 0.5] with the weights of each particle distributed according to a Gaussian funciton. \n", + "\n", + "The resampling picks $i \\in 1,...,N$ particles with probability to pick particle with index $i \\propto \\omega_i $, where $\\omega_i$ is the weight of that particle\n", + "\n", + "\n", + "To get the intuition of the resampling step we will look at a set of particles which are initialized with a given x location and weight. After the resampling the particles are more concetrated in the location where they had the highest weights. This is also indicated by the indices \n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbIAAAEkCAYAAABKTLRCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XecXVW5//HPN5OeAAlJDCSEhNADUod2UUFAmkDQixgUBYWL2L1iAfFqLCiWa7voDxFpogiCKKiIIISiFCcQICGEFEoSEtJICKTPPL8/9hrZDFNOMnPKznzfr9d5zT67rPWcNfuc5+y119lbEYGZmVlR9ah2AGZmZp3hRGZmZoXmRGZmZoXmRGZmZoXmRGZmZoXmRGZmZoXmRGY1TdKXJF1ewnpXSfpmF9e9yWVK6ifpVkkrJP2uK+PaXEg6XNK83PNpkg6vYkhWUD2rHYAVm6RngeFAI/AqcBvwiYh4ZRPKOhy4NiK2a54XEd/qmkgr7hSydhkSERuqHUwRRMQe1Y7BislHZNYVToyIgcB+QD3w5Y0tQNLm9qVqNPD0piSxzrbFZtiWZu1yIrMuExHzyY7I9gSQ9CFJ0yWtlDRH0kea123uVpL0RUkLgevStiMkvZIeIyRNlHRtbru3SPqnpOWS5ko6s7VYJJ0gaUpa75+S9sot+6Kk+SmuGZKObOdlDZV0R1r3Hkmjc+XslpYtS+WcmuZ/DfgK8N70Os6S1EPSlyU9J2mRpGskbZXWHyMp0nrPA3el+QfnXutj7XW7SXo2va7HgVcl9Uztd5OkxZKekfSp3PoHSmqQ9LKkFyX9ILfsd5IWpm7ReyXtkVt2laSfSbotvbZ/SNpG0o8kvSTpKUn7tojrAklPpuVXSurbzms4Kk1PlHRDaqeVqduxPrfufpIeTct+J+n6ru5atgKJCD/82OQH8CxwVJoeBUwDvpGevxPYERBwGLAK2C8tOxzYAHwH6AP0S/PmtSh/Ill3I2RHOSuB04BewBBgn7TsKuCbaXpfYBFwEFAHnJHi7APsCswFRqR1xwA7tvHarkr1vS1t+2Pg/rRsQCrnQ2Rd9PsCS4BxLeNOzz8MzALGAgOB3wO/ysUQwDWp3H7ASGApcDzZF853pOfD2vk/TEn/g35pm8lkCbV3qncOcExa/wHgA2l6IHBwi1i3SK/5R8CUFm2yBNgf6EuWdJ8BPpja+pvA3S3impri2hr4R+7/9Lr/N6/flyYCa9LrrwO+DTyYlvUGngM+TbYfvBtY11yuH93v4SMy6wp/kLQcuB+4B/gWQET8OSJmR+Ye4G/AW3PbNQFfjYi1EbG6hHreB9wZEddFxPqIWBoRU1pZ7xzg5xHxUEQ0RsTVwFrgYLJzeX2AcZJ6RcSzETG7nTr/HBH3RsRa4ELgEEmjgBOAZyPiyojYEBGPAjcB72mjnPcDP4iIOZGdP7wAmNCiG3BiRLya2uJ04C8R8ZeIaIqIO4AGsg/2tvwkIuam7Q8gS3pfj4h1ETEH+AUwIa27HthJ0tCIeCUiHmwuJCKuiIiV6TVPBPZuPnpMbo6IyRGxBrgZWBMR10REI3A9WVLPuyTFtQy4iOyLSCnuT6+/EfgVsHeafzDZl4efpP3g98DDJZZpmyEnMusKJ0fEoIgYHREfa05Kko6T9GDqeltO9iE8NLfd4vRhWKpRQHtJp9lo4LzUJbc81T2K7ChsFvAZsg/oRZJ+K2lEO2XNbZ5ICWgZMCLVcVCLOt4PbNNGOSPIjiKaPUf2YTy8tbpS+e9pUf5bgG1LiTVtP6LF9l/K1XcWsAvwlKR/SToBQFKdpIslzZb0MtlRErz+//Zibnp1K88HthPXc2RtUYqFuelVQN+U+EcA8yMif8XzuVi35ZPCVhaS+pAdoXwQ+GNErJf0B7JuxmYtb73Q0a0Y5gIHllD9XOCiiLiotYUR8RvgN5K2BH5O1r35gTbKGtU8IWkgWffYC6mOeyLiHSXEQ9pmdO759mRdqy8CzaM0W34w/yoi/qvE8lvb/pmI2LnVFSNmAqdJ6kHWNXejpCFpejxwFFkS2wp4idf/3zbWqNz09mRt0RkLgJGSlEtmpX7Jsc2Qj8isXHqTdeEtBjZIOg44uoNtXgSGtOjGyvs1cJSkU9NghiGS9mllvV8A50o6SJkBkt4paQtJu0o6IiXaNWRHEE3txHS8sgEmvYFvkJ2nmQv8CdhF0gck9UqPAyTt3kY51wH/LWmHlBC/BVwfbY9qvBY4UdIx6Sipr7IBMtu1sX5LDwMr0wCQfqmMPSUdACDpdEnDIqIJWJ62aSI7N7aW7Hxc/xRnZ31c0naStibrnr2+k+U9QNZF/Im0H4yntC84tplyIrOyiIiVwKeAG8i+0b8PuKWDbZ4i+8Cfk7rDRrRY/jxZ9+R5ZF18U3jtvEl+vQbgv4BLUt2zgDPT4j7AxWQDFhYCbyI7X9WW3wBfTfXtT3buqvn1HU12zumFVFbzwJXWXEF2nudessERa4BPtlVpSpbjyboDF5MdYX2eEt+z6bzSCcA+qb4lwOVkR1gAxwLTJL1CNohlQuoSvoas+28+8CTwIJ33G7Lzo3PIjpo6NbowItaRHTmeRZaETyf7YrG2c2FaUen13cxmZl1H2Q/mz46IO8tcz0PApRFxZTnrsdrkIzIzKxxJh6Xfr/WUdAawF/DXasdl1VH2RCZplKS70w8ip0n6dLnrNLPN3q7AY2Rdi+cBp0TEguqGZNVS9q5FSdsC20bEI5K2IPuR5skR8WRZKzYzs26h7EdkEbEgIh5J0yuB6WRXLTAzM+u0iv6OTNIYsl/9P9Ri/jlkV2NgwIAB+++2226VDMvMzKpg8uTJSyJiWGfLqdioxfTbmXvIfqj6+7bWq6+vj4aGhorEZGZm1SNpckTUd7xm+yoyalFSL7KrPPy6vSRmZma2sSoxalHAL4HpEfGDjtY3MzPbGJU4IjuU7Dp2Ryi7P9QUSe1dwdvMzKxkZR/sERH307kLjpqZmbXJV/YwM7NCcyIzM7NCcyIzM7NCcyIzM7NCcyIzM7NCcyIzM7NCcyIzM7NCcyIzM7NCcyIzM7NCcyIzM7NCcyIzM7NCcyIzM7NCcyIzM7NCcyIzM7NCcyIzM7NCcyIzM7NCcyIzM7NCcyIzM7NCK3sik3SFpEWSppa7LjMz634qcUR2FXBsBeoxM7NuqG7ixIllrWDixInPfe1rX+sDvG/ixIk/62j9yy67bOI555xT1pisFUuWwE9/CkOHwpVXbtrfSy6BSZNg2203vQyXXRtlFz3+jS17112hf//qvf+6qa997WsLJk6ceFmnC4qIsj+AMcDUdpafAzQADdtvv31YFXz3uxEQcfzxnfvbFWW47Noou+jxb8y63/1utd+B3RLQEF2QY3p2OhN2gYi4DLgMoL6+PqocTvf0oQ9lf086CQ4/fNP+HnBAVsZpp216GS67NsouevwbW3bz/m+FpCwplrkSaQzwp4jYs6N16+vro6GhoewxmZlZdUmaHBH1nS3Hw+/NzKzQKjH8/jrgAWBXSfMknVXuOs3MrPso+zmyiDit3HWYmVn35a5FMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrtIokMknHSpohaZak8ytRp5mZdQ9lT2SS6oCfAscB44DTJI0rd71mZtY9VOKI7EBgVkTMiYh1wG+B8RWo18zMuoGeFahjJDA393wecFB+BUnnAOekp2slTa1AXOUyFFhS7SA2UZFjh2LHX+TYodjxO/bq2bUrCqlEIutQRFwGXAYgqSEi6qsc0iYrcvxFjh2KHX+RY4dix+/Yq0dSQ1eUU4muxfnAqNzz7dI8MzOzTqtEIvsXsLOkHST1BiYAt1SgXjMz6wbK3rUYERskfQK4HagDroiIae1sclm5YyqzIsdf5Nih2PEXOXYodvyOvXq6JH5FRFeUY2ZmVhW+soeZmRWaE5mZmRVaVRKZpPdImiapSVKbQ0fburRVGjjyUJp/fRpEUjGStpZ0h6SZ6e/gVtZ5u6QpuccaSSenZVdJeia3bJ9aij2t15iL75bc/Kq1fYntvo+kB9L+9bik9+aWVaXdO7pEm6Q+qS1npbYdk1t2QZo/Q9IxlYi3RWwdxf5ZSU+mtv67pNG5Za3uQ5VUQvxnSlqci/Ps3LIz0r42U9IZlY28pNh/mIv7aUnLc8uq2vaSrpC0SG38JliZn6TX9rik/XLLNr7dI6LiD2B3sh/CTQLq21inDpgNjAV6A48B49KyG4AJafpS4KMVjv+7wPlp+nzgOx2svzWwDOifnl8FnFKlti8pduCVNuZXre1LiR3YBdg5TY8AFgCDqtXu7e3HuXU+BlyapicA16fpcWn9PsAOqZy6Gov97bn9+qPNsbe3D9VY/GcCl7Sy7dbAnPR3cJoeXEuxt1j/k2QD6Wql7d8G7AdMbWP58cBtgICDgYc60+5VOSKLiOkRMaOD1Vq9tJUkAUcAN6b1rgZOLl+0rRqf6i21/lOA2yJiVVmjKs3Gxv5vNdD2HcYeEU9HxMw0/QKwCBhWsQjfqJRLtOVf143AkamtxwO/jYi1EfEMMCuVVykdxh4Rd+f26wfJfidaKzpzebxjgDsiYllEvATcARxbpjhbs7GxnwZcV5HIShAR95J9eW/LeOCayDwIDJK0LZvY7rV8jqy1S1uNBIYAyyNiQ4v5lTQ8Ihak6YXA8A7Wn8Abd7KL0iH1DyX16fII21Zq7H0lNUh6sLlLlOq3/Ua1u6QDyb7Nzs7NrnS7t7Uft7pOatsVZG1dyrbltLH1n0X2LbtZa/tQJZUa/3+mfeJGSc0XbyhM26fu3B2Au3Kzq932HWnr9W1Su5ftd2SS7gS2aWXRhRHxx3LV21Xaiz//JCJCUpu/YUjfMt5M9ju6ZheQfRD3JvsdxReBr3c25lydXRH76IiYL2kscJekJ8g+YMuqi9v9V8AZEdGUZpe13bszSacD9cBhudlv2IciYnbrJVTNrcB1EbFW0kfIjoyPqHJMG2sCcGNENObmFaHtu0zZEllEHNXJItq6tNVSssPQnunba1kuedVe/JJelLRtRCxIH5iL2inqVODmiFifK7v5qGKtpCuBz3VJ0K+V3+nYI2J++jtH0iRgX+Amytz2XRG7pC2BP5N9aXowV3ZZ270NpVyirXmdeZJ6AluR7efVvrxbSfVLOorsi8ZhEbG2eX4b+1AlP0w7jD8iluaeXk52HrZ528NbbDupyyNs28b87ycAH8/PqIG270hbr2+T2r2WuxZbvbRVZGcE7yY77wRwBlDpI7xbUr2l1P+Gvuv0Idx8zulkoJJX++8wdkmDm7vdJA0FDgWerIG2LyX23sDNZP3vN7ZYVo12L+USbfnXdQpwV2rrW4AJykY17gDsDDxcgZibdRi7pH2BnwMnRcSi3PxW96GKRZ4pJf5tc09PAqan6duBo9PrGAwczet7VcqtpEv7SdqNbFDEA7l5tdD2HbkF+GAavXgwsCJ90dy0dq/kSJbciJV3kfV9rgVeBG5P80cAf2kxsuVpsm8SF+bmjyV7Q88Cfgf0qXD8Q4C/AzOBO4Gt0/x64PLcemPIvmH0aLH9XcATZB+k1wIDayl24D9SfI+lv2fVQtuXGPvpwHpgSu6xTzXbvbX9mKxL86Q03Te15azUtmNz216YtpsBHFfJ/bzE2O9M7+Hmtr6lo32oxuL/NjAtxXk3sFtu2w+n/8ks4EO1Fnt6PhG4uMV2VW97si/vC9J7cR7Z+dNzgXPTcpHdcHl2irE+t+1Gt7svUWVmZoVWy12LZmZmHXIiMzOzQnMiMzOzQnMiMzOzQnMiMzOzQnMiMzOzQnMiMzOzQnMiMzOzQnMiMzOzQnMiMzOzQnMis82epEsl/U+J614l6ZvljinVNUnS2Wn6/ZL+1oVlT5N0eJqeKOnaLiz7S5Iu76ryzDrLicxqjqQLJN3WYt7MNuZN6Ki8iDg3Ir7RRbGFpJ26oqy8iPh1RBxdQv0lJdqI2CMiJnU2LkmHS5rXouxvRcTZnS3brKs4kVktuhf4D0l18O9bbfQC9m0xb6e0riXpfmZm3YoTmdWif5Elrn3S87eS3WJjRot5syPiBcjuyyTpDknLJM2QdGpzYS2PYiR9QdICSS9IOruVo6zBkv4saaWkhyTtmLZrTpqPSXpF0nslDZX0J0nLU933SWr1fSXpHZKekrRC0iVkt7JoXnampPvTtCT9UNIiSS9LekLSnpLOAd4PfCHVf2ta/1lJX5T0OPCqpJ5pXv4mpX0lXZ9e0yOS9s7V/brX39xekgYAtwEjUn2vSBrRsqtS0kmpK3N56i7dPbfsWUmfk/R4et3XS+rbWvuYbSonMqs5EbEOeAh4W5r1NuA+4P4W8+4FSB+4dwC/Ad5EdhPCn0ka17JsSccCnwWOIjuiO7yVECYAXyO7YeEs4KIUV3Pde0fEwIi4HjiP7H5Lw4DhwJeAN9wbKd3g8PfAl4GhZPdhOrSNJjg6vb5dyO4WfSqwNCIuA34NfDfVf2Jum9OAdwKDIrt7d0vjye55tjVZO/1BUq826ie93leB44AXUn0Dm7845F7XLmT3nvpMaoO/ALemm0E2OxU4FtgB2As4s716zTaWE5nVqnt4LWm9lSyR3ddi3j1p+gTg2Yi4MiI2RMSjwE3Ae1op91TgyoiYFhGryG5M2NLNEfFwSgi/5rWjwNasB7YFRkfE+oi4L1q/yd/xwLSIuDEi1gM/Aha2U+YWwG6AImJ6ZHfPbc9PImJuRKxuY/nkXN0/ILuZ58EdlFmK9wJ/jog7UtnfB/qR3dwxH9sLEbEMuJX229NsozmRWa26F3iLpK2BYRExE/gn2bmzrYE9ee382GjgoNS1tVzScrIuuG1aKXcEMDf3fG4r6+QTzCpgYDtxfo/sqO1vkuZIOr+N9V5Xb0p2rdVNRNwFXEJ2B91Fki6TtGU7MdBWWa0tj4gmsqPIER1sU4oRwHMtyp4LjMytszHtabbRnMisVj1A1q32X8A/ACLiZeCFNO+FiHgmrTsXuCciBuUeAyPio62UuwDYLvd8VGeCjIiVEXFeRIwFTgI+K+nINur9d12S1F7dEfGTiNgfGEfWxfj55kVtbdJBqPm6e5C1QXM34Sqgf27d/BeAjsp9geyLRHPZza9rfgfbmXUZJzKrSamLrIHsfNZ9uUX3p3n50Yp/AnaR9AFJvdLjgPygg5wbgA9J2l1Sf6Ck35flvAiMbX4i6QRJO6UP8BVAI9DUynZ/BvaQ9O40svBTtH7ESIr9oHQO61VgTa7M19W/EfbP1f0ZYC3wYFo2BXifpLp0DvGwFq93iKSt2ij3BuCdko5M8Z6Xyv7nJsRotkmcyKyW3UM2eOP+3Lz70rx/J7KIWEk2QGIC2RHCQuA7QJ+WBUbEbcBPyEZBzuK1D/O1JcY0Ebg6dWGeCuwM3Am8QnYU+bOIuLuVepeQnbO7GFiatvtHG3VsCfwCeIms224pWRcmwC+Bcan+P5QYM8Afyc5nvQR8AHh3OqcF8GngRKC5S/bf5UbEU2SDOeakOl/XHRkRM4DTgf8DlqRyTkwDdswqQq2flzbrHtJR21SgTxuj/cysxvmIzLodSe+S1EfSYLIjt1udxMyKy4nMuqOPAIvIfsvVCLQ2KMTMCsJdi2ZmVmg+IjMzs0KruQuMDh06NMaMGVPtMMzMrMwmT568JCKGdbacmktkY8aMoaGhodphmJlZmUl6ruO1OlZzicysO1mzvpGHnlnGvU8v5p+zl7JF354ctssw3rrzUPYcsRU9eqjjQsy6OScysyqICH5x3xz+929Ps3ZDE7179uCAMYNZvmo937t9Bt+7fQbbbNmXH753Hw7ZcUi1wzWraU5kZhX26toNfOGmx/nz4wt4x7jhnH7waA4cszX9etcBsOSVtdw/cwn/d9dMTv/lQ1xw3G6c9ZYdyK6CZWYtOZGZVdCzS17lI7+azMxFKzn/uN34yNvGviFBDR3Yh5P3HcmRu7+J8254jG/+eTpPzF/Bxe/e69/Jzsxe40RmViEvLF/NKZc+wIamJq7+8IG8def2B2tt0bcXl56+Pz+bNIv/veNplr26jivPPICedf7VjFme3xFmFbB6XSPn/KqBNesb+d1HDukwiTXr0UN84oidufjdb+a+mUv49m1PlTlSs+LxEZlZmUUEX7jpcaa98DKXf7CenYdvsdFlvPeA7Zm+YCW/vP8ZdttmC95T36nbqJltVnxEZlZmP5s0m1sfe4HPH7MrR+4+fJPL+fI7d+fQnYZw4c1TmfzcS10YoVmxOZGZldF9Mxfz/b/N4KS9R/DRw3bsVFk963pwyWn7se2gvpx77WSWvlLqLdTMNm9OZGZlsnLNer544+OMHTqA756yV5cMnx88oDeXnr4/y1etY+KtT3ZBlGbF50RmViYX3/YUC15ew/feszd9e3XdsPndt92STx6xM7c+9gK3T1vYZeWaFZUTmVkZ/HP2En790POcdegO7Lf94C4v/6OH78ju227Jl/8wlRWr1nd5+WZF4kRm1sVWrdvA+Tc9wZgh/Tnv6F3LUkevuh5875S9WPbqOr7+J3cxWvfmRGbWxb53+wyeX7aK7/xnea/EsefIrfjoYTty0yPzuHvGorLVY1brnMjMutDU+Su4+p/P8oGDR3PQ2PJf7PeTR+7ETm8ayFf+OJU16xvLXp9ZLXIiM+siTU3BV/44la0H9OZzx5SnS7GlPj3r+Pr4PZi7bDWX3jO7InWa1ZqSEpmkYyXNkDRL0vmtLO8j6fq0/CFJY9L8MZJWS5qSHpd2bfhmtePGR+bxyPPLOf+43dmqX6+K1fsfOw7lxL1H8LNJs3l+6aqK1WtWKzpMZJLqgJ8CxwHjgNMkjWux2lnASxGxE/BD4Du5ZbMjYp/0OLeL4jarKStWrefi256ifvRg3r3vyIrXf+Hxu9Orh/jardMqXrdZtZVyRHYgMCsi5kTEOuC3wPgW64wHrk7TNwJHyjdPsm7k+3+bwfJV6/j6+D2rclfnbbbqy2eO2oW/P7WIO598seL1m1VTKYlsJDA393xemtfqOhGxAVgBNJ/p3kHSo5LukfTW1iqQdI6kBkkNixcv3qgXYFZtU+ev4NqHnuODh4xh3IgtqxbHmYeOYec3DWTirdM88MO6lXIP9lgAbB8R+wKfBX4j6Q3v9Ii4LCLqI6J+2LDSbm9hVguamoIv/2EqQwb05r/fsUtVY+lV14Ovj9+TeS+t5meTPPDDuo9SEtl8IH/PiO3SvFbXkdQT2ApYGhFrI2IpQERMBmYD1X23m3WhGyfPY8rcyg/waMshOw7hpL1HcOk9s3lu6avVDsesIkpJZP8Cdpa0g6TewATglhbr3AKckaZPAe6KiJA0LA0WQdJYYGdgTteEblZdy1et4+K/Vm+AR1sufGc28GPiLdOIiGqHY1Z2HSaydM7rE8DtwHTghoiYJunrkk5Kq/0SGCJpFlkXYvMQ/bcBj0uaQjYI5NyIWNbVL8KsGv73b09XdYBHW4ZvmQ38uHvGYu6c7it+2OZPtfaNrb6+PhoaGqodhlm7ps5fwYmX3M8Zh4xh4kl7VDucN1jf2MTxP76P1esbufOzh3Xp1ffNuoqkyRFR39lyfGUPs43UWEMDPNrSq64HXxu/B/NeWs0ld82qdjhmZeVEZraRrv7ns0yZu5wvv3NcTQzwaMt/7DiUd+87kkvvmc30BS9XOxyzsnEiM9sIc5et4nu3z+Dtuw5j/D4jqh1Oh/7nhCzZfvGmx9nQ2FTtcMzKwonMrEQRwZdufoIegm++680U4eI1gwf05qsn7cHj81Zw5T+erXY4ZmXhRGZWopsemc99M5fwxeN2Y+SgftUOp2Qn7rUtR+72Jv73jhn+bZltlpzIzEqwaOUavvGnJ6kfPZjTDxpd7XA2iiS++a496dmjBxf8/gmammprpLJZZzmRmXWgqSn47PWPsWZ9Ixf/51419ZuxUm27VT++dPzu/HP2Un5xn69JYJsXJzKzDvy/e2Zz/6wlTDxpD3Z608Bqh7PJTjtwFMftuQ3fu30Gjzz/UrXDMesyTmRm7Wh4dhk/uONpTtx7BBMOGNXxBjVMEhf/515ss1VfPnXdo6xYvb7aIZl1CScyszYsX7WOT133KNsN7se33rVnIUYpdmSrfr34v9P2ZeGKNZx/0+O+FqNtFpzIzFqxobGJ/75+CotfWcv/nbYvW/St3R8+b6x9tx/M54/ZldumLuSX9z9T7XDMOs2JzKyFiOB//jiVu2cs5qsn7sFe2w2qdkhd7r/eOpZj9hjORX+Zzp8ef6Ha4Zh1ihOZWQs/unMm1z08l4+/fUdOP7hYQ+1L1aOH+PGEfdl/+8F89vrHeGD20mqHZLbJnMjMcn790HP8+O8zOWX/7fjc0btWO5yy6turjsvPqGf7If0555oGX4/RCsuJzCy54V9z+Z8/TOXwXYfx7XcX4xJUnTWof2+u/vCBDOjTkw/88mGmzl9R7ZDMNpoTmXV7EcH3b5/BF256nEN3GsrP3r8fveq6z1tj5KB+XHv2gfSuE6f+/AHufso347Ri6T7vVrNWrN3QyH9fP4VL7p7FhANGccWZB9C/d89qh1VxO71pC27++KHsMHQAZ1/TwK8feq7aIZmVzInMuq2ZL67kvT9/kD9MeYHPH7Mr3373m7vVkVhLw7fsyw0fOYS37TyUC2+eyhdvfNw/mrZC6L7vWuu21m1o4kd3Ps3xP7mP55a+ys/evx8ff/tO3eKcWEcG9OnJLz5Yz7mH7cjvJs/lHT+4h79OXVjtsMzapVr7ZX99fX00NDRUOwzbDDU2BbdPW8iP7nyap198hZP2HsFXTxzHkIF9qh1aTXpi3gq+cNPjTF/wMu8YN5xPHrHTZvmbOqseSZMjor7T5TiR2eZu1boN3DR5Hpff/wzPLV3F6CH9+coJ4zhy9+HVDq3mrW9s4rJ753DppNmsXLuBQ8YO4ZzDxnLYzsMKeRcAqy0VTWSSjgV+DNQBl0fExS2W9wEqhDz0AAAXcElEQVSuAfYHlgLvjYhn07ILgLOARuBTEXF7e3U5kVlXWPLKWu6avog7pr/IfTMXs2Z9E/uMGsRH3jaWo/fYhjp/CG+Ul9es57cPP88V9z/LwpfXMGyLPhy1+5s4avfhHLrTUPr2qqt2iFZAFUtkkuqAp4F3APOAfwGnRcSTuXU+BuwVEedKmgC8KyLeK2kccB1wIDACuBPYJSIa26rPicw2xorV65n/0mrmL1/NzEUrmTp/BVPnv8zzy1YBMGKrvrxj3HBO3HsE+48e7PNgnbRuQxO3T1vIX6ct5J4Zi3ll7QZ69hC7DN+CN4/cij1HbsmYoQMYOagfIwb1c4KzdnVVIitlnPGBwKyImJMq/i0wHngyt854YGKavhG4RNknxnjgtxGxFnhG0qxU3gOdDbwty1et489PLChX8bYR8t+RIs2I3LKmiNf9bYygsSl7bGhsYl1jsL6xiTXrG1m9rpHV6xt5Ze0Glq9az/LV63jp1fW8snbD6+octXU/3jxyKyYcOIrDdhnGuG23dPLqQr179uDEvUdw4t4jWLuhkQfnLOPBOUuZOn8Ff3tyIdc3zH3d+oP792Jw/95s1b8Xg/r1on+fnvTrVUf/3nX06dmDXnU96FnXg149RF2d6CFRJ9GjhxDQQ9ntZyQQQPpf5v+j/vdWxwFjtmaX4VtUOwygtEQ2EsjvnfOAg9paJyI2SFoBDEnzH2yx7ciWFUg6BzgHYPvtty819la9+PJaLrx5aqfKsNrQu64HPetEv1519Otdl30A9unJkIG92XHYAAb1782IQX0ZOag/Iwf3Y8yQ/gzq37vaYXcbfXrWcdguwzhsl2FA9sPyBSvW8PyyVf8+Sn7x5TUsX72eFavWs/iVtaxatoo16xpZtb6Rteub2NDUxPrG2jpPb6X5xsl7FiqRlV1EXAZcBlnXYmfKGjtsAA9/6cguicu6gPKTuW/WQA9l38AR1PUQPZTN69lD1PWQj6QKRhIjUpfixogI1jcGTZE9GpuCpiYIckfsvHaE/9pxPeAcWDUD+9ZE+gBKS2TzgfytcbdL81pbZ56knsBWZIM+Stm2S/Wq68GbtuxbzirMrAtJondPf2mxTVfKYI+eZIM9jiRLQv8C3hcR03LrfBx4c26wx7sj4lRJewC/4bXBHn8Hdm5vsIekxUBXXB9nKLCkC8rZnLhNWud2aZ3b5Y3cJq3b1HYZHRHDOlt5h0dk6ZzXJ4DbyYbfXxER0yR9HWiIiFuAXwK/SoM5lgET0rbTJN1ANjBkA/Dx9pJY2qbTLwpAUkNXjIbZnLhNWud2aZ3b5Y3cJq2rdruU1MkZEX8B/tJi3ldy02uA97Sx7UXARZ2I0czMrE2+1qKZmRXa5pzILqt2ADXIbdI6t0vr3C5v5DZpXVXbpeautWhmZrYxNucjMjMz6wacyMzMrNA2u0Qm6VhJMyTNknR+teOpFkmjJN0t6UlJ0yR9Os3fWtIdkmamv4OrHWulSaqT9KikP6XnO0h6KO0z10vqdte5kjRI0o2SnpI0XdIh3ldA0n+n989USddJ6tsd9xdJV0haJGlqbl6r+4cyP0nt87ik/cod32aVyNKV+n8KHAeMA05LV+DvjjYA50XEOOBg4OOpLc4H/h4RO5P9QL07JvtPA9Nzz78D/DAidgJeIrvtUHfzY+CvEbEbsDdZ+3TrfUXSSOBTQH1E7En2O9oJdM/95Srg2Bbz2to/jgN2To9zgP9X7uA2q0RG7kr9EbEOaL5Sf7cTEQsi4pE0vZLsg2kkWXtcnVa7Gji5OhFWh6TtgHcCl6fnAo4gu2sDdM822Qp4G9mFDYiIdRGxnG6+ryQ9gX7pCkf9gQV0w/0lIu4lu9hFXlv7x3jgmsg8CAyStG0549vcEllrV+p/w9X2uxtJY4B9gYeA4RHRfJ+bhUB3u03yj4AvAE3p+RBgeUQ03w+mO+4zOwCLgStTl+vlkgbQzfeViJgPfB94niyBrQAm4/2lWVv7R8U/hze3RGYtSBoI3AR8JiJezi+L7LcX3eb3F5JOABZFxORqx1JjegL7Af8vIvYFXqVFN2J321cA0jmf8WSJfgQwgDd2rxnV3z82t0RW8avt1zJJvciS2K8j4vdp9ovNh/np76JqxVcFhwInSXqWrNv5CLJzQ4NS1xF0z31mHjAvIh5Kz28kS2zdeV8BOAp4JiIWR8R64Pdk+1B331+atbV/VPxzeHNLZP8Cdk6jinqTnZi9pcoxVUU69/NLYHpE/CC36BbgjDR9BvDHSsdWLRFxQURsFxFjyPaNuyLi/cDdwClptW7VJgARsRCYK2nXNOtIsgt9d9t9JXkeOFhS//R+am6Xbr2/5LS1f9wCfDCNXjwYWJHrgiyLze7KHpKOJzsP0nyl/m55wWJJbwHuA57gtfNBXyI7T3YDsD3Z7XJOjYiWJ3E3e5IOBz4XESdIGkt2hLY18ChwekSsrWZ8lSZpH7IBML2BOcCHyL7odut9RdLXgPeSjQJ+FDib7HxPt9pfJF0HHE52u5YXga8Cf6CV/SMl/UvIumFXAR+KiIayxre5JTIzM+teNreuRTMz62acyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyKxwJH1J0uUlrHeVpG9WIqZcne+SNFfSK5L2rWTdtUzS4ZLm5Z5PS9e7NOs0JzLrcpKelbQ6fZi/mBLKwE0s63UfgAAR8a2IOLtrou1y3wc+EREDI+LR1BZHVTuoWhMRe0TEpGrHYZsHJzIrlxMjYiDZfa3qgS9vbAG5ez4VyWhgWlcUlG6D0e57tKBtZNalnMisrNLt4m8D9gSQ9CFJ0yWtlDRH0kea120++pL0RUkLgevStiPS0d0rkkZImijp2tx2b5H0T0nLU7fema3FIukESVPSev+UtFdu2RclzU9xzZB0ZBtlvFPSo5JeTnVNTPP7SHqF7PZBj0maLelXZLe4uDXF/oW07sG5eB/Ld7FJmiTpIkn/ILsFxthWYng2xfs48KqknqldbpK0WNIzkj6VW/9ASQ0p5hcl/SC37HeSFkpaIeleSXvkll0l6WeSbkvx/0PSNpJ+JOklSU/lu09TXBdIejItv1JS3zba8d9Hqun/eYOka1L7T5NUn1t3v9TmK1O811e6y9hqXET44UeXPoBngaPS9CiyI5RvpOfvBHYEBBxG9mG9X1p2ONl9n74D9AH6pXnzWpQ/Ebg2TY8GVgKnAb2AIcA+adlVwDfT9L5kd7A9iCzZnJHi7APsCswFRqR1xwA7tvHaDgfeTPYlcC+yezOdnFsewE6ttUV6PhJYChyfynhHej4sLZ9EdkPHPYCeQK822ndKatt+qZzJwFfI7ic2luyeYsek9R8APpCmBwIH58r6MLBFaocfAVNyy64ClgD7A32Bu4BngA+mNvwmcHeLuKamuLYG/pFr/9f9H3n9PjIRWJPapA74NvBgWtab7F5Xn07/33cD65rL9cOPiPARmZXNHyQtB+4H7gG+BRARf46I2ZG5B/gb8Nbcdk3AVyNibUSsLqGe9wF3RsR1EbE+IpZGxJRW1jsH+HlEPBQRjRFxNbAWOBhoJPsgHyepV0Q8GxGzW6ssIiZFxBMR0RQRj5MdNR5WSoMkpwN/iYi/pDLuABrIPsSbXRUR0yJiQ0Ssb6Ocn0TE3NRGB5Alwq9HxLqImAP8guwu2ADrgZ0kDY2IVyLiwdzruSIiVkZ2Y8iJwN6StsrVc3NETI6INcDNwJqIuCYiGoHryb4g5F2S4loGXET2BaMU96c2aQR+Beyd5h9MltB/kv6/vwceLrFM6yacyKxcTo6IQRExOiI+1pyUJB0n6UFJy1KiO57srrPNFqcPzVKNAlpNOi2MBs5L3XnLU92jyI7CZgGfIfsgXyTpt5JGtFaIpIMk3Z268FYA57aIv5Q43tMijrcA2+bWmVtCOfl1RpN1v+bL/BIwPC0/C9gFeErSvySdkF5LnaSLUzfoy2RHSbR4PS/mple38rzlIJ58XM8BrbZjKxbmplcBfdP5vxHA/IjI3wG4lPaxbsSJzCpGUh/gJrKRfcMjYhDwF7JuxmYtb1ne0S3M55J1VXZkLnBRSq7Nj/4RcR1ARPwmIt5ClhSCrHuzNb8BbgFGRcRWwKUt4m+pZfxzgV+1iGNARFzczjYdlTsXeKZFmVtExPHptc2MiNOAN6XXdaOkAWRHs+OBo4CtyLpU6eD1dGRUbnp74IVOlAWwABgpKR/TqLZWtu7JicwqqTdZF95iYIOk44CjO9jmRWBIi+6uvF8DR0k6NQ16GCJpn1bW+wVwbjqikqQBaeDGFpJ2lXRESrRryI40mtqobwtgWUSskXQgWTLoKP78gI1rgRMlHZOOiPoqG+SyXQfltOdhYGUaANIvlbunpAMAJJ0uaVhENAHL0zZN6bWsJTtH15/U/dtJH5e0naStgQvJuh874wGyrt9PpP/veODAzgZpmxcnMquYiFgJfAq4AXiJLAnc0sE2T5Gdh5qTus1GtFj+PFn35HnAMrJBEHu3Uk4D8F/AJanuWcCZaXEf4GKygQ0LyY5cLmgjpI8BX5e0kmxwxQ3txU82cOHLKfbPRcRcsqOgL5El9LnA5+nEezGdVzoB2IdsMMYS4HKyoyyAY4FpykZV/hiYkLp6ryHr/psPPAk8SOf9huy85xyyLt9OjS6MiHVkAzzOIkvCpwN/IkvAZgDo9V3PZmabRtKzwNkRcWeZ63kIuDQirixnPVYcPiIzs5om6bD0+7Weks4g+9nDX6sdl9WOsicySaPSKK8n0w8dP13uOs1ss7Ir8BhZ1+J5wCkRsaC6IVktKXvXoqRtgW0j4hFJW5D9cPPkiHiyrBWbmVm3UPYjsohYEBGPpOmVwHSyqxuYmZl1WkUvOCppDNmVAB5qMf8csisvMGDAgP132223SoZlZmZVMHny5CURMayz5VRs1KKy23jcQ/aj1N+3tV59fX00NDRUJCYzM6seSZMjor7jNdtXkVGLknqRXdHh1+0lMTMzs41ViVGLAn4JTI+IH3S0vpmZ2caoxBHZocAHgCOU3QtqiqTjO9rIzMysFGUf7BER99O5i5CamZm1yVf2MDOzQnMiMzOzQnMiMzOzQnMiMzOzQnMiMzOzQnMiMzOzQnMiMzOzQnMiMzOzQnMiMzOzQnMiMzOzQnMiMzOzQnMiMzOzQnMiMzOzQnMiMzOzQnMiMzOzQnMiMzOzQnMiMzOzQnMiMzOzQit7IpN0haRFkqaWuy4zM+t+elagjquAS4BrKlCXWfk88ABMmgRDhsDSpW/8u3w5TJmS/Z01C4YPh2eegQ0boF8/WL0a6upg/XqQoKkJevSAiGy6WUTVXqJZEZU9kUXEvZLGlLses7J64AE48khYuzZLOlKWcJr/tmbZstemV67M/m7Y8Pp1GhvfuF17ZZrZG9TEOTJJ50hqkNSwePHiaodj9kaTJsG6da8dOTUnGiccs6qriUQWEZdFRH1E1A8bNqza4Zi90eGHQ+/eWVcgZEdN+b9mVjWVOEdmVnyHHAJ//7vPkZnVICcys1Idckj2MLOaUonh99cBDwC7Spon6axy12lmZt1HJUYtnlbuOszMrPuqicEeZmZmm8qJzMzMCs2JzMzMCs2JzMzMCs2JzMzMCs2JzMzMCs2JzMzMCs2JzMzMCs2JzMzMCs2JzMzMCs2JzMzMCs2JzMzMCs2JzMzMCs2JzMzMCs2JzMzMCs2JzMzMCs2JzMzMCs2JzMzMCq0iiUzSsZJmSJol6fxK1GlmZt1D2ROZpDrgp8BxwDjgNEnjyl2vmZl1D5U4IjsQmBURcyJiHfBbYHwF6jUzs26gZwXqGAnMzT2fBxyUX0HSOcA56elaSVMrEFe5DAWWVDuITVTk2KHY8Rc5dih2/I69enbtikIqkcg6FBGXAZcBSGqIiPoqh7TJihx/kWOHYsdf5Nih2PE79uqR1NAV5VSia3E+MCr3fLs0z8zMrNMqkcj+BewsaQdJvYEJwC0VqNfMzLqBsnctRsQGSZ8AbgfqgCsiYlo7m1xW7pjKrMjxFzl2KHb8RY4dih2/Y6+eLolfEdEV5ZiZmVWFr+xhZmaF5kRmZmaFVpVEJuk9kqZJapLU5tDRti5tlQaOPJTmX58GkVSMpK0l3SFpZvo7uJV13i5pSu6xRtLJadlVkp7JLdunlmJP6zXm4rslN79qbV9iu+8j6YG0fz0u6b25ZVVp944u0SapT2rLWaltx+SWXZDmz5B0TCXibRFbR7F/VtKTqa3/Lml0blmr+1AllRD/mZIW5+I8O7fsjLSvzZR0RmUjLyn2H+biflrS8tyyqra9pCskLVIbvwlW5ifptT0uab/cso1v94io+APYneyHcJOA+jbWqQNmA2OB3sBjwLi07AZgQpq+FPhoheP/LnB+mj4f+E4H628NLAP6p+dXAadUqe1Lih14pY35VWv7UmIHdgF2TtMjgAXAoGq1e3v7cW6djwGXpukJwPVpelxavw+wQyqnrsZif3tuv/5oc+zt7UM1Fv+ZwCWtbLs1MCf9HZymB9dS7C3W/yTZQLpaafu3AfsBU9tYfjxwGyDgYOChzrR7VY7IImJ6RMzoYLVWL20lScARwI1pvauBk8sXbavGp3pLrf8U4LaIWFXWqEqzsbH/Ww20fYexR8TTETEzTb8ALAKGVSzCNyrlEm3513UjcGRq6/HAbyNibUQ8A8xK5VVKh7FHxN25/fpBst+J1orOXB7vGOCOiFgWES8BdwDHlinO1mxs7KcB11UkshJExL1kX97bMh64JjIPAoMkbcsmtnstnyNr7dJWI4EhwPKI2NBifiUNj4gFaXohMLyD9Sfwxp3sonRI/UNJfbo8wraVGntfSQ2SHmzuEqX6bb9R7S7pQLJvs7Nzsyvd7m3tx62uk9p2BVlbl7JtOW1s/WeRfctu1to+VEmlxv+faZ+4UVLzxRsK0/apO3cH4K7c7Gq3fUfaen2b1O5l+x2ZpDuBbVpZdGFE/LFc9XaV9uLPP4mIkNTmbxjSt4w3k/2OrtkFZB/Evcl+R/FF4OudjTlXZ1fEPjoi5ksaC9wl6QmyD9iy6uJ2/xVwRkQ0pdllbffuTNLpQD1wWG72G/ahiJjdeglVcytwXUSslfQRsiPjI6oc08aaANwYEY25eUVo+y5TtkQWEUd1soi2Lm21lOwwtGf69lqWS161F7+kFyVtGxEL0gfmonaKOhW4OSLW58puPqpYK+lK4HNdEvRr5Xc69oiYn/7OkTQJ2Be4iTK3fVfELmlL4M9kX5oezJVd1nZvQymXaGteZ56knsBWZPt5tS/vVlL9ko4i+6JxWESsbZ7fxj5UyQ/TDuOPiKW5p5eTnYdt3vbwFttO6vII27Yx//sJwMfzM2qg7TvS1uvbpHav5a7FVi9tFdkZwbvJzjsBnAFU+gjvllRvKfW/oe86fQg3n3M6Gajk1f47jF3S4OZuN0lDgUOBJ2ug7UuJvTdwM1n/+40tllWj3Uu5RFv+dZ0C3JXa+hZggrJRjTsAOwMPVyDmZh3GLmlf4OfASRGxKDe/1X2oYpFnSol/29zTk4Dpafp24Oj0OgYDR/P6XpVyK+nSfpJ2IxsU8UBuXi20fUduAT6YRi8eDKxIXzQ3rd0rOZIlN2LlXWR9n2uBF4Hb0/wRwF9ajGx5muybxIW5+WPJ3tCzgN8BfSoc/xDg78BM4E5g6zS/Hrg8t94Ysm8YPVpsfxfwBNkH6bXAwFqKHfiPFN9j6e9ZtdD2JcZ+OrAemJJ77FPNdm9tPybr0jwpTfdNbTkrte3Y3LYXpu1mAMdVcj8vMfY703u4ua1v6WgfqrH4vw1MS3HeDeyW2/bD6X8yC/hQrcWenk8ELm6xXdXbnuzL+4L0XpxHdv70XODctFxkN1yenWKsz2270e3uS1SZmVmh1XLXopmZWYecyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrNCcyMzMrND+PwcpXEHjl/YTAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAF29JREFUeJzt3Xv8bfW87/HXW0Xk0m3pkVhWSHYcwhLh0Ymc2LQ3DuXYhfA4ySF1Hmcjtk0He+9w3B2XkHJX6CjsLnttKZud1mKlm+RRq13pTiK5VJ/zx/j+avr5Xea6zDmtOV7Px2M+5hhjjstnzPFb6z3H7TtSVUiS+usuky5AkjRZBoEk9ZxBIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HObTrqAYWy77ba1bNmySZchSRuVVatWXV9VSxYbb6MIgmXLlrFy5cpJlyFJG5Uklw0znoeGJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqec2ijuLpT9Xyw7/xsSWvebIZ01s2Zou7hFIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgST1nEEhSzxkEktRzBoEk9ZxBIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgST1nEEhSzxkEktRzIwuCJA9I8q0kFyQ5P8mhbfjWSU5LcnF732pUNUiSFjfKPYJbgf9VVbsATwBelWQX4HBgRVXtBKxo/ZKkCRlZEFTVVVX1g9b9K+BCYAfg2cCxbbRjgeeMqgZJ0uLGco4gyTLg0cBZwHZVdVX76Gpgu3HUIEma28iDIMk9ga8Ah1XVTYOfVVUBNc90ByVZmWTlddddN+oyJam3RhoESTajC4HPVdVX2+BrkmzfPt8euHauaavqqKpaXlXLlyxZMsoyJanXRnnVUIBPAhdW1XsGPjoReEnrfgnwtVHVIEla3KYjnPeTgBcB5yZZ3Ya9ETgSOC7Jy4HLgP1GWIMkaREjC4Kq+g6QeT7ea1TLlSStHe8slqSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ6b987iJFsvNGFV/XzDlyNJGreFmphYRddEdIClwC9a95bAfwA7jrw6SdLIzXtoqKp2rKoHAf8C/FVVbVtV2wD7AKeOq0BJ0mgNc47gCVX1zZmeqvpn4ImjK0mSNE7DtD76syRvAj7b+vcHfja6kiRJ4zTMHsELgSXACcBXW/cLR1mUJGl8Ft0jaFcHHZpki6q6eQw1SZLGaNE9giRPTHIBcGHrf1SSD4+8MknSWAxzaOi9wNOBGwCq6hxgj1EWJUkan6HuLK6qy2cNum0EtUiSJmCYq4YuT/JEoJJsBhxKO0wkSdr4DbNHcDDwKmAH4Epg19YvSZoCw1w1dD3dvQOSpCk0zFVD70xy7ySbJVmR5LokB4yjOEnS6A1zaGjvqrqJro2hNcBDgNeOsihJ0vgMEwQzh4+eBRxfVb8cYT2SpDEb5qqhryf5MXAL8MokS4DfjrYsSdK4LLpHUFWH07U2uryq/gDcDDx71IVJksZj0T2CJC8e6B786NOjKEiSNF7DHBp63ED35sBewA8wCCRpKgxzH8Ehg/1JtgS+OLKKJEljNVRbQ7PcjM8rlqSpMcw5gpPoHmIPXXDsAhw3yqIkSeMzzDmC/zPQfStwWVVdMaJ6JEljNsw5gm+PoxBJ0mSsyzkCSdIUMQgkqefWKgiSbJXkkUOOe3SSa5OcNzDsiCRXJlndXs9c24IlSRvWMM1Qn96aod6a7kayjyd5zxDzPgZ4xhzD31tVu7bXN9euXEnShjbMHsF9WjPU/xX4dFU9HnjaYhNV1RnAz9ezPknSiA3VDHWS7YH9gK9vgGW+OsmP2qGjrTbA/CRJ62GY+wjeCpwCfKeqzk7yIODidVzeR4C30d2g9jbg3cDL5hoxyUHAQQBLly5dx8WpL5Yd/o1JlyBttIZphvr4qnpkVf2P1n9JVT1vXRZWVddU1W1VdTvwcWC3BcY9qqqWV9XyJUuWrMviJElDmHePIMnrquqdST7InU1M3KGqXrO2C0uyfVVd1XqfC5y30PiSpNFb6NDQhe195brMOMkXgD2BbZNcAbwF2DPJrnTBsgZ4xbrMW5K04cwbBFV1Uns/dl1mXFUvnGPwJ9dlXpKk0fHOYknqOYNAknpunYIgyV03dCGSpMkYtomJZQP9uwFnj7AmSdIYDXND2T8BJyf5ALAD8JfAS0dalSRpbIZ5MM0pSQ4GTgOuBx5dVVePvDJJ0lgMc2jo74EPAnsARwCnJ3nWiOuSJI3JMIeGtgF2q6pbgO8lORn4BGDjLpI0BYY5NHQYQJJ7VNVvquoy4L+MvDJJ0lgMc2ho9yQXAD9u/Y9K8uGRVyZJGoth7iN4H/B04AaAqjqH7nyBJGkKDHVDWVVdPmvQbSOoRZI0AcOcLL48yROBSrIZcCh3tkwqSdrIDRMEBwPvp7uZ7ErgVOBVoyxKGyefEtYPk9zOa470yvVRGOaqoeuB/cdQiyRpAhZ6QtmcTyabsS5PKJMk/flZ6GTxSmAVsDnwGLoH1l8M7ArY+qgkTYmFnlB2LECSVwJPrqpbW/9HgTPHU54kadSGuXx0K+DeA/33bMMkSVNgmKuGjgR+mORbQLiz8TlJ0hQY5qqhTyX5Z+DxbdDrbYZakqbHsI+q3AS4DvgF8NAkNjEhSVNi0T2CJO8AXgCcD9zeBhdwxgjr0nrwxq5+cDtrQxnmHMFzgJ2r6nejLkaSNH7DHBq6BNhs1IVIkiZjmD2C3wCrk6wA7tgr8M5iSZoOwwTBie0lSZpCw1w+euw4CpEkTcZCjc4dV1X7JTmXORqfq6pHjrQySdJYLLRHcGh732cchUiSJmOhRueuau+Xja8cSdK4DXtnsSRpShkEktRzQwVBkrsn2XnUxUiSxm/RIEjyV8Bq4OTWv2sS7yuQpCkxzB7BEcBuwI0AVbUa2HGENUmSxmiYIPhDVf1y1rB5H2o/I8nRSa5Nct7AsK2TnJbk4vbuk84kacKGCYLzk/wNsEmSnZJ8EPjuENMdAzxj1rDDgRVVtROwovVLkiZomCA4BHg4XYNzXwBuAg5bbKKqOgP4+azBzwZmmqw4lq6Ja0nSBA3T1tBvgL9rr/W13cyNasDVwHYbYJ6SpPUwzBPKlgNvBJYNjr++bQ1VVSWZ91xDkoOAgwCWLl26PouSJC1gmGaoPwe8FjiXOx9Vua6uSbJ9VV2VZHvg2vlGrKqjgKMAli9fvujJaUnSuhkmCK6rqg1138CJwEuAI9v71zbQfCVJ62iYIHhLkk/QXeUz+ISyry40UZIvAHsC2ya5AngLXQAcl+TlwGXAfutYtyRpAxkmCF4KPIzuucUzh4YKWDAIquqF83y019DVSZJGbpggeFxV2c6QJE2pYe4j+G6SXUZeiSRpIobZI3gCsDrJpXTnCEJ39aePqpSkKTBMEMxuJkKSNEUWenj9vavqJuBXY6xHkjRmC+0RfJ7uwfWr6K4SysBnBTxohHVJksZkoYfX79PeffaAJE2xYZ5QtmKYYZKkjdNC5wg2B+5Bd2fwVtx5aOjewA5jqE2SNAYLnSN4Bd1zB+5Hd55gJghuAj404rokSWOy0DmC9wPvT3JIVX1wjDVJksZo0XMEhoAkTbdhmpiQJE2xeYMgyZPa+93GV44kadwW2iP4QHv/3jgKkSRNxkJXDf0hyVHADkk+MPvDqnrN6MqSJI3LQkGwD/A04Ol0l49KkqbQQpePXg98McmFVXXOGGuSJI3RMFcN3ZDkhCTXttdXktx/5JVJksZimCD4FHAi3R3G9wNOasMkSVNgmCC4b1V9qqpuba9jgCUjrkuSNCbDBMH1SQ5Iskl7HQDcMOrCJEnjMUwQvAzYD7gauAp4PvDSURYlSRqfRZ9ZXFWXAX89hlokSRNgW0OS1HMGgST1nEEgST03zDOL3zTQbUukkjRlFmqG+vVJdqe7SmiGLZFK0pRZ6KqhHwP7Ag9Kcmbr3ybJzlV10ViqkySN3EKHhm4E3gj8FNgTeH8bfniS7464LknSmCy0R/B04M3Ag4H3AD8Cbq4qbyaTpCky7x5BVb2xqvYC1gCfATYBliT5TpKTxlSfJGnEFr2zGDilqlYCK5O8sqqenGTbURcmSRqPRS8frarXDfQe2IZdP6qCJEnjtVY3lPmkMkmaPsMcGtrgkqwBfgXcBtxaVcsnUYckaUJB0DzFQ0ySNHm2NSRJPTepICjg1CSrkhw0oRokSUzu0NCTq+rKJPcFTkvy46o6Y3CEFhAHASxdunQSNUpSL0xkj6Cqrmzv1wInALvNMc5RVbW8qpYvWbJk3CVKUm+MPQiSbJHkXjPdwN7AeeOuQ5LUmcShoe2AE5LMLP/zVXXyBOqQJDGBIKiqS4BHjXu5kqS5efmoJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgST1nEEhSzxkEktRzBoEk9ZxBIEk9ZxBIUs8ZBJLUcwaBJPXcJB9ePxbLDv/GpEuQtIH08d/zmiOfNfJluEcgST1nEEhSzxkEktRzBoEk9ZxBIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgST1nEEhSzxkEktRzBoEk9ZxBIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPTSQIkjwjyUVJfprk8EnUIEnqjD0IkmwC/F/gL4FdgBcm2WXcdUiSOpPYI9gN+GlVXVJVvwe+CDx7AnVIkphMEOwAXD7Qf0UbJkmagE0nXcB8khwEHNR6f53koknWM4RtgesnXcSEuO791ef1H8u65x3rNfkDhxlpEkFwJfCAgf77t2F/pKqOAo4aV1HrK8nKqlo+6TomwXXv57pDv9d/mtZ9EoeGzgZ2SrJjkrsC/w04cQJ1SJKYwB5BVd2a5NXAKcAmwNFVdf6465AkdSZyjqCqvgl8cxLLHqGN5jDWCLju/dXn9Z+adU9VTboGSdIE2cSEJPWcQbCOkmyS5IdJvt76d0xyVms240vtRPhUmmPdj0lyaZLV7bXrpGsclSRrkpzb1nNlG7Z1ktOSXNzet5p0naMwz7ofkeTKgW3/zEnXOSpJtkzy5SQ/TnJhkt2nZdsbBOvuUODCgf53AO+tqocAvwBePpGqxmP2ugO8tqp2ba/VkyhqjJ7S1nPm0sHDgRVVtROwovVPq9nrDt3f/cy2n7Zzf4PeD5xcVQ8DHkX3b2Aqtr1BsA6S3B94FvCJ1h/gqcCX2yjHAs+ZTHWjNXvdBXRNpBzbuqd22/dZkvsAewCfBKiq31fVjUzJtjcI1s37gNcBt7f+bYAbq+rW1j/NzWbMXvcZ/5DkR0nem+RuE6hrXAo4Ncmqdvc7wHZVdVXrvhrYbjKljdxc6w7w6rbtj95YD40MYUfgOuBT7bDoJ5JswZRse4NgLSXZB7i2qlZNupZxW2Dd3wA8DHgcsDXw+nHXNkZPrqrH0LWe+6okewx+WN1leNN6Kd5c6/4R4MHArsBVwLsnWN8obQo8BvhIVT0auJlZh4E25m1vEKy9JwF/nWQNXcupT6U7drhlkpn7MuZsNmMK/Mm6J/lsVV1Vnd8Bn6JrYXYqVdWV7f1a4AS6db0myfYA7f3ayVU4OnOte1VdU1W3VdXtwMeZ3m1/BXBFVZ3V+r9MFwxTse0NgrVUVW+oqvtX1TK65jH+tar2B74FPL+N9hLgaxMqcWTmWfcDBv4hhO4Y6XkTLHNkkmyR5F4z3cDedOt6It02hynd9vOt+8y2b57LlG77qroauDzJzm3QXsAFTMm2/7NtfXQj9Hrgi0neDvyQdlKpJz6XZAkQYDVw8ITrGZXtgBO6vGNT4PNVdXKSs4HjkrwcuAzYb4I1jsp86/6ZdrlwAWuAV0yuxJE7hO5v/a7AJcBL6X5Mb/Tb3juLJannPDQkST1nEEhSzxkEktRzBoEk9ZxBIEk9ZxD0VJLbWmuR5yU5Psk91nL6wwanSfLNJFsuMP4RSf52PWu+X5IvLz7mH01zYJIPte6Dk7x4fWpYZDn3W8tpTk8yFc+8ndFaon3+4mPeMf6yJFN578HGxCDor1taa5GPAH7PWlz7n2QT4DDgjiCoqme2RrhGpqp+VlVD/yczx/QfrapPb8iaBhwIrFUQbCgDd7RL68QgEMCZwEMAkvy/1qjY+YMNiyX5dZJ3JzkH+Du6//S+leRb7fM1SbZt3S9ujZCdk+QzsxeW5MFJTm7LOTPJw9rwfdseyjlJzphjujt+PbZf4F9t87k4yTsHxntpkp8k+T5dsxgzw+/YK0nykCT/0pb1gyQPbsNfm+TsVv//bsO2SPKNNu55SV4wq67nA8vpbjZaneTuSfZqjZOd2xpjm68hvhcN7JntNrC8o5N8v83j2XN8F3u27+5EujtcSXJAm2Z1ko+le27EJu1X+nmtlv/Zxv3vbT3PSfKVmb27Nu5Hkvx7kkvaco5O1/7+MbP+Ht7b/k5WpLuhcHaNj03y7badT8mdd6A/ti33HOBV83wvGqeq8tXDF/Dr9r4p3W3xr2z9W7f3u9M1F7BN6y9gv4Hp1wDbzu4HHg78ZOazgfkdAfxt614B7NS6H0/XVAXAucAOrXvLOWpeBpzXug+ku7vzPsDmdHd1PgDYHvgPYAlwV+DfgA/NUcNZwHNb9+Z0ezd70z2HNnQ/kr5O1/Tw84CPD9RxnzlqOx1YPjC/y4GHtv5PA4fNM83HW/ceA+v2j8ABM99D+z63mDXtnnQNn+3Y+v8COAnYrPV/GHgx8FjgtIHptmzv2wwMeztwSOs+hq4dqdA1sXwT8J/a97EK2HXg72H/1v3mge/4GLqmVjYDvgssacNfABzdun8E7NG63zWz3r4m93KXsr/unmTmATJncmeTGK9J8tzW/QBgJ+AG4DbgK0PM96nA8VV1PUBV/XzwwyT3BJ4IHJ+uuQKAmV/L/wYck+Q44KtDLGtFVf2yzfcC4IF0YXR6VV3Xhn8JeOisGu5FFzgntBp/24bvTRcGP2yj3rOt/5nAu5O8A/h6VZ25SF07A5dW1U9a/7F0v3zfN8e4X2g1nJHk3unOs+xN17jfzDmVzYGl/OnDgL5fVZe27r3o/tM/u32vd6drAO0k4EFJPgh8Azi1jf+IdM2hbNnW85SB+Z5UVZXkXOCaqjq3fT/n04XxarpmyL/Uxv8sf7q9dgYeAZzW6tkEuKqt35ZVNbPH9xm61kw1QQZBf91SVX/0SMkkewJPA3avqt8kOZ3uPyGA31bVbRtguXehe3bDnzzOsqoOTvJ4ugffrEry2Kq6YYF5/W6g+zbW/+85wD9V1cf+5IPkMcAzgbcnWVFVb13PZc2Y3cZLtTqeV1UXLTLtzYMlAsdW1Rtmj5TkUcDT6c4D7Qe8jO6X+3Oq6pwkB9LtYcyY+V5v54+/49uZ/zuevR4Bzq+q3WfVMu8FBZoczxFo0H2AX7QQeBjwhAXG/RVwrzmG/yuwb5JtoHue7+CHVXUTcGmSfdvnaf9RkeTBVXVWVb2Z7iEgD1iHdTgL+M9JtkmyGbDv7BGq6lfAFUme05Z7t3aM/BTgZW2vhSQ7JLlvuquBflNVn6U7lPGYOZY7+H1cBCxL8pDW/yLg2/PU+4K2rCcDv2x7OKcAh6T9lE7y6CHWewXw/CT3bdNsneSB6c7b3KWqvgK8aaD2e9H9Qt8M2H+I+c92F+5sbfdvgO/M+vwiYEmS3Vs9myV5eHUXFNzY1pd1XLY2MPcINOhk4OAkF9L9Q/73BcY9Cjg5yc+q6ikzA6vq/CT/AHw7yW10h1kOnDXt/sBHkryJ7ljyF4FzgHcl2Ynu1+SKNmytVNVVSY4AvgfcSHcYYy4vAj6W5K3AH4B9q+rUJH8BfK/9H/xr4AC6E+nvSnJ7G/eVc8zvGOCjSW4BdqdrmfL4dFf0nA18dJ46fpvkh3Tfw8vasLfRHUb6UZK7AJcC+yyy3he07/PUNs0f6A5H3UL3VK2ZH30zewx/Txea17X3uUJ9ITcDu7VlXksLtIF6fp/uJPoH0j3mcdO2TufTfTdHJynuPFSlCbL1UUlrLcmvq+qek65DG4aHhiSp59wjkKSec49AknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ77/xD8j+K3VrtHAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# CODE SNIPPET #\n", + "def normalize_weight(particles):\n", + "\n", + " sumw = sum([p.w for p in particles])\n", + "\n", + " try:\n", + " for i in range(N_PARTICLE):\n", + " particles[i].w /= sumw\n", + " except ZeroDivisionError:\n", + " for i in range(N_PARTICLE):\n", + " particles[i].w = 1.0 / N_PARTICLE\n", + "\n", + " return particles\n", + "\n", + " return particles\n", + "\n", + "\n", + "def resampling(particles):\n", + " \"\"\"\n", + " low variance re-sampling\n", + " \"\"\"\n", + "\n", + " particles = normalize_weight(particles)\n", + "\n", + " pw = []\n", + " for i in range(N_PARTICLE):\n", + " pw.append(particles[i].w)\n", + "\n", + " pw = np.array(pw)\n", + "\n", + " Neff = 1.0 / (pw @ pw.T) # Effective particle number\n", + " # print(Neff)\n", + "\n", + " if Neff < NTH: # resampling\n", + " wcum = np.cumsum(pw)\n", + " base = np.cumsum(pw * 0.0 + 1 / N_PARTICLE) - 1 / N_PARTICLE\n", + " resampleid = base + np.random.rand(base.shape[0]) / N_PARTICLE\n", + "\n", + " inds = []\n", + " ind = 0\n", + " for ip in range(N_PARTICLE):\n", + " while ((ind < wcum.shape[0] - 1) and (resampleid[ip] > wcum[ind])):\n", + " ind += 1\n", + " inds.append(ind)\n", + "\n", + " tparticles = particles[:]\n", + " for i in range(len(inds)):\n", + " particles[i].x = tparticles[inds[i]].x\n", + " particles[i].y = tparticles[inds[i]].y\n", + " particles[i].yaw = tparticles[inds[i]].yaw\n", + " particles[i].w = 1.0 / N_PARTICLE\n", + "\n", + " return particles, inds\n", + "# END OF SNIPPET #\n", + "\n", + "\n", + "\n", + "def gaussian(x, mu, sig):\n", + " return np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))\n", + "N_PARTICLE = 100\n", + "particles = [Particle(N_LM) for i in range(N_PARTICLE)]\n", + "x_pos = []\n", + "w = []\n", + "for i in range(N_PARTICLE):\n", + " particles[i].x = np.linspace(-0.5,0.5,N_PARTICLE)[i]\n", + " x_pos.append(particles[i].x)\n", + " particles[i].w = gaussian(i, N_PARTICLE/2, N_PARTICLE/20)\n", + " w.append(particles[i].w)\n", + " \n", + "\n", + "# Normalize weights\n", + "sw = sum(w)\n", + "for i in range(N_PARTICLE):\n", + " w[i] /= sw\n", + "\n", + "particles, new_indices = resampling(particles)\n", + "x_pos2 = []\n", + "for i in range(N_PARTICLE):\n", + " x_pos2.append(particles[i].x)\n", + " \n", + "# Plot results\n", + "fig, ((ax1,ax2,ax3)) = plt.subplots(nrows=3, ncols=1)\n", + "fig.tight_layout()\n", + "ax1.plot(x_pos,np.ones((N_PARTICLE,1)), '.r', markersize=2)\n", + "ax1.set_title(\"Particles before resampling\")\n", + "ax1.axis((-1, 1, 0, 2))\n", + "ax2.plot(w)\n", + "ax2.set_title(\"Weights distribution\")\n", + "ax3.plot(x_pos2,np.ones((N_PARTICLE,1)), '.r')\n", + "ax3.set_title(\"Particles after resampling\")\n", + "ax3.axis((-1, 1, 0, 2))\n", + "fig.subplots_adjust(hspace=0.8)\n", + "plt.show()\n", + "\n", + "plt.figure()\n", + "plt.hist(new_indices)\n", + "plt.xlabel(\"Particles indices to be resampled\")\n", + "plt.ylabel(\"# of time index is used\")\n", + "plt.show()" ] }, { @@ -275,15 +631,10 @@ "source": [ "### References\n", "\n", - "http://www.probabilistic-robotics.org/" + "http://www.probabilistic-robotics.org/\n", + "\n", + "http://ais.informatik.uni-freiburg.de/teaching/ws12/mapping/pdf/slam10-fastslam.pdf" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": {