eval_SBI.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SBI eval "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 100k sim - 1 hr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "from petitRADTRANS import Radtrans\n",
    "import petitRADTRANS.retrieval_examples.emission.master_retrieval_model as rm\n",
    "from petitRADTRANS import nat_cst as nc\n",
    "import petitRADTRANS.rebin_give_width as rgw\n",
    "from scipy.interpolate import interp1d\n",
    "import sklearn\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "import pandas as pd\n",
    "\n",
    "from torch.distributions import Independent, Distribution\n",
    "from fab_priors import BoxUniform_New\n",
    "\n",
    "import torch\n",
    "from sbi.inference import SNRE_A, SNRE, prepare_for_sbi, simulate_for_sbi\n",
    "from sbi.utils.get_nn_models import posterior_nn\n",
    "from sbi import utils as utils\n",
    "from sbi.types import Array, OneOrMore, ScalarFloat\n",
    "\n",
    "#from Simulator import Simulator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "op= '/home/mvasist/petitRADTRANS/petitRADTRANS/retrieval_examples/emission/'\n",
    "observation_files = {}\n",
    "observation_files['NIRISS SOSS'] = op +'NIRISS_SOSS_flux.dat'\n",
    "observation_files['NIRSpec G395M'] = op +'NIRSpec_G395M_flux.dat'\n",
    "observation_files['MIRI LRS'] = op +'MIRI_LRS_flux.dat'\n",
    "\n",
    "# Wavelength range of observations, fixed parameters that will not be retrieved\n",
    "WLENGTH = [0.3, 15.0]\n",
    "# LOG_G =  2.58\n",
    "R_pl =   1.84*nc.r_jup_mean\n",
    "R_star = 1.81*nc.r_sun\n",
    "gamma = 1\n",
    "t_equ= 0\n",
    "\n",
    "# Get host star spectrum to calculate F_pl / F_star later.\n",
    "T_star = 6295.\n",
    "x = nc.get_PHOENIX_spec(T_star)\n",
    "fstar = interp1d(x[:,0], x[:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NIRISS SOSS\n",
      "NIRSpec G395M\n",
      "MIRI LRS\n"
     ]
    }
   ],
   "source": [
    "####################################################################################\n",
    "####################################################################################\n",
    "### READ IN OBSERVATION\n",
    "####################################################################################\n",
    "####################################################################################\n",
    "\n",
    "# Read in data, convert all to cgs! \n",
    "\n",
    "'''\n",
    "Im using only data_flux_nu_error['MIRI LRS'] from here to calculate the likelihood. \n",
    "'''\n",
    "\n",
    "data_wlen = {}\n",
    "data_flux_nu = {}\n",
    "data_flux_nu_error = {}\n",
    "data_wlen_bins = {}\n",
    "\n",
    "for name in observation_files.keys():\n",
    "    print(name)\n",
    "    dat_obs = np.genfromtxt(observation_files[name])\n",
    "    data_wlen[name] = dat_obs[:,0]*1e-4\n",
    "    data_flux_nu[name] = dat_obs[:,1]\n",
    "    data_flux_nu_error[name] = dat_obs[:,2]\n",
    "    \n",
    "    data_wlen_bins[name] = np.zeros_like(data_wlen[name])\n",
    "    data_wlen_bins[name][:-1] = np.diff(data_wlen[name])\n",
    "    data_wlen_bins[name][-1] = data_wlen_bins[name][-2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Simulator(params): \n",
    "\n",
    "    NaN_spectra = 0\n",
    "\n",
    "    atmosphere = Radtrans(line_species = ['H2O', 'CO_all_iso', \\\n",
    "                                         'CO2', 'CH4', \\\n",
    "                                          'Na', 'K'], \\\n",
    "          rayleigh_species = ['H2', 'He'], \\\n",
    "          continuum_opacities = ['H2-H2', 'H2-He'], \\\n",
    "          wlen_bords_micron = WLENGTH)#, mode='c-k')\n",
    "\n",
    "\n",
    "    pressures = np.logspace(-6, 2, 100)\n",
    "    atmosphere.setup_opa_structure(pressures)\n",
    "    temperature = 1200. * np.ones_like(pressures)\n",
    "\n",
    "    \n",
    "    t_int = params[0].numpy()                             #200.\n",
    "    log_kappa_IR = params[1].numpy()                      #-2\n",
    "    log_gravity = params[2].numpy()                       #params[5].numpy() 1e1**2.45 \n",
    "\n",
    "    gravity = np.exp(log_gravity)\n",
    "    kappa_IR = np.exp(log_kappa_IR)\n",
    "    \n",
    "    temperature = nc.guillot_global(pressures, kappa_IR, gamma, gravity, t_int, t_equ)\n",
    "    \n",
    "    abundances = {}\n",
    "    abundances['H2'] = 0.75 * np.ones_like(temperature) #0.74 * np.ones_like(temperature) (params[3].numpy())\n",
    "    abundances['He'] = 0.25 * np.ones_like(temperature)  #0.24 * np.ones_like(temperature) (params[4].numpy())\n",
    "    abundances['H2O'] = 0.001 * np.ones_like(temperature)\n",
    "    abundances['CO_all_iso'] = 0.01 * np.ones_like(temperature)\n",
    "    abundances['CO2'] = 0.00001 * np.ones_like(temperature)\n",
    "    abundances['CH4'] = 0.000001 * np.ones_like(temperature)\n",
    "    abundances['Na'] = 0.00001 * np.ones_like(temperature)\n",
    "    abundances['K'] = 0.000001 * np.ones_like(temperature)\n",
    "\n",
    "    MMW = rm.calc_MMW(abundances) * np.ones_like(temperature)\n",
    "    #print(MMW, abundances)\n",
    "    \n",
    "    atmosphere.calc_flux(temperature, abundances, gravity, MMW)\n",
    "\n",
    "    wlen, flux_nu = nc.c/atmosphere.freq, atmosphere.flux/1e-6\n",
    "\n",
    "\n",
    "    # Just to make sure that a long chain does not die\n",
    "    # unexpectedly:\n",
    "    # Return -inf if forward model returns NaN values\n",
    "    if np.sum(np.isnan(flux_nu)) > 0:\n",
    "        print(\"NaN spectrum encountered\")\n",
    "        NaN_spectra += 1\n",
    "        return torch.ones([1,371])* -np.inf\n",
    "\n",
    "    # Convert to observation for emission case\n",
    "    flux_star = fstar(wlen)\n",
    "    flux_sq   = flux_nu/flux_star*(R_pl/R_star)**2 \n",
    "\n",
    "    flux_rebinned = rgw.rebin_give_width(wlen, flux_sq, \\\n",
    "                data_wlen['MIRI LRS'], data_wlen_bins['MIRI LRS'])\n",
    "\n",
    "    #flux_rebinned = np.reshape(flux_rebinned, (371,1))    \n",
    "\n",
    "    FR= torch.Tensor(flux_rebinned)\n",
    "    \n",
    "    return FR    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prior= utils.BoxUniform(low =torch.tensor([0., 0, 2.]), \\\n",
    "#                              high =torch.tensor([2., 4000, 3.7]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "Prior= utils.BoxUniform(low=torch.tensor([0., -4 , 2 ]), high=torch.tensor([2000., 0, 3.7 ]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "  Read CIA opacities for H2-H2...\n",
      "  Read CIA opacities for H2-He...\n",
      " Done.\n",
      "\n",
      "\n",
      "  Read CIA opacities for H2-H2...\n",
      "  Read CIA opacities for H2-He...\n",
      " Done.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "simulator, prior = prepare_for_sbi(Simulator, Prior)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "inference = SNRE_A(prior= Prior, device= 'cpu', classifier='resnet')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#reading the csv file\n",
    "\n",
    "# Xpd= pd.read_csv('X_1M.csv')\n",
    "# Tpd= pd.read_csv('T_1M.csv')\n",
    "\n",
    "X=[]\n",
    "T=[]\n",
    "\n",
    "# for i in range(1,11):\n",
    "# #     print(i)\n",
    "#     dfX= pd.read_csv('/home/mvasist/simulations/3_params/X_'+ str(i) + '.csv')\n",
    "#     dfT= pd.read_csv('/home/mvasist/simulations/3_params/T_'+ str(i) + '.csv')\n",
    "#     X.append(dfX.values)\n",
    "#     T.append(dfT.values)\n",
    "    \n",
    "# for j in range(1,51):\n",
    "#     dfX_= pd.read_csv('/home/mvasist/simulations/3_params/_X_'+ str(i) + '.csv')\n",
    "#     dfT_= pd.read_csv('/home/mvasist/simulations/3_params/_T_'+ str(i) + '.csv')\n",
    "#     X.append(dfX_.values)\n",
    "#     T.append(dfT_.values)\n",
    "    \n",
    "# for k in range(1, 101):\n",
    "#     dfX_M= pd.read_csv('/home/mvasist/simulations/3_params/X_1Msim_'+ str(k) + '.csv', low_memory = False) #chunksize=1000, iterator=True, dtype={'col2': np.float32}\n",
    "#     dfT_M= pd.read_csv('/home/mvasist/simulations/3_params/T_1Msim_'+ str(k) + '.csv', low_memory = False)\n",
    "#     X.append(dfX_M.values) #for x,t in zip(dfX_M,dfT_M):\n",
    "#     T.append(dfT_M.values)\n",
    "\n",
    "\n",
    "#log gamma, Tequ and log gravity \n",
    "for k in range(1, 101):\n",
    "    if k==15: continue\n",
    "    dfX_M= pd.read_csv('/home/mvasist/simulations/3_params/1/X_100ksim_TintLkIRLg'+ str(k) + '.csv', low_memory = False) #chunksize=1000, iterator=True, dtype={'col2': np.float32}\n",
    "    dfT_M= pd.read_csv('/home/mvasist/simulations/3_params/1/T_100ksim_TintLkIRLg'+ str(k) + '.csv', low_memory = False)\n",
    "    X.append(dfX_M.values) #for x,t in zip(dfX_M,dfT_M):\n",
    "    T.append(dfT_M.values)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "99"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data = pd.read_csv(filepath, chunksize=1000, iterator=True) \n",
    "#Loop through the chunks and process the data\n",
    "    \n",
    "comb_np_array_X = np.vstack(X)\n",
    "comb_np_array_T = np.vstack(T)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xframe = pd.DataFrame(comb_np_array_X)\n",
    "Tframe = pd.DataFrame(comb_np_array_T)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "list_of_tensors_X = [torch.tensor(np.array(Xframe),dtype=torch.float32)]\n",
    "list_of_tensors_T = [torch.tensor(np.array(Tframe),dtype=torch.float32)] #torch.tensor\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "XX = torch.cat(list_of_tensors_X)[:, 1:]\n",
    "TT = torch.cat(list_of_tensors_T)[:,1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([98901, 371])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(XX)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([98901, 3])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(TT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "inference = inference.append_simulations(TT, XX)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 100k sim - epochs - 20 min"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Neural network successfully converged after 55 epochs.\n"
     ]
    }
   ],
   "source": [
    "density_estimator = inference.train()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "posterior = inference.build_posterior(density_estimator)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Posterior conditional density p(θ|x) (amortized).\n",
      "\n",
      "This RatioBasedPosterior-object was obtained with a SNRE_A-class method using a residualnet.\n",
      "It provides MCMC to .sample() from the posterior and can evaluate the posterior density with .log_prob().\n"
     ]
    }
   ],
   "source": [
    "print(posterior)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "observation = torch.load('3param_observation_TintkIRLg.pt') #log_gamma, T_int, T_equ = 1.5, 750, 2000\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 100k samples is around 5 hours for 100k sim\n",
    "# 100k samples is around 6 hours for 200k sim\n",
    "\n",
    "#TintLkIRLg 100k sim: 10k sampl-23 min, 100ksampl-\n",
    "\n",
    "#100k gmTequG- 5:41 hrs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Tuning bracket width...: 100%|██████████| 50/50 [01:41<00:00,  2.02s/it]  \n",
      "Generating samples: 100%|██████████| 20/20 [00:03<00:00,  6.49it/s]\n",
      "Generating samples:  40%|████      | 40250/100000 [1:37:53<2:19:28,  7.14it/s]"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "sampls= 100000\n",
    "\n",
    "samples = posterior.sample((sampls,), x=observation)\n",
    "log_probability = posterior.log_prob(samples, x= observation)\n",
    "end= time.time()\n",
    "print('it takes: '+ str((end-start)/60) + ' min')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([10000, 3])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "samples.size()\n",
    "#log_probability.size()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = torch.cat([torch.ones(1)*750.,torch.ones(1)*-2.,torch.ones(1)*2.58])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAKmCAYAAACixf+8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABOe0lEQVR4nO3deZxcVZn/8e/TS9KdhKQDYe0AEQVZFAFR0SBgFEVxRHEBR0R0Rsdxxfm5RJ1RXMaJOo4bg7uDC4qiMDKgKBIQjUaWGNk3ISxNIGTprJ3e6vn9cW/VXVLVXd1d3VXV5/N+verVp26de+tU9fbUec5i7i4AAIBQtNS7AQAAAFOJ4AcAAASF4AcAAASF4AcAAASF4AcAAASF4AcAAASF4AeAzOyVZna9ma0zsz4ze9DM/tfMTknVOcnM3MyGzOyQMtd4xMwuTN1fFNcv3obN7DEzu8jM9i9zvo9we+VkvXYA4WmrdwMA1JeZvUfSlyV9V9LnJW2X9GRJp0paIumq3Cmtkj4p6cwqn+I/JF0uaYak4yR9XNJhZvYcdx/M1b1Q0jfKXOPuKp8LAEZF8APg/ZL+193/IXVsuaRvmVm53uHfSHqdmf2Hu/+1iuvf7+4r4/L1ZtYu6dOSnilpZa5uT6ouAEwK0l4Adpf0WLkH3L1Q5vD5ktYqCmDGY1X89YBxng8AE0LwA+AGSW8ysw+UG8tTRp+iwOflZnbcOJ5vUfz1b2UeMzNry9/G8RwAUBHBD4C3S7pP0uck3W1m683sx2b24hHO+bak+yX9exXXb4mDmFlmtkTSRyX93N1vLlP3I5IG8zczWzCG1wMAIyL4AQLn7vdIOlrSiYqCmdWSXiXp12b2rxXOGZR0nqQlZvaiUZ7iG4qCmO2SrpH0uKSzKtT9rqRnlbn1Vvt6AGA0BD8A5O7D7n69u/+ru79I0kGSbpX0cTObX+G0iyTdrtHH/nxaUQBzoqLxQsdIuqBC3bXuflOZ29CYXxQAVEDwA2AX7v6ootRWm6SDK9QpSPo3Sc8xs9NGuNyDcQBzvbu/W9L3JL3ZzJ5d63YDQDUIfoDAmdm+FR46NP5adiaYJLn7ZZJulPQpSVblUy5VNGj649W2EQBqiVkUAG4zs99K+qWkByTNlfQyRQOhf+ruD41y/kcVrf1TFXd/zMz+W9L7zeyZuYHP3RVmkD3o7murfQ4AGAk9PwA+KqlT0arNv5H0E0nPVdRD88bRTnb3qyVdN8bn/KykrZI+ljt+jqQ/lbm9YYzXB4CKzN3r3QYAAIApQ88PAAAICsEPAAAICsEPAAAICsEPAAAICsEPAAAICsFPkzGzk8zMc7fHzOyX+RVzzey8+PFd1nMys0XxY/84da1vXKn39aQ6Pf/xZnahmd1mZkNmtqZCvXNS3/dddmA3sxNTj4+259Zobarbe2Jma8zswnGc90oz+5cyx+v6/R0rM/uMmf3GzDbE7T5nhLpvNbO7zKzfzO42s7eXqbPGzH44qY1uAKnfj0X1bgsaG8FP83qPorVYnifp3ZL2kvRbM3tSXVvVvFYpej9X1en5Xyjp+Yr2yrqzivpbVX4NnjfFj9VCvd+T8XilpF2CHzXfa3m3orWXrhipkpm9VdHGsT+XdIqkSyRdYGb/POktBJoYwU/zutPdV7r7n9z9EklnStpN0R9AVMnMWs2szd23xO/nllpdc4ynfcrdn+zuZ0j6axX1L5V0lpmVtpQws05Jr1H0j3DCavme1FsTvpZ57v58RduGlBX/jP27pB+4+0fd/Vp3/1dJF0r6lJm1T01TgeZD8DN9FP+oT9ofvDgt84iZHW1mvzezHWZ2b4Vu9ieZ2UVm9kTcHb/azF6VevyZcff08alj746PfTp17OD42KkjtKuY0nh13MZNZrYlfv49cnXdzP7dzJaa2QOSBiQ9vVxaxCLvi1MJA2a21szON7O51VxzLO9tvEnoWPxA0oGSjk8de5Wi3+mqgh8zO8TMLjOzdWa208weMrNLioFbhffkOjP7g5mdEn9P+8zsL2b2HDNri9M1a81sY/y9mJ06t2zqqZpUhZntaWbfMLN74p+7h83sR2bWnapzoaKer25LUn9rRngtY/n+ftrM3mNmD5jZVjP7nZkdkav3EjP7o5ltNrNt8XXzK1hXpcqfh+dK2lNSPp31A0l7KPuzkWFRgP7N+PekYnrUzPYxs++Z2aPx7/FaM7vCzPZK1fmEma2Kr7XezJZbbouS1Pv/yvj7uNHMes3sS3FbnhX/XG03s9vN7CW584t/e55nZjfGP69rzOzdVbxPMrO3mdlf4/PWm9l3zGz3XJ33mtmd8c/0JjO7yVJ/szC9sLdX82qJ/0mZpIWSPiNph6T/m+TnnSvpR5K+pGg7hDdL+pqZ3e3u10qSme0v6c+S1kl6n6QnJJ0h6edm9kp3v1zSXyT1Sloi6Q/xtZco2vBySer5lkgaknR9FW37kqTfSnq9op3IPyNpP0kvyNU7R9L9kt4vabukRyXNK3O9f5f0YUn/reh9PVzRJ/FnmNmJuX9Q5a5ZlpldJ2mRuy+q4jVV8qCi9+SNkn4fHztb0mWStlV5jSslbZL0z5LWS+pWtKfXaB+KniLp84ren22SPifp8vjWpui9OCyus07SB6tsz0h2l7RT0ffjCUXf1/8naYWZHeruOxV9b/aU9CxJr4jP6x/hmmP5/p4l6W5J75U0I35tv4ife8jMDlL0+n+m6PdiQNHP4EETfeEjKAZft+WO3x5/PVzStfmTLOoh/LGi4Okkdx8pFVgMsj8g6WFJeytK0c5K1emW9EVJj0iarei9ut6ifdtuzV3vS4p6Lc+QdIKkf5XUKulFit7TnvjYpWZ2oLuvT507V9HWK5+VdJ+i3u6vmNlWd7+w0gsws2WKfla+Er+ObkmflvQ0M3ueuw+b2RskfUHR9+73ilKORyr6ucN05O7cmugm6SRJXubWK+llubrnxY+1lbnOovixfxzDc18Yn/OC1LGZkjZI+mbq2HcU/YPaI3f+1ZJWp+7/QtK1cblF0kZFf4AGJc2Jj18saWWV78lVueNviI+/MHXMFQUmnRWucVJ8f3dF/zgvzNU7K673itGuOUJ7r5F03wiP/1DSmgqPnRM/31MkvUVR8NIhaV9FQeLJqdfyohGeY0H+dYzwvp6UOnZd/P05KHXsFXG93+bOv1TSAyNdL/eaFqWOrcm/97lzWiXtH5/3qtzP6COjvZZxfH/vldSeOvaa+PjzcvfnVvv7VOXPylPi655T5rGPxI915I63xcf/Lfd+/lDSfEUfNv4m6clVPP82Se8ZQ3tb4+e/W9KXy7z/383VXxUfPz517Mj42Jty31eXdGbu/KsVfRAobtWU+VlS9HduWNLHcuctjuu9Mr5/vqRVtfzecWvsG2mv5vVORZ9wn6Xo0/pVki4xsxMn+Xl3eNzDI0nu3i/pHkkHpOqcomiH8M1xKqQt7qX6taJP1cW0wnJJzzWzDklHSepS1IvQr2jwrxT12uzy6bWCn+buXyKpoOgTbtpV7t43yrWOU/QJP59SuFhRkJF/n6u5piTJ3V/o7k+ppu4oLlEUfP6dokDvMUWBVTU2KOqpWmbRbKGDx/C897j7/an7d8Vff52rd5ekhWbJuKSJMLN/jlMX2xR9D4q7zT91HJcb6/f3ancfTN0v9mgUf+5XKwoKLzaz16TTQg1kP0WBzyxFQdvfqjjnRkkfiFNCTy/3vTSzF5nZtWa2QdF7NyjpEJX/vvwqd/8uSdvd/Q+5Y1IU3KYNa9eU7sWKvgfdKu9kRR+sLsr9LfqzookBJ6Re51Fm9tX49cyqcD1MEwQ/zesed78pvv1KUarnfkXBQ9FQ/LW1zPmtuTrV2lTmWL+i3oeivRSlYAZzt8/HjxfH4Vyr6J/38xQFOX9198cV/YF+QTymYi9FQVI1Hk/fcfeBuL35P4xrq7hWsbs7U9fdhxQFDvnu8GquWVPuvlXS/ypKfZ0t6SKvcuyQu7uifww3SfoPSfeY2f1W3Syh/M/AwAjH21T+529M4rEdFyhKa54u6dmKAhgp+7NXrbF+fzfm7hfTaR3xefdJeomiv6k/kPSYma2c5A8jxfd7fu54se35Nh+pKBX2k/j3rBpnKErnfVDSLZJ6zOxjZtYiSWZ2jKIPOtsk/YOi78mzFA3aL/d9Kfcz0ps+EP/eqsz5m3IBqJT8zlcKfopB6H3a9e/Rbkr+Fn1fUfr3OYqC+I1mdqkxZX7aYszPNOHubmZ3KuoFKFoXf91P0gO5U/aLv1b7R3AsNijKm3+2wuPF8TC3KhprskTS0UqCnOWSXqdojMGApBVVPu/e6TtmNkPRP4aeXD2v4lrFfxz7KBlDUZxhs4d2/cdSzTUnw/cVjd1pURQAVy3uvTk7/jT/DEnvUjRNek0cUNfazvjrjNzxPfIVyzhT0jXu/v+KB2xiyzqM9fs7qrhH9Fozm6korfJJSVea2SLPjl2plWK7j1A2iDs8/npHrv5VioKSz5rZTnf/8mhP4O7rFPUyv9PMnqpoQPknFKW1vybp1Yo+QJ2eDkzMbL5yQU0NzDez9lwAVPydz/+OF22Iv75Y5T+4bZBKHwa+IekbcdtfrCgF/xNFARGmGXp+pon4k9gRiv4oFV0Xf311mVNereif0cpJaM5Vij5l3p7qnUrf+qXSH5zrFPVAPF/Z4OdoRbOXbnD3HVU+7+ty91+r6Gf8T+N4DSsVBV5n5o6foehDw3XjuOZkuFpRuu/r7n77aJXL8chqJevjPK1Gbct7sML1K87kS5ml6NN62pvL1OtXNFh1NJP2/XX3fndfrqgXdrakyVp760+KPjy8IXf8LEXB2y4fGtz984oG5X/JzN43lidz97vd/SOKgoji93CWonRUKfg3syXKpsFrpVW7/i07U1H6s1Lwc7Wi1PcBFf4W5T8Uyt03uftPFP1eTdbvAuqMnp/mdVg89kGKZricregTX2lmjbvfY2bflvQf8RiEaxV9z0+V9DZJn3D3cp+GJupjkm5QNOPjfEWDLecr+kNykLu/JVX3WkWzbYaVzFr6i6J8/AsUfXqu1hFm9j+KxgEcomg2z3XuXu04mBJ332hmX5D0YTPbrqhr/zBFs0T+oKi3ZVzM7BpJB6bH/ZjZnkrGmRwgaZaZvSa+f4e75z/FF9s5rDH2+MTPd6SkLyv6ZHufon8s5yj6FF9tmnFM3H2tmf1O0Xu6XlHP5FmqbkbUVZI+ZGYfUfSztUTRIOO8OyTtHqfvbpK003edcVTz769Fyz2cEF/nYUUDyj+sqJfztrjOIkU9sJ9w9/NGud6Jin6v94kPHVv8fXf3n8VfB83s3xT11vUoSgkuUTQQ/t2p9FH+tf+XmQ1L+qKZtbj7Fyq0YV58zYsUjcMZlHSaot/l38TVrpJ0rqQL49+9QyT9myoHIxOxVdLnzGyBogHor1c0S+yc+IPULtz9b2b2WUnnxz1Xv1P0oW9/RR+6vu3u15rZN+Pr/0nRz+UhitLJvyl3XTQ/gp/m9ZVUeZOi2RV/7+4/ztV7u6I/XOcoWhW6oOgfxFvd/TuT0TB3f8jMjlU02+wziv6Ib1D0T+B7uerFwcw3ebwAnUdTT3+naBZRtYOdpWga8isU/UNvVTR9+T3jfBmS9FFFPWlvl/QORa/h+5I+XO3YmgqKM2LSjlA0gDmteP8Tit7LWnpM0Sfmf1G0VMJORWnIl7v7zTV+rrSzFKVLvhI/53cVBRzfGuW8TyoaEP8+RWNBfqdojM39uXrfVjTu5DNx/QcVzfgpp5bf379Keqmi8VN7Kep5+YOkN6QGwhfXPHqsiut9QtlB1++Mb1K0vIUkyd2/bmauaCr3BxR9T9/l7heMdHF3/7KZDUn6qpm1uvvnylTbqWg21lsVTXcvKPo78wZ3/0V8nV+b2XsU/Ry9WtHv+NmKpqvX2hZFPT1fVrSO1uOS3uvu+b8pGe7+kXhIQPE9dEUB6jWKgigp6iV7s6KAZ56ioPWHkj5e+5eBRmAVAmagaVi0cN21kk5299/WtzVAeWb2NkW9kQeOIZULlRawfJG7L6x3WzA9MOYHAKbGiZK+SOAD1B9pLyie7TPidOR4CjCAcXL3/MBkAHVC2gsys3Mk/c9Iddy9JgvVAQBQbwQ/kEWbf444Hdfdb5qi5gAAMKkIfgAAQFAY8AwAAIJC8AMAAIJC8AMAAIJC8AMAAIJC8AMAAIJC8AMAAIJC8AMAAIJC8AMAAIJC8AMAAIJC8AMAAIJC8AMAAIJC8AMAAIJC8AMAAIJC8AMAAIJC8AMAAIJC8AMAAIJC8AMAAIJC8AMAAIJC8AMAAIJC8AMAAIJC8AMAAILSNoa6PmmtAMbO6t2Ak1tey+8EGsrVhUvq+nvB7wQaTaXfCXp+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+RnLuudENAABMG2NZ5yc8q1fXuwUAAKDG6PkBAABBIfgBAABBIfgBAABBIfgBAABBIfgBAABBIfgBAABBIfgBAABBYZ2fOli8bLl6evvU3dWpFUuX1Ls5AAAEhZ6fOujp7dOaZaeqp7ev3k0BACA4BD8AACAoBD8AACAoBD8AACAoBD8AACAoBD8AACAoBD8AACAoBD9TbPGy5eru6qx3MwAACBaLHE6x4ho/AACgPuj5AQAAQSH4AQAAQSH4AQAAQSH4AQAAQSH4AQAAQSH4AQAAQSH4AQAAQSH4qaPurk4tXra83s0AACAoBD91tGLpEvX09tW7GQAABIXgBwAABIXgBwAABIXgBwAABIXgZwqxozsAAPXHru5TiB3dAQCoP3p+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+psjiZcvV3dVZ72YAABC8tno3IBQ9vX1as+zUejcDAIDg0fMDAACCQvADAACCQvADAACCQvADAACCQvADAACCQvADAACCQvADAACCQvADAACCQvADAACCQvBTZ91dnVq8bHm9mwEAQDAIfupsxdIl6untq3czAAAIBsHPBC1etnzUnhs2NQUAoHGwsekEFXttigHQiqVLdrnPpqYAAEwRs1GrEPzUSD51RSoLAIDGRNoLAAAEhZ4foFG0tJaKrYc9pVS+60OzM9UO+/BjpfJQz6OT3y4AqLXU37uWjplJeff5pfLg/gsyp2zfPxk7u2PPpO9m+/HbM/XeceTvRn/66lsKAADQ/Ah+AABAUEh71UBxGjuDnDEhheFScfj2u0vlwz5+YKbalTf+slS+ri/5/PK5o4/P1BvesqXWLQQQipFmTFnyd6elsyPzUMv8rlJ583MWlsq9B7dm6vUdkfy/3GP3baXyDUdfUirfPpD9n3r19sNK5QtuObFU7rgxOzTg0otfXCr/y2XlXwLBTw0Up7cvWnplnVsCAABGQ9oLAAAEheBnEnV3dWrR0itZ3RkAgAZC2msSFdNhwEQMPfBg5v5L9juqVB5+wTGl8skrfp+p953bnlcqP/kf7yuVC9uz00IBTBPpcTru1Z0yMzXNfE4ydmbokP0z9TYcOatU3rooOT6492Cm3ouOuLNU7hzoL5UPaBnO1Lvp+kOTpi5Pxg297P2vTdrWn7124fEnSuWD+m9Pzh8a0ljR8wMAAIJC8DMBbFgKAEDzIe01AZU2LCUowlRpvXZVqbz86dnpngfvv6FUvuMLRyQPpHrDD/9ENqU29NjjtW0ggEmTTllJ2WnnNjv5e7D9qO5MvY2HtpfKfXsnfxCG5qXSR7mZ7u3rk/Lsh5MH5/4x24fS84VkervWbSwVd/Rn014H7bixVE6nrbK1Jg89PzVUHOAsMd4HAIBGRc9PDRHwAADQ+Ah+xonUFhrd0MOPlMqHvD0pt3bNK5XvOu+wzDmFjmQ16cM/m6TA8jPOANRAanaWzZhRKrfMmpWpVnhykrbasTB57PFjs6smtwwm19u5MJkp1dG1M1PP70jSXt3XJSmnzp5kpWXrG8ie81BPUk6lqXw4m6garnKWWb0R/IxTpfE+AACgsTHmBwAABIWeHyAww72bS+WnnLsy81jLbruVyvf/y9NL5Z377J2p96TLCqVy+29uqnUTgebWkqSjWma0Zx9KpZ13HpHMjFp3dDJza9tBuUX7Zia/b22dSTpqeH12ttecB5Pn3f+7Sb32R3sz9XxTksb2gaReYWcqPTbSxqZNktoaCT0/AAAgKAQ/AAAgKAQ/AAAgKIz5AVBS2Lq1VD7gE3+sWG/zG44rlQfe/txSed4D2Y0IZ/zm5uTONBgngMClx8FYtu8gPbanZd9kjNy2p2XHy208NPm3274t9TuRulzno9l/zd3X7SiV+/dIpsTPuXdTpp4/mJqOPpD8Lg4NZX8vM7+Llcb2TPPfV3p+AABAUAh+AABAUEh7ARizeRetLHu87UkHZu6vOydJj1mqF/3pv7hRp++4RUv671WnD6pPbbpGB+hnOkRrbc6ktBnIqJTC8kKmWnf7gE4fukMvHF4T/axau5bPPFg/O+Rkre3YvVRv+37Jiv+9T07+tQ5nZ6Orc33yi5BOe817ILXS8n3rM+dofbJBaNtgahPQHTuy9caTqprm6a1K6PkBMKWe98id+trGn+mUnXdqtg+qRdJsDemlekDf1NV6lq+tdxMBSdKzfK2+3n+FXjp8n2Yr/ln1QZ2y805947YL9Kzee+vdRIwTwQ+AKdO9Zb0+e+331KEhtSv7ibNdrg4N62NaqX19W4UrAFNjX9+mj/mf1KHh8j+rhUH9230/0b47N1a4AhoZaS8ANZPfAHX31H2bOVP/MHiD2oaH8qdltJrrVd29+vJhJ2rmH+8sHS/ku/iBMWrp6CiVC6nZUK3z5iaV9ttLkvTaR3+l1o0jp4RaNayXtNykz5zwanlqj9H2Lcl5uz2SvcaMzcnP/6x7niiVfVNv0rbtfZlzMpuHptNygaasaoGeHwBT5oXDD+zyKTqv3Qs6ee1fpqhFQHlLem9Tuwoj1mkvFPR3d9w8Yh00JoKfUax6qFeLly2vdzOAaaFTI/f6FM0aHhi9EjCJOgvV/QzOHuif5JZgMhD8jGDVQ72SpJ7evpErAqhKX5WZ9h2tM0avBEyivpbqfga3z5g5eiU0HMb8jGBgaFjHHbSHurs6tXjZcq1YuqTeTQKalvf36xodoJdq5NTXoEzLW5+kjhvu1cBzDk3Oz61E27YjGbNhN92R1BuqrncJ00j6ZyM1DsbaswGMpx5rnTO7VB4+ONldffsBsyRJV/U9U6+4/89q98qpr0Fr0ZUHP1PtW117rkpWR0+3oWVLbvzOw48mzzucXNsHUz1N+VWXGdtTc/T8VGHF0iX0/gA18DMdouFR/uwMq1WXzXjaFLUIKO/iQ07QUEvriHWGWlt10dNPmKIWoZYIfgBMmbU2R5/UcdqpVg0q++l2UKadatWnZ71Aa1vnVrgCMDV65izQR597tvpa2jWY28dr0FrU19auD5z8Jj0yb0GdWoiJIO01DouXLVd3V+foFQHs4kbbV2/zk/Vq3asX6UF1akh9atNvdaB+roO1tm+O1LdFktR67aqK12mdP79U7nvhUaWyDWdTBG19yTThtlvvL5WHt2yZ6EtBo6iQFspMEZfUOndecmfvJGgpzEh6eAbmJIHO9Qcfrtce8H6d9dfrdeo9N2v2YL92tM/UFYc8U5cuPF49sxZo1voodeUtSTDfum2EwdKtyXP5zp2lcjpF5/mNSFFzBD/j0NPbpzXLTq13M4Cmtdbm6HwdrfN1dL2bAozokXkLtOyE07XshNMzW7TMemLkafBobKS9Kli8bLlmtI2c7wUAAM2Hnp8Kenr7dMwBXZlji5ctV09vHykvoAEMb9pUKs/49U2VK6YGrRaOTmaPWUvy2a/Qnv0c2JpKjxW2bhWaSGqmVEtnR/axttS/vFSqbGh2crxjUzZVNtSR/GzM2JY8Nuve3LYW7ck1bMv25GlSqS1JKvSl7qfampnthUlH8FNGuTE9xfukuwAAaG4EP1JpBefiOj6lMT1XfbJUhzV+AACYHhjzo2QFZ7axAABg+qPnJ7Zi6RItWnplvZsBoNYKyTgNv/n2slUst5hdywHdyelHPCk5npvC3LI1GdsxlFq5N/2cGEWF1ZlHrJdbd8fay/8rs1mzsge6kvWjCh3J1PKZ65NFbPv2nZ05ZcbW5HvZMpCa4ZWfEPNYskN7esd45abbW2pKvA9VuXJzte8RqkbPDwAACErwwQ8LFgIAEJbg014sWAggn6YaWvNQcidVzi9rZ6lVptsW7Z/Um5tNt9iO/uTO5mTq/PATGyq2IRjjSOOkU0eSMqklm5nssm4z2rP1+pO0ZWH3JL3VMpR8Zzsf25G9dEfybzK9EnShM3vt9Eapvj43DT5TMZUuS6fERnofSHXVXNA9P/T6AAAQnmCDn/z09vRxAiIAAKavYNNeldJdpMEAVCu9yrTS5RGk0zJtqVll6ZSMJCmdskk9NrRufbZeqOmylPQGpi2pmV8+kN0g1FKzptrXJt+vwrwkZWUDQ7mLJymntq1J+rKld1u2Xv77V2pD9rgXUiks0ll1E2zwUw49PgAATH/Bpr3KWbF0SV1Wcu7u6mSBRQAApgg9Pw2ABRaBcHh/kjrJzCobgaU25Gzbe8/sg/nNOwPg+YUD25IUYTrVlX7fJMm3Jakqm5cseNiyMbV5bWtuk9vUJqWFLUm9wmA2pVbIbWCKxhZc8EMPCwAAYQsu+Cnu4wUAAMI0rcf8LF62nJ4eAACQMW17fopBDzu2A2h2PpRMvx5a+1gdW1JHI0wL98GBsmXt2FGmdqx3cy1ahSY1bXt+enr7MjO3ikEQ09kBAAjbtO35KacYDDGzCgCAcE3bnh8AAIBypmXwk96fq7urU4uWXkm6CwAASJqmaa/0/lzlVmwmEAIAIFzTMvgZTT22sAAAAI1hWqa9AAAAKpk2wQ8LGgIAgGpMm7QX21YAAIBqTJueHwAAgGo0fc8PqS4AADAWTR/85NNd6TV+AAAA8po++MlLr/EDAACQNy2CH3p6AABAtaZF8MOihQAAoFrM9gIAAEEh+AEAAEEh+AEAAEEh+AEAAEEh+AEAAEEh+AEAAEEh+AEAAEEh+AEAAEEh+AEAAEEh+AEAAEEh+AEAAEEh+GkQ3V2dWrxseb2bAQDAtEfw0yBWLF2int6+ejcDAIBpj+AHAAAEpamDn8XLlqu7q7PezQAAAE2krd4NmIie3j6tWXZqvZsBAACaSFP3/AAAAIwVwQ8AAAgKwQ8AAAgKwQ8AAAiKuXu92wAAADBl6PkBAABBIfgBAABBIfgBAABBIfgBAABBqXqFZzO7TdLOSWxLLSyQtL7ejRgFbayNDnd/Wr0bAQBoPmPZ3mKnux87aS2pATO7iTZOXLO0sd5tAAA0J9JeAAAgKAQ/AAAgKGMJfr45aa2oHdpYG7QRADBtscIzAAAICmkvAAAQFIIfAAAQlKqCHzM7xczuNrP7zGzpZDdqlLasMbNbzWx1cbqzme1uZleb2b3x1/nxcTOzr8TtvsXMjpmkNn3XzNbFayEVj425TWb2prj+vWb2pilo43lm1hO/l6vN7GWpxz4ct/FuM3tJ6vik/SyY2f5mdq2Z3WFmt5vZe+PjDfVeAgCanLuPeJPUKulvkg6SNEPSXyUdPtp5k3WTtEbSgtyxz0laGpeXSvpsXH6ZpF9JMknHSfrzJLXpBEnHSLptvG2StLuk++Ov8+Py/Elu43mS3l+m7uHx93mmpCfF3//Wyf5ZkLSvpGPi8m6S7onb0lDvJTdu3Lhxa+5bNT0/z5Z0n7vf7+4Dki6WdFoV502l0yR9Ly5/T9IrU8e/75GVkrrMbN9aP7m7Xy9p4wTb9BJJV7v7RnffJOlqSadMchsrOU3Sxe7e7+4PSLpP0c/BpP4suPtad18Vl7dKulNStxrsvQQANLdqgp9uSQ+n7j8SH6sXl/QbM7vZzN4WH9vb3dfG5cck7R2X69n2sbapXm19V5wy+m4xndQIbTSzRZKOlvRnNc97CQBoAs044Pl4dz9G0kslvdPMTkg/6O6uKEBqGI3YptjXJD1Z0lGS1kr6Ql1bEzOzOZJ+Lulcd9+SfqyB30sAQJOoJvjpkbR/6v7C+FhduHtP/HWdpMsUpWIeL6az4q/r4ur1bPtY2zTlbXX3x9192N0Lkr6l6L2saxvNrF1R4HORu18aH2749xIA0DyqCX5ulHSwmT3JzGZIOlPS5ZPbrPLMbLaZ7VYsS3qxpNvi9hRn9LxJ0i/i8uWSzo5nBR0naXMqfTLZxtqmX0t6sZnNj9NPL46PTZrc+KdXKXovi20808xmmtmTJB0s6QZN8s+CmZmk70i6093/K/VQw7+XAIDmMequ7u4+ZGbvUvTPo1XSd9399klvWXl7S7os+h+pNkk/cverzOxGST81s3+Q9KCk18X1f6loRtB9knZIevNkNMrMfizpJEkLzOwRSR+XtGwsbXL3jWb2KUUBhiR90t2rHaA83jaeZGZHKUojrZH0T3Fbbjezn0q6Q9KQpHe6+3B8ncn8WVgs6Y2SbjWz1fGxj6jB3ksAQHNjewsAABCUZhzwDAAAMG4EPwAAICgEPwAAICgEPwAAICgEPwAAICgEPzVgZnukdkZ/LLVT+jYzu2CUc7vM7B1T1VYAAELHVPcaM7PzJG1z9/+ssv4iSVe4+9Mms10AACBCz88kMrOTzOyKuHxevHnodWZ2v5m9J662TNKT456iz9evtQAAhGHUFZ5RU4dKeoGk3STdbWZfk7RU0tPc/ah6NgwAgFAQ/EytK929X1K/ma1TtF0HAACYQqS9plZ/qjwsgk8AAKYc/3zrb6uiNBiaD7MF0Eis3g04ueW1k/c7YZVfXsvMmaVyoT/5jNkya1b2EjNmJHf26ErKm7dVftoZ7aWyb03q+cBApl5h586K10D9XF24pOwPDj0/debuGyStMLPbGPAMAMDkY6o7MH788qCRTO+enyq1dHSUytbZmXnM+/qSx9K9QO1JEmSXc7ZsLfs8u/T8pHqcrLU1qTc0VEWrMVno+QEAABDBDwAACAzBDwAACArBT46ZbcvdP8fMzk/dP9fMzo7Lu5vZ1WZ2b/x1fpnrnWxmN5vZrfHXJanHnhkfv8/MvmIWTWeodF2LfCWuf4uZHVPhNfy2XFsAYLrzoaHk1teXuam1tXSz3eaUbmppLd18y9bMLX1Oob+/dEsfV2urWmbOLN0yzCrfUDcEP2NgZm2S3iLpR/GhpZKucfeDJV0T389bL+nv3P3pkt4k6Qepx74m6a2SDo5vp4xy3Zem6r4tPr+cH0his1Q0vnPPjW4AMIUIfsZmiaRV7l4cvn+apO/F5e9JemX+BHf/i7s/Gt+9XVKnmc00s30lzXX3lR5Nuft+6vxK1z1N0vc9slJSV3ydvMslvX58LxGYQqtXRzcAmEIscjg2iyXdnLq/t7uvjcuPafTtKl6tKHjqN7NuSY+kHntEUvco1+2W9HCZc9amjsndN8UB1h7xOkIA0NzSaaIRlmixdNppeDjzWHra+fD65E9jy9y5SaX27L9Fa08WOWxpSx7zreWnwO+C5WQaEj0/1Sn+9O4r6YmyFaLem4o/5WZ2hKTPSvqnMT3xKNcdwTpJ+43jPAAApjWCn131mVlq9SvtrmjcjiT1SepIPfZ4Me0Uf11X7oJmtlDSZZLOdve/xYd7JC1MVVsYHxvpuj2S9q9wTl5H3F4AAJBC8LOr30k6S5LMrFPS6yRdGz92p6SnpOpermgQs+Kvv4jPe7aZfT8ud0m6UtJSd19RPDFOa20xs+PiWV5nF8+vdN34+NnxrK/jJG1OpcdK4uvtI2nNOF4/ADQe99LN2toyt4zh4dLN3TO3ioaGSjff2Z+99Q+Ubjaro3RLz/wq9Pdnn2t4uHRjhldjIvjZ1XslnW5mqyWtlHSJu18fP/YrSSek6i6TdLKZ3SvpRfF9STpASa/LuxQFTB8zs9Xxba/4sXdI+rak+yT9Lb7+SNf9paT74/rfUmpGV9zeomdKWpkamA0AAGIMeM5x9x5JL6/w2INmtsHMDnb3e+PBxC8sU/U5kv47PufTkj5d4Xo3SXpameNlrxuP/3lnhWsdlbr7RkkXlKsHAEDoCH7Gbqmigc/3Vqrg7h+YuuaUdZu7X1PnNmAKLF62XD29USdjd1enVixdMsoZAACCnzFy97sl3V3vdozE3b9V7zZgavT09mnNslMlSYuWXlnn1gBTw3NT2NNjaQo7d5bKLbNnZ6t1JNPgffuO5PjM1ByX3Bii9GO+bXup3NrVlW1Tasd4pXd1H6GtTIOvH8b8AACAoNDzAzSZfKoLADA2BD9Ak0mnuoAgjbTCc3oV5qHshFffNlS+3nAhKe/YoUp8Z3+p3LLbnOyDqfRWob9faGykvQAAQFAIfgAAQFBIewEAGt84ZkmVXf25+NisWUl5Rnvq+O6ZU9IzvNLn+OBgtt7AQFIvNdsrPfMrX4+ZX/VDzw8AAAgKwQ8AAAgKaS8AQONLpYVs5szsQ6nZVZnH8gsM5lJQZfUPZO/Pn5c8z6zk2rZ2Q6aazUgthjg0wraKpLcaAj0/AAAgKAQ/AAAgKAQ/AAAgKIz5AQA0Fc+voJze2HR7MjU9v7Fpxett3ZZcalZuy5jNyWMtg6mxPK3ZvoNCenxRemPTXFvTY5J2eR2YMvT8AACAoBD8AACAoJD2AgA0tUorKuennLfMnZucMzOZmq7UCs/56fGFRfsk59/3SPnz88+Vn2KfNtJjmDL0/AAAgKAQ/AAAgKCQ9gIANL4qNwFNbxzaktqIVNo1VVU6Z+fO5E5qRWdJalnzWHJn3z1LxcL9D43U2tRzVl6NGvVDzw8AAAgKwQ8AAAgKwQ8AAAgKY34AAI2vBruhFzb1lsotXcnYHuvoSOp0ZMfopKfRFzqTKfGtC/aoeO10Sy09VkmSVzl2KduIcZyDEdHzAwAAgkLPDzBNdHd1atHSK0vlFUuX1LlFANCYCH6AaSId7BSDIKDhjSOlY23Zf135lZwrHW9JTzuvsNKy7ay8aWrL/Y8m186db6lp9Zaazl7YsSPXqNRrrPa1k+qqOdJeAAAgKAQ/AAAgKKS9AAD1M46Uzi5prlT6KD07qyW3urLak395PlxIyvumZm61ZPsEvDWV9ko9Ztuy6axC7+bksRnlV5LeBemsuqHnBwAABIXgBwAABIXgBwAABIUxPwCA5pYaO5MeD1TI76Ceut8yZ3apbEPJ+B/r68ucMrxgbupOMr09sxO8JEutDO3bk/FALfPnZ+oVNm0q21ZMLXp+AABAUAh+AABAUEh7AQCmJc+lvVrnJims9GamPjBYKg/u05U5p2391lK5MH9Ocq3cVHefk1rhuSWZbu9bt2brVVpZOrVqNemwyUfPDwAACArBDwAACAppLwDA9JTeOFTZ2V/pT/7WlqSp2jZnZ3sN7jcvOacvSUel01ySZH3JtTPrNqdWnI4eLL+qcybVlWs3K0HXHj0/AAAgKAQ/AAAgKKS9AADTUy5d1DIrlapKz65KbXhamJ3dDLV1R5KOsv5kVphmtGefKzf7q/ScXfMy99PprfRsNGZ7TS16fgAAQFAIfgAAQFAIfgAAQFAY8wMACEJhRzIup3XPBckDqfE6hb2zY3TSBndPxgN1PJ4b45Oe+r6ht1TMT1K31DT29GOM85la9PwAAICgEPwAAICgkPYCAExL6enjkuQDA0m5PynbrGST08E52XM61m5L7hSStNdwZ3aqe/sTvcm108+bS2cVcputoj7o+QEAAEEh+AGmoe6uTi1aeqUWL1te76YAQMMh7QU0gcXLlqunN9pwsburc9T6K5YukSQtWnrlpLYLaGQ+PJy5b6lNRn3r1uSBeXNKxc6/bcheJLWSs7el+gvyXQfpzUgLqefNpd7SbUhveuqkw6YUwQ/QBHp6+7Rm2an1bgYATAukvQAAQFAIfgAAQFBIewEAgpOect66OZnO7t17ZupZ34DKGW7P9R2kdob33fcqlVvWbcper9Lu7ekxQ3meXycaE0XPDwAACArBDwAACAppLwBAc0unjEZIEWWmvqfqWWvSD2DrN2fP6UxWdW7dkVohum1mpt7gfvNL5ZaB1PN0dmTqtaTuF9atT66X2nR1l9dQKSVGOmzc6PkBAABBoecHAIAK9h3YpNM33aglW+9QZ2FAfS0zdM0eR+qnTz5Rj87ao97NwzgR/AAAmlul9M8I6aPMrKsdfcnxubuVysf2rdFHN/xSrT6sdhUkSbMLA3rp+lU6eeNf9bGjz9Kf9zo0usbs5HozNg8m15s3O9uE1Oyxlq55yQOplJzn2p3ekJVUV22Q9gIAIGffoV59dMOV6vDBUuBT1O4FdQ4P6pN/+aH2276hwhXQyAh+AADIedWWVWr1woh12grDet0Dv5+iFqGWCH4AAMhZsuOuXXp88tq9oBc/umqKWoRaYswPACAINmNGqZzeRb2QmmbeEu+03unJuJ2RzBoakBVcg7slO7S3DCblwbnZMT8t/Z2lcuedqent6d3eczu8Z3ajT68KjXEj+AGmse6uTi1aemWpvGLpEknS4mXL1dPbt8txAJE+tWu2Rg+AdrTNGLUOGg/BDzCNpYOaYhAkST29fVqz7NRdjgOILJ/5FJ3Sf5faVXl21aC16Df7HT2FrUKtEPwAAIKQmTKemvaeXvm5mAL7mR+sk3WP2pVarTlnyFp12Z7PVfu2IbVvS9JR2/ZPVnFuGcoGT7s9uLVUHjho71K5dUdqheg1a7NPNJhcO70hq+fSY6geA54BAMhZ27KbPtX+fO1UmwaV3V5iUC3qa2nXJw99vdZ2stBhM6LnB0AJY4GAxI2t3Xp7+9/p9IE79MKh+9WpQfXZDC2f9VT99NAXEfg0MYIfACWMBcK0Vml15HQKLJdK6lGrvqqn66v2dLV0xOmsgtTy+LData5Ub3ifJG3V1pcMgh6cnU2w7Fg0t1Ru3Zmk1Fr6U/Xmz0ufosLDjyZNTa9MnU7jSaz+PAakvQAAQFDo+QECkZ/2DgChIvgBAlFp/E6loKjSGkHAtDNSuiiVEivs3Jmcsj67p1fLQLImUGfrPsk5B8zK1BuckyRc2rcns7i8PVnIsDAve07rtq7kzswkpWZbsmsMDff2JndIgY2I4AcIXKWgptIaQQDQ7BjzAwAAgkLwA2BUxRTY4mXL690UAJgw0l4ARlVMgZH+QnAqjJ3JT4kvbNlSKrf+LZnCvtuOvTL1th6aTGPfuUd7cs7O5HmGZnVkzpk7mFxvePbMUrm9Jdt/YakNWpVatTq9gnV0gPFA9PwAAICgEPwAAICgkPYCUDWmvyNoZhUfSqfBPLUKc8v2vky99q1zSuWhWcn09g1HJOcUcv+ZZ2xOztl8UJIq2/2u1my9ziQlpic2JeWhoUy9wtZkc9VMSiydDsu/1mmWKiP4AVA1pr8DmA5IewEAgKDQ8wMAQDWqTP0U0rOu1j6eeaxzW/KYd+9ZKg/OSjY83bIo2y/x6POTVJcNp9JRPjNTb25HkgabvXlbUi03My29OarNTK5R2L49de3ca02nwaZBCoyeHwAAEBSCHwAAEBTSXgDGhZlfAJoVwQ+AcWHmFzC6Qm68TUtrMi6nZXMyxmb2Y53JOW3Z3dqHkodUSIb/qH9+djr6xlnJg4UZC0vlOff2ZtuwcXNy7cfXlW/4NJ/qTtoLAAAEheAHAAAEhbQXAAC1NEKKKD0NvrDmoVJ5ZmoKfMdd7Zlzth/ZXSr3z0/SZtv3z6am2rcm99cdk/x733jogky9fW5IVoyesXuy0ao/2JOUc+m6XTZHLT3QnOkwen4AAEBQCH4AAEBQSHsBADBVKqSJhp94olROr8AsSZ3XJ7OzZu27V6k8c1M2nWWF5NqbDklmjG05JJuyemh+sqpz151JufWo3Uvl3R7cmTmn7a4kRVfYvKXsa5Akz22imjSusWaP0fMDAACCQvADAACCQvADAACCwpgfAAAaSH7cTOb+fQ+Uip2PrM3Ua+lKpq139Mwvldv65mfqbTkoKW86LCm3b036Q7btNytzzoI5yUkzNybT4FvXbszUK2zqTdqdmi5fcap8ndDzAwAAgkLwAwAAgkLaC8CUW7xsuXp6+3TpQ7065oCuejcHaEr5TVMLqU1Kbf2GUnnvrftm6i1YMLdUHpjfUSpv3T+ZHt/71Oxz9ZyYrDrdsT6pN3PTnEy9rvuTqfgz1qxPHhguZNuaal9hZ3Za/VQg+AEwJYoBjyR1d3VqzbJTtfJHS+vcKgAhIvgBMCV6evu0Ztmp9W4GABD8AADQlEZYJTk9Q2zowYczj7VuStJenXN3K5VnbE5WeJ71RGfmnHXPTMKFnXslz7vtwGwbNh6ZrBg9q2f/UnnuQ9m019w7kudtW78paXf/QKk83Nuriia4QjTBD4CaSqe30rq7Onc5NqOtVSvv36Azl16p7q5OrVi6ZCqaCCBwBD8AJqy7q1OLll5ZKleb3ioOdl6z7NTS+QAw2Qh+AEwYPTZAA0hvHjpCWmh4S2pj0q1bk3LPo6Xi7K6uzDkHPr6wVN50aDLDa+ce2RVzBhYn1/MDk1TXS8++OVPvO388oVTuujVJt+29MtnEtaV/78w59ngyQ2x4fWom2ThSYKzzAwAAgkLwAwAAgkLwAwAAgsKYHwAApoPxTP+ucI73ZWds2h33l8p73J9MZ7fUVHlJ2n5vssLzg69Ijl/Y99xMvec/465SufOYwVL5wH9KNkq9YdOizDkPXHZwqbz3n/crldue2JKpp83bNBp6fgAAQFAIfgAAQFBIewEAgIwRNxvdvj0pb9iYeajj4Z5S+bAb5pfKg4cuzNS7e9ERpXLfnskU/b7X3lYqP6Prkcw5n3vfpaXy73c8uVT+9IqXZ+rN+lt2I9dyCH6ABpXfCHS6Sy+UmD/OOkIAaongB2hQoW0EWinAqWbl52KgSKAEoBoEPwCaXjFQZIsMoL7SG6oOP/FEqdySKkvS/BuSGWML9lxQKj/6l6eUyvftfXjmnB8fmawKPdyRzFJ7/nPvyNT7/suuT917X9l2EvwAaGj5fcPo2QEwUQQ/ABpaOthZvGx51YFQfswUQROAIoIfAE1jpEAoLz1minQYgDSCHwBNqVxPTj5FBqAxeX9/qTyU2k2+7bHHS+V5M2dmzpn3v8OlcsseyU7waw96cqbec/dPxgr9+Yfln5/gB8C0USm1xTR6AGkEPwCmvYlMowcw/RD8AIAYIA3UTWpz1fRU+XQ5r5BOleVWmZ53Q2qzVtJeAFAZA6SBcLCxKQAACAo9PwBQAakwoPGNuAlrBQQ/AII12tT4dCpsLAssAmhsBD8AgjXa1Ph0QJSuO5bNVovXI1gCGgfBDwDkjBaoVLPfGAOogcZF8AMAYzTWXiAAjYXgBwAmYKTVo/N1SH8BjYHgBwAmoJpgplinVr1ExfFEBFPA+BD8AEAdTGRAdHE8ESk3YHwIfgCgDmoxIDqdciteC8DoCH4AoM4qjRvK18n3DpHyAsaH4AcApkilRRWrCWJIcQG1Q/ADAFNkIj01o61GDaB6BD8A0ARIcQG1w67uAAAgKAQ/AAAgKAQ/AAAgKAQ/AAAgKAx4BhpIftVfAEDtEfwADSS96i8AYHKQ9gIAAEExd693GwAAAKYMPT8AACAoBD8AACAoBD8AACAoBD8AACAoTHUHxsnMbpO0s97tGMUCSevr3YhR0Mba6HD3p9W7EUAzIPgBxm+nux9b70aMxMxuoo0T1yxtrHcbgGZB2gsAAASF4AcAAASF4AcYv2/WuwFVoI21QRuBaYQVngEAQFDo+QEAAEEh+AEAAEEh+AEmwMyeZWZDZvaaerclz8xOM7NbzGy1md1kZsfXu015ZvaGuI23mtkfzewZ9W5TnpkdamZ/MrN+M3t/vdtTjpmdYmZ3m9l9Zra03u0BGh1jfoBxMrNWSVcrWujwu+7+szo3KcPM5kja7u5uZkdK+qm7H1rvdqWZ2fMk3enum8zspZLOc/fn1LtdaWa2l6QDJb1S0iZ3/8/6tigr/jm8R9LJkh6RdKOk17v7HXVtGNDA6PkBxu/dkn4uaV29G1KOu2/z5NPNbEkN90nH3f/o7pviuyslLaxne8px93XufqOkwXq3pYJnS7rP3e939wFJF0s6rc5tAhoawQ8wDmbWLelVkr5W77aMxMxeZWZ3SbpS0lvq3Z5R/IOkX9W7EU2oW9LDqfuPxMcAVEDwA4zPlyR9yN0L9W7ISNz9sjjV9UpJn6pzcyoysxcoCn4+VO+2AJj+CH6AKpnZO+PBw6slHSvpYjNbI+k1ki4ws1fWsXmSsm00s/2Kx939ekkHmdmCOjZP0q5tjMcjfVvSae6+od7tkyq/jw2qR9L+qfsL42MAKmDAMzBBZnahpCsacMDzUyT9LR7wfIyk/5O00Bvol97MDpC0XNLZ7v7HerdnJGZ2nqRtDTjguU3RgOcXKgp6bpT09+5+e10bBjQwdnUHpq9XSzrbzAYl9Uk6o5ECn9jHJO2hqOdMkoYabfd0M9tH0k2S5koqmNm5kg539y11bVjM3YfM7F2Sfi2pVdHMQwIfYAT0/AAAgKAw5gcAAASF4AcAAASF4AcAAASF4AcAAASF4AcAAASF4AcAAASF4AcImJltq/H11jTCKtKSZGbnmNn5cfk8M+uJV2y+w8xeX+/2Aagfgh8Aofiiux+laMfzb5hZe53bA6BOCH4AyCKfN7PbzOxWMzsjPt5iZheY2V1mdrWZ/dLMXlPF9TrN7Fdm9lYzm2Nm15jZqvjap8V1FsXXvcjM7jSzn5nZrPixNWb2ubj+DfFWHTKzvzOzP5vZX8zst2a291hfq7vfK2mHpPljPRfA9EDwA0CSTpd0lKRnSHqRpM+b2b7x8UWSDpf0RknPreJacxTtI/Zjd/+WpJ2SXuXux0h6gaQvWLyXhaSnSrrA3Q+TtEXSO1LX2ezuT5d0vqQvxcf+IOk4dz9a0sWSPjjWFxrvc3avu68b67kApgeCHwCSdLyiYGXY3R+X9DtJz4qPX+LuBXd/TNK1VVzrF5L+x92/H983SZ8xs1sk/VZSt6Rij83D7r4iLv8wfr6iH6e+FoOuhZJ+bWa3SvqApCPG8BrfZ2a3S/qzpH8fw3kAphmCHwC1tkLSKanenTdI2lPSM+MxN49L6ogfy28u6KOUvyrp/LhH6J9S16nGF939CEUbvn7HzMZyLqrj3Lg10K0igh8AkvR7SWeYWauZ7SnpBEk3KApkXh2P/dlb0klVXOtjkjZJ+u/4/jxJ69x90MxeIOnAVN0DzKzYq/P3itJaRWekvv4pda2euPymal9cmrtfrmiX9nGdD6D5EfwAkKTLJN0i6a+Slkv6YJzm+rmkRyTdoSgttUrS5iqu915JnWb2OUkXSTo2TlWdLemuVL27Jb3TzO5UNAD5a6nH5sepsvdKel987DxJl5jZzZLWj+N1Fn1S0r+YGX8DgQCZ+4g9QwACZ2Zz3H2bme2hqDdocRwYTfS6iyRd4e5PK/PYGknHuvtEAhxMPf6hoJFYpQfaprIVAJrSFWbWJWmGpE/VIvABgHqi5wfAmJnZZZKelDv8IXf/dR3a8mZFqbG0Fe7+zqluC+j5QUOp2PND8AMAqBX+oaCRVAx+GOwHYFzy+4Kl99KK759rZmePcH5rvFLzFVU8l5nZV8zsPjO7JV6osFy935oZKzej9s49N7phWiD4AVBzZtYm6S2SfjRCtfdKurPKS75U0sHx7W3KzgpL+4Gyq0QDtbF6dXTDtEDwA2AyLJG0yt2Hyj1oZgslnSrp21Ve7zRJ3/fISkld8fYbeZdLYsd2ACMi+AEwGRZLunmEx7+kaF+uQpXX65b0cOr+I/GxDHffJGlmPC0fAMoi+AFQS8UBr/tKeqJcBTN7uaIVn0cKjiZinaT9JunaAKYB1vkBMF59ZjbD3Qfi+7srWXW5T/G+W2a2v6Jd3iXp64q2t3iFmb0srjPXzH7o7meN8Fw9kvZP3V+oZJuLvI74+QGgLHp+AIzX7ySdJUlm1inpdUp2fb9T0lMkyd0fdvej4tvX3f3D7r7Q3RdJOlPS8nKBj5m9ysz+I757uaSz41lfx0na7O5ry5xjkvaRtKaWLxTA9ELwA2C83ivpdDNbLWmlpEvc/fr4sV8p2hx1TMzsk2b2ivjukyVticu/lHS/pPskfUupGV3x8xc9U9LKSgOtAUBikUMAkyReBfqD7n7vOM//oaT3uXvZsUMVzvmypMvd/ZrxPCcmbPr+QznppOjrdddp8bLl6untU3dXp1YsXVLXZmFELHIIYMotVTTweVzc/ayxBD6x2wh8MNl6evu0Ztmp6ullaFmzYsAzgEnh7ndLunuKn/NbU/l8AJoTPT8AACAoBD8AACAopL0AABjFqod6NTA0rDOXXqnurs56NwcTRPADAMAoBoaGddxBe2jNslPr3RTUAGkvAAAQFIIfAAAQFIIfAAAQFIIfAAAQFIIfAAAQFIIfAAAQFIIfAAAQFNb5AQCgjOLu7ZJ0aVtrnVuDWiL4AQCgjOLu7ZKklZ+vb2NQU6S9AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUAh+AABAUNjVHQCAceju6tSipVeWyiuWLqlzi1Atgh8AAMYhHewUgyA0B9JeAAAgKAQ/AAAgKAQ/AAAgKAQ/AAAgKAQ/AAAgKAQ/AAAgKAQ/AAAgKKzzAwBAbPGy5erp7ZMULVyI6YngBwCAWE9vn9YsO7XezcAkI+0FAACCQvADAACCQvADAACCwpgfAAAmiB3emwvBDwAAE8QO782F4AcAADSU4pIDk9WLRvADAMAkya8bRDqsOsUlByarF43gBwCASZJeN4h0WONgthcAAAgKwQ8AAAgKaS8AADCpGm3sE8EPAAA1lF/zB4039ongBwCAGqp3rwZGx5gfAAAQFIIfAAAQFIIfAAAQFIIfAAAQFIIfAAAQFGZ7AQAwBfJT4JkVVj8EPwAATIF0sFOrtW4abfHAZkHwAwAIWj6AaKS2jBbMNNrigc2C4AcAELR0AFEP+YCnXDATQg/PaEFoLd8Dgh8AAOqomuArXWfxsuXTcvuM0d6HWvZyEfwAANBEpmOvT62ke4dGCqQIfgAAQF2Mluoa6yax1aYwCX4AAEDNVTOQfLRgZbJ6uQh+AABAzVUKbMbamzMZCH4AAMCUqXVvzniWKiD4AQAATWs8SxUQ/AAAMMUaIfUTMoIfAACmGNPV64td3QEAQFDo+QEAoAE1Y2psqvZJm+h7Q/ADAEADqlVqbCr3BZuqfdIm+hoIfgAAmMbY+X1XjPkBAABBIfgBAABBIe0FAMA0kB8EPFXT6adqkHMtEfwAADANpIOdqRzbM1WDnGuJtBcAAAgKPT8AAARoIlPgmzHVlUbwAwBAIPLjgsY7Bb4ZU11pBD8AAASCPcUiBD8AAKAkndKqpBlTXWkEPwCA4DT7mJXJ1OwprWoQ/AAAghPCP3hUxlR3AAAQFHp+AACYZvKzupBF8AMAwDTDrK6RkfYCAABBIfgBAABBIfgBAABBIfgBAABBIfgBAABBIfgBAABBIfgBAABBYZ0fAAACF9qiiAQ/AAAELrRFEUl7AQCAoBD8AACAoBD8AACAoBD8AACAoBD8AACAoBD8AACAoBD8AACAoBD8AACAoBD8AACAoBD8AACAoBD8AACAoBD8AACAoBD8AACAoBD8AACAoLTVuwEAAEyFxcuWq6e3T5LU3dVZ59agngh+AABB6Ont05plp9a7GWgApL0AAEBQCH4AAEBQCH4AAEBQCH4AAEBQCH4AAEBQCH4AAEBQCH4AAEBQCH4AAEBQCH4AAEBQCH4AAEBQ2N4CADBtsZ8XyiH4AQBMW+znhXJIewEAgKAQ/AAAgKAQ/AAAgKAQ/AAAgKAQ/AAAgKAQ/AAAgKAQ/AAAgKAQ/AAAgKAQ/AAAgKAQ/AAAgKCwvQUAYFphPy+MhuAHADCtsJ8XRkPaCwAABIXgBwAABIXgBwAABIXgBwAABIXgBwAABIXgBwAABIWp7gCAppRezyeNtX0wGoIfAEBTYj0fjBdpLwAAEBSCHwAAEBTSXgCApsG+XagFc/d6twEAAGDKkPYCAABBIfgBAABBIfgBAABBIfgBAABBYbYXAKAmzOw2STvr3Y4JWCBpfb0bMUHN/hpq2f717n5KuQcIfgAAtbLT3Y+tdyPGy8xuaub2S83/Gqaq/aS9AABAUAh+AABAUAh+AAC18s16N2CCmr39UvO/hilpPys8AwCAoNDzAwAAgkLwAwAAgkLwAwCYEDPb38yuNbM7zOx2M3tvvdtUSTVtNbOTzGyzma2Obx+rR1urZWYdZnaDmf01fk2fqHebRlJNe83sHDN7IvU9+MdatoF1fgAAEzUk6f+5+yoz203SzWZ2tbvfUe+GlVFtW3/v7i+vQ/vGo1/SEnffZmbtkv5gZr9y95X1blgF1bb3J+7+rsloAD0/AIAJcfe17r4qLm+VdKek7vq2qrxmamu1PLItvtse3xp2NlMjtJfgBwBQM2a2SNLRkv5c56aMapS2PjdOy/zKzI6Y2paNnZm1mtlqSeskXe3uDf3+V9neV5vZLWb2MzPbv5bPT/ADAKgJM5sj6eeSznX3LfVuz0hGaesqSQe6+zMkfVXS/05x88bM3Yfd/ShJCyU928yeVucmjaiK9v6fpEXufqSkqyV9r5bPT/ADAJiweOzGzyVd5O6X1rs9Ixmtre6+pZiWcfdfSmo3swVT3MxxcfdeSddKKruhZ6Op1F533+Du/fHdb0t6Zi2fl+AHADAhZmaSviPpTnf/r3q3ZyTVtNXM9onrycyereh/5Yapa+XYmNmeZtYVlzslnSzprro2agTVtNfM9k3dfYWisVk1w2wvAMBELZb0Rkm3xuM4JOkjca9JoynbVkkHSJK7f13SayT9s5kNSeqTdKY39nYI+0r6npm1KgrUfuruV9S5TSMp214z+6Skm9z9cknvMbNXKJqdt1HSObVsANtbAACAoJD2AgAAQSH4AQAAQSH4AQAAQSH4AQAAQSH4AQAAQSH4AQAAQSH4AQBgBGa2bfRa9Wdmx5rZV+LySWb2vHq3qVGxyCEAAA3GzNrcfWgs57j7TZJuiu+eJGmbpD/WuGnTAj0/AABUwSKfN7PbzOxWMzsjPt5iZheY2V1mdrWZ/dLMXjPCdV4W173ZzL5iZlfEx88zsx+Y2QpJPzCzRWb2ezNbFd+eF9e72MxOTV3vQjN7Tdzbc0W8W/3bJb3PzFab2fPN7IF4TzOZ2dz0/RDR8wMAQHVOl3SUpGdIWiDpRjO7XtGWGYskHS5pL0X7UH233AXMrEPSNySd4O4PmNmPc1UOl3S8u/eZ2SxJJ7v7TjM7WNKPJR0r6SeSXifpSjObIemFkv5Z0nMkyd3XmNnXJW1z9/+Mn/c6Sacq2qH+TEmXuvvgRN+QZkXPDwAA1Tle0o/dfdjdH5f0O0nPio9f4u4Fd39M0S7llRwq6X53fyC+nw9+Lnf3vrjcLulbZnarpEsUBUaS9CtJLzCzmZJeKun61DmVfFvSm+PymyX9zyj1pzV6fgAAaBzbU+X3SXpcUU9Ti6SdkhT3BF0n6SWSzpB08WgXdfcVcRrtJEmt7n5bbZvdXOj5AQCgOr+XdIaZtZrZnpJOkHSDpBWSXh2P/dlb0WDjSu6WdFA8LkeKgpdK5kla6+4FRTvRt6Ye+4miHpznS7qqzLlbJe2WO/Z9ST9S4L0+EsEPAADVukzSLZL+Kmm5pA/Gaa6fS3pE0h2SfihplaTN5S4Qp6feIekqM7tZUZBStq6kCyS9ycz+qihdlu4V+o2kEyX91t0Hypz7f5JeVRzwHB+7SNJ87ZpqC465e73bAABAUzOzOe6+zcz2UNQbtDgOjEaqa5L+W9K97v7FKWjjaySd5u5vnOznanSM+QEAYOKuMLMuSTMkfapS4BN7q5m9Ka77F0WzvyaVmX1V0eDol032czUDen4AAJgEZnaZpCflDn/I3X9dj/YgQfADAACCwoBnAABqIL8HmJmdY2bnp+6fa2Znx+XPx6s832Jml8Ups/z1DoxXdl5tZreb2dsrPO9P4jqrzWyNma2Ojz/dzC6s4UucNgh+AACYZGbWJuktiqaaS9LVkp7m7kdKukfSh8uctlbSc939KEWrNy81s/3yldz9DHc/Kq73c0mXxsdvlbTQzA6o8ctpegQ/AABMviWSVhU3K3X336Q2Ll0paWH+BHcfcPf++O5MjfI/O5499jplp7L/n6LtLJBC8AMAwORbLOnmCo+9RdGWFbsws/3N7BZJD0v6rLs/OsJzPF/S4+5+b+rYTfFxpBD8AAAweYqzivaV9ET+QTP7qKQhRQsQ7nqy+8NxauwpihY83HuE53q9dl3AcJ2kXVJloSP4AQCgNvriXdaLdpe0vviYpI50ZTM7R9LLJb3BR5l6Hff43KYKvTjxmKLTFW17kdYRPzdSCH4AAKiN30k6S5LMrFPR+JviDu93Kuq9Ufz4KZI+KOkV7r4jdbzbzK6Jywvj68jM5ivaPf7uCs/9Ikl3ufsjueOHKAqakELwAwBAbbxX0unxVPOVki5x9+vjx36laCPUovMVbTx6dTxF/evx8X0VpcEk6TBJf4739vqdpP+MZ3DJzL5tZsemrnemyu/Z9QJJV074lU0zLHIIAMAUiFd8/mBuQHK+zrskPeTul9fg+WYqCpqOT80sgwh+AACYEmb2VEl7p3qDJvv5DpbU7e7XTcXzNROCHwAAEBTG/AAAgKAQ/AAAgKAQ/AAAgKAQ/AAAgKAQ/AAAgKD8f+2SBOant7X+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x720 with 9 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = utils.pairplot(samples, points, limits=[[0,2000.],[-4.,0.],[2.0,3.7]], fig_size=(10,10), \\\n",
    "                           labels = ['Tint \\n U(0,2000.)', 'log_kappa_IR \\n U(-4.,0.)', 'log_gravity \\n U(2,3.7)'],\\\n",
    "                           title= 'SNRE \\n BU_new prior :' + str(1) + 'M simulations, '+ str(10) + 'k samples',\\\n",
    "                          points_colors = ['red'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Saving the samples file\n",
    "\n",
    "df_samples = pd.DataFrame(samples.numpy())\n",
    "df_samples.to_csv('/home/mvasist/samples/10ksamples__resnet_e113_SBI_1Msim_gammaTequG_modif.csv',mode='a', header=False)\n",
    "\n",
    "df_lnprob = pd.DataFrame(log_probability.numpy())\n",
    "df_lnprob.to_csv('/home/mvasist/samples/10klnprob__resnet_e113_SBI_1Msim_gammaTequG_modif.csv',mode='a', header=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# resnet , mlp, linear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# saved samples \n",
    "ss=[]\n",
    "\n",
    "#100klnprob_e76__SBI_200ksim.csv\n",
    "\n",
    "for i in ['10']:  #'100', '50', \n",
    "    s= pd.read_csv('/home/mvasist/samples/'+ i +'ksamples_e96__SBI_100ksim_gammaTequG.csv')   #waiting for this to be saved \n",
    "    ss.append(s.values)\n",
    "\n",
    "sss= np.vstack(ss)\n",
    "samples = pd.DataFrame(sss)\n",
    "\n",
    "#10ksamples_e96__SBI_100ksim_gammaTequG.csv\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = torch.tensor(np.array(samples),dtype=torch.float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([9999, 4])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "samples.size()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = torch.cat([torch.ones(1)*1.5,torch.ones(1)*2000.,torch.ones(1)*2.58])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAKmCAYAAAChcvGUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABIOUlEQVR4nO3deZxcVZ3+8edJOksTSJpdTGACgiiioCKCUWAiAgKyiQOOOiCOjuuAzgzGcX6O69iOM+q4zbghqAiKgIMg2xgWibLvuwEi0AQChCaEdJbu/v7+uLe6bhVV3Z2kqqu6z+f9etWrb90699ap6qr0k7Pc44gQAABAiia1ugIAAACtQhACAADJIggBAIBkEYQAAECyCEIAACBZBCEAAJAsghAA2T7K9tW2l9nus/1n27+2fUihzAG2w3a/7ZfWOMejtk8v3J+bly/dBmw/bvtM29vXOD6GuR3VrNcOIG0dra4AgNay/feS/kvSaZK+Kul5SS+RdJik+ZIuqTpksqTPSzp+lE/xZUkXSJoqaR9J/yrp5bZfHxHrqsqeLul7Nc5x3yifCwDWC0EIwD9K+nVEvK+wb6GkH9iu1Wp8maS/sv3liLhtFOd/MCKuzbevtj1F0hclvVbStVVlewplAaDp6BoDsIWkx2s9EBGDNXZ/W9JSZWFmQ9yc/9xhA48HgIYhCAG4XtIJtv+p1tifGvqUhaDDbe+zAc83N//5QI3HbLuj+rYBzwEAo0IQAvBBSYsl/buk+2w/Zfss2wcNc8wPJT0o6UujOP+kPNBsYnu+pE9LOjcibqpR9p8lrau+2d5qPV4PAIwaQQhIXETcL+nVkvZXFmxulXS0pEtt/0udY9ZJ+qyk+bYPHOEpvqcs0Dwv6XeSnpD07jplT5P0uhq33tG+HgBYHwQhAIqIgYi4OiL+JSIOlLSTpDsk/avtzescdqakuzTyWKEvKgsz+ysbX/QaSd+tU3ZpRNxY49a/3i8KAEaBIATgBSLiMWXdXx2SdqlTZlDS/5P0ettHDnO6P+dh5uqI+JikMyS91/beja43AKwvghCQONvb1XnoZfnPmjPKJCkizpd0g6QvSPIon3KBsgHX/zraOgJAszAbA8Cdtv9P0m8lPSRppqRDlQ2i/mVEPDzC8Z9Wdm2hUYmIx21/R9I/2n5t1aDp2XVmov05IpaO9jkAYLRoEQLwaUmdyq4WfZmkX0jaV1nLzXtGOjgiLpd05Xo+51ckPSfpM1X7T5T0xxq3d63n+QFgVBwRra4DAABAS9AiBAAAkkUQAgAAySIIAQCAZBGEAABAsghCAAAgWQQhoI3Z3s32abYfsr3a9krbt9n+mu2dN/CcV9q+snD/ANth+4B6Zdbj3EfZ/sSG1KtwjhPz+swt7Fti+/SNOe8G1OMA25+1zb+TwATGBRWBNmX7eGXLUdwpqVvSfZKmStpL0t9Kequkl2/AqT/cqDrWcJSkAyV9rcHnPVrSigafcyQHKLv69RclDY7xcwMYIwQhoA3Zfpmk0yX9RtJxETFQePgy21+VdNKGnDsi7t74Go6tiLhlpDK2p0XEmrGoD4CJgyZfoD2douz7+ZGqECRJioh1EfG90n3bc/PupBOL5Rrc7bWr7fNt99rus32t7UMKj58u6QRly2REflsywjl3sn2R7VW2n7T9X5Km1ShX0TVW6D7bz/Y5tnslXZc/1mH7U7bvtb3G9mO2/9P29KpzzrDdbfuBvNzjts+1va3tz6q8Ftq60utZ3/cMQPujRQhoT2+WdENEPNHqikiS7RdLukbZshgflfSspI9Iusj24RFxsbKFV7eW9DpJR+SH1m2hsT1V0uXKlvf4iKRlkv5O0jHrUbUzJZ0l6ViV/z37maS3KVvG4w/Kug+/IGmupLdXPfceyrodr5U0S9LBkjaX9ENJcyS9T9IbJb0gjAKYGAhCQHuaI+nG6p22J6uwyntE9I9RfT6hLCDsGxGL87r8VtLdkr4k6eKIeMD2k5LWRsS1ozjnCZJ2ys95bX7OiyXdsR71+lVEnFq6Y/tNko6TdEJE/CTf/X+2l0v6me09I+JWSe9Wtp7akRFxQfF8hXM9mm9eN4bvM4AxRtcYML48L2ld6bahM8c2wH6Sri2FIEnKu+zOkrSn7ZkbcM59JT1SDE0RMSjpl+txjvOr7h8iaa2kX+VdZB22O5QtJlt6HZJ0kKTHq0IQgAQRhID29KikHWrsf4OyrqfPjW11tIWkpTX2P66shWrzDTjndpJqdf2tT3dgdZ22UTazriIwKut2k6QtCz971uN5AExQdI0B7WmhpPfa3iYiSn/EFRE3S5Lt3avKr85/Tq3av6UaY7mkF9XY/yJJIemZDTjnUkmvqLF/2/U4R/UA5qeVvRdvqlP+sfznU5Kq30MACaJFCGhP31D2R/47+bigkTyhbGBy9R/3wxpUn6sk7VN1kcPJysbj3BIRpWv8rFE2+Hk0/ihpe9v7FM45SdJfbUQ9L5E0XdKsiLixxq0UhC6T9CLbbxvmXKWB3qN9PQDGIVqEgDYUEffYPknSjyVdZ/sHyi6oOFnSjpI+oKzLZ3VePmz/QtL7bN+flz1M2UUBG+Hrkk6UdLntf1V2ccMPS3qpKsPW3ZK2sP0hZYO9V0dEvcHPZ0haIOk82/+srPvqg5I2ZLyRJCkirrR9lrIxQl+TdL2yiyHOlXSopE9GxP3KZpa9X9JZtr+sbOr9ZspmjX0jIu7NX4sk/UM+iHsgIl4wgB3A+EYQAtpURJxp+3ZlM7b+WVk31DpJD0j6P0l/HRGPFg45WVkr72fzn7+U9DFJFzagLo/ZfqOyKen/rexaP7dKOiwiLikU/aGkfST9m6QuSX9WFkJqnXOt7bdI+rak7yob1/NzSRdJ+p+NqO67lb3ukyR9WlnLzhJJlyoffxQR62wfpOxaQR/Ifz4taZGybkApe9++qyzwfUbZWKihGXsAJgZHcI0wAACQJsYIAQCAZBGEAABAsghCAAAgWQQhAACQLIIQxpztOba/ZfuP+arjUbw+TVXZ6ba/antpvuL5H23vV6PcpHzF8SW2V9u+zfbba5QL219swmt6wSrvY6V6Zfb1OO4o25+osb9lr2VD2P4325fZfjqv94nDlH1/YVX6+2x/sEaZJbZ/1tRKtwHbJw733QNSQRBCK+ys7KJ5z0j6/Qhlf6Tsei+fkXS4sqsRX2p7z6pyX1A2bfzbkt6qbDXxc2wf2rBaD+9mZWtn3TxGz9cIRymbml9tvL2Wjym76OGwlwmw/X5J35N0rrI1yc6R9N38mkcAEsV1hNAKV0fEtpJk+2+VLYD5Arb3kPTXkk6KiB/n+66SdJekz0s6It+3jaR/lNQdEf+RH35FviBpt6TfNvG1SJLyKyuPZsX1tjcOX8usiBjMf99/U6tAvvDqlyT9NCI+ne++wvaLJX3B9g8jYt0Y1RdAG6FFCGMuX2F8NI5QdgHBXxSO7Zd0tqSDbU/Ldx+sbI2t6u6Mn0l6pe0d6z2B7U1s/ybvettjmHIvtX2+7WV519vDts/J/8DW7E6yfaXta2wfYvvWvGvvFtuvz1dF/7f8eZfbPt32jMKxNbunRtOdYXtr29+zfX/e9fiI7Z/bnl0oc7qkEyTNzs8XtpcM81ps++N5d9LavN7fdtWq86WuR9t/b/sh28/Zvsr2K6rKHWz7D7aftb0yP+9n6r2m4Yzy87SvpK31ws/IT5Wtx/bGegfanmz7+7ZX2D5wmHIvsn2G7cfyrrelti/Mg3qpzOds35yf6ynbC11YYiQvU3r/j8p/j8tt99r+Rl6X1+Wfq+dt32X74KrjT7f9qO032L4h/7wusf2xUbxPsv0BZ13Lq/M6/sj2FlVlTrZ9T/6Zfsb2jbaPHs35gXZDixDa2SskPRQRq6r236Us+Oycb79C2dWDF9coJ0m7SXqo+uT5P+4XStpK0hsi4gVlCi5S1pX3IWULds5WtmTDSP+Z2FnSV5W1RqyU9O+SLshvHcqWrXh5XmaZpFNHON9obKFs6Y1PSXpS0osl/YOkRbZfFhGrlXUlbq1sJfsj8uPW1DhXyZfy831H0m+UvadfkLSH7f2rwsi7lS3xcbKy39NXJf1v/tz9tndS9vp/paxlb62kXSTttLEvfBilIHZn1f7iZ+SK6oNsd0o6S1mQOqC06G0dP5X0F5L+SdIjyhaPfbOkTQplZitbruRRSTOUvVdX235tjaVIviHpPGXrue0n6V+ULbFyoLL3tCffd57tv4iIpwrHzlT2H4ivKPteHC/pm7afi4jT670A293KPivfzF/HbElflLS77TdExIDtd0n6T2W/u98r65Z8lbLPHTDuEITQzrZQ7VXNlxceL/3sjRdeJr263BDbOyhbcmGlpHkR8WS9StjeSlmgOTIiLig89PMRX0HW2vCGiHgwP9ckSf8raceIKLUuXOpsAPg71IAgFBGlEFKq/2RlS0c8rGz81PkR8YDtJyWtjYhhu8HywPgPks6IiI8W6vyksj/+hysLNiXrJB1e6mqyLWXjcfaW9AdJr1EWkD5UWKx14Ya/4lEpfQaqP0/DfUY2Vxb6tlP2O3xghOfYV9I/R8SZhX3nFAtExN8Wzj9Z2SKxd0n6WxV+Z7mFEVEaw3W57cMkfVTSmyLimvwcSyXdpmy9tzMKx24m6QMRcXZ+/5K8RfBzts+o8V1R3sr4T5I+FxGfL+y/X9I1kt4m6df567y9WEZj0P0MNAtdY0jRbsr+ID8i6S+HC0G5pyU9KKnb2ayjXdbjue4vhaDcvfnPS6vK3StpjvPUsLFsfyjv3lgpqV9ZCJKkXTfgdPuodtfj2fm596/af3nVeJtSS8cO+c9blYWls20fW+w6aiMvVvbHfxONLgRJ0g2S/invNnplrd+l7QNtX2H7aWXv3TplC9fW+r1cXHX/XknPl0JQYZ8kbV9VdkDZoPCis5X9Dmartrco+5twprOu2w5nXb/XSXpOWatU6XXu6Wzm54G2N6lzPmBcIAihnT0jafMa+0v/e19eKNdV4w9PdbmS/ZT9MfhRRKwcqRL5/57fomw19S9Lut/2gx7dbKPqFoi1w+zvUNb1sVHysSDfVbYw6zHKWmJK41Cmb8ApS+/j0uLOfLzW03pha0r1+13qcpueH7dY2biuScpalB63fa3t6kDVSKX3u/rzVO8z8iplgfkXEfHEKJ/jOGUtY6dKul1Sj+3P5K2Asv0aZS0nKyW9T9nv5HXKWnRq/V5qfUZ6izsiovR5qj7+mRqDv0uvo14QKgXSxcoCWvG2mbLWTUn6ibIu4tcrC/TLbZ9npuFjnKJrDO3sLklH296kapzQbsr+KCwulJsm6SWqHCe0W/7z7qrzfk/SLEk/td0fEdX/c36BvFXnb/KwtYeyLorv2l4SEdX/c2+E1fnPqVX7t6wuWMPxkn4XEf9Q2uFhBoyPQikkvEjlMTWlmVhb6oUhYkQRcYWyWVvTJM1TNt7kIttzq8a6NEqp3q9QZaCr9xm5RFlA+Yrt1RHxXyM9QUQsk/QRSR+xvauyweifUzZO678lvV1ZK9AxxZCSd8H1ru8LGsHmtqdUhaFt8589dY55Ov95kGp3ST8tDf3H4HuSvpfX/SBlY4Z+oSwcAeMKLUJoZ7+RNEXZ2BlJQ398j5N0WUSUWhouUfa/1ndVHf9uSXfWGAQd+ViX7yrrnnmHRikyt6p8/Z3dR3vsevpznfMfNopjN1H2fhS9t0a5NcoGuo7kWmXB8/iq/ccp+8/UlaM4R00RsSYiFiobRD5D0sYEtuH8Udkg91qfkeXKxlBV1+2ryi7L8A3bH1+fJ4uI+yLin5UFitLvcBNlXVZD43Nsz1e5y7CRJisLXkXHK+sirReELpc0KGmHiLixxu0Fkwki4pmI+IWkX6p53wWgqWgRQkvYPjbffG3+86354NsnI+IqSYqIW2z/QtkfoinKZn59SNkfy6E/aBGxzPbXJH3K9nPKLgR4nKT5Ks+IeoGIOMX2gKSf256U/4Neq66vkvRfyv7Hu1jZH5kTlf3vvimDfCNiqbNrJn3K9lPKZpS9W6ObWXWJpE/a/mdJ1yt7H46tUe5uSVvkXXw3SlpdY+aSImK57f/M6/K8su6dlyubTXSNshl1o+bsas775ed5RNmsvU9Jekz5rK68m+UhZQN3PzvC+fZXNgPuRfmuvfKxUYqIX+U/19n+f8pa8XqUdRvOl3SSpI8VupiqX/vX8s/I1/PPyH/WqcOs/JxnKhu3s07Skcq64i7Li10i6RRJp9v+sbKxQf9P9YPJxnhO0r/nA/3/JOmdymabnVhroLQk5QPovyLp23mL1lXKWia3V9Y1/MOIuML29/Pz/1HZ5/Klkt5TeJ3A+BIR3LiN+U3Z/4pr3a6sKtcp6WuSHlf2j/J1yqYxV59vsrKpxH9W1tJxu6Rj6zzvF6v2fUVZqHlnnbpuo2xGzv2SVilrQbhK0sGFMgfk5z6gsO9KSddUnWtuXu5vq/Z/Nt/fUdg3R1mrWG/++v9N2eyikDS3UG6JpNOr3rP/VtYl85yySwTsmB/32UK5Gcqmhj+TP7ZkmNdiSR9XNi1+rbLupe9ImjmK97f0mk/M7++rbObcI/nvaqmy2VW7Fo55RX7MB0fxWbqy3uepRtm/y3+Pa5QFhA/XKLNE0s+q9n1EWWvJqXXqME1Zd9FdysYArVA2qPivq8p9TFnA68sfPzCv/5WFMqX3/8CqY0+X9OhIn+lSOUlvyJ9jtbLvxd9XHXdi9Wcp3/8eZa2Az+ev5R5lV2yfkz9+Ql7nZfn7+JCySwLMrPXecOPW7jdH1PzPAQC0jO0PKLt20V/EC68jhWE4u1jmgRExp9V1AcYDxggBaEf7S/o6IQhAszFGCEDbiYjqQc0A0BR0jQEAgGTRNQYAAJJFEAIAAMkiCAEAgGQRhAAAQLIIQgAAIFkEIQAAkKz1uY4Q8+zRTtzqCrxl0jv4TqCtXD54Tku/F3wn0G5G852gRQgAACSLIAQAAJJFEAIAAMkiCAEAgGQRhAAAQLIIQgAAIFkEobF0yinZDQAAtIX1uY4QNtatt7a6BgAAoIAWIQAAkCyCEAAASBZBCAAAJIsgBAAAkkUQAgAAySIIAQCAZBGEAABAsghCLTKve6HmdS9sdTUAAEgaF1RskZ7evlZXAQCA5NEiBAAAkkUQAgAAySIItRjjhAAAaB2CUIsxVggAgNYhCAEAgGQRhAAAQLIIQgAAIFkEIQAAkCyCEAAASBZBCAAAJIsgBAAAkkUQAgAAySIItdjsrk6uLg0AQIuw+nwLLek+TJI0d8FFLa4JAABpokUIAAAkiyAEAACSRRACAADJIggBAIBkEYQAAECyCEIAACBZBCEAAJAsghAAAEgWQQgAACSLIAQAAJJFEAIAAMkiCLXAvO6Fmt3V2epqAACQPBZdbYGe3r6hBVcBAEDr0CIEAACSRRACAADJIgi1gdldnZrXvbDV1QAAIDkEoTawaMF89fT2tboaAAAkhyAEAACSRRACAADJIggBAIBkEYQAAECyCEIAACBZBCEAAJAsghAAAEgWQQgAACSLIAQAAJJFEAIAAMkiCAEAgGQRhAAAQLIIQgAAIFkEoTYxu6tT87oXtroaAIB67PINEwZBqE0sWjBfPb19ra4GAABJIQgBAIBkdbS6AgAmsGIXQkTr6gE0QuEz7I7KP5/R31++M1zXWZ3vQfF8Feca7nxV5xr2HHWM+phRfpfrnq/6NbTRvwe0CI2xmx/u1eyuzlZXAwAAiCA05tb2D2jRgvmtrgYAABBdYwA2kqdMHdqetOmMysemTxvajtVryttrytuDfVWTBNqoyRwodvVM2nzzoe1YvbqyYOHzPalrVnn/4EBluWLX0bMrys+z2WajOkaF7ibPmllRLJ55tvzYzE0LdV1TUU6Dhe/YVl3l7WVPl4+fPl31xKryd/YF3WkD5bpXvF/PPVd5joHBOicv7x9t997GokUIAAAkiyAEAACSRRACAADJYowQgPU3afLQpqdOKe/fbuuKYs/uVh4j8Ozc8jGb/6nc97/ZHcsqjhl8sjxOYbBqXAEw1orjXFbO27G8f03lWLZ1m5U/31GYKe6qIW9TVpbH0Dz/otp/gtd0VU41d2HI0GaPlu/07jy5otz0p7areb7BqZX3B6aXz99fmMS86SNbDW2vnVVZh2m95bE7U1eWt101nGnqM2uHtif1F8r1b1tZhxnlfzc6nimPOZq0sjz2qv/BJRoLtAgBAIBkEYQAAECy6BoDsP6K03sHBuoWW75ruen+bcf+YWj75Z2PDW1/7pojKo7Z+YxyV8SUu/9cfprlz1SenGn2GAOrXjd3aPuYL102tH19744V5R5d2TW0/aotyp/vSa6cJr799OVD279/epfyMbN6hra3nbKi4pila8vT8e99rtzFdMDMxyvKregvT3l/dFW5Ph948dUV5a5+bteh7W2mlp9r2drydPxf3bdnxTGTZpS7rNxR/s4//sgWFeW8pnDJjCnl7+i0ZZXdeJPWlrveZj1YPubJ15bLvOQflmgs0CIEAACSRYsQgDG32cOrtftpPXrJBU/pxOf/oOenT9P5+7xaPzj4TZK6Wl09oC3MeniVXn36w9r1N09o2qp+rdmkQwvfsqsuOH4PPT571sgnwKgQhABslMG164a2JxeuTCtJW9/eNbR9zZt2kiTtc8ODOvZLt0jrQs4nj222eo3++prrdNy1N+jkQ0/SNTu9XJK0/UOFKS0uXzVXkhT1u+SARllbmA2227Ry99UpOy6pKPe5J3cb2p7VsWpoe6eplbMif1/olvrbQpfV4/1dQ9s/XrKv9r3hAX35i1epo39QU/KrME9/vl9vvfAuHXLpXVr+vc31qZcfXnHupX3lcDR/q3uHtrsmraost6Zcbl2UX19PX7kOe23/SMUxty6dPbS9anW5W8vrqjqWNitcDXp1YSZdZc+YVu9cvtr16m0LUaQFPd50jQEYM7Mfe0b/9oVfy33lEFTS0R+atnpAX7/gdG3f+1RrKgi0gdmPPaMvf/HX6lzTPxSCStwvTeqTtvi7Z7T1o1xeohEIQgDGzDvPvUEd/XXWGMp1DA7oPTdeNUY1AtrPX5838vfE66QDz7x7jGo0sRGEAIyZQxbe/YL/4VabMjiot9190xjVCGg/hyy8a8TvifulfS9+aIxqNLExRgjAxilMpR9Y3lvx0Ixry/9Qr1sxR52r1mo0Zqxdoxf94XnF86tGLgw0UdctTw5tf+az7xva/sTsynaEKc+Xtz1Qf6BLR/kiyrps832HtgcLf40/sKpyuns9057v163f2rPmY+c/s8PQ9q+mVV4let0m5bpP6SsHrsmry/XuWFU5Bm/LWeUKbrK0/CIGpq+rKNe3Tfky1rPuLF8qYPWcmRXlil3jk9aWxwut22yKxhotQgDGTN/kaSMXktQ3eerIhYAJatWU0X1PRlsOwyMIARgz/7ftnlo3wj876zRJl2+1xxjVCGg/v93ltVrnEb4nnqTf7vTaYctgdOgaA9A4g5XN6QNPlrsVJl39lM6NLXSwrOEavwc0Sf+7ekd13PuwBlcVusZi+DETQDMMPli+uvmWfeUunC0HKz+P0V/u6/H0QkvN2squo9h0k/KdFSsLT1Tulvrtqu10REzSFNX/zA/EJF38+BxteWF5mnysKddPhfq5s7N4aEVdta5Qv0nl8OWOyngwpfBdLD5WHdc6CueOqeWW3al/qpomWvw+F0Lf5ML+sZpJT4sQgDGz1Jvq895Xq9Whdaoat6BJWu0OfXHTA7V08sw6ZwAmvqWTZuoLnX9Z53tirVaHvjiT70mj0CIEYEzd4O30wc636Zi1d+vN/Q+qU+vU56la2Lmrzt90Tz1WWC8JSNWNU+bog5OP1DFr7sq+J7FWfZ6q303dWedv8qosBA2y3l4jEIQAjI3CIqk9fZP1Lb1S3/Iry4/3SerryzdyLvxvmEVW0QLFbqT+R3vqlps0rbDY6HPlzpaK7l2p/me6xv4eWd/S7vp2x57lh9b2S2v7JS2vPKb6fEWrV9feP1x9hlHRtdaAcq1G1xgAAEgWQQgAACSLrjEA7YvuMLSTYT6Pg8Xup+ouq9GcY5hz1+1iasT3g+8YLUIAACBdBCEAAJAsghAAAEgWY4QAAGgkxt2MK7QIAQCAZBGEAABAsghCbWR2V6fmdS9sdTUAAEgGQaiNLFowXz29fSMXBAAADUEQAgAAySIIAQCAZBGEAABAsghCY+jmh3s1tWNyq6sBAAByBKExtLZ/QK/ZoavV1QAAADmCEAAASBZBCAAAJIsgBAAAkkUQAgAAySIIAQCAZBGEAABAsghCAAAgWQQhAACQLIIQAABIFkGozczu6tS87oWtrgYAAEkgCLWZRQvmq6e3r9XVAAAgCQQhAACQLIIQAABIFkEIAAAkiyAEAACSRRACAADJIggBAIBkEYQAAECyCEIAACBZBCEAAJAsghAAAEgWQQgAACSLIAQAAJJFEBoj87oXamrH5FZXAwAAFBCExkhPb59es0NXq6sBAAAKCEIAACBZBCEAAJAsghAAAEgWQQgAACSLIAQAAJJFEAIAAMkiCLWh2V2dmte9sNXVAABgwiMItaFFC+arp7ev1dUAAGDCIwgBAIBkEYQAAECyCEIAACBZBCEAAJAsghAAAEgWQQgAACSLIDQG5nUv1OyuzlZXAwAAVCEIjYGe3j4tWjC/1dUAAABVCEIAACBZBCEAAJAsglAbm9e9kDXHAABooo5WVwD1sd4YAADNRYsQAABIFkGoTc3u6mTKPQAATUbXWJsqTbefu+CiFtcEAICJixYhAACQLIIQAABIFkEIAAAkizFCTTSve6F6evsY9AwAQJsiCDVRT2+flnQf1upqAACAOugaAwAAySIIAQCAZBGEmmRe90LGBgEA0OYIQk3S09s3dFHEjTG7q5OFVwEAaBKCUJtbtGA+i68CANAkBKEmaPSyGLQKAQDQHI6IVtcBAACgJWgRAgAAySIIAQCAZBGEAABAsghCAAAgWaNea8z2nZJWN7EurbSVpKdaXYkmmaivbXpE7N7qSgAAxrf1WXR1dUTs1bSatJDtG3lt44vtG1tdBwDA+EfXGAAASBZBCAAAJGt9gtD3m1aL1uO1jT8T9XUBAMYQV5YGAADJomsMAAAkiyAEAACSNaogZPsQ2/fZXmx7QbMrNVZsn2Z7WX6NpAnD9va2r7B9t+27bJ/c6jo1iu3ptq+3fVv+2j7X6joBAMavEccI2Z4s6X5Jb5H0qKQbJL0zIu5ufvWay/Z+klZK+slEujif7e0kbRcRN9veTNJNko6aIL8zS5oRESttT5F0jaSTI+LaFlcNADAOjaZFaG9JiyPiwYhYK+lsSUc2t1pjIyKulrS81fVotIhYGhE359vPSbpH0uzW1qoxIrMyvzslvzHiHwCwQUYThGZLeqRw/1FNkD+qKbA9V9KrJV3X4qo0jO3Jtm+VtEzS5RExYV4bAGBsMVh6ArO9qaRzJZ0SEStaXZ9GiYiBiNhT0hxJe9ueMN2aAICxNZog1CNp+8L9Ofk+tLF8/My5ks6MiPNaXZ9miIheSVdIOqTFVQEAjFOjCUI3SNrF9o62p0o6XtIFza0WNkY+oPhHku6JiK+1uj6NZHtr2135dqeyQfz3trRSAIBxa8QgFBH9kj4q6VJlg25/GRF3NbtiY8H2WZL+KGlX24/afl+r69Qg8yS9R9J827fmt0NbXakG2U7SFbZvVxbSL4+IC1tcJwDAOMUSGwAAIFkMlgYAAMkiCAEAgGQRhAAAQLIIQgAAIFkEIQAAkCyCEAAASFbSQcj2ypFLAQCAiSrpIAQAANJGEFK2JIXtr9q+0/Ydto/L90+y/V3b99q+3PZvbR87zHkOzcveZPubti/M9+9t+4+2b7H9B9u75vtPtP3r/NxLbH/U9ifyctfa3iIvd6Xtr9u+0fY9tl9n+zzbf7L9xcLz/zp/7rtsf6C57xoAAONfR6sr0CaOkbSnpD0kbSXpBttXK1uqYq6k3SRto2yJkdNqncD2dEnfk7RfRDyUL99Rcq+kN0VEv+0DJf2bpLfnj+0u6dWSpktaLOmTEfFq21+X9DeSvpGXWxsRe9k+WdL/SnqtpOWSHrD99Yh4WtJJEbE8X4PrBtvn5vsBAEANBKHMGyWdFREDkp6wfZWk1+X7z4mIQUmP275imHO8TNKDEfFQfv8sSaVWmVmSzrC9i6SQNKVw3BUR8Zyk52w/K+k3+f47JL2qUO6Cwv67ImKpJNl+UNL2kp6W9Pe2j87LbS9pl3w/AACogSA0Nr6gLPAcbXuupCsLj60pbA8W7g+q8vezpkaZoXK2D5B0oKR9I2KV7SuVtTKheVioD+3Era7AWya9g+8E2srlg+eM+L1gjFDm95KOsz3Z9taS9pN0vaRFkt6ejxXaVtIBw5zjPkk75UFHko4rPDZLUk++fWID6100S9IzeQh6maR9mvQ8AABMGAShzPmSbpd0m6SFkk6NiMclnSvpUUl3S/qZpJslPVvrBBHRJ+nDki6xfZOk5wpl/13Sl23foua1wl2irGXoHkndkq5t0vMAADBhOIKWzOHY3jQiVtreUlkr0bw8JA1X1pK+I+lPEfH1sawvxhRfHrQTusaAKqPpGmOM0MgutN0laaqkL9QLQbn32z4hL3uLsllkAACgTRGERhARB1Tvs32+pB2rdn8yb/2hBQgAgHGCILQBIuLokUsBAIB2x2DpFsmvFr1Xvm3bC23PzO8fYvs+24ttLxjFud5l+/b8qth/sL1HnXJn59cyAtrDKadkNwBoEYLQBrDd6Ja0QyXdFhErbE9WNtD6rcquaP1O27uNcPxDkvaPiFcqu2bR9+uU+29JpzaozsDGu/XW7AYALZJs11i+8vx/STpcUp+kIyPiifw6QKcpW2rjSUnvjYiHbZ8uabWy5TAW5euA9eX3t5F0krIlMfaVdF1EnJg/z38ru0p1p6RfRcS/1qjOu1QOL3tLWhwRD+bHny3pSGVT+GuKiD8U7l4raU6dor+XdLrtjojor3c+AABSkXKL0AxJ10bEHpKulvT+fP+3JJ0REa+SdKakbxaOmSPpDRHxifz+5sqCz8eVLYHxdUmvkPRK23vmZT4dEXspWy5jf9vFZTNK5km6Kd+eLemRwmOP5vtG632SLq71QL5UyGJla6oBAJC8lIPQWkkX5ts3KVtcVcqCzc/z7Z8qW2+s5Jx8PbKS30R2IaY7JD0REXfkYeOuwvn+yvbNyqbTv0JZd1e1LfL1xjaK7b9UFoQ+OUyxZZJevLHPBQDARJBs15ikdVG+muSARvdePF91f6T1v3aU9I+SXhcRz+Tda7XW/+q3PSkPUT3KFkwtmaPy8hx15S1NP5T01hFWnJ+urEsPAIDkpdwiVM8fJB2fb79L2biaDTVTWXh6Nl+r7K11yt0naad8+wZJu9je0fbUvC4XVB9g+2jbX863d5B0nqT3RMT9I9TppZLuXO9XAgDABEQQeqGPSXqv7dslvUfSyRt6ooi4TVmX2L3KutsW1Sl6kfIFXfNBzB+VdKmkeyT9MiLukiTbn7d9RH7MSyStyLc/I2lLSd+1favtG0sntv1b2y/Ot7eV1DfC1bEBAEgGa421AdvbSfpJRLxlPY75maSPR8ST63HMxyWtiIgfbUA18UJ8eTbQvO6FkqRFl3w+23Hlla2rzMTBWmNAFdYaGyciYqntH9ieGRErRj5Cioh3b8BT9SobAA60VE8vw9QAtAeCUJuIiF+OwXP8uNnPAQDAeMIYIQAAkCyCEAAASBZBCAAAJIsgBKClbn64d2gWGQCMNYIQgJaY3dWpmx/u1dr+AWaRAWgZghCAlli0YL7W9g+MXBAAmojp8wBaZmrH5FZXAUDiCEIAWuY1O3S1ugoAEkfXGAAASBZBCAAAJIsgBAAAkkUQAgAAySIIAQCAZBGEAABAsghCAAAgWQQhAACQLIIQgDE1r3uhZnd1troaACCJK0sDGGM9vX1a0n1Yq6sBAJJoEQIAAAkjCAEYM3SLAWg3dI0BGDN0iwFoN7QIAQCAZBGEAABAsghCAAAgWQQhAACQLIIQAABIFkEIAAAkiyAEAACSRRACAADJIggBaLnZXZ2a172w1dUAkCCCEICWW7Rgvnp6+1pdDQAJIggBAIBkEYQAAECyCEIAACBZBCEAAJAsghAAAEgWQQgAACSLIAQAAJJFEAIAAMkiCAEYE/O6F2p2V2erqwEAFTpaXQEAaejp7dOS7sNaXQ0AqECLEAAASBZBCAAAJIsgBAAAksUYIQAAJgK7vB3RunqMM7QIAQCAZBGEAABAsugaAwBgfRW7oSR56tSh7Vi7tvxAdRdV1XE1y1WXoZurqWgRAgAAySIIAQCAZNE1BgDAKHjatPKdgYHKxyZPHtqe1NU1tB1r1lSUG+zrq3lMUVSdu153WrE7rtZz1T1+Y7vaJtjsNFqEAABAsghCANrC7K5Ozete2OpqAEgMQQhAW1i0YL56evtGLggADcQYIQAA6pi87TZD27HtFuXtaVMqC64rj+tZvVXn0PbUZ1ZXFJtUKOcnlpe3J5XbJWLTTSqOcV957E88v6q83d9fee5NyscVxwtVjG2SNNDbq1omdZbrXXEJAFWNWxpmXFDxuYa9jEC94zvKsaT69TULLUIAACBZBCEAAJAsusYAAKhjcIdth7bvf0+h62mTgVrFJUlTnip3m3U8N61uuXUzNx/anra8PCW9c1llN9KkQg+RB8uPTV5Xeb4ZPeVuuMEp5XaOSesGKwsO7lCu34ryMX525dB2zOisOES9K8rlppRf3+CzKyqKFbvDJm26aXl/3zDj/wqXESheUoCuMQATxrzuhZrd1TlyQQAYY7QIAWi6nt4+Lek+rNXVAMaVOSue0rvuvEqHPnCTNlm3RqumTNPFc1+jn++2v3o226rV1ZswCEIAANSx4iUzhrb/59AfDW0ftEllv9SVfeUOlikud+nc0LdTRbl7nt9uaPsvpi9XLZc9/nK9/voH9fnPX6CO/kFNGci6tjZdt0ZHPXidDn/4Rv3jgmN00TZ7VRwXHYUrTU8qd6HN3KayW+q5h2cObW/65/LrG+zYcmi7untu7azy7LmuxeXX17HqRRXlpj1R7l5Tf2GG3GBV91yBBwqPrS2/r4OrVtUo3Xh0jQEA0EZe/FivPv/5C9S5pn8oBJVMGRhU55p1+o/u87TDk0+1qIYTC0EIAIA2ctyvblRHf/0WFEnq6B/Q+xb+foxqNLERhAAAaCMH/e7uF7QEVZsyMKhjrr9pjGo0sTFGCACAOjr6yoHkkXXlMTRX9i2rKPfhW/56aHvWjPKYnC+/9PyKcqsGytPpd52+dGj7DdMfG9repK/yqs71bLpmjU576w+G7j89UJ6u3jtQnuq/fGBGxXHauby5y7QnhrZPf2ze0PYx21aGrC/eeujQ9ordy2ORJj8ztaLc4LSuoe2tbipM4e+vHHO0Zlb5sU2eKr/Hq7vKlxHY6vs9Ggu0CAFoGyy8CkhrN5k8ciFJAzP4E94IvIsA2gYLrwLSvW97kQY6PGyZ/g7r8aM2G6MaTWx0jQEAUMdmtz0+tP21nx4ztD1zSeUYnukzy8Fl+fblgPLhaz5YUW7N1uXjpi8rt0WsmlueNr7D7st16Xlf0yaqunR0wbrJHfqvow/UsifLU+FvfnzO0PZglOvz5h3ur3ue/WfcO7S99+ZLhrZPe/iNFeXeNPfBoe0nVpdf391LXlxRzpPLr+/JfctXoHZ/ZbCLjvIU/Li7HEVWl3sfNVZXSqJFCACANvLw1lvqQ3/3Hq2aOkVrJ1f+mV43eZL6pk3RPy44RsvmzKxzBqwPghAAAG3myt1fpoM/8wmd9cbX6/kZUzRo6fkZU3TuIa/WO771fi3aa+eRT4JRoWsMAIA6Bh4tz+aa++PCbK7JVQOaB8pXUY6ZhcVGOytnVVUcMqPw2FXlzb5ts5llg9pK39/qGP301eXlaZbP7pRuzFox7tAWFeebXl4/tWJx1ot3e11FuY6V5W6q32xbvjp1x3Pl/TOqJmzdvqZ8Rey+bcrlpleliCnPl7fXFSerVTW7THmuvB2Ft3Lyao05WoQAAECyCEIAACBZdI0BaCulawktWjC/1VUBFP3l2U2Dz64o7y90hUmq6Bqb9HxhsdDC8dU6Omr/Cd7s7spuN08vX4RxmzvK5/aMTSrKxeo15cc2nzW0vcX1lW0exUVOB2eWzzFpRf1FTge2LM8Um7Sq3EUYkytng/XP6hzaXrVdud4z719RUa7vxeV+s6nPlmfH9W0zTWONFiEAbYVrCQEYSwQhAACQLIIQAABIFmOEAAAYhcG+YbpsozxdfWDFivrlGmhSVC5kGmvLY3cmFcYsDa6qP/aneBmAgeK4p6rLA/jJp8vPM0y5qbPKF3mcclt5jrxnVV78ccZTz9aszpQHy+euP7qqsWgRAgAAySIIAQCAZNE1BgDAaFR1RbXa4HPPVe5weSr7wJo1GpXi9P7C8Rrm8gDFSwq46hIAA4XuuaJ4tPJS1Z5avqp2sUvP1VfsHgO0CAEAgGQRhAAAQLLoGgMAoBWKXVGN6HYrnmNDzj1MuahzhewX7C92qQ13vmLXXaGu9Z6nmWgRAgAAySIIAQCAZBGEAABAshgjBKCp5nUv1OyuzpELAqlp5nT8Vk3135DnbfFlCQhCAJqqp7dPS7oPa3U1AKAmusYAAECyCEIAACBZBCEAAJAsghAAAEgWQQgAACSLIAQAAJJFEALQNFxDCEC74zpCAJqGawgBaHe0CAEAgGQRhAAAQLIIQgAAIFkEIQAAkCyCEAAASBZBCAAAJIsgBAAAkkUQAgAAySIIAQCAZBGEALSd2V2dmte9sNXVAJAAghCAptiYdcYWLZivnt6+BtcIAF6ItcYANAXrjAEYD2gRAgAAySIIAQCAZBGEAABAsghCAAAgWQQhAACQLIIQAABIFkEIAAAkiyAEAACSRRAC0HAbc1XpEpbZADAWCEIAGq6nt0+LFszfqHOwzAaAsUAQAgAAySIIAQCAZBGEAABAsghCAAAgWQQhAACQLIIQAABIFkEIAAAkiyAEAACSRRAC0FBcDRrAeEIQAtBQjbwaNMtsAGg2ghCAhlvSfVhDzsMyGwCajSAEAACSRRACAADJIggBAIBkEYQAAECyCEIAACBZBCEAAJAsghCAtsa1hAA0E0EIQMPM616o2V2dDT0n1xIC0Ewdra4AgImjp7evYRdTBICxQIsQgIZoRmsQADQbQQjARiuN4Vm0YH5Tzs84IQDNQhACsNF6evuaFoIkxgkBaB6CEAAASBZBCMC4Ma97IV1kABqKWWMANti87oXq6e0bk0HSpeegiwxAIxGEAKy3UqvMWE6XL41Bmte9UHMXXKTZXZ1NHZcEIA2OiFbXAQAAoCUYIwQAAJJFEAIAAMkiCAEAgGQRhAAAQLKYNQZsINt3Slrd6nqMYCtJT7W6EiOgjo0xPSJ2b3UlgPGGIARsuNURsVerKzEc2zdSx403XurY6joA4xFdYwAAIFkEIQAAkCyCELDhvt/qCowCdWwM6ghMUFxZGgAAJIsWIQAAkCyCEAAASBZBCCiwfZrtZfk1gkr79rR9re1bbd9oe+98v21/0/Zi27fbfk3hmBNs/ym/nTAGddzD9h9t32H7N7ZnFh77VF7H+2wfXNh/SL5vse0FDa7j9ravsH237btsn5zv38L25fn7crntzfP9Y/peDlO/d+T3B23vVXVMO72PX7V9b/5enW+7q5X1BMa1iODGjVt+k7SfpNdIurOw7zJJb823D5V0ZWH7YkmWtI+k6/L9W0h6MP+5eb69eZPreIOk/fPtkyR9Id/eTdJtkqZJ2lHSA5Im57cHJO0kaWpeZrcG1nE7Sa/JtzeTdH9el3+XtCDfv0DSV1rxXg5Tv5dL2lXSlZL2KpRvt/fxIEkd+f6vFN7HltSTG7fxfKNFCCiIiKslLa/eLanUwjJL0mP59pGSfhKZayV12d5O0sGSLo+I5RHxjKTLJR3S5Dq+VNLV+fblkt5eqOPZEbEmIh6StFjS3vltcUQ8GBFrJZ2dl21UHZdGxM359nOS7pE0O3+OM/JiZ0g6qlDPMXsv69UvIu6JiPtqHNJW72NEXBYR/XmxayXNaWU9gfGMIASM7BRJX7X9iKT/kPSpfP9sSY8Uyj2a76u3v5nuUvkP2zskbd8udbQ9V9KrJV0naduIWJo/9LikbVtdz6r61dNu72PRScpa09qinsB4QxACRvYhSR+PiO0lfVzSj1pcn1pOkvRh2zcp60JZ2+L6SJJsbyrpXEmnRMSK4mMREcpa21pmuPq1k3r1tP1pSf2SzmxV3YDxjiAEjOwESefl2+co62aQpB6VW16krHuiZ5j9TRMR90bEQRHxWklnKRsP0tI62p6i7I/3mRFRev+eyLu8lP9c1qp61qlfPe32Psr2iZIOl/SuPFS2tJ7AeEUQAkb2mKT98+35kv6Ub18g6W/yGU/7SHo27/a5VNJBtjfPZ0UdlO9rGtvb5D8nSfoXSf9TqOPxtqfZ3lHSLpKuVza4ehfbO9qeKun4vGyj6mNlLWf3RMTXCg9doCxYKv/5v4X9Y/ZeDlO/etrqfbR9iKRTJR0REataXU9gXGv1aG1u3Nrppqw1ZamkdcrGUbxP0hsl3aRsps11kl6bl7Wk7yhrfblDlbOMTlI2UHWxpPeOQR1PVjaj6H5J3cqvGp+X/3Rex/uUz37L9x+al39A0qcbXMc3Kuv2ul3SrfntUElbSvqdsjD5f5K2aMV7OUz9js7f0zWSnpB0aZu+j4uVjfkp7fufVtaTG7fxfGOJDQAAkCy6xgAAQLIIQgAAIFkEIQAAkCyCEAAASBZBCAAAJKuj1RUAMHHYLk2Pl6QXSRqQ9GR+f+/I1rkCgLbB9HkATWH7s5JWRsR/tLouAFAPXWMAmsr2a21fZfsm25cWlth4re3b8ttXbd+Z7z/R9rcLx19o+4DW1B7AREcQAtBMlvQtScdGtg7aaZK+lD/2Y0kfi4g9WlU5AGCMEIBmmiZpd0mXZ8tmabKkpba7JHVFxNV5uZ9KemtLagggaQQhAM1kSXdFxL4VO7MgVE+/KlurpzehXgAgia4xAM21RtLWtveVJNtTbL8iInol9dp+Y17uXYVjlkja0/Yk29tL2nssK4yNEty4tdFtVGgRAtBMg5KOlfRN27OU/ZvzDUl3SXqvpNNsh6TLCscskvSQpLsl3SPp5rGsMIC0MH0eQMvZnivpwojYvdV1wUbhDwraiUdTiK4xAACQLFqEAACNwh8UtBNahAAAAIZDEAIwLNsrq+5XX/n5FNt/k29vYfty23/Kf24+yueYbPsW2xcW9u1o+zrbi23/wvbUfP+0/P7i/PG5hWM+le+/z/bB+b6ptq+2zeQQAC9AEAKwwfJwcZKkn+e7Fkj6XUTsomzx1QWjPNXJymaIFX1F0tcjYmdJz0h6X77/fZKeyfd/PS8n27tJOl7SKyQdIum7tifnC73+TtJx6/8KgTZwyinZDU1BEAKwMeZLujki+vP7R0o6I98+Q9JRI53A9hxJh0n6YWGf83P/qsa5is/xK0lvzssfKensiFgTEQ9JWqzyNYh+rcprFQHjx623Zjc0BUEIwMaYJ+mmwv1tI2Jpvv24pG1HcY5vSDpV2TWHSraU1FsIWI9Kmp1vz5b0iCTljz+blx/aX+OYOyW9bhR1AZAYghCADVGaHbSdpCdrFsimpA47i8j24ZKWRcRNw5XbWBExIGmt7c2a+TwAxh+CEICR9JUGKue2kPRU6TFVrgX2hO3tJCn/uWyEc8+TdITtJZLOljTf9s8kPS2pqzDAeY6knny7R9L2+XN0SJqVlx/aX+MYKVsAdvUI9QGQGIIQgJFcJendkmS7U9JfSboif+weSTsXyl4g6YR8+wRJ/1vrhLbvlaSI+FREzImIucoGOi+MiHfnrUlXKFueo/pcxec4Nj8m8v3H57PKdpS0i6Tr8+fbUtJTEbFug94BABMWQQjASE6WdIztWyVdK+mciLg6f+xiSfsVynZLeovtP0k6ML8v23vZ/mG+vZVGd6GzT0r6hO3FysYA/Sjf/yNJW+b7P6F8ZlpE3CXpl8rWKLtE0kfyLjFJ+ktJF63n6waQAK4sDWCj2D5f0qkR8adRlj9c0k4R8c3m1qziOc+TtCAi7h+r50wUf1Ca4YADsp9XXtnKWoxHo7qyNBcYA7CxFigbND2qIBQRF45cqnHy8U2/JgShVeZ1L1RPb59md3Vq0YL5ra4OqtA1BmCjRMR9ha6ythMRayPiJ62uB9LV09unJd2Hqae3r9VVQQ0EIQAAkCyCEAAA62le90LN6144Js9z7YNP6+aHe5v+XKlijBAAAOtprLq5enr7tM9OW+raB59+wWOMPWoMWoQAABiHGHvUGLQIAQDQIqVWnRJad8YeQQgAgBYpteqUzF3QnOt+0o1WH11jAAA00LzuhZq74KIxGUw9WnSj1UcQAgCggZodOkpBa3ZXZ1POnxq6xgAAGEequ9PWx+yuzqEQRRdZhiAEAEAiSuGnWWORxiO6xgAA2ACl1pVGjgVqxjkxPFqEAADYAM1oXWn0OYuzxVAbQQgAgAlqtOOJUp5eT9cYAACJS3l6PUEIAAAkiyAEAACSRRACAGAMzeteyKywNsJgaQAARrAhg4mLFy8sKp1nfa8OzQyw5iAIAQAwgtJg4vWZ1j5cYNqQmVkbc0XpeqrD1eyuTs3rXpjUzDG6xgAA2Ajj8SKIpTpL0pLuw4aCz6IF85ObOUaLEAAAo1SrS6sUIlq9GOr6dN+l1OIzEoIQAACj1OjurtGa2jF5xIBV6jprdSAbb+gaAwCgzb1mh65RB61FC+YPjSUiDI2MFiEAACYgur9GhxYhAADazIZMr8eGoUUIAIA2U9Gac8nnW1eRBNAiBAAAkkUQAgBgHKMbbePQNQYAwDjGoOiNQ4sQAABIFkEIAIAxQBdWe6JrDACAMTDeurDWZ8mO8YwWIQAAMKR6QdaJvggrLUIAAGDIRG79qYUWIQAAkCyCEAAASBZBCAAAJIsgBAAAkkUQAgAAdZVmkc3rXtjqqjQFs8YAAEBdpVlkpSn1Ew0tQgAAIFm0CAEAgI0ynq9CTYsQAADYKD29feP2KtQEIQAAkCyCEAAAGNFEnT3GGCEAADCiiTp7jBYhAACQLIIQAABIFkEIAAAkizFCAABgvRWvHTSe0SIEAADW23i+dlARQQgAACSLIAQAAJJFEAIAAMkiCAEAgGQRhAAAwKiVltoY77PFSpg+DwAARq201MZEQYsQAABIFkEIAAAkiyAEAACSRRACAADJYrA0AABoG8U1zMZiYDYtQgAAYIPN7ups6FT6sV7DjBYhAACwwWq12ox1q87GoEUIAAA0VKlVR5LmLrhI87oXtrhG9dEiBAAAmqLUGjR3wUUtrkl9tAgBAICGGI/Lb9AiBAAAGqLdxwPVQhACAABjrl0GVNM1BgAAxlRp8PRYTpOvhxYhAAAwpoqzylqNFiEAAJAsghAAAGiq0myydryeEF1jAACgqdr5ekK0CAEAgGQRhAAAQLIIQgAAIFkEIQAAkCyCEAAASBZBCAAAJIsgBAAAxkQ7rk7PdYQAAMCYaMfV6WkRAgAALdPqq07TIgQAAFqm1VedpkUIAAC0XKvGD9EiBAAAWq5V44doEQIAAMkiCAEAgGQRhAAAQLIIQgAAIFkEIQAAkCyCEAAASBZBCAAAJIsgBAAAkkUQAgAAySIIAQCAZBGEAABAsghCAAAgWQQhAACQLIIQAABIFkEIAAAkiyAEAACSRRACAADJIggBAIBkEYQAAECyCEIAACBZBCEAAJAsghAAAEgWQQgAACSLIAQAAJJFEAIAAMkiCAEAgGQRhAAAQLIIQgAAIFkEIQAAkCyCEAAASBZBCAAAJIsgBAAAkkUQAgAAySIIAQCAZBGEAABAsghCAAAgWQQhAACQLIIQAABIFkEIAAAkiyAEAACSRRACAADJIggBAIBkEYQAAECyCEIAACBZBCEAAJAsghAAAEgWQQgAACSLIAQAAJJFEAIAAMkiCAEAgGQRhAAAQLIIQgAAIFkEIQAAkCyCEAAASBZBCAAAJIsgBAAAkkUQAgAAySIIAQCAZBGEAABAsghCAAAgWQQhAACQLEdEq+sAAADQErQIAQCAZBGEAABAsghCAAAgWQQhAACQrI5WVwAAMDHYvlPS6lbXo0m2kvRUqyvRJBP1tU2PiN1HKkQQAgA0yuqI2KvVlWgG2zfy2sYX2zeOphxdYwAAIFkEIQAAkCyCEACgUb7f6go0Ea9t/BnV6+LK0gAAIFm0CAEAgGQRhAAAQLIIQgCAjWL7NNvL8usITRi2t7d9he27bd9l++RW16lRbE+3fb3t2/LX9rlW16nRbE+2fYvtC4crRxACAGys0yUd0upKNEG/pH+IiN0k7SPpI7Z3a3GdGmWNpPkRsYekPSUdYnuf1lap4U6WdM9IhQhCAICNEhFXS1re6no0WkQsjYib8+3nlP1Rnd3aWjVGZFbmd6fktwkze8r2HEmHSfrhSGUJQgAAjMD2XEmvlnRdi6vSMHnX0a2Slkm6PCImzGuT9A1Jp0oaHKkgQQgAgGHY3lTSuZJOiYgVra5Po0TEQETsKWmOpL1tj7gu13hg+3BJyyLiptGUJwgBAFCH7SnKQtCZEXFeq+vTDBHRK+kKTZxxXvMkHWF7iaSzJc23/bN6hQlCAADUYNuSfiTpnoj4Wqvr00i2t7bdlW93SnqLpHtbWqkGiYhPRcSciJgr6XhJCyPi3fXKE4QAABvF9lmS/ihpV9uP2n5fq+vUIPMkvUdZi8Kt+e3QVleqQbaTdIXt2yXdoGyM0LDTzCcqltgAAADJokUIAAAkiyAEAACSRRACAADJIggBAIBkEYQAAECyCEIAACBZBCEAAIZhe+XIpVrP9l62v5lvH2D7Da2u03jQ0eoKAACASrY7IqJ/fY6JiBsl3ZjfPUDSSkl/aHDVJhxahAAAGAVnvmr7Ttt32D4u3z/J9ndt32v7ctu/tX3sMOc5NC97k+1v2r4w3/9Z2z+1vUjST23Ptf172zfntzfk5c62fVjhfKfbPjZvBbrQ9lxJH5T08fxq2G+y/VC+bppszyzeTx0tQgAAjM4xkvaUtIekrSTdYPtqZUtxzJW0m6RtJN0j6bRaJ7A9XdL3JO0XEQ/ly5MU7SbpjRHRZ3sTSW+JiNW2d5F0lqS9JP1C0l9Jusj2VElvlvQhSa+XpIhYYvt/JK2MiP/In/dKSYdJ+rWy9bfOi4h1G/uGTAS0CAEAMDpvlHRWRAxExBOSrpL0unz/ORExGBGPK1vJvZ6XSXowIh7K71cHoQsioi/fniLpB7bvkHSOspAkSRdL+kvb0yS9VdLVhWPq+aGk9+bb75X04xHKJ4MWIQAA2sfzhe2PS3pCWQvUJEmrJSlvIbpS0sGSjpN09kgnjYhFeVfbAZImR8Sdja32+EWLEAAAo/N7ScfZnmx7a0n7Sbpe0iJJb8/HCm2rbKByPfdJ2ikfxyNlQaaeWZKWRsSgpPdImlx47BfKWnbeJOmSGsc+J2mzqn0/kfRz0RpUgSAEAMDonC/pdkm3SVoo6dS8K+xcSY9KulvSzyTdLOnZWifIu7A+LOkS2zcpCyw1y0r6rqQTbN+mrEut2Fp0maT9Jf1fRKytcexvJB1dGiyd7ztT0uZ6YXdc0hwRra4DAADjmu1NI2Kl7S2VtRLNy0PScGUt6TuS/hQRXx+DOh4r6ciIeE+zn2s8YYwQAAAb70LbXZKmSvpCvRCUe7/tE/KytyibRdZUtr+lbGD1oc1+rvGGFiEAAJrA9vmSdqza/cmIuLQV9UFtBCEAAJAsBksDALAeqtces32i7W8X7p9i+2/y7a/mV5G+3fb5efdZrXMeYvs+24ttLxhFHT6YX936VtvX2N6tRpld88dLtxW2T8kf+w/b89fvlU9MtAgBALAebK+MiE0L90+UtFdEfNR2h7JZY6+JiH7bB0lamG9/RZIi4pNV55ss6X5Jb1E2++wGSe+MiLuHqcPMiFiRbx8h6cMRccgw5SdL6pH0+oj4s+2/kPSDiDhoQ96DiYQWIQAAGme+pJtLC6ZGxGWFxVOvlTSnxjF7S1ocEQ/mU+HPlnTkcE9SCkG5GZJGatV4s6QHIuLP+fF/lrSl7ReN9IImOoIQAACNM0/STXUeO0nZ8hjVZkt6pHD/0XzfsGx/xPYDkv5d0t+PUPx4vfD6QTfn9U0aQQgAgI1XapHZTtKT1Q/a/rSkfmUXNWzME0Z8JyJeIumTkv6lXrl8YdYjlK1XVrRM0osbVZ/xiiAEAMD66cvDRckWkp4qPSZperFwPobocEnvitoDc3skbV+4PyffN1pnSzpqmMffqqy77omq/dPz+iaNIAQAwPq5StK7Jcl2p6S/UnnF+Xsk7VwqaPsQSadKOiIiVtU53w2SdrG9Yx6wjpd0QX78l20fXX2A7V0Kdw+T9Kdh6vtO1V5W46WSkl98lSAEAMD6OVnSMbZvVTYA+pyIuDp/7GJli7GWfFvZ4qeX51PY/0eSbL/Y9m8lKR9M/VFJlyoLUr+MiLvy418pqdZVqj9q+668Dp+QdEL1efP7M5TNRjuveLDtKcoC240b9A5MIEyfBwCggfIrSp8aEcO10oz2XJdGxMENqFb1eY9WNsX//zX63OMNQQgAgAayvaukbQutRG3H9jskXR4Rva2uS6sRhAAAQLIYIwQAAJJFEAIAAMkiCAEAgGQRhAAAQLIIQgAAIFn/H/GFxNlbvXCZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x720 with 9 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = utils.pairplot(samples[:,1:], points, limits=[[0,3],[1800, 2200],[1, 4]], fig_size=(10,10), \\\n",
    "                           labels = ['log_gamma \\n normal(0.,2.)', 'Tequ \\n U(0.,4000)', 'log_gravity \\n U(2.0, 3.7)'],\\\n",
    "                           title= 'SNRE \\n Guillot direct \\n' + str(100) + 'k simulations, '+ str(10) + 'k samples',\\\n",
    "                          points_colors = ['red'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# temp_params = {}\n",
    "# temp_params['log_delta'] = -5.5\n",
    "# temp_params['log_gamma'] = np.log10(0.4)\n",
    "# temp_params['t_int'] = 750.\n",
    "# temp_params['t_equ'] = 0.\n",
    "# temp_params['log_p_trans'] = -3.\n",
    "# temp_params['alpha'] = -0.25\n",
    "# pressures, temperature = nc.make_press_temp(temp_params)\n",
    "\n",
    "# # Create the Ratrans object here\n",
    "# rt_object = Radtrans(line_species=['H2O', 'CO_all_iso', \\\n",
    "#                                          'CO2', 'CH4', \\\n",
    "#                                           'Na', 'K'], \\\n",
    "#                     rayleigh_species=['H2','He'], \\\n",
    "#                     continuum_opacities = ['H2-H2','H2-He'], \\\n",
    "#                     mode='c-k', \\\n",
    "#                     wlen_bords_micron=WLENGTH)\n",
    "\n",
    "# # Create the RT arrays of appropriate lengths\n",
    "# rt_object.setup_opa_structure(pressures)\n",
    "\n",
    "# def Simulator(params): \n",
    "#     ##################\n",
    "\n",
    "#     NaN_spectra = 0\n",
    "    \n",
    "#     log_gamma = params[0].numpy()                            #log(0.4)\n",
    "#     t_equ = params[1].numpy()                                #1500.\n",
    "#     log_gravity = params[2].numpy()                                #params[5] 2.45 \n",
    "    \n",
    "#     # Make dictionary for modified Guillot parameters\n",
    "#     temp_params['log_gamma'] = log_gamma\n",
    "#     temp_params['t_equ'] = t_equ\n",
    "    \n",
    "#     log_P0 = -2.                                      #params[6] \n",
    "  \n",
    "#     gravity = np.exp(log_gravity)\n",
    "#     P0 = np.exp(log_P0)\n",
    "#     gamma = np.exp(log_gamma)\n",
    "\n",
    "    \n",
    "#     # Make dictionary for log 'metal' abundances    \n",
    "#     abundances_metals = {}\n",
    "        \n",
    "#     abundances_metals['H2O'] = -3              #np.exp(params[7]) * np.ones_like(temperature)\n",
    "#     abundances_metals['CO_all_iso'] = -2       #np.exp(params[8]) * np.ones_like(temperature)\n",
    "#     abundances_metals['CO2'] = -5              #np.exp(params[9]) * np.ones_like(temperature)\n",
    "#     abundances_metals['CH4'] = -6              #np.exp(params[10]) * np.ones_like(temperature)\n",
    "#     abundances_metals['Na'] = -5               #np.exp(params[11]) * np.ones_like(temperature)\n",
    "#     abundances_metals['K'] = -6                #np.exp(params[12]) * np.ones_like(temperature)\n",
    "    \n",
    "#     ####################################################################\n",
    "#     ####################################################################\n",
    "\n",
    "    \n",
    "#     wlen, flux_nu = \\\n",
    "#             rm.retrieval_model_plain(rt_object, temp_params, log_gravity, \\\n",
    "#                                          log_P0, R_pl, abundances_metals)\n",
    "    \n",
    "#     # Just to make sure that a long chain does not die\n",
    "#     # unexpectedly:\n",
    "#     # Return -inf if forward model returns NaN values\n",
    "#     if np.sum(np.isnan(flux_nu)) > 0:\n",
    "#         print(\"NaN spectrum encountered\")\n",
    "#         NaN_spectra += 1\n",
    "#         return -np.inf    \n",
    "\n",
    "#     # Convert to observation for emission case\n",
    "#     flux_star = fstar(wlen)\n",
    "#     flux_sq   = flux_nu/flux_star*(R_pl/R_star)**2 \n",
    "    \n",
    "#     # Rebin model to observation\n",
    "#     flux_rebinned = rgw.rebin_give_width(wlen, flux_sq, \\\n",
    "#                     data_wlen['MIRI LRS'], data_wlen_bins['MIRI LRS'])\n",
    "    \n",
    "#     FR= torch.Tensor(flux_rebinned)\n",
    "        \n",
    "#     return FR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:petitRT]",
   "language": "python",
   "name": "conda-env-petitRT-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}