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.40e-01 | test loss: 1.38e-01 | reg: 2.88e+01 : 100%|██| 20/20 [00:30<00:00,  1.50s/it]\n"
     ]
    }
   ],
   "source": [
    "from kan import KAN, create_dataset, SYMBOLIC_LIB, add_symbolic\n",
    "import torch\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,5,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.train(dataset, opt=\"LBFGS\", steps=20, lamb=0.01, lamb_entropy=10.);"
   ]
  },
  {
   "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": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAFICAYAAACcDrP3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAyXUlEQVR4nO3de3Ac9WEH8O/u3ekeOj1OJ+kkW5ZsyZLwg7eNIaQOSROU1gyZCWFa2mnTR9q0TDGBkqRO/20CzQCxcYZpm2kyZui0pDXjTGuCMy7haXCMbWKDDZEs2dLp/fDpdLo73WN//cPsdn93JyPbJ+2e9P3MMPGt7qSflN397u+tCCEEiIiIiki1ugBERLT8MFyIiKjoGC5ERFR0DBciIio6hgsRERUdw4WIiIqO4UJEREXHcCEioqJjuBARUdExXIiIqOgYLkREVHQMFyIiKjqGCxERFR3DhYiIio7hQkRERee0ugBEpUAIgcnJScRiMfj9fgSDQSiKYnWxiGyLNReiy4hEItizZw/a29tRV1eHdevWoa6uDu3t7dizZw8ikYjVRSSyJYU7URIVdujQIdx3332Ix+MALtVedHqtxefzYf/+/ejq6rKkjER2xXAhKuDQoUPYsWMHhBDQNG3e96mqCkVRcPDgQQYMkQnDhShHJBJBU1MTEonEZYNFp6oqvF4vwuEwqqurF7+ARCWAfS5EOfbt24d4PL6gYAEATdMQj8fx3HPPLXLJiEoHay5EJkIItLe3o7e3F1dyaSiKgtbWVnR3d3MUGREYLkSSiYkJ1NXVXdPng8FgEUtEVJrYLEZkEovFrunzMzMzRSoJUWljuBCZ+P3+a/p8RUVFkUpCVNoYLkQmwWAQbW1tV9xvoigK2traUFNTs0glIyotDBciE0VR8NBDD13VZ3fu3MnOfKKPsUOfKAfnuRBdO9ZciHJUV1dj//79UBQFqnr5S0Sfof/iiy8yWIhMGC5EBXR1deHgwYPwer1QFCWvuUs/5vV68dJLL+Huu++2qKRE9sRwIZpHV1cXwuEwdu/ejdbWVulrra2t2L17NwYHBxksRAWwz4VoAYQQOHbsGB5++GHs2bMHW7duZec90WWw5kK0AIqiIBAIwOPxIBAIMFiIPgHDhYiIio7hQkRERcdwISKiomO4EBFR0TFciIio6BguRERUdAwXIiIqOoYLEREVHcOFiIiKjuFCRERFx3AhIqKiY7gQEVHRMVyIiKjouOQ+0QLF43F89NFH6OzshM/ns7o4RLbGcCFaIE3TkEwm4fF4PnH7Y6KVjuFCRERFx8cvIiIqOqfVBSDSZbNZnDlzBrOzs1YXpeR1dHSgpqbG6mLQCsZwIdvIZDLYu3cvUqkUvF6v1cUpWeFwGH/3d3+HO++80+qi0ArGcCHbEEJAVVU8+OCDuOmmm6wuTsno7u7G2NgYPvWpT0HTNHz729+2ukhEDBeyH4fDgbKyMquLYXtCCPzmN7/B448/jkgkgnQ6jbvuuguKolhdNCKGC1EpEkKgu7sb//AP/4DBwUEAwDPPPAOXy2VxyYgu4WgxohIjhEBvby+++93vIhwOG8crKytRW1trYcmI/h9rLkQlRAiB/v5+fPe730V/f79xvKGhAX//93+P9vZ2C0tH9P9YcyEqEUIIjI6O4oknnkBfX59xvL6+Ht/5znewceNG9reQbTBciEqAEAKRSARPPvkkPvzwQ+N4MBjErl27sHnzZgYL2QrDhcjmhBBIJBLYu3cvTpw4YRyvqqrCt771Ldx4440MFrIdhguRzaXTafzkJz/Ba6+9Bn0pwPLycjz88MPYsmULg4VsieFCZGPZbBYHDhzAz372M2iaBgAoKyvD1772NWzfvp2rM5Nt8cwksikhBF5//XXs27cP6XQaAKCqKn7v934P99xzDxwOh8UlJJofw4XIhoQQOHXqFPbu3Yt4PA4AUBQFXV1d+IM/+AM4nZxFQPbGcCGyGX0uy1NPPYWLFy8ax2+99VZ8/etfh9vttrB0RAvDcCGyESEEpqam8OSTT2JgYMA4vn79ejz66KOorKxkBz6VBIYLkU0IITA7O4tnnnkG77//vnG8vr4ejz32GBoaGhgsVDIYLkQ2kU6n8eMf/xhvvvmmcayiogKPPPIIOjo6GCxUUhguRDaQzWbx05/+FP/93/9tDDl2u934+te/jttuu43BQiWH4UJkMU3TcOjQITz//PPIZDIALu1p88ADD6Crq4tzWagk8awlspAQAu+++y7++Z//GXNzcwAuzWXZsWMHfv/3f59DjqlkMVyILKJv+PX0008jGo0ax++88078xV/8BXfjpJLGcCGygBACY2NjePLJJzE6Omoc37RpE3bu3Iny8nL2s1BJY7gQWWB2dhZ79uxBd3e3caypqQmPPfYYgsEgg4VKHsOFaImlUin85Cc/wdGjR41j1dXVeOyxx9DS0sJgoWWB4UK0hDRNw0svvSQNOfZ4PHjwwQdx/fXXM1ho2WC4EC0RIQSOHz+OH//4x8Yqx/qQ48997nMcckzLCs9moiWgL0a5Z88ezMzMALi0yvHnPvc53H///Vw+n5YdhgvRIhNCYHp6Gk8//TQGBweN4xs3bsRf/dVfwePxWFg6osXBcCFaZKlUCv/yL/+C06dPG8caGhrw6KOPIhAIsJ+FliWGC9Eiymaz2L9/P37xi19ACAEAKC8vx86dO7Fu3ToGCy1bDBeiRSKEwBtvvIHnn38e2WwWAOB0OvEnf/InXIySlj2GC9EiEELg7Nmz+OEPf4hEIgHgUgf+7/7u7+Lee+9lBz4tewwXoiITQmB4eBhPPfUUJicnjeNbt27Fn//5n8PlcllYOqKlwXAhKiIhBKLRKJ5++mn09vYax9etW4eHH34YFRUVbA6jFYHhQlREyWQSzz77LE6cOGEcCwaDePTRR9HY2MhgoRWD4UJUJKlUCvv27cPhw4eNkWE+nw8PPfQQNm3axGChFYXhQlQE+pDjF1980VgzzOVy4U//9E/x6U9/msFCKw7DhegaZbNZHDx4EM8995yxZpiqqrjvvvs4MoxWLIYL0TXQNA2HDx/GP/3TPyGZTAK4NOS4q6sLf/zHf8zdJGnFYrgQXSVN0/Daa69h79690lyWT3/60/jrv/5rrhlGKxrDhegqCCHw9ttvY/fu3ZidnTWOb926FY888gj8fj/7WWhFY7gQXSEhBE6cOIGnnnoK0WjUOH7zzTfjm9/8JqqrqxkstOIxXIiugBACH3zwAb7//e/j4sWLxvFNmzbhW9/6FoLBIIOFCAwXogUTQqCnpwePP/44xsbGjOPr16/Hrl27EAqFGCxEH2O4EC2AEALnz5/H9773PQwNDRnHW1pa8J3vfAerVq1isBCZMFyIPoEQAkNDQ3j88cdx/vx543hjYyN27dqFtWvXMliIcjBciC5DCIHR0VF873vfQ3d3t3G8rq4Ou3btQkdHB4OFqACGC9E8hBCYmJjAE088gTNnzhjHA4EAvvnNb2Lz5s0MFqJ5MFyIChBCYHp6Gk8++SR+/etfG8crKyvxt3/7t9iyZQuDhegyGC5EOYQQmJ2dxZ49e3Ds2DHjuN/vxze+8Q3ccccdDBaiT8BwITIRQhh7srz++uvG0vlerxd/8zd/g8985jNQVV42RJ+EVwnRx4QQmJubw49+9CMcOnTIWDrf7XbjL//yL/H5z3+ewUK0QE6rC0BkB3qN5Uc/+hF+9rOfSXuy/NEf/RHuueceLp1PdAUYLrTi6fveP/vsszh8+LARLA6HA/fffz/uv/9+OJ28VIiuBK8YWtH0CZJPP/00Tp48afSxOBwOfPnLX+aeLERXieFCK5YQAufOncPjjz+O3t5e47jL5cL999/PYCG6BgwXWpGEEDhz5gyeeOIJhMNh47jP58Of/dmf4d5774XL5bKwhESljeFCK46+H8v3v/99aXXjmpoaPPzww7jzzjvZeU90jRgutGIIIZDNZvHaa6/hhz/8ISKRiPG1hoYGfPvb38aNN97ICZJERcBwoRVB0zQMDg7ihRdewOHDhzE3N2d8rbm5Gbt27UJnZyeDhahIGC60rAkhEIvFcODAARw4cABTU1PS19vb27lsPtEiYLjQsiWEQH9/P55++mmcPn3aGGYMAKqq4rbbbsPDDz/MHSSJFgHDhZYlIQQ+/PBD/OM//iMuXLggfa2+vh5f+cpXsGPHDni9XgYL0SJguNCyo2kajh8/jieffFIaDebz+fDFL34RX/nKV9DQ0MBQIVpEDBcqeUIIY22wkZERvPrqqzhw4ABmZmaM99TV1eGRRx7BbbfdBlVVGSxEi4zhQiVL76w/ceIEjhw5go8++gijo6NIJpPS+5qamrBr1y5s2LCBoUK0RBguVJKy2SyOHDmCffv2oa+vz1hs0kxRFGzYsAGPPfYYR4MRLTGGC5WcVCqF//iP/8C///u/59VSdFVVVfjCF76ABx54AIFAgMFCtMQYLlRSstks/vM//xPPP/880um0cdzlcmHVqlXo7OzE9ddfj5tvvhkNDQ1cxoXIIgwXKhmapuGVV16RgkVVVdxyyy34wz/8Q3R0dMDj8UBRFNZUiCzGcCFbE0JgZmYGo6Oj6Ovrw7PPPms0hSmKgq6uLjz44IMoLy9noBDZCMOFbEsIgZMnT2Lv3r0YGhpCJpOROu5vv/12BguRTTFcyLYmJibwgx/8QNpvRbdhwwZ84xvfYLAQ2RTDhWxJCIHXX38dg4OD0nFVVXHTTTfh0UcfRV1dHYOFyKYYLmRLmUwGb7/9trHYZHl5Oe655x5s2LABW7Zsgc/nY7AQ2RjDhWwpEong/PnzxusNGzbga1/7GpxOnrJEpUC1ugBEhQwODiIajRqvN2/ezDkrRCWE4UK2dP78eWkuy/r169kMRlRCGC5kSyMjI3C5XAAAr9eLpqYmi0tERFeCDdhkSw888AC2b9+O8+fPY2JiAnV1dVYXiYiuAMOFbEXTNJw5c8aYhV9VVYWqqiqcPHnS4pKVhmw2K/VVEVmF4UK2oaoq1q1bh6NHj+Lo0aNWF6dklZeXw+/3W10MWuEUoU8kILKYEALZbBY8Ja+dw+GAqrJLlazDcCEioqLjow0RERUdw4WIiIqO4UJEREXHcCEioqJjuBAtkKZpiMfj0oZlRFQYw4Vogc6dO4cdO3bg3LlzVheFyPYYLkREVHQMFyIiKjqGCxERFR3DhYiIio7hQkRERcdwISKiomO4EBFR0TFciIio6BguRERUdAwXIiIqOoYLEREVHcOFiIiKjuFCRERFx3AhIqKiY7gQLYAQAlNTU0gmk5iamoIQwuoiEdkaw4XoMiKRCPbs2YP29nbcfvvteOedd3D77bejvb0de/bsQSQSsbqIRLakCD6CERV06NAh3HfffYjH4wAg1VYURQEA+Hw+7N+/H11dXZaUkciuGC5EBRw6dAg7duyAEOKy2xqrqgpFUXDw4EEGDJEJw4UoRyQSQVNTExKJxGWDRaeqKrxeL8LhMKqrqxe/gEQlgH0uRDn27duHeDy+oGABAE3TEI/H8dxzzy1yyYhKB2suRCZCCLS3t6O3t/eKRoQpioLW1lZ0d3cb/TFEKxnDhchkYmICdXV11/T5YDBYxBIRlSY2ixGZRKPRa/r8zMxMkUpCVNqcVheAyErZbBbj4+MYHh7GyMgIzp07d03f76233kI0GkVzczM792lFY7jQipJKpTA2NoaRkRGMjIxgfHwcmqahrKwMDQ0N2L59O1paWnDhwoUr/t4NDQ2IRqM4cOAAAKCyshLNzc1oaWlBc3PzNTW3EZUa9rnQspZMJo0gGRkZMZZu8Xq9aGhoMP4LBAJGR/yePXvwyCOPXHGH/u7du7Fz504kEgkMDAzgwoUL6O/vx9DQEIQQ8Pl8aG5uNgInFApBVdkyTcsTw4WWlVgshtHRUSNM9OVZKioqpDCprKyc93sUe55LKpVCOBxGf38/+vv7EQ6HkclkUFZWhjVr1hg1m9WrV8PhcFztr05kKwwXKmnT09NSzSQWiwEAqqurpTApLy+/ou97pTP0X3rpJdx9990L+t7ZbBZDQ0Po7+83ajepVAoOhwNNTU1G7WbNmjUoKyu7onIT2QXDhUqGvjKxHiSjo6NIJBJQFAXBYNAIklAoBI/Hc80/b6Fri7344osLDpZChBAYHR01gubChQuIx+NQFAWNjY1Sv43X6722X4poiTBcyLY0TcP4+LgRJCMjI0in03A4HKirqzPCpL6+Hi6Xa1HKEIlE8Nxzz+GZZ56RRpK1tbVh586d+OpXv4qqqqqi/9yJiQmpZjM9PQ0AqKurM4KmpaUFFRUVRf/ZRMXAcCHbyGQy0kiusbExZLNZuFwu1NfXG2FSV1e35H0Teq1pZmYGFRUVqKmpWdKZ+NPT00bYXLhwAZOTkwCAQCAg1WxqamqWrExEl8NwIcvMzc1Jne+Tk5PQNA1ut1vqL6mpqeGoqhyzs7NSzWZkZAQA4Pf7jbBpaWlBXV0dl6MhSzBcaMnE43EpTKampgBc6rdobGw0+kuqq6t5Q7xCyWQSAwMDxoi0wcFBaJoGj8cjDX9ubGxkUNOSYLjQopmZmZFGculLq1RWVko1E/YbFF86ncbg4KBRuwmHw0in03C5XGhqajKa0ZqamuB0ci41FR/DhYpCCIFIJGJ0vg8PDxujrGpqaqSRXD6fz+LSrjyapmF4eFhqSksmk1BVFatXrzZqN83NzXC73VYXl5YBhgtdFU3TpGHBIyMjmJubg6qqqK2tRSgUMsKENyv7EUJgfHzcGCDQ39+PWCwGRVEQCoWkQQJXOkeICGC40ALpCzyaR3Lpw4LNI7nq6+vZzFKipqampJrNxYsXAQDBYNAYINDc3LwoQ69p+WG4UEHpdFrqfDcv8KjXShoaGlBbW8sO4mUqGo0aAwQuXLiA8fFxAEBVVZVUs6mtrbW4pGRHDBcCcGm0Ue6wYPMCj3qgLPX8DrKPRCIh1WyGh4eNBTn1oGlubkZDQwPPEWK4rFSzs7NSf4m+wKPf75dGcrEJhOaTSqWk4c/hcBjZbBZutxtr1qwxajerVq3igpwrEMNlhZienjZGcc23wGMoFILf77e4pFSqMpmMtCDnwMAAUqkUnE4nVq9eLQ1/5oKcyx/DZRkSQuDixYtSzURf4DF3WDAXQqTFomla3oKc+nm4atUqafgzz8Plh+GyDGiahomJCWm14FQqBVVVpQUeQ6HQoi3wSPRJhBB5C3LqE2vr6+ulQQKcWFv6GC4lyLzA4+joKEZHR5HNZuF0OhEKhRAKhdDY2GjJAo9EVyISiUgLcupLAgUCAWn150AgYHFJ6UoxXEpEOBw2+ksmJiaMBR7Nw4KDwSCHBVNJi8ViUs1mdHQUwKWdRPWg6ezsZM2mBDBcSsTJkyeRSqVQVVWFyspKVFVVwefzccgnLWvpdBrT09PGfzMzM9i0aRPn1pQAhkuJ0BcdJFrJstksFEVhDb0EMFyIiKjoGP9ERFR0XGHwY5qmGYsx0rWpra3lvIUSpWkawuEw5ubmrC5KyWtsbFzRk5IZLh/TNA3vvPOOsWc7XTkhBKLRKLZv346Wlhari0NXIZvN4uWXX0Ymk+Es+mswNTWFL33pS+js7LS6KJZhuJgoioJt27ahsbHR6qIUpGkaNE2z7ZL2mqbh0KFDVheDrpGiKLj77ruxdu1aq4tSMtLpNFRVhcPhgKZp+Ld/+zeri2Q5e96lLKSfIHYihEAsFsM777yDaDSKzs5ObNy4kSNmaNGoqgqn0wkhBIQQyGQycLlcHPqeI51O48iRI/jggw9QXl6OO+64A2vXruXfCQyXknHq1Cn09/cDAI4fP45QKIS6ujqLS0XLmb5cy6uvvorx8XF0dnbit37rt9hc9jEhBI4fP4433njDWIJpcHAQv/M7v2N10WyBj74lIJ1OY2RkxHidyWQwMjICjiKnxZTJZHDo0CGcPXsWExMTOHLkCE6ePMnz7mOzs7M4duwYNE0zjlVWVrI58WMMlxKQTCYRj8elY/oaTESLZWRkxKgtA5ee1E+cOMGRZLj0tzh//ryxDxJwqSnxjjvu4NI0H2O4lIBEIoFMJiMdi8Vi0hMTUbH19vbmnXdTU1MYHh62qET2IYRAd3e3UYtzuVz47d/+bWzevJn9LR9juJSARCKRFySJRALZbNaiEtFyp3fi+/1+6WaZzWZx/vz5Fd80lkwmEQ6Hjde1tbW49dZb4Xa7LSyVvbBD3+aEEJidnc07nkqlkE6n2blKi0JRFHz2s5/Fbbfdhv7+fvz85z83mmYHBgagaZrtRlUupcnJSczMzBivm5ubOT8uB2suJSC3vwW41NnK1QRoMamqioqKClx33XXS3K+JiYmC5+RKIYTA4OCg0WSoKAonDRfAcCkBhWou2WyWHau0JFRVxerVq43X8Xgck5OTFpbIeoODg8a/PR4P6uvr2deSg+Fic0KIgiGiaRrDhZbM6tWrjZunpmkYGhpasf0uqVTK2MQMAKqrqzlCrACGi83NFy4AGC60JBRFQW1tLTwej3FsJYdLNBqV+lsaGhpsuySTlRguNpfNZuftW0kmkyv2Aqel5ff7pX3sV/IK4uPj40ilUsbrVatWWVga+2K42JymafNexKy50FJxOp1oaGgwXkejUUxPT1tYImsIIaTVMRwOB0KhEPtbCmC42Fwmk8mbyKZjuNBSMj+hp9NpjI+Pr7iasx4uOp/Ph+rqausKZGMMF5vLZDLzzsRnuNBSURQFoVDI6FvQh+OuNHNzc9JIuUAgIPVF0f9juNhcOp2edyZ+KpVacU+OZJ1AICDtrDg8PLziliCKRqPS1IBQKLSiJ5NeDsPF5tLp9LwBcrmvERWbx+NBbW2t8XpqampFTabUtyAwd+ab+6FIxnCxubm5ucuGy0p7ciTrKIoizdSPx+MrbnVuc3+L0+nk5MnLYLjYXO5IMXMVnOFCS0lRFKxatcq4mWazWQwPD6+Y2nOhzvyqqioLS2RvDBcbE0JIVXBVVeHz+YzX2WyWKyPTkqqrq5NW/h0aGrKwNEtrbm5OqqkFAgF4vV4LS2RvDBebY7iQnfj9fmno7ejoqHSOLmfT09OIxWLG6/r6eqgqb6Hz4V/G5swXrqIoUrhomjbvHBiixeByuaRO7Onp6RUxmVIIgfHxcamZ2tz/RPkYLjZnDheHwyFVwxkuZIWmpibj3/oijiuh38W8Ayc78z8Zw8XmzE9KueEihFix6zuRNfQRY+aNsfr7+y0s0dLQNE3qzM9tHqR8DBcbyw0Pp9MpzQbO7fAnWgqBQACVlZXG68HBwWX/kJNIJKSZ+TU1NdzS+BMwXGysULjkntDL/aIm+3G73VJ/w8WLF5d9v0vuhNFVq1axM/8T8K9jY7krIuvhYj6pWXMhK6xdu9b499zcHMLh8LLtd9HXUdNHZiqKIm2eRoUxXGxM0zRpqLHL5YLL5ZJOaoYLLTX95lpWVmYc6+vrW9bhEg6Hjdf6tsZ0eQwXG8sdDaaHi7nmwvXFyAqBQADBYNB4HQ6HkUgkLCzR4kkmk9JIsWAwKC3gSYUxXGwsd5JkWVkZHA4Hm8XIci6XS2oai0aj0iZay4UQAqOjo9K2xmvWrOG2xgvAcLGx3L1cXC4XnE5n3vpiRFZoa2szHnQ0TUNPT4/FJVocfX19xkOeqqpSqNL8GC42ls1mpXApKyuDqqpSuHBPF7KCoihoaGiQFm7s7e1ddhvYpdNp9Pb2Gq/9fj8aGxvZmb8ADBcby+1P0cPFXCXnDH2yitfrxbp164zXU1NTy6ppTF/yZXx83Di2Zs0aaQkmmh/DxcZyw0XvzDeHC5fdJyt1dHQYTWPZbBbd3d0Wl6i4PvroI6PpWVEUdHZ2stayQAwXGytUc1EURVp6I7dfhmip6EOSzbP1z507t2wGmczNzeGjjz4yXldUVKClpYXhskAMFxvLXRFZn+PCcCG78Pl8Ugf3cmka0+e2mPdvaW1t5RDkK8BwsancpV/MoWIOl9yJlkRLrbOz02gay2Qyy6Zp7OzZs9IosY0bN1pcotLCcLGx3I3C9L4W88xobhhGVirUNNbT01PyTWOxWAznzp0zXgeDQTQ1NbFJ7AowXGxsIeHCPV3IauXl5WhpaTFeT01NYXR01MISXRshBM6fPy9NnOzo6OAqyFeI4WJj5mYx8/wWhgvZjXkUVSaTwblz50q230UIgbNnzxrld7lcuO666ywuVelhuNhY7i6UDocjr0NfCMFwIUvpTWMVFRXGsd7e3pI9L6PRKAYGBozXoVCIu05eBYaLTRXay0XvNM0Nl1Jv36bSV15ejtWrVxuvx8fHpZFWpUIIgQsXLkh7t3R0dHAtsavAcLGphYYLwPXFyHqKomD9+vXG61QqhfPnz5dc05gQAr/5zW+McpeVlaGtrY21lqvAcLGpQsvtm8OFy+6TnSiKgubmZni9XuPYuXPnSm4O1szMjNQkVl9fj9raWgtLVLoYLjZVKFzM/+aGYWQ3VVVVCIVCxuvh4WFpxJXd6U1is7OzxrH29nY2iV0lhotN5YaLvvQLIDeRAQwXsgeHw4G2tjbjdTwex8DAQMnUqoUQ+Oijj6RRYuvXr2eT2FViuNhUob1cdLkbhs3X5yKEwMzMDIaHh5FIJErmIqfSpCgK1q1bZ5yrQgj09PSUzHnHJrHiYn3PAtlsFkNDQ4hEImhsbEQwGMx7OspkMnm7UOpyNwwrVHPRq/hHjhxBMplEVVUV7rrrLtTU1PBJjBZNbW0tgsEgRkZGAAADAwNIJBIoLy+3uGSXp0+cNDeJcZTYtWHNpUiEENA0DUKIyz6paZqG06dP4/Dhw/jVr36Fl19+uWDTQW7Nxe12G6HgcDikk77QhmGJRALHjh0zaiyRSATHjx/nUjG0qFwuF1pbW43X0WgUg4ODtq+9aJomTZwsKytjk9g1YrjMwxwWnySTyeD999/Hyy+/jCNHjmB2dnbez0WjUbz//vtGcMzNzeHkyZN5TVuFltvXFdrTxfxeIQQGBwcRjUal7zk8PIxIJDLv7zs3N4fe3l709PQsuBktk8lwtBpJ1q9fb9SsNU3Dhx9+aHGJPtnU1BT6+/uN1w0NDairq7OwRKWP4VKAEALDw8N48803EY/HL3vjFELggw8+wLFjxzA8PIwPP/wQb7zxRsF+EL2pKncr2KmpKUxMTEjHcmsj5nDJnaWf+14hhNR2rMtkMhgeHi74+6RSKbz66qv45S9/iddeew3/+7//i0QicdnfO5PJ4N1338Wrr76KZDI573tp5VAUxWjq1fX19UnNTXajL/diPoc3btwoNT3TlWO4FDA9PY233noL3d3deOWVVxCLxeYNmHg8LlWngUs1hHA4nPcZTdMQDofzvoemaXk3/dy9XC4XLrn9M+l0Oi+sdKOjo3nl0ieOmcs2Ojqa93vlfub06dM4c+YM+vv78dprr0mzmmnl0puUdNFoFH19fbat3SYSCbz//vvGa7/fj/b2djaJXSOGS450Oo23337baFIaGxvD8ePHC04GE0JgZGQk76lMCIHe3t6C/SDzNUuNjY1J7zfXbnLDJDdscsNlZmZm3lrHxYsXCzbB9fT05L23r69v3sECExMTOH36tFHm4eHhgrUlWpmuu+46adSYuSnYToQQ6O7uxuTkpHFs/fr1qKqqsrBUywPDJYfD4UBra6t0875w4QLGx8cLPnkNDQ0V/D7j4+N5N/jp6em8JrFCX9P7P3SqqkrlASAt/53NZo05MXrnvXmOjPm9iUQir1wXL17E9PR0XplmZmZw8eLFvONCCJw5c0YKqeuuu056WqWVS1EUhEIhNDY2GscGBgakG7hdpNNpHD9+XJrbctNNN7HWUgQMlxyqqqKjowNbtmyRlhDv7u4uOKJrvuanRCKBixcvGp8RQmBqasp4rSgKampqjPcnk0mpBmQOl9zRYYAcGLkTLs2B4HA4pL02MpmM1NGv174KjSLTNC2vRgVcCh1zE1ogEMDNN9/MNmoyOJ1ObNq0yXidTCbx/vvv26ppTJ+HY35AXLt2LVatWsVwKQKGSwGKoqC1tVWqGg8ODuZ1WsfjccRiMeO13+83TkohBMbHx6X3m1eJdblc0k0/m80iGo0aF585XJxOpxQuiqLkhYvefKXXXHRut1u6WIQQmJ6elkLPvLGTqqrSBM3cGps+WMBcvvb2dm6kRBJFUdDR0SEtw3/mzBlb9culUim88847RnOdw+HA1q1b+ZBUJAyXeZSVlUk3/3g8ntcZPj09LTUNrVu3Dh6Px3g9MTFhvF/TNGmdJZ/Ph1AoJN3I9VDIbRbLnTQJ5Ndc9HDJZrNS4Pl8PgQCAanPxlxzSafTUhjV1NSgurpaKpO5VqRpmtS34na70dzczCc9ylNZWYnOzk7j9cWLF6UVh62kD2LJrbW0tLTwXC4Shss8FEXBmjVrjJu6PndEJ4SQmr1UVUVjYyP8fr/xnmg0atyY0+m01Ozl9/tRWVkp9aXoN31zWACXgs4cQgDyagp6GKVSKalPpaKiAuXl5dL7zTWk2dlZ6f25S17E43Hp67Ozs1LbeV1dnfR0SmR24403Gue4EKLgnC4rzM3N4ejRo0atxel0Ytu2bZyRX0QMl8sIBALSjXN0dFS6MMx9Gy6XC9XV1QgEAsaxeDxuNKUlk0mpNlJZWQmPxyMtUT4zMwNN06BpmvRzcpfYBy6Fi/nY3NwchBBIJBJSMFVWVsLpdErLb8zOzhqhl1szCQaD0u+QTqeNGpcQAmNjY9LvsXr1aj7pUUGKoqChoQFr1641jg0PD6O/v9/S2ou+QOXw8LBxbN26dVi7di3P5SJiuFyGy+WSlhCPxWJG7SKbzUojrHw+Hzwez7w35lgsJt3Eq6qq4HA4pJt+PB5HOp1GNpuVwsXj8eSd9LmBo4fY7Oys1DlfWVkJVVWlGpUedHrtS6eqqhGQ5j4acw3NfEE6nU40NjbygqR5qaqKm2++2ThXs9ksTpw4Yemw5GQyiaNHj0ojxG6//Xb2tRQZw+UTmIdTZjIZY/RUKpWSmrkqKyvhcDhQXV0t3ZgjkQiEEFJTlKIoqKysNP5XNzc3h2QyiXQ6LQWRuR9H53K5pItBDxdzf4qqqkbNy/xz0um0tOaYzu12o7y8HBUVFVIfjT4QIZ1OS4MU/H4/m8ToshRFwdq1a9HQ0GAc6+vrKziZdynow+jNg1ja2tqwZs0aPiQVGcPlMhRFQTAYlPpFxsbGAFyqIZibn/Snfb/fL7Xb6jdvcy3H6XQaNQnzzTmTySAejyOVSkm1j/nCxfxzksmkscS+zuFwwOfz5YWYpmmIxWLGCDWdz+eD2+2Gx+OBz+czjk9PTyObzWJ2dlYaLFBbW5u35TJRrrKyMtx8883GzXtubg7vvfeeJeGSSCTwq1/9Sqq1bNu2jbWWRcBw+QT6k7xucnIS6XTauOHq9BFWXq9XCgP9feabuH4Dz73p6+GQTCalZgOv15v3VOVwOKTQ0z9jDhf95wCXahnmZrRoNIpUKiUNDdVrX06nUyqXHqQTExNSc525yZBoPoqioLOzUxqFePbsWam5dSnoKwWY56Z1dHSgqamJtZZFwHD5BE6nUxo9FYvF8mauOxwOVFVVGcu0mMMoFoshmUxKT/zl5eXGE395eblUA4lGo3krEptrEeafaR4BNjc3l9dU5/P5jO/t8/nyhiPHYjGp9qVf/IqiSH1Hc3NzmJ2dlZoSnE4n6urqeFHSgpSXl+OGG24wXs/OzuK9995b0jLMzs7i2LFjxrXldruxbdu2vMEyVBz8qy6A+Qldn5VvDhe9rwK4dGM2T75MJpOIRCJ5w4P1E9rr9ebd9M0BoapqwWYxRVGkkWb6EGTzzzHXVtxutxRGsVgM09PTUg3J3F9kfsrMZrOYnJyUnvjY30JXQlEU3HDDDdI5c+rUKWnVisUkhMB7770nDaPv7OzkgJRFxHD5BHq/i7l2MTg4KPWhVFRUSE1U5htzJpPJG8JsDp+ysjKpZhKLxaRaTm4NxcwcLnpTnfnnmJu2Cg1HNnfOOxwO4/2KoqC6ulpqh+7v75ea9oLBIPtb6IpUV1fj+uuvN17PzMxIo7YWi3mzPJ3H42GtZZHxL7sAfr9fGsobDoelAKiurjZOUr3mYh4xNjAwINUQzOGiqmrecGRzcLlcrrxFK3XmUMpms5iYmMgbhqxTFEV6akwkElIzl9vtlr5feXm5VGMaHByUgss8+odoIRRFwa233iqdl6dOncK5c+cWNWCEEDh69Kh0XW3evBmhUIi1lkXEcFkAl8sl9buk02kpLMwbIwHIGzFmroo7HA5UVFQYJ3Vup34ymZTWIHO73QVnDSuKIoWSpml5a4SZ1zoD5FBLp9NS0155ebkUYm63W3q/+fdlfwtdrerqatx2223GuZNKpfDyyy9f09BkfdfYQjvH6g93v/71r41jFRUV2LZtG8/fRca1DhaosbGx4J4nDocDwWBQOlG9Xi+8Xm/BZS5yawgA8m7i5k52n8837zBJfZixfkGZ+0Rym8Fy+4Jy6ZM6ze+vra0tuKVAZWUl+1voqiiKgltuuQXd3d24cOECgEsrXfz0pz/FHXfcgYaGBmQyGcRiMczOzsLpdCIUCqGhoQFOp1O6zvTdUE+fPm1spbx+/Xps3rzZaDKenZ3F4cOHjVUlFEXB1q1bUVNTw3BZZAyXBdD3p3C73Xn7sfh8PqnmAVy6sVdVVeXtYQ/k1xD0mouqqgVnLZeXl897EeijwfQQM3/e4/Hk9dVUVFTA6XRKEzR15uX/dXqzQe7T4KpVq9jfQlfN7Xajq6sLL7zwgtFUFYlE8POf/xyqqkIIIZ1zTqcTa9euxRe+8AXU1tYa10M2m8Xhw4elzfx6enpw6tQp3HXXXQgEAnjllVekNQEbGxtxyy23MFiWAJvFFsjv96Ouri7veH19fV6fiD4IoJBAIJBXE8ltRjPLDS4zt9s9b3+M3+/PCwB9iZpcennNF5yiKKirq5P6moD/v9CJrpb+sPalL31JGvwCoGDTViaTQU9PD/7rv/7LGF2mz1kptJTM0NAQXnjhBfzrv/4rzp49axz3eDz4/Oc/Lw2EocXDcFkgVVXR2dkpjS5RVRWtra0Fn4Lm65MoFFC5M+LNzIMDcrlcrnk/Z14fTFdWVpZ3Mes/v1CTmcfjkZZMB4CWlhbp6ZHoaiiKgpaWFjzwwAPYsGFDwYek3JFc4+Pj+MUvfoFUKoWZmRm8+eabxgAWRVGkB7RsNivtv+R0OvHZz36WS+ovITaLLZC+BP+mTZuM9t3Ozs6CqwLru0x6PB5p3sl8HeFOpxOBQEBa50s/frl+En3yZu6mZEDhZi69NmLeRRK4FETz1Wg2bNiAVCqFoaEh1NbWSosQEl0LvV/vy1/+MiYmJjA8PGz0s1RUVMDn8yESieCXv/ylsfJET08P3n33XcTjcWngy3XXXYdbbrkFr7/+OsLhsFT7KS8vx2c+8xlpCRpafAyXK6CqKrZs2YKOjg4AMBafLMTn86GxsRG9vb3GsWAwWDAsFEVBfX09+vr6pON+v1/qlC+kUIg4HI55OyxXr16NU6dOSUOWm5ub5w0Ml8uFLVu2QNM0abg1UTEoigKHw4FQKFRwOSEhBMrKynDgwAFkMhkIIfD666/nrWCxfft21NfXo6mpCT09Pejr60MymURdXR02btzI2rYFGC5XQFEUY4LhQt57ww03YHx8HDMzM3C73bjxxhsLjvzS973I7WwPhUKX7TjXayK5gwHmmz2v9600NzcbQRYMBrFu3bp5Lzz9OBf2Iyvo65LddNNNePfddwEgbxTmli1bUF9fb2z/vXHjRmzcuDHv+9DSYrgsEr1p7Itf/CImJydRXV0tLa+SKxAIYM2aNcZNv6ysDO3t7Z/4c/Rtic1NBE1NTfN29DscDnzqU59CKBRCOp1GW1sbOzjJ1lRVxfbt2zE8PCyN/AIu1cS3bt2aNxiFrMdwWUT6MOPLjfjSqaqKbdu2wev1IhaLob29fUETFV0uF2666Sa89dZbmJubQyAQyHtqyy2Tx+PBpk2brvj3IbKCPmH43nvvxf/8z/9gcHAQQgg0NDTgnnvumXdQC1mL4WITiqLA5/Ph9ttvl44t5HNr165FdXU1ZmZmUFtbW3CJfqJSpnf+P/DAAxgaGoKmaVi1ahXPdRtjuJgIIfIWmSw1hUaOLRUhRN4kUyo9QgiEw2HbXwcDAwNWF6EgTdOkUaIrFcPFJBAIIBwO5w3VpYUrKyubt7+H7E8fudjT01NwuSNaGPNGfSuVIqzYa9SG9Fm//HNcO1VV2VRRovRFIOnarfTrgOFCRERFx6nWRERUdAwXIiIqOoYLEREVHcOFiIiKjuFSIrLZLGKxmLTgJNFKk81mMTMzw+ugBDBcSkQkEsELL7yQtyw/0UoyNjaGH/zgBxgbG7O6KPQJGC5ERFR0DBciIio6hgsRERUdw4WIiIqO4UJEREXHcCEioqJjuBARUdExXIiIqOgYLkREVHQMFyIiKjqGCxERFR3DhYiIio7hQkRERcdwISKiomO4lAAhBCYnJzExMYHJyUkIIawuEtGS06+DSCTC66AEKIL/D9lWJBLBvn37sHfvXpw7d8443tbWhoceeghf/epXUV1dbV0BiZYAr4PSxHCxqUOHDuG+++5DPB4HAOkpTVEUAIDP58P+/fvR1dVlSRmJFhuvg9LFcLGhQ4cOYceOHRBCQNO0ed+nqioURcHBgwd5YdGyw+ugtDFcbCYSiaCpqQmJROKyF5ROVVV4vV6Ew2E2DdCyweug9LFD32b27duHeDy+oAsKADRNQzwex3PPPbfIJSNaOrwOSh9rLjYihEB7ezt6e3uvaCSMoihobW1Fd3e30Q5NVKp4HSwPDBcbmZiYQF1d3TV9PhgMFrFEREuP18HywGYxG4nFYtf0+ZmZmSKVhMg6vA6WB4aLjfj9/mv6fEVFRZFKQmQdXgfLA8PFRoLBINra2q64vVhRFLS1taGmpmaRSka0dHgdLA8MFxtRFAUPPfTQVX12586d7MSkZYHXwfLADn2b4fh+Il4HywFrLjZTXV2N/fv3Q1EUqOrl/+/RZya/+OKLvKBoWeF1UPoYLjbU1dWFgwcPwuv1QlGUvGq+fszr9eKll17C3XffbVFJiRYPr4PSxnCxqa6uLoTDYezevRutra3S11pbW7F7924MDg7ygqJljddB6WKfSwkQQmBqagozMzOoqKhATU0NOy1pxeF1UFoYLkREVHRsFiMioqJjuBARUdExXIiIqOgYLkREVHQMFyIiKjqGCxERFR3DhYiIio7hQkRERcdwISKiomO4EBFR0TFciIio6BguRERUdAwXIiIqOoYLEREV3f8BmgHypOq/5yMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 500x400 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model = model.prune()\n",
    "model(dataset['train_input'])\n",
    "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 , r2\n",
      "gaussian , 0.7090268761989152\n",
      "1/x^2 , 0.21051195154680438\n",
      "sin , 0.1822506022370818\n",
      "abs , 0.12418544555819415\n",
      "tan , 0.10407480103502795\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "('gaussian',\n",
       " (<function kan.utils.<lambda>(x)>, <function kan.utils.<lambda>(x)>),\n",
       " 0.7090268761989152)"
      ]
     },
     "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', '1/x', '1/x^2', '1/x^3', '1/x^4', 'sqrt', '1/sqrt(x)', 'exp', 'log', 'abs', 'sin', 'tan', 'tanh', 'sigmoid', 'sgn', 'arcsin', 'arctan', 'arctanh', '0', 'gaussian', 'cosh'])"
      ]
     },
     "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)"
   ]
  },
  {
   "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 , r2\n",
      "gaussian , 0.7090268761989152\n",
      "J0 , 0.2681378679614782\n",
      "1/x^2 , 0.21051195154680438\n",
      "sin , 0.1822506022370818\n",
      "abs , 0.12418544555819415\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "('gaussian',\n",
       " (<function kan.utils.<lambda>(x)>, <function kan.utils.<lambda>(x)>),\n",
       " 0.7090268761989152)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# J0 shows up but not top 1, why?\n",
    "\n",
    "model.suggest_symbolic(0,0,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e78f4674",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "function , r2\n",
      "J0 , 0.9717763100936939\n",
      "gaussian , 0.7494106253678943\n",
      "sin , 0.49679878395526067\n",
      "1/x^2 , 0.21051195158162733\n",
      "abs , 0.12435207425739554\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "('J0',\n",
       " (<function torch._C._special.special_bessel_j0>, J0),\n",
       " 0.9717763100936939)"
      ]
     },
     "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.67e-02 | test loss: 1.80e-02 | reg: 2.87e+00 : 100%|██| 20/20 [00:08<00:00,  2.25it/s]\n"
     ]
    }
   ],
   "source": [
    "model.train(dataset, opt=\"LBFGS\", steps=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4773e989",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAFICAYAAACcDrP3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAx30lEQVR4nO3dfXAU930/8PfuPZ90kk6n0wMSj+KZGIyJa+rE1NgdSEz6YDtJ3TRuZpymbusxk8ZtY0/HdR9+mTjxtIXYrWmSZgpxa2canGls8EDTlLjjB/AzBowtEBgkoWeddNI9735/f+Dd7PfuJASsdHu692uGMbd6+grv3nu/T59VhBACRERENlJL3QAiIpp7GC5ERGQ7hgsREdmO4UJERLZjuBARke0YLkREZDuGCxER2Y7hQkREtmO4EBGR7RguRERkO4YLERHZjuFCRES2Y7gQEZHtGC5ERGQ7hgsREdnOXeoGEJUDIQSGhoYwPj6O6upqRCIRKIpS6mYRORZ7LkRTiMVi2LlzJ5YtW4ZoNIrFixcjGo1i2bJl2LlzJ2KxWKmbSORICp9ESVTcgQMHcOeddyKRSAC42HsxGL2WYDCIvXv3YuvWrSVpI5FTMVyIijhw4AC2bdsGIQR0XZ/081RVhaIo2LdvHwOGyILhQpQnFouhra0NyWRyymAxqKqKQCCArq4u1NXVzXwDicoA51yI8uzevRuJRGJawQIAuq4jkUhgz549M9wyovLBnguRhRACy5YtQ2dnJy7n0lAUBUuWLEFHRwdXkRGB4UIkGRwcRDQavaqvj0QiNraIqDxxWIzIYnx8/Kq+Ph6P29QSovLGcCGyqK6uvqqvD4VCNrWEqLwxXIgsIpEI2tvbL3veRFEUtLe3o76+foZaRlReGC5EFoqi4P7777+ir92+fTsn84k+wgl9ojzc50J09dhzIcpTV1eHvXv3QlEUqOrUl4ixQ//ZZ59lsBBZMFyIiti6dSv27duHQCAARVEKhruMY4FAAPv378eWLVtK1FIiZ2K4EE1i69at6Orqwo4dO7BkyRLpY0uWLMGOHTvQ3d3NYCEqgnMuRNMghMCLL76Iz3zmM3j++eexadMmTt4TTYE9F6JpUBQFdXV10n+JaHIMFyIish3DhYiIbMdwISIi2zFciIjIdgwXIiKyHcOFiIhsx3AhIiLbMVyIiMh2DBciIrIdw4WIiGzHcCEiItsxXIiIyHYMFyIish1L7hNNk6ZpGB8fR3V1NVwuV6mbQ+RoDBeiaRJCQNd189HGRDQ5hgsREdmOcy5ERGQ7d6kbQGTQNA1HjhxBPB4vdVPK3rXXXovGxsZSN4MqGIfFyDGSySS2bt2KVCqF6urqUjenLAkhcPr0aezatQu33XZbqZtDFYw9F3IMIQRUVcU3v/lNfPKTnyx1c8pCIpHASy+9hA0bNqC+vh7ZbBaf/exnwXtGKjWGCzmOy+WCz+crdTMcTwiBQ4cO4ctf/jLmzZuH22+/HXfffTdXspEjcEKfqExpmoann34aqVQKnZ2deOKJJ9DR0VHqZhEBYLgQla3Tp0/j0KFD5us1a9bg4x//eOkaRGTBcCEqQ0II/OQnP8HIyAgAQFEUfP7zn+dCCHIMhgtRGRoZGcGzzz5rvm5qauLqMHIUhgtRmRFC4Gc/+xlOnTplHtu6dSva2tpK2CoiGcOFqMykUik89dRT0DQNABAMBvGFL3wBqsrLmZyDZyNRGRFC4I033sCRI0fMYxs3bsS1117LJcjkKAwXojKiaRqeeuopJJNJAIDb7cbv//7vc18QOQ7DhaiMdHZ24uDBg+br5cuX4+abb2avhRyH4UJUJnRdxzPPPIOhoSEAF5cf/87v/A7q6upK2zCiIhguRGWit7cXP/7xj83Xzc3NuP3229lrIUdiuBCVAWPT5Llz58xjt99+O+bPn1/CVhFNjuFCVAaGh4exZ88es9pxXV0dvvjFL7LXQo7FcCFyOCEEnnvuOXzwwQfmsU9/+tNYsWIFw4Uci+FC5HCjo6P4wQ9+AF3XAQDV1dX48pe/DJfLVeKWEU2O4ULkYEII7N+/H8eOHTOP3XrrrVi3bh17LeRoDBciBxsbG8P3vvc9qdTLV77yFXg8nhK3jGhqDBcihxJCYN++fXjnnXfMY7fccgt+5Vd+hb0WcjyGC5FDxWIx7Nq1y+y1BAIB/OEf/iG8Xm+JW0Z0aQwXIgcy9rXkz7Vs3LiRvRYqCwwXIgcaHBzEv/zLv5i9lurqatx3333stVDZYLgQOYwQAk899RTef/9989inPvUpXH/99ey1UNlguBA5iBACZ8+exfe+9z1zN35tbS3uu+8+uN3uEreOaPoYLkQOous6/vmf/xnd3d3msc9//vPc10Jlh+FC5BBCCLz99tv40Y9+ZB5rbW3Fn/zJn/ARxlR2eMYSOUQmk8GOHTswOjoK4OLzWr7yla9g8eLF7LVQ2WG4EDmAEAI///nPpadMrl69GnfffTeDhcoSw4XIAWKxGL797W8jlUoBANxuN7Zv345IJFLilhFdGYYLUYnpuo5/+7d/w9tvv20e++QnP4nf+I3fYK+FyhbDhaiEhBA4efIknnzySbOkfk1NDf7iL/4CwWCwxK0junIMF6ISSqfTePTRR9HX12ce+93f/V3ccMMN7LVQWWO4EJWIEAL/9V//hRdeeME81t7eju3bt/NBYFT2GC5EJSCEwLlz5/Doo48ik8kAALxeL/78z/8cra2t7LVQ2WO4EJVANpvFt7/9bXR2dprHPv3pT+O3f/u3GSw0JzBciGaZEALPP/88fvzjH5vHWlpa8NBDD8Hv95ewZUT2YbgQzSKjMOXf/u3fSntaHnjgAaxcuZK9FpozGC5EsyiVSuHv/u7vcObMGfPYli1b8IUvfIHBQnMKw4Volui6jqeffho//elPzWNtbW145JFHuKeF5hyGC9EsMCoef/Ob30Q2mwVwcXXYgw8+iBUrVrDXQnMOw4VohgkhMDw8jIceegj9/f3m8TvuuAOf+9znGCw0JzFciGZYNpvFt771LRw+fNg8tmbNGjz88MPw+XwlbBnRzGG4EM0gIQT+8z//E7t375YeW/yNb3yDmyVpTmO4EM0QIQTeeOMN/M3f/I257NjlcuGrX/0qNm3axGChOY3hQjQDhBDo6enBAw88IBWl/MxnPoN7772XtcNozmO4ENlMCIHx8XE89NBD0jNaPvaxj+H//b//x2XHVBEYLkQ2y2QyePTRR/Hcc8+ZxxoaGvDYY49h/vz5HA6jisBwIbJRLpfDrl278N3vftd8+FcgEMBf//VfY+PGjQwWqhgMFyKbaJqGH/7wh1IZfZfLhfvuuw933XUXVJWXG1UOnu1ENtA0DU8//TQefvhhJBIJAICiKLjrrrvwta99DR6Pp8QtJJpdDBeiq6RpGp555hk8+OCDiMfjAC4Gy7Zt2/CNb3yDE/hUkRguRFdB0zT86Ec/koIFAG699Vb84z/+I+rq6jjPQhWJ4UJ0hTRNw3/8x3/g61//OsbGxszjmzdvxj/90z8hGo0yWKhiMVyIrkA2m8UPfvADPPjgg1Kw3HzzzXjyySfR1NTEYKGK5i51A4jKiRACiUQCO3bswHe+8x2zrAtwscfy5JNPorm5mcFCFY/hQjRNQgh0d3fjL//yL/Hcc89B0zQAFyfvf/3Xfx1PPPEEeyxEH2G4EE2Drus4fPgw/uzP/gzHjh0zj6uqit/8zd/E3//93yMSiTBYiD7COReiS8hms3jmmWfwxS9+UQoWn8+HP/7jP8bjjz/OYCHKw54L0SSEEBgbG8OOHTvw5JNPIplMmh+LRqP4q7/6K9x1113weDwMFqI8DBeiInRdx9GjR/HII4/gxRdfNOuEAcDatWvxD//wD9iwYQNLuhBNguFCZCGEQDwex7/+67/i8ccfx9DQkPkxVVVx22234Vvf+hafIkl0CQwXoo/oum4+OfKll16SeivBYBD33nsvvva1ryEUCjFYiC6B4UIVTwiBwcFB7Nq1C9///vcRi8Wkj69cuRKPPPIItmzZApfLxWAhmgaGC1UsIQRSqRReeOEFPPbYY3jvvfcghDA/HgwG8Xu/93t44IEHuDGS6DIxXKjiCCGQzWZx+PBh7NixA7/4xS+QzWbNjyuKgo997GN4+OGHceutt7K3QnQFGC5UMYQQ0DQNb731Fp544gkcPHjQfPaKoba2Fvfccw/uu+8+NDQ0MFSIrhDDheY8I1Teeecd7Nq1Cy+88IJUHh8APB4PNm/ejK9//etYv349XC5XiVpLNDcwXGjOEkIgk8ngrbfewne/+10cPHiwIFRUVcXatWvx1a9+FZ/61Kfg9/vZWyGyAcOF5hRjQn5oaAiHDh3Cv//7v+PVV18tGP5SVRXLly/HH/3RH+H2229HbW0tQ4XIRgwXmhOEEEgmkzh27Bh+8pOfYP/+/Th37py0VwW4OFm/dOlS/MEf/AE+97nPob6+nqFCNAMYLlS2hBBIp9Po6OjAwYMH8fzzz+PEiRPSM1YMqqpixYoVuOeee3DHHXew0CTRDGO4UNkQQpjlWT744AP87//+L/77v/8bJ06cwPj4eNGv8fv9WL9+Pe6++27cdtttfKY90SxhuJCjGU9+PH/+PN5880383//9H44cOYLz588jnU4X/RpVVdHW1oatW7fis5/9LNatW8eJeqJZxnAhR9J1HZ2dnfjhD3+In/3sZzh79iwmJiakHfRWiqKgvr4eN954I+644w7cdNNN5tAXQ4Vo9jFcyHFyuRwee+wxfP/730d/f/+kn6eqKiKRCDZs2IBt27Zh06ZNaGtr4456IgdguJDjqKqK7u7uosESDAaxYMEC3HDDDdi8eTOuv/56NDc3M1CIHIbhQo6jqiruvfdePPfccxgeHkZzczM+8YlP4KabbsLHP/5xLFq0CNXV1QDAQCFyKIYLOdLKlStxzz33IJvN4p577sGCBQs4f0JURhgu5Ci6ruO1115DIpHAunXr4HK5cOLECZw4caLUTSsLuVwOIyMjpW4GEcOFnENVVaxatQoHDx7EwYMHS92cshUKhVBXV1fqZlCFU8RkazuJZpkQArlcbtLlxjR9brcbqqqWuhlUwRguRERkO97aEBGR7RguRERkO4YLERHZjuFCRES2Y7gQTZMQApqmcTUb0TQwXIim6ejRowiHwzh69Gipm0LkeAwXIiKyHcOFiIhsx3AhIiLbMVyIiMh2DBciIrIdw4WIiGzHcCEiItsxXIiIyHYMFyIish3DhYiIbMdwISIi2zFciIjIdgwXIiKyHcOFiIhsx3AhmgYhBEZGRqT/EtHkGC5EU4jFYti5cyeWLVuGzZs3Y3x8HJs3b8ayZcuwc+dOxGKxUjeRyJEUwVswoqIOHDiAO++8E4lEAgCk3oqiKACAYDCIvXv3YuvWrSVpI5FTMVyIijhw4AC2bdsGIQR0XZ/081RVhaIo2LdvHwOGyILhQpQnFouhra0NyWRyymAxqKqKQCCArq4u1NXVzXwDicoA51yI8uzevRuJRGJawQIAuq4jkUhgz549M9wyovLBnguRhRACy5YtQ2dn52WtCFMUBUuWLEFHR4c5H0NUyRguRBaDg4OIRqNX9fWRSMTGFhGVJw6LEVmMjY1d1dfH43GbWkJU3tylbgBRKWmahqGhIfT396O/vx+dnZ1X9f3effddZDIZtLS0IBQK2dRKovLDcKGKks1mMTAwgIGBAfT392NoaAi6rsPr9SIajeLGG2/EokWLcPbs2cv+3vPmzUMqlcKhQ4cAAFVVVWhpaUFLSwuam5sRDoft/WWIHIxzLjSnpdNps1cyMDBglm7x+/1obGw0/9TW1poT8Tt37sSf/umfXvaE/o4dO7B9+3ak02n09vbiwoUL6O3txcDAgPkzm5ubzcCpr6+HqnJkmuYmhgvNKRMTE2avpL+/35xDqa6uRmNjI6LRKBobG6ccsrJ7n0s2m0V/f78ZNn19fdA0DR6PB01NTWbYRKNRuFyuK/7diZyE4UJlbWxsTOqZTExMAABqa2vNXkk0GkUwGLys73u5O/T379+PLVu2TOt7a5qGgYEBqXeTzWbhcrnQ2NhoDqM1NTXB4/FcVruJnILhQmVDCIFYLCaFSSqVgqIoCIfDUpj4fL6r/nnTrS327LPPTjtYihFCYGhoyAybCxcumL9XQ0ODGTbNzc3w+/1X90sRzRKGCzmWrusYHh6WwiSbzUJVVTQ0NJhDXA0NDTN2hx+LxbBnzx585zvfwenTp83j7e3t2L59O770pS+htrZ2Rn6u0au5cOECxsfHAQDhcNgMm5aWFlRVVdn+s4nswHAhx8jlchgcHDSDZHBwEJqmwe12m0ESjUYRiURmfW5CCIHh4WHE43GEQiHU19fP6k788fFxs1fT29trlvqvqamRFgnU1NTMWpuIpsJwoZLJZDLS5Pvw8DCEEPD5fGaYNDY2IhwOs6RKnmQyKc3ZDA4OArg4TGcNG/7bUakwXGjWJJNJKUyMu+9gMCiFSU1NDd8QL1Mmk0Fvb68ZOAMDA9B1HT6fz5yvaWlpQUNDA5c/06xguNCMGR8fN4e4+vv7zdIooVBIWhZcXV1d4pbOPblcDv39/WbY9PX1IZfLwe12m8ufm5ub0djYCLebe6nJfgwXsoUQomBZsLHKqq6uTlrJFQgEStzayqPrOgYHB6UVaZlMBqqqIhqNSivSvF5vqZtLcwDDha6IEAIjIyNSmKTTaSiKgkgkIk3A883KeYz/f9ZFAolEwvz/ZwyjNTc382aArgjDhabFKPBoDHENDAwgl8vB5XIVLAvmMEt5Ghsbk8LGqG5QV1cnLRLgMCZNB8OFispms+ayYGuBR4/HI02+sz7W3DUxMSHttRkZGQFwsZSOda8NH+1MxTBcCMAvCzwaPZP8Ao9GoNTV1XElV4VKpVLSirTBwUHzHLFWf45EIjxHiOFSqRKJhDRfMjo6CuBimXjr5Ds35dFkstks+vr6zN5Nf38/NE2D1+uVVqSxIGdlYrhUiHg8boZJf3+/WeCxpqZGChOWE6ErZRTkNOZt+vr6CgpytrS0oLGxkQU5KwDDZQ4yCjxaNyzmF3iMRqOIRqMshEgzxqgNZ10kYJyH0WhUWpFmR6FRchaGyxwwVYHHSCRi9kxmssAj0aUYNz3WvTZGD7q+vl5akXa5j0gg52G4lKFcLic9991a4LGhocEMk1IUeCS6HPF4XFqRZsz91dTUSCvSOPdXfhguZaKnp0daFiyEgNfrlVZyhcNhLgumspZIJKSCnENDQwAu1p8zejULFy7k3GAZYLiUiXfeeQeZTAY1NTXmn2AwyCWfNKdls1mMjY1hdHQUo6OjGB8fx+rVqxGJRErdNLoEhkuZMIoOElUyTdOgKAp76GWA4UJERLZj/BMRke04zvIRXdcxNDSEbDZb6qaUvXA4zEq6ZUrXdfT39yOTyZS6KWWvoaGhopdUM1w+ous63njjDXNJL12ZeDyOjRs3oq2trdRNoSug6zpefvllaJrGPVFXYXR0FDfffDMWLlxY6qaUDN9FLRRFwXXXXYfGxsZSN6UoIYSjw0/XdfziF78odTPoKimKgo0bN2LevHnQdR3JZBLJZJIFKadJ13W88MILqPTpbGe+S5WQoiiO23gohEAqlcKbb76JWCyGBQsWYPXq1Y5rJ80dqqri1KlTOHXqFIaHh+HxePBbv/VbHO4sIpPJYHBwEJFIxHwwHkOY4VI2Tp48iXPnzgEATpw4gfr6esybN48nMc2YWCyGnp4eABffQGOxGMMljxACPT09+PnPf47q6mosWrQIy5cvL3WzHIGrxcpALpfDhQsXzNe6rpsXPdFMaW5uNm9edF1HX19fxQ/15BNCoLOzE5qmYXR0FMeOHTPrpVU6hksZSCaTSCQS0rGRkRHoul6iFlEliEQiUtXs3t5ehkueRCIh3fiFw2FWD/gIw6UMTExMIJfLSceSySSXTdOMCgQC0iOMh4eHkUqlStcghxFCoL+/X7rxW7hwIVfZfYThUgYmJiYK7hgzmQzS6XSJWkSVQFVVNDc3m6+TySRGRkZK2CLn6erqMv/ucrm4BN+C4eJwQoiiY7iapvEukmaUoigF8y4XLlzg0NhHjMc8G2pqaqSeXqVjuJSB/PkW4GLoJJPJErSGKkl9fb20Qqy3t5dzfR8xqjQbmpqaOCRmwXBxuKlCJJlM8i6SZpTf75cmqEdGRore7FQaY77FOhfa0tJSwhY5D8PF4XRdn7TOE4fFaKYpiiK9aabTafNhdZVMCCGtEvN6vWhoaOC+MwuGi8NpmsZwoZIx5l2M56cYmwYrnbEr31BbW4vq6uoStsh5GC4Ol8vlCpYhG1KpVMXfQdLMy3/j7O3tnfScrBSjo6PSQpumpiaWY8rDcHG4qcIlk8kwXGjG+Xw+RKNR8/XY2Bji8XgJW1RaQghpYUP+0CFdxHBxuEwmM+nqnGw2y5U7NCvmzZtn/j2bzaK/v79ib2yMcDH4fD5WjC6C4eJw6XRauoitJ/BUvRoiuyiKUrDMtru7u2LDJZVKYWhoyHxdV1dX0Q8FmwzDxeHyJ/OtJ7GmadA0bbabRBUoFApJGwQHBgYq9mmVw8PD0vYA64IH+iX+iziYEEIq8aKqKkKhkPla0zTWF6NZ4XK5pHmFiYkJDA8Pl7BFpWGsljOGo1VVRUtLC4fEimC4OJz17lBRFGnVjq7rDBeaFYqiSM8PMh77UGlDY5qmSftbgsEg6uvrS9gi52K4OJy15+JyuaRhMSFExQ5N0OyLRCJSKZienp6KG5aNx+OIxWLm64aGBumxBPRLDBeHs4aHy+UqeBIgey40W/x+v7QkeWRkRKqtNdcZq8Ss12RbWxuHxCbBcHGw/J6Jx+OB3++XTmb2XGi2KIqC1tZW83Umk6moJclCCKnEvtfrlapGk4zh4mBCCKln4na74fP5pJUp3EhJs8UoBeN2u81j3d3dJWzR7Eomk+jv7zdfh8NhaYENyRguDqbrurSPxePxwOPxSOHCB4bRbKqpqUFtba35uq+vryLOQaMKsnUJcmtrK0u+TIHh4mD54eL1euF2u6UTmsNiNJvcbre0W39iYqJiqiSfO3fO/D2Np05ySGxyDBcHy+Vy0mocj8cDl8slhUs2m62IC5uco7W1VVqSXAlDY+l0Wir5UlNTg3A4XMIWOR/DxcHyd+B7vV64XC5pzJvhQrNJURQ0NDSgqqrKPNbT0zOnyxAJITA4OCitjGttbeVTJy+B4eJguVxOKkzp9XqhKIp0UjNcaLb5fD40Njaar0dGRjA2NlbCFs288+fPS7vy58+fX+IWOR/DxcHyg8Pj8RSEi6ZprIxMs0pRFOnNNZfLzelCltlsVnpAWnV1NZ86OQ0MFwfLD5diPZf83g3RTDOWJFt3pnd1dc3J81AIgZGREYyOjprHmpub4fP5Stiq8sBwcbD83fdGqHi9XvMYKyNTKVRVVSESiZiv8+ck5pLz58+b15iiKFiwYEGJW1QeGC4Olb8739pjsfZc8pcrE80GVVWlN9l0Oo0LFy7MuaExTdOkXfmBQACNjY0cEpsGhouDWXsuqqqaq8SsPReGC5WCUSXZei5++OGHcy5cYrEYRkZGzNdNTU0F9f2oOIaLg1l7LqqqmhP61gs6v0QM0WypqamRhsb6+/vn1NCYUUvMevO2YMEC9lqmieHiYPnhYmyetA6LMVyoVFwulzQ0lkql5tQzXjRNw/nz583Xfr+fDwa7DAwXB7OGhnVnvtGDKfZ5RLNFURS0tbVJNztzaWhsbGwMQ0ND5uvGxkZp8yhNjeHiUMUqIhsFK/N3BrMyMpVKTU0NGhoazNd9fX2Ix+MlbJE9jCEx6zW4cOFC9louA8OlRIQQ5p/JPj5ZuFj/DkzecxFCIJlMYmhoiAFEM8LlcmHhwoXm63Q6ja6urrI/1zRNw7lz58zXfr9feswzXZr70p9CdjPGcgcGBtDQ0IAFCxYUlO4uVm7fOLGNcDHW3herjGw8Ne+1115DMplEOBzGjTfeiKqqKl4gZBtjaMzr9Zrn4ZkzZ7By5cqyLkc/OjqKwcFB83VjYyOqq6tL2KLyw57LLBNC4MyZMzh8+DBOnTpl/jf/Tm+qcClWGTlfLpfD0aNHMTExAV3XMTQ0hPfee6/s7yjJeWpqatDU1GS+HhwclJ4zX26EEDh37hyHxK4Sw8UmlxrmMmQyGZw8edIslSGEQEdHR0HvI3/nvXWeJT9c8oe8hBAYHh4uuMB7enoq4sFONLtUVcXixYvN19lstqwn9nO5HD788EPzdSAQkB4zQNPDcLGB8Wb+7rvvorOzc9JKxUIIDA0NFewFGB8fx/DwsHQsvyCldW+LdUMlULznMjAwUFDrKZlMSjWSirVPCAFd18v2jYFmn6IoaG1tRTAYNI+dPXu2LFcxGtey9XpsaWnhKrErwHCZwnR6I0ZhuxdffBHHjx/HkSNHcPTo0Um/pq+vr+BjRuhYjxcrt2+w7nkBCgtcGt9vsrZOFnzj4+M4fPgwDh06hDNnzszJQoQ0M4LBoPSEypGREQwMDJTlTcqZM2ekWmKLFy9mr+UKMFyKMOp6HT9+HC+//DLOnTs36RutEAIffPABUqmUeezMmTNFx5yNuY9iYrGYdCFms9miz3IBcMmy+7lcbtLloJONhWcyGbzyyis4c+YM+vr68Prrr+P8+fNTrmZLpVLo7+/H6Ogoi2dWOEVRsGTJEnMVo67rOH36dIlbdfnS6bQ0JFZdXc2Nk1eIq8WKEELg+PHjOHnyJACgu7vbfHZ4/kmWTqfR19cnHTOe/xAOh6XPz2Qyk5bHiMfj0HXdvDiLPcvFkF8CJr+Xk0qlpLAr9nOsPR9jAtMafJqm4YMPPkBra6s0BGd8/uDgII4cOYJ4PA6Xy4WGhgZcc801RX8mzX1GGf6amhrzBqarqwsTExNls8pKCIGenh7pxmz+/PnSowVo+thzKSKVSkl3L5qmoaOjo2jvZWxsDMlksuB4X19fwedPTExMOqGeTCalSf3Jyu0Xe61pmrSyLJFITFrMMplMFnxM13Xp9zVM9oTBTCaD119/HWNjYxBCIJfLScs2qTJ5vV5pYn9iYqKs9rwIIXD69GmzvS6XC+3t7SVuVfliuBQxOjpacOc/PDxcECLG5F+xi2dsbEwKEiGE2WswhEIh8++5XE76mfnl9q09FaCw7L4xLCWEwMTEhNQm62RkJpMpCLjx8fGiw2WaphXMBQkh0N3dXfD5kUgEtbW1Bd+DKocxP2E9N0+dOlU2Q6ajo6PSEycjkQgikQiHxK4Qw6WIYoGRTqel0tvWzy0mnU4XDIFZewEulwvNzc3ma03TzPDK352vKIo0NJUfNvl7Yqw/1+VySc871zQNiUTCfG1M/k+2sie/RyKEkIr5uVwu1NfXY/HixVLVAKpMdXV10nk9MDBQcIPiREIIdHZ2Sjd17e3tBUPCNH18NyhCURSEQqGCHcb5F4mmaUWHjYCLJ2v+JL11LNfj8RTcFSUSCfPzrSe5y+UqOMmnCpeJiQnz7263WyqLLoSQfg5w8Q3A+rtbf1b+ZH0qlZJCNhKJ4JZbbuHT+QjAxZWMS5cuNc/rXC6Hjo6OErfq0tLptLQAIRgMcuPkVWK4FLFy5Ups2bIFt9xyi/RgoPweTTqdlobKwuGw9MZsHTrSNE160w8EAgiFQtLdvrVHkV9uPz9cJiu7b9QTM/h8PtTU1Eg/x9oOTdOksAiFQohGo1KbrMNoY2Nj0vBdY2Mj3G53WZf6IPsY5WBqamrMYx9++KGjn/NiFKm03iguWLCAe1uuEsOlCFVV4fV6UV9fL10k4+Pj0vDRxMSE9Lq5uVkKo7GxMXOOJZvNSm/6VVVVCAQCUmhYPz5ZuX1D/hyMEUaapklv/n6/H8FgUAoja7ikUinpdTgclsIlk8mYH8/fJ6MoilQRlwi4eEOzZMkS83UikcDZs2cdOzSmaRref/99s31utxvLli1jr+UqMVymoCgK6urqzNepVMrsXQghpPAAgPr6eulux7pqK5PJSL2R6upqeDwe6U0/mUyamzatn5tfBRkofKaL8fm5XE76Wr/fD5/PJ4WRdVgsHo9LQRYOh6WJeeP3ND7f2svxeDwIhUK8CEmiKAra29vh8/nMYx0dHY7csS+EQF9fn7SdoKmpCQ0NDTyvrxLD5RKs4aJpGuLxuPlGay2l4na7UVNTI/V00um02YtIJBJSEFVXV8Plcklv+ul02iy9Mlm5fYPH45GOGfXFstmsNP8SCATgcrmktfqpVAqaphXMCxlhWl1dXTDvYvz+1t85GAxyDwAVVVtbi7a2NvP18PCwI5clCyFw8uRJaUf+ihUrOMxrA4bLFBRFKZivMO7cjTt6g8/nM+dRDMY8S/7yYEVRzNL31jfnTCZj7ra3hot1d74hP3CM3ko6nZYm4IPBIFRVlYbr0um0GUDWeSGjJxIIBKTQM3ou+XNMNTU1vAipqPw36fw3cScwthJ0dXWZx+rr69HW1sZeiw0YLpdQVVUlde9HR0fNjYPWScqqqiq43W5pmMio1wXIy4ONN/v8cDF6HfmbIvOHwIDilZGBi70S692h8f2tw3XZbNYMIWtABgIB+Hw+uN1u6fMnJiaQy+UK5pisvToiK0VR0NTUJC2D7+3tLVpbr1SMwLMOI69YsaJgPpOuDMPlErxer/RGG4/HzUlz6yoqo4cTDAalN33jrt8aLh6Px3zTt/YoNE0zey/WO7xiJ3v+Ci1jWMw6mW8NL2vFWl3XkUqlkMlkpBVqxvJrVVWlHpjxu+bPMdXW1vIOjyblcrmwatUq8xzRNA0nTpxwREFUIQRGR0dx5swZ81goFGKRShsxXC5BVVVpHiWZTCKZTCIej0u9C+Mu3phAN4yPj0PXdelN3O/3m70Ra89F13VkMpkpKyJb22VdDGDUIrMOW6mqCp/PZw7DGYy9LolEQrprs4aFdVLf6KXlzzGVS80oKg1FUTB//nxpn1VXVxcGBwcd0Xs5efKkdDO2YsUK6WaPrg7DZRrC4bD592w2i7GxMXN4DPjl3IyxAdF6ghr1xKwncSAQMOdLrOFizGvkV0S2hpUhvzJysXCxLhjw+/0Fe13yy9EY4WL9fYx2jY6OSuFizDERTcXj8WDVqlXm61wuh+PHj5c8XEZHR6VNk9XV1Vx+bDOGyyUoioJwOCy9MQ8ODkplX7xer3kXr6qqdEefTqcRj8cLliEbJ7Hf75dOaGO4ynrxGb2P/Hblh0v+HhePx2Ou+vL7/dIKsImJCSksXC6XNBSWv2Is/yFnxlJqoqkoioJFixZJN2jnzp0r6bNehBB47733pBuxZcuWcdOkzRgu0xAKhaQeRm9vr7TfI3/S3/omnc1mMTQ0JA2hWcPH6/VKwWVUR7ZeeMWGxYqV3Tcm6q1fZ8zLeDwe6fMnJiaklWJer1fqieQP7/X39xesFONdHk2Hz+fD6tWrzdfZbBbHjh0rydyLsfz+1KlT5rGqqiqsWLGC57PNGC7T4PV6pTmIkZER6S4+HA6bb+LGkJJB1/WCFTLWnou1dwGgYKFAsYrIBuubv67rSKfTUg/J5/OZweVyuaTwyJ9DCQaD0s/Jn1NJpVLSm4H1TpRoKsaDxOrr681j586dK8nKMSEEjh07VjDXwvlD+zFcpsF4ENJkrMstgYt3QtbeiLUwpMvlkrrfxcLFeuLnT9xbWcPFqKpsXSpsHXIzVrIZjAl9Q/5+nvzqBFaqqnKlGF0Wn8+HNWvWSAUtjx49Oqv7XoQQGBgYkOZaQqEQey0zhOEyDYqioKWlpeibvN/vLygVEQgECp63YrAuQwZ+WcfMkB8uLper6M/NX2lmLHe2Dr/l756f6u6sWJBY7zStfD4fx6fpshjPerHWouvu7p7yUdp20zQNb7/9tnQDtnr1ap7LM4ThMk2hUKho76WlpUXqEQCF8xdW+cFjLBc25O89KVZu32D9OuNhZNY7wfw2WOeCrBRFKeiJGD2XYsFWW1tbdAUb0VS8Xi/Wrl1r9pB1Xcdbb7016dNZ7SSEwNmzZ6Xd+JFIBMuXL2evZYYwXKZJURSsWbNGCpKqqippk5ghf+WVVf5zYoqVgLFWKfZ4PJOWWMkvC2NdHg3ArAJg/Jz8Ev/Wn1GsV1NVVVX092hsbOQFSZdNURQsWLAA8+bNM48NDQ3hxIkTM9p7Mcovvfnmm+YogqqqWLduHW+SZhDDZZqMO/lNmzZh+fLlWLFiBTZt2jTpqqnJJryLDT/ll4Cx9lzyC1Ra5a80s07Q54cWUDhpbz1erKflcrnQ2toqHfN4PJg3bx7Dha6Iy+XC+vXrpfPw2LFj6O/vn7GA0XUdb775pnR9LFiwgA8Dm2EMl8tgBMx1112H9evXTzqprSgK6uvri+5NCYfDBcetb+xCCGloK38fjFV+ryb/AWP5d2WTzZXU19cX7R0Z4+RGIBql1K0r54guh6IoaGxsxIoVK8xj6XQar7zyivnICTsJIXDmzBlp6XEgEMB1113HoqszjOFymYwd7MafydTW1hb0HHw+X8Ebs6IoU+50nypc3G73pCvJ3G53QS9FVdWCh3sZixUmEwwGsWnTJlx//fX4xCc+gWuuuYZ3e3RVFEXB2rVrpbIwAwMDOHz4sLQg5WoZVY+PHDkildS/9tpri978kb0YLjPE7/cXvGk3NjYWDRLr3Eix7zOZ/OfBWHm93qILAebPny8dr6mpmXIOxahLtnTpUvNreVHS1TBuqDZu3Cj1rk+fPo233nrLluXJRhHXl19+WZrDXLhwIZcezxKGywxRFAWrVq1CXV2dWfxy9erVRU9qn883aRc9fyWaVbGhL4Pf7y/4noqiIBKJ4JprrkEoFEJ9fT02bNjASU2adUaPecOGDea8oRAC7777Lo4fP35Vu/eNR2K8+uqr6O3tNY/X1tbihhtumHT1JdmL/8ozxFidtXnzZoyPj6OqqmrSIS6jl1FsSGCqXs1UQ2r5hSoNqqpixYoVaG9vh6IocLlcvIujklAUBStXrsT4+DjeffddCCGg6zreeOMN+P3+KyokacxZvv7669JmSa/Xi1/91V/lY7lnEcNlBhkrti71KGCPxwOfzydtngQuBsGlvnayns1UG8Pyi14SlYqqqrjuuuuQTCbR0dEB4OLu/VdeeQUejweLFi2adhgYj/l+7bXX8N5775mLA1RVxYYNG9Da2spgmUUMFwcwnkxpXSoJFO7mL2ayEGGtJCoHxmMqNm7ciFQqhfPnzwO4uPLxxRdfRDqdxtKlSy/ZwxZCYGhoCEeOHEF3d7f0/a+55hqsWrVq0iX9NDMYLg6gKErRMPD5fFM+ctWYbFcURVrCaXw/3qVROVAUBT6fDzfddBP+53/+B319fQAuBsxLL72Ezs5OLF++HC0tLeYwsXHOG1WO33//fXR0dBQUfV2zZg2uu+46BksJMFwcwlpJ2VBdXX3JtfhVVVVwu91SvSS32816SVRWFEVBMBjEzTffjEOHDpkBo+s6uru70d3djWAwiMbGRjQ3NyMYDCKZTKK7uxu9vb3SHi/g4jWwbt06rF27lhP4JcJ/dQewPpAsv6z9pXoffr8fwWCwoHw+nxJJ5cZYBHPrrbfi1VdfxdmzZ6XrIZFI4OzZszh79uyU36e2thbXX389Fi5cyB5LCTFcHKKmpgbBYNB8ToyqqmhqarpkuLjdbkSjUSlcotEo79aoLBk9mF/7tV/DokWLcPz4cQwODk5r70tVVRWWL1+O1atXT7nKkmYH34EcwufzYdWqVXj77behaRoWLlwo7WCeytKlS9Hb22sueV66dOkMt5Zo5hiT/EuWLMHChQsxPDyMDz/8EF1dXRgdHUUulzPnGD0eD8LhMNrb27Fo0SJzDpJKj+HiEMbT+qLRKLLZLOrq6qZV+8iod7Z582bEYjHU1dXxAqM5wQiZxsZGRKNRrF+/HhMTE4jH48hkMvB4PAiFQubcJM95Z2G4WAghMDg4aGt9oytl3Vl8OWKxGGKxmL2NmSYhRMHEKpUfIQT6+voccR1MxuVyQdd1jI6OFizhLzXjkeOVjuHyEeOBWT09Pejp6Sl1c8pWsYKZVF7C4TDOnz9v7jmhy2dsjK5kipitZ4w6nLFmnv8cV09VVQ5RlCleB/ap9OuA4UJERLbjInAiIrIdw4WIiGzHcCEiItsxXIiIyHYMlzKhaRrGx8dteQQsUbnSNA3xeJzXQRlguJSJ0dFR/PSnP3XchjGi2TQyMoKnn34aIyMjpW4KXQLDhYiIbMdwISIi2zFciIjIdgwXIiKyHcOFiIhsx3AhIiLbMVyIiMh2DBciIrIdw4WIiGzHcCEiItsxXIiIyHYMFyIish3DhYiIbMdwISIi2zFcyoAQAkNDQxgYGMDQ0BCEEKVuEtGsM66DwcFBXgdlQBH8P+RYsVgMu3fvxuOPP47Tp0+bx9vb23H//ffjS1/6Eurq6krXQKJZwOugPDFcHOrAgQO48847kUgkAEC6S1MUBQAQDAaxd+9ebN26tSRtJJppvA7KF8PFgQ4cOIBt27ZBCAFd1yf9PFVVoSgK9u3bxwuL5hxeB+WN4eIwsVgMbW1tSCaTU15QBlVVEQgE0NXVxaEBmjN4HZQ/Tug7zO7du5FIJKZ1QQGArutIJBLYs2fPDLeMaPbwOih/7Lk4iBACy5YtQ2dn52WthFEUBUuWLEFHR4c5Dk1UrngdzA0MFwcZHBxENBq9qq+PRCI2toho9vE6mBs4LOYg4+PjV/X18XjcppYQlQ6vg7mB4eIg1dXVV/X1oVDIppYQlQ6vg7mB4eIgkUgE7e3tlz1erCgK2tvbUV9fP0MtI5o9vA7mBoaLgyiKgvvvv/+Kvnb79u2cxKQ5gdfB3MAJfYfh+n4iXgdzAXsuDlNXV4e9e/dCURSo6tT/e4ydyc8++ywvKJpTeB2UP4aLA23duhX79u1DIBCAoigF3XzjWCAQwP79+7Fly5YStZRo5vA6KG8MF4faunUrurq6sGPHDixZskT62JIlS7Bjxw50d3fzgqI5jddB+eKcSxkQQmB4eBjxeByhUAj19fWctKSKw+ugvDBciIjIdhwWIyIi2zFciIjIdgwXIiKyHcOFiIhsx3AhIiLbMVyIiMh2DBciIrIdw4WIiGzHcCEiItsxXIiIyHYMFyIish3DhYiIbMdwISIi2zFciIjIdv8fh8bVkXch0NMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 500x400 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "104199f4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "function , r2\n",
      "J0 , 0.9985560043309399\n",
      "gaussian , 0.6101756259771707\n",
      "sin , 0.5737221152646913\n",
      "tan , 0.08366297315238909\n",
      "1/x , 0.08315973336762218\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "('J0',\n",
       " (<function torch._C._special.special_bessel_j0>, J0),\n",
       " 0.9985560043309399)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.suggest_symbolic(0,0,0,a_range=(-40,40))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe1f857d",
   "metadata": {},
   "source": [
    "Finish the rest of symbolic regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "eb6c0f43",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "r2 is 0.9985560043309399\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "tensor(0.9986)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.fix_symbolic(0,0,0,'J0',a_range=(-40,40))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "11a27268",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "skipping (0,0,0) since already symbolic\n",
      "fixing (0,1,0) with x^2, r2=0.9999802186534139\n",
      "fixing (1,0,0) with sigmoid, r2=0.9999663092809886\n"
     ]
    }
   ],
   "source": [
    "model.auto_symbolic()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "5076005f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAFICAYAAACcDrP3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAsyElEQVR4nO3deZAcZ3038G/37Ox9r3Z1WrJ3tZKwsQ02BIXDsSGlBUQOY0wCgVAVQi4Kvw45SF6KkLPIUUkkQoKTEAqLJIYKMgm25VcOIcYUBtvYgG2MbXnXRl7J0mqP2Wtmd2ann/ePnx718/Qcu5J6Zrp3vp+qrdX2Xr32dH/7uX6Po5RSICIiCpFb6xMgIqL1h+FCREShY7gQEVHoGC5ERBQ6hgsREYWO4UJERKFjuBARUegYLkREFDqGCxERhY7hQkREoWO4EBFR6BguREQUOoYLERGFjuFCREShY7gQEVHoGmp9AkRxoJTC1NQUFhYW0N7ejr6+PjiOU+vTIoostlyIykilUjh48CCGh4fR39+Pyy67DP39/RgeHsbBgweRSqVqfYpEkeRwJ0qi4o4ePYqbbroJ6XQagLReNN1qaW1txeHDhzEyMlKTcySKKoYLURFHjx7F/v37oZSC53klv851XTiOg3vuuYcBQ2RguBAFpFIpbNu2DZlMpmywaK7roqWlBePj4+ju7q78CRLFAMdciAJuv/12pNPpNQULAHieh3Q6jUOHDlX4zIjigy0XIoNSCsPDwxgbG8P5XBqO42BwcBDHjh3jLDIiMFyILJOTk+jv7y/6uWYAuwA8C2CpzPf39fVV6OyI4oPdYkSGhYWFkp/bBeAbZ9+XMj8/H/YpEcUSw4XI0N7eflHf39HREdKZEMUbw4XI0NfXh6GhofMeN3EcB0NDQ+jt7a3QmRHFC8OFyOA4Dj70oQ9d0PfecsstHMwnOosD+kQBpda5XAUZc3kDgMeNr+c6F6JCbLkQBXR3d+Pw4cNwHAeuW/4S0Sv077zzTgYLkYHhQlTEyMgI7rnnHrS0tMBxnILuLn2spaUFR44cwb59+2p0pkTRxHAhKmFkZATj4+M4cOAABgcHrc8NDg7iwIEDOHHiBIOFqAiOuRCtgVIKsw88gM63vQ1zd9+Nruuu4+A9URncLIxoDRzHkTEV4z0RlcZuMSIiCh3DhYiIQsdwISKi0DFciIgodAwXIiIKHcOFiIhCx3AhIqLQMVyIiCh0DBciIgodw4WIiELHcCEiotAxXIiIKHQMFyIiCh1L7hOtVT4PLCwA7e1AIlHrsyGKNIYL0VopBXge4LosuU+0CoYLERGFjmMuREQUOu5ESZGRz+fx7MMPIzc/X+tTib1tr3gFegcGan0aVMfYLUaRkclk8I2RESwvLSHR3l7r04klpRSaRkcxcNttuOqtb6316VAdY8uFokMp5FwXOz7xCex5/etrfTbxkE4D3/wmcO21QG8vsrkcHnjHO8BnRqo1hgtFTkMigcamplqfRvQpBdx/P/D+9wNbtgA33gi8971wOJONIoAD+kRxlc8Dd9wBLC0BY2PApz4FHDtW67MiAsBwIYqv0VFpuWhXXAG86lU1Ox0iE8OFKI6UAr78ZWBmRj52HOCd75TqAUQRwHAhiqOZGeDOO/2PN24EODuMIoThQhQ3SgFf/Srw3HP+sZERYNu22p0TUQDDhShulpaAf/1XGdAHgNZW4N3vlppnRBHBVyNRnCgFPPoo8PDD/rG9e4FXvILFNClSGC5EcZLPS6slk5GPGxqAX/xFgOuCKGIYLkRxMjYG3Hef//GuXcD117PVQpHDcCGKC88DvvAFYGpKPnYc4Od+DujurulpERXDcCGKi1OngC99yf940yYp+cJWC0UQw4UoDvSiyePH/WM33ghcckntzomoDIYLURxMTwOHDknIANIV9p73sNVCkcVwIYo6pYC77gKefdY/9pa3ALt3M1woshguRFE3Owt89rMyoA9I/bD3vx9IJGp7XkRlMFyIokwp4MgR4Mkn/WNvehNw9dVstVCkMVyIomxuDvjnf7ZLvXzgA0AyWdvzIloFw4UoqpQC7rkH+P73/WNvfCPwYz/GVgtFHsOFKKpSKeC22/xWS0sL8Cu/AjQ21vS0iNaC4UIURXpdS3CsZe9etlooFhguRFE0OQn84z/6rZb2duCDH2SrhWKD4UIUNUpJ5eNnnvGPvfnNwKtfzVYLxQbDhShKlAJeeEFmiOnV+F1d0mppaKjpqRGdD4YLUZR4HvAP/wCcOOEfe+c7ua6FYofhQhQVSgHf+x7wxS/6x7ZuBX7jN7iFMcUOX7FEUZHNAgcOSLkXQFoqH/gAcNllbLVQ7DBciKJAKeBrX7N3mbz8cuC972WwUCwxXIiiIJUC/vIvgaUl+bihAbjlFqCvr6anRXShGC5EteZ5wOc+J+Mt2utfD/zUT7HVQrHFcCGqJaWAp58GPv1pv6R+Zyfwu78rRSqJYorhQlRLy8vAn/85cPq0f+xd7wJe8xq2WijWGC5EtaIU8F//Bdx7r39saEjGWrgRGMUcw4WoFpQCjh+XVks2K8caG4Hf+R1Z28JWC8Ucw4WoFnI5mR02NuYfe8tbgJ/9WQYLrQsMF6JqUwq4+27gS1/yj23eDPz+7wPNzbU7L6IQMVyIqkkXpvzjP7bXtPzWbwF79rDVQusGw4WompaWgD/5E+D55/1j+/YB7343g4XWFYYLUbV4HnDHHcBXvuIf27YN+PjHuaaF1h2GC1E16IrHn/iEDOYDMjvs934P2L2brRZadxguRJWmFDA9LQP2ExP+8be/Hbj5ZgYLrUsMF6JKy+WAv/gL4KGH/GNXXAF87GNAU1PtzouoghguRJWkFPAf/wHcfru9bfGf/RkXS9K6xnAhqhSlgEcfBf7oj/xpx4kEcOutwHXXMVhoXWO4EFWCUsDJk7J+xSxK+ba3Ab/6q6wdRusew4UobEoBCwsygG/u0fLylwN/+qecdkx1geFCFLZsVgpS3nWXf2zDBuCv/gq45BJ2h1FdYLgQhWllBbjtNuCf/snf/KulBfjDPwT27mWwUN1guBCFJZ8HPv95u4x+IgF88IPAz/884PJyo/rBVztRGPJ5Ke3ysY8B6bQccxwJlQ9/GEgma3t+RFXGcCG6WPk88IUvSCmX+Xk55jjA/v2ynoUD+FSHGC5EFyOfB774RTtYAOBNbwL+9m+B7m6Os1BdYrgQXah8Hvj3fwc+8hFgbs4/fsMNwN//PdDfz2ChusVwIboQuRzw2c9Ki8UMluuvBz79aWDjRgYL1bWGWp8AUawoJQP2Bw4An/ykX9YFkBbLpz8NbNrEYKG6x3AhWiulgBMngI9+VBZI5vNy3HGAn/xJ4FOfYouF6CyGC9FaeJ6UzP/t3waefNI/7rrAT/808Nd/DfT1MViIzuKYC9FqcjmZavye99jB0tQE/PqvA3/3dwwWogC2XIhKUUoG6w8ckLGUTMb/XH8/8Ad/IIskk0kGC1EAw4WoGM8DHn8c+PjHgQce8OuEAcBVVwF/8zfAtdeypAtRCQwXIpNSshjyX/5FurumpvzPuS7w1rfKlsXcRZKoLIYLkeZ5/s6R3/ym3VppbZVNvj78YaCjg8FCtAqGC5FSwOSklMr/zGeAVMr+/J490j22b59UOWawEK2K4UL1SylZBHnvvbKR1w9/KMe01lbgF35Btirmwkii88JwofqjlEwvfughmQn29a/Lx5rjyJbEH/uYFKBka4XovDFcqH4oJavqv/tdWU1/333+3itaVxfwS78kG3xt2MBQIbpADBda/3SofP/7Mq5y7712eXxA1qrccINUOH7lK6W1QkQXjOFC65dSst3wd78re9rfd19hqLiurFu59VbgzW8GmpvZWiEKAcOF1hc9ID81Bdx/P/Bv/wZ8+9uF3V+uC+zaBfzarwE33ijdYQwVotAwXGh9UErKszz5JPDlLwNHjgDHj9trVQAJkJ07gV/+ZeDmm4HeXoYKUQUwXCi+lAKWl4Fjx6TL6+67gaeesvdY0VwX2L1bBuvf/nYWmiSqMIYLxYdSfnmWZ58F/vd/gf/+bwmUhYXi39PcLAP0732vlG7hnvZEVcFwoWjTOz+++CLw2GPAN74BPPywfLy8XPx7XBfYtg0YGQHe8Q7g6qs5UE9UZQwXiibPA8bGgM9/HvjqV4EXXgAWF+0V9CbHkfGT175Wur3e8Aa/64uhQlR1DBeKnpUVKcfymc8AExOlv851JUCuvRbYvx+47jppsXBFPVHNMVwoelxX9qovFiytrcD27cBrXiOLHl/9aqn7xUAhihSGC0WP60p5+7vuAqanJTxe9zrp6nrVq4BLLwXa2+VrGShEkcRwoWjas0emDedy8n77do6fEMUIw4UiJeF5ePGRR5BLp2WWVyIhU42feqrWpxYL+ZUVrMzMgBFMtcZwochwXRfqZS8D7rsPJ++7r9anE1tNHR1IdnfX+jSozjlKlZrbSVRdSimsrKyAL8mL19DQANd1a30aVMcYLkREFDo+2hARUegYLkREFDqGCxERhY7hQkREoWO4EK2VUkA+X7p4JhGdw3AhWqvHHwd6euQ9EZXFcCEiotAxXIiIKHQMFyIiCh3DhYiIQsdwISKi0DFciIgodAwXIiIKHcOFiIhCx3AhIqLQMVyIiCh0DBciIgodw4WIiELHcCEiotAxXIiIKHQMF6I1UEphZmYG3tn3inu6EJXFcCEqI5VK4eDBgxgeHsb1N9yAhYUFXH/DDRgeHsbBgweRSqVqfYpEkeQoPoIRFXX06FHcdNNNSKfTAIArlcI3ALwBwBOOAwBobW3F4cOHMTIyUrsTJYogtlyIijh69Cj279+PTCYDpVRBN5g+lslksH//fhw9erRGZ0oUTWy5EAWkUils27YNmUwGnuedO34VcK7lYm507LouWlpaMD4+ju7u7uqeLFFEseVCFHD77bcjnU5bwVKO53lIp9M4dOhQhc+MKD7YciEyKKUwPDyMsbGxgq6wUi0XAHAcB4ODgzh27Bics+MxRPWMLRciw9TUFEZHR61g6QKwE8C7ADQA2AQgGfg+pRRGR0cxPT1dtXMlirKGWp8AUZQszM2hF8AAgP6z7xsBKABNkKex1wHYAWAKwEsATp19vwRgfn4efX19NThzomhhuFB9y+eBqSlgYgKYmEDv2BhGAOQBnAHwDIAJSJD8J4DPQYKkG8BmANsBvPzsj5oF0PPEE0A2C2zeDHR0VPVPIYoShgvVl1wOOHNG3iYmJFg8D2hsBPr70fHa1+LYpZfisRdeQHA4Pw9/rOUMgGNn/90CYAuAa7ZsQefSEnD//fKJtjYJmc2bgU2bgJ6eyv99RBHBcKH1bXn5XKsEZ84AMzOAUkBzMzAwAFxzjbzv6gIcBw6A99x6Kx79zd+Ur1uDDIAxx8EtH/kInJtvlt956hTw0kvyfnTU/52bNvmB09sLuBz2pPWJs8VofVlc9FslExPA3Jwcb2+XEOnvl/dluqxKrXMpZdV1LrmcnIsOm9OnpTsumQQ2bvTDpr8fSCQu8A8nihaGC8Xb3JzdMllclONdXRIiOlBaW8/rx+oV+kqpsgHjui4cx8GRI0ewb9++tf3wfF7O1Wzd5HISLAMDfjfaxo0SQEQxxHCh+FAKSKXsMFlaAhxHxjPMMGlquuhfF6wtZl4qjlFb7M4771x7sBSjlIz96LB56SX/79qwwQ+bTZuka40oBhguFF2eB0xP22GSy8k4xYYNfhfXhg0Ve8JPpVI4dOgQPvnJT2J0dPTc8aGhIdxyyy143/veh66urkr8Yr9V89JLwMKCHO/p8cNm82aZNEAUQQwXio6VFWBy0g+SyUnpQmpo8IOkvx/o66v62IRSCtPT05ifn0dHRwd6e3uruxJ/YcFv1Zw6JeEDAJ2d9iSBzs7qnRNRGQwXqp1s1h58n56WLqKmJj9MBgbkaZ0lVWyZjD1mMzkpx1tb7bDhfzuqEYYLVU8mY4eJfvpubbXDpLOTN8Tzlc1KyOjAOXNGuhWbmvzxms2bpQuR05+pChguVDkLC34X18QEMD8vxzs67GnB7e21Pc/1aGVF/pvrsDl9Wo41NPjTnzdtkv/+DVzuRuFjuFA4lCqcFnx2lhW6u+2ZXC0tNT3VuuR50nVmzkjLZqUV099vz0hrbKz12dI6wHChC6OUrHY3w2R5Wbqz+vrsAXjerKJH//8zJwmk0/7/P92NtmkTHwbogjBcaG10gUfdxXXmjHSzJBKF04LZzRJPc3N22OjqBt3d9iQBdmPSGjBcqLhczp8WbBZ4TCbtwXfWx1q/FhfttTYzM3K8vd1ea8OtnakIhgsJXeBRt0yCBR51oHR3cyZXvVpasmekTU76rxGz+nNfH18jxHCpW+m0PV4yOyvH29rswXcuyqNScjmZhaZbNxMT0n3a2GjPSGNBzrrEcKkX8/N+mExM+AUeOzvtMGE5EbpQuiCnHrc5fbqwIOfmzfJvFuRc9xgu65Eu8GguWAwWeOzvlzcWQqRK0bXhzEkC+nXY32/PSAuh0ChFC8NlPShX4LGvz2+ZVLDAI9Gq9EOPudZGt6B7e+0Zaee5RQJFD8MljlZWrH3frQKPGzb4YVKDAo9E52V+3p6Rpsf+OjvtGWkc+4sdhktcnDxpTwtWSgZOzZlcPT2cFkzxlk7bBTmnpuR4a6vfqtmxg2ODMcBwiYvvf1/KdXR2+m+trZzySetbLieLOWdn5W1hAbj8cmmVU6QxXOJCFx0kqmf5vDxQsYUeeQwXIiIKHeOfiIhCx36WszzPw/zUFLxcrtanEnutPT1oYiXdWPI8D6mJCahsttanEnttGzaguY6nVDNczvI8D4uPPop8Pg+HYxsXzJ2fB/buRdO2bbU+FboAnudh4cEHkc/n4XJN1AVRABKzs8D116N5x45an07N8C5q8BwH7ddcg86BgVqfSnFK+etZIsjzPJz++tfBQbx48xwHnXv3onvLFlmgm8nIGwtSronneTh5772o9+HsaN6lashxHCSitvBQKSmb8dhjssJ5+3aZjhm18wTAW8/64LguEs89Bzz3nFR/SCaBn/kZbhxWTDYrC5n7+vyN8RjCDJfYePpp4Phx+fdTT0m5jC1b+CKmykmlZPEuIDfQVIrhEqSU/Df62tdkn5tLLwV27ar1WUUCZ4vFwcqKrFjWPM+/6IkqZdMm/+HF86TKcZ139RRQChgbk+7q2VngySf9eml1juESB5mMlMUwzczIBU9UKX19dtXsU6cYLkHptP3g19PD6gFnMVziYHFRWi+mTEZKYxBVSkuLvYXx9LSM/ZFQSmr9mQ9+O3aw8vhZDJc4WFwsfGLMZmVrYqJKcV3pGtMyGWkxk2983P93IgFwCv45DJeoU6p4H24+z6dIqizHKRx3eekldo1peptnrbPTbunVOYZLHATHWwC5wDOZ6p8L1ZfeXnuG2KlTHOvTdJVmbeNGdokZGC5RVy5EMhk+RVJlNTfbA9QzM8UfduqNHm8xx0I3b67d+UQQwyXqPE/GV4phtxhVmuPYN83lZX+zunqmlD1LrLFRdoHlurNzGC5Rl88zXKh29LiL3j9FLxqsd3pVvtbVJYso6RyGS9StrBROQ9aWlvgESZUXvHGeOlX6NVkvZmftiTYbN0ayHFMtMVyirly4ZLMMF6q8piagv9//eG4OmJ+v3fnUmlL2xIZg1yEBYLhEXzZbenZOLseZO1QdW7b4/87lZDC7Xh9sdLhoTU2sGF0EwyXqlpfti9h8AZdr1RCFxXEKp9meOFG/4bK0JJMatO5uoI43BSuF4RJ1wcF880Wcz8sbUaV1dNgLBM+cKT3RZL2bnraXB5gTHugc/heJMqXsEi+uKxe5ls+zvhhVRyJhjyssLspNtt7o2XK6O9p15b8Lu8QKMFyiznw6dBx71o7nMVyoOhzH3j9Ib/tQb11j+by9vqW1VaoYUAGGS9SZLZdEwu4WU6p+uyao+vr67FIwJ0/WX7fs/LxsmqZt2GBvS0DnMFyizgyPRKJwJ0C2XKhampvtKckzM3ZtrfVOzxIzr8lt29glVgLDJcqCLZNkUi5w88XMlgtVi+MAW7f6H2ez9TUlWSm7xH5jo101miwMlyhTym6ZNDTInHpzZgoXUlK16FIwDQ3+sRMnanc+1ZbJSJhqPT32BBuyMFyizPPsdSzJpLyZ4cINw6iaOjulHIx2+nR9vAZ1FWRzCvLWrSz5UgbDJcqC4dLYKE+N5gua3WJUTQ0N9mr9xcX6qZJ8/Lj/d+pdJ9klVhLDJcpWVuzZOMmkvKjNcMnl6uPCpujYutWeklwPXWPLy3bJl85O6RajkhguURZcgd/YKMFi9nkzXKiaHEem37a1+cdOnlzfZYiUkvL65sy4rVu56+QqGC5RtrJiF6ZsbJSL23xRM1yo2pqagIEB/+OZGamUvJ69+KK9Kv+SS2p7PjHAcImyYHAkk4Xhks+zMjJVl+PYN9eVlfVdyDKXszdIa2/nrpNrwHCJsmC4FGu5BFs3RJWmpySbK9PHx9fn61ApaZnNzvrHNm2S1huVxXCJsuDqex0qjY3+MVZGplpoa5NyMFpwTGI9efFF/xpzHGD79tqeT0wwXKIquDrfbLGYLZfgdGWianBd+ya7vCwFHddb11g+b6/Kb2mR8SZ2ia2K4RJlZsvFdf1ZYmbLheFCtaCrJJuvxR/9aP2FSyol3WLaxo2F9f2oKIZLlJktF9f1B/TNCzpYIoaoWjo77a6xiYn11TWma4mZD2/bt7PVskYMlygLhotePGl2izFcqFYSCbtrbGlpfe3xks/LeIvW3MyNwc4DwyXKzNAwV+brFkyxryOqFseREijmw8566hqbm5PSNtrAgL14lMpiuERVsYrIumBlcGUwKyNTrXR2ypoP7fRp2VAr7nSXmHkN7tjBVst5YLjUilL+W6nPlwoX899A6ZaLUlLFdWqKAUSVkUjITVdbXpabctxfa/m8FKrUmpvtbZ5pVQ2rfwmFTvflnjkjT33btxeW7i5Wbl+/sHW46Ln3xSoj613zHnlEAqanB3jta6VZzwuEwqK7xhob/dfh888De/bEuxz97Kys3dEGBmRlPq0ZWy7VppRcfA89BDz3nP8++KRXLlyKVUYOWlkBHn9cSqJ7nrRefvjD+D9RUvR0dsoUXW1y0t5nPm6UklYLu8QuCsMlLKt1c2nZLPD0036pDKWAY8cKWx/BlffmOEswXIJdXkoB09OFF/jJk/WxsRNVl+sCl13mf5zLxXtgf2VFzl9rabG3GaA1YbiEQd/Mn3gCGBsrXalYKWlBBNcCLCzI95uCBSnNtS3mgkqgeMvlzJnCWk+ZjF0jqdj5KSXfF9cbA1Wf48jNt7XVP/bCC/GcxaivZfN63LyZs8QuAMOlnLW0RnRhuwceAH7wA+Dhh6U7qtT3nD5d+DkdOubxYuX2NXPNC1AYZvrnlTrXUsG3sCDddPffL11367EQIVVGa6u9Q+XMjDzgxPEh5fnn7Vpil13GVssFYLgUo+t6/eAHwIMPSv9rqRutUsCzz8oCMu3554v3Oeuxj2JSKftCzOWK7+UCrF52f2Wl9HTQUn3h2SzwrW/JuZ8+DXznOzLpoNxstqUlWZU9O8vimfXOcYDBQX8Wo+cBo6O1PacLsbxsd4m1t3Ph5AXibLFilJJgefpp+fjECX/v8OCLbHlZbsYmvf9DT4/99dls6fIY8/NyQeqLs9heLlqwBEywlbO0ZIddsd9jtnz0AKYZfPm8hObWrXYXnP76yUlppc3Py8/asAG48sriv5PWP12Gv7PTf4AZH5cJJXGZZaWUXLfmg9kll9hbC9CaseVSzNKS/fSSz8uge7HWy9ycjGUEnT5d+PWLi6UH1DMZe1C/VLn9Yh/n8/bMsnS6dDHLTKbwc55n/71aqR0Gs1lp2czNyQW5smJP26T61NhoD+wvLsZrzYtS0trS55tIAENDtT2nGGO4FDM7W/jkPz1dGCJ68K/YxTM3ZweJUn6rQevo8P+9smL/zmC5fbOlAhSW3dfdUkrJRW2ekzkYmc0WBtzCQvHusny+cCxIKWnJBb++rw/o6ir8GVQ/9PiE+dp87rn4dJnOzto7Tvb1yRu7xC4Iw6WYYoGxvGyX3ja/tpjl5cIuMLMVkEhIN4KWz/vhFVyd7zh211QwbIJrYszfm0jY+53n89Ky0fTgf6mZPcEWiVJ2Mb9EAujtlZuKy5dT3evutl/XZ84UPqBEkVIy09N8qBsaKuwSpjXj3aAYx5FWRXCFcfAiyeeLdxsB8nXBQXqzLzeZLHwqSqf9rzdf5IlE4Yu8XLgsLvr/bmiwy6IrZf8eQG4AWjDIgoP1S0t2yPb1AW98I3fnI+G6wM6d/ut6ZUW6lKNuedmegNDayoWTF4nhUsyePcC+fXLTNDcGCrZolpftrrKeHvvGbHYd5fP2Tb+lRQLMfNo3WxTBcvvBcClVdl/XE9OammSQ1fw95nnk83ZYdHQA/f32OZndaHNzdvfdwICcW5xLfVB4dDmYzk7/2I9+FO19XnSRSvNBcft2rm25SAyXYlxXWga9vfZFsrBgdx8tLtofb9pkh9HcnD/GksvZN/22NvlaMzTMz5cqt68Fx2B0GOXz9s2/uVmewswwMsNlacn+uKfHDpds1v98cJ2M49gVcYkAeaAZHPQ/TqdlUWVUu8byeeCZZ/zza2gAhofZarlIDJdyHEf6kLWlJb91oZQdHoCEkfm0Y87aymbt1kh7u9zwzZt+JuMv2jS/NlgFGSjc00V//cqK/b3NzXKxm2FkdovNz9tB1tNjD8zrv1N/vdnKSSalpcOLkEyOI+MVTU3+sWPHorliXymZ2WkuJ9i4UR6a+Lq+KAyX1Zjhks/LzVjfaM1SKg0N0soxWzrLy34rIp22g6i9XVoj5k1/edkvvVKq3L6WTNrHdH2xXM4ef2lpkd9jztVfWpK/JTgupMO0vb1w3EX//ebf3NrKNQBUXFeXdI9p09PRnJaslKxnM1fk797Nbt4QMFzKcZzC8Qr95K6f6LWmJn8cRdPjLMHpwY7jl743b87ZrL/a3gwXc3W+Fgwc3VpZXrYH4Ftb5evM7rrlZT+AzHEh3RJpabFDT7dcgmNMnZ28CKm44E06eBOPAr2UYHzcP9bbK6HIVstFY7ispq3Nbt7PzvoLB81ByrY2ueGb3US6Xhdgf62+2QfDRbc6gosig11gQPHKyIC0SsynQ/3zze66XM4PITMgW1rkb21osL9+cVHOJzjGZLbqiEyOI91L5jT4U6eK19arFR14Zjfy7t2F45l0QRguq2lstG+08/P+oLk5i0q3cFpb7Zu+fuo3wyWZ9G/6Zosin/dbL+YTXrEXe3CGlu4WMwfzzfAyK9Z6nnxdNmvPUNPTr13XboHpvzU4xtTVxSc8Ki2RAF72Mv81ks8DTz0VjYKoSsmD4vPP+8c6OlikMkQMl9W4rj2OksnI2/y83brQT/F6AF1bWJCLybyJNzf7rRGz5eJ5csMvVxHZPC9zMoCuRWZ2W7munIvuhtP0Wpd02n5qM8PCHNTXrbTgGFNcakZRbTiO1OYy11mNj8vC3Ci0Xp5+2n4Y273bftiji8JwWYueHv/fuZw8wevuMcAfm9ELEM0XqK4nZr6IW1r88RIzXPS4RrAishlWWrAycrFwMScMNDcXrnUJlqPR4WL+Pfq8ZmftcNFjTETlJJPSetFWVqQobK3DZXbWXjTZ3s7pxyFjuKzGcSRczBvz5KRd9qWx0X+Kd137iX55WW7iwWnI+kXc3Gy/oHV3lXnx6dZH8LyC4RJc45JM+rO+mpvtGWCLi3ZYJBJ2V1hwxlhwkzM9lZqoHMcBLr3UfkA7fry2e70oJVt+mw9iw8NcNBkyhstadHTYLYxTp+z1HsFBf/MmncvJjdnsQjPDp7HRDi5dHdm88Ip1ixUru68H6s3v0+MyyaT99YuL9kyxxka7JRLs3puYKJwpxqc8WoumJuDyy/2PczngySdrM/aip98/95x/rK1NusT4eg4Vw2UtGhvtMYiZGfspvqfHv4nrLiXN8wpnyJgtF7N1ARROFChWEVkzb/6eJ99ntpCamvzgSiTs8AiOobS22r8nOKaytGTfDMwnUaJy9EZivb3+sePHazNzTCkJtuBYC8cPQ8dwWQu9EVIp5nRLQJ6EzNaIWRgykbCb38XCxXzhBwfuTWa46KrK5lRhs8tNz2TT9IC+FlzPE6xOYHJdzhSj89PUBFxxhV3Q8vHHq7vuRSm5Fs2xlo4OtloqhOGyFo4jW50Wu8k3NxeWimhpKdxvRTOnIQN+HTMtGC6JRPHfG5xppqc7m91vwdXz5Z7OigWJ+aRpampi/zSdH73Xi1mL7sSJ8ltphy2fB773PfsB7PLL+VquEIbLWnV0FG+9bN5stwiAwvELUzB49HRhLbj2pFi5fc38Pr0ZmfkkGDwHcyzI5DiFLRHdcikWbF1dxWewEZXT2AhcdZXfQvY84LvfLb07a5iUkuKZ5mr8vj5g1y62WiqE4bJWjiPNejNI2trsRWJacOaVKbhPTLESMGaV4mSydImVYFkYc3o04FcB0L8nWOLf/B3FWjVtbcX/joEBXpB0/hxHStlv2eIfm5qShZWVbL3o8kuPPeb3IrgucPXVfEiqIIbLWukn+euuk6ed3bvl36VmTZUa8C7W/RQsAWO2XIIFKk3BmWbmAH0wtIDCQXvzeLGWViIBbN1qH0sm5ebAcKELkUgAr3yl/Tp88kmZjVipgPE8CRbz+ti+nZuBVRjD5XzogLnmGrlASg1qO46MVxRbm9LTU3jcvLErZXdtBdfBmIKtmuAGY8GnslJjJb29xVtHup9cB6IupW7OnCM6H44jLd/du/1jy8vAt77lbzkRJqWkxIs59bilRa5hFl2tKIbL+dIr2PVbKV1dhS2HpqbCG7PjlF/pXi5cGhpKzyRraChspbhu4eZeerJCKa2t0kJ79auB170OuPJKPu3RxXEcGXsxy8KcOQM89JA9IeVi6arHDz9sl9R/xSuKP/xRqBguldLcXHjTHhgoHiTm2Eixn1NKcD8YU2Nj8YkAl1xiH+/sLD+GouuS7dzpfy8vSroY+oFq7167dT06KgP8YUxP1kVcH3zQHsPcsYNTj6uE4VIpjiOD/d3dfvHLyy8v/qJuairdRA/ORDMV6/rSmpsLf6bjyNPilVfKQH1vL3DttRzUpOrTLeZrr/XHDZUCnnhCao9dzOp9vSXGt78t1TS0ri7gNa8pPfuSQsX/ypWiZ2fdcIOsP2lrK93FpVsZxboEyrVqynWpBQtVaq4rT25DQ/L9iQSf4qg2HAfYs0eujyeekFDwPODRR+X1eyGFJPWY5Xe+Yy+WbGwEfvzHuS13FTFcKknP2FptK+BkUloP5uJJQIJgte8t1bIptzAsWPSSqFZcVwbXMxng2DE5trIiA/zJpBS9XGsY6G2+H3lEClPqyQGuKy2krVsZLFXEcIkCvTOlOVUSKFzNX0ypEGGtJIoDvU3F3r3ycPXii3I8mwUeeEBmku3cuXoLWylZM/Pww7Ly3/z5V14pXdSlpvRTRTBcosBxiodBU1P5LVf1YLvj2FM49c/jUxrFgePIa/0NbwD+53+koCUgAfPNbwJjY7K2bPNmv5tYv+Z1leNnnpGWT7Do6xVXSMuIwVJ1DJeoMCspa+3tq8/Fb2uTJz+zXlJDA+slUbw4jnTxXn89cP/9fsB4nrRETpyQzw8MSBmm1lbpSjtxQgbtzTVegFwDV18tU545gF8T/K8eBeaGZMGy9qu1Ppqb5UILls/nLpEUN3oSzJveJDO9XnjBvh7SaTn2wgvlf05Xl6zL2rGDLZYaYrhERWenhILeJ8Z1gY0bVw+Xhgagv98Ol/5+Pq1RPOkWzE/8hAzm/+AHsvPrWta+tLVJ99nll5efZUlVwTtQVDQ1yaDj974nF9KOHfYK5nJ27pSuAT3leefOip4qUUXpQf7BQbkOpqeBH/1IKhrPzspsMj3GmExKC39oSMJIj0FSzTFcokLv1tffL+Mn3d1rq32k653dcIMMbHZ38wKj9UGHzMCAXBevfKWstp+flzGWZFK60fTYJF/zkcJwMThKYXFyEl6Y9Y0ulLmy+HykUvJWA0opeNkseInHm6MUFk6fjsZ1UEoiIeMxs7OFU/hrTHkevOXlur8OGC5nOXrDrJMnkTl5stanE1uJhga45aZPU6Q5ANDTA/Xii1jUa07ovCWSSTh1XlbJUapae4xGm1Lq3BtdHNd1JawpdngdhKferwOGCxERhY6TwImIKHQMFyIiCh3DhYiIQsdwISKi0DFc4iKflxX4YWwBSxRX+bwsouR1EHkMl7iYnQW+8pXILRgjqqqZGeCOO+Q9RRrDhYiIQsdwISKi0DFciIgodAwXIiIKHcOFiIhCx3AhIqLQMVyIiCh0DBciIgodw4WIiELHcCEiotAxXIiIKHQMFyIiCh3DhYiIQsdwISKi0DFcYkAphampKZw5cwZTU1NQStX6lIiqTl8Hk5OTvA5igOESYalUCgcPHsTw8DCGd+3C/7n1Vgzv2oXh4WEcPHgQqVSq1qdIVHHmdbB7zx78349+FLv37OF1EHGOYvxH0tGjR3HTTTchnU4DALqVwpsB/D8AKccBALS2tuLw4cMYGRmp3YkSVVDwOuhVCm8HcCeAaV4HkcaWSwQdPXoU+/fvRyaTgVKqoPmvj2UyGezfvx9Hjx6t0ZkSVQ6vg3hjyyViUqkUtm3bhkwmA8/zzh3vAc61XMwNXl3XRUtLC8bHx9Hd3V3dkyWqkFLXQR9wruUyZXw9r4PoYcslYm6//Xak02nrgirH8zyk02kcOnSowmdGVD28DuKPLZcIUUpheHgYY2NjBV0ApVouAOA4DgYHB3Hs2DE4Z/uhieKq3HVQquUC8DqIGrZcImRqagqjo6PnPcVSKYXR0VFMT09X6MyIqofXwfrAcImQhYWFi/r++fn5kM6EqHZ4HawPDJcIaW9vL/m5OUiX2FyZ7+/o6Aj7lIiqrtx1kIJ0iaXKfD+vg2hguERIX18fhoaGivYX5yFjLfki3+c4DoaGhtDb21vpUySquNWugynwOogDhkuEOI6DD33oQxf0vbfccgsHMWld4HWwPnC2WMSUmt9fCuf303rE6yD+2HKJmO7ubhw+fBiO48B1y//vcV0XjuPgzjvv5AVF6wqvg/hjuETQyMgI7rnnHrS0tMBxnIJmvj7W0tKCI0eOYN++fTU6U6LK4XUQbwyXiBoZGcH4+DgOHDiAwcFB63ODg4M4cOAATpw4wQuK1jVeB/HFMZcYUEphenoa8/Pz6OjoQG9vLwctqe7wOogXhgsREYWO3WJERBQ6hgsREYWO4UJERKFjuBARUegYLkREFDqGCxERhY7hQkREoWO4EBFR6BguREQUOoYLERGFjuFCREShY7gQEVHoGC5ERBQ6hgsREYXu/wOIupGGO/6M3wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 500x400 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "79816b25",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "train loss: 1.12e-03 | test loss: 1.17e-03 | reg: 4.76e+01 : 100%|██| 20/20 [00:08<00:00,  2.38it/s]\n"
     ]
    }
   ],
   "source": [
    "model.train(dataset, opt=\"LBFGS\", steps=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "ba171cc4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAFICAYAAACcDrP3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAu5ElEQVR4nO3de3BcZ3038O85q9vqYt0sXxVblqLgODj3BAjBBRJsEgOBuqEpmZAwMB2gE0OmfVsuvVCYt0P7FmKHkqYwA7H5ow3EJqTYVEnTQCEJhNxI4uDYkS3bsuKb5NV1Je3ued4/fn58nme1K8v2WZ2z2u9nxhP2SCsdmz3ne57b73GUUgpEREQBcsM+ASIimnsYLkREFDiGCxERBY7hQkREgWO4EBFR4BguREQUOIYLEREFjuFCRESBY7gQEVHgGC5ERBQ4hgsREQWO4UJERIFjuBARUeAYLkREFDiGCxERBa4s7BMgKgZKKfT392NkZAS1tbVobm6G4zhhnxZRZLHlQjSNRCKBzZs3o7OzEy0tLVixYgVaWlrQ2dmJzZs3I5FIhH2KRJHkcCdKoty6urqwYcMGjI2NAZDWi6ZbLdXV1di2bRvWrVsXyjkSRRXDhSiHrq4urF+/HkopeJ6X9/tc14XjONixYwcDhsjAcCHKkkgk0NraimQyOW2waK7rIh6Po7e3Fw0NDYU/QaIiwDEXoixbtmzB2NjYjIIFADzPw9jYGLZu3VrgMyMqHmy5EBmUUujs7MS+fftwNpeG4zhob2/H3r17OYuMCAwXIsuJEyfQ0tKS82txACsB7AaQnOb9zc3NBTo7ouLBbjEiw8jISN6vrQTwwqn/5jM8PBz0KREVJYYLkaG2tva83l9XVxfQmRAVN4YLkaG5uRkdHR1nPW7iOA46OjrQ1NRUoDMjKi4MFyKD4zi4++67z+m9Gzdu5GA+0Skc0CfKkm+dyxWQMZcrAbxofD/XuRBNxZYLUZaGhgZs27YNjuPAdae/RPQK/e3btzNYiAwMF6Ic1q1bhx07diAej8NxnCndXfpYPB7Hzp07sXbt2pDOlCiaGC5Eeaxbtw69vb3YtGkT2tvbra+1t7dj06ZNOHz4MIOFKAeOuRDNgFIKg08+iYYbbkDiiSdQ/573cPCeaBrcLIxoBhzHOT2m0tDQADBYiKbFbjEiIgocw4WIiALHcCEiosAxXIiIKHAMFyIiChzDhYiIAsdwISKiwDFciIgocAwXIiIKHMOFiIgCx3AhIqLAMVyIiChwDBciIgocS+4TEVHg2HIhIqLAMVyIiChwDBciIgocd6KkyMhkMtj9zDNIDQ2FfSpF74KrrkLzwoVhnwaVMA7oU2Qkk0k8ef31mEgmEaurC/t0ipJSCpV792LRli24/EMfCvt0qISx5ULRoRRSrosV996LVe9+d9hnUxxGR4Ff/hK45hqguRmTqRT+56abwGdGChvDhSKnrKwMFZWVYZ9G9CkFPP44cNttQGsr8NGPAp/4BBzHCfvMiDigT1S0MhngBz8AxseBN94AvvENYO/esM+KCADDhah4vfEG8MQT/utLLwXe9rbwzofIwHAhKkZKAQ8/DPT3y2vHAT72MaC2NtzzIjqF4UJUjAYHgYce8l8vXAjcckt450OUheFCVGyUAn7+c2D3bv/Y+98PLFsW2ikRZWO4EBWbVAr43veAdFpeV1UBd9wBuLycKTr4aSQqJkoBL78sLRft6qtlIJ9TkClCGC5ExUQpYMsWYHhYXrsucNddQHV1qKdFlI3hQlRMenuBH//Yf71iBfCBD7DVQpHDcCEqFnr6cV+ff+yP/xhYsCC8cyLKg+FCVCwSCekS03XDmpqA229nq4UiieFCVAyUAn72M2DXLv/YTTcBF10U3jkRTYPhQlQMxsaABx6QemIAEI8Df/qnQCwW7nkR5cFwIYo6vWjyN7/xj61ZA1x7LbvEKLIYLkRRNzkJ3H+//BcAKiqAz3wG4LYEFGEMF6IoUwp4+mngySf9Y9dcA9xwA1stFGkMF6IoS6WAb30LSCbldVkZ8OlPAzU14Z4X0RkwXIiiSing178GHnvMP3b55cAHP8hWC0Uew4UoqlIp4L77gNFReR2LAZ/9LDBvXrjnRTQDDBeiKFIKeOopWduirV4NfPjDbLVQUWC4EEXRxARw772yvgWQVsvddwMNDaGeFtFMMVyIokYp4IkngMcf949deSXwh3/IVgsVDYYLUdQMDwNf/zowPi6vy8qAz30OqK8P97yIzgLDhShKlAJ++EOZJaZdfz1wyy1stVBRYbgQRcnRo8A3vuFvYVxdDXzhC1zXQkWH4UIUFZ4HfOc7wOuv+8c+9CHgPe9hq4WKDsOFKAqUAvbulcrH5n4tf/mXQHl5uOdGdA4YLkRRkMkA3/wm8Oab/rG77gIuvZStFipKDBeisCkF/OpXwL//u39s+XJZ1+LyEqXixE8uUdiGh4G//3v5LyCBsnGjBAxbLVSkGC5EYVIKePBB4Je/9I9dey3wiU8wWKioMVyIwqIUsHs38E//5G9fXF0N/N3fscwLFT2GC1FYJiaAv/1b4PBh/9jHPsaNwGhOYLgQhUEp4KGHgEcf9Y91dABf+pKUeyEqcgwXotmm17R85SvA5KQcq6iQVkxbG1stNCcwXIhmWzIJfPGLQE+Pf+zDHwY++lEGC80ZDBei2eR5wPe/b3eHtbUBX/saUFkZ2mkRBY3hQjRblAKefx746lf9wpSVlRIsnZ1stdCcwnAhmg1KASdOAJ//PHDsmH/89tuBW29lsNCcw3Ahmg2Tk7J+5Zln/GOrV8vK/IqK8M6LqEAYLkSF5nnA1q3A977nVzxuaADuvRdYupStFpqTGC5EhaSLUn75y7JoEpB1LF/+MvDudzNYaM5iuBAVilLAvn3AZz8LHD/uH7/1VjkWi4V3bkQFxnAhKgSlgIEBCZFdu/zjV10ltcTi8fDOjWgWMFyIgqaULJT8q78CHn/cP75kCXD//RxnoZLAcCEKWioF/MM/AFu2+AP4tbXA5s3ANdcwWKgkMFyIgpROA//6r8A//7O/ULKiQqYhf+QjDBYqGQwXoqBkMsAPfmDPDHNdGXe5+24O4FNJYbgQBcHzgB/+ELjnHmB0VI45DnDbbVLuhQslqcQwXIjOl+cBP/qRtFAGB/3jN98M3HefjLewO4xKDMOF6HxkMsB//Afw6U8DiYR//L3vBb77XaCpicFCJYnhQnSu0mkpn/+Zz9jBsmYN8OCDwKJFDBYqWQwXonMxOQls2gR87nPA0JB/fM0aqSPW2spgoZLGzbqJzoZSMmD/1a/KuhW9TTEA3HCDtGQYLEQMF6IZU0r2Yvn852UAP5OR444DfPCDwL/9G7BwIYOFCAwXoplRCnj1VRlfeeop/3gsJht+3Xsv0NjIYCE6hWMuRGeSTgOPPCKtEzNYKiuBv/gLqRfGYCGysOVClI9SwPCwlHL55jf9xZGAbPb1f/8v8KlPcYEkUQ4MF6JcdDfYn/858MQTslBS6+gAvv1t4H3vk/IuRDQFw4XIpMvl/+AHMiOsr8//muvK4sh/+RfgoovYDUY0DYYLkeZ5wGuvAX/zN8BPf+pXNQaA6mrgz/4M+NKXgPp6BgvRGTBciJSSmmDf+Y6MrRw9an+9s1N2j/zAB2R2GIOF6IwYLlS6lJJFkI8/Dnzta8Bzz9ljKxUVwEc/Kl9bvpyhQnQWGC5UepSSBZAvvgj84z8CO3YA4+P291x4IfCVrwB/9EcSMgwWorPCcKHSoZQ/rrJ5s6yyN+uCAVIe/447gC98AbjgAoYK0TliuNDcp1sqr74qWxD/6EfAyZP298RiwPXXy3bE73oXx1aIzhPDheYupWS74WeflbpfO3bYm3kBEiAXXywr7W+9FaipYagQBYDhQnOLUvLfEyeAri7ge98Dnnlm6piK4wBtbVIr7K67gPnzGSpEAWK40NygWym7dgEPPQRs3w7s32/P/gL8UPnUp4A77wSWLGGoEBUAw4WKl1Ky0LGnB/jZz4Bt24Dnn7drgGmxGLByJfDJTwK33cZdIokKjOFCxUWvTenpAf7nf4Cf/ETGVLIH6LV4HHj72yVUbr5ZCk4yVIgKjuFC0aaUv4L+tdckUB57DHjllamD85rjAEuXyor6O+4ArrxSyuMzVIhmDcOFoml0FNi3D/j1r4Enn5TWyaFD9rbC2errpZVy223A2rV+1xdDhWjWMVwoepQCvvhFmemVa/zENG8ecPnlwC23SLdXeztQXs5AIQoZw4WiR8/oyhUsrivThq+6SsLkhhskUFiihShSGC4UTdddB1RVyfqU2loJkHe+E7jxRuDaa6XLi6voiSKL4ULRtHKl7J9y8cUyjtLWJnuqAAwUoiLAcKFIiXkeDjzzDFKjo8CaNXKwu1v+0Bml02mkBwbA+KWwMVwoMlzXhbrkEmDnTvTu3Bn26RStynnzUN7YGPZpUIlzlNLFmIjCpZRCOp0GP5Lnr6ysDK7rhn0aVMIYLkREFDg+2hARUeAYLkREFDiGCxERBY7hQkREgWO4EBFR4BguRDP1wgtSHeCFF8I+E6LIY7gQEVHgGC5ERBQ4hgsREQWO4UJERIFjuBARUeAYLkREFDiGCxERBY7hQkREgWO4EBFR4BguREQUOIYLEREFjuFCRESBY7gQEVHgGC5ERBQ4hgvRDCilcPLkSQDAyZMnoZQK+YyIoo3hQjSNRCKBzZs3o7OzEzfceCMA4IYbb0RnZyc2b96MRCIR7gkSRZSj+AhGlFNXVxc2bNiAsbExAMDlSuEFAFcCeMlxAADV1dXYtm0b1q1bF96JEkUQWy5EOXR1dWH9+vVIJpNQSk3pBtPHkskk1q9fj66urpDOlCia2HIhypJIJNDa2opkMgnP804fvwI43XJ50fh+13URj8fR29uLhoaG2T1Zoohiy4Uoy5YtWzA2NmYFy3Q8z8PY2Bi2bt1a4DMjKh5suRAZlFLo7OzEvn37pnSF5Wu5AIDjOGhvb8fevXvhnBqPISplbLkQGfr7+9Hd3W0FywIA1wL44qnX7QCqst6nlEJ3dzcGBgZm50SJIq4s7BMgipKRwUEsBbACwHIAbQDiADwAlae+5w8BXA+gD0C38WcEwPDwMJqbm2f5rImih+FCpS2VAnp7gZ4eYP9+tLz+Oj4DIAXgIICnAfQAOASgHMBKAMcALAXQAWAVgDWnftQxAM2PPQZccQXQ0QE0NgLsIqMSxTEXKi3j48DBgxImPT0SLJkMUFUFtLVBLV+ONR//OJ45cACZGf7IeZCguW7hQnzr85+H8+ab8oWGBgkZ/WfhQoYNlQy2XGhuGx31g6SnB3jzTUApoK4OaGsDbr5Z/nvqxu8A+KN77sFT99wj3zcDQ5BFlXd96UtwNm4ExsaAffuA7m758+KLgOcBNTV22CxdCrgc9qS5iS0XmlsSCeDAAWD/fgmT48fleGMjsGKFBElbG9DUlLcVkW+dSz5nXOcyMSHnosPmwAHpjquslHO68EIJm2XLgDI+79HcwHCh4qUU0N/vB0lPj4QLACxYICGyYgWwfDlQX39WP1qv0FdKTRswruvCcRzs3LkTa9eundkPT6ela85s3UxMSLAsXy5h094u515ZeeafRxRBDBcqHp4HHDlit0xGR6UFsmSJ3yppawOqq8/712XXFjMvFceoLbZ9+/aZB0sungf09flB090NjIxIl1lrq9+N1t4uXWtERYDhQtGVyQCHD0uQHDggYaKf8Ftb/SBZtqxgT/iJRAJbt27Ffffdh+7u7tPHOzo6sHHjRtx5552oP8tW0RkpBRw7ZofNqXL/WLTIb9lceOFZt8iIZgvDhaJjchI4dOj0tGAcOiRdSBUV0l2kw6S1ddbHJpRSGBgYwPDwMOrq6tDU1DS7K/EHBqQb7Y03JGyOHZPjzc122DQ3c0YaRQLDhcKTTPotkp4eaaV4nnRpmV1cixdzVlW24WF7zObwYWnxzJtnz0hbvJhhQ6FguNDsGR62pwUfPerfEPXge1sb0NLCG+LZSialtafD5uBB6VasrpZWjW7ZtLYCsVjYZ0slgOFChaGUjBOYYdLfL19rbvZnca1YIYsNGSbBmpyUgNHdaD09cqyiQgK8o0PCZvlyoLw87LOlOYjhQsFQStaUmNOCh4YkNBYutKcF19WFfLIlKJORagTmJIFkUloxy5b53WgrVgDxeNhnS3MAw4XOjefJanc9+H7ggKxMd11Zea7HS5Yv580qipSS///MsNEPA0uX2tOf+TBA54DhQjOTTvsFHnt6JEwmJ2XW1rJlfphccIF0vVBx0QtSdTdad7ffjblggd+N1t4u1Q2IzoDhQrlNTPgFHvfvtws8mtOCly7lAPFclUhIyOgp0EeOyPHGRntG2oIFHDOjKRguJMbG7MH3vj55mq2psWtyLVzIacGlanTUnv7c2yvdo7W1dtgsWcLPCDFcStbgoB0melFeQ4MdJlyUR/lMTNjTnw8ckO7Tqip/+jMLcpYshksp0P3pZpjociItLfZMrlxVfYlmIpXyC3K+8YYEz8SETHVevtxv2bS1sSBnCWC4zEWeJwsUzTAZGZEWyOLF9up3FkKkQvE8qRzQ3S1hs2+fdK25rkz8MGekBVBolKKF4TIX6AKP5kyu8XEZaDcLPC5fzidGCo9S8tBjTn/WWyQsWeJXEWhvZ0HOOYDhUoxSKSnqqBcsHjokxyoq7GnBra1cfU3Rpas46JZNd7e/udv8+fYkAY79FR2GS7HYs0fCZP9+v8BjPO63SFaskC4vTgumYjY05I/Z7Nvnz1qsr/dbNqtXs2VTBBguxeK552RwtKFB/tTXy3gJn+ZoLkulJHAGB+XP8DBwySXSkqFIY7gUi1RKpnMyTKiUZTJyDXAdTeQxXIiIKHCMfyIiChyXzZ7ieR6OHzoEb3w87FMpevVLl6K6tjbs06Bz4Hkeju7fj8zERNinUvQaW1tRM29e2KcRGobLKZlMBokdO5BJpeBwLci5UQqx/n5gwwZUX3xx2GdD5yCdTuPk9u1I8zo4L2XHjwO3346a1avDPpXQMFw0peA5DhpvvhktK1aEfTa5eZ7Ubiovj+TAfiaTwb6tW8FBvOKWcRw0feQjWNTZKQPoIyMyS4sFKXNLp+XfqaICcBxkMhnseeCBsM8qdAyXLK7roixqRfaUkumYO3fKCudLLwXWrIlkMcDoRR6di5jrouy552QK/JtvSmWHe+7hxmHZlAJ27ZJr84ILZB3OxRfzOgAH9IuDUsDPfy4f4hMngF/8QhaZcaIfFdKxY8DevdJySSTkwYam2rtX9rr57W+Bhx8GBgbCPqNIYLgUg/FxCRMtkwF+/3uGCxVWe7vfDZbJSHUIfuZsqZSUYNIaGqR0DTFcikIiIX3epiNH5IInKpQlS+xqxfv2ybgf+U6elN4E7YILZD8bYrgUhf5+GTQ0DQ4CyWQ450OlYd48YNEi/3Vfn5TMJ6GUFI01r8OOjkhOtgkDwyXqdOXYbBMTvNCpsGIxKYiqDQ9Li5l83d3+/9abohEAhktx0HtemNLpqV1lREFyHHkSN8dd9u3juIs2OWmPtzQ1cbzFwHCJOj0NeabHiYK0eLE9/bi7m2N92sCAdFlry5ZxMz4DwyXq9CK2XIaG+BRJhVVXJwGjvfkmW8yAXHc9PdI9rV14IcdbDAyXqMtk8g/c8yKnQnNd6RrTRkdls7pSp5S9PKCyUsZbGC6nMVyiLpWyn45Mo6NsuVBhOY6sd9E7nHoeF/ACsvbswAH/dUuLjLnQaQyXqBsfl4DJZWyM6w6o8BYtksWB2r59+T+TpeLYMXsW54oVUluMTmO4RN34+NQ1LloyyXChwquuBlpb/dfZN9ZSo5QErL4uHUfGW8jCcIm66QJkYiJ/8BAFxXGAzk7/dTIJHDxYul1jnif1xLTqalmZz/EWC8Ml6pJJ+yI2y3GkUuyeoMJzHOn2MafZ7tlTuuEyPGxPali8GKivD+98IorhEmVKybiK5rpAc7P/errBfqIgtbTYCwR7eqTLttQoBfT22jM1Ozv9CQ90GsMl6sxwicWAxkb/dSYjq4SJCq2iQmaNaQMDpVuCf+9ev6s6FuP6ljwYLlFnrnGJxexZO57H4pU0ey66yC8Fk06X5pTkyUm7nlhDg13ck05juESZUnZ4lJdLpVrz66XYNUGzz3FkxlhNjX9sz57SKwXT3y+z5bS2NnsclE5juERddrjU1dlNcIYLzZb6emDpUv91X59s/VAqlJJWiznOedFF7BLLg+ESZZ5nh0dFhTwlucb/bewWo9niusBb3uK/Hh2VVeql0jXmecDrr/uv43GZRcdwyYnhEmWeZw/YV1bKLnfmzJTsqcpEhaIXC5aXy2ulSmu77ZERWd+jLV7Mki/TYLhEWfZssMpK+VNW5h9jy4Vm04IFMi1Z27+/ND6DSkmwZE9BNq9FsjBcoiyTsRdJVlVJ15j5gR4fL50nRwpfZaVd6mRgQMrwl8JncPdufwpyWZl0EbJLLC+GS5SlUnZ5l6oq+VDrbgmA4UKzb+VKe3fK3bvDPZ/ZMD5ul9hvarL3uaEpGC5RNjk5NVxiMbv66sQEw4Vmj+NIHS1zSvzrr8/9MkRHjgAnTvivOzrkeqS8GC5RNjlpF62Mx+XiNms8ZX8PUaHV1sosKe3oUeD48fDOp9CUktaZftBzXWDVKnaJnQHDJcomJuzgqKqSD3Z2uJTaQjYKl+MAF1/sv56YmNuFLFMpu+uvrg5YtozhcgYMlyjLLkqpm+FmuKTTLLtPs8txpFvIXK3/2mtz9yHn+HGZtKC1tUnA0LQYLlFmjqfo7jDHsft6GS4UhoYGGXvRenvn5gZiSkmrzHzQW7XKXshMOfFfKKqy64aZYy1muGRPVyaaDbGY3GS1sTGpFjzXusYyGWmVaTU1rII8QwyXKDPDxXX9WWIMFwqb48g6D/Oz+Oqrc29ySX8/cOiQ/3rZMrsyOeXFcIkysykei8n6luxuMaW4YRiFo7lZKiVrBw7Mra4x3SVmViC45BJuDDZDDJcoMz/U5uJJc0Cf4UJhKSuTm602MjK3usYyGWmNafE4qyCfBYZLVGWHRlmZX/ZFD+zn+j6i2TLXu8b6+6U1pi1bZm8zTtNiuERVrnDRzXEzXACWgKHwLFhg7/Gyf//c6BpTSioPmL0Hb30ru8TOAsMlDEpJddW+vvzBkB0uFRX+9EfzfwP5NwxTChgakhXUk5MMIApeWZncdLXRUVlwWOyftUwGePll/3U8zkKVZ4n1omebUjL75OGH5ca/dClw662yy5/5wc3ey2W6cMnVLaYHIx99VKaJdnQAH/mIvfCN6Hzp1fqPPeY/5b/8MvD2txd3Ofpjx+xZYm1t7BI7S2y5zLZ0Gvjv/5ZS5em09Ok+88zUJ71c5fZ1oJSXTy27n21iQn7P4KBfvuLFF4v/iZKip6VFxiO0gwftIo/FRilg1y77urr0UnaJnSWGy2w7cUJWM5v27LFbKcDUcDErIZuD+8DUrjWlpMvt2DH7Z+7axdX8FLxYTG6+WjIpn7VifZCZnLS7xOrq2CV2DhguQVFKwiCTyX9RKSXBkh0kiYS0ZEz6Z2lVVf6HOxazwyVXt1hv79RaTydOSFfcdH8HXU6mWG8MNPscR/Z4qa31j7388tTPeTHQD2ZmLbGODi6cPAdF3CkaIamUdG397neyz8XatcCiRbmfdA4fzv3+Y8dk8yH9nuyaYeZ0z1x7unie322mL5Bs4+NShC9X37HnyTTSZ5+V97/73SxzQTPX2Cg34d/9Tl4fPix/Vqwors+QUsBLL/m9Bq4LXH45a4mdA/6L5aOf4lOp6Z/ilZKb8hNPyKysvXuBn/wk9zhIJjO1q0rLPp69T4sZLq5r70aZPRMsnc7d562UbHqU/ffRLapHHwV6emQc6JFHpramst+jJx14Hls6pc51gSuu8G/CqZTcpIvtczE2Jl16WlMTH7LOEcMlF6Vkvv6WLcB3vwv89rf5y4mnUlO/3tcn78++sCYmZIA9l/5++/tzhYv+gOfaMMz8/ePjslo6l+PHc4fLc8/ZgTg4KF0b+aZJj40BP/0p8MADwI9/LFOrqXQ5DtDZKTdj7dVXZWpysVBKHg77+/1jq1ZxhuU5YrjkMjAAbN8uAdHXB+zcKQuqct1oBwakNWDyvNzfPzIiN+VcEgk7ILK3LzbDJDtc0mn7vSMj+de+DAxMXUE9Pi4tlmx79uSeAOB5MhPt2WelxfXii8COHSygWeqqq+01LydP5r9uosjzgOef96+P8nJpjbHVck4YLtk8D3jqKbnZa+m0jKlk32jzDdADcjz7ZptI2D/D/NCOjNg/x5wBll2s8kx7ugwO2mFj/p6hoakTAE6cyN2iyjcBoL/frrkEAL//vV2anEqP48jNWHfZ6pt1sWwidvw48MYb/uvWVvnDcDknDJdsnieD5WbLAJDBSbO5DMjN/+DB3D9ncNC+MSslrQYzMFpa/K8nk3arxgwAx7EH8IGp4aKDTCn53TP9PUrJzJhcN4DxcRlHyp7mvHevXRYjFpMn1uXLp/4MKh2OAyxZYn8O9u3LPc4XNXog3/xcX3mlPbZJZ4Xhkq2sTGZ73XWXDORpk5My0J09cH70aO6fMzExNYzMmkvl5fbCs3TaHyfJ3ijMdad2i2Xv6WK2esxWV67fk90aMaddmnLNOvM8uWFoFRXABz8oq/85XZPKy4GrrvKf9sfHgRdeCPecZmJsTLp3tXnz5IGJrZZzxnDJxXWlOXzzzfZNPDtckkn7Rl5XZ08HNp/6lbK/t6pKph5rnmcPiufbKMx8v/le3dLRLRetslJKzOiLxPPkPPR5Zc9gi8ft35X91JlM2mNMS5fKVE0+4REgn7NVq+wHjZdeivaED10qybwOLrmED0vnieGSj+PIh6ux0T925Ig9jnLypB0CF11k3/SPHbNv4uZNv6YGmD/fLilhdmfl28tFy7dhmFL2TLGqqqm/x2xBjY/b57VkiVS61fr77VZRf7/98y+4gGUxyDZvHrB6tf96YCDaK/YzGeA3v7EH8q++mq2W88RwmU5Zmd26GBryu5SUkgFAPVbhOLJgzFyl3N/vfz2VsqdlzpsnLR1zpb35s83Qyi73AthTk81wSaft31NdLUUxzW61kyf9Cz17BtuiRXa4jIz4P0+vkzH/zuZOhESAfC6uucb/zCklN+8ortjXhWS7u/1jy5dLVzLD5bwwXKajByg1PY6ib8zHj/tfKysDFi605/mbM7PGxuzAaGiQG7/ZBTU8LD87u9x+dqFKIPeeLoCEmPl7amvl98Tj/rHBQf8pLXsGW0uL/D3Mv7PZnWeOz1RUyPfzIiSTvm7MMctDh2SsLmqtl+zgc13gbW9jN28AGC7TcRy50epuHz2zCpCbs9lHW1MjrRGztEoy6XchZU81bmyUD7B50x8Z8Ve7Z+/lkt31VFk5dU8X3eIxf8+8efJ75s3zjw0P+/XDzLB0XQlHMzA8z194mT0+U1cnf4iyxWLAO97hf27TaeDpp6M1LVkpmW7/yiv+sZYWGTPiA9N5Y7icSVOTPb6hB+lTKXvsoqFBbvjz5/vH0mn/qd9sLejvLyuTVoU2OioXX/ZeLtmtFGDqni56jCaZtMeF9CSD+nr7e/X3m2ViKiokhBob7RaVDpRk0v47NzdPnWhABPgr9s1dKl9/XdZ/RaX1olst5hji1VdzRX5AGC5nEo/bN+YTJ+TmbY5FAP6geVOTPWNMtwzMm3Is5m8OZo7RjI/71ZDNlkt2KwWY2lWmw0K3fjTdsjAnJkxOyrl7nl0/rLpaLqzaWvsCO37cn2VmTjTIV5yTCJDP7Tve4X9GJiZkgXJ2hYiwJBKyyFOrr7enUdN5YbicSVmZvQhxaEhuzImE3brQg+D19fbT/IkTU8OlokJu3o5jd1dNTkrAmIsiAQm47A98ebndL6y7xfS4DSDvqauT/5rhksnI3yOVsmeK6S60igp7GmYiITeGY8fs8RmGC03HcWSfF3OCyCuvRGNRpVJSE9C8Lq+80r5O6LwwXM7EceQmqo2PywdSP83r79HjFNmD5wMDcjM3B8Xjcf97zDGLdFoG/icn7Zu4+fO0WCx/uJjfo1sg9fV2i+rkSfldZuurqUne47p2oI6OSovIXN9SXm53ARLlUlMjWx5rySTwq1+F33oZHJQuMa22Vs6TpfUDw3/JM9GD+uYA95Ej9o22stJ/4qmosLvRBgflgjJv+nV1cnPWLQvN8+RGnl3l2Bzz0WIxe3rx+Li831x9X1Hhv1f/Tu3kSTk3s/vNDAvzaTOVkjA1qxHU1tp/T6JcHEe6mszP1osvSuWHsFovnifBYlbQuOIK+zNP543hMhPNzfYN/tAhe9ZUba0/dhKL2U3r0VG5kZstBLMVUVtrtyiGh6eGiznor2WXhJmYsEvIAHaNtOpq+++gd780W1/NzfJf3RIzz+vAAXt8pqlpav01olzmzQPe+U7/dTIJ/OIX4bRedKv96af9YzU1cn7s4g0Uw2Umsp/SDxywZ1nNn2+Ps5hPaZOTUvTSbCE0Nfkf5Joauyk+PCwXn/lUl2vMxXHs7rJUyh+o16qr/dZKRYU9eWBw0G59ZYdiU5P983fvtltfixaxC4FmxnFkFpbZMnj55dx7HhWaUsD//q/dTX311XbvBAWCd4eZKC+3LwwdAJq5PbHjSLjo15mMrP41WyLmWpiqKjuYhoZkLMQclM815pIrXLIrHpvBFYvZg/TDwzItVIvH7S66mhr7+/XYkWZOMSU6k7o6YM0a//M4MSF7As3mHkBKyYPes8/6x+rrgXe9i8FSAAyXmXAcqaGV72tmYUhAnvrNacLd3X5Y6OnK+vsrK+3upewpzq6bO1wAu7sslZq6V4vZ5aa7vbTRUbvicX293W2WXfrGVFnJmWJ0dvTYi1kuaM8e4He/m73WSzoNPPaY/QB23XV+dzAFiuEyEzpcskuwAHKDN8ulAPKUZt74zZlflZX29OPsVfp6ZpaWPbZiMteiZDLSl5y9Ot9kzgDTK+615uapVQDy7c/S1MSKsXT24nHgxhv96yiTAR5/3C7YWihKyTToXbv8Y4sWSbiwe7cg+K86U/Pn5556u2jR1BIoVVX558tnB4/r2iExNmavPdHrTrI5jv0+XcrCDAwzXHR3Xb4KxkuW2E9vOlBztZra21l7ic6eLse/apV/7OhR4IknCju4r7eh+K//8q+PWEyCjuWLCobhMlOVlfZFAcjFsnr11Bt2LGavjTG1tNg3Zteduko/e1OxfC2X6mq7MvKRI/bsL72AUtPFMrO5bu5ursZGqfRsqqiw90knOhtlZcD7329/5n/9aykNU6jWSyYDdHXZU+kvvhi47DJ2hxUQw+VsXH21PZDd0SGbCuX6gOYbo8m1J7c5E2183J5/X1mZuzsOkKAwm/TmxZMdWoC0QnK1vqqrcx93XRmE1S0gPesnu5VDNFOOI2N573mP/xmanAR+8hN7K4igKCUz08xB/Lo64Kab2PousDx3LZpCtwQ+9jHgtdfkhr9qVe4FjnqQPx6fuulXrn0izO4rz7NL5mcHiEkHj27qmwOV5eVTWymxGNDWJlNATYsW5S7Wp/8eH/84sHevtHze8hb2UdP5cV1ZV7J7t3yuAGl1/+d/An/yJ8EVQ9Wt+Ucf9cc9XRe44QY+IM0C3iXOhuNIK+Ptb/erp+b7gDY2ytiEafHiqYP/ur5Yvp+TvQ7GVFWV/+mrsjL3eMnKlXY3m+NIN1e+sRhd/ub66+X7Kip4UdL5q6oCbrnFHvN46SXg5z8Ppiy/XpD8wx/a3cyrVnEQf5bwX/hc6FXs03FdafrraY51dfLElOuprK4u/809e9zEVF6eu+UESChl/y4dFNddJ+91Xel7futbz/z3mcnfmWimdKt4/Xq/29fzZPbY88+f3wC/3g/pxz+WDcq0+fOBD3+Y20TMEnaLFYquSfbJT8pYSFOTvb7FpIPAnLKsZU8nNpWVyXvNagFavsByXeAP/kC6t1Ip6R7gxUZhcF3ZDrmvT1bNAzL+sn27PPxcdtnZtzB0sDzyiNQw0+JxYMMGe4EzFRTDpZB0l9d0AQFI66O62h4z0aYrDum6+adSNjTkvogcR0LJXMxGFJZYTAbXBwaAV1+VY8kk8NBDEhLXXCOf85kEglKyTuyRR4DnnvMnB5SVSQtp5UoGyyxit1gUlJfnDhHXnX48xnHyL2Y0V+MTRZUuY3TrrfYYZTIJPPwwsGOHXQ4pH8+TFtD3vy/7tJgVMW68keMsIeC/dhToveuz6W2Hp5MrRLLrmxFFmZ4oc8cdMr1fS6VkgeUDD8jqer1nkaaUv0Pq448D998PvPGG/3UdLO97H4MlBOwWi4pciy5raqbfz1vXC3NdewC0sjJ3WBFFld4t9c47pcXyyit+kBw8CDz4oFwjF10ksy4rKmQ2WE+PBIpZ1QKQa+Dmm2WWY751YlRQ/FePAseRgfWyMntQf8GCMw+2NzfLeI1Zj6yhgWUtqPjoMcrbb5cWyy9+4RdizWSkovHhw2f+OQsWyDTnVavYYgkRwyUq5s+XoDBX2V944Zm7tmprJZj27PGPtbdzBhgVJ8eRCS433SStlK4umU48k7Uv8bhMAHjve/NPaKFZw3CJiqoq2Vfipz+Vp7Vly/KXljHp1c59fTJTZv584NprZ+eciQrFdeXhavlyWcX/299KZYmREb917zjyENXUJK2Ua67xN/1isISO4RIVjgNceqk06YeHZapwriKTud63YoWspxkYkH7p6WaYERULHR6rVsli39FR+YwPDUnA6OrjjY1+S52f+8hguBgcpZA4eBAZc0+UsBw8eG7v6+uzNwGbRZ7nwUsmwcu7uDlKYaCnJxrXQT56fPL4cfkTIRnPg5drzVqJYbic4jgOsGAB1J49GDLHL+isxCor4eYrSUOR556qWqxeew2J114L+3SKVllVFWL5dpAtEY5Ss7XHaLQppeB5HvjPcf5c14XLWTpFiddBcEr9OmC4EBFR4Eo3VomIqGAYLkREFDiGCxERBY7hQkREgWO4FItMRhZXBrEFLFGx4nVQNBguxeLoUeDrX7drjxGVmjffBP76r+W/FGkMFyIiChzDhYiIAsdwISKiwDFciIgocAwXIiIKHMOFiIgCx3AhIqLAMVyIiChwDBciIgocw4WIiALHcCEiosAxXIiIKHAMFyIiChzDhYiIAsdwKQJKKQwMDCCRSGBgYABKqbBPiWjW6etgaGiI10ERYLhEWCKRwObNm9HZ2Ym3rl6Nb99/P966ejU6OzuxefNmJBKJsE+RqODM6+DSyy7D9x98EJdedhmvg4hzFOM/krq6urBhwwaMjY0BABYphT8D8G0ARxwHAFBdXY1t27Zh3bp14Z0oUQFlXwdLlML/AfD/APTxOog0tlwiqKurC+vXr0cymYRSakrzXx9LJpNYv349urq6QjpTosLhdVDc2HKJmEQigdbWViSTSXied/r4YuB0y8Xc4NV1XcTjcfT29qKhoWF2T5aoQPJdB0uB0y2Xw8b38zqIHrZcImbLli0YGxuzLqjpeJ6HsbExbN26tcBnRjR7eB0UP4ZLhCil8K1vfeuc3nvfffdx9gzNCbwO5gaGS4T09/eju7v7rC8OpRS6u7sxMDBQoDMjmj28DuYGhkuEjIyMnNf7h4eHAzoTovDwOpgbGC4RUltbm/drxyGD+ceneX9dXV3Qp0Q066a7Do5BBvOPTfN+XgfRwHCJkObmZnR0dMA5NX/flIbMEkvneJ/jOOjo6EBTU1OhT5Go4Ka7DlKQWWKpHO/jdRAtDJcIcRwHd9999zm9d+PGjTkvRqJiw+tgbuA6l4jJN78/H87vp7mI10HxY8slYhoaGrBt2zY4jgPXnf7/Htd14TgOtm/fzguK5hReB8WP4RJB69atw44dOxCPx+E4zpRmvj4Wj8exc+dOrF27NqQzJSocXgfFjeESUevWrUNvby82bdqE9vZ262vt7e3YtGkTDh8+zAuK5jReB8WLYy5FQO9jMTw8jLq6OjQ1NXHQkkoOr4PiwnAhIqLAsVuMiIgCx3AhIqLAMVyIiChwDBciIgocw4WIiALHcCEiosAxXIiIKHAMFyIiChzDhYiIAsdwISKiwDFciIgocAwXIiIKHMOFiIgCx3AhIqLA/X/4lPfT9Tg3VgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 500x400 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "e26b771f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "function , r2\n",
      "exp , 0.9999988610586863\n",
      "cosh , 0.9999699077016541\n",
      "sigmoid , 0.9999693609882967\n",
      "arctan , 0.9999174139339265\n",
      "gaussian , 0.9999096961395885\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "('exp',\n",
       " (<function kan.utils.<lambda>(x)>, <function kan.utils.<lambda>(x)>),\n",
       " 0.9999988610586863)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.suggest_symbolic(1,0,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "b939a769",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "r2 is 0.9999988610586863\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "tensor(1.0000, grad_fn=<SelectBackward0>)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.fix_symbolic(1,0,0,'exp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "a0e2813a",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "train loss: 8.09e-04 | test loss: 8.51e-04 | reg: 4.68e+01 : 100%|██| 20/20 [00:05<00:00,  3.96it/s]\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 1.0 e^{1.0 x_{2}^{2} + 1.0 J_{0}{\\left(- 20.0 x_{1} \\right)}}$"
      ],
      "text/plain": [
       "1.0*exp(1.0*x_2**2 + 1.0*J0(-20.0*x_1))"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# why can't we reach machine precision (because LBFGS early stops?)? The symbolic formula is correct though.\n",
    "model.train(dataset, opt=\"LBFGS\", steps=20);\n",
    "model.symbolic_formula()[0][0]"
   ]
  }
 ],
 "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
}