Examples/Burgers-SampleSize-Copy1.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",
    "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": [],
   "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",
    "K = 4000\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=K, # 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",
    "# center outputs (optional)\n",
    "# y_pre = y_pre - y_pre.mean(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "65777e7c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Set parameter Username\n",
      "Academic license - for non-commercial use only - expires 2024-06-04\n",
      "u_t = (0.101490 +0.000000i)u_11\n",
      "    + (-0.905064 +0.000000i)uu_1\n",
      "   \n"
     ]
    }
   ],
   "source": [
    "optimizer_name = \"miosr\" # \"l0bnb\", \"miosr\"\n",
    "assert (optimizer_name == \"l0bnb\") or (optimizer_name == \"miosr\")\n",
    "if optimizer_name == \"l0bnb\":\n",
    "    optimizer = L0BNB(max_nonzeros=2, lam=1e-3, is_normal=True, normalize_columns=False)\n",
    "elif optimizer_name == \"miosr\":\n",
    "    optimizer = ps.MIOSR(target_sparsity=2, normalize_columns=False, regression_timeout=100)\n",
    "    \n",
    "n_ensemble_models = 50\n",
    "model = ps.SINDy(feature_library=weak_pde_lib, optimizer=optimizer, \n",
    "                 differentiation_method=differentiation_method,\n",
    "                 cache=True,\n",
    "                 feature_names=['u'])\n",
    "model.fit(np.expand_dims(un, -1), t=dt, ensemble=True, \n",
    "          library_ensemble=True, n_candidates_to_drop=1, n_models=n_ensemble_models)\n",
    "\n",
    "# print found pde\n",
    "print_pde(model.get_coef_list()[np.argmin(np.sum((np.squeeze(np.tensordot(X_pre, np.array(model.get_coef_list()).T, axes=([-1], [0])), axis=1)-y_pre)**2, axis=0))].reshape(-1,1), model.get_feature_names())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ae80736a",
   "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": "3b30bdf4",
   "metadata": {},
   "source": [
    "#### Model selection by the UBIC (Algorithm 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8baf3445",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "threshold: 0.1840442876442533\n",
      "max_lam: 3.5738341379586243\n",
      "1 <---> 1 inf\n",
      "1 <---> 1 inf\n",
      "1 <---> 1 inf\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.0,\n",
       " array([-26242.31758607, -31080.27754417, -31057.76671574, -31047.68458305,\n",
       "        -31030.75643991, -30998.58497224, -30843.66312884, -30814.3812685 ]),\n",
       " 1,\n",
       " 1)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8bf2ca4",
   "metadata": {},
   "source": [
    "#### Percent coefficient error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "179c5e17",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.036833226680753, 0.5846604704856845)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "perf = evaluate_coefficients(np.linalg.lstsq(X_pre[:, best_subsets[best_bc]], \n",
    "                                             y_pre, rcond=None)[0].flatten())\n",
    "perf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d7b76b36",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQwAAADBCAYAAAA6jY5NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAu4klEQVR4nO2de5QcV3Wvvy3b2DgEtyWDeQSMWxDwI0EeSwkPJ1eSRzxubgKYkRwluZDEeOQLgTwwko1gRSwwZmRYQAgJIwVwwg3ElgjJhQskkm1BCCGxNDZYxiRBE2O4Jjw8bocEjI39u3+cXTM11dXd1Y/q6kd9a9Xq7tNV5+yqrtp9zj777G2SKCkpKcnCiqIFKCkpGR5KhVFSUpKZUmEUgJlVzWyfmR2IlVXM7ICZzRQpW16Y2XYzmzSzaTOretmEme1L2XfCzGbMbMq3CT8+um5HYt9tN7Nj/T+j8eT4ogUYRyTNm9ksMGtmVUnzkmpmNiPpYD9kMLMpoBa1F30GKsC8pLlOyhu0NRG1ZWaTwIKXrQWuS+w7CWyTtDlx/MWSdruSPV/S/tj3c9F17PiCLJehAky7zHsa7NP0/JPXd1QoexjFMgtM9btRfyC2JT5vknTQH8SZTsqbsBJYiD5IqvkDtuzBd2aBS+MFvm/qg2dmFX8oKy1kaIe1zb5sdf7J6ztKlAqjAPyGmgf2UMyNtRY4EPu8BYh362v+r95ueSMOA9V4gQ+9diTKqhAUSkods8kC/xePjmnYw2kXV0BpMkS0Ov/k9R0ZyiFJMayNDQWim61GeLAy4UpnS5Ndrk978Lytw0D8Bq+w/AFZIPQK2i1PxYdbc1Hb/nrMy+NDieg6pNURH25UXVlsIkWRRHR6jTJQocH5N7i+I0OpMIpnltDL2Bc9FGa2ndADAahK2p08yG/01PF1K/xBbbVbpUflUZuLQwoz2yZpm5/nHjOb9O/nGtXjQ4+af5yXtN/MDuI9jDQbRjfXqAMq8XYzXN+hpBySFM/1xP4FY8ay/T4+nvcHq2vMbJqlf+d1wCYfBtRY/qCuJCisdsuzyhD1Clb7Q70AS70I7xkkqbMrxGwhkBjy5EyNlPNvcn1HhrKHUTD+b7T4T+mvtdguNWBV8riYJb8Re5Ld7bjF38zWAQd8xuZ6lhvuKpLmzGy+nfImssRljupInmfENmAvsDlx3ELKvtH3VVJ6Jp1co4ykXi9CDylqe/H6dlD/wFIqjD7iN/YM4d84bvWfZemB2MPymzF1nO43et1QJaMcE4QxdsXM5l1pHPApTaL2XZllLve6jwCbGzwoV0ra4XXMp/UkfOp1ITYsqwELMSWziaV/cQgKdhsxBROrq0YH18jPaxNL1yeyNx0BLmx2/r5f3fVtV4ZBxcq1JINFZKiLegPROL/Df8JCMLOJFL+ECWBlwpYR+TIs9HKWoyQ/yh7G4DEJnO//pguEsXCN/hnvusLlrvtHTVMIKT4YJQNO2cMo6SmxGY+SEaRUGCUlJZkpp1VLSkoyUyqMkpKSzIyF0dPMvgL8u3+sUD//n1bWzr6jXDZo8gxS2aDJk6XscZKeQadIGvkNOBR7vyvl+7qydvYd5bJBk2eQygZNnixl8Wehk20chySHuty3m7JhbKNf7QxjG/1qp9dtdMxYzJKY2SFJ65t8v17SoT7IkXs75bkMXhv9aidLG62ehVaMYw9jETPbamZHgRvM7KiZbc2zvX7cmP1oo1/tjEobebdjxlYzjoJuMOOoGan3sZmtp8tAQ2Nh9EzDlcNVwCXA54ALgPebGZI+UqhwJSUZceWQch+DxLL7WNIhM6sljo/Ww2ySQkAjMyKX/aq03MN4XHoYFTPb5Ro2YifhIn8eeC9hrHeJl5eUDAs7gd8CngY8LHETGe9jMyaACYmDwIQZVVcWeFlcoQDjozBqknYluoVnAZ+T9EPgfwBnEDT0WQXIV1LSKWcDHyKsQTrFyzLdxxJzErvNokDGzBPWLkVrgeZJRA4b2yEJcAeh+3YTIY7BBHCml5eUDCRmGPBs4PkSvw/cBVwh8Rex3S4gdh97z3q9f6ykVLuWpRilye+XxWIZZ4VxFcFmcQlwK/AS4LmUQ5KSAcWMFwBvAk4D3m3GCuBK4K1mfIuYDYPYfew960OhDlufrFfioBmbY7aLhvFZx1ZhSPqIx118D6H79n1gujR4lgwSZpwK/DIhiNIJwNXAxyUe8l0+4uFDo/v4DmBn0uDZoO4Z4JgbNmsERXEzS72MKono5+Pih3Er8FcEL7dDKd8/iRDp+XEa8AsSi9pVBXYoRKiqEMLaQcjpsTK2z9Vevg64WSF47gThn2mxDq97AriYcNNAGMNOKiUIcUKmKS1PLDQNy0MCFoVfm8Uo7cOCGWcCrwV+FfgE8CqJ/+i+3iU/DDOqLIWG3CyFlBdmbMeH6VIiYlk3bqLDstHCHRYw4DvAE4qWNeP5TAOzibIJQoTx6PMUMJPY515C/Mm6OghGs30pdR5pJUuD8mM9Pucpl3Eq3MjZ9/Gyai/lyed3lYHWg04CvRB0Fain92SrZ6HVNi6zJE1RuJKR4bNrIocwM3uoHw5hMRoGyo3RaHzaVsYxaBjdOyLKQ9I1WTKtNdtHoXcxsJnIzDjRjF8HbgH+CHiyxKckdkrc3bt2bD1dOm6VCmOJniiMmEPYq4GT/PWqPiqNRnJNA1crJSBtuxnHYkzTWKFcRxje9IIsmdZa7VPpkSw9w4zHmPFown23FbgCOFfiX/JoT2E4XuumjrE1eqYwRxgvdstO4BJJN/nnm3wm5j3Q2hDVY6pR8F3/XGuwX9aMY0nWqYF9Q8FWkpYOsUL72cgqCfnSMq212qeWFpw4b9wTcydLBsmrgNuA3wFeCvyKxKeA5/dTrk4ZF4VRMbNdNDB6OnPAO5pVYmZZDaI3pmW+ynK8pF6mzFqWVdzMjpjZguqD72bNOJYZ79HMJzOSqXfZyCpt7nMzwcDXN4WRcNv+e2Ajocf248A7gZ+U+E6/5OkF4zIkqane0zPJPPBoM3tMox0kWasNuB3YmCjbCNye8fgszFP/D7syw4N9kDBbkjyvtjKOtcJTIxykd9npa7TOtJZ1n34SLT84Dfgi8Ezg14G7Jd7cb2VRLj7rIZJkZrcQuud/00VVcYewVEeablGYSk0a/ioZDp2gsU2i7YxjaW1GPQtXQvOe/Gd37PsK7Wcja5RprJ19qjQx4ObEWcAR4OvARcCNhGeukOUHSll81i6lwlhOZPjsWGGo3iHMHWl67hC22ZYnbY4nCKoSDI7xDGHrgDkt+WFEGcQmfWahYcaxJjIciA85vN3VWu5/sWxYog6ykal1prWW2chcro4yxXXBHQQP4v0SNwCYLXfbHjrynlsehI2Mc88Eo+e+PGUZpY3wr53qhzFoGwmflP60qa2gedAG0An+Og/aWuB1ONTN8eNiw8hKz3wxxgGFXkOlaDla4b2OZtPDefFR4D+BPwTuJ/Q6M7ltDyrjojDS4mGk8S/A6WZ2av4ijQaSdseGPQNHZMhVMQmR1wP/KXGOxHES5xapLErHrezU1HqWBEkPEazZa/oh1KigAc6RKqmm4taRvBT4y4LarkM9cNzqu8Iws4qZTZjZVGTp98/H3E/gSKx8yswmo8VM7ZR1QTksKekVpzNACqMXFNHD2EJYPbgfFqfhVkpaLel8wnqG2aibq6WVlJNZy7qUr1QYJT1B4sVSfSb7YabvCkPSHi1Nu1UJc/bxLmM0BZcWKixrWTeUCqOka8y4woxNRcsRZ6htGD5nvxBXFmY2HRsPVxKHrGqjLElk9Mxi+PwycIaZParFfiUlqXgYvVdC71aa9oKkDcOMihkTZkx5MJ2o/F4zDnhcjGXk4rjlw4Sk63KyJzElKbnkeBNL6wxqKXVkLUtSk7SrxT4ASHrQzG4nuPH+fZZjSgYPv9kPSv1bOxJjLfADwp/PILMFQGKPGevMmFaIvrVZSveKzUVhtLKae4Sm3f5+QtJcyjqGtFBhlYxl3XKEMCwpFcbwslrJaFH94xzgwxIDHb1Ny3OOxJ+dihnVNPtLEbMkk8BMNCPCUu9gJbF1C650qr5/RUuBUVqW9UDMgbVjeO6IfWYciXJGeNdyn2+VnNvvWf2d1NXG+a/08qm8r0kSiWsl3tTPNrvBQ/UtxHoVK4EFsxRnt6JcVAfZHZbQpfxi0XI3lk/ToESIPk2Acg9DB8rsCg5qGqavnbraPX/QPi+bBFX6+NucDXpr0ffIcplYD+zy7dYUmbc3OJcZ0FS8bFwct9rlKPA0Mzupk4OjXJdmPNQs12UOZAnR1y3thLo7v4d1ZWEBFrN1LQCTEgelvi5rnwIe2cf2WiLpkILj4i4SjltmTMmHbm4AnfaMaKmMi8LI6hoOgKT7CW7iP9VuQ7GgKYkQfX1TGpgx6ZbuCX8/m/h+yreJ6ObwG2XSjO2xrvwxL9vnZZOw+N1Eor1JM2ai7r9/f6SRPMm6/PXemDz70qz0GanCYp7Qbv1y2uUiwhqSgcevzYwP7yLzwPX+XZQycbk9sujuUp+6ZIc6OOYDwLb2j9NR0IZE2QbQ0d6dT8MueSX2+QBoIta1nIztN+vvK6BZL5uJl8XrSHT166KIx/afjOqJjm8mT7KuSJaorm7Ov5j7TKtAt4OOK1KO5jKWq1Xzos7wacYuMxTbzvdtsYyQ6/JzZtwdK3s7cJYZexLHP8GMX4x93tXjc4is3PfEyi7GreESNYVcFBcD9/i/ezxXxUqF/JtNvRW9johKm/LE2UfvAgf3HYl7CEF8H2q5cwEMtePWEFCnMCR2SVhsO+LbYhlh7v0CiSfEyi4H7pCYThx/t8THY593ZZStQYi+5ePT5OcmVIA5Vw5z0qKHYkObiFvWo/czMbtBQxrJE9WlYKWfdMV1uElVmc6/35jxduBxRcrQDPVj8ZmZrTGzi8xso5mt6aaxIeOLwDlm9og2j/MQfWww4wQzNhBC9F3VK8H8waomiisZD78OllyW/eHcl1LWiEgpVH3faeAeKR7xK/OU9LK6nIMEY2Wt0UFdnn8umPETwG8A3y1SjrxJddwyszMJFuxTCBrpMHAqcL6ZXUbIoDUr6c7+iNk1WaKGL0PSf5nZnYQhxq1ZG5I6z3XZJpvdKJgSoi8MLfz7/QRlsNqMwxJzbuSKEu/Oe1k1ZiBc8PdrY95/EbOuJOLt7kj0MFbGZJgm3D918gB/nqgLgkJLKoO2zr8gXgx8QuLBguXIlbrcqmZ2ISEm7o1ND8y43yAQzyfZ5nH/G7hR0gd6L1VJGmZMqBh37q4w48PAdRJ/XbQsjXAbxrskrem0jrQhyeEsSkDSDfi02QgzsB6fo4YtLX5qtS5oUPlV4ONFC9GMXGwYku6L3pvZo83sFU0EuK/RdyNCqTD6xwEf1jQzdg4kZvwC8AsSDxctS940smG8lLCSbYFigqcOCrcCP21mxymE7yvJCTVYHTkkvAr406KF6AfLehhm9jozu84/Xk8IzX5r36XqPW15ekYo5ND4d+DpOchUMgK4Z+sFwCcLFqUvJHsYByVdE30ws0vN7NgwGDZbUFPGeBgpRMOSQY9tUFIMFwCHJL5XtCCt6LnjlqRbEp/3EqZTx5nSjlHSEIlP0Jv8sbmTu+OWO2rd645ba7ppaIgpFUYKvkhtrK+LGSeb8SYYbd+LOI2MnucRclPWWOaMY6cQ7BrDPkRphzngPDNbIWnkreBtUGREq0HhBcBzpcGOrNVLGvUwzpT0PElbJF3m2xZJzwesnwIWjaTvEhRnFu/D3HGPzFlfjpz8ruJLxDtdFj4UEa0GiKFZyg75Lj6rNTlmGLVpR7MkMQZmWOIrR/cB8/EFYE4Uam3xn9+MYx3UfwA4HE11+rqOq4EdiTUeVbocEw8rZhxHcHH/q4JFyUzShmGNo4ZPedySusRgjRTGqWb2PjO73Mxe4dvlZvbHDEHy3RRqypAqsQltKQwzfsWMO8142F9/pcN2m3Ed2SJWtYp61Q6DENFqIPAl7Kslvlm0LF3gScXwpGJMxwLneGKw5QGIUhWGpI8SIhbdR5glOdXfXyGpq9RvaakSvbzjtIg9TpWYRmaF4cphL3AGYfh2BrC310rDf+TFH9N7G8viVrSKetVF80VGtBoIzHgZ8Nii5egGiT2xhYXR/dM0MVjDWRJJ90naK+ka3/b2yBW8LlWiR/yO8pbMRwrF5eh3qsQ05oAJM8tiv3krcHKi7GQv7zXzsZmKSko8jDn8x/d/jMOx97VOH3a/0WoSu4fcQ7MjzDgeeAeMhiu4LY8aXkl8vSwxWNsBdMxsY+eipadKJNzI+8xsgpAqcY7iUiWmyfxNwtTZkzLs/uQ2y7thlvYjVLWKelXSmguAuyTuLFqQLJjZ+ijzH+kmhalY1LQaTRYANppWbbTgzAhOKl1PqyZTJZrZLMGYF/1jVRKHdJ0q0d9njomRIBqW3NViv7sIw5C08p4icdBnTGazLgnPaG+Yh7q8oCuHcdl5TvwSQ5SV3e/3Q7A4U7KIJaKGk55AbJGGRk/fLLFBsGU0xW0K04kt2f1dTJXo3x2UtBqo+TCjRo9TJfbA8JnFgPh64PuJsu97eR4cpMeehoMY0WrAuAL4g6KF6Ja0qOFuG4uCKVWSQ85GqRL3EOwMN9Q3Yi3T16vNVInARPSZMH23hWJTJaYxBzRc6h8h8WG3dLyVMAy5C3i9xId7IYSPN3d4oJndhF5Z9N0UwT9iWiFfZsuoVx6Fq5bS1KBFtBoIzDgLqEr836Jl6RZXBqtTyqNnse43r4u4lTfem5hlaT54B+Fm3kK4OauRjcPMtuNDgZiCyVSWaLOjiFuJOs4AviDp8d3UUzLcuL/CgxJvKFqWdulFxK2+K4wi6JHCMEKA13PdCFoyZphhwL8CW4bVntPtszAuaQa69fREQbMOjMdnSSGcCxwH3NJqx1Glk2nVp3jqgUfnIVBOdOvpGVEqjPHmKCHnzOh3yxvQSQ9jv0fhWtWtT8YQUiqM8eblDHHekUIyn0la66//NmbL3KFUGGOLGT9JmMEb2tgXuQXQ8WjhwzTk6BfHCAvzTitakJK+cxHwsXGIDN6MRj2MTcBa6N4VfJTwADq3AOcVLUtJ3xmq2Bd5kYwavsbfHgAWzOxyRqML3vUsSYwjjMY1KWmPrcBnixaiaJKenps8d2rkan0Pvox5yOkmaniSOeAXe1RXyRBgxrOAbw973tQ8ooZfI+ky4FJJWwguxFu6aWAEKQ2f48cMcE7RQnRLnlHD3+YN3CDpT9IihpvZxjE1jP4z8HgPiFwy4phxOvBM8lujNFTUKQyPGL4u9nkNiZWLvs8VMH7Rljxl4peANQWLUtIfXgR8SuL+ogUZBNKSMd8CvM3MrnNlsY4Qci7OFQTbxr/lLmFv6KXRE8phyTjxUfILTdBXcnPc8uXpe4HdhAdj0mdMMLMzgdWSbklmShtgeuUaHtG1wvAVtiPNsJ+jGacAPyUNzR9jU9JsGB4h/ECi7F4zDqSlq2gUDyOKkbm4Ht7M5s3sdQT/jM3diT70zAGv7fRgD1bcNGbIiLDfzLanhRxoRSxWa5HX6ZcIwYkOFShDrkjsN6uLPr+5UazWzK7hHhT4GkkXSxoJjdsFXwbONLMfa/dAD024WlLLQETDjp/jKj/nzHg0+ZUFKwsYX2etSkrOGyDd6HlR1lrb2XeUkPQAQWk8s4PDt0FXIf6Hjatpw5fHY7veHAsUXQhmPArYCHy8SDkKYiWwkJaKom5IIukvzexSQoazg5LujH9vZk8huI6fIuntuYjbe6IgwJ0GAE4jsmN8vs3jJsahdxEhqZa1h+E2j2MD0LMAeAB4ocS9RQvSLW7sXO8fK632j3KVmFHzIMGLv0eqDUPSXvczmDazdQTlYcAphPnojw7ZsKSXnp4Rc8DPtHOAPzh1ysLH6ysJUdT3x8qqnYz/M8rSss0eyjVvZhOePqKRPBPAJknJaOVFcT4jEiinWdTwJFH810YRxVolMrrGkzBfHCVjlvT2bpRFk8xn2z0x0fZY2aBkPkujk5mSulyk/qDMEVIPxnOMXEyKcukFWdrssVzHaJ3Meh/ZUj/mjhknAp8ERt4x0aODr41SJALXe3mUMnFZb6+IEH1pmc+2E3oBUY6SQct8lsZtwNPN7KQ2jqlCXXLkmg9RLma5N+Ek+UXqztJmL+WqUe/8NxF7P0MIzDQoQ7VJ4DaJbxUtSN54btxTI8XgGe3mJPZL9banNKPn68zsZjO72u0VmNkpZnZhmot4+wKmZj5bhSf6dSYZoMxnaUj6ASEg7LltHlpJ1BPJPhmLlh4leap1KWYqWdrMQa7FY1xBHIkpjWmCcXRQGMnZkV44bqXZMA4SjJ2L4zfPqXqDK441HqKvK2I330HPdbLNzKI8CcfobeazvIiGJYcz7r9ASh4I7xXF/10X/8XNrCpp3veJ20AmJO1uVN5KkDba7EiuRHOVRD1XE5TEXp8VOZyXcuyQ9wFfK1qIXiPpkJnVuqkjTWGcqQYZ2iXdZ2atsozFjWVxomTLEYuZz/zfbIcrkWMs9Tp6lfmsF6kS02jXjpGWgjD+XcRmlnLN1sysAqF3ZmZHgAsJCja1PKrEr2c1cd0ztdmNXCntrCY2jPGZkx2E6eXZBscUghlPAr4h8e2iZRlEGnp6NqHSaocOMp9BuLH3m9k6/+eMt9Vt5rM8ZkkgKIxfy7qzpLm0KUbvZW2O7DGE5dSbiM1OAAf93/w6/zfe3KIcgpfilYS0l2232aVccaIE2/H295jZNoLCvT7lmKK4AvgGgzVEGhwkLduAtyXLEt+/rtn3rTZCt/YYIXLVEcIYGcLNPUXo0kb7bvf9t7dblmjzUDcyNzmXRxHypp7QxjGzQKXNdmYIynCfv1YJD1pqefJ653HuzeRK7FMBZhscXyEok1zka/98tAJ0N+jpRcuS4292qJvj6zKf2dLS9VdI+l6s/NGEBWlXqwc2jH5iPch81qTuO4CLJX0p4/5VYJukdrwfpwlDgwXCWp4Fhd5Yanni2KlkWa/I2P52BmsGpCFmPBvYK7VtyB4KrAepEtM8PW8xs73ALWYmwnj2VILz1mXDpiz6QGTHyKQwFAyF90RGw4zHxN2k51qVRzRyFOsVGduvDYOycL4GvLJoIfJCPTB6NlreflDSU4HnE6JvbZL0NKVkcx8Seh0PI07bDlwK9pupljt2iaR5NfGu7ANTKnhNSFY8b+opwN8VLcsgkzYk2agmCYpafT+I5Dwk2QC8WdIFedRf0h/MWEPwvXiqNLqpELt9FtJmSXab2c1NjllLLIRfCbcAzzSz4xTC95UMJxcBHx1xZbGeHBy35kkZj8bop2PUwKPgU/Bt4GnAV4qWp6RjLgJeUbQQedILG0aawrhaTULvmVlWr8ZxIrJjDKTCMLMzCNbxlxQtyyBixgrgHcA/FS3LoJMaBNhCbtWnRGVmdpGvLbmomTIZYwY2KLCZrQCuBf6hYFEGmccDf6Yxz5uahbTFZy8F7iQsDrrcQjCd1xPiYVxmZsPoAZfnLAkEhXF+TnV3y2uAEwj/oCXp/DXw80ULkTd52TDWSVrpDVwKSNLaWKNv66bBgqgpH9fwiFuA88xshULC5oHAzM4GdgLPKg2y6ZjxFOAMxmA6Nc2G4XEvtklLa5y8rAZU5dG3ItL8MBZnSCTtpT73SLMZlLFE0reB7wFnFi1LhJmdAPwZ8AZJyRgcY48ZW804SjDyH8eYRsJXIkBOLHCOr0penqwsrYeRnFa6p8X3JYHIjjEoD+dO4DvAUDhO9RMztgJXAZcQhtpPBa4yA4mPFCpc8awDrvP3UXyZxZXGaQrj9RbieEZUPV4FhIt7IZC6/H3MiRTGvqIF8d/vfwHnKemZVwJBmV5CUKjfkrjRjH8F3gNjrzAqic/L3Cga+WHEp06T06ht5ZgYEPKIGp5kDvitnOrOjJk9EvgQ8BpJdxctz4ByFiE37k5gKyEkwue8fOSw9qKG12gSXyZNYexQkyC/Zlbk2oROydvoCd7DMDMr+F/9auAWSde13HMM8TUj/wVcCjxL4qv+1QXAHYUJliNqI2o4wUZZ8fd18WXSjJ5Jm0Wy8WFKL9BP7gYeBn6iKAHMbCNhUduripJhkDHjfHf93g2cBDzJjBPM2AC8n2DXGCuSUcPdCFr18ooSKRPTehh7PXxaPCgvkv4jJ5lHAkny3tcE8PV+t28hj8wHCXFMFlrtP054FrN3A//NjLUSbzHjGMFmcRahZ7FzHA2erhBOTZRFMVnrQjumKQwjdEmWVeKxPE8hBAgulUc6RwgK468LaPvdwCclfbqAtgcW97M4CHwGOE/ie7A4GzJWCiIvx61LFaKEN2r0IoZvlqQfRk8IdozfyLH+VMzsJYQx+Jp+tz2omHE8wSHrLuCVEn9bsEiFk0sAnWbKYoipSdqVs7KAAtaUmNljgT8CXi7pP/vZ9qBixmqC5+YbJB4slUXvSFtLsrHFMXV5NdrFQsaySVueKrHjtIhpZQXxNeCRZva4fjRmIbT6HuBaSX/fjzYHHTNeAnwB+AuCr0VJD2kngM4qgutzVy60nlYgSsIT5SKZgMWw9tMetr7SaZka5+HIlZjh8zzgU31o8uWE3+TiVjuOOmasAn5A8LTdKHFbwSKNJO0E0JnvRUxPjzE550lw5j0o7jbq3VFXdVFWiMJwomFJrgrDY1zsJqQR+GGebQ06ZmwizBC9Rho6+1rfyMvo2dRxq4esZWndRSXxXTtpEdPKimSOnBcyeYyLDwLvyJreYBRxJ6y3ExJ8v1xiWINU94VcIm71QllYhlSJWp5Vq5ayfzdlSfJKlZjGHPlnzXo1cCLhYRlLzDhN4rtmHAGukih9T/pAJ6kSW9IscY4bOo95+Pka4WFPc0etdFGWpB+u4RFfBVaZ2co8HKjM7CzgDYxpjAsPp/c7wOvMOFviwwWLNFak5iXJmVkgyvpdkbTHFUw1Vnawm7ICzmkRD6BzK8Hw2VNiMS7eOI4xLsx4LPA3wEuB50jcW7BIQ0UvbBh1eUlGkTzzkjRo713A3Z6wqJf1/j7wbOCF47Zs3Ywf97eXAH8o8aMi5RlWun0WiuhhjAM9d+Ays7WENH6XjLqyiKJhmfGQGV824xBwrcT3JN5VKoviGBeFkXcQ4CQ9VRixGBe/Len/9areQSQWDevVhMC8pxCu5ceLlKskUA5J8mnveOA+4PG9WKhnZu/0un65a+EGHI+z+QGCB+uzgZOB/wDeM6pZ1ftJHqkSS7pE0o/M7DbCYrDPdlOX527dDDyzB6INLGY8A/hd4BzgV4FPSmHGy4wTGNFoWP2kF0bPcRmSFEHXwxKPcXEtYQVx08BGw4YHrtlkxh4znkAIaHMM+BfgcmlZFrmRjYbVT9z/qBYvM+NeMw6YsT1LHWUPIz/m6D45zruAT0nqx7qU3PGewo+A/05QhPOEoMn3S9wK3GrG14H3m3EJIc7mBYRoWDuLkHkM2JyMqtWMcVEY/YqHEWeO4GDUEWb2YoLCGeqhiCuJCwnDqhcBP0vIYbpW4mvJ/SU+YgaU0bD6RcWMqsR8611Lo2eebT6C0P07TdL32zz2scAXgalhXLZuxiOAjQSv218n+E7sAz4qcVeBoo0liajhL5a0Zuk7poHrgRmJba3qKm0YOSHpAcK/40+3c1wsxsWfDqKySPhIHPVp0Oi7CTOuBb4JvBE4HfiAxHMk3lkqi2KQdMgDSO0iYcOQ2CNRA2pRIOBmlAojXzoxfL6MsCbm93svTnckfCROAn4PeKcZN5nxY8CPE/LMPlPiuRJ3e5TukgHEjGmz9u7PcbFhFEVbWd09xsXbGdwYF28E3ixxkxn/E/gDQszMpwMm8RlCsN2S4eB6QkqBeIqBppQ9jHzJ3MNIxLj4YrsNNRsqtFnPOjM2m3G5WRjTmvFuM2oEI+Rr3UZxE8FnYi1wukQZT3TASfphSNQk5iT2S+zIUse4KIx+u4ZHfAl4hpmdmGHfKMbFNe02kjJUeDUhuXDcvnC8GWeY8fNmnGLGk83Ya8bfmvHPZlzuu76VkD7wiYSM9BAie60GbgdeLfGAxDck7qb0kRga0vwwOqlk5DfCdGpRbd9GiGHabJ+zgO8CT+2sDR0FbfD3zwJdCvoQaMHLXgn6IejroM+BzgWtAl0GeiHoLNDJGdrZCpoHbQCd4K/zoK1F/8bllvVe6e5ZGJceRpHUDUvM7GW+3iQe4+INkr6acnwWzgKe6u9fQFiDMU+YYzeCk9SjJJ4kcYHEUYl7JN4n8SmJOyRaTv0q+ELsJPhI3O+vpY9EwZjZVjM7amYP+WtHw9EslAojf5YpDLdVzBKGHwCvJ+SznW23YjN+yYzPEnK6PhtAYpfEbxKS737Z/xi+L/FgV2fhSHxE4lyJ4/y1VBZNyPth9vpShqM5KY2iu0jD0A3rsu2fA74Q+/wkQnAdCAbDbwFPzF6fVoNe6e9fBroI9GvlUGHwNoItaB7YAJzgr/NAz34X4KjXezxwgpdtAI6m7LseuLWb9sZlWrUI1/CIW4GfMrPjJf2IYDz8arsxLszYCOwg9FauNeM4iT+Lff8QpTv1oLETeA3wX4SH9STC8HO3mT3KPye3ExuUN9pWETxqVwBvI/RYP0fK6l7lETV8RKmpf0GAlyHpe2b2DcIPeBvB1nCMMBvxRUl/YWY/B7wFuFLS56NjzXgiIYT+uwjTYR8CXiRxf30745dceJAws1XA2YntHEIGtn8GFgh2nwcIM1A/65+j7QfAvYmytO2Hic+fISilG+TdCPKcuSqoqzbp20ysbAo4kNhvyvebbrcsUc+h9HJt9RmGh/w1ly484UF+ub9/G/AnwDeAnwE+BtxJiAGxwuX6GdDHQAug94J+rIjfaZQ3wnDhKPCQv7b87QEjuLtvAF4FvJfgj/ItQsCkfyCsrH0t8EKCotiYqCN1uNDleWQe9jR6FjK3V8APNQFs9/cHgGrsuwOx91OExVcA064QMpVluUj9nCIEXge8299/jPBv83+A7wDbgZNAjwH9HuhE0KRPjT6q37/PIGydPMwd1N/wIXPF8BPA8wgrjvcQuvkLBAP1Z4H3Ef7ZJ4En4As522mn39eLHtgwirwpKsBsoiyuMGZw/wX/UbZnLUtp61B92TLfhd8F3Qz6DOi7oMeDHg063x/kupuhzXO9EPg7f/9dQkyIdwOngc4BfRhUA30QdFqeN02/bs4u6++HofB5fg9WCfE5/pgwJPgCISTgN4EbgT8kBF9eDzw2TTEMwu/Shjx1z0I7W5E2jHiqxDQqic/dpEpMy3x2FuFfA+BPgc8Tbp4/J4wzVxOGDmcAJ5nxGok/MeO9wL8TMrXfLnHEjBUSDzc+lR2vh2suMDtOcKpg6p9gz78SpkNPBv4ReJU6zLMRm1pbFnTGzJDUM7tGN+34KtwTfDs+8Rp//2bCn8D9wHMJRsAPEgyFjwTi20ldvP8EwW5wH2G8/xVCwOHXAneoR0mo/LqMjG0pF4XRbqpEpWdKq6XUkbWsri7VGz3vINzwNymk2ftHM04m+C7cQ+h6nhfOJ+TE8KxbXwKeDGwiKJUjwMfNWENYiDUP/BpBIa2Gd1wB1z8HHkPIEf0dg0t+Ft5/Brz0E3DqbYSwdI9wI+cjWmwnppS9kTDEeS7hH3kFYTz9B56eYAVwnG/x98nPrb57FkHJX8XSAw5wrZm9heaKYAWhZ/Wgbz9KvEbvVxOGlg8SDHw/8Ncn+u8VGQijbcFf4+Wt3v8j8GpJN0U3g8dOvVADGFJgkBiUVIlpdJM+MQtXkTEUnLS4pgLSHaxeBDyeoEhOl5AZq4HL4IPPCdXOEp5nI9yzr3gcvGKe8DA80MaWtv/phHUfJxEetCgWx1ZC1/phQpc42pp9bvbd3wK/7fVHD7qAw8DzSVcI0fuH5H3iZpjZUeDylId5taTfbHV8FszMf3srwwC2SwFjqCpLsySzsfJJwhhyKla2nYRdImtZlnFbP2ZJYIXgAcEPBN8S3CuoKUyI9OyaHgU2JMp6ao3vVzsMmKFwlDaG2ejZ5wt1qLi2z34YbhQott0oOPvhHp5fPx+y8mEe4i35LICmfFYu1SWh7viiT6CIi9Tftt/yaTjTlcQD/nqm4C2f7vE5jsQsSbnlu8WfBVcW7pKgaVCdS0Ld8UWfQJ8u0q3ALmB9g+9Ty3vX/ls+HXoaKxRee6ss+nku/WxnVNoYpHNJKIwZkLskaBKUOqSPb+OyWrWmEAT1UIPv1+fZuLTzBdLtK+DhN0m3r5B2viDH5tbnWHe/2xmVNvrVTmobZrbeA0jtYrkbQiWxa5pLwjLGRWEs0k7UrbR9uykbxjb61c4wttGvdrptQ42jhtdo7ZKwvD7vpow0ZvYVgrMVBK1aS+ySVtbOvqNcNmjyDFLZoMmTpexxkp4BEA/+6/lJ5tUiC9pYKIySkpJ0PKfqHDAhsbvl/qXCKCkpycrY2TBKSko6Z+wVhplNmVlWd/JO6q+Y2YS3M5NXO97WpG+5tuNt5daGmd1rZgfMbHtebXg70e/SMkVgF/UfM7MjvuV5zab8t5/Oqw0oFQZqsu6lR2wB1kbt5PWDmtkEYZn/QWDCzKp5tONtTRJc/PNis6RNklqOqbvkSv9dVuZ0vVZKWi3pfOBSOgj0nAX/PaLFnfN+L+TCuIToKwyFRXYR7SyOa7edOWDOzCqEm2c+j3b8wcql7hgVM6vmdQ6wqLhv9nb2tDygAxRbnU0IFJXXn9Nh4IiZbfZ2ms50dMPY9zD6hT9oC3n+mE6rOCPdkuuD7KwEFswsl39kZzXBUWnBzGZd0eaCmU3n2ZOVVCP0XvbRRi7fTigVRv+YkrQt70ZcIa3OY1xuZpN9UHhI2uMPQS0v+4JzzNs5QojBkRebcqw7GpIclLSanK9ZqTD6gAcJ2u3vcxlfmtlMzD5So00PvowsuGFtCqjmcS5mNp3nGDzGzbH3FbrNOdqAPHsuMSZ8SApwNfn89kCpMCLtvDZHS/kkMBNZysnvx5wlGLwmgUoe43JJc97DWEn9OoRecT0sRm3LzSjt9Vb8eiVtTb1kJSEqWJ7scUU7CWzJ8VxKx62SkpLsjH0Po6SkJDulwigpKclMqTBKSkoyUyqMkpKSzJQKo6SkJDOla/gQEksUtUDwH8jNvbkNmSruBDUwNEmS1Wj/CmHdT+7OacNK2cMYMtzFfJ17Q+4nKI1KsVIBYZFdW5hZbi7snbhjRwovz4V7w06pMIaPKiGNI7C46CzvtR1Z6MTtPZd1D914V3rvIncX/mGlVBjDx2HgSjPbHv0TxpbOT0T/2mZWde/SSf886XEmFmNmeKyO1HI/JoqxMBWrZ8pjPEya2b6oDoLX5FSaW7fLMunyzUSyEtZwRO7g232f2cjFPfJe9O8qafU0YBrqY1NGbufRuTSh0uL78aUfORnKref5JyYIy+RFeOgqse8OxN5vByZjnxf39Tr2NSr313g6ytnYPgf8+2q87ibyzsSOnUrKGslIeFCPxOSYiZXPNqonpb19KWWTsXaaZvnydiaK/p0HcSt7GEOIwpqOTZKM8E+aOZKTfJyuMJSZbFK+jRAcNuJYbP+VLkPWodB1hHgNs6QMn7RkZNxLCDQDcDFwj/dEqr41racFC8A+P/b6qNB7L8le0c3kGyBoaCkVxpCR7E5L2kH/bu5o4VzDxVQNDIbzBHvFAULvJe24SUK8kDmvowLMuWKak7QpSz2NcEUY2UyujH21WksrPePU2ql/XCgVxvBRia+sjSJsxb6vxYx+60iMx2P2iSqxcX5KeTQsiVhHil0gRqRE0hTGtKSagq1lxtKNkjPADn8fDYsW40h4LyBLPVB/zpM+azLvbdwT+3ql20emYvVVyX+F6VBS+mEMKTGlUWXpQYMw1t9iZvOEm36bmR3Uko/EpJnVCA/jpbHjlpVLqkVGRtyGIClaPr/WH8C470dkrEwdKri884TwgbVoqOHHVAjG3LVevkrSjlj7+LnU1dPg8hyw5SH+Dvuxkd0lGSu0yvIexeqUfUool7ePFWZ2RCEgbabyYcV7SZNq4czmymgzIfLW7lj5jA/1ShKUQ5KSkcN7FpUMuy72zuLTz+QU3XsUKBXGmOAPQjQEaFk+7Eja3SqKmjx2qKTdkg5GNow2Zn/GjnJIUlJSkpmyh1FSUpKZUmGUlJRk5v8DTVCYbNIQI24AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 252x189 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "complexity_axis = [len(bs) for bs in best_subsets]\n",
    "with plt.style.context(['science']):\n",
    "    fig, ax = plt.subplots()\n",
    "    ax2 = ax.twinx()\n",
    "    ax.set_zorder(ax2.get_zorder()+1)\n",
    "    ax.patch.set_visible(False)\n",
    "    \n",
    "    l1, = ax.plot(complexity_axis, last_ubic, 'o-', c='black', markerfacecolor='none', label=f\"$\\lambda = {abs(last_lam)}$\")\n",
    "    ax.set_xticks(complexity_axis)\n",
    "    ax.set_ylabel(\"$\\\\textrm{UBIC}(\\\\xi^{k}, 10^{\\\\lambda})$\", fontsize=12)\n",
    "    ax.set_xlabel(\"Support sizes ($s_{k}$)\", fontsize=12)\n",
    "    ax.annotate('$(u_{xx}, uu_{x})\\checkmark$', fontsize=16, \n",
    "                 xy=(complexity_axis[best_bc], last_ubic[best_bc]), \n",
    "                 xytext=(complexity_axis[best_bc]+0., last_ubic[best_bc]+2000), \n",
    "                 arrowprops={'arrowstyle': '->'})\n",
    "    \n",
    "    l2, = ax2.plot(complexity_axis, b_uns, 'o--', c='blue', markerfacecolor='none', label=\"Uncertainty $\\\\textrm{U}^{k}$\")\n",
    "    s1 = ax2.scatter(complexity_axis[np.argmin(b_uns)], b_uns[np.argmin(b_uns)], c='blue')\n",
    "    ax2.tick_params(axis='y', labelcolor='blue')\n",
    "    if K == 500:\n",
    "        ax.set_title(\"$N_{\\Omega} = 500$, \" + f\"\\%CE = ${round(perf[0], 2)}$\")\n",
    "    elif K == 2000:\n",
    "        ax.set_title(\"$N_{\\Omega} = 2000$, \" + f\"\\%CE = ${round(perf[0], 2)}$\")\n",
    "    elif K == 4000:\n",
    "        ax.set_title(\"$N_{\\Omega} = 4000$, \" + f\"\\%CE = ${round(perf[0], 2)}$\")\n",
    "    elif K == 5000:\n",
    "        ax.set_title(\"$N_{\\Omega} = 5000$, \" + f\"\\%CE = ${round(perf[0], 2)}$\")\n",
    "    elif K == 10000:\n",
    "        ax.set_title(\"$N_{\\Omega} = 10000$, \" + f\"\\%CE = ${round(perf[0], 2)}$\")\n",
    "    elif K == 50000:\n",
    "        ax.set_title(\"$N_{\\Omega} = 50000$, \" + f\"\\%CE = ${round(perf[0], 2)}$\")\n",
    "    else: \n",
    "        assert 0\n",
    "#     ax.legend([l1, l2, s1], [f\"UBIC ($\\lambda = {abs(round(last_lam, 2))}$)\", \"Uncertainty $\\\\textrm{U}^{k}$\", \"Min $\\\\textrm{U}^{k}$\"], \n",
    "#               labelcolor='linecolor', loc='upper center', fontsize=12)\n",
    "    ax.legend([l1, l2, s1], [f\"UBIC ($\\lambda = {abs(round(last_lam, 2))}$)\", \"Uncertainty $\\\\textrm{U}^{k}$\", \"Min $\\\\textrm{U}^{k}$\"], \n",
    "              labelcolor='linecolor', loc='upper center', bbox_to_anchor=(0.385, 1.00), fontsize=12)\n",
    "    plt.savefig(f\"burgers_ubic_samplesize{K}.pdf\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5952fd7a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d918703",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}