Sparse-Identification-with-SINDy-on-pendulum-with-control.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test SINDy on a Simple Pendulum System\n",
    "\n",
    "This code demonstrates how to use the SparseNonLinearModel class to identify the governing equations of a non-linear dynamical system with a forced input.\n",
    "\n",
    "The system used in these examples is a [simple pendulum system](https://services.math.duke.edu/education/ccp/materials/diffeq/pendulum/pend1.html) of ordinary differential equations.\n",
    "\n",
    "Contents:\n",
    "1. Generate data by simulating pendulum\n",
    "2. Identify system dynamics using sparse regression\n",
    " 1. Using SparseNonLinearModel class\n",
    " 2. Using identified model for prediction\n",
    " 3. Estimation from 10 data points\n",
    "\n",
    "References:\n",
    "\n",
    "1. Brunton, S. L., & Kutz, J. N. (2019). Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control (1st ed.). Cambridge University Press. http://www.databookuw.com"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import partial\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import odeint\n",
    "from scipy.interpolate import interp1d\n",
    "\n",
    "# Local modules\n",
    "from dynopt.models.models import NonLinearModel, SparseNonLinearModel\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Generate data by simulating pendulum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# System equations\n",
    "\n",
    "def pendulum_dydt(t, y, u, m=1, L=2, g=-10, b=1):\n",
    "    \"\"\"Simulates the non-linear dynamics of a simple pendulum system.\n",
    "    These non-linear ordinary differential equations (ODEs) return the\n",
    "    time-derivative at the current time given the current state of the\n",
    "    system.\n",
    "    \n",
    "    Args:\n",
    "        t (float): Time variable - not used here but included for\n",
    "            compatibility with solvers like scipy.integrate.solve_ivp.\n",
    "        y (np.array): State vector. This should be an array of\n",
    "            shape (2, ) containing the current state of the system.\n",
    "            y[0] is the angle (theta) in radians of the pendulum from\n",
    "            the vertical down position and y[1] is the angular \n",
    "            velocity (dtheta/dt).\n",
    "        u (function): u(t) is a function which should return the \n",
    "            torque applied to shaft.\n",
    "        m (float): Mass of pendulum.\n",
    "        L (float): Length of pendulum.\n",
    "        g (float): Acceleration due to gravity.\n",
    "        b (float): Damping coefficient for joint friction (proportional\n",
    "            to angular velocity).\n",
    "    Returns:\n",
    "        dydt (np.array): The time derivate of the state (dy/dt) as a\n",
    "            shape (2, ) array.\n",
    "    \"\"\"\n",
    "\n",
    "    dydt = np.empty_like(y)\n",
    "    dydt[0] = y[1]\n",
    "    dydt[1] = -(b/m) * y[1] + (g/L) * np.sin(y[0]) + u(t) / L\n",
    "\n",
    "    return dydt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = 0.01\n",
    "T = 10\n",
    "t = np.arange(0, T + dt, dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cart-pendulum system parameters (constants)\n",
    "params = {\n",
    "    'm': 1,  # Mass of pendulum\n",
    "    'L': 2,  # Length of pendulum\n",
    "    'g': -10,  # Acceleration due to gravity\n",
    "    'b': 1  # Damping coefficient for joint friction\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initial condition\n",
    "x0 = [np.pi, 0]  # Pendulum up position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array(-0.), array(-1.30057568))"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Control input\n",
    "u = -2*np.sin(t*(1+t/20))  # Sine function with varying frequency\n",
    "\n",
    "uf = interp1d(t, u, kind='previous')\n",
    "u = u[:-1]\n",
    "assert u.shape == (1000,), u.shape\n",
    "\n",
    "# u is a function defined in this range\n",
    "uf(t[0]), uf(t[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.000000e+00, -6.123234e-16])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dYdt = pendulum_dydt(t[0], x0, u=uf, **params)\n",
    "dYdt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAp4AAAEWCAYAAADPfZe7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABV5UlEQVR4nO3dd3hUZdrH8e/09IQkpBJCqKGT0DsqIKIoggKCHTu4Ivva1lWBVVld6y6KYsNVUVRQUFEBRXonIXQIEBIIpJKezExmzvvHJFmRIoHMnJnJ/bkur8sMQ84vORnmzlPuR6MoioIQQgghhBBOplU7gBBCCCGEaByk8BRCCCGEEC4hhacQQgghhHAJKTyFEEIIIYRLSOEphBBCCCFcQgpPIYQQQgjhElJ4CiGEEEIIl5DCUwghhBBCuIQUnkIIIYQQwiWk8BRCCCGEEC4hhacQQgghhHAJKTyFEF5t/vz5aDQatm3bpnYUADZs2MCMGTMoKipy2jVmzZpFhw4dsNvtdY8pisKsWbNYvXr1Gc/94IMPiI2Npby83Gl5hBCilhSeQgjhQhs2bGDmzJlOKzyzs7N5+eWXmTVrFlrt//6JP3jwIM899xwnT5484/l33HEH/v7+vPzyy07JI4QQvyeFpxBCeJE333yTkJAQxowZc8bj27dvB6B79+5nPK7X67n//vt58803qaiocFlOIUTjJIWnEKJRmTFjBhqNhj179nDLLbcQHBxMZGQkd999N8XFxed8bkpKCmPGjCEoKIjg4GBuvfVW8vLyznjunXfeSYsWLc57vdr/f+yxxwBISEhAo9Gg0Wj47bffLpg5IiKChx566KzHe/XqxciRI+s+tlgsfPDBB0ycOPGM0c7u3bszadIkANq2bYtGoyEwMBBFUQCYNGkSJSUlfPHFFxfMIYQQl0sKTyFEozR27Fjatm3LokWLePLJJ1mwYAGPPvroOZ9744030rp1a77++mtmzJjBt99+y9VXX43Vaq3XNe+55x4efvhhABYvXszGjRvZuHEjycnJ5/07x44dIy8v76znVFdXk5aWdsbjmzdvpqCggCuuuOKM586bN48ePXqQnJxcd81169bVFcRRUVEkJibyww8/1OvrEUKI+tKrHUAIIdQwefLkutHHoUOHkp6ezocffsgHH3xQV5DVGjNmTN0ayOHDhxMZGcmkSZP48ssv60YSL0azZs1o3rw5AElJSeccIf2j2k1Rf5wi3717N2az+YzCc+PGjQBnFandu3cnMzOTCRMm0KdPn3NeJzk5mZUrV1701yKEEJdCRjyFEI3S9ddff8bHXbp0oaqqitzc3LOe+8ficty4cej1elatWuXUjOAoPI1GI506dTrj8do1m78vMrOzs9FoNISHh5/x3KysLHJzc88qXn8vIiKC3NxcqqurGzC9EEKcSQpPIUSjFBYWdsbHJpMJgMrKyrOeGxUVdcbHer2esLAwCgoKnBewxtatW+ncuTMGg+GMx7dv305oaOgZo6aVlZUYDAZ0Ot1Zz4WzR01/z8fHB0VRqKqqarjwQgjxB1J4CiHEnzh16tQZH1dXV1NQUHBG8erj44PZbD7r7+bn51/WtVNSUs65BnTFihVnPR4eHo7FYjmrJ+f27dvx8/MjMTHxvNcpLCzEZDIREBBwWXmFEOJCpPAUQog/8dlnn53x8Zdffkl1dTVDhgype6xFixbk5uaSk5NT95jFYuHnn38+4+9eaGT1jywWC4WFhcTGxp7x+KJFi0hPTz+r8KwtLA8fPnzG42lpaSQmJp41Evp7R44coUOHDn+aSQghLocUnkII8ScWL17M448/zooVK3jjjTe4//776dq1K+PGjat7zvjx49HpdEyYMIFly5axePFihg8fjs1mO+Nzde7cGXD029y4cSPbtm2jtLT0nNc1Go20bNmSJUuWcPjwYbKzs3n77bf561//Cpy9iai2EN60adMZj4eEhHD48GGWLFnCpk2bOHbs2Bl/brfb2bJly1m74YUQoqFJ4SmEEH9i8eLF7N+/nzFjxvDss88yatQoli9fjtForHtOQkICS5YsoaioiJtuuonHHnuMm2++mdtvv/2MzzVkyBCeeuopvvvuOwYMGEDPnj3r1mCey8cff4zNZqNjx47079+fjIwMnnrqKeDswjMuLo6BAweyZMmSMx5/7rnn6NatG5MmTaJv374sWrTojD//7bffKC4urtcOfSGEuBQapbaDsBBCiDPMmDGDmTNnkpeXd9ZOcXe1aNEixo8fz7Fjx86aoj+f2267jSNHjrB+/XonpxNCNHYy4imEEF5kzJgx9OzZk9mzZ1/U8w8fPszChQt56aWXnJxMCCGk8BRCCK+i0Wh47733iImJwW63/+nzMzMzmTNnDgMGDHBBOiFEYydT7UIIIYQQwiVkxFMIIYQQQriEFJ5CCCGEEMIlpPAUQgghhBAuoVc7wIXY7Xays7MJDAxEo9GoHUcIIYQQQvyBoiiUlpYSExODVnvhMU23Ljyzs7OJi4tTO4YQQgghhPgTWVlZNGvW7ILPcevCMzAwEHB8IUFBQU6/ntVqZfny5QwfPhyDweD064mGJ/fQ88k99HxyDz2b3D/P5+p7WFJSQlxcXF3ddiFuXXjWTq8HBQW5rPD08/MjKChIXmweSu6h55N76PnkHno2uX+eT617eDHLImVzkRBCCCGEcAkpPIUQQgghhEs4tfCcPXs2PXv2JDAwkIiICEaPHs2BAweceUkhhBBCCOGmnFp4rl69milTprBp0yZWrFhBdXU1w4cPp7y83JmXFUIIIYQQbsipm4t++umnMz7+6KOPiIiIYPv27QwaNMiZlxZCeABFUbDY7FTbFAw6LUa9rP4RwhUURaHcYsOg02DUaaVXtnAZl+5qLy4uBiA0NPScf242mzGbzXUfl5SUAI7dWVar1en5aq/himsJ55B76L7sdoV9p0rZcKSAHceKOFpQQdbpSizV9rrnNPEzEB1sIsSupWxrJle2jyQ8wKRianEp5HXofux2hfWHC/j1QB5bjp4m83QFVVbHa8/PqCMxKpCe8U24tnMUrcIcrzm5f57L1a/B+lxHoyiK4sQsdRRF4YYbbuD06dOsXbv2nM+ZMWMGM2fOPOvxBQsW4Ofn5+yIQggnKKiCzblatuRpOG2p36iKBoU2wQp9IxS6hinoZFBGiHqx2WFjroZfsrUUmi/uBZQQqDCimZ3EEJeUB8ILVFRUMHHiRIqLi/+0/aXLCs8pU6bwww8/sG7duvN2tT/XiGdcXBz5+fku6+O5YsUKhg0bJr3LPJTcQ/dxNL+cuWuOsnTnSWx2xz8z/kYdvRKa0CchlLaRgcSH+RLia0Cv1VJptZFXaiY9t4Sl63eRowSz52Rp3edr1sSXewe0YFz3WPQ6mZJ3Z/I6dA87Mot46pvdHMmvACDIR8+1naMY0DqMxKhAwgOM2OwKp0rM7D5Rwsr9uaw6kIfV5ni9XtUunBnXdyAqyEfNL0NcAle/BktKSggPD7+owtMlU+0PP/wwS5cuZc2aNRc8SslkMmEynT2tZjAYXPqPl6uvJxqe3EP1lFRZeWPFIT7emFFXcPZvHcb4ns0Z3iESH4PunH8vCIgM8addVCCa4zsZObIv2SUWvkk5wScbj3H8dCXPfbePBVuO89z1HejXKtyFX5W4FPI6VIfNrvDGyoO8tSoduwKh/kb+cmVrJvRqfs7XX2igHx1imzCuVzy5JVW8veoQ/910jF8O5LNtzkbenNCNIe0iVPhKxOVy1WuwPtdwauGpKAoPP/ww33zzDb/99hsJCQnOvJwQQmXL95zi6W93k1fqmLm4MjGCR65qQ9e4kEv6fPFh/kwb2pb7B7Vi4dZM3vjlEAdySpn43mYm9W7O36/tgK/x3IWsEI1RaZWVaV+k8sv+XADGJMXy3PUdCfa9uMIgIsiHp0cmEll+hGUFoew6UcJd87fy1DWJ3DeolTOji0bCqYXnlClTWLBgAUuWLCEwMJBTp04BEBwcjK+vrzMvLYRwoQpLNc//sI8FmzMBSAj3Z+b1HRnUtmmDfH5fo447+ycwOimWV5cf5JNNx/hscyYbjxQw55ZkOsQ4fymOEO6uuMLKbR9uJu14MSa9lpfGdmF0Uuwlfa4Yf/j8xl688ONBPt+SyYvL9nO6wsrjV7eTHfDisjh1odTcuXMpLi5myJAhREdH1/23cOFCZ15WCOFCWYUVjHl7Q13Ref+glvw0bWCDFZ2/F+Jn5B+jO/Hp5N5EBpk4klfO2Lkb+Gn3yQa/lhCe5HS5hYnvbyLteDGh/kYW3t/3kovOWia9ltljOvPUNYkAzP3tMP/8aX9DxBWNmFMLT0VRzvnfnXfe6czLCiFcZGtGIaPfWs/+U6WEB5j47J7ePDWyPSa9c6e/B7QJ56dHBjGwTTiVVhsPfLqDt1al46K9kkK4lXJzNbd9uJk92SWE+Rv5/N4+dLvE5S3ncv/gVrxwYycA3l19hHdWH26wzy0aH9kaKoS4JEtSTzDxvU0UlFvoGBPE0qn96d/adRt+mvgb+ejOntzVvwUA//r5AC8u2yfFp2hUqm12Hv48hd0nHEXnF/f1oV1UYINfZ1LveP420jHy+c8f9/PjLpllEJdGCk8hRL19viWTaQtTsdoURnaO4qsH+hIT4vp123qdludGdeTZ6zoA8N7ao/ztm13Y7VJ8isZh1vd7+XV/Lia9lvfv6EGbyIYvOmvdN6gVkwc4Ngn/9audHDhV+id/Q4izSeEphKiXD9cd5anFu1AUuL1vPHNuScbP6NJD0M5y94AEXh7bBa0GPt+Sxczv9sjIp/B6X23L4r8bj6HRwBvju5HUvInTr/nUNYn0bx1GhcXGvf/dRnGFnG4k6kcKTyHERftw3VFmfb8XcGwimnl9R7Ra99jhOq5nHK+O64pGAx9vPMYryw+oHUkIpzmYU8ozS3YD8OjQtlzTOdol19XrtMy5JZlmTXzJLKzgb9/ukl/yRL1I4SmEuCiLdxyvKzr/cmVrnrwm0e3aqtyY1Ix/3ODYBPHWqsPMX39U5URCNLwKSzUPfbaDKqudgW3CmXJFa5dev4m/kbcnJaPXavgh7STfpp5w6fWFZ5PCUwjxp1buzeGxr9MAuLt/Ao8Oa+t2RWetW/vE89jV7QDH+rffDuSqnEiIhvXisn2k55YRGWTi9fHd0Kkw69ClWQiPXNUGgGe/3cPx0xUuzyA8kxSeQogL2pF5mikLdmCzK4xNbsbfr23vtkVnrYeGtGJcj2bYFZi6IEU2QQivsT49n083OXrmvjauG+EBZx8z7SoPDmlF9/gmlJqra9Z9y5S7+HNSeAohzuv46Qru++82zNV2rkqM4KWxnd1mTeeFaDQanh/dmd4JoZSZq5n88VZOl1vUjiXEZSmtsvJ4zczDbX3iXdq+7Fz0Oi2v3twVo17L2kP5LN2ZrWoe4Rmk8BRCnFOZuZp7Pt5GfpmF9tFB/PuWJPQ6z/knw6jX8s6t3YkP8+P46Uqmf5kqbZaER3vpp/2cKKokLtSXJ2tOE1Jbi3B//nKlY43pP77fS1GF/IInLsxz3kWEEC5jtytM+yK17kSi9+/ogb9J3ZZJl6KJv5G5k7pj0mtZdSCPd9bIiSvCM6UdL+KzmmNpXxrbxa1ej/cNakXriADyyyy8JEdqij8hhacQ4ixzVx9m5b4cjHot827vTqwKzeEbSoeYIGZe3xGAV5cfZPORApUTCVE/drvCM9/uRlHgxqRY+rVSd4r9j4x6LS+MdnST+GJrFnuyi1VOJNyZFJ5CiDNsOJzPqzU9MP9xQ0eSXdCU2tnG94zjxqRYbHaFv3yRItOBwqN8sTWLnceLCTTpeWqke0yx/1HvlmGM6hqDosDz38vRteL8pPAUQtTJLaniL5+nYldgbHIzxvWIUztSg3BsNupEy6b+5JSYeWbJHrUjCXFRCsstvPyzY/r60WFtiQj0UTnR+T1+dTuMei0bjxSwcp+0MRPnJoWnEAKAapudhz9PIb/MTLvIQJ4f3cnt2ybVh79Jz+vjHD0Pv9uZLTtwhUd4Y+VBiiqsJEYFcnvfeLXjXFBcqF/dWe4vLtuH1WZXOZFwR1J4CiEAmLMqnc1HC/E36nj71mR8jTq1IzW4rnEhTK055eWZb3eTU1KlciIhzu9ofjkLajYUPTuqg0d0lXhoSCvC/I0czS/n6+3H1Y4j3JD7/xQLIZwuJfM0//k1HYAXx3SmVdMAlRM5z9QrW9M5NpjiSitPLEqTtWjCbf3r5/1U2xWuaNfU7TYUnU+gj4EHh7QC4D+/HMJcbVM5kXA3UngK0chVWKqZ/uVObHaF67vGcEO3WLUjOZVBp+X18V0x6rT8diCP79JOqh1JiLOkZJ5m2a5TaDTwhJv07LxYt/aJJzLIRHZxFQu3ZqkdR7gZKTyFaOSe/2EfR/PLiQ724R83dFI7jku0jghkak3T61nf7ZFd7sKtKIrC7GWODUU3JTcjMSpI5UT142PQ1S1pmfNrOlVWGfUU/yOFpxCN2C/7curWkL1yc1eC/QwqJ3KdBwa3ok1N0+vaN3kh3MFvB/LYklGISa9l+vC2ase5JON6xhEb4ktuqZlPNx1TO45wI1J4CtFIFVdYeXLxLgAmD0hQ/dxnVzPqtcwe0xmAhduy2HhYGssL9SmKwusrDwJwR78WRAd75uENJr2Oh2tmFd5ZfVhGPUUdKTyFaKReXLaPvFIzLcP9eezqdmrHUUWPFqFM6t0cgGeX7Jb2L0J1qw7kkna8GF+DjvsGtVQ7zmUZ270ZsSG+5JdZ+Ep2uIsaUngK0QhtSM9n4TbHov+XbuqCj8H7WiddrMevTiTU38ih3DI+2ShTgkI9iqLwxspDANzeN57wAJPKiS6PQaetK57fW3OEavnFTiCFpxCNTqXFVjfFflufeHq2CFU5kbqC/Qx1I76vrzxIfplZ5USisfr9aOe9Hj7aWWtcjzhC/Y1kFlawbPcpteMINyCFpxCNzOsrD5JZWEF0sA+Pj2icU+x/NK5HHJ1igyitquaVnw+oHUc0Qt422lnL16jjzn4tAJj722Hpmyuk8BSiMUk7XsT7a48A8MKNnQj0aTy72C9Ep9UwY1RHwLHRKO14kbqBRKOz5lC+14121rq9bzx+Rh37Tpaw5lC+2nGEyqTwFKKRsNkVnlq8C7sCN3SL4crESLUjuZUeLUK5MSkWRYHnlu6RkRnhUu/8dhiACb3ivGa0s1aIn5Fbejk28dV+naLxksJTiEbis83H2JNdQpCPnmeu66B2HLf05DWJ+Bl1pGQWsWyXrEcTrrEzq4iNRwrQazXcM9C7RjtrTR6QgE6rYeORAvadLFE7jlCRFJ5CNAL5ZWb+VbN28bERiV43otJQIoN8uLfmjf9fP++X9krCJd5d4xgFvL5rDLEhntm388/EhPgyomMUAPPXZ6gbRqhKCk8hGoHZy/ZTWlVN59hgJtZMeYlzu3dQS8IDjGQUVPD5lky14wgvdzS/nB9rdnvfP7iVymmc667+LQD4NvUEheVyTG1jJYWnEF5ua0Yhi3YcR6OBf4zuhE6rUTuSWwsw6XlkqOOYwjdXHqK0yqpyIuHN5q05gqLAlYkRtIsKVDuOU3WPb0Ln2GDM1Xb5pa4Rk8JTCC9WbbPzzLe7AZjQM45ucSHqBvIQE3rG0TLcn4JyC++tOaJ2HOGlckurWLTDcaLPA14+2gmg0WjqRj0/2XhMlrI0UlJ4CuHFPtl0jP2nSgnxM/DY1Ylqx/EYBp22rqn8e2uPkltSpXIi4Y3+u+EYlmo7Sc1D6NmiidpxXOLaLtGEB5g4VVJVt8RANC5OLTzXrFnDqFGjiImJQaPR8O233zrzckKI3zldbuH1FQcBeOzqdoT6G1VO5FlGdIoiqXkIlVYbb61KVzuO8DJVVhsLaqab7xvYEo2mcSyBMel1TOrtWGf+qRxR2yg5tfAsLy+na9euzJkzx5mXEUKcwxsrD1JSVU376CAm9JQNRfWl0WjqRj0/35JFdlGlyomEN1m6M5vCcguxIb4M69C4eure0qs5Oq2GLRmFHMopVTuOcDGnFp7XXHMNzz//PGPGjHHmZYQQf5CeW8qnmx2jKc9c1142FF2ifq3C6dMyFIvNzhwZ9RQNRFEUPqppKXRb33j0usa16i0q2IerEiMA+GyzbDJqbBrXT7sQjcTzP+zDZlcY1iGSfq3C1Y7j0aYPc4x6frk1i6zCCpXTCG+w+Wgh+06W4GPQMqFnnNpxVDGxZrp98Y7jVFpsKqcRrqRXO8Dvmc1mzGZz3cclJY7TDaxWK1ar81ua1F7DFdcSziH3ENYeyue3A3kYdBoeH97a474X7nYPk5oF0r9VGOsPF/DmyoPMvrGj2pHcnrvdQ3fz4TpHp4TR3WLwN2jc7vvkivvXt0UIzUJ8OF5UxZKULMYmxzrtWo2Rq1+D9bmORnHRgcQajYZvvvmG0aNHn/c5M2bMYObMmWc9vmDBAvz8/JyYTgjvYFPg5Z06TlVqGBJt58YW0q6kIWSUwuu79WhR+Fs3G02983AZ4QIFVfCPFB0KGp7sWk10I35rW3FCw/eZOloEKDzaWUY9PVlFRQUTJ06kuLiYoKCgCz7XrUY8n3rqKaZPn173cUlJCXFxcQwfPvxPv5CGYLVaWbFiBcOGDcNgMDj9eqLhNfZ7+NnmTE5V7qeJn4FX7hpAsK/nfQ/c9R7usOxg9cF89hDHyyM7qx3HrbnrPXQHL/18EIUM+rUKZfJNPdSOc06uun+9ysz89K81ZJRBQtJA2kd7dwN9V3L1a7B2hvpiuFXhaTKZMJnOPkPaYDC49B8vV19PNLzGeA+LK638e5VjCu/RYW0JD/LsoRR3u4d/Hd6O1QfzWbLzJFOvakurpgFqR3J77nYP1VZhqebLbY6G8Xf3b+n23xtn37/oJgau7hjFD7tO8uWOEzw/Wn6ha2iueg3W5xpO3VxUVlZGamoqqampABw9epTU1FQyM2UXmxAN7d3Vhykst9Cqqb+cx+4EXZqFMLR9JHYF3l51WO04wgN9m5JNSVU18WF+XFmzq7uxq91k9G1KNhWWapXTCFdwauG5bds2kpKSSEpKAmD69OkkJSXx7LPPOvOyQjQ6OSVVfLj+KABPXtO+0bVncZWHr2wNwLepJ2SHu6gXRVH4ZJOjYfptfeLRSoszAPq1CqN5qB9l5mp+kpOMGgWnvjsNGTIERVHO+m/+/PnOvKwQjc4bKw9RZbXTI74JQ9vLSIqzdI0LYWCbcGx2hXlyhruoh9SsIvadLMGo13JT92Zqx3EbGo2m7vvxVc0yBOHdZFhECA+XnlvGl9uyAHjimsRGc/SeWh4a4hj1XLgtS85wFxettlH6dV2iCfGT42t/b2z3Zmg0sPFIgcwkNAJSeArh4V75+QA2u8LQ9pH0bBGqdhyv16dlKN3jm2CptvPBuqNqxxEeoLjCync7swGY1Dte5TTuJzbEl/41B118vV1GPb2dFJ5CeLAdmaf5ac8ptBp4fEQ7teM0ChqNhilXtALg003HKKqwqJxIuLtFO45jrrbTPjqI5OYhasdxSzf3cEy3f739OHa7S9qLC5VI4SmEh1IUhX/+uB+AscnNaBspPfBc5Yp2EbSPDqLcYqs7c1uIc1EUhc82OzYVTerdXJbCnMfVHaMI9NFzoqiSTUcK1I4jnEgKTyE81G8H8thytBCjXsujw9qqHadR+f2o5/wNGZSZpQ2MOLdNRwo5nFeOv1HH6CQ5FvJ8fAw6RnWNAeArmW73alJ4CuGBbHaFl35yjHbe1a8FMSFyhqOrXdMpmpZN/SmutPJZTZscIf6odrTzhqRYAkxudWaL27m5Znf7j7tPUlLlXufXi4YjhacQHujblBPsP1VKkI+eB4e0UjtOo6TTanhwsON7//66o5ir5axpcaa8UjM/73H0ppRDHf5ct7gQWkcEUGW180PaSbXjCCeRwlMID2OptvPaioMAPDiktbRmUdEN3WKJCvIhr9TMktRsteMIN/PV9iysNoVucSF0ig1WO47b02g0daOetS3ihPeRwlMID/PltixOFFXSNNDEnf1aqB2nUTPqtdw9oAUA7605IrtxRR27XeHzLY7enbf2kRZKF+vG5Fi0GkjJLCIjv1ztOMIJpPAUwoNUWW28tSodgClDWuFr1KmcSEzo1ZwAk55DuWWsPpindhzhJhzN0CsJ9NFzbedoteN4jIhAHwa0aQo4jqYV3kcKTyE8yBdbMjlZXEV0sA8TZM2YWwjyMXBLrzgAOUZT1KmdKr6hW4z8glhPNyY5drd/m3ICRZFZBG8jhacQHqLKauOt3w4DMOWK1vgY5M3MXdzVPwG9VsPGIwXsOl6sdhyhsuIKKz/udmwqGtcjTuU0nmd4hyh8DToyCipIzSpSO45oYFJ4CuEhPt10jLxSM7EhvvJm5mZiQnzrehDOWyujno3d0rRsLNV2EqMC6SybiurN36Tn6o6RgGPUU3gXKTyF8AAVlmreWe0Y7fzLVa0x6uWl627uGZgAwLJdJ8kqrFA5jVDTVzXT7Df3iJOTii5RbbP979JOYrXZVU4jGpK8ewnhAf678Rj5ZRaah/oxJrmZ2nHEOXSMCWZA63BsdoUP1x9VO45Qyb6TJaQdL8ag0zC6W4zacTzWgNbhhAcYKSy3sPaQbNrzJlJ4CuHmyszVvFs32tkGg05etu7q3kEtAVi4NYviCjl5pTGq3VQ0tH0kYQEmldN4Lr1OW7d85ZsU6ZHrTeQdTAg3N3/9UU5XWGkZ7i8jKG5uUJtwEqMCqbDY+GyLHKPZ2JirbXVrEsf1lHXYl+vGmun25XtOUSpHaHoNKTyFcGMlVda6Fj2PDG2DXkY73ZpGo+HegY5Rz4/WZ2CplrVpjckv+3I5XWElKsiHQTW9KMWl6xwbTMum/pir7fy8J0ftOKKByLuYEG7sg7VHKamqpk1EANd1kdFOTzCqawwRgSbySs38sEumCBuT2mn2sd1j0WllU9Hl0mg03NjNMeopu9u9hxSeQripogoLH65zbFKZNrStvJF5CKNey201RyR+tD5DGmA3EieLK1lTc3LVzd1lmr2h3FBTeK4/nE9OSZXKaURDkMJTCDf1/tqjlJqrSYwK5JpOUWrHEfUwsXdzjHotaceL2X7stNpxhAss3nECuwK9EkJpEe6vdhyv0TzMjx7xTVAUWJoqMwjeQApPIdxQYbmFj2pa8jw6rC1aGe30KGEBpropwo/WZ6gbRjidoih10+zj5XCHBndDzabK79Kk8PQGUngK4YbeXXOYcouNTrFBDO8QqXYccQnuGtACgB93n+T4aWko7802Hy3kWEEFASY913SW2YmGdk3naHRaDWnHizmaX652HHGZpPAUws3klZr57wZHK57pw9rKySceKjEqiH6twrAr8MlGaa3kzWpHO0d1jcbPqFc5jfcJDzDRr1UYAN/vlFFPTyeFpxBu5p3Vh6m02ugaF8IV7SLUjiMuw939Hcdofr4lkwpLtcpphDOUVllZtusk4DgiUzjH9TXN5JfuzJYNex5OCk8h3EhOSRWfbpLRTm9xZWIE8WF+lFRVs2iHtIPxRt+nnaTKaqd1RABJcSFqx/FaV3eKwqjTcii3jP2nStWOIy6DFJ5CuJG3V6VjrrbTI74Jg9qEqx1HXCatVsOd/VoA8NH6o9jtMlLjbWqn2cf1aCa/KDpRkI+BIe0cTfm/k+l2jyaFpxBuIruoks+3ON7EZLTTe9zcI45Ak54jeeWsOZSndhzRgA7llJKSWYReq+HGpGZqx/F61/9ud7tMt3suKTyFcBNvrUrHYrPTp2Uo/VrLaKe3CDDp687t/lBaK3mV2tHOKxMjaBpoUjmN97sqMRI/o46swkpSs4rUjiMukRSeQriBrMKKujexR4e2VTmNaGh39G2BRgNrDuaRnivr07yB1WZncc263XGyqcglfI06htW0l1sq0+0eSwpPIdzAnF/TsdoUBrQOp3fLMLXjiAbWPMyPYe0db5jSUN47/Lo/l4JyC00DTXVrD4Xz1e5u/yHtJDZZM+2RpPAUQmUZ+eV8veM44DilSHinuwc4Wist2nGcogqLymnE5fqqZoZiTHIsep28lbrKwDZNCfY1kFtqZvPRArXjiEsgrxYhVPbvXw9hsysMadeU7vFN1I4jnKR3Qijto4OostrrNpEJz5RbUsWqA46NYjd3l2l2VzLqtVzTyXE6lOxu90wuKTzffvttEhIS8PHxoXv37qxdu9YVlxXC7R3OK+PbFMc6MVnb6d00Gg13928BwCcbM6i22dUNJC7Z4pQT2OwKPeKb0DoiQO04jU7tdPuPu09hqZbXkadxeuG5cOFCpk2bxtNPP01KSgoDBw7kmmuuITMz09mXFsLtvbnyEHYFhraPpKs0n/Z6o7rGEOZvJLu4iuV7c9SOIy6Boii/690po51q6N0yjKaBJooqrKxLlxZlnsbphedrr73G5MmTueeee2jfvj1vvPEGcXFxzJ0719mXFsKtHcwp5bs0x1TRo8PaqJxGuIKPQcek3s0BR0N54Xm2HzvNkbxy/Iw6RnaJVjtOo6TTari2s+N7/93OkyqnEfWld+Ynt1gsbN++nSeffPKMx4cPH86GDRvOer7ZbMZsNtd9XFJSAoDVasVqtTozKgApxwqYs0dLUt9yopv4O/16ouHV/py44uflcr2+/ACKAld3iKBtUz+PyOwKnnQPL8X4HrG8/dthtmacJvVYAR1jgtSO1OC8+R5+scUxW3dNp0hMWsUrv0ZPuH8jO0Uwf0MGy/ecoqS8Cl+jTu1IbsXV97A+13Fq4Zmfn4/NZiMyMvKMxyMjIzl16tRZz589ezYzZ8486/Hly5fj5+fntJwAigKv7dKRWa7lkY/XcmdbWTfiyVasWKF2hAs6UQ4/7tGjQaGbIZtly2SR/B+5+z28HF1DtWzP1/LC1xu4tbX3/lvjbffQbIOlqTpAQ6w5k2XLvHvJmDvfP0WBUJOOQrON1xYuJylMWiudi6vuYUVFxUU/16mFZ60/Hv2nKMo5jwN86qmnmD59et3HJSUlxMXFMXz4cIKCnD8q0KxzIePe20pKgZZ7WyQzrEOE068pGpbVamXFihUMGzYMg8GgdpzzemhBKpDLyE7R3HNTF7XjuBVPuYeXI/Z4MTe9u5nUQh3/HnQF4QHedeqNt97Dr3ecwLJlDwlhfkwd399rj7X1lPu3z3CId9ceJVsfzdMju6kdx624+h7WzlBfDKcWnuHh4eh0urNGN3Nzc88aBQUwmUyYTGf/A2wwGFzyjevaPJQrYxVWntAw4/t99G8bQbCv+77oxPm56mfmUuw6XsyKfbloNfDo8HZum1Nt7nwPL1ePhHCSmoeQklnEwu3ZTPPSjgbedg8XpzhmJsb1bI7RaFQ5jfO5+/27IakZ7649ym8H86m0QZCP+2ZVi6vuYX2u4dTNRUajke7du5811LtixQr69evnzEtfshHN7CSE+ZFbamb2sn1qxxFe6PWVBwG4oVustGJpxO7q72go/+mmTGkJ4wGO5JWxNeM0Oq2GscmxascRQPvoQFpHBGCptrNij3SJ8BRO39U+ffp03n//fT788EP27dvHo48+SmZmJg888ICzL31JDFp4YXRHAL7YmsWG9HyVEwlvsiPzNL/uz0Wn1fCXq2Qne2N2TacoIoNM5JeZ+WGXrPF1d19td5wuNqRtUyKCfFROI8CxjK+2p6ec3e45nF54jh8/njfeeINZs2bRrVs31qxZw7Jly4iPj3f2pS9ZzxZNuK2PI9+Ti3dRabGpnEh4i9eWO0Y7xybHkhAunRMaM4NOW/fvzEfrM1AU2RzhrqptdhbVFJ43S+9OtzKqpvBcl55PQZn5T54t3IFLTi566KGHyMjIwGw2s337dgYNGuSKy16Wx0e0IybYh8zCCl5bcUDtOMILbD5SwLr0fAw6DQ9fKaOdAm7p1RyjXkva8WJ2ZJ5WO444j98O5JFbaibM38iVibLp1J0khPvTpVkwNrvCsl3S09MTyFnt5xHoY+CFGzsD8MG6o6RmFakbSHg0RVF4tWa0c3zPOOJCndseTHiGsAATo7s5Rmw+XJ+hbhhxXgtrTioakxyLUS9vm+5Gpts9i7yCLuCKxAhGd4vBrsATX6fJBgBxydal57MloxCjXsvUK2S0U/xP7Sajn3afIruoUuU04o9yS6v4dX8u4PilUbif67rEoNHA1ozTnJDXkNuTwvNPPDuqI6H+Rg7klPLO6sNqxxEe6PejnZN6NycqWDYmiP9pHx1En5ah2OwKn246pnYc8QeLd5zAZldIbh5C64hAteOIc4gK9qF3QigA38uop9uTwvNPhPobmXG9Y5f7f349xKGcUpUTCU+z6kAuqVlF+Bp0PDikldpxhBuqHfX8fEsmVVbZzOguFEXhy62OaXYZ7XRv13d1tLhakiqFp7uTwvMijOoSzdD2EVhtCo8vSsNml92n4uL8frTz9n7xRATKaKc429D2kTRr4svpCivfppxQO46ose3YaY7kl+Nn1HFtlxi144gLuKZTFHqthr0nS0jPLVM7jrgAKTwvgkaj4R+jOxFg0pOSWcTHGzLUjiQ8xM97TrEnuwR/o477B8lopzg3nVbDHX1bANJayZ0srBntvK5LNAEml5wwLS5RE38jg9o2BWSTkbuTwvMiRQf78tTIRAD+9fMBsgorVE4k3J3drvD6ikMATB6QQKi/9x+xJy7duJ5x+Bl1HMgpZePhArXjNHqlVVZ+SHO055Fpds9wQ02HiKWpJ+SXNzcmhWc93NKzOb0TQqm02vjbN7vkB1tc0Pe7TnIgp5QgHz2TB7ZUO45wc8G+BsYmNwOktZI7+D7tJJVWGy2b+pPcvInaccRFGNo+Eh+DloyCCnadKFY7jjgPKTzrQavV8M+xXTDptaw9lM/XNSdZCPFH1TY7b6xwrO28d2BLgn0NKicSnuDO/i0A+GV/DpkFMquiptpp9vE94tBoNCqnERfD36RnaPtIAJbKJiO3JYVnPSWE+zN9WFsA/vH9XnJLq1ROJNzRt6nZHMkvp4mfgbsGJKgdR3iIVk0DGNy2KYoCH2/MUDtOo3Uwp5TUrCL0Wg1jakahhWeobSb/XVq2bAR2U1J4XoLJAxLoHBtMSVU1zy3Zo3Yc4WasNjv//sWxtvOBwa1kU4Kol9pRzy+3ZlFmrlY3TCNV20LpysQImgaaVE4j6mNwu6YE+ejJKTGz5Wih2nHEOUjheQn0Oi0vje2CXqvhx92n+Gm3nA8r/ueLrVlkFlYQHmDi9pqdykJcrMFtmtIy3J9SczWLZDmPy1mq7SyuaWklm4o8j0mv45pO0YDsbndXUnheog4xQTww2NEe55kleyiusKqcSLiDCkt13WjnI1e1xteoUzmR8DRaraZu1HP+hgzsMl3oUr/sy6Gw3EJEoInBNe15hGe5vmZ3+4+7T8pR125ICs/LMPXK1rRq6k9eqZkXlu1VO45wAx+tzyCv1EzzUD/G92yudhzhocYmNyPQR8/R/HJWH8xTO06jsnCbY5r9pu7N0OvkLdIT9WkZRtNAE0UVVtYektePu5FX1WXwMeh4aWwXNBr4cttx1h3KVzuSUNHpcgvv/HYYgL8Ob4tRLy8vcWn8TXrG93BM8364/qjKaRqP7KJK1tQU+uN6yDS7p9JpNVzXRabb3ZW8M16mHi1Cub1PPABPLk6jwiKbARqruasPU2qupn10EKPkeD1xme7o1wKNBtYeyic9t1TtOI3CF1uzsCvQp2UoLcL91Y4jLkPt7vYVe3OotNhUTiN+TwrPBvDYiERiQ3w5frqSV34+qHYcoYKTxZXMrzlK9fER7dBqpe+fuDxxoX51PQnnyzG9Tldts7NwayYAE3vHq5xGXK5ucSE0D/WjwmJj5b4cteOI35HCswEEmPS8cGMnAD7acJQdmadVTiRc7c2Vh7BU2+mVEMoQ2ZAgGshdNZuMFm0/IRsYnezX/bnklJgJ8zdydcdIteOIy6TRaBjV1THdvkSaybsVKTwbyJB2EYxJjkVR4Imv0zBXy9B+Y5GeW8aXNRsSnhiRKKeciAbTt2UYiVGBVFptLNyWqXYcr/bZZsf396YezTDppRuFN7ihWywAqw/myi9ubkQKzwb0zLUdCPM3cii3jLdXHVY7jnCRV5cfwK7AsA6RdI+XM51Fw9FoNHWjnh9vOEa1TVrDOENWYQVranY/3yLdKLxG28hAEqMCsdoUftoj/bbdhRSeDaiJv5GZN3QE4O3f0tl/qkTlRMLZdmYV8ePuU2g08NjV7dSOI7zQDd1iaeJn4ERRpaxVc5IvtmaiKDCgdbhsKvIyo2o2GX1TcyiAUJ8Ung3s2s7RDOsQidWm8MTXaXJWrBdTFIV//rgfgDFJzWgbGahyIuGNfAw6JvZ2jMJ9uC5D3TBeyGqz8+U2xwlRk3rLaKe3GZ3kmG7fdKSQ46crVE4jQArPBqfRaHh+dCcCffTsPF7MR9KDz2v9si+XjUcKMOq1PDqsjdpxhBe7rU8LDDoNWzIKSZHNiw1q5d4c8krNNA00MbSDbCryNrEhvvRtGQbIJiN3IYWnE0QG+fC3ke0BeGX5ATIL5Lcsb2O12Xnxx30ATB6QQLMmfionEt4sKtiH67s6Rm7mrTmichrvUrupaFyPZhjkpCKvNCbZ8dpZtP04iiKzkGqTV5mTTOgZR9+WYVRZ7Ty5OE1+2L3M51syOZJXTpi/kYeGtFI7jmgE7hvUEoCf9pwiI79c5TTeISO/nHXp+Wg0MEE2FXmtazpH42vQcSS/nNSsIrXjNHpSeDqJRqNh9pjO+Bi0bDhcwIIt0grFW5RUWXlj5SEApg1tQ6CPQeVEojFoFxXIFe2aoijw/joZ9WwIn9c0jB/UpilxoTJr4a0CTHpGdIoCYPEO2WSkNik8nahFuD//N9yx0/mFH/bJlLuXeGtVOoXlFlo19eeWXjJKIlznvkGO0fWvth2noMyschrPVmW1sXCro//urX3kpCJvVzvd/l1atvTZVpkUnk52d/8EeiWEUmGx8X9f7ZRd7h4uq7CCj2p2Fv9tZHv0siZMuFCflqF0aRaMudrOxxuPqR3Hoy1JPUFRhZVmTXy5MjFC7TjCyfq1CicyyERRhZVV+/PUjtOoybumk2m1Gl69uSt+Rh1bMgpll7uHe/nnA1hsdvq1CpM3K+FyGo2mbq3nJxszqLTIyM2lUBSF+RschfttfeLRaeW0MW+n02rqWist3nFc5TSNmxSeLhAX6sffr+0AOAqXQzmlKicSlyIl8zTf7cxGo4Gnr20vR2MKVYzoGEVcqC+nK6x8tT1L7TgeaWvGafadLMHHoGV8zzi14wgXGZPUDIBVB3IpLLeonKbxksLTRW7pFcfgtk2xVNv561c7scrRdx7Fblf4x/d7ARib3IyOMcEqJxKNlV6n5d6BjlHP99celWM0L8HHGzIAGN0tlhA/o7phhMu0iwqkU2wQVpvC92nS01MtUni6iEaj4aWxXQjy0ZN2vFjOcvcw36ScYEdmEX5GXd2GMSHUcnP3OJr4GcgsrOCnPafUjuNRThZX1n3P7ujXQt0wwuVqRz0Xye521Ujh6UJRwT78Y3QnAP7z6yF2nyhWOZG4GKVVVmbXHI059crWRAX7qJxINHa+Rh239W0BwLurj0if4Hr4bFMmNrtCr4RQ2kcHqR1HuNj13WLQazXszCqSZW8qcWrh+cILL9CvXz/8/PwICQlx5qU8xvVdY7imUxTVdoXpX6ZSZZXNAe7u378cIr/MTEK4P5MHJKgdRwgA7ugbj49By64Txaw5lK92HI9QZbXxeU1P5TtltLNRCg8wcUXNxtDadlrCtZxaeFosFm6++WYefPBBZ17Go9Se5R4eYORgThmvLj+gdiRxAem5pXy0PgOAZ0d1wKTXqRtIiBphASYm9nL0n3zr13SV03iGH9JOUlBuITrYh+FyLnujdUsvx4ayRTuOS09PFTi18Jw5cyaPPvoonTt3duZlPE5YgInZY7oA8N7ao6w9JD3F3JGiKMxYupdqu8LQ9hFc0U7aJwn3ct+glhh1WrZkFLL5SIHacdyaoih8sM7Rzu7WPvHSg7cRG9w2gqggH05XWFm+J0ftOI2OXu0Av2c2mzGb/3caR0lJCQBWqxWr1er069dewxXXGtImlIm9mrFgy3H++uVOvpvSl1B/2V15uRryHv68J4d16fkY9VqeHNHWJT8XwrWvQ08X5qdjbHIMn289zr9/OcT8OPdYs+iO93DD4QL2nizB16BlXHKMW2VzN+54/xra2OQY3vrtCJ9vOcaIDk3VjtPgXH0P63MdjeKCVenz589n2rRpFBUVXfB5M2bMYObMmWc9vmDBAvz8vO8cXYsNXtmlI6dSQ6cmdu5pZ0daQ7oHiw1m79RRaNYwPNbOtc2lZY1wTwVV8HyKDjsapneqJj5Q7UTu6Z19WvYVaRkYZeemBHk9N3YFVfCPFB0KGp5JqiZc9oxeloqKCiZOnEhxcTFBQRf+BbjeI57nKw5/b+vWrfTo0aO+n5qnnnqK6dOn131cUlJCXFwcw4cP/9MvpCFYrVZWrFjBsGHDMBgMTr8eQLsepYx9dxO7T2spbtqRib2kmfHlaKh7+NqKQxSajxId7MO/7u6Hn9GtJge8mhqvQ0+3m90sTslmZ3U0D45MUjuO293Dgzml7Nu4Ea0GnrtlEPGh3jeQ0ZDc7f45y6+l21mXXkB+YBtuH9ZG7TgNytX3sHaG+mLU+9106tSpTJgw4YLPadGiRX0/LQAmkwmTyXTW4waDwaU//K68XpfmoTwxIpHnf9jH7J8O0K91U9pEypDF5bqce3jgVCnv1ZzH/tyojgT7+zZgMnGxXP2692RTrmzDN6nZ/LI/j0N5lXSIcY8pd3e5hx9tdOxeHtEpitaRcvjDxXKX++csE3vHsy69gEUp2fz16kSvXPfrqntYn2vUu/AMDw8nPDy8vn9NXMDd/RNYfTCPtYfy+csXqXw7pZ/snlaJ3a7wt292UW1XGNYhkhGdotSOJMSfatU0gGs7R/N92kneWpXOW5OS1Y7kNnJKqliS6mgWXnvikxAAQ9tHEuZvJLfUzKoDeQyTTgcu4dTyPjMzk9TUVDIzM7HZbKSmppKamkpZWZkzL+txtFoNr47rSpi/kX0nS3j5J2mxpJbPt2ay/dhp/I06Zl7fUe04Qly0qVe2BuCHXSfZd/Lip7283fwNGVhtCj1bNCGpeRO14wg3YtRrGdvdcZLRFzX9XYXzObXwfPbZZ0lKSuK5556jrKyMpKQkkpKS2LZtmzMv65EiAn14+SZHi6UP1h3ll33S4sHVckuq+GfNCUX/d3U7YkJkil14jsSoIK7rEg3Aq8sPqpzGPZRUWfl00zFARjvFuY3v6dhXsepALtlFlSqnaRycWnjOnz8fRVHO+m/IkCHOvKzHuqp9ZN1pGtO/3Mnx0xXqBmpkZn2/l9Kqaro0C+b2muMIhfAk04a2RauBlftySM0qUjuO6j7ZeIzSqmpaNfVnaHuZRhVna9U0gD4tQ7Er1P2SIpzL+1bSeri/jWxP17gQiiutTFmQgqVa2n64wsq9OXyfdhKtBl68sTM6rfS1Ep6ndUQAY5IdU4eN/VS0Cks17689AjiWIWjlNS3Oo3bA54utWXKMtQtI4elmjHotc25JItjXwM6sIl5ctk/tSF7vdLmFp77ZBTim4zrFyq5X4bkeuaoNeq2GtYfyG/VpRp9tyuR0hZX4MD9GdYlRO45wY0PbRxIT7ENhuYXvdmarHcfrSeHphuJC/XhtXFfAsTB+2a6TKifybjO+20NeqZlWTf15dFhbteMIcVniQv3q1q29uvwgLjgjxO1UWW3MqxntnDKktVe2yRENR6/TcmvfeAA+3pjRKF8zriSvRjd1VftIHhjcCoDHv04jI79c5UTe6afdJ1mSmo1WA6+O64aPQdpYCc839crWGPWOM9zXHspXO47LLdyaRV6pmdgQX0YnxaodR3iACT2bY9Rr2X2ihB2ZRWrH8WpSeLqx/xvell4tQikzV/PgZzuotMjak4ZUUGbm6W92A/DA4FZ0iwtRN5AQDSQ62Jfb+jhGcF7+eT92e+MZwbFU23ln9WEAHhjSCqNe3ubEnwv1N3JDV8eSjI83ZKgbxsvJK9KN6XVa/n1LUl1/zycXp8kUQAN6dukeCsottI0M4JGh3nVcmhAPDWlFgEnP7hMlLNl5Qu04LvPltixOFlcREWji5poejUJcjDtqNhkt23WS3JIqdcNcpkcW7mTRUS05bvh1SOHp5qKCfZgzMRmdVsOS1Gzeq1m3JC7PNynH+SHtJDqthldv7iYnRQmvExZg4qErHMt1/vXTgUaxW7fSYuPfvxwCYMoVrWXpjKiXTrHB9IhvQrVd4RMPbq10KKeUZbtzWHtKg9kNO+NI4ekB+rYK47lRHQD454/7WX0wT+VEnu1YQTl/r5lif/jK1nRuJrvYhXe6u38CsSG+ZBdX8cG6o2rHcbqPN2aQW2qmWRNfbunVXO04wgNNHpAAwH83HqPcXK1ymkvz4XrHa71zqELzUD+V05xNCk8PcVufeMb3iMOuwMMLdnBUNhtdEqvNzl++SKXcYqNniyZMvaK12pGEcBofg47Hrm4HwNzfDpNfZlY5kfMUV1qZ+5tjbeejQ9vK2k5xSYZ3jCIh3J/iSiufe+AxmgVlZhbtcCytGRLtfqOdIIWnx9BoNMwa3ZHk5iGUVFVz73+3UVplVTuWx3ltxUF2ZhUR5KPnjQlJ0mZFeL3ru8bQpVkwZeZqXl/hvUdpvrfmCMWVVtpEBMhOdnHJdFoN9w1yHK/6wbqjHneIy6ebMrFU2+kSG0TLQLXTnJu863oQk17HO7d2JyrIh/TcMh75IpVqm2e9KNS0Pj2/brfrP8d2IVbOYheNgFar4emR7QHHySz7T5WonKjh5ZZU1U0v/nV4Ozl5TFyWG5NiaRpo4mRxFUs9qKF8ldXGJ5syALizXzwaN30ZSOHpYSKCfHj3tu6Y9Fp+3Z/LjO/2yE73i3CyuJK/fJ6CosCEnnGM7BytdiQhXKZ3yzBGdo7CZld49lvv+zfjXz8foMJio1tcCFd3lDPZxeXxMei4u79jree7qw97TDuypTuzyS+zEB3swwg3fh1I4emBusaF8OaEJDQax7D6vDWy0/1CzNU2Hvx0BwXlFjpEBzHj+o5qRxLC5f5+bQd8DTq2ZBSyeIf3tFfadbyYr3ccB+DZUR3QuOswj/Aok/o0J9Ck51BuGSv25agd50/Z7Qrv1dQCd/ZrgcGNl5G5bzJxQSM6RfHMtY6d7rN/3C/ny17A89/vI7VmXec7t3aXFiuiUYoJ8eUvVzn61c7+cR/FlZ6/RlxRFGZ9vwdFgdHdYkhu3kTtSMJLBPkYuL2f4xCG11ccdPtRz5/3nOJQbhmBJj0T3LyjgxSeHuzuAQnc1b8FAH/9cidbjhaqG8gNLdp+vK4f25sTkmge5n6tJYRwlckDEmjV1J/8MguvLT+gdpzL9sOuk2zNOI2PQcvjIxLVjiO8zL0DWxJo0rP/VCnLdp9UO855KYrCf35NB+DO/i0I9jWonOjCpPD0cH+/tgNXd4zEYrMzef5Wdp8oVjuS29iReZqnvtkFwCNXteGKxAiVEwmhLqNeyz9u6ATAfzcdY/ux0yonunTl5mpmL9sPOI68jZHNgqKBhfgZmTzQsdbzjZWHsLnpqOcv+3LZe7IEf+P/1qa6Myk8PZxOq+GN8Un0Sgil1FzN7R9uIT23VO1Yqss6XcG9H2/DUm1naPsIHrlKjsQUAqBf63DGJMeiKPD41zs99kSj11cc5ERRJbEhvtw/qJXacYSXuntAAsG+BtJzy1jqhkfPKorCf1Y5Rjtv7RtPE3+jyon+nBSeXsDXqOODO3rQpVkwheUWJr2/mazCCrVjqaaiGu79JKVuM9GbE5LQSnsVIeo8e10HmgaaOJxXzps1R0x6kl3Hi+vaJz0/uhO+Rlm3LZwjyMdQ19fzzZWHsLpZC8MVe3PYmVWEj0HLvQNbqh3nokjh6SUCfQx8fFcv2kYGkFNiZuL7mzhVXKV2LJez2uzMP6jlcF45kUEmPrizB/4mvdqxhHArIX5GXhjtmHKft+YIaceL1A1UD9U2O099k4ZdgVFdY2QJjXC6O/u1IMzfSEZBBZ+50Rnu1TY7L//sWKs9eUAC4QEmlRNdHCk8vUgTfyOfTu5NfJgfWYWVTJi3keyiSrVjuYzdrvDE4t0cKNbia9DywR09iQ6WdV9CnMvwjlGM6hqDza4w/cudVFo8Y8r9o/UZ7D5RQpCPnmev66B2HNEI+Jv0TB/eFoDXVx6iqMKiciKHr7cfJz23jBA/A/cP9pzlJlJ4epmIIB8+u6c3caG+ZBRUMO7djY1i2l1RFJ5dupvv0k6h1Si8Ob4rnWKD1Y4lhFubeX1HIgJNpOeWMev7vWrH+VMHc0p5pWY3/tPXtqdpoGeM8AjPN75HHO0iAymutLrF8pRKi43XVzqOwJ16RWuCfNx7J/vvSeHphZo18WPhfX1JCPfn+OlKxr27kSN5ZWrHchpFUfjXzwf4dFMmGg3c2trOFe2aqh1LCLcX6m/kjfHd0Gjg8y2Z/JDmvi1jzNU2pn2RirnazpB2TRnXI07tSKIR0eu0/P06x9Gzn2w8pvom3rmrD5NTYiY2xJfb+sarmqW+pPD0UjEhviy8rw+tIwI4WVzFuHc3etQ6roulKAqvLD/A2785zmCfOao93cPds+WFEO6oX+twHhrimKZ7cnGa286QvLbiIHtPlhDqb+Tlm7rICUXC5Qa2acrQ9hFU2xX+tni3ak3lj+aX807Ne97fRrbHpPeszXVSeHqxiCAfvrivD+2jg8gvszBh3iZW7c9VO1aDURSFF5ft461Vjhfg369tzy09ZRREiPqaNrQtSc1DKK2q5sHPtrvdes9f9+fw7mrHcYCzx3QmItBH5USisZpxfUf8jI6jZz/fmuny6yuKwnNL92Cx2RnYJpyRnaNcnuFySeHp5cIDTHx5fx8GtgmnwmLjnv9u4/Mtrn+xNDSbXeHZJXt4b62jpcrM6ztyj4e0khDC3Rh0Wv5zSxKh/kZ2nyjhycVpKIp7zBxkFlQw7YtUAG7rE8/VHT3vjVZ4j2ZN/Pi/4e0A+Oey/eSUuLZ7zLJdp1hzMA+jTsusGzp55Mi/FJ6NQKCPgQ/v7MnY5GbY7ApPLd7FrO/2ul0/sotVabHx4Kfb+WTTMTQaePHGztzRr4XasYTwaM2a+PH2pGT0Wg1LUrN5p2aEUU2VFhsPfLqdkqpqkpqH8IzsYhdu4I5+LegaF0KpuZonFqW5bMo9t7SKZ5bsBuCBwS1JCPd3yXUbmhSejYRBp+WVm7vUneDz4fqjTHp/M7mlntXrM7/MzC3vbWL53hyMOi1vTkhiYu/mascSwiv0aRnGc6Mcxd1LP+1n8Y7jqmWx2RUe/jyFvSdLCPM38vakZIx6ecsS6tNpNbw8tgsmvZbfDuTVHWbgTIqi8MTXaRSWW2gfHcSUK1s7/ZrOIq/iRkSj0fDosLa8e1t3Akx6thwtZNR/1rE1o1DtaBclNauIG+asJzWriBA/A5/d25vru8aoHUsIr3Jrn3gmD3Cc9/z412msOuD6deGOdWy7WbkvB6Ney7u3dZeevMKttIsK5O/X/e+XtF3Hi516vc82Z7LqQB5GvZY3J3TzuA1FvyeFZyN0dccolkztT5sIxylH49/dyL9+3o+l2j2n3hVF4ZONGdz8zgZOFFWSEO7Pogf70bNFqNrRhPA6Go2Gp0e258akWKrtCg9+up21h/Jcdn1FUXh95aG69mhvju9GD3mtCzd0a+/mXN0xEqtN4aEF2yksd05j+dSsoro+u0+MSKRtZKBTruMqUng2Uq2aBvDtlP6MSY7FrsBbqw5z49vr2Ztdona0M+SVmnnw0x08s2QPVpvCiJqiuVXTALWjCeG1tFoNL9/UhasSI6iy2pn88TZ+2Zfj9OvWtkf7d02D7mev68A1naOdfl0hLoVGo+GlsV2IC/Ulq7CSBz7Zjrm6YTtC5JZUcf8n27BU2xnaPpK7vGA/gxSejZi/Sc9r47oxd1IyTfwM7MkuYdScdfzj+72UmatVzaYoCt/tzGb466v5ac8p9FrHKMzcW5M96oQGITyVQadl7q3dGdExCku1nfs/2c6XW7Ocdj2bXWHW93vPaI92V/8Ep11PiIYQ4mfkwzt6EmjSsyWjkEcXplLdQBt3iyus3PHRVnJKzLSOCOD18V3Raj1vF/sfSeEpuKZzND9PG8TIzlHY7AofrDvKVa/+xhdbMlXZ+b7vZAm3fbCFhz9P4XSFlQ7RQSydOoB7B7X0yNYRQngqo17LnIlJjO4WQ7Vd4fFFacxYuqfB3lhrlVRZuXv+Vj5anwHAjFEdpD2a8BhtIgN5+9ZkDDoNy3ad4q9f7bzs10hJlZW7P97KvpMlhAeY+OCOHgR6yaCLFJ4CcDSbf3tSd+bf1ZP4MD9ySsw8uXgXw15bzTcpx11SgGbkl/PE12lc+++1rEvPx6jTMm1oG5ZM7U+HmCCnX18IcTa9Tsvr47sxfVhbAOZvyOCmdzaSntswx/DWbhpcfTAPH4Oj0L1TRjqFhxnYpilvTfxfO7J7/7vtkmcOc0qqGPfORrYfO02Qj55PJvciPswzWyedi9MKz4yMDCZPnkxCQgK+vr60atWK5557DovFOYtvRcMY0i6Cn6cN4pnrOhDmbySjoIJHF+5kwEu/MufXQ+SVmhv0eoqisCPzNFM+28GVr/7Gwm1Z2BW4tks0v/x1MNOGtsWgk9+PhFCTRqPhL1e14Z1buxPooyc1q4iR/17LnF8PUWG5tDfXkiorLy7bx9i5GziaX050sA9fP9CP67pIpwrhmYZ3jOLtScn4GLSsOpDHmEvYN7EhPZ/r56xj/6lSwgNMLLjXcfqgN9E76xPv378fu93Ou+++S+vWrdm9ezf33nsv5eXlvPLKK866rGgAPgYdkwckMKFnHPM3ZPDR+gxySsy8svwgr688RJ+WoVzbOYYrEyOICq7/0XWKorD3ZAm/7Mvl25QTHMkvr/uzIe2aMvWK1rKLVQg3NKJTFF3jgnli0S7WHMzjleUHmb/hGA8MbsnN3eMI9vvzqcDT5RY+23yMD9dn1O0Cvr5rDLNu6EiIn9HZX4IQTjW8YxRf3NeXez7exsGcMm54ax2TB7TkgcEtL/jznVtSxesrD/HF1kwUBVo19eejO3vRPMzPheldw2mF54gRIxgxYkTdxy1btuTAgQPMnTtXCk8P4W/SM+WK1tw7sCXLdp1k/oYMUrOKWJ9ewPr0AgDiQn3p2SKUtpGBxIf6ERfqR5CPAT+TDptdocpqo6DcQnZRJUfyytl1opiUzCLyy/43cupj0HJt5xjuHZRAYpR3/WYnhLeJDvbl47t6siQ1m1dXHCCrsJLnf9jHyz8fYFCbpgxoHUbrpn4UVDkOfLDYLZwqrmLXiWLWpeezPj0fq81x0kvLpv78/dr2XJkYqfJXJUTD6RYXws/TBvLU4l0s35vDO6sP8+mmY4zoFMVViRG0iQzE36SjrKqavSdLWL43h5V7czDXtDS8pVccz1zXAT+j00o0Vbn0qyouLiY09PwjWWazGbP5fwVJSYljiNpqtWK1Wp2er/YarriWJ9EA13aK4NpOERwrrOCn3Tks35vD7uwSsgorySo8Ue/P6WfU0bdlKMPaR3B1x0gCTI4fxcv93ss99HxyDz3DtZ0iGJYYzuKUbD7bnMn+nDJW7sthZV3bJT2zUlaf8+92jAnkrn4tGNkpEoNOK/fazchr8PIFmbTMmdCFX/fn8fov6RzIKePr7cf5evv5TwNLigvm8avb0iO+CaBc1vff1fewPtfRKIrikkNGDx8+THJyMq+++ir33HPPOZ8zY8YMZs6cedbjCxYswM/P+4abPV1VNRwt05BRqiGvCvKrNJw2Q5UNLHYNGhQMWvDTQ6gJwkwKsf4KcQEKLQJATr8TwjsoCpyogH1FGo6UaMip1FBiBatdg0GrEGSAGD+F+ECFzk0UouSfc9GI2BU4Wgo78rUcK9OQUwk2xfEeGO0LLQIVeoTbaeYPntq4paKigokTJ1JcXExQ0IVnLutdeJ6vOPy9rVu30qNHj7qPs7OzGTx4MIMHD+b9998/798714hnXFwc+fn5f/qFNASr1cqKFSsYNmwYBoN3tC1Qi92uoNHg8vZHcg89n9xDz1d7D4cOHYrRKOs2PY28Bj2fq+9hSUkJ4eHhF1V41nuqferUqUyYMOGCz2nRokXd/2dnZ3PFFVfQt29f5s2bd8G/ZzKZMJlMZz1uMBhc+sPv6uuJhif30PPJPfR8RqNR7qEHk9eg53PVPazPNepdeIaHhxMeHn5Rzz1x4gRXXHEF3bt356OPPkKrlblVIYQQQojGymmbi7KzsxkyZAjNmzfnlVdeIS8vr+7PoqKinHVZIYQQQgjhppxWeC5fvpz09HTS09Np1qzZGX/mov1MQgghhBDCjTht7vvOO+9EUZRz/ieEEEIIIRofWXQphBBCCCFcwq3b4teOjtY2knc2q9VKRUUFJSUlspPPQ8k99HxyDz2f3EPPJvfP87n6HtbWaRczq+3WhWdpaSkAcXFxKicRQgghhBAXUlpaSnBw8AWf47KTiy6F3W4nOzubwMBAlzQir21Yn5WV5ZKG9aLhyT30fHIPPZ/cQ88m98/zufoeKopCaWkpMTExf9o6061HPLVa7Vk74l0hKChIXmweTu6h55N76PnkHno2uX+ez5X38M9GOmvJ5iIhhBBCCOESUngKIYQQQgiXkMLzd0wmE88999w5z4sXnkHuoeeTe+j55B56Nrl/ns+d76Fbby4SQgghhBDeQ0Y8hRBCCCGES0jhKYQQQgghXEIKTyGEEEII4RJSeAohhBBCCJeQwvN33n77bRISEvDx8aF79+6sXbtW7UjiIs2ePZuePXsSGBhIREQEo0eP5sCBA2rHEpdo9uzZaDQapk2bpnYUUQ8nTpzg1ltvJSwsDD8/P7p168b27dvVjiUuUnV1NX//+99JSEjA19eXli1bMmvWLOx2u9rRxHmsWbOGUaNGERMTg0aj4dtvvz3jzxVFYcaMGcTExODr68uQIUPYs2ePOmFrSOFZY+HChUybNo2nn36alJQUBg4cyDXXXENmZqba0cRFWL16NVOmTGHTpk2sWLGC6upqhg8fTnl5udrRRD1t3bqVefPm0aVLF7WjiHo4ffo0/fv3x2Aw8OOPP7J3715effVVQkJC1I4mLtJLL73EO++8w5w5c9i3bx8vv/wy//rXv/jPf/6jdjRxHuXl5XTt2pU5c+ac889ffvllXnvtNebMmcPWrVuJiopi2LBhlJaWujjp/0g7pRq9e/cmOTmZuXPn1j3Wvn17Ro8ezezZs1VMJi5FXl4eERERrF69mkGDBqkdR1yksrIykpOTefvtt3n++efp1q0bb7zxhtqxxEV48sknWb9+vcwUebDrrruOyMhIPvjgg7rHxo4di5+fH5988omKycTF0Gg0fPPNN4wePRpwjHbGxMQwbdo0nnjiCQDMZjORkZG89NJL3H///arklBFPwGKxsH37doYPH37G48OHD2fDhg0qpRKXo7i4GIDQ0FCVk4j6mDJlCtdeey1Dhw5VO4qop6VLl9KjRw9uvvlmIiIiSEpK4r333lM7lqiHAQMG8Msvv3Dw4EEAdu7cybp16xg5cqTKycSlOHr0KKdOnTqjtjGZTAwePFjV2kav2pXdSH5+PjabjcjIyDMej4yM5NSpUyqlEpdKURSmT5/OgAED6NSpk9pxxEX64osv2LFjB1u3blU7irgER44cYe7cuUyfPp2//e1vbNmyhb/85S+YTCZuv/12teOJi/DEE09QXFxMYmIiOp0Om83GCy+8wC233KJ2NHEJauuXc9U2x44dUyMSIIXnGTQazRkfK4py1mPC/U2dOpW0tDTWrVundhRxkbKysnjkkUdYvnw5Pj4+ascRl8But9OjRw9efPFFAJKSktizZw9z586VwtNDLFy4kE8//ZQFCxbQsWNHUlNTmTZtGjExMdxxxx1qxxOXyN1qGyk8gfDwcHQ63Vmjm7m5uWf9piDc28MPP8zSpUtZs2YNzZo1UzuOuEjbt28nNzeX7t271z1ms9lYs2YNc+bMwWw2o9PpVEwo/kx0dDQdOnQ447H27duzaNEilRKJ+nrsscd48sknmTBhAgCdO3fm2LFjzJ49WwpPDxQVFQU4Rj6jo6PrHle7tpE1noDRaKR79+6sWLHijMdXrFhBv379VEol6kNRFKZOncrixYv59ddfSUhIUDuSqIerrrqKXbt2kZqaWvdfjx49mDRpEqmpqVJ0eoD+/fuf1cLs4MGDxMfHq5RI1FdFRQVa7ZllgU6nk3ZKHiohIYGoqKgzahuLxcLq1atVrW1kxLPG9OnTue222+jRowd9+/Zl3rx5ZGZm8sADD6gdTVyEKVOmsGDBApYsWUJgYGDd6HVwcDC+vr4qpxN/JjAw8Kz1uP7+/oSFhck6XQ/x6KOP0q9fP1588UXGjRvHli1bmDdvHvPmzVM7mrhIo0aN4oUXXqB58+Z07NiRlJQUXnvtNe6++261o4nzKCsrIz09ve7jo0ePkpqaSmhoKM2bN2fatGm8+OKLtGnThjZt2vDiiy/i5+fHxIkT1QutiDpvvfWWEh8frxiNRiU5OVlZvXq12pHERQLO+d9HH32kdjRxiQYPHqw88sgjascQ9fDdd98pnTp1Ukwmk5KYmKjMmzdP7UiiHkpKSpRHHnlEad68ueLj46O0bNlSefrppxWz2ax2NHEeq1atOud73x133KEoiqLY7XblueeeU6KiohSTyaQMGjRI2bVrl6qZpY+nEEIIIYRwCVnjKYQQQgghXEIKTyGEEEII4RJSeAohhBBCCJeQwlMIIYQQQriEFJ5CCCGEEMIlpPAUQgghhBAuIYWnEEIIIYRwCSk8hRBCCCGES0jhKYQQLjRt2jRGjx6tdgwhhFCFFJ5CCOFCW7dupVevXmrHEEIIVciRmUII4QJWqxV/f3+sVmvdY7169WLz5s0qphJCCNfSqx1ACCEaA51Ox7p16+jduzepqalERkbi4+OjdiwhhHApKTyFEMIFtFot2dnZhIWF0bVrV7XjCCGEKmSNpxBCuEhKSooUnUKIRk0KTyGEcJHU1FQpPIUQjZoUnkII4SK7du2iS5cuascQQgjVSOEphBAuYrfbSUtLIzs7m+LiYrXjCCGEy0nhKYQQLvL888+zcOFCYmNjmTVrltpxhBDC5aSPpxBCCCGEcAkZ8RRCCCGEEC4hhacQQgghhHAJKTyFEEIIIYRLSOEphBBCCCFcQgpPIYQQQgjhElJ4CiGEEEIIl5DCUwghhBBCuIQUnkIIIYQQwiWk8BRCCCGEEC4hhacQQgghhHAJKTyFEEIIIYRLSOEphBBCCCFc4v8BActQYXTiiVkAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x250 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,2.5))\n",
    "plt.plot(t, uf(t))\n",
    "plt.title('Input $u(t)$')\n",
    "plt.xlabel('$t$')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3.54 s, sys: 37.1 ms, total: 3.58 s\n",
      "Wall time: 3.6 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "# Simulate using scipy.integrate.odeint method\n",
    "# Produces same results as Matlab\n",
    "fun = partial(pendulum_dydt, u=uf, **params)\n",
    "rtol = 10e-12\n",
    "atol = 10e-12 * np.ones_like(x0)\n",
    "X = odeint(fun, x0, t[:-1], tfirst=True, rtol=rtol, atol=atol)\n",
    "assert X.shape == (1000, 2), X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 3.142e+00,  0.000e+00],\n",
       "       [ 3.142e+00, -0.000e+00],\n",
       "       [ 3.142e+00, -0.000e+00],\n",
       "       [ 3.142e+00, -0.000e+00],\n",
       "       [ 3.142e+00, -1.000e-03]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X[:5].round(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAp4AAAD/CAYAAABYQdZFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABehUlEQVR4nO3dd3xT9frA8U+SpnvRlg5oC2W2bFr2RqYDZago4MCJC7levfpT73VcFbe4QHELoogiAiICsjcUCmVv2tKWLrpXmpzfH18ochm20OQ05Xm/XueV5DTJedqTNE++6zFomqYhhBBCCCGEnRn1DkAIIYQQQlwdJPEUQgghhBAOIYmnEEIIIYRwCEk8hRBCCCGEQ0jiKYQQQgghHEISTyGEEEII4RCSeAohhBBCCIdw0TuAS7HZbKSmpuLj44PBYNA7HCGEEEII8T80TaOgoIAGDRpgNF66TbNWJ56pqalEREToHYYQQgghhPgbycnJhIeHX/I+tTrx9PHxAdQv4uvra/fjWSwWlixZwuDBgzGbzXY/nqh5cg6dn5xD5ybnz/nJOXR+jj6H+fn5REREVOZtl1KrE88z3eu+vr4OSzw9PT3x9fWVN5uTknPo/OQcOjc5f85PzqHz0+scVmVYpEwuEkIIIYQQDiGJpxBCCCGEcAhJPIUQQgghhEPU6jGeQgghhBD2YrVasVgseodR4ywWCy4uLpSWlmK1Wq/4+cxmMyaTqQYik8TzHIVlFRRZIK/EgqsVDIDRYMBoMGAwgMGgbp/Zr/bJ+qJCCCGEM9E0jfT0dHJzc/UOxS40TSM0NJTk5OQay1P8/f0JDQ294ueTxPMvXlywl193uPDs1hXVepzRcG4iajIYcDMbcXMx4uZiUpdmI+4uptP7Tbibjfi4mfH1cMHX3YyPuwu+HmZ83M34ursQ6O1GfR83fN1dJLkVQgghatCZpDM4OBhPT8869zlrs9koLCzE29v7bxd0/zuaplFcXExGRgYAYWFhV/R8knj+hU3TLvNxf32suiyxXHnTNoCri5H6p5PQ+j5uBPu4EV7Pk8gAtUUEeODnYa5zbxohhBDCHqxWa2XSGRgYqHc4dmGz2SgvL8fd3f2KE08ADw8PADIyMggODr6ibndJPP/i7VFt6e+ZwtChQzG5uKBpoJ1OKm2ahgZoNiqvn9mPdjb51ACrVaPcaqXUYqOswkaZxaouK05fWmyUWKwUlFrIL60gv8RCQWkF+X+5nVVYRkFpBeUVNk7klnAit+Sicfu4u9A40IvmId60DPGhRagPLUJ8aODnLgmpEEII8RdnxnR6enrqHIlzOfP3slgsknjWFKPRgMkAZpMRs0vNDKK9EqUWK5kFZWQWlqnLgjJO5peScqqEpJxiknKKySxQCWriiTwST+Sd83gfNxdaNfClY2Q9OkT40zHSnxBfd51+GyGEEKL2kIaZ6qmpv5cknrWYu9lERIAnEQEX/1ZWUm4l5VQxhzOLOHCyoHI7kllEQVkFm47msOloTuX9G/i506lxAD2bBdKzWRDh9eQbnxBCCCEcQxJPJ+fhaqJ5iA/NQ3wY2ia0cn95hY3DmYXsTMklITmX7Um5HDhZQGpeKfN3pDJ/RyoAjQI96dksiEExIXRvGoi7Wf+WXiGEEELUTZJ41lGuLkZiwnyJCfNldOdIAIrKKtiRnMuGI9msO5TFjpQ8jmcXczw7iVmbkvByNdGvZTCDW4fQPzoYX3ep0SuEEEI4gxdffJGXXnoJgHfffZfx48dX+bFnutH9/PzsvsSUJJ5XES83F3o0C6JHsyD+ObglBaUWNh/NYeX+TJbsSedkfhm/JabxW2Iabi5GBrUKYVRsOL2bB+FikiJXQgghRG3WunVrli1bhre3NxUVFQCUlZXx5JNP8v3331NSUsKAAQOYOnUq4eHhlY9LS0tj9uzZvPDCC3aPURLPq5iPu5kBMSEMiAnhpRtbk3gijyV70vl9VzpHMotYuDONhTvTCPJ2ZVRsOOO6NbrkeFMhhBBC6MfFxYXQ0FBsNhv5+fkATJo0iQULFvDDDz8QGBjIP//5T2644Qbi4+MrZ6eHhobi5+fnmBgdchRR6xmNBtpH+NM+wp8nB7dk14l8ft6WwoIdqWQVlvPp6iNMX3OEAdHB3NG9Mb2bBWE0yoxAIYQQzk/TtBpbf7s6PMymKs8Wz8zMpG3btkycOJFnn30WgE2bNtG7d28WLlx4wcfk5eXxxRdfMGPGDAYOHAjAzJkziYiIYNmyZQwZMqRmfpFqkMRTnMdgMNA23I+24X48d30MK/dnMmPjcVYfyGTZ3gyW7c2gebA3j17TjOvbhkk3vBBCCKdWYrHS6j9/OPy4e14egqdr1VKx+vXr8+WXXzJ8+HAGDx5MdHQ048aN4+GHH2bw4MGsX7/+vMfEx8djsVgYPHhw5b4GDRrQpk0b1q9fL4mnqH3MJjXWc1CrEI5kFjJj43F+2prCwYxCHv8hgfeWHuChfk0Z0TEcVxdJQIUQQgh7ue6667j//vsZO3YsnTt3xt3dnddff/2i909PT8fV1ZV69eqdsz8kJIT09HR7h3tBkniKKmtS35sXhrVm0sAWzNhwjC/WHuVYdjFP/5zI1JWHeXpoNNe2CZVFeYUQQjgVD7OJPS87vvXP4zKWMHz77bdp06YNP/74I1u3bsXdvfqFYTRN0+2zWhJPUW1+HmYevaY543tG8f3mJD5ZdYTj2cU8/N02YiP9ee76GOIaBegdphBCCFElBoOhyl3eejty5AipqanYbDaOHz9Ou3btLnrf0NBQysvLOXXq1DmtnhkZGfTo0cMR4Z5H+kbFZfNyc+G+3k1Y9VQ/Jg5ojofZxLakXEZN28DE77eTUVCqd4hCCCFEnVFeXs7YsWMZPXo0r7zyCvfeey8nT5686P3j4uIwm80sXbq0cl9aWhq7du2SxFM4Ly83F54Y1IKVT/Xjts4RGA0wf0cqA95ZxXebjmOzaXqHKIQQQji95557jry8PD744AP+9a9/ERMTw7333nvR+/v5+XHvvffyz3/+kz///JPt27czbtw42rZtWznL3dEk8RQ1JsTXnddHtePXR3rRtqEfBaUVPPfLLm75dAPHsor0Dk8IIYRwWitXrmTKlCnMmDEDX19fjEYjM2bMYO3atUybNu2ij3vvvfcYPnw4t956Kz179sTT05MFCxZUruHpaJJ4ihrXNtyPeY/05IVhrfByNRF//BTXvr+GmRuPo2nS+imEEEJUV79+/bBYLPTq1atyX2RkJLm5uTz00EMXfZy7uzsffvgh2dnZFBcXs2DBAiIiIhwR8gVJ4inswmQ0ML5nFEue6Ev3JoGUWKw8P28Xd3+1hcyCMr3DE0IIIeqcxMREvL29L9kCeiHe3t5MmDDBTlGdyzmmcAmn1dDfg+/u68pX64/xxuJ9rDqQyfUfrOGjMbF0iZKZ70IIIURNmDhxIuPGjQMgMDCwWo9NSEgAcEj3uySewu6MRgP39oqid/MgHvluGwczCrn9s408PbQl9/duIut+CiGEEFcoICCAgADVoPPXWu1V0axZM3uFdR7pahcO0yLEh3mP9OSmDg2w2jReW7SPh7/bRkm54+vjCiGEEMLxJPEUDuXl5sKU0R14ZXgbXE1Gft+VzujpG8jIlzU/hRBCiLpOEk/hcAaDgXHdGvHd/V2p52lmZ0oeN328jj2pVe8WEEIIIYTzkcRT6KZz4wDmPdKTpvW9SMsr5eZP1rPuUJbeYQkhhBDCTiTxFLpqFOjF3Id60qNpIMXlVsZ/tYXFu9L1DksIIYQQdiCJp9Cdn6eZr8Z3ZmjrUMqtNh7+Lp6f4lP0DksIIYQQNUwST1EruLmY+GhMR26JC8emwZNzdjBjwzG9wxJCCCGcwosvvojBYMBgMPD+++9X67FnHufv72+f4P7Cronn5MmT6dy5Mz4+PgQHBzN8+HD2799vz0MKJ+ZiMvLmze24r1cUAP/+dTezNiXpHJUQQgjhHFq3bk1aWhr3339/5b7p06fTr18/fH19MRgM5Obmnve4tLQ0pkyZ4pAY7Zp4rlq1ikceeYSNGzeydOlSKioqGDx4MEVFRfY8rHBiBoOB566P4YE+TQB49pdEftyarHNUQgghRO3n4uJCaGgonp6elfuKi4sZOnQozz777EUfFxoaip+fnyNCtG/losWLF59z+6uvviI4OJj4+Hj69Oljz0MLJ2YwGPi/a6Mpr7Dx9fpjPP3zTswmAyM6husdmhBCiLpI08BS7Pjjmj2hitX7MjMzadu2LRMnTqxMIjdt2kTv3r1ZuHDhRR83adIkAFauXHml0dYIh5bMzMvLA6gs6fS/ysrKKCsrq7x9ptyTxWLBYrHYPb4zx3DEscTfe3Zoc8osFXy/JYUn5+zE29VIvxb1L/kYOYfOT86hc5Pz5/zq+jm0WCxomobNZsNms6md5UUYX3d844btmRRw9arSfQMDA/n8888ZOXIkAwcOJDo6mnHjxvHQQw8xcOBA1q1bp57TZkPTNIDK3/PM/jOXlb/3X2P5n/td6OeapmGxWM6r6V6d14rDEk9N03jiiSfo1asXbdq0ueB9Jk+ezEsvvXTe/iVLlpzTbGxvS5cuddixxKV1McGh+ka2ZBp55LttPNrKSiOfv3+cnEN9uVtOEZKXQFDBXnzKUnG3nMJoq6DC5E6JOYBczyiyfGLI8G2H1eh2weeQc+jc5Pw5v7p6Ds90RxcWFlJeXq52Worx1yGW/IICMFe9bHSvXr248847GTt2LB07dsRsNvN///d/5OfnU1ZWhtVqPadGe0FBQeX14uLiyn1G4/kjLUtLS9E07aI13svLyykpKWH16tVUVFSc87Mzz10VDks8H330UXbu3MnatWsvep//+7//44knnqi8nZ+fT0REBIMHD8bX19fuMVosFpYuXcqgQYMwm812P56omiFWGw/O3M6aQ9l8dcSD2fd3ISrowt8Q5Rzqy5CyGePGjzEc+B2Ddv63ZrOtBA/LKQKKD9Mkaxmaqze21qOw9ZgI/o0AOYfOTs6f86vr57C0tJTk5GS8vb1xd3dXOzUf1froYL7V6Go/4/3336ddu3bMmzePzZs3ExwcDICbmxsmkwlfX180TaOgoAAfHx8Mp5//TAOej4/PBXMqd3d3DAbDRfOt0tJSPDw86NOnz9m/22kXS1YvxCGJ52OPPcb8+fNZvXo14eEXb8p2c3PDze381g+z2ezQF7+jjycuzWyGT+7oxO2fbWRnSh73fLuNuQ/3INjH/RKPkXPoULlJ8MezsHfB2X0N46D5EAhrD/4RYHKDsjw4dQySN8P+RRhykzBt/wZTwkzoNB6ueR7M3oCcQ2cn58/51dVzaLVaMRgMGI3Gc1v+TFXoTqsFjh07RmpqKjabjeTkZDp06ABQmWAajcbK7vIzv+eZ/WcuL9Ti+b/3u9DPDQbDBV8X1Xmd2DXx1DSNxx57jF9++YWVK1cSFRVlz8OJOszLzYUv7+7MqGnrOZ5dzIMz4vn+/m64m01//2BhP5oGCbPg96ehvAAMRugwFro/CsHRF35MwzhoMwqGTIbja2HtFDj8J2z5HHbPw3D9FEf+BkII4TTKy8sZO3Yso0ePJjo6mnvvvZfExERCQkL0Dq3K7Lqc0iOPPMLMmTOZNWsWPj4+pKenk56eTklJiT0PK+qoIG83vrq7M77uLmxPyuXZXxIrB1ALHVgtsHAS/PqwSjojusKEdXDTRxdPOv/KaISoPnDHXLhrIQS1hOIsXOaMo03KTKgo+/vnEEKIq8hzzz1HXl4eH3zwAf/617+IiYnh3nvvveRj0tPTSUhI4NChQwAkJiaSkJBATk6OI0I+j10Tz2nTppGXl0e/fv0ICwur3GbPnm3Pw4o6rEl9bz4eG4vJaGDuthNMX31E75CuTsU5MGMExH8NGOCaf8P43yGk1eU9X1RvmLAWuj0CQNPMJZi+G6mOI4QQgpUrVzJlyhRmzJiBr68vRqORGTNmsHbtWqZNm3bRx33yySd07NixclH5Pn360LFjR+bPn++o0M9h18RT07QLbnfffbc9DyvquN7N6/OfG1SC8/riffy596TOEV1linPg2xvh2Bpw9Ybbf4A+T4LxCoc9uLjC0NeouPU7yk2eGFM2wecDIUe+XAghRL9+/bBYLPTq1atyX2RkJLm5uTz00EMXfdyLL75Yq3IxqdUunNKd3Rtxe5dINA0mzU4gKVuHhX+vRmeSzvRE8KoP9/wBLYfW6CG05kNY2/zfaH4RkHMYvhgCGftq9BhCCFEXJSYm4u3tfckW0Avx9vZmwoQJdorqXJJ4CqdkMBh46cbWxEb6U1BawUPfxVNqqfpaaOIylBfBzJFnk867FkLohdfkvVIFHg2puHsxhLSFogz4Zpgkn0IIcQkTJ05k//79JCQkMGbMmGo9NiEhgV27drF9+3Y7RXeWJJ7Cabm6GPloTCwBXq7sTs3npQV79A6p7rJZ4ef7IHU7eASopLMqE4iuhHcI3DUfQs8knzdA9mH7HlMIIZxUQEAAzZo1o1mzZtWuu37mcY5YfUgST+HUGvh7MGV0BwwG+H5zEr9sT9U7pLrpj+dg/yK1FuftP9g/6TzDMwDuPJN8ZqoW18IMxxxbCFGnyaoo1VNTfy9JPIXT69OiPo8PaA7AfxbsIbVI54Dqmh0/wKbT44VGfgqRXR17fM8AGDcX6jVWi89/dwuUFTo2BiFEnXFmsfPqlHkUZ/9eV1pUwGElM4Wwp8euaU788VOsOZjF1wdNjC231smKGw53cjcsmKSu930GWo/QJw7vYJV8fjEI0hJgzl1w+2wwyb8wIUT1mEwm/P39ychQvSeenp6VVX/qCpvNRnl5OaWlpRetRFRVmqZRXFxMRkYG/v7+mExXtoKJ/NcWdYLJaGDK6A4MnbKak4XlvLnkAK+MaKd3WM6trAB+vBMqSqDpNdD3X/rGE9gUxsxRYz0PLYNlL8CQV/WNSQjhlEJDQwEqk8+6RtM0SkpK8PDwqLGk2t/fv/LvdiUk8RR1RqC3G2+MbMM9325j5qZk+keHMCDGecqI1Tp/PAvZh8C3IYz87MrX6awJ4XEwfJpq8dzwkRr72f42vaMSQjgZg8FAWFgYwcHBWCwWvcOpcRaLhdWrV9OnT58a6f0zm81X3NJ5hiSeok7p3TyIvmE2VqUZ+ddPO/l9Um+Cfdz1Dsv57F8M274FDDByOngF6R3RWa2HQ/qTsOZtmD8RglpAw1i9oxJCOCGTyVRjCVVtYjKZqKiowN3dvdYNO5PJRaLOGRZpIzrEm+yicp6as1NmLlZXUTbMf0xd7/4INO516fvrof9z0GIoWMtg9h1SWlMIIZyEJJ6izjEb4d1b2uHmYmTVgUy+WX9M75Cch6bBb/9Q62bWj1Y12Gsjo1G1xAY0hfwUmPeQil0IIUStJomnqJOah3jz7HUxgKrnfixL1liqkn0LYc+vYHSBEZ+CuRYPU3D3g1u+VmuLHlisxnwKIYSo1STxFHXWHd0a0aNpIKUWG0/9tAObTVrELqmsABadnrne83Fo0EHXcKokrB1c+7q6vuxFSN6sazhCCCEuTRJPUWcZjQbeGNUOL1cTW46d4ivpcr+0Fa9BQapaqL3PU3pHU3Vx46HNKLBVwJzxUHJK74iEEEJchCSeok6LCPDk/053ub/1xz6OSpf7haXtgE2fqOvXvQNmD33jqQ6DAW6Ycna858InZLynEELUUpJ4ijpvbNdIejY73eU+ZwdW6XI/l80GC/8Bmk1VJmo+UO+Iqs/dF0Z9psam7p4LiXP0jkgIIcQFSOIp6jyD4WyX+9bjp/hq3VG9Q6pdEn+EE/Hg6g1DJusdzeVrGKfKegL89k/ITdI3HiGEEOeRxFNcFcLrefLs9We63PfLLPczyovUpByA3v8E3zBdw7livf4B4V2gLB9+eQhsVr0jEkII8ReSeIqrxpguqsu9rMLGs78kysLyAOs+gII08I+Ebg/rHc2VM7nAyE9V6+3xtbD+Q70jEkII8ReSeIqrhsFg4LURbXE3G1l/OJs58Sl6h6SvvBOw7n11fdB/a/eandUR0ASGnl5iafkrkLZT33iEEEJUksRTXFUaBXoxaWALAF79bS+ZBWU6R6SjP1+GihKI7AGtbtI7mprVcRxE3wA2C8x9ACqu4vMshBC1iCSe4qpzX68oWoX5kldi4eWFe/QORx8nd8PO2er6kFfVkkR1icEAwz4Ar/qQuRdWvaF3REIIIZDEU1yFXExG3hjVDqMBFuxIZcW+DL1DcrzlrwKaaulsGKt3NPbhFQjXv6uur52iZu4LIYTQlYveAQihh7bhftzTM4rP1x7l+Xm7WPKPPni5XSVvh5StsP83MBih//N6R2NfrW6ENjfDrp9g3sPwwKq6M5ZV6K80D7IOQV4SFGVBcbaqnGW1gGZVqyq4uIGbr1pr1t0f/CPAL1JdOlOhBiFqyFXySSvE+Z4Y3ILFu9NJOVXCO0sO8J9hrfQOyTH+fEldth8D9VvoG4sjXPcWHF0Nmftg1esw8MUL3k3TNNLzS0nOKSGzoIzMglKKyq1YrDYqrBquLka83FzwcXMhxM+dhv7uNPD3wNNV/o1eFYpz1Je25E2QskW9ngpPXtlz+kdCaDsIbQthHSCyG3j410S0QtRa8h9TXLU8XV14dURb7vpyM1+vP8pNHRrQPsJf77Ds68hKlYSZXKHf03pH4xieAXDDezB7rJrFH30DhHciu7CMLcdy2HQ0h4TkXA6dLKSgrKLaT98o0JPWDXxpFeZLbKN6xDWqh5uLyQ6/iHAoTVOlZA/8AQd+h9TtF76fdwjUiwKvILV5BKhWToNRbRVlal3Z0nzVIpqXrIoblBeqy9wk2Lfw9JMZVBLauBdE9YUmfaVVVNQ5kniKq1rfFvUZ3qEB8xJSefrnnSx4rBdmUx0d+qxpaiY7QKd7VGvL1SLmBmh7KyT+yKlZ93K/5xS2phSfdzeT0UB4PQ+CfdwI8nbDx90FF5MRs9FAWYWNonIr+SUWTuaXciK3hILSCo5nF3M8u5hFiekAuJuNdG4cQO/mQVzbJoyIAE9H/7biSmQfhh3fw47Zqgv9rwKaqlbJ8M4Q1g4Cm4G7X/WPoWmqBTVjN6Qnqi15M+QchvSdats4Fcye0PQaaHkdtLxWfYkSwslJ4imuev++oRWrDmSyL72AL9YeZULfpnqHZB8Hl6gJNmZPVaXoKlFcXsHCnWn8ln4zb2lLCC4+xsD8z9nKGFqG+NAlKoBOjesRHepLVJAXri5V/+KRU1TO3rR8dqfmkXgin41HssksKGPNwSzWHMzitUX7aNvQj2vbhjKiY0PC/KT1qlaqKIc9v8KWz1RX+hlmT2jSH1oOheaDwSe0Zo5nMKjJb1F91HZGfhocXwfH1sLBpZCfolpD9y0Eo4uKod1oaDFUxioLpyWJp7jqBXq78dz1rXhyzg6mLDvAtW1CaRTopXdYNUvTzi4p1Pk+8A7WNx4HyCos46t1R/l2w3EKSlUX+vOme5lufocHXBZx67gJBET3vqJjBHi50rNZED2bBQFqnOihjELWHspi2d6TbDicTeKJPBJP5PH2H/u5JjqY27tE0q9lMCZjHVvCyhkV50D8V7D5M1XBC1T3eNMB0P421dLo6sAWa98waHuz2jRNtXzuWwT7foOTibB/kdrc/aD1CIi9q+6uSiHqLEk8hQBGxTZk7rYU1h/O5vl5u/j2ni4Y6tLaloeXq9ZOFw/o8Zje0djVyfxSpq08zA9bkii12AA1DnN05whujhsAS09g3PkDAUsnQdO1NTqGzmAw0DzEh+YhPozvGUV2YRl/7D7JvIQTbD6aw7K9GSzbm0GYnzt39WjMmK6R+Lqba+z4oopKclVX9oapUF6g9nmHQOf7IfaOmmvZvBIGA4S1V1v//4OMfbDzB9j5I+SfgPiv1dYwTsXdeoS0ggqnIImnEJwtpzlkymrWHMxiXsIJRnQM1zusmqFpsPotdT3u7jrb2llcXsGnq44wffURSixWANpH+PNwv6YMignBeKaF8drX1SSr7EOqpOaQV+0WU6C3G2O6RjKmaySHMwv5YXMSP8WnkJZXyuu/7+Oj5YcY2zWS8T2jCPWTpMHuygpg4yew4UO1FBJASBv1Zaz1SHBx1Te+SwmOVisyXPMfOL4Wts2APfPUF8oT8bDkOeh4B3S5H/zqyP8uUSfV0VkUQlRf4yAvJg5oDsB/F+4lp6hc54hqyLG1kLRBzWTvOVHvaGqcpmn8FJ9Cv7dW8v6fBymxWIlrVI/v7uvKvId7MKR16NmkE8CjHgw7XaN+w8eQtNEhcTat781z17di47MDeOvmdjQP9qawrIJPVx+h95vLeX5eIifzSx0Sy1XHZoPt38GHcbDiFZV01o+BW7+FB9eobvXanHT+ldGoxoWO+gz+sQeu+Tf4hqsZ8+umwPvt4ZcJkLFX70iFuCBJPIX4iwf6NKFliA85ReW8+lsd+ce9+k11GXsn+DbQN5YadiyriLGfb+LJOTvIKCgjMsCTqWNj+WlCd3o2C7r4cImWQ6HDWECDeQ9BeZHDYnZzMXFLpwj+mNSHL+/uRJfGAVisGjM3JtHnzRW8+tueuvOlpzZI2QpfDIRfH1brbgY0gVFfwEPrVeUuoxN/DHrXhz5PwuM74LZZ0Lg32CrUrPyp3WDWbQ77YiVEVTnxO06Immc2GZk8qi0GA/y8LYV1h7L0DunKJG1S63YazdBzkt7R1JgKq41PVh1myJTVrD+cjbvZyNNDo1n6RB+uaxtWtfG5Q14DnwaQc+TsMlMOZDQauCY6hB8ndOeHB7rRqVE9yipsfLbmKL3fWM57Sw9QXF79dUXFaSW5MP8x+HyA6op29YFB/4WHN6nJO86ccP4vkwtEXw93L4T7lkPMMMCg1h/9cgh8OVSN89Y0vSMVQsZ4CvG/YiPrcWe3Rnyz4TjP/ZLI4kl9cDc76YLga0/XKu9wuyrRVwck5xTzj9kJbD1+CoCezQJ5bUTb6q9E4OEPN30IM0fBpk/UwvJRVzbL/XJ1axLInAndWXkgk7f/2M/u1Hze//MgP2xJ4umh0Qzv0PDc4QLi0vYuhN/+CYVqbVU6jIMB/wGfEH3jcoTwOBg9E7IOwvoPYMcPaqjNjBEQ0U1NVIrqqyYv/Y3i8gpy88vJKiynqKyC4nIrJZbTl+VWNE09jcFgwAAYDeDhasLLzQUvNxe83VzwcnUhwMuVIG9XXOrqGsmiWiTxFOICnhzSkj92n+RYdjEfLj/IU0Oi9Q6p+jL2wYHFgAF6PK53NDXi14QTPP/LLgrKKvB2c+E/w1pxS1z45a9A0GygWpJm2zeqK/ahDeDmXbNBV5HBYKB/y2D6Nq/P77vSmfz7XlJOlfDEjzv4Zv0x/jOsFXGNZAHxSyrMhEVPqkk3AIHN4cYPoVF3XcPSRdDp373fs6pi19YvIXkjfHsTRPaAfs+QE9yNYznFJGUXcyy7iKTsYpJyiskoKCU910T5huU1Fo7BAIFertT3cSfYx42G9TxoHOhJo0AvGgd6ERngiYerk37BF9UiiacQF+Djbualm1rz4Ix4Pl11hGHtGxAd6qt3WNWz4UN1GX09BDXTN5YrVGqx8u95u5gTnwJAbKQ/U0Z3JDKwBtZYHPIqHF6hShcu/bcqr6kjo9HA9e3CGBATzFfrjvHR8oPsSMlj1LQNDGvfgP+7NpoG/rIQ/XkOLFFfHooywWCCno9D36dliSHfMGxDJnM8+l60NVOIPPojLknr4dsbOWiL5l3LLWzSYi7wQPVlzt1sJNBLVfHydDXh6eqCh6sJD7MJowE0VA++BthsGqUWK4VlFRSWVVB0+vJUsQWrTSOrULWe7k27cKgN/T2ICfMhJsyX6FBfYsJ8aBToJWve1jGSeApxEUNahzKkdQh/7D7JMz8n8vNDPZznH2B+mir5B+oD2Ikl5xQzYWY8u1PzMRrg0WuaM/GaZjXXbefmAzd9BN/eqFqFYoapMoU6czebeKhfU0bFNeTdJQeYvTWZBTtSWbbnJI8PbM49PaOqVWWpzrKUwNIXYPOn6nZwKxjxiVr/8ipVXF5BQlIuW4+fYuvxU2xPOnW6iMJ1hNCNh11+5TbTCroa9zHb7b8kuMWxKvwhXMI7EhngSZCXC7u2buDmGwbj7+V+xWsaW20aOUXlZBSUklFQRkZ+KSmnSjiWXczx7CKOZRWRX1rBidwSTuSWsGxvRuVjPcwm2ob7ERtZj9hIf2Ib1SPI2+0K/0JCT5J4CnEJL93YhnWHsklIzuW7Tce5s3tjvUOqmk2fgM2ixnRFdNE7msu26kAmj/+wndxiCwFernx4e8fKKkE1qklftQj3ls/g18fg4fWXV4PbDoJ93Hl9VDvu6N6IF+fvZsuxU7z++z5+ik/hvze1oXvTQL1D1M/J3fDzfZCxR93u+pBa6/Iqa+WssNpISM49Xao1kx0peVht504k8jCbTrcmRuLSoA8HfIuIPjgdc8IMOpTF0+HwfeA+Eto9j8U3kozd4O3mUiOFNExGA/V93Kjv40bri9znVFE5B04WsDctn71pBexNz2d/egElFiubj+aw+WhO5X0jAzyJa1SP7k0C6dEskPB6DqwuJa6YJJ5CXEKonztPD23Jv3/dzZuL9zO4VWjtX+i7rAC2fqWuO+m6nZoGn64+yjvLDqJp0D7cj6nj4mhozy7mQS/BoWVw6ij88Szc9LH9jnUZWjfw48cHu/PzthNMXrSXQxmF3P7ZRoZ3aMCz18cQ7FPLX5c1bdsMNYHIWgZewTB8KjQfpHdUDpOeV8qyvSdZfSCTDYezKSg7dwWEMD934hrVo1OjenRqHEB0qM/5vQQxU6DXRFgxGRLnwO65sOdXjB3G4W5xbCnOel6udG0SSNcmZ79IWW0aRzIL2Z6Uy7akU2xLOsWBk4Uk5aixqL9sPwGoRLRH00C6n96uuveCk5HEU4i/MbZrI37ZfoJtSbm8MH8Xn97RSe+QLi3+GyjLUxMrWlyrdzTVVl5h4/vDRjZlHgTg9i4RvDCstf1XFnD1UsnLV9fB9pkQcxO0GGzfY1aTwWDg5rhwBsWE8NaSfXy3KYl5Can8uTeDJ4e0ZFy3Rs4zHORyWUrh96dg27fqdrNBMHyaWtOyjjuUUcgfu9NZsuckO5Jzz/mZv6eZns2C6NM8iJ7NgqreChjQRC1G33Mi/PlfOPgHpu3fMNAwC3wOQ58nVNEFHZiMZ0vQ3tpZrcqRV2IhITmXLUdzWH84ix0peZWJ6A9bkgFo3cCX/i2D6R9dnw4R9er+e8LJSOIpxN8wGg1MHtmO6z9Ywx+7T7J4VzpD29SCWs4XYrWoGtSgygA62VqFecUWHpwRz6ZMI0YDvDCsNXf1aOy4ABr1gG4Pw8aPYcFEeHiDbh+6l+LnaeaV4W25JS6Cf/+6i50pebwwfzc/bk3mleFt6BhZ+2KuEaeOwY93QtoOwADXPAe9/ul0r/Pq2J9ewK8JJ/hjdzqHM88WOjAYoGOEPwNiQujVLIg2Df2uLMEKbQtjf4Tj67EtfRFTyibY8IFa8aHX49B1gvpypjM/DzN9W9Snb4v6QEsKyyoqk9D1h7PZk5bP7lS1fbTiEP6e6v79WwbTt0V96nk5SYWqOkwSTyGqoGWoDw/2bcLHKw7zwvxd9GwWiI+7We+wzrdrLuSfUF2P7UbrHU21HM8uYvzXWziSWYSbUeOjsbEMaq1DpaUB/4aDSyD7IPz+DIz81PExVFH7CH9+ebgnszYn8dbifexOzWfE1PXc1jmCfw2NJqAufcgeXKrGc5bmgkcA3PxFrZgEZg+5xeXM35HKnK0pJJ7Iq9xvNhno0TSIIa1DGdgq2D5dyo16YL1zIZt/mEzXwsUYMvaoAgubPoU+T6nlx2pReVFvNxf6RwfTPzoYgKzCMlbtz2TF/gxWH8gkt9jCrwmp/JqQitGg1mlWE0dDa2ZVDFFtDkk8p06dyltvvUVaWhqtW7dmypQp9O6tz0LNQlyux65pzm870ziWXczk3/fx2oi2eod0Lk0729rZ9QGnmmCRkJzLPV9vIaeonDA/d+5oVEi/Fjp1nZo9VNftl4Nh5w9qlnvMDfrEUgUmo4E7ujXi2jahTF60j5+3pfDDlmQW707nX0Oiua1zhHMvPq9pah3KZS8CGjSIVTXW60hBhDMqrDbWHMzip/gUlu45SbnVBqhk85roYK5rG0b/6GB8HfGF12DgpF8HKkY/jXn/fFj+CuQeV2ukbvgI+j8HbWpn9acgbzdGxYUzKi6cCquNbUm5rNifwYp9GexLL6ic6f/qor1Eh/owuHUog1uF0LqBb41MpBJ/z+6J5+zZs5k0aRJTp06lZ8+efPrpp1x77bXs2bOHyMhIex9eiBrjbjYxeWQ7bv9sI7M2JTG0dSh99EqOLiR5M6QlgMkN4sbrHU2VrT6QyYSZ8RSXW2nb0I9PxnZg65o/9Q0qorNahmrte6rsYngn8KmlwytOC/J2451b23Nblwj+PW8X+9ILePaXRGZvSeK/w9vQLtxf7xCrr6IMFkyCHbPU7bjxcO0b4FJ3ltM5lFHIT/EpzN2WQkZBWeX+VmG+3NIpnJs6NNSv5dpogna3Qqvhqst91ZtquMPc+2HtFNU70GJolaog6cHFZKRLVABdogJ4emg0J3JLWHp6jOymoznsSy9gX3oBH/x5kIb+HgxuHcKQ1qF0alRPqizZkUHT7Fu8tWvXrsTGxjJt2rTKfTExMQwfPpzJkydf8rH5+fn4+fmRl5eHr6/9F++u2PUrB9fNo0V0a0xmN1Xf2mgCk6ta68/dF9z81DIr7n7gGahq5Ipaw2KxsGjRIq677jrMZvu0DLw4fzdfrz9GmJ87iyf1wc+jlnS5zxmvZqV2GAfDa9eM7IuZvyOVf/6YgMWq0bt5EJ+Mi8PVqNn9HFZJRRl8PhDSd6ou3bE/18oWngupsNr4dsNx3l16gMKyCgwGGNMlkqeGtMTf075JTI29BwszYfY4VW3HYFIJZ5f7ay5QHeWXWliwI5Wf4lPYnpRbuT/Ay5WbOjTg5rhwWjfQbzmvi57D8iK1VNva99UERoCIrqocaeNe+gR7mU4VlbN8XwZ/7E5n9cFMSi22yp/V8zQzMCaEwa1D6d08yClLJjvis/CvqpOv2TVrKi8vJz4+nmeeeeac/YMHD2b9+vXn3b+srIyysrPf+PLz8wH1B7RYLPYMFQDDnnnEpM2FtLlVur9mMIJnEPiEonmHgHcomm8DtHpRENAELaBprVkL8Gpx5nViz9fLPwY0YcW+DI7nFPPS/F28MbKN3Y5VZfmpuOz5FQNgibsXHPB+uVIzNibx30X70DS4vk0ob45qg6tRc8g5rBoj3PQJLl8MwHB4OdYNH2PrMkHnmKrujq7hDGlVnzf/OMCvO9L4blMSixLTeGpwc0Z1tF/t9xo5fxl7cPlxLIa8ZDQ3X6wjv0Rr0s8pXtcXY7VpbDiSw8/bTrB0bwZlFSrRMRkN9G0exKjYBvRrUb+yKICer/+LnkODK3SbCO3vwLjhQ4xbPsOQvAm+vh5bkwFY+z8Hoe10iLj6vF0N3NguhBvbhVBSbmXtoWyW7stgxb5MThVbmBOfwpz4FDzMRno1C2JQTDD9WgZRz85f3GpKRcp2IrNWYrE4Zomx6rxe7drimZqaSsOGDVm3bh09evSo3P/aa6/xzTffsH///nPu/+KLL/LSSy+d9zyzZs3C09P+g4Ajs1dRr+gIBq0Co2bFoFkxYsNos+BiK8FsLcFsLcbFWozZWoKBv//Tlbn4UOQWQr57Q/I9GpHnEUG+RyQVJil558yO5MMHu01oGLivpZW2AXbtOPhbMalzaHFyAVneLVnX/DldY/k7mga/pxj5I0V9wPYOsTEyykZtHYbYOGs57ZO/xmpwYXXLF8n3cL4hQofyYM5RE+kl6o/cyFvj5igrkfqUpb+kkLwEOh2bioutlEK3EDY1+QeF7jpMMqshmSWwOdPI5kwDueVnX+ShHhpdg210CtLwdY5c5jzullO0SP+VRlmrMGIF4IR/F/aGjaLIPUzn6C6PVYMj+QZ25hhIzDFw6i/nzIhGU1+NtgFqC6iFIz7MFQXEpP1M46wVaAYjy6MnU+Ru/2FCxcXFjBkzpkotng5JPNevX0/37t0r97/66qvMmDGDffv2nXP/C7V4RkREkJWV5ZCudovFwtKlSxk0aNDfN03brFCcBQXpGApPQmE6hoJ0DHkpcOoIhpwjGIoyLvpwrV4UWlh7tPAu2MK7QEgbMEq3/ZWq1jm8Qq8v3s8X644T5O3Kosd66PdN2FKCy4ftMZTkUDHqK7ToYfrEUQVWm8bLv+1l1mZVc33iNU15tF+Tcwb1O/IcVommYZozDuPBP9DqR1MxfqmagORkLFYbMzYm8cHywxSVqyThxnZhPDGoWY0uzH8l58+47WuMi/+FQbNha9wb68gva+VyVn+nsKyC33edZO72E2w9nlu539fdhWHtwhgV24A2tXgyS7XP4amjmFa9jmH3XAxoaAYTWvvbsfb+F/g675cGTdPYk1bAsr0ZLNubwb6Thef8PCbUh0ExwQyMCSY61Fvf82mrwLjtW4yrXsNQmgtAin83Am+fikuA/b8s5+fnExQUpH9Xe1BQECaTifT09HP2Z2RkEBISct793dzccHM7/yuE2Wx26AdQ1Y5nBrdwqBd+8buUFUDOEcg+BBl7IT1RbfknMJw6iuHUUdgzDxOA2QvC46BRTzWerGGcGl8qLosjXjNPDY1h1cFsDmUU8p/5+5g2Llaffzy7ZkNJDvhF4NLqxlo77riswso/f9rBb4lpGAzw8k1tuKNbo4ve39Hv+0saPhWmdseQuQ/zylfgujf1jqjazGZ4sF9zhsdG8Mbv+5i7/QTzd6axeM9J7ukZxcP9m9bojOlqnT9Ng+X/hTXvqNsd78B4w3sYTbXk/FeBzaax6WgOc+KT+T0xnRKLSu6NBujdvD63dApnYEyIU40XrPI5DG4Bt3wJvZ+A5f/FcGAxhoSZGBPnQOd7ocdE8HXOFtAOjQLp0CiQJ4fGkJRdzJI9anLS1mM57E0vYG96AR+sOEx4PQ8Gtwqlf3R9OjcOcNx51jS1/NuSf0PW6V7kkDZUDH6N+F25XBcQ6ZD/o9U5hl0/oVxdXYmLi2Pp0qWMGDGicv/SpUu56aab7Hno2sHNB8Laq+2virLVhIWUrZC8Sc1GLsuDo6vVtnKyGhvapB80HQDNBoJfQ11+BXFx7mYT793agZHT1rF4dzrfb05mTFcHd8NqmhrsD9D5vlqbdBaWVfDgjK2sO5SN2WRgyuiOXN/OiT6IvILUEkvfjYLNn0KzAdBiiN5RXZYQX3feHd2Be3pF8cpve9h4JIdPVh3mx63JTLymGbd3jcTNxYHJUUW5Wjlg5w/qdr9noe+/au1M6f+VnFPMz9tS+HlbCsk5JZX7mwR5cXOncEZ2DK/9ZXZrSmgbGDMbkjbCspcgab1a4m3LFxB3F/Sc5NSfZZGBntzXuwn39W5CdmEZf+7LYMnuk6w5mEnKqRK+XHeUL9cdxd1spHuTQPq2qE+fFvWJCvKyT6NEeiL88RwcXaVuewRA/2chbjyaTYNdi2r+mDXA7p9STzzxBHfccQedOnWie/fuTJ8+naSkJCZMcJ5B+jXOKxCa9lcbgM0GmfvU7M0jK9VWmgd7flUbQMNO0OpGiLkRAqL0ilz8j7bhfvxrSDSvLtrLywt307lxPZqH+DgugOPr1T8fFw+IvdNxx62G7MIyxn+9hZ0peXi5mvj0jk70ah6kd1jV13wgdH0INk2DXybAhDXgd4kej1quTUM/vr+/G8v3ZfDaor0czizixQV7+GzNUR7p34yb48IrJ7rYTWk+/HiH+p9nMMGNH0DHcfY9Zg0oKqtg8a50fopPYcOR7Mr93m4uDGsfxs1xEcRG+tfarnS7i+wG4xfB4eWw6g3VwLJ5OsR/DR3vgF7/cPp1WAO93bi1UwS3doqguLyC1QeyWLb3JKsPZJJRUMaK/Zms2J8JQESAB31b1KdXsyA6Nw4g0PsKB4fmp6m1VRO+AzS18k7XB6H3k+Dhr+5jq70T8eyeeI4ePZrs7Gxefvll0tLSaNOmDYsWLaJRo4t3sV11jEYIaaW2TveAtQJSt6k37aFlqmX0xOlt6X9UabO2t6rKND7nD1kQjnVvryhWH8xkzcEsHvt+O/Me6em4bpYzrZ3tR4NngGOOWQ0pp4q584vNHMkqIsDLla/Hd3bO9STPGPQSJG1Q66XOGa8+XJ2oO/h/GQwGBsSE0LdFfX7YksyHyw9yIreEZ39J5OMVh3jsmmaMigvHbI81DfNT4btb4OQucPWGW79RvTu1lNWmseFwNnO3pbB4dzrFp8fJGgzQo2kgt8RFMKR1KB6uztOVblcGg+oZaHqNapFb+YZqAd36BWz7FjqOVQlovcZ6R3rFPF1dGNomlKFtQtE0jX3pBaw+kMmqA5lsOZZDck4JMzcmMXNjEgDNgr3pEhVA19NrjIb5VXGMdckpWP8hbJwGlmK1r/VIGPiCU/0d7b6O55Vw9Dqejl73qsoK0mHfQtgzH46tBU39w8NgUv+oO4yBltfWqUWVL5de5zCjoJTr3l9DVmE5d3ZvxMs3OWCJpfxUeK+Nej08tEF9calFDpws4I4vNnEyv4yG/h58e28Xmtb/+2nUtfZ9eEbOUfi0rxoe0/1RGPKq3hHVmFKLle83JzF15WEyTy9m3sDPnbt7Nua2LpFVGgNapfOXsRdm3gz5KeAdAmN+hAYdavA3qTkHTxbw87YTzNt+gvT80sr9jQM9GRkbzsjYhoTXq1ulF+32Hjy2Fla+DsfWqNsGE7Qert5HDWNr7ji1SFFZBRuPZLPqQCabTy9a/78a+nvQtqEfbcP9aNPQj7YN/c4tGlCar5LNDR+fXT81vAsMeU0Vu7iAq3YdT1FDfELV+L3O96nxoXvnw47vVffFwT/U5hGgxtB0vs+pu/+cVbCPO2/d0p7xX23h2w3HiY2sx/COdh7LtP07lXRG9qh1SWf88VPc8/UW8kosNA/2Zsa9XevOOLeAKDXZaPZYVT6wUQ+Ivl7vqGqEu9nE+J5R3N4lkpkbj/PJqiOk5pXy2qJ9vL/sIKM7RzK+Z2MiAq4g0Tq2Fn4Yo4YTBTaHcT9DvdrVA5acU8yixDQW7kw7p1a6r7sLw9o3YGRs+NXdlX65GveCuxeqIUKr3oQjK2DXz2pr3Bt6PAbNBjlNoYaq8HJzYUBMCANiVO/kqaJyth4/xeaj2Ww+msOu1HxO5JZwIreExbvPTsRu6O9Bu2AXbrH9Ts+T3+FmUa9DLbgVhv7Pqf85Tvr6k8TT2XgFQqfxass6pErJ7fgB8k+o8n7rPlB1pbtOgMjuTvvCdEb9WwbzaP9mfLTiEM/M3UmLEB9aNbBTS73NqrqrAOLuts8xLtOK/Rk8NDOeUouN2Eh/vry7s92r5ThczA2qlWbDR/DLQ/Dgqjo19trdbOK+3k0Y160Rvyac4PM1RzmYUciX647y9fqj9GsZzOjOEVwTHVy9bvhdP6vxsdZyiOgGt39fa4aIJGUX81tiGosSz002XYwG+rUMZlRsQ66JCXbsxKu6qlEPuHMepO1U76FdP6tW0GNrIKgldH8E2t4CrnWrJRmgnpcrg1qFMKiVSkQLyyrYdSKPxJQ8Ek/ksetEHrlZadxaOIe7SpbgbygC4JCtAVMqRrEivQchizwJW7eJUF8Pwvzcqefliq+7C74eZnzcXfBxM2My2MgqVVXMalvHkSSeziyomSpV1v852P+7Gu93bM3ZSUkNO0Gfp9TsW0lAHeIfg1qw80ReZf3xBY/2ws/TDu/6wysgLwnc/dWks1pi3vYTPDlnBxU2jX4t6zN1bCyernX038zAF9WKFCmbYc7dcM9ip1zf81LczSZGd47k1k4RrDqQyedrjrL2UBbL92WwfF8GQd5ujIpryLB2DWj9d+tSbvgY/nhWXY+5EUZO1/XvZbVpJCSfYsW+TJbvy2BPWn7lz4wG6BoVyHVtQ7mubdiVTwYRFxbWTr0OBvxHfX5t/VotCbRgIiz9N7QfoxpZ6rfUO1K78XZzoVuTQLo1CYRTx2D9d2jbZ2KoUCskZLuF85P3GL4r6kJyXjlaucaRzCKOZBZV4dlduHaQBQ/32vX6raOfCFcZo0m1wMTcACd3w6ZPYedsNRnp+9GqhFmfpyD6hjrVhVEbmYwGPritA8M+WktSTjGPz97OF3d1xlTTZXniv1KX7W+rNcnOl2uP8vLCPQDc1KEBb9/S3j6TUmoLkxlu+Qo+6a0mGy14HEZ8Wie/5BkMqtWvX8tgjmQWMntrMj/Hp5BVWManq47w6aojNAr05No2YQyOCcL215kDNptKIjZ8pG53naDGpumwTvGJ3BI2Hclm5f5MVh/MJLf47MxfowG6Nw3kurZhDG4VSn2f2vVhXaf5hcPgV9Tn1LZvYfNnkHtcrSCxaZrqho+7W32GmevIkJ0zbDY4/Cds+RwO/AFoGAAadISekwiMGcaDRhMPosZgp+eVkpZXSnp+ibrMK+VUsYWCUgv5JRbySysoLK2grMJKcWk5bvZemeIyyOSiv6j1kxqqozBDzX7b8gVYTn8zCmsPg/4LTfrqG5sd1ZZzuDs1j5FT11NWYeO+XlE8f0MNjsEsSIf3WoOtolZMKrLaNF79bS9frjsKwN09GvOfG1pddi3w2nIOq+zIKpgxQo23HfRf6DlR74gcwmK18efeDH5NOMGK/RmUWmyVP/Ny0egXHUa/Zn5cd+S/eB2Yp34w6GW1mLgDkvMKq40jWUUkJOWy8fR4upRTJefcx9fdhb4tg+nfsj59W9SXls3TdH8P2mxqVZetX8KB30E7/dpy84VWN6kVXRr1dO6GlMIMSJilfsfc42f3N71GrXca1eeK3icyuUg4nncwDP6vegFvmgYbP4G0HfDtjWrw9qCXIKS13lHWWa0b+PHWLe2Z+P12Pl97lIgAT+7q0bhmnjzhO5V0RnTVPeksLq/g8R8SWLrnJABPDWnJw/2aXl2TLpr0haGvw+9PwbIXIDgGmg/SOyq7M5uMlUvIFJVVsHJ/JosS01h5IIOiMiurdx1h3P538TLtwYILXwU9RXHJtbTem0FUkBeRAZ41sk6ozaZxsqCUo1lFHMsqZm9aPrtS89ibln9OMgyqR6JNQz96Ng3kmuhgOkT441KXW+WdldGo1s1tPhDyUiD+GzWhNi8Zts9Qm29DaHuzGrbRINY5ktCyAti7EBJ/VGvXnkmo3f2gw1i1nGJQc11DdARJPOs6r0C45nnVvbXqTbWG2qGlan3QTuPV2BonrIXsDG5s34DknGLe+mM/Ly3YTQN/j8oB5ZfNZlP/hEH3SUWZBWXc980WdqTk4Woy8vat7bmxvfPWZb4iXe5X61Fu+wZ+uhfu//Oq+AA5w8vNhevbhXF9uzCKS8uY+f0sRmR+SP3iwxRqHjxomcS6lLaQcrDyMUYDNKznQbi/J4HergR6uRLg5YaPuwtmkwGzyYjZZMSqaZRZrJRabJRYrOQWW8gqLCOzoIzMwjJSThWfl2Ce4elqok0DPzpH1aNrVCBxjerh5SYfe07FLxyueQ76/Z9aQ3fnbNg9T02oXfe+2ryCoeVQaHkdRPWtXZOSCk6qlWf2L1atuBV/aXVv2El9DrceWbtitjN5B14tvIJUfemuD8KfL8OeeaqJf+8CNbam3eg6OTZNbw/3a0rKqWK+35zMY99v45vxXejaJPDyn/DoStUt4+YHrYbXVJjVtj+9gHu/2ULKqRL8Pc18dmcnOjeuHbOTdWEwwHVvQ9YB9eE4azTct6zWzNh2JHPOQcalvoynJQe8QzCP/pF/ao0YlJzLzpQ89p8s4GhWEcXlVpJzSs4pM3m5TEYDEfU8aBzkRYsQH1o38KVNQz+iAr0ue8iHqGWMRmjcU23Xvqnqk++eC4f+hKIMNTZ027dgNEN4JzUutHEvdd3Vy3FxluZB0ia1WP6RVaoYzF8FNlMFYNreDIFNHRdXLSKJ59UmsKmqEHJsLSx8Qs0g/OVB2DYDhr2vZsqLGmMwGHj5pjak55WyYn8m93y9hRn3dSU28jJbmeO/VpftR+v2DXlRYhpPztlBcbmVxoGefDW+C1FBDvzHXlu5uMKtM+CzayDnsEo+75pfayZ/OcTxDbh8PxqzJQ8tsBmGcXNxq9eIWDjnNa9pGpkFZRzNKiI9v5TswnKyi8rIKSqnsMxKhdWGxWqjrMKGyWjAw2zC3WzC3WzE18NMfW836vu4EeTtRgN/D8LredTtiWziXGZ3tZpHqxuhohyOr1Uru+z/XXXHJ21Q2+o3wWCEoBYQ1kHNcwiOhoAm4BsOpitIgTQNijIhY4+a1HtyjxrOdnIX8D9TZxrEqiIvLYaoyb5XeSOPJJ5Xq8a9YMJaNdN01ZvqjftJLzX4v/N9zjFexkmYTUamjYvjnq+3sP5wNnd9uZnv7+9Gm4Z+1XuiwgzY95u6HntXzQf6N6w2jbeX7GfaysOAKhP40ZjYcytsXO2868PYOfDlYLXM0k/3wugZuszgdrhdc+GXCRisZeR4NcPnzkWY/S48tMRgMBDs606wbx2boSwcz8VVTchpeo1qCc05ohpWzmwFqZC5T207fzj7OKMZ/CPBJ0z1THgGqs3FXb1fjS7q0lIC5YVQXqQqCBWkQt4JVTmu4iKt9QFN1FqlkT1U2VCfUMf8LZyEJJ5XMxdX6P0EtBkF8x9T9XR/fwr2/wY3fSwVkGqQu9nE53d14q4vN7Pl2CnGfLaRL+/uTKfqdE8nzFKTihp2glAHlOT8i1NF5Tw+O4HVBzIBuL93FE8PjZaJGRcSHA23/wDfDlfvpUVPwfXv1N1WDk2D1W/BClU61NbiWtZ5jGLoVTjMQOjMYFC9eoFNVSU/UKuApCao1si0HZB9SK2XaS1TPRM5h6/kgKpwREhrCG6tLiO6SKL5NyTxFKpc3R3z1DpiS/+jZttN7Q43fgCtR+gdXZ3h6erCl3d35u6vthB//BRjP9/E1LGxlaXULslmUxNXwOGTitYfzuIfsxM4mV+Gu9nImzdfxZOIqqpRDxj1Gfx4l5rQ5xmoJkjUNRVl6kvrztnqdrdHsPb/D7bFf+gblxBn+ISenng09Ow+m021XOYcVeNDi3OgOFttFWWqMpytQm1mD3D1BjdvNVbUJ0zNqPdrqC5dZAmu6pLEUyhGI3R9QHVXzJsAKVtUNZbj69XkI3lz1QgfdzMz7+3KI7O2sXxfBg/MiOfV4W24rUvkpR94bI3qQnL1gTYjHRKrxWrjvaUHmLbqMJoGTep78dHtsfYrA1rXtLpJdf39/pQaa2Zyhb5P6R1VzSnKgh/GQvJGMJjg+rfVcjAWy98/Vgg9GY2qR0969XQh/WTiXEHNYPxi6PWEur15Onw5RHVNiBrh4Wri0zviGNmxIVabxjNzE3n2l0TKKqwXf9CZSUXtbnXIDM1dJ/IYMXUdU1eqpPP2LhEsfKyXJJ3V1fUBNW4aYMUrsPY9feOpKRl74fMBKul084NxP6mkUwgh/oYknuJ8JhcY+AKMmaPW+EzdDtP7wdE1ekdWZ5hNRt6+pT1PDGqBwQCzNiUx+tONnMi9wGD1oiy17BWcHbdkJyXlVl5btJebPl7HrhP5+HmYmTo2lskj29Xdmuv21vNxuObf6vqyF9W6g85s11z4bID6MurfCO5bqnpKhBCiCiTxFBfXYjA8uEYtBVFyCmYMh61f6R1VnWE0Gpg4oDlf3t0ZPw8zCcm5DHlvNbO3JHFOJdsd34PNomr3hrW3SyyapvHbzjQGT1nF9NVHsNo0rm8XxtIn+nBd2zC7HPOq0udJ6PuMur70P7DsJTUpx5lYK+CP5+Cn8aoMb1RfuH8F1G+pd2RCCCciiae4NP8IGL9IzXy3VcDCSbDoX2rwtagR/VsGs+DRXsRG+lNYVsHTPycy7otN7EvPV8nJmW52O0wq0jSNjUeyGTF1PY/M2kZyTglhfu58cVcnPh4TS7CPLHdTY/o9c7blc+27sOBx53kfFZxUXzw3fKRu95wE4+aqymhCCFEN0ncm/p7ZA0Z9AfVj1Di1zZ+qGYEjP1cL+YorFhnoyZwJPfhi7RHeXnKAdYeyue79NTwdncWD2YfUrMo2o2rseDabxtK9J/l01WG2JeUCqrzgA32acH/vJlJW0B4MBtXy6RkIvz2hVikoPAkjPwP3Wjx29sAfMO9hKM5Sr8PhU9XEKSGEuAzy6SKqxmBQM3KDmsHcB9SYw+9uhtu+A/dqLoQuLshkNPBAn6YMbR3G5N/38vuudEIO/QAmWOvRD5LK6N7UG9MVlAA8mlXEL9tSmLv9BCmn1HhSVxcjt3YKZ+KA5tLC6QidxqsFq3++Hw4shs8Hwu3f177yeZZSWPYCbPpE3Q5pAzd/KV3rQogrIomnqJ7WI1SLzfdj1BI/X10P434GnyqsRSmqJDLQk2nj4ti27zBtf9gCwOsZ3dj1xSZCfd25JiaYXs2CiI2sR4ivG4ZLLEyeWVDGrtQ8Nh7JZtX+TPalF1T+zMfdhTu6NeLuno0l4XS0Vjepkn2zx6qytZ9dAyM+UWX1aoOkTWp9zqz96na3h2HAC9LDIYS4YpJ4iuqL6gPjf4OZo+BkInx9Hdy1EHxlEkpNij21GLBQEtiGDhH9SN6RRnp+KbM2JTFrUxIAvu4uNAr0wt/TjI+7CzYblFVYySwsIy23lOyi8nOe02iA3s3rMzK2IYNbheLhehWUcqytwuPggZUwe5xaN/f726DTvWrdXFdPfWIqzYc/X1bFJNDAK1hVMWsxWJ94hBB1jiSe4vKEtYd7/oBvb1IlyL6+Hu5eCL5S0aZGaBrEq0pFHt3u4ZXObfn3Da1YezCLNQez2HA4m0OZheSXVpB4Iu+iT2MwQNP63rQP96dPiyB6N68vtdVrE59QuPs3lext+EhVOTq6GoZNgca9HBeHzaaqD/35shq/DdBxHAz6rxoWIIQQNUQST3H5ApuqZPPrYare7dfXq5ZPv4Z6R+b8kjaqbk6zJ7S9BQA3FxMDYkIqS2yWVVg5nFFEWl4Jp4otFJZaMJmMuJoM1PdxI8TXncaBXjJRqLZzcYMhr0KzgTDvIcg+qN5LbW5WrZ/27EnQNDiyApa+AOk71b56jWHY+9Ckn/2OK4S4asknkrgy9RqfTj5vUCUdv7lBtYR6B+sdmXM7s4RSm1EXnfHs5mKiVQNfqSZUVzTtDw+th+WvwNYvYddPsH+RqgjUY2LNjqO2WdUEwXVTVIEIADdf6P1P6DpBxnIKIexG1vEUV65eIzXm0z9SJZ8zRkJJrt5ROa/iHNj9i7oeN17fWIRjeQbADe+qsZ/hXcBSrLrgp7RVs+CPrVXd4pcr66DqTp/SDubcpZJOFw+VbE7cDr0mSdIphLArafEUNcM/Eu6YB18OVROOvr9NLTCt1yQJZ7bzR7CWQUhbaBirdzRCDw06wL1L4NCfsOoNSNkMiT+qzTsEmg9Wk/zCOkBgMzBeoA1B01RZy/RESNoAh5dD5r6zP/eoB53vh64PgleQg34xIcTVThJPUXMCm8Idc9USS0kbVIvKbbPAZNY7MudxTqWiu9TsIHF1Mhig+UBoNgBSt6nJZrvmqkXnt89QG4DRrJJRr0AwmECzQUkOFGZARen/PKdJjSXtcDu0uFZaN4UQDieJp6hZoW1hzGyYMQIOLoGF/4AbP5QEqqpStkDmXtX9eXpSkbjKGQzQME5t170Fx9fBwaWQslW1ZlaUQH6K2v6XyRWCY6BBR2jSX7WSyix1IYSOJPEUNa9Rd7jla/jhdtUqE9hMjR0Tf69yUtFI8PDXMxJRG7m4QdNr1AZgrYDCdFVLvThLtZgbDOARAN71wbeh9DgIIWoVSTyFfbQcCkMmw+KnVdm9gCip7/x3SnJVVypA3N16RiKchckF/MLVJoQQTkBmtQv76TYBujygrs99AFLi9Y2ntkuco7pNg1tBeGe9oxFCCCFqnCSewr6GTFYzcCtKVV3qgpN6R1Q7nTOp6G4ZEyuEEKJOksRT2JfJBW7+EupHQ0GamuleUf73j7vanNgGJ3eBizu0u1XvaIQQQgi7kMRT2J+bD4z+TlVGSdoAS57TO6LaJ/4rddlquFpfUQghhKiDJPEUjhHUDEZ8qq5vng4J3+sbT21Smg+7flbXZVKREEKIOkwST+E40ddB36fV9YWTIH2XruHUGolzVGnEoJYQ2U3vaIQQQgi7kcRTOFbfZ85ONvppPJQX6R2RvjTtbDe7TCoSQghRx0niKRzLaIThn4BPGGQdgEVP6R2RvlK3q+ozJjdof5ve0QghhBB2JYmncDyvQBj1ORiMkPAd7Jitd0T62faNumx1o5QyFEIIUedJ4in00bjXX8Z7/gOyDukbjx5K82HnHHVdJhUJIYS4CkjiKfTT5ylo3BssRfDT3VBRpndEjpU4R/3uQS2gUU+9oxFCCCHsThJPoR+jCUZ+Bp6Bapzjysl6R+Q4mgZbz0wqGi+TioQQQlwVJPEU+vINg2Hvq+vr3oekjfrG4ygn4uFkoqpUJJOKhBBCXCXslngeO3aMe++9l6ioKDw8PGjatCkvvPAC5eVSLlH8j5hh0H4MaDb45UEoK9Q7Ivvb+qW6bD1CJhUJIYS4atgt8dy3bx82m41PP/2U3bt389577/HJJ5/w7LPP2uuQwpld+zr4RcCpY3W/pGbJqbOVijrdo28sQgghhAO52OuJhw4dytChQytvN2nShP379zNt2jTefvvtCz6mrKyMsrKzE0zy8/MBsFgsWCwWe4Va6cwxHHEs8T9Mnhhu+ACX70ZA/NdUNBuC1mxQtZ/GGc6hcfssTBWlaMGtqQjpALU4Vj04wzkUFyfnz/nJOXR+jj6H1TmOQdM0zY6xnOP5559n8eLFbN269YI/f/HFF3nppZfO2z9r1iw8PT3tHZ6oBdqkfEfTzD8odfFjecxrWFx89A6pZmka1+z7P3xKU9kRfifH6g/UOyIhhBDiihQXFzNmzBjy8vLw9fW95H0dlngePnyY2NhY3nnnHe67774L3udCLZ4RERFkZWX97S9SEywWC0uXLmXQoEGYzWa7H09cgKUEly8HYMg6gK3VCKwjPqvew2v5OTQkrcdlxo1oZi8qHt8FbnUssa4Btf0cikuT8+f85Bw6P0efw/z8fIKCgqqUeFa7q/1irZJ/tWXLFjp16lR5OzU1laFDh3LLLbdcNOkEcHNzw83N7bz9ZrPZoS9+Rx9P/IXZDCM+hc8HYtzzC8Z2t0D09ZfxNLX0HCbMAMDQdhRmb5lUdCm19hyKKpHz5/zkHDo/R53D6hyj2onno48+ym23XXr5l8aNG1deT01NpX///nTv3p3p06dX93DiatQwFno8BuumwMIn1OLqHv56R3XlirJhz6/qukwqEkIIcRWqduIZFBREUFBQle574sQJ+vfvT1xcHF999RVGoywbKqqo3zOwbyFkH1Kz3G/6WO+Irty2b8BaDmEdoEFHvaMRQgghHM5umWBqair9+vUjIiKCt99+m8zMTNLT00lPT7fXIUVdYvY4nWwaYPtMOPSn3hFdGasFtnyurnd9UN9YhBBCCJ3YLfFcsmQJhw4dYvny5YSHhxMWFla5CVElkd2gywPq+oJJUFagazhXZO8CyD8BXvWhzSi9oxFCCCF0YbfE8+6770bTtAtuQlTZgP+AfyTkJcGyS09qq9U2faouO90DLudPoBNCCCGuBjLoUtRubt4w7AN1fctncHy9vvFcjtTtkLwRjGaZVCSEEOKqJomnqP2a9oeOd6jrvz4KlhJ946mujZ+oy9YjwCdU31iEEEIIHUniKZzD4FfAJwxyDsOK1/SOpuoKTp6ty95tgr6xCCGEEDqTxFM4Bw9/uP5ddX3DR3AiXtdwqiz+K7BZILwLNIzTOxohhBBCV5J4CucRfR20vQU0G8x7BCrK/v4xeiovOjupSFo7hRBCCEk8hZMZ+gZ4BkHmXljzjt7RXNr2mVCSA/UaQ8xNekcjhBBC6E4ST+FcvALhurfU9TXvQPoufeO5GKsF1n+orveYCKZqFwkTQggh6hxJPIXzaT0Com8AWwX8+ghYK/SO6HyJP0FeMngFQ4exekcjhBBC1AqSeArnYzDA9e+Aux+kJcCGD/WO6Fw2G6x9T13v/jCY3fWNRwghhKglJPEUzsknFIZMVtdXTIasg/rG81f7F0HWfnDzlQXjhRBCiL+QxFM4rw5joOkAsJapLnebVe+IVGvnytMJcef7VKusEEIIIQBJPIUzMxhg2Pvg6g3Jm2DzZ3pHBHvmwcldqrWzx2N6RyOEEELUKpJ4CufmHwGDXlLX/3wJsg/pF4u14mxVpe6PgmeAfrEIIYQQtZAknsL5xd0DUX3AUozp1wkYNJ1mue+cDdkHwSMAuj2kTwxCCCFELSaJp3B+RiMM/wTc/TGmJRCd9ovjYygvghWvquu9JoG7r+NjEEIIIWo5STxF3eDXEIZNAaD5yYUYkjY49vhrp0D+CfCLhC4POPbYQgghhJOQxFPUHa1HYGt3OwY0TPMfhpJcxxz31HFY/4G6PuQVMHs45rhCCCGEk5HEU9Qp1sGvUeQajCEvGeY9rJY3srclz0NFKTTuDTE32v94QgghhJOSxFPULW4+bI16GM3kCvt/g7Xv2vd4+36DvfPBYIRr31BLPAkhhBDigiTxFHVOrmcTrEPeUDeWvwKHltnnQMU5sPAf6nqPiRDS2j7HEUIIIeoISTxFnaR1vANi7wQ0+OleyKrh9T01DRY9BYUnIagF9Pu/mn1+IYQQog6SxFPUXde+BQ3joDQXZo6AgvSae+6tX8Kun8Bggpumgtm95p5bCCGEqKMk8RR1l9kdbp8N9aIgNwlm3gyleVf+vClbYfEz6vrAFyGi85U/pxBCCHEVkMRT1G3e9eGOueBVH04mwrfD1djMy5V1CGbdCtZyaHmd1GMXQgghqkEST1H3BTSBcXNVKcvUbfD1DZCbXP3nyTmiuuyLsyGsA4ycLrPYhRBCiGqQxFNcHcLawd2/gXcIZOyG6f3g6OqqP/5EPHwxWHXZBzSBsXPAzcdu4QohhBB1kSSe4uoR0gruWwahbaE4C74ZBgsmQVH2xR9TUQYr31BJZ1Gmeuz4xeAd7LCwhRBCiLrCRe8AhHAo/0i4Z4maHLTtG4j/Cnb8AO1ugeZDVGLp4qa64o+sgPhvID9FPTbmRrjpY3D31fd3EEIIIZyUJJ7i6uPqCTd+AO1uVeUuU7fDtm/VdiE+YTDov9D2ZhnTKYQQQlwBSTzF1atxL7h/BRxbC3vmwdE1cOoo2CrULPgGsdDqRmg9UtbpFEIIIWqAJJ7i6mYwQFRvtYGqSHRmvxBCCCFqlCSeQvyVJJxCCCGE3cisdiGEEEII4RCSeAohhBBCCIeQxFMIIYQQQjiEJJ5CCCGEEMIhJPEUQgghhBAOUatntWunl7bJz893yPEsFgvFxcXk5+djNpsdckxRs+QcOj85h85Nzp/zk3Po/Bx9Ds/kaWfytkup1YlnQUEBABERETpHIoQQQgghLqWgoAA/P79L3segVSU91YnNZiM1NRUfHx8MDlhfMT8/n4iICJKTk/H1lXrczkjOofOTc+jc5Pw5PzmHzs/R51DTNAoKCmjQoAFG46VHcdbqFk+j0Uh4eLjDj+vr6ytvNicn59D5yTl0bnL+nJ+cQ+fnyHP4dy2dZ8jkIiGEEEII4RCSeAohhBBCCIeQxPMv3NzceOGFF3Bzc9M7FHGZ5Bw6PzmHzk3On/OTc+j8avM5rNWTi4QQQgghRN0hLZ5CCCGEEMIhJPEUQgghhBAOIYmnEEIIIYRwCEk8hRBCCCGEQ0jiedrUqVOJiorC3d2duLg41qxZo3dIooomT55M586d8fHxITg4mOHDh7N//369wxJXYPLkyRgMBiZNmqR3KKIaTpw4wbhx4wgMDMTT05MOHToQHx+vd1iiiioqKnj++eeJiorCw8ODJk2a8PLLL2Oz2fQOTVzE6tWrGTZsGA0aNMBgMDBv3rxzfq5pGi+++CINGjTAw8ODfv36sXv3bn2CPU0ST2D27NlMmjSJ5557ju3bt9O7d2+uvfZakpKS9A5NVMGqVat45JFH2LhxI0uXLqWiooLBgwdTVFSkd2jiMmzZsoXp06fTrl07vUMR1XDq1Cl69uyJ2Wzm999/Z8+ePbzzzjv4+/vrHZqoojfeeINPPvmEjz76iL179/Lmm2/y1ltv8eGHH+odmriIoqIi2rdvz0cffXTBn7/55pu8++67fPTRR2zZsoXQ0FAGDRpEQUGBgyM9S5ZTArp27UpsbCzTpk2r3BcTE8Pw4cOZPHmyjpGJy5GZmUlwcDCrVq2iT58+eocjqqGwsJDY2FimTp3KK6+8QocOHZgyZYreYYkqeOaZZ1i3bp30FjmxG264gZCQEL744ovKfaNGjcLT05MZM2boGJmoCoPBwC+//MLw4cMB1drZoEEDJk2axNNPPw1AWVkZISEhvPHGGzz44IO6xHnVt3iWl5cTHx/P4MGDz9k/ePBg1q9fr1NU4krk5eUBEBAQoHMkoroeeeQRrr/+egYOHKh3KKKa5s+fT6dOnbjlllsIDg6mY8eOfPbZZ3qHJaqhV69e/Pnnnxw4cACAHTt2sHbtWq677jqdIxOX4+jRo6Snp5+T37i5udG3b19d8xsX3Y5cS2RlZWG1WgkJCTlnf0hICOnp6TpFJS6Xpmk88cQT9OrVizZt2ugdjqiGH374gW3btrFlyxa9QxGX4ciRI0ybNo0nnniCZ599ls2bNzNx4kTc3Ny488479Q5PVMHTTz9NXl4e0dHRmEwmrFYrr776KrfffrveoYnLcCaHuVB+c/z4cT1CAiTxrGQwGM65rWnaeftE7ffoo4+yc+dO1q5dq3coohqSk5N5/PHHWbJkCe7u7nqHIy6DzWajU6dOvPbaawB07NiR3bt3M23aNEk8ncTs2bOZOXMms2bNonXr1iQkJDBp0iQaNGjAXXfdpXd44jLVtvzmqk88g4KCMJlM57VuZmRknPctQdRujz32GPPnz2f16tWEh4frHY6ohvj4eDIyMoiLi6vcZ7VaWb16NR999BFlZWWYTCYdIxR/JywsjFatWp2zLyYmhp9//lmniER1PfXUUzzzzDPcdtttALRt25bjx48zefJkSTydUGhoKKBaPsPCwir3653fXPVjPF1dXYmLi2Pp0qXn7F+6dCk9evTQKSpRHZqm8eijjzJ37lyWL19OVFSU3iGJahowYACJiYkkJCRUbp06dWLs2LEkJCRI0ukEevbsed4yZgcOHKBRo0Y6RSSqq7i4GKPx3LTAZDLJckpOKioqitDQ0HPym/LyclatWqVrfnPVt3gCPPHEE9xxxx106tSJ7t27M336dJKSkpgwYYLeoYkqeOSRR5g1axa//vorPj4+la3Xfn5+eHh46BydqAofH5/zxuR6eXkRGBgoY3WdxD/+8Q969OjBa6+9xq233srmzZuZPn0606dP1zs0UUXDhg3j1VdfJTIyktatW7N9+3beffdd7rnnHr1DExdRWFjIoUOHKm8fPXqUhIQEAgICiIyMZNKkSbz22ms0b96c5s2b89prr+Hp6cmYMWP0C1oTmqZp2scff6w1atRIc3V11WJjY7VVq1bpHZKoIuCC21dffaV3aOIK9O3bV3v88cf1DkNUw4IFC7Q2bdpobm5uWnR0tDZ9+nS9QxLVkJ+frz3++ONaZGSk5u7urjVp0kR77rnntLKyMr1DExexYsWKC37+3XXXXZqmaZrNZtNeeOEFLTQ0VHNzc9P69OmjJSYm6hqzrOMphBBCCCEc4qof4ymEEEIIIRxDEk8hhBBCCOEQkngKIYQQQgiHkMRTCCGEEEI4hCSeQgghhBDCISTxFEIIIYQQDiGJpxBCCCGEcAhJPIUQQgghhENI4imEEEIIIRxCEk8hhHCgfv36MWnSJL3DEEIIXUjiKYQQQgghHEJqtQshhIPcfffdfPPNN+fsO3r0KI0bN9YnICGEcDBJPIUQwkHy8vK49tpradOmDS+//DIA9evXx2Qy6RyZEEI4hoveAQghxNXCz88PV1dXPD09CQ0N1TscIYRwOBnjKYQQQgghHEISTyGEEEII4RCSeAohhAO5urpitVr1DkMIIXQhiacQQjhQ48aN2bRpE8eOHSMrKwubzaZ3SEII4TCSeAohhAM9+eSTmEwmWrVqRf369UlKStI7JCGEcBhZTkkIIYQQQjiEtHgKIYQQQgiHkMRTCCGEEEI4hCSeQgghhBDCISTxFEIIIYQQDiGJpxBCCCGEcAhJPIUQQgghhENI4imEEEIIIRxCEk8hhBBCCOEQkngKIYQQQgiHkMRTCCGEEEI4hCSeQgghhBDCIf4f5WvE0PqOulQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x250 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,2.5))\n",
    "for i in range(X.shape[1]):\n",
    "    plt.plot(t[1:], X[:,i], label=f'x[{i}]')\n",
    "plt.xlabel('t')\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate derivatives (vectorized version)\n",
    "dXdt = pendulum_dydt(t[:-1], X.T, uf, **params).T\n",
    "assert dXdt.shape == (1000, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Identify system dynamics using sparse regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A. Using SparseNonLinearModel class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SparseNonLinearModel(['theta', 'theta_dot', 'u'], ['theta_dot', 'theta_ddot'], estimator=LinearRegression(), custom_features=['sin(x0)', 'cos(x0)', 'sin(x0)**2', 'cos(x0)**2', 'sin(x0)*cos(x0)'], poly_order=2)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Labels for inputs and outputs\n",
    "x_names = ['theta', 'theta_dot', 'u']\n",
    "y_names = ['theta_dot', 'theta_ddot']\n",
    "\n",
    "# Choose additional features for model identification\n",
    "custom_features = ['sin(x0)', 'cos(x0)', 'sin(x0)**2', \n",
    "                   'cos(x0)**2', 'sin(x0)*cos(x0)']\n",
    "\n",
    "poly_order = 2\n",
    "model = SparseNonLinearModel(x_names, y_names, \n",
    "                             custom_features=custom_features,\n",
    "                             poly_order=poly_order)\n",
    "model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>theta_dot</th>\n",
       "      <th>theta_ddot</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>-6.123234e-16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-1.348867e-11</td>\n",
       "      <td>-1.000483e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-9.955803e-05</td>\n",
       "      <td>-1.992160e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-2.978226e-04</td>\n",
       "      <td>-2.975509e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-5.939895e-04</td>\n",
       "      <td>-3.951002e-02</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      theta_dot    theta_ddot\n",
       "0  0.000000e+00 -6.123234e-16\n",
       "1 -1.348867e-11 -1.000483e-02\n",
       "2 -9.955803e-05 -1.992160e-02\n",
       "3 -2.978226e-04 -2.975509e-02\n",
       "4 -5.939895e-04 -3.951002e-02"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dXdt = pd.DataFrame(dXdt, columns=y_names)\n",
    "dXdt.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>theta</th>\n",
       "      <th>theta_dot</th>\n",
       "      <th>u</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3.141593</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>-0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3.141593</td>\n",
       "      <td>-1.348867e-11</td>\n",
       "      <td>-0.020010</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3.141592</td>\n",
       "      <td>-9.955803e-05</td>\n",
       "      <td>-0.040037</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3.141590</td>\n",
       "      <td>-2.978226e-04</td>\n",
       "      <td>-0.060081</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3.141586</td>\n",
       "      <td>-5.939895e-04</td>\n",
       "      <td>-0.080139</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      theta     theta_dot         u\n",
       "0  3.141593  0.000000e+00 -0.000000\n",
       "1  3.141593 -1.348867e-11 -0.020010\n",
       "2  3.141592 -9.955803e-05 -0.040037\n",
       "3  3.141590 -2.978226e-04 -0.060081\n",
       "4  3.141586 -5.939895e-04 -0.080139"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_u = pd.DataFrame(np.hstack([X, u.reshape(-1, 1)]), columns=x_names)\n",
    "X_u.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "threshold = 0.2\n",
    "model.fit(X_u, dXdt, threshold=threshold)\n",
    "\n",
    "model.n_params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['x1', 'x2', 'sin(x0)']"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.x_features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/billtubbs/anaconda3/envs/torch/lib/python3.10/site-packages/sklearn/base.py:458: UserWarning: X has feature names, but LinearRegression was fitted without feature names\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# R-squared\n",
    "model.score(X_u, dXdt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "y0    0.0\n",
       "y1    0.0\n",
       "dtype: float64"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Display coefficients\n",
    "model.intercept_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x1</th>\n",
       "      <th>x2</th>\n",
       "      <th>sin(x0)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>y0</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>y1</th>\n",
       "      <td>-1.0</td>\n",
       "      <td>0.5</td>\n",
       "      <td>-5.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     x1   x2  sin(x0)\n",
       "y0  1.0  0.0      0.0\n",
       "y1 -1.0  0.5     -5.0"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# True model equations \n",
    "# y0 = dxdt[0] = x[1]\n",
    "# y1 = dxdt[1] = -(b/m) * x[1] + (g/L) * np.sin(x[0]) + x[2] / L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.0"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "-params['b'] / params['m']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-5.0"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params['g'] / params['L']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1 / params['L']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>theta_dot</th>\n",
       "      <th>theta_ddot</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>-6.123234e-16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-1.348867e-11</td>\n",
       "      <td>-1.000483e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-9.955803e-05</td>\n",
       "      <td>-1.992160e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-2.978226e-04</td>\n",
       "      <td>-2.975509e-02</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      theta_dot    theta_ddot\n",
       "0  0.000000e+00 -6.123234e-16\n",
       "1 -1.348867e-11 -1.000483e-02\n",
       "2 -9.955803e-05 -1.992160e-02\n",
       "3 -2.978226e-04 -2.975509e-02"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dXdt.loc[0:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## B. Using identified model for prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>theta</th>\n",
       "      <th>theta_dot</th>\n",
       "      <th>u</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3.141593</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>-0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3.141593</td>\n",
       "      <td>-1.348867e-11</td>\n",
       "      <td>-0.020010</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3.141592</td>\n",
       "      <td>-9.955803e-05</td>\n",
       "      <td>-0.040037</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3.141590</td>\n",
       "      <td>-2.978226e-04</td>\n",
       "      <td>-0.060081</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3.141586</td>\n",
       "      <td>-5.939895e-04</td>\n",
       "      <td>-0.080139</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>3.141578</td>\n",
       "      <td>-9.873016e-04</td>\n",
       "      <td>-0.100208</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      theta     theta_dot         u\n",
       "0  3.141593  0.000000e+00 -0.000000\n",
       "1  3.141593 -1.348867e-11 -0.020010\n",
       "2  3.141592 -9.955803e-05 -0.040037\n",
       "3  3.141590 -2.978226e-04 -0.060081\n",
       "4  3.141586 -5.939895e-04 -0.080139\n",
       "5  3.141578 -9.873016e-04 -0.100208"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_u_sample = X_u.loc[0:5]\n",
    "X_u_sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/billtubbs/anaconda3/envs/torch/lib/python3.10/site-packages/sklearn/base.py:458: UserWarning: X has feature names, but LinearRegression was fitted without feature names\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>theta_dot</th>\n",
       "      <th>theta_ddot</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>-6.123234e-16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-1.348867e-11</td>\n",
       "      <td>-1.000483e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-9.955803e-05</td>\n",
       "      <td>-1.992160e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-2.978226e-04</td>\n",
       "      <td>-2.975509e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-5.939895e-04</td>\n",
       "      <td>-3.951002e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>-9.873016e-04</td>\n",
       "      <td>-4.919100e-02</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      theta_dot    theta_ddot\n",
       "0  0.000000e+00 -6.123234e-16\n",
       "1 -1.348867e-11 -1.000483e-02\n",
       "2 -9.955803e-05 -1.992160e-02\n",
       "3 -2.978226e-04 -2.975509e-02\n",
       "4 -5.939895e-04 -3.951002e-02\n",
       "5 -9.873016e-04 -4.919100e-02"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.predict(X_u_sample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'theta_dot': 0.0, 'theta_ddot': -6.123233995736765e-16}"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# For one-step ahead prediction use dicts or Series\n",
    "x_u = {\n",
    "    'theta': np.pi, \n",
    "    'theta_dot': 0, \n",
    "    'u': 0\n",
    "}\n",
    "model.predict(x_u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## C. Estimation from 10 data points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of parameters: 4\n",
      "Coefficients\n",
      "          x1        x2   sin(x0)\n",
      "y0  1.000000  0.000000  0.000000\n",
      "y1 -0.999997  0.499999 -4.999998\n",
      "Intercepts\n",
      "y0    0.0\n",
      "y1    0.0\n",
      "dtype: float64\n",
      "R-squared: 1.000\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/billtubbs/anaconda3/envs/torch/lib/python3.10/site-packages/sklearn/base.py:458: UserWarning: X has feature names, but LinearRegression was fitted without feature names\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "# The following data points were randomly sampled\n",
    "# from a time-series generated by simulating a \n",
    "# simple pendulum with the following equations\n",
    "# Pendulum system dynamics:\n",
    "# y[0] = dXdt[0] = x[1]\n",
    "# y[1] = dXdt[1] = -(b/m) * x[1] + (g/L) * np.sin(x[0]) + u(t) / L\n",
    "# with m=1, L=2, g=-10, b=1 and X0 = [pi, 0]\n",
    "\n",
    "X_u_data = np.array([\n",
    "    [-0.31337, -0.35066, -1.95775],\n",
    "    [ 0.43377,  0.1971 ,  0.12873],\n",
    "    [ 0.02058, -0.97899, -0.13313],\n",
    "    [-0.64254, -2.22694, -0.583  ],\n",
    "    [ 0.13985, -0.35462, -0.13526],\n",
    "    [ 0.29423,  0.4599 ,  1.98026],\n",
    "    [ 0.29177, -0.34668,  0.43937],\n",
    "    [-0.17303, -0.78595, -0.89596],\n",
    "    [ 0.3605 , -0.48347, -1.54148],\n",
    "    [ 0.33364, -0.20928,  0.92452]\n",
    "])\n",
    "\n",
    "dXdt_data = np.array([\n",
    "    [-0.35066,  0.91312],\n",
    "    [ 0.1971 , -2.23421],\n",
    "    [-0.97899,  0.80954],\n",
    "    [-2.22694,  4.93158],\n",
    "    [-0.35462, -0.40998],\n",
    "    [ 0.4599 , -0.91979],\n",
    "    [-0.34668, -0.87185],\n",
    "    [-0.78595,  1.19881],\n",
    "    [-0.48347, -2.05099],\n",
    "    [-0.20928, -0.9659 ]\n",
    "])\n",
    "\n",
    "x_names = ['theta', 'theta_dot', 'u']\n",
    "y_names = ['theta_dot', 'theta_ddot']\n",
    "X_u = pd.DataFrame(X_u_data, columns=x_names)\n",
    "dXdt = pd.DataFrame(dXdt_data, columns=y_names)\n",
    "\n",
    "custom_features = ['sin(x0)', 'cos(x0)', 'sin(x0)**2', \n",
    "                   'cos(x0)**2', 'sin(x0)*cos(x0)']\n",
    "\n",
    "model = SparseNonLinearModel(x_names, y_names, \n",
    "                             custom_features=custom_features,\n",
    "                             poly_order=2)\n",
    "\n",
    "model.fit(X_u, dXdt, threshold=0.2)\n",
    "print(f\"Number of parameters: {model.n_params}\")\n",
    "print(f\"Coefficients\\n{model.coef_}\")\n",
    "print(f\"Intercepts\\n{model.intercept_}\")\n",
    "print(f\"R-squared: {model.score(X_u, dXdt):.3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "torch",
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}