gempy3_first_model_generative copy 2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import all the libraries and packages "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/jt925938/anaconda3/envs/gempy_pyro/lib/python3.10/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n",
      "  from .autonotebook import tqdm as notebook_tqdm\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Setting Backend To: AvailableBackends.numpy\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import joblib\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import torch\n",
    "import torch.nn.functional as F\n",
    "\n",
    "import pyro\n",
    "import pyro.distributions as dist\n",
    "from pyro.infer import MCMC, NUTS, Predictive, EmpiricalMarginal\n",
    "from pyro.infer.autoguide import init_to_mean, init_to_median, init_to_value\n",
    "from pyro.infer.inspect import get_dependencies\n",
    "from pyro.infer import SVI, TraceEnum_ELBO, config_enumerate, infer_discrete\n",
    "\n",
    "import gempy as gp\n",
    "import gempy_engine\n",
    "import gempy_viewer as gpv\n",
    "from gempy_engine.core.backend_tensor import BackendTensor\n",
    "import arviz as az\n",
    "from gempy_probability.plot_posterior import default_red, default_blue, PlotPosterior\n",
    "\n",
    "from sklearn.mixture import GaussianMixture\n",
    "from sklearn.mixture import BayesianGaussianMixture\n",
    "from sklearn.cluster import KMeans\n",
    "from scipy.stats import multivariate_normal, norm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get the path of all the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename_a = './Fw__Hyperspectral_datasets_from_the_KSL_cores/CuSp131.pkl'\n",
    "filename_b = './Fw__Hyperspectral_datasets_from_the_KSL_cores/CuSp133.pkl'\n",
    "filename_c = './Fw__Hyperspectral_datasets_from_the_KSL_cores/CuSp136.pkl'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load and get the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['XYZ', 'MWL_SWIRKaol', 'MWL_SWIRClay', 'BR_Fe3_Fe2_FeOH', 'BR_Anhydrite', 'BR_Gypsum', 'BR_Qtz_Fsp_Cal', 'MWL_MWIRCarbonate', 'BR_Kaol_AlOH_H2O', 'MWL_Carbonate', 'BR_Fe2_FeOH'])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with open(filename_a, 'rb') as myfile:\n",
    "    a= joblib.load(myfile)\n",
    "a.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Description about Hyperspectral data\n",
    "## Data is obtained by scannig the core data from different boreholes using hyperspectral sensors. There were around 450 channels for each pixels initially. It was preprocessed and seperated based on 10 different types of rocks. In each of the file we have \"X\", \"Y\", \"Z\" coordinates points corresponding to sensors and corresponding to each rock type we have a transformed RGB correspondence information. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(235, 34)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "column_name =[]\n",
    "for keys, _ in a.items():\n",
    "    if keys=='XYZ':\n",
    "        column_name.append(\"Borehole_id\")\n",
    "        column_name.append(\"X\")\n",
    "        column_name.append(\"Y\")\n",
    "        column_name.append(\"Z\")\n",
    "    else:\n",
    "        column_name.append(keys+\"_R\")\n",
    "        column_name.append(keys+\"_G\")\n",
    "        column_name.append(keys+\"_B\")\n",
    "data_a =[]\n",
    "for keys, values in a.items():\n",
    "    if keys=='XYZ':\n",
    "        label_a = np.ones((235,1))\n",
    "        data_a.append(label_a)\n",
    "        data_a.append(values)\n",
    "    else:\n",
    "        data_a.append(values/255)\n",
    "\n",
    "# Concatenate the arrays horizontally to create an array of size 5x30\n",
    "concatenated_array_a = np.hstack(data_a)\n",
    "# sort the data based on the depth\n",
    "sorted_indices = np.argsort(-concatenated_array_a[:, 3])\n",
    "concatenated_array_a = concatenated_array_a[sorted_indices]\n",
    "concatenated_array_a.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "dataframe_KSL = pd.DataFrame(concatenated_array_a,columns=column_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Borehole_id</th>\n",
       "      <th>X</th>\n",
       "      <th>Y</th>\n",
       "      <th>Z</th>\n",
       "      <th>MWL_SWIRKaol_R</th>\n",
       "      <th>MWL_SWIRKaol_G</th>\n",
       "      <th>MWL_SWIRKaol_B</th>\n",
       "      <th>MWL_SWIRClay_R</th>\n",
       "      <th>MWL_SWIRClay_G</th>\n",
       "      <th>MWL_SWIRClay_B</th>\n",
       "      <th>...</th>\n",
       "      <th>MWL_MWIRCarbonate_B</th>\n",
       "      <th>BR_Kaol_AlOH_H2O_R</th>\n",
       "      <th>BR_Kaol_AlOH_H2O_G</th>\n",
       "      <th>BR_Kaol_AlOH_H2O_B</th>\n",
       "      <th>MWL_Carbonate_R</th>\n",
       "      <th>MWL_Carbonate_G</th>\n",
       "      <th>MWL_Carbonate_B</th>\n",
       "      <th>BR_Fe2_FeOH_R</th>\n",
       "      <th>BR_Fe2_FeOH_G</th>\n",
       "      <th>BR_Fe2_FeOH_B</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.0</td>\n",
       "      <td>101726.067469</td>\n",
       "      <td>1.256648e+06</td>\n",
       "      <td>-192.063946</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.007843</td>\n",
       "      <td>...</td>\n",
       "      <td>0.007843</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.007843</td>\n",
       "      <td>0.015686</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.003922</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0</td>\n",
       "      <td>462396.723231</td>\n",
       "      <td>5.712055e+06</td>\n",
       "      <td>-757.994770</td>\n",
       "      <td>0.027451</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.019608</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>...</td>\n",
       "      <td>0.203922</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.956863</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>0.011765</td>\n",
       "      <td>0.007843</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0</td>\n",
       "      <td>462396.707804</td>\n",
       "      <td>5.712055e+06</td>\n",
       "      <td>-758.491319</td>\n",
       "      <td>0.019608</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.007843</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>...</td>\n",
       "      <td>0.180392</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.917647</td>\n",
       "      <td>0.035294</td>\n",
       "      <td>0.019608</td>\n",
       "      <td>0.011765</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.0</td>\n",
       "      <td>462396.691591</td>\n",
       "      <td>5.712055e+06</td>\n",
       "      <td>-758.987619</td>\n",
       "      <td>0.027451</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.019608</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>...</td>\n",
       "      <td>0.262745</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.066667</td>\n",
       "      <td>0.023529</td>\n",
       "      <td>0.015686</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.0</td>\n",
       "      <td>462396.675420</td>\n",
       "      <td>5.712055e+06</td>\n",
       "      <td>-759.483931</td>\n",
       "      <td>0.015686</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>...</td>\n",
       "      <td>0.188235</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.686275</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>0.039216</td>\n",
       "      <td>0.019608</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>230</th>\n",
       "      <td>1.0</td>\n",
       "      <td>462391.365850</td>\n",
       "      <td>5.712036e+06</td>\n",
       "      <td>-870.771163</td>\n",
       "      <td>0.007843</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.015686</td>\n",
       "      <td>0.062745</td>\n",
       "      <td>0.082353</td>\n",
       "      <td>...</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.082353</td>\n",
       "      <td>0.156863</td>\n",
       "      <td>0.011765</td>\n",
       "      <td>0.007843</td>\n",
       "      <td>0.011765</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.019608</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>231</th>\n",
       "      <td>1.0</td>\n",
       "      <td>462391.333756</td>\n",
       "      <td>5.712036e+06</td>\n",
       "      <td>-871.258504</td>\n",
       "      <td>0.015686</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.035294</td>\n",
       "      <td>0.121569</td>\n",
       "      <td>0.180392</td>\n",
       "      <td>...</td>\n",
       "      <td>0.133333</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.172549</td>\n",
       "      <td>0.321569</td>\n",
       "      <td>0.027451</td>\n",
       "      <td>0.019608</td>\n",
       "      <td>0.023529</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.058824</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>232</th>\n",
       "      <td>1.0</td>\n",
       "      <td>462391.301701</td>\n",
       "      <td>5.712036e+06</td>\n",
       "      <td>-871.745859</td>\n",
       "      <td>0.019608</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.027451</td>\n",
       "      <td>0.054902</td>\n",
       "      <td>0.133333</td>\n",
       "      <td>...</td>\n",
       "      <td>0.113725</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.160784</td>\n",
       "      <td>0.286275</td>\n",
       "      <td>0.027451</td>\n",
       "      <td>0.023529</td>\n",
       "      <td>0.027451</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.062745</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>233</th>\n",
       "      <td>1.0</td>\n",
       "      <td>462391.268635</td>\n",
       "      <td>5.712036e+06</td>\n",
       "      <td>-872.233237</td>\n",
       "      <td>0.015686</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.027451</td>\n",
       "      <td>0.039216</td>\n",
       "      <td>0.117647</td>\n",
       "      <td>...</td>\n",
       "      <td>0.090196</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.149020</td>\n",
       "      <td>0.278431</td>\n",
       "      <td>0.027451</td>\n",
       "      <td>0.023529</td>\n",
       "      <td>0.027451</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.058824</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>234</th>\n",
       "      <td>1.0</td>\n",
       "      <td>462391.235636</td>\n",
       "      <td>5.712036e+06</td>\n",
       "      <td>-872.720620</td>\n",
       "      <td>0.015686</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.023529</td>\n",
       "      <td>0.035294</td>\n",
       "      <td>0.113725</td>\n",
       "      <td>...</td>\n",
       "      <td>0.101961</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.125490</td>\n",
       "      <td>0.207843</td>\n",
       "      <td>0.023529</td>\n",
       "      <td>0.019608</td>\n",
       "      <td>0.023529</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.027451</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>235 rows × 34 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     Borehole_id              X             Y           Z  MWL_SWIRKaol_R  \\\n",
       "0            1.0  101726.067469  1.256648e+06 -192.063946        0.000000   \n",
       "1            1.0  462396.723231  5.712055e+06 -757.994770        0.027451   \n",
       "2            1.0  462396.707804  5.712055e+06 -758.491319        0.019608   \n",
       "3            1.0  462396.691591  5.712055e+06 -758.987619        0.027451   \n",
       "4            1.0  462396.675420  5.712055e+06 -759.483931        0.015686   \n",
       "..           ...            ...           ...         ...             ...   \n",
       "230          1.0  462391.365850  5.712036e+06 -870.771163        0.007843   \n",
       "231          1.0  462391.333756  5.712036e+06 -871.258504        0.015686   \n",
       "232          1.0  462391.301701  5.712036e+06 -871.745859        0.019608   \n",
       "233          1.0  462391.268635  5.712036e+06 -872.233237        0.015686   \n",
       "234          1.0  462391.235636  5.712036e+06 -872.720620        0.015686   \n",
       "\n",
       "     MWL_SWIRKaol_G  MWL_SWIRKaol_B  MWL_SWIRClay_R  MWL_SWIRClay_G  \\\n",
       "0          0.000000        0.000000        0.000000        0.000000   \n",
       "1          0.003922        0.003922        0.019608        0.003922   \n",
       "2          0.003922        0.003922        0.007843        0.000000   \n",
       "3          0.003922        0.003922        0.019608        0.003922   \n",
       "4          0.000000        0.000000        0.003922        0.000000   \n",
       "..              ...             ...             ...             ...   \n",
       "230        0.000000        0.000000        0.015686        0.062745   \n",
       "231        0.000000        0.000000        0.035294        0.121569   \n",
       "232        0.003922        0.003922        0.027451        0.054902   \n",
       "233        0.003922        0.003922        0.027451        0.039216   \n",
       "234        0.003922        0.003922        0.023529        0.035294   \n",
       "\n",
       "     MWL_SWIRClay_B  ...  MWL_MWIRCarbonate_B  BR_Kaol_AlOH_H2O_R  \\\n",
       "0          0.007843  ...             0.007843                 0.0   \n",
       "1          0.003922  ...             0.203922                 0.0   \n",
       "2          0.000000  ...             0.180392                 0.0   \n",
       "3          0.003922  ...             0.262745                 0.0   \n",
       "4          0.000000  ...             0.188235                 0.0   \n",
       "..              ...  ...                  ...                 ...   \n",
       "230        0.082353  ...             0.058824                 0.0   \n",
       "231        0.180392  ...             0.133333                 0.0   \n",
       "232        0.133333  ...             0.113725                 0.0   \n",
       "233        0.117647  ...             0.090196                 0.0   \n",
       "234        0.113725  ...             0.101961                 0.0   \n",
       "\n",
       "     BR_Kaol_AlOH_H2O_G  BR_Kaol_AlOH_H2O_B  MWL_Carbonate_R  MWL_Carbonate_G  \\\n",
       "0              0.007843            0.015686         0.000000         0.000000   \n",
       "1              0.000000            0.956863         0.058824         0.011765   \n",
       "2              0.000000            0.917647         0.035294         0.019608   \n",
       "3              0.000000            1.000000         0.066667         0.023529   \n",
       "4              0.000000            0.686275         0.058824         0.039216   \n",
       "..                  ...                 ...              ...              ...   \n",
       "230            0.082353            0.156863         0.011765         0.007843   \n",
       "231            0.172549            0.321569         0.027451         0.019608   \n",
       "232            0.160784            0.286275         0.027451         0.023529   \n",
       "233            0.149020            0.278431         0.027451         0.023529   \n",
       "234            0.125490            0.207843         0.023529         0.019608   \n",
       "\n",
       "     MWL_Carbonate_B  BR_Fe2_FeOH_R  BR_Fe2_FeOH_G  BR_Fe2_FeOH_B  \n",
       "0           0.000000            0.0            0.0       0.003922  \n",
       "1           0.007843            0.0            0.0       0.000000  \n",
       "2           0.011765            0.0            0.0       0.000000  \n",
       "3           0.015686            0.0            0.0       0.000000  \n",
       "4           0.019608            0.0            0.0       0.000000  \n",
       "..               ...            ...            ...            ...  \n",
       "230         0.011765            0.0            0.0       0.019608  \n",
       "231         0.023529            0.0            0.0       0.058824  \n",
       "232         0.027451            0.0            0.0       0.062745  \n",
       "233         0.027451            0.0            0.0       0.058824  \n",
       "234         0.023529            0.0            0.0       0.027451  \n",
       "\n",
       "[235 rows x 34 columns]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dataframe_KSL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>MWL_SWIRKaol_R</th>\n",
       "      <th>MWL_SWIRKaol_G</th>\n",
       "      <th>MWL_SWIRKaol_B</th>\n",
       "      <th>MWL_SWIRClay_R</th>\n",
       "      <th>MWL_SWIRClay_G</th>\n",
       "      <th>MWL_SWIRClay_B</th>\n",
       "      <th>BR_Fe3_Fe2_FeOH_R</th>\n",
       "      <th>BR_Fe3_Fe2_FeOH_G</th>\n",
       "      <th>BR_Fe3_Fe2_FeOH_B</th>\n",
       "      <th>BR_Anhydrite_R</th>\n",
       "      <th>...</th>\n",
       "      <th>MWL_MWIRCarbonate_B</th>\n",
       "      <th>BR_Kaol_AlOH_H2O_R</th>\n",
       "      <th>BR_Kaol_AlOH_H2O_G</th>\n",
       "      <th>BR_Kaol_AlOH_H2O_B</th>\n",
       "      <th>MWL_Carbonate_R</th>\n",
       "      <th>MWL_Carbonate_G</th>\n",
       "      <th>MWL_Carbonate_B</th>\n",
       "      <th>BR_Fe2_FeOH_R</th>\n",
       "      <th>BR_Fe2_FeOH_G</th>\n",
       "      <th>BR_Fe2_FeOH_B</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.007843</td>\n",
       "      <td>0.023529</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>...</td>\n",
       "      <td>0.007843</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.007843</td>\n",
       "      <td>0.015686</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.003922</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.027451</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.019608</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.949020</td>\n",
       "      <td>...</td>\n",
       "      <td>0.203922</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.956863</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>0.011765</td>\n",
       "      <td>0.007843</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.019608</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.007843</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.909804</td>\n",
       "      <td>...</td>\n",
       "      <td>0.180392</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.917647</td>\n",
       "      <td>0.035294</td>\n",
       "      <td>0.019608</td>\n",
       "      <td>0.011765</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.027451</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.019608</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.992157</td>\n",
       "      <td>...</td>\n",
       "      <td>0.262745</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.066667</td>\n",
       "      <td>0.023529</td>\n",
       "      <td>0.015686</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.015686</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.694118</td>\n",
       "      <td>...</td>\n",
       "      <td>0.188235</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.686275</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>0.039216</td>\n",
       "      <td>0.019608</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>230</th>\n",
       "      <td>0.007843</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.015686</td>\n",
       "      <td>0.062745</td>\n",
       "      <td>0.082353</td>\n",
       "      <td>0.188235</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.019608</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>...</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.082353</td>\n",
       "      <td>0.156863</td>\n",
       "      <td>0.011765</td>\n",
       "      <td>0.007843</td>\n",
       "      <td>0.011765</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.019608</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>231</th>\n",
       "      <td>0.015686</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.035294</td>\n",
       "      <td>0.121569</td>\n",
       "      <td>0.180392</td>\n",
       "      <td>0.341176</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>...</td>\n",
       "      <td>0.133333</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.172549</td>\n",
       "      <td>0.321569</td>\n",
       "      <td>0.027451</td>\n",
       "      <td>0.019608</td>\n",
       "      <td>0.023529</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.058824</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>232</th>\n",
       "      <td>0.019608</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.027451</td>\n",
       "      <td>0.054902</td>\n",
       "      <td>0.133333</td>\n",
       "      <td>0.278431</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.062745</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>...</td>\n",
       "      <td>0.113725</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.160784</td>\n",
       "      <td>0.286275</td>\n",
       "      <td>0.027451</td>\n",
       "      <td>0.023529</td>\n",
       "      <td>0.027451</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.062745</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>233</th>\n",
       "      <td>0.015686</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.027451</td>\n",
       "      <td>0.039216</td>\n",
       "      <td>0.117647</td>\n",
       "      <td>0.270588</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>...</td>\n",
       "      <td>0.090196</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.149020</td>\n",
       "      <td>0.278431</td>\n",
       "      <td>0.027451</td>\n",
       "      <td>0.023529</td>\n",
       "      <td>0.027451</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.058824</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>234</th>\n",
       "      <td>0.015686</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.003922</td>\n",
       "      <td>0.023529</td>\n",
       "      <td>0.035294</td>\n",
       "      <td>0.113725</td>\n",
       "      <td>0.247059</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.027451</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>...</td>\n",
       "      <td>0.101961</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.125490</td>\n",
       "      <td>0.207843</td>\n",
       "      <td>0.023529</td>\n",
       "      <td>0.019608</td>\n",
       "      <td>0.023529</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.027451</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>235 rows × 30 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     MWL_SWIRKaol_R  MWL_SWIRKaol_G  MWL_SWIRKaol_B  MWL_SWIRClay_R  \\\n",
       "0          0.000000        0.000000        0.000000        0.000000   \n",
       "1          0.027451        0.003922        0.003922        0.019608   \n",
       "2          0.019608        0.003922        0.003922        0.007843   \n",
       "3          0.027451        0.003922        0.003922        0.019608   \n",
       "4          0.015686        0.000000        0.000000        0.003922   \n",
       "..              ...             ...             ...             ...   \n",
       "230        0.007843        0.000000        0.000000        0.015686   \n",
       "231        0.015686        0.000000        0.000000        0.035294   \n",
       "232        0.019608        0.003922        0.003922        0.027451   \n",
       "233        0.015686        0.003922        0.003922        0.027451   \n",
       "234        0.015686        0.003922        0.003922        0.023529   \n",
       "\n",
       "     MWL_SWIRClay_G  MWL_SWIRClay_B  BR_Fe3_Fe2_FeOH_R  BR_Fe3_Fe2_FeOH_G  \\\n",
       "0          0.000000        0.007843           0.023529                0.0   \n",
       "1          0.003922        0.003922           0.000000                0.0   \n",
       "2          0.000000        0.000000           0.000000                0.0   \n",
       "3          0.003922        0.003922           0.000000                0.0   \n",
       "4          0.000000        0.000000           0.000000                0.0   \n",
       "..              ...             ...                ...                ...   \n",
       "230        0.062745        0.082353           0.188235                0.0   \n",
       "231        0.121569        0.180392           0.341176                0.0   \n",
       "232        0.054902        0.133333           0.278431                0.0   \n",
       "233        0.039216        0.117647           0.270588                0.0   \n",
       "234        0.035294        0.113725           0.247059                0.0   \n",
       "\n",
       "     BR_Fe3_Fe2_FeOH_B  BR_Anhydrite_R  ...  MWL_MWIRCarbonate_B  \\\n",
       "0             0.003922        0.000000  ...             0.007843   \n",
       "1             0.000000        0.949020  ...             0.203922   \n",
       "2             0.000000        0.909804  ...             0.180392   \n",
       "3             0.000000        0.992157  ...             0.262745   \n",
       "4             0.000000        0.694118  ...             0.188235   \n",
       "..                 ...             ...  ...                  ...   \n",
       "230           0.019608        0.000000  ...             0.058824   \n",
       "231           0.058824        0.000000  ...             0.133333   \n",
       "232           0.062745        0.000000  ...             0.113725   \n",
       "233           0.058824        0.000000  ...             0.090196   \n",
       "234           0.027451        0.000000  ...             0.101961   \n",
       "\n",
       "     BR_Kaol_AlOH_H2O_R  BR_Kaol_AlOH_H2O_G  BR_Kaol_AlOH_H2O_B  \\\n",
       "0                   0.0            0.007843            0.015686   \n",
       "1                   0.0            0.000000            0.956863   \n",
       "2                   0.0            0.000000            0.917647   \n",
       "3                   0.0            0.000000            1.000000   \n",
       "4                   0.0            0.000000            0.686275   \n",
       "..                  ...                 ...                 ...   \n",
       "230                 0.0            0.082353            0.156863   \n",
       "231                 0.0            0.172549            0.321569   \n",
       "232                 0.0            0.160784            0.286275   \n",
       "233                 0.0            0.149020            0.278431   \n",
       "234                 0.0            0.125490            0.207843   \n",
       "\n",
       "     MWL_Carbonate_R  MWL_Carbonate_G  MWL_Carbonate_B  BR_Fe2_FeOH_R  \\\n",
       "0           0.000000         0.000000         0.000000            0.0   \n",
       "1           0.058824         0.011765         0.007843            0.0   \n",
       "2           0.035294         0.019608         0.011765            0.0   \n",
       "3           0.066667         0.023529         0.015686            0.0   \n",
       "4           0.058824         0.039216         0.019608            0.0   \n",
       "..               ...              ...              ...            ...   \n",
       "230         0.011765         0.007843         0.011765            0.0   \n",
       "231         0.027451         0.019608         0.023529            0.0   \n",
       "232         0.027451         0.023529         0.027451            0.0   \n",
       "233         0.027451         0.023529         0.027451            0.0   \n",
       "234         0.023529         0.019608         0.023529            0.0   \n",
       "\n",
       "     BR_Fe2_FeOH_G  BR_Fe2_FeOH_B  \n",
       "0              0.0       0.003922  \n",
       "1              0.0       0.000000  \n",
       "2              0.0       0.000000  \n",
       "3              0.0       0.000000  \n",
       "4              0.0       0.000000  \n",
       "..             ...            ...  \n",
       "230            0.0       0.019608  \n",
       "231            0.0       0.058824  \n",
       "232            0.0       0.062745  \n",
       "233            0.0       0.058824  \n",
       "234            0.0       0.027451  \n",
       "\n",
       "[235 rows x 30 columns]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spectral_data = dataframe_KSL.iloc[:,4:]\n",
    "spectral_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABagAAAMWCAYAAAAUNGFuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABv+0lEQVR4nO39abgk51kf/nd1n+6zzqqZ0WjfbMn7vuINLxgcEi8oOPyCcYydBJNg9pD8zRYcB0ycmPhnQiAkEAFmMxhDCGBs49gYvODYkiwvsi1LM9JoFs1+zpylu6vq94L/Vd1PozPqGXXXM1P1+by6a+665vq+qKuq+zlP3Z3keZ43AAAAAACgZM3YAQAAAAAAqCcL1AAAAAAARGGBGgAAAACAKCxQAwAAAAAQhQVqAAAAAACisEANAAAAAEAUFqgBAAAAAIjCAjUAAAAAAFFYoAYAAAAAIIqZcU9stS8f67y0d/95h4FRZ7vumsng7yu97n1lxKEGdm55ZFGfWj8T9JIkKep+90Bpmai+4Xvd8L2t0Wg02q3Bo3p17Z6yIlEDw9ddp9UOejvml4r6/hOfLy0T1TY3d3VR99L+puf5PsEkbVu8oahXumubnue6Y1LO9rmu1Rwcr6/vLy0T1WfNjhjOdt0Nf48d535nBzUAAAAAAFFYoAYAAAAAIIqxR3xsm1ss6tMbq1MJA6OGXznuZeGroHkjLzsONTA61mNYnrvmmI7VL7ynqLN7Phf0vu0Nf1x2HGro7bufFxy//v99XKQkVNnZxnrAtBx82zcW9df/9B1Bb/fM0ujpMFFZngXHK5/7nUhJqLpdC9uK+sT6StDzPZYYzvVznx3UAAAAAABEYYEaAAAAAIAoLFADAAAAABDF2DOoe1k6zRzwoPpD1525SUBV5aeODOrTx4Nef2R2IUzD35kQt3wqRgxqbKkzHzsCFdX/4t1F3c3Du91abi465Rr+zAeTtN7vFnUrCfeipg3fJ7jw2UENAAAAAEAUFqgBAAAAAIhi7BEf7WarqNemEgX+rpmh666XGfEBVFPrmicUdT4y0mMu+UTZcaihrSNvfibX3BQnCLW10J6NHYGK6t17pqj7eTi28obW1rLjUHPDn/lgkhY7c0XdWz8T9JqNpOw4cM7soAYAAAAAIAoL1AAAAAAARDH2iI/1fq+o89yoBcrRywa/rO26A6pq/ad/oKh795wKevf3vZLH9H263Q2OX/62Xynq+Re8ruQ01NGlc9tjR6CiPv2py4p6pf+loHdPulx2HGpu7Se/p6jbv/CnEZNQNWk2mNfWaYVLfRtD63lwobKDGgAAAACAKCxQAwAAAAAQxdgjPlpNa9mUb6bZKup+lp7lTICLV742GK+QrYfjjLLRk2EKenl4pfVWfO6jXJ87dk/sCFRUrzEYlTU6MvCS1nzZcai7vk92TEd7aO3kxPpKxCRwfnz7AAAAAAAgCgvUAAAAAABEYYEaAAAAAIAoxp5BDTFkI3Piho3OkAO4WCXznaJuzq0FvebQ7EyYlnYS7lloL5mRCVTDbGPwnSFJwmfqQqM1ejpMVTLfjh2BilrprRd1lvscx8XHDmoAAAAAAKKwQA0AAAAAQBRjj/jop+k0c8CDGn41xUgPoKo+/JuLRX2ytSXoPdD9YtlxqKEv9k8Ex7//uauL+rvKDkMtLXTmYkegor7YGYxU2NjoBb2PnPla2XGoub/47cHnvFe8M2IQKicdWjtJjAjkImQHNQAAAAAAUVigBgAAAAAgirFHfMy3O0Xdy/pTCQOjgldTvKUCVNQ33PJ1RZ3vuyvo/eZbj5cdhxr6e81Lg+PvePMlkZJQV0/cfm3sCFTUP/3pa4r6t//1gaD3hPausuNQcy/5nlbsCFTUwsxsUe+Z3xb0Dq+eLDkNnDs7qAEAAAAAiMICNQAAAAAAUYw94iPL82nmgAfVag7+htLP0ohJAKYn2XVlUednloNeq/HlsuNQQ+3Rj3lbtkTJQX0d6y0/9ElwHvKv3V3UG3kv6M0m9mtRruSmG2NHoAaWe2vBcZplkZLA+DyRAQAAAACIwgI1AAAAAABRWKAGAAAAACCKsWdQp/lgZk1uHjUlGZ6V5LoDqqr37l8d1PecCnrH06TsONTQl1vd4Hj1Nz5a1AvfXnYa6mi1vxE7AhX1lf+xUtSneqtB73DHdUe5Nn73L4p64dU/ETEJVdPL+kW93g/n7Q+v58GFyg5qAAAAAACisEANAAAAAEAUY4/46KfpNHPAg8obxnoA1Xfkg+tFferktqC3kh4pOw41dChbC47vuu2Sot5VdhhqqTv0ajJM0mfSrUW9moYjPdZz33Ep1+Hb54t6R8QcVM/weNTtc4tB7/jactlx4JzZQQ0AAAAAQBQWqAEAAAAAiGLsER+zM+2i7qbhL4ImSTK5RDAkaQyuLeM+gKq64le/u6gv2/fFoHf193+s7DjU0PPzcLTM4992VaQk1NWNi5fHjkBFPWP2RFEvbMwGva9rbB09Habq2t/857EjUFFzM52inm22g167NfbSH0RjBzUAAAAAAFFYoAYAAAAAIAoL1AAAAAAARDH2IJpu2t+0l+dmAzMd5k4DdbD+jv93UB/Igt4DffdBpu/zzfXg+MQvfKKoF15bdhrqqJensSNQUXedGcyZ3sj2B72vJd2y41Bz6//5XUU99xuvipiEqsnywXeI1f5G0OunnrFc+OygBgAAAAAgCgvUAAAAAABEMfaID2M8AGA6PvqhvUV9shX+7fh474tlx6GGvpqeCo7/+L4rivq7yg5DLX3+1P6HPgnOw+/ODUYYLZ9ZC3p/cPKOov7F0hJRZx/9wKVF/fKIOaietf5gZNFGvxf0jE7lYmAHNQAAAAAAUVigBgAAAAAgirFHfMzOtIu6l/WnEgZGJY2kqL2WAlTVy+54a1H377k16P3hzbeUnIY6+o780uD4//nNl0RKQl3NzXRiR6CifmJuMNbjb+a3Bb0tM/Nlx6Hmhj/zwSTNDz1Hu2m4Zje8rgIXKjuoAQAAAACIwgI1AAAAAABRjD3io5l4JYDytVuDS7Sb9s5yJsDFq3/7h4o63/eVoHcmd+9j+o61wuP8q3cMDl5Qbhbq6ejqqdgRqKgjp5aKupcdCnot+7Uo2fBnvvaLb4iYhKrZMTu41y13185yJlyYPJEBAAAAAIjCAjUAAAAAAFFYoAYAAAAAIIqxZ1BDDK3m4G8oSRrOQc8bedlxAKYi/+JtRZ3dezDorWVmUDN9R5M0OE5v/1KkJACTdVdrrqh7G/2gN5O0Rk+Hqcq/cOvg4MXRYlBBw/ezpOE35Lj42EENAAAAAEAUFqgBAAAAAIhi7BEf6/3BK8Z5brQC5dgYvu5GRnq4DoGq+Mrb9xf1yfXZoPdAfrrsONTQF7LwOvvEby8V9Te8s+w0AJPziZmNol7tbwS9jdwYLcr1lf94b1E/8U0Rg1A5a2m3qI1D5WJkBzUAAAAAAFFYoAYAAAAAIIqxR3xs6cwXdTf1KhTlmJ/pFPVavxv0vLYCVMVjPv2fi7p/z61B79E331JuGGrpVdmO4Pi5f/SKSEkAJuuHl04V9V+sbg16L5y9suw41NzwZz6YpCcvXlXUJzZWgl4zScqOA+fMDmoAAAAAAKKwQA0AAAAAQBQWqAEAAAAAiGLsGdQQQ2JWEiUbvuby3JxzytG/8+NFne/7UtBby/tlx6GGjrXC4/yOTw4Onvmt5YYBmKDjJxaKOs0fCHpZ2WGoveHPfO3dN0RMQtWsDn1nGP0em/pey0XADmoAAAAAAKKwQA0AAAAAQBRGfAAMMdaDGGZuenZR92fng958cnvZcaihS9LwOHnE4+IEAZiwvZctF/XM/nCe0c68NXo6TNXwZz6YpFZj8/GoRqdyMbCDGgAAAACAKCxQAwAAAAAQxdgjPlb7G0XtFXjKcqa7XtRZ7ne2gWpa/uevL+rVA+Hrxvt6vbLjUEMfnVkLjl/4pncX9XW3v67kNNRRu2XyINPx/kOXFfVy/1DQ+6PuvUX9Y6Ulos6Wv2vwmW/nez8SMQlV8+nTdxf1Wr8bMQmcHzuoAQAAAACIwgI1AAAAAABRjP0u3VyrXdRrycZZzoTJmWkOXnXvZUbLANW05b/9SlHP33Nr0Lvm5ltKTkMdPb8/Hxzvfdc/iJSEuuql/dgRqKhv3HuwqN9x70LQm2u2R0+HqdryS7/y0CfBecgbg/USY3m5GNlBDQAAAABAFBaoAQAAAACIwgI1AAAAAABRjD2D+sT6SlGbZ0NZumkvdgSAqfu2p/5AUa/n4RzWT5z8atlxqKGfWrs9OP7d1xwr6j+793UlpwGYnG89OHiu3nP6SNC7t3m07DjU3D9+6g8W9Xv2/2HEJFTNyfUzRT08j7rRsIbHxcEOagAAAAAAorBADQAAAABAFGOP+Oi02kW90e9OJQyMaiaDv6FkeRYxCcD0/OrfGzxXu/tWg94333ZZ2XGooZduuTE4fsvVXnunXNvmFmNHoKLenF5e1N+zcDzoXTe/p+w41Nz/+Kb12BGoqO1Dz9GVbnidWUvhYmAHNQAAAAAAUVigBgAAAAAgirFHfEAMo78+C1BFzZ1LRT1zeiPotcoOQy3NJeGV1tnp+Uu59s7viB2BGtrZWogdgZpp7jDOiOlrNcO9qFlqxAcXPjuoAQAAAACIwgI1AAAAAABRjD3iY26mXdTdtDeVMDCqmQz+hpI1wtdS8tzrx0A1zP7AzxZ1655bg95lH7ul5DTU0ZN77eB44Qe/M1IS6irNvX7MdFzZXi3qThp+/b2macQH5Zr94bfHjkBFbW0P7mcb/XDNbnhdBS5UrlIAAAAAAKKwQA0AAAAAQBQWqAEAAAAAiGLsGdQQQzNJijobGTmdDPUALmbpobsGB0f2B71unpachjo6PbJlIT+8/8FPhCm569TB2BGoqNlOv6iT9fD7w5a8VXYcam74M1979w0Rk1A1p7pninq1txExCZwfO6gBAAAAAIjCAjUAAAAAAFGMPeKjlVjLpnxJwxgPoPpaewevePbXl4NeJ/H6MdO3NQuPk0uvjhOE2srz/KFPgvNw7Mx8Uad5eLOb812Dkg1/5oNJmm21i9o4VC5GVp0BAAAAAIjCAjUAAAAAAFGMPeJjpbte1F7Boyz9LC1q1x1QVZ994g8V9Yl0Nuh9IT9Ydhxq6Hebx4LjXa9+f1G//PDrSk5DHc3OdGJHoKJ+c27wqvup42eC3rvXvlzUP1ZaIursM0/84aJ+5v1/EDEJVbPW7xb16Ije0fFGcCGygxoAAAAAgCgsUAMAAAAAEIUFagAAAAAAohh7BrWZNcSQue6AGphtD+btz6Vp0EsayejpMHHtkT0L8w3PX8q11JmLHYGKmhuaxTr6TF3urZUdh5qbnUkf+iQ4D8MzqId/ywsuFnZQAwAAAAAQhQVqAAAAAACiGHvEx+xMu6hXu14XoBxJMngNL8/ziEkApuem3399UWf7vhj0rvv+1bLjUEMvbuwIjp/z9qsiJaGubly6PHYEKuofrg++u/5RZyHoXTG3s+w41Nyj3vudsSNQUTvnlor62Npy0Msb1lK48NlBDQAAAABAFBaoAQAAAACIYuwRH2nm19wp3/AvbXstBaiq/PA9g4PDB4Nev+H5y/StJuEzNr//4CZnwnSsphuxI1BRx/PBqMp+Ho6q7CRjfx2GiQg+88EEDY9HhYuRHdQAAAAAAERhgRoAAAAAgCjGfqep0xqc2k17UwkDAHU084yXF3V/z61BbynZV3Ia6mhvGr4W2nz218cJQm0d2jgZOwIV1RoaE9hshPe6a1tbyo5DzQ1/5oNJ6meDEUZpbkQgFx87qAEAAAAAiMICNQAAAAAAUVigBgAAAAAgirFnUA/PsMnz/CxnwuQkyWBOXDIyMw6gKvq3f6io831fCXpncr/7wPQdaoWf7fLbPzU4eMHryg1DLZ1cPxM7AhW12hzsycoa4b1utZGOng5T1f/4HxR1++U/HDEJVbMx9Ftx1uy4GNlBDQAAAABAFBaoAQAAAACIYuwRH61ksJY9PHYBAHh4km27izrfdjTotYw3ogQL+ch1tnNnnCDUVtP3C6bkxvZKUc/22kHv+sZC2XGoufz4A7EjAFyQ7KAGAAAAACAKC9QAAAAAAEQx9oiP1d5GUftFUMqS5VlRu+6Aqnrrze8p6pUkC3qfXdlXdhxq6Fe7dwXHd/3o4Jfg3/XtZaehjtaGvmvAJH1fulrU9585HvR+L7mzqH+6tETU2dt+8t6ifsvr4uWgela660U9vI4CFws7qAEAAAAAiMICNQAAAAAAUVigBgAAAAAgirFnUM80W0WdZulUwgBAHb1m9kRRn1ieD3of7CyVHYcaurqzMzi+ec3sQsq10JmLHYGK+sbWpUX9pc6BoJc0krLjUHOvSM/EjkBFbekMvkMsd9eCnt/z4mJgBzUAAAAAAFFYoAYAAAAAIIqxR3wstmeLupv2phIGRrWbg0u0l/UjJgGYnhv+6ueLun/PrUHvpptvKTkNdfSqbEdw/MzffUmkJNTVanc9dgQq6pWdwRitX2mHY7R2to3RolxP+Ow7Ykegor5h52OK+s+Pfz5iEjg/dlADAAAAABCFBWoAAAAAAKIYe8RHL0unmQMeVJpnsSMATF3/zo8Xdb7vS0FvLTfeiOk71gqP86/eMTh4QblZqKfZmU7sCFRUqzX4PtFMkqC3rTU/ejpM1fBnvvbuGyImoWoO9leK2vodFyM7qAEAAAAAiMICNQAAAAAAUYw94sOoBQCYjmRusajzuYWwV3YYaqmdj/zD3FyUHNTXVUu7YkegotbX20Wd5eHNbmtitAwl21iNnYCKOtg9UdT91IgPLj52UAMAAAAAEIUFagAAAAAAorBADQAAAABAFGPPoG4l1rIBAGrB7EJKttxfix2BilrtDmZQj/6u0g4zqClZsm137AjUQN4Y/XERuPBZdQYAAAAAIAoL1AAAAAAARDH2iI/h16Hy3OsCADAp2Sc+OKjvPRj0VrJu2XGooYPNcKRH/xO3Dg5eW24W6unwyonYEaioextzRd3Pw3vdSqNfdhxqLvv4BwYHT3tVvCBUzrH15aJOs+wsZ8KFyQ5qAAAAAACisEANAAAAAEAUY4/4WOttTDMHPKjMaBmgBv7Z2+8v6vWR14+/sHxv2XGood8786Xg+Av/6/Ki/l+/UHYagMn599ndRf3A6qmg9+H8y2XHoea+6z8eKupff1PEIFTO6Y3Voh5eR4GLhR3UAAAAAABEYYEaAAAAAIAoxh7xMd+eLWrjPihLuzm4RHuZX9kGquk/P+pYUZ853Al6N/d2lx2HGnrW4rXB8U/MrcUJAjBhb2zfUNT/bv500Nvd2Vp2HGru5x5zNHYEKmrL7HxRL2+En+PyhnGpXPjsoAYAAAAAIAoL1AAAAAAARGGBGgAAAACAKMaeQd1utoraVELKkiRJ7AgAU7f9N361qJfuuTXoXXvzLSWnoY5e3F8Mjq/+b6+IlARgsm6+6d6iftfnwnvd7vaWsuNQc9t+7Vcf+iQ4D8O/3zUqaVhX4cJnBzUAAAAAAFFYoAYAAAAAIIqxR3xkeT7NHABQW72P/k5R5/vvCXqns42S01BH986En/Oyj35wcPCC15Ubhloy1o1p+cQdlxf1mf7ng16vnZUdh5ob/szXvvnNEZNQNXnDmh0XNzuoAQAAAACIwgI1AAAAAABRjD3io91sTTMHPKhOa3CJ9tJ+xCQA09N+/j8q6v49twa9rT93S8lpqKOr+uF4hebzXxIpCcBkPf7So0U9d7AT9B7b3ll2HGpu+DMfTNLM0JqdsVlcjOygBgAAAAAgCgvUAAAAAABEYYEaAAAAAIAoxp5BnebZNHPAg+qeZe503shLTAIwPemhuwYHR/YHvW6elpyGOjo9smUhP7z/wU+EKclzn+uYjpnZwXO0OTKXdXc+9tdhmIjhz3zt3TdETELVpNlgzc4zlYuRHdQAAAAAAERhgRoAAAAAgCjGfqeplw1ejfK6AGXJhkbLGOkBVNX6z/zMoD4QjtQ63HfvY/o+31wPjk+848NFvfDqnyg7DsDE/NHBy4v6ZPczQe8Lc2fKjkPNrf/s24p67jdeGjEJVZMZy8tFzg5qAAAAAACisEANAAAAAEAUY4/4aCWDtexk5NePYVpmmq2iHv5V2kbDyA+gOjovfmpRz+y7N+gt3XWi7DjU0KWNTnC87UWXREoCMFnPmzlZ1Iszc0Hv8uZcA8rUeeFTYkegoppnWbMzppeLgR3UAAAAAABEYYEaAAAAAIAoxh7xMdtqF/VKY20qYWDU8GiZUUnDqBmgGma+4TVFne67Peht/Z/vKTsONXRVGj5vW89/QaQkAJO1bcfgu2vrTHive2TaHj0dpmr4Mx9M0vzMYFzb6IBAY3q5GNhBDQAAAABAFBaoAQAAAACIwgI1AAAAAABRjD2DerW/Mc0c8KDW+t2izht5xCQA03Ps27+/qJcfmAt6+7vLJaehjj7aOhMcP/df/llRP/auN5UdB2Bivnz/rqLeyI4GvV/u3lXU31NaIurs+Hd8X1Hv/chHIiahak6srxR1lmcRk8D5sYMaAAAAAIAoLFADAAAAABDF2CM+FmZmi3qtZ9wH5Zif6RT1qusOqKgdv/Dmot6+//NB74o3frDsONTQ07PF4Pi6n3lypCQAk/XMm08X9fb3hve6x85eWnYcam77z/9o7AhU1I65paLeSHsRk8D5sYMaAAAAAIAoLFADAAAAABDF2CM+IIaskceOADB1rb03FHV/fTnodZJW2XGooa0jP/aeXHp1nCAAk5Zt/n1iV9LZtAfTMPyZDyapm/WLOsuto3DxsYMaAAAAAIAoLFADAAAAABDF2CM+ekOvC0BZeungusuN+wAqqn/nx4s63/eloLeWe/4yfcdGJsnkX71jcPCCcrNQT0mSxI5ARS3fulHUG1kv6K3madlxqLnhz3zt3cZ9MDnB2okRH1yE7KAGAAAAACAKC9QAAAAAAERhgRoAAAAAgCjGnkG90l0vavNsKMvwHCWAqnrnP/7Toj6VZEHv9pX9Zcehhn5j467g+OC/Gzx/3/GGstNQR1tnF2JHoKLedP9SUd93Zl/Q+0iejZ4OU/VzQ5/53rz/NRGTAFxY7KAGAAAAACAKC9QAAAAAAEQx9oiPhfZsUa9016YSBkY1k8HfUPKG0TJANX3fb76sqPN9Xwp6t/2bU2XHoYZeM3tDcPwvf3xvpCTU1Xq/FzsCFfWWmcEYj9vntwe9ly8+suQ01N0P/Po3xI5ARS3MDNbsTiWrEZPA+bGDGgAAAACAKCxQAwAAAAAQxdgjPppJMs0c8KBazcHfUPpZGjEJAADT0k2N+GA6duwZvOreXmkFvT1Za/R0mK4TR2InoKLWPUe5yNlBDQAAAABAFBaoAQAAAACIwgI1AAAAAABRjD2Der0/mGeT5/lUwsCoNMuK2nUHVNW+17+7qE8szwe9/dmxsuNQQx/NTwTHT3/zRlF//RvKTkMd+ZzHtLz/3iuK+lTvs0Hvj9oHi/pHSktEne37kQ8X9aNe6apjcs701os6y7OznAkXJjuoAQAAAACIwgI1AAAAAABRjD3iY7E9W9TdtHeWM2FyZpqtou5lXv0EqumGv/r5ou7fc2vQu+nmW0pOQx29KtsRHD/zd18SKQnAZH3L8w4U9c99eEvQ29acKzsONTf8mQ8mae/C4LPcfStHIyaB82MHNQAAAAAAUVigBgAAAAAgirFHfPSydJo54EHlDWM9gOpL991e1Pn+O4Peet4vOw41dHpky0J+31fjBKG2Lt9ySewIVNTq/kHdz8PvtItJu+Q01F3/s+8v6vZL/0XEJFRNL9v8O0OWW1fhwmcHNQAAAAAAUVigBgAAAAAgirFHfKz2Noo693oAJekPjZZx3QFV9dab31PUK0kW9D67sq/sONTQr3bvCo7v+tFeUb/r28tOQx3dv3wsdgQq6p/eO1/U950Jn6nLvbWy41BzP/PPPlbUb9lnxAeTc3j1ZFH3UiMCufjYQQ0AAAAAQBQWqAEAAAAAiMICNQAAAAAAUYw9g3p2pl3Ua73sLGfC5DSTwd9QsobrDqimN//09UWd77836N3xX64oOw419C1z1wfHP/RaswuBavjB7lJRf2FuS9B7/tINZceh5v7NW6+LHYGK2jq7UNQn1lYiJoHzYwc1AAAAAABRWKAGAAAAACCKsUd8tBJr2QAwDe2XfmdR9++5Neht+6+3lJyGOrqqnwTHrW/+h5GSAEzWZUtninrmVCvoXdeYKzsONdf+ptfHjkBF9bM0dgR4WKw6AwAAAAAQhQVqAAAAAACiGHvEB8TQTAavHGd5xCAAU9T/1B8Vdb7va0FvJe+WHYcaOtQKH7LZX/354OCZ31pyGoDJSZLNv0R0G75gUK7e//mtom5/649FTELVZPngfpa7t3ERsoMaAAAAAIAoLFADAAAAABCFBWoAAAAAAKIYewb18CxgAGByZp7x8qLu77k16C0l+0pOQx3tTcPPec1nvTBSEoDJ2nPTSlHPfbod9K7IWmXHoeaaj31O7AhUVJ6bO83FzQ5qAAAAAACisEANAAAAAEAUY4/46GXpNHPAg8q8pgLUwNFveUNRLz8wF/Tu7i6XHYca+tDMmeD4md/x7qJ+7F3/uOw4ABPziU9fXtTLvS8Evd+eub+ov6e0RNTZ8X/6lqK+/K+N02JyrNlxsbODGgAAAACAKCxQAwAAAAAQxdgjPiAGv0QL1MG2739pUW/df0/Qu+Q/7Ss5DXV0Y2MhOL72dYuRklBXszOd2BGoqOc871BRb/3L8F7XSXwdplzb3/DU2BGoqGu27C7qfcsPREwC58cOagAAAAAAorBADQAAAABAFGO/05RmWVEbu0AMrjugqs78tz8b1IfD19yPp+tlx6GGvtxYDY7v/63B8SN/vOw01FGWZw99EpyHP/7rK4r6RPe2oNdK7NeiXAf+0+eK+pFviBiEyjm2sVzU/SyNmATOjycyAAAAAABRWKAGAAAAACAKC9QAAAAAAEQx9gxqiCFvmDsNVF9/dfD34u56+Gg2f58Yuhs+IlKubbOLsSNQUbNDz9EkSYLejhnXHeU6fXoudgQqaqbZih0BHhY7qAEAAAAAiMICNQAAAAAAUYz9/manNTi1m/amEgZGJY1k+ACgknb8wpuLevv+zwe9K974wbLjUEPPzsLX3K//z8+NlIS66mdp7AhU1Dd+w+Gi/vfvXwp6z2/tLjsONbdr90rsCFRUL+3HjgAPix3UAAAAAABEYYEaAAAAAIAo/EQ7AETW2ntDUffXl4NeJ/GL3EzftjQPjpMrHhEpCXW1YYQgU7L6tcH4mG4evgI/n5shSLkW9rjXUb48zx/6JIjMDmoAAAAAAKKwQA0AAAAAQBRGfABAZP3bP1TU+b6vBL0zuVdBmb57Rz4RZn/1/sHBU15RbhhqqZ+lD30SnIfVU52izvKVoNc34YOSzWyxRxDgwbg7AgAAAAAQhQVqAAAAAACisEANAAAAAEAUY8+gbiYGdFG+ZOi6y/M8YhKA6Zl5wouLur/1kqC3mNxZdhxqaG8afs5LnvCMSEmoqzTLYkegoi79e0tFvXjLmaB3Y9f3C8rVvnFX7AhUVGa9hIucHdQAAAAAAERhgRoAAAAAgCjGHvGx3u8VtVELlCXLB697uu6Aqjr6LW8o6uUH5oLe3d3lsuNQQx+aCV97f+br31PUj73rdSWnoY7m2p3YEaioT//PVlEf768EvZ/t7C/qV5eWiDr7s18Z3Ov+4U9HDELltJqb7z9NjOzlImAHNQAAAAAAUVigBgAAAAAgirFHfPSzdJo54EEZ6wHUwe/cc2VRn2qG970HugfKjkMNfan7QHD8W42ri/qtZYehltZ73dgRqKh3zw2eq6ePrwa9tb7rjnJ9YHYwOvUfRsxB9ax014t6eFQqXCzsoAYAAAAAIAoL1AAAAAAARGGBGgAAAACAKMaeQT070y7qtZ55NpSjmQz+hpI3zKMGqumNf3BzUWf7vhj0PvX9J8qOQw19x8y1wfE//+lrH/Q8mJad81tiR6CifmLP8aL+2Nr2oHdZZ0fJaai7t3/DqdgRqKg9C9uK+tAZ3x+4+NhBDQAAAABAFBaoAQAAAACIYuwRHwDAdOSnHhgcnApfyUuNN6IEvWTkH5aXo+Sgvq5d3BM7AhV1333bi3o9PR70tjQ7Jaeh7h745KDeGi8GFXT57M6iPrbucxwXHzuoAQAAAACIwgI1AAAAAABRGPEBAJHl93y5qLN99wa9bp6WHYcaOpVkwXG2/0CkJNTVvWtHY0egolpD97dmEs4zWkzaZceh5o4dWyjqGyLmoHoObgzGBPbSfsQkcH7soAYAAAAAIAoL1AAAAAAARDH2iI9WYi0bAKah+YyXFnWy9/NBb/FXP1h2HGroqjT8nNd6/gsiJaGuFlqzsSNQUY987smiXvhIeJ09urHQgGlKRsbKdGayTc6Eh2d4hFFzZP0uz/Oy48A5s+oMAAAAAEAUFqgBAAAAAIjCAjUAAAAAAFGMPYM6zc1KAoCp2DhTlPn6WtDKGmbGMX1r4YjMRn7scJwg1NbR9dOxI1BRBz8zmDO9lobX2WrDd1yma3T278JCN1ISqu7k0PeJfpZGTALnxw5qAAAAAACisEANAAAAAEAUY4/4WOttFPXoayowLZnRMkANvPXm9xT1ShLe9z67sq/sONTQr3bvCo7v+tFeUb/r28tOQx2tDn3XgEn6V0OX1uG1k0HvD5uDe9/PlJSHevveM4OZWn8WMQfVs9wdjAkcXUexhsfFwA5qAAAAAACisEANAAAAAEAUY4/4aDcHp3bT3lnOhMlJksErUF5LAarqBeuDX9o+0WwFvT+ZmSs7DjW0p701OH7OxtgfEWEiZkbufTApz0t2FPWnZ2aD3mpqtAzlenxze+wIVFQrGew/HR3xMbyuAhcqO6gBAAAAAIjCAjUAAAAAAFFYoAYAAAAAIIqxBwwudQYzME+s96cSBkYNzyPsZ2nQSxrmKAHV8PV3/HRR9++5Nej9zs23lJyGOnpttis4fuX7vyNSEupqsT370CfBefjeNw/ub3/wlj1Bb3trvuw41Ny/3PFA7AhU1MLQc3S5m53lTLgw2UENAAAAAEAUFqgBAAAAAIhi7BEf62lvmjngQaXZ5q+m5I28xCQA09P9xZ8s6nTfkaB3ODVWi+m7ox0+b1/16784OHj7N5Wchjp6xNJlsSNQUff95zuL+kRvJei17NdiyppJeI199tBgzMy1JWeh2oZHfJzprUdMAufHExkAAAAAgCgsUAMAAAAAEMXYIz4ghiRJijrPjfQAqqn5dc8v6uSKrwS9pfffUXYcamhX3gqOm496RKQk1NUXTt0bOwIV1ZkbjMpKGknQ2z2zWHYcaibLwxFaT9rzQKQkVF13aCxgZu2Ei5Ad1AAAAAAARGGBGgAAAACAKMYe8dFKrGVTvrON9TDyA6iM2YWhejZoNUdeR4ZStNuxE1AzM83WQ58EE7Y3mY8dgZrZdkM3dgQqKm8M1keslXAxsuoMAAAAAEAUFqgBAAAAAIjCAjUAAAAAAFGMPYN6o98ravNsKIs5SkAd7Hv9u4v6xHI4D3N/dqzsONTQxxunguNn/8haUT/3tWWnoY5Ob6zGjkBF/e+je4t6uX8w6B3K10ZPh6n6809cWdTfFjEH1ZNmWewI8LDYQQ0AAAAAQBQWqAEAAAAAiGLsER8QQ9JIinp43AdAldy+vKOoTzXDvx2v9bplx6GGjqRnguNb29uK+rllh6GWstyryUzHJ1uD8TGrvY2gd6i/XHYcau5YK3nok+A8bKRDY3mtnXARsoMaAAAAAIAoLFADAAAAABDF2CM+ZmfaRd3L+lMJA2eTJF6HAqrpm//F4O/F6b4jQe+WP90xejpM3JPbu4Pj1zzrQKQk1FW7ZfIg0/HGoakefzm7FPSe0dlbchrq7luuvS92BCrqCduvLerPndoX9PLcyA8ufHZQAwAAAAAQhQVqAAAAAACiGPtdum46GOvh9QDKMvzrs647oKpO//ngdc/lB+aC3mq6MXo6TNzxvBscH/ibhaLeXnIW6ilpGOXGdHy0M1/Uq2vhM/W23tGy41Bz77vnyqL+7og5qJ4vnL63qNd73bOcCRcmO6gBAAAAAIjCAjUAAAAAAFFYoAYAAAAAIIqxZ1B3WoNTu2lvKmFgVDCP0GhCoKJ2/Nz3FfW2fV8Mepd9/8fKjkMNPT1bDI6ve8tjIyWhroa/a8AkfecNg995+LXPb4mYBBqNb338/tgRqKjZVruou61+0MvyrOw4cM7soAYAAAAAIAoL1AAAAAAARDH2u3SpVwKIIEkGcz3yPI+YBGB68sP3DA4OHwx6/YbnL9O3mow8Y48fjxOE2soaPucxHYe+trWou9ly0FvPjK6kXIe/PBgzsytiDqonH3qOjo70yKylcBGwgxoAAAAAgCgsUAMAAAAAEMXYIz5aibVsAJiG/NCBos4OHQl6PSO2gBqYa7VjR6Citu1YK+rWaviddmtrruw41MzwyMpGo9E4vLJQ1I8tOwyVNus5ykXOqjMAAAAAAFFYoAYAAAAAIAoL1AAAAAAARDH2DOqV7mB2V57nUwkDo7Kh2auuO6CqXv6Ttxd1N0+D3m0n7yk5DXX088u3Bscf+NnLi/rDbyo5DLV0fG05dgQq6u8fHFxb+5bD33k4tHqi7DjUzOh32O/O7i7qO8sOQ6UdPnOyqNMs3fxEuEDZQQ0AAAAAQBQWqAEAAAAAiGLsER/z7dmiXuttTCUMjEoaSVHnDSM+gGr6o3d+fVHn++8Jeq/8T2M/quG8fdfWJwXHP/zd7ThBqK3ZmU7sCFTUb87vLupvSXtB72VLjyw7DjX36ddeETsCFfX4ndcU9Z2nDgS9zLhULgJ2UAMAAAAAEIUFagAAAAAAohj7veFWYi2b8iXJYMTH8LgPAACqI8uz2BGoqPXu4Cvv6MjAtu8XlGz546eLekvEHFTP4fWTRd1N+/GCwHmy6gwAAAAAQBQWqAEAAAAAiGLsER8AwJQcHPzSdn7oaNDq5mnZaaih1UY4XiE/djJOEGqrn7nXMR37s4Wi7mbha+/LudfgKdeX7tpV1JdHzEH1rPTWizrP87OcCRcmO6gBAAAAAIjCAjUAAAAAAFFYoAYAAAAAIIqxZ1CnefbQJwEA5+zM//rcoD7cCXor2fro6TBx9zXC6+zg/14r6ke+vew01JF5mUxLlmzeu733QHlBoNFoXL5tJXYEKmrX3NaiXu93IyaB82MHNQAAAAAAUVigBgAAAAAgirFHfLSbraJeO8t5MEnNZPBOXubNT6Ciln7kNUW9uO8rQW/3T95Rdhxq6En5YnB8xfdeHSkJwGT9vdesFvXbfn1r0NvanCs7DjW3bY/VFKYjb1gw4eJmBzUAAAAAAFFYoAYAAAAAIIqxR3xkflmbCJrJ8N9Q0mg5AKZq/cxQvR60vK4HAOfvzt8aPEdP9VeD3pbWfNlxqLlt3/bY2BGoqG7WL2rrd1yM7KAGAAAAACAKC9QAAAAAAEQx9oiPjX6vqHOvC1CSXjp4TcV1B1TVp9/wV0V9PG8HvXuyo2XHoYY+lB8Ljh/3k4Pn78veWHYagMn5tZm5oj65cSbozSStsuNQc7/0s6eK+ge+J2IQKufY2nJRZ3kWMQmcHzuoAQAAAACIwgI1AAAAAABRWKAGAAAAACCKsWdQz7QG87l6Wf8sZ8LkJEkyqBvJWc4EuHg95Ye3FXW2/1DQu/T3t42eDhP3tOaO4Ph5b+hGSgIwWc9bH3yPfW97Nug9eeHKsuNQM8PfZxuNRuNpvfVISai6VnNo/2kaLwecLzuoAQAAAACIwgI1AAAAAABRjD3iAwCAamqPjNFK5uciJQGYrOc/9kBRb/ncfNC7PHGvY7razXDJ5f2znaJ+XtlhAC5gdlADAAAAABCFBWoAAAAAAKIw4gMAYlvbKMp8tRe0skZedhpqaLWRBcf56bVISQAm664vXFLUa+nXgt5qnpYdh5rJRz7HHW50IyWh6ppD49qSkdFtcDGwgxoAAAAAgCgsUAMAAAAAEIUFagAAAAAAohh7BnW72SpqUwkpS6c1uEQ3+uFc1tF5XgAXq9bLbi7qZN8Xg96293+s7DjU0FVZKzhuft2zIiUBmKzHvXow83fb7ywEvWf0Z8uOQ8300n5w/M6fuiFSEqruSduvK+rPnPzaWc6EC5Md1AAAAAAARGGBGgAAAACAKMYe8bE+NF4hz41WoBz9LC1qIz2Aqrr95bcU9Yk0fN34a/kDZcehhv64cSw4vvqNtxX1y19VdhqAyfm1399e1IfXw9fe3zH7laJ+fVmBqLVP/evBNfj1r4uXg+o5la4WdZplEZPA+bGDGgAAAACAKCxQAwAAAAAQxdgjPjqtwandtHeWMwGAc/H4X35BUef7wtePr/6ZU2XHoYaen+wMjl/yY9siJQGYrH/8rPuK+hf/cmvQ29KaLzsONfe4pxyOHYGKWu6vFXWWG/HBxccOagAAAAAAorBADQAAAABAFGOP+Ei9IkAESSOJHQFg6rI7vzCo9x0MeutZv+w41NCpJPycl351f6QkAJO1dv/g+0Q/T4NeM/Fdg3J1rluMHYGK2tFeKuqDyYmg517HxcAOagAAAAAAorBADQAAAABAFBaoAQAAAACIYuwZ1K3EWjblyxt5UY/Oox7uAVzMmjc9pqiTubmgN9f8StlxqKFdeSs4bj36hkhJACarvWUwd3r0+8SVM1vLjkPNtJrh8/XX/+iSov6Xv1B2GqrsgY1TRZ2N/IZcZumEi4BVZwAAAAAAorBADQAAAABAFGOP+IAYsty7KAAAwPkZedM9sL3RLi8ItTQ6KvXPmyeL+l+WnIVqe2DtdFGn2VlufHCBsoMaAAAAAIAoLFADAAAAABDF2CM+mkny0CdBiUZ/hRvgYjXzjJcXdX/PrUFvKdlXchrqaG8/HKnVfMaLIiUBmKyFZ19W1PN3HAh6j8tmy45DzbSa4Z7A65pLkZJQdcPXWi+NGATOkx3UAAAAAABEYYEaAAAAAIAoLFADAAAAABDF2DOo4UKTN/KHPgngItC/8+NFne/7UtBby/tlx6GGjrXC33XI7vjE4OApryg5DcDknP7gkaJeTTfCXtv3CaYry8Nr7O5sJVISqi7NstgR4GGxgxoAAAAAgCgsUAMAAAAAEMXYIz5GX00BACZkY3WoDl8/zowzogS9ZOQfVlcf9DyAi83sJWlRt+4K92dtzUdvfjBZ/SwNjucTU1aZjmbifsbFzQ5qAAAAAACisEANAAAAAEAUY79fkuZ+EZS4cq+5AxXVe897i7q7LxytcCL1uh7Td3ezGxyv/uFni3rhjWWnAZicz33q0qJe6X856L23faSov6e0RNRJOjLi45+sz0VKQtVdMrelqA+vnowXBM6THdQAAAAAAERhgRoAAAAAgCjGHvHRSqxlU748N9YDqL7WI64o6tnOoaA3e/t62XGooS15Kziee9zOSEkAJuvGG44W9fxXO0Gvm/fLjkPNJEk4qu2XZs8U9cvKDkOlHVtfLuo0M6KXi49VZwAAAAAAorBADQAAAABAFBaoAQAAAACIYuwZ1M2R2UlQhuGZXeZRA1XVfPIzizrZ+ZWgt/AHd5Qdhxq6KgtnUDe/7lmRkgBMVmdbWtSjv6u0p7VYdhxqrm2PIMCDcncEAAAAACAKC9QAAAAAAEQx9oiPzHgFAJiO9TND9XrQyhuev0xff/Qfet0YMaixZmLfDNORdQcjA0e/0+5IZsuOQ82M3tvO5L1ISai6NMtiR4CHxSdBAAAAAACisEANAAAAAEAUY4/4gNiSRhIce+0dqIy5xaF6LmiN3vtgGkY/ECZbd0TJQX21mvbNMB0zWwd1MwmfqQtJq+Q01M3oNbe12YmUhKpLEt8ZuLj5JAgAAAAAQBQWqAEAAAAAiMKIDy5oeW6MB1ADswtD9WzQahrxQQnao4/b9uyDngdwsXngC/NF3c1Wgt5cbr8W0zU6qq2XZ5GSUHXWTrjYeSIDAAAAABCFBWoAAAAAAKKwQA0AAAAAQBRmUHNBS5LBzC4zlYDKOnF4qD4RtNKGWYVM39rIqPP80P44QaitXtqPHYGK2vP49aKe/ct20NtqvxYlm01asSNQUcNrJ6PyhrUULnyeyAAAAAAARGGBGgAAAACAKMYe8dE8y+sCEEPScE0C1TDzrFcWdX/vrUFvKbm33DDU0t5++Opn88kvjpQEYLJmX/T4ol7469uD3nPWu2XHoWZGRytcnsxGSkLVtZqD/ae9NOxZO+FiYAc1AAAAAABRWKAGAAAAACCKsUd8ZLlf/aR8+dB155dngarq/sbPFnV23+GgdzRdLTsONXRHJ3zG/v3f/sXBwVteWHIa6qjVbMWOQEXd8467ivpYdzno/deF7UX9orICUStplgXHX8vOREpC1Y1ea8OspXAxsIMaAAAAAIAoLFADAAAAABCFBWoAAAAAAKIYewZ1M0mmmQMAaqvzmn9d1P17bg16u95zS8lpqKPr++GehdY3/f1ISairLN98diY8HNf8k51Fvf0Xls9yJkze6L1te9KJlIQ6SxrW87jw2UENAAAAAEAUFqgBAAAAAIhi7BEfWZ5PMwc8qGRotEzuGgQqKt13e1Hn++8Meut5v+w41NDqyJaF/PD+OEGoLZ/zmJaj7z1Y1Mv9taB3emaj7DjU3NPT+dgRqKhOa7C810t9f+DiYwc1AAAAAABRWKAGAAAAACCKsUd8pH5ZGwAAgIvI+mq7qPNGOEpmKemUHYeaGx2pBZNiVBYXO7dHAAAAAACisEANAAAAAEAUY4/4aCXWsgFgGvJTDwwOTp0IemnD63pM31oy8g/Hj0bJATBpnzu1s6jX0gNBbz3vlx2HmtuRxk5AVSXJ6Ic5uLhYdQYAAAAAIAoL1AAAAAAARGGBGgAAAACAKMaeQQ0x5LnZq0D1zTzhxUXd33pJ0FtM7iw7DjV0VS983jaf/c2RklBXW2YXYkegok62BnuyRr9bbGvOlh2HmhmdC3yo5fst05HmWewI8LDYQQ0AAAAAQBQWqAEAAAAAiMKIDy5ow69EGfcBVFX/zo8Xdb7vS0FvLe+XHYcaOjQTvoKc3fGXg4PHvLDkNNTRam8jdgQq6lHpWlG3m+HX3/lGq+w41Mzod9jre5GCUHnNRvLQJ8EFzA5qAAAAAACisEANAAAAAEAUY4/4aCZeFyCuZOSVlbxh5AdQDTM3Pbuo+7PzQW8+ub3sONTQ1pEffk8uvTpOEGorzdLYEaioK/eeKurOgfDr76PyubLjUHNHZ6yrADwYO6gBAAAAAIjCAjUAAAAAAFFYoAYAAAAAIIqxZ1Cv93vTzAEPKsuzTXt5bgY1UA3L//z1Rb16oBX09vU8f5m+j86sBcff+K9vKeorPvG6ktMATM6vnNhT1Ce69wW93565p6j/VVmBqLXbmuuxI1BR3bRf1H6vi4uRHdQAAAAAAERhgRoAAAAAgCjGHvExN9Mu6m7qdWPK0UwGf0MZHfeRJEnZcQCmYv57vr2o5/Z9Jejt+sk7yo5DDT0pmw+Od7zpOZGSUFdXbtkVOwIVdVNv8J2h0xz76y9MxZHciA+mb6YZjgzMjEflImAHNQAAAAAAUVigBgAAAAAgirHfcfJKAABMR/eW3ynq3gPhGK1+Pld2HGqoNzo1K02j5KC+XrDlEbEjUFH7hr7xpiMjA/uZex3TNTqWckcyGykJVTd8rW3phKPbljfWyo4D58wOagAAAAAAorBADQAAAABAFGOP+NjoD145zo37IIL5dvg61PA1CXAx+7GP7irqpTz823Evf6DsONTQ3c1ucLz6W39V1AuvLTsNdfTOpx+PHYGK+p2NrxX1XKsT9F44f03ZcaiZpZFRCwtJK1ISqq7TGizvPXbLVUHvjnx/2XHgnNlBDQAAAABAFBaoAQAAAACIwgI1AAAAAABRjD2DOkmSB62hLGu9jeB4pml+F5PXabWLetfC1qB3//KxsuNQE7986ONFfdO2K4LeFZ0dZcehJprJYJ/CjWk4l3XhHz6j7DjUXOc5N8aOQEU90D1d1C/d+qig98nuobLjUDPLG6vB8Y/u9dsiTMfs0PfYx8yE3x++OuNex3Q8e8+jNu19/MiXzun/soMaAAAAAIAoLFADAAAAABDF2CM+FtuzRd1Ne1MJA6OaQ+NksjxiEGrj0sXtRf0znccFvdcsf7jkNNRFL+1v2ruutVRiEuoky7Oi/jsfCHfuKjULNJ/+/NgRqKhHLOwt6v/6yy8Oer/2+r8uOw41t+O1j3vok+A8XLu4p6j/wz+bDXpHfun6suNQE3/+8y/dtLfl1UZ8AAAAAABwEbBADQAAAABAFGOP+Fhqzxf1qZFfooVpyfLBXI88D2d89LO07DjUwOnu4P5232JyljNhOpIkvO5mGq5Dpu/vDJlZPhUjBnXm+wUlyM+cDo53pNkmZ8J05BsbsSNQUXNJu6ibT3pK0NvT/HjZcaiJ7BOfmtj/ZQc1AAAAAABRWKAGAAAAACCKsUd8bGS9oh7+1XeYpubQq+7ZyGvvidfemYK1Xreov5p4BY/yDT9vG41G42C2HikJdXK0GX62y770lUhJqKvuLb9T1PMv/ucRk1A1x/srRX3iZ98f9H6zvVTUry4tEXXW+xvPV6bjZP9MUS+/431B70i2vdww1MaB901uRJsd1AAAAAAARGGBGgAAAACAKCxQAwAAAAAQxdgzqNPM3GkuLHkjjx2BCupl/dgRqLnuyDV4Ju9tciZMzuFGNzg+9RfHinpp9GSYgjzzuY7pONMf/KbIbXdfF/TWO6fKjkPdWVZhSnp5WtRrx8OlvszaCRcBO6gBAAAAAIjCAjUAAAAAAFGMPeKj3WwVdTOxrk05hq+1ZOR9qHZr7MsXxrZ1dqGor2p0Iiahrmab7eB4Z3MuUhLq5LFZeJ3teNNzIiWhrjovemrsCFTU1pn5on76k+4Pet/4havLjkPNdZ77mNgRqKgtQ/e6S16xJ+jt/FUjA5mOy154lubXzu3/stIMAAAAAEAUFqgBAAAAAIhi7BkJW9qD1wWOJctTCQOj8nzwa7P5yC/P9tJ+2XGogVPrZ4r67sbGWc6E6eg0w0fz1sQ4I6Zv9CpLtu6IkoP6Sh7/rNgRqKhT/dWi/shtjw1672keKOofKC0RdZZcc33sCFRUOxkay3tDeJ3tOtdZCzCm/X/eeuiTxmQHNQAAAAAAUVigBgAAAAAgCgvUAAAAAABEMfZgy/W0V9Rpnk0lDIwavtaG51E3Go1GIyk5DLXTbbjXUb5+ngbHayPHMA3HmuEzNrv1s4ODf1ByGGop33/n4ODp8XJQPRtD32M/vxR+gTizvl52HGou+/znBgeer0xQsE53/FjQW018r2U6ltdmJ/Z/2UENAAAAAEAUFqgBAAAAAIhi7BEfq/2Nos6M+CCCJAlfyUvM+GAKhq+zff1TEZPA3/LEBWrhzHLsBFRUNjQm8EQSjs062TtTdhxqbv2vvlbUixFzUD1rWbeou5+8M+gdzbeXnIa6+FBnfmL/lx3UAAAAAABEYYEaAAAAAIAozmvER57nZzkTJmd4nMzfue5M+GAKmsng73YP9E5HTEKdzM50inp0fNFy3is7DjV0NOkHx73PH4mUhLrKvvq1hz4JzkOrOfhstydrBb0bFy4rOw4196lPXV7UL4uYg+pZSwcjPvb/9c6gd2/f91qm47Z8ZWL/lx3UAAAAAABEYYEaAAAAAIAoxh7x0U/Thz4J4CKXZoN73Up/PWIS6qSbDsZ4nOiGr0nd2+yMng4T95U0fPXz4x++tKi9gkwZ/u8vDZ6/L3hLxCBUTpoNRgauJOHIwGP9yb2aDOP4k/nBNej5yiQt99eK+v3Z9qB3tGd0G9NxMF2e2P9lBzUAAAAAAFFYoAYAAAAAIAoL1AAAAAAARDH2DGqIIWkkwweBZuLvK0xekgwutNMbqxGTUCdLnfmi3tKeD3qdpkc107c1CWedXz53JlIS6mrvdrOAmY6FmdmiftlGN+jd2t5WdhyAqRiet388CX9Drpv1y45DTdy7dnRi/5cVPgAAAAAAorBADQAAAABAFGO/N7zUmSvqE+teD6Acw+MWGnnYayYjMz9gAhbbg3vdDVv2Br3bjt1ddhxq4h/vekpR787DR/PdjfWy41BDX5dvCY4f+R9vjJSEurrm3d8VOwIVde387qJ+4jedDnq/e+UlZceh5n5o+7HYEaioy+Z3FPWPPP1g0Lvtk3tHT4eJ+Oatj9q090vLHzun/8sOagAAAAAAorBADQAAAABAFGOP+Ohl6UOfBBOW5/mmvewsPThfaT7068e9lYhJqJPh+9lykgW9I+la2XGoobtaveC4+94PF/XCq3+i7DjUUPq/3zM4eNqr4gWhck4PPUfX724FvebBA0W9VFoi6uzTh/YU9bXxYlBBa2m3qO/99M6gdyL1vZbp+FL/5MT+LzuoAQAAAACIwgI1AAAAAABRjD3iY/i197ONXYCyNJMkdgQqqJv2Y0eghn7tyN8U9aO3Xxn0nti5tOw41NBlWfiRsP2MR0RKQl01n/yU2BGoqMXmbFEvvfGFQe/Iv/8/Rb2rrEDUWuYrLCVotcI1u42st8mZcOGwgxoAAAAAgCgsUAMAAAAAEIUFagAAAAAAohh7BnUrGaxlJ2b/UpKzXWtJw3UIVEPvLLPPO/6WTAkW8vCZmuwyjZVyJZffEDsCFTX8PbZ53WODXpJ8uOw41Mzo99m2n/NiSoZ/o6vVyoLebLNddhw4Z771AgAAAAAQhQVqAAAAAACiGHvEB1xo0jx76JPgHKVZWtSdplsk5VjszBX1bBK+gjeb+FsyAJyvYCxgK3zG3nP/jqK+pqxA1FozN+OD6ZhJWkW9bc9a0Ouc8b2W6djS7Ezs//KtFwAAAACAKCxQAwAAAAAQhX3+XNDyoVeg8kb4OlTwuh5MSKvZeuiTYMIuX9hZ1Ftas0Hv/my97DjU0MFmGhz3P3Hr4OC15WahnrKPf2Bw8LRXxQtC5aykg+do+oH3Br1uw+c+pisfGenxwIxrjukYHoHaPRNeZ928X3YcOGd2UAMAAAAAEIUFagAAAAAAorBADQAAAABAFGPPoG6by8oFZnQmNUxCmg3msM43OxGTUCd3nz5c1J3t4aN5z+zesuNQQ+2R33VoLs5uciZMycJC7ARUVNoYzGVNv3Z/0EsaW8uOQ83t7qcPfRKch05z8B1i6epw5nT7iJ+fYzp2JXMT+7/soAYAAAAAIAoL1AAAAAAARDH2Pv+FmcGrnqeS1amEgVFJMnjlOM+N9GD6ljrzRX11Z2fQu6Oxr+w41EQ37RX16L1ua+KVPKZvWx7uWWg+6hGRklBXyU1Pih2BippN2kU984zHh7333V12HGrucXuOxY5ARQ2Pp5x73g1h7zMHyo5DTbyku/lYwP95jv+XHdQAAAAAAERhgRoAAAAAgCjGfm94a3vwy9pHklNTCQOjWsngbyhZngW9pJGMng4P22Jn8Cu0z0m2B70/KTkL9dFpDV4/nmt1znImTMf86BStLVuj5KC+km17Ykegomabg6+8zSd8XdDbs+1zZceh5nY/23dYpmNxaMRH87kvDXqXvOu3y45DTfz9V59lffgXzu3/soMaAAAAAIAoLFADAAAAABDF2CM+Dq2dKOpe1p9KGBiVjoz1GJYkXo9i8nrp4P72J+mhiEmok+Hn6nraDXrrZ7kPwqScTsIZH/l990ZKQl3l+z8/OHj8Szc/Ec7RetYr6uwLnwp6nzl9SVHfVFoi6uyzf7ilqF/w3yIGoXLOZIPvENnH/jzorebW8JiOd75v28T+LzuoAQAAAACIwgI1AAAAAABRWKAGAAAAACCKsWdQ753fUdQr3fWphIFRrWTwN5Q8D+djNs2gZsqe0t4dHP9V44uRklB1w/e3Xp6GvUY+ejpMXG/0kdrtPeh5MC35yWOxI1BRw79pkx85EvTubZedhrp79BOOPPRJcB6G106SS/cGvS3JidHTYSJef8X9m/Z+cvPWg7KDGgAAAACAKCxQAwAAAAAQxdgjPh4/d1lR33vm6FTCwKj5dqeoR0fLJA0jPpi8y+d3FvU/2ugGvXeVHYZaWmjNBseXNjqbnAmTc1U/fKY2X/jNkZJQV+2XfmfsCFTUjpmFou78kzcHvce//UfLjkPNLX37s2NHoKK2NueKuv1Nrw96u3/cvY7p2P2H/33z5uXPO6f/yw5qAAAAAACisEANAAAAAEAUY4/4ONBfLupu2p9KGBi11huMWMiGfoG70Wg08jwvOw41cNfyoaL+893XR0xCnexa2FbUzZHxRbdlp8qOQw19uh2ONHr5vx0MNZr/028tOw41tP7W7yvq9jv/OGISquZw73RRr//Um4Jet7Ft9HSYqtXf+XhRL7wuXg6q58jQmt36234w6O3LrJ0wHav/6rsm9n/ZQQ0AAAAAQBQWqAEAAAAAiMICNQAAAAAAUYw9g7qXp9PMAQ8pGZnLmiTJJmfC+btqaVdRP2k9O8uZMDm757YW9fWdnRGTUFdP63WC48V/8/pISairuR97Z+wIVNSl7cEzdu4nwuvsuff+07LjUHPzL31U7AhU1J6ZLUU992/eEfS2/N6PlB2Hmuh852s2b777zef0f9lBDQAAAABAFBaoAQAAAACIYuwRH8vpWlHneT6VMDAqzQcjFrI8HLcwOvIDJuErJ+8v6v9+6Z6ISaiTr50+XNS9LBypdfXsJWXHoYY+OrMWHL/wTe8u6utuf13Jaaij5e8ajJXZ+d6PRExC1Rzonijqk9/xhqD3G1+4uqh/oLRE1NkH37Zc1K9w0TFBw/e606//zqB3XzpXdhxq4n3f9oGJ/V92UAMAAAAAEIUFagAAAAAAohh7xMex7uBVlF7Wn0oYOJskCUd6GPHBtO3vHo8dgZp4267nFvVCOM2o8cfNUyWnoY6ua4Svfu5+5dgfEWEi5m5+XuwIVNRCc7ao27vDe9vTe2ujp8NUPf2GQ7EjUFEzzVZRty+bD3qdW1ujp8NEfLCzMbH/yw5qAAAAAACisEANAAAAAEAUY7+/udwdvP6U53nQGx29AJMyfK0Z6UEZhq+5e1aORExCndzWXC/qpVb4Ct5Kf3KvTcFmDie94HjjtsNFvVR2GGop/dStg4NvjxaDClpOB99j/+dHbgp6VzWy0dNhoppJuCfwD/dfUdTfXXYYKu3IxmAs4Lv/5Nqgd9fGl0pOQ118cPnLE/u/7KAGAAAAACAKC9QAAAAAAERhgRoAAAAAgCjGnkG90e9t2hudSQ2TkjfMoCaelaHZ+zBNH1u9p6i3zMwHvSw3H5Pp+1q6HBz/n08OZmTeXHYYaukzvz1X1C94Z8QgVM6xjcH97cB8GvR+rzWYt+9exzS0muGewHfn9xe1GdRM0rG1wb3uXTPhXOADZ46VHYeauG/56MT+LzuoAQAAAACIwgI1AAAAAABRjD3io5UM1rLTRnqWM2E6hsd9AFTJye6Zoh5+3jYajcZCa7bsONBoG99GyeZb/dgRqKiZZquoL89aQW8+aZcdh5rppeG9bWtzbpMz4eFpD93rrujsCHqHWifKjgPnzA5qAAAAAACisEANAAAAAEAUY4/4SPNsmjngIeVeNwYq6pLZLUXdTsLXjx/YOFV2HGroYPdkcPyB+a1F/fKSs1BPf9xeKOpnRMxB9Wz0e0X9peZG0DvVXSs7DjV3oGvUAtPRywajeA/2Tga99aH7IFyo7KAGAAAAACAKC9QAAAAAAEQx9oiP4V8/TodeHYCyJEkSOwLAVCy1Br/ovmNmIehtn1ksOw419Li5vcHx9y6cjBOE2voX1x6IHYGK2j47eI6+oRe+5v6uzq6y41AzzSTcE/isuSsjJaHq5mbaRX3T7KVB70jHyEAufHZQAwAAAAAQhQVqAAAAAACisEANAAAAAEAUY8+gXmzPFnU37Z3lTJic4ZldWZ5t2oNJ2bWwragfv+WqoPfhw3eUHYeauGf1SFHPL4WzCa+Z2Vp2HGrosXk4+/yK774sUhLqatuPvDJ2BCpqT2fw2e5xP3V90Hvljx8qOw41M/od9q2POhwpCVV3+cLOon7XlctB79Xp5WXHoSZ+bu+LNu39wKG/OKf/ywofAAAAAABRWKAGAAAAACCKsUd8QAzNJCnqPE+CXqfl8mXyhq+ry1pLEZNQJ8dWTxf1oc7JoHdJa77kNNTRapIHx9nBI5ucCdORff5zg4NvjpeD6kmHRizkp04FvTs7yejpMFXzzzRCi+mYbbaLevER4V7Uzv5W2XGoiaf3Vyf2f9lBDQAAAABAFBaoAQAAAACIYuwZCc3EWjblSxqbv3aX5fmmPThf3bRf1NuT9lnOhOloN8NX8HYks5GSUHWd1uAetysLP+c1b7i25DTUXXLDI2JHoKJmksFzNbn6mqB3sPn5suNQc809u2JHoKLmm52ibn/9U4Le1v/zpbLjUBNXXH5q8+bRc/u/rDoDAAAAABCFBWoAAAAAAKKwQA0AAAAAQBRJno83yLcze1VRp1m66Xlp7/6Hnwr+/2Y6VxT12S5V1x2T0mpfPtZ5rjkmqd25sqh3zC8FvR2dwfGdD3y6tExU3/D97mm7Hxn0XjEzuCbfvP/dpWWi2s72jP35S19U1N9932+UEYeauGTLjUX9vB03Bb071w8V9ReP/E1pmai2s93rrt66p6jvPnZrCWmoi7m5q4t66+xC0Du9sVrU6+v7S8tE9U1y/cQOagAAAAAAorBADQAAAABAFDPjnrjQni3qle7aVMLAqKSRDB9Aqa7aujs4vvf0A5GSUHXDr+HdsLg36O1qLZYdhxoaHunRaDQa3/f9C5ucCdPxul96auwIVNSNWwavH7/7DVuC3od/wTOW6UqS8EvsL87ctMmZ8PBcv3XwHeJ9O3cGvdecWh09HSbiBZc+btPeRw7fcU7/lx3UAAAAAABEYYEaAAAAAIAoxh7xsdHvFXWe51MJA6PyxuBaG73uRl+XgkkYvq7O9NYjJqFOWs3B34uPdE8FvUP5yZLTUEe/39sXHK++c/BL8G/9wbLTUEe/+Mb/W9Q/cG/EIFTO184cLur/8j+uCHp/0Rk8c/9+aYmok9HvsLfMdYv6G8sOQ6Ud3Thd1L+4fn3Qe6B7qOw41MSh7smJ/V92UAMAAAAAEIUFagAAAAAAohh7xMfW2cGvuR9bO32WM2FyZpqtok6zLOgNvxIPk7Jjbqmov23nk4LeLxz4y5LTUBdf/dbBOIXe4XC0zGtvXRo9HSbu5vY1wfH3//C2SEmoqzf+4lNjR6Ciblq6vKjf9F+eEvT+0U/8SdlxqLkf7BuXynQ8eunKov73//XZQW/pjdZOmI4/u2rLpr3rTpzb/+UqBQAAAAAgCgvUAAAAAABEYYEaAAAAAIAoxp5BvdxdK+o8NzeJcozOnR6WNJISk1AXJ9ZXivrdxz4TMQl18uMf3LFpbzk7WmIS6urLSTj7/Oj/PFTUV7+p7DTU0an/8L6inv/mH4gXhMo52lsu6tPv+KOg98PHthf1e8oKRK39enu2qJ8WMQfVc2DjeFHf/ob/E/T+IjlV1G8pKQ/18KJ95zho+izsoAYAAAAAIAoL1AAAAAAARDH2iI9Wc7CWnSRGK1CO4Wut6bqjBMMjjFa662c5Eybn0mzwOB79y/F9rcVyw1BL1zXmguNLXnlZpCTU1ZYXXRE7AhU112wX9fxN80Gvd9fm4wRhErbMLgTHX0pPbXImPDyne6tF/UuL4TeKfcsPlB2Hmji8dnJi/5cd1AAAAAAARGGBGgAAAACAKMYe8ZFmg9efhl+Bh2lqJYO/oeyYXwp6vbRfdhxqJs3S2BGoiX97+KNFfcO2vUHvMXNGLTB9vcboZztjtZi864bub3efOhT08vWNsuNQEyd6Z4r6o39wVdD7YvbVsuNQM1ct7gqOl1P3OqZjeWOtqN977Lag1/e9lilZneBYVDuoAQAAAACIwgI1AAAAAABRWKAGAAAAACCKsWdQm1lDDGk+mH1+ZmS2TfZ35mUCXJyG552f6q0GvYMzy2XHoYa+nJ8Jjr/63wfP3yf+eNlpqKq5VnvT3id+pVXUL/mZMtJQF6e6g/vb2+eOBr3jy56xTNdyfy043sh6kZJQdcNrJ6s9s865+NhBDQAAAABAFBaoAQAAAACIYuwRH9nQ6wJQljwfjPHopv2ISQDKsd7vBsdHe6cjJaFO7tw4EhzfMn9NUb+j7DBU1n1njm3a++jc4GvJS8oIQ20Mf4c4k4Wvva+nxi0wXSc3whFas2cZdQQPx/DayVy7E/Q2+u51XPjsoAYAAAAAIAoL1AAAAAAARDH2iI/59mxRt5JwXTs1/oMp2Tq7UNQv2HFT0EsaSdlxAKZuS3s+OL569pJISaiTJ89dHhz/0N7DkZJQZY/YcllRf3bjrqC3K7Nvhum4eml3Ud+9Et7brlrcVXYcamChM1fUT9t+fdA7bHQbU9JqDp6jaWaNjouPT4IAAAAAAERhgRoAAAAAgCjGHvFxzdCrUf+/1iOC3tuzuyeXCIbcuGXwyvEtL0+DXvPKPWXHoWaSJBwj026OfcuE87ajsxQcP3ZmR6Qk1Mnz++FomZ3v/BeRklBlz54dfK77bCMc8fFvT36yqL+3tETUwSe+ZWdRv/APO0Gvk/hsx+Q9a8cji/p7ejuDXi8xVobp2DE3+A5xcv1M0Bsd0wuT0m5t/hxdbM9t2nswrlIAAAAAAKKwQA0AAAAAQBQWqAEAAAAAiGLsoVszSauob0rCeTYLzdnJJYIhSSPZvLe0WGIS6ijP8+C4l/UjJaFO+nk4b3915BimYXV0y8KJQ1FyUG0HstVNe6u9jRKTUCcn/6ZX1KMzp7+hffno6fCwbR1aH+k0wu8TW3LfJ5iOS+e2F/Vydy1eEGqll25+TzvTWD+n/8sOagAAAAAAorBADQAAAABAFEk++g77Jhbmry3qyxZ3BL2DZ04U9eraPRMJBo1Go7Fl8fqifvTWK4PeXLNd1B898KHSMlFtrTFf9Ux79085CXUyfN3Nt8OxWYvtuaI+fOqLpWWi+oavu2u3XRr0nrpwdVH/7v4/LC0T1bZzyyOL+tT6mU3P84xlkp562XOL+osn7wt6zWQwTnBl9e7SMlFt25ceUdS75rYGvbnW4DvsHYc/WVomqq8ze1VRp9nmIwI9Y5mkSa6f2EENAAAAAEAUFqgBAAAAAIhi5qFP+VvXbNld1P8puSHo/fTOI5NLBENu2npFUf+vRyVBb+EpS2XHoeY6Q6/kwbRcvrgzOH7i/JWbnAmT84qlm4Ljf/uS45GSUGXP2j4Y8fH+Q7fGC0Kt/PaOxaL+toZnKtP35O3XFfWrknCE1uaDF+DhuXbrnqK+d/loxCTwt7bPnduanR3UAAAAAABEYYEaAAAAAIAoxh7xsdCaLepH7z4W9FrHrHMzHdtb80W95dseE/SSJz2n7DjUXKvpXsf0DT9vG41G46rm/CZnwuRcn4YfCWe++UWRklBlVzQXYkeghq5856uKett3vjfo+WTHNOwautc9q78S9FrNvOw41MRSa66oZ1qtiEngb3VaYy85NxoNz2QAAAAAACKxQA0AAAAAQBQWqAEAAAAAiGLsgSBJIynq+R3doJceM9+G6ZhJBtdWct2NQa+58/Ky41Bz5zpDCc7H8H2v0Wg05nJ/S2b6tqXhTMzmNY+OlIQqu7zRiR2BGmre8JSibifvixeE2mgng89uSwsbQa/b9X2C6UiS5KFPghItzsw99ElDfOsFAAAAACAKC9QAAAAAAEQx9vsl/Twt6lOHFoNeq9GfXCIYspoNxslkn/zroJefPDo4eORzyopEja31ug99EjxMG1kvOD7W6G1yJkzOsVb4Wmj21dsGB095RclpqKqvNdY37Xk1mWnp/9a7ino5DcctNF13TMF6PlgfeWB5Iei1knz0dJiIbja47vLcdUZ8JzZWzul8O6gBAAAAAIjCAjUAAAAAAFGMPeJjoTlb1HufHY70WPyQX+RmOpaGrrvmM78u6DUf8bSy41Bzvcw4I6ZvttkOji9ptDc5EyZnbz98FXTm614ZJwiVttBobdrzOjLT0voH317Ui//154KeER9Mw55krqivvfx40Jvd6vsE09FsuJ9xYTnT23y024OxgxoAAAAAgCgsUAMAAAAAEIUFagAAAAAAokjyMQe+bVm8vqhv3Hp50Ltr+VBRn1z56oSiQaNxyZYbi/obdz426M0MzVj6tX3vLS0T1dZqX/7QJzUajbR3/5STUCfD191SZz7o7ZzbUtR3H7u1rEjUwPB1d9OOK4Pec+avLupf3vd7pWWi2i7fMfgsd3jlxKbnecYySS+44iVFfdupezY9z/dYJmXPtkcV9bULe4JeMjT3/JP3f6S0TFTf0sJ1Rb3R7216Xq97XxlxqIlJrp/YQQ0AAAAAQBQWqAEAAAAAiGJm3BMvW9hR1P+6cXXQe8vi5q8PwMNxzeLuov4Pl54Keks3jp4NkzX8Ch6U5dKF7cHxE+evfPATYYJePPRaaKPRaPzUUw5tciacv5sWB6+Bnm3EB0zSf8gHo7K+e+nSiEmoi0cvDT67vaaxN+it2SLIlOyYWyrqo6unIyahTtqtsZeVH5LbIwAAAAAAUVigBgAAAAAgirH3Ys8220X9mPlw1MLi+tzkEsGQbTMLRb3j264Nes2bHlNyGupmbqYTHHfTfqQk1MlCazY4vrzpGcv03ZC2g+O5m58XKQlVdmVry0OfBBN26Z7lol56IHymthL7tZi8S1rzRf3UPBy1MNPKyo5DTWxvLxb1ydaZiEmok+HRMqNWeuvn9H95IgMAAAAAEIUFagAAAAAAohh7xMf1s7uK+qZPvjPo3fi0fzW5RDBkPhm8ctx8evi68cyjnlN2HGrgiZdcV9QPdMNX8k6sr5QdhxqaSVrB8WLe2uRMmJytI28cJ9cZo8Xk7U46m/Ym+SvwMOzSX/7uop7/ll8OekkjKTsONfAd3cE4o5u+O+wluy8pOQ11cXVnZ1EfaB6LmIQ6uWFx76a9u84cOqf/yw5qAAAAAACisEANAAAAAEAUFqgBAAAAAIhi7GFvq1mvqPu3fiDoncw2JpcIhqSNwVDMfP+Xw94lVxV1e/cNpWWi2l7dubao/9Pqp4NeK/E3Paavn6fB8XqSbXImTM7q6O3t2MEoOai2MyP3t2Fp5l7HdCRD883T3HXG9J1oDR6q/bvDGazNB06WnIa6WMm6Re2ZSlmO9zb/na5O89x+X8RqCwAAAAAAUVigBgAAAAAgirH3Wx/snSzq9732L4Pe3dnRiQWCYcvpYHxM90//Ouh1lpcHB495YVmRqLgf+r9vKeofvfx5QW/b3GLZcaihmaQVHM/l/pbM9F3dDUcvzDzrlXGCUGkrjf6mva2zCyUmoU6SLbuK2rg2yvDhmdWivv5PdwS9maHRbc/9+dIiUQNfWR2MZ+tlm4/Ugkm688R9m/au27b3nP4vT2gAAAAAAKKwQA0AAAAAQBRJnuf5OCd2Zq8q6plm+PpxN+0Vdb97YELRoNFYmL+2qK/ZsjvstWaL+v8e/FhZkai4V1/9iqL+/YN/E/SSJClq9zomqdW+vKiXOvNBb/vsYLTMvuO3l5aJ6hu+7m7acWXQe8781UX9y/t+r7RMVNv2pUcU9fLGatDzjGVannrZc4v6CyfuDXrD193q2j1lRaLi5uYGz9CrR77DDo+Z+eKR8LsGPBztzuCzXN7YfJnPM5ZJGv4+MardGkyVXl/f/5D/lx3UAAAAAABEYYEaAAAAAIAoZh76lL9reKQHTFNz6LW7bhb+8vvWGb/2zuQ9IlnctNduntctE87J3EwnON41uzVSEurkis6O4Pj5vdlNzoTzd7bJgqPjjWBSnjt7RVHf0z4SMQl10UsH31t3tJeC3p7W0ujpMBHDI4vOMuEDSrPYnjun8+2gBgAAAAAgCgvUAAAAAABEYYEaAAAAAIAoxh6oOjvTLup+mk4lDIwanlmzqxPOYX1R+/Ky41ADb/629aL+2XdEDEJtXbWwKzh+aeeqSEmok9dm4XX3rX/86khJqLLd89uKeqW7FvTazVbZcaiJn/7OwbX1V794acQk1MXwLODndS4Lek/bcK9jOmaGnqNn+80HuFDZQQ0AAAAAQBQWqAEAAAAAiGLsER9plhV1L+tPJQyMWu93i/r+9eNB72+as2XHoQ5am//dbtfC1k17MCmjz9jlxFgtpu9UKwmOs31fHBw85RUlp6Gqrpq7pKjvPnUo6K2nvbLjUBPpvUeLup97pjJ9w+MV7kyXg95cZ9vo6TBxeSMc8ZE0kk3OhOl59NYrz+l8O6gBAAAAAIjCAjUAAAAAAFGMPeKj1Rxay/ZmFCUZ/gXk7e3FoLfViA+moPMvfmJw8PZvCnrb2gslp6GO2s3w0bwl92vvTN+ufvgqaOtxz4uUhCo722e39V530x48HM2924u6lawGPa+9M217m/PB8bV9ewSZjubQ2kmeh/e24XUVKMs1M+c20sjdEQAAAACAKCxQAwAAAAAQhQVqAAAAAACiGHsGdZpl08wBD6qfDQaeH+suB72FlhnUTF7vf7xt014365eYhLrq5+EPPZxJ/PAD03eqFc4mzO7/yuDgkc8pOQ1VdTrb2LS32JkrMQl18tVfWy/qU73Vs5wJk3cgC6+5z7X9tgjTkeX5Q58EJVrJe+d0vh3UAAAAAABEYYEaAAAAAIAoxh7xMSz36gAlGX5NpdkIXz9eT7tlx6EGjv3R4U17ezvbywsCUKJeMvIPa8sPeh48HDub85v21vs+1zEdf5VtLerR7w+txH4tpitthKNSW43RBy5MxtLQqKzljbWg12md19IfPCxPT7ad0/meyAAAAAAARGGBGgAAAACAKMbe599qDtayk8RrKZRvpbceHG9tL0RKQpVd8q1XDQ5u/2LQO95fKTkNdTSThL/uvpj7tXemb1c/HN/WvOEpkZJQZVuS9qa9zAhBpuSZ+eDz28LMbNBLjFtgyrp5Ghyvj4z8gElZO8uorG7aLzEJ/K0vNlbP6Xw7qAEAAAAAiMICNQAAAAAAUYw94qOfDl5Nyb2CR0myfPAKVD8LX486018fPR0etg//3ObX1b0rR0tMQl2d7oevQn2lfSZSEurkU7PhM/Zl/+4ninr21z5Udhwqal96etNe0whBpuS323NFvbYWvgLfNOKDKTvaWw6O70qMbmM6htfs0jwcJeNeRwwfW77rnM63gxoAAAAAgCgsUAMAAAAAEIUFagAAAAAAohh7BjXEMNMczOhaHJof12g0Gpd0tpQdhxp46o2HBgefCHutpr/pMX3zzU5wvLc5t8mZMDlXp+FHwvY3PjtSEqpsKek89EkwYdf1B/e32WY76LUSn+2Yrp0zS8HxU5o7IiWh6ubbg2fsWi+ct99pWfqjfO3muV13nsgAAAAAAERhgRoAAAAAgCjs8+eC1s/Sol7urgW9Q82TJaehDk4f2XycwuirUjANaSMLjtfzdJMzYXJWkzw4zg8fjpSEKlvP+5v22l4/Zkq+PNMr6v7IMzVv5KOnw0Qd6y8HxwdntkVKQtWt9wf3utF7Wzfd/PkL03JiY+WczreDGgAAAACAKCxQAwAAAAAQxdjv0vUyrwRQvjQbvOqejbx+vNrbKDsONfDW9YVNe62mv+kxfQ+snw6OP9M8FCkJdfLh7GhwfNV/3FXU3/GDZaehqvZtHN20t+ZzHVNyMFsv6pMbZ4JeM0nKjkPNjH6u+1R+b6QkVN3weNQ8N76I+DrnOL7NagsAAAAAAFFYoAYAAAAAIAoL1AAAAAAARDH2QJB2c3BqN+0FvcTsLqZkptkq6tHrbFtn81nBcL6eng2uq98c6e2YWyo3DLV0yeyW4PiJnUsjJaFOntvaFRx/y3f5bMfk7Wh7jlK+H0kH97PPzG0NemZQM2175rYFx4+a9bmO6Ri+n2Vn6UFZbljce07n20ENAAAAAEAUFqgBAAAAAIhi7BEfZ5Pn+ST+GzirVjP8e0oz8fcVJm93f/P72e7O1k17MCmtkXvbXNLa5EyYnEuy8NXP5HGPi5SEKmudZW+M7xNMy/zcYDxls+s1d8o13+wExx3fYZmSpDG4vzWTzXtQlie0dz30SUPcHQEAAAAAiMICNQAAAAAAUYw94iPxq59E5tVPynBVvr5pL81Hfw8ZJu/vjPjwt2RKsC0Nn7HNG54UJwiV1j7LyCLfNZiW3dcvF/XMbeE12HTdMWVZI3y+7k5mIyWh6oafo9ZOuBBsnOP6iW+9AAAAAABEYYEaAAAAAIAoxh7x4RUBYhh+TWV+JvwF5MUZr0cxeTuX1gYHx8Jeuzn2LRPO2+jzttfw/GX6TrXC19zz+78yOHjCN5Wchqra2tz8s1vSMGqB6Zh/3Painvvcsc1PhCnY0poLjpfzfqQkVN3wyKLM1wcuAFuSc1s/sYMaAAAAAIAoLFADAAAAABCFBWoAAAAAAKIYeyDI8CxggKpa39j8ttjLzIxj+kaft21zWSnBQjbyDzv2RslBtWVnmanvuwbTkp/Z2LTnumPa0jx8wM74XEcE7nXEcK4z9+2gBgAAAAAgCgvUAAAAAABEMfaID4gty8PXQlv+vsIUnOzObtqba3VKTAJ/Kz3LK/EwKaujj9QTh6LkoNpWsu6mvTRLS0xCnTQvWSrqJDkTMQl1lIyM9DjbqCOAKtnVaJ/T+Vb4AAAAAACIwgI1AAAAAABRjD3io++1OyLIh8Z6bKS9oLeSrpcdhxp4/9zmt8XDGyfLC0JtraYbwfH9mXsd03dnK3zGrv233yvq+W96U9lxqKjj/ZXYEaihX//1+aI+uH4iYhLqaPS+d6i5+ThBeDjSLCvq0fGoWW49j/L9TXrsnM63gxoAAAAAgCgsUAMAAAAAEMXYIz4ghiQZ/OrxTLMV9OZbnbLjUAOv7A7GK/zHkV4365cbhlqab4b3tr3NuUhJqJMb0vBXtude+XWRklBls83Nf829NfI5Dybl8f21oj7bNQjTsNQKP8dd3VqMlIQ6aw6tq0BZznUsrx3UAAAAAABEYYEaAAAAAIAoLFADAAAAABDF2DOoW8lgLTttpFMJA6OGZyUttcP5XXNNM6iZvIVOb9PeVfO7SkwCUJ7RD4TJDvc7Jq/Z2HwGZpr5fsF07FwazKBunQz3Z5nLyrQtNWeD47PdB+Hh6LQGn+Y20vA7bTOxN5XyvbBzxTmd7yoFAAAAACAKC9QAAAAAAEQx9ogPiK018lqK16OYhiPrC5v2Wv6mRwmSkdeNXXWUoZ2P/MPClig5qLZ549mIYHVtcN3ljfBml+ajNz+YLGNkuBC4Dolh9hzX7HzvBQAAAAAgCgvUAAAAAABEYcQHF421tBscZw2v5DF5N115dHBwPOydTlfLDUMtpXkWHK+PHMM0nB7ZspDf8+U4Qai00fEKw9otX0uYjrV0cG1lRnpQstHPdUu5PYJMRy9LN+259xHDXY21czrf3REAAAAAgCgsUAMAAAAAEIUFagAAAAAAojDsjQva8KykbtoPesv9c5tnA+P4/H27h46+GPRa/qZHCUZnFfYaZlAzfb3R2cCrZu4zeRt5f9NeL928Bw/Hibxd1JnfdaBkp9P14HhlxjXIdDSTpKgzI6e5AMw1Wud0vtUWAAAAAACisEANAAAAAEAURnxwQRt+TWW4bjQajaWZubLjUAOXdjZ/rf2B7ukSk1BX/TwNjpfP8ko8TMrRZnjd9f7mzkhJqLL1tLtpr93ytYTp+LqnHCjq7Z9ZjJiEOspGRmiloyO1YELSbPPxMcYbEcPnukfO6Xw7qAEAAAAAiMICNQAAAAAAUXiXjgtaEvwSbfg61Exybr8ICuPY6G9+WzRWhjJsnVkIjh+ZeB2Z6XtyN3ymzn7/90ZKQpXtbC9t2uulxhkxHQt/73GD+rYvBL2kkYyeDhO1pRV+f9iVW4JhOlrNwf7T0XEfSeJeR/kua289p/PtoAYAAAAAIAoL1AAAAAAARDH2+yVeCSCG5tBrd5fMbQl6C81O2XGogUe/+OTg4HeixaDGNrJecHys0dvkTJicY63wc1721dsGB095RclpqKpenm7am53xuY4paQ1GGI2O9JhtGrfAdLVGrrml3LoK0zE6EnVY03oeFwE7qAEAAAAAiMICNQAAAAAAUVigBgAAAAAgCkO3uKBljcEcpbW0G/RWs+7o6fCwbRzYfD5mK/E3PaZv9Dqbc91RgvmRsYXJpVfHCUKldbP+5r3UvH2mI9t3oKjTRhb0fJ9g2tJG+IDtGgXMlLSag+8MaRbe6zotS3+Ur9NoPfRJQ3zrBQAAAAAgCgvUAAAAAABEYZ8/F7ThV1M2+uGrn6vpRtlxqIF9d+7YtLdjZqnEJPC3Wg3vgjJ97ZERH425xSg5qLatrblNe3k+ehHCZKzffrSoz/j+QGQHE+OMmL7RZ+p633VH+R6dnNv3CTuoAQAAAACIwgI1AAAAAABRGPHBRSNvePWT6etmm//S7LpfeieC0V9/h1KkXgVl8tzPiGH92OArb5avB71mYowW07WR9YPjvnFGTIlRWVxojib9hz5piB3UAAAAAABEYYEaAAAAAIAoLFADAAAAABDF2DOoW83BWnZiVhclOdtcuHay+axgOF+P+/rjg4PfD3un+qvlhqGWZkbubYu5ex3Tt5CFcwubl90YKQlV1mqc5XNdy0/jMB0nj84Xdd5YDnqZka1M2ej32b2NdqQk1Jl5+8TwtXT5oU8aYgc1AAAAAABRWKAGAAAAACCKsd+l66dpUee5d6EoR5plRd1t9IPeid6ZsuNQA7/811cMHd0Z9LpZeA3CNJzorQTHn2+fjpSEOvnobC84/obv/6Ginv3Dj5Qdh4o60D2xaW/4Mx9M0gd6O4r6ZPfuiEmoowe64ee421qLkZJQdcPP0TQPn6m5eUZEsJJtnNP5dlADAAAAABCFBWoAAAAAAKIYe8RHMvSrn4lfAKUknaFfdN81vzXoXT23q+w41MBLk1NF/W9GejvbS+WGoZYWWrPB8RXNhUhJqJPHpJ3geP51L4uUhCrbPrP5q+1ZbsQH0/HP/km3qN/9Py6JmIQ62tvZHhxf3zTig+mYnWkX9UY/HN3WatqbSvlumtl5Tue7SgEAAAAAiMICNQAAAAAAUYw94iPP8wetYZqyxuBaayXh31My1yFTsPuKlcHB/WGvm/fLDUMtrabhrx0fyFYjJaFO7myFr4L2/uADg4NXjQ48gvMz/LkOytL78pGi3sh6ZzkTJm8lWw+OV/M0UhKqrpsOvqvmI8/bfua6o3yn83N75tpBDQAAAABAFBaoAQAAAACIwgI1AAAAAABRjD2DOs2zaeaAB9UbmqO03FsLet3MPGAm7w/2XzF0dGfQu//M8XLDUEuj97p7No5FSkKd3NEP72+//6HLivr1ZYehsg6sbX4/u27b3hKTUCdv++TgfnZo7TMRk1BHx7rLwfEdzU6kJFTd8Jzp0d+NS5Kk7DjQ+NTK3ed0vh3UAAAAAABEYYEaAAAAAIAoxh7xMdNsFXU69OoATNNsq13Ue2a3Bb097a1lx6EGvuObjhT1994S9q5a2lVyGuro0rntwfGz5q6ME4RauTm5NDj+f/77MyMlocquXdhT1AdXwrEySzNzZcehJn70ZSeK+v3vuyRiEuro+vnw+fq8lu8TTEenNVjeS7NwRK8RH8Tw3C03nNP5dlADAAAAABCFBWoAAAAAAKIYe8RHqzlYy/Z6AGVJ88GrKQfXTwS94bEzMCmH/3rz+1vL3/QoQSsJr7O5xHXH9LXzkX+Y3xIlB9W22Oxs2juycarEJNTJmTvWi3oj60VMQh3183A86qkk2+RMmJzRNbumNTwimGuc25qdb70AAAAAAERhgRoAAAAAgCjGHvEx+iugULZu2g+Ps/4mZ8L563U3fw0lbbgPUr5Wwyt5TN/fGfExtxglB9XWHXnVfdg187tLTEKd/N+v7C3q1f4Xgp7RlUxbb+S+lzZGH7gA1bTe2Pxz34OxgxoAAAAAgCgsUAMAAAAAEIUFagAAAAAAohh7BjXEkOWDGV15Hs7rys3vYgpOnZnbtHd47WR5Qait0fn6x/NepCTUyddmwhn76e+9e3DwzG8tOQ1VlZxlpv5lM1tKTEKdnGgN/b7IyCN1Jtn8t0dgEtaybnB8b7YaKQlVd7a1E7+kRAyfXTtwTufbQQ0AAAAAQBQWqAEAAAAAiCLJR/f+AwAAAABACeygBgAAAAAgCgvUAAAAAABEYYEaAAAAAIAoLFADAAAAABCFBWoAAAAAAKKwQA0AAAAAQBQWqAEAAAAAiMICNQAAAAAAUVigBgAAAAAgiv8PGS1Hb9p3ZJUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 2000x1000 with 10 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 10, figsize=(20,10))\n",
    "for i in range(10):\n",
    "    colum_start = 3 *i\n",
    "    colum_end = 3 *i +3\n",
    " \n",
    "    ax[i].imshow(spectral_data.iloc[:,colum_start:colum_end]) \n",
    "    \n",
    "# Remove axis labels for clarity\n",
    "for ax in ax:\n",
    "    ax.axis('off')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>BR_Anhydrite_R</th>\n",
       "      <th>BR_Anhydrite_G</th>\n",
       "      <th>BR_Anhydrite_B</th>\n",
       "      <th>BR_Qtz_Fsp_Cal_R</th>\n",
       "      <th>BR_Qtz_Fsp_Cal_G</th>\n",
       "      <th>BR_Qtz_Fsp_Cal_B</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.007843</td>\n",
       "      <td>0.023529</td>\n",
       "      <td>0.015686</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.949020</td>\n",
       "      <td>0.866667</td>\n",
       "      <td>0.203922</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.909804</td>\n",
       "      <td>0.831373</td>\n",
       "      <td>0.196078</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.035294</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.992157</td>\n",
       "      <td>0.905882</td>\n",
       "      <td>0.215686</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.062745</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.694118</td>\n",
       "      <td>0.631373</td>\n",
       "      <td>0.149020</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.125490</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>230</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.054902</td>\n",
       "      <td>0.121569</td>\n",
       "      <td>0.070588</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.094118</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>231</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.117647</td>\n",
       "      <td>0.262745</td>\n",
       "      <td>0.156863</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.125490</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>232</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.125490</td>\n",
       "      <td>0.282353</td>\n",
       "      <td>0.184314</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.066667</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>233</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.125490</td>\n",
       "      <td>0.290196</td>\n",
       "      <td>0.168627</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.149020</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>234</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.109804</td>\n",
       "      <td>0.247059</td>\n",
       "      <td>0.172549</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.211765</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>235 rows × 6 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     BR_Anhydrite_R  BR_Anhydrite_G  BR_Anhydrite_B  BR_Qtz_Fsp_Cal_R  \\\n",
       "0          0.000000        0.007843        0.023529          0.015686   \n",
       "1          0.949020        0.866667        0.203922          0.000000   \n",
       "2          0.909804        0.831373        0.196078          0.000000   \n",
       "3          0.992157        0.905882        0.215686          0.000000   \n",
       "4          0.694118        0.631373        0.149020          0.000000   \n",
       "..              ...             ...             ...               ...   \n",
       "230        0.000000        0.054902        0.121569          0.070588   \n",
       "231        0.000000        0.117647        0.262745          0.156863   \n",
       "232        0.000000        0.125490        0.282353          0.184314   \n",
       "233        0.000000        0.125490        0.290196          0.168627   \n",
       "234        0.000000        0.109804        0.247059          0.172549   \n",
       "\n",
       "     BR_Qtz_Fsp_Cal_G  BR_Qtz_Fsp_Cal_B  \n",
       "0                 0.0          0.000000  \n",
       "1                 0.0          0.000000  \n",
       "2                 0.0          0.035294  \n",
       "3                 0.0          0.062745  \n",
       "4                 0.0          0.125490  \n",
       "..                ...               ...  \n",
       "230               0.0          0.094118  \n",
       "231               0.0          0.125490  \n",
       "232               0.0          0.066667  \n",
       "233               0.0          0.149020  \n",
       "234               0.0          0.211765  \n",
       "\n",
       "[235 rows x 6 columns]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "column_list =['BR_Anhydrite_R',\n",
    " 'BR_Anhydrite_G',\n",
    " 'BR_Anhydrite_B',\n",
    " 'BR_Qtz_Fsp_Cal_R',\n",
    " 'BR_Qtz_Fsp_Cal_G',\n",
    " 'BR_Qtz_Fsp_Cal_B'\n",
    " ]\n",
    "hsi_data = dataframe_KSL[column_list]\n",
    "hsi_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Since it is difficult to classify a hyperspectral data in general. We can apply different classical clustering methods to have some starting guess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#gm = KMeans(n_clusters=3, random_state=42).fit(X)\n",
    "#gm = GaussianMixture(n_components=3, random_state=0).fit(X)\n",
    "gm = BayesianGaussianMixture(n_components=3,covariance_type='full', random_state=0).fit(hsi_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[7.71968839e-01, 7.05429203e-01, 1.67939396e-01, 4.09525013e-04,\n",
       "         0.00000000e+00, 2.91311307e-02],\n",
       "        [9.40121984e-02, 1.71913397e-01, 2.09966969e-01, 3.27839958e-03,\n",
       "         0.00000000e+00, 7.28620311e-01],\n",
       "        [8.19402576e-03, 1.15901741e-01, 2.47505076e-01, 1.65117002e-01,\n",
       "         0.00000000e+00, 4.19518439e-02]]),\n",
       " array([[[ 4.33002341e-02,  3.92883055e-02,  8.77528201e-03,\n",
       "          -3.05456806e-04,  0.00000000e+00,  5.32015149e-05],\n",
       "         [ 3.92883055e-02,  3.56846591e-02,  8.03772688e-03,\n",
       "          -2.44771188e-04,  0.00000000e+00,  1.57119257e-04],\n",
       "         [ 8.77528201e-03,  8.03772688e-03,  1.95029488e-03,\n",
       "           1.20617224e-05,  0.00000000e+00,  2.54546847e-04],\n",
       "         [-3.05456806e-04, -2.44771188e-04,  1.20617224e-05,\n",
       "           5.65984466e-05,  0.00000000e+00, -3.76308121e-07],\n",
       "         [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
       "           0.00000000e+00,  9.57443183e-07,  0.00000000e+00],\n",
       "         [ 5.32015149e-05,  1.57119257e-04,  2.54546847e-04,\n",
       "          -3.76308121e-07,  0.00000000e+00,  2.63701229e-03]],\n",
       " \n",
       "        [[ 2.21480994e-02,  1.67385445e-02, -1.74426993e-03,\n",
       "          -2.57068059e-04,  0.00000000e+00, -5.40908620e-03],\n",
       "         [ 1.67385445e-02,  1.32208920e-02, -1.86351496e-04,\n",
       "          -2.04969136e-04,  0.00000000e+00, -5.79384057e-04],\n",
       "         [-1.74426993e-03, -1.86351496e-04,  2.71290581e-03,\n",
       "           2.33609363e-05,  0.00000000e+00,  7.99424700e-03],\n",
       "         [-2.57068059e-04, -2.04969136e-04,  2.33609363e-05,\n",
       "           3.29858306e-04,  0.00000000e+00, -6.90047843e-04],\n",
       "         [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
       "           0.00000000e+00,  8.80783534e-07,  0.00000000e+00],\n",
       "         [-5.40908620e-03, -5.79384057e-04,  7.99424700e-03,\n",
       "          -6.90047843e-04,  0.00000000e+00,  3.43544736e-02]],\n",
       " \n",
       "        [[ 6.02135857e-03,  4.69543985e-03, -5.13864917e-04,\n",
       "          -1.24450858e-03,  0.00000000e+00,  1.01950816e-04],\n",
       "         [ 4.69543985e-03,  4.40669430e-03,  1.28831391e-03,\n",
       "           2.73946461e-04,  0.00000000e+00, -6.11968745e-04],\n",
       "         [-5.13864917e-04,  1.28831391e-03,  3.89066851e-03,\n",
       "           2.93462709e-03,  0.00000000e+00, -1.59436184e-03],\n",
       "         [-1.24450858e-03,  2.73946461e-04,  2.93462709e-03,\n",
       "           3.86334631e-03,  0.00000000e+00, -1.44604084e-03],\n",
       "         [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
       "           0.00000000e+00,  9.02729088e-07,  0.00000000e+00],\n",
       "         [ 1.01950816e-04, -6.11968745e-04, -1.59436184e-03,\n",
       "          -1.44604084e-03,  0.00000000e+00,  8.76907498e-03]]]))"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gm.means_ , gm.covariances_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "loc_mean = torch.tensor(gm.means_)\n",
    "loc_cov  = torch.tensor(gm.covariances_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1\n",
      " 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2\n",
      " 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2\n",
      " 2 2 2 2 2 2 2 2 2 2 2 2 2]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "tensor([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "        1., 1., 1., 1., 1., 1., 1., 1., 1., 2., 2., 2., 2., 2., 2., 2., 2., 2.,\n",
       "        2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2.,\n",
       "        2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 3.,\n",
       "        3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3.,\n",
       "        3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3.,\n",
       "        3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3.,\n",
       "        3.])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_gmm_label = gm.predict(hsi_data)\n",
    "print(y_gmm_label)\n",
    "\n",
    "y_gmm_label_arranged = torch.Tensor(y_gmm_label+1)\n",
    "y_gmm_label_arranged "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# We use the normalised hsi data y_obs_label information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "y_obs_label = hsi_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/8AAAK6CAYAAABv4Ow9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABcjklEQVR4nO3dd5RU5eH/8c+ysBaqSJOICkRRbNhiAzWCBQVbJEENmtjQKLEXFHtFf1b0i9hijZ0kogRji4oajcaS+EVFUAFBUFFAFBeW/f3hYb7ZYEEFlr15vc7JOc6dZ2aei5PB99x7nymrrq6uDgAAAFBY9Wp7AgAAAMCSJf4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAAUn/gH4r3HyySdn++23r7GtU6dOGTJkSC3NaPHo1KlTzj777NqexmL33HPPpVOnTnnuuee+dWy/fv3Sr1+/JTKPZfE9sizOCYBlm/gHYJk1YcKEnH766enevXvWX3/9bLzxxunbt29uvvnmzJkzZ6nNY9KkSenUqVM6deqUhx56aKH7hwwZkk6dOmX69OlLbU58P9tvv33p3+Xaa6+dTTfdNL17985pp52WV155pbanV8MTTzwh8AFYbOrX9gQA4Kv89a9/zVFHHZWKiorsvvvuWWuttTJ37ty8+OKLufjii/PWW2/lnHPOWerzuvrqq7PjjjumrKxsqb/2f5vNNtssr776aho0aLBYn3edddbJr3/96yTJ7NmzM378+IwaNSp33313fvWrX2XgwIE1xr/66qspLy9frHNYFE888URuv/32DBgwYKH7amtOANRd4h+AZc7EiRNzzDHHpG3btrn55pvTqlWr0n377bdf3n333fz1r39d6vNaZ511MmbMmDz88MPZcccdl/rr16YvvvgiDRo0SL16S++kwXr16mW55ZZb7M/bunXr7L777jW2HX/88TnuuONy0003ZfXVV8++++5bum9R5vDZZ59lxRVXXOxz/TpL4s8FgGJz2j8Ay5zrr78+n332Wc4777wa4b/A6quvngMOOKDGtj/96U/Za6+9ssEGG+QnP/lJjjnmmEyZMmWxzmuXXXbJGmuskauvvjrV1dXfOHb77bfPySefvND2/7w2fcF17SNHjsxVV12Vbt26ZaONNspvf/vbzJo1K5WVlTnvvPOy5ZZbZqONNsrAgQNTWVn5la95//33Z6eddsr666+fvfbaK3//+98XGjN16tQMHDgwW221VdZbb73suuuuuffee2uMWTCnBx98MJdddlm6deuWDTfcMJ9++mnmzp2bcePGZdq0ad/653XyySdno402yuTJk9O/f/9stNFG6datW26//fYkyRtvvJH9998/Xbp0yU9/+tOMGDHiK+fxn9f833XXXenRo0c22GCD7L333nnhhRe+dS7fZvnll89FF12UZs2a5Zprrqnx7/c/r69fcJnHW2+9leOOOy6bbbZZjS8LFvW9+Morr+SQQw7JZpttli5duqR37965+eabS392C/6cFlym0KlTp6+dU5L87//+bw4++OBsvPHG2WijjXLAAQfk5ZdfrjFm+PDh6dSpU1588cVccMEF2WKLLdKlS5ccccQRC122MmvWrIwbNy6zZs36xj+7k046KZtvvnnmzp270H0HHnhgdtppp298PABLh/gHYJnz+OOPp127dtl4440XafzQoUNz0kknZfXVV8/JJ5+c/fffP88++2z222+/zJw5c7HNq7y8PIcffnhef/31PPzww4vteZPk2muvzejRo3PooYfmZz/7Wf7yl7/kjDPOyMCBA/POO+/kyCOPzA477JDhw4fn2muvXejxf//733P++ednt912y29/+9t88sknOfjgg/Pmm2+Wxnz44Yf5+c9/XvqzOfXUU7Paaqvl1FNPzU033bTQc/7P//xPnnjiiRx00EE59thj06BBg0ydOjW77LJLLr300kXar6qqqhxyyCFp06ZNjj/++PzoRz/K2WefneHDh+fggw/Oeuutl+OPPz4NGzbMSSedlIkTJ37j891zzz05/fTT06JFi5xwwgnZeOONc/jhhy+WL3oaNmyYHj16ZOrUqXnrrbe+dfxRRx2Vzz//PMccc0z69OmTZNHfi08//XT222+/jBs3Lvvvv38poBec0fKLX/wiW2+9dZLkoosuKv3v64wdOzb77bdfXn/99Rx88ME5/PDDM2nSpPTr1+8r1zI499xz8/rrr+fII4/MPvvsk8cff3yhRSMffvjh7LLLLt/6Xt99993zySefZPTo0TW2f/DBB/nb3/6W3Xbb7RsfD8DS4bR/AJYpn376aaZOnZru3bsv0vj33nsvQ4YMydFHH53DDjustH3HHXfMnnvumd///vc1tv9QvXv3ztChQ3P11Vdnhx12WGzX/ldVVeXWW28tXd/+8ccf58EHH0y3bt1y3XXXJfnykocJEyZk+PDhOfLII2s8/s0338x9992X9dZbL0my6667Zuedd86VV16Zq666Kkly2WWXpaqqKiNGjMhKK62UJNlnn31y7LHH5qqrrkrfvn2z/PLLl57ziy++yH333Vdj23f1xRdfZLfddkv//v2TfPnn161bt5xyyim59NJLs8suuyRJttpqq/Ts2TN//OMfv/Ia9ySZO3duLrvssqyzzjq55ZZbUlFRkST58Y9/nNNOOy2rrLLK957nAmuuuWaSLxebXPDPX2fttdfOJZdcUrq9qO/FqqqqnH766WnVqlX++Mc/pkmTJqWxC8442GijjbLGGmvk6aefXugSha9y+eWXZ+7cubnjjjvSrl27JMkee+yRnXfeORdffHFuu+22GuObNWuWG2+8sfT+nT9/fm699dbMmjUrjRs3/tbX+3dbbLFF2rRpk/vvvz8//elPS9sffPDBzJ8/X/wDLCMc+QdgmfLpp58m+fIo7KJ4+OGHM3/+/PTs2TPTp08v/a9FixZZffXVF+ln4r6Lfz/6/8gjjyy25919991rLGy3wQYbpLq6Oj/72c9qjNtggw0yZcqUzJs3r8b2jTbaqBT+SdK2bdt07949o0ePTlVVVaqrq/OXv/wl22+/faqrq2v8WXXt2jWzZs3Ka6+9VuM599hjj4XCf9VVV80bb7yRCy+8cJH3bcFR8SRp0qRJ2rdvnxVWWCE9e/Ysbe/QoUOaNGnyjUf+//Wvf+Wjjz5K3759S+GfJHvuued3Dtavs+B9N3v27G8d27dv3xq3F/W9+L//+7+ZNGlS9t9//xrhn+R7fZlUVVWVp59+Oj169CiFf5K0atUqvXr1yosvvlj6/9UCP//5z2u81qabbpqqqqq89957pW177bVX3njjjey1117f+Pr16tVL796989hjj9V4nfvvvz8bbbRRjTkBUHsc+QdgmdKoUaMkixZfSfLOO++kurr6axfgq19/8f9V17t37/zP//xPrr766vTo0WOxPGfbtm1r3F4Qs/95NLtx48aZP39+Zs2aVTp6n3y5DsJ/WmONNfL5559n+vTpqVevXmbOnJm77rord91111fO4T+v+V511VW/1778u+WWWy7NmzdfaB/atGmzUOg2btz4Gy/TmDx5cpKF97VBgwaLLTAXvO8W5cun//zzWdT34oIvONZaa60fMtWS6dOn5/PPP0/79u0Xuq9jx46ZP39+pkyZUuNMhv98vy34EuL7Xiazxx575LrrrssjjzySPfbYI+PHj89rr72Ws84663s9HwCLn/gHYJnSqFGjtGrVKmPHjl2k8fPnz09ZWVmuu+66r/zpsyWxAvuCo/8nn3xyHn300e/02Kqqqq+c59etov91279twcH/NH/+/CTJbrvtlj333PMrx/z7gnJJftDp/gt83c/Rfd3277pfi9uC991XfZnyn/5zxf3aeC9+X4vrfbXAj3/846y77rq5//77s8cee+T+++9PgwYNapzdAUDtEv8ALHN++tOf5q677spLL72UjTba6BvHrrbaaqmurs6qq676lUc+l5TddtstQ4cOzVVXXZXtt99+ofubNm36lUdRJ0+evEROg3733XcX2vbOO+9khRVWKB15b9iwYebPn5+tttpqsb/+0rDgaPW7776bLbfcsrR97ty5mTRpUtZee+0f9PyzZ8/OI488klVWWSUdO3b8zo9f1Pfign//b7755jf+u1jUSwCaN2+eFVZYIW+//fZC940fPz716tVbLOshfJs99tgjF154YaZNm5YHHngg2223XZo2bbrEXxeAReOafwCWOQcffHBWXHHFDBo0KB9++OFC90+YMKH0k2g77rhjysvLc9VVVy101LK6ujoff/zxEpnjgqP/Y8aMyWOPPbbQ/e3atcsrr7xS42f5Hn/88cX+84MLvPTSSzWu2Z8yZUoeffTRbL311ikvL095eXl22mmnPPTQQzV+AWCB/zzl/+t8l5/6W9zWW2+9NG/ePHfeeWeNP9c//OEPP/hXHebMmZMTTzwxn3zySQ477LDvde39or4X11133ay66qq55ZZbFpr3vz9uhRVWSPLtp+KXl5dn6623zqOPPppJkyaVtn/44Yd54IEHsskmm5Qup/kuFvWn/hbo1atXysrKct5552XixIkW+gNYxjjyD8AyZ7XVVsv/+3//L8ccc0x22WWX7L777llrrbVSWVmZl156KaNGjSotQrbaaqvl6KOPziWXXJL33nsvPXr0SMOGDTNp0qQ88sgj+fnPf56DDjpoicxzwbX/Y8aMWei+Pn365KGHHsrBBx+cnj17ZsKECRkxYkRWW221JTKXtdZaKwcddFD69euXioqK3HHHHUlSY+X84447Ls8991x+/vOfp0+fPvnxj3+cGTNm5LXXXsuzzz6b559//ltfZ8FP/e25557fadG/xaFBgwY5+uijc/rpp+eAAw7ILrvskkmTJmX48OHf6WyKqVOn5k9/+lOS5LPPPsu4ceMyatSofPDBBznwwAMXWshvUS3qe7FevXo588wzc/jhh2ePPfbIXnvtlZYtW2b8+PF56623csMNNyT58kuC5Muf5evatWvKy8uz6667fuVrH3300XnmmWey7777Zt999015eXnuuuuuVFZW5oQTTvhe+/Pwww9n4MCBueCCC7510b/kyzMQunXrllGjRqVJkybZbrvtvtfrArBkiH8Alkndu3fP/fffnxtuuCGPPvpo7rjjjlRUVKRTp045+eST8/Of/7w09tBDD80aa6yRm266KVdffXWSpE2bNtl6662/8pT8xaV+/fo5/PDDM3DgwIXu69atW04++eT87ne/y/nnn5/11lsv11xzTQYPHrxE5rLZZpulS5cuufrqqzN58uT8+Mc/zgUXXFDjVPgWLVrknnvuydVXX52HH344d9xxR5o1a5Yf//jHOf7445fIvBa3X/ziF6mqqsoNN9yQiy66KGuttVaGDh2aK664YpGfY8yYMTnxxBNTVlaWhg0bZpVVVslPf/rT9OnTJxtssMEPmt+ivhe7deuWm2++OVdffXVuvPHGVFdXp127djXe1zvuuGP69euXBx98MPfff3+qq6u/Nv7XXHPN3H777bnkkksybNiwVFdXZ4MNNsjFF1+cDTfc8Aft03ex++675/HHH0/Pnj1r/CIDALWvrLq2V9YBAKAQHnnkkRxxxBG5/fbbs+mmm9b2dAD4N675BwBgsbjnnnvSrl27bLLJJrU9FQD+Q50+7f8/f5Lo3z311FNp1apVkqSysjKXXXZZ7r///syePTsbbbRRTjvttHTo0KHGY8aNG5dzzz03L730Uho2bJjdd989Rx99tNPWAAC+wYMPPpg33ngjf/3rX3Pqqad+rwUTAViy6vRp/y+//PJC20466aSssMIK+eMf/1jadvrpp2fkyJE5+eST07p161xzzTWZOHFiHnzwwTRu3DhJMmPGjOy6665ZY4010r9//0ydOjUXXnhhdtttt5x++ulLaY8AAOqeTp06ZcUVV8wuu+ySs846K/Xr1+njSwCFVKc/mbt06VLj9qRJk/LOO+/UWNX2/fffz7333pszzjgje++9d5Jk/fXXz09/+tPceeedOeSQQ5Ikd955Z2bPnp2rrroqzZo1S5JUVVXlrLPOSv/+/dO6deulsk8AAHXNG2+8UdtTAOBbFOqa/wceeCBlZWXp1atXadvo0aMzf/787LzzzqVtzZo1y9Zbb50nn3yytO3JJ5/MlltuWQr/JOnZs2fmz5+fp59+eqnMHwAAAJaEOn3k/z89+OCD2WyzzdKmTZvStvHjx2fllVdO06ZNa4zt2LFj7r333hrjfvazn9UY06RJk9Lv7i4Jm266aSorK9OyZcsl8vwAAAAsWz744INUVFTkhRdeWKqvW5gj/6+//nrefPPNGkf9k2TmzJml6/r/XZMmTTJjxowa45o0abLQuKZNm9YYtzh98cUXmTdv3hJ5bgAAAJY98+bNyxdffLHUX3eZOvI/a9asTJs27VvHtWvXbqEV+EeMGJEGDRpkp512WlLTW+wW/BrBo48+WsszAQAAYGno3r17rbzuMhX/o0aNyqBBg7513MiRI9OxY8fS7erq6owcOTLdunWrcc1+8uUR/k8//XSh55g5c2aNSwGaNGmSWbNmLTRuxowZC10yAAAAAHXJMhX/ffr0SZ8+fb7z41588cVMnjy5xir/C3To0CEffvjhQhE/fvz4dOjQoca4/7y2f9asWfnggw9qjAMAAIC6phDX/I8YMSIrrrhitt9++4Xu69q1a+rVq5e//OUvpW0zZszI6NGjs80225S2bbPNNnnmmWcyc+bM0rZRo0alXr162XrrrZfsDgAAAMAStEwd+f8+5s2bl4ceeig9evTI8ssvv9D9bdq0yd57752LLroo9erVS+vWrTNs2LA0btw4ffv2LY3r27dvbr311hxxxBHp379/pk6dmosuuih9+/ZN69atl+YuAQAAwGJV5+N/9OjR+fjjjxda5f/fDRo0KA0bNswll1yS2bNnZ+ONN87vfve7Gr8C0LRp09x8880555xzcsQRR6Rhw4bZe++9c8wxxyyN3QAAAIAlpqy6urq6tifx32rBKo9W+wcAAFj2VVVVZe7cuV97f4MGDVJeXv6Nz1FbHVjnj/wDAADAklRdXZ0pU6Zk4sSJqayszFcdQy8rK0tFRUXatWuXVVZZJWVlZbUw068n/gEAAOAbTJkyJWPGjMkf/vCHTJgwIVVVVQuNKS8vz2qrrZY999wzSdK2bdulPc1vJP4BAADga1RVVWXixIn5wx/+kL/97W/fOHby5MlJkkaNGqV169bfegnA0lSIn/oDAACAJWHu3LmprKzMhAkTFmn8hAkTUllZ+Y1rA9QG8Q8AAADfoLq6+itP9f8qVVVVX7kmQG0T/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAN+grKxskX+2r7y8PGVlZUt4Rt+d+AcAAICv0aBBg1RUVGS11VZbpPGrrbZaKioq0qBBgyU8s++mfm1PAAAAAJZV5eXladeuXfbcc88kyYQJE77yZ//Ky8uz2mqrZc8990y7du0W+UyBpUX8AwAAwDdYZZVVkiSNGjVKZWVlqqurFxpTVlaWioqKtGvXrjR+WSL+AQAA4BuUlZWlbdu2ad26debOnfu14xo0aLDMHfFfQPwDAADAIigvL19m4/7bWPAPAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAAUn/gEAAKDgxD8AAAAUnPgHAACAghP/AAAAUHDiHwAAAApO/AMAAEDBiX8AAAAoOPEPAAAABSf+AQAAoODEPwAAABSc+AcAAICCE/8AAABQcOIfAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAAUn/gEAAKDgxD8AAAAUnPgHAACAghP/AAAAUHDiHwAAAApO/AMAAEDBiX8AAAAoOPEPAAAABSf+AQAAoODq1/YEfohOnTp97X1PPfVUWrVqlfHjx+e2227L3/72t7z33ntZeeWV061btxx11FFp3rx5afzw4cMzcODAhZ7nkEMOyfHHH79E5g8AAABLQ52O/7vuumuhbSeddFJWWGGFtGrVKknyzDPP5IUXXsgvfvGLrL322pk8eXKuvPLKPP/88/nTn/6UioqKGo+//vrr07hx49Lt1q1bL9mdAAAAgCWsTsd/ly5datyeNGlS3nnnnZxwwgmlbbvuumv222+/lJWVlbatvvrq2WefffL4449np512qvEc6667bo0zAgAAAKCuq9Px/58eeOCBlJWVpVevXqVtK6200kLjOnfunCSZNm3aUpsbAAAA1JZCLfj34IMPZrPNNkubNm2+cdyLL76YJOnYseNC9/Xq1SvrrLNOunfvnmHDhqWqqmqJzBUAAACWlsIc+X/99dfz5ptv5uyzz/7GcV988UUGDx6czp07Z8sttyxtb9myZQYMGJANN9wwZWVleeyxx3L55Zdn6tSpOf3005f09AEAAGCJWabif9asWYt0Kn67du0WWqhvxIgRadCgwULX8P+nM844I5MmTcqdd95ZYx2Abt26pVu3bqXbXbt2zXLLLZebb745hx12WGkBQQAAAKhrlqn4HzVqVAYNGvSt40aOHFnjlP3q6uqMHDky3bp1S7Nmzb72cZdddllGjBiRa665Jmuttda3vk7Pnj1z4403ZsyYMeIfAACAOmuZiv8+ffqkT58+3/lxL774YiZPnlxjlf//dOutt2bYsGG58MILaxzhBwAAgKIrxIJ/I0aMyIorrpjtt9/+K+9/4IEHct555+XYY4/NHnvsscjPO3LkyJSXl5d+HQAAAADqomXqyP/3MW/evDz00EPp0aNHll9++YXuf/7553PyySdniy22yE9+8pO8/PLLpfvatGlT+mWAgw46KJtvvnk6deqUJHn00Udz9913Z//990/Lli2Xyr4AAADAklDn43/06NH5+OOP06tXr6+8/7nnnsvcuXPz7LPP5tlnn61x35FHHpkBAwYkSdq3b5/77rsv77//fubPn5811lgjp5xySvr167fE9wEAAACWpLLq6urq2p7Ef6vu3bsn+fIsAwAAAIqvtjqwENf8AwAAAF9P/AMAAEDBiX8AAAAoOPEPAAAABSf+AQAAoODEPwAAABSc+AcAAICCE/8AAABQcOIfAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAAUn/gEAAKDgxD8AAAAUnPgHAACAghP/AAAAUHDiHwAAAApO/AMAAEDBiX8AAAAoOPEPAAAABSf+AQAAoODEPwAAABSc+AcAAICCE/8AAABQcOIfAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAAUn/gEAAKDgxD8AAAAUnPgHAACAghP/AAAAUHDiHwAAAApO/AMAAEDBiX8AAAAoOPEPAAAABSf+AQAAoODEPwAAABSc+AcAAICCE/8AAABQcOIfAAAACk78AwAAQMGJfwAAACi4+rU9gR+iU6dOX3vfU089lVatWn3tuBYtWuTpp5+usW3cuHE599xz89JLL6Vhw4bZfffdc/TRR6eiomLxThwAAACWojod/3fddddC20466aSssMIKpfBfoF+/funVq1fpdoMGDWrcP2PGjBxwwAFZY401MmTIkEydOjUXXnhh5syZk9NPP33J7AAAAAAsBXU6/rt06VLj9qRJk/LOO+/khBNOWGjsKqusstD4f3fnnXdm9uzZueqqq9KsWbMkSVVVVc4666z0798/rVu3XowzBwAAgKWnUNf8P/DAAykrK6txhH9RPfnkk9lyyy1L4Z8kPXv2zPz58xe6PAAAAADqkkLF/4MPPpjNNtssbdq0Wei+a6+9Nuuuu2423XTTHH300Zk8eXKN+8ePH58OHTrU2NakSZO0bNky48ePX6LzBgAAgCWpTp/2/+9ef/31vPnmmzn77LMXum+PPfbIdtttlxYtWuTNN9/M0KFDs+++++ZPf/pTmjZtmiSZOXNmmjRpstBjmzZtmhkzZizx+QMAAMCSskzF/6xZszJt2rRvHdeuXbuFVuAfMWJEGjRokJ122mmh8YMHDy7982abbZZNNtkke+21V+6+++4ccsghP3ziAAAAsAxbpuJ/1KhRGTRo0LeOGzlyZDp27Fi6XV1dnZEjR6Zbt241rtn/OmuvvXbat2+f1157rbStSZMmmTVr1kJjZ8yYUTo7AAAAAOqiZSr++/Tpkz59+nznx7344ouZPHnyV67yv6g6dOiw0LX9s2bNygcffLDQWgAAAABQlxRiwb8RI0ZkxRVXzPbbb79I48eMGZO3334766+/fmnbNttsk2eeeSYzZ84sbRs1alTq1auXrbfeerHPGQAAAJaWZerI//cxb968PPTQQ+nRo0eWX375he6/4YYbMmHChGy++eZp3rx5xo4dm2uuuSZt2rSpcZZB3759c+utt+aII45I//79M3Xq1Fx00UXp27dvWrduvTR3CQAAABarOh//o0ePzscff5xevXp95f3t27fPX/7yl/z5z3/O7Nmzs9JKK2XbbbfN0UcfXWN1/6ZNm+bmm2/OOeeckyOOOCINGzbM3nvvnWOOOWZp7QoAAAAsEWXV1dXVtT2J/1bdu3dPkjz66KO1PBMAAACWhtrqwEJc8w8AAAB8PfEPAAAABSf+AQAAoODEPwAAABSc+AcAAICCE/8AAABQcOIfAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAAUn/gEAAKDgxD8AAAAUnPgHAACAghP/AAAAUHDiHwAAAApO/AMAAEDBiX8AAAAoOPEPAAAABSf+AQAAoODEPwAAABSc+AcAAICCE/8AAABQcOIfAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAAUn/gEAAKDgxD8AAAAUnPgHAACAghP/AAAAUHDiHwAAAApO/AMAAEDBiX8AAAAoOPEPAAAABSf+AQAAoODEPwAAABSc+AcAAICCE/8AAABQcOIfAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFFz92p7AD9GpU6evve+pp55Kq1atMnz48AwcOPArx3Tt2jU33HBDknztuEMOOSTHH3/84pkwAAAA1II6Hf933XXXQttOOumkrLDCCmnVqlWSZLvttlto3DvvvJOTTjop22yzzUKPv/7669O4cePS7datWy/mWQMAAMDSVafjv0uXLjVuT5o0Ke+8805OOOGE0rbmzZunefPmNcY99dRTKS8vzy677LLQc6677roLjQcAAIC6rFDX/D/wwAMpKytLr169vnXcFltskZYtWy6lmQEAAEDtKVT8P/jgg9lss83Spk2brx3zz3/+M++8887XfkHQq1evrLPOOunevXuGDRuWqqqqJTVdAAAAWCrq9Gn//+7111/Pm2++mbPPPvsbxz3wwANZbrnlsuOOO9bY3rJlywwYMCAbbrhhysrK8thjj+Xyyy/P1KlTc/rppy/JqQMAAMAStUzF/6xZszJt2rRvHdeuXbtUVFTU2DZixIg0aNAgO+2009c+bv78+XnwwQez3XbbpVGjRjXu69atW7p161a63bVr1yy33HK5+eabc9hhh5UWEAQAAIC6ZpmK/1GjRmXQoEHfOm7kyJHp2LFj6XZ1dXVGjhyZbt26pVmzZl/7uOeeey4ffPBBevfuvUjz6dmzZ2688caMGTNG/AMAAFBnLVPx36dPn/Tp0+c7P+7FF1/M5MmTa6zy/1VGjBiRJk2aZNttt/2+UwQAAIA6pxAL/o0YMSIrrrhitt9++68dU1lZmYcffjg77LDDQpcMfJ2RI0emvLw8nTt3XlxTBQAAgKVumTry/33MmzcvDz30UHr06JHll1/+a8c98cQTmTlz5tee8n/QQQdl8803T6dOnZIkjz76aO6+++7sv//+fhIQAACAOq3Ox//o0aPz8ccff+1P9y0wYsSItGzZMptvvvlX3t++ffvcd999ef/99zN//vysscYaOeWUU9KvX78lMW0AAABYasqqq6ura3sS/626d++e5MuzDAAAACi+2urAQlzzDwAAAHw98Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAAUn/gEAAKDgxD8AAAAUnPgHAACAghP/AAAAUHDiHwAAAApO/AMAAEDBiX8AAAAoOPEPAAAABSf+AQAAoODEPwAAABSc+AcAAICCE/8AAABQcOIfAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAAUn/gEAAKDgxD8AAAAUnPgHAACAghP/AAAAUHDiHwAAAApO/AMAAEDBiX8AAAAoOPEPAAAABSf+AQAAoODEPwAAABSc+AcAAICCE/8AAABQcOIfAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAAUn/gEAAKDg6nT8V1VV5brrrsvOO++cDTfcMN27d8/gwYMze/bsGuMqKyszePDgbL311unSpUt+/etfZ/z48Qs937hx4/LrX/86Xbp0ydZbb52LLroolZWVS2t3AAAAYImo0/E/dOjQXH755dlrr70ybNiw/OpXv8qdd96Z008/vca4c889N/fcc0+OOeaYDBkyJJWVlfnVr36VWbNmlcbMmDEjBxxwQObOnZshQ4bkmGOOyd13350LL7xwae8WAAAALFb1v8vgq666Kl26dEnXrl2/dezLL7+cu+66KxdccMH3nty3eeCBB9K7d+8ceuihSZItttgiH3/8ca677rrMmzcv9evXz/vvv5977703Z5xxRvbee+8kyfrrr5+f/vSnufPOO3PIIYckSe68887Mnj07V111VZo1a5bkyzMLzjrrrPTv3z+tW7deYvsBAAAAS9J3OvJ/1VVX5ZBDDsn555+fuXPnfuPYCRMm5I9//OMPmdu3mjdvXho1alRjW+PGjVNdXV26PXr06MyfPz8777xzaVuzZs2y9dZb58knnyxte/LJJ7PllluWwj9Jevbsmfnz5+fpp59ecjsBAAAAS9h3Pu2/Xbt2ueWWW9K3b99MmDBhScxpkfXp0yf3339/nn322cyePTuvvvpqbr311vTt2zf16395UsP48eOz8sorp2nTpjUe27FjxxrX/Y8fPz4dOnSoMaZJkyZp2bLlV64PAAAAAHXFdzrtP0kGDBiQOXPm5Lzzzsuee+6ZM844I7vtttuSmNu36t+/fyorK/PrX/+6dLR/t912yymnnFIaM3PmzDRu3HihxzZp0iQzZsyoMa5JkyYLjWvatGmNcQAAAFDXfOf4T7484t6lS5ccffTROemkk/Lss8/m9NNPzworrPCDJjNr1qxMmzbtW8e1a9cuFRUVue2223LLLbdk4MCB6dy5c8aOHZsrrrgi55xzTs4444wfNBcAAAAoiu8V/0my5ppr5r777svZZ5+d4cOH55VXXsmll16atdde+3tPZtSoURk0aNC3jhs5cmSaN2+ewYMH58QTT0y/fv2SJJtttlkaNWqUE044Ifvvv3/at2+fJk2a5NNPP13oOWbOnFnjUoAmTZrUWP1/gRkzZix0yQAAAADUJd87/pNk+eWXz/nnn58tt9wyZ555Zn7xi1/kxBNPzH777fe9nq9Pnz7p06fPIo199dVXU1lZmXXWWafG9s6dOyf5csHB9u3bp0OHDvnwww8Xivj/vMa/Q4cOC13bP2vWrHzwwQcLrQUAAAAAdcl3XvDvq/Tu3Tv33XdfOnTokHPPPTdHHnnkEr9Ovm3btkmS1157rcb2f/3rX0mSVVddNUnStWvX1KtXL3/5y19KY2bMmJHRo0dnm222KW3bZptt8swzz2TmzJmlbaNGjUq9evWy9dZbL7H9AAAAgCXtBx35/3drrLFG7rrrrgwePDi33357/vrXvy6up/5KLVq0SI8ePXLFFVekqqoqnTt3zltvvZUhQ4Zkq622SseOHZMkbdq0yd57752LLroo9erVS+vWrTNs2LA0btw4ffv2LT1f3759c+utt+aII45I//79M3Xq1Fx00UXp27dvWrduvUT3BQAAAJakxRb/SVJRUZHTTjstW2yxRU499dSvvIZ+cRo8eHCuvvrq3HHHHZk6dWpatmyZ3r17Z8CAATXGDRo0KA0bNswll1yS2bNnZ+ONN87vfve7Gr8C0LRp09x8880555xzcsQRR6Rhw4bZe++9c8wxxyzRfQAAAIAlrax6wW/kLWYffPBB3n777fzkJz9ZEk9fCN27d0+SPProo7U8EwAAAJaG2urAxXrk/9+1bNkyLVu2XFJPDwAAACyixbLgHwAAALDsEv8AAABQcOIfAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAAUn/gEAAKDgxD8AAAAUnPgHAACAghP/AAAAUHDiHwAAAApO/AMAAEDBiX8AAAAoOPEPAAAABSf+AQAAoODEPwAAABSc+AcAAICCE/8AAABQcOIfAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAAUn/gEAAKDgxD8AAAAUnPgHAACAghP/AAAAUHDiHwAAAApO/AMAAEDBiX8AAAAoOPEPAAAABSf+AQAAoODEPwAAABSc+AcAAICCE/8AAABQcOIfAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKrn5tT+CHqKqqyo033pj77rsvU6ZMSYsWLbLjjjvmyCOPTMOGDZMk48ePz2233Za//e1vee+997LyyiunW7duOeqoo9K8efPScw0fPjwDBw5c6DUOOeSQHH/88UttnwAAAGBxq9PxP3To0AwdOjRHHXVUNthgg4wdOzaXXnpppk2blksuuSRJ8swzz+SFF17IL37xi6y99tqZPHlyrrzyyjz//PP505/+lIqKihrPef3116dx48al261bt16q+wQAAACLW52O/wceeCC9e/fOoYcemiTZYost8vHHH+e6667LvHnzUr9+/ey6667Zb7/9UlZWVnrc6quvnn322SePP/54dtpppxrPue6669Y4IwAAAADqujod//PmzUujRo1qbGvcuHGqq6tLt1daaaWFHte5c+ckybRp05bsBAEAAGAZUKcX/OvTp0/uv//+PPvss5k9e3ZeffXV3Hrrrenbt2/q1//67zVefPHFJEnHjh0Xuq9Xr15ZZ5110r179wwbNixVVVVLbP4AAACwNNTpI//9+/dPZWVlfv3rX5eO9u+222455ZRTvvYxX3zxRQYPHpzOnTtnyy23LG1v2bJlBgwYkA033DBlZWV57LHHcvnll2fq1Kk5/fTTl/i+AAAAwJKyTMX/rFmzFulU/Hbt2qWioiK33XZbbrnllgwcODCdO3fO2LFjc8UVV+Scc87JGWec8ZWPPeOMMzJp0qTceeedNdYB6NatW7p161a63bVr1yy33HK5+eabc9hhh6VVq1Y/fAcBAACgFixT8T9q1KgMGjToW8eNHDkyzZs3z+DBg3PiiSemX79+SZLNNtssjRo1ygknnJD9998/7du3r/G4yy67LCNGjMg111yTtdZa61tfp2fPnrnxxhszZswY8Q8AAECdtUzFf58+fdKnT59FGvvqq6+msrIy66yzTo3tCxbzmzBhQo34v/XWWzNs2LBceOGFNY7wAwAAQNHV2QX/2rZtmyR57bXXamz/17/+lSRZddVVS9seeOCBnHfeeTn22GOzxx57LPJrjBw5MuXl5aUvFAAAAKAuWqaO/H8XLVq0SI8ePXLFFVekqqoqnTt3zltvvZUhQ4Zkq622Kq3k//zzz+fkk0/OFltskZ/85Cd5+eWXS8/Rpk2btGnTJkly0EEHZfPNN0+nTp2SJI8++mjuvvvu7L///mnZsuVS3z8AAABYXOps/CfJ4MGDc/XVV+eOO+7I1KlT07Jly/Tu3TsDBgwojXnuuecyd+7cPPvss3n22WdrPP7II48sjW3fvn3uu+++vP/++5k/f37WWGONnHLKKaX1BAAAAKCuKqte8Bt5LHXdu3dP8uVZBgAAABRfbXVgnb3mHwAAAFg04h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAAUn/gEAAKDgxD8AAAAUnPgHAACAghP/AAAAUHDiHwAAAApO/AMAAEDBiX8AAAAoOPEPAAAABSf+AQAAoODEPwAAABSc+AcAAICCE/8AAABQcOIfAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAAUn/gEAAKDgxD8AAAAUnPgHAACAghP/AAAAUHDiHwAAAApO/AMAAEDBiX8AAAAoOPEPAAAABSf+AQAAoODEPwAAABSc+AcAAICCE/8AAABQcOIfAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAAUn/gEAAKDgxD8AAAAUnPgHAACAghP/AAAAUHDiHwAAAApO/AMAAEDB1a/tCfwQVVVVufHGG3PfffdlypQpadGiRXbccccceeSRadiwYWlcp06dFnpsixYt8vTTT9fYNm7cuJx77rl56aWX0rBhw+y+++45+uijU1FRscT3BQAAAJaUOh3/Q4cOzdChQ3PUUUdlgw02yNixY3PppZdm2rRpueSSS2qM7devX3r16lW63aBBgxr3z5gxIwcccEDWWGONDBkyJFOnTs2FF16YOXPm5PTTT18q+wMAAABLQp2O/wceeCC9e/fOoYcemiTZYost8vHHH+e6667LvHnzUr/+/+3eKquski5dunztc915552ZPXt2rrrqqjRr1izJl2cWnHXWWenfv39at269JHcFAAAAlpg6fc3/vHnz0qhRoxrbGjdunOrq6u/8XE8++WS23HLLUvgnSc+ePTN//vyFLg8AAACAuqROx3+fPn1y//3359lnn83s2bPz6quv5tZbb03fvn1rHPVPkmuvvTbrrrtuNt100xx99NGZPHlyjfvHjx+fDh061NjWpEmTtGzZMuPHj1/i+wIAAABLSp0+7b9///6prKzMr3/969LR/t122y2nnHJKjXF77LFHtttuu7Ro0SJvvvlmhg4dmn333Td/+tOf0rRp0yTJzJkz06RJk4Veo2nTppkxY8aS3xkAAABYQpap+J81a1amTZv2rePatWuXioqK3HbbbbnlllsycODAdO7cOWPHjs0VV1yRc845J2eccUZp/ODBg0v/vNlmm2WTTTbJXnvtlbvvvjuHHHLIEtkXAAAAWFYsU/E/atSoDBo06FvHjRw5Ms2bN8/gwYNz4oknpl+/fkm+DPtGjRrlhBNOyP7775/27dt/5ePXXnvttG/fPq+99lppW5MmTTJr1qyFxs6YMaN0dgAAAADURctU/Pfp0yd9+vRZpLGvvvpqKisrs84669TY3rlz5yTJhAkTvjb+v0qHDh0WurZ/1qxZ+eCDDxZaCwAAAADqkjq74F/btm2TpMbR+yT517/+lSRZddVVv/axY8aMydtvv53111+/tG2bbbbJM888k5kzZ5a2jRo1KvXq1cvWW2+9OKcOAAAAS9UydeT/u2jRokV69OiRK664IlVVVencuXPeeuutDBkyJFtttVU6duyYJLnhhhsyYcKEbL755mnevHnGjh2ba665Jm3atKlxlkHfvn1z66235ogjjkj//v0zderUXHTRRenbt29at25dW7sJAAAAP1idjf/ky4X8rr766txxxx2ZOnVqWrZsmd69e2fAgAGlMe3bt89f/vKX/PnPf87s2bOz0korZdttt83RRx9dY3X/pk2b5uabb84555yTI444Ig0bNszee++dY445pjZ2DQAAABabsuoFv5HHUte9e/ckyaOPPlrLMwEAAGBpqK0OrLPX/AMAAACLRvwDAABAwYl/AAAAKDjxDwAAAAUn/gEAAKDgxD8AAAAUnPgHAACAghP/AAAAUHDiHwAAAApO/AMAAEDBiX8AAAAoOPEPAAAABSf+AQAAoODEPwAAABSc+AcAAICCE/8AAABQcOIfAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAAUn/gEAAKDgxD8AAAAUnPgHAACAghP/AAAAUHDiHwAAAApO/AMAAEDBiX8AAAAoOPEPAAAABSf+AQAAoODEPwAAABSc+AcAAICCE/8AAABQcOIfAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAAUn/gEAAKDgxD8AAAAUnPgHAACAghP/AAAAUHDiHwAAAApO/AMAAEDBiX8AAAAoOPEPAAAABVe/tifwQ1RVVeXGG2/MfffdlylTpqRFixbZcccdc+SRR6Zhw4ZJkuHDh2fgwIFf+fiuXbvmhhtu+MZxhxxySI4//vgltxMAAACwhNXp+B86dGiGDh2ao446KhtssEHGjh2bSy+9NNOmTcsll1ySJNluu+1y11131XjcO++8k5NOOinbbLPNQs95/fXXp3HjxqXbrVu3XrI7AQAAAEtYnY7/Bx54IL17986hhx6aJNliiy3y8ccf57rrrsu8efNSv379NG/ePM2bN6/xuKeeeirl5eXZZZddFnrOddddd6HxAAAAUJfV6Wv+582bl0aNGtXY1rhx41RXV3/j4x544IFsscUWadmy5ZKcHgAAACwT6nT89+nTJ/fff3+effbZzJ49O6+++mpuvfXW9O3bN/Xrf/VJDf/85z/zzjvvpFevXl95f69evbLOOuuke/fuGTZsWKqqqpbkLgAAAMASV6dP++/fv38qKyvz61//unS0f7fddsspp5zytY954IEHstxyy2XHHXessb1ly5YZMGBANtxww5SVleWxxx7L5ZdfnqlTp+b0009fovsBAAAAS9IyFf+zZs3KtGnTvnVcu3btUlFRkdtuuy233HJLBg4cmM6dO2fs2LG54oorcs455+SMM85Y6HHz58/Pgw8+mO22226hywW6deuWbt26lW537do1yy23XG6++eYcdthhadWq1Q/fQQAAAKgFy1T8jxo1KoMGDfrWcSNHjkzz5s0zePDgnHjiienXr1+SZLPNNkujRo1ywgknZP/990/79u1rPO65557LBx98kN69ey/SfHr27Jkbb7wxY8aMEf8AAADUWctU/Pfp0yd9+vRZpLGvvvpqKisrs84669TY3rlz5yTJhAkTFor/ESNGpEmTJtl2220Xz4QBAACgDqizC/61bds2SfLaa6/V2P6vf/0rSbLqqqvW2F5ZWZmHH344O+ywQyoqKhbpNUaOHJny8vLSFwoAAABQFy1TR/6/ixYtWqRHjx654oorUlVVlc6dO+ett97KkCFDstVWW6Vjx441xj/xxBOZOXPm157yf9BBB2XzzTdPp06dkiSPPvpo7r777uy///5+EhAAAIA6rc7Gf5IMHjw4V199de64445MnTo1LVu2TO/evTNgwICFxo4YMSItW7bM5ptv/pXP1b59+9x33315//33M3/+/Kyxxho55ZRTSusJAAAAQF1VVr3gN/JY6rp3757ky7MMAAAAKL7a6sA6e80/AAAAsGjEPwAAABSc+AcAAICCE/8AAABQcOIfAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAAUn/gEAAKDgxD8AAAAUnPgHAACAghP/AAAAUHDiHwAAAApO/AMAAEDBiX8AAAAoOPEPAAAABSf+AQAAoODEPwAAABSc+AcAAICCE/8AAABQcOIfAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAAUn/gEAAKDgxD8AAAAUnPgHAACAghP/AAAAUHDiHwAAAApO/AMAAEDBiX8AAAAoOPEPAAAABSf+AQAAoODEPwAAABSc+AcAAICCE/8AAABQcOIfAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIKr0/FfXV2d6667Lttvv33WW2+99OrVKyNHjvzKcddee2222267bLDBBvnFL36Rl19+eaFxU6dOzYABA7LRRhvlJz/5SU499dR8+umnS2FPAAAAYMmp0/F//fXX5/LLL89ee+2Va665Jj/5yU9y7LHH5rHHHqsx7rrrrsuVV16ZX/3qVxk2bFhatmyZAw88MBMnTiyNmTt3bg4++OC88847ueSSS3LmmWdm9OjROe6445b2bgEAAMBiVb+2J/B9VVZWZujQoenXr1+OPPLIJEnXrl0zefLkXH755dl+++2TJF988UWGDRuWAw88ML/61a+SJJtsskl23nnn3HDDDTnzzDOTJA899FDGjh2bkSNHpkOHDkmSJk2a5KCDDsqrr76aDTbYYKnvIwAAACwOdfbI/8SJEzN79uxsvfXWNbZ37do1b7zxRiZPnpwk+cc//pFPP/00PXv2LI2pqKjIDjvskCeffLK07cknn0ynTp1K4Z8kW2+9dZo1a5YnnnhiCe8NAAAALDl19sj/F198keTLkP93C26PGzcubdu2zfjx45OkRtQnSceOHXPzzTdnzpw5WX755TN+/PiFxpSVlaV9+/al51jcpk2blqqqqnTv3n2JPD8AAADLlilTpqS8vHypv26dPfK/2mqrpaysLK+++mqN7QsW8psxY0aSZObMmamoqMhyyy1XY1yTJk1SXV1dY1zjxo0Xep2mTZuWxixuyy23XOrXr7PfvwAAAPAd1a9ff6E+XSqvu9Rf8RvMmjUr06ZN+9Zx7dq1S6NGjbLbbrvl+uuvz1prrZUuXbrk8ccfz4MPPpjky6P2y7oXXnihtqcAAADAf4FlKv5HjRqVQYMGfeu4kSNHpmPHjhk4cGA+/PDDHHrooUmSlVZaKUcddVQGDx6cli1bJvnyCH9lZWW++OKLGt+uzJw5M2VlZWnatGlp3Ff9rN+MGTOyyiqrLI7dAwAAgFqxTMV/nz590qdPn0Uev9JKK+XGG2/M1KlTM2PGjKyxxhp59NFH06BBg3Tu3DnJ/13r//bbb2fttdcuPXb8+PFp27Ztll9++dK4N998s8bzV1dX5+23315oUUEAAACoS+rsNf//rnXr1llrrbVSXl6eO+64I7vssksaNWqUJNl4443TqFGj/PnPfy6Nnzt3bv7yl79km222KW3bZptt8vrrr+edd94pbXv22WfzySefZNttt11q+wIAAACLW1l1dXV1bU/i+7r//vvzxRdfZLXVVsu0adNy1113ZdKkSbn77rvTqlWr0rhrr702Q4YMyfHHH5+11lord9xxR0aPHp0//elPadeuXZIvvxDYa6+9kiTHHntsPv/881x00UXp1KlThg0bViv7BwAAAItDnY7/P/3pT7nmmmsyadKkrLjiitl2221z7LHHpk2bNjXGVVdX59prr83vf//7TJ8+Peuss04GDhyYjTbaqMa4qVOn5txzz83o0aNTv3797LDDDjnllFNKZxEAAABAXVSn4x8AAAD4doW45h8AAAD4euIfAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKDjxDwAAAAUn/gEAAKDgxD8AAAAUnPivBVVVVbU9BVhqqqura3sKAADwX69+bU/gv8lnn32WIUOGZPbs2VlnnXWy1VZbZfXVV6/tacFSU1VVlfLy8tqeBnwvTz31VJo3b54WLVqkdevWtT0d+N4qKytTUVGRJJk/f37q1XMsiLrrqaeeykorrZSWLVv6bKZwqqurU1ZWttier6zaYbml4oUXXsjRRx+dH/3oR2nWrFnGjBmTBg0aZOjQoVlzzTUX679UWBY888wzueGGG9KkSZP86Ec/yuGHH56GDRvW9rTgO/vnP/+ZE044IbNmzUpZWVmqq6tz5plnZuutt86KK6642P9ihiXlmWeeyf/8z/9khRVWSIsWLfLb3/42bdq08f6lTvq2z2ao62677ba899576dChQzbZZJN06NDhB39hK/6XkgsuuCCvvfZaLr/88jRt2jTvv/9+fvOb36SioiIDBw7MpptuWttThB9sQQTddNNNGTJkSHbaaafMmTMnTz/9dDp06JDf/va32XLLLcUSdUZlZWUOO+ywNG3aNAMGDMisWbNy22235e9//3v222+/HHLIIbU9RfhGCz5v77333lx88cXZaaedUlFRkYcffjiNGzfOwQcfnD322MPnMnXKN30277vvvjn00ENre4rwvU2cODG//e1v88knn6RTp075xz/+kUaNGuWyyy7Lhhtu+IOe23leS8lDDz2UTTbZJC1atEiDBg3Srl27XHHFFZk2bVruvvvufPTRR7U9RfjBysrKUlVVlYceeih77rlnzj///Fx66aW58847M3v27Fx11VV59913S9/Qw7Lu9ddfzz//+c907949HTp0yIYbbpjzzz8/m2yySe69996MHj06yZenTsOyaEHQjxo1Kl27ds1ZZ52VQYMGZfjw4Vl11VUzZMiQvPHGGz6XqVO+6bP5vvvuy1NPPVXbU4Tv7cknn8y8efMybNiwDB06NLfcckt+/OMf55hjjsmLL774g55b/C8F8+bNS7t27TJhwoTStqqqqnTo0CG//OUv8+KLL+bhhx+uxRnC4vP222/njTfeyNZbb53kyyhq3759Dj/88MyaNStXXnllkjjCxDLt32N+9uzZWXfddZMkX3zxRRo0aJADDjggrVq1ylVXXZUkrplmmTZ58uSMGzcuG220USnyV1555Rx88MFZaaWVcsYZZyTxucyyb1E/m6+++uramiJ8J1/1petjjz2WVq1aZa211kpZWVnWXnvtXHrppamsrMwNN9yQiRMnfu/X818rS0H9+vWz9tprZ+rUqXnrrbeS/N+/6F//+tdp2rRpnnrqqUydOrU2pwnf2fTp0/PZZ5/V2Na6detUVFTk7bffTvLll19JsuOOO6Z79+557rnn8sQTTyRxtJRlx6uvvpobb7wxo0aNyowZM0oxv9JKK6VVq1b54x//mOTLz/Mk2WCDDbLzzjvn/fffz913353E+5llw9SpUzN9+vQa29q2bZuysrK89957SZK5c+cmSTbZZJPss88+efXVV/Pggw8m8T5m2bI4Ppuhrpg3b14qKyvTrFmzVFVVpbKyMsmXn9mNGjXKaaedltGjR+evf/3r9/6sFv9LycEHH5zXXnstf//731NdXZ369etn7ty5qaioyD777JNnn30206ZNq+1pwiL55JNPcswxx2S//fbLL3/5y1x33XWlM1sqKyuzySab5M9//nNpRel58+alvLw8O+64Y1ZbbbXce++9SRwtpfZVVlbmrLPOyi9/+cv8+c9/znHHHZf+/ftn+PDhSZKGDRtmww03zOjRozN16tSUl5eXwmnbbbdNhw4d8vjjj+ezzz7zfqZWzZgxI0cffXT23Xff7LPPPjn33HPz2muvJfnyfb7tttvmgQceqPG5XFZWlq222irbbLNNrrvuuiQ+l1k2LOpn89NPP/2Vn83t27cvfTbDsuiZZ57JQQcdlGOPPTb/7//9v8yePTv169dPRUVFVlpppbz//vsZO3ZskqRBgwZJkp122ildunTJ/fffn8mTJ3+v1/UJv5S0bt06e+yxR+6555688cYbSf7vL9gFC+28/vrrSfwuOsu2jz76KAMGDMj777+fgw46KOutt17uuOOOHHXUUZkzZ05WXnnlbLTRRvnkk09y//33J/m/9/o666yTDTfcMBMnTsyYMWNqczcgSfLSSy/l2WefzSWXXJLf/e53+cMf/pBWrVrljDPOyMsvv5zmzZuna9eumTFjRulLqwYNGqS6ujqrrrpq1ltvvbz77rvWbaFWzZ49OyeddFLee++9nHjiienZs2f+8Y9/5LDDDsvEiRNTUVGRzTffPEly/fXXJ/m/U/xXWWWVdOvWLZ9//nlpDQuobYv62fzJJ5985Wfz+uuvn3fffTcffvhhLe8J/J8FjXfTTTdlwIABad26dcrKynLPPffkkEMOKa1Vsd9+++Wdd97Js88+W+Pof5Icd9xx+ec//1k69f+7dqP4X4pOPvnkjB07Nn/4wx/yySeflH7vfNasWWnVqlXpND3X3LEse+211zJmzJgMHDgwe++9d84+++ycdtpp+fjjj3PiiScmSXbdddesvPLKGTFiRKZOnZp69erliy++SPLlt5ZvvfVW6XIAqE0vv/xyZs6cma5du6ZRo0ZZa621cvzxx2fTTTfNoEGDMmvWrOy222758Y9/nIcffjivvPJKkpT+Mt55550zfvz42twFyNixY/P888/nqKOOyk477ZSjjz46559/flZaaaWcfPLJmT59erbaaqtssskm+eMf/5iJEyemvLy89D7eYost8tFHH5U+p6G2LAiZH/rZ3LNnz4wfP95/U7NM+bqFse+66658+umnueaaazJu3Li0b98+vXv3zu9///vSJeMNGjRIVVVV1llnnay55pp55JFHSs/5XYj/pWiFFVbIxRdfnMcee6z0zXtVVVXGjx+fzz77LJtsskktzxC+3dSpU9O4ceO0atWqtK1r16458cQT8/DDD+ePf/xjVlllleyxxx6ZMmVK6VTS5ZZbLsmX1+nVq1fPt/EsdWPHjs0FF1yQ999/v7Tto48+Sps2bTJ79uzSttVWWy3HH3983nvvvfzud7/Lcsstl1/+8pcpKyvLRRddlOT/3s8vvvhi2rRpk3r16jlri6Xm7bffzogRI0q3P//88zRo0CBt2rQpbVt77bVzzjnn5KWXXsodd9yRJk2aZM8990y9evVy4YUXJkkqKiqSJK1atcrcuXMzc+bMpbsjkGTOnDmlf14QMt/1s3nw4MFJ/u+z+YUXXvDZzDLpqxbGXmONNXL44Ydn5syZpcUqjzvuuHz88ce5/fbbS+vCLThw3LBhw5SVlX2v6/7F/1K20047pXfv3rnpppuy9957Z9CgQTnmmGOy3nrrpWPHjj6gWGZ89tlnueGGG3LTTTflz3/+c6ZMmZLky7+YP/nkk3z00Uel92uDBg3So0eP7LTTThk6dGg+++yz7LXXXtl+++1zzz335K677sqsWbMyf/78PPTQQ1l99dWz8cYb1+bu8V9m/vz5Oeigg3LzzTdn1KhRNRY8e/PNN/PBBx/UGL/uuutmv/32y+23357p06ena9euOeCAAzJu3LgceuihufPOO/Pkk0/m3nvvzRZbbJEf/ehHjjCxxH3++ecZOHBgfvGLX+SGG27Iq6++muTLo53z5s3LO++8U/pcrq6uzoYbbph99903t99+e8aNG5euXbtmv/32y+OPP54rrrgikyZNyty5c3P//fdnlVVWKV0aAEvDZ599lgsvvDBHHnlkjj/++Dz55JOl+77rZ/P48eN9NrPM+a4LY/fo0SPPP/98Hn/88bRp0ybHH398/vznP+emm24qPX7q1Kn54IMP0qlTp++1RktZtdqsFU899VSeeuqpTJgwIV27ds0vf/nL2p4SlNx999255JJLssoqq6S6ujrvvvtuOnTokOHDh+fTTz/NT3/60xxwwAE58sgjM3/+/NKHzz//+c/st99+OeOMM/Kzn/0sU6ZMye9///vcfPPN+dGPfpQf/ehHee655/Kb3/wmhx12WBKXubDkzZ8/P5WVlRkwYECmT5+eiRMn5o477kjHjh0zffr09OnTJ1tuuWXOPffcVFdXl96T48aNyz777JODDz44hx56aObOnZtXXnklF198cT766KN8/vnn2XLLLXPOOedkhRVWqOW9pOg+/fTTnHrqqZkyZUqOOOKIrLzyyllrrbVKR++33XbbbLPNNjnttNNSUVFRei9Pnz49O+64Y371q1/lyCOPzJw5c3LfffflyiuvTIMGDdKxY8e8+OKLOeigg3LUUUdZ8I+l4oknnsipp56a1VdfPR06dMhLL72UqVOnZujQodl0000zffr0/PznP88WW2zhs5k655NPPslZZ52V119/PSussEJ69uyZnXbaKauttlo++uijnH766Zk2bVpuv/320gKs9evXz5gxY3LOOedkpZVWKp0BcP7552fkyJFZaaWVssUWW+SFF17I/Pnzc+WVV2b11Vf/znMT/7Xs3z/MoLZVVVVl+PDhueWWW7L77rtnt912S8OGDfPYY4/ltNNOy8EHH5wjjzwyZ511VkaOHJmHHnoozZo1K72PP/300wwcODBTpkwpLcCTJI888kjGjRuXqVOnZu+9907nzp1rcS/5bzR79uzssMMOOfXUUzN06NCss846ueiiizJv3rzceOONueKKKzJ8+PCsvfbapffzzJkzM2jQoHzyySf53e9+VzrdrrKyMh9++GHKysqyyiqr1PKe8d/iueeey3HHHZeLL744P/nJT0rvxwXv1+uuuy5XX311br/99tJvny+475xzzslf//rXPProo6Xne/nllzNmzJhMmTIlu+yyS9Zee+1a2S/++8yZMycnnHBCKioqMnDgwDRr1izz5s3L3nvvnS5duuTcc8/NnDlzcvPNN/tsps756KOPcvTRR2fevHn52c9+lldffTWjR49Os2bN8vvf/z7LL798rrvuutx9993p379/9t5771RVVZXex4MHD86zzz6bc845J+uvv37mzJmTMWPG5O677860adPSoUOHnHjiiaVfAPiufL1by4Q/y5J69erl0UcfzQYbbJCf/exnadWqVRo2bJhtttkm2267bR577LFUVVVlzz33TP369XPZZZcl+b/3caNGjdKhQ4dUV1dn+vTppdNPe/Tokf79++f0008X/ix1VVVVKSsry5prrplmzZrll7/8ZR588ME888wzadCgQXbYYYd06tQp5513XubMmVN6Pzdp0iQNGjRIRUVFjetGKyoq0rZtW/9xyVL1j3/8I506dcqWW26Z8vLyTJw4MbfcckuefPLJfPrpp9lzzz3Trl27XHXVVZk1a1aNx3bu3DnV1dV55513Stu6dOmSffbZJ8cee6zwZ6n66KOP8vjjj2fTTTdNixYtUl1dneWWWy6rr7566fr/5ZdfPj169Mjaa6/ts5k65esWxp4+fXppYexevXrVWBi7vLy8tODqzjvvXFrkL/nyfb3RRhvlggsuyNChQ3Pqqad+7/BPxD/wb8rKynLmmWdm0KBBWWmllZJ8eeSoadOmqaioSP369VNeXp6OHTvmoIMOyt13351HHnmkxsr9c+bMySeffFJajARqW3l5eaqrqzN27NisvPLK2WOPPbLBBhtkyJAhmTNnTjp06JCjjjoqr7zySq688sp88sknSb48W2D69Olp1apVysrKvJ+pVVVVVZk0aVKS5NJLL82OO+6Ym266Kf37988vf/nLvPjiizn++OPz+OOPZ/jw4Zk9e3bpPfvhhx/m888/T9OmTWtzFyDJl4uV/fjHP86//vWvJF+uG/Txxx9nypQp+clPflL6bfMOHTrkt7/9bV566aVcddVVPptZZv37wntftzD2CSecsEgLYzdr1qzGwtj/finWgsu8foj6P/gZgEJZsFr0gmv5F5xmV11dXfqmsWHDhvn5z3+el19+OWeddVY++OCD7Ljjjvnoo4/yyiuvpFevXqUPMahtC675X3XVVTNr1qxUVFTksMMOy29+85uMGDEin3/+edq2bZuzzjor559/fv7+979n1113zRtvvJH//d//zcEHH1zbuwBZY401svzyy+d3v/tdXnjhhVxxxRVZd911M2fOnBx33HH53e9+l7POOiv9+/fPVVddlQ8++CD77rtvZs2aleeffz477LBDmjRpUtu7AWnWrFl69uyZ66+/Ph9++GHatWuX++67L82aNcu1116bKVOmZP/998+BBx6Y7bbbLscee2yuueaaPPfccz6bWWZ8/vnnufvuu7PvvvuWfoavvLw8DRo0KC2M3bp165SVlZXOMlywMPaOO+6YvfbaK++8807uuOOOrLnmmtlll13SsGHDJb4wtmv+gUXSu3fvdOvWLSeeeGLpi4HKysqcdNJJ+etf/5q2bdvmo48+Stu2bXPZZZd9r0VIYEn57LPPsu222+aGG27IBhtskCTp169f/vGPf2TFFVfMoEGDsvvuu+eRRx7J8OHDM2PGjFRVVeXEE0/0yxQsEz7++OPsu++++eKLL9K2bdsMGzYsK664YsrKyvLEE0/kwgsvzHbbbZeTTjopl19+eX7/+9+ncePG+fTTT7PKKqvkkksuSceOHWt7NyDJlyucP/PMM3n33XczdOjQ7LPPPunbt28aNGiQP/7xj7nyyitz7LHHlhbE9tnMsuThhx/OWWedlQ8//DDHH398Dj744NKifTNnzswOO+yQfv36LZMLY4t/4FtNmzYtO+20Uy6++OL06NEjSUrfcH7++ed56623Mn78+DRs2LB0PyxLPvjggwwYMCDnnXdeqqurc8455+Tvf/976tevn9133720QvoCH3/8cenSF1hWXHvttbn00kvTp0+fnHPOOZk7d27pjKyDDz44lZWVueWWWzJ//vy8++67mThxYqqrq7PtttvW8szhq9166625++67c++996aioqIUOgvWsLj88strnPbss5na9re//S2nnXZaVl555dSrVy+ffPJJbrjhhqyyyiqZP39+5s2bl4suuigjRoxYJhfGds0/8K3GjRuXuXPnpn379km+PI36zTffzPTp07PCCitk/fXXz+677y78WWaVlZXljTfeyJlnnpndd989SXL99dfnwAMPzAMPPFC6xnQB/3HJsuhXv/pVVllllTz44IOZPHlyGjRoULrWtGnTpvnoo4+SfHmNaPv27UuLtcKyasKECWnbtm3Ky8tTVlZWWkOodevWmTBhwkI/Pemzmdo2b968tGnTJmeffXYOPPDAVFVV5Zprrkny5WdvRUVFevfuvcwujC3+ga+14D8qX3nllbRo0SIdOnTI5MmTc8EFF2TPPffMqFGjanmGsGhWXnnl9OjRIx999FHOO++8DBkyJFtttVX69OmTFVdcMc8//3xtTxG+VUVFRQYPHpzlllsuZ555ZqqqqvLFF1/k9ddfz2uvvZY+ffrU9hThO+nSpUvGjBmTZ599trTt5Zdfzrhx4/Kzn/2sFmcGX23zzTfPrbfemrXWWivdunXL5ptvnieeeCIvvPBCaUynTp1y4IEHLpMLYzvtH/hWhx9+eOrVq5ctt9wyV199dRo1apRBgwY5okSdMm3atNI39guOJlVXV+ejjz5KixYtanl2sOjuv//+nH322WnSpEk6duyYMWPGZJVVVsmll16adu3a1fb0YJFVV1end+/e+eKLL7LFFlukadOmGTlyZNZYY42ce+65adu2bW1PEb5SZWVlKioq8ve//z0XXHBBWrZsmWHDhpXu//TTT3PKKafkpZdeym9+85vSwtinn356Nt988xxzzDG1Mm/xD3yjOXPmpFevXpk0aVKWX375DBgwIAcddFBtTwvgv9qCo6Vvv/121ltvvfziF7+o7SnB9zJ16tRcccUVefvtt1NVVZXddtuttNAf1AVDhgzJPffck2OPPTZ77LFHafG/ZXFhbPEPfKsTTjghzZs3z3HHHbdYfmMUAODfzZgxI40aNUp5eXltTwUWyYKV/MePH5+zzz47s2bNym233ZYVVlihdGbAnDlz8uabb+btt99eJhbGFv/At/r3nykBAAD+zx133JHrrrsue+21V4488sg89thjadu2bTp16lQr1/Z/Hf81D3wr4Q8AADUtWBx71113zYYbbpjhw4fn0EMPzW9+85uMHj26lme3sPq1PQEAAACoa+rVq5fKysp8+umnmTdvXiZPnpxGjRrl2muvzTbbbFPb01uI0/4BAADgexg3blx++ctfZs6cOTnxxBOzzz771PaUvpb4BwAAgO+hsrIyN998c/r165fll1++tqfzjcQ/AAAAFJxVvAAAAKDgxD8AAAAUnPgHAACAghP/AAAAUHDiHwAAAApO/AMAAEDBiX8AAAAoOPEPAAAABSf+AYCl4vjjj8/666+ft99+e6H7rr322nTq1CmPP/54LcwMAIqvrLq6urq2JwEAFN9HH32Unj17Zu21184tt9xS2j5x4sT06tUr2267ba688spanCEAFJcj/wDAUrHyyivn+OOPz3PPPZc//OEPpe1nnXVW6tevn1NPPbUWZwcAxSb+AYClpk+fPtl4440zePDgfPzxx3nwwQfz1FNP5eijj07r1q1re3oAUFhO+wcAlqqxY8dmzz33TPfu3fPiiy+mdevWueeee1KvnmMSALCk1K/tCQAA/13WXHPNHHjggRk2bFjKy8szbNgw4Q8AS5i/aQGApW6llVZKkrRq1SprrrlmLc8GAIpP/AMAS9WUKVNy5ZVXZq211sqUKVNy/fXX1/aUAKDwxD8AsFSdffbZSZLrrrsuO++8c6655ppMnDixlmcFAMUm/gGApebhhx/OY489lqOOOipt2rTJKaeckgYNGuSss86q7akBQKGJfwBgqfj0009z7rnnpnPnzunXr1+SpHXr1jnqqKPy1FNP5c9//nMtzxAAikv8AwBLxeWXX55p06blrLPOSnl5eWn7fvvtl3XXXTfnn39+Pv3001qcIQAUl/gHAJa4f/3rX/n973+ffffdNxtssEGN+8rLy3PmmWfmww8/zOWXX147EwSAgiurrq6uru1JAAAAAEuOI/8AAABQcOIfAAAACk78AwAAQMGJfwAAACg48Q8AAAAFJ/4BAACg4MQ/AAAAFJz4BwAAgIIT/wAAAFBw4h8AAAAKTvwDAABAwYl/AAAAKLj/Dzc0d9UfzswgAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1040x715 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "geo_model_test = gp.create_geomodel(\n",
    "    project_name='Gempy_abc_Test',\n",
    "    extent=[0, 1000, -500, 500, -900, -700],\n",
    "    resolution=[100,100,100],\n",
    "    refinement=3,\n",
    "    structural_frame= gp.data.StructuralFrame.initialize_default_structure()\n",
    "    )\n",
    "p2d = gpv.plot_2d(geo_model_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([False, False, False, False, False])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "geo_model_test.grid.active_grids_bool"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'brk1' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[17], line 5\u001b[0m\n\u001b[1;32m      1\u001b[0m gp\u001b[38;5;241m.\u001b[39madd_surface_points(\n\u001b[1;32m      2\u001b[0m     geo_model\u001b[38;5;241m=\u001b[39mgeo_model_test,\n\u001b[1;32m      3\u001b[0m     x\u001b[38;5;241m=\u001b[39m[\u001b[38;5;241m100.0\u001b[39m, \u001b[38;5;241m900.0\u001b[39m],\n\u001b[1;32m      4\u001b[0m     y\u001b[38;5;241m=\u001b[39m[\u001b[38;5;241m0.0\u001b[39m, \u001b[38;5;241m0.0\u001b[39m],\n\u001b[0;32m----> 5\u001b[0m     z\u001b[38;5;241m=\u001b[39m[\u001b[43mbrk1\u001b[49m, brk1],\n\u001b[1;32m      6\u001b[0m     elements_names\u001b[38;5;241m=\u001b[39m[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124msurface1\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124msurface1\u001b[39m\u001b[38;5;124m'\u001b[39m]\n\u001b[1;32m      7\u001b[0m )\n\u001b[1;32m      8\u001b[0m gpv\u001b[38;5;241m.\u001b[39mplot_2d(geo_model_test, cell_number\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m11\u001b[39m)\n",
      "\u001b[0;31mNameError\u001b[0m: name 'brk1' is not defined"
     ]
    }
   ],
   "source": [
    "gp.add_surface_points(\n",
    "    geo_model=geo_model_test,\n",
    "    x=[100.0, 900.0],\n",
    "    y=[0.0, 0.0],\n",
    "    z=[brk1, brk1],\n",
    "    elements_names=['surface1', 'surface1']\n",
    ")\n",
    "gpv.plot_2d(geo_model_test, cell_number=11)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gp.add_orientations(\n",
    "    geo_model=geo_model_test,\n",
    "    x=[500],\n",
    "    y=[0.0],\n",
    "    z=[brk1],\n",
    "    elements_names=['surface1'],\n",
    "    pole_vector=[[0, 0, 1]]\n",
    ")\n",
    "\n",
    "gpv.plot_2d(geo_model_test, cell_number=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model_test.update_transform(gp.data.GlobalAnisotropy.NONE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gp.compute_model(geo_model_test, engine_config=gp.data.GemPyEngineConfig())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model_test.interpolation_options.kernel_options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gpv.plot_2d(geo_model_test, cell_number=[5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#gpv.plot_3d(geo_model_test, show_surfaces=True, image=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model_test.structural_frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "brk2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "element2 = gp.data.StructuralElement(\n",
    "    name='surface2',\n",
    "    color=next(geo_model_test.structural_frame.color_generator),\n",
    "    surface_points=gp.data.SurfacePointsTable.from_arrays(\n",
    "        x=np.array([100.0, 900.0]),\n",
    "        y=np.array([0.0, 0.0]),\n",
    "        z=np.array([brk2, brk2]),\n",
    "        names='surface2'\n",
    "    ),\n",
    "    orientations=gp.data.OrientationsTable.initialize_empty()\n",
    ")\n",
    "\n",
    "geo_model_test.structural_frame.structural_groups[0].append_element(element2)\n",
    "# Compute and visualize the updated model:\n",
    "gp.compute_model(geo_model_test)\n",
    "gpv.plot_2d(geo_model_test, cell_number=5, legend='force')\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#gpv.plot_3d(geo_model_test, image=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model_test.structural_frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model_test.structural_frame.structural_groups[0].elements[0], geo_model_test.structural_frame.structural_groups[0].elements[1] = \\\n",
    "geo_model_test.structural_frame.structural_groups[0].elements[1], geo_model_test.structural_frame.structural_groups[0].elements[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model_test.structural_frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gpv.plot_2d(geo_model_test, cell_number=5, legend='force')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create a custome grid where the observed data information is available "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_loc = 300\n",
    "y_loc = 0\n",
    "z_loc = dataframe_KSL['Z']\n",
    "xyz_coord = np.array([[x_loc, y_loc, z] for z in z_loc])\n",
    "gp.set_custom_grid(geo_model_test.grid, xyz_coord=xyz_coord)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp_coords_copy_test = geo_model_test.interpolation_input.surface_points.sp_coords.copy()\n",
    "geo_model_test.transform.apply_inverse(sp_coords_copy_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gp.add_surface_points(\n",
    "    geo_model=geo_model_test,\n",
    "    x=[x_loc, x_loc],\n",
    "    y=[0.0, 0.0],\n",
    "    z=[brk1, brk2],\n",
    "    elements_names=['surface1', 'surface2']\n",
    ")\n",
    "gpv.plot_2d(geo_model_test, cell_number=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model_test.surface_points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gp.compute_model(geo_model_test)\n",
    "gpv.plot_2d(geo_model_test, cell_number=5, legend='force')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp_coords_copy_test = geo_model_test.interpolation_input.surface_points.sp_coords.copy()\n",
    "geo_model_test.transform.apply_inverse(sp_coords_copy_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Change the backend to PyTorch for probabilistic modeling\n",
    "BackendTensor.change_backend_gempy(engine_backend=gp.data.AvailableBackends.PYTORCH)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model_test.interpolation_options.uni_degree = 0\n",
    "geo_model_test.interpolation_options.mesh_extraction = False\n",
    "geo_model_test.interpolation_options.sigmoid_slope = 1100."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#geo_model_test.solutions.octrees_output[0].last_output_center"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "custom_grid_values = geo_model_test.solutions.octrees_output[0].last_output_center.custom_grid_values\n",
    "custom_grid_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_obs_label.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model_test.surface_points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model_test.orientations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model_test.transform.apply_inverse(sp_coords_copy_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_obs_label = torch.tensor(y_obs_label.to_numpy(),dtype=torch.float64)\n",
    "y_obs_label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "@config_enumerate\n",
    "def model_test(y_obs_label):\n",
    "    \"\"\"\n",
    "    This Pyro model represents the probabilistic aspects of the geological model.\n",
    "    It defines a prior distribution for the top layer's location and\n",
    "    computes the thickness of the geological layer as an observed variable.\n",
    "    \"\"\"\n",
    "    # Define prior for the top layer's location\n",
    "    prior_mean_surface_1 = sp_coords_copy_test[2, 2]\n",
    "    prior_mean_surface_2 = sp_coords_copy_test[5, 2]\n",
    "    \n",
    "    \n",
    "\n",
    "    mu_surface_1 = pyro.sample('mu_1', dist.Normal(prior_mean_surface_1, torch.tensor(0.02, dtype=torch.float64)))\n",
    "    mu_surface_2 = pyro.sample('mu_2', dist.Normal(prior_mean_surface_2, torch.tensor(0.02, dtype=torch.float64)))    \n",
    "    \n",
    "    \n",
    "    # Ensure that mu_surface_1 is greater than mu_surface_2\n",
    "    pyro.sample('mu_1 > mu_2', dist.Delta(torch.tensor(1.0, dtype=torch.float64)), obs=(mu_surface_1 > mu_surface_2))\n",
    "    # Update the model with the new top layer's location\n",
    "    interpolation_input = geo_model_test.interpolation_input\n",
    "    \n",
    "    \n",
    "    interpolation_input.surface_points.sp_coords = torch.index_put(\n",
    "        interpolation_input.surface_points.sp_coords,\n",
    "        (torch.tensor([2]), torch.tensor([2])),\n",
    "        mu_surface_1\n",
    "    )\n",
    "    interpolation_input.surface_points.sp_coords = torch.index_put(\n",
    "        interpolation_input.surface_points.sp_coords,\n",
    "        (torch.tensor([5]), torch.tensor([2])),\n",
    "        mu_surface_2\n",
    "    )\n",
    "    \n",
    "    \n",
    "    # # Compute the geological model\n",
    "    geo_model_test.solutions = gempy_engine.compute_model(\n",
    "        interpolation_input=interpolation_input,\n",
    "        options=geo_model_test.interpolation_options,\n",
    "        data_descriptor=geo_model_test.input_data_descriptor,\n",
    "        geophysics_input=geo_model_test.geophysics_input,\n",
    "    )\n",
    "    \n",
    "    # Compute and observe the thickness of the geological layer\n",
    "    \n",
    "    custom_grid_values = geo_model_test.solutions.octrees_output[0].last_output_center.custom_grid_values\n",
    "    \n",
    "    lambda_ = 20\n",
    "    # #class_label = F.softmax(-lambda_* (torch.tensor([1,2,3], dtype=torch.float64) - custom_grid_values.reshape(-1,1))**2, dim=1)\n",
    "    z_nk = F.softmax(-lambda_* (torch.tensor([1,2,3], dtype=torch.float64) - custom_grid_values.reshape(-1,1))**2, dim=1)\n",
    "    # #class_label = torch.mean(F.softmax(-lambda_* (torch.tensor([1,2,3,4,5,6], dtype=torch.float64) - custom_grid_values.reshape(-1,1))**2, dim=1),dim=0)\n",
    "    \n",
    "    N_k = torch.sum(z_nk,axis=0)\n",
    "    N = len(custom_grid_values)\n",
    "    pi_k = N_k /N\n",
    "    #print(pi_k)\n",
    "    mean = []\n",
    "    cov = []\n",
    "    for i in range(z_nk.shape[1]):\n",
    "        mean_k = torch.sum( z_nk[:,i][:,None] * y_obs_label, axis=0)/ N_k[i]\n",
    "        #cov_k = torch.sum( (normalised_hsi - mean_k.reshape((-1,1))) (normalised_hsi - mean_k).T )\n",
    "        cov_k = torch.zeros((mean_k.shape[0],mean_k.shape[0]),dtype=torch.float64)\n",
    "        for j in range(z_nk.shape[0]):\n",
    "                cov_k +=  z_nk[j,i]* torch.matmul((y_obs_label[j,:] - mean_k).reshape((-1,1)) ,(y_obs_label[j,:] - mean_k).reshape((1,-1)))\n",
    "        mean.append(mean_k)\n",
    "        cov_k=cov_k/N_k[i] + 1e-5 * torch.diag(torch.ones(cov_k.shape[0],dtype=torch.float64))\n",
    "        cov.append(cov_k)\n",
    "    mean_tensor = torch.stack(mean, dim=0)\n",
    "    cov_tensor = torch.stack(cov,dim=0)\n",
    "    \n",
    "    \n",
    "    with pyro.plate('N='+str(y_obs_label.shape[0]), y_obs_label.shape[0]):\n",
    "        assignment = pyro.sample(\"assignment\", dist.Categorical(pi_k))\n",
    "        \n",
    "        obs = pyro.sample(\"obs\", dist.MultivariateNormal(loc=mean_tensor[assignment],covariance_matrix=cov_tensor[assignment]), obs=y_obs_label)\n",
    "        \n",
    "    return obs\n",
    "dependencies = get_dependencies(model_test, model_args=(y_obs_label,))\n",
    "pyro.render_model(model_test, model_args=(y_obs_label,),render_distributions=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print( y_obs_label.shape)\n",
    "model_test(y_obs_label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior = Predictive(model_test, num_samples=100)(y_obs_label)\n",
    "\n",
    "# Key to avoid\n",
    "avoid_key = 'mu_1 > mu_2'\n",
    "\n",
    "# Create sub-dictionary without the avoid_key\n",
    "prior = dict((key, value) for key, value in prior.items() if key != avoid_key)\n",
    "\n",
    "data = az.from_pyro(prior=prior)\n",
    "az.plot_trace(data.prior)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pyro.primitives.enable_validation(is_validate=True)\n",
    "nuts_kernel = NUTS(model_test, step_size=0.0085, adapt_step_size=True, target_accept_prob=0.9, max_tree_depth=10, init_strategy=init_to_mean)\n",
    "#nuts_kernel = NUTS(model_test, step_size=0.00085, adapt_step_size=True, target_accept_prob=0.9, max_tree_depth=10)\n",
    "#nuts_kernel = NUTS(model_test)\n",
    "initial_values = {'mu_1': torch.tensor(0.01, dtype=torch.float64),'mu_2': torch.tensor(0.01, dtype=torch.float64) }\n",
    "#mcmc = MCMC(nuts_kernel, num_samples=200, warmup_steps=50, disable_validation=False, initial_params=initial_values)\n",
    "mcmc = MCMC(nuts_kernel, num_samples=200, warmup_steps=50, disable_validation=False)\n",
    "mcmc.run(y_obs_label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "posterior_samples = mcmc.get_samples()\n",
    "posterior_predictive = Predictive(model_test, posterior_samples)(y_obs_label)\n",
    "data = az.from_pyro(posterior=mcmc, prior=prior, posterior_predictive=posterior_predictive)\n",
    "az.plot_trace(data)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "posterior_samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loc_mean_posterior, loc_cov_posterior =[], []\n",
    "for key, values in posterior_samples.items():\n",
    "    print(key)\n",
    "    if key == \"sample_data\":\n",
    "        mean = values.mean(dim=0)\n",
    "        cov = values.std(dim=0)\n",
    "        print(\"mean\\n\",mean)\n",
    "        print(\"cov\\n\", cov)\n",
    "        loc_mean_posterior.append(mean.detach().numpy())\n",
    "        loc_cov_posterior.append(cov.detach().numpy())\n",
    "    elif key == \"sigma_data\":\n",
    "        print(\"mean\\n\",values.mean(dim=0), \"\\nstd\\n\", values.std(dim=0))\n",
    "    else:\n",
    "        print(\"mean\\n\",values.mean(), \"\\nstd\\n\", values.std())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loc_mean_posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loc_cov_posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 6))\n",
    "x = np.linspace(-0.5,1.5,8000)\n",
    "  # Combine x and y grids into a (100, 100, 2) array\n",
    "for i in range(3):\n",
    "    # Create a multivariate normal distribution\n",
    "    rv = norm(loc_mean_posterior[0][i], loc_cov_posterior[0][i])\n",
    "\n",
    "    # Calculate PDF values for each point in the grid\n",
    "    pdf_values = rv.pdf(x)\n",
    "\n",
    "    # Plot the Gaussian distribution using contour plot\n",
    "    \n",
    "    plt.scatter(x, pdf_values, s=2.5, label='Gaussian_'+str(i+1))\n",
    "#plt.colorbar(label='Probability Density')\n",
    "plt.xlabel('hsi_spectra')\n",
    "plt.ylabel('Unnormalized Probability Density')\n",
    "plt.title('Gaussian Distribution')\n",
    "#plt.scatter(loc_mean[i], color='red', label='Mean')\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "az.plot_density(\n",
    "    data=[data.posterior, data.prior],\n",
    "    shade=.9,\n",
    "    var_names=['mu_1'],\n",
    "    data_labels=[\"Posterior Predictive\", \"Prior Predictive\"],\n",
    "    colors=[default_red, default_blue],\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "az.plot_density(\n",
    "    data=[data.posterior, data.prior],\n",
    "    shade=.9,\n",
    "    var_names=['mu_2'],\n",
    "    data_labels=[\"Posterior Predictive\", \"Prior Predictive\"],\n",
    "    colors=[default_red, default_blue],\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Update the model with the new top layer's location\n",
    "interpolation_input = geo_model_test.interpolation_input\n",
    "interpolation_input.surface_points.sp_coords = torch.index_put(\n",
    "    interpolation_input.surface_points.sp_coords,\n",
    "    (torch.tensor([2]), torch.tensor([2])),\n",
    "    posterior_samples[\"mu_1\"].mean()\n",
    ")\n",
    "interpolation_input.surface_points.sp_coords = torch.index_put(\n",
    "    interpolation_input.surface_points.sp_coords,\n",
    "    (torch.tensor([5]), torch.tensor([2])),\n",
    "    posterior_samples[\"mu_2\"].mean()\n",
    ")\n",
    "\n",
    "#print(\"interpolation_input\",interpolation_input.surface_points.sp_coords)\n",
    "\n",
    "# # Compute the geological model\n",
    "geo_model_test.solutions = gempy_engine.compute_model(\n",
    "    interpolation_input=interpolation_input,\n",
    "    options=geo_model_test.interpolation_options,\n",
    "    data_descriptor=geo_model_test.input_data_descriptor,\n",
    "    geophysics_input=geo_model_test.geophysics_input,\n",
    ")\n",
    "\n",
    "gpv.plot_2d(geo_model_test, cell_number=5, legend='force')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp_coords_copy_test2 =interpolation_input.surface_points.sp_coords\n",
    "sp_coords_copy_test2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp_cord= geo_model_test.transform.apply_inverse(sp_coords_copy_test2.detach().numpy())\n",
    "sp_cord"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model_post = gp.create_geomodel(\n",
    "    project_name='Gempy_post_Test',\n",
    "    extent=[0, 1000, -10, 10, -900, -400],\n",
    "    resolution=[100,10,100],\n",
    "    refinement=7,\n",
    "    structural_frame= gp.data.StructuralFrame.initialize_default_structure()\n",
    "    )\n",
    "gp.add_surface_points(\n",
    "    geo_model=geo_model_post,\n",
    "    x=sp_cord[3:,0],\n",
    "    y=sp_cord[3:,1],\n",
    "    z=sp_cord[3:,2],\n",
    "    elements_names=['surface1', 'surface1','surface1']\n",
    ")\n",
    "gp.add_orientations(\n",
    "    geo_model=geo_model_post,\n",
    "    x=[500],\n",
    "    y=[0.0],\n",
    "    z=[brk1],\n",
    "    elements_names=['surface1'],\n",
    "    pole_vector=[[0, 0, 1.0]]\n",
    ")\n",
    "gpv.plot_2d(geo_model_post, cell_number=5, legend='force')\n",
    "geo_model_test.update_transform(gp.data.GlobalAnisotropy.NONE)  \n",
    "gp.compute_model(geo_model_post, engine_config=gp.data.GemPyEngineConfig())\n",
    "gpv.plot_2d(geo_model_post, cell_number=[5])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "element2 = gp.data.StructuralElement(\n",
    "    name='surface2',\n",
    "    color=next(geo_model_post.structural_frame.color_generator),\n",
    "    surface_points=gp.data.SurfacePointsTable.from_arrays(\n",
    "        x=sp_cord[:3,0],\n",
    "        y=sp_cord[:3,1],\n",
    "        z=sp_cord[:3,2],\n",
    "        names='surface2'\n",
    "    ),\n",
    "    orientations=gp.data.OrientationsTable.initialize_empty()\n",
    ")\n",
    "geo_model_post.structural_frame.structural_groups[0].append_element(element2)\n",
    "gp.compute_model(geo_model_post)\n",
    "#gpv.plot_2d(geo_model_post, cell_number=5, legend='force')\n",
    "geo_model_post.structural_frame.structural_groups[0].elements[0], geo_model_post.structural_frame.structural_groups[0].elements[1] = \\\n",
    "geo_model_post.structural_frame.structural_groups[0].elements[1], geo_model_post.structural_frame.structural_groups[0].elements[0]\n",
    "gpv.plot_2d(geo_model_post, cell_number=5, legend='force')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model_post.surface_points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Creating marginal distribution plots\n",
    "# p = PlotPosterior(data)\n",
    "# p.create_figure(figsize=(9, 5), joyplot=False, marginal=True, likelihood=False)\n",
    "# p.plot_marginal(\n",
    "#     var_names=['mu_1', 'mu_2'],\n",
    "#     plot_trace=False,\n",
    "#     credible_interval=.70,\n",
    "#     kind='kde',\n",
    "#     marginal_kwargs={\"bw\": 1}\n",
    "# )\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Visualizing the posterior distributions\n",
    "# p = PlotPosterior(data)\n",
    "# p.create_figure(figsize=(9, 6), joyplot=True)\n",
    "# iteration = 99\n",
    "# p.plot_posterior(\n",
    "#     prior_var=['mu_1', 'mu_2'],\n",
    "#     like_var=['mu_1', 'mu_2'],\n",
    "#     obs='obs',\n",
    "#     iteration=iteration,\n",
    "#     marginal_kwargs={\n",
    "#         \"credible_interval\": 0.94,\n",
    "#         'marginal_kwargs': {\"bw\": 1},\n",
    "#         'joint_kwargs': {\"bw\": 1}\n",
    "#     }\n",
    "# )\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Creating a pair plot for selected parameters\n",
    "# az.plot_pair(data, divergences=False, var_names=['mu_1', 'mu_2'])\n",
    "# plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gempy_pyro",
   "language": "python",
   "name": "python3"
  },
  "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}