{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bistable gene-regulatory networks\n",
    "\n",
    "## Model definition\n",
    "\n",
    "We consider the piecewise linear system defined in Filippov sense as in [Augier and Yabo (2021)][id1]. In synthetic biology, this model can represent a _genetic toggle switch_, a synthetic regulatory network designed by [Gardner et al. (2000)][id2] in the bacteria E. coli through two genes lacI and tetR mutually repressing each other. Dynamics is given by\n",
    "\n",
    "[id1]: https://hal.archives-ouvertes.fr/hal-03099681\n",
    "[id2]: https://pubmed.ncbi.nlm.nih.gov/10659857/\n",
    "\n",
    "$$\n",
    "\\left\\{ \\begin{array}{l}\n",
    "\\dot x_1 = -\\gamma_1 x_1 + u(t) k_1 s^{-}(x_2, \\theta_2), \\\\\n",
    "\\dot x_2 = -\\gamma_1 x_2 + u(t)k_2 s^{-}(x_1, \\theta_1),\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "where for $j\\in \\{1,2\\}$, $x_j\\in \\mathbb{R}$, and for $\\theta\\in \\mathbb{R}$, $s^{-}(\\cdot,\\theta):\\mathbb{R}\\to \\mathbb{R}$ is such that\n",
    "\n",
    "$$\n",
    "    s^{-}(x,\\theta)= \\left\\{ \\begin{array}{ll}\n",
    "    1 & \\textit{if } x < \\theta, \\\\\n",
    "    0 & \\textit{if } x > \\theta,\n",
    "    \\end{array} \\right.\n",
    "$$\n",
    "\n",
    "and the control $u(\\cdot) \\in L^\\infty ([0,t_f],[u_{\\text{min}},u_{\\text{max}}])$, with $0<u_{\\text{min}}<1 \\leq u_{\\text{max}}$. Such system is piecewise linear in $\\mathbb{R}^2$ and, in particular, is linear in the regular domains\n",
    "\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "&B_{00}=\\left\\{(x_1,x_2)\\in \\mathbb{R}^2 \\mid 0<x_1<\\theta_1, \\ 0<x_2<\\theta_2\\right\\},\\\\\n",
    "& B_{01}=\\left\\{(x_1,x_2)\\in \\mathbb{R}^2 \\mid 0<x_1<\\theta_1, \\ \\theta_2<x_2<\\frac{k_2}{\\gamma_2}\\right\\},\\\\\n",
    "& B_{10}=\\left\\{(x_1,x_2)\\in \\mathbb{R}^2 \\mid \\theta_1<x_1<\\frac{k_1}{\\gamma_1}, \\ 0<x_2<\\theta_2\\right\\},\\\\\n",
    "&B_{11}=\\left\\{(x_1,x_2)\\in \\mathbb{R}^2 \\mid \\theta_1<x_1<\\frac{k_1}{\\gamma_1}, \\ \\theta_2<x_2<\\frac{k_2}{\\gamma_2}\\right\\}.\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "The following image shows the dynamics for the open loop system (i.e. $u \\equiv 1$):\n",
    "\n",
    "<img src=\"openloop.png\" width=436 height=337>\n",
    "\n",
    "## Time-optimal control problem\n",
    "\n",
    "The optimal control problem involves minimizing the time of a state transfer from $B_{10}$ to $B_{01}$, and so is defined for the state $x(t) = (x_1(t),x_2(t))$ as \n",
    "\n",
    "$$\n",
    "\\left\\{\\begin{array}{l}\n",
    "\t\\textit{minimize} \\ t_f  \\\\\n",
    "\t\\dot x = f(x),\\\\\n",
    " x(0) = x_0 \\in B_{10}, \\\\ x_1(t_f)\\in [0,\\theta_1), \\ x_2(t_f)=x_2^f > \\theta_2, \\\\ u(\\cdot) \\in [u_{\\text{min}},u_{\\text{max}}].\n",
    "\\end{array}\\right.\n",
    "$$\n",
    "\n",
    "with $f(x)$ being the right-hand side of the system previously introduced.\n",
    "\n",
    "## Problem regularization\n",
    "\n",
    "Bocop requires $s^-$ to be regularized to a smooth function. We then define, for $x\\in \\mathbb{R}$ and $k\\in \\mathbb{N}$, the Hill function\n",
    "\n",
    "$$\n",
    "\\delta(x_i,\\theta_i,k) = \\frac{\\theta_i^k}{x_i^k + \\theta_i^k},\n",
    "$$\n",
    "\n",
    "which can approximate $s^-$ for large values of $k$ and, when $k \\rightarrow \\infty$, meets\n",
    "\n",
    "$$\n",
    "\t\\lim_{k \\rightarrow \\infty} \\delta(x_i,\\theta_i,k) = \\left\\{ \\begin{array}{ll}\n",
    "\t1 & x_i < \\theta_i, \\\\\n",
    "\t0 & x_i > \\theta_i, \\\\\n",
    "\t1/2 & x_i = \\theta_i.\n",
    "\t\\end{array} \\right.\n",
    "$$\n",
    "\n",
    "System parameters are fixed to $\\gamma_1 = 1.2$, $\\gamma_2 = 2$, $\\theta_1 = 0.6$, $\\theta_2 = 0.4$ and $k_1 = k_2 = 1$, and control bounds are set to $u_{min}=0.5$ and $u_{max}=1.5$. Initial conditions are set to $x_1(0) = 0.8$, $x_2(0)=0.3$ and $x_2(t_f) = 0.7$. The parameter $k$ of the Hill function is set to $k = 200$.\n",
    "\n",
    "<font color=white>[Thumbnail](thumbnail.png)</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[37m// +++DRAFT+++ This class implements the OCP functions\u001b[39;49;00m\r\n",
      "\u001b[37m// It derives from the generic class bocop3OCPBase\u001b[39;49;00m\r\n",
      "\u001b[37m// OCP functions are defined with templates since they will be called\u001b[39;49;00m\r\n",
      "\u001b[37m// from both the NLP solver (double arguments) and AD tool (ad_double arguments)\u001b[39;49;00m\r\n",
      "\u001b[37m//#pragma once\u001b[39;49;00m\r\n",
      "\r\n",
      "\u001b[36m#\u001b[39;49;00m\u001b[36minclude\u001b[39;49;00m \u001b[37m<OCP.h>\u001b[39;49;00m\u001b[36m\u001b[39;49;00m\r\n",
      "\r\n",
      "\u001b[34mtemplate\u001b[39;49;00m <\u001b[34mtypename\u001b[39;49;00m \u001b[04m\u001b[32mVariable\u001b[39;49;00m>\r\n",
      "\u001b[36mvoid\u001b[39;49;00m OCP::finalCost(\u001b[36mdouble\u001b[39;49;00m initial_time, \u001b[36mdouble\u001b[39;49;00m final_time, \u001b[34mconst\u001b[39;49;00m Variable *initial_state, \u001b[34mconst\u001b[39;49;00m Variable *final_state, \u001b[34mconst\u001b[39;49;00m Variable *parameters, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *constants, Variable &final_cost)\r\n",
      "{\r\n",
      "  \u001b[37m// Minimize final time\u001b[39;49;00m\r\n",
      "  final_cost = initial_state[\u001b[34m2\u001b[39;49;00m];\r\n",
      "}\r\n",
      "\r\n",
      "\u001b[34mtemplate\u001b[39;49;00m <\u001b[34mtypename\u001b[39;49;00m \u001b[04m\u001b[32mVariable\u001b[39;49;00m>\r\n",
      "\u001b[36mvoid\u001b[39;49;00m OCP::dynamics(\u001b[36mdouble\u001b[39;49;00m time, \u001b[34mconst\u001b[39;49;00m Variable *state, \u001b[34mconst\u001b[39;49;00m Variable *control, \u001b[34mconst\u001b[39;49;00m Variable *parameters, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *constants, Variable *state_dynamics)\r\n",
      "{\r\n",
      "\tVariable x1 = state[\u001b[34m0\u001b[39;49;00m];\r\n",
      "\tVariable x2 = state[\u001b[34m1\u001b[39;49;00m];\r\n",
      "\tVariable tf = state[\u001b[34m2\u001b[39;49;00m];\r\n",
      "\tVariable u  = control[\u001b[34m0\u001b[39;49;00m];\r\n",
      "\t\r\n",
      "\t\u001b[36mdouble\u001b[39;49;00m g1 = constants[\u001b[34m0\u001b[39;49;00m];\r\n",
      "\t\u001b[36mdouble\u001b[39;49;00m g2 = constants[\u001b[34m1\u001b[39;49;00m];\r\n",
      "\t\u001b[36mdouble\u001b[39;49;00m t1 = constants[\u001b[34m2\u001b[39;49;00m];\r\n",
      "\t\u001b[36mdouble\u001b[39;49;00m t2 = constants[\u001b[34m3\u001b[39;49;00m];\r\n",
      "\t\u001b[36mdouble\u001b[39;49;00m k  = constants[\u001b[34m4\u001b[39;49;00m];\r\n",
      "\t\r\n",
      "\tVariable delta1 = pow(t2,k)/(pow(t2,k) + pow(x2,k)); \u001b[37m// delta(x2,t2,k)\u001b[39;49;00m\r\n",
      "\tVariable delta2 = pow(t1,k)/(pow(t1,k) + pow(x1,k)); \u001b[37m// delta(x1,t1,k)\u001b[39;49;00m\r\n",
      "\r\n",
      "\tstate_dynamics[\u001b[34m0\u001b[39;49;00m] = -g1*x1 + u*delta1;\t\u001b[37m// x1dot\u001b[39;49;00m\r\n",
      "\tstate_dynamics[\u001b[34m1\u001b[39;49;00m] = -g2*x2 + u*delta2;\t\u001b[37m// x2dot\u001b[39;49;00m\r\n",
      "\tstate_dynamics[\u001b[34m2\u001b[39;49;00m] = \u001b[34m0\u001b[39;49;00m;\t\t\t\u001b[37m// tf (constant)\u001b[39;49;00m\r\n",
      "\t\r\n",
      "\t\u001b[37m// Rescale the state\u001b[39;49;00m\r\n",
      "\t\u001b[34mfor\u001b[39;49;00m (\u001b[36msize_t\u001b[39;49;00m i=\u001b[34m0\u001b[39;49;00m; i<\u001b[34m2\u001b[39;49;00m; i++)\r\n",
      "\t\tstate_dynamics[i] *= tf;\r\n",
      "}\r\n",
      "\r\n",
      "\u001b[34mtemplate\u001b[39;49;00m <\u001b[34mtypename\u001b[39;49;00m \u001b[04m\u001b[32mVariable\u001b[39;49;00m>\r\n",
      "\u001b[36mvoid\u001b[39;49;00m OCP::boundaryConditions(\u001b[36mdouble\u001b[39;49;00m initial_time, \u001b[36mdouble\u001b[39;49;00m final_time, \u001b[34mconst\u001b[39;49;00m Variable *initial_state, \u001b[34mconst\u001b[39;49;00m Variable *final_state, \u001b[34mconst\u001b[39;49;00m Variable *parameters, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *constants, Variable *boundary_conditions)\r\n",
      "{\r\n",
      "    boundary_conditions[\u001b[34m0\u001b[39;49;00m] = initial_state[\u001b[34m0\u001b[39;49;00m]; \u001b[37m// x1(0)\u001b[39;49;00m\r\n",
      "    boundary_conditions[\u001b[34m1\u001b[39;49;00m] = initial_state[\u001b[34m1\u001b[39;49;00m]; \u001b[37m// x2(0)\u001b[39;49;00m\r\n",
      "    boundary_conditions[\u001b[34m2\u001b[39;49;00m] = final_state[\u001b[34m0\u001b[39;49;00m];   \u001b[37m// x1(tf)\u001b[39;49;00m\r\n",
      "    boundary_conditions[\u001b[34m3\u001b[39;49;00m] = final_state[\u001b[34m1\u001b[39;49;00m];   \u001b[37m// x2(tf)\u001b[39;49;00m\r\n",
      "}\r\n",
      "\r\n",
      "\u001b[34mtemplate\u001b[39;49;00m <\u001b[34mtypename\u001b[39;49;00m \u001b[04m\u001b[32mVariable\u001b[39;49;00m>\r\n",
      "\u001b[36mvoid\u001b[39;49;00m OCP::pathConstraints(\u001b[36mdouble\u001b[39;49;00m time, \u001b[34mconst\u001b[39;49;00m Variable *state, \u001b[34mconst\u001b[39;49;00m Variable *control, \u001b[34mconst\u001b[39;49;00m Variable *parameters, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *constants, Variable *path_constraints)\r\n",
      "{}\r\n",
      "\r\n",
      "\u001b[36mvoid\u001b[39;49;00m OCP::preProcessing()\r\n",
      "{}\r\n",
      "\r\n",
      "\u001b[37m// ///////////////////////////////////////////////////////////////////\u001b[39;49;00m\r\n",
      "\u001b[37m// explicit template instanciation for template functions, with double and double_ad \u001b[39;49;00m\r\n",
      "\u001b[37m// +++ could be in an included separate file ? \u001b[39;49;00m\r\n",
      "\u001b[37m// but needs to be done for aux functions too ? APPARENTLY NOT !\u001b[39;49;00m\r\n",
      "\u001b[34mtemplate\u001b[39;49;00m \u001b[36mvoid\u001b[39;49;00m OCP::finalCost<\u001b[36mdouble\u001b[39;49;00m>(\u001b[36mdouble\u001b[39;49;00m initial_time, \u001b[36mdouble\u001b[39;49;00m final_time, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *initial_state, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *final_state, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *parameters, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *constants, \u001b[36mdouble\u001b[39;49;00m &final_cost);\r\n",
      "\u001b[34mtemplate\u001b[39;49;00m \u001b[36mvoid\u001b[39;49;00m OCP::dynamics<\u001b[36mdouble\u001b[39;49;00m>(\u001b[36mdouble\u001b[39;49;00m time, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *state, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *control, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *parameters, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *constants, \u001b[36mdouble\u001b[39;49;00m *state_dynamics);\r\n",
      "\u001b[34mtemplate\u001b[39;49;00m \u001b[36mvoid\u001b[39;49;00m OCP::boundaryConditions<\u001b[36mdouble\u001b[39;49;00m>(\u001b[36mdouble\u001b[39;49;00m initial_time, \u001b[36mdouble\u001b[39;49;00m final_time, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *initial_state, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *final_state, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *parameters, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *constants, \u001b[36mdouble\u001b[39;49;00m *boundary_conditions);\r\n",
      "\u001b[34mtemplate\u001b[39;49;00m \u001b[36mvoid\u001b[39;49;00m OCP::pathConstraints<\u001b[36mdouble\u001b[39;49;00m>(\u001b[36mdouble\u001b[39;49;00m time, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *state, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *control, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *parameters, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *constants, \u001b[36mdouble\u001b[39;49;00m *path_constraints);\r\n",
      "\r\n",
      "\u001b[34mtemplate\u001b[39;49;00m \u001b[36mvoid\u001b[39;49;00m OCP::finalCost<double_ad>(\u001b[36mdouble\u001b[39;49;00m initial_time, \u001b[36mdouble\u001b[39;49;00m final_time, \u001b[34mconst\u001b[39;49;00m double_ad *initial_state, \u001b[34mconst\u001b[39;49;00m double_ad *final_state, \u001b[34mconst\u001b[39;49;00m double_ad *parameters, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *constants, double_ad &final_cost);\r\n",
      "\u001b[34mtemplate\u001b[39;49;00m \u001b[36mvoid\u001b[39;49;00m OCP::dynamics<double_ad>(\u001b[36mdouble\u001b[39;49;00m time, \u001b[34mconst\u001b[39;49;00m double_ad *state, \u001b[34mconst\u001b[39;49;00m double_ad *control, \u001b[34mconst\u001b[39;49;00m double_ad *parameters, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *constants, double_ad *state_dynamics);\r\n",
      "\u001b[34mtemplate\u001b[39;49;00m \u001b[36mvoid\u001b[39;49;00m OCP::boundaryConditions<double_ad>(\u001b[36mdouble\u001b[39;49;00m initial_time, \u001b[36mdouble\u001b[39;49;00m final_time, \u001b[34mconst\u001b[39;49;00m double_ad *initial_state, \u001b[34mconst\u001b[39;49;00m double_ad *final_state, \u001b[34mconst\u001b[39;49;00m double_ad *parameters, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *constants, double_ad *boundary_conditions);\r\n",
      "\u001b[34mtemplate\u001b[39;49;00m \u001b[36mvoid\u001b[39;49;00m OCP::pathConstraints<double_ad>(\u001b[36mdouble\u001b[39;49;00m time, \u001b[34mconst\u001b[39;49;00m double_ad *state, \u001b[34mconst\u001b[39;49;00m double_ad *control, \u001b[34mconst\u001b[39;49;00m double_ad *parameters, \u001b[34mconst\u001b[39;49;00m \u001b[36mdouble\u001b[39;49;00m *constants, double_ad *path_constraints);\r\n"
     ]
    }
   ],
   "source": [
    "!pygmentize problem.cpp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import bocop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[EXEC] > ['cmake -DCMAKE_BUILD_TYPE=Debug -DPROBLEM_DIR=/opt/ct-gallery/gallery/examples/bistable -DCPP_FILE=problem.cpp  -DCMAKE_CXX_COMPILER=g++  /opt/ct-gallery/env/lib/python3.7/site-packages/bocop']\n",
      ">\t-- The C compiler identification is GNU 7.5.0\n",
      ">\t-- The CXX compiler identification is GNU 10.2.0\n",
      ">\t-- Detecting C compiler ABI info\n",
      ">\t-- Detecting C compiler ABI info - done\n",
      ">\t-- Check for working C compiler: /opt/ct-gallery/env/bin/x86_64-conda-linux-gnu-cc - skipped\n",
      ">\t-- Detecting C compile features\n",
      ">\t-- Detecting C compile features - done\n",
      ">\t-- Detecting CXX compiler ABI info\n",
      ">\t-- Detecting CXX compiler ABI info - done\n",
      ">\t-- Check for working CXX compiler: /usr/bin/g++ - skipped\n",
      ">\t-- Detecting CXX compile features\n",
      ">\t-- Detecting CXX compile features - done\n",
      ">\t-- Problem path: /opt/ct-gallery/gallery/examples/bistable\n",
      ">\t-- Using CPPAD found at /opt/ct-gallery/env/include/cppad/..\n",
      ">\t-- Using IPOPT found at /opt/ct-gallery/env/lib/libipopt.so\n",
      ">\t-- Found Python3: /opt/ct-gallery/env/bin/python3.7 (found version \"3.7.8\") found components: Interpreter Development Development.Module Development.Embed \n",
      ">\t-- Found SWIG: /opt/ct-gallery/env/bin/swig (found suitable version \"4.0.2\", minimum required is \"4\")  \n",
      ">\t-- Build type: Debug\n",
      ">\t-- Configuring done\n",
      ">\t-- Generating done\n",
      ">\t-- Build files have been written to: /opt/ct-gallery/gallery/examples/bistable/build\n",
      "[DONE] > ['cmake -DCMAKE_BUILD_TYPE=Debug -DPROBLEM_DIR=/opt/ct-gallery/gallery/examples/bistable -DCPP_FILE=problem.cpp  -DCMAKE_CXX_COMPILER=g++  /opt/ct-gallery/env/lib/python3.7/site-packages/bocop']\n",
      "[EXEC] > make\n",
      ">\tScanning dependencies of target bocop\n",
      ">\t[  3%] Building CXX object src/CMakeFiles/bocop.dir/AD/dOCPCppAD.cpp.o\n",
      ">\t[  7%] Building CXX object src/CMakeFiles/bocop.dir/DOCP/dOCP.cpp.o\n",
      ">\t[ 10%] Building CXX object src/CMakeFiles/bocop.dir/DOCP/dODE.cpp.o\n",
      ">\t[ 14%] Building CXX object src/CMakeFiles/bocop.dir/DOCP/dControl.cpp.o\n",
      ">\t[ 17%] Building CXX object src/CMakeFiles/bocop.dir/DOCP/dState.cpp.o\n",
      ">\t[ 21%] Building CXX object src/CMakeFiles/bocop.dir/DOCP/solution.cpp.o\n",
      ">\t[ 25%] Building CXX object src/CMakeFiles/bocop.dir/NLP/NLPSolverIpopt.cpp.o\n",
      ">\t[ 28%] Building CXX object src/CMakeFiles/bocop.dir/OCP/OCP.cpp.o\n",
      ">\t[ 32%] Building CXX object src/CMakeFiles/bocop.dir/tools/tools.cpp.o\n",
      ">\t[ 35%] Building CXX object src/CMakeFiles/bocop.dir/tools/tools_interpolation.cpp.o\n",
      ">\t[ 39%] Building CXX object src/CMakeFiles/bocop.dir/opt/ct-gallery/gallery/examples/bistable/problem.cpp.o\n",
      ">\t[ 42%] Linking CXX shared library /opt/ct-gallery/gallery/examples/bistable/libbocop.so\n",
      ">\t[ 42%] Built target bocop\n",
      ">\tScanning dependencies of target bocopwrapper_swig_compilation\n",
      ">\t[ 46%] Swig compile /opt/ct-gallery/env/lib/python3.7/site-packages/bocop/src/bocopwrapper.i for python\n",
      ">\tLanguage subdirectory: python\n",
      ">\tSearch paths:\n",
      ">\t   ./\n",
      ">\t   /opt/ct-gallery/env/include/python3.7m/\n",
      ">\t   AD/\n",
      ">\t   DOCP/\n",
      ">\t   NLP/\n",
      ">\t   OCP/\n",
      ">\t   tools/\n",
      ">\t   /opt/ct-gallery/env/include/cppad/../\n",
      ">\t   /opt/ct-gallery/env/include/coin/\n",
      ">\t   /opt/ct-gallery/env/include/python3.7m/\n",
      ">\t   /opt/ct-gallery/env/lib/python3.7/site-packages/bocop/src/AD/\n",
      ">\t   /opt/ct-gallery/env/lib/python3.7/site-packages/bocop/src/DOCP/\n",
      ">\t   /opt/ct-gallery/env/lib/python3.7/site-packages/bocop/src/NLP/\n",
      ">\t   /opt/ct-gallery/env/lib/python3.7/site-packages/bocop/src/OCP/\n",
      ">\t   /opt/ct-gallery/env/lib/python3.7/site-packages/bocop/src/tools/\n",
      ">\t   ./swig_lib/python/\n",
      ">\t   /opt/ct-gallery/env/share/swig/4.0.2/python/\n",
      ">\t   ./swig_lib/\n",
      ">\t   /opt/ct-gallery/env/share/swig/4.0.2/\n",
      ">\tPreprocessing...\n",
      ">\tStarting language-specific parse...\n",
      ">\tProcessing types...\n",
      ">\tC++ analysis...\n",
      ">\tProcessing nested classes...\n",
      ">\tGenerating wrappers...\n",
      ">\t[ 46%] Built target bocopwrapper_swig_compilation\n",
      ">\tScanning dependencies of target bocopwrapper\n",
      ">\t[ 50%] Building CXX object src/CMakeFiles/bocopwrapper.dir/CMakeFiles/bocopwrapper.dir/bocopwrapperPYTHON_wrap.cxx.o\n",
      ">\t[ 53%] Linking CXX shared module /opt/ct-gallery/gallery/examples/bistable/_bocopwrapper.so\n",
      ">\t-- Moving python modules to /opt/ct-gallery/gallery/examples/bistable\n",
      ">\t[ 53%] Built target bocopwrapper\n",
      ">\tScanning dependencies of target bocopApp\n",
      ">\t[ 57%] Building CXX object src/CMakeFiles/bocopApp.dir/AD/dOCPCppAD.cpp.o\n",
      ">\t[ 60%] Building CXX object src/CMakeFiles/bocopApp.dir/DOCP/dOCP.cpp.o\n",
      ">\t[ 64%] Building CXX object src/CMakeFiles/bocopApp.dir/DOCP/dODE.cpp.o\n",
      ">\t[ 67%] Building CXX object src/CMakeFiles/bocopApp.dir/DOCP/dControl.cpp.o\n",
      ">\t[ 71%] Building CXX object src/CMakeFiles/bocopApp.dir/DOCP/dState.cpp.o\n",
      ">\t[ 75%] Building CXX object src/CMakeFiles/bocopApp.dir/DOCP/solution.cpp.o\n",
      ">\t[ 78%] Building CXX object src/CMakeFiles/bocopApp.dir/NLP/NLPSolverIpopt.cpp.o\n",
      ">\t[ 82%] Building CXX object src/CMakeFiles/bocopApp.dir/OCP/OCP.cpp.o\n",
      ">\t[ 85%] Building CXX object src/CMakeFiles/bocopApp.dir/tools/tools.cpp.o\n",
      ">\t[ 89%] Building CXX object src/CMakeFiles/bocopApp.dir/tools/tools_interpolation.cpp.o\n",
      ">\t[ 92%] Building CXX object src/CMakeFiles/bocopApp.dir/opt/ct-gallery/gallery/examples/bistable/problem.cpp.o\n",
      ">\t[ 96%] Building CXX object src/CMakeFiles/bocopApp.dir/main.cpp.o\n",
      ">\t[100%] Linking CXX executable /opt/ct-gallery/gallery/examples/bistable/bocopApp\n",
      ">\t[100%] Built target bocopApp\n",
      "[DONE] > make\n",
      "Done\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "problem_path = \".\" # using local problem definition\n",
    "bocop.build(problem_path, cmake_options = '-DCMAKE_CXX_COMPILER=g++')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "State 0 and State 1 correspond to $x_1$ and $x_2$ respectively. State 2 represents the final time to minimize."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6d62a4c4958c4c9994eab32d1b45b6be",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=113, continuous_update=False, description='iteration', max=113), Output(…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Done\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAAEkCAYAAACLy/7uAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZQcZ33v//dn9n0fjTSSRiPJsmx5wTbCCxhjAgHLYBvIcg3kYsA5BgK54QfkYuL7cwgnJIEcciAB4pjgg5OATXJig8mF2PzYbLwiyVotydrX2fd9/f7+qJpRa9Q9M9JMT9e0vq9z+nT1U1Vd33m6pr9dVU89j8wM55xzLkoyUh2Ac845N5UnJ+ecc5Hjyck551zkeHJyzjkXOZ6cnHPORY4nJ+ecc5Hjyck551zkeHJy7jxIulHSc5K6JLVLelbS68J5H5T063N4r3pJJilrDvFcJWmLpP7w+arzfS/nosCTk3PnSFIJ8F/APwAVwHLgL4ChFMWTA/wQ+DegHHgY+GFY7tyi5MnJuXN3MYCZPWJmY2Y2YGZPmdkOSZcCDwA3SOqV1Akg6R2SXpbULem4pM/HvN/T4XNnuM4N4ToflrRHUoekJyWtShDPzUAW8FUzGzKzvwcE/Na8/+XOLRBPTs6du1eBMUkPS9okqXxihpntAT4KPG9mRWZWFs7qAz4AlAHvAD4m6V3hvJvC57JwnefDeX8GvAeoBp4BHkkQz2XADjuzL7IdYblzi5InJ+fOkZl1AzcCBnwLaJH0hKSaadb5pZntNLNxM9tBkGjeNM1mPgL8tZntMbNR4K+AqxIcPRUBXVPKuoDi2f9VzkWLJyfnzkOYND5oZiuAy4Fa4KuJlpd0naRfSGqR1EVwdFU1zSZWAV+T1BmeGmwnOFW3PM6yvUDJlLISoGf2f5Fz0eLJybk5MrO9wHcIkhQER1RTfQ94AlhpZqUE16U0zfLHgY+YWVnMI9/Mnouz7G7gSkmKKbsyLHduUfLk5Nw5knSJpE9LWhG+Xgm8F3ghXKQJWDGltVwx0G5mg5KuBd4XM68FGAfWxJQ9AHxO0mXhNkol/V6CkH4JjAH/S1KupE+E5T8/7z/SuRTz5OTcuesBrgNelNRHkJR2AZ8O5/+c4KilUVJrWPZHwBck9QD3A/8+8WZm1g98EXg2PI13vZk9DnwJeFRSd/j+m+IFY2bDwLsIGlx0Ah8G3hWWO7coyQcbdM45FzV+5OSccy5yPDk555yLHE9OzjnnIseTk3POucg5716Qk6mqqsrq6+tTHYZzzrkk2rJlS6uZVcebF8nkVF9fz+bNm897/Xv/cwdlBTlcsbyUy5eXUFdRwJn3JzrnnEs1SUcTzZtVcpJ0C/A1IBP4ZzP7mynz/xR4f8x7XgpUm1m7pCME94WMAaNmtvGc/4JzMDZu7Gno5pWGbkbGgmbyxXlZXF5byhUrSrmstoQrlpdSX1lIRoYnLOeci6IZ73OSlEnQC/NvAyeA3wDvNbNXEix/G/D/mNlvha+PABvNrDXe8vFs3LjR5nLkBDA0Osarjb3sOtXFzpNd7D7ZxZ7GHoZHxwEoys1iQ20Jl9cGR1eX1ZayprqQ7Ey/DOeccwtB0pZEByyzOXK6FjhgZofCN3sUuAOIm5wIunFJ1LX/gsnNyuSKFcHR0nvDspGxcfY3BQlr18ng8b2XjjI4EiSs7Exx0ZJiLl1azCXLirlkaQmXLCumuijXTws659wCmk1yWk7QCeWEEwRdt5xFUgFwC/CJmGIDnpJkwD+Z2YMJ1r0HuAegrq5uFmGdu+zMDDbUlrChtoTf37gSgNGxcQ629LG3sZs9DT3sbezmuYNtPPbyycn1KgtzJpPV+qXFXLq0hHU1ReRlZyYlTuecu9DNJjnFO2RIdC7wNuBZM2uPKXuDmZ2StAT4qaS9Zvb01BXDpPUgBKf1ZhHXvMjKzGD90mLWLy3mjqtOl3f0DbO3MUhWe8Ok9d0XTx9lSbCyvIB1S4q4KHysqylmbXUhxXnZCxW+c86lpdkkpxPAypjXK4BTCZa9kymn9MzsVPjcLOlxgtOEZyWnqCkvzOGGtZXcsLZysmxs3Dja1seehh72N/ewv7mXA029PLO/leGx8cnllpXmnU5YS4rD5yLKC3Pibco559wUs0lOvwHWSVoNnCRIQO+bupCkUoKRPf8gpqwQyDCznnD6bcAX5iPwVMjMEGuqi1hTXQQsmywfHRvnWHt/kKzCx/7mHh596TgDI2OTy1UV5bCmqoj6qgLqqwpZU1VIfVUh9ZWFforQOedizJiczGw0HB/mSYKm5A+Z2W5JHw3nPxAu+m7gKTPri1m9Bng8bEyQBXzPzP57Pv+AKMjKzJhMWm+/7HT5+LhxsnOAAy3BEdb+5h4Ot/bx870ttPaeOOM9akvzWF0dJKrVVcGjvqqQleUF5GR5C0Ln3IUlkkNmzEdT8qjrGRzhSGs/h1p7OdLaz+HWXg639XO4pZfuwdHJ5TIzxIryfOoqClhZURA8l4fPFfmU5md7S0Ln3KI016bkLgmK87Inm7rHMjM6+kc43NrH4dY+joTPxzv62bWzgY7+kSnvk3VGsqqrKGBFmMSWl+X76ULn3KLkySliJFFRmENFYQ6vXVV+1vyewRGOtw9wrL2fEx39HGvv53h7P/ube/jFvmaGRsfPWH5pSR4ryvOpLZt45FFbGkwvL8unJD/Lj7ycc5HjyWmRKc7LZkNtNhtqS86aNz5utPQOcbx9ImkFSexkZz/bjnfyk10Nk106TSjIyTyduErzzkpiS0vz/OjLzYpZsP8dbO5jeGycmpJc1lQV+TVTd148OaWRjAxRU5JHTUkeG+srzpo/Pm609g1xqnOQU50D4SOc7hrglVNdtPYOn7VeRWEOS4pzqSnJY2lJHjUluSwJtzPxurIol0zvq/CCMjo2zp6GHjYfbWfzkQ42H22nqXvojGXysjPYuKqCG9ZWcv2aSq5cUepdhLlZ8eR0AcnIEEuK81hSnMdVK8viLjM4MkZj10TCCp6bugfDxxB7Grpp7R1ifEo7mgxBdZjAasKEVVMcTpfmUVWUQ3VRLhWFOWT5l9Oi1DM4wsvHOtl8tIPNR9rZdryT/uHgVonlZflct7qSq+vKuGhJ0HtKQ9cgW4928MKhNv72yX0AFOZksrG+guvXBPcQXl5b4vuDi8tb67lzNjo2TlvfMI1dYdLqGaI5TGCN3aenpzbegKBnjfKCHKqKcqgqyj39KA5eV8e8rizM9VNCKdTeN8xLh9t44VA7Lx1uZ29jN+MW/BC5dFkJG1eV89r6CjauKqe2LH/a92rrHeKlw+08f6iN5w+2sb+5F4Di3Cxet7qCG9YER1Ybakv8CPwCMl1rvVklp1kMmXEz8EPgcFj0mJl9YTbrxuPJKT0MjozR0jNEU/cgLT1DtPYO0dI7TGvvEK3h69bw9cQv8KlK87PPSGTlhdlUFORQHjYaKS8InwtzqCjIIT/Hr4+dr46+YV483M4Lh9p44VAbext7AMjPzuSaVWVsXFXBxvpyrq4rpyh3biddWnqGJrfz/KE2DrUEt0eW5GVx7ergqOqGNZVcsrTYh7ZJY3NKTrMZMiNMTp8xs3ee67rxeHK68PQPj9LaM0xL70TSGqK1Z/j0dO8Qbb3DdPQP0zkwQqLdNi87I37yKsihojB7MomVFeRQWpBNaX42hTmZF2SLxc7+IBk9f/DsZLSxvpzr1yzcdaKm7sEgUYWxHGnrB6CsIJvrwiOrG9ZWcXFN0QX5WaWrhR4yY77WdReQgpws6iqzqKssmHHZsXGja2CE9r4gWbX3DdPRN0x7f/jcNzJZfqy9n/a+YXpibmyeKitDlOYHiaokfC7Nz6asIPuM8rKJeTHl+dmLJ7ENj46z5WgHvz7QwjP7W9l5sguz040WPvO2ZdywtpIrlpct+OnUmpI87rhqOXdctRyAU50Dk8nq+UNtPLm7CQhGCLh+TSXXrangmrpyLlla7Nes0tR8Dplxg6TtBJ3CfsbMdp/DugsyZIZLD5kZp+8Fm62RsfHgqKt/ZDKZdQ2MxH109A9zpK2ProERugdGzmr8ESsnMyNMaFkU52VTnJcVPHKD6aK8sDw364zXRblZlITTedkZSUlwZsbBll6efrWVZ/a38OLhdvqHx8jKEFfXlfHJt1zMGy6q5MoVC5+MZlJbls97rlnBe65ZAcDx9v7JU4AvHGzj/+5sAIKjvCtXlHLNqnKuqSvnmroyKotyUxm6myfzNWTGVmCVmfVKuhX4AbBulusGhSkaMsNdGLIzMyZbKp6L8XGjZ2iU7ikJrLP/zNfdAyN0D47QMzjKqc4BeodG6RkcTXgtLVZmhoLElRsvkWVRlJtNUW4mBTnBMoW5WRTmZp6ezgle52Rl0NQ9yO5T3fx6fyvP7G+lsXsQgDVVhfzua1fwxnXVXL+mYtEN67Iy7L7r9zauxCzos3LrsU62Hu3g5WMdfOvpQ4yGvyJWVRZMJqqr/ehq0ZqXITPMrDtm+seSvimpajbrOhdlGTGn/FbOvPhZRsfG6Rsao2coSFw9g6P0xkwHj5HJZDbxurF7kJ7m0/Om3jw9k9L8bG68qIob11Vx40VVrKyY+XTpYiGJFeUFrCgv4PbX1AJB45udJ7vYerSDrcc6+PWBVh4PBwwtyMnkiuWlvGZlGVcsL+XKFaXUVRQsmtOxF6p5GTJD0lKgycxM0rVABtAGdM60rnPpLCszg9KCDEoL5nakMjw6Tt/QKL1Do/QNj4bTY/QNjU4+hkbHqSnJY2110QXXJDsvO5PX1VfwuvDmczPjRMcAW4918PKxTl4+1sF3nj0yOe5aSV4WV64o44oVpVy5POjjcnlZviesCJmvITN+F/iYpFFgALjTgmaAcddN0t/iXNrKycogJyvHB6ycJUmTpwInGlkMj47zalMPO092seNEFztPdp5xOrCiMCc4wlpRyuXLS7l0WQkryj1hpYrfhOucu2ANjoyxr7GHHSe72Hmikx0nutjf3MtYmLBK8rK4ZFkJG8LHpctKWFdT5P1NzhMfMsM55+LIy87kNSvLeM3KMmAVECSsVxq62TP56OHfNx+fbNySmSHWVhdyaUzCunRZCVVFOX6UNY88OTnnXIy87Mywtd/pIWvGx41j7f2TSeuVU9385nA7P9x2un1XeUE262qKubimiItrilm3JJj2pu3nx5OTc87NICND1FcVUl9VyK1XLJss7+wf5pWGbvY29LC/uYdXm3r54bZTZ9z0XVmYw7qJhFVTzLolwfS53Kd3IfLk5Jxz56msIIfXr63i9WurJsvMjKbuIV5t6uHVph72N/XyanMPj209Se/Q6aRVVZTDRUuKWFNdxJqqQtZWF7G6qpAV5fl+XxaenJxzbl5JYmlpHktL87jp4urJcjOjoWtwMmHtb+5hf3MvP97ZQGdMD/7ZmWJVZSGrqwpZU13I2qoiVlcXsqaqkIrCC+e6licn55xbAJImR5q+ef2SM+a19w1zuLWXgy19HGrp43BrL4da+vjVvpbJe7MguLl6dVWQuOoqClhVWUBdRQF1lQVUF+WmVeKaVXKaxZAZ7wc+G77sBT5mZtvDeUeAHmAMGE3UbNA55y5UQV+RFbx21ZkjWI+NGyc7BjgYJquJpPXS4XZ+sO3kGb3zF+RkUhfe27VqInFVBklseVl+5PpPnMmMySkc9uIbxAx7IemJKcNeHAbeZGYdkjYR9JEX28Hrm82sdR7jds65tJeZIeoqgyOjN68/c97Q6BgnOgY41tbPsfZ+jrb1c6y9jyOtfTz9agtDo6ePuDIUdKa7ojyf5WUFLC/PZ0VZPsvL81lels+ysjxys6J179a8DJlhZs/FLP8CQR96zjnnkiQ3K5O11UWsrS46a974uNHSO8TRtn6OtvVxvL2fo+39nOgY4NkDrTT1DJ5x1CXBkuJclpfls7y8IHw+M4EVznGAyXM1n0NmTLgb+EnMawOekmTAP4W9jzvnnEuSjAxRU5JHTUke166uOGv+8Og4jV2DnOjs52THACc7Byafd5zo5L93NZzV2XBZQTbLSvOpLc1jWVket1y2jBvXVZ313vNlvobMCBaU3kyQnG6MKX6DmZ2StAT4qaS9ZvZ0nHV9PCfnnFsAOVkZk6cL45k48joRJqwTHf2c6hygoXOQhq5Bth7roL6yMOXJaVbDXki6EvhnYJOZtU2Um9mp8LlZ0uMEpwnPSk4+npNzzkVD7JHXa1eVx10m2f2yztjxq6Qs4FXgLQTDXvwGeF9s7+KS6oCfAx+Ivf4kqRDIMLOecPqnwBfM7L9n2GYLcPT8/qRJVcBiaYThsSaHx5ociylWWFzxXmixrjKz6ngz5mvIjPuBSuCbYTv7iSbjNcDjYVkW8L2ZElP4nnGDPReSNi+WZusea3J4rMmxmGKFxRWvx3rarJpfmNmPgR9PKXsgZvoPgT+Ms94h4DVzjNE559wFZnHdleWcc+6CkM7JaTE1WfdYk8NjTY7FFCssrng91lAkR8J1zjl3YUvnIyfnnHOLlCcn55xzkZN2yUnSLZL2STog6d5UxxNL0kpJv5C0R9JuSX8Sln9e0klJ28LHramOFYIe5SXtDGPaHJZVSPqppP3hc/w79BaYpPUx9bdNUrekT0albiU9JKlZ0q6YsoR1Kelz4T68T9LbIxDr30raK2mHpMcllYXl9ZIGYur3gcTvvGCxJvzMI1iv34+J84ikbWF5qus10XfVwu2zZpY2D4L7sA4Ca4AcYDuwIdVxxcS3DLgmnC4muLl5A/B54DOpji9OvEeAqillXwbuDafvBb6U6jgT7AeNwKqo1C1wE3ANsGumugz3ie1ALrA63KczUxzr24CscPpLMbHWxy4XkXqN+5lHsV6nzP8KcH9E6jXRd9WC7bPpduQ02YO6mQ0DEz2oR4KZNZjZ1nC6B9hD0LHuYnIH8HA4/TDwrhTGkshbgINmNtdeRuaNBf1Jtk8pTlSXdwCPmtmQmR0GDhDs2wsiXqxm9pSZTYwxHpmRBxLUayKRq9cJCnoq+H3gkYWKZzrTfFct2D6bbskpXg/qkfzyl1QPXA28GBZ9Ijxl8lBUTpVxukf5LWHHvAA1ZtYAwQ4MLEm4durcyZn/5FGsW0hcl1Hfjz/MmSMPrJb0sqRfSXpjqoKaIt5nHuV6fSPQZGb7Y8oiUa9TvqsWbJ9Nt+Q06x7UU0lSEfCfwCfNrBv4R2AtcBXQQHB4HwVvMLNrgE3AxyXdlOqAZiIpB7gd+I+wKKp1O53I7seS7gNGge+GRQ1AnZldDXwK+J6kklTFF0r0mUe2XoH3cuYPqkjUa5zvqoSLximbU92mW3KaVQ/qqSQpm+DD/q6ZPQZgZk1mNmZm48C3WMBTDdOxmB7lgYke5ZskLQMIn5tTF2Fcm4CtZtYE0a3bUKK6jOR+LOku4J3A+y280BCexmkLp7cQXGu4OHVRTvuZR7Ves4D3AN+fKItCvcb7rmIB99l0S06/AdZJWh3+gr4TeCLFMU0Kzyt/G9hjZn8XU74sZrF3A7umrrvQJBVKKp6YJrggvougPu8KF7sL+GFqIkzojF+gUazbGInq8gngTkm5klYD64CXUhDfJEm3AJ8Fbjez/pjyakmZ4fQaglgPpSbKyZgSfeaRq9fQW4G9ZnZioiDV9Zrou4qF3GdT1RokWQ/gVoKWJQeB+1Idz5TYbiQ41N0BbAsftwL/CuwMy58AlkUg1jUErW+2A7sn6pKg9/mfAfvD54pUxxoTcwHQBpTGlEWibgkSZgMwQvAr8+7p6hK4L9yH9xGMkZbqWA8QXFOY2G8fCJf9nXD/2A5sBW6LQKwJP/Oo1WtY/h3go1OWTXW9JvquWrB91rsvcs45FznpdlrPuQUh6UZJz0nqktQu6VlJrwvnfVDSr8/hveolWXjt4XzjeTC8+XFc0gfP932ciwpPTs6do7DV1H8B/wBUEDSZ/QtgKIVhbQf+iOAUkHOLnicn587dxQBm9ogFrcIGLLhJdYekS4EHgBsk9UrqBJD0jvCelW5JxyV9Pub9ng6fO8N1bgjX+XDYfUyHpCclrUoUkJl9w8x+Bgwm4w92bqF5cnLu3L0KjEl6WNKm2Bt7zWwP8FHgeTMrMrOycFYf8AGgDHgH8DFJE3fXT9w/Vhau83w4788ImhhXA88Qkd4DnFsInpycO0cW3Iw40ZrpW0CLpCck1Uyzzi/NbKeZjZvZDoJE86ZpNvMR4K/NbI8F3Qb9FXDVdEdPzqUTT07OnYcwaXzQzFYAlwO1wFcTLS/purCX5xZJXQRHV1XTbGIV8DVJneGpwXaCu/Cj0t2Oc0nlycm5OTKzvQT3qlw+URRnse8R3HOz0sxKCa5LaZrljwMfMbOymEe+mT03v9E7F02enJw7R5IukfRpSSvC1ysJeqZ4IVykCVgR9lIyoRhoN7NBSdcC74uZ1wKME9z4POEB4HOSLgu3USrp96aJKUdSHkHCy5aUJ8n/v92i5Tuvc+euB7gOeFFSH0FS2gV8Opz/c4K7+xsltYZlfwR8QVIPcD/w7xNvZkF3QF8Eng1P411vZo8TjJv0qKTu8P03TRPTU8AA8HrgwXA68h31OpeI9xDhnHMucvzIyTnnXOR4cnLOORc5npycc85Fjicn55xzkePJyTnnXOScdxf9yVRVVWX19fWpDsM551xodNwYHTPGxscZHTdyszLJy57b8c2WLVtazaw63rxIJqf6+no2b96c6jCccy5tjY6N0943THPPEC29Q7SGzy09waO1d4i23mHa+obp6B9m6l1Hf/zbF/O/3rJuTjFIOppoXiSTk3POufMzMDxGY/cgDV0DNHefTjQtU5JPe5yEA1CUm0V1cS5VRTlctKSIawtzqCzMoaIwh8qi3GC6KIdlpflJ/Ts8OTnn3CJgZnT0j9DYNUhT9yANXYM0dg/S1DVIw8Rz1wDdg6NnrZuTlcGS4lyqinJZWVHANavKqS7Kpao4l+qiXKqLcyfn5+dkpuCvO5snJ+eci4CRsXEauwY50THAyc4BTnT0c7JjgBMdA5zqGqCha5Dh0fEz1pGgqiiXZaV51FUWcN2aCmpK8lhaksey0jyWlOSxpCSX4twsJCXYcjR5cnLOuQUwMjY+mWxOdPSHCWggLOunsXuQ8Smn2ZYU57KiPJ8rlpfy9suWUhMmnYnn6uJcsjPTs9G1JyfnnJsngyNjHGvv52hbP0fb+jjS1hdOB8loLCb7ZAiWleazvCyf69dUsrw8nxXl+SwvK2B5eT7LSvPIy47GKbZU8OTknHPnYGRsnKNt/Rxs6eVgSy9HW/s50tbHsfZ+GroGz1i2JC+L+qpCXrOyjNtfU0tdZQF1FQUsL8tnaWle2h71zAdPTs45F0f34AgHm3s52NIXJKLmMBm19TMacwRUVZRLfWUBN6ytpL6ykFWVBayqLKS+soCygpxptuCm48nJOXdB6x4cYV9jD3sbe3i1sYcDYRJq7hmaXCYrQ9RXFXLRkiLeftlSLlpSxNrqItZUF1Kcl53C6NPXrJKTpFuArwGZwD+b2d9MmX8z8EPgcFj0mJl9YTbrOufcQhgZG+dQSx97G7vZ29jDvvBxsnNgcpni3CwuqinipourWVtdxNrqICGtrCjwU3ALbMbkJCkT+Abw28AJ4DeSnjCzV6Ys+oyZvfM813XOuXnT0TfMrlNd7D7Vzd6GIBkdbOllZCw4HZeVIdZWF/HaVeW8//o6LllazPqlJdSW5i26JtfpajZHTtcCB8zsEICkR4E7gNkkmLms65xzM2rtHWLXya7w0c3Ok11nHA3VluaxfmkxN69fwqXLilm/tJg1VUXkZPmRUJTNJjktB47HvD4BXBdnuRskbQdOAZ8xs93nsK5zzs2os3+Yl493suN4FztPdrH7VNcZLeRWVxVydV0ZH7hhFVcsL+Wy2lJKC/ya0GI0m+QU7xh3ao9MW4FVZtYr6VbgB8C6Wa4bbES6B7gHoK6ubhZhOefS2ejYOHsbe3j5eCcvH+tg27FODrX2AUHPCGuqCrl2dcVkErpseQkl3jghbcwmOZ0AVsa8XkFwdDTJzLpjpn8s6ZuSqmazbsx6DwIPAmzcuDFuAnPOpa/m7kG2Huvg5WOdvHy8k50nuhgYGQOgqiiHq1aW8zuvXcHVdWVcuaKMolxvbJzOZvPp/gZYJ2k1cBK4E3hf7AKSlgJNZmaSriUYxLAN6JxpXefchcfMONExwIuH23npcBsvHW7nSFs/ANmZ4rLaUu68diVX15Vz9coyVpTne0OFC8yMycnMRiV9AniSoDn4Q2a2W9JHw/kPAL8LfEzSKDAA3GlmBsRdN0l/i3MuosyMA829YTIKHo3dwbWi0vxsXldfwfuuq2NjfQUblpVc0N32uIAs3oAeKbZx40bzwQadW7xaeobYdaqLV051s/14J5uPdtDeNwxAdXEu162u4LrVFbxudQUXLykmI8OPii5EkraY2cZ48/ykrXNuTsbGjVdOdfPrA6385kg7u0910dR9uneFVZUFvHn9Eq5bXcG1qytYVVngp+jcjDw5OefO2fH2fp7Z38qvD7Tw3ME2OvtHAFhbXcjr11ZxWW0Jl9WWsqG2hNJ8b0Hnzp0nJ+fcjIZHx9l8tJ1f7mvh53ubOdDcC8Cy0jzeemkNN15UxevXVrKkJC/Fkbp04cnJORdXc/cgv9zXwi/2NfPM/lZ6h0bJyczgujUVvPfaOt50cTVrqwv9FJ1LCk9OzjkgaFG3t7GHp3Y38dM9jew6Gdy+uKw0j9teU8ub11fzhouqKPT7i9wCSMu9bHzcvPWPc7MwNm68fKyDJ3c38uTuJo619yPBa+vK+d+3rOe3LlnC+ppiPzpyCy4tk9Obv/JL8rMzWVdTzPqaovC5mJUVBWR60nIXuKHRMZ472MZTuxv56StNtPYOk5OZwesvquRjN6/lrZfWUF2cm+ow3QVuvsZzej/w2fBlL/AxM9sezjsC9ABjwGiiNu3zZXzc2HT5Ml5t6mHr0Q5+tP10b0l52RlctKSIi2uKuThMWBcvLfZu8l3aGxkb57mDbfxo+yme3N1Iz+AohTmZ3HzJEt5+2VLevL7aB81zkTJf4zkdBjP5oGcAAB8dSURBVN5kZh2SNhH0kRfb+/ibzax1HuNOKCND3LvpksnXvUOj7G/qYX9TL/uaeni1qYdnD7Ty2NaTk8sU5WaxtrqQNdVFrK4qZE11IWuqgun8HL9T3S1O4+PGS0fa+dH2U/xkVyPtfcMU52bx9suXcusVS3n92irvicFF1ryM52Rmz8Us/wJBB6+RUJSbFfTPVVd+RnlX/wivNgfJ6tXGHg629PHioTYef/nkGcvVluadmbSqi1hTVUhtWb6fInSRMz5u7DjZxRPbTvF/d56iqXuI/OxM3rqhhtuuXMZNF1d7QnKLwnyO5zThbuAnMa8NeEqSAf8U9j5+loUeMqO0IOjP63X1FWeUDwyPcbi1j0OtvRxu6eNQa/D4wbaT9AyOTi6XnSlWlBewsqKAuop86ioKqKuYeF3gp0hcXGbG0Og4gyNjDI2OMzQyztBoOD06Fr6OKTtj/jhDIzHTU5bvHhxlz6luesIm3zevr+a219TylkuXUJCTlpeXXRqbr/GcggWlNxMkpxtjit9gZqckLQF+KmmvmT191htGZMiM/JxMNtSWsKG25IxyM6O1dzhIXC29HG3v51h7P8fb+9lxonPyDvkJ5QXZQcKqLJxMXivLC1hRXkBNaS65Wf7rNWrMjIGRMfqGxugbGmVgZIyBkTEGh8dOT4+Mn1U2MDzGYMz0wEiQNOLOHxljrt1Z5mVnkJuVSU5WBrmTj0wKcjK5/apaXruqnLdcWuM9M7hFbV7GcwKQdCXwz8AmM2ubKDezU+Fzs6THCU4TnpWcok4S1cW5VBfncu3qirPmdw2McDxMWMemJK6f7GxgdPzMb6Tq4lxqy/JZXpbHstL8yenasnyWleZTVZTjjTRmMDZu9A6N0jM4Qt/QGL1Do/QPj9I3NErf0Bj9w6P0hs8TCadvYv7wmeX9w2P0DY+ec+LIycwgLzuD/JxM8rMzycvOJD8nk7ysTKqLc2PKMsjPDpcJ5+fnZJKblREmmczTiSY7wXRWJtmZ8v3CXRDmazynOuAx4H+a2asx5YVAhpn1hNNvA74wX8FHSWl+NqXLS7l8eelZ80bHxmnoGuR4ez8nOwc41TnIqc4BTnUNsK+xh1/sbZkcVG1CTlYGtaWnk9XS0lyWFOdRU5LLkpI8akryqC7KJScrY6H+xHk1Pm70hImlZ3A0fIxMPncPjtKdYN5EWe/Q6MwbIhg1tTAni8LcTApzsigIn5cU51FQmUlRbhYFOVkU5WZSkJtFYU4mBTlZFOQEiWQiqZydgDLIylyc9e9c1M3XeE73A5XAN8NfdRNNxmuAx8OyLOB7ZvbfSflLIiwrM4OV4fWoeMyMzv6RMHEN0NAVJK+J188dbKW5Z4ix8bN/1lcU5rCkOExYxbnUlAQJrDpMZDUleVQlMYmNjRs9gyN09o/Q0T9M58AIXf0jdPYP09E/QtdAMN05ECwzMd01MDLjUUpOZgbFeVnhI5vivCyqqgonp4vzsikJ5xflZgfJJzdIKhMJpzA3SCh+tOHc4uLjOS0S4+NGW98wzT2DNHcP0dQ9SHNP8NzUPURLT/jcGz+JFedlUVWUS0VhDhWFOVQVBc+VhblUFuWQlx2cMsrOzGDcOOtIpXtg5HSCmUg4/SN0D06fZErysigryKGsIJvS/GzKw+my/GxK8rMpiUk0UxORtypzLr35eE5pICPj9DWvy2oTLzc2brT1DdHcPUTzRMLqGaK9b5i2vmHaeoc43t7Py8c66egfjpvIztq2gib55YU5lOVnU1qQQ31lweR0WX425YXZlOXnUBomnrKCHErzs725vXPuvHhySjOZGWJJcR5LivOAs69/xRofN7oGRmjrG2ZodIyRMWN0bByJM06dFeb4aTHn3MKK5Gk9SS3A0Tm+TRWwIL1SzAOPNTk81uRYTLHC4or3Qot1lZlVx5sRyeQ0HyRtTnY/fvPFY00OjzU5FlOssLji9VhP83awzjnnIseTk3POuchJ5+QUtw+/iPJYk8NjTY7FFCssrng91lDaXnNyzjm3eKXzkZNzzrlFypOTc865yEm75CTpFkn7JB2QdG+q44klaaWkX0jaI2m3pD8Jyz8v6aSkbeHj1lTHCiDpiKSdYUybw7IKST+VtD98Lp/pfRaCpPUx9bdNUrekT0albiU9JKlZ0q6YsoR1Kelz4T68T9LbIxDr30raK2mHpMcllYXl9ZIGYur3gQjEmvAzj2C9fj8mziOStoXlqa7XRN9VC7fPmlnaPAg6pj0IrAFygO3AhlTHFRPfMuCacLoYeBXYAHwe+Eyq44sT7xGgakrZl4F7w+l7gS+lOs4E+0EjsCoqdQvcBFwD7JqpLsN9YjuQC6wO9+nMFMf6NiArnP5STKz1sctFpF7jfuZRrNcp878C3B+Rek30XbVg+2y6HTlNDilvZsPAxJDykWBmDWa2NZzuAfYQjDS8mNwBPBxOPwy8K4WxJPIW4KCZzbWXkXljwQCb7VOKE9XlHcCjZjZkZoeBAwT79oKIF6uZPWVmE2OUvEAwrlvKJajXRCJXrxMU9A/2+8AjCxXPdKb5rlqwfTbdklO8IeUj+eUvqR64GngxLPpEeMrkoaicKiMY8fgpSVsk3ROW1ZhZAwQ7MLAkZdEldidn/pNHsW4hcV1GfT/+MPCTmNerJb0s6VeS3piqoKaI95lHuV7fCDSZ2f6YskjU65TvqgXbZ9MtOc16SPlUklQE/CfwSTPrBv4RWAtcBTQQHN5HwRvM7BpgE/BxSTelOqCZSMoBbgf+IyyKat1OJ7L7saT7gFHgu2FRA1BnZlcDnwK+J6kkVfGFEn3mka1X4L2c+YMqEvUa57sq4aJxyuZUt+mWnGY1pHwqScom+LC/a2aPAZhZk5mNmdk48C0W8FTDdMzsVPjcDDxOEFeTpGUA4XNz6iKMaxOw1cyaILp1G0pUl5HcjyXdBbwTeL+FFxrC0zht4fQWgmsNF6cuymk/86jWaxbwHuD7E2VRqNd431Us4D6bbslpckj58Bf0ncATKY5pUnhe+dvAHjP7u5jyZTGLvRvYNXXdhSapUFLxxDTBBfFdBPV5V7jYXcAPUxNhQmf8Ao1i3cZIVJdPAHdKypW0GlgHvJSC+CZJugX4LHC7mfXHlFdLygyn1xDEeig1UU7GlOgzj1y9ht4K7DWzExMFqa7XRN9VLOQ+m6rWIMl6ALcStCw5CNyX6nimxHYjwaHuDmBb+LgV+FdgZ1j+BLAsArGuIWh9sx3YPVGXQCXwM2B/+FyR6lhjYi4A2oDSmLJI1C1BwmwARgh+Zd49XV0C94X78D5gUwRiPUBwTWFiv30gXPZ3wv1jO7AVuC0CsSb8zKNWr2H5d4CPTlk21fWa6LtqwfZZ777IOedc5KTbaT3nnHNpwJOTc+dB0o2SnpPUJald0rOSXhfO+6CkX5/De9VLsvDC+PnEcrGkH0pqCWN5UtL683kv56LCk5Nz5yhs0vtfwD8AFQT3c/wFMJSikMoIrq2sB2oILkRHraGKc+fErzk5d44kbQT+PzMrizPvUuBlIBsYAEbNrEzSO4C/JLj/pgv4tpl9PlznGEEz3L7wbX7bzJ6X9GHgT4GlBAnnHptFrxeSKggahlRZ2BzZucXGj5ycO3evAmOSHpa0KbbXCTPbA3wUeN7MimISWB/wAYKjnHcAH5M00fXLxM3NZeE6z4fz/ozg/pdq4Blm37XNTUCjJya3mHlycu4cWXCn/ERT228BLZKekFQzzTq/NLOdZjZuZjsIEs2bptnMR4C/NrM9FvRp91fAVZJWTRebpBXANwh6FXBu0fLk5Nx5CJPGB81sBXA5UAt8NdHykq4LhyBokdRFcHRVNc0mVgFfk9QpqZOgw1AxTX9lkqqBp4BvmlkkOhB17nx5cnJujsxsL8GNlJdPFMVZ7HsEjRZWmlkp8ACn+yOLt/xx4CNmVhbzyDez5+LFEJ5afAp4wsy+eP5/jXPR4MnJuXMk6RJJnw5PoSFpJUG3SS+EizQBK8IutCYUA+1mNijpWuB9MfNagHGCXjkmPAB8TtJl4TZKJf1egnhKgCeBZ80sUgNsOne+PDk5d+56gOuAFyX1ESSlXcCnw/k/J+h6plFSa1j2R8AXJPUA9wP/PvFmFvRV90Xg2fA03vVm9jjBoH6PSuoO339TgnjeDbwO+JCk3phH3Tz+zc4tKG9K7pxzLnL8yMk551zkeHJyzjkXOZ6cnHPORY4nJ+ecc5FzXr0gJ1tVVZXV19enOgznnHNJtGXLllYzq443L5LJqb6+ns2bN6c6DOecc0kkKWFHxn5azznnXOR4cnLOORc5npycc85Fjicn55xzkePJyTnnXOR4cnLOORc5npycc85Fjicn55xzkePJyTnnXOR4cnLOORc5npycc85FzpySk6SHJDVL2pVgfrmkxyXtkPSSpMvnsj3nnHMXhrkeOX0HuGWa+X8GbDOzK4EPAF+b4/acc85dAOaUnMzsaaB9mkU2AD8Ll90L1Euqmcs2nXPOpb9kX3PaDrwHQNK1wCpgRZK36ZxzbpFLdnL6G6Bc0jbgj4GXgdF4C0q6R9JmSZtbWlqSHJZzzrkoS+pgg2bWDXwIQJKAw+Ej3rIPAg8CbNy40eay3b/40W5eOdU9l7dwzjk3jQ21Jfz5bZcl7f2TeuQkqUxSTvjyD4Gnw4TlnHPOJTSnIydJjwA3A1WSTgB/DmQDmNkDwKXAv0gaA14B7p5TtLOUzGzunHMu+eaUnMzsvTPMfx5YN5dtOOecu/B4DxHOOecix5OTc865yPHk5JxzLnI8OTnnnIscT07OOecix5OTc865yPHk5JxzLnI8OTnnnIscT07OOeciJ9kj4ZZK+pGk7ZJ2S/rQXLbnnHPuwpDskXA/DrxiZq8h6IPvKzEdwTrnnHNxJXskXAOKw+EyisJl447n5Jxzzk1I6nhOwNeBJ4BTQDHwP8xsPMnbdM45t8glu0HE24FtQC1wFfB1SSXxFvSRcJ1zzk1IdnL6EPCYBQ4QjIJ7SbwFzexBM9toZhurq6uTHJZzzrkoS/ZpvWPAW4BnJNUA64FDM620ZcuWVklH57jtKqB1ju+xUDzW5PBYk2MxxQqLK94LLdZViWbIzM77XWNHwgWamDISrqRaghZ9ywABf2Nm/3beGzy32Dab2caF2NZceazJ4bEmx2KKFRZXvB7rackeCfcU8La5bMM559yFx3uIcM45FznpnJweTHUA58BjTQ6PNTkWU6ywuOL1WENzuubknHPOJUM6Hzk555xbpDw5Oeeci5y0S06SbpG0T9IBSfemOp5YklZK+oWkPWEv7X8Sln9e0klJ28LHramOFUDSEUk7w5g2h2UVkn4qaX/4XJ7qOAEkrY+pv22SuiV9Mip1G68H/+nqUtLnwn14n6S3RyDWv5W0V9IOSY9LKgvL6yUNxNTvAxGINeFnHsF6/X5MnEckbQvLU12vib6rFm6fNbO0eQCZwEFgDZADbAc2pDqumPiWAdeE08XAq8AG4PPAZ1IdX5x4jwBVU8q+DNwbTt8LfCnVcSbYDxoJbvCLRN0CNwHXALtmqstwn9gO5AKrw306M8Wxvg3ICqe/FBNrfexyEanXuJ95FOt1yvyvAPdHpF4TfVct2D6bbkdO1wIHzOyQmQ0DjwJ3pDimSWbWYGZbw+keYA+wPLVRnbM7gIfD6YeBd6UwlkTeAhw0s7n2MjJvLH4P/onq8g7gUTMbMrPDwAGCfXtBxIvVzJ4ys4kRBV4AVixUPNNJUK+JRK5eJ4QjN/w+8MhCxTOdab6rFmyfTbfktBw4HvP6BBH98pdUD1wNvBgWfSI8ZfJQVE6VEQx58pSkLZLuCctqzKwBgh0YWJKy6BK7kzP/yaNYt5C4LqO+H38Y+EnM69WSXpb0K0lvTFVQU8T7zKNcr28Emsxsf0xZJOp1ynfVgu2z6ZacFKcscm3lJRUB/wl80sy6gX8E1hL03N5AcHgfBW8ws2uATcDHJd2U6oBmomAwy9uB/wiLolq304nsfizpPoIx2b4bFjUAdWZ2NfAp4HtKMPLAAkr0mUe2XoH3cuYPqkjUa5zvqoSLximbU92mW3I6AayMeb2CYCypyJCUTfBhf9fMHgMwsyYzG7NgrKtvsYCnGqZjQfdTmFkz8DhBXE2SlgGEz82pizCuTcBWM2uC6NZtKFFdRnI/lnQX8E7g/RZeaAhP47SF01sIrjVcnLoop/3Mo1qvWcB7gO9PlEWhXuN9V7GA+2y6JaffAOskrQ5/Qd9JMNhhJITnlb8N7DGzv4spXxaz2LuBXVPXXWiSCiUVT0wTXBDfRVCfd4WL3QX8MDURJnTGL9Ao1m2MRHX5BHCnpFxJq4F1wEspiG+SpFuAzwK3m1l/THm1pMxweg1BrDOOPJBM03zmkavX0FuBvWZ2YqIg1fWa6LuKhdxnU9UaJFkP4FaCliUHgftSHc+U2G4kONTdQTAI47Yw3n8FdoblTwDLIhDrGoLWN9uB3RN1CVQCPwP2h88VqY41JuYCoA0ojSmLRN0SJMwGYITgV+bd09UlcF+4D+8DNkUg1gME1xQm9tsHwmV/J9w/tgNbgdsiEGvCzzxq9RqWfwf46JRlU12vib6rFmyf9e6LnHPORU66ndZz7oIi6TuS/jLVcTg33zw5OTcDSe+TtFlSr6QGST+RdOM8vG/SE0sY+1FJfZJ+IKkimdtzbr54cnJuGpI+BXwV+CugBqgDvskC3NwdtuKay/qXAf8E/E+C2PsJYncu8jw5OZeApFLgC8DHzewxM+szsxEz+5GZ/Wm4TK6kr0o6FT6+Kik3nHezpBOSPh32qdYg6UPhvHuA9wP/Ozwi+1FYfkTSZyXtAPokZUm6VNIvJXWG/ZzdPss/4f3Aj8zsaTPrBf5f4D0TrTCdizJPTs4ldgOQR3CPVyL3AdcT3PD5GoJ7av5PzPylQCnB3fJ3A9+QVG5mDxLcyPplMysys9ti1nkv8A6gjODmxh8BTxHcjf/HwHclrZ9F/JcRtPYCwMwOAsOk+D4k52bDk5NziVUCrXa6T7l43g98wcyazawF+AuC02gTRsL5I2b2Y6AXmCmx/L2ZHTezAYLEVwT8jZkNm9nPgf8iSGAzKQK6ppR1EXTk6VykeXJyLrE2oGqGaz+1QGwHs0fDssn3mJLc+gmSxnRi+yirBY5b0NtB7DZm029ZLzC1y5sSoGcW6zqXUp6cnEvseWCQ6XteP0UwNMeEOmbfbUuimwxjy08BKyXF/q/WASdn8f67CU41ApM9DeQS3KTuXKR5cnIuATPrAu4nuE70LkkFkrIlbZL05XCxR4D/E3Y3UxUu/2+z3EQTQU8c03kR6CNoOJEt6WbgNoLhYGbyXeA2SW8Mu6D6AvCYBUMgOBdpnpycm4YF/Yp9iqCRQwvBKbdPAD8IF/lLYDNBNy87Cbqame29S98GNoSt8H4QbwELxiW7naBD21aCpuAfMLO9s4h9N/BRgiTVTHCt6Y9mGZtzKeXdFznnnIscP3JyzjkXOZ6cnHPORY4nJ+ecc5Ezp+Qk6aGwW5a4A7iF3bd0SdoWPu6fy/acc85dGObUsSTBIFlfB/5lmmWeMbN3nsubVlVVWX19/RzCcs45F3VbtmxpNbPqePPmlJzM7GlJ9XN5j3jq6+vZvHnzfL+tcxeGn/8cvvhFGB+feVnnztfdd8Mf/MGc3kLS0UTz5nrkNBs3SNpOcKf7Z8J7L84S9tJ8D0BdXd0ChOVcmnrkEXj2WbjuulRH4tJZkm9DSnZy2gqsMrNeSbcS3Li4Lt6CYS/NDwJs3LjRb75y7nydOAGXXw6/+lWqI3HuvCW1tZ6ZdYfjyBD2yJwddvHinEuW48dh5cpUR+HcnCQ1OUlaKknh9LXh9tqSuU3nLngnTsCKFamOwrk5mesw0I8ANxMMK3AC+HMgG8DMHgB+F/iYpFFgALjTvL8k55Knpwe6uvzIyS16c22tN+2AZ2b2dYKm5s65hXDiRPDsR05ukfMeIpxLJ56cXJrw5ORcOjkeDqLrp/XcIufJybl0MnHkVFs7/XLORZwnJ+fSyYkTUFMDubmpjsS5OfHk5Fw6OX7crze5tODJybl04vc4uTThycm5dHLihDeGcGnBk5Nz6aK3Fzo7/cjJpYVkDzYoSX8v6YCkHZKumcv2nHPT8HucXBqZ65HTd4Bbppm/iaAX8nUEw2H84xy355xLxJOTSyNzSk5m9jTQPs0idwD/YoEXgDJJy+ayTedcAg0NwbPf4+TSQLKvOS0Hjse8PhGWnUXSPZI2S9rc0tKS5LCcS0ONjcHz0qWpjcO5eZDs5KQ4ZXF7JTezB81so5ltrK6OO6S8c246jY1QUABFRamOxLk5S3ZyOgHEtmtdQTBcu3NuvjU2BkdNiveb0LnFJdnJ6QngA2GrveuBLjNrSPI2nbswTSQn59JAsgcb/DFwK3AA6Ac+NJftOeem0dgIl1yS6iicmxfJHmzQgI/PZRvOuVlqbISbb051FM7NC+8hwrl0MDQE7e1+Ws+lDU9OzqWD5ubg2ZOTSxOenJxLB36Pk0sznpycSweenFya8eTkXDrw5OTSjCcn59LBRHJasiS1cTg3Tzw5OZcOGhuhogJyc1MdiXPzwpOTc+nAe4dwaWaugw3eImlfOJjgvXHm3yypS9K28HH/XLbnnEvAk5NLM+fdQ4SkTOAbwG8TdPD6G0lPmNkrUxZ9xszeOYcYnXMzaWyE669PdRTOzZu5dF90LXDAzA4BSHqUYHDBqclp4d19N/T1pToK5xbO8ePw7nenOgrn5s1cklO8gQSvi7PcDZK2EwyV8Rkz2x3vzSTdQzCUO3V1dXMIC9i1C7q65vYezi0m69bB296W6iicmzdzSU6zGUhwK7DKzHol3Qr8AFgX783M7EHgQYCNGzfGHZBw1l58cU6rO+ecS625JKcZBxI0s+6Y6R9L+qakKjNrne6Nt2zZ0irp6BxiA6gCpt1OhHisybOY4vVYk2cxxXshxboq0Yy5JKffAOskrQZOAncC74tdQNJSoMnMTNK1BK0D22Z6YzOb8zjtkjab2ca5vs9C8FiTZzHF67Emz2KK12MNnHdyMrNRSZ8AngQygYfMbLekj4bzHwB+F/iYpFFgALgzHOPJOeecS2iugw3+mGC029iyB2Kmvw58fS7bcM45d+FJ5x4iHkx1AOfAY02exRSvx5o8iylejxWQn2VzzjkXNel85OScc26R8uTknHMuctIuOc3UGW2qSVop6ReS9kjaLelPwvLPSzoZ00nuramOFUDSEUk7w5g2h2UVkn4qaX/4XB6BONfH1N02Sd2SPhmlepX0kKRmSbtiyhLWpaTPhfvxPklvj0Csfytpr6Qdkh6XVBaW10saiKnjBxK/84LFmvBzT2W9ThPv92NiPSJpW1ie6rpN9H2V/P3WzNLmQdCk/SCwBsgBtgMbUh3XlBiXAdeE08XAq8AG4PME3TulPMYp8R4BqqaUfRm4N5y+F/hSquOMsx80EtzgF5l6BW4CrgF2zVSX4T6xHcgFVof7dWaKY30bkBVOfykm1vrY5SJSr3E/91TXa6J4p8z/CnB/ROo20fdV0vfbdDtymuyM1syGgYnOaCPDzBrMbGs43QPsIeincDG5A3g4nH4YeFcKY4nnLcBBM5trLyPzysyeBtqnFCeqyzuAR81syMwOAwcI9u8FES9WM3vKzEbDly8Q9AqTcgnqNZGU1itMH68kAb8PPLKQMSUyzfdV0vfbdEtO8TqjjewXv6R64GpgojPAT4SnTB6KwqmykAFPSdoSds4LUGNmDRDsvEDUxga/kzP/uaNYrxMS1WXU9+UPAz+Jeb1a0suSfiXpjakKaop4n3vU6/WNBL3q7I8pi0TdTvm+Svp+m27JaTad0UaCpCLgP4FPWtAH4T8Ca4GrgAaCQ/soeIOZXQNsAj4u6aZUBzQdSTnA7cB/hEVRrdeZRHZflnQfMAp8NyxqAOrM7GrgU8D3JJWkKr5Qos89svUaei9n/rCKRN3G+b5KuGicsvOq33RLTjN2RhsFkrIJPujvmtljAGbWZGZjZjYOfIsFPtWQiJmdCp+bgccJ4mqStAwgfG5OXYRn2QRsNbMmiG69xkhUl5HclyXdBbwTeL+FFxnCUzht4fQWgusMF6cuymk/90jWK4CkLOA9wPcnyqJQt/G+r1iA/TbdktNkZ7ThL+g7gSdSHNMZwnPK3wb2mNnfxZQvi1ns3cCuqesuNEmFkoonpgkuiO8iqNO7wsXuAn6YmgjjOuOXZxTrdYpEdfkEcKekXAWdK68DXkpBfJMk3QJ8FrjdzPpjyqsVjIyNpDUEsR5KTZSTMSX63CNXrzHeCuw1sxMTBamu20TfVyzEfpuqViBJbF1yK0GLkoPAfamOJ058NxIc5u4AtoWPW4F/BXaG5U8AyyIQ6xqCljfbgd0T9QlUAj8D9ofPFamONYyrgKDX+9KYssjUK0HSbABGCH5h3j1dXQL3hfvxPmBTBGI9QHA9YWK/fSBc9nfC/WM7wRhut0Ug1oSfeyrrNVG8Yfl3gI9OWTbVdZvo+yrp+613X+Sccy5y0u20nnPOuTTgyck551zkeHJyzjkXOZ6cnHPORY4nJ+ecc5Hjyck551zkeHJyzjkXOf8/JpBSJ20pxuQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "bocop.run(problem_path)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading solution:  ./problem.sol\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl8AAAEMCAYAAAD6TdPZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxV1bn/8c9KSAiTzDKFSRHDGCQooKg4D9WCirU4V8S5tna8vbetbe1wb6vW2loHFFsVh1at1qlqf+IsCAkJEaTIZAggmIEhZD5n/f7IoGefQKZzzjr72c/79eIlCSf7PF8XKzxZe++1jbUWpZRSSimVGCmuC1BKKaWUChJtvpRSSimlEkibL6WUUkqpBNLmSymllFIqgbT5UkoppZRKIG2+lFJKKaUSqIvrAtqqX79+9rDDDnNdhnKotraW9PR012WoA8jNzS2x1g50XUeQDBgwwI4aNeqgr5E+b6TnA/kZpeY72PdE3zRfI0eOZOXKla7LUA6Vl5fTt29f12WoAzDGfOq6hqAZNWpUq98Xpc8b6flAfkap+Q72PdE3px3D4bDrEpRjlZWVrktQynekzxvp+UB+Run5WuKb5kuplBT966pUe0mfN9LzgfyM0vO1xDeJjTGuS1COpaWluS5BKd+RPm+k5wP5GaXna4lvrvnS046qoqKCAQMGOHv/uro6iouLqa6udlZDMsjIyCAzMzOQ3zD9yPW8iTfp+UB+Run5WuKb5qtLF9+UquLE9eQsLi6mV69ejBo1KrArsdZaSktLKS4uZvTo0a7LUW3get7Em/R8ID+j9Hwt8c1px9raWtclKMeKi4udvn91dTX9+/cPbOMFDaf/+/fvH/jVPz9xPW/iTXo+kJ9Rer6W+GY5qWvXrq5LUI6NGTPGdQmBbrya6P+DlhljlgFft9ZuMcYMA5631k5zXVcyzJt4kp4PWs5YFwpTWiFjUaLnoZl8tsc/P9D165FOepfOrV35pvnSn7TVmjVryM7Odl2GUlFMQ0c6Amja12cyUOiuoi9InzfS80HLGb/5+Cr+teYzRxUF2ws3zWJSZu9OHcM3zVe3bt1cl6Ack/4Ntq2qqqo488wzeeONN9i1axe33HILGzZsYO/evQwbNoylS5d26vi1tbWceuqpvPHGG3qtZduNATZba23jx3Frvowx1wDXAAwdOpSSkhLq6uoIh8N0796dsrIyhg4dyubNm8nKyiIUCgGQm5tLTk4OeXl5TJkyhXXr1jF69Gi2b99Ov379qKysJCUlhbS0tOYLoIuLixkzZkzzP/5Nx2j6b2FhIWPHjqWoqIhBgwaxe/duMjIygIYfmPv06cPOnTsZMWIE69evZ9KkSVHHKCgoYMKECWzYsIHMzExKSkro2bPnQTPl5+czderUwGQaO3YsGzdujMi0tXQfo/p25cqZIynfXc7AgQPZsX07mcOH8+mnnzJy5Mjm/27bto1Bhw6irKyUXoccQlVVFWldGm6Yqauvo1u3buzbu5d+/fqzc9dOhg0bFnWM4q1bGTJ0KJ9//jl9+/SlomIfXTMyCIVC2LAlPT2d/ZX76d27N6UlJQweMoStRVsZMXJE8zGKPi1i+IjhfLZjB/0HDGDPnj306N6D2tpaTIohNTWVmupqevbsldSZhvbJaNPfvYOy1vri17hx46wKtpUrVzp9/7Vr1zp9/yZ/+tOf7F133WWttfaUU06xTz75ZPOfrV69Oibv8bOf/cw+9thjB/zzlv5fACttEnyvcPELOB+470sf/x04FRgMvAX8APgrcC3wNDCx8XW3An8Afg4MBB4GMoHFQFpr75uTk3PAMWriet7Em/R81rac8Zy737FXLF7uoJrYkzqGB/ue6JsL7rt37+66BOVYTk6O6xKSwpIlS5gzZw6hUIg333yTE088sfnPJk2aFJP3mDt3LkuWLInJsQKiH1AFYIwZB3yFhpWvo4BnrbW/BXoDi2hozEY2XheWBuwGZlhrPweKgDuAm621dbEoTPq8kZ4PWs4YCltShVx/GYQx9PJN8xXExw+oSLm5ua5LcK62tpZNmzYxatQoUlNTOfXUU8nOzubaa6/lvffei9n7TJw4kRUrVsTseAHwKnCKMeZvwIVAqbV2JzAFeNUYk9b4uTAwkYbG7Dbg/2hYEdtmjOkJHAbUW2srYlWY9HkjPR+0nDFsLSkpMpqvIIyhl2+aL135UkH86cirpKSEPn36NH/8yiuv8Mwzz9C7d2/OPPNMnnvuuZi8T2pqKunp6ezbty8mx5POWrvVWjvZWvs1a+0vrLXDG/9oDLAemAB83Pi5UdbaImAN8D3gO8Aq4G7gx0C+MWZ2rGqTPm+k5wNd+ZLIN1fTVlVVuS5BOVZYWBiz02qx8Oabb/LWW281f7xw4UIAFi1a1Py5E088kdmzZ3PHHXdQUdGwmDFkyBCuueYaXnjhBfLy8ppf+53vfIdevXod9D27desWceevMYZZs2Yxa9YsysvLWb16NXPnzqWwsJBXX32V733ve9xwww3cfvvt7f4BpqampvlCY9Ux1toFjb/Nb/yFtfayxv/ecYAv+10sa0i2eRNr0vNByxlD1pIqZOUrCGPo5ZvmS/8RUGPHjnVdQoTZs2cze/bsqM/feuutUZ/77ne/G/W5c889l3PPPbdd79m3b19CoRDV1dW89dZbnHTSSaSnp7Nr1y7effddFi9eDMCqVaua7w6trKxsd+NVWlrKwIED9RFCAiTbvIk16fmg5YzhsJzTjkEYQy/fnHbUHe5VUVGR6xKSwumnn867777L008/zbhx48jOzuacc87htttuY+bMmUDDT5KTJ09m7969HdoUdenSpZx99tmxLl05IH3eSM8HLWcMWUuqjN4rEGPo5ZuVL91vSA0aNMh1CUnhpptu4s477+TRRx894GuysrK4/fbb6dKlC1lZWTz00EOEQiFWrVrFtGnTGDx4MJMmTeKll14iPT2dESNGUFZWRmlpKTfccAOPP/44v/nNbxKYSsWL9HkjPR+0nDEcRszKVxDG0Ms3HU3TRoEquHbv3s0hhxziugznjjrqKE466SRCoRCpqaktvmbBggURH99xxx1897vf5Xe/+x1nnXUW77zzDuvXrycjI4Nt27bxta99jTPOOIP33nuP2tpa5s6dy5FHHpmIOCrOpM8b6fmg5YySLrgPwhh6+ea0Y0qKb0pVcaLX/X3hqquuOmDj1ZLU1FRqa2vp2bMnH330EZMmTWreSqJfv37cddddfP/732flypWkp6dz+eWXx6t0lWDS5430fNByRkkX3AdhDL18s/KllOq4b3/72wBcf/31zZ8bP368q3KUUp0k6YL7IPLNclI4HHZdgnJMH66uVPtJnzfS80HLGRsuuJfRfAVhDL1803y15xSLkunLm4sqpdpG+ryRng9azhgKyzntGIQx9PJN81VfX++6BOXYzp07XZeglO9InzfS80HLGcNhS4qQla8gjKGXb5qv9PR01yUox0aMGOG6BKV8R/q8kZ4PWs7YcMG9g2LiIAhj6OWboQviOWEVaf369a5LwFrrugTn9P+BvyTDvIkn6fmg5YyS9vkKwhh6+ab56tatm+sSlGOun/2VkZFBaWlpoJsPay2lpaWBvDXcr1zPm3iTng9azijpgvsgjKGXb7aaqKysdF2Cciw3N5ecnBxn75+ZmUlxcTGff/65sxqSQUZGBpmZma7LUG3ket7Em/R80HJGSRfcB2EMvXzTfLX3wcBKHteTMy0tjdGjRzutQan2cj1v4k16PojOGA43rL5LueA+CGPo5ZvTjrrypXJzc12XoJTvSJ830vNBdMZQ46UPkla+gsY3zZeufKkg/nSkVGdJnzfS80F0xlBYVvMVhDH08k3zVVVV5boE5VhBQYHrEpTyHenzRno+iM4YtrJOOwZhDL1803zp3VVqwoQJrktQynekzxvp+SA64xcrXy6qib0gjKGXb4aupqbGdQnKsQ0bNrguQSnfkT5vpOeD6IxNjzqWsvIVhDH08k3zpTvcK93eQKn2kz5vpOeD6IzSLrgPwhh6+ab50mc7qpKSEtclKOU70ueN9HwQnVHaBfdBGEMv3zRfKSm+KVXFSc+ePV2XoJTvSJ830vNBdEZpF9wHYQy9fNPRBPmRLqpBXV2d6xKU8h3p80Z6PojOKG3lKwhj6OWb5kupcNNVpkqpNpM+b6Tng+iMzc2XkJWvIIyhl2+aLz3tqHSjXaXaT/q8kZ4PojM2n3YUsvIVhDH08k1Hoxfcq7KyMtclKOU70ueN9HwQnVHaPl9BGEMv3wxdWlqa6xKUY0OHDnVdglK+I33eSM8H0RmlXXAfhDH08k3zVVtb67oE5djmzZtdl6CU70ifN9LzQXTGUOMlUlIuuA/CGHo5bb6MMT8yxtxujMlq7bX6eCGVldXqXxOllIf0eSM9H0RnlHbBfRDG0MtZ82WMmQ7MBzZba9e19vrKysr4F5UApRU1FGzdTWmFPi6pvfLz812XoJTvSJ830vNBdEZpF9wHYQy9ujh87/XAm9bae9ryYgl3Qzyfv40fPrOatJQU6sJhfnvBZL46ZZjrsnxj6tSprktQynekzxvp+SA6Y1PzJWXlKwhj6JWQlS9jzFJjzGmNv/+lMeZuYApQ0NZj+H3lq7Sihh8+s5rqujD7auqprgvzg2dW6wpYO+Tm5rouQSnfkT5vpOeD6IzSNlkNwhh6Jeq0463A/xhjLgGOAm4BsoGDrjUaY64xxqw0xqzcs2cPJSUl7Nixg23btlFeXs7GjRupqqpi7dq1hMNh8vLygC8GMi8vj3A4zNq1a6mqqmLjxo2Ul5ezbds2duzYQUlJCVu2bKGiooJ169ZRX19PQUFBxDGa/ltYWEhNTQ2ffPIJe/fupaioiF27drFr1y6KiorYu3cvn3zyCTU1NRQWFkYdo7i8ihQ8u/Rby0ebd/g2E0BBQQH19fWsW7eOiooKtmzZErdxOvLII8VlkjROKjnl5OS4LiGupOeD6IzSTjsGYQy9TKIe22OMeQvoCcy21u5r79ePHz/erl27NvaFJUhpRQ3H/d8bVNdF7uT7vdPHcuNJYzBClo/jKS8vL5DL035hjMm11k5zXUeQTJs2za5cufKgr5E+b6Tng+iMH24u42v3f8BjC6Yz64gBDiuLDaljeLDviYk67TgJGALUdKTxAv9f89W/Z1d+e8FkMtJS6NW1C127GCYNO4TbX1vPD55eTU19yHWJSW/KlCmuS1DKd6TPG+n5IDpj02lHKQ9+CcIYesV96IwxQ4AlwBxgvzHmjMbPzzXGLDLGPG+MOb2141RXV8e50vj76pRhvPfDk3ns6um8/1+n8PyNs7j5lCP4e24xFy9aTole/3VQ69a1elOsUspD+ryRng+iM0q74D4IY+gV1+bLGNMdeBb4rrX2Y+A24GcA1trnrLULgSuBi1o7Vnp6evwKTaD+PbuSPbwP/Xt2JSXF8J3TxnLPxVP5aNse5t7zHut3dmhhMBBGjx7tugSlfEf6vJGeD6IzSrvgPghj6BXX5staW2mtnWmtfb3x47ettTM9L/sx0Op2E3V1dfEoMSl8ZfIQ/nbtTGrqw5z/5/dZ+p9drktKStu3b3ddglK+I33eSM8H0RlDwi64D8IYerncZNUYY/4PeMVam9fa67t0cbklWfxlD+/D8zcex4h+3VnwlxU8/N5mEnUzhF/069fPdQlK+Y70eSM9H0RnDAvb4T4IY+jl8nK9bwKnAvOMMde19uJwONzaS3xvaJ9u/P26mZw6bhA/f2EtP37uI+pC8nO3ld/3elPKBenzRno+iM4o7bRjEMbQy9lykrX2buBuV++frHp07cJ9l+bw21f/w31vbeTT0kruuXgqvbunuS7NuRQpt/YolUDS5430fBCdsXmfLyErX0EYQy/fJA7SPlgpKYb/OiuL382bzPLNpZx373tsKdnvuizn0tK0AVWqvaTPG+n5IDpj0wkRKStfQRhDL980X0E47eh14bThPLZgOuX7a5n75/f4cHOZ65KcqqiocF2CUr4jfd5IzwfRGZsuuE/1zb/gBxeEMfTyzdBJv+D+QKYf1p/nbjyOfj3SufSh5by4Onh3hTQZMMD/OzkrlWjS5430fBCdsemCeymnHYMwhl6+ab5qa2tdl+DMyP49ePb6Y8nO7M1Nj6/igbc3BvJOyOLiYtclKOU70ueN9HwQnVHaBfdBGEMv3zRfXbt2dV2CU326p/Pogul8ZfIQfv3yOm7955rmCRgUY8aMcV2CUr4jfd5IzwfRGUPCLrgPwhh6+ab5kvB4oc7KSEvlj18/imtPOIxHPviUax9dSWVtveuyEmbNmjWuS1DKd6TPG+n5IDpjWNjKVxDG0Ms3zVe3bt1cl5AUUlIMPzp7HL+YM4E31u1i/gPL+HxfMJ4JmZ2d7boEpXxH+ryRng+iM35xwb2M5isIY+jlm+YriJuwHczlM0dx/2XT+M/OfZx/73ts/Fz+3SK5ubmuS1DKd6TPG+n5IDqjtAvugzCGXr5pvrp37+66hKRz2vhBPHnNTKpqQ1xw7/us2CJ7K4qcnBzXJSjlO9LnjfR8EJ1R2gX3QRhDL980X7ry1bIpw/vw7PXH0a97Opc8uJxXCne4LilugvjTkVKdJX3eSM8H0RlDjfdaSXm2YxDG0Ms3zZeufB3YiP7deeb6Y5k0rDc3PJ7HY8s+dV1SXATxpyOlOkv6vJGeD6IzNp929M2/4AcXhDH08s3QVVVVuS4hqfXtkc5jC6Zz8pGH8uPnPuL3r68XtxdYYWGh6xKU8h3p80Z6PojOKO2C+yCMoZdvmq+MjAzXJSS9bump3H9ZDhfmZPKH//cJ//PcR6L2Ahs7dqzrEpTyHenzRno+iM4YEnbBfRDG0Ms3zVeQd7hvjy6pKfx23mRumH04jy8v4sYleVTXhVyXFRNFRUWuS1DKd6TPG+n5IDqjtH2+gjCGXr5pvoL6bMeOMMbwgzOz+Ok54/nXms+4YvGH7K2uc11Wpw0aNMh1CUr5jvR5Iz0fRGdsPu0oZOUrCGPo5ZvmKxSSsXqTSFfNGs0fvj6FvKJyLrp/Gbv2+vspAbt373ZdglK+I33eSM8H0Rm/uOBeRvMVhDH08k3zlSLlto4EmzNlGIuvPJpPS/dz/r3vs7lkv+uSOkyv+1Oq/aTPG+n5IDpjyFoxpxwhGGPopR1NABx/xECevGYGlbUh5t37PquLg/dThlJKSREKyznlGFS+ab7C4bDrEnxtcmYfnrn+WLqlp3LxouUs21TquqR204erK9V+0ueN9HwQnTFsrZg9viAYY+jlm+FLTU11XYLvjR7Qg6evO5YhvTO4YvGHvLFup+uS2qVPnz6uS1DKd6TPG+n5IDpjKGxFrXwFYQy9fNN81dfXuy5BhMG9M3jq2pmMHdSLax7J5Z8F212X1GY7d/qrWVQqGUifN9LzQXTGUNiKudgegjGGXr5pvtLT012XIEa/Huk8vnA6U0f25VtPruLx5f7YY2XEiBGuS1DKd6TPG+n5IDpjWNgF90EYQy/fNF9BPCccT70y0njkqmOYPXYg//2PQu5/a6Prklq1fv161yUo5TvS5430fBCdUdppxyCMoZdvmq9u3bq5LkGcjLRU7r9sGudMHsJvXlnH715dl9TPg5w0aZLrEpTyHenzRno+iM7YcMG9nOYrCGPo5Zvmq7Ky0nUJIqV3SeEPXz+K+ccM556lG7n1n2uaN/BLNrm5ua5LUMp3pM8b6fkgOqO0la8gjKGXb57Z0717d9cliJWaYvj1eZPolZHGA29vYl91Pb+bN5kuqcnVm+fk5LguQSnfkT5vpOeD6IyhsJznOkIwxtAruf51PQhd+YovYww/OiuL759xJP9YtY3rk/CB3EH86UipzpI+b6Tng+iM0vb5CsIYevlm+HTlK/6MMdx40hh+/tUJvL52Jwv+uoLK2uTZ4iOIPx0p1VnS5430fNDSypes045BGEMv3zRfVVVVrksIjCuOHcUdF2bzwcZSrlj8Ifuq61yXBEBBQYHrEpTyHenzRno+iM4YEnbBfRDG0Ms3zVcQH7zp0gU5mfxx/lRWFe3m0oc+ZE+l+wZswoQJrktQynekzxvp+SA6Y1jYylcQxtDLN81XTU2N6xIC5yuTh/DnS6by8fa9XPzgMsr21zqtZ8OGDU7fXyk/kj5vpOeD6IyhsKxNVoMwhl6+ab50h3s3Tp8wmAcuz2HDrgrmP7CMz/e5a4IzMzOdvbdSfiV93kjPB9EZw9aSImjlKwhj6OWb5kuf7ejO7CMPZfGVR1NUVslFD3zAZ3vcPG2gpKTEyfsq5WfS5430fBCdUdrKVxDG0Ms3zVeKpPtqfei4MQP461XHsGtvDRc98AHF5Ynf+qNnz54Jf0+l/E76vJGeD6IzhiyiLrgPwhh6+aajSebH3gTFMaP78eiCYyjbX8tF9y/j09L9CX3/ujr3F/0r5TfS5430fBCdseGCe0fFxEEQxtDLN82XSg5HjejLEwtnsL+2novuX8bGzysS9t7hcDhh76WUFNLnjfR8EJ1R2mnHIIyhl2+aLz3tmDwmDuvNk9fMoD4c5qL7l/Gfz/Yl5H11o12l2k/6vJGeD6IzhoRdcB+EMfTyTUejF9wnl6zBh/DkNTNJMfD1Bz5gzfY9cX/PsrKyuL+HUtJInzfS80F0xrCwla8gjKGXb5qvtLQ01yUojzGH9uRv186kW1oqlz64nI937I3r+w0dOjSux1dKIunzRno+iM4YsrKaryCMoZdvmq/aWrcbfKqWjRrQgyeumUHXLqlc8uDyuJ6C3Lx5c9yOrZRU0ueN9HwQnTEclnXaMQhj6OWb5ksfL5S8RvZvaMDSUg0XL1rGJzvj04BlZWXF5bhKSSZ93kjPB9EZpa18BWEMvXzTfFVWJn5fKdV2owf04PGFM0hJMcxftJwNu2J/F2R+fn7Mj6mUdNLnjfR8EJ0xFEbUylcQxtDLN81XEO+G8JvDB/bkiYUzAJi/KPbbUEydOjWmx1MqCKTPG+n5IDpjwwX3joqJgyCMoZdvhk9XvvxhzKE9eWLhdKy1zH9gGZtLYrcRa25ubsyOpVRQSJ830vNBdEZppx2DMIZevmm+dOXLP44Y1IslV8+gPtzQgMVqJ/ycnJyYHEepIJE+b6Tng+iM0i64D8IYevmm+dKVL385cnAvllw9nZr6EPMfWEZRaefHLy8vLwaVKRUs0ueN9HwQnTEsbOUrCGPo5ZvmS1e+/GfckEN47OrpVNaFmL9oGVvLOteATZkyJUaVKRUc0ueN9HwQnVHaDvdBGEMv3zRf1dXVrktQHTBhaG8eWzCdfdV1zF+0jOLyjjdg69ati2FlSgWD9HkjPR9EZwwLu9sxCGPo1ebmyxhzmjFmkTFmSuPH18SvrGjp6emJfDsVQxOH9WbJ1TPYU9XQgG3fXdWh44wePTrGlSmVPIwxi40xu4wxH8XyuNLnjfR8EJ0xbC2CzjoGYgy92rPydQPwfeBSY8zJQELXCevq6hL5dirGJmU2rIDt3t/QgO3Y0/4GbPv27XGoTKmk8RfgzFgfVPq8kZ4PojNKu+YrCGPo1Z7m63Nr7W5r7feA04Gj41RTi7p06ZLIt1NxkD28D48sOIbSilouXrScnXvbdyq5X79+capMqQbGmBXGmIeMMd82xpxsjBmYqPe21r4NxPwJw9LnjfR8EJ0xFAYj6LRjEMbQqz3N10tNv7HW/hfwSOzLObBwOJzIt1NxctSIvvz1qqPZtbea+Q8sY1c7GjC941UlwBzg70A6cB2wxRjzqduSIhljrjHGrDTGrNyxYwclJSXs2LGDbdu2UV5ezsaNG6mqqmLt2rWEw2HWrFkDfLGXUl5eHuFwmLVr11JVVcXGjRspLy9n27ZtNB1vy5YtVFRUsG7dOurr6ykoKIg4RtN/CwsLqamp4ZNPPmHv3r0UFRWxa9cudu3aRVFREXv37uWTTz6hpqaGwsLCFo9RUFBAfX0969ato6Kigi1btrSaqenuuNzcXCorK8Vl8o7T7t27IzKFQiHqaqp9nenL41RcXCxinLx/9w46j621B3+BMXcBt9jWXhhn2dnZtmkQlP+t2FLGFYs/ZEjvDJ68ZiYDe3Vt9Wt27NjBkCFDElCd6ghjTK61dprrOmLJGDMOmGetvS1B7zcKeNFaO7Etr582bZpduXLlQV8jfd5IzwfRGaf84jW+mj2UX8xp01+TpCd1DA/2PbEtK18VwD+NMd0bD3a6Mea9WBbYFpKWWBUcPaofD195NNt3V3PxomWUVNS0+jVpaWkJqEwFmTFmxJc/ttZ+DExwVE5MSJ830vNBdMaQsE1WgzCGXq02X9baHwNPAG8ZY94Fvgv8V7wL89LTjvJMP6w/i688mq3llVyyaDmlrTRgFRWxf1i3Uh5PGWOKjTHvGGP+bIy5E8hyXVRnSJ830vNBdEZrZW01EYQx9Gq1+TLGnAIsBPYDA4GbrbXvxLswL73gXqaZh/fnoSuOZkvpfi55cDll+2sP+NoBAwYksDIVRNbamdbaTOAbwOvAGuCcRLy3MeYJ4APgyMYGcEEsjit93kjPB9EZQ8IerB2EMfRqy/D9D/ATa+1sYB4NPxmeHNeqWlBbe+B/lJW/HTdmAA9eMY1NJfu59MHl7K5seayLi4sTXJkKKmvtBmvtP6y1D1lrE/IXz1o731o7xFqbZq3NtNY+FIvjSp830vNBdMawsB3ugzCGXm057Xiytfbdxt8XAmcBv4x3YV5du7Z+Qbbyr+OPGMiiy6ex4fMKLnlwOXsqo/d1GzNmjIPKlPI36fNGej6Izhi2lhRB+3wFYQy92r1waa3dAZwSh1oOSh8vJN+JYwdy/2U5fLKzgksfWs6eqsgGrOmWeaVU20mfN9LzQXTGsIVUQStfQRhDrw6dNbbWduz5MJ3QrVu3RL+lcuCkIw/l3kunsu6zvVy++EP2Vn/RgGVnZzusTCl/kj5vpOeD6IwNdzs6KiYOgjCGXr65ZE832AyOU8YN4p6Lp7Jm2x6uWPwh+xobsKZN7ZRSbSd93kjPB5EZm7bclHTaMQhj6JWQ5ssYk2qM+YMxZo0xptAYc1h7j9G9e/d4lMhnxbUAACAASURBVKaS1OkTBvOni6dSWLyHKx9eQUVNPTk5Oa7LUsp3pM8b6fkgMmMo3Nh8CTrtGIQx9ErUytePgE3W2gnA3TQ8pLtddOUreM6cOJg/zj+K/K27ueWp/ED+dKRUZ0mfN9LzQWTGxt5L1IO1gzCGXnFvvowxPYDzrLV/aPzUZqDdtzboylcwnTVpCDeeNIbX1+4kbdDhrstRynekrypIzweRGcONpx0FLXwFYgy9ErHydSow3BiTb4zJBxYDZW35wi8/QLa4uDimD7yU+BBPiZk2FRXzYn4xfTJSKSn6REQmieO0Y8eOg01l5VDT3wmppOeDyIxNzZekux2DMIZerT5Yu9NvYMytwE5r7X2NHz8IrAaKgK8AhwL3WGtfO9hxcnJybBCXJoPMWst3/1bAP/K38ehV0zl6RC/d7y2JSXywdrJry4O1a2pqRM8b6fkgMuO+6jom/ew1/ufscSw8od2XTyclqWPY2Qdrd1ZfoLKxkC7A6cAL1trnrLULgSuBi1o7iO5wHzxPrdjKs6u2sXDGMGYdMYCioiLXJSnlO9LnjfR8EJmx6THHku52DMIYeiWi+VoPzGj8/S3AS9bazV/68x8D97R2EH22Y7Cs3b6Xnz5XyKiMar5z5ngABg0a5LgqpfxH+ryRng8iMzaddhTUewViDL0S0Xw9AUw1xmwAJgPfATAN/g94xVqb19pBQqFQfKtUSWNfdR03LMkljXoe+MaxZHRNB2D37t2OK1PKf6TPG+n5IDJjqOmaL0HdVxDG0Cvuy0nW2nK+WPn6sm/ScDF+b2PMmKZrwg4kJcU3+8GqTrDW8oOnC9haXsVjVx3L2JEDm/8sIyPDYWVK+ZP0eSM9H0Rm/OJuRznNVxDG0MvZuTxr7d007PmlVLNHPtjCKx/t5OIJ3Zk5ZmCrr1dKqSBpuuZL0t2OQeSb5aRw0984JVbB1t3c9uJaDsuo5KcXHhv15/pwdaXaT/q8kZ4PIjNKvOYrCGPo5ZvmKzU11XUJKo72Vtdxw5KVZFDLfd84joyM6NuO+/Tp46AypfxN+ryRng8iMzY/XkhQ9xWEMfTyTfNVX1/vugQVJ9ZafvRsIZ/treUPF2UzduTQFl+3c+fOBFemlP9JnzfS80FkxqatOSWddgzCGHr5pvlKT093XYKKk6dWbOWl1TtYcMwgTsk+8KaBI0aMSGBVSskgfd5IzweRGZvudpR0D1oQxtDLN8MXxHPCQbB+5z5uff4jRqRXcssZEw7+2vXrE1SVUnJInzfS80Fkxi+u+ZKz8hWEMfTyTfPVrVs31yWoGKuuC3H9oytICddxz2XH0K3bwW83njRpUoIqU0oO6fNGej6IzBgOy2u+gjCGXr5pviorK12XoGLsFy+uZWNJFT86eRiTjhjZ6uv12Z5KtZ/0eSM9H0RmlLjJahDG0Ms3z+zp3r276xJUDL20egePLy9iwbEjuPy0tv3Uk5OTE+eqlJJH+ryRng8iMzY/21FO7xWIMfTSlS+VcFvLKvnB0/kM6lLFN2ePavPXBfGnI6U6S/q8kZ4PIjNKvOYrCGPo5ZvmS1e+ZKgLhblpyUrq6ur4/YWT6HNIrzZ/bRB/OlKqs6TPG+n5wLPyJbD5CsIYevmm+aqqqnJdgoqBO19fT8G2fSzM7sGx2Ue262sLCgriVJVSckmfN9LzQWTGpk1WJV3zFYQx9PLNNV9BfPCmNG+v/5x739zI/GNG8P3z2393y4QJB9+KQikVTfq8kZ4PIjM29l4IWvgKxBh6+Wblq6amxnUJqhN27avmW0/k0S+1mu+fOrpDx9iwYUOMq1JKPunzRno+iMwYFni3YxDG0Ms3K1+6w71/hcOWW55cxb6qWu46ZzT9DunZoeNkZmbGuDKl5JM+b6Tng8iMEvf5CsIYevlm5Uuf7ehf97+9ifc2lnHBYZZzZk3t8HFKSkpiWFX8lFbUULB1N6UVulqr3PPLvOko6fkgMmNI4AX3QRhDL9+sfKVIepBVgKwqKueO1/7D2RMHc9vXOreLcc+eHVsxS6Tn87fxg6cL6JKSQshafnvBZL46ZZjrslSA+WHedIb0fBCZsenB2oLOOgZiDL1809HYpr9xyjcqaur55uN5dE+p45dzxnf61HFdXV2MKouP0ooafvB0ATX1lv21IarrwvzgmdW6AqacSvZ501nS80FkRol3OwZhDL1803wp/7n1+Y/YtruKb0/vS79enX82Z7hpa+ck9WlpJSFPiWkpKRSX6zYpyp1knzedJT0fRGZs3udLUPMVhDH08k3zpacd/eWFgu08k7eNEwdW841zjo/JMZN5o93a+jB/WvoJ9eHIFdq6cJjMvvpQeOVOMs+bWJCeDyIzStxkNQhj6OWbjkYvuPeP4vJK/vsfhUwY3J3fX306JkbfJMrKymJynFirrgtx7aMreWPd55x/1DAy0lLo1bULGWkp/PaCyfTv2bXT77F//36uuOIKFi5cyJIlSxL+9cq/knXexIr0fBCZsWl1PVVQ8xWEMfTyTfOVlpbmugTVBqGw5dtPrKKurp4/X3YMfXsfErNjDx06NGbHipXK2noW/HUFb67/nF+fN4k7L5rCez88mceuns57Pzy5TRfbV1VVceKJJxIKhQ74mmeffZZ58+axaNEi/vnPf7a7zpa+vra2lhNOOEF/sBEuGedNLEnPB5EZm1a+BPVegRhDL980X7W1ta5LUG1wz9INrCzazdcOh5H9e8T02Js3b47p8TprX3UdVyz+kA82lnLHhdlcPH0EAP17diV7eJ82r3gtXryY888/n9TU1AO+pri4mOHDhwMc9HXt+fr09HROOeUUnnrqqXYfT/lHss2bWJOeDyIzhgVecB+EMfTyTfOljxdKfnlF5dz17/WM617Bjy89LebHz8rKivkxO2p3ZS2XPricVUW7+eP8qZw/teObBC5ZsoQ5c+YAcNJJJ/H6668D8OMf/5ibb74ZaNiEsLi4GGj94tSWjnGgr587d66ehhQumeZNPEjPB5EZw81bTchpvoIwhl6+2eersrLSdQnqIPZV1/HtJ/Ppl2H44xWz4vJEgvz8fKZO7fgmrbFSUlHDpQ8uZ9Pn+7nv0hxOHT+ow8eqra1l06ZNjBo1CoCf//zn/PSnP2XXrl2sWrWq+RTh+eefz0033cRLL73Eueeee9BjtnSM6urqFr9+4sSJrFixosP1q+SXLPMmXqTng8iMoebHC7msKLaCMIZevmm+gng3hJ/c+vwaissr+du1Mxkzsl9c3iMZJufOvdVcvGgZ23ZX8dCV0zj+iIGdOl5JSQl9+vRp/viEE07AWsudd97Jm2++2XyKsEePHjz88MNtOmZLxzjQ16emppKens6+ffvo1atXp7Ko5JQM8yaepOeDyIy2+ZovOStfQRhDL9/0zrrylbz+WbCdZ1c1bCsxbVR8Gi+A3NzcuB27LYrLK/na/R/w2Z5q/vqNYzrdeAF069aN6urq5o8LCwvZsWMHXbt27XAz1N5j1NTU6Gl9wVzPm3iTng8iMzZvsiqo+QrCGHr5pvnSla/kVFxeyX8/U8DgLpXcuSD213l9WU5OTlyPfzCbS/Zz0f3LKN9fy2NXT2f6Yf1jcty+ffsSCoWorq5mx44dXHLJJTz//PP06NGDV199FYBNmzaxYMEC5s2bF/X1p5xyCtu2bWv++EDHeO6551i4cCFz5szhtddea359aWkpAwcO1LuJBXM5bxJBej6IzCjxmq8gjKGXb5ovXflKPvWhMLc8lU9dfYjffPXImG4r0ZK8vLy4Hv9A1mzfw4X3vU9VXYjHF87gqBF9Y3r8008/nddee43zzz+fO+64g3HjxvGTn/yEn/3sZwAcdthhPPTQQ1FfFw6H2bBhA/36Naw2VlZWHvAYc+fOZdGiRfzlL3+JuLtx6dKlnH322THNo5KLq3mTKNLzQWTGprsdJe07HoQx9NJrvlSH/fnNjazYUs7t8yZz0rQRcX+/KVOmxP09vFZsKeOqh1fQK6MLjyyYzphDY/8A2Jtuuok777yTDz74oPlzJ5xwQsTHLVm7di0XXHAB3bo17KDfvXv3Vo/xy1/+khtvvLH548cff5zf/OY3sYihkpSLeZNI0vNBZEaJO9wHYQy9fNM7f/m6GOVe07YSU/rWMy8BjRfAunXrEvI+TZau28VlDy1nYK+u/P36Y+PSeAEcddRRnHTSSQfdZLUlEydO5M4772zTa621/PCHP+Sss85qvri1traWuXPncuSRR7a7ZuUfiZ43iSY9H0Rm/OJuRznNVxDG0Ms3zVc8ti5QHbOvuo6bH8+jh6nlf+cl7ieW0aNHJ+y9ns/fxsJHVjLm0J787bqZDOsT3+czXnXVVQfcPLW0tJTrrruOVatWdXiV6o9//CP//ve/efrpp7nvvvuAhjl1+eWXd7hm5Q+JnDcuSM8HkRmbTjsKWvgKxBh6+ea0Y11dnesSVKPbXlzLtt1V/Pf0Q8g6fGTC3nf79u0cfvjhcX+fR5d9yk+f/4ijR/XjoSum0SvD7cXo/fv3b26YOurmm29u3rBVBUui5o0r0vNBZMamC+4l3e0YhDH08k3z1aWLb0oVrWDrbv62spj5UwawYM4xCX3vpgvL48Vayz1LN3D7a+s5JetQ7rlkKhlp7X+Uj1LJJN7zxjXp+SAyY0jg44WCMIZevjnt2NojVVRiLP2oCIDPqlIoqUjs8zbjecertZZfvfQxt7+2nvOOGsZ9l+Vo46VEkH6nuPR8EJkxLHCT1SCMoZdvmi/lXigUouunH3Bldi/e31jK6Xe9zQsF2xP2/ilxure6PhTmB0+v5sF3N3PlsaO448Js0iQ9u0MFWrzmTbKQng8iM4YFXnAfhDH08k1iSV2+Xy1dupS+fXpz69eP56Wbj2dk/x5884lV3Ph4HmX7478KFo+NQKvrQtywJI+/5xbzrVOO4NZzx5Mi6JuaUtI30JWeDyIzfrHJqqNi4iAIY+jlm+ZLTzu6Za2ltraWc889F2MMYw7tyTPXzeT7ZxzJa2s+4/Tfv82/1+6Maw0VFRUxPd7uyloufXA5r3+8k1vPHc8tp43VJl+JE+t5k2yk54PIjE3XfEna5ysIY+jlm+ZLL7h3p6amhrKyMs4++2x69OjR/PkuqSnceNIYnr9xFgN6pnP1Iyv5/t8LqKipj0sdAwYMiNmxtu+u4sL7PmB18R7+OP8ovnFc8G51VsEQy3mTjKTng8iMVuAmq0EYQy/fNF+1tYm9uFt94V//+hfLly8/4J+PH3oI/7xpFjeedDjP5BVzzt3vsLp4d8zrKC4ujslx/vPZPs7/8/sND8i+6hjOmTw0JsdVKhnFat4kK+n5IDJjqPEkkKRrvoIwhl6+ab66du3quoRA+vjjj/n000859dRTD/q69C4pfP+MLJ68Zia19WEuuPd9Hnh7Y/OGgLEwZsyYTh9j+aZS5t33PmFr+dt1M5l5eGwekK1UsorFvElm0vNBZMYvHi/kqprYC8IYevmm+dLHCyVeOBxm6dKlnHfeeW1+wsAxo/vx8reO55SsQfz65XVc8fCH7NoXm7Fbs2ZNp77+lcIdXLb4Qw7t1ZVnbziWcUPi+yBwpZJBZ+dNspOeDyIzhq3FGFk3oQVhDL1803w1PTxYJYa1FmMMCxcuZPjw4e362j7d07n30qn86ryJfLi5jLP/8A5L/7Or0zVlZ2d3+Gsf+WALNzyex8Shh/D0dceS2Vcf1K6CoTPzxg+k54PIjKGwFXW9FwRjDL1803wFcRM2l/Ly8nj11Vc7fAuwMYZLpo/khW/OYkDPrnzj4RXc9uJa6kIdv2s1Nze33V9jreV3r67jp8+v4ZSsQSy5egZ9e+hzQlVwdGTe+In0fBCZMWxlPVoIgjGGXr5pvrp315WKRCkvL+eNN95g6tSpnT7W2EG9eO7G47h85kgeencz8x9Yxq69HTsNmZOT067X14XCfP/p1dyzdCPzjxnBfZdOpVu67lqvgqW988ZvpOeDyIxNpx0lCcIYevmm+dKVr8QIh8M899xzzJo1i0MPPTQmx8xIS+UXcyZy9/yjWLN9L2ff/S7LN5W2+zjt+eloT1UdVz78IU/nFvPtU4/g1+dNpIvuWq8CSPqqgvR84Fn5CltRdzpCMMbQyzf/GunKV2IYY5g5cyYzZsyI+bG/mj2U5286jl4ZXbj4weU8+M6m5j1r2qKtPx1tLatk3r3v8+HmMm6/MJtvn6qbp6rgkr6qID0fRGYMWSvutGMQxtDLN81XVVWV6xLE27lzJ2vWrCErKytuzcrYQb14/qbjOCXrUH750sd884lV7G/jpqyFhYWtvqZg627O+/P77NzbsIfXvJzMzpaslK+1Zd74mfR8EJnRWsSddgzCGHr5pvnKyMhwXYJooVCIf/zjH9TV1cX9vQ7JSOP+y3L44ZlZvFy4g7n3vMfGz1t/vMTYsWMP+uf/+ugzLnrgA7qlp/DsDcdy7OHB2zVZKa/W5o3fSc8HkRlDAk87BmEMvXzTfOkO9/H15ptv0qdPH6ZMmZKQ9zPGcP3sw3l0wXRK99cy50/v8a+Pdhz0a4qKilr8vLWWRW9v4voluWQNPoR/3HAcYw7tFY+ylfKdA80bKaTng8iMYStvq4kgjKGXb5ovfbZj/ITDYcrKyjjnnHMSfm3UcWMG8OI3Z3H4oT257rE8/veVddQfYDuKQYMGRX2uPhTmJ89/xK9e/pizJg7myWtmMKCnPg1BqSYtzRtJpOeDyIxha0kRtvIVhDH08k3zFQqFXJcgUl1dHVVVVVx44YX07NnTSQ1D+3Tjb9fO4JLpI7jvrY1c9tCHlFTURL1u9+7I50VW1NRz9SMreWxZEdeeeBh/mj+VjDTdSkKpL/POG2mk54PIjOGwrEcLQTDG0Ms3zVdKim9K9ZXXX3+dN99803UZdO2Syq/Om8Tv5k0mr6icc//4LquKyiNe8+Xr/pruaHznkxJ+fd4kfnTWOHE/DSoVC9Kvl5WeDyIzSrzbMQhj6KUdTYBt2rSJ//znP5x88smuS2l24bThPHP9sXRJNXzt/g94bNmnUdtRLNtUypx73mPb7ioWX3k0F08f4ahapZRKrHDjo9+Uv/mm+QqHO/5YGhUtHA7z8ssvc+655ybdczMnDuvNCzfN4rgxA/jxcx9x85P57KlsOD36l/c2c+mDy+nTPY3nbzyOE8cOdF2uUkmtujo2D7ZPVtLzQWRGiZusBmEMvXxzFXtqql7LE0spKSlceeWVzq7zak2f7uksvuJo/vzmBu769ye8vf5zenVNpXh3NSdnHcpdX5/CIRkde+6kUkHSp08f1yXElfR8EJkxZOVd8xWEMfTyzcpXfX3bNuJUrVu/fj1vvPFG0jZeTVJSDDedfATP3XgcJ2cdSmavVH5/UTYPXTFNGy+l2mjnzp2uS4gr6fkgMqPEux2DMIZevln5Sk9Pd12CCJWVlbz44oucf/75rktps4nDevP7i6ZQU1ND1666jYRS7TFihOxrIqXng8iM4bC8fb6CMIZevln5CuI54Xh45ZVXGD9+PKNGjXJdSrutX7/edQlK+Y70eSM9H0RmDAu82zEIY+jlm5WvZLso3I+stQwbNsy3DzGdNGmS6xKU8h3p80Z6PojMGArLe7ZjEMbQyzcrX5WVla5L8LWKigq2bNnCjBkzSEvz5/VSubm5rktQynekzxvp+SAyo7Xy7nYMwhh6+ab56t69u+sSfMtay0svvcTGjRtdl9Ipfl2xU8ol6fNGej6IzBgS2HwFYQy9fNN86cpXxxUWFlJaWsrs2bNdl9IpQfzpSKnOkj5vpOeDyIxhi7hNVoMwhl6+ab505atjrLWsWrWKuXPn+v7h5EH86UipzpI+b6Tng8iM4bAlVVbvFYgx9PJN81VVVeW6BN+x1hIKhbj88ssZOnSo63I6raCgwHUJSvmO9HkjPR9EZgxbeVtNBGEMvXzTfAXxwZudVVBQwLPPPitmiXrChAmuS1DKd6TPG+n5IDJjKCxvk9UgjKGXb5qvmpoa1yX4yt69e3n99dc54YQTXJcSMxs2bHBdglK+I33eSM8HkRkbVr4cFhMHQRhDL980X7rDffu89NJLHH300QwePNh1KTGTmZnpugSlfEf6vJGeDyIzhi3i7nYMwhh6+ab50mc7ts/xxx/P8ccf77qMmCopKXFdglK+I33eSM8HkRlDAh8vFIQx9PJN85WS4ptSndq7dy/vvPMOmZmZpKamui4nppL9QeBKJSPp80Z6PojMaAVecB+EMfTyTUdjrXVdQtKz1vLCCy8QDoddlxIXdXV1rktQynekzxvp+SAyY0jgNV9BGEMv3zRfqnX5+flUVFQwa9Ys16XEhdSmUql4kj5vpOeDyIzhsLxrvoIwhl6+ab70tGPrSkpKmDNnjrjTjU10o12l2k/6vJGeDyIzhq0Vs31QkyCMoZdvOhq94P7ArLWUl5dz2mmnibq70ausrMx1CUr5jvR5Iz0fRGYMW0uqsOYrCGPo5ZvmKy0tzXUJSaugoIBnnnlG/HVxEnbpVyrRpM8b6fkgMmPDJqsOi4mDIIyhl2+GsLa21nUJSalpM9VzzjlH3FK01+bNm12XoJTvSJ830vNBZEZrEXe3YxDG0Ms3zZc+Xqhlb7zxBtOmTRN9urFJVlaW6xKU8h3p80Z6PojMGLJW3AX3QRhDL980X5WVla5LSEpnnHGGqEcIHUx+fr7rEpTyHenzRno+iMwocZPVIIyhl2+aryDeDXEw+/bt46mnnqJr165i7270mjp1qusSlPId6fNGej6IzCjxtGMQxtCri+sC2kpXvr5greXll19mwIABgdqCIzc3l5ycHNdlKHVQxphlwNettVuMMcOA562101zV88CL79OlzxBXbx93W4u3MjxzuOsy4urLGfdW14nbZDWI39t903zpytcX1qxZQ2lpKRdccIHrUhIqaJNT+Y9puOtlBPBp46cmA4Uxfo9rgGug4S6xkpIS6urqCIfDdO/enbKyMoYOHcrmzZvJysrizS1VvF+8NpYlJJ984fkgImPvtBAlJSVUVFQwYMAAiouLGTNmDGvWrCE7O7u5mWn6b2FhIWPHjqWoqIhBgwaxe/fu5uuoq6ur6dOnDzt37mTEiBGsX7+eSZMmRR2joKCACRMmsGHDBjIzMykpKaFnz54H/buXn5/P1KlTm4+Rl5fHlClTWLduHaNHj2b79u3069ePwYMHs2PHDtLS0sRkam3ByPhle4Lx48fbtWsDMMHa4OOPP+aQQw5h2LBhrktJqLy8vEAuT/uFMSbX5QpPMjDGHAH8xVp7XOPHPwTqrLV3tvJ1t1lrf9Le95s2bZpduXLlQV/z3vKVTJyU3d5D+0bB6gKyJ8vNB56MBg7J6CLq7nap39sP9j1RV758ZuPGjWRlZYmaeG01ZcoU1yUo1ZpJRK50TQPuN8YMBp4CXgImAO8DpwE/A0qALsaYTOBR4J/ADGvtRbEoaObRU0VfnjDrGNn5QH7GIH5v981oVldXuy7BuY8//piXX345sLv9r1u3znUJSrWmH1AFYIwZB3yFhmbsKOBZa+1vgd7AIuDvwMjGP8sHsoHnrLW/B2I2yaXPG+n5QH5G6fla4pvmKz093XUJTlVVVfHKK68wZ86cwO72P3r0aNclKNWaV4FTjDF/Ay4ESq21O4EpwKvGmLTGz4WBiTQ0ZlP4ovl6tfE4MbseRPq8kZ4P5GeUnq8lvmm+6urqXJfgVH5+PuPGjWPEiBGuS3Fm+/btrktQ6qCstVuttZOttV+z1v7CWtt0G94YYD0Npxw/bvzcKGttUeOffdL0GmPMAOCzWNUkfd5IzwfyM0rP1xLfXPPVpYtvSo05ay0zZswgHA67LsWpfv36uS5BqQ6x1i5o/G1+4y+stZd5/uyqxv+WAN+L1XtLnzfS84H8jNLztcTZypcx5kxjzH+MMRuMMf/V2uuD2njU1NRw//33U1lZGZjNVA9E93pTqv2kzxvp+UB+Run5WuKk+TLGpAL3AGcB44H5xpjxLmpJdq+//jpDhw6lR48erktxTvLdPkrFi/R5Iz0fyM8oPV9LXCU+Bthgrd1kra0FngTmHOwLgri1wpYtW1i/fj2nn36661KSQlBvNFCqM6TPG+n5QH5G6fla4mSTVWPMPOBMa+3VjR9fBky31t7keV3zTs403Bn0UUILVclmAA3Xw6jkNNJaO9B1EUFijPmcL3bTPxDp80Z6PpCfUWq+A35PdHUVe0vLWFFdoLX2AeABAGPMyqDvnh10+ndAqUhtaXalzxvp+UB+Run5WuLqtGMx8OUnoWYCwbvXVCmllFKB46r5WgEcYYwZbYxJB75OwyM1lFJKKaVEc3La0Vpbb4y5iYbdnFOBxdbaNa182QPxr0wlOf07oFT7SZ830vOB/IzS80VxcsG9UkoppVRQBW9zDaWUUkoph7T5UkoppZRKoKRvvtr7GCKllAqi1r5XmgZ3N/75amPMVBd1dlQb8s02xuwxxuQ3/vqpizo7yhiz2BizyxjT4n6WAsavtXy+Hr/2SurmSx9DFGzGmFRjzB+MMWuMMYXGmMNc16RUMmrj98qzgCMaf10D3JvQIjuhHf8WvGOtndL46xcJLbLz/gKceZA/9+34NfoLB88H/h6/dknq5osOPIZIifIjYJO1dgJwN3CD43qUSlZt+V45B3jENlgG9DHGDEl0oR0k/t8Ca+3bQNlBXuLn8WtLvkBJ9uZrGLD1Sx8XN35OCWeM6QGcZ639Q+OnNgNjHJakVDJry/dKP38/bWvtM40xBcaYV4wxExJTWsL4efzaSvL4RXD1eKG2atNjiJRIpwLDjTH5jR/3A/7tsB6lkllbvlf6+ftpW2rPo+FZehXGmLOB52g4RSeFn8evLaSPX4RkX/nSxxAF1xTgp03n/4HXgHxjzGHGmIeMMU87rk+pZNKW75V+jGcSbgAAAy5JREFU/n7aau3W2r3W2orG378MpBljBiSuxLjz8/i1KgDjFyHZmy99DFFw9QUqAYwxXYDTgRcar/lY4LQypZJPW75X/hO4vPGuuRnAHmvtjkQX2kGt5jPGDDbGmMbfH0PDv2+lCa80fvw8fq0KwPhFSOrTjh18DJGSYT0wA3gEuAV4yVq72W1JSiWnA32vNMZc1/jn9wEvA2cDG2j4weYbruptrzbmmwdcb4ypB6qAr1sfPcLFGPMEMBsYYIwpBm4F0sD/4wdtyufr8WsvfbyQSkrGmL7AK8AA4APgGmtt1Zf+/Glr7TxX9SmllFIdpc2X8hVjTH/gV8BpwIPW2t84LkkppZRqF22+lFJKKaUSKNkvuFdKKaWUEkWbL6WUUkqpBNLmSymllFIqgZJ6qwmllFLKjxpvDvp/jR8OBkLA540fV1prj3VSmEoKesG9UkopFUfGmJ8BFdba213XopKDnnZUSimlEsgYU9H439nGmLeMMX8zxqw3xvyvMeYSY8yHxphCY8zhja8baIx5xhizovHXcW4TqM7S5ksppZRyJxv4FjAJuAwYa609BngQ+Gbja/4A/N5aezRwQeOfKR/Ta76UUkopd1Y0PaPRGLMReK3x84XASY2/PxUY3/joQ4BDjDG9rLX7Elqpihld+VIJYYxZaow5rfH3vzTG3O26JqWUSgI1X/p9+Esfh/ligSQFmGmtndL4a5g2Xv6mzZdKlFuB/zHGXAIcRcPDspVSSrXuNeCmpg+MMVMc1qJiQJsvlRDW2rcBA3yHhqfVh4wxhxljHjLGPO24PKWUSmY3A9OMMauNMWuB61wXpDpHt5pQCWGMmQQ8A5R497cxxjxtrZ3npjKllFIqsXTlS8WdMWYIsASYA+w3xpzhuCSllFLKGW2+VFwZY7oDzwLftdZ+DNwG/MxpUUoppZRDetpROdP4+I1fAacBD1prf+O4JKWUUirutPlSSimllEogPe2olFJKKZVA2nwppZRSSiWQNl9KKaWUUgmkzZdSSimlVAJp86WUUkoplUDafCmllFJKJZA2X0oppZRSCaTNl1JKKaVUAmnzpZRSSimVQP8fZBQQe27nrewAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Final time is 1.88756560904673\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "solution = bocop.readSolution(problem_path + \"/problem.sol\")\n",
    "tf = solution.state[2][0]; t = solution.time_steps*tf\n",
    "x1 = solution.state[0]; x2 = solution.state[1]\n",
    "u = solution.control[0]\n",
    "\n",
    "umin = 0.5\n",
    "umax = 1.5\n",
    "g1 = 1.4\n",
    "g2 = 1.8\n",
    "theta1 = 0.6\n",
    "theta2 = 0.4\n",
    "k1 = 1\n",
    "k2 = 1\n",
    "x1max = 0.9\n",
    "x2max = 0.8\n",
    "\n",
    "plt.figure(0, figsize=(10,4))\n",
    "plt.subplots_adjust(hspace=0.4, wspace=0.4)\n",
    "plt.subplot(121)\n",
    "ax = plt.gca()\n",
    "plt.xlabel('$x_1$'); plt.ylabel('$x_2$')\n",
    "ax.set_ylim(0.,x1max); ax.set_xlim(0.,x2max)\n",
    "xv = np.arange(0,theta1,0.01)\n",
    "sumax = k2/g2*umax - (k2/g2*umax - theta2)*((k1/g1*umax - xv)/((k1/g1*umax - theta1)))**(g2/g1)\n",
    "plt.plot(x1,x2,label='',zorder=3)\n",
    "plt.plot(xv,sumax,'--', label='$(S_{u_{max}})$',linewidth=1, color='grey')\n",
    "plt.legend()\n",
    "plt.xlim([0,theta1*1.8])\n",
    "ax.text(x1[0]+0.04, x2[0]-0.05, '$(x_1^0, x_2^0)$')\n",
    "plt.scatter([x1[0],x1[-1]],[x2[0],x2[-1]], zorder=3, s=20)\n",
    "plt.xticks([0, theta1], [0,'$\\\\theta_1$'])\n",
    "plt.yticks([0, theta2, x2[-1]], [0,'$\\\\theta_2$','$x_2^f$'])\n",
    "plt.grid(linestyle='dotted')\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.plot(t[1:],u)\n",
    "plt.xlabel('Time'); plt.ylabel('$u$')\n",
    "plt.yticks([umin,1,umax], ['$u_{min}$', 1, '$u_{max}$'])\n",
    "plt.grid(linestyle='dotted')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "print(\"Final time is \" + str(tf))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bocop returns status 0 with objective 1.888 and constraint violation 1.486e-14\n",
      "Costate at first time stage (t0+h/2):  [-2.21241861069235, 0.680215961376846, 0.99503700918767]\n",
      "Multipliers for initial conditions:  [-2.18337787e+00  6.68757718e-01  1.38521948e-12]\n"
     ]
    }
   ],
   "source": [
    "print(\"Bocop returns status {} with objective {:2.4g} and constraint violation {:2.4g}\".format(solution.status,solution.objective,solution.constraints))\n",
    "p0 = []\n",
    "for i in range(solution.dim_state):\n",
    "    p0.append(solution.costate[i][0])\n",
    "print(\"Costate at first time stage (t0+h/2): \",p0)\n",
    "print(\"Multipliers for initial conditions: \",solution.boundarycond_multipliers[0:solution.dim_state])"
   ]
  }
 ],
 "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}