Examples/Reaction-Diffusion-2D-big-nPIML.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"
   ]
  },
  {
   "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": [],
   "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==1])\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])"
   ]
  },
  {
   "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.179 v + -0.815 uuu + 0.769 vvv + -0.813 uvv + 0.767 uuv + 0.083 u_22 + 0.084 u_11\n",
      "(v)' = -0.151 u + 0.855 v + -0.799 uuu + -0.818 vvv + -0.797 uvv + -0.818 uuv + 0.086 v_22 + 0.086 v_11 + 0.009 vvvv_22 + -0.010 uvvv_22 + -0.011 uuvv_22 + 0.036 uu_1 + -0.085 uuuu_1 + 0.040 vvvv_1 + -0.090 uvvv_1 + -0.010 vvvv_11 + 0.011 uvvv_11 + 0.012 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=False)\n",
    "model.print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e1038a14-5752-4f4d-ba5a-541abed6c517",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "u_t = (0.067813 +0.000000i)u\n",
      "    + (1.103205 +0.000000i)v\n",
      "    + (-0.073863 +0.000000i)uuu\n",
      "    + (-0.205475 +0.000000i)vvv\n",
      "    + (-0.072619 +0.000000i)uvv\n",
      "    + (-0.210551 +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": [
      "u_t = (-1.100886 +0.000000i)u\n",
      "    + (0.037555 +0.000000i)v\n",
      "    + (0.203575 +0.000000i)uuu\n",
      "    + (-0.041613 +0.000000i)vvv\n",
      "    + (0.206424 +0.000000i)uvv\n",
      "    + (-0.039739 +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)}\n",
      "stridge_v ...\n",
      "{(0, 1, 2, 3, 4, 5)}\n",
      "{(0, 1, 2, 3, 4, 5, 9, 15, 35, 37, 39, 40, 44, 47, 49, 71, 73, 75), (0, 1, 2, 3, 4, 5)}\n"
     ]
    }
   ],
   "source": [
    "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),\n",
       " (0, 2, 4, 8, 9, 14, 15)}"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "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, 14, 15),\n",
       " (0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 65, 70),\n",
       " (0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 70),\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, 1, 2, 3, 4, 5, 8, 9, 14, 15, 65, 70)\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 70)\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15, 70)\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",
      "(1, 2, 4, 15)\n",
      "(2, 4, 15)\n",
      "(2, 4)\n",
      "(2,)\n",
      "--------------------------------------------------\n",
      "(0, 2, 4, 8)\n",
      "(0, 2, 4)\n",
      "(0, 4)\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",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5)\n",
      "(0, 1, 2, 3, 4)\n",
      "(0, 2, 3, 4)\n",
      "(0, 2, 4)\n",
      "--------------------------------------------------\n",
      "(0, 1, 9, 15)\n",
      "(0, 1, 15)\n",
      "(0, 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",
      "--------------------------------------------------\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, 8, 9, 14, 15)\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 70)\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, 39, 40, 44, 47, 49)\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)\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, 1, 2, 3, 4, 5, 8, 9, 14, 15, 65, 70)\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 70)\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15, 70)\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",
      "(1, 2, 4, 15)\n",
      "(2, 4, 15)\n",
      "(2, 4)\n",
      "(2,)\n",
      "--------------------------------------------------\n",
      "(0, 2, 4, 8)\n",
      "(0, 2, 4)\n",
      "(0, 4)\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",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5)\n",
      "(0, 1, 2, 3, 4)\n",
      "(0, 2, 3, 4)\n",
      "(0, 2, 4)\n",
      "--------------------------------------------------\n",
      "(0, 1, 9, 15)\n",
      "(0, 1, 15)\n",
      "(0, 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",
      "--------------------------------------------------\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, 8, 9, 14, 15)\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 70)\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, 39, 40, 44, 47, 49)\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)\n",
      "--------------------------------------------------\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{1: ((0,), -79286.70104043813), 2: ((0, 15), -79782.41500301847), 3: ((0, 1, 15), -80148.37059246466), 4: ((0, 1, 9, 15), -82433.05404289016), 5: ((0, 1, 8, 9, 15), -82510.35281231243), 6: ((0, 1, 8, 9, 14, 15), -82591.15240150668), 7: ((1, 2, 3, 4, 5, 9, 15), -88655.88617596863), 8: ((0, 1, 2, 3, 4, 5, 9, 15), -88860.20954906201), 9: ((0, 1, 2, 3, 4, 5, 9, 14, 15), -88857.16238963649), 10: ((0, 1, 2, 3, 4, 5, 9, 14, 15, 70), -88852.8420606986), 11: ((0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 70), -88847.59336399431), 12: ((0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 65, 70), -88841.99953940697), 13: ((0, 1, 2, 3, 4, 5, 9, 15, 39, 40, 44, 47, 49), -88828.9240292045), 14: ((0, 1, 2, 3, 4, 5, 9, 15, 39, 40, 44, 47, 49, 75), -88822.6988259611), 15: ((0, 1, 2, 3, 4, 5, 9, 15, 39, 40, 44, 47, 49, 73, 75), -88813.64405214202), 16: ((0, 1, 2, 3, 4, 5, 9, 15, 39, 40, 44, 47, 49, 71, 73, 75), -88804.7869661112), 17: ((0, 1, 2, 3, 4, 5, 9, 15, 37, 39, 40, 44, 47, 49, 71, 73, 75), -88796.20936907997), 18: ((0, 1, 2, 3, 4, 5, 9, 15, 35, 37, 39, 40, 44, 47, 49, 71, 73, 75), -88796.10294784103)}"
      ]
     },
     "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,) -79286.70104043813\n",
      "(0, 15) -79782.41500301847\n",
      "(0, 1, 15) -80148.37059246466\n",
      "(0, 1, 9, 15) -82433.05404289016\n",
      "(0, 1, 8, 9, 15) -82510.35281231243\n",
      "(0, 1, 8, 9, 14, 15) -82591.15240150668\n",
      "(1, 2, 3, 4, 5, 9, 15) -88655.88617596863\n",
      "(0, 1, 2, 3, 4, 5, 9, 15) -88860.20954906201\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15) -88857.16238963649\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15, 70) -88852.8420606986\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 70) -88847.59336399431\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 65, 70) -88841.99953940697\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/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAo/klEQVR4nO3df3RT9f0/8OdNBfRMmzSosNUfaQLyS6dNGhIPbnrWFJxn9EybUnrOnEghAXVzHqUpx03PNqUm6BS3HdqUFvfj0LUJPQfcDmJS5pzOQmlgKgq0uQ1s3dRJmxR04CT3+wff3E/T5qZJm+TmJq/HOTltkpt7X1ebF+973+/3681wHMeBEEKIIJnYARBCSLajREkIIZOgREkIIZOgREkIIZPIWKL0+XzQ6XQTXmdZFg6HA263Gw6HA8FgMFMhEUJIQphM9Hq73W6o1WrodDqMP5xOp0NfXx+AS0nTZrPB5XKlOyRCCElYRhIlfzCGiUqULMuiurqaT5QAUFRUhJGRkUyFRAghk7pMzIN7vV4olcqo15RKJXw+H7RabczPLFmyBBqNJuZ7Q0NDKC4uTvj46d4+G2Oi7cU/Bm0v/jHibe/3+3Hs2LHoF7kMGn84u93OmUymqNfUajXn8XgE96HVarmVK1fyj127dvHvrVy5Mql40r19Jo5B26d2+0wcg7YX/xhjt9+1a1dUTtFqtRO2F7VFKSReh05xcTH27t0b873a2tqkjpPs9lOR7pjSfQ75Fv9UPpNt5yD1+KfymemcQ21tbdTnKysrJ26UVBqepvGHa25unpC9FQpF3BblVP51ShUxj50qUj8HqcfPcdI/B6nHz3HxzyHWe6KOozSZTDFfLysry3AkiclECzTdpH4OUo8fkP45SD1+IPlzyHiv98jICBQKBf/a+OFBVqsVHo9HcB+VlZUTLr37+/vR1taGQCAAlUqFtWvXYv78+Wk5B0JIbouVYzJyj9Lr9fLJr7GxEXq9HmazGQDgcrlgs9mg1+vR29ub9BjKnTt3Yt26dfzQI4Zh4HA40NraijVr1qT6VAgheSijLcpUGJvt+/v7sXDhQoTD4QnbyWQynDhxAvPmzct0iIQQCYvVopT0XO+2tjYwDBPzPYZh0NramuGICCG5SNKJMhAITJgSGREOh3Hy5MkMR0QIyUWSS5RDQ0OorKxEe3s7VCqVYIuS4zjs2bMHNTU12LdvH7788ssMR0oIkZL29nZUVlZiaGhowns5e4+SYRjYbDa8+uqrOHbsGL72ta/h/vvvx5o1a7Bw4cJMh00IkYicu0c5f/58tLa2QiaToaCgIOpnW1sbGhsb8d5776G3txf33nsvnE4nFi1ahNtvvx3Nzc1U0o0QkhBJtygjBgYG0Nrayo+jrKuri9nbfeHCBezduxevvPIKXnvtNcycORP33nsv1qxZg/LychQUFGTqNAghWSpWjsmJRDkV//rXv/D73/8eO3fuxPHjx3Hdddfh+9//Ph544AHcdNNNUdvSgHZC8kfMHJPiKZRpl+p5puFwmOvp6eE2bNjAyeVyDgC3bNkyrqWlhQuFQlxbWxsnk8m4goKCqJ87d+5MaRyEkOyQdXO9swHDMDAYDNi+fTv+/e9/o729HVdeeSUsFgvmzJmDtWvXIhwO4+LFi1E/6+rqMDAwIHb4hJAMyPtEOdYVV1yB1atX47XXXsPp06dhMBgEt6UB7YTkD0qUAq677jp89atfFRynefHiRRw4cCDmmCtCSG6hRBmHSqWCTCb8n+jQoUO47rrrcMstt+CJJ56Ax+PB+fPnMxghISQTJJcox87MSbe1a9cKTpGUyWQ4ePAg2tvbUVZWhl27dmH58uVQKpX49re/jZdeegkffvih4OcJIdklZ2fmZMIrr7yCurq6qDJuHMdNKOPGcRzef/997N+/H/v378ebb76JL774AjfccANWrFiBFStWoLy8PKoW51g0BImQ7EDjKKco0QHtY3322Wf4y1/+wifOEydOoKCgAAaDgU+cZWVlKCgoiFlTM1YyJoSkHyVKEZ06dYpPmt3d3QiFQlAqlTAajdi3b1/MS3SqqUlI5uXcXG8pufHGG2GxWLB79258+umneOutt/Dwww/jyJEjgvcxaQgSIdkhK5erzXWXXXYZli1bhmXLlqG/vx8dHR0xkyXHcQgEApkPkBAShVqUIos3BIlhGKhUqswGRAiZgBKlyOINQeI4DnV1dRmOiBAyHiVKkY2vqTl2JlBrayt15BCSBSTX663T6VBcXIza2tqcWIg9YuwQpKGhIfT09ODYsWM0lpKQDGlvb0d7ezuGhobQ19cX9Z7kEqVUhwcl47///S8WLVqEr3/96zl/roRkGxoeJBFXXHEFtm7dildffRWvv/662OEQkvcoUWYps9mMb3zjG3jsscdoBUlCRJYViZJlWTidTrjdbjgcDrAsK3ZIomMYBtu2bcOHH36IpqYmscMhJK9lxYBzt9uN+vp6/rnVakVzc7OIEWWH0tJSrFu3Dk899RRqa2sxe/ZssUMiJC9lRYuyo6ND7BCy1jPPPIOLFy/i6aefFjsUQvJWViRKpVIJnU4HlmXh9XpRUVEhdkhZ49prr8VTTz2FpqYmvP/++2KHQ0hemjRRBgIBbN26FcuXL4der+cfK1aswPPPP5+SucgulwsAoNFo4HK5YDabp73PXPKDH/wAarUajz32GBUCJkQEccdRNjQ0gGEYrFq1CqWlpRPeP3LkCDo6OsAwDBobG6cchNvthkKhAMuysFqtsFgsgvcoIwPOI3Jt4LmQP/7xj1i5ciX27NmDyspKscMhJGdEBppHxBpwLriut8Ph4ILBYELr4AaDQa6hoSGhbcfz+/1cfX191HOFQsH5/f6Y26d6XW+pCIfD3IoVKziNRsOdP39e7HAIyVlJreu9adMmyOXyhDKyXC6fcovS5/NBr9fzz9VqNTZv3oxgMDil/eUqhmHwi1/8AoFAANu2bRM7HELyiuidOVqtFr29vVGvnTlzBlqtVqSIstfixYvx8MMP45lnnsFHH30kdjiE5I0JiXJwcBAbN27E0aNHEQgEsHz5csyfPx9Hjx5NSwBqtRoVFRVwOBxwOp1wOp2wWq1pOVYuePrppzFjxgw8+eSTYodCSN6YMOC8ubkZJpMJTU1NKCoqgt1uh1wuR0tLC2677ba0BGEymWAymdKy71yjVCrx85//HI888ggeeugh6HQ6sUMiJOdNaFFWVFSgqqoKTU1NMJlMKC0thVqtpkSWRSwWC5YsWYJHH32UhgsRkgETEmUwGMTmzZsBAOXl5QCA559/HqFQKLOREUGXXXYZXnrpJbz99ts0q4mQDJiQKKuqqtDQ0BD1mlwupxZllikvL8d3v/td1NfX4/PPPxc7HEJyWsxe7/HDgvR6PQ4fPoyuri4cOHAAo6OjGQkulqGhIVRWVkYNEM1Xzz//PD7++GNs3bpV7FAIkbz29nZUVlZiaGhownuCM3NCoRDWr18Pr9eL2bNn8/fChoeHwTAMysrK4HK5UFhYmN7ox8mHCufJaGhowMsvv4wTJ07g+uuvFzscQiQvqQrnLS0tsNvtGB4eRn9/PwYGBjAwMIDh4WGcOXMG27dvn9a0RZIaTz75JAoLC2Gz2cQOhZCcJZgoS0pKUFJSIvhBtVqNsrKytARFEnfVVVehsbER7e3tePvtt8UOh5CcJJgoE6kyPn5GDRHHAw88AJ1Oh0cffRThcFjscAjJOYIVzs1mM8rKysAwDJRKJf/68PAw/3ukPBoRl0wmw7Zt23DHHXfgN7/5DR588EGxQyIkpwgmypKSEhw+fBiDg4Pwer18kQqFQgGTyRT3spxk3rJly1BbW4vNmzejqqoq451shOSySdfMKSkpwfr16zMRC5kmu92OBQsWYMuWLXjuuefEDoeQnDGt6kFdXV2pioOkwPXXXw+bzYYXX3wRfr9f7HAIyRlTTpSjo6OidObQgPP4Nm3ahDlz5uCJJ54QOxRCJGVKA87Lyspw5MgRwZ1yHAeGYXDx4sXURZoAGnA+uT/84Q+ora2F1+vl5+sTQhKT9IDzw4cP4+LFizEf4XAYmzZtSnvQJHk1NTVYtmwZfvSjH+HLL78UOxxCJE8wUZaWlmJwcDDuh2lZ2ezEMAy2bduGY8eOwel0ih0OIZIX9x7lfffdF/fDdFmXvXQ6HR588EH85Cc/iRr7SghJnuhr5pD0efbZZ/G///0PP/3pT8UOhRBJo0SZw+bOnYsf//jH+PWvf40PPvhA7HAIkayEEuWOHTuSep1kj0cffRQqlQpWqxUNDQ387J3+/n6xQyNEMhJKlDqdDhs3box6bePGjVQ9SAJmzZqFe+65B2+99Ra2bt2Kzs5ObN26FQsXLsQrr7ySsuP09/dj8+bNlIhJbuISxLIst2HDBo7jOM5qtXI+ny/Rj6bUypUrRTmuVJ08eZJjGIYDMOEhk8m4/v7+aR+jra2Nk8lkXEFBQdTPnTt3Tv8ECMmwWDlGcMB5LEeOHIHJZEJ3d3falq6djE6nQ3FxMWpra1FbWytKDFKyefNmbN26NebEAIZhcNddd2H16tWQy+UxH1/5ylfAMIzg/vv7+7Fw4cKY5d1kMhlOnDiBefPmTfs8+vv70dbWhkAgAJVKhbVr12L+/PnT3i8hEe3t7Whvb8fQ0BD6+vqi30wm01ZXV3Ner5dbtWpVinJ38qhFmZzVq1dzMpksZosSADdjxgzBFicArqCggCsqKuJUKhV36623ct/85je5lStXct/73ve4hx9+mLv99tvjtlgff/xxLhwOT+scqMVKMilWjpm0ehBwaf0ci8WClpYWyOVy6HQ61NTUoKWlhcp5ZTmVSiXYIiwoKMDjjz+OZ599FufOnUMoFEIoFEIwGOR/F3oEAgGEQiH4/X7BtcXD4TBeeOEF/PKXv4RcLodCoYBCoUjo98jzM2fOoK6uLuYx6urqcMcdd6SkxUpIPAklysbGRjidTn51RoVCAafTicbGRlo3J8utXbsWDocj5nscx6Gurg4ymQyFhYUoLCxMeoGyeJf2MpkM99xzD+6++24Eg0E+AUd+P3XqVNTzL774IqljMwyD1tZW+hskaZdQooxV21Aul6f0D9Tr9YJlWajVagCgdcRTZP78+WhtbUVdXR0YhuGLmXAch9bW1mm3xuIlYgB48cUXEz7G+fPnJyTUn/3sZ3jnnXditig5jkMgEJhq6IQkLKFEmW5erxculwvNzc1gWRYVFRVUTzGF1qxZgzvuuAOtra18Z0hdXV1KLllTmYgvv/xyzJ07F3PnzuVfe+ONN3Dw4EHBziiVSjXtcyBkMjF7vUOhEEZGRpL6IxwdHQWAKd2z1Gg06Ovrg0KhAIColuV4VGYtOw0MDKQlEWeqV52QiFg5JmaLUi6Xo7OzE7Nnz560MAYA7N69GyMjI1i3bl3SQbEsi+HhYSgUCvh8PqjVasEkSbLXvHnz0nKvcHyLNdKylMlkKbl1QEgiBC+9169fjyNHjmDVqlXQaDTQ6/VQq9VQKBQIBoNgWRaHDh3C4OAgrFYrqqqqphSAz+eDUqmE2+2GyWSC0+mEWq2G2Wye8kmR3DL21sGbb76Jv/3tb3j77bdhNBrFDo3kiYQGnIdCIXR2dsLv9yMYDEKhUECj0aRkNUan0wmr1YqRkRE+CRcVFQkOOYkMOI+ggef55fTp07jxxhuxZ88eVFZWih0OyQGRgeYRsQacJzUzJx28Xi+qq6sxMjLCv8YwDPr6+qDVaidsT/co8xvHcSguLsaaNWuwZcsWscMhOSippSAyhe5HkmQwDAODwYCDBw+KHQrJI1mRKMvKyhAMBgH8X493rNYkIQBgNBpx6NChjC9sR/JXVoyjdLlcsNls0Ol06Ovrg8fjETskksUMBgPOnTuHDz74ALfccovY4ZA8kBWJUqFQoLm5WewwiESUlZVBJpPh4MGDlChJRgheegcCgahHRGRgOSFiufLKK3HzzTejp6dH7FBInhBMlB6PB2q1Gna7HSzL8q+fOXMGu3fvRldXV0YCJCQWo9FIHTokYwQvvdVqNfr6+lBaWhr1eklJCUpKShAKhdDV1ZXQzB1CUs1oNKKlpQWjo6NU6o+knWCLMhQKTUiSY8nlcsFB4ek0NDSEysrKqAGiJP8YDAZwHIfe3l6xQyE5or29HZWVlRgaGprwnmCiHB4ennTHYweJZ0pxcTH27t1Ls3Hy3MKFC1FYWEiX3yRlamtrsXfv3qiZfxGCiTKRMmdUCo2IRSaTYenSpdShQzJCMFGq1eq463bv2LEDSqUyLUERkgij0Yienh5RbgGR/BK3etDy5cvhcrmwYcMGvvgFy7L8mMf9+/dnJkpCYjAajXjmmWcQCASmXZyFkHjiTmF8/fXXUV5ejrq6Ouh0Omi1WpjNZmi1WkqSRHRLly4FALr8Jmk36Vzv+vp6DA8PY2BgAH6/H+FwmBZzIlnhmmuugUajoQ4dknaCifLAgQNRzyPjJ8c6evRoWoIiJFEGg4FalCTt4s7MOXv2LEZHRwUfHR0dmYyVkAmMRiOOHDmCCxcuiB0KyWGChXtlMhkYhhH8YGS1vUyXuopUOKfK5gQAent7+WFCBoNB7HCIhEUqnceqcC7YorRYLBgYGMDw8HDMx8DAANavX5/24MejAedkrFtvvRWzZs2iy28ybfEGnAsOD6quro475EIul6O6ujo1ERIyRTNnzoRWq6UOHZJWgi3K8vLyST+cyDaEpBt16JB0izs86MCBA9i4cSP+/Oc/A7hUo3L58uWYPXs2HnrooYwESMhkjEYjBgcH8cknn4gdCslRgoly9+7dMJlM6O3thdlsxtGjR2GxWKBQKNDQ0ICBgQFs3Lgxk7ESElNkfW+6/CbpIniPsrOzEyMjI5DL5QgGg7BYLLDZbPzl9qZNm7Bq1aqMBUqIkBtuuAFz5sxBT08PVq5cKXY4JAcJtijLysogl8sBXFrTpqKiYsI9Sb1en97oCEkAwzBU8ZyklWCiLCoqinoea/3t8dsQIhaDwUBL2JK0iVuP8uzZs/zzUCgUc5tMowrnJBaj0YizZ8/iww8/FDsUIlHxKpwnPDMnMhNn/PNM/wteWVmJvXv3ZvSYJPudPXsWCoUCTqcTdXV1YodDJCxWjok7M2fsTJyRkZEJs3PEmJlDSCxXXXUVlixZQuMpSVoI9npbrVa+MyfeNoRki0jFc0JSTbBFGW8FxmS2ISRTDAYD3n///ah764SkwqSFezPNZrMhGAyKHQaRIKPRSEvYkrTIqkTp8/ngcDjEDoNI1KJFi2gJW5IWWZUoWZaNOV6TkETIZDLo9Xq6T0lSLmsSpdvthtlsFjsMInGRGTq0hC1JpaQS5datW1FTUwMA6O7uxujoaEqCCAaDUCgUKdkXyW8GgwEff/wxTp06JXYoJIcIDg8ar6GhARqNBiaTCcClWpRdXV247777ph1EZ2cnLBZLQttGZuZE0JIQZKzIchA9PT1QqVTiBkMkIbIERESsmTkJJ0q9Xo+qqip0d3enJrr/z+v1JlWFKLIUBCGxXHvttVCr1Th48CBWr14tdjhEAsY3tsY2xCISTpSDg4MAEDWNsbe3N2UtygiWZdHY2Iiamhpotdpp75vkH6p4TlIt4URZWlqKsrIyzJ49Gx6PB16vF3a7fdoBRC7lI6xWK6xWK/V+kykzGo3o6urChQsXMGvWLLHDITkg4c6c8vJydHZ2orS0FBzHwel04lvf+lbKAgkGg/wYSrvdDp/Pl7J9k/xiMBhw4cIF/P3vfxc7FJIjEm5Rbty4Edu3b8dzzz2XlkAUCgXq6+tRX1+flv2T/HHbbbdh5syZ6OnpwdKlS8UOh+SAhFuUHo8HO3bswIEDB9IZDyHTNmvWLFrClqRUwi3Kvr4+yOVyhEIh7NixA0qlEiaTCYWFhemMj5ApMRgMePXVV8UOg+SIhFuUkZJrcrkcHMehvr6e6lGSrGU0GsGyLP7zn/+IHQrJAQknypqaGmzcuBFKpRI+nw8ejwcdHR3pjC0mWgqCJCIy8Jwuv0mi4i0FkXCi7Ovrg06nw/DwMLZv346SkpKUBpmoyIBzmo1D4lGpVLj22mtpPCVJWG1tLfbu3Yvi4uIJ7yV8j9Jut6OqqiqlgRGSLrSELUmlhFuUsZJkIBBIZSyEpFRkCdtwOCx2KETiBFuUXV1dUb3aO3bsiHo/GAzC4/Fg//796Y2QkCkyGo0YHR3F8ePHsXjxYrHDIRIm2KLcsmULDh8+zD9vamrCyMgI/+A4DmfOnMlIkIRMRVlZGRiGofuUZNoEW5RjkyQAtLS0TFhMbPw8bUKySWFhIb+E7dq1a8UOh0hYwvcoi4qK+N9DoRB2794d9Roh2Yg6dEgqJJwovV4v/7tcLkdVVVXUa4Rko8gStufOnRM7FCJhcYcHhUIhdHZ2gmEYeDyeCe/39fVh3bp1aQsulsiAc6psThJhNBoRDodx+PBh3HXXXWKHQ7JYpNJ50hXO5XI5TCYT7HY7/H7/hEHmYlT6oQrnJBmLFi3ClVdeiZ6eHkqUJK5I42tKFc5LSkrQ1NSE7u5ulJeXpyVAQtKloKAAS5cupZ5vMi1JFe4d68CBA+jq6kp5QISkGi1hS6Yr4SmMwKVB6CzLAgA4jsPhw4dTsmYOIelkMBiwZcsWnD59GjfeeKPY4RAJSmq52mAwiOHhYajVagSDQVit1nTGRkhKjK0kRImSTEXCiVKj0WD9+vUYHBwEwzBQqVRU7ZxIwpw5c6BSqdDT05PU0siERCR8j1KtVuPUqVMoKSmB2+1OZ0yEpJzRaKQOHTJlCSfKYDAItVqN0dFRfPrpp1ixYgWam5vTGRshKWM0GuHz+fDFF1+IHQqRIIabYldgd3c3ysrK+CUiMkWn06G4uJgGnJOk9PT04Pbbb8ehQ4eg1+vFDodkobEDzvv6+qLeS7hFOV55eTkGBwenHVyyqMI5mYrS0lLMnDmT5n0TQVOqcD6+/uR4IyMj6OzsRG9v7/QjJCTNZs2ahdtuuw09PT145JFHxA6HSIxgomxqakJNTU3cD9MAXiIlRqMRf/rTn8QOg0iQYKK02+2TTlmkepRESoxGI15++WV8+umnuPrqq8UOh0iI4D3KWEly69atfCuzu7sbGo0mJUH4fD44HA44HA5UV1cjGAymZL+EjEVL2JKpSrgzp6GhAQqFgm9FlpeXp6wepdfrRX19Perr66HX66n4BkmLkpISXHPNNZQoSdISTpR6vR7r16+HWq1OaQA+nw+NjY38c7PZDJ/Px88pJyRVGIaBwWCggeckaQknyshQIIZh+NdS0eOt1WrR0tLCP49cdiuVymnvm5DxIpWEaAlbkoyE53qXlpairKwMs2fPhsfjgdfrhd1uT0kQZrOZ/72jowMmkwkKhSLmtpEK5xE08JwkI7KE7YkTJ7Bo0SKxwyFZIDLQPCJWhfOkZuYMDg7y0xZramomrMo4XcFgEDqdDn19fYKJsrKykiqckykbHR2FQqFAa2srHnzwQbHDIVkoVo5JamZOSUkJnnvuOTz33HMAJh+UniybzQaPxyOYJAmZrsLCQixevJg6dEhS4ibK0dFRdHV1IRAITHivtLQ0pcN4HA4HbDYbX+uShgiRdKEOHZIswUQ5ODgIlUoFs9kMjUaDQCCA0dFRbNiwAfPnz0dBQUHKpi+63W5otVo+SXZ2dlKrkqSN0WjEe++9R0vYkoQJJkq73Q6Xy4VwOIze3l40NTXx4xurqqrQ1NSEjo6OaQfAsiyqq6tRUVEBhmFQVFQEm8027f0SIiSyhO34CjGECBHs9dZoNHxi1Gq16O7uhtPpTHkHjlqtpjnjJKMWL17ML2F75513ih0OkQDBFuXY8ZLApWQ5PkkePXo0LUERkk4FBQXQ6/XUoUMSJpgoz5w5g7Nnz2J0dBSjo6MIhUL875FHqqYwEpJpBoMB77zzDl3NkITEvUepUChQVFSEoqIimM1m/vfI61SLkkiV0WjERx99hH/84x9ih0IkQDBRWiwWDA8P84+RkZGo38+cOYOSkpJMxgrg/2bmjB1JT0iyqJIQGa+9vR2VlZXJzcw5cuTIpB03iWyTajQzh6SKSqVCVVUVXnjhBbFDIVkkqZk5iSTATCdJQlIpUiCDkMlMeXExQqTOYDCgr6+PlrAlk6JESfKW0WjE+fPn8e6774odCslylChJ3iotLcWMGTPo8ptMihIlyVuXX345v4QtIfFQoiR5zWg0UqIkk6JESfKawWDAwMAAzpw5I3YoJItJLlHSgHOSSkajEQANPCfxB5xLLlEWFxdj7969tE4OSQm1Wo2rr76aEiVBbW0t9u7di+Li4gnvSS5REpJKtIQtSQQlSpL3aAlbMhlKlCTvGQwGhEIhnDx5UuxQSJaiREny3tKlS8EwDF1+E0GUKEnek8vlWLRoEXXoEEGUKAkBLWFL4qNESQgudei8++67+Oyzz8QOhWQhySVKGnBO0sFgMNAStnluShXOsxVVOCfp8OWXX0Iul2Pp0qWYO3cuVCoV1q5di/nz54sdGsmwWDlGcF1vQvLJ7373O3z++ed44403IJPJwDAMHA4HWltbsWbNGrHDIyKT3KU3IanW39+PdevW8c/D4TAuXryIcDiMuro6DAwMiBgdyQaUKEnea2trA8MwMd9jGAatra0Zjohkm6y49GZZFm63G2q1GizLwmKxQKFQiB0WyROBQABCt+ovXryIPXv2YNGiRViwYAEWLFgw5b/N/v5+tLW1IRAIpOUeqNT3n4ljTHn/XBbQarX8736/nzObzYLbrly5MhMhkTzS0NDAFRQUcABiPgoLC6Oez507l7vrrrs4q9XKvfjii9y+ffs4lmW5ixcvCh6jra2Nk8lkXEFBQdTPnTt3puQcpL7/TBwj0f3HyjGi93qzLIvq6uqoYRlFRUUYGRmJuT31epNU6+/vx8KFC2MWxZDJZDhx4gTmzp2LkydP4vjx4zhx4gSOHz+O48eP4+TJkzh//jyAS0tL3HTTTViwYAEWLlzIP2QyGXQ6Xdz9z5s3L63xZ/P+M3GMZPaflb3eXq8XSqUy6jWlUgmfzwetVitSVCSfzJ8/H62trairqwPDMOA4jv/Z2trKf4G0Wu2Ev8lwOIzTp0/ziTOSRP/617/io48+mvTYHMfh/vvvx5133slf/o/9Geu18e+99dZbgrcOOI7DqlWrYDAY+M+MfYzdl9BrR44cibv/73znO7j55pujPhMOh2MeT+gxMDAgWL0pHA7jrrvuwg033JD0fiPbf/LJJ4L7j9yHbmxsFPz/JHqiDAaDMV8fHh6O+XpkwHlEbW0tFfEl07ZmzRrccccdaG1t5e9f1dXVTdqKkclkUKlUUKlUuPvuu6PeC4VCOHHiBH74wx8KziPnOA7vvfcePvnkEzAMw3cqjf0Z67Wx7/3zn/+Mm8gCgQC/7fjPjn/Eeu+TTz6Ju/9QKIRz584J7jMy3Creo7+/P+5/5xkzZmDhwoWT7kfouPv27RNc7iMcDuO3v/0tjh07BgAxB5yLfo/SbrdzJpMp6jW1Ws25XK6Y29M9SiI18e6BFhQUcA0NDXm9/2w7h1g5RvThQQqFYkLrcXh4mHq9Sc5Yu3Zt3BZZXV1dXu8/E8eY7v5FT5Qmkynm62VlZRmOhJD0iNwDlclkKCgoiPo59h5ovu5fCucgeq83AOh0Or7Xm2VZWK1WeDyemNtSrzeRqoGBgaTvgebT/jNxjET2HzPHTOvCP0X8fj9XX1/PuVwurr6+nhsZGRHcVsx7lLt27RLt2Kki9XOQevwcJ/1zkHr8HBf/HLLyHiVwaclQu90Os9kMu92etfcnc6G0m9TPQerxA9I/B6nHDyR/DlmRKFMl2ZPPxP/wdMeU7nPIt/in8plsOwepxz+Vz6T7HChRplm+/ZFLPf6pfCbbzkHq8U/lM+k+h6zozEnGkiVLoNFoYr43NDSE4uLihPeV7u2zMSbaXvxj0PbiHyPe9n6/nx98HiG5REkIIZmWU5fehBCSDpQoCSFkEpQoCSFkEqJXD5ICn88Hr9cLAOjt7UVLS0vWjvWcjM1mw+bNmyUZv9frBcuyUKvVAISnv2YjlmX5koIsy8JsNvPnka18Ph/Wr18/YQlfKa1IIHQOSX+n0zXyPZfY7fao38dWZJeSvr4+DkDcmU/ZyuPxcBaLheO4SzO51Gq1yBElZ+zfEMdx/LlkK5fLxf+9jJfMigRiincOyX6nKVFOoq+vj1MoFPxzv9/PAeD8fr+IUU2Ny+Xi1Gq1JBPl+Lil9t9//Bcx2xNlxPgk4/f7J5zL2O9HNhp/DlP5TtM9yklotVq0tLTwzyOFhsdXZc92brcbZrNZ7DCmhGVZvvSez+dDMBjM+svW8ZRKJXQ6HX8JXlFRIXZIUxJvRQKpmMp3mhJlAsYmmI6ODphMpqy9JxNLMBiUVLzj+Xw+KJVK/r6Y0+mE2+0WO6ykuFwuAIBGo4HL5ZLsP1rJrkiQrZL9TlNnThKCwSDcbveEG8PZrrOzExaLRewwpmx4eBgsy/J/zBaLBUVFRYKFWLOR1+uF3W7nywgCQHNzs8hRpY5QAs12iX6nqUWZBJvNBo/HI6nWmdfrxapVq8QOY1rUajUUCgX/3z3yUyqXeyzLore3FyaTCRaLBX6/H52dnWBZVuzQkpZrKxIk+p2mFmWCHA4HbDYb1Go1/6+nVP44Ojs7+d9ZlkVjYyNqamoks8ql1O5Hjufz+aDX6/nnarUamzdvlmQrzGQyxWwJS3FFgqS+0+nqacolLpeL83g8HMdx3MjICNfc3CxyRFMHifbYm0wmvtdbasODIoWpxxr/PFshxnCy8cODxi8OmG1inUOy32kqijEJlmUnVCtSKBQYGRkRKaKpCQaDcDqdsNlssFgssFqtkmlRApfit9ls/LIhkZaAVHi9Xvh8Pr7FYjKZsjp+r9cLj8cDh8OB+vp66PV6vgOEZVk0NzdDr9ejt7c3aycwCJ3DVL7TlCgJIWQS1JlDCCGToERJCCGToERJCCGToERJCCGToERJCCGToERJCCGToERJCCGToERJ0i4y0D1S9cfr9cLpdGZ8rrPX64VOp4PT6ZzWfnQ6neSqF5HpobneJK0qKipQXV0Nu93Ov+bz+VBRUQG/35/RWEwmE2pqaqa9H7vdPmFus9RL2ZH4qEVJ0sbhcADAhBJvWq1W0mXfxtcuZFk2qvAIyT2UKEnaNDY28rUXx6uurs5wNOkztrVMchMlSpIWLMvGXbJhbFEIn88Hh8MBt9sNh8PB37sce08xcl+zuroaLMvynxmbiL1eLzQaDaxWK9xuN5xOJ6xWa9xyZl6vlz+2zWYDcGnZDJ1OxxdOcDgc0Gg0cDgc8Pl8Ufc5vV4vDh8+DI/HA6fTiV/96lfQaDT8PwTBYBAajYbfN5EmukdJRMWyLF88NUKn06G7uxsmkwkmkwkej4dfSsHlcsHtdqO+vh5arRYajYa/P2gymWA2mzF79my+0o3b7UZ1dXXU/scfO1Ldenh4mK80YzKZUF5ezu+7r6+Pv9wee58zEqNGo+FvJ8ycOZPfp0Kh4Cs2EemiFiVJi0hrUahnO9LibG5unlDuTa1W8/f8Zs+eHVX0VqFQRLVSY1XcHnv/0Gw2w+v1xmxVNjc3Q6lUwuv1Rq3xHNlHS0sLdDodysrKkuqosVgsfPwsy0qyqC2JRomSpE19fb3gujBjazMmK5W9y1qtlm8VWiwWvuUaOY5Wq0VHR0fC+4skZIvFAqfTCZ/PJ6m6nyQ2SpQkbSKdHOPHLY5t3dXU1PCtuQifzxd3nZ/JllAY+77b7Z7QSx15P9axI8+DwSC8Xi9cLhdYlp0wblIohsjnrVYrdfLkELpHSdIqUmHaZrNBo9HwaydH7iFqtVrY7XY4HA6o1Wr09vbC5XLxa3hHWnMmk4nvxIl8zuv1gmVZ2O122O12Phn6/X7+cjuyPwAT9hc5ts1m4y/vTSYTnE4n7HY731Gk1+uxfv16fiXIjo4OKJVKmM1mqNVqPik6nU6YTCYAl24faLVayS5LS6JRhXOSUyIJORs6T9xuNyXKHEEtSkJSyGq1orq6Gkqlku5N5hC6R0lyRqT32uVyibbmd3V1NYLBIIaHh7N68TCSHLr0JoSQSVCLkhBCJkGJkhBCJkGJkhBCJvH/AHJe92mxZqzsAAAAAElFTkSuQmCC\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, 7.343079250158649, 11.816843193856668)"
      ]
     },
     "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.80229304597699, 2.4003721135393445)"
      ]
     },
     "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.57350851e+00, 9.07779455e+00, 8.71183896e+00, 6.42715551e+00,\n",
       "        6.34985674e+00, 6.26905715e+00, 2.04323373e-01, 0.00000000e+00,\n",
       "        3.04715943e-03, 7.36748836e-03, 1.26161851e-02, 1.82100097e-02]),\n",
       " 'complexities': array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12]),\n",
       " 'min_bic_score': -88860.20954906201,\n",
       " 'found_coeff': array([ 0.93623127, -0.95866111, -0.89608068, -0.95636529, -0.89636987,\n",
       "         0.09433103,  0.0936821 ]),\n",
       " 'percent_coefficient_error_mean': 6.80229304597699,\n",
       " 'percent_coefficient_error_std': 2.4003721135393445}"
      ]
     },
     "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_V2.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}_V2.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, 9, 14, 15, 70)\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 70)\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 65, 70)\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, 15)\n",
      "(2, 4, 15)\n",
      "(2, 4)\n",
      "(2,)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 65, 70)\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 70)\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15, 70)\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, 15)\n",
      "(0,)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15, 70)\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, 70)\n",
      "--------------------------------------------------\n",
      "(0, 1, 8, 9, 15)\n",
      "(0, 1, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 9, 15)\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, 15)\n",
      "(2, 4, 15)\n",
      "(2, 4)\n",
      "(2,)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 65, 70)\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 70)\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15, 70)\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, 15)\n",
      "(0,)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15, 70)\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, 70)\n",
      "--------------------------------------------------\n",
      "(0, 1, 8, 9, 15)\n",
      "(0, 1, 9, 15)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 9, 15)\n",
      "--------------------------------------------------\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{1: ((0,), -79286.70104043813), 2: ((0, 15), -79782.41500301847), 3: ((0, 1, 15), -80148.37059246466), 4: ((0, 1, 9, 15), -82433.05404289016), 5: ((0, 1, 8, 9, 15), -82510.35281231243), 6: ((0, 1, 3, 5, 9, 15), -84201.47229495611), 7: ((1, 2, 3, 4, 5, 9, 15), -88655.88617596863), 8: ((0, 1, 2, 3, 4, 5, 9, 15), -88860.20954906201), 9: ((0, 1, 2, 3, 4, 5, 9, 14, 15), -88857.16238963649), 10: ((0, 1, 2, 3, 4, 5, 9, 14, 15, 70), -88852.8420606986), 11: ((0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 70), -88847.59336399431), 12: ((0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 65, 70), -88841.99953940697)}"
      ]
     },
     "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,) -79286.70104043813\n",
      "(0, 15) -79782.41500301847\n",
      "(0, 1, 15) -80148.37059246466\n",
      "(0, 1, 9, 15) -82433.05404289016\n",
      "(0, 1, 8, 9, 15) -82510.35281231243\n",
      "(0, 1, 3, 5, 9, 15) -84201.47229495611\n",
      "(1, 2, 3, 4, 5, 9, 15) -88655.88617596863\n",
      "(0, 1, 2, 3, 4, 5, 9, 15) -88860.20954906201\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15) -88857.16238963649\n",
      "(0, 1, 2, 3, 4, 5, 9, 14, 15, 70) -88852.8420606986\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 70) -88847.59336399431\n",
      "(0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 65, 70) -88841.99953940697\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/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAqZElEQVR4nO3df3Rbdf0/8OdNNkDQJU0Vp4XZ3rRQPqjHJqkTN3CHphtsWY+wdLUe0dKuycY+ih6h6cQfx4PaJeMIU5E2XTr88Vltk/XzYWwMlnSgAuK6ZPscnGxrcxvAyo+ztUkBAVlzv3/sm/tp2vxsk9zc5PU4J6dNcnPv60Lz2vu+7/v9ejM8z/MghBASl0zsAAghJN9RoiSEkCQoURJCSBKUKAkhJImcJUqfzwetVjvvdY7jYLPZ4HK5YLPZEAwGcxUSIYSkhMnFXW+XywWWZaHVajH3cFqtFl6vF8DFpGmxWOB0OrMdEiGEpCwniVI4GMNEJUqO49DY2CgkSgAoKSnB1NRUrkIihJCkloh5cI/HA5VKFfWaSqWCz+eDRqOJ+Znrr78earU65nsTExMoKytL+fjZ3j4fY6LtxT8GbS/+MRJt7/f7cerUqegX+Ryaezir1crr9fqo11iW5d1ud9x9aDQafuPGjcJj3759wnsbN25MK55sb5+LY9D2md0+F8eg7cU/xuzt9+3bF5VTNBrNvO1FbVHGk+iGTllZGQ4cOBDzvebm5rSOk+72C5HtmLJ9DsUW/0I+k2/nIPX4F/KZxZxDc3Nz1OcbGhrmb5RWGl6kuYfr6emZl72VSmXCFuVC/nXKFDGPnSlSPwepx8/z0j8HqcfP84nPIdZ7oo6j1Ov1MV/X6XQ5jiQ1uWiBZpvUz0Hq8QPSPwepxw+kfw45v+s9NTUFpVIpvDZ3eJDZbIbb7Y67j4aGhnmX3qOjo+jr60MgEEB5eTlaW1tRVVWVlXMghBS2WDkmJ32UHo9HSH5dXV2ora2F0WgEADidTlgsFtTW1mJkZCTtMZR79+7Fli1bhKFHDMPAZrPB4XCgpaUl06dCCClCOW1RZsLsbD86Oorq6mqEw+F528lkMpw5cwaVlZW5DpEQImGxWpSSnuvd19cHhmFivscwDBwOR44jIoQUIkknykAgMG9KZEQ4HMbZs2dzHBEhpBBJLlFOTEygoaEB/f39KC8vj9ui5Hkejz32GJqamnD48GFcuHAhx5ESQqSkv78fDQ0NmJiYmPdewfZRMgwDi8WCxx9/HKdOncInP/lJ3HHHHWhpaUF1dXWuwyaESETB9VFWVVXB4XBAJpNBLpdH/ezr60NXVxdefPFFjIyM4LbbboPdbsd1112HG264AT09PVTSjRCSEkm3KCPGxsbgcDiEcZRtbW0x73a///77OHDgAB599FE8+eSTuOSSS3DbbbehpaUFdXV1kMvluToNQkieipVjCiJRLsQ///lP/P73v8fevXtx+vRpXHXVVfj617+Ob3zjG7jmmmuitqUB7YQUj5g5JsNTKLMu0/NMw+Ew/8ILL/Bbt27lFQoFD4BftWoV39vby4dCIb6vr4+XyWS8XC6P+rl3796MxkEIyQ95N9c7HzAMg5UrV+KRRx7Ba6+9hv7+fnz4wx+GyWTCxz/+cbS2tiIcDmNmZibqZ1tbG8bGxsQOnxCSA0WfKGf70Ic+hK985St48skn8corr2DlypVxt6UB7YQUD0qUcVx11VX4xCc+EXec5szMDI4ePRpzzBUhpLBQokygvLwcMln8/0THjh3DVVddhc985jO455574Ha78d577+UwQkJILkguUc6emZNtra2tcadIymQy/PWvf0V/fz90Oh327duHtWvXQqVS4dZbb8VDDz2El156Ke7nCSH5pWBn5uTCo48+ira2tqgybjzPzyvjxvM8/va3v+Gpp57CU089hT/96U/497//jRUrVmDdunVYt24d6urqompxzkZDkAjJDzSOcoFSHdA+2zvvvIM//vGPQuI8c+YM5HI5Vq5cKSROnU4HuVwes6ZmrGRMCMk+SpQievnll4WkOTw8jFAoBJVKhS984Qs4fPhwzEt0qqlJSO4V3FxvKfnUpz4Fk8mE/fv349y5c3j22Wexfft2nDhxIm4/Jg1BIiQ/5OVytYVuyZIlWLVqFVatWoXR0VEMDAzETJY8zyMQCOQ+QEJIFGpRiizRECSGYVBeXp7bgAgh81CiFFmiIUg8z6OtrS3HERFC5qJEKbK5NTVnzwRyOBx0I4eQPCC5u95arRZlZWVobm4uiIXYI2YPQZqYmMALL7yAU6dO0VhKQnKkv78f/f39mJiYgNfrjXpPcolSqsOD0vHuu+/iuuuuw2c/+9mCP1dC8g0ND5KID33oQ9i1axcef/xxHDlyROxwCCl6lCjzlNFoxI033ojvfOc7tIIkISLLi0TJcRzsdjtcLhdsNhs4jhM7JNExDIPdu3fjpZdeQnd3t9jhEFLU8mLAucvlQkdHh/DcbDajp6dHxIjyQ01NDbZs2YIf/vCHaG5uRmlpqdghEVKU8qJFOTAwIHYIeesnP/kJZmZm8KMf/UjsUAgpWnmRKFUqFbRaLTiOg8fjQX19vdgh5Y0rr7wSP/zhD9Hd3Y2//e1vYodDSFFKmigDgQB27dqFtWvXora2VnisW7cODzzwQEbmIjudTgCAWq2G0+mE0Whc9D4LyTe/+U2wLIvvfOc7VAiYEBEkHEfZ2dkJhmGwefNm1NTUzHv/xIkTGBgYAMMw6OrqWnAQLpcLSqUSHMfBbDbDZDLF7aOMDDiPKLSB5/EcPHgQGzduxGOPPYaGhgaxwyGkYEQGmkfEGnAed11vm83GB4PBlNbBDQaDfGdnZ0rbzuX3+/mOjo6o50qlkvf7/TG3z/S63lIRDof5devW8Wq1mn/vvffEDoeQgpXWut733nsvFApFShlZoVAsuEXp8/lQW1srPGdZFjt27EAwGFzQ/goVwzD4+c9/jkAggN27d4sdDiFFRfSbORqNBiMjI1GvnT9/HhqNRqSI8td//Md/YPv27fjJT36C119/XexwCCka8xLl+Pg4tm3bhpMnTyIQCGDt2rWoqqrCyZMnsxIAy7Kor6+HzWaD3W6H3W6H2WzOyrEKwY9+9CMsXboU9913n9ihEFI05g047+npgV6vR3d3N0pKSmC1WqFQKNDb24vPfe5zWQlCr9dDr9dnZd+FRqVS4f7778d//ud/4q677oJWqxU7JEIK3rwWZX19PTZt2oTu7m7o9XrU1NSAZVlKZHnEZDLh+uuvx913303DhQjJgXmJMhgMYseOHQCAuro6AMADDzyAUCiU28hIXEuWLMFDDz2E5557jmY1EZID8xLlpk2b0NnZGfWaQqGgFmWeqaurw5e//GV0dHTgX//6l9jhEFLQYt71njssqLa2FsePH8fQ0BCOHj2K6enpnAQXy8TEBBoaGqIGiBarBx54AG+88QZ27doldiiESF5/fz8aGhowMTEx7724M3NCoRDa29vh8XhQWloq9IVNTk6CYRjodDo4nU4sW7Ysu9HPUQwVztPR2dmJX/ziFzhz5gyuvvpqscMhRPLSqnDe29sLq9WKyclJjI6OYmxsDGNjY5icnMT58+fxyCOPLGraIsmM++67D8uWLYPFYhE7FEIKVtxEWVFRgYqKirgfZFkWOp0uK0GR1H3kIx9BV1cX+vv78dxzz4kdDiEFKW6iTKXK+NwZNUQc3/jGN6DVanH33XcjHA6LHQ4hBSduhXOj0QidTgeGYaBSqYTXJycnhd8j5dGIuGQyGXbv3o3Vq1fjN7/5De68806xQyKkoMRNlBUVFTh+/DjGx8fh8XiEIhVKpRJ6vT7hZTnJvVWrVqG5uRk7duzApk2bcn6TjZBClnTNnIqKCrS3t+ciFrJIVqsV1157LX72s59h586dYodDSMFYVPWgoaGhTMVBMuDqq6+GxWLBgw8+CL/fL3Y4hBSMBSfK6elpUW7m0IDzxO699158/OMfxz333CN2KIRIyoIGnOt0Opw4cSLuTnmeB8MwmJmZyVykKaAB58n94Q9/QHNzMzwejzBfnxCSmrQHnB8/fhwzMzMxH+FwGPfee2/Wgybpa2pqwqpVq/Dtb38bFy5cEDscQiQvbqKsqanB+Ph4wg/TsrL5iWEY7N69G6dOnYLdbhc7HEIkL2Ef5e23357ww3RZl7+0Wi3uvPNO/OAHP4ga+0oISZ/oa+aQ7PnpT3+KDz74AD/+8Y/FDoUQSaNEWcCWL1+O73//+3j44Yfx97//XexwCJGslBLlnj170nqd5I+7774b5eXlMJvN6OzsFGbvjI6Oih0aIZKRUqLUarXYtm1b1Gvbtm2j6kEScOmll2L9+vV49tlnsWvXLgwODmLXrl2orq7Go48+mrHjjI6OYseOHZSISWHiU8RxHL9161ae53nebDbzPp8v1Y9m1MaNG0U5rlSdPXuWZxiGBzDvIZPJ+NHR0UUfo6+vj5fJZLxcLo/6uXfv3sWfACE5FivHxB1wHsuJEyeg1+sxPDyctaVrk9FqtSgrK0NzczOam5tFiUFKduzYgV27dsWcGMAwDNasWYOvfOUrUCgUMR9XXHEFGIaJu//R0VFUV1fHLO8mk8lw5swZVFZWZvScCMmG/v5+9Pf3Y2JiAl6vN+q9pEUxZuvq6sLg4CC6urpEW/2vrKyMZuakIRAIxF3Slud5PPvss3jmmWfibiOXy7Fs2bK4idTn88X9LMMwcDgcVAmfSEKk8dXQ0DDvvZQSZSgUgslkQm9vLxQKBbRaLZqamtDb20vlvPJceXl53BahXC7Hd7/7Xfz0pz/F22+/jVAolPQRDAYRCoUQCAQQCoXg9/sTJuJAIJDFsyMkR1K5ZrdYLHwwGIx6LRgM8p2dnYvrDFgA6qNMz9mzZ3mZTJa1PsrOzk5eLpfH3L9cLhflb4SQxYiVY1K6671z5855S9gqFIqMXlJ5PB7Y7XZ4PB54PJ6M7bfYVVVVweFwQCaTQS6XR/10OByL7j9sbW2N26IMh8Noa2tb1P4JyQdp9VFmi8fjgdPpRE9PDziOQ319PdVTzKCWlhasXr0aDocDgUAA5eXlaGtry8hNlkgibmtrA8MwQtIMh8NoamqiGzmkIMS86x0KhTA1NYXy8vKUdzQ9PQ0AC+qzVKvV8Hq9UCqVAC4ubMaybMxtqcxafhobG4tKxOfOncPevXtx9OhR3HTTTWKHR0jKYuWYmC1KhUKBwcFBlJaWJi2MAQD79+/H1NQUtmzZknZQHMdhcnISSqUSPp8PLMvGTZIkf1VWVkZ1xVy4cAFjY2PYvHkzfD4fPvnJT4oYHSGLE/fSu729HSdOnMDmzZuhVqtRW1sLlmWhVCoRDAbBcRyOHTuG8fFxmM1mbNq0aUEB+Hw+qFQquFwu6PV62O12sCwLo9G44JMi4luyZAn+8Ic/QKPRYPPmzXj66aexdOlSscMiZEFSGnAeCoUwODgIv9+PYDAIpVIJtVqdkdUY7XY7zGYzpqamhCRcUlIS9wZBZMB5BA08z2/PP/88vvSlL2H79u146KGHxA6HkHkiA80jYg04T3kKY7a43W5eqVRGvQaA93q9Mben4UHS88tf/pIHwPf394sdCiFJLXh4UDZRf2Th2759O7761a9iy5YtOHXqlNjhEJK2vEiUOp0OwWAQwP/d8dZoNOIGRjKGYRjY7XZUVFTg9ttvF0ZIECIVeTGO0ul0wmKxQKvVwuv1wu12ix0SybArrrgCQ0ND0Ol0uPPOO+FyuRIW2yAkn+RFolQqlejp6RE7DJJlVVVV+M1vfoPbbrsNDzzwAK3iSSQj7qV3IBCIekTQZRNZjC9/+cvo7OxEZ2cnnnnmGbHDISQlcROl2+0Gy7KwWq3gOE54/fz589i/fz+GhoZyEiApPPfffz/WrFmDpqYmTExMiB0OIUnFvfRmWRZerxc1NTVRr1dUVKCiogKhUAhDQ0MpzdwhZLYlS5agv78fWq0WjY2NeOaZZ3DJJZeIHRYhccVtUYZCoXlJcjaFQhF3UHg2TUxMoKGhIWqAKJGeK6+8Ei6XC8ePH8c999wjdjiEoL+/Hw0NDTGvcuImysnJyaQ7npqaWlxkCxCpcE6zcaRv5cqV2L17N375y1/iv/7rv8QOhxS55uZmHDhwIGrmX0TcRJlKmTMqhUYWa+vWrbjjjjtgMpnw4osvih0OITHFTZQsyyZct3vPnj1QqVRZCYoUD4Zh0N3djcrKStx+++0IhUJih0TIPAmrB61duxZOpxNbt24Vil9wHCeMeXzqqadyEyUpaJdffjn2798PnU6HlpYWDA0N0WB0klcSTmE8cuQI6urq0NbWBq1WC41GA6PRCI1GQ0mSZFRlZSV+97vf4X/+539gs9nEDoeQKEnnend0dGBychJjY2Pw+/0Ih8O0/CjJio0bN+K+++7D9773PQwPD4sdDiGCuIny6NGjUc8j4ydnO3nyZFaCIsXrxz/+Merq6tDc3IxXX31V7HAIAZBkZs5bb72F6enpuI+BgYFcxkqKgFwux759+3DZZZehsbER77//vtghERK/wrlMJkvYoc7zPBiGwczMTNaCiyVS4Zwqmxe2Y8eO4cYbb8SWLVvw8MMPix0OKQKRSuexKpzHbVGaTCaMjY1hcnIy5mNsbAzt7e1ZD34uGnBeHD7/+c/jF7/4BX7961/jd7/7ndjhkCKQaMB53OFBjY2NCdfDUSgUaGxszEyEhMRgMpnwwgsvwGw2Q6lU4vnnnxeWw21tbUVVVZXYIZIikdLiYvmE1vUuLu+++y6qq6vxyiuvQCa7eAHEMAx4nofD4UBLS4u4AZKCEyvHJBwedPToUWzbtg1PP/00gIs1KteuXYvS0lLcdddd2YuUkP/vH//4h3D3OxwOIxwOY2ZmBuFwGG1tbRgbGxM5QlIM4ibK/fv3Q6/XY2RkBEajESdPnoTJZIJSqURnZyfGxsawbdu2XMZKilBfX5/QkpyLYRg4HI4cR0SKUdw+ysHBQUxNTUGhUCAYDMJkMsFisaCurg4AcO+992Lz5s05C5QUp0AgELecH8/zUdX3CcmWuC1KnU4HhUIB4OKaNvX19UKSjKitrc1udKTolZeXxx2mxjAMysvLcxsQKUpxE2VJSUnU81jrb8/dhpBMa21tTdiibGtry3FEpBglrEf51ltvCc9jlb8Sox4lVTgvLlVVVXA4HJDJZJDL5ULrMtI/WVlZKXKEpFAkqnCe8sycyEycuc9zPTOHhgcVp7GxMTgcDvzv//4vDh8+DIfDgdbWVrHDIgUoVo6JezPHZDLBarXG3RnP8+js7MxcdIQkUFlZia6uLvA8D5ZlceLECbFDIkUk7qW32WyGQqGI+1AqlTCbzbmMlRAwDIMNGzbg4MGDoixuR4pT3ESZaAXGdLYhJNMMBgMCgQD+/ve/ix0KKRJJC/fmmsViQTAYFDsMksfWrFmDyy+/HIcOHRI7FFIk8ipR+nw+WgaAJHXZZZdBr9fj4MGDYodCikReJUqO42KO1yRkLoPBgOeeey6l9ecJWay8SZQulwtGo1HsMIhErF+/HuFwmBa5IzmRVqLctWsXmpqaAADDw8OYnp7OSBDBYBBKpTIj+yLFoaysDDU1NXT5TXIi7jjKuTo7O6FWq6HX6wEAdXV1GBoawu23377oIAYHB2EymVLaNjIzJ4KWhCheBoMBv/rVr3DhwgUsWZLynzIhUSJLQETEmpmT8l9XbW0tNm3alPFlRD0eT1pViCJLQRBiMBhw//334y9/+QtuvPFGscMhEjW3sTW7IRaR8qX3+Pg4AERNYxwZGVlMfILBwUHY7XbY7XZwHIeuri74fL6M7JsULp1OhyuvvJKGCZGsS7lFWVNTA51Oh9LSUrjdbng8noRTHFMVuZSPMJvNMJvNdPebJCWTybB+/XocPHgQO3fuFDscUsBSblHW1dVhcHAQNTU14HkedrsdN998c8YCCQaDwhhKq9VKLUqSEoPBgFOnTlEBX5JVKS8utm3bNjzyyCPZjicpqh5EZpuensZHP/pRPPjgg9i+fbvY4ZACkPbiYrO53W7s2bMHR48ezXhghCzUsmXLcNNNN9EwIZJVKfdRer1eKBQKhEIh7NmzByqVCnq9HsuWLctmfIQkZTAY0NnZiXfeeQdXXHGF2OGQApRyizKyfo5CoQDP8+jo6EB7e3vWAiMkVQaDAe+//37Gh64REpFyomxqasK2bdugUqng8/ngdrsxMDCQzdhioqUgyFyVlZW45ppr6PKbLEqipSBSTpRerxdarRaTk5N45JFHUFFRkdEgUxUZcE6zcchsBoOBivmSRWlubsaBAwdQVlY2772UE6XVasWWLVsyGhghmWIwGPDaa6/REhEkK1JOlJs2bZr3Go1dI/li9erVWLZsGc3SIVkR96730NBQ1F3tPXv2RL0fDAbhdrupzBXJC0uXLsW6detw8OBB/OAHPxA7HFJg4rYof/azn+H48ePC8+7ubkxNTQkPnudx/vz5nARJSCoMBgOOHTuGN954Q+xQSIGJ26KcnSQBoLe3d95iYnPnaRMipltvvRUMw+Dw4cNoaWkROxxSQFLuoywpKRF+D4VC2L9/f9RrhIjtYx/7GFauXEnDhEjGpZwoPR6P8LtCocCmTZuiXiMkHxgMBhw5cgT//ve/xQ6FFJCEUxhDoRAGBwfBMAzcbve8971eb86HDEUGnFNlcxKLwWDA97//ffz5z39GXV2d2OEQCYlUOo814Dxp9aDx8XFYrVYcP348Zu3IXA88p+pBJBGe57FixQoYjUY8+OCDYodDJChWjklaFKOiogLd3d0YHh6mf6FJ3mMYBhs2bMDjjz+On//851EV+QlZqLQK98529OhRDA0NZTwgQhbLYDDA7/fj7NmzYodCCkRaS9cNDQ2B4zgAFy9xjh8/npFVGAnJpJtvvhmXXXYZDh06hGuvvVbscEgBSLlF2dnZiSNHjuDYsWM4d+4c/H4/zGZzNmMjZEEuv/xy3HzzzTRMiGRMyi1KtVqN9vZ2jI+Pg2EYlJeXU7VzkrcMBgO+9a1vIRQKCbVUCVmolFuULMvi5ZdfRkVFBVwuVzZjImTRNmzYgAsXLuDIkSNih0IKQMqJMhgMgmVZTE9P49y5c1i3bh16enqyGRshC7ZixQp85jOfoctvkhEpr8I41/DwMHQ6Xc4va7RaLcrKymjAOUnqe9/7Hnp7e/H6669DLpeLHQ7Jc7MHnHu93qj3Um5RzlVXV4fx8fFFB5cuqnBOUmUwGHDu3DkcO3ZM7FCIBCSqcB73Zs7c+pNzTU1NYXBwECMjI4uPkJAsWLlyJUpLS3Ho0CHccMMNYodDJCxui3Ju/cm5DwC0PgnJa3K5HLfeeiv1U5JFi9uitFqtSacsUj1Kku8MBgN+//vf49VXX8XVV18tdjhEouK2KGMlyV27dqGpqQnAxZs5arU6I0H4fD7YbDbYbDY0NjYiGAxmZL+ErFu3DnK5HE888YTYoRAJS2tmjlKpFFqRdXV1GatH6fF40NHRgY6ODtTW1lLxDZIxSqUSq1evpstvsigpJ8ra2lq0t7eDZdmMBuDz+dDV1SU8NxqN8Pl8wpxyQhbLYDBgeHgY7777rtihEIlKOVFGhgLNLluViTveGo0Gvb29wvPIZbdKpVr0vgkBLibKd999F08//bTYoRCJSnmud01NDXQ6HUpLS+F2u+HxeGC1WjMShNFoFH4fGBiAXq+HUqmMuW2kwnkEDTwnyVx77bVgWRYHDx7E+vXrxQ6H5JnIQPOIBVU4n218fFyYttjU1DRvVcbFCgaD0Gq18Hq9cRMlVTgnC3H33Xfjv//7v/Hyyy9TMV+SUKwck9bMnIqKCuzcuRM7d+4EkHxQerosFgvcbnfcJEnIQhkMBrz66qt48cUXxQ6FSFDCRDk9PY2hoSEEAoF579XU1GR0GI/NZoPFYgHLsggGgzREiGTUTTfdhCuuuAKHDh0SOxQiQXET5fj4OMrLy2E0GqFWqxEIBDA9PY2tW7eiqqoKcrk8Y9MXXS4XNBqNkCQHBwepVUky6tJLL8XatWtpmBBZkLiJ0mq1wul0IhwOY2RkBN3d3cL4xk2bNqG7uxsDAwOLDoDjODQ2NqK+vh4Mw6CkpAQWi2XR+yVkLoPBgL/85S84d+6c2KEQiYl711utVguJUaPRYHh4GHa7PeM3cFiWpTnjJCfWr18Pnufx5JNP4mtf+5rY4RAJiduinHtnUKPRzEuSJ0+ezEpQhGTD8uXLodPp6PKbpC1uojx//jzeeustTE9PY3p6GqFQSPg98sjUFEZCcsVgMODJJ5/EBx98IHYoREIS9lEqlUqUlJSgpKQERqNR+D3yOtWiJFKzYcMGhEIhPP/882KHQiQkbqI0mUyYnJwUHlNTU1G/nz9/HhUVFbmMFcD/zcyZPZKekFRpNBosX76cLr/JPP39/WhoaEhvZs6JEyeS3rhJZZtMo5k5ZLG2bNmC5557Di+99JLYoZA8lNbMnFQSYK6TJCGZYDAYcPr0afj9frFDIRKx4MXFCJEqvV6PSy65hGbpkJRRoiRF58Mf/jDWrFlD/ZQkZZQoSVEyGAx45pln8NZbb4kdCpEASpSkKG3YsAEffPABjQUmKaFESYoSy7K47rrr6PKbpIQSJSlaBoMBhw4dQjgcFjsUkucklyhpwDnJlA0bNuCNN96Az+cTOxSSBxINOJdcoiwrK8OBAwdonRyyaF/84hehVCrp8psAuLj+1oEDB1BWVjbvPcklSkIyZenSpbjlllsoUZKkKFGSomYwGOD1evHaa6+JHQrJY5QoSVG75ZZbIJPJ8MQTT4gdCsljlChJUSstLcUNN9xAl98kIUqUpOgZDAa43W68//77YodC8hQlSlL0NmzYgHfeeQd//OMfxQ6F5ClKlKToffrTn8aKFSvo8pvEJblESQPOSaYxDAODwYCDBw/SiqBFbEEVzvMVVTgn2fDEE09gw4YN2LJlC95++22Ul5ejtbUVVVVVYodGcixWjom7rjchxeTVV18FADgcDjAMA4ZhYLPZ4HA40NLSIm5wRHSSu/QmJNNGR0dx1113AQB4nkc4HMbMzAzC4TDa2towNjYmcoREbJQoSdHr6+sDwzAx32MYBg6HI8cRkXyTF5feHMfB5XKBZVlwHAeTyQSlUil2WKRIBAKBuDdxZmZm8Nhjj+G6665DdXU1rr32WigUigUdZ3R0FH19fQgEAlnpA5X6/nNxjAXvn88DGo1G+N3v9/NGozHuths3bsxFSKSIdHZ28nK5nAcQ87Fs2bKo58uXL+fXrFnDm81m/sEHH+QPHz7Mj4+P8zMzM3GP0dfXx8tkMl4ul0f93Lt3b0bOQer7z8UxUt1/rBwj+l1vjuPQ2NgIr9crvFZSUoKpqamY29Ndb5Jpo6OjqK6ujlnAVyaT4cyZM1i+fDnOnj2L06dP48yZMzh9+jROnz6Ns2fP4r333gMAXHbZZbjmmmuElmd1dTWqq6shk8mg1WoT7r+ysjKr8efz/nNxjHT2n5d3vT0eD1QqVdRrKpUKPp8PGo1GpKhIMamqqoLD4UBbWxsYhgHP88JPh8MhfIE0Gs28v8lwOIxXXnlFSJyRJPqnP/0Jr7/+etJj8zyPO+64A1/60peEy//ZP2O9Nve9Z599Nm7XAc/z2Lx5M1auXCl8ZvZj9r7ivXbixImE+zcYDPj0pz8d9ZlwOBzzePEeY2NjcSvNh8NhrFmzBitWrEh7v5Ht33zzzbj7j/RDd3V1xf3/JHqiDAaDMV+fnJyM+XpkwHlEc3MzFfEli9bS0oLVq1fD4XAI/VdtbW1JWzEymQzl5eUoLy/HLbfcEvVeKBTCmTNn8K1vfQt//etfY36e53m8+OKLePPNN4VhSQCifsZ6bfZ7//jHPxImskAgIGw797NzH7Hee/PNNxPuPxQK4e233467T5lMFve9yGN0dDThf+elS5eiuro66X7iHffw4cM4f/58zH2Hw2H89re/xalTpwAg5oBz0fsorVYrr9fro15jWZZ3Op0xt6c+SiI1ifpA5XI539nZWdT7z7dziJVjRB8epFQq57UeJycn6a43KRitra0JW2RtbW1Fvf9cHGOx+xc9Uer1+piv63S6HEdCSHZE+kBlMhnkcnnUz9l9oMW6fymcg+h3vQFAq9UKd705joPZbIbb7Y65Ld31JlI1NjaWdh9oMe0/F8dIZf8xc8yiLvwzxO/38x0dHbzT6eQ7Ojr4qampuNuK2Ue5b98+0Y6dKVI/B6nHz/PSPwepx8/zic8hL/soAYBlWVitVhiNRlit1rztnyyE0m5SPwepxw9I/xykHj+Q/jnkRaLMlHRPPhf/w7MdU7bPodjiX8hn8u0cpB7/Qj6T7XOgRJllxfZHLvX4F/KZfDsHqce/kM9k+xzy4mZOOq6//nqo1eqY701MTKCsrCzlfWV7+3yMibYX/xi0vfjHSLS93+8XBp9HSC5REkJIrhXUpTchhGQDJUpCCEmCEiUhhCQhevUgKfD5fPB4PACAkZER9Pb25u1Yz2QsFgt27Nghyfg9Hg84jgPLsgDiT3/NRxzHCSUFOY6D0WgUziNf+Xw+tLe3R9WKBaS1IkG8c0j7O52tke+FxGq1Rv0+uyK7lHi9Xh5AwplP+crtdvMmk4nn+YszuViWFTmi9Mz+G+J5XjiXfOV0OoW/l7nSWZFATInOId3vNCXKJLxeL69UKoXnfr+fB8D7/X4Ro1oYp9PJsywryUQ5N26p/fef+0XM90QZMTfJ+P3+eecy+/uRj+aew0K+09RHmYRGo0Fvb6/wPFJoeG5V9nzncrlgNBrFDmNBOI4TSu/5fD4Eg8G8v2ydS6VSQavVCpfg9fX1Yoe0IIlWJJCKhXynKVGmYHaCGRgYgF6vz9s+mViCwaCk4p3L5/NBpVIJ/WJ2ux0ul0vssNLidDoBAGq1Gk6nU7L/aKW7IkG+Svc7TTdz0hAMBuFyueZ1DOe7wcFBmEwmscNYsMnJSXAcJ/wxm0wmlJSUxC3Emo88Hg+sVqtQRhAAenp6RI4qc+Il0HyX6neaWpRpsFgscLvdkmqdeTwebN68WewwFoVlWSiVSuG/e+SnVC73OI7DyMgI9Ho9TCYT/H4/BgcHwXGc2KGlrdBWJEj1O00tyhTZbDZYLBawLCv86ymVP47BwUHhd47j0NXVhaamJsmscim1/si5fD4famtrhecsy2LHjh2SbIXp9fqYLWEprkiQ1nc6W3eaConT6eTdbjfP8zw/NTXF9/T0iBzRwkGid+z1er1w11tqw4Mihalnm/s8XyHGcLK5w4PmLg6Yb2KdQ7rfaSqKkQTHcfOqFSmVSkxNTYkU0cIEg0HY7XZYLBaYTCaYzWbJtCiBi/FbLBZh2ZBIS0AqPB4PfD6f0GLR6/V5Hb/H44Hb7YbNZkNHRwdqa2uFGyAcx6Gnpwe1tbUYGRnJ2wkM8c5hId9pSpSEEJIE3cwhhJAkKFESQkgSlCgJISQJSpSEEJIEJUpCCEmCEiUhhCRBiZIQQpKgREmyLjLQPVL1x+PxwG6353yus8fjgVarhd1uX9R+tFqt5KoXkcWhud4kq+rr69HY2Air1Sq85vP5UF9fD7/fn9NY9Ho9mpqaFr0fq9U6b26z1EvZkcSoRUmyxmazAcC8Em8ajUbSZd/m1i7kOC6q8AgpPJQoSdZ0dXUJtRfnamxszHE02TO7tUwKEyVKkhUcxyVcsmF2UQifzwebzQaXywWbzSb0Xc7uU4z0azY2NoLjOOEzsxOxx+OBWq2G2WyGy+WC3W6H2WxOWM7M4/EIx7ZYLAAuLpuh1WqFwgk2mw1qtRo2mw0+ny+qn9Pj8eD48eNwu92w2+341a9+BbVaLfxDEAwGoVarhX0TaaI+SiIqjuOE4qkRWq0Ww8PD0Ov10Ov1cLvdwlIKTqcTLpcLHR0d0Gg0UKvVQv+gXq+H0WhEaWmpUOnG5XKhsbExav9zjx2pbj05OSlUmtHr9airqxP27fV6hcvt2f2ckRjVarXQnXDJJZcI+1QqlULFJiJd1KIkWRFpLca7sx1pcfb09Mwr98ayrNDnV1paGlX0VqlURrVSY1Xcnt1/aDQa4fF4YrYqe3p6oFKp4PF4otZ4juyjt7cXWq0WOp0urRs1JpNJiJ/jOEkWtSXRKFGSrOno6Ii7Lszs2ozpyuTdZY1GI7QKTSaT0HKNHEej0WBgYCDl/UUSsslkgt1uh8/nk1TdTxIbJUqSNZGbHHPHLc5u3TU1NQmtuQifz5dwnZ9kSyjMft/lcs27Sx15P9axI8+DwSA8Hg+cTic4jps3bjJeDJHPm81muslTQKiPkmRVpMK0xWKBWq0W1k6O9CFqNBpYrVbYbDawLIuRkRE4nU5hDe9Ia06v1ws3cSKf83g84DgOVqsVVqtVSIZ+v1+43I7sD8C8/UWObbFYhMt7vV4Pu90Oq9Uq3Ciqra1Fe3u7sBLkwMAAVCoVjEYjWJYVkqLdboderwdwsftAo9FIdllaEo0qnJOCEknI+XDzxOVyUaIsENSiJCSDzGYzGhsboVKpqG+ygFAfJSkYkbvXTqdTtDW/GxsbEQwGMTk5mdeLh5H00KU3IYQkQS1KQghJghIlIYQkQYmSEEKS+H9DiWk3h/dAaAAAAABJRU5ErkJggg==\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.2901852658927275, 11.816843193856668)"
      ]
     },
     "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.57350851e+00, 9.07779455e+00, 8.71183896e+00, 6.42715551e+00,\n",
       "        6.34985674e+00, 4.65873725e+00, 2.04323373e-01, 0.00000000e+00,\n",
       "        3.04715943e-03, 7.36748836e-03, 1.26161851e-02, 1.82100097e-02]),\n",
       " 'complexities': array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12]),\n",
       " 'min_bic_score': -88860.20954906201,\n",
       " 'found_coeff': array([ 0.93623127, -0.95866111, -0.89608068, -0.95636529, -0.89636987,\n",
       "         0.09433103,  0.0936821 ]),\n",
       " 'percent_coefficient_error_mean': 6.80229304597699,\n",
       " 'percent_coefficient_error_std': 2.4003721135393445}"
      ]
     },
     "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_V2.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}_V2.mat\"\n",
    "    print(file_name)\n",
    "    sio.savemat(file_name, metric_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e77fa604-381a-4851-96a9-50d06ba8ead7",
   "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
}