robot models.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: TkAgg\n"
     ]
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import random\n",
    "import math\n",
    "from matplotlib import patches\n",
    "%matplotlib\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figures for the motion models\n",
    "Using velocity motion model drawing some motions using particles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_vel_motion_model(control,priorpose,alphas):\n",
    "    delta_t = 1\n",
    "    \n",
    "    transvel = random.gauss(control[0],math.sqrt(alphas[0]*control[0]**2+alphas[1]*control[1]**2))\n",
    "    angvel   = random.gauss(control[1],math.sqrt(alphas[2]*control[0]**2+alphas[3]*control[1]**2))\n",
    "    gamma    = random.gauss(0,math.sqrt(alphas[4]*control[0]**2+alphas[5]*control[1]**2))\n",
    "    #print(transvel,angvel,gamma)\n",
    "    new_x = priorpose[0] - transvel/angvel * math.sin(priorpose[2]) + transvel/angvel * math.sin(priorpose[2]+angvel*delta_t)\n",
    "    new_y = priorpose[1] + transvel/angvel * math.cos(priorpose[2]) - transvel/angvel * math.cos(priorpose[2]+angvel*delta_t)\n",
    "    new_theta = priorpose[2] + angvel*delta_t + gamma*delta_t\n",
    "    \n",
    "    return [new_x,new_y,new_theta]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "particles = []\n",
    "high_transunc_particles = []\n",
    "high_rotunc_particles = []\n",
    "control = [3,1]\n",
    "alphas_transunc = [0.01,0.01,0.00001,0.0001,0.0001,0.0001]#high alpha 0,1 value for high translational error\n",
    "alphas_rotunc = [0.0001,0.0001,0.01,0.01,0.0001,0.0001] #high alpha 3,4 value for high rotational error\n",
    "for i in range(100):\n",
    "    particles.append([0,0,0])\n",
    "    high_transunc_particles.append(sample_vel_motion_model(control,particles[i],alphas_transunc))\n",
    "    high_rotunc_particles.append(sample_vel_motion_model(control,particles[i],alphas_rotunc))\n",
    "particles = np.array(particles)\n",
    "high_transunc_particles = np.array(high_transunc_particles)\n",
    "high_rotunc_particles = np.array(high_rotunc_particles)\n",
    "\n",
    "arrowlenght = 0.1\n",
    "high_transunc_particles_arrows = np.array([arrowlenght*np.cos(high_transunc_particles[:,2]),\n",
    "                                    arrowlenght*np.sin(high_transunc_particles[:,2])])\n",
    "high_rotunc_particles_arrows = np.array([arrowlenght*np.cos(high_rotunc_particles[:,2]),\n",
    "                                    arrowlenght*np.sin(high_rotunc_particles[:,2])])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2.5244129544236893, 1.3790930823955807, 1.0]\n"
     ]
    }
   ],
   "source": [
    "#calculating the ideal new pose (no errors!)\n",
    "ideal_new_pos = sample_vel_motion_model([3,1],[0,0,0],[0,0,0,0,0,0])\n",
    "print(ideal_new_pos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using odometry motion model drawing some motions using particles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_odo_motion_model(control,priorpos,alphas):\n",
    "    #convert control\n",
    "    delta_rot1 = math.atan2(control[1][1]-control[0][1],control[1][0]-control[0][0]) - control[0][2]\n",
    "    delta_trans = math.sqrt((control[1][0]-control[0][0])**2+(control[1][1]-control[0][1])**2)\n",
    "    delta_rot2 = control[1][2] - control[0][2] - delta_rot1\n",
    "    #sample deltas\n",
    "    sampled_delta_rot1 = random.gauss(delta_rot1,math.sqrt(alphas[0]*delta_rot1**2+alphas[1]*delta_trans**2))\n",
    "    sampled_delta_trans = random.gauss(delta_trans,math.sqrt(alphas[2]*delta_trans**2+alphas[3]*(delta_rot1**2+delta_rot2**2)))\n",
    "    sampled_delta_rot2 = random.gauss(delta_rot2,math.sqrt(alphas[0]*delta_rot2**2+alphas[1]*delta_trans**2))\n",
    "    #calc new pos\n",
    "    new_x = priorpos[0] + sampled_delta_trans*math.cos(priorpos[2]+sampled_delta_rot1)\n",
    "    new_y = priorpos[1] + sampled_delta_trans*math.sin(priorpos[2]+sampled_delta_rot1)\n",
    "    new_theta = priorpos[2] + delta_rot1 + delta_rot2\n",
    "    \n",
    "    return [new_x,new_y,new_theta]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "particles_2 = []\n",
    "high_transunc_particles_2 = []\n",
    "high_rotunc_particles_2 = []\n",
    "control_2 = [[0,0,0],[2.5244129544236893, 1.3790930823955807, 1.0]]\n",
    "alphas_transunc_2 = [0.0001,0.00001,0.01,0.01]#high alpha 2,3 value for high translational error\n",
    "alphas_rotunc_2 = [0.01,0.001,0.0001,0.0001] #high alpha 1,2 value for high rotational error\n",
    "for i in range(100):\n",
    "    particles_2.append([0,0,0])\n",
    "    high_transunc_particles_2.append(sample_odo_motion_model(control_2,particles_2[i],alphas_transunc_2))\n",
    "    high_rotunc_particles_2.append(sample_odo_motion_model(control_2,particles_2[i],alphas_rotunc_2))\n",
    "\n",
    "particles_2 = np.array(particles_2)\n",
    "high_transunc_particles_2 = np.array(high_transunc_particles_2)\n",
    "high_rotunc_particles_2 = np.array(high_rotunc_particles_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: Qt5Agg\n"
     ]
    }
   ],
   "source": [
    "#setting up the plots\n",
    "%matplotlib\n",
    "plt.rcParams.update({'font.size': 13})\n",
    "plt.rcParams['figure.figsize'] = 10, 10\n",
    "\n",
    "fig, ((ax1, ax2),(ax3,ax4)) = plt.subplots(nrows=2,ncols=2)\n",
    "ax1.set_aspect('equal')\n",
    "ax2.set_aspect('equal')\n",
    "ax3.set_aspect('equal')\n",
    "ax4.set_aspect('equal')\n",
    "\n",
    "#print(ax1)\n",
    "\n",
    "#for i in range(len(particles)):\n",
    "#    ax1.arrow(high_transunc_particles[i,0],high_transunc_particles[i,1],\n",
    "#              high_transunc_particles_arrows[0,i],high_transunc_particles_arrows[1,i],head_width = 0.02)\n",
    "#print(high_transunc_particles_arrows[0,:])\n",
    "#print(high_transunc_particles_arrows[1,:])\n",
    "\n",
    "ax1.plot(high_transunc_particles[:,0],high_transunc_particles[:,1],'k.')\n",
    "ax1.set_xlim(-0.5,3.5)\n",
    "ax1.set_ylim(-0.5,3.5)\n",
    "ax1.set_xlabel('x')\n",
    "ax1.set_ylabel('y')\n",
    "ax2.plot(high_rotunc_particles[:,0],high_rotunc_particles[:,1],'k.')\n",
    "ax2.set_xlim(-0.5,3.5)\n",
    "ax2.set_ylim(-0.5,3.5)\n",
    "ax2.set_xlabel('x')\n",
    "ax2.set_ylabel('y')\n",
    "ax3.plot(high_transunc_particles_2[:,0],high_transunc_particles_2[:,1],'k.')\n",
    "ax3.set_xlim(-0.5,3.5)\n",
    "ax3.set_ylim(-0.5,3.5)\n",
    "ax3.set_xlabel('x')\n",
    "ax3.set_ylabel('y')\n",
    "ax4.plot(high_rotunc_particles_2[:,0],high_rotunc_particles_2[:,1],'k.')\n",
    "ax4.set_xlim(-0.5,3.5)\n",
    "ax4.set_ylim(-0.5,3.5)\n",
    "ax4.set_xlabel('x')\n",
    "ax4.set_ylabel('y')\n",
    "#print(high_transunc_particles)\n",
    "ax1.set_title(\"Velocity Motion Model\\n high translatory uncertainty\")\n",
    "ax2.set_title(\"Velocity Motion Model\\n high rotatory uncertainty\")\n",
    "ax3.set_title(\"Odometry Motion Model\\n high translatory uncertainty\")\n",
    "ax4.set_title(\"Odometry Motion Model\\n high rotatory uncertainty\")\n",
    "\n",
    "#draw arc for velocity motion model\n",
    "#center of circle and radius as well as angles\n",
    "x_center = particles[0][0] - control[0]/control[1] * math.sin(particles[0][2])\n",
    "y_center = particles[0][1] + control[0]/control[1] * math.cos(particles[0][2])\n",
    "radius = abs(control[0]/control[1])\n",
    "circle_startangle = -np.pi/2 * 180/np.pi\n",
    "circle_endangle = (-np.pi/2 + control[1])* 180/np.pi\n",
    "#draw arc\n",
    "arc1 = patches.Arc([x_center,y_center],width =2*radius,height=2*radius,theta1=circle_startangle,theta2=circle_endangle,linewidth=0.5)\n",
    "ax1.add_patch(arc1)\n",
    "arc2 = patches.Arc([x_center,y_center],width =2*radius,height=2*radius,theta1=circle_startangle,theta2=circle_endangle,linewidth=0.5)\n",
    "ax2.add_patch(arc2)\n",
    "#draw line for odometry motion model\n",
    "ax3.plot([0,ideal_new_pos[0]],[0,ideal_new_pos[1]],\"k\",linewidth=0.5)\n",
    "ax4.plot([0,ideal_new_pos[0]],[0,ideal_new_pos[1]],\"k\",linewidth=0.5)\n",
    "\n",
    "#draw arrow for start and end position\n",
    "arrowlenght = 0.6\n",
    "startarrow = [0,0,arrowlenght,0]\n",
    "endarrow = [ideal_new_pos[0],ideal_new_pos[1],arrowlenght*math.cos(ideal_new_pos[2]),arrowlenght*math.sin(ideal_new_pos[2])]\n",
    "ax1.arrow(startarrow[0]-startarrow[2]/2,startarrow[1]-startarrow[3]/2,startarrow[2],startarrow[3],color = 'blue',head_width=0.1,alpha=0.5)\n",
    "ax1.arrow(endarrow[0]-endarrow[2]/2,endarrow[1]-endarrow[3]/2,endarrow[2],endarrow[3],color='blue',head_width=0.1,alpha=0.5)\n",
    "ax2.arrow(startarrow[0]-startarrow[2]/2,startarrow[1]-startarrow[3]/2,startarrow[2],startarrow[3],color = 'blue',head_width=0.1,alpha=0.5)\n",
    "ax2.arrow(endarrow[0]-endarrow[2]/2,endarrow[1]-endarrow[3]/2,endarrow[2],endarrow[3],color='blue',head_width=0.1,alpha=0.5)\n",
    "ax3.arrow(startarrow[0]-startarrow[2]/2,startarrow[1]-startarrow[3]/2,startarrow[2],startarrow[3],color = 'blue',head_width=0.1,alpha=0.5)\n",
    "ax3.arrow(endarrow[0]-endarrow[2]/2,endarrow[1]-endarrow[3]/2,endarrow[2],endarrow[3],color='blue',head_width=0.1,alpha=0.5)\n",
    "ax4.arrow(startarrow[0]-startarrow[2]/2,startarrow[1]-startarrow[3]/2,startarrow[2],startarrow[3],color = 'blue',head_width=0.1,alpha=0.5)\n",
    "ax4.arrow(endarrow[0]-endarrow[2]/2,endarrow[1]-endarrow[3]/2,endarrow[2],endarrow[3],color='blue',head_width=0.1,alpha=0.5)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figures for the measurement models\n",
    "basic functions needed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#gaussian function with mean,std input\n",
    "def gauss(x,mu,sigma):\n",
    "    return 1./(sigma*math.sqrt(2*math.pi))*np.exp(-(x-mu)**2/(2*sigma**2))\n",
    "#the resampling algorithm of the particle filter:\n",
    "#weights need to be normalized!\n",
    "def resample(particles,weights):\n",
    "    new_particles = []\n",
    "    \n",
    "    for i in range(len(particles)):\n",
    "        randomnumber = random.random()\n",
    "        j=0\n",
    "        accumulated_weight = weights[0]\n",
    "        while accumulated_weight < randomnumber:\n",
    "            j +=1\n",
    "            accumulated_weight += weights[j]\n",
    "        new_particles.append(particles[j])\n",
    "        \n",
    "    return new_particles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "feature based:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#algorithm return probability of measurement given a feature the robots pose and the used variances (distance and angle)\n",
    "def feature_based_measurement_prob(measurement,featurepos,robotpos,variances):\n",
    "    #calculate the expected measurement values\n",
    "    expected_distance = math.sqrt((featurepos[0]-robotpos[0])**2+(featurepos[1]-robotpos[1])**2)\n",
    "    expected_angle = math.atan2(featurepos[1]-robotpos[1],featurepos[0]-robotpos[0])-robotpos[2]\n",
    "    #print(math.atan2(featurepos[1]-robotpos[1],featurepos[0]-robotpos[0]))\n",
    "    #print(expected_distance,expected_angle)\n",
    "    #use these values to calculate the (Gaussian) probabilities\n",
    "    p1 = gauss(measurement[0],expected_distance,math.sqrt(variances[0]))\n",
    "    p2 = gauss(measurement[1],expected_angle,math.sqrt(variances[1]))\n",
    "    #return prob product\n",
    "    return p1*p2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_position = [5,5] #x,y pos of the feature as given\n",
    "measurement = [3,np.pi/4] # measurement of feature (distance and direction from robot)\n",
    "variances = [0.3**2,0.5**2] #uncertainty +- 0.3 for distances and +-0.1rad for angles\n",
    "\n",
    "#generating gaussian distributed particles centered at [3,3,0]\n",
    "#(for the used feature this pose corresponds to an expected measurement [2.82,np.pi/4])\n",
    "particles = []\n",
    "parcount = 300\n",
    "weights = []\n",
    "for i in range(parcount):\n",
    "    #generate particle\n",
    "    x = random.gauss(3,1)\n",
    "    y = random.gauss(3,1)\n",
    "    phi = random.gauss(0,np.pi/4)\n",
    "    particles.append([x,y,phi])\n",
    "    \n",
    "    #weigh particle\n",
    "    weights.append(feature_based_measurement_prob(measurement,feature_position,particles[i],variances))\n",
    "\n",
    "weights /= sum(weights)\n",
    "    \n",
    "#resampling\n",
    "new_particles = np.array(resample(particles,weights))\n",
    "particles = np.array(particles)\n",
    "#print(new_particles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlwAAAFDCAYAAADvdNBjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOyde5xVZbnHf88MMzCD4IACxgRCpFCIgmDgoQuYREnapOUlrCyTOnVORy0S0gJLgg6mdc7JSrudkgxvzdFQERNLTUkRECnwkggOCggM12GYy3v+eNca9qy93nXZe6+19tr79/189gf22uvyrLX3eua3nvd5nleUUiCEEEIIIdFRkbQBhBBCCCGlDgUXIYQQQkjEUHARQgghhEQMBRchhBBCSMRQcBFCCCGERAwFFyGEEEJIxFBwFREiMl9E2nPYbpiIKBG5NGPZZSLymQIeI6ftosA63/kiMtTlMyUi14Xc32YR+XnhLMwPEamzzu/UpG0hpFQQkcdE5JEC7esEy9dcVoj95YvlL97vsjz0OYvIr0Xk5cJZlz8icqWInJe0HfnSI2kDSDd+DuDBHLZ7A8CZADJvkssAtAP4TYGOUUwMAzAPwCMAtjg+OxPA1rgNKjB10Of3MoDnE7aFEFL8zIP2939xLP8ygFJotnkltL+/L2lD8oGCq4hQSr0O4PUctmsF8HSUxygGRKQCPr9ZpVSg60AIyUZEBEBPpdThpG0h/ohIjVKqxfS5UurvcdpDvOGQYhHhHLbLGCr8vIgsFpEdIrJLRH4nIv1c1rvUev8YgA8A+KC1XInIr92OYS1bICJrRWS/iGwXkWUi8q4c7LftuEJEbhGR3SKyV0R+LiK9M9arEZH/EpF/iMhBEXldRH4vIvWO/T0mIo+IyMUi8ncArQD+BcBKa5XHM85vmLVN1pCiiJwhIg+ISLN1vDUiconPubxTRO60zqFFRP4qIu8NcA1miMgq61ruE5H1InKFY52zReQvli17ReRuERliX0MAr1qr/jbj/Kb4HZuQsNjDRyLyQRFZDeAwgIutz6pE5NvW50esofdvWqLM3r7e8kdvishhEdkqIndZD0cQkUHW/f+ydR+9KiI/y/Rf1nqbrfW+bP3/oHXPHm8d437rnnpZRC52bGv7iYtEZJNlxyoROSPA+Q8SkV9Y9rdafvBjLutdIyJNll0PQUfZg1xf+7y+KiKvWdfgERF5p2O9/7BsbrZ8zmMicqZjnfki0i4iY0XkcRFpAXCNiNgRrO9m+IvLMq+NYz9vs875DeucXxGR7/mcRx8Rudn6flut63yF1zbWdu+yvru3Mr7//3Gs4+lrRWQzgBMBXJ5xfvP9jl2MMMKVDuzhs08DGAFgMbRj/Lxh/S8DuB06xPxVa9lOj/2fYO2zCcCxAL4I4CkRGaWUejNHe/8C4FMARgH4HoBeAOwcsxrr/TwA2wEMAnAVgCetY2Y+XY8BMB/AdwDsAvBPAF8B8GPLTnvI7Q03QyyntRLAWgBfAvAWgFOhb2BXRORE6IjhqwD+FcA+a9tHRGSSUmqtYbsRAP4AYCmAbwHoBPBu6Gtqr3NexjrfB9DbOr+VInKadR7nA7jXWr7c2pRPqiQqjgdwG4AFADYD2GYt/x2AD1nLnwPwHgDfBtAXwBxrnd8CqIe+f7cBGAxgBgBblB0HYD+Ab0DfeycCmAvgAejh/0ymAxgO4N8sm34E4JfQ/uEuAP8FfR/eLiLPKKVeydh2NIAbLPsOA7gWwAoReadS6i23kxaRYwE8Yb29BvreuwjAH0Tko0qpB6z1vgRgEYD/BvBHAJOh79+gfBjAWOhhsV7Q1/Nhy9cdsdY5EcDPoK9/L2jR+2cRGa+UWp9pNoC7LVvmAdgL4CEAT1nb/9paL/PaZJ5zfwB/BdAT2qduso492WS8iFRB+6Hh1jYvQf8ufioilUqpn3qc+x+hv/crAOyxjjUpY99BfO3HoX8vfwOw0No0laM0UErxVSQv6D+w7Rnvh0GPvz/oWO+/ARxyWe/SjGWPAXjE7xgun1dC34y7AFwVdDuHHascy78MoAPASI9jDoQWKB93nEM7gJMc60+xjvNel30pANdlvH8C2kFUe9i9GcDPM97/EvqGPsZh4z8A3Omxn09Yx+/rsc4rABpdrtsRAF82fZ988RXFC/oPtAJwlmP5+63lDY7l10ILmn7W+wMAvhrieD0ATLD2PS5j+WYAbwKozVh2o7Xe1RnL+lk+4ZqMZY9Z652Ssex4AAcBLHCs90jG+29b9r/dYeMKAH+z/l9h+YI7HevYtl3mc76bARwCcHzGstGWr7vCsE2ldZ3WAfhRxvL51jE/77JNN7/ncc43AGgDcLLPb+LljPefsew93bHebdZ3VmHYz/GWXed5HCuQr4XDR6f1xSHFdPCA4/0LAGpE5PhC7FxEPmKFqPdAO7PDAPoDGJnjLhsd7++GdlzvyTjmpSLyrIjss465HfrpzXnMTUqpl3IxQkRqoZ+ib1dHnySDMB36yeywiPQQkR6WbY8A8BpWXAd9LneIyHnW02SmPScBeAeApfZ+rX2/DmCjz74JiYqDSqlHHcumQwuWBxy/1YehH8js4bpnAcwWka+ISJa/EM1/iMgLInIQ+o/9M9bHzvWfVEodynj/ovXvCnuBUmoPgB0Ahji2fVEp9ULGem8B+DMyoikuTAfwOIA3Hee4HMB4EellHaceOiqdyZ0e+3XyuMqIsimlNkDf75mRnkki8pCI7IT2IW3QkXg3H3x/iGM7ORvAE0qpF33XPMp06EjY8y7XaRCAdxq22wUtlBaKrpp3G1XI1demEgqudLDH8b7V+rdXvju28hzugw77fhbARGhn+noe+9/heL8T+glpsHXMj0MPRawFcGHGMTtcjpnLkKZNP+jfeFPI7QZCD1e2OV7/Bj1E4oolDM+BHjK9G8AOK19jTMZ+AT1U49z3GK99ExIh212WDYQe7m5F99/p36zP7d/qRdBVz/MAbLRydL6YsZ8rAfwQ+o/qx6AfumZYnznvdaefO+Kx3Lmt0+cA+rwGuyy3GQg93Oe8FxdD+41+AN5m2L/bNTPhaZslRB4GUAs9nHYmtD98Btnn2amU8koP8eM45OYPRyH7Ot2Vsc8slA5NTYNO+/gRgM0issHy/5n7Du1r0wpzuEgDdETrk0qpdqCrUimfH/tAx/sB0A7MzrP6BHTk6gv2CiJyAnQo2Uk+Jc17oIVevd+KDnYB+BOAm8Pao5RaAZ070gvAB6GHHpYBGGrtFwC+huzybUDnLxASN26/6V3Qv8cPGrb5JwAopbYDmAVglvVgcRV0bs/LSqk/Qd/ry5VSds4XRMQr6pQrTp8D6OiLa26nxS7oFjLfMHz+FoBqw/4HFcA2u6L6QwD6ALhQZeTMikhf6OHITPJt8fAWcvOHm3A0B9fJRtOGSqmXAVwiuohiAnRu610i8m4rypazr00jFFylSyv0E6oftdAh7Mwf9yegozS50oCjyY32/hSOPh3X4ujTq82nQ+w/UIRPKXVIRP4K4FIRWRhiWHE5gNMBPB9yKDLz2IcBLLMS6X9kOc9NAF6DzmW7yWPzgkUwCcmR5dCJ5FVKqaeCbKCUWi8iVwL4HHSe0p+g7/VdjlXD3OtBOVlERlvDdbDSLT4AHVkxsRw6uvKKNVSZhYhshS4G+DiAOzI+ujCEbe8TkePtYUURGQ0dMfqh9Xmt9W+Xr7FGHkYieITfLernxiMA5ojIySGGFZdDRyd3KaVe9VvZDaVUJ4C/ici1AD4Kff4vIrivbUUJ+EMKrtLlHwCusMK3WwG8pZTa7LLecuiw/89F5HYApwD4OrKdZBjqReR30E1XR0JXKf5eKWU/CS0H8BMRWQSdn/Ev0MOZQTvZvwg9/Hi5iByCdjamG/Ya6CrFP4vIj6DD+6cAqFFKfd+w/29Dh/NXisgt0MOrxwMYD6BDKfUtt42soZT3QefcNUEPGfw7dBHBPmudfwdwr4jUQFciNlvrTQXwkFLqHujhht3QT4abALRARwT3B7w+hOSFUmqliNwJ4D4RuRHAagBV0FXSH4MeijsG+v69HUejHJ+FHhJ6zHq/HMDXRWQ2gDXQQ+7TIzD5TejqQrtK8Tpof+IWObG5Gboa8HERuRm6oOVY6NypeqXUl5RSnSKyAMCPLf+xDLqi75MhbNsN4CERuQFHqxS34GhF4Z+g/dlvrGMMhR6iDdPA+R8AzhORR6Ejk68qpdx8+M3QgneliHwX+iHw7QDep5SaZdj37QAut7a5EcAG6If5UQAmKaU+4baR6Jkyboau6HwF+ty/Cl1ZucpaLaiv/QeAKSLyYego3TallF1Nmx6Sztrn6+gL5irFSx3rXWYtf7tpPeg/4g9DV+EoAL92O4a17CroyEsLdMnwRGgH+muTbQb7bTuugC5R3gN98/8S2VUoi6DD/QehnfLJ0I5yfsZ6j8Gl0tL67MuWzR3WMYdZy7OqdazzeRi6PP0AtOO/KOPzzXBUwECXL/8vtCM/Au0I/g/AhzzO/0zohNYm6Cey1wH8AsAJjvXeZ53zXuuav2ytNzJjnQughaUdfZyS9O+Tr9J7wVGR5visEvrha4P1e94NHaWeB53Y3BPArdB/DA9CPzz8JfMegY7e3Ar9R3IfdN7PWDgq/Az34GXI8HOmdW0/Ad2G5iXL1r8BmOjYLsufQKdO/A+0ADpi+aSHAVzsWG8udKTrkPX5mc5zMFzDzdCze1wFLaAOQwuskx3rfcq63w9D57aeA93u4bGMdebD4IOhq0rXWefQZZfhnAdb3/sO63gvA7jB6zdhfY8LrHWPWNs+DuArHuc+EPqh+2VoP/cWdL7fGY71fH0tdHudv1r2KmT8nUjTS6yTISRv5GjTzk8rpW5P1hpCSDkgutFzu1Lq7KRtcWI17XxEZeSrkvKFVYqEEEIIIRETi+CySkEPZLxaRLfnPz2O4xNCSD7QhxFC8iWRIUUrCbFBKTU69oMTQkie0IcRQsIS+5Ci1Un289BJ1YQQkirowwghuZBEDlcDdOntbxI4NiGE5At9GCEkNLEPKYrICgCvK6U+57HOLOjuxejdu/f4UaNGxWUeMdHRATQ1AUOGACJJW0NKmNWrV7+llBqQtB0m/HwY/VeR8sYbQF0dUJNPT2dC/DH5sFgFl9V1+yUAZyqlVvmtDwATJkxQzz77bLSGEW86OoCrrwa+/W3guJKb3ooUGSKyWik1IWk73Ajrw+i/ioQ77gCOOQY499ykLSFlgMmHxT2k+EUA64KKLVIkLFoEzJpFsUUIfVj6eOYZ4M03KbZI4sQmuESkGrpz8E/jOiYpAL/7HXDaacBoFmOR8oY+LIVs2wYsXQpceWXSlhASa4TrfOgJkZfEeEySD3/7G7BjB/DRjyZtCSHFAH1Ymjh8GFiwAPjOd5h3SoqC2ASXUur3Sqk+SqkDcR2T5MG2bcBddwH/8R9JW0JIUUAfliKU0jmnc+cCtbVJW0MIAKBH0gaQIqSlRT8Z3ngjnwwJIenjf/4H+PjHgbe/PWlLCOmCcymS7igFzJ+vnwxZPk0ISRvLl+v2D2eembQlhHSDgot05667gLPP5pMhISR97NoFPPoo8OlPJ20JIVlQcJGjbN8OPP88MG1a0pYQQkh4vv99YM6cpK0gxBUKLqJRCvjP/wS+8Y2kLSGEkPA0NgIf+ADQr1/SlhDiCgUX0dx5J/CRjwB9+yZtCSGEhGPXLuDpp4EZM5K2hBAjFFxEDyWuX69ztwghJG1wKJGkAAquckcp7ayuuSZpSwghJDx/+AMwZYquTCSkiKHgKneWLtVh+D59kraEEELC8dZbekaMc85J2hJCfKHgKmfefBPYsAH44AeTtoQQQsLD6DxJERRc5cyNN9JZEULSyf/9HzB1KocSSWqg4CpXnngCOOMM4JhjkraEEELC0dYG/PnPHEokqYKCqxxRCrj7buCTn0zaEkIICc+vfw1cdlnSVhASCgqucuSPf9SJ8hX8+gkhKePQIeCVV4BTT03aEkJCwb+45UZHB/CnP3H6HkJIOvnZz4AvfjFpKwgJDQVXubFkCTBzZtJWEEJIeHbvBvbuBYYPT9oSQkJDwVVOHD6sJ6c+44ykLSGEkPDccgvw5S8nbQUhOUHBVU78/OfArFlJW0EIIeF5/XWguhoYODBpSwjJCQqucmHfPt3o9OSTk7aEEELC85OfAP/6r0lbQUjOUHCVCwzFE0LSysaNQH09pyAjqYaCqxzYt0/nbw0enLQlhBASnt/8Brj88qStICQvKLjKgd/+FvjMZ5K2ghBCwrN9O9CvH9CzZ9KWEJIXFFylTmcnsHkz8I53JG0JIYSEh13lSYlAwVXqPPAA5xsjhKST1ladEjFgQNKWEJI3FFylzmOPAVOmJG0FIYSE5667gAsvTNoKQgoCBVcp89JLwMiRgEjSlhBCSHiefx447bSkrSCkIMQquETkbBF5WkQOiMhbInJLnMcvO37/e+Dii5O2gpCSgT4sRl54ARg9OmkrCCkYPeI6kIhMAXA3gC8AuB+AAHh3XMcvO44c0a0g2LeGkIJAHxYzd90FXHNN0lYQUjBiE1wAFgL4qVLq7oxlz8V4/PLi/vuB885L2gpCSgn6sLhoadEV1rW1SVtCSMGIZUhRRHoDeA+AwyLynBWKf0xEJsRx/LJk1SrgPe9J2gpCSgL6sJi5917ggguStoKQghJXDlc/61hXALgMwGAADwN4QETqnCuLyCwReVZEnt25c2dMJpYQu3YBgwYxWZ6QwhHYh9F/FYAXXgDGjk3aCkIKSlyCa7/176+UUs8rpY5Ah+erAPyLc2Wl1K1KqQlKqQkD2H8lPA88AMyYkbQVhJQSgX0Y/VeeHDrEoURSksQiuJRSewFsBqDcPo7DhrLiH//Q7SAIIQWBPixGHnkEOPvspK0gpODE2RbiFgCfE5F3i0gPALMBHAbw1xhtKH1aW4Hqag4nElJ46MPigPmnpESJs0rxRgB9ADwKoBeANQA+Yj05kkLx5z+zszwh0UAfFjWdnfrfyspk7SAkAmKLcCnNt5VSJyil6pRSU5VSa+M6ftnwxBPA5MlJW0FIyUEfFgPPPAOccUbSVhASCZzap5RQCmhvB6qqkraEEELCs2IFMG1a0lYQEgkUXKUE5x0jhKSZQ4eA3r2TtoKQSKDgKiUefBD48IeTtoIQQsLz6qvAiScmbQUhkUHBVUrs3Qsce2zSVhBCSHj++Efgox9N2gpCIoOCq1Robgb690/aCkIIyY3t24H6+qStICQyKLhKheeeA8aPT9oKQgjJDfYOJCUOBVep8NxzwOmnJ20FIYSE58039fyvhJQwFFylwt69QF3WPOCEEFL8rF7NCD0peSi4CCGEJMu6dWxpQ0oeCq5SYO9eoG/fpK0ghJDcaGkBamuTtoKQSKHgKgXWrGH+FiEkvTBhnpQBFFylABPmCSFpZccOYMCApK0gJHIouEqB5magX7+krSCEkPAwYZ6UCRRchBBCkmPdOmDs2KStICRyKLjSzoEDnOyVEJJeDh1iwjwpCyi40s7WrcDQoUlbQQghucGEeVImUHClnTfeAAYPTtoKQgghhHjQI2kDSJ688QZwxhkAgMY1TVi8fBO2NbdgcF0NZk8fiYZxnAyWEFKkdHQAFUef+2fe9hSefGV31/vJI/pjyRVnJmEZIQWHEa6088YbwNvehsY1TZh773o0NbdAAWhqbsHce9ejcU1T0hYSQog7b73V1RLCKbYA4MlXdmPmbU8lYRkhBYeCK+3s3w/06YPFyzehpa2j20ctbR1YvHxTQoYRQogP27YBb3sbAGSJLRvTckLSBocUS4RtzS2hlhcjHBIlpMxgDiopIxjhSjtWhc/guhrXj03Liw0OiRJShlgpEYSUAxRcaUcpAMDs6SNRU1XZ7aOaqkrMnj4yCatCwyFRQsqQN98ETjgBgE6Qd8O0nJC0QcFVIjSMq8fC88egvq4GAqC+rgYLzx+TmiG5UhgSJYSE5MgRoLoaALDkijOzxBWrFEkpwRyuNHPgAHDMMV1vG8bVp0ZgORlcV4MmF3GVliFRQkgOOJqeUlyRUoYRrjTz1lvA8ccnbUVBSPuQKCEkB6yUCELKAUa40kxbG1BVlbQVBcGOzJmqFFnBSEgJwml9SBkRi+ASkV8DmAmgNWPxN5RSt8Rx/JKlvT204Cpm4WIaErUrGO2keruC0d6GkKihDysOTp33EPa1Hi2u6duzEs9f/+EELSIkOHEOKf6vUuqYjBcdVb60twM9gmtmr9YLjWuaMHnRoxg+ZxkmL3q0qNoxlEwF45IlwLBheiqTYcP0e5Im6MMSxCm2AGBfawdOnfdQQhaVGfRfecMhxTQTUnCZhMv8+zagtb2zaCNIJVHBuGQJMGsWcOiQfv/aa/o9AMycmZxdhKQEp9jyW04KCP1XQYgzwnWBiOwWkRdFZLGIHOO/CfGkvR2orPRfz8KtChAAmlvaijqClPamrgCAa6896qxsDh3Sy0laoA8rNAVKmr+ucX1B9kMM0H8VhLgE138DGAXgeAAfB/ABALeZVhaRWSLyrIg8u3PnzphMTCGVlUBnZ/DVQyaoFksEKYoKxtiHULdsCbecFBuBfRj9VwgKlDR/x6qtBdkPMUD/VRBiEVxKqdVKqe1KqU6l1AYAVwH4hIj0NKx/q1JqglJqwgBrJnniQo8eOsoVkI6QT5PFEkEqdFPXRKYRGjo03HJSVITxYfRf0dC3pzmaH9a3kZDQfxWEpHK47LAMa4LzIaTgqjc0F62rqeqWwwVE1wMr1yrJQjZ19UrCjyxnbcGC7jkQAFBbq5eTNEIfFjPPX/9hDJuzzPWzsNF7EhL6r4IQS4RLRC4WkTrr/ycB+AGA+5RSh+M4fsnSo4fuxRUQ09Dc/PNGxzItULFMUJ1IEv7MmcCttwInnqiHUU48Ub9nwmkqoA+LiJCRqUsnuUdULpk4pBDWEBP0XwUhrgjXlwDcYoXfdwD4A4D5MR27dAkpuPyai0ZdkZhIZMmFQk4jFCpiN3MmHVR6oQ8rAm5oGANA52x1KIVKEVwycUjXchIh9F95E4vgUkpNieM4xUJszUXr6oB9+0JtkuR8i8XS3mH29JHdGqkCuQ2hsiFr+VBuPiw2chgKvKFhDAUWSSWcS7HAxDps1q8fsHt34fcbEcXS3qFQSfgl05CVEEJI5LDxaYGJctjMNXKW1x7jpVCRpUJgR/rsa3rV0rVYvHxTqGhksUTsCEkD0256DC/tONj1/qSBvbFCBOjoCNVPkJC0QsFVYIL+EQ477GgcvurslxrR5ZdDFjf5DgkWMheMkDRxXeP6UHlUTrEFAC/tOIhplWOw4q23gEGDojaZkMSh4CowQf4I5/KH3hg5w5DUCC4g2RwyJ/lGI4spYkdIXFzXuB63P3204WWHUl3vTaLLKba6lndUA9u2UXCRsoA5XAUmSFf0XHJ/jJEzVOdhbfESRyf4fIcEC92QlZA0YOrqnlu3dwHeeCM/gwhJCYxwFZggw2a5/KE3Rs6kDWhtBXq6Nu1PJXFV/xViSLCYInaExIGpq3vO3d4puEiZQMEVAX5/hHP5Q28cvjquBXjzTTTu7lE0uVGZ5NIiI65+XRwSJCQ8lSKu4sqr2/tJA3u7DiuedHwt8Oabrh3kNy+akZ+hhBQZHFJMgFwmYzYOX43qj8ZVrxZFB3cnubbIiKv6j0OChITH1NXdq9v7iqun4KSBvbstO2lgb6z4+lQM23+q6zamaXwISSuMcMWAW5Rn4fljQkd+XCNn697C4r9uRUtb96fLJDq4O8k1UhVn9R+HBAkJR67d3ldcPcXwCedBJOUBBVfEmPKRFp4/Bk/OOSv/A9TXY9sR92TVpPtB5Rqp4lAfIcUNu70TEh4OKUZM5N3Ijz8eg3HE9aOk+0GZjl8h4ll9yKE+QgghpQYjXBETRz7S7IotmFs5qugiQm6RKuBoNZNX9SGH+ggpFxQ4rEjKAUa4IiaO+QMbag9g4bkjiyIilNk/a/HyTbhgfH2XXW5VTLlE++Lo0UUIiYfNX3gXgM7s5axSJCUGI1wRE0s+0qmnogE70VCInLA8cMtXu2d1U5f4G26oOgoz7VHjmibMvmsd2jqPRslm37UOQGF7dBFCYmLECGyu+l/gu99N2hJCIoWCK2JimT9w/HjgvvuA97yncPv0wBZETc0tXT156utqcOhIu2dVYiGmPZp/34YusWXT1qkw/74NFFyEpBGP/l2ElBIUXDEQeT7S4MF6PrIYcAqizHwsE3YEK0i0z6+VRHNLm+sxTMsJISmgshJobwd68E8SKV2Yw1UKxPiE6CaI/LAjWEGqD+NqekoIKSJGjQI2bkzaCkIihY8TpULPnrHMqRhW+DgjWPlOe9Svtgp7DmVHs/rVVoWyixBSRIwfDzz1FHDKKUlbQkhkMMJVKowZA6xfH/lh/Kor62qqcqqWtCsPm5pbsgrEM0XbvHNHo6qy+xpVlYJ5544OcxqEkGJixAjglVeStoKQSGGEq1QYPx5YtgyYMCHSw5h6awFaGM0/b3TofDVnXphzWtxeVUefC4IUIeQyYTYhJEEqKoDO7NYQhJQSFFylQn098NprnqsUQohkCh5nlaLX/ryO7ZcXtudQW7dKRa9hSb8qR0JIkVJTA+zfD/TpY1yFD1MkzVBwlQoiQFUVcOQIUF2d9XEhhUjYqku/YwfJCws6GXeuE2YTQhLmQx8CVqwAzj/f9WM+TJG0wxyuIiZ0R/X3vQ94/HHXjyKf09EDv2MH7bofRJixypGQlHL66cCzzxo/TtKHEVIIylZwFfv0MPbTXFNzCxSOPs152vn+9wN/+YvrR6Y+WU3NLZGfu58Imj19ZKAfYhBhFsdUSsX+2yEklVRU6FeHe3qByY949QAkpJgoS8GVk5iJmZye5qqrgbY2QHVPO29c0+Q5NWzU5+4ngp59bbfLTGrdCTod0uzpI1FTVZnTtkFIw2+HkNRy5pnA00+7fuT10HRdY/QV2oTkS1kKrjSEpnMeGhs9Gvj737stWrx8U1blXyYtbR24culajL3+4VDCIWikx00ECYCpowYAAO5YtdXzOGHaSwRprpoPafjtEJJazt0MhgkAACAASURBVDoLeOQR14+8Hppuf3oLhs1ZhmFzluGdc93nbCUkaWJPmheRCgBPADgTwBCl1Otx25CGPJ8g8w668pGPAL/4BTB6dLc5D4PQ3NIWeCLoMAmsDePq8exru7Hk6S1dwk8BuGd1Eyac2L9reiA36utq8GTISbmjnEopDb8dEh3F4L9Kmpoa4PBhAO4ViUFoV8A75y7DywtnRGkpIaFJokrxKgCHEjhuFzmLmRhx63eVGRUy0r8/sGdPliAKSlun6lbRZyrDDlsNuHLjzqwom72+3VrCDVvIFEs5eBp+OyRSEvdfaeK6xvW4Y9VWdCiFShFcMnEIbmgY473RSSeh8cHVmPvXnVkPdCJZGROutKs8jk9IRMQquETkZABfBnABgDVxHjuTIJMoJ41fVMiz0ecx78DiB/4eWmzZZIocUxQrbKTHa/nMSUNx+9NbXD8fXFdTVOXgafjtkGgoFv+VFq5rXN/tvu5Qquu9U/SMuvYBHO6wPd0g4MU3AEfmaUtbB2qrKnCoLViD1DDHJyQOYsvhskLxvwQwG0BzXMd1I+o8n0LhFRWycU3iPjgYTftbjfut8Jnr2o7WeEWxwlYDei2/oWEMJo/on/WZLWSKKW8qLb8dUliKyX+lBVNupnN5d7HlTVCxFeb4hMRFnBGu/wDwplLqXhEZ5rWiiMwCMAsAhg4dGokxUeb5FIog0SJXMdKhUKkUOsRdWXX6+DZ72NLr+DdfNDZUpMcvMrTkijONw4ZXLV1rtCNX8hmiTMNvhxScovJfacCUJuBc7i62fJ4KfeghQHvA4xMSF7EILhF5J4CvAQg00Z9S6lYAtwLAhAkTyvbuCJIvZBIdHSKo6SFoaQ9/+VZu3Ol7/CBzGmYSZH2TkDHZcWxNVehzA9ixmoSD/is3TLmZlYYHwWwUchFePQR4eeEMjJj7QJ7HJ6SwxBXhei+AAQBeEP1jt4cynxeR65RSt0R58GJJuA5LkHwhkxip79sTs9texuJjRmJbcwuOralCc0tboONmNiT1Or7fnIZu19zrupu2mT19JK5eujarV9fBI+1oXNMUep+c/oeEJFH/lVYumTjENTfzkolDgu1AAfX9emFb82HUVlfi4BFzTurmRdkViXkfn5ACE1cO150ARgAYa73OsZZ/CMBvojxwmhtVBskXcutxBQAH2xTQcghPfv39eHXRDPTuGVxb2xG0XPOVcrnmftu4ZW60dSjPPC7TPk1tMtjagRhIzH+lmRsaxuDSSUO7IkqVIrh00tCshPVelYaIkwim9m7Fq4tm4LBP7tbM257K+fiExIWoBMazrRyIVxGgj82ECRPUsx7za/kxedGj7hGgHPo7hSWuyFrjmiZcf/8G7DnUPYJVUwEsfKdCw+c/imFzgjUDrKmqzDsJPJdr7rUNYJ6+QwC86vJ067VP01BHHL8J4o+IrFZKBRq+S4I4/Ve5YE6cV7h00onGKuZMJo/ojyVXnFl44wgJicmHJdGHC0qpzcg3KzIgSTWqjDNPyO6L5RRcLZ3A4k2tgDW1j5+0rq+rwdRRA7B4+SZctXRtziIxlznPcv2evPK4jPltSqGmqpKtHUhOxOm/yoWNC84x5FwJ7li1xbNXn82Tr+wGULjeW2lNRSHFS8lP7RPHZMZu5NPKIJfJkY2CRXpi8R9f8BVbAj08ec/qptBDgU5bva7tMMM5eX1PXvuz87hM27phD42ytQMhxYNXVWPQvCu795a9L7v3Vti5FtOcikKKl5IXXFFPZmwi14hNrje6UbAc2wvbDvonyw+uqwktEk22Th01wPPx3+2c3LaxvydTnhrgncfl9d03jKvHk3POwquLZuDJOWdRbBGSMKbqwUoo3PCRk7vlY5koVO+tYur9R0qHkhdcSTWqzDWyZrrR59+3wTPqZRQXH34XBssRz2PaIsQ05GdaPv++Da623rFqq29ELdN5Na5pwj2rm7ptIwAuGF/fVdm48HzzkIBJxLJJKSHpwRTFuuSdfYDf/x43NIzBKwvPcW2SDOgcrqC9vwDvkQTOmUqiIJEcrrhJolFl2Clg/Caabm5p62rr4JYP5tnn6q0RmPvoFrSoo/razumqz1jva3euC9y3pnFNk7HNRNDGgk3NLRg+ZxkqXPIzFI72A7PPz3R9vEQsm5QSkg7sPCvX/Kuvf11PoiiCJVeciZm3PdWVswUcTZgP2nvLL8eWc6aSKCgLwZUEYRqD5jLRtFvfKJO4aJh2GvDICiyuG4tte1uNtoR5OixUaF15HNf5NOklYpngSkj6uaFhjHuC+7RpwIMPAufojhymasSgvbf8evFNHTWg2zy2AAtrSP5QcEVI0OiK280fhDDh7YZZH0fD/fcDc6/M+swWKybqXZ7q4gitO58mTSIWQCwVoRR1hJiJ9P740IeAq64CPvxhoCI7E2baTY/hpR0Hs5abqhS9hgz9UhwIyRUKriIgV/ESKrw9YgSwaxfQ3AzU1XUtDhJds+dWdB7bq81DvmQ+Tfo58smLHo28czynAyLETOT3hwhw8cXAHXcAM2d2+8gktgAdPV+5cWfWjBReQ4ZuD8DOFAdCcqHkk+aB3NosxEkueQFu4W3f8/zKV4Af/7jboiDRtaXPbA2UpF9VkV9rokqRrOT2IFWbcSS4smqJEDOFuD8mLliBYXOWdb0mLljRfYVJk4A1a4Aj3YuATGLLxs1neFUwM2GeREXJC6409FMJmxdQKZJVbRfoPE84QYfjt23rWhTEibi1XnCrAFz8ydNQl+OE0jVVlfjBhadltWkI4sjj6LVGJ0yImXzvj4kLVmD7/u5Cavv+I9mi6/LLgV/8IrR9Tp/hVcGcVO9GUvoEFlwiksoZP9MQmWgYV28UKm69qX5w4WlZYfrA5/mVrwC3HJ1rN6gTcXOczl5WueLVriGII4+j1xqdcLpJq/9KC/neH06xZVz+rncBW7cC+/eHsg/I9iVuvfhm3vaU61AjE+ZJIQgT4fqniDwgIg0i4t6FsghJS2Ri/nmjs0SD3brBLmnOV5gAAPr21ZGuF18EYJ782omf47QjbKZWESbs+QtNeR5BHHkc/baSaqBLCkYq/VdaiPX++Nd/BX7yk663Jw3sHWgzPx/mbDVhU1UB9u8jBSGM4BoFYC2A/wHwuogsFJER0ZhVOEw3mQKKKp8rUzQA6Db3oT33n1fVT6gnzC98Afj5z12P60aF+A97mnLB6mqqPPftJ3yDOvKoO8eziWrqSaX/Sgux3h9DhgCHDwM7dgAAVlw9xVd0BRF/bmILANo6WRhDCkNgwaWUekUp9U0AQwHMgnZgG0TkERG5SESKsuLRK4Ljl88Vd7K9LRrq62qyOrX7DYOGesLs1Qs45RTg2WezjuuGAnDV0rWe18AknPa2tHnu2++ps2FcPS4YX98V5asUSaw8275ON180FoD/NSHFQ1r9V5rI56FnUJ/qUMvxb//WrQBoxdVTPPff0taBK5euxbA5y3DqvIcC20VIIQmdNK+U6gSwAsAfAGwC8B4ACwG8KiLTC2te/vhFcExCJslk+1yGQUM/Yc6cCfzmN0Bnp+/+lYLvNfCLsLkJQoF7y4lM7J44mZPR3rO6KTGRk4YiDGImbf6rXFh17bQscTWoTzVWXTvNfYP+/XV7m03hc3H3tXbkLLqKveKdFDehnupEZDyALwC4GMArAG4BsEQpdUBELgXwSwBFF3u1G5AOn7PMdY4/N6Hh14k4DGEbApp6xNTVVmHyokeN+wnSaLWbLcd9BLNvXoKGr326a/97DnnnYJmugd9URg3j6vHsa7u7dW9WAO5Z3YQJJ/Y32p3P9xBFI8ZC/i78YKPVwpJW/1Uq+P2ejeLKwdGKxpOBX76EQX23BN7WZl9rdvrD5BH9XYcV7bkb2YuP5EuYKsW1AFZCi7RpSqkJSqmfKaUOAIBS6nYAufUEiIkweU6FSrbPJSLi2uOqUnDgcHtekZUsW1o6MXdnHRofXoPGNU04cLg90H7cxGCQCNvKjTtDD5WGnVDbJqpIVFxFGIykFZZS8F9pplC/56z2ESLYvr8Vw+csy9vGJVecmTUxtj1HI5COindS3ISJcP0U+mnQWI+rlBqYv0nREWZC6UJNXppLRMRtCpuDre1ZFYBhIyuutqACix95Bei/G22dwSaddpvM2rbby5YwYqVxTRPm37chtA02put+/f0b8noajWtS2zgjaWVC6v1XminU79m9fYS4jlzkgmmORiA9Fe+keAksuJRSP43SkDgIM6F0GHHmRdCb1C3cntnbapjhCS7M9DpGWypqgOYWZHf9cqdDKQyfswzH1lRBBGg+1JbXUKlTrDSuacLsu9Z5CkDThNc2pnPdc6gta5qPMBTqd+EHnXthKQX/lWZyjVQHRik9/U8InD71pIG9PZPv43rYIqVLyXeadxK0kqZQZc5BhjH9wu1eYXe/SE8QWxTC/xAUgOaWNuw51JbXUKmbWFm8fJNvtM2r1QTg7QTzGQKIq/ydjVZJKWHyU2H8lycF2M9LOw5i2k2PGT9nLz6SLyyFdlDoROUgERG/3AA7MdONzEiPn+1utnTtB5L1lFhVIYDoqX2CkMtQqdv19YviBHFys6ePxJVL17p+Ztp/0O8+SHFCvsQVSSMkDkwRab9ItZNBfaqNXekLgde8jGFGSAhxg4IrgyiqUILcpF7DR36TS9uRniC2Z9riGsoXQSUUOiFddnqub7DZxiRg/K6lKXRvE6QPV8O4esy/b4Nr53u3KFGxVSDRuZNSot5wT/tFqp2sunaa67yLbpw0sDf2tbTlJdCua1yPO1ZtRYdSqBTBJROH5DSNGSuOCUDB1Y2oEpX9RIZXboBXtCcz4hHUdr8WGR3onsnlt76bzUB+Amb29JGeOVx2KwnAW5DMP2904ChRMSapxxFJIyQOChmxtVtAmPJadcKD4KUdBzGoTzU2L5rhsa6Z6xrX4/ant3S971Cq6/0NDWMC76fYHuZIcpRdDpcXSSUqe+UGmHJ2KkW65Q6ZIkK27c6GfccaJsu2K36cOVlBcoeCCMArA3RnbxhXj8WfPM04obddbehXZh4m34pJ6oRER7xTYx19ZNy+/wgmLlgROFcsc4qgO1ZtdV3n9qe3YOZtTwW2hu0kiA0jXBkkVYXiN3zk9mSY6awa1zR1m3sxEwVg7PUP4+CR9q5crKbmFlRVCqoqxDM5PTPCY8r/6l1diUNHOrJs9hIqQZ7w/CJrbg1avSJ6fuT63ZuGCtyWA+7fMYcbSDkQJmLrnEg6sx9WWLbvP4JLJw3tFq1yw1ml6JVf9uQruzHztqcC2WTyhU3NLRgx94Fuw5VhImckfVBwZZBkorLJGQXJ5Vm8fJPncJ9bHpMtvmqrKtDS1mnc3nYWYXOK/PKwMp/wvPbptx8nuZaZ5/Ldm4YKnn1tN+5Z3dRt+ey71nUrQPBal8MNpJxxii3gqMD55IShXf7ClBfmhi1knPlYXgKnUiSQ6PIThl4+LHPKslyGK0m6EBWySiRuJkyYoJ61JlmOgzRGG4LmV5m4dNJQrNy409UpVALoBHyvhfO6TR01oJuQCIpb9M5NCB1u74DbT7dSBK8sPCfUMU3n4PfdT170qPs183HUmVQI4BZkrK+rySk5txQQkdVKqQlJ21EI4vZfpUCYfKtMf+G13eZFMzz345YcD8A3KuaGU3S5+TAT+fgvUjyYfFhsES4RWQDgUwCOA3AYwF8AXK2UCv+LjpA0JiqHjQI5uWPVVvzgwtNcnYL9zivy4hbpuWd1Ey4YX28Uciacw4LOyspKEU/HFbbMPJOw371pqCCMDaYRXeaOFR9p8WHlRKa/MLWMGCTe88OakuMvnTQ00FCkE2d0zm10wC/iRUqTOJPmfwtgrFKqL4BhALYA+H2Mxy9Z3JLuw9ChVFZSq1uSqSnR05QUunLjTjw55yz88KKxoexzig07h6ymqtLXIYUtM88Hr4KGqPZNEoU+rAix/cWqa6dhUJ/qbp8N6lONVaP2AsuXG7c3JcffsWorbmgYkzW/Yi44G25H3giWFCWxCS6l1Eal1F7rrUCPVJVVF0dnpWChJiJ2iqW6mipUVQa/ce2b3HYKN1801ihs3CIvfhV+Tvv8cBMbfv3IgPgbg5qqSy+ZOCR78vGK4N8HG5wWJ/Rh8RBW4GT6i1XXTsMPLxrb5Wt6VFai8dSzgaeeAl580XV7v6aspkmt88Eesgy6nJQGsSbNi8inAPwEQF8A7QCujvP4SRJ1LxbncJidjxRkOC/zJrftNOEmhkwhcgV0VeHUZ+REmXKfALPY8BpiE/jnmEWBVyHBhBP7Zy03NWKtq6lC7549UpU3WK6Usw+LiyVXnOmaOO+G018Y/ezHPo+G/14MfPe7QF1dt3145Vw688J6VQo2LtA5ViYbg4ixXJL4SfpJJGleRE4AcDmAJ5VSj7l8PgvALAAYOnTo+Ndeey1eAyPAJDKiTo72Sth0u8k9xVAFsPCEfWj46iWBj9FteyvBFchudQEA/WqrMO/c0a5iwytBvVOpVAgVUwFAdP2I0kkakua9fFgp+q8kMQmb6krBf37itG73jqefnTUWmDcP+MEPgOqjQ4/OHC4/vERXPu0rSOlg8mGJVSmKyEAA/wQwVCllfJRJc5VPZtWb6SoLgFcNFTRhejn52ZEZXelXW4UZp74NKzfuzNqPV8XjDy8ai4YXHgV69gQuvNDVVr+Imi0ww1YEBnGKSYqXoOeTxirYuEmD4AKC+bA0+68ocKsGDBLVcU7nM6hPdVfH+UxM/qvLz27eDPzoR8B//idQdbSxcqZdQfCrenSD9375UIyCazCAJgBjlFIvmNZLq8MKGvUxRbjctnebTDqIyAhiiwCY6dUeQgQ/uNB6mvzFL4DjjgMaGrLWC1LSnYuzGnv9w67DcU7CRgwL4QSva1yPJU9v6eboGbnKnRQJLl8fllb/FQWmh6ZLJw31FF1htvOK0Het//LLwE9+okVXZfdcy6AtKcL6sCij2xRyxYfJh8WSNC8iFSLyb9YTIUTk7QB+DGAzgI1x2BA3QZK8qyoEh460uybRu23f1qm6iS0g2BQRQWxRAJY8vQVTRw1wrSjsUOro1DmXXw688QawLNs5BamyyaVYIIjYAnSuV9DiBNsJek0P5EfjmqYssQVw6o5Soxx9WKHxqgYs1HZeFdu3P70F1zWuB975TuALXwC++U2gs9PHanfC+rCopvcphA8j8RFnW4hzALwgIgcBrAJwCMDZSqn2GG2IDb8k77qaKkD0FDVuN0qYPkx+6wbdlwKw7Pk3sPD8MZ5tIRrXNGHy3pEY/rjC5HnLut3cQULyUQqRutqqwA6oEE7Qq8s/e2mVHGXlwwqNXzVgrttNu+kxDJuzDMPmLMOVS9eiby9zCxpbpE2893UMk/dj2DcfwLA5yzBxwQoAwasPw/qwqOZq5TyN6SIWwaWU6lRKnaOUGqiU6q2UqldKzVRKvRLH8YNSyLYNpj5K9XU1eHXRDPTu2cMzWhWmD5PfumH2Zc9R2GlwcraA0YJG0NQKzL1rTde1CtIHKxcn06/WNNn2UWqqKqEUAjsgrznOgn73XufCXlqlQ1p8WDGTa+8pr+2m3fQYXtpxsNtyt+anNh1KOfLBpGubiQtWuLaAcCOsDzP5gnx9RFRCjkRDnBGuoqbQoVm30LYAmDpqAADvP/aTFz3qOrRXVSFZ/bWC9GwK2xh18fJNnk09swRNp2DxH54LfCznvoMI3Xnnjs4698oKQV1NFQRa6C08fwz2GoYem5pbsvbv5eyCfvemfQjAXlqEZJBr7ymv7Zxiy49KEaMgs5cvueJM9O0Z3Idd17geI+bqSNmIuQ/oYUsHpp59+fqIqIQciQbOpWgRRdsGr2Rqv4q+mqrKrqlxcqlSdCZSDjuuJlBfG0CLhZsvGuua5GnKBRMovHrqPuBTn/KsWKyqFPSu7oG9LW3GORcFeniz3nF+QZJDvZJmM8/Dqz2FjfO7N1WNOvdhFyCwp05upCVpPghMmu9OrlWKpu3CzLsIwHuqHqWw+TvTdBU2gFPnPYR9rd65r6bphNwS+qNIbmermSJEKUhFRXFVKQYlLoflW06cA14ibvb0kb6Vg7mKPTehFwZT64apowYY91tfV4Mn3/EW0NQEXHklYA0DZO6jrrYKBw63oy1jAkFbXJkI6zzCVoc2rmnClUvXuq6T+d17OTYgfKsOYuCllyAnn0zBRQLhJbgunTQ0S6RNOLG/8X4HgM0HHtS9ugYM6FoWpom0TZyTULNKsYg4cgS4/nrI976X7OTVxY6pW3o+oVmv8XXnpMxBtg9yY5mq5oKSGeZ2dq+fvOhRoyidPX0kMO4s4JlngNmzge98B6it7baPyYse7coRs/Gz0zmZtYlMpxhkOrLMaYdM30Hmd++VnPrknLPo4ArBww8DTzyRtBUkRZw0sLfrsOJJA3vjhoYxrm0jTAzqUw1cuxD41reAz38eOPVUAEf9YJhoWq6TUOcinpx+miTEzp3A9dcDX/sa8L3vua7CHC6LKMbY/cbX7bkLTYnmmdsHzTHzqppz49JJQ3XFpEWvKvNPwiQgFTKmJzrjDODqq7Xoev31QNv7Yee1BWnvAACZvs6kvSpEuvYX5LtncmqEKKX7IjU1aaFOSEBWXD0FJw3s3W1ZpQBfmXqS6/pe9+uqa6cBffoAN90EPPAA0NiYs12Vln8JU4TFFg8pZt06YMECYOFCYPhw42oUXBbOCZbtJOx8nhyCirgg6wUt/w0rACac2B+t7Ud70ew51Ga8yb0SxLutP3gwcOONwH/9l5401qIuQKWhiabmFsy+a13g9g42JvGZ2VcsyHfP5NSIaG3V/ZBOOQX43OeStoakkK9MPamb/+xQ5qIXr/v11HkP6f9UVABz5ujf5o03dn+CC8ikd/TD1+5a1008fc3gv2xMPv7aP5jntiVFwL33Ag8+qIV6nz6eq1JwZWBHnF5dNKMgw0RBRVyQ9YJGWMIIgLqaqlB9XGZPH+kaMVJw6UtTUwN8//t6iPG3v9Xr5Zku2NapMP++DVnLg4jMChfDM8/T77uPqsoIKGw7klSxfbuOhn7pS8D73pe0NSSlhPVhJrIS5C+6CJg6VQ8RHQxeDXnppKFYu3UvOjq7O7yOTuUpnkx+7OCRDtfKR5IwnZ3A4sVAR4cW6BX+coo5XBETdHzdb72gOWamZHxncnpVhWD+eaNxlSGB1O3mbxhXb0w4dXUWIsBXv6rD89/5Dva2jHfdNpOqCsExvXpk5XrZuHWcN12bTDoNYi9oRDAz5y7KKiN7GCHzmCXJc88Bt9+uRfkxxyRtDUkxYYb7vXyYK+PHA/X1wDe+AfQJlgS/6p+7cPCIe8TdtBzw9mN3rNrKqudi4uBBnet36aXA6acH3oyCKyW4CanMvl42JmHgtixowngm/WqrXMVQ5nBhduLnaWj45HAM/tUGNFV4R+DaOhVqq82Cy40gFZ8mwkQEo0hO9Xo6L1nBdeedwGuv6eGaAE+FhHgRpuBp2k2P+e7P2Uz1pIG9seKmm4B5K2DOCj1K2N5gNrOnjzSKwVyT8EkEbN6sI1vf/jYwaFCoTSm4ioAglSkN4+rx7Gu7u1UgKgD3rG7ChBP7d1vfJAzclrmJFa+hMtN9by83RmzOH4PZ552KuX/ciBZ4NxVsam5B7+pK16dBt47zQSo+62qq0NreGfg8/ShUKXZZJePv2qUjWlOm6KIKQgpAGB/mJYb69qx07Vz/0o6DmPbjp9CrsgKHOzoRRHSZ8NqyYVw9rlq61jXvNMgctWGYedtT3foyTh7RH0uuOLOgxyg5lAJ+8QudCnHTTV392sLAx8sEaVzThHd960FcuXRtoMqUlRt3Bp4kOWheUNhiAdMk0vZyz4jNe0/GwovGo76ngl9DiCPtnah0JF5VVQrmnTvaeB5PzjkLP7xorGuu1fzzRhesKKKQ1URlk4z/f/+nI1pz5wLnxNOfiJQHhSp4ev76DxsF2Us7DmLjgnPQq9L6k5ljxGnmpKE5fe7XjT+UDQ6xBQBPvrIbM297yrAFwWuv6f6SY8YA116bk9gCGOFKjMY1TZh917puTUBtWto6cOXStVi8fFO3yInXdEDD5yzrirQ4I2FNzS24aulaPPvabtc8gDBDZZUixvD2dY3rjRGmzL5XADD3nufRklEd6aStU6Gupgq9e/YI3ZMGMOdaFWKYrpDDgGEjjKlj925g0SLg/e/XJdOEREBcvag2LtAPCzoSdgBBo11Bu+rbn+fSjT8ophlHgs5EUlYoBfzyl8Abb+jofK9eee2OgishFi/f5Cq2MnEmUHslVdqRltl3r8uaFNv+fMnTW7KGH8PilUtgnDIDLk1EPcSWzd6WNqyd96FwBiLb+drRvkIluxdyGDCqZPyi4L77gCef1FGtfv2StoYQT8JEqHUkzF9s5TJTiVvDVpIAW7YAP/gBcMklwOWXF2SXFFwJEfSPc2bkJEhyuJvYsrHbNzSMq885B6k+QEWgk6BNRJ0UYlgtiirAQs9KUHKdonfv1k+D732v/peQFGD7BTeczVWDUnKpAeWAUsCvfqUbMS9apFscFQgKLg+inKMqSCsDG+dwnG1TLlkE25pbPEVI5v7dzjmXikC3JqJBJpguREJ7hcsQaL5VgH7DgGU7t5lSOlfrqaeAa64B+vdP2iJShnjdf6bJpgEYfdpJA3tjxdVTQtuRqw+LOqF98oj+rsOHk0fwfsU//wn86EfAxRfr6Z0KDCevNhD1LOxeOVxOTJNYmybH9tsXANft+tVW4XBbdiWf85wb1zTha3euC1yqvNkRUne7tlWVgt7VPbC3pQ2DpRWza3ei4aqZaHxpbyjxEnTy6jChfjcHDrgL06h/N0XLE08Ad98NfPSjwNln57UrEeHk1SQngtx/Exes6Ca6BvWpxo79R4zzxLr5CbdqRkBPK9SpFAarBOkAhgAAGbNJREFUw5j9viFo+OhE1+Otunaaq/1uCe1OelVKVy5ZUJw+rLa6opv9ZV+luH07cMstOvXhi1/MO6pl8mGMcBmIuj+SvY/5923oqvDrXV2JI+2d3USY8ykp88apcZn3sKpSPIcVp44agCWGXCu33ldu5+zWosKEXxsHo5DauhWNN9yKuXISWpQ+zyDDgV7T/GQSdJohUzRw4fljXEVw2fXVWr8e+PWvgYkTdak0+2qRBAly/9li57rG9bhj1VZjxAswDwmuuHqKe78uOxJ25Ajwy19i4pzXsR3dE6237z+CiQtWuIquIInrhzsURl37QGDR5ebDaqoq8cOLxpamTwrDvn16Htf2dl2FGHGuKQWXgTj6I7nl7niFw503zqG27onnAuCiM4Z4Jq/fs7oJdYbmpSac59y4pgn3rG7KElsV0r2ju18bB8+bfcgQLK4bixbHsf3ES9DvJ2hgN6yAKpu+Wps3Az/9KTBihM5zqMp9nkxCCkXQ+++6xvWefhLwHxL0HGasrga+9CVsn7PM9WMvkReEwx4P1U7K7iEwCK2twG236erDL39ZzyYQA2UnuILm1xQ6MTooXkLEL3qjoHt1eSW2t7R1oGePCtRUVWaF3Xv2qDBOnRPEjr69wrdx8CIX8RJUTO419BPL14akfjexsXMn8OMf60lav/UtoHduycSEREHQ+++OVVuN+xBr/VLJvSybh8AgdHQAS5boyPwXvgCMjLf9TlkJrjAVa8XYHynIDbKtuQU3XzTWM49pb0sbbr5orGteUpBzNtmRaxsHEybnWaeOAI8/rqvgMjowN65pwoHD7YH3nY8Npu2L8XdTEPbs0U+Era3Av/87cNxxSVtESBZB7z+v/NOwbRxyQwEtLVm5QqaE9nwo+YfAIHR26jY1jz0GzJwJfOYziZhRVoIrTGi1GPsjBYneDK6r6bLRlNhur5M5VGmfp85tUmixhit7ueSJxXUDz54+0rWv2IHKajRu3I2Gxq8DY8cCF14I9OwZqLcZoCfIDiqAwgqoYvzd5MU//qGfCHv1Ai67DHj725O2iBAjQe8/UwPnQk+hY6qKHNSrAliwQAuuz362675acsWZgRPng1KyD4FB2LcP+M1vdKf4GTOAm2/u9pAeN2UluMKGVoupP1KQ6E3mTdTV0d3nRnNG/ZyCbs+htqwoYFw3cMO4+m5FBTZtncDiXX3Q8IMfAGvXAtdfDxx7LLY1nxJsxyHut1wEVDH9bnKisxNYvhz405+AUaOAb34TqK1N2ipCAhHk/rtkonuuayGn0AF0gr65SvEc3bPuf/8XePNNoKEBmDSpW7Vg45om1wmtF33itMA2lNxDYBBefhm4/XYtrj79aeAd70jaIgBl1hbC1EbB1HahmDDZLgJAmXMOnNErpfTQn72+14TPmTivkSkXLpceVF7bDJ+zLFi59o4dmPzDv6KpM3jydn05OJ4wNDXpaNaePcC0acDUqYk8DbItBIkDu0oxqil0QtHeroe8Vq0Chg4FPvUpoF+/QK1/yr6lg01rK/CHPwDPPaeLeT71KZ1rmgAmH1ZWgivNPZICCw8DpnMP2sA0yHFyub5u2wj0JK43NIwJJZL1vp7vGg4NQlq+/8hobQUeflj30Ro8WDupAQMSNYmCi5Q1mzcDd9wB7NuH4eq9gRpcl63oUgp44QXg3nu1aG1oAMaPT9oq9uEC0h1azTdvypS/5jUZddjj5FJ+7LZN5ryPYYYvM7/foA1hy7I8eudO4IEHgE2b9Kz3Z5+tWzskmNtACLEYNkzPP9rejsHfXY6mVv9Nymri6bY24C9/0a/2duDd7wa+9jXgmGOStsyXshJcQHrza0xT6hw60o7GNU2+52TKU+tQCtaopJGg+Vm5lB+bPrPnfbSjWEFFsv39Bu0472dfSaAUsHGjFlm7dgHHH68TSD/72aQtI4SY6NEDsxu8K87Lhj17gAcf1O0cqquB978fuPZa/f8UEYvgEpHvA/gogCEADgBYBuAapVQZyfL8cOtMD+ik9iuXrsX192/AvHNHG4WIKUJWX1fT1X0+U3TZIsyU5+SWd2U6RoWIURR6zauYOYdkWJHsjHZ5RfIG92gHtm7VlUKlEuU5fBh4+mlg5Ur9RDhqlBZYxx+ftGWphD6MJEGo+XPb24EeJRJDUQp48UXgoYd0QUG/fsA55wCXXJJqHx1LDpeIfA/AXQBeAFAH4DcAjiilPua3bZpzIKKYxNgridKZj+RMmD9wuD1r2iB7/TC2mvKu/mVEfzy3Za/r05gpV6pxTROuWrrW1ZEUupjBFPXq17MS8/rsQMO+l/WC6mrglFN0LkAaRFhLC/D888Dq1cC2bdrenj2BM84ApkzR/08RxZjDlasPS7P/Kgc8p+cpQkwtIyb36cSS9ue06AKA4cO1/3r3u4t/FgilgFde0f7rH//QVdKATnyfPh044YRk7cuBokqaF5EZAH6nlDrWb920OqyoEvRNyfM2tkjxnSA6DwForJiETna3K39Mtjm5rnG9a4TNTpwvJI1rmlxbTXT7blpbdeh69Wrg9df1Cg4R1rjujWRyAQ8ePGrb9u16Wc+ewKmnatsGD47ehogpRsHlJKgPS6v/KgdME1CnTXRlJcwrpRPvV68G/v73oyLsxBOB008H3v1uDJv3SNZ+N8fR8LW9HfjnP7VtGzdqcSWixdX48ToSXwJRumJLmv8ggOcTOnYsRDV/ldcQHICuz9yO39ah0Ltnj7y7wXvlXa3cuBOdBhHf1NyC4XOWZQkUW1Rlii4FPe/jhBP7F1TINIyrx+Llm7IEV7fvpmdPYMIE/bJpbQU2bACWL0fjay2Y2z4MLajsOq+5S1cDK1eiYWQ/4G1v06+BA4HKSn+jlAIOHNDzetmvbduAvXuPriOi16ut1cLvggtS+eRXQpS8Dyt13MRW5nJn24hJ7+iHzbtaEi+48q1GFNERruHDjy5TSjf/XL0aw5a+Dv1Imxm5Vxg2535sHn9YP7TZPqxXr2AR/rY2PfSX6b927tSCyv57IKIntx82TIurCy8M5h9LiNgFl4hcAOAKAB/wWGcWgFkAMHTo0JgsKyxRzV9lSp63sTslRzl/ll/eldfnCu5TKq3cuDMrchdVBWFO16ZnT/10ePrpWLzo0exJtVGJxS0D0TCwn84H+9vfgB07jobHbZTKcmCN7f2xuHMItnVWYXCvCsw+YwQaLv8I0Ldv8Q9nliF+PqwU/Fe545zcukOpblElr2nhihIRLXSGDQOecZtQ2xJg73qXFkyPP67/bWnJ9kGZPswWU1VV+gHwbW8DhgwBJk7U7WXKTFD5EavgEpFPAvgZgPOUUs+Z1lNK3QrgVkCH5GMyr6BENf2NfXO7dR8Gjs4RFsXx7Twvrwib/eRnysuycYqpOCZYte032RX02hht3d+m86ZC2jT33vVo6bTm9zysMPfpXcDgA2gY5zviTmImiA8rBf9V7nhNbm2TREuZKPKCu3HKKfpFIiF7oryIEJHPQTuqc5VSK+M6blLMnj4SNVXd1X2Y6W8a1zRh8qJHMXzOMkxe9Cga1zR1fdYwrh71BnFgLw97fK/j2Z/PvXe9p9iy998wrj5Qs75M4WISO7kIRLdz8bM/zHdTSFu9hp5JcVFuPqzUOWlgb+PyIL0JgXhbymT6sMyRAqevJsVLLIJLRL4K4EYA05VST8ZxzKRpGFePheePQX1dDQRaCAVNmHe7sa5cuhbjvvNw183lJ6gaxtXjgvH1XUOMlSK4YLx7e4UgN7KbMLD363Z+JkGYSaZAcTufqkrBwdZ2owh0w3Qu19+/wTgMG+a7Mdma61yScUT2SP6Uow8rdVZcPSVLdNkJ80Ensc53xCIMfDhLP3ENKf4IQDuAlZLxQ1ZKFX9r2DzItcmqSdy4TSRtCi83rmnCPaubup7UOpQyJqEHSfA3CYBOpVyn/PHLNRNrHRvn+dhtLOzk9qA5E6Zz8bIjbOuJQs5YENXQMyk4ZenDSh1TNaJpcutMcn3IypVCPZxtXjQDw+Zk53HFUqVY5sQiuJRSzPwNgdcNlCmEvARd0CrJxjVNvo1HgfDCwD7G9fdvwJ5D3SsC7ZYPTtszz2fyokeztguSMxHW+eQqbAo1Y0GYqYtIctCHlRd25fTvVm1Bp8voYr/aKs9G01FQyIcziqtkSH/DixLEr/VDEFER5GnIHn7zssMmF2GQOc2OWzTIuXzqqAFYuXGnZ0dlv3M3Xbu6miq0tndmNWttam7B5EWPJlbineb5PQkpZW5oGIMbGsZEn6geED6cpR8KriLEbzguyBNNkKch09AlcPRGznQ2x9ZUoVdVBZoPhWuc6hYNcjZmbWpu8Q3hO+13w+SU5p83GsDRqX4y549MusQ7rfN7ElIOFMv9yYez9EPBVYSY5k0Egj/RBHka8ooWLTxfh9Qz99Hc0oaqCkFdbRW2Nbd0JWsWMk/NiyDn7ueUGsbVu3bKT6LEmxBCwlAs4o/kBgVXxOQajvYbjguyPeD9NOQ1obUtTLK61XeqrtyqfCJDYXKtBMg7ohbk2KwMJIQQEhUUXBHiNmwWVqB4iQc/MeYnPNyiYJl5TV55ZDa5Rob88tRsgk5gHUaYsjKQEFJMTFywAtv3H+l6P6hPNVZdOy1Bi0gUxNb4tByJsm9KIZrgZfYKA5CV1xS0LCuXyJBbLysnQYdPw1yLxjVNONjanvOxCCGkkDjFFgBs338EExesyFrXr0E1KW4ouCIkyqGrfMWcfeNeZU0RVFdTlVUZqIBAoitIZMjpKABkNYa9dNLQnBrFBr0WtjBz5sX1q60K1fiUEEIKhVNsmZaz03z64ZBihEQ5dJWPmHMb6jShoMWPaR1nA9Ogx5t773osPH9M6KajbgS9FqZE/drqHhRbhJCiJmhvRVK8MMIVIYWcAsZJPvP5hakQtHOoTJEuBf98tKinpAh6LZgsTwhJK/Rf6YeCK0LymU/Rj3zEXNAbNHN/JlFTKeIb0o7aUQS9FoWcdJoQQgrBoD7VgZbTf6UfDilGTFR9U/Jpgmca6uxXW4Xa6h7dOr8vXr4JVy1di7raKlRVCNoc81x0KOVbeRl2aDVsK4yg14Kdmgkhxcaqa6cFqlKk/0o/FFwpJZ/pJkw3bubcYM68qz2H2lBVKd0qGW388gjCOIpcW2kEEbbs1EwIKUaCtICg/0o/FFwxUOi5uPLt7xXkxnXLu2rrMM1w6D08GMZRRJ0Yyk7NhJC0Qv+Vbii4IqYQzU+dFEKU5NqN3YRfHkFQR8HEUEIIIaUIk+YjJooKvThEiUlA9autiqzy0uu4ufT6Yn8aQgghxQIFV8REIY7iqFYxVf7NO3d0ZJWXXsedOmqAp5hiU0BCSLHCh0ECcEgxcqJofpprtUqYXDK/vKuo8gjcjjt11ADcs7rJc1g2yaaAhc7RI4SUDlGklZB0QsEVMVGU8uZSrZLLTZ9UgqbzuJMXPeorpsJGEgslkuhMCSFesEM8saHgipioSnnDiqEwN32uYiSqSE8QMRUmklhIkURnSgjxgoVAxIaCKwaKoZQ36E2fqxiJMtITREyFiSQWUiTRmRJCvIhyTl2SLpg0XyYETbTPtaoyyvkSg0zdE2YapUKKJE63QQjxIso5dUm6YISrTAgaAcpVjEQZ6Qk6LBs0kljIJ05Ot0EI8YId4okNBVcJ4pVL5XfT5ypGogybFzo3rJAiic6UEOJHMaSVkOSh4Cox/HKp/G76XMVIVJGeKHLDCi2S6EwJIYT4QcGVIoJEevJNCM9VjEQV6cnnfLyuF0USIYSQOKHgSglBIz2FyKUqJjGS6/k0rmnC7LvXdU243dTcgtl3rwPA/liEEELih1WKKSFoFWBSVXNRTa2T6/lcf/+GLrFl09ahcP39G/KyhxBCCMmF2ASXiFwsIo+LyD4RaY/ruKVC0EhPoUqQ/eb+cn4+/74NkbSFyPV89hxqC7WcEC/ovwgh+RLnkOIeALcAqAFwa4zHLQmCVgGacqkAPUVOkPwqv+FLt89N5NsWIorcsMY1TRxWJGGh/yKE5EVsgksptRwARGRKXMcsJcJUATpzsMJW+vklqrt9bqIQQ5m55JTV1VShucU9msW5DklY6L8IIfnCHK6UEKaTupOwXeD9hi+DRq2SbAA6/7zRqKoQ188K1QGfEEIICUpRVimKyCwAswBg6NChCVsTPUEbe+ZaPRi20s9v+NL0eb/aKtRW9yiKBqD2ca9cutb1c851SKKi3PwXISQYRRnhUkrdqpSaoJSaMGDAgKTNiZSoqvsyCVvp55eobvp83rmj8eScs/Dqohl4cs5ZiQ/ZNYyrRz3nOiQxU07+ixASnKIUXOVElJM+24St9PMbvsxneDNuOHEsIYSQYiC2IUURqQRQBaDaet/L+qhVKaWMG5Y4UU76bJNLpZ/f8GUxNUf1gnMdkkJA/0UIyZc4c7g+DeBXGe9tRTEcwOYY7Sgqopz0OZO0CKQoKOdzJwWD/osQkhexDSkqpX6tlBKX1+a4bChGOORFSPFD/0UIyZeirFIsJzjkZSZo9SYhhBBS7FBwFQEc8jqKLbKamlsgAOzkGL9mrYQQQkgxwypFUjRktsgAjootGzYsJYQQklYouEjREGTKIDYsJYQQkkYouEjREERMsWEpIYSQNELBRYoGPzHF6k1CCCFphYKLFA1uLTLs6aeLuZs9IYQQ4gerFEnRwBYZhBBCShUKLlJUsEUGIYSQUoRDioQQQgghEUPBRQghhBASMRRchBBCCCERQ8FFCCGEEBIxFFyEEEIIIRFDwUUIIYQQEjEUXIQQQgghEUPBRQghhBASMRRchBBCCCERQ8FFCCGEEBIxFFyEEEIIIRFDwUUIIYQQEjEUXIQQQgghEUPBRQghhBASMRRchBBCCCERQ8FFCCGEEBIxFFyEEEIIIRETm+ASkUoRWSwiO0Vkv4jcIyLHx3V8QgjJB/owQkg+xBnhmgPgYwAmAni7tey3MR6fEELygT6MEJIzPWI81iwA31FK/RMAROQbAF4WkWFKqc0x2kEIIblAH0YIyZlYIlwiciyAoQBW28uUUq8A2Afg1DhsIISQXKEPI4TkS1wRrr7Wv3sdy5szPutCRGZBP00CQKuIvBChbWnneABvJW1EkcNr5E2xXZ8TkzbAhcA+jP4rFMX22ytGeI28Kcbr4+rD4hJc+61/j3Usr4N+QuyGUupWALcCgIg8q5SaEK156YXXxx9eI294fQIR2IfRfwWH18cfXiNv0nR9YhlSVEo1A9gC4HR7mYi8A/rJ8Pk4bCCEkFyhDyOE5EucVYq3ArhGRIaLSF8A3wewnMmmhJCUQB9GCMmZOKsUFwHoB+AZAD0BrABwaYDtbo3SqBKA18cfXiNveH2CkYsP47X1htfHH14jb1JzfUQplbQNhBBCCCElDaf2IYQQQgiJGAouQgghhJCIKUrBxTnLvBGR74vIBhHZJyLbROQ2EemftF3FiIhUiMhfRUSJyNv9tygvRORsEXlaRA6IyFsickvSNpUC9GFm6L+CQ//lTdr8V1EKLnDOMj86oJN1jwNwGvQ1+lWiFhUv/9/e/YRaOsdxHH9/RkqYYqGmaWYhf7K0MIoUUrqaSE02RkLNRjKhKWXK7GxmYyku2VhpikRZzGU1IqOmUZKFKWkUC8Uk0dfi3Ks7454/99bj9zvnvl+786w+i3M+fXvO73m+zwEXWofoUZJ7gHeB44y+S3uAN1pmWiB22Hj21+zsrzHmsb+6PDSf5ByjnWXLq59vAL4DrvcR7P9Ksh94p6oufSnjtpbkZuAj4ADwFbC3qn5om6ofSU4Bn1bVi62zLBo7bHb218bsr8nmsb+6u8PlzrItuQ9fvniRJDuAN4EjjNavaJ0kVwG3A38kOb16O/6TJHPxxuae2WGbZn9dwv6abF77q7uBi03uXdzukhwADgGHW2fpzGHgfFWdaB2kU9cy+v0fAp4AdgMfAx8muaZhrkVgh83I/hrL/ppsLvurx4FrU3sXt7MkjwCvAw9V1enWeXqR5EbgBeCZ1lk6tvY7e6uqzlTVn8ArwOXAne1iLQQ7bAb218bsr5nMZX91N3C5s2w2SZ4EXgMerKqV1nk6cxdwHXA2yc/AWpmfSfJ0u1j9qKpfge+BjQ5x9newc47YYdPZXxPZX1PMa3/1emj+JeBxYAn4BVgGdlbVUtNgnUjyLPAysFRVX7TO05skVwLrHzPfA5wC9gHfVNVvTYJ1JskRRn9d3A98CzzP6KmoW1YLTVtkh41nf01mf81mHvvr/9yluBlb3bu4XbwK/AWsJPn3YlVd3SxRR6rqAusepU6y9j0/b1ld5DiwEzgJXMHoSagHei2rOWOHjWd/TWB/zWzu+qvLO1ySJEmLpLszXJIkSYvGgUuSJGlgDlySJEkDc+CSJEkamAOXJEnSwBy4JEmSBubAJUmSNDAHLkmSpIE5cEmSJA3MgUvdSLI7yU9JHlt3bTnJSpLLWmaTpEnsL03T6y5FbUNV9WOSg8CJJF8yWta6H7i1qv5um06SxrO/NI27FNWdJMeAR4FdwMNVdbJtIkmajf2lcRy41J0ku4BzwGdVdXfrPJI0K/tL43iGS11JsgN4G/gAuCnJU40jSdJM7C9N4hku9eYosJfR+Yd9wPtJPq+qs21jSdJU9pfG8i9FdSPJvcB7wB1V9fXqtaPAQeC2qvq9ZT5JGsf+0jQOXJIkSQPzDJckSdLAHLgkSZIG5sAlSZI0MAcuSZKkgTlwSZIkDcyBS5IkaWAOXJIkSQNz4JIkSRqYA5ckSdLA/gEVSKHDyg6ZmgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%matplotlib\n",
    "plt.rcParams.update({'font.size': 13})\n",
    "plt.rcParams['figure.figsize'] = 10, 5\n",
    "\n",
    "fig, ((ax1, ax2)) = plt.subplots(nrows=1,ncols=2)\n",
    "ax1.set_aspect('equal')\n",
    "ax2.set_aspect('equal')\n",
    "#print(ax1,ax2)\n",
    "\n",
    "ax1.plot(particles[:,0],particles[:,1],'o')\n",
    "ax2.plot(new_particles[:,0],new_particles[:,1],'o')\n",
    "\n",
    "ax2.plot(feature_position[0],feature_position[1],'ro')\n",
    "arc1 = patches.Arc(feature_position,width =measurement[0]*2,height=2*measurement[0],linewidth=0.5, color='red')\n",
    "ax2.add_patch(arc1)\n",
    "ax1.plot(feature_position[0],feature_position[1],'ro')\n",
    "arc1 = patches.Arc(feature_position,width =measurement[0]*2,height=2*measurement[0],linewidth=0.5, color='red')\n",
    "ax1.add_patch(arc1)\n",
    "\n",
    "ax1.set_xlim(0,7)\n",
    "ax1.set_ylim(0,7)\n",
    "ax1.set_title(\"initial particle set\")\n",
    "ax1.set_xlabel(\"x\")\n",
    "ax1.set_ylabel(\"y\")\n",
    "ax2.set_xlim(0,7)\n",
    "ax2.set_ylim(0,7)\n",
    "ax2.set_title(\"resampled particle set\")\n",
    "ax2.set_xlabel(\"x\")\n",
    "ax2.set_ylabel(\"y\")\n",
    "plt.show()\n",
    "#plt.savefig(\"featurebasedex.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "a map matching example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "gridsize = 1\n",
    "#convert a map value into a simple 1 occupied or 0 not occupied\n",
    "def mapconvert(value):\n",
    "    if value > 0:\n",
    "        return 1\n",
    "    else:\n",
    "        return 0\n",
    "\n",
    "#this map_weigh function is a direct copy from the one used in the later implementation\n",
    "#therefor it may seem a bit more complicated then needed\n",
    "def map_weigh(pose,localmap,globalmap_lib):\n",
    "        #calculate probability of localmap corresponding to global map at particle position \n",
    "        \n",
    "        #find overlapping cells calc m_bar\n",
    "        N = len(localmap)\n",
    "        overlap = []\n",
    "        m_bar = 0\n",
    "        for i in range(N):\n",
    "            for j in range(N):\n",
    "                #only use localmap entries that are actually measured (!=0)\n",
    "                if localmap[i][j] != 0:\n",
    "                    local_x = (i+1/2)*gridsize - N*gridsize/2\n",
    "                    local_y = (j+1/2)*gridsize - N*gridsize/2\n",
    "                    globalmap_x = pose[0] + local_x * math.cos(pose[2]) - local_y * math.sin(pose[2])\n",
    "                    globalmap_y = pose[1] + local_x * math.sin(pose[2]) + local_y * math.cos(pose[2])\n",
    "                    globalgrid = (int(math.floor(globalmap_x/gridsize)),int(math.floor(globalmap_y/gridsize)))\n",
    "                    if globalgrid in globalmap_lib:\n",
    "                        m_bar += mapconvert(globalmap_lib[globalgrid]) + mapconvert(localmap[i][j]) \n",
    "                        overlap.append(((i,j),globalgrid))\n",
    "                    #else:\n",
    "                    #    m_bar += 0 + mapconvert(localmap[i][j]) \n",
    "                    #    overlap.append(((i,j),globalgrid))\n",
    "        #calculate correlation between the two maps (only if there is an overlap)\n",
    "        if len(overlap) > 0:\n",
    "            m_bar = m_bar/(2*len(overlap))            \n",
    "\n",
    "            sum_1 = 0\n",
    "            sum_2 = 0\n",
    "            sum_3 = 0\n",
    "            for ((i,j),(k,l)) in overlap:\n",
    "                    sum_1 += (mapconvert(globalmap_lib[(k,l)])-m_bar)*(mapconvert(localmap[i][j])-m_bar)\n",
    "                    sum_2 += (mapconvert(globalmap_lib[(k,l)])-m_bar)**2\n",
    "                    sum_3 += (mapconvert(localmap[i][j])-m_bar)**2\n",
    "            correlation = sum_1/math.sqrt(sum_2*sum_3)\n",
    "\n",
    "            weight = max(0,correlation)\n",
    "            \n",
    "        return weight"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "#building the global map\n",
    "global_map = np.zeros((40,40))\n",
    "globaldim = len(global_map)\n",
    "for i in range(1,globaldim-1):\n",
    "    global_map[2][i] = 1\n",
    "    global_map[1][i] = 1\n",
    "    global_map[globaldim-3][i] = 1\n",
    "    global_map[globaldim-2][i] = 1\n",
    "    global_map[i][2] = 1\n",
    "    global_map[i][1] = 1\n",
    "    global_map[i][globaldim-3] = 1\n",
    "    global_map[i][globaldim-2] = 1\n",
    "\n",
    "#plt.imshow(global_map,origin = \"lower\",cmap=\"Greys\")\n",
    "#plt.show()\n",
    "\n",
    "globalmap_lib = {}\n",
    "for i in range(1,globaldim-1):\n",
    "    for j in range(1,globaldim-1):\n",
    "        globalmap_lib[(i,j)] = global_map[i,j]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAD4CAYAAADl7fPiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAOMElEQVR4nO3df4xlZX3H8fenoE3KkMouP0R+rMYSEmrK1mywhrSBWikQItrYFtK0tKVZK0tSkzaBtokY+49NY00qq3bVDdgo0l+rJK7IhjZBE38tZBGoWLZkhXEJC2LBjTZm9ds/5oyZZ/be3Zl77p25M/t+JZN7fjz3nOdk2A/n3PvM801VIUnzfma1OyBpuhgKkhqGgqSGoSCpYShIapy82h0YZGZmpjZu3Lja3VgznnrqqYkc9/zzz19y2zPOOGMifdBkHDhwgOeffz6D9k1lKGzcuJFbbrlltbuxZmzbtm0ix13O7+Cmm26aSB80GVu2bBm6z8cHSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSY2pHNGo5dm+ffuS2zryUMfjnYKkhqEgqWEoSGoYCpIahoKkhqEgqWEoSGocd5xCkp3ANcChqnpdt+1u4MKuySuA/62qzQPeewD4PvBj4EhVDZ/uRdJUWMrgpTuA24FPzG+oqt+dX07yfuDFY7z/8qp6ftQOSlpZxw2FqnogyasH7UsS4HeAXx9vtyStlr7DnH8VeLaqnhiyv4D7khTwj1W1Y9iBkmwFtgJs2LChZ7dOLMuZuHU5ba0zemLqGwrXA3cdY/+lVXUwyZnAniSPV9UDgxp2gbEDYNOmTf7XKK2Skb99SHIy8FvA3cPaVNXB7vUQsAu4ZNTzSVoZfb6S/A3g8aqaHbQzySlJTp1fBq4AHu1xPkkr4LihkOQu4MvAhUlmk9zY7bqORY8OSV6VZHe3ehbwpSQPA18DPldV946v65ImYSnfPlw/ZPsfDth2ELi6W34SuLhn/yStMEc0SmoYCpIahoKkhqEgqWEoSGo4m/M64GzOGifvFCQ1DAVJDUNBUsNQkNQwFCQ1DAVJDUNBUsNQkNQwFCQ1DAVJDYc5rwPLmaF5ORwSfWLyTkFSYylzNO5McijJowu2vSfJd5Ls636uHvLeK5N8K8n+JLeOs+OSJmMpdwp3AFcO2P6Bqtrc/exevDPJScB24CrgIuD6JBf16aykyTtuKHTFW14Y4diXAPur6smq+hHwaeDaEY4jaQX1+Uzh5iTf6B4vThuw/xzg6QXrs922gZJsTbI3yd7Dhw/36JakPkYNhQ8DrwU2A88A7x/QJgO2DS0HV1U7qmpLVW2ZmZkZsVuS+hopFKrq2ar6cVX9BPgog8vBzQLnLVg/Fzg4yvkkrZyRQiHJ2QtW38bgcnBfBy5I8pokL2euotQ9o5xP0so57uClrmzcZcDpSWaB24DLkmxm7nHgAPCOru2rgI9V1dVVdSTJzcAXgJOAnVX12ESuQtLYjFo27uND2v60bFy3vhs46utKSdPLEY2SGoaCpIahIKlhKEhqGAqSGoaCpIahIKlhKEhqGAqSGoaCpIYTt55gnIxVx+OdgqSGoSCpYShIahgKkhqGgqSGoSCpYShIaoxaNu7vkjze1X3YleQVQ957IMkjXWm5vePsuKTJGLVs3B7gdVX1S8B/A395jPdf3pWW2zJaFyWtpJHKxlXVfVV1pFv9CnM1HSStA+P4TOGPgc8P2VfAfUkeTLL1WAexbJw0HXr97UOSvwaOAJ8c0uTSqjqY5ExgT5LHuzuPo1TVDmAHwKZNm4aWl5M0WSPfKSS5AbgG+L2qGviPuKsDQVUdAnYxuLycpCkyatm4K4FbgLdU1Q+GtDklyanzy8AVDC4vJ2mKLOUrybuALwMXJplNciNwO3Aqc48E+5J8pGv7qiTzFaHOAr6U5GHga8DnqureiVyFpLGZWNm4qnoSuLhX7yStOEc0SmoYCpIahoKkhqEgqWEoSGo4m/OU2rZt22p3QSco7xQkNQwFSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSQ2HOU+p7du3L7mtQ6I1Tt4pSGosKRSGlI7bkGRPkie619OGvPeGrs0T3QzQkqbYUu8U7uDo0nG3AvdX1QXA/d16I8kG4DbgDcxN737bsPCQNB2WFAqDSscB1wJ3dst3Am8d8NbfBPZU1QtV9T3malAuDhdJU6TPZwpnVdUzAN3rmQPanAM8vWB9ttsmaUpN+oPGDNg2sJqUtSSl6dAnFJ5NcjZA93poQJtZ4LwF6+cCBwcdrKp2VNWWqtoyMzPTo1uS+ugTCvcA898m3AB8dkCbLwBXJDmt+4Dxim6bpCm11K8kB5WOex/w5iRPAG/u1kmyJcnHAKrqBeBvgK93P+/ttkmaUksa0TikdBzAmwa03Qv8yYL1ncDOkXonacU5zHlKOXRZ8z70oQ8tue1NN93U+3wOc5bUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNRzmrBPScoYOL9Vyhhgv5/zLGfLuMGdJY2coSGoYCpIahoKkhqEgqWEoSGoYCpIaI4dCkguT7Fvw81KSdy1qc1mSFxe0eXf/LkuapJEHL1XVt4DNAElOAr4D7BrQ9ItVdc2o55G0ssb1+PAm4H+q6ttjOp6kVTKuYc7XAXcN2ffGJA8zVxnqL6rqsUGNkmwFtgJs2LBhTN1au7Zv377ktssZBpsMquSnlTap2bqXOnz6ueeeG7qv951CkpcDbwH+ZcDuh4BNVXUx8EHgM8OOY9k4aTqM4/HhKuChqnp28Y6qeqmqDnfLu4GXJTl9DOeUNCHjCIXrGfLokOSV6e5Xk1zSne+7YzinpAnp9ZlCkp9jro7kOxZs+1OAqvoI8HbgnUmOAD8ErquqgaXoJU2HXqFQVT8ANi7a9pEFy7cDt/c5h6SV5YhGSQ1DQVLDUJDUMBQkNQwFSQ1nc14HJjUkelLHXY7l9GEtmebr8k5BUsNQkNQwFCQ1DAVJDUNBUsNQkNQwFCQ1DAVJDUNBUsNQkNRwmPMJZlLDa6d52K6WxzsFSY1xTPF+IMkjXVm4vQP2J8k/JNmf5BtJXt/3nJImZ1yPD5dX1fND9l0FXND9vAH4cPcqaQqtxOPDtcAnas5XgFckOXsFzitpBOMIhQLuS/JgV/ptsXOApxesz3bbGkm2JtmbZO/hw4fH0C1JoxjH48OlVXUwyZnAniSPV9UDC/YPKl54VO2HqtoB7ADYtGmTtSGkVdL7TqGqDnavh5grRX/JoiazwHkL1s9lrtispCnUKxSSnJLk1Pll4Arg0UXN7gH+oPsW4leAF6vqmT7nlTQ5fR8fzgJ2deUiTwY+VVX3Liodtxu4GtgP/AD4o57nlDRBfcvGPQlcPGD7wtJxBUxmVk9JY+eIRkkNQ0FSw1CQ1DAUJDUMBUkNQ0FSw1CQ1DAUJDUMBUkNQ0FSw1CQ1DAUJDUMBUkNQ0FSw1CQ1DAUJDUMBUkNQ0FSY+RQSHJekv9M8s0kjyX5swFtLkvyYldSbl+Sd/frrqRJ6zNH4xHgz6vqoW5G5weT7Kmq/1rU7otVdU2P80haQSPfKVTVM1X1ULf8feCbDKj8JGltGctnCkleDfwy8NUBu9+Y5OEkn0/yi8c4hmXjpCkwjlL0M8C/Ae+qqpcW7X4I2FRVFwMfBD4z7DhVtaOqtlTVlpmZmb7dkjSivhWiXsZcIHyyqv598f6qeqmqDnfLu4GXJTm9zzklTVafbx8CfBz4ZlX9/ZA2r+zakeSS7nzfHfWckiavz7cPlwK/DzySZF+37a+A8+GnVaLeDrwzyRHgh8B1XcUoSVNq5FCoqi8xuMz8wja3A7ePeg5JK88RjZIahoKkhqEgqWEoSGoYCpIahoKkhqEgqWEoSGoYCpIahoKkhqEgqWEoSGoYCpIahoKkhqEgqWEoSGoYCpIahoKkRt/ZnK9M8q0k+5PcOmD/zya5u9v/1a4+hKQp1mc255OA7cBVwEXA9UkuWtTsRuB7VfULwAeAvx31fJJWRp87hUuA/VX1ZFX9CPg0cO2iNtcCd3bL/wq8aX7Kd0nTqU8onAM8vWB9lqNrSf60TVUdAV4ENg46mGXjpOnQJxQG/R9/cU2HpbSZ22jZOGkq9AmFWeC8BevnAgeHtUlyMvDzwAs9zilpwvqEwteBC5K8JsnLgeuAexa1uQe4oVt+O/AfVoiSplufClFHktwMfAE4CdhZVY8leS+wt6ruYa7W5D8l2c/cHcJ14+i0pMnpU0tyvpL07kXb3r1g+f+A3+5zDkkrK9N4N5/kOeDbizafDjy/Ct2ZtPV6XbB+r209XNemqjpj0I6pDIVBkuytqi2r3Y9xW6/XBev32tbrdc3zbx8kNQwFSY21FAo7VrsDE7JerwvW77Wt1+sC1tBnCpJWxlq6U5C0AgwFSY01EQrHm8xlrUpyIMkjSfYl2bva/ekjyc4kh5I8umDbhiR7kjzRvZ62mn0cxZDrek+S73S/t31Jrl7NPo7b1IfCEidzWcsur6rN6+B77zuAKxdtuxW4v6ouAO7v1teaOzj6ugA+0P3eNncje9eNqQ8FljaZi1ZZVT3A0X8Bu3CSnTuBt65op8ZgyHWta2shFJYymctaVcB9SR5MsnW1OzMBZ1XVMwDd65mr3J9xujnJN7rHizX3WHQsayEUljxRyxp0aVW9nrlHo21Jfm21O6Ql+TDwWmAz8Azw/tXtznithVBYymQua1JVHexeDwG7mHtUWk+eTXI2QPd6aJX7MxZV9WxV/biqfgJ8lHX2e1sLobCUyVzWnCSnJDl1fhm4Anj02O9acxZOsnMD8NlV7MvYzAdd522ss99br/kUVsKwyVxWuVvjcBawq5vc+mTgU1V17+p2aXRJ7gIuA05PMgvcBrwP+OckNwJPsQbn1hhyXZcl2czcY+wB4B2r1sEJcJizpMZaeHyQtIIMBUkNQ0FSw1CQ1DAUJDUMBUkNQ0FS4/8B4XgWB0BTtnoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#a localmap example that has been generated using the simulate_alldirrobot.py file (robot was at pose=[5,5,0])\n",
    "localmap = np.array([[0.,0.,0.,0.,0.,0.,0.,0.,0.,-0.,-0.,0.,0.,0.,0.,0.,0.,0.,0.,0.],\n",
    " [0.,0.,0.,0.,0.,0.,0.,0.,0.,-0.,-0.,0.,0.,0.,0.,0.,0.,0.,0.,0.],\n",
    " [0.,0.,0.,0.,0.,0.,0.,0.,0.,-0.,-0.,0.,0.,0.,0.,0.,0.,0.,0.,0.],\n",
    " [0.,0.,0.,0.,0.,0.,0.,0.,0.,-0.,-0.,0.,0.,0.,0.,0.,0.,0.,0.,0.],\n",
    " [0.,0.,0.,0.,0.,0.,0.,0.,0.,-0.,-0.,0.,0.,0.,0.,0.,0.,0.,0.,0.],\n",
    " [0.,0.,0.,0.,0.,0.,0.,0.,0.,-0.,-0.,0.,0.,0.,0.,0.,0.,0.,0.,0.],\n",
    " [0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,1.,1.,0.,0.,0.,0.,0.,0.,0.,0.],\n",
    " [0.,0.,0.,0.,0.,0.,0.,3.,3.,4.,3.,3.,3.,2.,2.,1.,0.,1.,1.,0.],\n",
    " [0.,0.,0.,0.,0.,0.,1.,3.,-3.5,-3.,-3.5,-3.5,-1.5,-1.5,0.,0.,-1.,0.,0.5,-0.5],\n",
    " [-0.,-0.,-0.,-0.,-0.,-0.,1.,3.,-3.,-9.5,-10.,-3.,-2.,-1.5,-1.5,-1.5,-1.5,-1.5,0.,-1.],\n",
    " [-0.,-0.,-0.,-0.,-0.,-0.,0.,4.,-3.5,-10.,-10.5,-3.5,-1.5,-1.5,-2.,-1.5,-1.5,-1.5,-1.5,-1.5],\n",
    " [0.,0.,0.,0.,0.,0.,1.,3.,-3.5,-3.,-3.5,-3.5,-1.5,-1.5,-1.5,-2.,-1.5,-1.5,-1.5,-1.5],\n",
    " [0.,0.,0.,0.,0.,0.,0.,3.,-1.5,-1.5,-2.,-1.5,-1.5,-1.5,-2.,-2.,-1.5,-1.5,-1.5,-1.5],\n",
    " [0.,0.,0.,0.,0.,0.,0.,2.,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-2.,-1.5,-1.5,-1.5,-1.5,-2.],\n",
    " [0.,0.,0.,0.,0.,0.,0.,2.,0.,-2.,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5],\n",
    " [0.,0.,0.,0.,0.,0.,0.,1.,0.,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5],\n",
    " [0.,0.,0.,0.,0.,0.,0.,0.,-1.,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5],\n",
    " [0.,0.,0.,0.,0.,0.,0.,1.,0.5,-1.5,-1.5,-2.,-1.5,-2.,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5],\n",
    " [0.,0.,0.,0.,0.,0.,0.,0.,-1.,-1.5,-1.5,-1.5,-2.,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5],\n",
    " [0.,0.,0.,0.,0.,0.,0.,1.,0.5,0,-1.5,-1.5,-1.5,-1.5,-2.,-1.5,-2.,-2.,-1.5,-1.5]])\n",
    "\n",
    "plt.imshow(localmap,origin = \"lower\",cmap=\"Greys\", vmin=-0.5, vmax=0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "300\n"
     ]
    }
   ],
   "source": [
    "particles = []\n",
    "parcount = 300\n",
    "#for i in range(3,globaldim-3):\n",
    "#    for j in range(3,globaldim-3):\n",
    "#        particles.append([i,j,0])\n",
    "        #particles.append([i+1/2,j+1/2,0])\n",
    "for i in range(parcount):\n",
    "    x = random.gauss(7,5)\n",
    "    while x < 3:\n",
    "        x = random.gauss(7,5)\n",
    "    y = random.gauss(7,5)\n",
    "    while y < 3:\n",
    "        y = random.gauss(7,5)\n",
    "    phi = 0\n",
    "    \n",
    "    particles.append([x,y,phi])\n",
    "    \n",
    "#print(particles.index([10,5,0]))\n",
    "particles = np.array(particles)\n",
    "print(len(particles))\n",
    "#plt.plot(particles[:,0],particles[:,1],'.')\n",
    "#plt.imshow(global_map,origin = \"lower\",cmap=\"Greys\")\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "300\n"
     ]
    }
   ],
   "source": [
    "weights = []\n",
    "for i in range(len(particles)):\n",
    "    weights.append(map_weigh(particles[i],localmap,globalmap_lib))\n",
    "weights = np.array(weights)\n",
    "weights /= sum(weights)\n",
    "print(len(weights))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(particles[:,0],particles[:,1],s = weights*10000)\n",
    "\n",
    "plt.imshow(global_map,origin = \"lower\",cmap=\"Greys\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_particles = resample(particles,weights)\n",
    "new_particles = np.array(new_particles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(new_particles[:,0],new_particles[:,1],'.')\n",
    "plt.imshow(global_map,origin = \"lower\",cmap=\"Greys\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(particles[:,0],particles[:,1],'.')\n",
    "plt.imshow(global_map,origin = \"lower\",cmap=\"Greys\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: Qt5Agg\n"
     ]
    }
   ],
   "source": [
    "%matplotlib\n",
    "plt.rcParams.update({'font.size': 13})\n",
    "plt.rcParams['figure.figsize'] = 10, 5\n",
    "\n",
    "#setting up the axes\n",
    "fig, ((ax1, ax2)) = plt.subplots(nrows=1,ncols=2)\n",
    "ax1.set_aspect('equal')\n",
    "ax2.set_aspect('equal')\n",
    "\n",
    "from mpl_toolkits.axes_grid.inset_locator import inset_axes\n",
    "inset_axes = inset_axes(ax2, \n",
    "                    width=\"40%\", # width = 30% of parent_bbox\n",
    "                    height=\"40%\", # height = 30% of parent_bbox\n",
    "                    loc=1)\n",
    "\n",
    "\n",
    "inset_axes.set_xticks([])\n",
    "inset_axes.set_yticks([])\n",
    "inset_axes.set_aspect('equal')\n",
    "\n",
    "#plot localmap\n",
    "inset_axes.imshow(localmap,origin = \"lower\",cmap=\"Greys\", vmin=-0.5, vmax=0.5)\n",
    "#plot globalmap\n",
    "ax1.imshow(global_map,origin = \"lower\",cmap=\"Greys\")\n",
    "ax2.imshow(global_map,origin = \"lower\",cmap=\"Greys\")\n",
    "\n",
    "#plot the old and new particles\n",
    "ax1.plot(particles[:,0],particles[:,1],'.')\n",
    "ax2.plot(new_particles[:,0],new_particles[:,1],'.')\n",
    "\n",
    "#labelling\n",
    "ax1.set_title(\"Initial set of particles\")\n",
    "ax2.set_title(\"Resampled particles\")\n",
    "inset_axes.set_xlabel(\"local map\")\n",
    "ax1.set_xlabel(\"x\")\n",
    "ax1.set_ylabel(\"y\")\n",
    "ax2.set_xlabel(\"x\")\n",
    "ax2.set_ylabel(\"y\")\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "location based ray casting model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#returns probability for gaussian hit distribution\n",
    "#normalizer changed due to cut-offs at max and min ranges\n",
    "def prob_zhit(exp_meas,real_meas,max_range):\n",
    "    std_hit = 0.5 #needs to be inflated to at least half a grid cell\n",
    "    norm_hit = 1./(0.5*np.sqrt(2*np.pi)*std_hit* (math.erf((max_range-exp_meas)/(np.sqrt(2)*std_hit))-math.erf(-exp_meas/(np.sqrt(2)*std_hit))))\n",
    "    return norm_hit * np.exp(-(exp_meas-real_meas)**2/(2*std_hit**2))\n",
    "\n",
    "\n",
    "def range_weigh_particle(measurement,expected_measurement):\n",
    "        #range weigh particles, using beam endpoint model for z_hit\n",
    "        z_hit = 0.4 #likelihood to get correct measurement\n",
    "        z_rand = 0.2 #likelihood to get random measurement\n",
    "        z_short = 0.3 #likelihood meas to short, only needed if unexpected objects are involved\n",
    "        z_max = 0.1 #likelihood to get max reading due to failure of measurement\n",
    "\n",
    "        max_range = 20\n",
    "        gridsize = 1\n",
    "        error = 0.5 #gaussian meas error, always inflate to half the grid size\n",
    "\n",
    "        prob = z_rand/max_range # random prob\n",
    "        if measurement >= max_range: #max range prob\n",
    "            prob += z_max\n",
    "        #to short prob\n",
    "        lambda_short = 0.2\n",
    "        if measurement < expected_measurement:\n",
    "            prob += z_short*lambda_short*np.exp(-lambda_short*measurement)/(1-np.exp(-lambda_short*expected_measurement))\n",
    "        #gaussian error prob\n",
    "        prob += z_hit*prob_zhit(expected_measurement,measurement,max_range)\n",
    "        \n",
    "        return prob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAngAAAFJCAYAAAAIWVCvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd3xcd53v/9dnVEZdsizJknuJHdtxnDiEhJBNARJKFkLb/bFLW1g6LG1hw+aywO6ywFIv7QYSdiHUS5bABUILgZAEUkjsJHZiO3Hcq2TJVu/l8/vjnJFlRbI10syc0fj9fDzmoZkzZ8585Dxy/Pa3mrsjIiIiIrkjFnUBIiIiIpJaCngiIiIiOUYBT0RERCTHKOCJiIiI5BgFPBEREZEckx91AdmkpqbGly5dGnUZIiIiIqe1adOmFnevneg9Bbwxli5dysaNG6MuQ0REROS0zGzfZO+pi1ZEREQkxyjgiYiIiOQYBTwRERGRHKOAJyIiIpJjFPBEREREcowCnoiIiEiOUcATERERyTEKeCIiIiI5RgFPREREJMco4ImIiIjkGAU8ERERkRyjgCciIiKSYxTwRERERHKMAp6IyBTtPNpFV/9Q1GWIiJyWAp6IyBR85fdPcdUX7uY/f7096lJERE5LAU9E5DS2H+ng83fsAKCtZzDiakRETi8/6gJERLLdTffsprQwj7KifIaGPepyREROSy14IiKn0DswzG8eb+TlFyygujTO0IgCnohkPwU8EZFTuHvHUXoHh7lmXQMFecbQyEjUJYmInJYCnojIKfx++1Eqiwu4aFk1+TFTF62IzAoKeCIip/Dg3uNBuMuLkZ8XY2BILXgikv0U8EREJnG0o499x3q4aGk1ACWFeTy49zgP7D4WcWUiIqemgCciMokH9x4H4JnLgoD30RevZXlNKR/4n830DQ5HWZqIyCkp4ImITOKR/W3E82OcM78CgOW1ZXz8Zes41NbLzzcfjrg6EZHJKeCJiExi2+EOVteXU5B34lb57BVzOauujFseOhBhZSIip6aAJyIyAXdne2MHa8PWuwQz4yXr5/Pw/lZauvojqk5E5NQU8EREJtDY0UdbzyBrGiqe9t5Va+twhzufOBpBZSIip6eAJyIyge1HOgAmDHhrGyqoKSvUbFoRyVoKeCIiE9h2OAh4q+vLn/aemXHhkmo27m3NdFkiIlOigCciMoEdTV0sqCqmvKhgwvcvXDqH/cd7aOroy3BlIiKnp4AnIjKBPS3dLK8tnfT9ZyyZA8DD+9SKJyLZRwFPRGQcd2dPSzcrassmPWdNQwV5MRsdqycikk0U8ERExmnu6qerf4hlNZO34BUV5LG8ppRtCngikoUU8ERExtnT3A1wyoAHsHZ+xehkDBGRbKKAJyIyzp6WKQa8hgoOt/fR1jOQibJERKZMAU9EZJw9Ld0U5seYX1V8yvMSu1yom1ZEso0CnojIOLtbulk6t4S8mJ3yvLPqgkkYu8IuXRGRbKGAJyIyzp6W7tN2zwLUVxRRUpjH7uauDFQlIjJ1CngiImO4OweO97C4uuS055oZy2pK2a0WPBHJMgp4IiJjtHQN0D80wsI5pw94AMtry9ilFjwRyTKRBzwze6eZ7TGzPjPbZGaXneLcK8zsPjM7Zma9ZvaEmX1wgvNeaWbbzKw//Pny9P4WIpIrDrb2ALBwzqknWCQsrynlUFsvfYPD6SxLRCQpkQY8M3sV8CXgk8AG4D7g12a2eJKPdAFfBi4H1gL/Afybmb1zzDUvAW4Bvg+cH/78kZldnK7fQ0Ryx8HWXoAkWvBKcYe9x9RNKyLZI+oWvH8Ebnb3b7j7dnd/N3AEeMdEJ7v7Jnf/obtvdfc97v494HZgbKvf+4A/uPsnwmt+ArgrPC4ickqH2oKAt2CKLXiJ7cx2HVXAE5HsEVnAM7NC4BnAb8e99Vvg2VO8xobw3LvHHL5kgmvePtVrisiZ7WBrD1UlBZTF86d0/pK5QUvf/uM96SxLRCQpUbbg1QB5QNO4401A/ak+aGYHzawf2Ajc4O5fH/N2fTLXNLO3mtlGM9vY3NycTP0ikoMOtvZOefwdQHlRAXNKCjjQqoAnItkj6i5aAB/32iY4Nt5lwIXA24H3mdnrpntNd7/J3S909wtra2unWLKI5KqDrb0sOM0OFuMtqi7hgFrwRCSLTK0PIj1agGGe3rJWx9Nb4E7i7nvCp4+Z2TzgX4Hvhscap3NNERF352BrD1esSu4fe4vmlGi7MhHJKpG14Ln7ALAJuHrcW1cTzKadqhgQH/P6/hRcU0TOQMe7B+gbHEmqixZgYXUxh1p7GRk5XeeDiEhmRNmCB/AF4Ltm9iBwL0GX63zg6wBm9h0Ad399+PrdwB7gyfDzlwMfBG4Yc80vAfeY2fXA/wNeDjwH+It0/zIiMrslu0RKwqI5JQwMj9DU2UdDZXLhUEQkHSINeO5+i5nNBf4FaAAeB65x933hKePXw8sDPg0sBYaAXcA/EwbC8Jr3mdnfEK6RF57zKnf/cxp/FRHJAYmAN50xeAAHjvcq4IlIVoi6BQ93v4GTW+DGvnfluNdfBL44hWveCtyaivpE5MxxpH2aAS/s0j1wvIeLllWnvC4RkWRlwyxaEZGscKS9j6KCGBXFyf3bd8GcYszQUikikjUU8EREQo0dwRg6M0vqc/H8POaVF3HgeG+aKhMRSY4CnohIqKm9j3kV8dOfOIGFc4o53KaAJyLZQQFPRCR0pL2P+oqiaX22vrJodAyfiEjUFPBERICREedoZx/105wFO7+qmCPtfbhrLTwRiZ4CnogIcLxngMFhp36aXbT1FUX0D43Q2jOY4spERJKngCciAjS29wFBV+t0zK8KPqdxeCKSDRTwREQYG/Cm10WbWOD4SHgdEZEoKeCJiBAskQJMe5JFQ9jy16iJFiKSBRTwRESApo4+YgY1ZYXT+nxNWZyCPOOwWvBEJAso4ImIEHSt1pUXkZ83vdtiLGbMqyjiiMbgiUgWSHovWjOLA/OBYqDZ3ZtTXpWISIY1dfQxb5oTLBIaKos0Bk9EssKU/qlqZuVm9g4zuwdoB3YCjwONZnbAzL5hZs9MZ6EiIunU2N437SVSEhoqixXwRCQrnDbgmdn7gb3A3wN3AC8FzgdWAZcAHyNoCbzDzH5jZivTVq2ISJok9qGdiYaqIhrb+xgZ0WLHIhKtqXTRPhu4wt0fn+T9B4FvmtnbgTcBVwBPpag+EZG06+4forNviHnTnEGb0FBRxMDwCMd7Bqgpm1lroIjITJw24Ln7X0/lQu7eD9ww44pERDJsdImUyhl20VaFa+G19SngiUikNItWRM54TeG4uZm24CU+39ShcXgiEq2kZtGa2ULgH4EGYA/wCPCIu+9MQ20iIhlxtLMfSEXAi590PRGRqCTbgvdjgkkWvcBFBF2yO8ysPZxhKyIy6zSHgayufGbdqjVlcczUgici0Ut2Hbx1wLPc/bHEgbBVbwOwPpWFiYhkSnNXP0UFMcriSS8NepKCvBhzSwvVgicikUv2brYRKBt7wN0PAgeB21JVlIhIJjV39lNbHsfMZnyt2vIijqoFT0QilmwX7QeAj5tZVTqKERGJQnNnP7UpmvU6ryKuFjwRiVyyAa8TKAWeNLOvmdnrzexcM8tLQ20iIhmRaMFLhXnlRRqDJyKRSzbg3QJUhT8bgI8Dm4FOM3swxbWJiGTE0c6+lAW8uoo4LV39DGs3CxGJULJj8FYCF7n71sQBM5sDXECwfZmIyKwyMDRCa88gtWUzWyIloa6iiBGHY1391M1w2RURkelKNuA9AFSPPeDurcDvw4eIyKxyrDsYL5eyFrzyE2vhKeCJSFSS7aK9Efg3M6tJRzEiIpmWWAMvZWPwtJuFiGSBZFvwfhj+3GFmvyBo0XsE2OzuPSmtTEQkA1Id8Ma24ImIRCXZgLeIYKzdeeHP9wIrADezp9x9bYrrExFJq1QHvMR11IInIlFKKuC5+yHgEPDLxDEzKyEIfNrJQkRmnUTAqykrTMn1ErtZNHWoBU9EojOzfXmAsGv2/vAhIjKrNHf1U1lcQDw/dct51lUU0dypFjwRiU6ykyxERHJKKhc5Tqgrj6sFT0QipYAnIme05s7+0YkRqRJsV6YWPBGJTsoCnpktNjMFRhGZVZq70tGCV0Rzp3azEJHopDKQ7QU2m9nlKbymiEhaNXf2U1uW+ha8ET+xiLKISKalMuD9PfAT4LMpvKaISNp09w/RMzCc8ha82vJgseOjGocnIhGZ0SxaMyt3904Ad785PPyxmRYlIpIJqV4DLyFxveYuBTwRicZMW/DuHn/AzM6Z4TVFRDIiEcDSMYsWTgRIEZFMm1YLnpldC6wDSs1ssbvvH/P2DwgWPhYRyWqJLtRUB7yacExfi1rwRCQi0+2i3QzUAzXAzWa2DDgCHAYGU1SbiEhaJRYjTvUki+LCPMri+WrBE5HIJBXwzKzW3ZvdfR9wk5k94e73hO8tINirdlsa6hQRSbnmrn7yYsacktRsUzZWTVkhLV0DKb+uiMhUJDsG7z4zW554kQh34fND7v6Au3ekrDoRkTRq7uynpqyQWMxSfu3a8ri2KxORyCQb8H5FEPIuGHvQzC43s3tTV5aISPqlY5uyhJqyuFrwRCQySQU8d38v8DngD2b2fDM738x+A/wB2H/qT4uIZJfmrtQvcpwQtOBpDJ6IRCPpSRbu/jkzywN+ARjwU2C9u29NdXEiIunU3NnP2oaKtFy7tixOe+8g/UPDxPPz0vIdIiKTSaoFz8wWmdmNwL8DDwH9wC8V7kRkthkZcVq6BtLXRRte95i6aUUkAsmOwXsK2AC82N0vBa4F/reZfTjllYmIpFFrzwDDI05duK1YqiW6ftVNKyJRSLaL9rXufmvihbvfaWZXAr80swXu/s6UVicikibp2sUiIdGCp8WORSQKyU6yuHWCY5uBS4ErU1STiEjapWsf2oRabVcmIhGa6V60AIQLH1+aimuJiGTCaMBL0yzauaWFJ32PiEgmpSTgAbh7a6quJSKSbuluwSsqyKOiKF9dtCISidOOwTOzPYBP49pfdPcvT+NzIiJp19zZT0lhHqXx6W7JfXo15fHRsX4iIpk0lTvbG6Z57b3T/JyISNo1d6VvF4uE2rI4LZ1aJkVEMu+0Ac/d757sPTMrd/fO1JYkIpJ+zZ3p28UiobY8ztbD2p5bRDJvpmPwnhb+zOycGV5TRCTtjqZxH9qEmrI4LZpkISIRmFbAM7Nrzex/AaVmtnjc2z+YeVkiIunVnIGAV1sep7N/iN6B4bR+j4jIeNMdXbwZqAdqgJvNbBlwBDgMDKaoNhGRtOgfGqa9dzD9XbRlJxY7XlRdktbvEhEZa1oBL1z37iYze8Ld7wEwswXAImBbCusTEUm5lnB/2Ey04EEwoUMBT0QyKakuWjP7VzNbnnidCHfh80Pu/oC7a0SxiGS1dK+Bl1Cj/WhFJCLJjsH7KPBHMztr7EEzi5vZc1JXlohI+iQCV00GZtGC9qMVkcybziSLHwB/MLMVY45VAb9LTUkiIumVCHh1FekNeHPLtF2ZiEQj2TF4DnwWaAHuMrPL3X1P+J6ltDIRkTRJBK65pekNeAV5MeaUFKgFT0QybrqTLD5tZjGCkHcF0Mv0tjMTEcm4lq5+5pQUUJifsu24J1VbHlcLnohkXLIBb7SVzt0/lQh5wN+ksigRkXTKxBp4CTVlCngiknnJBrwPAd2JF+7+iTDk3ZbSqkRE0igT+9Am1JbHeWR/W0a+S0QkIan+CXf/rLt3jzv2ceCLgPakFZFZIRP70CaoBU9EonDagBfuUnFK7v4Jd6+ywKLUlCYiknruTnNnf9qXSEmoLY/TOzhMd/9QRr5PRASm1oJ3v5n9t5ldMtkJZjbHzN5BsIvFS1NWnYhIinUPDNM7OJzRMXigpVJEJLOmMgZvNfBh4JdmNgxsIth3tg+YA6wF1gAPAu9z99vTVKuIyIxlaheLhLHblS2tKc3Id4qInLYFz93b3P2fgAXAO4AnCBY2XgYMAd8GNrj7pQp3IpLtMh3wasLFjlvUgiciGTTlWbTu3gvcGj5ERGalxKLDUbTgiYhkStILHZtZnKDbFuBJd+9LbUkiIukz2oKXoUkWc0vjxEwteCKSWUktk2JmVwMHgEfCxzEz+y8zq0tHcSIiqdbc2U9ezJhTUpiR78uLGdWlcbXgiUhGJbtPz1eBu4GzgMXA64CVwMNm1pDi2kREUq65s5+5pYXEYpnbPlvblYlIpiXbRbsY+Et33x2+Pgj8xMy+SxD+XpnK4kREUi2Tu1gkKOCJSKYl24L3GLBwguP/Abxg5uWIiKRXJvehTagti3NUAU9EMijZgPcD4Btmdta447VAa2pKEhFJn0xuU5ZQWx6npaufkRHP6PeKyJkr2S7aL4Q/t5nZbQQTLWLA3wIfSGVhIiKpNjLiHOuOpot2cNhp7x1kTmlmJneIyJkt2YBXC2wAzg8ffwOcHb73MTN7BbAF2OLuv0hZlSIiKdDeO8jgsEcS8CAY/6eAJyKZkFTAc/djwO/CBzC6Lt65nAh91wDXEex2ISKSNZozvMhxQl35if1oV80rz+h3i8iZKemFjsdz935gY/gQEclaiZmsNRGMwRv7/SIi6ZbsJAsRkVkr0/vQJijgiUimKeCJyBkjqoBXHs8nnh/TbhYikjEKeCJyxmju6ieeH6M8PuPRKUkxMy12LCIZpYAnImeMlnCRY7PMbVOWoIAnIpmUVMAzs5+a2YvNTMFQRGadKLYpS6gtU8ATkcxJNqh1A7cAB83sk2a2Mg01iYikRRS7WCTUlsc1Bk9EMiapgOfurwEagI8DVwFPmtk9ZvZ6MytOR4EiIqnS3NlPTVQteOVxjncPMDg8Esn3i8iZJemuVnfvcPevuftFBAscbwJuBBrN7EYzW5PqIkVEZmpweITjPQORtuABHOsaiOT7ReTMMu2xdGY2H3gp8GJgCLgVWARsMbMPpqY8EZHUON49gHvml0hJSARLjcMTkUxIdpJFgZn9lZn9CtgHvAz4DNDg7m9y92uA1wD/kvpSRUSmL6o18BJO7EfbF8n3i8iZJdnFoI4ABvwA+Gd33zLBOXcArTMtTEQklaLahzZBu1mISCYlG/C+BHze3XvGHrRgUalF7r7f3VuBZakqUEQkFZo7woAX0Ri8GnXRikgGJTsG71+BsgmOVwN7ZlyNiEiaNHUEXaN1FdEEvKKCPCqK8hXwRCQjkg14ky3/XgZoYImIZK2jnf1UlRQQz8+LrAathScimTKlLloz+3L41IFPmtnYLto84CLg0RTXJiKSMk0dfcwrL4q0hrryIrXgiUhGTHUM3rnhTwPWAGMXchoAHgY+l8K6RERSqqmzP7Lu2YTa8jhbDrZFWoOInBmmFPDc/TkAZvYt4L3u3pHWqkREUuxoRx8r62oiraG2XPvRikhmJLtV2RsV7kRkthkZcY529jMvC1rwugeG6e4firQOEcl9p23BM7OfA691947w+aTc/dqUVSYikiLHugcYHnHqIh6Dl1iipaWrn9J4sqtUiYhM3VTuMMcIJlcknouIzCpHO4NJ/tnQggfBjN4lc0sjrUVEcttpA567v3Gi5yIis8XRcJHjuoqIW/C0m4WIZEiy6+ClnJm908z2mFmfmW0ys8tOcW6Dmf3AzJ4ws2Ezu3mCc95gZj7BI9o7u4hEJrHI8TwFPBE5Q0x1DN6UJDsGz8xeRbD92TuBP4U/f21ma919/wQfiQMtwH8Cbz3FpXuAFeNq00LMImeopoi3KUuYU1JIXswU8EQk7aY6Bi9d/hG42d2/Eb5+t5m9EHgHcP34k919L/AeADP7q1Nc1929McW1isgs1dTZx9zSQgrzo+20yIsZc0sLFfBEJO2SGoOXSmZWCDyDpy+Q/Fvg2TO8fLGZ7SPYZeNR4CPu/sgMrykis9TRjr7R7tGoabsyEcmEKP85W0MQwJrGHW8C6mdw3SeBvwdeCvwtwR6595rZyolONrO3mtlGM9vY3Nw8g68VkWwVrIGXHcNwtdixiGRCNqyD5+Ne2wTHpn4x9/uB+0cvZnYfQSveuwm7d8edfxNwE8CFF1447e8VkezV1NHH6vryqMsAgnGATxzpjLoMEclxUa6D1wIM8/TWujqe3qo3be4+bGYbgQlb8EQktw2POM1Z1oLX0tXPyIgTi1nU5YhIjopsHTx3HzCzTcDVwI/GvHU18ONUfY+ZGbAe2Jyqa4rI7HGsq58Rj34NvIR5FUUMjTjHewaoiXhWr4jkrmnvlWNmZQDu3jWD7/8C8F0zexC4F3g7MB/4evgd3wm/4/Vjvvf88GkFMBK+HnD3beH7HwMeAJ4Kz3kPQcB7xwzqFJFZKrFESl2WTLJI7KbR2N6ngCciaZN0wDOz9xEsb7IgfH2YIKh90d2TGsPm7reY2VzgX4AG4HHgGnffF56yeIKPjZ8N+xJgH7A0fF1FMKauHmgPz7/c3R9MpjYRyQ0ntinLnhY8SNRVGW0xIpKzkgp4ZvYZggWGP8uJiQyXAB8lCGjXJVuAu98A3DDJe1dOcOyUg1bc/f3A+5OtQ0RyU6IFL+p9aBPqK4OA19iumbQikj7JtuC9GXizu9865tidZvYkcCPTCHgiIunU1NGHGVnTHVpTFscMGju0uY6IpM901sHbMsmxyPe1FREZ70h7LzVlcQrysuMWVZAXo6YsTlO7Ap6IpE+yd7zvAO+a4Pg7gO/OvBwRkdQ60t7H/MrsGH+XUF9RRFOnAp6IpM9UFjr+8rjzX2tmLyCYqQpwMcHM1++nvjwRkZk50t7HWbVlUZdxknkVcQ629kZdhojksKmMwTt33OtN4c8l4c/G8LE6VUWJiKSCu3OkrZe/OKsm6lJOMq+iiE37WqMuQ0Ry2FQWOn5OJgoREUm1zv4hugeGmV+VfV20rT2D9A8NE8/Pi7ocEclB01kHLx+4iGCNusIxb7m7axyeiGSNI23BOLeGyuKIKznZvHBM4NGOfhZVl0RcjYjkomTXwVsN3AYsA4xgL9l8YBDoRxMtRCSLHG4Pxrk1ZNkki8Rix40dfQp4IpIWyc6i/SLBGLxKoAdYA1wIPAq8MrWliYjMTGO4FElDVXa14NWHAa9Ja+GJSJok20X7TOAKd+82sxEg390fNrPrgK8Q7PkqIpIVjrT1YpY9+9AmJAJeo9bCE5E0SbYFzwha7gCaCfejBQ4CZ6WqKBGRVDjS3kddefYscpxQUZxPPD+mFjwRSZtkW/AeB84DdgMPAh8ys2HgLcDOFNcmIjIjR9r7sm6CBYCZUV9ZNLpProhIqiX7z9pPELTiAfwLsAj4A/B84D0prEtEZMYOt/dm3QSLhHkVRdqPVkTSJqkWPHe/fczz3cBaM6sGWt3dU12ciMh0uTuN7X1cuaou6lImVF9RxOaDbVGXISI5Kul18BLMrAzA3Y+nrhwRkdTo6B2iZ2A4i1vw4jR19OHumNnpPyAikoSkRx6b2fvMbD/QDrSb2QEze7/pDiUiWWR0Dbws28UiYV5FEX2DI3T0DkVdiojkoGQXOv4M8Fbgs8D94eFLgI8CDcB1Ka1ORGSaRtfAy8JJFgD1lScWO64sKYi4GhHJNcl20b4ZeLO73zrm2J1m9iRwIwp4IpIlDrVl5y4WCYm6jrT3cnZ9ecTViEiumc7iUFsmOZZdC02JyBntYGsvBXk2ui1Ytpkf7q5xuE0zaUUk9ZINZd8B3jXB8XegfWhFJIscbO1hflUxebHsHB5cV15EXsw4HLY0ioik0mm7aM3sy+POf62ZvQB4IDx2MTAf+H7qyxMRmZ6Drb0smlMSdRmTyosZ9RVFo13JIiKpNJUxeOeOe70p/Lkk/NkYPlanqigRkZk62NrDVWvmRV3GKS2oKlbAE5G0OG3Ac/fnZKIQEZFU6R0YpqVrgIVzsnMGbcL8qiI27muNugwRmaH/eegAt29tfNrxxXNL+OiL10ay1qUmRohIzjnU1gPAwizuogVYMKeYxvY+hke0EZDIbPZ/H9rPn/ccp6mzb/Tx+OF2vnXvXvqHRiKpKemdLMxsHsFEi7WAA9uAG9y9KcW1iYhMy4HjQbfnoupsb8ErZmjEae7sH10XT0Rmpw2Lq/jumy4eff21u3bx6d88EVk9SbXgmdmlwE7g1UAv0Ae8BnjKzC5JfXkiIsk72Do7WvASS6VoHJ6IpFqyXbSfA/4vsMrdX+furwNWAT8EPp/q4kREpuNgay+FeTFqy+JRl3JKC0bXwlPAE5HUSraL9nzgDe4+2qHs7iNm9gXgkZRWJiIyTQdbe1kwp5hYlq6Bl5DYzUIBT0RSLdkWvHZg2QTHlwFtMy9HRGTmDrT2ZP0MWoDyogIqivLVRSsiKZdswPsh8N9m9hozW2ZmS83stcA3CLpuRUQid7C1N+vH3yXMrypWC56IpFyyXbTXAQZ8c8xnB4GvAf+cwrpERKalu3+I493ZvwZeQrDYsfajFZHUSirgufsA8F4zux5YQRD2drp7TzqKExFJ1t5j3QAsnVsacSVTM7+qWIsdi0jKTbmL1swKzOzPZna2u/e4+2PuvkXhTkSyyd6W4Ja0tGb2dNG29w7S1T8UdSkikkOmHPDcfZBgMoWWXBeRrJVowVsyS1rwEl3JibX7RERSIdlJFt8G3pKOQkREUmFvSze15XHK4klv1BOJxdVBS+P+Ywp4IpI6yd4BS4HXmNnVwCage+yb7v6eVBUmIjIde491s2yWtN4BLJkbBrzjCngikjrJBrw1wMPh8+Xj3lPXrYhEbk9LD89dXRt1GVNWWVxAeVG+Ap6IpFSys2ifk3huZmXhsa5UFyUiMh2dfYO0dPXPmvF3AGbGkrklCngiklLJjsHDzN5nZvsJdrVoN7MDZvZ+M8vuPYFEJOftC8exLauZPQEPgnF4GoMnIqmUVAuemX0GeCvwWeD+8PAlwEeBBoKFkEVEIjHb1sBLWFRdwu+2HWV4xMnL8v1zRWR2SHYM3puBN7v7rWOO3WlmTwI3ooAnIhHa2xIGvFmyBl7CkupSBoZHaOzoY0HV7NiBQ0SyW9JdtMCWSY5N51oiIimzu6WbeRVxSgtkMa0AABw1SURBVApnxxIpCVoqRURSLdlQ9h3gXRMcfwfw3ZmXIyIyfTuPdrGyrjzqMpJ2YqmU7tOcKSIyNcn+MzcOvNrMXgA8EB67GJgPfN/Mvpw4UWviiUgmjYw4O4928apnLoq6lKQ1VBaRHzPNpBWRlEk24K3mxDp4S8KfjeFjzZjztCaeiGTUobZeegaGZ2ULXn5ejAVzikdnAYuIzNS018HLRY3tfexq7mJFbVnUpYhIknYeDZbkXDlvdv7/u7haa+GJSOpoYsQYzV39PO/zd/PKr93HLQ/tp6t/KOqSRGSKnjraCcDKutkZ8JbOLWVPSzfu6gARkZlTwBtjTUMF179oNW09A3zox4/xzP/4HR/4n838efcx3XRFstxTTV3UlsepKimMupRpWVFbSmffEM1d/VGXIiI5YHatJZBm+THjbVes4K2XL+fh/W3cuukAt20+wo8fPsiSuSX89TMW8ooLFjJf61SJZJ0dR7tmbesdwIqw9l1Hu6krL4q4GhGZ7dSCNwEz4xlL5vCpV6znwQ8/jy/8f+fRUFnE5367g0s/fSd/e9MD3PLQftp7B6MuVUQAd2dnU+fsDnjh2N9dzdreW0RmTi14p1FSmM8rLgha7vYf6+EnjxzkZ48e5kM/foyP/HQrz11dx8s2zOfKs+soKsiLulyRM9LB1l66B4ZZOW/2zaBNqK8oorggj93NWgtPRGZOAS8Ji+eW8L6rVvHe561ky8F2fvroIW7bfITfbG2kvCifa9Y18NLz53Px8rnaT1Ikg7Yebgdg3YLKiCuZvljMWF5bqhY8EUkJBbxpMDPOW1TFeYuq+PA1a7hv1zF++ughfrHlMLdsPEB9RREvOreevzy3gQsWzyGmsCeSVlsPd5AXM1bXz94WPAi6aR/e3xp1GSKSAxTwZig/L8blq2q5fFUtvS8b5nfbm/jZo4f5/gP7+da9e6mvKOKF6+q55twGLlyisCeSDlsPd7CitnTWD5NYXlvKbVsO0zc4POt/FxGJlgJeChUX5vGS8+bzkvPm09k3yO+3H+VXjx3hBw/u5+b79lJXHudFibC3tFrduCIpsvVwO89eURN1GTO2orYMd9jT0s2ahoqoyxGRWUwBL03Kiwp42YYFvGzDArr6h/j99iZ+9dgRfvjQAb59/z5qy+O84Jx5PH9tPc9aPpfCfE1oFpmOlq5+mjr6OWf+7A9EY2fSKuCJyEwo4GVAWTyfl56/gJeev4Du/iHufCJo2bt100G+98B+yuL5XHF2Lc9fO48rz66jsrgg6pJFZo2thzsAWJsDAW95bSkxgx2NnbA+6mpEZDZTwMuw0nj+aDdu3+Awf3qqhTu2NfH7J5r45ZYj5MeMi5dXc/WaeVy1dh4L55REXbJIVnv8UDCD9pyG2TuDNqGoII9lNaVsb+yMuhQRmeUU8CJUVJDHVWuDIDcy4jxyoI07tjXxu+1N/Ott2/jX27axpqGCq9fUceXqOs5bWKVxeyLjPLyvlRW1pVSW5EbL9+qGCrYcbIu6DBGZ5RTwskQsFuye8Ywlc/jnF61mT0s3d2xr5I5tTXz1Dzv58p07mVNSwGUra7ny7GDWbk1ZPOqyRSLl7jy8v5Wr1syLupSUWdtQwS+3HKGzb5DyotwIrSKSeQp4WWpZTSlvvXwFb718BW09A9zzVAt3PXmUe3Y08/PNhwFYv7CSK1YFge/8RXPUuidnnD0t3bT2DPKMJXOiLiVlEmv5PdnYyYVLqyOuRkRmKwW8WaCqpJBrz5vPtefNZ2TE2Xakg7uePMpdTzbzf/6wk6/cuZPK4gIuW1nDZStruPSsGo3dkzPCw/uDrswLcijgJWbPbj/SoYAnItOmgDfLxGLGugWVrFtQyT88dyXtPYP8cWczdz3ZzN07mvnFliMALJlbwqVn1XDpihouWTGX6tLCiCsXSb1N+1opL8rnrHB5kVzQUFlERVG+JlqIyIwo4M1ylSUFvHj9fF68fj7uzs6jXfxpZwv37jzGzx89zA/+vB+zYFzPpWcFrXvPXDqHkkL9p5fZb+Pe4zm3HaCZsaahgu1HOqIuRURmMf0tn0PMjJXzylk5r5w3XrqMoeERthxq596nWrh3Vws337uXm+7ZTUGesWHxHC5eVs3Fy+ZywZIqBT6ZdZo6+njqaBd/9YyFUZeScmsaKrjloQMMj7jG1orItOhv9RyWnxfjgsVzuGDxHN79vJX0Dgzz0N7j3Lurhft3HeOGu3bxlTt3kh92+168rJqLllVz4ZLqnFlyQnLXfbtaALj0rNm/Rdl45y2q5Ob79rKjqVM7WojItCjgnUGKC/O4fFWwxApAV/8Qm/a18uCeYzy45zjfuncvN96zGzNYXV8xGvieubSa2nItySLZ5d6dx6gqKWBtDgag8xcFk0Y2H2hTwBORaVHAO4OVxfO5YlUtV4SBr29wmEcPtPHgnuM8uOc4tzx0gJvv2wvA4uoSLlhcxQVLghbB1fXl5Odp/1yJhrtz784Wnr1ibk6Nv0tYOreEyuICHj3Qxt9ctDjqckRkFlLAk1FFBXk8a/lcnrV8LgCDwyM8fqidh/Ye5+F9bdy36xg/fTRYg6+4II/1CytHA9+GxVVaeFky5qmjXRxp7+PdZ9VGXUpamBnnLari0QPa0UJEpkcBTyZVkBdjw+I5bFgcdBe5O4faenl4fxsP72vlkf2tfOOe3QyNOBAszbJhURUbFs/h3IWVrG2ooKggL8pfQXLUb7c2AnDVmrqIK0mf8xdV8dU7n6K7f4jSuG7VIpIc3TVkysyMhXNKWDinhGvPmw8E3bqPHWrn4X2tPLy/lXvHtPLlx4xV88pZv7CScxdWsn5BFWfXl1OYr65dmZnfbmtiw+Iq6iqKoi4lbTYsqmLE4bFD7aOt6iIiU6WAJzNSVJDHM5cGEzEgaOVr7Ohj84F2HjvUxpaD7fxmayM/fOgAAIV5MVY3BKFv/YIqzl1Yycq6Mo3nkyk73NbLloPtXPfCs6MuJa3OX1QFBGv9KeCJSLIU8CSlzIyGymIaKot54bp6IAh9B1uDv5S3HGpjy4F2fvbIYb73wH4ACvNjnD2vnLUNFaxpKGft/EpWN5RToY3WZQK/eizYreUF59RHXEl6zSktZHV9OffvPsY/PHdl1OWIyCyjgCdpZ2Ysqi5hUXUJf7m+AYCREWfvsW62HGxn6+F2th/p5I7tTdyy8cDo5xZVF7OmvoK18ytY01DB2oYKFs4pxiz3Zk3K1Lg7t246yHmLqliRQ9uTTeaSFXP5wZ/30z80TDxf41lFZOoU8CQSsZixvLaM5bVlvGzDAiD4y/toZz/bDnew7Ujw2H6kgzu2N+HBPA7Ki/JZ01DB2fPKWTWvjJXzylk1r1x77Z4hth7u4InGTj7+snVRl5IRlyyfy7fu3cuj+9u4WN20IpIEBTzJGmbGvIoi5lUU8ZzVJ2ZH9gwM8WRj52jg23a4g58+cojO/qHRc2rKClkVhr2V88qC53Xl2pEjx/zwof0U5sV4SdgSnOsuXjYXM7h/9zEFPBFJigKeZL2SwvyTlmuBE5M5nmzs5KmmLnY0dbLjaBc/2niA7oHh0fPqyuOjoW9lXTnLakpZUVtKbXlcXb2zTGv3ALduOshLz59PVcmZ0WJbWVLAOfMruG/nMd53VdTViMhsooAns9LYyRxXnn2itW9kxDnc3nsi9DV18dTRTn744AF6B08Ev/J4PstqS1leUxp2FZeyvKaMZTWlFBdqrFM2+v6f99E3OMKbL1sedSkZdeWqOr529y7aegbOmGArIjOngCc5JRY7sVbf2G7eRPDb3dzN7uYu9rR0s7ulm4f2to6u25ewoKqYZTWlYegrZWlNKUvmlrKgqlhr+EWkvXeQ//rTHq48u5az68ujLiejnremjq/+YSd3Pdk8Ol5VROR0FPDkjDA2+F2+6uTtrXoHhsPA1zUaAHe3dPOThw/RNWacX8xgflUxS+aWsLi6lCVzS1hSXcLiuSUsmVtKmXYbSJuv372L9t5BPvj83F77biLnLQy2Afzd9iYFPBGZMv2NJGe84sI81s4PlmMZy91p7uxn3/Ee9h3rYf+x7tHnt29t5Hj3wEnnzy0tDMJedQmL55aypLqEBXOKWVBVTENlkRZznqadRzv57z/t4eXnL2Ddgsqoy8m4WMx43uo6fvXYEQaGRtSKLCJTooAnMgkzo66iiLqKotGdOsbq6Btk/7Eg8O073j36/KG9rfxs8+HRpV0A8mJGfUURC6qKR0Pfwjknns+vKta+vRMYGh7hAz/aQmlhHtdfsybqciLz/HPmccvGA9yzo5mr1s6LuhwRmQUU8ESmqaKogHULKidsVeofGuZQay+H2no51NrLwTHPH9xznCPtvYz4yZ+pKYuPhr6FYeibV1FEQ2XwmFsWJy92Zs38/fdfbGPzgTa++uoN1JbHoy4nMpevqqW6tJCfPHJQAU9EpkQBTyQN4vl5ows5T2RweITG9r7R0DcaBNt62HqonTu2NjEwPHLSZ/JjRl15nPrKIhoqi6mvLKK+oih8XTS6hmAudOG5O1+9cyffuX8fb7t8OS9ePz/qkiJVkBfj2vPm84MH99PeM6j1HUXktBTwRCJQkBcb3b5tIiMjzrHuAZo6+jjS3kdje2/ws6OPxvY+th/p4M4njp609EtCTVl8NPDVlsepK49TGz4Sz2vK4lnbJdw3OMwnfrmd7z6wj1dsWMB1L1wddUlZ4eUbFnDzfXv5+ZbDvO5ZS6IuR0SynAKeSBaKxWw0lE02scDd6egbonE0+AUhMBEKD7b28Mj+Vo73DJw0HjChsrgg+I6yOHUVwc/a8sTzE+GwsriAWAa6ht2du3Y086lfbWdHUxdvuWwZ179oTUa+ezZYv7CSdQsq+PZ9e3nNRYv15yIip6SAJzJLmRmVxQVUFheccm24weERjncPcLSjn+auPpo7+8Pn/cHzzn4e2d/G0c4++gZHnvb5vJgxp6SA6tLCcY841SUFVJfFmVtayJySQuaWBT+n2k3cOzDM44fbeWDXMX62+TA7j3axuLqEb73hmSetYyjBf+83/cUy3n/LZu5+qpnnnK0/HxGZnAKeSI4ryIuNjs+DyZcZcXe6+odGQ1/i5/Hufo53D4w+nmjspLV7gLbewQlbBgHK4vmUFOZRUphHcWHwvKggxtCwMzA8Qv/gCE0dfRwbs9TMM5bM4dOvPJeXb1iYE+MI0+Evz53Pf/76Cb5+1y6uXFWr7fZEZFIKeCICBC1E5UUFlBcVTDo5ZKyh4RHaegdPCn/HugdoDZ/3DgzTOzhMz8AwvYND9A2OkB8zyuL5zC2Ncd6iKhZUFbFqXjkXLq2mulTbcJ1OYX6Md155Fh/7+Vbu2qFWPBGZnAKeiExLfl6MmrJgwoZkzt9etJhv3ruH//zVE/zFWTUUaAFtEZmA7gwiIrNIYX6M61+0hiebOrnx7l1RlyMiWUoBT0RklnnhunpevL6BL/3+KbYcbIu6HBHJQgp4IiKz0L+/dB115UW85TsbaWzvi7ocEckyCngiIrNQdWkh//V3F9LVN8Sr/+sBDrf1Rl2SiGQRBTwRkVlqTUMF33rjRTR39POKG+7jgd3Hoi5JRLKEAp6IyCx20bJqfvi2Z1FUEOPV33iA63+yRV22IqJlUkREZrtz5lfyi/dcxud/+yTfe2Af/7PxIM9dXceL1tXz7BU11FcWRV2iiGSYAp6ISA4oi+fzsZecw99fuozvPbCPHz98iDu2NQFQUxbnrLpSFleXUF0abC1XWVxAvCBGPD9GYX6MeH4eeTHDCBa9NiN8DjD2tREzCM4UEYCe/mHK4tkVqbKrGhERmZFF1SVcf80aPvTC1Wxv7OD+XcfY0dTJzqNd3L2jmePdAwwOT7LHnIhM29KakqhLOIkCnohIDorFjHPmV3LO/JP3H3Z3OvuHaO8ZZGB4hIGhEfqHRugfHGZ4xHHAHRwPfwafcYDw+MhIBL+QSJZbv3Dyvb6joIAnInIGMTMqigqoKCqIuhQRSSPNohURERHJMQp4IiIiIjkm8oBnZu80sz1m1mdmm8zsstOcf0V4Xp+Z7Tazt8/0miIiIiK5JNKAZ2avAr4EfBLYANwH/NrMFk9y/jLgV+F5G4BPAV8xs1dO95oiIiIiuSbqFrx/BG5292+4+3Z3fzdwBHjHJOe/HTjs7u8Oz/8G8G3ggzO4poiIiEhOiWwWrZkVAs8APjfurd8Cz57kY5eE7491O/B3ZlZAsA5nUtc0s7cCbw1f9pvZ41P6BUTkTFUDtERdhIhktdH7RPGn0/o9SyZ7I8plUmqAPKBp3PEm4KpJPlMP/G6C8/PD61my13T3m4CbAMxso7tfOMX6ReQMpPuEiJxONtwnou6ihWAdzbFsgmOnO3/88WSvKSIiIpIzomzBawGGCVrlxqrj6S1wCY2TnD8EHCMIcsleU0RERCSnRNaC5+4DwCbg6nFvXU0w83Ui9/P0rtargY3uPjjNa4510xTOEZEzm+4TInI6kd8nzD26nstwSZPvAu8E7iWYJfsm4Bx332dm3wFw99eH5y8DHge+AdwIXArcAPytu/94KtfM3G8nIiIiEo1I96J191vMbC7wL0ADQXi7ZkwQWzzu/D1mdg3wvwmWPTkMvCcR7qZ4TREREZGcFmkLnoiIiIikXjbMohURERGRFFLAC2n/WpHZzcyuMzOf4PHvUdcmIrPTbL6vqIuW0YkZ3yOYmPGn8OcbgbXuvj/K2kRkasysHCgdc+iDwGuAy9x9ZzRVichsNpvvKwp4gJn9Gdji7m8Zc+wp4FZ3vz66ykRkOszsQ8B7gOe6+5MR1bAUOM/df5bk524HPuDu2jZRJItkw30lGWd8F+2YPXHH73F7qj1xRSRLmdn1BDfh50R8E34BsG4an1sNPJHiWkRkBrLovjJlZ3zA49R74o7fEUNEspiZfZhgiMUV7r5jzPGVZvYLM9toZn80s/rw+L1m9qzw+U1m9sHw+U/N7BYze8jMdpvZRae5ziIz+5mZPWJmW83sCuBTwBvM7FEzqzzFZ9ea2Z/MbLOZvR847u5DGfxjE5FTOMV95TYz+6KZPWBmT5nZRWb2/8xs75h7yavN7M9m9piZ3W5mJeHxCe89qaSAd4L2rxWZxczsI8DbgCvHjo0xszjwdeCd4ebf3wbeFb79b8D/MrP3AAXu/rnw+Hpgs7s/E/gw8E+TXSfsBfgV8EV330DQ8v8AsAV4vrufD/RN8tl84JvA2939POASYHM6/nxEJHmT3VdC64Dt7v4s4A/Al4A3AM8h2GAB4HZ3v9jdzwV2cWKnrcnuPSkT6ULHWWI6e+KKSBYJ/4X9XuBaoDvROga0AS8F1gA/NzOAQuCHAO7+WzP7D+AVwPPDa5UBceAz4TW2EtysXzbJdV4OPOjufwiv2R5eZymwN7zGZJ99BcFWi4nxdk8Cran4MxGRmTnNfaUAiLn7jeGxXuAr7t4ettJ1hMffZGZ/TfD//CLgZzDxvSfVzviA5+4DZpbYv/ZHY966GvjxxJ8SkWxhQWK6Dqgg2J5wrKuAc4F/cvfvT/DZCwiGaTwe7mUNQevd1jHdpBsIWuMmvE54k35w3LGFQKOfmMV2qs+ObbF7BvD50/7SIpJWU7iv9HDy//fnAl8On68DHjOzvwPWApe7e6+Z7QC2hdef6N6TUuqiDXyBYKzMm81sjZl9CZhP0KUiIlnMA5XubhM8fg80AleHN2zM7Nzw53yC7tGrgAYzOz+85HpgqZkVmNkcgu7cL0x2HYKW/nXhsZgFWyUuIthKMWGyzx4b89nnAS8kCJMiEqEp3FfWcfL/q8uA3eHzc4HHwnPuD8Pd24Badz9wintPSingEexfC7yPYP/aR4G/QPvXiuSKbwGVwHYzexR4o5kVE7TYv9fddwOfBD4Snr+eoPX+T8AfgX9z94MTXSc8/2ZguZltBR4m6IrdBiwJB1afe4rPfhd4tpk9QjB256C7N6fnj0FEUuhcwoBnZg2c3GK/jiDgfRf4ZzO7B2ggaNU71b0npbQOnojIGGb2J+B17r4n6lpERKZLLXgiIidbzInJESIis5Ja8ERERERyjFrwRERERHKMAp6IiIhIjlHAExEREckxCngiIiIiOUYBT0RERCTHKOCJiIiI5BgFPBEREZEco4AnIiIikmMU8EREzjBmdpeZfTXqOkQkfRTwRERSxMwKo65BRAQU8EQkC4QtSl8zs8+b2XEzazaz95pZ3Mz+j5m1mdl+M3vdmM+YmV1nZrvMrNfMHjOz14677gvN7I9m1hpe93YzWzPm/cvN7AEz6zKzdjP7s5mtG1fXV8dd82Yz+8W4uj9nZs3AvUnUlvTvPJVrh9e9wcw+aWYtZnY0rC+WqB+4AniXmXn4WDrJf5frxpwz9vHvU/nvKiLRUcATkWzxGqATuBj4T+CLwE+BHcCFwLeB/zKz+eH5/wG8CXgXsBb4FHCjmf3lmGuWhte5CLgSaAduM7NCM8sHfgb8CTgv/N4vAcNJ1v1awIDLgNcnUdt0fuepXvs1wBDwbOAfgPcBrwrfey9wP/AtoCF8HJjkd/vamHMagM8DjcB3Tv1HIiJRM3ePugYROcOZ2V1A3N0vCV8bcBS4392vDY8VAN3Aq4FfAy3A8939j2Ou80VglbtfM8n3lAIdBC1Y24BjwJXufvcp6nrc3f9hzLGbgRp3f3H4frW7rx/3HaetLdnf2d1vncq1x183fP8OYJ+7v3my3+t0zOxDwHuA57r7k2Gr33nu/rOpXkNEMic/6gJEREJbEk/c3c3sKPDYmGODZtYK1BG0XBUBvzGzsf9KLQD2Jl6Y2Qrg4wQtZLUEvRYxYLG7/ykMa7eb2e+B3wM/cvfJWrMms2nc6ynVNo3fOZlrb+Fkh8dcI2lmdj1BS+Bz3H1HePgFQA1BK6iIZBkFPBHJFoPjXvskxxIhDeAlwP5TXOc24BDwtvDnEEHLXSGAu78xbP16IXAt8Akze5m73x5+foSg+3WsgnGvu8e9nmptE70+1e+czLVPdY2kmNmHgbcDV7j7zvDYFQRdw8fM7K/D99qnc30RSQ8FPBGZjbYB/cASd79zohPMbC6wBniXu/8hPHYB4+577r4Z2Ax82sx+DfwdkAh4zQRjz8Y6j6e3xCVV2wyk6toDQN7pTjKzjwBvIejG3pU47u53m9kW4I3uvmcGdYhImijgicis4+6dZvY54HPh2LV7gDLgWcCIu98EtBKMV3uLmR0AFgCfJWjFw8yWEbTs/ZygdW85sJ5gYkHCncAXzexa4Mnw/EWcIuBNsbZ0/t5TsRe4KBxH1wUcd/eRsSeELXfvJWjZ7Daz+vCtNnfvA5Zy6qArIhFSwBOR2eojQBPwQYJQ1gE8CnwGwN1HzOxVwJeBx4GdwAeAH4ef7wFWAT8iGEvWBHwf+PSY7/gmQej7Zvj6BuD/hedPu7YZSsW1P0cwQ3cbUAws4+SxiwZcB1QQLv0yxlVm9iTQ6JqlJ5K1NItWRESSYmaXAP/k7q+IuhYRmZjWwRMRkWRtA5aEiyyfG3UxIvJ0asETERERyTFqwRMRERHJMQp4IiIiIjlGAU9EREQkxyjgiYiIiOQYBTwRERGRHKOAJyIiIpJjFPBEREREcowCnoiIiEiO+f8BkiXjOenrTUYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.arange(0,19.901,0.001)\n",
    "probs = []\n",
    "for i in range(len(x)):\n",
    "    probs.append(range_weigh_particle(x[i],10))\n",
    "max_range_prob = range_weigh_particle(40,10)\n",
    "x = np.append(x,[19.9,20,20.1,20.1])\n",
    "probs = np.append(probs,[max_range_prob,max_range_prob,max_range_prob,0])\n",
    "\n",
    "plt.rcParams.update({'font.size': 14})\n",
    "plt.rcParams['figure.figsize'] = 10, 5\n",
    "\n",
    "plt.xlim(0,21)\n",
    "plt.ylim(0,0.36)\n",
    "plt.xlabel(\"measurement $z_t$\")\n",
    "plt.ylabel(\"probability $p(z_t|x_t,m)$\")\n",
    "plt.xticks([0,10,20],(0,\"$z_{expected}$\",\"$z_{max}$\"))\n",
    "plt.yticks((0,0.15,0.3))\n",
    "plt.plot(x,probs)\n",
    "#plt.savefig(\"measurementprob.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#the incoming measurement:\n",
    "measurement = 8\n",
    "#create some particles\n",
    "particles = []\n",
    "for i in range(100):\n",
    "    x = random.gauss(0,5)\n",
    "    y = random.gauss(-10,5)\n",
    "    phi = np.pi/2\n",
    "    particles.append([x,y,phi])\n",
    "\n",
    "    \n",
    "    \n",
    "weights = []\n",
    "for i in range(len(particles)):\n",
    "    #this is a very simplified ray casting, normally this involves actual geometry calculating the first occlusion in the sensors path\n",
    "    #in this case however all particles face in the same direction and there is only one wall\n",
    "    expected_measurement = -particles[i][1] \n",
    "    if expected_measurement < 0: #particle beyond (or in) wall = not possible\n",
    "        weights.append(0)\n",
    "    else:\n",
    "        weights.append(range_weigh_particle(measurement,expected_measurement))\n",
    "        \n",
    "particles = np.array(particles)\n",
    "weights = np.array(weights)\n",
    "weights /= sum(weights)\n",
    "new_particles = np.array(resample(particles,weights))\n",
    "\n",
    "#map for plotting with wall at 0\n",
    "map = np.zeros((30,30))\n",
    "for i in range(len(map)):\n",
    "    map[i][len(map)-1] = 1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmcAAAESCAYAAAC4tXL+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3df5wddX3v8fcnmyXZ8CMbBNsHe+WnPoJyEZDtA7zBC0QlIAUj/qRQpe0DrFYtFCOhoASLEok/bh9tbYXWBpVaUHAvNghKg1SpUhNDoCjRoCHcxR+A2QTIgpvN9/4xczazszNz5pwzc+bHeT0fj30kZ87MnO+cPeezn/n+NOecAAAAUA6zii4AAAAA9iA5AwAAKBGSMwAAgBIhOQMAACgRkjMAAIASITkDAAAoEZKzHJnZajPbnOH5DjOzb5nZDjNzZrY0q3OjXszsAjN7Z9HlQG+Liln+zweKLhuKZWaHmtkKMzu46LKUEclZtVwn6WWS3ibp1ZLuLbY4KLELJJGcoWhRMWupJJIzHCrpKkkkZxFmF10AtOTlkv7DOXdnpycyswHn3HgGZao13iegIzNilpll+gJm1i9pt3NuMtMTVwQxqp6oOesCM1tsZg+Y2fNm9t9mdnrEPq8zs/8ws+fMbLuZfdXMXuI/d6iZOUlHSfpDv3lgS+DY3zezH5jZuJk9ZWZfNLPfDZ3fmdlVZvZXZvaEpOcCz73UzG4xs9/45/hPMzspxXVtMbN/NLP3+v9/zszuMLMDzGzIzL5uZs+Y2WYze0cr1xzY52P+e/eMmf3KzNaY2ctD+7zcf62n/PL/3Mz+NvD8CjPbFfH63zazu8P7mdmxZvYdMxuXdFng+feZ2cNm9oKZ/cLMPmVmcwLPX+C/z6/2y/mcmW01swv8599tZj/zm3huM7P9Q+XZ18w+Y2aP+6+xycwuDO2z2n8/TzSz+81sp5n9yMzOCV6XpJMlvdYvjzOz1XG/R6BVZnaCmX3NzEb9z+DDZrbMzPr85yNjlv85fJekIwKfzW8Hzts0FgXizgfM7GeSnpc0LW4E9u3GdzLxvQjst8zMfmLe34Gnzew+M/tfgeedmV0ZOuYUf/tJof3ajuWWc9xu/O7N7I/NbJWZ/dq/3n8xswWN65J0j3/IdwKfhUOjfo89yTnHT04/klZL2iZpi6Q/kfT7kr4p6beSXhHY72xJk5L+RdKZ8poAfiRps6S9Jc2RdKJ/njX+/4/zjz1T0m5JI/7//0jSLyX9WNK8wGs4SU9IusMvx1v97YdIekrSDyS9XdIZkv6vvIB3bJPr2yLpcUn/7p/zAknbJd0u6X5JH5T0ekm3Stol6Yi01xzY758knSfpFElv9Ms/Jul3A/s86r/em/z93iXp7wPPr5C0K6L835Z0d2i/Sb8Mfy5psaTj/edWSXpB0jX+Nf2ZX45/DRx/gf8+bwpc+1f9389KSWslneXvt0PS6sCx/ZL+U9IvJL1H0uvkNQlNSvrTiM/Uj/zf9RJJd0qaaLy/kl4h6YeS/kveZ+XE4HvPDz+d/vjf17/0v/enSLrE/z6s9J+PjFmSjvAfPx74bL7CPyZVLPLP+YT/nT/H/07NjylnN76Tie+Fv8875cXAK/x9zpJ0taQzA/s4SVeGyn+Kv/2k0H5tx3LlHLflNVc6SY/Ji99LJL1XXhL5eX+f/fxtTtJFgc/CnKI/22X5KbwAdf6R94fUSXpDYNtekv6fpJsC2x6VNBI69lB5Sdx7A9seCQYPf9t6SRslWWDbq/3XfXdgm5P0M0mzQ8d/3i/PPoFtffKSu1uaXN8WeYlgMAn8pP9afxHYtsD/kl/W6jWHnu+TF/SflnSJv+0A//XOTijnCqVPzpykPw7td5gfkC4ObT/P37/xx+WCiGsf9K89/D59StLOwON3yvuD8arQa9zgHzsr9Jk6IbDPgX75PhR3bfzwk9ePJJPXReZSeTcOwVgUFbNWS9occZ5UsciPO9slDaYoW+7fyTTvhaS/lfTDJmVtJTlrO5Yr57itPcnZN0L7/U3o/Z1xbfzs+aFZM3/PS/pG44Fz7rfy7lBOlCQze5mkwyXdbGazGz/yvmSPSIptXjSzfeTdjd7i/E+7/xrfk3fX8r9Dh3zDORdu3lsi6d8kPR94bZN0d9JrB9znnNsZePwT/99vBcqzTdKv5Tc9tHLNZnaGeU2M2+QFiucl7S9pob/L0/KCzbV+E8YhKcrczNdDj18nrwvALaHyftN/Pvw+Ba99TN61h9+nTZIGzOwA//ESf9uDode4S9LvSHpp4NjfOOfuD7zGk/5r0LEWXWFmC/zmvi3y/jBPyPsDPyjv89qOVmLRd/3vVlq5fSdTvhfrJB1rZn9tZq8xs7ktlD1Kp7E817jtuyP0+L81/f1FApKz/D0VTJx8v5J0kP//F/v//ou8L3Xw52hJL0o496C8L98vI577pbwkJrwt7MWS3h3x2u9r8toN20KPf5uwvRGQUl2zmf2evET2KXlNlSdI+j15wWCuJPnv7eslPSjpryVt8ft8vClF2aPs9pOdoEZ5R0Nl/bW/Pfw+RV173PsUfE+O1Mz34ysRrxE+l+Q1uXYa8IG0VsurlfqkpNfK+15+xH+u3c9hK7EoKpYlyfM7uVrN34sbJb1fXm3RvZKeMrN/Dvdxa0GnsTy3uJ3wGi/4/xKnUmC0Zv4OMDMLJWi/I68fg+TV/EheNfh/RBy/I+HcY/KqhaPuVH9HXhV0UDhJbLz+v0v6TMRzUftnIe01L5VXU/bWxl2imZlCQcA5t1nSuWY2S9KwpA9L+oqZvcI59xP/HLPMbHbobvNF8hLlaadLKO8pCnS+DRiN2Naqp+XdpZ8f8/wjGbwG0DG/1udMeU1wwYE3r+vw1K3EorxiU1DT72Ta98KP/38n6e/8mqM3yWtG7dOeKW9ekNftJSjuBrmIWN7J3yq0iOQsf3Pldcy8Q5LMbC95nSobH+5N8pogFzrnPt3KiZ1zz5rZBklvMbNrGwmgmZ0gr91/ZYrT3CXpVZIe9JtcuyHtNc+T15QZDCxvkTQQtbNzbrek/zKzK+R1dD1SXnX9Vnk1jEfKq1qXmR0mr2k0nJxFuVte35ODnHNfTrF/O+6SN+DhaefczzM43wvyBpMAWZsjL6mYihf+jdEfpDw+rpa3iFiUpOl30szmq8X3wjn3lKQbzOwMSf8z8NRWeYN5gs5osbx5vn9t/62KQU1aApKz/I1J+qyZfVReIvABebVaH5O8Oyoze7+k28xsQNJt/jEHSTpV0p3OuVsTzv8ReX2kbjWzf5RX9fxxeV+kL6Yo30fkje65x8w+K6/J8ABJx0uadM59uMXrbaqFa75L0sWS/tHMviQvkH1Qe+7gZGavlHeneLO8msK58t7j7fJGHklen78d8gLiR+U1B18ur7k0TXk3m9kn/eOPkvRdeR3wD5V31/znzrnH2ns3pnxJ3ojee/zXelhecnWkpBOdc29p8Xw/lnSh37z7uLzm9S0dlhGQc267md0vaZmZ/VLed+s9Sn8z0PhsXiDvc77DObdJBcSiJpp+J9O+F2Z2vbyY9D158euV8vqI/U1gt5slLTezD0naIOl0eSPG08r1/cvgb1XYT+TF0T8xs53yEtyyJOaFIznL39OS/lTSp+XV1GyW9Ebn3MONHZxzXzezxZKulNc3YS95TWX3yq/pieOcW2Nmb5Q30/JtknbKq6X7YKjDZ9zxj/l9uz4qr5p9f3l9qdbLq4bPRZprds7daWZ/IS9Be4e8gPUWf/+GX/rHfUjSkLxmxx9Ier1z7lf+ebaZ2dnykrhb5Y10WqbAHGYpynuZmW2SN/z7UnmBZIu8aSyeTjg07fknzOw0eUPtL5bXuX9MXpL9r22c8jp5d+FflPeH4kZ5/WKALPyBpH+QdL2kZyR9Qd7n9Aspjv0nSYvk9RHdT953/pSiYlGcFr6Tad6L++Qlen8kaR95tWTXyb9J931cXjPmpfJqJ78q70YzPEAprry5v3+d/K2KONfT5i3jdZm8KTlmyRsZvyWLslZdY5gvAAAASoDRmgAAACWSeXJmZu/w56XaYdFL5pzuT3Uwbt5SRqdlXQYAaAfxC0AZ5FFztk3SZ+W10U9jZofL6xd1raT5/r9fM9bTAlAOxC8Ahcutz5l5C5ve7ZybHdh2taTFzrnXBLZ9x9/v6lwKAgAtIn4BKFK3+5wdI2/kSNAP/e0AUGbELwBd0e2pNPaVN9dL0Jiko6J2NrOL5K1Yr7333vv4I488Mt/SASiN9evXP+WcO7DocgQQvwCk1kkM63Zy9oy8vhpBg4pZ9sE5d728uWM0PDzs1q1bl2/pAJSGmXU6sW/WiF8AUuskhnW7WXOjvOUlgo7ztwNAmRG/AHRFHlNp9PmLwe7lP57r/5i8GZOHzexcM+s3s3PlLS1xY8IpAaAriF8AyiCPZs0/lPTPgcfj/r+HOeceNbNz5C0t8Xl5y+i8iTX/AJQE8QtA4TJPzpxzqyWtTnj+TnnrEQJAqRC/AJQByzcBAACUCMkZAABAiZCcAQAAlAjJGQAAQImQnAEAAJQIyRkAAECJkJwBAACUCMkZAABAiZCcAQAAlAjJGQAAQImQnAEAAJQIyRkAAECJkJwBAACUCMkZAABAiZCcAQAAlAjJGQAAQImQnAEAAJQIyRkAAECJkJwBAACUCMkZAABAiZCcAQAAlAjJGQAAQImQnAEAAJQIyRkAAECJkJwBAACUCMkZAABAiZCcAQAAlMjsoguA+hrZMKpVd23SE2PjOmhwQMuWLNTS44aKLhYAAKVGcoZcjGwY1eW3PaTxiUlJ0ujYuC6/7SFJIkEDACABzZrIxaq7Nk0lZg3jE5NaddemgkoEAEA1kJwhF0+Mjbe0HQAAeEjOkIuDBgda2g4AADwkZ8jFsiULNdDfN23bQH+fli1ZWFCJAACoBgYEIBeNTv+M1gQAoDUkZ8jN0uOGSMYAAGgRzZoAAAAlQnIGAABQIiRnAAAAJUJyBgAAUCIkZwAAACVCcgYAAFAiJGcAAAAlwjxnaMvIhlEmmAUAIAckZ2jZyIZRXX7bQxqfmJQkjY6N6/LbHpIkEjQAADpEsyZatuquTVOJWcP4xKRW3bWpoBIBAFAfJGdo2RNj4y1tBwAA6ZGcoWUHDQ60tB0AAKRHcoaWLVuyUAP9fdO2DfT3admShQWVCACA+mBAQEmVeTRkoxxlLR8AAFVGclZCVRgNufS4odKUBQCAOiE5K6Gk0ZBJCVGZa9sAAEA6JGcl1M5oyCrUtgEAgOYYEFBC7YyGrPPcYyMbRrVo5VodtnyNFq1cq5ENo0UXCQCA3JCclVA7oyHrOvdYo0ZwdGxcTntqBEnQAAB1Zc657r6gWZ+klZIukDRX0jclvds591ST47pb0ILNe/nJWnDyu9S33wGa3PGUtt17o3b++N7Y/Yf+9POaPf/FM7bv2v5rjf7DH+dZ1FzV9bqQjnPOii5DWDsxrNfiFwBJ0nrn3HA7BxbR52y5pDdKOkHS05I+L+mLks4ooCyltfPH9yYmY2Hb7r1RLzrj/ZrVP3dq2+6J57Xt3htTHd9qMtgtffsd0NJ2oAuIYQByVURydpGkjzrnfiZJZvYhSZvN7FDn3JYCylMLjUSqnQRr3stPnpbYzZ7/Yr3ojPdPO29RJnc8FVlzNrkjsaIVyBMxDECuupqcmdl8SQdLWt/Y5px71Mx2SHqlpC2h/S+SFwiRQqu1bQ0LTn7XtBo3SZrVP1cLTn5X4clZpzWCQJZaiWHELwDt6nbN2X7+v9tD28cCz01xzl0v6XpJGh4eduvWrcu3dD3qsOVrFNUhpn/+i9XtPolRZs7fdqyWfurNRRcrFvPNZWN4uK2uGnlLHcOIX6gi4ld2zNrvMtvt5OwZ/9/5oe2DknZ0uSzwHTQ4oNGIUZ1lWci8SqsRMN9c7RHDUFvEr/Lo6lQazrkxSVslvaqxzcwOl3fH+WA3y4I9qrCQeVXmOqvzfHMghqHeiF/lUcSAgOslXWZm98gb6fQJSXfRkbY4ZV/IvEp3c3Wdbw7TEMNQS8Sv8igiOVspaYGkH0iaI+lbks4voBwI6KTpMO8+Cu2uNVqEsjcRIxPEMNQS8as8ur5CgHNu0jn3QefcAc65fZ1z5zSbgBbl1Y0Z/Kt0N1eFJmJ0hhiGuiJ+lQcLnyOVuNqxbtRqVeluruxNxAAQh/hVHiRnaCqpz1c3arWWLVk47fWlct/NVWl0KQAEEb/KgeQMTSXVjnWjViuLuznm7gEAVAXJGZpKqh37zNuP7UqtVqcDFqoy2hMAgK4PCED1xNWCHTQ4oKXHDenac47W0OCATNLQ4ICuPefoUiU9zN0DAKgSas7QVLM+X1n2Ucij+bFKoz0BVBfdJ5AVkjM01a0RPHk1P1ZptCeAaqL7BLJEcoZUujGCJ69pOao22hNA9VRpsmyUH8kZSiOv5kfm7gGQN7pPIEskZyiNPJsfmbsHQJ7oPoEsMVoTpcHSIQCqiviFLFFzhlS6MQqpiOZHRlcByALxC1kiOUNT3RyFlEXzY9qAxegqAFnqZvcJ4le90ayJpqo0iWsjYI2OjctpT8Aa2TA6Y9+467r0lo2R+wNAWVQpLqN1JGdoqkqjkFoJWHHln3QuNqEDgDKoUlxG60jO0FTS8k1FGNkwqkUr1+qw5Wu0aOXaaUlUKwErqfzcgQIos7LFZWSL5AxNlWkUUrNmy1YCVtR1BXEHCqCsyhSXkT2SMzRVpsXNmzVbthKwGtfVZxb5WtyBAiirMsVlZI/RmkilLJO4Nmu2bHU4e2M7yzsBqJqyxGVkj+QMHevmXDtpZuFuNWCxvBMAoExIzpBaVBImqatz7eS1iDl3oACAsiA5QypxEx7OmT0rtg9YHskOtVwAgLojOUMqcR3xw9sa8hzpSC0XAKDOGK2JVFpNthjpCABAe6g5QypxHfEXzOvX8xO7e2akIwsNAwDyRs0ZUombP+yqs47qmbl2Wlm3EwCAdlFz1gOyqO1p1hG/jslYWNIEuL1w/UARrhx5SF++/3FNOqc+M517wkt0zdKjiy4WkCuSs5qLG2UptZ5Q5dURvypNhSw0DHTXlSMP6Uvf3zr1eNK5qcckaKgzmjVrrtlyR0WrUlMhCw0D3fXl+x9vaTtQF9Sc1Vy7tT3dqs3Ksqmw0zI3Oz6vCXABRJt0rqXtQTSHospIzmouzXJHYVk2hYbPG05+kpLHVpKtTsvc7PhGWcYnJtVnpknnNFTiJligDhrftajtSfJqDiXhQ7fQrFlzcaMsk2p70jSFjmwY1aKVa3XY8jVatHJt02bIuObL+QP9kfsPzutvqbmz0+bbpOODZZe8QN94D0nMgPyce8JLWtrekEdzaCPhaySLjYTvypGH2j4nEIfkrOaWHjeka885Wgvm7UmC5sxO/rU3awptp59YXPJjpsjk0Tm1lGylab5NSiiTji97vz2grq5ZerTOP/FgzQpUlA30z9LwIfsnHtdJc2icMvR/a/WmGNVFctYjnp/YPfX/sfGJxGSqWcf3dpKVuORn284Jze3f8zEcHOjXteccre3jEy2dp1mZmyWUScczShMozvAh+2vO7D03cOMTu5veDMY1ezZrDk2SR8LXiioNnkLnSM56QFwyteL2hyP3b9YU2k6yEpf8mLwEreGFXbsT94/b3qzMzRLKpOMZpQkUJ+67e+ktG2MTk3abQ5PkkfC1ghr83kJy1gPikqax8YnI4NZoCo2b9b+dZCUq+TFJ4XvORrBpta9cszI3SyiTjm+n3x6AbMR9dyedi605ajSHNhKnPjOdf+LBHXXezyPhawU1+L2F0Zo9IG7EpqTYKSuSJpxtZ0qJqBUG4sr0xNh40xUJwpqN7EwzajXumuPKIkmLVq4t/eS5QJUlxYqkaXeuWXp0piMpG+fKY7TmeTd8T/c9+pupx4uO2F83Xfjqafu0M/K+09dEcUjOaqyRsMQFNqm9u65WE6fgccF9Fq1cmxhs0q5IkGYajJ2/3TXjuFZqv8JlyWu6EQCe13/62/rpr59rul83a46yTvikmUmSJN336G903g3fm0qWrhx5SE9sn3md7dbgp3lNFIvkrKbCyUOcdu+6sljKqdNJXZOSz2BfjKj3YXCgXyvOPqrta2CdTSA/aRMzqbp9P6MSpKDGc+E52xrm9c/SxwNdN1oR97pJ5UF3kZzVVFTyEFZ0v6moGrhTjzxQq+7apEtufiCxRi5N8hk3DYYk7T1ndkdJFP0/gPykTcyKjmFhL718jXYFOtLONmnztWfO2K9ZYhYUN1XHC7scN4I1RnJWUp0uRZSUJJhUmj5SwRq4VpoK0ySfg/P6E/u1dSLr/h9AnaRNUtox5E9vU5YY1hC+Zkna5bzt4WtPm5gdecUdhU/hgWKQnJVQFv2Z4pKHocEB3bd8cXaFzVDSlB/h626WXPX3mZ59fmY/s4ZOkyjW2QSitZKktKOs8St8zc22p/H8ZPzBnUzhseiI/SMTxEVHJE/ui+5hKo0SymI+m1OPPFDhr27Zk4e0U36MbBjVrITANDQ4oL33mq2J3dGBLYv3odnUHUCvyiNJaXjZi/fu/CQFOHT5mpb2T5MkdTKFx00XvnrGazBas1yoOSuhTvszjWwY1a3rR6fNIWaS3nx8553485Rmyo9GrWJUlf5Af99UgnRYQjDMKonKYlAEgOleedWdsc996y9O6V5BcvLSy+NjUzBBSkroOp2zTRKJWMlRc1ZCnc5IH1Xz5iTd88iTnRYtV0m1WY3ENK6vWZ9ZqolyhwYHSKiAEtvxQnJf0rKanbKVMakGMW3ClPV0HigfkrMS6nRG+qqOJFx63NC0BdqDGslW3DXsdtNHLjGrP1CMuCQlbfJSVVkNeGiY2xf9hsVtR72QnJVQp/2ZqrwW5FVnHZWYVKW9NvqEAcXYfO2ZMxKxLEdr9opHPvaGGYnY3D7TIx97Q0ElQjfR56ykOunPVOWRhM1WH2jl2ugTBhSjk0Rsvzl9kU2b+83pi9i7XLasPDOyr9iWlXvej9kW3bQZVbNIIta7SM5qqN3llcoiKamq+rUBSPbg1afrlVfdOS1B229Onx68+vQCS5VeMBGLsvnaM3OdBw71YK4iE9kNDw+7devWFV0MtKHTCXXzVOay9brh4WGtW7euFh1siF9A7zGz9c654XaOpeYMuSrzAuFlLhsAoHcxIAC56mRC3ZENo1q0cq0OW75Gi1aunTYRbdFlAwAgL9ScIVftTusxsmFUy76ycWqW/9GxcS37ykZJ2dVqVXXKEQDlR5cJdCLTmjMz+4CZ3W9mO81sc8w+7zSzR/197jez47MsA8ql3Wk9Vtz+8IzllyZ2O624/eHCy4Z6In4hK40uE6Nj43La02Ui69p/1FfWzZpPSLpO0seinjSzkyT9vaT3SFog6VZJd5jZfhmXAwUJN0WeeuSBLU8GO7JhVGPjE5HPxW1vBxPVIoT4hUysuP1hukygI5kmZ865rzrnbpUUd3twoaTbnHPfdM69IGmVpBckvSnLcpRR3v2nyiDqbvHW9aN68/FDqSeDbZyj2etkIeuJanvhd1xnxC9kIenmki4TSKvbfc6OkbS68cA558xsg7+9tnplVGBcB/t/2/gLPXDVaW2fI+ySmx/Qusd+M219uXb7d2Q1UW2v/I57XE/GL7QmqXYsrssE/dMQlqrmzMxWm5lL+Lkm5evtK2l7aNuYpMhmATO7yMzWmdm6J58s96LdSXplVGDcXeHY+ETqWqQ0d5ZO0k3f3zp1zjL07+iV33EVEb/QTUkxLKrLRBniF8onbbPm+yQdmPDz8ZTneUbS/NC2QUk7onZ2zl3vnBt2zg0feOCBKV+ifHplVGBSR/q0SUrazvgucM5uJ0ZRzZe98juuKOIXuiYuhi2Y1x9ZG8aNHaKkatZ0zj0r6dkMXm+jpFc1HpiZSTpW0m0ZnLu0Dhoc0GjEH+m8RwV2WlXe6vHLlizUxTc/EPlc2iQlau3MOI1zdjMximu+nD/QH9nPhJGfxSN+VVMWTX1FNBfGrf971VlHRe7PjR2iZD2Vxmwzmyup33toc/3HDTdIOsfMXmtme0m6VNJcSV/Lshxl02xUYB4dyTutKm/n+KXHDWnBvP7I59ImKeFO+nHnC56zm1NixN3lmomRnxVH/CqPLJr6imoubHWgEVP6IErWU2lcKWlc0vWSDvf/P5X+O+e+K+m98oLcdklvk/QG51xks0BdJH1Z8wognVaVt3v8VWcd1XGSsvS4Id23fLF+vvJMbfjIaTr/xIMVXmAxeM40U2JklQDH9qvbOZHpyE8UgvhVElk09RXZXBiMYfctX9y0xaFb8QvVkeloTefcCkkrmuzzBUlfyPJ1qyBuVGBSAOnkD3unVeXtHt8oc1JTQqtNDcOH7K81D/5C23Z6zYaDA/1acfZRU8c0e80sR1ImNVFnNfITxSB+lUcWTX1VaS7sZvxCdbB8U8HyCiCd9nPr5PikJKXVQBPeX5Je2LW7pdfMMgGO609C8yWQnSz66RbV17cd3YpfqA4WPi9YXv0NOp39Pq/Z81ttamh1/7xHUmY9cS2AmbKIP3VZAaQqNYDIFjVnBcurJiZN82Inx7c7CirqTjZpeyuBKapW7uKbH5CZvLk3QtpNgGm+BPLVafzK6hxlUKUaQGSH5KxgeQaQTpOIuOM76QPRZ6ZJNzNT6jObOnfwvRic1z/V1ywoKjDFrS4Q8XKVvIMGekkWN0HdvpHKY+oOulL0JpKzEqhaTczVX49f1DcumWsErIg8SZI06Vxk0tc/y9TfZ5qY3HNkXGBqVs3fZ6bdzlX2DhpA5/Ka+yyvjvt1qQFEa0jOkNrIhlFd/fWHI2uypHRNjXGGBgcia74mdjsNDvRr7zmzmwamuOr/ht3O6ecrz0wsB2vcAfV15chDuun7W6duErMc+Zhnx/2q3cCjcyRnSCVNktVKU2NQoybskpjVBbaPT6RaOL3Z6gLN+mgwZB2or5ENo9MSs4asEig67iNLjNZEKmmSrHaaGhfM658a7RiXPA0mrBIQ1BhJGbWqQJo+GqxxB9TXqrs2xXaryCKBYqZ/ZImaM6TSLHgNDuxZ1DfYNDgrZgBAw0nA/PMAAA5NSURBVI7xXVP/X7ZkoZZ9deO0/mWStG3nhA5dvkZDKZoZG9X/7TRPtnPnSzMoUA1J3+NWEqi473wVO+4Tv8qL5AypJPXnGujv04qzvUV9w02DSYlZ4/lg0+GK2x+OXEBc2tPMuO6x3+ieR55MDChRfTSaBaJWh6zTDApUR9z32xRd6x8lzXc+TbLTTlJ05chD+vL9j2vSOfWZ6dwTXqJrlh6dqtztXguKY67JH8+yGB4eduvWrSu6GD0rrs9ZeCmlRSvXJnbKjzM0OKD7li/WYcvXxDY9NISnLeufZdpn7mxt2zkxNVVHuJYtqvwD/X3TJpCN2+fNxw9FJoNx19q4FnRmeHhY69atCy+rWknEr+JFfb9N0nknHpw6ycniOx9VjkYMG9s5EZmsXTnykL70/a0zznV+C2UPI37lz8zWO+eG2zmWmjOkkvauMKnpYKC/L7bfWuO4ZiMupZnzyU7sdlMjSBs1deG7wDQjqaKu8dQjD9St60cj7y7pAAxURxZTUmTxnY8bld6IYVE1WF++//HIc335/sfbTs6IX+VGcobU0gznjkuuGjVZl96yMbKps9F02GzEZSuCyVfaQBS+xkUr18YmdczcDVRLp1NSZPGdT5P8hG8c47qHNOs2koT4VW6M1kSmktazW3rckD71tmMS17sLrl0pec0OQa22cQVr5KI0C0RJSV1d1u4DkE4W3/m0yU8w9jRWUAmL254G8avcSM6QqWYLg4eTrz6zqbvEkQ2jU/vct3yxtqw8U595+7HTznXeiQfPCChJgjVy7QSipKSORdCB3hL+zg8O9Gtu/yxdcvMDWrRy7VQMSxIVi6IEY8+5J7wkcp+47WkQv8qNAQEoRCedc4MjneYP9Ou53+6aMf2GFN3hv9X+JmkGEiAfDAhAmXUSG5rFsKjzZD1aE/nrZEAAyRkKETdSyCR95u3HtpT4NALd6Nh47GjNTjAXUDFIzlBmC6/8hl7YtXvG9nZGOxJj6onRmqicuBGZTmp5KRXWnQPQTefd8L3IxExqb7QjMQxhJGfIXXjB9MGBfplJcZW2ZRrKzUSNAKTptVtJ7U2MdkQWSM6Qq5ENozOWZIpbAaDhoMGByGp+qbM5itqRZn40APUWNwl3lGVLFpYmfqG6SM6QmaiAtOquTZGd9eOYpFOPPHBGbdWyr2yUTFPn6lYNFhM1Aoi6SUtSlviF6mIqDWSicWc56lf5N4JPK0s5NUZr3vPIk5EzaIeTvEYNVp7anR8NQH2kvRlbdMT+sSsAFBG/UF0kZ8hEXPNf0iSJgwP90+bY+czbj9U1S49uqVYq7xosJmoEkOZmbNER++umC19dqviF6qJZE5mICzKTzqm/z2bcNfbPsmkLpgelWV8zuG+esliPD0C1RS0rFzenWZniF6qLmjNkIi7IDA0OaNVbjtGCef1T2wYH+rXqrcfEJjhRtVX9s0z9fdNr4bpVg9VYseDnK8/UfcsXk5gBPaaV2fTLFr9QTdScIRNxd5aNWqZW5y2TZtZWRW0jUQLQDWnjGPELWSA5Qyaybv6LC4RZBzNm5gaQtW7FL4kYVlckZ8hM1Wa5ZoJZAFVGDKsv+pyhZyVNMAsAZUcMqy9qzpCozlXmTDAL1Fud45dEDKszas4QK25i2ZENo0UXLRNMMAvUV93jl0QMqzOSM8Sqe5U5E8wC9VX3+CURw+qMZk3EqnuVORPMAvVV9/glEcPqjOQMseJmuq5TlXnVRpgCSKcX4pdEDKsrmjURiypzAFVF/EKVUXOGWFSZA6gq4heqjOQMU+KGnRPMAFQR8QtVRXIGScw0DQBAWdDnDJJ6Y9g5AABVQHIGSb0x7BwAgCogOYMkZpoGAKAsSM4gqfhh5yMbRrVo5VodtnyNFq1cW6slVgAAaAUDAiCp2GHnDEYAAGAPkjNMKWrYedJgBJIzAECvoVkThWMwAgAAe5CcoXAMRgAAYA+SMxSu6MEIAACUCX3OUDjWwAMAYA+SM5QCa+ABAOChWRMAAKBESM4AAABKhGZNlNrIhlH6ogGoJOIX2kVyhtJi5QAAVUX8Qicya9Y0szlm9jkz+6mZPWNmW81slZnNDe23zMxGzew5M7vbzA7PqgxpsIZjdSStHABkqSrxC9VB/EInsuxzNlvSU5LOkjQo6TWSFkv6RGMHMztP0jJ/nwMl/UjS7WbWN+NsOWjcyYyOjctpz50MCVo5sXIAuqj08QvVQvxCJzJLzpxzzznnrnDOPeKcm3TOPSbp85JOCex2kaTPOed+6JzbKekvJR0u6aSsypGEO5lqYeUAdEsV4heqhfiFTuQ9WvO1kh4MPD5G0vrGA+fcs5J+6m/PHXcy1cLKAShYqeIXqoX4hU6kSs7MbLWZuYSfayKOuVjeHeUVgc37Stoe2nVM0n4xr3uRma0zs3VPPvlkuitKwJ1MtSw9bkjXnnO0hgYHZJKGBgd07TlH05m2R6xfv359872aq0v8QrUQv9CJtKM13yfpgwnP7ww+MLNLJF0mabFzbmvgqWckzQ8dOyhpR9RJnXPXS7pekoaHh13KssZatmThtNEzEncyZcfKAchALeIXqof4hXalSs786vtn0+xrZh+W9G5JJzvnwp25Nkp6laQRf999JL3M35471nAEek9d4heA3pHpPGdmtkrS2+QFtkcjdrle0qfN7GuSHpF0jaSfS/puluVIwp0MgChViF8AekNmyZmZHSKv6eC3kjaaWeOpx5xzR0mSc+4mMxuStEZec8D3JJ3tnJuMOCUAdAXxC0CZZJac+UPPLcV+10m6LqvXBYJYLgXtIH4BKBOWb0JtsFwKAKAO8p7nDOgaJhkGANQByRlqg0mGAQB1QLMmauOgwQGNRiRiaScZpr8agKoiftULNWeojU6WS2n0VxsdG5fTnv5qIxtGcyotAGSD+FU/JGeojU6WS6G/GoCqIn7VD82aqJV2JxmmvxqAqiJ+1Q81Z4Di+6Wl7a8GAEUhftUPyRmgzvqrAUCRiF/1Q7MmoD2T1DLaCUDVEL/qh+QM8LXbXw0Aikb8qheaNQEAAEqE5AwAAKBESM4AAABKhOQMAACgREjOAAAASoTkDAAAoERIzgAAAEqE5AwAAKBESM4AAABKhOQMAACgREjOAAAASoTkDAAAoERIzgAAAEqE5AwAAKBESM4AAABKxJxzRZchFTN7RtKmostRkAMkPVV0IQrUy9ffy9e+0Dm3b9GFyALxq2c/w1JvX38vX7vUQQybnXVJcrTJOTdcdCGKYGbrevXapd6+/l6/9qLLkCHiV4/q5evv5WuXOothNGsCAACUCMkZAABAiVQpObu+6AIUqJevXert6+fa66FO19KqXr52qbevv5evXerg+iszIAAAAKAXVKnmDAAAoPZIzgAAAEqktMmZmc0xs8+Z2U/N7Bkz22pmq8xsbmi/ZWY2ambPmdndZnZ4UWXOmpl9wMzuN7OdZrY5Zp93mtmj/j73m9nx3S5nHsysz/99P+n//m81swOKLlcezOwdZvYdM9thZrsinj/dzB42s3Ez+28zO62IcubBzD7hX9sOM3vCzG4ws/1D+1TyM97rMayX45fUOzGM+JVP/CptciZvDranJJ0laVDSayQtlvSJxg5mdp6kZf4+B0r6kaTbzayv66XNxxOSrpP0sagnzewkSX8v6T2SFki6VdIdZrZf10qYn+WS3ijpBEn/w9/2xeKKk6ttkj4r6eLwE/4f6tskXStpvv/v18zs0C6WL0+Tks6X9CJJx8j7Xf9z48mKf8Z7PYb1cvySeieGEb/yiF/Oucr8SPozSRsDj++V9FeBx/tI2inp5KLLmvF1XyBpc8T2GyV9MfDYJG2V9K6iy5zBNT8m6U8Cj4+Q5CQdWnTZcrzmUyTtCm27WtJ3Qtu+I+mqosub03twpqTtgce1+oz3YgzrxfjlX09PxTDiV7bxq8w1Z1FeK+nBwONjJK1vPHDOPSvpp/72XhC+fidpgyp+/WY2X9LBmn5tj0raIemVRZWrINN+x74fquK/4wTNvuNV/4wTw/ao2+92CjFsCvGrzc94IcmZma02M5fwc03EMRdLOknSFYHN+0raHtp1TFKpq8Xbuf4Ylbz+FBrlr+O1taquv+MZzOzNki6U9OeBzaW8/l6OYcSvVIhhnjr/jqfJOn4Vtbbm+yR9MOH5ncEHZnaJpMskLXbObQ089Yy8duygQXl3J2XW0vUniLv+R9spVIk84/9bxd9t1qr6GW+Jmb1V0uckne2c+2HgqbJ+xns5hhG/miOGear4+W5ZHvGrkOTMr7p/Ns2+ZvZhSe+W1wdjU+jpjZJeJWnE33cfSS/zt5dWK9ffROP6JUlmZpKOldcBs7Kcc2NmtlXetT0gTXUs3U/Tq4x7wUZJp4a2HSfp3wsoSy7M7I8kfUrSWc65+0JPl/Iz3ssxjPjVHDFsCvGr3c940R3omnSuWyWvU+URMc+fJ+lX8n7ZA5L+j6SHJfUVXfaMrn+2pLnyqkof9f8/N/D8SfKC5Gsl7SXvbvZXkvYruuwZXPsVkjZJOkxeQPuKpDuLLldO19rn/25Pk7Sr8XuW13n0CHk1EedK6vf/fU416VQs6QOSnpb0ezHPV/oz3ssxrJfjl399PRHDiF/5xK/CLy7hog+RN7LlBf/iGj8Ph/b7kLwh2zvlZeORQbCKP5JW+O/BtJ/QPu+U9DNJ45L+S9LxRZc7o2vvk/RJeVMRPCPvTuOAosuV07VeEPV7bgQwSaf7f7DH/X9PK7rMGV67kzQR+o4/G9qnkp/xXo9hvRy//GvriRhG/MonfrG2JgAAQIlUbSoNAACAWiM5AwAAKBGSMwAAgBIhOQMAACgRkjMAAIASITkDAAAoEZIzAACAEiE5AwAAKBGSMwAAgBL5/+JsQu9ynG7eAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "plt.rcParams.update({'font.size': 13})\n",
    "plt.rcParams['figure.figsize'] = 10, 5\n",
    "\n",
    "#setting up the axes\n",
    "fig, ((ax1, ax2)) = plt.subplots(nrows=1,ncols=2)\n",
    "ax1.set_aspect('equal')\n",
    "ax2.set_aspect('equal')\n",
    "\n",
    "\n",
    "ax1.imshow(map.T,origin='lower',extent = (-20,20,-20,0.5),cmap='Greys')\n",
    "ax1.plot(particles[:,0],particles[:,1],'o')\n",
    "ax1.set_xlim(-20,20)\n",
    "ax1.set_ylim(-25,10)\n",
    "ax1.set_title(\"before measurement\")\n",
    "ax1.set_yticks((-20,-10,0,10))\n",
    "ax2.imshow(map.T,origin='lower',extent = (-20,20,-20,0.5),cmap='Greys')\n",
    "ax2.plot(new_particles[:,0],new_particles[:,1],'o')\n",
    "ax2.set_xlim(-20,20)\n",
    "ax2.set_ylim(-25,10)\n",
    "ax2.set_title(\"after measurement\")\n",
    "ax2.set_yticks((-20,-10,0,10))\n",
    "#plt.savefig(\"locationbasedmeasure.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mapping with known localization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "feature mapping:\n",
    "\n",
    "example of 2 consecutive measurements of the same feature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import random\n",
    "import math\n",
    "from matplotlib import patches\n",
    "from scipy.stats import multivariate_normal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [],
   "source": [
    "#this is one kalman filter update step, though since a feature will perform actions we only need measurement updates\n",
    "#all inputs need to be numpy arrays! since alot of np matrix multiplication are used\n",
    "def update_feature(pose,measure,feature_posmean,feature_posvariance):\n",
    "    #range and bearing conditionally independent --> var = [[var_rad,0],[0,var_phi]]\n",
    "    #radial var =0.09 (std=0.3); angle(phi) var = 0.01 (std=0.1)\n",
    "    measure_variance = np.array([[0.09,0],[0,0.01]]) \n",
    "    \n",
    "    if type(feature_posmean) == type(None):\n",
    "        #the feature has not been observed before, initialize it with the standard errors\n",
    "        new_featuremean = np.array([pose[0]+measure[0]*np.cos(measure[1]+pose[2])\n",
    "                                    ,pose[1]+measure[0]*np.sin(measure[1]+pose[2])])\n",
    "        \n",
    "        #the jacobi matrix for previously unobserved features\n",
    "        dr_dmx = (new_featuremean[0]-pose[0])/measure[0]\n",
    "        dr_dmy = (new_featuremean[1]-pose[1])/measure[0]\n",
    "        dphi_dmx = -(new_featuremean[1]-pose[1])/(measure[0])**2\n",
    "        dphi_dmy = (new_featuremean[0]-pose[0])/(measure[0])**2\n",
    "        jacobi = np.array([[dr_dmx,dr_dmy],[dphi_dmx,dphi_dmy]])\n",
    "        new_featurevar = np.linalg.inv(jacobi)@measure_variance@np.linalg.inv(jacobi).T\n",
    "    else: #feature has been observed before update it using kalman filter measurement step   \n",
    "        #the measurements one would expected given the current feature position mean and robot pose\n",
    "        expected_r = math.sqrt((feature_posmean[0]-pose[0])**2+(feature_posmean[1]-pose[1])**2)\n",
    "        expected_phi = math.atan2(feature_posmean[1]-pose[1],feature_posmean[0]-pose[0])-pose[2]\n",
    "        expected_measure = np.array([expected_r,expected_phi])\n",
    "        #the jacobi matrix for kalman filter [[dr/dmx,dr/dmy],[dphi/dmx,dphi/dmy]]\n",
    "        dr_dmx = (feature_posmean[0]-pose[0])/expected_r\n",
    "        dr_dmy = (feature_posmean[1]-pose[1])/expected_r\n",
    "        dphi_dmx = -(feature_posmean[1]-pose[1])/(expected_r)**2\n",
    "        dphi_dmy = (feature_posmean[0]-pose[0])/(expected_r)**2\n",
    "        jacobi = np.array([[dr_dmx,dr_dmy],[dphi_dmx,dphi_dmy]])\n",
    "\n",
    "        #calc the kalman gain and apply measurement update\n",
    "        kalman_gain = feature_posvariance@jacobi.T@np.linalg.inv(jacobi@feature_posvariance@jacobi.T+measure_variance)\n",
    "        new_featuremean = feature_posmean + kalman_gain@(measure-expected_measure)\n",
    "        new_featurevar = (np.identity(2)-kalman_gain@jacobi)@feature_posvariance\n",
    "    \n",
    "    return new_featuremean,new_featurevar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [],
   "source": [
    "#feature at x,y = 2,2\n",
    "#definine 2 poses and measurements\n",
    "pos1 = np.array([0,0,0]) #first pose\n",
    "meas1 = np.array([3,np.pi/4]) #measurement at first position leads to feature at (2.12,2.12)\n",
    "pos2 = np.array([5,0,0]) # second pose\n",
    "meas2 = np.array([3.5,2.55]) #measurement at second position leads to feature at (2.01,1.95)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [
    {