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": "stderr",
     "output_type": "stream",
     "text": [
      "train loss: 5.00e-03 | test loss: 3.94e-02 | reg: 2.72e+00 : 100%|██| 20/20 [00:11<00:00,  1.78it/s]\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,1,1], grid=20, k=3, seed=0)\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])\n",
    "\n",
    "# train the model\n",
    "model.train(dataset, opt=\"LBFGS\", steps=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3f95fcdd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAFICAYAAACcDrP3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAvQklEQVR4nO3d6XNUV34+8Od2q1vqTUtrXxBakJCQABvM6g3PEMPYhInHs9jJULiSmVQqiV12/oG8ynuIp1yZSX7jMVMpMuPALA4eMF4BLxixGLSCFrQLqSW11OpF6u57fi+Ye6uvNiR0pd6eT5Wr5l51w7FHt58+53sWSQghQEREpCNDtBtARESJh+FCRES6Y7gQEZHuGC5ERKQ7hgsREemO4UJERLpjuBARke4YLkREpDuGCxER6Y7hQkREumO4EBGR7hguRESkO4YLERHpjuFCRES6Y7gQEZHuUqLdAKJ4IITA6OgopqamYLfbkZ2dDUmSot0sopjFngvRItxuN44fP46qqirk5uaivLwcubm5qKqqwvHjx+F2u6PdRKKYJPEkSqL5nTt3Di+++CJ8Ph+A+70XhdJrsVqtOHXqFA4cOBCVNhLFKoYL0TzOnTuH559/HkIIyLK84OsMBgMkScKZM2cYMEQRGC5Es7jdbpSUlMDv9y8aLAqDwQCLxYK+vj5kZmaufgOJ4gBrLkSzvPPOO/D5fEsKFgCQZRk+nw8nTpxY5ZYRxQ/2XIgiCCFQVVWFzs5OLOfRkCQJFRUVuHPnDmeREYHhQqThcrmQm5u7ovdnZ2fr2CKi+MRhMaIIU1NTK3q/x+PRqSVE8Y3hQhTBbrev6P0Oh0OnlhDFN4YLUYTs7GxUVlYuu24iSRIqKyvhdDpXqWVE8YXhQhRBkiS8+uqrD/Xe1157jcV8oj9jQZ9oFq5zIVo59lyIZsnMzMSpU6cgSRIMhsUfEWWF/unTpxksRBEYLkTzOHDgAM6cOQOLxQJJkuYMdyn3LBYL3n//fTz77LNRailRbGK4EC3gwIED6Ovrw7Fjx1BRUaH5WUVFBY4dO4b+/n4GC9E8WHMhWgIhBNra2nDy5Em8/PLL2LhxI4v3RItgz4VoCSRJgtPpRHZ2NpxOJ4OF6AEYLkREpDuGCxER6Y7hQkREumO4EBGR7hguRESkO4YLERHpjuFCRES6Y7gQEZHuGC5ERKQ7hgsREemO4UJERLpjuBARke4YLkREpDuGC9ESWSwWbNmyBRaLJdpNIYp5PM+FaInC4TB8Ph+sViuMRmO0m0MU0xguRESkOw6LERGR7lKi3QAihSzLGBoawszMTLSbEvfy8vJgtVqj3QxKYgwXihmyLOPChQsIhUIwm83Rbk5cEkJgYmIC+/fvR0VFRbSbQ0mM4UIxQwgBSZLw5JNPoqSkZN7XhMNhCCGQksJf3fmEw2H88Y9/jHYziBguFHsMBsOc2VhCCAwPD+PixYsIh8PYs2cP1q1bB0mSotTK2MX/JhQLWNCnuCCEwLVr1zA4OIjh4WFcvnwZ4XA42s0iogUwXCguhMNhuN1u9XpiYgJ+vz96DSKiRTFcKC4Eg0FNmExPT2NiYiKKLSKixTBcKC4EAgHNFGVZluFyucA1wESxieFCccHr9SIUCmnujYyMRKk1RPQgDBeKC1NTU3N6KW63G7IsR6lFRLQYhgvFPCEEPB7PnPs+nw/BYDAKLSKiB2G4UFyYnJycc29mZoZbxRDFKIYLxTxZljE1NTXnfjAYRCAQiEKLiOhBGC4U88LhMLxe75z7sizPe5+Ioo/hQjFvZmZm3h6KEAJer5fTkYliEPcWo5g3e41LpPkK/WtJCIGpqSl4PB6kpaUhIyMDBoOB+3tR0mO4UMzzer0L7iMWzXAJhUK4fv06bt68Cb/fD5PJhPLycuzduxc2m40BQ0mNw2IU05RpyAsNfXm93qisdZFlGQ0NDbh8+TJ8Ph+EEJiZmUFbWxs++ugjTE9Pr3mbiGIJw4ViXmTvxGAwaA4S8/v9ax4uQgj09/fjxo0b84Zeb28vvvnmG9aCKKkxXCjmRYaL2WyG0+lUrwOBwJovpAyFQmhoaFD/XkmSUFJSApPJBOB++Ny6dQvj4+Nr2i6iWMJwoZgmy7ImXCwWiyZcFppJtlqEEBgYGMDg4KB6Lz8/HwcPHsTmzZvVe36/H42Njey9UNJiuFBMC4fD8Pl86rXVatWESygUWtO1LkIINDc3qxMMDAYDtm3bhrS0NGzZsgUOh0N9bUdHx7yLP4mSAcOFYtr09LSmOO5wOJCRkaHOxFpo37HVMjk5ib6+PvU6JydHPW7Zbrejurpa/dnU1BS6u7vZe6GkxHChmOb3+zU1FYfDAbvdDqPRqN6bnJxckw9wIQS6u7s1w3DV1dVqrQUAqqqqNBMO2tvbuXMzJSWGC8W02Wtc0tPTYbFYNB/o821quRpkWUZnZ6d6nZqairKyMrUXJUkSnE4nCgoK1NcMDw/zxExKSgwXilmz17goQ0+pqalIS0tTXzc1NbUmvYPJyUnNAWX5+flIT0/XvMZgMKC8vFy9np6eRl9fH4fGKOkwXCimRdZTUlJSYLPZYDQaYbVa1fuLreDXixACg4ODmvpPWVkZDAbtI6RMS05NTVXv9fT0MFwo6TBcKKZFDnmZTCakpaXBYDBoZmUFAoFVXxEvhEBPT4+mLcXFxfNu8ZKeno7s7Gz1enh4WDPjjSgZMFwoZs0+x8VisajF8sjhqGAwCL/fv6ptmZ6exr1799Rrp9OJjIyMeV9rNBpRUlKiXvt8PoyMjLD3QkmF4UIxa3ZoKENikiRpwmWhw8T0NDY2pvk7ioqKkJIy/76vytCYMmSmbBdDlEwYLhSzpqenNdN+HQ6HOgxlt9vnrHVZrZ6BsipfmTQgSRKKi4sXfY/T6YTNZlOvBwcHV70uRBRLGC4Us/x+P0KhkHodWWex2WyansNqTkeWZRkDAwPqtcViQU5OzqJb6qempiIvL0+9Hh8f56mZlFQYLhSzZm+nn56ern6gp6WlaWZkreZCykAggLGxMfXa6XRqZqvNR5IkFBUVqdczMzNwuVysu1DSYLhQTBJCaALDYDDAbrerPzeZTLBYLOr1aq51GR8f18z2KigomDMFeTZJkpCfn6/uJKBMZSZKFgwXilmz17hE9haMRqMmbHw+36psvS+EwL179zT1lsLCwiW9NzMzU1N3ifxziBIdw4ViktJzUZjNZs2qfEmSNFOBA4HAqtQ0Zvc40tLSkJ2dvaQjjFNTUzXrXWb3gIgSGcOFYpIsy5qwmL2fGADNB3c4HIbb7da9pjE9Pa2pt2RkZGiG4xYjSZJmn7FAIMADxChpMFwoJgWDQc23fLvdrqlzSJKErKwsTU1jdHRU93ZMTk5qQi6yjvIgkiQhLy9PM2X63r17LOpTUmC4UEwKBAKYmZlRryPXuETei9zefnR0VNcPbiEEXC6XZjp0ZE9kKZxOp2Y4j+FCyYLhQjHJ5/MtuMZFkZaWprk/Pj6ueY8ehoaG1P9tMpmWXG9RWCwWZGZmqtdjY2Oa0CRKVAwXikmzpxZHrnFRGI1GTd1lampK16J+OByGy+VSrx0Ox7whtxiDwaBZTOn1etfs/BmiaGK4UEyK/ACePe04UuQHdzAY1BTfV8rn82na4XQ6F9xPbDGRbQyFQlxMSUmB4UIxKS8vD/X19Vi/fj3y8/M160UUkiQhJydHs0GknrsPj4+Pa7byz8/PX9aQWGQbI2e6Re6uvFI+nw8ulwvT09NcQ0MxZflfw4jWQGlpKUpLSyGEgCzLC87QSk9PR1pamjqzTAmX5YbAbEIIDA8Pa3YIyM3Nfag/1+FwwGazwe12AwBcLhfC4fBD9YJmt/H27dv48ssvYbPZkJ2djd27d6/ozyTSC3suFJMkSYIkSTAYDEhJSVnwQ91isWgWU46Pj+tSMFfCRZGamqopzC+HyWSC0+lUrycmJnQ7f+bevXsIhUKYmJhAf3//ikOVSC8MF4prBoMBOTk56rXX69XlbJfZ9ZvlLJ6cTdlnTKHXYspgMKhZ25Oenj7v8CFRNDBcKO7NLpiPjY2tuO7i8Xg0M89ycnKWvHhyoTZGLqbUozbk9Xo1+69lZ2eveKiNSC8MF4prkiTN+VCNHM56GMpq/8iNMCN7Hg/TxszMTM0RAZH1nIc1NjamWxuJ9MZwobinFPUVLpdrxTOnImd0paSkPPBwsAexWq2a2tDs8Fqu2VvJGI1G5ObmPvSfR6Q3hgvFPbPZjKysLPV6pUV9WZYxMjKiXlut1mUvnpzNaDRqakNTU1Mrqg3NbuPsiQ1E0cZwobinTBNW+P3+Fa2CDwQCmJiYUK+zsrI0Q1oPK3LYSinGP+zQ2MzMjGbCwexhN6JoY7hQQogMF2Xblof94Ha73ZqpwpHF+IelLKaMrA2tZDHl7OnMeXl5Dzwdk2gt8beR4p5S1NdjFbwykyvy5MnI2WgrkZ6erpnOHPn3LLeNykJMvdtIpBeGCyUEu92u2X9sZGRE/fBdrsjZZmazGU6nU5fFiampqZrakNvtRiAQeKg/a6W7NROtNoYLJQTlA1YxOTn5UEcKB4PBOTsh67UwcfZiSr/fr6ntLJWy+aVidrASxQKGCyWE2R/cs48nXqqpqSnNwsTZdZKViqzfKDO+llsbmr1t/+yNMYliAcOFEoJSd4jcIflhFirOXpi43JMnH9RGp9OpOT1zubUhIcScA8e4eJJiEcOFEkZmZuaKjhQWQmBoaGjOwkQ9axk2mw3p6enq9cjIyLIXU0a2UTmMjPUWijUMF0oY8x0pHHkey4PIsqzpScxeVa8HZbW/YrmLKWVZ1kw4sFqtD71bM9FqYrhQwph9pPDskyQfxO/3a3Yrzs7OXpWFiZFDbcoEgqX2sAKBgKaWpNcCTyK9MVwooUTWH0Kh0JIL5spmlZFTgwsKCnQfblJqQ8okASEEBgcHl/z+8fFxzeLJgoICLp6kmMTfSkoYyir4yIL5wMDAkt8/u5bxMMcaL8Xsc1fu3bu3pDU5Sk0ocoHnagQgkR4YLpRQ7Ha7pmA+PDy8pE0sZVnW9CAsFovm9Eg9mc1mTd3F7XZrzo5ZyOxeTlpammZtD1EsYbhQQjGZTJqahsfjWdKpj36/X3Oqo9PpfOiTJx9EkiQUFRWp1zMzM0savgsEAprFk1lZWbBaravSRqKVYrhQwikpKVGHikKhEAYHBxf94FbqLZG1jKKiolWrZSjDWcrJlkKIBw7fCSEwPj6u2XWA9RaKZfzNpISiFMwjZ1D19fU9sFfQ39+vqbcUFhauai0jMzNTs2XL4OAgQqHQou8ZGBjQ1FuKi4tXrX1EK8VwoYRjt9s1tYjh4eFF9xkLh8Po7+9Xr20226rXMsxms2batNvt1mw7M5ssy5o2WiyWFZ+OSbSaGC6UcAwGA0pKStRrv9+vmQk22+TkpKYuk5eXp1npvxokSdK0cWZmZtHhO6/Xq6m35OTkrFpNiEgPDBdKOJIkobS0VLOWpKura97XKvWOyJX869atW/UegSRJKCws1Eyb7unpmTdclCnIkWtwSkpKWG+hmMbfTkpITqdTsy3KwMCApmCvEEKgu7tbvTabzSgqKlqT4ab09HTN+S5DQ0PzthEAuru71eBJSUlBcXExh8QopjFcKCGZTCaUlpaq1x6PZ95hp6mpKc3BW9nZ2brvJ7aQlJQUrFu3Tr32er3zDt/5/X5NvSUzM3PV1uAQ6YXhQgmrrKxMM923vb1d88EthEBvb6+m2L9+/Xr1PatNkiTN3yeEQEdHx5w2Dg4Oaja3LC0t5fktFPMYLpSQJElCbm6u5ht+X1+f5kNalmV0dHSo12azGeXl5Ws63JSTk6NpY29vr6aNQgjcuXNHcwxAWVnZmrWP6GExXChhmUwmVFZWqtc+nw+dnZ0QQqgLJyO3U8nLy9PUQNaqjRUVFZo2dnV1qW10u93o7e1Vf56dna37GTNEq4HhQglLkiRUVlZqFlS2tLRgenoaQgg0Njaq+45JkoSqqqo1n4E1Xxubm5vV2WstLS2aWWJVVVUcEqO4wHChhJaZmYn169er16Ojo2hvb8fIyAja29vV++np6aioqIhKjyArK2tOG5ubmzEyMoKWlhb1vs1mw4YNG9hrobiQEu0GEK0mSZJQX1+Pzs5OhEIhCCHw9ddfIy0tTbO2ZdOmTVFblChJEjZv3oyuri4Eg0EIIdDQ0IDGxkbN1OTa2lo4HI6otJFoudhzoYSmbBIZWXvxer1zdkCura2NWo9AaePGjRvVe9PT05iYmFCvs7KysHnzZvZaKG4wXCjhGQwG7Ny5c96z5s1mM/bu3Rv1reslScLOnTtRWFg452dpaWl48sknNQeMEcU6hgslPEmSkJGRgYMHDyIvLw+SJEGSJDgcDjzzzDNYv3591HsEkiTBarXiwIED2LhxI9LS0tSzaQ4cOIDS0tKot5FoOVhzoaSgHIH83e9+F8PDw5BlGTk5ObDZbDHzoS1JEux2O/bv3w+v14twOAybzYaUlJSYaSPRUjFcKKYomzQGg8FV/7uGh4dX/e9YqbGxsWW9XgihmbpMFC0MF4oZkiTB6XTi7t27uHv3brSbE7fMZrNm3QxRNEjiQUf0Ea0RIYR60iKtjMFg4FAaRRXDhYiIdMfZYkREpDuGCxER6Y7hQkREumO4EBGR7hguREsUDofh8XgQDoej3RSimMdwIVqi0dFRvP3225pNL4lofgwXIiLSHcOFiIh0x3AhIiLdMVyIiEh3DBciItIdw4WIiHTHcCEiIt0xXIiISHcMFyIi0h3DhYiIdMdwISIi3TFciIhIdwwXIiLSHcOFiIh0x3AhWgIhBEZHR9V/hBDRbhJRTJMEnxKiBbndbrzzzjt488030dHRod6vrKzEq6++iqNHjyIzMzN6DSSKUQwXogWcO3cOL774Inw+HwBoeiuSJAEArFYrTp06hQMHDkSljUSxiuFCNI9z587h+eefhxACsiwv+DqDwQBJknDmzBkGDFEEhgvRLG63GyUlJfD7/YsGi8JgMMBisaCvr49DZER/xoI+0SzvvPMOfD7fkoIFAGRZhs/nw4kTJ1a5ZUTxgz0XoghCCFRVVaGzs3NZM8IkSUJFRQXu3Lmj1mOIkhnDhSiCy+VCbm7uit6fnZ2tY4uI4hOHxYgiTExMrOj9Ho9Hp5YQxbeUaDeAKJpmZmbQ3t6O5uZmNDc34+bNmyv6806ePIldu3ahrq4OeXl5HCKjpMVwoaTi8/nQ1taGlpYWNDc3o729HaFQCDabDTU1Nfjxj3+Mr7/+Gr29vcv+s3NyctDX14fPP/9cva6vr0ddXR3q6uqwbt06hg0lDdZcKKFNTk6qQdLS0oKuri4IIZCZmYlNmzahtrYWmzZtQmlpqfrBf/z4cbzxxhvLLugfO3YMr732GjweD1paWtDY2IjGxka0t7dDlmWkp6erQVNfX4/y8nIYjcbV+lcniiqGCyUUl8ulhklzczP6+/sBAHl5edi0aZMaKAUFBQv2IvRe5xIIBNDa2orGxkY0NTWhra0NMzMzsFgsqK2tRX19Perr61FVVQWTybSif3+iWMFwobglhMDg4KAaJC0tLRgZGQEAlJSUaMJkuTO4lrtC//3338ezzz67pD87GAzizp07aGpqUv/x+/0wmUzYuHGjOpRWU1MDi8WyrHYTxQqGC8UNWZbR3d2tGeaamJhQ15goQ1y1tbVwOBwr/vuWurfY6dOnlxws85FlGV1dXWhqalJ7NxMTEzAYDKisrFSH0erq6nT59yJaCwwXilmhUAgdHR1qkLS0tKjf8Dds2KD2TKqrq1ftG77b7caJEyfw7//+73N2RX7ttddw9OhRZGRk6Pp3CiHQ39+vBk1jY6PaIystLVWDpr6+nmtqKGYxXChmTE9P4/bt2+ow1+3btxEMBpGWloaNGzeqYbJhw4Y1r00IITA2NgaPxwOHwwGn07mmM7+Gh4c1PZu+vj4AQEFBgSZsFqslEa0lhgtFjdfrVXskzc3N6OzsRDgchsPhQG1trTrMVVZWxllVs7jdbk3NRtmuxul0amakrV+/nmFDUcFwoTUzPj6uqZf09PSoH4iRxfeSkhJ+IC6TEtRK7+bOnTsIhUKw2+3YtGmTGjaVlZVISeHyNlp9DBdaFUIIDA8Pa6YFDw0NAQAKCws1a0xyc3MZJjpThhiVtTatra2Ynp5Gamoqampq1OnPGzduhNlsjnZzKQExXEgXQgj09fWpvZLm5maMjY1BkiSUlpaqPZOamhpkZWVFu7lJR5kcEVm38Xq9SElJQVVVlVq3qa2thc1mi3ZzKQEwXOihhMNh3L17V+2VtLa2wuPxwGg0orKyUu2V1NTU8MMqBgkh0N3drZmRNj4+rk7rVuo2dXV1PACNHgrDhZZEWfin9Era2toQCATUhX9KmFRXVyM1NTXazaVlEkJgaGhIHUZrampShzFLSko0a23y8vKi3FqKBwwXmpff70dbW5vaM1E2eLRaraipqVGHuSoqKlggTlAul0udjdbY2Iienh4AQG5urmb6c3FxMWtmNAfDhQDcP4cksl6ibPCYkZExZ4NHg4HHACUjj8ejqdl0dHRAlmVkZGRoejbl5eX8HSGGS7KK3OCxpaVFXZSXm5urmRZcWFjIb6U0L7/fP2dDzmAwCKvVqk5/rqur44acSYrhkgSUDR4jw2R4eBjA/fH0yD25cnJyotxailczMzPqhpyNjY1obm5GIBCA2WxGdXW1Ov25pqYGaWlp0W4urTKGSwKSZRk9PT2aPbncbjckSUJ5ebkaJDU1Nbrvi0WkCIfD6Orq0kwS8Hg8MBgM2LBhgzqMtmnTJm7ImYAYLgkgFAqhs7NTEyY+nw8pKSnqBo9KmHALd4oWIQR6e3s1kwRcLhcAoKysTFO3cTqdUW4trRTDJQ4pq6+VIFEOn1JWXyvDXBs2bODqa4pZyi4OStA0NjZiYGAAwP1dHCLDhhtyxh+GS5y4fv06mpqa0NzcjI6ODoTDYdhsNs1MLh6bS/FufHxcMyPt7t27mg056+vrsXv3bh41EAcYLnHi5s2bCAaDSE9Ph8PhQHp6OiwWC7/NUUILhUKYmpqCx+OBx+OB1+tFVVUVdw2IAwyXOBEKhWA0GhkmlNTC4TAkSeI6mjjAcCEiIt0x/omISHfcFOrPwuEwbt++Db/fH+2mxL2KigqOicepcDiM1tZW+Hy+aDcl7m3YsCGpj5dguPxZKBTC22+/jZmZGa4eXoHBwUH88z//M3bs2BHtptBDCIVC+MUvfoGZmRmuiXpIQggMDAzgX/7lX7Br165oNydqGC4RJEnCK6+8grq6umg3JW74/X7cunULNTU1sFqt+Ld/+zewjBe/hBAwGAz46U9/is2bN0e7OXFDWXuWn58Pp9OJf/3Xf03654DhMovBYOAme4sQQiAUCsHlcuHWrVv4+OOP0dXVhZdeegmHDh3ibLYEwedgcUIITE9PY2BgADdu3MClS5fQ1dWFv/qrv8Jf//Vf8zkAw4WWQHmQenp6cP36ddy6dQu9vb3w+Xzqt7MPP/wQjz/+eJRbSrQ6lN9zn8+H9vZ2NDQ04JtvvkF/fz+mp6fV133xxRc4dOhQtJoZUxguNC8hBGZmZtDb24srV67g6tWr6O/vRzAYnPf14+PjaG1tXeNWEq2eyEDp6OjAV199hWvXrmFoaAihUGje9wwPD6OtrW0tmxmzGC6kEkIgGAyiv78fDQ0NaGhoQG9vL2ZmZhZ8j8ViQW1tLZ5//nlUV1fj/Pnza9hiIv0JIeD3+9HR0YHLly/j6tWrGBwcXDBQJEmC3W5HXV0dvv3tb6Ourg5/+MMf1rjVsYfhkuSEEJBlGS6XC9euXcMXX3yBrq4uBAKBeV8vSRIyMjJQUVGBzZs3Y8uWLSgqKkJKSsqCDx9RrFN66l1dXfjqq69w5coV9Pf3L/g7bTAY4HQ6UVtbi+3bt6Ourg65ubkwGo18Dv6M4ZKkhBDq0caXLl1CU1MTpqam5p3hIkkSMjMzUVtbi927d2Pjxo3IzMyEwWBg4ZLilvLFanR0FJcvX8Znn32Gzs7OBXvqBoMBubm52Lp1K/bs2YOqqiqkp6dDkiQ+B/NguCQRZaZXd3c3Pv/8c3z99dcYGRmBLMtzXitJEtLT01FTU4Ndu3ahrq6OgUIJQRn2amtrwyeffIJr165hYmJi3i9WBoMB2dnZaqDU1NTA4XAwUJaA4ZIElF7KjRs38Mknn+D27duaGS6RrFYrqqursXfvXmzevBlOp5OBQnFP6aUMDw/j888/x4ULF9DT0zPvEJYkScjOzsbmzZuxd+9e1NbWsofyEBguCSryYbp48SIuXryIoaGheb+dmUwmrFu3Drt378aOHTtQWFjIHZgpISi1lLa2Npw/fx5Xr17F5OTkvK9Vzkfat28ftmzZgoyMDAbKCjBcEowy9NXV1YWPPvoIX3/9NTwez5zXKd/Otm/fjr1796KiogJpaWl8kCghCCHg9XrR0NCAs2fPoq2tbd5p9EajEcXFxXj88cfxxBNPoLi4mF+sdMJwSRDKQsempiacO3cOTU1N8w59mc1mVFdX4+mnn8YjjzzCb2eUUIQQcLvduHTpEs6ePYve3t55a4oOhwNbt27FM888g/r6elitVj4DOmO4xDkhBAKBAK5fv473338f7e3t844jp6enY8eOHXjmmWdQUVEBk8nEh4kShhACo6Oj+OSTT/DBBx/MOwRsMBhQUlKCp556Ck8++SQKCgpYT1xFDJc4pYwlX79+He+99x7a29sRDoc1r5EkCQUFBerDlJuby4eJEooQAi6XCx9++CE++OADuFyuOaFiNptRW1uLgwcP4tFHH4XNZuMzsAYYLnFGqak0NTXh97//PVpaWuaEitFoREVFBf7iL/4Cjz32mDp1kihRKMNf58+fx5/+9CeMjIzMeY3VasX27dvx3HPPYePGjeytrzGGS5xQZn91dXXhd7/7Ha5duzanQGkymVBTU4PvfOc72Lx5Mwv0lHCUNSoXL17E6dOnMTAwMKenYrfbsXfvXjz33HMoKytjgT5KGC5xQAiB4eFh/N///R8+++yzOacEGo1G1NbW4vDhw6irq4PZbObDRAlFCIFwOIxbt27h5MmTaG1tnVOodzgceOqpp/Dcc8+hpKSEQ8BRxnCJYcq3tE8//RR//OMf4XK5ND83GAzYsGEDvvvd7+KRRx5hqFBCEkJgaGgIv/nNb3Dx4sU5syBtNhueeuop/OVf/iWKi4sZKjGC4RKDlCGw27dvz/stTZIkFBcX4/Dhw9i9ezcsFgsfJko4yvT6jz/+GO++++6cuorZbMauXbvw4osvory8nKESYxguMUbZquW9997DuXPn5gyBZWVl4eDBg9i/f7+6JQVRolHqiydOnMCNGzc0k1YkSUJNTQ1eeuklbNmyBSkpKXwOYhDDJYbIsozW1lacOHECHR0dmkKl2WzG3r178cILL6CwsBAGgyGKLSVaHUpv5ezZs3j33XcxMTGh+XlOTg5eeOEF7N+/nwsfYxzDJQYoCyHPnj2L3//+9/B6verPJElCeXm5+i2NM18oUQkh0NfXh1/+8pe4evWqZijYZDLhiSeewEsvvYSioiI+A3GA4RJlykywX/3qV3MeKKvVioMHD+LQoUNcq0IJLRwO48svv8Qvf/lLDA8Pa35WUlKCI0eOYOfOnRwCiyMMlyiSZRmNjY34f//v/6G/v1+9L0kSKioqcOTIEdTW1sJoNEaxlUSrR5kR+b//+7/4wx/+oJkJZjab8cwzz+Dll19GdnY2QyXOMFyiJBgM4qOPPsLJkyc1w2Amkwn79+/H97//fRbsKaEpW7f8/Oc/x+XLlzU1xvz8fBw9ehR79uxhbyVOMVzWmPJN7be//S3OnTunWWWflZWFI0eOYO/evaytUEITQqCzsxNvvvkm2tvb1fuSJGHbtm34yU9+gpKSEj4DcYzhsoaEEJiYmMB//dd/zfmmVl1djZ/85CcoLy/nA0UJTZZl3Lx5E2+++Sbu3bun3jebzTh8+DB+8IMfcCZYAmC4rBEhBEZGRvDWW2+hsbFRvW8wGLB3714cPXoUmZmZfKAoocmyjC+//BJvvfWWZppxRkYG/vZv/xZPP/00UlL4sZQI+P/iGlC2rzh+/LhmCMBkMuHw4cN44YUXkJqaymChhBYOh3HhwgX8x3/8h6bOWFhYiFdffRX19fVcv5VAGC6rTAiBgYEBHD9+HJ2dnep9i8WCI0eO4Fvf+ha/qVHCC4fD+PTTT/Hzn/9cs+vEhg0b8Prrr2P9+vX8cpVg+Km2ipQ1LLODxeFw4Kc//Sl27drFacaU8GRZxqVLl/CLX/xCEyz19fV4/fXXkZ+fz2BJQAyXVSKEwPj4OH72s59pgiUjIwP/+I//iEcffZRDAJTwhBC4cuXKnKGwrVu34o033uD6lQTGcFkFQgj4fD7853/+J1paWtT7GRkZ+Kd/+ic8+uijfKAo4Qkh0NTUhJ/97GfweDzq/S1btjBYkgC/Oq+CUCiEkydPoqGhQb1ns9nw93//93jkkUf4QFHCU/YJe/PNNzE+Pq7er62tZbAkCYaLzsLhMD744AOcP39eXceSmpqKo0ePYseOHRwKo6QwOTmJt956S7OtUVlZGd544w3k5OQwWJIAP+l0JITArVu38Nvf/lY9f8JoNOJ73/senn76aQYLJYVgMIj//u//1qznys3Nxeuvv47CwkIGS5Lgp51OlJlhb7/9tlq4lCQJTz31FA4dOsRZYZQUhBC4cOGCpudutVrxD//wD6isrGSwJBGGi05mZmbw61//WjMMsHHjRvz4xz+G2WyOYsuI1oYQAr29vfj1r3+t7plnNBrxwx/+EI899hiDJckwXHQgyzI+/PBDXLlyRb3ndDrxd3/3d9zZmJLG9PQ0fvWrX8Hlcqn39uzZw557kmK4rJAQAj09Pfjd736n1llMJhP+5m/+BmVlZQwWSgpCCHz88ce4evWqeq+wsBCvvPIKUlNTo9gyihaGywrNzMzgf/7nf+B2u9V7Tz31FPbu3ctgoaSgbHEUOZHFbDbjyJEjXH2fxBguKyCEwOXLl3H9+nX1XlFREX7wgx9wvzBKGqFQCL/5zW80w2FPPPEE9uzZw2BJYgyXFXC73Th9+rT6bS0lJQU/+tGPuECMkoYQAjdu3MClS5fUe3l5eXj55Zf5BSvJMVweklLEj5wdtn37duzYsYPBQknD7/fj3XffxczMDID75xN9//vfR0FBAZ+DJMdweQjKwV8ffvihOpffbrfje9/7HkwmU5RbR7Q2hBD4/PPP0draqt6rq6vDvn37GCzEcHkYQgicP38eo6Oj6r2nn36as8MoqUxNTeG9996DLMsA7m9z9MMf/hAWiyXKLaNYwHBZJmUl/meffabey8rKwsGDB7m9CyUNIQS++uor3L17V7332GOPob6+nl+wCADDZdmU+fyRO73u27ePUy4pqfj9fvzpT39Sey1paWk4fPgwi/ikYrgs09jYGC5cuKBeZ2Vl4dvf/jZ7LZQ0lA1aIw/Be+SRR1BdXc0vWKTiJ+IyCCHw5Zdfamotjz/+OPLy8qLYKqK1FQ6H8dFHHyEUCgG4PwX/4MGD7LWQBsNlGXw+Hz799FPNDLFvfetb/LZGSWVwcBA3b95UrysrK1FXV8fngDQYLkukHNna19en3tu2bRuKi4v5UFHSUHrvU1NTAO4fK7Fv3z6kpaVFuWUUaxguSyTLMi5cuKDZnHLfvn2stVBSCQQC+OKLL9TrjIwMLhymefGTcYmGh4fR1NSkXpeVlbGASUlFCIHOzk709PSo97Zu3Yrc3NwotopiFcNlCYQQ+Oabb9ShAOD+ORXcSpySzeXLlzVbvTz++OPsvdO8+FuxBKFQCFeuXFEL+TabDdu2bWOvhZKKz+fTnNeSm5uLTZs28TmgeTFclsDlcmnm9FdWViI/Pz+KLSJae3fv3sXAwIB6vXXrVqSnp0exRRTLGC4PIIRAa2urZkhs+/btnNNPSUUIgatXryIYDAIAjEYjdu3axV4LLYjh8gCyLOObb75Rh8QsFgv3T6KkMz09rTkUz+l0ckILLYrh8gA+nw/t7e3qdVFREYfEKOkMDAygt7dXvd60aROHxGhRDJcH6O/v12z3UlNTw1lilFSUBcSBQADA/YWT27Zt4ywxWhR/OxYhhMCdO3fUcWZJkrBp06Yot4pobcmyrNnuxWazoaamhkNitCiGyyJkWUZbW5t6bbfbeSAYJZ2pqSnNbMmSkhIunKQHYrgsIhAIoLu7W70uKChAZmZm9BpEFAUDAwMYGxtTr2tqanicNz0Qw2URIyMjmoeqoqICZrM5ii0iWltCCLS3t88ZGmbvnR6E4bIAIQS6u7vVrS4kSUJVVVWUW0W0tpS6o8JqtWL9+vVRbBHFC4bLIjo7O9X1LWazGevXr+c3NkoqMzMzmo0qc3Nz4XQ6o9giihcMlwWEQiHcvXtXvc7MzEROTk70GkQUBR6PByMjI+r1unXrOBWfloThsgCfz4ehoSH1urCwEFarNYotIlp7IyMj8Hq96jV777RUDJcFuFwueDwe9Xr9+vVcNEZJRQiBvr4+hEIhAPfrjqWlpQwXWhJ+Ws5DCIH+/n7NDJmysrLoNoooCiK3fElNTUVhYWEUW0PxhOGygJ6eHrWYbzKZUFxczG9slFRkWdZssW+325GVlRXFFlE8YbjMQ5ZlzTc2u93OGTKUdILBoKaYn5WVxbojLRnDZR7T09O4d++eep2bmwubzRbFFhGtPZ/PB7fbrV7n5eXxHCNaMobLPCYnJzE+Pq5eFxUV8aGipDM5OQmfz6de5+fnc2iYlozhMg+Xy6VuLw7c36iPKNmMjY2pO1QA9/fWI1oqhss8BgYGEA6HAQAGg4HFfEpKw8PDkGUZwP3nIDc3l88BLRnDZR4jIyMwGo0A7m/7wpMnKRmNjo6qz4HJZEJ2dnaUW0TxhIWEeRw6dAjbt2/HwMAAxsfH+VBRUvrOd76DLVu2YGRkBG63G3l5edFuEsURhksEIQRu376N6elpAIDD4YDD4UBTU1OUWxYfwuEwpqamot0MWiFZltHa2qo+BzabDTabjc/BEoXDYc3uHsmK4fJnkiRh3bp1uHbtGq5duxbt5sQti8UCu90e7WbQQzIYDCgtLUVDQwMaGhqi3Zy4ZbVak/45kISyDD3JCSEQDofB/xwrZzQauQ9bnOJzoJ9kfw4YLkREpLvkjVUiIlo1DBciItIdw4WIiHTHcCEiIt0xXOKEsoZE2ZaGKBmFw2F4vV4+B3GA4RInuru78corr6C7uzvaTSGKmrt37+JHP/oR7t69G+2m0AMwXIiISHcMFyIi0h3DhYiIdMdwISIi3TFciIhIdwwXIiLSHcOFiIh0x3AhIiLdMVyIiEh3DBciItIdw4WIiHTHcCEiIt0xXIiISHcMFyIi0h3DJQ4IITA2Ngav14uxsTEIIaLdJKI1pzwHfr+fz0EcYLjEMLfbjePHj6Oqqgo7duzA2bNnsWPHDlRVVeH48eNwu93RbiLRqot8Dnbu3ImPP/4YO3fu5HMQ4yTB+I9J586dw4svvgifzwcAmm9pkiQBAKxWK06dOoUDBw5EpY1Eq43PQfxiuMSgc+fO4fnnn4cQArIsL/g6g8EASZJw5swZPliUcPgcxDeGS4xxu90oKSmB3+9f9IFSGAwGWCwW9PX1ITMzc/UbSLQG+BzEP9ZcYsw777wDn8+3pAcKAGRZhs/nw4kTJ1a5ZURrh89B/GPPJYYIIVBVVYXOzs5lzYSRJAkVFRW4c+eOOg5NFK/4HCQGhksMcblcyM3NXdH7s7OzdWwR0drjc5AYOCwWQ6amplb0fo/Ho1NLiKKHz0FiYLjEELvdvqL3OxwOnVpCFD18DhIDwyWGZGdno7KyctnjxZIkobKyEk6nc5VaRrR2+BwkBoZLDJEkCa+++upDvfe1115jEZMSAp+DxMCCfozh/H4iPgeJgD2XGJOZmYlTp05BkiQYDIv/36OsTD59+jQfKEoofA7iH8MlBh04cABnzpyBxWKBJElzuvnKPYvFgvfffx/PPvtslFpKtHr4HMQ3hkuMOnDgAPr6+nDs2DFUVFRoflZRUYFjx46hv7+fDxQlND4H8Ys1lzignGPh8XjgcDjgdDpZtKSkw+cgvjBciIhIdxwWIyIi3TFciIhIdwwXIiLSHcOFiIh0x3AhIiLdMVyIiEh3DBciItIdw4WIiHTHcCEiIt0xXIiISHcMFyIi0h3DhYiIdMdwISIi3TFciIhId/8fXgRsa3C2dpEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 500x400 with 4 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.999988712412588\n",
      "r2 is 0.9999928603717329\n",
      "r2 is 0.9968394556850537\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "tensor(0.9968, grad_fn=<SelectBackward0>)"
      ]
     },
     "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.46e-15 | test loss: 6.78e-16 | reg: 2.72e+00 : 100%|██| 20/20 [00:02<00:00,  8.21it/s]\n"
     ]
    }
   ],
   "source": [
    "model.train(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(2.01 x_{1} \\right)} + 2.0 \\log{\\left(0.62 x_{2} \\right)} + 5.85 \\right)}$"
      ],
      "text/plain": [
       "1.0*sin(2.0*log(2.01*x_1) + 2.0*log(0.62*x_2) + 5.85)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.symbolic_formula()[0][0]"
   ]
  },
  {
   "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": "stderr",
     "output_type": "stream",
     "text": [
      "train loss: 2.08e-01 | test loss: 2.18e-01 | reg: 2.75e+01 : 100%|██| 20/20 [00:38<00:00,  1.92s/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=5, k=3, seed=1)\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.train(dataset, opt=\"LBFGS\", steps=20, lamb=0.01, lamb_entropy=10.);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3a69ec41",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAFICAYAAACcDrP3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABIR0lEQVR4nO3de3QU53k/8O/salfalVaXXQmBEEjouhJISICxDQHLgM2tbhM7bdM2Td30NI3b2DltnDRt2ubXy8ml6TkxjpuLnbq1/0nS2OT0NBBszEUCDMExCAnQSkJCQgIh0N6099vM7w86m0Xosrua2Xln9XzO0T+wu/PuuzPzzHt7Xk4QBAGEEEKIhDRKF4AQQkj2oeBCCCFEchRcCCGESI6CCyGEEMlRcCGEECI5Ci6EEEIkR8GFEEKI5Ci4EEIIkRwFF0IIIZKj4EIIIURyFFwIIYRIjoILIYQQyVFwIYQQIjkKLoQQQiRHwYUQQojkcpQuACFqIAgC7HY7vF4vCgoKYLFYwHGc0sUihFnUciFkHi6XCwcOHEB9fT3KysqwZs0alJWVob6+HgcOHIDL5VK6iIQwiaOdKAmZ3TvvvINnnnkGfr8fwL3Wi0hstRiNRrz99tvYvXu3ImUkhFUUXAiZxTvvvIP9+/dDEATwPD/n6zQaDTiOw6FDhyjAEJKAggshM7hcLlRWViIQCMwbWEQajQYGgwHj4+MoLi6Wv4CEqACNuRAywxtvvAG/359UYAEAnufh9/vx5ptvylwyQtSDWi6EJBAEAfX19RgeHkYqlwbHcaipqcHg4CDNIiMEFFwIuc/U1BTKysoW9X6LxSJhiQhRJ+oWIySB1+td1Ps9Ho9EJSFE3Si4EJKgoKBgUe83mUwSlYQQdaPgQsj/EQQBZrMZNTU1Kb+X4zjU1tbCbDbLUDJC1IeCC1nSBEGI/wH3phW/8MILaX3OCy+8QIP5hPwfCi5kyZkZUDiOi/8JgoBPfepTMBqN0GiSvzzElgsh5B4KLmRJmBlQgF8HlcTXAEBJSQkOHjwIjuMWDDAajQYajQZbtmzBs88+i3/8x39EOByW50sQoiIUXEhWmyugzOy+SmzFAMDu3btx6NAhGAyGWV8v/pvBYMAvfvELdHV14atf/Sq+973vYd++fbh+/brM34wQtlFwIVlnvm6vuV4rvi7R7t27MT4+jpdeeumBQf6amhq89NJLuHnzJp588kloNBp87nOfw5EjR+B2u9HR0YGf/OQnKS3EJCSb0CJKkhVmnsbJDKzPFVTmeq3D4YDH44HJZILZbJ7zfV6vF1/+8pfx4x//GM888wz+7d/+jaYokyWHggtRrdlO3VSCSrKvT9fBgwfxhS98AWazGT/4wQ+wadMm2Y5FCGuoW4yoTrLjKHO9N/E9cnr66adx8uRJlJWVYf/+/fj2t7+NWCwm6zEJYQW1XIgqpNPtNddnZHotSiQSwbe+9S18+9vfxpYtW/D9738fK1asyGgZCMk0Ci6EWVIElJmfo+QixzNnzuCzn/0sgsEgXn75Zezdu1exshAiN+oWI0xJZaZXsp+X+DlK2rp1Kzo7O/Hoo4/iD//wD/HFL34RwWBQ0TIRIhdquRAmSNVKme3zlA4qMwmCgDfeeANf+cpXsGbNGrz66qtobm5WuliESIpaLkQxUrdSEj838fNYw3Ecnn32WRw7dgwajQa7du3Cf/zHf9CaGJJVqOVCMkrqFspcn89iUJlNKBTCV7/6Vfzwhz/E7t278fLLL9NmYyQrUHAhspM7oMw8hloCS6J33nkHzz//PHQ6Hb73ve9h+/btSheJkEWhbjEii2QSRUp5LDk/PxN2796Nrq4uNDY24plnnsE///M/IxKJKF0sQtJGLRciqUy0UmY7nlqDykw8z+OVV17B1772NbS2tuLVV19FdXW10sUiJGXUciGLJtfAfDLHFI+XLcTNyg4fPgyHw4GOjg789Kc/VbpYhKSMggtJixIBJfHYicfMRhs2bMDJkyexb98+PPfcc/jsZz8Lj8ejdLEISRp1i5GkpZsoUo7jZ2tQmc1bb72FF198ERaLBa+++io2btyodJEIWRC1XMiCFpMoUsoyJB57Kfn4xz+Ozs5OWCwW7N+/Hy+99BIlwCTMo5YLmVWmB+bnk41jK+mIRCL45je/iQMHDmDr1q343ve+RwkwCbMouJA4lgIKsHS7wRZy6tQpPPfccwiHw3j55ZexZ88epYtEyAOoW2yJU3JgfqFyJZaH/Nq2bdvQ1dWFhx9+GJ/85Cfx13/915QAkzCHWi5LFGutFBG1VpInCAL+8z//E3//93+PmpoavPrqq2hqalK6WIQAoJbLksJqK0VErZXUcByHT3/603jvvfcAALt27cLrr79OCTAJEyi4ZDnWA4qIBu3T19TUhKNHj+KTn/wkvvSlL+FTn/oUHA6H0sUiSxx1i2UhVru8ZkPdYNI6cuQIXnjhBej1enz/+9/HRz7yEaWLRJYoarlkkUwlipQKdYNJb8+ePejs7ERdXR0+9rGP4V/+5V8oASZRBLVcVE5NrZRE1A0mr1gshldeeQVf//rXsX79erz66quoqqpSulhkCaGWiwqpZRxlNtmacJI1Wq0Wn//853H48GFMTU3hsccew1tvvaV0scgSQsFFJdQcUETUDZZ5YgLMvXv34rOf/Sz+/M//HF6vV+likSWAusUYpnSiSKnQoD0bfvrTn+LFF19EWVkZXnvtNbS3tytdJJLFqOXCIBYSRUqFWivs+O3f/m10dnbCbDZj7969ePnll8HzvNLFIlmKWi6MUOvA/HxobIVNkUgE3/jGN/Dyyy9j27Zt+O53v4vly5crXSySZSi4KCgbAwpA3WBq0dXVheeeew6RSATf+c53sHv3bqWLRLIIdYtlWDYMzM+HusHUY/v27ejq6sJDDz2EP/iDP8CXv/xlhEIhpYtFsgS1XDIkW1spiagbTJ0EQcDrr7+Ov//7v0ddXR1effVVWK1WpYtFVI5aLjLK9laKiNauqBvHcfiTP/kTvPfee+B5Hjt37sR//dd/UQJMsigUXCS2VAKKiLrBskdzczPee+89/P7v/z5efPFF/NEf/RElwCRpo24xCSyFLq+ZaNA+ux0+fBif//znkZeXh+9///vYunWr0kUiKkMtl0VQW6JIqVBrJfvt27cPnZ2dqKmpwUc/+lF87WtfowSYJCXUcknRUmylJKKxlaUlFovh5Zdfxje+8Q20t7fjBz/4ASXAJEmhlksSlto4ymxo0H5p0mq1+Mu//EscOnQId+/eRUdHBw4ePKh0sYgKUHCZAwWUX6NuMLJp0yacOHECTz75JD7zmc/gc5/7HCXAJPOibrEE2ZIoUio0aE9mEgQB//3f/40vfelLWLZsGV577TW0tbUpXSzCIGq5ILsSRUqFWitkNhzH4Xd/93dx8uRJFBUVYc+ePXjllVcoASZ5wJJtuSz1gfn50NgKSUY4HMbXv/51fOc730FHRwf+/d//HeXl5UoXizBiSQUXCijzo24wko7Ozk4899xz8a2Vn3jiCaWLRBiQ9d1iNDCfHOoGI+l67LHHcOrUKWzYsAG/93u/h7/5m7+hBJgke1su1EpJHnWDESkIgoAf/vCH+OpXv4r6+nq89tpraGhoULpYRCFZ1XKhVkpqaO0KkRLHcfjTP/1THD16FJFIBDt27MCbb75JCTCXKNW3XKiFkh4KKkROgUAAf/d3f4c33ngDv/Ebv4Hvfve7MBqNSheLZFBWBBe6QaaO6o2kK5VbRjgcRigUQkFBQcrnG52f6pajdAEWi07A9FC9kXT967/+K9atWyfb5/v9fvh8Pjz77LOyHYPIT/XBhRCSWYODg3jxxRdl+ewTJ07gW9/6Fqqrqym4qFxWDegTQuTHcRy0Wq3kf36/H5/61Kfwla98RemvSCRAwYUQojhBEPCxj30MTz31FHbt2qV0cYgEqFuMEKK4t956C4ODgzhy5AiNB2YJCi6EEEW53W587nOfwy9+8Qvk5NAtKVtQtxghRDE8z+Opp57C7/zO76C9vV3p4hAJUXAhhCiC53l84QtfQDAYxEsvvUTdYVmGggshJOOcTieee+45dHZ24siRI9BqtUoXiUiMggshJGMikQjeeustPPHEE5iensbRo0dhNpuVLhaRAY2eEUIy4tq1a/irv/or3L59G3/7t3+L3/qt36IWSxajlgshRFaCIOCdd97Bnj170N7ejvfeew9PP/00BZYsRy0XQohsBEHAG2+8gX/4h3/AK6+8gqeeeooG7pcICi6EEFnEYjF885vfxOuvv46f/exn2LBhAwWWJYSCCyFEcna7HS+++CIuX76Md999FzU1NUoXiWQYjbkQQiQjCAKOHDmCXbt2ged5HD16lALLEkUtF0KIJMLhML7yla/g5z//Of7f//t/+PjHP06D9ksYBRdCyKJFo1F84hOfwPT0NN59911UVlbS+MoSR8GFELIogiDgj//4j+HxePC///u/MBgMSheJMICCCyFkUd544w2cOXMGvb29FFhIHAUXQkjaxsbG8KUvfQknT55Efn6+0sUhDKHZYoSQtASDQezZswdf/OIX0dzcrHRxCGMouBBCUnb79m18/OMfh9VqxYsvvqh0cQiDKLgQQlK2d+9elJaW4kc/+hHNCiOzojEXQkhKdDodPvnJT2Lt2rU4efKkLMdYtWqVLJ9LMocTBEFQuhCEEPUIBAKyH0Or1UKv18t+HCIfCi6EENkIgoBYLAatVkvdZ0tMVoy5UHxMHdUZyYTe3l4sX74cvb29SheFZFhWBBdCCCFsoeBCCCFEchRcCCGESI6CCyGEEMlRcCGEECI5Ci6EEEIkR8GFEEKI5Ci4EEIIkRwFF0IIIZKj4EIIIURyFFwIIYRIjoILIYQQyVFwIYQQIjkKLoQQQiRHwYUQQojkKLgQQgiRHAUXQgghkqPgQgghRHIUXAghhEiOggshhBDJUXAhhBAiOQouhBBCJEfBhRBCiOQouBBCCJEcBRdCCCGSo+BCCCFEchRcCCGESI6CCyGEEMlRcCGEECI5Ci6EEEIkR8GFEEKI5Ci4EEIIkRwFF0IIIZKj4EIIIURyFFwIIYRIjoILIYQQyVFwIYQQIjkKLoQQQiRHwYUQQojkKLgQQgiRHAUXQgghkqPgQgghRHIUXAghhEiOggshhBDJUXAhhBAiOQouhBBCJEfBhRBCiORylC7AYgiCALvdDo/HA5PJBIvFAo7jlC4W06jO0iPWm9frRUFBAdVbEgRBgNPpBM/zcDqdEASB6iwJ2XKuqbLl4nK5cODAAdTX16OsrAw1NTUoKytDfX09Dhw4AJfLpXQRmUN1lp6Z9bZmzRqqtwUk1llHRwecTic6OjqozhaQdeeaoDJHjhwR8vPzBY7jBI7jBADxP/Hf8vPzhSNHjihdVGZQnaWH6i11VGfpycZ6U1VwOXLkiKDVagWNRnNf5c/802g0glarVdUPIReqs/RQvaWO6iw92VpvnCAIgtStITm4XC5UVlYiEAiA5/kFX6/RaGAwGDA+Po7i4mL5C8ggqrP0UL2ljuosPdlcb6oZc3njjTfg9/uT+gEAgOd5+P1+vPnmmzKXjF1UZ+mheksd1Vl6srneVNFyEQQB9fX1GB4eRirF5TgONTU1GBwcVOVsi8WgOksP1VvqqM7Sk+31porgMjU1hbKyskW932KxSFgi9lGdpYfqLXVUZ+nJ9npTRbeY1+td1Ps9Ho9EJVEPqrP0UL2ljuosPdleb6oILgUFBYq+X40W+51NJpNEJVEXqrfUUZ2lJ9vrTRXBxWKxoLa2Nq3+xYqKCgwODmJwcDC+SngpWEydLV++HC6XK+lBxmxisVhQU1OT8vs4jkNtbS3MZrMMpWJXNBqFTqdDVVVVyu9dqnUm3FsCArPZjJqampSvUbXUmyqCC8dxeP7559N63+c///n4zfLy5cs4d+4c+vv7Ybfbs/rmmW6dAcBHP/pRdHV14Sc/+Qk+/PDDRTff1UQQBPzFX/xFWu97/vnnmR5glVIwGITD4cDdu3cRDAbx3HPPpfU5L7zwwpKoMzGgiA+3HMdBo9HghRdeSOvz1FBvqhjQB6SZD+7z+TA1NQW73Q6fzweNRgOz2QyLxQKz2YycHFWnWnuAy+XCypUrEQgEkmqxJdYZAPT392NwcBDRaBSVlZVoampCRUUF8yd1Oniej59X09PTWL16ddLnGsdxMBgMGB0dRUlJCTQaTdbWkd/vh9/vRywWg16vh9FoRF5eHtxud0rXp1hnN2/eZH69Rjpmu95mOydonQsDiouL8fbbb4PjuAUvXPHiPnjw4H0/QH5+PqqqqrBhwwZs2rQJVVVVCIVC6O/vx7lz53D58mXcvn0bkUhE5m+TGTk5OXjttdfiT0nzmVlnxcXFePjhh/GJT3wCW7Zsgc/nw7vvvou3334bvb29CAaDGfoW8hIEAbFYDDzPQ6PRICcnB2azOX6uJVNvGo0GP/vZz+Izd8TPyxbhcBhOpxOTk5Pwer3Izc1FaWkpLBYLDAYDOI677/pcqM7Ea/jrX/86DAZDhr6F/Ga2ToBff9e57lmp1Ntc9zVWqablInrnnXfw9NNPw+/3P/B/4g9oNBpx8OBBPPnkk0l9ZigUgt1ux9TUFNxuNwCgqKgofgHl5uZK9wUyxOPxxLOqvv/++3jmmWfidTbz5AeSq7O7d++ir68P169fB8dxWLNmDaxW66KmUypFEATwPB/P1Dtba+Odd95Jq94SW0Fi8FEbQRDirZRoNIqcnBwYjUYYjcZ5H+6SrbMf/ehHKC0tBQC0t7cjLy9Pxm8jn2RbKAuZr97Ez0z1vqY4ufLKyKmrq0v4sz/7M6Gmpua+3Du1tbXCgQMHBJfLlfZnh8NhYWJiQujt7RVOnToldHV1CRcvXhTGxsYEv98v4beQz/T0tHDr1i3B4/HE/83pdAoHDhwQamtrF11ngUBA6OnpEX76058Kr7/+uvA///M/Qn9/vxCJROT4OpKLxWJCJBIRotGowPP8vK9dTL2lchxWhMNhweVyCRMTE8KtW7cEh8MhhEKhlD4j2ToLBALC2bNnhTNnzgg+n0+OryMLnucf+JPCXPVWU1Oz6PuaElTXchEEAYcPH8bq1auxbt06OByO+N4kZrNZ0r7uaDQKh8MBu90Oh8MBnueRn58Pi8WC0tJS5OfnS3YsqUxPT8Pn86GwsHDW8gmCIFmdCYKAmzdvwmazYWxsDHq9HnV1dbBarSgqKlrsV5HcYloU6dabkEQLSWmCICAYDMLv9yMcDkOr1cZbKYtpdSVTZ+FwGN3d3YhEImhra2PymgKka6Eke6zEeisuLmbyvFmI6oLL5OQkzpw5gx07dmS03zEWi8HpdMJut8NutyMWi8FgMMQDDQtzzt1uN/x+P4qKimA0GjN6bK/Xi/7+fgwMDCAYDGLFihWwWq1YvXq14t1CiTd4pbqphP8b2wHY6SqLxWLxri+e55GbmxsfoM+kSCSC7u5uhEIhtLW1MbMubeatUambu/B/4zgsnDOpUF1w+fDDD2G32xXtdxQEAS6XKz7zLBKJQK/Xo7S0FKWlpSgsLMz4iehyuRAIBFBcXKzoIGksFsPo6ChsNhsmJydhNBrR0NCAxsbGjAc8FlsNLIzHhEIh+Hw+hEKh+Owjo9Go6GzJaDSKS5cuwe/3Y/369SgsLFSkHKwElEQUXDKA53kcOnQIdXV1aGpqUro4AO798NPT0/FAEwqFoNPpYLFYYLFY4k1aOY/vcrkQDAZRUlLC1MCo0+mEzWbD0NAQotEoVq9ejaamJqxYsUL2YyfexLVaLRM3iUSxWCyjQW/mNGKdTgej0Rif7cWCaDSKnp4eeL1etLa2ZqxngsWAMhPP80nNlGWJqoLLrVu3cO7cOTzxxBNMdEPNxuPxxANNIBCAVquNB5qSkhJotVrJjiX83x7l4XAYxcXFTAWWRJFIBENDQ+jr64PL5UJRUREaGxtRX18PvV4v6bFY7H6ay8yWlZTnhigcDsPv9yMQCMTXlhiNRuh0OsmPJYVYLIbe3l643W60traipKREluOoIaAkUmPrRVXB5fz58/B4PNi5c6fSRUmKz+eLT3EWF22WlJSgtLR00Ys2EwNLSUmJaqZLT05Ooq+vD6Ojo/E0FlarddHZXTNxo5aL1F1lgiAgEAjA5/PdN43YYDCo4ubE8zwuX74Mp9OJdevWSZL5N5MD8nKg4CKjaDSKQ4cOwWq1orGxUenipCwQCMQDjcfjiS88EwNNKk/w4mySSCSS8ntZEQgEMDg4CJvNBp/Ph7KyMlitVqxZsyblwJDpLia5LLYrLxqNwufzxTMy5OXlwWg0qubBIxHP87h69Sqmpqawdu3atNZSqT2gzKS2rjHVBJexsTF88MEH2LNnT8YHhqUmLtq02+1wuVwA7i3aFGeezXcz4HkeDocD0WhUtYElkSAIGBsbg81mw82bN5Gbm4v6+no0NjYuOKjLwuC41FKdhDBzGrFGo4lPI1ZT6202giCgr68Pd+7cQVNTE8rLy5N6z0xquRkvRG2tF9UEl7NnzyIUCqGjo0PpokgqEonEA42YtbmgoCA+8yxx5pcYWGKxGMxmM7P95umanp6O5zMLhUJYuXIlrFYrKisr77ugWJhaLLeFvuPMacR6vR75+fnIzc3NmpspcK8e+vv7MTExAavVOutkkGwOKIkouMggHA7j8OHDaGlpQW1trdLFkU00GoXT6cTU1FR80abRaIx3nQWDQfA8D4vFknVJNhPFYjFcv34dNpsNd+/eRX5+fnwCQG5ublZ0gSVrZussEonA7/cjGAzGU4IoPY1YboIgYHBwEDdv3kRDQwNWrly5ZALKTGrqGlNFcBkZGcHFixexd+9eZmdESY3n+XiguXv3LpxOJ/R6PWpqalBeXq7YOoBMs9vtsNlsGBwcBM/zqKqqQnNzc0amM7OC53n4fD54PJ74Ysf8/HymphFnwuDgIMbGxlBbW4vVq1cvqe8uUlPrRRXB5fTp0xAEAdu2bVO6KBkXi8XiCTV5nofL5Yov2hTHaIqKirL2QhOf3MPhMIaHhzEwMAC3243i4mJYrVbU1dVlXfegKBKJwOfzxTNQ5+bmIi8vDzqdbsm03GbenkZGRjAyMoI1a9agurpamUIpTMzgzTrmg0swGMQvfvELtLe3L7mTScxtxnEczGYztFptfNGmOPMsFAohJycnHmjkXrSZKfONOUxMTMBms2F0dBRarTaez0yuNRGZJE4j9vv9iEQis+b5yvYxp4XWoIyOjmJ4eBhVVVVp7RqqdmrpGmM+uAwNDaG3txf79u1T/cyoVESjUdjtdmg0GlgsljlvIF6vF1NTU5iamoov2jSbzSgtLZV80WYmpDJbyu/3x/OZ+f1+LFu2DE1NTaiqqlLd945Go/HFjoldXwvNHMyW2XKpLmocGxvDtWvXUFlZifr6ejmLxhy1dI0xH1xOnjwJvV6PLVu2KF2UjIlEInA4HPFAkexJ5Pf744FG6kWbmZDuOg+e53Hjxg3YbDZMTEwgLy8vns+MlSSIcxGnEYt5vtKZRqzGdT5SDMjfvHkTAwMDqKioQENDgyq+t1TU0HphOrj4/X4cOXIEDz30EFatWqV0cTIi3cAyUzAYjKehmZ6eji/aFFPRsNQKlPIJ3O12w2az4dq1awiHw1i1ahWsVitWrlzJzIUYi8XiXV8ztwtezPYHrCXpnEmOGV63b9+GzWZDeXk5rFYrc99ZLhRcFmlgYAB9fX3Yv38/80/dUgiHw3A4HNDpdJLuTRMOh+/baVMQBBQWFsZ32lRqBp6cKVui0SiGh4dhs9lgt9thMpni05mV+r7hcDg+QC9Xni/WcqtlYsrwnTt3cPXqVZSVlaG5uZnpG65U1NA1xnRwOXbsGEwmEzZv3qx0UWQXCoXgdDolDywziS2jqampBxZtWiyWjGU/yGRXzt27d2Gz2XD9+nUIghDfnnnZsmWyHVPE83y8lSLm+crENGIls0IrsQZlamoKly9fhsViwdq1a5m+6UqF9dYLs8HF4/Hg6NGjeOSRR1BRUaF0cWQlBha9Xo+SkpKMnSyxWOy+QBOLxWA0GuMzz+QYr1ByEDoUCsXzmXk8HpjNZlitVtTU1Eg+nVlc7JiY5ys/Pz+j3ZGZ7CpjYVGjw+FAb28viouLsW7dOtVN6kgVBZc09fX14dq1a9i3b19WnyTBYBBOpxN5eXkoLi5W7EQRF22KqWii0Sjy8vLigWaxizZZ6q4RBAG3bt2Kb8+ck5MTn868mD1ExDxfPp9vzmnESpBr6jKLaeudTid6e3thMpnQ2tqa1fcO1rvGmA0u7777LiwWCzZu3Kh0UWQTCATgcrlgMBgyumXzQsQNyMRAEw6H0160yfpA88ztmZcvX46mpqaUtmeOxWLxbMRKbhe8EClajSwGlJncbjcuXbqEgoICtLa2ZvV4LcutFyaDi8vlwvHjx7F169akMqGqkd/vh9vtZi6wzCQIQnwDtJmLNsUN0Oa6SYk3MxaDykw8z2N0dBR9fX2YnJyEwWBAY2MjGhoakJ+fP+t7Zk4jZmG74GSk+ruoIaDMND09jUuXLsFgMGD9+vVZm8WBgkuKLl++jJGREezfv5/JSlssMbAYjUYUFRUpXZyUiIs27XY7/H5/fMq0xWKJZxFQ++K+2bZnFjPyCoLwwHbB+fn5i5pGrIT5WpQsjJ9Iwev1oru7G3q9Hm1tbUxNv5cKy11jzAUXQRBw5MgRrFixAm1tbUoXR3I+nw/T09PIz89XffJJv98fn+Ls9XrBcRwKCwuT2pdGDcTtmcXszLm5uVi9ejVqampQVFTE9HbByRKDjJivKvEmpcaAMpPf78fFixeRk5ODtrY21Z+Ts2G19cJccLHb7ejs7MT27dtRWlqqdHEk5fV64fF4UFBQAJPJpHRxJCM+zYstGq/XC41Gg6KiovgUZzU+NSZuFzwxMYHr169jcnISer0+PgFAzedo4qU/s6uMxSfhdAUCAXR3d4PjOLS1tTE3FrZY4u9IwWUB3d3dmJiYwJ49e5irrMXweDzwer0wmUzMpyRJxWz99+KiTXGnTXHRptiiYf3inm+74JnbM5eWlsa3Z2Z9rAVYuMtLTeNkqQgGg+ju7gbP82hvb79vEz61Y7VrjKngIggCDh8+jNWrV6OlpUXp4khmenoaPp8PhYWFcw4Oq02y4ypiAk5xp02e55Gfnx/faZOVLatT3S5YEASMj4/DZrNhfHwcer0e9fX1sFqtzHV3pjqGwvoMv3SFQiF0d3cjGo2ira0ta65FgM2uMaaCy507d3D69Gk8/vjjWZE+Hbg3LdLv98f76NVuMSlbYrHYfTttxmIxGAyGeKBRokU313bBqbSuPB5PfDpzKBRCRUUFmpqaHtieOZOkmOHF0tokqYTDYXR3dyMcDqOtrS1rehFYbL0wFVw+/PBD2O12PPnkk0oXRRIulwuBQADFxcWqb4ZL/TQrbnwmjtNEo1Hk5ubGx2gKCwtlfQoLhUKSbxcci8UwMjKCvr6++PbMDQ0NaGhoyMiDhVxThtU++2+mSCSCS5cuIRAIYP369cy1NNNBwWUePM/j0KFDqK2tRXNzs9LFWTSn04lgMJgVgUXuPFWCIMDtdscDTTgchk6nu28DNCmOKeb58vl88WnERqNRljxf4vbMw8PD8e2ZrVYrli9fLulxMrkGRY2p/ecSjUbR09MDr9eL9evXq25JwGxY6xpjJrhMTEzg7Nmz2LVrl6qfJMTV7cFgECUlJcwPXs9HqW6RxJ02g8EgcnJy7tsALdVyzNwuWMzzlYlpxOFwGENDQ+jr67tve+ba2tq0Z9ApuahRzkzWmRaLxdDT0wOPx4OWlhbVd8Wz1nphJricP38eHo8HO3fuVLooaRMEAU6nE+FwGCUlJaqdU8/SgK7P54tnB/D7/dBoNPFAIy7anM1s2wWL2YiVuvjE7Zlv3LgBjUaD2tpaWK1WmM3med/H4qLGbOkq43kevb29cLlcWLduHSwWi9JFShsFl1lEo1EcOnQIVqsVjY2NShcnLYIgwOFwIBKJwGw2q3JdB8D2VNRAIBAPNOKizcSdNnU6XXy7YL/f/8A0Ylb4/X4MDAygv78/vj2z1WpFdXV1PFiyGFBmo2Rqf6nwPI8rV67Abrdj7dq1KCsrU7pIaWOpa4yJ4DI2NoYPPvgAu3fvVuX0wGwILGp7Eg2FQvExGrfbjVAohNzc3PjeNMXFxSlvF5xpPM9jbGwMNpsNt27dQm5uLurr69HY2BhfZMvCTWIhLLV00yUIAq5evYq7d++iubk5I3v9yIGl1gsTweXs2bMIhULo6OhQuigp43k+Pq1WfHpWE7nSsWeCOI3Y5XLB4XDA4/EgFApBp9PFd9pkfdGmePm53W709/djcHAQkUgElZWVsFqtqKysVM2NWs3nEnCv/DabDbdv30ZTU5Pkky8ygYJLgkgkgkOHDmHdunWoq6tTsigp43kedrsdPM+rLrCo+Wlz5jRig8GA/Px85OTkIBqN3rcBWuKiTYvFwkTLeL4ur2g0iuvXr6Ovrw92ux0FBQXx7MwsB8lEamsFJxIEAQMDA7h16xYaGhqwcuVKpYuUMla6xhQPLiMjI7hw4QL27dunmosH+PUujjzPw2KxqCL1h4jlcZW5pLNdsLhoU8wQIC7aFKc4ZzK/WzozvKampuLTmYUMb88sBTWeZ6LBwUGMj4+jrq4Oq1atUro4KWGl9aJ4cDl9+jQEQcC2bduULEZKYrEY7HY7AMBsNqsmsKjxiXLmdsHinimpjmuJizbFQBOJRJCbm3vfTptyrN9JlO7nh0IhXLt2DTabDdPT0ygpKUFTU5Ms2zNLTc0t5OHhYYyOjqKmpgZVVVVKFyclYpZrJSkaXEKhEA4fPoz29nZUV1crVYyUiN0uAGCxWJgeMBaprS9czPMl13bBgiBgeno6PvNMykWbcq5BEWTanjkT1JpKZmRkBNevX0dVVRVqamqULk7SWOgaUzS4DA0NoaenB/v371fFDCsxCaO41oL1wKK2p0ZxGnGmtwtO3GkzcdGmuNPmQr+zEosafT5fPJ9ZIBDA8uXLYbVaUVVVxfSNW42t5xs3bmBoaAirVq1SzbgwC11jigaXzs5O6HQ6bNmyRakiJC0SicDhcECj0cBisTB/Uahp/QFL2wWLizbtdjt8Pl/8QULcaTMnJ4epNSji9sziLCeDwYCGhgY0NjYyMXlhLmpLJXPz5k0MDAxg5cqVqK+vZ768gPKtF8WCi9/vx5EjR7Bp0yasXr1aiSIkTQws4pa+LAcWtXQ/8DzP/HbBgUAgnoZmenoaHMehuLiY2Q3QXC4XbDYbrl27hmg0ilWrVsFqtaKiooKZOk2ktlQyYoYFsZXIYp0mWrLBZWBgAH19fdi/fz/TA+LhcBgOhwM6nQ5ms5nZE0otF2o4HI53fYnTiFndLjjx0giFQvHJAG63GwDu22mTpQwAkUgEw8PD6Ovrg9PpRGFhIRobG1FfX89UOUVqeSACgMnJSfT19aGsrAzNzc3M3g8A5bvGFAsux48fR35+Ph5++GElDp8UtQQW1rsYxDxfPp8vPo1YzEbM2o0kmS6vSCQSb9GIO22aTKZ4oGEpC/adO3dgs9lw/fp1cByHmpoaZrdnVktX7t27d3HlyhWUlpaiubmZuXM4kZKtF0WCi8fjwdGjR/HII4+goqIi04dPSigUgtPphF6vR0lJCZMnOuuDozOnEYvZiFnrTlrMGIo4e9But8fXPeXn58dnnrEy7hEMBuPbM3u9XlgsFjQ1NTG3PbNaJqHY7XZcvnwZJSUlWLduHXPXnmjJBZe+vj4MDg5i//79THbfBINBOJ1O5OXlSbaXiJRYnlosTiMWtwsWpxEbDAamfms5ZnjxPH/fTpvRaBR5eXnxNDSZXLQ5F0El2zOzfI6LnE4nent7UVhYiJaWFqbOb5GSXWOKBJejR4+ipKQEmzZtyvShFxQIBOByuZCXl8fc/g4sP9XN3C44U9OIU5HJKcPivj7izLNIJAK9Xh/vOisqKlL8txO3Zx4cHEQwGERFRQWsVitWrVrFzM2c9da5y+VCT08PCgoK0NraylQrUKRU6yXjwcXtduPYsWPYunUrysvLM3noBYmBxWAwMLcojdVUGuLK91AoJNl2wVJScmOtxDKIizbtdns8uWZpaSkT6ybE7ZltNhvu3LkDo9GI9vZ2NDQ0KF20OFbPf+De5naXLl2C0WhEa2src5NTVB1c+vr6km5Se71ejI+Po7GxMekvy/M8otEo1qxZk3LZvF5v0je6YDCIcDicUveA2OxMZ4908Yksldem8+SWzns8Hk9SF4n4hJ6Xl5fSNOLF1FsqxFZeJty4cSOp7xMIBOB2u8HzfEpjjuKsqnSy9XZ3dyfVEvf5fJiYmEBBQUFKx4lGo4jFYmkFpGSvA7Hlnm5wSfU6GBkZSXrMLBAI4NatW/ftybMQnucRi8XSGndO5bYtvjadOlvMtSPJ4+X09HRKm3zV19cn/VqXy4VLly7BYDCkFVyi0WjSJ0iq0zTFPnZxXCFVqfSFptOfG4vF0j45otEoCgoKknptqje7WCwGl8uVdr2lIpNPa36/H5WVlUm9NtVsu4FAAIODg8jNzU0ruLhcLrS2tib1WqvVmtJnT05Oore3F/n5+WkFl2SvA7HVkqp0rwOfz5fSGrxkf3vg3sPbtWvXkJubK/ukJqVaeZIEl3R/9IUIgoAzZ85g7dq1uHPnTtqfI1flut3uRQ/iyVU2KXo75Sqby+WCTqdLqeWmFnJdB1euXEFFRQU8Hk/anyNH2SKRCI4dO4ZHH30UN27cSPtzWL0O5Po9e3t7UVNTE89TmA5W60zE1ujYDL29vdDr9Wm1WOQmCALC4TCKioqULsqsWL1xC4KASCTC1Mwk1jmdTkSjUeb2FhEEAceOHUN5eTmzSR1ZvA7Gxsag0WiwYsUKpYsiK2aDi9frxbVr17B9+3amBu9E4opiFssGpNbllkniwD+r9cYacXdEFtONDA4Owm634/HHH2eubCLWrgOe5zE8PIz169czW2dSYafWEwiCgK6uLjQ3NzO12jmR0+lEfn4+0ycIi2VzuVzUaknBjRs3oNPpmJsW7/f7cfbsWezevZuZmYFzYek6uH79ejzlUbZjLrgIgoCenh5otdqUJglkkjhrJ9kBb3KPOEOM1QcG1kSjUYyNjaGlpYWpG6QgCDhy5AisVissFovSxVENQRAwNjaG1tZWpn5PuTAXXJxOJ4aHh5ntDgPu5Ryjrp3UUZdY8sRBfNaSYgJAT08PIpEIHnroIfotUyBmRGBpYbGcmAousVgMp0+fRltbG9NPt263m4lUHnNZzLx2OVGXWPLcbje8Xm9K68Eywefz4eLFi9i7dy9TYxmzYek6EAQBw8PDzLVC5cRMZ6kgCDh37hyKi4uZ3vJYXMjFcvBjEXWJJY/neVy9ehUNDQ1M3cDF7rC1a9fSQ0KKxCnkUnelJ04bZi1oMXPmjoyMwG63Y8uWLcxVUiKfz8d0OnBAunnqUqKuxOQIgoCBgQEYjUbm0uLbbDaEw2Fs3LhR6aIkhaXr4PLly6itrZU9nx1LmAguHo8H3d3d+MhHPsL0zBNBEOD1epnLOzaTmEtosWKxGCKRiCQnrcvlYrorkRUulwt2ux3r1q1jKhCHQiH88pe/xJNPPslUa2o+YqoYpcViMYTD4YysU0onLYxcFL+Tx2IxdHV1wWq1MjfdciZxbQvLAVAkxUU1NTUFjuOg0+lQVFSU9meKXYlLYfrlYkQiEVy9ehWNjY1MnWOCIODo0aNYs2YNzGaz0sVJCQsBemhoCCaTSfKypBscEt8nZ+49RcO6OM5iMpmYXCQ2kxrWtkippKQEZrMZPM/D4XCkfTJTl9jCxJQgZWVlzE3vHR0dhdPpxLZt2+g3TJEgCLh16xbWrl0r2zEWe20tZrO8+SgWXMS+ZafTyfw4C/DrDKasr22Rsqmr1+uRk5MTDzDT09NpfT7NEpufeC0IgoD6+nqmroVIJIKuri7s3LmTiS6mZLEyU8zr9YLjuIxOJ5/vGl3o+pXy/qHY2WK323H16lVs376dqS6A2QiCEN+ZUumTNVlSlpPjOFgsFvj9fkQikZTeS7Pr5icurHM6ncwtrhNzh1VUVGR9Hiy5XLlyBdXV1Uz9riK5y6RIcAkEAjhz5gw2btyoiidacWBbjiSV4hTdZJ4YknmtXIN0Wq0WRUVFKXeP+f1+5jZ3YoUgCLh9+zbGx8eZ3GRqeHgYU1NT6OjoUN3vx0LCSp7nEQwGsWrVqowcb6HfaL5py3LcNzIeXKLRKDo7O1FdXZ2xSl8MQRDgcDhkG2tJ3Lp4oaAhds2Jr5+NnDNkDAYDtFpt0t1j4g6MrE/UUIIgCJiYmMD169exbt065iY7+Hw+nD59Gk888QRzPQvJnntKd+MNDw8jPz9f1nLMdU9KNljI+dCQ0drneR6nT5+GyWRipgsgsTUw8wcRBAEejwccx8k21iJOGxa3cZ2tDGIuM+DXs8BisdicJ5Bc9cpxHMxmM/x+P6LR6IKvF8vM2hO50gRBwI0bNzA6Oop169YxN0U7Fovh8OHDaGpqwrJlyzJ6bPEhauZf4r8v9EAmx3hL4jGTff34+DhaWlokK8NM6X6/TN13M/ZIwvM8zp49C57n8cgjjygaWMQTJRQKIRKJxJ/2dToddDodNBoNBEGA3+9HIBBAaWmpbOUVNxvTarXxVknik87MrV05josHm5nTCDOxiEqr1cJkMsHhcGDZsmXzPjnJ2eJTq2g0iqGhIbjdbrS2tia9S2qm8DyP48ePw2g0Zjx3WOJD1FzndeK/L2bL41TK5HK5MDk5iWAwCIPBgIqKChQUFMx73MnJSeTk5GQ8L5x4f5gp2daelDISXCKRCH75y18iHA5j27Zti969cTEEQYDP54Pf74dWq43PiBL7R30+X/y1Wq0WFotF1vImnqBigBG3ZRV/7MTAIr5HfO1s0xDlviHk5+fH94EvKiqa9XhiFx5rT+VKEc+7wcFBaLVarF+/nrmElDzP48yZM/B4PHjqqacUCSypBAvxwWvmNSBly2VwcBB3795FeXk5li9fjunpaVy6dAl1dXUoLy+f9RiCIKC/v1/xPVtmW8MyW3nmCkiLJWtwEVe0nz9/HjqdDtu2bVO0iyQWi8HtdiMWi6GoqAh6vR7A/ZUrNr3FG3gmTw7xmGJzf2ZQmQ3P8/Hgl+5e4emU02w2486dO8jNzX1gFp3YajEajUu+1SIIAqLRKCYmJnDr1i0sX74cq1evVnw8IJEgCAgGg3j//fcxPT2NvXv3Zvw6TXxQSuWc4TgOsVjsvmtVynFHn8+H9vZ2GAwGcByH8vJyWCwWXLlyBYIgYPny5Q+c+6Ojo/GFxyxQKi2MLMFF3AL4xo0b6O/vx6pVq7Bu3TrFWiziBe50OqHX61FcXPzAySeeIFqtVtGWVWLrZL6LLLH1knjyZOqmpdVqYTab4XA4YDabodfr40Ha7/cjFoupYiagXMRrYGpqChMTE9DpdGhqakJhYSEzAVcs49jYGLq7u1FSUoJ9+/Yp0qLiOC7lLi7xPbN1J0tVx62trffdD8QHq5aWFvT29gJAPMCIY7QjIyPYvHmzYr9z4sNyOoFFqnJLFlxcLhei0Sj8fj8cDgcmJyeh1+uxadOmOZuPmRIKheB2u5Gfn6+KMYB0ntyUeBLW6/UoKiqC0+lEQUEBcnNzEQqF4PF4YLFYmK9nOXi9Xni9XrhcLng8Huj1eqxevVr27tVk2O32+IOWz+eLBz4AWL9+PWpraxVrUaU7djLzIUvsCpLq3JvtN+M4DsXFxWhpacGVK1fg8/mwbNky+P1+DA4Ooq6ujsl1XQs9rDI75nL+/HlwHBe/4axfvx6lpaXxH0fJjJ1utxuFhYXxTXpYyh662LJoNJp4V57UdZ3M54jTk71eb3xNi9lshk6nY6qeM8VmsyE3NxeFhYUPDPwqXR+dnZ0A7p0zubm5KC4uRnt7O1asWBHvBlOyjIs5tngdAIh3LctdpuLiYrS1tWF4eBh9fX3QarVoaGiIz7CTuy5T/fxUZrpJQZLgotfrUVdXB41Gg5ycnPjT9OTkpBQfDwBpZyLWaDTx1kooFJKsPIkW0z8txQ+ZOItstv9Lh0ajSam+jEZjvGtCEISk3sva+onF0ul0WLlyZfwaCIfDcDgckh4j3dlleXl58fxW4nUqtlLE1stiLSappRTXgfh9pLoOdDodpqamFnzd8uXLEY1GodVqodFoknqPiPV0UovBCRL8qplYDZtuUzdTT2JUtvRkU9cZy3WWqRXr6XSrsVpvdF9b3PUpSXBJhd/vx+joKGpqapibiikuDmR1IHq2qYVKE9cBFBQUMLlYksU6AwCHwwG/34+KigrmZo6JKeLLy8uVLs6sWPxNg8EgxsfHUVVVxdx1IMei0mRk/KwOBoMYGRnBu+++i5s3b2b68PPiOA4+nw/BYFDpojwgceEkSzweD0KhEHMXe7L52pQSiUQwOjqKixcv3re2Smkcx2FiYgJHjx6F1+tVujgPYPE64Hkely9fht1uV3zSxkxKXgcZDy5msxm7du3CsmXL8Mtf/hIXLlyIr8pVmsFgQF5eHtxuNxOJ71gXCoXg8/lgMpmYGj9JfFJjLeiJysvL0d7eDo7j0N3djVu3bildpLiHH34Yer0ep06dYuomzqrR0VF4vV40Nzcz0wpNTFej1HWgSE3odDo8/PDD2LBhA8bGxnD8+HG43W4livIAceGTy+VStiCM43kebrcbubm5TKUwUaoLIB1GoxFtbW1Yvnw5hoaGcOXKlZS3NJCDXq/H9u3bcfv2bVy5ckXp4jBtenoao6OjqK6uZiYbReKUbCUzkisaZqurq7Fjxw5oNBqcOHEC165dU7I4AO4NSBYVFSEUCsHv9ytdHGZNT0+D53mmViGrKbCINBoNamtrsXbtWng8Hly4cIGJB5vly5dj7dq1+PDDD+F0OpUuDpNisRiuXr0Kk8mEqqoqpYsDAIq3VhIp3oYzmUzo6OhATU0Nenp68P7778s2ZThZeXl5MBqNmJ6eZqbLjiXBYBCBQABFRUVM9DGroRtsIWazGRs2bIDRaERvby9GRkYU75IS91vq6uqibuJZDA0NIRQKoampSfHzTuwGA9JfkCo1xYMLcG/RU2trK7Zs2QKn04ljx47hzp07ipapsLAQGo2GiadIloj52fLy8phYhazG1spc9Ho9WlpasGbNGoyPj+PSpUuKTi7RarXYvn07XC4XLl68qFg5WORwOHDz5k3U1dUpvhfPzG4wVrBTEtxriu/cuROFhYU4ffo0ent7FXtiElM8hMNhJmfNKEUcG1O6O0yt3WDJqKysxPr16xGJRHDhwgVFH7QsFgva29vR29ur+AMfKyKRCPr6+mA2m7Fy5UpFy8JSN9hMTAUX4F6X1NatW9HS0oKhoSGcPHlSsZu7Xq9HQUEBPB4PEwOtSvP7/QiFQrMm/sykbOgGW4jJZMKGDRtgsVjQ39+P/v5+xbpoW1paUFZWhq6uLroOAAwMDEAQBFitVsXKwGI32EzMBRfg3k2jvr4eHR0diEajOH78OEZHRxUpi7g40OVyKd4HrqRoNIrp6WkYjUbFFr9mc2tlNlqtFo2NjWhsbITdbseFCxfg8XgyXg6O47B9+3YEAgGcP38+48dnyeTkJO7cuYOGhgbFrwPWusFmYrdkuJdPbMeOHVi5ciU+/PBDnD9/PuNPTmL3WCwWU+TCZoXL5YJWq1Use8FSaK3MZdmyZdiwYQN0Oh0uXbqEsbGxjD/omEwmbN68GQMDAxgbG8vosVkRCoUwMDCA8vLyjG//DLCxdiUVTAcX4F5yw40bN2Lz5s2YnJzEsWPHJE8GmEwZTCYTfD4fwuFwRo/NAq/Xi0gkguLiYkVO6KXUWplLXl4e1q9fj8rKSoyMjODy5csZPxcbGxuxatUqnDlzhsksFnISBOG+zMdKHJ+FtSupYD64iCorK7Fz507k5eWhs7MTNpsto09v+fn50Ov1S657LBKJwOPxKJI7bKl1gy2E4zhUV1ejpaUFfr8fFy5cgN1uz2gZtm7dCp7n8f7772f0uEq7efMmnE4nmpqaMp6NQk2tlUSqCS7AvRXNjz32GBobG9HX14dTp04hEAhk7PjFxcXgeR7T09MZO6aSxKSUOp0u46uPl3I32EKKi4uxYcMGFBYW4urVqxgaGsrYrEqDwYAtW7ZgdHSUiUXPmeD3+zE0NITKykqUlJRk7LhqGLSfj6qCC3DvZtPc3Ixt27bB5/Phvffey1heJnHMwe/3L4luAY/Hg1gslvZeOumg1kpydDodmpubUVdXh9u3b6O7uztjCTCrq6tRV1eHc+fOMZV0Uw6CIODq1avIy8tDTU1NRo+rhkH7+aiz1ABKS0uxc+dOlJWV4dy5c7h48WJGpmoajcYlkdxSiaSU1FpJ3YoVK9DW1gZBENDd3S3Zxl8LEZNbdnV1ZXU38cjISDwpZaayUai1G2wm1QYX4N46lEceeQTt7e24ceMGTpw4kZEEmOICQlaSbUpNiaSU1FpJX35+Ptrb27F8+XJcu3YNV69elX1WpV6vx7Zt23D79m1cvXpV1mMpJdNJKdXeDTaTqoOLaM2aNXj88cfBcRxOnDiBoaEhWY8nJrcMBoNZmdwyk0kpqRtMGmICzObmZkxPT2ckAeaKFSviyS2zLU1SppNSJrZW1NoNNlN2fAvcywXW0dGBNWvW4NKlS7InwMzW5JaZTEpJ3WDSs1gsGU2AuXHjRphMJnR2dmZVN3GmklJmW2slUdYEF+DegPv69euxZcsWOBwO2RNgZltyS7E7LBNJKam1Ih+9Xo9169ahurpa9gSY2ZjcMlNJKbNh0H4+2feNcC8B5q5du+IJMC9fvizLU1W2JbcUg6Sc3WHUDZYZHMdh1apVGUmAabFY0NbWlhXJLTOVlDJbBu3nk5XBBfh1Asx169ZhcHAQnZ2dsgQAMbmluIpdrTKRlJK6wTIvUwkwsyW55cDAAHiely0pZTZ3g82UtcEFuHcTa2hoQEdHByKRCI4fP44bN25IfpyCggLk5OSodvW+3EkpqbWirEwkwNRoNNi2bRsCgQA++OADST87U8SklI2NjbJeB9naDTZT9n9DACUlJfEEmL/61a/wwQcfSPp0xXEcioqKVJvcUs6klNRaYcfMBJjj4+OSfn5hYSE2b96M/v5+1SW3lDMppdoSTkplSQQX4NcJMB966CHcvn0bx48flzQBpk6nQ0FBgeqSW8qZlJJaK+xJTIB5/fp19Pb2Snq+NjY2orKyUlXJLeVMSqnGhJNSWTLBRbRq1Srs2LEDubm56OzsRH9/v2RdWQUFBapKbilXUkrqBmPbbAkwpXzQEpNbnj17VrLPlJOYlNJqtUqajWIptlYSLbngAtxb0bx9+3Y0NjbiypUrOH36tGQJMNWS3DIxKWVBQYGknwtQN5gaiAkwTSYTrly5IlkCTKPRiC1btmBkZET2Bc2LlZiU0mw2S/KZS2nQfj5LMrgA9350MQGm1+vFsWPHJMnLpJbklolJKaU6+am1oj46nQ5r165FbW1tPAGmFFknqqurUVtbi7NnzzKb3FKOpJRLbdB+Pkv72wMoKyvDzp07UVpairNnz6K7u3vRUzVZT24ZDoclTUpJ3WDqV1FREU+AefHiRUketB555BHo9XqcOnWKyW5iMSllU1OTJNkolno32ExLPrgA9yfAHB0dxYkTJxbdrcVqckue5+FyuaDX6yVJSkndYNlDTIBZXl4uSQJMMbnlxMQEc8ktE5NSLnaWJHWDzY6CSwIxASYAHD9+HMPDw2l/FqvJLcWklIvdo4VaK9lJo9Ggrq4Ozc3NcLvduHDhwqIekFhMbhmLxdDX1ydJUspsTDgpFaqNGQoLC/H4449jzZo16O7uxtmzZ9OeqslackupklJSayX7iQkwDQYDenp6FpUAU0xu2dXVxUQ38dDQEILB4KKSUlJrZWEUXGYhJsB89NFHYbfb8d5776WdM4mV5JZSJaWk1srSkZubi5aWlkUnwBSTWzqdTnR3d0tf0BSISSlra2vTTkpJg/bJoZqZx4oVK+5LgHnlypWUn7wSk1sqOWtmsUkpqRtsaZotAebdu3dT/hwxuWVPT49iyS0jkQhsNtuiklLSoH3yKLgsIDEB5sDAADo7O1MOEuLgucfjUSSp32KTUlI3GDGZTGhvb4fFYoHNZsPAwEDKXb0tLS0oLS3FqVOnEI1GZSrp3MQyW63WlM9j6gZLHQWXJCQmwAyHwzh27FjKCTBNJhO0Wi3cbndGp2UuJikltVZIopycnHgCzKmpKVy8eDGlTOMajQbbt2+H3+/PeHLLO3fupJ2UkrrB0kM1lYKSkhLs3LkTFRUVKSfAFLvHIpFIRvd+cbvd0Gg0KU+3pNYKmcuyZcvQ3t4OrVaL7u7ulBJgFhYW4qGHHoLNZpM8ceZcQqEQ+vv7sWzZspSSUi7VhJNSoeCSopycHGzatAkPPfQQJiYmUkqAqdPpYDKZ4PV6M5LcUjxOqqvwqbVCFmIwGNDW1oaVK1fi+vXruHz5ctLntNVqRWVlJU6fPi17FgtBEGCz2VJOSrmUE05KhYJLmlatWoWdO3dCr9enlAAzU8ktE5NS6vX6pN5D3WAkFRzHYc2aNWhpaYHP50spAWamklvevHkTDocDVqs16eSs1FqRBgWXRcjPz8djjz2GhoaGlBJgyp3cMp2klNQNRtKVTgLMTCS3FJNSrly5MqmklDRoLy0KLouk0Wiwdu1abNu2DR6PJ6kEmHInt0w1KSW1VshipZMAU0xuee7cOcmn6Yt7tOTl5aG2tjap19OgvbSoFiVSVlaGXbt2wWKxJJUAU67klqkkpaRuMCK1mQkwb9++Pe/rH3nkEeh0OsmTW46OjsLj8SSVlJK6weRBwUVCer0ejz76KNra2pJKgCl1ckuxOyyZpJTUDUbkkpgAc3BwEH19fXOua0lMbtnX1yfJ8aenpzEyMoKqqqp5Z0lSN5i8KLjIoKamJp4A88SJE3MmwExMbinFZmViK2i+pJTUWiGZICbAbGpqgsvlmjcB5ooVK9Dc3Ixf/epXi06TJCalLCgoQHV19Zyvo4ST8qNalYmYALOqqmreBJhiri+3272o5JbJJKWk1grJtNLSUmzYsAF5eXno6enB6OjorN1fmzZtkiS55fDwMILBIJqbm2c9x6m1kjkUXGSk1WrR1tZ2XwLM2fIyFRUVLSq5ZTJJKam1QpQiJsCsqqrC2NgYenp6HpjIIkVyS4fDgfHx8TmTUtKgfWZRDWfAihUrsHPnTphMJpw6deqBBJiLTW45X1JK6gYjLOA4DqtXr0ZrayvC4fCsCTAtFgvWr1+Pnp6elJNjRqNR2Gw2lJSUzJqUkgbtM4+CS4YYDAZ85CMfwdq1azEwMICurq77AklicstUkvrNl5SSusEIawoLC9He3g6z2TxrAszW1laUlpaiq6srpetA/JyZe7RQN5hyKLhkEMdxaGxsREdHB0KhEI4dO4axsbH4/4vJLZNdvR+LxWZNSkmtFcKynJwcWK1WNDQ0PJAAU0xu6fP5kk5ueefOHUxOTqKhoWHW64C6wZRBNa6AkpIS7NixAxUVFfjggw/wq1/9CtFodM7kloIgYGpqCiMjI5iamooHDpfL9UBSSmqtELUoLy+fNQFmYWEhNm/eDJvNhps3b8ZfP9t1kJiUsry8PP466gZTHidkMv87ecCNGzfQ3d2N3NxcbN68GSUlJfB6vfB4PNBqtfjRj36E73znO/elyKitrcVnPvMZ/OZv/iZqamriucOotULUSBAEjIyMYHx8HCUlJWhoaIBer8e7774Lp9OJjo4O/PjHP571Onj66aexc+dO7NixAzqd7r7WCl0HyqLgwgCfz4fz58/D5XJh7dq1qK+vx1tvvYVnn30Wfr8fHMfN2k1mNBpx8OBBPPnkk/F/owuKqJXT6cTAwAAEQUBDQwPy8vLwT//0Tzhw4ABCoRAAzHsdPPHEEwCo1c4KCi6M4HkefX196O/vx/Xr1/HFL34RwOwXk0i8iA4dOoQ9e/ZkqqiEyCYSiWBgYAAOhwNDQ0P49Kc/DSC56+DnP/859u7dm6mikgVQcGHM4OAgWlpaEA6HkxrU5zgORqMR4+Pj867MJ0RNrl69io0bNyIUCtF1oFI0oM+Yw4cPJ31BAfee6Px+P958802ZS0ZI5hw9ehTBYJCuAxWjlgtDBEFAfX09hoeHU8oQy3EcampqMDg4SH3NRPXoOsgOFFwYMjU1hbKyskW932KxSFgiQjKProPsQN1iDElc25IOj8cjUUkIUQ5dB9mBggtDkt2SeC4mk0mikhCiHLoOsgMFF4ZYLBbU1tam3F/McRxqa2uT2iecENbRdZAdKLgwhOM4PP/882m994UXXqBBTJIV6DrIDjSgzxiXy4XKykoEAoGkNk3SaDQwGAw0v59kFboO1I9aLowpLi7G22+/nVQmVzGF+MGDB+mCIlmFrgP1o+DCoN27d+PQoUMwGAyz5kkS/81gMODw4cP35RYjJFvQdaBuFFwYtXv3boyPj+Oll15CTU3Nff9XU1ODl156CTdv3qQLimQ1ug7Ui8ZcVEAQBDgcDng8HphMJpjNZhq0JEsOXQfqQsGFEEKI5KhbjBBCiOQouBBCCJEcBRdCCCGSo+BCCCFEchRcCCGESI6CCyGEEMlRcCGEECI5Ci6EEEIkR8GFEEKI5Ci4EEIIkRwFF0IIIZKj4EIIIURyFFwIIYRIjoILIYQQyf1/KU+FG93cIgEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 500x400 with 16 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ad2e8d6f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAFICAYAAACcDrP3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAmWUlEQVR4nO3dfWxc1Z3G8efOq2fssccvcZyQGuzUQKB0kaoVsC9dEAUjsarK8s+qS1sB21YUSNsVkDZQkrShbNvtkhQVSquiJStUlVUQK23opqJVS6Vtq7YSJUAbUlIS8uLEdjwztmfGM3Pv2T/SeztjO4kdX3vevh8Jqb7GcBru9TPnnN/5XcsYYwQAgI8C1R4AAKDxEC4AAN8RLgAA3xEuAADfES4AAN8RLgAA3xEuAADfES4AAN8RLgAA3xEuAADfES4AAN8RLgAA3xEuAADfES4AAN8RLgAA34WqPQCgHhhjND4+rqmpKbW1tam7u1uWZVV7WEDNYuYCnEUqldLOnTs1NDSkVatWaWBgQKtWrdLQ0JB27typVCpV7SECNcniTZTA/Pbu3atbb71V2WxW0unZi8udtcTjce3evVvDw8NVGSNQqwgXYB579+7VzTffLGOMHMc5498XCARkWZb27NlDwABlCBdgllQqpXXr1imXy501WFyBQECxWExHjhxRMplc/gECdYA9F2CWZ555RtlsdkHBIkmO4yibzWrXrl3LPDKgfjBzAcoYYzQ0NKSDBw9qMY+GZVkaHBzUgQMHqCIDRLgAFcbGxrRq1aol/Xx3d7ePIwLqE8tiQJmpqakl/fzk5KRPIwHqG+EClGlra1vSzycSCZ9GAtQ3wgUo093drfXr1y9638SyLK1fv15dXV3LNDKgvhAuQBnLsnTvvfee189u3LiRzXzgT9jQB2bhnAuwdMxcgFmSyaR2794ty7IUCJz9EXFP6D///PMEC1CGcAHmMTw8rD179igWi8myrDnLXe61WCymF198UTfeeGOVRgrUJsIFOIPh4WEdOXJEO3bs0ODgYMX3BgcHtWPHDh09epRgAebBnguwAMYY/exnP9OHPvQhvfDCC/rbv/1bNu+Bs2DmAiyAZVlKJpMKBAJKJpMEC3AOhAsAwHeECwDAd4QLAMB3hAsAwHeECwDAd4QLAMB3hAsAwHeECwDAd4QLAMB3hAsAwHeECwDAd4QLAMB3hAsAwHe03AcWqFQqKZ1Oq6OjQ6FQqNrDAWoa4QIskDFGtm0rGAzSch84B8IFAOA79lwAAL5j4Rg1w7Zt/epXv9Lk5GS1h1L3/uIv/kK9vb3VHgaaGMtiqBm5XE4333yz8vm82traqj2cumSM0cGDB/XNb35TN910U7WHgybGzAU1wxijQCCg7du366//+q+rPZy6YIzRD37wAw0MDGjDhg0qFov6x3/8R/GZEdVGuKDmhEIhRaPRag+j5hlj9Morr+jzn/+8gsGgNm3apFtuuYVKNtQEwgWoU+l0Wlu2bNHY2Jgk6XOf+5xs267yqIDTqBYD6pBt23riiSf0i1/8wrs2NDSk66+/voqjAv6McAHqjDFGL7/8sr797W/LcRxJUltbm7Zu3aq+vr4qjw44jXAB6ogxRidOnNC2bds0NTUlSQoEAvr4xz+u97///VUeHfBnhAtQR4rFor761a/qtdde865dffXV+tSnPqVgMFjFkQGVCBegThhj9D//8z967rnnvGs9PT3aunWrkslk9QYGzINwAeqAMUbHjh3TV77yFeXzeUmnS7Y/85nP6Morr6T8GDWHcAHqQKlU0mOPPaa33nrLu/aBD3xAH/nIRxQI8Bij9nBXAjXOGKOf/OQnFcthvb292rx5s+LxeBVHBpwZ4QLUuPHxcT366KPKZrOSpGAwqLvvvluXXnopy2GoWYQLUMNs29ZTTz2lffv2edeuueYaffSjH2U5DDWNuxOoUcYY/frXv9bTTz/tNaJMJpPavHkzXaNR8wgXoEZlMhk98sgjSqfTkiTLsnTnnXfqfe97H8thqHmEC1CDHMfRf/zHf1T0Dnvve9+rT3ziEyyHoS5wlwI1xhij3/72t3riiSe83mGtra168MEH1dXVxawFdYFwAWpMJpPRtm3bND4+Lun0cthtt92m97///QQL6gbhAtQQ27b1rW99S//3f//nXbvssst077330jsMdYVwAWqEMUY//elP9a1vfatiOeyhhx7S6tWrmbWgrhAuQA0wxujQoUN66KGHNDk5Ken0ctjtt9+u6667jmBB3SFcgCozxiiTyWjz5s168803vevXXHONNm7cqFCIt5Gj/hAuQJXl83lt375dL730kndtzZo12r59uzo7O6s4MuD8ES5AFeXzeX3961/Xf/7nf3r7LPF4XNu2bdMVV1zBchjqFvNtoAqMMcrlcvrqV7+qp556SqVSSdKf39HywQ9+kGBBXSNcgBVmjNH4+Li2bdum5557TrZtSzrd7fj222/XXXfdxT4L6h53MLCCjDHat2+fNm/erF/+8pdeQ8pgMKjbbrtNDz74oFpaWqo8SmDpCBdgBRhjlM1m9b3vfU///u//rpMnT3rfi0Qi+ud//mdt2rRJ8Xic5TA0BMIFWEbGGDmOo9/+9rf613/9V7388sve/ooktbe36/7779cdd9yhaDRaxZEC/iJcgGVijNHIyIiefPJJPfvss17rfNfQ0JAeeeQR/d3f/R2tXdBwCBfAZ+4S2AsvvKCdO3fq4MGDFd9vaWnRLbfcok2bNumCCy5gGQwNiXABfOQ4jl555RU9+uij+tnPflaxBGZZli699FI98MADGh4eVjgcJljQsAgXwAfGGE1NTek73/mOnnzySU1MTFR8v7u7W3fccYfuuOMO9fT0ECpoeIQLsETGGP3hD3/Qgw8+qJ/85CfeSXtJikajuummm/Qv//Iv2rBhA2+RRNMgXIAlsG1bP/7xj/W5z31Ohw4d8q5blqUNGzbogQce0A033KBIJMJsBU2FcAHOU6FQ0K5du/TlL39ZmUzGux6Px/WRj3xEn/70p7Vq1SpCBU2JcAEWya0G+/rXv66nnnpKMzMz3vf6+/v1xS9+UTfddJOCwSDBgqZFuACLYIxROp3Www8/rO9///teXzDLsnT11Vfra1/7mi655BJCBU2PcAEWyBij0dFR3X///XrxxRcr+oLdcsst+tKXvkQlGPAnhAuwAMYYHT16VJ/+9Kf18ssve8Hi9gV74IEH1NraSrAAf0K4AOdgjNFbb72lu+++W7/5zW+867FYTPfff78++clPUg0GzEK4AGfhtsi/++679bvf/c673tbWpq1bt+q2227j3SvAPHgqgDMwxujnP/+5Nm7cqLffftu73tXVpUcffVQf+tCHaDgJnAHhAszDcRz97//+r+67776Kd6/09fXpscce0/XXX89pe+AseDqAWQqFgp555hndc889FcFy0UUX6dvf/rY+8IEPECzAOTBzAf7EGKPJyUn927/9m7773e9WHI68/PLL9fjjj+uKK65g4x5YAMIF0J8rwr7whS/oRz/6UUXzyb/6q7/Sjh07NDAwQLAAC0S4oKkZY1QsFrVnzx5t3769ovlkMBjUBz/4QT3yyCP0CAMWiXBB0zLG6MSJE/rKV76i//qv/1I+n/e+F4/Hddddd+nee+/lcCRwHggXNKVSqaSXX35ZW7du1RtvvFHxvXXr1mnr1q36+7//e86wAOeJJwdNxRijVCqlxx9/XN/97nc1PT3tfS8YDOr666/X1q1bNTQ0xGwFWALCBU3DGKNXX31VDz30kH7xi194/cGk068h/sxnPqOPfvSjisfjBAuwRIQLmkKpVNILL7ygLVu26MSJE971QCCgq666Stu2bdOVV17J+RXAJ4QLGpoxRrlcTt/4xjf0zW9+U7lczvteW1ubPvGJT+hTn/qUOjo6mK0APiJc0LCMMcpkMtqyZYu+973veS/2kqShoSFt375d1157Lf3BgGVAuKAhuRv3DzzwgF544QVvfyUQCOiGG27Ql7/8ZfX39zNbAZYJ4YKG47ZxmR0s4XBYd955pzZt2qS2tjaCBVhGhAsaTj6f17Zt2/Tf//3fXrC0tLTovvvu01133aVoNFrlEQKNj3BBQymVSnriiSf07LPPev3BWlpatHnzZn384x9XOByu8giB5kDdJRqGMUY/+MEPtHPnTpVKJUmnl8I++9nPEizACiNc0BCMMdq/f78efvhhZbNZSac37z/2sY/p7rvvJliAFUa4oCFMTk7qC1/4gt555x3v2nXXXafPf/7zamlpqeLIgOZEuKDu2batp556Sj/96U+9awMDA3rkkUfU3t5exZEBzYtwQV0zxujnP/+5nnzySW8DPx6Pa9u2bVq/fj3lxkCVEC6oW8YYjY6OasuWLcpkMpIky7J0++2368YbbyRYgCoiXFC3SqWSHnvsMb366qvetfe9733auHEj72EBqoxwQV0yxuiHP/yhnn32We+gZGdnp7Zt26aurq4qjw4A4YK6Y4zRoUOH9MUvftErOw4Gg7rnnnv0l3/5lyyHATWAcEHdyWaz+tKXvqS33nrLu3bdddfpzjvv5H0sQI3gSURdsW1bTz/9tPbs2eNdW7t2rbZs2aLW1tYqjgxAOcIFdcMYox//+Md67LHHvPYu0WhUmzdv1qWXXspyGFBDCBfUBWOM9u3bp02bNlWUHX/4wx/WP/zDPxAsQI0hXFDzjDE6cOCA7rnnHh0+fNi7fvXVV2vTpk30DQNqEOGCmuY4jn7961/rjjvu0BtvvOFdHxwc1Ne+9jV1d3czawFqECfNUJOKxaLefvttPf/883r66ac1Pj7ufa+vr087duzQJZdcQrAANYpwQc2xbVsPP/ywnnvuOaXT6Yrv9fX16fHHH9c111xDsAA1jGUx1JxgMKiurq45wXL55Zfr6aef1rXXXkuwADWOcEFNGh4e9s6tdHd365Of/KS+//3vcwIfqBMsi6EmXXzxxfqnf/onXXTRRbrhhhvU39+vQCBAsAB1gnBBTXEcR7/61a+UzWb1N3/zNwoEAtq/f7/2799f7aHVhVKppImJiWoPAyBcUDsCgYAuvfRSvfTSS3rppZeqPZy6lUgk1NHRUe1hoMlZxu1XDlSZMUalUknckksXCoVo4omqIlwAAL7jow0AwHeECwDAd4QLAMB3hAsAwHeEC7BAVLMBC0e4AAu0b98+9fX1ad++fdUeClDzCBcAgO8IFwCA7wgXAIDvCBcAgO8IFwCA7wgXAIDvCBcAgO8IFwCA7wgXAIDvCBcAgO8IFwCA7wgXAIDvCBcAgO8IFwCA7wgXYAGMMZqYmJDjOJqYmOCdLsA5EC7AWaRSKe3cuVNDQ0O69tprNTExoWuvvVZDQ0PauXOnUqlUtYcI1CTL8BEMmNfevXt16623KpvNSlLFbMWyLElSPB7X7t27NTw8XJUxArWKcAHmsXfvXt18880yxshxnDP+fYFAQJZlac+ePQQMUIZwAWZJpVJat26dcrncWYPFFQgEFIvFdOTIESWTyeUfIFAH2HMBZnnmmWeUzWYXFCyS5DiOstmsdu3atcwjA+oHMxegjDFGQ0NDOnjw4KIqwizL0uDgoA4cOODtxwDNjHAByoyNjWnVqlVL+vnu7m4fRwTUJ5bFgDJTU1NL+vnJyUmfRgLUt1C1BwBUkzFGxWJRhUJBhUJB09PTS/rnjYyMqKWlRclkUi0tLT6NEqg/hAuaijHGC5JCoaBisShjjAKBgMLhsC688EINDg7q4MGDi/5nv+td75JlWfrd734nSYpGo0omk0omk+ro6FBra6vf/3eAmkW4oKE5jjMnTKTT5cORSESJREKRSEThcNj7mY0bN+qzn/3sojf077vvPl111VUqFotKp9NKp9NKpVI6efKkjDEKh8Ne0CSTSbW1tbH5j4bFhj4aim3bFWFSKpUkScFgUJFIxPsrFDrz5yq/z7nYtq1MJqNUKqVUKqVMJiPHcRQMBr2gSSaTSiQSCgTYBkVjIFxQ10qlUkWY2LYtSQqFQhVhEgwGF/XPXewJ/RdffFE33njjgv7ZjuNocnJSqVTKm93Ytq1AIKD29nZvdtPR0bHocQO1gnBBXSnffC8UCt4v/nA4XBEmfswAFtpb7Pnnn19wsMzHGKPp6WlvZpNKpVQsFmVZlhKJhDe76ejoqFi+A2oZ4YKaNbuSq1AoyBgjy7LmhMly7V2kUint2rVL3/jGN/TWW29519evX6+NGzfqYx/7mDo6Onz/92azWS9o0um08vm8JKm1tbVi3yYajfr+7wb8QLigZpypksuyrIogCYfDK74RbozRqVOnNDk5qUQioa6urhUdQz6f95bQUqmUN5uKxWIV+zaxWGzFxgScDeGCqjlXJVd5mKBSoVCoCBv38GckEvGCJplMKh6PU5GGqiBcsGLcMJmZmVGxWPTCZDGVXJhfqVSqKH/OZDIyxigUClUsoyUSCcIGK4JwwbI5U1nwUiu5cG62bXsVaW7Y2LatYDDoVaQlk0m1t7dT/oxlQbjAN2cqC16OSi4sjjFmTvlzqVSSZVlzyp+ZOcIPhAvO27nKgqPRqMLhMGFSg4wxFRVpqVRKhUJBlmV5FWlu4EQikWoPF3WIcMGC1EJZMJZXLperKH/O5XKSTp/lKQ8bGnJiIQgXzKuWy4KxMmZmZiqW0dyO0S0tLRXlz/F4vMojRS0iXCCJsmCcm9uQ0w2cycnJioac7syGhpyQCJem5UeDRzQ327YrwsZtyBkKhbziABpyNi/CpUksV4NHwOU4jjKZTEXgzG7I6ZY/c581PsKlQa1kg0dgPsYYTU1NVRQJlDfkLC9/Zrm18RAuDYBKLtQDt/y5vG3NzMyMJM0pf6YhZ/0jXOqQGyYzMzNUcqGu5fP5irM2bvlzLBaraFtDQ876Q7jUCTdI3L5cEpVcaDyFQqGi/NltyBmNRr2g6enpYWZTBwiXOjE+Pi7HcbwQoZILzcC2bU1NTXl/5XI5DQwMLMs7dOAvwqVOOI7D5juanm3bsiyLZ6EOEC4AAN8R/wAA37Fo/yduBZZ7HgTnLxwOc0iuThljvMOPWJpEItHUHaUJlzLpdFqSKN9dglKppGQySbjUKcdxdODAATmOw3/DJcjlctqwYYN6enqqPZSqIVzKuCeHm/nTxmIYYzQzM+Od9DfG6NSpU9UeFpbIsiy9+93vVmdnZ7WHUheMMRobG1Nra6vi8bgcx9Grr76qZt/OJlzmwczl3NyW/KlUSsFgULFYjINuDYSKrIVx3/D55ptvKhAI6KKLLlJvb2+1h1UTCBecN/c0tXsWgUOcaDa2beuPf/yj11X8wIEDVR5R7eCjCc5LqVTy+kJJf26ICTQLY4yOHz+uiYkJ71pra6u6u7urOKraQbhg0YwxyufzFWvKsViM5UQ0DbcJ5+HDh71rgUBAg4ODfMj6E8IFi+aGiysQCNDrCU3FGKNDhw6pUCh41/r6+iiCKEO4YFHcjXx3jVk63VSQslU0C7cqcnR01LvW0tKi/v5+Zu9lCBcsmruR76JKDM3Etm0dOnTIO3BtWZYuvPBCtbS0EC5lCBcsim3bFUsB4XCY98agaRhjdOLECWUyGe9aR0eHent7eQZmIVywYO6hyfIWOXxaQ7Nwl4Tfeecd75p7toVl4bkIFyxK+Ua+ZVmEC5rKsWPHKpaFe3t71dHRwTMwD8IFC1Yqlby3YEpSJBLhExuagjFGuVxOx44d866Fw2G9613vIljOgHDBgsx3tqWlpaWKIwJW1tGjR+eUHre2thIuZ0C4YEHOdLaFBwuNzj0weeLECe9aNBrVBRdcUMVR1T7CBedkjFGpVKo42+J2QgaawdGjRyuWhNesWcN+4znw2wELUj5rkVgSQ3NwZy0nT570rrW0tGjNmjUEyzkQLjgntwTZFQwGFYlEeLjQFI4dOzZn1kK7o3MjXHBW7uufWRJDs3ErxGbvtfT19fHBagH4DYFzKp+1SCyJoXkcP36cWct5IlxwVvMtidHuBc1gZmamYtYSiUSYtSwC4YIzYkkMzcrtIVb+wWr16tXM2heB3xI4K5bE0IyKxaJGRka8r0OhkNasWVPFEdUfwgVnxJIYmpExRuPj48pms961VatWKR6Pc+8vAuGCebEkhmblOE5FD7FgMMis5TzwmwJnxJIYmo0xRqlUSpOTk961ZDKpRCLBrGWRCBfMiyUxNCNjjI4fP+41aLUsS2vXruW+Pw+EC+ZFLzE0o+npaZ06dcr7uq2tTclkknA5D/y2wByzZy2SODiGhueWH9u27V1bs2YN7yw6T4QL5pgdLoFAgF5iaHiFQkGjo6Pe19FoVD09Pdz354lwwRzzvXGSJTE0Mrf8uLz7d29vryKRSBVHVd/4jYEK8y2JUSWGRuc4TsWhyWAwqNWrV1dxRPWPcMEcLImhmRhjlMlk5pQf8wrjpSFcUIEqMTSjEydOyHEcSafLj2lQuXT81oDHXRJza/wlqsTQ+GZmZjQ+Pu59HYvF1NnZSbgsEeGCCuUbmiyJodEZYzQ2NqZCoeBd6+3tVSgUquKoGgN/gmfgfnp3HEfGGFmW5S0PNeov29lLYuFwmBp/NDTHcSre2RIKhbRq1aqGfcZXEuEyj1KppHw+r3w+7/2ytSxLoVBIsVhMLS0tsiyroW7A+ZbEqBJDIzPGaHJyUlNTU961ZDKpeDxexVE1DsJlHpOTkxXLQ9LpG7FQKKhQKCiXyymRSDRcry2qxNBsyjfypdMvBOOe9wd7LvM41yZ2oVDQxMSE8vl8xSf9ejb74CRLYmh0hUKBjfxlxMxlHtFoVMFgUI7jKBKJKBgMqlQqVWz6OY6jdDotY4xisVhd35AsiaHZuCfyy2frq1atYiPfR/xJziMQCKijo0PBYLDi0/vMzIwmJye9fRj38JVlWd4+TL2iSgzNZPZGfjAYVG9vL/e8jwiXeViWNe/SWDQaVSgUUjqd9mYxxhil0+m6/oVMlRiaiTFG09PTFSfyOzo61NraWsVRNR72XBbBsiwFg0Elk8mKhnZuwNi2XXd7MCyJoRmdPHmyorU+G/n+I1wWyT3vkkwmFQ6Hveu2bSuTydRduEgsiaG5FItFjY2NeV+3tLSwkb8MCJfz4AZMR0dHRd+tmZkZTU9P11XAFIvFOe31WRJDozLGaGJiQrlczrvW09NDa/1lQLicJ/dQZXt7e8UnnunpaRUKhboIGGPMnPM8LImhkblvm3QFAgH19vZWcUSNi3BZArdKrPxEr1tB5raNqWW8cRLNJpvNKpVKeV8nEgm1tbVxzy8DwmWJLMtSW1tbxbS6VCpVtJSoRcYYFYvFiiqxaDRKe300LGPMvBv53PPLgz9VH1iWpUQiUXGTZrPZOVVYtaZ83VliSQyNrVQqaXR01Ps6Eomou7ubWcsyIVx8YFmWwuHwnDr5ycnJml0ecxynYkksFAqxJIaGZYxRKpWq+EDV3d3N+4qWEeHiE8uyFI/H62J5zN1rKW/YF41GCRY0LGOMRkZGvA96gUBAq1evrvKoGhvh4iPLstTe3l6xPJbL5Wpyeay8SswtTAAakTFmzkZ+W1vbnEpP+Itw8ZFbnly+POa+M6KWlsdmN+EMh8MN9/oAoNzJkycrilf6+vrYyF9m/On67EzLY9PT01Uc1Z+5Z1vKgy4Wi1VxRMDyKhaLOnnypPd1NBplI38FEC7LwK0eK795s9lsTRyuNMZUbGoGAgH2W9CwjDE6derUnBP5bOQvP8JlGcxXPeYuj1UzXNyN/PI6/5aWFpYH0LAcx9Hx48e9r9nIXzn8VlkmlmWptbW1orllsViseu+x2WdbWBJDo3K7ZWQyGe9aMpmcs6qA5UG4LKP5lsemp6dVLBZXPGDcE/nlG/mRSISNfDQsY4yOHz/uldxblqU1a9Zwv68QwmUZWZalSCQyp/dYtZbHcrkcG/loCm758fj4uHctHo/TWn8F8SbKFdDa2qqZmRmvFLJQKCibzaq1tXXFbnTbtivOtgSDQTbycd6MMTLGqFQqVdzb4XDYmxFLqur9NTIyMqf8OBTiV95K4U96mbnvfkkkEkqlUt7MYWpqasWWpdwKsfIT+bFYjI18LJpbFDI+Pq7R0dE5y7yBQEDhcFjt7e1avXq1Ojs7FQwGVzxk8vl8RWv9aDSq3t5ePkytIMJlBViWpWg0qng87p13cTcbu7q6lv2GdxxnTvlxLBbjQcOCGWNUKBR09OhRjYyMVPSlK+f2rBsdHdXY2Jja29t10UUXKZlMyrKsFbnn3FYv5fuLfX19lB+vMD66rhC3eqx8Wl4sFjU1NbWs+y/urGV2+TFvm8RCGGPkOI5GR0f1yiuv6NChQ2cMlvl+Np1O67XXXtMf//hH2ba9InuNMzMzFeXH4XBYfX19fJhaYcxcVlAgEFB7e7smJia8hyybzSoSiSzb/ofjOMpms97XbgcB4FzcCsO33367ourKFQwGlUgk1N7e7s0Kcrmc0ul0xYcm27Z1+PBhTU9P6+KLL17WvT5jjI4dO1YRgKtXr6Z4pQoIlxXkVo+1trZ63ZLLl8f8Xpt2K2bKZy2xWEyhUIhPcTinbDar/fv3K51OV1wPhULq6+vTmjVrFI/HK5a73JlOKpXS22+/XXHGZHx8XK+//ro2bNiwLMuy7ix99qxl7dq1vv57sDAsi60wd3msvPeYbdu+lycbY2TbNrMWnBfbtucNlu7ubl155ZV697vfrdbWVgUCgYqQsCxLwWBQXV1deu9736v+/v6KwpFMJqPXX399WQ4TG2P0zjvvzNlrcQMQK4twqYL5WvPn83nfH7jp6ek5FWLMWrAQwWBQ69at8+7RYDCowcFBXX755d475892H7ktkAYGBnTJJZdU7DVOTU3pjTfe8PV+d1cAZleIXXDBBdzvVUK4VIHbmj+RSFRcn56entOx+Hy4lT2zK8RmvykTOJuenh719/crFovpsssuU39//6KXbt1eXhs2bKiYrU9PT/saMLZte0UDrnXr1vGeoioiXKrEsizFYrE5zS0zmcySuycbY+ZUobW2tlblvAHqVyAQUH9/v6688soltai3LEvd3d1nDJhsNrvk+/3o0aMVLwNLJBK0eqkywqWKLMtSW1tbRf294zhKp9MqlUrn9cAZYzQ9PT3nZWCsO+N8BINBtbS0LPnesSxLnZ2duuyyy+YEzOuvv37eAWOMUSqV0uHDh71rgUBAAwMDnMavMsKlytz9l/IHwbZtpVKpRQeMe3q6/MVk8zXPBKrBsiwlk8l5ZzDns8nvVoe9+eabFW1e1q5dSw+xGkC4VJlbXZNMJisONpZKJU1MTCy4g7Lb5ymTyVT8/e5bMXnQUAvcGcx8AfPaa68pnU4v+H6fmZnR73//+4qKyPb2dl144YXc7zWAcKkB7gZ/R0dHRQWZbduamJjwNvnP9NC5wZJKpSo2NCORiFfZA9SKMwVMLpfT66+/rhMnTshxnLPe77lcTm+88UZFqXQkEtHQ0BCvkagRhEuNcA9YJpPJioBxD6Sl02mvfYb70Ln/O5/Pa2JiomJpIBgMqr29nYcMNal8D6Z8z7FQKGj//v168803vVdElN/vtm1rdHRUr776akWwhEIhXXzxxSwB1xB2vGqIGzCdnZ3epr4rl8spn88rGo0qEokoEAjItm3NzMxUbN5Lpzc0k8kkZ1pQ09w9mPe85z36/e9/7+0Vuq8mHhsbU3d3tzo6OhQMBpXL5XTq1Kk5S7+hUEhDQ0Pq6enhfq8hhEuNcQ+fdXZ2KpPJVPRIcmcp5e9lmS0YDKqjo4OlAdQFt+Dkiiuu0B/+8AeNjY153ysWixoZGdHIyMgZfz4ajRIsNYpwqUHlm/y5XE5TU1NzmgbOJxqNqr29nfMsqCuWZamlpUUbNmzQ8ePHdfjw4Tmz8fl+prOzU+vXr1/Rl+5h4QiXGuW214jH44pGo96y2Oy25e5MJx6Pe6eRedBQb9yilnXr1qmnp0fHjx/X6Oio8vl8xQerUCik9vZ2rV27Vl1dXXN6m6F2EC41zn3o2tra1NraKtu2vYAJBAIKBoNeCTMPGeqdO4sZGBhQf3+/crmcZmZm5DiOwuGwYrGYV1rP/V7bCJcybk+ulXih0VJZluWVIJdv/FfT2cqlUT/cNkQLWYpdKW4Bi/uqilrmPpfNjnApEw6HNTMzs+A37WEuPlHWv9bWVo2Pj2t8fLzaQ6lbwWBQ4XC42sOoKsvwUVOS+MTtMwKmPjH79E+zf9AiXAAAvuOEPgDAd4QLAMB3hAsAwHeECwDAd4RLnTDGLPjdLkCjsm1b2Wy24tUSqE2ES50olUoaGxvjcBaaWjab1S9/+cuKF4ShNhEuAADfES4AAN8RLgAA3xEuAADfES4AAN8RLgAA3xEuAADfES4AAN8RLgAA3xEuAADfES4AAN8RLgAA3xEuAADfES4AAN8RLnXAGKOxsTG98847Ghsb450uaErGGI2Pj2tkZETj4+M8BzWOcKlhqVRKO3fu1NDQkNauXaurrrpKa9eu1dDQkHbu3KlUKlXtIQLLrvw5GBgY0Ic//GENDAzwHNQ4yxD/NWnv3r269dZbvZcilf9nsixLkhSPx7V7924NDw9XZYzAcuM5qF+ESw3au3evbr75Zhlj5DjOGf++QCAgy7K0Z88eHiw0HJ6D+ka41JhUKqV169Ypl8ud9YFyBQIBxWIxHTlyRMlkcvkHCKwAnoP6x55LjXnmmWeUzWYX9EBJkuM4ymaz2rVr1zKPDFg5PAf1j5lLDTHGaGhoSAcPHlxUJYxlWRocHNSBAwe8dWigXvEcNAbCpYaMjY1p1apVS/r57u5uH0cErDyeg8bAslgNmZqaWtLPT05O+jQSoHp4DhoD4VJD2tralvTziUTCp5EA1cNz0BgIlxrS3d2t9evXL3q92LIsrV+/Xl1dXcs0MmDl8Bw0BsKlhliWpXvvvfe8fnbjxo1sYqIh8Bw0Bjb0awz1/QDPQSNg5lJjksmkdu/eLcuyFAic/T+PezL5+eef54FCQ+E5qH+ESw0aHh7Wnj17FIvFZFnWnGm+ey0Wi+nFF1/UjTfeWKWRAsuH56C+ES41anh4WEeOHNGOHTs0ODhY8b3BwUHt2LFDR48e5YFCQ+M5qF/sudQBY4xOnTqlyclJJRIJdXV1sWmJpsNzUF8IFwCA71gWAwD4jnABAPiOcAEA+I5wAQD4jnABAPiOcAEA+I5wAQD4jnABAPiOcAEA+I5wAQD4jnABAPiOcAEA+I5wAQD4jnABAPju/wFmaj6x3Z6ioAAAAABJRU5ErkJggg==\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": "code",
   "execution_count": 9,
   "id": "f2dc6ceb",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "train loss: 3.72e-03 | test loss: 3.45e-03 | reg: 3.42e+00 : 100%|██| 20/20 [00:03<00:00,  5.13it/s]\n"
     ]
    }
   ],
   "source": [
    "model.train(dataset, opt=\"LBFGS\", steps=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2a302028",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAFICAYAAACcDrP3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAqyElEQVR4nO3de4xc1X0H8O+Z587szuzszr5sr8H7sknsgBMrpWrU0qRgI5y0URxFgkaqWtGmCcUFxQkFlKCShBZDUpsoSpo0MbZAiaoYWopBVsIjqhSiCBEwj9RevH6t7d31zM5zZ3Ze9/QPcy/3zK7tfdzZe2fm+5EseWa99vHaZ7/3vH5HSCkliIiILOSyuwFERNR4GC5ERGQ5hgsREVmO4UJERJZjuBARkeUYLkREZDmGCxERWY7hQkRElmO4EBGR5RguRERkOYYLERFZjuFCRESWY7gQEZHlGC5ERGQ5hgsREVnOY3cDiOqBlBLxeBzZbBZtbW2IRqMQQtjdLCLH4siF6DKSyST27t2LkZERdHd3Y2BgAN3d3RgZGcHevXuRTCbtbiKRIwneREk0v8OHD2PHjh3I5XIALo5edPqoJRgM4uDBg9i2bZstbSRyKoYL0TwOHz6M7du3Q0oJTdMu+etcLheEEDh06BADhsiE4UJUJZlMor+/H/l8/rLBonO5XAgEAhgfH0ckEql9A4nqANdciKrs378fuVxuQcECAJqmIZfL4cCBAzVuGVH94MiFyERKiZGREYyNjWExXUMIgcHBQYyOjnIXGREYLkSKWCyG7u7uZX1+NBq1sEVE9YnTYkQm2Wx2WZ+fyWQsaglRfWO4EJm0tbUt6/NDoZBFLSGqbwwXIpNoNIqhoaFFr5sIITA0NITOzs4atYyovjBciEyEELjzzjuX9Lk7d+7kYj7Re7igT1SF51yIlo8jF6IqkUgEBw8ehBACLtflu4h+Qv+pp55isBCZMFyI5rFt2zYcOnQIgUAAQog50136e4FAAM899xy2bt1qU0uJnInhQnQJ27Ztw/j4OPbs2YPBwUHlY4ODg9izZw/Onj3LYCGaB9dciBZASomXXnoJf/Znf4YXXngBH//4x7l4T3QZHLkQLYAQwlhTiUQiDBaiK2C4EBGR5RguRERkOYYLERFZjuFCRESWY7gQEZHlGC5ERGQ5hgsREVmO4UJERJZjuBARkeUYLkREZDmGCxERWY7hQkRElmO4EBGR5Vhyn2iBpJSQUs57eRgRqThyIVoEhgrRwnjsbgBRvWCwEC0cRy5ERGQ5jlzIMSqVCl555RWk02m7m1L3tmzZgt7eXrubQU2MC/rkGPl8Hn/8x3+MfD6PUChkd3PqkpQSo6OjePzxx/Hnf/7ndjeHmhhHLuQYUkq4XC585zvfwZ/+6Z/a3Zy6cfr0abz00ku47bbbAAC33HIL+MxIdmO4kON4PB74/X67m+F4mqbhlVdewZe+9CWMjo7C6/Xitttu48YDcgSGC1GdkVKiWCziiSeewP3334+pqSkAwL333ovBwUGbW0d0EcOFqI5IKXH27Fk88MADePLJJ1EsFo2PCSEQj8dtbB3R+7gVmagOSClRKBTw85//HDfddBP27dunBMvmzZvx9NNPY/v27Ta2kuh9HLkQOZymaXjnnXfw0EMP4emnn0ahUDA+5vF4sGPHDuzevRv9/f0olUo2tpTofQwXIoeSUuL8+fP493//d/zgBz/AhQsXlI/39PTg/vvvx+23346WlhYu5JOjMFyIHEZKiVQqhSeffBJ79+7F8ePHla3FHo8HN910E771rW/h2muvhcvF2W1yHoYLkUNIKZHP5/Hss89i9+7deP3116FpmvJrhoaGcM899+DWW29FMBjkaIUci+FCZDMpJcrlMv73f/8X//qv/4pf/epXc9ZOOjs78Td/8zfYuXMn1qxZw1Ahx2O4ENlI0zS8/fbb2L17N55++mnkcjnl462trfj0pz+NXbt2YdOmTXC73Ta1lGhxGC5ENtAX67/3ve/hRz/6EWKxmPJxn8+HT3ziE7jnnnvwR3/0R/B4PBytUF1huBCtICklcrkc/vM//xMPP/wwRkdHlcV6l8uFLVu24Ktf/SpuueUW7gKjusVwIVoBUkqjFtiDDz6Il19+GeVyWfk1Q0NDuPvuu3Hbbbehvb2doUJ1jeFCVGN6yZZHH30Ujz/++Jz7aqLRKG6//Xb8wz/8A1avXs1QoYbAcCGqESklZmdn8fOf/xzf/OY3MTo6qnzc7/fjU5/6FO677z6eV6GGw3AhqgF9F9gDDzyAQ4cOKVuLhRDYvHkzvv71r+Pmm2+Gz+fjaIUaDsOFyEJSSmSzWfzoRz/CI488gsnJSeXjPT09uPPOO/H3f//36OzsZKhQw2K4EFlE0zS89tpruPfee/Hyyy+jUqkYH/N6vfjkJz+JBx54AJs2beIUGDU8hgvRMkkpkclk8P3vfx+PPvronDtVRkZG8LWvfQ2f/exn4ff7OVqhpsBwIVoGTdNw5MgR3HPPPXjhhReUWmCBQACf//zncf/992Pt2rUMFWoqDBeiJdB3gu3btw/f+MY35qytbNq0Cd/61rdw880383Q9NSWGC9EiSSlx4sQJ/NM//RP+67/+SzkMGQwG8dd//de477770NfXx1ChpsVwIVqEcrmMQ4cO4atf/eqccyvXXHMN/uVf/gW33HILRyvU9BguRAugX+C1e/dufPe738XMzIzxMZ/Ph1tvvRUPPvgg+vv7GSpEYLgQXZGUEr///e9x11134cUXX1QW7VevXo1vfOMbuO2223gYksiE4UJ0GeVyGf/93/+NXbt24dSpU8b7LpcLn/jEJ/Dtb38bmzZtYqgQVWG4EM1DP2n/yCOP4N/+7d+UabDW1lb84z/+I77yla8gHA4zWIjmwXAhqiKlxOnTp3H33XfjmWeeUabBBgYG8J3vfAfbt2+Hx8PuQ3Qp7B1EJlJK/PrXv8aXvvQlvPnmm8b7LpcLN910E/bs2YP169dztEJ0BSxwRPSecrmMJ554Ajt27FCCpaWlBV/+8pfxs5/9jMFCtEAcuVDTk1Iin89j9+7dePTRR5HL5YyP9fb24uGHH8att97KsytEi8BwoaYmpUQ8HseuXbvw5JNPKpWMr7vuOnz/+9/H9ddfz1AhWiSGCzUtKSVOnjyJL3zhC3jhhRcgpQRwcX1l+/bteOyxx3DVVVcxWIiWgGsu1JSklHj99dfxmc98Br/85S+NYPH5fLjzzjtx4MABBgvRMnDkQk1HSomXXnoJt99+O06ePGm8HwqF8M///M/44he/yNP2RMvEcKGmomkannnmGXzxi19UyuT39vbisccew2c+8xm43W4bW0jUGBgu1DQqlQp+9rOfYefOnUgkEsb7Q0ND+I//+A/8yZ/8CUcrRBbhmgs1hUqlgv379+OOO+5QgmXz5s04ePAgg4XIYhy5UMPTg+Xuu+9GJpMx3v/Yxz6Gxx9/HIODgwwWIotx5EINTdM0PPnkk3OC5cYbb8RPf/pTBgtRjTBcqGFpmoann34ad911lxIs27Ztw/79+7FmzRoGC1GNcFqMGpKUEi+++CLuuOMOJJNJ4/2tW7di37596O3tZbAQ1RBHLtRw9AOSf/d3f4epqSnj/RtuuAE//vGPGSxEK4DhQg1FSokzZ87gb//2b5UDklu2bMFPfvITrF69msFCtAIYLtQwpJRIp9PYuXMnXnvtNeP94eFh7Nu3D+vWrWOwEK0Qhgs1jHK5jIceegjPPvus8V5PTw9++MMfYuPGjQwWohXEcKGGIKXET3/6U3z3u981riUOBoN45JFHeECSyAYMF6p7Ukq8+uqruPfeezE7OwsAcLvd+PKXv4xbb70VLhf/mxOtNPY6qmtSSsRiMdx11104f/688f5f/MVf4Ctf+Qo8Hu62J7IDw4XqWqlUwoMPPojf/OY3xnubNm3Ct7/9bbS2ttrYMqLmxnChuiWlxMGDB/HjH//YuOwrEolgz549vOiLyGYMF6pLUkocO3YM999/v7LOcs899+CGG25gsBDZjOFCdSmXy+G+++5TDkp+8pOfxB133MHLvogcgOFCdUfTNOzbtw//8z//Y7w3MDCAhx9+mOssRA7BcKG6IqXEW2+9hYceegjlchkA4Pf78c1vfhMjIyOcDiNyCIYL1ZXZ2Vk88MADmJiYMN77/Oc/jx07djBYiByE4UJ1Q0qJZ555Bs8//7zx3oYNG/C1r30NXq/XxpYRUTWGC9WNZDKJRx55BMViEQDg8/nw9a9/HWvXruWohchhGC5UF/QzLW+88Ybx3tatW/HpT3+awULkQAwXqguZTAY//OEPUalUAACtra3YtWsXWlpabG4ZEc2H4UKOJ6XEr371K2XUcvPNN+MP//APOWohciiGCzlepVLBE088gVKpBODi1uPbb7+di/hEDsZwIcc7deoUXn75ZeP1ddddh4997GMctRA5GMOFHE1KiV/84heIxWLGe5/97Gd5Ep/I4Rgu5GjlchnPPvusUvX4lltu4aiFyOEYLuRo586dw6uvvmq8/vCHP4yhoSEbW0REC8FwIceSUuK3v/2tMiW2detW+Hw+G1tFRAvBcCHHklLixRdfhKZpAIBAIICPf/zjnBIjqgMMF3KsXC6nXF+8bt06bNiwwcYWEdFCMVzIsU6dOoUTJ04Yr//gD/4AoVDIxhYR0UIxXMiRpJR4/fXXkclkAABCCJ5tIaojDBdyrN/+9rfGFuRAIIAtW7YwXIjqBMOFHKlYLOK1114zXq9atQrr1q2zr0FEtCgMF3KkWCyG48ePG683btyIcDhsY4uIaDEYLuRIY2NjmJ6eNl5v3rwZLhf/uxLVC/ZWcqQjR44YN066XC58+MMftrlFRLQYDBdypBMnTsDv9wMAQqEQrrnmGi7mE9URj90NIJrPfffdh8997nM4cuQIzp49izVr1tjdJCJaBIYLOYqmaXjllVcwMzMDAOjp6UFPTw9efPFFm1tWH8rlsrJWRWQXhgs5hsvlwsaNG/H888/j+eeft7s5dSscDqOjo8PuZlCTE1I/pUZkMyklyuUy+F9y+TweD3fXka0YLkREZDk+2hARkeUYLkREZDmGCxERWY7hQkREluNWZKIFMu99YbUAosvjyIVogX73u9/B5XLhd7/7nd1NIXI8hgsREVmO4UJERJZjuBARkeUYLkREZDmGCxERWY7hQkRElmO4EBGR5RguRERkOYYLERFZjuFCRESWY7gQEZHlGC5ERGQ5hgsREVmO4UJERJZjuBAtgJQSiUQCAJBIJJS7XYhoLoYL0WUkk0ns3bsXIyMjuPHGGwEAN954I0ZGRrB3714kk0l7G0jkUELyEYxoXocPH8aOHTuQy+UAzH8TZTAYxMGDB7Ft2zZb2kjkVAwXonkcPnwY27dvh5QSmqZd8te5XC4IIXDo0CEGDJEJw4WoSjKZRH9/P/L5/GWDRedyuRAIBDA+Po5IJFL7BhLVAa65EFXZv38/crncgoIFADRNQy6Xw4EDB2rcMqL6wZELkYmUEiMjIxgbG1vUjjAhBAYHBzE6OmqsxxA1M4YLkUksFkN3d/eyPj8ajVrYIqL6xGkxIpN0Or2sz89kMha1hKi+eexuAJGdKpUKEokE4vE4YrEYTp48uazfb2JiAi0tLYhEImhpabGmkUR1iOFCTaVUKmF6ehrxeBzxeByJRAKapsHr9SIajeKjH/0o1q1bt6SQWbt2LYQQeOeddwAAfr8fkUjE+NHa2mrx34bIuRgu1NAKhYIRJPF4HKlUClJKtLS0IBqNYtOmTejq6kIoFDIW4u+66y7cfffdi17Q37VrF66//nqUSiWkUikkk0kkk0lMTU1BSgmv16uETVtbGxf/qWFxQZ8aSj6fRywWM8JEXwMJBoPo6upCNBpFNBpFW1vbJX8Pq8+5VCoVI2xSqRRSqRQ0TYPb7UZ7e7sRNuFwGC4Xl0GpMTBcqK5ls1klTPRSLaFQSAmTQCCwqN93sSf0n3vuOWzdunVBv7emachkMsbIJplMolKpwOVyIRwOG2HT3t4Ot9u9qHYTOQXDheqGlBLpdFoJk0KhACEE2tvbEY1GjUDx+XzL/vMWWlvsqaeeWnCwzEdKiWw2q4RNqVSCEAKhUEgJG6/Xu7y/FNEKYbiQY2mahmQyqYRJuVyGy+VCR0eHESSdnZ3weGqzfJhMJnHgwAE89thjOH78uPH+0NAQdu7cib/6q79Ce3u75X9uLpdTwmZ2dhYA0Nraqqzb+P1+y/9sIiswXMgxKpWKsZMrFoshkUigUqnA4/Ggs7PTGJlEIpEVny6SUmJ6ehqZTAahUAidnZ0ruhg/OzurhI0+mgoEAkrYLHb6j6hWGC5km1KpZIxIYrEYkskkpJTw+XzGWklXVxfa29u5q6pKsVhEKpVCIpFAKpUyNi74fL4525/5tSM7MFxoxczOziphop+Gb2lpMaa4urq6uEV3CcrlsrL9OZ1OQ0oJj8ejhI15yzVRLTFcqGZmZmaUMJmZmQEAtLW1KSOTYDBoc0sbT6VSQTqdVsKmUqnA7XbP2ZHG7c9UCwwXsoS+40lffI/FYsYidDgcVrYFsyzKypNSztn+XC6XIYSYEza12hxBzYXhQksipUQqlVJ2chWLRQgh0NHRYQRJNBrl9lkHklJiZmZGCRv936+trQ3t7e3o6OhAe3u7Jdu6qfkwXGhBKpWKsi14enoa5XIZbrd7zrZgHvyrT/l8XgmbfD4P4OJZHvO6DUeetBAMF5pXuVzG9PS0ESbmAo+dnZ1GmEQiEc7ZN6hCoaCEjb5mpld91n9wzYzmw3AhABe3tuprJeYCj36/X1l8D4fD3G3UpKoLcmYyGRbkpEtiuDSpfD6vhIm5wKO5jMrlCjxSczMX5NR3pGmaBo/HoxTkDIVCHN02IYZLk8hms0qYmAs8mkcmPOFNS6VpmrL9OZVKsSBnE2O4NCC9wKP5jEl1gUf9B2tTUa2wIGdzY7g0AL3Ao/lSrFKpZBR4NIcJzzCQXaSUcwpyFgoFADC2P7MgZ+NguNQhc4FHfVuwfvraHCQdHR2cfiBHY0HOxsVwqROTk5PKtmB9l4558Z2lPKjeFYtFJWyy2SwAwO/3G0HT1dXFkU0dYLjUibfeegvFYhGhUAjhcBihUAiBQIBbPqmhVSoVZLNZ40cul8Pg4GBN7tAhazFc6kS5XOZ6CTW9SqUCIQRH6HWA4UJERJZj/BMRkeU4z/Ie/RrbUqlkd1PqHosb1i+92nWlUrG7KXUvFAo1dUVphst7NE3DkSNHUKlUeKBrifQy7h/5yEewatUqu5tDS6BpGkZHR6FpGrexL0M+n8cHPvABdHV12d0U2zBc3iOlhBACH/rQh9DV1QUpJQqFAgqFAnemLJCUEr/+9a/tbgYtkxACw8PD6OjogJQSxWIRxWKRBSkvQUoJKaWxyUB/UG325WyGSxUhBMbHx3HmzBlkMhm43W7ccMMN3Fc/j2QyifPnzyMajSIcDsPr9fKbT4MQQmBychITExOYmZmB2+3Gli1bmnqaZz76dPqpU6eM+ml8GL2I4VJFCGFc16u/zmQyDJcqUkpMTk7i6NGjEELA7/fjox/9qN3NIgvppVqAi1vhc7kcw2UeiUQC6XQa6XQa586dw3XXXWd3kxyBu8XmEY1GjSdw/cmk2Ye41aSUiMfjxs+llFzEbzDmJ3B9oZ/9QKVXgtb5fD6WqnkPw2Ue+hSPLh6Ps1NVKRaLSqcKh8Mc3TWYtrY2pR8kk0n2gyqFQsGohwZc3CHGw84XMVzm4ff7EQqFjNfpdJpblKuk02mjoi0AdHZ28tR0g/H5fGhtbTVez8zMoFwu29giZ9GvFDB/TSKRCNcd38PvBvNwuVzo7Ow0Xs/Ozho3NdL7U2L6U6wQAtFo1OZWkdX0+390xWIRMzMzNrbIefQ1KQDGpWh0EcNlHkIIdHZ2ct3lEvSvh656pEeNQQ8Xcz/gusv7NE1DKpUyXvv9fq63mDBcLqH6djyuu7yvUChwvaVJcN3l0mZnZ5HP543X4XCY6y0mDJdLaGlpUYa4qVQKxWLRxhY5g/70av5amHfXUWPxer1z1l3YD96/StxcJicSidjXIAdiuFyCPjWmq35ab2bV6y3mKURqLEII5ZtmqVQyLvBqdub1FrfbjXA4zH5gwnC5BCEEurq6lPlmTo1dnGfWz7cAF6+j5SJm45pv3YVTYxcPlZofNgOBAM95VWG4XEb1WgLD5WJBPvPOufb2dp7abnCtra1KP0ilUtA0zcYW2S+fz2N2dtZ4HQ6HWeizCsPlMnw+n7IVM51OK/+hmo3+1Go+89PMVV+bhdfrVXYD5nI55YxTs9H7gTlgOzo6bGyRMzFcLkOfGtMVi8WmnxK4cOGC8XO3283F/CZhXnfRp4SatR/o4aLzeDwIhULsB1UYLpehHw40nzw3f3NtNuVyWTnfEgwGlZ1E1Jj0dRdzP0gkEja2yF6lUkmZGm5ra+NW/HkwXK6gra0NwWDQeD09Pd20t/Rls1mljlJnZyf39TeJQCCg9INUKtW0pWAymYwyNcySL/NjuFyB1+tV5lOz2WxTlsCQUiIWiynB2t3dbWOLaCW53W5l/bG6YGOzqN4tp2/VZrjMxXBZAPM30Uql0pS7xvRw0fl8PnR0dLBTNRHzQ5amaU25/qj/vXV+v59Tw5fAcLkC/ZCguQTGhQsXmq5TFQoFpVOFw2Hu628iQgiEQiGlHyQSiabrB/l8XhmxVV/PQe9juCxAIBBQtmImEommKoEhpZzzd+7u7maJ/Sbj8/mUfpDNZptqS7Je+sg8NWyu4kEqfndYAJfLpWxJ1p/im+mpzTxa078enBJrLkIIZWqsVCo11Zbk6mrgHo9HqV5AKobLAggh0N3drZTAmJqasrlVK6dcLivrLcFgkCX2m5C+eG0esZq/2TY6PUx11ZULSMVwWaD29nZlK2Y8Hm+aLcmZTEbZIReNRjnP3KSCwWBTbknWqyCbp4Y7Ojo4NXwZ/MoskNfrVeZXs9lsU9xOKaXEhQsXlFIXPT09NraI7OR2u5XT+rOzs8hms00xNWYepelThJwSuzSGyyL09vYaP69UKojFYg3fqTRNU6oS+P1+dqomZ/73r16HaFSVSkXZLRkIBLgF+QoYLgukP6mYKwBPTU01fHXYXC6nXOUaiUS4BbmJCSEQDoeVfpBIJBq+H2SzWaVobSQSYXWKK2C4LEIgEFCmBJLJpHLNaaPR77Axl7ro7e3lqKXJeb1e5bR+Lpdr6NP6+lZ8c4ByC/KVMVwWQQihTI2VSqWGnhqTUmJyctJ47fF4uAWZAFzc1KGrVCoNfaBS0zRl6s/v9/PWyQVguCyCviXZPByenJxs2E41OzurdKpwOMx5ZjKqJJt3DDZySaRcLqfslqyeFqT5MVwWqbW1VbnWd3p6uiEvENMXas0nsHt7e7n1kgBcfHqvPq3fqP0gkUgoxw54Qd7C8DvFIrlcLmVqrFAoYHp6uiGf2s6fP2/83O12o6enh1MBBOD9u4505XK5IafG9HVHnb7exH5wZQyXRRJCoKenR7kv2/xNuFHMzs4qnSoUCvFUPhn03ZPmKeJGXH/M5/PIZrPG63A4zFP5C8RwWYLqb7TxeLyhpgTmmxKrDlSi6oKu6XS6IfuBuQIBr/VeOIbLErjdbmVqTH/Kb6SntvPnzyuFKvv6+tipSNHoU2PVdxh5PB4eIF4EhssSCCHQ19fXsFNjs7OzSqcKh8PKJgYi4P27jsz9oJGmxvL5vFLiiXcYLQ7DZYmqp8ZisVhDHKjUFzDN0xu9vb2cEqN5NerUmN4PqneJcdSycAyXJXK73ejr6zNeFwqFhnhqk1Li3Llzxmv978lORfMRQihbc8vlckNMEWuapozevV4vp8QWieGyRPNNjZ07d67uO1U+n1d2iYXDYe4So0vS112qd43Ve62xXC43Z5cYp8QWh+GyDKFQSKk1Fo/H67rGkl5e37xLbNWqVZwSo8tqaWlRao1lMpm67wexWEyZEjNfFkgLw3BZBpfLhdWrVxuvS6USJiYm6nb0omkazp49a7z2eDwsVElXpJdF0tX7dRR6+3U+n49TYkvAcFkGvZCluc7QuXPn6vaGymw2i0QiYbzu7OxEW1ubjS2iejDfdRTVT/71Qko55+bVjo4OHpxcAobLMgWDQWWvfyqVQjqdrrunNiklzp8/rxwYW716NWuJ0YLol8jpZmZm6rIfABfvadLbrVfkoMXjd45lEkJgzZo1xpC5UqkoU0v1olwuK7vEWlpaWEuMFsX8/0VKiampKZtbtHjFYlGpBB4IBFhLbIkYLsukb8UMBALGexMTEygWiza2anH0MhfmA2M9PT3K34nocvQy/Ob/M9UlhJxOr4BsPqdTvROOFo7hYgG/36+cecnlcrhw4UJdTQmcOXNGmQpYs2aNzS2ieuPxeJSF/XqrGF59OZ7L5eKU2DIwXCzS399vbNmVUirfrJ0ul8spUxjhcBidnZ2cCqBFma9i+OTkZF2ceZFSYmZmBqlUynhPvxyP/WBpGC4W0KcEzAua8Xi8LhY09RP55mm8NWvWcCqAliQYDCpnXtLpNLLZrOP7AXBxId+8w42X4y0Pv3IWcblc6O/vN16Xy2WcOXPGxhYtTLlcxvj4uPHa5/Nh9erVfFqjJdErV+g0TcPExISNLVqYUqmECxcuGK/9fj/L6y8Tw8UieqcKBoPGe+fPn3d0ET/9JHI6nTbe6+npUf4ORIuhn3kxL+zHYjFHL+zrRSrNhWe7u7uVczu0eAwXC/n9fuXEfi6XU+5FcRopJU6dOqXc23LVVVfxaY2Wxev1KgvhxWJROTviNNWjq+r7mmhpGC4W6+/vV9YrTp8+rRxMdAopJVKplFLmIhKJcCGflk2vXGHuB5OTk47tB+l0Whm9t7e3o62tjf1gmRguFhJCIBwOK09t+jdwJz61nTp1SunwV111FYtUkiUCgQA6OzuN1zMzM468pVLf0KLvaBNCYNWqVQwWCzBcLCaEwNVXX23sMpFS4sSJE47ajqlvuzSfyG9tbeW9LWQZIYRSPqj6m7gTSCmRy+WUE/mtra0sUmkRhovF9PstqkvxO+0w2enTp5Xtx2vXrmVxPrKMPoo3X4+dSqWQSqUc1Q/OnTunjN77+vq4Dd8iDJcacLvdGBgYUOqNnThxwjGdKp/PK9uk/X4/1q5dy6c1spTL5VLq7ulXOjipH5gPD7OenrUYLjWgL2ian9omJyeRTCZt71h69QDztsv+/n5uPybLCSHmXNswPT3tiMPFehXwUqlkvNfX18ftxxZiuNSI1+vFwMCA8bpSqeD48eO2d6rZ2VmcOnXKeO3z+XD11VfzaY1qwu12K3XqNE3D+Pi4I/qBefuxz+fjmqPFGC41oi9omkcvExMTtq696OdazFfQrlmzBqFQyJb2UOPTq4abRy/xeNzWtRcpJc6ePausOfb19aGlpcWW9jQqhksNzTd6effdd23ZMaPvjDl58qTxns/nU9pHVAsej2fO6OXMmTO29oPqUYv58DNZg+FSQ3rpenMhv6mpKdvK8R8/flwpR7N27VqEQiFOBVBNCSHQ3d2tjJATiYQt5170NUfzWsuqVavQ0tLCfmAxhkuNeb1eDA8PKztmjh07tqKnlaWUSCaTyg6xlpYWDA4OrlgbqLl5PB5lR6KmaTh16pRShbjW9NP41TvEWKi1NhguNaaf+I1Go8Z7iURiRe970TQNR48eVZ7WBgYGEAwG2aloRejnv6rL8U9MTKxYP6hUKjh58qQSaP39/TzfVSMMlxXgdruxYcMG5TKx0dFRzMzM1Lxj6SejzTfshUIhrFu3jsFCK8rtdivVK4CLh3nz+fyK9IOpqSkkk0njvVAoxB1iNcRwWQH6U5t5UTOfz+Po0aM1XdSUUhp/jvkK4w0bNnA/P604IQQikYhSe69QKCiVuWtB7wfVFcDXrVvH0/g1xHBZIUIIrF+/Xrnn4uzZszUtyS+lxNGjR5HNZo33+vr6WJiPbKPX3jNPRU1NTdW0uKumaThx4oSymaW7u5sVwGuM4bJChBBobW3Fhg0blEXNd955pybXwOp7+avLvFxzzTW8upVsI4RAIBBQDu5qmobjx4/XZHpMSonJyck5t0xyWrj2+F1mBQkh0N/fr1wDm8vl8NZbb1m6e0zfFfP2228rpcTXr1+PcDjMTkW20m9tNZfkn52dxfHjxy3dPab3g7GxMWVaeGBgAIFAgP2gxhguK8ztduODH/ygUstrcnLSsvUXKSUKhQLeeOMNZRqgr6+PZV7IMVwuF4aGhpTpsVgshtOnT1vaD44dO6bskuzp6WFxyhXCcFlhQgi0tbVh48aNysVcY2Njy66cLKVEuVzGkSNHlDsq2trasGnTJl4ERo4hhEAwGMTQ0JAyTXvmzJllb0/W+0H1emNraysGBwc5LbxC+FW2gX72ZWhoaM76y+nTp5fUsfQO9eabbyqXgHm9Xlx77bU800KOo5/c7+/vN97TNA3vvvsuJicnl9UPjh07pjxgeb1erF+/Hn6/n/1ghTBcbKKvgZi3J1cqFRw5cgRjY2PQNG3BnUtKiWKxiDfeeAOnT5823ne5XNi4cSO6u7vZociR9N1j3d3dxnuVSgXHjh3D+fPnl9QP/u///k85he9yuTA8PIz29nb2gxXETd42EULA7XbjQx/6EIrFotEZKpUK3n77bczMzBjnUS7XIaSUSKVSePPNNxGPx5Xff8OGDbjqqqvYocix9H4wMjKCcrmMRCIB4GI/GB0dRS6Xw9VXXw2Px3PFfpBOpzE6OopMJmO873K5MDg4yHUWGzBcbCSEgM/nw0c+8hG89tprRsBomoaxsTHE43GsX78ePT09xmEvIYTxJKffzTI2NqaUD3e5XFi/fj2Gh4c5v0yOp/eDa665Br///e+NU/R69eRkMomrr74aHR0dcLvdRkiY+8G5c+fmXFnscrkwMDCANWvWsB/YgOFiMyEE/H4/tmzZgjfffFO5BjaVSuHVV19FKBRCd3c3wuEwPB4PCoUCpqenEYvFlB1hwMUCgR/4wAcwMDDADkV1Q+8HH/zgBzE6OqqcS8lkMnj77bfR2tqKjo4OtLa2wu12o1AoIJVKIZlMKjvCgIv9YGhoCH19fewHNmG4OID+5LZ582aEQiG8++67RmfRh/vpdPqKv4++K6y3t5dTAFR3zCOYYDCI8fFx49yLlBLZbFbZ/XUpbW1tGB4eRiQSYT+wEcPFIYQQ8Hg8WL9+Pbq6unD06NEF3/vi9Xqxdu1aDA8P83AY1TW9HwwMDKCjowMnT55c8K2VXq8Xq1atQn9//xXXKqn2GC4OI4RAZ2cnrr/+esRiMZw5cwbxeByFQkE5XOZ2uxEMBtHX12dc+qV/PlG904tcXnvttUgkEpiYmEAqlUKpVFKCxu12IxAIoKurC729vXy4chCGiwPpO2j008TFYhG5XA6zs7PQNA0ejweBQADBYNA4GMkORY1G7wfRaBTRaBSlUgn5fB7FYhFSSng8HrS0tMDv98PlcrEPOAzDxURKienp6RW9JXKxXC4XNE3DzMwMZmZm7G6OQj9nQPVNX+ez4477hRBCQAgBTdOQy+WQy+XsbpJCP8jZ7Bgu7xFCIBwOY2JiAhMTE3Y3p255PB54vV67m0HL0Nraing8rpybosVxu91N3w+EXKk7Rh1OSrli1602Ov3JkuoP+4F1mr0fMFyIiMhyPF1ERESWY7gQEZHlGC5ERGQ5hgsREVmO4VInKpUKZmZmLL1jnKjeVCoV5HI59oM6wHCpE5lMBr/4xS+UuyqImk0ul8NvfvMbxx2cpLkYLkREZDmGCxERWY7hQkRElmO4EBGR5RguRERkOYYLERFZjuFCRESWY7gQEZHlGC5ERGQ5hgsREVmO4UJERJZjuBARkeUYLkREZDmGCxERWY7hUgeklIjH45iamkI8HoeU0u4mEa04vR9MTEywH9QBhouDJZNJ7N27FyMjIxgeHsYXvvAFDA8PY2RkBHv37kUymbS7iUQ1Z+4HAwMD+Mu//EsMDAywHzickIx/Rzp8+DB27NhhXIpk/mcSQgAAgsEgDh48iG3bttnSRqJaYz+oXwwXBzp8+DC2b98OKSU0Tbvkr3O5XBBC4NChQ+xY1HDYD+obw8Vhkskk+vv7kc/nL9uhdC6XC4FAAOPj44hEIrVvINEKYD+of1xzcZj9+/cjl8stqEMBgKZpyOVyOHDgQI1bRrRy2A/qH0cuDiKlxMjICMbGxha1E0YIgcHBQYyOjhrz0ET1iv2gMTBcHCQWi6G7u3tZnx+NRi1sEdHKYz9oDJwWc5BsNrusz89kMha1hMg+7AeNgeHiIG1tbcv6/FAoZFFLiOzDftAYGC4OEo1GMTQ0tOj5YiEEhoaG0NnZWaOWEa0c9oPGwHBxECEE7rzzziV97s6dO7mISQ2B/aAxcEHfYbi/n4j9oBFw5OIwkUgEBw8ehBACLtfl/3n0k8lPPfUUOxQ1FPaD+sdwcaBt27bh0KFDCAQCEELMGebr7wUCATz33HPYunWrTS0lqh32g/rGcHGobdu2YXx8HHv27MHg4KDyscHBQezZswdnz55lh6KGxn5Qv7jmUgeklJienkYmk0EoFEJnZycXLanpsB/UF4YLERFZjtNiRERkOYYLERFZjuFCRESWY7gQEZHlGC5ERGQ5hgsREVmO4UJERJZjuBARkeUYLkREZDmGCxERWY7hQkRElmO4EBGR5RguRERkOYYLERFZ7v8BaZ5U6btgEu0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 500x400 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "7cbec774",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fixing (0,0,0) with x^2, r2=0.9999921393183026\n",
      "fixing (0,1,0) with x^2, r2=0.9999940727994734\n",
      "fixing (1,0,0) with sqrt, r2=0.9998914314178492\n"
     ]
    }
   ],
   "source": [
    "model.auto_symbolic()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e14000d8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 1.01 \\sqrt{x_{1}^{2} + 1.0 x_{2}^{2}} - 0.01$"
      ],
      "text/plain": [
       "1.01*sqrt(x_1**2 + 1.0*x_2**2) - 0.e-2"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.symbolic_formula()[0][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "031fabd6",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "train loss: nan | test loss: nan | reg: nan : 100%|█████████████████| 20/20 [00:03<00:00,  5.70it/s]\n"
     ]
    }
   ],
   "source": [
    "# will give nan, it's a bug that should be resolved later. \n",
    "# But happy to see the above already give a formula that is close enough to ground truth\n",
    "model.train(dataset, opt=\"LBFGS\", steps=20, lr=1e-3, update_grid=False);"
   ]
  }
 ],
 "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
}