GitHub_collection_pykan/tutorials/Example_5_special_functions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "134e7f9d",
   "metadata": {},
   "source": [
    "# Example 5: Special functions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2571d531",
   "metadata": {},
   "source": [
    "Let's construct a dataset which contains special functions $f(x,y)={\\rm exp}(J_0(20x)+y^2)$, where $J_0(x)$ is the Bessel function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2075ef56",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "train loss: 1.31e-02 | test loss: 9.92e-02 | reg: 3.78e+00 : 100%|██| 20/20 [00:06<00:00,  3.20it/s]\n"
     ]
    }
   ],
   "source": [
    "from kan import *\n",
    "\n",
    "# create a KAN: 2D inputs, 1D output, and 5 hidden neurons. cubic spline (k=3), 5 grid intervals (grid=5).\n",
    "model = KAN(width=[2,1,1], grid=20, k=3, seed=0)\n",
    "f = lambda x: torch.exp(torch.special.bessel_j0(20*x[:,[0]]) + x[:,[1]]**2)\n",
    "dataset = create_dataset(f, n_var=2)\n",
    "\n",
    "# train the model\n",
    "model.fit(dataset, opt=\"LBFGS\", steps=20);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f30c3ab",
   "metadata": {},
   "source": [
    "Plot trained KAN, the bessel function shows up in the bettom left"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3f95fcdd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAFICAYAAACcDrP3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAxS0lEQVR4nO3de1zUdb4/8NdnBgZmGGC4Kooat0zwUipeEARzAy/ZRSs77Z7Wsj092spzPG3b5dfuduy23VatrFO555xsO1kbVutiupmKcvGGmpcURURFREAYGJj7zPf3B873gJqX/A4zA6/n47GPFr/MzBvlM6/5XL9CkiQJREREClL5ugAiIup9GC5ERKQ4hgsRESmO4UJERIpjuBARkeIYLkREpDiGCxERKY7hQkREimO4EBGR4hguRESkOIYLEREpjuFCRESKY7gQEZHiGC5ERKQ4hgsRESkuyNcFEAUCSZJw9uxZtLe3Q6/XIyYmBkIIX5dF5LfYcyG6BKPRiKVLlyItLQ1xcXFISkpCXFwc0tLSsHTpUhiNRl+XSOSXBO9ESXRx69atw5w5c2A2mwF09l48PL0WnU6HwsJCFBQU+KRGIn/FcCG6iHXr1mHmzJmQJAlut/tHv0+lUkEIgaKiIgYMURcMF6LzGI1GJCYmwmKxXDJYPFQqFbRaLWpra2EwGLxfIFEA4JwL0Xk++ugjmM3mKwoWAHC73TCbzVixYoWXKyMKHOy5EHUhSRLS0tJQXV2Nq2kaQggkJyfjyJEjXEVGBIYLUTdNTU2Ii4u7psfHxMQoWBFRYOKwGFEX7e3t1/R4k8mkUCVEgY3hQtSFXq+/pseHh4crVAlRYGO4EHURExODlJSUq543EUIgJSUF0dHRXqqMKLAwXIi6EELg8ccf/0mPXbBgASfzic7hhD7RebjPhejasedCdB6DwYDCwkIIIaBSXbqJeHbor1q1isFC1AXDhegiCgoKUFRUBK1WCyHEBcNdnj/TarVYs2YN8vPzfVQpkX9iuBD9iIKCAtTW1mLJkiVITk7udi05ORlLlizBqVOnGCxEF8E5F6IrIEkSNm7ciKlTp+K7777DlClTOHlPdAnsuRBdASGEPKdiMBgYLESXwXAhIiLFMVyIiEhxDBciIlIcw4WIiBTHcCEiIsUxXIiISHEMFyIiUhzDhYiIFMdwISIixTFciIhIcQwXIiJSHMOFiIgUx3AhIiLFMVyIiEhxDBciIlIcw4WIiBTHcCG6DIfDgVOnTuHgwYMAgKNHj6K5uRlut9vHlRH5L97mmOhHGI1GFBYW4pNPPsGBAwdgMplgt9sRGhqKuLg45OTkYP78+Zg0aRKCgoJ8XS6RX2G4EF1EeXk5Fi5ciL179yIzMxMzZ87EyJEjodfrYTQaUVFRgdWrV6Oqqgpz587Fiy++iLi4OF+XTeQ3GC5E5/nHP/6BefPmQa/X45VXXsGMGTNgt9uxcuVK2Gw2RERE4N5774XD4cDKlSvx/PPPIyMjAx9//DH69evn6/KJ/ALDhaiLw4cPY9q0aQgLC8PKlSuRnp4OIQSqq6sxevRotLa2IikpCRUVFYiKioIkSSgpKcF9992HvLw8LF++HCEhIb7+MYh8jhP6ROe4XC68/PLLaGlpwTvvvCMHy6UIIZCdnY3XXnsNX3/9NdauXdtD1RL5N4YL0TlVVVVYvXo1Zs+ejezs7MsGi4cQAnfccQcmTJiADz/8EE6n08uVEvk/LnEhOqesrAzt7e2YM2cOampq0NHRIV+rra2Fy+UCANjtdhw4cAARERHy9QEDBmD27Nl4/vnnUV9fj8TExB6vn8ifMFyIzjl06BB0Oh2Sk5Px8MMPo7S0VL4mSRJsNhsAoK6uDrfccot8TQiBN998EyNGjIDZbEZdXR3Dhfo8hgvRORaLBUFBQQgJCYHNZoPVar3o90mSdME1p9MJrVbbLYSI+jKGC9E58fHxsFgsMBqNGD9+PMLCwuRrFosFZWVlcohkZWXJGyeFEBg8eDAaGhqgUqkQFRXlqx+ByG8wXIjOGTNmDBwOB7Zv345XX32127Xq6mpkZmaitbUV/fr1w2effQaDwSBfF0Lg2WefRf/+/TkkRgSuFiOSjRs3DsnJyfjoo4/Q0dEBtVrd7X8eQgioVCr5z1UqFU6fPo0vvvgCM2fORGRkpA9/CiL/wHAhOicmJgaPPfYYdu3ahbfeeuuKlxTbbDa88MILsFgsePjhh694CTNRb8ZhMaIu5s2bh82bN+PVV1+FTqfDI488gtDQUABAUFAQgoKC5F6MJEkwmUx46aWXsHLlSixevBhDhw71ZflEfoPHvxCdp7GxEY8++ij+/ve/o6CgAAsXLsSwYcNQWVkJt9sNjUaD1NRUbN++HW+88Qb27NmDRYsW4ZFHHuk2fEbUlzFciC6io6MDH374Id566y2cOXMGycnJSEtLQ3h4OFpaWlBZWYm6ujqMGTMGf/jDH5CbmwuViqPMRB4MF6JLqK+vx3fffYfi4mJ8//332L59O3JycjBp0iTk5+dj/Pjx0Ol0vi6TyO8wXIiu0I4dOzBu3Djs2LEDY8eO9XU5RH6N/XiiK6RWq+VlyER0aWwlRESkOIYLEREpjuFCRESKY7gQEZHiGC5ERKQ4hgsRESmO4UJERIpjuBARkeIYLkREpDiGCxERKY7hQkREimO4EBGR4hguRESkOIYLEREpjvdzIbpCkiTB7XZDpVJBCOHrcoj8GnsuRFeB93IhujJBvi6ASCmSJOHIkSM4e/asr0u5JiqVCsOHD0dYWJivSyH6yTgsRr2G2+3Go48+ikGDBkGv1wMAXC4X1Gq1jyu7Olu2bMHvfvc7jBw50telEP1k7LlQrxISEoL58+cjPj4e33zzDZYtW4ZHHnkEM2bMCIghLUmS0N7eDn7mo0Dn/62N6Cew2Wx44403sHbtWjzwwAMoLCzkGzZRD2K4UK9UV1eHffv2AQCam5vx1FNP4fTp0z6uiqjvYLhQr3Ty5EmYTCaEhoZCr9fj5MmT+Pbbb9l7IeohDBfqlaqqquBwODB06FDceeedcLvdWLFiBaxWq69LI+oTGC7U60iShMrKSgBASkoKHnroIeh0Omzbtg07duxg74WoBzBcqNeRJAknTpwAACQlJWHs2LHIysqCxWLBihUr4Ha7fVwhUe/HcKFex+VyyRsp+/Xrh5CQEMybNw9qtRpFRUU4duyYjysk6v0YLtTrOBwONDc3AwDi4+MhhEB+fj7S0tLQ0NCAL774gkNjRF7GcKFex2q1orW1FSqVCrGxsQCA6OhozJ07FwDw6aefwmg0+rBCot6P4UK9TkdHB0wmE4KDgxETEwMAEELgnnvuQWxsLA4dOoT169ez90LkRQwX6nXa2tpgtVqh0WhgMBjkP09NTUVBQQGcTifeeecdmEwm3xVJ1MsxXKjXMZlMcDqd0Ol03U4WVqvVeOyxx2AwGFBeXo6VK1ey90LkJQwX6nWuv/56fPbZZ3j77bflYTGgc2hszJgxuP/+++FyubB48WLU1dX5sFKi3ovhQr1OdHQ0br31VsyZMwehoaHdrqlUKixYsABDhgzBkSNH8P7773PfC5EXMFyoTxFC4LrrrsPjjz8OIQSWL1+OQ4cOcXiMSGEMF+pzhBC4//77MWLECJw5cwZ/+tOf4HK5fF0WUa/CcKE+KTo6Gr/5zW8QHByML774AuXl5ey9ECmI4UJ9khACt99+O3Jzc2EymfDyyy/DbDb7uiyiXoPhQn2WTqfDM888g/DwcGzYsAGfffYZey9ECmG4UJ8lhEB2djZ+8YtfwOl04qWXXsLRo0cZMEQKYLhQn6ZWq/HUU0/hhhtuQE1NDX7/+9/DZrP5uiyigMdwoT5NCIHExEQsWrQIWq0WX375JT755BP2XoiuEcOF+jwhBGbNmoX7778fdrsdzz//PPbt28eAIboGDBciAMHBwXjuuedw4403oq6uDk8++SRaW1t9XRZRwGK4EKGz95KQkIDXX38dBoMBGzZswB//+Ec4HA5fl0YUkBguROcIIZCbm4tnnnkGKpUK77zzDj7//HMOjxH9BAwXoi7UajV+/etf46677oLFYsGTTz6JkpISBgzRVWK4EJ1Hq9Xi9ddfx4QJE3DmzBk8/PDD+OGHHxgwRFeB4UJ0Hs/8y/vvv4+0tDRUVlbigQceQHV1NQOG6AoxXIguQgiBjIwMfPjhhxg4cCAqKiowb9481NTUMGCIrgDDhehHeI6H+eCDD9CvXz+UlZXhn//5n9mDIboCDBeiSxBCoKCgAMuXL0e/fv1QXl6Of/qnf8KBAwcYMESXwHAhugwhBKZPn47//u//lofI7r77bmzevJm3SCb6EQwXoisghEB+fj4+/fRTXH/99Th8+DDmzp2L//mf/4Hdbmcvhug8DBeiKySEQFZWFgoLC5GTk4OmpiY89thjeOKJJ9DU1MSAIeqC4UJ0FYQQGDZsGD7//HM8+OCDkCQJ7733Hu644w6Ul5dzmIzoHIYL0VUSQiAuLg5vv/02li5diri4OGzduhW33XYbXnzxRfZiiMBwIfpJhBAICQnBQw89hNWrVyMvLw+tra1YtGgRpk+fjq+++gpWq5UhQ30Ww4XoGqhUKowZMwarVq3Cyy+/jPj4eOzatQs///nPce+996K8vBwOh4MhQ30Ow4XoGgkhEBkZiSeeeALr16/HL37xC6jVaqxevRrTp0/HvHnzUFpaCpvNxpChPoPhQqQQlUqFYcOGYfny5fj6669xyy23wG63Y+XKlZg+fTpuv/12fPbZZ/KcDIOGejOGC5GChBAIDg7GlClTsGrVKnz22WfIz8+HJEn49ttvcf/992PSpEn4zW9+g5KSErS1tTFoqFdiuBB5gRACOp0Os2bNwpdffol169Zh/vz5iI+Px9GjR7FkyRIUFBQgOzsb//7v/47169fjzJkzsNlsvi6dSBFBvi6AqDcTQiA0NBRZWVmYOHEiamtrsXbtWhQWFmLnzp344YcfcODAAbz33nuIi4vD1KlTkZyc7Ouyia4Zw4WoBwghIITA4MGD8atf/Qrz5s3DsWPHsGHDBqxZswY7d+5EQ0MDLBYL1Gq1r8slumYMF+pVJElCS0sLgoODfV3KZcXFxWHu3LmYM2cOTp8+jerqauj1emzZssXXpRFdM4YL9RpCCAwZMgRvv/12QH/6t1gsiIyM9HUZRNdESFymQr1Eb1p15RlGIwpUDBciIlIclyITEZHiOOdCdIW6dvI5ZEV0aey5EF2h3bt3Q61WY/fu3b4uhcjvMVyIiEhxDBciIlIcw4WIiBTHcCEiIsUxXIiISHEMFyIiUhzDhYiIFMdwISIixTFciIhIcQwXIiJSHMOFiIgUx3AhIiLFMVyIiEhxDBciIlIcw4WIiBTHcCG6ApIkoaWlBQDQ0tIC3h2c6NIYLkSXYDQasXTpUqSlpeFnP/sZJEnCz372M6SlpWHp0qUwGo2+LpHILwmJH8GILmrdunWYM2cOzGYzgIvf5lin06GwsBAFBQU+qZHIXzFciC5i3bp1mDlzJiRJgtvt/tHvU6lUEEKgqKiIAUPUBcOF6DxGoxGJiYmwWCyXDBYPlUoFrVaL2tpaGAwG7xdIFAA450J0no8++ghms/mKggUA3G43zGYzVqxY4eXKiAIHey5EXUiShLS0NFRXV1/VijAhBJKTk3HkyBF5PoaoL2O4EHXR1NSEuLi4a3p8TEyMghURBSYOixF10d7efk2PN5lMClVCFNgYLkRd6PX6a3p8eHi4QpUQBTaGC1EXMTExSElJuep5EyEEUlJSEB0d7aXKiAILw4WoCyEEHn/88Z/02AULFnAyn+gcTugTnYf7XIiuHXsuROcxGAwoLCyEEAIq1aWbiGeH/qpVqxgsRF0wXIguoqCgAEVFRdBqtRBCXDDc5fkzrVaLNWvWID8/30eVEvknhgvRjygoKEBtbS2WLFmC5OTkbteSk5OxZMkSnDp1isFCdBGccyG6ApIkYePGjZg6dSq+++47TJkyhZP3RJfAngvRFRBCyHMqBoOBwUJ0GQwXIiJSHMOFiIgUx3AhIiLFMVyIiEhxDBciIlIcw4WIiBTHcCEiIsUxXIiISHEMFyIiUhzDhYiIFMdwISIixTFciIhIcQwXIiJSHMOFiIgUx3AhIiLFMVyIiEhxDBeiy3A4HDh16hQOHjwIADh69Ciam5vhdrt9XBmR/+Jtjol+hNFoRGFhIT755BMcOHAAJpMJdrsdoaGhiIuLQ05ODubPn49JkyYhKCjI1+US+RWGC9FFlJeXY+HChdi7dy8yMzMxc+ZMjBw5Enq9HkajERUVFVi9ejWqqqowd+5cvPjii4iLi/N12UR+g+FCdJ5//OMfmDdvHvR6PV555RXMmDEDdrsdK1euhM1mQ0REBO699144HA6sXLkSzz//PDIyMvDxxx+jX79+vi6fyC8wXIi6OHz4MKZNm4awsDCsXLkS6enpEEKguroao0ePRmtrK5KSklBRUYGoqChIkoSSkhLcd999yMvLw/LlyxESEuLrH4PI5zihT3SOy+XCyy+/jJaWFrzzzjtysFyKEALZ2dl47bXX8PXXX2Pt2rU9VC2Rf2O4EJ1TVVWF1atXY/bs2cjOzr5ssHgIIXDHHXdgwoQJ+PDDD+F0Or1cKZH/4xIXonPKysrQ3t6OOXPmoKamBh0dHfK12tpauFwuAIDdbseBAwcQEREhXx8wYABmz56N559/HvX19UhMTOzx+on8CcOF6JxDhw5Bp9MhOTkZDz/8MEpLS+VrkiTBZrMBAOrq6nDLLbfI14QQePPNNzFixAiYzWbU1dUxXKjPY7gQnWOxWBAUFISQkBDYbDZYrdaLfp8kSRdcczqd0Gq13UKIqC9juFCfV1NTg02bNqGkpARmsxlGoxHjx49HWFiY/D0WiwVlZWVyiGRlZckbJ4UQGDx4MBoaGuB0OlFVVYXMzEyEhob66kci8jkuRaY+5+TJkyguLkZxcTE2bdqEEydOQAiBIUOG4OjRo1i2bBkeeuihbo+prq5GZmYmWltbcd1112Hnzp0wGAzydSEEnn32WbzxxhtQq9UIDQ3F+PHjMXnyZOTl5SEzM5NLlKlPYbhQr1dXV4dNmzZh8+bN2LRpE44dOwYAGDlypPzmn52dDbfbjezsbERFRWHt2rXdJux/bJ8L0DlMVldXh9zcXMyaNQvz5s3D5s2bsXnzZmzZsgVGoxFarRYTJkyQX2/MmDHQaDQ++fsg6gkMF+p1zpw50y1MqqqqAAAZGRnym3tOTg5iYmIueOyyZcvwxBNP4LnnnsPTTz8tD31dKlysViv+7d/+DatXr8aGDRswdOhQ+flcLhf27t2LzZs3o7i4GCUlJWhra4NOp8PEiRORm5uLvLw83HTTTQgODu6Bvx2insFwoYDX2NiI4uJiOUwqKysBADfccANyc3ORm5uLyZMnX9HZXx0dHXjwwQexZs0a/Md//AceeeQRhIaG4tixYxg3bpw8LLZ9+3YYDAaYTCa89NJLeP/997F48WI88MADl3x+p9OJPXv2yPWWlpaivb0der0eWVlZcr033ngjD8OkgMZwoYDT1NSELVu2yPMmP/zwAwAgLS1NfnPOzc39yed8NTY24tFHH8Xf//53FBQUYOHChRg2bBgqKyvhdruh0WiQmpqK7du344033sCePXuwaNEiPPLII1Cr1Vf1Wg6HA7t375Z7WmVlZTCbzYiIiMCkSZPkn2XkyJFX/dxEvsRwIb/X0tKCLVu2yG/A+/btAwAkJyd3C5MBAwYo9podHR348MMP8dZbb+HMmTNITk5GWloawsPD0dLSgsrKStTV1WHMmDH4wx/+gNzcXKhU137ghd1uR0VFhRyc5eXlsFqtMBgMyM7Oln/W4cOHK/J6RN7CcCG/09raipKSEjlMvv/+e0iShCFDhiAvL0+eN+mJjYr19fX47rvvUFxcjOrqalitVkRFRWH48OHIz8/H+PHjodPpvPb6NpsNO3bskMNm27ZtsNlsiI6ORk5Ojhw2V3IOGlFPYriQz5lMJpSWlspvoLt374bb7cbAgQORl5cnT3oPGTLEp3W6XC5IkgSVSuWzXoPVasW2bdvk+aXt27fD4XAgNjYWkydPlsNm6NChDBvyKYYL9bj29naUl5fLPZOKigq4XC4kJCR065kkJSXxDfIyzGYztm7dKofNzp074XQ6ER8f323IMDU1lX+X1KMYLuR1njdAT5js2LFDfgPsGiZ8A7x27e3t2Lp1q7xBdNeuXXJwdw0bBjd5G8OFFOcZuvGEybZt2+BwOBAXF4fJkyfLYcKhG+9ra2tDeXn5BUOOiYmJ8r9Dbm6uz4ccqfdhuNA1s9ls2L59+0UnnbvOA3DS2feMRiNKS0vlTZ1dF0t0DRue6kzXiuFCV81ut2Pnzp3y0MvWrVvl5bJdVzBxuaz/a25u7raYwrPMOykpSf53zMvLQ0JCgo8rpUDDcKHLcjgc2LVrlxwm5eXl8ka/7Oxsed6EG/0C39mzZy+6QTU1NVWRDarUdzBc6AJdjyjZtGmTfIdGvV7fbdc4jyjp/RobG+UhtOLiYvlonaFDh3YLm9jYWB9XSv6G4ULy4YqeMCktLZUPV8zKypJ7JqNHj+bhin1cfX19t7DxHAqanp7e7Ry36OhoH1dKvsZw6YPcbjf2798vh0lJSQmMRiNCQ0MxceJEOUzGjh3LY+Hpkurq6uSgKS4ulm9nMGLECDlsPLcxoL6F4dIHSJKEH374QQ6TLVu2oLm5GSEhIRg/frwcJuPGjeMNreianDx5Ut7QWVxcLN+IbdSoUXLYTJo0CZGRkb4ulbyM4dILSZKEyspKuYFv2bIFjY2NCA4Oxvjx4+VGPn78eN6Kl7yqpqamW9icOnUKKpUKN910k/x7mJWVhfDwcF+XSgpjuPQCkiShqqpKbsDFxcVoaGhAUFAQMjMz5UY8YcIErx6ySHQpkiTh2LFj3W4xXV9fD7VajdGjR8t7bCZOnIiwsDBfl0vXiOESgDyNtGuYnD59Gmq1GmPGjJH3JkyYMAF6vd7X5RJdlOdDUdc5G8+HorFjx8ph4+2Tp8k7GC4B4vjx493CpLa2ttvwQl5eHiZOnNjtvu9EgcQznNs1bM6ePQuNRoPMzEx5bpDDuYGB4RIgRowYgSNHjsgTo3l5ecjKyoLBYPB1aURe4Xa7cfDgQTloPAtRPv74Y9x9992+Lo8ug+ESINxuN4QQPJuL+ixJkiBJEttBgGC4EBGR4nh2h0I8k5Nnz571dSnXRKVSISMjg6t16KqxDVBXDBeFSJKEpUuXIjExUbEVWp5O5aWGAFwuF1QqlWLDBCUlJfh//+//YcSIEYo8H/Ud3mgD5/PcalqtVnttaIxtQBkMFwWFhITggQceUOTE2DNnzuC5555DQkICnn766QuWYkqShF27duGFF17ArFmz8MADD1zz8faSJKG9vR0cKaWfSsk2AHT+TtpsNmzbtg1FRUXYt28fHA4HBg8ejPz8fNxyyy2Ijo5WLGjYBpTDcPFDkiShsLAQf/nLXxAcHIzc3FxMnTq12/e43W68+eabWLt2LQ4cOIAZM2bwnhvUq7jdbuzZsweLFi3Cpk2bYLVau11fuXIl0tPT8fTTT+O2225DcHAwJ/r9CO/k5IckScK2bdsgSRLsdju2b99+wSeptrY2fP/99wA6T6r1nE5L1Bu4XC789a9/xezZs7F27Vqo1WoUFBTglVdewVtvvYVf/vKXiI+Px/79+zF//nw8+eSTaGlpYY/Dj7Dn4ofsdjuOHj0qf3348GF5CaZHU1MTmpqaAHTezKu6uhrZ2dn85EYBz+1245NPPsETTzwBk8mEzMxMvPjii5g4caJ8y4eHHnoIR48exauvvorPP/8cH3zwAU6ePIn33nsP8fHxbAd+gD0XH5IkCUajEW1tbd0+cbW3t6OhoUH++uTJk3A6nd0eW19fD7PZLH9dXV3d7brL5UJTUxPsdjs/zVHAkCQJ69atw29/+1u0t7fj1ltvxeeff47JkydDo9HIe1xUKhVSU1OxbNkyvPbaa4iIiMCaNWuwYMGCC9oT+QbDxUckScKhQ4eQn5+P2267DbW1tfI1o9GI1tZW+euGhoYLxptPnz7dLXBOnjzZ7bn/8pe/ICsrC7/73e8uCCYif+RpEwsXLoTRaMTUqVPx7rvvon///hftiQghEBISgl/96ldYvHgxwsPDsXr1arz55ptwuVw++AmoK4aLl0mSBKvVCqvVesGnqRUrVmDv3r3YunUrVq9eLV9vbGyExWKRv6+1tRXt7e3dnvP06dPdnq+urk4OkY6ODrz77rs4ceIEPvroowvmY1wuF1paWuBwOPgJj/xGe3s7nn76adTU1GDYsGF46623EBcXd9khLrVajblz5+LJJ5+EEALvvvsuNmzYwN9tH2O4eJHnIL5Zs2bhrrvu6ta7sFqtKC0tlb8uKyuTG0NDQwOcTifUajVUKhU6OjrQ1tbW7bnr6+sBQL6HvWcIDOjsxXiGyVpbW7Fjxw75ud1uN/785z9j4sSJeOqpp2Cz2bz00xNdOc/v5fr16xEREYHXX38dSUlJVzx3olar8etf/xoFBQVob2/H73//ezQ2Nnq5aroUhosXud1u/OlPf8KWLVuwfv16rFixQn6Tb2pqkm8JCwCVlZVyb6WhoQFutxsDBgxAeHg47HY7Wlpauj23p+EkJydDrVajtbVVnoOpqqpCR0cHgM6Aq6iokB/X1NSExYsXo6amBv/1X/8lr0oj8hXPcNiSJUvgdrvx8MMPY8qUKVc9KR8WFoZFixYhISEB33//Pd5++20Oj/kQw8WLGhsbsXHjRvnrjRs3yr2LEydOdJtXOXPmjLyU0tMrGTJkCAwGA5xOJ5qbm+Xvdbvd8hEbw4YNg0ajgclkkjd/HT58GG63W/7+gwcPykNmBw4ckOd3rFYrioqKvPTTE10Zh8OB119/HfX19Rg1ahQWLFgAtVp91c8jhEB6ejoWLlwIlUqF5cuXo6Kigh+efITh4kUHDhzA6dOn5a+rqqrQ2Ngon8Fkt9sRFRUFrVaLtrY2OVQ8jxk0aBCioqK6hQkAOJ1OuSeTlpYGrVYLq9Uq/9mRI0cAAP3794dKpcKJEydgMpnkXf0Oh0P+VFhWVnbBYgGiniJJEkpLS/H1118jJCQETz31FGJjY3/y86lUKvzyl7/ExIkT0dLSgtdee41Dvz7CcPESSZKwY8cOOJ1OJCUlITIyEi0tLTh+/DiAzr0rADBy5EgkJCTAbrejtrYWbrdbXoacmJgo36/l7Nmz8icwu90Oo9EIIQSuu+466PV6OBwONDU1wel0yq+Rl5cHrVaLs2fPyqG2f/9+AJ33h9FoNKiurpZDjain2Ww2LF26FGazGTfffDOmTZt2zXtUIiIi5COTvv32W07u+wjDxUvcbjd2794NALj55psxZMgQ2O12ecjKEy6jRo3CwIED4Xa7UVNTA6fTKc+nDBw4ENHR0QAgb5gEOoez2tvboVarMWDAAERFRcHlcuHMmTOwWq2or6+HEALjxo1DdHQ0zGYzjh8/DpvNJr/urFmzEB8fj9bWVhw6dKgn/2qIAPxfr2XTpk3QarX413/9V0XuMCmEwOTJkzF9+nTYbDa8/fbb3VZfUs9guHhJe3s7Dh48CCEEsrKykJKSAqBz/sNqtaKmpgYAkJGRgcGDBwMAampqYLPZ0NzcDCEE+vXrh6ioKADoNudiNpthsVgQFBSE2NhYxMTEAOict2lra8PZs2cRHByM9PR0DBw4EC6XC1VVVTAajairq0NQUBAmTJiAlJQUuFwu7N69m5/sqMc5HA588MEHsFqtmDJlCiZOnKjYzvrg4GA8+uijCAsLQ1lZGUpLS/k73sMYLl5SV1eH06dPQ6fTYcSIERg2bBiAzlVhTU1NOHPmDDQaDVJSUnDdddcBAI4fP462tja0tbUhKCgI8fHxcs+l67lJZrMZNpsNwcHBiIyMlE+gra+vR2NjI0wmE3Q6HQYNGoTk5GT5devq6tDS0gK9Xo+UlBSMGjUKALBnz55uCwCIvE2SJOzduxcbNmyARqPBv/zLv0Cj0Sj2/EIIjB07FlOmTIHNZsPy5cu5mbiHMVy8QJIk/PDDD+jo6EC/fv2QmJiIoUOHQgiBmpoaHDlyBK2trQgPD0diYiKGDBkCoDOQGhsb0dHRgZCQEMTExMjhYjQa5WWVJpMJdrsdoaGhCAsL6xYup06dgs1mQ2RkJKKionD99dcD6Jzkr6yshM1mQ3x8PGJjY3HTTTdBCIFDhw5126RJ5G2SJGHFihUwmUwYO3YscnJyFD8PLDg4GPPnz4dGo8HGjRtx4MAB9l56EMPFSyoqKuB2u3HDDTcgIiICqampCAkJQUNDA0pLS2G325GQkIDY2FgMHjwYQUFBaGpqkofG9Hq9HBBAZ6B4wqWtrQ0ulwtarRZarVY+ar+hoQHV1dVwu93o378/wsLCMGzYMKhUKhw7dkw+Xfm6666DTqdDRkYGdDodTp8+jbq6Op/9XVHfc+LECaxevVpe3XX+/YqUIIRATk4ObrzxRrS1teHTTz9luPQghosXmM1mbNu2DQAwduxYqNVqeXK+ra0NX375JSRJwtChQ6HT6dC/f39otVq0trZi//79cDqdMBgM0Ov1MBgMEELIvRWgc4jM7XYjLCwMISEhSEhIgBACjY2NOHjwIIDOZczBwcG4/vrrERYWhvr6enz77bcAgPT0dKjVagwaNAj9+/dHe3u7vPiAyNskScKqVatQX1+P5ORkTJ8+3WunGIeFheHnP/85hBD46quvuDKyBzFcFFZRUYG77roL27dvh0ajwaRJkyCEQFRUFJKSkuByueTVWWPGjIEQAjExMTAYDLBYLNi6dSskSUJcXBxCQ0MRGRkJtVotz7MAnUNkkiQhPDwcwcHBSEhIkHs+e/fuBdC5c18IgQEDBmDgwIGwWCw4evQohBAYOXIkgM4lm6NHj4YkSVi0aBHefPNNOcCIvKW1tVXuRdxzzz2Ii4vz2msJITBz5kwkJibi5MmTWLt2LXsvPYThojAhBHbu3AmHw4HU1FT5jVyj0WDs2LHy94WGhsqrY8LDw5GQkACXy4Xy8nIAkIfKwsPDoVarYbFYYLFY5GP6AcjXEhISoNPp0NzcLO9jSU1Nlb8nMzNTfl29Xo9Ro0bJx5bPmDEDarUaNTU1KC0t5X0wyKskScLmzZtx8OBBREdH45577vH679yAAQMwc+ZMuN1ufPrpp9w03EMYLgobMWIEHnvsMYwaNQq/+93v5E2QADB9+nRotVoAnftbhg8fDqAzeNLS0gBAPh/Ms3RZr9cjODgYdrtdXqvvCZfIyEioVCpER0cjOjoaLpcLVqtVfj7PvS/uvPNOeSXOjTfeKK8gE0JgxowZuPvuuzFq1Cg88cQT8s2YiLzB6XTif//3f+F0OnHzzTfLH4K8SQiBuXPnQqvVoqKiAnv37mXvpQcwXBQWHByM5557Dhs3bsQdd9whfyrz7Hd55plnMG3aNLzyyivQ6/XytRtvvFF+DpVKhYyMDAghoNPpEBISAofDIa/o6hounp5PUlKS/PjY2Fh5ebMQAjfffDOeeuopTJs2DS+88EK3jWoRERH44IMP8N1332HSpEle/Jsh6ryp3ebNmxEcHIz77rvvJ50hdrU87Wv06NEwm80oLCz0+msSb3PsFWq1+qKrXzQaDZ588km4XC6o1epuwTNhwgTodDqYzWZERUVhxIgRADqHz0JDQ+V7ukiSJB946ekVeTZFbtiwAUBn76nr+UwhISF49tlnL3hdz2trNBpoNBp+miOvkiQJRUVFaG5uRnp6OrKysnpsGDY0NBRz5sxBWVkZ1qxZg9/+9rfXdIYZXR57Lj1MCIGgoKALGtXw4cMxffp0hIaG4s4778SgQYMAdAaDVquF2+1Ga2ur/F/g/8JFCIHZs2dj0KBBiIyMxIMPPijf58Vz/cdel6inWCwW/O1vfwMA3HrrrYiMjOyx1xZCYNq0aYiNjUVNTY28cIa8h+HiJ0JDQ7Fs2TKsX78ef/zjH+Vw0Gg00Ov1cLvd8kZKk8kEAIiKipLDIiMjA9988w3Wr1+PW2+9lSFCfufQoUPYt28fdDqdT35HBw8ejKysLDidTnz11Vc8lcLLGC5+QggBg8GAsWPHQq/Xyw3Ps2IM6FzC6XQ6YTKZIITo9slPCIHU1FSMGDGiR8axia6GJElYt24dOjo6MGzYMGRkZPR4DWq1GrfffjtUKhWKi4t5p0ovY7j4OZVKhYiICACdE/l2ux0dHR0XhAuRP7NarVi3bh0AYNq0aV7ZkX85nh37/fv3R11dHcrLyzk05kUMFz+nUqnkEDEajbBarbBarVCr1XKPhsjfHTlyBPv374dWq0V+fr7Phm0HDBiArKwsuFwuFBUVcWjMixguAcATLq2trbBYLLBard2Gy4j8mSRJ2LBhA9rb23H99dcjPT3dZ7V4Ng6rVCps3ry5232SSFkMlwDgObzSaDTKZ4xpNBqEhYX5uDKiy7Pb7Vi7di0AYOrUqfL+Ll8QQiA7Oxvx8fGoq6vDjh07ODTmJQyXAOBZctza2ipP6oeGhsq7/Yn82fHjx/H9999Do9GgoKDA1+VgwIAByMzMhNPpxDfffMNw8RKGi5/zrCIDOu9u2dLSApfLJW+uJPJnnrPEjEYjkpKSMHLkSJ8vk1er1fJJzJ7aSHkMlwBgMBigUqlgMpnQ2NgIt9sNvV6PkJAQX5dGdEkul0vuHeTm5vrFCkfPqjGDwYDjx49j3759vi6pV2K4BADPAZUWiwWnT5+GJEmIiIjotgufyB/V1dWhoqICQUFBKCgo8HmvxWPw4MEYOXIk7HY71q9fz6ExL2C4BABPkFgsFpw6darbnxH5K0mSsH37djQ0NCAhIUG+f5E/0Gg0mDp1KgBg48aN8onjpByGSwAIDw+HRqOB1WrFiRMnAADR0dFQqfjPR/5LkiR8++23cLvdGDdunFdvCna1hBCYMmUKtFotKisrUV1d7euSeh2+OwUAnU4HrVYLh8OBmpoaAP+3PJnIX7W2tqKsrAxCCNxyyy1+92Fo6NChSE5OhslkQklJCYfGFOZf/9p0UVqtFjqdDg6HAydPngTQ2XMh8meHDh3CiRMnEBkZKd911Z+Eh4fL9zDasGEDXC6XjyvqXRguASA0NFReZeNwOADAr4YYiM4nSRK2bNkCm82GG264AYMHD/Z1SRc1ZcoUqFQq7N69G2fPnvV1Ob0KwyUABAUFdbtdshACsbGxfvdJkMjD6XRi8+bNAIDs7Gy/3JMlhMDo0aMRExOD+vp67N+/39cl9SoMlwCgVqu79VTUajXvokd+zfNmHRQUhMmTJ/u6nB+VkJCAjIwMOBwObNmyxdfl9CoMlwAghEBCQoL8dUhICMOF/Nq+ffvQ2NiI+Ph4DB8+3G972cHBwcjJyQEAlJSUwG63+7ii3oPhEiAGDhwo/3+9Xt9tmIzIn0iShJKSErhcLgwfPtyv5wc9B1lqNBocPHgQdXV1vi6p12C4BAAhBIYMGSIv5YyJifGLYzSILsZms2Hr1q0AgEmTJvn9Zt9hw4ahf//+aGlpwe7du31dTq/BcAkQKSkp8inIqampPBGZ/FZ9fT0OHz4MjUbjl0uQzxcTE4ObbroJbrdbXoRA147hEiCSk5Mxfvx4hIWF4bbbbvO7DWlEHtXV1bDb7ejfvz9uuOEGX5dzWSqVCjk5OVCr1Th06BDnXRTi3/3VACNJEoxGI4KDg73y/IsXL0ZdXR3S09PR0tLildew2WxeeV7qGyRJQnp6Ov7617+ioaEBQgg0Nzf7uqzLysnJwX/+539i5MiRWLVqla/L6RUYLgoRQmDw4MFYtmwZ1Gq1V1/rb3/7m9ee22KxICIiwmvPT72Xpw38+c9/lttAWVmZj6u6Ort27WIbUIiQeKCOIiRJ6jVnEwkh/H6cnPwP2wB1xXAhIiLFcVaYiIgUx3AJEJIkwe1295phB6Kfgu0gcDBcAsSePXug1WqxZ88eX5dC5DN79uyBTqdjOwgADBciIlIcw4WIiBTHcCEiIsUxXIiISHEMFyIiUhzDhYiIFMdwISIixTFciIhIcQwXIiJSHMOFiIgUx3AhIiLFMVyIiEhxDBciIlIcw4WIiBTHcCEiIsUxXAKAJEloaWkBALS0tPBGSdQnedpB1/+S/2K4+DGj0YilS5ciLS0NU6dOhd1ux9SpU5GWloalS5fCaDT6ukQir2M7CExCYvz7pXXr1mHOnDkwm80A0O1TmhACAKDT6VBYWIiCggKf1EjkbWwHgYvh4ofWrVuHmTNnyvcL/zEqlQpCCBQVFbFhUa/DdhDYGC5+xmg0IjExERaL5ZINykOlUkGr1aK2thYGg8H7BRL1ALaDwMc5Fz/z0UcfwWw2X1GDAgC32w2z2YwVK1Z4uTKinsN2EPjYc/EjkiQhLS0N1dXVV7USRgiB5ORkHDlyRB6HJgpUbAe9A8PFjzQ1NSEuLu6aHh8TE6NgRUQ9j+2gd+CwmB9pb2+/psebTCaFKiHyHbaD3oHh4kf0ev01PT48PFyhSoh8h+2gd2C4+JGYmBikpKRc9XixEAIpKSmIjo72UmVEPYftoHdguPgRIQQef/zxn/TYBQsWcBKTegW2g96BE/p+huv7idgOegP2XPyMwWBAYWEhhBBQqS79z+PZmbxq1So2KOpV2A4CH8PFDxUUFKCoqAharRZCiAu6+Z4/02q1WLNmDfLz831UKZH3sB0ENoaLnyooKEBtbS2WLFmC5OTkbteSk5OxZMkSnDp1ig2KejW2g8DFOZcAIEkSmpubYTKZEB4ejujoaE5aUp/DdhBYGC5ERKQ4DosREZHiGC5ERKQ4hgsRESmO4UJERIpjuBARkeIYLkREpDiGCxERKY7hQkREimO4EBGR4hguRESkOIYLEREpjuFCRESKY7gQEZHiGC5ERKS4/w98MZZEafTE+gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 500x400 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "733a2a41",
   "metadata": {},
   "source": [
    "suggest_symbolic does not return anything that matches with it, since Bessel function isn't included in the default SYMBOLIC_LIB. We want to add Bessel to it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "031db28f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  function  fitting r2   r2 loss  complexity  complexity loss  total loss\n",
      "0        0    0.000000  0.000014           0                0    0.000003\n",
      "1        x    0.001598 -0.002293           1                1    0.799541\n",
      "2      cos    0.159266 -0.250262           2                2    1.549948\n",
      "3      sin    0.159266 -0.250262           2                2    1.549948\n",
      "4    1/x^2    0.098715 -0.149929           2                2    1.570014\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "('0',\n",
       " (<function kan.utils.<lambda>(x)>,\n",
       "  <function kan.utils.<lambda>(x)>,\n",
       "  0,\n",
       "  <function kan.utils.<lambda>(x, y_th)>),\n",
       " 0.0,\n",
       " 0)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.suggest_symbolic(0,0,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4b8549a2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['x', 'x^2', 'x^3', 'x^4', 'x^5', '1/x', '1/x^2', '1/x^3', '1/x^4', '1/x^5', 'sqrt', 'x^0.5', 'x^1.5', '1/sqrt(x)', '1/x^0.5', 'exp', 'log', 'abs', 'sin', 'cos', 'tan', 'tanh', 'sgn', 'arcsin', 'arccos', 'arctan', 'arctanh', '0', 'gaussian'])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SYMBOLIC_LIB.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "cbde1924",
   "metadata": {},
   "outputs": [],
   "source": [
    "# add bessel function J0 to the symbolic library\n",
    "# we should include a name and a pytorch implementation\n",
    "add_symbolic('J0', torch.special.bessel_j0, c=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bda24c6d",
   "metadata": {},
   "source": [
    "After adding Bessel, we check suggest_symbolic again"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "83e5cfdd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  function  fitting r2   r2 loss  complexity  complexity loss  total loss\n",
      "0        0    0.000000  0.000014           0                0    0.000003\n",
      "1        x    0.001598 -0.002293           1                1    0.799541\n",
      "2      cos    0.159266 -0.250262           2                2    1.549948\n",
      "3      sin    0.159266 -0.250262           2                2    1.549948\n",
      "4    1/x^2    0.098715 -0.149929           2                2    1.570014\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "('0',\n",
       " (<function kan.utils.<lambda>(x)>,\n",
       "  <function kan.utils.<lambda>(x)>,\n",
       "  0,\n",
       "  <function kan.utils.<lambda>(x, y_th)>),\n",
       " 0.0,\n",
       " 0)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# J0 fitting is not very good\n",
    "model.suggest_symbolic(0,0,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e78f4674",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  function  fitting r2    r2 loss  complexity  complexity loss  total loss\n",
      "0        0    0.000000   0.000014           0                0    0.000003\n",
      "1       J0    0.999225 -10.314291           3                3    0.337142\n",
      "2        x    0.001598  -0.002293           1                1    0.799541\n",
      "3      cos    0.585822  -1.271642           2                2    1.345672\n",
      "4      sin    0.585822  -1.271642           2                2    1.345672\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "('0',\n",
       " (<function kan.utils.<lambda>(x)>,\n",
       "  <function kan.utils.<lambda>(x)>,\n",
       "  0,\n",
       "  <function kan.utils.<lambda>(x, y_th)>),\n",
       " 0.0,\n",
       " 0)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# This is because the ground truth is J0(20x) which involves 20 which is too large.\n",
    "# our default search is in (-10,10)\n",
    "# so we need to set the search range bigger in order to include 20\n",
    "# now J0 appears at the top of the list\n",
    "\n",
    "model.suggest_symbolic(0,0,0,a_range=(-40,40))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "47fb0d09",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "train loss: 1.31e-02 | test loss: 1.00e-01 | reg: 3.78e+00 : 100%|██| 20/20 [00:03<00:00,  5.16it/s]\n"
     ]
    }
   ],
   "source": [
    "model.fit(dataset, opt=\"LBFGS\", steps=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4773e989",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAFICAYAAACcDrP3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAxbElEQVR4nO3de1yUdb4H8M/vGZhhhgEGEBVvJEgq5v2CAgalSWa2pVm77dk93dvq2Dm+2q1sa2u7bRdL7XJsT3XOmqd92Tlqp7xhpYng/YLllbiIioCAMDjD3Od5zh82z0KZaT7DDPB5v1772mBuX5DffOZ3fYSiKAqIiIg0JIW6ACIi6noYLkREpDmGCxERaY7hQkREmmO4EBGR5hguRESkOYYLERFpjuFCRESaY7gQEZHmGC5ERKQ5hgsREWmO4UJERJpjuBARkeYYLkREpDmGCxERaS4i1AUQdQaKouDMmTOw2+0wm81ITEyEECLUZRGFLfZciC7AarVi8eLFSE9PR1JSEgYOHIikpCSkp6dj8eLFsFqtoS6RKCwJXomS6Pw2bNiA2bNnw+FwADjXewkI9FpMJhNWrlyJ/Pz8kNRIFK4YLkTnsWHDBsyYMQOKokCW5R+9nyRJEEJg7dq1DBiiNhguRN9jtVrRr18/OJ3OCwZLgCRJMBqNqK6uhsViCX6BRJ0A51yIvmfp0qVwOBwXFSwAIMsyHA4HPvzwwyBXRtR5sOdC1IaiKEhPT0dlZSUupWkIIZCamoqysjKuIiMCw4WoncbGRiQlJV3W4xMTEzWsiKhz4rAYURt2u/2yHm+z2TSqhKhzY7gQtWE2my/r8TExMRpVQtS5MVyI2khMTERaWtolz5sIIZCWloaEhIQgVUbUuTBciNoQQmDu3Lk/67GPPPIIJ/OJvsMJfaLv4T4XosvHngvR91gsFqxcuRJCCEjShZtIYIf+qlWrGCxEbTBciM4jPz8fa9euhdFohBDiB8Ndge8ZjUasW7cO06ZNC1GlROGJ4UL0I/Lz81FdXY1FixYhNTW13W2pqalYtGgRTp06xWAhOg/OuRBdBEVR8NVXX2HKlCnYuHEjrrnmGk7eE10Aey5EF0EIoc6pWCwWBgvRT2C4EBGR5hguRESkOYYLERFpjuFCRESaY7gQEZHmGC5ERKQ5hgsREWmO4UJERJpjuBARkeYYLkREpDmGCxERaY7hQkREmmO4EBGR5hguRESkOYYLERFpjuFCRESaY7gQ/QSv14tTp07hyJEjAICKigo0NTVBluUQV0YUvniZY6IfYbVasXLlSnz00Uc4dOgQbDYbPB4PoqKikJSUhMmTJ+Oee+5BdnY2IiIiQl0uUVhhuBCdx/bt2zFv3jx88803GD9+PGbMmIERI0bAbDbDarVi7969WL16NcrLy3H77bfjhRdeQFJSUqjLJgobDBei7/n8889x5513wmw24y9/+QtuuOEGeDweLF++HG63G7GxsfjlL38Jr9eL5cuX49lnn8WwYcOwbNky9OrVK9TlE4UFhgtRG99++y2uv/56REdHY/ny5cjIyIAQApWVlRgzZgxaWlowcOBA7N27F/Hx8VAUBcXFxbjjjjuQl5eH999/HwaDIdQ/BlHIcUKf6Dt+vx8vvfQSmpub8fbbb6vBciFCCOTk5ODVV1/Fp59+ioKCgg6qlii8MVyIvlNeXo7Vq1dj1qxZyMnJ+clgCRBC4Oabb8bEiRPx3nvvwefzBblSovDHJS5E39m2bRvsdjtmz56NqqoqtLa2qrdVV1fD7/cDADweDw4dOoTY2Fj19j59+mDWrFl49tlnUVdXh379+nV4/UThhOFC9J2jR4/CZDIhNTUVDzzwALZu3arepigK3G43AKCmpgbXXXedepsQAq+//jqGDx8Oh8OBmpoahgt1ewwXou84nU5ERETAYDDA7XbD5XKd936KovzgNp/PB6PR2C6EiLozhgvRd3r27Amn0wmr1YrMzExER0ertzmdTmzbtk0NkaysLHXjpBACAwYMQH19PSRJQnx8fKh+BKKwwXAh+s7YsWPh9Xqxa9cuvPLKK+1uq6ysxPjx49HS0oJevXrh448/hsViUW8XQuDJJ59E7969OSRGBK4WI1JNmDABqampWLp0KVpbW6HT6dr9L0AIAUmS1O9LkoTa2lqsWLECM2bMQFxcXAh/CqLwwHAh+k5iYiL+5V/+Bfv27cObb7550UuK3W43nn/+eTidTjzwwAMXvYSZqCvjsBhRG3feeSe2bNmCV155BSaTCQ8++CCioqIAABEREYiIiFB7MYqiwGaz4cUXX8Ty5cuxcOFCDB48OJTlE4UNHv9C9D0NDQ14+OGHsWbNGuTn52PevHkYOnQoSktLIcsy9Ho9Bg0ahF27dmHBggXYv38/nnvuOTz44IPths+IujOGC9F5tLa24r333sObb76J06dPIzU1Fenp6YiJiUFzczNKS0tRU1ODsWPH4plnnkFubi4kiaPMRAEMF6ILqKurw8aNG1FYWIivv/4au3btwuTJk5GdnY1p06YhMzMTJpMp1GUShR2GC9FF2r17NyZMmIDdu3dj3LhxoS6HKKyxH090kXQ6nboMmYgujK2EiIg0x3AhIiLNMVyIiEhzDBciItIcw4WIiDTHcCEiIs0xXIiISHMMFyIi0hzDhYiINMdwISIizTFciIhIcwwXIiLSHMOFiIg0x3AhIiLN8XouRBdJURTIsgxJkiCECHU5RGGNPReiS8BruRBdnIhQF0CkFUVRUFZWhjNnzoS6lMsiSRKuuuoqREdHh7oUop+Nw2LUZciyjIcffhj9+/eH2WwGAPj9fuh0uhBXdmmKiorw9NNPY8SIEaEuhehnY8+FuhSDwYB77rkHPXv2xIYNG/DWW29h7ty5yM/P7xTzJIqiwG63g5/5qLPjADJ1SW63G6+99hoKCgrwwAMPYNeuXXzDJupADBfqkk6fPo0DBw4AAKqrqzFv3jzYbLYQV0XUfTBcqEuqrq6GzWaDTqeDTqdDSUkJey9EHYjhQl1SVVUVvF4v0tLSkJ2dDY/Hg//+7/+GLMuhLo2oW2C4UJf07bffQlEUpKen4/7770dERATWr1+P8vLyUJdG1C0wXKjLURQFx48fBwCkpKQgPz8fQ4YMQWNjIz7++GMOjRF1AIYLdTmyLKsbKXv16oX4+HjccccdAIDly5ejsbExlOURdQsMF+pyfD4fmpqaAABJSUkAgDlz5qBXr14oLy9HQUEBey9EQcZwoS7H4/HAarVCCIEePXpACIGUlBTMnDkTfr8ff/vb3+ByuUJdJlGXxnChLsfpdKrLkBMSEgCcO6/rt7/9LUwmE3bu3MllyURBxnChLsdut8PhcCAyMhLx8fEAACEExo4di+zsbDidTixZsgRerzfElRJ1XQwX6nJsNhvcbjcMBgNiY2PV7xsMBjz44IMwGAxYs2YNvvzyS/ZeiIKE4UJdjtvtRlxcHHr06NHu2HohBPLz83HDDTfA6XTilVde4ZEwREHCcKEuZ9SoUdi5cycKCgqQmJjY7jaDwYDHH38cFosFO3bswIoVK9h7IQoChgt1OXq9Hn379sXAgQMREdH+qhJCCIwZMwa//vWv4fP58MYbb6C+vj5ElRJ1XQwX6nYkScIjjzyCvn374ujRo3j//fd55hiRxhgu1O0IIZCamorf/e53AIAlS5agtLSUw2NEGmK4ULckSRLuu+8+DB8+HLW1tXjllVe4NJlIQwwX6rZ69OiB+fPnQ6/XY8WKFfjiiy/YeyHSCMOFui0hBGbOnIkbb7wRTqcTzz77LBoaGkJdFlGXwHChbi0qKgp/+tOf0Lt3b5SUlGDBggXw+XyhLouo02O4ULcmhMCwYcPw+OOPQ5Ik/PWvf+XwGJEGGC7U7UmShLvvvhvTp0+H3W7HY489hpMnTzJgiC4Dw4UIQHR0NF5++WWkpKTg8OHDmD9/PpxOZ6jLIuq0GC5EODc8NmTIELz88sswGo1YsWIF3n33Xfj9/lCXRtQpMVyIviOEwC233IKHHnoIsizjxRdf5PwL0c/EcCFqIzIyEvPnz8eUKVNgtVrxyCOP4PDhwwwYokvEcCH6HovFgjfffBODBw9GRUUF7r//fk7wE10ihgvR9wghkJ6ejnfffRe9e/fGjh07cN9996Guro4BQ3SRGC5E5yGEwOTJk/HOO+8gPj4eX375Je677z7U19czYIguAsOF6EcIIXDTTTfh7bffRlxcHAoKCnDPPfewB0N0ERguRBcgSRJuu+02vPXWW4iLi8P69evxm9/8BlVVVQwYogtguBD9BEmS8Ktf/QpLlixBYmIivvrqK9x6663Yt28fA4boRzBciC6CJEm49dZb8be//Q39+vXD/v37MWvWLKxYsYIHXRKdB8OF6CJJkoTp06djxYoVGD16NKqrq3H33XfjT3/6E6xWK3sxRG0wXIgugRAC48aNwyeffII5c+bA4/Hgtddew5w5c1BSUgJZlkNdIlFYYLgQXSIhBPr164f//M//xEsvvYS4uDhs2rQJM2bMwBtvvMFeDBEYLkQ/ixACJpMJ8+bNwyeffILMzEw0NDRg/vz5uPHGG7F+/Xq43W6GDHVbDBeiyyBJEnJycrB69Wo8/fTTiI+Px/bt2zFnzhz85je/we7du+Hz+Rgy1O0wXIgukxACiYmJePrpp/H555/j1ltvhRACK1euRH5+Pu69917s2LGDPRnqVhguRBqRJAkjR47Ehx9+iBUrViAvLw9OpxPLli1Dfn4+Zs2ahf/93//FmTNnoCgKg4a6NIYLkYaEEDAYDMjPz8dnn32Gjz76CNdeey1kWUZBQQH+6Z/+CdnZ2fjDH/6A4uJinD17lkFDXRLDhSgIhBCIjo7GrFmzsHr1aqxfvx533XUXkpKSUF5ejoULFyI/Px/Z2dmYN28ePv/8c9TV1cHj8YS6dCJNRIS6AKKuTAiBqKgo5OTkICsrC9XV1SgoKMCqVauwZ88eHDlyBIcPH8aSJUvQs2dPTJkyBQMHDgx12USXjeFC1AGEENDpdEhJScH999+Pu+66C8eOHcPmzZuxbt067N69G/X19XC5XIiIYLOkzo9/xdSlKIqC5uZmREZGhrqUn5SUlIQ5c+bglltuQV1dHSoqKmA2m1FUVBTq0oguG8OFugwhBFJSUvDWW29Bp9OFupyfzel0Ii4uLtRlEF0WoXCZCnURXWnVlRACQohQl0H0szFciIhIc1yKTEREmuOcC9FFatvJ55AV0YWx50J0kUpKSqDT6VBSUhLqUojCHsOFiIg0x3AhIiLNMVyIiEhzDBciItIcw4WIiDTHcCEiIs0xXIiISHMMFyIi0hzDhYiINMdwISIizTFciIhIcwwXIiLSHMOFiIg0x3AhIiLNMVyIiEhzDBeii6AoCpqbmwEAzc3N4NXBiS6M4UJ0AVarFYsXL0Z6ejqmTp0KRVEwdepUpKenY/HixbBaraEukSgsCYUfwYjOa8OGDZg9ezYcDgeA81/m2GQyYeXKlcjPzw9JjUThiuFCdB4bNmzAjBkzoCgKZFn+0ftJkgQhBNauXcuAIWqD4UL0PVarFf369YPT6bxgsARIkgSj0Yjq6mpYLJbgF0jUCXDOheh7li5dCofDcVHBAgCyLMPhcODDDz8McmVEnQd7LkRtKIqC9PR0VFZWXtKKMCEEUlNTUVZWps7HEHVnDBeiNhobG5GUlHRZj09MTNSwIqLOicNiRG3Y7fbLerzNZtOoEqLOjeFC1IbZbL6sx8fExGhUCVHnxnAhaiMxMRFpaWmXPG8ihEBaWhoSEhKCVBlR58JwIWpDCIG5c+f+rMc+8sgjnMwn+g4n9Im+h/tciC4fey5E32OxWLBy5UoIISBJF24igR36q1atYrAQtcFwITqP/Px8rF27FkajEUKIHwx3Bb5nNBqxbt06TJs2LUSVEoUnhgvRj8jPz0d1dTUWLVqE1NTUdrelpqZi0aJFOHXqFIOF6Dw450J0ERRFwVdffYUpU6Zg48aNuOaaazh5T3QB7LkQXQQhhDqnYrFYGCxEP4HhQkREmmO4EBGR5hguRESkOYYLERFpjuFCRESaY7gQEZHmGC5ERKQ5hgsREWmO4UJERJpjuBARkeYYLkREpDmGCxERaY7hQkREmmO4EBGR5hguRESkOYYLERFpjuFC9BO8Xi9OnTqFI0eOAAAqKirQ1NQEWZZDXBlR+OJljol+hNVqxcqVK/HRRx/h0KFDsNls8Hg8iIqKQlJSEiZPnox77rkH2dnZiIiICHW5RGGF4UJ0Htu3b8e8efPwzTffYPz48ZgxYwZGjBgBs9kMq9WKvXv3YvXq1SgvL8ftt9+OF154AUlJSaEumyhsMFyIvufzzz/HnXfeCbPZjL/85S+44YYb4PF4sHz5crjdbsTGxuKXv/wlvF4vli9fjmeffRbDhg3DsmXL0KtXr1CXTxQWGC5EbXz77be4/vrrER0djeXLlyMjIwNCCFRWVmLMmDFoaWnBwIEDsXfvXsTHx0NRFBQXF+OOO+5AXl4e3n//fRgMhlD/GEQhxwl9ou/4/X689NJLaG5uxttvv60Gy4UIIZCTk4NXX30Vn376KQoKCjqoWqLwxnAh+k55eTlWr16NWbNmIScn5yeDJUAIgZtvvhkTJ07Ee++9B5/PF+RKicIfl7gQfWfbtm2w2+2YPXs2qqqq0Nraqt5WXV0Nv98PAPB4PDh06BBiY2PV2/v06YNZs2bh2WefRV1dHfr169fh9ROFE4YL0XeOHj0Kk8mE1NRUPPDAA9i6dat6m6IocLvdAICamhpcd9116m1CCLz++usYPnw4HA4HampqGC7U7TFciL7jdDoREREBg8EAt9sNl8t13vspivKD23w+H4xGY7sQIurOGC7U7VVVVWHz5s0oLi6Gw+GA1WpFZmYmoqOj1fs4nU5s27ZNDZGsrCx146QQAgMGDEB9fT18Ph/Ky8sxfvx4REVFhepHIgo5LkWmbufkyZMoLCxEYWEhNm/ejBMnTkAIgZSUFFRUVOCdd97Bvffe2+4xlZWVGD9+PFpaWnDFFVdgz549sFgs6u1CCDz55JNYsGABdDodoqKikJmZiauvvhp5eXkYP348lyhTt8JwoS6vpqYGmzdvxpYtW7B582YcO3YMADBixAj1zT8nJweyLCMnJwfx8fEoKChoN2H/Y/tcgHPDZDU1NcjNzcXMmTNx5513YsuWLdiyZQuKiopgtVphNBoxceJE9fXGjh0LvV4fkt8HUUdguFCXc/r06XZhUl5eDgAYNmyY+uY+efJkJCYm/uCx77zzDh599FE89dRTeOKJJ9ShrwuFi8vlwr/9279h9erV2LRpEwYPHqw+n9/vxzfffIMtW7agsLAQxcXFOHv2LEwmEyZNmoTc3Fzk5eVh9OjRiIyM7IDfDlHHYLhQp9fQ0IDCwkI1TEpLSwEAQ4YMQW5uLnJzc3H11Vdf1Nlfra2tuPvuu7Fu3Tr8+c9/xoMPPoioqCgcO3YMEyZMUIfFdu3aBYvFApvNhhdffBF//etfsXDhQtx1110XfH6fz4f9+/er9W7duhV2ux1msxlZWVlqvaNGjeJhmNSpMVyo02lsbERRUZE6b3L48GEAQHp6uvrmnJub+7PP+WpoaMDDDz+MNWvWID8/H/PmzcPQoUNRWloKWZah1+sxaNAg7Nq1CwsWLMD+/fvx3HPP4cEHH4ROp7uk1/J6vSgpKVF7Wtu2bYPD4UBsbCyys7PVn2XEiBGX/NxEocRwobDX3NyMoqIi9Q34wIEDAIDU1NR2YdKnTx/NXrO1tRXvvfce3nzzTZw+fRqpqalIT09HTEwMmpubUVpaipqaGowdOxbPPPMMcnNzIUmXf+CFx+PB3r171eDcvn07XC4XLBYLcnJy1J/1qquu0uT1iIKF4UJhp6WlBcXFxWqYfP3111AUBSkpKcjLy1PnTTpio2JdXR02btyIwsJCVFZWwuVyIT4+HldddRWmTZuGzMxMmEymoL2+2+3G7t271bDZuXMn3G43EhISMHnyZDVsLuYcNKKOxHChkLPZbNi6dav6BlpSUgJZltG3b1/k5eWpk94pKSkhrdPv90NRFEiSFLJeg8vlws6dO9X5pV27dsHr9aJHjx64+uqr1bAZPHgww4ZCiuFCHc5ut2P79u1qz2Tv3r3w+/1ITk5u1zMZOHAg3yB/gsPhwI4dO9Sw2bNnD3w+H3r27NluyHDQoEH8XVKHYrhQ0AXeAANhsnv3bvUNsG2Y8A3w8tntduzYsUPdILpv3z41uNuGDYObgo3hQpoLDN0EwmTnzp3wer1ISkrC1VdfrYYJh26C7+zZs9i+ffsPhhz79eun/jvk5uaGfMiRuh6GC102t9uNXbt2nXfSue08ACedQ89qtWLr1q3qps62iyXahg1PdabLxXChS+bxeLBnzx516GXHjh3qctm2K5i4XDb8NTU1tVtMEVjmPXDgQPXfMS8vD8nJySGulDobhgv9JK/Xi3379qlhsn37dnWjX05Ojjpvwo1+nd+ZM2fOu0F10KBBmmxQpe6D4UI/0PaIks2bN6tXaDSbze12jfOIkq6voaFBHUIrLCxUj9YZPHhwu7Dp0aNHiCulcMNwIfVwxUCYbN26VT1cMSsrS+2ZjBkzhocrdnN1dXXtwiZwKGhGRka7c9wSEhJCXCmFGsOlG5JlGQcPHlTDpLi4GFarFVFRUZg0aZIaJuPGjeOx8HRBNTU1atAUFhaqlzMYPny4GjaByxhQ98Jw6QYURcHhw4fVMCkqKkJTUxMMBgMyMzPVMJkwYQIvaEWX5eTJk+qGzsLCQvVCbCNHjlTDJjs7G3FxcaEulYKM4dIFKYqC0tJStYEXFRWhoaEBkZGRyMzMVBt5ZmYmL8VLQVVVVdUubE6dOgVJkjB69Gj17zArKwsxMTGhLpU0xnDpAhRFQXl5udqACwsLUV9fj4iICIwfP15txBMnTgzqIYtEF6IoCo4dO9buEtN1dXXQ6XQYM2aMusdm0qRJiI6ODnW5dJkYLp1QoJG2DZPa2lrodDqMHTtW3ZswceJEmM3mUJdLdF6BD0Vt52wCH4rGjRunhk2wT56m4GC4dBLHjx9vFybV1dXthhfy8vIwadKkdtd9J+pMAsO5bcPmzJkz0Ov1GD9+vDo3yOHczoHh0kkMHz4cZWVl6sRoXl4esrKyYLFYQl0aUVDIsowjR46oQRNYiLJs2TLMmTMn1OXRT2C4dBKyLEMIwbO5qNtSFAWKorAddBIMFyIi0hzP7tBIYHLyzJkzoS7lskiShGHDhnG1Dl0ytgFqi+GiEUVRsHjxYvTr10+zFVoXMxTm9/shSZJmwwTFxcX44x//iOHDh2vyfNR9BKMNfJ/f74csy4iIiAja0BjbgDYYLhoyGAy46667NDkxtq6uDk899RR69+6N+fPn/+BTlKIo2LdvH55//nnMnDkTd91112Ufb68oCux2OzhSSj+Xlm0AOPc36XK5sGvXLqxZswYHDx6Ex+PBgAEDMHXqVOTn5yMxMVGzoGEb0A7DJQwpioJVq1bho48+QkREBPLy8jB16tR295FlGa+//joKCgpw6NAh3HDDDbzmBnUpsiyjpKQEf/7zn1FYWAi3263etnXrVnz88ccYMmQInnjiCdx8882IjIzkRH8YYbiEIUVRsHPnTiiKAq/Xi927d2PKlCntGs7Zs2fx9ddfAzjXyykvL2e4UJfh9/vx8ccfY/78+Th9+jSio6NxzTXXIDc3F9HR0SgpKUFBQQEOHz6M++67D8XFxXjmmWeQkJDAgAkTDJcw5PF4UFFRoX5dWlqqLsEMaGxsRGNjI4BzF/OqrKxETk4OGxZ1en6/H8uWLcPvf/97tLa2IjMzE88//zwmTpyoXvIhcErFyy+/jI8//hjvvfceTp48iXfffRc9e/ZkOwgDvAZtCCmKAqvVirNnz7Yb47Xb7aivr1e/PnnyJHw+X7vH1tXVweFwqF9XVla2u93v96OxsREej4fjx9RpKIqCdevW4fHHH0drayt+8Ytf4H/+538wefJk6PV6dYGLJElITU3F22+/jddeew2xsbFYv3495s6d+4P2RKHBcAkRRVFw5MgRXHfddbjppptQXV2t3ma1WtHS0qJ+XV9fD5fL1e7xtbW17QLn5MmT7Z572bJlyMrKwlNPPfWDYCIKR4E28eijj6KlpQXTpk3D22+//aM9ESEEDAYD7r33XixcuBAxMTFYs2YNFixYAL/fH4KfgNpiuARZYLWLy+X6waepZcuW4cCBA9ixYwdWr16t3t7Q0ACn06ne7+zZs7Db7e2es7a2tt3z1dTUqCFit9uxZMkSnDhxAkuXLkVZWVm71/X7/WhqaoLX6+UnPAobNpsNjz/+OE6cOIGMjAwsXrz4olaC6XQ63H777XjssccghMCSJUuwceNG/m2HGMMliBRFwdGjRzFz5kzMnj27Xe/C5XJh69at6tfbtm1TG0N9fT18Ph90Oh0kSUJrayvOnj3b7rnr6uoAQL2GfWAIDACqq6vVYbKzZ89i9+7d6nPLsoz3338fkyZNwmOPPdZuBQ5RqMiyjA8++ACbNm1CbGwsFixYgJSUlIueO9HpdHjooYdw/fXXw2634+mnn243tEwdj+ESRLIs44033kBRURE2btyIpUuXqm/yjY2N6iVhgXOT9oHeyunTpyHLMvr06YPY2Fh4PB40Nze3e+6GhgYAQGpqKnQ6HVpaWtQ5mLKyMrS2tgL4x36YgMbGRixcuBDHjx/Hf/3Xf6mr0ohCJfAhbPHixZBlGb/73e+Qm5t7yZPyJpMJzz33HJKTk3HgwAG8+eabHB4LIYZLEDU0NGDz5s3q15s3b1Z7F8ePH283r3L69Gk0NzdDURScPn0aAJCSkgKLxQKfz4empib1vrIsq0dsDB06FHq9HjabTd38VVZWBlmW1fsfOXIEXq8XAHDw4EGcOnUKwLne09q1a4PzwxNdJK/Xi1dffRV1dXUYOXIk5s6dC51Od8nPI4RARkYGHn30UUiShA8++AB79uzhh6cQYbgE0aFDh1BbW6t+XV5ejoaGBiiKgoqKCng8HsTHx8NoNOLs2bPqUFfgMf3794fFYmkXJgDg8/nUnkx6ejqMRiNcLpf6vcAcS+/evSFJEk6cOKEGT0lJCbxer/qpcPv27T9YLEDUURRFQXFxMT777DMYDAY88cQT6NGjx89+PiEEfvvb3yIrKwtWqxWvvPIK/75DhOESJIqiYPfu3fD5fBg4cCDi4uLQ3NyM48ePAzg3DAYAI0aMQHJyMjweD6qrqyHLstpz6devn3q9ljNnzqifwDweD6xWK4QQuOKKK2A2m+H1etHY2Aifz6e+Rl5eHoxGI86cOaOG2oEDB9TX1ev1qKioUEONqKO53W4sWrQIDocDU6ZMwfXXX3/Ze1RiYmIwf/58mEwmbNy4kZP7IcJwCZLA0RUAcO211yIlJQUejwfffvstZFlWexejRo1C3759Icsyqqqq4PP51M2Rffv2RUJCAgCo3wPODWfZ7XbodDr06dMH8fHx8Pv9OH36NFwuF2prayGEwIQJE5CQkACHw4Hjx4/D7Xarrztz5kz07NkTLS0tOHr0aEf+aogA/KPXsmXLFphMJvzrv/4rDAbDZT+vEAI5OTmYMWMG3G433nrrrXarL6ljMFyCxG6348iRIxBCICsrC4MGDQJwbv7D5XKhqqoKAJCRkYEBAwYAAKqqquByudDU1AQhBHr37o34+HgAaDfn4nA44HQ6ERERgR49eiAxMRHAuXmblpYWNDU1ITIyEhkZGejbty/8fj/KyspgtVpRU1ODiIgIZGZmIi0tDX6/H/v37+cnO+pwXq8X//Ef/wGXy4VrrrkGEydO1GxnfWRkJB566CGYzWbs2LEDRUVF/BvvYAyXIKmpqUFtbS1MJhOGDx+OIUOGADg3HNbY2IjTp09Dr9cjLS0NV1xxBYBzk/w2mw1nz55FREQEkpKS1J5LYLIfOBcubrcbkZGRsFgs6gm0dXV1aGxshM1mg8lkQv/+/ZGamqq+bk1NDZqbm2E2m5GWloaRI0cCAPbv399uAQBRsCmKgq+//hqbNm2CwWDA/fffD71er9nzCyEwduxYXHvttXC73fjggw/URS3UMRguQaAoCg4fPozW1lb06tUL/fv3x5AhQyCEwLFjx1BWVoaWlhbExMSgX79+SElJAXAukBoaGtDa2gqDwYDExEQ1XKxWq7qs0mazwePxICoqCiaTCb179wZwLlxOnToFt9sNi8WC+Ph4XHnllQDOTfKXlpbC7XajZ8+e6NGjB0aPHg0hBI4cOdJukyZRsMmyjKVLl8Jut2PcuHFBORcvMjIS99xzD/R6Pb766iscOnSIvZcOxHAJkr1790KWZQwZMgQxMTFIS0uDwWBAfX09iouL4Xa7kZycjB49emDAgAGIiIhAY2Mjqqqq4Ha7YTabERcXpw6LnT17Vg2XwH8bjUYYjUY1XOrr61FRUQFZltGrVy9ER0dj6NChkCQJVVVV2LVrFxRFwRVXXAGTyYSMjAwYjUbU1taipqYmZL8r6n5OnDiBNWvWQJIk/PM//zNMJpPmryGEQHZ2NkaPHg2bzYa///3vDJcOxHAJgtbWVuzcuRMAMH78eOh0OnVy3maz4f/+7/8AAEOGDFF7HkajES0tLTh48CB8Ph8sFgvMZjMsFguEELDb7eoemebmZsiyjOjoaBgMBiQnJ0MIgYaGBnVyvn///oiMjMSVV16J6Oho1NXV4fPPPwcADBs2DDqdDgMGDECfPn1gt9vbbbQkCiZFUbBy5UqcPn0aaWlpmD59etBOMY6Ojsavf/1rCCHw6aefttsaQMHFcNHYnj17MHv2bOzatQt6vR5ZWVkQQiA+Ph4DBw6E3+9XA2DMmDEQQiAxMREWiwVOpxM7duyAoihISkpCVFQU4uLioNPp1HkW4NwQmaIoiImJQWRkJPr06aP2fL755hsA53buCyHQt29f9O3bF06nE5WVlRBCqJdvjY2NxZgxY6AoCl544QUsXrxYDTCiYLFarVi+fDkURcFtt912WftafooQAjNmzED//v1RXV2NgoIC9l46CMNFY5IkYd++ffB6vRg0aBBGjBgBANDr9Rg3bpx6v6ioKEyaNAlCCMTExCA5ORl+vx/bt28HAHWoLCYmBjqdDk6nE06nUz2mH4B6W+/evWEymdDU1ISDBw8CgLo6zWw2Y/z48errms1mjBw5Uj22fMaMGdDpdDh27Bg2b97M62BQUCmKgsLCQhw9ehSJiYm47bbbgv43l5ycjBtvvBGyLGP58uXcVNlBGC4aGz58OObOnYtRo0bh6aefVjdBAsD06dNhNBoBACNHjlR7EHq9Hunp6QCgng+WlpYG4FwYREZGwuPxqGv1A+ESFxcHSZKQkJCAhIQE+P1+uFwuGAwGpKenq9e+uOWWW9SVOKNGjVJXkAkhMH36dPzqV7/C6NGj8Yc//EG9GBNRMPh8Pvz973+Hz+fDlClT1L/zYBJC4LbbboPJZMLevXu59L6DMFw0FhkZiT/+8Y/YtGkTbr75ZvVTWWC/y5NPPonp06fj5ZdfRnR0tHrbqFGj1OeQJAnDhg2DEAImkwkGgwFer1dd0dU2XAI9n4EDB6qPT0xMVJc3CyFw7bXX4vHHH8f111+PF154AVFRUep9Y2JisGTJEnz55ZfIysoK4m+GCKioqEBRUREiIyNxxx13/KwzxC6VEAIjR47E2LFj4XA4sGrVKoZLB+BljoNAp9Odd/WLXq/H73//e8iyDEmS2gXPxIkTYTKZ4HA4EB8fr/ZqoqKiEBUVhZaWFvV8sMCBl4GVZBEREZg4cSI2bdoE4NzRLm3HsQ0GA5588kn4/X7odLp2wxBCCERGRiIyMpINjoJKURSsXbsWzc3NGDZsmKabJn9KVFQUZs+ejeLiYvVKl8Gc6yH2XDqcEOIHb/AAcNVVV2H69OmIiorCLbfcgv79+wM4FwxGoxGyLKOlpUX9f+BczyXwnLNnz8aAAQNgsVhw9913q9d5CdwuhEBERATnVChknE4nPvvsMwDAjTfeqP79dgQhBPLz85GUlITjx4+3u34SBQd7LmEiKioK//7v/46ysjIMGTJEDQe9Xg+z2QxZltWNlDabDcC5nksgLDIyMrB+/Xo4HA5kZGQwRCjsHD58GAcPHoTJZMKNN97Y4X+j/fv3R3Z2Nj755BN8+umn6mIWCg72XMKEEAJxcXEYN24czGaz2vACK8YAoKWlBT6fDzabTb1/28enpaVh+PDhbDAUdhRFwYYNG9Da2ophw4YhIyOjw2vQ6XT4xS9+AUmSsGXLFl6pMsgYLmFOkiTExsYCODeR7/F40Nra+oNwIQpnTqcTGzZsAADk5+cHZUf+TwmclpycnIyamhoOjQUZwyXMSZKkhojVaoXL5YLL5YJOp1N7NEThrqysDIcPH4bJZEJ+fn7Ihm2Tk5ORlZUFv9+PtWvX8sDWIGK4dAKBcGlpaYHT6YTL5Wo3XEYUzhRFwcaNG2G323HllVdi6NChIaslsHFYkiQUFxejoaEhZLV0dQyXTiCw5NhqtaonIuv1enWfDFE483g8KCgoAABMnTo1pH+3gf1mPXv2xKlTp7B7924OjQUJw6UTCOzyb2lpUSf1o6Ki1N3+ROGsqqoK33zzDQwGA6ZNmxbqctCnTx9MmDABfr8f69evZ7gECcMlzAkh1HCx2+1obm6Gz+eD0Whst9OeKBwFzhJraWnBwIEDMWLEiJAvk9fpdOpJzEVFReqJF6QthksnYLFYIEkSbDYbGhoaoCiKetw+UTjz+XzqScS5ubnqysdQCqwas1gsOH78OA4cOBDqkrokhksnEDig0ul0ora2FoqiIDY2tt0ufKJwVFtbi7179yIiIiKkq8S+b8CAARgxYgQ8Hg+++OILDo0FAcOlEwgEidPpxKlTpwCcCxyGC4UzRVGwc+dONDQ0IDk5Wb1+UTjQ6/WYOnUqAGDz5s3qieOkHYZLJxATEwO9Xg+Xy4UTJ04AOLeCTJL4z0fhS1EUfPHFF5BlGZmZmUhKSgp1SSohBPLy8mA0GlFaWoqKiopQl9Tl8N2pEzCZTDAajfB6vaiqqgLwj+XJROGqpaUF27ZtgxAC1113Xdh9GBo8eDBSU1Nhs9lQXFzMoTGNhde/Np2X0WiEyWSC1+vFyZMnAQAJCQkhrorowg4fPoyTJ08iLi5OvepqOImJiUFOTg4AYNOmTfD7/SGuqGthuHQCUVFR6i59r9cLAGE1xED0fYqiYMuWLXC73cjIyMCAAQNCXdJ5XXvttdDpdCgpKUFjY2Ooy+lSGC6dQERERLvLJQshkJSUFHafBIkCvF4vCgsLAQCTJ08Oy2XzQgiMHj0aPXr0QF1dHZcka4zh0gnodLp2PRWdTofExMQQVkR0YbW1tTh8+DAiIyNx9dVXh7qcH9W7d29cddVV8Pl8ahiSNhgunYAQAsnJyerXBoOBl2ilsPb111/jzJkz6NWrF4YNGxa2vey24bd161a43e4QV9R1MFw6ib59+6r/bTab2w2TEYUTRVGwdetW+P1+jBgxIqw/CAUOstTr9SgtLVX3kdHlY7h0AkIIpKSkqEs5ExMTeaEwCltutxs7d+4EAGRnZ4f9lVGHDh2KPn36wGq1oqSkJNTldBkMl04iLS1NPQV50KBBPBGZwlZtbS3Kysqg1+uRmZkZtkNiAfHx8Rg1ahRkWUZxcXGoy+kyGC6dRGpqKjIzMxEdHY2bbrop7DakEQUcO3YMPp8PvXv3xpVXXhnqcn6SJEnIycmBTqfD0aNH4fF4Ql1Sl8DDqTSkKAqsVisiIyOD8vwLFy5EbW0thg4diubm5qC8Bic06XIoioJhw4Zh5cqVqK+vhxACTU1NoS7rJ02ePBnvvvsuRowYgVWrVoW6nC6B4aIRIQQGDBiAd955J+hjzJ9++mnQntvpdIbFsejU+QTawPvvv6+2gaKiohBXdWn27dvHNqARofBAHU0oitJlziYSQoT9ODmFH7YBaovhQkREmuOsMBERaY7h0kkoigJZlrvMsAPRz8F20HkwXDqJ/fv3w2g0Yv/+/aEuhShk9u/fD5PJxHbQCTBciIhIcwwXIiLSHMOFiIg0x3AhIiLNMVyIiEhzDBciItIcw4WIiDTHcCEiIs0xXIiISHMMFyIi0hzDhYiINMdwISIizTFciIhIcwwXIiLSHMOFiIg0x3DpBBRFQXNzMwCgubmZF0qibinQDtr+P4UvhksYs1qtWLx4MdLT0zFlyhR4PB5MmTIF6enpWLx4MaxWa6hLJAo6toPOSSiM/7C0YcMGzJ49Gw6HAwDafUoTQgAATCYTVq5cifz8/JDUSBRsbAedF8MlDG3YsAEzZsxQrxf+YyRJghACa9euZcOiLoftoHNjuIQZq9WKfv36wel0XrBBBUiSBKPRiOrqalgsluAXSNQB2A46P865hJmlS5fC4XBcVIMCAFmW4XA48OGHHwa5MqKOw3bQ+bHnEkYURUF6ejoqKysvaSWMEAKpqakoKytTx6GJOiu2g66B4RJGGhsbkZSUdFmPT0xM1LAioo7HdtA1cFgsjNjt9st6vM1m06gSotBhO+gaGC5hxGw2X9bjY2JiNKqEKHTYDroGhksYSUxMRFpa2iWPFwshkJaWhoSEhCBVRtRx2A66BoZLGBFCYO7cuT/rsY888ggnMalLYDvoGjihH2a4vp+I7aArYM8lzFgsFqxcuRJCCEjShf95AjuTV61axQZFXQrbQefHcAlD+fn5WLt2LYxGI4QQP+jmB75nNBqxbt06TJs2LUSVEgUP20HnxnAJU/n5+aiursaiRYuQmpra7rbU1FQsWrQIp06dYoOiLo3toPPinEsnoCgKmpqaYLPZEBMTg4SEBE5aUrfDdtC5MFyIiEhzHBYjIiLNMVyIiEhzDBciItIcw4WIiDTHcCEiIs0xXIiISHMMFyIi0hzDhYiINMdwISIizTFciIhIcwwXIiLSHMOFiIg0x3AhIiLNMVyIiEhz/w/mJpPJqef7VAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 500x400 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.plot()"
   ]
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}