GitHub_collection_pykan/tutorials/Example_9_singularity.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "134e7f9d",
   "metadata": {},
   "source": [
    "# Example 9: Singularity"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2571d531",
   "metadata": {},
   "source": [
    "Let's construct a dataset which contains singularity $f(x,y)=sin(log(x)+log(y))\n",
    " (x>0,y>0)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2075ef56",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cuda\n",
      "checkpoint directory created: ./model\n",
      "saving model version 0.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "| train_loss: 1.14e-01 | test_loss: 1.29e-01 | reg: 6.34e+00 | : 100%|█| 20/20 [00:03<00:00,  5.03it"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "saving model version 0.1\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "from kan import *\n",
    "import torch\n",
    "\n",
    "device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n",
    "print(device)\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=5, k=3, seed=2, device=device)\n",
    "f = lambda x: torch.sin(2*(torch.log(x[:,[0]])+torch.log(x[:,[1]])))\n",
    "dataset = create_dataset(f, n_var=2, ranges=[0.2,5], device=device)\n",
    "\n",
    "# train the model\n",
    "model.fit(dataset, opt=\"LBFGS\", steps=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3f95fcdd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAFICAYAAACcDrP3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAtEElEQVR4nO3deVBUV74H8O9pmn0RRFwQdWjouGMGF1AERlEwMTqOZqLJxHrGOJNRg5Xt6cSYp/FpzExMiVuciXmT0SQvasRoFNxGfYC7QXGLoogoS0DAbqRplqb7vj8iXWBcUC709v1UpaaK27fvD8bT3z7n3HuOkCRJAhERkYwUli6AiIjsD8OFiIhkx3AhIiLZMVyIiEh2DBciIpIdw4WIiGTHcCEiItkxXIiISHYMFyIikh3DhYiIZMdwISIi2TFciIhIdgwXIiKSHcOFiIhkx3AhIiLZKS1dAJEtkCQJ5eXl0Ol08PLygr+/P4QQli6LyGqx50L0EFqtFitXroRarUZAQACCg4MREBAAtVqNlStXQqvVWrpEIqskuBMl0f3t3bsXkyZNgl6vB/Bz76VBQ6/Fw8MDycnJSEhIsEiNRNaK4UJ0H3v37sXYsWMhSRJMJtMDX6dQKCCEQEpKCgOGqBGGC9E9tFotgoKCUF1d/dBgaaBQKODu7o6CggL4+vq2foFENoBzLkT32LBhA/R6fbOCBQBMJhP0ej02btzYypUR2Q72XIgakSQJarUaubm5eJymIYSASqXC1atXeRcZERguRE2UlZUhICCgRef7+/vLWBGRbeKwGFEjOp2uRedXVlbKVAmRbWO4EDXi5eXVovO9vb1lqoTItjFciBrx9/dHSEjIY8+bCCEQEhKC9u3bt1JlRLaF4ULUiBACiYmJT3TunDlzOJlPdBcn9InuwedciFqOPReie/j6+iI5ORlCCCgUD28iDU/ob9u2jcFC1AjDheg+EhISkJKSAnd3dwghfjHc1fAzd3d3pKamIj4+3kKVElknhgvRAyQkJKCgoABJSUlQqVRNjqlUKiQlJaGwsJDBQnQfnHMhagZJknDo0CHExcXhwIEDGDFiBCfviR6CPReiZhBCmOdUfH19GSxEj8BwISIi2TFciIhIdgwXIiKSHcOFiIhkx3AhIiLZMVyIiEh2DBciIpIdw4WIiGTHcCEiItkxXIiISHYMFyIikh3DhYiIZMdwISIi2TFciIhIdgwXIiKSHcOFiIhkx3AhegSDwYDCwkJcunQJAHDt2jXcvn0bJpPJwpURWS9uc0z0AFqtFsnJyfj6669x8eJFVFZWoq6uDm5ubggICEB0dDReffVVREVFQalUWrpcIqvCcCG6j2PHjuHNN9/EuXPnMHjwYIwdOxZhYWHw8vKCVqtFZmYmdu7ciZycHEyePBlLlixBQECApcsmshoMF6J77Nu3D9OmTYOXlxeWLVuGZ599FnV1ddi0aRNqa2vh4+ODKVOmwGAwYNOmTVi0aBH69u2LL7/8Ep06dbJ0+URWgeFC1MiVK1cwZswYeHp6YtOmTejTpw+EEMjNzUV4eDgqKioQHByMzMxM+Pn5QZIkHD58GC+99BJ+85vf4PPPP4erq6ulfw0ii+OEPtFdRqMRH374ITQaDdasWWMOlocRQmD48OH429/+hh07dmDPnj1tVC2RdWO4EN2Vk5ODnTt3YuLEiRg+fPgjg6WBEAITJkxAZGQk1q9fj/r6+laulMj68RYXoruOHj0KnU6HSZMmIS8vD1VVVeZjBQUFMBqNAIC6ujpcvHgRPj4+5uOBgYGYOHEiFi1ahOLiYgQFBbV5/UTWhOFCdNfly5fh4eEBlUqF1157DUeOHDEfkyQJtbW1AICioiKMHj3afEwIgU8++QT9+/eHXq9HUVERw4UcHsOF6K7q6moolUq4urqitrYWNTU1932dJEm/OFZfXw93d/cmIUTkyBguRHd17NgR1dXV0Gq1iIiIgKenp/lYdXU1jh49ag6RYcOGmR+cFEKge/fuuHXrFhQKBfz8/Cz1KxBZDYYL0V0DBw6EwWDAyZMn8de//rXJsdzcXAwePBgVFRXo1KkTNm/eDF9fX/NxIQTmz5+Pzp07c0iMCLxbjMhsyJAhUKlU2LBhA6qqquDk5NTkvwZCCCgUCvPPFQoFfvrpJ2zduhVjx45Fu3btLPhbEFkHhgvRXf7+/nj99ddx+vRprFq1qtm3FNfW1uK///u/UV1djddee63ZtzAT2TMOixE1Mm3aNKSnp+Ovf/0rPDw8MHPmTLi5uQEAlEollEqluRcjSRIqKyuxdOlSbNq0CStWrEDPnj0tWT6R1eDyL0T3KC0txezZs7Fr1y4kJCTgzTffRO/evZGdnQ2TyQQXFxeEhobi5MmTWL58ObKysrB48WLMnDmzyfAZkSNjuBDdR1VVFdavX49Vq1ahpKQEKpUKarUa3t7e0Gg0yM7ORlFREQYOHIiFCxciNjYWCgVHmYkaMFyIHqK4uBgHDhxAWloazp49i5MnTyI6OhpRUVGIj49HREQEPDw8LF0mkdVhuBA106lTpzBkyBCcOnUKgwYNsnQ5RFaN/XiiZnJycjLfhkxED8dWQkREsmO4EBGR7BguREQkO4YLERHJjuFCRESyY7gQEZHsGC5ERCQ7hgsREcmO4UJERLJjuBARkewYLkREJDuGCxERyY7hQkREsmO4EBGR7LifC1EzSZIEk8kEhUIBIYSlyyGyauy5ED0G7uVC1DxKSxdAJAeDwYCbN2/CZDJZupQWE0Kge/fucHFxsXQpRE+M4UJ2oaCgALNmzcLAgQMtXcoTkSTJPNSWmZmJTz/9FCEhIRauiujJMVzILkiShLCwMCxdutTSpTyRdevW4cSJE4iPj0d9fT04FUq2juFCdsfWJttNJhMOHDiA7du3w9PTE15eXpYuiajFODtJZGF1dXW4efMmAEClUtlcOBLdD8OF6D4kSYJOp8PZs2dRXl7eqsNUOp0Ot27dghACwcHBrXYdorbEcCG6j9u3b2PKlCkYNmwYRo8ejcuXL7dawJSXl0Or1cLZ2Rndu3dvlWsQtTWGC9E9JEnCp59+ij179qCmpgZnz57FggULYDAYWuV6N2/eRHV1Nby9vdGlS5dWuQZRW2O4EN2jrKwMGzduhCRJiIyMhFKpxL59+3DmzBnZryVJEnJycmA0GtGpUyf4+fnJfg0iS2C4EDUiSRIOHz6MGzduwNfXF0lJSejbty/0ej2+/fbbVhkay87OBgD06NEDbm5usr8/kSUwXIgakSQJqampMBqNGDRoEJ5++mlMmDABALB3717odDpZr2cymZCTkwMACA0N5fIyZDf4L5moEZ1OhyNHjgAAEhISoFQqMWbMGLi5uSE3NxeXLl2S9Xo1NTXIy8sDAPTq1UvW9yayJIYLUSPXrl3DzZs34erqiuHDh0MIgV69eiE4OBg1NTXIyMiQdWhMo9Hg1q1bcHJy4nIvZFcYLkR3SZKEzMxMVFdXo2vXrlCr1QAAb29vREREAADS0tJkXRyzqKgIlZWVcHd3523IZFcYLkSNHDt2DADQv39/+Pj4mH8eGxsLIQTOnTsHjUYj2/Vyc3NhMBjQvn17BAQEyPa+RJbGcCG6q+GZFgAYMmSIeXJdCIGBAwfCw8MDJSUluHr1qizXkyQJV65cgSRJCAoK4ppiZFcYLkR33bp1Czdu3ICTkxPCw8ObrPHVvXt3dOvWDXV1dTh16pRs8y4NtyGrVCo4OzvL8p5E1oDhQnTX1atXUVFRAR8fHzz11FNNjnl6emLAgAEAgOPHj8sSLnV1dbh+/ToAoGfPni1+PyJrwnAhws9DVOfOnYPRaES3bt3QsWPHJseFEOZJ/fPnz0Ov17f4mjqdDoWFhRBCQK1WczVksisMF6K7srKyAPz8vMm9T8o3zLs4OzsjPz8fhYWFLb7erVu3oNFo4OLiwtWQye4wXIgA1NbW4vLlywCAAQMG3LcXERoaCn9/f+h0Oly8eLHF17xx4waqq6vh4+ODwMDAFr8fkTVhuBDh5yX28/PzoVAo0K9fv/uGS/v27REaGgpJklo8qS9JErKzs2EymdC5c2f4+vq2oHoi68NwIcLPy95rtVp4eHggNDT0vq9xdnZGeHg4ACAzMxNGo7FF1zx//jwAQK1Wc8FKsjsMF3J4Db2Iuro6dOzYEZ07d37ga4cMGQIhBLKzs6HVap/4mgaDwTwM17dvX07mk91huBABuHDhAgAgODj4gQ8zCiHQr18/uLu7o7S0FLm5uU98PY1Gg7y8PAgh0L9/f4YL2R2GCzk8k8lkXu24d+/ecHJyeuBru3fvjsDAQNTW1iIrK+uJ511u3LiB27dvw8PDg6shk11iuJDD0+v1uHbtGgCgX79+D32tl5cX+vbtCwA4ceLEE11PkiT8+OOPqK2tRadOndC1a9cneh8ia8ZwIYdXWlqKkpISKJVK9OzZ86FDVAqFwvww5ZkzZ1BdXf1E1zx9+jSAnyfzuaYY2SOGCzm8vLw86HQ6eHt7o0ePHg99rRACgwYNglKpRF5eHoqKih77evX19eY5nrCwsIcOwxHZKoYLOTRJknD58mUYjUZ06dIFHTp0eOQ5vXv3hr+/PyorK3Hu3LnHnnfRarXIycmBEOIXC2QS2QuGC9kdSZIe6wO/oRcREhLSrOdNAgIC0KtXL5hMJvOWyI8jJycHpaWl8PDwMM/fENkbhgvZFb1ej3/+85/QaDTNCpj6+nrz8yZ9+vQx7+HyMEqlEkOHDgXw8wrJdXV1za6vYbfLuro6BAUFcfdJslsMF7IbBoMBc+fOxaxZszB79mxUVVU98hydToe8vDwAzX+YUQiB4cOHw8nJCdnZ2Y+1iKUkSebdLgcMGABPT89mn0tkSxguZDecnJwQEhIChUKBrVu3Yu3atY/c7764uBhlZWVwcXH5xR4uDxMWFgZ/f39otVpkZmY2exiuqqrKvNvl0KFDOd9CdovhQnZDoVDg9ddfxx//+EeYTCYkJSU9ckviq1evorq6Gn5+fujWrVuzr9WxY0eEhYVBkiQcOHCg2efl5eXh5s2bcHFxwaBBgxguZLcYLmRXnJ2dMW/ePKhUKty6dQufffbZA3svDQ8zmkwmdOvW7bFWJlYqlRg5ciQA4MiRI9DpdI88R5IknDhxAlVVVejSpQt3nyS7xnAhuxMYGIhXX30VALB161b89NNP931dw+6TwM/bDLu6ujb7GkIIjBgxAm5ubsjNzTUvH/MwkiTh4MGDAIDw8HAus092jeFCdkcIgcmTJ6Njx44oLCxESkrKfedEampqmmwQ9rh69+4NlUqFmpoa7N+//5HzLhqNxrxkTFxcXLPuTCOyVfzXTXape/fuiI+PhyRJ2Lx5831vFy4rK0NBQQGcnJweuEHYw3h5eSEuLg4AkJqaipqamge+VpIknD17FoWFhfDy8kJMTAznW8iuMVzILikUCkyZMgVKpRI//PCDuYfSWE5ODioqKuDt7Q21Wv1E1xk3bhycnZ1x/vz5Rw6N7d69GwaDAb169YJKpXqi6xHZCoYL2SUhBIYOHYqQkBDodDrs2LGjybCVJEnIyspCfX09goKC0KlTpye6xqBBg6BWq1FVVYXvvvvugUNjd+7cwZ49ewAAzzzzDHeeJLvHcCG71a5dOzz33HMAgO3btzd5qFKSJJw8eRIAzBuAPQkfHx+MHz8eAJCcnAyNRvOL1zRM5GdnZ8PLywvjx4/nkBjZPYYL2S0hBCZOnAg3NzdcunSpycOOOp3OfKdYZGTkE3/YCyEwZcoU+Pj4ICcnB/v27ftF76W+vh5ffPEFjEYjYmJiHrlnDJE9YLiQXQsLC8OAAQNQV1eHzZs3mz/4r1+/jvz8fLi4uGDw4MEt6kn06tULcXFxMBqNWLduHfR6vfmYJEm4cOEC0tLS4OTkhFdeeQXOzs4t/r2IrB3Dheyau7s7fv/73wMAUlJSUFxcDEmScPToUej1enTt2vWxln25H6VSiVmzZsHd3R3Hjx9HcnKyOcQMBgNWrlwJnU6HXr16YeTIkRwSI4fAcCG7JoTAhAkTzM+8JCcnw2g0Yvfu3QCAiIiIFj/M2LCQ5bhx41BfX49Fixbhxx9/hCRJ2LFjB7Zs2QInJye8/vrraNeunQy/FZH1Y7iQ3evRowd+97vfQZIkfPrpp8jIyMDRo0ehUChkm1x3dnbGokWL0L17d9y4cQPPP/883njjDSQmJqK2thbx8fF4+eWX2Wshh8FwIbsnhMDs2bPRoUMHXLlyBX/4wx+g0WigUqlkG6YSQuCpp57CunXr0LlzZ1y5cgVr1qxBaWkp+vfvj08++eSJ70gjskUMF7J7Qgj07t0b7777LlxcXFBSUgIXFxe88847zdrW+HGuk5CQgNTUVMyYMQPx8fGYO3cudu7ciaeeeoq9FnIoSksXQCS3+z3IKITArFmz0LFjR+zfvx9xcXF44YUXHvj6lggLC8Pf//53mEymJuuHyX0dImsmJP6LJztw/fp1zJw5E5GRkZYupcWOHz+OtWvXIiQkxNKlED0xhgvZhbq6OuTm5sJoNFq6lBZTKBQICQmBi4uLpUshemIMFyIikh3nXIiaqfH3ME7OEz0c7xYjaqYzZ87AyckJZ86csXQpRFaP4UJERLJjuBARkewYLkREJDuGCxERyY7hQkREsmO4EBGR7BguREQkO4YLERHJjuFCRESyY7gQEZHsGC5ERCQ7hgsREcmO4UJERLJjuBARkewYLkREJDuGC1EzSJIEjUYDANBoNOAGrkQPx3AhegitVouVK1dCrVZj1KhRkCQJo0aNglqtxsqVK6HVai1dIpFVEhK/ghHd1969ezFp0iTo9XoA99/m2MPDA8nJyUhISLBIjUTWiuFCdB979+7F2LFjIUkSTCbTA1+nUCgghEBKSgoDhqgRhgvRPbRaLYKCglBdXf3QYGmgUCjg7u6OgoIC+Pr6tn6BRDaAcy5E99iwYQP0en2zggUATCYT9Ho9Nm7c2MqVEdkO9lyIGpEkCWq1Grm5uY91R5gQAiqVClevXjXPxxA5MoYLUSNlZWUICAho0fn+/v4yVkRkmzgsRtSITqdr0fmVlZUyVUJk2xguRI14eXm16Hxvb2+ZKiGybQwXokb8/f0REhLy2PMmQgiEhISgffv2rVQZkW1huBA1IoRAYmLiE507Z84cTuYT3cUJfaJ78DkXopZjz4XoHr6+vkhOToYQAgrFw5tIwxP627ZtY7AQNcJwIbqPhIQEpKSkwN3dHUKIXwx3NfzM3d0dqampiI+Pt1ClRNaJ4UL0AAkJCSgoKEBSUhJUKlWTYyqVCklJSSgsLGSwEN0H51yImkGSJBw6dAhxcXE4cOAARowYwcl7oodgz4WoGYQQ5jkVX19fBgvRIzBciIhIdgwXIiKSHcOFiIhkx3AhIiLZMVyIiEh2DBciIpIdw4WIiGTHcCEiItkxXIiISHYMFyIikh3DhYiIZMdwISIi2TFciIhIdgwXIiKSHcOFiIhkx3AhIiLZMVyIHsFgMKCwsBCXLl0CAFy7dg23b9+GyWSycGVE1ovbHBM9gFarRXJyMr7++mtcvHgRlZWVqKurg5ubGwICAhAdHY1XX30VUVFRUCqVli6XyKowXIju49ixY3jzzTdx7tw5DB48GGPHjkVYWBi8vLyg1WqRmZmJnTt3IicnB5MnT8aSJUsQEBBg6bKJrAbDhege+/btw7Rp0+Dl5YVly5bh2WefRV1dHTZt2oTa2lr4+PhgypQpMBgM2LRpExYtWoS+ffviyy+/RKdOnSxdPpFVYLgQNXLlyhWMGTMGnp6e2LRpE/r06QMhBHJzcxEeHo6KigoEBwcjMzMTfn5+kCQJhw8fxksvvYTf/OY3+Pzzz+Hq6mrpX4PI4jihT3SX0WjEhx9+CI1GgzVr1piD5WGEEBg+fDj+9re/YceOHdizZ08bVUtk3RguRHfl5ORg586dmDhxIoYPH/7IYGkghMCECRMQGRmJ9evXo76+vpUrJbJ+vMWF6K6jR49Cp9Nh0qRJyMvLQ1VVlflYQUEBjEYjAKCurg4XL16Ej4+P+XhgYCAmTpyIRYsWobi4GEFBQW1eP5E1YbgQ3XX58mV4eHhApVLhtddew5EjR8zHJElCbW0tAKCoqAijR482HxNC4JNPPkH//v2h1+tRVFTEcCGHx3Ahuqu6uhpKpRKurq6ora1FTU3NfV8nSdIvjtXX18Pd3b1JCBE5MoYLOby8vDykpaXhyJEj0Ov10Gq1iIiIgKenp/k11dXVOHr0qDlEhg0bZn5wUgiB7t2749atW6ivr0dOTg4GDx4MNzc3S/1KRBbHW5HJ4eTn5yM9PR1paWlIS0vDzZs3zQGRk5ODtWvXYsaMGU3Oyc3NxeDBg1FRUYFf/epX+OGHH+Dr62s+LoTA/PnzsXz5ciiVSri5uSEiIgIxMTGIjY3F4MGDeYsyORSGC9m9oqIipKWlmQPl+vXrAICwsDDExMQgJiYG0dHRMBqNGD58OPz8/LBnz54mE/YPes4F+HmYrKioCLGxsRg3bhymTZuGjIwMpKenIz09HRUVFXB3d0dkZKQ5bAYOHAgXFxeL/D2I2gLDhexOSUlJkzDJyckBAPTt29f84R4dHY327dv/4ty1a9fi7bffxoIFC/CXv/zFPPT1sHCpqanBG2+8gZ07d+LgwYPo2bOn+f2MRiPOnz9vruXw4cO4c+cOPDw8MHToUMTGxiI2Nha//vWv4ezs3AZ/HaK2wXAhm1daWmruJaSlpSE7OxsA0KtXryZh0py1v6qqqjB9+nSkpqbigw8+wMyZM+Hm5obr169jyJAh5mGxkydPwtfXF5WVlVi6dCn+8Y9/YMWKFXjllVce+v719fU4e/asOfyOHDkCnU4HLy8vDBs2zBw2AwYM4GKYZNMYLmRzysvLzcNOaWlp+PHHHwEAarXaHCYxMTFPvM5XaWkpZs+ejV27diEhIQFvvvkmevfujezsbJhMJri4uCA0NBQnT57E8uXLkZWVhcWLF2PmzJlwcnJ6rGsZDAacOXPGHDZHjx6FXq+Hj48PoqKizGHTv3//x35vIktiuJDV02q1yMjIMH8Anz9/HgCgUqmahElgYKBs16yqqsL69euxatUqlJSUQKVSQa1Ww9vbGxqNBtnZ2SgqKsLAgQOxcOFCxMbGQqFo+YIXdXV1yMzMNAfnsWPHUFNTg3bt2iE6Otr8u/br10+W6xG1FoYLWZ2KigocPnzYHCbnzp2DJEno0aOH+cM1Nja2TR5ULC4uxoEDB5CWlobc3FzU1NTAz88P/fr1Q3x8PCIiIuDh4dFq16+trcWpU6fMYXPixAnU1tbCz8/PHDaxsbHNWgeNqC0xXMjiKisrceTIEfMHaFZWFkwmE7p27Wr+8IyNjUWPHj0sWqfRaIQkSVAoFBbrNdTU1ODEiRPmv9XJkydhMBjQoUOHJmHTs2dPhg1ZFMOF2pxOp8OxY8fMPZPTp0/DaDSiS5cuTXomwcHB/IB8BL1ej+PHj5vD5ocffkB9fT06duxo/jvGxsYiNDSUf0tqUwwXanUNH4ANYdL4A7BxmPADsOV0Op35b52WltYkuBuHDYObWhvDhWTXMHTTECaNh24aHlrk0E3buHPnjrmX2HjIMSgoqEnYWHrIkewPw4VarGHSuSFMGk86Nw4TTjpbXsPNEg3DaI1vlmgcNlzVmVqK4UKPreF22YZvw8ePHzffLtsQJrxd1jZoNJomYdNwm3dwcHCTmym6dOli4UrJ1jBc6JEaP+jX8OxFw4N+w4cPN8+b8EE/21deXm6+DbzxA6qhoaHmoGnJA6rkOBgu9AuNlyhJS0sz79Do5eWFqKgo8/AJlyixf6WlpeYHWBsvrdOzZ88mYdOhQwcLV0rWhuFC5sUVGz5Ajhw5Yl5csWG9q5iYGC6uSCguLm4SNg2Lgvbp08ccNg9aFJQcC8PFAZlMJly4cME8zp6RkYGKigq4ubmZV+qNiYnhsvD0SEVFRU32xmnYzqB///7mf0fR0dFN9r4hx8BwcQCSJOHHH380fwBkZGRAo9HA1dUVERER5g8BbmhFLfWgjdjCwsLMPZuoqCi0a9fO0qVSK2O42CFJkpCdnd0kTMrKyuDs7IwhQ4aYwyQiIoJb8VKrysvLaxI2hYWFUCgUePrpp81hM2zYMHh7e1u6VJIZw8UOSJKEnJwccwNOT0/HrVu3oFQqMWjQIHOYREZGtuoii0QPI0kSrl+/bv53mpaWhuLiYjg5OSE8PNwcNkOHDoWnp6ely6UWYrjYoHsbaXp6On766adfNNLIyEh4eXlZulyi+2rOl6LY2NhWX3maWgfDxUbcuHGjSSMsKCj4xfDC0KFDm+z7TmRL7h3OTU9PR3l5OVxcXDB48GAO59oYhouNGDBgAK5evcqJUXIYJpMJly5d+sWNKF9++SWef/55S5dHj8BwsREmkwlCCK7NRQ5LkiRIksR2YCMYLkREJDuu3SEDg8GA/Px8mEwmS5fSYkIIdOvWjQ9P0mNhG6B7MVxkUFhYiDlz5iA8PNzcdbfV1YBPnz6NVatWQaVSWboUsiGFhYVITEy0mzawevVqtoEWYrjIQJIk9OvXD+Hh4di2bRtGjBiBadOmWbqsJ7JgwQJwpJQeV0MbGDBgAHbt2oWoqCjMmDHD0mU9kffee49tQAYMFxnt3bsXW7duRUVFBV5++WWbW+SRDYpaQgiB//u//8PmzZtRUVGB6dOn29wWDGwD8rHNfqsVEkIgLi4OQgicPXsWJSUlli6JqM0NGTIEAHDp0iVUVFRYuBqyJIaLjIYMGYL27dujrKwMp06d4rcgcjgDBgyAu7s7SkpKkJuba+lyyIIYLjLq0qULnn76aZhMJuzbt4/hQg6nR48eCAwMRE1NDU6fPs024MAYLjJSKpUYPXo0AJj3SCFyJN7e3ggLCwMAHDlyhOHiwBguMhsxYgQ8PT2Rn5+Ps2fPWrocojalUCgQFRUF4OdbeisrKy1cEVkKw0VmoaGh6N27NwwGA1JTU/nNjRxOZGQk3NzcUFBQYN4GmRwPw0Vmbm5uGDNmDADg3//+N+7cuWPhiojaVmhoKLp164aamhocPnyYX7AcFMNFZkIIjBkzBh4eHrh27RrOnDnDxkUOxdvbG5GRkQCAgwcPwmg0WrgisgSGSyvo3bs3+vXrB4PBgB07dli6HKI2JYTAqFGjIIRAVlYWiouLLV0SWQDDpRW4ublh/PjxAH5+ar+8vNzCFRG1HSEEIiMj0aFDB5SVleHYsWPsvTsghksrEELg2WefRbt27ZCfn4+MjAw2LnIogYGBGDRoEEwmE3bt2mUXqyXT42G4tJKQkBAMGzYMRqMRmzZt4rgzORQnJyc899xzEEIgIyODyyE5IIZLK1EqlZgyZQoUCgUyMjJ4SyY5lIa19vz9/VFcXIxDhw6x9+5gGC6tRAiBESNGIDg4GFqtFt9++y0bFzmUoKAgREdHQ5IkbN68GQaDwdIlURtiuLQif39/TJo0CQCwZcsWlJWVWbgiorajUCjw4osvQqlU4tixY7h48aKlS6I2xHBpRUIITJkyBX5+frh+/Tp27drF3gs5DCEEYmJi0LNnT+h0Onz55Zec2HcgDJdWplar8cwzz8BkMuHzzz+HTqezdElEbcbHxwdTp06FEALbtm1DXl6epUuiNsJwaWUKhQJ/+tOf4OnpiXPnzmH37t3svZDDEELg97//Pbp164aSkhKsX7+evRcHwXBpZUII/PrXv8bo0aNhNBqxevVqrhRLDqVLly6YPn06hBD46quvcOnSJX7BcgAMlzagVCoxZ84ceHp6IisrC1u3bmXjIochhMC0adOgVqtRVlaGDz/8EHV1dZYui1oZw6UNCCEwaNAgTJgwAUajEUlJSVxviRxKx44dMW/ePDg7O2Pnzp28Nd8BMFzaiJOTE9566y0EBAQgJycHq1at4lP75DCEEJg4cSKee+45GAwGLFy4EOfPn2fA2DGGSxsRQqBnz56YNWsWhBD45z//yTXHyKG4urpi6dKlCAkJQVFRERITE1FcXMw2YKcYLm2o4c6xwYMHo7KyEvPnz0dJSQkbFzkEIQR+9atfYcWKFWjXrh1OnjyJxMREaDQatgE7xHBpY76+vli2bBn8/PyQlZWFhQsXora21tJlEbWJhjXHli5dCjc3N6SmpuL1119nwNghhksbE0IgIiICf/nLX6BUKvHNN99gzZo1qK+vt3RpRG1CoVDgP/7jP/D+++/DxcUF27dvx/Tp01FYWMiAsSMMFwtQKBT44x//iJdffhlGoxHLli3DV199xYfLyGEolUokJibigw8+gJubG/bu3YtJkybhxIkTbAd2guFiIa6urliyZAni4+NRXV2NefPm4auvvuIdZOQwnJ2dMXv2bKxevRodOnTAuXPn8Pzzz2PNmjWoqqpiL8bGMVwsRAgBPz8/fPrpp4iOjkZlZSXefvttrFu3jg+YkcNQKpV46aWXsGXLFoSFhaG8vBzz58/H5MmTkZmZyV6MDWO4WJAQAp07d8a//vUvxMXFoaqqCgsWLMD8+fNRUVHBb27kEIQQiIyMxI4dO/CnP/0JLi4uOHDgAMaNG4eFCxfyjkobxXCxsIaA+eKLLzB58mQYjUasW7cOf/jDH3DlyhU2KnIIDe3gk08+wf/+7//i6aefRkVFBZYvX46EhAR8++23qKmpYXuwIQwXKyCEgL+/P9auXYt58+bB3d0dBw8exPjx4/Hdd9/BYDCwUZFDUCqVSEhIQEpKCv7rv/4LAQEByM7OxowZM/Dyyy8jKyuLQ2U2guFiJYQQ8PDwwLvvvosvvvgCPXr0QH5+PmbMmIH//M//5NAAOQwhBNq3b4+5c+di9+7d+N3vfgchBFJSUjB27FgsWLAABQUFbA9WjuFiZZycnPDcc8/h+++/xzPPPAODwYDPPvsM48aNQ2pqKnsx5DAUCgX69OmDf/3rX9i4cSP69+8PrVaLFStWYNSoUVi1ahVu3brF9mClGC5WSAiB0NBQbNy4EUuWLEGHDh1w8eJFTJ06FX/+85+RnZ3NoQFyCEIIuLi4YPz48di9ezc++OADBAYG4saNG3j33XcxcuRIfPzxx7h+/TqMRiODxoowXKyUEAKenp5ITEzE999/j9GjR6O+vh7ffPMN4uPj8f777+PatWswmUxsUGT3GuYl33nnHfz73/9GYmIi/P39ce3aNSxcuBAxMTGYPn06du/ezaVkrATDxcopFAoMGDDAvEyMWq1GeXk5VqxYgREjRmD27NlIS0vDnTt3IEkSGxXZNSEEgoOD8dFHH+HQoUN455130KNHD9y+fRtbtmzBCy+8gOjoaMybNw/Hjx+HXq9nm7AQhosNaJjsnzp1Kvbv34/FixcjODgY5eXl2LBhA377298iNjYWb7/9Nvbv34/S0lL2aMiuKRQKhIaGYvHixUhPT8f69esRHx8PT09PXLt2DatXr8aYMWMQHx+P1atXIz8/n0PJbYzhYkOEEOjYsSPeeustHDp0CKtXr0ZUVBRcXFyQnZ2Nv//975g0aRKioqIwffp0bN++HWVlZQwZslsNbeLFF1/E1q1bkZ6ejiVLlmDQoEEQQiAzMxPz5s1DTEwM5s2bh8uXL/OLVxthuNighgb1yiuvYNeuXTh48CCWLl2KmJgYeHt7o7CwEFu2bMHUqVMRExODhQsXmudniOyREALOzs7o2bMn3nrrLezbtw979uzBrFmz0LVrV5SUlGDNmjWIi4vD3LlzkZeXx4BpZQwXGyaEgKurK/r374833ngD33//PQ4fPox169Zh7Nix8PX1RV5eHpYvX464uDh88MEH+Omnn9ioyK4JIeDu7o6IiAh8/PHHSE9Px0cffYSnnnoKGo0Ga9euxahRo7By5UpotVq2h1bCcLETDbdsqlQqTJ06Fd988w3S0tKwcOFCBAcHo7S0FB9//DFGjRqF//mf/zHfAEBkzxQKBQIDA5GYmIgDBw7go48+Qo8ePVBUVIT33nsPv/3tb5GRkcHVyFsBw8UOCSGgVCoREhKCuXPn4uDBg3j//ffRqVMnXL9+HW+88QaeeeYZfP3115z8J4cghECHDh2QmJhovpXZ09MTJ0+exPPPP4/FixdDq9Vauky7wnCxc0IIdOrUCfPmzcP+/fvxyiuvwNPTE2fOnMGf//xnDB8+HC+++CLeeustbNu2zdLlErUqIQSCgoKwbNkybN++HZGRkdDpdPj444/xwgsv4MKFC/yiJROGi4MQQiAkJASrVq1CamoqXnzxRfj6+iI/Px87d+7EP/7xD6Snp1u6TKI24eTkhKFDh+K7777D3Llz4enpiYyMDLz33nvcT0kmSksXYG+s/VuPQqFAeHg4PvvsM9y4cQMZGRk4c+YMysvLER4ejitXrli6RLJx1t4GGmvXrh3ef/99DB06FB9++CHmz5+P77//3tJl2QWGiwyEELhw4QKWLl1q6VKeSIcOHeDv74+bN2/i3LlzEEJYuiSyMQ1tYMmSJZYu5YlFRkZi3759bAMyEZItfc2wUnV1deaF82ydQqGASqWCi4uLpUshG8I2QPdiuBARkew4oW8jJEniLcPk8NgObAfDxUacPXsWnp6eOHv2rKVLIbIYtgPbwXAhIiLZMVyIiEh2DBciIpIdw4WIiGTHcCEiItkxXIiISHYMFyIikh3DhYiIZMdwISIi2TFciIhIdgwXIiKSHcOFiIhkx3AhIiLZMVyIiEh2DBciIpIdw8UGSJIEjUbT5H+JHA3bgW1huFgxrVaLlStXQq1WY+TIkaitrcXIkSOhVquxcuVKaLVaS5dI1OrYDmyTkBj/Vmnv3r2YNGkS9Ho9ADT5liaEAAB4eHggOTkZCQkJFqmRqLWxHdguhosV2rt3L8aOHWveL/xBFAoFhBBISUlhwyK7w3Zg2xguVkar1SIoKAjV1dUPbVANFAoF3N3dUVBQAF9f39YvkKgNsB3YPs65WJkNGzZAr9c3q0EBgMlkgl6vx8aNG1u5MqK2w3Zg+9hzsSKSJEGtViM3N/ex7oQRQkClUuHq1avmcWgiW8V2YB8YLlakrKwMAQEBLTrf399fxoqI2h7bgX3gsJgV0el0LTq/srJSpkqILIftwD4wXKyIl5dXi8739vaWqRIiy2E7sA8MFyvi7++PkJCQxx4vFkIgJCQE7du3b6XKiNoO24F9YLhYESEEEhMTn+jcOXPmcBKT7ALbgX3ghL6V4f39RGwH9oA9Fyvj6+uL5ORkCCGgUDz8/56GJ5O3bdvGBkV2he3A9jFcrFBCQgJSUlLg7u4OIcQvuvkNP3N3d0dqairi4+MtVClR62E7sG0MFyuVkJCAgoICJCUlQaVSNTmmUqmQlJSEwsJCNiiya2wHtotzLjZAkiTcvn0blZWV8Pb2Rvv27TlpSQ6H7cC2MFyIiEh2HBYjIiLZMVyIiEh2DBciIpIdw4WIiGTHcCEiItkxXIiISHYMFyIikh3DhYiIZMdwISIi2TFciIhIdgwXIiKSHcOFiIhkx3AhIiLZMVyIiEh2/w+lxIsMjvXaAAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 500x400 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ccb7ec43",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "r2 is 0.9974619150161743\n",
      "saving model version 0.2\n",
      "r2 is 0.997527003288269\n",
      "saving model version 0.3\n",
      "r2 is 0.9740613698959351\n",
      "saving model version 0.4\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "tensor(0.9741, device='cuda:0')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.fix_symbolic(0,0,0,'log')\n",
    "model.fix_symbolic(0,1,0,'log')\n",
    "model.fix_symbolic(1,0,0,'sin')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "0937db67",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "| train_loss: 2.66e-07 | test_loss: 2.75e-07 | reg: 0.00e+00 | : 100%|█| 20/20 [00:01<00:00, 15.69it"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "saving model version 0.5\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "model.fit(dataset, opt=\"LBFGS\", steps=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e959cda3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - 1.0 \\sin{\\left(2.0 \\log{\\left(5.017 x_{1} \\right)} + 2.0 \\log{\\left(1.512 x_{2} \\right)} - 7.194 \\right)}$"
      ],
      "text/plain": [
       "-1.0*sin(2.0*log(5.017*x_1) + 2.0*log(1.512*x_2) - 7.194)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ex_round(model.symbolic_formula()[0][0], 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16e4da06",
   "metadata": {},
   "source": [
    "We were lucky -- singularity does not seem to be a problem in this case. But let's instead consider $f(x,y)=\\sqrt{x^2+y^2}$. $x=y=0$ is a singularity point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "1ce52cec",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "checkpoint directory created: ./model\n",
      "saving model version 0.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "| train_loss: 3.65e-03 | test_loss: 3.97e-03 | reg: 4.84e+00 | : 100%|█| 20/20 [00:03<00:00,  5.36it\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "saving model version 0.1\n"
     ]
    }
   ],
   "source": [
    "from kan import *\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,1,1], grid=5, k=3, seed=0)\n",
    "f = lambda x: torch.sqrt(x[:,[0]]**2+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": "code",
   "execution_count": 7,
   "id": "3a69ec41",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAFICAYAAACcDrP3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAsi0lEQVR4nO3de1TUdf4/8Od7ZlAYUQdwUAFFRllvaXlBlItgXnDla+tlTbdvmta2pqbVbu3Z3e+2Vmu1lZuXLEvtoq6uXajU8EheEFQ0FCUUL4ngZcALgsP9/vn8/kjmB+YF5TPzmRmej3M8neM4zAviNc95Xz7vj5BlWQYREZGCNGoXQERErofhQkREimO4EBGR4hguRESkOIYLEREpjuFCRESKY7gQEZHiGC5ERKQ4hgsRESmO4UJERIpjuBARkeIYLkREpDiGCxERKY7hQkREimO4EBGR4nRqF0DkDGRZRkFBAUpLS+Hp6QkfHx8IIdQui8hhceRCdAcWiwXLli1DcHAwjEYjgoKCYDQaERwcjGXLlsFisahdIpFDErwTJdGtJSQkYPLkySgvLwfw8+ilXv2oRa/XIy4uDjExMarUSOSoGC5Et5CQkIDY2FjIsgxJkm777zQaDYQQiI+PZ8AQNcBwIbqJxWJBQEAAKioq7hgs9TQaDTw8PGA2m2EwGGxfIJET4JoL0U3Wrl2L8vLyJgULAEiShPLycqxbt87GlRE5D45ciBqQZRnBwcHIzs7GvbSGEAImkwlnzpzhLjIiMFyIGrl27RqMRmOznu/j46NgRUTOidNiRA2UlpY26/klJSUKVULk3BguRA14eno26/lt27ZVqBIi58ZwIWrAx8cH3bt3v+d1EyEEunfvDm9vbxtVRuRcGC5EDQghMH/+/Pt67oIFC7iYT3QDF/SJbsLrXIiajyMXopsYDAbExcVBCAGN5s4tUn+F/tdff81gIWqA4UJ0CzExMYiPj4eHhweEEL+Y7qr/Ow8PD2zbtg1jxoxRqVIix8RwIbqNmJgYmM1mLF26FCaTqdFjJpMJS5cuRW5uLoOF6Ba45kLUBLIsIzExESNHjsSuXbswYsQILt4T3QFHLkRNIISwrqkYDAYGC9FdMFyIiEhxDBciIlIcw4WIiBTHcCEiIsUxXIiISHEMFyIiUhzDhYiIFMdwISIixTFciIhIcQwXIiJSHMOFiIgUx3AhIiLFMVyIiEhxDBciIlIcw4WIiBTHcCEiIsUxXIjuoqamBrm5uTh58iQA4OzZsygsLIQkSSpXRuS4eJtjotuwWCyIi4vDhg0bkJmZiZKSElRXV8Pd3R1GoxGRkZF46qmnEB4eDp1Op3a5RA6F4UJ0CwcOHMALL7yAjIwMhISEIDY2Fv3794enpycsFgvS0tKwdetWZGVlYerUqVi0aBGMRqPaZRM5DIYL0U2+//57zJw5E56ennjzzTcxbtw4VFdXY9OmTaiqqkK7du0wbdo01NTUYNOmTXjllVfQt29frF+/Hh07dlS7fCKHwHAhauCnn37C2LFj0aZNG2zatAl9+vSBEALZ2dkYOHAgioqKEBQUhLS0NHh5eUGWZezbtw+PPfYYoqOjsWbNGrRu3Vrtb4NIdVzQJ7qhrq4Ob7zxBq5fv44VK1ZYg+VOhBCIiIjA22+/jc2bN2P79u12qpbIsTFciG7IysrC1q1bMWnSJERERNw1WOoJITBhwgQMHToUq1evRm1trY0rJXJ83OJCdENKSgpKS0sxefJknDt3DmVlZdbHzGYz6urqAADV1dXIzMxEu3btrI/7+flh0qRJeOWVV3D58mUEBATYvX4iR8JwIbrh1KlT0Ov1MJlMmD17Nvbv3299TJZlVFVVAQDy8vIwevRo62NCCPz73/9Gv379UF5ejry8PIYLtXgMF6IbKioqoNPp0Lp1a1RVVaGysvKW/06W5V88VltbCw8Pj0YhRNSSMVyIbvD19UVFRQUsFgtCQ0PRpk0b62MVFRVISUmxhkhYWJj1wkkhBLp27YqrV69Co9HAy8tLrW+ByGEwXIhuGDRoEGpqapCamoq33nqr0WPZ2dkICQlBUVEROnbsiM8//xwGg8H6uBACf/vb39CpUydOiRGBu8WIrIYMGQKTyYS1a9eirKwMWq220Z96QghoNBrr32s0Gly6dAlfffUVYmNj0b59exW/CyLHwHAhusHHxwfPPvssjhw5guXLlzd5S3FVVRX++c9/oqKiArNnz27yFmYiV8ZpMaIGZs6cieTkZLz11lvQ6/WYM2cO3N3dAQA6nQ46nc46ipFlGSUlJXj99dexadMmLFmyBD179lSzfCKHweNfiG6Sn5+PefPm4bvvvkNMTAxeeOEF9O7dG6dPn4YkSWjVqhV69OiB1NRULF68GOnp6XjttdcwZ86cRtNnRC0Zw4XoFsrKyrB69WosX74cV65cgclkQnBwMNq2bYvr16/j9OnTyMvLw6BBg7Bw4UJERUVBo+EsM1E9hgvRHVy+fBm7du1CUlISfvzxR6SmpiIyMhLh4eEYM2YMQkNDodfr1S6TyOEwXIia6NChQxgyZAgOHTqEwYMHq10OkUPjOJ6oibRarXUbMhHdGbuEiIgUx3AhIiLFMVyIiEhxDBciIlIcw4WIiBTHcCEiIsUxXIiISHEMFyIiUhzDhYiIFMdwISIixTFciIhIcQwXIiJSHMOFiIgUx3AhIiLF8X4uRE0kyzIkSYJGo4EQQu1yiBwaRy5E94D3ciFqGp3aBRApoaamBhcuXIAkSWqX0mxCCHTt2hWtWrVSuxSi+8ZwIZdgNpsxd+5cDBo0SO1Smi0tLQ0ffPABunfvrnYpRPeN4UIuQZZl9O/fH3369MHu3bvh6+uLiRMnYsCAAXBzc1O7vHvyl7/8BVwKJWfHcCGXkpqairVr1wIA3n//fURERGD+/PkYPXo0dDqdwy/EM1TIVXB1klzKsGHD8MQTT2DQoEGoqanB999/jylTpuCpp55CVlYW37yJ7IThQi5l2rRp+Pjjj5GYmIhvvvkGY8aMQW1tLTZs2ICHH34Yy5cvR3FxMUOGyMYYLuRyhBBo06YNxo4di7i4OHz00UcICgpCXl4eXnzxRYwbNw67d+9GbW0tQ4bIRhgu5LKEENDr9XjiiSeQmJiIefPmQa/X48CBA5gwYQLmzZuHc+fOMWCIbIDhQi5PCIGAgAAsWbIE3333HaKiolBVVYU1a9ZgxIgRWLlyJUpKShgyRApiuFCLIISAVqtFREQEtmzZguXLl6Nr1664ePEinnvuOTzyyCPYt28f6urq1C6VyCUwXKhFEULA09MTs2fPRmJiIp5++mm4u7sjOTkZ48ePx0svvYRLly5xFEPUTAwXapGEEAgMDMSKFSuwefNmDBs2DGVlZVi2bBlGjRqFuLg4VFdXM2SI7hPDhVosIQR0Oh1GjBiB+Ph4vPHGG/D19cWpU6cwY8YMPP300zh//jwDhug+MFyoxRNCoH379njxxRexY8cOPPLII6irq8N//vMfjBo1Cp9//jlHMUT3iOFCdIMQAn379sXGjRvxwQcfICAgADk5OZg1axbmzJmDvLw8BgxREzFciBoQQsDDwwNPPvkkduzYgYkTJ0KSJHz22WcYO3Ysdu7cyR1lRE3AcCG6BSEEgoODsX79eixZsgRGoxEnTpzAlClT8NZbb6GsrIyjGKI7YLgQ3Ub9KGbOnDnYtm0bwsPDUVpaioULF+KJJ57AxYsXGTBEt8FwIboLIQQGDBiAb7/9Fs8++yzc3NzwzTff4H/+53+QkpLCgCG6BYYLURMIIeDt7Y133nkHH374IXx9fZGZmYnJkydjw4YNqK2tVbtEIofCcCG6B25ubpg+fTq++eYb9OvXD/n5+XjmmWfwxhtvoKKigqMYohsYLkT3SAiB0NBQbNmyBbGxsaisrMSiRYswf/58WCwWBgwRGC5E90UIgS5dumD9+vWYPXs2hBD47LPPMGPGDF4PQwSGC9F9E0KgXbt2ePfdd/Hqq6/C3d0d27Ztw9SpU3H27FkGDLVoDBeiZhBCoHXr1njxxRexYsUKtG/fHgcOHMCjjz6KkydPMmCoxWK4EClAp9NhxowZ+OSTT9ChQwf8+OOPmDZtGk6cOMGAoRaJ4UKkEI1Gg9/85jdYu3YtOnbsiMzMTPzv//4vsrKyGDDU4jBciBQkhEBMTAw+/fRTGI1GHDt2DDNnzuQiP7U4DBcihQkhMGbMGHz00UcwGAw4ePAg5s6di+LiYrVLI7IbhguRDQghMH78eCxevBju7u6Ij4/Hyy+/jOrqarVLI7ILhguRjWg0GkyfPh0vvfQSNBoNVq9ejU8//RSSJKldGpHNMVyIbEin0+HPf/4zJk+ejOrqarz88ss4cOAA11/I5TFciGzMw8MDixcvxoMPPoiCggI899xzuHr1KgOGXBrDhcjGhBDw8/PDsmXL4OXlhaNHj2LhwoU8SZlcGsOFyA6EEAgPD8df//pXaLVarFu3Dt9++y1HL+SyGC5EdqLRaPDMM8/g17/+NaqqqvD3v/8dZrOZAUMuieFCZEd6vR5vvvkmOnXqhKysLLzzzjvcPUYuieFCZEdCCPTu3Rt//OMfIYTAf/7zHxw7doyjF3I5DBciOxNCYObMmejTpw+KioqwcuVKhgu5HIYLkQq8vb2tNxnbsmULzp8/r3ZJRIpiuBCpQAiBiRMnws/PD1evXsV3333H0Qu5FIYLkUo6deqEsWPHAgC+/vprnjtGLoXhQqQSIQQmTJgAnU6H9PR05OTkqF0SkWIYLkQqEUIgJCQEfn5+KC4uxt69e9UuiUgxDBciFfn4+CAkJAQAsGvXLq67kMtguBCpSAiBhx9+GABw5MgR3lCMXAbDhUhF9VNj7u7uyM3NRXZ2ttolESmC4UKkMpPJhM6dO6OyshLp6elql0OkCIYLkcratWuHPn36AABSU1O57kIuQad2AURKc7Y3Z41Gg7Fjx8Ld3R3h4eHIyMhQuySiZmO4kEsQQuDYsWN49dVX1S7lvtTV1aFXr144e/Ysjh07BiGE2iURNYuQne1jHtEtVFdXIzs7G3V1dWqX0mwajQbdu3dHq1at1C6F6L4xXIiISHGcFiNqooafwzhtRXRn3C1G1ERHjx6FVqvF0aNH1S6FyOExXIiISHEMFyIiUhzDhYiIFMdwISIixTFciIhIcQwXIiJSHMOFiIgUx3AhIiLFMVyIiEhxDBciIlIcw4WIiBTHcCEiIsUxXIiISHEMFyIiUhzDhYiIFMdwIWoCWZZx/fp1AMD169fBG7gS3RnDhegOLBYLli1bhuDgYIwaNQqyLGPUqFEIDg7GsmXLYLFY1C6RyCEJmR/BiG4pISEBkydPRnl5OYBb3+ZYr9cjLi4OMTExqtRI5KgYLkS3kJCQgNjYWMiyDEmSbvvvNBoNhBCIj49nwBA1wHAhuonFYkFAQAAqKiruGCz1NBoNPDw8YDabYTAYbF8gkRPgmgvRTdauXYvy8vImBQsASJKE8vJyrFu3zsaVETkPjlyIGpBlGcHBwcjOzr6nHWFCCJhMJpw5c8a6HkPUkjFciBq4du0ajEZjs57v4+OjYEVEzonTYkQNlJaWNuv5JSUlClVC5NwYLkQNeHp6Nuv5bdu2VagSIufGcCFqwMfHB927d7/ndRMhBLp37w5vb28bVUbkXBguRA0IITB//vz7eu6CBQu4mE90Axf0iW7C61yImo8jF6KbGAwGxMXFQQgBjebOLVJ/hf7XX3/NYCFqgOFCdAsxMTGIj4+Hh4cHhBC/mO6q/zsPDw9s27YNY8aMUalSIsfEcCG6jZiYGJjNZixduhQmk6nRYyaTCUuXLkVubi6DhegWuOZC1ASyLCMxMREjR47Erl27MGLECC7eE90BRy5ETSCEsK6pGAwGBgvRXTBciIhIcQwXIiJSHMOFiIgUx3AhIiLFMVyIiEhxDBciIlIcw4WIiBTHcCEiIsUxXIiISHEMFyIiUhzDhYiIFMdwISIixTFciIhIcQwXIiJSHMOFiIgUx3AhIiLFMVyI7qKmpga5ubk4efIkAODs2bMoLCyEJEkqV0bkuHibY6LbsFgsiIuLw4YNG5CZmYmSkhJUV1fD3d0dRqMRkZGReOqppxAeHg6dTqd2uUQOheFCdAsHDhzACy+8gIyMDISEhCA2Nhb9+/eHp6cnLBYL0tLSsHXrVmRlZWHq1KlYtGgRjEaj2mUTOQyGC9FNvv/+e8ycOROenp548803MW7cOFRXV2PTpk2oqqpCu3btMG3aNNTU1GDTpk145ZVX0LdvX6xfvx4dO3ZUu3wih8BwIWrgp59+wtixY9GmTRts2rQJffr0gRAC2dnZGDhwIIqKihAUFIS0tDR4eXlBlmXs27cPjz32GKKjo7FmzRq0bt1a7W+DSHVc0Ce6oa6uDm+88QauX7+OFStWWIPlToQQiIiIwNtvv43Nmzdj+/btdqqWyLExXIhuyMrKwtatWzFp0iRERETcNVjqCSEwYcIEDB06FKtXr0Ztba2NKyVyfNziQnRDSkoKSktLMXnyZJw7dw5lZWXWx8xmM+rq6gAA1dXVyMzMRLt27ayP+/n5YdKkSXjllVdw+fJlBAQE2L1+IkfCcCG64dSpU9Dr9TCZTJg9ezb2799vfUyWZVRVVQEA8vLyMHr0aOtjQgj8+9//Rr9+/VBeXo68vDyGC7V4DBeiGyoqKqDT6dC6dWtUVVWhsrLylv9OluVfPFZbWwsPD49GIUTUkjFcqMU7d+4ckpKSsG/fPpSXl8NisSA0NBRt2rSx/puKigqkpKRYQyQsLMx64aQQAl27dsXVq1dRW1uLrKwshISEwN3dXa1viUh13IpMLc7FixeRnJyMPXv2ICkpCRcuXIAQAoGBgcjKysL777+P3//+942ek52djZCQEBQVFaFbt244fPgwDAaD9XEhBP72t79h8eLF0Gq1cHd3R2hoKIYPH47o6GiEhIRwizK1KAwXcnl5eXlISkqy/snJyQEA9O/fH8OHD0dUVBQiIiIgSRIiIiLg5eWF7du3N1qwv911LsDP02R5eXmIiorC+PHjMXPmTCQnJyM5ORl79+6FxWKBh4cHhg4dag2bQYMGoVWrVqr8PIjsgeFCLufKlSvYs2cPkpOTkZSUhKysLABA3759rWEyfPhweHt7/+K577//Pv70pz/h73//O/7yl79Yp77uFC6VlZV4/vnnsXXrVuzevRs9e/a0fr26ujpkZGRYa9m3bx+Ki4uh1+sxbNgwREVFITo6GgMGDICbm5sdfjpE9sFwIaeXn5+PpKQk6xv46dOnAQC9evVqFCZNOfurrKwMTz75JLZt24ZXX30Vc+bMgbu7O3JycjBkyBDrtFhqaioMBgNKSkrw+uuv46OPPsKSJUswa9asO3792tpapKenW+vdv38/SktL4enpibCwMERFRSEqKgoPPfQQD8Mkp8ZwIadTUFBgnXZKSkrCiRMnAADBwcHWMImKirrvc77y8/Mxb948fPfdd4iJicELL7yA3r174/Tp05AkCa1atUKPHj2QmpqKxYsXIz09Ha+99hrmzJkDrVZ7T69VU1ODo0ePWkdaKSkpKC8vR7t27RAeHm79Xvr373/PX5tITQwXcnjXr1/H3r17rZ/2jx07BgAwmUyNwsTPz0+x1ywrK8Pq1auxfPlyXLlyBSaTCcHBwWjbti2uX7+O06dPIy8vD4MGDcLChQsRFRUFjab5B15UV1cjLS3Nuj504MABVFZWwmAwICIiwvq9PvDAA4q8HpGtMFzI4RQVFWHfvn3WMPnxxx8hyzICAwOtb65RUVF2uVDx8uXL2LVrF5KSkpCdnY3Kykp4eXnhgQcewJgxYxAaGgq9Xm+z16+qqsKhQ4esYfPDDz+gqqoK3t7eiIyMtP4smnIOGpE9MVxIdSUlJdi/f791muvo0aOQJAn+/v7WBe+oqCgEBgaqWmddXR1kWYZGo1Ft1FBZWYkffvjBupU6NTUVNTU16NChQ6NRXM+ePRk2pCqGC9ldaWkpDhw4YP00fuTIEdTV1aFz586Ijo62vkkGBQXxDfIuysvLcfDgQWvYHD58GLW1tfD19W00yuvRowd/lmRXDBeyufo3wPowafgG2DBM+AbYfKWlpdaf9Z49exoFd8OwYXCTrTFcSHH1Uzf1YdJw6qZ+WzCnbuyjuLi40SixfsoxICDAekGnI0w5kuthuFCz1S8612+nbbjoPHz4cGuYcNFZfRaLpdH6VsPNEg3Dhqc6U3MxXOie1W+XrT+b6+DBg9btspGRkdYw4XZZx1dYWIj9+/dbRzb127yDgoKsU2jR0dHo3LmzypWSs2G40F01vNCv/tqL+gv96q+9GD58OC/0cwEFBQXWa4oaXqDao0ePRms293uBKrUcDBf6hYZHlCQlJVnv0Ojp6Ynw8HDryIRHlLi+/Px86xRaw6N1evbs2ShsOnTooHKl5GgYLmQ9XLH+DWT//v3WwxUbnnfFwxXp8uXLjcKm/lDQPn36WH9PbncoKLUsDJcWSJIkHD9+3PoGsW/fPlgsFri7u1tP6o2KiuKx8HRXt7udQb9+/ay/R/W3MaCWheHSAsiyjBMnTljfAPbu3YvCwkK0bt0aoaGh1jcB3tCKmut2N2J78MEHrb9n4eHhaN++vdqlko0xXFyQLMs4ffq0NUySk5Nx7do1uLm5YciQIdYLF0NDQ3krXrKpc+fONQqb3NxcaDQaDBgwwBo2YWFhaNu2rdqlksIYLi5AlmVkZWU1mp64evUqdDodBg8ebA2ToUOH2vSQRaI7kWUZOTk51t/RPXv24PLly9BqtRg4cKD1Gpthw4ahTZs2apdLzcRwcUI3N2lSUhIuXbr0iyYdOnQoPD091S6X6Jaa8qEoKirK5idPk20wXJzE+fPnrZ/2kpOTYTabfzG9MGzYsEb3fSdyJjdP5yYlJaGgoACtWrVCSEgIp3OdDMPFSfTv3x9nzpxptDAaFhYGg8GgdmlENiFJEk6ePPmLjSjr16/HlClT1C6P7oLh4iQkSYIQgmdzUYslyzJkWWYfOAmGCxERKY5ndyigpqYGFy9ehCRJapfSbEIIdOnShRdP0j1hD9DNGC4KyM3Nxfz58zFw4EC1S2m2I0eO4L333oPJZFK7FHIi7AG6GcNFAbIso1+/fnjttdcU/9r5+fnYu3cvDh06hMLCQnh5eWHw4MGIjIy0ycm0//d//wfOlNK9snUP7Nu3D4cOHUJBQQEMBoO1Bzp16qT467EHlMFwUZgSC42yLKO0tBSffPIJVq5c+YvpBiEEAgIC8Ic//AFPP/002rVrp9jrEjWXUr+LJSUl+PTTT2/bA/7+/tYeaN++PXvAwTBcHIwsyzh//jwWLFiAXbt2QZIkdOjQAQ899BA6d+6MK1euID09HRcvXsQ//vEP7Ny5EytWrED37t25g4ZcQv1FwgsWLEBiYiIkSYKPjw8GDBiAzp074+rVqzh69CjMZjMWLlyIXbt2sQccEMPFgciyjJ9++gkzZsxARkYGPDw8MGvWLMybNw9du3aFVquFJEkwm81YuXIl1qxZg6SkJDz66KNYt24d+vbty+YipybLMk6ePIkZM2YgMzMTer0es2bNwty5cxv1QG5uLlauXInVq1ezBxwU70HrIGRZhtlsxqxZs5CRkYEOHTpg1apVeOutt2AymaDT6SCEgFarRWBgIF5//XV88skn6NixI06ePImZM2ciJyeHw3pyWrIs49y5c5g5cyYyMzPh6+uL1atX41//+tcveqBr165YtGhRox6YNWsWzp8/zx5wEAwXB1FaWornnnsOR48ehbe3N1avXo3Jkyff9k6PWq0W48ePxyeffAKj0YjMzEzMnz8fRUVFdq6cSBlFRUWYP38+jh07BqPRiI8//hgTJ068aw98/PHHMBqNOH78OBYsWIDi4mI7V063wnBxAHV1dVi6dCkSEhLg7u6Od955B2PGjLnr8F4IgREjRmDJkiXQ6/VITEzE22+/jdraWjtVTqSM2tpavP3229i9ezf0ej3effddjBw5skk98PDDD2Px4sXQ6/XYuXMn3n33XdTV1dmpcrodhovKZFnG/v378d5770GWZTzzzDN49NFHodE07X+NEAITJkzAggULIITAqlWrkJiYyKkBchqyLGP37t1YtWoVhBBYsGABJk6c2OS1EyEEJk2ahLlz5wIAPvjgAyQnJ7MHVMZwUVlRURFefvlllJSUICQkBC+99NJtpwFuR6vV4vnnn0dYWBjKysqwcOFCFBYW2qhiImUVFhbiH//4B8rKyhAeHo7nn38eWq32nr6GTqfDH//4R4SGhqK0tBQLFy7kFLHKGC4qkmUZa9euxeHDh+Hp6YlFixbd973G27Vrh3/+859o37490tPT8fHHH/OTGzk8SZKwatUqZGRkwGAwYNGiRfd924j657dt2xZpaWn47LPP2AMqYrio6MKFC3jvvfcgSRIef/xxhIWF3fc2SiEEQkJCMHPmTMiyjA8//BDZ2dkKV0ykHFmWcfbsWXz00UeQZRmzZs3C4MGDm9UDQ4cOxfTp0yFJElasWIELFy4oXDU1FcNFJZIk4cMPP0Rubi78/f3x3HPP3fNUwM20Wi2effZZBAYG4tKlS3j//fe5sEkOqz4Arly5gqCgIMybN6/Ja423o9VqsWDBAvj7+1uvhXGFwzSdEcNFBfW3d92wYQOEEPjDH/6AwMBARb62v78/5syZAyEEPv/8c5w+fZpTA+RwZFnGiRMn8MUXX0AIgblz58LPz0+Rr921a1fMnj0bQghs3LgRZ8+eVeTr0r1huKhAlmWsWbMG+fn56NatG6ZPn67YVcVCCDz22GMIDg5GYWEhVq1axXAhh1M/crdYLOjZsyd+97vfKdoD06dPR7du3ZCfn4/Vq1dz9KIChosKLly4YP3E9uSTTyp+smuHDh3w+9//HkIIxMXFce2FHM6ZM2fw7bffWkfu3t7ein79jh074sknn4QQAl9++SXXXlTAcLEzWZaxceNGXLlyBX5+fpg2bZriZyEJITBlyhQEBgbi2rVrWL9+PUcv5DBkWca6detQWFiIoKAgTJ482SY98Lvf/Q7+/v64fPkyNm7cyB6wM4aLnRUUFGDjxo0AgKlTp8Lf398mr9OxY0c89thjAIDPP/8cV65cscnrEN2rS5cu4csvvwQAPP744zAajTZ5HT8/P0ydOhUAsHHjRhQUFNjkdejWGC52JMsytm/fjpycHHh5eeHxxx+32Qmu9WsvHTp0wIULFxAfH89PbqQ6WZaxZcsW5ObmwtfX1yYj93pCCDz++OPw8vJCTk4Otm/fzh6wI4aLHVVVVWHdunWQJAkjR45EcHCwTV+vW7duiImJgSzLWL9+PSoqKmz6ekR3U15ejg0bNkCWZYwbN06xXZK306NHD4waNQqSJGHdunWoqqqy6evR/8dwsRNZlpGeno7Dhw/Dzc0NTzzxRLOva7kbjUaDGTNmoHXr1khPT0daWho/uZFqZFlGamoqMjIy4O7uruguydvRarWYMWMG3NzccPjwYaSnp7MH7IThYieyLOOLL75ARUUF+vbti6FDh9q8sYQQGDx4MB588EFUVVXhv//9LxuLVCPLMv773/+iuroaAwYMwIABA+zSA0OHDsUDDzyAiooKfPHFF+wBO2G42Mm1a9cQHx8PAJgyZQratGljl9f18PDAo48+CiEEEhIScPnyZbu8LtHN8vLysHPnTgghMHXqVLi7u9vlddu0aYMpU6YAAOLj45Gfn2+X123pGC52IMsy9uzZA7PZDC8vL4wfP95ut2IVQiA2NhYdOnTApUuXsHPnTn5yI7uTZRk7duzA5cuX4evri3Hjxtm1B8aPHw9vb2+YzWbs2bOHPWAHDBc7qKurw1dffQVJkhAREYGgoCC7vn6XLl0wfPhwyLKMr776ijcTI7urqanBV199BVmWER0drdhRL03VrVs3REREQJIkfPnllzxzzw4YLnZw4cIF7N+/HxqNBr/97W9tvpB/M41GgylTpkCr1SI1NZVnLZHdZWVl4fDhw9Bqtfjtb3/b7AMq75VWq8WUKVOg0WiQkpKC8+fP2/X1WyKGi43VTwcUFhaic+fOGD58uN2mA+oJIRAeHo4uXbqgqKiI+/3Jruqv7youLkZgYCCGDRumSg9ERkbCz88P169fR0JCAnvAxhguNlZTU4PNmzcDAKKjo+Hr66tKHd7e3hg1ahQAYMuWLdzvT3ZTVVWFLVu2AABGjx6t+DliTWU0GjFixAgAwObNm1FTU6NKHS0Fw8XGzp07hyNHjkCr1WLChAl2/8RWTwiB3/zmN9DpdDh27BjOnDmjSh3U8vz00084fvw43Nzc8Mgjj6hWhxACEyZMgFarRXp6OqeHbYzhYkOyLGPnzp0oKiqCv78/QkNDVQ2XQYMGoVu3bigtLcX333/PaQGyOVmWkZCQgLKyMnTr1g0DBw5UtQeGDBmCLl26oLi4GDt27GAP2BDDxYZqa2ut17ZERUXBx8dH1XoMBgMefvhhAMC2bdtQXV2taj3k+qqqqrBt2zYAwKhRo9C+fXtV6/H29kZ0dDSAn3uAU2O2w3CxoYsXL+Lo0aPQaDSIjY1V7RNbQ7GxsdapsZycHLXLIReXnZ2NzMxM6HQ6jBs3Tu1yrNd9abVa/Pjjj9w1ZkMMFxuRZRnJycmwWCzo1KmTqlNi9YQQGDhwIPz9/VFSUoLExEROC5DNyLKM3bt3o6SkBAEBAXY57uVu6qfGOnfuDIvFgqSkJPaAjTBcbESSJOuW36FDh9rsnhX3ytvbGxEREQCA7du382Iyspna2lokJCQAACIjI+Hl5aVyRT/r0KEDhg0bBoA9YEsMFxu5evUqUlNTIYTAr3/9a7tfNHY7Des5cuQILl26pHZJ5KLy8vKs08Jjx45VfdRSTwhhPX7m8OHDvJGejTjGO56LkWUZaWlpuHr1KgwGA8LCwhyqsUJDQ+Hj44PCwkL88MMPnBYgxcmyjIMHD6KwsBBGo9EhpoXrCSEwbNgweHt7Iz8/H4cOHWIP2ADDxUZ27NiBuro69OvXDwEBAWqX00inTp0wYMAASJLELclkE7IsW3+3Bg4ciI4dO6pdUiN+fn548MEHIUkSduzYoXY5LonhYgOlpaXYu3cvgJ+3X7q5ualcUWNardZ6tX5KSgqKi4tVrohcjcViwYEDBwD8fFW+vc/TuxudTmftgX379qGkpETlilwPw8UGTp8+jZycHLRu3RpRUVEOMx1QTwiB4cOHQ6/X4+LFizh58qTaJZGLyczMhNlsRps2bRAZGemwPeDu7o7z58/j1KlTapfkchguCpNlGXv37kVlZSUCAwPRq1cvtUu6pR49esBkMqG6uprbMUlRsiwjMTERNTU1CA4OhslkUrukW/rVr36FoKAgVFVVITk5mT2gMIaLwurq6pCYmAgACAsLQ9u2bVWu6Nb0er11S/KePXt4jxdSTHV1Nfbs2QPg5y3IHh4e6hZ0G56enggLCwMAJCYmsgcUxnBR2NWrV5GRkQEhhPWoFUcVHR0NjUaDzMxM3v6YFJObm4tTp05Bq9UiOjra4abEGhoxYgQ0Gg2OHTvGLckKY7goLD09HdeuXYPBYMDgwYMdtrHqD7L09vZGYWEh0tLS1C6JXMThw4dhsVhgNBrx0EMPqV3ObQkhMHjwYHh5eaGgoADp6elql+RSGC4KS0pKQl1dHXr16mX3W7neq44dO6Jv376QJAnJyclql0MuoH69RZZl9OvXT7X7FzVV586d0bt3b0iShKSkJLXLcSkMFwVVVVUhJSUFABAREYFWrVqpXNGd6XQ6REZGAgAOHDiAiooKlSsiZ1dWVoYffvgBwM8ngTvaFuSbubm5WXsgJSUFlZWVKlfkOhguCjKbzThz5gx0Op0qtzO+V0IIREREQKfT4ezZs7h48aLaJZGTy8nJwfnz59GqVSuEh4c7RQ9ERkZCp9MhKysLZrNZ7ZJcBsNFQUeOHEFxcTGMRiP69u2rdjlN0rt3b/j6+qK0tJTrLtRsBw8eRHl5Ofz9/fGrX/1K7XKapE+fPvD19UVJSQl7QEEMFwXV75V/4IEHHOYU5Lvx9vZG//79rbcI4F5/ul+yLFunhQcOHAiDwaBuQU3k4+ODfv36sQcUplO7AFchyzIiIyNRUFCAESNGQKvVOsUvqVarxSOPPIL27dtj9OjROHTokNolkZMSQuBPf/oTBg8ejJ49e0II4TQ9EBsbC3d3dwwfPpy7xhTCcFGAEAKZmZnw9PREr169kJeXh0WLFqldVpNJkgSTyYQTJ07g+PHjDj9PTo5HCIHjx4/j22+/BfDz9NjBgwfVLeoe1NXVoXfv3jhz5gx7QCFCdoaPFg6uuroaOTk5LnHTIY1GA5PJ5PA73cixsAfoZgwXIiJSHBf0nYQsy5AkySnmsIlshX3gPBguTiI9PR16vZ6LjdSisQ+cB8OFiIgUx3AhIiLFMVyIiEhxDBciIlIcw4WIiBTHcCEiIsUxXIiISHEMFyIiUhzDhYiIFMdwISIixTFciIhIcQwXIiJSHMOFiIgUx3AhIiLFMVyIiEhxDBcnIMsyrl+/3ui/RC0N+8C5MFwcmMViwbJlyxAcHIyRI0eiuroaI0eORHBwMJYtWwaLxaJ2iUQ2xz5wTkJm/DukhIQETJ48GeXl5QDQ6FOaEAIAoNfrERcXh5iYGFVqJLI19oHzYrg4oISEBMTGxlrvF347Go0GQgjEx8ezscjlsA+cG8PFwVgsFgQEBKCiouKODVVPo9HAw8MDZrMZBoPB9gUS2QH7wPlxzcXBrF27FuXl5U1qKACQJAnl5eVYt26djSsjsh/2gfPjyMWByLKM4OBgZGdn39NOGCEETCYTzpw5Y52HJnJW7APXwHBxINeuXYPRaGzW8318fBSsiMj+2AeugdNiDqS0tLRZzy8pKVGoEiL1sA9cA8PFgXh6ejbr+W3btlWoEiL1sA9cA8PFgfj4+KB79+73PF8shED37t3h7e1to8qI7Id94BoYLg5ECIH58+ff13MXLFjARUxyCewD18AFfQfD/f1E7ANXwJGLgzEYDIiLi4MQAhrNnf/31F+Z/PXXX7OhyKWwD5wfw8UBxcTEID4+Hh4eHhBC/GKYX/93Hh4e2LZtG8aMGaNSpUS2wz5wbgwXBxUTEwOz2YylS5fCZDI1esxkMmHp0qXIzc1lQ5FLYx84L665OAFZllFYWIiSkhK0bdsW3t7eXLSkFod94FwYLkREpDhOixERkeIYLkREpDiGCxERKY7hQkREimO4EBGR4hguRESkOIYLEREpjuFCRESKY7gQEZHiGC5ERKQ4hgsRESmO4UJERIpjuBARkeIYLkREpLj/BxNiOy3uy0HpAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 500x400 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "abef7aa9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "r2 is 0.9999973773956299\n",
      "saving model version 0.2\n",
      "r2 is 0.9999948740005493\n",
      "saving model version 0.3\n",
      "r2 is 0.9998846650123596\n",
      "saving model version 0.4\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "tensor(0.9999)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.fix_symbolic(0,0,0,'x^2')\n",
    "model.fix_symbolic(0,1,0,'x^2')\n",
    "model.fix_symbolic(1,0,0,'sqrt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "aa71848c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rewind to model version 0.4, renamed as 1.4\n"
     ]
    }
   ],
   "source": [
    "model = model.rewind('0.4')\n",
    "model.get_act(dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e14000d8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 1.00775534257195 \\sqrt{0.999962771771901 \\left(6.10769914067904 \\cdot 10^{-5} - x_{1}\\right)^{2} + \\left(9.20887777110479 \\cdot 10^{-5} - x_{2}\\right)^{2} + 0.00441348508007971} - 0.00955450534820557$"
      ],
      "text/plain": [
       "1.00775534257195*sqrt(0.999962771771901*(6.10769914067904e-5 - x_1)**2 + (9.20887777110479e-5 - x_2)**2 + 0.00441348508007971) - 0.00955450534820557"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula = model.symbolic_formula()[0][0]\n",
    "formula"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "c56ee3d5",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 1.01 \\sqrt{1.0 x_{1}^{2} + x_{2}^{2}} - 0.01$"
      ],
      "text/plain": [
       "1.01*sqrt(1.0*x_1**2 + x_2**2) - 0.e-2"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ex_round(formula, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fd57d41",
   "metadata": {},
   "source": [
    "w/ singularity avoiding (LBFGS may still get nan because of line search, but Adam won't get nan)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "de708f21",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "| train_loss: 5.11e-04 | test_loss: 5.64e-04 | reg: 0.00e+00 | : 100%|█| 1000/1000 [00:14<00:00, 70.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "saving model version 1.5\n"
     ]
    }
   ],
   "source": [
    "model.fit(dataset, opt=\"Adam\", steps=1000, lr=1e-3, update_grid=False, singularity_avoiding=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fd34c4c",
   "metadata": {},
   "source": [
    "w/o singularity avoiding, nan may appear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "031fabd6",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "| train_loss: nan | test_loss: nan | reg: nan | : 100%|█████████| 1000/1000 [00:17<00:00, 57.55it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "saving model version 1.6\n"
     ]
    }
   ],
   "source": [
    "model.fit(dataset, opt=\"Adam\", steps=1000, lr=1e-3, update_grid=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "124c9ca4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}