Examples/GS-reproducible.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5c8736f7-3fcd-49e0-9367-6a3627c14b26",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sklearn's version: 1.2.2\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import warnings; warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import os\n",
    "import sys; sys.path.append('../')\n",
    "import numpy as np\n",
    "from numpy.random import default_rng\n",
    "\n",
    "import scipy.io as sio\n",
    "import pysindy as ps\n",
    "from PDE_FIND import build_linear_system, print_pde, TrainSTRidge, measure_pce\n",
    "from best_subset import *\n",
    "from frols import frols\n",
    "from UBIC import *\n",
    "from solvel0 import solvel0\n",
    "from findiff import FinDiff\n",
    "import sgolay2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ddb3bc6",
   "metadata": {},
   "source": [
    "#### Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0e366b04-cdf3-40e8-9cfc-b0e0c25ea165",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 128\n",
    "# data = sio.loadmat(f\"../../PDE-FIND/Datasets/reaction_diffusion_3d_{n}.mat\")\n",
    "data = sio.loadmat(f\"../Datasets/reaction_diffusion_3d_{n}.mat\")\n",
    "\n",
    "u_sol = (data['usol']).real\n",
    "v_sol = (data['vsol']).real\n",
    "x = (data['x'][0]).real\n",
    "y = (data['y'][0]).real\n",
    "z = (data['z'][0]).real\n",
    "t = (data['t'][0]).real\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",
    "\n",
    "# Ground truth\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])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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",
    "\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": "markdown",
   "id": "c7189a62",
   "metadata": {},
   "source": [
    "#### Add noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7ad8b85d-f21a-491d-ad42-d1c604f86f66",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(100)\n",
    "noise_lv = 0.1\n",
    "domain_noise = 0.01*np.abs(noise_lv)*np.std(u)*np.random.randn(*u.shape)\n",
    "u = u + domain_noise"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "874e4593",
   "metadata": {},
   "source": [
    "#### Denoise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "13716122",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|████████████████████████████████████████████████████████████████| 100/100 [00:15<00:00,  6.54it/s]\n",
      "100%|████████████████████████████████████████████████████████████████| 100/100 [00:12<00:00,  8.01it/s]\n"
     ]
    }
   ],
   "source": [
    "denoise = True\n",
    "if denoise: \n",
    "    un = u[:, :, :, :, 0]\n",
    "    vn = u[:, :, :, :, 1]\n",
    "\n",
    "    div = 30\n",
    "    ws = max(un.shape[:-2])//div; po = 5\n",
    "    if ws%2 == 0: ws -=1\n",
    "\n",
    "    nun = np.zeros_like(un)\n",
    "    for i in trange(un.shape[-1]):\n",
    "        for j in range(un.shape[-2]):\n",
    "            nun[:, :, j, i] = sgolay2.SGolayFilter2(window_size=ws, poly_order=po)(un[:, :, j, i])\n",
    "    un = nun.copy()\n",
    "    del nun\n",
    "    \n",
    "    nvn = np.zeros_like(vn)\n",
    "    for i in trange(vn.shape[-1]):\n",
    "        for j in range(vn.shape[-2]):\n",
    "            nvn[:, :, j, i] = sgolay2.SGolayFilter2(window_size=ws, poly_order=po)(vn[:, :, j, i])\n",
    "    vn = nvn.copy()\n",
    "    del nvn\n",
    "    \n",
    "    dim = 10\n",
    "\n",
    "    un = un.reshape(-1, len(t))\n",
    "    uun, sigmaun, vun = np.linalg.svd(un, full_matrices=False); vun = vun.T\n",
    "    un = uun[:,0: dim].dot(np.diag(sigmaun[0:dim]).dot(vun[:,0:dim].T))\n",
    "    un = un.reshape(len(x), len(y), len(z), len(t))\n",
    "\n",
    "    vn = vn.reshape(-1, len(t))\n",
    "    uvn, sigmavn, vvn = np.linalg.svd(vn, full_matrices=False); vvn = vvn.T\n",
    "    vn = uvn[:,0: dim].dot(np.diag(sigmavn[0:dim]).dot(vvn[:,0:dim].T))\n",
    "    vn = vn.reshape(len(x), len(y), len(z), len(t))\n",
    "\n",
    "    u = np.stack([un, vn], axis=-1)\n",
    "    del un, vn"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4add67a5",
   "metadata": {},
   "source": [
    "#### Applying best-subset regression on the weak formulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "65df608d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 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()\n",
    "# X_pre, y_pre = model.feature_library.cached_xp_full[0], model.cached_x_dot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "1f872310",
   "metadata": {},
   "outputs": [],
   "source": [
    "### To reproduce the exact %CE reported in the paper ###\n",
    "X_pre = np.load(\"./denoised_files/GS_X_pre_noise1e-1_denoised.npy\")\n",
    "y_pre = np.load(\"./denoised_files/GS_y_pre_noise1e-1_denoised.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "19450be1",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████████████████████████████████████████████████████████████| 12/12 [00:04<00:00,  2.89it/s]\n",
      "100%|██████████████████████████████████████████████████████████████████| 12/12 [00:04<00:00,  2.40it/s]\n"
     ]
    }
   ],
   "source": [
    "max_complexity = 12\n",
    "\n",
    "target_index = 0\n",
    "frols_coeffs,selected_inds_list,_,_,_ = frols(y_pre[:, target_index:target_index+1], X_pre, \n",
    "                                              max_nonzeros=max_complexity)\n",
    "selected_inds_list = set(selected_inds_list)\n",
    "ever = reduce(np.union1d, selected_inds_list)\n",
    "_, bs = brute_force_all_subsets(X_pre[:, ever], \n",
    "                                y_pre[:, target_index:target_index+1], \n",
    "                                max_support_size=max_complexity)\n",
    "best_subsets_u = [tuple(ever[list(b)]) for b in bs]\n",
    "del bs\n",
    "\n",
    "target_index = 1\n",
    "frols_coeffs,selected_inds_list,_,_,_ = frols(y_pre[:, target_index:target_index+1], X_pre, \n",
    "                                              max_nonzeros=max_complexity)\n",
    "selected_inds_list = set(selected_inds_list)\n",
    "ever = reduce(np.union1d, selected_inds_list)\n",
    "_, bs = brute_force_all_subsets(X_pre[:, ever], \n",
    "                                y_pre[:, target_index:target_index+1], \n",
    "                                max_support_size=max_complexity)\n",
    "best_subsets_v = [tuple(ever[list(b)]) for b in bs]\n",
    "del bs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d563c662",
   "metadata": {},
   "outputs": [],
   "source": [
    "u_pre, v_pre = y_pre[:, 0:1], y_pre[:, 1:2]\n",
    "del y_pre"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6db6f1aa",
   "metadata": {},
   "source": [
    "#### Model selection by the UBIC (Algorithm 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "fca396c4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.13696923174845235\n",
      "5 <---> 5 inf\n",
      "5 <---> 5 inf\n",
      "5 <---> 5 inf\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.0,\n",
       " array([-199409.01571746, -206249.85744539, -219768.49973045,\n",
       "        -250412.62102061, -271318.6918666 , -313925.64324027,\n",
       "        -313922.93529288, -313915.76704176, -313908.28846223,\n",
       "        -313900.68000578, -313893.37111871, -313878.47810631]),\n",
       " 5)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tau = 3; per = 75 # 85\n",
    "scale = np.log(len(u_pre))\n",
    "# scale = 1\n",
    "\n",
    "post_means, b_bics, b_uns = baye_uncertainties(best_subsets_u, (X_pre, u_pre), u_type='cv1', take_sqrt=True)\n",
    "predictions = X_pre@post_means\n",
    "\n",
    "b_bics = np.array(b_bics)\n",
    "complexities = np.array([len(bs) for bs in best_subsets_u])\n",
    "d_complexities = complexities[decreasing_values_indices(b_bics)]\n",
    "d_bics = b_bics[decreasing_values_indices(b_bics)]\n",
    "thres = np.percentile(np.abs(np.diff(d_bics)/(np.diff(d_complexities)*d_bics[:-1])), per)\n",
    "# thres = 0.02\n",
    "print(thres)\n",
    "\n",
    "lower_bounds = []\n",
    "for k, efi in enumerate(best_subsets_u):\n",
    "    com = len(efi)\n",
    "    assert com == np.count_nonzero(post_means[:, k:k+1])\n",
    "    # lower_bound = 2*log_like_value(predictions[:, k:k+1], u_pre)/np.log(len(u_pre))-com # old version without scale\n",
    "    lower_bound = 2*np.abs(log_like_value(predictions[:, k:k+1], u_pre))-np.log(len(u_pre))*com\n",
    "    lower_bounds.append(lower_bound)\n",
    "    \n",
    "last_lam = np.log10(max(lower_bounds/(b_uns*scale)))\n",
    "delta = last_lam/tau\n",
    "now_lam = last_lam-delta\n",
    "last_ubic = UBIC(b_bics, b_uns, len(u_pre), hyp=10**last_lam, scale=scale)\n",
    "last_bc = np.argmin(last_ubic)\n",
    "while now_lam >= 0:\n",
    "    now_ubic = UBIC(b_bics, b_uns, len(u_pre), hyp=10**now_lam, scale=scale)\n",
    "    now_bc = np.argmin(now_ubic)\n",
    "    \n",
    "    diff_com = now_bc-last_bc\n",
    "    diff_bic = b_bics[now_bc]-b_bics[last_bc]\n",
    "    # better use imp = abs(diff_bic/(b_bics[min(last_bc, now_bc)]*diff_com))\n",
    "    # if generally now_bc < last_bc\n",
    "    imp = abs(diff_bic/(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 (diff_bic > 0 or imp < thres)) or \\\n",
    "        (diff_com < 0 and diff_bic > 0 and imp > thres):\n",
    "        break\n",
    "    \n",
    "    last_lam = now_lam\n",
    "    now_lam = round(last_lam-delta, 8)\n",
    "    last_ubic = now_ubic\n",
    "    last_bc = now_bc\n",
    "    \n",
    "best_bc = last_bc\n",
    "if abs((b_bics[last_bc]-b_bics[last_bc-1])/b_bics[last_bc-1]) < thres:\n",
    "    best_bc = best_bc - 1\n",
    "\n",
    "last_lam = round(last_lam, 10)\n",
    "last_lam_d_u = last_lam\n",
    "last_ubic_d_u = last_ubic\n",
    "last_bc_d_u = last_bc\n",
    "uns_u = b_uns\n",
    "last_lam, last_ubic, last_bc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "bcd05278",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.10598674457151297\n",
      "4 <---> 4 inf\n",
      "4 <---> 4 inf\n",
      "4 <---> 7 1.465320735032407e-05\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(1.51625528,\n",
       " array([-157262.09616402, -130309.21118815, -188006.82684939,\n",
       "        -245195.50010284, -325565.48189545, -325564.20781034,\n",
       "        -325558.74969488, -325553.04892461, -325504.96014557,\n",
       "        -325460.46445174, -325441.57916254, -325238.58852882]),\n",
       " 4)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tau = 3; per = 75 # 85\n",
    "scale = np.log(len(v_pre))\n",
    "# scale = 1\n",
    "\n",
    "post_means, b_bics, b_uns = baye_uncertainties(best_subsets_v, (X_pre, v_pre), u_type='cv1', take_sqrt=True)\n",
    "predictions = X_pre@post_means\n",
    "\n",
    "b_bics = np.array(b_bics)\n",
    "complexities = np.array([len(bs) for bs in best_subsets_v])\n",
    "d_complexities = complexities[decreasing_values_indices(b_bics)]\n",
    "d_bics = b_bics[decreasing_values_indices(b_bics)]\n",
    "thres = np.percentile(np.abs(np.diff(d_bics)/(np.diff(d_complexities)*d_bics[:-1])), per)\n",
    "# thres = 0.02\n",
    "print(thres)\n",
    "\n",
    "lower_bounds = []\n",
    "for k, efi in enumerate(best_subsets_v):\n",
    "    com = len(efi)\n",
    "    assert com == np.count_nonzero(post_means[:, k:k+1])\n",
    "    # lower_bound = 2*log_like_value(predictions[:, k:k+1], v_pre)/np.log(len(v_pre))-com # old version without scale\n",
    "    lower_bound = 2*np.abs(log_like_value(predictions[:, k:k+1], v_pre))-np.log(len(v_pre))*com\n",
    "    lower_bounds.append(lower_bound)\n",
    "\n",
    "last_lam = np.log10(max(lower_bounds/(b_uns*scale)))\n",
    "delta = last_lam/tau\n",
    "now_lam = last_lam-delta\n",
    "last_ubic = UBIC(b_bics, b_uns, len(v_pre), hyp=10**last_lam, scale=scale)\n",
    "last_bc = np.argmin(last_ubic)\n",
    "while now_lam >= 0:\n",
    "    now_ubic = UBIC(b_bics, b_uns, len(v_pre), hyp=10**now_lam, scale=scale)\n",
    "    now_bc = np.argmin(now_ubic)\n",
    "    \n",
    "    diff_com = now_bc-last_bc\n",
    "    diff_bic = b_bics[now_bc]-b_bics[last_bc]\n",
    "    # better use imp = abs(diff_bic/(b_bics[min(last_bc, now_bc)]*diff_com))\n",
    "    # if generally now_bc < last_bc\n",
    "    imp = abs(diff_bic/(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 (diff_bic > 0 or imp < thres)) or \\\n",
    "        (diff_com < 0 and diff_bic > 0 and imp > thres):\n",
    "        break\n",
    "    \n",
    "    last_lam = now_lam\n",
    "    now_lam = round(last_lam-delta, 8)\n",
    "    last_ubic = now_ubic\n",
    "    last_bc = now_bc\n",
    "    \n",
    "best_bc = last_bc\n",
    "if abs((b_bics[last_bc]-b_bics[last_bc-1])/b_bics[last_bc-1]) < thres:\n",
    "    best_bc = best_bc - 1\n",
    "\n",
    "last_lam = round(last_lam, 10)\n",
    "last_lam_d_v = last_lam\n",
    "last_ubic_d_v = last_ubic\n",
    "last_bc_d_v = last_bc\n",
    "uns_v = b_uns\n",
    "last_lam, last_ubic, last_bc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e421b70a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQQAAACyCAYAAACk7LHWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAA630lEQVR4nO1de3gU1fl+TyKBYAJLggECmmVDA4SgmIsWawUkaC+iFgNKf1rpgyRqsbW2TYrWAm0thnprbdEEbWsrWEiqVVDQBBJFRSVZ0AQjgawBQgKBJBsSEliS/X5/nNnNXmYvM7PXZN7nmWdnzznv932zM3P2XL+PERFUqFChAgAigm2AChUqQgdqhaBChQorLgm2AcHCfffdR5MmTUJjYyO0Wi0AuDx3lxdszlDRqdrpP51r1659mYjuAwAQ0ZA8Vq9eTURElk/H8xUrVpAtXJVzde4ow1uOKxv8xXeX5y++O44v74F6D73jA1hDwnsx5LsMc+fOFT03mUxelXN17ijDW44rG/zFd5fnL747ji/vgXoPvedbQUH4dw6Fw7GGdcTbb7/tNt8bKJWh8oPLDwUbAsGHTQthyI4heMLIkSNl8bTPaXG08+hAwj7+kTQ6CY0PNwbEBpXvG34o2BBovlohuEB8fLws3tHOo6DVfG1HTU0NZs6cCQBga1nAbAgHvslkQkNDA3p6elyWiY2NRXV1tWz9SvmhYIMv+SNHjkRycjKioqJcllcrBBeIjY2Vze0z9+Gpj5/CLeNuCZoNoc5vaGiARqPB1KlTEREx5Iey/A6z2YxTp06hoaEB06dPd1luyN6JxsZGrFmzBpWVlS7z5SJvWx5W7VoF/WG9bBlKbQh1fk9PD8aNG6dWBgFCREQExo0b59QiE55/reX7kG0haLVarFmzxmW+u1rUE2paa9C9qhvGNiOICIxJ7y4otSEc+GplEFiI/d7CbEOjtUzArAkzyPl37Df3AwDe/uHbuDTqUuS8noOKxoqA2jCY+CoCD7VCcIHe3l5J5bce3Irbt9yOpNFJSHgqAWwtwydnPsH8f80HW8uQNDrJ7zYMNr4tXnsNSEsDIiP552uv+Uy0bOTl5cFoNPpVR2lpKcrLy1FcXCwrXyrUCsEFMjIyvC5b8XUFVr6zEn+Y9wc0PtwIWk2g1YTmlc2Ii47D0YePSp5ylGrDYORb8NprwGOPAc8/D5w/zz8feyz4lcKCBQuwdetWv8kvLS0FAGRnZwMAysvLJeXLgVohuIC3Uz0dvR1Y+t+l2JKzBVeNv8ou71DNIazPXo+ei66n1nxhw2DlW/DEE8DLLwPz5gHDhvHPl1/m6XJhMBhQXFwMvZ4P/K5fv16yjJycHJSUlMg3wgP27dsHnU4HANDpdFZbvc2XgyE7qGiZZZg7d67oMs6EhASPMs73nceY6DHQ5+mRGJsoKmNu6lx0m7pxoe8Chl8yXJKN3tgwWPiexl1vvFE6j9y4+rA09dvb2wHwl8uCgoICFBYWujfIBgaDATqdDgaDAYsXL0ZJSQni4uJQUFCAgoIC60tr0euqVZGbmytqowVtbW2S8r2BOssgwNMsg7vFGwBw+txp3PDPG7AlZwuuHHelWxn3vHEPfjDtB/jRVT+SZKMnGwYT393Lm5bGuwnz5g2kVVQADz0E1NbKsy09PR3r1q1DSUkJ9Ho9FixYAIC/3JbD9kUWQ2lpKRYvXoyioiIUFhZCp9NZD4CPMTjK0Gg0Ti++K2g0GmuFJSffGzjOMgzZCsETmpqaMGXKFNG8c6ZzuOW1W7Bo2iKXlYGtjOVXL8fvP/i95ArBnQ1DgW/BY48By5fzbsL11wMffsi/K+ky2KK8vBzZ2dkwGAyIi4tDVlaWV5WBRqNBTk6OpLESKS2ErKwsayvAYDBYKy1v82WB/LiBKJQPT5ubTp8+7TLvkZ2P0I//92Mym81eyejr76PJz02mT5s+dVteig3hzq+qqpIka/NmohkziCIi+OfmzYpMIyKiwsJCKisro/z8fCoqKqKOjg4qKyujsrIy6ujooI6ODkpPT3fiVVdXU1FRkZ2c6upqIiLKycmxppeUlPjMxsLCQmtadna223x3EPvdYbO5KSgvI4AcAGUOaR0AygDkO5TLBpDrizTbw1OF8NFHHzml9Zv7qfN8J3Vd6CJTn8kt31HGzsM76XDbYY8cTzYMFr7UCiFQqK6uppKSEmpoaLB+lwJLhWIrI5TgqUIISpeBiEoZY3kOyYuJyDpvwhjLEcqWM8ZyGWPZADRy02xle4P+/n6ntFXlq3Dq3Cn88/Z/SpZx85Sb0dzVDON5IzQjNLJtkIJw5wcD6enpSE9Pl823TAGGK0Jp2lHDGLPtuGUBMAjnBgDpCtMkwbJL0YLnPnkOb9W/hadvelq2jN9W/BYv7HtBNl8qwp0fClBSOYQjQmlQMQ5AO2OsiIjyIPzL2yBeYZodamtrsWzZMowdOxY6nQ6pqanIyMhAdXU1EhISsH//fkycOBFpaWl4//P3UVhViPeWvoeaz2owYcIEAEBLSwtmzZqFmpoaREZGIiUlBbW1tZg0aRJMJhPee+89LF++HNXV1YiOjsadujtx9/a7sXjiYvSe60VbW5tVZ0xMDBITE1FfX4/k5GS0tbWhsrLSytdoNIiPj0dDQwNSUlLQ3NyM7u5uKz8+Ph6xsbFobGzE9OnT0djYiL1791r5CQkJiIqKQlNTE9LS0lBfX4/+/n7MnDkTBw4cEL2mzz//HHfddZfdNbW2tlp1RkdHQ6vVoq6uDlqtFl1dXXbXdOjQISxcuNDumoxGo88WLKmQh5qaGuzatQunTp1CS0sLYDPt6M8xglyHI9uhTJkLbqHAL7RwwMcCCpWkOerxNIbw1VdfERHRya6TRETUY+rx3EFzIcMWs1+aTW/UvSGbr1R/qPBDdQxhsCMoYwhEVCqlPGMsF0AVEdkutdqHgX96HfiAo0ZBmkc4eTsSMGnUJBz/+XFvRHjEhu9vwGUjL/OJLBUqfI2gjCEIA3+ZloFDAFuFdMtAYqlQqegsg4REVK4kzRu7LN6OaDVh0zWbkBibiP8u+S+azjbJuk6hOWaHWeNnwXjeiGOdx2TxleoPJ34oI1Abm9ytLRgzZgwWLFhgXXZtNBqh1+tRWlqKgoICeUpJQldgMB1iXQasgfV8yX+W0Av7XnBKl4KOjg7R9Md3P04r314pm69UfyjwpXQZkp5NIqyB05H0bJIi+5SgpKTEbi2Cv2C75sARZWVldt+LioqsNlnWVjjCU5chlGYZQgr3j78f92fer0hGTU2NaHpeRh421WzC2QtnZfGV6g8XvgW2LTfbQ6x75y1sNzcZDAbJ/6j+3tjkDYxGIwwGg/V7bm6udbWjN0uvxaBWCC4wfJi0jUhiiIyMFE2fOGoi5uvm49+f/1sWX6n+cOH7E7abm3Q6HZKTkwFAcsVgeSENBgMyMjJgMBhgNBqRl5dn97JadBYXF4sectDe3o64uDjk5dkv6bEswZa1JoL80BwPh+Pee++l1atXU0VFxUDTyaZrYLvsVm6Xwd3S3WPGY9R5vlM2X6n+YPMdm66rK1bbdQeqTlRR1Ykqu7TVFauJiGjCUxOsaelFfGnxirdW2JU9cfaER/tslxlXV1dTQ0MD5eTkeLXC0NJlyM/PdynPF3DXZbAgPz/fbpm0u2XMjr97RUUFAfgnCe9FKK1DCCjEdjsmjU4SdZcux9sRwNc6uIqQc/noy7H3+F6MuGQErp5wtWS+Uv2hxl8zdw3WzF3jlE6rCWwts7q2B4DmXzQD4K7tq3O5z4XihcUoXijvn1av1yM9PR1GozFkNjZ5QnFxMTIzM50WTpWWliI/Px/AwHW5g7rb0Q1svRodOXJE8U69SZMmuc3/4tQX2NmwE2/c+YYsvlL9oc73N7KyslBeXo64uDgAQFVVlbViAID58+c7OXnR6/Vob29HTg6fILvzzjtFXzyDweCUJmXrM8B3YVZVVaG0tNSqb8GCBSgrK8OSJUtgMBisXpNycnJQXl6OgoICrFu3DgAk+XSwQK0QXMAxLqA/ZPzflf+HR3c/imOdx3DF6Ct8bkO48y3wdcvNAss/qQVxcXEwGAzWcYWNGzc6cRz3OtjKyMvLQ3l5OYxGo0+WPGdnZ6Ojo8MurayML6nRaDRWWyyVRXZ2NhoaGhTpVCsEF2htbUVqaqpfZcRExeCeK+/BKwdeweNzHve5DeHOt0COP0o5GOobmwC1QnAJX6y390bG2rlrMXKYePy9YDs5DTY/2BhqG5uAITzt6Clyky8chHojY/SI0fjw2IfYeWSnz20Id74K/0OyT0XG2Czw/QBGAO1EdMAfhgUannwqRkdHK9bhrYxuUzf+sOcP+M6U7/jUhlDnm81mNXpTAGE2m53SvIrcxBibzBh7kjH2AoA7ATAAYwAsYIy9yBhbxxjT+t7k0IFWqw2YjO9943toPdeKz0585lMbQpk/cuRInDp1SvQhVeF7WIK9egoP79RCYIzNB0BE9Gt3RMbYfMaYjoh2KzM1NFFXV4dx48YFREZkRCQezHwQb9S9gWsmXuMzG0KZn5ycjIaGBjQ3N8uWr0IaLOHg3UGsy1BFRJ2ehBPRLsbYaLnGhToC2UIAgIe/+TAuibC/HaH8D6+UHxUV5TEYbGNjoyIblPJDwYZA8526DLaVAWNsFGPsPldkbyqOcEVXV1dAZQyLHIYyQxle0r/kMxtUfmDv4WDguxpDuIMxtgXc81CVIovCFHKi4CiVMXbkWPzu/d+hz9znExtUfuDvYbjzGdHAGnHG2K8AZEJwWAKgmogaFVkUoli2bBlptVqXody6uroQGxurSIccGde9fB3yv5WP26fdrtgGlR+cexhO/MrKSsybN+8VIloGOLcQyonoTiL6LxH9F3xWwUVUvfCGZdrR1eabQK1DcMRD1zyET5s+9YkNKj849zCc+G43NxHRfofvGxljdyiyKEwRExMTFBmrdq3C0c6jePKjJ3nC+/wjaXSS5CW8Sq9hqPNDwYZA890uTBIWJXUwxhYBMAyWRUneIDHROZpzIGQc7TyKd+9+Fx8c/QAPTnvQKkNsc48/9Kv80LIh0HxXg4pXM8beA/AogCUAbgLwKGPs3cHahXBEfX190GTMuGwGNuzbgP1f7vdc2A/6VX7o2BBovqsWwmQiukksQ1i4NOjhaQGHP2VMHDUR2bps6Pv1+D6+H3D9Kj90bAg039VCcqMbDrnJGzQI9pTVT6/9KY61e3bV7i/9Kj80bAg031ULYQxj7EUARzBQOWgAJMPLoCehDstuR1fTjr7wua9ExvVXXI8+Qx+ICIxJHz9Qql/lh4YN/uY77na0W4dgC2FZ8hIMREAyAtg6WFYnrlmzhtztdgzWHLar6FFyZhlCfQ481PmhYEMg+IyxtUS0BnDjD4GIOoloIxH9STg2DpbKwBsEaw678eFGa9yBHdfvQHx0PBp+2iDLa1Coz4GHOj8UbAg0X/JmdKWzDIwxDWMsnTGWwxgrtEnPYYxlC3Ee/ZLmDV57DUhLA+bPn4O0NP5dLjQajXwygPHx4/HjWT/Ghn0bgqJ/qPNDwYZA80XHENxsaGLgkZmVbHleAgBEVMwYyxJe2HYhrZwxlmuJyejLNPIivuNrrwGPPQa8/DKQlNSEo0cvx/LlPG/pUukXGh/vFIVeMv/BpAexqWaTbL5S/UOZHwo2BJrvqoUwRjiYwwEAiroNRFRMRBYH+joABgBZwieEz3Q/pHnEE0/wyuBb3wLuuSca+/cDjz/O0+VAqQfchoYGTB4zGb+54TfoudgTFP1DmR8KNgSa72qWoRhAJhHtcsxgjBlEyksGY0wH7pKtnDG22CE7HgODmb5Ks0NtbS2WLVuGsWPHQqfTITU1FXV1c9Df/wHq6i7DLbcwVFaaoNePw4kTEbj11pNYsmQUxozRY8qUBAA8uvGsWbNQU1ODyMhIpKSkoLa2FpMmTYLJZMKZM2fQ1dWF6upqREdHQ6vVoq6uDlqtFl1dXWhra0NGRgaqq6sRExODxMRE1NfXIzk5GW1tbVb+5srNeObIMyi7owwGgwEpKSlobm5Gd3e3lR8fH4/Y2Fg0NjZi+vTpaGxstNOfkJCAqKgoNDU1IS0tDfX19ejv78fMmTNx4MABTJgwwemaOjs7cebMGbtram1tter0dE29vb1obm62uyaj0WjN12g0iI+PR0NDg+g19ff3o7Gx0e6aent7rfmerikiIgKHDh3yeJ/cXVN/fz8qKyvd3id319Tb22vHF7tP7q6ps7MTe/bscXuf3F3TmTNn8Omnnzrdp127duHUqVOWCN1a64tBPgyPZjnAuxW5Dke2Q5l8m/NCSz6AbOG7T9McbRSL/jxjBtHu3fx83759RES0axfRlClEzz5LdPPNRDExRDfcQPTHPxLp9UT9/c7hsjZv5rIiIsw0Ywb/LgcWG8xmM0376zSq/LpSFl8uhjo/FGwIBB820Z/94oadiErd5TPGcohovXCeDmAfBv7VdeBrHTQ+TvOIxx4Dli/n3Yb+/nOoqADuu493GZYuBR5+GOjpAT74ANi5k6cZjcDNNwPf+Q6wYAFQVjYwDtHf/wEiI+fIHofo7u4GADDGsDJrJZ7/7HnM0c6RzJeLoc4PBRsCzicFLQE5B/g/dgOAauGw/JPnC3m2LQefptkeYi0EIun/7gYD0YYNRLfeSjRqFFF0NNE99xDt3UvU2XmWiHirY8YMjxW1E86ePTtwfv4sPbf3Odl8ORjq/FCwIRB82LQQ5LzQWgCzAIySyg2lw1WFYIFtVGhvceECUUQE0S9+QfSNbxBddVUH7dtHZDLxdKlwtMFsNtMx4zHZfKX6hxo/FGwIBN+2QpDjFL+U+Dbo+MG881HOdE9UFDB9OvD97wNffgncdls3br0V+N73ADlxYx1t+Nr4NTI3ZuJ833lZfKX6hxo/FGwIlWlHlyCiTOHzaxqkLtgByF4uahmH2LMHuPvuPmzcCFRVASdOAI8+Cpw9K98G3Rgdrh5/NbYeFA8p7okvFUOdHwo2BJrvyh/CKMbYKEWWhDg8hXJrbGyUJXfpUj4I+dBDwLRpSSgoADZsAA4dAlpagJQU4IUXgL4+72x0xMpr+OAib+lJ50vBUOeHgg3+5jtubnI1TnAHgBuF8xvFyoT74WkM4eTJkx77Xp4gJmP/fqL584mmTSPato3IbJbG7zf305baLdRvFpnv9IIvBUOdHwo2BIIPV2MIgss0gE/TtTPGfgkvV/kNNvjr32XWLD41+dRTQH4+MH8+sN+FYyQxfgSLwKLpi1DV7Nk7fqj/O4U6PxRsCDTfscuwQPCDsBHcfVoy+KrFIYfe3l6/yWCMDzx+8QWwZAkfdFy2DGhq4vmWDVbXXXeN6AarrgtduPnVm3Gy+6Rfr0EuP29bHvK25flcf+mXpSj90u0SF7d8X9gw6PnkplkNYD6A+9yVCdfDU5chkHPYnZ1Ejz5KFBdHdPvtRFotX7vQ1naWdu8mmjzZeT1E7lu5tLZyrU/0+5pf1lBGZQ1lPtWf/14+FVUVyeb7wobByoenaUfLgCIR7SKil8TKDHYEci/9qFF8IPLAAaCiAujsBOrrgf37qzFvHl/16LjBauU1K/Fi1Ysw9ZsU6/c1P1uXjWxdts/0523LQ9bELORmSNrJrvpDkMF3Ne1YaPvFZmzBNu3GcJ6J8DTLkJCQoFiHVBmXXw50dQE7dvBKYNOmNADA9dcDdXX2ZWeOm4k37nzDKUCsEv2hyF//0XokxyUjJzUn4Pp9ISPU+Y6zDGLh4K8G3z5s+30ygAMOab8G8CKA15UYHCxYIje5QlRUlGIdcmRMn873S+zYAcyaNRpbtgAJCTzdEbPGz0Lpl6VYMmOJz/Tb4nD3Yfzs3z+DZoQGeRl5AAB9ix6aERruFbpFj/bedlQ3V6NoYZE1v6C8AADwwuwXrN9tZZQ1lCE5LtnjP/7h7sMoM5Sh7B55bjyDdQ/Die8YucnV2EEOgC3gS5RXAGhzyN8C4GoAV4vxw+Hwx9JlX8jYvJmPGezeTbRhQxWNGkU0caL4noqL/Rfp8mcuJ32z3mf6HfklB0tI92cdNbQ3EBFRR28HYQ2orKHMWi77X9lUcrDE+r2soYyy/5Vt1e9KhickPplo5ci1XynCYemxUj48jSEQ3624EcB68GnHbGEKEoyxyQCSiWg/OYR+G0xIS0sLigzbhU0rV6YjNha4eJHvpnTEJRGX4IHMB/D8Z8/7TL8jXzNCA4CvkgTg9N2S1t7b7lK/KxnG80a78voWvfW8oKwAd0y/w06PHPuVwhe/YTjx3TlZLSeim4joAeHF3yhEh34SgKNDk0GHYEb9WboUqK0F9uzZi6Ym4K67gLvvBsxm57IrMlbgja/eQOd5Z0dWvor6I/ZSxkXHDZyPiHPKd9Tv6cUuKCtARnGGtVIo1hfj+6PlB6lx1B8sGeHG93ovA3EvzH8iHh36a8mWhRn6+/uDLsPCf+opPtgoNuQxduRY1D5Qi9EjRvtNfyD4q769CpoRGqzYtgLF1cXITMzEyIiRAdPvLxnhxneqEITArl5BStlww8yZM4Muw8IfNgzYuhX4xz+AN990Ljc+ZjwKPyxEn7lPlK9Uv1h3wBbt58XzbfV7kqEZoUFhdiH0LXrkbc/D4tTFPrM/mDLCje9UIRDR64yxFYyx+xhjWsd8xphWyP8lEYXlDAPgedrxwIEDinUolWHLHz8eKCnhHpwOHbIvFxkRiTcPvYlth7b5VP/WPVtR+FEhDB0GFFcXw3jeiIIyPoNQUFZgTde36FFUXYRyQ7n1vKq5Cr9/5/fQt+g9yrAgNyMX6RP4SvklM5b49PcLloxQ53u1uYkPPGI0gF8B2Ao+q7AVwLsAfgkeDDboMwVKDk+zDF999ZWnwVmPUCpDjF9URDR9OpHjArTNX2ymef+c53f9/uZ39HZYZxaCbX8o2BAIPrzxqUg8StOf5NdNKvyB3Fxg3z6+96G0lO+LAIA7Uu9AQXkBms42YdKoSUG1UQk0IzTWWQgVgYdXg4qMscmMsV+KdSEGKwT31EGV4Yr/17/yjVBPPjmQFhUZhfqH6u0qA3/pHyr8ULAh4HzyonkNYJY35cLp8NRl6Ojo8NjU8gSlMtzxjx8nmjCBaOfOgTSz2Uw/feen1NHb4Xf9Q4EfCjYEgg+FPhUBiO9vGEyoqakJugx3/EmT+LboH/0IMAjjcowxnO45jX/s/4ff9Q8FfijYEGi+txUCY4xVMcbeFY73AJRIti6MEBkZGXQZnvhz5nA/jYsW8f0PAPDQNQ/hb/v+BjOZ/a5/sPNDwYZA8xlvMXgoxNhkcliMJJYWTli2bBlptVrMnTvXssHDDmfOnMHYsWMV6VAqwxs+EXDPPfzz1VeByX/W4mjnUadySaOTJIeUD4T9ocwPBRv8za+srMS8efNeIaJlgJctBLEXP5wrA2Bgt6NYZQDw2I9KoVSGN3zGgOJi7vb9z38GjnYehfm3ZtBqQsWcCtBqAq0m0UrCF/pdgYjw1idvyeYr1e8LfijY4G++425HsZWKv2KM7WOMrbPMKjDGRjPG5g/2cQNbTJqkfOpOqQxv+SNHAq+/Dqxbx78zxnDfW/fh4qiLAdEvhj3H9uAvjX8Jmn5f8EPBhkDzxVoI5QByiWgVETUC1n0MuwB8rbRSYIxpGGPpjLEcxlihTXoHY6yMMZZvk5bDGMtmjOX6Ik0KTCbXnogCJUMKf/Jk4N//5ufHj/P5/NePKFtIqsT+bYe24bpx1wVNvy/4oWBDoPliFcJkcrGtmfhiJfGtbd5jCXio+VIAsHlhFxPRAhoIApsj6CwXvmcrSZNqZGtrq/wr9JEMqfybbuKfc+cCpU/chhd3fyLqpNVf+m2xrX4brhpxlWy+O/15efyQy3eF0lJ+KJGh1IZg8+VMO2pkcKwgomIisnhy1gGwLGbXMMZs98hm2eQZwP0yKEmThIyMDKkUn8uQyz95Eph26WxcfVUECp/twmOPBVh/90l0m7px15y75Cn2oH/xYn7I5YuhoABobwdyHDy1BeseBosvtnT5Grh3i5YsSYMLCC9/u+VfHLzl0c4YKyKiPDhXPPEK0+xQW1uLZcuWYezYsdDpdEhNTUVGRgaqq6uRkJCA/fv3Y+LEiUhLS0N9fT36+/sxc+ZMHDhwABMmTADAV4HNmjULNTU1iIyMREpKCmprazFp0iSYTCa89957WL58OaqrqxEdHQ2tVou6ujpotVp0dXWhra3NqjMmJgaJiYmor69HcnIy2traUFlZaeVrNBrEx8ejoaEBKSkpaG5uRnd3t5UfHx+P2NhYRHZfgZ58hncBoBW4pXUUcC9wSfcVOHLkCJqamry+ps8//xx33XWX3TW1trZadbq6psunXY5nZzyL/7z2HyxcuNDumoxGo5Xv6ZoMBgNuvPFGNDY2Yvr06WhsbERvby+uvZbnf/llAqKiolxe07Fjx3Dttdd6vE8FBRosWTIGl1/+KT791P6adu/eDZ1O5/Y+ubumbdu2YerUqU73yfGabJ8922v6z3/+g6uuukrys2e5Ty+//DJmz57tdJ927dqFU6dOWVYyaq0vBjms4AN3jbYFQKxD+igh3eOqRXAXbLkOR7ZDGdEw7eAOXnOET0uo+Gzhu+w0Rz2eVip+8sknHld4eYJSGXL4ERE82nRTE5E2Yz/NfORXdO6cvOjTcu2f9895VN5QHpTrl8ovLORHMG0INh/uNjcR0X7G2EYA+xljBMAIYAwAAnA/8cjPbkHC+IArMMZyaGCsIB1AJoAqItLbFNuHgX96HXg0KY2CNEnQarVSKT6XIYc/fTrw4YfAvHnA/16NRebmFzF7zm+RkhITEP0dvR2oaq7C7Mtno+vSLrdl9XreVNdoBsYE9Hr+PTsbqKmZis8/B6qrgaIiew7AI2CJySgrA5KTgdtuc2+/Xs/Llrl5OoJxD4PJd+VTsZyIpgC4Gdxl2gIi+gbxmQZFEAb4Chlj1YyxavCuwlYhzzIYWCpUKjqhvEawSXaaVDvrHP2ey4BSGXL4lujTFRXA6aYmXBn/TRw2v4euLkBqVDA5+t9teBc3JN2AkcNGeuSnp/OXWK8HdDpeCeTm8jSDAZgy5QByc/m5ZbAvPX2gQnAlY9UqnuZJ/+LFAxWNKwTjHgaTL+aG/UYSwrwTkQEDg3NO+XIgvJxi4xB64Si1KbteOC33RZoUhGsLYelS/vnQQ0Bd3Q1IuKUG3132OW4wL8J113GPS1lZ7mUo0V9uKMfClIVe8zUa/qnTOX+PiNBa09rdOFxyJSMuzt6Po17PKxCAVyo5OQMcVwi3f3ilfLFBxfWMsX1uOJmwidswWNHV5b65GwgZcvlLl/KjpqYWaWk/AROcJmi1PI5kcTHwgx/4R/9D1zyE1MtSJfHFXsq4OOD48S7ruRwZtvoLCoD163n3Iz2d/wZfe7HWNlj3MFh8sS6DAQP/1mJHWC9Z9hZtbW1Bl+ELPmMMf/3sr9C36HHbbTwAzMqVwNNP8/0PvtT/4bEP8fC7D2NY5DBZfKX6HdHR0WE9X7WKtxxWrOCVQWbmQEvCnzaEG1+shbCO3MRbYIx5jkMeBrD4VHS1uSmc1yE48k92n8TWg1uRPiEdmZnA3r08+nRDA/CXvwCXuPCbJVX//776H+YkzZHMd9UdsPDddRfcybjqqoGFURoNUFg4sKjJ09iBow1yEep8r3wqgk8xam2+LwKwDsAisfLheIRq5CZ/8D85/gml/i3VLs9oJLrpJqLvftfZP6Nc/VOfn0pVJ6q85ldXE2VnE2k03FdkRwdRfj4RQJSbS7Rp014qKiLS6YjS04nKyjgnJ2eA407GLbecoAaHwE/p6TzPW78joXIP/cmHu2lHxtgd4FGbiDG2DkAngDzwwbn7GWNZRLRKUbUVBoiJkT5N52sZvuJnTcxC14UutHS1YEIsX9gyejSwfTvwk5/wYLJvv82drsjV32hsxNkLZ3H1hKu95qenO0/5FRbyAwCqqi5BZiafebBFiYMnDlcyqqqaodMl2uXt2sVbE950F7y5hsHGFxtDyCKiOCKKB68MiIgyiejXRHQTAKbIwjBBYmKi50J+luErfgSLQMNPG6yVgQXDhvGm8913A7NnA/v3i/O9wdiRY/Hu3e8igg08UqFy/bbQaDzPLPjbhlDmi1UI1hkGItoI50FEdzMQgwaDLQxYe287nv74aacyjAG/+hXw7LN8c9T27fL037H1DnSbul3ql4Ng80PBhlAI5eY49uw4TOnZxdIgQHKy8i0bSmX4kh8TFYO17691CrBqQU4OsG0bH4VftgxISwPmz5/j1W7Jjt4O7D2+F1eNt9/dGErXH642BJovViE8KjhHWSeMIayy+f4kgEE/fgAMnmlHCy6NuhQ3JN2AnUd2uiz/zW/y+frNm4Fp04CPP/4czz/PVz+6qxR2HtmJOdo5GDnMPhZjKF1/uNoQaL6rdQhVNsdWm/N9GCTrEDyFcjMajYp1KJXha/7ClIVuKwQAeOklvky4rQ146ql4zJsHvPwyD1HvCjuO7MAt37jFo36pCDY/FGzwN9/jtCM8hGnzlB8uh6dpx7Ou5uIkQKkMX/N7TD10oe+CW45lt+Tp00Tx8f108CD/7m635IdHP6SuC10e9UuFLf/8eaKLF+XzfWHDYOXDQ1wGt20MCnPnqt6iuro66DJ8zY8eFo1dhl2oOeXaV79lt+TYscBddxnwi1/w79Oni5f/tOlT/PnTPyMmynl6y1f2m0zc5fyOHfL4vrBhqPDFKoSNQoTnUbaHIqvCEBpvJ6r9KMMf/M9OfIZ/HPiHS47tbsn/+79zqK0FfvhDuPS6VPJliXXvgjf6pcDCf/RR4LLLgFuceyVe8X1hw1Dhiy1aZeC+BMbYJTIWBx4RupyIzsozL3wQH+/kZCngMvzBv3XqrVhSugRP3/S0ddOTLex3S16JSZOAiAhn12IWbKvfhs2LNnutXwri4+OxfTuwdStfIyFirke+UoTiPfQnX6yFsIKIDhDRfodjFxG9Du6BaNCjoaEh6DL8wZ81fhbO953HobZDLnlLlwK1tcCuXe+jsRGYMQN48UXncofbDqPb1I30CeIuK5Xav3dvE5Yv5zMcct6LwXoP/coniYNxGCT7Ge69915avXq1y7XeJ06c8DgY4wlKZfiLf7zzOJnNZq/5NTVEl11G1NZmn9/W00Y7D+8UYbrX7w1MJqLMzAu0bp1sEYP6HvqKX1FRQQD+Sa4GFRljN3qoQ3ziZDXY8BS5qbm5WbEOpTL8xR8zYgy2HtzqNT8tjXcZ1q61z1/5zkokaZIk6/cGv/0tcMklvcjP91zWH/p9JSPU+Y6Rm6Q4SIkHMBmAFw6wwx/d3d2eC/lZhr/4kRGRyN2ei/m6+Rg70nXcP1v+2rVAaipw//18xqGjtwPb67fjpVtfkqzfE3bu5EFnnn++FhER35IlQ4l+X8oIN74UBylFRJRFQjSnwY7B5A/BESMuGYH5k+fjncPveM2/7DLuZOSXv+TfXa1O9Ea/O5w4wZdOb9oEZGdfKZmvVL+vZYQbX6xCKCCijSKHYger4YTBPod969Rb8dYh98FYHfkrVwL19fwffFv9NqvvRDn6xdDXx6c4f/ITvu4g2L+fL2SEG1/MDfuQWHjkCYN9ymrR9EW4IekGSfyoKO567ZFHgIJN38Ot074vW78Yfvc77r3p0Ufl8ZXq94eMcOPLCeU2JBAbGxt0Gf7kjxo+Cuf7zuNw22FJ/IULgZhv6LGp7CDGRI8RYXmn3xHl5XwfxaZNQGSkdL5S/f6SEW78IVsheNrc1Cg1iIEfZPib/+ZXb+Ivn7oO2S7GZwyYtmQTPn5/BGx8mMrSb8HJk8CPfgT861/A+PHS+Ur1+1NGqPO98qk4FA5Pm5tOnjzpNt8bKJXhb37tqVq64tkrXK5JcMX/xl++QYt+Uk0PP6xMPxFRXx/R/PlEv/mNPL5S/f6WEQ58eNjcpAJD498l9bJURLJI1LSKb3YS4x86cwg9F3uw4fGr8e9/A4dcL3j0yv4//hG4eBFYvVoe3x2Gwj30Nd+FA27/QgixBvAQcQVCWg54HEkdCeHifZ0mBb29vfIuzocy/M1njGFLzhYkjRZfXCTGHxY5DH+c/0eMG8fw61/zacht2+Tpf/99YMMGoKpK3BV8sH+/ULAh4HwKcFMdQDqEyM/gQVh14NGec4S0XPD9Ej5Nc7RjKPpDEENffx/tPb7Xa/4jOx+hlq4WIuI+CpKTid59V7r+1laiiROJduxwbVuwf79QsCEU/CH4FUSkJ6L1jDENAAPx+JFZGIghaQCvNHydJglDZQ7bTGZ8d9N30dLV4pHf3tuOjfqNGD18NABg+HDgqaf4NGRfn3v9nZ02Os3APffw4zvfUWa/OwyVe+hLfjDHEDIBWLZiaRzy4v2QJgkJCQlSKT6XEQj+sMhhuDn5Zrx9+G2P/J1HdmKudi6ih0Vb0267DUhIADZudK3/44/tt0//6U9AVxdfd6DUfn/yQ8GGQPP9MoYg9N8dQ3QayCYsOxGVM8YW2/T1Hcv7Os0OtbW1WLZsGcaOHQudTofU1FRkZGSguroaCQkJaGlpQWtrK9LS0lBfX4/+/n7MnDkTBw4cwIQJPL5BS0sLZs2ahZqaGkRGRiIlJQW1tbWYNGkSTCYT9u/fj8svvxzV1dWIjo6GVqtFXV0dtFoturq60NbWZtUZExODxMRE1NfXIzk5GW1tbTh48KCVr9FoEB8fj4aGBqSkpKC5uRnd3d1Wfnx8PGJjY9HY2Ijp06ejsbERR44csfITEhIQFRWFpqYmp2tK7k/Ga/rX8O1Lv213TS0tLUhISLBe06v7XkVqdCq6urrsrumHPzyGX/86HVde+RUuXjxttamjowMajQZbtpzDuHETcODAEezZY8Yzz8zC889/ioMHo91e07lz5zBy5Ei7a+rt7bW7T66u6cCBAzCZTIiMjPR4n1pbW60yHe/T119/bZcvdp+MRqM13/E+HTp0yI4vdp/cXVNtbS3a2tokP3sWnfv370dXV5fTs7dr1y6cOnUKLS0tQDCnHQEUAsi1PUcIjiEMplBuntDR20E7Djt35h35qytWU2t3q6iMFSuIHnlEnL9iBdHf/kZ05gzRFVcQbdvmlVlB//1CwYZAh3ILRpehCIBBmGnQEFExEZUC0Nmklfs6TaqRaWlpii9UqYxA8TUjNJg9abbTOIIt/2DrQUyMnYjLLr1MVMbvfw+88grf6+DI/+ILvoV62TJg8WLuCq2nB/jsM9/Y7y9+KNgQaH4wBhUNwktbTkR5NunrhbT1/kqTgqEW9ecl/UtYU7nGJf/v+/+Olm7ngUcLxo0D8vN5FChbvtkMHDwI7NkDnD4NPP44jxKVnMyXKvvKfn/wQ8GGUIjcpAJAf39/0GUEkr9w6kJsP7wdZjKL8r3Z3fizn1lcrw3wv/4aiIkBnnsOuPlm7kvhgw+4B+ViD6tDgv37hYINgearFYILzJw5M+gyAslPiU9BbFQs9C16J75ldeKs8bPcyhg+nM8g/PznQH8/53/8MdDayqcaP/+cR5l+4w1glntRku33Bz8UbAg0f8hWCJ42Nx04cECxDqUyAs1/5uZnEB/NZ2hXvLUCb33M/SV0XujEz7/5c1EvzY74wQ/4moQrrgDi40fjkUf4hqX33gP+9z/g6qs9isBrr/Exh/j40V7FlvQ1PxRsCBRf3dzk5SzDV1995TbfGyiVEQz+qe5TRER09+t30xPvPEFERM98/Iw1MtPF/otuHbRu3sxXII4ZQ7R3bz3t3k00eTJP9wabN/Pyu3cT1dR8FXB+KNgQaD5sZhmCspdBRWjCTGbM2DAD+lw90i5LQ11THdp727Hm/TW4P/N+nOo+hTn/nINNizYhI1HcNdcTT3B/iJs2AT/5yUTMmAFMmQI88IB3kZfeegu45hrgH/8AOjvHY9QoPgD5wAPAO+49vgHg+yqysoC//x3o7JyAUaMAnY7z33ZeeyWK7duBzEw+6Hn2rL2M7dsHypGLOOjvvANkZPAxEgtfq+X8t9w7qQLAf6f0dKCoCOjqmoDYWCApifuzfPNN+7JiNuzcyfkvvghce20sHnmEx+d86KGBuBuuoFYILtDS0oKpU6cGVUag+REsAt+d8l1sq9+GtIQ0lOpLsePwDszVzsX5vvO46dWbsDRtqcvKAADq6oDrr+cvxJNPHsW0adPR18cdoGR7EdHj1VeBu+7im53q6pqRmjoafX18oPKmmzzzN23ibtg4/wRSU0ehrw/Yvdv9MmlbbN4M3H23RUYTUlNTrTK+9z37smK9qP/8B7j3Xmd+ZSVw662e9W/ZAvz4x878998Hbr/dubyjDSUlwH33cX539wkAibj+en5vPIIC2EwPpcNTl6Gjo8NtvjdQKiMY/K21W+k7r36H2nra6NGdj9KdJXfShs820OyXZtPPd/7cYzyHGTN4U9VW/+7dPN0bBJsfCjYEmg+bLkPQX8xgHZ4qhA8++MBtvjdQKiMY/M7znfTs3meJiGj5K8vpji130CfHP6E/f/Jnr4K72PZfd+/eo6j/Gwx+KNgQaL5aIXhRIXz00Udu872BUhnB4pvNZjpnOkcpT6fQFc9eQfnv5Uvib97M/40iIsw0Y4a0lzEU+KFgQyD5thUC49+HHpYtW0ZarRZz584Vjd505swZjB3rOoiJN1AqI1j8l/UvY8+xPfjy1Jc4e/EsvnjgC0RFRgVMf6jwQ8EGOXztc1oc7TzqlJ40OgmNDzfapVVWVmLevHmvENEyYAgPKlpCublCbW2tyzBv3kKpjGDwxR6m4X8YLvowecsHxB/GUOSHgg1K+Uc7j4JW8z/6yspK6zPA1jqPgHoTyk0FgEmTJgVdRjD4lofpyheuxG9m/QZLZi8BIP4wueMDwJEjRzBlypSw4oeCDWJ8IkLE7/g6QlO/Cf3mfpjJDDOZETs8Fuf7zqPb1G1den6h7wIiIyIxceJEL6+aY8iuVPQEk8kUdBnB5Bd8qwDmPjP+vv/vSHw6EQCQ+HQiKhsrcbzzOBKfTrQej+16DABw4ys3Wstds/EamEwmrNuzzo5fd7oOVc1VdnyLK/jUv6UCACY8PQG3bL4FJpMJP9vxM0x4eoI1vb23HW/Xv40JT0+wHltqt8DUb8L4p7gP9/FPjcfyN5fDZDLhztI7Me6pcQCAic/wl2Nj9UaMXT8W8evjEb8+HhVfV+BY5zGMKRyDMYU81kR+GY8ym7MjB6OfHI1R60ZZf5s1lWsQ88cY63Gw9SA+afoEI58YiegnuPOY9R/xPXXX/fc6jPjDCIz4wwgr/7637sPwPwy3Hm09bXi97nVE/T4KUb/nXbN/ff4vXOi7gKmbpoKtZdbKAAAWbVmEuPVxSHgqAUnPJVnLT/vrNKRt4Lsb9zbtRUtXC3ou9Ei78eTDgbpwOlR/COLAGtjxuy500YmzJwhrQCfOnqDei710sf8inTh7wnoYe41ExFc5Wsqd7DpJFRUV1Hm+045v6jPRhb4Ldvyz57nfv5auFsIaUPPZZjp97jRVVFRQR28HNZ9ttqb3m/upx9RDzWebrcc50zkym81WfktXC7X3tFNFRQW19bTRya6ThDWwrsLsMfXQ6XOn6fS503Tm3Bm60HeB+vr7qL2nndp72glrQOdM54iI6O2yt8nYayRjr9H625y/eJ66LnRZj77+Purr76NzpnN0znSOsAZk6jMREdG7u96l3ou91Hux18o39Zno/MXz1sNsNlO/uZ8u9F2gC30XCGtAff19ZDabqXx3uXV2x/beSLmHYul25dWVip6hBgrl/JioGMRExQAAEmMTrXm25xYkXJpglzcyYyRih8di1PBRXvHHx/B/+AmxE6z6Y0fEQjNCY5cePSzazo2bI9/ymZGRgdjogchFFvtc8W0jUVmC2H772m8jdoR99KPhlwzHcAx34o+MGAh8OyxyGABgdtZsjLhkhF05S54tGJjdwG1kBA9fdU3mNV7tIXEFXwR7VQHVQafKDw0b5PCTRieBrWVgaxlGPTPKeu7K3b4thmwLwbLb0dW045dffql4lkGpjGDwLQ+TFe8PpA8FfijYoJRvOxOxYcMGPPjggy7LqrsdvRxD+MUvfuE23xsolaHyg8sPBRsCwYcaym0Atv4QbM/PnDnjVTlX544yvOW4ssFffHd5/uK74/jyHqj30Hu+BUN2pSJj7CUATeDNpUYh2dW5u7xgc4aKTtVO/+mcRET3AUO4QlChQoUzhnyXQYUKFQNQKwQVKlRYoVYIImCM5TDGymRyNYyxdEFGoUwZ2cIhi28jRzafMdbBGCtjjOXL5Ft+gxzPpUW5DYyxauGQfB2C7mzGWK5UrsDPF/iSrl/s2ZFiixu+V8+jY1mpz6NaIYiAeOQnuVgCINMiQ+oDyRhLB5BOPNpUOmNMJ8cIIWqVLK6AxUS0gGQGugGwSvgN4mRcQxwRJRNRBoAV4NG+vIZw7ZZYogbhN5XCzwdgFPiQwnd8diwVoo0st47kxJ49Kc+jSFlJz6NaIfgYxEPTWUKQ6DAQlt5bvp6I1jPGNOAPtSQ+AAgvoGSeAzQKKqNcAPsYYzrh95D6G9iG3tPJ+A2qAJQIL7KOiPSeCA6IB9Bu890Lb5AukYWBe2EAIKlyUgqpz6NaIfgJwsvUTjLiSgrIBNAgkyvnJXJEHIB2xpikf2cByRBeKsZYkVC5SQZjLFdOa42IjOCtihIAcjZ0FAHIEuxOlsG3hcbhe7xCebLg7fOoVgj+Qw7ZxK6UCuHGJUvtgzPGshVUQrb6i4UXyyhnHABAg8CvBo/ALQcL5JCEZnk5ESVDhv3E448WgFeKDVDW2jIKcoINr55HtULwAxhjOZa+t4z+a6FNP88I6Q9TuzCAlQMeAVtyE5UxliuHZ4N9Nuca8OuQaoNGgf50m27COkj8DS2DcEIrK0vhmNI+DLQSdABkDVYrgZTnUa0QRCD8w2TKHCHPBlBoGSGH9Be6CHwgzBLK3kNIVHsIYxDlgl6NRN0WbAXsBsQkvRBCeY1lAE3qNQiIg30/XgqKhUotG8ASOb8hYL3+dVK4js+O8FvobO6n29ab2LMn5Xl0LCv1eVRXKqpQocIKtYWgQoUKK9QKQYUKFVaoFYIKFSqsUCsEFSpUWKFWCCpUqLBiyPpUHGwQppksU3VG8NWKcqb7fGmTRlicFDIQ5uS9nkYV1kNk+mKxVzhAbSEMAgjLUrOE1YWl4JWCJrhWAeAbaySBMSZ3ubY3siUvhbZUaHL3dYQb1AphcEAHoM3yRVhYo3Qvgy8gZ+m28oAYIlCy8lFoHchehh5OUCuEwYEqAKuEPfw6YGB1ocW3gHCuE1asZQvfswW/B1b/C8L+edF0gWPZ259jIydH8F+QzRgrscgAX62YI7ZcVrAlW7Cv0GIr+N4Hy/Jpi0+CIstybssKRCFPIybHBXIBODX7Lcu0PW1LRmi0uPwPCpDbc/Xw7wG+rbYMAIG/VBqbvDKb83wA2TbfrWUFGSWu0oXPfBtukU2ZMiFfZyvbjb2FNtwcR1stNoK/iNU2dhTapBe5kiOir0QkLdtGT66H37cQfI9E0O+1Pw+1hTBIQHwPwwIiYuD/hF57GSKhn0y8q5HtJj0PgK1vgQab8nGCDd52VbYAqBa2VztxaGAQbyO4kxQAuBNAm8XPgXC4leMB7eB+E4og7N8ArN6SHFs1+6DM4UxYQK0QBgEcm7vEt+4G6uG1bJZxuRHJxYCcAXy8oAy89SHGywbfw68XZGgA6IWKR09EC7yR4wpCRWcZs1hlk5VM4k5VjFLkhyPUCmFwQOOwO04D+39Lo82gWhYc+sM24wM62PSzRdIt3QYLsiDSL7eBpZIQqxByichIfKyjkIkP+hUCKBDOLd0Wq48E4V/cGzmA8zVnC7MOBkFHm012nDA+kWMjTwf5uy/DBuo6hEEEm0pBh4EXCeB97SWMMQP4Q53HGCungTUC2YwxI/jLtsKGZ5dOREbLIB6EPjwRWbZqZwovmO3aB8tgoGhTXrDXAO4qzmjpCggcDfhgaaaQHk9EBTb6IVyLkxwXP08Z4y7dLLZUCVzLuIej70gd7FsEySJlBh3U7c9DHIyxauLOTL1KD1cIrZxs8rBYS6hsFoN7fFpvk14odMUGNdQug4ohAaFloPGiqLV1ZTs9C4men8MVaoUwhCE86JYmusf0cAdxb9ZuvQ6R4EuSiNYTUbllDEHC7ElYQ+0yqFChwgq1haBChQor1ApBhQoVVvw/GaCzu7NPdpkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 252x189 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "assert list(map(len, best_subsets_u)) == list(map(len, best_subsets_v))\n",
    "complexities = list(map(len, best_subsets_u))\n",
    "with plt.style.context(['science', 'grid']):\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.plot(complexities, last_ubic_d_u, 'o-', c='blue', markerfacecolor='none', label=\"$u_t,\\, \\lambda_{\\\\textrm{U}}=$ \"+str(round(last_lam_d_u, 2)))\n",
    "    ax.plot(complexities, last_ubic_d_v, 's--', c='green', markerfacecolor='none', label=\"$v_t,\\, \\lambda_{\\\\textrm{U}}=$ \"+str(round(last_lam_d_v, 2)))\n",
    "    ax.set_xticks(complexities)\n",
    "    ax.set_ylabel(\"$\\\\textrm{UBIC}_{\\\\Gamma}(\\\\xi^{k}, 10^{\\\\lambda})$\", fontsize=12)\n",
    "    ax.set_xlabel(\"Support sizes ($s_{k}$)\", fontsize=12)\n",
    "\n",
    "    plt.annotate('min$\\checkmark$', fontsize=16, c='blue', \n",
    "             xy=(complexities[last_bc_d_u], last_ubic_d_u[last_bc_d_u]), \n",
    "             xytext=(complexities[last_bc_d_u]+0.5, last_ubic_d_u[last_bc_d_u]+0.5e5), \n",
    "             arrowprops={'arrowstyle': '->', 'linestyle':'-', 'color':'blue'})\n",
    "    \n",
    "    plt.annotate('min$\\checkmark$', fontsize=16, c='green', \n",
    "             xy=(complexities[last_bc_d_v], last_ubic_d_v[last_bc_d_v]), \n",
    "             xytext=(complexities[last_bc_d_v], last_ubic_d_v[last_bc_d_v]+0.925e5), \n",
    "             arrowprops={'arrowstyle': '->', 'linestyle':'--', 'color':'green'})\n",
    "    \n",
    "    plt.legend()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52bd7b8e",
   "metadata": {},
   "source": [
    "#### Percent coefficient error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "fafbbe39",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.055415449636699976 0.08751867999585022\n",
      "0.010965596710386558 0.010720272536196155\n"
     ]
    }
   ],
   "source": [
    "assert np.alltrue(best_subsets_v[last_bc_d_v] == ground_indices_v) and np.alltrue(best_subsets_u[last_bc_d_u] == ground_indices_u)\n",
    "errs_u = measure_pce(np.linalg.lstsq(X_pre[:, ground_indices_u], u_pre, rcond=None)[0].flatten(), \n",
    "                     ground_coeff_u)\n",
    "errs_v = measure_pce(np.linalg.lstsq(X_pre[:, ground_indices_v], v_pre, rcond=None)[0].flatten(), \n",
    "                     ground_coeff_v)\n",
    "print(errs_u.mean(), errs_u.std())\n",
    "print(errs_v.mean(), errs_v.std())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:pysr]",
   "language": "python",
   "name": "conda-env-pysr-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.10"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}