roughWork.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "747250f5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Defaulting to user installation because normal site-packages is not writeable\n",
      "Collecting deeptime\n",
      "  Downloading deeptime-0.4.4-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (2.1 MB)\n",
      "\u001b[2K     \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m2.1/2.1 MB\u001b[0m \u001b[31m1.7 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m[36m0:00:01\u001b[0m[36m0:00:01\u001b[0m:01\u001b[0m\n",
      "\u001b[?25hCollecting scipy>=1.9.0\n",
      "  Using cached scipy-1.10.1-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (34.4 MB)\n",
      "Requirement already satisfied: numpy>=1.20 in /usr/lib/python3/dist-packages (from deeptime) (1.21.5)\n",
      "Collecting scikit-learn>=1.1\n",
      "  Using cached scikit_learn-1.2.2-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (9.6 MB)\n",
      "Requirement already satisfied: threadpoolctl>=3.1.0 in /home/past/.local/lib/python3.10/site-packages (from deeptime) (3.1.0)\n",
      "Requirement already satisfied: joblib>=1.1.1 in /home/past/.local/lib/python3.10/site-packages (from scikit-learn>=1.1->deeptime) (1.2.0)\n",
      "Installing collected packages: scipy, scikit-learn, deeptime\n",
      "\u001b[31mERROR: pip's dependency resolver does not currently take into account all the packages that are installed. This behaviour is the source of the following dependency conflicts.\n",
      "tf-models-official 2.11.3 requires opencv-python-headless, which is not installed.\n",
      "easyocr 1.6.2 requires opencv-python-headless<=4.5.4.60, which is not installed.\u001b[0m\u001b[31m\n",
      "\u001b[0mSuccessfully installed deeptime-0.4.4 scikit-learn-1.2.2 scipy-1.10.1\n"
     ]
    }
   ],
   "source": [
    "!pip install deeptime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e2e31d8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d84e377a",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(0, 1, 100)\n",
    "x = 3 * np.exp(-2 * t)\n",
    "y = 0.5 * np.exp(t)\n",
    "X = np.stack((x, y), axis=-1)  # First column is x, second is y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "1229c913",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "from deeptime.sindy import STLSQ\n",
    "\n",
    "library = PolynomialFeatures(degree=3)\n",
    "\n",
    "optimizer = STLSQ(threshold=0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "875ae5b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from deeptime.sindy import SINDy\n",
    "\n",
    "# Instantiate the estimator\n",
    "estimator = SINDy(\n",
    "    library=library,\n",
    "    optimizer=optimizer,\n",
    "    input_features=[\"x\", \"y\"]  # The feature names are just for printing\n",
    ")\n",
    "\n",
    "# Fit the estimator to data\n",
    "estimator.fit(X, t=t);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "dcb7f199",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = estimator.fetch_model();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9c4cc722",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x' = -1.999 x\n",
      "y' = 0.999 y\n"
     ]
    }
   ],
   "source": [
    "model.print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "3bbe3a36",
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = 6\n",
    "y0 = -0.1\n",
    "\n",
    "t_test = np.linspace(0, 1, 100)\n",
    "x_test = x0 * np.exp(-2 * t_test)\n",
    "y_test = y0 * np.exp(t_test)\n",
    "\n",
    "sim = model.simulate([x0, y0], t=t_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "733d14fa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV8AAAE9CAYAAACsiWHQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA8Y0lEQVR4nO3dd3hUxfrA8e+bEGroNfSWUAIkQEAB6b1XkSKCqCjXK/YriAVUlOv12n5cVAQpihQpgqgoTQFBIPQSIXQCMYXQQ4Ak8/tjlyUJG9Jzssn7eZ59dufMnLPvSXlzMjtnRowxKKWUyl5uVgeglFJ5kSZfpZSygCZfpZSygCZfpZSygCZfpZSygCZfpZSyQD6rA8gJypQpY6pXr251GEqpXGbnzp2Rxpiyzuo0+QLVq1cnMDDQ6jCUUrmMiJxKrk67HZRSygKafJVSygKafJVSygLa55uMW7duERISQkxMjNWhqHQqWLAglStXxsPDw+pQlLqLJt9khISEULRoUapXr46IWB2OSiNjDOfPnyckJIQaNWpYHY5Sd9Fuh2TExMRQunRpTbwuSkQoXbq0/ueiciy98r0HTbyuTb9/KlPcuAG//ALLl0NYGJQvD/37Q9euUKBAug9ryZWviJQSkTUiEmx/LplMu69EJFxEDqR2fxGZICJHReSwiHTN6nPJSp6enim2efzxxzl06BAA7777bqK6li1bZsp7ZJZRo0axZMkSIONxK5UtbtyA11+H//0Pbt7EVKqMuXHTVn79dVt9OlnV7TAeWGeM8QbW2cvOzAG6pXZ/EakPDAF87ftNFxH3zA3diRs3YOVKePRR6NHD9rxyZYa+Mak1c+ZM6tevD9ydxLZs2ZLl759erhq3ylvM6l+ICjzOfs/7+SmqMV8fbcHf7hWhWjXYu9d2RZxOViXfvsBc++u5QD9njYwxG4GoNOzfF1hojLlhjDkBHAWaZ07IyUjyl5HKlW3PmfCX8bbffvuNdu3aMWjQIOrWrcvw4cO5vQJJu3btCAwMZPz48Vy/fh1/f3+GDx8O3LmqvXr1Kh07dqRJkyY0bNiQFStWpPie8+bNo1GjRvj5+TFixAgATp06RceOHWnUqBEdO3bk9OnTgO2Kdty4cbRs2ZKaNWs6rm6NMfzzn/+kfv369OzZk/DwcMfxUxO3MYaXX36ZBg0a0LBhQxYtWpTi10OpjDAGoqJg/37Dl18eYtSLi+iy9xStfviYeUf2EhPnQWh0cRCBMmVsXRHpfzOT7Q/gYpLyhXu0rQ4cSM3+wDTg4QTbZwGDkjnuGCAQCKxatapJ6tChQ3dtc2rFCmO6dDHmiSeMGTPmzuOJJ2zbV6xI3XGcKFKkiDHGmA0bNphixYqZM2fOmLi4OHP//febTZs2GWOMadu2rdmxY0ei9kn3v3Xrlrl06ZIxxpiIiAhTq1YtEx8f73QfY4w5cOCA8fHxMREREcYYY86fP2+MMaZXr15mzpw5xhhjZs2aZfr27WuMMWbkyJFm0KBBJi4uzhw8eNDUqlXLGGPM0qVLTadOnUxsbKw5e/asKV68uPnuu+9SHfeSJUsc+//999+mSpUq5ty5c/f8eiSV6u+jyrMuXDDmwIF489VXf5lHHplumjZ90BQtWs4AiR7+pQPMF62/Nj90/T/b7/jjjxvTo8c9jw0EmmRyW5Z94CYia4EKTqomZtV7As4+YXF6SWSMmQHMAAgICEj/ZdPy5ba/gEk/3En4l7FPn3Qf/rbmzZtTuXJlAPz9/Tl58iQPPPBAqvY1xvDqq6+yceNG3NzcOHv2LGFhYVSo4OzbA+vXr2fQoEGUKVMGgFKlSgGwdetWli1bBsCIESP417/+5dinX79+uLm5Ub9+fcLCwgDYuHEjQ4cOxd3dnYoVK9KhQ4c0nfPmzZsd+5cvX562bduyY8cOihUrlqGvh8rbrl6Fs2dtj82bA/nxx485cmQDFy+eu+d+56JDMMZw4UZhjAGJjoZy5dIdR5YlX2NMp+TqRCRMRLyMMaEi4gWEJ9c2GcntHwJUSdCuMnDvr2hGhYXZuhqcKVwYzmXO2xdI8Kmqu7s7sbGxqd53/vz5REREsHPnTjw8PKhevfo9h2AZY1I1UiBhm4TxmQRdABkZcZDwOEll5Ouh8pYbN2y/hkeOXObixWJcunSnLirqAtu3z3e6X5EipfAu50/960V5oFYdWpQvRGXP3ZQucA3BQGQkPP10uuOyqs93JTDS/nokkHInZOr2XwkMEZECIlID8Aa2ZzDWeytfHq5dc16Xwb+MaeXh4cGtW7fu2n7p0iXKlSuHh4cHGzZs4NSpZCdaAqBjx44sXryY8+fPAxAVZet2b9myJQsXLgRsCT2lK802bdqwcOFC4uLiCA0NZcOGDWmKu02bNixatIi4uDgiIiLYuHEjzZtnbRe+cn3x8fD33/DHHzG89dY6+vZ9hXbtmtKrV1nCwq4kalu7divc3W13QBYqVBw/vz4MHfox06fvYfPmCNZs+olvHvZhbJE9+BcMpkz+K8i1q3DqFPj52YabpZNV43ynAotF5DHgNPAggIhUBGYaY3rYywuAdkAZEQkB3jTGzEpuf2PMQRFZDBwCYoGnjTFxWXom/fvbPlwrUiRx14PJ+F/GtBozZgyNGjWiSZMmzJ9/56/58OHD6d27NwEBAfj7+1O3bt17HsfX15eJEyfStm1b3N3dady4MXPmzOHTTz9l9OjR/Oc//6Fs2bLMnj37nsfp378/69evp2HDhvj4+NC2bds0xd2/f3+2bt2Kn58fIsL7779PhQoV+Ouvv9LwVVF5weXLcOaMYdOmA/z22xr27/+V4OCN3Lp1PVG74OCNNGzY01HOn78wo0fPxte3Dq1bN6ZKFfckvYgF4J2374zzPXfOdkH19NMZHucr9/rXLq8ICAgwSefzDQoKol69einvfHu0w969tj7ewoVtV7yRkba/jG+/naFvkMqYVH8flUuJjbXlwZAQ+P77H9iy5XsOHlx9z35bNzd3Bg78gC5dnqNcOVtvYcWKtlzqlkV9ACKy0xgT4KxO73DLqAIFbAk2C/4yKqXuuHwZTp+GM2dsv2Zx9v9p162by65dS53uU768D/XqdaF588507twOH59ieHlBTphrSZNvZihQwDaiIRNGNSilbG733QYFXeann9awbdtPFClShoED/52ona9vN0fyLVy4BHXrdsLPrwudO3emadPqVK5s+4c0p9Hkq5TKMWJibFe2W7Yc5ZdfVrF79yqOHPmd+HjbaJZixSowYMDURKNofH2706vXG7Rp05X27ZtTo0Y+SpW6e/RnTqPJVyllqcuX4dixOH75ZSsbNqxg794fCAs7nEzbvzlzZjdVqzbB0xOqVoUqVSrxwguTc0RXQlpo8lVKZStjICICTp60jdgKDb3Ka6/V4sqV5If7V6nSmEaNetKuXU86dvSjRg0oUSLbQs4SmnyVUlnudv/t3r1RnDplcHMr7agrWNCT0qWrJUq+Hh6FqFevE02b9qZ79x4EBFSiUiXIn9+K6LOGTqaeg7m7u+Pv7+94TJ06NdOOvWfPHn766adMOVZqpqX8+OOPiY6OdpR79OjBxYsXM+X9Vc4UH28bnbB06TlGjvwfHTp0pFevcvz447S72vr59aVYsQq0bj2G8eNXsW7dedasWcmsWU/w4IOVqFEjdyVe0HG+QAbH+WYhT09Prl69miXHnjNnDoGBgUybdvcvQlqlJs7q1asTGBjomC8iu+SE72NeEhdn+8Dszz9Ps2LFUnbsWMqxY38kalO5sh+vv77HURaBUqWu4+1dgBo13ChaNJuDzkL3GuerV74u5tKlS9SpU4fDh20fSAwdOpQvv/wSgLFjxxIQEICvry9vvvmmY58dO3bQsmVL/Pz8aN68OZcuXeKNN95g0aJF+Pv7O6ZqvO3gwYM0b94cf39/GjVqRHBwMAAffvghDRo0oEGDBnz88cd3xfbbb7/Rq1cvR/mf//yn4864c+fO0b59e9q3bw/YknFkZGSyxz158iT16tXjiSeewNfXly5dunD9+vW73lNZLy7OdoW7cOFphg37L+3b38fQodVYuPCFuxIvgIdHQeLirlGjBrRvD488AgMHFqJRo9yVeFOifb4pmDEja48/Zkzydbfnub1twoQJPPTQQ0ybNo1Ro0bx7LPPcuHCBZ544gkApkyZQqlSpYiLi6Njx47s27ePunXr8tBDD7Fo0SKaNWvG5cuXKVy4MG+99VayV76ff/45zz77LMOHD+fmzZvExcWxc+dOZs+ezbZt2zDGcN9999G2bVsaN26c4jmOGzeODz/8kA0bNtx15ZvccUuWLElwcDALFizgyy+/ZPDgwSxdupSHH344dV9YlaXi4203Ohw9avvg7Pjxfbz9tp/TtiJu+Pi0o3nzgfTu3ZdmzSpRuTLky+PZJ4+ffs5WqFAh9uzZc9f2zp0789133/H000+zd+9ex/bFixczY8YMYmNjCQ0N5dChQ4gIXl5eNGvWDIBixYql+L4tWrRgypQphISEMGDAALy9vdm8eTP9+/enSJEiAAwYMIBNmzalKvneS3LH7dOnDzVq1HD88WnatCknT57M0HupjDHGNolfYGAkoaGFMKaIo65ixQaULFmFCxfOAODmlo+6dTty//2D6NevL02alMXLK+tu43VFmnxdUHx8PEFBQRQqVIioqCgqV67MiRMn+OCDD9ixYwclS5Zk1KhRxMTEpHp6yISGDRvGfffdx48//kjXrl2ZOXNmqlaKyJcvH/Hx8Y5yalYOTsu0kdrtYI2oKNi37yqLF3/Ppk3fcujQGh5++AtatRrtaOPm5kbz5sM4c2YPLVsOZsCAfvj7l6JChZx/s4NV9O+QC/roo4+oV68eCxYsYPTo0dy6dYvLly9TpEgRihcvTlhYGD///DMAdevW5dy5c+zYsQOAK1euEBsbS9GiRbly5YrT4x8/fpyaNWsybtw4+vTpw759+2jTpg3ff/890dHRXLt2jeXLl9O6detE+1WrVo1Dhw5x48YNLl26xLp16xx1yb1fao6rsl90NOzaFcvEiT/TvfvDdO1ans8+G8GBAz8THx/Ljh0LErUvVAgmTnyPtWtXM2PGaLp3L4WXlybee9Er3xTcq082qyXt8+3WrRujR49m5syZbN++naJFi9KmTRveeecdJk+eTOPGjfH19aVmzZq0atUKgPz587No0SKeeeYZrl+/TqFChVi7di3t27dn6tSp+Pv7O/qSb1u0aBHffPMNHh4eVKhQgTfeeINSpUoxatQox3y6jz/++F1dDlWqVGHw4ME0atQIb2/vRPVjxoyhe/fueHl5JZrXt0mTJk6Pq10M2S82Fk6cMKxevZcVK+ayffu3yd74cOtWDG5ut/Dx8aBWLexdCppp00KHmpFzh5qpjNPvY8r+/huOHIHNm/9i+vTBnD2732k7L6/6tGjxMAMGDOWBB2wT1rhn/drgLk2nlFRKJRIdbUu4hw/jWFanaNEqREaeSNSueHEvmjcfRu/eD9O1qx81aojLzaGQU2nyVSqPiI+33QCxevVBvvtuJk2aDKJ27VaO+gIFitCkySACAxfRuPEAunR5hIEDO+Lj454jp2R0dZp8lcrlLl+GXbuuMn/+Ytat+5ITJ/4E4Nq184mSL8DQoVN4++1PaNy4GPZFq1UW0eSrVC4UH2+bMWzVqj0sWfI527bN58aNxLeA79y5hCFDplGkSDGqVYM6daBy5Yo6FjebaPJVKhe5dg12745m3rxFrF37BSdObLurjbu7B/7+/ejW7THatStCnTpQsKAFweZxmnyVygXOnYODB+HgwQtMnFiL6OgLd7UpX74Obds+wfDhI2jRohxly1oQqHLQfzBysClTpuDr60ujRo3w9/dn2zbbVUy7du24PTSuevXqDBw40LHPkiVLGDVqFGCbuaxs2bI0btwYb29vunbtypYtW7L9PG6bNGkSH3zwQYbbKJtbt+DQIVi8GFatghMnoHDhklSr1tTRJl++/DRvPox33vmdDRuC+OabF+nTRxNvTqBXvjnU1q1bWbVqFbt27aJAgQJERkZy8+ZNp20DAwM5ePAgvr6+d9XdnogHYMOGDQwYMIANGzbo2FcXdvky/PnnRWbP/op8+YrSuvUTierbtXuayMgTtG//FA8/PIoWLcqQzTN5qlTQK98cKjQ0lDJlyjjmNyhTpgwVK1Z02vall17i3XffTfGY7du3Z8yYMcyYMYNjx47RpEkTR11wcDBNmza9a5927drx/PPP06ZNG+rVq8eOHTsck+289tprjnbJTTc5ZcoU6tSpQ6dOnRzTYAIcO3aMbt260bRpU1q3bs1ff/2VYvx5XWgozJp1mD59/kHv3pVYuPBFVq2aTFzcrUTt2rbtw/r1R/jmm5fo3VsTb06lV745VJcuXXjrrbfw8fGhU6dOPPTQQ7Rt29Zp28GDBzN9+nSOHj2a4nGbNGnCF198Qa1atShevDh79uzB39+f2bNnO7orksqfPz8bN27kk08+oW/fvuzcuZNSpUpRq1Ytnn/+eU6ePOl0Wsj4+HgWLlzI7t27iY2NpUmTJo4EP2bMGD7//HO8vb3Ztm0b//jHP1i/fn26v165VXw8HD1q+Pbb9Sxb9hH79/+YqP7ixbPs2fM9zZo9SPXq4OsLXl56TeUK9LuUSpMmTUJEUvUY42RCiDFjxiRqM2nSpHu+n6enJzt37mTGjBmULVuWhx56iDlz5jht6+7uzssvv8x7772X4nkkvJ388ccfZ/bs2cTFxbFo0SKGDRvmdJ8+ffoA0LBhQ3x9ffHy8qJAgQLUrFmTM2fOJJoW0tPT0zEt5KZNm+jfvz+FCxemWLFijuNcvXqVLVu28OCDD+Lv78+TTz5JaGhoirHnJTdvws6dt3j66W/o2LExkyd3uivxVqrUkEcf/ZIRI3oydCh06mSbY0G5Br3yzcHc3d1p164d7dq1o2HDhsydOzfZq9MRI0bw3nvvOe33TWj37t2O/t6BAwcyefJkOnToQNOmTSldurTTfW53fbi5uSWa5tHNzY3Y2Nh7TgvpbDrL+Ph4SpQo4XSu4rwuOhr274cvvpjF999P4sKFkET1IkLDhr0YMOB5hgxpR+3aovMruChLrnxFpJSIrBGRYPtzyWTafSUi4SJyIMn2/4jIXyKyT0SWi0gJ+/bqInJdRPbYH59nw+lkicOHDzuW7wHbgpfVqlVLtr2HhwfPP/+80+V9bvv999+ZMWOGY+WLggUL0rVrV8aOHcujjz6a7liTmxayTZs2LF++nOvXr3PlyhV++OEHwDahe40aNfjuu+8A29V4wknh86LLl2HTJvj2W9i7Fy5cCE+UePPnL0zbtv/gyy//4uefV/Lmm+2pU0cTryuz6sp3PLDOGDNVRMbby684aTcHmAbMS7J9DTDBGBMrIv8GJiTY/5gxxj+zA540aVKKXQX3MmPGDGakYU2iq1ev8swzz3Dx4kXy5ctH7dq1U9z/scce45133km0bdGiRWzevJno6Ghq1KjB0qVLE410GD58OMuWLaNLly5pO6EEkpsWEmyjLfz9/alWrVqieXrnz5/P2LFjeeedd7h16xZDhgzBz8/5MjS5WVQUbNgQQVRUWRL+A9GmzVP8/PMU8ucvQseOz/D442N54IHSFC9uXawqc1kypaSIHAbaGWNCRcQL+M0YUyeZttWBVcaYBsnU9wcGGWOGp9Q2OXl5SskPPviAS5cu8fbbb1sdSpbIqd/HiAj44YcjzJz5Htu3z+e11/ZQsWL9RG1CQrbSvXtjGjcuqHeguaicOKVkeWNMKIA9AZfLwLFGAwmX360hIruBy8BrxphNznYSkTHAGICqVatm4O1dV//+/Tl27JiOMshG4eGwdOkBZs+eQmDgIkd/+S+//JtHH50LQNGi4OcHPj4t8vwik7lZln1rRWQtUMFJ1cRMfI+JQCww374pFKhqjDkvIk2B70XE1xhzOem+xpgZwAywXflmVkyuZPny5VaHkGfcTrqzZr3Nrl3f3fUh5aVLoRQvHkfTpu7UrKkLTeYFWZZ8jTGdkqsTkTAR8UrQ7eB8rZJ7EJGRQC+go7H/JBtjbgA37K93isgxwAcITPZASmWhyEhYuvQQX3452WnSbdCgO0OHTmTYsFZUq6ZrnuUlVv1TsxIYCUy1P69Iy84i0g3bB2xtjTHRCbaXBaKMMXEiUhPwBo6nN8j0rPyrcg4rl8iKioI//ohmypQn2b59/l2x+Pn1YcSI13nooQAqV7YoSGUpq5LvVGCxiDwGnAYeBBCRisBMY0wPe3kB0A4oIyIhwJvGmFnYRkAUANbYk+OfxpingDbAWyISC8QBTxljotITYMGCBTl//jylS5fWBOyCjDGcP3+egtn8SdWVKxAYCMHBYEwhwsODEyVeP78+PPLIGwwZ0pRk7hZXeYQuoInz0Q63bt0iJCSEmJgYi6JSGVWwYEEqV66MRzYsOnb9OuzaZQgKEuLj72wPClrLxx93pkGD7owa9RZDhwZo0s1D7jXaQZMvzpOvUqkRGwvbtkXz0UefsGfPKl56aSNubnfufDDGcO3aPgYP9tPuhTwoJw41U8qlGQOHDsXxySfz+O6717l48SwA27Z9Q4sWIwEoWRKaNROqV897N4+olGnyVSqNzp6FL75Yy1dfvcDZs/sT1W3fvoDOnUcSEADe3jp6QSVPk69SqXTxIixeHMS0aS/fNcNYsWLlGTBgMuPGPUajRuicCypFmnyVSsHNm/D77xf5978nsWHDNOLj4xx1+fMXpmvXl3nhhZdo2dKT/PktDFS5FE2+SiXDGDhyBDZtiubll325ePGco05EaNFiFM899w7du1fE09PCQJVL0uSrlBORkbB5s+22YChMQMAQ1q79EAAfn7aMHfsRw4c31oUoVbpp8lUqgRs3YOvWWxw5knhscM+er3P06Cb69PkX48YNpFYt/SRNZYxO36GU3eHD8YwdO4OBA2tz4cLZRHXFipXgu++28d//DtLEqzKFXvmqPO/SJZg3bx8ff/wUx49vBWDp0pd4/PEFAFSvDi1bgqenJl2VeTT5qjwrPh62bbvOW29N5tdfP0g0iuHEie14eFyic+fiemeayhKafFWeFBEB06at57PPxhARccyx3d3dg+7dX2HChFe5775COl5XZRlNvipPiY2F33+/xKRJL7F588xEdd7ebXjhhS8YOrSurpWmspwmX5Vn/P03fPrpr3z22WjHXAwAhQoVZ8iQ//DKK49Rp45+Bq2yhyZflevFxcGOHbBvH0RFXUmUeBs3HsCrr06jVy8vXaRSZStNvipXi4iADRts8zIANGkykCZNBhEcvJFHH/0fL700iCpVLA1R5VGafFWuFB8PW7de47ffzlG2rHeiumHDpuPjI3TtWkbnYlCW0eSrcp3Ll+Hzz7fz4YcjMMbw2mu7KVCgCACentCzZ1kqVbI4SJXn6acLKlcJCopj1KgpTJjQkrCwI4SHB7N06csA1KkDgwahiVflCHrlq3KFmzdh2bIzvPHGwwQHb3RsL1iwKLVrN6NrV6hWzcIAlUpCk69yeRERMHXq93z++Wiioy84tteq1YrXX/+awYNrUKiQhQEq5YQmX+XSAgNj+Ne/XmbDhmmObSJu9O79Bu+8M5GGDfVHXOVM+pOpXNLNm7BgwXHeeONBTp/e5dheqlRVnn/+W55+uhUlS1oYoFIp0OSrXE5kJKxdC8uWzU+UeP39+/Huu7Po0qWUzsmgcjwd7aBcyuHDsGKFbThZ9+6vUrduR9zdPRg27BOWLl1G9+6aeJVr0Ctf5RLi4mDLFggKurPNzc2dxx77lvj4U4wd24yiRa2LT6m0suTKV0RKicgaEQm2PzvtnRORr0QkXEQOJNk+SUTOisge+6NHgroJInJURA6LSNesPheV9aKj4cMPdzBhwmOJ5twFuP/+crz4oiZe5Xqs6nYYD6wzxngD6+xlZ+YA3ZKp+8gY429//AQgIvWBIYCvfb/pIqL/hLqw8HB44YV5vPpqa/744yt++GESAPnyQYcO8MADaDeDcklWJd++wFz767lAP2eNjDEbgag0HnehMeaGMeYEcBRonoE4lYUOHYpl+PAX+OKLkcTG3gDg99+nI3Kefv2gdm1r41MqI6xKvuWNMaEA9udy6TjGP0Vkn71r4na3RSXgTII2IfZtyoUYA+vXX+LBB3uxdu1Hju0VK/ry6afbeeSR0pQqZWGASmWCLEu+IrJWRA44efTNhMN/BtQC/IFQ4L+339ZJW5NMfGNEJFBEAiMiIjIhJJUZYmNh7txjPPxwCw4d+sWx3d+/H19/vZVHH61FgQIWBqhUJsmy0Q7GmE7J1YlImIh4GWNCRcQLCE/jscMSHOtLYJW9GAIknJ21MnAumWPMAGYABAQEOE3QKntFR8PHH2/h3Xf7cO3aecf2Xr1e47//nYyPj46MVLmHVT/NK4GR9tcjgRVp2dmesG/rD9weDbESGCIiBUSkBuANbM9grCobREXBK68s4s03OzgSb758BXjqqfnMmvW2Jl6V61g1zncqsFhEHgNOAw8CiEhFYKYxpoe9vABoB5QRkRDgTWPMLOB9EfHH1qVwEngSwBhzUEQWA4eAWOBpY0zisUkqxzl3Dlatus6CBeMdH6wVLVqW8eNXMG5cCzw9LQ5QqSwgxuh/3AEBASYwMNDqMPKkY8dsy/zEx0NoaBDvv9+KYsXK8+67P/LwwzXx8LA6QqXST0R2GmMCnNXpHW7KMgcO2O5au83Lqx7PPbeGBx6oSdeuJXHTngaVi2nyVZbYsOEyv/4aRI0a9yXaPnhwU/z8LApKqWyk1xYqWxkDK1eGM2pUez76qCMnT+4AwM3NdseaJl6VV2jyVdkmPh4WLjzDmDFtOH16FzduXOP//q8HN29epFs3vWNN5S3a7aCyRVwczJt3lJde6khU1GnAtuLE4MFTGTy4BOXSc4+jUi5Mk6/KcrGxMHNmEK+80oHLl/8GwN3dg6efXsDkyQMpUcLa+JSygiZflaViY+Gzz/YxcWInrlyx3cbt4VGIF15YxsSJ3XQqSJVnafJVWSY2FqZP38urr3Z03LVWoIAnEyb8yMsvt6FwYYsDVMpCmnxVlriTeDtw7ZptVtCCBYvx2murefHFFhQsaHGASllMk6/KdLGx8PPPsGbNckfiLVSoOJMmrWHcuGaaeJVCk6/KZHFx8OuvEBoKvXq9yY0bV9m8eSaTJq3l2WcDdDpIpex0bgd0bofMEh9vS7ynT9/ZZvv5OsOoUVU18ao8R+d2UFnOGPj++0jCw0vhlmBShjJlhF69NPEqlZTe4aYyxYoVYfzjH62YM2ckcXG3AChZEnr2RBOvUk7ola/KsPXrLzFuXDfCwo4QFnaE2NgbvPTSYnr2RD9cUyoZeuWrMiQwMIaxY/ty5swewHbL8AMPDKFXL3Qcr1L3oMlXpdvhw3E88cQwjhz53bFt9OiZTJo0QFefUCoFmnxVupw9a3jqqWfYs2e5Y9uDD77Pe+89qnM1KJUKmnxVmkVFwXPPvc9vv33m2Na58wt89NHLlC1rYWBKuRBNvipNoqPh9dcXsGTJeMe2Zs2G8L///YdKlSwMTCkXo8lXpdqtW/DRR5uZMWOUY5uPT1umT5+Dt7f+KCmVFvobo1LFGFi/Hq5fL4KnZxnAtuDlxx8vJyBAB/IqlVaafFWq7NgBp05B1aqNGT9+G/Xrd+Hdd3+kW7eSVoemlEvSmyxUioKDYc+eO+WSJSvz1lu/0KcPiFgWllIuTa981T2Fh8PPP19ItK1QIejaFTw8LApKqVxAk69KVnQ0TJo0n1dfrU1Q0FoA3N1tiVdvolAqYzT5Kqfi4+Hzz3cyc+bjXLsWxSefdGXv3pW0aYOuNKxUJrAk+YpIKRFZIyLB9menn9qIyFciEi4iB5JsXyQie+yPkyKyx769uohcT1D3eTacTq60enUkU6YM4NatGAAqVKjDgAHt8Pa2ODClcgmrrnzHA+uMMd7AOnvZmTlAt6QbjTEPGWP8jTH+wFJgWYLqY7frjDFPZW7YecPhw3G89NIwoqJss6IXLFiMSZO+p337YhZHplTuYVXy7QvMtb+eC/Rz1sgYsxGISu4gIiLAYGBBJseXZ124AC+/PJmgoDWObf/4x9c88ogPbtpJpVSmserXqbwxJhTA/pzeXsTWQJgxJjjBthoisltEfheR1hkNNC+5dQv+859fWLXqHce2nj0nMmFCH52XV6lMlmXjfEVkLVDBSdXETHyboSS+6g0FqhpjzotIU+B7EfE1xlx2Et8YYAxA1apVMzEk17V8eQiffjqc2+v61a3bkX//ezJlylgcmFK5UJYlX2NMp+TqRCRMRLyMMaEi4gWEp/X4IpIPGAA0TfCeN4Ab9tc7ReQY4APctTqmMWYGMANsC2im9f1zm0OHYpk4cRjXrp0HoHhxL6ZMmY+vr7vFkSmVO1nV7bASGGl/PRJYkY5jdAL+MsaE3N4gImVFxN3+uibgDRzPYKy53sWLMG3aCo4e3QTYVqN4/vmF9OlT3trAlMrFrEq+U4HOIhIMdLaXEZGKIvLT7UYisgDYCtQRkRAReSzBMYZw9wdtbYB9IrIXWAI8ZYxJ9gM7BXFxsG4d+PsPZNSoOeTPX5i+fSczblwb8unN50plGbndv5eXBQQEmMDAu3om8oStW2H//jvliIjj9OtXTbsblMoEIrLTGBPgrE6vbfKwkJDEiRfg/vtr4utrTTxK5SU6cjOPiomBmTP/5ObN645tnp7QWgfnKZUtNPnmUd99d4J//7sL777blNOndyECHTpA/vxWR6ZU3qDJNw86ciSet98eRUzMFUJDg5g9+xEaNYqngrNR2UqpLKHJN4+xTRP5KcHBGwFwc3PnmWdm0ayZ/igolZ30Ny6PWbDgCEuWTHCUe/SYwJNP3qfzNiiVzfRXLg/566843nvvUcc0kZUr+zFp0uuU1GXYlMp2mnzziOvX4e23p3Hs2BYA3Nzy8eKLc2nSRD9hU8oKmnzziCVLTrBkyauOcs+eExk1yk8XwFTKIikmXxH5Z3IrTSjXcPKk4b33nuTmzWgAKlZswBtvvEqJEtbGpVRelpor3wrADhFZLCLd7BOYKxdx6xZ89NEyx+ToIm6MGzdLuxuUsliKydcY8xq22cFmAaOAYBF5V0RqZXFsKhPs2AE+Pr3o3Xsy+fLlp0OHcTz+eHMd3aCUxVI1t4MxxojI38DfQCxQElgiImuMMf/KygBV+kVGwsGD4OFRgF693iAgYDCtW1emdGmrI1NKpZh8RWQctjl3I4GZwMvGmFsi4gYEA5p8cyBjYNMm2/Nt3t51de4GpXKI1Fz5lgEGGGNOJdxojIkXkV5ZE5bKqKAgQ0RE4u75Vq3QOXqVyiFS0+f7RtLEm6AuKPNDUhkVEwNTpnzKtGm9CA8/CkCNGqBL1SmVc+jHLrnQ6tV/s2zZG+zf/yOTJ/ty/PgmWra0OiqlVEKafHOZyEh4//0JxMTYFmwuXboGAwfeR5EiFgemlEpEk28uM3fuDrZuneMojx79KU2b6phepXIaTb65yNGjhv/9b5yj7OfXl7Fju+iYXqVyIP21zCViY+Hjjxdy4sSfAOTLl59x4z6gShWLA1NKOaXJN5cIDLzOggXjHeWOHZ9l4MDaFkaklLoXTb65wPXr8MEHnxAVdRoAT88yvPjiRIoXtzgwpVSyNPnmAuvXR/Ljj+85ygMGvEWbNpp5lcrJNPm6uMuXYdGiZY6hZeXL1+HZZ5+gQAGLA1NK3ZPebOritm+HBx4YQ7lydVi69GUGDJiIn59+W5XK6fS31IVFRsLx47bXPj5tGT9+Gx06gLu7tXEppVJmSbeDiJQSkTUiEmx/vmulDBGpIiIbRCRIRA6KyLOp2V9EJojIURE5LCJds+ucrLB9e+Jy2bJC7do6171SrsCqPt/xwDpjjDewzl5OKhZ40RhTD7gfeFpE6t9rf3v9EMAX6AZMF5FceR0YGgqnT8cl2tasGbomm1Iuwqrk2xeYa389F+iXtIExJtQYs8v++goQBFRKYf++wEJjzA1jzAngKNA8C+K33Pz523nttdr88cdXxMXF4uWF3lChlAuxKvmWN8aEgi3JAuXu1VhEqgONgW0p7F8JOJNg1xDuJOxcIyQEZs9+nfPnTzJv3mMsWfISzZpZHZVSKi2y7AM3EVmLbfHNpCam8TiewFLgOWPM5ZSaO9lmnGxDRMYAYwCquthEt19//QeHDv0K2BbEfOihsVRw9pVWSuVYWZZ8jTGdkqsTkTAR8TLGhIqIFxCeTDsPbIl3vjFmWYKq5PYPARL+810ZOJdMfDOAGQABAQFOE3ROFBICX389yVG+776H6devjnUBKaXSxapuh5XY1oXD/rwiaQP7EvWzgCBjzIep3H8lMERECohIDWyrLicZE+DavvnmD4KC1gLg5ubOE0+8Trl7dtoopXIiq5LvVKCziAQDne1lRKSiiPxkb9MKGAF0EJE99kePe+1vjDkILAYOAauBp40xiYcEuLCzZ+Hrr99ylJs3H07v3jp5jlKuyJKbLIwx54GOTrafA3rYX2/GeR9usvvb66YAUzIt2Bxk0aIdifp6H3tsImXLWhyUUipddG4HFxEWBvPm3fmb0qzZEPr08bEwIqVURmjydRHLlx9g7947XeOPPDJB+3qVcmGafF1AVBTMm/eBo+zv34/+/RtYGJFSKqN0Yh0XsGcP9O8/leLFvfjtt+kMHTqeihWtjkoplRGafHO4q1fh2DEoXrwC/fu/R48er9G7t64Dr5Sr026HHG7/fjAJbgEpX74I1apZF49SKnNo8s3Bbt6Ev/5KvM3PT2cuUyo30OSbg+3de4PFiycQGXkCgEKFoLbeU6FUrqDJN4eKj4eZMxeyevVUXnutNosWPUeDBrpKhVK5hSbfHOr4ccNPP30EgDHxlCxZgfr1U9hJKeUyNPnmUAsX/k5IyF4APDwKMXr0GF2RWKlcRJNvDhQRAUuXfuoot2gxkpYtS1kYkVIqs2nyzYHWrj2V6FbiYcOeoUQJ6+JRSmU+Tb45zPXrMG/eZxgTD0C9ep3o3l07e5XKbTT55jB798awefMsR7lHj2d0YUylciFNvjmIMTB37ndcvRoJQOnS1Rg+vKfeVKFULqTJNwc5fRp++WW6o9yu3VPUr68De5XKjTT55iDbt1/k+vWLALi7e/Dww6N1eJlSuZTOapZDXL0KFy6UYNKkQwQHb+Ts2f20aqWzpSuVW2nyzSFuT6AjIvj4tKVly7a6PptSuZh2O+QA8fFw+HDibfXqWROLUip7aPLNAUJC4Nq1O+V8+XT2MqVyO02+OcDq1QeZOrUFmzfPIibmKrVrg4eH1VEppbKS9vla7Pp1WLLkK06c+JMTJ/7kr7/WMWTIt1aHpZTKYnrla7FDh27x55/fOMqdOo3QJeGVygM0+Vps8eKfuXIlHIASJSry4INdLI5IKZUdNPlaKDISfv11jqN8//2PUKeO3tGmVF5gSfIVkVIiskZEgu3PJZ20qSIiG0QkSEQOisizCer+IyJ/icg+EVkuIiXs26uLyHUR2WN/fJ6Np5VmO3acZ//+VY7ywIGjKFzYwoCUUtnGqivf8cA6Y4w3sM5eTioWeNEYUw+4H3haRG7PrbgGaGCMaQQcASYk2O+YMcbf/ngq604hY+LjYeHCRcTF3QKgRo376NChjsVRKaWyi1XJty8w1/56LtAvaQNjTKgxZpf99RUgCKhkL/9qjIm1N/0TqJzVAWe2M2dg06Z5jnKrViOoVs3CgJRS2cqq5FveGBMKtiQL3PPzfRGpDjQGtjmpHg38nKBcQ0R2i8jvItI6k+LNdBs2BHPihO103NzyMXjwQ7oysVJ5SJaN8xWRtUAFJ1UT03gcT2Ap8Jwx5nKSuonYuifm2zeFAlWNMedFpCnwvYj4Jt3Pvu8YYAxA1apV0xJSht28CcuX3xnL27BhD5o1K5OtMSilrJVlydcY0ym5OhEJExEvY0yoiHgB4cm088CWeOcbY5YlqRsJ9AI6GmOM/T1vADfsr3eKyDHABwh0Et8MYAZAQECASccpptvJkxAfbyhYsBgxMZdp02Y45ctnZwRKKatZ1e2wEhhpfz0SWJG0gYgIMAsIMsZ8mKSuG/AK0McYE51ge1kRcbe/rgl4A8ez5Awy4Ngx6N17Eh98EMaTTy5h0KDeulqFUnmMVcl3KtBZRIKBzvYyIlJRRH6yt2kFjAA6JBg61sNeNw0oCqxJMqSsDbBPRPYCS4CnjDFR2XROqRITY5tIB8DDoyBNmgykQYNC1gallMp2lsztYIw5D3R0sv0c0MP+ejPg9HrQGON0zi9jzFJs3RQ51okTtrXabitZEkqVsi4epZQ19A63bHbsWOKyTh2pVN6ks5plo+homD59IufPnyIgYDD163ehVq2CVoellLKAXvlmo+PH49myZTbbt89n+vS+hIdvplgxq6NSSllBk282WrXqDy5dCgWgaNGy9O7dztqAlFKW0eSbTa5fhzVrljjKjRsPwMdHe32Uyqs0+WaT48fj2b37zn0ibdsO1C4HpfIwTb7Z5JdfdnDhgm2Ab5EipejZs521ASmlLKXJNxvcuAG//nrnqrdRoz7UqaMrZCqVl2nyzQanThl2717uKLduPYASJayLRyllPU2+2eD33w8RHh4MQIECRejVq7PFESmlrKbJN4vFxsLPP3/vKPv6dqduXb2xQqm8TpNvFjt3Ds6fD0HE9qW+775+lC1rcVBKKcvpQNMsdvIkDB/+GX37vs3+/T/Ru3ePFPdRSuV+mnyzkDFw+rTttadnGVq0eISGDa2NSSmVM2i3QxaKjLRNpnObhwd4eVkXj1Iq59Dkm4VOnUpcrlwZXSRTKQVo8s1Sv/yymeXLX+XYsS3Ex8fp0vBKKQdNvlkkOhrWrPmW1avf4/33W7Fy5RtUqWJ1VEqpnEKTbxY5fdpw4MDPjnKLFp0opEu1KaXsNPlmkT/++Ivz508CULBgUbp0ecDagJRSOYom3ywQHw9r1vzkKNev34WaNXUiHaXUHZp8s0B4OOzbt9pR9vfvTpkyFgaklMpxNPlmgSNHrhEcvNFR7tq1KyIWBqSUynE0+WaBX375jdjYmwBUrOhL48aVLY5IKZXTaPLNZDExsHXrL46yr283KmvuVUolock3k4WEwKFDvzrK99/fRYeYKaXuosk3k+3efYawsMMA5MtXgE6dWlsckVIqJ7Ik+YpIKRFZIyLB9ueSTtpUEZENIhIkIgdF5NkEdZNE5KyI7LE/eiSomyAiR0XksIh0za5zui0mphKvvbaHAQPep1On56ldWy97lVJ3s2pKyfHAOmPMVBEZby+/kqRNLPCiMWaXiBQFdorIGmPMIXv9R8aYDxLuICL1gSGAL1ARWCsiPsaYuCw9G7uLF+H6dTeqVPGjShU/8uWD8uWz452VUq7Gqm6HvsBc++u5QL+kDYwxocaYXfbXV4AgoFIqjrvQGHPDGHMCOAo0z6ygU3L2bOJyhQo6i5lSyjmrkm95Y0wo2JIsUO5ejUWkOtAY2JZg8z9FZJ+IfJWg26IScCZBmxBSTtiZJmnyrZRt76yUcjVZlnxFZK2IHHDy6JvG43gCS4HnjDGX7Zs/A2oB/kAo8N/bzZ0cwiRz3DEiEigigREREWkJySljIDAwiAsXQhzbNPkqpZKTZX2+xphOydWJSJiIeBljQkXECwhPpp0HtsQ73xizLMGxwxK0+RJYZS+GAAknbqwMnEsmvhnADICAgACnCTotzp+HRYteYd++HyhXrjaPPTaL0qXbZPSwSqlcyqpuh5XASPvrkcCKpA1ERIBZQJAx5sMkdQkX4+kPHEhw3CEiUkBEagDewPZMjt2p06djOXLkdwDCw4/i7V1BbylWSiXLquQ7FegsIsFAZ3sZEakoIrenA2sFjAA6OBlS9r6I7BeRfUB74HkAY8xBYDFwCFgNPJ1dIx02bdpNTIytV6REiYo0a+adHW+rlHJRlgw1M8acBzo62X4O6GF/vRnnfbgYY0bc49hTgCmZE2nqxMfDH3/85ij7+LSnUiW97FVKJU/vcMsEkZEQFPSbo9ygQXtK3nXbiFJK3aHJNxOcORPL0aObHOW2bdtZF4xSyiVo8s0Ef/yxh5iYKwCUKFGJJk1qWhyRUiqn0+SbQcbA1q13Jk738WlLxYra36uUujdNvhkUFQVBQXeSb716bbS/VymVIk2+GXT2bDxHj252lB94oI2O71VKpUiTbwadOhVN06aD8PKqT9GiZWnWrK7VISmlXIBVU0rmGpcueTJ8+OcAxMRcxctLL3uVUinTK98MuHwZoqPvlD09PXWJeKVUqmjyzYC//05cLlcO3PQrqpRKBU0VGRAWlrhcoYI1cSilXI8m3wz4+eefmDatFz/99C6nT+/WJYOUUqmmH7il082bsGPHOvbv/5H9+3/kxo0rlC/f2OqwlFIuQq980yk8HI4d2+Io+/m1JH9+CwNSSrkUTb7pdObMDc6c2eUot2zZwsJolFKuRpNvOm3dupvY2JsAlCvnTZ06OsZMKZV6mnzTwRjYsWOro1yz5v36YZtSKk00+abDpUsQHPyno+ztfT/Fi1sYkFLK5WjyTYewMDh5cpuj3LTpfTqZjlIqTTT5psPhw2GcP38KAA+Pgtx3XyOLI1JKuRpNvumwdeud1eirVm1KxYoeFkajlHJFmnzTKDYW9u270+VQvXpzypWzMCCllEvSO9zS6Px56NZtAr6+3ThxYjsNG95PwYJWR6WUcjWafNMoIgIKFChC7doPULv2A9SqZXVESilXpN0OaRQRkbhctqw1cSilXJsm3zTS5KuUygyafNPg5k24eDHxNl25QimVHpYkXxEpJSJrRCTY/nzXYusiUkVENohIkIgcFJFnE9QtEpE99sdJEdlj315dRK4nqPs8M+OOjExcLlkSPHSUmVIqHaz6wG08sM4YM1VExtvLryRpEwu8aIzZJSJFgZ0issYYc8gY89DtRiLyX+BSgv2OGWP8syLoChVg0CBb10NEBBQunBXvopTKC6xKvn2BdvbXc4HfSJJ8jTGhQKj99RURCQIqAYdutxERAQYDHbI8Ymzrs5UqZXvUqZMd76iUyq2s6vMtb0+ut5PsPW9TEJHqQGNgW5Kq1kCYMSY4wbYaIrJbRH4XkdaZGLNSSmWaLLvyFZG1gLMlJSem8TiewFLgOWPM5STVQ4EFCcqhQFVjzHkRaQp8LyK+TvZDRMYAYwCqVq2alpCUUirDsiz5GmM6JVcnImEi4mWMCRURLyA8mXYe2BLvfGPMsiR1+YABQNME73kDuGF/vVNEjgE+QKCT+GYAMwACAgJMGk9PKaUyxKpuh5XASPvrkcCKpA3s/bmzgCBjzIdOjtEJ+MsYE5Jgn7Ii4m5/XRPwBo5ncuxKKZVhViXfqUBnEQkGOtvLiEhFEfnJ3qYVMALokGDoWI8ExxhC4i4HgDbAPhHZCywBnjLGRGXliSilVHqIMfofd0BAgAkMvKtnQimlMkREdhpjApzV6R1uSillAU2+SillAU2+SillAe3zBUQkAjiVBYcuA0Sm2Cpn03PIGfQccoa0nkM1Y4zTuQ81+WYhEQlMrrPdVeg55Ax6DjlDZp6DdjsopZQFNPkqpZQFNPlmrRlWB5AJ9BxyBj2HnCHTzkH7fJVSygJ65auUUhbQ5JsFROQrEQkXkQNWx5Je91rGyVWISEER2S4ie+3nMNnqmNJDRNztc1SvsjqW9LAv9bXfPj+LS97HLyIlRGSJiPxl/51okeFjardD5hORNsBVYJ4xpoHV8aSHfapPr4TLOAH9jDGHUtg1x7DPjFfEGHPVPj3pZuBZY8yfFoeWJiLyAhAAFDPG9LI6nrQSkZNAgDHGZcf4ishcYJMxZqaI5AcKG2MuZuSYeuWbBYwxGwGXnk3NGBNqjNllf30FuL2Mk8swNlftRQ/7w6WuNkSkMtATmGl1LHmViBTDNmPiLABjzM2MJl7Q5KtS4R7LOOV49n/Z92CbsH+NMcbVzuFj4F9AvMVxZIQBfhWRnfYVZFxNTSACmG3v/pkpIkUyelBNvuqeUljGKcczxsTZV7OuDDQXEZfpBhKRXkC4MWan1bFkUCtjTBOgO/C0vVvOleQDmgCfGWMaA9ewrbieIZp8VbLutYyTq7H/m/gb0M3aSNKkFdDH3me6ENvCAt9YG1LaGWPO2Z/DgeVAc2sjSrMQICTBf01LsCXjDNHkq5xKxTJOOZ59WakS9teFsC89ZWlQaWCMmWCMqWyMqY5t5Zb1xpiHLQ4rTUSkiP0DW+z/qncBXGoUkDHmb+CMiNSxb+oIZPiD5yxbQDMvE5EFQDugjIiEAG8aY2ZZG1Wa3V7Gab+9zxTgVWPMT8nvkuN4AXPt6/q5AYuNMS45XMuFlQeW2/6Wkw/41hiz2tqQ0uUZYL59pMNx4NGMHlCHmimllAW020EppSygyVcppSygyVcppSygyVcppSygyVcppSygyVcppSygyVcppSygyVepVBCRZiKyzz5HcBH7/MAuM0+Eynn0JgulUklE3gEKAoWw3ev/nsUhKRemyVepVLLfWroDiAFaGmPiLA5JuTDtdlAq9UoBnkBRbFfASqWbXvkqlUoishLb1I41sC2x9E+LQ1IuTGc1UyoVROQRINYY8619lrQtItLBGLPe6tiUa9IrX6WUsoD2+SqllAU0+SqllAU0+SqllAU0+SqllAU0+SqllAU0+SqllAU0+SqllAU0+SqllAX+Hz6cs1SwLDdGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(5, 5))\n",
    "ax.plot(x0, y0, \"ro\", label=\"Initial condition\", alpha=0.6, markersize=8)\n",
    "ax.plot(x_test, y_test, \"b\", label=\"Exact solution\", alpha=0.4, linewidth=4)\n",
    "ax.plot(sim[:, 0], sim[:, 1], \"k--\", label=\"SINDy model\", linewidth=3)\n",
    "ax.set(xlabel=\"x\", ylabel=\"y\")\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14a5ec90",
   "metadata": {},
   "outputs": [],
   "source": [
    "#  import matplotlib.pyplot as plt\n",
    "\n",
    "#  import numpy as np\n",
    "\n",
    "#  from scipy.integrate import odeint\n",
    "\n",
    "#  from sklearn.preprocessing import PolynomialFeatures\n",
    "\n",
    "\n",
    "#  from deeptime.sindy import SINDy, STLSQ\n",
    "\n",
    "\n",
    "\n",
    "#  # Generate measurements of the Rossler system\n",
    "\n",
    "#  a = 0.1\n",
    "\n",
    "#  b = 0.1\n",
    "\n",
    "#  c = 14\n",
    "\n",
    "\n",
    "\n",
    "#  def rossler(z, t):\n",
    "\n",
    "#      return [-z[1] - z[2], z[0] + a * z[1], b + z[2] * (z[0] - c)]\n",
    "\n",
    "\n",
    "\n",
    "#  dt = 0.01\n",
    "\n",
    "#  t_train = np.arange(0, 150, dt)\n",
    "\n",
    "#  x0_train = [-1, -1, 0]\n",
    "\n",
    "#  x_train = odeint(rossler, x0_train, t_train)\n",
    "\n",
    "#  x_dot_train = np.array([rossler(xi, 1) for xi in x_train])\n",
    "\n",
    "\n",
    "#  # Plot training data\n",
    "\n",
    "#  fig = plt.figure(figsize=(8, 6))\n",
    "\n",
    "#  ax = fig.add_subplot(111, projection=\"3d\")\n",
    "\n",
    "#  ax.plot(x_train[:, 0], x_train[:, 1], x_train[:, 2], color=\"firebrick\", alpha=0.7)\n",
    "\n",
    "#  ax.set(xlabel=\"x\", ylabel=\"y\", zlabel=\"z\", title=\"Training data (Rossler system)\")\n",
    "\n",
    "\n",
    "#  # Instantiate and fit an estimator to the data\n",
    "\n",
    "#  estimator = SINDy(\n",
    "\n",
    "#      library=PolynomialFeatures(degree=3),\n",
    "\n",
    "#      optimizer=STLSQ(threshold=0.05),\n",
    "\n",
    "#  )\n",
    "\n",
    "#  estimator.fit(x_train, y=x_dot_train)\n",
    "\n",
    "\n",
    "#  # Get the underlying ODE model\n",
    "\n",
    "#  model = estimator.fetch_model()\n",
    "\n",
    "#  model.print(lhs=[\"x\", \"y\", \"z\"])\n",
    "\n",
    "\n",
    "#  # Simulate from novel initial conditions\n",
    "\n",
    "#  t_test = t_train\n",
    "\n",
    "#  x0_test = [2, 3, 0]\n",
    "\n",
    "#  x_test = odeint(rossler, x0_test, t_test)\n",
    "\n",
    "\n",
    "#  # Plot test data\n",
    "\n",
    "#  fig = plt.figure(figsize=(8, 6))\n",
    "\n",
    "#  ax = fig.add_subplot(111, projection=\"3d\")\n",
    "\n",
    "#  ax.plot(x_test[:, 0], x_test[:, 1], x_test[:, 2], label=\"True solution\", color=\"firebrick\", alpha=0.7)\n",
    "\n",
    "#  ax.set(xlabel=\"x\", ylabel=\"y\", zlabel=\"z\", title=\"Testing data (Rossler system)\")\n",
    "\n",
    "\n",
    "#  # Simulate data with SINDy model and plot\n",
    "\n",
    "#  x_sim = model.simulate(x0_test, t_test)\n",
    "\n",
    "#  ax.plot(x_sim[:, 0], x_sim[:, 1], x_sim[:, 2], label=\"Model simulation\", color=\"royalblue\", linestyle=\"dashed\")\n",
    "\n",
    "#  ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6b9de26",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}