KS_exp/IPI-Mac.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6e4730cb",
   "metadata": {},
   "source": [
    "# PDE-FIND for the Kuramoto Sivashinsky Equation\n",
    "\n",
    "This notebook demonstrates PDE-FIND on the Kuramoto Sivashinsky equation.\n",
    "$$\n",
    "u_t + u_{xxxx} + uu_x + u_{xx} = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "dc63748a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sklearn's version: 1.2.2\n",
      "mrmr is not installed in the env you are using. This may cause an error in future if you try to use the (missing) lib.\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import warnings; warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import numpy as np\n",
    "from numpy.random import default_rng\n",
    "import pandas as pd\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, wiener\n",
    "from scipy import integrate\n",
    "from scipy.integrate import trapz\n",
    "from scipy.integrate import simpson as simp\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",
    "import pysindy as ps\n",
    "\n",
    "from os.path import join as join_path\n",
    "FIGURE_EXPORT_PATH = \"/Users/pongpisit/Documents/figure_export/\"\n",
    "\n",
    "def evaluate_coefficients(prediected_coeffs):\n",
    "    GROUND = np.array([0.1, -1])\n",
    "    errs = 100*np.abs(GROUND-np.array(prediected_coeffs).flatten())/np.abs(GROUND)\n",
    "    return errs.mean(), errs.std()\n",
    "\n",
    "import torch, sympytorch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "from sympy import symbols, simplify, lambdify\n",
    "from mathparser import math_eval\n",
    "from varname import nameof\n",
    "from misc import h5file\n",
    "\n",
    "import derivative\n",
    "from tvregdiff import TVRegDiff, tvregdiff, numdiff, pysindydiff, savgol_denoise\n",
    "from functools import partial\n",
    "from findiff import FinDiff\n",
    "\n",
    "from tqdm import trange\n",
    "import heapq\n",
    "\n",
    "def percent_coeff(pred, ground=np.array([-1., -1., -1.])): \n",
    "    return 100*np.abs(pred.flatten()-ground)/np.abs(ground)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9983aec0",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = sio.loadmat('../Datasets/kuramoto_sivishinky.mat')\n",
    "u = data['uu']\n",
    "x = data['x'][:,0]\n",
    "t = data['tt'][0,:]\n",
    "dt = t[1]-t[0]\n",
    "dx = x[2]-x[1]\n",
    "X, T = np.meshgrid(x, t)\n",
    "XT = np.asarray([X, T]).T\n",
    "\n",
    "assert np.max(x)/32 == np.pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b219ce45",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "noise_lv = 30\n",
    "un = u + 0.01*noise_lv*u.std()*np.random.randn(u.shape[0], u.shape[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "dcbb14dd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "STLSQ model: \n",
      "(x0)' = -0.994 x0_11 + -0.996 x0_1111 + -0.992 x0x0_1\n"
     ]
    }
   ],
   "source": [
    "library_functions = [lambda x: x, lambda x: x * x]\n",
    "library_function_names = [lambda x: x, lambda x: x + x]\n",
    "pde_lib = ps.PDELibrary(\n",
    "    library_functions=library_functions,\n",
    "    function_names=library_function_names,\n",
    "    derivative_order=4,\n",
    "    spatial_grid=x,\n",
    "    is_uniform=True,\n",
    ")\n",
    "\n",
    "print('STLSQ model: ')\n",
    "optimizer = ps.STLSQ(threshold=0.1, alpha=1e-5, normalize_columns=False)\n",
    "model = ps.SINDy(feature_library=pde_lib, optimizer=optimizer)\n",
    "model.fit(np.expand_dims(u, -1), t=dt)\n",
    "model.print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a1b8b8d8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MIOSR model: \n",
      "Set parameter Username\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2023-06-05 15:55:36,367 [INFO] Set parameter Username\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Academic license - for non-commercial use only - expires 2024-06-04\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2023-06-05 15:55:36,370 [INFO] Academic license - for non-commercial use only - expires 2024-06-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(x0)' = -0.994 x0_11 + -0.996 x0_1111 + -0.992 x0x0_1\n"
     ]
    }
   ],
   "source": [
    "print('MIOSR model: ')\n",
    "optimizer = ps.MIOSR(target_sparsity=3, normalize_columns=False, alpha=1e-1, regression_timeout=100)\n",
    "model = ps.SINDy(feature_library=pde_lib, optimizer=optimizer)\n",
    "model.fit(np.expand_dims(u, -1), t=dt)\n",
    "model.print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "007b854e",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Differentiator(ps.BaseDifferentiation):\n",
    "    def __init__(self, diff_name, diff_func_kwargs, \n",
    "                 d=1, axis=1, \n",
    "                 is_uniform=True, periodic=False):\n",
    "        super(Differentiator, self).__init__()\n",
    "        self.diff_name = diff_name\n",
    "        self.diff_func_kwargs = diff_func_kwargs\n",
    "        self.diff_func = getattr(derivative, self.diff_name)(**diff_func_kwargs)\n",
    "        \n",
    "        self.d = d\n",
    "        self.diff = partial(pysindydiff, **{\"diff_method\":self.diff_func, \"order\":self.d})\n",
    "        \n",
    "        # Other info...\n",
    "        self.axis = axis\n",
    "        self.is_uniform = is_uniform\n",
    "        self.periodic = periodic\n",
    "        self.transform = np.vectorize(composite_function(self.diff, lambda _: _, left2right=True), signature=\"(m),(m)->(m)\")\n",
    "\n",
    "    def _differentiate(self, x, t):\n",
    "        in_shape = x.shape\n",
    "        if len(in_shape) == 2: x = np.expand_dims(x, -1) # x should now be 3-dimensional\n",
    "        if isinstance(t, float) and self.is_uniform: \n",
    "            t = np.linspace(0, stop=t*(x.shape[self.axis]-1), num=x.shape[self.axis])\n",
    "        out = []\n",
    "        # wrt to x var\n",
    "        if self.axis == 0:\n",
    "            for i in range(x.shape[-1]):\n",
    "                diff = self.transform(x[:,:,i].T, t).T\n",
    "                out.append(np.expand_dims(diff, axis=-1))\n",
    "        # wrt to time var\n",
    "        elif self.axis == 1:\n",
    "            for i in range(x.shape[-1]):\n",
    "                diff = self.transform(x[:,:,i], t)\n",
    "                out.append(np.expand_dims(diff, axis=-1))\n",
    "        return np.concatenate(out, axis=-1).reshape(in_shape)\n",
    "\n",
    "class FiniteDifferentiator(ps.BaseDifferentiation):\n",
    "    def __init__(self, acc=2, d=1, axis=1, is_uniform=True, periodic=False):\n",
    "        super(FiniteDifferentiator, self).__init__()        \n",
    "        self.acc = 2*(acc//2)\n",
    "        self.d = d\n",
    "        self.axis = axis\n",
    "        self.is_uniform = is_uniform\n",
    "        self.periodic = periodic\n",
    "    def _differentiate(self, x, t):\n",
    "        if self.is_uniform:\n",
    "            if isinstance(t, float): \n",
    "                dt = t\n",
    "            else: \n",
    "                dt = t[1]-t[0]\n",
    "            self.diff = FinDiff(self.axis, dt, self.d, acc=self.acc)\n",
    "        else:\n",
    "            raise NotImplementedError(\"is_uniform=False\")\n",
    "            self.diff = FinDiff(self.axis, self.d, acc=self.acc)\n",
    "        return self.diff(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "57631251",
   "metadata": {},
   "outputs": [],
   "source": [
    "### originally from derivative ###\n",
    "## slow ##\n",
    "# diff_name, diff_func_kwargs = \"Kalman\", {'alpha':kalpha}\n",
    "## fast ##\n",
    "# diff_name, diff_func_kwargs = \"FiniteDifference\", {'k':1}\n",
    "# differentiation_method = Differentiator\n",
    "# differentiation_kwargs = {\"diff_name\": diff_name, \"diff_func_kwargs\": diff_func_kwargs}\n",
    "\n",
    "# originally from pysindy\n",
    "# differentiation_method, differentiation_kwargs = ps.SmoothedFiniteDifference, {}\n",
    "\n",
    "# originally from findiff (accurate and fast)\n",
    "differentiation_method, differentiation_kwargs = FiniteDifferentiator, {'acc':2}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ac2a2c65",
   "metadata": {},
   "outputs": [],
   "source": [
    "diff_order = 4; spatial_weights_p = 4\n",
    "weak_pde_lib = ps.WeakPDELibrary(library_functions=[lambda x: x, lambda x: x * x], \n",
    "                                 function_names=[lambda x: x, lambda x: x + x], \n",
    "                                 derivative_order=diff_order, p=spatial_weights_p, \n",
    "                                 spatiotemporal_grid=XT, \n",
    "                                 include_bias=False, is_uniform=True, K=10000, # new random K points in every calls to the ps.WeakPDELibrary\n",
    "                                 differentiation_method=differentiation_method, \n",
    "                                 differentiation_kwargs=differentiation_kwargs, \n",
    "                                 cache=True\n",
    "                                )\n",
    "kwargs = {'fit_intercept':False, 'copy_X':True, 'normalize_columns':False}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4dbe2eee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # optimizer = ps.MIOSR(target_sparsity=3, normalize_columns=False, alpha=1e-1, regression_timeout=100)\n",
    "# optimizer = L0BNB(max_nonzeros=3, lam=5e-2, is_normal=True, normalize_columns=False, threshold=1e-4) # tune lam จาก (X_pre, y_pre) | ตอนนี้ยังไม่ได้ tune\n",
    "# # optimizer = BruteForceRegressor(3)\n",
    "\n",
    "# # if feature_library=weak_pde_lib, then just differentiation_method=None is fine.\n",
    "# n_ensemble_models = 50\n",
    "# model = ps.SINDy(feature_library=weak_pde_lib, optimizer=optimizer, \n",
    "#                  differentiation_method=differentiation_method(**differentiation_kwargs), \n",
    "#                  cache=True,\n",
    "#                  feature_names=['u'])\n",
    "\n",
    "# model.fit(np.expand_dims(un, -1), t=dt, ensemble=True, library_ensemble=True, n_candidates_to_drop=1, n_models=n_ensemble_models)\n",
    "# print(model.get_feature_names())\n",
    "# print(model.print())\n",
    "\n",
    "# X_pre, y_pre = np.squeeze(model.feature_library.cached_xp_full), model.cached_x_dot\n",
    "# print_pde(model.get_coef_list()[np.argmin(np.sum((np.squeeze(np.tensordot(X_pre, np.array(model.get_coef_list()).T, axes=([-1], [0])), axis=1)-y_pre)**2, axis=0))].reshape(-1,1), model.get_feature_names())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "14c6d5ea",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['u', 'uu', 'u_1', 'u_11', 'u_111', 'u_1111', 'uu_1', 'uuu_1', 'uu_11', 'uuu_11', 'uu_111', 'uuu_111', 'uu_1111', 'uuu_1111']\n",
      "(u)' = -0.993 u_11 + -0.995 u_1111 + -0.998 uu_1\n",
      "None\n",
      "u_t = (-0.990625 +0.000000i)u_11\n",
      "    + (-0.993308 +0.000000i)u_1111\n",
      "    + (-0.995604 +0.000000i)uu_1\n",
      "   \n"
     ]
    }
   ],
   "source": [
    "optimizer = ps.MIOSR(target_sparsity=3, normalize_columns=False, alpha=5e-2, regression_timeout=100)\n",
    "n_ensemble_models = 50\n",
    "model = ps.SINDy(feature_library=weak_pde_lib, optimizer=optimizer, \n",
    "                 differentiation_method=differentiation_method(**differentiation_kwargs), \n",
    "                 cache=True,\n",
    "                 feature_names=['u'])\n",
    "\n",
    "model.fit(np.expand_dims(un, -1), t=dt, ensemble=True, library_ensemble=True, n_candidates_to_drop=1, n_models=n_ensemble_models)\n",
    "print(model.get_feature_names())\n",
    "print(model.print())\n",
    "\n",
    "X_pre, y_pre = np.squeeze(model.feature_library.cached_xp_full), model.cached_x_dot\n",
    "print_pde(model.get_coef_list()[np.argmin(np.sum((np.squeeze(np.tensordot(X_pre, np.array(model.get_coef_list()).T, axes=([-1], [0])), axis=1)-y_pre)**2, axis=0))].reshape(-1,1), model.get_feature_names())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "1c7fa09b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "u_t = (-0.992887 +0.000000i)u_11\n",
      "    + (-0.994861 +0.000000i)u_1111\n",
      "    + (-0.998153 +0.000000i)uu_1\n",
      "   \n",
      "0.46999004572032915\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{3: ((3, 5, 6),\n",
       "  0.01617495345907047,\n",
       "  array([[-0.99288685],\n",
       "         [-0.9948606 ],\n",
       "         [-0.99815284]]))}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_results = {}\n",
    "for effective_indices in set([tuple(np.nonzero(model.get_coef_list()[i][0])[0]) \n",
    "                              for i in range(len(model.get_coef_list()))]):\n",
    "    coeff = np.linalg.lstsq(X_pre[:, effective_indices], y_pre, rcond=None)[0]\n",
    "    mse = ((y_pre-X_pre[:, effective_indices]@coeff)**2).mean()\n",
    "    n_terms = len(effective_indices)\n",
    "    if (n_terms > 0) and (n_terms not in model_results or mse < model_results[n_terms][1]):\n",
    "        model_results[n_terms] = effective_indices, mse, coeff\n",
    "        \n",
    "for com in model_results:\n",
    "    effective_indices, _, coeff = model_results[com]\n",
    "    print_pde(coeff, np.array(model.get_feature_names())[list(effective_indices)])\n",
    "    \n",
    "print(percent_coeff(model_results[3][-1]).mean())\n",
    "model_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "d257da53",
   "metadata": {},
   "outputs": [],
   "source": [
    "# optimizer = L0BNB(max_nonzeros=5, lam=1e-1, is_normal=False, normalize_columns=False, threshold=1e-4)\n",
    "# optimizer.fit(X_pre, y_pre).coef_\n",
    "\n",
    "# mses = []\n",
    "# all_possible_models = list(itertools.combinations(range(14), 3))\n",
    "# for effective_indices in all_possible_models:\n",
    "#     coeff = np.linalg.lstsq(X_pre[:, effective_indices], y_pre, rcond=None)[0]\n",
    "#     mses.append(((y_pre-X_pre[:, effective_indices]@coeff)**2).mean())\n",
    "# all_possible_models[np.argmin(mses)], np.min(mses)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebb1ccd9",
   "metadata": {},
   "source": [
    "#### Evidence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "bfb21c1e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "u_t = (-0.333791 +0.000000i)uu_1\n",
      "   \n",
      "u_t = (-0.223216 +0.000000i)u_11\n",
      "    + (-0.574863 +0.000000i)uu_1\n",
      "   \n",
      "u_t = (-0.992887 +0.000000i)u_11\n",
      "    + (-0.994861 +0.000000i)u_1111\n",
      "    + (-0.998153 +0.000000i)uu_1\n",
      "   \n",
      "u_t = (-0.993005 +0.000000i)u_11\n",
      "    + (0.001229 +0.000000i)u_111\n",
      "    + (-0.994965 +0.000000i)u_1111\n",
      "    + (-0.998308 +0.000000i)uu_1\n",
      "   \n",
      "u_t = (0.003505 +0.000000i)u\n",
      "    + (-0.993824 +0.000000i)u_11\n",
      "    + (-0.995146 +0.000000i)u_1111\n",
      "    + (-0.998524 +0.000000i)uu_1\n",
      "    + (0.000302 +0.000000i)uuu_11\n",
      "   \n",
      "u_t = (0.010872 +0.000000i)u_1\n",
      "    + (-0.993135 +0.000000i)u_11\n",
      "    + (-0.994363 +0.000000i)u_1111\n",
      "    + (-0.998453 +0.000000i)uu_1\n",
      "    + (0.000462 +0.000000i)uuu_11\n",
      "    + (0.000465 +0.000000i)uuu_1111\n",
      "   \n",
      "u_t = (0.013149 +0.000000i)u_1\n",
      "    + (-0.993254 +0.000000i)u_11\n",
      "    + (-0.994077 +0.000000i)u_1111\n",
      "    + (-0.998535 +0.000000i)uu_1\n",
      "    + (0.003518 +0.000000i)uuu_1\n",
      "    + (0.000717 +0.000000i)uuu_11\n",
      "    + (0.000720 +0.000000i)uuu_1111\n",
      "   \n",
      "u_t = (0.002925 +0.000000i)u\n",
      "    + (0.012903 +0.000000i)u_1\n",
      "    + (-0.994150 +0.000000i)u_11\n",
      "    + (-0.994419 +0.000000i)u_1111\n",
      "    + (-0.998819 +0.000000i)uu_1\n",
      "    + (0.003328 +0.000000i)uuu_1\n",
      "    + (0.000953 +0.000000i)uuu_11\n",
      "    + (0.000701 +0.000000i)uuu_1111\n",
      "   \n"
     ]
    }
   ],
   "source": [
    "all_ols_models = []\n",
    "all_eff_rhs = []\n",
    "n_ensemble_models = 50; miosr_alpha = 5e-2\n",
    "for mnz in range(1, 9):\n",
    "    optimizer = ps.MIOSR(target_sparsity=mnz, normalize_columns=False, alpha=miosr_alpha, regression_timeout=100)    \n",
    "    model = ps.SINDy(feature_library=weak_pde_lib, optimizer=optimizer, \n",
    "                     differentiation_method=differentiation_method(**differentiation_kwargs), \n",
    "                     cache=True,\n",
    "                     feature_names=['u'])\n",
    "    model.fit(np.expand_dims(un, -1), t=dt, ensemble=True, \n",
    "              library_ensemble=True, n_candidates_to_drop=1, n_models=n_ensemble_models)\n",
    "    ols_model, eff_rhs_indices = model.validate_coef_list()[0]\n",
    "    all_ols_models.append(ols_model)\n",
    "    all_eff_rhs.append(eff_rhs_indices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "617a9a57",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_eff_rhs.append(tuple(i for i in range(X_pre.shape[1])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "93a0192c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RFE_PLinearRegression\n",
      "(3, 6)\n",
      "(6,)\n",
      "--------------------------------------------------\n",
      "(3, 5, 6)\n",
      "(3, 6)\n",
      "--------------------------------------------------\n",
      "(3, 4, 5, 6)\n",
      "(3, 5, 6)\n",
      "--------------------------------------------------\n",
      "(0, 3, 5, 6, 9)\n",
      "(0, 3, 5, 6)\n",
      "(3, 5, 6)\n",
      "--------------------------------------------------\n",
      "(2, 3, 5, 6, 9, 13)\n",
      "(3, 5, 6, 9, 13)\n",
      "(3, 5, 6, 13)\n",
      "(3, 5, 6)\n",
      "--------------------------------------------------\n",
      "(2, 3, 5, 6, 7, 9, 13)\n",
      "(2, 3, 5, 6, 9, 13)\n",
      "--------------------------------------------------\n",
      "(0, 2, 3, 5, 6, 7, 9, 13)\n",
      "(2, 3, 5, 6, 7, 9, 13)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13)\n",
      "(0, 1, 2, 3, 4, 5, 6, 7, 9, 10, 11, 12, 13)\n",
      "(0, 1, 2, 3, 4, 5, 6, 7, 10, 11, 12, 13)\n",
      "(0, 2, 3, 4, 5, 6, 7, 10, 11, 12, 13)\n",
      "(0, 2, 3, 5, 6, 7, 10, 11, 12, 13)\n",
      "(0, 2, 3, 5, 6, 7, 11, 12, 13)\n",
      "(0, 2, 3, 5, 6, 7, 11, 13)\n",
      "(2, 3, 5, 6, 7, 11, 13)\n",
      "(2, 3, 5, 6, 11, 13)\n",
      "(3, 5, 6, 11, 13)\n",
      "(3, 5, 6, 13)\n",
      "--------------------------------------------------\n",
      "RFE_PLinearRegression\n",
      "(3, 6)\n",
      "(6,)\n",
      "--------------------------------------------------\n",
      "(3, 5, 6)\n",
      "(3, 6)\n",
      "--------------------------------------------------\n",
      "(3, 4, 5, 6)\n",
      "(3, 5, 6)\n",
      "--------------------------------------------------\n",
      "(0, 3, 5, 6, 9)\n",
      "(0, 3, 5, 6)\n",
      "(3, 5, 6)\n",
      "--------------------------------------------------\n",
      "(2, 3, 5, 6, 9, 13)\n",
      "(3, 5, 6, 9, 13)\n",
      "(3, 5, 6, 13)\n",
      "(3, 5, 6)\n",
      "--------------------------------------------------\n",
      "(2, 3, 5, 6, 7, 9, 13)\n",
      "(2, 3, 5, 6, 9, 13)\n",
      "--------------------------------------------------\n",
      "(0, 2, 3, 5, 6, 7, 9, 13)\n",
      "(2, 3, 5, 6, 7, 9, 13)\n",
      "--------------------------------------------------\n",
      "(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13)\n",
      "(0, 1, 2, 3, 4, 5, 6, 7, 9, 10, 11, 12, 13)\n",
      "(0, 1, 2, 3, 4, 5, 6, 7, 10, 11, 12, 13)\n",
      "(0, 2, 3, 4, 5, 6, 7, 10, 11, 12, 13)\n",
      "(0, 2, 3, 5, 6, 7, 10, 11, 12, 13)\n",
      "(0, 2, 3, 5, 6, 7, 11, 12, 13)\n",
      "(0, 2, 3, 5, 6, 7, 11, 13)\n",
      "(2, 3, 5, 6, 7, 11, 13)\n",
      "(2, 3, 5, 6, 11, 13)\n",
      "(3, 5, 6, 11, 13)\n",
      "(3, 5, 6, 13)\n",
      "--------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "st = backward_refinement(all_eff_rhs, (X_pre, y_pre), mode='rfe', verbose=True)\n",
    "st += backward_refinement(all_eff_rhs, (normalize(X_pre, axis=0), y_pre), mode='rfe', verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "b649d56e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(6,)\n",
      "(3, 6)\n",
      "(3, 5, 6)\n",
      "(3, 4, 5, 6)\n",
      "(0, 3, 5, 6, 9)\n",
      "(2, 3, 5, 6, 9, 13)\n",
      "(2, 3, 5, 6, 7, 9, 13)\n",
      "(0, 2, 3, 5, 6, 7, 9, 13)\n",
      "(0, 2, 3, 5, 6, 7, 11, 12, 13)\n",
      "(0, 2, 3, 5, 6, 7, 10, 11, 12, 13)\n",
      "(0, 2, 3, 4, 5, 6, 7, 10, 11, 12, 13)\n",
      "(0, 1, 2, 3, 4, 5, 6, 7, 10, 11, 12, 13)\n",
      "(0, 1, 2, 3, 4, 5, 6, 7, 9, 10, 11, 12, 13)\n",
      "(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13)\n"
     ]
    }
   ],
   "source": [
    "decreasing_aics = []\n",
    "complexities = []\n",
    "eff_feature_names = []\n",
    "now_aic = 1e6\n",
    "for com in sorted(st.track.keys()):\n",
    "    indices, aic = st.track[com]\n",
    "    print(indices)\n",
    "    if aic < now_aic:\n",
    "        complexities.append(len(indices))\n",
    "        decreasing_aics.append(aic)\n",
    "        eff_feature_names.append([convert2latex(f) for i,f in enumerate(model.get_feature_names()) if i in indices])\n",
    "        now_aic = aic\n",
    "decreasing_aics = np.array(decreasing_aics)\n",
    "complexities = np.array(complexities)\n",
    "relative_decreasing_aics = decreasing_aics-np.min(decreasing_aics)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "cb26bcbe",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxYAAAGICAYAAADcTXa9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAADXZUlEQVR4nOzdd3hM2RvA8e9MmkQighDRWdHLrt57z+rdYrEsVl/8Fmv1vnpZq7N6r6uzevTeWZ2oaUif8/tjdm4yBAlhUt7P8+SRe++ZmXeOTHLee5pOKaUQQgghhBBCiE+gt3QAQgghhBBCiPhPEgshhBBCCCHEJ5PEQgghhBBCCPHJJLEQQgghhBBCfDJJLIQQQgghhBCfTBILIYQQQgghxCeTxEIIIYQQQgjxySSxEEIIIYQQQnwySSyEEEIIIYQQn0wSCyGEEEIIIcQnk8RCCCGEEEII8ckksRAilhkMBkuHIIQQQgjxxUliIcQnuH//Pq9fv+b+/ftcv34dAL1ePlZCCCGESHysLR2AEPHV0aNHGTJkCIcPHyZ//vy8fv2aFy9eUKdOHb777jsKFSoEgFIKAJ1OZ8lwhRBCCCE+K50ytXqEEDFy9uxZnj59Ss6cOTl58iQGg4FLly6xatUqzp8/T40aNRg7diy5c+e2dKhCCCGEEJ+dJBZCxDKlFPv27WPkyJE8evSIgQMH0rhxY0uHJYQQQgjxWclQKCE+klIKnU731lAnnU5H+fLlyZw5M7Nnz2bYsGE4ODjg6elpyXCFEEIIIT4r6bEQIpaYEo3w8HCsrKy082PHjmXbtm3MmzePzJkzWy5AIYQQQojPSJavESKWmHosFi5cSHBwsHb+559/pn79+ty5c8dSoQkhPjMPDw8KFCjwzuv+/v64urrSoEEDs/O3b9/GysqK169fA/D8+XO6du1KxowZcXJyIkeOHIwaNYrw8PDPGr8QQsQGSSyEiEUvXrxg48aNbN68GUDrvejSpQtFixa1cHRCiM/h6NGj6HQ6fH19OX/+fJRlRo0aRVhY2FvXz549S/bs2XFwcODp06eUKFGCoKAgjhw5QkBAAKtXr+bw4cOyqpwQIl6QoVBCxBLTUKh9+/Zx8OBB+vfvrzU2kidPbunwhBCfSdeuXXFxcSEwMBC9Xs+YMWPMrj948AAPDw/+/PNPWrduTUBAAA4ODgAMHTqUixcvsmLFCm2Rh5UrV37x9yCEELFBeiyE+AQGg0EbohAcHMzixYtZuHAhs2bNonbt2rRq1YqffvqJ9evXWyS+BQsWoNPp0Ol0/PPPP29dV0rx1VdfaRPOY8vgwYM/+g6rKebbt2/HWjyWfs3Dhw8zePBgfH19v+jrvo/p/+jZs2df/LU/xaf8bH0OYWFhrFy5kiZNmtC0aVOWLl2KwWAwK/Prr79StWpVmjVrho2NDZcuXdKunT17lgIFCnDjxg1Wr17N0KFDv/RbEEKIWCOJhRAxZOrk8/b2Rq/XaxO1T506xaZNm2jQoAH169fH29ub+vXr06lTJ7799ltLhoyTkxNz58596/y+ffu4efMmTk5OFogq8Th8+DBDhgx5K7GoVasWR44cIW3atJYJTHyy7du3kypVKvLkyUOhQoVIkiQJ+/bt066fP3+eJUuWMHz4cKysrPDw8DAbDmVKLPbs2YOHhwc5c+a0xNsQQohYIYmFEDGk0+m4fv067u7u9OjRg7t37wJQsmRJVqxYQa1atZg4cSK5c+cmadKklC5d2myVKEto0qQJa9aswd/f3+z83LlzKVGiBBkzZrRQZImbq6srxYsXx87OztKhiI+0ePFimjRpoh03adKExYsXa8d9+/alSZMm5MmTB4BcuXJpicXLly/5999/KVCgAE+ePCF9+vRfNnghhIhlklgI8RGyZ8/OX3/9xaVLl+jevTvr1q3j5cuXANqKUD169OCvv/4iNDTUkqEC0KxZMwCWLVumnfPz82PNmjW0bds2ysccPHiQSpUq4eTkhIODAyVLlmTLli1vlduyZQsFCxbEzs6OLFmy8Pvvv78zjuvXr9O8eXNSp06NnZ0duXLlYvr06R/1np4+fUqHDh3IkCEDdnZ2uLq6UqpUKXbt2hUrrxmTx125coVmzZqRJk0a7OzsyJgxI61atSI4OJjBgwfTp08fALJkyWI2NO1dQ6GiU/emIUEXL16kWbNmODs7kyZNGtq2bYufn1+06/HevXvUr1+fZMmS4ezszHfffcfTp0/Nyty4cYM2bdpok4zTpUvHt99++9ZE5Nj+P4nJz1ZUPnc9vnz5ko0bN76VWKxZs4agoCD27NnDtm3b2LRpE25ubri5ubF582at3s6dO4eLiwvp06cnY8aM3Lt3L0bvTwgh4hpJLIT4CEopmjRpQv/+/bG3t6d///5MnDgRpZR29/mrr74iU6ZM2NjYWDhaSJYsGQ0bNmTevHnauWXLlqHX680aRSb79u2jYsWK+Pn5MXfuXJYtW4aTkxPffvstK1as0Mrt3r2bOnXq4OTkxPLlyxk3bhwrV65k/vz5bz3npUuXKFKkCBcuXGD8+PFs3ryZWrVq0a1bN4YMGRLj99SyZUvWr1/Pb7/9xo4dO5gzZw6VK1fm+fPnn/yaMXnc2bNnKVKkCF5eXgwdOpStW7cyatQogoODCQkJ4YcffqBr164ArF27liNHjnDkyBG++eabKF87unVv0qBBAzw8PFizZg2//PILS5cupWfPntGux3r16vHVV1+xevVqBg8ezPr166lWrZpZQvzw4UNSpkzJ6NGj2bZtG9OnT8fa2ppixYpx9epVrVxs/p/E5GfLUvW4Zs0asmbNSqZMmQgKCiIoKIjs2bOTNm1aNm3aRN++fenZsydXrlzhzJkznDlzhj///FNLLEzDoAA8PT3x9fVlxIgRvH79msDAQA4fPkz79u2j9X6FECJOUEKIT2IwGNS0adOUo6Ojyps3r5o9e7YyGAxvlbGE+fPnK0AdP35c7d27VwHqwoULSimlihQpor7//nullFJ58uRR5cqV0x5XvHhxlTp1ahUQEKCdCwsLU3nz5lXp06fX3k+xYsWUu7u7CgwM1Mr5+/urFClSqDd/vVSrVk2lT59e+fn5mZ3v0qWLSpIkiXrx4oVZzLdu3Xrve3N0dFQ9evR4b5mPfc3oPk4ppSpWrKiSJ0+unjx58s44xo0bF+V7iuq9RrfuBw0apAA1duxYs+fs3LmzSpIkyQd/5kyP79mzp9n5JUuWKEAtXrz4nY8NCwtTISEhKnv27GaPj83/k5j8bEXlS9RjlSpVFBDll7W1tcqSJYt69eqV2WNu376tAPXkyRP1448/mtXXuXPnVLVq1VTKlCmVs7OzKly4sJo9e/YH36sQQsQV0mMhRAzNnDmToKAg7Vin0/HTTz8xdOhQrl27xrFjx9DpdGYrw8SFVWzKlStHtmzZmDdvHufPn+f48eNRDoN69eoVR48epWHDhjg6OmrnraysaNmyJffv3+fq1au8evWK48ePU79+fZIkSaKVM90VjiwoKIjdu3dTr149HBwcCAsL075q1qxJUFAQXl5eMXo/RYsWZcGCBQwfPhwvL6+3hpx97GvG5HGvX79m3759NG7cGFdX1xjFH5Xo1n1ktWvXNjvOnz8/QUFBPHnyJFqv2aJFC7Pjxo0bY21tzd69e7VzYWFhjBw5kty5c2Nra4u1tTW2trZcv36dy5cva+Vi6/8kJj9bUflS9bhjxw6UUlF+hYaG8u+//2rLyppkypQJpRSurq7MnDmTiRMnatfy5cvHtm3bePbsGb6+vhw/fpwffvjhg+9XCCHiCkkshIiBO3fuMGDAAGrVqsWqVat49eqVdq1Fixa0a9cuRsNQviSdTkebNm1YvHgxM2fOxMPDgzJlyrxVzsfHB6VUlCsVubu7A8bdgX18fDAYDLi5ub1V7s1zz58/JywsjKlTp2JjY2P2VbNmTYAYL3u6YsUKWrduzZw5cyhRogQpUqSgVatWeHt7f9JrxuRxPj4+hIeHx9qk2+jWfWQpU6Y0OzYNxQsMDIzWa775f2VtbU3KlCnNXqdXr14MHDiQunXrsmnTJo4ePcrx48cpUKCA2evE1v9JTH62omKJehRCCAHWlg5AiPgkffr0bNiwgQULFjBt2jTu3LlDp06dSJo0KSEhIZw5c0ZbLlKvj3t5+/fff89vv/3GzJkzGTFiRJRlXFxc0Ov1PHr06K1rDx8+BCBVqlS4uLig0+m0RmNkb55zcXHR7hb/9NNPUb5ulixZYvReUqVKxaRJk5g0aRJ3795l48aN/PLLLzx58oRt27Z99GvG5HEpUqTAysqK+/fvxyj2d4lu3ccmb29v0qVLpx2HhYXx/Plzs4b24sWLadWqFSNHjjR77LNnz8w2f4yt/5MkSZJE+2crKpaoRyGEEJJYCPFB6r8dtcE4lKJ06dJky5aNWbNmMXz4cJYsWULbtm3ZsmUL2bNnR6fTER4ebvElZqOSLl06+vTpw5UrV2jdunWUZZImTUqxYsVYu3Ytv//+O/b29oBxM8DFixeTPn16PDw80Ol0FC1alLVr1zJu3DhtyEpAQACbNm0ye04HBwcqVKjA6dOnyZ8/P7a2trH6vjJmzEiXLl3YvXs3hw4d+qTXjMnj7O3tKVeuHKtWrWLEiBHvbKxG9+53dOs+Ni1ZsoRChQppxytXriQsLMxsw0SdTvfWkrhbtmzhwYMHfPXVV1E+76f+n0T3ZysqX6oeLTHEUf23j44QQsRFklgIEU2DBw/mypUrTJ8+nbRp0zJo0CC6devGgAEDGD16NJUqVWLUqFFA3JhT8S6jR4/+YJlRo0ZRpUoVKlSoQO/evbG1tWXGjBlcuHCBZcuWae9v2LBhVK9enSpVqvDzzz8THh7OmDFjSJo0KS9evDB7zsmTJ1O6dGnKlClDp06dyJw5MwEBAdy4cYNNmzaxZ8+eaL8HPz8/KlSoQPPmzcmZMydOTk4cP36cbdu2Ub9+/U9+zZg8bsKECZQuXZpixYrxyy+/8NVXX/H48WM2btzIn3/+iZOTE/ny5dOet3Xr1tjY2JAjR45PqvvYsnbtWqytralSpQoXL15k4MCBFChQgMaNG2tlPD09WbBgATlz5iR//vycPHmScePGmQ0Bi+3/k5j8bEXlS9RjdBv5Hh4e2Nvbc/bs2U9+TSGEiNMsMmVciHho06ZNqmDBgqpMmTJqzpw56unTp2bXI69eE1dEXhXqfd5cFUoppQ4cOKAqVqyokiZNquzt7VXx4sXVpk2b3nrsxo0bVf78+ZWtra3KmDGjGj16tLbSzptu3bql2rZtq9KlS6dsbGyUq6urKlmypBo+fPhbMb9vVaigoCDVsWNHlT9/fpUsWTJlb2+vcuTIoQYNGvTWKjwf+5rReZzJpUuXVKNGjVTKlCm1evj+++9VUFCQVqZfv37K3d1d6fV6Bai9e/e+871Gp+5Ndfzmz2F0V9UyPf7kyZPq22+/VY6OjsrJyUk1a9ZMPX782Kysj4+PateunUqdOrVycHBQpUuXVgcOHFDlypXTfm5i+/9EqZj9bEXlS9Tjh3h5eSkPDw+VMWNGde7cuU96LiGEiOt0Skm/qhDRdebMGebNm8epU6eoWLEiPXv2xMXFxdJhCSHiqK5du+Li4kJgYCB6vZ4xY8Zo165cuUKXLl04deoUSim+++47pk6dasFohRDi00hiIUQ0qf/mWrx8+ZK//vqLwYMHkyJFCkaMGGE21EMIIcA4ET5dunTs2bOHoKAg6taty507d7SFHQoVKsTPP/9Ms2bNCAgI4Nq1axQuXNjCUQshxMeLe8vWCBFHmPah8Pf3ByLmTTg6OtKpUydOnDjB06dPadiwIUuXLrVYnEKIuGn79u2kSpWKPHnyUKhQIZIkScK+ffu06//++y/h4eEYDAaSJUsmSYUQIt6TxEKIdzDdVRwxYgQuLi4sXbqU169fA8akI3369PTo0YP+/ftrPRbSASiEMFm8eDFNmjTRjps0acLixYu142XLljFr1ixttbY3NxQUQoj4RoZCCfEB165dY9q0aVy4cIEyZcrQtm1bMmXKBMCsWbPw9vbmt99+M1uWVgiRuL18+ZI0adJw6tQpbQWw8+fPU6ZMGby9vc12FL9z5w5ly5bljz/+0DYKFEKI+Eh6LIT4AA8PD/r370/NmjVZvHgx9evXZ/v27Tx8+JAhQ4Zou1dLji6EMFmzZg1Zs2YlU6ZMBAUFERQURPbs2UmbNi2bNm1i7dq13Lp1CzDuFB4cHEyuXLksHLUQQnwa6bEQIgaePXvGL7/8wrx580idOjXFihVjw4YNlg5LCBHHVK1alZ07d0Z5rXbt2mTOnJkVK1bw8uVLsmbNyvDhw6ldu/YXjlIIIWKXJBZCRINSCqWUNu/i8ePH3Llzh6xZs5IqVSoMBoN2TQghhBAiMZLEQgghhBBCCPHJ5BarEEIIIYQQ4pNZWzqAhMRgMPDw4UOcnJxkdaAEzmAwoJTCysrK0qEIIYQQX5RSioCAANzd3WUYsDAjiUUsevjwIRkyZLB0GEIIIYQQn929e/dInz69pcMQcYgkFrHIyckJMH7QkiVLFrMHHz0KVavCzp1QtGisxPPixQuKFCnCnj17tH0XPqeWLVtSrFgxunTp8tlfy5KqV69O0qRJmT17NilSpLB0OEIIIcQX5e/vT4YMGbR2jxAmkljEItPwp2TJksU8sdi8OeLfypVjJZ6hQ4dSu3Zt8uXLFyvP9yHDhg2jQoUKdOnSJebvPw4ybXj39OlTUqRIgU6nQylFq1at2LRpE3v27KFt27aWDlMIIYSwCBn2Ld4kA+PiAoMBVq40fr9yJcTCQl2BgYHMnTuXH3744ZOfK7ry589P5syZWbJkyRd7zc/J9Atz8uTJXL9+Hb1ej5WVFR06dKBNmzZcvXrVwhEKIYQQQsQdkljEBcePg7e38ftHj4zHn2jr1q1YW1tTokQJs/Pbtm3D3t6esLAw7dzly5fR6XQ8e/bsg9c/pHbt2ixbtuyT448LlFIYDAaSJEnCokWLtPN6vZ6GDRsyZswYrZwQQgghRGIniUVcsHo1WP83Ks3a2nj8ifbv30/hwoXfOn/mzBny5MmDtbW12bl06dKRKlWqD17/kKJFi3Ls2DGCg4M/+T1YkmkYlF6vp3r16jx//pwtW7bw5MkTTp8+zYEDB7Sy0hUshBBCCCFzLCxPKVixAkw9BGFhxuMxY+ATGqy3b9/G3d39rfNnz56lYMGCZudOnz5NgQIFonX9Q9KlS0dwcDDe3t5fZML456LT6QgICKBNmzY4ODhw6NAh1q5dS6ZMmciWLRvfffcdEJGACCHiJqUUYWFhhIeHWzoUIeINKysrrK2t5e+biDFJLCztzBm4d8/83N27cPYsvNHAj4nAwECSJEkSxcudoXPnzm+dM/VufOj6h9jb2wPw+vXrjwk7TkmSJAkNGzakUqVKPH78mFGjRjF58mR0Oh0pU6YEpLdCiLgsJCSER48eJYjfR0J8aQ4ODqRNmxZbW1tLhyLiEUksLG3NGrCygsh306ysjOc/IbFIlSoVPj4+ZucCAwO5fv26WY+EwWDg1KlTtGvX7oPXTYoVK8b06dMpXLgwrVu3pnjx4nTq1AkwLnEL4Orq+tGxxxV79uyhadOmgPH95M2blyFDhjB16lQLRyaE+BCDwcCtW7ewsrLC3d0dW1tbuREgRDQopQgJCeHp06fcunWL7NmzyyZ4CZDBYPgs/6+SWHxOU6fC4MHvX+Xp5UvzpAKMx2PGwPTp736cTgdDhsA79oz4+uuvWbx4sdm5mzdvEh4eTo4cObRz27dv5/nz5xQoUOCD100GDhzIyJEjKVWqFI6OjlpSAXDhwgXSp08frfkYcdn169eZOHEiWbNmJXv27AD069ePkiVL8ujRI9KmTWvhCIUQ7xMSEoLBYCBDhgw4ODhYOhwh4hV7e3tsbGy4c+cOISEhUY6AEHGfKXl4/vw5gYGBJEuWjFu3bpErVy6znqjYHNYtKejnVLIk2NiAj8+7v0JDAXjKGw3x0ND3P87WFt5Y8SmyatWqcfHiRbNei5QpU6LT6Th27BgAXl5edOnSBXt7e7Jnz/7B6yaenp78+++/bN++ncmTJ5u97oEDB6hateonVZslRF7ZKTw8nOzZs1O/fn2OHDlidn7x4sWSVAgRj8idViE+jnx24j/T/2Hv3r0pWbIkxYsXp0ePHmTJkgVPT09WrVpFcHCwllTExiqX0mPxORUqBBcvQps2sGnTO4sFkoTsXMeDa7RmIU1ZTkpevPt5a9eGefPgv3H+UcmXLx+FCxdm5cqV/PjjjwCkTZuWYcOG0apVKxwdHSlfvjyNGjVi9+7dWFlZffC6ybFjx/D19cXDw8Ns9aigoCDWrVvH9u3bY1BJlhMeHo6VldVbmfrLly+5evUqJ0+e5PDhw1y+fJkUKVJw9epVunXrZsGIhRBCCCGi79WrV1SrVo3Jkydz5swZfHx88PX1ZcOGDbRu3ZrUqVMzbNgwWrZsGTu9FkrEGj8/PwUoPz8/8wsGg1LTpytlY6OUlZVSxsFR2tcKGpmdsiFY1WONWk9tFYyN8aSVlfHxM2YYny8atmzZonLlyqXCw8Nj7T3ev39f5cuXT/3777+qYMGC6tKlS9q1adOmqSpVqsTaa30pBQsWVKNGjVJhYWFKKaUWLVqkPD091YQJE1TRokVVkSJF1G+//abWr19v4UiFENEVGBioLl26pAIDAy0dSrzTunVrVadOHUuH8ckyZcqkJk6caOkwYhWg1q1bF+3yn/J/+b7P0DvbOyJeuXv3rurZs6fKli2b6tWrl3r9+vUnP6f0c30JOh107gwnT0K2bPBG9+JCWpsdh2LLOupTlw2k4wHdmMKJ9HVRJ05Cp07RXoa2Zs2a/Pjjjzx48CBW3kZgYCANGzZk2rRpZMmShb59+zJ8+HDtuo2NTbyb2GwwGGjSpAkrVqygcePGbN++nZYtW7JixQp69uzJ/v37yZUrFz169KBOnTqWDlcIkUgcPnwYKysrqlev/ta127dvo9PpOHPmjNn5NWvWUL58eZydnXF0dCR//vwMHTpUW1QjsTl+/DgdOnSwdBhCWJSKNLxJvTHUKUOGDAwdOpShQ4eyb98+BgwYgMFg+KTXk8TiS8qXD06fhh9+MB7/lyCM4X/0ZhxuPHrrIc9wZSpdKXJnNXmb5WPMGIhJntC9e3cyZMgQG9Fjb2/PkSNHKFu2LADNmjVjyZIl2vUOHTqYTfyOD/R6PX369GHkyJFYW1szfPhwZs+erU32tLOzI1++fCxYsMCygQohEpV58+bRtWtXDh48yN27dz9YfsCAATRp0oQiRYqwdetWLly4wPjx4zl79ix//fXXF4g4+kJCQr7I67i6usrEfZHoRR7eFNVQJ0dHR5o3b86MGTPw8vJi1apVn/R6klh8aQ4O8OefkDu3diovFxlHX+6Rga1UpynLSELgWw+9dAl++QUyZICqVWHJEnj16ksGnzBZWVlRo0YNJk2aRIkSJejSpQs1atRg//79ANSoUYMrV65YOEohRGLx6tUrVq5cSadOnfD09PzgjY1jx44xcuRIxo8fz7hx4yhZsiSZM2emSpUqrFmzhtatW7/zseHh4fTq1YvkyZOTMmVK+vbt+9ZdTaUUY8eOJWvWrNjb21OgQAFWr15tVubixYvUqlWLZMmS4eTkRJkyZbh58yYA33//PXXr1mXUqFG4u7vj4eEBwIMHD2jSpAkuLi6kTJmSOnXqcPv2be05jx8/TpUqVUiVKhXOzs6UK1eOU6dOmb3u4MGDyZgxI3Z2dri7u5vNg8ucOTOTJk3SjnU6HXPmzKFevXo4ODiQPXt2Nm7caPZ8GzduJHv27Njb21OhQgUWLlyITqfD19f3nXWo0+n4888/8fT0xMHBgVy5cnHkyBFu3LhB+fLlSZo0KSVKlNDqw+SPP/4gW7Zs2NrakiNHjrcSwOvXr1O2bFmSJElC7ty52blz51uv/aE6fNPq1avJly8f9vb2pEyZksqVK/NKGhKJjqlX4saNGzx8+BCAokWL8ttvv+Hn5/dpvRafPJhKaKI95vD27bfmWbz55UsyNZt2qgz73lvU0VGpNm2U2rtXqVicSpGghYaGKqWUunr1apTXV6xYoXQ6ncqaNat6/Pix2bXYnK8ihPh83hwf7uur1IEDlvvy9Y1+7HPnzlWFCxdWSim1adMmlTlzZmWINLfu1q1bClCnT59WSinVrVs35ejoqEJCQmJcT2PGjFHOzs5q9erV6tKlS6pdu3bKycnJbFx+//79Vc6cOdW2bdvUzZs31fz585WdnZ36559/lFLGuXcpUqRQ9evXV8ePH1dXr15V8+bNU1euXFFKGcf5Ozo6qpYtW6oLFy6o8+fPq1evXqns2bOrtm3bqnPnzqlLly6p5s2bqxw5cqjg4GCllFK7d+9Wf/31l7p06ZIWW5o0aZS/v79SSqlVq1apZMmSqb///lvduXNHHT16VM2aNUuL+805FoBKnz69Wrp0qbp+/bpWb8+fP9fq1cbGRvXu3VtduXJFLVu2TKVLl04BysfH5511CKh06dKpFStWqKtXr6q6deuqzJkzq4oVK6pt27apS5cuqeLFi6vq1atrj1m7dq2ysbFR06dPV1evXlXjx49XVlZWas+ePUop49+avHnzqvLly6vTp0+rffv2qa+//tpsjkV06jDyHIuHDx8qa2trNWHCBHXr1i117tw5NX36dBUQEBDl+5I5Fglf37591fz585VSSptn+vLly096TkksYlG0P2gTJiil138wuVCglF6vbgyYpwYNUipLlvcXzZxZqYEDlbp27Yu83XgvY8aMqkyZMmr37t1mCcPTp0/VDz/8oFatWqWUikhEhBDxx5uNogMHovcr93N9HTgQ/dhLliypJk2apJQy/v5JlSqV2rlzp3b9zcSiRo0aKn/+/B9VT2nTplWjR4/WjkNDQ1X69Om1xujLly9VkiRJ1OHDh80e165dO9WsWTOllFL9+vVTWbJkeWdi07p1a5UmTRqtsauUMXnKkSOHWcIUHBys7O3t1fbt26N8nrCwMOXk5KQ2bdqklFJq/PjxysPD452vG1Vi8euvv2rHL1++VDqdTm3dulUppdT//vc/lTdvXrPnGDBgQLQSi8jPe+TIEQWouXPnaueWLVumkiRJoh2XLFlStW/f3ux5GjVqpGrWrKmUUmr79u3KyspK3bt3T7u+detWs8QiOnUYObE4efKkAtTt27ff+V4ik8Qi4TL9zBw4cEB169YtVp9bhkJZwooVb2+aZ1rONdKyrgAoRbbdsxg8GG7ehP37jVM0kiV7+2lv34Zhw8DDA0qVMo64ek/vbaKmlGL+/PlkyJCBMWPGMG/ePK07OFWqVPj5+ZExY0YAs6V2hRDic7p69SrHjh2jadOmAFhbW9OkSRPmzZv3zseoaGxudffuXRwdHbWvkSNH4ufnx6NHjygRaU8ka2trChcurB1funSJoKAgqlSpYvb4RYsWaUN7zpw5Q5kyZbCxsXnn6+fLl89sQ66TJ09y48YNnJyctOdMkSIFQUFB2vM+efKEjh074uHhgbOzM87Ozrx8+VKbc9KoUSMCAwPJmjUr7du3Z926dYSFhb23HvLnz699nzRpUpycnHjy5AlgrPsiRYqYlS9atOh7ny+q502TJo32niOfCwoKwt/fH4DLly9TqlQps+coVaoUly9f1q5nzJiR9OnTa9dLvLF3VXTqMLICBQpQqVIl8uXLR6NGjZg9e7bZXlcicTAYDOh0Ovz8/Hj9+jW3b9+mR48ebNiwgdWrVzNhwoRPen7Zx+JLe/gQjh59+7ynJ/z+O/TuDRs2RJxXylj+0SN0adNSpgyUKQNTphiLLVwIO3bAm8PhDh82fnXvbtz2onVr47yM9/zeT1R0Oh0VK1bE3d2dGTNm8Msvv7Bp0yZ+++03rl27xv79+1m5cqVWVgghvoS5c+cSFhZGunTptHNKKWxsbPDx8cHFxeWtx3h4eHDw4EFCQ0Pf2bh3d3c3W0UqRYoU0YrHNNZ6y5YtZjGBcXELMC7s8SFJkyZ963kLFSpktgCIiaurK2Ccm/H06VMmTZpEpkyZsLOzo0SJEtrk7wwZMnD16lV27tzJrl276Ny5M+PGjWPfvn3vrIc3z+t0Ou09RpWgqWhuGBb5eU3PEdW5yGPXo3qt921U9mb56NRhZFZWVuzcuZPDhw+zY8cOpk6dyoABAzh69ChZsmT54HsUCYNeryc4OJiqVatSqlQpwsLCmDZtGnfv3iVdunQ0btwYiNi1O6YksfjS1q0zrgallLF3wsrKmCV06GA8v26dsauhe3djtmC6+7JunXHJ2v/Y20PTpsavR4+ME7kXLoQLF8xfLjgYVq0yfqVODS1aQKtWULDgl3vLccHs2bM5deoUU6dOxdramrCwMKytrcmZMydTpkyhffv29OrViyJFipA7d24GDhwIoJUTQsRv+fLBgQOWff0PCQsLY9GiRYwfP56qVauaXWvQoAFLliyhS5cubz2uefPmTJkyhRkzZtC9e/e3rvv6+pI8eXK++uqrt66lTZsWLy8vbbW/sLAwTp48yTfffANA7ty5sbOz4+7du5QrVy7KuPPnz8/ChQvfm9i86ZtvvmHFihWkTp2aZFF1wQMHDhxgxowZ1KxZE4B79+7x7NkzszL29vbUrl2b2rVr89NPP5EzZ07Onz+vxR8TOXPm5O+//zY7d+LEiRg/T3TkypWLgwcP0qpVK+3c4cOHyZUrF2Cs97t37/Lw4UPc3d0BOHLkiNlzRKcO36TT6ShVqhSlSpXit99+I1OmTKxbt45evXrF0jsTcZkpebWzs2PDhg24uLhgZ2fH8OHDqVChAiVLltQS2I/eeT1WB1YlctEac1i2rDZ3QuXMqdSFC1GXO39eqRw5IuZilCv3wdc3GJQ6dUqpHj2UcnV9/3jf/PmV+v13pR49+rj3Gt+cP39eVa1aVS1evFg7ZzAYVEhIiDa/wmAwqIcPH6qrV6+anRNCxD/xcYO8devWKVtbW+UbxUzv/v37q4IFCyql3p5joZRxEqaVlZXq06ePOnz4sLp9+7batWuXatiwoTZfIyqjR49WLi4uau3atery5cuqffv2b03eHjBggEqZMqVasGCBunHjhjp16pSaNm2aWrBggVJKqWfPnqmUKVNqk7evXbumFi1aZDZ5+81N2kwTj8uXL6/279+v/v33X/XPP/+obt26afMKChYsqKpUqaIuXbqkvLy8VJkyZZS9vb02b2L+/Plqzpw56vz58+rmzZtqwIAByt7eXj179kwpFfUcizc3l3N2dtYmr/7777/KxsZG9e3bV129elWtWLFCpU+fXgFR/p+863mj+v/Zu3ev2VyNdevWKRsbG/XHH3+oa9euaZO39+7dq5QyTt7OnTu3qlSpkjpz5ozav3+/KlSoUJSTt99Xh5Hr3svLS40YMUIdP35c3blzR61cuVLZ2tqqv//+O8r3JXMsEq4JEyaYHR86dEjVrFkzVuaUSmIRiz74QXvyJCJR6NRJqQ/tcPj6tVIdO0YkIk+fRjuWkBClNm5UqmFDpWxt351gWFkpVbOmUsuXKxWP/v5+lD///FM5OzurH374Qd25c8fsWlhYmKz4JEQCEh8TC09PT23y7ptME29PnjwZZcNVKeOKdmXLllVOTk4qadKkKn/+/Gro0KHvnXgcGhqqunfvrpIlS6aSJ0+uevXqpVq1amWWCBgMBjV58mSVI0cOZWNjo1xdXVW1atXUvn37tDJnz55VVatWVQ4ODsrJyUmVKVNG3bx5Uyn17t2fHz16pFq1aqVSpUql7OzsVNasWVX79u21v6GnTp1ShQsXVnZ2dip79uxq1apVZsnCunXrVLFixVSyZMlU0qRJVfHixdWuXbu0549pYqGUUhs2bFBfffWVsrOzU+XLl1d//PGHAt77c/QxiYVSSs2YMUNlzZpV2djYKA8PD7Vo0SKz57169aoqXbq0srW1VR4eHmrbtm1vvdaH6jBy3V+6dElVq1ZNubq6Kjs7O+Xh4aGmTp36zvcliUXC9OzZM1WvXj21dOlSs/MNGjRQ586d++Tn1ykVzQGE4oP8/f1xdnbGz88v6m7JW7eMEx6GDoV69aL/xOvWwW+/wcaN8BHjIF+8gJUrjUOlvLzeXc7ZGRo3Ns7HKFky2ht8xyuHDx+mf//+JEmShOrVq1OlShXy5Mlj6bCEELEsKCiIW7dukSVLFpIkSWLpcEQ8NWLECGbOnMm9e/csHcoX977P0AfbO7FMRWOBAhF9+/fvZ+/evfTq1QsnJyeeP3+Oo6OjNm/qU0hiEYu+9AftY1y7BosWwV9/wfs2c82WzTgXo2XLj8pl4hxl7J1Dr9fzzz//sGDBAi5duqRtKpQ9e3aqVKnCv//+i52dHWXKlNHGtQoh4h9JLMTHmDFjBkWKFCFlypQcOnSIrl270qVLF4YPH27p0L44SyUW586d48GDB6RMmZI7d+6QL18+cubMqV2XJOPjvXjxggULFnD06FEOHDhAjhw5SJo0Kba2tnTo0IHq1at/8mtIYhGL4kNiYWIwwL59xiRj1ar37+BdtqwxyWjUKOplbuMjg8HAjh07+Oeff/Dx8eHixYvcv3+fYsWK4ejoSOfOnSlUqJClwxRCfCRJLMTH6NmzJytWrODFixdkzJiRli1b0q9fv0S5iIclEou7d+/yww8/cO/ePaytrXF1deXEiRO4urrSokULOnXqRNq0aWPt9RKb8+fP8+eff/Ldd9+xc+dO/vnnH0aMGIGDg4PZksmfQhKLWBSfEovIXr2CtWuNScbu3W9vsWGSJIlxBFfr1lC58ttbbsQXby6h9vLlSxwdHQEICAjAycnJUqEJIWKJJBZCfJovnViYRhZcunSJvHnz8ujRI548ecLr16/Ztm0bixYt4vHjx/Tu3Zuff/4ZZ2dn6b2IoSdPnpA6dWrtuGnTprRo0YJvv/021l5DEotYFF8Ti8ju3YtYuvbKlXeXS5sWvvvOmGTE1ykKkRMM+eUkRMIiiYUQn8aScyyi+pv86tUrFixYwKJFiyhZsiRjx46N9vLGAry8vJg9ezZz584lODgYOzs7rly5wq+//srSpUvNNrD8FHFm5+1Ro0ah0+no0aOHdk4pxeDBg3F3d8fe3p7y5ctz8eJFs8cFBwfTtWtXUqVKRdKkSalduzb37983K+Pj40PLli21nTtbtmyJ7xtbUt+9e5dvv/2WpEmTkipVKrp166ZtwpOYZMgAv/wCly7BsWPw008Q1T5Kjx7BuHGQNy8UKmTciuPp0y8f76eI3GshSYUQQggRN0T1Nzlp0qS0bduW/v37s3XrVlq2bImfn58Foot/lFJ4eHhgZ2eHl5eXNkk7W7ZseHp6YmtrG+3NID8kTiQWx48fZ9asWW+N7xo7diwTJkxg2rRpHD9+HDc3N6pUqUJAQIBWpkePHqxbt47ly5dz8OBBXr58iaenJ+Hh4VqZ5s2bc+bMGbZt28a2bds4c+YMLVu21K6Hh4dTq1YtXr16xcGDB1m+fDlr1qzh559//vxvPo7S6aBIEZg2zZhErF0LdetCVMNMT50y7ufn7g516sCaNcaN+b6o981EF0IkWtIpL8THiYufHXt7e+rUqcPevXsJCQnh8OHDlg4pTjP9H4aHh2NlZYWLiwuTJk3i77//5uzZs0ybNo26desCsXiD9ZMXrP1EAQEBKnv27Grnzp2qXLlyqnv37kop47rZbm5uavTo0VrZoKAg5ezsrGbOnKmUUsrX11fZ2Nio5cuXa2UePHig9Hq92rZtm1LKuG4zoLy8vLQyR44cUYC2ec/ff/+t9Hq9evDggVZm2bJlys7OLkZrNCeGdZ2fPlVq6lSlChd+/wZ8Li7GrTq8vIwb931Whw4ZX/Tw4c/8QkKI+CIsLExdunRJ2yxNCBEzz549U5cuXVJhYWFvXfuS7R3THlM7duxQ3t7e2vl9+/aptWvXfvbXTwh+++031aBBA9W0aVPl4OCgChcurCpXrqzGjBkT669l8WUOfvrpJ2rVqkXlypXNlnO7desW3t7eVK1aVTtnZ2dHuXLlOHz4MD/++CMnT54kNDTUrIy7uzt58+bl8OHDVKtWjSNHjuDs7EyxYsW0MsWLF8fZ2ZnDhw+TI0cOjhw5Qt68ec2WF61WrRrBwcGcPHmSChUqfOZaiD9SpYIuXYxfly5FLF378KF5OR8f+OMP41eOHBFL12bI8BmCWrEi4t8SJT7DCwgh4hsrKyuSJ0/OkydPAHBwcJAhj0JEg1KK169f8+TJE5InT46VhVdq0ev1KKX466+/qFChAm3atEEpRdmyZS0aV3zSt29f7t+/T44cOdi5cydHjx7l119/JTQ0FIjdeaYWTSyWL1/OqVOnOH78+FvXvL29AUiTJo3Z+TRp0nDnzh2tjK2tLS4uLm+VMT3e29vbbAa8SerUqc3KvPk6Li4u2NraamWiEhwcTHCkMT/+/v7vLJsQ5c4No0fDiBGwZ49xwvfatRAYaF7u6lUYMAB+/RUqVDBO+K5fH/5biOnTGAzG3f/A+O/EiQlzZz8hRIy5ubkBaMmFECL6kidPrn2GLMm00Eq3bt3YuHEjYBy2E5uN4YTIVG8HDhzAwcGBQoUKoZSidOnSjB07lurVq1O4cGEgdueZWiyxuHfvHt27d2fHjh3vXbHjzTcbnR+kN8tEVf5jyrxp1KhRDBky5L2xJAZWVlClivErIABWrzb2ZPzzj3k5pYwJyJ490LkzNGhg7MmoUAH0Hzvb5/hxMCV/jx4Zj4sW/ZS3I4RIIHQ6HWnTpiV16tTanTkhxIfZ2NhYvKfCRK/X4+Xlxd69e9m0aRNXr14lZ86cBAYG4u7ubrboj4hgWqDmyJEjXLt2jTlz5gDGeSodO3ZkwYIFWmIRmyyWWJw8eZInT56YbUIWHh7O/v37mTZtGlevXgWMvQmRN0N58uSJ1rvg5uZGSEgIPj4+Zr0WT548oWTJklqZx48fv/X6T58+NXueo0ePml338fEhNDT0rZ6MyPr160evXr20Y39/fzJ8lrE+8YeTE7RpY/y6fRsWLzb2ZNy4YV7u1Stj8rFokXF4lGnp2hw5YviCq1cbZ5SHhRn/Xb1aEgshhBkrK6s400gSQsTM3bt36d27N82bN6do0aKsXbuW/PnzkyVLFmrUqAHIkvGRhYeHo9frtfro1q0bI0eO5OLFi+T5b3+AGjVqaLtsv7m316ey2D4WAQEB2pAmkzZt2pAzZ07+97//kSdPHtzd3enZsyd9+/YFICQkhNSpUzNmzBh+/PFH/Pz8cHV1ZfHixTRu3BiAR48ekT59ev7++2+qVavG5cuXyZ07N0ePHqXofw3Oo0ePUrx4ca5cuUKOHDnYunUrnp6e3L9/X0tiVqxYQevWrXny5Em012hOCPtYfA5KgZeXMcFYsQLeWOnXTNGixgSjadOol7l964kzZTJuvmGSMaMxo5FfMEIIIcRn8SXaO6ZkISwsDCsrK62h/Ouvv1KyZElq1qz5WV43PjHVkb+/P5s3b6Z58+Zm19evX8+VK1dYuXIlQUFBFC1aFCsrK2xsbJg5c+ZnCyrOiLwqlFJKjR49Wjk7O6u1a9eq8+fPq2bNmqm0adMqf39/rUzHjh1V+vTp1a5du9SpU6dUxYoVVYECBcxWMahevbrKnz+/OnLkiDpy5IjKly+f8vT01K6HhYWpvHnzqkqVKqlTp06pXbt2qfTp06suXbrEKP7EsCrUpwoMVGrlSqVq1VLKyurdq0rZ2ChVv75SGzYoFRLyjic7dSrqB58+/SXfkhBCCJGofKn2TkBAgOrRo4dSyrg6lMFgUGfPnlUNGzZUQUFBn/W145MdO3YonU6natasqXbu3KmtpNWiRQu1YMECNW/ePOXh4aFmzJihNm3apK2CaioXmyy+KtT79O3bl8DAQDp37oyPjw/FihVjx44dODk5aWUmTpyItbU1jRs3JjAwkEqVKrFgwQKzbu8lS5bQrVs3bfWo2rVrM23aNO26lZUVW7ZsoXPnzpQqVQp7e3uaN2/O77///uXebCKRJAk0amT8evwYli41Doc6c8a8XGiocSL42rXGlaiaNzfOx/jmm0idEWvWGCd4RNqzBCsr4/mCBb/QOxJCCCFEbFORVqdat24d9erVAyB//vy4urpaOLq4pUqVKhw6dIiZM2cyfvx4/v33X3744QcWL16sDXVKlSoVW7ZsoVOnTtrjYnMIlInFhkIlRDIU6uOdO2dMMBYvNiYc75JHf5lWdsv5znYV7q9vGDOQN9nYvH/JKZ0OhgwxrpkrhBBCiBj5ku2dlStX4uXlRYcOHciaNSsnT54kb968ODk5ydyK/5iSh2vXrjFz5kz++usvSpQoweTJk8mSJYtWrlGjRvz6668UKFDgs8UiiUUsksTi04WFwc6dsHDiC9bvdCCYqFcM0xNOZXbRmoXUZT0OBEZZLkpubrB5M0RaOEAIIYQQ0fMl2jv//PMPCxcu5MmTJxw+fJiMGTNiZWVF7ty5GT58OJkzZ/4sr5sQnDx5knr16vHgwQP69OnD6NGjAfjjjz+oVKkSHh4eny0pi/0+ECE+gbU11KgBy3ekwPvma2bln0YpDr5VzoAVO6hGC5bihjc/8zthRGPVl9q14cIFSSqEEEKIOCxTpkxUrFiRxYsXM2XKFEqXLs2+ffuYNGmSJBX/uX37NuXLl+fkyZPaOaUUhQoVokuXLhQpUsRsg+hOnTrh4eEBxO7eFZFJYiHirORZU9D+zE8cnH6OG9Y5+U03jMzceqtcAMmYwM8s4Puon8jKyjg8asYMWL8eUqb8rHELIYQQIuZMg2h8fHx4/fo1LVu2xMXFhZYtW3Ljxg0OHDhAqlSpLByl5RgMBrNj0ybNP//8M5MmTcLb21tLGFq1akXevHmpUqVKlI/9XCSxEHGbTgedO5Pt1CqGZF/MTV129lGWdszBCfOdzpfR7O3H6/Xw1Vdw8iR06iTL0AohhBBxlKlRvH37dsaMGQNENJ5Hjx7NrFmzLBZbXGCabB3+36I1OXLkYObMmRQtWpQFCxbQr18/Hjx4AMCqVat48eIFjo6OKKU+y0TtKGP8Iq8ixKfKlw9On0bfvh1lOcAcXQe8caMpy7Qi/1Cep/x3J8OUQLRvD6dOGR8vhBBCiDivfv36pEuXjgsXLmBnZwdA3rx5GTp0KBDRsE4sdu/eTZs2bZg+fTrBwcFmK5/mzJmTsWPHMmrUKE6dOkX27NmpUKECY8aMoU2bNsCX660Ambwdq2Ty9heSJw9cvgxKsYcKVGKPdulPOtCB2cbEIlcuuHjRgoEKIYQQCc/naO+Eh4djZWXFo0ePuHr1KtOmTePp06fUq1ePpEmTcvLkSSZOnIi9vX2svF58ceLECTp16kSHDh3YvHkzd+/eZfv27aROnRqIqDeTlStXsnfvXsqVK0fTpk2/eLySWMQiSSy+gDt3INKkrTCsSMsjnmFc07oKO9hBNfPyGTN+4SCFEEKIhOtztncaN26Mq6srjo6OzJ49m0qVKpElSxaqVKlClSpVEt0Ss+3btydHjhz07t0bgP79++Pk5ES/fv0ICAhgy5YtNG7cGKWUWYJhKTIUSsQva9ca5038x5pw6rNWO95DRZ6Twnig1xvLCyGEECLOMt3jfvDgAdWqVWP69OmMGTOGzZs3kz17dsaOHatNQk5MSYWpN6JkyZLacevWrTlw4AAA8+bNY/Pmzej1erM5FJbsM5DEQsQvK1bAGx+Yhvp12vfhWLOeusYDpYzlhRBCCBFnmZKF69evs2bNGv79918Avv76a27dusWZM2csGJ3lWFlZMWbMGHLkyKEd58iRA3d3dyZMmMCmTZvo1asXYJ5MWDL5ksRCxB8PH8LRo28lFuVrJSVF8oiJXKtpaPxGKWP5R4++ZJRCCCGE+Ajly5enVq1a7Nmzh5CQEOzt7XF2dsbLywuw7J14S3F2diblG8vkDxw4kN9//52AgAC++eYbgC+26tOHWFs6ACGibd0646RspYx7U1hZwZQp2HToQL32OubONRbbRWV8rFLhEv4s4nGdO1subiGEEEK8xWAwYDAYsLa2ZtOmTZw8eZJdu3Zx8eJF1q5di1KKzJkz07JlSyBxDYN6n0yZMtGvXz/Sp08PQFhYGNbWcaNJHzfSGyGiY+VKY1Kh10P27MZlZH/8EXQ6GjaMKBaGDRtcfzCWU8r4OCGEEEJYjKm3Yd++ffTs2ZOwsDD0er3WIH706BFubm5MmTIFT09PSpYsyejRoxkwYABJkyZNlL0VJpHrrkePHoSGhtK1a1fq1asHECcmbZtIYiHih6dP4eBB4/c//mhMKvLk0S5XrAjJk0cUX11wOHToYDw4cACePftysQohhBDCjKm3ITAwkOXLl1OhQgWWLFnCy5cvAejQoQMdO3bkm2++YdKkSWzatAl3d3ftrnxi7q2IXHcrV66kYsWKLF26lICAALPrcYEkFiJ+ePkScuc2rvI0Ywa8sY61rS3UqRNxvGO3FX6j/zCWz50b/vvwCSGEEMJyqlevzrZt2yhQoABz585l6tSpPIt08y80NJSUKVPSuHFjFi1aZMFI457q1auzdetWChQowJw5c5g2bRovXrzQrseFXh3ZxyIWyT4WlrVlC3h6Rhz/9Rd8953l4hFCCCESoo9p75j2nzAYDOj1ep49e8b8+fMZO3YsGTJkYOLEiZQrV04rf/HiRR4/fkzFihW1xyRWMa07S5LEIhZJYmFZwcGQOjX4+xuPa9eGDRssG5MQQgiR0MRme+fBgweUKVOG27dvM3DgQIYMGRJLUSZ8cbHuEm/6JxIcOztjMmGyfXtEkiGEEEIIy1mzZg2urq6sX7+e8HDjEvEGg4F06dLRqVMnWrRoQYsWLbTzIkJ8qjtJLESC0qhRxPfBwcbhUUIIIYT4st5s4BYsWBBPT08mT57MiBEjuH37tja8qWLFijg4OODh4QHEnT0Z4or4VHfyPycSlKpVwdEx4njVKsvFIoQQQiQ2phH2pgauKcHIli0bw4YNo0qVKixbtow2bdpw7NgxAIYOHYqbm5tZeRHBVHdVq1aN83UniYVIUJIkgW+/jTjeutW4oJQQQggRW2R6atTCwsKYMGECBQoUoEePHty7dw+9Xq/VV/r06enfvz8bN27ExcWF4sWLkz59eq5cuULXrl2BuLV0alySPn16+vXrF+frTiZvxyKZvB03rF0LDRpEHK9YAY0bWy4eIYQQ8Ytp9Z2XL19y5coV3NzcuHPnDo6OjhQoUMCsrGnFnsTkXe2dyZMnc/DgQYYOHcqff/6Jl5cXO3bs0MoEBQVhZ2en1dfFixc5c+YM+fPnJ1++fIl+9af3MRgM6HS6OF93kljEIkks4obAQHB1hVevjMeNGsnm20IIIWJu0KBBbNu2jRs3bpA/f37u3r1LYGAgDRo0oF27dhQsWNDSIVrEu9o7pUuXZujQoVSsWBGAjh07kjp1aoYOHcqRI0c4ffo0nTt3TpTJWGyJ63Vn+dRGiFhmbw+1akUcb9kCr19bLh4hhBDxU40aNdixYwfXrl1j6NChzJs3j169enHo0CG++eYbGjduzI0bNywdZpwQEBBAtWrVyJIlizbWf9CgQZw6dQqDwcCMGTMIDQ0FzIeSyf3tmImcVMTFupMei1gkPRZxx6pV5sOfVq82Hx4lhBBCRMe77hBv3bqVcePGATBu3DgKFSr0pUOzmHe1d16/fo1SiqRJkxIWFoa1tTV9+vRBp9Nx9OhRdu/ejbW1tQUjF5+b9FiIBKlmTWPPhcnq1ZaLRQghRPz1rmEnNWrUYOLEiaRLl47vv/+eEydOfOHI4h4HBweSJk0KoCUQ3333HRMmTMDDwwNra+s4sXJRQhAWFiY9Fgmd9FjELQ0bwpo1xu+TJoWnT82TDSGEECI29OjRg7CwMMaPH4+dnZ2lw/ns3tfeMfXwRO7pWb9+Pd988w0ZM2aMM5OM46Pw8HCsrKwsHcZ7yf+sSLAaNoz4/tUr407cQgghxKd6/fo1+/bt04779u1Lrly5ePbsmQWjiluGDx9O+fLluXnzJnXr1iVjxoyA5Tdwi89MSUWWLFkYMWKEtgt3XCL/uyLBqlXLuK+FiQyHEkIIERt2797NnDlzMBgMhIaG4u7uzk8//YS7u7ulQ7M4Uy9F5cqVAeNQqKlTp/L48WPA2KMhg2U+3osXL/D09GTfvn2cPXvW0uG8RRILkWA5OUH16hHHGzdCUJDl4hFCCJEw1KxZkwwZMnDq1ClsbGy083F5GdAvrUSJEvz555+UL1+eZcuWMXz4cO7fv2+2F4OIvmXLlgGQIkUKRo4cSbFixTh69KiFo3qbJBYiQYs8HCogAHbutFwsQggh4i/TpOMnT55w8uRJ7ty5Q69evZgyZQoLFy6kQ4cO+Pj4WDhKyzP1RiilyJEjB8OHD6dnz55s3bqVokWL8ueffxISEmLhKOOXa9eusXDhQs6fPw+Ak5MTw4YN48cffwTi1rKzkliIBM3TE2xtI45lOJQQQoiPYZobMHjwYGbOnEmGDBm4dOkSO3fu5PTp09SpUwcXFxcLR2k5psTLz8+PV69eab0SVlZWNGrUiOPHj2NnZ0enTp1YvHixJUONVwwGAx4eHnz33XdmPRQBAQHaz2Rc6gGSxYRFgubsDNWqwaZNxuMNGyAkxDzZEEIIIaLDtNGbr68vyZMnp1GjRqxcuZIxY8bEqbvGlmBq5E6ZMoXhw4czadIkmjZtSooUKTAYDLi4uNCrVy+8vb1p2rQpEPd3kbYkU934+flx/PhxduzYwZEjRzh16hROTk48fPiQ7t27U7hwYUuHakZ6LESCF3k4lJ8f7NpluViEEELEP6akYfHixRw8eJDkyZMTEhJC/vz5uXXrFkeOHJEG8n/atGlD3759Wb16NSNGjODKlSta0uHo6IjBYMDBwUGSimiaP38+c+fOpVy5cmTNmpWLFy+SMWNGvv/++ziXVIAkFiIRqF0bIs2tk+FQQgghYsTUAH769CkLFiwAwNbWFhsbGxo2bMiBAweAuDXW3VIyZMhAnz59aN68OVu3bqVWrVr89ddf3Lx5k4EDB1KuXDlA6upDIu8FsnDhQtq1a8fWrVtJkyYNtWvXplKlSpYOMUqSWIgEL3lyqFIl4nj9eggNtVQ0Qggh4pPIDeC2bdvi6urKnj17CAkJITAwkOLFi9OnTx8LRhj3ODs788MPP3Du3DmaNGlCx44dKVasGCVKlKD6f8s1yn4WH+bj48OJEydYtWoVYKyzatWqsToO3yGV/1WRKEQeDuXjA3v2WC4WIYQQ8Yept+K3335j4MCBrFu3jubNm1OtWjUaN27MvHnztDIytMdIKUV4eDjW1taMHDmSx48fs3//fmbNmgVETPQWEaLqwUmRIgUtWrTg9OnTXLlyBaUUdnZ2HD9+/J2PsTSZvC0ShTp1wNoawsKMx6tXGyd1CyGEEB8SGhqKg4MDtWrVomnTpgwdOpRBgwbh6upKhgwZAGNjWe7CG+l0Om2XaDDOrcidO7d2LPVkFHmeienf8PBwrKysuHDhAlOnTsXb25vDhw+zc+dObG1tKVq0KL/88stbj48rdCoupjvxlL+/P87Ozvj5+ZEsWTJLhyPeUK0a7Nhh/D5lSvD2NiYbQgghRFRMDbdz586RP39+7fzy5ctZt24dK1assGB0liPtndjj7+9PpUqVGDduHOXLl9fOX7t2jV27dtGgQQP27NnDtm3bmDBhAkFBQaRLl85yAX+ApIwi0WjUKOL758/hn38sFooQQoh4QKfTERQUxLBhwzh9+jQAISEhNG3alFevXnHixAkLRxg/mIZGibeFh4eTIUMGfvrpJ3r16sWlS5cA8PDwoHPnzqRJk4ZmzZrx8uVLdu7cGaeTCpDEQiQidetCpJ5ZWR1KCCHEe4WHh5MkSRJq1qzJzp07AeNqUAADBw4kX758lgwvTjMlEgEBAW8NjRIRXFxcmDdvHl27duXChQsMHjyYQ4cOaddN9Th+/Hj++uuvOJ+gSWIhEo1UqSBSLyNr10Ic/3wKIYSwICsrK549e8arV6/YunUro0ePZv369YwfP55Xr15hZ2dn6RDjLFMiUadOHQoXLsyDBw8sHFHcYpqJYDAYSJ48OR07dmTw4MH4+vpSqVIl+vfvz+PHj7GysiI8PJx06dLRq1cvrKys4vTkd0ksRKISeTjU06fw39LjQgghxFuOHTtGrVq1uHfvHq9evWLChAksW7aM4OBgcuTIAcTNlXniCoPBQIcOHciUKRPLly+3dDhxhsFgQKfT8fTpU7OJ7CVLlmTHjh20b9+e0aNHazuUW1lZYWNjo+1dEZcnv8vUVZGo1K0LnTuDKdlftcq8F0MIIYQwKVq0KKtXryZDhgyEh4fTvXt3Bg0aRIoUKbQ78nFtVZ64ICQkBFtbW/R6PU2bNiV16tTs+W+d97i4ktGXptfrefz4MWnTpqVr165069aNbNmyaddbtmzJ06dP6dy5MxCxUlR8EHdTHiE+gzRpoGzZiGMZDiWEECIqz58/Z8qUKVpSYWVlRc6cORk7dmy8aeRZwvPnz/nzzz+BiN6cihUrMnz4cEuGZXFv9mwlS5aMP/74g2vXrjFw4EBtDg8YE9rXr19rS/TGp583SSxEohN5szxvbzh82HKxCCGEiJtu377Nli1buHDhgtaw++GHH7TVoUTUbt++zebNm7lw4QI6ne6t+QCJrbfCYDDw4sWLt963vb09bdu2pUePHvj4+NCqVSuGDh2Kj48Pw4YNw8fHh1SpUsW7oXaSWIhEp359iPz5XrXKcrEIIYSwjDcbbAaDwawRXKhQIZo0acKBAwe0lXiePXvG1q1btfLibW/Wm16v59q1a0DinI8yduxY6tSpQ/HixVmwYAFgrAelFDY2NlSrVo2tW7cyatQo5s6dS8qUKVm1ahW//vorEP9+zmSDvFgkG8bEH2XLRkzcdneHe/cgDs+FEkIIEctMY/0HDx5M69atyZIli3Zt37592m7HZ8+epWTJkhgMBuzs7Pjzzz9xdXW1YOSW92Z7x7Tr+JEjRzhw4ADbtm0zqzcHBwemT59O6tSpLR36F3X27FlatmzJiRMn2Lp1K2PGjGHEiBFUqFABMNajg4MD1v/t1vvy5UuuXLmCg4OD2U7l8YlM3haJUsOGEYnFw4fg5QUlS1o2JiGEEF+OTqcjICCADRs2MH/+fH744Qfatm1LunTpuHv3LjqdjjFjxjBjxgxcXFxo0aIFKVOmxNXVVSYgv8G0StGGDRtwcXFh5MiRzJo1i+TJk2v1ljp16kRXb/PmzePbb7/F1taWOnXqoNfrGThwIHv27CEsLIzu3btr81EAHB0dKVy4sAUj/nSSWIhEqUED6N494nj1akkshBAisXFycmLv3r3MnTuXbdu28fDhQzp37kzLli21MjNmzKBy5cr06NGDjBkzAolvnkB03L9/HysrK/73v/8BkC9fPqpWrUrXrl3JnDkzkPjqrXz58qRKlQowDmn69ttvWb58OYcPH+b69eu8evUKW1tbrccnIUgY70KIGEqXzjyRWL06YglaIYQQiUN4eDjJkyenS5cu/Pjjjxw7dowqVaowc+ZMrYy9vT1t2rRh/vz5Fow0bjMYDDg7O+Pj48PGjRsBSJo0KR06dGDDhg0Wjs5y6tWrR8GCBYGIpKp169YsWbKEuXPn0q9fPwtG93lIYiESrcirQ927B8ePWy4WIYQQX0Z4pDXGTas92dnZ0bBhQ06ePEmjRo3o3Lkz+fPn5/HjxwCUK1eOfPnyAfFvMu3n8GYd6PV6HB0d8fDwYNmyZWzbto2rV69y9uxZQkNDo3xMYuHk5AREJBbly5fn2rVrPH/+nK+//hqI2xvexZRF38kff/xB/vz5SZYsGcmSJaNEiRLaagtgnFg1ePBg3N3dsbe3p3z58ly8eNHsOYKDg+natSupUqUiadKk1K5dm/v375uV8fHxoWXLljg7O+Ps7EzLli3x9fU1K3P37l2+/fZbkiZNSqpUqejWrRshISGf7b0Ly2vQwPx49WrLxCGEEOLLMSUTRYoUYdiwYQQEBJhd/+GHHyhdujRt2rQhTZo0AGTPnp369esDCasR+DEiD9sxJWmLFi2ibt26eHl5sX37dnr27EmbNm1IlSoVP/30E5D4hkFFJTw8HFtbW1auXKntRB4WFmbhqGKXRT8d6dOnZ/To0Zw4cYITJ05QsWJF6tSpoyUPY8eOZcKECUybNo3jx4/j5uZGlSpVzH4J9OjRg3Xr1rF8+XIOHjzIy5cv8fT0NLsj0bx5c86cOcO2bdvYtm0bZ86cMRs/GR4eTq1atXj16hUHDx5k+fLlrFmzhp9//vnLVYb44jJmhGLFIo5XrQJZI00IIRKeN++WBwcH4+npyfbt2+nevTuHDh1CKUV4eDgFChQgb968FPvvD4Qsnhnh0KFD2NnZ8ddffwERSVqdOnUYMGAAy5cvZ9OmTXz77bccPnyYHj164ODgAEhiYaq72bNnkzJlSgoVKgSgrQiVYKg4xsXFRc2ZM0cZDAbl5uamRo8erV0LCgpSzs7OaubMmUoppXx9fZWNjY1avny5VubBgwdKr9erbdu2KaWUunTpkgKUl5eXVubIkSMKUFeuXFFKKfX3338rvV6vHjx4oJVZtmyZsrOzU35+ftGO3c/PTwExeoywrHHjlDKmE8av48ctHZEQQojPxWAwKIPBoJRSKiAgQC1dulRVrFhRFSxYUC1atEgppdTly5dVhgwZVFhYmCVDjZNevXql+vXrpwoXLqwAtXv37rfqyWAwqDp16qhNmzZZKMq4yVR3RYsWVT/88IM6cuSICg8PV0qpBPWzFmf688LDw1m+fDmvXr2iRIkS3Lp1C29vb6pWraqVsbOzo1y5chz+b6vkkydPEhoaalbG3d2dvHnzamWOHDmCs7OzducBoHjx4jg7O5uVyZs3L+7u7lqZatWqERwczMmTJz/r+xaWJcOhhBAi4Tp06BBdunShUaNG7N+/H51Op905d3R0pFmzZqxZs4bq1avTsWNHUqZMSePGjalRowZWVlZmox8EODg40L9/fzp27AjAr7/++tbkbJ1Oxy+//CKT3d9gqrvu3btz9+5d+vTpw/r164GInp+EwOKJxfnz53F0dMTOzo6OHTuybt06cufOjbe3N4A2vtEkTZo02jVvb29sbW1xcXF5b5moNmRJnTq1WZk3X8fFxQVbW1utTFSCg4Px9/c3+xLxS5YsEHnJaBkOJYQQCUeHDh2oVKkSBQsWpEOHDtq4djCObTcYDCRPnpxRo0Zx/vx5Bg8eTJ8+fRg7diwgw3dMVKQ/jI6OjjRq1Agwzj1p2LAhjRs35siRI4DxRnHx4sWZOnWqdpyYvVl3zZs3Z+7cueTMmfOtuksILJ5Y5MiRgzNnzuDl5UWnTp1o3bo1ly5d0q6/+aFW0dhc5c0yUZX/mDJvGjVqlDYh3NnZmQwZMrw3LhE3RV4d6t9/4cwZi4UihBAilqxfvx53d3fq1avHgAEDWLJkCbNnz9ZWepoyZQpPnz7VymfNmpWuXbtqi72ATNQ20el07N27963zU6dOpWPHjqxevZrFixcDEXVmGgWSkO7Gf4yo6i59+vTMnj37rbpLCCz+ibG1teWrr76icOHCjBo1igIFCjB58mTc3NwA3uoxePLkida74ObmRkhICD4+Pu8tY/olEtnTp0/Nyrz5Oj4+PoSGhr7VkxFZv3798PPz077u3bsXw3cv4oLIiQXIcCghhEgIXr9+TZ06dQAICQmhUKFCZM2alRMnTnDv3j1mzpz51t94JV3WUTp+/DiVKlWiVatWeHl5mdVT+/bt+fHHH+nTpw+QeJeVfZfEVncWTyzepJQiODiYLFmy4Obmxs6dO7VrISEh7Nu3j5L/7WxWqFAhbGxszMo8evSICxcuaGVKlCiBn58fx44d08ocPXoUPz8/szIXLlzg0aNHWpkdO3ZgZ2enzdqPip2dnbZUrulLxD/ZssF/S0kDMhxKCCESgubNm9O4cWPAeBMToEGDBqxfv54BAwbQoUMHwHyojgx9MmdqBBcoUICVK1fy8uVLRowYYTavIkmSJNy5c0fbXTux91C86c26W7NmjXbtzbpLED1kFpo0rpRSql+/fmr//v3q1q1b6ty5c6p///5Kr9erHTt2KKWUGj16tHJ2dlZr165V58+fV82aNVNp06ZV/v7+2nN07NhRpU+fXu3atUudOnVKVaxYURUoUMBshn316tVV/vz51ZEjR9SRI0dUvnz5lKenp3Y9LCxM5c2bV1WqVEmdOnVK7dq1S6VPn1516dIlRu9HVoWKv0aMMF8d6uxZS0ckhBDic2jSpInS6XTq6dOnSimlrRIllPL391dHjx41a2dFdurUKdWiRQuVJEkSBaiNGzeqChUqqB49eiilEtbqRrHt9OnTWt01b95c7d69O0HWnUUTi7Zt26pMmTIpW1tb5erqqipVqqQlFUoZP+yDBg1Sbm5uys7OTpUtW1adP3/e7DkCAwNVly5dVIoUKZS9vb3y9PRUd+/eNSvz/Plz1aJFC+Xk5KScnJxUixYtlI+Pj1mZO3fuqFq1ail7e3uVIkUK1aVLFxUUFBSj9yOJRfx19ap5YjFwoKUjEkIIERtMS3qa3LhxQ02ePFkplXAac7EhLCxMNW7cWJUoUULlzZtX9evXT4WGhqrw8PC3kq+dO3cqQNnb26t69eqpZ8+eKaUkSYvKm3Vy+PBhVbZs2QRbdzqlZNBHbPH398fZ2Rk/Pz8ZFhUP5c8P588bv8+VCyKtISCEECIeCw8Pp1y5cgwaNIgqVapo5yPvIp3YzZo1i127drFy5UouX77MgAED6NWrF6VLlwbg3r17uLm5YW1tTUBAAM7Ozty7dw8nJyecnZ2jtbhOYmYwGMyWO/b390cpleDqTj5NQvznv9XzALh8Gf7bAF4IIUQ89+rVK9zc3OjWrRs9evTgwoULgHFMu9xfNdq0aRP16tUDIFeuXDRo0IDffvsNMO731atXL2xsbMwawMmSJdNW0EooDePPRa/XJ4q6k8RCiP/I6lBCCJEwJUuWjPnz59OrVy8uX77MkCFDOHDgAGBs1CX25EIpRb9+/ShSpIh2rkWLFjg4OLBr1y4WLFigLXiT2PelEO8nQ6FikQyFiv/y5IkYApU3b8TQKCGEEPFX5CFPR44cYejQoezdu5cePXrQs2fP9y4tn5gEBwdjZ2dHeHg4VlZWHDlyhGHDhvHo0SO8vLyws7MDpL3zKUw/i6Y6Tmikx0KISCIPh7pwAa5csVwsQgghPp7pzvrz58/N5lGUKFGCrVu30rVrV8aOHUvt2rUJCwuzVJhxiilxMDV4CxQoQFBQEPnz59cSDvFxTPfxQ0JCeP36dYJMKkASCyHMyHAoIYRIGKysrAgICMDd3Z3OnTtz48YNs+stWrSgYcOGDBs2DGtr6wSxOdmnMjV+g4KCuH//Pg4ODmzZsoVJkyYBCWsuwJdmqrsZM2bg6OjI+vXrLRvQZxKjxMLf3z/KD154eDj+/v6xFpQQlpInD+TIEXEsiYUQQsRftra2TJs2jVu3bvHrr7+yfft2rfFcsGBBQkNDyZUrF5BANif7RKbG75IlS8iYMSO//vorfn5+uLi4AFJHsaFatWr069ePhQsX4uPjY+lwYl20f0LWrVtH4cKFCQoKeutacHAwRYoUYdOmTbEanBBfmk5n3mtx9ixcv265eIQQQnw8Ozs72rRpQ69evXj58iVt2rRhyJAh+Pr6Mm3aNG7fvk2GDBksHWac07ZtW2bNmoWXlxc///wzu3bt0oZBybCx6DMlsb6+vtpxnjx56N27N8WLF9cStoQk2onFH3/8Qd++fXFwcHjrmoODA//73/+YNm1arAYnhCVEnmcB0mshhBDxmbW1NVWqVGHz5s2MHj2aBQsWkCJFCv744w/69esHyEpHb9LpdLRs2ZJevXrh5+dH586dmTp1KmCsTxE9ph6giRMn8vLlS20FMhcXF/73v/9ZOLrPI9qJxYULFyhfvvw7r5ctW5bzsoSOSADy54evvoo4lsRCCCHiN1Pi0KpVKy5cuMCRI0dYunQpjRs3BkiwE2k/hZ2dHTVr1mTz5s388MMP9O7dmwIFCrBy5UpLhxavvHz5kgcPHrBo0SIg4S9vHO3EwsfH573dX6GhoQlyrJhIfN4cDnXqFPz7r+XiEUII8WkiJw6Ojo4UK1aMAgUKWDCiuMXU0L1//z6jRo0iNDTU7Hrfvn1p3bo1L1684NWrV5YIMc6KKkkwGAzanGRHR0dat27NjRs3uHLlCkopLl++zIMHD750qF9EtBOLzJkzc+LEiXdeP3HiBJkyZYqVoISwNBkOJYQQIrEwDdm5fPkyv//+O02aNGHLli1m82q7du1KgwYNtF4eYWRKLEaMGKGtnqXX69Hr9dy/f5+ffvqJyZMns2zZMurXr0+xYsUYPHgw165ds2DUn0+0E4v69eszYMAAHj9+/NY1b29vfv31Vxo0aBCrwQlhKV9/DVmyRBxLYiGEECKhq1ChAsuWLcPOzo6xY8cyf/58reF84sQJXrx4QdKkSRP0UJ6YMq2UpdfrGTt2LPXq1WPXrl2EhISglCJDhgxMmDCBKVOmkCtXLjZs2MD48eOpUKGChSP/PKK983ZAQAAlSpTg7t27fPfdd+TIkQOdTsfly5dZsmQJGTJkwMvLCycnp88dc5wlO1EmLH37wrhxEce3b4N0ygkhRPymlMJgMMi8ivfw9vZm8uTJTJo0idy5c1O7dm1mzpzJH3/8Qd26dfH19cXFxUXaO2/YuXMn8+fP5/nz5zRv3pzWrVubXW/Tpg3Fixfnxx9/tFCEn1+0EwsAPz8/+vXrx4oVK7T5FC4uLjRp0oSRI0eSPHnyzxVnvCCJRcJy7BgUKxZxPH489OpluXiEEEJ8nPDwcKysrAgMDMTe3t7S4cRZSimzTfDu37/PsGHDOHnyJPXr16d///6AtHfeFLnebt68yeTJk5k9ezaVK1dm7Nix2l4pL168oGnTpmzYsCHB/hzGKLEwUUrx7NkzlFK4urrKToz/kQ9awqIUZM4Md+8aj4sXhyNHLBqSECKRMBgM6HQ6+fsaSwwGA3q9nubNm3Pu3Dk2bdpElsjjXYUZg8GAUsqsV8dUh0opAgICEn17x1Qf70pW9+7dS6VKlbC2tmbjxo1Ur16d8PBwtm3bRq1atd5K4hKKj9pCUafT4erqSurUqRNkpQgBb68O5eUF9+5ZLh4hRMLl7+8PGEcGvHz5Er1eL39fY4mpARgeHk79+vXJkSMHS5YssXRYcZper9eSCtP9Z9NcAvm5NDLVR5cuXShUqBBeXl5m14sXL84PP/xA3759qVixImBcnaxWrVpAwq3HaO9y8vXXX0erEk6dOvVJAQkRlzRsCBMmRByvXQvdu1suHiFEwmG6Y3nu3DlGjhyJl5cXBQsW5Pnz5zg4ONCsWTNq1KhBmjRpgIgGsogeU32tXr2aEiVKkCFDBho2bIibmxsrVqywdHginnqzp6F169b4+voyaNAg6tevT+PGjXFxccHe3p68efMCYGtrm2B7KN4U7aFQQ4YMidYTDho06JMCis9kKFTCYzAYJ2zfv288LlUKDh60bExCiIRl+/bteHt7kzZtWi5fvoy1tTWHDx9m69atGAwGOnXqxKBBg0iSJImlQ413DAYDHTt2JFu2bFHudJxYGnuf6s16SoztndDQUGxsbLTjyIn++fPnmTVrFps2baJkyZKMGjWKVKlSkTNnTtasWUPRokUTzc/aR82xEFFLjB+0xKBHD5g8OeL4/n1Il85i4QghEjhTA8THx4fly5czadIkkiVLxtixY6lQoUKiaaDElps3bzJ//nx69epFihQpCA4ORqfTYWtra+nQ4qywsDCsra05d+4cwcHBFClSxOx6YmvvXLx4kYULF7Jx40ZatWpF7969o+yFOHjwIL179+bYsWNkzZqVzJkzs2vXLgtG/uXFSmLh4+PD4sWLmTt3LmfOnImFsOKnxPZBSywOHoQyZSKOp06FLl0sF48QImExNU7eNdTp9OnTzJw5k5s3bzJ9+nRy5MhhgSjjp9WrV+Pl5cXGjRtxdXUla9asWFlZkTdvXnr37m3p8OK8rl27Mn36dDZu3Iinp6d2PrG1d2rXrk316tX55ptvGD9+PE5OTsybNw8wrjgWHh6uJarh4eGcOHGCR48eUbRoUdzd3bVVyRKDT0osdu3axdy5c1m/fj2pUqWifv36TI58azeRSWwftMTCYID06eHRI+NxuXLwzz8WDUkIkcg8ffqUn3/+mfv377N06VLc3NwsHVKcZUrUzp8/T+/evfnxxx+5efMms2fP1pbGL1asGE5OTtL7E0nknjI7OzuSJEnCnTt3WL16NTt37mT9+vU4ODgAiau98/jxY8qXL8/58+extjZOTS5dujR9+vShTp06TJkyhTx58lCpUiX5eSIGk7dN7t69y/z585k/fz4vX77Ex8eHlStXyq7bIsHS66FBA5g2zXi8fz94e4P8XRdCfC6mBsrVq1fx8PDA1dWV2bNnM27cOP79919JLN7D1LDz8vJi0KBBlCxZEoDXr1/z4MEDGkZa7i+xNwIjM9XF7NmzKVSoEJUqVSJLliz07t2bAgUK4ODgkCgbzv7+/nz//feEhIQAYG1tzS+//MKmTZuoU6cOs2bNYu3atRaOMu6I9vISK1eupGrVquTKlYsLFy4wefJkHj58iF6v1zb+ECKhirzsrFKwbp3lYhFCJHw6nQ5fX1/GjRtHaGgoBoMBOzs7+vTpI39zo8lgMLBw4ULtuHv37pw8eVLb4FdELVWqVKxatUo71ul0VK1aVfs+scmePTsdO3bE1tYWa2trwsLCKF++PGFhYbRv357ChQvj4eGRKJOuqEQ7sWjevDmFCxfG29ubVatWUadOHZn4JBKN0qXhvxUfAVi92nKxCCESPlMi4eTkxMqVK7W5F3Z2dri4uFg4urhFKaXttWAwGLTzFSpUwMrKioULF3L79m0eP37MpUuXSJ48uYUijVveVW+lS5dGr9ezcOFC7ty5w9WrV9myZYulwowTnJ2dtWFQ1tbWODo6UqpUKebOnUu9evWAiP0+ErtoD4Vq27YtM2bMYN++fbRs2ZImTZrILzeRaFhZQf368McfxuN//oEnTyB1aouGJYRIoPR6vbYO/tatW3F1dcXDw4ODBw9SpEgRcubMaekQLco0GTbyEqBKKS0Ba9y4MU5OThw7dowFCxaQP39+MmXKRL9+/dDpdIlqMm1kH1NvadOmpW7dulrZxHpX/s26a9u2LXZ2dtSpUwdA9pj5T4wmbwcGBrJy5UrmzZvH0aNHqVatGlu2bOHMmTPaJiCJWWKazJQY7dkDlSpFHP/5J3ToYLl4hBAJ199//8306dNxcnJi7969ODo6kjJlSkqWLMkvv/wicyz+U6FCBXLmzMmoUaPMeiJWr15NqVKlUEoxcOBAfvnlF1KmTImzs3OiTCjeFJN6c3V1JVmyZOj1ei2xSMztnQoVKuDh4cGYMWOk9ysKH70q1PXr15k3bx6LFi3i5cuX1KpVi4YNG1K/fv3YjjHeSMwftMQgLAzSpoVnz4zHVarAjh2WjUkIkTAFBgZy7NgxypUrx7Fjx1i8eDFTpkzh9evX2so8Av766y9mzJhBsmTJaNeuHbVq1SJp0qRAxAZmK1asYOnSpWzYsMHC0cYdn1pvibm9Y6o7JycnfvjhB2rWrImjo+M7l4tObD66BrJnz86oUaO4d+8eixcv5vXr1zRr1iw2YxMiTrG2Ng6HMtmzJyLJEEKIT2W6z/fo0SMOHTpEuXLlAChatChPnz5l8eLFklS8oWXLlvz++++kT5+eiRMnMnXqVAIDA4GIoSlNmjQhMDCQs2fPWjLUOEXq7eOZ6i5DhgxMnDiRadOmERQUJEmFiYpFjx8/js2ni3f8/PwUoPz8/CwdivhMduxQyrgulPFrzhxLRySESGhWr16t6tSpo/z9/VV4eLhSSqnr16+rypUrWziyuMvf319NmDBBOTs7q2+++UatWbNGKaVUWFiYUkqpa9euWTK8OOtj603aO++uO4PBYOHILCtWdt4WRom5azCxCA017l/x4oXxuHp12LrVsjEJIeK38PBw9Hq92aTYUaNGUbJkSa3XIjg4GF9fX9KkSZNoJx4DhIWFYW1tzdWrV6PcgfzIkSOUKlWK1KlTs3fvXrOleVUinngc2/WWmNo7n1J3iZH02wgRAzY28N/iGADs2gWyJLoQ4mNs2rQJACsrK63hdubMGZYuXcrevXv58ccf6dGjB3369KFVq1Y4OTlp5RMr05Kf7dq1w8PDgzVr1hAUFAQYG8Bff/01PXr0YMSIEW818BJrUgFSb5/iU+ouMZLEQogYatQo4vuwMJD5gEKImLp16xZ16tShbNmy/P3339o+AvPmzePChQs0atQInU5HaGgoOXLkoH///trOxwKmTJlCxYoVmT59OlOnTuXp06fodDqSJElCQEAAqVKlAmRvgTdJvX08qbvokaFQsSgxdQ0mZiEhxs3yfH2Nx7VqwebNFg1JCBEPnTt3jmnTpnHr1i1q1apFmzZtcHZ2JiQkBFtbW65du8aQIUNYtGhRou2l8Pb2xt7eHmdn57eu3b9/n3nz5vHHH3+QOXNmBg8eDEDDhg15/vx5ot7E93PXW2Jt7zx48IC5c+fKz9x7fFSPha+vL3PmzKFfv368+G+w+alTp3jw4EGsBidEXGRrC//thwMYl5z187NcPEKI+Cc8PJz8+fPTr18/ihUrxvjx42ncuDGXL1/WGiceHh44Ozuzfft2C0drGQsWLKBbt27kyJGDtm3bcvfuXSBil+j06dPz22+/cfr0abJly0aNGjXo1q0bvXv3xtbWlvDwcEuGbzFSb59PunTppO4+IMY9FufOnaNy5co4Oztz+/Ztrl69StasWRk4cCB37txh0aJFnyvWOC+xZvCJ0ZYt4OkZcbxoEbRsabl4hBBx3/smD1+8eJHvvvuOc+fO0bp1a2bMmEGSJElYu3YtgYGBtGjRItFNPs6XLx+rV6/G3d2dIUOGEB4ezsSJEwHj31tra2uSJEmiLfPp6+vLs2fPyJgxI7a2tomuvky+RL0l1vaOKTmTn7l3i3GPRa9evfj++++5fv06SZIk0c7XqFGD/fv3x2pwQsRVlStD5N+lq1dbLhYhRPxganA0adKEQ4cOaeeVUuTJk4f+/fuTJUsWsmTJov19rV+/Pi1atDB7fGKwe/dunJ2dyZEjB05OTgwZMoS9e/eyc+dOADp27MjFixfN9g5Injw5X331ldbjk5jqy0Tq7fPS6/VSdx8Q48Ti+PHj/Pjjj2+dT5cuHd7e3rESlBBxnZ0d1K4dcbx9O/j7Wy4eIUTcZLrDafLs2TOePXtGr169GDFiBHfv3kWn06GUomHDhlSuXJmGDRsCxuFSiXUaZNasWenTpw+hoaEopUiaNCm//vor69at49KlS1y4cIEiRYpYOsw4R+pNWFqME4skSZLgH0UL6urVq7i6usZKUELEB//97QcgOFgmcIuEK3LjOLE2dD+W6e6maex1qlSpmDt3LjVr1mTDhg0MGDCA69evo9PpuHTpEkeOHNGWrIy8DG1ikyVLFqpUqYKNjY22Ola9evXQ6/W0a9eO1q1bA8iY9jdIvX1Z4eHhb908SOxinFjUqVOHoUOHEhoaChi7fe7evcsvv/xCgwYNYj1AIeKqatXA0THiWIZDiYQqctd/Ym3oxtTly5cZNmwY06dPx8/Pz2xVp8yZMzNo0CBGjhzJtWvXKFCgALVq1aJVq1aULVsWkIYfgIODg/a9jY0NVlZW5MmTh6NHj9KsWTPA/GdTGEm9fV6mz6avry9WVlZSl2+I8eRtf39/atasycWLFwkICMDd3R1vb29KlCjB33//TdKkST9XrHFeYp3MlJg1bw7Llhm/T5IEnj41TzaEiM8CAgJYuXIl3t7epEyZkidPnlCiRAnKlSunjSk2GAzodDpJON5QsGBBGjduzKFDh7h16xbLli2jQIECgLHOIjdGVq5cyYoVKyhUqBBdunQhWbJkMgn0P6Zdjx88eEC6dOkAOHnyJIUKFUrUO5B/yOeut8Te3lFKYW9vT8uWLZk6darZnOPE7qP3sdizZw+nTp3CYDDwzTffULly5diOLd5J7B+0xGjtWojcUbdiBTRubLl4hIhNc+bMYefOnfj4+GBtbY1er2ffvn3Y2NjQqFEj+vfvT6ZMmSwdZpyzYMECNm3axJo1awCYPXs2ly5d4vfff8fKyor169dTtWpVs5V5xPs5OTlRpkwZxowZQ758+SwdTrzxueotsbd3AgICmDdvHtu2baNHjx5Uq1bN0iHFGTFOLG7fvk3mzJk/UzjxW2L/oCVGgYHg6gqvXhmPGzaEVassG5MQsSk0NBQbGxuCgoLw9fXl9evX7N69m+nTp3P16lW6detGnz59SJUq1Vt34hOrX3/9lfTp09OxY0dCQkLw9fWlRYsWLFy4kBcvXvDdd99x5swZs8dID8X77du3j8mTJ+Pv70+TJk1o2LAhLi4uhIeHo9frpe7e4XPVW2Jt70T+nIaFhbFgwQKSJUtG48aN5TP8nxgnFnq9npIlS9KyZUsaNWpEihQpPlds8U5i/aAldk2awMqVxu8dHODJE0jEIwJFAhRVwvDq1SsWL17MypUrqVixIgMGDLBQdHHPy5cv8fX1JX369No509zEhw8fUrBgQbp27SpDeaLJ1GA7d+4cc+fOZefOnZQtW5bhw4eTKlUqS4cXZ33OektM7R3T778DBw7w7Nkz6tWrJ5/d94hxYnHq1CmWLVvG8uXLefr0KdWqVeO7776jdu3a2NnZfa4444XE9EETEVatMh/+tHq1+fAoIRKq4OBg1q9fT7du3ahXrx4TJkwwmziamL1599JgMFChQgXOnDmDt7c39vb2Fowuftu0aRPdu3cnMDCQzp0707t3b5IkSSJ3iz8gNustMbZ3Fi5cyObNm/nrr79IkiSJJBfvEOM+62+++YZx48Zx9+5dtm7dSurUqfnxxx9JnTo1bdu2/RwxChGn1awJkdsIMhRKJFSm+1A3b95EKYWdnR1NmjRh//79hIWFce/ePQtHGHe82VjT6/W0b9+erl27Ym9vL6s+fYCpfv766y8ePHhgdu3bb79l7ty5+Pn5sXfvXlk8IBKpt9gTHh5u9jlt3bo1BQoUYPfu3RgMBqysrLh8+bIFI4ybPnrydmSnTp2iXbt2nDt3LlH/skyMGbwwatgQ/punSdKkxtWh5IakSIiCgoLo1q0bEyZMwNHRkZCQEGxsbDh06BCpUqUiZ86clg4xTnn58iXDhw9n0KBBZr0UMh/lw0JDQylevDhBQUH89NNPNGrUCFdXV23fgP/9739UrFiRGjVqSH1G8iXqLaG2d0yf18GDB5ut9DRr1izu3LnD33//zYMHDyhatCjW1ta4ubkxc+ZMC0Yc93z0p/DevXuMHTuWggULUqRIEZImTcq0adNiMzYh4o3Im+W9emXciVuIhMa0tGyyZMlY9t86y7a2tuh0OkqXLi1JRRRu3brF3LlzKV26NPPmzcPX1xeQfQSiw8bGhqVLl+Lp6cmyZcsYO3Ys9+/fR6/Xo9fr8fLyIkuWLIDsrxKZ1NvHM31eS5UqxZw5c7TP64MHDyhUqBBTpkyhYMGCNG3alEGDBjFixAgA2SQvkhj3WMyaNYslS5Zw6NAhcuTIQYsWLWjevLmsFEXCzeDFhwUEGFeHCg42HrdoAYsXWzYmIT6XZcuW8ffff/Pdd9+RK1cujh49Ss6cOWUZ0He4ceMGs2bN4sSJE5QpU4YOHTpoewvInfYPe/36NcuWLWP48OGEhYXRsWNHfH192bhxI1evXrV0eHHW56y3hNzeMX1eT548SenSpenatavZZPejR48ycOBAtm3bJp/dKMQ4sciQIQNNmzalRYsWFCxY8DOFFT8l5A+a+LB69WD9euP3Tk7G1aFkzxyREJg22zpw4ABjx47F0dGRffv2YWtrS6pUqShatCh9+vTR7oKKCKYJnj4+PixevJixY8eSLFkyJk6cSNWqVS0dXpz35iT4YcOGMXv2bPLkyUOfPn2oWLGi9vMpInzuekuo7Z2oPq+Ojo5MnDiR6tWra+Xat29P3bp1qVWrlgWjjZtinFjIOr3vllA/aCJ6liyB776LON64Eb791nLxCPEpXr9+zdChQ83mBvj7+3P48GGqVq3K6dOnmTt3LlOnTiUwMBBH2XJe876/k35+ftSqVYvDhw/TpEkTbUiZeDellDZZ1iQwMFBW1vqAz1lvCam98zGf1507d5I5c2ayZ88u7eI3RKsP59y5c9r4sfPnz3Pu3Ll3fgmRWHl6gq1txPHq1ZaLRYhPFdVY42TJklG9enX0ej2FChUiMDCQOXPmSFLxBp1Ox8OHD0mfPj3Lli3TFjUxGAw4Ozvz008/UalSJVq2bAmQqBc9iQ6dTvfWsp6SVHyY1Fv0xPTzajAYqFKlCtmzZ9ceLyJEK7EoWLAgz549077/+uuvKViwoPZlOv76669j9OKjRo2iSJEiODk5kTp1aurWrfvWuD+lFIMHD8bd3R17e3vKly/PxYsXzcoEBwdrY+CSJk1K7dq1uX//vlkZHx8fWrZsibOzM87OzrRs2VKblGNy9+5dvv32W5ImTUqqVKno1q0bISEhMXpPIvFydoZq1SKON2yImHMhRHyTJ08evLy8qFKlCsuWLWP8+PHcvXvXrMzw4cNZvny5hSKM26ytralWrRp//PEHAwYM4MqVK9p47Lp165I6dWpKlCgBIGvhC2FhVlZW0f68SiLxftFKLG7duoWrq6v2/b///sutW7e0L9Pxv//+G6MX37dvHz/99BNeXl7s3LmTsLAwqlatyqtXr7QyY8eOZcKECUybNo3jx4/j5uZGlSpVCAgI0Mr06NGDdevWsXz5cg4ePMjLly/x9PQ0uwvUvHlzzpw5w7Zt29i2bRtnzpzRsk8w3jGqVasWr1694uDBgyxfvpw1a9bw888/x+g9icQt8upQfn6we7flYhHiU4SHh5MtWzZ++eUX6tWrx8KFC6lRowY7duzQyqRLl47Zs2dr5UWE1KlTM3LkSBo0aMCOHTvo0KED+/fvB2D27Nk8ffoUFxcXYmHF90Qjcl3JKjzRJ/X2YWnSpGHIkCHUrVv3g59XSSw+QMXQvn37VGho6FvnQ0ND1b59+2L6dGaePHmiAO15DAaDcnNzU6NHj9bKBAUFKWdnZzVz5kyllFK+vr7KxsZGLV++XCvz4MEDpdfr1bZt25RSSl26dEkBysvLSytz5MgRBagrV64opZT6+++/lV6vVw8ePNDKLFu2TNnZ2Sk/P79oxe/n56eAaJcXCc+LF0rZ2CgFxq82bSwdkRDRZzAY3nnNz89PlS1bVul0OtWoUaMvGFX8d/PmTdW0aVNlZWWlMmTIoDJnzqzWrVunlFIqLCzMssHFE6afzUOHDlk4kvjlc9VbQmjvmNqyly9fVoGBgdr569evq8aNG8vn9SPFeJ2sChUq8OLFi7fO+/n5UaFChU9Kcvz8/ABIkSIFYOwd8fb2Nls5w87OjnLlynH48GEATp48SWhoqFkZd3d38ubNq5U5cuQIzs7OFCtWTCtTvHhxnJ2dzcrkzZsXd3d3rUy1atUIDg7m5MmTn/S+ROLh4gKVK0ccr18PoaEWC0eIGNHpdHh7e5M+fXqWLl1qNtY4WbJkdOrUiYoVK9KqVStAeik+RCmFUoqsWbOybNkyTp48Sa9evfjzzz+pW7cuIMOgokun03Hr1i1Kly5N27ZtZZhyNEm9vZtpRay+ffuSIkUKZsyYgZ+fH1999RUrVqyQz+tHinFiod7RDfT8+XOSJk360YEopejVqxelS5cmb968AHh7ewPGLqrI0qRJo13z9vbG1tYWFxeX95ZJnTr1W6+ZOnVqszJvvo6Liwu2trZamTcFBwfj7+9v9iVEo0YR3/v4wJ49lotFiJjS6/VUr16dmTNn0r9/fy5fvmw21jhNmjSULFkSkD+yH6LT6dDpdNrwkwIFCtCjRw9ZZvYjpUqVihUrVnDnzh12yzjTaJN6e79Fixbxv//9jxUrVjB48GBOnz4NyOf1Y0V7AeP69esDxl+U33//PXZ2dtq18PBwzp07p/2x+RhdunTh3LlzHDx48K1rbyYy70pu3lcmqvIfUyayUaNGMWTIkPfGIRKfOnXA2hrCwozHq1ebT+oWIi5LnTo1Q4YMYeXKlSxevJgjR44wfPhwypYty5w5c3j8+DEpUqSQscYxEHkTLam36DHVk2nOpYODA05OTjRq1IgUKVKQKVMms3LCSOotZgwGA8mTJ6d79+5kyJCBmTNn8s8//zBo0CCtl0LETLR7LEyrKSmlcHJy0o6dnZ1xc3OjQ4cOLP7IrYa7du3Kxo0b2bt3L+nTp9fOu7m5AbzVY/DkyROtd8HNzY2QkBB8fHzeW+bx48dvve7Tp0/Nyrz5Oj4+PoSGhr7Vk2HSr18//Pz8tK979+7F5G2LBCpFCqhYMeJ43bqIJEOIuMg0pOnq1asEBgaSLl06evbsycqVK0mXLh0VK1YkQ4YMTJgwgY4dOwIyCfRjSWMuZhYvXszu3bvR6XTaMJ5KlSqRLVs2QOrzXaTeoseU9CdPnpy2bduyc+dOChUqRP369alQoQJnz561cITxUEwnZQwePFi9fPkyViZ4GAwG9dNPPyl3d3d17dq1KK+7ubmpMWPGaOeCg4OjnLy9YsUKrczDhw+jnLx99OhRrYyXl1eUk7cfPnyolVm+fLlM3hYfZfbsiAncoNTOnZaOSIgPq1u3rkqSJImaPn268vX11c6fPn1aTZgwQW3dutWC0cV/YWFh750gn9hFVTdLlixRrVq1skA08Ycl6i0+t3dME7CXLl2q5s6d+9Z1X19f9c0336iSJUuqy5cvf+nw4r0Y77wdmzp37szSpUvZsGEDOXLk0M47Oztrm7iMGTOGUaNGMX/+fLJnz87IkSP5559/uHr1Kk5OTgB06tSJzZs3s2DBAlKkSEHv3r15/vw5J0+e1MYA16hRg4cPH/Lnn38C0KFDBzJlysSmTZsA4x27ggULkiZNGsaNG8eLFy/4/vvvqVu3LlOnTo3W+0lIO1GKT/PsGbi5gWlu648/wsyZlo1JiA/x8/Nj8uTJ7N69m6+//ppWrVrxzTffWDqseM1gMJgNhRLvp5QiLCwMGxsbAB49esTYsWPJly8fVatWRa/Xc+rUKTw9PS0cadzypestIbR3hgwZwqxZs/j2229p27YthQsXRimFlZUV8+bNIyAggO7du1s6zPjnY7KRVatWqUaNGqlixYqpr7/+2uwrJoAov+bPn6+VMRgMatCgQcrNzU3Z2dmpsmXLqvPnz5s9T2BgoOrSpYtKkSKFsre3V56enuru3btmZZ4/f65atGihnJyclJOTk2rRooXy8fExK3Pnzh1Vq1YtZW9vr1KkSKG6dOmigoKCov1+4nMGL2JfpUoRPRaurkrJKnUiLgsPD1dKKeXj46Pmzp2rihQpogoUKKDWrl1r4cjiN9Pd0QYNGqiePXuaLWspjPbt26fKli2rrl69ana+devW6ocfflAFChRQSZIkUSVKlFA1a9ZUkydPtlCkcYsl6y0htHf8/f3VwoULVdWqVVX16tXVpk2btGv16tVTf/31l1Iq4nejiJ4YJxaTJ09Wjo6O6qefflK2trbqxx9/VJUrV1bOzs6qf//+nyPGeCMhfNBE7PnjD/PhUHv2WDoiIaLPz89PtW/fXul0OlWuXDl16tQpS4cUb4WEhKi+ffuqUqVKqS1btlg6nDjn8uXLqlSpUurrr79Wo0ePVnfu3FFKGYeqXLt2Td2/f1+1a9dOnTx5Ur148UKFhIQopd6/70piYMl6S0jtnQsXLqhmzZopvV6vKlWqpHr16qWSJUsW5Z5t4sNiPBQqZ86cDBo0iGbNmuHk5MTZs2fJmjUrv/32Gy9evGDatGmx26USjySErkERex4/Bnd3MM1x7dwZpk+3bExCRGYaprNy5UpevnxJ27Ztza4HBARQoUIFbGxsmDdvHrly5bJQpPGL+m/Fna1bt1KuXDkcHBx49eoV8+bN49y5c9pu5SLC3bt3WbRoETt37iRXrlz07duXrFmzate3bdvGuHHjZLnUN1iq3hJCe0e9sTKWl5cXv//+O8HBwXz//fc0aNCA8PBwWVY7pmKaidjb26vbt28rpZRydXVVZ86cUUopde3aNZUiRYrYS3nioYSUwYvYUb58RI+Fm5sMhxJx09ChQ1XatGlVhw4d1NGjR1V4eLg2hGfevHlq4sSJlg0wHnr16pWqX7++toiISXBwsFJKhlcoFXHn3FQX4eHhat26dSpv3rzKxcVFTZs2zWw4cr169dQ///xjkVjjkrhQbwmpvRMeHi6fx1gU4xllbm5uPH/+HIBMmTLh5eUFGHfJVpabBy5EnNSwYcT33t7w30bvQsQpPXv2ZPTo0dy5c4dBgwbx999/a3fp/v77b20DUlliNnrCw8NxcHCgffv2nDx5UjsfFhaGra0tgEzoxrjk6evXr7W60Ov11K1bl3PnzpEvXz66du1K7969tfKjR4+mTJkylgo3zpB6i116vV6rS2nHfroY/2arWLGitpJSu3bt6NmzJ1WqVKFJkybUq1cv1gMUIj6rXx8iLxe+apXlYhHiXRwdHWnVqhUTJkwgZcqU1KlTh8qVK9O7d2+2b99O06ZNAWkMv094eLjWKLGysuLMmTN4eXmxatUqOnbsyNixY/nll19YunSphSONO3bu3ImjoyO///47T58+BSKGp3Tu3Jnvv/9eW5UnLCwMDw8P2X8BqTcRt8X4r8SsWbMYMGAAAB07dmTBggXkypWLIUOG8Mcff8R6gELEZ2nTQunSEcdr1kTMuRAirjA1iHPnzs3ixYs5dOgQyZMn5/Lly8ybNw87OzttEz1hZOq92blzJ2BMJkyNt+PHj9OlSxdcXFzImjUrmzZtIiAggHz58lGjRg1A7owCfP311wwePJitW7cyaNAgTp06pdVh/vz5CQwM5KuvvgLA2toakI3dQOrtc9LpdPLZ/EQW3ccioUkIk5lE7JsyBSIvhX3oEJQsabl4EgoVaeKd7BUQO0yN5ajqUr0x0VHA06dPcXNzo2jRogwYMICaNWtqdRcYGKjtx9SlSxfat29PgQIFLBlunOTr68uGDRuYNm0az58/57fffqNBgwY0b96cnDlzMm7cOPl8R8HS9ZaQ2jumCdorVqygZMmSZMiQwdIhxWvRSizOnTsX7SfMnz//JwUUnyWkD5qIPffvQ+TfUz16wMSJFgsn3vPx8UGv12Nra8u9e/fw8PCwdEgJkiQS0XPp0iWmTJnCjRs3qF69Ot999x1ubm4ABAcHY2dnx9atW1m/fj0zZ85M9HVqasQtXLiQokWLaiuNhYeHM2zYMCZOnEhISAh58+Zl8+bNpEmTRn4WiXv1ltDaO0op8uTJg1KKLVu2mK2sJWImWomFXq+PVveQTqdL1N3lCe2DJmJPqVIRE7fTp4c7d0BuwMWcUooOHTqwdetW0qRJg7OzM3fu3KFs2bK0bNmSihUrmpVN7I2RTyH192Gmxt7t27dZsGABCxYsIGfOnIwfP548efJo5QIDA2ndujUrV660YLRxS7Vq1Th69CidO3emXbt2ZMuWDYAXL15w5coV0qZNS5YsWeTn8A1xpd4SWnvHYDBw/vx5pkyZgouLC7///rulQ4q3opVY3LlzJ9pPmClTpk8KKD5LaB80EXsmToRevSKOvbygWDHLxRNfvXjxgrVr11KjRg1OnDjB69evefz4MevWrePQoUPky5ePMWPGULVqVUuHmiBIo+5t76uTK1eu8P3333P8+HGaN2/OjBkzcHJyMisjw3qMgoKCmD9/PuvXryd9+vS0a9eOkjJG9IPiSr0llPbOjRs3tPkoYNy75+bNmxQsWFD2sPhYX2BJ20QjIa3rLGLXnTvmu3D//LOlI0p4Lly4oFq2bKkyZcqkfvvtN9k19SOY6mzFihXaLr4iau3atVP79u3Tjk17C2zdulVlypRJ9ezZU7sma+SbM+2REhwcrDZs2KBKly6t0qRJoyZNmmThyOK2uFRvCaG9s2vXLtWnTx+llLFOTZ/TxL6j+6f6qNsmf/31F6VKlcLd3V3rzZg0aRIbNmyIxZRHiIQjY0bzHorVq40phogZFanS1BsVmCdPHqZOnUr//v3ZvHmzdGV/BGtra5RSDB06lKpVq3Ljxg1LhxQn+fv7c/v2bbp3786QIUO4desWOp0Og8FA9erV+fbbb2nWrBkgPRQmkfdAMd0FtrW1pXbt2hw4cIBOnTrRs2dPsmXLFqN5nQmd1Nvnky1bNp48ecLVq1extbXVPqfSS/tpYvzb7o8//qBXr17UrFkTX19fbU5F8uTJmTRpUmzHJ0SCEXmzvDt3INK+WSKaIv/Cj+qXv7OzMx06dGDIkCGsWLGCjRs3fsnw4h1TcrZjxw5evnypnV+6dCllypRh+vTplgotTjI18hwdHfn7779p1qwZW7ZsoX///ly+fBm9Xs+TJ0/YvXu3tpCJJBVGps9ru3bt6Nu3L/7+/mbXO3ToQJkyZahVq5Y2MVlIvcU202fYYDCQPHly0qZNy/jx4zlw4AA3btxg4sSJBAQEWDjK+M06pg+YOnUqs2fPpm7duowePVo7X7hwYbOdHoUQ5ho0gD59Io5Xr4bChS0XT0ISFhaGtbW19q+npyfBwcE8efJE5gm8h06nIzg4mDlz5hAaGkqtWrXQ6XTkz5+f33//XeuxkLHGRnq9Hn9/f7Jly8aYMWPo27cvFStWpFevXnzzzTdUqVKFgIAAcufOre39IfVmZFoAJk+ePCxfvpwHDx7Qpk0bypcvD0DatGmpXr06GTNmxMbGRj63/5F6i12muunXrx/e3t68ePGCbdu2cf78eezt7alevfpb86JEzMT4VsqtW7f4+uuv3zpvZ2fHq1evYiUoIRKiLFnME4lVq2Q4VGyZPXs2ELEZFECDBg1o3Lix/JF9j/DwcOzs7Gjfvj2nT5/WzhsMBpydnSlUqBCANI4jsbe3p1OnTkybNo2WLVvi4ODA/v37WbhwIaGhoRQvXpwpU6YAMqTiTaadofv374+vry//+9//mDt3rva5nTdvXqJesv5dpN5ij+lmSurUqRk5ciSbNm1i3rx5fPfdd+zZs4cePXoAsoHlp4hxYpElSxbOnDnz1vmtW7eSO3fu2IhJiAQr8nCof/+FKD5KIoYePnzI5s2b2bp1K2D8g2D6oxCfVyv5EqysrDh37hzHjh1j7dq1dOzYkd9//50BAwawYMECS4cXJ9nY2DBgwAAGDhyIv78/7dq1Y/HixTRu3JitW7cyatQo3N3dARkGFZUkSZJQt25dli5diqenJz///DNubm6ULVsWFxcX8uXLJ3fdoyD1FnsOHz7MsWPHsLW1RSlFtWrVWLJkCdevX8fW1haQmwKfIsa/9fr06cNPP/3EihUrUEpx7NgxRowYQf/+/ekTeZyHEOItDRqYH69ebZk44jOllDa3SymFu7s7Xbt25fz584DxD8Lz588tGWK8ceHCBX788UccHBzImjUr69ev58WLF3z11Vd4enoCcufOJHI92NnZUa9ePWbMmEHhwoX5/vvvqVu3Lvv27XurrDBnSvydnZ0ZMmQIFy9epGfPntSoUYMlS5YA5hOWhZHUW+ypUKECRYoU4dixY+h0OlKnTo2np6csQBRLorWPxZtmz57N8OHDuXfvHgDp0qVj8ODBtGvXLtYDjE8SyrrO4vP6+uuInors2eHqVZCbI+9nmjvx5MkTUqdOrZ1//PgxGzduZNeuXRw8eJBChQphb2+PUorOnTtr45DFuwUGBmJvbw9A9+7dad26Nd98842Fo4pbTCs7HT9+nCJFirx1/ffff6dv376UK1eOvXv3WiDC+MfU9JA7wzETV+otvrV3TJ/hW7ducenSJebMmcOVK1fw9PTEzs6Omzdv8uuvv5InTx5Zye0TfVTNtW/fnjt37vDkyRO8vb25d+8e7dq148GDB7EdnxAJTqNGEd9fvw7/3WgX72EaS1ysWDEqV67M2bNnAbh58yZHjx6lffv21K9fn9evX9O2bVv69u0rScUHPH36lAEDBmBvb09YWBgAtWrVYtq0aXLH/Q16vZ7nz59TqlQp6taty4EDB8zqqEGDBrRp04YRI0YAaD1q4t10Op3WOJaft+iTevs4pkRh8eLFbNy4kWrVqmFjY8Pt27dJmTIlXbp0IU+ePGZlxcf5pNpLlSoVqVOnxtvbm65du5rtXiiEiFrkeRYgw6FiYsmSJTg7O9OjRw9mzJhByZIlmTNnDpUrV2bq1KkkS5YMJycnCstyWx/k7+/PxYsX2b9/v5a4Va5cmZCQEIvfDY2LkidPzrp167Czs2PUqFEsWbKE0NBQwDj38M6dO2TPnh2Qye4xJUN4Po7UW8yEhoYycOBAJk+eTMeOHdm5cyc2NjZ06tSJUqVKWTq8BCPaiYWvry8tWrTA1dUVd3d3pkyZgsFg4LfffiNr1qx4eXkxb968zxmrEAmChwfkyxdxLIlF9CilKFmyJEOHDuWbb75hxowZtG3bFm9vb61Mr169WLBggfzBjYZs2bJRv359du3apfU2nz59mrlz5wJyN/RNVlZW1KpVi4EDB5IuXTo6d+5MixYtOHr0KH379iUwMBBXV1ept2iKXE+SiEWf1FvMmeps6dKlLF++nCRJkhAYGEiaNGlIly4da9eutXCECUu0E4v+/fuzf/9+WrduTYoUKejZsyeenp4cPHiQrVu3cvz4cW2nUSHE+0UeDnX5Mly8aLlY4gvTXfQ8efIwfvx4xo0bx6FDh8ifPz9DhgwhICAADw8P/t/efYdFcX0NHP8uoBRB7AVBUGMXuzH23mKJMSrGWKImVuwxUaOvxoaa2GKsJKImMfqzxpIYW9TYBSXWYG9B7KJYQNj7/rHZhRVEWNHZhfN5Hh53CsvhsOPcM3PvncKFC2NnZycNvCQYC65t27YxZMgQli5dyowZM2jbti3vvPMO06ZNIzIyEtC+D7e1KlOmDIGBgezcuZPbt29TrVo19u7dy6RJkwC5ipxSxuPzk08+YceOHRpHYzskb6ln/L/MxcWF9evXExcXZxpX1q5dO9PFKTlnpI0UD9729vbmhx9+oGHDhly4cIG33nqLAQMGyNO2E7C1wUxCO6dPQ8LZmceOhTFjNAvHahkfMLZhwwYqVKiAp6en2fYjR45Qu3ZtihYtyubNm8mbN69pm0y9+GJTp04lT548NG7cmGXLlnHjxg0+//xzANNVd8ld0vR6vVk/99u3bxMTE2OaYlak3L1792jWrBmOjo4sXboUb29vrUOyCdaQN1tp7yQciH3v3j0CAgIoVqwYzZo1w8nJibCwMKpXr65xlOlLiu9YhIeHm55TUbhwYZycnPjkk09eW2BCpGclS5oXFtIdKmnGW/3Tp0+nVq1afPvtt9y4cQMwnDAqVqzIoEGDGDFiBHnz5jUbNCsN46Rdu3aNixcv0rp1azw8PBg4cCB79uzh4sWL5M6dG5DcJcfOzs4sP7ly5ZKiIpX2798PGMatLF++HG9vb9OUqeLFJG+pZywq+vTpwxdffMHWrVvx9/enffv2dOzYkZ07d2obYDqU4sJCr9eTKVMm07K9vT1ZsmR5LUEJkREkHMR94gT88492sVi7BQsW0LFjR1atWsWUKVO4dOmS6YRx9OhR8ufPD8hsHi8TFxeHp6cnBQoUYOPGjYDhgW9jxoyRfsbijbh69SrTp0/nzp076HQ6fHx8WLp0KX369AGkK9mLSN5eTZMmTRg4cCAbN26kc+fOTJkyhWXLlvHZZ58B0g0qLTmkdEelFB9//DGOjo4APH36lN69eycqLuTkJETKtGsH48bFL69aBaNGaRePNStWrBgjR45kxYoVBAQEsGLFCvz9/bG3t+fw4cPUqlULkCvtL6PT6bh9+zZKKRYvXsyjR4/w9PRk27ZtFCxYEEDmcE8hY55iY2Oxt7eXz14KxMXF4eXlRcWKFfn999/p1KmTaVv27NkBuTiQFMmbZYzH6M6dO2ndurVpfadOnRgxYgS7d+82Hbdy/KadFBcWXbt2NVtO+MEWQqRe6dJQvLjhAXkghUVSTp48Sd++fZk/fz4lS5ake/fudO/enYkTJ/Ldd99RqFAhZs+eDcQ/RE+82PLly/nhhx/w9fXlwoULTJs2jTJlytCwYUPTA06lgZIydnZ2PH36FCcnJ61DsRnGwjY2NtZU2Hp7e3Pw4EFatWpFhQoVtA7RKkneLGMs/OfNm8f9+/dp3bo1MTEx1KlTBy8vL3777TeaN2+udZjpjkVP3hZJs5XBTMJ6jBoF/z1TC4AzZwxP4xYGFy5coEuXLkRFRfHBBx/QsWNHihQpYtoeFRWFq6urhhHanosXL1KoUCGuX7/OxIkTmTVrVqKuruLl/vrrL+rUqcPcuXPp3bu31uHYhGXLlpkK2/Xr1+Pg4ICvry81a9akT58+UqS9gDXmzRbaO8bJP3bs2MGePXv4v//7P9O2iIgIsmfPbuqFI9KOFBZpyBYONGFd/v4bypePX540CUaM0Cwcq/Tvv/+ydOlS/vjjD4oVK8bgwYMpWbKk1mHZDOMMT6dOnSI0NJSOHTua7u4MGjSIUqVK0bNnT63DtDmXLl3ihx9+4OjRo8yZM0dmNEqhpApbpZTcbXwJa8ubtbZ3jMWE0cmTJ9m4cSMrV66kUqVKFC9enDt37lCyZEnpefOayD1vITRUtiwkfGC9zA5lYLzeodfrKVCgAMOHD2fo0KEcOnSIGjVqMHv2bJ4+fapxlLbB2Hf4zJkzLFu2jDt37pgaI/7+/mzatEnL8GyC8fP45MkTwPC59PHx4YsvvqBKlSrSKE6GMXcnT55k2bJlFCpUiNjYWPLnz4+DgwOLFi2S/CVB8pYyxpkAjbM7JSwqTpw4QY8ePXBwcKBQoUKsWbOGGzduULBgQRo2bAjIoO3XQQoLITSk05nPDnXkCFy4oF081kKn0/H48WNTf3+dTkfLli0JDQ2latWqDBw4UK6yp1Lr1q2pW7cuu3btMq1zdXVl5cqVgMwqkxJBQUGcO3cOOzs79Ho9rq6ujBkzhgIFCmgdmtUyFrZnz55NsrA1zk4mzEneUsbe3p6HDx/SrFkzKleubHoAHhgeZrljxw6GDh3KypUr6dKlC23btqVXr17ky5cPkEHbr4MUFkJoLGFhAXLXAiA4OBhXV1cmT57MzZs3gfgrS71796Zjx44MGzYMwOzZFSJeXFwcsbGxAKxbt45JkyaxevVqevfuTZs2bWjTpg2DBw82ja2QQdsvptPpiIuL48KFCwQFBQGGfMlnL+WksLWM5O3l3NzcCA0N5e2332bWrFl88803hIeHA4anbUdHRwPQvHlz5syZIzl7zWSMRRqy1j6HwropBUWKwMWLhuUqVeDQIW1j0tq9e/eYM2cOO3fu5K233qJHjx5UqVIFMDysc8CAAaxcuVKuNiVg7Ft87do1smfPbjYV+IQJE8iTJw/FixcnICCAhg0bUqNGDQoWLEiBAgVkitlkJHwK+YULF5g/fz5dunShTJky3Lt3jwcPHsj4ihcwfq7WrFnD6dOn2bhxI+fPn6dmzZoAODo6smzZMjmOn2MLebOm9o7x/74rV66wePFilixZQpEiRZg2bRq+vr5m+7Vv355Vq1bJZ+41ksIiDVnTgSZsy+efw9dfxy9fugQZva3y6NEj1q9fb3ra9siRI+nYsSN9+/bFzs6ORYsWSYM4CZ9++imLFi1i9uzZdO7cGTc3N7Pt58+f55NPPuG3337D2dlZoyhtT0BAAOfPn2fHjh3Y2dnh5eVFjhw5aNCgAX379tU6PKs2ffp0XF1dpbBNJWvOm9btnYQF//POnDlD9+7d2bdvHx06dGDu3Llky5Ytxd8vXpESaSYyMlIBKjIyUutQhI05eFApw70Lw9c332gdkfWIjY1V48ePV+7u7ipz5syqbNmy6p9//lFKKRUXF6dxdNbn4cOHasKECap27dqqX79+6vDhw0qv1yu9Xq+ePHmilFLq888/Vz/++KPGkVo/vV6vlFLqyJEj6uOPP1bBwcFq3bp1qnbt2urIkSPq1KlTKjo62mxfYWA8Nq9fv262/ty5c6pu3brq8ePHWoRl9Wwlb1q3d4zH24ABA9Sff/6ZaP3OnTuVl5eX6tmzp2mbnC/eDLlEIIQVqFIF/nvwMSDjLMBwRUmv12Nvb8+oUaM4f/4827ZtY+XKlRQvXhyllFzlfI5SCldXV/z9/enRowdHjx7lk08+YfXq1eh0OtN89y1atDB135H+xi9mvKK5ZcsW/Pz8qFSpEu+99x6NGzdm2bJllCxZksyZM5vtKwyMx+asWbNYvnw5AI8fP6ZIkSJUrlzZNEZAmJO8pYxxgo9z584xaNAgxowZw4ULF0zjoerUqUPbtm1NU8rKXbE3R7IshBV4fnaoAwfg6lXt4rEGOp3O7ESQM2dOatWqRbFixUzbhTljTtzd3enSpQubN2+mRo0atG/fnrp163LgwAEAatWqRa1atQAZtJ0c9V9P4VKlSrFhwwbT+i+//JIrV64QERGhVWhWz5i72rVrs2fPHsAwkBagS5cupuNYCltzkreUc3FxYcGCBfj5+bF161ZGjRrFiRMnsLe358mTJ2zevJnSpUsD8v/cmySTIAthJdq2henT45dXr4ZBgzQLR9gI45W4HTt2cPHiRXr06GHa5ubmxsyZMzl27Bi3b982NUbk6t3LGXMUExNDjhw5uHPnDmPHjqVOnTo8e/aMiIgI05SVwlzC3GXNmpXbt28zduxY6tevT0xMDLdu3eLDDz8EpMGXkOQt9Tw9PRkxYgTvvPMOo0ePpkqVKjRp0oSYmBg8PT3JkSOH/H/3hklhIYSVqFoVPD3h2jXD8qpVUlg8z/jEaBHPeMI8fvw448ePZ9++fXzyySdUrVoVgEyZMjFw4EDCwsKoXr262feIFzNOJ9ugQQN8fX25fv06q1atYu/evRQsWJAJEyYAiZ/0K16cu3379lGgQAHat28PyADa50neUs+Yi3r16rFnzx5WrVrFokWLKFasGP379zftI94cmRUqDWk9S4KwfYMGwaxZ8cvXroE8e0ukxOPHj1m3bh0///wzcXFx9OrVi/fffx+Afv36UahQIT777DO5epdKZ86cwcPDA1dXV7788kvatWtH6dKlTc//EC+WVO7KlCkjFwdewhbyZm3tHSnwrYecXYSwIs8/LG/tWm3isBbG6x7z58/H29ub0NBQbQOyYi4uLnTs2JFp06aRN29e2rZtS6NGjRg7diyLFy+mY8eOgIxNSY2FCxdSrFgxXF1dAahbty6jRo2SoiIFXpQ7a2ocWyPJm2USFhVyvVxbUlgIYUWqV4f8+eOXM/oEIMZG8FtvvUX58uVZvHgxz5490zgq62Q8mZYoUYIlS5Zw4MABsmbNyp9//sm3336Lh4cHer1eCosUun37Nhs3bmTdunWAIb+NGjUCkAL3JSR3lpG8ifRAukKlIWu7NShsU//+8N13htc6HYSHQ0YZI2rsLxseHm4aGGvstnPs2DHGjRvHL7/8IleMk2EcoJ2wu5Mxr9I3O2WMedqxYwcHDx5k+PDh6HQ6oqKicHBwME3bKxKT3FnG1vJm7e0d+b9OO1JYpCFrP9CEbdi1C+rWjV+eOxf69NEsHE0MHjyYIUOG4OXllWTfWTlppIzkyTJRUVGsWrWKHTt2sHPnTsqVK0fWrFl59uwZ3bp1o1mzZlqHaFXi4uIAQ3eUqKgoVq9ezfbt2yV3L2HLebPW9s6sWbPo2bMnzs7OWoeSYUlXKCGsTM2akDdv/HJG6g6l1+uJi4vD1dWVH3/8ESDJAXnSWBav06lTp9i2bRsdO3akTZs23Llzhw4dOjB48GCaNm2qdXhWY+LEiTx69Ah7e3vTcXr06FG2bNkiuUuG5O31OH78OJ999hm9e/cmJiZG63AyLLljkYastYIXtqdvX5g3z/Dazg6uX4c8ebSN6U3avn07y5Yto23btlSoUIGzZ8/i5ORElSpVtA7N5sgsUKkXFRVlGjwL0LFjRz799FPq1aunYVTW5dq1a1SuXJmcOXPi7+9Phw4dyJ49e6L9JHfm0kverKm9c/78eYoUKcLjx485dOgQo0eP5osvvqBFixaaxpVRydlGCCuUcHYovR7+G8uX7l29epWOHTvyww8/sH37drp160bz5s1ZuHAh0dHRWodnE2JjYwH4448/uHz5shQVqRQSEsLw4cMBePr0KWDomhcUFGTKrTA8mOzw4cO8//77rFq1iqlTp3LhwgXTduMV48GDB7N06VLJ3X8kb2krJCSE2bNnA4aZ8erWrcuPP/5IzZo1AZkhSgtyxhHCCtWuDblyxS+vWqVdLG+CccDxyZMnqVWrFosWLWLDhg00bdqUXbt2MXv2bNOJQiTPOC3lV199RY0aNThx4oTGEdkOpRT58uUjKiqK4OBg04BZX19fypUrJ1N+JhAXF4eXlxfDhg2jffv2/O9//+Pdd99l06ZNAGTOnBmlFMWKFcPHx0dy9x/JW9oxHq/3798nODjYtN7Hx4ds2bIB0m1WC1JYCGGFHBzgv2ebAbBjB9y+rV08r5vxqvratWvx9PTEyckJX19f8uTJw+TJk00nCZGYsSj766+/uHv3rmn9ggULaNOmDTNnztQoMutnzF3Cq5p58+blrbfeYuHChezdu5crV66wYsUKevXqpVWYVsOYL4gf++Tu7k6vXr04fvw4RYsWpWXLljRp0oSnT5+i0+lwd3dnzJgxQMa9eix5S1svO17nzp1r9n+heLOksBDCSrVrF/86Lg5+/VW7WF4344m3Y8eOHDx40LR+8uTJBAcH8+jRI61Cs3rGouz7779n27ZtgOHE6+vry4QJE+jSpQtg3rgRBsbcGXMTEBBA586dOXToEEuXLmXgwIF07NiRc+fOmY25yKiM+SpTpgxBQUGmWY2UUri4uODv788777xD+/btcXJyQq/XmzUCM+rVY8lb2jLmY9KkSXTq1CnR8Xrjxg1y5MihcZQZl9xjE8JK1a0LOXKA8cLLqlXQo4emIb02dnZ2hIWFcezYMbZt28bNmzcpWbIkT58+xcfHhyxZssgg5Bcw5uWTTz5hy5YttG/f3nTizZo1K7Vr1waQ3CXw448/MnjwYH788UeaNWtmuor84Ycf0rRpUypWrMiKFSu4fPkyn3/+uWmshUzfC5GRkdSsWZPvv/+eY8eO0b17d3x9fVFK0aRJE5YvX276zOl0ugyfLyPJW9rr3r074eHhVKpUSY5XKyJnGiGsVKZM0Lp1/PK2bfFFRnpz+fJlOnXqxOPHj/H29uZ///sfV65cwdnZmf/7v//TOjyrZizKQkND2bp1Kz179mT27NlMmDCB74xPWhRmmjRpQuvWrRkzZgzDhg3j1KlTKKUoVKgQFSpUAKBFixasX7+eY8eOmcZaSCPF0IVnwoQJdOrUif3799OvXz+2bduGTqdj+/btBAcHU7RoUUDylZDkLW0Y7y4ePnyYe/fuUalSJcBwvP76669yvFoBTQuL3bt307JlSzw8PNDpdKbH2BsppRg7diweHh44OztTt25dTp48abZPdHQ0/fv3J1euXGTJkoVWrVpx7do1s33u3btH586dcXd3x93dnc6dO3P//n2zfa5cuULLli3JkiULuXLlYsCAATIPstBcwu5QsbGwfr12sbwOxtv93t7e7N69my+++IIVK1bQp08f2rZty6BBg/Dw8ADkivuLJFWUXbhwATc3N1q1agVIH+3n5cmTh4CAALp27crBgwf5/PPP2bdvHxDfGMmSJQu9e/dm0aJFWoZqlXLlykWfPn345Zdf8PHxoVmzZhQqVIh+/frRoUMHIP7hbyKe5O3VGc8DW7ZsYc6cOYBhJrwsWbLg7+9PYGCgluEJAKWh3377TX355Zdq9erVClBr16412z558mTl5uamVq9erY4fP678/PxU/vz51YMHD0z79O7dWxUoUEBt3bpVHTlyRNWrV0+VK1dOxcbGmvZp2rSpKlOmjNq3b5/at2+fKlOmjGrRooVpe2xsrCpTpoyqV6+eOnLkiNq6davy8PBQ/v7+qfp9IiMjFaAiIyMtS4gQz4mOVipbNqXA8NW8udYRpb0nT56okSNHKqWUiouLU0op9ddff6kOHTqYHcfCnF6vN71+/Pix6fXw4cPVnj17tAjJ6hlzljB3x44dU82bN1f29vaqX79+6syZM0qp+M9iRESE2bIw5C9hDo8ePaomTJigfvnlFxUTE2PaR5hLT3l7k+0dvV6vYmNjzXLz4MEDNXLkSHXp0iXTuqdPn6q7d+8qpeR41ZKmhUVCzxcWer1e5cuXT02ePNm07unTp8rd3V3Nnz9fKaXU/fv3VaZMmdTy5ctN+/z777/Kzs5Obd68WSml1KlTpxSgDhw4YNpn//79ClD//POPUspQ4NjZ2al///3XtM8vv/yiHB0dU3XQSGEhXoeuXeMLi0yZlLp3T+uI0tbVq1fV+++/rzZt2mS2vkOHDioqKkqjqGyDFGUpY2yQPHnyRB09ejTJfZo2bap0Op0aN26c2feIF5PGm2XSQ95eZ3vHeOw9fvxYbd++PdH6nTt3qu+++069/fbbqnLlymrgwIFq8ODBqk+fPmkei0g9q+1bcPHiRSIiImjcuLFpnaOjI3Xq1DHdsg4JCeHZs2dm+3h4eFCmTBnTPvv378fd3Z2qVaua9nnnnXdwd3c326dMmTKmLhdg6IMbHR1NSEjIC2OMjo7mwYMHZl9CpLWED8t79gw2bNAultfB09OTTp06sW/fPsLDwwHD8yyCgoLIkiWLdONJxu3btzl9+jS//fabqYuA8XkfxgGMIr5706ZNm6hYsSJDhgzh1KlTQHw3sYEDBzJgwAD69u1rtl68WMLuiZKvlJO8Jc94vO7YsYOGDRvSoUMHDh06ZFo/ffp09Ho9fn5+hIeHkzdvXsqXL89nn30GyAx4WrPawiIiIgIwzFGcUN68eU3bIiIiyJw5M9mzZ092nzx58iR6/zx58pjt8/zPyZ49O5kzZzbtk5SAgADTuA13d3e8vLxS+VsK8XKNGkHWrPHL6elheQcPHmTQoEEEBgYyZ84c2rRpQ7Vq1Rg7dixXrlwBZABecqQoS53GjRszf/58Tp8+zfjx49mxY4fp82VnZ8edO3fImTOnaVmknDTmLCN5e7HmzZuzZcsW9Ho948eP5+effwbg119/pW/fvgwZMoRvv/2WBw8e0KVLFwoXLgzIsas1q8/+840KlYLpw57fJ6n9LdnneSNGjCAyMtL0dfXq1WTjEsISjo7w3xhcAP74A9LLzTFnZ2eKFCnCggULGDNmDGXLlmX16tXMmjWLYsWKaR2eVZOiLPXc3Nz49NNP+fzzz7l37x6tWrVi6NChHDhwgMGDB1OtWjVABtBawjhlr0gdyduLKaVo2LAhI0aMIF++fAwdOpSPP/6Yu3fvmvLWunVrwsLCOHPmjMbRCiOrLSzy5csHkOiOwc2bN013F/Lly0dMTAz37t1Ldp8bN24kev9bt26Z7fP8z7l37x7Pnj1LdCcjIUdHR7JmzWr2JcTrkLA7VHQ0bNyoXSyvyngVPSoqCnd3d/r370/BggUZNGgQ586d48yZM2bdEkXSpChLPePFonr16rF582YWL17M77//TuPGjSlTpoypG5Q09l7OeBxfvXqVNm3acOHCBY0jsg2St5QzXhipUKECgYGBLF68mE2bNlGoUCEmTpwIGI7VKlWqSNdPK2K1hUWhQoXIly8fW7duNa2LiYlh165dVK9eHYBKlSqRKVMms32uX7/OiRMnTPtUq1aNyMhIDh06ZNrn4MGDREZGmu1z4sQJrl+/btpny5YtODo6muZIFkJLTZpAwgf/2nJ3qIT93Y0nB+NJISAggPnz52sWm7WTouzVGD97xjsSbdu25dSpU5w/f940TaV0TUkZYy4vX77M8ePH+frrr3n27JnGUVk/ydvLhYeH4+fnx/nz503rlFI0bdqUfv36UbBgQXx8fEzbRowYQdmyZTWIVCTpjQ8XT+Dhw4fq6NGj6ujRowpQ06dPV0ePHlWXL19WShmmm3V3d1dr1qxRx48fVx9++GGS0816enqqbdu2qSNHjqj69esnOd1s2bJl1f79+9X+/fuVr69vktPNNmjQQB05ckRt27ZNeXp6ynSzwqp8+GH87FBOTko9fKh1RK/m/v37atCgQSosLMy07unTp2rXrl1KqfQxc8rrsnz5cvXpp58qpQwzHSml1IEDB5Sfn5+WYYl0zjgrz4MHD1RwcLBSSqlnz56p4OBgVb16ddM6YS495u11tnf+/vtv5evrq+rXr68WLlyo7ty5Y9p29epV1atXL9PPldnvrI+mhcWff/6pgERfXbt2VUoZDsYxY8aofPnyKUdHR1W7dm11/Phxs/d48uSJ8vf3Vzly5FDOzs6qRYsW6sqVK2b73LlzR3300UfKzc1Nubm5qY8++kjde27OzsuXL6vmzZsrZ2dnlSNHDuXv76+ePn2aqt9HCgvxOq1eHV9YgFIJZlm2Cca5yJUyHJNHjx5VnTp1Uu+++64KCgpSa9asUZ999pnZMxlE0qQoE1pavny5aWrPpJ4NIpKWnvL2uts7R48eVf7+/qpixYpq8ODBpudTfP/996p169ZKKdvNXXqnU0qmDUkrDx48wN3dncjISBlvIdLckyeQOzc8emRYbtsWVq7UNqbkqCQmP9Dr9djZ2dGpUyecnJzIlCkTS5cupW7duuTPn593332XNm3apGiShoxEKYVer8fe3p67d+9y5coVpk2bxt27d2nXrp1p+uxx48bh7Oysdbg2w/g5i42NxcHBQetwbEZkZCRjx46lT58+Mp4nFdJT3t5Ee0ev1/Pzzz8zatQoHj16RIMGDdi7dy/Tpk3Dz8+PuLg4GQ9lhaSwSENSWIjXzc8P/vc/w2sXF7h5E7Jk0TamFzE22saOHcuZM2f47rvvyJEjB2A4YTx8+BB3d3e2bt1KSEgIw4cPl4LiP1KUvV7GBsn9+/fJli2b1uHYBGPOXlTY7ty5k0mTJpHFWv9D0kh6zdvrbu8k/H9Mr9czb948duzYQePGjenVq1ea/zyRdqSwSENSWIjXbeVKaN8+fnnVKvjgA+3iSYmNGzcyevRo3Nzc6N69O61atSJHjhzo9XqUUjx+/Bg/Pz+mT59OiRIltA7XKkhR9mb07t2bpUuX8ueff5o9RFW8WKdOnXB0dCRz5sxmhW3z5s15//335XP4Auktb2+qvfOiuxK2lq+MRAqLNCSFhXjdHj0ydId68sSw7OcHy5drG1NKhIaGsmjRIo4cOUL9+vUZNGiQqaEMMGPGDBwcHOjfv7+GUVofKcrS3rNnz8iUKRNxcXEcPHiQn3/+mZiYGObMmUPmzJm1Ds/qKaW4d+8eOXLkkMI2FdJb3rRo79hqrjIaq51uVgiRWJYs8O678csbN8YXGdZKKUX58uWZNGkSH330EQsWLKBGjRqsWbPGtM/gwYPp16+faX9h0KJFC4KCgihfvjzff/89M2fO5O7du9jZ2WFvb4+bmxuNGjUym3JbJM04jezy5cv5+++/sbe3p3r16gwbNoyCBQtKUZEMY+6WLFnCqVOnyJEjB7Gxsbzzzjvs2rWLf/75Rxp8SZC8pS3JlW2QwkIIG5PwYXmPHhmexG1NjCfTB/89Htx4MnB1daVPnz4EBwdz69Yt2rZty7Jly0zfZ2dnZ7a/kKIsLRk/X2fOnCEoKMi03sfHh9GjR2sVlk0w5u7cuXMsXLgQAAcHB9zc3GjVqhVbtmzRMjyrJXlLO/J8GdshhYUQNqZ5c3B0jF+2tpmhjCfTiRMnkj17dpYtW8bjx48Bw8nB09OTQYMGMXLkSNq0aQNIg9hIirLXI+Hny9/fH2dnZ0JDQwF49OgRly9f1igy6xYXF2d6mCAYcufi4mLK3dOnT2nfvr2pC6McxwaSt7SRMC/G/+OE9ZP59YSwMW5u0KwZrFtnWN6wAZ4+BScnTcNKpEePHjx58oTvv/+esLAwunfvjre3NwB58uQhIiICJycn6TebQMKibOHChcyZM4fWrVvj4uJiVpQ9ffrUrCiT/CXPmJ/Jkydz7do1tmzZwooVKyhatChubm5UrFiRkSNHahyltoyzjgUHB7Njxw4+//xzs0GzU6ZM4erVq2a5y5IlCxUrVmTUqFFAxixsJW+vT1xcHA4ODowcOZI6derQpEkTrUMSKSCFhRA2qG3b+MLi4UPYuhVattQ0pESKFSvGyJEj+emnn5g3bx4bN25k0qRJ+Pr68tVXX/HTTz8B0jBOihRlae/x48fcvXuXDh060LhxYwICAvjyyy9xd3enaNGiQMb+LBqL2kuXLjF27Fi2bdtG//79adSoEUopwsPDJXdJkLy9Pg4ODkRHR7N37142b95MwYIFKVmypNZhiZeQWaHSkMwKJd6UyEjIkwdiYgzLnTvD0qXaxpSc27dvM3z4cBYtWkSePHmoWrUqv/76q9ZhWbWIiAhTUZYtWzZTUValShV++ukn6tWrZ7paKl7M2HALCwujePHipvXz5s3j+PHjzJ07V8PorNPBgwcJDAzkypUrtGrVCn9/f7Pt8+fP59ixY5K752SkvL2u9o7xeA0JCaF8+fLY2dlx9epVpk6dir29PbNmzUqznyVekzR+kneG9rofcS9EQi1aKAWGL3d3pZ4+1TqixPR6vYqLizMtR0REqIMHD6pbt24ppZTZNpG0W7duqR49eiidTqfy5s2rWrVqpXVINufhw4eqbdu26uLFi0oppWJiYpRSSjVp0kSFhYVpGJn10ev1SimlLl26pEaMGKHc3d1VnTp11N9//22237vvvqvOnDmjRYhWKaPl7XW2dx4+fKg6dOigLl26ZLb+33//VUrJecPayaUuIWxUu3bxryMjYft27WJ5EZ1OZ3ZFPW/evLz99tvkypULkAF5yVFKodfryZUrF99//z3Xr19n/fr1/PDDD4DMkpKchLmJi4vD1dWV6tWrs3HjRgAyZcoEwKBBgyhSpIgmMVorY7ccb29vJk2axL59+wgLC6N8+fJMnDjRtF/Hjh0pVKiQVmFaHcmb5ZI6Xt9++202bNhgtp+Hhwcg5w1rJ38dIWxUy5bwX/sIMDyFW6QfUpSlnHEGnsjISMA8N/b29jx69AhHR0c2bdrEggUL2LlzJ3PmzCF37txJPtU3I9qwYQPFixdn165dpnV6vZ5SpUrRq1cv2rRpQ6tWrUzbPvroIxwcZJim5C31LDlev/32W3bv3q1JvCJ15MwkhI3Knh0aNoxfXrcOnj3TLJwU0ev1ZtMwCpEWjMVB69atqVatmmlaT4CdO3fStGlTTp8+zaVLl/jqq6+YOnUqN2/eNF0BFVC4cGFKlSrFV199xTfffEN4eLipwffee++RO3dufH19NY7S+kjeUs+S4/Xu3bsUK1ZMo4hFasjg7TQkg7fFmxYUBN27xy9v3gzWNiNfXFwc9vb2PH36FCdrmxPXhuj1epRScoU9GcHBwQQEBHDz5k3atm1Lu3bt8PDw4Pjx4/j6+nLnzh3+7//+j2+++QYHBwdTlyhhcPbsWYKCgli7di2lSpViypQpvPXWW3z44YcULlyYiRMnyoQBSciIeUuL9o4cr+mTFBZpSAoL8abdvQt580JsrGH5k08gMFDbmF6kVq1auLq68tNPP5EzZ06tw7EJUpSlnEow+9PixYvZvHkz5cuXZ9y4cXh5eZn2+7//+z+cnZ0ZMWKEhtFat4MHDzJs2DD27NlD8eLFefz4Mfv378fDw0OmR01GRsrbq7Z35HhNv6SwSENSWAgtNGkCW7YYXufMCdevm4+9eNOMJ4xbt26RI0cOdDodSikCAwPZsGEDH3zwAd0T3mYRLyVFWer98ccfDBkyhIiICAYMGMDgwYPJmjUrt27don379vz5559ah2h1nr+qvnv3bo4cOUK1atWoWrVqurvqnlYyYt7Sur0jx2v6IYVFGpLCQmghMBB69oxf3rrVfOyFVkaNGkWnTp0oUaIEYDj5rlmzhsOHDzNlyhSNo7NOUpSljvGOzoEDByhQoIDZlU4wPLSscuXKuLi4sH79esqXLw/Aw4cPcXNzS5cNvrQgebFMRsqbJe0dOV4zBvkLCWHj3n8fEna713p2KOM0qU5OTixN8NQ+Ozs72rZtayoq5JpGYsauErNmzeLs2bPY2dlhb29Pz5496datG2FhYRpHaF2M401mzpxJhQoVmD59OhEREYChkefj48Nnn33GyJEjTY0UpRRubm6AzKz1IgnzIsdpyknekifHa8YgfyUhbFyuXFC3bvzymjWg1cRLxivudnZ2NG3alDt37rBp0yZu3rzJ0aNH+euvv0z7pqf+xmlFijLLTJ06lf79+7N+/XomTZpEWFiYqRFy7tw50xXV9NbP/XUzThggUkfyljw5XtM36QqVhqQrlNDK/PnQp0/88o4dUK+eNrE8fPiQbt264eLiwt69e3nw4AHe3t4UKVKETp060bJlSzlhJCFhToKDgwkMDKRVq1ZUqVKFf//9l6ioKGrVqqVxlNbryZMnrF69moCAAO7evcuQIUPw9PSkR48e3Lp1iyxZsmgdok2QY9MyGS1vr9rekeM1/ZLCIg1JYSG0cuMGeHiA8QGmffvCnDnaxPLs2TNWr15NgwYNuHHjBgEBAcyaNQudTicDj19CirK0MXXqVGbMmEGuXLn4+OOPGTp0qKl/t0iesR/7119/zblz55g2bRqurq5ah2X1Mlre0rK9I8dr+iKFRRqSwkJoqV492LnT8DpfPrh2zXzsxZvyxx9/0CTBwzQCAgIIDw9n9uzZbz4YGyNF2at5vjFy79493NzccHBwkIIshYx5mjp1KkuXLmXo0KF069ZN67CsXkbLW1q0d+R4TZ9kjIUQ6UTbtvGvIyJg7943H8PZs2eZMWMGZ8+eNa0bMWIEISEhXL9+/c0HZGN27NhBhw4dyJ07N2XKlKFMmTJ89dVXUlSk0PNXOLNnz46DgwMgY3pSypgnf39/Ro8ezc8//6xxRLZB8pZ6crymT1JYCJFOtGkDCf8vfhOzQyW84RkXF0fRokVp06YN+/fvN1v/008/kT9//tcfkA2Tokxo7ddffzW9dnFxwc/Pj40bNwKGrj4iaZI3IeJJYSFEOpE/P9SsGb+8enX8mIu0FPfflFPP36qOiori0KFDhISE8PXXXzNixAi+/vprevXqRVRUVNoHYuOkKBPW5OzZsyxYsIBz584BEBsba5qhDGSqzxeRvAlhTj7xQqQjCbtDhYdDgjZqmjHevq5YsSKTJ082FRrr169n/PjxlChRAhcXF7Zv305UVBQtW7akbNmyaR+IDZKiLG3FJTGvsvEKsVJKpvxMRsLcGQvb1q1bs2/fPgAcHBySzG9GJ3mznByvGYMUFkKkI23amC+/ru5Qer0ePz8/VqxYQfv27fnjjz/o3LkzK1asYPDgwezevZuSJUsyaNAg3nvvvdcThA2SouzVxMbGAnDkyBEgPp8JGyV2dnb8/fff6HQ66aedwItyB/Do0SMOHjyYqLDt27cvISEhmsRrLSRvlpPjNWOSwkKIdMTTE6pXj19eter1dIeys7Nj2LBhTJo0CQcHByZMmEBgYCAuLi4AODo64uvry+LFi9P+h9s4KcosZxzY+f7771OiRAnWrVtHXFycWaPkt99+o0KFCqaHCQqDF+UOYPXq1UycOJGSJUuaFbYtWrSgUqVKWoatOcmb5eR4zaCUSDORkZEKUJGRkVqHIjKw6dOVgvivAwde788LDw9Xw4YNU5kzZ1ZNmzZVu3btUkopdeLECdWzZ8/X+8NtVGxsrPrtt99U+/btVc2aNdXChQvNtn/99ddq+vTpGkVn/S5evKj69++v6tSpo8aNG6cuXbqklFJKr9era9euqQULFqjatWurU6dOaRyp9XlR7oznrejoaNWlSxd19+5dLcO0OpI3c6lp78jxmrFIYZGGpLAQ1uDyZfPCYujQtHnfZ8+eKaWUCgsLS3L7ihUrlE6nU4ULF1Y3btww2xYXF5c2QaQzUpSlnvGzdO3aNTV58mRVokQJVb9+fXX48GGz/ZYvX64ePXqkRYhWK6nc1atXTx08eNBsPylszUneEktpe0eO14xHCos0JIWFsBZVq8YXFt7eSun1affeBQsWVLVq1VLbt283Kxhu3bqlPvnkE7Vy5UqlVHwhIqQoe1X6ZD7A586dU40aNVI6nU61a9dOXbt27Q1GZv1Sk7vw8HCllFKnTp3K8IWt5C15ybV35HjN2GSMhRDpUMLZoS5fhrQaR6iUIigoCC8vL6ZMmcKiRYt49OgRALly5SIyMpKCBQsCiR9+lJEZ+xo3atSI2rVrs2PHDrP57evXr0+PHj2YMmUKefLkMQ16BJmuEuKn5v3666/ZunWr2bYiRYowc+ZMPD09efLkCdmzZ9ciRKuVmty5u7sDULJkSRYsWABk3OcwSN4sJ8drxiZnLCHSoQ8+MF9euTJt3len01G/fn1Gjx5N8eLFGT58OB07diQkJIRffvmF3bt38/bbb5v2FfGkKLOcnZ0dSimOHDlC3759GTlyJGFhYYBhCstSpUrRo0cPunXrhouLi0xbmUBqc6fX683yl1ELW8mb5eR4zeC0uE2SXklXKGFNKlWK7w5VuLBl3aEWLlyoevfuberK83z3pmPHjqmGDRsqnU6nSpcurb777rsk9xPxTp8+rfr3769y5sypWrVqpYKDg9WyZctU3rx5tQ7NKiXsVhEREaGmTZumatWqpfz8/FRwcLBp29tvv62OHj2a6HsyMsmdZSRvL/ei9o7kTuiUklIxrTx48AB3d3ciIyPJmjWr1uGIDG7KFBg+PH75yBGoUCF173HixAmGDh1Kly5d+OijjwDDlffY2Fjs7e1NV6YiIiJ4+PAhb731lmldRr9jERgYyJEjR5g9ezYODg7ExsaaukQBHD9+nCFDhrB9+3ZKlSpFnz596NevX6L9MhrjZ+fRo0c8fPiQfPnyJdpn9+7djBkzhv3799OyZUvc3d3ZtWsXZ8+e1SBi6yG5s4zkLfWM7Z379+/j7u4uuRMmGfdenRDp3PPdoSx5WF6ZMmX44IMP6NevH59++ilXrlxBp9ORKVMm7OzsiIuLQylF/vz5KVasmOn2f0YvKgCqVavGhQsXWLFiBWAYZ6GU4tmzZ+j1enx9fdmyZQv//vsva9asoU+fPoB0gzJ+dn788UeqVavGyZMnAcPDtoxjT2rXrs2ff/7JL7/8wt27d3n8+LHpmSkJx6dkNJI7y0jeLCe5E8+TOxZpSO5YCGtToQKEhhpeFy0KYWFgSZt/3759jBw5EicnJ5o2bUqjRo0oXbp0msaaHi1cuJDPP/+cdu3aMXr0aNMYCsD0oKiM3Bf7ZT777DMiIiKYNWsWOXPmBDAVZ5kyZZICNhmSO8tI3lImqfaO5E6A3LEQIl1r1y7+9dmzcPx46r5fKYVer6d69eqMHTuWfPnysWzZMoYMGULHjh0ZM2YMe/bsYenSpaxYsYLw8PC0/QVsXM+ePfntt984e/YsPXv2ZObMmaYresauZOLFevXqRXh4OGXLluXbb7/lxo0b6HQ6MmfOjE6nIzY21jTwMyPPwpMUyZ1lJG+Wk9wJkDsWaUruWAhrc+YMFC8evzx6NIwbZ/n76fV6tmzZws6dO7l37x4nT57k2rVrVK1aFVdXV/r27UulSpVePfB0QBmeE4SdnR07d+5k8eLFnDp1iuzZs5MzZ06KFi1Ko0aNuHDhAo6OjtSqVQsPDw+tw7Y6z549Y/z48fz5558UKFCAmjVr8vbbb5tmH4uOjsbR0VHjKK2T5M4ykreXe1F7R3InpLBIQ1JYCGtUtmz8nYoSJeDUKcu6Q+n1erMr7FFRUbi6ugLw8OFD3Nzc0iLcdEuKstQzfuYiIyNZt24da9as4datWyiluHnzJpUrVyZLlizcunWL4cOHU6NGDa1DthqSO8tI3lImqfaO5E6AFBZpSgoLYY3GjYMxY+KXT5yAVxkekbDAkNmfUkaKsrSh1+vZu3cvV69eJSYmhtDQUIoVK4a3tze+vr5mY1iEOcmdZSRvSUtJe0dylzFJYZGGpLAQ1uj0aShVKn557FjzQkO8OVKUWSZhtzKROpI7y0jekpdce0dyl7HJX12IdK5kSfPCwuwp3FeuvPF40g0LcpfwRJthiwoL8pZw9izjtTDjxAIZhoXHquQO+cxZSj5zwgJSWAiRAbRtG//65EnDXQz27QNvb9i/X7O4bJbkzjJpkDdjQZahpupNo8+b5M4ykjfLZcjcZXDyVxYiA0g47SzA6tXAfw9uM/0rUk5yZxnJm2Ukb5aT3FlG8iYsJGMs0pCMsRDWSilDl6iwMMNyuXKK0BseEBEB+fPDv/9aNlVURqTXQ4ECkrvUkrxZRvJmOcmdZVKQN2nviBeROxZCZAA6nXl3qL//1nE2wjArEdevw+HD2gRmiw4fNpxwQXKXGpI3y0jeLCe5s4zkTbwCKSyEyCASFhYAq+zaG144OMCqVW8+IFu1apUhZyC5Sw3Jm2Ukb5aT3FlG8iZegXSFSkNya1BYM6WgWDE4d86wXJEQQqhsWChYEC5dkm4CL6OUYUDj1avx6yR3Lyd5s4zkzXKSO8ukMG/S3hEvIncshMggnu8OdYRKXKCQYeHKFfj7b20CsyWhoeYnXJDcpYTkzTKSN8tJ7iwjeROvSAoLITKQRN2h+G+Fvf1/U0WJZK1ebchVQpK7l5O8WUbyZjnJnWUkb+IVSVeo58ydO5evv/6a69evU7p0aWbOnEmtWrVS9L1ya1BobvZsw6O1X3BYKwVF7gdzkcIAVOEQh6hq2JgpE7i6vvi9dTr46ivw90/joK3ES3IHQFQUPHuWeH1Gzp3kzTKSN8tJ7iyThnl7oBTu9+9Le0ck4qB1ANZkxYoVDBo0iLlz51KjRg0WLFhAs2bNOHXqFAULFtQ6PCFernp1wwngxo0kN+uAtqziaz4H4DBvcwlvfLhsOJncu/fi986XD6pVew1BW4mX5C5ZGTl3kjfLSN4sJ7mzzOvMmxD/kTsWCVStWpWKFSsyb94807qSJUvSunVrAgICXvr9csdCWIU7d6BbN9iwIcnNh6hCVQ6ZlqtyAC+uJrmvkc4jP1SpApkdE297yTjIV9n+Ot87ye3R0YYnzv57DR3J/9f4su0mnl5Qozo4OqVsf1sU/RT27oNryX+OIHHeEi6bvS7oBbVqgaOT2d8pJa9Tut+rvE6T93ryBDb/ju78OfP3fVFOknidaF3RYtCyBToXlxTFkpq4rem9dI8fwc8/w8kTqcrRi+hQ4FsWOnd+4R2NlIz5tvp9oqIgKAj+DjVsf0lusvAIP/6XaP0DwB2kvSMSkcLiPzExMbi4uLBy5Uref/990/qBAwcSGhrKrl27En1PdHQ00dHRpuUHDx7g5eUlB5rQnlIwbx4MGmR42FFcXPwmwIdLXMFbs/CEEEJYPy+uJHmueGBnh7teL+0dkYgM3v7P7du3iYuLI2/evGbr8+bNS4TxQTHPCQgIwN3d3fTl5eX1JkIV4uV0OujbF0JCoEgRsIs/1HXAR/ysXWxCCCFsl50dFC6sdRTCSskYi+fonrtvqJRKtM5oxIgRDBkyxLRsvGMhhNXw9YWjR2HwYFi40FBwKMUoJnCXHBzi7aS/L1culIenWUHyvJfd67Tm7Sn+XqU39Ee+f998O8n0QciWDfLkAV0Gvm6j9HDjJkTeN1/9XN4SLqus2SB3LtDZmf19Uvvaku+xlp8BCvX4CcTE/LeUID9JvFboDA8wc3REqQTbX/KzLI0xuXyIDOC/8weffmoYBJ4/v9YRCSskhcV/cuXKhb29faK7Ezdv3kx0F8PI0dERR8fEfc6FsCouLrBgAezZA6dPG1bxhPn0SbyvTgclS8LJk284SGtlB+SH0g0NuUuuNSW5S8AOyAelG0jeUkUHuEDpKjabtzdRyCT7ukoV1D9h8N/YgaQvAuigeHE4dCjRe7xIetzHbHutWnDmDDw35sIOvfk3lSoF8+fDgwcvD0BkSDLGIoGqVatSqVIl5s6da1pXqlQp3nvvPRm8LWzb5cvg45O6/WUmNAPJnWUkb5aRvFlOcmcZC/L2IFs2ae+IJGXge/WJDRkyhO+//55FixZx+vRpBg8ezJUrV+jdu7fWoQnxatasSbZbkxk7O8P+wkByZxnJm2Ukb5aT3FlG8ibSkBQWCfj5+TFz5kzGjRtH+fLl2b17N7/99hve3jJ7jrBxK1Ykvi9ufLrq809ZVcqwvzCQ3FlG8mYZyZvlJHeWkbyJNCSFxXP69u3LpUuXiI6OJiQkhNq1a2sdkhCvJjwcDh5MfOJo0QLOnjX8m5BShv2vX39zMVoryZ1lJG+WkbxZTnJnGUvz9oLZMoWQwkKI9G7t2vinJNnbQ+bMhsF3a9fCW28Z/p03z7DewcH8+zI6yZ1lJG+WkbxZTnJnGUvz9oIHsAohhYUQ6d3//me4ymRnB0WLwpEj0KuX+SNue/c2f+aFUobvy+gkd5aRvFlG8mY5yZ1lLM1bRi/IxIspkWYiIyMVoCIjI7UORQiDmzeVsrNTCpTq00epx4+T3//xY6V69zbsb2en1K1bbyZOayS5s4zkzTKSN8tJ7izzCnmL1OmkvSOSJHcshEjPoqIM846vWQNz54Kzc/L7OzsbbnuvWWP4vocP30yc1khyZxnJm2Ukb5aT3FnmVfJWosSbiVHYHHmORRqS51gIIYQQIr2T9o54EbljIYQQQgghhHhlDi/fRaSU8ebPA3nUvRBCCCHSKWM7Rzq9iOdJYZGGHv7XT9PLy0vjSIQQQgghXq+HDx/i7u6udRjCisgYizSk1+sJDw/Hzc0NnXGqthR68OABXl5eXL16VforppLkzjKSN8tJ7iwjebOM5M1ykjvLvCxvSikePnyIh4cHdnbSq17EkzsWacjOzg5PT89Xeo+sWbPKf34WktxZRvJmOcmdZSRvlpG8WU5yZ5nk8iZ3KkRSpMwUQgghhBBCvDIpLIQQQgghhBCvTAoLK+Ho6MiYMWNwdHTUOhSbI7mzjOTNcpI7y0jeLCN5s5zkzjKSN2EpGbwthBBCCCGEeGVyx0IIIYQQQgjxyqSwEEIIIYQQQrwyKSyEEEIIIYQQr0wKC43t3r2bli1b4uHhgU6nY926dVqHZBMCAgKoUqUKbm5u5MmTh9atWxMWFqZ1WDZh3rx5lC1b1jQ/ebVq1fj999+1DsvmBAQEoNPpGDRokNahWL2xY8ei0+nMvvLly6d1WDbh33//pVOnTuTMmRMXFxfKly9PSEiI1mFZNR8fn0SfN51OR79+/bQOzerFxsYyatQoChUqhLOzM4ULF2bcuHHo9XqtQxM2Qh6Qp7FHjx5Rrlw5unXrxgcffKB1ODZj165d9OvXjypVqhAbG8uXX35J48aNOXXqFFmyZNE6PKvm6enJ5MmTeeuttwBYsmQJ7733HkePHqV06dIaR2cbDh8+zMKFCylbtqzWodiM0qVLs23bNtOyvb29htHYhnv37lGjRg3q1avH77//Tp48eTh//jzZsmXTOjSrdvjwYeLi4kzLJ06coFGjRrRr107DqGzDlClTmD9/PkuWLKF06dIEBwfTrVs33N3dGThwoNbhCRsghYXGmjVrRrNmzbQOw+Zs3rzZbDkoKIg8efIQEhJC7dq1NYrKNrRs2dJseeLEicybN48DBw5IYZECUVFRfPTRRwQGBjJhwgStw7EZDg4OcpcilaZMmYKXlxdBQUGmdT4+PtoFZCNy585ttjx58mSKFClCnTp1NIrIduzfv5/33nuP5s2bA4bP2y+//EJwcLDGkQlbIV2hRLoQGRkJQI4cOTSOxLbExcWxfPlyHj16RLVq1bQOxyb069eP5s2b07BhQ61DsSlnz57Fw8ODQoUK0aFDBy5cuKB1SFZv/fr1VK5cmXbt2pEnTx4qVKhAYGCg1mHZlJiYGH766Se6d++OTqfTOhyrV7NmTbZv386ZM2cA+Pvvv9mzZw/vvvuuxpEJWyF3LITNU0oxZMgQatasSZkyZbQOxyYcP36catWq8fTpU1xdXVm7di2lSpXSOiyrt3z5co4cOcLhw4e1DsWmVK1alaVLl1KsWDFu3LjBhAkTqF69OidPniRnzpxah2e1Lly4wLx58xgyZAgjR47k0KFDDBgwAEdHR7p06aJ1eDZh3bp13L9/n48//ljrUGzCF198QWRkJCVKlMDe3p64uDgmTpzIhx9+qHVowkZIYSFsnr+/P8eOHWPPnj1ah2IzihcvTmhoKPfv32f16tV07dqVXbt2SXGRjKtXrzJw4EC2bNmCk5OT1uHYlITdPX19falWrRpFihRhyZIlDBkyRMPIrJter6dy5cpMmjQJgAoVKnDy5EnmzZsnhUUK/fDDDzRr1gwPDw+tQ7EJK1as4KeffmLZsmWULl2a0NBQBg0ahIeHB127dtU6PGEDpLAQNq1///6sX7+e3bt34+npqXU4NiNz5symwduVK1fm8OHDzJo1iwULFmgcmfUKCQnh5s2bVKpUybQuLi6O3bt389133xEdHS0DklMoS5Ys+Pr6cvbsWa1DsWr58+dPVOyXLFmS1atXaxSRbbl8+TLbtm1jzZo1WodiM4YNG8bw4cPp0KEDYLgQcPnyZQICAqSwECkihYWwSUop+vfvz9q1a9m5cyeFChXSOiSbppQiOjpa6zCsWoMGDTh+/LjZum7dulGiRAm++OILKSpSITo6mtOnT1OrVi2tQ7FqNWrUSDSN9pkzZ/D29tYoIttinNTDOBBZvNzjx4+xszMffmtvby/TzYoUk8JCY1FRUZw7d860fPHiRUJDQ8mRIwcFCxbUMDLr1q9fP5YtW8avv/6Km5sbERERALi7u+Ps7KxxdNZt5MiRNGvWDC8vLx4+fMjy5cvZuXNnopm2hDk3N7dEY3iyZMlCzpw5ZWzPS3z22We0bNmSggULcvPmTSZMmMCDBw/kCuhLDB48mOrVqzNp0iTat2/PoUOHWLhwIQsXLtQ6NKun1+sJCgqia9euODhIUyelWrZsycSJEylYsCClS5fm6NGjTJ8+ne7du2sdmrAVSmjqzz//VECir65du2odmlVLKmeACgoK0jo0q9e9e3fl7e2tMmfOrHLnzq0aNGigtmzZonVYNqlOnTpq4MCBWodh9fz8/FT+/PlVpkyZlIeHh2rTpo06efKk1mHZhA0bNqgyZcooR0dHVaJECbVw4UKtQ7IJf/zxhwJUWFiY1qHYlAcPHqiBAweqggULKicnJ1W4cGH15ZdfqujoaK1DEzZCp5RS2pQ0QgghhBBCiPRCnmMhhBBCCCGEeGVSWAghhBBCCCFemRQWQgghhBBCiFcmhYUQQgghhBDilUlhIYQQQgghhHhlUlgIIYQQQgghXpkUFkIIIYQQQohXJoWFEEIIIYQQ4pVJYSGEyHAuXbqETqcjNDRU61BM/vnnH9555x2cnJwoX7681uFYLR8fH2bOnKl1GEIIIZIghYUQ4o37+OOP0el0TJ482Wz9unXr0Ol0GkWlrTFjxpAlSxbCwsLYvn271uG8FmlRFBw+fJiePXumTUBCCCHSlBQWQghNODk5MWXKFO7du6d1KGkmJibG4u89f/48NWvWxNvbm5w5c6ZhVNp7lbw8L3fu3Li4uKTZ+wkhhEg7UlgIITTRsGFD8uXLR0BAwAv3GTt2bKJuQTNnzsTHx8e0/PHHH9O6dWsmTZpE3rx5yZYtG1999RWxsbEMGzaMHDly4OnpyaJFixK9/z///EP16tVxcnKidOnS7Ny502z7qVOnePfdd3F1dSVv3rx07tyZ27dvm7bXrVsXf39/hgwZQq5cuWjUqFGSv4der2fcuHF4enri6OhI+fLl2bx5s2m7TqcjJCSEcePGodPpGDt2bJLvs2rVKnx9fXF2diZnzpw0bNiQR48emWIZNGiQ2f6tW7fm448/Ni37+Pgwfvx4OnbsiKurKx4eHsyePdvse3Q6HfPmzaNZs2Y4OztTqFAhVq5cabbP8ePHqV+/vimOnj17EhUVZdpu/JsEBATg4eFBsWLFqFu3LpcvX2bw4MHodLpk70yNHTuWggUL4ujoiIeHBwMGDDD7HYx3PRYvXmx6r4RfCfMXFBREyZIlcXJyokSJEsydO9e0LSYmBn9/f/Lnz4+TkxM+Pj7Jfh6FEEIkTwoLIYQm7O3tmTRpErNnz+batWuv9F47duwgPDyc3bt3M336dMaOHUuLFi3Inj07Bw8epHfv3vTu3ZurV6+afd+wYcMYOnQoR48epXr16rRq1Yo7d+4AcP36derUqUP58uUJDg5m8+bN3Lhxg/bt25u9x5IlS3BwcGDv3r0sWLAgyfhmzZrFtGnT+Oabbzh27BhNmjShVatWnD171vSzSpcuzdChQ7l+/TqfffZZove4fv06H374Id27d+f06dPs3LmTNm3aoJRKVa6+/vprypYty5EjRxgxYgSDBw9m69atZvuMHj2aDz74gL///ptOnTrx4Ycfcvr0aQAeP35M06ZNyZ49O4cPH2blypVs27YNf39/s/fYvn07p0+fZuvWrWzcuJE1a9bg6enJuHHjuH79OtevX08yvlWrVjFjxgwWLFjA2bNnWbduHb6+vknu6+fnZ3qv69ev88svv+Dg4ECNGjUACAwM5Msvv2TixImcPn2aSZMmMXr0aJYsWQLAt99+y/r16/nf//5HWFgYP/30k1nRKoQQIpWUEEK8YV27dlXvvfeeUkqpd955R3Xv3l0ppdTatWtVwv+WxowZo8qVK2f2vTNmzFDe3t5m7+Xt7a3i4uJM64oXL65q1aplWo6NjVVZsmRRv/zyi1JKqYsXLypATZ482bTPs2fPlKenp5oyZYpSSqnRo0erxo0bm/3sq1evKkCFhYUppZSqU6eOKl++/Et/Xw8PDzVx4kSzdVWqVFF9+/Y1LZcrV06NGTPmhe8REhKiAHXp0qUkt9epU0cNHDjQbN17772nunbtalr29vZWTZs2NdvHz89PNWvWzLQMqN69e5vtU7VqVdWnTx+llFILFy5U2bNnV1FRUabtmzZtUnZ2dioiIkIpZfib5M2bV0VHR5u9j7e3t5oxY8YLf0ellJo2bZoqVqyYiomJSXL7i97j3LlzKmfOnGrq1KmmdV5eXmrZsmVm+40fP15Vq1ZNKaVU//79Vf369ZVer082JiGEECkjdyyEEJqaMmUKS5Ys4dSpUxa/R+nSpbGzi//vLG/evGZXue3t7cmZMyc3b940+75q1aqZXjs4OFC5cmXTlfmQkBD+/PNPXF1dTV8lSpQADOMhjCpXrpxsbA8ePCA8PNx0Fd2oRo0app+VEuXKlaNBgwb4+vrSrl07AgMDLRqfkvB3Ni4/H0dy+5w+fZpy5cqRJUsW0/YaNWqg1+sJCwszrfP19SVz5sypjq9du3Y8efKEwoUL8+mnn7J27VpiY2OT/Z7IyEhatGhBs2bNGDZsGAC3bt3i6tWr9OjRw+xvOGHCBNPf7+OPPyY0NJTixYszYMAAtmzZkup4hRBCxJPCQgihqdq1a9OkSRNGjhyZaJudnV2irj7Pnj1LtF+mTJnMlnU6XZLr9Hr9S+Mx9v3X6/W0bNmS0NBQs6+zZ89Su3Zt0/4JG9gpeV8jpVSqZsCyt7dn69at/P7775QqVYrZs2dTvHhxLl68CKQ8VymJLbl9kos74fqU5uV5Xl5ehIWFMWfOHJydnenbty+1a9d+4e8SFxeHn58fWbNmJTAw0LTe+LcODAw0+/udOHGCAwcOAFCxYkUuXrzI+PHjefLkCe3bt6dt27YWxS2EEEIKCyGEFZg8eTIbNmxg3759Zutz585NRESEWYM5LZ89YWxgAsTGxhISEmK6K1GxYkVOnjyJj48Pb731ltlXahrNWbNmxcPDgz179pit37dvHyVLlkxVvDqdjho1avDVV19x9OhRMmfOzNq1awFDrhKOW4iLi+PEiRPJ/s7GZePvnJJ9SpUqRWhoqGnQOMDevXuxs7OjWLFiycafOXNm4uLiXvp7Ojs706pVK7799lt27tzJ/v37OX78eJL7Dh48mOPHj7N27VqcnJxM6/PmzUuBAgW4cOFCor9foUKFTPtlzZoVPz8/AgMDWbFiBatXr+bu3bsvjVEIIURiDloHIIQQvr6+fPTRR4lmKKpbty63bt1i6tSptG3bls2bN/P777+TNWvWNPm5c+bMoWjRopQsWZIZM2Zw7949unfvDkC/fv0IDAzkww8/ZNiwYeTKlYtz586xfPlyAgMDsbe3T/HPGTZsGGPGjKFIkSKUL1+eoKAgQkND+fnnn1P8HgcPHmT79u00btyYPHnycPDgQW7dumUqTurXr8+QIUPYtGkTRYoUYcaMGdy/fz/R++zdu5epU6fSunVrtm7dysqVK9m0aZPZPitXrqRy5crUrFmTn3/+mUOHDvHDDz8A8NFHHzFmzBi6du3K2LFjuXXrFv3796dz587kzZs32d/Bx8eH3bt306FDBxwdHcmVK1eifRYvXkxcXBxVq1bFxcWFH3/8EWdnZ7y9vRPtGxQUxNy5c1m7di12dnZEREQAmLo9jR07lgEDBpA1a1aaNWtGdHQ0wcHB3Lt3jyFDhjBjxgzy589P+fLlsbOzY+XKleTLl49s2bKl5E8ihBDiOXLHQghhFcaPH5+oK0/JkiWZO3cuc+bMoVy5chw6dCjJGZMsNXnyZKZMmUK5cuX466+/+PXXX02NXQ8PD/bu3UtcXBxNmjShTJkyDBw4EHd3d7PxHCkxYMAAhg4dytChQ/H19WXz5s2sX7+eokWLpvg9smbNyu7du3n33XcpVqwYo0aNYtq0aTRr1gyA7t2707VrV7p06UKdOnUoVKgQ9erVS/Q+Q4cOJSQkhAoVKjB+/HimTZtGkyZNzPb56quvWL58OWXLlmXJkiX8/PPPlCpVCgAXFxf++OMP7t69S5UqVWjbti0NGjTgu+++e+nvMG7cOC5dukSRIkXInTt3kvtky5aNwMBAatSoQdmyZdm+fTsbNmxI8tkeu3btIi4ujlatWpE/f37T1zfffAPAJ598wvfff8/ixYvx9fWlTp06LF682HTHwtXVlSlTplC5cmWqVKnCpUuX+O2331L99xVCCGGgU8+fyYUQQqRLPj4+DBo0KNHzLhLS6XSsXbuW1q1bv7G4hBBCpA9yWUYIIYQQQgjxyqSwEEIIIYQQQrwy6QolhBBCCCGEeGVyx0IIIYQQQgjxyqSwEEIIIYQQQrwyKSyEEEIIIYQQr0wKCyGEEEIIIcQrk8JCCCGEEEII8cqksBBCCCGEEEK8MikshBBCCCGEEK9MCgshhBBCCCHEK5PCQgghhBBCCPHK/h/dXzrG576mjQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 800x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,4))\n",
    "ax.plot(complexities, relative_decreasing_aics, color='blue', linewidth=2.5, label=\"AIC-decreasing models\")\n",
    "ax.scatter(complexities, relative_decreasing_aics, marker='*', c='red', s=300)\n",
    "for com, name, rel_aic in zip(complexities, eff_feature_names, relative_decreasing_aics):\n",
    "    rotation = (90+45)/2\n",
    "    if com == 1: rotation = 0\n",
    "    ax.annotate(f\"    $({','.join(name)})$\", (com, rel_aic), rotation=rotation)\n",
    "plt.xlim(0.8,8.1)\n",
    "plt.title(r\"Model selection based on $\\frac{\\Delta AIC}{\\Delta s}$\")\n",
    "plt.ylabel(\"Relative AIC\")\n",
    "plt.xlabel(\"Number of support sizes\")\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "2feae1e6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 30275.8339728   25599.5127601  -12858.14249527 -12864.25869979\n",
      " -12865.05464797 -12880.66375586 -12898.78130837 -12904.31448525]\n",
      "1 ['uu_{x}']\n",
      "2 ['u_{xx}', 'uu_{x}']\n",
      "3 ['u_{xx}', 'u_{xxxx}', 'uu_{x}']\n",
      "4 ['u_{xx}', 'u_{xxx}', 'u_{xxxx}', 'uu_{x}']\n",
      "5 ['u', 'u_{xx}', 'u_{xxxx}', 'uu_{x}', 'u^2u_{xx}']\n",
      "6 ['u_{x}', 'u_{xx}', 'u_{xxxx}', 'uu_{x}', 'u^2u_{xx}', 'u^2u_{xxxx}']\n",
      "7 ['u_{x}', 'u_{xx}', 'u_{xxxx}', 'uu_{x}', 'u^2u_{x}', 'u^2u_{xx}', 'u^2u_{xxxx}']\n",
      "8 ['u', 'u_{x}', 'u_{xx}', 'u_{xxxx}', 'uu_{x}', 'u^2u_{x}', 'u^2u_{xx}', 'u^2u_{xxxx}']\n",
      "3 improves 2\n",
      "Percent improve: 1.502280751035022\n",
      "7 does not improve 3\n",
      "Percent improve: 0.0007901377107136002\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(decreasing_aics)\n",
    "for c in complexities: print(c, eff_feature_names[c-1])\n",
    "find_transition_V5(decreasing_aics, complexities, percent_ok=0.09)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "2941b177",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.46999004572032915, 0.21721774824355908)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "errs = percent_coeff(model_results[3][-1])\n",
    "errs.mean(), errs.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51f8c73b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f13cade",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "fa4b0771",
   "metadata": {},
   "source": [
    "## Calculate a weak form of each canidate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "d42d44e6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.01072770358302814\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[-0.99983724],\n",
       "       [-0.99991587],\n",
       "       [-0.99992506]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "acc=4\n",
    "# Define the derivative:\n",
    "d_dx = FinDiff(0, dx, 1, acc=acc)\n",
    "d_dxx = FinDiff(0, dx, 2, acc=acc)\n",
    "d_dxxxx = FinDiff(0, dx, 4, acc=acc)\n",
    "d_dt = FinDiff(1, dt, 1, acc=acc)\n",
    "\n",
    "coeff = np.linalg.lstsq(np.stack([u*d_dx(u), d_dxx(u), d_dxxxx(u)]).reshape(3, -1).T, \n",
    "                        d_dt(u).reshape(-1,1), rcond=None)[0]\n",
    "print((100*np.abs(coeff+1)).mean())\n",
    "coeff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "d5666ce5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.008675862373541676\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[-1.00009514],\n",
       "       [-1.00010748],\n",
       "       [-1.00005765]])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Xf = np.stack([integrate.trapz(d_dxx(u), t, dt, axis=1), \n",
    "               integrate.trapz(d_dxxxx(u), t, dt, axis=1), \n",
    "               integrate.trapz(u*d_dx(u), t, dt, axis=1)]).T\n",
    "yf = integrate.trapz(d_dt(u), t, dt, axis=1).reshape(-1,1)\n",
    "coeff = np.linalg.lstsq(Xf, yf, rcond=None)[0]\n",
    "print((100*np.abs(coeff+1)).mean()) # 0.00867587799540844\n",
    "coeff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "19bbfa1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# sio.savemat(\"../../PDE_Discovery_Weak_Formulation/Datasets/KS_data_tmp.mat\", \n",
    "#             {'uu': un, 'xx':x, 'tt':t, 'dx':dx, 'dt': dt})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "d5c6f331",
   "metadata": {},
   "outputs": [],
   "source": [
    "def weight_poly(x, m, k):\n",
    "    \"\"\"\n",
    "    Polynomial piece of weighting function used to satisfy BC A = d^k/dx^k[(x^2 - 1)^m]\n",
    "    x: independent variable\n",
    "    m: power of base function\n",
    "    k: order of derivative\n",
    "    \"\"\"\n",
    "    a = np.zeros((2*m+1, 1)) # initial coefficient vector\n",
    "    for l in range(m+1):\n",
    "        a[2*l] = ((-1)**(m-l))*np.math.comb(m, l) # set polynomial coefficients\n",
    "    \n",
    "    c = np.zeros((2*m+1, 1)) # final coefficient vector\n",
    "    for n in range(2*m-k+1):\n",
    "        c[n] = a[n+k]*np.math.factorial(n+k)/np.math.factorial(n)\n",
    "    \n",
    "    p = 0\n",
    "    for n in range(2*m-k+1):\n",
    "        p += c[n]*(x**n) # final windowing function\n",
    "\n",
    "    return p\n",
    "\n",
    "def weight_full(k, p, x, t):\n",
    "    \"\"\"\n",
    "    Assemble the 1D weight functions into the full weight\n",
    "    k = [kx,ky,kt]: order of derivative(s)\n",
    "    p = [px,py,pt]: exponents of weight polynomials\n",
    "    \"\"\"\n",
    "    if len(k) == 3:\n",
    "        wx = weight_poly(x, p[0], k[0])\n",
    "        wy = weight_poly(x, p[1], k[1])\n",
    "        wt = weight_poly(t, p[2], k[2])\n",
    "        wX, wY, wT = np.meshgrid(wx, wy, wt)\n",
    "        W = wX * wY * wT\n",
    "    elif len(k) == 2:\n",
    "        wx = weight_poly(x, p[0], k[0])\n",
    "        wt = weight_poly(t, p[1], k[1])\n",
    "        wT, wX = np.meshgrid(wt, wx)\n",
    "        W = wX * wT\n",
    "    \n",
    "    return W, wx, wt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "ecec1660",
   "metadata": {},
   "outputs": [],
   "source": [
    "u_weak = un.copy()\n",
    "n_domain = 10000\n",
    "x_size, t_size = 100, 100\n",
    "\n",
    "xsup = np.linspace(-1,1,x_size+1)\n",
    "tsup = np.linspace(-1,1,t_size+1)\n",
    "S_x = 2/(dx*x_size)\n",
    "S_t = 2/(dt*t_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "280cc095",
   "metadata": {},
   "source": [
    "#### Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "3f253123",
   "metadata": {},
   "outputs": [],
   "source": [
    "# p = [4,3]\n",
    "# dA01, _, _ = weight_full([0,1], p, xsup, tsup)\n",
    "# dA10, _, _ = weight_full([1,0], p, xsup, tsup)\n",
    "# dA20, _, _ = weight_full([2,0], p, xsup, tsup)\n",
    "# dA40, _, _ = weight_full([4,0], p, xsup, tsup)\n",
    "\n",
    "# X_weak = []\n",
    "# y_weak = []\n",
    "# for _ in range(n_domain):\n",
    "#     ix = np.random.randint(0, len(x)-x_size, size=1)[0] + np.arange(0, x_size+1)\n",
    "#     it = np.random.randint(0, len(t)-t_size, size=1)[0] + np.arange(0, t_size+1)\n",
    "    \n",
    "#     usup = u_weak[ix, :][:, it]\n",
    "    \n",
    "#     f1 = -(1/2)*(usup**2)*dA10*S_x\n",
    "#     f1 = trapz(trapz(f1, tsup, 1), xsup, 0)\n",
    "    \n",
    "#     f2 = usup*dA20*(S_x**2)\n",
    "#     f2 = trapz(trapz(f2, tsup, 1), xsup, 0)\n",
    "    \n",
    "#     f3 = usup*dA40*(S_x**4)\n",
    "#     f3 = trapz(trapz(f3, tsup, 1), xsup, 0)\n",
    "    \n",
    "#     X_weak.append([f1, f2, f3])\n",
    "    \n",
    "#     f4 = -usup*dA01*S_t\n",
    "#     f4 = trapz(trapz(f4, tsup, 1), xsup, 0)\n",
    "#     y_weak.append(f4)\n",
    "    \n",
    "# X_weak = np.array(X_weak)\n",
    "# y_weak = np.array(y_weak)\n",
    "\n",
    "# np.linalg.lstsq(X_weak, y_weak, rcond=None)[0]\n",
    "# coeff = np.linalg.lstsq(X_weak, y_weak, rcond=None)[0]\n",
    "# coeff, (np.abs(coeff+1)*100).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "bf842b39",
   "metadata": {},
   "outputs": [],
   "source": [
    "p, q = 4, 4 # p, q = 5, 4 also good\n",
    "wx, wt = np.meshgrid((xsup**2-1)**p, (tsup**2-1)**q)\n",
    "W = wx*wt\n",
    "\n",
    "# Define the derivative:\n",
    "d_dx = FinDiff(0, xsup[1]-xsup[0], 1)\n",
    "d_dxx = FinDiff(0, xsup[1]-xsup[0], 2)\n",
    "d_dxxxx = FinDiff(0, xsup[1]-xsup[0], 4)\n",
    "d_dt = FinDiff(1, tsup[1]-tsup[0], 1)\n",
    "\n",
    "dA01 = d_dt(W)\n",
    "dA10 = d_dx(W)\n",
    "dA20 = d_dxx(W)\n",
    "dA40 = d_dxxxx(W)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "06636cbf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "savgol\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 15/15 [06:53<00:00, 27.55s/it]\n"
     ]
    }
   ],
   "source": [
    "# tuning denoising parameters based on min norm(X_weak@coeff-y_weak) !?\n",
    "losses = []; heapq.heapify(losses)\n",
    "denoiser = 'savgol'\n",
    "step = 2 if denoiser == 'savgol' else 1\n",
    "print(denoiser)\n",
    "for ws in trange(3, 33, step):\n",
    "    X_weak = []\n",
    "    y_weak = []\n",
    "    for _ in range(n_domain):\n",
    "        ix = np.random.randint(0, len(x)-x_size, size=1)[0] + np.arange(0, x_size+1)\n",
    "        it = np.random.randint(0, len(t)-t_size, size=1)[0] + np.arange(0, t_size+1)\n",
    "        \n",
    "        usup = savgol_filter(u_weak[ix, :][:, it], ws, polyorder=2)\n",
    "        usup = wiener(u_weak[ix, :][:, it], ws, noise=None)\n",
    "\n",
    "        f1 = -(1/2)*(usup**2)*dA10*S_x\n",
    "        f1 = trapz(trapz(f1, tsup, 1), xsup, 0)\n",
    "\n",
    "        f2 = usup*dA20*(S_x**2)\n",
    "        f2 = trapz(trapz(f2, tsup, 1), xsup, 0)\n",
    "\n",
    "        f3 = usup*dA40*(S_x**4)\n",
    "        f3 = trapz(trapz(f3, tsup, 1), xsup, 0)\n",
    "\n",
    "        X_weak.append([f1, f2, f3])\n",
    "\n",
    "        f4 = -usup*dA01*S_t\n",
    "        f4 = trapz(trapz(f4, tsup, 1), xsup, 0)\n",
    "        y_weak.append(f4)\n",
    "\n",
    "    X_weak = np.array(X_weak)\n",
    "    y_weak = np.array(y_weak)\n",
    "\n",
    "    coeff = np.linalg.lstsq(X_weak, y_weak, rcond=None)[0]\n",
    "    loss = ((X_weak@coeff-y_weak)**2).mean()\n",
    "    \n",
    "    losses.append((loss, X_weak, y_weak, coeff))\n",
    "    losses = heapq.nsmallest(10, losses)\n",
    "\n",
    "losses = sorted(losses)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "4267e12c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.021832395563001\n",
      "1.8977996022640904\n",
      "2.5133935538041627\n",
      "3.0511071077478005\n",
      "3.440400704507855\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([-1.02511986, -1.02981218, -1.04827998])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# print(losses[:n_tops])\n",
    "n_tops = 5 # same significant number A: A.xxxe-y\n",
    "avg_denoised_weak_coeff = 0\n",
    "for i in range(n_tops):\n",
    "    avg_denoised_weak_coeff += losses[i][-1]\n",
    "    print(percent_coeff(avg_denoised_weak_coeff/(i+1)).mean())\n",
    "avg_denoised_weak_coeff = avg_denoised_weak_coeff/(n_tops)\n",
    "avg_denoised_weak_coeff\n",
    "\n",
    "# X_weak, y_weak = 0., 0.\n",
    "# for i in range(5):\n",
    "#     X_weak += losses[i][1]\n",
    "#     y_weak += losses[i][2]\n",
    "# X_weak /= n_tops\n",
    "# y_weak /= n_tops\n",
    "# percent_coeff(np.linalg.lstsq(X_weak, y_weak, rcond=None)[0]).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "1c8ae6e6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(3.9218733651876165e-06,\n",
       "  array([[ 0.08257882, -0.03687394, -0.04376773],\n",
       "         [-0.01337275, -0.0072121 ,  0.00680215],\n",
       "         [ 0.00109509,  0.00905165, -0.00603873],\n",
       "         ...,\n",
       "         [-0.00825402,  0.14312572, -0.07744333],\n",
       "         [ 0.05028671, -0.07998783,  0.02141321],\n",
       "         [ 0.00127049,  0.01686838, -0.01173849]]),\n",
       "  array([-1.93476859e-05,  1.58406514e-02, -3.14225527e-03, ...,\n",
       "         -5.67749348e-02,  9.56400248e-03, -6.29575293e-03]),\n",
       "  array([-1.00578273, -1.01026803, -1.01460421])),\n",
       " (6.913193480969802e-06,\n",
       "  array([[ 0.01695572, -0.13478587,  0.10936995],\n",
       "         [-0.00868196, -0.00893185,  0.01384142],\n",
       "         [ 0.02706135,  0.0393157 , -0.0141495 ],\n",
       "         ...,\n",
       "         [-0.08337341, -0.02903938,  0.10867105],\n",
       "         [ 0.11850393, -0.18407172,  0.05798834],\n",
       "         [-0.07793961,  0.26757709, -0.18453369]]),\n",
       "  array([ 0.00827385,  0.00599532, -0.05524916, ...,  0.00071997,\n",
       "          0.00773445, -0.0057381 ]),\n",
       "  array([-1.02133969, -1.02641914, -1.03545418])),\n",
       " (1.3006029566075381e-05,\n",
       "  array([[ 0.07773   , -0.26363009,  0.16161822],\n",
       "         [ 0.00735836,  0.05017395, -0.06156669],\n",
       "         [-0.07059055,  0.2396383 , -0.16452653],\n",
       "         ...,\n",
       "         [-0.02484593,  0.02485609, -0.04266426],\n",
       "         [ 0.07053598, -0.21911286,  0.15725725],\n",
       "         [ 0.15662912, -0.29966616,  0.14238113]]),\n",
       "  array([ 2.12366851e-02,  5.18028423e-03,  5.03492332e-05, ...,\n",
       "          3.90759910e-02, -1.92747361e-02,  2.10224582e-04]),\n",
       "  array([-1.0287893 , -1.03349409, -1.05005406])),\n",
       " (2.1799545339912492e-05,\n",
       "  array([[ 0.0096619 , -0.00184803, -0.00393732],\n",
       "         [ 0.03294723, -0.14002594,  0.09884441],\n",
       "         [-0.12480687,  0.21583851, -0.10153252],\n",
       "         ...,\n",
       "         [ 0.00585412, -0.17221191,  0.12154552],\n",
       "         [-0.09145148,  0.15135691, -0.06222279],\n",
       "         [-0.03965582,  0.04293985, -0.00494842]]),\n",
       "  array([-0.01113958,  0.00852876,  0.01081167, ...,  0.04291615,\n",
       "          0.00544632,  0.00210292]),\n",
       "  array([-1.0352138 , -1.03942867, -1.06528496])),\n",
       " (3.5664366079346574e-05,\n",
       "  array([[ 0.00389666, -0.03548249,  0.01511611],\n",
       "         [-0.04553414,  0.09208309, -0.03883753],\n",
       "         [ 0.02088138, -0.01236522,  0.0028568 ],\n",
       "         ...,\n",
       "         [-0.02061217,  0.07222077, -0.04180216],\n",
       "         [ 0.00088965,  0.00773475, -0.00410008],\n",
       "         [-0.04836059,  0.04219193,  0.0118996 ]]),\n",
       "  array([ 0.02668554, -0.0059649 , -0.00312116, ..., -0.00798618,\n",
       "         -0.00564119, -0.00820369]),\n",
       "  array([-1.0344738 , -1.03945096, -1.0760025 ])),\n",
       " (5.232175770661518e-05,\n",
       "  array([[ 0.00649777, -0.0009705 , -0.00821371],\n",
       "         [ 0.01687908, -0.07773553,  0.02000276],\n",
       "         [-0.12414218,  0.18235998, -0.06482215],\n",
       "         ...,\n",
       "         [-0.0451568 ,  0.08479695, -0.05043195],\n",
       "         [ 0.04826807, -0.07145146,  0.00396652],\n",
       "         [ 0.01257037, -0.09968219,  0.0726386 ]]),\n",
       "  array([-0.00697565,  0.03209204,  0.01411594, ...,  0.0107783 ,\n",
       "          0.02328503,  0.01088965]),\n",
       "  array([-1.01765957, -1.01899329, -1.06302934])),\n",
       " (7.144055175563794e-05,\n",
       "  array([[-0.06386146,  0.05371681, -0.02447554],\n",
       "         [ 0.1280545 , -0.09956064, -0.01823773],\n",
       "         [ 0.08352246, -0.14586773,  0.05316688],\n",
       "         ...,\n",
       "         [-0.02588993, -0.112117  ,  0.1214644 ],\n",
       "         [-0.0314522 ,  0.07380332, -0.04514729],\n",
       "         [-0.0227159 , -0.04633548,  0.05527025]]),\n",
       "  array([ 0.04299293,  0.00075224, -0.01608345, ...,  0.00644818,\n",
       "          0.00016175,  0.00515888]),\n",
       "  array([-0.98999804, -0.99407257, -1.04994903])),\n",
       " (9.127798845302397e-05,\n",
       "  array([[-0.04812328, -0.07441925,  0.10458453],\n",
       "         [-0.009995  ,  0.03886031, -0.01052762],\n",
       "         [ 0.14072811, -0.12030428, -0.00067809],\n",
       "         ...,\n",
       "         [-0.03157333,  0.0293129 , -0.02059509],\n",
       "         [ 0.00124449, -0.01276824,  0.01144961],\n",
       "         [ 0.04287749,  0.00666978, -0.04330066]]),\n",
       "  array([ 0.00152203, -0.00957844, -0.01180018, ...,  0.01782602,\n",
       "          0.00634506,  0.00302596]),\n",
       "  array([-0.95694539, -0.96126841, -1.02632432])),\n",
       " (0.00011318182784530615,\n",
       "  array([[-6.03978849e-03, -6.40922320e-02,  4.30425514e-02],\n",
       "         [-3.77226963e-03,  5.48299994e-03, -1.39737575e-04],\n",
       "         [-5.44915261e-02, -8.15003787e-02,  9.21264078e-02],\n",
       "         ...,\n",
       "         [ 3.76035065e-02, -1.41036760e-01,  6.06529455e-02],\n",
       "         [ 2.40577340e-02, -5.18689467e-02,  1.50893347e-02],\n",
       "         [ 8.95324006e-02, -2.53283760e-01,  1.68953629e-01]]),\n",
       "  array([ 0.04192653,  0.00314559,  0.03891514, ...,  0.03142542,\n",
       "          0.01528115, -0.02703658]),\n",
       "  array([-0.90255228, -0.90455609, -0.97251326])),\n",
       " (0.00013337660661913748,\n",
       "  array([[ 0.00987576,  0.05491446, -0.05325866],\n",
       "         [-0.06020636,  0.29115839, -0.21087157],\n",
       "         [ 0.00376262, -0.02406232,  0.00986413],\n",
       "         ...,\n",
       "         [-0.00167452,  0.01607946, -0.01217901],\n",
       "         [-0.0025002 ,  0.01162738, -0.00358885],\n",
       "         [ 0.0667074 , -0.29440329,  0.20128965]]),\n",
       "  array([ 0.01199576, -0.00474138,  0.00934204, ..., -0.0057292 ,\n",
       "         -0.00321924,  0.01041669]),\n",
       "  array([-0.85315631, -0.85545265, -0.92848689]))]"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "losses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c00c7c19",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:pysr]",
   "language": "python",
   "name": "conda-env-pysr-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}