3D_exp/IPI-v.ipynb

{
 "cells": [
  {
   "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",
      "Sklearn's version: 1.0.2\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/\"\n",
    "\n",
    "np.random.seed(100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f428bf63-f8aa-4de3-97c8-dd17f45a8c55",
   "metadata": {},
   "outputs": [],
   "source": [
    "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([-1.25      , -1.23046875, -1.2109375 , -1.19140625, -1.171875  ,\n",
       "        -1.15234375, -1.1328125 , -1.11328125, -1.09375   , -1.07421875,\n",
       "        -1.0546875 , -1.03515625, -1.015625  , -0.99609375, -0.9765625 ,\n",
       "        -0.95703125, -0.9375    , -0.91796875, -0.8984375 , -0.87890625,\n",
       "        -0.859375  , -0.83984375, -0.8203125 , -0.80078125, -0.78125   ,\n",
       "        -0.76171875, -0.7421875 , -0.72265625, -0.703125  , -0.68359375,\n",
       "        -0.6640625 , -0.64453125, -0.625     , -0.60546875, -0.5859375 ,\n",
       "        -0.56640625, -0.546875  , -0.52734375, -0.5078125 , -0.48828125,\n",
       "        -0.46875   , -0.44921875, -0.4296875 , -0.41015625, -0.390625  ,\n",
       "        -0.37109375, -0.3515625 , -0.33203125, -0.3125    , -0.29296875,\n",
       "        -0.2734375 , -0.25390625, -0.234375  , -0.21484375, -0.1953125 ,\n",
       "        -0.17578125, -0.15625   , -0.13671875, -0.1171875 , -0.09765625,\n",
       "        -0.078125  , -0.05859375, -0.0390625 , -0.01953125,  0.        ,\n",
       "         0.01953125,  0.0390625 ,  0.05859375,  0.078125  ,  0.09765625,\n",
       "         0.1171875 ,  0.13671875,  0.15625   ,  0.17578125,  0.1953125 ,\n",
       "         0.21484375,  0.234375  ,  0.25390625,  0.2734375 ,  0.29296875,\n",
       "         0.3125    ,  0.33203125,  0.3515625 ,  0.37109375,  0.390625  ,\n",
       "         0.41015625,  0.4296875 ,  0.44921875,  0.46875   ,  0.48828125,\n",
       "         0.5078125 ,  0.52734375,  0.546875  ,  0.56640625,  0.5859375 ,\n",
       "         0.60546875,  0.625     ,  0.64453125,  0.6640625 ,  0.68359375,\n",
       "         0.703125  ,  0.72265625,  0.7421875 ,  0.76171875,  0.78125   ,\n",
       "         0.80078125,  0.8203125 ,  0.83984375,  0.859375  ,  0.87890625,\n",
       "         0.8984375 ,  0.91796875,  0.9375    ,  0.95703125,  0.9765625 ,\n",
       "         0.99609375,  1.015625  ,  1.03515625,  1.0546875 ,  1.07421875,\n",
       "         1.09375   ,  1.11328125,  1.1328125 ,  1.15234375,  1.171875  ,\n",
       "         1.19140625,  1.2109375 ,  1.23046875]),\n",
       " array([-1.25      , -1.23046875, -1.2109375 , -1.19140625, -1.171875  ,\n",
       "        -1.15234375, -1.1328125 , -1.11328125, -1.09375   , -1.07421875,\n",
       "        -1.0546875 , -1.03515625, -1.015625  , -0.99609375, -0.9765625 ,\n",
       "        -0.95703125, -0.9375    , -0.91796875, -0.8984375 , -0.87890625,\n",
       "        -0.859375  , -0.83984375, -0.8203125 , -0.80078125, -0.78125   ,\n",
       "        -0.76171875, -0.7421875 , -0.72265625, -0.703125  , -0.68359375,\n",
       "        -0.6640625 , -0.64453125, -0.625     , -0.60546875, -0.5859375 ,\n",
       "        -0.56640625, -0.546875  , -0.52734375, -0.5078125 , -0.48828125,\n",
       "        -0.46875   , -0.44921875, -0.4296875 , -0.41015625, -0.390625  ,\n",
       "        -0.37109375, -0.3515625 , -0.33203125, -0.3125    , -0.29296875,\n",
       "        -0.2734375 , -0.25390625, -0.234375  , -0.21484375, -0.1953125 ,\n",
       "        -0.17578125, -0.15625   , -0.13671875, -0.1171875 , -0.09765625,\n",
       "        -0.078125  , -0.05859375, -0.0390625 , -0.01953125,  0.        ,\n",
       "         0.01953125,  0.0390625 ,  0.05859375,  0.078125  ,  0.09765625,\n",
       "         0.1171875 ,  0.13671875,  0.15625   ,  0.17578125,  0.1953125 ,\n",
       "         0.21484375,  0.234375  ,  0.25390625,  0.2734375 ,  0.29296875,\n",
       "         0.3125    ,  0.33203125,  0.3515625 ,  0.37109375,  0.390625  ,\n",
       "         0.41015625,  0.4296875 ,  0.44921875,  0.46875   ,  0.48828125,\n",
       "         0.5078125 ,  0.52734375,  0.546875  ,  0.56640625,  0.5859375 ,\n",
       "         0.60546875,  0.625     ,  0.64453125,  0.6640625 ,  0.68359375,\n",
       "         0.703125  ,  0.72265625,  0.7421875 ,  0.76171875,  0.78125   ,\n",
       "         0.80078125,  0.8203125 ,  0.83984375,  0.859375  ,  0.87890625,\n",
       "         0.8984375 ,  0.91796875,  0.9375    ,  0.95703125,  0.9765625 ,\n",
       "         0.99609375,  1.015625  ,  1.03515625,  1.0546875 ,  1.07421875,\n",
       "         1.09375   ,  1.11328125,  1.1328125 ,  1.15234375,  1.171875  ,\n",
       "         1.19140625,  1.2109375 ,  1.23046875]),\n",
       " array([-1.25      , -1.23046875, -1.2109375 , -1.19140625, -1.171875  ,\n",
       "        -1.15234375, -1.1328125 , -1.11328125, -1.09375   , -1.07421875,\n",
       "        -1.0546875 , -1.03515625, -1.015625  , -0.99609375, -0.9765625 ,\n",
       "        -0.95703125, -0.9375    , -0.91796875, -0.8984375 , -0.87890625,\n",
       "        -0.859375  , -0.83984375, -0.8203125 , -0.80078125, -0.78125   ,\n",
       "        -0.76171875, -0.7421875 , -0.72265625, -0.703125  , -0.68359375,\n",
       "        -0.6640625 , -0.64453125, -0.625     , -0.60546875, -0.5859375 ,\n",
       "        -0.56640625, -0.546875  , -0.52734375, -0.5078125 , -0.48828125,\n",
       "        -0.46875   , -0.44921875, -0.4296875 , -0.41015625, -0.390625  ,\n",
       "        -0.37109375, -0.3515625 , -0.33203125, -0.3125    , -0.29296875,\n",
       "        -0.2734375 , -0.25390625, -0.234375  , -0.21484375, -0.1953125 ,\n",
       "        -0.17578125, -0.15625   , -0.13671875, -0.1171875 , -0.09765625,\n",
       "        -0.078125  , -0.05859375, -0.0390625 , -0.01953125,  0.        ,\n",
       "         0.01953125,  0.0390625 ,  0.05859375,  0.078125  ,  0.09765625,\n",
       "         0.1171875 ,  0.13671875,  0.15625   ,  0.17578125,  0.1953125 ,\n",
       "         0.21484375,  0.234375  ,  0.25390625,  0.2734375 ,  0.29296875,\n",
       "         0.3125    ,  0.33203125,  0.3515625 ,  0.37109375,  0.390625  ,\n",
       "         0.41015625,  0.4296875 ,  0.44921875,  0.46875   ,  0.48828125,\n",
       "         0.5078125 ,  0.52734375,  0.546875  ,  0.56640625,  0.5859375 ,\n",
       "         0.60546875,  0.625     ,  0.64453125,  0.6640625 ,  0.68359375,\n",
       "         0.703125  ,  0.72265625,  0.7421875 ,  0.76171875,  0.78125   ,\n",
       "         0.80078125,  0.8203125 ,  0.83984375,  0.859375  ,  0.87890625,\n",
       "         0.8984375 ,  0.91796875,  0.9375    ,  0.95703125,  0.9765625 ,\n",
       "         0.99609375,  1.015625  ,  1.03515625,  1.0546875 ,  1.07421875,\n",
       "         1.09375   ,  1.11328125,  1.1328125 ,  1.15234375,  1.171875  ,\n",
       "         1.19140625,  1.2109375 ,  1.23046875]),\n",
       " array([ 0.        ,  0.1010101 ,  0.2020202 ,  0.3030303 ,  0.4040404 ,\n",
       "         0.50505051,  0.60606061,  0.70707071,  0.80808081,  0.90909091,\n",
       "         1.01010101,  1.11111111,  1.21212121,  1.31313131,  1.41414141,\n",
       "         1.51515152,  1.61616162,  1.71717172,  1.81818182,  1.91919192,\n",
       "         2.02020202,  2.12121212,  2.22222222,  2.32323232,  2.42424242,\n",
       "         2.52525253,  2.62626263,  2.72727273,  2.82828283,  2.92929293,\n",
       "         3.03030303,  3.13131313,  3.23232323,  3.33333333,  3.43434343,\n",
       "         3.53535354,  3.63636364,  3.73737374,  3.83838384,  3.93939394,\n",
       "         4.04040404,  4.14141414,  4.24242424,  4.34343434,  4.44444444,\n",
       "         4.54545455,  4.64646465,  4.74747475,  4.84848485,  4.94949495,\n",
       "         5.05050505,  5.15151515,  5.25252525,  5.35353535,  5.45454545,\n",
       "         5.55555556,  5.65656566,  5.75757576,  5.85858586,  5.95959596,\n",
       "         6.06060606,  6.16161616,  6.26262626,  6.36363636,  6.46464646,\n",
       "         6.56565657,  6.66666667,  6.76767677,  6.86868687,  6.96969697,\n",
       "         7.07070707,  7.17171717,  7.27272727,  7.37373737,  7.47474747,\n",
       "         7.57575758,  7.67676768,  7.77777778,  7.87878788,  7.97979798,\n",
       "         8.08080808,  8.18181818,  8.28282828,  8.38383838,  8.48484848,\n",
       "         8.58585859,  8.68686869,  8.78787879,  8.88888889,  8.98989899,\n",
       "         9.09090909,  9.19191919,  9.29292929,  9.39393939,  9.49494949,\n",
       "         9.5959596 ,  9.6969697 ,  9.7979798 ,  9.8989899 , 10.        ]))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 128\n",
    "data = sio.loadmat(f\"../Datasets/reaction_diffusion_3d_{n}.mat\")\n",
    "u_sol = real(data['usol'])\n",
    "v_sol = real(data['vsol'])\n",
    "x = real(data['x'][0])\n",
    "y = real(data['y'][0])\n",
    "z = real(data['z'][0])\n",
    "t = real(data['t'][0])\n",
    "\n",
    "del data\n",
    "\n",
    "dt = t[1] - t[0]\n",
    "dx = x[1] - x[0]\n",
    "dy = y[1] - y[0]\n",
    "dz = z[1] - z[0]\n",
    "# For WeakPDELibrary\n",
    "ground_indices_u = (0, 1, 5, 9, 15, 23)\n",
    "ground_coeff_u = np.array([0.014, -0.014, -1.000, 0.020, 0.020, 0.020])\n",
    "ground_indices_v = (2, 5, 10, 16, 24)\n",
    "ground_coeff_v = np.array([-0.067, 1.0, 0.01, 0.01, 0.01])\n",
    "\n",
    "x, y, z, t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "51a0a7f1-5302-46c7-9974-0466b1de2282",
   "metadata": {},
   "outputs": [],
   "source": [
    "time = t\n",
    "u = np.zeros((n, n, n, len(time), 2))\n",
    "u[:, :, :, :, 0] = u_sol\n",
    "u[:, :, :, :, 1] = v_sol\n",
    "# u_dot = ps.FiniteDifference(axis=3)._differentiate(u, dt)\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",
    "X, Y, Z, T = np.meshgrid(x, y, z, time, indexing=\"ij\")\n",
    "spatiotemporal_grid = np.asarray([X, Y, Z, T])\n",
    "spatiotemporal_grid = np.transpose(spatiotemporal_grid, axes=[1, 2, 3, 4, 0])\n",
    "weak_lib = ps.WeakPDELibrary(\n",
    "    library_functions=library_functions,\n",
    "    function_names=library_function_names,\n",
    "    derivative_order=2,\n",
    "    spatiotemporal_grid=spatiotemporal_grid,\n",
    "    is_uniform=True,\n",
    "    include_interaction=False,\n",
    "    include_bias=True,\n",
    "    periodic=True,\n",
    "    K=10000,\n",
    "    cache=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7ad8b85d-f21a-491d-ad42-d1c604f86f66",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(100)\n",
    "noise_lvs = [0, 0.04, 0.1]\n",
    "noise_lv = noise_lvs[-1]\n",
    "domain_noise = 0.01*np.abs(noise_lv)*np.std(u)*np.random.randn(*u.shape)\n",
    "u = u + domain_noise\n",
    "# np.save(\"u_noise1e-1.npy\", u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "1e8c3c5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "denoise = True\n",
    "if denoise: \n",
    "    u = np.load(\"u_noise1e-1_denoised.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "186a6ed9-e4ce-4734-a664-f27af8c1fdd1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(x0)' = 0.014 1 + -0.014 x0 + -1.000 x0x1x1 + 0.020 x0_33 + 0.020 x0_22 + 0.020 x0_11\n",
      "(x1)' = -0.067 x1 + 1.000 x0x1x1 + 0.010 x1_33 + 0.010 x1_22 + 0.010 x1_11\n"
     ]
    }
   ],
   "source": [
    "include_sr3 = True\n",
    "# Fit the 3D reaction-diffusion equation\n",
    "optimizer = ps.SR3(threshold=1e-5, normalize_columns=True, max_iter=100)\n",
    "model = ps.SINDy(feature_library=weak_lib, optimizer=optimizer, cache=True)\n",
    "model.fit(u, ensemble=False, library_ensemble=True)\n",
    "model.print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "cb0bae31-7357-4fce-8174-273997affaaa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hasattr threshold in optimizer...\n",
      "x0_t = (0.013997 +0.000000i)1\n",
      "    + (-0.013965 +0.000000i)x0\n",
      "    + (-1.000266 +0.000000i)x0x1x1\n",
      "    + (0.019998 +0.000000i)x0_33\n",
      "    + (0.020000 +0.000000i)x0_22\n",
      "    + (0.020004 +0.000000i)x0_11\n",
      "   \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0, 1, 5, 9, 15, 23)"
      ]
     },
     "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": "e00d2a8f-1d40-41c9-9e18-4e1d0b512c1e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hasattr threshold in optimizer...\n",
      "x1_t = (-0.067002 +0.000000i)x1\n",
      "    + (1.000066 +0.000000i)x0x1x1\n",
      "    + (0.009997 +0.000000i)x1_33\n",
      "    + (0.009999 +0.000000i)x1_22\n",
      "    + (0.010000 +0.000000i)x1_11\n",
      "   \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(2, 5, 10, 16, 24)"
      ]
     },
     "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": "41ce20d8-be68-4bf6-b9df-91fcf1e226bf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "({(0, 1, 5, 9, 15, 23),\n",
       "  (0, 3, 5, 9, 15, 23),\n",
       "  (1, 2, 3, 4, 5, 6, 9, 10, 15, 16, 21, 22, 23, 24)},\n",
       " {(0, 1, 2, 3, 4, 5, 6, 7, 9, 10, 15, 16, 17, 20, 21, 22, 23),\n",
       "  (1, 2, 5, 6, 9, 16, 23, 24),\n",
       "  (1, 3, 4, 5, 6, 9, 10, 15, 16, 22, 23, 24),\n",
       "  (2, 5, 10, 16, 24)})"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# pub 32 v\n",
    "# ({(0, 1, 2, 4, 6, 9, 10, 15, 16, 22, 23, 24),\n",
    "#   (0, 1, 5, 9, 15, 23),\n",
    "#   (0, 1, 5, 10, 15, 23),\n",
    "#   (0, 3, 5, 9, 15, 23)},\n",
    "#  {(0, 1, 2, 3, 4, 5, 9, 10, 17, 21, 24),\n",
    "#   (0, 1, 2, 4, 6, 10, 15, 16, 22, 23, 24),\n",
    "#   (1, 2, 3, 4, 5, 6, 10, 15, 16, 21, 22),\n",
    "#   (1, 4, 5, 6, 10, 16, 22, 24),\n",
    "#   (2, 5, 10, 16, 24)})\n",
    "\n",
    "sr3_u = set()\n",
    "sr3_u.add(ensemble_indices_u)\n",
    "for i in range(len(model.get_coef_list())):\n",
    "    sr3_u.add(nonz(model.get_coef_list()[i][0]))\n",
    "sr3_v = set()\n",
    "sr3_v.add(ensemble_indices_v)\n",
    "for i in range(len(model.get_coef_list())):\n",
    "    sr3_v.add(nonz(model.get_coef_list()[i][1]))\n",
    "sr3_u, sr3_v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "8f196144-5590-49d3-9cda-677a696e2da7",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_pre, y_pre = model.feature_library.cached_xp_full[0], model.cached_x_dot\n",
    "X_pre = np.load(\"./IPI_output_files/X_pre_noise1e-1_denoised.npy\")\n",
    "y_pre = np.load(\"./IPI_output_files/y_pre_noise1e-1_denoised.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "480c677f-954a-4bd7-a89b-dec21c1c956a",
   "metadata": {},
   "outputs": [],
   "source": [
    "max_complexity = 12; target_index = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "b379cbbb-331e-425f-8b36-192b7c1498e3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{(2,), (2, 5, 10, 16, 24), (2, 5, 16, 24), (5, 16, 24)}"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# pub 32 v\n",
    "# {(2,), (2, 5, 10, 16, 24), (2, 5, 16, 24), (5, 16, 24)}\n",
    "\n",
    "nor_axis = None\n",
    "Xbnb = X_pre.copy()\n",
    "if nor_axis is not None: Xbnb = normalize(Xbnb, axis=nor_axis)\n",
    "\n",
    "thres = 1e-4\n",
    "bnb_sols = bnb(Xbnb, y_pre[:, target_index:target_index+1].ravel(), max_nonzeros=max_complexity, \n",
    "               lam=1e-2, threshold=thres, normalize=True, corrected_coefficients=True)\n",
    "bnb_feature_hierarchy = 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": 15,
   "id": "e79134d8-a02a-451c-9760-40070e2f0c23",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{(2,),\n",
       " (2, 3, 5, 8, 9, 10, 11, 16, 17, 18, 23, 24),\n",
       " (2, 5, 8, 9, 10, 11, 16, 17, 18, 23, 24),\n",
       " (2, 5, 8, 9, 10, 11, 16, 18, 23, 24),\n",
       " (2, 5, 8, 10, 11, 16, 18, 23, 24),\n",
       " (2, 5, 10, 11, 16, 18, 23, 24),\n",
       " (2, 5, 10, 11, 16, 23, 24),\n",
       " (2, 5, 10, 16, 23, 24),\n",
       " (2, 5, 10, 16, 24),\n",
       " (2, 5, 16, 24),\n",
       " (2, 5, 24),\n",
       " (2, 24)}"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nor_axis = None\n",
    "Xfrols = X_pre.copy()\n",
    "if nor_axis is not None: Xfrols = normalize(Xfrols, axis=nor_axis)\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": 16,
   "id": "9f940b54-aa89-45a2-8286-fa0355578344",
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_hierarchy = bnb_feature_hierarchy.union(selected_inds_list)\n",
    "if include_sr3:\n",
    "    if target_index == 0:\n",
    "        feature_hierarchy = feature_hierarchy.union(sr3_u)\n",
    "    elif target_index == 1:\n",
    "        feature_hierarchy = feature_hierarchy.union(sr3_v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "56e5dee8-db5e-4c3e-a9e6-74240ff0db54",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RFE_PLinearRegression\n",
      "(2, 5, 24)\n",
      "(2, 24)\n",
      "(2,)\n",
      "--------------------------------------------------\n",
      "(2, 5, 8, 9, 10, 11, 16, 17, 18, 23, 24)\n",
      "(2, 5, 8, 9, 10, 11, 16, 18, 23, 24)\n",
      "(2, 5, 8, 10, 11, 16, 18, 23, 24)\n",
      "(2, 5, 10, 11, 16, 18, 23, 24)\n",
      "(2, 5, 10, 11, 16, 23, 24)\n",
      "(2, 5, 10, 16, 23, 24)\n",
      "(2, 5, 10, 16, 24)\n",
      "(2, 5, 16, 24)\n",
      "(5, 16, 24)\n",
      "(16, 24)\n",
      "(24,)\n",
      "--------------------------------------------------\n",
      "(1, 2, 5, 6, 9, 16, 23, 24)\n",
      "(2, 5, 6, 9, 16, 23, 24)\n",
      "(2, 5, 9, 16, 23, 24)\n",
      "(2, 5, 16, 23, 24)\n",
      "(2, 5, 16, 24)\n",
      "--------------------------------------------------\n",
      "(2, 5, 8, 9, 10, 11, 16, 18, 23, 24)\n",
      "--------------------------------------------------\n",
      "(2, 3, 5, 8, 9, 10, 11, 16, 17, 18, 23, 24)\n",
      "(2, 5, 8, 9, 10, 11, 16, 17, 18, 23, 24)\n",
      "--------------------------------------------------\n",
      "(2, 5, 10, 16, 24)\n",
      "--------------------------------------------------\n",
      "(2, 5, 16, 24)\n",
      "--------------------------------------------------\n",
      "(2, 5, 8, 10, 11, 16, 18, 23, 24)\n",
      "--------------------------------------------------\n",
      "(2, 5, 10, 16, 23, 24)\n",
      "--------------------------------------------------\n",
      "(2, 5, 10, 11, 16, 23, 24)\n",
      "--------------------------------------------------\n",
      "(1, 3, 4, 5, 6, 9, 10, 15, 16, 22, 23, 24)\n",
      "(1, 3, 4, 5, 6, 9, 10, 16, 22, 23, 24)\n",
      "(1, 4, 5, 6, 9, 10, 16, 22, 23, 24)\n",
      "(1, 4, 5, 6, 10, 16, 22, 23, 24)\n",
      "(1, 4, 5, 6, 10, 16, 22, 24)\n",
      "(1, 4, 5, 6, 10, 16, 24)\n",
      "(1, 4, 5, 6, 16, 24)\n",
      "(4, 5, 6, 16, 24)\n",
      "(4, 5, 16, 24)\n",
      "(5, 16, 24)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 6, 7, 9, 10, 15, 16, 17, 20, 21, 22, 23)\n",
      "(0, 1, 2, 3, 4, 5, 6, 7, 9, 10, 15, 16, 17, 20, 21, 22)\n",
      "(0, 1, 2, 3, 4, 5, 6, 7, 10, 15, 16, 17, 20, 21, 22)\n",
      "(0, 1, 2, 3, 4, 5, 6, 10, 15, 16, 17, 20, 21, 22)\n",
      "(1, 2, 3, 4, 5, 6, 10, 15, 16, 17, 20, 21, 22)\n",
      "(1, 2, 3, 4, 5, 6, 10, 15, 16, 20, 21, 22)\n",
      "(1, 2, 3, 4, 5, 6, 10, 15, 16, 21, 22)\n",
      "(1, 2, 3, 4, 5, 6, 10, 16, 21, 22)\n",
      "(1, 2, 4, 5, 6, 10, 16, 21, 22)\n",
      "(1, 2, 4, 5, 10, 16, 21, 22)\n",
      "(1, 2, 5, 10, 16, 21, 22)\n",
      "(1, 2, 5, 10, 16, 22)\n",
      "(1, 2, 5, 10, 16)\n",
      "(1, 2, 5, 16)\n",
      "(2, 5, 16)\n",
      "(2, 5)\n",
      "(2,)\n",
      "--------------------------------------------------\n",
      "(5, 16, 24)\n",
      "--------------------------------------------------\n",
      "(2, 24)\n",
      "--------------------------------------------------\n",
      "(2, 5, 10, 11, 16, 18, 23, 24)\n",
      "--------------------------------------------------\n",
      "SelectKBest\n",
      "(2, 5, 24)\n",
      "(2, 24)\n",
      "(2,)\n",
      "--------------------------------------------------\n",
      "(2, 5, 8, 9, 10, 11, 16, 17, 18, 23, 24)\n",
      "(2, 5, 8, 9, 10, 11, 16, 18, 23, 24)\n",
      "(2, 5, 8, 10, 11, 16, 18, 23, 24)\n",
      "(2, 5, 10, 11, 16, 18, 23, 24)\n",
      "(2, 5, 10, 11, 16, 23, 24)\n",
      "(2, 5, 10, 16, 23, 24)\n",
      "(2, 5, 10, 16, 24)\n",
      "(2, 5, 16, 24)\n",
      "(5, 16, 24)\n",
      "(16, 24)\n",
      "(24,)\n",
      "--------------------------------------------------\n",
      "(1, 2, 5, 6, 9, 16, 23, 24)\n",
      "(2, 5, 6, 9, 16, 23, 24)\n",
      "(2, 5, 9, 16, 23, 24)\n",
      "(2, 5, 16, 23, 24)\n",
      "(2, 5, 16, 24)\n",
      "--------------------------------------------------\n",
      "(2, 5, 8, 9, 10, 11, 16, 18, 23, 24)\n",
      "--------------------------------------------------\n",
      "(2, 3, 5, 8, 9, 10, 11, 16, 17, 18, 23, 24)\n",
      "(2, 5, 8, 9, 10, 11, 16, 17, 18, 23, 24)\n",
      "--------------------------------------------------\n",
      "(2, 5, 10, 16, 24)\n",
      "--------------------------------------------------\n",
      "(2, 5, 16, 24)\n",
      "--------------------------------------------------\n",
      "(2, 5, 8, 10, 11, 16, 18, 23, 24)\n",
      "--------------------------------------------------\n",
      "(2, 5, 10, 16, 23, 24)\n",
      "--------------------------------------------------\n",
      "(2, 5, 10, 11, 16, 23, 24)\n",
      "--------------------------------------------------\n",
      "(1, 3, 4, 5, 6, 9, 10, 15, 16, 22, 23, 24)\n",
      "(1, 3, 4, 5, 6, 9, 10, 16, 22, 23, 24)\n",
      "(1, 4, 5, 6, 9, 10, 16, 22, 23, 24)\n",
      "(1, 4, 5, 6, 10, 16, 22, 23, 24)\n",
      "(1, 4, 5, 6, 10, 16, 22, 24)\n",
      "(1, 4, 5, 6, 10, 16, 24)\n",
      "(1, 4, 5, 6, 16, 24)\n",
      "(4, 5, 6, 16, 24)\n",
      "(4, 5, 16, 24)\n",
      "(5, 16, 24)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 6, 7, 9, 10, 15, 16, 17, 20, 21, 22, 23)\n",
      "(0, 1, 2, 3, 4, 5, 6, 7, 9, 10, 15, 16, 17, 20, 21, 22)\n",
      "(0, 1, 2, 3, 4, 5, 6, 7, 10, 15, 16, 17, 20, 21, 22)\n",
      "(0, 1, 2, 3, 4, 5, 6, 10, 15, 16, 17, 20, 21, 22)\n",
      "(1, 2, 3, 4, 5, 6, 10, 15, 16, 17, 20, 21, 22)\n",
      "(1, 2, 3, 4, 5, 6, 10, 15, 16, 20, 21, 22)\n",
      "(1, 2, 3, 4, 5, 6, 10, 15, 16, 21, 22)\n",
      "(1, 2, 3, 4, 5, 6, 10, 16, 21, 22)\n",
      "(1, 2, 4, 5, 6, 10, 16, 21, 22)\n",
      "(1, 2, 4, 5, 10, 16, 21, 22)\n",
      "(1, 2, 5, 10, 16, 21, 22)\n",
      "(1, 2, 5, 10, 16, 22)\n",
      "(1, 2, 5, 10, 16)\n",
      "(1, 2, 5, 16)\n",
      "(2, 5, 16)\n",
      "(2, 5)\n",
      "(2,)\n",
      "--------------------------------------------------\n",
      "(5, 16, 24)\n",
      "--------------------------------------------------\n",
      "(2, 24)\n",
      "--------------------------------------------------\n",
      "(2, 5, 10, 11, 16, 18, 23, 24)\n",
      "--------------------------------------------------\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{1: ((2,), -205294.33276816175), 2: ((16, 24), -208548.38635326718), 3: ((5, 16, 24), -225673.94037683957), 4: ((2, 5, 16, 24), -254978.96449502304), 5: ((2, 5, 10, 16, 24), -325867.846506729), 6: ((2, 5, 10, 16, 23, 24), -325874.43433331477), 7: ((2, 5, 10, 11, 16, 23, 24), -325880.3525063712), 8: ((2, 5, 10, 11, 16, 18, 23, 24), -325882.1715341), 9: ((2, 5, 8, 10, 11, 16, 18, 23, 24), -325877.40728853544), 10: ((2, 5, 8, 9, 10, 11, 16, 18, 23, 24), -325872.07208175457), 11: ((2, 5, 8, 9, 10, 11, 16, 17, 18, 23, 24), -325865.9927872089), 12: ((2, 3, 5, 8, 9, 10, 11, 16, 17, 18, 23, 24), -325859.614889683), 13: ((1, 2, 3, 4, 5, 6, 10, 15, 16, 17, 20, 21, 22), -214355.2010601697), 14: ((0, 1, 2, 3, 4, 5, 6, 10, 15, 16, 17, 20, 21, 22), -214354.82395822284), 15: ((0, 1, 2, 3, 4, 5, 6, 7, 10, 15, 16, 17, 20, 21, 22), -214351.3346953505), 16: ((0, 1, 2, 3, 4, 5, 6, 7, 9, 10, 15, 16, 17, 20, 21, 22), -214345.46636742094), 17: ((0, 1, 2, 3, 4, 5, 6, 7, 9, 10, 15, 16, 17, 20, 21, 22, 23), -214339.54408288823)}"
      ]
     },
     "execution_count": 17,
     "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": 18,
   "id": "1d053d3e-7f83-4370-82dd-c4c53df3587a",
   "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"
     ]
    }
   ],
   "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": 19,
   "id": "80daf1ca-5b16-45d6-98e6-09e8d397a55e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2,) -205294.33276816175\n",
      "(16, 24) -208548.38635326718\n",
      "(5, 16, 24) -225673.94037683957\n",
      "(2, 5, 16, 24) -254978.96449502304\n",
      "(2, 5, 10, 16, 24) -325867.846506729\n",
      "(2, 5, 10, 16, 23, 24) -325874.43433331477\n",
      "(2, 5, 10, 11, 16, 23, 24) -325880.3525063712\n",
      "(2, 5, 10, 11, 16, 18, 23, 24) -325882.1715341\n",
      "(2, 5, 8, 10, 11, 16, 18, 23, 24) -325877.40728853544\n",
      "(2, 5, 8, 9, 10, 11, 16, 18, 23, 24) -325872.07208175457\n",
      "(2, 5, 8, 9, 10, 11, 16, 17, 18, 23, 24) -325865.9927872089\n",
      "(2, 3, 5, 8, 9, 10, 11, 16, 17, 18, 23, 24) -325859.614889683\n"
     ]
    }
   ],
   "source": [
    "# pub 32 v\n",
    "# (2,) -225753.91368340503\n",
    "# (16, 24) -228763.57914524354\n",
    "# (5, 16, 24) -244641.6767333153\n",
    "# (2, 5, 16, 24) -267260.97739774553\n",
    "# (2, 5, 10, 16, 24) -272890.6729010754\n",
    "# (2, 5, 10, 16, 17, 24) -272914.162795415\n",
    "# (2, 3, 5, 10, 16, 17, 24) -272924.30835923535\n",
    "# (2, 3, 5, 7, 10, 16, 17, 24) -272932.6275170798\n",
    "# (0, 2, 3, 5, 7, 10, 16, 17, 24) -272938.3679217328\n",
    "# (0, 2, 3, 5, 7, 10, 16, 17, 23, 24) -272943.1201767921\n",
    "# (0, 2, 3, 4, 5, 7, 10, 16, 17, 23, 24) -272939.5732545059\n",
    "# (0, 2, 3, 4, 5, 7, 10, 16, 17, 22, 23, 24) -272934.92400972766\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": 20,
   "id": "2780a396",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_subsets = [br.track[com][0] for com in complexities]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "d502d5a2",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████████████| 12/12 [00:02<00:00,  5.57it/s]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[(2,),\n",
       " (16, 24),\n",
       " (5, 16, 24),\n",
       " (2, 5, 16, 24),\n",
       " (2, 5, 10, 16, 24),\n",
       " (2, 5, 10, 16, 23, 24),\n",
       " (2, 5, 10, 11, 16, 23, 24),\n",
       " (2, 5, 10, 11, 16, 18, 23, 24),\n",
       " (2, 5, 8, 10, 11, 16, 18, 23, 24),\n",
       " (2, 5, 8, 9, 10, 11, 16, 18, 23, 24),\n",
       " (2, 5, 8, 9, 10, 11, 16, 17, 18, 23, 24),\n",
       " (2, 3, 5, 8, 9, 10, 11, 16, 17, 18, 23, 24)]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ever = reduce(np.union1d, best_subsets)\n",
    "_, bs = brute_force_all_subsets(X_pre[:, ever], y_pre[:, target_index:target_index+1], max_support_size=len(best_subsets))\n",
    "best_subsets = [tuple(ever[list(b)]) for b in bs]\n",
    "del bs\n",
    "best_subsets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "7a1b877c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.save(\"./IPI_output_files/X_pre_noise1e-1.npy\", X_pre)\n",
    "# np.save(\"./IPI_output_files/y_pre_noise1e-1.npy\", y_pre)\n",
    "# np.save(\"./IPI_output_files/feature_names.npy\", model.get_feature_names(), allow_pickle=True)\n",
    "# np.save(\"./IPI_output_files/best_subsets.npy\", best_subsets, allow_pickle=True)\n",
    "\n",
    "np.save(\"./IPI_output_files/v/feature_names.npy\", model.get_feature_names(), allow_pickle=True)\n",
    "np.save(\"./IPI_output_files/v/best_subsets_denoised.npy\", best_subsets, allow_pickle=True)\n",
    "\n",
    "# np.load(\"./IPI_output_files/X_pre_noise1e-1.npy\")\n",
    "# np.load(\"./IPI_output_files/y_pre_noise1e-1.npy\")\n",
    "# np.load(\"./IPI_output_files/feature_names.npy\", allow_pickle=True)\n",
    "# np.load(\"./IPI_output_files/best_subsets.npy\", allow_pickle=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "1c73c56f",
   "metadata": {},
   "outputs": [],
   "source": [
    "BICs = []\n",
    "uncertainties = []\n",
    "for com in complexities:\n",
    "    effective_indices, _ = br.track[com]\n",
    "    fit_res = sm.OLS(y_pre[:, 1:2], X_pre[:, effective_indices]).fit()    \n",
    "    BICs.append(fit_res.bic)\n",
    "    uncertainties.append(sum(fit_res.bse))\n",
    "BICs = np.array(BICs)\n",
    "uncertainties = np.array(uncertainties)/min(uncertainties)\n",
    "lam = 4\n",
    "UBICs = BICs+(10**lam)*log(len(y_pre)*uncertainties)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "b7444cf4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.argmin(BICs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "532d955e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABAAAAAKTCAYAAACHP3pUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABwwElEQVR4nO3de1yUdd7/8fdwGlABUVTEI9pBSd0UWgMjtFY0a29t28Tc2Lw7bJSuIttvXbXa8r7TLNe2KLE27bDtqm2sm23aDZZ5WMnURVLDTfMAiqiYMh45zu8PnJERRFDgmsPr+XjMQ+fiM9d8mPXucX/f1/f6jMlqtVoFAAAAAADcmpfRDQAAAAAAgOZHAAAAAAAAgAcgAAAAAAAAwAMQAAAAAAAA4AEIAAAAAAAA8AAEAAAAAAAAeAACAAAAAAAAPICP0Q24m6qqKhUWFiowMFAmk8nodgAAAAAAbs5qterUqVMKDw+Xl9flr/MTADSxwsJCdevWzeg2AAAAAAAepqCgQF27dr3szwkAmlhgYKCk6g8+KCjI4G4AAAAAAO7OYrGoW7du9vXo5RAANDHbtv+goCACAAAAAABAi7nSbegMAQQAAAAAwAMQAAAAAAAA4AEIAAAAAAAA8AAEAAAAAAAAeAACAAAAAAAAPAABAAAAAAAAHoAAAAAAAAAAD0AAAAAAAACAByAAAAAAAADAAxAAAAAAAADgAQgAAAAAAADwAAQAAAAAAAB4AAIAAAAAAAA8AAEAAAAAAAAegAAAAAAAAAAPQAAAAAAAAIAHIAAAAAAAAMADEAAAAAAAAOABCAAAAAAAAPAABAAAAAAAAHgAQwOA7777TqNHj1ZoaKiCgoI0ZMgQrVmzxqEmPz9fP/3pT9W6dWuFhoZq8uTJKisrc6jZvn274uPjFRAQoC5dumjWrFmyWq0ONWvXrlVUVJT8/f3Vq1cvLVy4sFY/GRkZioyMlNlsVmRkpJYvX970vzQAAAAAAAYwNAC4++67VVFRoS+++EJbt27VzTffrHvuuUdFRUWSpMrKSt199906c+aMNmzYoKVLlyojI0O/+c1v7OewWCwaPny4wsPDtXnzZqWlpWnevHmaP3++vWbfvn0aNWqU4uLilJOToxkzZmjy5MnKyMiw12RnZysxMVFJSUnKzc1VUlKSxo4dq02bNrXcBwIAAAAAQDMxWS+9VN5CiouL1aFDB61bt05xcXGSpFOnTikoKEirV6/WnXfeqVWrVumee+5RQUGBwsPDJUlLly7VhAkTdPToUQUFBSk9PV3Tp0/XkSNHZDabJUkvvvii0tLSdPDgQZlMJk2bNk0rVqxQXl6e/f2Tk5OVm5ur7OxsSVJiYqIsFotWrVplrxk5cqRCQkK0ZMmSy/4epaWlKi0ttT+3WCzq1q2bSkpKFBQU1HQfGAAAAAAAdbBYLAoODr7iOtSwHQDt27dX37599f777+vMmTOqqKjQm2++qU6dOikqKkpS9VX5fv362Rf/kjRixAiVlpZq69at9pr4+Hj74t9WU1hYqP3799trEhISHN5/xIgR2rJli8rLy+ut2bhxY72/x5w5cxQcHGx/dOvW7eo+kBZ2rqyy1m0SAAAAAAD3ZVgAYDKZlJWVpZycHAUGBsrf31+vvPKKPvvsM7Vt21aSVFRUpE6dOjm8LiQkRH5+fvbbBOqqsT2/Uk1FRYWKi4vrrbGd43KmT5+ukpIS+6OgoKARn4IxTpdWaNyfvtL//DNPVVWEAAAAAADgCZo8AHjuuedkMpnqfWzZskVWq1VPPvmkOnbsqPXr1+vrr7/W6NGjdc899+jw4cP285lMplrvYbVaHY5fWmO7st0UNXW9f01ms1lBQUEOD2e3Yfcx5Rac1OJ/7dNvM75RRWWV0S0BAAAAAJqZT1OfcNKkSRo3bly9NT179tQXX3yhf/7znzpx4oR90bxgwQJlZWXpvffe0+9+9zuFhYXVGsJ34sQJlZeX26/Wh4WF1bpKf/ToUUm6Yo2Pj4/at29fb82luwLcwch+nfWH+3+k32Z8o4+2HpTlXLlee2Cg/H29jW4NAAAAANBMmnwHQGhoqPr06VPvw9/fX2fPnq1uwMuxBS8vL1VVVV+RjomJ0Y4dOxx2BGRmZspsNtvnBMTExGjdunUOXw2YmZmp8PBw9ezZ016TlZXl8D6ZmZmKjo6Wr69vvTWxsbFN8Kk4n/uiuir9F4Pk5+OlzG+P6OF3N+t0aYXRbQEAAAAAmolhMwBiYmIUEhKihx56SLm5ufruu+/0//7f/9O+fft09913S5ISEhIUGRmppKQk5eTk6PPPP9dTTz2lxx57zL5rYPz48TKbzZowYYJ27Nih5cuXa/bs2UpNTbVv309OTtaBAweUmpqqvLw8LV68WIsWLdJTTz1l72fKlCnKzMzU3LlztWvXLs2dO1erV69WSkpKi382LSXhpjC9+9+3qLWftzZ+f1y/+NNXOnGm7MovBAAAAAC4HMMCgNDQUH322Wc6ffq07rjjDkVHR2vDhg36+OOP9aMf/UiS5O3trU8//VT+/v4aMmSIxo4dqzFjxmjevHn28wQHBysrK0sHDx5UdHS0nnzySaWmpio1NdVeExERoZUrV+rLL7/UzTffrP/5n//Ra6+9pvvuu89eExsbq6VLl+qdd97RgAED9O6772rZsmUaPHhwy30oBojtHaq/PnarQlr5Kvdgica+ma2ikvNGtwUAAAAAaGImK98F16Qa+v2LzmbP0VN68O2vVWQ5ry5tA/SXRwerZ2hro9sCAAAAAFxBQ9ehhu0AgHO5rmOgPnoiRj3bt9Khk+f084XZ+rbQYnRbAAAAAIAmQgAAu64hrfS35Fj17Ryk4tOlGvdWtrYe+MHotgAAAAAATYAAAA46BJq19Fe3KrpHiCznK/SLtzdp7XfHjG4LAAAAAHCNCABQS3CAr/78yGDF39BB58ur9Oh7m/XPbwqNbgsAAAAAcA0IAFCnAD9v/emX0bpnQGeVV1r16yU5WvJ1vtFtAQAAAACuEgEALsvPx0uvjhuo8YO7y2qVpv99uxau/d7otgAAAAAAV4EAAPXy9jLphTH99MTQ3pKkF1ft0ourdolvjwQAAAAA10IAgCsymUyaNrKPfndXH0nSwrXfa8byHaqsIgQAAAAAAFdBAIAGS47vrTk/6y+TSVrydb4mL81RWUWV0W0BAAAAABqAAACN8sCPu+v1BwbJ19ukT785rMfe36JzZZVGtwUAAAAAuAICADTa3QM66+2HblGAr7fWfndMDy7apJKz5Ua3BQAAAACoBwEArkr8DR30waM/VpC/j7YeOKHEt7J19NR5o9sCAAAAAFwGAQCuWlSPdlr2eIxC25i1q+iUxi7MVsEPZ41uCwAAAABQBwIAXJO+nYP0UXKMuoYEaP/xs/r5wo3afeSU0W3BYBWVVXr9i90aPHu1/ralwOh2AAAAAIgAAE2gZ2hrfZQcq+s7ttERS6nufzNbuQUnjW4LBtl95JR+lr5R8zK/0xFLqV5YmafTpRVGtwUAAAB4PAIANImwYH99+HiMftStrU6eLdf4P32ljd8XG90WWlBllVUL136vu1/boG8OlijI30dhQf46ebZc72fvN7o9AAAAwOMRAKDJhLT2018eHazY3u11pqxSE97ZrMydRUa3hRbw/bHT+vnCjXpx1S6VVVZp2I0dlJUar9+OvFGS9Kd1e3WGXQAAAACAoQgA0KTamH20eMItSojspLKKKj3xl38rY+tBo9tCM6mssurt9Xs16tX1ysk/qUCzj176+QAtnnCLOgX5679+FK6e7VvpxNlyffDVAaPbBQAAADwaAQCanL+vtxb8YpB+HtVVlVVW/eZvuVq8YZ/RbaGJ7S8+o3FvZet/P81TaUWV4q4P1f9NvV1jo7vJZDJJkny8vTRx2HWSpLfW7dXZMnYBAAAAAEYhAECz8PH20kv3DdDDQyIkSbP++a1eyfpOVqvV4M5wraqqrHpv437d9ep6bd5/Qq39vDX73v56/+EfK7xtQK36MQO7qHu7Vjp+pkx/3ZRvQMcAAAAAJAIANCMvL5OeuaevfjP8BknSq5/v1vOffKuqKkIAV1Xww1n94u1N+v2KnTpXXqmYXu31WcrtGj+4u/2q/6V8vb00cVhvSdLCtXt1rqyyJVsGAAAAcAEBAJqVyWTSr++8Xs//102SpHc37tdTf8tVeWWVwZ2hMaxWq/66KV8j/7hO2XuPK8DXW8//1036y6OD1a1dqyu+/meDuqpL2wAVny7Vkq/ZBQAAAAAYgQAALeKh2J76Y+LN8vYy6e85h/TEB//W+XKuBLuCwpPn9MvFX2vG8u06U1apW3qGaNWUOD0U21NeXnVf9b+Ub41ZAAvXfs//9gAAAIABCADQYsYM7KI3H4yS2cdLq/OOaMI7X+vU+XKj28JlWK1Wfbi5QCNeWaf1u4tl9vHS03f31dJfxahnaOtGn+/nUV0VHuyvo6dKtWxzQTN0DAAAAKA+BABoUT+J7KT3Hv6x2ph99NXeHzT+T5t0/HSp0W3hEkcs5/Xwu5v124xvdKq0QgO7t9XKKXF6NK6XvBt41f9Sfj5eeuLCLoD0L79XaQW7AAAAAICWRACAFndrr/Za+qtb1a61n7YfKtHYN7NVePKc0W1B1Vf9l+cc1PD5a7XmP8fk5+2l393VRx8lx6p3hzbXfP6x0V0VFuSvIst5fbjlYBN0DAAAAKChCABgiH5dgvXh4zEKD/bX98fO6P6F2dp77LTRbXm0Y6dK9as/b9XUZbmynK/QgK7B+ufk25Qc3/uqr/pfyuzjrSeGVn8jQPqaPewCAAAAAFoQAQAMc13HNvrbE7HqFdpah06e09g3s7WzsMTotjzSJ7mFSnhlrbK+PSJfb5OeSrhBf38iVjd0Cmzy90q8pZs6BppVWHJeGVsPNfn5AQAAANSNAACG6tI2QB8mx+im8CAVny7TuDe/0ub9Pxjdlsc4frpUE//yb/16SY5OnC1XZOcgrZh0mybdcb18vJvnPw/+vt5Kjq/eBfDGmj0qq+ArIQEAAICWQAAAw4W2MWvJr27Vj3u206nSCiUt2qQ1u44a3Zbb+2zHYSW8sk6fbj8sHy+Tptx5vf4xcYj6dg5q9vceP7i7QtuYdejkOS3PYRYAAAAA0BIIAOAUgvx99d7DP9awGzvofHmVHnt/i1bkFhrdlls6caZMU5bmKPmDf+v4mTLd2ClQ/5g4RFOH3yA/n5b5T0L1LoBekqTX1+xReSW7AAAAAIDmRgAApxHg5623fhmt//pRuCqqrJqyNEcffHXA6Lbcyupvjyjhj+v08bZCeZmkJ4f21opfD1G/LsEt3ssvBvdQaBs/FfxwTv/IYRYAAAAA0NwIAOBUfL299MfEm/Xgrd1ltUpP/2OH3lizR1ar1ejWXFrJuXL95sNcPfr+Fh07VareHVrr708O0W9H9pHZx9uQngL8vPVY3MVdABXsAgAAAACaFQEAnI6Xl0n/M7qfJg27TpL08v/9Ry+u2kUIcJW+/M9RjXhlnTL+fVAmk/Sr23vp08lxurlbW6Nb04O39lC71n46cPwst3wAAAAAzYwAAE7JZDLpqRE3auaovpKkN9ft1e8ytquyihCgoU6dL9fvMr7RhHc2q8hyXhGhrfW3x2M0Y1Rf+fsac9X/Uq3NPno0LkKS9PoXe/jfFwAAAGhGBABwao/d3ksv3TdAXiZp2ZYC/XrJv1VaUWl0W07vX3uKNfKP67V0c4Ek6b+H9NTKyXGK7tnO4M5q+2VMT7Vt5au9xWf0z2/YBQAAAAA0FwIAOL2xt3TTG+MHyc/bSyu3F+nR97bobFmF0W05pTOlFXr6H9v1i7c36dDJc+rWLkBLf3Wrfv/TmxTg5xxX/S/VxuyjR2+r3gXw2ue72QUAAAAANBMCALiEu/p31qIJ0Wrl5631u4v14NubVHK23Oi2nMpXe49r5Kvr9MFX+ZKkpFt76LMpt+vWXu0N7uzKHortqSB/H31/7IxWbj9sdDsAAACAWyIAgMuIu76DPnh0sIIDfPXv/JNKfCtbRy3njW7LcOfKKvX8Jzs17q2vVPDDOXVpG6APHhms/xnTT63NPka31yCB/r565LbqbwRI+2K3qtgFAAAAADQ5AgC4lEHdQ/Th4zHqGGjWrqJT+vnCbOUfP2t0W4bZsv8H3fXqOr3zr/2SpHG3dNNnKXG67fpQYxu7ChOG9FSgv4++O3Jan+0sMrodAAAAwO0QAMDl3BgWqI+SY9W9XSvl/3BWP1+4Uf8pOmV0Wy3qfHmlXvj0W93/Zrb2Hz+rsCB/vfvft+jF+wYo0N/X6PauSnCAr/57yMVZAOwCAAAAAJoWAQBcUvf2rfRRcoxu7BSoo6dKNfbNbOXknzC6rRaRk39Cd7+2Xn9av09Wq3TfoK76v6m3a+iNHY1u7Zo9PKSn2ph9tKvolDK/PWJ0OwAAAIBbIQCAy+oY5K9lj9+qgd3bquRcuX7x9iZt2F1sdFvNprSiUi99tkv3pW/U98fOqEOgWX/6ZbT+MPZHCg5wzav+l2rbyk8TYntKqt4FYLWyCwAAAABoKgQAcGltW/npg0cGK+76UJ0tq9TD727WZzvcb4r8jkMl+q+0f2nBl9+ryiqNvjlcmSm3a3hkJ6Nba3KP3Bah1n7e+vawRavzjhrdDgAAAOA2CADg8lqbffT2Q9G6q1+Yyiqr9ORf/q0PtxQY3VaTKKuo0vys7zT6jX/pP0dOqX1rPy18cJBeHTdQIa39jG6vWYS09tMv2QUAAAAANDkCALgFs4+30h4YqLHRXVVllX770Td6e/1eo9u6Jt8WWjTmjX/ptc93q7LKqlH9w5Q59XaN7NfZ6Naa3aO3RSjA11vbD5VozX/YBQAAAAA0BQIAuA0fby/NvW+AHourniT/v5/m6Q+Z/3G5K8jllVVK+3y3Rr+xQd8etqhtK1+lPTBQC34RpfZtzEa31yLatzHrlzE9JEmvfr7H5f43BAAAAJwRAQDcislk0oxRffX/RtwoSUr7Yo9+v2Kny3yl3HdHTulnCzbqD1nfqbzSquGRnZQ59Xb99EfhRrfW4h6N6yV/Xy/lFpzU2u+OGd0OAAAA4PIIAOB2TCaTJg67Tv8zpp9MJun97AOa+uE2lVdWGd3aZVVWWZX+5fe657UN2n6oREH+Pnol8Ud6KylKHQP9jW7PEB0CzXpwsG0XALMAAAAAgGtFAAC3lXRrD/0x8Wb5eJn08bZCPf7nrTpfXml0W7V8f+y0fr5wo+Z+tktllVUadmMHZaXG696BXWUymYxuz1C/iu8ls4+XcvJPasMe9/2KRwAAAKAlEADArY2+uYve+mWUzD5e+mLXUf1y8deynC83ui1J1Vf9316/V6NeXa+c/JMKNPvopZ8P0OIJt6hTkGde9b9Ux0B/jR/cXZL06mp2AQAAAADXggAAbu+OPp3050cGK9Dso6/3/aDxf/pKx0+XGtrT/uIzGvdWtv730zyVVlQp7vpQfTb1do2N7ubxV/0vlRzfW34+Xtpy4ISyvz9udDsAAACAyyIAgEf4cUQ7LfnVrWrf2k87Dll0/5vZOnTyXIv3UVVl1Xsb9+uuV9dr8/4Tau3nrdn39tf7D/9YXdoGtHg/rqBTkL8euKWbpOpZAAAAAACuDgEAPEa/LsH6W3KMurQN0N5jZ3R/+kZ9f+x0i71/wQ9n9Yu3N+n3K3bqXHmlYnq112cpt2v84O5c9b+C5KG95eftpU37ftBXe9kFAAAAAFwNAgB4lF4d2uhvyTHq3aG1CkvO6/6F2dpxqKRZ39Nqteovmw5o5B/XKXvvcQX4euv5/7pJf3l0sLq1a9Ws7+0uOgcHaOwtXSVJr7ELAAAAALgqBADwOOFtA/Th4zHq3yVYP5wp07i3vtKmZrqqXHjynH65+GvNXL5DZ8oqdUvPEK2aEqeHYnvKy4ur/o3xxNDr5Ott0sbvj2vz/h+MbgcAAABwOQQA8Ejt25j118cGa3BEO50urdAvF3+tL3YdabLzW61Wfbi5QCNeWaf1u4tl9vHS03f31dJfxahnaOsmex9P0qVtgH4eVT0LgF0AAAAAQOMRAMBjBfr76r2Hf6yf9O2o0ooq/er9rfp426FrPu8Ry3k9/O5m/TbjG50qrdDA7m21ckqcHo3rJW+u+l+TJ4f2lo+XSet3F2vrgRNGtwMAAAC4FAIAeDR/X2+lPxilewd2UUWVVSnLtunP2fuv6lxWq1XLcw5q+Py1WvOfY/Lz9tLv7uqjj5Jj1btDm6Zt3EN1a9dK9w1iFgAAAABwNQgA4PF8vb30h/t/pIdieshqlZ75eKfSPt8tq9Xa4HMcO1WqX/15q6Yuy5XlfIX6dwnWPyffpuT43lz1b2ITh10nby+T1n53TNsKThrdDgAAAOAyCAAASV5eJj33Xzdp8p3XS5L+kPWdXvg0r0EhwCe5hUp4Za2yvj0iX2+Tnkq4QX9/MlY3dAps7rY9Uvf2rXTvwC6S2AUAAAAANAYBAHCByWRS6vAb9Mw9kZKktzfs028/+kYVlVV11h8/XaqJf/m3fr0kRyfOliuyc5BWTLpNk+64Xr7e/J9Wc5o47Dp5maQvdh3VNwdPGt0OAAAA4BKadZXywgsvKDY2Vq1atVLbtm3rrMnPz9dPf/pTtW7dWqGhoZo8ebLKysocarZv3674+HgFBASoS5cumjVrVq0rs2vXrlVUVJT8/f3Vq1cvLVy4sNZ7ZWRkKDIyUmazWZGRkVq+fHmtmgULFigiIkL+/v6KiorS+vXrr/4DgEt65LYIvfzzAfIySX/belCT/pqj0opKh5rPdhxWwivr9On2w/LxMmnKndfrHxOHqG/nIIO69iwRoa015mbbLoA9BncDAAAAuIZmDQDKysp0//3364knnqjz55WVlbr77rt15swZbdiwQUuXLlVGRoZ+85vf2GssFouGDx+u8PBwbd68WWlpaZo3b57mz59vr9m3b59GjRqluLg45eTkaMaMGZo8ebIyMjLsNdnZ2UpMTFRSUpJyc3OVlJSksWPHatOmTfaaZcuWKSUlRTNnzlROTo7i4uJ01113KT8/vxk+HTiz+6O7Kf3BKPl5e+mznUV6+N3NOlNaoRNnyjR5SY6SP/i3jp8p042dAvWPiUM0dfgN8vPhqn9LmnhH9S6A1XlHtONQidHtAAAAAE7PZG3MpLOr9O677yolJUUnT550OL5q1Srdc889KigoUHh4uCRp6dKlmjBhgo4ePaqgoCClp6dr+vTpOnLkiMxmsyTpxRdfVFpamg4ePCiTyaRp06ZpxYoVysvLs587OTlZubm5ys7OliQlJibKYrFo1apV9pqRI0cqJCRES5YskSQNHjxYgwYNUnp6ur2mb9++GjNmjObMmVPn71ZaWqrS0lL7c4vFom7duqmkpERBQVwNdnX/2lOsX72/RWfKKnVTeJCOnirVsVOl8jJJyfG9NeUn18vs4210mx5rytIcfbytUCNu6qQ3k6KNbgcAAAAwhMViUXBw8BXXoYZesszOzla/fv3si39JGjFihEpLS7V161Z7TXx8vH3xb6spLCzU/v377TUJCQkO5x4xYoS2bNmi8vLyems2btwoqXq3wtatW2vVJCQk2GvqMmfOHAUHB9sf3bp1a+SnAGc25LpQ/eWxW9W2la92Flp07FSpendorb8/OUS/HdmHxb/Bfn3HdTKZpP/beUR5hy1GtwMAAAA4NUMDgKKiInXq1MnhWEhIiPz8/FRUVHTZGtvzK9VUVFSouLi43hrbOYqLi1VZWVlvTV2mT5+ukpIS+6OgoKBBvztcx83d2upvj8cotnd7PTG0tz6dHKebu7U1ui1Iuq5joO7u31mSlPYF3wgAAAAA1KfRAcBzzz0nk8lU72PLli0NPp/JVPs70q1Wq8PxS2tsdy00Rc2lxxpSU5PZbFZQUJDDA+7n+k6B+utjt2rayD7y9+WqvzP59R3VX924cnuR/lN0yuBuAAAAAOfl09gXTJo0SePGjau3pmfPng06V1hYmMMQPkk6ceKEysvL7Vfiw8LCal2BP3r0qCRdscbHx0ft27evt8Z2jtDQUHl7e9dbA8D53BgWqFH9w7Rye5HSvtit18cPMrolAAAAwCk1egdAaGio+vTpU+/D39+/QeeKiYnRjh07dPjwYfuxzMxMmc1mRUVF2WvWrVvn8NWAmZmZCg8PtwcNMTExysrKcjh3ZmamoqOj5evrW29NbGysJMnPz09RUVG1arKysuw1AJzTpGHVuwA+3X5Yu4+wCwAAAACoS7POAMjPz9e2bduUn5+vyspKbdu2Tdu2bdPp06clVQ/Yi4yMVFJSknJycvT555/rqaee0mOPPWbfSj9+/HiZzWZNmDBBO3bs0PLlyzV79mylpqbat+YnJyfrwIEDSk1NVV5enhYvXqxFixbpqaeesvcyZcoUZWZmau7cudq1a5fmzp2r1atXKyUlxV6Tmpqqt99+W4sXL1ZeXp6mTp2q/Px8JScnN+fHBOAaRYYHKSGyk6xW6fU1e4xuBwAAAHBKzfo1gBMmTNB7771X6/iaNWs0dOhQSdUhwZNPPqkvvvhCAQEBGj9+vObNm+cw9X/79u2aOHGivv76a4WEhCg5OVnPPvusw735a9eu1dSpU7Vz506Fh4dr2rRptRbuH330kZ5++mnt3btXvXv31gsvvKCf/exnDjULFizQSy+9pMOHD6tfv3565ZVXdPvttzf4d27o1y8AaFo7DpXonrQN8jJJWanx6t2hjdEtAQAAAC2ioevQZg0APBEBAGCcR9/brNV5R/WzgV00P/Fmo9sBAAAAWkRD16GGfg0gADSlyXdWzwL4x7ZD2ld8xuBuAAAAAOdCAADAbQzo2lbDbuygKqv0BrMAAAAAAAcEAADcim0XwPKcQzpwnF0AAAAAgA0BAAC3MrB7iG6/oYMqq6xasOZ7o9sBAAAAnAYBAAC3M+XCLoCMfx9UwQ9nDe4GAAAAcA4EAADcTlSPEN12Xagqqqxa8CW7AAAAAACJAACAm5ryk+pdAB9tLdChk+cM7gYAAAAwHgEAALd0S892iunVXuWVVqV/yTcCAAAAAAQAANyWbRfAh5sP6nAJuwAAAADg2QgAALitW3u1148j2qmsskoLmQUAAAAAD0cAAMCtpVz4RoAlmwt0xHLe4G4AAAAA4xAAAHBrMb3bK7pHiMoqqrRwLbsAAAAA4LkIAAC4NZPJZJ8F8NdN+Tp6il0AAAAA8EwEAADc3m3XhWpg97YqrajSW2v3Gt0OAAAAYAgCAABuz2QyacqFWQAfbDqg4tOlBncEAAAAtDwCAAAeIf6GDvpR12CdL6/Sn9axCwAAAACehwAAgEeoOQvg/ewDOs4uAAAAAHgYAgAAHmPYjR3Vv0uwzpVX6u0N+4xuBwAAAGhRBAAAPIbJZNLkC7MA3t+4XyfOlBncEQAAANByCAAAeJSf9O2oyM5BOlNWqUXsAgAAAIAHIQAA4FFq7gJ4d+N+lZwtN7gjAAAAoGUQAADwOAmRndQnLFCnSyu06F/sAgAAAIBnIAAA4HG8vC7uAnjnX/tUco5dAAAAAHB/BAAAPNLIm8J0Q6c2OnW+Qu/+a7/R7QAAAADNjgAAgEfy8jLp13dU7wJYtGGvTp1nFwAAAADcGwEAAI81qn9n9e7QWpbzFXpv436j2wEAAACaFQEAAI/lXWMWwNsb9ul0aYXBHQEAAADNhwAAgEe7Z0C4eoW21smz5Xo/e7/R7QAAAADNhgAAgEfz9jJp0h3XSZLeXr9PZ9gFAAAAADdFAADA4/3Xj8LVo30r/XCmTB98dcDodgAAAIBmQQAAwOP5eHtp4rDqXQBvrdurc2WVBncEAAAAND0CAACQdO/ALurWLkDHz5TpL5vYBQAAAAD3QwAAAJJ8vb00cWj1LoCFa/fqfDm7AAAAAOBeCAAA4IKfDeqqLm0DVHy6VH/dlG90OwAAAECTIgAAgAv8fLz05LDekqSFa79nFwAAAADcCgEAANTw86iuCg/219FTpVq2ucDodgAAAIAmQwAAADWYfbz1xNDqXQDpX36v0gp2AQAAAMA9EAAAwCXG3tJNYUH+KrKc14dbDhrdDgAAANAkCAAA4BJmH28lx/eSJKWv2aOyiiqDOwIAAACuHQEAANRh3I+7q2OgWYUl55Xxb3YBAAAAwPURAABAHfx9vfV4fPUsgDfW7FF5JbsAAAAA4NoIAADgMsb/uLtC25h18MQ5Lf/3IaPbAQAAAK4JAQAAXEaAn7cev716FsDr7AIAAACAiyMAAIB6/OLW7mrf2k/5P5zVx9sKjW4HAAAAuGoEAABQj1Z+PnrMtgvgi92qYBcAAAAAXBQBAABcQdKtPRTSylf7j5/VJ9+wCwAAAACuiQAAAK6gtdlHj8ZV7wJI+2KPKqusBncEAAAANB4BAAA0wEOxPdW2la/2Hjujf7ILAAAAAC6IAAAAGqCN2UePDImQxC4AAAAAuCYCAABooIeG9FSQv4/2HD2tVTsOG90OAAAA0CgEAADQQEH+vnr4tupdAK99vltV7AIAAACACyEAAIBG+O8hEQo0++i7I6f1fzuLjG4HAAAAaDACAABohOAAX/33kJ6SpFfZBQAAAAAXQgAAAI308G0RamP20a6iU8rKO2J0OwAAAECDEAAAQCO1beWnh2J7SKqeBWC1sgsAAAAAzo8AAACuwqO39VIrP2/tLLTo87yjRrcDAAAAXBEBAABchZDWfvplTE9J1bMA2AUAAAAAZ0cAAABX6bG4CAX4emv7oRJ9+Z9jRrcDAAAA1IsAAACuUvs2ZiXFVM8C+CO7AAAAAODkCAAA4Bo8FtdL/r5eyi04qXW7i41uBwAAALgsAgAAuAYdAs36xeDqXQCvrv6OXQAAAABwWgQAAHCNHr+9l8w+Xvp3/kn9a89xo9sBAAAA6tSsAcALL7yg2NhYtWrVSm3btq3189zcXD3wwAPq1q2bAgIC1LdvX7366qu16rZv3674+HgFBASoS5cumjVrVq2rbGvXrlVUVJT8/f3Vq1cvLVy4sNZ5MjIyFBkZKbPZrMjISC1fvrxWzYIFCxQRESF/f39FRUVp/fr1V/8BAPAIHYP89cCPu0uSXv2cXQAAAABwTs0aAJSVlen+++/XE088UefPt27dqg4dOuiDDz7Qzp07NXPmTE2fPl2vv/66vcZisWj48OEKDw/X5s2blZaWpnnz5mn+/Pn2mn379mnUqFGKi4tTTk6OZsyYocmTJysjI8Nek52drcTERCUlJSk3N1dJSUkaO3asNm3aZK9ZtmyZUlJSNHPmTOXk5CguLk533XWX8vPzm+HTAeBOkuN7y8/bS5v3n1D2XnYBAAAAwPmYrC1wqerdd99VSkqKTp48ecXaiRMnKi8vT1988YUkKT09XdOnT9eRI0dkNpslSS+++KLS0tJ08OBBmUwmTZs2TStWrFBeXp79PMnJycrNzVV2drYkKTExURaLRatWrbLXjBw5UiEhIVqyZIkkafDgwRo0aJDS09PtNX379tWYMWM0Z86cBv2uFotFwcHBKikpUVBQUINeA8A9PPvxDr2ffUCDI9pp2eMxRrcDAAAAD9HQdajTzQAoKSlRu3bt7M+zs7MVHx9vX/xL0ogRI1RYWKj9+/fbaxISEhzOM2LECG3ZskXl5eX11mzcuFFS9W6FrVu31qpJSEiw19SltLRUFovF4QHAMyXH95avt0mb9v2gr9gFAAAAACfjVAFAdna2PvzwQz3++OP2Y0VFRerUqZNDne15UVFRvTUVFRUqLi6ut8Z2juLiYlVWVtZbU5c5c+YoODjY/ujWrVtjfmUAbiS8bYDGRlf/N+C1z3cb3A0AAADgqNEBwHPPPSeTyVTvY8uWLY1uZOfOnRo9erSeffZZDR8+3OFnJpPJ4bntroWax6+25tJjDampafr06SopKbE/CgoKLlsLwP09MbR6F8DG749r8/4fjG4HAAAAsPNp7AsmTZqkcePG1VvTs2fPRp3z22+/1R133KHHHntMTz/9tMPPwsLCal2BP3r0qKSLOwEuV+Pj46P27dvXW2M7R2hoqLy9veutqYvZbHa4PQGAZ+sa0ko/j+qqJV8X6LXPd+vPjww2uiUAAABA0lXsAAgNDVWfPn3qffj7+zf4fDt37tSwYcP00EMP6YUXXqj185iYGK1bt05lZWX2Y5mZmQoPD7cHDTExMcrKynJ4XWZmpqKjo+Xr61tvTWxsrCTJz89PUVFRtWqysrLsNQDQEE8OvU4+Xiat312srQdOGN0OAAAAIKmZZwDk5+dr27Ztys/PV2VlpbZt26Zt27bp9OnTki4u/ocPH67U1FQVFRWpqKhIx44ds59j/PjxMpvNmjBhgnbs2KHly5dr9uzZSk1NtW/NT05O1oEDB5Samqq8vDwtXrxYixYt0lNPPWU/z5QpU5SZmam5c+dq165dmjt3rlavXq2UlBR7TWpqqt5++20tXrxYeXl5mjp1qvLz85WcnNycHxMAN9OtXSv9bFAXScwCAAAAgPNo1q8BnDBhgt57771ax9esWaOhQ4fqueee0/PPP1/r5z169LBP+Jek7du3a+LEifr6668VEhKi5ORkPfvssw735q9du1ZTp07Vzp07FR4ermnTptVauH/00Ud6+umntXfvXvXu3VsvvPCCfvaznznULFiwQC+99JIOHz6sfv366ZVXXtHtt9/e4N+ZrwEEIEkHjp/RHX9Yq8oqq/4xcYhu7tbW6JYAAADgphq6Dm3WAMATEQAAsPnNh7nK+PdB3dGnoxZPuMXodgAAAOCmGroOdaqvAQQAdzLpjuvkZZK+2HVU2w+WGN0OAAAAPBwBAAA0k4jQ1hp9c/UsgFeZBQAAAACDEQAAQDOaOOw6mUzS6rwj2nGIXQAAAAAwDgEAADSj6zq20U8HhEuSXv9ij8HdAAAAwJMRAABAM/v1HdW7AD7bWaS8wxaj2wEAAICHIgAAgGZ2fadAjerfWRK7AAAAAGAcAgAAaAG/vuM6SdLKHYf13ZFTBncDAAAAT0QAAAAtoE9YkO7qFyarVUpjFwAAAAAMQAAAAC3k13dcL0n65zeF2nOUXQAAAABoWQQAANBCIsODlBDZSVYrswAAAADQ8ggAAKAFTb6zehfAitxC7T122uBuAAAA4EkIAACgBfXrEqyf9O2oKqv0+hp2AQAAAKDlEAAAQAuz7QL4eFuh9hefMbgbAAAAeAoCAABoYQO6ttWwGzuossqqN9gFAAAAgBZCAAAABrDtAvh7ziHlHz9rcDcAAADwBAQAAGCAgd1DdPsN1bsAFnzJLgAAAAA0PwIAADDIlDuvkyR9tPWgCn5gFwAAAACaFwEAABgkqkc73XZdqCqqrEpf+73R7QAAAMDNEQAAgIFsswD+tqVAh06eM7gbAAAAuDMCAAAw0I8j2unWXu1UXmnVwi/ZBQAAAIDmQwAAAAabcucNkqRlmwt0uIRdAAAAAGgeBAAAYLBbe7XTj3u2U1llld5cu9fodgAAAOCmCAAAwGAmk0lTflI9C+CvX+friOW8wR0BAADAHREAAIATiO3dXlE9QlRWUaWPth40uh0AAAC4IQIAAHACJpNJo/p3liRtKzhpbDMAAABwSwQAAOAk+ncJliTtOFRicCcAAABwRwQAAOAkbgoPkskkHS45r2OnSo1uBwAAAG6GAAAAnERrs496hbaWxC4AAAAAND0CAABwIrbbALYTAAAAAKCJEQAAgBPpRwAAAACAZkIAAABOZEDXtpK4BQAAAABNjwAAAJxIzUGAxacZBAgAAICmQwAAAE6k5iBAbgMAAABAUyIAAAAnYx8EeJAAAAAAAE2HAAAAnAyDAAEAANAcCAAAwMnYdgAwCBAAAABNiQAAAJzMTV2CGQQIAACAJkcAAABOpg2DAAEAANAMCAAAwAnZbwNgECAAAACaCAEAADgh2yDAb9gBAAAAgCZCAAAATohBgAAAAGhqBAAA4IQYBAgAAICmRgAAAE6ojdlHEQwCBAAAQBMiAAAAJ8UgQAAAADQlAgAAcFK2AIAdAAAAAGgKBAAA4KQYBAgAAICmRAAAAE7KNgiwkEGAAAAAaAIEAADgpBgECAAAgKZEAAAAToxBgAAAAGgqBAAA4MQYBAgAAICmQgAAAE6MQYAAAABoKgQAAODEag4CPM4gQAAAAFwDAgAAcGIMAgQAAEBTIQAAACdnnwPAIEAAAABcAwIAAHByDAIEAABAUyAAAAAn149BgAAAAGgCBAAA4ORuCg+SxCBAAAAAXBsCAABwcoH+vurVgUGAAAAAuDYEAADgAvpzGwAAAACuEQEAALgAWwDwDd8EAAAAgKtEAAAALoBBgAAAALhWBAAA4AIYBAgAAIBr1awBwAsvvKDY2Fi1atVKbdu2rbf2+PHj6tq1q0wmk06ePOnws+3btys+Pl4BAQHq0qWLZs2aJavV6lCzdu1aRUVFyd/fX7169dLChQtrvUdGRoYiIyNlNpsVGRmp5cuX16pZsGCBIiIi5O/vr6ioKK1fv77RvzcANLVAf1/1CmUQIAAAAK5eswYAZWVluv/++/XEE09csfaRRx7RgAEDah23WCwaPny4wsPDtXnzZqWlpWnevHmaP3++vWbfvn0aNWqU4uLilJOToxkzZmjy5MnKyMiw12RnZysxMVFJSUnKzc1VUlKSxo4dq02bNtlrli1bppSUFM2cOVM5OTmKi4vTXXfdpfz8/Gv8JADg2nEbAAAAAK6FyXrppfRm8O677yolJaXWlX2b9PR0LVu2TM8++6zuvPNOnThxwr5jID09XdOnT9eRI0dkNpslSS+++KLS0tJ08OBBmUwmTZs2TStWrFBeXp79nMnJycrNzVV2drYkKTExURaLRatWrbLXjBw5UiEhIVqyZIkkafDgwRo0aJDS09PtNX379tWYMWM0Z86cBv2uFotFwcHBKikpUVBQUIM/IwC4krfX79X/fpqnETd10ptJ0Ua3AwAAACfR0HWo4TMAvv32W82aNUvvv/++vLxqt5Odna34+Hj74l+SRowYocLCQu3fv99ek5CQ4PC6ESNGaMuWLSovL6+3ZuPGjZKqdyts3bq1Vk1CQoK9pi6lpaWyWCwODwBoDhd3APDfGQAAADSeoQFAaWmpHnjgAb388svq3r17nTVFRUXq1KmTwzHb86KionprKioqVFxcXG+N7RzFxcWqrKyst6Yuc+bMUXBwsP3RrVu3K/3aAHBVbIMAD508xyBAAAAANFqjA4DnnntOJpOp3seWLVsadK7p06erb9++evDBB+utM5lMDs9tdy3UPH61NZcea0jNpb9DSUmJ/VFQUFDv7wIAV4tBgAAAALgWPo19waRJkzRu3Lh6a3r27Nmgc33xxRfavn27PvroI0kXF+2hoaGaOXOmnn/+eYWFhdW6An/06FFJF3cCXK7Gx8dH7du3r7fGdo7Q0FB5e3vXW1MXs9nscHsCADSnfl2Ctbf4jHYcKtHQGzsa3Q4AAABcSKMDgNDQUIWGhjbJm2dkZOjcuXP255s3b9bDDz+s9evXq3fv3pKkmJgYzZgxQ2VlZfLz85MkZWZmKjw83B40xMTE6JNPPnE4d2ZmpqKjo+Xr62uvycrK0tSpUx1qYmNjJUl+fn6KiopSVlaW7r33XntNVlaWRo8e3SS/LwBcq/5dgrUit5AdAAAAAGi0RgcAjZGfn68ffvhB+fn5qqys1LZt2yRJ1113ndq0aWNf5NvY7tfv27ev/VsAxo8fr+eff14TJkzQjBkztHv3bs2ePVvPPvusfWt+cnKyXn/9daWmpuqxxx5Tdna2Fi1aZJ/uL0lTpkzR7bffrrlz52r06NH6+OOPtXr1am3YsMFek5qaqqSkJEVHRysmJkZvvfWW8vPzlZyc3IyfEgA0XP+uDAIEAADA1WnWAODZZ5/Ve++9Z38+cOBASdKaNWs0dOjQBp0jODhYWVlZmjhxoqKjoxUSEqLU1FSlpqbaayIiIrRy5UpNnTpVb7zxhsLDw/Xaa6/pvvvus9fExsZq6dKlevrpp/XMM8+od+/eWrZsmQYPHmyvSUxM1PHjxzVr1iwdPnxY/fr108qVK9WjR49r/CQAoGnUHAT4w5kytWvtZ3BHAAAAcBUmq+3GezSJhn7/IgBcrTvmfam9xWf03sM/VvwNHYxuBwAAAAZr6DrU0K8BBAA0Xr8u1bcBbD940thGAAAA4FIIAADAxfS3BQAMAgQAAEAjEAAAgIux7QBgECAAAAAagwAAAFzMTV0cBwECAAAADUEAAAAuJsjfV71CW0viNgAAAAA0HAEAALigi7cBEAAAAACgYQgAAMAF2QYBfsM3AQAAAKCBCAAAwAUxCBAAAACNRQAAAC6IQYAAAABoLAIAAHBBQf6+imAQIAAAABqBAAAAXBSDAAEAANAYBAAA4KIGXAgAth8kAAAAAMCVEQAAgIuy7QDgFgAAAAA0BAEAALgoBgECAACgMQgAAMBFMQgQAAAAjUEAAAAujEGAAAAAaCgCAABwYf0v3AbAIEAAAABcCQEAALiw/l3aSuIWAAAAAFwZAQAAuDAGAQIAAKChCAAAwIUxCBAAAAANRQAAAC6OQYAAAABoCAIAAHBxDAIEAABAQxAAAICLs+0A4BYAAAAA1IcAAABcnC0AOHTynE4wCBAAAACXQQAAAC6OQYAAAABoCAIAAHAD3AYAAACAKyEAAAA3wCBAAAAAXAkBAAC4AXYAAAAA4EoIAADADTAIEAAAAFdCAAAAbiDI31c927eSxC4AAAAA1I0AAADcBLcBAAAAoD4EAADgJgZ0rQ4AdhAAAAAAoA4EAADgJtgBAAAAgPoQAACAm7AFAAdPMAgQAAAAtREAAICbYBAgAAAA6kMAAABuhNsAAAAAcDkEAADgRvp3YRAgAAAA6kYAAABupH9XdgAAAACgbgQAAOBGGAQIAACAyyEAAAA3wiBAAAAAXA4BAAC4GQYBAgAAoC4EAADgZhgECAAAgLoQAACAm+nPDgAAAADUgQAAANzMTQwCBAAAQB0IAADAzQQHXBwEuKOQXQAAAACoRgAAAG6IQYAAAAC4FAEAALgh+xyAgwQAAAAAqEYAAABuiEGAAAAAuBQBAAC4IQYBAgAA4FIEAADghoIDfNWDQYAAAACogQAAANwUgwABAABQEwEAALipARcCgB0EAAAAABABAAC4LQYBAgAAoCYCAABwU7ZBgAU/MAgQAAAABAAA4LYYBAgAAICaCAAAwI0xCBAAAAA2BAAA4Mb6MwgQAAAAFxAAAIAbYxAgAAAAbAgAAMCN9Qu/OAjw5FkGAQIAAHgyAgAAcGPBrS4OAmQXAAAAgGcjAAAAN8cgQAAAAEgEAADg9hgECAAAAKmZA4AXXnhBsbGxatWqldq2bXvZunfffVcDBgyQv7+/wsLCNGnSJIefb9++XfHx8QoICFCXLl00a9YsWa1Wh5q1a9cqKipK/v7+6tWrlxYuXFjrfTIyMhQZGSmz2azIyEgtX768Vs2CBQsUEREhf39/RUVFaf369Vf3ywOAk2AQIAAAAKRmDgDKysp0//3364knnrhszfz58zVz5kz97ne/086dO/X5559rxIgR9p9bLBYNHz5c4eHh2rx5s9LS0jRv3jzNnz/fXrNv3z6NGjVKcXFxysnJ0YwZMzR58mRlZGTYa7Kzs5WYmKikpCTl5uYqKSlJY8eO1aZNm+w1y5YtU0pKimbOnKmcnBzFxcXprrvuUn5+fhN/MgDQchgECAAAAEkyWS+9lN4M3n33XaWkpOjkyZMOx0+cOKEuXbrok08+0Z133lnna9PT0zV9+nQdOXJEZrNZkvTiiy8qLS1NBw8elMlk0rRp07RixQrl5eXZX5ecnKzc3FxlZ2dLkhITE2WxWLRq1Sp7zciRIxUSEqIlS5ZIkgYPHqxBgwYpPT3dXtO3b1+NGTNGc+bMadDvarFYFBwcrJKSEgUFBTXoNQDQ3OJfXqMDx8/qg0cG67brQ41uBwAAAE2ooetQQ2cAZGVlqaqqSocOHVLfvn3VtWtXjR07VgUFBfaa7OxsxcfH2xf/kjRixAgVFhZq//799pqEhASHc48YMUJbtmxReXl5vTUbN26UVL1bYevWrbVqEhIS7DV1KS0tlcVicXgAgLNhECAAAAAMDQD27t2rqqoqzZ49W3/84x/10Ucf6YcfftDw4cNVVla9TbWoqEidOnVyeJ3teVFRUb01FRUVKi4urrfGdo7i4mJVVlbWW1OXOXPmKDg42P7o1q1bYz8GAGh2F+cAnDS2EQAAABim0QHAc889J5PJVO9jy5YtDTpXVVWVysvL9dprr2nEiBG69dZbtWTJEu3evVtr1qyx15lMJofX2e5aqHn8amsuPdaQmpqmT5+ukpIS+6Pm7gUAcBYMAgQAAIBPY18wadIkjRs3rt6anj17NuhcnTt3liRFRkbaj3Xo0EGhoaH2wXthYWG1rsAfPXpU0sWdAJer8fHxUfv27eutsZ0jNDRU3t7e9dbUxWw2O9yeAADO6NJBgG1b+RncEQAAAFpao3cAhIaGqk+fPvU+/P39G3SuIUOGSJL+85//2I/98MMPKi4uVo8ePSRJMTExWrdunf2WAEnKzMxUeHi4PWiIiYlRVlaWw7kzMzMVHR0tX1/femtiY2MlSX5+foqKiqpVk5WVZa8BAFcV3MpX3du1kiTtOMSsEgAAAE/UrDMA8vPztW3bNuXn56uyslLbtm3Ttm3bdPr0aUnSDTfcoNGjR2vKlCnauHGjduzYoYceekh9+vTRsGHDJEnjx4+X2WzWhAkTtGPHDi1fvlyzZ89WamqqfWt+cnKyDhw4oNTUVOXl5Wnx4sVatGiRnnrqKXsvU6ZMUWZmpubOnatdu3Zp7ty5Wr16tVJSUuw1qampevvtt7V48WLl5eVp6tSpys/PV3JycnN+TADQIrgNAAAAwLM1+haAxnj22Wf13nvv2Z8PHDhQkrRmzRoNHTpUkvT+++9r6tSpuvvuu+Xl5aX4+Hh99tln9iv3wcHBysrK0sSJExUdHa2QkBClpqYqNTXVft6IiAitXLlSU6dO1RtvvKHw8HC99tpruu++++w1sbGxWrp0qZ5++mk988wz6t27t5YtW6bBgwfbaxITE3X8+HHNmjVLhw8fVr9+/bRy5Ur7bgQAcGX9uwbr0+2HtYMAAAAAwCOZrLZpeWgSDf3+RQBoaf/aU6xfvL1J3du10rrfDjO6HQAAADSRhq5DDf0aQABAy7ENAsz/4axOni27QjUAAADcDQEAAHgIBgECAAB4NgIAAPAgDAIEAADwXAQAAOBB+l0IABgECAAA4HkIAADAg7ADAAAAwHMRAACAB7EFAPk/nFXJ2XKDuwEAAEBLIgAAAA9ScxAguwAAAAA8CwEAAHgYbgMAAADwTAQAAOBhGAQIAADgmQgAAMDDsAMAAADAMxEAAICH6dclSBKDAAEAADwNAQAAeJi2rfzsgwB3FLILAAAAwFMQAACAB+I2AAAAAM9DAAAAHsg2CHD7QQIAAAAAT0EAAAAeiB0AAAAAnocAAAA8EIMAAQAAPA8BAAB4oLat/NStXYAkBgECAAB4CgIAAPBQ3AYAAADgWQgAAMBD9e/SVhIBAAAAgKcgAAAAD2XbAbCDAAAAAMAjEAAAgIeyDQI8cJxBgAAAAJ6AAAAAPBSDAAEAADwLAQAAeDAGAQIAAHgOAgAA8GD9CAAAAAA8BgEAAHgwBgECAAB4DgIAAPBgtgCAQYAAAADujwAAADwYgwABAAA8BwEAAHg4BgECAAB4BgIAAPBwDAIEAADwDAQAAODhGAQIAADgGQgAAMDD9QuvMQjwHIMAAQAA3BUBAAB4uJDWFwcB7mQXAAAAgNsiAAAAMAgQAADAAxAAAADsgwC/IQAAAABwWwQAAAAGAQIAAHgAAgAAAIMAAQAAPAABAABAIa391DWEQYAAAADujAAAACCJQYAAAADujgAAACBJ6t+VAAAAAMCdEQAAACSxAwAAAMDdEQAAACQxCBAAAMDdEQAAACQxCBAAAMDdEQAAAOy4DQAAAMB9EQAAAOz6EQAAAAC4LQIAAICdbQfADgIAAAAAt0MAAACwswUA+xkECAAA4HYIAAAAdgwCBAAAcF8EAAAABwwCBAAAcE8EAAAABwwCBAAAcE8EAAAABwwCBAAAcE8EAAAABzUHAVrOMwgQAADAXRAAAAAc1BwEyC4AAAAA90EAAACohdsAAAAA3A8BAACgFtsgwG8OEgAAAAC4CwIAAEAt7AAAAABwPwQAAIBaGAQIAADgfggAAAC1hLT2U5e2DAIEAABwJwQAAIA6cRsAAACAeyEAAADUqX/X6gBg+yGLwZ0AAACgKRAAAADqZNsBsP3gSWMbAQAAQJMgAAAA1IlBgAAAAO6FAAAAUCcGAQIAALiXZg0AXnjhBcXGxqpVq1Zq27ZtnTWbN2/WnXfeqbZt2yokJEQJCQnatm2bQ8327dsVHx+vgIAAdenSRbNmzZLVanWoWbt2raKiouTv769evXpp4cKFtd4rIyNDkZGRMpvNioyM1PLly2vVLFiwQBEREfL391dUVJTWr19/1b8/ALg6BgECAAC4j2YNAMrKynT//ffriSeeqPPnp06d0ogRI9S9e3dt2rRJGzZsUFBQkEaMGKHy8urtphaLRcOHD1d4eLg2b96stLQ0zZs3T/Pnz7efZ9++fRo1apTi4uKUk5OjGTNmaPLkycrIyLDXZGdnKzExUUlJScrNzVVSUpLGjh2rTZs22WuWLVumlJQUzZw5Uzk5OYqLi9Ndd92l/Pz8ZvqEAMC5MQgQAADAfZisl15KbwbvvvuuUlJSdPLkSYfjW7Zs0S233KL8/Hx169ZNUvXV/gEDBmjPnj3q3bu30tPTNX36dB05ckRms1mS9OKLLyotLU0HDx6UyWTStGnTtGLFCuXl5dnPnZycrNzcXGVnZ0uSEhMTZbFYtGrVKnvNyJEjFRISoiVLlkiSBg8erEGDBik9Pd1e07dvX40ZM0Zz5sxp0O9qsVgUHByskpISBQUFNf7DAgAnsva7Y3po8deKCG2tNU8NNbodAAAA1KGh61BDZwDceOONCg0N1aJFi1RWVqZz585p0aJFuummm9SjRw9J1Vfu4+Pj7Yt/SRoxYoQKCwu1f/9+e01CQoLDuUeMGKEtW7bYdxJcrmbjxo2SqncrbN26tVZNQkKCvaYupaWlslgsDg8AcBe2WwD2FZ9hECAAAICLMzQACAwM1JdffqkPPvhAAQEBatOmjf7v//5PK1eulI+PjySpqKhInTp1cnid7XlRUVG9NRUVFSouLq63xnaO4uJiVVZW1ltTlzlz5ig4ONj+sO1kAAB30I5BgAAAAG6j0QHAc889J5PJVO9jy5YtDTrXuXPn9PDDD2vIkCH66quv9K9//Us33XSTRo0apXPnztnrTCaTw+tsdy3UPH61NZcea0hNTdOnT1dJSYn9UVBQcNlaAHBFDAIEAABwDz6NfcGkSZM0bty4emt69uzZoHP99a9/1f79+5WdnS0vLy/7sZCQEH388ccaN26cwsLCal2BP3r0qKSLOwEuV+Pj46P27dvXW2M7R2hoqLy9veutqYvZbHa4PQEA3E3/rsH6bGcRgwABAABcXKMDgNDQUIWGhjbJm589e1ZeXl4OV9htz6uqqiRJMTExmjFjhsrKyuTn5ydJyszMVHh4uD1oiImJ0SeffOJw7szMTEVHR8vX19dek5WVpalTpzrUxMbGSpL8/PwUFRWlrKws3XvvvfaarKwsjR49ukl+XwBwRf3YAQAAAOAWmnUGQH5+vrZt26b8/HxVVlZq27Zt2rZtm06fPi1JGj58uE6cOKGJEycqLy9PO3fu1H//93/Lx8dHw4YNkySNHz9eZrNZEyZM0I4dO7R8+XLNnj1bqamp9uAgOTlZBw4cUGpqqvLy8rR48WItWrRITz31lL2XKVOmKDMzU3PnztWuXbs0d+5crV69WikpKfaa1NRUvf3221q8eLHy8vI0depU5efnKzk5uTk/JgBwagwCBAAAcA+N3gHQGM8++6zee+89+/OBAwdKktasWaOhQ4eqT58++uSTT/T8888rJiZGXl5eGjhwoD777DN17txZkhQcHKysrCxNnDhR0dHRCgkJUWpqqlJTU+3njYiI0MqVKzV16lS98cYbCg8P12uvvab77rvPXhMbG6ulS5fq6aef1jPPPKPevXtr2bJlGjx4sL0mMTFRx48f16xZs3T48GH169dPK1eutH8jAQB4ItsgwEMnz2nnIYtierc3uiUAAABcBZPVNi0PTaKh378IAK4k+c9b9dnOIs0c1VeP3d7L6HYAAABQQ0PXoYZ+DSAAwDX071p9G8A3zAEAAABwWQQAAIArYhAgAACA6yMAAABcEYMAAQAAXB8BAADgimyDACVp5yGLwd0AAADgahAAAAAapF+X6oEy3AYAAADgmggAAAANMqBrW0nSdgIAAAAAl0QAAABoENsgQAIAAAAA10QAAABoEAYBAgAAuDYCAABAgzAIEAAAwLURAAAAGoxBgAAAAK6LAAAA0GD9mQMAAADgsggAAAANZhsEyA4AAAAA10MAAABoMNsOgL3FZ3SKQYAAAAAuhQAAANBg7duY7YMAdzAIEAAAwKUQAAAAGoVBgAAAAK6JAAAA0CgMAgQAAHBNBAAAgEZhECAAAIBrIgAAADQKgwABAABcEwEAAKBRag4C3FnIIEAAAABXQQAAAGg0BgECAAC4HgIAAECj2W4D+OYgAQAAAICrIAAAADQagwABAABcDwEAAKDRGAQIAADgeggAAACN1r6NWeHB/pIYBAgAAOAqCAAAAFeF2wAAAABcCwEAAOCqDOhaHQBsJwAAAABwCQQAAICrYtsBsJ1vAgAAAHAJBAAAgKvCIEAAAADXQgAAALgqDAIEAABwLQQAAICrxiBAAAAA10EAAAC4arbbABgECAAA4PwIAAAAV60f3wQAAADgMggAAABXzbYDYB+DAAEAAJweAQAA4KqFXhgEaLUyCBAAAMDZEQAAAK4JgwABAABcAwEAAOCaMAgQAADANRAAAACuCYMAAQAAXAMBAADgmtQcBHi6tMLgbgAAAHA5BAAAgGsS2saszrZBgOwCAAAAcFoEAACAa8YcAAAAAOdHAAAAuGYEAAAAAM6PAAAAcM0YBAgAAOD8CAAAANeMQYAAAADOjwAAAHDNGAQIAADg/AgAAABNoh9zAAAAAJwaAQAAoEkMuBAA7CAAAAAAcEoEAACAJmEbBPgNAQAAAIBTIgAAADQJBgECAAA4NwIAAECTYBAgAACAcyMAAAA0GQYBAgAAOC8CAABAk+nPIEAAAACnRQAAAGgy/dkBAAAA4LQIAAAATcZ2C8BeBgECAAA4HQIAAECT6RDIIEAAAABnRQAAAGhSDAIEAABwTgQAAIAmxSBAAAAA50QAAABoUgwCBAAAcE4EAACAJsUgQAAAAOdEAAAAaFIdAs0KC6oeBPhtocXodgAAAHABAQAAoMn178ptAAAAAM6GAAAA0OTscwAOnjS2EQAAANgRAAAAmhyDAAEAAJxPswUA+/fv1yOPPKKIiAgFBASod+/e+v3vf6+ysjKHuvz8fP30pz9V69atFRoaqsmTJ9eq2b59u+Lj4xUQEKAuXbpo1qxZslqtDjVr165VVFSU/P391atXLy1cuLBWTxkZGYqMjJTZbFZkZKSWL19eq2bBggWKiIiQv7+/oqKitH79+ib4NADAszAIEAAAwPk0WwCwa9cuVVVV6c0339TOnTv1yiuvaOHChZoxY4a9prKyUnfffbfOnDmjDRs2aOnSpcrIyNBvfvMbe43FYtHw4cMVHh6uzZs3Ky0tTfPmzdP8+fPtNfv27dOoUaMUFxennJwczZgxQ5MnT1ZGRoa9Jjs7W4mJiUpKSlJubq6SkpI0duxYbdq0yV6zbNkypaSkaObMmcrJyVFcXJzuuusu5efnN9fHBABuiUGAAAAAzsdkvfRSejN6+eWXlZ6err1790qSVq1apXvuuUcFBQUKDw+XJC1dulQTJkzQ0aNHFRQUpPT0dE2fPl1HjhyR2WyWJL344otKS0vTwYMHZTKZNG3aNK1YsUJ5eXn290pOTlZubq6ys7MlSYmJibJYLFq1apW9ZuTIkQoJCdGSJUskSYMHD9agQYOUnp5ur+nbt6/GjBmjOXPmNOh3tFgsCg4OVklJiYKCgq7h0wIA1/boe1u0Ou+InrknUo/cFmF0OwAAAG6roevQFp0BUFJSonbt2tmfZ2dnq1+/fvbFvySNGDFCpaWl2rp1q70mPj7evvi31RQWFmr//v32moSEBIf3GjFihLZs2aLy8vJ6azZu3ChJKisr09atW2vVJCQk2GvqUlpaKovF4vAAAEgDLnwTwA7mAAAAADiFFgsAvv/+e6WlpSk5Odl+rKioSJ06dXKoCwkJkZ+fn4qKii5bY3t+pZqKigoVFxfXW2M7R3FxsSorK+utqcucOXMUHBxsf3Tr1q3+DwIAPIRtEOA3fBMAAACAU2h0APDcc8/JZDLV+9iyZYvDawoLCzVy5Ejdf//9evTRRx1+ZjKZar2H1Wp1OH5pje2uhaaoufRYQ2pqmj59ukpKSuyPgoKCy9YCgCdhECAAAIBz8WnsCyZNmqRx48bVW9OzZ0/73wsLCzVs2DDFxMTorbfecqgLCwtzGMInSSdOnFB5ebn9SnxYWFitK/BHjx6VpCvW+Pj4qH379vXW2M4RGhoqb2/vemvqYjabHW5PAABUsw0CLLKc17eFFv04ot2VXwQAAIBm0+gdAKGhoerTp0+9D39/f0nSoUOHNHToUA0aNEjvvPOOvLwc3y4mJkY7duzQ4cOH7ccyMzNlNpsVFRVlr1m3bp3DVwNmZmYqPDzcHjTExMQoKyvL4dyZmZmKjo6Wr69vvTWxsbGSJD8/P0VFRdWqycrKstcAABrHtgtgO3MAAAAADNdsMwAKCws1dOhQdevWTfPmzdOxY8dUVFTkcIU9ISFBkZGRSkpKUk5Ojj7//HM99dRTeuyxx+yTC8ePHy+z2awJEyZox44dWr58uWbPnq3U1FT71vzk5GQdOHBAqampysvL0+LFi7Vo0SI99dRT9veaMmWKMjMzNXfuXO3atUtz587V6tWrlZKSYq9JTU3V22+/rcWLFysvL09Tp05Vfn6+w9wCAEDD2eYAMAgQAADAeI2+BaChMjMztWfPHu3Zs0ddu3Z1+Jnt/nxvb299+umnevLJJzVkyBAFBARo/Pjxmjdvnr02ODhYWVlZmjhxoqKjoxUSEqLU1FSlpqbaayIiIrRy5UpNnTpVb7zxhsLDw/Xaa6/pvvvus9fExsZq6dKlevrpp/XMM8+od+/eWrZsmQYPHmyvSUxM1PHjxzVr1iwdPnxY/fr108qVK9WjR4/m+pgAwK3171od5rIDAAAAwHgmq201jibR0O9fBABPcPTUef34hc9lMkk7nhuh1uZmy50BAAA8VkPXoS32NYAAAM/TMdBfYUH+slqlnYUWo9sBAADwaAQAAIBmxSBAAAAA50AAAABoVgwCBAAAcA4EAACAZsUgQAAAAOdAAAAAaFa2WwC+P3ZaZ0orDO4GAADAcxEAAACaVcdAf3UKMstqlb49zCBAAAAAoxAAAACaXf8ubSVJ3xzkNgAAAACjEAAAAJodgwABAACMRwAAAGh2DAIEAAAwHgEAAKDZMQgQAADAeAQAAIBmxyBAAAAA4xEAAABahG0OwHYGAQIAABiCAAAA0CJs3wTAIEAAAABjEAAAAFqEbRDgNwQAAAAAhiAAAAC0CAYBAgAAGIsAAADQIhgECAAAYCwCAABAi2EQIAAAgHEIAAAALcZ2GwCDAAEAAFoeAQAAoMXYdwAQAAAAALQ4AgAAQIvpzyBAAAAAwxAAAABaTMeg6kGAVQwCBAAAaHEEAACAFsUgQAAAAGMQAAAAWhSDAAEAAIxBAAAAaFEMAgQAADAGAQAAoEXVHAR4toxBgAAAAC2FAAAA0KI6BvmrY+CFQYCFDAIEAABoKQQAAIAWN6Br9S6AbxgECAAA0GIIAAAALY5BgAAAAC2PAAAA0OIYBAgAANDyCAAAAC2OQYAAAAAtjwAAANDiGAQIAADQ8ggAAACG4DYAAACAlkUAAAAwRP+uBAAAAAAtiQAAAGAI+w4AvgoQAACgRRAAAAAMwSBAAACAlkUAAAAwBIMAAQAAWhYBAADAMAwCBAAAaDkEAAAAw/QjAAAAAGgxBAAAAMPYdgDsIAAAAABodgQAAADD2L4KcM9RBgECAAA0NwIAAIBhOjEIEAAAoMUQAAAADMUgQAAAgJZBAAAAMBSDAAEAAFoGAQAAwFAMAgQAAGgZBAAAAEMxCBAAAKBlEAAAAAzVKchfHS4MAsw7zCBAAACA5kIAAAAw3IALtwF8c5DbAAAAAJoLAQAAwHAMAgQAAGh+BAAAAMMxCBAAAKD5EQAAAAzHIEAAAIDmRwAAADAcgwABAACaHwEAAMAp2G4D2M4gQAAAgGZBAAAAcAr2AOAQOwAAAACaAwEAAMApXAwAThrbCAAAgJsiAAAAOAUGAQIAADQvAgAAgFNgECAAAEDzIgAAADgNBgECAAA0HwIAAIDT6McgQAAAgGZDAAAAcBq2HQA7DrEDAAAAoKkRAAAAnMaAC4MAdx89pXNllQZ3AwAA4F4IAAAATqPmIMBvD7MLAAAAoCkRAAAAnAqDAAEAAJpHswUA+/fv1yOPPKKIiAgFBASod+/e+v3vf6+ysjJ7TW5urh544AF169ZNAQEB6tu3r1599dVa59q+fbvi4+MVEBCgLl26aNasWbJarQ41a9euVVRUlPz9/dWrVy8tXLiw1nkyMjIUGRkps9msyMhILV++vFbNggULFBERIX9/f0VFRWn9+vVN8GkAABqKQYAAAADNo9kCgF27dqmqqkpvvvmmdu7cqVdeeUULFy7UjBkz7DVbt25Vhw4d9MEHH2jnzp2aOXOmpk+frtdff91eY7FYNHz4cIWHh2vz5s1KS0vTvHnzNH/+fHvNvn37NGrUKMXFxSknJ0czZszQ5MmTlZGRYa/Jzs5WYmKikpKSlJubq6SkJI0dO1abNm2y1yxbtkwpKSmaOXOmcnJyFBcXp7vuukv5+fnN9TEBAC7BIEAAAIDmYbJeeim9Gb388stKT0/X3r17L1szceJE5eXl6YsvvpAkpaena/r06Tpy5IjMZrMk6cUXX1RaWpoOHjwok8mkadOmacWKFcrLy7OfJzk5Wbm5ucrOzpYkJSYmymKxaNWqVfaakSNHKiQkREuWLJEkDR48WIMGDVJ6erq9pm/fvhozZozmzJnToN/RYrEoODhYJSUlCgoKauAnAwCwKSo5r1vnfC4vk7Tz+ZEK8PM2uiUAAACn1tB1aIvOACgpKVG7du0aVZOdna34+Hj74l+SRowYocLCQu3fv99ek5CQ4HCeESNGaMuWLSovL6+3ZuPGjZKksrIybd26tVZNQkKCvaYupaWlslgsDg8AwNXrFGRWaBvbIED+mwoAANBUWiwA+P7775WWlqbk5OTL1mRnZ+vDDz/U448/bj9WVFSkTp06OdTZnhcVFdVbU1FRoeLi4nprbOcoLi5WZWVlvTV1mTNnjoKDg+2Pbt26XbYWAHBlJpPJ/nWA2w+eNLYZAAAAN9LoAOC5556TyWSq97FlyxaH1xQWFmrkyJG6//779eijj9Z53p07d2r06NF69tlnNXz4cIefmUwmh+e2uxZqHr/amkuPNaSmpunTp6ukpMT+KCgouGwtAKBhGAQIAADQ9Hwa+4JJkyZp3Lhx9db07NnT/vfCwkINGzZMMTExeuutt+qs//bbb3XHHXfoscce09NPP+3ws7CwsFpX4I8ePSrp4k6Ay9X4+Pioffv29dbYzhEaGipvb+96a+piNpsdbk8AAFw7BgECAAA0vUYHAKGhoQoNDW1Q7aFDhzRs2DBFRUXpnXfekZdX7Q0HO3fu1B133KGHHnpIL7zwQq2fx8TEaMaMGSorK5Ofn58kKTMzU+Hh4fagISYmRp988onD6zIzMxUdHS1fX197TVZWlqZOnepQExsbK0ny8/NTVFSUsrKydO+999prsrKyNHr06Ab9vgCApmELAHYfPaVzZZUMAgQAAGgCzTYDoLCwUEOHDlW3bt00b948HTt2TEVFRQ5X2Hfu3Klhw4Zp+PDhSk1Ntf/82LFj9prx48fLbDZrwoQJ2rFjh5YvX67Zs2crNTXVvjU/OTlZBw4cUGpqqvLy8rR48WItWrRITz31lP08U6ZMUWZmpubOnatdu3Zp7ty5Wr16tVJSUuw1qampevvtt7V48WLl5eVp6tSpys/Pr3duAQCg6TEIEAAAoOk1egdAQ2VmZmrPnj3as2ePunbt6vAz2/35f/vb33Ts2DH95S9/0V/+8hf7z3v06GGf8B8cHKysrCxNnDhR0dHRCgkJUWpqqlJTU+31ERERWrlypaZOnao33nhD4eHheu2113TffffZa2JjY7V06VI9/fTTeuaZZ9S7d28tW7ZMgwcPttckJibq+PHjmjVrlg4fPqx+/fpp5cqV6tGjR3N8RACAyzCZTOrfJUhr/nNMOw6VKKpHiNEtAQAAuDyT1bYaR5No6PcvAgDqNz/zP3rtiz36eVRXzbv/R0a3AwAA4LQaug5tsa8BBACgMfp3bStJ2n6QQYAAAABNgQAAAOCULh0ECAAAgGtDAAAAcEoMAgQAAGhaBAAAAKdkGwQoSTsOcRsAAADAtSIAAAA4LdttANsJAAAAAK4ZAQAAwGn1uxAAsAMAAADg2hEAAACc1oAL3wSw++hpBgECAABcIwIAAIDTsg0CrKyyMggQAADgGhEAAACcFoMAAQAAmg4BAADAqTEIEAAAoGkQAAAAnBqDAAEAAJoGAQAAwKn171odAOw+elrnyxkECAAAcLUIAAAATi0syF+hbfwYBAgAAHCNCAAAAE6tehDghTkAB7kNAAAA4GoRAAAAnB6DAAEAAK4dAQAAwOkxCBAAALS0isoq7T5ySqfOlxvdSpPxMboBAACu5NJBgP6+3gZ3BAAA3MmJM2XKO2xRXtEp5R22aFeRRd8dOa2yiir96ZfRGh7ZyegWmwQBAADA6dkGARafLtO3hy0a1D3E6JYAAIALKq+s0r7iM9WL/cMXF/tHLKV11rfy89aJM2Ut3GXzIQAAADg9k8mkfl2C9eV/jmnHoRICAAAAcEXHT5dq14Ur+nmHT2lXkUW7j5xWWWVVnfXd27VS386B6hMWpL6dg9S3c6C6hbSSl5ephTtvPgQAAACX0P9CAMA3AQAAgJrKK6v0/bHT2nX4lPKKLiz2D1t09FTdV/Vb+3mrz4UFfvViP1A3hgWpjdn9l8fu/xsCANwC3wQAAACKT5dWb9uvsdjfc/SUyiutddb3bN/KfkW/T+dA9Q0LUteQALe6qt8YBAAAAJfAIEAAADxHWcWFq/pFF+/Vzzt8SsWn676qH2j2UZ8LV/T7dA5U385BurFToFp7wFX9xuDTAAC4BAYBAgDgno6dKrUP47Mt9r8/drrOq/omk9SzfWv79v0+YdWL/a4hATKZPPOqfmMQAAAAXAKDAAEAcG1lFVXac/S0w2J/V5FFxafrnrIf6O+jvjWu6PcJC9SNYYFq5ccy9mrxyQEAXAaDAAEAcH5Wq7X6qn5R9TC+6gX/Ke05eloVVXVf1Y9o39q+yLfdr9+lLVf1mxoBAADAZfRjECAAAE6ltKJSu4+c1i7bYr+oekDf8TN1X9UP8vdRn85Biqyx2L+hU6AC/Jjt0xIIAAAALmMAgwABADCE1WrV0Qv36tu27lffq39GlXVc1fcySRGhrWst9jsH+3NV30AEAAAAl8EgQAAAmt/58kr7vfq2xf6uolP64TJX9YMDfO1D+SIvbN+/viNX9Z0RAQAAwGUwCBAAgKZjtVp1xHLhqr5tKN9hi/YWX/6qfq8Obez36tsW+2FBXNV3FQQAAACXwiBAAACqWa1WVVmliqoqVVZZVVFlVWWlVZVWq8PziqoqVVmtOl1aqe+KTtnv088rsujk2fI6z922la/6hgXZB/JFdg7SdR3bcPudiyMAAAC4FAYBAgCsVqvKKqtUUXlhkVt1YZFbdcliuMaj+nn1a2ovkK2qsl5SU1WjrvJifZXV9rzK4T0qLnmfyio51FRUWVVV5dhv7f7qr6u45M+6rtI3lreXSb07tFafGov9vmFB6hRk5qq+GyIAAAC4lP5dGAQIAK6mtKJSZ0srdaasQmfLKnW6tKLG8wqdLq3U2dIKnSmz/VmhM6WVOnvhz+rn1a+1/VnX18nhIm8vk7y9TPKp8aefj5eu69jm4mI/LJCr+h6GAAAA4FI6B/urfWs/HT9TprzDFg1kDgAANKmyiqrqhfeFxfjpSxbe1c8vLtBP11io2153prTi4oK+rELllS2zWPf1rl7sepsuLHq9veyLXy+TST7eNRfFXvL2kry9vBwWybUXzl7yumQh7V2rrvocl6urWVP7tdV9eZmq38v23PH3qNGzySRvb8f38bq0X5O4eo86EQAAAFyKyWRS/64X5gAcKiEAAODRyiurLn8lvazikivttivvF39+pvTiwv5MWXVtWWVVs/Vr9vFSa7OPWpu91drPR638vKuf+/mole3YhT+rj3ur1YU/66rz8/a6sHCuXgQDqB8BAADA5TAIEICrslqtOlNWKcu5cp06X6FT56v/tJwvd9jyXn0lvfYC/Uyp45X2sormW6z7+XhVL8D9fNTGXGPh7ed9yfMLC3rzhQW9bfFu9nb8ma+3fLy9mq1fAFdGAAAAcDkMAgRgBKvVqrNllbKcv7h4t5yvcFjI2xf09gV+hUP96dIKNcet637eXg4LdIcFuP35ZRboteqrF/e+LNYBt0MAAABwOQwCBNBYVqtV58or7Ytzi8Oi3XEBb7nM4v50aUWTTF2XJB8vkwL9fRTo73vhTx+1MfvaF+Ftalw9t1+B97t4lb36CryP/Qq9nw+LdQBXRgAAAHA5DAIEPIvVatX58ir7otxyydX2mgt5y/lyWc7VOFZ68WdNtXj3ti/efRRo9lVQwMWFfFCNBX3tY74KuvCnv68XQ9oAtDgCAACAyzGZTOrXJVhrvzumHQwCBJya1WpVaUWVfdFec2t87avutr/Xrmmqr3zzMqmehXnNK/IXagJstRePBfh6s3gH4JIIAAAALqn/hQCAOQBoSVVVVlVaraqyWlVVJVVaraqssqqqqvpYZY3j9mP2P2X/+6XHqy7UV9Y4XvM8tuNWa/U5ah6vsl7oy+H9ZP+743s6Hrdabee+9HezqvKS89pfd+G44+9Qfdxa49y27fanzpc32VfAeZmkNuYLC/Y6FuaXLuCDLvkz0L96+zyLdwCeigAAAOCS+netngPwDd8EYKiqKqvKq6pUWWVVRZVVFZVWVdieV1Yfq6yqqvGz6uflldYar6m6cNyq8krHc13ptXWdy/73qqoL57Cq3NZHpe1njs8vfe+6njfV9nFPZbqweA+qZ6t8fVfdA/191ZrFOwBcEwIAAIBLaq5BgJdeha2+uqtax+yPS6/c1lFb+5wXrwRX2H5+6VVd+7FLeqr1nrp8Tw61cjhWe7Fd9yLeYXFdVaXKyhqL6wtXpHGRySR5m6q/j9zbZJK3l6n62IXntuNeJlX/vcZxL5PkdeE11a8zyfvCa23HveznuHjc/rMLx71s9bZar+qeTDXOXf0aObz+4vkuec8LvTXo9/K6+DsE+Ho7XJlv7efD97QDgMEIAAAALqnmIMB7F2yUr7fJYbuzbRt0zWMVVZcu1h23bHOFt2mYTJKvl5e8vUzy8TLJx9skby8v+999Liw2fb1r1lz8e10/u3jcdOH4hWMXzlfzuf29vWscr/G8Zl+25w7nraPHmgviSxe6F/8UV6cBAE6NAAAA4JJMJpNu7d1en35zWHmHLS32vjUXgBevsFYvGC9dENoWhbWPOS4kbcd8Lvl5zauttiu2l14NvvSKbnW9V/XV3Vqvv3hO2yK39uL6wkL9Mgt3nxqLa2+vC4vtSxbhXOUFAMA5EQAAAFzWiz/rr3tv7iKrZN/O7LAl+TKL7Uu3Ml9pYV7zNQAAAK6KAAAA4LIC/X31k8hORrcBAADgEryMbgAAAAAAADQ/AgAAAAAAADwAAQAAAAAAAB6AAAAAAAAAAA9AAAAAAAAAgAcgAAAAAAAAwAMQAAAAAAAA4AEIAAAAAAAA8AAEAAAAAAAAeAACAAAAAAAAPAABAAAAAAAAHoAAAAAAAAAAD0AAAAAAAACAByAAAAAAAADAAxAAAAAAAADgAQgAAAAAAADwAM0WAOzfv1+PPPKIIiIiFBAQoN69e+v3v/+9ysrK6qw/fvy4unbtKpPJpJMnTzr8bPv27YqPj1dAQIC6dOmiWbNmyWq1OtSsXbtWUVFR8vf3V69evbRw4cJa75GRkaHIyEiZzWZFRkZq+fLltWoWLFigiIgI+fv7KyoqSuvXr7/6DwEAAAAAACfRbAHArl27VFVVpTfffFM7d+7UK6+8ooULF2rGjBl11j/yyCMaMGBAreMWi0XDhw9XeHi4Nm/erLS0NM2bN0/z58+31+zbt0+jRo1SXFyccnJyNGPGDE2ePFkZGRn2muzsbCUmJiopKUm5ublKSkrS2LFjtWnTJnvNsmXLlJKSopkzZyonJ0dxcXG66667lJ+f34SfDAAAAAAALc9kvfRSejN6+eWXlZ6err179zocT09P17Jly/Tss8/qzjvv1IkTJ9S2bVv7z6ZPn64jR47IbDZLkl588UWlpaXp4MGDMplMmjZtmlasWKG8vDz7OZOTk5Wbm6vs7GxJUmJioiwWi1atWmWvGTlypEJCQrRkyRJJ0uDBgzVo0CClp6fba/r27asxY8Zozpw5DfodLRaLgoODVVJSoqCgoMZ/SAAAAAAANEJD16EtOgOgpKRE7dq1czj27bffatasWXr//ffl5VW7nezsbMXHx9sX/5I0YsQIFRYWav/+/faahIQEh9eNGDFCW7ZsUXl5eb01GzdulCSVlZVp69attWoSEhLsNXUpLS2VxWJxeAAAAAAA4GxaLAD4/vvvlZaWpuTkZPux0tJSPfDAA3r55ZfVvXv3Ol9XVFSkTp06ORyzPS8qKqq3pqKiQsXFxfXW2M5RXFysysrKemvqMmfOHAUHB9sf3bp1u2wtAAAAAABGaXQA8Nxzz8lkMtX72LJli8NrCgsLNXLkSN1///169NFH7cenT5+uvn376sEHH6z3PU0mk8Nz210LNY9fbc2lxxpSU9P06dNVUlJifxQUFNT7uwAAAAAAYASfxr5g0qRJGjduXL01PXv2tP+9sLBQw4YNU0xMjN566y2Hui+++ELbt2/XRx99JOnioj00NFQzZ87U888/r7CwsFpX4I8ePSrp4k6Ay9X4+Pioffv29dbYzhEaGipvb+96a+piNpsdbk8AAAAAAMAZNToACA0NVWhoaINqDx06pGHDhikqKkrvvPNOrXv8MzIydO7cOfvzzZs36+GHH9b69evVu3dvSVJMTIxmzJihsrIy+fn5SZIyMzMVHh5uDxpiYmL0ySefOJw7MzNT0dHR8vX1tddkZWVp6tSpDjWxsbGSJD8/P0VFRSkrK0v33nuvvSYrK0ujR49u0O8LAAAAAICzanQA0FCFhYUaOnSounfvrnnz5unYsWP2n4WFhUmSfZFvY7tfv2/fvvZvARg/fryef/55TZgwQTNmzNDu3bs1e/ZsPfvss/at+cnJyXr99deVmpqqxx57TNnZ2Vq0aJF9ur8kTZkyRbfffrvmzp2r0aNH6+OPP9bq1au1YcMGe01qaqqSkpIUHR1t37GQn5/vMLfgSmy7GBgGCAAAAABoCbb15xW/5M/aTN555x2rpDofl7NmzRqrJOuJEyccjn/zzTfWuLg4q9lstoaFhVmfe+45a1VVlUPNl19+aR04cKDVz8/P2rNnT2t6enqt8//tb3+z3njjjVZfX19rnz59rBkZGbVq3njjDWuPHj2sfn5+1kGDBlnXrl3bqN+7oKDgsr83Dx48ePDgwYMHDx48ePDg0VyPgoKCeterJqv1ShEBGqOqqkqFhYUKDAysd3ig0SwWi7p166aCgoJ6vycSuBb8O0Nz498YWgL/ztAS+HeG5sa/MfdmtVp16tQphYeH17r1vqZmuwXAU3l5ealr165Gt9FgQUFB/AcAzY5/Z2hu/BtDS+DfGVoC/87Q3Pg35r6Cg4OvWNPorwEEAAAAAACuhwAAAAAAAAAPQADgocxms37/+9/LbDYb3QrcGP/O0Nz4N4aWwL8ztAT+naG58W8MksQQQAAAAAAAPAA7AAAAAAAA8AAEAAAAAAAAeAACAAAAAAAAPAABAAAAAAAAHoAAAAAAAAAAD0AA4IEWLFigiIgI+fv7KyoqSuvXrze6JbiROXPm6JZbblFgYKA6duyoMWPG6D//+Y/RbcHNzZkzRyaTSSkpKUa3Ajdz6NAhPfjgg2rfvr1atWqlm2++WVu3bjW6LbiJiooKPf3004qIiFBAQIB69eqlWbNmqaqqyujW4MLWrVunn/70pwoPD5fJZNI//vEPh59brVY999xzCg8PV0BAgIYOHaqdO3ca0yxaHAGAh1m2bJlSUlI0c+ZM5eTkKC4uTnfddZfy8/ONbg1uYu3atZo4caK++uorZWVlqaKiQgkJCTpz5ozRrcFNbd68WW+99ZYGDBhgdCtwMydOnNCQIUPk6+urVatW6dtvv9Uf/vAHtW3b1ujW4Cbmzp2rhQsX6vXXX1deXp5eeuklvfzyy0pLSzO6NbiwM2fO6Ec/+pFef/31On/+0ksvaf78+Xr99de1efNmhYWFafjw4Tp16lQLdwojmKxWq9XoJtByBg8erEGDBik9Pd1+rG/fvhozZozmzJljYGdwV8eOHVPHjh21du1a3X777Ua3Azdz+vRpDRo0SAsWLND//u//6uabb9Yf//hHo9uCm/jd736nf/3rX+yUQ7O555571KlTJy1atMh+7L777lOrVq305z//2cDO4C5MJpOWL1+uMWPGSKq++h8eHq6UlBRNmzZNklRaWqpOnTpp7ty5evzxxw3sFi2BHQAepKysTFu3blVCQoLD8YSEBG3cuNGgruDuSkpKJEnt2rUzuBO4o4kTJ+ruu+/WT37yE6NbgRtasWKFoqOjdf/996tjx44aOHCg/vSnPxndFtzIbbfdps8//1zfffedJCk3N1cbNmzQqFGjDO4M7mrfvn0qKipyWA+YzWbFx8ezHvAQPkY3gJZTXFysyspKderUyeF4p06dVFRUZFBXcGdWq1Wpqam67bbb1K9fP6PbgZtZunSp/v3vf2vz5s1GtwI3tXfvXqWnpys1NVUzZszQ119/rcmTJ8tsNuuXv/yl0e3BDUybNk0lJSXq06ePvL29VVlZqRdeeEEPPPCA0a3BTdn+f/661gMHDhwwoiW0MAIAD2QymRyeW63WWseApjBp0iR988032rBhg9GtwM0UFBRoypQpyszMlL+/v9HtwE1VVVUpOjpas2fPliQNHDhQO3fuVHp6OgEAmsSyZcv0wQcf6K9//atuuukmbdu2TSkpKQoPD9dDDz1kdHtwY6wHPBcBgAcJDQ2Vt7d3rav9R48erZUCAtfq17/+tVasWKF169apa9euRrcDN7N161YdPXpUUVFR9mOVlZVat26dXn/9dZWWlsrb29vADuEOOnfurMjISIdjffv2VUZGhkEdwd38v//3//S73/1O48aNkyT1799fBw4c0Jw5cwgA0CzCwsIkVe8E6Ny5s/046wHPwQwAD+Ln56eoqChlZWU5HM/KylJsbKxBXcHdWK1WTZo0SX//+9/1xRdfKCIiwuiW4IbuvPNObd++Xdu2bbM/oqOj9Ytf/ELbtm1j8Y8mMWTIkFpfY/rdd9+pR48eBnUEd3P27Fl5eTn+v+Pe3t58DSCaTUREhMLCwhzWA2VlZVq7di3rAQ/BDgAPk5qaqqSkJEVHRysmJkZvvfWW8vPzlZycbHRrcBMTJ07UX//6V3388ccKDAy07zgJDg5WQECAwd3BXQQGBtaaK9G6dWu1b9+eeRNoMlOnTlVsbKxmz56tsWPH6uuvv9Zbb72lt956y+jW4CZ++tOf6oUXXlD37t110003KScnR/Pnz9fDDz9sdGtwYadPn9aePXvsz/ft26dt27apXbt26t69u1JSUjR79mxdf/31uv766zV79my1atVK48ePN7BrtBS+BtADLViwQC+99JIOHz6sfv366ZVXXuHr2dBkLnf/2DvvvKMJEya0bDPwKEOHDuVrANHk/vnPf2r69OnavXu3IiIilJqaqscee8zotuAmTp06pWeeeUbLly/X0aNHFR4ergceeEDPPvus/Pz8jG4PLurLL7/UsGHDah1/6KGH9O6778pqter555/Xm2++qRMnTmjw4MF64403CNA9BAEAAAAAAAAegBkAAAAAAAB4AAIAAAAAAAA8AAEAAAAAAAAegAAAAAAAAAAPQAAAAAAAAIAHIAAAAAAAAMADEAAAAAAAAOABCAAAAAAAAPAABAAAAAAAAHgAAgAAAAAAADwAAQAAAAAAAB7g/wPpjHVVi5sCCgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1200x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(UBICs); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "cf14c59b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from UBIC import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "d9f763b0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n_lams = 3\n",
      "Using baye_uncertainties\n",
      "[158.85535149 258.75771253 124.57513916  32.35651272   1.\n",
      "   1.02600143   1.06362583   1.0884958    1.23178153   1.36129565\n",
      "   1.40364847   2.05389896]\n",
      "0 ---> 8\n",
      "1.3333333333333333 ---> 6\n",
      "2.6666666666666665 ---> 5\n",
      "4.0 ---> 5\n",
      "5.333333333333333 ---> 5\n",
      "[ 473822.03308536  897656.69719111  306892.41629262 -116652.88906709\n",
      " -321592.7852043  -321488.21533194 -321333.28686819 -321228.78526099\n",
      " -320611.46574409 -320052.44972782 -319865.30952148 -317079.07092362]\n",
      "The optimal complexity is currently at the support sizes of 5.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([ 473822.03308536,  897656.69719111,  306892.41629262,\n",
       "        -116652.88906709, -321592.7852043 , -321488.21533194,\n",
       "        -321333.28686819, -321228.78526099, -320611.46574409,\n",
       "        -320052.44972782, -319865.30952148, -317079.07092362]),\n",
       " 4,\n",
       " 2.6666666666666665)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_lams = 3\n",
    "ics, bc, lam = UBICs(best_subsets, \n",
    "                     (X_pre, y_pre[:, 1:2]), \n",
    "                     u_type='cv', take_sqrt=True, use_baye=True, \n",
    "                     delta=1+1/3, n_lams=n_lams)\n",
    "ics[-n_lams], bc, lam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "0d69ff4f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.247674455990428"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "_, _, b_uns = baye_uncertainties(best_subsets, (X_pre, y_pre[:, 1:2]), u_type='cv', take_sqrt=True)\n",
    "lower_bounds = []\n",
    "for efi in best_subsets:\n",
    "    com = len(efi)\n",
    "    w = np.linalg.lstsq(X_pre[:, efi], y_pre[:, 1:2], rcond=None)[0]\n",
    "    pred = X_pre[:, efi]@w\n",
    "    lower_bound = log_like_value(pred, y_pre[:, 1:2])/np.log(len(y_pre[:, 1:2]))-com\n",
    "    lower_bounds.append(lower_bound)\n",
    "max_lam = max(np.log10(lower_bounds/b_uns))\n",
    "delta = np.linspace(0, max_lam, 4)[1]\n",
    "max_lam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "84cef451",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n_lams = 3\n",
      "Using baye_uncertainties\n",
      "[158.85535149 258.75771253 124.57513916  32.35651272   1.\n",
      "   1.02600143   1.06362583   1.0884958    1.23178153   1.36129565\n",
      "   1.40364847   2.05389896]\n",
      "0 ---> 8\n",
      "1.4158914853301425 ---> 6\n",
      "2.831782970660285 ---> 5\n",
      "4.247674455990428 ---> 5\n",
      "[ 787960.07394269 1409352.65270575  553240.99559745  -52667.5626646\n",
      " -319615.2752375  -319459.28728199 -319229.95618283 -319076.27396736\n",
      " -318175.60549507 -317360.47401245 -317089.58067781 -313017.465258  ]\n",
      "The optimal complexity is currently at the support sizes of 5.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([-167172.9175316 , -146452.83771546, -195778.94116253,\n",
       "        -247214.18952482, -325627.87086299, -325628.21898013,\n",
       "        -325625.10821244, -325620.95905375, -325581.80972334,\n",
       "        -325545.39428185, -325529.1513415 , -325366.7291644 ]),\n",
       " 4,\n",
       " 2.831782970660285)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_lams = 3\n",
    "ics, bc, lam = UBICs(best_subsets, \n",
    "                     (X_pre, y_pre[:, 1:2]), \n",
    "                     u_type='cv', take_sqrt=True, use_baye=True, \n",
    "                     max_lam=max_lam, \n",
    "                     delta=delta, \n",
    "                     n_lams=n_lams)\n",
    "ics[-n_lams], bc, lam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "3e652b6b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.831782970660285"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5abd37e7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eac9b889",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "c1bff1e5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4 <---> 4 inf\n",
      "4 <---> 5 2.0216252251847544e-05\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(2.8317829707,\n",
       " array([ 787960.07394269, 1409352.65270575,  553240.99559745,\n",
       "         -52667.5626646 , -319615.2752375 , -319459.28728199,\n",
       "        -319229.95618283, -319076.27396736, -318175.60549507,\n",
       "        -317360.47401245, -317089.58067781, -313017.465258  ]),\n",
       " 4)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tau = 3\n",
    "thres = 0.02\n",
    "\n",
    "post_means, b_bics, b_uns = baye_uncertainties(best_subsets, (X_pre, y_pre[:, 1:2]), u_type='cv1', take_sqrt=True)\n",
    "preddictions = X_pre@post_means\n",
    "lower_bounds = []\n",
    "for k, efi in enumerate(best_subsets):\n",
    "    com = len(efi)\n",
    "    lower_bound = log_like_value(preddictions[:, k:k+1], y_pre[:, 1:2])/np.log(len(y_pre[:, 1:2]))-com\n",
    "    lower_bounds.append(lower_bound)\n",
    "    \n",
    "last_lam = np.log10(max(lower_bounds/b_uns))\n",
    "delta = last_lam/tau\n",
    "now_lam = max_lam-delta\n",
    "last_ubic = UBIC(b_bics, b_uns, len(y_pre[:, 1:2]), hyp=10**last_lam)\n",
    "last_bc = np.argmin(last_ubic)\n",
    "while now_lam > 0:\n",
    "    now_ubic = UBIC(b_bics, b_uns, len(y_pre[:, 1:2]), hyp=10**now_lam)\n",
    "    now_bc = np.argmin(now_ubic)\n",
    "    \n",
    "    diff_com = now_bc-last_bc\n",
    "    imp = abs((b_bics[now_bc]-b_bics[last_bc])/(b_bics[last_bc]*diff_com))\n",
    "    print(min(last_bc, now_bc), '<--->', max(last_bc, now_bc), np.nan_to_num(imp, nan=np.inf))\n",
    "    \n",
    "    if (diff_com > 0 and imp < thres) or (diff_com < 0 and imp > thres):\n",
    "        break\n",
    "    \n",
    "    last_lam = now_lam\n",
    "    now_lam = last_lam-delta\n",
    "    last_ubic = now_ubic\n",
    "    last_bc = now_bc\n",
    "\n",
    "last_lam = round(last_lam, 10)\n",
    "last_lam, last_ubic, last_bc"
   ]
  }
 ],
 "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
}