0_introduction_to_numerical_robotics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to numerical robotics\n",
    "\n",
    "This notebook is a very general introduction to Pinocchio. It presents the main method to manipulate the geometry model of a manipulator robot: set the configuration, compute the position of an effector, check the collision or the distance to the obstacle. The main idea is to give a brief introduction of the general topic: how to discover and learn a robot movement constrained by the environment, using iterative optimization methods.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NB: as for all the tutorials, a magic command %do_not_load is introduced to hide the solutions to some questions. Change it for %load if you want to see (and execute) the solution.\n"
     ]
    }
   ],
   "source": [
    "import magic_donotload"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up\n",
    "\n",
    "For this class, we need the UR5 robot model (urdf, in the robotpkg_example_robot_data), the pinocchio python software, the python optimizers from scipy and for the display of the results: the gepetto viewer and the python matplotlib. We have that with this set of imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pinocchio as pin\n",
    "import time\n",
    "import numpy as np\n",
    "from numpy.linalg import inv,norm,pinv,svd,eig\n",
    "import matplotlib.pylab as plt\n",
    "from scipy.optimize import fmin_bfgs, fmin_slsqp\n",
    "from tp0.load_environment import createRobotWithObstacles,Target"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first load the robot model and display it. You should first start Gepetto Viewer: in a shell (CTRL-ALT-T to open a terminal), run the command: gepetto-gui. A new GUI widow will open (with the logo of our team): this is the viewer. You are not suppose to close this window during the rest of the class, and if you do close it, remember to reopen it before trying to display anything from Python."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this class, I wrapped the methods to load the robot model and create the obstacle field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "robot = createRobotWithObstacles()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next few lines initialize a 3D viewer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "You can open the visualizer by visiting the following URL:\n",
      "http://127.0.0.1:7000/static/\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "            <div style=\"height: 400px; width: 100%; overflow-x: auto; overflow-y: hidden; resize: both\">\n",
       "            <iframe src=\"http://127.0.0.1:7000/static/\" style=\"width: 100%; height: 100%; border: none\"></iframe>\n",
       "            </div>\n",
       "            "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Viewer = pin.visualize.GepettoVisualizer\n",
    "Viewer = pin.visualize.MeshcatVisualizer\n",
    "viz = Viewer(robot.model, robot.collision_model, robot.visual_model)\n",
    "viz.initViewer(loadModel=True)\n",
    "viz.display(robot.q0)\n",
    "hasattr(viz.viewer, 'jupyter_cell') and viz.viewer.jupyter_cell()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The scene should not be displayed in the viewer (check it). The robot and the red obstacles are encoded in the robot object (we will not need to see in depth what is inside this object). The object Target is the green dot that the robot should reach. You can change the target position by editing target.position, and display the new position with target.display().\n",
    "\n",
    "You can display a new configuration of the robot with viz.display (take a numpy.array of dimension 6 in input):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "viz.display(np.array([0,1.,-1.5,0,0,0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also set up a target with is visualized as a green dot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "target = Target(viz,position = np.array([.5,.5]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the robot model\n",
    "The robot is a 6 degrees-of-freedom (dof) manipulator. To make the example simple, we will only use the joints 1 and 2. We will then typically represent a configuration of the robot by a vector q2 of dimension 2. Yet the robot model is dimension 6, so let's define 2 helper functions to pass from 2 to 6 and from 6 to 2.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load -r 17-27 tp0/solution.py\n",
    "\n",
    "def q2_to_q6(q2):\n",
    "     '''\n",
    "     Transform a vector 2 into a vector 6, corresponding to locking 4 joints of the 6-dof arm. \n",
    "     '''\n",
    "     q6 = np.zeros(6)\n",
    "     q6.flat[[1,2]] = q2\n",
    "     return q6\n",
    "\n",
    "def q6_to_q2(q6):\n",
    "     return q6[1:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function computes the position of the end effector (in 2d):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load -r 28-32 tp0/solution.py\n",
    "\n",
    "def endef(q):\n",
    "     '''Return the 2d position of the end effector.'''\n",
    "     pin.forwardKinematics(robot.model,robot.data,q)\n",
    "     return robot.data.oMi[-1].translation[[0,2]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next function computes the distance between the end effector and the target."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load -r 33-37 tp0/solution.py\n",
    "\n",
    "def dist(q):\n",
    "     '''Return the distance between the end effector end the target (2d).'''\n",
    "     return norm(endef(q)-target.position)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last function checks if the robot is in collision, and returns True if a collision is detected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load -r 37-42 tp0/solution.py\n",
    "\n",
    "def coll(q):\n",
    "     '''Return true if in collision, false otherwise.'''\n",
    "     pin.updateGeometryPlacements(robot.model,robot.data,robot.collision_model,robot.collision_data,q)\n",
    "     return pin.computeCollisions(robot.collision_model,robot.collision_data,False)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random search of a valid configuration\n",
    "The free space is difficult to represent explicitely. We can sample the configuration space until a free configuration is found:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qrand(check=False):\n",
    "    '''\n",
    "    Return a random configuration. If check is True, this\n",
    "    configuration is not in collision\n",
    "    '''\n",
    "    while True:\n",
    "        q = q2_to_q6(np.random.rand(2)*6-3)  # sample between -3 and +3.\n",
    "        if not check or not coll(q): return q\n",
    "\n",
    "viz.display(qrand(check=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now find a valid configuration that is arbitrarily close to the target: sample until dist is small enough and coll is false (you may want to display the random trials inside the loop)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load -r 63-74 tp0/solution.py\n",
    "\n",
    "# Sample a random free configuration where dist is small enough.\n",
    "def qrandTarget(threshold=5e-2, display=True):\n",
    "     while True:\n",
    "          q = qrand()\n",
    "          if display:\n",
    "               viz.display(q)\n",
    "               time.sleep(1e-3)\n",
    "          if not coll(q) and dist(q)<threshold:\n",
    "               return q\n",
    "viz.display(qrandTarget())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## From a random configuration to the target\n",
    "Let' s now start from a random configuration. How can we find a path that bring the robot toward the target without touching the obstacles. Any idea?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load -r 82-89 tp0/solution.py\n",
    "q = qrand(check=True)\n",
    "for i in range(100):\n",
    "     dq = qrand()*.1                           # Choose a random step ...\n",
    "     qtry = q+dq                               # ... apply\n",
    "     if dist(q)>dist(q+dq) and not coll(q+dq): # If distance decrease without collision ...\n",
    "          q = q+dq                             # ... keep the step\n",
    "          viz.display(q)                       # ... display it\n",
    "          time.sleep(5e-3)                     # ... and sleep for a short while"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Configuration space\n",
    "Let's try to have a better look of the configuration space. In this case, it is easy, as it is dimension 2: we can sample it exhaustively and plot it in 2d. For that, let's introduce another function to compute the distance to collision:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load -r 51-59 tp0/solution.py\n",
    "\n",
    "def collisionDistance(q):\n",
    "     '''Return the minimal distance between robot and environment. '''\n",
    "     threshold = 1e-2\n",
    "     pin.updateGeometryPlacements(robot.model,robot.data,robot.collision_model,robot.collision_data,q)\n",
    "     if pin.computeCollisions(robot.collision_model,robot.collision_data,False): return -threshold\n",
    "     idx = pin.computeDistances(robot.collision_model,robot.collision_data)\n",
    "     return pin.computeDistance(robot.collision_model,robot.collision_data,idx).min_distance - threshold\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's sample the configuration space and plot the distance-to-target and the distance-to-obstacle field (I put 500 samples to spare your CPU, but we need at least 10x more for obtaining a good picture)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEICAYAAABLdt/UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzsnXd4VFX6xz/vTCaT3gshlNA70kHAgoCCUmxrZxV17evPteuqa1u7rl3U1cWOBbFhRelNeg29hARCQnrPlPP7415CEtImmZJyP88zT+bce8p7Zybfe+4p7ytKKQwMDAwMWg8mXxtgYGBgYOBeDGE3MDAwaGUYwm5gYGDQyjCE3cDAwKCVYQi7gYGBQSvDEHYDAwODVoYh7AYuISKzRORhX9vRXBCR2SLypK/tMDCojCHsBhWIyAERKRGRAhHJFZEVInKTiFT8TpRSNymlnmhgXRM8a7FriMiZIpLahPLXiMgyd9pUrf5FInK9p+pvQPuPisjHvmrfwH0Ywm5QnalKqVCgM/AMcB/wnm9NMmgIIuLnaxsMmglKKeNlvFBKARwAJlQ7NgJwAv319GzgSf19DPADkAtkA0vROgsf6WVKgELgXj3/l0A6kAcsAfpVamc28AYwHygAVgPdKp3vB/ymt3MUeFA/bgLuB/YCWcAXQFQN1xas2+PUbSoE2gNW4GXgsP56GbDWUL4PUAo49LK5DbS7dyW7dwKX1PLZ/1uvu1Sv/3X9+CvAISAfWAecVqnMo8BXwMf6+euBQOADIAdIBu4FUiuVaQ/MBTKB/cDt+vFJQDlg09vf5Ovfo/Fqwv+yrw0wXs3nVZOw68dTgJv195WF/WlgFmDRX6cBUltdwLVAaCUx3Vjp3Gxd/EYAfsAnwBz9XChwBLgLCNDTI/VzdwCrgA56vW8Dn9VyfWdWFjn92ON6+TggFlgBPFFL+WuAZdWO1WV3sC7KM/VzQ4BjVLqhVatrEXB9tWNXAdF6+bvQbowB+rlHdSE+H+0GF4j2lLUYiNQ/k83Hr1nPsw54BPAHugL7gHMq1fexr3+HxqvpL2MoxqAhHAaiajhuAxKAzkopm1JqqdIVoiaUUu8rpQqUUmVoInKKiIRXyvK1UupPpZQdTSAH6cenAOlKqReVUqV6Hav1czcC/1RKpVaq92IXhiWuBB5XSmUopTKBx4AZDSzbELsPKKX+p5SyK6XWo/WWL25oxUqpj5VSWXr5F9FuXr0qZVmplPpGKeVUSpUAlwBPKaVylFKpwKuV8g4HYpVSjyulypVS+4B3gctcvF6DZo4h7AYNIRGtV1qd54E9wK8isk9E7q+tAhExi8gzIrJXRPLRevSgDeccJ73S+2IgRH/fEW2opSY6A/P0yd5ctOEHBxBfzzUdpz1wsFL6oH7MFWqzuzMw8rhtun1XAu0aWrGI3CUiySKSp5cPp+pndqhakfbVjlV+3xloX82eB2n4Z2XQQjAmWwzqRESGown7SatBlFIFaMMDd4lIP2ChiKxRSv0OVO+5XwFMByagiXo42jiwNMCMQ8DldZy7Vim1vAH11PQ0cRhN8Lbp6U76sYaWr4tDwGKl1MQG5q9Sv4ichjZ5PR7YppRyikj1z6y6TUfQhmC26+mO1ezZr5Tq0ZD2DVouRo/doEZEJExEpgBz0MZdt9SQZ4qIdBcRQZu8c+gv0CY4u1bKHgqUoU1wBgFPuWDOD0A7EblDRKwiEioiI/Vzs4B/i0hn3aZYEZleSz1Hgehqwz+fAQ/p5WLQxp9rW/J3FOggIv4u2N1TRGaIiEV/DReRPnXUX/0zs6NNdPqJyCNAWD1tfgE8ICKRIpII3Fbp3J9AvojcJyKB+lNUf/3mfbz9pMrLWw1aJsYXaFCd70WkAK1390/gJbTJv5roASxAW0WxEnhTKbVIP/c0mmDmisjdwIdowxxpaL3JVQ01SH8ymAhMRRv22A2M00+/AnyHNhxUoNc7spZ6dqAJ+T7drvbAk8BatEnGLcB6/VhN/IHWs08XkWMNtPtstDHsw7rtz6KNk9fEK2jzAzki8irwC/ATsAvtsyvl5KGX6jwOpKKteFmAtmqmTLfHgfYZDtLPHwP+i/b0BNqqJYAsEVlf3/UZNF+Or2AwMDBohYjIzcBlSqkzfG2LgfcweuwGBq0IEUkQkTEiYhKRXmhzIPN8bZeBdzEmTw0MWhf+aGv5u6BtHJsDvOlTiwy8jjEUY2BgYNDKMIZiDAwMDFoZPhmKiYmJUUlJSb5o2sDAwKDFsm7dumNKqdj68vlE2JOSkli7dq0vmjYwMDBosYjIwfpzGUMxBgZNxuF0+toEA4MqGKtiDAwaidOpePjTX/hx3Q6iw4J48ZopnNLFVTczBgbux+ixG7Qq8ktKeePnFbz0/RIOZeV6tK3565L5YW0yTqXIzCvioU9+8Wh7BgYNxeixtwHKbXY+WbSBzLwiJg/rxYCkBF+b5BGcTsXf3prL9tQMAL5bm8zX98wgKiTII+3lFJZUTReV1JLTwMC7GD32FkBGTiEL1u1id2pmo8o/9NEvvPLtMj5dtIHrXv6S3Wn1ujlxO0u27eOdX1azfm+jQ47WS0Z+YYWoA2QVFLM15ajH2jt7UM8qN41Lx57isbYMDFzB6LE3c/YdzuK65z4nv7gMs0l4bOYkJo/s7VIdS7buq3hfbnewelcKPRJj6ijhXr5cvpknv/gdAJMIr/xtGqf361pPKdeJDA4kIjiA3KJSAPxMJjrGhNdTqvG0iwzl83uuZHnyAeLCQxjTJ8ljbRkYuEKL7bF/sXgT5z70Hpc8+RGb9tXmPrtuNuxK5b3vV7Fs0776M/uIr5duIb+4DACHU/HRr+tcrqNLfFSdaU/z07odFe+dSvHLhl0eacdq8eON689nYOcEeibE8PRVk+gS59lrjQsP4YJR/Q1RN2hWtMge+7YD6TzzxR8c94Zwx6zvWPDMDZhNDb9PLd+8jztf+RanXsn9M8Zz0bjm9ygdHOBfLW1xuY4XrpvCk3N+JzOvkOmj+jGmb5KbrGsY7SJDq6QTqqXdycDOCXzyf0akN4O2TYsU9sPZ+VR2cZNbWEJxaTmhQQENruPX1TsrRB3gl9U7mqWwXzVxKKuTU9iy7wixEcHcc9m4+gtVIzEmnLduu9AD1jWMuy84g+zCEnakZjC8R0eumzjCZ7YY1ExeWSnr0tNICAmjT3S9GxtbNTaHg2OFxcSEBGExm31tTqNokcI+pHsi0WFBZOUXAzCqdyeXRB2gXXRYnenmQmiQldn3X0ZeUSmhgVZMpoZEkmteRIUEMetm391YDOomo6iQC77+hLTCAgR4/LQJzOg/qN5yrZGD2bnM/GguaXn5dIoM538zLqJDhOfmaTxFixT26LBgPrznMr5fnUxIgD8XnzbQ5TpmnjeClKM5rE0+RK/Ocfzj0obHIXA6FR9+sZJ1mw7SvWscN/31dKxW14dIXCE82LUbl4FBQ5m7axtphQWAFvT0jfWr2qywv7xwOWl5+QCk5OTx+uJVPDP9HB9b5TotUtgB2keHc+O5oxpdPsBq4embpzQ4/+ZtqWzbkUbvHgnsOZDJ+59qsZM3bUtFORV33Dih0bYYNG/S8vPZk51Fn9hY4oJDfG2O2wn0s9SZbkuU2uxV0iU2m48saRotVti9ydKVu3nkqW9wOhUi0L9fhyrnd+313FppA9+y4lAK13/7DaV2O2FWK59dfAl9Yj0zBu1Uiq/2bOFoSSGTO/Wke4R3lqRe1mcAv+zbzcrDhwj19+eJ09tuJ+W6U4eyav8him02gv39uXbUUF+b1CgMYW8Avy3chtOpTbQqBYWFpVXODx7QyRdmNYjyMjtbtxwiJCSAnr1b545TTzJrzRpK7VovLr+sjPfWr+OFcyZ5pK1/rvqFz3Zv0trduprvzvsr3cKjPdLWcexOJ18e3MApXSK5YdgQTo3vQoBf65SFgvIyjpUU0yE0DIup5knRYZ07MP+Wv7I7I4te8TG0C/PcCi5P0uRvUEQ6okWgbwc4gXeUUq80td7mREx01S+3b88ELp4+jPWbU+jeJZZLzx/uI8vqpqzMxt23fsSO7do6/ytnjuWav53pW6NaGFY/c7W050Tv+wPJFe+L7OUsTNvrcWF/bOOPfHFgPQAf71vDnDNm0j+y9TkyW3E4hb/9No9CWzl9o2KZc95lhFtrnrdqHx5G+/DmuZiiobhjg5IduEsp1QcYBdwqIn3dUG+NlFUbA/MG1141llHDuhIc5M/QQZ254ZozmHL2QB65ewpXXDQSs7l57vNatXx3hagDfPbBcsrLvf/5tURK7TZeWr0cTBAeYAUgKSKC20aM9FibicFh1dKeX43x+5GdFe9tTgdLj+71eJu+4InVCym0lQOwPTuT2dvW+9giz9Lk7odS6ghwRH9fICLJQCKwval1VyYzv5Bb3v2GHYcz6dU+ljevP5+4cO9MZP32+zbs5XYmnN6HG689k+Bgq1fabSr+/lW/XrOfyaVNXG2Z+37/le92aztmTSK8eu55nNOtu0fXNb92+jTuWv4jR4sLuahbfyZ37uWxto7TOSSKrLKiKunWiM3hqJIudzpqyQn787J5Z+saBLhp4Eg6hUZ42Dr349bnShFJAgYDq91ZL8BrP61gx2HNCdbOw5m88uNy/n2555chLVySzCtv/gbAug0HKSws45EHpnm8XXcwcnQPTj+rD0v+SMbPz8Ttd0/G7OdbYS8rt/Pq50tI3p/OwB6J3PaXsfj51SyWNoeDp79ayIodB+ieEMPjV5ztMU+N1VmRllLx3qkUe3OzsJg9K7Q9I2L5/ryrTzqebythQ/ZBYgNC6Rue6NY2Xxh+AQ+t/57DxflM6dCPczv0c2v9zYU7h47h9oU/YHM6SQwJY0afmpdz5peXccmPn5FRot3sFqXuZ8GF1xJk8a8xf3PFbcIuIiHAXOAOpVR+DedvAG4A6NTJ9cnG/JKyaunSWnK6lx0706umdx3xSrvuwGQSHn7yIjIz8gkM9Cck1LW18A67A5PZhIj7NkW9/fVyvlywEYCte9MJCfTn+vNPrTHvJ4s3MHflFgCO5BTw9Fd/8Pw1DV+i2hT6xsSxJOXAiXRsnFfarU52WRF/XfE2qcU5ANzVZxIzuo5xW/2JQRH8b+wMt9XXXDm3Sy/6x8STVpBPv5h4wvxrfurek5tVIeoAaUX5HMjPpW+09v07leLeZT/z7d7ttAsO5c1x0xgQ084r1+AKbum+iYgFTdQ/UUp9XVMepdQ7SqlhSqlhsY1YLnbJqQMrHoP9zCYuHe2d7f8D+1dd2njKgI5eadedxMaFuSzqbz4zn6kjHucvpz/Nn0vd57Rrd0pVl8G7D9XuivjQsaqBMlKz8txmR338Z+JkLujVlxHtO/D4GeOZ2KV7vWUOFebywJ8/cs+qH9id1zgXy9X56fDmClEHeHfPYrfU2xbpFBrBqe071SrqAJ1DIwit1DsP9w8gMeTE3MfcPVv5cvcWyp0OUgpyuXPJjx61ubG4Y1WMAO8ByUqpl5puUs2M7tWZL+68gu2pGfRJjKNHgnfW+I45tQf/vHcKy1fupkNiFDOuGO2Vdn3JmmW7+G6ONppWWFDKsw98xVdLH3BLz31E/06s3nYiHu/wvp1rzXvWgO58vXJrhU+fCQN7NLn9hhIVGMRLEyc3OH+p3cYVf3xCWpF281l4eDe/nXcjkdamDR0FmqtuFgoyt6whgZZGdGAQ/zv7Yv6zfjkmEe4aMrbK6pljJcVV8h8rLapeRbPAHUMxY4AZwBYR2agfe1Ap5fZbWfd2MXRv5z0/4seZeFY/Jp7VssYeHXYHGSnHiIgPJ9BFdwQFeVUjARUXlWG3O7BYmv5zuWryMIIC/Nm+L51TeiYy7fT+teYd0yeJd265iFW7UujWLppzh9bvhz63WBuii2ig76DdWVl8sGkDAX5+3Dh0OLHBwQ27kGqkFuVViDpAdlkJu/OOMSKuaXscpnYYxO/p21meuZsgsz8PD5jepPoM6md4fAc+nXxpjefOTerJW5tXk1eu/c4u79n8HAcCiKrsJtFLDBs2TK1du9br7XqTPbvTWbNqL4kdojh9XB+vtp2fVcA9Ex5j36aDhEYG88T3D9BvdMMn/gryirn9qnc4cigbgHMvHsbtDzX/CeM3/ljJ6wtXAXDLmSP5+/i6n64yi4qY9PEH5JRq/6Q9o6OZf8WMRq0cKrKVc+b3b5JVpvXoQvz8+X3KTcQGumflVnZZIcF+Vqzmtrvdv7lwqCCPhan7aB8cyoRO9Q/RuRMRWaeUGlZfvta5xczH7Niexp23fYStXFtS9dfrTmfGzNO81v7XL89n3yZtuKMgp4hZd33AayufanD50PAgXv34BlYv3UVoWCAjT/f8srumkpKdWyHqAG8uWs20QX3pHF37UrVtmRkVog6wKyuLjKIiEkJd320YbPHng3GX89LmxdidTm7rP8Ztog4QZW19PmpaKh1Dw/lrn8G+NqNODGH3AEsX7agQdYDff9nqVWG3ldnqTDeE0PAgJkxpOR7+qjtv0o7Vfd1dIyPxN5sp19c4xwYFExUY2Ggb+kbG898zLml0eQMDd2HsVvEAsXFVdxDGxnt3e/LUW84hun0kABZ/P656+GKvtu8LesRFM7HvicfiCX260zO+7vmYTuERvHneVIYmtGdMx0787/wLPOoywMDAWxhj7B7AYXfy8vM/smLZLhI7RPHAv84nob13d6/lZxewZ/1+ErrGk9A13qtt+wqnU/HngUMAjEjq2CKDktREqaOMhRkrsDsdnBk3ilCL+4ZlbM4yDhZtJcAcTIcg14KkG3ifho6xG8LuA9au3c8fC7YRExPCFVeNIaARcUwN2gYO5eRfW19kZ4EWcL19YDzPDLyfQHPTA6+UO0v5YN8DHCnV/MOMibmI8e1O3vlq0HwwJk+bKcnb03jg3jkVboBTUrJ49PGLXKrD4XAyb9YCUnamM3xCP06b1jJ9RhvUT0bpsQpRBzhccpR9hSn0C+/Z5Lp3F6ytEHWA5ce+5oy4K/AzGR2Nlo4h7F5m48aUClEHWL/ugMt1vP/Y13z91gIAfvtsBQ9/YGb0uS1notOg4YRagvE3WSh3ahPBJkxE+btnWM8iVXdg+okFkxjTbq0B41v0Mt27x9eZbgjrFydXSW9YlFxLTgN3cCg/l1WHD1FQXlZ/ZjcT4hfMP3peT6w1ikhLODd1u5KEQPf4rekROpT+4VqsX7P4MSXxVkziOe+VBt7D6LG7SFpOHlvSjtIjPppusa4HQRg+oit33DmJBb9uJSY2lFtvm+hyHV36JnJge9qJdDV/Ngbu46d9O7n9jxNeAb+efiXxXo57OixqIMOiXA/YXh8iJi7seBdnJ1yLRaxYzZ7znJlXuokDef9FxEK3yL8TbOnisbYMDGF3iS1p6Vzzv68oLrdhMZt45bKpjOvV1eV6pk4bwtRpQxptx63PXY6fn5mU3ekMn9CfyTPGAlCUX8LaRcmERAQx9HRjhYM7eHHtcmxOJwBphfl8vH0jdw0f62Or3EuIX6RH6y+zZ7A+/XrsqhCAvNKNjO7wE2ZTy4hr0BIxhN0FPl29ieJybazT5nDy4Yr1jRL2phIcGsidr1VdvVCYV8w/zn+Z1L0ZAEy75jRudnFS1lXKy2xkHc4mun0U/tbWOeHmV829gKWZRstqzhTa9laIOkCp4whljqMEmZpvrOCWjvErdYEQa1XPesEBzcfT3pqFyRWiDvDDh8uweTAM3qGdaVzT4+/8tdttXNfn/ziy/6jH2vIlj5w6jhDdjWuf6Fiu7tf4J622SoilO35yYvgqwJxAgF/z82HemjB67C5w0xkjWXswjR3pmXSIDOPus93jJuDAkWze/341SimumTKCbomue7AMDqu6Fd4a5O/RaEmzH/mczNQsANIPZPLR419y7/9u81h7vmJ0YmdWXHkTx4qL6BgWXmt0e4PasfrFMiThfQ7kvodJ/OkWeSsmaT6dotaIIewuEB0SxLxbriKvpJSwAKtb/JMXlZRz83NfcixX8+v85/YUvnzqGsJcdLU74qy+TL7iVH76dCWBwVbuefkqTB6Mb3qSP5pS1/3R+JrC0jKS0zNpHx5KYmTtgaPD/K11BmcwqJ9w6wBOiX/Z12a0GQxhbwThgU3f9Xec1MzcClEHyM4vJiU9h/7dElyu6/ZnLuXGRy/Az+KH2cNjwZfcM52Nf2ylpLCU4PAgLr6r+bv1rUxGfiFXvPs5abn5WMwmnr94Muf0b/qmHwODytidpZQ6sgjyi8ck3pNbQ9h9TGJMOBEhgeQWasEtwoKtdIhv/AYUq5fG/fuP6c1721/m4PZUugzoRHSCZ1dWuJs5azaTlquF5rU5nLz2x8o2L+zF9qOY8CPAz/VlvM0FpRSFBS9TWrYQi18vwsIfxWRqXPCUppJTtoslh/9OqSObUEsnxiXOItDP9bCgjcEQdh2ncnKwKItgPytxAd7zxhgSZOWNey7m3W9X4HTC9dNHERHSeNex3iS2QzSxHVqmCFRf7VI93dZYk/Fv9uV/AwgDo2+lT2TL9BlTXPwRBQXPA2ArX4vCTmTkf3xiy+as1yh1aMFqCmwpbM+ZzdDYe7zStiHsgN3p4M51n7IscycmhHv6nselSaO81n7PTrE8/3cj5Jk3uXLUIH7bvocd6ZkEW/25b/IZvjbJZxwr3aKLOoBic9YbdAmbRoC5ZT2FAdjKt1VJ223basnpeZzKVi1d7rW2DWEHlmbsZFnmTgCcKF5M/pELOw03VkC0ADZsTmHz1lR69Yhn1PBuDS4XHhjAFzddTmpOPjEhQYQGtN3J0ZMFR6GU55bKNgaHIxt7yTxEgrAETkZMNQ9XWq1jKC7+qCLtb/XdZrI+kTPJKt2KQ5XhbwqnZ/hlXmvbEHbAWc11sVMpfOHO2MA1lq3czUNPzuP4V3X37ecwdVLDgwtbzGa6xLS8Xqm7iQkYRELQaI4UrwCge/jFXhsLbgj20sUUZ18DaJGuygrfJCR2PmI6ecg0MGgaCjtlpYvxs/QkJOQm7xpbiXZBI5nc6UsKbAeJ8O9JgF+U19o2hB04Pb4XQ6OSWJd9AIBbe03A39y4j+ZocQEf7tyAWYRreg8lKsBz/jfaOn8s2UHl++8fi5NdEnYDDZOYGZvwElmlWzCJheiAfr42qQKnbSelOTdjwoFTP6YcB7GXLcUSeF6NZYKCLiQo6ELvGVkHwZYEgi2ur3BrKm4RdhGZBLwCmIH/KqWecUe93sJi8uOtETNJzj9MqF8ASSGN660U2sq46OdPSC3MA+Cngzv5Yco1WBt5kzCom3bVQhDGxXo3BGFrwiRmYgObl+tn5SyiNPtKhCIQwaRUhbiLyXjSqosmK46ImIE3gIlAKrBGRL5TSm1vat3exM9kZkBExybVsTMns0LUAXbnZXGwIIeeEc3nsbY1MePyU0k7ksumrYfo2b0dt1w/ztcmGbgR5UhDOTMr0iICSmEJmoGfdbQPLWv+uKMrOQLYo5TaByAic4DpQIsSdnfQPjgMf5OZcqc2Fhjs509cYM0uXgvKy/hw20bK7DYu7T2QxFCjt+kqgQH+PPagsZqotSJ+HRFTAsp5RD8QSnDMj5gtjXceppQTHKlgCm3VvX53CHsicKhSOhUY6YZ6WxwJwWG8ecb5vLBhCWaTiQeHnkmE9eQ16U6lmDH/KzZmaD/YL3Zu5eeLryYioGWsXzfwLsszdpNSlMXImK50DXVPkI2WgEggAdFzsBW+ilJ2/ENuxtQkUbehcm6C8qWABcKfQgJbZ8fAHcJek8OUk5aUiMgNwA0AnTq1XnedEzp2Z0LH7nXmSS8qqBB1LV3Ixox0zuxkBB8wqMrsvct4eccvAASYLcw+9Xp6h7f3sVXew+SXhDXiJfdUVvqzLuoANlT+o61W2N2x3S4VqDw43QE4XD2TUuodpdQwpdSw2Ni2PeYcaQ2s4lTKLEIHYyimXpxOJ9/NXspbj8xl/ZKdvjbHK3xzaF3F+1KHjZ8Pb/GhNS2dao7qlL3VLmt2h7CvAXqISBcR8QcuA75zQ72tlkCLhf9OuoA+0bF0DY/khTMn0z2yZW7N9ybvPvEtbz0yl+9mL+Xhv85iw7JdvjbJ40RbQ+pMG7iAdRL4nVjKKaF31OuhVTmOoUp/RdlqjivsLPsTe+GbOMuW1njeVzR5KEYpZReR24Bf0JY7vq+U8t0+3hbCiIQO/HRxy/THAVBSXM4T93/JhjX76dI9jkefv5S4drW7vnUHa/44MR/vdCrWLkpm8NjW7bjrkQHTuXv9HFKKsjgjvjeXJbXJ6Su3IKYgiJ4Dtk0gkYilR535lT0FlX0pOLMAE4T9Gwk6EZXMUfob9pybAScOwC/8ecxBno1a1lDcssBaKfUj8KM76jJoGcyZvYy1q/YCsGdnOm+99Av/eu4Sj7bZoVscaftPLH/r2K31TyR2Donhy9NbXwATXyFiBf8RDcqrSr7SRR3AiSp6p4qwO0u+g4qV9eAo+bbZCHvbdmln0GhysouqpQtryek+7njuMkZPGkhSrwQuvW0i51zmPUdtBm0PkWq7xqWq+18xJ9SYLi+cRXHm2ZRm/xWnPdWjNtaGsSXSoFFMPG8gv/+0GZvNgQhMnu75WKARMaE8/M61Hm/HwACAoBlQthhs67Shm7B/VTltDrkdZT+Is/xPxDIAv7D7sZcuwFbwNAAO+07Kcm8jMOabmmr3KIawGzSKAYM78/oH17NlQwpduscxYHBnX5tk0IxxKju7cj8gr2wHMYHD6Rbu2WE7dyCmYCT6M5QzGyQMqRYBSUwhWKLernLMad9TLb3X43bWhCHsBo2mS/d4unSP97UZBi2A5OxZ7Mx9D4C0ogWYxEyXsOYxHl0fYmq4V0azdSy2An+gXE/7xs2FIewGbsPpVJhMTQ/wbdD6OFa6vmq6ZH2LEXZXMFv6ExA9B3vJD4g5Hkuwb4YODWE3aDJrNxzgied+oKCwlAumDObvN473tUkGzYxIa1+ySjdUSjcf18Duxuw/FLP/UJ/aYAi7QZN5/NnvycvXgnF/9e06hg9JcimakUHrp1/07Yj4VRpjv9yt9Zfb08ks/BKTKZC4kCswm9p2HIRhhfPvAAAgAElEQVRWJexK2bE5c7GYohAxVnJ6A7vDSUFhaZVjuXklPrLGoLliFn8GRN/hkbrtjly2pV9EuUPzZJJT/Ct94j+vd1dpa6bVqF9R+W5WpY5j5aHRrDt8PuV6dHADz+JnNnHeOQMr0vGxYYwa3tWHFhm0NQrK1leIupZeg82RWUeJ1k+r6bHvzX6GcsdRAIpsOziU9zbdoh7wsVVtg7tuO5uRQ7uSX1DC6JHdiAhv24/BBt7F6peI1kfVdoGaJRQ/s2fdWzR3Wo2wO1Rx1bSzuJacnufPrLWsz9lIXEAcU9tPxmKy+MwWbyAinDa6br8brYHVi3bw4au/YfYzcf3dkxk4wngyqYkfDv/Bl6k/4m+ycGPXKxgWNcCj7QX596JL1FMcznsdkymQpKjHMYm1/oKtGPGF28phw4aptWvXurXOY8W/sz3jdhQ2/ExhnNLuY0L8e9eav8SeQ4EtjXD/TljN7nOZuyl3Cy/sfKUifWbsaVzXteU6+2pu2MrtvH7nh2xctJ0u/Tty55vXERbleY+Hmel5XDfpRWzldgACg618+Pu9hIS1zOAoSin+SNlHga2Mszp1q+JGuikcKErlrk1PVaStJn/eH/4sAea2JbQOZyH5pX9iMccQYh1Yf4EGIiLrlFLD6svXanrsMUHjGZb4PcXlewm1DsDq167WvJkl2/kl7U7KnYX4m8LpEDIOs/jTL2IqUdamBbvYnr+jzrRB0/jy5R/5efZiANIPZDIrNIB7373R4+0eTcupEHWAkqIysjMLWqyw37fkF77Ypfl27x4RzbzpVxLqBnHPKc+rki5zllNkL2lTwm535LEt/UJK7fsA6BBxN4nht3rVhlYzeQoQZOlKTPDEOkUdYGP2bMqdhSgFBY5CtufNZ0vuPOal3E6hrWmTLklBVaNDdQ6uO1rUvNQVXLHiWW5Y8yq7CtKa1LY72ZKazvxNOziSW+BrU6pwZH9GlfThfRm15HQvXXq1Iy4hoiLdqVsc7TrWviOxqKiMt17+hccf+JJli5rXzb2wvLxC1AH25GaxLO2AW+ruFdqVhADN66ZTgVXa8dquH1mbvdst9bcEsovnV4g6wOG8171uQ6vpsbuCVLqfqUrvy5yFHC1NJsTS+AhPp8aMJNeWx7qcDcQHxHNFp9p9YmzNO8B/ds7TEiXwwKb/MXfsQ41u2118tXYL//p2AUpBWICVT264lG5xzSMQyOgpQ/jt42UVkW/GTPXORpDgkABe+ORGvv90JWaziekzRuPvX/u/z7P/mscqPRDI8kXJPP/m1Qx0gz+dgxk5PPnpAo7lFzH91H5cM3G4y3VY/cwE+lkosZ+IKBReQ2zexhDkF8jTA+5m6bG1fLx/NSlFWaQUbWRxxhb+O/J2uoUk1F9JC8dUzSukSbz/VNcmhX1w9HVklG6lxJ6LICg9RKtgItK/Yz2l62dywtlMTji73nyHi7OqpDPL8ihz2LCaGzfZuihtH6uPHqJ/VDznJdU+v1AfH6xYz/Gpl/zSMuat38bdk05vdH3u5NTzhvDvb+5i4+JkuvbvyLhLTvVa27Htwrn2zkkNyrt1U0rFe6Vg+5ZDbhH2e/87n11p2lPlK98so3v7GMb2c2340GIy858zz+WuxT9SYrczs98QRrd3XxziUEsIk9qdztPbfqo4ZlMOtuYebBPCHh08hazi+eSWLMAkgXSJfsbrNrRJYY8O6MHFSZ9TaD9CubOM1cfew+YsYVDkJU0eY3eFUyK7EeYXRL5dW8EzIqpno0X9x4M7uGXxtxXpzJIirunTuN5sWEDV8dDQgIBG1eMpho4fwNDxnl1p0VR69W3PutUnHsd79XFPAOqUzJyq6YwccFHYASZ16cnZST2wO534m81usa0yJjHRI7R9xfCiCaFnaKLb22mOiPjRK+5dbI5jmCUEk8n7/z9tUtgB/M3BRJm7AzC9Y81R0EsdJQSYm/YYtTJrIxtzttMxKIFzE87AVGlHbHxABG8Nv42fj6wlxC+QCzuMaXQ7v6RUHcP8OWVXo4X94anjufmjbziaX8jILh2ZMXpwo+1qqzzw+IX89/XfOZaZz1nnDGCwmzZtnTGgG7+s0wJ5B1j8GNm78U8BJhGPiPpxnj3lGt7YPZ9cWxHTE0fSJ7zpT8MtCYs5xmdtt5rlju6k2F7A+/ufJqV4N1H+cVzX9Z/EWl3vca3K2sizO96tSF+YeDYzkqa709QKXtq4lFc3r6hIX9p9IM+Ontzo+pRSlNjsBPm37jX4LQ2b3cHnSzaRmVvIpGG96NPJcJvclmhzyx3dycKMb0gp1nrA2eUZfJ/2Add2dX0X64acqpHNN+YmMwPXhL3QVsq9a79mbdZB+ke258VhFxNpDT4p3y0DTuVwUT6r9DH2B4c1zQ+0iBii3kh2bzzIx89/D8BV906lxynuC0Ji8TNz1Vmej1Zl0LIxhL0GShxV43mWVks3lM7BVSeKOgW5PnH0xo7FLDqqra5Ylbmf57f9xlNDzj8pX4DZjxfGnNcoOw3cR0FOEQ/+5WUKcrTfzPY/9/L+micJjTj5Zmxg4CmaJOwi8jwwFS1cyF5gplIq1x2G+ZJR0RPZmLuccmcpJkyMjW2cYE5qdzrZ5XlsyEmmU1AC13d1PRzY0dL8KumMkvxacho0B44cyKwQddCEPv3AMUIH1S7s5eV2tqzZT1CIlej4cA7uyaBLz3hi4ty3I9qgbdHUHvtvwANKKbuIPAs8ANzXdLN8S4egbtzZ8wUOFu8iPqAj7QOTGlWPSUxc1Xk6V3Vu/Lj6lA4D+DVtO059SebUju7bnmzgfhK7xREZF0ZOhnYDjowLo323uFrzl5fZuHfme+zYdAgAs78fdiAo2MrTb/2V3gM6eMNsg1ZGk4RdKfVrpeQq4OKmmdN8iLLGE2X1/cTUWQm9+fC0mWzITqFfRHtGxRqOp1zBVm7nz0U7MJtNDD+zN06Hkx9mLyHraD5nTB9Cj4HuW78NEBwWxHPf3sXnL/8MwKV3TCY4tPaVVeuW76kQdQBHuR0sZoqLyvjqoxU89FzzCvqcU1RCUVk5iZFhbdrfeXPHnWPs1wKfu7E+A50h0Z0YEu1eAWoLOOwOHr7+fTat0iLFjxjXh0A/WPzNOgB+mL2El+ffTVJv96wxP07HHgnc/cbMBuW1BlSdoK68Rs1qbV5TYF+v3cpj837H7nRyZp+uvHrVVMymVuWVpNVQ77ciIgtEZGsNr+mV8vwTsAOf1FHPDSKyVkTWZma2bSf4bZEt+4+wKvkg5TZ7/ZndxJ5thytEHeDPhcks/i0ZAgNAhLJSGxuW+NaPy+BTuzF+2iAAzGaT5qlShHaJkcy4yTcR7mvC4XTy5Ld/YHdqPs8XJe/j9+176yll4Cvq7RIopSbUdV5ErgamAONVHYvilVLvAO+Ato7dRTsNWjD/+XoJHy7QeskDktrx7j/+gtXi+d5oUEgtHgXNZgiwQkkpiV1rH//2BiLCPU//hevunIQ1wIJ/gB85xwqJignFz+K5zUOuohQ4nFX/bW12h4+sMaiPJj1HicgktMnSaUop30W2MGi2FJeWV4g6wJYD6azYdsArbXfsFsfV/zgHk0lOjAfrf01+ZmY+OI0RE/p7xZb6iIoNJTg0AIvFj7iEiGYl6qCFQLx1wgm/PAM7tmN8v+4+tKjl4VROSh22+jO6gaZ2m14HrMBv+j/OKqXUTU22yqDVYDab8DObsDucFcesdXhFdDeX3XwWF157Gil7M7jnyrcpLS4HYOqMMVxy60Sv2dEauGHcCM7s05X8klIGdmyHv1/zmgNwhfkHd/BLyi46h0Zw64DRBJg9ey3LMnby4KbPKbKXMS1xCI8MuNCjk89NXRVj3LIN6sRq8eOBy87i6c+08dnzRvTm1D7u24nZEPytFrr3TeTlL25l5e/biUuIYJw+rt2a2bw1la3b0+jVox1DG+BZMq+whIV/7ibAamHCqF74mU9+oO/Zznf+T9zForR93Lrkm4r04aJ8XhwzxaNtPrz5S4rsZQB8l7ae0+J6c1a7fh5rr+XeclsxR4sLKLbbSAqNbBVLyi4cM4CJQ3pSVm4nJtx3OzA794incw/fL2Gtj7JSG7Ne+pmd2w7Tb1BHbrjjHCwuDs0sX7Wbh5/4Bqc+Lv7g3edx9vjahaSouIy/PfoZKUc075EL/9zFs//wjF8jX7P6aEqV9KpqaXfjUM4KUT9Oob3Uo20aa5WaGf/dtoZRX77JuHnvcuvib3E2wknb1zu2c+ZH73HOp7NZkerZH21DCQ20+lTUWxKz3/qDH+etZ++udL77Yg2fvrfE5ToWLNxeIeoAv/6+rc7863ekVog6wOK1e8gtaJ3TZv2jK0VYc0J7azgZRYUea88sJi5PGl2R7hAUxelxfTzWHhg99mZFsa2cp9YtrFjL/OPBnSw9vJ8zEhu+KWlvTjb3/v4zDv2GcOP8b1k180aC/f09YLGBJzi4L7POdEOIiQ6tko6NCa0lp0ZUeNWoPwFWPwKtdf9mvl66hbmLNxMREsh9l4+jU3yky3b6gvM69+bY8InM3b2VHWk5rDmQzrhP3ueDqRcxLMG9PuP35B7jtS0rUUpxf+8LCbNaODW2B2EWz0ZVMnrszQgn6qQeusPFHvvhgvwqZQpt5eSUlrjFPgPvMOzU7tXS3VyuY+ZVYxg9shvBwVaGDu7MjdedUWf+ft0SuPEvYwiw+hEZFsjjt55X5yT3ht1p/PujBexIyWDV9oPc+eZ3LtvoS67uPZTugXGUO7Qlm0U2G/9evpgDuTn1lGw4BeVlXP7bHL7dv53vDiTzzJoVjIrp6XFRB6PH3qwIsVi5feDoCr/qYxI6c1r7JJfqOCW+HYmhYaQVaL5KBsUnkBBSd2/NoHlx4RWjCA6xsmNbGgMGdeKsya77BwoKsvLUoxe5VGbm+aOYef6oBuXdf6RqWMeD6Tk4nQqTqeXMCflV2zW74cgRJsyezaPjzuKqU05pcv0HC3LJLDnhEC6rtJj9+dkMinHvTueaMIS9mXHn4NOYktSbQns5A6MTTvrx1UeYNYC5F13O58lb8Debuar/IGPbdwvknGmDOWeaeyNXlZbZyMkvIS46pMm/iaG9OhLg70dpubaTeGTfTi1K1AH+PmwUy1NTtE6QAnGAUymeW7bULcLeKTSCKGsg2WXaE3OEfwBJod4ZrjIiKBkYtAE27Ujl3ue+oaCojF5d4njlob8QFtK0WJzbD6Tzw8pkIkICmHH2MAKtLS8wS6ndznsb1vHCsmUI2o0p3Gplwy23uqX+5JwMXtm0HKdS3D5wDP2jm7Yqy4igZFAvdqeT535dwvK9B+kVH8OjU8YT1swCVxu4h1c+WERBkbbkbuf+DL78aT3X/WV0PaXqpm9SO/omtas/YzMmwM+P6wcPZdmBg6xOTcViMvGvcWe5rf4+kXHMOvMCt9XXUAxhb8N8sHI9H67eAMDeY9lY/Mw8e/4kH1tl4AnKbY46020Zq58fH190MQdycwkPCCAmKKj+Qs0cY/C1DbPvWHaV9P5j7lsRYNC8uPaiUZj1naQxkSGcP7HpY8itCbPJRLeoqFYh6mD02Ns0Z/bsytyNJzaunNmjiw+tMfAkZ53ai55d4jickUfvrvGEhXh+yZ2B7zCEvQ0zsU933rp8Oiv3pdAzLoaLhzQPT4cGnqFDu0g6tGsZm4gMmoYh7G2ccT27Mq6nEW7PwKA1YYyxGxgYGLQyfLKOXUQygYONLB4DHHOjOb6ktVxLa7kOMK6lOdJargOafi2dlVKx9WXyibA3BRFZ25AF+i2B1nItreU6wLiW5khruQ7w3rUYQzEGBgYGrQxD2A0MDAxaGS1R2N/xtQFupLVcS2u5DjCupTnSWq4DvHQtLW6M3cAziMgsIE0p9YSvbWkOiMhsIFUp9ZCH23kU6K6UusqT7dRjwyLgY6XUf31lg4F7aYk9dgMXEZEDIlIiIgUikisiK0TkJhGp+P6VUjc1RNT1uiZ41mLXEJEzRSS1CeWvEZFl7rTJ04jIoyLysa/tMGieGMLedpiqlAoFOgPPAPcB7/nWJAMDA0/QIoVdRJ4Qkc0islFEfhURz4ck8QAi8ryI7NCvZZ6IRHi6TaVUnlLqO+BS4GoR6a/bMltEntTfx4jID3rvPltEloqISUQ+AjoB34tIoYjcq+f/UkRyRMQhIkpELql0jbNF5A0Rma8/MawWkW6VzvcTkd/0do6KyIP6cZOI3C8ie0UkS0S+EJGo6tcjIsHAT0B73aZCEWkvIlYReVlEDuuvl0XEWkP5PsAs4FT9qcYhInuAAUBkHXb3rmT3zsrXXEMb7UXkOz3vHhH5W7UsASLyud7OehE5pVLZ+0QkTT+3U0TGi8gk4EHgUv16N+l5Z4pIsojY9OtIq2bHdP1/Jl//XGt05Ski1+r15IjILyLSWT8uIvIfEckQkTz9d+sxPxQi0lFEFuq2bBOR//NUW55GRAJE5E8R2aRfy2MebVAp1eJeQFil97cDs3xtUyOv42zAT3//LPCsh9o5AEyo4XgKcLP+fjbwpP7+aTSxs+iv0zgxH3NSXcC1wFCgP5AK7Kx0bjaQDYxAc2HxCTBHPxcKHAHuAgL09Ej93B3AKqADYAXeBj6r5frORBsPr3zscb18HBALrACeqKX8NcAyYC/QFfDXbc6txe5g4BAwUz83BG3TSb9a6l8MvKlf4yAgExivn3sUsAEX65/13cB+/X0vvZ32et4koFulch9Xa+c8oBtwOnA94ASG6OdGAHnARLQOXSLQWz+3CLhef38+sAfoo1/bQ8AK/dw5wDogAhA9T4IH/z8SKtkfCuwC+vr6/7aR1yJAiP7eAqwGRnmqvRbZY1dK5VdKBgMtcgZYKfWrUsquJ4+LmDc5DJzUC0YTmgS0XW42pdRSpf8ia0Ip9b5Sap1Saiua8PcUkfBKWb5WSv2pX+snaOIGMAVIV0q9qJQqVUoVKKVW6+duBP6plEpVSpWhCdnFItJQ/0ZXAo8rpTKUUpnAY8CMOvKHAnuUUvuUUuVo4rq7DrsPKKX+p5SyK6XWA3PRxLkKItIRGAvcp1/jRuC/1WxZp5T6SillA15CuwGMAhxoN7W+ImJRSh1QSu2t7QKUUvOVUnuVUkuABUAh2k0Z4DrgfaXUb0opp1IqTSm1o4ZqbgSeVkol69f9FDBI77Xb9M+pN9qNPlkpdaSOz7RJKKWO6J8tSqkCIBnthtTiUBqFevJ4h8ljutUihR1ARP4tIofQ/oEf8bU9buBatCEFb5KI1jOtzvNovbZfRWSfiNxfWwUiYhaRZ/RH+3w0QQJt6/Rx0iu9LwZC9Pcd0XrJNdEZmKcPB+Wi/VM7gIbGFmtPVbcVB/VjteGP1js+TlG185Xt7gyMPG6bbt+VQE3hhNoD2bowVbalskBVtKuUcqI99bRXSu1Be3J5FMgQkTlSx7CjiEwWkVUikg1s1u09/j3U9VlXpjPwSqXrykbrbSYqpf4AXgfeAI6KyDsiEtaAOpuMiCQBg9F6ui0S/X9lI5AB/FapE+N2mq2wi8gCEdlaw2s6gFLqn0qpjmg9qdt8a23t1Hcdep5/Asd7hd6yaziauJy0GkTvOd+llOoKTAXu1Md2F6AJ1bvHrwOth34lMAEIR3vyAGhIZONDaEMHtZ2brJSKqPQKUEql1ZC3pp7PYTSROk4n/VhNuNpzOgQsrmZbiFLq5lrsiBKR0Gq2VL6OjsffiLZSqcNxW5VSnyqlxurXotCG7E6yWZ8/mAu8gHbzG4jWYz/+PdT1WVe/thurXVugUmqFbs+rSqmhQD+gJ3BPA+psEiISgnZtd1R7Wm9RKKUcSqlBaN/vCE/OTzRbYVdKTVBK9a/h9W21rJ8CF/nCxoZQ33WIyNVoj/ZX1jXc4S5EJExEpgBz0MZot9SQZ4qIdBcRAfLResoOpdQEYD3ao3p/pVR/tPH4LP0VBLgSreMHoJ2I3CHaZGeoiIzUz80C/l1p4i628s2wGkeB6GrDP58BD+nlYtCe6mpbHngUbdy4U6VjwUBBzdn5AW24aYaIWPTXcNEmYquglDqENr7/tD6BNhBtWKTyTXyoiFyoDzPdAZQBq0Skl4icpYt2KVCC9l0ctzlJTixZ9UcbtslE6yScwYknDNBWQM3Ub9AmEUkUkd41XNss4AER6QcgIuEi8hf9/XARGSkiFrQnmtJK9ngEva25wCdKqa892Za3UErlos1reCwOZbMV9roQkR6VktOAmsYKmz36qoT7gGlKqWIPN/e9iBSg9cj+iTaWO7OWvD04MUa7EnhTKbVIP/c0mmDmisjdwIdoQwtpwHa0G0GD0IcnJqI9FaQDu4Fx+ulXgO/QhoMK0J4ERtZSzw40Id+n29UeeBJYizYksQXthvRkLab8AWwAxou2csUf7Qa1qw67zwYuQ+tZp6P1pE9adaNzOdrE52FgHvAvpdRvlc5/i7ZKKQdt7P1CfbzdirY09ZjeRhzaahiAL/W/WSKyXrfpduALvZ7pVPoulFJ/on3f/0GbRF1M1Sea4/nm6dcyRx9a2wpM1k+HAe/q9R9Eu5m/UMs1Nxm9Y/EekKyUeslT7XgDvYMRob8PRHvC9ZhutcidpyIyF23FgBPtB3ZTLY/ozRrRltVZ0f5BAFYppW7yoUmNRkQuAF5DW4GSC2xUSp3jW6tcQ0TOBV4GzGgTjf/2sUmNQkQ+Q1spFIPWs/+XUqrF7VkQkbHAUrQbs1M//KBS6kffWdU49Ce1D9B+WybgC6XU4x5rryUKu4GBgYFB7bTIoRgDAwMDg9oxhN3AwMCglWEIu4GBgUEro6G7+NxKTEyMSkpK8kXTBgYGBi2WdevWHVMNiHnqE2FPSkpi7dq1vmjawMDAoMUiIgfrz+UjYTcwaC2kZxewYvN+YiNCOG1QV1+bY2AAGMJuYNBojhzL5+onPiGnoASAq88dzm0Xn1ZPKQMDz2NMnhq0KpZu3Mfkf7zNWbe9wae/rvNoW3+s210h6gBfL9rs0fYMDBqKIextgG3707ng/vcYe9OrPPfx77TWTWklZTYenPUDx/KKKCgu4z9zFrMzJcNj7UWGBVZNhwbWktPAwLsYwt4C+HjJem6YNZdn5i2kpNzmcvmH3/6R1Iw8ysrtfPnHJhasqdEFisfIKyzh7te+Zcrd7/Loez9TZrPXX6gRFBaXUVpete5judW977qPSSP7MHVsP8xmE/FRoTx2/eT6CxkYeIEWO8auVDmUrwEJQfxPqb9AC+XbNdt49pvFAKzclUJxmY3HLzvbpTqy8orqTHualz9fwuINmivw+cu3kxAdxo3nj3Z7OzERwYzq35lVW7WFA53iIzilh+eiJppMwiPXnsPDM89G81dlYNA8aJHCrlQZKnsG2DZq6aDrMIXd53odeY9A+Wqw9EHCn0ZMHg856jJbU45WSW9JSa8lZ+1MPa0/ny/YAEB4SABnDO7uFtsaSlpmbrV0nkfaERFeuv18flqZTGm5jXNG9iYksDaHi+5t18CgOdEihZ2yJRWiDkDxe6iQ2xBTcIOrUIVvQ+k8vb7DqPynkYhn6y7kA4Z0TWTO8k0V6aFdXY8MdvcV4xjSM5HM3CJOH9SNhBivBL2p4KyhPdiwS3O+KQLjhnjuxmLxMzPtNI/FLzAwaBG0TGGXgGoHLCAW1+pwHKo73UyYPLgXpeU2libvp0t8NDdOHNGoes4a1tPNljWcyyYOITYyhF0pmQzp1YGR/U5yA27gY35dvYNfVu8kITqUmy4c45UnnebKsk372LrvCAO7t2f0AFfixjQfWqaw+4+FgGlQ+h3gh4Q9ihYboeFIwNmo0u84HmFMAlxzHa5UOThSwdQOMQW5VNZVLhjZnwtGtuxe6PhhPRnvw5uLQe2s3naQf759wsV5enYBL/y9tmBVrZvvl23j8fd/qUg/ecO5nDOqpkBTzZsWKewigkS8gHLcBxKAmELrL1S9joCJEPk+qnw1YumLBDQ8SpVyHNHG+B0pIJEQ9T5i6eeyDQYGzYEte49USW/eU1to2NbPgjU7T0obwu5lxFyvL5y6y1vHINYx9eZTqhSV/yiUrwPLQG0oyJGin8xBFbyIRL3fJFsMmicOp5MP5q9hR8pRhvTswGUTh/jaJLczoFtClXT/aum2RPuY8Krp2PBacjZvWrSwewtV+BqU6HF0HQfBXH3yzzPrsg18zzvfrOT9H1YDsHDdHsxmE385a5BH2jpwJJvH3vuZzNwipozpy00X1N/pcAcj+3Xm8b9N5tc/d9IuKpRbLxrrlXabI7dePJas/CK27j3CKd0TPbIs1xsYwt4Q7NUcqpnCQUWBMxskEAm+xTd2NYB1Ow8xd+FmQoOs3DD9VKLDG75yyAA27q4aSnfjrjSPCfuDb/3A7tRjALz3/Wp6dYpj3NAe9ZRqGnvTjvHQ2z+SnlXAOSN7cc+VZ2Eytb7lm06n4pslW0jNyGXsKV0Z0qtDjflCAq08d+s0L1vnfpos7CLSES1SfTu0gLPvKKVeaWq9zQkJGI8q+/VEOnAqBEwG2y7w64yY2/nQuto5cCSb21/8mnK7A4Bt+9L5+NGrfGxVy6J3pzjW70ytSPfr6rnv+vCx/CppT633r8y/3v2ZPfrNZO6izfTvlsCUMa1vvui1r5bw8c+a76BPf13Hm/f8pVZxbw24o8duB+5SSq0XkVBgnYj8ppTa7oa6q6DKN4FtE1gGIv6e6TXVhNM6nSNOfyL9txAUPBAJPFc7YR3pNRsaw/YD6RWiDrAzJYOSMhuBVheXhrZBtu49wt2vf0dWXhEd4yNIiApjWN9OXO7BMfbxw3rw3bJtAARaLYwZ6Pmldpm5hVXTOYW15GzZLF6/t+K9w6lYtmlfncJudzgRAbOpZXpdabKwK6WOAEf09wUikgwkAm4VdlX6Oyr3VrSHAhNEvOLyEsXGkF9cyg2vzyU5NYPQwDBe+c7CHvYAACAASURBVNtAhnp342aj6dkxDj+zCbvDCUCXhKhmIeqFJWUcOZZPYmw4QQF1L1PNLShhx8GjdIyPJNGLE1lPzv61wvXCoaO5XD5xiMeGYI7z4DUT6dc1gczcQsYP60GX9tEebQ/gvNF9+ehnLehNkNXCmR7cPOZLOsZHcCgjt0q6Nj75ZR2vf7UUEeGOS8/gkvHe60S6C7eOsYtIEjAYWF3DuRuAGwA6derkct2qZC6aqAM4USXzvCLsny3ZSHKq5iGwoKSM575ezOf3Xunxdt1B9w4xvPD3aXzx+0ZCgwO4rRlMiu0+lMmtL3xFTkEJsRHBvHXvJXRuF1lj3oPpOfzt6TnkFJTg72fmmVumei2YRUFxWdV0UVktOd2H2WTiwjMHnnT897W7mL9iO7ERIdxy4RjCQ9znRfL2S06nb5d2pGflM/aUriQlRLmt7ubEwzPP5snZv3HoaA5nDO7O+acPqDFfakYur3yxmOMOUF/8dCFnDO5GfNSJJdWZOYWs3n6Q+KhQhvdxXcu8gduEXURCgLnAHUqp/OrnlVLvAO8ADBs2zHW/saZqSxtNnu/NANgqDWUAlNtb1gqYMQO7Mmag62KYlVfE4g17CQ2yMn5YT7dNqL3zzYoKH+aZuUW8//0qHvtbzV4Rv/h9Q0XecruD/81f7TVhv2ziEF79YgmgueNt6Fpmu92BUyn8Le7519qwK5UH3vqhQmjSMvN4/a6L3FL3cSYMb/0bx2IiQnj5jgvqzZdfVEplr9ZOpSgsKSMeTdiPHMvnmic/JTu/GICbLhjNdVNHecTmpuCWX5+IWNBE/ROl1NfuqPOkNkLvQDn2Qfl6sAxCQu/yRDMncdHoAXz353aO5hbiZzJx46Tm9yW6m5z8Yq558lPSswoAmDLmAP+6zj1PR05VPV37Pd7fYq6StviZa8npfmZMGkb/ru04kpXP8N6diI0MqbfMl39s5KXPFuF0Km44/1S3/MNv3ZdeRWi2tOHNQ96gZ6c4BvdMrPBtNLJf5ypPMQvW7KwQdYDPF2xoncIummu794BkpdRLTTeplnZMkUjUR56qvlYSosL46v4ZbEs5SvuoMDrH1Txs0JpYufVAhagDzF+xjQf+Ot4tvdDrpo5kw65UCorLiAwN5Opza/d9M2PScFZsOcC+tCwiQwP5v0tOb3L7NeF0KkRO9tI4uGcHBjewjmN5RbzwycKKG9WseSv+v73zjo+yyvr490xJT0gjoSYQSmhCwIAFpEhVRBB1ZWWx8O4irmV11Xctq+6uq2t519eyq7voWvC1Y0PBAjZApYMI0ktCKAmQXieZue8fE4YkpDC95H4/Hz4f7jPPc+95MjO/uc+5557DuGF9yOjq3pPloIxOiOAQ98G9vZeGWAMmo4F/3HE5327ai8EgjMnq1WgBNS66cZ6qpu1AwRMz9pHAHOAnETmZcvFepdTSVq4JKuKiIjivX3Alrtq4M4/vf9pPj86JToevJcY1zn0TFxXhsdnygJ6deP9vc8nNLyK9U0Kr/uLEuChe/9McCorKSIyLIiKs9YXf0opqPlzxE0opZow+64x80S8uXs1Ln6whzGTk3msnMOkc17aPV1ZbTnv6KKty3y8/tG83Hr3xEpZ8v52OCXYfu8a7hJlNTByR2exrU0cOZPW2HL5cv4vEuGgemOv9dT5XEH+UScvOzlbr16/3+bi+wmZTPPX2t3y7aS9pqfE8MHfyGT3Ke4r123O5+e/vYa33e8y95BxudFIQnn13BW8v30RsVAR//s1FjBgQmItEJ6mprePav7zO3kMnAMjomsTCB2YT3spTxo6cfOb8+XVHO8xk5IunbyQ60rmEcgBKKe58djErNtvD6ob27cpzd16ByYfuI43vqK2z+tQ1eBIR2aCUym7rPL3z1At8uOIn3ly2EYDDx0t46OUveOb3M302/reb9jpEHeDrjXucFvZbrhzNLVd6x/XhDQ4cLnSIOsC+QyfYf/gE/dJTW7ymuEEharAv0FZWW1wSdhHh8ZunserH/VitVkYNydCiHsL4Q9SdQQu7F8graFoxqLiFM71Dt5TGMbrdUwKvMpSnSYqPJsxkdGzIMpuMJLeRPiGrT1f6dEt2bOMfO7SXW09WRoOBMUN7uXy9RuMptLB7gdFZvXjjiw2OWfPYYd7N99GUKy4cwsH8Ilb+uJ8enRO459oJPh3fHyR3iOavN1zMM++uQCm49coLSI5vXaQjws28cM8svtm4h/AwE+OCZeeZRtMG2sfuJTbvOsSqLftIS01g2qiBp0Vc+MtHpwk+1v6cy/++9Q21dVZumHF+iwt7rpBTvolNRUsIN0YzquM1xJp9sz9E4xrax+5nsvp2Javv6fVJrTYbDyz4lGXrdpIQG8XjN01jSB/n65iCPZ+FyRicuSw0Z0Z5ZQ13PfsRlTW1ADzwwqf075F6mrvNFU7U5LLo4ANYlX3TXX7VHq7r9U+3+9X4H60KPmbJdz/zxdqdKAWFpZWNynCdKQfyC5n+l1cYftvT3PiP96mq/9JrQo/CskqHqIP9xzy/sKyVK86co1V7HKIOUFCzD4ut2iN9a/yLFnYfU1rR+ItTUuH8F+lv73xNTkERSsEPO3JY+NUGT5mnCTC6JMWRmZZyqp0cR2Z6SitXnDmpkb0wyqmH9uTwHoQZAnPDjcY5tCvGx0wYkclrn613bEu+avyZ7m08RXFFVattjeew1NbxypK1HCwoZuyw3j4vyG0yGXn+v69g0Vc/UltnZebYwcREhnuk7+TwdGZ2/xMbCz8mwhjD6JTrPNJvSyhbEWByqUaxxjn04qkTKKV4Zek6Nuw4SGZaR26Ycb5L2+yPF5ezZlsuKYkxLmWHe//7n3jozeUoBZFhZl66/Ur6d285XlvjOg+99LkjRzrAU7fNcCmpWnvHVvoIVL4CGJDYPyDR1/vbpKBEL556gbeWb+K591YBsGZbDlab4rarxjjdT3J8DFNHDnDZjpnnn0VGpyT2Hy3k7D7dSOtoX0j7PHcnr+3aSFxYBPcOG0e3mNCPX/c2GxpUTwLYsCNPC7uTqNqt9aIOYEOVPQaR0xFDaKYIDgS0sDvBz/uPttr2JVkZXcjKOJUQ6ufCfG5a+QHW+iew3SXHWDZtnldtOFLxFUU1W0mMyKJTVPDsUnWGfukpjUrU9fOQf7tdoZq6Cm2gvJ/fvj2jhd0Jsvp05bPVOxq1PUV1TS0KXK5wtL24wCHqAHtKTlBjrSPc6J23+EDpe/x4/CFHe1jHh+keO9UrY/mT+66bSExkOAcLihkztLfLScLaNeahEHYeWH6wtyMuQ4yd/WtTiKOF3QkuHzcEq03V+9hTuObi4R7p97VP1/GPRatQKObPGMncac7XUs1K6kK40USN1R6+lpXcxWuiDvbZetN2KAp7bFQEf7x+kr/NCGpETJDwH7CsAQlDwjzzvdG0jBZ2J/nF+CyP1kA8cryUZxetdOTbfv6D75gwoi9pqc7lfe/VIYnXxs/izd2b6RAWwa2DvZveNcrcFRo8YUebg6/i+5frd/HNxj10T4nnuqkjPFb1SHM6IiYI1ymHfYX+JPuZimoLTQOTKqstLvU1PKU7w1O6e8CqthmQeCsWazFFNT+RFJFFZsJ8n4zrKb7bsp+7n/vE0c4vLOP+AM2trXGOmto68gqKSU2IJSbKM6GhrvLF2p3syi1geP80zhnou5oOWtj9TEaXJC7IymDl5n0AnDsonT7dO7Zxlf8xG2IYnvq4v81wmY07DzZqN41+aW9U1R0hr2wRBjGTFjcbc5DGmh8rKmf+4++Qm19MbFQ4/3vbZQzxU9Wp1z/fwFNvfwvAwk/X8cTNlzJmqG8SzWlhr0fV7kZVLUIMcRB1HWJoPeWrpzAYhCduvpTVW3NQSnHuoB6NSnFpvEPD3ZzNtdsTtdZS1h6eTbXVHuVVUPEV53Z5E5HgS1L3+hcbyM23p8kuq6zhn4tWsuDuq/xiy1cbdjv+rxR8s3GvFnZfoqyHUYWzQJWhACyrfVpf1WgwMHJwT5+Np4FJ5/TjeEkl32zcTfeUeG6b5fx+hFChxLLVIeoApZatVNcdJdLsuagvX2G12hq165q0fUn3lHi2NCg+3rVjB5+NrYUdwLIeVIPESpY1KFslYohq+RqN36m21PLwK8vYuDOP/j1SeWDuZKeKC189aRhXTxrmRQuDg0hTFwQTCntElckQi9kYOEXblbJSU/YU+RWLsEkkSXG/Jy7qkmbPnTVxGF+u38Wx4grCw0zcMON8H1t7ittnjaG8qoZduccY3r8711zU5oZRj6FTCgCqdgvqxJVA/d/C0BlDyrd+tUnTNs+9t4qXl6x1tKdfMEiHJrrIkfKl7C1+DoOY6Zd4L4mRgROSaC2+l5yyhVSf/H5iICP1CyLC+jd7fnllDXsPHadLcgef1hr2BTqlgBOIeTDEPYKqfAUkFol70OW+VOUiVPn/AEYk9j4k8mKP2alpTMMdoc21NWdO55iL6RwTeJ9VVf4cVC+ihzmKQmstBTYLYKOmbleLwh4TFe5yjYNQwSOrdCIyRUR2isgeEbnbE336Gom6HEPyxxiS3kDMrlWoUXW5qNL7wVYItmOokrtQ1hNtX6hxifHZfWhYmOrCbN+WINR4F2U9iip/ipNvcaLRTDiCQcKJDNMutNZwe8Yu9qXzfwITgTxgnYgsVkr97G7fQYftBGBtcKAWVBGgy415gwuz+/LM7TPZtOsQ/dJTGHe2FvaQQp2+nyPO3JeYxKcIM/lmv0aw4glXzAhgj1JqH4CIvAVMB9qfsJv7gykT6nbWt4eCsUezp36Tt48/fr+Mamsdv8s6nzn9nc/LroFzB/Xg3EE9/G2GxguIKQ0VMQOqP7QfCB9HcvzziLjuaNiTd5zPVm8nPiaSX4zPCtndxp64q65Aw90eecBpyU5EZB4wDyAtzfkc5MGASAQkvgFVH4KYIHKGfSt1EyprLdz41UdU1tlLnt3/wzJGdOpGZkLgb0zS+JbC0kr+tnA5B/OLGDO0N/MvO/+0wuihjCH+cZTlKsAK5my3RD2voJhfP/IWFfU7u3/cfYgnbpnuIUsDC08Ie3OfstNCbZRSC4AFYI+K8cC4AYkYYiF6TqvnlFhqHKIO9j/W0YpyLeya03j4lWWs2LwXgL2HTtAlOY7po8/ys1W+RcLO9kg/a3/OdYg6wIrN+7DZFAZD6P1QemLxNA9o6PDqBhxu4VwN0CkqhnM7nfqTpcXGMzTFP9ueg4mco0X88oGFjLrhae791xJq66xtXxTkHDjSePH9wNFCP1kS/KSlNi4807Vjh5AUdfDMjH0d0EdEegKHgFnA1R7oN2QREV6ddAXv7N5KdV0tM3sPJC7Mv8mKXEHV7QPLJjD3sYeMeplHXl3GnrzjACxbu5OBPTsxe7JnZnOBysjBGeQu2wiAQYTzBukdyq6S3T+N22eN4b2vt5AQG8k9105o85r3v9nCt5vsGUB/e/kooiLCHK9VW2p5dOGXbNp1iAE9Urnv+okeq0frLm4Lu1KqTkRuBj4HjMBLSqltbVzW7okwmbkmiBdMlWUzqvAaoBowQIcnkMhpXh3zRElFq+1Q5LarxtC1YwcO5hdxQVYvRgwIzfUpX3H1pLO5etKZTQa+Wr+Lvy1c7mgXl1fx1xtO1Rz4z8drWPK9PUbk8PESOsREcPc1bf9Y+AKPLAkrpZYCSz3RlyY4UFXvYhd1ABuq8nWvC/v0CwbxzLsrAYgIMzG5HVQzMhiEqyYE7wQgmNm672ir7byC4kbtg/Xtmto6duUWkBgX7dP8MA0JzVgfjfcxxLfe9gJzLhpORtdkDuYXcc7AdHp20fsDNN5jSJ8uvPZZg3aT9L9jh/Zm+bpdp9rDelNeWcO8R99md95xjAbhj9dP4pKRA31lsgMt7BqXkOj5KMtmqF0Hxgwk9l6fjGvPgqn9zMGIUgpUKWLwzyzWWcYM7c2f/muKPQNoagLzpp/X6PXJ5/YjOjLM7mPvmcr47L68tWwju+vXgaw2xbPvrtTCrgkexBCLJL2OUhZEwtq+QNOuUXUHUEVzwZqHMvVHEl9GDIn+NqtNpo4cwNSRA1p8fdSQDEYNyXC0DU1qKRj8tOdAV3TQuIUWdc2ZoMoeA2t9laq67fbkXiHItFEDGZTRCQCzycgdV4/zix16xq5xG6UU+8q/p9paSs+Y84gyed/frgkyVJMIJlXuHzu8TGS4mRfuvoqc/CISYqNIjPNPTQct7Bq3+Sr/abYW2wtDx5pSmNXjOS3umkZI1HUoywagFiQaifylv03yGiaTkV5dk/1qg3bFaNzCpqxsKz4V6VpWV8CBirWtXKFpj0jEhUjyR0j800jSx0jYEI/2v/SF5czucSNzB9zGpq9+8mjfwUhIzdhV9Zf2zIph53gsv4SmdQxiJNwYQ7W11HEs0hjnR4s0gYqYeoPJ88Wc92/N5an5CzhZDe5PM5/gnSMvEB4gu0D9QcjM2FXFQlTxjajyp1CFs1E1q/xtUrthSpf7iDTGY8DIkIQZ9Iw5198madoRBTnHaFjis7K0ivLiSj9a5H9CZsauqj9p0LKhqpci4aP8Z4+qtqfxbQekR5/NvD6LUMrmVlrVQEcpxc8H8jEZDGSmp/jbnIBFqWqo+QYIh/CxXk8zPOD8TFLSkinItcePnzW6P4md2vcaT8gIO8bOULvZ0RSjf7IlKlWNKroJLCtRhk5Iwr8Qc8txsKGEr0S9qLSSXQePkd4pgU5JvnH7KKW4919LHDsNZ44dzD0BkhckkFCqBlU4G2rr/dwR05D4v3t1zNiEGJ754RE+f/lrIqLDufg3E9pVzvrmCBlhl9j7UbYSqNsBYedD9G9aPFcpBVVvoGq3IuZsJOpylFKe+TBUvAYWez4TbEdRpQ8iSe+6368GgP2HTzDv0XcoLq8i3Gzk77fO4JyB6V4fd/uB/Ebbx9//ZgtzpmTTLSU4Z4ZHqo7yRu4iqqxVTOk0kezELM90bFl3StQBqj9GWe9GjN6tNZDUOYGr753p1TGc4dtNe1m8ciuJcVHcdPko4mMjfTp+6Ai7MRlJfOXMTq5YgCq3zyJU1XuosodAWSFmHhJzi1t2KFXa+ICttPkTHefXUGvZiBjiMJt9v/U42Hhr+SaKy6sAqKm18tIna3wi7Ebj6U8jpmaONaTOaqOy2kJcdGC55JRSPLbzaY7V2F0Xu8r38sig++kW1dX9ziW2yQETiuar8YQq2/Yf5Q//XIzVZvf75xwtZMHdV/nUhtB1iLaCsnzf5EAlUIMqfxZl2ehW3xJ5GcipXBgSfU3LdqgqCo9dTuHxyzlRMJHy0ifdGttTlFfV8Miry5j/+Du88cUGf5vTCHMTMQ0zGX0ybmZaCjPHnso5f93UEa26gdZvz2XS755n/C3PcftTH3i0KEjO0SI27cqj2lLb9snNUGmtcog6gFVZOVjlmdo4EjYEoudhl3ITYIBjI7EV34VSNo+MEehs23fUIeoAW/Ye8bkNITNjdwpTJlh+aP4124nmj58hYsqA5MVgWQvGdCSs5Ufc6qpl1Nae+iEpL3uS6Nib/b5N/9GFy/l8jb0g94YdeSTERnHRef39atNJrr14BD9szSE3v4iE2EhuvuICn419zzUTmDMlG5PR0KZv/+FXl1FWWQPAqi37+eS7bVw2xv1iJIu+/pHH/+9LlILe3ZJ54Z6rnC7uEG2KokdUGgcqcwGIMITTK7qH27adxBB7JyrmFlTBKFAl9oPVH0HEBIiY7LFxApWBPVMxGsQh7mdldPa5De1S2CX2DhR1dl+grRSs++0vGHtAmPuhemLsDJFtF8kVMTc5YsZeq8R5Kmst/H7VElbn5zIoMZWnR19KUoRr25m3Hyho1N6Rkx8wwt4xIYa3HrqG/MIykjtEExHe9G/oXc7Up15RZWm17SrPvbeKk5F9e/KO8/nqHVw+zvnNPnf3u40PDy+lylrF+JQxpER42gdurn8SboAtNNMINGVgRmceu2kai1duIzEuit/OHOlzG9qnsEs4EvcAAEpZofpTey6LiMn2YtQepqWF2fCISYRHXERN9aeAibj4RxBxTdif3fIDn+XaF/dWHcnh4fVf8+SoqW1c1Txn9+tGbn6Roz0ss3srZ/ses8kY8IuW11w0nKffWQFASkIMU871TFGQpn79tvz8LRFrjmVOuvf8viIGVPR1UPGC/YAxHSLGe228QGPM0N6MGer5zVhnSrsU9oaIGCHyktOOq9qtdneKKRMJd+0XV9kqUcW3gmUVytQbiX8OMZ0qbSZiJD7xRWzWQ4hEYTC6nsb0cEVpq21nuGv2hSR3iCY3v4jRWb0YM7SXy321V341JZusPl0pKCpjWGZ3j0VF3DX7Qh588TNq66wMy+zGlAB5kmoOQ+xdqLDRoAohbCRi0DuSfUW7F/bmUDWr7bmjqbMfiPsLEjXL+X4qXgSLfdZG3S5U2UNIwguNzhERjKZurfZz8Hgx2/Ly6du5IxmpzYv/pT378/GB7djqn9Nn9HQ9dt5sMjJvxvkuX6+B4yUVdE6OY1Avz/pXJ47IZPiANEorqunasQNGQ2DHP0j4Of42oV2ihb0ZVPViHKIOqKoPXBJ2bIWN21bnF2Y37z/MvAXvUWWpw2w08szcSxnVr8dp543v3pu3J1/N2vyDDExKZWzXjNM70/iEfyxayatL1wFw/dQR/PZyz+6Ajo+JJD7Gt3HRmuAisH/u/YUhtUnbte3jEnkZcCpiQaKc92m+8d1mqiz2H5laq5X/W7mpxXOHp3bjpsHnaVH3IzlHixyiDvDykrWN1is0Gl/g1oxdRJ4ApgEWYC9wvVKquPWrAh+J/g2qboc9JNKUicTd51o/YUMg+QOwrAdTLyRsuNN9xEQ0Dn2MCdcViwKZ2rq6045ZaluPYd+0K4+3l28iKiKMpA5RHDhSxODeXfjV5LPb/dZ4jWu464pZBtyjlKoTkceAe4A/uG+WfxFDFJLwvGf6cjNV6W8nncfmA0fYfeQ4acnx3DbVf4nNNG3Tq2syE4b3daQfmDgik15dk1o8P6+gmFuefJ8aS+MfhG827gFgzpRs7xmrCVncEnal1BcNmquBK9wzR9OU5Lho3r9zDiWV1cRFhusZnJNs2XOY/3y8GoPBwPz6BeH7FyzlWHEF00YN5Pe/HOvR8USER+ZP5arxQwEY0qdLq+/ZztyC00T9JD/uPhRQwl5dW8fCjZsoqa7hsoH96Z3c8g+Wxr94cvF0LvB2Sy+KyDxgHkBaWlpLp2laoENUYOUbCQYKSyu59X/fd2wO+mnPYYxGA4Wl9o0zby7byODeXZgwvK9HxxURsvqeWd6VPt07EmYyYmkm5cDAnp08ape73PThx3y77wAAb27awsfX/4quHXQIYyDSprCLyHKguU/YfUqpj+rPuQ97GMnrLfWjlFoALADIzs5WLZ2nCT0Kisr49wffU1Ft4RfjhzIss/XwTk9xML+o0Y7PkorqZmzz727ItNQEnvzdDN5avonIcDOJcZHkHClicJ8uXHOx82sy3qKmrs4h6gClNTWsyT3IzLN04rpApE1hV0q1mnRaRK4FLgHGq4ZlTDQa7Ltub/77++w/bA/1XLVlP2/+eQ7dUxO8PnbPLkkkxEZSVFbV7Otx0eGMzvJ/BNE5A9N9kqHSHcJNJrrExXK4tMxxLD0hsHf/tmfcCncUkSnYF0svVappYgiNBsoqaxyiDlBjqWNHTkErV3iOuOgI/vXfv2DqyAGc1WSjUHxMJK/ePzvgUxMEEv++fDpDu3QmIzGBByeO4+xuHkjz60cqqy34ai5qsyk+W72dt5ZtJL+wrO0L3ETcuTER2YM9UPvkN3e1Ump+W9dlZ2er9evXuzyuJnhQSnHFva84YrnDTEZe//McenR2PX2CK9TU1nHHMx+xZlsOMZHhPH7zNIb312s97ZFqSy13PruYNdtySOoQzZO3TmeAl9cz/vLS53y8ahsASR2iee2B2XRMiHG6HxHZoJRqc0XdLWF3FS3s7Ysjx0v55/urqKiyMGvCUL+6HYpKK4mJCsfsozzurlJaUU1eQTFpqQnERDmXlhfsP2SPLvySzbvz6N+jE/deO6HN9L7L1+1i4dJ1RISbuOOX40K2ruv/fbbekaANoH96KgsfnO218aw2GyPnPd0oR/sfr5vI9NFnOd3XmQq7TikQYJTV1vDq7rVU1Fm4qudQesQ6P7O1Kis7SndjNpjpG+v/BF6dk+P467yL/W0GAAlxrqUy9iU7cwq46X8WUVJRTWJcFM/fdSUZrcTCN8eLi1fzyXf2GWJeQQlx0RHcPafl7Ir7D5/gjwuWYrXai2H87qkP+OR/fuNy9shA5mSe/JOUV9W0cKZnMBoMJMRGcbykwnEsOd752bozhN67FuT8euWbPL3tW17c+QO/+OoVCqqc88dZlZXHdjzDX7c/yYPbHuPFfa95yVKNt3hx8Q+OCJ7C0kpeXrLG6T7y8htvAD/YRlqDnKNFDlEHOFFSQWkzUUQNqauzsv1APoeOlThtnz+5ZOQAOsScCh+OiQrn9c83YLN5x3tRWW3hT7+eTJfkDkSFm7nmouGMHNzTK2OdRM/YA4jCmko2nshztIsslWw6cYjJ3c48l/fO0j38VLLd0f6yYCUzu11CYpj3o1A0gcOYYb1Yvv5U8e22coMPyuhEh5gISsrtYt4/PZWEVlIN19TWcdMTi/hxz2EMItw5exxXXuihgthepntqAm/8eQ7//uB7Fq/axvYD+Ww/kE+1pZb/muZ+oZ2GvLJkLc+9vwqAG2ac7/H+W0ILewARZ44gOTya4zX2RzYDQnqMc4JsNjSuKCQIJtFvczDx6+nnsWn3IUrKq0nqEM3cS5xPfTvl3P5ERYSxedchBvTs1OYmrOT4GF68ZxaLvv6RqHAzv5qS3eqO2S/X7+LHPfY6qTalePqdFVwxbkjQ7IxOSYiluc2Z5gAABglJREFUsqZxzdh3v/qRqpparp96DtGR7udkOnSshH++t8rR/tcH3zNxRCZpPgj11d/4AMJkMPDCBbP4y6bPqait4Tf9zqNffGrbFzagT2wGE1PHsiz/GwRhdvoVxJk9XxVK4z0y01L44NG5HDpWQreUeKdrmp5kdFYvRmed+RpLj86J3Hn1uDM6V5Am7eAjMy3FkdMH7O6nV5euY1fuMZ75/Uy3+6+qOb3YeFW1awXInUULe4AxKKEz71x4nVt9zO15NTO7TsUoRmLN3l2k0XiH2KgI+qV7Lo2EUopPvvuZg/lFjBzckyF93ItBH5/dhw9X/MTGnXkYDcLts8YGzWz9JHMuyqaqppYv1u4gr+DUOsH6HQc90n+vrkmMzurFis17ARg5uCd9unu6tmzz6HDHdo6ylULdLjCmIcbQDG/T2Itgv7xkLQBGo4Hn77qCoX3dS+1QZ7Wx//AJ4qIjSE0M3qfCjTvzmP/4O44i4QN6pvLq/Z4Jf7TabKz9ORebTXHuoHS3K17pcEdNm6i6HFThbLAVgERB/L+QcN8s7mh8y5frdzv+b7Xa+HbTXreF3WQ0+GwG6k2GZXbjgbmT+WjFVpI6RHH7rLEe69toMHDeoB4e6+9M0cLejlGVL9tF3d5AlT+jhT1E6ZbSoVElJ51KoTGXjBzIJSNDJ6GZjmNv1zT5XdfRMyHLH6+bxDkD0+mcFMdV47OYOWawv03SeBH9TW7HSPSvUTVfg/UgSAck5g5/m6TxEh0TYvjHHZf72wyNj9DC3o4RYydIXmoXdkMnxKAjaDSaUEALeztHJNytmqwajSbw0D52jUajCTH8EscuIseAHBcvTwaOe9AcfxIq9xIq9wH6XgKRULkPcP9e0pVSbcaY+kXY3UFE1p9JgH4wECr3Eir3AfpeApFQuQ/w3b1oV4xGo9GEGFrYNRqNJsQIRmFf4G8DPEio3Euo3AfoewlEQuU+wEf3EnQ+do1Go9G0TjDO2DUajUbTClrYNRqNJsQISmEXkYdEZIuIbBaRL0Ski79tcgUReUJEdtTfywciErQp90TkShHZJiI2EQnK0DQRmSIiO0Vkj4jc7W97XEVEXhKRAhHZ6m9b3EFEuovI1yKyvf6z9Tt/2+QqIhIhImtF5Mf6e/mzV8cLRh+7iMQppUrr/38rMEApNd/PZjmNiEwCvlJK1YnIYwBKqT/42SyXEJH+gA34N3CnUiqoKqmIiBHYBUwE8oB1wC+VUj/71TAXEJHRQDmwUCk1yN/2uIqIdAY6K6U2ikgssAGYEaTviQDRSqlyETEDq4DfKaVWe2O8oJyxnxT1eqKB4Pt1ApRSXyil6uqbqwH3Kh/4EaXUdqXUTn/b4QYjgD1KqX1KKQvwFjDdzza5hFJqBVDobzvcRSl1RCm1sf7/ZcB2wL2afn5C2Smvb5rr/3lNt4JS2AFE5GEROQjMBh7wtz0eYC7wqb+NaMd0BRoWu8wjSEUkFBGRHsBQYI1/LXEdETGKyGagAFimlPLavQSssIvIchHZ2sy/6QBKqfuUUt2B14Gb/Wtty7R1H/Xn3AfUYb+XgOVM7iWIaa4Sc1A+CYYaIhIDvAfc1uRpPahQSlmVUlnYn8xHiIjX3GQBm7ZXKTXhDE99A1gCPOhFc1ymrfsQkWuBS4DxKsAXPJx4T4KRPKB7g3Y34LCfbNHUU++Pfg94XSn1vr/t8QRKqWIR+QaYAnhlgTtgZ+ytISJ9GjQvBXb4yxZ3EJEpwB+AS5VSlf62p52zDugjIj1FJAyYBSz2s03tmvoFx/8A25VST/rbHncQkY4no95EJBKYgBd1K1ijYt4DMrFHYeQA85VSh/xrlfOIyB4gHDhRf2h1MEb3AIjIZcCzQEegGNislJrsX6ucQ0QuBp4CjMBLSqmH/WySS4jIm8BY7Cli84EHlVL/8atRLiAio4CVwE/Yv+sA9yqllvrPKtcQkcHAq9g/WwbgHaXUX7w2XjAKu0aj0WhaJihdMRqNRqNpGS3sGo1GE2JoYddoNJoQQwu7RqPRhBha2DUajSbE0MKu0Wg0IYYWdo1Gowkx/h9eiUbcIEfGbAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f9618ebc6d8>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# %load -r 95-129 tp0/solution.py\n",
    "\n",
    "def sampleSpace(nbSamples=500):\n",
    "     '''\n",
    "     Sample nbSamples configurations and store them in two lists depending\n",
    "     if the configuration is in free space (hfree) or in collision (hcol), along\n",
    "     with the distance to the target and the distance to the obstacles.\n",
    "     '''\n",
    "     hcol = []\n",
    "     hfree = []\n",
    "     for i in range(nbSamples):\n",
    "          q = qrand(False)\n",
    "          if not coll(q):\n",
    "               hfree.append( list(q6_to_q2(q).flat) + [ dist(q), collisionDistance(q) ])\n",
    "          else:\n",
    "               hcol.append(  list(q6_to_q2(q).flat) + [ dist(q), 1e-2 ])\n",
    "     return hcol,hfree\n",
    "\n",
    "def plotConfigurationSpace(hcol,hfree):\n",
    "     '''\n",
    "     Plot 2 \"scatter\" plots: the first one plot the distance to the target for \n",
    "     each configuration, the second plots the distance to the obstacles (axis q1,q2, \n",
    "     distance in the color space).\n",
    "     '''\n",
    "     htotal = hcol + hfree\n",
    "     h=np.array(htotal)\n",
    "     plt.subplot(2,1,1)\n",
    "     plt.scatter(h[:,0],h[:,1],c=h[:,2],s=20,lw=0)\n",
    "     plt.title(\"Distance to the target\")\n",
    "     plt.subplot(2,1,2)\n",
    "     plt.scatter(h[:,0],h[:,1],c=h[:,3],s=20,lw=0)\n",
    "     plt.title(\"Distance to the obstacles\")\n",
    "\n",
    "hcol,hfree = sampleSpace(100)\n",
    "plotConfigurationSpace(hcol,hfree)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimize the distance under non-collision constraint\n",
    "Finally, let's use one of the classic python solvers (from scipy) to search for a robot configuration that optimizes the distance to the target, under the constraint that the distance to collision is positive.\n",
    "For that, we need to define a cost function (taking the robot configuration and return a scalar) and a constraint function (taking again the robot configuration and returning a scalar or a vector of scalar that should be positive). We additionally use the \"callback\" functionnality of the solver, to render the robot configuration corresponding to the current value of the decision variable inside the solver algorithm.\n",
    "We are going to use the \"SLSQP\" solver from scipy, which implements a \"sequential quadratic program\" algorithm and accepts constraints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully    (Exit mode 0)\n",
      "            Current function value: 7.270182709684149e-07\n",
      "            Iterations: 44\n",
      "            Function evaluations: 334\n",
      "            Gradient evaluations: 44\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([-1.34173219,  1.11253745, -0.43407717,  0.        ,  0.        ,\n",
       "        0.        ])"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# %load  -r 133- tp0/solution.py\n",
    "\n",
    "def cost(q):\n",
    "     '''\n",
    "     Cost function: distance to the target\n",
    "     '''\n",
    "     return dist(q2_to_q6(q))\n",
    "     \n",
    "def constraint(q):\n",
    "     '''\n",
    "     Constraint function: distance to the obstacle should be positive.\n",
    "     '''\n",
    "     return collisionDistance(q2_to_q6(q))\n",
    "     \n",
    "def callback(q):\n",
    "     '''\n",
    "     At each optimization step, display the robot configuration in gepetto-viewer.\n",
    "     '''\n",
    "     viz.display(q2_to_q6(q))\n",
    "     time.sleep(.01)\n",
    "\n",
    "def optimize():\n",
    "     '''\n",
    "     Optimize from an initial random configuration to discover a collision-free\n",
    "     configuration as close as possible to the target. \n",
    "     '''\n",
    "     return fmin_slsqp(x0=qrand(check=True),\n",
    "                       func=cost,\n",
    "                       f_ieqcons=constraint,callback=callback)\n",
    "\n",
    "optimize()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the output of the solver. It always returns a variable value, but sometimes the algorithm fails being traped in an unfeasible region. Most of the time, the solver converges to a local minimum where the final distance to the target is nonzero"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}