drake_symbolic_and_autodiff_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams['axes.grid'] = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Symbolic Computation and Optimization\n",
    "\n",
    "Core idea: make your program keep track of what computations you have performed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "z:  (x + pow(y, 2))\n",
      "formula:  ((x + pow(y, 2)) = 1)\n"
     ]
    }
   ],
   "source": [
    "from pydrake.all import (\n",
    "    Variable\n",
    ")\n",
    "# Declare some Variables\n",
    "x = Variable(\"x\")\n",
    "y = Variable(\"y\")\n",
    "# Declare an Expression\n",
    "z = x + y**2\n",
    "print \"z: \", z\n",
    "# Declare a Formula using z\n",
    "z_equals_one = (z == 1)\n",
    "print \"formula: \", z_equals_one"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can do simple things, like do evaluations or substitutions..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using x = 3.000000, y = 4.000000...\n",
      "Partial evaluation (for only x): z =  (3 + pow(y, 2))\n",
      "Partial evaluation (for only y): z =  (16 + x)\n",
      "Full evaluation: z =  19.0\n",
      "Substituting x = y^2: z =  (2 * pow(y, 2))\n"
     ]
    }
   ],
   "source": [
    "xi = 3.\n",
    "yi = 4.\n",
    "print \"Using x = %f, y = %f...\" % (xi, yi)\n",
    "# Partial evaluation (specifying *some* variables)\n",
    "print \"Partial evaluation (for only x): z = \", z.EvaluatePartial({x:xi})\n",
    "print \"Partial evaluation (for only y): z = \", z.EvaluatePartial({y:yi})\n",
    "# Full evaluation (specifying *all* variables)\n",
    "print \"Full evaluation: z = \", z.Evaluate({x:xi, y:yi})\n",
    "print \"Substituting x = y^2: z = \", z.Substitute({x:y**2})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And you can ask for derivatives..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Derivative of z w.r.t. x:  1\n",
      "Derivative of z w.r.t. y:  (2 * y)\n",
      "2nd derivative of z w.r.t. y:  2\n",
      "You can ask for Jacobians if you want:  [<Expression \"1\"> <Expression \"(2 * y)\">]\n"
     ]
    }
   ],
   "source": [
    "print \"Derivative of z w.r.t. x: \", z.Differentiate(x)\n",
    "print \"Derivative of z w.r.t. y: \", z.Differentiate(y)\n",
    "print \"2nd derivative of z w.r.t. y: \", z.Differentiate(y).Differentiate(y)\n",
    "print \"You can ask for Jacobians if you want: \", z.Jacobian([x, y])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mathematical Programming\n",
    "\n",
    "As you might guess, you can build this up to cool stuff! The primary thing we use this for is setting up and solving optimizations. By specifying an optimization (i.e. an objective + a list of constraints) symbolically, Drake figures out what class of optimization it is and dispatches to an appropriate solver.\n",
    "\n",
    "http://drake.mit.edu/doxygen_cxx/group__solvers.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result:  SolutionResult.kSolutionFound\n",
      "Solver used:  OSQP\n",
      "Optimizing x:  0.0\n"
     ]
    }
   ],
   "source": [
    "from pydrake.all import (MathematicalProgram)\n",
    "prog = MathematicalProgram()\n",
    "# NewContinuousVariables spawns a 1x1 array of decision\n",
    "# variables. We're more concerned with individual variables\n",
    "# for this example, so take out the first element of that\n",
    "# array to get the actual Variable.\n",
    "x = prog.NewContinuousVariables(1, \"x\")[0]\n",
    "prog.AddQuadraticCost(x**2)\n",
    "prog.AddBoundingBoxConstraint(0., 1., x)\n",
    "result = prog.Solve()\n",
    "print \"Result: \", result\n",
    "print \"Solver used: \", prog.GetSolverId().name()\n",
    "print \"Optimizing x: \", prog.GetSolution(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result:  SolutionResult.kSolutionFound\n",
      "Solver used:  SNOPT\n",
      "Optimizing x:  0.0\n"
     ]
    }
   ],
   "source": [
    "prog = MathematicalProgram()\n",
    "x = prog.NewContinuousVariables(1, \"x\")[0]\n",
    "prog.AddCost(x**4)\n",
    "prog.AddBoundingBoxConstraint(0., 1., x)\n",
    "result = prog.Solve()\n",
    "print \"Result: \", result\n",
    "print \"Solver used: \", prog.GetSolverId().name()\n",
    "print \"Optimizing x: \", prog.GetSolution(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SOS Stability Analysis Example\n",
    "\n",
    "As you might guess, you can build this up to cool stuff! For example, after defining some dynamics, you can pass symbolic variables through it, and then do stability analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xd8FHX+x/HXJ50QINSQQq8iPYBSRCKigAULKJ6iYMGGeud5lvPO84pnO3uvgDUqiiCCSAmoIEJQOgKhSQhFOkEIBD6/P3bwt8Yk7G6yO5vweT4e88jOzHdm3swu+9npoqoYY4wxvopwO4AxxpiKxQqHMcYYv1jhMMYY4xcrHMYYY/xihcMYY4xfrHAYY4zxixUO4zoReVBE3gmDHHVFZJWIxLmdJVRE5BwR+TQEy3lSRG7y6o8VkR9FpF6wl23KnxUOU6GIyAYROTtIs78XGK2qh4I0/3D0X+CRss5ERNqISLaI7Ha66SLSxqvJ48D9IhIDoKoFwJvAPWVdtgk9KxzG4PkFDFwDuL7lEyoi0hWooarzymF2ecBgoBZQB5gIZB4fqapbgB+BC72meQ+4xln3pgKxwmHKnYg0E5FdItLZ6U8RkR0i0sfpbyIis0Vkv4hMw/NF4z39hSKyXET2iMgsETnFGf420BD4TETyReTucox9GrBHVXMDmVhEYkRkkYjc5vRHisgcEXkgwPkNd6Z/TkT2Ort1+nqNTxGRic56zhGRG5zhcSJyUETqOP1/E5FCEanu9P9HRJ52ZjMAmO01zx7O+9TA6e/gvAetT5RXVfeo6gb13IpCgKNA8yLNZgHneU2TC+wGTvd3/Rh3WeEw5U5V1+LZBfGuiMQDo4ExqjrLafIesBBPwfg3nl/6AIhIS+B94I9AXWAynkIRo6rDgJ+AC1Q1QVUfK7psEWnofNmV1P2hhNjtgFVl+DcfBq4C/uUUunuBSOChQOeJp5itw7Oe/gF8IiK1nHHvA7lACp5f+v8Vkb7ObrYFwJlOu97ARqCnV//xYvGbf7OqzgVeAcaKSBXgbeBvqvqjr4FFZA9wCHgOz24wbyuBDj4MM2HOCocJClV9DVgDfAckA/eD54sd6Ar8XVULVPUr4DOvSS8HPlfVaap6BPgfUAXo4eNyf1LVxFK690qYNBHYH8i/1WvZy4D/AOOBu4Bhqnq0DLPcDjytqkdU9QM8X/LnOVsEvYB7VPWQqi4CXgeGOdPNBs4UkSigPfCs0x+HZ91/7bQr7t/8IFADmI9n99ML/gRW1URn+lHAD0VG73eWeaJhJsxZ4TDB9BrQFnjOORgKnl/Iu1X1gFe7jV6vU7z7VfUYsAlIDXLW3UC1kkaKyJXO7rF8EZlSynzGAo2Byaq6poR5NfSaV34p89qsv70L6UY86ycF2KWq+4uMO76OZgN9gM7AUmAani2Q04EcVd3htPvdv9kp1mPwvG9PFFm+T5z39mXgrSJnTVUD9hRpXtwwE+ascJigEJEE4GngDeBBr10sW4CaIlLVq3lDr9d5QCOv+QjQANjsDCr1i6zol3Ix3ZUlTLoEaFnSfFX1XWf3WIKqDiglwovAJOBcEelVwrx+8ppXQinzSnX+/cc1xLN+8oBaIlKtyLjj62gu0Aq4GJitqiuc8efhdUyDYv7NIpKKZ7fYaOCJMhy4jgDi+W3BPwVYXKRdccNMuFNV66wr9w5PwfjQef3q8ddO/zw8u6Bi8Oxy2Qe844xrBRwA+gLReHb5rANivKYdGYS8McDPQGoZ5jEMWAskAH84/jrAeQ0HCoE7nPUwxFlPtZ3xXwPPA3F4dkdtA/p5TT/XaX+G0/+R0z/Eq01nYLVXv+DZOnnUeT0VeMxr/Bg8x6qKy9sP6ITnuE51PLvH8oA4rzZfApd59acCO4FYtz+v1vn5+XQ7gHWVrwMG4fn1W8vpTwBygCud/qbOF1++80X1/PHC4Yy/GFgB7MXzC/nUIvP+Cc/ujbvKOffjeI4bBDJtQ+dLsKfXsA+A1wKc33BgjrNu9gKrgXO8xqfh2bLZ5RSom4pM/zBw8PiXMp5jDgokFWm3ADjNeX0Hnq2Q40U6BU8xPV58ZgA3lJB3CJ7TbfOdaSYD7b3GJ+M5mB/jNewvwJNuf16t878T5w005qQnInXxFLROqnrQ5SzDgetVtdjdXeW4nHOAW1T1ohO0i8GzS6m9eo6D+LucJ4C1qvqi0x/rzK+3qm73P7lxkxUOY8JQqAqHMYGwg+PGGGP8Ylscxhhj/GJbHMYYY/wS5XaAYKhTp442btw4oGkPHDhA1apVT9wwxCyXfyyXfyyXfypjroULF+5Q1bo+NXb7tK5gdOnp6RqorKysgKcNJsvlH8vlH8vln8qYC8hWH79jbVeVMcYYv1jhMMYY4xcrHMYYY/xihcMYY4xfrHAYY4zxi2uFQ0QaiEiWiKx0HhN6RzFtRESedR6NueT4o0iNMca4x83rOAqBP6vq985zBRaKyDT1PDvguAFAC6c7DXjJ+WuMMcYlrm1xqOoWVf3eeb0fz7OHiz7lbRDwlnOa8TwgUUSSg5Hn2DHlhawcNuwty5M+jTGm8guLe1WJSGPgK6Ctqu7zGj4JeERVv3H6Z+B5XkJ2MfMYCYwESEpKSs/MzPQrw4Ejyt/nHET1GP/qWZVqMXLiiUIoPz+fhITSHhbnDsvlH8vlH8vln7LkysjIWKiqXXxq7OuVgsHq8DzkZyFwSTHjPgd6efXPANJPNM9ArxxfvGm3Nrtvkv7htW/1SOHRgOYRLJXxStVgslz+sVz+qYy5qChXjotINPAx8K6qflJMk1w8z5s+Lg3P4yiDon1aIle3iWFOzk4e/3JVsBZjjDEVmptnVQme51KvVNUnS2g2EbjaObvqdGCvqm4JZq7eadFceVpDXpm9js+XBHVRxhhTIbl5VlVPYBiwVEQWOcP+iufZzajqy3ieWzwQz/OqfwFGhCLYAxe0YcWWffxl3GJaJCXQMqlaKBZrjDEVgmuFQz0HvEs9Au3sd7s1NIn+X2xUJC9flc55z37DjW8vZMKonlSPiw51DGOMCUt25XgJkqrH8eKVndm06xfu/GARx465f/aZMcaEAyscpejWpBZ/P78N01du5392sNwYY4BK+gTA8nR190as2rafF2etpXm9BC7pnOZ2JGOMcZVtcZyAiPDPC0+le9Pa3PvxUhZu3OV2JGOMcZUVDh9ER0bw4pWdSUmMY+RbC8nd/YvbkYwxxjVWOHxUs2oMr1/TlcNHj3H92GzyCwrdjmSMMa6wwuGH5vUSePHKzqzZns8fMxdx1M60MsachKxw+OmMFnX5xwVtmL5yGw9PXul2HGOMCTk7qyoAV3dvzLqfD/D6N+tJSazCtb2auB3JGGNCxgpHgP5+fhu27j3Evz9fQf0acQxsF5THhBhjTNixXVUBiowQnh7akc4Na/LHDxYxf72dpmuMOTlY4SiDuOhIXr+6C2mJVbjhrWxytu93O5IxxgSdFY4yqlk1hrHXdiM6MoJr3lzA9n2H3I5kjDFBZYWjHDSoFc/o4V3Z/cthho9ewL5DR9yOZIwxQWOFo5y0S6vBi1d2ZvW2/Vw/JpuDh4+6HckYY4LCCkc56tOqHk9d3pEFG3dx87sLOVx4zO1IxhhT7tx+5vibIrJdRJaVML6PiOwVkUVO90CoM/rrgg4p/Pfidsxa9TN3fmhXlxtjKh+3r+MYAzwPvFVKm69V9fzQxCkfV3RryL6DR3h4yo9Ui4vmvxe3xfOIdWOMqfhcLRyq+pWINHYzQ7DceGYz9h48wouz1lK9ShT3DTjF7UjGGFMuxPNYbxcDeArHJFVtW8y4PsDHQC6QB9ylqstLmM9IYCRAUlJSemZmZkB58vPzSUhICGjaolSVt1ccZuamQi5pEc2FzWICnld55ipPlss/lss/lss/ZcmVkZGxUFW7+NRYVV3tgMbAshLGVQcSnNcDgTW+zDM9PV0DlZWVFfC0xTl69Jj+6YMftNE9k/T5mWsCnk955yovlss/lss/lss/ZckFZKuP39tuH+Molaru83o9WUReFJE6qrrDzVz+iIgQHh/cgWPHlMenrkIEbunT3O1YxhgTsLAuHCJSH9imqioi3fCcBbbT5Vh+i4wQnrisIwo89sUqIkW48cxmbscyxpiAuFo4ROR9oA9QR0RygX8A0QCq+jIwGLhZRAqBg8BQZ5OqwomMEJ4Y0oFjCg9P+ZEIEW7o3dTtWMYY4ze3z6q64gTjn8dzum6lEBUZwVOXdeCYKg9NXokIXH+GFQ9jTMUS1ruqKqOoyAievrwjqsp/Pl9JQeExbs2wYx7GmIrDCocLoiMjeHZoJ2IiF/P41FXkFxRy97mt7CJBY0yFYIXDJVGRETx5WUfiY6N4adZafiko5B8XnEpEhBUPY0x4s8LhoogI4aGL2lI1JpLXvl7PgcNHefTS9kRa8TDGhDErHC4TEf468BSqxkbx9PQ1HDx8lKcu70hMlN242BgTnqxwhAER4Y9nt6RqTBQPTV7J/oJCXrqyM1Vj7e0xxoQf+1kbRm7o3ZRHL23HnJwdXPHaPHbkF7gdyRhjfscKR5i5vGtDXh2Wzupt+7n0pbls3HnA7UjGGPMbVjjCUN9TknjvhtPZd/AIl740l6W5e92OZIwxv7LCEaY6N6zJuJt7EBsVyeWvfsuyHYVuRzLGGMAKR1hrVjeB8bf0oFHtqjy5sIB3v9vodiRjjLHCEe7qVY/jwxtPp22dSO4fv4x/fbbCnmNujHGVne9ZAVSLi+aPnWP5Jj+JN+esZ8POAzx7RScS7HRdY4wLbIujgogQ4YEL2vCfi9oye/XPDH5pLrm7f3E7ljHmJGSFo4K56vRGjB3Rjc17DnLRC3PI3rDL7UjGmJOMFY4KqFeLOoy/pScJsVFc8do83v52AxX0+VbGmArI1cIhIm+KyHYRWVbCeBGRZ0UkR0SWiEjnUGcMV83rJTBhVC96Na/D3ycs5y/jlnDoyFG3YxljTgJub3GMAfqXMn4A0MLpRgIvhSBThVGjSjRvXNOV2/u2YNzCXIa8/C2b9xx0O5YxppJztXCo6ldAaTvpBwFvqcc8IFFEkkOTrmKIiBDu7NeS167uwoYdB7jguW+Ym7PD7VjGmEpM3N43LiKNgUmq2raYcZOAR1T1G6d/BnCPqmYX03Yknq0SkpKS0jMzMwPKk5+fT0JCQkDTBpMvubbkH+O5Hw6x5YBySYtozmsaTUSQnypYkdeXGyyXfyyXf8qSKyMjY6GqdvGpsaq62gGNgWUljPsc6OXVPwNIP9E809PTNVBZWVkBTxtMvubaf+iI3vbe99ronkl61evzdPu+Q2GRK9Qsl38sl38qYy4gW3383nb7GMeJ5AINvPrTgDyXslQICbFRPDO0I49c0o7563cx4JmvmWO7rowx5SjcC8dE4Grn7KrTgb2qusXtUOFORBjarSETR/UiMT6aq974jie/XEXh0WNuRzPGVAJun477PvAt0EpEckXkOhG5SURucppMBtYBOcBrwC0uRa2QWtWvxsRRPRncOY1nZ+bwh9e+s6vNjTFl5urNjlT1ihOMV+DWEMWplOJjonh8SAe6N6vNAxOWM+Dpr3nwwlO5pHMqEuQD58aYyincd1WZcnJJ5zSm3HEGpyRX588fLeaWd79n14HDbscyxlRAVjhOIg1qxfP+yNO5d0Brpq/cxrlPf0XWqu1uxzLGVDBWOE4ykRHCTWc2Y8KtvagVH8OI0Qv46/il7D90xO1oxpgKwgrHSapNSnUmjOrJyN5NyZz/E+c89RVZP9rWhzHmxKxwnMTioiP568BT+PjmHiTERjFizAL+mPmDHfswxpTKCoehU8OaTLq9F3f0bcGkJVs4+8nZTFycZ7dqN8YUywqHASA2KpI/9WvJpNt7kVazCre//wPXj81m0y677sMY81tWOMxvtK5fnU9u7sH9A09h7tqd9HtqNs/PXENBoT3rwxjjYYXD/E5UZAQ39G7K9D+fSUarevzvy9UMePprvllj97wyxljhMKVITazCS1elM2ZEV46qctUb33Hre9+zde8ht6MZY1xkhcOcUJ9W9Zj6x9786eyWTFuxjbOemMVzM9bYo2qNOUlZ4TA+iYuO5I6zWzD9T2fSu0Vdnpi2mrP+N4u5eYUcO2ZnXxlzMrHCYfzSsHY8Lw9LJ3Pk6dRKiOHVJQVc/NJcFm4s7QnAxpjKxAqHCcjpTWsz8dZeXN8uhq17D3LpS99y63vfs3HnAbejGWOCzAqHCVhEhNArNZqsu/pwR98WzFi5jb5PzOb+8UvZts8OoBtTWbn9IKf+IrJKRHJE5N5ixg8XkZ9FZJHTXe9GTlO6+Jgo/tSvJV/9JYMrujXkgwWb6P1YFg9PXsluu32JMZWOa4VDRCKBF4ABQBvgChFpU0zTD1S1o9O9HtKQxi/1qsfx74vaMvPPfTivXTKvfr2O3o9l8eyMNeQXFLodzxhTTtzc4ugG5KjqOlU9DGQCg1zMY8pJw9rxPHl5R764ozfdm9XmyWmr6f1YFi9k5djt242pBMStG9mJyGCgv6pe7/QPA05T1VFebYYDDwM/A6uBP6nqphLmNxIYCZCUlJSemZkZUK78/HwSEhICmjaYKnKutXuO8mnOEZbuOErVaOjXKJp+jaKpGh28R9dW5PXlBsvln8qYKyMjY6GqdvGpsaq60gFDgNe9+ocBzxVpUxuIdV7fBMz0Zd7p6ekaqKysrICnDabKkGvRT7v1ujELtNE9k7TtA1/oY1+s1J35Ba7nCiXL5R/L5Z+y5AKy1cfvbzd3VeUCDbz604A87waqulNVC5ze14D0EGUzQdChQSKvX9OFybefwRkt6/DirLX0enQm/5m0grw9B92OZ4zxUZSLy14AtBCRJsBmYCjwB+8GIpKsqluc3guBlaGNaIKhTUp1XrwynTXb9vN8Vg6j525gzNwNnN8+mRt6N+XUlBpuRzTGlMK1wqGqhSIyCpgKRAJvqupyEfkXnk2micDtInIhUAjsAoa7ldeUvxZJ1XhmaCf+cm4rRs/ZQOb8n/h0UR69mtfhht5N6d2iDiLBOw5ijAmMm1scqOpkYHKRYQ94vb4PuC/UuUxopdWM5+/nt+H2vi1477ufGD1nPde8OZ/W9atxXa8mXNAhhbjoSLdjGmMcduW4CRs1qkRzc59mfHPPWfxvSAdU4S/jltDjkZk89sWPbLbjIMaEBVe3OIwpTkxUBIPT07i0cypzcnby1rcbeHn2Wl6evZZ+bZK4pntjujerbbuxjHGJFQ4TtkSEXi3q0KtFHTbvOci78zaSuWATU5dvo3m9BK7u3ohLOqeREGsfY2NCyXZVmQohNbEKd/dvzdx7z+KJIR2oGhPJAxOW0+2h6dwzbgnf/7T7+LU/xpggs59qpkKJi47k0vQ0Lk1PY9GmPWTO/4mJi/P4IHsTLZMSuLxrQy7plOp2TGMqNSscpsLq2CCRjg0S+dv5bZi0OI/MBZv496QVPDrlRzrVFaLTdtC9aW0iIuxYiDHlyQqHqfASYqMY2q0hQ7s1ZOWWfXywYBMfLdjAla9/R4NaVbi4UxqXdEqlcZ2qbkc1plKwwmEqlVOSq/PghafSo+p2DtZuybiFuTw3cw3PzlhD54aJXNw5jQvaJ5MYH+N2VGMqrBMWDhGJ9bpfVInDjAknMZHCOR1TGdQxla17DzFh0WY+/j6Xv3+6jH9/toKzWtfj4s6pZLSqR0yUnSNijD982eL4FujswzBjwlL9GnHceGYzRvZuyvK8fYz/YTMTFm3mi+VbqRkfzfntU7ioUwqdGtS04yHG+KDEwiEi9YFUoIqIdAKO/4+qDsSHIJsx5UpEaJtag7apNbhvQGu+XrODT37YzIfZm3h73kZSE6twXvtkLmifQtvU6naBoTElKG2L41w8NxVMA570Gr4f+GsQMxkTdFGREWS0rkdG63rsP3SE6Su38dniLbz5zXpe/WodjWvHc0GHFC7okELLpGpuxzUmrJRYOFR1LDBWRC5V1Y9DmMmYkKoWF83FndK4uFMae345zBfLtvLZkjxeyMrhuZk5tEqqxvntkzm/QwpN7MwsY3w6xjFDRJ4Eejv9s4F/qere4MUyxh2J8TG/ntr78/4CpizbwmeL83hi2mqemLaadqk1OK99MgPbJtOwtu2xNScnXwrHG8Ay4DKnfxgwGrgkWKGMCQd1q8VydffGXN29MXl7DvL5ki18tiSPR6b8yCNTfuTUlOoMbJfMgLb1aVo3/J4/bUyw+FI4mqnqpV79/xSRRcEKZEw4Skmswg29m3JD76Zs2vULXyzbyuRlW3h86ioen7qK1vWrMaBtMgPb1aeFHRMxlZwvheOgiPRS1W8ARKQnUC4PRhCR/sAzeJ4A+LqqPlJkfCzwFp5nje8ELlfVDeWxbGMC1aBW/K9FJG/PQb5YtpUpy7bw9IzVPDV9Nc3qVnW2RJI5JbmanZ1lKh1fCsfNeA6S18BzSm65PMJVRCKBF4B+QC6wQEQmquoKr2bXAbtVtbmIDAUeBS4v67KNKS8piVW4tlcTru3VhO37DjF1+VYmL93664H1xrXjGdDOc0ykbWp1t+MaUy5OWDhUdRHQQUSqO/37ymnZ3YAcVV0HICKZwCDAu3AMAh50Xo8DnhcRUbt/tglD9arHMax7Y4Z1b8yO/AK+XL6NKcu28OpX63hp1lrSalahbY0jVG+6m45piXaxoamwpKTvYBG5s7QJVfXJ0safcMEig4H+qnq90z8MOE1VR3m1Wea0yXX61zptdhQzv5HASICkpKT0zMzMgHLl5+eTkBB+Bzotl3/CKVf+YeWH7YUs2HaU5TsKOapCrTghPSmSbvWjaJYYQYTLu7PCaX15s1z+KUuujIyMharaxZe2pW1xHD/C1wroCkx0+i8Avgoo2W8V9z+laBXzpY1noOqrwKsAXbp00T59+gQUatasWQQ6bTBZLv+EW67znb+fT8uioHYLJi/dyuzVPzNt4yHqV4+jf9v6nNc+mfSG7tz2JNzW13GWyz+hylXaBYD/BBCRL4HOqrrf6X8Q+Kgclp0LNPDqTwPySmiTKyJRQA08x1iMqZCqRgvndU7jks5p7D90hBkrtzN56Rbem/8TY+ZuoF61WAa0rc/Adsl0aVyLSNudZcKQLwfHGwKHvfoPA43LYdkLgBYi0gTYDAwF/lCkzUTgGjw3VRwMzLTjG6ayqBYXzUWdUrmoUyr5BYXMWLmNyUu3kLlgE2O/3UhdryLS1YqICSO+FI63gfkiMh7PbqKLgbFlXbCqForIKGAqntNx31TV5SLyLyBbVSfiufjwbRHJwbOlMbSsyzUmHCXERjHIuQ38gYJCZv7o2RL5MHsTb327kToJsfRvm8TAdsmc1qS2FRHjKl/OqnpIRKYAZziDRqjqD+WxcFWdDEwuMuwBr9eHgCHlsSxjKoqqsVG/3mDxQEEhWau2M2XpVj5euJl35v1EnYQYzj21Pue1S6Zbk1pERdrzRExo+fQEQFX9Hvg+yFmMMUVUjY3i/PYpnN8+hV8OFzJr1c98vnQLn3y/mXe/+4naVWM4v30yF3ZMpXPDRLvY0ISEPTrWmAoiPiaKge2SGdgumYOHjzJr1XYmLfn/YyINalVhUIdUBnVMsduemKCywmFMBVQlJpIB7ZIZ0C6Z/YeO8OXybXy6aDMvzsrh+awcTkmuzkUdPbu7UhKruB3XVDK+PHO8TZHbgCAifVR1VtBSGWN8Vi0umkvT07g0PY2f9xcwaUkeExbl8fCUH3l4yo90a1KLQR1TOK9dMonxMW7HNZWAL1scH4rI28BjQJzztwvQPZjBjDH+q1stlhE9mzCiZxM27jzAxEV5fLpoM/ePX8Y/J66gX5skBndJo3eLunZmlgmYL4XjNDw3F5yL52ryd4GewQxljCm7RrWrclvfFow6qznL8/YxbmEuExZt5vOlW0iqHsvFndIY0iWNZvYsEeMnXwrHETy3Ua+CZ4tjvaoeC2oqY0y5ERHaptagbWoN7hvYmpkrt/PRwlxe+3odL89eS+eGibSvdoT0Q0eoFhftdlxTAfhSOBYAE/Dcr6o28IqIDFbVwUFNZowpd7FR/39Qffu+Q4z/YTMfLcxlzE+HyVw9nYFtk7ny9IZ0bljTTu01JfKlcFynqtnO663AIOdOtsaYCqxe9ThuPLMZI3s35c0JM1mr9Zi4KI9PfthM6/rVuPL0RlzUMcW2QszvnPCSU6+i4T3s7eDEMcaEmojQLDGS/17cju/+2peHL2lHZITw90+Xcdp/Z3DfJ0tZtnmv2zFNGLHrOIwxv6oaG8UV3RoytGsDFufu5d15Gxn/Qy7vz/+JDg0Sueq0hlzQIYW46Ei3oxoX2U1ujDG/IyJ0bJDI40M68N19Z/OPC9pwoKCQv4xbQq9HZ/LUtNX8vL/A7ZjGJbbFYYwpVY34aEb0bMLwHo2Zu3Ynb3yznmdmrOGlWWu5sGMK1/ZsQpsUe576ycQKhzHGJyJCz+Z16Nm8Dut+zmf0nA2MW5jLuIW5dG9am+t6NeGs1vXsWeonAdtVZYzxW9O6Cfz7orbMu68v9w1ozcadB7j+rWzOfmo2H2Zv4nChXepVmVnhMMYErEZ8NDee2YzZd2fw3BWdiIuK5O5xSzjz8Sze/GY9vxwudDuiCQJXCoeI1BKRaSKyxvlbs4R2R0VkkdNNDHVOY4xvoiMjuKBDCp/f3osxI7rSoFY8/5q0gp6PzOTZGWvY+8sRtyOacuTWFse9wAxVbQHMcPqLc1BVOzrdhaGLZ4wJhIjQp1U9PryxO+Nu6k6nhjV5ctpqejwyg0em/MjuA4fdjmjKgVuFYxD//9zyscBFLuUwxgRJl8a1eHN4V6bccQZnnZLEK1+tpdejM3niy1W2BVLBiaqGfqEie1Q10at/t6r+bneViBQCi4BC4BFV/bSUeY4ERgIkJSWlZ2ZmBpQtPz+fhITwu1uo5fKP5fJPKHJt3n+MT9ceZsHWo1SJgv6No+nXKJr46JLPwjqZ11cgypIrIyNjoap28amxqgalA6YDy4rpBgF7irTdXcI8Upy/TYENQDNflp2enq6BysrKCnjaYLJc/rFc/gllruWb9+oNYxeN4O78AAATW0lEQVRoo3smafsHp+rzM9fogYIjrufyR2XMBWSrj9/vQbuOQ1XPLmmciGwTkWRV3SIiycD2EuaR5/xdJyKzgE7A2mDkNcaERpuU6rx6dReW5u7lqemreXzqKsbM3cCd/VoyJD2NqEg72TPcufUOTQSucV5fg+e27b8hIjVFJNZ5XQfPw6NWFG1njKmY2qXV4M3hXfn45u40rBXPfZ8sZcAzXzNj5bbjexxMmHKrcDwC9BORNUA/px8R6SIirzttTgGyRWQxkIXnGIcVDmMqmfRGtRh3U3deviqdwmPKdWOzGfrqPBZv2uN2NFMCV245oqo7gb7FDM8GrndezwXahTiaMcYFIkL/tvXpe0o9Muf/xNPT1zDohTmcnhxJq04HSa5Rxe2IxovtTDTGhI3oyAiGdW/M7LszuO2s5mRvO8pZ/5vNC1k5FBQedTuecVjhMMaEnYTYKP58Tise7lWF3i3r8PjUVZzz1Fd2/CNMWOEwxoStuvERvDKsC29f142oCOG6sdmMGLOAdT/nux3tpGaFwxgT9s5oUZcv/tibv513CtkbdtP/6a95Zvoa233lEiscxpgKIToyguvPaMrMu87k3Lb1eWr6as579hvmr9/ldrSTjhUOY0yFUq9aHM9d0YnRI7py8PBRLnvlW+77ZInd/yqErHAYYyqkjFb1mHZnb244owkfLNhE3ydn89niPDt4HgJWOIwxFVZ8TBT3n9eGiaN6kVwjjtve/4Fb3v2enfkFbker1KxwGGMqvLapNRh/Sw/u7t+KGSu3c85TXzFl6Ra3Y1VaVjiMMZVCVGQEt/Rpzme39SI5MY6b3/2e29//wR4eFQRWOIwxlUqr+tUYf0tP7uzXkslLt9DPuXDQlB8rHMaYSic6MoLb+7Zgwqie1EmI4bqx2TwwYRmHjth1H+XBCocxptI6NaUGE0b15LpeTXjr240Men4Oq7budztWhWeFwxhTqcVGRfL389swZkRXdh4o4ILnv+GtbzfYabtlYIXDGHNS6NOqHlPu6E2PZrV5YMJyrh+bzS47cB4QVwqHiAwRkeUickxESnw4uoj0F5FVIpIjIveGMqMxpvKpWy2W0cO78sD5bfh6zQ7Oe/Zrvv9pt9uxKhy3tjiWAZcAX5XUQEQigReAAUAb4AoRaROaeMaYykpEuLZXEz65pQdRkcLlr3zLmDnrbdeVH1wpHKq6UlVXnaBZNyBHVdep6mEgExgU/HTGmJNB29QaTBp1Bme2rMuDn63gtvd/IL+g0O1YFYK4WWVFZBZwl/PI2KLjBgP9VfV6p38YcJqqjiphXiOBkQBJSUnpmZmZAWXKz88nISEhoGmDyXL5x3L552TOdUyVKeuPMG71EZKqCrd1jCO1Wum/qSvj+srIyFioqiUeOvgNVQ1KB0zHs0uqaDfIq80soEsJ0w8BXvfqHwY858uy09PTNVBZWVkBTxtMlss/lss/lkt1bs4OTf/3NG39tyk6aXFeqW0r4/oCstXH7/eogEqTbwXp7DLOIhdo4NWfBuSVcZ7GGFOs7s1qM/n2Xtz87vfc+t73/Li1OX86uyUREeJ2tLATzqfjLgBaiEgTEYkBhgITXc5kjKnE6lWP470bTuPyLg14bmYON76z0I57FMOt03EvFpFcoDvwuYhMdYaniMhkAFUtBEYBU4GVwIequtyNvMaYk0dsVCSPXNqOBy9ow8wft3PJi3P4aecvbscKK26dVTVeVdNUNVZVk1T1XGd4nqoO9Go3WVVbqmozVX3IjazGmJOPiDC8ZxPGjujGtn0FXPjCN8xdu8PtWGEjnHdVGWOMq3q1qMOEW3tSNyGWq9+Yz7iFuW5HCgtWOIwxphSN61Tl41t6cHrT2tz10WKemrb6pL9Y0AqHMcacQPW4aN4c3pXB6Wk8M2MNry89zOHCY27Hco0VDmOM8UFMVASPD27Pnf1aMievkOGj57P34BG3Y7nCCocxxvhIRLi9bwtuaBfDgg27GPzSXPL2HHQ7VshZ4TDGGD/1TI1m7LXd2Lr3EINfmkvO9ny3I4WUFQ5jjAlAj2Z1yLzxdA4fPcZlr3zLktw9bkcKGSscxhgToFNTajDuph7Ex0RyxavzmJtzclzrYYXDGGPKoHGdqnx8cw9Sa1Zh+OgFfLFsq9uRgs4KhzHGlFFS9Tg+vLE7p6ZW55Z3F/Jh9ia3IwWVFQ5jjCkHifExvHv9afRsXoe7xy3h3e82uh0paKxwGGNMOYmPieK1q7twVut63D9+GWPnbnA7UlBY4TDGmHIUFx3Jy1el069NEv+YuJzXv17ndqRyZ4XDGGPKWUxUBC9e2ZmB7erzn89X8tKstW5HKldBewKgMcaczKIjI3h2aCciIxbz6Bc/Unj0GLf1beF2rHJhhcMYY4IkKjKCpy7rQFSE8MS01URECLdmNHc7Vpm5UjhEZAjwIHAK0E1Vs0totwHYDxwFClW1S6gyGmNMeYiKjOB/Qzqgqjw+dRVx0ZFc16uJ27HKxK0tjmXAJcArPrTNUNWT43JMY0ylFBkh/G9IBwoKj/HvSSuIjYrgqtMbuR0rYK4UDlVdCZ47TRpjzMkgKjKCZ4Z2ouCdhfzt02XERUcyOD3N7VgBETefZCUis4C7StlVtR7YDSjwiqq+Wsq8RgIjAZKSktIzMzMDypSfn09CQkJA0waT5fKP5fKP5fJPWXIdPqo88/0hVuw8xk0dYjktufx+v5clV0ZGxkKfDweoalA6YDqeXVJFu0FebWYBXUqZR4rztx6wGOjty7LT09M1UFlZWQFPG0yWyz+Wyz+Wyz9lzXWg4IgOeWmuNr3vc526bEv5hNKy5QKy1cfv96Bdx6GqZ6tq22K6CX7MI8/5ux0YD3QLVl5jjAmV+Jgo3hjehbapNRj1/g/MW7fT7Uh+CdsLAEWkqohUO/4aOAfPFosxxlR41eKiGTO8Kw1qVuGGsdmsyNvndiSfuVI4RORiEckFugOfi8hUZ3iKiEx2miUB34jIYmA+8LmqfuFGXmOMCYaaVWN467rTSIiL4prR89m06xe3I/nElcKhquNVNU1VY1U1SVXPdYbnqepA5/U6Ve3gdKeq6kNuZDXGmGBKTazCW9d243DhMYa98R078gvcjnRCYburyhhjThYtkqrx5vCubN13iBGjF5BfUOh2pFJZ4TDGmDCQ3qgmL17ZmRVb9nHj29kcLjzmdqQSWeEwxpgwcVbrJB69tD1zcnZy7ydLjl+WEHbsJofGGBNGBqensXn3QZ6avprGtatyexjeUdcKhzHGhJnb+zZn484DPDltNQ1rxXNRp1S3I/2GFQ5jjAkzIsLDl7Zj856D3D1uCSmJVejWpJbbsX5lxziMMSYMxUZF8sqwdNJqVmHk29ms33HA7Ui/ssJhjDFhKjE+htEjuhIhwojR89l94LDbkQArHMYYE9Ya1a7Kq8PSydt7iBvfWciRo+6fpmuFwxhjwlyXxrV47NL2zF+/i39+ttztOHZw3BhjKoKLOqWycss+XvlqHackV+fK09x7gqBtcRhjTAVxd//WnNmyLv+YsJz563e5lsMKhzHGVBCREcKzV3SiYa14bn5nIZv3HHQlhxUOY4ypQGpUiebVq7twuPAYI9/K5uDhoyHPYIXDGGMqmOb1Enj2ik6s2LKPuz8O/T2trHAYY0wFlNG6Hn85txWfLc7jzTkbQrpst54A+LiI/CgiS0RkvIgkltCuv4isEpEcEbk31DmNMSac3XxmM849NYn/Tl4Z0oPlbm1xTAPaqmp7YDVwX9EGIhIJvAAMANoAV4hIm5CmNMaYMCYiPD6kA41qxXPre9+z51BoLg5069GxX6rq8UdczQPSimnWDchxHiF7GMgEBoUqozHGVATV46J5eVg6+YcKeWFRQUiuLBe3HxQiIp8BH6jqO0WGDwb6q+r1Tv8w4DRVHVXCfEYCIwGSkpLSMzMzA8qTn59PQkJCQNMGk+Xyj+Xyj+XyTzjmmrelkCVbDzGiQ1WiI8Tv6TMyMhaqahefGqtqUDpgOrCsmG6QV5v7gfE4BazI9EOA1736hwHP+bLs9PR0DVRWVlbA0waT5fKP5fKP5fJPZcwFZKuP3+9Bu+WIqp5d2ngRuQY4H+jrhC4qF2jg1Z8G5JVfQmOMMYFw66yq/sA9wIWq+ksJzRYALUSkiYjEAEOBiaHKaIwxpnhunVX1PFANmCYii0TkZQARSRGRyQDqOXg+CpgKrAQ+VFX3bwtpjDEnOVfujquqzUsYngcM9OqfDEwOVS5jjDEnZleOG2OM8YsVDmOMMX6xwmGMMcYvVjiMMcb4xfUrx4NBRH4GNgY4eR1gRznGKS+Wyz+Wyz+Wyz+VMVcjVa3rS8NKWTjKQkSy1dfL7kPIcvnHcvnHcvnnZM9lu6qMMcb4xQqHMcYYv1jh+L1X3Q5QAsvlH8vlH8vln5M6lx3jMMYY4xfb4jDGGOMXKxzGGGP8clIWDhEZIiLLReSYiJR46pqI9BeRVSKSIyL3eg1vIiLficgaEfnAue17eeSqJSLTnPlOE5GaxbTJcO4ofLw7JCIXOePGiMh6r3EdQ5XLaXfUa9kTvYa7ub46isi3zvu9REQu9xpXruurpM+L1/hY59+f46yPxl7j7nOGrxKRc8uSI4Bcd4rICmf9zBCRRl7jin1PQ5RruIj87LX8673GXeO872ucZ/uEMtdTXplWi8ger3FBWV8i8qaIbBeRZSWMFxF51sm8REQ6e40r/3Xl6xOfKlMHnAK0AmYBXUpoEwmsBZoCMcBioI0z7kNgqPP6ZeDmcsr1GHCv8/pe4NETtK8F7ALinf4xwOAgrC+fcgH5JQx3bX0BLYEWzusUYAuQWN7rq7TPi1ebW4CXnddD8TwyGaCN0z4WaOLMJzKEuTK8PkM3H89V2nsaolzDgeeLmbYWsM75W9N5XTNUuYq0vw14MwTrqzfQGVhWwviBwBRAgNOB74K5rk7KLQ5VXamqq07QrBuQo6rrVPUwkAkMEhEBzgLGOe3GAheVU7RBzvx8ne9gYIqW/DCs8uJvrl+5vb5UdbWqrnFe5wHbAZ+ujvVTsZ+XUvKOA/o662cQkKmqBaq6Hshx5heSXKqa5fUZmofnaZvB5sv6Ksm5wDRV3aWqu4FpQH+Xcl0BvF9Oyy6Rqn6F50diSQYBb6nHPCBRRJIJ0ro6KQuHj1KBTV79uc6w2sAe9Txoynt4eUhS1S0Azt96J2g/lN9/aB9yNlWfEpHYEOeKE5FsEZl3fPcZYbS+RKQbnl+Ra70Gl9f6KunzUmwbZ33sxbN+fJk2mLm8XYfnl+txxb2nocx1qfP+jBOR44+SDov15ezSawLM9BocrPV1IiXlDsq6cuVBTqEgItOB+sWMul9VJ/gyi2KGaSnDy5zL13k480kG2uF5QuJx9wFb8Xw5vorn8bz/CmGuhqqaJyJNgZkishTYV0w7t9bX28A1qnrMGRzw+ipuEcUMK/rvDMpn6gR8nreIXAV0Ac70Gvy791RV1xY3fRByfQa8r6oFInITnq21s3ycNpi5jhsKjFPVo17DgrW+TiSkn61KWzhU9ewyziIXaODVnwbk4bmBWKKIRDm/Go8PL3MuEdkmIsmqusX5otteyqwuA8ar6hGveW9xXhaIyGjgrlDmcnYFoarrRGQW0An4GJfXl4hUBz4H/uZsxh+fd8DrqxglfV6Ka5MrIlFADTy7H3yZNpi5EJGz8RTjM1W14PjwEt7T8vgiPGEuVd3p1fsa8KjXtH2KTDurHDL5lMvLUOBW7wFBXF8nUlLuoKwr21VVsgVAC/GcERSD50MyUT1HnLLwHF8AuAbwZQvGFxOd+fky39/tW3W+PI8fV7gIKPYMjGDkEpGax3f1iEgdoCewwu315bx34/Hs//2oyLjyXF/Ffl5KyTsYmOmsn4nAUPGcddUEaAHML0MWv3KJSCfgFeBCVd3uNbzY9zSEuZK9ei8EVjqvpwLnOPlqAufw2y3voOZysrXCc7D5W69hwVxfJzIRuNo5u+p0YK/zwyg46yoYZwCEewdcjKcSFwDbgKnO8BRgsle7gcBqPL8Y7vca3hTPf+wc4CMgtpxy1QZmAGucv7Wc4V2A173aNQY2AxFFpp8JLMXzBfgOkBCqXEAPZ9mLnb/XhcP6Aq4CjgCLvLqOwVhfxX1e8Oz6utB5Hef8+3Oc9dHUa9r7nelWAQPK+fN+olzTnf8Hx9fPxBO9pyHK9TCw3Fl+FtDaa9prnfWYA4wIZS6n/0HgkSLTBW194fmRuMX5LOfiORZ1E3CTM16AF5zMS/E6WzQY68puOWKMMcYvtqvKGGOMX6xwGGOM8YsVDmOMMX6xwmGMMcYvVjiMMcb4xQqHMcYYv1jhMMYY4xcrHMYEmYh0dW7UFyciVcXzbJC2bucyJlB2AaAxISAi/8Fz5XgVIFdVH3Y5kjEBs8JhTAg49z1aABwCeuhv76hqTIViu6qMCY1aQAJQDc+WhzEVlm1xGBMC4nn+dCaeB/8kq+oolyMZE7BK+zwOY8KFiFwNFKrqeyISCcwVkbNUdeaJpjUmHNkWhzHGGL/YMQ5jjDF+scJhjDHGL1Y4jDHG+MUKhzHGGL9Y4TDGGOMXKxzGGGP8YoXDGGOMX/4P2rLKfyI/+soAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f93a5111f50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define dynamics\n",
    "def dynamics(x):\n",
    "    return -x**3 - x\n",
    "\n",
    "# Do some visual inspection of this system\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from pydrake.all import sin\n",
    "x_sample = np.linspace(-1, 1, 1000)\n",
    "f = np.array([dynamics(x_i) for x_i in x_sample])\n",
    "plt.plot(x_sample, f)\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"x dot\")\n",
    "plt.title(\"xdot = \" + str(dynamics(Variable(\"x\"))));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Verified globally stable!\n",
      "Resulting Lyapunov function:  0.010768550023664367*x(0)^2 + -8.3573990004666846e-19*x(0)^3 + 0.012087454329122366*x(0)^4 + -2.5625928704208748e-14*1 + -2.7488648052396447e-18*x(0)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEWCAYAAABmE+CbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzsvXd8XNWZ//9+1IvVm2VZ7g3bgImN6aCAMbApkAQSshtCFghJvmF/+91sCdlsypeQhLTNbnrYhFA2CRASwEkcHBsQvdhgjDEukrtkyeq9a57fH3PHDGIkjayZuVOe9+s1L82999w7nxmdmc895znnPKKqGIZhGEaS2wIMwzCM6MAMwTAMwwDMEAzDMAwHMwTDMAwDMEMwDMMwHMwQDMMwDMAMIe4QkfUi8kiQZV8WkRXh1mQYwSAinxCRZ6d5jd+KyFVBlDtNRJ6fzmvFI2YIMYSIbBKR2wLsv1JEGkUkBfgGcEeQl/wu8I7rGUYoCLK+nuy1q0XkpjH7TgNOBx6d7HxVfR3oEJH3nayGeMQMIba4G7hORGTM/uuAXwNnAHmq+mKQ19sAvFtEykMn0TBOcDcT1FdVHQnx633KuW6ws21/7ZxjOJghxBaPAIXABb4dIlIAvBe4F7gCeMrv2Lki0iIilc726SLSISLLAFR1AHgFWB+5t2AkEBPWVxEpEpENItIlIi8DC/1PdurvVhHpdP6e6+z/unPNH4lIj4j8yDllbP3/qYg85Lf9LRF53M+gqoFLRCQ91G88VjFDiCFUtR94EPi43+4PA3tUdQdwKrDXr/zzwM+Be0QkE7gP+A9V3eN3/m68zWzDCClB1NcfAwNAOXCD8wBARAqBPwM/AIqA/wT+LCJFqvpF4BngFlWdoaq3iEg2MB+/+g/8M3CaE5u4ALgRuN7XglDVemAYWBr6dx+bmCHEHvcA1zg/8OD9st3jPM8HuseU/yqQB7wMHMP7JfSn2znPMMJBwPoqIsnAh4Avq2qvqr7BW/UY4D1Ajarep6ojqvpbYA8wXp+/rw6fqP+q2gd8DK+Z/C/wD6paN+Y8q/9+mCHEGKr6LNAMXCkiC4Azgd84h9uBnDHlh/H25a4EvhegfzUH6AinZiNxmaC+lgApwFG/4of9ns8as+07XjHOS/nq8Nj6/zJwABC8rZWxWP33wwwhNrkX753WdcBfVfW4s/91YIl/QRGpAL4C/Ar4XoD+0lOAHeGVayQ4geprMzACVPqVm+P3/Bgwd8x15gD1zvO33dioai+wn3fW/88C6c71/m3MsVlAGm/vZkpozBBik3uBdcAneXszeyNwkW/DCZ7dDfwSb/9pA/A1v+PpwGpgc9gVG4nMO+qrqo4CfwC+KiJZIrIcuN7vnI3AEhH5WxFJEZGPAMuBPznHjwMLxrzO2Pq/BLgdb7fRdcC/icgqv/JVwBOqOhiSdxkPqKo9YvCBd4REO5A+Zv9W4Czn+T/ibTWkOduz8N6ZXeBsXwP8we33Yo/4fwSqr3i7jf4EdOGNcX0NeNbv+Pl4R8F1On/P9zt2DrDPueYPnH0rgV14u4dSnGve6nfOZ4CdPg14g9bvd/uziaaHOB+MESeIyHrg/6hqMLM1XwJuVG9AzzBiHhH5DfCgqk44W19ETgXuVNVzIqMsNjBDMAzDMACLIRiGYRgOZgiGYRgGYIZgGIZhOJz0aoNuUFxcrPPmzQt4rLe3l+zs7MgKimIdED1aokUHTKzllVdeaVHVkghLiol6DdGjJVp0QGxoeeWVV7qAF1T18kkv4vYwp6k8Vq9erePx5JNPjnsskkSLDtXo0RItOlQn1gJsU6vX4xItWqJFh2psaJlKvbYuI8MwDAOwGIJhGIbhYIZgGIZhAGYIhmEYhoMZgmEYhgGYIRhGQETkchHZKyK1InJrgOPpIvKAc/wlEZnnd+wLzv69InJZJHUbxnQwQzCMMTjZvH6MN0fvcuCjzvLM/twItKvqIuD7wLecc5cD1wIrgMuBnzjXM4yoJ6YmpgWitWeQe184TPGAx20pRpTz65cOc/jYCFWTF10L1KrqAQARuR+4EnjTr8yVeNOTAjyEN+G7OPvvV+8a+wdFpNa53guheRfxjarS0TdMU/cgx7sG6Owfpn9olN6hEfqGRhn1KIcODbFLaxGBtOQkcjNTyc1IJTczhaLsdGblZ5CTker2W4lJYt4QhkY9/PfjNVy/PM1tKUaUc+/zh8nSkWCKVvD21I51wFnjlVHVERHpxJsMvgJ4ccy570j7KCI3AzcDlJWVUV1dHVBIT0/PuMciTai1tA94qOnwcLTbQ0OPh2M9Hpr6lZFg7u1qJk5ylpUCRZlJlGcLc3KTmJOTxLzcZHLTJTTiHeLt/xPzhlCak0FKktA6YMt4GxNzrKOfs8qC+kEIVGhsBRuvTDDnoqp3AncCrFmzRquqqgIKqa6uZrxjkWa6WnoGR3h6XzNbdh/npQNt1Hf0A5AkMK8omxVzZ/De4mzKcjMoy02nNCeDgqxUstJTyEpNJjMtmdTkJJ6srubCCy/Co8rgsIeugWHvo3+Elp5B6jv6qW/vp669j5qmHl7e139Cw9KyHM5dVMSFi0s4b1ExaSnT6zWPp/8PxIEhJCcJM/MyaOkfcluKEcV0DQzTPThCUWZQLck63p7rdzbenLyBytSJSAqQB7QFeW7CMDLq4al9zdy/9ShP7W1maNRDflYq5y4s4obz57NmbgHLynNITwk+zJKSJCd+yDNSk8nLmrh7qLNvmN2NXbx6pJ3na1v5zUtH+NVzh8jNSOGKleV88F0VrJ1fiLfHL7GJeUMAqMjPpK3D0qIa43PMuRstygzqS78VWCwi8/Emdb8W+NsxZTbgzQH8AnA13ty8KiIbgN+IyH/iTVm6GG8qx4Sib2iEX794hLueO0hD5wDFM9K57py5rF9exuq5BaQkR248S15WKmcvKOLsBUX8n6pFDI6M8lxtC3/a0cCfXj/GA9uOsrw8lxvOn8+Vq2aRGkFt0UZ8GEJBJjUN7W7LMKKY+nbHEDImNwQnJnALsAlIBu5S1V0ichvehcI2AL8E7nOCxm14TQOn3IN4A9AjwGfVm1A+IRge9XDP84f4SfV+2nqHOHdhEV953wouOaU0an5o01OSuXhZGRcvK6N/aJRHXqvnrmcP8i+/28FPnqzlXy9byuUrZyZkiyE+DCE/k/YBZWTUE9E7DyN2ONFCCMIQAFR1I7BxzL4v+z0fAK4Z59yvA18/Sakxy4sHWvnSI29Q09TDBYuL+b/rFrN6bqHbsiYkMy2Zj66dw7VnVrL5zeN8e9NePvPrVzlrfiHf+tBpzCuOjqWtI0Vc/HpW5GeiQGPXgNtSjCilvmOA1GQJ+SgTA4ZGPHxz426uvfNFBkZG+Z+Pr+HeG9ZGvRn4IyKsXzGTx/7xAr7+gZW8eayLy//7ae5+7iCaQHnn46OFUJAJeLsFZhdkuazGiEaOdfRTnpdJUgJ2A4STYx39fOq+V9hZ38nfnTWH/3jPcjLTYnceXkpyEn931lwuWVbGvz+8k6/+8U22Hm7n2x86jez0uPi5nJC4eIez8r2GcKyzf5KSRqJS39HPrPwMwAYfhIqddZ3ceM9W+odG+fl1q7lsxUy3JYWMmXkZ/PL6Nfz86QN8+7E91B7v4Z4b1jIzL8NtaWElbrqM4K3AoWGM5VhHPxX51noMFc/XtvDhn79AanISD33m3LgyAx8iwqcvWsi9N5xFfUc/H/rp8xxs6XVbVliJC0PISE0mJ40TE10Mw5/hUQ/HuwaoyI/vu7tI8fz+Fm64ZytzCrN4+LPnsnRmjtuSwsr5i4v57SfPpn94lGt+9gKH4tgU4sIQAIozkqjvsKCy8U4aOwfw6Ftdi8bJ8+qRdm64eyuVBVn8+pNnUZqTGCZ76uw8HvzU2XhU+dgvX6KxMz5/a+LGEIoyhfr2PrdlGFGIb8ipGcL0aO7z8Ml7tlGak8FvPnk2xTPS3ZYUURaV5nD3359Je+8Q19/1Mr2DQa2LFVPEjyFkCPUd/Qk1RMwIDt9gAzOEk6dncIT/enWA4VEPd33iTEpyEssMfJw2O5+fXbeamqZu/vWhHXH3exM/hpCZxMCwh/a+YbelGFHGMacrscIM4aT58iNvcKxH+enHVrOodIbbclzlgsUl3HrFMjbubGTjwfj6vYkjQ/COL7eRRsZY6tr7KchKjenx8W7y+1fq+MP2eq5alMp5i4rdlhMVfPKCBbzntHJ+XzPM63UdbssJGfFjCM6SBPUdFkcw3k5dex+VhTbk9GQ42tbHlx59g7PmF/K+hZZ0xoeI8I0PnEpemvBPD7zGwHB8LFcVN4ZQnOl9KzbSyBhLXXs/lTaDfcqoKl985A0E+P5HVtks7zHkZaZy46np7G/u5TubJk7YEyvEjSFkp0JWWrJ1GRlvw+NRZ0kTix9MlQ07jvH0vmb+7fJlFpAfh5XFyXzs7Dn86rmD7DrW6bacaRM3hiAizMrPPDHE0DAAmroHGRr1MNu6jKZEZ/8wt/3xTVZV5vOxs+e6LSeq+df1yyjISuMrj+7C44ntUUdxYwjgHUVis5UNf446c1OshTA1flJdS1vfEF//wEqSk6yraCLyslL5/OXL2Ha4nYe317stZ1rElSHMMkMwxlDnGILFEIKnvqOfXz13iA+cUcGKWXluy4kJrl49m9Nn5/G9v+6N6QBzUIYgIpeLyF4RqRWRWwMcTxeRB5zjL4nIPGf/pSLyiojsdP5e7HfOamd/rYj8QEKQnmh2QSZtvUNxOYPQODmOtnlvEKyFEDzfcwKk/7J+qctKYoekJOHfLl/Gsc4BfvPSEbflnDSTGoKIJAM/Bq4AlgMfFZHlY4rdCLSr6iLg+8C3nP0twPtU9VS8+Wfv8zvnp8DNeHPOLgYun8b7ADgxtLDOAsuGQ117HyU56WSk2hyEYKht6ubh1+r5+/PmWSB5ipy3qJhzFxbx4ydr6YnRm9JgWghrgVpVPaCqQ8D9wJVjylwJ3OM8fwi4REREVber6jFn/y4gw2lNlAO5qvqCeud+3wtcNd03M8cxhCNtNhfB8HK0rZ9Kax0EzU+q95ORksynLlzotpSY5F8vW0pr7xD3vXDYbSknRTAJciqAo37bdcBZ45VxEpR3AkV4Wwg+PgRsV9VBEalwruN/zYpALy4iN+NtSVBWVkZ1dXVAkT09PXTt3g7Aky+/TmqTO5Noenp6xtUYaaJFi5s6ahv6WJifdOL1J9MiIoXAA8A84BDwYVVtD1DueuA/nM3bVfUeEckCfgcsBEaBP6rqO7pYo5WjbX08+toxrj9nHoXZaW7LiUnOmFPA+YuKueu5g9xw/jzSU2KrZRqMIQTq2x87tmrCMiKyAm830vopXNO7U/VO4E6ANWvWaFVVVUCR1dXVXHTRRXzh2U2kFc6iqmpFwHLhprq6mvE0Rppo0eKWjpFRD21/fYwPL5tHVdWyYLXcCjyuqnc48bJbgc/7F3BM4yvAGrz19hUR2YA3Hdt3VfVJEUkDHheRK1T1L6F+b+HgzqcPkCTwyQvnuy0lpvn0RQv52C9f4pHt9XzkzDluy5kSwXQZ1QGVftuzgWPjlRGRFCAPaHO2ZwMPAx9X1f1+5WdPcs0pIyJUFmadGFliJDaNXQOMenSqI4z8uz/vIXBX5mXAZlVtc1oPm4HLVbVPVZ8EcLpXX+Xt9Txqae8d4sFtR/ngGbMpz7Mutulw3qIiVszK5edPH4i5eQnBGMJWYLGIzHfueq4FNowpswFv0BjgauAJVVURyQf+DHxBVZ/zFVbVBqBbRM52Rhd9HHh0mu8F8AaWLYZggP8IoykZQplTP331tDRAmUDdqG/r8nTq/vuAx6fy4m7x4LajDI54uOF8ax1MFxHh5gsXcKC5l+p9TW7LmRKTdhk5MYFbgE1AMnCXqu4SkduAbaq6AfglcJ+I1OJtGVzrnH4LsAj4koh8ydm3XlWbgM8AdwOZwF+cx7SZU5jFszUtqCohGMlqxDAn5iAUvv2Od926dTQ2NgY6JT/IS0/WRZoC/Bb4gaoeCHiBKcTGwh1/8ajyi6f7WVqQRMOeV2jYE7icxaTeyXhasj1Kbprwgz9vJ6kxMlnlQvG5BBNDQFU3AhvH7Puy3/MB4JoA590O3D7ONbcBK6ciNhgqCzLpHx6lpWcoYZN4GF6OtveTJLyjC2TLli0By4tIBzAqIuWq2uCMhgt0i1cHVPltzwaq/bbvBGpU9b/G0zaV2Fi44y9P7mmiuX8rX75qFVWnzxq3XKLHpAIxkZbrhvfw0+r9LF51VkRycYTic4mrmcoAc4ps6Knhpa6tj5m5GaSlTKma+3d/Xk/grsxNwHoRKRCRAryDJTYBiMjteGNo//ekhUeY+148TElOOpetmOm2lLji2jPnoMADL8fORLX4M4QTk9PMEBKduvb+k1nU7g7gUhGpAS51thGRNSLyCwBVbQO+hje+thW4TVXbnAEUX8Q7gfNVEXlNRG4KzbsJD8e7Bqje28RH1lRO1TiNSagszKJqSQn3bz3KyKjHbTlBEVSXUSzhCyAeaTVDSHSOtvdxzsKiKZ2jqq3AJQH2bwNu8tu+C7hrTJk6AscXopZHX6vHo/DBdwWcBmRMk4+cWcmn//dVntvfykVLStyWMylxd0uQkZpMaU76iVUujcRkaMRDY9eALWo3CX94tZ5VlfksKEnsPMnhomppKbkZKTwaI6ugxp0hgA09NbwrdqpiqTMnYHdDF3sau611EEYyUpN5z2nlPLarkb6h6F/fKC4NYU5h1okx6EZicqi1F4B5RWYI4/Hw9npSkoT3njb+yCJj+ly1qoK+oVE2v3ncbSmTEpeGUFmYRUNnP8MxEsgxQo8vhjS3KNtlJdGJx6P8cccxqpaW2LpFYebMeYVU5GfGRPKc+DSEgkw8iqXTTGAOtfaSnZZM8Qz7sQvE6/WdNHQO8DenlrstJe5JShLec1o5z9W20DUw7LacCYlLQ7BlsI3DrX3MKcq22erj8NgbjaQkCZcsK3NbSkJw2YoyhkeVJ/dE91IWcWkIlWYICc/h1l6LH4yDqrJpVyPnLCwiL8udZeITjTMqCyjJSWfTroBLpkQNcWkIvtmpNhchMRn1KEfb+k/MWjfeTk1TDwdbem1mcgRJShLWLy+jem9zVOdcjktDSEoS5hRmnRhpYiQWDZ39DI16mGcB5YA89kYjIrB+uXUXRZLLVsykb2iUZ2paJi/sEnFpCADzirI51GIthETkrRFG1kIIxOY3j3NGZT6luZFZhdPwcvaCInIyUtj8ZvR2G8WxIXhbCLGWoMKYPodsyOm4tPQMsrO+k4uXBUrzYISTtJQkLlxcwlP7mvGmko8+4tcQirMZHPFwvHvAbSlGhDnc2ktaShLldgf8Dp6paQbgoiVmCG5w0ZISjncNsvd4t9tSAhK3hjC/2Ht3eLDF4giJxuHWPioLMklKsiGnY3lqbzNF2WmsmJXrtpSE5EJngbun9ja7rCQwcWsIvv5jiyMkHodaey2gHACPR3mmpoULFhebWbrEzLwMls3M4al9ZggRZVZeJmkpSTbSKMFQVY609Vn8IAC7jnXR2jvERUujfxnmeOaiJSVsPdRG72D0LXYXt4aQlCTMLczikHUZJRTNPYP0DY3aCKMAPOUkfL9gsRmCm1y0pIThUeX5/a1uS3kHcWsI4A0sWwshsThsQ07H5ZmaFlbMyqV4huUad5PV8wpIT0niBTOEyDKvKIvDrX029DSBOGxDTgMyMDzK9qMdnLNgahnkjNCTnpLM6rkFvHjADCGi+IaeNnTZ0NNE4XBrL8lJQkV+pttSoorXjnYwNOLhLDOEqODsBUXsbuyio2/IbSlvI64NYb5zl3jY4ggJw4GWXmYXZFrC+DG8dKANEVg7r9BtKQZeQ1CFlw+2uS3lbcT1t2auby6CxREShgPNvSwoPvnuIhEpFJHNIlLj/C0Yp9z1TpkaEbk+wPENIvLGSQsJMS8dbGXZzFxb3TRKOL0yj/SUJF48YIYQMcpzM0hPSbKRRgmCx6McbOmZbsL4W4HHVXUx8Liz/TZEpBD4CnAWsBb4ir9xiMgHgZ7piAglQyMeXj3SzlnzrXUQLURrHCGuDSEpSZhblMVBm5yWEDR0DTAw7GFBybQCylcC9zjP7wGuClDmMmCzqrapajuwGbgcQERmAJ8Dbp+OiFDyel0HA8Mezl5ghhBNRGMcIcVtAeFmblG2tRAShAPN3pvy+dPoMgLKVLUBQFUbRCTQoj8VwFG/7TpnH8DXgO8BE96FiMjNwM0AZWVlVFdXByzX09Mz7rFg+eN+7w/O8LE9VLfsPenrhEJLKIgWHTA9LWmdo6jC3X96mlWl0/8pDsXnEveGML84m6f2NTPqUZJtun5cc6DZa/wLJ+kyWrduHY2NAZcgzg/ypQJVJBWRVcAiVf0nEZk30QVU9U7gToA1a9ZoVVVVwHLV1dWMdyxY7j74MotL+3nf+oumdZ1QaAkF0aIDpqdl7dAI39n2V0byKqmqWuqqFh9xbwgLS7IZGvFQ324ZtOKdA809ZKclU5oz8cSrLVu2BNwvIh3AqIiUO62DciBQEtw6oMpvezZQDZwDrBaRQ3i/W6UiUq2qVbiEqvLa0Q4uW27Z0aKNrLQUTinP4dUj7W5LOUFcxxDgrbvF/c1RE+MzwsSBll4WlMxAZFotwQ2Ab9TQ9cCjAcpsAtaLSIETTF4PbFLVn6rqLFWdB5wP7HPTDMCbG6Kjb5gz5gTb+DEiybvmFLDjaAejUTJ5NmEMobbJDCHeOdDcO92AMsAdwKUiUgNc6mwjImtE5BcAqtqGN1aw1Xnc5uyLOrY7d59nzAk4etZwmXfNKaB3aJS9jdGRHyHuu4wKstMoyk6zFkKc0z80Sn1HPx8urpzWdVS1FbgkwP5twE1+23cBd01wnUPAymmJCQHbj3QwIz2FRaXTGoprhIl3OUb96pF2lkdBjoq4byGAt5VghhDf+BIhhaCFEFdsP9rO6ZV5NqAiSqkszKQoOy1q4giJYQil2dZlFOccaPH+f80Q3qJ/aJQ9Dd2sqrT4QbQiIpwxp4DtRzrclgIEaQgicrmI7BWRWhEJNHMzXUQecI6/5BtyJyJFIvKkiPSIyI/GnFPtXPM15xG2JK8LS2bQ3jdMW2/0TAAxQotvyOk05yDEFW8c62TEo5xRafGDaOaMOfkcbOmls2/YbSmTG4KIJAM/Bq4AlgMfFZHlY4rdCLSr6iLg+8C3nP0DwJeAfxnn8n+nqqucR6DhfSFhYamNNIp3DjT3MCsvg6y0uA+LBY0voLzKRhhFNadW5AFeA3ebYFoIa4FaVT2gqkPA/Xin9/vjP93/IeASERFV7VXVZ/Eag2ss8g09tW6juMU35NR4ix11ncwuyLSEOFHOCUOod98QgrmdCjRN/6zxyqjqiIh0AkVAyyTX/pWIjAK/B25X1XcMxg3FFH+PKqlJ8OSru5nZd2ASSdMjXqbVx5IOVWVfQx/nVaRM+jrR8plEgjePdbFyVp7bMoxJKMhOoyI/k50xYggBp+mfRJmx/J2q1otIDl5DuA649x0XCdEU/0WvP8NQejpVVWsnkTU94mVafSzpONbRz8CmJ6g6YylV58xzVUu00D0wzMGWXj54RsXkhQ3XObUiLypaCMF0GdUB/oO7ZwPHxisjIilAHjDhRB1VrXf+dgO/wds1FTYWlmSzv9kWuYtHapyuwMVlOS4riR52N3gnOq2ocH9suzE5KytyOdTaR9eAu4HlYAxhK7BYROaLSBpwLd7p/f74T/e/GngiUPePDxFJEZFi53kq8F4grMlEFpXO4Gh7HwPDo+F8GcMFao57f/wW2+SrE+xyApQrrMsoJljpxBF21Xe5qmNSQ1DVEeAWvOu37AYeVNVdInKbiLzfKfZLoEhEavGuBX9iaKqz0Nd/Ap8QkTpnhFI6sElEXgdeA+qB/wnd23onC0tmoPrWBCYjfth3vJui7DSKLHh6gl3HuiiekTbpQn9GdBAtgeWgxuip6kZg45h9X/Z7PgBcM86588a57OrgJIYG/0XuTim3ZnQ8UdPUY0szjGHXsS5WzMqb7kJ/RoQompHOrLwM1wPLCTFTGbwzWEWg5rgNPY0nVJXa4z0ssfjBCQZHRqk53s2KKFgbxwieFRV5J7r63CJhDCEjNZl5RdnUNEXHqoJGaGjsGqB7cIQlZdZC8LGvsYcRj1r8IMZYNjOHQ63uxjkTxhAAlpTNYE+ULDNrhIZ9TotvUam1EHy8FVC2FkIssWxmLqMedXXdtYQyhKVlORxq6bWRRnGEb4SRtRDeYndDF9lpycwptAyBscTSmd6bGjdzIySUISyZmYNHbU2jeKLmeA+FNsLobew93s2SmTkk2ZLXMcW8oizSUpLYe9wMISIsdQKP+1z8wI3QUtPUbfMPxlBzvOdEXTdih5TkJBaXututnVCGMK84m7TkJIsjxAmqSs3xHhZbd9EJWnoGae0dslnbMcrSmTnsbXRvclpCGUJqchILSrLZZ4YQFxzvGnRGGNmPnw9f3bYWQmxyysxcjncN0u5S7paEMgTwOvA+m4sQF/i6/kI5KU1ECkVks4jUOH8DZpcRkeudMjUicr3f/jQRuVNE9onIHhH5UMjEBYHvM1ky01pNsYgvsOxWL0bCGcKSshzqO/rpdnkRKWP67HGa1qfMDOnwyluBx1V1MfA4fsuw+BCRQuAreJeBXwt8xc84vgg0qeoSvAmlngqluMnYe7yH/KxUSizIHpMsOzHSyJ1uo4QzBN8HboHl2Gd3QzczczMoyE4L5WX9kz3dA1wVoMxlwGZVbVPVdmAzcLlz7AbgmwCq6lHVyXKChJR9x7tZUpZjS1bEKCU56RRkpbo20ijh8g36+pv3Nvawem6hy2qM6bC7oYtl5SHvKy9T1QYAVW0YJ9d3oKRRFSLiy1X5NRGpAvYDt6jq8bEXCEXip7GoKm/W93HOrMkTBZ0M0ZJcKFp0QHi0lKaPsnVfPdXVE2YQCIuWhDOEivxMstOSrYUQ4wyNeKht6uHdywL9Xk/MunXraGxsDHQo2OTD4yWESsGbL+Q5Vf2ciHwO+C7e5E9vLxyixE/+NHT207/pCd5YYvZ4AAAgAElEQVQdRKKgkyFakgtFiw4Ij5bN7Tv50+sNXHTRRVNq6YVCS8IZQlKSsLgsx9XZgMb0qW3yrtdzMivXbtmyJeB+EekARkWk3GkdlANNAYrWAVV+27OBaqAV6AMedvb/DrhxygJPEl+dtlFXsc3Ckhl09g/T1jsU8QmXCRdDAO+QvL3Hu5kgh48R5bwVUA75j59/sqfrgUcDlNkErBeRAieYvB7Y5CSF+iNvmcUlwJuhFjgevpV8zRBim4WlvqX6I5+7JSEN4ZTyHNp6h2jqHnRbinGS7G7oIi0lifnF2aG+9B3ApSJSA1zqbCMia0TkFwCq2gZ8DW82wa3Abc4+gM8DX3WSP10H/HOoBY7H/uYeirLTQh1kNyLMwhJvnXZjkbuE6zICWO4sC7zrWCdluRkuqzFOhj2N3Swpm0FKcmjvaVS1Fe+d/dj924Cb/LbvAu4KUO4wcGFIRQXJgZbecBikEWFm5WWSkZrkypprCdtCAPfzlxonz+6GrlDPP4h5DjT3sqDEDCHWSUoSFhTPMEOIFDkZqcwrymLXMTOEWKSpe4CWniFLhepH18AwLT2DLCixGcrxwMJSM4SIsmJWHrsa3E1XZ5wcexq8o2nCMAchZjngBCAXWJdRXLCwJJu69v6I525JWENYPiuXo239dPbbEhaxxu4Gb8tuubUQTnCwxXs3aV1G8cHCkhmowsGWyI40SmhDgLd+XIzYYXdDF+V5GeRn2WgaHweae0lOEuYUmiHEAwtLfENPI9ttlLCG4Ms3a3GE2GNnfaclkB/DgeZeKgsySUtJ2K90XDG/OBsR2N9kLYSIUJqTQUlO+omE5EZs0DM4woGWXk6tMEPwZ39zjwWU44jMtGQq8jOthRBJVszK5U1rIcQUbx7rQhVOnW3xAx8ej3Ko1eYgxBvzi7M53GothIixYlYuNU09EY/kGyfPznpvi26ltRBO0NA1wMCwxwLKccbcoiwOtfZF9DUT2hCWl+cx6lFb+TSGeKO+k7LcdEpzbIa5jwNOt8KCYusyiifmFmbT2T9MR1/k0mkmtCFYYDn22FnfafGDMfiGJloLIb6YW5QFwOEIthIS2hDmFGaRk5FyohvCiG56BkfY39xj3UVjONLaR0ZqEqU5ljYznpjnxIQORTCOkNCGkJQknD47nx1HO9yWYgTBiYCyGcLbONzWx5zCLEubGWfMKbQWQsQ5vTKPPY3dFliOAXwtOTOEt3O0rc8mpMUhGanJzMzNMEOIJKfPzmfUozYfIQZ4o76T0px0Sm3J8hOoKkecFoIRf8wtyoro0FMzhEpvGt0dR80Qoh0LKL+Tlp4h+oZGTwQgjfgi0kNPE94QynIzmJmbwY46iyNEM76A8qmzzRD8OdLmvXu0FkJ8Mrcom5aeQXoHRyLyekEZgohcLiJ7RaRWRG4NcDxdRB5wjr8kIvOc/UUi8qSI9IjIj8acs1pEdjrn/EBcjIidXplngeUoZ8fRDlThjDkFbkuJKo60ee8e51gLIS6ZV+SNDUUqjjCpIYhIMvBj4ApgOfBREVk+ptiNQLuqLgK+D3zL2T8AfAn4lwCX/ilwM7DYeVx+Mm8gFJxemc+h1r6ITgAxpsb2I+0ArJqd77KS6OJIaz8iMLsg020pRhjwdQX6WoLhJpgWwlqgVlUPqOoQcD9w5ZgyVwL3OM8fAi4REVHVXlV9Fq8xnEBEyoFcVX1BVRW4F7hqOm9kOpzu/Mi8XmdxhGhl+5EOFpZkk5eVGtbXEZFCEdksIjXO34BNEhG53ilTIyLX++3/qNPyfV1EHhOR4nDqPdzWS3luBukpyeF8GcMlfC2/SMURUoIoUwEc9duuA84ar4yqjohIJ1AEtExwzbox16wIVFBEbsbbkqCsrIzq6uqAF+zp6Rn32GT0DSsAjzyzHc+x6a2xPx0doSZatExXh6ry8v4+VpWmTPv9BKHlVuBxVb3D6R69Ffi8fwERKQS+AqwBFHhFRDYA3cB/A8tVtUVEvg3cAnx1WqIn4GhbH5UWP4hbcjNSKcxOi9hIo2AMIVDfvp5EmZMqr6p3AncCrFmzRquqqgJesLq6mvGOBcPC16vpSsmmqurMk75GKHSEkmjRMl0dh1t76d5UzRVrT6HqrDnh1nIl4CtwD1DNGEMALgM2q2obgIhsxtvl+RDeup0tIq1ALlA7LcGTcLi1j6qlJeF8CcNlKgsyqWvvj8hrBWMIdUCl3/Zs4Ng4ZepEJAXIA9omuebsSa4ZUVZVFlC9twlVtRmfUcarTvzgjDkRiR+UqWoDgKo2iEhpgDKBWs0VqjosIp8BdgK9QA3w2UAvEoqW7+Co0tQ9yEjn8Yi0BOOlxRlKIqElbWSAmnrPpK8TCi3BGMJWYLGIzAfqgWuBvx1TZgNwPfACcDXwhBMbCIjzResWkbOBl4CPAz88Cf0hY828An7/ah0HW3ot0UiUsf1IB9lpySwpywnJ9datW0djY2OgQ8E6TsAWroikAp8BzgAO4K3TXwBuf0fhELR89x3vhs1Pc9HqFVStCtjjGlLipcUZSiKh5YW+3ex4/hAXXngRSUnj36yGQsukhuDEBG4BNgHJwF2quktEbgO2qeoG4JfAfSJSi7dlcK3vfBE5hLfpnCYiVwHrVfVNvF+cu4FM4C/OwzXWzPXGDrcdbjdDiDK2H+ng9Mp8kif4MkyFLVu2BNwvIh3AqIiUOzct5UBTgKJ1vNWtBN4WbjWwCkBV9zvXexBvDCIsHHECjTYHIb6pKMhkaMRDS89g2GfpB9NCQFU3AhvH7Puy3/MB4Jpxzp03zv5twMpghYabhSUzyMtM5ZVD7Xx4TeXkJxgRoX9olN0NXXzqogWReklfa/cO5++jAcpsAr7hNwJpPd6WQAawXERKVLUZuBTYHS6hde1eQ7CgcnzjG1J8tL0/7IaQ8DOVfSQlCavnFrDt8EShDyPS7KzvZMSjnFEZsQlpdwCXikgN3h/0OwBEZI2I/ALACSZ/DW936lbgNlVtU9VjwP8DnhaR1/G2GL4RLqH1Hf1kpCZRlD29kXFGdDO7wGv49R3hDywH1UJIFFbPLeCJPU209w5RYF+yqMBn0BEKKKOqrcAlAfZvA27y274LuCtAuZ8BPwunRh/1Hf3Mys+0QRBxTkW+t4XgaxGGE2sh+OGLI7xyuN1lJYaPlw+2sbh0BkUzLPnLWOo7Bk78WBjxS3Z6CgVZqREZemqG4MfplfmkJgvbzBCiglGPsu1QO2vnF7otJSqpb+83Q0gQZhdkmSFEmozUZFbMyuMViyNEBbsbuugZHDFDCMDA8CgtPYNmCAnC7IJM6zJygzVzC9hR18ngiGVQc5uXDnqN2QzhnRxzAowVtqhdQjC7IJP69n4mmN4VEswQxrBmXiFDIx5b6C4K2HqwjcrCTMrz7EdvLMc6vOtFWgshMZhdkMXgiIeWnvCuyGyGMIaz5hciAi/sb3VbSkKjqrx8qI2184rclhKV1Hd4uw9mmSEkBJEaaWSGMIaC7DROmZlrhuAy+5t7aOsd4izrLgpIfXs/SQIz8yy/dCIwu9BnCOENLJshBOCchUW8cqSdgWGLI7iFxQ8mpq6jn5m5GaQm21c4EXirhWCGEHHOXVjE0IjnxCqbRuR5+WAbJTnpljx+HI519FtAOYHIyUglLzP1RFdhuDBDCMCZ8wtJEnjRuo1cQVV5rraVcxcW2SzccfDNUjYSh/K8DBo7ByYvOA3MEAKQm5HKqRV5vHDADMEN9h7vpqVnkPMWhTX7ZMwy6lEabJZywlGel0GDGYI7nL2wiNeOdtA3NOK2lITj2Rpv5tXzzRAC0tQ9wIhHrcsowZiZl2ktBLc4d2Exw6PepROMyPJsbQsLSrKtS2QcTkxKs88noSjPy6C1dyisg13MEMbhzHkFpCQJz+1vcVtKQjE04uGlA23WOpgA36Q0M8zEotwZYny8K3ytBDOEcchKS2HNvAKe3meGEEm2H2mnf3jU4gcT4Os2KAtzshQjuvDN2A9nHMEMYQIuWlLK7oausDqy8XaerW0hSbxzQYzANHYNkJWWTG6GpTNJJHyTEMMZRzBDmICLlpQA8PS+ZpeVJA7P1rZwemU+uRmpbkuJWhq7BpiZm2FDchMMX5fRsc7wTU4zQ5iAU8pzKMlJ5ykzhIjQ3jvEjqMdXGDdRRPS2Dlg3UUJSHZ6CrkZKdZCcAsR4aIlJTxT08LIqMdtOXHPU/ua8Si8e1mp21KimsbOgRN3i0ZiUZ6XaTEEN6laWkJn/zA7bDnssPPEniaKstM4fXZk8ifHIh6P0tQ9QJkZQkJSnp9Bg3UZucf5i4pJEqzbKMyMjHp4al8zFy0tISnJvb5xESkUkc0iUuP8LRin3GMi0iEifxqzf76IvOSc/4CIpIVSX2vvEMOjykzrMkpIwr18hRnCJORnpbGqMp+n9ja5LSWu2X60g87+YS52v7voVuBxVV0MPO5sB+I7wHUB9n8L+L5zfjtwYyjF+Ua82bLXiUl5XiYtPUNhy+hohhAEFy8rZUddJ002/DRsPLGnieQk4YLFJW5LuRK4x3l+D3BVoEKq+jjQ7b9PvMN+LgYemuz8k8XXf2wthMTEdyNwvHMwLNe3gcxBsH7FTL7713389c3jfOzsuW7LiUue3NPEmrkF5GW6Pty0TFUbAFS1QUSm0mQpAjpU1bcAVh1QEaigiNwM3AxQVlZGdXV1wAv29PS87dgzR4YBOLDrVdr3R/Z+bqwWt4gWHRB5LU0t3pbBX556gaWFySHXYoYQBItLZzCvKMsMIUzUd/Szp7GbL1yxLCKvt27dOhobGwMdmm40O1DwI2BWdFW9E7gTYM2aNVpVVRXwgtXV1fgf27ZpL8l79vP+9e8mOcKxlrFa3CJadEDktcxu6uG7255i5oJlVK16+71GKLSYIQSBiLB+xUx+9dxBugaGbdJUiHnsDe+P8/oVMyPyelu2bAm4X0Q6gFERKXdaB+XAVIJHLUC+iKQ4rYTZwLFpC/ajoXOA0pz0iJuBER2cmJzWEZ7ua4shBMn65WUMjyrVe220Uah57I0Gls3MYX5xtttSADYA1zvPrwceDfZEVVXgSeDqkzk/GI532aS0RCY7PYUZ6SlhW07HDCFIzphTQPGMNP66K2BXg3GSNHUPsO1wO5evjEzrIAjuAC4VkRrgUmcbEVkjIr/wFRKRZ4DfAZeISJ2IXOYc+jzwORGpxRtT+GUoxTV22aS0RKc0J53mbgsqu0pykrDulDL+9HoDgyOjpKckT36SMSmbdh1HFa5YWe62FABUtRW4JMD+bcBNftsXjHP+AWBtuPQ1dg7Y0uAJTmluurUQooH1K8roGRw5kdHLmD6PvdHAguJslpTNcFtK1NMzOELP4IjNQUhwSnMyaApTC8EMYQqcv6iEvMxUNuwIaZwwYWnvHeLFA21cvnKmrdwZBL4ZqtZllNiU5abT1D2AN2QVWswQpkBaShJ/c+pMNr95nP6h8KWxSxQ27Wpk1KNR010U7VhiHAO8LYSBYQ9dA6HP926GMEXed/os+oZG2bL7uNtSYp5HXqtnQXE2Kyty3ZYSEzR22SxlwxtDAGjuDn0cIShDEJHLRWSviNSKyDvWdhGRdGchr1pnYa95fse+4Ozf6zcSAxE5JCI7ReQ1EdkWijcTCc6aX0RpTrp1G02TYx39vHigjStXVVh3UZD4RpaU5KS7rMRwk9IcX27l0McRJjUEEUkGfgxcASwHPioiy8cUuxFoV9VFwPfxLvCFU+5aYAVwOfAT53o+3q2qq1R1zbTfSYRIThLee9osntrbTGf/sNtyYhafoV51xiyXlcQOzd2DZKclk51ugwMTGV8LocmlFsJaoFZVD6jqEHA/3gXA/PFfEOwhvGOzxdl/v6oOqupBoJYwDsmLFO9fNYuhUQ+PvdHgtpSY5ZHt9bxrTj5zi6JiMlpM0NQ9QKl1FyU8pU4LsSkMLYRgbjUqgKN+23XAWeOVUdUREenEOymnAnhxzLm+BTgU+KuIKPBzZ22Xd3Cyi4CFE1WlLEu464ldlPUecE3HZESLlrE6jnZ72NPYz3XL0yKuL1o+k5OhuXuQkhnWXZTozEhPISstOSxdRsEYQjALdo1XZqJzz1PVY85qkptFZI+qPv2Owie5CFi4+Ti1fGfTXuasWMOCkrfG0CfywlvB6vjGxt2kJB3knz50EYXZIc0fM2UtsURzzyCnzLQAfKIjIpTmpLvWZVQHVPptB1qw60QZEUkB8oC2ic5VVd/fJuBhYqwr6erVs0lOEh7cVue2lJhicGSUh16pY90pZRE3g1inuXvQAsoGAKW5GWHpMgrGELYCi53UgGl4g8QbxpTxXxDsauAJZ6GvDcC1ziik+cBi4GURyRaRHAARyQbWA29M/+1EjrLcDN69tJTfv1rH8KjHbTkxw+Y3j9PWO8S1aysnL2ycYGB4lO6BETMEA8C9FoKzjO8twCZgN/Cgqu4SkdtE5P1OsV8CRc6CXp/DSTuoqruAB4E3gceAz6rqKFAGPCsiO4CXgT+r6mOhfWvh5yNnVtLcPciTeyy9ZrDc//JRKvIzoyEzWkxhQ04Nf3zLV4R6tnJQ49dUdSOwccy+L/s9HwCuGefcrwNfH7PvAHD6VMVGG+9eWkJpTjoPbjsasbX8Y5nDrb08W9vC5y5dYuv5TxHf2jWlZggG3uUr+oZG6RkcISeE+VlspvI0SElO4urVs3liTxN17X1uy4l67t96lCSBD6+x7qKp4puVai0EA/znIoQ2jmCGME0+dvZcRIT7XjjstpSoZmB4lAe2HuXiZWW2WudJYF1Ghj9vzVYObRzBDGGazMrP5PIVM/nty0foGwr9YlPxwiPb62nrHeKG8+e5LSUmae4eJEmgKNsMwfB2GQEhT5RjhhAC/v68eXQNjPCHV+vdlhKVqCp3PXeQU8pzOWdBkdtyYpKm7kGKZlguZcNLidNCCPXQUzOEELB6bgGnVuRx9/OHwrJGeayzq3WUfcd7uPH8+baQ3Ulis5QNf3IzUkhLTqKl1wwh6hAR/v68edQ29bCj2fIkjGXTIe/4+fedbnkPTpbmHpuUZryFiFA0I42W7qGQXtcMIUS87/RZVORn8sf9w9ZK8OON+k52toxy/TlzYyIPtYgUishmEalx/haMU+4xEekQkT+N2f9rZ6n3N0TkLhEJyZjA5u5BG3JqvI3iGem09FgLISpJTU7iM1UL2d/p4YX9rW7LiRp+9EQtmSlw3Tnz3JYSLLcCj6vqYuBxZzsQ3wGuC7D/18Ay4FQgE7hpuoI8HrVlK4x3UDwjzQwhmrl69Wzy04UfPlHrtpSoYE9jF4/tamT93FTyMkM3eSbM+C/lfg9wVaBCqvo40B1g/0Z1wDsLf/Z0BXX0DzPiUTME420Uz0intSe0XUaWaSOEZKQmc8X8VH67p5Wth9o4c16h25Jc5YdP1DIjPYVL58aMGQCUqWoDgKo2OKvxThmnq+g64B/HOR70su4bn3gWgOaj+6mudm++S7QsHR4tOsBdLb1tQzR3D/Pkk08iIiHRYoYQYqpmp7C5TvjOY3t54FNnJ+yomr2N3Wzc2cBnqxYxIy26EgmtW7eOxsbGQIfyQ/gyPwGeVtVnAh2cyrLueRUr4LmXqTrrXayd795NRrQsHR4tOsBdLbXJB9h4cDdnnHUe+VlpIdFihhBi0lOEf1y3mC898gZP7GniklPK3JbkCt/8y25y0lO46YL5vPZydBnCli1bAu4XkQ5gVETKndZBOTDllQtF5CtACfCpaQl1sFnKRiB89aGlZ4j8rNAsJW8xhDBw7ZmVzC/O5luP7WHUk3gjjp6taaF6bzP/cPHikFXUCOK/lPv1wKNTOVlEbgIuAz6qqiFZF90MwQhE8QyfIYQusGyGEAZSk5P418uWsu94D79/JbES6Hg8yjc27mZ2QSYfP3eu23JOhjuAS0WkBrjU2UZE1ojIL3yFROQZ4Hd484fXichlzqGf4V3e/QUReU1Evsw0aesdIj0liey06B+2a0SOcBiCdRmFiStWzuSMOfl8e9MeLls5M5ZG2UyLh16t482GLn7w0TNiYt7BWFS1FbgkwP5t+A0hVdULxjk/5N+plp4himekJ2w8yghM0Qxv67slhOsZWQshTIgIX7tyJW29Q3zvr3vdlhMR2nqH+ObG3ayeW8B7T7VZyaGitXfQ0o0a76AgK40kgdbe0A09NUMIIysr8vj4OfO478XD7KzrdFtO2Pnmxt10D4zwjQ+cSpItwhYy2nqHTtwNGoaP5CShMDu0s5XNEMLM59YvoSg7nX9/eGdc515+6UArv3uljpsuWMDSmTluy4krWnuGbNlrIyDFM9JoDuF6RmYIYSY3I5XbrlzBzvpOfvLkfrflhIXewRE+//vXmV2QyT9esthtOXGFqtLaO2gtBCMgoV7PyAwhAvzNqeVcuWoWP3yiJi67jm7/824Ot/XxvWtOJ9NGwoSUwVEYGPZQZDEEIwDFM9JoDeES2GYIEeK296+kaEYan3vwNfqH4meJ7M1vHue3Lx/hUxcu5CxLfhNyuoa881iKLBeCEYDiGekhXQLbDCFC5GWl8t1rTqe2uYcvPrwzLpbIrmvv498e2sHy8lw+d+kSt+XEJd0+Q7AWghGA4px0+odH6R0MTfpeM4QIcsHiEv7vJUv4w/Z6fv3SEbflTIuB4VE+/b+vMOJRfvx37yItxapSOHirhWCGYLwT3+S0UK16at/iCPMPFy+iamkJ/++Pu9h6qM1tOSeFqvLvD+/kjfou/usjq5hfnO22pLjFZwg2D8EIhO9GoTlEgWUzhAiTlCT810dWUVmQxU33bKO2qcdtSVPmvx+v4Q+v1vNP65Yk7OJ9keKtLiOLIRjvpCTEy1eYIbhAflYa99ywltRk4fq7Xqapa8BtSUHzvy8e5r+21HDN6tn8f5cscltO3NM9qGSnJdvoLSMgvhaCdRnFOJWFWfzqE2tp7xvi2v95keMxYAoPb6/jS4++wcXLSvnmB0+1tXUiQNeQUmjxA2McfF2JbSEaemqG4CKnzs7j7r9fy/HOAa6980UaO6PXFH7z0hE+9+AOzllQxI//9l2kJFvViQTdQ9ZdZIxPekoyM9JTQraekX2rXWbt/ELuvXEtzd2DfOinz7OnscttSW9DVflJdS3//vBOqpaUcNcnzrTuiwjSNaQUWwvBmIDC7DTazBDih9VzC7n/5rMZ8Xi4+qcv8OTeKSfpCgsDw6N87sEdfPuxvbzv9Fn8/Lo1ZKSaGUSS7iG1EUbGhJghxCErK/J45LPnMacwixvu3sodf9nj6mJ4tU09XP2z53l4ez3/fOkSfnDtKptrEGFUla4htVnKxoQUZadZUDkeKc/L5PefOZdrz5zDz57az9U/fZ7dDZHtQhr1KPe+cIj3/vAZ6tv7+Z+Pr+EfLllsAWQX6BoYYVRtlrIxMaFsIVjGtCgjMy2Zb37wVC5YXMx/PPIG7/3hs3zi3Hn847rF5GaEN+vaK4fb+eqGXeys7+TCJSV89+rTKM3NCOtrGuPj+5LbLGVjIgpneA1BdfrfVTOEKOVvTi3n3IVFfHvTXu567iAPbjvKDefN54bz5pOXFVpjeOVwGz96opYn9zZTlpvODz56Bu87rdxaBS7T6kw2slFGxkQUZacxNOphIARrZgbVZSQil4vIXhGpFZFbAxxPF5EHnOMvicg8v2NfcPbv9UtEPuk1De8Etm984FT+9A/nc+7CIv778RrWfmMLn3vgNV7Y38rINGIMHX1D3PfCIa780bN86Kcv8NrRDv5l/RKe+Ocq3n/6rIQ1AxEpFJHNIlLj/C0Yp9xjItIhIn8a5/gPRWRa09BbnH5hCyobE1Ho3DD4ZrVPh0lbCCKSDPwYuBSoA7aKyAZVfdOv2I1Au6ouEpFrgW8BHxGR5cC1wApgFrBFRHzLYk52TcNhxaw8fn7dGnY3dPG/Lx5mw2vH+MP2evIyU7lgcTFr5hawoiKPxaUzyMtMfceP+cioh4bOAXY3dPFGfSfP1Law42gHHoVlM3P4yvuW85EzK8lKswYjcCvwuKre4dyo3Ap8PkC57wBZwKfGHhCRNUD+dIX4uoyKLahsTEBhtrfHoCsShgCsBWpV9QCAiNwPXAn4/3hfCXzVef4Q8CPx/ipdCdyvqoPAQRGpda5HENc0xnBKeS5f/8Cp/Md7lvPEniae3NvEU/ua+dPrDSfKpKUkUZydRnKyMNA/wOgzm2nvG8K32naSwGmz87nl4sWsX17Gyoo8l95N1HIlUOU8vweoJoAhqOrjIlI1dr9zA/Ud4G+BD0xHiK/LyFoIxkREtIUAVABH/bbrgLPGK6OqIyLSCRQ5+18cc26F83yyawIgIjcDNwOUlZVRXV0dUGRPT8+4xyJJpHRkA+8tgfcUJ9M+mMnhLg9NfUrHoNI1OIIHZTjZQ06Gh9yZqRRkCJU5SVTMSCIjZRg4RkvNMaprwi41av43EJSWMlVtAFDVBhEpneJL3AJscM4dt1Aw9bq3aYR3FSvPP/v0FCWEh2j5P0aLDogOLe0DHlaXJZM8MjBtLcEYQqBaPdaKxisz3v5AsYuA9qaqdwJ3AqxZs0arqqoCiqyurma8Y5EkWnRA9GiJFh3g1XL77bfT2NgY6PC0unlEZBZwDW+1MMYlmHpdRfR9dtGgJVp0QPRo+cDlodESjCHUAZV+27OBY+OUqRORFCAPaJvk3MmuaRhhYcuWLQH3i0gHMCoi5c4dfjkwlWnjZwCLgFqndZAlIrWqasvCGjFBMKOMtgKLRWS+iKThDRJvGFNmA3C98/xq4An15ojcAFzrjEKaDywGXg7ymobhBv51+Xrg0WBPVNU/q+pMVZ2nqvOAPjMDI5aYtIXgxARuATYBycBdqrpLRG4DtqnqBuCXwH1O0LgN7w88TrkH8QaLR4DPquooQKBrhv7tGcaUuQN4UERuBI7g7QLyjYYqKd8AAAU9SURBVBz6tKre5Gw/AywDZohIHXCjqm5ySbNhhISgxhmq6kZg45h9X/Z7PoDzxQlw7teBrwdzTcNwG1VtBS4JsH8bcJPf9gVBXGtGaNUZRnixtYwMwzAMwAzBMAzDcDBDMAzDMAAzBMMwDMNBVKc/3TlSiEgzcHicw8VASwTljEe06IDo0RItOmBiLXNVtSSSYiBm6jVEj5Zo0QGxoWUx8IKqXj7ZBWLKECZCRLap6hrT8RbRoiVadEB0aQmGaNIbLVqiRQfEnxbrMjIMwzAAMwTDMAzDIZ4M4U63BThEiw6IHi3RogOiS0swRJPeaNESLTogzrTETQzBMAzDmB7x1EIwDMMwpoEZgmEYhgHEsCGIyDUisktEPM5KlOOVu1xE9opIrZMjN9Q6gk3KPioirzmPkC71Pdl7dJYff8A5/pKIzAvl609BxydEpNnvc7gp0HVCoOMuEWkSkTfGOS4i8gNH5+si8q5w6DgZoqVeO6/hat2OlnodpJb4qNuqGpMP4BRgKd6ct2vGKZMM7AcWAGnADmB5iHV8G7jVeX4r8K1xyvWE6XOY9D0C/wf4mfP8WuABl3R8AvhRBOrGhcC7gDfGOf43wF/wZvQ7G3gp3JqmoD0q6rXzOq7V7Wip11PQEhd1O2ZbCKq6W1X3TlJsLVCrqgdUdQi4H28S9VByJd5k7Dh/rwrx9ScjmPfor/Eh4BKRCRL+hk9HRFDVp/Hm5RiPK4F71cuLQL6THc11oqheg7t1O1rqdbBaIkK463bMGkKQVABH/bbrnH2h5G1J2YHxkrJniMg2EXlRREL5xQrmPZ4oo6ojQCdQFEINweoA+JDTlH1IRCoDHI8EkagX4SRS+t2s29FSr4PVAnFQt4NKkOMWIrIFmBng0BdVNZjUhoHuFqY8znYiHVO4zBxVPSYiC4AnRGSnqu6fqpZA8gLsG/seQ/I5hEDHH4HfquqgiHwa793dxSHWEQyR+DzGf/EoqdeTaZnCZcJRt6OlXgf7OnFRt6PaEFR13TQvUQf4O/Vs4FgodYjIcQkiKbuqHnP+HhCRarwJ2UNhCMG8R1+ZOhFJAfKYuNkZFh3qzUbm43+Ab4VYQ7CEpF6cLNFSryfT4nLdjpZ6HZSWeKnb8d5ltBVYLCLzRSQNb+AppCN8CCIpu4gUiEi687wYOA9vnulQEMx79Nd4NfCEOhGoEDKpjjF9me8HdodYQ7BsAD7ujMg4G+j0dY3ECJGo1+Bu3Y6Weh2Ulrip2+GOiocx2v4BvG44CBwHNjn7ZwEbx0Td9+G9Y/liGHQUAY8DNc7fQmf/GuAXzvNzgZ14RyfsxJuQPZQa3vEegduA9zvPM4DfAbXAy8CCMP1PJtPxTWCX8zk8CSwLk47fAg3AsFNHbgQ+DXzaOS7Ajx2dOxlnNE8i1+toqNvRUq8TqW7b0hWGYRgGEP9dRoZhGEaQmCEYhmEYgBmCYRiG4WCGYBiGYQBmCIZhGIaDGYJhGIYBmCEYhmEYDmYIcYaInOkssJUhItnO2vor3dZlGNPF6nb4sYlpcYiI3I53FmcmUKeq33RZkmGEBKvb4cUMIQ5x1lvZCvz/7d27bUJBEAXQOyI1CRkhbTikECqgLkQVhK7JAdY4YBtAsHro6ZwKJrjS3U8wv0m+u/tv4ZHgLWR7Lk9G67RL8pVkm8dpCtZCtidyQ1ihsdf2muSQZN/d54VHgreQ7bk+eh8Cz6uqU5J7d1+qapPkp6qO3X1bejZ4hWzP54YAQBJ/CAAMCgGAJAoBgEEhAJBEIQAwKAQAkigEAIZ/1Re5yVjRnmwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f93a5111350>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from pydrake.all import (\n",
    "    MathematicalProgram,\n",
    "    Polynomial,\n",
    "    SolutionResult,\n",
    "    Variables\n",
    ")\n",
    "prog = MathematicalProgram()\n",
    "x = prog.NewIndeterminates(1, \"x\")[0]\n",
    "\n",
    "# Search over Lyapunov candidates\n",
    "(V, constraint) = prog.NewSosPolynomial(Variables([x]), 4)\n",
    "# V(0) = 0\n",
    "prog.AddConstraint(V.ToExpression().EvaluatePartial({x: 0.}) == 0.)\n",
    "# V radially unbounded\n",
    "prog.AddSosConstraint(V - 0.01*Polynomial(x**2))\n",
    "\n",
    "# Calculate Vdot = dV/dx dx/dt\n",
    "# Because V is a polynomial, its Jacobian is a Polynomial.\n",
    "# To keep Vdot as a polynomial, we convert the dynamics to a\n",
    "# polynomial as well, which will fail if you specify non-polynomial\n",
    "# dynamics.\n",
    "Vdot = V.Jacobian([x]).dot(Polynomial(dynamics(x)))[0]\n",
    "# Vdot negative semidef\n",
    "# (assume dynamics(0) = 0) so skip that constraint\n",
    "prog.AddSosConstraint(-1.*Vdot)\n",
    "\n",
    "result = prog.Solve()\n",
    "if result == SolutionResult.kSolutionFound:\n",
    "    print \"Verified globally stable!\"\n",
    "    V_result = prog.SubstituteSolution(V)\n",
    "    Vd_result = prog.SubstituteSolution(Vdot)\n",
    "    print \"Resulting Lyapunov function: \", V_result\n",
    "    x_sample = np.linspace(-1, 1, 1000)\n",
    "    V_eval = np.array([V_result.Evaluate({x: x_sample_i}) for x_sample_i in x_sample])\n",
    "    Vd_eval = np.array([Vd_result.Evaluate({x: x_sample_i}) for x_sample_i in x_sample])\n",
    "    plt.subplot(1, 2, 1)\n",
    "    plt.plot(x_sample, V_eval)\n",
    "    plt.xlabel(\"x\")\n",
    "    plt.title(\"V(x)\")\n",
    "    plt.subplot(1, 2, 2)\n",
    "    plt.plot(x_sample, Vd_eval)\n",
    "    plt.xlabel(\"x\")\n",
    "    plt.title(\"Vdot(x)\");\n",
    "else:\n",
    "    print \"Could not prove globally stable.\"\n",
    "    print \"Result: \", result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Automatic Differentiation (\"autodiff\", specifically forward-mode)\n",
    "\n",
    "Core idea: Don't bother storing the complete symbolic expression (which can be huge), but instead remember your derivatives with respect to a list of variables of interest.\n",
    "\n",
    "Every Autodiff object has a value, and a list of derivatives\n",
    "\n",
    "$$\n",
    "\\left[ [f], [\\dfrac{df}{dx_1}, \\dfrac{df}{dx_2}, ...] \\right]\n",
    "$$\n",
    "\n",
    "and these objects can be combined by overloading every relevant operator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Input (-1.000000, -1.000000), value 2.0, derivative [-2. -2.]\n",
      "Input (-1.000000, 0.000000), value 1.0, derivative [-2.  0.]\n",
      "Input (-1.000000, 1.000000), value 2.0, derivative [-2.  2.]\n",
      "Input (0.000000, -1.000000), value 1.0, derivative [ 0. -2.]\n",
      "Input (0.000000, 0.000000), value 0.0, derivative [0. 0.]\n",
      "Input (0.000000, 1.000000), value 1.0, derivative [0. 2.]\n",
      "Input (1.000000, -1.000000), value 2.0, derivative [ 2. -2.]\n",
      "Input (1.000000, 0.000000), value 1.0, derivative [2. 0.]\n",
      "Input (1.000000, 1.000000), value 2.0, derivative [2. 2.]\n"
     ]
    }
   ],
   "source": [
    "from pydrake.all import AutoDiffXd\n",
    "def inner_product(v):\n",
    "    # Expects a [n x 1] array input\n",
    "    return v.dot(v)\n",
    "for x in np.linspace(-1, 1., 3):\n",
    "    for y in np.linspace(-1, 1., 3):\n",
    "        x_ad = AutoDiffXd(x, [1., 0.])\n",
    "        y_ad= AutoDiffXd(y, [0., 1.])\n",
    "        f_ad = inner_product(np.array([x_ad, y_ad]))\n",
    "        print \"Input (%f, %f), value %s, derivative %s\" % (x, y, str(f_ad.value()), str(f_ad.derivatives()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Way more systems in Drake currently support Autodiff than Symbolic, though we hope to close that gap in the future!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}