2_trajopt_geom.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# First trajectory-optimization problems\n",
    "This notebook drives you through simple trajectory optimization problems, considering only \"kinematic\" version i.e only motion variables (no forces or torques yet). \n",
    "We will first look at a problem where the trajectory is described by a sequence of configurations $[q_0 ... q_T]$, with the only constraints is to have them sufficiently closed to each other.\n",
    "The we will look at a real integration, where each timestep is described by a configuration $q_t$, velocity $v_t$ and acceleration $a_t$ (all in the configuration space. \n",
    "This will allow us to formulate a proper contact constraint, even if only partially described from the motion point of view. \n",
    "The next notebook will drive us to complete version of the trajectory optimization problem, including forces, torques and complete contacts.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import magic_donotload"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up\n",
    "We will use the following tools:\n",
    "- the ur10 model (loaded by example-robot-data)\n",
    "- pinocchio.casadi for writing the problem and computing its derivatives\n",
    "- the IpOpt solver wrapped in casadi\n",
    "- the meshcat viewer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_kine_imports\n",
    "import time\n",
    "import unittest\n",
    "import numpy as np\n",
    "import casadi\n",
    "import pinocchio as pin\n",
    "from pinocchio import casadi as cpin\n",
    "import example_robot_data as robex\n",
    "from scipy.optimize import fmin_bfgs\n",
    "from numpy.linalg import norm\n",
    "\n",
    "from utils.meshcat_viewer_wrapper import MeshcatVisualizer\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The notebook is written first for a UR10 robot model configured in fixed basis, then with a Talos pair of legs configured in floating mode. Replace the robot by your own URDF..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "robot = robex.load(\"ur10\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The robot will start from an arbitrary initial configuration (make the solver life simpler by avoiding the singular 0 configuration) and aim at an object located a placement $M_{target} \\in SE(3)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_kine_configurations\n",
    "Mtarget = pin.SE3(pin.utils.rotate(\"y\", 3), np.array([-0.5, 0.1, 0.2]))  # x,y,z\n",
    "q0 = np.array([0, -3.14 / 2, 0, 0, 0, 0])\n",
    "endEffectorFrameName = \"tool0\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open a viewer, display the robot and markers to visualize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_kine_viewer\n",
    "# Open the viewer\n",
    "viz = MeshcatVisualizer(robot)\n",
    "viz.display(robot.q0)\n",
    "# Add a vizualization for the target\n",
    "boxID = \"world/box\"\n",
    "viz.addBox(boxID, [0.05, 0.1, 0.2], [1.0, 0.2, 0.2, 0.5])\n",
    "# Add a vizualisation for the tip of the arm.\n",
    "tipID = \"world/blue\"\n",
    "viz.addBox(tipID, [0.08] * 3, [0.2, 0.2, 1.0, 0.5])\n",
    "\n",
    "\n",
    "def displayScene(q, dt=1e-1):\n",
    "    \"\"\"\n",
    "    Given the robot configuration, display:\n",
    "    - the robot\n",
    "    - a box representing endEffector_ID\n",
    "    - a box representing Mtarget\n",
    "    \"\"\"\n",
    "    pin.framesForwardKinematics(model, data, q)\n",
    "    M = data.oMf[endEffector_ID]\n",
    "    viz.applyConfiguration(boxID, Mtarget)\n",
    "    viz.applyConfiguration(tipID, M)\n",
    "    viz.display(q)\n",
    "    time.sleep(dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viz.viewer.jupyter_cell()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are mostly interested by the robot model. We also need to access some particular frames."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_kine_modeldata\n",
    "robot.q0 = q0\n",
    "\n",
    "# The pinocchio model is what we are really interested by.\n",
    "model = robot.model\n",
    "data = model.createData()\n",
    "endEffector_ID = model.getFrameId(endEffectorFrameName)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we will need the same helpers for getting around the Casadi SX/MX issue."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_kine_casadi\n",
    "# --- Casadi helpers\n",
    "cmodel = cpin.Model(model)\n",
    "cdata = cmodel.createData()\n",
    "\n",
    "cq = casadi.SX.sym(\"x\", model.nq, 1)\n",
    "cpin.framesForwardKinematics(cmodel, cdata, cq)\n",
    "\n",
    "error3_tool = casadi.Function(\n",
    "    \"etool3\", [cq], [cdata.oMf[endEffector_ID].translation - Mtarget.translation]\n",
    ")\n",
    "error6_tool = casadi.Function(\n",
    "    \"etool6\",\n",
    "    [cq],\n",
    "    [cpin.log6(cdata.oMf[endEffector_ID].inverse() * cpin.SE3(Mtarget)).vector],\n",
    ")\n",
    "error_tool = error3_tool"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are done for the set up, let's noew write the problems we are interested by."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimizing a sequence of configurations\n",
    "We will first consider a problem where the trajectory is represented by a sequence of configurations.\n",
    "\n",
    "Decide: $Q = [ q_0 ... q_T ] \\in R^{n_q \\times (T+1)}$\n",
    "\n",
    "Minimizing:   $sum_{t=0}^{T-1}  || q_t - q_{t+1} ||^2 + || e(q_T) ||^2$\n",
    "\n",
    "Subject to: $q_0 = robot.q0$\n",
    "\n",
    "The hyperparameters for this problems are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_kine_hyper\n",
    "T = 10\n",
    "w_run = 0.1\n",
    "w_term = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For that, we first define an optimization problem in Casadi where the decision variables are the T+1 configurations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_kine_casadi_opti\n",
    "opti = casadi.Opti()\n",
    "var_qs = [opti.variable(model.nq) for t in range(T + 1)]\n",
    "totalcost = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now implement the running cost (the sum of configuration differences):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%do_not_load tp2/generated/trajopt_kine_casadi_runcost"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now implement the terminal cost. You can use either the 3d cost (position only) or 6d cost (full placement)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%do_not_load tp2/generated/trajopt_kine_casadi_termcost"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now implement the initial constraint."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%do_not_load tp2/generated/trajopt_kine_casadi_q0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we are good for optimizing the problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_kine_solve\n",
    "opti.minimize(totalcost)\n",
    "opti.solver(\"ipopt\")  # set numerical backend\n",
    "opti.callback(lambda i: displayScene(opti.debug.value(var_qs[-1])))\n",
    "\n",
    "# Caution: in case the solver does not converge, we are picking the candidate values\n",
    "# at the last iteration in opti.debug, and they are NO guarantee of what they mean.\n",
    "try:\n",
    "    sol = opti.solve_limited()\n",
    "    sol_qs = [opti.value(var_q) for var_q in var_qs]\n",
    "except:\n",
    "    print(\"ERROR in convergence, plotting debug info.\")\n",
    "    sol_qs = [opti.debug.value(var_q) for var_q in var_qs]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the progress of the solver is shown by displaying the terminal configuration in the viewer at each new step. We can also visualize the full trajectory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_kine_disptraj\n",
    "def displayTraj(qs, dt=1e-2):\n",
    "    for q in qs[1:]:\n",
    "        displayScene(q, dt=dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "displayTraj(sol_qs)"
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0YAAAGNCAYAAAA1lGFGAAAABHNCSVQICAgIfAhkiAAAIABJREFUeF7snQfcHUW5/wcFBASV3oIQSihCQgmEQEgIqBRpQVCaFBEuF5UQkaJ+KBFBxEJolxvLDZEiIF00UkMIIgmBEFoIVToixT9FQNH3f77jfd672eyec/ac3T1bfvP5nM/7ni2zM9+d3TO/eZ55ZqG+RnJKIiACIiACIiACIiACIiACIlBjAh+qcd1VdREQAREQAREQAREQAREQARHwBCSM1BBEQAREQAREQAREQAREQARqT0DCqPZNQABEQAREQAREQAREQAREQAQkjNQGREAEREAEREAEREAEREAEak9Awqj2TUAAREAEREAEREAEREAEREAEJIzUBkRABERABERABERABERABGpPQMKo9k1AAERABERABERABERABERABCSM1AZEQAREQAREQAREQAREQARqT0DCqPZNQABEQAREQAREQAREQAREQAQkjNQGREAEREAEREAEREAEREAEak9Awqj2TUAAREAEREAEREAEREAEREAEJIzUBkRABERABERABERABERABGpPQMKo9k1AAERABERABERABERABERABCSM1AZEQAREQAREQAREQAREQARqT0DCqPZNQABEQAREQAREQAREQAREQAQkjNQGREAEREAEREAEREAEREAEak9Awqj2TUAAREAEREAEREAEREAEREAEJIzUBkRABERABERABERABERABGpPQMKo9k1AAERABERABERABERABERABCSM1AZEQAREQAREQAREQAREQARqT0DCqPZNQABEQAREQAREQAREQAREQAQkjNQGREAEREAEREAEREAEREAEak9Awqj2TUAAREAEREAEREAEREAEREAEJIzUBkRABERABERABERABERABGpPQMKo9k1AAERABERABERABERABERABCSM1AZEQAREQAREQAREQAREQARqT0DCqPZNQABEQAREQAREQAREQAREQAQkjNQGREAEREAEREAEREAEREAEak9Awqj2TUAAREAEREAEREAEREAEREAEJIzUBkRABERABERABERABERABGpPQMKo9k1AAERABERABERABERABERABCSM1AZEQAREQAREQAREQAREQARqT0DCqPZNQABEQAREQAREQAREQAREQAQkjNQGREAEREAEREAEREAEREAEak9Awqj2TUAAREAEREAEREAEREAEREAEJIzUBkRABERABERABERABERABGpPQMKo9k1AAERABERABERABERABERABCSM1AZEQAREQAREQAREQAREQARqT0DCqPZNQABEQAREQAREQAREQAREQAQkjNQGREAEREAEREAEREAEREAEak9Awqj2TUAAREAEREAEREAEREAEREAEJIzUBkRABERABERABERABERABGpPQMKo9k1AAERABERABERABERABERABCSM1AZEQAREQAREQAREQAREQARqT0DCqPZNQABEQAREQAREQAREQAREQAQkjNQGREAEREAEREAEREAEREAEak9Awqj2TUAAREAEREAEREAEREAEREAEJIzUBkRABERABERABERABERABGpPQMKo9k1AAERABERABERABERABERABCSM1AZEQAREQAREQAREQAREQARqT0DCqPZNQABEQAREQAREQAREQAREQAQkjNQGREAEREAEREAEREAEREAEak9Awqj2TUAAREAEREAEREAEREAEREAEJIzUBkRABERABERABERABERABGpPQMKo9k1AAERABERABERABERABERABCSM1AZEQAREQAREQAREQAREQARqT0DCqPZNQABEQAREQAREQAREQAREQAQkjNQGREAEREAEREAEREAEREAEak9Awqj2TUAAREAEREAEREAEREAEREAEJIzUBkRABERABERABERABERABGpPQMKo9k1AAERABERABERABERABERABCSM1AZEQAREQAREQAREQAREQARqT0DCqPZNQABEQAREQAREQAREQAREQAQkjNQGREAEREAEREAEREAEREAEak9Awqj2TUAAREAEREAEREAEREAEREAEJIzUBkRABERABERABERABERABGpPQMKo9k1AAERABERABERABERABERABCSM1AZEQAREQAREQAREQAREQARqT0DCqPZNQABEQAREQAREQAREQAREQAQkjNQGREAEREAEREAEREAEREAEak9Awqj2TUAAREAEREAEREAEREAEREAEJIzUBkRABERABERABERABERABGpPQMKo9k1AAERABERABERABERABERABCSM1AZEQAREQAREQAREQAREQARqT0DCqPZNQABEQAREQAREQAREQAREQAQkjNQGREAEREAEREAEREAEREAEak9Awqj2TUAAREAEREAEREAEREAEREAEJIzUBkRABERABERABERABERABGpPQMKo9k1AAERABERABERABERABERABCSM1AZEQAREQAREQAREQAREQARqT0DCqPZNQABEQAREQAREQAREQAREQAQkjNQGREAEREAEREAEREAEREAEak9Awqj2TUAAREAEREAEREAEREAEREAEJIzUBkRABERABERABERABERABGpPQMKo9k1AAERABERABERABERABERABCSM1AZEQAREQAREQAREQAREQARqT0DCqPZNQABEQAREQAREQAREQAREQAQkjNQGREAEREAEREAEREAEREAEak9Awqj2TUAAREAEREAEREAEREAEREAEJIzUBkRABERABERABERABERABGpPQMKo9k1AAERABERABERABERABERABBYWAhEQAREQAREQAREQAREQAecuvPBC99e//tWjOProo93tt9/uP/b9E5/4hDBVmMBCfY1U4fqpaiIgAiIgAiIgAiIgAiLQksDBBx/sJkyY4BA//L/QQgu5IUOG+P832WQTN3bsWC+WlKpLQK501b23qpkIiIAIiIAIiIAIiEAbBBBEp5xyihdFJP5OnTrVCyEsSKuvvrrbdtttI3P605/+1G9lijxAG0tDQMKoNLdKBRUBERABERABERABEciCwMYbb+zWWGON/qzvv/9+t8cee/jvbMedjmPCCdE0ZswYx/FK5ScgYVT+e6gaiIAIiIAIiIAIiIAIdEEgbA2aNm1arIUoeBnmJGlWShfgC3aqhFHBboiKIwIiIAIiIAIiIAIi0DsCFmwhLJbCJbr22mu9FUkBGcJkyvtdwqi8904lFwEREAEREAEREAERSIGARaIjK4QRQReCgof5R8HEvCI+rcRTCkVTFjkSkDDKEbYuJQIiIAIiIAIiIAIiUCwCzA9aeumlHRYgEu5xQVHE/uD8IztGEeqKdR/TKI2EURoUlYcIiIAIiIAIiIAIiEBpCRB1jvDciB0i1GFBwnI0efJk/5eQ3ZbYH/xe2kqr4AsQ0DpGCyDRBhEQAREQAREQAREQgToRMNc4mzOEMMJSFDWHCFEUtCAhlIhgh1udBFO5W42EUbnvn0ovAiIgAiIgAiIgAiLQQwJYmljzSPONengTUrr0winlo2xEQAREQAREQAREQAREoDYEcLEjrDcJlzusTrIYlfv2y2JU7vun0ouACIiACIiACIiACIiACKRAQMEXUoCoLERABERABERABERABERABMpNQMKo3PdPpRcBERABERABERABERABEUiBgIRRChCVhQiIgAiIgAiIgAiIgAiIQLkJSBiV+/6p9CIgAiIgAiIgAiIgAiIgAikQkDBKAaKyEAEREAEREAEREAEREAERKDcBCaNy3z+VXgREQAREQAREQAREQAREIAUCWscoBYjKQgREQAREQAREQAREoHoELrzwQr9GEekTn/iE22OPPdxBBx1UvYqqRp6A1jFSQxABERABERABERABERCBAIG//vWvbsyYMY5FXMNp2223dVOnTg1v1vcKEJArXQVuoqogAiIgAiIgAiIgAiKQHoFTTjklUhRxBcQS+5WqR0AWo+rdU9VIBERABERABERABESgCwILLbRQ07Nxq3vjjTeaHqOd5SMgi1H57plKLAIiIAIiIAIiIAIikBGBKPe58KVwtWvnuPB5+l5sAgq+UOz7o9KJgAjUgMC1117r+JHlY5N7J0yY4DbeeGN3//33y2WjBm1AVRQBESgOAd7D7aQ11lijncN0TIkISBiV6GapqCIgAtUjYKLo4IMP9pXDfYNtREJCHI0fP17CqHq3XTUSAREoKIF3333X3Xrrrf5d3NfXF1vK1Vdf3UkYxeIp7Q4Jo9LeOhVcBESgCgRwxUAABRMRjxixtL/BfX/605/cM8884zeNGjWqCghUBxEQARHoOYH333/fXXDBBe4HP/iBe/nll/07GCt+XAq/t+OO0/ZyEVDwhXLdL5VWBESgQgTMfc5GHXGb22STTdzTTz8dORKJFYlzzMWO73z4riQCIiACIpCcwJtvvunOP/98d84553hBtOKKK7qjjz7aDR8+3M2cOdMdd9xxC2Q6adIkZ1b+BXZqQ6kJSBiV+vap8CIgAlUiwAgkH6xCUQkLEm525v/Oj/e0adPc7Nmzow7XNhEQAREQgRgCiKAf/ehH7qc//al766233PLLL++OPfZYd8ghh7gXX3zRbbTRRt6djvexuTzbHFC50MVArcBmCaMK3ERVQQREoBoEWFGdH16sQJYsIAPfsQwFLUT8WLMAYTM/+GqQUS1EQAREIB0CTzzxhDvttNPcJZdc4v7xj3+4gQMHum984xteEP3zn/90b7/9tltllVXSuZhyKR0Bhesu3S1TgUVABKpCANGD+5wl5hsF3eLYh/ixxPfgfvZpnlFVWoPqIQIikCWBO++80+22225u0KBBfoBp8ODB7rLLLnMIpSOOOMK7KRN4QaIoy7tQ/LwVfKH490glFAERqCgBfNQRQ/wgm5UoLHziVldHJHFuUFhVFJOqJQIiIAIdEfj73//uLr/8cnfuuee6e+65x+ex4447+nlDo0eP9t+xGhHQ5iMf+YifX6RUbwISRvW+/6q9CIhADwngOkeaPHmy/4vQYd7QuHHj/Lwh/o9KiCH28bfd9Tai8tE2ERABEagiAeYPnXfeee5nP/uZe+WVV9wiiyziDjzwQHf88ce7DTbYoL/KiCLmEC255JJu5ZVXriIK1SkhAc0xSghMh4uACIhALwkghrAuWahY/ld0pF7eEV1bBESgKATuuOMOL4iuueYa98EHH3ixc/jhh3tXuZVWWmm+YrIfAbXUUku5j3/840WpgsrRYwISRj2+Abq8CIiACLRLAFHEgq9BSxJhY80Nr918dJwIiIAIVIXAe++95y6++GIviObMmeOrNWLECPfVr37Vff7zn/fWomAiWA0udogihNOiiy5aFRSqRwoEJIxSgKgsREAERCAPArjN/b//9//muxTBF3DBUxIBERCBOhF48MEH3cSJE310OeZpLrHEEm7//ff3A0dBd7kgk3/961+OdYueeuopv2Yc4biVRCBIQMJI7UEEREAEREAEREAERKDwBFhvCCH0i1/8ws2aNcuXd+2113ZHHnmkO/TQQ93HPvax2Dogil577TX3wgsvaFHsWEraIWGkNiACIiACIiACIiACIlBYAoTaRgxdccUV7m9/+5uPILfnnnt6MbTddtu1tPwgirAqIYzWWWedwtZTBes9AQmj3t8DlUAEREAEREAEREAERCBA4NVXX/URO3/+85+7Rx991O8ZMmSIF0Nf+tKX2o7IyaKtb7zxhj9/ueWWE2MRaEpAwqgpHu0UAREQAREQAREQARHIg8D777/vrr/+eh9MYcqUKX6NIdzj9t13X/eVr3zFDR06NFExOP/ZZ5/1gkiR5xKhq+3BWseotrdeFRcBERABERABERCB3hIgStz06dPdRRdd5K688krv8kbaZpttvHXoC1/4glt88cUTFxKRNW/ePDdgwACJosT06nuCLEb1vfequQiIgAiIgAiIgAj0hMATTzzhlxogmAKLrJKY/0NkORZjHThwYMflIhz33LlzfR7NAjJ0fAGdWFkCEkaVvbWqmAiIQFUIMILKZ4011qhKlVQPERCBGhJ4/fXXvRDCOnTPPfd4Assvv7z74he/6A444AA3bNiwrqhgfWJdI9YoWnXVVbVGUVc063myhFE977tqLQIiUAIC1157rRs3blz/aCrrGLFGx8knn1yC0quIIiACIuD82mtXXXWVjyh36623ug8++MCvObT77rt769COO+7oPvzhD3eNytYoIqT3aqut1nV+yqCeBCSM6nnfVWsREIGCE8DF5JBDDoks5cEHH+wmTZoUuU8bRUAERKDXBN5++2139dVXezF00003+SAKiB9Ca2MZ+vznP+8++tGPplZMIs+98sorDlE0aNCg1PJVRvUjIGFUv3uuGouACBScAP72rMpuk5CjioswQiApiYAIiEARCLzzzjsOKzdi6MYbb3QEP/jQhz7kRowY4QMo4C6XRbhsLFAvvviiv57WKCpCSyh3GRSVrtz3T6UXARGoIAE6F81EEVU+88wz3SOPPOJHSIOfRRdd1LuphD9LLbWUW3bZZf2HzskyyyzjP/jhK4mACIhAJwRYbJXw2r/61a/8X0tbb721twohhlZZZZVOsm7rHCxFrHeENUqiqC1kOqgFAVmMWgDSbhEQARHImgAjnY8//nj/5/LLL3f33Xdf1pftz3+FFVZwa6+9tv+stdZa/f/zHfGkJAIiIAJGgDlDiCBc5XCTQxyRCJyAZWifffbJVAxZObAU4T7HoA8fJRFIg4CEURoUlYcIiIAItEEAC89TTz3lV3FHCBGu9rHHHnPPP//8fGezZse7777bNMfRo0e7n/3sZ75D8JGPfMR/FltsMX8eHZXwB59/OhGvvfbaAp+nn37alysqLb300n61+Q022MB96lOf8p+NNtpIgikKlraJQEUJvPTSS+6aa67xrnI333xzfy0333xzt9dee3nL0Oqrr55L7Yk8Rzhu3p9rrrlmR2sc5VJQXaSUBCSMSnnbVGgREIGiE8Dn/Q9/+INfuPCPf/yjmzVr1nxFZhX29ddf31tocAHhw/9sI6RtqzU8zj//fLftttv6iE8bbrihD0uLOOomIdaefPJJL9bs/4cfftivHB9OhNhFMNExGj58uNtyyy192F0lERCBahDgXUA0ueuuu87ddddd/ZUaOXKk22OPPdyee+6ZmxiyiyOKGOThvbTxxhv7wSAlEUiTgIRRmjSVlwiIQG0J4PpG5wExxOe5557rZ8F8H4SDiQhcTlZaaaWmrAjLffbZZ0ceM2rUKHf77bc7C0+78MIL+44C10FQ8T3NDgOTqh988EF/DT5YvviO+AsmQuRuttlmPnDE0KFDHR2oJZdcMrIO2igCIlA8ArNnz+63DPGMW9p5553d3nvv7XbZZZdMAii0Q8JEEYKNQZmFFlqondN0jAgkIiBhlAiXDhYBERAB50cs7777bnfnnXd6axBCCPFgCXGAMNhiiy3838GDB3eEbcKECe6kk07ywRUsffnLX3ZYi7AOBTsGiCQrA2XD+oQbHCFxF1lkkUw6EbjlIQZhMHPmTL9g45tvvjlfXQnPu9NOO/m1SrBsKYmACBSLwO9+9zt3ww03OP4+88wzvnC46PLc4iaHKEoztHYntUcU8R7k/TJgwIBOstA5ItAWAQmjtjDpIBEQgboTmDt3rvvtb3/rfvOb37g77rijHwcWGgQQViD+Eo0pzTR16lS/BgiucriOYBUiChMWJ0LhRiVEEnONEEqMAGPBQUSxQCwiKcvEXCWsZ4gkOloPPfRQ/+WIgIdAYn7UZz7zGUfQB6VyELj//vu9SxXWSlw4CSk/efJkt8Yaa7iDDjqoHJVQKT0B5gsxV4jnkwVXbT4jrrBjxozxbnKIoqIkIs/95S9/8e+wFVdcsSjFUjkqSkDCqKI3VtUSARHonsAtt9zihRCjqRacgJFUOvW4lNDJX3nllbu/UEwOdAiYQ4TlJ2gdoiyMoCLK4sQR+0n8RVi9/PLLfv4QLih8R5TEnRtTnI4241KIoKQTdtttt81nWSOIwzbbbOMXfcTtTnOUOkKc+Um4ddJZRhDx4X+E0SmnnOJFLsKIY5SKSwCrLu+xKVOm+MESS7wPsAjtuuuufq5g0RLvKuY4YiGXpahod6ea5ZEwquZ9Va1EQAQ6IIAlhlFxOvGEocVljkQH/nOf+5zbYYcdfMcwr0Tnk5C0uMUFhRFiZ86cOX57O2FqTSRhSUJsIZCoGxGdGC1mnZE8RBLccD9EIBHZiv+DifDg22+/vV8QEtYSSnm1tPjrYFnA0mjtHjGEiyfrbGFFQiTxnb9KxSGAy9nvf/97L4QY3MHt1RLPGAM7rDPEvMCiJiLPMQhEoJosB6CKWn+VqzcEJIx6w11XFQERKAgBRk8ZScWqwTwZRMTHPvYx9+lPf9q7kyCIevWjTFhvRAvlCSc6DW+88Ya3/CSdhIxAIjEaSyAFAjXYNbhe0vzCZWv3O6KMOUq4Jk6bNs3NmDHDvffee/2nI0ixznEvcOHCjVApXwKInqA1CIGE+9yFF14YWRCOJyGcEPaciwuoUrYEeG/hwnrjjTf6D88Vgyok1iIzqxB/yxAQhbLj8kdZsZgriUBeBCSM8iKt64iACBSGAOtfMD/ikksucazhQ8KlBCsFQgi3rl4nxAsdg2ZChf3md9+JmAm62zEfiTDdzF1CnGCJwmLTSb7dsEMgYVFi7gNBLYLJxCqClbDmSvkTwHqE+Dn44IMXuDjWpfHjx/e7avH9kEMO8QJeKX0CWIGwvNraQoT5t8R6Y7j67rbbboV4nyWpPQte//nPf/bzibpdgiDJdXWsCEBAwkjtQAREoBYE6DRcfvnl7qKLLvJR1PKcK9QJ4BdeeMGHwybCXZybG8IGKxcigfp0K2JMKMGKRWexlNH5QizRIe42/6QccPcj+AQj4HQAEW6WcANk9BuXIAStUvYEcJ0jkAeDCViNwgkLEZ10szDxPwIK65FSOgQYOOB5YP4jAU4sETUO6yrPAgMHeS22mk6t/p2LLT9AZDysxXHvvTSvqbxEIExAwihMRN9FQAQqQ4CRx+uvv96LIfzt6fgzispk8d133z3zCG3dgGS+E6OlreYQ4Y6GiEIopNWRgBMfhBBBGwiTSzleeeUVP6+pV644dMiZL4HbIx1E7i+J8nBfuadMImdOglL6BLAU8UEAWUL8xM0vwu0OYRRlXUq/dNXMkQEKXH0RQ1hSLRw+z+amm27qhRBtf6uttnIf/vCHSwuBuY8MwjAYJNfL0t7GShRcwqgSt1GVEAERCBJgzgpucpdeeqkPMoBbiXXQlltuuVLAIsz1oEGDfJjuZgkBQ+Q6OkrME0rbqhN0tyPCHOILMcZn3XXXXWA9pWZlTXPf3/72N29F4kOgjMcff7w/ezrkuBAhlBCMSp0RQABhIZo0aZIXP3RYsRyyuDCJ/cw1IiBDMLGf7RzLPv4qtUeAsNRYSRFBWIUIlGLJ5gp99rOf9dbSZZddtr1MC34UliKsitSdd4qSCPSSgIRRL+nr2iIgAqkRYN4QHTgEEa4YdBr23XdfL4hwRytTYlQYYbTlllu2bQW69957vZBqZWHqhoNZkujIYKXiWrAmKAKuVYi4tIVZu+UlpC9WQYQS85NsXgsdLeaN8SE0eNbrOLVb3jIcZ1HnEDomhrAY8T8d2aDbXFR9OM5Ce0scRRFy3hprQggxxHNvgxELL7ywtwSZVcjWI4vOqZxbsRTRlqhzWQatyklapW6XgIRRu6R0nAiIQCEJXHPNNe6ss85y06dP9+WjA3zooYf6hQrLmrCGID6SuKxhGXv44Ye9CKRDlXWyzhvRo3C3o6xEuGOtEeYk9VIkUfe7777bCyVCr9tcDOZhEMABSxLzMCinUnMCJooYYED8mpUI61HYhY5jOYaPJYTyySefvIBVqflVq7sXSyth6hHv8CISYzAhflgbipDaBIFJ8g4oGzWiYuIqiCVM7q9lu3vVLa+EUXXvrWomApUlQHCAiRMnugsuuMDh3sVEXeYN8anCqOOsWbP8XJ6kYWptkjudjDwtNyaSiGaHWx9zHR588EHvwoh73+KLL97Ttsh8LVtkFrc747TBBhs43JKYtI77ncKBd3ebYIhlCEuSJdohllyEVR0TnX/ED9YgPgR+IdS+JQKnIIRY5Ji/iIQ6JOYHzps3z6+jlPQ9Vwc+qmPvCEgY9Y69riwCIpCQABaJH//4x+5//ud//Jn77LOP+9rXvua23nrrhDkV93BcS3BhYkS+k8nUjOhjCSFwQ57iCKImkPgfqxedH1vPBksWbmxFEB/MQcPljgntwchedOwtqhfh25WSEWBeEffbotIxv8gsTslyKu/RuMYRZh4BRDuj/sGENc2EOGKojlZLhCEug+uss06mrr/lbUUqeS8JSBj1kr6uLQIi0JIAQuHqq6925557rneXY3TxsMMOc2PHjvVr/FQtISj4MEeqE2FjLnVDhw7tSFilxTMokigTIo+J5NQNSw1BHHBt63XCwsUkdyxJfBCWJFwCmeCOaybud0UQdL1m1c71zdWOY6Pc7drJo0zHEEUNAYR7HILogQce8G6wlgh5bxYhBNEnP/nJMlUv1bLyTsCqzPprq666qtYoSpWuMkuLgIRRWiSVjwiIQKoEcJfDVY4PE/3XXnttd8wxx3iXnMUWWyzVaxUps9mzZ3t3GjpQnQgj6gIvEp2ytEJ4p8GIDiPCiHphFUPk0kEqiiWJOs6dO9fPTcL1jnkglhBHuN3xl3kgSvUkMGfOHC+AEEN33XWXd+UNJlzjmBtE4BQ+6623Xj1BhWptaxQxSMKgg5IIFJWAhFFR74zKJQI1JUDHA3c51h4i4dqEIGK0tQ6JKG9YwrqJnsbILGsOITzIp1OBlSVvOkpMRCd4AxH1mFOFBYlPr+ckWb3feecd73KMjYp3AAAgAElEQVSHSMKaROQ7EsLVRBKR7ogGqFQ9AjyLtE0CeTDvj8WUaRPBhAgaPny4j3hIBDnNl1mwHWD1JxQ37HjOlUSgyAQkjIp8d1Q2EagJgbC7HJ1jAikgiOq0Dg1uJnS+sY51a+nBMvPUU0/1u60VsSkh4BBt3H8mqTM/hbkHWGSYh4AA6XV0uyC3xx57zLvdIZYIsYwbHokJ5HSQ7SMrQRFbW/MyYaEmSALih3lnfBhcCKYVV1zRix/mNPKxeXPNc673XgY+iDzHMy5RVO+2UJbaSxiV5U6pnCJQUQK/+MUv3Gmnneaefvpp72Ixbtw495WvfMVHM6tbYg4OHQksEN1aeRAdMOXvwIEDuxZaWd8Lm5PEXzpRuAMyfwORgfsNVrQ8wpAnqSeWBAQSc99wr2L9KRKREQm3jNvdjjvuKNehJFBzOBbLBcIHKxBBErAKYR0KJsJkI3y22GIL7xKHWOc5UmqfAM8x4pJBjtVXX739E3WkCPSQgIRRD+Hr0iJQZwJBQUSkpvHjx7sDDzywzki8BQI3MqwkaSREBhaoIqwrlLQ+lJ0PQpHFe/nLiD2R7piX1EnEvqRlSHo8Hexp06b5SGT8NaHEIrNEvBs2bJgbMWKEj8allA8B5owROp41vpjXhkUSS2o4mQDadNNNHR+WAFDqnADPK+5zWP/rOMjVOTmd2WsCEka9vgO6vgjUjABrmpx66qnemoHLGCF9999//5pRWLC6dCSYzI2LTpqdflzUWMcHcdStFWrBUuezxaJ8MfLM2idLLbWUt4DBCZHUrdthVrXAGkEAB6xK3FtcJUlYlLjPRG0jeiB/NSG9u7uAhRHxwxxFE0L33XdfZKbwxhq0+eab+7/cA6V0CDCYwTuH9zsBZIoyXzCd2imXOhCQMKrDXVYdS0OAEU0WSLSV4wl9y7wLfsjLnhBE3/ve9/xoLa5iJ510kgRR4KYyx4EPfvhpCxiYl8Wlrlk7NysSfHCvw/UQYfTaa6/5Se8rrLBC6uyalSfpPuav4HKHSOJ/C+ZAPizKy3O+4YYb+g8Wi8GDBxdunRfeSVh3Law5ZSZ0vr2zkjJJejzCmA+WoEcffdT/z/pmNt8rmB9iE47GlAWHEUJK2RDg+cRNEYFKu2AtNSURKBsBCaOy3TGVt5IEbO0P3G1YJJG/CCJEEmKJNGHChFLWHUH03e9+13ekEETf+c53au8yF3Uj33jjDR9BjrkNaSc6LIQBh38W+add3nbzw5KESEJQ4rZDsAYsB8xnQCilLTDbLVe7x3HPsWpwb1j/hg4+7njBRGAHE0l07JlzRUe/FyHreTcdcsghC1SP99RZZ53lQ+mnkVgkFdGD+DERZEIoKn9EJQvy0hknXDZ8+I5lUSkfArxjuG8EKMEVsahW3Hxo6CplJiBhVOa7p7JXhgArxZvwYX7CmDFj+jsajMTyg3/ttdeWpr5MuiXcNhYiRvXpzJ144olun3320Q9mzF1k7ZxRo0b5Dm8WHXrc0BAQzNPJIv+YauWyOWhJevnll/08JOqICyFisEwdZMQezwzuYMyH4cP/jz/+uA9KQaJuCCaeK+Yv8dc+WS16zHuJhUrjEuIIgdfKcoRFAUsZkcriPrTTcML6gDWV+nJPgx8shUq9I2CiCGGEBVdJBMpMQMKozHdPZa8EAUZh6UxgJSIhgBBGjCbT2cDSwl8+wYSQQjDZeUWBgSDCQsSEeToxzCFCECnFE6AjT/QmXH+yFC2IBvIvustZPKnWe+ikWSKqHfOQsCYhOGiPdLCzZNy6hJ0fgXUJC4q5kPGMhS1MzOnAwoQAXn755f18Jv7yYZ4ZljTbnkQwYg2aPHly08Lvsssubo899vAC3D68x/ife4EQovMcl5iozzPA3BQCVCB+sAAxF7FOYfvj+BRxO2KdAQh7rxSxjCqTCCQhIGGUhJaOFYEcCGA9wn2OEdqohJBCLCGg6IQgPIqQgoKIzgxziCSI2rszf/7zn33HPesACYgG5rZssMEGpbKitEdxwaNMJMEWUYgbIe5quNzhbldmkRSsLRYYXJjM9Yz/qS+CkL/NEkIKDuEPlku2WYh0QluHFzdtlm94HwM7iB4+WBXsf0QQVi6+hwd/wnnoe7EIEGSBtrfEEku4lVdeuViFU2lEoEMCC3d4nk4TARHIiACCCMETl8yPP044xZ2X1fabbrrJHXXUUX4kG0F0ySWXuP322y+ry1UyXyI4MUKetSWD/JnQ/9xzz/kwulWfB2A8sRrR+eY7kcgY4SZCHOsPYUXCqkLnrqwJccGHdZOiEvMVqTNCKfiXRYARO1gswx/4sM3c99qJlIirGy6zWKUQn/ZXneaou1LubYgi3DsRs7q/5b6XKv38BCSM1CJEoEAE6MAQbjYYaAHrEKmV737e1eBHkWhUU6ZM8T+MrEvExOysO/d51zPr62HVwFWIjmQeCUsArnS4NLG+SF3ul9UTC4jNsxo5cqQXBkQ0u/POO71rKgKpSgEqaFPmiks7S5IQRbi/YXE79thj3VVXXdX09K997WvuoIMOanqMdpafAEsL4BqJhTuv91b5qakGZSHwobIUVOUUgaoSQAThHkeyAAvB8NzsL5IoYr7A17/+de+OhdWKKHOIpC9/+cu16WSn2RYZxceVLjg3Js38w3khEOgoE8K7G9eocL5l+m4iCVcxrBp08BBJiEZcz2jXLM6KeMzrvhSFHyKI+W7MXyIMOhwQ0j/84Q99SPG4hCtcnMUq7hxtLx8BgrhIFJXvvqnE7ROQMGqflY4UgUwIjBs3zs8ZwlrE3CI6H/xPstDdmVw4YaaMEp599tneunHeeee5vffe289rIPIcbllKnRGg881k+LwtN9xHop2Zq1RnpS//WXDngyjCgoY72BZbbOFDpxMdjoAHiAOsSlUVSdSL+uEOixhCGGEJoF0SIpzne+DAgV4wMjcrnAiNfd1113lmRNFTqh4BBDNthEEwRLDWKKrePVaN/k1AwRfUEkSgxwSwEtnCrswfQiQhiLASNYs6RzQ6PnkEX/jtb3/rjjnmGN9xGjZsmBdI/FXqngDWIua4tDOHo/urzZ8DEcMQBYjxvIVZ2nVJMz8TQPxlhJxOIcEHLNKbrTdVdmbvvvuunzxvgzG4EDIXi3oz/yyufryfzMWXd5TNiYQTFmUW38UK14u1ltJsB8rr3wS4r1i2X3rppUosNq77KgLNCEgYNaOjfSJQYAJ5CCPCAjNv4NZbb/WTu8844wy37777FphKuYrGCCxrv2yzzTY9EUbQomNMh3/RRReN7QiXi2q6pQ2KJCacE7CAxVixuCEeCHldpjlJ1IHgGwhxrIV0erEGmQjqNiAH+cGMTjSiUmG2022PeedmYhdhRHAdJRGoOgG50lX9Dqt+lSVg7nZZVJBR369+9at+BfkZM2a4U0891VuLJIrSpU0HEsbddka7KRUuU7iL0QFSWpCAudpxjxCPzEkaMWKEn3dDwtqL+xlig0huRUzc2xdeeMHPG6KDSzkpP1ZphAsiifql0Q7Jw6IAGp+quiAW8V6nWSaEM1ZlkkRRmmSVV5EJyGJU5LujsolABIGgG4vtxgWPTk63iXlE55xzjhdCzKsgwhRWIq0s3y3Z6POZw4K1ptcWByKPkXCjSqNzHF3b6m2lw88H4YHljWeG+RcMWqy22mo9tcJRLsrBfCHEHIKIsmHhYl8erpt0rJlDx4cy9LqdV68FZlcjLIsEImHeXbOgG9mVQDmLQG8ISBj1hruuKgKFI3D33Xf7yHIsEjlq1Cg/j4hJ1UrZEECE3nbbbW777bfPpZParBZ0lJnrhPWI+TNxc0ua5VH3fWYV4b4yQZ2w4CyeyndEki2UmiUnyoCrHyJt2WWX9fN9ECQMbCDeEL1531sTj7jWUTbW61IqNgHWryKwDu7TEkXFvlcqXfoEJIzSZ6ocRaBUBJjnwholP//5z30Uqp/85Cdu//33L1UdylhYOh9Mfid0dhESZWGhWRY8zcOaUIQ6Z1UGRAgCBJc1rIKIEwISIJRWXXXVVPkiPJjLgxgif66JCEOMkXohhqK4Yj3CnQ+3QxYZVhuLotT7bbyXiFY5aNAgby1SEoG6EZAwqtsdV31FIEDgsssuc0cffbR3t8FahCjSj2E+TYTAFozGEr0r71H8qBrSwWZtI8qCW6Zc6qIoJdtm1hKYIoxwWcR6w1/uOwMRnXK2RTb5ixDCGsN8IbNMdZpvshomOxrBaAvq4tIna0QyflkeTVtFFCGwCcmucNxZ0lbeRSYgYVTku6OyiUBGBJ555hl3+OGHu5tuuslbCH75y1/6tVuU8iNw7733+sALReqA0DnCasRcI4VaTr8tmCWJuT9Mamd9IMJe416G5bCVmOH8V1991c8d4ngGNIgoh6WIZIEi0i95ujliPcJ1Ewask6TUWwI898wDszlyvS2Nri4CvSUgYdRb/rq6CORKgA7Jj3/8Yzd+/Hgfqvfb3/62/+Qx/yHXihb8Yozy0ylk/aIiWIuCuHDLwt2pKJasgt/KjotHZ5SEQEDwYPGBO2J5iSWW6G8X1mnFDY32gjDC6sR8MPa1ElMdFzDjE6kzk/tffPFF71pHxD+l/AnwO4DAxlqURgCf/GugK4pAugQkjNLlqdxEoLAEWKASdzn8x0eOHOkmT56sH8Ie3S0mNtMhwVpXxI4tLnWk4Po2PUJV+cuaQKKiCAUsiMaf/wmigBseYqmdxVfLBAxxRPQzxB5WL8SeUn4EGKDBrZO/rMulJAIi0JiXKQgiIALlJ8BaKmPGjOl3paFDi1WIhIvE17/+dTds2DC/6OKkSZPctGnTJIp6eNvpBBLxqWjWIkNC+yGiGR1ypWwJmPsbf7EIEQwFMWRWXKJEkszVjqAFRW03SUkxKIClaOWVV/YWCxbOVcqHgFmKaE8SRfkw11XKQUAWo3LcJ5VSBGIJIIpGjx7t5x2EE2KIEUHccFiTaMKECYWJghYua12+0yHBYrTeeusVuoMrl7p8WiRWE9zpCErAnCMEKUEUbO4ZFiXEEmLoj3/8o3c7QyRVbU0gOLDWEvMfWSagSHPv8mkJ+V0FCxGssUIqAEZ+3HWlchCQMCrHfVIpRSCWAH7hdCbiEqPQ1113ndtqq63iDtH2HAlgtWNexaabblpoYQQS5h6QiJ5WFStFjrc69lKInfDiq4TxppOKQIgKZc05dGiZi4QVGJdYAqZgWSKSZBXuD3Vk4ADXOurD4sdK6RGwNsTvBe3NgnakdwXlJALlJyBhVP57qBrUmMDtt9/urUXNEqOvWJWUikGATi2dWaK+Fb0zS0eKhX+JHKYw7t21H1iyVhRR6IKLr9L5RwwlWW+IvJibg1B67rnnvMWJ55x8sCYVvV21Ikk9EOUMIlCvIs7Da1WHou2nzeAa+8gjj7iNN95YFrmi3SCVpzAEFi5MSVQQERCBxAQQRq3SnDlzfIcM1xQ6HCRGpJN2xlpdR/vbI8BIfxmsRdSGDjadKNY2wXVLHdT27rEdZQIGfosssoi3hiCIseLyMTEUZSFqdiXuC3NzyJOFOJkTRt7MRyJP1qHh+WZh2TLeM8pM5D3qOHv2bLf++ut7ty+lzgjQDgnFjQvv5ptvXso20VnNdZYIJCcgYZScmc4QgcIQwI2uVcJl4uabb3ajRo3yHVxG/ukw8UPJaDU/miQ6WbZYJNvofJV95LkVm7z3M48EgVGm+RN05AndzWgzi3KqTbRuNYgS5vUxTwt+PFcE28BSmOZzZfcCAcEzu8kmm/hrEcgAKzFBDWhv7C+bJQlxxHuK9xWh7XE9JCqfUjICJorghyhSEgERaE5ArnTN+WivCBSaACOAuDnRGYpL++67r7v44ot9xwk3LlxwEEHz5s3z0YiefPJJ/x2RhZvPWmut5fNjFJuOMMdbh4482G7CKe6a2h5NAHFBp7mMbml0tGkvVZv0H32nkm/lvrIOEcET6NAzuR3LjVk60hRE7ZSO8mBF4pnHfYpFZBFMlI17WCaBy/sGniw+jLWV95VSawLcf9ojiXmCSiIgAq0JSBi1ZqQjRKCQBK6++mp36KGH+o4PP4BRiQ44kaxwS0HcIHSi3Hbee+89717BKDcdJ0Zo6eQRHctcWcifjhUWKCaAW+eKa9NRoSNGZ6uMrjtR7LLYdscdd7gNNtjAzzEpU8cUFrSzhx9+2I86R7WhLHgVPU8bbAguvkoHlPvL81AETpTRPsxFIvAHgyBYjLEElsVFDZ4MLBA4gPcSEfyU4gkwoAUr3vkKYhHPSXtEIExAwihMRN9zIUBnHbcu/Ow32mgj33FXao8ALhFHHnmk+9WvfuV97y+//HJ36623ulNOOcWH9bWE69yFF17oO0EIn/vuu8+7oqy22mr+kHY7bYggEnkQLYo8sB7Q+SOqEZ2tDTfc0I/oIsTosNCJsTlNJpaK0lHsB5TzP9SfhTvp1JVVPDL6zP0sm1tW2reaZ4FOJ4MDNn+PAQMEiM0bSvuaaeQXFElYhxnoYFFV3hu8J3C5K7Jgp/wMxPAOhLu9y9JgU6U8EEVPPPGEH+RC/CqJgAi0T0DCqH1WOjIFAlgfDj/8cDdr1izfgaCzSDrwwAPd5MmTU7hCtbNAAB1wwAE+YtO4cePc6aef7jszJDoLCBbc4Og07LLLLgt0chAxdC7gjijF0tNNR4gfYDpVjEoSHQuxhSBChBEenNFpOo9ci/KZAOaaHEsnxwQa27opS9HvPJz40Jkrcz0Z0GAEuuid6LTbA88VdSfRZpk/hJiw9ls2sWsiiWcQd1oGOfifD8EcbIHZtDmmkR/vL0Q67sC41ins9P9R5Z3Mu5jBC0SvkgiIQDICCr6QjJeO7oIAI1hYMVh746677nLDhw/34Vjp4OP/rhRPgI7Ad77zHfeDH/zAT+KeOnWqGzly5Hwn8EO47bbb+m2PPvqonxtkoskONJcKfPURqSwWyTGdjnIjrBA9pODq6ViVKDP5MirN/4gkrmWdS+45o9Z0Lk0gVXlOE8yr0FGhDgjwoUOHtm11tPZXtr+0WwYTcCPk+cJShGggoIJZh8pWJyuvDUTwjK677rp+M88qzySutNQbCzCWsHaty3mxoMxYrAk2QXkpYxnn7aXNi3f+yy+/7K36tFElERCB5ARkMUrOTGd0SOCzn/2stxrMnDmzwxzqeRojo3vuuaebNm2a22mnnbzrHAKiWcJigwsFUanirBP469PRw1WOzi5zDfIa9bY5TYglOjRYUrCCEXqY8OKILDpjlA0LC39x0aN8ZZ3TxJwtOtdF62Q2a0dx+xC33As6YHm1mbiypL0dwcN7ikEbRD/3jeeo2eKraZehV/lRd0u8H3gmqTcWJQJJMH+qaO2X9wHvD8QcgWPqmBDwiHdEkb0768hBdRaBNAhIGKVBUXm0JHDFFVe4L37xi2769OluxIgRLY/XAf8mgIhEFPGDd9ppp7njjz++LTS40uHqQ6emVaJTgRghKAA+6QipXnZ2KTsdNEavGbUeMGCAY+0fBBSCkInuzEujw8p3BBPHY42igxB0yetlPYLc6WQSQZAR7jih2uo+FWk/vLk3COpu3TGLUC/qg/hm3hB1IjAB7a2TxVeLUJ+0ymBCCaEIH74jipnbCJ+iiCSee8r4wAMP+OAgdQrMgDDkXcggE2uOKYmACHRHQMKoO346u00C+IETVvqiiy5q8wwddu6557pjjjnGWxmuvfZaP2en3cR6ORZBrN2OOJYpOvAczwh5EecY4CpCB8gEEmXF0nXPPfd410xGuHEhoUNLfZjThDiiQ0cHnk4E341Ju2za5R53HGLO1rLJ65pxZUlrOwIWIYFLZFE6yEnqRpuwyF0W/pltuKqauK7KvUrCJe5YE0kEWYELnXGeRQYpeOZ63QbsfiLcsGy1MygUV9eybOd9xnuOe8LAlpIIiED3BCSMumeoHFoQmDhxojviiCN8R10v7xawGrvpcBKM4sorr/Ri6KqrruooshDudEOGDEnUYWHklTVPsL7gOkNnp4gCKYoiZWekH0sTYol60MnF2oYoYZTbXILsfDpPWNaop3WC+d/mR0Vdp5NtdF4IQlEF64rVn44oUfbgxjyxoljnWt0fOpMIVdoJ94Q2Q/nt3kgMNSdoAom/iCOCrcydO9cHQMCNi++9FEn2HuC5L6tob34H/r2Xdoz7IO82RZ5rh5iOEYH2CEgYtcdJR3VBgE7Hbrvt5s4555wucqnHqXTYdt55Z+8S8o1vfMOdeeaZHXcycHWi49eJWwkWJwQa4bm5f3R2ytxhZE4TnbXnn3/euwBhGYM1Qh3WTKinY8x2xFSac5rouNx0001uhx126PheFrX10zkmqETRJ3tbJ572zKR9/uKiac+GWYiKyrmo5YIr7BAjdNJ5hmbMmOFD0mO1hW8v3huUC8GLizDuZVULzMA7BeHH+x2PAiUREIH0CEgYpcdSOUUQOOOMM9zJJ5/sIwfhnqUUT4AOOgEq6GAQunyfffaJP7iNPcxpYVQRS0mniR9frAKsU8SIMOKiLJaBduuMAKRjh1sVE+4J9oB1k/oy54rQt0Tv475gicKND65Y1Wz03KxMUSPliC3mcSEeetFJbJdDp8dhcYMNbotFqh/3BvbmWkUZmT9EQAWzCBapvJ3yL8p5QZFEIADaBSJps802888N4iRP3pQHAcEACNeuQkRI7jWWTt7tWPSrJviK0pZVjnoTkDCq9/3PtPZ0RLA2HH300W78+PGZXqvsmd94440+yAKWmSlTprhhw4Z1XSV+QC2iVrcdEkZf6WDSuWAkODhPp+uCFjgDOla4weFyx+RmBBGdkbvvvtvPaYIvroaIfto7cxvCc5qw3JnI6vY+FBUVViPqRmet13W0kNN0xhG73B+sQyTK1uvyFfUepl0uhAkCFGstzwmDDTwzPEO8R/K6D5QBgcx8OMKSh5cwSLveWeUHT8QmgWhwEazDHKqsWCpfEWhGQMKoGR3t64rAcccd5375y196a5HWVIhHiYshazkhIm+++WbvhpJWYiFdAl+kYeWhw8n6SFg+cJGhg1GW+Udp8QzmQwePDhcWJ6xIDz30kO+A46ZFZxA3PaxtCErW7sLNCAsUnRu40TGks4OVyUbbsyhnHnlSfgJgYJ1sFUo+i/IgYLHs0QlmcIF7wvNkFrw02n8W5a5LnrQPPghVngkGGwjGw3cseFmLJNoF7rFY5Xkuy+Z+BjveKQxAEJKcNq4kAiKQDQEJo2y41j5Xfvx4gU+YMMGNHTu29jyiAPBjTVCKn/3sZz6E+fXXX+87zmkmhBETotPsCPADjfUEoWSTm7Pu2KTJJI+8bE4ToohODR0Z/qcziIBiQUrEJS52WFnoyNs6UgguRBTtA65l6dTTHohSlpdLHVyxxsEQCwSRuRBDNggj61AeLT3ZNUwgmUjiWeDZIHgA27C4Zvkusbk5tBFbmDpZDfI/mvcAcz6x2vPeUBIBEciWgIRRtnxrm/thhx3m7rzzTh+tSGlBAoxe4jqHCx1zibCsWcjgBY/ufAsWDRKTdNNOdPTp0GNBonNTpYhrabGCEUKHDrt1+BBNdNDYTjAIhJFFAkRU4PLDnC46/HTgcN2jc2TrNFE2szIVTTRZKGcCHGTRwaXzHFx8FcsDrp0IfxhFzfFK614qn3QJmEjivvEeQSTxLuEe4nqalbudCY158+a59dZbL5N3Y1qkeEewBAHvC+Y+KomACGRPQMIoe8a1u8K9997rhg4d6gMIEHZaaX4CiJXtttvOYc1h7tWJJ56YGSJ+VHExYgJ0VomOKi52+O/j9160znpW9W4nX+YDEOWu3fkANqcJdzSsIVhhEBnMaeKZYtSYRCcJlz1G2E0kIUw5P+jemIU4aVZvOrtMuMcylqZLnS2+igDi+cHaxrwurifLULM7Up59CBYSc/WeeOIJ73HA4ACWnbQDN9BuGKBAjHGdIgZmMPdQ/mKZVxIBEciHgIRRPpxrdZW99trLR/f6wx/+UKt6t1NZOraf/vSnvSi68MILMxeOjDgiWvCrz7KTTGeV6E9YPLAe1X3+EW2BzhfiJa05FHSQEAiwRiAgunCtQZjSrlhoE2sTggGRxL1HkHEe4on7T+cT4WSCop02m/QYEzGDBg3q2IJD+Sg38xMpO2XmuwV3kBhKelfKc7xZkigxvyMWQAZrJPe/k+UH4mpPm+IZpY0xyFCURLmoOywYWFESARHIj4CEUX6sa3ElggcQcpr5Mrvuumst6txuJem8jho1yq+tkYcosnLRucTdCPeUrJO5xDBXirWB6MDW1YKEUGQuFta6rBnYnCY6U7jfEeCBidq45CGMmcOBBQfXMzqXWKI4zgRH2MrUjYjGLZAIZOSBQI7Ki7bIAEE4IdzghrhC/CH6KK/lE5VXOA99rw6BoEjiPWYupQwOYOmJa19JCJgrHwMNBA9p17qb5BpJjkUUYelnzmGac0OTlEHHikCdCUgY1fnuZ1D37bff3nfAb7jhhgxyL2+WvRJFECMSE64ozHPJKzG3BssRI7Fcu47zj1gDimeh1x0txCrWIzqZWJQQGqyDgshgTseTTz7pBRSdTTqadMYsAATnkLiP7c7hufLKK/3cuWYJ4czixSSugVWAziDWNcQbljDmxZllS4KoGc167LO2SCQ7RBJCH+GMqOfd1s0cTfImX1zrsND0IjADZUAUmeUdYaQkAiKQP4GF87+krlhVAldccYW77bbbfNAFpf8jEBRFl156adcLtyZli386nWPrZCY9v5Pj6VzY6L8tEGtCqZP8yngO9535Qb1OQUshAohkc87oiGHNMtc7cyu66667fJh3LDeIJMQU4oX60BFlGyKKv5xvwiWJgGERUDq4uF/CClcmOqTmekheSfLrNWddP1sC1hZodzbnBoGNVRZRz3uGNoTYD1pA2ykVeZPvkCFD/BwnBjVol1lbeq1svJuJPIeldeONN+xpQZ0AACAASURBVFY47nZumo4RgYwIyGKUEdg6ZrvJJpv4keeLLrqojtWPrHOvRRGFsknGdIZ70dFkAjUdFzrfZj2qevQw6nz//ff7MOxlrSuiCUGNgKEDymg6nUXuJZ1HOpFYBm1eGSIJaxMDJLfcckvk82AbGe3/+te/7jueFm1LQqgpMu2MIGBWJHYhLGiXWFoQGASDod0mtVaTJ+9MXJ5Z6iBr6xHX432BFZd3dFnfFxG3R5tEoJQEZDEq5W0rXqEnTpzoO4KXXHJJ8QrXoxLx4/qZz3zGT5LvhaXIqm3rujDvpBcLA+JKxocO9h//+EcfCMImOuc1Ipt3E6AzVvZODqPuCBiLLrf55pv3Y0QwmdANhj9mcIRoeq0SnT+sUKRO2gDvGtynspyDgZsWFi1G8JWKSSA40IPrJaKIbYh2W28N8T548GC/rx2RxPlYWYkCyXw5UlbiyEQR7qMsAq0kAiLQewIf6n0RVIIqEPj+97/vjjrqKN/pVXLexWjMmDFu5syZDtHIWkW9TIx8YgHoZaKDPXr0aO9ih8sWUaZ6XaaseBBumDZQ1UTHEeGEtYfAHkTC23bbbb0AxoLUKtH5RBCZKEKEnHLKKe7ggw/2YqRZYv+1116bqSji+ogurtOqPM3Kqn35EjC3TmuHiG9Cx5OwyPA+xrJE8JGgtSmqlLRvfs/Ik3k/DCylmXg/IIiY27TmmmummbXyEgER6IKAhFEX8HTqvwkgioiGdcIJJwjJ/xKgg/f73//ed/YOPfTQnnPhB5jOehESnWjmqjDZnkUW6XBUSURQF+bmpLmOTxHuW7tlSDq/g3yPPvpoL6zaESE8U3zySHleK4/61Oka5prJuwaRi8ghhDyC6L777vPWIEQSwiROJNGWccfDYsQ5/M6lkbgmgVC4bpHChKdRN+UhAmUnIGFU9jvY4/Lj033GGWe4b33rW37UWMm5k046yV188cVeEJ188smFQEKUr15HRwuCwJWKdXdsHR5GZNsZxS0EzBaFwAq23nrrdeQi1iLryu6+/fbbvTDiLy5ycQkXurxd27ge11UqNwEsP4gcXO6GDRvmBgwY4C3WU6dO9SIFq2WUSLLADLjG8nvHgE43iYGgxx9/3LvrZeWi1035dK4I1J2AhFHdW0CX9cdaxAseYaTkvCA69dRT3W677eZ++tOfFgYJc4uY48MPf5ESI7k2eZ8FgZnYT2clbgS3SGWPKwsj0UymLnMd4uqWxXY6pO0m1v9CQIUTeZh4weqUppDhelxXqToEEEk295H1tBAoDMywDh8uvkRgxOXXEuKI+UkEdOD37rnnnvPWpqTJLEUMImogMSk9HS8C+RBQ8IV8OFfyKqx/8sMf/tBNmDChJ5P6iwb1nnvucV/+8pf9iDaRuTqZVJ5lnRippMOOa0lw0nKW12w3b+apWNhmJkuzuj2dFZsz0G4+RTiOEWlCBheNcV5sEOHcv2bJ9tt8IaxEWIt4dpoFVIiyGNl8I8TR2Wef7Q466CAvZPbYYw//6TbJYtQtwWKfzzuGZ5Y2ueOOO3rXXlzmiEpHyHqCOFjobizdHEsbZz/vLZYmiLI0hWuN0EJ0YakKhtAPH6fvIiACvSWgcN295V/qqx922GGOUX6L3FPqynRZeH5IiYRE4geziKOBjHBijSmSS10UdkZuCbdrYgmB0cm8lai8s95GJ4q2QEe/rsIoKWObL2R/ETiIJP5irQm61iGawhYmhBECiHOYq4R4CosstnFuMK9gOcmDxDFRFqmo6yatp44vBwGz9CJkmC/I3Ews2WuttZYXTSuttJKvCPtYO4n36Y9//GPvZtcsbbnllu6QQw4p3IBZszJrnwjUkYAsRnW86ynUedasWe7nP/+5++Uvf5lCbp1lQWcHC8ioUaN8BvYdgdJs1Lmzq8WfRejinXfe2S9SOX369EKKIkrPjzo/5La4Z3yNeruH0VhGaumMMMLK34EDB3pxVHSxwVwFOlR0rope1t7e5f+7OlYjgpVYQtwgkhAxfNhvKWp+kVmFeP7t/6C4IS+7RpQwMjc5ymDlCIujqOsWhZ/KkS4Be26xDvFhviCJ9yftiN8cAqtgaSKaHAM57UTXREAVzYsgXXLKTQSqQUBzjKpxH3OvBXOLttpqK/elL30p92tbaF8uPHv2bD8KZ2F1GbUjJHSeiVDcdMpY2Da41kueZWjnWohFfvTLMvfFXNLw7Z87d66P9NaOy0o7LLI6htFkRpbVAWqfMM9OULAgjPgetgyRY1Akha+AlcgETfA4E1nh4+0YxFRQCEUNqlBGpXoS4FnmYwM2RLZjoWpcfgnEwLpd3QyCIMzNYtmMMM8Dz0bWKfj7mvW1lL8IFJGAhFER70oHZaLD+OKLL/rRdSbZM6qVVWKC6tVXX92zgAvMaaKzwyguPxT8sNhoMfvy7PifeeaZ7vrrr3fHHnus23vvvbNCnkq+/Hjj4oeFq0wJaxGhdmnjCCR+uPO8x+2ywrWG8pXF7a/demV9HAIlKIxsTg/WG7PmWBmCx7GN5553ASkojPi/WWK/5W1CiHcH+UdZh6K2Nctf+6pHgPcnH0QSFiMGwXifIpI6XW7ARDttPWogIEjRfveyJmvupO2ItazLovxFoBcEJIx6QT3la15wwQU+DCkv6ssvv9wLlrFjx6Z8lf/LDmvR5z73ObfLLrtkdo24jOkIBUd3rQNkLjB8D4/u8oPCS97+xuWddPvdd9/tWY8cOdKHLC9DwqKGgC5bQmzYHJGnn37adyKKtv4RbIs+f6to9z0oZoJlQ4jwzNJxDFp/2B4UPeyjLVinkb/sbxV0gXeICSquiyBikIW8wmIsroxFY6ny5EvAAsMQzGbRRRft6OK0V9oibSzKUmmZ2jPQ7JiOChBzEmWSMIqBo82VJ6A5RhW4xXTOiYLGPBfSHXfckVlYZq7Dug933nlnT8iFR24RQaxuHh5JtsLR2Qm6yliHqVXHqVXlcKHYfffd/UTca665pjSuU1hfmAtV1jkwjM4y/wi3ukcffdRHiOL+Wyel1X3Lcj+uNvDtxq0my/IVMW86fMGBDspoHUUrb9DVjufZOpPs5znu9FlGBNHRJE8GVni3mPCygRauQQcxDxemIt4flak1gU6fd9q+/Z6Ff9fCV0XEB4V8eH8W321wotPnK4syKU8RyIOALEZ5UG5c44YbbvCd6Msuu8x35JZYYgm3//77z7cWwtChQx2uWfvuu68PHfqjH/3Il44AB0wAJcQni2LixkbCbY58cZvDz/mWW27xH0bSLSra7373O3+tefPm9dd011139RafThLWIuYVbb311p2cnvo5zUZz6fgQvjfY8aLDQ8eqm4TbBKKItS5+85vfuGWWWaab7HI9l/k6zz//fCFd0ZKAoE0PHjzYu7BgAWu1gn2SvDs9lgiNRXTx67Q+WZ5nHT1ET1CEcE06YtOmTXPXXXfdAlHibACkldtRsOwIG/KbPHlyvzWZ94a9B3g/kB/HhEWQiae4gZcsGSnv6hKwdkb74llo1Z6DgwNGhW20Yyyc/M9f2nTQwtoNwfAARTd56VwRKBMBWYxyvFu8/M4//3w/J4URb9bbOP744915553XXwpcstjGnBVGx6dMmeK++c1vOtzlGCmfMWOGP2+11VbzayuY28eVV17p/Z5J9913nw+MwD6sSPvtt58/h44bnQP+PvTQQ4lr/t///d/+BXzJJZckPjerE/hhCI6kmWsNo11hlzrKwAgxHaBu0gknnOBwozvrrLP8PSlTwiWNqGmElUZclD2ts846PiIU0esee+wxP4Cw2GKL5W7BY94Wk7KxGnU6glz2e5Gk/HS6bMQ87B5kLm3kxwBEOPG8m+tceF/U9yirEte3ARMbEed9EXYf4jqyFkVR1bZuCNDmaG+0waj2ZYLcrhF+Rthugwr8BjK4YPl1YuEJX4/8435Du6m3zhWBMhCQxSjHu0RYz//6r//y690MHz7cr31AyGs6qZY+85nPeGFEh3vttdd2p59+uhs/frz7/Oc/7910iIDGSxArEpYnE0bkhTWKz2c/+9n5avWTn/zEr60ybtw4/zn33HP9ApRJE9aio446yk+E71Wio0IUOhI/BFjLgm4IdGxsdDdKGHVb7t/+9rfekofFLeoHrdv88zgf8VAlywZiDwspzwPWMFwFedY6nRDdyT1goAMrr0RRe/ToEPIei3MhoiMY1Rkkd3N/S2tk3EocLgudRZ7xuHK0V1MdJQLRBPh9Crc5jkT4t/rtou2bAEIYBf/n94+2y2+l7WuWX/D86JJqqwjUi4AsRjneb3OFs0sOGzbMd+BY+2Tdddf1m7fYYov5SoRl58gjj/SfYEriCoflaeLEiW6nnXZyu+22m3fhS5oQRbjrMZ+plwnhw48JL35EIWsYmRsC+5p1qLotN51u3BwJx4wALWtCOGBhifpRLmudKDeToPngYopllYEFFonNY/4Rixyvv/76ZcZXqrJn+ZwbCAmiUjWJnhaWiKT8ljdLAwYMmG+3WXiCG9kWtT2cb9C1k+PDg3SIIo4xqyj7EUDhxG+n/X6G9+m7CNSVgIRRjnc+PJocNWrPqHMw0cm76qqr3J577tlVSXGfY7E61l4g7DFzTdpNr7/+ukMY4UJmq363e27ax9koGCLIXOj4H5HED0HQjSDNjj/uWljtYIcrJJP+y5pwwcTCUdYADK24c29GjBjhBxyIFEe7Z1tWC8TCESscoiz8jLcqq/aLgAiUn0DS9euixAjbzL00TCR8vAkafvMQPOYCyl9+A4NCiWNNSJmHCdZaEsfzv823C143fM1wmfRdBKpKQMIoxztLZ/TBBx/0ARRIM2fO9HMScAGKS7jdMc+oG2HEnCOCOuAGhuXpe9/7nnfPazcx7wnB1mtrEeVF7IQFj73kw/UJH8d+XvZR8xbC54a/n3TSSf5+4QrZS1fCcLk6+Y4opjOPyOs0zGwn183zHKxEuJ5SxyeffNI999xz3qLDdoRSmomBAwJxLLfccmlmq7yaEMhj7o+5H0W9R5oUTbtEoCWBKDc6EylRbm9h66VZhLiQiR5Elf1vx/N7hxgyQWQWJM4zERVX2HYsV3HnarsIlJmA5hjlePcQQQiTOXPm+Mn7xxxzjDv00EP9ZPG49O1vf9u/1LDYEFmOlxWT/gmb3U4iQh2BF772ta+5HXbYwU2aNMmvuTN79ux2Tne///3v3Q9/+EN34oknehFXpsSPA+s5BUfDOulQEZocZgiq//zP/ywTgtiysggw7nRRVsvYk0q4AxGINQeXOsQLdcblJc16k6/WL8qvcfA+DHbwsroygojReDqcSiKQJoGw6OA7v1fW3vgebHdYgcwqRDnMIsRzELQYBUV8UBSZ9QfhZOKJbVyHfO3/YB1NqKVZb+UlAqUg0OggKOVAoBHWua/ROeu7+OKL+xoR5foa8436GoEU+hod1P6rb7bZZn2N6HMLlObXv/51X8Ny1NcY3e9bfvnl+xpzhfruuecefxznNxpa39y5c/vP++pXv9rXEEP+e8PK09eIltXXsFb172+8VPsaVqu+Rgexfxv/NEREX2POTl8japvP/7TTTvN5N+bVzHdc2b40xGDfySef3NcQlH2NNZgSFb/Rke5ruA/2NX5M+hpzcxKdW+SDaQ+NH8MiFzGTsjXWn+prdCR83RvukX2NCH1dX6chjPoagR66zkcZtCbQWNzXv6fyTI0BkTwvp2tVmEBjzTv/G0Sboi1HJX6rOC6YGi7Bido9xzdEU1/D3dx/7Jkh33DelIff/XBSuw8T0fe6EFiIipZCwZW8kKw39JWvfMW9/PLLha3Jmmuu6Rov6/7yMRflP/7jP9x3v/vdwpY57YJhUWL9FEbhSA0h6uepzJo1y0cTrEpi7hoLpDbEeO3mxeBeh5WHhCsh1p5O3euIioj1dZtttuk4j6q0qTzqwUg5Hxv1ruo186iXrpE/Ady+cfHlnRMOmEBpsBKZpSho3WEfFhysSu1YS/n9st8wzuUcPsG82R60FJG3WZyYv0tZ83zOKI+SCBSBgIRRTnehDMIIFLjcsdYSrn5EzatT4oeAdZ7CiWiCLJTbzg9S+Nwif2fOFGHhCUpQx0SYfNxaCShiEaOSCqRXX33Vr6OEq54CLyRvRXTUcOehA0bHrFVHjGcw2OHjinTu2MZfnmFzOepkPZeoGpAfZbRgL1HHaJsItEuAdhqeM9TuubRFnoFOz2/nOrR1kubWtUNLx1SRgOYYVfGudlGnc845x49mTZ8+vYtcyncqPzhRooiasHjnmDFjKhfWlLDjrULMlu9Otl9i5vYh/olY13Cx82sgMScviRGdgCqEw1dKTsDEBgLGRquD8yrazZE5R5aHjYwjkDrJK+qadBCtsxi1X9tEIAmBbkSNtfMk10t6LO1doigpNR1fJQISRjndzV122aXQbnSGgahdBxxwgLvkkktyIlOMy1jUnrjS2Kh03P4ybsedjJDWdU+E2WaBWAQRHWCEMFagVgKJY4gySVASWYuStyLEi1lh6CyaxQjrDwMVuLWa9Yfco4QO98ssQ+y30XSeV7M+kZ8tCh1VSpuAzrWirELkQ35KIiACIiAC1SdQTx+a6t/XrmvYzahW1xfvQQZh95yoItx8881u++239wvdYiXABY25OszFYtFULBBswwqz1FJL+b9ERcM9i0520VzWPvnJTzpcwSibOvb/DnuLWx1WoBdffNHPL+M+MljA3LtGgIX5mgVzlVZZZZX+eXnLLLOM46PUmgBiJCxeEDUk5l6YhQZRYsInSpwER7Z5hqPmbZBHcO5kuHQIIq5BXogxrq0R8zAlfS8CgTzaZh7XKAJLlUEE4ghIGMWRqfF21jvab7/9akWAzhPzTZolOr10khERdJIbEc185w7Xw5deesmLJfbbmjmEV0dgIo4ef/xxt8UWW7gHHnjALwTKulB01pjjwzo7nItL1wsvvODWXXdd3zFfYokl+t28GpEMfRAIrBNs57qsm4N1A8HFdubMsI9yUQ62I3r4PypxDsKIzr2E0b8JIW5N4D700EP+XnB/6TAjmJql3Xbbze26667NDtG+/yVA+w2LGhNGJopMsBi0VmLFLE0cH+zcRc1b4lhz3zMrkQmv8HWiBJlupAjkTcDm0YXbZ9rl4HnhmbDnIu38lZ8IFJ2AhFHR71DO5WOeDS4p+++/f85X7u3lGDFuJYwOPPBALzyI3mdCgkn3pMGDB/u/bMcSw9/hw4f3V4r5PAgUOn/mooUgISoaC49yPPsRUXTMEU8mbjiGc3DdsmuwHg9lIcohFimE1WOPPeY78ggztiGcqBM/pAgttnHNRmh3v0gtC5OynY4/+VEv1jfiWpzLPlzMEF9soyOJOOOvlYUyIMhsjSvKnDSAQT+kAv0Dz6222spbAnE3RGzmlXj+ojrzaV4/j2s0K2+wfogYXFlNEHEeooiOoIklyyssUrAQ0R45zvLkmPBxUWUJHsO1w0LMzmFfuBxR+WmbCGRFgLZq7qVZXcPyteeJZxJ3VyURqBsBRaWr2x1vUV8Ww2ysaeAmTpzY4shq7eaHB9GACIhKJ5xwgkMYYSmCUZpucSaUoqw2ts/KZMeY+1vw3KhtdOjNwmXHmpDB1Y+OP6KGoAOIAawiiB5cyAhGgPBDLHEsrmVPPPGEj+BGmGpYYclCfCH82MYHcTZjxgwvvt555x3/YRvzcQhVy7XYxv9PPfWUjwpHObk29wDrG9Y5ysux/E++WNmsLggy9iEkzYoXxS/qXibdRjmOOuooLwCbpTiLER1rs1C06mhwLKIlrYhqceXlGnS0olzP4s5JeztMSDx7lIMykYKiDbHDcSZ6bF6SfWefWYfIh3MRMeFRdbbZ9fxFAil4PQsFHuTCyLm52oXP1XcRyINAL9ogbT44xy+PeuoaIlAEArIYFeEuFKQMjQVQHe5fU6ZMKUiJ8isG7oN09On8Yz2xZOtNNBbd86IIKwodd4TR0ksv7a073aZmHfq4fbY9uD9qm7nRBYUc/3MsooI6P/vss96lDwGAOx+J/xGAHGcihe1Dhw712xAzJP7fbrvt/P+W2EawERLMEFpca5NNNvEWMubmIGpgB2/CobOND4nj2c9frFrsp4xYsdhGBDmEFq5uzO+irFh1WJMJkbbqqqv66+GWiCUPMUcdEICcu8466/h9CEHKwzUQePylTBzH3DHOQayRV5Bzf0Xb/IcOho2+NhNGdOw5Dre9rBPCgpFhxEIvrCEmYKgnnS8rAyKH+sMhKtFBZJ+5+XBeJ+U36xQCChFEnjav0v7a9TnWrhdVJm0TgU4IBAV5q/N5TqPaoFk9abPkx99w+22Vd9x+ng0bdIg7RttFoIoEJIyqeFc7qBMd5DPOOMMdf/zxfiS/TonO8pFHHuk7wix6SifcfnD4cbAfGgQFbnJYN4joxlpPG264od+fpgUpa/ZBAYXLHsEkzNpk1w4LgWbfw/vIwwSZuQeyDSHCsYgPOwfrk10b4UNC9Nh+RA4J0WMJ6xT7R44c2b/NRNzo0aP9Nu4RboCInEGDBvXPt8LNDzFEWcz9j+ubxYq/CDSsZZQTwcVxHNNJomOBCKHzTselWaJznrWlKHh9RBqfToRFs3q02sezxXNlFh7qbUKI7ZSHjqCJJPhZsv85t90OIB1K8kQAwdfuhZ2PEON6JPIPilfK1UurWiuW2l8+ArRda5O0cRYUnzRpUmxF7LcofADtlueFNspf2m7Q8ho+Pul38mTQgr9KIlAnAhJGdbrbTer6/e9/37tLfetb32pyVDV3feMb3/AuZddcc40XR3yaJTr7WDromHMeHV7Y0aEvk0CijjY/ifKnYf1qxi0oyILHxW3nmCjRZecG99n/Qf5m6cHqZwkXQZIJLv5HHJKwJFk+tmCrdZ5N6PkDE6QkI669cG2jfknKmKDqsYdyTTpz1rELu+uYMIkTbHQq7RN7kcCOqJH2oAClQxkUX3YqXChrXDnaubaOEYEwAdqeCSN+O8Lt09qdtUm+h9sgQp58OIb/EU8cExTxnMfAzpAhQyLbt5WLc6OuQduPE2XhOum7CFSJwIeqVBnVpTMCuIchjBBFwU5kZ7mV6yzmwrBm084775x4tJ7OMrz4ccICQdADrE9x85SKSob5OwRxqFoKCi7+jxNZti+4PzhvKe68Vrys828jxM2ONzeY8DF0WOj40OFBQATdzsLHdvKd0WCzlnRyfqfn0KmjTuEOX7v5hTuT7Z6X5DjY5GnBS1I2HVtOAjzPPOvW7nn2gs8A7woGBlpZlzknLJwQMiaM7Jnefffd/bZmzzjPUtx+yqskAnUjIGFUtzseUV9EEXMscKOrWzriiCO81eTss8/uqOp0mrFSwG/11Vf3LmQzZ8707mkWua2jjHM8CUsJc4s6dRfLsailuhQdHOusBDvydELCHQ46QmGXFY4zi0XQshPuSE2bNq0pF64VdwydK40KN8WnnSKQGgGe6aDY5nvweWbgI6kYD+dBYYNWWPKLEz5sD793gpVtti81KMpIBApGQMKoYDck7+Iwgf28885zRF2rWyI0OZ1GIo6tvfbaXVcfUcQ8mREjRvgJ+1iQ4Ft0gUS5cblQSp8A7cs6FwgQ3McQI3RcWlk9OM5Gha1DxXc7j7yxSuEqw7aoUWauh6gizHqwA5Z+TXuTYx6iLo9r9Iaerpo3AZ5FPiTaFQMW9lyGRZOVjfdHeCCF9wfbwsLGBJG5qZIHx9g1GKyxvOx9Ye+YMAvKZ2UN79N3EagyAc0xqvLdbaNuWIuISHb44Ye3cXR1DiHa2LHHHusjy40fPz7VimGBIl/c64h+xkKuzFsixDQR1oqWCC5gkeAsIEHRytjr8tBWWq1lFHZDpeMR7HTwnU4Jri0IGTodzcRR8NygwDIWdH4sHzo+dIrCgQIQTmwjLwsyUBSBBA8EOZbKTspE3Tgv684b1zHrX6/boa5fbgK0I55HPBQYrOCdQfs1kRLVltkWFuc86wSasXXlEFj8zpB/MNkzYs8X3y3xTmlmTbL95Sau0otAcgKyGCVnVpkzpk+f7i699NJaBlz43ve+58NyE4nPJuSneWPNxQ73Oouids8993ihRAjroiUCRxTdstVLZkTPIyJhs09YGIXFDCO/1jEJ7wt2jqyejPSamAq6tNiocFh4cTyJzo6NCttfttO5CrvGhDtceTMmGlewjO1en7rCLG60u9182jmOzmYzAdtOHjpGBIwAAxRjx471bdcEiz2HtGs+PBPBZ5Pjgs8J75GDDjqof+4haw+GRRHvCZ53toefMRNidi2+2zYrZ/gdpTsoAnUhIGFUlzsdUU+sRdtvv73ba6+9IvZWdxPhuH/0ox/5BUgPO+ywzCuKqxoWpK233tpHfiPQwcMPP9y/bk/mBWjjAow8Pvnkk20cqUPaJRDXsaDDQ+coOA+ADky4Y8IxHGudGvazzTpT4XJYR8o6RMH9jFLTmQqPSJN3XH7h/NP+TseQOnVyfYRK0rkYnZYfZtwfE56d5qPz6k2Atj5mzJh+CGaJZAPti+eAtkYwnHDAF3t+2yVI3uPGjfMfLEvWdhFJ/M+zx/PD36jrcUxez1e7ddJxIpAbgYYbg1INCTTWTmBhlr5bbrmldrXfb7/9fN0b7gc9qXvDMtPX+JHse+SRR/oagqSvEcWuJ+UIXvS9997rawijnpejSgVoRITqmz179nxV4n43FlL226ZOnTrfvsYIcF/DHaZ/G8faMfwfzov87XiOayxCPF9+9oXrcRx5BPNnP3n0MjWsbIkvD4e4uibOrM0TYMf9URKBTgnw7DUsRf6Z5pkMP/+t8uX4xpISrQ7rej9tPe/nq+tCKwMRSJHAQuSVmwrThQpDYPjw4T6S2hVXXFGYMuVRkHnz5nlL0U477eRuuOGGPC4Zew0ePUbrGh09t+aaa/r70el6ObEXSbCDQBEslFq2R1GKvwAAIABJREFUtZgSVDGXQxmVtXlEQasQFh1Gb81qg+tdcD8jyoz0tuu2ZW4wXAtLko0A2/XNjWbOnDn+mly/0bnqvz7nc82wC04ukBoX4frU1eZFYb1qp+4cw+h60NJE3cgv6H7E/zYynkaduB7XUBKBXhGIs0KnWZ48rpFmeZWXCKRNQMEX0iZagvzofN19993u/PPPL0Fp0y0iazUxif60005LN+MOcsNdgs4xC8W+/vrr3sWOzhwR8hZddNEOcuzuFAJS4FJHZL2wK0d3OdfvbOv0h2veGCnu3xR2a6MTz4dzg53+cB72nWM4HnFlooh9iCTLOyx6bLsJiXaESNz1u91OPSmPuexQ7nY6ZZwXDjLBNsuHttsYnZ9PdHZbVjsfbuH7llbeykcEWhEIzxFsdXwn+/O4Rifl0jkikBcBWYzyIl2g62AVYLLmxIkTC1Sq7Ivy4IMPusGDB7s999zTXXXVVdlfMOEVCMrw/PPPu7feestbj4gQx5ykvBLCiAAMRNNT6h2BdoVRNyVEgCCqetnJR9ghhqwjZlYt6gUDrFmUkWOC86OiLDcmWDgnGHzCGFFfs85FiUHO5xqUhb/kEZ5jwXXNWtUNe50rAiIgAiJQXAIKvlDce5NJyRixxp0My0ndki1gWwRrURR7wnwTxe5Tn/qUt9jce++9fuS7VZjoqLw62fb3v//dPf74452cqnNSJNCOtajbyyEAeimKKH/QOhQOAmFCBxb8z/5myeoSJyptcntcHlyDD9fjE3e9XjOLK7+214cAbZMBgCwTz0LQzTfLaylvESgaAQmjot2RDMuDRYBIdAgERmLrlGbNmuWmTJni9t13Xz+PpsgJUUQI8S233NJ99KMfda+88op79NFH3fvvv59psS3ctKYdZopZmTcIhAWMWWno9NEpw1qDmGkEiemP2GXgeHeFhYt9D+bL/+QVl+hcmkseefI/68GQV9hVz/KQm1EcTW1vlwDtC4FNm0+aOJd2nfXvd3CgIWkZdbwIlJ2AhFHZ72CC8rNmzz/+8Y9aWosIW0pgA9YvKkNCHGFBWmGFFdzyyy/vBRIhvgnWkJVA4pqrrrqqa0SoKwOiypYxj5HaPK7R7AaFhREixSw11imj84cFlX3BEXLESfA7eSFkTFRxXQQRxzSz8JB/cA4Wxw4ZMsQXu9d8mrHTvnIToP3S9mjvSRNupXGiPWlerY5ncALx1mxwoVUe2i8CZSQgYVTGu9ZBmVlYFGsRLnThhSg7yK5Up9x2223uzjvvdAceeKCP/la2hEBizpGNVtMBfOKJJzJxsSMIBPOclHpDgE5PJx2mpKXlGnl1sKLKRqcrKEos8ILN67EOGW2eT3B0nWOCwoW6kFfQ/YdnJJh/VBnYZs8U4spG8smf8gQT+/K4L3Hl1PbqEDC3zaRWH9po3hbL8LNWnbugmohAPAEJo3g2ldqDKKJzfdxxx1WqXu1U5rvf/a63vvC3zAmLF4ERNttss34XOwRSmhYkhCPXkDtd/i3FOj7NrBxplYpr0Mnimr1IXDtYTzqJdMKiXHiw/AQ7hBzLccGRbESQCRfyiRIx4ZFv8jVxGLY2hc9HmPVSSPbiHuma2RAwayltDgEetH42uyKDAeF2yfHkxz7ad/D/Znm1u4/ryXraLi0dVxUCEkZVuZNN6sH6NITmrmPAhfvuu8/PG8BahDAse8LdjXWGVl55Ze9it/jii7sHHnjAvfnmm47gCd2mRRZZxGFdlDDqlmTy8+l8t2PlSJ5z9BlcK2wZiT4y/62Uyzp5lNMsSVYSosNFRZeLK6lZg+jo8T8p6EqH6LGOJX+DXDg+LOTirqPtItCKgA1GIO6TWIitHQbzJy+zqtK2yRMhk5aYIT97XlrVS/tFoCoEtI5RVe5kk3pgLdp0003dYYcd1uSoau4aP368j/B2wgknVK6CWMGYE4RIevvtt71AWnLJJd2gQYM6XiiWPBFY7777rrdKKXVOgA42i6syb6WVFShu1NjyYK4NHXk6Qvzlk1biGq3Kl9a12s2H8oTFUPBc9iOYojqLUdeIc0EKbo+7HnzyFKxR5de26hDgGWYQxJ7hZ555ZoHK0eYQN8F2F/eO4FmwfbRnBgwsb65l1lDyavbeiJu/1OycBQquDSJQAQKyGFXgJjarwh133OEuvfTSWlqLCD39m9/8xnewEAtVTbjYEcWOH0Xc4IhiR5jvTl3sCBcui1F3rYUOCR0U1gvj/ygXmOAV6LyEOyDWOSIPOjVmxUjTpYs2U9YRYRsp7+5OtT671b1rnYOOEIF/E+CZ5mOCnGdv9913XwAP746wEIoSUNY2eceYsLf3CHnz4X3B+yMossIXTNPKFM5b30WgbAQkjMp2xxKWl0h0n/70p91ee+2V8MzyH37qqaf6Dv63v/3t8lemRQ3MxW6VVVbxLnaLLrqoI0Q5UeySutixwCxzl5Q6I0DHh86ICZlwR8hGjIO5R1k+OC7YmaETZCPCdi4dGo4Ld6JsP9s5hhFqyhGVOF9JBEQgewLhQRKeSwRN8PllW5QYZ4AkLgXztYEOLEkmkvhr23kPBJ/5qEGZuOtouwjUgYCEUYXvMmuAsHZPHecWvfTSS+5Xv/qV/4EZOnRohe/yglUzF7vhw4f7nXPnznWPPPJI21HsEFaLLbaYrEYLom1rC4Io2LGxTosJGPaZeApmGHZnC7p20XFCFHGMjTYzqsz/5Mc1w8KH79bx4phgmey6YStVWxXUQSIgAh0RCAsjvvOcm1CxAZLwu4CL8awGBRTnMHBiQseeZcuL7/YO4f0QHKgJiqSoQZmOKqeTRKAiBCSMKnIjo6rB3KK9997bbbfddlG7M9nGSxZBZi9w/k6ePDl3d53TTz/dffDBB5WcW9TujbModrjG4WL3l7/8xT377LMt1ykiAAMBHljzSik5gXBHwzpDQeESdmGjI2SdlfAVwyO69mwF5ylwPt9JUfnQSQoLJ44lL4mjMHF9F4HsCAQHKBAu9r7g+eS3Eys/cxNxnQs+y5wX/G6DJLxf+Jg7XNhljvcCx9p2nndzx2XfwIEDfYAirsvfYAq/y7KjopxFoDgEJIyKcy9SLQkvvBkzZuRqLbJRa0z+NgrGC5vJ56NHj47ssKVa6f/N7NVXX3W/+MUv3IYbbuh22GGHLC5RmjzNxW7AgAFuueWW8y52tAt+BJsJH/YjpDTXKPmtpgNjI760f/7nwzNhI7hsD84VQiiFhYu52PAsW2eK8ywF5xyw3863+Qlck/9JnBdlMeIcCaPk91hniEAnBPiNDD5vPJ88l2YF5reTd659gtfgfRCMNsc5vEPYHvzf3j2cy/HsQxQFLUY2uMK1g9cLv+9596Q5p7ETZjpHBPImoKh0eRPP6XpYiw4//HC3ySab5HJFXvC8QO2lzF++M+LEC5kXbvCFnWWhuB5R1b7zne9keZnS5Y2L3UorreRWWGEFH8Xuscce8/dl/fXX9+s8EaTjD3/4Q9N6MYdJXJsi8p0QOhQ2/ycsSNhO58hEErlxjImYYO48P9YhspHfKCETFFXBzhN5kYfNMwqXnA5SuHzhY/RdBEQgewL8PvIs8izzzNqAil2Z72ZpNiHVrFQ82yaGOM7OYTvvA943ts3eScH3iFmTo943za6rfSJQdgISRmW/gxHlP+uss3yn98Ybb4zYm82m8EucF7t19BBISUadeDmfffbZjhDFZv5vt9QIItZsYqHSffbZp93TanWcRbFbd9113YsvvuhefvllL5Bg1ypQQ6v9tQIZU1kbFGA3HY6gALIOCB2S8D46IEHXtqDAiRIvwcnYQZeX8P/WQQpfj+15DVbEoNJmERCBEIHg+yMMh99Rfld5V7R6djkmaGG2vHj3hIVV0JptxzEQEzVYEy6TvotA1QjIla5id/Rvf/ubw1rEuj15jvQEhQ8dLlyxojpzrXDzMuZzzTXX+E5i0sS5CKuxY8cmPbVWx5uL3Sc/+UlvRWJeES6ISt0RQHyYmEeg0IZNGLFvzJgx3ooL/3CnhU5IUAy1KgmjwZzDvASSXYdtPANcGxdWrsn1woMTPCvBEeVW19N+ERCB3hPg+W4lirotJe8PiaJuKer8shJYqDFS3FfWwqvcCxI48cQT3Xnnnec7RR//+McXPCCHLXS46IQFzfJJL4uo4pPk5UxTXmuttXwH/89//rNbfPHFk1621sdz39pxpWPRXKVoAiaGrO0jVpJ0YhBGHN/JoEJ0iaK3mgALWrOij9RWERABERABEagPAbnSVehev/DCC451i0477bSeiSJwRk30DrvxZIH9d7/7nV/Y9JhjjpEo6gAwVgWl7ghEuakkyRGh0omlNMk1OJZy5mlRTlo+HS8CIiACIiACvSAgYdQL6hldExc6XKOOO+64jK4QnS2j5OPGjXOTJk3ynS1EUNBtJ+gyhHiZOXPmfBmtt956qcwHYm4V82coi1I8AeYJIYKYX0SUOqLTzZs3z88xUuo9gTwESx7X6D1JlUAEREAEREAEkhGQMErGq7BHP//88z7owE9/+tPcy4gQeuONN7wLEK5viCIb9eYvwshc4ggBTWCIYGKNnW4THftbb73Vr9u06qqrdptdJc5///333TvvvOM+8pGPeObLLrusd298/fXX3eDBg91bb73lVlxxRbfUUku5ESNGeKGkJAIiIAIiIAIiIAJ1JSBhVJE7f9lll3lhcthhh+VeI4QQ1zZ3Odx0bGI3I9PBeUIHHXSQ45N2OvPMM32WdbMW/fOf//QR5RA7CKEllljCLw7IXCsCYJC4B0svvbRfxwjRSGhu0jLLLOP/KmVDAEsq7JPMMbKSEBQhHCwh7VIyaEEZNc8obbLKTwREQAREoKwEJIzKeudC5Z47d67bcsste1IbOn7hTlzSMNvdFBxr1cUXX+w233xzN3z48G6yKuy5CCCsPQifZ5991oubhRde2D3wwAP+vsPgox/9qJ9bBYcll1zSu8tZSsMqV1g4BSyYBSBBeCQVRoiirIMvgAzRxmBGMLx3AVGqSCIgAiIgAiKQGwEJo9xQZ3sh5tZgNSh7YhT7mWee8dYO/k6bNs2vZ0QnLi6dc845fv2db37zm3GHlGb7Bx984N3fEDq4viF8sPbY/d1ggw28OELoIIRggwBi0dZu09ChQ71rXbOE251SawJYYbCUNmu3UbnQ/vmE1xmJOjaNbQxoUNYkYcLTuK7yEAEREAEREIEiEpAwKuJd6aBMBxxwgB9lnjFjhhs2bFgHORTnFFzDCKRAahVNHiGBMMJFbK+99ipOJZqUhDoh/Ah8gJh98803vSCZPn26n/vDfuqFVehTn/qUnwPE/2xHBCW1QDQpit9lUQS5Fp+4ZKGo4/ZnuR0LTJ5WyG7rEhWZsZ08e1FPRFin5W2nTjpGBERABERABMpCQOsYleVOtVFORvxxuZoyZYpftLMO6de//rX7whe+4E4//XT3rW99q1BVJtobAoc1lZjfQyhxvmMJeuqpp7wIwdLF4qpYhfibd8JaYCGig3PBwuXAjc8WFA3vy+O7LThYlgVJscRgLeKDBYhBi3asQByHSAkm6o5Fh3zIl/2IVJvb1y1/8uPT7P53ew2dLwIiIAIiIAJlIPChMhRSZWyPwPXXX++tCr2aa9ReKdM96oILLvBzbQ4//PB0M06QG65vb7/9tl9YFtc33OBuu+02x7wvttGxxdIzYMAAhyvcaqut5kaNGuWWX355b+nCDa4XogjrhHXeW03AR5D00mKDlcwsGwluTc8ONXEDV7iF5+AlKZhZkRBNJmS5H90soBy8fpQYS1I+HSsCIiACIiACVSEgV7qq3MlGPVZZZRU/J+emm26qUK3iq/Lkk0+6qVOnehc6rDBZpn/961/eGof1hzk+iJ9XXnnFR3978MEH/V/CYmP5Ye7PyJEjvWAjIYRIbC9Sskn3dIxbJawViKheJhMZ7ZS3l+VEsPAxMQTnYEI0sd+YNhOlHGP7ycfqHhRFiCSzTCHCwq6WQRdI/pdlqJetQ9cWAREQAREoMgEJoyLfnQ7K9vGPf9yv5VPVhEvR5MmTfcfyueee89U84ogjUqsurm2EvcbCg+vbyiuv7INAsG2dddbxwohOKAKIeUGIpO22226+CHCpFSajjGB39tlnu+uuu87PdaIj3axzTkc+ShTRybYOOnnwv1l2sii6XSOLvNPME15B8QanIF+EjAU7gCvHhsWMlSfInXzDlifLh/wRUYiesLuhRcgjL/LoxTymNPkqLxEQAREQARHIioCEUVZklW/qBBgNRxSFE9HosBzFdS7Dx/MdC9Brr73mXQ/5nwVimfPDQqjkg3sbFh4+G264oQ+UgFjC/a3sifqdfPLJvhM9adKkltUJd/Q5gW3kQ6eeD/cGYcT2LCOclSFQQJgXQsTmCCFOguImSZsNsrVrBK1I5G3WKfaTuDeIpihh2/LG6wAREAEREAERqBkBzTGq2Q0va3XpxEeJIupDZzA8ks52BM9bb73lrT2PPvqoe/75573wITgFogiLExYi1v7ZaKONvBWIdZCYB0SwhDXXXNPvw0UuuCZQWRkGy411gXlOlhA1uGEecsghbVfPggmYRYR70EwUcU2sVAiFuERe48eP7+/gh49LIiTC5+b1PSyMuC7lDooVtpklJ1inoIscx9Du+di5bINjXDJhRJ6WL+LIzuPc8Fwx8m4nMETcNbVdBERABERABKpCQBajqtzJitejWWeaqiOaiEpHNL5HHnnEz7dC/BAZbt111/XChrDXCB5c4qpi/en0ttOBDlsRWBMJ98F2UrCzTT7tChZcPWfPnh17CfLBihcWCLEnFHAHwiPI1iLJBQVJ0AXOXBCpCly5NyZUzEWOfQgYngPyNv7h6tt5YaFDnpzLdYPX43wTtuG89F0EREAEREAE6kZAFqO63fGS1pd5Pq3SrFmz/Ho/BEJAGG266aZuiy228J12xBGd8sUWW6xy1p9WXKL2Bzvf7IdRWChFnRfeFrSO0OFuZs2Iyz84JybumPB1i/w9bL1EEAXnGCFQaM8DBw70c72CiXODPBA4nIsQ4n/yCooiE1LkAXvbx30xKxP3moGDMWPGeKtg2KoXvIdF5qqyiYAIiIAIiEDWBCSMsias/HMjQOjrJZdccr7w11VzgUsLZhL3KTrk5qLF9RFAuB2S6GSbtSj4P/mPHj26/9Os3EmsQ0nK3eyavdyHgCGqIfPisI4FrW0IwzDvZmVFCHE+7PmERRnncu8YFOB6fILCCpGGq56SCIiACIiACIiAc3KlUysoBYEhQ4a4OXPmxJaVjl+wwxd7oHZ4AkHrQiskcA0HDBg7dqzviGPdoHNtgiU4r6Xd+5GkY07+QSHRquxF3N+qvmY1CrvDxdUlPGeI44Lso/ZzDIIJlu1eJ+762i4CIiACIiACVSEgYVSVO1nxetCZxBUoLkWNlMcdW+ftJmCSdIZNiATnpgQ79606+q14I6ziOu/BczmuWVjxVtcp0/482jNtIEk7KBM/lVUEREAEREAEOiEgV7pOqOmc3AnQIY4LLY31otvOee4V6sEFEUVwQmBEdbzZhwjiL8cGE9vC6+MkrQJ5YmXCUhG8X0E3PcqGNYvj+D+YOK4dAZW0XDpeBERABERABERABCCwUMPHvU8oRKAsBKzTbPNS6ChXYcJ+XvxtHlC7bm7BciFM4N7Jud3WD1HGvS67G123HHS+CIiACIiACIhAdgQkjLJjq5xFQAREQAREQAREQAREQARKQkCudCW5USqmCIiACIiACIiACIiACIhAdgQkjLJjq5xFoHIEkoTWTrPyvbpumnVQXiIgAiIgAiIgAsUmIGFU7Puj0olAYQgQDKFXAoU5Rr26dmFugAoiAiIgAiIgAiKQKQEJo0zxKnMRqAYBIsoR+KBXgS6IoheMZFcNqqqFCIiACIiACIhAkQhIGBXpbqgsIlBQAoiSXq4hZAu7hsOIFxSXiiUCIiACIiACIlBCAhJGJbxpKrII5EkAMdJLUWR1JVx3eG2jPDnoWiIgAiIgAiIgAtUmsHC1q6faiYAIhAnYWlAbb7yxn7eDNabZ2kQIo/B+zmN7eN4PeWUlonDjo+xKIiACIiACIiACIpAFAQmjLKgqTxEoKAGEDMLF5gxhhQmLnnDROZY5PsEUtCIttNBC7umnn3a2eGz4fH0XAREQAREQAREQgTIQkCtdGe6SyigCKRFA4CCGsOyQsMAgjBBMCJs4a48db8XA2mTnDxkyxAdlIO+48zn2/vvv9wEUmgVRQHCRN3/jkqxGcWS0XQREQAREQAREoBsCEkbd0NO5IlAyAgiOoIUIsYKoMRe4sGtcXPUsOl04v7jj2Y7gaWWdYn9YhIXz7FVkvHA59F0EREAEREAERKBaBCSMqnU/VRsRaEoA4WPWHixErYQKmSFUwlYaBBUpKIz4v11hZYUk37CbXtMKaKcIiIAIiIAIiIAIZERAwigjsMpWBIpIADc6RA0ihkVT2xFGHGNCiDpxLmKGbSaE+IvIMWsP0eM4jw/XjEtYf5rtD5+XVHiFz9d3ERABERABERABEYgjIGEUR0bbRaCCBBBDJHNHa0cYMW8I65IlEzs2L4ntiKSgwPn/7Z09dBpZtu/3XeuuJTLRkZhITGQmMhOZjowj48g4ajpqHDWOjKOmM2WDo8FR42hwNDgaFA1KXqN1g0Eveei+4KFo0E0GRYMyiOr996kqqCqq+BCSQOp/9WK1VR/n43fO2Wfv87GP/lsNKP0tc7HtzmAtwx2MY9n7fE4CJEACJEACJEAC6xCgYbQOLb5LAg+YgBo3alx4HScEjZKwGRl3D5L3mRo+rlGlhtMqBpaiC4a/zlI6NbDWmV16wEXFpJMACZAACZAACWyBwH9YuLYQL6MkARK4ZwJqwLx48cLMFv3666/GyHANIzVYvMvl9L7XCYI+V29y7ozTTZLuXXqnYbtxe401nWFyL02nO7Ol99WIomF0E/L8hgRIgARIgARIYBUCNIxWocR3SIAEpoZTcJbpPtAsciV+H/EzDhIgARIgARIggcdPgIbR4y9j5pAESIAESIAESIAESIAESGAJAe4xWgKIj0mABEiABEiABEiABEiABB4/ARpGj7+MmUMSuFUC3r1ItxrwksC2Fe+SZPExCZAACZAACZDAIyFAw+iRFCSzQQKrEFAnCnoGkTpiuMm1zPX2TcJc9Rt12LCJ84dV4+F7JEACJEACJEACv00CNIx+m+XOXP9GCahxcVOjSL3C6W8bzhe0uFy34V7Pdb/RYmS2SYAESIAESIAE7oAADaM7gMogSWCXCeiStFXPHfLmQ911b9tdtsa/zVmrXS5Xpo0ESIAESIAESGAzAjSMNuPHr0ngwRHQGRc1jHRZXfDA1UWZ0Xfdc4W87+kskhuO99+Lwtrkmcal8fAiARIgARIgARIggdskQMPoNmkyLBJ4AATcw1J15khngfSMoGWXfhNmFOnsjRopOpOj+3/0PT1E9i4vXcpHRwx3SZhhkwAJkAAJkMBvkwANo99muTPXv1ECalCogaP7jHTWSP8dNDLU0AmbkdH9Sd7LPXRVw9Fn+tN/v3//findZfuEdFYomC43UBpGS/HyBRIgARIgARIggRsQoGF0A2j8hAQeKgE1SLz7hNT48DpT0OfqtS5oGIUZMmpcucaSu29JDa1F+5DU4FGDSmeqoi4NS2efomayggZaVDi8TwIkQAIkQAIkQALrEKBhtA4tvksCD5yAGjg6q+Ne7t/uHiF9FuZ1zvtNEIG7V8ldardo39IqXvE0/kXxBePn3yRAAiRAAiRAAiRwGwRoGN0GRYZBAg+EgHeGSGdk1ABRYyVsRiiYpaDB4+4vcsPR990ZIfffp6en4v6C33vD13QEZ6mC8bt/a1rDjLeo93mfBEiABEiABEiABFYh8J+rvMR3SIAEHj6B4LI5d4ZHl63p8rlFlxouwcNV1TD65ptvzGdq9KgB9M9//nO6lE7v/frrr9NgDw8Pp0vvgnGtYpi532g+Fi3XC4bNv0mABEiABEiABEhgFQI0jFahxHdI4BEQ0FkW774d/XudM4F0ZklndVyDymvMuIbK8+fPp6T0vUV7iTZB6qZhkzD4LQmQAAmQAAmQAAl4CXApHesDCZDAlIAaO+fn52b2J7i0bVXX3otw6izS8fGxXF5emjjc5XXepXQarz7TdAQ906kht2x2a1H8fEYCJEACJEACJEACUQT+w8IV9ZD3SYAESMBLQA0Tna3ZhnMENZLUcKNhxDpJAiRAAiRAAiRwFwRoGN0FVYZJAo+YQHCv0n1ldVvx3lf+GA8JkAAJkAAJkMB2CdAw2i5/xk4CJEACJEACJEACJEACJLADBLjHaAcKgUkgARIgARIgARIgARIgARLYLgEaRtvlz9hJgARIgARIgARIgARIgAR2gAANox0oBCaBBEiABEiABEiABEiABEhguwRoGG2XP2MnARIgARIgARIgARIgARLYAQI0jHagEJgEEiABEiABEiABEiABEiCB7RKgYbRd/oydBEiABEiABEiABEiABEhgBwjQMNqBQmASSIAESIAESIAESIAESIAEtkuAhtF2+TN2EiABEiABEiABEiABEiCBHSBAw2gHCoFJIAESIAESIAESIAESIAES2C4BGkbb5c/YSYAESIAESIAESIAESIAEdoAADaMdKAQmgQRIgARIgARIgARIgARIYLsEaBhtlz9jJwESIAESIAESIAESIAES2AECNIx2oBCYBBIgARIgARIgARIgARIgge0SoGG0Xf6MnQRIgARIgARIgARIgARIYAcI0DDagUJgEkiABEiABEiABEiABEiABLZLgIbRdvkzdhIgARIgARIgARIgARIggR0gQMNoBwqBSSABEiABEiABEiABEiABEtguARpG2+XP2EmABEiABEiABEiABEiABHaAAA2jHSgEJoE7MTx+AAAgAElEQVQESIAESIAESIAESIAESGC7BGgYbZc/YycBEiABEiABEiABEiABEtgBAjSMdqAQmAQSIAESIAESIAESIAESIIHtEqBhtF3+jJ0ESIAESIAESIAESIAESGAHCNAw2oFCYBJIgARIgARIgARIgARIgAS2S4CG0Xb5M3YSIAESIAESIAESIAESIIEdIEDDaAcKgUkgARIgARIgARIgARIgARLYLgEaRtvlz9hJgARIgARIgARIgARIgAR2gAANox0oBCaBBEiABEiABEiABEiABEhguwRoGG2XP2MnARIgARIgARIgARIgARLYAQI0jHagEJgEEiABEiABEiABEiABEiCB7RKgYbRd/oydBEiABEiABEiABEiABEhgBwjQMNqBQtg0CZPLMzk5vZSJL6CJXJ2fyMmF/+6mcfH7WyJwfSEnXz7J55OLQLmtEv5ELs9O5JRluwosvhNJ4H5lxPXFqZycXUWmZvrg+lJOj8/kcvmbd/sG2ujpyZmskOK7TccOhb5yGe5Qmm87KeH97W3HwvBI4EouLq4fHIbJ5YVcPnC1824Mo/Mv8u77N/LqxQt58eKVvHr1Cv9/Id9++63v9+LtF6fzu5azz+/k3afTR9UJTS5P5NO7t/L5YoO6feGwfPXCcNTfm+/fys8nV3J5/EHefPs7iSczkjvqyNhEcy2nH7+XF7+PSyKdk0prsEHkO/7p5EK+fHgHFltXodYDdflVvs/mpFQ5kvJRW4Yrf30lxz+/kW9RtslMTo46o5W/5Iu7QGBX5NxiGXF1+lHeffgs57fUJ198eSuv/viNxFNZKdR70QUBWff21R/ld/GkZAs1WfBmdBi38GRy/hnp+IN8E09JttiQ/i2EuSwIleXfv0E/iT7y3fEtgV8W6RrPVy7DNcJc9Or16Sd5qzoE+j3VHb5fsROdnP480zvQT374ent9Q3h/uygX2392+fmF/PHD2Q0G37af9vtOwW3Lvc3SfynHb/NS6aKPn5zJp7ffyxunLbht4sULvz797bcv5O0t1vfw9F/J2Zefnbb5Rt5++CRfAx3FuF+TQv6zPOhxW+sOr0EtYwGulSi0rMFoZI3HY2s0GlqDftdq13JWMtewRhr/uGUVYoJ3M1ZtcIcJupegR1a/XbfK+YyVjMeQp4RV6m4escsynm9ZY29w47bNLuuwdJ91y1YS7NPV/uaRLw1hjDId+NO19Js1Xxj2rb6pLJ6rW7ISyKOkKlZvzeC29vqoY5XSOasx1BSswG08tPqDQMbbRSuOfGfrJpDtX2Fp3FaqdiktQQa7JudCZcTQamRVFsesfNMnaYK5We/vYd3Kos7Gi+0l343s+GMFq7Xkzdt7HNYOB1YtAw6JktW5vYgWhjRGeaQkZVV2VZitXIYLs7nGw5HVyseNDiGJotVeWh21zLTPvcM+Iaq/XSNX9/dqz6qkwCKOthTsO+8vETsYU1h7vyO5d6Pcj6xOOW2lICt9xebIa0mVrQ50AtWnx9Crh8OB1e+2rWouaeWad1jQo65VzWh7jFuJZNJKxLWf0F/CytW9eubY6lXSVrLQtHZEQ1m7FO5mxsgxNROJhMTw7xj+f7i/L3t7e7K/fyCHT57Jy0JeUri/r+/u5aTabki91ZDSofPxA/3fxaeCFBsjyVbb0jlK31ouXJbxZEL2vKHuxSSmkIMXbobdDr628d/XJ/LhDwlBI15j5mO9WC8+v5LfYVasGhxCfnYkrWZdms2KPF0vyK29fXZUlPo4K5kDTcIe2sKhvzy9KTv7KN8mk5KvD/zpva+yXYVSVBpX+fa239mltITlbdfkXGg9OsCsTlvqjbbU8j5JE5aj1e+FxhX2OeRW/F4klx15pPyKkKthSb6le+N+XwbxjGRTtxTgbQezchneVsT7kkjGJa71YdiUamPx7M/ktCq1vtPvJVOSvK1keMOJ6m/vIq4Nw5yc1KQ5ALtRS2pNLgg1OCPb+x3JvRuU4fVJWQrNlFRrL2392A0DdRoDwVCok5I8tPXpPejVBweH8uTZS8zS4HncaNR3c/Wgo48K0h7+W/71z3/Kv/49kl49j3Y2lHa5LF+nVWxPnh7VJd8rSeHz4jZ7NwndPNQ7NYwWauZxGEOVjJMDKIjPf5AfXz+JVhI3z+u9hPDk/d/lH3/9SV4/2Q83WG6ainvUFdZL4kgwo4FBg4TdaNf7eKW3VWEYjhG+kQre60CeffejfPf0DoVBMMqN/r6QTmeAAZcVjdZRX/rDMZSDuYxvlIpb/XiX0rhLaQmF/DDk3N6Tl/LjD8/l8BbtolAcO3Hz7uXXatmcSLfTxaKJnGR+E9xXowItUDLFkqSxULxTrcpZ5N6FK2liMDJTKdyNQbRqcnfmPfCo9STfqEomBna1mpzvTNq2mZDo9r4bcu9cquWGxAoVeR2i1ixSAxO5qkxV6jtAfD0cSfqoKi/NoK5e+/L0xwYG0LB+ZdyVtnfgeu+ZlMsp6R5V5GT3VgYvpXO3htGC6K9HcXkK42F2XcvF2UXEHiPdJAwHAycncnZxFbledoLNsvp8/tKwz+fDxibfs3Pnff33yen8uvrJpZyfwomBpi1SKM/HuLt3sHH//Dx0c9zk6lzOQ3bNXcO5w+nJMfif4XmglqOlamONYTRxvj9HuV2cib1/UMtQN18HnURcIz24fwxnAiiLUMRmSiwmiTCp4C3DIPTJlV12i+oN3tF6Z+cK7yOfx1oPblDY18ir4XQKvqHCYCSQLatfHraLPjKbgY+P5RRsQ6N1P75RXdb6ouUTka+FafS2afwb7ejUs5l07fbq5gNlbrdJtGlvhVmUFq0L03L20vTWUe/9xWk3b96IZ7ic87LQdniKOnt67m8r7v1oGWjLSRex3W4Xy8zwehXFxGRa7HqONoJKHtpe3UC1bWqZR7Xr8MgDdxHfpfKw4/M+nFxfydXVtVxf68/+t/0c31xdee4vSOVC+eWNbVXZsAYfX0570u6MJZ3L+EeJV2Lkf2m5HNL3V6jf5rU1ynBpe1ilPwjJcLosFVW+Bg05ipr5OK9JbVSSoxzeW3BNrqAjoG1pXxa5r30VvWAuDq1zl3J5iXo398zTx4XIIK2rV6YO42f+7dbXQD2e3p+LYP7GhY7uF6X8HX7KpA8F9nRBO0AIRgahvWr/FypjjAzFc+2rQ/WhdXUL7/tghLBPELe/XJz7kX0qEu4tL+3PVV8z74fkd2F7XyT3wGdadxboB8pwqlOq8xa731xVndBZvkYfq0QK662BUVk4PnwqPpXaUyvcvmBWbuAK2bnOFUtXpJoPWmv7ksokEYwD1hPgYb4o6WFLqlFtdp3I7/vdtRffrfHBuJW3gMtKloObbPpWtVizzKrEfsMqZtMWlF4LZnJgXfnY6jfLVj5fsqr1hlUrY19SLGYlskWrUqlYlXLZKte71n83ilY2lTBxxQqeNeyDplXKuWHnLXfJ/LBVtnJp+33B3pxBr2ZlE/Z6yVzDWaM5Hlitct7KadyNhlUvZ7GmMmnlsZZy6VJnh9GwrnuswvcY9esIL5bAOv7VVmFGs+xYRV3rGdxj1Ktgvbpnj9GwbVXMvifNZ9qabT0aWd1a3sok7bXcKe8Cd11TmstYuXLdajabVr1atNLYNzXdtwRGnXre7GWKpXJWsVi0CoWiVf/fPateyFopd49VByyLKVM+kixbdm2wy7aQy1rZbMZKJey14QnsO/OuVh31mlbZrGtNWBmEXSwWrEKlbf3fQBn6V9YOrXYlZ6VRT7TeNKolK5uMWfFUwar37DfHvbpVyKYsTN6YfQTtQQt7f+IWlm141s32VirrcQ91OJOx8uWa1WjUrUohjXATVhbpnJXuyOo1S1Za44unUa/yqNd5q9qNqE3DrtUopu32g3yYfBdqlnm9Y++tyqIudrFXLxFHus0ePS0H5CVYpW5Yl8dom4VM1irWGih/tL+S1tnkbM9cVBr/O9Cmhx2rYsrQbp/9NdurK3LGg7ZVK2atdDpjZdJJs89K4tivpfsSI9LyX/9VQx5Qzr66pyH662jR3UgSlEeBtJu03IRnMFwnU34WLci7gpFxbnkmsjWrh7Xk/Yb3fsxK5utW3606465VM+1N64DKm6HVrWTsPGu+lXsiZ1U6gYoRlBHjQLsNiO1Rt2YVc2jnlTrqedUqaHvJlOfqm5ZTtQjZmUOdzSNdaNtYTm3Ss3yP0dhq5tEGtS9Q9tmE0x41H5D9OeVhCsHqNUpWxpHb8UzJqjW7zpr8IeRSxcoltW1nrXKj51+r77DXcgyVX2aPz9CqZ9eTDavycaP3/X9Qg2xIWmV3f9G4b7XqNatWq1q11mobb1eSQ8F6GFa/le46ZbisPQTrVWh/EErF6pZT9r6Jnu6/QnmE7ifFXqRCCn0pZHvf7vfE3b/sBot81oo59DXZmeyIpaxyZ9ZzrKQXYMfZfH+LuljLou9AWks1axYk9rJA5mexp7QEGdqoV9FPQO9IZKxS09Uj0FYb5an+kcpXrXrLfTbAXmz7WSyVtyrT++GsZnfH6CIcHlqWTn8RzzfD2wHqehv6TSZfsero55toO/lU3IpBp3AlxqBVtDKZAnQhPG/UrDL2s8R075K2xXV1i+D7Xbt/wBJaW0dAuZS0E8N+XO039L7KMO03y+1ZeQ0COsAQ+l4B/bx51/ziVroIfWIqJyPae1BfCcg9C3pQDbpTpgA+0AVr0AtTkNFJlJVb1v0mdFBHp4yB86BTgWya9csSz1hwpLCk4MZWu4h+UmVfmFrg7PGTdM0KSoROuWi5qqsvEqS9rvU+l4f+kLMyKFcsxbPS0Jd9ut6SlC163ClBRsfcfdPeN3tWGTI4lq0/uL1GsijDmz5zlfl4Fp1WvW4qVR3/r1VQsdDhu4U7huCsplGRA4bRsJmH0pezvLbDsJEzHWwsA6W3VrOaprJhE1qvaqX1vtcw0gyggzGbaGMzw8jcHmAzsCqTqawxvGoqtFCxC6ZGDqx6DhvM4OhgpkpA+EFIyxqbYxcZRirwtaNf1TnCxoaRyfTQahZUQfUaRuaBNcKGfnVk4G0s3TKEX6CDGXdgqDoeMoZdlGetYDqsGBSTBsq30WhZxvYYYwNhKWmEU0YbZQWdPIyVZLpqGz5o5Lk0nCZMBcAQnZsqQQmr6O6yNcoBHFlkNc0wSqsaPn7oWE2q3TL0GYV2vLEM4vEKFwjjIhqpxLNW3al4YwgNzPYaQZGBoe7qjSMYyhhkw/2s4yTBRBd+9VGP4lpvvELP6SjBJVnq2J0RFPdmA8qahqsGWqtltVptm9VcyGpEob2UbMMoVVRFFPmG4mfqo9vRobMso4O1Y0Ynis5QO4SEz9vHTeuytpuYlQl4Q+lXEafpOBancdamYQQW8E2tblXzKStl0naD9tqtwnjPwpD0cNYyBYMm6toiXmN1eKFlPzXKHeD4zm6HqHOuYaSpm8qjsLTflKc3XO8A0IwF3EhahWrHGph6O0Bb1faADjiT89/P2/dztgcPOzOe9pbOYeCo2bOGGg7avDHITX2Gc5u+k3f9X9Aw0nvTcPwDOiPUuVQMcmPWYKEQNUx9judmm2wHMOxSqN9NOxNO2lA3HUNtZcMIcjZbKFl1tHUt8REUdde41jblNu1h0+4P5gff1LCBgtUJ0zDsZC2UX65htKJsWJWPh77vnyP0a7G442AA9bpagUI3stuguDIz6mO9v6ocwquL6zdq3lpluGJ7WNYfRORtahihPBo5VXrjAVmLD2FUZrAhvWc4hBtGtqHg6c2hL2S0/8/4lczFeoFGEDSM0N+gbidDBghc/cUM3EyvgdUw7dc/KDrAQKkOHGbmHOpgEDmTXs8hx6gJw8brrMJxwhDh3GqsfT+MHJ9zixEcYhUch06mncNo9xkMaF95xOHma03dYqaLxK0k9MOGM2A56lamg4dZ6GVz9wMyfFpeyYxVKNWtDjw0QaJCXsAQzqv+gf7QcWKwsL1HyD0ollYlHTP9lleSjOAIQQc5MYsy1WHgjc3C5hAMxqgcbztyXAeAMail95e2Y6ecQgwfg9k1jHSQC/1pHYO+Rp+uwnhO+vVbfd3oMemCPXA4rX+or47umQ+1vqYvrvYPOCMpJsAHlWNe0jqDXMG6tVrIW33rfgwjKM1qxNScXxUzPylYkTN54RklnOJwvAIFK9MIjVQrJDpjX0GMm/b9oGGE8FpmBDJYcRwBF8/PPLaMBvCehwoFgwzOrqEsBYra8YI23wmHl+Eiw0iVw+FgGFKZwsO6FcMIQfdg7MwbRnjQtw3LmWHkeGkJ8oex0vVaHCiPLBp9mMLjetJLV1xtzJvngdWDEPNdTqcWNBbVQAv3WBjspDSDdscY5rHNLlekFXXELllnVFjrQKCoexU1MtAJTyV/WLmMMLoNoy3Uc5XbGaU9yigMMVXQA51xWMjmXrtgOss5TzOOYZQKehxEZ5YPdPY3r8voGBGWVwk1aRp0/cZcVBpB2Iz8w6CdKacja+AqzOu017EaNhC+03o0I4Z9B84M5AJekDRmcCRoGOGTIRQSY4x7DCNtm1FpvzlPTXOYnNPbjuwKjnI77SE4OOEqf3F/oqEf6gx1zPIZTA4qd0Ap4TEqQg0jLWITjtcwcupyJjjy58gIHS3UpmwUqFi4t0SV2xFywkmi8z+XPdqNT6HEY1X41MBzDQj9AnWjmNCyLdmzqW5gqqT4+hh/LNO/IuXXOrJhRT4RSdB60caAlZbzEEZRrQovruZdzLRX0HcGZ/rmwllXDkXX73XLcJ32sLg/mMuUuTEzjLSoS46C6aw2MW9gdgSrEaZ1PsIwGvV7Hn3D/s7WCwIGgWv4hOgFdgo9fQ5myppYHZApueVlv2EnCzIL9TIOj7xzCiMG9IxR5pVHKru1bgf7Wxh9OXj3CvSUnojm/9mvZqx0wLWhyjntS8J0l341bQYB/TYZBpu6TgM0ek9w8A/NvecM1DlJWF23sD+w3w8aXGNbX0O/4TfE3PsY3PHJBac8YBjPLb7ATLoZ+IRBOM1bZHsPk3u4Z7gFB5Lt+tMpBgao4HXU9L9zMySOLHHl5HyR2Xem3wc8DLvvewyjslefhoyYM4xGOhCMgU13FNgTZ3ifF5WoRffHpn3Goc94x8u8X3R0BizCIF8U8raf/Scs6ju/4pmi/PT+mSceqNJvO0vidfZjjIcy8r7peMZJpJIhe1oigoREiLxw1k/GXTa5fyiHePGk3UGcMRn12nIicbN/RsZjGcPJgAY1HA4jg1v9wZ4cHE53sa3+2aZvmv06q1wHkkonRWoVSf++Lfl8TrDkDb/n8uzJKt/rviN9L4H9xO4H3jwfYo8ZHuua4G5Xuj38usod38HhwE2vi3Yb544k4Nlpnu1+Lo+NqC2s5W/hjJSXMq2RsYTEAxukUukUynqAYtd1yvO7p0z6Jh1p6VlC6bSk5hL8VPJYy1Pt96TVvZb3UYt/575b/QaW0Plf3k8KmgW8ENltRuty98Z1OS0ZDIk16zlJ9nKSx7lLdvk/k6cROEJTHktLbrqbfF8ONVHLrkAVnXTgeXAQk1x2vuI9/+mnZaEtfr6oOYSk/QTnhd2ZbIg7XjrdFCfSZiP5IOo+ZJL3stsbzrgKqfsH+QLqflvacGQykucy3zpmIc2JiMuOdPp4nu5LB3sAEkYGQx7iv/5YZeMQMlHkqtuQzighpfSi0BcXx/RpDF7FgnVlH/UQQ7Kt9kDU34vJxN5zKRdT0qi2pN6pybOXduU8qzehd7aNPN/oWkU2rMgnevMQ9hfhvJJEqgvvYXkp//TaKZ9Def2nX+S1k4HLr2+l1Aj0PUl4uazFbiaHQur31Zf1ynAd+bK4P1heSnvPK1LONKTcrUr1pCR/0bK+wr6jblaO6ovr3P6Tpwa/7rXodnvSQ59j6jRaw0ibUVCmhegFvhRi70Q5W5ZWvCb9v7rl5Xmj15QWigpujeeClsOcFLC8Rfu7zoXY/SnqdrmQFKyQwF6gsvzy3E7QOTxEpsqt6KoTxDY5lVp9iHGCvnz9ajLoXPBkhn/1GzU5Pvqrb2N/Ev1XYgyPwKk/SMv081nze/bMaT3JrKTjkMHwetbP5SSXw/MMfs+e+fM2JziCiQv8bd7H3mFfN7bnOBvC0HTofUidMPUAsjIVLEM4ACigD671B9JTFIurSIizrCvpNLsyjuXsftWX/D3JFLISbzSl0+7K5IfX8LpsvxCDvPb3zNClsIdBsIfQpD24TccNdzwyz2Po1xd1S5JAXXn/3ifbMqO3cFCFgJw0XLe1z1RZHC0BF8YRKKqwP6/RBoudLGTP+0idQPOibWwDlS4s6ju/tyXnCykplHNLvJiloAgloRXATefMD6BcNhrSTRShr0/V2g0hBT2EXaOT115X+/uhDKBI9Hs96atnNMlLDY0A+zo2jPNhfP6sCsUf5ZQcdeDhpiKlfFZS8W/k2w8nIRtMF+QptAXiALMPr+QP6YLAeY6kcmUIdPU+ZMgvCGzxo2HfVhxCQ9iHsqXtFNLHLuEFYYWmOfA+DBCjoIVGhrkSNSz18WhpbAsSssmjTeryobxvYeMkOkPpQemswvVmDh1o4g/y9ngdF5wAGeyw1szSEB3bdggG074JzzUzra9rHQyrh1H3F0Wxb3fW8UQy0Gkv+sh5NhwYV/wYKpLRAPIQxn5P5eJgKOlKSzrdlpSeYMAI5WSaQliaV4hm+Stw34zNV9AcfErT01IZRt9QWnAxbrYTT06k3s1g03mUBrI8poVvBPO3Ip/IMC+hIA9x5EEShLt1KX/4Mu8ECK6Deh0tBf+VwGDV4Y3lEDISaJvrleEN20OQXzBTkX8fSrGSR/0dwgNdwxwOf1atybhUEceOiPzy+uyzfP/H30sWRsYQLtELcCd8lEsY2R0uvoN6QTDoMfQE9CNtuFb+dD7nhGSCthEtsxKSNJ6EvArjnjwva/83wGCUW49PpTGAofxsdQF63apJG4pzNo62in5n+kPIec2vuu5u+R1U7b2sS7tRwjdwuwzDqVLEACI8of7uxUfbC+DBd9Js16SIaa4+jlWplosY7IIc+f0b+QLD7rYvMxAdckXdD3l1eiturCttW4veino2NF5h9QqrI3tJe+BKLeuw575QI/IUfMe8poPwUUmKuJ+Gu2wUz/QawBI0YYSjjAhl9duT84+SrwiO2flFni8QszhrySTCVPcHdN3LjNE8D7itNdMFiy4IilpbMD4k1VJKfl/D2S8pCCsYRZ3ed5HeN+arMEb816plMXhTtksxmf9Rfrwt+2tRVnf12R5GLP/8d/zU2xU6b8y01GsYUayhTIr/T/709KYJn8jZB4z+1mM4m+j/yE9uVVjPSUpo5PYZKBjBHiCwqBZ7W2dcoJ6o2IVGiO4sakAKQgGdzHauDevywUv56W/4mVk9lH+rCeO1JY1SVYqvIRBvJVPBxjnfXu0OcYRBCqhCL6NHwJYnJxgXvgi5FR3OhjyjA777J1eYZRnHJYdRztXVLCdZMEZMfx3PSv7HHyIHXu3RQTWc0PaCLvTBeS3UoUSujDIaC7q0PizACKpIoVWX1tVryWOGcVDA7NHaGQ2NdPnNFflEBXTd6Ugfo9LV2p/l5d6VfHmVkmwR89V/+84zuHwgr3+BLA4L5Pr25NB6ZXj/7WH/NQboUhgs7WDW6Dgmw3ZaKt0lMuHqCwZ1StLLtczsjqvHDcNYrnovUZBGsyxH2bzAMZPkxh1p//R02rb2nDqB5fIwmmb3p8EbFQMz/N6u4Yl6kDuSIowbU4+7TRkXq2vMel5Ko9aXXLUp78MsxexY2qkK+u+GXPzwk8w0MJw788Mv8it+ZkYN9bFZh3e0zpFUGkX59ccD2X/2Xv7yD/zU6xqe6/mB9SYGRCrHUvgbZktW5XbP75mBUng1zC2pIuHJwgCMEXwYGNLKEhGGni25wDYIDzrs7h4WLEKh0IFUlZXrhGnPiMLL38VQEk/chOrAfogsDot7jXuTi8+SL3al1PqbvF4yC2cGhTEzHVzcskZ0W3l1SzNGK+Z1L4kpXQiPAmaJ/vE3+etf/ip/+VOUUWR3DmPM9vj064n6V1+nS96TtDlhD6MnzbO5kaAVU77Sa9dwl6gLte7zcg0HM9PhuS47ugTNe03k5OMn5+wDLIF7+lK+e/8LRpSKGH8Jmc4OQRxyy4mgj+UwiC2B6XivfQzpE91ZRYfmTXUaI6hxpK/b6szPal1BaUO+UzhceJlZ7qcT8dc+lham8WzQllbIyNmgO8A2pZwUsxt2G6tlPSSRm9TlM/n48dQOc+9Anjz/Tn7889+kiUXbkNxiVhh6r7XTuHp7PchkzFLFbqO+2lkcc2lRBQ4BQDYMfem+xlLDXgi3qFub8IwK837uXzQb0oe71crr5XVxDp9ZSqMF0JS2f7DZl/hkBjOKaHsdLGMLzilemeWrK15zCXC+u2pLs5fEaH8xYJztS76cx3IgHErbOJFGIyalYoQWE5WEqDij3vfeX5FPeFA4vwirEGbnFx1gZS4oDuxZOoFS//HLAuga6C3KofXKcL32sAniGbunUj7KQcYPpFEoSR/D1mHnvXhZTzCj2YHcT2MJ2EzZnM12hZfL8rsxDBz90m1LOY06X0Ff9tEzc4SleFiJLCPU+25QViLt/T4M/GxBcj7FEoeMmnrcQT0+lkYrKaWXq6vHkzN1WZ6XUphRpNmB4VXKIlE9vOdx3X319aN8cRrs/uEzeQmj6S9YlpqPYUZMR/yhCH88drSq/Sfy7PWP8qe/tqWaQZjaFzioVtctXLZRNSLq/vIy8b2Bw1wx+QZ5UZrv70OimL91KJlMEkGqvjLfBieYOR/AuM0XFMRtXCnBCn70U07bDwlyPo3el7o4G7RtdKgUZLHqQe1Ga+6YGnsWJxD45Bpu4kMiDN66/CrFQkvyzaZ8t1TEXgn2FOvSGUmvXo2DMW7l7zs1jFzLd2wW8S668NMtvfwAABAoSURBVNxMueHnESITrLXNwCgawNwc4XyAK/Xxb34hBsWecwBotyZHp/b5H9cXJ/KxpN9r3BqHNw3uH/PTloc4XAsbIWVQhyDBsiGvXLs+h1/6ECU4LHcjk2+EP5qTjHL+8x+wLAl7YT6tFlg0y4h8ONOxkFu+S+PEOJu0qrpcA+nCbMDZl3dS6UCH19ROGWGpAIR+LyAPRuiwx6mCFN2ZNCxrUaVz1O/IGcK70jOS3BhnSQvgcZbCwKBonGprxHenn+UDDugLJNd8Z88YDDBjhbLQM2ymBx3M531PRxTRI41aFTkKnAZ4gWUC3WRZ6rrux1xOnZhl2rk/+9+CR3jpiVnakYDgrJe/+JXB62MsWRhLDgfrzTpux6D01cO5KGc3nBmpPg5+vILgupiea2PX2VABFwju5nV5KD0obP7aCWUCFnUS5xNMtw0tSGNYm54mb532+hQjs/mE6dALb4/9Z3BNJrP2GZkW1DczFAfFuXZmzj6aXJ/L1w8laY4RbphsCJFHmvab89Svtb5p4ePnFQl6z9we+wdKou67EIPvm/to25WPcjpttxhFPH6H0b00FK2fMHbtufC9iTbY6PSmXu7/97JSwax9DPzKBSyv8XagOLvm5PjUdL66B6SC9RzwXCm5t1/t80hQb8+PfxasYIIKYWd/Xho68TmRmno9xt68r+ezwY1rGOrFusSxiqAeovjtPS8LthpJ7wjLnZMlmTtuwxuF99+R8mtWJlGfmqLUa0U+4eHo/iL/+UVmlNWZ1b5o9TB+tGRY9gZyKLQealbWLMO12kOwXoUD8d01/V5AfzjIY1YlqdUzI5Wyr0ajI3L6c8/ypj1nz0YXytyFVj5t+z9XcF5MWAJmiXT/5X8r8Hz/ufwZMyhqJHTVOMIyc9P0DgpypIkcNqRSv/C168lpTeqDjBzVinOTEHsvMfNp6nFR2pmCv72GJXd670paRw34mFr0DQyvUtb0/w3MKrkiAh7cMHgc0Iox0zWMZaSoUy04OLuNDVm+djuBcYf9hLlibmpsrq5b2Imetp9AvqLuL8w+9p2Vv3g4Q158Kh7JGIZCw2s5R7Z3hB5SP59VKpKFgdipBg8q1dm5tsQxM33kyqOJE0BkBpz+PzIjB5LFcnUZYqlymJHizCRFfT45bUoXeyITeGEPWxNUHo7bJcl9cPtM6BDHH6VcD1b8S/mMZZa6TP7DWVTouH91Im9zRenGk9KvV+Tdu3fy7u1befv99/L992/kzRv0O95KMumZfXwJ7EdTe+9BXXfi/aGPs2dKBXN2DGDAixC8xhThEnvugBXEru/CzzpWyeFd+J2Hi0Z1QWwuuBvMmLM5wn44J6bqnlthvz5owL23Ccf+xdMlqwU3c03j5hNerRB2BWdCDFoVq5THWTP6HtyxqltYPevC6z1Gz7nJp+z0x1M4yycPX/BpeKAp1CM9cNipwBkaOHugBJ/x9rki8KAHN8MF5L/qOY+iD89PMCngTSfoeilQIgtY+vOhjMtwb/k/iL+Es1vUYwryh7iL3rMV9GwAnD8yZYr852o45wPnRCT1fZyxUMC5RV14n6vCo0ssaZ/BYM6xgf/+DFxr+o/VcN1s28yT8MIzhDeYOvzqu+UfT+eRhqrvvJMR3H5jFNqkIwb3lsUG3E6r21UnzSW4u3Sdq6orduhbTl3KGJfB4Xl3StA5k0jPPShU6lar3UJ68lY2Pzt3QM/80HqXdOsd6kDVUz+nrubhDlrTEl1K8N6Ds4Q0nATY1Jpww93U816ynrMq9KiHKupExoLQsuscyqq+7FwDeNwqO3VQvYTl6v/LauNMpkLGLj8tm2K54amPric8tI2C1gWbx43qsuMBJ6Fu8et6doWez5Sxst5zIbSqzqURbTfQplNwH12qNG037Z7qvWp7NZ+Ys1Ccs7DgOlZdwOdxxlYKbbztumwKS4sT3xiyJOuRJTGcb6Pn+riesuJwa1qqdazhCmm/Ec8IOddrlFFX3HylrbxpJ3pGT/B+xXLFouvFTut3vlSe3re9YMITljn3A25wIa9SSbyj+fKJFZVRITLiyt9uEzi3ZFb39awTt/4mcEZYzsrh/LE0zmgpe+W6nlGSdc6Ig8yNxeJWpoL44VEup23NpLkSkCHexKmLarjqxvlztnxGPlIpKwXPpk2f/31fhswfdlmGeZCaf3d2Z15+/c+NZMOKfIJJQTvL4NwWrxOxMVwWZ+D+udmq4dyY1c5SM+7zV5BDK7XNNctwaXtYoT8IYsFeK1sv0DarfVKgzqi3sOTUu6h+Df44Zy8/7dv03Du3X7ddZGP8DrI3bqVy6AfQxsyZMapzOLJymV4Q3uc4KYd8qqHeq56RzDtn8aEPbRiZhTZQrFoNHNHQRLvLZWfn6QXzrX9rPY6tclSE87F9flYS8Ti6hu8sJTeGvtUse9q8pjOnbUpPf0C6UQfNOXw4x0jP/cvhiAD33DO7L1ROeo6PnnOEo03Q9vM1v/5lzh1aRbf4H3hehH6YUW+SSHMS5/QZ/Us9Mgbul+rOcRfIhu0p1n5fZbUt9h2vdChXlRexRAoyCfIiq3Jxzh+gqSf2sSDR+opf7qFl4UyiLNKq/W251sQxGw30hVkrh3MKXb1Az9ss4jgKI7Og7+aLHnmNWA1jwxx5nQrykNJXHQh5sY+NcZ/b8jrv6nQq21PQRTM4Wwk/k98kvFoi7oTXU6nqQdNvcDYe2lEJxzjoGZrG6+jUEyuOcVEX8uY4hpA0mVvO0QGav8gfjtbxuk+ES/OkiSesHKLi2Y37/6HJQEZv+cIo7mRv6qXDBK4ju3DbMb+QI/ju7O/JJaZCKzrKX5ACaiYK3gxvDvsYNcQ62e4wDa9fnj0q+lxnFPrwG5VIypNDZ/4OcbuJUYN2z4wye9IX/NskWC+Mtp51pIvRExhQWGL3fMHepulHdl694WtImoaQe3uuKxPP5/5/BvmYwGyWQaYR9zVufzx6wjP2xYxjGJh8IgdOodgptMM3vPSka3jv6WMJ0hjjEMl0Rp6HelezwxtiD0Py6aEZQQrGGfzb5FH3ruga4OTsxGY3DZ4is3GgXC90/V8iJU80wVF5n8Lzlh02X2LU4qmbUTtyfx0NhqfJw8+gmf5jGvj8P7D2+rSDqXWTRHh0g+s+/+yxU+fcL0Pimw9U7+ip6XA+gFGap1qfA9/NcY1M9Pp1WU/67mITp3FGgvhTWNLmZ+imOJBGALPLz23twb89OV2lvfpeV69S6owBs0ApLD1w6tvslWBaPB9j9uJSvbKhLaeeHDiyyFO4tnBYMe3r8gy2Y+dvZeWRi9PyDNaPub+Rr0DlvPr8LTzHD6TU7ctRfODUxafiikEPich6FKxP8+3QqeeQHYkUlpI+dTn6QjenxPe1McRT8D5pt4RpWMG8+BMm1xgt3d/fM7PDfciHMepeelpe/nh8f138LN9iqWfr79H7oMK/DsqvQFmFpDeymblyYAmfWToQ9+VYDoKFpPJuFJ/1YeEJn7+7ihxC4mf9QXTbXK8MF7eH+XoV7JOCWVlSBug7LscH/ro9V1m9YUygGqDPA9MUZIYtmaa9ja0PaDkv0guC9WDub8iXUUwOff2M6wlPnTepx8iMZCLazIwAljVhOufgwN+DBAl565BX3wqydt/T+9qvT3WwafpVZnbhUEUdRuieWKQRHum8sV9fQO72scRLmzT64HQWcmUakDdlq+kW/qKalZM/7YE64OmIZ9+fyttvstJIN2TYysI5DOSFT76HUQu29zB9Zaoy2gG4utAAs0KJJJbYPZure966M8fam3aEGIrOxHQtx2+SUow1ZODZDxfs983fwayF6tda/y5kgL2ZifRTo+tdfX6BfqIvpc6/4P1wFoi/TIKB24z0bpjOavIbeHb2Ac5OumXpY2/a4XxwO33njgyjW8jz5Wd5gU1zsHyw+W8e6/WXV5LABrBC+99w23kL8TEIEiABEnjABGaGETq836DTmItPr6SaatlunB9wOTLpJEACqxKYGUajX39Yy2HBqjHc+3uXn6D7NiXf+4e8n1d9N05OlGG0ccDeALCV4A0cfSSxFPjP9+YF5/ZycKd7jDZJ5gRnCHVHsMyz4TXDbPSDtwvsMeNFAiRAAiTwWyaA81vq8FBWfhk9FvtbxsO8kwAJPBACh++xfymBvdBf5hwnPIwcXMtJuSLDMrxHPkCjSBnvrGG0ZzafYgN46C60S7iT7EoSm+uKy/alPoyaxFSSAAmQwEYEIvf8bhTqLn9sLwPTta5nR0fSL5TX2Ky+y/li2kiABNYl4PpOWPe7XXz/8IemNOFevvT5YheTtzBNV8c4kzKG87Q87usXfrCDD3fWMJKnR/CckpYeDhh999X2JIVF59jzcywf3+AMgRTO08ECSY4P7mCtYpJIgATumcA1zqzAfj3sEhj2r+bXn99zau4luvMjuBH/vfzxjymcXQRvdD9wlOxeuDMSEtgRAhPdN6YWETy56dGFj+fC2VLvm3BJn7zlLF2bfZvaTwwGd9NPHGRr0vzl5YNe1ri7e4yc6jC5OJVmC2cBwKHCWE88x+bFfD4vz8J3/t1yJWJwJEACJLDjBK5P5VOlCafxOIhVPdRg6miM87Mqf3r94Da9rkN6cvEVLo87OD4gJ6XSa3nCUbJ18PFdEnjQBC6+fpBaB+dB2ULPeKBPF2vy4wNdvnXnhXHxRT7UunCk9dvqJ27CdecNo5tkit+QAAmQAAmQAAmQAAmQAAmQwDoEdncp3Tq54LskQAIkQAIkQAIkQAIkQAIksAEBGkYbwOOnJEACJEACJEACJEACJEACj4MADaPHUY7MBQmQAAmQAAmQAAmQAAmQwAYEaBhtAI+fkgAJkAAJkAAJkAAJkAAJPA4CNIweRzkyFyRAAiRAAiRAAiRAAiRAAhsQoGG0ATx+SgIkQAIkQAIkQAIkQAIk8DgI0DB6HOXIXJAACZAACZAACZAACZAACWxAgIbRBvD4KQmQAAmQAAmQAAmQAAmQwOMgQMPocZQjc0ECJEACJEACJEACJEACJLABARpGG8DjpyRAAiRAAiRAAiRAAiRAAo+DAA2jx1GOzAUJkAAJkAAJkAAJkAAJkMAGBGgYbQCPn5IACZAACZAACZAACZAACTwOAjSMHkc5MhckQAIkQAIkQAIkQAIkQAIbEKBhtAE8fkoCJEACJEACJEACJEACJPA4CNAwehzlyFyQAAmQAAmQAAmQAAmQAAlsQICG0Qbw+CkJkAAJkAAJkAAJkAAJkMDjIEDD6HGUI3NBAiRAAiRAAiRAAiRAAiSwAQEaRhvA46ckQAIkQAIkQAIkQAIkQAKPgwANo8dRjswFCZAACZAACZAACZAACZDABgRoGG0Aj5+SAAmQAAmQAAmQAAmQAAk8DgI0jB5HOTIXJEACJEACJEACJEACJEACGxCgYbQBPH5KAiRAAiRAAiRAAiRAAiTwOAjQMHoc5chckAAJkAAJkAAJkAAJkAAJbECAhtEG8PgpCZAACZAACZAACZAACZDA4yDw/wEvvqni037hrAAAAABJRU5ErkJggg=="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimization with proper integration\n",
    "So far, we represented the trajectory with a collection of configurations. The regularization $|| q-q'||^2$ can be interpreted as a regularization on the velocity, and extended to higher level differenciation. This is the approach explored in eg the excellent tutorial by [Marc Toussaint (2017)](https://argmin.lis.tu-berlin.de/papers/17-toussaint-Newton.pdf). Going in that direction, you might want to represent your trajectory by a proper polynomials, which would make your derivatives more stable, then leading to collocation.\n",
    "\n",
    "We will go another way. The drawback of $|| q-q'||^2$ and further derivatives is that it is not Markovian any more (you might formulated a 2nd-order Markovian). For Casadi and the backend solver IpOpt behind it, it makes no importance as they are not able to exploit the Markovian structure. But other efficient solvers, in particular based on Belman recursivity, are able to exploit this structure for improved computation speed.\n",
    "\n",
    "![image.png](attachment:image.png)\n",
    "\n",
    "\n",
    "Marc Toussaint (2017). A tutorial on Newton methods for constrained trajectory optimization and relations to SLAM, Gaussian Process smoothing, optimal control, and probabilistic inference. Geometric and numerical foundations of movements, 361-392."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We will then represent the trajectory by a list $X=[x_0...x_T]$ of state points combining configuration and velocity $x=[q,v]$, and a list of accelerations $A = [a_0 ... a_{T-1}$ (all of them in configuration space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_acc_helpers\n",
    "# --- Casadi helpers\n",
    "cmodel = cpin.Model(model)\n",
    "cdata = cmodel.createData()\n",
    "\n",
    "nq = model.nq\n",
    "nv = model.nv\n",
    "nx = nq + nv\n",
    "ndx = 2 * nv\n",
    "cx = casadi.SX.sym(\"x\", nx, 1)\n",
    "cdx = casadi.SX.sym(\"dx\", nv * 2, 1)\n",
    "cq = cx[:nq]\n",
    "cv = cx[nq:]\n",
    "caq = casadi.SX.sym(\"a\", nv, 1)\n",
    "\n",
    "# Compute kinematics casadi graphs\n",
    "cpin.forwardKinematics(cmodel, cdata, cq, cv, caq)\n",
    "cpin.updateFramePlacements(cmodel, cdata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to redefine the forward geometry helper for $x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_acc_error\n",
    "# Sym graph for the operational error\n",
    "error_tool = casadi.Function(\n",
    "    \"etool3\", [cx], [cdata.oMf[endEffector_ID].translation - Mtarget.translation]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We must then enforce these variables to be derivatives of each other. This is done by introducing a numerical integration, for example:\n",
    "$$ v_{t+1} = v_t + a_t \\Delta T$$\n",
    "$$ q_{t+1} = q_t \\oplus v_{t+1} \\Delta T$$\n",
    "\n",
    "For now, the configuration space is simple, and we don't have to worry on the sum in $\\mathcal{Q}$, more on this later: for now, it is a regular plus. \n",
    "\n",
    "Here we use the basic Euler integration scheme, but any other scheme (Runge-Kutta, etc) can be chosen. We will suppose that it is explicit, ie $x'$ is a litteral of $x,a$. Roboticists don't pay much attention to the integration scheme, in part because they are lazy, but also because nobody shown a convincing proof of the superiority of better integrators for our setups. \n",
    "\n",
    "The optimization problem then is\n",
    "\n",
    "\n",
    "Decide: \n",
    "- $X = [x_0 ... x_T]$ the state trajectory\n",
    "- $A = [a_0 ... a_{T-1}]$ the acceleration\n",
    "\n",
    "Minimizing:  \\sum_{t=0}^{T-1} v_t^2 + || e(q) ||^2\n",
    "\n",
    "Subject to:\n",
    "- $q_0 = robot.q0$\n",
    "- $\\forall t=0..T-1, $\n",
    "$$\\quad\\quad\\quad\\quad v_{t+1} = v_t + a_t \\Delta T$$\n",
    "$$\\quad\\quad\\quad\\quad q_{t+1} = q_t \\oplus v_{t+1} \\Delta T$$\n",
    "\n",
    "The following hyperparameters are useful:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_acc_hyper\n",
    "T = 50\n",
    "DT = 0.002\n",
    "w_vel = 0.1\n",
    "w_conf = 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "This integration scheme is summarized with the following casadi helper.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_acc_cnext\n",
    "# Sym graph for the integration operation x' = [ q+vDT+aDT**2, v+aDT ]\n",
    "cnext = casadi.Function(\n",
    "    \"next\",\n",
    "    [cx, caq],\n",
    "    [\n",
    "        casadi.vertcat(\n",
    "            cpin.integrate(cmodel, cx[:nq], cx[nq:] * DT + caq * DT**2),\n",
    "            cx[nq:] + caq * DT,\n",
    "        )\n",
    "    ],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From that, you can define the trajectory optimization problem. Follow the steps:\n",
    "\n",
    "1. Define the variables $X$ and $A$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opti = casadi.Opti()\n",
    "var_xs = [opti.variable(nx) for t in range(T + 1)]\n",
    "var_as = [opti.variable(nv) for t in range(T)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Compute the cost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%do_not_load tp2/generated/trajopt_acc_ocp2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Enforce the initial condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%do_not_load tp2/generated/trajopt_acc_ocp3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Enforce the integration scheme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%do_not_load tp2/generated/trajopt_acc_integration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. Solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_acc_ocp5\n",
    "### SOLVE\n",
    "opti.minimize(totalcost)\n",
    "opti.solver(\"ipopt\")  # set numerical backend\n",
    "opti.callback(lambda i: displayScene(opti.debug.value(var_xs[-1][:nq])))\n",
    "\n",
    "# Caution: in case the solver does not converge, we are picking the candidate values\n",
    "# at the last iteration in opti.debug, and they are NO guarantee of what they mean.\n",
    "try:\n",
    "    sol = opti.solve_limited()\n",
    "    sol_xs = [opti.value(var_x) for var_x in var_xs]\n",
    "    sol_as = [opti.value(var_a) for var_a in var_as]\n",
    "except:\n",
    "    print(\"ERROR in convergence, plotting debug info.\")\n",
    "    sol_xs = [opti.debug.value(var_x) for var_x in var_xs]\n",
    "    sol_as = [opti.debug.value(var_a) for var_a in var_as]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_acc_ocp6\n",
    "print(\"***** Display the resulting trajectory ...\")\n",
    "displayScene(robot.q0, 1)\n",
    "displayTraj([x[:nq] for x in sol_xs], DT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Enforcing a contact constraint\n",
    "\n",
    "The second problem we defined is actually very similar to the first one. It is a simple reformulation, but the optimum is the same. But now we can nicely enforce acceleration constraints, or set acceleration costs, which is very useful for dynamics.\n",
    "\n",
    "As a first step, let's enforce a contact constraint. \n",
    "\n",
    "Rigid (bilateral) contacts lead to a motion constraint imposing a geometry quantity to remain constant: either the position of a contact point in a 3D contact, or the relative placements of the contact bodies for 6D contact, etc. \n",
    "\n",
    "### Working with a floating robot (ie on a Lie group)\n",
    "\n",
    "But first, for that to be meaningul, let's do it on a humanoid robot.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_acc_talos\n",
    "robot = robex.load(\"talos_legs\")\n",
    "# Open the viewer\n",
    "viz = MeshcatVisualizer(robot)\n",
    "viz.display(robot.q0)\n",
    "\n",
    "# The pinocchio model is what we are really interested by.\n",
    "model = robot.model\n",
    "data = model.createData()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This robot has a more complex configuration structure, with the rotation of its basis represented by a quaternion (dimension 4) subject to normalization constraint, and the velocity of its basis with an angular velocity (dimension 3). Then $n_q>n_v$, and we cannot trivally sum $q \\oplus v \\Delta t$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_acc_helpers\n",
    "# --- Casadi helpers\n",
    "cmodel = cpin.Model(model)\n",
    "cdata = cmodel.createData()\n",
    "\n",
    "nq = model.nq\n",
    "nv = model.nv\n",
    "nx = nq + nv\n",
    "ndx = 2 * nv\n",
    "cx = casadi.SX.sym(\"x\", nx, 1)\n",
    "cdx = casadi.SX.sym(\"dx\", nv * 2, 1)\n",
    "cq = cx[:nq]\n",
    "cv = cx[nq:]\n",
    "caq = casadi.SX.sym(\"a\", nv, 1)\n",
    "\n",
    "# Compute kinematics casadi graphs\n",
    "cpin.forwardKinematics(cmodel, cdata, cq, cv, caq)\n",
    "cpin.updateFramePlacements(cmodel, cdata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This $\\oplus$ operation is the integration on the configuration space, formally a Lie group, and is implemented in pinocchio with pin.integrate(model,q,dq). The following Casadi helper synthesizes this information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_acc_integrate\n",
    "# Sym graph for the integration operation x,dx -> x(+)dx = [model.integrate(q,dq),v+dv]\n",
    "cintegrate = casadi.Function(\n",
    "    \"integrate\",\n",
    "    [cx, cdx],\n",
    "    [casadi.vertcat(cpin.integrate(cmodel, cx[:nq], cdx[:nv]), cx[nq:] + cdx[nv:])],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can set up the optimization variables to contain a configuration $q$, but in that case you must enforce the normalization constraint $q[4:7] == 1$. More efficiently, we represent the decision variable to be $dq \\in R^{n_v}$ a small displacement from the robot reference configuration robot.q0. Then $q = robot.q0 \\oplus dq$ becomes a function of the decision variable, and no constraint need to be additionally enforce. The nice aspect of Casadi is that you can define e.g var_dq to be your decision variable, then call var_q the integration of $dq$ from robot.q0, and use var_q as if it was you variable for the simple robot case. More about that in practice later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The optimization problem can now be written:\n",
    "\n",
    "Decide:\n",
    "- $X = [ x_0 ... x_T ]$ the state trajectory (with $x = (robot.q0 + dq, v)$)\n",
    "- $A = [a_0 ... a_{T-1} ]$ the accelerations\n",
    "\n",
    "Minimizing: $\\sum_{t=0}^{T-1} v_t^2 + e(q_T)^2$\n",
    "\n",
    "Subject to:\n",
    "- $x_0 = (robot.q0, 0)$\n",
    "- $\\forall t=0..T-1$\n",
    "$$ v_{t+1} = v_t + a_t \\Delta t$$\n",
    "$$ q_{t+1} = q_t \\oplus v_{t+1} \\Delta t$$\n",
    "$$ \\forall c  \\textrm{ contacts }, \\alpha_c{t} (q_t,v_t,a_t) = \\kappa$$\n",
    "\n",
    "where $\\alpha$ is the acceleration of the contact interfaces, and $\\kappa$ is a numerical trick we will explain below.\n",
    "\n",
    "Before that, let's write the helpers needed to go around the SX/MX Casadi typing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_acc_error\n",
    "# Sym graph for the operational error\n",
    "error_tool = casadi.Function(\n",
    "    \"etool3\", [cx], [cdata.oMf[endEffector_ID].translation - Mtarget.translation]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_acc_cnext\n",
    "# Sym graph for the integration operation x' = [ q+vDT+aDT**2, v+aDT ]\n",
    "cnext = casadi.Function(\n",
    "    \"next\",\n",
    "    [cx, caq],\n",
    "    [\n",
    "        casadi.vertcat(\n",
    "            cpin.integrate(cmodel, cx[:nq], cx[nq:] * DT + caq * DT**2),\n",
    "            cx[nq:] + caq * DT,\n",
    "        )\n",
    "    ],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Contact model\n",
    "\n",
    "We will consider 2 contact types: punctual 3d contacts, and flat 6d contacts. They are defined by their type, and the frame of the robot where the contact occurs. We will store this information using the SimpleNamespace type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from types import SimpleNamespace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_acc_frames\n",
    "Mtarget = pin.SE3(pin.utils.rotate(\"y\", 3), np.array([-0.1, 0.2, 0.45094]))  # x,y,z\n",
    "contacts = [SimpleNamespace(name=\"left_sole_link\", type=pin.ContactType.CONTACT_6D)]\n",
    "endEffectorFrameName = \"right_sole_link\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_acc_framesId\n",
    "endEffector_ID = model.getFrameId(endEffectorFrameName)\n",
    "for c in contacts:\n",
    "    c.id = model.getFrameId(c.name)\n",
    "    assert c.id < len(model.frames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can display them in the viewer for more understanding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_acc_viz\n",
    "# --- Add box to represent target\n",
    "# Add a vizualization for the target\n",
    "boxID = \"world/box\"\n",
    "viz.addBox(boxID, [0.05, 0.1, 0.2], [1.0, 0.2, 0.2, 0.5])\n",
    "# Add a vizualisation for the tip of the arm.\n",
    "tipID = \"world/blue\"\n",
    "viz.addBox(tipID, [0.08] * 3, [0.2, 0.2, 1.0, 0.5])\n",
    "for c in contacts:\n",
    "    c.viz = f\"world/contact_{c.name}\"\n",
    "    viz.addSphere(c.viz, [0.07], [0.8, 0.8, 0.2, 0.5])\n",
    "\n",
    "\n",
    "def displayScene(q, dt=1e-1):\n",
    "    \"\"\"\n",
    "    Given the robot configuration, display:\n",
    "    - the robot\n",
    "    - a box representing endEffector_ID\n",
    "    - a box representing Mtarget\n",
    "    \"\"\"\n",
    "    pin.framesForwardKinematics(model, data, q)\n",
    "    M = data.oMf[endEffector_ID]\n",
    "    viz.applyConfiguration(boxID, Mtarget)\n",
    "    viz.applyConfiguration(tipID, M)\n",
    "    for c in contacts:\n",
    "        viz.applyConfiguration(c.viz, data.oMf[c.id])\n",
    "    viz.display(q)\n",
    "    time.sleep(dt)\n",
    "\n",
    "\n",
    "def displayTraj(qs, dt=1e-2):\n",
    "    for q in qs[1:]:\n",
    "        displayScene(q, dt=dt)\n",
    "\n",
    "\n",
    "displayScene(robot.q0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Contacts are then mathematically defined by a forward geometry function to remain constant: the world position of a point for 3D contacts, the relative placement of two frames for 6D contacts. \n",
    "\n",
    "As we formulate the problem in accelerations, we will rather differentiate twice this constraint to have it at the acceleration level. This is more generic, and assuming you start from a state with a valid contact configuration, then it is leading to respecting the initial constraint as well.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_acc_helper_contact\n",
    "# Sym graph for the contact constraint and Baugart correction terms\n",
    "# Works for both 3D and 6D contacts.\n",
    "# Uses the contact list <contacts> where each item must have a <name>, an <id> and a <type> field.\n",
    "dpcontacts = {}  # Error in contact position\n",
    "vcontacts = {}  # Error in contact velocity\n",
    "acontacts = {}  # Contact acceleration\n",
    "\n",
    "for c in contacts:\n",
    "    if c.type == pin.ContactType.CONTACT_3D:\n",
    "        p0 = data.oMf[c.id].translation.copy()\n",
    "        dpcontacts[c.name] = casadi.Function(\n",
    "            f\"dpcontact_{c.name}\",\n",
    "            [cx],\n",
    "            [-(cdata.oMf[c.id].inverse().act(casadi.SX(p0)))],\n",
    "        )\n",
    "        vcontacts[c.name] = casadi.Function(\n",
    "            f\"vcontact_{c.name}\",\n",
    "            [cx],\n",
    "            [cpin.getFrameVelocity(cmodel, cdata, c.id, pin.LOCAL).linear],\n",
    "        )\n",
    "        acontacts[c.name] = casadi.Function(\n",
    "            f\"acontact_{c.name}\",\n",
    "            [cx, caq],\n",
    "            [cpin.getFrameClassicalAcceleration(cmodel, cdata, c.id, pin.LOCAL).linear],\n",
    "        )\n",
    "    elif c.type == pin.ContactType.CONTACT_6D:\n",
    "        p0 = data.oMf[c.id]\n",
    "        dpcontacts[c.name] = casadi.Function(f\"dpcontact_{c.name}\", [cx], [np.zeros(6)])\n",
    "        vcontacts[c.name] = casadi.Function(\n",
    "            f\"vcontact_{c.name}\",\n",
    "            [cx],\n",
    "            [cpin.getFrameVelocity(cmodel, cdata, c.id, pin.LOCAL).vector],\n",
    "        )\n",
    "        acontacts[c.name] = casadi.Function(\n",
    "            f\"acontact_{c.name}\",\n",
    "            [cx, caq],\n",
    "            [cpin.getFrameAcceleration(cmodel, cdata, c.id, pin.LOCAL).vector],\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yet in practice, we barely get perfect 0 contact acceleration (because nonexactness of the numerical optimization), and we don't have a perfect numerical integrator. Both lead to some drift in the contact interface. We then have to correct this nonphysical drift with a nonphysical corrector, asking the acceleration to follow the corrector $\\kappa$. The simplest solution is the so-called Baumgart correction term, which is a simple linear feedback in both contact position and velocity.\n",
    "\n",
    "$$ \\kappa(q,v) = -K_p (p(q) \\ominus p_0) - K_v \\nu$$\n",
    "\n",
    "where $p(q)$ is the contact location, that should be constant to $p_0$, and $\\nu is the contact velocity, $K_p,K_v$ are the feedback gains typically chosen so that $K_p = 2\\sqrt{K_v}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_acc_Baumgart\n",
    "Kp = 200\n",
    "Kv = 2 * np.sqrt(Kp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_acc_contact_placement\n",
    "# Get initial contact position (for Baumgart correction)\n",
    "pin.framesForwardKinematics(model, data, robot.q0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_acc_corrector\n",
    "cbaumgart = {\n",
    "    c.name: casadi.Function(\n",
    "        f\"K_{c.name}\", [cx], [Kp * dpcontacts[c.name](cx) + Kv * vcontacts[c.name](cx)]\n",
    "    )\n",
    "    for c in contacts\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Writing the final problem\n",
    "You can now write the full problem, following the steps.\n",
    "\n",
    "1. Declare the $X$ and $A$ variables. Remember $X$ is written in the tangent, first as $var_dxs$ where each step is of size $n_{dx} = 2 n_v$, then as $var_xs$ by integrating (cintegrate) from the state $robot.q0,0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%do_not_load tp2/generated/trajopt_acc_ocp1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Compute the cost by summing the integration terms and the terminal cost."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%do_not_load tp2/generated/trajopt_acc_ocp2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Enforce the initial constraint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%do_not_load tp2/generated/trajopt_acc_ocp3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Enforce the integration constraint "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%do_not_load tp2/generated/trajopt_acc_integration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. Enforce the contact constraint (first without correction, then with correction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%do_not_load tp2/generated/trajopt_acc_ocp4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "6. Solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_acc_ocp5\n",
    "### SOLVE\n",
    "opti.minimize(totalcost)\n",
    "opti.solver(\"ipopt\")  # set numerical backend\n",
    "opti.callback(lambda i: displayScene(opti.debug.value(var_xs[-1][:nq])))\n",
    "\n",
    "# Caution: in case the solver does not converge, we are picking the candidate values\n",
    "# at the last iteration in opti.debug, and they are NO guarantee of what they mean.\n",
    "try:\n",
    "    sol = opti.solve_limited()\n",
    "    sol_xs = [opti.value(var_x) for var_x in var_xs]\n",
    "    sol_as = [opti.value(var_a) for var_a in var_as]\n",
    "except:\n",
    "    print(\"ERROR in convergence, plotting debug info.\")\n",
    "    sol_xs = [opti.debug.value(var_x) for var_x in var_xs]\n",
    "    sol_as = [opti.debug.value(var_a) for var_a in var_as]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally display the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp2/generated/trajopt_acc_ocp6\n",
    "print(\"***** Display the resulting trajectory ...\")\n",
    "displayScene(robot.q0, 1)\n",
    "displayTraj([x[:nq] for x in sol_xs], DT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}