Examples/Burgers-Traditional-Model-Selection.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e1412bc1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sklearn's version: 1.2.2\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import warnings; warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import os\n",
    "import sys; sys.path.append('../')\n",
    "import numpy as np\n",
    "from numpy.random import default_rng\n",
    "\n",
    "import scipy.io as sio\n",
    "import pysindy as ps\n",
    "from PDE_FIND import build_linear_system, print_pde, TrainSTRidge\n",
    "from best_subset import *\n",
    "from UBIC import *\n",
    "\n",
    "from sklearn import feature_selection as fs\n",
    "import shap\n",
    "\n",
    "def evaluate_coefficients(prediected_coeffs):\n",
    "    GROUND = np.array([0.1, -1])\n",
    "    errs = 100*np.abs(GROUND-np.array(prediected_coeffs).flatten())/np.abs(GROUND)\n",
    "    return errs.mean(), errs.std()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec2cb006",
   "metadata": {},
   "source": [
    "#### Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0e1d253e",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = sio.loadmat('../Datasets/burgers.mat')\n",
    "u = (data['usol']).real\n",
    "x = (data['x'][0]).real\n",
    "t = (data['t'][:,0]).real\n",
    "dt = t[1]-t[0]\n",
    "dx = x[2]-x[1]\n",
    "X, T = np.meshgrid(x, t)\n",
    "XT = np.asarray([X, T]).T"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a52b3e7d",
   "metadata": {},
   "source": [
    "#### Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "00fab164",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOkAAAC9CAYAAABWKFoCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAA48UlEQVR4nO19T4wlx3nf76vu92Zmd0kOhxL119JqV0rsKEAcWnQcXwwYVC4JkIMl+ZIz6UuCHAILQoBcHSnwIUcpOQeIJOcSwAhCxgcDgWNIpJ04kuBYXK4pySJFcXbI5e78ed315VD1VX1VXf1ez87Mm7fL+gEzXV1VXVXdr371/anqamJmVFRUbC7MZTegoqJiOSpJKyo2HJWkFRUbjkrSiooNRyVpRcWGo5L0EQIR7RLRc0T0u+r865fdroqzoZL00cJzAL4L4Fl//jkAr15ecyrOA5WkjxCY+dtwRH3RR30ewEuX16KK80Al6aOH30Yk5jPM/MplNqbi7KgkffSwy8y3fHjP26U3LrVFFWcC1WWBjxaI6HkA+wAO4NTd73g1uOIhRSVpRcWGo6q7FRUbjkrSioo1gIi+QEQvrkh/zpsrCSpJKyrWgGV+ASL6gs/zkj9/TqdXklZUXD6eBSAe+VsAntGJ7dqbc0p80HyEH+MnR9MPcQ87uHop6ZdZ91nTa9vG01/DD/6Qmf/xaAYAv/evP8BHx6nT9faPFrj9owUA4I//5Og/MPNAdR3Bbnb+lD65cJIS0TMAvgU3JQAALzHzl7M8d+CWs73IzF/TaR13uEmfHS3/Vf7epaVfZt1nTa9tG09/jX8wzmCPw2OLf/Ov9kbTZx959W9WlaFwAGC0sHVI0j1mvgkEwh4U8nxR9PHT4kl88ELTN7Xus5Z/mXWftf7Lfu4A0LPFgvszl+PxHURpegNxWSeANdikGfluqNUwGqOrYrZxZWn5e/T0haYv+0Evs+6zln+Zda9Kv8y6PW6vymDBWHA/+pfDO4M+J04iH/ciEJxKN3ye3Vxgrc0mJaLnmfkbI8l7APaJ6OvM/EKe+Cp/D4D78SY84HPFuuurdV9O3fv8M9zBW5Pz92AsYCfn98R7Mov7vAqLmTfQKNfpOPo8gCJJhbxEdEBEX9Du6h1cXWo/VFScB/boaezBDQ6v8Q9ur8rPYFisZ7XeWqZgiGh3Sdrz3latqHhoYBlYMI/+nSfWNU+6B7foO0CtvvimP5cJ3boYvGLj0YNxwuN/54m1qLveWfRCFvd5fzwA8Ir/qwSteCjQg7AAFdP+9E+OAeD6edW18SuODnEPr/L3sM8/u+ymVLwP4PvZ9VX5nLpLxb9nfm0bmOAhnoqNX3FUHUcV68QePQ3wlCkYQj8iSc8bG0/SiopNhAVhwetRRCtJKyoeAD0TFtyspa6NJ6nYpJexkKHi/YepNmkPg5NKUodqk1asE6exSReoJK2o2FhYJvQjNun//l/3gHOcgqkkrah4AFiYUUn6d37tceD9NAVTUbGJ6GGw4PXQZ+NJWh1HFevEZMcRU3UcCarjqGKdmO44qpK0omKj4RxHdcVRRcXGokpShWqTVqwTp7FJx1Yc/eWfHkwqYyo2nqTVJq1YJ05nk5ZJeuNXnwIetimYZVt2+vQvwO0ieGPJPkgVFRuDnten7q7rfdIvMvPnlxB0dIv9iopNBMO9+D32d55YF0mXfch26Rb7FRWbiN6ru2N/OZZ9kGlV+lr3OCKirxfSdrPzZIt9cRzV3RkqLhL7/LPQzzBpZwaDhW1H/zRWaYv+/JZPv5VvzLeuPY5Gt+zEii32q+OoYh047Zaep3yf9FkA/9mHRVvU++t+F8DLRPRFOL/MejfH9uL7u8z8ykiWpVvsV1RsIkre3de/8xZe/27YYPu6StrNLk+0RWY+8Frmt1DYHHsd6m5xy86pW+xXVGwimN2qI/338c89jV//nc/i13/ns0A6BXOAJdqi7/sv+W8mHehPUQBrkKRjW3ZO3WK/LmaoWCdOszPDKdTdVdriM4oDvwfgSzqxLmaoqFCYvJiBCZ2dRlJm/jYR/W6uLRLRi15YfcObhbdQWCuw8SStqNhEWD7dboElbTHbIH50EU8laUUF4FYnnGINgoVBV98nrai4IJQIecpFQpZRX1UTVMdRxQNDk1GHNbfInTC74D6/BUxczDBmk77x8k8nlTEVG0/S6jiqWA0GmBz5xtRWQyo3jQrO02zp2Y3YpB945mPAw/YWTEXF+WCMgRSjFRlFShZyJ5cABNJ8m/AB79N4d8+KStKKhwhjBKWgrq66nIqZaMLFKRgEWz/Y5FBt0vcBxmzHAAqHSEYCMwfSEQ2yptcBUckNeRPjFABwZ/KnDwmdrR9sAlBt0kcWY44ckzN0nIxGnD6gED9GxhiktEJCQtY9+giwmLqYoZK04mHGKumoDMNljpwxMqpCgmRN7EwiRVYKfiVokhICgZkoljMBlscdR+eNStKK88MEz6oDDc8CwaiQJSNjyE/ZxYqlpOxMovgmifEkpYyo3nNUtlmHWCZJ7/z5j4D30xRMxUOEKMa8ZOI8YZA3SLeE4DlRI/lIu2EpIyUAmJhOEtaEDQMGeXLqdCoMKGUwHFFLeOLvfQJ4P03BVMfRJWOlUwcIaqty6pCK11MdchjYjpltqMnImkwhv1Hqq+TJyeiJKMQjE+sZ/Ln4t49/AkyQgrzkq2rnjXW89L0L93rODQDPMvOXC3lGdxOsjqNLwIhTJ9ptGCZCm3sj+UYcOcEm1YSR/J6MZAhhMjPkIS852dmtxri2GwmzO/pytBTlJpbFjQExwI3Bk1c/Dbw5zXHUP0KOoy8BbgsVInqWiJ4vbNv5xfqy92VArdTR0E6dgXSMeYY2JA2uD1cm6alNKBIwJWqUbiIBKaijPt1QWpYx4ZyNcxyxiWlO8Pk041RyyW9DeROf3KNE0oyQY9uj7BLRDWa+VUiruDDETslEKekkh3ACeQfWhJOgtvfkqFVURUZNtkBIo2zEEdXU6HjjSJgRln3ZrIlovIfXp1up37h7h6FA4ilwknQio8+ItdmkfkvP/RGJGXYTZOYXdILYpACqXXpajNiQwR5UaZSeqIMio0A5V4qOnCSvsh0DMVLVlbRE1H9eJWVfBpOJcURAS04TMARuvNrspSGTATfSdgNuAbDLB+MkIRpHTGbCwbuv4eDgltzv9eUP1juOpordM2KdjqMv5AQULNtNsNqkE7H09Ss1T6jpqHia2IXJ9eQFnolV5GT0UhBAajsK+eAIR9peBNzR52GRnsaXFaSb8aTzElLye/XVqbQmSEQ2UVqKZLRNbK8rw1VvvY3NhnDtsc/g2sc/AxjC66/9j9srH/cSdfe9/3MbeNimYDzxvubDz+idAyfsJlhRRGZP5qoopZHpKalzT91sOmMwB6m8p+Ey08TCjYFIyWA7JmSMxGW5rjFgIbhBuIaNGxBgHLkl7CSjI5kltyg+IbGXnEk4qLPucYlqC+MfoUFUc0nCE54+E+yIunvl734KeJimYPy+Ll8loq/4qC/7eNnf5ZtwuwUmuwlWrEJKqoEYJSUzB6LTRwihpIxASpGeuQqKRMqFdE3GJoaFQDBG2YZeDdX2o5BKbEgj8fB2ZZSUMAQrhBZSeUnMfpwIxPTasEtHPPeEDTapPBLy4Ql4pLy73ga9WYjX+7sMdhOsKEOrpYM3P4RYSYQ+qnRNPsDZjsrJEsiq7EVS3lOtdgaiGi8ZvRrLEt+Q6/yNic4cr64yEaghL+Wc7RhVW6TkAwGNkoi+XPakdRIXkWyegCxS1z+KIC39rQuRddzK34ExKknPG3UxwyZgkoPHq6WSKFySExoeR506QrpANpPGB4+pmtJoTEpI76hB49VNY6LNKNMd3nkjaqo4bSwMuBVJ58ko0rNxN24bT0jGUHUNEtQ/J5NJSkBJS0SC53n00Yff+8kPgUn2JDnn00Ss+nKg/7TEDWCoTW48SR9Zx1FxwYCTBuSJlnQB38fGJ/6VHRkkJCm1FIgT+uQlJ6K9qMgYpjZ8PHvbUBw5srqHWy85GwNrol1pxWkjktUTVSQgN5Gc1qurEiaSfFDqrn9GiniJHQmVLmEfT5qU8uhFUfC/AVP8Pa5c/zQwwZ6UzbGnQH8LhoieJ6LnCrMcX2HmL/r0ZDpy40n68KO8YMBFpVIyhLQ01BdSXMmT2Iu5cycQD6mqmqumgCNf48QXmwZojVcxo+3IDTmiybSGV1m5FQnpJR95FdZ7U8Wr6vJ6aWYA2dCAPQkTCUcANy7REd0RIqyQzyUhxbKEpJKsd1QpaipZHI3kK4GZwNPV3aXfgvHO0+94cg6kbCXphSOVdiFO9aZoJ0JHhvwiWRMbEchsRBNJ2qSEDM6Ypolxjay8adzEvonxLq/xDpqYl40nns9vlUQMElLiKPWsBh9LkzpzEjXV34oVcvl4UtIyhNXjFSaGsmJU6XEXMdVhFPKfzibdzc6fys7FZyNfHvyy99UAeAhIuvE26cgC9GhPlh08xmQRWhICibo6cO4IARtnmJFWTUUyBueMk4rUeLL5BQGOTMZJzibOQVoJN15tbeDiiRzBGviwI5kVaQfJ6zt8E502EhaysNOo3Ulw7kRiagkJifPXW0nKiJcrK4yUeHk4oZdKvPf6RJu0IEkPv38LR99/TU51GQdY8i0Yj1f9h5teBvA8gLCGfeNJulE26aoFA0YJzSARxxw8qYQlTcjcm6qmMiiQz4ACCRtHTpkCCSqrJ5Mx4NY5e6x4UBsDawjUApYMbOM8oL0nHzwZhXhCUlFjqXF7z7Jvos1tR3+boobqZxXyIT2ylpIFVTTYoypO7xnGWf5sGAyZxuQfA7j6yWk26aBCADu/dAM7v3QDAHDwX/5Il7HqWzDfQSTxLhypAzaepJeOifu1BlqKRFDnYR1quDTONeopDZG8aJro8Glk4r8Bt02cvmgap062jfO6Nt6REyRi48jSGNjWL39rvQrbeIlJEk7VVfhpjyAlKdqObJQKqtRUIWIgksSJPSn6p6iwiPm0appoHOonKEUGqVr4efJ00gWJeoyUvKcBMybbpKu+BZOl5+vdK0kjBrqqC5KorSkh4Q9RIKZScuDgySRlUFFBQQIGe7JpPDFaJwE9Wbk1wbZ0tqGbynDTHE4aipPH+r9AwkaI50ncIpDOKmLKlxOEjIlklDlIJaaSsImPYOCsyfqzJqomkwVABcIMCCXXc4zKpWSJgCEfp+mlOpfidI6jpd+CGUsXVJLGBakxTtuLqjfFl/5VD1POoPTlYp9H1FZvP7p5Rb9dh5eMZAy4bbxKGckoYW7IScygulIgJLf+jY7WEc/65XO2hVoQ4ONbZS8KCcM0COIUhyevu2lFuIKKWiJkIrAovSTJUyAG5fFK4nFGwpKSk+dJCJgTU9U1uG4VxJBeAzaepOfmOBpx8GibMglLyEiMkpSAIqF/y0MT1BBI1FQhaFBj/fSHSExDsLM2LIXjWevUzrZRZFTEbF2+viHw3Dtv/MS/NTHMJpLUEsBtJA03iqyeZHKflmI4SM1RlT97xl7bkKTcXgw8Z/UoObMlS8QteYAKRyqdF/LS2LU4neMIdcWRw5kcR2MLBkQDTRPUaeGNEL0dh16946UjaY+rcuoEYoo92Tbg1pNyZmAbR1ieEXpx7szctbb1fwbAzBPTAJgBvfa2mkhM6z2nLNKQhLAAgaMEBRKHjgWHgYclDPa2Nft4Fx52TfKPmwIBRG3NiZJ4ujVpCClbtXtWEZ8lvyJlTuKcrKGFjLRcfY2Pu/bxiYsZEAeki8bGk3QaODVyJJaiBBzbYWB8jtKfk+ShMJfoJKdMdcAdm8Z1orb13lXjVVQnSe2WE2XcNui3vM3ZGkfIhtDPHCHRAP3ceIlHsDMEb2s/82RovB1pUrWVGw73nUx7GE8vkvTsswqBgD7M6QvgkiZxqarKgRgUmBifMZfsPV6tloam6bxLiCl5E3LGsSaNyweFnMhTwBiVpEd/+SrwsL2qdvFQBNXST4hLI9Mg4RBFK3kVViQnyUQ/vFQUB48R29GrrU3jJaLYkQY8a51t2DaOjK1TV7uZtydnhD6osQgOHmdPIi4OMI5o3LBf7cORoCYwZDDHqDUGJnZykZwiaj3xONLThVkPYhQkl+Zk0snzeZGMKFQKA0OylSRelmdQ1hKJGMqgYVxoP7K8WfQyEAM08s2Ync/cBB6mV9XOioFNmjzgKPlK3wJxMx+5hKUgKYtvhPgpjmSOUryvTRPmHLn1nldPRmczNk6FbQk8b2BbR0KeeUnpiWlbxLCB87Q2DGvgnDvGqaVurtIRMkhG43rZYIpDHgf5R0RD7VFIqbUJgs5IidQaFBDiOVaU5C/EoSAhS3UUVONh+Q9wXUak0JbMYJZy3/vx+9QmnfAGwGh6apNS7JByHsqAF5K5pExJSpnnNZBS7EmRot6OhHHS00nD1ntRW2DWoPc2Jebkwq2BnQF966RmP/dknEdiYgb0hr09yV56KsnYuN4jiwjg7UhCVGPTgUruO38qNIzznZeVWpo+RQ+Of0VJqKRnEqeuGVNdk/IlbkTyDuJwiuv0feWScmzQAPDYRycuZmBM+vraeWCUpET0OIA9Zr59lgpWvQGw+g0BpbMpFdVJhdKWH6nqW1zJIxJSeV1l2w+etc7h07bRppw3sF6NtfMmEK/basCNI2W/5Ym45eclW0I/56C22hnAxjoyemkpzh7xsDj7ksOA49RVDgvMy9MZjhyivsbOqHTUZCnPSEf1Hc6pthRJsJSsiES3hXQ8GDGLtmWpDBSuU3Gl+yyReTAQTIG694vGMkn6274p/5GIfgvAyw9I2KVvAKxKJy/xwphP8ShS073epPNQNEJFdQWCo4fEjjTkJOWs9R7SBnbbk3TWop8Zp77OG/RzoJ81sHOnuvKM0G07dZVnQD9z5LMzJxm5deorWk/U1vcaw0DLUU1VNqW7tWgjBu0gEM8bgiwuGqeeupMCKUVgsrsmdn4aksNKp45kJ8vlvJqYhc4/sAFVOM87ml6Qcsuka3HgyctakvdBSDpmk543lpH0JQDPENHjzPwHRPSbeDBjeDc7z98AWJ5O7k2NqL1Kz1USlshxMkhMZU82xqm48haIPzqbkoC2gZ07+5JnDfq587jauXF25Nw4gm756ZBt5+CxLWDn7B1DQD9jJ+1mvgd5goK8utrYKPVJ9RQ/yIhRTZ5woc/o0VoIFsjBvhMq9VPsJOvzMhKpSANVlaLapjuez8vw+awib06GMUkr0lmFR/MuK0+njbRhQDxG9Ozm5WbxxSmgFSDg9MR+QCwj6W8BeA3A14joUwC++oB1HGD5GwBL0w9xH6/a/wsA2Gs+hL3mI1Ft1d/6kB0GTFzhw97pY/0UiVsg4MiIWeNsxy1HzH5G4Llx0rEl9FuEbgbwnNBvOVLyDLAzhm0YmAO2tZGMDTupaACQBTUA4FVXw6GZ+ghAzbVFlZSDJIyEikd26dadBlKyJptKD1JS2YxAQuIoRbN0JjGLU7JjSLxcIiaqd0aUkvQsXp/d+0rpOULYnIRaA3B5GHd/+kO899NXJeY6VmGJ4+jwh1N3d5iGZSR9DU7F/QMiegLArzxgHVPeABhN32kew6fnf9/N6yU7C4haGz2xLORsWu/s8SSdNbCtAc8b9LPGeVu3nerazQn9toHdIvRzQjcXMgL9FoNbC54D3FjwDMDMul95TiCygOlBxoXJuN5gjO4tGhT7iicgMwVSMhxpUxWWwryly+uqZ62+WlLkipJUJOVAkmrSJXnd9WOqa1IOgAHJlklCW8ibSzSocIl8K6TpgNyDAYUTcsbrGY9/6CYe/9BNEAhv/Nl/v41V0PeeYefmKd6kmYBRknpy/rI/fR7AHQB/dNoKTvkGwG6+rQS1LWhrBpB7JxJAeAuESO29Ix5YeVVr5sI8M+i3Wke6eeMk5ZzQ7Xj1dQvotr3qOmP0WwzMGDxj0MwCjXXENIymYZjGkdGA3TGMGW7GEUi5af3kpZBPCBrO/Y/t3qpwaUK8kLdXEtVLSdJ2oVZxJV4TysbraaDaUgwTylK1IG1DHVIWZ+WSuq6kXkp7MUK2FQQc9R4z0oSQV50rrcXF0elVV/0sLxhLp2CY+c/98d+dpZKzvAGAtgVtbwdiguCOYbF5E98C2WrdYoC5V2HnjpD9tkE3Iyctt52Tp99m2C2GnVvwjnXq6oyBeQ80Fm3LaIxF01g0hkFkYQzDGOdHDU4eP0cLuAVzluE3qCJY37Gt1aSksGdreN3JkxY2StXwHphIOQmDisQUaRHJ6e1c3wG1NCNl1ybEK5ENw3yjTqJC3Mr0PC07H1V9wzkPr/MaSBpHad1JWpn0y7BsMcN546FYzPD/jl7B3tVP4slr192USOumPqy3Me2MwLMG3Zafr5wT+h23kKDfBrodZ1P2O46YPLfAlgVmPWjOmM06NG2PtmW0TY+GLGaNBRHDELs1A+TMTUtxBasNZCP0SoIxuykRa50kZZUvbLvhX3ViRUb36hMF4jED1KuR3ypvrnbkiERMSKpU20FelKUkMJSkujyJRyFuBflGybaM0EhV0iEZx8qlgTpduk6WUmnivvvm1MUMeLQWM5wFW1efwqc++hvgxqBrCLY1QOvJuNXANkC35aRmv+1X9mwD3Tajn1vwNsHu9OA5g2Y9sNVhvt2jbTrM2h7ztsfM9Ggbi5YYBOtmahhgw2gI6Jk82dx2lJaN3y0OXsXlcG4t0FshMFxeq0hqhZyORCJJIwkj2Yw4g5TUNNrZo6SmOIhSYvmBo0d0Vg2IHMMgl1cwlm+g9no1uXhdlhdQ6rDOC12PyqxJaFldR2m7szqK5C1KXVef5H3iA9MXM1RJ6tHvtDjZ3Q5viXDrFqDbubMpFzsIHtjFNsBbFnYb4CudsyW3emxtL9A2PXa2Omy1J5gbi62mg2ksZsRoyHoHMUeVlYwjGhOMNegIbsdy6/pQz+4jsn1PsGzQWQPrydn3Jqi1fS/qrNiWFGxPthSkIzGAnkCBcNqO1M6gIUlz23JAIrskLHVIp9VOpJJ01OESwQrhpemeOMCQoK4OZT/qdcX6/kKcupalbErSQpm6jgeQiAQkEvgisfkk3TY43pvBthSmR+zcq7Detuy3LXjLAjs9sNWj2bK4snOCebPAznyBq/NjzKnHdtthZiwa6tEav4qVACJ2JEMTPh/Qs3cEMdAz0FtHzJPe52GDrm+clLUGfW+C+uukpoHt4SWjkNtLzT56aZ03VREzSK6cjBRJZYWMKg8QO/8y6SeOnJJNihHi5qTTpFqWV86BoK7G61ilK6kozrec5IM4HpwneUBInGgDZ1FGaD1oTIF/liXcf219UzAbgcUVg/sfbMBzoNsC+jnDbvu/uQXtWDQ7C7TzDtd2TrAzO8aV2QLX5ifYahbYMR22TBeIGZyOREFSLriFJQZbC8uNJyXhxBIWtkFnGyysgbUGi96gYwO2BgvrCGpFivqjSEgoFZcStdZJx0jMqOJGtTWqsHqqZYyMidpaIuvIdUm8hFEoBwXVFKXzSL40D0MvpAiSH1BxJQnMIC0RocuK8aTT9bMgzkio2gPdHsFEqZoPRgpXJ26wPRUbT9LDk3385Hv/DdufuYn5J2+AZ56Y28fY3urw2JUTPDY/xJXZAk/Mj7DTnOBKs8A2naA1Fm5/ZwsiQucloIXBwhoccwvLjli9NVjYFse2waIndDzDcU/ouEHfO5J21qDvG3ReUva9IyYzwfZOYjqbEwiT3YqEQV3tfbo4eyReVFFOCZvMbSJVWUs25oBYeqpjlfTLSVoMc1HdLdYhP6T2tjKUHelPEnUUsQBO523ddf5ZkbpuICUdiclm5FXh+P5rTHzn57eACVKwencVzMd28fg/+U1gm2G272Fry0nM3a17uLa1wO7sPh5vD3G1dcScmx4z6kHGeScsnOTr0YBs48KWsLAtFrbBwjY4si1O+gYLbnHUteh8vpOuQc8NOmvQdV4N9qqtU2cJbN0mQKwJ2fsFCBZAT27BkUhCTVhInEsnNU9KOfmAUadPQkqtEqtOudTmzPMNCMcjRGaVNyVDWSLqdCES4Ea1zF7kXA3GUEUV1upBQzciSEy5VtunpEgqLSPs7t0AHpa3YDYFzbUFdj50DztbJ/jA1fu4NjvC3vw+npzdw2PNCR5rjzA3J5jDeWkZAJNBx46AJwws2NmOR7bBcT/DghvctzMcdTMsbIujvsWJNehsg+NuhoUlR9S+QW8JXd9Edbb336XkBtwjqrJ9lIyOmJn0zNVWIVgyxYKEsNqOBDIpmaWHBQM6r06fIDWLZFTqa74IIFFfWRMnHkiVm1yX24uSR0vUUI6rYJjXS8skr4S1Ohsla7lN/gZOuR/KaRxHq17XVPm+ysxf1nEbT9Inr97H33r6DTw5P8KHt9/DNbqPJ2aHeMwcYtssMHe9HCCnyi7YYMEtjrnFggnH3OKwn+OQZzjq5zi0DU7sDEfdHIe9U3GPujmOLcH2TqJ2lmDVsbfOM8uWwL1sIERAJ7amszGZyU1hiIrWQzmFlPqaOIiQpgOBgHpaBVBqrlwHde6vA2J6CBc8vrnKqiVjus43s9vYq4kJsePeDvF6kXLqXNdNHLQGUCTdwP4UmzQMEHrQsMV4AEl7Ivm8LUuqLfq5kFKBV+A06u7EDzbBr7q7kcdvPEl39t/AO/9pH9d/7Sp++dcNdsjZnA31MMSwnpgn3MLyHAtucGRnuG+3cNS7491uC8e2xaGd47Brcdy3OO7nOOocoRddg5PeoLcNuk48tQbcE2xv4lSJJZCfUnHSEwhOH+vVWvHYKimaO4AoJ6SWREpKalImxITKn9ubiVRzGYsqrEiahKS5pBRyFGxQKUt7Y60QQCSTjyd9rTh9opSLJNRkUp1ATaOQskM1OQceWgYSm1VLypyHShW/sz/NJj2lurvqdU0Q0Q2fNsDGk/QXP3mCf/bPr2G3fQ87dOw9tIwOBj0bHHHridnirt3Ce90WDnkb7/bbuNfNcGy3cLebO5L2MxwuWucg6lucdMartS26zjuDOhMWHDjb0jjp2Hti2ig90cN7aaMDSEgoKm3obL2SaIM5TqSk86QJ0lFJylFHjl8GOH2qhFIS5pIySCgbNYMgpTgpP6/D5SlJvtheeEmakDjPq/MHqWzjBQOy6faRKo+TOkrOKYl48olPAQ+4mOHe6z+ULUGBlOi72dX565qAU4NfSj8C7bDxJP1gcw+f2TLYJoCJ0QPoGFgwcMIG93iOd7sd3LPbOLA7uNtv49Bu453FNo76Fvf7Oe4t5ljYBofdDMddi84anHQtbG/Q9Qa2d9Mo6Ancmej06UxUZ3tlN/Z+9W4PR2ACqBNbKyOjlpqy8qePcfmL0yVnkVZ3IxkyYnnnS0kqph7XnJA2rQ/KXrQMTR5dtyDOQeZxmuiclZ9JTHWeqKu6voRkGUEzskZHkZKcekBJ2gYgvNeHySipu9c+/mnZEhQ//5/JmzQHWPI65pj6K9h4kj5uCB9sCQzGMROOeuA+E+7227hrt3HHXsVBv4P7/TYOuiu4u9jyKu4cR4sWR/3cOYh6ck6hvkHfe4kp3lkhZu//RJUVqZlNk4Rwj0QiRpvTN17Hl6SmxGkSErKFBpxd5zPLYgYdTrytWiopacFp3lyNzdXCYMNZpNJP4MkQBoZAII73GWxOqLx5GQUJGNL1dTpN3Udi40o+3Z5UsqZ81PVOZGpG/BVY9brmvrwFBuAGET3DzK9I4saT9Cc/tvj937+LX/2HM3z2H8zxLjc4sDt4q7+Gd+0O7nTXsL/Ywb1+G3cXO7i7mOPIOul53DU46Wc4WTRY9Aa2b2A7wPYN0Lmt3VkT03tmg/TsSJHPq69+bSsJGXMVNqwGQqqe6lVBOYmBKOG0VHOLhYeOHC99tEQhVp1R1F9NlESVFOlJQ2JoUrO/CSgpmtSpzpk8Ufy5ImY6H6nAmhKa3JzmUXyLBEzbIsSLNOM0PXEySZt0ugvceWe6TTrVcTThdc1X/PnzGKrGm0/ST/xCi3/xL6/iHQbe6Ldwp9/Gz/sn8NbiMdzttvF2dxXvnuzgXjfDe902DhcNTnonRbvOL91bNLC9cdLR25zUE6gzYaqEekRVVhPTE9Foz6omqPzoospq6ZhLUi0FdUcX0uk5Qjm3ccjWpEhImHQWTusS8osTRUtPUR0tp+lCBE0KPeAoR1DMy+V4geX4gFQwmpgpMaWsQGKrEjXZpSy5Vz0gJA4rXwhTqC+o3+qavWvXgQk2KQGpRrECq17X9OffALD+L30T0S6ciL8B4Nl8DsjnuQPguwBeVDcDALhve7zRE96y2/jxYhcH3TX8rHsMPz++hnvdNu6cbOFet43jxQz3Fg0WXYO+b9GdeGJaA5xEYqLzc5i9l5SepKZP1Vk5Gul/vZKOmc0JKBslEA5RCgze6/TniHkZANl8WiSzF1l1DLkuSD4oKaHaIQ4XssFM0+pmlJqR6QN1Vs9HFju/HxTy+DAQyHWqnYobof0FTTdIV9EQCnWHdmpVNZ9vlTiCH/hUBacgm77sUVpx9CXAjRJE9CwRPV+YzP3imOG8b2f4q5M9vNnt4m8WT+LOyTbeXlzFO8fbuN/P8d7xFo4WMyc1j1s/ZeKICa/KypHEG2sJ6EQ6UpCMWkJCEXOg1gbbMdqLjqCuozjpGTukSEC9YwGx9ppyLAP6GkTC5/GSBlKqbp7u69US0qqOGyTzsgUBkkenxUwlx1F5id+wXU55VyQp5lWJehBK0gqNz+sQ9COMZC3eJ0C0gDXgwkmaEbJkNAPALhHdYObBPNGbi8fxF/c/hre6x/Hm0eN492SO9xbbuHe8jeMFYbGYgxfG/XXGLTDoyXlmey8lF8qe9Cos9VH6uQUI6tz/AIl09FIuTLOA1bpatRomEFdJMSGJdHSJcw8oJaFSA4NdqqRMauNmIiUhcYG4OQFziYlMYmJ4HTFn8apxYRVRlkc7dUJDYxvEN5W3u9w2dxShmFeTqOgg94OxLiO/OX9fYJzmlbVkei3D5F3wJ2JtNqmfrN0fkZh7cB6urzPzCzrhr/66xV//+zdx1O+j/aWbmP3ih3ByMkN37BxBtGiAhSOm6QyogyNiB5heEVI8sRlBk7dPGKBeqarBnszmIVmVJTZWyMcJiQGl8kJdC5Vf4v0xXM+phMqnQBIySf3SEcMUBKmytVQiVYZ0YJsST7cz/8U4uw5KbU0y+VrUFFGo2LfHJNcAAw4l0ylOjUnb7facCmwXspGNYV1Goa79u7exf++2PJHrw4zD68bU3cm74E/EuZDUL3vK54FuZYT8Qk5AgUhbIjogoi8w87cl7fCJD+PqP/0N9CcznBzNsNgHeNGCjt3L2NQRzMJLzI6CVDSKmNQhkMNY93sZL/HYxwWSCXEJ0ZPLAKzTVcVuZFb2qrYPFQmdSqvKCGH1SUFNZKVWennhSU/BERMkBuSeSpJN6siII+lEoW1BiinJ7pxIVkk5JYVCOkM3hvK2JbajFnGqrrQThHITL20uYbVNqR1c0Jk5uQ+BlqPR1+Tqe+rqJ/DU1U8AAG69+ce3hw1MUZonvSicC0k1qUrwxPuaDydzQN7t/F0dp3H/3hbe/fkV8KIBHRlQb0ALIaS3M/3RyKoeFunpO1AfH6jpAXhyGGVPxrdHODp4RA2GqKxKaoY4CSOGlWNCe2+B2PG0oyWRpkE4cFoPhuHERsvUYrD/lAyV0mIZRqeFcqPxrO2/dJGCakfuiS2SEADCxzCQNIjz68RUKJSTlKmeSyndllmk80YVOX8IK8B4dGxSPzf0VSL6io/6so9/0bugvwk3gSuLkFPC329Ab22h8RKTenJk7Dy5eobpvbTxqqrES6dPFxQo4vpOF+1MKAdQJK44jILaJx2SAdl3J13VoySY7suyR4/u0NrTWOroSghRTjZ1nXhOyw4gpB1el1fKLxEjJCGlOmrbMJV+qv15uUmbR4g2pp4mxSQyEe6px99omFpqJ6AHh8l42CTpMniV92Yh/vP+eADgFf83kMj8N3dw97++iKsf/TQe/8hNQIjYu04f1I7eOs+tK9RPqXjyCjETMnL4oWWPoOgMctdpUjr3P6WkCs4bpK+LDexITn9/i7hCRnfanGw50QYdX/JSQfVL26Gh3zRxKnUuTOP9ZxcOyieVlA4QhdFk4LwhVS8PiIUkZ+H+JDVUk0r+7DvJxXBsl8Pb918Hprz0jWyQu0Bs/GKGK7Mn8eG//Y+c3fiOjQvVvTOHLAfnj9hRokKmc5VKEoW4zJOZ58lIQQMSyL90sbpO086geF3esWikA2IoAaHKTVRC6eE0lJwleFstfto0GzR0pezJJGxeZlsW6kiW3Im3Veze7D6N+i6OxKdOL93GkXAYGEXFwUhY3YcfIJ/a/gWg7hZ4OrSHjJ1978Gxyk6zcUWOVrGC40LSvdoWV/NwOIo7v5gueXQH83apvi7tsFnHyjt+EC5KZEh/FOktZVGsT3fYQXypQyNWqQeDtKNm6VpSq/sfkAmFfMlzjY8h3FsYQGQAHBmQkrLj71ROXxFmfRwL5/dbiBvBMsfRu288pFMwD4rmqMf8YAHx6rkOHdXOOMou6aicvbOoJWeulkknGjhDVLZMuiUjvdbmJFDsZ8POkEqMAoRMoQ7OEvI69PnqzpeolCWhWmrQ5LzT2nAu1wgs4sBw3lhC0iee3sApmItEc9Rh9s5xlEDy6UMZpnNSyonaSDmRWgBQsgdLKHSQfEVP/AKMqn7ViFwiEKlSxgSElD+QzIUKJnfMB+jEpyHOKaTT6PVnIeqYhgBgOKIOrNVRPHRTMBeJo/tv47XXX8Leziewd+UX4D6vhuHzzF+WTTp+yIRUjaJh3pwwyu4LjqGBOgekavDUu3sAqbUiS9rvHlBy6XlQQW7+ITtP7iMf/DjLVGpsKayuH5Q1Fj5luclvzHj76MfApMUMfLbB4xTYeJLu4Ao+s/MrruMcHrtIIhGJ3s2mbTNyRysrEnyaUT+Q3oIjzyvhkBfDzli6LrcX8wvDXIWz8YKDp9SePK+gUGwaLgw2XLoQq/vyaB2FvH4ZZFGiR9VGFVAamXTYotw4uS89CJTKYPUs5Pll4eRZRdf8U1sfA6rj6JRY9OCjIwTyBYICIELYbqL0hd7iufxbNcrqvBqnHakL6aOSG8vT87xrGsknYbQtY0ScohaM5J1628mAm5Wn0wba02rod4svGhtPUrY9+OQEIBM3oVJf+GbyEyWEQFgGqb5MQeAGCNmBEbWYpcDTN/hBrnkUUSKBPgJpeBXO47ocnA8Ep/jtGMnKsovExpP0sL+LHx7+GfbMh7DXfBiyqJZ9B3DCzhGK4YgaxkQJk3pZSUlH9mkpgfOfKzI8iUMmJ/OBYCrOTOpcApTOl6SJmrw07wPWodX+bP7YhacQI79OSKfVZlVWcAoC6WL7LD68/I0k7u2TnwBTd7BfkyKz8STd4Su4aT4LgGD7zktLT9BM9QWcFA3PXxFACAzoS8rSMt2xjcKvwYqt2hwTh1EMq+5DlJJe+mwsdjRczAvV5OR+cqIkBmMhfJq8D3rdqrzLDN+p12V1s4pbFWZVno97avYR4Iw26cFbrwLvp3lSBoPV2yR6Coa1Hht7ricZRxLKkj4PS8Nl3gMHEBDDubSVa8AgMmkd8skEVRbpskJzYwtIdTxKOqZOz8Nx4Il7TOu2K+k+Niio61i1k4mW5E2aNr0OIjBbZZLExfYhzABTDMfBsZyehGHPp46JIOZRdffJpyZ+qmIiNp6kYKg5T4K1DAovIKqVQmQBGEdqlxVhTasGGeg39Un/y/kB1cGypKjI2VGJXALFCk8H7Yn0Nnkgvw7r1Ud6UNBFFcIEDidkdTg+l3Cd6pu0KhyuU6uwgEnhtI3TrjtTHaf5WZZI0vPG5pNU6OGfaZirhPwkLp2ZQCTL1/womVwk6NUmAfKeplahICUoovvciqUifTVRhlCUFrvZtyshfMItRS6lzrIiI/d9DCPrW/petfNmzPGSeItRCHP+aM6OTXWunXaRxthWLOeMjSfpIe7hVf4ensQHsUdPq13nRL1z+YRsDOnoohC7rpx0C3mJOnRAlRq0RYK8EOyIoCStzxikjzaRkrLUj8j5IBClvDQlhqOcBgBLUWOwck+Uyg0W9V7VE1RgpGQeEFvHFVT7Iig+n1EP6jLv7tRpplJ5JailnKODUJKX0nOfd7/7KXABjqNlH2xatVnfWki6bDdAnz56Azu4ipv0WZ25XImNT811dHIh77WzST/gtBhS0ivMo6t3R0PejAS+7DAYsA957VfWIcT+qSvldH4/EDxqB9J+Ck4PIZA4SGI+8veRqrvINBA/uGnvZrg+GpVWzAjShFb36US7v4VoJ3MeDrYh1M4NnMQtC6snNRxgSoNJULHHB8ckXg8M/rjXfBiY5Dgat0lzTPhg09LN+tYlSUd3A5z6xakiaCQMAGC3RUiQpKQ6bMwcBmBNPvn9Ik8QF2rLtZEorrsIidgL6EgC1wxK6tDFxMtl8KA4kMh2m7pByYV5h5R6KdE6wmJzMMKHlfKH4AcOysOs26aqdi1VtxRYkmkHaiDJbn9ZeDTdj3o2OOvioKDNhZXh5Dcfxi2DW7s7WZQu/WDTqs361kXS3bHdADHhi1MPjFylCZ0s2mqavMmUhgyugQt+hzrSRIidMq8OUJtgBSkVSk/tvKCiufrYeyRYv8Xh66V8hX9ylF3k/f1KvboNANw7nWbYhrGwDBbWX5tXncwP5eH8WT0gSN8HvIbDcY8oH5/s6p/4MlReFQ7FS95Juj6cJM1s0jt3buHOQeji11XSbnZ16YNNGNusb10kHd0NECtuQGxSANEufRAkw/EIefWEu/Tn4rQJIonVCG2ZgzRhTxjSF6R7V6bg2ElSr6Emdd6R8h7lpEtMJXDYyc9Jm0RlHgWlA46QXiR+3jyb5ZXBSLd5FUbvScqNqn1x7jOpvDBvO3iXNm3Ufv8G9u2bcnp9ZXsZ2Sbb7ots/qtsuH37j26rpAMs+WCTQnGzvrXsFrhsN0CsuIGBTXreWDVy2ozEujORUg0h3T+VxqF/6wio3QJ1WaGooEsiiI1A1sxCy04T+xmc3F60EFUdQVLqZ8GpWafMBd1ybfMmKJCABup6JrZYBXSbkId1Osf0pNxBobHcHP7aPfow9toPAwBudX9xu5AzvYzHx9sCVn2wCcs267vw3QJX7QaICTewdqwykgAldRNmIIrD2LOG35wUaSY2UuyHA6dW6JBe7/adM9hmoVNT1taM3KNEp3iu5qNjU2Vhg7tPddXSKvIwhy0YI8mY/QIC8veNeA95PKuBkEW1lSktUX1XNaL4LNS5HpBXoaDujmdd/sGmsc36BOtQd4u7AaovShVvYCNRIm9RdeZsgbmEs07iCUEJJwrSRUlX/ZYbW1ZVcLGK8kCRjwL+GKKVZ4hLYd05Ke3oJUEW1N+UHFQKW8kaC2JFnjQcn7FuQqphZJJXaQvSuKDlTCWoXH6KBfbLPtg0tlmfYB27BR6gsBug/qJU6QYeSiQ8oEJY65Cqw2tpLAxMOkw8KRa7pIrg3LVWzThEBxQzQCbK9UwZzrSEUi8eUy+z6MSmRBomBvTefsv6/tgD4Ow4Fh6U7e34wTNfDrcscD1Ljh6+xQyPCkY7m6jRiOTlfJSn5FA4KVSRE51imDMhxwgGl6vah6M1ilz9HfbxTP3O1ehBeAkzS3b4qstOsaopaiAkzqPrKy+y09Xds2LjSXrhjqNNBhXImEwlKQk8IHIopBC/itB5zuF/UW2TqdikOUoVh6zc0iQfaVqp7TnhVgl3X+s0+xSK6Iw9+iAwdXH8yJzq/t3bwPvpLZiKDMvZlCL1SGWdfUz1lH8Dxfd0zUmaNiT5WNPSZnK5mcmiDWkyxQxhwFomuc8Iy6Nrd/eufBJ4X70FU/HgGGNNyUBLPVKKNRZh87dwTU7o82taUUsd2JaUNZmD4yyu1NLSvCCyCzZo8a2psfZXmzTikbVJNw1LWZNLDG03j+wVJIQHVtrODwo9J6yrHaZjeAtZnNjhk21S5tEPQp03Np6k72ub9GHAKgcYEI1V5ZAqerDHbOo1QJq+Rx8CeOLODNVxVPHIoOQAKznDxmznJT6fSG5lk060p88Eed1uDagkrbh8rDJKJ5FbDNSSPa3Jmw8Yo+xfDmagr+puRcVylMg9ak+rYz7vChQk+AossUnfPpz2+cSpMKuzXC7EcbTPP7vsplQ8KihKaxfn+9n1lWVYBvq++PfUfOIu+BOx8ZK0Oo4q1ok9enqi46iquxUVm43qOKqo2HCscZ602qQVFQqTbVIet0nRn++XnDZeklabtGKdmGyT2vXZpBcuSYnoGSJ6lYhe9n9fLeS5Q0QvEtHvXnR7KirOBWzBfV/8e/voR8BD9hbMHjPfBBxh4fY0yjG65WdFxUYivMg+xF478aNPE3HhkjQj39i2nrt+O8OKiocDbB89mzTflTvD6Jaf57alZ0XFEuzzz3AHb8np9VX52XL4Js9FYy1benp8HkCRpCu2/LxUx9E+/+zSBoZa9/qwR09jD67O1/gHqy8QSToRyz6lsir9XNRdZv42M38j+wsE9R+kGWv8895WLeII98+jiQ8MNbrWut8HdXtcX5mDvTQd+dPQn1Lx58+dJn1d86R7APazhsn+ut/058mWn1Oxav70rOmbWvdZy7/Mus9a/2U/dwBgO+7dLajBz8J9QgWIn1KZnL4Wm9Q7i17I4mTP0QMUtvwUHODnf/gSf/s7S4q/vmJe68zpr/EPxtIvs+6zln+Zda9s2yXWDTjSLMVr+P5XXrPf3x6UG6XwR1X8bpYv/xbM0vSNX8zAzP/4sttQUZGDmf/tKbIfYPm3YJamb/yywIqKRwCrPqWyNL2StKLiguH9LDdK34JZli4gXtPrNg+CVW7rC6x3F0s+j77Gdnz1Mur23vYbwOkdeedQ99p/c1/nC/rTJ5fV90rYWEm6yi19wfgSgM+pj0s9v8a64et8Dp4ol4Cv+HvfW+dKMH/PMr9+a9nU3HkiH4guue8NsLEkxWq39YVB5nr96Q3VjrXAE2Otdaq6nwfwHf9l9m+MLOO8KHwXwLdEki/5XOZF49L6XgmbTNLd7Lz4CfOLxNjn0deAsTXO68BNuGe9T0RfX7YQ5bzhp+O+DuBbAH5lXfUWsJudr73vaWwySQ8w7RPmF4ni59EvEkT03Aa8EfSqJ8zLANam6nu18iX/1tSBqJ2XgANcft8L2GSSXuoXwPPPo6+x6n0ies530Btrrhtwz12wi/KrhRcF/Rn638PlEWWjvj6/sSRd5Za+SKjPo79MRC9jjZ2FmV/x97qHodq1jvq/Dffq4HP+fJ2ezW/4tdzPAfjSGr27zwH4XLY09VL6XrF9mzwFU1FRscGStKKiwqGStKJiw1FJWlGx4agkrajYcFSSVlRsOCpJH1EQ0Y3LWHNccf6oJH108RzcWtiKhxyVpI8g/CqlF3B5b9FUnCMqSR9B+KV1t9b9LmjFxaCS9BGEf3Nlf1W+iocDlaSPJj4H4MX66Y5HA5WkjyZu4ZIW6FecP+oC+4qKDUeVpBUVG45K0oqKDUclaUXFhqOStKJiw/H/ATp3kTGaNnNrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 252x189 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "with plt.style.context(['science']):\n",
    "    ufig, uax = plt.subplots()\n",
    "    contour = uax.pcolormesh(t, x, u, shading='gouraud')\n",
    "    uax.set_xlabel(\"$t$\")\n",
    "    uax.set_ylabel(\"$x$\")\n",
    "    uax.set_title(\"$u$\")\n",
    "    cbar = ufig.colorbar(contour)\n",
    "    ufig.savefig(\"../Figures/burgers.png\", dpi=300)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f25fc2a6",
   "metadata": {},
   "source": [
    "#### Add noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "41704e4d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----- Clean -----\n",
      "STLSQ model: \n",
      "(x0)' = 0.100 x0_11 + -1.001 x0x0_1\n",
      "l0_penalty: 1.6123253173828125\n",
      "PDE derived using STRidge\n",
      "u_t = (-1.000403 +0.000000i)uu_{x}\n",
      "    + (0.100145 +0.000000i)u_{xx}\n",
      "   \n",
      "\n",
      "----- Add 30% noise -----\n",
      "30.0\n",
      "l0_penalty: 1.93969921875\n",
      "PDE derived using STRidge\n",
      "u_t = (-0.775461 +0.000000i)uu_{x}\n",
      "   \n"
     ]
    }
   ],
   "source": [
    "print(\"----- Clean -----\")\n",
    "library_functions = [lambda x: x, lambda x: x * x]\n",
    "library_function_names = [lambda x: x, lambda x: x + x]\n",
    "pde_lib = ps.PDELibrary(\n",
    "    library_functions=library_functions,\n",
    "    function_names=library_function_names,\n",
    "    derivative_order=2,\n",
    "    spatial_grid=x,\n",
    "    is_uniform=True,\n",
    ")\n",
    "\n",
    "print('STLSQ model: ')\n",
    "optimizer = ps.STLSQ(threshold=0.1, alpha=1e-5, normalize_columns=False)\n",
    "model = ps.SINDy(feature_library=pde_lib, optimizer=optimizer)\n",
    "model.fit(np.expand_dims(u, -1), t=dt)\n",
    "model.print()\n",
    "\n",
    "Ut, R, rhs_des = build_linear_system(u, dt, dx, D=3, P=3, time_diff = 'FD', space_diff = 'FD')\n",
    "rhs_des = rhs_des[1:]\n",
    "R = R[:, 1:].real; Ut = Ut.real\n",
    "w = TrainSTRidge(R,Ut,10**-5,0.1)\n",
    "print(\"PDE derived using STRidge\")\n",
    "print_pde(w, rhs_des)\n",
    "\n",
    "print()\n",
    "\n",
    "print(\"----- Add 30% noise -----\")\n",
    "np.random.seed(0)\n",
    "noise_lv = float(30)\n",
    "print(noise_lv)\n",
    "noise = 0.01*np.abs(noise_lv)*(u.std())*np.random.randn(u.shape[0],u.shape[1])\n",
    "un = u + noise\n",
    "np.mean(noise), np.std(noise), np.max(noise), np.min(noise)\n",
    "\n",
    "Utn, Rn, rhs_des = build_linear_system(un, dt, dx, D=3, P=3, time_diff = 'poly',\n",
    "                                       deg_x = 4, deg_t = 4, \n",
    "                                       width_x = 10, width_t = 10)\n",
    "rhs_des = rhs_des[1:]\n",
    "Rn = Rn[:, 1:].real; Utn = Utn.real\n",
    "# Solve with STRidge\n",
    "w = TrainSTRidge(Rn,Utn,10**-5,1)\n",
    "print(\"PDE derived using STRidge\")\n",
    "print_pde(w, rhs_des)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa4d5361",
   "metadata": {},
   "source": [
    "#### Use the precomputed denoised variable given by the K-SVD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0c7322fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "denoised = True\n",
    "if denoised: un = np.load(\"./denoised_files/burgers_dictlearn_denoised30_components32_8x8.npy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "737b71e1",
   "metadata": {},
   "source": [
    "#### Applying best-subset regression on the weak formulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8228253c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['x0',\n",
       " 'x0x0',\n",
       " 'x0_{x}',\n",
       " 'x0_{xx}',\n",
       " 'x0x0_{x}',\n",
       " 'x0x0x0_{x}',\n",
       " 'x0x0_{xx}',\n",
       " 'x0x0x0_{xx}']"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "poly_deg = None; diff_order = 2\n",
    "differentiation_method = ps.FiniteDifference; differentiation_kwargs = {}\n",
    "kwargs = {'fit_intercept':False, 'copy_X':True, 'normalize_columns':False}\n",
    "\n",
    "weak_pde_lib = ps.WeakPDELibrary(library_functions=[lambda x: x, lambda x: x * x], \n",
    "                                 function_names=[lambda x: x, lambda x: x + x], \n",
    "                                 derivative_order=diff_order, p=diff_order, \n",
    "                                 spatiotemporal_grid=XT, \n",
    "                                 include_bias=False, is_uniform=True, K=10000, # new random K points in every calls to the ps.WeakPDELibrary\n",
    "                                 differentiation_method=differentiation_method, \n",
    "                                 differentiation_kwargs=differentiation_kwargs, \n",
    "                                 cache=True\n",
    "                                )\n",
    "\n",
    "X_pre, y_pre, fns = ps_features(un, t, weak_pde_lib, kwargs)\n",
    "# y_pre = y_pre - y_pre.mean(axis=0)\n",
    "fns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2c750119",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_subsets = [tuple(np.where(abs(brute_force(X_pre, y_pre, support_size=i).flatten())>0)[0]) \n",
    "                for i in range(1, X_pre.shape[1]+1)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20140d21",
   "metadata": {},
   "source": [
    "### F-test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6d3cee96",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F-statistic: 11.349618489137109\n",
      "P-value: 0.0007574864881534299\n",
      "The additional feature is significant in predicting the target variable.\n"
     ]
    }
   ],
   "source": [
    "# Is u_{x} significant?\n",
    "alpha = 0.05\n",
    "model0 = sm.OLS(y_pre, X_pre[:, best_subsets[1]], hasconst=False).fit()\n",
    "model1 = sm.OLS(y_pre, X_pre[:, best_subsets[2]], hasconst=False).fit()\n",
    "f_test_res = sm.stats.anova_lm(model0, model1).iloc[-1]\n",
    "f_statistic, p_value = f_test_res['F'], f_test_res['Pr(>F)']\n",
    "print(f\"F-statistic: {f_statistic}\")\n",
    "print(f\"P-value: {p_value}\")\n",
    "if p_value < alpha:\n",
    "    print(\"The additional feature is significant in predicting the target variable.\")\n",
    "else:\n",
    "    print(\"The additional feature is not significant in predicting the target variable.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5db19a5a",
   "metadata": {},
   "source": [
    "#### Mutual information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9c0cafed",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.97556593, 0.94677428, 0.87256523, 0.6650905 , 0.6494026 ,\n",
       "        0.62173837, 0.48995862, 0.43541525]),\n",
       " array(['x0x0x0_{x}', 'x0x0_{x}', 'x0_{x}', 'x0x0', 'x0', 'x0x0x0_{xx}',\n",
       "        'x0x0_{xx}', 'x0_{xx}'], dtype='<U11'))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mutual_info = fs.mutual_info_regression(X_pre, y_pre)\n",
    "sorting_indices = np.argsort(mutual_info)[::-1]\n",
    "mutual_info[sorting_indices], np.array(fns)[sorting_indices]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "95a35c8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "shaps = np.array([0.00107794, 0.00099889, 0.00715736, 0.00148111, 0.00082081, 0.00439339, 0.00050938, 0.00097005], dtype=np.float32)\n",
    "fns = ['$u$','$u^{2}$','$u_{x}$','$u_{xx}$','$uu_{x}$','$u^{2}u_{x}$','$uu_{xx}$','$u^{2}u_{xx}$']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "6a9a50bc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPEAAAC0CAYAAABFXsrWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAPxUlEQVR4nO3du3Mcx7kF8HNcDlR1y8QSvn8AuLRzkJRcClUCIAfMxIdcDhxRr1wSH4mAiCZLuUlCiW5wq0iKzuwqGUsJsfkAIzuQCSG0A0JL3VJm13eD7lkMFtOzs8DMbPfu+VWxuDvs6Wks90PP9JNmBhFJ108mXQARORoFsUjifjrpAozj7Nmz9uOPP+KVV17B66+/vu/fdnZ2sLCwUMuxOvOapbKlVt5YyjFueXd2dvDFF1/82czOAgDMLJk/n376qZmZZX/n1Xms6fyntWxFx2IpR2plKzqWfw9g1XxcJHk7/cYbb1Q6dthzm87/KNdoOv82rhHLzxDzZ1SULnieTaBGPeyfot9WZb755pux0o+r6fzbuEbq+bdxjRh/BuRqYlpCXUyrq6u2uro66WKITBzJNTNbBdQ6LZK8pFqnQ372u/8dK/3//c9vGyqJSPtUE4skLqkg3tnZwerqKjY3NyddFJGJ8d//hex9UrfTCwsLUMOWzDrf1bSTvU8qiCdFz9wSs6Rup0XkIAWxSOIUxCKJKwxiknMkb5H8yr8/RXKx1ZKJSCWhhq13AVwG0AUAM9si+WZrpRKRykK301tm9hJAVAOr1U8sUr2fuEtyGcALkvMAVgA8B/B1w+UrpX5ikYP9xIU1sZmtA7gH4L8BXASwYWafN188ERlXYU1M8hSA98zsw+w9yQUz22mzcCIyWlkX05XshZltwTdyiUhcQs/EJ3zg5nVGZUbyPIC+T7ttZk8L0nQAvAegb2Z3csdvA7jt3y6b2c1R15sW4w7rBDS0U/aEgpgk/wBgw79fAfCkLCMfnCtm9r5/v+HPG/ZqIIsugPsAelkeIjJaYRCb2QOS2wDe8YfuFNTMwy7CtWBn+iRPD9fGZtYjWXRrfgPAYwDz1YouIkDJLCYftIPAJXnMzH4oyasDdyud2cV4AdkFsA2gQ/KGmV0eTpD1EwOumb3qqoEi02BzczM/RmIhexEM4oIRWhcAfDjmdTtVEw49H6+T7JrZdj6N+ollluUrrrW1tZ3seKiL6RZcrfoid/jkiGv0sT9o5+Fq1pF8g1g315i1W+U8EQnXxPfN7GH+AMkvR+R1D+65NtMpap0O2Mb+gJ8froVFpFion/gFyUWSx3LHzpVlZGZ9ABskl/2QzUFAk3ziW6/h/20FwAX/Gj7YuyTPk/wE7tZdRCoI1cRfA3gE19VkAAjgBIDPyjIzs8La2szO5F73APSqnisi5UJBfKHgdvpUC+URkTGFJkAMB/AigONtFEhExhNqnT4G4BqAuewQ3NziiU5FzPqJ1Ucss6zqfOJr2BtyGU13j/qJRSrOJ4abP/zQ31Yf96O3dDstEqHgVEQ/AQIA3iL5EdTtIxKlYMNWtiCAmV0B8BLAzEwNFElJpW1c/HI9IhKh0LrTi7nXJ0guaclakTiFnokH833N7DvfwKXleUQiNAhiv+vDuyQfA7hB8lHuz1eIoKtJ606LlPQT+8Xi1wGsk1waHrUVA/UTi1TvJ37UQllEpAahIP6c5EKbBRGRwwkF8W0MLa1D8lLjpRGRsYX6iW8DeE7y53BL9BDAGQDaykUkMqEgvmxmD/IHSC61UB4RGVPZutPHACzDLXjXi7G1WkTCI7aW4MZKn4SbvXQzhhFb6icWqT6fuGNmH+QPkCxdKK8N6icWOdhPHAri7ysek0SMu2mbNmxLRyiIj/s5xNm60Wfggniiy/OIyEFlDVun4LYgBdxi8gpgCVJNPzllC+VdhWudBoB5ko9HbKgmIhMQGrF1DcB1M5s3s3kAV+CCWkQiU7ZQ3mBbUzP7DgW7NojI5IUatiy/H7G/vc7WoAbJS2bW+hBMrTs92/Tc7VTtJ/4Sbux09p4AQPIq9vZlaj2I1U8sUr2f+MBeTHkaRy0Sj0p7MQFAfn6xxlGLxCPUxbQA4AO4/ZcAdwu9BOC1doolIlWFbqevALiP/YvjdUZlRvI8gL5Pu+03Dx9O04EbRNI3szvjnCsiBwUXBch3MQEAyU/KMvLBuWJm7/v3GwBWCpK+eoRzRWRIqJ/4e5IfkXyb5Jt+GuKNEXldBPA8975P8vRwIjPrwdW4Y58rIgeV3U734buWvJMj8upgf3Duwi0oUEWlc7N+YgDqK5aZs7m5mZ9Lv5C9CAXx/eEWaJJfHuK6nUOcEzxX/cQyy/IV19ra2k52PHQ7/YLkoh+plRm1KEAf+wNvHsB2xfId5VyRmRaqib+GW0CeJA17o7Q+K8nrHvY/N3fGaGE+yrkiM63yiC0/vzjIzPokN0hm0xcHQUnyCYAln2YZruW5Q3LbzHpl54pIudCiAEUjss4A2Co4nj+v8LnZzM7kXvdQMCMqdK6IlNsXxCS/hQvWp3BdPoQbtTWxSQ8iUm5fEJvZLwGA5IWCwR6lt9Nt0FREkYpTEYcDOHSsbepiEqm+tamIJCLUOi0yc1JdOUQ1sUjiQnsxzZG8RfIr//4UycVWSyYilYRup98FcBlAF3CNWjFsqCYiB4Vup7fM7CX2VvYQkUiFgrhL8jqAZT+f+Dp8rTxJ2tpUpHo/8bof3PEOgAuIZC8m9ROLVFyyluSiH9yxRfIEXM38ZgyBLCL7hRq2ugCeAYMtXL4jeQna2lTkSJroix4EMck5uLWu3gcw53d7yOwCuD3W1UWkFYMg9q3R6wDWSS5pgXiRNFTaAcIv1aN+YpEIlW0yfg17OyFm84r1TCwSmVDD1jUAG/71biBN6zSfWKT61qYb2S111rUUw+20+olFxphPTPIP/uVbJD+CG/QhIpEJNmyZ2Yf+9RUALwHcbLNgIlJNpUUBzGy96YKIyOHkB3ssjkh71czeabY4IjKufE38OVyLNOGGXe5ib5OzDvZvriYikcgH8QU/ThpFI7ZiWLJWRA4aNGxlAezNFaQ93nxxymk+sUj1fmL6hQAe+fevwe0IMdERW+onFqk4n9jMHpDchlsUAADuxrB4vIgcFOxiyhYFyN6TvGRm2otJJDKjNlQb/BO0oZpIlGrdUI3kebhuqQ6A7aKNwkNpSN7G3sIDy2amEWIiFdS2oRrJDoAVM3vfv9+A20y8apougPsAetm/i8hodY7Yuoj9t+B9kqeHauOyNDcAPAYwX6nkIgKg3hFbnVx6+POHA7IsTRfANoAOyRtmdnn4Alk/MQDNKZaZs7m5mR8jsZC9aHrEVqdqGjO7k7vWOsmumW3nE6qfWGZZvuJaW1vbyY7XOWKrj/1BOw9Xs45MQ/I8yU9yx6NZTUQkdnWO2LoH91yb6RS0ThemIQnsD/j54VpYRIrVNmLLzPokN0gu+0ODYCX5BMBSKI0P5PMku3DPxlpFRKSicUZsHTOzH8oyM7MvA8fPVEhTeFxEygWDuGBhvAsAPmy2OCIyrtC607fgGqFe5A6fbKNAIjKeUE18v6CLaeK3u1p3WqT6fOIXfgTXdu45+ByAzxos20jqJxapOJ8YrivpEVxXk2FvFtNEg1hEDgoF8QWtsSWShsq7IiKCNbZE5CDtiiiSuKR2RRSRg0Ibqm34/ZgeAjjuR2/pdlokQkntiqh1p0Uq9hP7Gvihf32F5LuIYFdE9ROLlOxPTPJjko9IXie54I/N+THUj4bmG4tIJPI1cQ9ukbrBzCUzewngax/Mi2b2rO0Ciki5/DPxidCcYR/MWsBOJELBhq0CnaYKISKHlw/iX41Iq6mIIhHKB/FdkndJ/iyfgOQxknexN/hDRCKSX+1yC8A6gC2S3/qW6n8AeALgTgyNWuonFhnRT2xmPQC/IHkCwGkAT2PqWlI/sUj1/Ym/AxBN8IpI2Dit0yISIQWxSOIUxCKJUxCLJC6pIFYXk0j1JWujpC4mkZKpiCKSJgWxSOIUxCKJqzWI/R7Dy/7v0+OkqXLuuP79z7/Vkc3E8m/jGqnn38Y1Yv8Zagtikh0AK2bW83sN36iapsq5h/Gff/29jmwmln8b10g9/zauEfvPUGdNfBHA89z7fkGNGkpT5dyBoi6mqr/JitIN59d0/ke5RtP5t3GNWH6GmD+jonShrlWaWaUMRyH5CYC+md3x72/DbZHaG5UGbsZU6bn++J8A/BfcKiPPkGtmh+s3y78/yrE685qlshUdi6UcqZWt6NjrAF7xr380s7NA8/3EnSOkOXA8K7SI7KnzdrqP/YE3D2C7Ypoq54pIgTqD+B72r8PVMbOnFdNUOVdECtT2TAy4biK4WhXAYKUQkHwCYMnM+iVpCo/Hzm/7iqaWL2o6/zaQPNH0CjGp/z8cJf9agzhGLXz45+Ba1vsAumZW6/avdeRPcgnAZQBPzOyqX37pPTO7WmdZc9e7BeAfcJsRPPNb5S6b2R+buJ6/ZvT/D03ln9QEiEM6CeC535qm9g8f7vn9HQB/9e9jzP+0mb1F8hTJ6wB+D6BbVwELXAbwKoDfkLwJt6NmD0BjQYw0/h8ayX8WgrjpD3/XzC6SPAW3EXvd6si/R/KYX9F0y//W79RWwiF+x5DBpnzA4G6gSSn8PzSS/ywEcdMffs/n3YWrcZ5FmP823C+ydQAwswck+3UVcBjJj+Fq4g24W+odAHNNXc/r+UenJv8fosx/Fp6J5+A+mC7chumf15x/o402LTUK1XoNX+s+BrAM98tjDm7wTq2f/dA1F4E4G56azn8WgrjuL2ijjTZtNAq18DPMwfVGNPkMPHzNaBuems5/6qYikrxF8qPsNxuAFyTfrvESlwFswTXa/AXuue+1hPJv/Bpm9rLNAPbmAVwDcB5uGO/M5D91NbGvBV4FsAL3YRyHq3Ea6U7x11wys4ejU8aZf1vXaBLJc/5Z/xQAq/u2N+b8p65hq+mW0aYbbdpoFJpQw1PTUmhgbCT/qQviFr6gTwHcgWu0uelr/vsJ5d/WNdo2n3WhzVr+03g73WjLaNONNm00Ck2i4aluqTcw1pn/NAZxq19Qkgu+tpcWNd32kVL+yQfxLIwLlmpSb2A8bP7T0MV02szeAnDPjwveRfPjgpvuApIRSH5M8i7JS35cPFBzA2Mq+U9Dw9YsjAuWg1JvYKwt/2m4nZ4DcNHM1nPHGrvtKWr9Jvm2bqfblXoDY535J3877UcHrQ8dy9eSJ2q+5FMA7wH4Hu436Fdwo22kQSSXSP7FPzIB7jNv7DFmeNRZ7pa3kfyPIvkgHtbCsMvHcL9BH5jZRTP7dZMD+2Wg1baPgu/Rbs3fo9okfzs9bBLDLqV5fjTTczP7wb8/B9cL8euGrpfM92jqgrhI6uOCpf22j0AZovweTePtdKNdAzIZbbd9pPQ9mroghhqeZkILbR/JfI+m7nZ6GsYFy2gtDYtM4ns0dUEssyvWZ9amTcOILZlBUzon+lCm8ZlYZkMyz6xN0+20JCmlZ9amKYhFEqfbaZHEKYhFEvf/r0rgWJe2c0wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 252x189 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "with plt.style.context('science'):\n",
    "    mi = pd.Series(mutual_info)\n",
    "    mi = mi/mi.sum()\n",
    "    mi.index = fns\n",
    "    mi.sort_values(ascending=False, inplace=True)\n",
    "    # (?) Mutual information for each candidate\n",
    "    # plt.title('Ranking by mutual information', fontsize=12)\n",
    "    mi.plot.bar(fontsize=12)\n",
    "    plt.ylabel(\"Candidate importance\", fontsize=12)\n",
    "    plt.savefig(\"mutual_info.pdf\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "136bf798",
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.load(\"X_weak_burgers.npy\")\n",
    "# np.load(\"y_weak_burgers.npy\")\n",
    "# np.load(\"best_subsets_burgers.npy\", allow_pickle=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "73b60131",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Exact explainer: 10001it [00:33, 226.01it/s]                                                                                                    \n"
     ]
    }
   ],
   "source": [
    "linear_model = SkLinearRegression(fit_intercept=False).fit(X_pre, y_pre)\n",
    "explainer = shap.Explainer(linear_model.predict, X_pre)\n",
    "shap_values = explainer(X_pre)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "62de53f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "candidate_importance = np.abs(shap_values.values).mean(axis=0)\n",
    "candidate_importance = candidate_importance/candidate_importance.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "2693482c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOsAAAC0CAYAAAB10tv3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQrElEQVR4nO2dP2wcxxXGvy8IEAeByBOTIiVJJU2aUJQduDRMWinURaIcuFBly1YfR38QwHQTx4Z7S6YbuxBg0UqqBHB4kllHMqkgQFIkplkFSSGasuEiQIKXYmbF5XH/zN7t7N3efT+A0O1wbmaWq7cz8+bNNzQzCCFGn28MuwFCiDBkrEK0hG8OuwFZnDlzxr7++ms88cQTePrppw/9bnd3F7Ozs9HSYpc/rm1rW3tHpR1FaQDw/vvv/8HMzgAAzGzkfl577TUzM0v+TRM7bRh1jkPbstJGpR1ta1s6DcCqebsY6WHwM888E5Q2yHd702KXP0gdscuvu45h/I1C00blb1TpuzaEnrPsJ+tNU8Qnn3xSKX9VYpffRB26h+GX308dSPWstBFculldXbXV1dVhN0OIoUPydTNbBUbUG7y7u4vV1VVsbm4OuylCDA3//382uR5Jb/Ds7CyyetZjF25WKuerD16oqUVCNI+ft+4m1yPZswohjiJjFaIlyFiFaAkjaaxyMAnRcgeTEJOEHExCtBQZqxAtIdNYSU6TvE7yY399kuRCoy0TQhwir2d9CcBlAFcAwMy2Acw01Sg5mIQIdzBtm9kjkkMJHJaDSYijDqY8Y50nuQzgIckZAM8B+AzA3aLCSZ4DsA+gA2DHzLYy8twAcMNfLpvZW1VuQIhJJdNYzWyN5EkAzwM4D+CWmZUZagfAc2b2sr/egDPyXuYBrAPoJnmFEOVkGqs31Itmdim5JjlrZrsFZZ2H630T9kkuZvSubwK4jwbnwEKMA0VBEVeSD2a2TfJZpMbPGXTghsAJe8g2yHkAOwA6JN80s8u9GRIHE+DG7aE77oUYBzY3N9PO1dnkQ56xznkPcJpOH/Ue+Y6ZvZt8JrlGct7MdtJ55GASk0y6g3r99dd3k/Q8YyXJdwBs+OvnAHxaUsc+DhvnDFwPmi70HID5lFNpr7TlQggA+Q6m2yR34BxMAPBuRk/byy24+WhCJ2O+uoPDBjzT26sKIbLJnbN643xsoCSnzOzLgvz7JDf8kg+QMlySnwJYMrMtkudIzsPNXVcGvgMhJoRcY/UOpTQrAC4VFWZmH+WknyrLkyZxMMm5JCaZoAgmktfh5qAPU8kn4jXrMHIwCREewbRuZnfSCSRLe0QhRDzyAvkfklwgOZVKO9tEg4QQ2eT1rHcB3INbwjEABDAH4O0mGqU5qxDhu25WMobBJ6O1qgfNWYUIlHXJMNQFAMfjNUsIUUaeN3gKwDUA00kSAEPJFjkhRDzyhsHXcBBqqJBAIUaAPG/whpnd8cPh4z6aqbFhsGRdhKigG0zyHb+f9TTJRbigiEaGwXIwCVHBwZRsPDezKwAeAZD8ihBDJEiR38zWYjdECFFMnm7wQurzHMmljMB+IUSD5DmY5pMPZva5dzTN5+StHTmYhChwMJGchhM9exnANMmrqe/t4UA+NDpyMAlRsOvGzB4BWAOwRnKpN4ophBDd4N68ZtatWo8Qk0jeMPhe1YJSusFdv8H8zZK80gwWogJ5xvoeydmKZWXqBufkfRIHEVJCiADylm5uoEdGlOSLZvZeQVkdBOgGewO+DyDPkKUbLCaaqrrBNwB8RvK7cNIuBHAKQJGxZtHJSvTiarlfkoNJTDJVdYMvm9ntdALJpZI69lGuG3wRwJ5XN3wKwHdJ7kiOVIhyinSDpwAswxldN8A7XKob3KPG/xTchgEZqhAB5EUwLcHFAp+A223zVlkEk5ntA9gguey1gw/pBnsPcHK9CDdnXfG9rBCihLxhcMfMXkknkCwVTAvRDfbXW8g+DhKANJiEAMK3yH0RmBYFOZiECNcNPk7yFwCSOecpOGOVrIsQQ6LIwXQSwEWftF528rkQIi5FgmlX4bzBADBD8n7RwVRCiLjkhRteA/CGmc2Y2QzcKehXc/LWjrbICRHuYNpIn8dqZp+TbGx3jBxMQgRqMAGw9Dk3/vN06vrFOM0TQuSR17N+BBcbnFwTAPyG9OTcm6pxwkKIAQg+6yZNQJywEKJmgs66AYD0/tZ+VCSqIAeTEOEnn88CeAXufBvADX2X4HbKREcOJiHCI5iuAFjH4XNuOpHaJIQIIHfzeXrpBgBI/rKB9gghcsgN5PexwTs4kGpZAXCpiUYJIY6St856BcD34PaznvI/J5pqlBxMQoRHMK1nnH6euVc1BnIwCRHuYHroz7vZSQXvnwXwdlHhISLfPs8OnBzpIakXIUQ+ecZ6F07omyQNB1FLucaaEvl+2V9voEcNwue56pUjtkh+AUDGKkQAwRFMfn9rEZki3+ne1es0nfLlzcPpBwshAsjbfJ4VoXQKwHZGekIHASLfAOAF1RbhPMxHkMi3mGSCRL5J/h3OKLfgeknCRTH1G7zfyUo0s67fJHARGSeqy8EkJpkgkW8z+yEAkFzJCIooGwbvo0Tku6euLsl1kt2i0+aEEI68QP4jw92stB5u4fBa7BGRb5IXSaaFwNPhjEKIAvIcTJXx59ds+Pko0CPyDbcR4BaAJ1Nz1hvqVYUIozZjBYJFvrs9/x5BIt9ChEcwDRU5mIQI1GAiOU3yOsmP/fVJH9EkhBgSeT3rSwAuA5gHnHOp7GAqIURc8nbdbJvZIxwoRQghhkyesc6TfAPAMsln/efGjmbUFjkhAh1MZrbmgyCehwsJbPSsGzmYhAjcIkdywQdBbJOcg+tpn9XhVEIMjzwH0zyAB4A7OgPA516Ff6yN9diFm5Xyf/XBC5FaIsRRHs9Z/XLNSyTvA3iT5L3Uz8doMDRQc1YhCuas3vu7BmCN5FJsIe8iNGcVIjAoImPj+YLWWYUYLkWHKV/Dwclxyb7WsZ6zCjHK5DmYrgHY8J+1jU2IESAvKGLDzO744fBxv4xzvKlGycEkRIVdNyTfMbNLAE6TXITbWN7IMFgOJiECgyJ8j3rHf75C8iVkaCX1UkE3eAZO62ndzHL3tQohDgjaz2pma2V5AnWDF+GM+CN//QUaHF4L0WYeG2vAftWrZvZ8we9LdYPhetQVAC/7672MPEKIDNI963twHmDChRvu4UAHuOPTi+igRDfYD3nTw96ZLEOVbrCYZEJ0g1d8HDCyIpgCpEiz6OT9wqscZop8y8EkJplS3eDEUD3TOErZ3HIfgbrB3sn0oYa/QoST52Ci33B+z18/BTcfLVq6uYWU/CgydIN9wYsAtsxsx593AzPLFQMXQjjylm5uk9yB23wOuF6wUOQ7UDd4Hm5JaM8fnzFjZvIGCxFA7tJNsvk8uSb5opkVnnUToBu8hYClGukGC1ESwZRxMNXjX6G/g6n6Qg4mIUoimAY8mEoIEZE6D6YSQkSkzggmIURE6oxgqg05mIQodjDFiGDqCzmYhCjQYKohgkkIEZE6I5iEEBGpLYJJCBGXKhFMU2b2ZRONkoNJiGoaTL06wSsALsVoVC9yMAkRfjDVdbhlm4ep5BPxmiWEKCOvZ13PWLrJDNIXQjRDnm7wQ39kxlQq7WwTDRJCZJPXs96FW7YhScPBrpu3m2iUHExChDuYVhTBJMRwqSLy/Rgf5F8awRQo8t0BcBHAvpm9G9ZsIURtp8iFiHx7nhygvUJMLHWeIhci8g0z6yZCaUKIcPKMdSMZCpN81szuBhym3EGJyHcoEvkWk0yIyPchajpFrlMxv2udHExigskT+c6TdbnjDRVmdgXAI5SfIrePQJFvIUR1HhsryVdJ3iP5BslZnzbth7/3eva7ZnELh0MSM0W+hRD9kR4GdwF001vhzOwRgLveaBfM7EFeQSEi3z7PMpyXuENyZ9LOZz124Wal/F998EKkloi2kTbWOTP7bVYmM3tEstRZFCDynXWS3BEUwSREhS1yGXTqbUo+cjAJUaDBBOAnJd/VFjkhhkjaWD8k+SHJY+kMJKdIfoiDIAkhxBBIn8+6TXINwLbfabMPFw9sAF4pci4JIeLTe9ZNF8APSM4BSM5RLVuyqR05mPpH3ubxIcjB5A20cSNNkINJiGIHkxBihJGxCtESqqyzNobmrKNL1TkxoHlxvwwSFNEYmrMKoTmrEK1FxipES5CxCtESRtJYEwdTStpCiIlDDiYhWkJUBxPJcySX/b+L/eapyn//9dc6ihla+U3UMQ73EHuk1cRIbpA6autZQ3SDK2gLV+J///4bvvn9Hw1azNDKb6KONt1D3lruf/58G9/68T+PpNe1jru5uRl9XX+QOursWTN1g/vI85ist1Do2zsrX1Z5vWmxyx+kjtjl113HMP5GIXUcu3AT3z79Kxy7cPPQz69/95cjad8+/au+7yG0Fw39Ls0sqMAySP4SqSMxSN6AOzqyWyWPT/89gO/AqVM8QGrcDjfhTl/XnRa7/HFtW1baqLSjjW1bgNum+rWZnQHiO5g6/eRJGieEOKDOYfA+ynWDQ/IIITKo01hDdIOlLSxEn9Q2ZwUOHfkI4LHyRK9ucGaeGupe8OU9qKO8pstvApJzMZU/mvgbTfJzrnWd1cw+MrNu8pNKP2Vm+0V5auAEAJCcDThEayjlk1wi+UeSb/jrueRzDEheJ/mL5D8ggIckfxarPsR/Bk3UEf0eSC6knkkwIxnB1CczAJ4H8Cd/XfUQrSbKXzSz0yRPeiP9DYCYx19ehjsP9+ck34ITwOsCyBRzr4HYz6CJOpq4hxMAPvPH1MybWVAd42Sse2Z2nuRJOEXGUSy/S3LKH1GyTfIsIoqn++NP7vgfAK53j1Uf4j+DJupo4h76eiGMk7F2/R94Hq4HeRCh/IUBy9+Be0hrAGBmt0nu19S+I5B8Fa5n3YA7x2gXB6fZxyD2M0jqWIhYR+zygT5fCLU6mIbJmDhPar0H34veB7AM95KYhgtCea+uOnrqi/oMfB0LQLsdTCSn4V4G8wCOhz6P1horyesA/gHXYzwgOQVgOe9wrRrqOwsXKrmPCvOMkjKj3oP/T7EU8W/S6DPwddb+HJos39fR10ttJPezBnIZwDac8+SPcPOypyLWNwPgGoBzcALodRD1HszsUUzDQfPPAIjzHKKWX5dXvrU9axYkl8zsTnnOvso+6+eYJwFYxGFYtHtogtjtj/0cYpTvRzhPwu0wW4T3ypvZ1SrltNbBNA7OkyHcQ60Mqf1NOBJrLb8ur3xrjRXAFoB34Zwnb/m313rE+maSJZcay2z6HupmGO2P8Ryill/XS621w+BxcJ7EvofYNNH+BpxwTTznWrzyrTXWXkjO+jdWXeXVMs8QgxH7OTTxnOt6qbXWWIe0bNBq58+40IATa6DyfU96GcCnZnaV7gjVi4O+ALR0kwPJV+lOgn/Rx3ACLXL+jAuxn0Ok8hfN7DSAWz4GfA81xIC31sHUQNxr250/40Ls5xCj/Cgx4G0eBh/xsJH8WVuif0QYDTgSay/fl3nezNZSaQMP3ds8DN4CcBHAF3BvxI/hok9qoYHoH5EBe/b8wj3TaFFRvc85NRQetMy1nrT0CHCun3LbbKz34d6It83svJn9NFaAumiUKPO9PDJCAff6CQWsWIfCDUX78dFDn5nZl/76LJwn9aeR6mtq6WbgOmSsYqSINd+r2Ibo9fVTR5uHwWIMiTXfy6OJJbq66pCxipGmrvleAVEdlXXWoWGwGGkaCjeMHd882eGGYnKZ1LDP1kYwicmg7Xt+60RzVjHqNDGnbAUaBouRRmGfB8hYhWgJGgYL0RJkrEK0hP8DGc5q8MutvcoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 252x189 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "with plt.style.context('science'):\n",
    "    mi = pd.Series(candidate_importance)\n",
    "    mi = mi/mi.sum()\n",
    "    mi.index = fns\n",
    "    mi.sort_values(ascending=False, inplace=True)\n",
    "    # (?) Mutual information for each candidate\n",
    "    # plt.title('Ranking by SHAP value', fontsize=12)\n",
    "    mi.plot.bar(fontsize=12)\n",
    "    plt.ylabel(\"Candidate importance\", fontsize=12)\n",
    "    plt.savefig(\"shap.pdf\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "12270b74",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([39597.08633088, 33249.0402039 , 26044.00310333,   644.04920612,\n",
       "          605.75237157,   543.5305861 ,    39.8190605 ,    39.77767518]),\n",
       " array(['$uu_{x}$', '$u_{x}$', '$u^{2}u_{x}$', '$u^{2}$', '$u_{xx}$',\n",
       "        '$u$', '$u^{2}u_{xx}$', '$uu_{xx}$'], dtype='<U13'))"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f_scores, _ = fs.f_regression(X_pre, y_pre)\n",
    "sorting_indices = np.argsort(f_scores)[::-1]\n",
    "f_scores[sorting_indices], np.array(fns)[sorting_indices]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "d5a15cfb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOsAAAC0CAYAAAB10tv3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQTElEQVR4nO3dMXAbxxkF4PcynkmRMQkzKVK4IKmkSROKslJrTMku0kWinHGRcUFTdh9LlhqLlRONe0uiU8RFZiRa6VLYhGTWsUUqTVIkolhlkkIU5EyKzCTzp9g98ADcAQdwbw8HvG9GI9zxuLsAuNjF7t7+NDOIyPj7VtUFEJFiVFlFauKFqguQZW1tzV5++WUcHBxgfn6+42dln6siz0koW93KOy7lGHRuY2PjN2a2BgAws2D/AFwAcNb/v1zk2qyfffDBB2ZmlvyfVva5KvKchLJlnRuXctStbOlzAK6brxfBusEkGwDOmVnTzD4D8OsB114alOaZM2cKnTvO73afKzv94+RRdvqh86jiNSp6blxeo6F+1wK0qO4DAOsALqeOt5DTusK1vpcxoGUt6ssvvxzq+mGVnX6MPPQcqk9/lDyQallpgaZuSF4G0DKz2/74FoAtM2t2XbcMYN9X7t3unwPAW2+9ZUkf/syZM4U/oUQmwc7ODnZ2dgAAGxsbvzWzt4DyB5gaWSfNrEUy95fm5+dx/fr1kookMt7SDdTGxsZBcj5kZW2hs3LOwbWgbSTXARySXARwGsB3Se6bWcd1eV78xe+GKtC/Pn1zqOtFxlnIynoXnYNKDTPbTV+QdJEBgORpANtFK6rItAtWWX3XdpvkWX+qXXFJPgSwYmYtf7wMYBlAY5iWVWSaBf3Oam7KJuv8qa7jXQDnQuYtMum03FCkJsaysh4cHOD69evt4WuRaeT//ueT47FcG6ypG5H2KqaD5HgsW1YR6ZVZWUnOkrxJ8nN/fJLkUtSSiUiHvG7w2wCuAFgEADPbI/lqtFKJSI+8bvCemT0HoD1fRMZEXsu66Bc3PCU5Bzcn+hjAg2glE5EOmS2rmW3CLR/8HoCLcMsCP4lVKE3diBScuiF5EsC6mb2bHJOcN7OD8ouoqRsRYLipm/eTB2a2Bz/YJCLVyKusC36AKa1RcllEpI+8ASaS/BjAtj8+B+BhnCKJSJbMympm90juA3jDn7rtu8IiUpHctcG+crYrKMkZM/smSqlEpEfuABPJV9P/0Gdr0dA0dSNSfOrmJtyeSk9Tp0+UV6xOmroR6Z26yesGb5nZ/fQJkpm7QIhIHHmV9am/y2Y/9T31PICPopSqIto9UcZZXmV9AOAruCkcA0AAC5jwyioyzvIq62pGN/hkhPKISI68hfzdFXUJwEsxCiQi2fJGg2cAXAMwm5yCu7c1yi1yydSN4tzINCu6Ydo1HC01PCy1RBk0dSNS/K6bbTO777vDL/nVTOoGi1So3wqmj/3D10j+EsBqnCKJSJa8hfz3Adz3j98n+TaAG4MSI3kBR9Hk9rsDU6Wu2Qfwik//dvc1ItKr0CbffpuXvkg2AJwzs0v+eBtd8Wz8NVd97Jtdks8AqLKKFJC3b/BS6vECyZUCW5FehNtULdHy0eLazKyVBKnyMVq/HqnUIlMod3dDAI8AwMyeAHhCcg39p24acF3gxCFcQOUefufEZeR8D06mbgBo+kamzs7OTvqOs/nkQbuykpyFax0vAZgleTX1+4cAbo2QbyPrpJk1SQLAOjK+C2vqRqZZuoHa2Ng4SM63K6vfc2kTwCbJle5VTAW00Fk55+AGkjL5CrtFspk1ECUinfKmbr4aIa276LzntdFdCUmuk0zfxB59wYVIXeVV1k9Izg+TkJm1AGyTPOu/k7YrJcmHfiT4buqaywBuqVUVKSZvgOkWur5vklwbtCu/mWXeoJ6MAHvNrv9FpIB+lfUxye/Cbe1CAKcARAmhoYX8IsUX8l8xs3vpEyRXSitVF40GixTcg8nvGzwD4CzcqG5zhNFhEQkobwXTCtz85wm4u21uKJiySLXyusENM3snfYLk+QjlEZEceVM3zwqeE5FI8lrWl/w9rMkc6Cm4yqrI5yIVydsw7R7c/ayr/t9DRT4Xiato+IwZAFfhRoMBYI7k17ECU03y1I02Epeiiu7BdA3Ah2Y2Z2ZzcFHQr+ZcKyIR9NswrR3u0d/TquWBIhXKq6zmu8IA2t3i2dTxWtkFE5FOeaPBn8GtDU6OCQD+hvQk7k20AScRGSLWTVrMdcIi4hSKdQMA6ftby14nrKkbkeJTN/MA3oGLbwO4ru8KgNNlFi4xyVM3IkUVjXz+PoAtdG670iipTCJSQO7N5+mpGwDw27CISEXyKuszvzZ4H0d7Aa8CeDdGoUSkV79ucAt+ysY7kX2piMSQV1m3MqKfZ26GJiJx5K1gekpyKb2KCUC0m881dSNSfMO0B3AbfZOk4WjV0kflFs/R1I1I8ambnhVMJE+WVioRGajwCia43SJEpCIdLSvJv8JVyl24WKuEW8VUaPH+EJHP53w+W2amW+9ECuiorGb2QwAguZqxKKJvN7hg5PNluEr8mT9+BrfVqYgMkNcN3ityrsvAyOdwLeql1PFhxjUikiFvgGkUDQyIfO67vOlu71xWV1mRz2WaDYx8XpJG3g98nNbVrJ9p6kam2cDI5wG0UDDyuR9kuqPYrCLF5cW6mSV5k+Tn/vgkyaUBaQ2MfO7TWgawa2a7JBdJLo5YdpGpkteyvg3gCoBFwA0uDQpMZWYtkts+6jnQFfkc7ub1RbjNww/9/k5zZqbRYJEC8irrnpk990sNCysQ+XwXmqoRGUleZV30LeRTknNw86WPoVg3IpXJm2fdhPsO+j24EdttxboRiavohmlLfhHEHskFuJb2VTOL0rJq6kakeKyb9gitmT3xC/s1aitSoXbLSnIWbsngJQCzfvf9xCGAW5HLJiIp7cpqZs8BbALYJLlS9kbeIjKcQvez+i1e+s6ziki5+gVTvoajyHHJfa2auhGpSN486zUA2/7xYc41pUmmbnTHjUyzohumbSdd4WTKJmY3WFM3IsWnbkDyY//wNb87f+btbCISR+4Ak5m96x+/D+A5gBsxCyYinQrdz+qXH4pIhdKLIpYGXHvVzN4otzgikifdsn4CNwJMuKWFhzjaU6mBziBVIhJZurKumtkTAMhawRRzR35N3Yj0mbpJKqo3i17RbhrX1I1I8Vg3JPkhXHAqADgN3XwuUqnMympm90juA0gGlO4U2ORbREqUO3WT3HyeHJNci7lbhIh0GhSYqv0jFAhMJSLlCRaYSkTKFTIwVTDaME2kz9TNOK1g0tSNSP+pG61gEhljY7mCSUR6BV3B5KPDteBa4v2cwFQNAOsAWmZ2e7jiikyvYCuYfCU8Z2aX/PE2XNiNbq+MXFqRKRZyBdNFdM7Ntkgud7euZtZUmEeR4Q2zgmnGzL7pk1YDRwNSgBugmhulUMnUDQDdeSNTZ2dnJz1tOZ88yK2sGRukrQJ4d8h8G0NeD0BTNzLd0g3UxsbGQXI+b9/gm3Ct5NPU6RNZ16a00Fk55wDsD1lOEcmR17JuZUzdZAZKTrmLVLRzAI2s0WARGU3eVqRPfciMmdS58/0SMrMWgG2SZ30g5nbFJfnQjxbD/+wcgFX/WEQKyGtZH8BN25Ck4eium4/6JWZmma2vmZ1KPW4CaI5UWpEplldZV7WCSWS8FI4ih4h7MOmuG5GCsW6qjiKnqRuR4humVRpFTkR65Y0Gb/t4N/cBvORXM0XrBotIL0WRE6kJRZETqYl2ZSX5HsmvSH5Ict6fm/VrhL/qut9VRCJLt6xNAOtmdtXMDgDAzJ6b2QMATwrs0RSMpm5E+k/dLJjZ77N+ycyekxzpdrdRaOpGpHfqJneAKUMjbFFEZBjpyvqTAdcOukVOREqUrqx3SN4h+WL6ApIzJO/gaJGEiFQgvbvhHslNAHv+TpsW3EIIA/COmT2qpIQiAqA31k0TwA9ILgBYBrBbxZSNIp+LFFzI7ytoZfOqGg0WOd5osIhUSJVVpCZUWUVqQpVVpCZUWUVqYiwrqxbyixScuqmapm5ENHUjUluqrCI1EbSykrzgw2dcILk86jXD+u8//hwimcrSj5FHjO//ZedR9/SPm0ew76xFIp8PER19KP/751/wwvd/dNxkKks/ZB4v/uJ3mef/86d7+PaP/95z/l+fvhkk/Rh5hEo/z87OTulr0Y+TR8iWNTPy+QjXtGV9ChVtgbKuy0qv+1zZ6R8nj7LTD51HFa9R0TyKtnBFy1v03HHyoJkVSnAQkpcBtMzstj++BRc6sjnMNf78HwB8B253ikdIjYjBDWWnj0OfKzv9SS1b1rlxKUcdy7YEd5vqv83sp0D5UzeNUa5JCiciR0J2g1sYHPm8yDUikiFkZb2Lzn2asiKfF7lGRDIE+84KuGkZuNYTQHvnCZB8CGDFzFp51wTIe8mn9yhEerHTnwQxXiOSC2XuXhLpOYyUR9B5VjP7zMyayb/U+VNm1up3TQAnAIDkvI8iEFrZ6YPk0nE2Uye5QvILkh/644XkcSTBXyOSN0n+MvW6PCX5sxBp5yj9fR41j7FcGzyiOQBvAPijPw4dS7bs9AH3Jj724UsWfTSEYSyb2WskT/pK+isAi6EL2UcZr9EVAK8A+DnJG3Cb+DUBZG5IH0CM93mkPCapsh6a2UWSJ+F2ZKxb+sDx/1CaJGd8iM49kucRd3P24K+RmT0HcN//A+B6ECHSzhHjfR4pj0mqrE3/5BfhPn0f1Sx94Ph/KPtwlX0TAMzsHslWwPIN0vTd1WCvEcn34FrWbQBNH4dp9rjp9hH8OYTKI+gAU5UiDDyUmr7PYxbuDVyEC2L9SeD0azc441vRrwGchfsgmoVbSBP0tUnltwSM5wBTbSsryZsA/gb3afuI5AyAs3nBtcYt/Zw8g1am2M/Bd7sfw432j/KdOyvNWbiZhNJe9678gj+HUHnU+Ra5KwD24AYevoD7TnO6RunHGOks/Tl0mQNwDcAFuE3ij82HHY1SUb3gzyFUHrVtWbOQXDGz+4OvHI/0favxCtydR8vwI51mdjVUHhl5lvYakTzvvyefBGB1nJOO8RxGzaO2A0xlDzzEGNgoe6SzosGZsgfhyhbjOYyUR20rK4BdALfhBh5u+FZqq0bpx6hMpT+HLnPJtFGJeZQtxnMYKY/adoPLHniIMbBR9khnhNco+iBcaDGeQ6g8altZu5Gc9y1TbcQe6Qytiu/cocV4DqHyqG1lreOnum9JrwB4aGZX6UJrrtfpj3uQsgf5YojxHEbJQ1M3cS2b2WsA7vq1u4eIu3Y3KJLvkbxDcs2vZwbKHcAKLsZzCJVHbQeYKlgzGkLVa3dDiz2AVYYYzyFIHnXuBveMpJL82Zh3g2cBXDSzzdS52nYb6/6dG4g2kBgkjzp3g3cBrAN4Bvdp9TncypCx5VfjbHadS/cMFuKXanQVrC46Nnbd8wv3N1Pq16fu1ynVFR5KnSvr13CfVvfM7KKZvV7W4u6yRFhuKL2ijxtkvM+Ho7zPte0GT4JJmPqoG79y6LGZfeOPz8ONyL9eYp7TPXUzqer8HbYOxmXcYNqmbmpvEqY+6qaKcYNQ77Mqa7VqN0g2aSKNGwR5n9UNrtAkTH3UXcTlhsd+n1VZRbqM67hBbVcwiYRQwT2/I9N3Vpl2tRk3UDdYplqdxg1UWUVqQt1gkZpQZRWpif8Ds4dxQeeyYy4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 252x189 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "with plt.style.context('science'):\n",
    "    mi = pd.Series(f_scores)\n",
    "    mi = mi/mi.sum()\n",
    "    mi.index = fns\n",
    "    mi.sort_values(ascending=False, inplace=True)\n",
    "    # (?) Mutual information for each candidate\n",
    "    # plt.title('Ranking by mutual information', fontsize=12)\n",
    "    mi.plot.bar(fontsize=12)\n",
    "    plt.ylabel(\"Candidate importance\", fontsize=12)\n",
    "    plt.savefig(\"f_regression.pdf\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "c2a9842e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.79839749, 0.7687904 , 0.7226005 , 0.05708643, 0.04711015,\n",
       "        0.03377256, 0.00378066, 0.00329149]),\n",
       " array(['$uu_{x}$', '$u_{x}$', '$u^{2}u_{x}$', '$u_{xx}$', '$u^{2}$',\n",
       "        '$u$', '$u^{2}u_{xx}$', '$uu_{xx}$'], dtype='<U13'))"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r_scores = np.array([SkLinearRegression(fit_intercept=False).fit(X_pre[:, i:i+1], y_pre).score(X_pre[:, i:i+1], y_pre) for i in range(X_pre.shape[-1])])\n",
    "sorting_indices = np.argsort(r_scores)[::-1]\n",
    "r_scores[sorting_indices], np.array(fns)[sorting_indices]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "59c0b061",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPEAAAC0CAYAAABFXsrWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAReUlEQVR4nO2dPXAcR3qG39fl4KpcIpbrVAG49OXgj64csgiIFzgTCcl1wTGBKCmXDiQTgpFMlvIjCbnKdcFVkRQvc6AD+BOfSIDRXWCTRGgHhJa6cmbX56B7gMHs/PTszszO7L5PFQo7s909vT/fTn9ff/02zQxCiO7yN9PugBBiMmTEQnScv512B8qwtrZm77//Pvb29rC4uHjkuSrP1d3+rPata/1tSz/G6e/Nmzf/1czWAABm1pm/y5cv240bN+zy5cuW5MaNG5Wdq7Kteepb2rm29KNrfUs7Fx0/ffrUAPybebvo1J14cXERGxsbePbs2chz586dC2ojrVzyXEiZSdqf5Bp1t9/ENdryGtr8HqWVi479/72DJ2wKd9Rx/9J+rfJ4+vRpqfJlqbv9Jq7R9fabuEYbXwOADfN2QevQFNPGxoZtbGxMuxtCTB2SN81sA1B0WojO0ymfeG9vDxsbGzh37twRf+G9X/++VDt//d2vKu6ZEM3hY0KL0XGnjDgKbAkxzyQDWxpOC9FxZMRCdJxOGXHkE6fNEwsxL8gnHgMFzkSbkE8sxIwhIxai46QaMckFkndIfu+PT5FcarRnKcgnFiLcJ/4UwDqAAQCY2S7J80WNk7wEYAigB+C1me1klHkN4Kxv+15oXc0TCzHqE2cZ8a6ZvSMZnFhNsgfgQzP7zB9vAfgwpcw1MzsDYIfkjwDuhdQVQqSTZcQDkisA3pLswxnUKwBPctr62JeJGJI8Hb+jmtkQwBkAIDkA8Dy0rhAinVQjNrNNkqcAfAJnYA/MLM+AATcMHsaO9wH00wr6H4jTAFbL1I18YgAj+dNCzDrPnj2Lx4MWowepRuwN+IqZfREdk1w0s72S1+2lnTSzbZIAcAXA7dC6s+oTl52HBjQXPY/Eb1w3b97ci87nTTFdjR6Y2S58kCuHIY4aXh8ugJWKmW0DuEbydNm6QohDsoz4hJm9S5zrFbT1AMDJePmkT0vyCslbsVP7oXWFEOlkBbZI8rcAtvzxhwBe5DVkZkOSW97fBYADYyX5AsAynLGejfnEdyNjzaobJ2s9sRDzRNA8sZk9IvkaLrAFAPf8kDoXM/su4/yZ2OF24n9u3Tiz6hMLUYbQeeLIDz4wXJLHzOynGvsmhBiDTCNOydBaBfBFvd0RQpQlK3f6DoALcIkZ0d/JtLJNotxpIcJzpx+a2eP4CZKFPmvdyCcWInw98VuSSySPxc5drK9bQohxyboTPwHwA9xUkwEggBMAvmmqY0KIMLKMeDVlOH2qgf7konliIcLniZMGvATgeH3dCkM+sRCB88TeF74OYCE6BcCQvxRRCDEFsobT13GYcrmfUUYI0QKyotNbZvbYD6uP++ytqQ+nhRCjZC5F9AsgAOACyS9xuIB/aijZQ4hyga3H/vFVkp8ie/H+ASWE8vpwWWAP/bpikLwL4K4vtmJmI9dTYEuIEgsg4pjZZlGZQKG803DG/Z0//hGHw/QBgIcAtqM2hBDFZOVOL8UenyC5HCBZmyp2lyjTBxA30P1YmVtwd+fUtcRCiHQy1S4BvAQAM3sD4A3JNeRPMfVQIHbnh87xdcT92JB7ACfJ0yN5y8zWkxeQUJ6YZwqF8kguwN1NPwOwQPJarP4+Dv3VMvSynvAyPQfBskhE3j+3SXJgZkd0tuQTi3kmSyjvwIi9ptYmgE2Sy8msrQCGCBS788Gt+zFpnksABrFgluamhQgkazj9wxhtPcBRfzZV7M77wDtm9toLyAPO2OMG30/ehYUQ6WTNE39LcrFMQ353hy2SK17w7ohQHsmeN+DHvtwrAC/MLJqKGpC8RPI3yJiT1jyxEOGiAHeR8GdJrpnZt3mNBwjl7SAj80tCeUKEETpPfBfAK5J/D+At3AKIMwByjVgI0TxZRrxuZo/iJ0guN9AfIURJUn1irzt9jORHJNf8Pkxlo9WVI59YiECf2N91V3GYgXWb5J2AnRFrRT6xEOE+cc/MPo+fICmhPCFaSNYU04+B54QQUybLiI+T/JLkef/3FYq3Nq0d+cRClNtQ7RTcJuCAW/c7dX0t+cRClBPKuwYg2mq0T/K5NlQTon1kDaevA/jazPpm1gdwFc6ohRAtI08o72BbU7+meDujbGPIJxYiPHfa4vsR++F1pEEdlEddB/KJhQifJ/4OLnc6OiYAeKGAaF+mESOeUCivsK4QYpTgvZjipOVRTyKUF1JXCJFO0F5MAODzp/eynkeGUF7ijtqHS+eMxPIiobyzAXWlsSXmmkKNrTheEOBzuP2XADeEXgbwQc41ehhTKM+LCOTWBeQTi/mmUGMrwVU4Dei41lVvjOtm1kkK5ZWpK4Q4JFMUID7FBABeNiePIcYUyitTVwhxlMwFED53+qMofxrFou4PAJyMHRcJ5e2QHHixvKC6QohR8obTQ/ipJc/J9KIOMxuS3PL+LZAQyoPzqQdwQnn7fvqqb2bHfZnUunGiwJaCWmKeCU32eJiMQJMsFLKTUJ4Q9ZNM9sgaTr8lueQztSIkCiBEC8m6Ez+BE5AnScNhltY3TXVMCBFGcMaWX188VeQTC1Fuk/EkZwDsppxvDPnEQhQsgCD5H3DGugOXBkm4rK3MRQ9CiOlyxIjN7OcAQHI1Jdlj6sNpIcQoWeLxI8PmtHNNI1EAIcLniVuJfGIhwueJhRAdQUYsRMdJNWKSCyTvkPzeH58iudRoz1KQTyxEuE/8KYB1+F0fzGzXr2TKJVBjqwcnSj80s3ux83fh9kUGgBUzu52sK59YiHChvF0ze+dTLoMooZN1NqOJAZwQwXbUhhCimCwjHvhlgW9J9uGM8RVcTnUWIRpbMLNtv4Y4yS0Az5EiyyOEyCYr7XLTJ3d8AiehE7IXUw8BOlk5DODUPHokb5nZerKAhPLEPFNWKG/JJ3fskjwBd2c+P8amar3Qggn/eJPkwMyOSPTIJxbzTJZQXtYU08Fw18ze+AURRVubDjGmThbJSwkNr/3MwkKIIxzciUkuwPm1nwFY8Ls9ROzjMHKcxQMcldUpo5P1GkcNvp+8Cwsh0jkwYjN7B2ATwCbJ5bwdINII0djyZVbgAmU9kq/NbNuL5l3yAa8BMqRstZ5YiDHXE/tEj36RTxygsZUmIJ9bN458YiHKbTJ+HYc7IUbrissGtoQQNZM1T3wdwJZ/rCCTEC0mb5Pxx35YfdxPN6VKzTaJcqeFKLGemORvzewLABf8rg0nMeXhtHxiIQLXE/u78Bf+8VUA7wCMLEgQQkyfIGUPM9usuyNCiPGIJ3ssFZS9Zmaf1NsdIURZ4nfib+Ei0oRLuNjH4YKGHo5urjYVlOwhRH5ga9XM3gBAWsZWGyRrFdgSIiewFRmwZwGjTH2KSQgxSlZgiyS/httUDQA+QLEogBBiCmRNMT2CW5X0C//3wMymvoWLkj2EKJHsEYkCRMck14oMeUKhvMK68omFKL+h2sFTKNhQbRKhvBJ1hRAJqtxQbRKhvKC6QohRstYTj7OhWg/jC+UF1ZVQnphnCoXyasrY6pUsn1tXPrGYZ7KE8qrM2BpiTKG8CesKMddUmbE1iVDeJHWFmGsqy9gysyGALZIrXgzviFCej0AjJpS3Gonq5dWNo3liIcLnicfK2JJQnhD1EySUZ2aPSL6G28YFAO4HRKeFEFOgTMbWMTP7qZFeCSGCyRLKA8nz8T9k+KlNIp9YiECfmOQduGmft7HTJ+vrVhjyiYUI32T8YcoUU2HgSQjRPFnD6bckl/xOEBEXm+iQEKIcWXfiJ3DTSyRpOFzF9E1THUtDGltChM8Tr0pjS4h2EiweHz/2iyOksSVEC9GuiEJ0nE7tiiifWIhwn3grGlKTPG9mT3zCx1SRTyxE+DzxWLsiBgrlpZYheRfAXV9sxcy0gZsQAWQtgHgM4LF/fJXkpyjYFTFE7K6gzADAQwDb0fNCiGIOotMkvyL5A8mvSS76cwt+GP1DYr1xGqlidyXK3IJT2px6jrYQXSJ+J96GuwserFwys3cAnnhjXjKzlzlt9VAsdpdXZgAnydMjecvM1pMXkFCemGcKhfIAnDCzP6RVNrN3JEOVK+P0QsskhOQ3SQ7M7IjOlgJbYp7JEsrLXIqYQq/g+SGKxe5Sy5C8RPI3sfOtmdYSou3EjfgXBWWLliI+SJRJE7vLKvMaRyV7+sm7sBAinbgR3yd5n+R78QIkj5G8j8Pkj1RChPKyynhDHsTuyKtp15AogBA5yR5mtktyE8CuX7k0hMuXNgCfFwS1ojZChPKyykgoT4gAcpM9vBLlP5A8AeA0gJ2AqSUhxBTJSvZ4A0DGK0QHKBOdnjryiYUoscl4G5FPLESgKIAQojt06k4sxue9X/++VPm//u5XNfVEVE2njFiiAO1FPxLNIZ9YiI4jn1iIGaNTd2Ix32jInk6n7sSaJxZi1Ceu1Ij9AoYV/z+p6pFbJqRu5BOHBrX+97/+PM7LCKbu9pu4Rtfbb+IaTdw0ylyjNp84pp+17RczjMjsZJUJqTsO//fff6mimam138Q1ut5+E9domxEnqfJOPInGVkjdA9JecOivcVq5ZHt1tz/JNepuv4lrtOU1hF4zjdC6VV4jqx7NLKjBIvw64GEks+MlaB/6lVG5ZeBWTOXW9ef/HcDfwamDvERsSAHnI8SPJzlXZVvz1Le0c23pR9f6lnbuHwH8zD/+HzP7J6D+6HRvgjIj56NOCyEOqXI4PcSYGluBdYUQKVRpxJNobIXUFUKkUJlPDBzZogXAgVIISL4AsGxmw5wyqecr6NOSb+9lFe013X4TzMhrOFGnCk2bv0eVGnEbIXkRLvI9BDAws0q3Z627fX+NJaDWL9BEr4HkMoB1AC/M7JqXd7piZtcq7+zhNe8A+E+4DQ9e+u14V7K00yu4Xmu/R/OQdtkH8AmAP/njqo2s7vYB52q88tvr1PFDMelrOG1mF0ieIvk1gH+B29GjTtYBnAXwzyRvw4k6bgOoxYjR4u/RPBjxvpl9TPIUnHJn19oH6v8CTfoatkke81sA7fq7Sq/SHibwWwwdbPwHHIwI6qK136N5MOJt/8YM4H6tX3asfaD+L9C2H7KP+xpew/3IbAKAmT0iOaywfyOQ/AruTrwFN6TeA7BQ4yUnfY9qa38efOK6Ax61tu+vsQD34Q4AHDezbytufwmo1eeu/D3yd93nAFbgfkAW4BKEKn1vYtdbAhTYaoS6Ax5NB1T8Nev+Iao0aNPEe+R/2JbrfN8T12ttYKtTSxEDWQewCxfw+COcz/RBh9oHyTskv4x+nQG8JflRlddI0AdwHcAluBTYSan9PTKzd00ZsKfq96iy9mfuTpwGyWUze1xcsh3t+7vMWQAfwn2gx+HuarVM2ZC86P3YUwCsjiFj3Z9B3dT9Hk3S/swFtuoOeDQRUJlC5LXS4NwUgk5N0NoA6cwZMYAdAPfgAh63/V3tYYfan4YR9KPpoYraq/09mgJVv0eVtT9zw+m6Ax5NBFTqjrw2EPxrNOhUB10KkM6cESchuejvZJ2hoR+ixnzuLlL3e1Rl+zNnxNOYApqUaeQep/Why4GnJmhrgFRTTO3gtJldAPDA5x7vo8bcY5JfkbxPcs3nYwPdDzxVSt3vUZXtz1xgawqR3SpoOvd4FgNPVdOZAOksDqdHIrskP2r5cHoBwMdmthk7V9vQbRYCT3XTpQDpLA6ndwBcAfAj3C/c93DZMK3FZx9tJs7FRxInarieDDgGyWWSf/TuDOC+M7W5YcnPIDakLs0sGvFzuF+4R2b2sZn9sq6k+LqYQtqlaD4ukfyM98f9jGduOD0LaAqoeXy21Csz+8kfX4SbIfhlTdfTFNO8oSmgemk6LpHRB00xzQqaAmqepuMSVX7GMuJ20rng3KzRQFyiss9Yw+kWoimg6dNQ2mUln7GMWIhA2hqXmLmMLSGqoEtrouUTC5FOZ+ISGk4LkUKX4hIyYiE6jobTQnQcGbEQHef/AbEr8061nLxuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 252x189 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "with plt.style.context('science'):\n",
    "    mi = pd.Series(r_scores)\n",
    "    mi = mi/mi.sum()\n",
    "    mi.index = fns\n",
    "    mi.sort_values(ascending=False, inplace=True)\n",
    "    # (?) Mutual information for each candidate\n",
    "    # plt.title('Ranking by mutual information', fontsize=12)\n",
    "    mi.plot.bar(fontsize=12)\n",
    "    plt.ylabel(\"Candidate importance\", fontsize=12)\n",
    "    plt.savefig(\"r_scores.pdf\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b729a5b1",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "956d98a8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "0ed1dd87",
   "metadata": {},
   "outputs": [],
   "source": [
    "# y_pre = y_pre-y_pre.mean(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "117c7b99",
   "metadata": {},
   "outputs": [],
   "source": [
    "# tau = 3\n",
    "# verbose = True\n",
    "# scale = np.log(len(y_pre))\n",
    "# # scale = 1 <- generalized UBIC\n",
    "# per = 75\n",
    "\n",
    "# post_means, b_bics, b_uns = baye_uncertainties(best_subsets, (X_pre, y_pre), u_type='cv1', take_sqrt=True)\n",
    "# predictions = X_pre@post_means\n",
    "\n",
    "# b_bics = np.array(b_bics)\n",
    "# max_complexity = len(b_bics)\n",
    "# complexities = np.arange(max_complexity)+1\n",
    "# d_complexities = complexities[decreasing_values_indices(b_bics)]\n",
    "# d_bics = b_bics[decreasing_values_indices(b_bics)]\n",
    "# slopes = np.diff(b_bics)/(np.diff(complexities)*b_bics[:-1])\n",
    "# thres = np.percentile(np.abs(np.diff(d_bics)/(np.diff(d_complexities)*d_bics[:-1])), per)\n",
    "# # thres = 0.02\n",
    "# print(\"threshold:\", thres)\n",
    "\n",
    "# lower_bounds = []\n",
    "# for k, efi in enumerate(best_subsets):\n",
    "#     assert len(efi) == np.count_nonzero(post_means[:, k:k+1])\n",
    "#     com = len(efi)\n",
    "#     lower_bound = 2*log_like_value(predictions[:, k:k+1], y_pre)-np.log(len(y_pre))*com\n",
    "#     lower_bounds.append(lower_bound)\n",
    "\n",
    "# last_lam = np.log10(max(lower_bounds/(b_uns*scale)))\n",
    "# print(\"max_lam:\", last_lam)\n",
    "# delta = last_lam/tau\n",
    "# now_lam = last_lam-delta\n",
    "# last_ubic = UBIC(b_bics, b_uns, len(y_pre), hyp=10**last_lam, scale=scale)\n",
    "# last_bc = np.argmin(last_ubic)\n",
    "# while now_lam >= 0:\n",
    "#     now_ubic = UBIC(b_bics, b_uns, len(y_pre), hyp=10**now_lam, scale=scale)\n",
    "#     now_bc = np.argmin(now_ubic)\n",
    "    \n",
    "#     diff_com = now_bc-last_bc\n",
    "#     diff_bic = b_bics[now_bc]-b_bics[last_bc]\n",
    "#     imp = abs(diff_bic/(b_bics[last_bc]*diff_com))\n",
    "    \n",
    "#     if verbose:\n",
    "#         print(min(last_bc, now_bc), '<--->', max(last_bc, now_bc), \n",
    "#               np.nan_to_num(imp, nan=np.inf))\n",
    "    \n",
    "#     if (diff_com > 0 and (diff_bic > 0 or imp < thres)) or \\\n",
    "#         (diff_com < 0 and diff_bic > 0 and imp > thres):\n",
    "#         break\n",
    "    \n",
    "#     last_lam = now_lam\n",
    "#     now_lam = round(last_lam-delta, 8)\n",
    "#     last_ubic = now_ubic\n",
    "#     last_bc = now_bc\n",
    "\n",
    "# best_bc = last_bc\n",
    "# if abs((b_bics[last_bc]-b_bics[last_bc-1])/b_bics[last_bc-1]) < thres:\n",
    "#     best_bc = best_bc - 1\n",
    "    \n",
    "# last_lam = round(last_lam, 10)\n",
    "# last_lam, last_ubic, last_bc, best_bc"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:pysr]",
   "language": "python",
   "name": "conda-env-pysr-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}