Free Final Time Problem.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e683bd95",
   "metadata": {},
   "source": [
    "# Setup\n",
    "\n",
    "Dependencies:\n",
    "- System: python3\n",
    "- Python: jupyter, numpy, scipy, jax, matplotlib\n",
    "\n",
    "Example setup for a Ubuntu system (Mac users, maybe `brew` instead of `sudo apt`; Windows users, learn to love [WSL](https://docs.microsoft.com/en-us/windows/wsl/install-win10)):\n",
    "```\n",
    "/usr/bin/python3 -m pip install --upgrade pip\n",
    "pip install --upgrade jupyter numpy scipy matplotlib jax jaxlib\n",
    "jupyter nbextension enable --py widgetsnbextension  # necessary for interactive sliders to show up\n",
    "jupyter notebook  # from the directory of this notebook\n",
    "```\n",
    "Alternatively, view this notebook on [Google Colab](https://colab.research.google.com/github/StanfordASL/AA203-Examples/blob/master/Lecture-5/Free%20Final%20Time%20Problem.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "ca1a0327",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import minimize, Bounds\n",
    "from scipy.integrate import solve_bvp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "212a9cb4",
   "metadata": {},
   "source": [
    "## [Closed-form solution](https://www.wolframcloud.com/obj/schmrlng/Published/Optimal%20Control%20Example%20%28Free%20Final%20Time%29.nb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "93e37a67",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.477694926940431"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def free_final_time_analytical(a=1.0, b=1.0):\n",
    "    # Closed-form solution for optimal final time.\n",
    "    # See https://www.wolframcloud.com/obj/schmrlng/Published/Optimal%20Control%20Example%20%28Free%20Final%20Time%29.nb\n",
    "    return (1800 * b / a)**(1 / 5)\n",
    "\n",
    "\n",
    "free_final_time_analytical()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "ef87432c",
   "metadata": {},
   "source": [
    "## Simple example solving the two-point boundary value problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "id": "04b8b9d6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "       message: The algorithm converged to the desired accuracy.\n",
      "       success: True\n",
      "        status: 0\n",
      "             x: [ 0.000e+00  2.500e-01 ...  3.667e+00  4.000e+00]\n",
      "           sol: <scipy.interpolate._interpolate.PPoly object at 0x30c4c77e0>\n",
      "             p: None\n",
      "             y: [[ 0.000e+00  5.114e-01 ... -1.136e+00 -2.000e+00]\n",
      "                 [ 2.067e+00  2.003e+00 ... -2.358e+00 -2.876e+00]]\n",
      "            yp: [[ 2.067e+00  2.003e+00 ... -2.358e+00 -2.876e+00]\n",
      "                 [-0.000e+00 -5.114e-01 ... -1.136e+00 -2.000e+00]]\n",
      " rms_residuals: [ 1.504e-04  1.021e-04 ...  2.797e-04  2.387e-04]\n",
      "         niter: 7\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABcGElEQVR4nO3dd3gU5cLG4d/spkMKIYEACb1JDV0Q6YoNBVQQAcVe0GP3YINjRc/nsaKCqIAI0kRAQQRRQJEaCL33TmjpdXe+P0ZRpGUhm8kmz31dezkJk+wzbsI+zLzvO4ZpmiYiIiIiNnDYHUBERERKLhURERERsY2KiIiIiNhGRURERERsoyIiIiIitlEREREREduoiIiIiIhtVERERETENn52Bzgft9vNgQMHCA0NxTAMu+OIiIhIPpimSWpqKhUrVsThOP85jyJdRA4cOEBcXJzdMUREROQi7N27l9jY2PPuU6SLSGhoKGAdSFhYmM1pREREJD9SUlKIi4s79T5+PkW6iPx5OSYsLExFRERExMfkZ1iFBquKiIiIbVRERERExDYqIiIiImKbIj1GREREpKhzuVzk5ubaHaPQ+fv743Q6L/n7qIiIiIhcpLS0NPbt24dpmnZHKXSGYRAbG0vp0qUv6fuoiIiIiFwEl8vFvn37CAkJITo6ukQtvGmaJklJSezbt49atWpd0pkRFREREZGLkJubi2maREdHExwcbHecQhcdHc2uXbvIzc29pCKiwaoiIiKXoCSdCfm7gjpuFRERERGxjYqIiIiI2EZFRERERGyjIiIiIiK2KZmzZo5uhWWfgjPgr4dfAPgFQ3AZCCkLIZHWo3QMBITYnVhKENM0OZGRy+GULE5m5JKcmUtyZg7Jmblk5rhxud3kuk3yXG5cbgjwcxDo5yDQ30Ggn5PQQD8iSwUQWTqAsqUCKFs6kNKBJfNXXUTO7o033uCFF1444/Pvvvsujz/+eKFmKZl/O53cbRWR/CoVDRGVIaIKlKkC0XWhXD2IrgN+gd7LKcXa8fQcNh9KZeuRVLYcTmXv8Uz2n8zkwMlMMnJcBfpc4cH+VI4MIS4ymLjIEKpHlaJuTBi1y4cSHHDpKyOKiPWPiMzcgv3dza9gf6dHs1geffRR7r777lMfDx48mDlz5nDLLbd4I955ebWIDB06lKlTp7Jp0yaCg4Np06YNb731FnXq1PHm015YRFVo9wy4ciAvx/qvKxtyMyHjOGQeh4wTkHEUcjMgPcl67E84/fsYTihbEyo0grhWENsCyjcAZ8nsd3JuKVm5JO45yao9J1m55wTrD6RwNC37vF8TWSqAiBB/IoL9iQgJIDzYnyB/J/5OAz+HA3+ngWEY5LrcZOe5yM51k5XnJiUzlxMZORxLy+F4eg6ZuS6SM3NZuz+ZtfuTT3sOhwFVo0pxWUwY8XERNKtahgYVwwnw01VbEU9l5rqoN/hHW557wytdCQnI/3tPaGgooaGhALz00kvMmTOH+fPnExsbS48ePZg/fz6dO3dmypQp3op8ilffMRcsWMDAgQNp0aIFeXl5PP/881x99dVs2LCBUqVKefOpzy+qJnR6MX/7Zp6Ek3v+eOyG4zvgyEY4vB6yTsLRzdZj7WRrf/8QqNQMqrWHmp2gQjw49C/OkiYtO4+lO47x69aj/L79KFuPpHG2FaDjIoOpUz6UWuVDqVo2hEoRIVSMCKJiRDBB/gXzc5OWnce+ExnsOZbB3hOZ7DmWzvakdDYeTOFYeg47ktLZkZTOzLUHAQj0c9A4LoJW1SJpVzuaJnER+DlVTESKo8GDBzN27Fjmz59P1apVAXjssce4++67GTNmTKFkMMxCXCA/KSmJcuXKsWDBAtq1a3fB/VNSUggPDyc5OZmwsLBCSOgB04TUg1Yh2b8S9i6FfcshO+X0/YLLQPUOUPMqqHOtNe5EiqVtR9L4cf0hFmxOYuWeE+S5T//VqhwZQpPKETStXIZGseHULh9KKRvHbpimSVJaNhsPprL+QDIrd58gYfcJTmScfvOu0CA/rqgRRfs60XSuW45yYUE2JRYpWrKysti5cyfVqlUjKCjIpy7NAAwZMoQxY8acVkL+NH/+fIYNG3beMyL/PP6/8+T9u1D/FkxOtk4LR0ae/c04Ozub7Oy/TlenpKScdb8iwTAgrKL1qHWV9Tm32zo7snsRbP8Fdi6EzBOw/lvrYTihaluodyPU7Qah5e09Brkkpmmydn8yP64/xOx1h9ielH7an1eODKFtrSiurBlFi2qRRJUuWuOJDMOgXGgQ5UKDaF87GrCOaXtSOit2Hef37cf4dWsSJzJymb3+ELPXH8IwoHmVMlzboALXNoyhQnjJW9Za5FwMw/Do8oidzldCCluhnRFxu93ceOONnDx5kt9+++2s+/znP//h5ZdfPuPzRfKMSH648mD/Ctg2Dzb/AIfX/u0PDah2JTS+3SomATZeqhKP7D2ewbRV+5m6aj87j/5VPvydBm1qRHFVvfK0qxVN5bK+P9vK5TZZtz+ZBVuS+HnTERL3njztz5tWjqBn01i6Na5IeLC/PSFFbHK+MwJF2WuvvcYHH3zAjBkzTishZcqUITDQ+gdTYZ4RKbQi8tBDD/HDDz/w22+/ERsbe9Z9znZGJC4uzneLyD8d2w4bv7Me+1f89Xn/UlDvJmjSF6pcYZ1tkSIlM8fFd2sOMCVhH8t2Hj/1+WB/Jx3rRtO1fgwd65YjLKh4vxkfOJnJ7HWH+GHdQVbsPnFq3Eugn4Ou9WO4tXksV9SIwuHQz7AUf75YREzTJCIi4qxXHJYtW0aLFi2AYnhp5pFHHuH7779n4cKF5ywhAIGBgafaWLFUtga0fdx6nNwDqydC4jg4sRNWj7ce0XWhxb3Q+DYIDLU7cYm3IymNcUv3MHnFXlKy8gCrJ7apUZaeTWK5pkGMreM8ClvFiGDubluNu9tW43BKFjMSDzA5YS9bDqcxY/UBZqw+QFxkMP0vr0Lv5pUJDynexUzE1xiGcWqYRFHh1TMipmny6KOP8u233zJ//nxq1arl0dcX6cGqBcU0rYGuieNg7TeQ+8ep/oBQiO8DLR+wZvlIoTFNk/mbk/hi0U5+3Xr01OfjIoO5rUVlejSpRMUIjY34059jZSav2Mf0xP2nCluQv4MeTSpxZ5uq1I0ppr+/UqL54hmR/OjSpQurV68mPT2dyMhIJk+eTOvWrc/YzycuzTz88MOMHz+e6dOnn7Z2SHh4OMHBF/6LvEQUkb/LSobVE2DZSDi29Y9PGnBZN+ssSqVmdqYr9nLy3Hy3+gCfLtzB5sOpgHX2o2OdcvRvXYX2taJ1yeECMnNcTE/cz+jfd7HpUOqpz7erHc3DHWrQqlpkib1luhQ/xbWI5JdPFJFz/YUzatQoBgwYcMGvL3FF5E+mCTvmw9LhsGX2X5+veiW0fQJqdNI4kgKUmeNi3NLdfPbrTg6lZAFQOtCPPi3juKN1VeIifX/QaWEzTZNlO48zZvEuZq87xJ8zmZtWjuDhDjXpVLecSp34PBURHygil6rEFpG/O7IRFr1vLZjmtk55E9vSWpCtent7s/m4rFwXXy3ZzfAFO06tchodGsjdV1Tj9laVNQukgOw5lsGIhduZnLCPnDw3AHVjQnnyqtpcVa+8zpCIz1IRUREpWU7uhcUfQcIoyLP+1U7VK6HjC1DlzGt3cm5ZuS4mLNvDx/O3cyTVKiBxkcEM7FCTHk0rEeinlXC94UhqFl/8touvluwmLdsq1Y3jInj66tq0rRmlQiI+R0VERaRkSj0Ev75jFRJXjvW5mlfB1a9CucvszVbEud0mM1Yf4P9+3Mz+k5kAVIoI5l+da9KzaSz+Wsa8UJzMyOHThTsYtWjXqVUoW1WLZNC1dWlSuYzN6UTyT0VERaRkO7kXfn0bVn1lXbIxHND0DusMSelydqcrcn7ffpQ3Zm1k3X5r7nxMWBCPdq7Jrc3idIM3mySlZvPJ/O18tXT3qUs2N8VX5Nlr6lJJs5LEB6iIqIgIWIuk/TTEWiQNIKC0NcOm9SPgr7/MdySl8drMjfy86QhgDUJ9qEMN7mlbrcBuKieX5sDJTN6Zu4VvVu7DNK3F0e67sjoPdahRotZoEd+jIqIiIn+3+3f48QU4sNL6uExVuPa/ULurrbHskpGTx7Cft/HZrzvJcbnxcxjc3qoy/+pcq8jd80Us6/Yn8+r3G1j6x8q15UIDefGGenRrVEHjR6RIUhFREZF/crth/VSY8xKkHrA+V+c6uGaoVUxKANM0+WHdIV77fgMHkq1Bve1rRzO4Wz1qRJe2OZ1ciGmazNlwmDdmbWT3sQzAWsX2lZsaULOcXj8pWlREVETkXLLTYOF/rVk27jzwC4J2T0Obx8AvwO50XrP3eAYvTFvHwi1JAMSWCWbwDfU0RdQHZee5+HTBDob9so3sPDf+ToP721XnkY61CA7QJTUpGlREVETkQo5sgllPw65frY/L1YebPix2K7S63CajFu3kf3O2kJnrIsDPwYPta/BwhxoaB+Lj9hzLYMiMdfyy2SqXVcqG8NbNjbi8elmbk4moiBRUEdF0geKsXF248zvoORJCysKR9fBZF2ssSU6G3ekKxIYDKfT8eBGvzdxIZq6LVtUi+fHxdjx5VW2VkGKgctkQvhjQguH9mhETFsTuYxnc9ukSXvh2LalZuXbHE5ECoCJS3BkGNOoFA5dBw1vBdMPiYfBJa9i1yO50Fy3X5ea9n7Zw47DfWL0vmdAgP97s2ZCv77ucalGl7I4nBcgwDK5pEMOcJ9vRp2VlAMYt3UPXdxfyy+YjNqcT8U1vvPEGhmGc8XjvvfcKPYsuzZQ0W36E75+AlP2AAa0HQqeXwN93TituPZzKk5NWs3a/dSvra+rH8MpN9SkX5jvHIBfv9+1HGfTNWvYct87q9WlZmRevv0xTfaXQ+fKlmdTUVNLT0099PHjwYObMmcNvv/1GbGxsvr6HxojIxctKgR+fh1VjrY/L1YMeI6BCI3tzXYDbbfLFop3898fN5OS5CQ/259XuDbixcUW7o0khy8xx8faczXyxaCemaY0deadXPM2qaGVWKTxnvBGbJuTadNnbP+Sib4b60ksvMXbsWObPn4/T6aR///4cOXIEPz8/XnrpJW699dazfp2KiFy6zT/AjEchPQkc/tDxObjicXAUvbEVB5MzeWJiIkt2WGtMtK8dzX9vaUR5nQUp0X7ffpSnJ63mQHIWDgMe7lCTf3WupdVypVCc8Uackw5v2PQPo+cPQIDnl6UHDx7Ml19+yfz586latSoHDx7k8OHDxMfHc+jQIZo1a8aWLVsoVerM763BqnLp6lwLDy+BujeAOxfmvQJf3gQpB+1Odpof1x/i2vd/ZcmO44QEOHmjR0NG39VCJURoUyOKHx5vR88mlXCbMOyXbdwy/Hd2H0u/8BeLlHBDhgw5rYQAVKhQgfj4eABiYmKIiori+PHjXs2hi6olXako6P0VJI6HWc9YU32HXwHdP7F9VdasXBevz9zI2CW7AWhYKZwP+jTRYFQ5TXiwP+/0jqdLvfI8N3Uta/Ylc/0HvzG0Z0O66bKdFCb/EOvMhF3P7YEhQ4YwZsyY00rIPyUkJOByuYiLiyuAgOemIiLWdcUmfSGuJUy5Cw6thfG94PKB0GUI+BX+kujbjqQycNwqNh9OBeD+dtV5+uo6OuUu53RdwwrEx0Xw2IRVLN91gke/XsWibUcZ0q2+FkGTwmEYF3V5pLC99tprfPLJJ8yYMYOgoCAOHToEQJkyZQgMtP6+P378OHfccQcjR470eh6NEZHT5WXD3MGwdLj1ccWm0GsMRFQutAjTE/fz3NS1ZOS4iCodwP96xdO+dnShPb/4tjyXm/fnbWXYL9swTahVrjSf9GtKzXKhdkeTYsYXZ82YpklERAQpKSln/NmyZcto0aIF2dnZXHXVVdx3333079//nN9LY0TEO/wC4dq3oM8ECC5j3URvRDvY+pPXnzo7z8VL09bx2IREMnJctKlRllmPXakSIh7xczp46uo6jLunFdGhgWw9ksaNwxbx/RqbTpmLFCGGYZCcnIxpmmc8WrRogWmaDBgwgE6dOp23hBQkFRE5uzrXwgMLoWITyDwB426BX94At8srT7fvRAa9hi8+NR7k0U41GXtPK8qF+sa/MqToaVMziln/upLW1cuSkePikfGrePm79eS63HZHEymyFi1axMSJE5k2bRrx8fHEx8ezdu1arz6nLs3I+eVlw+znYMXn1sc1OsHNn0NIZIE9xW9bj/LI1ys5mZFLRIg/7/aKp2PdcgX2/aVky3O5eXvOFoYv2A5A8yplGHZ7U2LCVXLl0vjipZmCpEszUjj8AuGGd6DHp+AXDNt/hpEd4fCGS/7Wpmny2a87uOOLpZzMyKVRbDjfP9pWJUQKlJ/TwaBr6zKifzNCA/1YsfsE3Yb9RsLuE3ZHExFURCS/GveG++ZBRBU4scu6ed6GGRf97bJyXTw5aTWvzdyI24Sbm8Yy6YHWxJbxbAqaSH51rR/Dd4+2pU75UJJSs+nz6RImrdhrdyyREk9FRPKvfH24fz5Uawe56TCp/x/jRjy75n7gZCa3Dl/Mt6v243QYDOlWj7dvbaS75YrXVY0qxdSH29C1fnlyXG6enbKG/8xYT57GjYjYRkVEPBMSCf2+hcsftj5e8JZVSHLyt5Jl4t6T3PTRItbuT6ZMiD9j72nJXVdUw7jIeySIeKpUoB+f9G3GE11qAzD6913c8cUyTmbk2JxMpGRSERHPOf3gmqHW6qvOANj0PYy6FlLOPz3y+zUH6D1iMUmp2dSNCWXGI21pUyOqkEKL/MXhMHisSy2G92tGSICT37cfo8fHv7PzqJaGFylsKiJy8eJvhzu/h5AoOLgaRnaCA4ln7GaaJh/M28oj41eRneemU91yTHmoDXGRGg8i9rqmQQxTH25DpYhgdh5Np8fHi1iy45jdscTHFOHJp15VUMet6bty6U7sgvG9IWmTdb+DniPhshsAa5Gyf09Zw7RE62zJvW2r8dx1l+F06FKMFB1HUrO478sEVu89ib/TYGjPRtzSLNbuWFLEuVwutm7dSkhICNHR0SXqErNpmiQlJZGRkUGtWrVwOk8f4+fJ+7eKiBSMrGSYPMCa3osBXV8nufH93Dd2Bct2HsfPYfBq9wb0aVl4S8WLeCIr18VTk1Yzc6119+lHOtbkqatrl6g3F/FcWloa+/btK5FnRQzDIDY2ltKlS5/xZyoiYg9XHvzw7KnFz74JuJGnU3pRKjCA4f2a0baWxoNI0eZ2m7wzdwvDftkGQM8mlXjz5ka62aKcl8vlIjc31+4Yhc7f3/+MMyF/8uT9W3fflYLj9IPr/8chRzlilg3l5pwZRAQnUenuL6kbpxIiRZ/DYfB01zpULhvCc1PXMnXVfo6kZvNJv6aEBvnbHU+KKKfTec43ZLkw1XwpUAu3HqXzknj+lTOQXPzobC6m7tw7IOO43dFE8q1X8zg+u7M5IQFOftt2lF4jlnA4JcvuWCLFkoqIFJjpifu5e/Ry0nNcHKt+Izl9JkNgGOxZDF9cA8n77I4okm8d65Rjwv2XE1U6gI0HU+j58e9sO5JqdyyRYkdFRArEqEU7eWxCInlukxsbV2TUgJaUqtMJ7v4RQivC0c3weVdI2mJ3VJF8axQbwdSHrqBaVCn2/7Ei8Oq9J+2OJVKsqIjIJTFNk//N2czL31k3wRvQpirv9Y7/a3Bf+XpwzxwoWwtS9sEXXWF/go2JRTxTuWwI3zzUhsax4ZzIyKXPyCX8tvWo3bFEig0VEbloLrfJC9PW8eHP1gyDp6+uzZBu9XD8c42QiDi4ezZUbAKZx2F0N9j+iw2JRS5OZKkAxt13OVfULEtGjou7Ry/nhz+m+YrIpVERkYuSk+fmsQmrGL90Dw4D3ujRkEc61Tr3mgulouDO76Bae+uGeeNuhQ3TCze0yCUoHejHFwNacG2DGHJcbgaOX8nXy/bYHUvE56mIiMeycl08+FUC3685iL/TYNjtTbm9VT4WKgsMhb6Tod5N4M61FkBL/NrreUUKSqCfk2G3N6VPy8q4TXhu6lo++3WH3bFEfJqKiHgkLTuPu0Yt5+dNRwj0czDyjuZc17BC/r+BXyDcMgqa9APTDdMehOWfeS+wSAFzOgze6NGAB9vXAOC1mRt5/6etJXJlTZGCoCIi+ZackUu/z5ayeMcxSgf68eXdLelQp5zn38jhhG4fQqsHrY9nPgWL3i/YsCJeZBgG/76mDk9fXRuAd3/awpuzN6mMiFwEFRHJl+PpOdw2cgmJe08SEeLPuHtb0ap62Yv/hg4HXPMmXPmU9fHcwfDz66C/yMVHGIbBI51q8dIN9QAYsWAHQ2asx+3Wz7CIJ1RE5IKSUrPp8+kSNh5MIap0IBPvb03juIhL/8aGAZ0HQ+ch1scL/ws//UdlRHzKPW2rMbRnQwwDvly8m0FT16iMiHhARUTO60hKFrd9upjNh1MpHxbIxAcup05MaME+yZVPwjVvWduL3oM5L6qMiE/p07Iy7/aKx+kwmLRiH89MWYNLZUQkX1RE5JwOJmfS+9MlbE9Kp0J4EBPvb02N6DNv91wgLn8Qrnvb2l48DGY/pzIiPqV7k0p8cFsTnA6Db1bu4+nJq1VGRPJBRUTOav/JTHqPWMLOo+lUighm0gOtqRpVyrtP2vI+uOE9a3vpJzDraXC7vfucIgXo+kYVGNanCX4Og29X7eeJiYnkufQzLHI+KiJyhgMnM+nz6RL2HM+gcmQIEx+4nLjIkMJ58uZ3wY3DAMOa1jvrKZ0ZEZ9ybcMKfNS3Kf5OgxmrD/CYyojIeamIyGkOJWfRZ+RfJWTC/ZcTW6aQSsifmvaH7p8ABqz4AmY9ozIiPqVr/Rg+6dsMf6fBzDUHeWLSapURkXNQEZFTDqdYJWT3sQziIoP5+v7LqRgRbE+Y+D5w00dYZ0ZGwuxBKiPiU7rUK8/wflYZ+W71AY0ZETkHFREBrNkxfT61xoTElgnm6/sup5JdJeRPTfrCjR9a20uHw48vqIyIT+l8WXmG3d4UP4fBtMQDPDtFU3tF/klFRDialk2fkUvY8cfA1K/vs+FyzLk07Q/d/lh1dclH1sJnKiPiQ7rWj+HDPn/Npnlu6lqVEZG/UREp4U6k59Dvs6VsT0qnYngQE+4vxIGp+dVsAFz/jrX9+wcwf6itcUQ8dW3DCrzXOx6HARNX7OWl6eu0HLzIH1RESrDkzFzu+GIZmw6lUi40kPH3FcES8qcW98C1/7W2F7wFv75jbx4RD3VrXJF3e8djGDBu6R5em7lRZUQEFZESKz07j7tGLWPt/mQiSwUw7t5W3l8n5FK1egC6vGxtz3sZlgy3N4+Ih26Kr8RbPRsB8PlvO3ln7habE4nYT0WkBMrMcXHPmOWs3HOS8GB/vrqnFbXKF/Cy7d7S9nFoP8janv1vSBhtZxoRj/VqEcfLN9YH4MOft/Hx/G02JxKxl4pICZOT5+ahcQks2XGc0oF+fHl3S+pVDLM7lmc6DII2j1rb3z0OaybZGkfEU3e2qcq/r6kLwH9nb2bUop02JxKxj4pICeJymzwxMZH5m5MI8ncw6q4WBXMX3cJmGHDVq9DiPsCEbx+ETbPsTiXikYc61OBfnWoC8PJ3G5i0Yq/NiUTsoSJSQpimyfNT1zJz7UH8nQYj+jenRdVIu2NdPMOwBq827gOmCyYPgJ0L7U4l4pEnrqrNvW2rATDomzXMXnfQ5kQihU9FpAQwTZPXZ25k4oq9OAz44LYmtK8dbXesS+dwWPelqXsDuLLh6z6wL8HuVCL5ZhgGL1x/Gb2ax+I24V9fJ/Lb1qN2xxIpVF4tIgsXLqRbt25UrFgRwzCYNm2aN59OzuHDn7fx2W/WNeg3b27EtQ0r2JyoADn94ObPoVo7yEmDcTfD4Q12pxLJN8MwGNqzEdc2iCHH5eb+sStYueeE3bFECo1Xi0h6ejqNGzfmo48+8ubTyHmMXbzr1BTBwTfUo1fzOJsTeYF/ENw2Hio1g8wTMLYHnNhldyqRfHM6DN67LZ4ra0WRkePirlHL2XQoxe5YIoXCq0Xk2muv5bXXXqNHjx7efBo5h+/XHGDwjPUA/KtzLe7+41p0sRQYCn2nQLl6kHbIKiNpR+xOJZJvgX5ORvRvRtPKEdZig58vY+/xDLtjiXhdkRojkp2dTUpKymkPuTgLtyTxxMRETBP6X16FJ7rUsjuS94VEQr+pEF4Zju+Ar26GLP0Mie8ICfBj1ICW1I0J5UhqNnd8sYyjadl2xxLxqiJVRIYOHUp4ePipR1xcMbyMUAgS957kwa8SyHWZ3NCoAv+5sT6GYdgdq3CEVYD+30JIFBxaAxNuh9wsu1OJ5Ft4iD9j7m5JbJlgdh5N565Ry0nLzrM7lojXFKki8txzz5GcnHzqsXev5tV7atuRNO4atYyMHBdX1orinV7xOB0lpIT8Kaom9JsCAaVh168w9V5wu+xOJZJv5cOC+PLulkSWCmDt/mQeGLuC7Dz9DEvxVKSKSGBgIGFhYac9JP8OJWdx5xfLOJGRS+O4CIb3a0aAX5F6iQtPxSbWAFZnAGz8DmY+CbrBmPiQ6tGlGX1XC0oFOFm07RhPTlyNy62fYSl+Sui7VPGTnJnLgFHL2H8yk+pRpRg1oAWlAv3sjmWv6u2h50jAsO5Js+C/dicS8Uij2AhG9G+Ov9Ng5tqDvPLdet2xV4odrxaRtLQ0EhMTSUxMBGDnzp0kJiayZ88ebz5tiZOV6+L+L1ew6VAq0aGBjPnjlK4A9bvDdf9nbc9/AxLG2BpHxFNta0Xxbu94DAPGLN7N8AU77I4kUqC8WkRWrFhBkyZNaNKkCQBPPvkkTZo0YfDgwd582hLF5TZ5clIiS3ceJzTQjzF3tSQuMsTuWEVLy/vgyqes7e+fgM2z7c0j4qEbGlXkpevrAfDW7E1MXbnP5kQiBccwi/B5vpSUFMLDw0lOTtZ4kbMwTZP/zFjPmMW7CXA6GH13C9rUiLI7VtFkmjB9ICSOA79guPM7iGthdyoRj7wxayOfLtyBn8PgiwEtaFccbtUgxZIn798aI+LDRizcwZjFuwH4X6/GKiHnYxjQ7X2oeRXkZcL4XnB0m92pRDwy6Jq63BRfkTy3yUNfJbBuf7LdkUQumYqIj5qeuJ83f9gEwEs31KNb44o2J/IBTn+4dTRUbAqZx6370qQl2Z1KJN8cDoP/u6UxV9QsS3qOiwGjlmv1VfF5KiI+aPH2Yzw9eTUAd19RjXuK89LtBS2wNNw+CcpUte5HM74X5KTbnUok3wL8HAzv14zLKoRxNC2bO0ct42RGjt2xRC6aioiP2XI4lfvHriDXZXJdwxhevP4yuyP5ntLR0PcbCI6EAythyj3g0sqV4jtCg/wZNaAFFcKD2JGUzn1friArVwueiW9SEfEhh1OyGPDFMlKz8mhepQzv9IrHUdJWTS0oUTWhzwTwC4ItP8APz2rBM/EpMeFBjLqrBaGBfizfdYKnJ6/GrQXPxAepiPiItOw87hq1nAPJWVSPLsXIO5oT5O+0O5Zvq9zqrwXPVnwOi963O5GIR+rGhDGifzP8nQbfrznIW7M32R1JxGMqIj4gz+XmkfEr2XAwhajSAYy5qyVltGBZwah3I1wz1Nr+aQism2pvHhEPtakZxVs3NwKsmXRjF++yN5CIh1REijjTNBkyYz3zNycR5O/gsztbaMGygnb5Q9DqIWv72wdhz1J784h4qGfTWJ6+ujYAQ2as5+dNh21OJJJ/KiJF3MhfdzBu6R4MA97r3YT4uAi7IxVPXV+HOteDKxsm9IFj2+1OJOKRgR1r0qt5LG4THhm/ivUHtMaI+AYVkSJs1tqDvDHLuub7wnWXcU2DGJsTFWMOJ9w80rprb8YxGHcrZBy3O5VIvhmGwes9GnJFzbJk5Li4e/RyDiZn2h1L5IJURIqoVXtO8MTERADuaF1Fa4UUhoBS0GcihMfB8e0woS/kZdudSiTf/J0OPu7bjFrlSnM4JZu7R68gLVtT06VoUxEpgvYez+C+L1eQneemU91yDL6hHoahabqFIrQ89J0MgWGw53eY/oim9YpPCQ/254sBLYgqHcjGgyk8Mn4leS633bFEzklFpIhJycrlnjHLOZqWQ70KYXzYpwl+Tr1MharcZdDrSzCcsHYSLPiv3YlEPBIXGcLndzYnyN/B/M1JvPr9BrsjiZyT3uGKkDyXm4HjVrLlcBrlQgP5fEBzSgX62R2rZKrREW54x9qe/wasnWJvHhEPNY6L4L3e8QCMWbyb0Yt22htI5BxURIqIP6fp/rr1KMH+Tj6/swUVwoPtjlWyNRsArR+xtqc9rGm94nOuaVCBQdfWBeCV7zfwy6YjNicSOZOKSBHxxaJdp6bpvn9bPA1jw+2OJABXvfK3ab23WzfKE/EhD7Sr/rdpvSvZeDDF7kgip1ERKQLmbTzMazOta7gvXHcZV9fXNN0i489pvTGNIOMojOsFWVqfQXyHYRi81r0hrauXJT3HxT2jl3MkNcvuWCKnqIjYbOPBFP719SpME/q0rKxpukVRQCm4fSKEVoSjm2HyXbpbr/iUAD8Hw/s1o3pUKQ4kZ3Hflwm6W68UGSoiNkpKzebeMStIz3HRpkZZXrmpvqbpFlVhFaHP1+AfAtvnwY/P251IxCPhIda03ogQf1bvPckzU9Zgamq6FAEqIjbJynXxwNgV7D+ZSbWoUnzctyn+mqZbtFWMhx4jrO1lI2D5Z7bGEfFU1ahSfNK3GX4Og+9WH+CDedvsjiSiImIH0zT59zdrWLnnJOHB/nx+Z3MiQnQ3XZ9Q70boPNjanvUsbP/Z3jwiHmpdoyyv92gAwLs/beH7NQdsTiQlnYqIDYb9vI3piQfwcxh80rcp1aNL2x1JPNH2SWjcB0wXTBoASVvsTiTikd4tKnPvH+PRnpq0mtV7T9obSEo0FZFCNnvdQf4313rjeuWmBrSpGWVzIvGYYUC39yHucshOhq976wZ54nOeu+4yOtctR3aem/u+XKEb5IltVEQK0br9yTwxcTUAA9pU5fZWlW1OJBfNLxBuGwfhleH4Dph8J7hy7U4lkm9Oh8H7fZpQNyaUI6nZ3P9lApk5mkkjhU9FpJAcSc3ivi9XkJnrol3taF68/jK7I8mlKhUFt0+AgNKwcyHMHmR3IhGPlA70Y+QdzYksFcDa/ck8PXm1ZtJIoVMRKQRZuS7u/zKBg8lZVI8upRvZFSfl60PPkYBhzaJZNtLuRCIeiYsMYXi/Zvg7DWauPaiZNFLo9G7oZaZp8tzUtSTu/XOGTAvCg/3tjiUFqe510GWItf3Dv2H7L/bmEfFQy2qRvNb9r5k0M9cctDmRlCQqIl42YuEOvl21H6fD4OO+TakWVcruSOINVzwOjW6zZtJMvhOObbc7kYhHerf4a2XnpyYnsnafbmUghUNFxIvmbTzMW7M3AfCfbvW4QjNkiq8/Z9LEtrDuRfP1bbonjfic566tS/va0WTlurl/7Ardk0YKhYqIl2w9nMpjExIxTbi9VWX6XV7F7kjibf5B0HschFWCo1vgm3vBrVkI4jv8nA4+vL0JNaJLcTA5iwfHJpCdp59h8S4VES84kZ7DvV+uIC07j1bVIvlPN91DpsQILW9N6/ULgq1zYN7LdicS8UhYkD8j72hOWJAfK/ec5IVv12kmjXiVikgBy3W5GTh+JbuPZRBbJpiP+zYlwE//m0uUik2g+8fW9qL3YfVEe/OIeKh6dGmG3d4UhwFTEvbx+W877Y4kxZjeIQvY6zM38vv2Y4QEOBl5R3PKlg60O5LYocHNcOXT1vaMR2HfCnvziHioXe1oXri+HgBvzNrIgi1JNieS4kpFpABNXL6H0b/vAuDd3vFcViHM3kBir44vQJ3rwZUNE/pCiqZEim+5+4qq9Goei9uER8avZHtSmt2RpBhSESkgK3Yd58Vp6wB48qradK0fY3MisZ3DAT1HQPRlkHYIJvaFXM1CEN9hGAavdm9AsyplSM3K474vV5CSpVsZSMFSESkAB05m8uBXCeS6TK5rGMOjnWraHUmKisBQ6DMegsvA/gT4/nHQwD/xIYF+Tob3a0aF8CB2JKXz2NercLn1MywFR0XkEmXmuLh/7AqOpuVwWYUw3r61sWbIyOkiq8Oto8FwwuqvYfFHdicS8Uh0aCCf9m9OoJ+DXzYn8faczXZHkmJEReQSmKbJoKlrWLc/hchSAYy8oxkhAX52x5KiqHoH6PqGtT33Jdj2k61xRDzVMDac/97SCIBP5m9neuJ+mxNJcaEicglGLNzB9MQD+DkMPunblNgyIXZHkqKs1QMQ3w9MN0y5W8vAi8+5Kb4SD7avAcCzU9ZoGXgpECoiF2n+5iOnlm8fcmN9WlUva3MiKfIMA254B2Jb/rEMfB/ISrE7lYhHnulah451osnOs5aBT0rNtjuS+DgVkYuwIymNR79ehWlCn5aV6deqst2RxFf4BULvsRBaAY5uhm8fALfb7lQi+eZ0GLzfpwnV/1gG/uFxCeTk6WdYLp6KiIdSs3K578sVpGbl0bxKGV6+Ucu3i4dCY6x70jgDYfMsWPCm3YlEPPLnMvChgX4s33WCl79bb3ck8WEqIh5wu00en5DI9qR0YsKC+Liflm+XixTbzLpbL8CCt2DDDHvziHioRnRp3u8Tj2HAuKV7GL90j92RxEfpXdQD78zdwrxNRwj0c/DpHc0oFxpkdyTxZfF94PKB1va3D8Jh/atSfEunuuV5+uo6AAyZsY4Vu47bnEh8kYpIPs1ae5Bhv2wDYGjPhjSKjbA3kBQPV71iTe3NTbcGr2boL3LxLQ93qMH1DSuQ6zJ58KsEDpzMtDuS+BgVkXzYdCiFpyevBuDettXo2TTW5kRSbDj94JZRUKYqnNxtTet15dmdSiTfDMPg/25txGUVwjialsODXyWQleuyO5b4EBWRCziZkcP9XyaQkeOibc0oBl1b1+5IUtyERMJt48E/BHb8Aj8NsTuRiEdCAvz4tH8zyoT4s2ZfMs9/uxZTtzKQfFIROY88l5tHv17FnuMZxEUG82GfJvg59b9MvKB8fej+ibW9eBismWRvHhEPxUWGMOz2pjgdBlNX7mfUol12RxIfoXfV83hr9iZ+3XqUYH8nn/ZvTplSAXZHkuKsfne48ilre8ajcGCVrXFEPHVFzSiev+4yAF6ftZHftx21OZH4AhWRc5ieuJ+Rv+4E4H+9GnNZhTCbE0mJ0PEFqNUV8rJgQj9IS7I7kYhH7r6iKj2bVsLlNhk4fiV7j2fYHUmKOBWRs1i3P5lnp6wBrBHh1zWsYHMiKTEcTrh5JJStCSn7YPKd4Mq1O5VIvhmGwRs9GtIoNpwTGbncPzaBzBwNXpVzUxH5h2Np2TwwNoHsPDcd60Tz1B9z5EUKTVC4NXg1IBR2L4Ifn7c7kYhHgvydDO/XjKjSAWw8mMKz36zR4FU5JxWRv8l1uRk4fiX7T2ZSLaoU793WBKdDy7eLDaLrQM9Pre1ln8Kqr+zNI+KhihHBfNy3GX4Og+9WH+DThTvsjiRFlIrI37w+cyNLdhyndKAfI+9oRniwv92RpCSrex10eM7a/v4J2LfC3jwiHmpZLZIhN9YHrMH/C7dozJOcSUXkD1MS9jH6910AvNOrMTXLhdobSASg3bNQ53pw5cDEfpB62O5EIh7p16oyvZvH4Tbh0a9XsftYut2RpIhREQFW7z3J89+uBeBfnWtxdf0YmxOJ/MHhgB7DIaoOpB6ESXdAXo7dqUTyzTAMXulenyaVI0jOzOX+LxNIz9bqwfKXEl9EklKzefCrBHLy3HS5rDyPd65ldySR0wWFWYNXA8Nh7xKY/W+7E4l4JNDPGrwaHRrI5sOpPDNltQavyimFUkQ++ugjqlatSlBQEK1atWLZsmWF8bQXlJPnZuC4lRxMzqJGdCne7d0YhwanSlEUVdOa1osBK76AhNF2JxLxSPmwIIb3a4a/02DW2kN8PH+73ZGkiPB6EZk4cSJPPvkkQ4YMYeXKlTRu3JiuXbty5MgRbz/1Bb02cwPLdh0nNNCPT+9oTmiQBqdKEVa7K3R6wdqe+TTsLRqFXiS/mlUpwys3NQDg7Tmb+WWT/e8DYj+vF5F33nmH++67j7vuuot69eoxfPhwQkJC+OKLL7z91Oc1aflevly8G4B3e8dTI7q0rXlE8uXKp+GyG8GdCxP7Q8pBuxOJeKRPy8rc3qoypgn/mrCKnUc1eLWk82oRycnJISEhgS5duvz1hA4HXbp0YfHixWfsn52dTUpKymkPb1i15wQvTlsHwJNX1aZLvfJeeR6RAmcY1s3xoi+DtEMwqT/kZdudSsQj/+lWn+ZVypCalcd9X64gTYNXSzSvFpGjR4/icrkoX/70N/ry5ctz6NChM/YfOnQo4eHhpx5xcXFeyXUyM5cAPwdX1yvPIx1reuU5RLwmsDTcNs5agXXfcpj1NGjgn/iQAD8HH/drSvmwQLYdSePJiYm43foZLqmK1KyZ5557juTk5FOPvXv3euV5OtYpx/RHruCd3vEanCq+qWwNuPkLwICVX1oDWEV8SLlQa/BqgNPBnA2HGfbLNrsjiU28WkSioqJwOp0cPnz6IkyHDx8mJubMtToCAwMJCws77eEtNaJLUzrQz2vfX8TranWBzoOt7R/+DXuW2JtHxENNKpfhte7W4NV3f9rCvI1asK8k8moRCQgIoFmzZsybN+/U59xuN/PmzaN169befGqRkqHtE1Cv+98Grx6wO5GIR3q1iKP/5VUwTXh8QiLbk9LsjiSFzOuXZp588klGjhzJmDFj2LhxIw899BDp6encdddd3n5qkeLPMKD7x1CuPqQfsZaBz82yO5WIR166oR4tqpYhNTuP+79cQWpWrt2RpBB5vYj07t2bt99+m8GDBxMfH09iYiKzZ88+YwCriFykgFJ/DF6NgP0JMPMpDV4VnxLg5+Djvs2ICQtie1I6T0xcrcGrJYhhFuF1dlNSUggPDyc5Odmr40VEioXtP8NXN4Pphuvehpb32Z1IxCOJe0/Sa8RicvLcPN6lFo93qW13JLlInrx/F6lZMyJyCWp0gi4vW9uzB8Gu3+zNI+Kh+LiIU4NX3/tpK3M3aPBqSaAiIlKctHkUGtwC7jyYdCec9M4UeBFv6dU8jjtaVwHgiYmJbDuiwavFnYqISHFiGHDjhxDTEDKOwsS+kJtpdyoRj7x0Qz1aVo0kLTuPB8Zq8GpxpyIiUtwEhMBt4yGkLBxcDd89psGr4lP8nQ4+6tuUCuEavFoSqIiIFEcRleHW0WA4Yc1EWPyR3YlEPBIdGmitvOrn4KeNh/ng5612RxIvURERKa6qtYOur1vbc1+C7b/Ym0fEQ43jInhdg1eLPRURkeKs1YPQ+HZrSu+Uu+D4TrsTiXjk1r8NXn1yolZeLY5URESKM8OAG96Fik0h8wRM6As56XanEvHIn4NXtfJq8aQiIlLc+QdB76+gVDk4sh6mPazBq+JT/hy8+ufKq09N0uDV4kRFRKQkCK9klRGHP2yYBr/+z+5EIh6JDg1keP9mBDgdzNlwmI/nb7M7khQQFRGRkqJyK7j+bWv759dgy4/25hHx0N9XXv3f3C38sumIzYmkIKiIiJQkzQZA83sAE765F5K22J1IxCO9WsTR7/LKmCb8a8Iqdh3VmCdfpyIiUtJc8yZUbgPZKTChD2SetDuRiEcG31Cf5lXKkJqVxwNjE0jPzrM7klwCFRGRksYvAHp9CWGxcGwbTL0P3C67U4nkW4Cfg4/7NiU6NJDNh1N59ps1FOEbycsFqIiIlESlo+G2ceAXBFvnwM+v2p1IxCPlwoL4pG9T/BwGM9cc5NOFO+yOJBdJRUSkpKoYDzf9sfT7b+/C2im2xhHxVPOqkQzpVg+At2Zv4retR21OJBdDRUSkJGt4C1zxmLU9/RE4kGhrHBFP9bu8Crc0i8VtwqNfr2Tv8Qy7I4mHVERESrrOQ6DmVZCXCRNuhzRNiRTfYRgGr3VvQMNK4ZzIyOWhcQlk5WrMky9REREp6RxOuPkzKFsLUvbDxP6Ql213KpF8C/J3Mrx/MyJLBbBufwovTlunwas+REVERCA4Avp8DYHhsHcJzHpay8CLT6kUEcyHfZrgMGBKwj7GLd1jdyTJJxUREbFE1YJbPgfDASu/hGUj7U4k4pErakbx7DV1AXj5u/Ws3HPC5kSSHyoiIvKXWldBl5et7dmDYMcCe/OIeOiBdtW5tkEMuS6Th75KIClVlxmLOhURETldm0eh0W1gumDynXBc6zOI7zAMg/+7tTE1y5XmcEo2A8evJNfltjuWnIeKiIiczjCg2/tQqTlknoDxt0FWit2pRPKtdKAfI/o3o3SgH8t2Hue/szfZHUnOQ0VERM7kH2StvBpaEY5utm6Qp2XgxYfUiC7N27c2BmDkrzuZueagzYnkXFREROTsQmP+tgz8jzDvFbsTiXjkmgYxPNi+BgDPTFnNtiOpNieSs1EREZFzq9T0r2XgF70HqyfaGkfEU09fXZs2NcqSkePi/rEJpGbl2h1J/kFFRETOr+EtcOVT1vaMR2HvcnvziHjAz+nggz5NiAkLYkdSOs9O0Z16ixoVERG5sI4vQt0bwJVtLQN/cq/diUTyLap0IB/3a4q/0+CHdYd0p94iRkVERC7M4YAeI6B8Q0g/AhP6QE663alE8q1p5TIM7lYfgP/+uJklO47ZnEj+pCIiIvkTWBr6jIdS0XBoLUy9H9xan0F8R79WlenRpBIut8kj41dxJCXL7kiCioiIeCKiMvQeB84A2PQ9/PK63YlE8s0wDF7v0YA65UM5mqbFzooKFRER8UzlVtDtA2v717dhzSR784h4ICTAj0/6NaV0oB/Ld53QYmdFgIqIiHguvg9c8Zi1PX0g7Flqbx4RD1SPLs3btzYCrMXOZq3VYmd2UhERkYvT+T9/zKTJsWbSnNhtdyKRfLumQQXub1cdgGenrGFHUprNiUouFRERuTh/zqSJaQgZR+Fr3ZNGfMuzXevQslokadl5PDxuJZk5uo2BHVREROTiBZaGPhOhdAwc2QDf3KN70ojP8HM6GNanCVGlA9l0KJXB09fZHalEUhERkUsTXgn6fA1+wbB1Dvz4gt2JRPKtXFgQH/SJx2HA5IR9TFquxfoKm4qIiFy6Sk2hx3Bre+knsGykvXlEPNCmRhRPXV0HgJemr2P9gWSbE5UsKiIiUjDqd4fOg63tH56FrT/ZGkfEEw+1r0HHOtFk57l5eNxKUnRzvEKjIiIiBaftkxDfF0w3TB4Ah9fbnUgkXxwOg3d6xVMpIpjdxzJ4drJujldYVEREpOAYBtzwHlS9EnJSYVwvSD1sdyqRfClTKoCP+lo3x5u9/hCjf99ld6QSQUVERAqWXwD0Hgtla0HKPvi6N+Rk2J1KJF/i4yJ47trLAHhj1kYS9560N1AJoCIiIgUvuAz0nQTBkXBgFUy9T9N6xWfcdUVVrqkfQ67LZOC4lSRnaLyIN6mIiIh3RFa3pvU6A60b5M15ye5EIvliGAb/vbURlSND2H8yk6cmr9Z4ES9SERER76l8OfT4xNpe8hEsHWFvHpF8Cgvy5+O+TQlwOvhp42E++3Wn3ZGKLRUREfGuBjdD5yHW9uxBsGmWvXlE8qlBpXBe6lYPgLdmb2LlnhM2JyqeVERExPvaPgFN77Sm9X5zD+xfaXcikXzp16oyNzSqQJ7b5NHxqzRexAtURETE+wwDrv8f1OgMuRkwvjec2GV3KpELMgyDoT0bUqWsNV7kmSkaL1LQVEREpHA4/eHW0VC+AaQfga9ugYzjdqcSuaDQIH+G9bHGi8zZcFjrixQwFRERKTxBYdB3MoRVgmNb4es+kJtldyqRC2oYG87z19UFrPVF1uw7aW+gYkRFREQKV1hF6DsFAsNh75I/1hhx251K5ILubFOVrvXLk+syeWT8Kt2PpoCoiIhI4StfD24bB84A2DgD5rxgdyKRCzIMg//e3JhKEcHsOZ7B81PXarxIAVARERF7VLsSuv+5xsjH8PuH9uYRyYfwEH+G3d4EP4fB92sOMmnFXrsj+TwVERGxT8Nb4KpXrO05L8KayfbmEcmHJpXL8NTVdQAYMmM9Ww+n2pzIt6mIiIi92vwLLn/Y2p72EGz/2d48IvnwQLvqXFkriqxcN49+vYqsXN1L6WKpiIiIvQwDrn7dWoHVnQsT+1s3yhMpwhwOg//1akxU6QA2HUrl9Zkb7Y7ks7xWRF5//XXatGlDSEgIERER3noaESkOHA5rvEi19pCTBuNuheM77E4lcl7lQoN4p1c8AGOX7Gb2uoP2BvJRXisiOTk53HrrrTz00EPeegoRKU78AqH3VxDTCNKTYGwPSD1sdyqR82pXO5oH2lcH4Nkpa9h/MtPmRL7Ha0Xk5Zdf5oknnqBhw4beegoRKW6Cwqw1RspUtZaA/+pmyDxpcyiR83v66jo0josgJSuPJyYm4nJrSq8nitQYkezsbFJSUk57iEgJE1oe+n8LpcrB4bV/rL6qf2VK0eXvdPDBbfGUCnCybOdxPv5lm92RfEqRKiJDhw4lPDz81CMuLs7uSCJih8jq0H+qtfrqnt9h8l3gyrM7lcg5VSlbile7NwDgvXlbSdh9wuZEvsOjIjJo0CAMwzjvY9OmTRcd5rnnniM5OfnUY+9eLRQjUmLFNITbJ4BfEGz5AWY8qqXgpUjr2TSW7vEVcblNHpugJeDzy8+TnZ966ikGDBhw3n2qV69+0WECAwMJDAy86K8XkWKmShvrjr0T+sLq8RAcAV3fsKb8ihRBr3RvQMKeE+w9nsngaet477Ymdkcq8jwqItHR0URHR3sri4jImepcCzcNsxY7W/IxBIVDh0F2pxI5q7Agf97r3YReIxYzLfEA7WpH07NprN2xijSvjRHZs2cPiYmJ7NmzB5fLRWJiIomJiaSlpXnrKUWkuIq/Ha55y9qePxQWf2xvHpHzaFalDI93rgXA4Onr2Xs8w+ZERZvXisjgwYNp0qQJQ4YMIS0tjSZNmtCkSRNWrFjhracUkeLs8geh44vW9o/Pwcqx9uYROY+HO9akRdUypGXn8fjERPJcGt90Ll4rIqNHj8Y0zTMeHTp08NZTikhx1+5paPOotf3dv2D9t/bmETkHp8PgnV7xhAb6kbD7BB/P3253pCKrSE3fFRE5L8OAq16FZgPAdMM398Lm2XanEjmruMgQXuleH4D3521l1R5N6T0bFRER8S2GAde/Aw1uAXceTOoP2+bZnUrkrLrHV6JbY2tK7xMTE0nP1no4/6QiIiK+x+GEHsOh7g3gyrGm9+76ze5UImcwDIPXbmpAxfAgdh3L4JXvNtgdqchRERER3+T0h1tGQa2rIS8TxvWCPUvtTiVyhvAQf/7XKx7DgIkr9vLj+kN2RypSVERExHf5BUCvsVC9A+Smw7hbYP9Ku1OJnKF1jbLc385a8PO5qWtJSs22OVHRoSIiIr7NPwhuGw+V20B2CoztAQcS7U4lcoYnr6pN3ZhQjqfn8NzUNZim7tILKiIiUhwElIK+kyC2JWSdhC9vgoOr7U4lcppAPyfv9o4nwOngp41HmLRC91MDFRERKS4CQ6HfNxDb4m9lZI3dqUROc1mFMJ66ujYAr3y3gT3HtOqqioiIFB9BYVYZqdQcMk/AlzfCobV2pxI5zb1XVqdl1UjSc1w8OSkRl7tkX6JRERGR4iUoHPpP/auMjLlRZ0akSHE6DP7XqzGlA/1YsfsEny7cYXckW6mIiEjxc6qMNIPM4zCmm2bTSJESFxnC4G71AHhn7mY2HUqxOZF9VEREpHgKCof+354+gHXvcrtTiZxya7NYulxWnlyXyVOTVpOTVzJvjKciIiLF159nRk5N7e0OuxfbnUoEsFZdfaNnAyJC/Fl/IIVhv2yzO5ItVEREpHgLDIV+U6BaO8hJg696ws6FdqcSAaBcaBCvdW8AwEe/bGPNvpP2BrKBioiIFH8BpeD2SVCjE+RmwLhbYcscu1OJAHBDo4rc0KgCLrfJk5NWk5XrsjtSoVIREZGSwT8Ybvsa6lwHeVkwoQ+sm2p3KhEAXr2pAVGlA9l2JI13526xO06hUhERkZLDPwh6fQkNbwV3HnxzD6z80u5UIpQpFcDQng0B+PTXHazYddzmRIVHRUREShanP/QYAc0GgOmGGY/C4o/tTiXCVfXKc0uzWEwTnpmypsRcolEREZGSx+GEG96DNv+yPv7xOfj5ddBNyMRmL91Qj/Jhgew8ms7/5my2O06hUBERkZLJMOCqV6DTi9bHC/8L3z8O7pLxr1ApmsKD/U9dovn8t50k7D5hcyLvUxERkZLLMKDdM3D9O4ABCaNh8p2Qm2V3MinBOtUtT8+mlXCb8MyU4j+LRkVERKTFPdBrDDgDYON3MO4WyEq2O5WUYENuqE+50EB2JKXz7k/FexaNioiICEC9m6w79waEwq5fYdT1kHLQ7lRSQoWH+PNGD+sSzciFO1i1p/heolERERH5U7V2cNdMKFUODq+Fz7rA4Q12p5ISqku98vRo8uclmjVk5xXPSzQqIiIif1ehMdw7F8rWgpR98EVX2LHA7lRSQg3pVu/UQmcf/Vw870WjIiIi8k9lqsI9c6Bya+tmeV/dDKsn2J1KSqCIkABevak+AB/P387Ggyk2Jyp4KiIiImcTEgn9p0H9nuDOhW8fgF+Gaq0RKXTXNqzANfVjyHOb/PubNeS53HZHKlAqIiIi5+IfBDd/Dlc8Zn284E2YcjfkZtqbS0qcV26qT1iQH2v2JfPFop12xylQKiIiIufjcFgLn934ITj8YP1UGHUdpB6yO5mUIOXCgnjxhnoA/G/OFnYdTbc5UcFRERERyY+md1iXaoLLwIGVMLITHFxtdyopQW5tFkvbmlFk57n59zdrcLuLx2VCFRERkfyqdiXcOw+iakPKfvjiGlg31e5UUkIYhsHQng0J9neydOdxJizfa3ekAqEiIiLiibI14J65UKMT5GbAlLtg7mDdo0YKRVxkCE93rQPA0B82ciTF929HoCIiIuKp4Ai4ffJfd+9d9L61LHzGcVtjSckwoE1VGsWGk5qVx8vf+f6CeyoiIiIXw+kHV79qzarxC4btP8PIjnBond3JpJhzOqxLNE6Hwcy1B5m38bDdkS6JioiIyKVoeIu1EmtEZTixCz7rDKu+sjuVFHP1K4Zz75XVAHhp2jrSs/NsTnTxVERERC5VTEO4fwHU7AJ5WTB9IEwbCDkZdieTYuzxzrWJiwzmQHIW/5vju3foVRERESkIIZHWuJFOL4LhgMSvrJvmHS2e9wcR+wUHOHmtu3WH3tG/72T13pP2BrpIKiIiIgXF4YB2z8Ad0607+B5ZD5+2131qxGva146me3xF3CYMmrrWJ5d/VxERESlo1drBg79ClbaQk2bdp+ab+yCr+N2wTOz34g31iAjxZ+PBFEb/vsvuOB5TERER8YbQGLhzBnR8EQwnrJ0Ew9vC3uV2J5NiJqp0IM9dWxeAd+Zu4cBJ37oXkoqIiIi3OJzQ/hm4e7Y1q+bkbviiK8x/C1y5dqeTYuTWZnE0r1KGjBwXr/jY2iIqIiIi3hbXEh78DRrcAqYL5r9hDWQ9stHuZFJMOBwGr/VogJ/DYPb6Qz61toiKiIhIYQgKh5s/g56fQVAEHEyEEe2sVVm1PLwUgLoxYdzT1lpbZPD09WTm+MbPlYqIiEhhMQxodCs8vARqXQ2uHOs+NaOuhaTNdqeTYuCxLrWoFBHM/pOZfPDzVrvj5IuKiIhIYQurALdPghuHQUAo7F1qDWSd/ybkZdudTnxYSIAf/7mxPgAjF+5gy+FUmxNdmIqIiIgdDAOa9oeHF/91dmT+UBh+JexebHc68WFX1SvPVfXKk+c2eXHaOkzTtDvSeamIiIjYKSLOOjtyyxdQKhqOboZR18CMf0H6MbvTiY8a0q0eQf4Olu08zvTEA3bHOS8VERERuxkGNLgZBi6DJv2tz60cAx82hWUjNZhVPBZbJoRHO9UC4PVZG0nJKrrTxVVERESKipBIuGkY3PUDlG8AWSdh1tMwoj3s/t3udOJj7r2yGtWiSpGUms17c4vuwFUVERGRoqZKG+tuvte9bU37PbzWmlkz6Q44tt3udOIjAv2cpwaujlm8i02HiuYtBlRERESKIqcftLwPHl0JTe8EDNgwHT5qCbOegfSjdicUH9C+djTXNojB5TYZPG19kRy4qiIiIlKUlYqCGz+AhxZBzavAnQfLPoX342HBf3UjPbmgF2+oR7C/k2W7jjMtcb/dcc6gIiIi4gvK14d+U+COGVChMeSkwi+vw3sNYeH/qZDIOVWKCObRzjUBeH3mpiI3cFVFRETEl1RvD/fNh5s/h6ja1oDWn1+D9xvBwrchK9nuhFIE3du2OtWjSnE0LZsP5xWtgasqIiIivsbhgIa3WEvF9/wMytaCzBPw86vwTn348QVI3md3SilCAvwcDO5WD4BRi3ax7UiazYn+oiIiIuKrHE7r3jUDl1qFJLqudclm8TB4vzFMvR8OrrY7pRQRHeqUo3PdcuS5TV75fkORGbiqIiIi4uv+LCQPL4HbJ0PVK61BrWsmWnf4/awLJH4NuVl2JxWbvXRDPQKcDhZuSWLexiN2xwG8WER27drFPffcQ7Vq1QgODqZGjRoMGTKEnJwcbz2liEjJZhhQ+2oY8D3c94u1WqvDD/Yth2kPwjt1rcs2RzbZnVRsUjWqFPdcWQ2AV2duIDvP/lV7vVZENm3ahNvtZsSIEaxfv553332X4cOH8/zzz3vrKUVE5E+Vmlr3r3liA3R6EcLjrHEki4fBx62sm+st/ghSD9mdVArZwI41KRcayO5jGXz+206742CYhXiR6P/+7//45JNP2LFjR772T0lJITw8nOTkZMLCwrycTkSkGHO7YOtc6x42W+dYl24ADAdUaweXdYM610NYBXtzSqH4dtU+npi4mpAAJz8/1YGY8KAC/f6evH8X6hiR5ORkIiMjz/nn2dnZpKSknPYQEZEC4HBCnWugz9fw1BZr+fjYFmC6Ycd8mPmUdelmZGf49R04tA7cbrtTi5d0j69E08oRZOS4ePOHjbZmKbQism3bNj788EMeeOCBc+4zdOhQwsPDTz3i4uIKK56ISMlRqqy1fPy9P1lLyHf5j1VKAPavgHkvw/Ar4O1aMOVuSBgDx3dCEZllIZfOMAz+c2N9HIa17XLb99p6fGlm0KBBvPXWW+fdZ+PGjdStW/fUx/v376d9+/Z06NCBzz777Jxfl52dTXZ29qmPU1JSiIuL06UZEZHCkHoINs+CzT/ArkWQm376n4dEQWxzqNTcGoNSobG1BL34rN3H0qlStlSBf19PLs14XESSkpI4duzYefepXr06AQEBABw4cIAOHTpw+eWXM3r0aByO/J+E0RgRERGb5OVYZ0d2LLAu3exPAPdZlgYPKWutXxJdx1rpNTwOwitBWKxVUgyj0KOL/bxaRDyxf/9+OnbsSLNmzfjqq69wOp0efb2KiIhIEZGbBYfWWuVk3wqrmJzYBZznLcQZCKWiITgCgsKtR2AoOP2tacV/PpwB4B9sPfz++G9IWShdziozpcpBQCmVGh9SJIrI/v376dChA1WqVGHMmDGnlZCYmJh8fQ8VERGRIiwnA45thaTNkLQJjm6FlP2QvB/SDnPekuKpgNIQWQ3KVIPI6taj3GXWzQADCv7SglwaT96//bwVYu7cuWzbto1t27YRGxt72p8VlWVlRUTkEgSEWONEKjQ+88/yciD1IGQctW7Ed+qRYl3icbusKcTuPMjLhrwsyM20HjnpkHEM0o9AWhLkZUJOmnVG5tDafzyRAWVrQoVGVo7KraFCPPgFFMb/ASkAhbqOiKd0RkRERMhOs860HN8JJ3bC8R1wbDscXg9pZ1mQzS/YGlRbpQ1Uaw9xrcDptX93y1kUiUszBUFFREREzivtCBxaY50p2bcC9iy2zqb8XVA41OgEtbpCzS5QOtqerCWIioiIiJRMpglHt8Du32HXb7D9Z8g8/rcdDKja1roPT72bIOTci2zKxVMRERERAWssyv4E2PKjtbT9oTV//ZnDD6p3hEa9rCXu/YPty1nMqIiIiIiczYldsP5bWPfN6QNfgyKgcR9odqc1G0cuiYqIiIjIhRzdCmunQOI4SN771+fjWkHL+6Fedw1yvUgqIiIiIvnldlljSRJGW8vbmy7r8+Fx0OpBaHoHBOk9yBMqIiIiIhcj9TAkjIJlI601UAACw6DZAGjzL824yScVERERkUuRmwlrJsLij6xZOAD+IdDiXrjiMd3s7wJURERERAqC223NtlnwFhxYaX3OvxS0vM8qJJr+e1aevH/n/1a4IiIiJY3DAXWugft+hj4TreXjc9Nh0XvwQTz8Psxaol4umoqIiIjIhRiGVUjunw99JkC5+ta9c+a8AB+1hPXTrMXUxGMqIiIiIvllGFDnWnjwV7jxQyhd3lqbZPKdMOpaOLTO7oQ+R0VERETEUw6nNa330ZXQ/t/Wjfb2LIYR7eDHF6wb9Um+qIiIiIhcrMDS0PF5eDQBLrvRWoNk8TDrcs2GGbpckw8qIiIiIpcqvBL0Hgt9p0BEFUjZD5P6w4TbrbVJ5JxURERERApKratg4FJo9ww4/GHzLPi4lbWUvM6OnJWKiIiISEHyD4ZOL8IDC6BCY8g8Ad/cAxP7QdoRu9MVOSoiIiIi3lC+Ptw7Dzq+AA4/2PQ9fNQKNs2yO1mRoiIiIiLiLU5/aP+stf5I+YaQeRwm9IFZz0Jult3pigQVEREREW+LaWitztr6EevjZSPgs86QtNneXEWAioiIiEhh8AuArq9bM2tCouDwOhjRHlZ9ZXcyW6mIiIiIFKZaV8FDi6B6B8jLhOkD4bvHSuw9a1RERERECltoDPT7Fjq+CBiQMNpaIj55n93JCp2KiIiIiB0cDmj/jHWpJigC9idYl2p2LrQ7WaFSEREREbFTrS7WmiMxDSHjKHzZHZZ+aneqQqMiIiIiYrcyVeHuOdDoNut+NT88A7OeAVee3cm8TkVERESkKAgIgR7Doct/rI+XfQpf94asFFtjeZuKiIiISFFhGND2Ceg1FvyCYdtP8PnVcGK33cm8RkVERESkqKl3I9w1C0rHQNJG+KwLHFxtdyqvUBEREREpiio1tVZjLd8Q0o/AqOuL5YwaFREREZGiKrwS3DUTqrSFnFT46mbYMN3uVAVKRURERKQoCwqHft/AZd3AlQOT7oTln9udqsCoiIiIiBR1/kFw6xhoNgAwYeaT8Ov/7E5VIFREREREfIHDCTe8B+2etT6e9wr88gaYpq2xLpWKiIiIiK8wDOj0wl9rjSx4C376j0+XERURERERX9P2Ceg61Npe9B78+LzPlhEVEREREV/U+mG4/o9xIks+hllP+2QZURERERHxVS3uhRs/BAxY/hnMfs7nyoiKiIiIiC9regfc9JG1vfQTmPeyT5URFRERERFf16TvX5dpfnsXFvzX3jweUBEREREpDlrcC1e/bm3PfwMWvW9vnnxSERERESku2jwCnV60tucOtsaNFHEqIiIiIsVJu2fgyqet7ZlPw/pv7c1zASoiIiIixU2nF6H53YAJU++HHQvsTnROKiIiIiLFjWHAdW//daO8CX3h4Gq7U52VioiIiEhx5HBCz8+gSlvISYWvboHjO+xOdQYVERERkeLKPwj6jIfyDSH9CIztCelH7U51GhURERGR4iwoHPpNgYjKcGInTLgdcrPsTnWKioiIiEhxFxoDfadAYDjsXQrTBxaZ1VdVREREREqC6DrQeyw4/GDdFJg/1O5EgIqIiIhIyVG9PdzwrrW94C1YPdHePKiIiIiIlCxN74ArHre2ZzwCu3+3NY6KiIiISEnTeQjUu8laY2Tq/ZCXY1sUP9ueWUREROzhcECPEeB2QftnwS/AtigqIiIiIiWRfzDcNs7uFLo0IyIiIvZRERERERHbeLWI3HjjjVSuXJmgoCAqVKhA//79OXDggDefUkRERHyIV4tIx44dmTRpEps3b+abb75h+/bt3HLLLd58ShEREfEhhmkW3hqvM2bMoHv37mRnZ+Pv73/B/VNSUggPDyc5OZmwsLBCSCgiIiKXypP370KbNXP8+HHGjRtHmzZtzllCsrOzyc7OPvVxSkpKYcUTERERG3h9sOq///1vSpUqRdmyZdmzZw/Tp08/575Dhw4lPDz81CMuLs7b8URERMRGHheRQYMGYRjGeR+bNm06tf8zzzzDqlWrmDNnDk6nkzvuuINzXQ167rnnSE5OPvXYu3fvxR+ZiIiIFHkejxFJSkri2LFj592nevXqBAScuUrbvn37iIuL4/fff6d169YXfC6NEREREfE9Xh0jEh0dTXR09EUFc7vdAKeNAxEREZGSy2uDVZcuXcry5ctp27YtZcqUYfv27bz00kvUqFEjX2dDREREpPjz2mDVkJAQpk6dSufOnalTpw733HMPjRo1YsGCBQQGBnrraUVERMSHeO2MSMOGDfn555+99e1FRESkGCjSd9/9cxyt1hMRERHxHX++b+dnPkyRLiKpqakAWk9ERETEB6WmphIeHn7efQp1iXdPud1uDhw4QGhoKIZhFOj3TklJIS4ujr179xbLqcE6Pt9X3I+xuB8fFP9j1PH5Pm8do2mapKamUrFiRRyO8w9HLdJnRBwOB7GxsV59jrCwsGL7AwY6vuKguB9jcT8+KP7HqOPzfd44xgudCfmT15d4FxERETkXFRERERGxTYktIoGBgQwZMqTYrmmi4/N9xf0Yi/vxQfE/Rh2f7ysKx1ikB6uKiIhI8VZiz4iIiIiI/VRERERExDYqIiIiImIbFRERERGxTbEuIh999BFVq1YlKCiIVq1asWzZsvPuP3nyZOrWrUtQUBANGzZk1qxZhZT04nhyfKNHj8YwjNMeQUFBhZjWMwsXLqRbt25UrFgRwzCYNm3aBb9m/vz5NG3alMDAQGrWrMno0aO9nvNieXp88+fPP+P1MwyDQ4cOFU5gDw0dOpQWLVoQGhpKuXLl6N69O5s3b77g1/nS7+DFHKMv/R5+8sknNGrU6NRCV61bt+aHH34479f40uvn6fH50mt3Nm+++SaGYfD444+fdz87XsNiW0QmTpzIk08+yZAhQ1i5ciWNGzema9euHDly5Kz7//777/Tp04d77rmHVatW0b17d7p37866desKOXn+eHp8YK2cd/DgwVOP3bt3F2Jiz6Snp9O4cWM++uijfO2/c+dOrr/+ejp27EhiYiKPP/449957Lz/++KOXk14cT4/vT5s3bz7tNSxXrpyXEl6aBQsWMHDgQJYsWcLcuXPJzc3l6quvJj09/Zxf42u/gxdzjOA7v4exsbG8+eabJCQksGLFCjp16sRNN93E+vXrz7q/r71+nh4f+M5r90/Lly9nxIgRNGrU6Lz72fYamsVUy5YtzYEDB5762OVymRUrVjSHDh161v179eplXn/99ad9rlWrVuYDDzzg1ZwXy9PjGzVqlBkeHl5I6QoWYH777bfn3efZZ58169evf9rnevfubXbt2tWLyQpGfo7vl19+MQHzxIkThZKpoB05csQEzAULFpxzH1/7Hfyn/ByjL/8emqZplilTxvzss8/O+me+/vqZ5vmPz1dfu9TUVLNWrVrm3Llzzfbt25uPPfbYOfe16zUslmdEcnJySEhIoEuXLqc+53A46NKlC4sXLz7r1yxevPi0/QG6du16zv3tdDHHB5CWlkaVKlWIi4u7YPP3Nb70+l2K+Ph4KlSowFVXXcWiRYvsjpNvycnJAERGRp5zH19/DfNzjOCbv4cul4sJEyaQnp5O69atz7qPL79++Tk+8M3XbuDAgVx//fVnvDZnY9drWCyLyNGjR3G5XJQvX/60z5cvX/6c19QPHTrk0f52upjjq1OnDl988QXTp0/nq6++wu1206ZNG/bt21cYkb3uXK9fSkoKmZmZNqUqOBUqVGD48OF88803fPPNN8TFxdGhQwdWrlxpd7QLcrvdPP7441xxxRU0aNDgnPv50u/gP+X3GH3t93Dt2rWULl2awMBAHnzwQb799lvq1at31n198fXz5Ph87bUDmDBhAitXrmTo0KH52t+u17BI331XCk7r1q1Pa/pt2rThsssuY8SIEbz66qs2JpP8qFOnDnXq1Dn1cZs2bdi+fTvvvvsuY8eOtTHZhQ0cOJB169bx22+/2R3Fa/J7jL72e1inTh0SExNJTk5mypQp3HnnnSxYsOCcb9a+xpPj87XXbu/evTz22GPMnTu3yA+qLZZFJCoqCqfTyeHDh0/7/OHDh4mJiTnr18TExHi0v50u5vj+yd/fnyZNmrBt2zZvRCx053r9wsLCCA4OtimVd7Vs2bLIv7k/8sgjfP/99yxcuJDY2Njz7utLv4N/58kx/lNR/z0MCAigZs2aADRr1ozly5fz/vvvM2LEiDP29cXXz5Pj+6ei/tolJCRw5MgRmjZteupzLpeLhQsXMmzYMLKzs3E6nad9jV2vYbG8NBMQEECzZs2YN2/eqc+53W7mzZt3zut/rVu3Pm1/gLlz5573eqFdLub4/snlcrF27VoqVKjgrZiFypdev4KSmJhYZF8/0zR55JFH+Pbbb/n555+pVq3aBb/G117DiznGf/K130O32012dvZZ/8zXXr+zOd/x/VNRf+06d+7M2rVrSUxMPPVo3rw5ffv2JTEx8YwSAja+hl4dCmujCRMmmIGBgebo0aPNDRs2mPfff78ZERFhHjp0yDRN0+zfv785aNCgU/svWrTI9PPzM99++21z48aN5pAhQ0x/f39z7dq1dh3CeXl6fC+//LL5448/mtu3bzcTEhLM2267zQwKCjLXr19v1yGcV2pqqrlq1Spz1apVJmC+88475qpVq8zdu3ebpmmagwYNMvv3739q/x07dpghISHmM888Y27cuNH86KOPTKfTac6ePduuQzgvT4/v3XffNadNm2Zu3brVXLt2rfnYY4+ZDofD/Omnn+w6hPN66KGHzPDwcHP+/PnmwYMHTz0yMjJO7ePrv4MXc4y+9Hs4aNAgc8GCBebOnTvNNWvWmIMGDTINwzDnzJljmqbvv36eHp8vvXbn8s9ZM0XlNSy2RcQ0TfPDDz80K1eubAYEBJgtW7Y0lyxZcurP2rdvb955552n7T9p0iSzdu3aZkBAgFm/fn1z5syZhZzYM54c3+OPP35q3/Lly5vXXXeduXLlShtS58+f01X/+fjzmO68806zffv2Z3xNfHy8GRAQYFavXt0cNWpUoefOL0+P76233jJr1KhhBgUFmZGRkWaHDh3Mn3/+2Z7w+XC2YwNOe018/XfwYo7Rl34P7777brNKlSpmQECAGR0dbXbu3PnUm7Rp+v7r5+nx+dJrdy7/LCJF5TU0TNM0vXvORUREROTsiuUYEREREfENKiIiIiJiGxURERERsY2KiIiIiNhGRURERERsoyIiIiIitlEREREREduoiIiIiIhtVERERETENioiIiIiYhsVEREREbGNioiIiIjY5v8BrRqDYy2mrNgAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import solve_bvp\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def ode(t, z):\n",
    "    return np.vstack((z[1], -np.abs(z[0])))\n",
    "\n",
    "def bc(za, zb):\n",
    "    return np.array([za[0], zb[0] + 2])\n",
    "\n",
    "t_mesh = np.linspace(0, 4, 5)\n",
    "z_guess = np.zeros((2, len(t_mesh)))\n",
    "\n",
    "sol = solve_bvp(ode, bc, t_mesh, z_guess)\n",
    "print(sol)\n",
    "\n",
    "x_plot = np.linspace(0, 4, 100)\n",
    "z_plot = sol.sol(x_plot)\n",
    "\n",
    "plt.plot(x_plot, z_plot[0], label=r\"$z_1$\")\n",
    "plt.plot(x_plot, z_plot[1], label=r\"$z_2$\")\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "471a1bc2",
   "metadata": {},
   "source": [
    "## Two-point boundary value problem solver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "id": "81293271",
   "metadata": {},
   "outputs": [],
   "source": [
    "def free_final_time_2pbvp(a=1.0, b=1.0, N=20):\n",
    "    # Indirect method (solving a two-point boundary value problem).\n",
    "\n",
    "    def ode(t, x_p_tf):\n",
    "        x1, x2, p1, p2, tf = x_p_tf\n",
    "        return tf * np.array([x2, -p2 / b, np.zeros_like(t), -p1, np.zeros_like(t)])\n",
    "\n",
    "    def boundary_conditions(x_p_tf_0, x_p_tf_N):\n",
    "        x1_0, x2_0, p1_0, p2_0, tf_0 = x_p_tf_0\n",
    "        x1_N, x2_N, p1_N, p2_N, tf_N = x_p_tf_N\n",
    "        return np.array([x1_0 - 10, x2_0, x1_N, x2_N, a * tf_N - p2_N**2 / (2 * b)])\n",
    "\n",
    "    return solve_bvp(\n",
    "        fun=ode, \n",
    "        bc=boundary_conditions, \n",
    "        x=np.linspace(0, 1, N + 1),\n",
    "        y=np.array([np.linspace(10, 0, N + 1),\n",
    "                  np.zeros(N + 1),\n",
    "                  np.zeros(N + 1),\n",
    "                  np.zeros(N + 1),\n",
    "                  np.ones(N + 1)]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "id": "abb6af3e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "       message: The algorithm converged to the desired accuracy.\n",
       "       success: True\n",
       "        status: 0\n",
       "             x: [ 0.000e+00  1.667e-02 ...  9.833e-01  1.000e+00]\n",
       "           sol: <scipy.interpolate._interpolate.PPoly object at 0x30d4f4d60>\n",
       "             p: None\n",
       "             y: [[ 1.000e+01  9.992e+00 ...  8.241e-03  0.000e+00]\n",
       "                 [-1.185e-19 -2.196e-01 ... -2.196e-01  0.000e+00]\n",
       "                 ...\n",
       "                 [ 2.993e+00  2.893e+00 ... -2.893e+00 -2.993e+00]\n",
       "                 [ 4.478e+00  4.478e+00 ...  4.478e+00  4.478e+00]]\n",
       "            yp: [[-5.306e-19 -9.833e-01 ... -9.833e-01  0.000e+00]\n",
       "                 [-1.340e+01 -1.295e+01 ...  1.295e+01  1.340e+01]\n",
       "                 ...\n",
       "                 [-5.985e+00 -5.985e+00 ... -5.985e+00 -5.985e+00]\n",
       "                 [ 0.000e+00  0.000e+00 ...  0.000e+00  0.000e+00]]\n",
       " rms_residuals: [ 3.002e-05  3.003e-05 ...  3.003e-05  3.002e-05]\n",
       "         niter: 2"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "free_final_time_2pbvp()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "id": "3fbef2e4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4.47772199, 4.47772199, 4.47772199, 4.47772199, 4.47772199,\n",
       "       4.47772199, 4.47772199, 4.47772199, 4.47772199, 4.47772199,\n",
       "       4.47772199, 4.47772199, 4.47772199, 4.47772199, 4.47772199,\n",
       "       4.47772199, 4.47772199, 4.47772199, 4.47772199, 4.47772199,\n",
       "       4.47772199, 4.47772199, 4.47772199, 4.47772199, 4.47772199,\n",
       "       4.47772199, 4.47772199, 4.47772199, 4.47772199, 4.47772199,\n",
       "       4.47772199, 4.47772199, 4.47772199, 4.47772199, 4.47772199,\n",
       "       4.47772199, 4.47772199, 4.47772199, 4.47772199, 4.47772199,\n",
       "       4.47772199, 4.47772199, 4.47772199, 4.47772199, 4.47772199,\n",
       "       4.47772199, 4.47772199, 4.47772199, 4.47772199, 4.47772199,\n",
       "       4.47772199, 4.47772199, 4.47772199, 4.47772199, 4.47772199,\n",
       "       4.47772199, 4.47772199, 4.47772199, 4.47772199])"
      ]
     },
     "execution_count": 123,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "free_final_time_2pbvp().y[-1]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "63861975",
   "metadata": {},
   "source": [
    "## Exploring indirect single shooting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "dff52c04",
   "metadata": {},
   "outputs": [],
   "source": [
    "import jax\n",
    "import jax.numpy as jnp\n",
    "import matplotlib.pyplot as plt\n",
    "from jax.experimental.ode import odeint\n",
    "\n",
    "from ipywidgets import interact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "id": "3e7b0bdb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def indirect_single_shooting_error(initial_costate_and_final_time, a=1.0, b=1.0):\n",
    "    initial_state = np.array([10., 0.])\n",
    "    initial_costate = initial_costate_and_final_time[:-1]\n",
    "    final_time = initial_costate_and_final_time[-1]\n",
    "\n",
    "    def shooting_ode(state_costate, t):\n",
    "        (x1, x2), (p1, p2) = state_costate\n",
    "        return jnp.array([x2, -p2 / b, 0, -p1])\n",
    "\n",
    "    states, costates = odeint(shooting_ode, (initial_state, initial_costate), jnp.array([0, final_time]))\n",
    "    return jnp.linalg.norm(jnp.append(states[-1], -costates[-1, 1]**2 / (2 * b) +  a * final_time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "id": "dcf30800",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ac32a8e473024cdfa9e84d94567604bb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=4.05, description='final_time', max=8.0, min=0.1), Output()), _dom_cla…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "vv_indirect_single_shooting_error = jax.jit(jax.vmap(jax.vmap(indirect_single_shooting_error)))\n",
    "X, Y = np.meshgrid(np.linspace(-4, 8, 100), np.linspace(-4, 8, 100))\n",
    "\n",
    "\n",
    "@interact\n",
    "def plot_slice(final_time=(0.1, 8.0)):\n",
    "    plt.figure(figsize=(12, 10))\n",
    "    plt.contourf(X,\n",
    "                 Y,\n",
    "                #  vv_indirect_single_shooting_error(np.stack([X, Y, final_time * np.ones_like(X)], -1)))\n",
    "                 np.log(1 + vv_indirect_single_shooting_error(np.stack([X, Y, final_time * np.ones_like(X)], -1))),\n",
    "                 levels=40,\n",
    "                 vmin=0,\n",
    "                 vmax=10)\n",
    "    plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b0ec4f0",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.10"
  },
  "vscode": {
   "interpreter": {
    "hash": "bda33b16be7e844498c7c2d368d72665b4f1d165582b9547ed22a0249a29ca2e"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}