Examples/Reaction-Diffusion-2D-big-nPIML-pub.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "823da98c-2a31-4f1c-b696-317d0d110853",
   "metadata": {},
   "source": [
    "## ReactionDiffusion 2D big\n",
    "    - noise_lv = 15\n",
    "    - Now need bruteforcing\n",
    "    - V2 means เอา Odd indices -> Look at issub variable | Odd or even not different\n",
    "    - Choose even"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5c8736f7-3fcd-49e0-9367-6a3627c14b26",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%pylab inline\n",
    "pylab.rcParams['figure.figsize'] = (12, 8)\n",
    "import warnings; warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from sympy import symbols\n",
    "\n",
    "import sys; sys.path.append('../')\n",
    "from PDE_FIND import *\n",
    "\n",
    "import sys; sys.path.append('../../parametric-discovery/')\n",
    "from best_subset import *\n",
    "from frols import frols\n",
    "from p_linear_regression import PLinearRegression\n",
    "from r_pca import R_pca\n",
    "from pde_diff_new import RPCA\n",
    "from RobustPCA.rpca import RobustPCA\n",
    "\n",
    "import scipy.io as sio\n",
    "from scipy.signal import savgol_filter\n",
    "import itertools\n",
    "\n",
    "from sklearn.preprocessing import normalize\n",
    "from sklearn.linear_model import Ridge, BayesianRidge, LinearRegression as SkLinearRegression\n",
    "from sklearn.base import BaseEstimator\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from abess.linear import LinearRegression\n",
    "from bess import PdasLm\n",
    "import pysindy as ps\n",
    "\n",
    "from tqdm import trange\n",
    "from os.path import join as join_path\n",
    "FIGURE_EXPORT_PATH = \"/Users/pongpisit/Documents/figure_export/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f428bf63-f8aa-4de3-97c8-dd17f45a8c55",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nonzero(wei):\n",
    "    return tuple(np.nonzero(wei)[0])\n",
    "\n",
    "def check_percent(bic_scores, complexities):\n",
    "    slope = (bic_scores[1:]-bic_scores[:-1])/(complexities[1:]-complexities[:-1])\n",
    "    slope_index = np.argmin(slope)\n",
    "    percent_improve = 100*np.abs(bic_scores[slope_index+1]-bic_scores[slope_index])/np.abs(bic_scores[slope_index])\n",
    "    percent_from_1 = 100*np.abs(bic_scores[slope_index+1]-bic_scores[0])/np.abs(bic_scores[0])\n",
    "    return slope_index+1, percent_improve, percent_from_1\n",
    "\n",
    "def percent_coeff(pred, ground):\n",
    "    return 100*(np.abs(pred-ground)/np.abs(ground))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0e366b04-cdf3-40e8-9cfc-b0e0c25ea165",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([-10.      ,  -9.921875,  -9.84375 ,  -9.765625,  -9.6875  ,\n",
       "         -9.609375,  -9.53125 ,  -9.453125,  -9.375   ,  -9.296875,\n",
       "         -9.21875 ,  -9.140625,  -9.0625  ,  -8.984375,  -8.90625 ,\n",
       "         -8.828125,  -8.75    ,  -8.671875,  -8.59375 ,  -8.515625,\n",
       "         -8.4375  ,  -8.359375,  -8.28125 ,  -8.203125,  -8.125   ,\n",
       "         -8.046875,  -7.96875 ,  -7.890625,  -7.8125  ,  -7.734375,\n",
       "         -7.65625 ,  -7.578125,  -7.5     ,  -7.421875,  -7.34375 ,\n",
       "         -7.265625,  -7.1875  ,  -7.109375,  -7.03125 ,  -6.953125,\n",
       "         -6.875   ,  -6.796875,  -6.71875 ,  -6.640625,  -6.5625  ,\n",
       "         -6.484375,  -6.40625 ,  -6.328125,  -6.25    ,  -6.171875,\n",
       "         -6.09375 ,  -6.015625,  -5.9375  ,  -5.859375,  -5.78125 ,\n",
       "         -5.703125,  -5.625   ,  -5.546875,  -5.46875 ,  -5.390625,\n",
       "         -5.3125  ,  -5.234375,  -5.15625 ,  -5.078125,  -5.      ,\n",
       "         -4.921875,  -4.84375 ,  -4.765625,  -4.6875  ,  -4.609375,\n",
       "         -4.53125 ,  -4.453125,  -4.375   ,  -4.296875,  -4.21875 ,\n",
       "         -4.140625,  -4.0625  ,  -3.984375,  -3.90625 ,  -3.828125,\n",
       "         -3.75    ,  -3.671875,  -3.59375 ,  -3.515625,  -3.4375  ,\n",
       "         -3.359375,  -3.28125 ,  -3.203125,  -3.125   ,  -3.046875,\n",
       "         -2.96875 ,  -2.890625,  -2.8125  ,  -2.734375,  -2.65625 ,\n",
       "         -2.578125,  -2.5     ,  -2.421875,  -2.34375 ,  -2.265625,\n",
       "         -2.1875  ,  -2.109375,  -2.03125 ,  -1.953125,  -1.875   ,\n",
       "         -1.796875,  -1.71875 ,  -1.640625,  -1.5625  ,  -1.484375,\n",
       "         -1.40625 ,  -1.328125,  -1.25    ,  -1.171875,  -1.09375 ,\n",
       "         -1.015625,  -0.9375  ,  -0.859375,  -0.78125 ,  -0.703125,\n",
       "         -0.625   ,  -0.546875,  -0.46875 ,  -0.390625,  -0.3125  ,\n",
       "         -0.234375,  -0.15625 ,  -0.078125,   0.      ,   0.078125,\n",
       "          0.15625 ,   0.234375,   0.3125  ,   0.390625,   0.46875 ,\n",
       "          0.546875,   0.625   ,   0.703125,   0.78125 ,   0.859375,\n",
       "          0.9375  ,   1.015625,   1.09375 ,   1.171875,   1.25    ,\n",
       "          1.328125,   1.40625 ,   1.484375,   1.5625  ,   1.640625,\n",
       "          1.71875 ,   1.796875,   1.875   ,   1.953125,   2.03125 ,\n",
       "          2.109375,   2.1875  ,   2.265625,   2.34375 ,   2.421875,\n",
       "          2.5     ,   2.578125,   2.65625 ,   2.734375,   2.8125  ,\n",
       "          2.890625,   2.96875 ,   3.046875,   3.125   ,   3.203125,\n",
       "          3.28125 ,   3.359375,   3.4375  ,   3.515625,   3.59375 ,\n",
       "          3.671875,   3.75    ,   3.828125,   3.90625 ,   3.984375,\n",
       "          4.0625  ,   4.140625,   4.21875 ,   4.296875,   4.375   ,\n",
       "          4.453125,   4.53125 ,   4.609375,   4.6875  ,   4.765625,\n",
       "          4.84375 ,   4.921875,   5.      ,   5.078125,   5.15625 ,\n",
       "          5.234375,   5.3125  ,   5.390625,   5.46875 ,   5.546875,\n",
       "          5.625   ,   5.703125,   5.78125 ,   5.859375,   5.9375  ,\n",
       "          6.015625,   6.09375 ,   6.171875,   6.25    ,   6.328125,\n",
       "          6.40625 ,   6.484375,   6.5625  ,   6.640625,   6.71875 ,\n",
       "          6.796875,   6.875   ,   6.953125,   7.03125 ,   7.109375,\n",
       "          7.1875  ,   7.265625,   7.34375 ,   7.421875,   7.5     ,\n",
       "          7.578125,   7.65625 ,   7.734375,   7.8125  ,   7.890625,\n",
       "          7.96875 ,   8.046875,   8.125   ,   8.203125,   8.28125 ,\n",
       "          8.359375,   8.4375  ,   8.515625,   8.59375 ,   8.671875,\n",
       "          8.75    ,   8.828125,   8.90625 ,   8.984375,   9.0625  ,\n",
       "          9.140625,   9.21875 ,   9.296875,   9.375   ,   9.453125,\n",
       "          9.53125 ,   9.609375,   9.6875  ,   9.765625,   9.84375 ,\n",
       "          9.921875]),\n",
       " array([-10.      ,  -9.921875,  -9.84375 ,  -9.765625,  -9.6875  ,\n",
       "         -9.609375,  -9.53125 ,  -9.453125,  -9.375   ,  -9.296875,\n",
       "         -9.21875 ,  -9.140625,  -9.0625  ,  -8.984375,  -8.90625 ,\n",
       "         -8.828125,  -8.75    ,  -8.671875,  -8.59375 ,  -8.515625,\n",
       "         -8.4375  ,  -8.359375,  -8.28125 ,  -8.203125,  -8.125   ,\n",
       "         -8.046875,  -7.96875 ,  -7.890625,  -7.8125  ,  -7.734375,\n",
       "         -7.65625 ,  -7.578125,  -7.5     ,  -7.421875,  -7.34375 ,\n",
       "         -7.265625,  -7.1875  ,  -7.109375,  -7.03125 ,  -6.953125,\n",
       "         -6.875   ,  -6.796875,  -6.71875 ,  -6.640625,  -6.5625  ,\n",
       "         -6.484375,  -6.40625 ,  -6.328125,  -6.25    ,  -6.171875,\n",
       "         -6.09375 ,  -6.015625,  -5.9375  ,  -5.859375,  -5.78125 ,\n",
       "         -5.703125,  -5.625   ,  -5.546875,  -5.46875 ,  -5.390625,\n",
       "         -5.3125  ,  -5.234375,  -5.15625 ,  -5.078125,  -5.      ,\n",
       "         -4.921875,  -4.84375 ,  -4.765625,  -4.6875  ,  -4.609375,\n",
       "         -4.53125 ,  -4.453125,  -4.375   ,  -4.296875,  -4.21875 ,\n",
       "         -4.140625,  -4.0625  ,  -3.984375,  -3.90625 ,  -3.828125,\n",
       "         -3.75    ,  -3.671875,  -3.59375 ,  -3.515625,  -3.4375  ,\n",
       "         -3.359375,  -3.28125 ,  -3.203125,  -3.125   ,  -3.046875,\n",
       "         -2.96875 ,  -2.890625,  -2.8125  ,  -2.734375,  -2.65625 ,\n",
       "         -2.578125,  -2.5     ,  -2.421875,  -2.34375 ,  -2.265625,\n",
       "         -2.1875  ,  -2.109375,  -2.03125 ,  -1.953125,  -1.875   ,\n",
       "         -1.796875,  -1.71875 ,  -1.640625,  -1.5625  ,  -1.484375,\n",
       "         -1.40625 ,  -1.328125,  -1.25    ,  -1.171875,  -1.09375 ,\n",
       "         -1.015625,  -0.9375  ,  -0.859375,  -0.78125 ,  -0.703125,\n",
       "         -0.625   ,  -0.546875,  -0.46875 ,  -0.390625,  -0.3125  ,\n",
       "         -0.234375,  -0.15625 ,  -0.078125,   0.      ,   0.078125,\n",
       "          0.15625 ,   0.234375,   0.3125  ,   0.390625,   0.46875 ,\n",
       "          0.546875,   0.625   ,   0.703125,   0.78125 ,   0.859375,\n",
       "          0.9375  ,   1.015625,   1.09375 ,   1.171875,   1.25    ,\n",
       "          1.328125,   1.40625 ,   1.484375,   1.5625  ,   1.640625,\n",
       "          1.71875 ,   1.796875,   1.875   ,   1.953125,   2.03125 ,\n",
       "          2.109375,   2.1875  ,   2.265625,   2.34375 ,   2.421875,\n",
       "          2.5     ,   2.578125,   2.65625 ,   2.734375,   2.8125  ,\n",
       "          2.890625,   2.96875 ,   3.046875,   3.125   ,   3.203125,\n",
       "          3.28125 ,   3.359375,   3.4375  ,   3.515625,   3.59375 ,\n",
       "          3.671875,   3.75    ,   3.828125,   3.90625 ,   3.984375,\n",
       "          4.0625  ,   4.140625,   4.21875 ,   4.296875,   4.375   ,\n",
       "          4.453125,   4.53125 ,   4.609375,   4.6875  ,   4.765625,\n",
       "          4.84375 ,   4.921875,   5.      ,   5.078125,   5.15625 ,\n",
       "          5.234375,   5.3125  ,   5.390625,   5.46875 ,   5.546875,\n",
       "          5.625   ,   5.703125,   5.78125 ,   5.859375,   5.9375  ,\n",
       "          6.015625,   6.09375 ,   6.171875,   6.25    ,   6.328125,\n",
       "          6.40625 ,   6.484375,   6.5625  ,   6.640625,   6.71875 ,\n",
       "          6.796875,   6.875   ,   6.953125,   7.03125 ,   7.109375,\n",
       "          7.1875  ,   7.265625,   7.34375 ,   7.421875,   7.5     ,\n",
       "          7.578125,   7.65625 ,   7.734375,   7.8125  ,   7.890625,\n",
       "          7.96875 ,   8.046875,   8.125   ,   8.203125,   8.28125 ,\n",
       "          8.359375,   8.4375  ,   8.515625,   8.59375 ,   8.671875,\n",
       "          8.75    ,   8.828125,   8.90625 ,   8.984375,   9.0625  ,\n",
       "          9.140625,   9.21875 ,   9.296875,   9.375   ,   9.453125,\n",
       "          9.53125 ,   9.609375,   9.6875  ,   9.765625,   9.84375 ,\n",
       "          9.921875]),\n",
       " array([ 0.  ,  0.05,  0.1 ,  0.15,  0.2 ,  0.25,  0.3 ,  0.35,  0.4 ,\n",
       "         0.45,  0.5 ,  0.55,  0.6 ,  0.65,  0.7 ,  0.75,  0.8 ,  0.85,\n",
       "         0.9 ,  0.95,  1.  ,  1.05,  1.1 ,  1.15,  1.2 ,  1.25,  1.3 ,\n",
       "         1.35,  1.4 ,  1.45,  1.5 ,  1.55,  1.6 ,  1.65,  1.7 ,  1.75,\n",
       "         1.8 ,  1.85,  1.9 ,  1.95,  2.  ,  2.05,  2.1 ,  2.15,  2.2 ,\n",
       "         2.25,  2.3 ,  2.35,  2.4 ,  2.45,  2.5 ,  2.55,  2.6 ,  2.65,\n",
       "         2.7 ,  2.75,  2.8 ,  2.85,  2.9 ,  2.95,  3.  ,  3.05,  3.1 ,\n",
       "         3.15,  3.2 ,  3.25,  3.3 ,  3.35,  3.4 ,  3.45,  3.5 ,  3.55,\n",
       "         3.6 ,  3.65,  3.7 ,  3.75,  3.8 ,  3.85,  3.9 ,  3.95,  4.  ,\n",
       "         4.05,  4.1 ,  4.15,  4.2 ,  4.25,  4.3 ,  4.35,  4.4 ,  4.45,\n",
       "         4.5 ,  4.55,  4.6 ,  4.65,  4.7 ,  4.75,  4.8 ,  4.85,  4.9 ,\n",
       "         4.95,  5.  ,  5.05,  5.1 ,  5.15,  5.2 ,  5.25,  5.3 ,  5.35,\n",
       "         5.4 ,  5.45,  5.5 ,  5.55,  5.6 ,  5.65,  5.7 ,  5.75,  5.8 ,\n",
       "         5.85,  5.9 ,  5.95,  6.  ,  6.05,  6.1 ,  6.15,  6.2 ,  6.25,\n",
       "         6.3 ,  6.35,  6.4 ,  6.45,  6.5 ,  6.55,  6.6 ,  6.65,  6.7 ,\n",
       "         6.75,  6.8 ,  6.85,  6.9 ,  6.95,  7.  ,  7.05,  7.1 ,  7.15,\n",
       "         7.2 ,  7.25,  7.3 ,  7.35,  7.4 ,  7.45,  7.5 ,  7.55,  7.6 ,\n",
       "         7.65,  7.7 ,  7.75,  7.8 ,  7.85,  7.9 ,  7.95,  8.  ,  8.05,\n",
       "         8.1 ,  8.15,  8.2 ,  8.25,  8.3 ,  8.35,  8.4 ,  8.45,  8.5 ,\n",
       "         8.55,  8.6 ,  8.65,  8.7 ,  8.75,  8.8 ,  8.85,  8.9 ,  8.95,\n",
       "         9.  ,  9.05,  9.1 ,  9.15,  9.2 ,  9.25,  9.3 ,  9.35,  9.4 ,\n",
       "         9.45,  9.5 ,  9.55,  9.6 ,  9.65,  9.7 ,  9.75,  9.8 ,  9.85,\n",
       "         9.9 ,  9.95, 10.  ]))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = sio.loadmat('../Datasets/reaction_diffusion_2d_big.mat')\n",
    "u_sol = real(data['u'])\n",
    "v_sol = real(data['v'])\n",
    "x = real(data['x'][0])\n",
    "y = real(data['y'][0])\n",
    "t = real(data['t'][:,0])\n",
    "\n",
    "n = 512; issub = 2\n",
    "if issub > 1:\n",
    "    spatial_sub_indices = np.array([i for i in range(n) if i%issub==0])\n",
    "    u_sol = u_sol[spatial_sub_indices, :, :][:, spatial_sub_indices, :]\n",
    "    v_sol = v_sol[spatial_sub_indices, :, :][:, spatial_sub_indices, :]\n",
    "    x = x[spatial_sub_indices]\n",
    "    y = y[spatial_sub_indices]\n",
    "m = 201; issub = 1\n",
    "if issub > 1:\n",
    "    time_sub_indices = np.array([i for i in range(m) if i%issub==0])\n",
    "    u_sol = u_sol[:, :, time_sub_indices]\n",
    "    v_sol = v_sol[:, :, time_sub_indices]\n",
    "    t = t[time_sub_indices]\n",
    "\n",
    "dt = t[1] - t[0]\n",
    "dx = x[1] - x[0]\n",
    "dy = y[1] - y[0]\n",
    "\n",
    "# For WeakPDELibrary\n",
    "ground_indices_u = (0, 2, 3, 4, 5, 8, 14)\n",
    "ground_coeff_u = np.array([1.000,-1.000,1.000,-1.000,1.000,0.100,0.100])\n",
    "ground_indices_v = (1, 2, 3, 4, 5, 9, 15)\n",
    "ground_coeff_v = np.array([1.000,-1.000,-1.000,-1.000,-1.000,0.100,0.100])\n",
    "\n",
    "x, y, t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1accd704-7384-41ab-bbbb-b145ab5732aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "u = np.zeros((x.shape[0], y.shape[0], len(t), 2))\n",
    "u[:, :, :, 0] = u_sol\n",
    "u[:, :, :, 1] = v_sol\n",
    "\n",
    "# Odd polynomial terms in (u, v), up to second order derivatives in (u, v)\n",
    "library_functions = [\n",
    "    lambda x: x,\n",
    "    lambda x: x * x * x,\n",
    "    lambda x, y: x * y * y,\n",
    "    lambda x, y: x * x * y,\n",
    "]\n",
    "library_function_names = [\n",
    "    lambda x: x,\n",
    "    lambda x: x + x + x,\n",
    "    lambda x, y: x + y + y,\n",
    "    lambda x, y: x + x + y,\n",
    "]\n",
    "\n",
    "# Need to define the 2D spatial grid before calling the library\n",
    "np.random.seed(100)\n",
    "X, Y, T = np.meshgrid(x, y, t, indexing='ij')\n",
    "XYT = np.transpose([X, Y, T], [1, 2, 3, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1010e2b8-1185-4063-963f-da0475f7ec76",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Need to increase the weak form mesh resolution a bit if data is noisy\n",
    "weak_lib = ps.WeakPDELibrary(\n",
    "    library_functions=library_functions,\n",
    "    function_names=library_function_names,\n",
    "    derivative_order=2,\n",
    "    spatiotemporal_grid=XYT,\n",
    "    K=10000,\n",
    "    is_uniform=True,\n",
    "    periodic=True,\n",
    "    include_interaction=True, \n",
    "    cache=True\n",
    ")\n",
    "\n",
    "# Noiseless\n",
    "u_noisy = u.copy()\n",
    "noise_lv = 15\n",
    "domain_noise = 0.01*np.abs(noise_lv)*np.std(u_noisy)*np.random.randn(u_noisy.shape[0], \n",
    "                                                                     u_noisy.shape[1], \n",
    "                                                                     u_noisy.shape[2], \n",
    "                                                                     u_noisy.shape[3],)\n",
    "u_noisy = u_noisy + domain_noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "bdb4aec0-0bd1-410e-9d35-3523d2a440bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "thres_stlsq = None\n",
    "thres_stlsq = 1e-1\n",
    "if thres_stlsq is not None:\n",
    "    optimizer = ps.STLSQ(threshold=thres_stlsq, alpha=1e-5, normalize_columns=True)\n",
    "\n",
    "# thres = 5e-4\n",
    "# optimizer = L0BNB(max_nonzeros=12, lam=1e-1, threshold=thres, \n",
    "#                   is_normal=False, normalize_columns=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "37acbc57-7a13-496f-b1df-5d011439dd95",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(u)' = 0.849 u + 0.180 v + -0.815 uuu + 0.768 vvv + -0.814 uvv + 0.766 uuv + 0.083 u_22 + 0.084 u_11\n",
      "(v)' = -0.151 u + 0.854 v + -0.798 uuu + -0.817 vvv + -0.796 uvv + -0.818 uuv + 0.086 v_22 + 0.085 v_11 + 0.009 vvvv_22 + -0.010 uvvv_22 + -0.011 uuvv_22 + 0.035 uu_1 + -0.083 uuuu_1 + 0.040 vvvv_1 + -0.089 uvvv_1 + -0.009 vvvv_11 + 0.010 uvvv_11 + 0.011 uuvv_11\n"
     ]
    }
   ],
   "source": [
    "model = ps.SINDy(feature_library=weak_lib, optimizer=optimizer, \n",
    "                 cache=True, feature_names=['u', 'v'])\n",
    "model.fit(u_noisy, ensemble=False, library_ensemble=True)\n",
    "model.print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e1038a14-5752-4f4d-ba5a-541abed6c517",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hasattr threshold in optimizer...\n",
      "u_t = (0.068705 +0.000000i)u\n",
      "    + (1.101967 +0.000000i)v\n",
      "    + (-0.074769 +0.000000i)uuu\n",
      "    + (-0.204172 +0.000000i)vvv\n",
      "    + (-0.073626 +0.000000i)uvv\n",
      "    + (-0.209162 +0.000000i)uuv\n",
      "   \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0, 1, 2, 3, 4, 5)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "_, ensemble_indices_u = model.validate_coef_list(target_index=0, criterion='bic')[0]\n",
    "ensemble_indices_u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ad6a44b6-4d84-443c-b60b-c655de9a0adc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hasattr threshold in optimizer...\n",
      "v_t = (-1.101935 +0.000000i)u\n",
      "    + (0.036918 +0.000000i)v\n",
      "    + (0.204715 +0.000000i)uuu\n",
      "    + (-0.040901 +0.000000i)vvv\n",
      "    + (0.207337 +0.000000i)uvv\n",
      "    + (-0.039198 +0.000000i)uuv\n",
      "   \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0, 1, 2, 3, 4, 5)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "_, ensemble_indices_v = model.validate_coef_list(target_index=1, criterion='bic')[0]\n",
    "ensemble_indices_v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "148a5f2b-5bf4-4e7c-9bf9-6dcc9c740f72",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "stridge_u ...\n",
      "{(0, 1, 2, 3, 4, 5)}\n",
      "{(0, 1, 2, 3, 4, 5, 8, 14), (0, 1, 2, 3, 4, 5), (0, 1, 3, 5, 8, 9, 14, 15, 16, 18, 20, 23, 25, 27, 40, 42, 44, 47, 52, 54, 56, 59, 61, 63)}\n",
      "stridge_v ...\n",
      "{(0, 1, 2, 3, 4, 5)}\n",
      "{(0, 1, 2, 3, 4, 5, 9, 15, 35, 37, 39, 71, 73, 75), (0, 1, 3, 5, 8, 9, 14, 15, 40, 44, 45, 46, 47, 48, 49), (0, 1, 2, 3, 4, 5), (0, 1, 2, 3, 4, 5, 9, 15, 40, 44, 47, 49), (0, 1, 2, 3, 4, 5, 9, 15, 35, 37, 39, 40, 44, 47, 49, 71, 73, 75)}\n"
     ]
    }
   ],
   "source": [
    "# pub\n",
    "# stridge_u ...\n",
    "# {(0, 1, 2, 3, 4, 5)}\n",
    "# {(0, 1, 2, 3, 4, 5, 8, 14), (0, 1, 2, 3, 4, 5), (0, 1, 3, 5, 8, 9, 14, 15, 16, 18, 20, 23, 25, 27, 40, 42, 44, 47, 52, 54, 56, 59, 61, 63)}\n",
    "# stridge_v ...\n",
    "# {(0, 1, 2, 3, 4, 5)}\n",
    "# {(0, 1, 2, 3, 4, 5, 9, 15, 35, 37, 39, 71, 73, 75), (0, 1, 3, 5, 8, 9, 14, 15, 40, 44, 45, 46, 47, 48, 49), (0, 1, 2, 3, 4, 5), (0, 1, 2, 3, 4, 5, 9, 15, 40, 44, 47, 49), (0, 1, 2, 3, 4, 5, 9, 15, 35, 37, 39, 40, 44, 47, 49, 71, 73, 75)}\n",
    "\n",
    "print(\"stridge_u ...\")\n",
    "stridge_u = set()\n",
    "stridge_u.add(ensemble_indices_u)\n",
    "print(stridge_u)\n",
    "for i in range(len(model.get_coef_list())):\n",
    "    stridge_u.add(nonzero(model.get_coef_list()[i][0]))\n",
    "print(stridge_u)\n",
    "\n",
    "print(\"stridge_v ...\")\n",
    "stridge_v = set()\n",
    "stridge_v.add(ensemble_indices_v)\n",
    "print(stridge_v)\n",
    "for i in range(len(model.get_coef_list())):\n",
    "    stridge_v.add(nonzero(model.get_coef_list()[i][1]))\n",
    "print(stridge_v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "83d7409c-12ff-4487-963b-970087cb975e",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_pre, y_pre = model.feature_library.cached_xp_full[0], model.cached_x_dot\n",
    "# nonzero(brute_force(X_pre[:, :16], y_pre[:, 0:1], support_size=7))\n",
    "# nonzero(brute_force(X_pre[:, :16], y_pre[:, 1:2], support_size=7))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "d4de109e-335c-4072-bfa7-c8f23df3b28e",
   "metadata": {},
   "outputs": [],
   "source": [
    "max_complexity = 12; target_index = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "f4dedacd-7ba5-4f87-872b-b545d4d33bf4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{(0,),\n",
       " (0, 2, 3, 4, 5, 8, 9, 14, 15),\n",
       " (0, 2, 4),\n",
       " (0, 2, 4, 8, 9, 14, 15),\n",
       " (0, 2, 4, 8, 14)}"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# pub\n",
    "# {(0,),\n",
    "#  (0, 2, 3, 4, 5, 8, 9, 14, 15),\n",
    "#  (0, 2, 4),\n",
    "#  (0, 2, 4, 8, 9, 14, 15),\n",
    "#  (0, 2, 4, 8, 14)}\n",
    "\n",
    "thres = 1e-4\n",
    "nor_axis = None\n",
    "if nor_axis is not None: Xbnb = normalize(X_pre, axis=nor_axis)\n",
    "else: Xbnb = X_pre.copy()\n",
    "bnb_feature_hierarchy = set()\n",
    "for lam in (1e-1, ): # (1e-1, ) and (1e-2, ) are fine.\n",
    "    bnb_sols = bnb(Xbnb, y_pre[:, target_index:target_index+1].ravel(), max_nonzeros=max_complexity, \n",
    "                   lam=lam, threshold=thres, normalize=True, corrected_coefficients=True)\n",
    "    bnb_feature_hierarchy = bnb_feature_hierarchy.union(set([tuple(np.nonzero(bnb_sols[:, i:i+1].flatten())[0].tolist()) \n",
    "                                                             for i in range(bnb_sols.shape[1])]))\n",
    "bnb_feature_hierarchy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "30b26968-a395-4d07-b510-923bf80de791",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{(0,),\n",
       " (0, 1, 2, 3, 4, 5, 8, 9, 12, 14, 15, 72),\n",
       " (0, 1, 2, 3, 4, 5, 8, 9, 14, 15),\n",
       " (0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 72),\n",
       " (0, 1, 2, 3, 4, 8, 9, 14, 15),\n",
       " (0, 1, 2, 4, 8, 9, 14, 15),\n",
       " (0, 1, 2, 8, 9, 14, 15),\n",
       " (0, 1, 8, 9, 14, 15),\n",
       " (0, 1, 8, 9, 15),\n",
       " (0, 1, 9, 15),\n",
       " (0, 1, 15),\n",
       " (0, 15)}"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nor_axis = 0\n",
    "if nor_axis is not None: Xfrols = normalize(X_pre, axis=nor_axis)\n",
    "else: Xfrols = X_pre.copy()\n",
    "frols_coeffs,selected_inds_list,_,_,_ = frols(y_pre[:, target_index:target_index+1], Xfrols, max_nonzeros=max_complexity)\n",
    "selected_inds_list = set(selected_inds_list)\n",
    "selected_inds_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "977899a8-368b-4936-95b9-97f9780fb639",
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_hierarchy_u = bnb_feature_hierarchy.union(selected_inds_list)\n",
    "if thres_stlsq is not None:\n",
    "    if target_index == 0:\n",
    "        feature_hierarchy_u = stridge_u.union((bnb_feature_hierarchy).union(selected_inds_list))\n",
    "    elif target_index == 1:\n",
    "        feature_hierarchy_u = stridge_v.union((bnb_feature_hierarchy).union(selected_inds_list))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "81da2d5c-7940-4dd1-90c6-72710bd2f387",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RFE_PLinearRegression\n",
      "(0, 2, 4, 8, 9, 14, 15)\n",
      "(0, 4, 8, 9, 14, 15)\n",
      "(0, 8, 9, 14, 15)\n",
      "(0, 8, 9, 15)\n",
      "(0, 9, 15)\n",
      "(0, 15)\n",
      "(0,)\n",
      "--------------------------------------------------\n",
      "(0, 2, 3, 4, 5, 8, 9, 14, 15)\n",
      "(2, 3, 4, 5, 8, 9, 14, 15)\n",
      "(2, 3, 4, 8, 9, 14, 15)\n",
      "(2, 3, 4, 8, 9, 15)\n",
      "(2, 3, 4, 9, 15)\n",
      "(2, 3, 4, 15)\n",
      "(2, 4, 15)\n",
      "(2, 4)\n",
      "(2,)\n",
      "--------------------------------------------------\n",
      "(0, 1, 3, 5, 8, 9, 14, 15, 40, 44, 45, 46, 47, 48, 49)\n",
      "(0, 1, 3, 5, 8, 9, 14, 15, 40, 44, 45, 47, 48, 49)\n",
      "(0, 1, 3, 5, 8, 9, 14, 15, 40, 44, 45, 47, 49)\n",
      "(0, 1, 3, 5, 8, 9, 14, 15, 40, 44, 47, 49)\n",
      "(0, 1, 3, 5, 8, 9, 14, 15, 40, 47, 49)\n",
      "(0, 1, 3, 5, 8, 9, 14, 15, 40, 47)\n",
      "(0, 1, 3, 5, 8, 9, 14, 15, 47)\n",
      "(0, 1, 3, 5, 8, 9, 14, 15)\n",
      "(0, 1, 3, 5, 9, 14, 15)\n",
      "(0, 1, 3, 5, 9, 15)\n",
      "(0, 1, 3, 9, 15)\n",
      "(0, 1, 9, 15)\n",
      "(0, 1, 15)\n",
      "(0, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5)\n",
      "(0, 1, 2, 3, 4)\n",
      "(0, 2, 3, 4)\n",
      "(0, 2, 4)\n",
      "(0, 4)\n",
      "(0,)\n",
      "--------------------------------------------------\n",
      "(0, 1, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 8, 9, 15)\n",
      "(0, 1, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 8, 9, 14, 15)\n",
      "(0, 1, 2, 4, 8, 9, 14, 15)\n",
      "(1, 2, 4, 8, 9, 14, 15)\n",
      "(1, 2, 4, 8, 9, 15)\n",
      "(1, 2, 4, 9, 15)\n",
      "(1, 2, 4, 9)\n",
      "(2, 4, 9)\n",
      "(2, 4)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 12, 14, 15, 72)\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 72)\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15)\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15)\n",
      "(1, 2, 3, 4, 5, 9, 15)\n",
      "(1, 2, 3, 4, 9, 15)\n",
      "(1, 2, 4, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 4, 8, 9, 14, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 8, 9, 14, 15)\n",
      "(0, 1, 8, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 8, 9, 14, 15)\n",
      "(0, 1, 8, 9, 14, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 35, 37, 39, 71, 73, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 37, 39, 71, 73, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 39, 71, 73, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 39, 73, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 39, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 72)\n",
      "--------------------------------------------------\n",
      "(0, 1, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 40, 44, 47, 49)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 40, 47, 49)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 40, 47)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 47)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 2, 4, 8, 14)\n",
      "(0, 2, 4, 8)\n",
      "(0, 2, 4)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 35, 37, 39, 40, 44, 47, 49, 71, 73, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 37, 39, 40, 44, 47, 49, 71, 73, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 39, 40, 44, 47, 49, 71, 73, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 39, 40, 44, 47, 49, 73, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 39, 40, 44, 47, 49, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 40, 44, 47, 49, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 40, 44, 47, 49)\n",
      "--------------------------------------------------\n",
      "(0, 2, 4)\n",
      "--------------------------------------------------\n",
      "SelectKBest\n",
      "(0, 2, 4, 8, 9, 14, 15)\n",
      "(0, 4, 8, 9, 14, 15)\n",
      "(0, 8, 9, 14, 15)\n",
      "(0, 8, 9, 15)\n",
      "(0, 9, 15)\n",
      "(0, 15)\n",
      "(0,)\n",
      "--------------------------------------------------\n",
      "(0, 2, 3, 4, 5, 8, 9, 14, 15)\n",
      "(2, 3, 4, 5, 8, 9, 14, 15)\n",
      "(2, 3, 4, 8, 9, 14, 15)\n",
      "(2, 3, 4, 8, 9, 15)\n",
      "(2, 3, 4, 9, 15)\n",
      "(2, 3, 4, 15)\n",
      "(2, 4, 15)\n",
      "(2, 4)\n",
      "(2,)\n",
      "--------------------------------------------------\n",
      "(0, 1, 3, 5, 8, 9, 14, 15, 40, 44, 45, 46, 47, 48, 49)\n",
      "(0, 1, 3, 5, 8, 9, 14, 15, 40, 44, 45, 47, 48, 49)\n",
      "(0, 1, 3, 5, 8, 9, 14, 15, 40, 44, 45, 47, 49)\n",
      "(0, 1, 3, 5, 8, 9, 14, 15, 40, 44, 47, 49)\n",
      "(0, 1, 3, 5, 8, 9, 14, 15, 40, 47, 49)\n",
      "(0, 1, 3, 5, 8, 9, 14, 15, 40, 47)\n",
      "(0, 1, 3, 5, 8, 9, 14, 15, 47)\n",
      "(0, 1, 3, 5, 8, 9, 14, 15)\n",
      "(0, 1, 3, 5, 9, 14, 15)\n",
      "(0, 1, 3, 5, 9, 15)\n",
      "(0, 1, 3, 9, 15)\n",
      "(0, 1, 9, 15)\n",
      "(0, 1, 15)\n",
      "(0, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5)\n",
      "(0, 1, 2, 3, 4)\n",
      "(0, 2, 3, 4)\n",
      "(0, 2, 4)\n",
      "(0, 4)\n",
      "(0,)\n",
      "--------------------------------------------------\n",
      "(0, 1, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 8, 9, 15)\n",
      "(0, 1, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 8, 9, 14, 15)\n",
      "(0, 1, 2, 4, 8, 9, 14, 15)\n",
      "(1, 2, 4, 8, 9, 14, 15)\n",
      "(1, 2, 4, 8, 9, 15)\n",
      "(1, 2, 4, 9, 15)\n",
      "(1, 2, 4, 9)\n",
      "(2, 4, 9)\n",
      "(2, 4)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 12, 14, 15, 72)\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 72)\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15)\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15)\n",
      "(1, 2, 3, 4, 5, 9, 15)\n",
      "(1, 2, 3, 4, 9, 15)\n",
      "(1, 2, 4, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 4, 8, 9, 14, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 8, 9, 14, 15)\n",
      "(0, 1, 8, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 8, 9, 14, 15)\n",
      "(0, 1, 8, 9, 14, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 35, 37, 39, 71, 73, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 37, 39, 71, 73, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 39, 71, 73, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 39, 73, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 39, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 72)\n",
      "--------------------------------------------------\n",
      "(0, 1, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 40, 44, 47, 49)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 40, 47, 49)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 40, 47)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 47)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 2, 4, 8, 14)\n",
      "(0, 2, 4, 8)\n",
      "(0, 2, 4)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 35, 37, 39, 40, 44, 47, 49, 71, 73, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 37, 39, 40, 44, 47, 49, 71, 73, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 39, 40, 44, 47, 49, 71, 73, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 39, 40, 44, 47, 49, 73, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 39, 40, 44, 47, 49, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 40, 44, 47, 49, 75)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15, 40, 44, 47, 49)\n",
      "--------------------------------------------------\n",
      "(0, 2, 4)\n",
      "--------------------------------------------------\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{1: ((0,), -79320.11373424961), 2: ((0, 15), -79806.76114696664), 3: ((0, 1, 15), -80165.2164044002), 4: ((0, 1, 9, 15), -82462.18168864735), 5: ((0, 1, 8, 9, 15), -82543.64479604733), 6: ((0, 1, 3, 5, 9, 15), -84243.34223096781), 7: ((1, 2, 3, 4, 5, 9, 15), -88659.65261778954), 8: ((0, 1, 2, 3, 4, 5, 9, 15), -88865.28261169813), 9: ((0, 1, 2, 3, 4, 5, 9, 14, 15), -88861.87067724303), 10: ((0, 1, 2, 3, 4, 5, 8, 9, 14, 15), -88856.31070591009), 11: ((0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 72), -88850.51899501913), 12: ((0, 1, 2, 3, 4, 5, 8, 9, 12, 14, 15, 72), -88844.6483186886), 13: ((0, 1, 2, 3, 4, 5, 9, 15, 40, 44, 47, 49, 75), -88833.72763134443), 14: ((0, 1, 2, 3, 4, 5, 9, 15, 39, 40, 44, 47, 49, 75), -88827.29038786446), 15: ((0, 1, 2, 3, 4, 5, 9, 15, 39, 40, 44, 47, 49, 73, 75), -88818.17606672437), 16: ((0, 1, 2, 3, 4, 5, 9, 15, 39, 40, 44, 47, 49, 71, 73, 75), -88809.2414487594), 17: ((0, 1, 2, 3, 4, 5, 9, 15, 37, 39, 40, 44, 47, 49, 71, 73, 75), -88800.62358299637), 18: ((0, 1, 2, 3, 4, 5, 9, 15, 35, 37, 39, 40, 44, 47, 49, 71, 73, 75), -88799.60664017875)}"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "br1 = backward_refinement(feature_hierarchy_u, (X_pre,y_pre[:, target_index:target_index+1]), mode='rfe', ic_type='bic', verbose=True)\n",
    "br2 = backward_refinement(feature_hierarchy_u, (X_pre,y_pre[:, target_index:target_index+1]), mode='k best', ic_type='bic', verbose=True)\n",
    "br_u=br1+br2; br_u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "3a723af4-722a-4f8b-990c-8188d2c6defd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ploting less than max_complexity\n",
      "Delete 13 complexity\n",
      "Delete 14 complexity\n",
      "Delete 15 complexity\n",
      "Delete 16 complexity\n",
      "Delete 17 complexity\n",
      "Delete 18 complexity\n"
     ]
    }
   ],
   "source": [
    "print(\"Ploting less than max_complexity\")\n",
    "for c in range(max_complexity+1, max(list(br_u.track.keys()))+1):\n",
    "    del_c = br_u.delete(c)\n",
    "    if del_c is not None: print(f\"Delete {c} complexity\")\n",
    "    else: print(f\"{c} not exist\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "9aa75fe0-d6c5-4f23-883a-c27488f6d1d3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0,) -79320.11373424961\n",
      "(0, 15) -79806.76114696664\n",
      "(0, 1, 15) -80165.2164044002\n",
      "(0, 1, 9, 15) -82462.18168864735\n",
      "(0, 1, 8, 9, 15) -82543.64479604733\n",
      "(0, 1, 3, 5, 9, 15) -84243.34223096781\n",
      "(1, 2, 3, 4, 5, 9, 15) -88659.65261778954\n",
      "(0, 1, 2, 3, 4, 5, 9, 15) -88865.28261169813\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15) -88861.87067724303\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15) -88856.31070591009\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 72) -88850.51899501913\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 12, 14, 15, 72) -88844.6483186886\n"
     ]
    }
   ],
   "source": [
    "# For reproducibility of 2% noise experiment\n",
    "# u\n",
    "# (1,) -90432.11492322256\n",
    "# (1, 14) -90900.56314841735\n",
    "# (0, 1, 14) -91096.59801155927\n",
    "# (0, 1, 8, 14) -92631.07213881976\n",
    "# (0, 1, 8, 14, 15) -92711.06074776711\n",
    "# (0, 1, 8, 9, 14, 15) -92912.16926463085\n",
    "# (0, 2, 3, 4, 5, 8, 14) -99273.17863634705\n",
    "# (0, 1, 2, 3, 4, 5, 8, 14) -99629.82290987358\n",
    "# (0, 1, 2, 3, 4, 5, 8, 14, 67) -99627.38471701351\n",
    "# (0, 1, 2, 3, 4, 5, 8, 14, 67, 72) -99626.97018095759\n",
    "# (0, 1, 2, 3, 4, 5, 8, 9, 14, 67, 72) -99619.72292591527\n",
    "# (0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 67, 72) -99611.13907227767\n",
    "# + stlsq before\n",
    "# (1,) -90432.11492322256\n",
    "# (1, 14) -90900.56314841735\n",
    "# (0, 1, 14) -91096.59801155927\n",
    "# (0, 1, 8, 14) -92631.07213881976\n",
    "# (0, 1, 8, 14, 15) -92711.06074776711\n",
    "# (0, 1, 2, 4, 8, 14) -94460.72867681735\n",
    "# (0, 2, 3, 4, 5, 8, 14) -99273.17863634705\n",
    "# (0, 1, 2, 3, 4, 5, 8, 14) -99629.82290987358\n",
    "# (0, 1, 2, 3, 4, 5, 8, 14, 67) -99627.38471701351\n",
    "# (0, 1, 2, 3, 4, 5, 8, 14, 67, 72) -99626.97018095759\n",
    "# (0, 1, 2, 3, 4, 5, 8, 9, 14, 67, 72) -99619.72292591527\n",
    "# (0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 67, 72) -99611.13907227767\n",
    "\n",
    "# v\n",
    "# (0,) -89704.75044145095\n",
    "# (0, 15) -90094.80085394751\n",
    "# (0, 1, 15) -90433.5949355086\n",
    "# (0, 1, 9, 15) -92311.22556061689\n",
    "# (0, 1, 8, 9, 15) -92462.0683253755\n",
    "# (0, 1, 3, 5, 9, 15) -94229.07163281288\n",
    "# (1, 2, 3, 4, 5, 9, 15) -98884.21234254984\n",
    "# (0, 1, 2, 3, 4, 5, 9, 15) -99223.00785465668\n",
    "# (0, 1, 2, 3, 4, 5, 9, 14, 15) -99226.63012262661\n",
    "# (0, 1, 2, 3, 4, 5, 9, 14, 15, 69) -99230.05524578446\n",
    "# (0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 69) -99224.97738853545\n",
    "# (0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 69, 72) -99216.5636109259\n",
    "# + stlsq before\n",
    "# (0,) -89704.75044145095\n",
    "# (0, 15) -90094.80085394751\n",
    "# (0, 1, 15) -90433.5949355086\n",
    "# (0, 1, 9, 15) -92311.22556061689\n",
    "# (0, 1, 8, 9, 15) -92462.0683253755\n",
    "# (0, 1, 3, 5, 9, 15) -94229.07163281288\n",
    "# (1, 2, 3, 4, 5, 9, 15) -98884.21234254984\n",
    "# (0, 1, 2, 3, 4, 5, 9, 15) -99223.00785465668\n",
    "# (0, 1, 2, 3, 4, 5, 9, 14, 15) -99226.63012262661\n",
    "# (0, 1, 2, 3, 4, 5, 9, 14, 15, 69) -99230.05524578446\n",
    "# (0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 69) -99224.97738853545\n",
    "# (0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 69, 72) -99216.5636109259\n",
    "\n",
    "complexities = np.array(sorted(list(br_u.track.keys())))\n",
    "bic_scores = []\n",
    "eff_feature_names = []\n",
    "for com in complexities:\n",
    "    indices, bic = br_u.track[com]\n",
    "    bic_scores.append(bic)\n",
    "    print(indices, bic)\n",
    "bic_scores = np.array(bic_scores)\n",
    "min_bic_score = np.min(bic_scores)\n",
    "div_scale = 1e3\n",
    "rel_bic_scores = (bic_scores-min_bic_score)/div_scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "c8de7ab2-ca6d-422a-8ff3-0fa870c89efc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUoAAADyCAYAAAArx4ypAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAqS0lEQVR4nO3df1Rb9f0/8OcN1Vl/NCHt5g96JtwUW9eeKQn4a7pOCXWnDbhJAKlnR1aEtN12dNMSul/+2CwmVY+eOVtCoW5zMpK0x0PbaUno1OlxgybUH9tEyCXdyqZbGwKdZ2pL7vePfnM/BPITktzc5PU4Jwdy7829r6vk1fe99/1+vRme53kQQgiJSCZ2AIQQkukoURJCSAyUKAkhJAZKlIQQEkPaEqXb7YZGo5mznOM4mM1m2O12mM1m+P3+dIVECCFxYdLx1Ntut4NlWWg0Gsw+nEajgcvlAnAuaRqNRthstlSHRAghcUtLohQOxjAhiZLjONTU1AiJEgDy8/MxMTGRrpAIISQmUe9ROp1OKJXKkGVKpRJut1ukiAghZK5FYh480v1In88X8TPLly/H2bNnhfcFBQUoKCgAAIyPjwu/xyPV22diTLS9+Meg7cU/xsztx8fHMT4+LqxbtGgRTpw4EfoBPo1mH85kMvFarTZkGcuyvM1mi7iPysrKea0TY/t0HIO2T+726TgGbS/+MRLNI6JeeisUijmtR5/PB4VCMa/91dfXp3T7+Uh1TKk+h1yLfz6fybRzkHr88/lMyr/LCaXhBZp9OI/Hw6vV6pBlCoWCn5iYiLiP+fzrlCxiHjtZpH4OUo+f56V/DlKPn+cl0KKceV+SZdmQdRzHobS0dN4tylRLRws01aR+DlKPH5D+OUg9fiDxc0hL9yCn0wmHwwGz2YyWlhaUlZVBr9cDOJcc29vbUVZWhsHBQWzfvj1qoqyqqkJvb2/IspGREXR1dcHr9aKwsBCbNm1CcXFxKk+JEJKlwuWYtPajTIbZJ7F3717ce++9Qh/N4M/Ozk40NDSIFyghRJLCJUpJj/UeGRnBvffei0AggOnp6ZCfjY2NGB0dFTtEQkgWkHSi7OrqAsMwEdd3dnamMRpCSLaSXKIcHx9HVVUVuru74fV654wdDwoEAujt7cXw8HCaIySESFF3dzeqqqpCOp8HSS5RFhQUoLe3F/X19SgsLIzYomQYBh6PB6tWrcINN9yAXbt2RR3xQwjJbfX19ejt7Q07wkdyiXKmTZs2RWxRMgwDl8sFq9WKZcuW4Xvf+x4uv/xyVFdXo7e3F2fOnElztIQQqZJ0oiwuLkZnZydkMhny8vJCfnZ2dmL16tWoqanBwYMHceLECTz++OPgOA533HEHrrjiCtx3331wu90Rky0hhABZ0D0IAEZHR9HZ2Sn0o2xsbMSKFSsi7uOdd97Br371K/z2t7/FRx99hNWrV+Oee+7B3XffjSuuuCLVp0AIyWBZ2Y9yIc6ePYu+vj78+te/xksvvYQzZ86goqIC99xzD+644w5ceOGFAKhDOyG5JFyOEbXMmtgWLVqE9evXY/369fD7/bDZbPjVr36FjRs34pJLLkFtbS2WLVsGs9kMmUwmdGg3m83UoZ2QHJLTLcpIRkdH8Zvf/AZdXV1z69L9fzKZDMPDw1Ev8Qkh0pN1I3NSZcWKFXjkkUdw9913QyaL/J/oueeeS2NUhBCxUKKM4vjx4xHXBQIBPP3001i7di0ee+wxDA4OYnp6Oo3REULSRXKJcubInFSL1qFdJpOhoqIC+fn5MJlMuO666/D5z38etbW16OjoiJpkCSGZJ9rIHLpHGcXIyAhWrVqFQCAwZ93Me5RnzpzBn//8ZzgcDvT19WFgYACBQABXXXUV1q1bh4qKCnzta1/DkiVLoh6LnqwTIr6wOSZ5NYPTI93Vlffu3cvLZDI+Ly8v5OfevXsjfsbn8/H79u3jDQYDX1RUxAPgFy1axN9yyy38o48+yv/pT3/iz549K2zf1dWV8DEIIakRLsdQizIOiXZon83j8aCvrw99fX04cuQIpqamoFAoUF5ejmuuuQYPPfRQ2NFB9GSdkPSjDucZ4OzZsxgYGEBfXx8cDgfeeuutiEMo8/LysG3bNrS1taU5SkJyF3UPygCLFi3CTTfdhIcffhhvvvkmvvnNb0Z8YMTzPLxeb3oDJITMQYlSZFdddVXEvpoMw6CwsDC9ARFC5qBEKbJopeJ4nkdjY2OaIyKEzEaJUmSzS8UFL8MZhkFnZyc9yCEkA0juYY5Go0FBQQHq6+uzYn7hoJlP1t99913861//AsdxkMvlYodGSE7o7u5Gd3c3xsfH4XK5QtZJLlFK/al3PE6cOIGVK1di69at2Llzp9jhEJJT6Km3RCxfvhxGoxHPPPMMTblLSAagRJmhHnzwQVx22WV48MEHxQ6FkJyXEYV7OY6D0+mEUqkEx3HQ6/VgWVbssER14YUXwmw2o76+Hv39/SgvLxc7JEJyVkYkSrvdjpaWFuG9wWBAe3u7iBFlhrq6Ojz77LO4//77MTQ0hEWLMuJ/FyE5JyMuvXt6esQOISMxDIOnn34a7733Hjo6OsQOh5CclRGJUqlUQqPRCJfgFRUVYoeUMUpLS3HPPffgJz/5Cfx+v9jhEJKTYiZKr9eLnTt3Yt26dSgrKxNet99+O5544omkjEW22WwAAJVKBZvNBr1ev+B9ZpMdO3bgk08+waOPPip2KITkpKj9KFtbW8EwDGpra1FSUjJn/dDQEHp6esAwzIIq3NjtdigUCnAcB4PBgObm5oj3KIMdzoOyreN5JDt27MBDDz2E9957DytXrhQ7HEKyRrCjeVC4DucRC/eazWbe7/fHVejS7/fzra2tCZTG/D8ej4dvaWkJea9QKHiPxxN2+3QX7s0U//vf//jCwkJ+w4YNYodCSFYLl2MiXnpv27Yt7uFzcrl83i1Kt9uNsrIy4T3Lsti+fTvdj5vlggsuwM6dO3Ho0CEcPnxY7HAIySmiP8xRq9UYHBwMWXbq1Cmo1WqRIspc1dXVuOWWW/D9738fZ8+eFTscQnLGnEQ5NjaGLVu24NixY/B6vVi3bh2Ki4tx7NixlATAsiwqKipgNpthsVhgsVhgMBhSciypYxgGzzzzDN5//33s3r1b7HAIyRlzejC3t7dDq9Vi9+7dwlSscrkcHR0duPbaa1MShFarhVarTcm+s01JSQk2bdqEhx56CBs3boRSqRQ7JEKy3pwWZUVFBaqrq7F7925otVqUlJSAZVlKZBnksccew5kzZ/Dwww+LHQohOWFOovT7/di+fTsACOOLn3jiCUxOTqY3MhLRpZdeih//+Md47rnn8Ne//lXscAjJenMSZXV1NVpbW0OWyeVyalFmmPvuuw9XXnklHnjgAbFDISTrhX3qPbtbUFlZGY4ePYr9+/cL81KLZXx8HFVVVSEdRHPR5z73OTzxxBN45ZVX8Pvf/17scAiRvO7ublRVVWF8fHzOuogjcyYnJ9HU1ASn04mlS5cKE2D5fD4wDIPS0lLYbDYsWbIktdHPkgsVzuPF8zzKy8vxz3/+E++++y7OO+88sUMiRPISqnDe0dEBk8kEn8+HkZERjI6OYnR0FD6fD6dOncKuXbsWNGyRLFywutDIyAh++ctfih0OIVkrYqIsKipCUVFRxA+yLIvS0tKUBEXi9+UvfxlNTU145JFHcPLkSbHDISQrRUyUHMfF/PDsETVEHI8++igCgQAeeughsUMhJCtFLJmt1+tRWloKhmFCOjX7fD7h92B5NCKuL3zhC/jpT3+KlpYWbNmyBWvWrBE7JEKySszpasfGxuB0OoUiFQqFAlqtNupleSrRw5zwPvvsM6xZswZXXnkl+vr6wDCM2CERIknhckzMSViKiorQ1NSUsqBIcpx//vl48sknUVVVhQMHDqCqqkrskAjJGguqHrR///5kxUGSQKfToaKiAg888AA+/fRTscMhJGvMO1FOTU2J8jCHOpxHxjAMnnrqKXAch2effVbscAiRlHl1OC8tLcXQ0FDEnfI8D4ZhMD09nbxI40D3KGP7zne+gxdeeAEjIyP4whe+IHY4hEhKwh3Ojx49iunp6bCvQCCAbdu2pTxokrhHHnkEMpkMP/nJT8QOhZCsEDFRlpSUYGxsLOqHaVrZzLRs2TI8/PDD2LNnD95++22xwyFE8qLeo7zzzjujfjhYho1knq1bt6K4uBj3338/YvQAI4TEIPqcOSQ1zjvvPDz11FN49dVX8dJLL4kdDiGSRokyi61fvx5f//rX8eCDD1J3IUIWIK5EuWfPnoSWk8zx1FNP4fjx4/jxj3+M7du3o76+Htu3b8fIyIjYoREiGXElSo1Ggy1btoQs27JlC1UPkoCrr74at956K5544gmYzWZYrVbs3LkTq1atwvPPPy92eIRIQlyJsqSkRCi4AACbN29Gc3NzymZlJMkzMjKC/v5+AEAgEEAgEBC6dzU2NmJ0dDRpx6EWK8lafALcbjevVCr5oaGhRD6WVGq1mq+srORffPFF0WKQktbWVj4vL48HMOclk8l4o9G44GN0dXXxMpmMz8vLC/m5d+/ehZ8AIWny4osv8pWVlbxarZ6zLmb1oJlqa2thMBhgsVjQ09OTotQdHY3MSUx9fT2sVisCgUDY9TKZDEuXLoVCoQh5yeXyOcvCrfvnP/+Jq6++Ouz+ZTIZhoeHsWLFilSfJiFJM6/qQcC5+XOam5vR0dEBuVwOjUaDuro6dHR0pH3OHJKYwsLCiCXXGIZBeXk51q5dC7/fj8nJSfj9fvj9fni9XuF3v9+PM2fORNxHpH9rGYZBZ2cnTRlCJC+uRNnW1gaLxSLMzqhQKGCxWNDW1kZfggy3adMmmM3msOsYhsFzzz0Xs8XH8zw++eSTkMQZfD311FNwuVxhkyXP8/B6vck4DUJEFVeifPzxx+csk8vlSU2STqcTHMeBZVkAoHnEk6S4uBidnZ1obGwUWn/Bn52dnXFdFjMMg8WLF2Px4sW4/PLLQ9a98847GBoaClschWEYFBYWJutUCBFNXIky1ZxOJ2w2G9rb28FxHCoqKuDxeMQOK2s0NDTg5ptvRmdnJ7xeLwoLC9HY2JiUe4fRWqzBJ+uESF3YRDk5OYmJiYmEWgNTU1MAMK97lgaDAS6XC8C52R0dDkfC+yDRrVixIiW3ScK1WIFzSZJlWXzxi19M+jEJSbew/SjlcjkcDkfcFcz37dsHq9U6ryTJcRx8Ph8UCgXcbjf8fr9w+U2koaGhAcPDw9i2bRtqa2vR0tICq9WKv//972hpaRE7PEIWLOKld1NTE4aGhlBbWwuVSoWysjKwLAuFQgG/3w+O4zAwMICxsTEYDAZUV1fPKwC32w2lUgm73Q6tVguLxQKWZaHX6+d9UiT9wrVYP/roI3zve9/DDTfcgLvuukukyAhZuLj6UU5OTsJqtcLj8cDv90OhUEClUiVlNkaLxQKDwYCJiQkhCefn50fscqLRaFBQUCC8r6+vR319/YJiIKnB8zzuvvtu9Pb2YmBgAF/60pfEDomQObq7u0OmlhkfHxduBQYl1OE8FZxOJ2pqajAxMSEsYxgGLpcLarV6zvbU4VxaPv74Y1x//fWYnp7GwMAALrnkErFDIiSqhKaCSBe6H5ndLrroIuzbtw8nTpxAY2MjFREmkpQRibK0tBR+vx8AhL6U4VqTRJpWrlyJvXv3wmaz4ZlnnhE7HEISlhH9KG02G4xGIzQaDVwuF3UPykJ6vR4PPPAAtm3bhtLSUtx8881ih0RI3DIiUSoUCrS3t4sdBkmxtrY2DAwMoLa2Fm63G5dddpnYIRESl4iX3l6vN+QVFOxYTkiizjvvPPT09IDnedx11104e/as2CEREpeIidLhcIBlWZhMJnAcJyw/deoU9u3bF3dndEJmuvzyy9HT04M33ngDP/zhD8UOh5C4RLz0ZlkWLpcLJSUlIcuLiopQVFSEyclJ7N+/P+aUtoTM9tWvfhUmkwkPPvggbrjhBvobIhkvYotycnJyTpKcSS6Xi9LVY3x8HFVVVSEdRIn0/OAHP0B1dTUaGhrwwQcfiB0OIeju7kZVVRXGx8fnrIuYKH0+X8wdz+wkni4FBQXo7e2l0TgSxzAMurq6cPnll6O6uhoff/yx2CGRHFdfX4/e3t6QkX9BERNlPGXOqBQaWYglS5Zg//794DgOmzdvps7oJGNFTJQsy0adt3vPnj1QKpUpCYrkjtWrV6OjowMvvPACdu3aJXY4hIQVtXrQunXrYLPZsHnzZqH4BcdxQp/Hw4cPpydKktU2btyIt956C/fffz80Gg2uv/56sUMiJETUIYx9fX0oLy9HY2MjNBoN1Go19Ho91Go1JUmSVE8++SQ0Gg30ej3+85//iB0OISFijvVuaWmBz+fD6OgoPB4PAoEATShGku7888+HzWbDJ598go0bN4adg4cQsURMlEeOHAl5H+w/OdOxY8dSEhTJTcuXL8fvfvc7HDlyBA8//LDY4RAiiDoy5/Tp05iamor46unpSWesJAeUl5fj5z//OX7+85/j4MGDYodDCIAohXtlMhkYhon4weC0p+m+RApWOKfK5tkrEAjgG9/4Bv74xz/C5XJRzVKSFsFK5wlVON+8eTOMRmPELkCnTp2C2WzG7t27kx9xFFThPDf4/X5oNBrI5XK8+eabWLx4sdghkRyRUIXzmpoaFBUVQS6Xh32xLIuampqUB01yk0KhwL59+/C3v/0N3/3ud8UOh+S4iP0oy8vLY344nm0Ima9rr70Wu3btwre//W2oVCqcPn0aXq8XhYWF2LRpE4qLi8UOkeSIqN2Djhw5gi1btuAPf/gDgHM1KtetW4elS5di69ataQmQ5LaGhgasXbsWP/rRj2A2m2G1WrFz506sWrUKzz//vNjhkRwRMVHu27cPWq0Wg4OD0Ov1OHbsGJqbm6FQKNDa2orR0VFs2bIlnbGSHDQyMoLXX38dwLmHPIFAANPT0wgEAmhsbMTo6KjIEZJcEDFRWq1WTExM4OjRo/B4PNixYweMRiOsViu2bduGvr4+nDp1Kp2xkhzU1dUFmSz8nynDMOjs7ExzRCQXRUyUpaWlkMvlAM7dWK+oqJhzT7KsrCy10ZGc5/V6I1YV4nk+ZJoSQlIlYqLMz88PeR+uL9vsbQhJtsLCwoj9eRmGQWFhYXoDIjkpaj3K06dPC+8nJyfDbpNuVOE8t2zatClqi7KxsTHNEZFsFa3Cedwjc4IjcWa/T/fIHOpwnnuef/55NDY2gmEYBAIB4W+vq6sLDQ0NYodHskxCHc6bm5vh8/mE18TERMh7n8+HpqamlAdNSENDA4aHh7Ft2zbccccdAIAdO3ZQkiRpE7HDucFgEB7mRNuGkHRYsWKFUN7vpptuwsDAgMgRkVwSsUUZbQbGRLYhJNkqKyvR19eHTz/9VOxQSI6IWbg33YxGI/x+v9hhkAym0+nw8ccf47XXXhM7FJIjMipRut1umM1mscMgGW7NmjX44he/iAMHDogdCskRGZUoOY6j2oMkJoZhUFlZiYMHD9IUtyQtMiZR2u126PV6scMgEqHT6eD1evHXv/5V7FBIDkgoUe7cuRN1dXUAgP7+fkxNTSUlCL/fD4VCkZR9kdzwta99DRdddBFdfpO0iDtRtra2QqFQQKvVAjhXi9LpdCYlCKvVKuw3luDInOCLRujkpgsuuAAVFRU0rw5ZsOCInOAr3MiciP0oZysrK0N1dTX6+/uTGqTT6URtbW3c2xcUFNDIHALg3OV3c3MzTp48iWXLlokdDpGo2fNvVVVVzdkm7hbl2NgYAIQMYxwcHFxIfAKr1QqLxQKLxQKO49DW1ga3252UfZPstWHDBgQCAbz88stih0KyXNwtypKSEpSWlmLp0qVwOBxwOp0wmUwLDmD2JbfBYIDBYKCn3ySmyy67DGVlZTh48CC+9a1viR0OyWJxtyjLy8thtVpRUlICnudhsVhw2223JS0Qv98v9KE0mUzUoiRx0el0eOWVV/DZZ5+JHQrJYhGrB822ZcsW7Nq1K9XxxETVg8hMQ0NDUKvV6O/vT+o/3CR3JVQ9aDaHw4E9e/bgyJEjSQ+MkPm69tprccUVV9DTb5JScd+jdLlckMvlmJycxJ49e6BUKqHVarFkyZJUxkdIVAzDQKfT4cCBA3jyyScjVkMnZCHiblEGS67J5XLwPI+WlhaqR0kyQmVlJUZHR/HBBx+IHQrJUnEnyrq6OmzZsgVKpRJutxsOhwM9PT2pjC0smgqCzHbbbbfhggsuoMtvsiDzmgpithUrVqC1tRX33ntv0gNMBD3MIeFUVlbi9OnTePXVV8UOhUjcgh7mmEwm0ZMkIZHodDq88cYbmJiYEDsUkoXiTpTV1dVzltGcyiRTbNiwAdPT0zh8+LDYoZAsFPGp9/79+0Oeau/Zsydkvd/vh8PhoD9MkhGWL1+OkpISHDhwAHfddZfY4ZAsE7FFuWPHDhw9elR4v3v3bkxMTAgvnudx6tSptARJSDx0Oh1efvllnD17VuxQSJaJ2KKcmSQBoKOjY85kYvGWRiMkHXQ6HX72s5/hrbfewi233CJ2OCSLxH2PMj8/X/h9cnIS+/btC1lGiNhKS0tx6aWXUjFfknRxJ8qZRXrlcjmqq6uTVriXkGSQyWTYsGED9ackSRd1COPk5CSsVisYhoHD4Ziz3uVypb3LULDD+exim4QA5y6/u7q64PF4oFKpxA6HSEh3dze6u7vDdjiP2qKUy+XQarU4evQoPB4PRkdHQ14tLS0pCzqSYIVzSpIknIqKCpx//vnUqiQJq6+vR29vLwoKCuasi1kUo6ioCLt370Z/fz/Ky8tTEiAhyXLxxRfj1ltvxcGDB3HfffeJHQ7JEgkV7p3pyJEj2L9/f9IDImShKisr8dprryVtllBC4i6zBpzrhM5xHACA53kcPXoUd955Z0oCI2S+NmzYgO9+97vo6+ujueJJUsSdKFtbW+H3++Hz+cCyLPx+PwwGQypjI2ReCgsLsWbNGhw8eJASJUmKuBOlSqVCU1MTxsbGwDAMCgsLqdo5yViVlZXo6OjA9PQ08vLyxA6HSFzc9yhZlsXx48dRVFQEu92eypgIWTCdToeTJ09iYGBA7FBIFog7Ufr9frAsi6mpKZw8eRK333472tvbUxkbIfN2/fXXY9myZdRNiCRF3IV7Z+vv70dpaakwRUS6aDQaFBQUUIdzEtM999yDoaEhvPPOO2KHQiRgZodzl8sVsi7uFuVs5eXlGBsbW3BwiaIO5yReOp0O7777Lo4fPy52KEQC5tXhfHb9ydkmJiZgtVoxODi48AgJSYF169Zh0aJFOHToELZu3Sp2OETCIrYoZ9efnP0CzvWlJCRTyeVyrF27lqoJkQWL2KI0mUwxhyxSPUqS6XQ6HYxGI/773//i4osvFjscIlERW5ThkuTOnTtRV1cH4NzDnGRVZ3G73TCbzTCbzaipqYHf70/KfgmprKzEZ599RiUByYLE/TCntbUVCoVCaEWWl5cn7Y/P6XSipaUFLS0tKCsro+IbJGlUKhVWrVpF3YTIgsSdKMvKytDU1ASWZZMagNvtRltbm/Ber9fD7XYLY8oJWSidTodDhw4hEAiIHQqRqLgTZbArEMMwwrJkPPFWq9Xo6OgQ3gcvu5VK5YL3TQhw7vL7ww8/nNM3jpB4xT3Wu6SkBKWlpVi6dCkcDgecTidMJlNSgphZuKCnpwdarRYKhSLstsEK50HU8ZzEctNNN0GhUODgwYMoKysTOxySYYIdzYPCVThPaGTO2NiYMGyxrq5uzqyMC+X3+6HRaOByuSImyqqqKvT29ib1uCT7bdy4EcPDw9SqJDGFyzEJjcwpKirC448/jscffxxA7E7piTIajXA4HBGTJCHzVVlZCbfbHba1QEgsURPl1NQU9u/fD6/XO2ddSUlJUrvxmM1mGI1GodYldREiyXT77bcjLy8Phw4dEjsUIkERE+XY2BgKCwuh1+uhUqng9XoxNTWFzZs3o7i4GHl5eUkbvmi326FWq4UkabVaqVVJkkqpVOIrX/kKdRMi8xIxUZpMJthsNgQCAQwODmL37t1C/8bq6mrs3r0bPT09Cw6A4zjU1NSgoqICDMMgPz8fRqNxwfslZLbKyko4nU7873//EzsUIjERn3qrVCohMarVavT398NisST9AQ7LsjRmnKSFTqfDtm3bcOTIEWzYsEHscIiERGxRzuwvCZxLlrOT5LFjx1ISFCGpsHLlSqxYsYIuv0nCIibKU6dO4fTp05iamsLU1BQmJyeF34MvGj9LpIRhGOh0Ohw8eJCuYkhCot6jVCgUyM/PR35+PvR6vfB7cDnVoiRSo9PpcOLECbz99ttih0IkJGKibG5uhs/nE14TExMhv586dQpFRUXpjBXA/43MmdmTnpB43XLLLViyZAnVqCRzdHd3o6qqKrGROUNDQzEf3MSzTbLRyByyULW1tTh+/Dj+/Oc/ix0KyUAJjcyJJwGmO0kSkgw6nQ4DAwP46KOPxA6FSMS8JxcjRKrWr18PhmFolA6JGyVKknOWLVuGG2+8kboJkbhRoiQ5qbKyEn19ffj000/FDoVIACVKkpN0Oh0+/vhjvPrqq2KHQiSAEiXJSatXr8aVV15Jl98kLpQoSU5iGAaVlZU0SofERXKJkjqck2TR6XTwer34y1/+InYoJANE63AuuURZUFCA3t5emieHLNjatWtx0UUX0eU3AXBu/q3e3l4UFBTMWSe5RElIslxwwQVYt24dJUoSEyVKktN0Oh3eeustnDx5UuxQSAajREly2vr16xEIBPDyyy+LHQrJYJQoSU677LLLcN1111E1IRIVJUqS83Q6HQ4fPozPPvtM7FBIhqJESXJeZWUlpqam8MYbb4gdCslQlChJzrvmmmtQUFBAl98kIsklSupwTpItOJfOgQMHaJRODptXhfNMRRXOSSocOnQIOp0OTU1NOH36NAoLC7Fp0yYUFxeLHRpJs3A5JuK83oTkkhMnTgAA9uzZA4ZhwDAMzGYzOjs70dDQIG5wRHSSu/QmJNlGRkawdetWAADP8wgEApienkYgEEBjYyNGR0dFjpCIjRIlyXldXV1gGCbsOoZh0NnZmeaISKbJiEtvjuNgt9vBsiw4jkNzczMUCoXYYZEc4fV6Iz7EmZ6eRm9vL1avXo1Vq1Zh5cqVuOSSS+Z1nJGREXR1dcHr9abkHqjU95+OY8x7/3wGUKvVwu8ej4fX6/URt62srExHSCSHtLa28nl5eTyAsK9LLrkk5H1BQQFfXl7Of+c73+F/8Ytf8A6Hg//HP/7BBwKBiMfo6uriZTIZn5eXF/Jz7969STkHqe8/HceId//hcozoT705jkNNTQ1cLpewLD8/HxMTE2G3p6feJNlGRkawatUqBAKBOetkMhmGh4dx6aWXYnh4GO+//z7ef/994fcPPvhAGNFz0UUXYdWqVSGvlStXgmEYXHPNNVH3v2LFipTGn8n7T8cxEtl/Rj71djqdUCqVIcuUSiXcbjfUarVIUZFcUlxcjM7OTjQ2NoJhGPA8L/zs7OwUvkClpaUoLS0N+ez09DS8Xq+QQIOvw4cPx1WRiOd5NDQ04LbbbhPeB9su0X6fuez111+PeOuA53nU1dXhxhtvFLaP9goEAnOWDQ0NRdx/IBDAhg0bsHr16oT2Ofs1OjoaNokFj/HVr34Vy5cvj2u/4ZafPHky4v6D96Hb2toi/n8SPVH6/f6wy30+X9jlwQ7nQfX19VTElyxYQ0MDbr75ZnR2dgr3rxobG2O2YvLy8qBSqaBSqbBhw4aQdSdPnsTw8DC+//3vY3BwMOzneZ6H2+3G8ePHhQdKwe5J8f7+4YcfRk2UHo8HZ8+eFT7DMAxkMlnI+3Cv4Db//ve/o3bEP336ND755JOwn4/3OBzHRf3vvHjxYlx77bVR44y2/wMHDkS8Sg0EAvj1r38tVLoP1+Fc9HuUJpOJ12q1IctYluVtNlvY7ekeJZGaaPdA8/Ly+NbW1pzef6adQ7gcI3r3IIVCMaf16PP56Kk3yRqbNm2K2uJrbGzM6f2n4xgL3b/oiVKr1YZdPvteECFSFbwHKpPJkJeXF/Jz5j3QXN2/FM5B9KfeAKDRaISn3hzHwWAwwOFwhN2WnnoTqRodHU34Hmgu7T8dx4hn/2FzzIIu/JPE4/HwLS0tvM1m41taWviJiYmI24p5j/LFF18U7djJIvVzkHr8PC/9c5B6/Dwf/Rwy8h4lALAsC5PJBL1eD5PJlLH3J7OhtJvUz0Hq8QPSPwepxw8kfg4ZkSiTJdGTT8f/8FTHlOpzyLX45/OZTDsHqcc/n8+k+hwoUaZYrv2RSz3++Xwm085B6vHP5zOpPoeMeJiTiNWrV0OlUoVdNz4+joKCgrj3lertMzEm2l78Y9D24h8j2vYej0fofB4kuURJCCHpllWX3oQQkgqUKAkhJAZKlIQQEoPo1YOkwO12w+l0AgAGBwfR0dGRsX09YzEajdi+fbsk43c6neA4DizLAog8/DUTcRwnlBTkOA56vV44j0zldrvR1NQUUisWkNaMBJHOIeHvdKp6vmcTk8kU8vvMiuxS4nK5eABRRz5lKofDwTc3N/M8f24kF8uyIkeUmJl/QzzPC+eSqWw2m/D3MlsiMxKIKdo5JPqdpkQZg8vl4hUKhfDe4/HwAHiPxyNiVPNjs9l4lmUlmShnxy21//6zv4iZniiDZicZj8cz51xmfj8y0exzmM93mu5RxqBWq9HR0SG8DxYanl2VPdPZ7Xbo9Xqxw5gXjuOE0ntutxt+vz/jL1tnUyqV0Gg0wiV4RUWF2CHNS7QZCaRiPt9pSpRxmJlgenp6oNVqM/aeTDh+v19S8c7mdruhVCqF+2IWiwV2u13ssBJis9kAACqVCjabTbL/aCU6I0GmSvQ7TQ9zEuD3+2G32+fcGM50VqsVzc3NYocxbz6fDxzHCX/Mzc3NyM/Pjzo9QaZxOp0wmUxCGUEAaG9vFzmq5ImUQDNdvN9palEmwGg0wuFwSKp15nQ6UVtbK3YYC8KyLBQKhfDfPfhTKpd7HMdhcHAQWq0Wzc3N8Hg8sFqtMeeJyUTZNiNBvN9palHGyWw2w2g0gmVZ4V9PqfxxWK1W4XeO49DW1oa6ujrJzHIptfuRs7ndbpSVlQnvWZbF9u3bJdkK02q1YVvCUpyRIKHvdKqeNGUTm83GOxwOnud5fmJigm9vbxc5ovmDRJ/Ya7Va4am31LoHBQtTzzT7faZCmO5ks7sHzZ4cMNOEO4dEv9NUFCMGjuPmVCtSKBQRp77MVH6/HxaLBUajEc3NzTAYDJJpUQLn4jcajcK0IcGWgFQ4nU643W6hxaLVajM6fqfTCYfDAbPZjJaWFpSVlQkPQDiOQ3t7O8rKyjA4OJixAxgincN8vtOUKAkhJAZ6mEMIITFQoiSEkBgoURJCSAyUKAkhJAZKlIQQEgMlSkIIiYESJSGExECJkqRcsKN7sOqP0+mExWJJ+1hnp9MJjUYDi8WyoP1oNBrJVS8iC0NjvUlKVVRUoKamBiaTSVjmdrtRUVEBj8eT1li0Wi3q6uoWvB+TyTRnbLPUS9mR6KhFSVLGbDYDwJwSb2q1WtJl32bXLuQ4LqTwCMk+lChJyrS1tQm1F2erqalJczSpM7O1TLITJUqSEhzHRZ2yYWZRCLfbDbPZDLvdDrPZLNy7nHlPMXhfs6amBhzHCZ+ZmYidTidUKhUMBgPsdjssFgsMBkPUcmZOp1M4ttFoBHBu2gyNRiMUTjCbzVCpVDCbzXC73SH3OZ1OJ44ePQqHwwGLxYJnn30WKpVK+IfA7/dDpVIJ+ybSRPcoiag4jhOKpwZpNBr09/dDq9VCq9XC4XAIUynYbDbY7Xa0tLRArVZDpVIJ9we1Wi30ej2WLl0qVLqx2+2oqakJ2f/sYwerW/t8PqHSjFarRXl5ubBvl8slXG7PvM8ZjFGlUgm3E84//3xhnwqFQqjYRKSLWpQkJYKtxUhPtoMtzvb29jnl3liWFe75LV26NKTorUKhCGmlhqu4PfP+oV6vh9PpDNuqbG9vh1KphNPpDJnjObiPjo4OaDQalJaWJvSgprm5WYif4zhJFrUloShRkpRpaWmJOC/MzNqMiUrm02W1Wi20Cpubm4WWa/A4arUaPT09ce8vmJCbm5thsVjgdrslVfeThEeJkqRM8CHH7H6LM1t3dXV1QmsuyO12R53nJ9YUCjPX2+32OU+pg+vDHTv43u/3w+l0wmazgeO4Of0mI8UQ/LzBYKCHPFmE7lGSlApWmDYajVCpVMLcycF7iGq1GiaTCWazGSzLYnBwEDabTZjDO9ia02q1wkOc4OecTic4joPJZILJZBKSocfjES63g/sDMGd/wWMbjUbh8l6r1cJiscBkMgkPisrKytDU1CTMBNnT0wOlUgm9Xg+WZYWkaLFYoNVqAZy7faBWqyU7LS0JRRXOSVYJJuRMeHhit9spUWYJalESkkQGgwE1NTVQKpV0bzKL0D1KkjWCT69tNptoc37X1NTA7/fD5/Nl9ORhJDF06U0IITFQi5IQQmKgREkIITFQoiSEkBj+HyWboyHF822hAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 350x262.5 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with plt.style.context('science'):\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.plot(complexities, rel_bic_scores, color='black', linewidth=2.5//2)\n",
    "    ax.scatter(complexities, rel_bic_scores, marker='o', c='black', s=25)\n",
    "    plt.xlim(min(complexities)-0.3,max(complexities)+0.3)\n",
    "    plt.ylabel(r\"Relative BIC ($\\times 10^{3}$)\")\n",
    "    plt.xlabel(\"Complexity\")\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "af69617f-09fb-43c7-9f48-70b1e4c50fde",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6, 5.24232570772613, 11.774490030146312)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "optimal_slope_index, im1, im2 = check_percent(bic_scores, complexities)\n",
    "optimal_slope_index, im1, im2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "505cc741-c917-4a9f-8f19-0ee41d499d91",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Correctly discover? True\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(6.850210733008831, 2.4196357674983653)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "if target_index == 0: \n",
    "    ground = ground_coeff_u\n",
    "    ground_indices = ground_indices_u\n",
    "elif target_index ==1 : \n",
    "    ground = ground_coeff_v\n",
    "    ground_indices = ground_indices_v\n",
    "correctly_discover = False\n",
    "if br_u.track[optimal_slope_index+1][0] == ground_indices: \n",
    "    correctly_discover = True\n",
    "    print(\"Correctly discover?\", correctly_discover)\n",
    "pred = np.linalg.lstsq(X_pre[:, br_u.track[optimal_slope_index+1][0]], \n",
    "                       y_pre[:, target_index:target_index+1], rcond=None)[0].flatten()\n",
    "errs = percent_coeff(pred, ground)\n",
    "errs.mean(), errs.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "a79c6916-4a16-4b97-aaa2-aa533a155741",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'rel_bic_scores': array([9.54516888e+00, 9.05852146e+00, 8.70006621e+00, 6.40310092e+00,\n",
       "        6.32163782e+00, 4.62194038e+00, 2.05629994e-01, 0.00000000e+00,\n",
       "        3.41193446e-03, 8.97190579e-03, 1.47636167e-02, 2.06342930e-02]),\n",
       " 'complexities': array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12]),\n",
       " 'min_bic_score': -88865.28261169813,\n",
       " 'found_coeff': array([ 0.93556164, -0.95866287, -0.89542675, -0.95640873, -0.895763  ,\n",
       "         0.0942679 ,  0.09359832]),\n",
       " 'percent_coefficient_error_mean': 6.850210733008831,\n",
       " 'percent_coefficient_error_std': 2.4196357674983653}"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "metric_dict = {'rel_bic_scores': rel_bic_scores, \n",
    "               'complexities': complexities, \n",
    "               'min_bic_score': min_bic_score, \n",
    "               'found_coeff': pred, \n",
    "               'percent_coefficient_error_mean': errs.mean(), \n",
    "               'percent_coefficient_error_std': errs.std()}\n",
    "metric_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "58e3d84d-2566-49d5-8ea9-597cb1b8bae5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "./Visualization/reaction_diffustion_2d_big_v_stlsq_bnb_frols_noise_lv15.mat\n"
     ]
    }
   ],
   "source": [
    "if correctly_discover:\n",
    "    target_index2name = ['u', 'v']\n",
    "    file_name = f\"./Visualization/reaction_diffustion_2d_big_{target_index2name[target_index]}_stlsq_bnb_frols_noise_lv{noise_lv}.mat\"\n",
    "    print(file_name)\n",
    "    sio.savemat(file_name, metric_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e659555e-e4f2-4049-bb8e-70ca87c03ff7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "4440d15e-2f78-487c-adc3-380bf3e2c5ce",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0,)\n",
      "(0, 15)\n",
      "(0, 1, 15)\n",
      "(0, 1, 9, 15)\n",
      "(0, 1, 8, 9, 15)\n",
      "(0, 1, 3, 5, 9, 15)\n",
      "(1, 2, 3, 4, 5, 9, 15)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15)\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15)\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15)\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 72)\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 12, 14, 15, 72)\n"
     ]
    }
   ],
   "source": [
    "all_possible_candidates = sorted(set().union(*[br_u.track[mc][0] for mc in br_u.track.keys()]))\n",
    "out = np.hstack([brute_force(X_pre[:, all_possible_candidates], y_pre[:, target_index:target_index+1], \n",
    "                             support_size=i+1) for i in range(len(all_possible_candidates))]).T\n",
    "for i in range(len(out)):\n",
    "    print(tuple(np.array(all_possible_candidates)[list(nonzero(out[i]))]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "47adf1b9-5803-4141-b18d-1daf65deec86",
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_hierarchy = set(tuple(np.array(all_possible_candidates)[list(nonzero(out[i]))]) for i in range(len(out)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "020b75fc-62bf-4e0e-b850-fd0acf484098",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RFE_PLinearRegression\n",
      "(1, 2, 3, 4, 5, 9, 15)\n",
      "(1, 2, 3, 4, 9, 15)\n",
      "(1, 2, 4, 9, 15)\n",
      "(1, 2, 4, 9)\n",
      "(2, 4, 9)\n",
      "(2, 4)\n",
      "(2,)\n",
      "--------------------------------------------------\n",
      "(0, 15)\n",
      "(0,)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15)\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15)\n",
      "(1, 2, 3, 4, 5, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 3, 5, 9, 15)\n",
      "(0, 1, 3, 9, 15)\n",
      "(0, 1, 9, 15)\n",
      "(0, 1, 15)\n",
      "(0, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 72)\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 8, 9, 15)\n",
      "(0, 1, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 12, 14, 15, 72)\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 72)\n",
      "--------------------------------------------------\n",
      "SelectKBest\n",
      "(1, 2, 3, 4, 5, 9, 15)\n",
      "(1, 2, 3, 4, 9, 15)\n",
      "(1, 2, 4, 9, 15)\n",
      "(1, 2, 4, 9)\n",
      "(2, 4, 9)\n",
      "(2, 4)\n",
      "(2,)\n",
      "--------------------------------------------------\n",
      "(0, 15)\n",
      "(0,)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15)\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15)\n",
      "(0, 1, 2, 3, 4, 5, 9, 15)\n",
      "(1, 2, 3, 4, 5, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 3, 5, 9, 15)\n",
      "(0, 1, 3, 9, 15)\n",
      "(0, 1, 9, 15)\n",
      "(0, 1, 15)\n",
      "(0, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 72)\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 8, 9, 15)\n",
      "(0, 1, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 12, 14, 15, 72)\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 72)\n",
      "--------------------------------------------------\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{1: ((0,), -79320.11373424961), 2: ((0, 15), -79806.76114696664), 3: ((0, 1, 15), -80165.2164044002), 4: ((0, 1, 9, 15), -82462.18168864735), 5: ((0, 1, 8, 9, 15), -82543.64479604733), 6: ((0, 1, 3, 5, 9, 15), -84243.34223096781), 7: ((1, 2, 3, 4, 5, 9, 15), -88659.65261778954), 8: ((0, 1, 2, 3, 4, 5, 9, 15), -88865.28261169813), 9: ((0, 1, 2, 3, 4, 5, 9, 14, 15), -88861.87067724303), 10: ((0, 1, 2, 3, 4, 5, 8, 9, 14, 15), -88856.31070591009), 11: ((0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 72), -88850.51899501913), 12: ((0, 1, 2, 3, 4, 5, 8, 9, 12, 14, 15, 72), -88844.6483186886)}"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "br1 = backward_refinement(feature_hierarchy, (X_pre,y_pre[:, target_index:target_index+1]), mode='rfe', ic_type='bic', verbose=True)\n",
    "br2 = backward_refinement(feature_hierarchy, (X_pre,y_pre[:, target_index:target_index+1]), mode='k best', ic_type='bic', verbose=True)\n",
    "br=br1+br2; br"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "4f216e7a-9282-455a-a894-29a5bd699496",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ploting less than max_complexity\n"
     ]
    }
   ],
   "source": [
    "print(\"Ploting less than max_complexity\")\n",
    "for c in range(max_complexity+1, max(list(br.track.keys()))+1):\n",
    "    del_c = br.delete(c)\n",
    "    if del_c is not None: print(f\"Delete {c} complexity\")\n",
    "    else: print(f\"{c} not exist\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "e965e4e2-db07-486b-b53b-96cacf6b1372",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0,) -79320.11373424961\n",
      "(0, 15) -79806.76114696664\n",
      "(0, 1, 15) -80165.2164044002\n",
      "(0, 1, 9, 15) -82462.18168864735\n",
      "(0, 1, 8, 9, 15) -82543.64479604733\n",
      "(0, 1, 3, 5, 9, 15) -84243.34223096781\n",
      "(1, 2, 3, 4, 5, 9, 15) -88659.65261778954\n",
      "(0, 1, 2, 3, 4, 5, 9, 15) -88865.28261169813\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15) -88861.87067724303\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15) -88856.31070591009\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 72) -88850.51899501913\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 12, 14, 15, 72) -88844.6483186886\n"
     ]
    }
   ],
   "source": [
    "# For reproducibility of 2% noise experiment\n",
    "# u\n",
    "# (1,) -90432.11492322256\n",
    "# (1, 14) -90900.56314841735\n",
    "# (0, 1, 14) -91096.59801155927\n",
    "# (0, 1, 8, 14) -92631.07213881976\n",
    "# (0, 1, 8, 14, 15) -92711.06074776711\n",
    "# (0, 1, 2, 4, 8, 14) -94460.72867681735\n",
    "# (0, 2, 3, 4, 5, 8, 14) -99273.17863634705\n",
    "# (0, 1, 2, 3, 4, 5, 8, 14) -99629.82290987358\n",
    "# (0, 1, 2, 3, 4, 5, 8, 14, 67) -99627.38471701351\n",
    "# (0, 1, 2, 3, 4, 5, 8, 14, 67, 72) -99626.97018095759\n",
    "# (0, 1, 2, 3, 4, 5, 8, 9, 14, 67, 72) -99619.72292591527\n",
    "# (0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 67, 72) -99611.13907227767\n",
    "# + stlsq before\n",
    "# (1,) -90432.11492322256\n",
    "# (1, 14) -90900.56314841735\n",
    "# (0, 1, 14) -91096.59801155927\n",
    "# (0, 1, 8, 14) -92631.07213881976\n",
    "# (0, 1, 8, 14, 15) -92711.06074776711\n",
    "# (0, 1, 2, 4, 8, 14) -94460.72867681735\n",
    "# (0, 2, 3, 4, 5, 8, 14) -99273.17863634705\n",
    "# (0, 1, 2, 3, 4, 5, 8, 14) -99629.82290987358\n",
    "# (0, 1, 2, 3, 4, 5, 8, 14, 67) -99627.38471701351\n",
    "# (0, 1, 2, 3, 4, 5, 8, 14, 67, 72) -99626.97018095759\n",
    "# (0, 1, 2, 3, 4, 5, 8, 9, 14, 67, 72) -99619.72292591527\n",
    "# (0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 67, 72) -99611.13907227767\n",
    "\n",
    "# v\n",
    "# (0,) -89704.75044145095\n",
    "# (0, 15) -90094.80085394751\n",
    "# (0, 1, 15) -90433.5949355086\n",
    "# (0, 1, 9, 15) -92311.22556061689\n",
    "# (0, 1, 8, 9, 15) -92462.0683253755\n",
    "# (0, 1, 3, 5, 9, 15) -94229.07163281288\n",
    "# (1, 2, 3, 4, 5, 9, 15) -98884.21234254984\n",
    "# (0, 1, 2, 3, 4, 5, 9, 15) -99223.00785465668\n",
    "# (0, 1, 2, 3, 4, 5, 9, 14, 15) -99226.63012262661\n",
    "# (0, 1, 2, 3, 4, 5, 9, 14, 15, 69) -99230.05524578446\n",
    "# (0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 69) -99224.97738853545\n",
    "# (0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 69, 72) -99216.5636109259\n",
    "# + stlsq before\n",
    "# (0,) -89704.75044145095\n",
    "# (0, 15) -90094.80085394751\n",
    "# (0, 1, 15) -90433.5949355086\n",
    "# (0, 1, 9, 15) -92311.22556061689\n",
    "# (0, 1, 8, 9, 15) -92462.0683253755\n",
    "# (0, 1, 3, 5, 9, 15) -94229.07163281288\n",
    "# (1, 2, 3, 4, 5, 9, 15) -98884.21234254984\n",
    "# (0, 1, 2, 3, 4, 5, 9, 15) -99223.00785465668\n",
    "# (0, 1, 2, 3, 4, 5, 9, 14, 15) -99226.63012262661\n",
    "# (0, 1, 2, 3, 4, 5, 9, 14, 15, 69) -99230.05524578446\n",
    "# (0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 69) -99224.97738853545\n",
    "# (0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 69, 72) -99216.5636109259\n",
    "\n",
    "complexities = np.array(sorted(list(br.track.keys())))\n",
    "bic_scores = []\n",
    "eff_feature_names = []\n",
    "for com in complexities:\n",
    "    indices, bic = br.track[com]\n",
    "    bic_scores.append(bic)\n",
    "    print(indices, bic)\n",
    "bic_scores = np.array(bic_scores)\n",
    "min_bic_score = np.min(bic_scores)\n",
    "div_scale = 1e3\n",
    "rel_bic_scores = (bic_scores-min_bic_score)/div_scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "d5719e52-107d-4840-a35f-d782fdbbf9ad",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUoAAADyCAYAAAArx4ypAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAqS0lEQVR4nO3df1Rb9f0/8OcN1Vl/NCHt5g96JtwUW9eeKQn4a7pOCXWnDbhJAKlnR1aEtN12dNMSul/+2CwmVY+eOVtCoW5zMpK0x0PbaUno1OlxgybUH9tEyCXdyqZbGwKdZ2pL7vePfnM/BPITktzc5PU4Jwdy7829r6vk1fe99/1+vRme53kQQgiJSCZ2AIQQkukoURJCSAyUKAkhJAZKlIQQEkPaEqXb7YZGo5mznOM4mM1m2O12mM1m+P3+dIVECCFxYdLx1Ntut4NlWWg0Gsw+nEajgcvlAnAuaRqNRthstlSHRAghcUtLohQOxjAhiZLjONTU1AiJEgDy8/MxMTGRrpAIISQmUe9ROp1OKJXKkGVKpRJut1ukiAghZK5FYh480v1In88X8TPLly/H2bNnhfcFBQUoKCgAAIyPjwu/xyPV22diTLS9+Meg7cU/xsztx8fHMT4+LqxbtGgRTpw4EfoBPo1mH85kMvFarTZkGcuyvM1mi7iPysrKea0TY/t0HIO2T+726TgGbS/+MRLNI6JeeisUijmtR5/PB4VCMa/91dfXp3T7+Uh1TKk+h1yLfz6fybRzkHr88/lMyr/LCaXhBZp9OI/Hw6vV6pBlCoWCn5iYiLiP+fzrlCxiHjtZpH4OUo+f56V/DlKPn+cl0KKceV+SZdmQdRzHobS0dN4tylRLRws01aR+DlKPH5D+OUg9fiDxc0hL9yCn0wmHwwGz2YyWlhaUlZVBr9cDOJcc29vbUVZWhsHBQWzfvj1qoqyqqkJvb2/IspGREXR1dcHr9aKwsBCbNm1CcXFxKk+JEJKlwuWYtPajTIbZJ7F3717ce++9Qh/N4M/Ozk40NDSIFyghRJLCJUpJj/UeGRnBvffei0AggOnp6ZCfjY2NGB0dFTtEQkgWkHSi7OrqAsMwEdd3dnamMRpCSLaSXKIcHx9HVVUVuru74fV654wdDwoEAujt7cXw8HCaIySESFF3dzeqqqpCOp8HSS5RFhQUoLe3F/X19SgsLIzYomQYBh6PB6tWrcINN9yAXbt2RR3xQwjJbfX19ejt7Q07wkdyiXKmTZs2RWxRMgwDl8sFq9WKZcuW4Xvf+x4uv/xyVFdXo7e3F2fOnElztIQQqZJ0oiwuLkZnZydkMhny8vJCfnZ2dmL16tWoqanBwYMHceLECTz++OPgOA533HEHrrjiCtx3331wu90Rky0hhABZ0D0IAEZHR9HZ2Sn0o2xsbMSKFSsi7uOdd97Br371K/z2t7/FRx99hNWrV+Oee+7B3XffjSuuuCLVp0AIyWBZ2Y9yIc6ePYu+vj78+te/xksvvYQzZ86goqIC99xzD+644w5ceOGFAKhDOyG5JFyOEbXMmtgWLVqE9evXY/369fD7/bDZbPjVr36FjRs34pJLLkFtbS2WLVsGs9kMmUwmdGg3m83UoZ2QHJLTLcpIRkdH8Zvf/AZdXV1z69L9fzKZDMPDw1Ev8Qkh0pN1I3NSZcWKFXjkkUdw9913QyaL/J/oueeeS2NUhBCxUKKM4vjx4xHXBQIBPP3001i7di0ee+wxDA4OYnp6Oo3REULSRXKJcubInFSL1qFdJpOhoqIC+fn5MJlMuO666/D5z38etbW16OjoiJpkCSGZJ9rIHLpHGcXIyAhWrVqFQCAwZ93Me5RnzpzBn//8ZzgcDvT19WFgYACBQABXXXUV1q1bh4qKCnzta1/DkiVLoh6LnqwTIr6wOSZ5NYPTI93Vlffu3cvLZDI+Ly8v5OfevXsjfsbn8/H79u3jDQYDX1RUxAPgFy1axN9yyy38o48+yv/pT3/iz549K2zf1dWV8DEIIakRLsdQizIOiXZon83j8aCvrw99fX04cuQIpqamoFAoUF5ejmuuuQYPPfRQ2NFB9GSdkPSjDucZ4OzZsxgYGEBfXx8cDgfeeuutiEMo8/LysG3bNrS1taU5SkJyF3UPygCLFi3CTTfdhIcffhhvvvkmvvnNb0Z8YMTzPLxeb3oDJITMQYlSZFdddVXEvpoMw6CwsDC9ARFC5qBEKbJopeJ4nkdjY2OaIyKEzEaJUmSzS8UFL8MZhkFnZyc9yCEkA0juYY5Go0FBQQHq6+uzYn7hoJlP1t99913861//AsdxkMvlYodGSE7o7u5Gd3c3xsfH4XK5QtZJLlFK/al3PE6cOIGVK1di69at2Llzp9jhEJJT6Km3RCxfvhxGoxHPPPMMTblLSAagRJmhHnzwQVx22WV48MEHxQ6FkJyXEYV7OY6D0+mEUqkEx3HQ6/VgWVbssER14YUXwmw2o76+Hv39/SgvLxc7JEJyVkYkSrvdjpaWFuG9wWBAe3u7iBFlhrq6Ojz77LO4//77MTQ0hEWLMuJ/FyE5JyMuvXt6esQOISMxDIOnn34a7733Hjo6OsQOh5CclRGJUqlUQqPRCJfgFRUVYoeUMUpLS3HPPffgJz/5Cfx+v9jhEJKTYiZKr9eLnTt3Yt26dSgrKxNet99+O5544omkjEW22WwAAJVKBZvNBr1ev+B9ZpMdO3bgk08+waOPPip2KITkpKj9KFtbW8EwDGpra1FSUjJn/dDQEHp6esAwzIIq3NjtdigUCnAcB4PBgObm5oj3KIMdzoOyreN5JDt27MBDDz2E9957DytXrhQ7HEKyRrCjeVC4DucRC/eazWbe7/fHVejS7/fzra2tCZTG/D8ej4dvaWkJea9QKHiPxxN2+3QX7s0U//vf//jCwkJ+w4YNYodCSFYLl2MiXnpv27Yt7uFzcrl83i1Kt9uNsrIy4T3Lsti+fTvdj5vlggsuwM6dO3Ho0CEcPnxY7HAIySmiP8xRq9UYHBwMWXbq1Cmo1WqRIspc1dXVuOWWW/D9738fZ8+eFTscQnLGnEQ5NjaGLVu24NixY/B6vVi3bh2Ki4tx7NixlATAsiwqKipgNpthsVhgsVhgMBhSciypYxgGzzzzDN5//33s3r1b7HAIyRlzejC3t7dDq9Vi9+7dwlSscrkcHR0duPbaa1MShFarhVarTcm+s01JSQk2bdqEhx56CBs3boRSqRQ7JEKy3pwWZUVFBaqrq7F7925otVqUlJSAZVlKZBnksccew5kzZ/Dwww+LHQohOWFOovT7/di+fTsACOOLn3jiCUxOTqY3MhLRpZdeih//+Md47rnn8Ne//lXscAjJenMSZXV1NVpbW0OWyeVyalFmmPvuuw9XXnklHnjgAbFDISTrhX3qPbtbUFlZGY4ePYr9+/cL81KLZXx8HFVVVSEdRHPR5z73OTzxxBN45ZVX8Pvf/17scAiRvO7ublRVVWF8fHzOuogjcyYnJ9HU1ASn04mlS5cKE2D5fD4wDIPS0lLYbDYsWbIktdHPkgsVzuPF8zzKy8vxz3/+E++++y7OO+88sUMiRPISqnDe0dEBk8kEn8+HkZERjI6OYnR0FD6fD6dOncKuXbsWNGyRLFywutDIyAh++ctfih0OIVkrYqIsKipCUVFRxA+yLIvS0tKUBEXi9+UvfxlNTU145JFHcPLkSbHDISQrRUyUHMfF/PDsETVEHI8++igCgQAeeughsUMhJCtFLJmt1+tRWloKhmFCOjX7fD7h92B5NCKuL3zhC/jpT3+KlpYWbNmyBWvWrBE7JEKySszpasfGxuB0OoUiFQqFAlqtNupleSrRw5zwPvvsM6xZswZXXnkl+vr6wDCM2CERIknhckzMSViKiorQ1NSUsqBIcpx//vl48sknUVVVhQMHDqCqqkrskAjJGguqHrR///5kxUGSQKfToaKiAg888AA+/fRTscMhJGvMO1FOTU2J8jCHOpxHxjAMnnrqKXAch2effVbscAiRlHl1OC8tLcXQ0FDEnfI8D4ZhMD09nbxI40D3KGP7zne+gxdeeAEjIyP4whe+IHY4hEhKwh3Ojx49iunp6bCvQCCAbdu2pTxokrhHHnkEMpkMP/nJT8QOhZCsEDFRlpSUYGxsLOqHaVrZzLRs2TI8/PDD2LNnD95++22xwyFE8qLeo7zzzjujfjhYho1knq1bt6K4uBj3338/YvQAI4TEIPqcOSQ1zjvvPDz11FN49dVX8dJLL4kdDiGSRokyi61fvx5f//rX8eCDD1J3IUIWIK5EuWfPnoSWk8zx1FNP4fjx4/jxj3+M7du3o76+Htu3b8fIyIjYoREiGXElSo1Ggy1btoQs27JlC1UPkoCrr74at956K5544gmYzWZYrVbs3LkTq1atwvPPPy92eIRIQlyJsqSkRCi4AACbN29Gc3NzymZlJMkzMjKC/v5+AEAgEEAgEBC6dzU2NmJ0dDRpx6EWK8lafALcbjevVCr5oaGhRD6WVGq1mq+srORffPFF0WKQktbWVj4vL48HMOclk8l4o9G44GN0dXXxMpmMz8vLC/m5d+/ehZ8AIWny4osv8pWVlbxarZ6zLmb1oJlqa2thMBhgsVjQ09OTotQdHY3MSUx9fT2sVisCgUDY9TKZDEuXLoVCoQh5yeXyOcvCrfvnP/+Jq6++Ouz+ZTIZhoeHsWLFilSfJiFJM6/qQcC5+XOam5vR0dEBuVwOjUaDuro6dHR0pH3OHJKYwsLCiCXXGIZBeXk51q5dC7/fj8nJSfj9fvj9fni9XuF3v9+PM2fORNxHpH9rGYZBZ2cnTRlCJC+uRNnW1gaLxSLMzqhQKGCxWNDW1kZfggy3adMmmM3msOsYhsFzzz0Xs8XH8zw++eSTkMQZfD311FNwuVxhkyXP8/B6vck4DUJEFVeifPzxx+csk8vlSU2STqcTHMeBZVkAoHnEk6S4uBidnZ1obGwUWn/Bn52dnXFdFjMMg8WLF2Px4sW4/PLLQ9a98847GBoaClschWEYFBYWJutUCBFNXIky1ZxOJ2w2G9rb28FxHCoqKuDxeMQOK2s0NDTg5ptvRmdnJ7xeLwoLC9HY2JiUe4fRWqzBJ+uESF3YRDk5OYmJiYmEWgNTU1MAMK97lgaDAS6XC8C52R0dDkfC+yDRrVixIiW3ScK1WIFzSZJlWXzxi19M+jEJSbew/SjlcjkcDkfcFcz37dsHq9U6ryTJcRx8Ph8UCgXcbjf8fr9w+U2koaGhAcPDw9i2bRtqa2vR0tICq9WKv//972hpaRE7PEIWLOKld1NTE4aGhlBbWwuVSoWysjKwLAuFQgG/3w+O4zAwMICxsTEYDAZUV1fPKwC32w2lUgm73Q6tVguLxQKWZaHX6+d9UiT9wrVYP/roI3zve9/DDTfcgLvuukukyAhZuLj6UU5OTsJqtcLj8cDv90OhUEClUiVlNkaLxQKDwYCJiQkhCefn50fscqLRaFBQUCC8r6+vR319/YJiIKnB8zzuvvtu9Pb2YmBgAF/60pfEDomQObq7u0OmlhkfHxduBQYl1OE8FZxOJ2pqajAxMSEsYxgGLpcLarV6zvbU4VxaPv74Y1x//fWYnp7GwMAALrnkErFDIiSqhKaCSBe6H5ndLrroIuzbtw8nTpxAY2MjFREmkpQRibK0tBR+vx8AhL6U4VqTRJpWrlyJvXv3wmaz4ZlnnhE7HEISlhH9KG02G4xGIzQaDVwuF3UPykJ6vR4PPPAAtm3bhtLSUtx8881ih0RI3DIiUSoUCrS3t4sdBkmxtrY2DAwMoLa2Fm63G5dddpnYIRESl4iX3l6vN+QVFOxYTkiizjvvPPT09IDnedx11104e/as2CEREpeIidLhcIBlWZhMJnAcJyw/deoU9u3bF3dndEJmuvzyy9HT04M33ngDP/zhD8UOh5C4RLz0ZlkWLpcLJSUlIcuLiopQVFSEyclJ7N+/P+aUtoTM9tWvfhUmkwkPPvggbrjhBvobIhkvYotycnJyTpKcSS6Xi9LVY3x8HFVVVSEdRIn0/OAHP0B1dTUaGhrwwQcfiB0OIeju7kZVVRXGx8fnrIuYKH0+X8wdz+wkni4FBQXo7e2l0TgSxzAMurq6cPnll6O6uhoff/yx2CGRHFdfX4/e3t6QkX9BERNlPGXOqBQaWYglS5Zg//794DgOmzdvps7oJGNFTJQsy0adt3vPnj1QKpUpCYrkjtWrV6OjowMvvPACdu3aJXY4hIQVtXrQunXrYLPZsHnzZqH4BcdxQp/Hw4cPpydKktU2btyIt956C/fffz80Gg2uv/56sUMiJETUIYx9fX0oLy9HY2MjNBoN1Go19Ho91Go1JUmSVE8++SQ0Gg30ej3+85//iB0OISFijvVuaWmBz+fD6OgoPB4PAoEATShGku7888+HzWbDJ598go0bN4adg4cQsURMlEeOHAl5H+w/OdOxY8dSEhTJTcuXL8fvfvc7HDlyBA8//LDY4RAiiDoy5/Tp05iamor46unpSWesJAeUl5fj5z//OX7+85/j4MGDYodDCIAohXtlMhkYhon4weC0p+m+RApWOKfK5tkrEAjgG9/4Bv74xz/C5XJRzVKSFsFK5wlVON+8eTOMRmPELkCnTp2C2WzG7t27kx9xFFThPDf4/X5oNBrI5XK8+eabWLx4sdghkRyRUIXzmpoaFBUVQS6Xh32xLIuampqUB01yk0KhwL59+/C3v/0N3/3ud8UOh+S4iP0oy8vLY344nm0Ima9rr70Wu3btwre//W2oVCqcPn0aXq8XhYWF2LRpE4qLi8UOkeSIqN2Djhw5gi1btuAPf/gDgHM1KtetW4elS5di69ataQmQ5LaGhgasXbsWP/rRj2A2m2G1WrFz506sWrUKzz//vNjhkRwRMVHu27cPWq0Wg4OD0Ov1OHbsGJqbm6FQKNDa2orR0VFs2bIlnbGSHDQyMoLXX38dwLmHPIFAANPT0wgEAmhsbMTo6KjIEZJcEDFRWq1WTExM4OjRo/B4PNixYweMRiOsViu2bduGvr4+nDp1Kp2xkhzU1dUFmSz8nynDMOjs7ExzRCQXRUyUpaWlkMvlAM7dWK+oqJhzT7KsrCy10ZGc5/V6I1YV4nk+ZJoSQlIlYqLMz88PeR+uL9vsbQhJtsLCwoj9eRmGQWFhYXoDIjkpaj3K06dPC+8nJyfDbpNuVOE8t2zatClqi7KxsTHNEZFsFa3Cedwjc4IjcWa/T/fIHOpwnnuef/55NDY2gmEYBAIB4W+vq6sLDQ0NYodHskxCHc6bm5vh8/mE18TERMh7n8+HpqamlAdNSENDA4aHh7Ft2zbccccdAIAdO3ZQkiRpE7HDucFgEB7mRNuGkHRYsWKFUN7vpptuwsDAgMgRkVwSsUUZbQbGRLYhJNkqKyvR19eHTz/9VOxQSI6IWbg33YxGI/x+v9hhkAym0+nw8ccf47XXXhM7FJIjMipRut1umM1mscMgGW7NmjX44he/iAMHDogdCskRGZUoOY6j2oMkJoZhUFlZiYMHD9IUtyQtMiZR2u126PV6scMgEqHT6eD1evHXv/5V7FBIDkgoUe7cuRN1dXUAgP7+fkxNTSUlCL/fD4VCkZR9kdzwta99DRdddBFdfpO0iDtRtra2QqFQQKvVAjhXi9LpdCYlCKvVKuw3luDInOCLRujkpgsuuAAVFRU0rw5ZsOCInOAr3MiciP0oZysrK0N1dTX6+/uTGqTT6URtbW3c2xcUFNDIHALg3OV3c3MzTp48iWXLlokdDpGo2fNvVVVVzdkm7hbl2NgYAIQMYxwcHFxIfAKr1QqLxQKLxQKO49DW1ga3252UfZPstWHDBgQCAbz88stih0KyXNwtypKSEpSWlmLp0qVwOBxwOp0wmUwLDmD2JbfBYIDBYKCn3ySmyy67DGVlZTh48CC+9a1viR0OyWJxtyjLy8thtVpRUlICnudhsVhw2223JS0Qv98v9KE0mUzUoiRx0el0eOWVV/DZZ5+JHQrJYhGrB822ZcsW7Nq1K9XxxETVg8hMQ0NDUKvV6O/vT+o/3CR3JVQ9aDaHw4E9e/bgyJEjSQ+MkPm69tprccUVV9DTb5JScd+jdLlckMvlmJycxJ49e6BUKqHVarFkyZJUxkdIVAzDQKfT4cCBA3jyyScjVkMnZCHiblEGS67J5XLwPI+WlhaqR0kyQmVlJUZHR/HBBx+IHQrJUnEnyrq6OmzZsgVKpRJutxsOhwM9PT2pjC0smgqCzHbbbbfhggsuoMtvsiDzmgpithUrVqC1tRX33ntv0gNMBD3MIeFUVlbi9OnTePXVV8UOhUjcgh7mmEwm0ZMkIZHodDq88cYbmJiYEDsUkoXiTpTV1dVzltGcyiRTbNiwAdPT0zh8+LDYoZAsFPGp9/79+0Oeau/Zsydkvd/vh8PhoD9MkhGWL1+OkpISHDhwAHfddZfY4ZAsE7FFuWPHDhw9elR4v3v3bkxMTAgvnudx6tSptARJSDx0Oh1efvllnD17VuxQSJaJ2KKcmSQBoKOjY85kYvGWRiMkHXQ6HX72s5/hrbfewi233CJ2OCSLxH2PMj8/X/h9cnIS+/btC1lGiNhKS0tx6aWXUjFfknRxJ8qZRXrlcjmqq6uTVriXkGSQyWTYsGED9ackSRd1COPk5CSsVisYhoHD4Ziz3uVypb3LULDD+exim4QA5y6/u7q64PF4oFKpxA6HSEh3dze6u7vDdjiP2qKUy+XQarU4evQoPB4PRkdHQ14tLS0pCzqSYIVzSpIknIqKCpx//vnUqiQJq6+vR29vLwoKCuasi1kUo6ioCLt370Z/fz/Ky8tTEiAhyXLxxRfj1ltvxcGDB3HfffeJHQ7JEgkV7p3pyJEj2L9/f9IDImShKisr8dprryVtllBC4i6zBpzrhM5xHACA53kcPXoUd955Z0oCI2S+NmzYgO9+97vo6+ujueJJUsSdKFtbW+H3++Hz+cCyLPx+PwwGQypjI2ReCgsLsWbNGhw8eJASJUmKuBOlSqVCU1MTxsbGwDAMCgsLqdo5yViVlZXo6OjA9PQ08vLyxA6HSFzc9yhZlsXx48dRVFQEu92eypgIWTCdToeTJ09iYGBA7FBIFog7Ufr9frAsi6mpKZw8eRK333472tvbUxkbIfN2/fXXY9myZdRNiCRF3IV7Z+vv70dpaakwRUS6aDQaFBQUUIdzEtM999yDoaEhvPPOO2KHQiRgZodzl8sVsi7uFuVs5eXlGBsbW3BwiaIO5yReOp0O7777Lo4fPy52KEQC5tXhfHb9ydkmJiZgtVoxODi48AgJSYF169Zh0aJFOHToELZu3Sp2OETCIrYoZ9efnP0CzvWlJCRTyeVyrF27lqoJkQWL2KI0mUwxhyxSPUqS6XQ6HYxGI/773//i4osvFjscIlERW5ThkuTOnTtRV1cH4NzDnGRVZ3G73TCbzTCbzaipqYHf70/KfgmprKzEZ599RiUByYLE/TCntbUVCoVCaEWWl5cn7Y/P6XSipaUFLS0tKCsro+IbJGlUKhVWrVpF3YTIgsSdKMvKytDU1ASWZZMagNvtRltbm/Ber9fD7XYLY8oJWSidTodDhw4hEAiIHQqRqLgTZbArEMMwwrJkPPFWq9Xo6OgQ3gcvu5VK5YL3TQhw7vL7ww8/nNM3jpB4xT3Wu6SkBKWlpVi6dCkcDgecTidMJlNSgphZuKCnpwdarRYKhSLstsEK50HU8ZzEctNNN0GhUODgwYMoKysTOxySYYIdzYPCVThPaGTO2NiYMGyxrq5uzqyMC+X3+6HRaOByuSImyqqqKvT29ib1uCT7bdy4EcPDw9SqJDGFyzEJjcwpKirC448/jscffxxA7E7piTIajXA4HBGTJCHzVVlZCbfbHba1QEgsURPl1NQU9u/fD6/XO2ddSUlJUrvxmM1mGI1GodYldREiyXT77bcjLy8Phw4dEjsUIkERE+XY2BgKCwuh1+uhUqng9XoxNTWFzZs3o7i4GHl5eUkbvmi326FWq4UkabVaqVVJkkqpVOIrX/kKdRMi8xIxUZpMJthsNgQCAQwODmL37t1C/8bq6mrs3r0bPT09Cw6A4zjU1NSgoqICDMMgPz8fRqNxwfslZLbKyko4nU7873//EzsUIjERn3qrVCohMarVavT398NisST9AQ7LsjRmnKSFTqfDtm3bcOTIEWzYsEHscIiERGxRzuwvCZxLlrOT5LFjx1ISFCGpsHLlSqxYsYIuv0nCIibKU6dO4fTp05iamsLU1BQmJyeF34MvGj9LpIRhGOh0Ohw8eJCuYkhCot6jVCgUyM/PR35+PvR6vfB7cDnVoiRSo9PpcOLECbz99ttih0IkJGKibG5uhs/nE14TExMhv586dQpFRUXpjBXA/43MmdmTnpB43XLLLViyZAnVqCRzdHd3o6qqKrGROUNDQzEf3MSzTbLRyByyULW1tTh+/Dj+/Oc/ix0KyUAJjcyJJwGmO0kSkgw6nQ4DAwP46KOPxA6FSMS8JxcjRKrWr18PhmFolA6JGyVKknOWLVuGG2+8kboJkbhRoiQ5qbKyEn19ffj000/FDoVIACVKkpN0Oh0+/vhjvPrqq2KHQiSAEiXJSatXr8aVV15Jl98kLpQoSU5iGAaVlZU0SofERXKJkjqck2TR6XTwer34y1/+InYoJANE63AuuURZUFCA3t5emieHLNjatWtx0UUX0eU3AXBu/q3e3l4UFBTMWSe5RElIslxwwQVYt24dJUoSEyVKktN0Oh3eeustnDx5UuxQSAajREly2vr16xEIBPDyyy+LHQrJYJQoSU677LLLcN1111E1IRIVJUqS83Q6HQ4fPozPPvtM7FBIhqJESXJeZWUlpqam8MYbb4gdCslQlChJzrvmmmtQUFBAl98kIsklSupwTpItOJfOgQMHaJRODptXhfNMRRXOSSocOnQIOp0OTU1NOH36NAoLC7Fp0yYUFxeLHRpJs3A5JuK83oTkkhMnTgAA9uzZA4ZhwDAMzGYzOjs70dDQIG5wRHSSu/QmJNlGRkawdetWAADP8wgEApienkYgEEBjYyNGR0dFjpCIjRIlyXldXV1gGCbsOoZh0NnZmeaISKbJiEtvjuNgt9vBsiw4jkNzczMUCoXYYZEc4fV6Iz7EmZ6eRm9vL1avXo1Vq1Zh5cqVuOSSS+Z1nJGREXR1dcHr9abkHqjU95+OY8x7/3wGUKvVwu8ej4fX6/URt62srExHSCSHtLa28nl5eTyAsK9LLrkk5H1BQQFfXl7Of+c73+F/8Ytf8A6Hg//HP/7BBwKBiMfo6uriZTIZn5eXF/Jz7969STkHqe8/HceId//hcozoT705jkNNTQ1cLpewLD8/HxMTE2G3p6feJNlGRkawatUqBAKBOetkMhmGh4dx6aWXYnh4GO+//z7ef/994fcPPvhAGNFz0UUXYdWqVSGvlStXgmEYXHPNNVH3v2LFipTGn8n7T8cxEtl/Rj71djqdUCqVIcuUSiXcbjfUarVIUZFcUlxcjM7OTjQ2NoJhGPA8L/zs7OwUvkClpaUoLS0N+ez09DS8Xq+QQIOvw4cPx1WRiOd5NDQ04LbbbhPeB9su0X6fuez111+PeOuA53nU1dXhxhtvFLaP9goEAnOWDQ0NRdx/IBDAhg0bsHr16oT2Ofs1OjoaNokFj/HVr34Vy5cvj2u/4ZafPHky4v6D96Hb2toi/n8SPVH6/f6wy30+X9jlwQ7nQfX19VTElyxYQ0MDbr75ZnR2dgr3rxobG2O2YvLy8qBSqaBSqbBhw4aQdSdPnsTw8DC+//3vY3BwMOzneZ6H2+3G8ePHhQdKwe5J8f7+4YcfRk2UHo8HZ8+eFT7DMAxkMlnI+3Cv4Db//ve/o3bEP336ND755JOwn4/3OBzHRf3vvHjxYlx77bVR44y2/wMHDkS8Sg0EAvj1r38tVLoP1+Fc9HuUJpOJ12q1IctYluVtNlvY7ekeJZGaaPdA8/Ly+NbW1pzef6adQ7gcI3r3IIVCMaf16PP56Kk3yRqbNm2K2uJrbGzM6f2n4xgL3b/oiVKr1YZdPvteECFSFbwHKpPJkJeXF/Jz5j3QXN2/FM5B9KfeAKDRaISn3hzHwWAwwOFwhN2WnnoTqRodHU34Hmgu7T8dx4hn/2FzzIIu/JPE4/HwLS0tvM1m41taWviJiYmI24p5j/LFF18U7djJIvVzkHr8PC/9c5B6/Dwf/Rwy8h4lALAsC5PJBL1eD5PJlLH3J7OhtJvUz0Hq8QPSPwepxw8kfg4ZkSiTJdGTT8f/8FTHlOpzyLX45/OZTDsHqcc/n8+k+hwoUaZYrv2RSz3++Xwm085B6vHP5zOpPoeMeJiTiNWrV0OlUoVdNz4+joKCgrj3lertMzEm2l78Y9D24h8j2vYej0fofB4kuURJCCHpllWX3oQQkgqUKAkhJAZKlIQQEoPo1YOkwO12w+l0AgAGBwfR0dGRsX09YzEajdi+fbsk43c6neA4DizLAog8/DUTcRwnlBTkOA56vV44j0zldrvR1NQUUisWkNaMBJHOIeHvdKp6vmcTk8kU8vvMiuxS4nK5eABRRz5lKofDwTc3N/M8f24kF8uyIkeUmJl/QzzPC+eSqWw2m/D3MlsiMxKIKdo5JPqdpkQZg8vl4hUKhfDe4/HwAHiPxyNiVPNjs9l4lmUlmShnxy21//6zv4iZniiDZicZj8cz51xmfj8y0exzmM93mu5RxqBWq9HR0SG8DxYanl2VPdPZ7Xbo9Xqxw5gXjuOE0ntutxt+vz/jL1tnUyqV0Gg0wiV4RUWF2CHNS7QZCaRiPt9pSpRxmJlgenp6oNVqM/aeTDh+v19S8c7mdruhVCqF+2IWiwV2u13ssBJis9kAACqVCjabTbL/aCU6I0GmSvQ7TQ9zEuD3+2G32+fcGM50VqsVzc3NYocxbz6fDxzHCX/Mzc3NyM/Pjzo9QaZxOp0wmUxCGUEAaG9vFzmq5ImUQDNdvN9palEmwGg0wuFwSKp15nQ6UVtbK3YYC8KyLBQKhfDfPfhTKpd7HMdhcHAQWq0Wzc3N8Hg8sFqtMeeJyUTZNiNBvN9palHGyWw2w2g0gmVZ4V9PqfxxWK1W4XeO49DW1oa6ujrJzHIptfuRs7ndbpSVlQnvWZbF9u3bJdkK02q1YVvCUpyRIKHvdKqeNGUTm83GOxwOnud5fmJigm9vbxc5ovmDRJ/Ya7Va4am31LoHBQtTzzT7faZCmO5ks7sHzZ4cMNOEO4dEv9NUFCMGjuPmVCtSKBQRp77MVH6/HxaLBUajEc3NzTAYDJJpUQLn4jcajcK0IcGWgFQ4nU643W6hxaLVajM6fqfTCYfDAbPZjJaWFpSVlQkPQDiOQ3t7O8rKyjA4OJixAxgincN8vtOUKAkhJAZ6mEMIITFQoiSEkBgoURJCSAyUKAkhJAZKlIQQEgMlSkIIiYESJSGExECJkqRcsKN7sOqP0+mExWJJ+1hnp9MJjUYDi8WyoP1oNBrJVS8iC0NjvUlKVVRUoKamBiaTSVjmdrtRUVEBj8eT1li0Wi3q6uoWvB+TyTRnbLPUS9mR6KhFSVLGbDYDwJwSb2q1WtJl32bXLuQ4LqTwCMk+lChJyrS1tQm1F2erqalJczSpM7O1TLITJUqSEhzHRZ2yYWZRCLfbDbPZDLvdDrPZLNy7nHlPMXhfs6amBhzHCZ+ZmYidTidUKhUMBgPsdjssFgsMBkPUcmZOp1M4ttFoBHBu2gyNRiMUTjCbzVCpVDCbzXC73SH3OZ1OJ44ePQqHwwGLxYJnn30WKpVK+IfA7/dDpVIJ+ybSRPcoiag4jhOKpwZpNBr09/dDq9VCq9XC4XAIUynYbDbY7Xa0tLRArVZDpVIJ9we1Wi30ej2WLl0qVLqx2+2oqakJ2f/sYwerW/t8PqHSjFarRXl5ubBvl8slXG7PvM8ZjFGlUgm3E84//3xhnwqFQqjYRKSLWpQkJYKtxUhPtoMtzvb29jnl3liWFe75LV26NKTorUKhCGmlhqu4PfP+oV6vh9PpDNuqbG9vh1KphNPpDJnjObiPjo4OaDQalJaWJvSgprm5WYif4zhJFrUloShRkpRpaWmJOC/MzNqMiUrm02W1Wi20Cpubm4WWa/A4arUaPT09ce8vmJCbm5thsVjgdrslVfeThEeJkqRM8CHH7H6LM1t3dXV1QmsuyO12R53nJ9YUCjPX2+32OU+pg+vDHTv43u/3w+l0wmazgeO4Of0mI8UQ/LzBYKCHPFmE7lGSlApWmDYajVCpVMLcycF7iGq1GiaTCWazGSzLYnBwEDabTZjDO9ia02q1wkOc4OecTic4joPJZILJZBKSocfjES63g/sDMGd/wWMbjUbh8l6r1cJiscBkMgkPisrKytDU1CTMBNnT0wOlUgm9Xg+WZYWkaLFYoNVqAZy7faBWqyU7LS0JRRXOSVYJJuRMeHhit9spUWYJalESkkQGgwE1NTVQKpV0bzKL0D1KkjWCT69tNptoc37X1NTA7/fD5/Nl9ORhJDF06U0IITFQi5IQQmKgREkIITFQoiSEkBj+HyWboyHF822hAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 350x262.5 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with plt.style.context('science'):\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.plot(complexities, rel_bic_scores, color='black', linewidth=2.5//2)\n",
    "    ax.scatter(complexities, rel_bic_scores, marker='o', c='black', s=25)\n",
    "    plt.xlim(min(complexities)-0.3,max(complexities)+0.3)\n",
    "    plt.ylabel(r\"Relative BIC ($\\times 10^{3}$)\")\n",
    "    plt.xlabel(\"Complexity\")\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "bb8b0ec2-6650-43b9-b37b-954041d4c9fb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6, 5.24232570772613, 11.774490030146312)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "check_percent(bic_scores, complexities)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "9d2fbfd9-8a10-4ac7-89ca-b5312a9c24ea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'rel_bic_scores': array([9.54516888e+00, 9.05852146e+00, 8.70006621e+00, 6.40310092e+00,\n",
       "        6.32163782e+00, 4.62194038e+00, 2.05629994e-01, 0.00000000e+00,\n",
       "        3.41193446e-03, 8.97190579e-03, 1.47636167e-02, 2.06342930e-02]),\n",
       " 'complexities': array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12]),\n",
       " 'min_bic_score': -88865.28261169813,\n",
       " 'found_coeff': array([ 0.93556164, -0.95866287, -0.89542675, -0.95640873, -0.895763  ,\n",
       "         0.0942679 ,  0.09359832]),\n",
       " 'percent_coefficient_error_mean': 6.850210733008831,\n",
       " 'percent_coefficient_error_std': 2.4196357674983653}"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "metric_dict = {'rel_bic_scores': rel_bic_scores, \n",
    "               'complexities': complexities, \n",
    "               'min_bic_score': min_bic_score, \n",
    "               'found_coeff': pred, \n",
    "               'percent_coefficient_error_mean': errs.mean(), \n",
    "               'percent_coefficient_error_std': errs.std()}\n",
    "metric_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "51a0a7f1-5302-46c7-9974-0466b1de2282",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "./Visualization/reaction_diffustion_2d_big_v_stlsq_bnb_frols_bruteforce_noise_lv15.mat\n"
     ]
    }
   ],
   "source": [
    "if correctly_discover:\n",
    "    target_index2name = ['u', 'v']\n",
    "    file_name = f\"./Visualization/reaction_diffustion_2d_big_{target_index2name[target_index]}_stlsq_bnb_frols_bruteforce_noise_lv{noise_lv}.mat\"\n",
    "    print(file_name)\n",
    "    sio.savemat(file_name, metric_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "1e653741-9957-441d-8960-9cd87baa084f",
   "metadata": {},
   "outputs": [],
   "source": [
    "d_bic_scores, d_complexities, _ = get_decreasing_vals(bic_scores, complexities)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "e77fa604-381a-4851-96a9-50d06ba8ead7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7 improves 1\n",
      "Percent improve: 0.11774490030146312\n",
      "8 does not improve 7\n",
      "Percent improve: 0.0023193187412436115\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "7"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# also tested with percent_ok = 0.09\n",
    "find_transition_V2(d_bic_scores, d_complexities, 0.095)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27bc0653-e9aa-4a09-9f11-1537e575688e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:py3.7] *",
   "language": "python",
   "name": "conda-env-py3.7-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.7.12"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}