GeostatsPy_bootstrap.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "\n",
    "<p align=\"center\">\n",
    "    <img src=\"https://github.com/GeostatsGuy/GeostatsPy/blob/master/TCG_color_logo.png?raw=true\" width=\"220\" height=\"240\" />\n",
    "\n",
    "</p>\n",
    "\n",
    "## Subsurface Data Analytics \n",
    "\n",
    "### Bootstrap for Subsurface Data Analytics in Python \n",
    "\n",
    "#### Michael Pyrcz, Associate Professor, University of Texas at Austin \n",
    "\n",
    "##### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Bootstrap\n",
    "\n",
    "Uncertainty in the sample statistics\n",
    "* one source of uncertainty is the paucity of data.\n",
    "* do 200 or even less wells provide a precise (and accurate estimate) of the mean? standard deviation? skew? P13?\n",
    "\n",
    "Would it be useful to know the uncertainty in these statistics due to limited sampling?\n",
    "* what is the impact of uncertainty in the mean porosity e.g. 20%+/-2%?\n",
    "\n",
    "**Bootstrap** is a method to assess the uncertainty in a sample statistic by repeated random sampling with replacement.\n",
    "\n",
    "Assumptions\n",
    "* sufficient, representative sampling, identical, idependent samples\n",
    "\n",
    "Limitations\n",
    "1. assumes the samples are representative \n",
    "2. assumes stationarity\n",
    "3. only accounts for uncertainty due to too few samples, e.g. no uncertainty due to changes away from data\n",
    "4. does not account for boundary of area of interest \n",
    "5. assumes the samples are independent\n",
    "6. does not account for other local information sources\n",
    "\n",
    "The Bootstrap Approach (Efron, 1982)\n",
    "\n",
    "Statistical resampling procedure to calculate uncertainty in a calculated statistic from the data itself.\n",
    "* Does this work?  Prove it to yourself, for uncertainty in the mean solution is standard error: \n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma^2_\\overline{x} = \\frac{\\sigma^2_s}{n}\n",
    "\\end{equation}\n",
    "\n",
    "Extremely powerful - could calculate uncertainty in any statistic!  e.g. P13, skew etc.\n",
    "* Would not be possible access general uncertainty in any statistic without bootstrap.\n",
    "* Advanced forms account for spatial information and sampling strategy (game theory and Journel’s spatial bootstrap (1993).\n",
    "\n",
    "Steps: \n",
    "\n",
    "1. assemble a sample set, must be representative, reasonable to assume independence between samples\n",
    "\n",
    "2. optional: build a cumulative distribution function (CDF)\n",
    "    * may account for declustering weights, tail extrapolation\n",
    "    * could use analogous data to support\n",
    "\n",
    "3. For $\\ell = 1, \\ldots, L$ realizations, do the following:\n",
    "\n",
    "    * For $i = \\alpha, \\ldots, n$ data, do the following:\n",
    "\n",
    "        * Draw a random sample with replacement from the sample set or Monte Carlo simulate from the CDF (if available). \n",
    "\n",
    "6. Calculate a realization of the sammary statistic of interest from the $n$ samples, e.g. $m^\\ell$, $\\sigma^2_{\\ell}$. Return to 3 for another realization.\n",
    "\n",
    "7. Compile and summarize the $L$ realizations of the statistic of interest.\n",
    "\n",
    "This is a very powerful method.  Let's try it out.\n",
    "\n",
    "#### Objective \n",
    "\n",
    "I want to provide hands-on experience with building subsurface modeling workflows. Python provides an excellent vehicle to accomplish this. I have coded a package called GeostatsPy with GSLIB: Geostatistical Library (Deutsch and Journel, 1998) functionality that provides basic building blocks for building subsurface modeling workflows. \n",
    "\n",
    "The objective is to remove the hurdles of subsurface modeling workflow construction by providing building blocks and sufficient examples. This is not a coding class per se, but we need the ability to 'script' workflows working with numerical methods.    \n",
    "\n",
    "#### Getting Started\n",
    "\n",
    "Here's the steps to get setup in Python with the GeostatsPy package:\n",
    "\n",
    "1. Install Anaconda 3 on your machine (https://www.anaconda.com/download/). \n",
    "2. From Anaconda Navigator (within Anaconda3 group), go to the environment tab, click on base (root) green arrow and open a terminal. \n",
    "3. In the terminal type: pip install geostatspy. \n",
    "4. Open Jupyter and in the top block get started by copy and pasting the code block below from this Jupyter Notebook to start using the geostatspy functionality. \n",
    "\n",
    "You will need to copy the data file to your working directory.  They are available here:\n",
    "\n",
    "* Tabular data - sample_data_biased.csv at https://git.io/fh0CW\n",
    "\n",
    "There are exampled below with these functions. You can go here to see a list of the available functions, https://git.io/fh4eX, other example workflows and source code. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geostatspy.GSLIB as GSLIB          # GSLIB utilies, visualization and wrapper\n",
    "import geostatspy.geostats as geostats    # GSLIB methods convert to Python        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also need some standard packages. These should have been installed with Anaconda 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np                        # ndarrys for gridded data\n",
    "import pandas as pd                       # DataFrames for tabular data\n",
    "import os                                 # set working directory, run executables\n",
    "import matplotlib.pyplot as plt           # for plotting\n",
    "from scipy import stats                   # summary statistics\n",
    "import math                               # trig etc.\n",
    "import scipy.signal as signal             # kernel for moving window calculation\n",
    "import random"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set the working directory\n",
    "\n",
    "I always like to do this so I don't lose files and to simplify subsequent read and writes (avoid including the full address each time). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#os.chdir(\"c:/PGE383\")                     # set the working directory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Loading Tabular Data\n",
    "\n",
    "Here's the command to load our comma delimited data file in to a Pandas' DataFrame object.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(r'https://raw.githubusercontent.com/GeostatsGuy/GeoDataSets/master/sample_data_biased.csv') # load our data form my GitHub account\n",
    "#df = pd.read_csv('sample_data_biased.csv') # load our data table locally"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's drop some samples so that we increase the variations in bootstrap samples for our demonstration below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using 58 number of samples\n"
     ]
    }
   ],
   "source": [
    "df = df.sample(frac = 0.2)                  # extract 50 random samples to reduce the size of the dataset   \n",
    "print('Using ' + str(len(df)) + ' number of samples')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualizing the DataFrame would be useful and we already learned about these methods in this demo (https://git.io/fNgRW). \n",
    "\n",
    "We can preview the DataFrame by printing a slice or by utilizing the 'head' DataFrame member function (with a nice and clean format, see below). With the slice we could look at any subset of the data table and with the head command, add parameter 'n=13' to see the first 13 rows of the dataset.  "
   ]
  },
  {
   "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>X</th>\n",
       "      <th>Y</th>\n",
       "      <th>Facies</th>\n",
       "      <th>Porosity</th>\n",
       "      <th>Perm</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>74</th>\n",
       "      <td>590</td>\n",
       "      <td>579</td>\n",
       "      <td>0</td>\n",
       "      <td>0.095128</td>\n",
       "      <td>1.650407</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>228</th>\n",
       "      <td>40</td>\n",
       "      <td>229</td>\n",
       "      <td>0</td>\n",
       "      <td>0.105188</td>\n",
       "      <td>3.080116</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>261</th>\n",
       "      <td>330</td>\n",
       "      <td>989</td>\n",
       "      <td>1</td>\n",
       "      <td>0.180307</td>\n",
       "      <td>207.432804</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>153</th>\n",
       "      <td>360</td>\n",
       "      <td>279</td>\n",
       "      <td>1</td>\n",
       "      <td>0.134607</td>\n",
       "      <td>9.776276</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>112</th>\n",
       "      <td>500</td>\n",
       "      <td>259</td>\n",
       "      <td>1</td>\n",
       "      <td>0.116832</td>\n",
       "      <td>3.513810</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>224</th>\n",
       "      <td>120</td>\n",
       "      <td>219</td>\n",
       "      <td>0</td>\n",
       "      <td>0.079483</td>\n",
       "      <td>0.870833</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>400</td>\n",
       "      <td>500</td>\n",
       "      <td>1</td>\n",
       "      <td>0.101716</td>\n",
       "      <td>7.029332</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>206</th>\n",
       "      <td>400</td>\n",
       "      <td>99</td>\n",
       "      <td>1</td>\n",
       "      <td>0.123665</td>\n",
       "      <td>4.412897</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>251</th>\n",
       "      <td>600</td>\n",
       "      <td>319</td>\n",
       "      <td>1</td>\n",
       "      <td>0.116017</td>\n",
       "      <td>9.836013</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>234</th>\n",
       "      <td>490</td>\n",
       "      <td>289</td>\n",
       "      <td>1</td>\n",
       "      <td>0.144525</td>\n",
       "      <td>16.251107</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>163</th>\n",
       "      <td>450</td>\n",
       "      <td>279</td>\n",
       "      <td>1</td>\n",
       "      <td>0.155484</td>\n",
       "      <td>26.561393</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>199</th>\n",
       "      <td>920</td>\n",
       "      <td>329</td>\n",
       "      <td>1</td>\n",
       "      <td>0.104723</td>\n",
       "      <td>6.333944</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>117</th>\n",
       "      <td>130</td>\n",
       "      <td>769</td>\n",
       "      <td>1</td>\n",
       "      <td>0.128778</td>\n",
       "      <td>16.159792</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       X    Y  Facies  Porosity        Perm\n",
       "74   590  579       0  0.095128    1.650407\n",
       "228   40  229       0  0.105188    3.080116\n",
       "261  330  989       1  0.180307  207.432804\n",
       "153  360  279       1  0.134607    9.776276\n",
       "112  500  259       1  0.116832    3.513810\n",
       "224  120  219       0  0.079483    0.870833\n",
       "18   400  500       1  0.101716    7.029332\n",
       "206  400   99       1  0.123665    4.412897\n",
       "251  600  319       1  0.116017    9.836013\n",
       "234  490  289       1  0.144525   16.251107\n",
       "163  450  279       1  0.155484   26.561393\n",
       "199  920  329       1  0.104723    6.333944\n",
       "117  130  769       1  0.128778   16.159792"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#print(df.iloc[0:5,:])                   # display first 4 samples in the table as a preview\n",
    "df.head(n=13)                           # we could also use this command for a table preview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Summary Statistics for Tabular Data\n",
    "\n",
    "The table includes X and Y coordinates (meters), Facies 1 and 0 (1 is sandstone and 0 interbedded sand and mudstone), Porosity (fraction), and permeability as Perm (mDarcy). \n",
    "\n",
    "There are a lot of efficient methods to calculate summary statistics from tabular data in DataFrames. The describe command provides count, mean, minimum, maximum, and quartiles all in a nice data table. We use transpose just to flip the table so that features are on the rows and the statistics are on the columns."
   ]
  },
  {
   "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>count</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>min</th>\n",
       "      <th>25%</th>\n",
       "      <th>50%</th>\n",
       "      <th>75%</th>\n",
       "      <th>max</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>X</th>\n",
       "      <td>58.0</td>\n",
       "      <td>437.586207</td>\n",
       "      <td>244.255477</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>280.000000</td>\n",
       "      <td>400.000000</td>\n",
       "      <td>607.500000</td>\n",
       "      <td>990.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Y</th>\n",
       "      <td>58.0</td>\n",
       "      <td>551.793103</td>\n",
       "      <td>290.084063</td>\n",
       "      <td>29.000000</td>\n",
       "      <td>300.000000</td>\n",
       "      <td>559.000000</td>\n",
       "      <td>834.000000</td>\n",
       "      <td>999.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Facies</th>\n",
       "      <td>58.0</td>\n",
       "      <td>0.793103</td>\n",
       "      <td>0.408619</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Porosity</th>\n",
       "      <td>58.0</td>\n",
       "      <td>0.133789</td>\n",
       "      <td>0.039522</td>\n",
       "      <td>0.065626</td>\n",
       "      <td>0.104839</td>\n",
       "      <td>0.128789</td>\n",
       "      <td>0.145473</td>\n",
       "      <td>0.221402</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Perm</th>\n",
       "      <td>58.0</td>\n",
       "      <td>255.523637</td>\n",
       "      <td>634.498889</td>\n",
       "      <td>0.239697</td>\n",
       "      <td>4.000799</td>\n",
       "      <td>9.806145</td>\n",
       "      <td>37.327095</td>\n",
       "      <td>3178.630458</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          count        mean         std        min         25%         50%  \\\n",
       "X          58.0  437.586207  244.255477   0.000000  280.000000  400.000000   \n",
       "Y          58.0  551.793103  290.084063  29.000000  300.000000  559.000000   \n",
       "Facies     58.0    0.793103    0.408619   0.000000    1.000000    1.000000   \n",
       "Porosity   58.0    0.133789    0.039522   0.065626    0.104839    0.128789   \n",
       "Perm       58.0  255.523637  634.498889   0.239697    4.000799    9.806145   \n",
       "\n",
       "                 75%          max  \n",
       "X         607.500000   990.000000  \n",
       "Y         834.000000   999.000000  \n",
       "Facies      1.000000     1.000000  \n",
       "Porosity    0.145473     0.221402  \n",
       "Perm       37.327095  3178.630458  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.describe().transpose()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualizing Tabular Data with Location Maps \n",
    "\n",
    "It is natural to set the x and y coordinate and feature ranges manually. e.g. do you want your color bar to go from 0.05887 to 0.24230 exactly? Also, let's pick a color map for display. I heard that plasma is known to be friendly to the color blind as the color and intensity vary together (hope I got that right, it was an interesting Twitter conversation started by Matt Hall from Agile if I recall correctly). We will assume a study area of 0 to 1,000m in x and y and omit any data outside this area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "xmin = 0.0; xmax = 1000.0               # range of x values\n",
    "ymin = 0.0; ymax = 1000.0               # range of y values\n",
    "pormin = 0.05; pormax = 0.25;           # range of porosity values\n",
    "nx = 100; ny = 100; csize = 10.0\n",
    "cmap = plt.cm.inferno                   # color map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try out locmap. This is a reimplementation of GSLIB's locmap program that uses matplotlib. I hope you find it simpler than matplotlib, if you want to get more advanced and build custom plots lock at the source. If you improve it, send me the new code. Any help is appreciated. To see the parameters, just type the command name:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<function geostatspy.GSLIB.locmap(df, xcol, ycol, vcol, xmin, xmax, ymin, ymax, vmin, vmax, title, xlabel, ylabel, vlabel, cmap, fig_name)>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "GSLIB.locmap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can populate the plotting parameters and visualize the porosity data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfkAAAF6CAYAAAAakiGSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3hUZfbA8e9JL6TQIfReRcQAiorY0Z+9gg3b6tpW195WxYJ1VVQsyCL23hCxYFkEFZcgCkpRei8hkADpmfP7497AJIQkDDO5w+R8nuc+ybxz33vPIHLmfe9bRFUxxhhjTOSJ8joAY4wxxoSGJXljjDEmQlmSN8YYYyKUJXljjDEmQlmSN8YYYyKUJXljjDEmQlmSN/WGiEwQkQfc34eIyCqvY9rXicgfIjLE6ziMMVWzJG/ClojcLiKTK5X9tZuyYUG+t4rIdhHZJiKbROQbETlnD+qH/EuEiCwTkQI3xvUi8rKINAjlPStT1V6q+l83nntF5PW6vL8xpnqW5E04+x44RESiAUSkBRAL9KtU1tk9N9j2V9UGQDdgAvCsiNwTgvvsjZPcGPsB/YG79vQCIhIT9KiMMWHBkrwJZzNxknpf9/Vg4DtgYaWyxaq6BkBEuovIFBHJEZGFInL23gahqtmq+hpwJXC7iDR273WxiMwXka0iskRErnDLk4HPgQy3lb1NRDJEZICI/CQiW0RkrYg8KyJxexufG+Nq95693RgyRGSi++ewSET+Vn6u2+J+X0ReF5E84CIRiReRp0RkjXs8JSLx7vlNRGSSG3eOiEwTkSj3vWUicrSIDAXuAM5xP+9vInKWiMzyj1NEbhSRj4PxmY0xNbMkb8KWqhYDP+Mkctyf04Dplcq+hx3JdQrwJtAMGA48JyK9ghTSJ0AMMMB9vQE4EUgFLgaeFJF+qrodOB5Yo6oN3GMNUAb8E2gCHAwcBVwVjMBEpA1wAjDbLXoLWAVkAGcCo0TkKL8qpwDvA+nAG8CdwEE4X572dz9jea/Aje61mgLNcZJ5hfWwVfULYBTwjvt59wcmAh1EpIffqecDrwXhIxtjasGSvAl3U9mZ0A/DSfLTKpVNdX8/EVimqi+raqmq/gJ8gJPk9pqqlgDZQCP39WequlgdU4Gv3Hh2V3+Wqs5wY1sGvAgcvpdhfSwiW3C++EzFSeZtgEOBW1W1UFV/BcYBF/jV+0lVP1ZVn6oWAOcB96nqBlXdCIz0O78EaAm0U9USVZ2mtdj0QlWLgHdwEjvul632wKS9/MzGmFqyJG/C3ffAoSLSEGiqqn8BPwKD3LLe7Hwe3w4Y6HYrb3GT33lAi2AEIiKxOK3ZHPf18SIyw+3C3oLTkm5STf2ubrf3OrebfNTuzheRF/y6+u+oJqxTVTVdVdup6lVuws4AclR1q995y4FWfq9XVrpOhnuO//kZ7u+PAYuAr9zHErdVE09lrwDniojgfGl4103+xpg6YEnehLufgDTgcuAHAFXNA9a4ZWtUdal77kpgqpv0yo8GqnplkGI5BSgF/uc+r/4AeBxorqrpwGRA3HOrauk+DywAuqhqKk63t1RxHqr6d7+u/lF7GOcaoJGIpPiVtQVW+9+iijrtKp2/xo1lq6reqKodgZOAGyp1/e/umqjqDKAYp4fjXKyr3pg6ZUnehDW3ZZoF3IDTTV9uulvmP6p+EtBVRC4QkVj36F/pmfAeE5FGInIeMAZ4RFU3AXFAPLARKBWR44Fj/aqtBxqLSJpfWQqQB2wTke44A/mCTlVX4vR2PCQiCSLSB7gU59n77rwF3CUiTUWkCXA38DqAiJwoIp3d1ngeztiCsiqusR5oXz4oz8+rwLNAqapO35vPZozZM5bkzb5gKs5AOv8EMc0t25Hk3e7pY4FhOK3QdcAjOMk4EL+JyDacrurLgH+q6t1+9/oH8C6wGaeVOtEvlgU4iXOJ++ggA7jJPW8r8BLO8+pQGY7z/HsN8BFwj6pOqeb8B3C+TM0B5gK/uGUAXYCvgW04PSvPlc+Nr+Q99+cmEfnFr/w1nMcq1oo3po5JLcbPGGNMwEQkEWcmQj93TIUxpo5YS94YE2pXAjMtwRtT90KW5EVkvIhsEJHf/coauQuV/OX+bOiWi4g87S7aMUdE+vnVGeGe/5eIjAhVvMaY4BORZcB1OHPtjfGciAwVZ6GsRVXNFBGRG0RknpuLvhGRdn7vlYnIr+4x0a+8g4j87Oapd8Rd5MpdZOod914/i0j7uviM/kLZkp8ADK1Udhvwjap2Ab5xX4OzcEgX97gcZxQyItIIuAcYiLM4xz3lXwyMMeFPVdu70/tm13y2MaElznLYY3ByTk9guIj0rHTabCBTVfvgLBj1qN97Bara1z1O9it/BHjSzW2bcQa64v7crKqdgSfd8+pUyJK8qn6PO5/Yzyk482Zxf57qV/6qu6jIDCBdRFoCxwFTVDVHVTfjrGZW+YuDMcYYUxsDgEWqusRdUfNtnPyzg6p+p6r57ssZQOvqLujOOjkS5wsB7JrbynPe+8BR7vl1pq6fyTdX1bUA7s9mbnkrKi7Oscot2125McYYs6f2NKdcirMnRLkEEclyF8EqT+SNgS2qWlrFNXfcz30/1z2/zoTL7lNVfbPRasp3vYDI5Thd/SQnJx/YvXv34EVn9sqWLVso2riRtPiKM9lyi4qIb9qU9PT0CuV5eXnkbNyAr6yU2PgEmjZvQUJCQl2GbIypwaxZs7JVtWmwr3vc0EzdlJ0bUN1Zsxb9ART6FY1V1bF+r/ckp5wPZFJx6em2qrpGRDoC34rIXJy1I3Z3zVrfL1TqOsmvF5GWqrrW7Y7f4JavAtr4ndcaZ37vKmBIpfL/VnVh9z/kWIDMzEzNysoKbuQmYMuXL+df51/A7b37UN5Tpao88sdc7nllAh06dNhx7meffMwXY0Zx/Rk9aJQUz6LsrTw+J4dHXn6H9u3be/QJjDGVicjyms/ac9nZOfw8c08XeXTERg0vVNXMak7ZXa6pQESOxtm06XD/ZZjLd7tU1SUi8l/gAJyVL9NFJMZtrftfs/x+q8TZ0jmNXR9jh1Rdd9dPBMpHyI/A2dWrvPxCd5T9QUCu253/JXCsiDR0B9wd65aZfUi7du3ofsQQXlowj5V5uazMy2Xcgnl0GXxYhQQP8Nrzo7llQGsaJTmt/s5NUvhb1xRef+l5L0I3xnhAfWUBHbUwE+jijoaPw1k4a6L/CSJyAM7mUSer6ga/8obit/0ycAgwz92s6Tt2boRVObeV57wzgW9rs7lTMIWsJS8ib+G0wpuIyCqcUfIPA++KyKXACuAs9/TJOJt7LALycbbtRFVzROR+nP8w4OySVaffgkztqSrffvMNH7/xDoUFhRx72omcevrpxMbGcvvIkXw9ZQqfv/c+oAy97VaOOfbYCvXLysqIKs6nQXzFHsAezVJ5a/68Ovwkpi7l5+fz5mtv8fUX39KocUMuvmIE/fv39zos4xVVdj7eDvaltVRErsFpLEYD41X1DxG5D8hS1Yk4GzI1AN5zex5XuCPpewAviogPp4H8sKqW/8N0K/C2iDyAMzr/P275f4DXRGQRTgt+WEg+WDUicsU76673xjOPP8GvH3zF0U27Ehcdw4/ZSyjq0pRn/jN2Rzd9TYYffySPHJhGSnzsjrL/rdjErBYDuOP+h0IVuvFIcXEx555xPqXL42nToDMFJdtZVPQbV9x6CcPPO8fr8Ew1RGRWDV3jATnwwPb6w893B1Q3MfbSkMS0L7MV70xQbNq0iW/e/5Rz2/enRYN0GiU24MQ2fSj+cw2zZs2q9XUuuvYGRs1Yxfqthagq89bn8p9F2zn/byHZy8V4bPJnkylaEU3Xhn1IjE2iUVJTDkwbwgtPjaW4uNjr8IwnooHUAA9TWbiMrjf7uAULFtA5Pp2oSi32rnENmf2/mWRm1u7L9XEnnEhCUjJPvvgsuZs20L5rDx4adzNt27YNRdjGYzOm/Y8msRkVyqKjYkgilZUrV9KpUyePIjNeUUqxp7LBY0neBEWzZs3YWFawS/nGsgJ6tGlTRY3dO3zIERw+5IhghWbCWNsObVj87Sya0qJCeYFvO40b1+l0YhNGQvVMvj6y7noTFF26dCG+XTNmbVhK+TiPZbkbWRi1jaOOOdrj6Ey4OmvYmWyMXUZe4WYAVH38lTuXAYP77bJ+gqkvFPWVBnSYXVmSN0Hz7xfGkNO7Of9eOpUnl0xlRsMCnnl1nC1kY3aradOmPPfK02xqvoSft37Jz/lf0u+UbjzwyH1eh1av5eXlMXnyZCZPnkxubmAL0wRMFdWygA6zK+uuN0GTkpLCqKcep6ysjLKyMuLi4rwOyewDevbsyXsT36a4uJiYmBiioqzt4aVvvv6WO28dSUxRI0AojXuEkaPu5Lihx9ZYNygkmiixXpxgsSRvgi46Opro6GivwzD7GPtS6L1t27Zx160jaRU/gNgk579HaVkJ99zxIAcdPJC0tLTQB6Gl+HzZob9PPWFfmY0xxgDw448/El3ckNjonV+4YqJjiSluzPTp0+ssDtXSgA6zK2vJG2OMAdjtolWqWusFrfaegq+kju4V+awlb4wxBoBBgwZRFr+ZktIde7JQUlZMWUIOhx56aJ3EoO6yttaSDw5ryRtjjAEgOTmZhx6/j9tvupvobQ0BKI3L4cFH7yU1tW5WlBOJIUpsjYRgsSRvjDFmhyFDDufrqZP58ccfUVUGDRpEgwYN6uz+qiX4ytbX2f0inSV5Y4wxFSQnJ3PMMcd4c3NVqN22saYWLMkbY4wJL/Z8PWgsyRtjjAkbgiK2RG3QWJI3xhgTRmKIimrmdRARw5K8McaY8KEl+ErXeB1FxLAkb4wxJowoWHd90FiSN/VGTk4OhYWFtGzZsg5X7zLG7CmxHeWCxpK8iXjZ2dmMvPlaCtYvIikmihxJ4fZRo+nVq5fXoRljKrMpdEFlSd5ENFXl9msvY0SbDfQ70FlFa21uIXf881Jeev+rOlvFyxhTWza6PpgsyZuI9ueff9KwaA392jbaUdYyLYGhbbbx1Refc+bZ53gYnTFmFxJLVExLr6OIGJbkTUTbsmULzZJ2LW+eJKzOtqUzjQk7WoKveJXXUUQMS/IhUFpaysSPJ/L5x5NJTE7inAvP4ZBDDvE6rHqpZ8+ePLnOR0mZj9jonZsufrdSueBvQ7wLzBhTNVVb8S6ILMkHmc/n45q/XUPe3M10S+1KcVkxD107ipMuPZkrrr7C6/DqnZSUFE6/5Dpuee1JzusVR1J8DJP+LCCp22Hst99+XodnjNmFIjbwLmgsyQdZVlYWG39fz+BmO1vuRyUdzvsT3mP4BcNtoJcHzj73Anrt349J779B4fZtDLnqNAYffniV0+gWLFjAq+OeYMWSP+nQpRcXXX4DnTp18iBqY+orG10fTJbkg+x/P/xMy6gWFcqiJIqmUU2YP38+AwcO9Ciy+q1Xr1706jWq2nPmzJnDY3ddyE0nxdLruCR+W/wzd19/Nvf8+y26du1aR5EaU7+JxBEV2zqE15ehwGggGhinqg9Xev8G4DKgFNgIXKKqy0WkL/A8kAqUAQ+q6jtunWlAinuJZsD/VPVUEUkDXgfa4uTbx1X15ZB9uCpYkg+ylq0zyPL9b5fybbqNpk2behCRqa1xzz7IPWfF06FlIgB9OzfgtlO3M27Mwzw6erzH0RlTP6ivGC1aEZJri0g0MAY4BlgFzBSRiao6z++02UCmquaLyJXAo8A5QD5woar+JSIZwCwR+VJVt6jqYX73+AD4xH15NTBPVU8SkabAQhF5Q1WLQ/IBq2BJPsiGnjCUsU+9SJvt2TRJboKqsnjLYhp3bkrHjh29Ds8AeXl5fPLRu/z5Rxat23ejQ+dezPzxK+bMmkazE1pVOLdnu2RWfbrQo0iNqY9C2l0/AFikqksARORt4BRgR5JX1e/8zp8BnO+W/+l3zhoR2QA0BbaUl4tICnAkcHH5qUCKOM8GGwA5OD0EdSaq5lPMnkhOTub5115gSeOVTN74FZOyPycuM4nRL4z2OjQDbNq0iStHnEzskjFcuN9vJC57ivtuPINeDabSPK2YRUuWsX7d2h3nr9xQSKMmGR5GbEz9IgriKwvoqIVWwEq/16vcst25FPh8lxhFBgBxwOJKb50GfKOqee7rZ4EewBpgLnCdqvpqE2iwWEs+BDp27MgbH75Bfn4+MTExxMXFeR2ScU0Y9zQXDsrj6H4NKSkpIapLEX3bJjNh+jquOLML47+Yz6VH5NCwUWO2FgoPfVjAJTff6HXYxtQje9WSbyIiWX6vx6rqWL/XVW1aoVVdSETOBzKBwyuVtwReA0ZUkbCHA+P8Xh8H/IrTuu8ETBGRaX5fAkLOknwIJSVVsQqL8UROTg5fT/mSiR+8w8nXNgQgPz+fBglCu+YxrFm/ncH9mlBW1o07Xl+IxG2geatOXHztAwwceJDH0RtTj0gcUXFtAq2draqZ1by/CvC/eGucVnbFEESOBu4EDlfVIr/yVOAz4C5VnVGpTmOcxwGn+RVfDDysqgosEpGlQHdg14FbIWJJ3kS86dOmMubRGzjlwDIuODiHB19dw2H7t+CswekUlEFBkSJRzpOrI/o3Y95K6HP8wxxxxBEeR25MPaTFaOGyUF19JtBFRDoAq4FhwLn+J4jIAcCLwFBV3eBXHgd8BLyqqu9Vce2zgEmqWuhXtgI4CpgmIs2BbsCSIH6eGlmSN57bunUrY5/9N1nTvyI6JpYjTjiLCy+5gtjY2L2+dnFxMc88egsvXJFManIMW/NiOWq/FYz6YD2D92+EFgoTvtvO0EPbArBo5XZmLErgSluh0BhvhHAXOlUtFZFrgC9xptCNV9U/ROQ+IEtVJwKP4QySe89dS2OFqp4MnA0MBhqLyEXuJS9S1V/d34cBFabjAfcDE0RkLs6jgltVNTskH243LMkbT5WVlXHd34ZzatvlXH1OKiWlJbw980VG3j6fBx4fs9fXnzNnDv3al5KanAxASmoKJSUtGNhlJTeMXUd8UlO2bC+jc3E8Py/aQHF0cx564mkbR2GMh8QXurFpqjoZmFyp7G6/34/eTb3Xcea87+66Q6ooWwMcG2iswWBJ3nhq+vTpdE9cyfF90gCIjhMuOiSN6z/4kZUrV9KmTcDP5gCIi4ujqKTiWJtGjRvTsCkce+aFXHb5VSQmJpKdnY3P56NZs2Z7dT9jzF6y/eSDypK88dSSP+exX7Nd/4fu3dzH0qVL9zrJ9+7dm0fWN2DF+kLaNk8AIL+wjI9mxvDEuAtITHQWvmnSpMle3ccYEyRR8UhC+wAr/xTMSCKCJXnjqY5de/Lzj9FU7h/7fX0UJ3bosNfXj4qK4t5HX+LOmy+jd0YeSfHKjEVRXPaPB63Vbkw48hWhBXU6Ni2iWZI3njr00EN55fk2fD5nOcf2TqWkVHl75lYadhq81634cl26dOH1D7/j119/pbCwkCsOOGBHC94YE240pM/k6xtL8sZT0dHRjH7pLcY++2/eeqd8dP0V3Hvp3/f4WuvWrWPypI/YmpvNwYcdR//+/XfsNBcVFUW/fv2CHb4xJhTUnskHiyV547mUlBRuvP1e4N6Ar/HDD9N47tF/cs6h0DUlmk9f+phPPzqUkQ8+SVSUrd5szD5DFawlHzSW5M0+z+fz8cyjt/PclamkNXDm1h/cG+57bTozZsxg0KBBHkdojKm1qHgkIdDxOL/WfEo9Y0ne7FPKysr4+eefyc7O5oADDqBNmzYsXryYTs1KSGuQXOHc/xsQx3+/mxR2SX7ZsmXMmTOH5s2b079/f+tpMMafrwjNr7zviwmUJXmzz1izZg23XD+CPu3yyWjk44G3fXTvdyJnD7+E3Pxd95jYsrWM5JR0DyKtms/nY9Td/2LpD1Pp1yCen0uVp2MSeWr8BJo2bep1eMaEB9Xa7ihnasGSvNlnPHjvDdx0FvTu3BKAc45TRo77jMWLD8eX2J6sBSvJ7J4KOHPh35hayv2jh3kZcgVfffEF23+ezt19u+0YEPjb+mxG3XE7T740robaxtQj9kw+aCzJm31CXl4ehbnL6N259Y4yEWHY0am8M/kd7nv4Oe66+QremrqSxqnC7yuFv107irZt23oYdUVffPAe57RruSPBA+zfvAmvzpxHcXGxLaVrDOBsNWtJPlgsyZt9hlS1E7SrSZMmvPDyB6xYsYKtW7dya5cu4Zc0VavczVoEnJ0ojTEoluSDyJK82SekpqYSl9KWeUu20LOj0yWvqrw9JZdjhp2947xwarlXdtwZZ/Hp6Ee5olfnHa35uRs2kdG1O/Hx8R5HZ0yYiIqHpE4BVv4rqKFEAkvyZp9x18inuPm6C+nbYS2tGvv4fo6PLvv/H4MHD/Y6tFo57vjjyfphOvfPmE6/5DjWlCpLohN4anzl3SlNKBQXFwOEXw+PqchXCNv+9DqKiGFJ3uwzMjIyePXtL5kxYwbZ2dncMfyAsG65VxYVFcW/Rj3EkiVLmDt3Lgc0a8bAgQNtCl2IbdiwgZG3/oul8/5EFbrs34N7Hr6fxo0bex2aqYIoiM8eXwWLJXmzT4mOjuaQQw7xOoy90rFjRzp27Oh1GPVCWVkZV15wKQNL23Jka2cbpD8Xr+LqEZfz5sT37AtWuLJn8kFjf8ONMRHrp59+otHWWDqlZ+wo69qwNYmby5g9e7aHkZndc0fXB3KYXVhL3hgTsdavX0+6JuxSnuZLZP369R5EZGoUlQANOgdYeWVQQ4kEluSNMRGrT58+vK0vcZBfmaqy3LeJ3r17exaXqUZZIeTZwLtgse56U6fKysooK7MlK8NVSUlJRM3Z79KlC10P2Z/Jq2eyqSCPjflbmLT6f/Q7ZtA+NWiz3lEN7DC78KQlLyL/BC7DWfZgLnAx0BJ4G2gE/AJcoKrFIhIPvAocCGwCzlHVZV7EbQKXk5PDY/ffxtI/ZgLQvkcmN9/9sI1wDhNZM2fy2L2jyN+US1kUnDTsdK645qqIGJh2/2MP8dmkSUx692MkKopzrr2SoccP9Toss1u24l0w1XmSF5FWwD+AnqpaICLvAsOAE4AnVfVtEXkBuBR43v25WVU7i8gw4BHgnLqO2wTO5/Nx49/P4+Ieazn4wgYAzFg0kxv/fh7j35kcEYlkX7ZkyRLuvfoWhrUYQHrbBpT6yvjq7W94rrSMa264zuvw9lpUVBQnnXwyJ518stehmNqwFe+Cyqt/XWOARBGJAZKAtcCRwPvu+68Ap7q/n+K+xn3/KJHqFjg14SYrK4tOCWs5uEuDHWUHdW5Al8R1zJw508PIDMDr4yYwJKUL6QnOf5+YqGiOa9WXye99TGlpqcfRmXonOgFSugV2mF3UeUteVVeLyOPACqAA+AqYBWxR1fJ/UVYBrdzfW+EOmVTVUhHJBRoD2f7XFZHLgcshvJc2rY/WrVtH+7TiXco7pBWxbt06DyIy/lYvW8Fhyc0rlEVLFEkSS35+PqmpqR5FZuqlskLIXeB1FBGjzlvyItIQp3XeAcgAkoHjqzi1fBRFVa32XUZYqOpYVc1U1Uzbmzu89O7dm59X7bo2+4zVCTbCOQwceMhAFm5ZU6Esv6SQ0vgoUlJSPIrK1F8K6gvsMLvworv+aGCpqm5U1RLgQ2AQkO523wO0Bsr/1VkFtAFw308Dcuo2ZLM3OnbsSJMeR/DElC2s3VLMui3FPPV1Lo26HU6nToFuRGGC5Zzzz+X3uE3MXPcX+SWFLM9dz5urfuKa22/CnoyZOqeATwM7zC68GF2/AjhIRJJwuuuPArKA74AzcUbYjwA+cc+f6L7+yX3/W42kOT71xL8eeJxJEz/hsYlvAnD0qcM46eTTPI7KAKSlpTH+vTd4bfwrfDrtR1p0yeD+K0az33777Thn8eLFzJs3j1atWnHAAQdY8jehZQk7aMSLfCkiI3FGyJcCs3Gm07Vi5xS62cD5qlokIgnAa8ABOC34Yaq6pLrrZ2ZmalZWVgg/gTH1Q1lZGTdedwu//jCfxOKGlMblk942nv+8Ntae1ddzIjJLVTODfd3MLg305yf6BFQ35uSfaoxJRIYCo4FoYJyqPlzp/RtwclIpsBG4RFWXu++NAO5yT31AVV9xy/+LMw28wH3vWFXd4L53NnAvTh/Fb6p6bkAfLkCezJNX1XuAeyoVLwEGVHFuIXBWXcRljKnonbffZd7UFeyXfviOstXLFvPAvQ/x6BMPeRiZiVRaWohuCc3AOxGJBsYAx+A8Cp4pIhNVdZ7fabOBTFXNF5ErgUeBc0SkEU7eysRJ2LPcupvdeuepaoXWpYh0AW4HDlHVzSLSLCQfrBo2QbkK+fn5zJs3j5ycunn0v3btWubPn79jv2tjwsX7b31E+wa9KpRlpHTkh6k/eRSRqRdC90x+ALBIVZeoajFO7/Ep/ieo6neqmu++nIEzRgzgOGCKqua4iX0KUNOqSn8DxpR/EShv3dclW7u+khefG8tb498mLaoRW8vy6D+4H/c/ch9xcXFBv1deXh63XHsj6xeuJCU6kQ2+PK65/XpOPPmkoN/LmECoT6lqgotUOenFmCBQIPCB8k1ExL81PVZVx/q93jEl27UKGFjN9S4FPq+mbiu/1y+LSBnwAU5XvgJdAUTkB5zHA/eq6hd78Hn2miV5P19//TUfvjiRwxsPRcTp5Pjjv7/x1OOjueWOm4N+vztvuJVmy4TDM5yu0OKyEl68fzRdunWlWzdb2MF47/Rhp/L6Y5/QNf3AHWVrti7loCG7PFkzJngCT/LZNTyTr9WUbAAROR+na778WVV1dc9z14BJwUnyF+Asxx4DdAGG4PQITBOR3qq6paYPEizWXe/n1bGv0SvlgB0JHqBHw/34/JPgf/HKzc1l6dy/6NG4w46yuOhYMpO68N4b7wT9fsYEYvi559D54Ob8sW0ai7LnsmDrDHytN3HXyNu9Ds1EqugESO8e2FGzHVOyXf7TtXcQkaOBO4GTVbWoprqqutr9uRV4k53jy1YBn6hqiaouBRbiJP06Yy15P3l5ebSJSaxQFiVR+EoVVQ3qtKH8/HwSo3ddIKZBXCLrNm2uooYxdS8mJoYxY59h/vz5LFiwgFatWpGZmcsId2sAACAASURBVGn7DZjQKSuEzSFb8W4m0EVEOgCrcfZNqTDaXUQOAF4EhlZ6hv4lMMpd0A3gWOB2d/2WdFXNFpFY4ETga/ecj4HhwAQRaYLTfV/t7LBgsyTv5+jjj2b6hJl0b7RzFbaN29fTuUfHoM8LbtGiBcUJSl7RdlLjk3eU/7F1GWefeHlQ72XM3urRowc9evTwOgxTH+zdM/nqL+0sjX4NTsKOBsar6h8ich+QpaoTgceABsB77r/7K1T1ZFXNEZH7cb4oANznliUDX7oJPhonwb/knvMlcKyIzAPKgJtVdVNoPl3VPJknH2qBzpPfvn07F549Al0bS9OY5uSVbmFjwhrGvTmW9u3bBz3OWVlZ3Hn1LfSJaUNaXAP+LFhDcs9mPPPSc0RHRwf9fsYYEyyhmid/YLsYnXFHYGswxP19c0hi2pdFZJJPS0nTy0ZcwnW3/HOPN6spLi7myy+/5JefZ9OhcwdOPf2UkC76sWHDBj56/0M2rd/IIUcO5rDDDrOuUGNM2Atpkr8twCR/lSX5yiIyybdPb60j+pzOz6XzmfDha7Ro0cLrkIwxJqKELMm3jdEZt6YFVDfumhxL8pVE7DP5dumtyN9UyGvjX+XmO27xOhxjjDG1EZNQ25HyVfgxqKFEgohN8gBtUlrw629/eB2GMcaY2iotRHMWeh1FxIjoJL962wa6DLJFZYwxZt8h4LMVFYMlYkd4rd66nl9L/2LEZRd5HYoxxpjaUkAlsMPsIiJb8ptL81jXoYDn73qJjIwMr8MxxhizB9Ra8kETkUm+S/eujBn/vNdhGGOM2VMxCUijQAfe1fkmb2EvIpO8McaYfVRpIZptA++CxZK8McaY8KHYwLsgsiRvjDEmjAhqg+iCxpK8McaY8OKL2Ilfdc6SvDHGmPARnQCNA13fZHVQQ4kEluSNMcaEDS0tRDf+6XUYEcOSvDHGmPBiz+SDxpK8McaYMCK2GE4QWZI3xhgTPhRQG3gXLJbkjTHGhI+YBKRJ1wArLwtmJBHBkrwxxpjwUVqIb8NfXkcRMSzJG2OMCSu2GE7wWJI3xhgTPlRsMZwgsiRvjDEmrNjo+uCxJG+MMSZ8xMQjzboEWNme5VdmSd4YY0wYse76YLIkb4wxJmxoSRG+9Yu8DiNiWJI3xhgTVuyZfPBYn4gxxpiwoioBHbUhIkNFZKGILBKR26p4f7CI/CIipSJyZqX3HhWRP0Rkvog8LSLiln8hIr+5770gItFu+f0iMkdEfhWRr0QkIwh/PHvEkrwxxpjwoTjP5AM5auAm3zHA8UBPYLiI9Kx02grgIuDNSnUHAYcAfYDeQH/gcPfts1V1f7e8KXCWW/6YqvZR1b7AJODuPf3j2FvWXW9MJfPnz2fp0qV07tyZrl0DXV5z31ZSUsKMGTMoKChgwIABpKenex2SqS9iE4hq3jnAyr/XdMIAYJGqLgEQkbeBU4B55Seo6jL3PV+lugokAHGAALHAerdOnntOjPu+VioHSC4vr0uW5I1xFRYWcv3lV7Nt0VpaSAqv61aa9GrH42NGExcX53V4dWbBggVceck/iNmeivhi2B7zENfc/HeGn3uO16GZ+qCkiLK1i0N19VbASr/Xq4CBtamoqj+JyHfAWpwk/6yqzi9/X0S+xPkS8Tnwvl/5g8CFQC5wxN5+gD1l3fXGuJ4b/SwNlxZwTuuDObxVb4a3PpjY+Tm8/NJ/vA6tzvh8Pq694gZal/Sjc8qBdErbn16JR/LsI2NZvny51+GZekAB1cAOoImIZPkdl1e6fFUP7mvVuhaRzkAPoDXOl4UjRWTwjouoHge0BOKBI/3K71TVNsAbwDW1/oMIEkvyxrj+O/krBjav2D0/qHk3vvxokkcR1b0FCxagW+NoEJ+2oyw6KoaGvnZ8+slnHkZm6g0F9UUFdADZqprpd4ytdPVVQBu/162BNbWM7DRghqpuU9VtOC32gyqErloITMR5BFDZm8AZtbxX0FiSN8alqiAVv+iLCD5f5Udzkcv5rLs2dqKIoqysrO4DMvVSCEfXzwS6iEgHEYkDhuEk5dpYARwuIjEiEosz6G6+iDQQkZYAIhIDnAAscF/7L913cnl5XbJn8sa4Djl6CL98s5D+zXf+f/m/dX9x9BnHexhV3erZsye+xHwKireRGNcAAJ/6yIlazokn3+5xdKY+kNgEYloEOvDu12rfVdVSEbkG+BKIBsar6h8ich+QpaoTRaQ/8BHQEDhJREaqai+c5+xHAnNxuvi/UNVPRaQ5MFFE4t1rfgu84N7yYRHpBviA5cDfA/xgARPVOh/sF3KZmZmalZXldRhmH7N9+3auHHEZsWu304IGrNY8ots3Ysz4sSQkJHgdXp357bffuOayf5JY0BTxxbAtdi0XX3Uel11xqdehmTAiIrNUNTPY1+3bLEm/Pqt7QHWbPjc7JDHty6wlb4wrOTmZV957k6ysLJYuWcKZXbvSt29fROrX6lv7778/X06dxPfff8+2bds49NBDadGihddhmXpDbMW7ILIkb4wfEaF///7079/f61A8lZSUxNChQ70Ow9RHSq1XrzM1syRvjDEmbDhT6GxMeLBYkjfGGBM2JDae6BadAqw9M6ixRAJL8sYYY8KGlhRRumaJ12FEDEvyxhhjwkjtd5SrT0QkD2cRCwWSgAJ2rtaXrKrRVdWzJG+MMSZ82MC7KqlqavnvIvKLqvbzf727epbkjTHGhBWbQlczEYlR1VL3ZezuzrMkb4wxJnzExhOTEejAu5+CGkoY+x54X0Q+Bw7Db6vcyizJG2OMCRtaUkTJ6qVehxHubgIuBfYDsti5jO4uLMkbY4wJHypg3fVVEpEDgME4A+6mq+qLNdWxFQeMMcaElRDuQrfPEpHrgZdxNs5pBLwsIjfUVM+TJC8i6SLyvogsEJH5InKwiDQSkSki8pf7s6F7rojI0yKySETmiEi/mq5vjDFm32VJvkqXAgep6r2qei8wELi4pkpeteRH42zT1x3YH5gP3AZ8o6pdgG/c1wDHA13c43Lg+boP1xhjTF1RnwR0RDjF2cq2XDQ758nvVp0/kxeRVJxnChcBqGoxUCwipwBD3NNeAf4L3AqcAryqzp64M9xegJaquraOQzfGGBNiEhtPTOsOAdaeGtRYwsxLODnwI/f16W5ZtbwYeNcR2IjzPGF/YBZwHdC8PHGr6loRaeae3wpY6Vd/lVtmSd4YYyKMFhdRsnKZ12GEHVV9RkS+x5kyB3Ceqv5WUz0vuutjgH7A86p6ALCdnV3zVamqD2aXLgoRuVxEskQka+PGjcGJ1BhjTJ1SBNWogI5IJiLtgC3Ap+6xxS2rlhct+VXAKlX92X39Pk6SX1/eDS8iLYENfue38avfGlhT+aKqOhYYC5CZmVnjcwpjjDHhyRf5g+gC8Sk7165PBtoBfwI9q6tU5199VHUdsFJEurlFR+Gs1jMRGOGWjQA+cX+fCFzojrI/CMi15/HGBIeq4vP5vA4jZMrKyrwOwQTARtfvSlX7qOp+7s9OwMHAdzXV82oxnGuBN0QkDliCMw0gCnhXRC4FVgBnuedOBk4AFgH51GLKgDGmeiUlJTz176eZ+MFnaCm0bt+Sex/6F927d/c6tKD4bNJkRj/2LPlbi4hPiuXqf17B6Wec5nVYphYkNp7Y1u29DiPsqepMEdntSnflPEnyqvorkFnFW0dVca4CV4c8KGPqkfvufoDZn/1Fv7RjiJIoclfncMUFV/PeZ2/RrFmzmi8QxqZNm8ZDtz9Ft+SDiUtKoKSsmCfveZHExESOP2Go1+GZGmhJEcUrl3sdRlgSkRY48+MVmAmcKiLi5skqRfZIBWPMLrZu3crUr36gc1ofosT5JyAtoRGNi9rz7lvveRzd3nvmiefpmHAgcTEJAMRGx9E5uT/PP13jbCMTBlQDPyKZiAzD2YHndPeYARxaXYIHW7vemHpn06ZNJEgyIhWfYabEpbNk0TJvggqijeuz6RbXu0JZQmwSeVvyPIrI7BmJ+JHyAboDOFBVcwBEpBHOejJvVVfJ/iSNqWcyMjIoitpGma+0Qvmm4rUcdGh/j6IKnp69u5OTv65C2ZaCbNp3qnG2kQkTPpWAjnog1+/3LbWpYC15Y+qZuLg4/n7d5Tz/8Hg6JexPYmwyq7cuIapVASedfJLX4e21G2+7ngvPvoySvGIaJ7ckJ38966PmM+6u57wOzdSCxMUTZwPvqvIZ8IWIlLfcz8MZmF4tS/KmRgUFBXzw3od8//U0mrdqzgUXn0fXrl29DiskNmzYwNuvvMZff8yne98+DLvgPBo3bux1WEE3/LxzaN+xLeNfeIV12cs4dvgRXHjR+SQmJtZ5LMXFxXz80Sd8M/kbGjdtzPmXnEfPntVO/a1Wx44deeujV3hxzDh+/20u3QZ15d/XjKdDh0CXSvXeihUrePPlV1mxZDmZhx7EWcPOJiUlxeuwQkKLiihaYQPvKlPV20XkJOBwt2i0qk6sqZ7U8Mx+n5SZmalZWVlehxERCgoKOP/sC2FlPK2S2rO1OI+lZfO45993csSRR3gdXlAtXbqU6y+4lMEJLeiY2pRFuRv4oXQjY958hVatWnkdXkQqKSnhkvMuoXhRKR2TO7K9ZDvziv7gunv/wUmn7Pu9CsEw+5dfuOvKmzg4sSMtkhvzV+4q/krMZfx7r5Oenu5ZXCIyS1WrmiW1V3qlp+m7hx8aUN3eEyfXGJOIDMXZJC0aGKeqD1d6fzDwFNAHGKaq7/u91xYYh7NAmwInqOoyEbkGuB7oBDRV1Wz3/O4428P2A+5U1ccD+mDOtfoAq1V1k4ik4ywR/6uqVrvQhT2TN9X6+MOP0ZXx9GjUl9SEdFqltqV/yhAeuueRiFtE5emHHuO0tM4MaNGRJkkpHNSyE8cntuH5J0Z7HVrEmvLVFAoXF9OvyYGkJzakVWprhjQ6ktEPjaakpMTr8MLCY/eO4vQm/enVpAONE1M5qEVPehQ04vWXX/U6tBAJbCGc2iyGIyLRwBic3U17AsNFpHK30QqcDdTerOISrwKPqWoPYAA7V2b9ATgaqNwFkQP8Awg4uft5GWcztySc6XMPA/+pqZIleVOtad/+QEZi2wpl8TEJSGE02dnZHkUVGkvmL6RTw4pzxHs0zuD3WbM9iijyTf92Oq3jWlcoi42OpYGmsGLFCo+iCh+lpaXkrc+hUWJqhfKejdsx47/TPYoq9EI48G4AsEhVl7g7oL6Ns9PpDqq6TFXnABVaMe6XgRhVneKet01V893fZ6vqsso3U9UNqjoTCMY31ihV3QocB0xW1WOBA2uqZM/kTbVaZDRn8ex1pCc22lGmqhRpIQ0aNPAwsuBLbJDE9pIikmPjd5TlFuWTmp7mYVSRrWXrlvxaMpcWtKxQvt23nYYNG3oUVfiIjo6mNEop85URHbVzK/FNBXk069jcw8hCR+LiiG8TspkQVe1qOrCWdbvibArzIdAB+Bq4TVXrau1kdXduvRinNwJqsZ+8teRNtc6/+DyW6QIKSvIBJ8H/uWUuhx19CElJSR5HF1znXDKCD1fOodTn/D9bUlbGB6t/57wrLvM4ssh11vCzWMpithVvA5y/Xws3L6BnZk8aNWpUQ+3IJyKcdPZpfLv2N3zuo9fC0mK+zfmdC/9+qcfRhYYWF1O4YkVAB9CkfDdS97i80uVrtavpbsTgbPN6E9Af55n4RQF+zEDchtNlvwH4SkRSgfE1VbKWvKlWx44defDpkTz4r4co26oU+QoZcsJg7rznDq9DC7pTzzyD3C25PDHhddKiYsnTUoZdeTHHHW9LoYZKixYtePSFRxh5232UbCmh2FdM/8H9ueeBu70OLWxced01jC4u4T8ffUaD6HgKYnxcNfJG+vbt63VoIeGsXhfwnPfsGgbe1WpX02rqzlbVJQAi8jFwELV4Lh4MqvoV8JVfUR7OAMJq2eh6UyuqyubNm0lKSiIhIcHrcEKqtLSU3Nxc0tLSiImx78F1ofzvV2JioifT+PYFxcXFbN26lYYNGxIV5X0nbKhG1/dMS9fXDxkSUN0DP/+k2phEJAZne9ajgNU4A9jOVdU/qjh3AjCpfHS9O2jvF+BoVd0oIi8DWao6xq/OMiCzfHS9X/m9wLa9HF0/nip6IlT1YhEZqar3VFXP/gUztSIi9ab7NCYmJiLnxoez+vT3K1BxcXH15O9l6LaNVdVSd7rblzhT6Mar6h8ich9Owp4oIv2Bj4CGwEluAu2lqmUichPwjThrQs8CXgIQkX8AtwAtgDkiMllVL3M3lMkCUgGfiFwP9FTVQNZYnlTNe1N394YleWOMMdXatm0b69atIyMjI+RjcSQujoS2bWs+MUCqOplKK8Wp6t1+v8/E6cavqu4UnPnzlcufBp6uonzd7q61p1T1w8plInKt+963u6tnSd4YY0yVfD4fjz/0GF9/8hXp0alsLsvl1PNO56p/hG73by0upmD5yppPrGfcHoi/4fQKlMsQkRuAp1S1yufzluSNMcZU6bUJr/LrRzM5tdlxiAg+9fHdq1No1SYojdMqKXs18C6SXY2ziE95V7/i7EJ3BJC/u0q1SvIikokzdSADKAB+B74u3/LOGGNM5Png9fc5qsnBO7YljpIoDm7cjzf+83robrp3o+sj2ZrKC+6ISHZNebjaIZoicpGI/ALcDiQCC3Hm6B0KTBGRV9y1fI0xxkSYosIiYqNiK5QlxCSwbeu2EN41sNXuInWrWRFJA1DVoyqVJwFv1FS/ppZ8MnCIqhbs5uZ9gS44a/0aY4yJIL369mbFb6tol75zavlfm5dy0JEH8cW0L0N2X2vJV/A/oFv5CxEZAFwKHAl8WlPlapO8//y/3bz/a+1iNMYYs6+56a6buXTYJWzOzqV5fFPWFq1nbXIOE/75APeOGhmSe0pcHInt2tR8Yv0xW0Q+BaYDZ+Msy/sycFVtltSt7TP5DsC1QHv/Oqp6cgABG2OM2QdkZGTwzqR3+ej9D/nzjz854oDjOeW0U0hOTg7ZPX3FxeQvWxWy6+9rVHWYiByBM7K+CTAFWFDbNfNrO7r+Y5yl+z6l0s48xhhjIldqaiojLrmoTu8Zqc/XA6Wq3wHfuc/nhwOvikgp8LKqjquubm2TfKE72d8YY4wJKa1yHxmjqrnAC8ALItILZ0e6atU2yY8WkXtwFscv8rvhL4EE6oWFCxfyzoTXyNmYzaHHHslJp5xCfHx8zRVNxPl+6vc88dCTbMnOIyE5niv+8TdOO+M0r8MyxgBo6Ja1jSTuevs31XRebZP8fsAFOKP5yrvr1X0d9r764kteuOchjk7vRJeEZGY++xaT3v2QsW++SlxcnNfhmTo0a9Ys7r7ufvqlHEpSajLFpUWMGfkSUdHRnHKqDTExxmsSH0ti+9AttlPf1DbJnwZ0VNXiUAYTCmVlZTw76jH+3u4gEmOchD40OZ3PVs3li8mfc/Kpp3gcoalLzz35Ar0SM0mKdQYOxcXE0zf9YF56ZpwleWPCgK+ohO1LV3sdRsSo7X6FvwHpoQwkVFavXk1j4nYk+HL7pbZkxne73bjHRKhVK1eTnlBxt7OEmES25W33KCJjTGXqdtnv6WF2VduWfHNggYjMpOIz+bBv+qSnp7OlpHCX8k2F22jeqqMHERkv7de3N2t/WEVG6s55uHlFubRo1dzDqIwx/mx0fUXulrWXAv+HszBODLAZ+AF4w91Zr0q1TfJVbka/L0hNTaXbgL78MPcvBjXvjIiQW5TP1G0rGHPuKK/DM3Xs2huu5qKfLsWXW0ZGShuy8zewsORX/n3XI16HZozB3aDGRtfvICJn4syRfxc4X1WXuOXJwIHAKSJyY+Vlb8tVm+RFRNSx237t8nMC/gR14J6HH2TUv0by1PRpNIiJpzgphjueepjWrW1wR33Trl07Xnl/PC888yJzZv9ApwM78eJ1Y+jevbvXoRljgKi4OJICHXg3LbixhIkvVPX9yoWquh34HvheRBrsrnJNLfnvROQD4BNV3bE+vYjE4WxSMwL4DpgQQOB1JjExkfsff5j8/Hzy8/Np3Ljxjl2VTP3Ttm1bRj32oNdhGGOq4CsuZputeLeDqu52NyAReV5Vr6zunJqS/FDgEuAtd2nbLUACEI0zZ/7JfWn9+qSkJJKSkrwOwxhjzG7ZILqqiMgSqPAcQ4FWIjIEeH53C9bVlOQ/BK5W1edEJBZn3dwCVd0ShJiNMcaYClRt4N1uZFZ6rcB/gZOBz4GAkvwE4EsRmQA8pqpr9ypEY4wxpgbWkt+VquZULhOR8aq6XERm7a5eTVvNvisinwF3A1ki8hp+G9So6hN7EbMxxhhTQVR8LMkdWgVW+cfgxhJO3KXlq6SqF+zuvdpMoSsBtgPxQAq2C50xxpgQ8RWVsNVWvKvKVr/fk3G66WvcP6amKXRDgSeAiUA/Vc3fmwiNMcaYmlh3/a4q95yLyCM4279Xq6aW/J3AWe5uN8YYY0xIKZbka0NVi0Vko4hEq2rZ7s6r6Zn8YcEPzRhjjNkdsRXvqiAinYGrgFzgSaAYuLO6BA+1X9bW1DNTpkzhpTEvs3nTZg4+bCDXXH8VzZo18zosY0yEi46PpUGgA+9mBDeWMPMBzoy3lsCzwGXAa8Dg6ipZkje7eP3VNxj76Ot0SupHk9jezP10GedOu5APJr1DWlqa1+EZYyJYWVEJeUvWhOz67liz0TiLuo1T1YcrvT8YeAroAwwrX1JWRNrhrB0TDcQCz6jqC+57w4E7cJ42rMFZYz5bRPYHXgAaAMuA81Q1L8DQt6vqk+795rjd9TWu7lbbrWZNPVFSUsLYZ8fTM/UQkuJSEIkiI60jCVua8/ab73gdnjEmwjmL4QR21EREooExwPFAT2C4iPSsdNoK4CLgzUrla4FBqtoXGAjcJiIZIhKD86XhCFXtA8wBrnHrjANuU9X9gI+Am/f8T2SH70TkYvczlLnd9zWyJG8q2LRpE7G+RKKjKnbyNIpvweyZczyKyhhTn6j7XH5Pj1oYACxS1SWqWgy8DZxS4d6qy1R1DpWmi6tqsaqWb7Uez878Ke6RLM6mKKk4rXlwtoX93v19CnDGHv1BVHQtzpeGAqAL8BY7v0zslnXXmwoaNWpEMfmo+hDZ+R1wS1E2B/Ts42Fkxpj6IoTL2rYCVvq9XoXTKq8VEWkDfAZ0Bm5W1TVu+ZXAXJw1Zf4Crnar/I4zn/0T4CygTaCBq2pqIPWsJW8qiIuL48xzT2Vh3kxKyooByMnfQG7ics67cLjH0Rlj6gPVwA6giYhk+R2XV7p0Vd8ear1VuqqudLvkOwMjRKS5u6/LlcABQAZOd/3tbpVLgKvdZWdTcEbEB0xEjhORf4vI4yJyfG3qWEve7OLa668hvWE6r7/8NkUFRXTu1ZGXR4610fXGmJCLio8lpWNGYJWzyFbVyhu5+FtFxdZ0a3Z2rdeaqq4RkT+Aw4DlbtliABF5F7jNLVsAHOuWdwX+b0/vVU5EbsHpFRiP82XlThHZT1Ufra6eJXmzi6ioKEZcfCEjLr7Q61CMMfVMWVEJuUtCthfaTKCLu3X6amAYcG5tKopIa2CTqhaISEPgEJwVYTcBPUWkqapuBI4B5rt1mqnqBnGefd6FM9I+UCOATFUtcK/9uvt5qk3y1l1vjDEmjDj7yQdy1ERVS3EGq32Jk4jfVdU/ROQ+ETkZQET6i8gqnGfoL7otdoAewM8i8hswFXhcVee6z+VHAt+LyBygLzDKrTNcRP4EFuD0GLy8F38wReUJ3v0sRUC1C+GAteQjUmFhIVlZWYgI/fv3Jy4uzuuQTIgtXLiQ1atX061bN1q1CnAhkQizfPlyFi9eTPv27enYsaPX4Zg9EMr95FV1MjC5Utndfr/PxOnGr1xvCs7c+aqu+QJVtNJVdTTO9LpgmCwiDVV1M4CIpOPsI18tS/IR5scffuSOf/6LlJImILA1NptHnh7FwIG1HkBq9iHbt2/nH3+7mm1Ls2kiKawq28T+Rwzg3lH3ERVVPzvqSktLue36m1mWNZ+M6DTW+/Jo3rsDj495yr7w7gOUPRgJV4+o6l2VXm+pbh/5cpbkI8j27du5/bq7OCBhCAlJiQAUlORzyzW388X3n5GYmOhxhCbYHh/1KE1WCke3PAQAVWXKdzP56IMPOeOsMz2Ozhuvjp9A0S+rGd5q59YbP8z7g+efeY7rbrzew8hMbUTHx5Ia6MC7Gjde3XeJyOk4z+XLp9Ip0F9ErgEmqOorVdWzJB9Bpk2bRnpJCxIa7EzmibFJpG1vxo8//shRRx3lYXQmFH74+ntGZBy947WIMKhpbz584716m+Q/ffdjzmp6YIWyAc268cbEzy3J7wPKCkvYsjhkA+/2ZaOAvwPly+Iqzqp8N+EMIqySJfkIUlpaClU8yxJflPOeiUhSaepvTFQ0pSUlHkXjPZ+vbJdHFdESha+sxjFKJgwolZaaM+XyVfW//gUiUqCq1XbZ18+HdhFq0KBB5ESvpdS3M6GXlJWwOXYdBx98sIeRmVDp078vizavrFCWtXEBJ5xxkkcRee+I449h9sZFFcrmZC9h0FHVbtZlwkioRtfv4wbVsqwCz1ry7iL7WcBqVT3Rnbf4NtAI58nKBe4uO/HAq8CBOPMR/7+9O4+PqjobOP57si9kIQlhh4CAgIiAUVGrUkRFtNJaq9QVpaK1FEuVaq1va7VaWnGtK4p1qaKorUWrooi2KsqqqAgCsmMggRBCCFnnef+4NxCSQMIwkzvL8/VzP8ycu8wz1wvP3HPOPeciVV3nUdghLSsrixv+byL33vEg2bVdEGBb7CZ+c/sNpKf7NSJixPruu+9YtWoVXbt2Deue1zff9jvGX3IV6wuKyJY2bNZi2vTO5aeXXuJ1aJ4Z/4truXb+OP7z3UI6Sjpb2cXu7Fget6r6sBEFCdsfb7pj4zf07nrliAAAIABJREFUfRGZpqoNR/cDvK2uvx7nOcW67PMX4D5VfVFEHgPGAY+6f+5Q1V4iMsbd7iIvAg4HPzz/h5x8ysm8N2cuMTHC6SNOJzs72+uwQobP5+O23/6eLz5YQKe4LIp8peT268o9jzxAUlKS1+EdspycHGb+5598+OGHbFy3gTHHHM3gwYNp+t+C6JCSksLTM59nwYIFrFz+DaP69OLEE0+M2qcNwk1sYjyZR3Twb+elgY0lxNx4kHX3HmiFqLb+wwruyEHPAHcCvwZ+ABQBHVS1RkROBG5T1bNEZLb7+hN3Sr8tQDs9SOD5+fm6aNGi4H8RE3ZefH4Gcx58mTM6H7e3bMHWr+kwagCTf3eTh5EZE15EZHEzQ8j6JS8lV3/Xx7/7uPFLHwpKTKFCRDrgTKijwEJVbbaHolc/be8HfsO+/hXZQIk7GhE44wvXjeixd9Ygd/1Od3tjDtm/XniZk3OP3q8sP7cvc99816OIjDH1KYJP/VsimVuT/Qlwvrt8KiLNzhrW6tX1InIuUKiqi0VkWF1xE5tqC9bVP+54YDxAt27dAhCpiURVVdXEpcTuVyYIvlrrz2tMqLDBcJp0C3CsqhYDiEgW8AHOvPIH5MWd/MnAeSKyDqej3XCcO/tMtzoe9p8ZaO+sQe76DKC44UFVdZqq5qtqfrt27YL7DUzYOvO8kSwu/Ga/shXF6zn2pOMOsIcxprVZ7/oD2lnvdUlLdmj1O3lV/S3uXLvunfyNqnqJiLwMXICT+K8A/u3uMst9/4m7fu7B2uONOZixP7uKn3/4KW9unE/nmCyKfLsozqjiiVv/7HVoxhjcjne9/Ox492VgYwkxbwJvi0jdnfslNBiDvymhNBjOTcCLIvIn4DNguls+HXhORFbj3MGP8Sg+EwGSk5N56sVnWLBgAcu++IrTe+Zx2rBhxMWF0l8FY6JXTWU1xau3eB1GyFHVm93m7mFu0QOqOqu5/Tz9l80dvecD9/Ua4PgmtqnAmfLPmICIiYlh6NChDB061OtQjDENqT0nfyCq+gbwxqHsY7cvxhhjQorPGmQbEZFS9nVEjwcSgN2qmnaw/SzJG2OMCSmW4xtT1f2GLRWRUYTysLbGGGNMQ7GJ8bQ9oqN/Oy8LbCyhTFXfFJG7gFsPtp0leWOMMSHF7uQbE5Ef13sbizOfS3lz+1mSN8YYEzJqKqvZZr3rm3JOvdc1wDpgdHM7WZI3xhgTUmwklMZU9Sp/9rNpmYwxxoQMxZnUxJ8lEonI6SIy6CDrk0XkFwdab3fyUUZVKSgoICEhgZycHK/DMcaYRrTJKUui1irgDhEZAHwMrAQqgFycdvkewIMH2tmSfBRZtmwZkyfeTPVOpVZr6HRELvc+NJXc3FyvQzPGGADiEuPJ9ndY2+WBjSUUqOoG4Ap3QpoRwAAgGWd+lz+r6kHnVbckHyV27drFdVdOpH/MibRJdR633Prtd1x71S949fWZiNgvZ2OM92oqqylaZR3vGnJnn5vpLi1mbfJR4q233iazsiNtEveNp9C+TSd2F1SxYsUKDyMzxph9VJ0R7/xZTGOW5KPE9qLtJGhSo/J4XxI7duzwICJjgmfunLlccNZ5jDphOOefcS5v/afZybpMyBDUz8U0Zkk+Snzv1JPZEVuwX1mtr5bSmG0cffTRHkVlTOB9Mm8eD958F+cmHM24Lqfzw+TBPPmH+5nz7hyvQzMtpOrf0hIiMlJEvhGR1SJycxPrTxWRJSJSIyIXNFj3toiUiMgbDcpFRO4UkZUislxEJrrlk0Xkc3f5SkRq3bb1VmNt8lFiwIAB5A8fyJI5H9Mpric1vmo2+lZyzQ1XkZZ20PkNjAkr0+5/hFHthpCWkAJAanwS53TIZ/qDjzHijBEeR2eaE5cYR46/He9WHny1iMQCDwNnAJuAhSIyS1W/rrfZBmAscGMTh7gbSAGuaVA+FugK9FVVn4jkAqjq3e4+iMgPgElu23qrsSQfJUSEP0+9i48++ojX//kfkpMzueXS6+jfv7/XoRkTUNu2FJLdrt9+ZRmJqezautOjiILrgw8+YPrDT7F923ZOPPVErp1wDe3atfM6LL9VV9ZQGLwR744HVrtTmyMiL+KMGrc3yavqOnddo0fvVfU9ERnWxHF/Dlysqj53u8ImtvkpMOMw4z9kluSjiIhwyimncMopp3gdijFBc0TfPmxcu5Wu6e33lm3dXUzHbp09jCo4Zs6YyRNTnmJg2nH0TBjAmjfWcOn7lzPz9RfJyMjwOjy/BXHEu87AxnrvNwEnBOC4RwAXiciPgCJgoqquqlspIinASGBCAD7rkFibvDEmovxi8vW8U/ol3+7YTK36WFdSwH+KP2Pib2/wOrSAqqmp4fEHpjE0axjpiRnESAx5mb3IKevISzMO6SmrkHMYI97liMiiesv4BoduqndeIH5SJAIVqpoPPAE81WD9D4CPW7uqHizJG2OAPXv28OILL3L9NRP5651/YcOGDV6H5LfevXvz0IzplBzThhd3zWdLv3ju+8djEdfBtLi4mARfMnEx+1fI5iZ35LP5n3kUVWAcRu/6baqaX2+Z1uDQm3Dazut0wRlU5nBtAl51X/8LGNhg/Rg8qKoHq643JuqVl5dzxYWXk7I1gbzUbmxeuo5xr13JXY9M4bjjjvM6PL/07NmTu+77q9dhBFVmZiYVlONTHzGy736teM82BvY90sPIDk9cYrz/He9WNbvFQqC3iPQANuMk34v9+7D9vAYMx7mDP416XQBFJMMtuzQAn3PILMkbE+Veeell0rYmM7idc/ORk5JN+8pc/nTLHbz2zr9tNMQQlZCQwA/HjGbus//l6Mx84mPj2V5exMa4b/nLFX/0Ojy/VVdWUxikEe9UtUZEJgCzceZkf0pVl4nI7cAiVZ0lIsfh3I23BX4gIn9U1aMARORDoC/QRkQ2AeNUdTYwBXheRCYBZcDP6n3sj4B3VHV3UL5UMyzJGxPl/vvOfzkivcd+ZWmJbajaVsmuXbtIT08/wJ7Ga7/81QQyMtKZ8fSLVO2uJq93dx677RE6dPDzTjhEBHOqWVV9E3izQdnv671eiFON39S+TfZaVtUS9p/vvf66p4Gn/Yv28FmSNybKteuQy84NpaQn7hsvwac+qqkmKanxKIkmdMTExDB23FjGjhvrcSSBFanTxnrBOt4ZE+Uuv/pylu75isqaKsCZjnjp9q8Ydvb3SUhI8Dg6E238He0umHf/4czu5E3UKSsr47VXXuWrRZ/Rs18ffjzmIrKzs70OyzP9+/fn13dO5r477yGxJpHy2j2cNOIkfvO7mzyJZ926dbz6wkts21rIySOGcdbIkcTHx3sSi2l98Unx5Pb2s7lhTWBjiQSiEfjzJz8/XxctOugUuyZKFRcXc/VFlzCwpg2909qxoayYT6qLuP/ZJ+nRo0fzB4hgPp+PLVu2kJmZSUpKiicx/O+D/zL1pj9waloPMhNT+XLnZnZ3T+PRp6dbog8xIrLYfS48oDokdtDLuvjXEX3qmnuCElM4s+p6E1WeeOgRTtYshnc+kq7pWZzcqRcXZPbivjv+7HVonouJiaFTp06eJXifz8c9t93JlV1PZGC77nRLz+GcrseQun4Xs99+25OYjDdsqtnAsSRvosqC/33MkNzu+5XlZbRj7YpmZrYIoJqaGjZu3EhZWVmrfWY4KCgoIKM2jtT4xP3KB2Z05qN35noUlWltehiLacza5E1UycjMYGdlOVnJbfaWVdbWEJvQOlXBb8yaxeN330+2JFJSXcHg007it7f/wTq4AWlpaeyqqWxUXlK5m5z2fTyIyHjF7soDx+7kTVT56fireO27ZdT4nId0fKq8uWkZoy++KOifvXTpUp698z4mdDqesV3zub7HydR8+g33TYnskdlaKj09nSOGHMXCwn29p3ZXV/K/XWv58SVjPIzMtDbrXR84didvosoZZ57Jlk3fcd/0Z8iNS6GoejfDzhvFFT+7KuifPWP605yd04fEOKfWQET4fqe+3PvWu9xwy83Exdlfx9vvnsLvJ9/Mo0s+JC0uiZ2x1Uya8oeo7xQZTeIS42nfx8/e9esDG0sksH9VTNhauXIlf/3DnXy3dgMxCXGMHnMB4669mpiYg1dQXXbVWC669GK2bNlCTk6OXx3NNm/ezNTbbmft18shNoZh54ziukm/Omi1+/bCbZRUKHd//Aa7KiqQGDg1ry8JxFBVVWVJHkhJSWHqww+yc+dOSktL6dy5c7P/P01kqa6spmBl0OaTjzr2r4oJS1u3bmXS2Gs5N/0YftC1N9W1Nbz3j9mUle5i0s03Nrt/QkIC3bp18+uzd+/ezcTLruCizA5c2T+fGp+P2e99yG0FW7jr/nsPuF+fYwbwj6mPc2mnk8nOTqPSV83b65eyNbvWsx7toSojIyOs50M3h8eq3gPHfiKbsDTz+RkcF9eNTmnOIDbxsXGc2WkQ7772HyoqKoL62W++/gb5ccn0zc4FIC4mhnO692bd4iUUFhYecL/S4hJOyupHPHFU1dZSW+sjP6M38THx+Hw2kKcxdQ5jPnnTgCV5E5bWfvMtndrsP0pdjMTQNi6V4uLioH72htWr6Zac1qi8c2IKW7YcuJpx87qNHHvkQFJys9iTIEh6Ckf07kXb5PSg/zAxJmzYsLYBZdX1JiwNPH4Iy5e9S/vUtnvLqn017PDtoV27dkH97AHHDuHD9z/m6Nx9nYN8qqzZU0ZeXt4B9xt4/BDWvP4Fg3N707ZtFgB7qiupTYwhOTk5qDEbEy7ikuLp4O+wthsDG0sksCRvwtKPL7yAy2e8QnLhao7O6cGOijLeKVzK2ElXB3340+EjRvCPx59gzoZvObVzd3ZVVfHPDas4/YLzDzot66VXXs6Vr19MwrY4+mZ3p2h3Ce9s/5zr/nijzdlujKu6wjreBZJV15uwlJaWxvSZz9HmrH68ULqQ+RmF/PKeW7nw4uA/Tx0fH89jzz9HyjkjmFrwLc9X7WDUzTfw819df9D9cnJyePLl55BTu/FC6Xy+bF/KmBvH8+lHnzL5lzcyZ84ca5s3BmuTDySboMYYj7z0wks8PXU6RyUdSWJsIqv2rKHL0Dzu+ds9dmdvQl6wJqhpl9BBR+f6N0HN9M02QU1DVl1vjAfKy8t58v5pnJM7grgY569hx7T2zJ3/IZ999hlDhgzxOEJjvGHj0AeWJXljPPD111+TG5O9N8HX6RrTiY8++MiSvIla8UnxdPR3xLvNgY0lEliSN8YDbdu2pVz3NCov1z3k5OZ4EJExoaG6oprN31jHu0CxjnfGeOCII44gpUsa63fue+ZnZ0Up62I2cc4PzvEwMmO8Z8/JB44leRNQtbW1LF68mHnz5tkAL8342xMPUdKznDe3zeGd4v+yMHYp9z5xvw3naqKaYr3rA8mq603ArFy5kklXTySzMpU4YilgOzfdcQunn3G616GFpKysLJ547kl27NhBZWUl7du3t171xqgzuJQJDEvyJiB8Ph+Txl/PaQlDyMrIBKCypoopt9zJMYOPISfH2pkPpG3bts1vZEyUiE+Ko/ORfna8s6b8RizJu0pKSqioqLC7KT99/vnnZFSkkJWeubcsMS6BPrFdmf3W21xymX/PvRpjokt1RQ2brONdwER9ki8pKWHy9Texdtla4iWB2HThT/fcwaBBg7wOLaxUVVURJ7GNyuOIpbKi0oOIjDHhSAGf1dYHTNR3vJt4zfXo1/GcknEmQ9OH0a/qWCZdfUPQZzKLNIMGDWKLbqeiZl9C96mPVbUbOf3MER5GZowJL+r3f6axqE7yGzZsoPDb7XRNz9tb1iahDR1rujLrtVneBRaGkpKS+O1dt/Lvorks3voVS7cu518F73H+uAvp3r271+EZY8KIT/1bWkJERorINyKyWkRubmJ9ooi85K6fLyJ5bnmCiPxdRL4UkaUiMqyJfWeJyFdNlN8oIioird45Kaqr63fs2EGyNJ7iMzm2DYVbioLymevXr+e2/7uDFV+vJi4+lgsu/CETJl5HbGzjqu5w8/3h3+eY2cfwztuzqaioZPKI4XTr1s3rsIwxYSQ+KZ4uff3seFd48NUiEgs8DJwBbAIWisgsVf263mbjgB2q2ktExgB/AS4CrgZQ1aNFJBd4S0SOU1Wfe+zzgbImPrOr+3kb/PtShyeqk3yfPn0okWJqfDX7DS+61beZq4cHvqNYSUkJl465iuTKPLqkDsXnq2Xm3+dQVLSNP931x4B/nheysrIYc/FPvQ7DGBOmqiuq2bgiaB3vjgdWq+oaABF5ERgN1E/yo4Hb3NevAA+J0xu7P/AegKoWikgJkA8sEJE2wK+B8cDMBp95H/Ab4N/B+ELNierq+uTkZCZMvo55O+ayqXQDRbsLWVL8CXn5XTj++OMD/nmvvvIv2N2WjNR2AMTExNIxvS9z3v4vpaWlAf88Y4wJN07HO/VrAXJEZFG9ZXyDw3cGNtZ7v8kta3IbVa0BdgLZwFJgtIjEiUgP4Figq7vPHcA9QHn9A4nIecBmVV16mKfFb1F9Jw/w45/8mH5H9eOlf7zErp1l/GL0NYwYMYKYmMD//vlm+UqS4/YfzUxEiI9JZevWraSnpwf8M40xJtwcRhe6bc1MNdvU89ENP+5A2zwF9AMWAeuBeUCNiAwCeqnqpLr2ewARSQF+B5zZ4uiDIOqTPED//v35YytUlx839FjmvfMcGbTbW6bqo4YyOndu+GPSGGOiUxAHvNvEvrtvgC7AdwfYZpOIxAEZQLGqKjCpbiMRmQesAk4DjhWRdTg5NVdEPgB+CfQAlrpjr3QBlojI8araagMBRHV1fWs799xzaNOhlq0la/D5aqmo2s2GkiVcMe4SUlJSvA7PGGM854xdr34tLbAQ6C0iPUQkARgDNHyUahZwhfv6AmCuqqqIpIhIKoCInAHUqOrXqvqoqnZS1Tzge8BKVR2mql+qaq6q5rnrNgFDWjPBgwd38m5Pw2eBDjhzCkxT1QdEJAt4CcgD1gEXquoOt8PDA8AonPaOsaq6pLXjDoTk5GReeuV5Hn/0Cd59Zy4Zuencce1kzjrrrIB9hqqycOFCXpv5GiCcP+ZH5OcfrPbKGGNCR0JiPF38HdZ228FXq2qNiEwAZgOxwFOqukxEbgcWqeosYDrwnIisBopxfggA5AKzRcSHM3P9Zf4F2bpEW3kiABHpCHRU1SUikgYsBn4IjMWpEpniPrvYVlVvEpFRONUeo4ATgAdU9YSDfUZ+fr4uWrQoqN8jVE2dMpW5Mz+gd2JvVGFV5UpGXnoW198w0evQjDERREQWN9P+7ZfMuPZ6SuaY5jdswhvbHwxKTOGs1avrVbWg7k5cVXcBy3F6M44GnnE3ewYn8eOWP6uOT4FM94eCaWDjxo288/I7nJp9Gp3SOtM5vTOn5pzG6y+8TkFBgdfhGWNMi6ifi2nM0zZ5tyfiYGA+0F5VC8D5IYBTNQIte+QBERlf99hEUVFwBrIJdYsXL6ZdbYf9JtiJkRja+zqwZElYtnAYY6JMkNvko45nSd4dPOBV4FeqerCHxFvyyAOqOk1V81U1v127dk3sEvnatm1LVWxFo/LKuEoyMzOb2MMYY0KNfwneknzTPHmETkTicRL886r6T7d4q4h0VNUCtzq+boDCljzyYICTTjqJP6dOYVv5NnJSnCGSi3YXUd6mjKFDh3ocnTHGNC8hKZ5u/na8mxfYWCKBF73rBaf34nJVvbfeqrrHFqa4f/67XvkEd/jBE4CdddX6Zn/x8fE89uyjTJ4wmWWFX+BTJaNTBo/97bGIGBvfGBP5KiuqWRe8YW2jjhd38ifjPHrwpYh87pbdgpPcZ4rIOJyB/H/irnsTp2f9apxH6K5s3XDDS15eHi+/8TKFhU5FSG5ubjN7GGNMaLFpYwOn1ZO8qn5E0+3sAKc3sb0CvwhqUBHIkrsxJlxZ+3rg2LC2xhhjQkZd73oTGJbkjTHGhIzEpHjy/J1P3jreNWJJ3hhjTMiorKhm7QrrWx0oluSNMcaEEMUnVl0fKJbkjTHGhAynTd7ndRgRw5K8McaYkGKP0AWOJXljjDEhw+l4598cZPM+DnAwEcCSvDHGmJBRWVHF2hWbvQ4jYliSN8YYEzoE63gXQJbkjTHGhAzreBdYluSNMcaEELUkH0CW5I0xxoSMxKQEevTt5Ne+iz8KcDARwJK8McaYkFFRUcW3KzZ5HUbEsCRvjDEmpKhYdX2gWJI3xhgTQqxNPpAsyRtjjAkZ1rs+sCzJG2OMCSGKUut1EBHDkrwxxpiQkZiUQM++Xfzad1lRgIOJAJbkjTHGhIzKiipWr9jodRgRI8brAIwxxpg6iuKj1q+lJURkpIh8IyKrReTmJtYnishL7vr5IpLnll8iIp/XW3wiMkhEUkTkPyKyQkSWiciUesf6tYh8LSJfiMh7ItI9QKepxSzJG2OMCSmKz6+lOSISCzwMnA30B34qIv0bbDYO2KGqvYD7gL8AqOrzqjpIVQcBlwHrVPVzd5+pqtoXGAycLCJnu+WfAfmqOhB4Bfjr4ZwXf1iSN8YYE0I0aEkeOB5YraprVLUKeBEY3WCb0cAz7utXgNNFRBps81NgBoCqlqvq++7rKmAJ0MV9/76qlrv7fFpX3pqsTd4YY0zISExKoFffrn7tu6r5jnedgfoN/puAEw60jarWiMhOIBvYVm+bi2j84wARyQR+ADzQxGePA95qNsIAsyRvos6C+fO5/8/3smNrMSkZqVw98VpGjhrpdVjGGJyOdytXrPd39xwRWVTv/TRVnVbvfcM7cnAezael24jICUC5qn61304icTh39w+q6poG6y4F8oHTmv8KgWVJ3kSVL774gt9fdwsjs4eS2S6d3VXlPPb7+wEs0RsTAvTwnpPfpqr5B1m/CahfTdAF+O4A22xyE3cGUFxv/RjcqvoGpgGrVPX++oUiMgL4HXCaqla26FsEkLXJm6jy+H2PMCxjMJlJ6QCkJqRwRs4JPPHAox5HZoypE8Q2+YVAbxHpISIJOAl7VoNtZgFXuK8vAOaqqgKISAzwE5y2/L1E5E84PwZ+1aB8MPA4cJ6qFh7SSQgQS/ImqmxYt572qTn7laUmpLB7526PIjLG7E/xaa1fS7NHVq0BJgCzgeXATFVdJiK3i8h57mbTgWwRWQ38Gqj/mN2pwKb61fEi0gXnTr0/sMR9vO5n7uq7gTbAy255wx8UQWfV9SaqHHlUPzZ+WUC3jH3zVZdUlNI2N8vDqIwxdRKTEund17/Hyde3YD55VX0TeLNB2e/rva7AuVtvat8PgKENyjbRdDs+qjqi+YiCy5K8iSrX/XoC1475Gd8DuqZ3ZOvu7XywczF/fPgur0MzxgCVFZWsXLHW6zAihlXXR7gdO3Yw7bEnmDzxNzz3zD8oKyvzOiRP9ezZk8defJLSo2P5Z9n/2NCtlL/+/X5OGDq0+Z2NMUFX1/HOn8U0Jm5/goiSn5+vixYtan7DCLdhwwbGjfkZuRWdyU5sR2HlFkozi/nHK8+SlWXV08YY/4nI4mZ6svslIa6N5mYO9Gvfzds/CUpM4czu5CPYXbf9md61Azgy6yhyUnPpnzWQ3J2deeRB60lujAldQexdH3WsTT6CrfxqJcPSz96vLC/zCD7+4ANvAjLGmGYkJiXQ58g8v/b9bt78wAYTASzJR7D4xHiqa6uJj43fW1ZeXU5623QPozLGmAOrrKjkmxVrmt/QtIhV10ewCy+7kK9KllDX78KnPr7cuZix11zRzJ7GGOMNBXx+/mcaszv5CHblz8ayrbCIt1+bTZvYNMq0lIt/PoaRZ9vwrSZwCgoKWLduHd27d6dTp07N72DMQam1rweQJfkIFhMTw0233sSESRMoKiqiQ4cOJCUleR2WiRC1tbX8329u5YuPlpAb25ai2h0cddIx/Onuu4iLs39ajP+0BaPXmZaxv4lRIDU1ldTUVK/DMBHmmelPs+XDNZzfYfjesnnzPuOpadMZf901HkZmwllSUiJH9u3p175Fn3wR4GjCn7XJG2P88u+X/sVxOQOo8dVSWllGja+W43IG8MYrrT48t4kg6lbX2yN0gWF38sYYv1RVVbFo25cs2rSMFEmlXHczpHM/qtOrvA7NhLGKikpWLF/tdRgRw5K8McYvuV068MXbXzM880xiJRaf+vho44ccMbyX16GZsKbWJh9AluSNMX4p2VZK37SjKa+pIF7iqNEa+qYNoLhkm9ehmTBnVe+BY0neGOOXPeV7OKHPSZSVlVGxp4L05CS6tWnD++VzvA7NhDMFVUvygWJJ3hjjl4HHHs3m+ZvpmtGV9HRnFMXNpZs4anB/jyMz4SwpKZG+ff1r8vl0wbcBjib8WZI3xvhl0k2TuPLCqygv3k375PYU7ilkY8J6pv/2Sa9DM2GsoqKC5StWeh1GxLBH6EKUqjJnzhwuuuhizj77XB566BHKy8u9DsuYvTp37syMWS9w7KWDKO29g8GXDGTGrBfo1q2b16GZMKbYLHSBZPPJh6iH/vYwTzz+HClJucTFxVNaVkT3ntm8+upM4uPjmz+AMcYEUbDmk4+NTdA2qR382rd010abT74Bu5MPQWVlZfz978/RNr07SYkpxMXGk5XRiY3rtjFnjnVqMsZEOp+fi2nI2uRD0Jo1axASEZH9ymNI4tNPF3L22WcfYE9jjAlvTse73n7tu2Dh5gBHE/4syYegDh06oNp41LBaraR3b//GdDbGmHBQUVHJ8uXfeB1GxLDq+hCUm5vLiScfx47S7/bOBb+7fCepacLo0aM9js4YY4LJxq4PJEvyIeree6dy/gVnsLtqAzvL19K7Xw4vzHiWtLQ0r0Mzxpggszb5QLHe9WFAVRu1zxtjjJeC17s+TlOSM/3at2z3dutd34C1yYcBS/DGmGiRlJRE335H+rXvokXzAhxN+As+X0nFAAAIG0lEQVSbJC8iI4EHgFjgSVWd4nFIxhhjAqyiooLly5cH7fjN5RIRSQSeBY4FtgMXqeo6EckDlgN1vQI/VdVr3X2OBZ4GkoE3getVVUUkC3gJyAPWAReq6o6gfbkmhEWbvIjEAg8DZwP9gZ+KiA2QbYwxESk4bfItzCXjgB2q2gu4D/hLvXXfquogd7m2XvmjwHigt7uMdMtvBt5T1d7Ae+77VhUWSR44HlitqmvUebbsRcC6mRtjTMRRUD+X5rUkl4wGnnFfvwKcLgdpMxWRjkC6qn6iTie3Z4EfNnGsZ+qVt5pwqa7vDGys934TcIJHsRhjjAkSp02+r1/7Ll68oLlNWpJL9m6jqjUishPIdtf1EJHPgFLgVlX90N1+U4NjdnZft1fVAvdYBSKSe2jf6PCFS5Jv6lfUfj/bRGQ8TnUJQKWIfBX0qMJLDrDN6yBCkJ2XxuycNGbnpDH/esc1o7y8fPbixQty/Nw9SUTqP1o1TVWn1XvfbC45yDYFQDdV3e62wb8mIke18JieCZckvwnoWu99F+C7+hu4/yOnAYjIInuMYn92Tppm56UxOyeN2TlprEEyDRhVHdn8Vn5rNpfU22aTiMQBGUCxWxVf6ca4WES+Bfq423c5wDG3ikhH9y6+I1AY6C/UnHBpk18I9BaRHiKSAIwBZnkckzHGmPDSklwyC7jCfX0BMNftKd/O7biHiPTE6WC3xq2O3yUiQ922+8uBfzdxrCvqlbeasLiTd9tFJgCzcR57eEpVl3kcljHGmDByoFwiIrcDi1R1FjAdeE5EVgPFOD8EAE4FbheRGqAWuFZVi911P2ffI3RvuQvAFGCmiIwDNgA/CfZ3bCgiR7wTkfEN2mGinp2Tptl5aczOSWN2ThqzcxIeIjLJG2OMMSZ82uSNMcYYc4giLsmLyEgR+UZEVotIq48u5BUR6Soi74vIchFZJiLXu+VZIvKuiKxy/2zrlouIPOiepy9EZIi33yB4RCRWRD4TkTfc9z1EZL57Tl5yO+AgIonu+9Xu+jwv4w4WEckUkVdEZIV7vZwY7deJiExy/958JSIzRCQp2q4TEXlKRArrP37sz3UhIle4268SkSua+izTeiIqyUt0D39bA9ygqv2AocAv3O9+oGEVz2bfEIzjcYZljFTX44w5XecvwH3uOdmBM4wlHHw4y0jyAPC2qvYFjsE5N1F7nYhIZ2AikK+qA3A6ZI0h+q6Tp9k3HGudQ7ouxBmr/Q84A8wcD/yh7oeB8UZEJXmiePhbVS1Q1SXu6104/3B35sDDKo4GnlXHp0Cm+xxnRBGRLsA5wJPuewGG4wxXCY3PSYuHswxHIpKO00t4OoCqVqlqCVF+neA8aZQsznPRKTgDn0TVdaKq/8PpTV7foV4XZwHvqmqxOxHLuzT+4WBaUaQl+aaGLOx8gG0jllt9OBiYT4NhFYG6YRWj5VzdD/yGfbNXZAMlqlrjvq//vfcbzhKoP5xlpOgJFAF/d5swnhSRVKL4OlHVzcBUnEecCnD+vy8muq+TOod6XUT89RJuIi3Jh/Twgq1BRNoArwK/UtXSg23aRFlEnSsRORcoVNXF9Yub2FRbsC5SxAFDgEdVdTCwm4PPjBXx58StTh4N9AA6Aak41dENRdN10pwDnQM7NyEm0pJ8S4YsjFgiEo+T4J9X1X+6xVvrqldl/2EVo+FcnQycJyLrcJpuhuPc2We61bKw//fee06k3nCWrRlwK9gEbFLV+e77V3CSfjRfJyOAtapapKrVwD+Bk4ju66TOoV4X0XC9hJVIS/JRO/yt2yY4HViuqvfWW3WgYRVnAZe7vWSHAjvrquUihar+VlW7qGoezrUwV1UvAd7HGa4SGp+TRsNZtmLIQaeqW4CNIlI3ucjpwNdE8XWCU00/VERS3L9Hdeckaq+Teg71upgNnCkibd0akjPdMuMVVY2oBRgFrAS+BX7ndTyt+L2/h1Mt9gXwubuMwmkrfA9Y5f6Z5W4vOE8ifAt8idOz2PPvEcTzMwx4w33dE1gArAZeBhLd8iT3/Wp3fU+v4w7SuRgELHKvldeAttF+nQB/BFYAXwHPAYnRdp0AM3D6JFTj3JGP8+e6AK5yz81q4Eqvv1e0LzbinTHGGBOhIq263hhjjDEuS/LGGGNMhLIkb4wxxkQoS/LGGGNMhLIkb4wxxkQoS/LGBIE4swKudSfswH1ueK2IdBeRjuLOiHcIx5sqIsODE60xJlJZkjcmCFR1I87MXFPcoinANFVdD/waeOIQD/k3Dj78rDHGNGLPyRsTJO4ww4uBp4CrgcGqWiUia4B+qlopImNxZvaKBQYA9wAJwGVAJTBKVYvd4y0GzlFn1DpjjGmW3ckbEyTqjIM+GWfO8V+5Cb4HzlzklfU2HQBcjDNV8p1AuTqTx3wCXF5vuyU44/EbY0yLWJI3JrjOxhkqdID7viPOVK/1va+qu1S1CGfa0tfd8i+BvHrbFeLMkmaMMS1iSd6YIBGRQcAZwFBgkjuL1x6csc/rq39X76v33oczNWydJHd/Y4xpEUvyxgSBO5vZozjV9BuAu4GpOJMn5fl52D44E6gYY0yLWJI3JjiuBjao6rvu+0eAvkA+8K2I9DqUg7md+HrhzB5njDEtYr3rjWllIvIj4FhVvfUQ9xmiqv8XvMiMMZEmrvlNjDGBpKr/EpHsQ9wtDufxOmOMaTG7kzfGGGMilLXJG2OMMRHKkrwxxhgToSzJG2OMMRHKkrwxxhgToSzJG2OMMRHKkrwxxhgTof4fDlkvpOg6lngAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.subplot(111)\n",
    "GSLIB.locmap_st(df,'X','Y','Porosity',xmin,xmax,ymin,ymax,pormin,pormax,'Well Data - Porosity','X(m)','Y(m)','Porosity (fraction)',cmap)\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.1, wspace=0.2, hspace=0.2); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A Very Simple Bootstrap Method in Python                       \n",
    "\n",
    "If you are new to bootstrap and Python, here's the most simple code possible for bootstrap.\n",
    "\n",
    "* specify the number of bootstrap realizations, $L$\n",
    "* declare a list to store the bootstrap realizations of the statistic of interest\n",
    "* loop over L bootstrap realizations\n",
    "    * n MCS, random samples with replacement for a new realization of the data\n",
    "    * calculate the realization of the statistic from the realization of the data\n",
    "* summarize the resulting uncertainty model, histogram, summary statistics etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeYAAAF6CAYAAADMN/v3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZicVZn38e9tQEAWUQkIARJERMHBiNEXd9yRV0RHEdxARNAZ13ccZ0SdcWVAR1xxQ0UWRQEZERRHEEVUFAiIQAQVJQghhIAKiAgE7vePc4pUmuqkOt3VddL9/VxXXV116qlT96mnun71LPU8kZlIkqQ23G/YBUiSpOUMZkmSGmIwS5LUEINZkqSGGMySJDXEYJYkqSEGs9ZoEfGuiPjSsOvomOh6IuKvEfGwev3oiPjQBPb9+Yj4j4nqr6vfiIivRMSfI+L8ie5fE6v7PaY2GMwiIjIiHj6i7X0R8dUh1HJ2RLyu3+kz878ys6/pxzumWtvfI+LWiLglIi6MiHdGxDpjraffcWbmBpn5h9Wtuev5XhMRPx3R9xsy84Pj7buHpwDPAbbMzCdMVKcRsX4NkdMnqs9hq/97t9VxLYqIj0XEjMmsofs9NtFf/rR6DGY1oS5lrQnvxzdl5obA5sDbgX2A0yMiJvJJImKtiexvks0GFmbmbWN94CrG/VLgDuC5EbH56hY3jucflMdk5gbAs4BXAAeOtYPJDnMNWGZ6meYXIIGHj2h7H/DVen1X4FpKEN0ALAb275p2PeBw4GrgZuCnwHr1vl2Ac4G/AL8Cdu163NnAIcDPgNuBrwF3A38H/gocUaf7JHANcAtwIfDUUeqcU8eyH/BH4Ebg3fW+3YA7gbtq378C9gIuHDHutwOnjPI6nQ28bkTb1sDfgBf0qGdd4KvATXX8FwCb1TH3GmcCbwR+B1w1ct4ARwOfB84EbgV+DMweMfa1RtYLPKo+1931+f7S1d+HuqY/ELgS+BNwKrDFiPfIG2ptfwY+A0SP1+iAEc/1/j77XmHco7z+P6yv3UXAv3a1vxP45ohpPwl8ql5/IPBlyvt2EfAhYEa97zWU99/Ha20fAratz3UT5T30NWDjrr53Bn5Z58FJwAkjXscXABfXeX4usFO//3u1v8774VF1Hv4FWAC8sGu6o4HPAacDtwHPruM8FlhK+V98D3C/Ov3D6/vl5jqmE0bWABxE+f+4s86704B3ACePqPnTwCeG/bk1lS9DL8DL8C8jPxxq2/tYMZiXAR8A1gZ2p4TRg+r9n6kfILOAGcCTgHXq7Zvq9PejrN68CZhZH3c2JUB3BNaqfZ/NfcPvVcBD6jRvB64H1u1R55w6li9Sviw8hrKE9aiR09bb61A+jB/V1fZL4CWjvE73qa22nwN8uEc9r68fbg+or8vjgI1G66vWfibwYJZ/sRkZzLcCT6u1fxL46Yix3yeY6/XXdKbtuv9oaqAAz6R8YO9c+/40cM6I2r4DbEz5MrIU2G2U12mF5+qz7xXG3aPPrYF7gB3qe+CSrvtmU96Pndd2BiWEd6m3TwG+AKwPbAqcD7y+q9ZlwJsp76/1KCH1nFrrzDp/P1Gnvz8l9N5Keb/+IyXIOq/jzpQvr/+n1rEfsBBYZ1X/e3Vs11O+3KxN+SLzrvqcz6zzfvuueXcz8GTK/9a6lFD+NrBhfT/8FjigTv914N1d0z5llBrufU/U25tTgn/jenutOr7HDftzaypf1oRVh2rDXcAHMvOuzDyd8o16+7r6+bXAWzNzUWbenZnnZuYdlEA9PTNPz8x7MvNMYD4lqDuOzswFmbksM+/q9cSZ+dXMvKlOczjlA3P7ldT6/sy8PTN/RVkyfswo/d5BWdp5FUBE7Ej5QPtOvy9KdR0lVEa6i/KF4uH1dbkwM29ZRV+HZuafMvP2Ue7/bmaeU2t/N/DEiNhqjPX28krgqMy8qPZ9cO17Ttc0h2XmXzLzj8CPgLkT2Peqxr0vJYx/TQmZHSPisQCZeTVlKfpFddpnAn/LzF9ExGbA84G3ZeZtmXkDZel4n66+r8vMT9f31+2ZeWVmnpmZd2TmUuBjwNPrtLtQwulT9X/hfyhB33Eg8IXMPK/O82MoXw53Wcnrc1FE/JnyJe5LwFfq9BtQXvM7M/OHlPfly7se9+3M/Flm3kN5r+0NHJyZt2bmQsparFfXae+ifIHZIjP/npkr7G8wmsxcTPlisldt2g24MTMv7OfxWj0Gs6Csdlx7RNvalH/mjpsyc1nX7b9RPjg2oXwD/32PfmcDe0XEXzoXyo5B3dsHr1lVcRHx9oi4PCJurn08sD7vaK7vUedojgFeUbcRvxo4sYbHWMyiLHmPdBzwfeAbEXFdRHwkIka+ziOt6vW49/7M/Gt93i3GUuwotqAsCXb3fRNlbB1jeV3H2veqxr0vZZUymXkdZbXsfl33H8/y0HpFvQ3lPbg2sLjrPfgFypJzz+eOiE0j4ht1Z6xbKJsjOu+3LYBFmZmjPH428PYR7/mtWPk82jkzH5SZ22bme2rQbgFcU693XM3or9kmLF+a7zX9vwEBnB8RCyLitSupZ6RjqF9e69/jxvBYrQaDWVBWJ88Z0bYNK/6Tj+ZGyjbFbXvcdw1wXGZu3HVZPzMP65pm5OnNVrgdEU8F/h14GWXV+caUVXirs7PVfU6llpm/oKyKfCrlA31MHzp1afVxwE969H1XZr4/M3egrN5/ASVgetayivaOe5eOI2IDypL6dZTVjVBWm3c8dAz9XkcJlU7f61OW9het4nH96KfvUeuLiCcB2wEHR8T1EXE9ZVXxy7t21joJ2DUitgRezPJgvoayxLpJ13two8zccSXPfWht2ykzN6KEUef9thiYNWJnv+41FtcAh4x4zz8gM78+2vhGcR2w1YgdIrdm9NfsRpYvFd9n+sy8PjMPzMwtKJtYPjvylxg9+uw4BdgpIh5NeQ9/bYxj0RgZzIKyOvc9EbFlRNwvIp4N7AF8c1UPrN/ojwI+FhFbRMSMiHhi/QnRV4E9IuJ5tX3diOh8eI5mCdD9m8oNKdsAlwJrRcR/Ahut3jBZAszpsff3scARwLJ+V/FFxAMi4umUbXrnU3bCGTnNMyLiH+oes7dQPjjv7qpldX47untEPCUi7g98EDgvM6+pq1wXAa+qr/VrWfHL0hJgy/q4Xo4H9o+IuXXe/Vfte+Fq1DjRfe9H2Qa9A2X1+Vzg0ZQvIc8HqOM/m7Ia+KrMvLy2LwbOAA6PiI3q+3vbOu9GsyF1J7mImEXZAarj55R5+KaIWCsi9gS6fxL2ReANEfF/6i8N1o+I/xsRG/Y51o7zKF+2/i0i1o6IXSn/k9/oNXFm3g2cCBwSERtGxGzgXyj/g0TEXl3/d3+mBPDdPbq6z/syM/9O+Sw4Hji/bsrQABnMgrJT17mUvan/DHwEeGVmXtbn4/8VuJSy1/GfgA9T9ga9BtiTsgPLUsrSxDtY+fvuk8BLoxyc4lOUVcHfo+zIcjVl6XyVq79HcVL9e1NEXNTVfhzlg76fpeUjIuJWygfYJ4CTKTtB3dNj2odSPtBuAS6nrH7t/I565Dj7dTzwXsrr/DjK9tuOAymv702UHerO7brvh5Q9e6+PiBtHdpqZZwH/UcezmBLq+4ycbnWMp++IWJeytuTTdamvc7mKMr9Grs5+NsuXljv2pazm/TXl/f1NVtycMtL7KTtx3Qx8F/ifrrHcSdnh6wDK3tKvomz7vaPeP58yH46oz3UlZQezManP80LKF48bgc8C+2bmFSt52JspYf4Hyv/y8ZQvzQCPB86LiL9S9op/a30NR/oysENdDX9KV/sxwD/gauxJEStuKpGmn4hYj7Kn6c6Z+bth16M1S0ScB3w+M78y7FoGJSK2Bq4AHtrHDowaJ5eYJfgn4AJDWf2IiKdHxEPrquz9gJ2A/x12XYNSN/38C/ANQ3lyrMlHF5LGLSIWUnbsedEqJpU6tqdsz92A8muEl9Zt2VNO3VFvCWUz0m5DLmfacFW2JEkNcVW2JEkNMZglSWrIGr2NeZNNNsk5c+YMuwxJksbkwgsvvDEzZ/a6b2DBXI+IdCzlt5z3AEdm5icj4n2U3/ktrZO+K8uxl4mIgym/D7wbeEtmfn9lzzFnzhzmz58/oBFIkjQYETHqkRUHucS8DHh7Zl5Uj3pzYUScWe/7eGZ+dESRO1AOOrAj5TixP4iIR9Qj2kiSNC0MbBtzZi7OzIvq9VspRz6atZKH7En5ndwd9Yg0V7Lioe4kSZryJmXnr3p6t8dSjv8K5Tizl0TEURHxoNo2ixUPtXgtKw9ySZKmnIEHcz0DzsmU86HeAnyOcqzcuZTj5h7embTHw+/zI+uIOCgi5kfE/KVLl/Z4iCRJa66BBnM99+zJwNfqCcXJzCX1BOL3UM7E0lldfS0rnj5tS8qpz1aQmUdm5rzMnDdzZs8d2iRJWmMNLJjr+Uq/DFyemR/rau8+q8uLgc4ZjE4F9omIdSJiG8r5V88fVH2SJLVokHtlPxl4NXBpRFxc295FObn5XMpq6oWUk3aTmQsi4kTKqdmWAW90j2xJ0nQzsGCuJ5zvtd34PieU73rMIcAhg6pJkqTWeUhOSZIaYjBLktQQg1mSpIYYzJIkNcRgliSpIWv0aR+lNd0ez30aixeNepKZodp81mxOO+OcYZchTTsGszREixddzfyD2zyC3bxD2/zCIE11rsqWJKkhBrMkSQ0xmCVJaojBLElSQ9z5S1Ney3s+37BkCdDmzl+ShsNg1pTX8p7Ps968aNglSGqMq7IlSWqIwSxJUkMMZkmSGmIwS5LUEINZkqSGGMySJDXEYJYkqSEGsyRJDTGYJUlqiMEsSVJDDGZJkhpiMEuS1BCDWZKkhhjMkiQ1xGCWJKkhBrMkSQ0xmCVJaojBLElSQwxmSZIaYjBLktQQg1mSpIYYzJIkNcRgliSpIQazJEkNMZglSWqIwSxJUkMMZkmSGmIwS5LUEINZkqSGGMySJDXEYJYkqSEGsyRJDTGYJUlqiMEsSVJDDGZJkhpiMEuS1BCDWZKkhhjMkiQ1xGCWJKkhBrMkSQ0xmCVJasjAgjkitoqIH0XE5RGxICLeWtsfHBFnRsTv6t8H1faIiE9FxJURcUlE7Dyo2iRJatUgl5iXAW/PzEcBuwBvjIgdgHcCZ2XmdsBZ9TbA84Ht6uUg4HMDrE2SpCYNLJgzc3FmXlSv3wpcDswC9gSOqZMdA7yoXt8TODaLXwAbR8Tmg6pPkqQWTco25oiYAzwWOA/YLDMXQwlvYNM62Szgmq6HXVvbRvZ1UETMj4j5S5cuHWTZkiRNuoEHc0RsAJwMvC0zb1nZpD3a8j4NmUdm5rzMnDdz5syJKlOSpCYMNJgjYm1KKH8tM/+nNi/prKKuf2+o7dcCW3U9fEvgukHWJ0lSawa5V3YAXwYuz8yPdd11KrBfvb4f8O2u9n3r3tm7ADd3VnlLkjRdrDXAvp8MvBq4NCIurm3vAg4DToyIA4A/AnvV+04HdgeuBP4G7D/A2iRJatLAgjkzf0rv7cYAz+oxfQJvHFQ9kiStCTzylyRJDRnkqmxJa7Drr1/CvB1nD7uMnjafNZvTzjhn2GVIA2EwS+op71nG/IPb/EnivEOvHnYJ0sC4KluSpIYYzJIkNcRgliSpIQazJEkNMZglSWqIwSxJUkMMZkmSGmIwS5LUEINZkqSGGMySJDXEYJYkqSEGsyRJDTGYJUlqiMEsSVJDDGZJkhpiMEuS1BCDWZKkhhjMkiQ1xGCWJKkhBrMkSQ0xmCVJaojBLElSQwxmSZIaYjBLktQQg1mSpIYYzJIkNcRgliSpIQazJEkNMZglSWqIwSxJUkMMZkmSGmIwS5LUEINZkqSGGMySJDXEYJYkqSEGsyRJDTGYJUlqiMEsSVJDDGZJkhpiMEuS1BCDWZKkhhjMkiQ1xGCWJKkhBrMkSQ0xmCVJaojBLElSQwxmSZIaYjBLktQQg1mSpIYYzJIkNWRgwRwRR0XEDRFxWVfb+yJiUURcXC+7d913cERcGRG/iYjnDaouSZJaNsgl5qOB3Xq0fzwz59bL6QARsQOwD7BjfcxnI2LGAGuTJKlJAwvmzDwH+FOfk+8JfCMz78jMq4ArgScMqjZJklo1jG3Mb4qIS+qq7gfVtlnANV3TXFvbJEmaViY7mD8HbAvMBRYDh9f26DFt9uogIg6KiPkRMX/p0qWDqVKSpCGZ1GDOzCWZeXdm3gN8keWrq68FtuqadEvgulH6ODIz52XmvJkzZw62YEmSJtmkBnNEbN5188VAZ4/tU4F9ImKdiNgG2A44fzJrkySpBWsNquOI+DqwK7BJRFwLvBfYNSLmUlZTLwReD5CZCyLiRODXwDLgjZl596BqkySpVX0Fc0Q8OjMvW/WUy2Xmy3s0f3kl0x8CHDKW55Akaarpd1X25yPi/Ij454jYeKAVSZI0jfUVzJn5FOCVlB205kfE8RHxnIFWJknSNNT3zl+Z+TvgPcC/A08HPhURV0TEPw6qOEmSppu+gjkidoqIjwOXA88E9sjMR9XrHx9gfZIkTSv97pV9BOV3x+/KzNs7jZl5XUS8ZyCVSZI0DfUbzLsDt3d+whQR9wPWzcy/ZeZxA6tOkqRppt9tzD8A1uu6/YDaJkmSJlC/wbxuZv61c6Nef8BgSpIkafrqN5hvi4idOzci4nHA7SuZXpIkrYZ+tzG/DTgpIjonltgc2HswJUmSNH31FcyZeUFEPBLYnnKKxisy866BViZJ0jQ0lpNYPB6YUx/z2IggM48dSFWSJE1T/Z7E4jhgW+BioHPWpwQMZkmSJlC/S8zzgB0yMwdZjCRJ012/e2VfBjx0kIVIkqT+l5g3AX4dEecDd3QaM/OFA6lKkqRpqt9gft8gi5AkSUW/P5f6cUTMBrbLzB9ExAOAGYMtTZKk6aff0z4eCHwT+EJtmgWcMqiiJEmarvrd+euNwJOBWwAy83fApoMqSpKk6arfYL4jM+/s3IiItSi/Y5YkSROo32D+cUS8C1gvIp4DnAScNriyJEmanvoN5ncCS4FLgdcDpwPvGVRRkiRNV/3ulX0P8MV6kVawx3OfxuJFVw+7jFHdsGQJMHPYZUhSX/o9VvZV9NimnJkPm/CKtMZZvOhq5h/cbvDNevOiYZcgSX0by7GyO9YF9gIePPHlSJI0vfW1jTkzb+q6LMrMTwDPHHBtkiRNO/2uyt656+b9KEvQGw6kIkmSprF+V2Uf3nV9GbAQeNmEVyNJ0jTX717Zzxh0IZIkqf9V2f+ysvsz82MTU44kSdPbWPbKfjxwar29B3AOcM0gipIkabrqN5g3AXbOzFsBIuJ9wEmZ+bpBFSZJ0nTU7yE5twbu7Lp9JzBnwquRJGma63eJ+Tjg/Ij4FuUIYC8Gjh1YVZIkTVP97pV9SER8D3hqbdo/M385uLIkSZqe+l2VDfAA4JbM/CRwbURsM6CaJEmatvoK5oh4L/DvwMG1aW3gq4MqSpKk6arfJeYXAy8EbgPIzOvwkJySJE24foP5zsxM6qkfI2L9wZUkSdL01W8wnxgRXwA2jogDgR8AXxxcWZIkTU/97pX90Yh4DnALsD3wn5l55kArkyRpGlplMEfEDOD7mflswDCWJGmAVrkqOzPvBv4WEQ+chHokSZrW+j3y19+BSyPiTOqe2QCZ+ZaBVCVJ0jTVbzB/t14kSdIArTSYI2LrzPxjZh4zWQVJkjSdrWob8ymdKxFx8oBrkSRp2ltVMEfX9YcNshBJkrTqYM5RrkuSpAFY1c5fj4mIWyhLzuvV69TbmZkbDbQ6SZKmmZUGc2bOmKxCJEnS2M7HLEmSBmxgwRwRR0XEDRFxWVfbgyPizIj4Xf37oNoeEfGpiLgyIi6JiJ0HVZckSS0b5BLz0cBuI9reCZyVmdsBZ9XbAM8HtquXg4DPDbAuSZKaNbBgzsxzgD+NaN4T6Bys5BjgRV3tx2bxC8rpJTcfVG2SJLVqsrcxb5aZiwHq301r+yzgmq7prq1tkiRNK63s/BU92nr+bjoiDoqI+RExf+nSpQMuS5KkyTXZwbyks4q6/r2htl8LbNU13ZbAdb06yMwjM3NeZs6bOXPmQIuVJGmyTXYwnwrsV6/vB3y7q33funf2LsDNnVXekiRNJ/2e9nHMIuLrwK7AJhFxLfBe4DDgxIg4APgjsFed/HRgd+BK4G/A/oOqS5Kklg0smDPz5aPc9awe0ybwxkHVIknSmqKVnb8kSRIGsyRJTTGYJUlqyMC2MUvSoFx//RLm7Th72GWMavNZszntjHOGXYbWUAazpDVO3rOM+Qe3exyDeYdePewStAZzVbYkSQ0xmCVJaojBLElSQwxmSZIaYjBLktQQg1mSpIYYzJIkNcRgliSpIQazJEkNMZglSWqIwSxJUkMMZkmSGmIwS5LUEINZkqSGGMySJDXEYJYkqSEGsyRJDTGYJUlqiMEsSVJDDGZJkhpiMEuS1BCDWZKkhhjMkiQ1xGCWJKkhBrMkSQ0xmCVJaojBLElSQwxmSZIaYjBLktQQg1mSpIYYzJIkNcRgliSpIQazJEkNMZglSWqIwSxJUkMMZkmSGmIwS5LUEINZkqSGGMySJDXEYJYkqSEGsyRJDTGYJUlqiMEsSVJDDGZJkhpiMEuS1BCDWZKkhhjMkiQ1xGCWJKkhBrMkSQ1ZaxhPGhELgVuBu4FlmTkvIh4MnADMARYCL8vMPw+jPkmShmWYS8zPyMy5mTmv3n4ncFZmbgecVW9LkjSttLQqe0/gmHr9GOBFQ6xFkqShGFYwJ3BGRFwYEQfVts0yczFA/btprwdGxEERMT8i5i9dunSSypUkaXIMZRsz8OTMvC4iNgXOjIgr+n1gZh4JHAkwb968HFSBkiQNw1CWmDPzuvr3BuBbwBOAJRGxOUD9e8MwapMkaZgmPZgjYv2I2LBzHXgucBlwKrBfnWw/4NuTXZskScM2jFXZmwHfiojO8x+fmf8bERcAJ0bEAcAfgb2GUJskSUM16cGcmX8AHtOj/SbgWZNdjyRNtOuvX8K8HWcPu4yeNp81m9POOGfYZWglhrXzlyRNWXnPMuYfPHPYZfQ079Crh12CVqGl3zFLkjTtGcySJDXEYJYkqSEGsyRJDTGYJUlqiMEsSVJDDGZJkhpiMEuS1BCDWZKkhhjMkiQ1xENyriH2eO7TWLyozUPp3bBkCdDm4QclaU1jMK8hFi+6utlj785686JhlyBJU4arsiVJaojBLElSQwxmSZIaYjBLktQQg1mSpIYYzJIkNcRgliSpIQazJEkNMZglSWqIwSxJUkMMZkmSGmIwS5LUEINZkqSGGMySJDXEYJYkqSEGsyRJDTGYJUlqiMEsSVJDDGZJkhpiMEuS1BCDWZKkhhjMkiQ1xGCWJKkhBrMkSQ0xmCVJaojBLElSQwxmSZIaYjBLktQQg1mSpIYYzJIkNcRgliSpIWsNu4CW7PHcp7F40dXDLqOnG5YsAWYOuwxJ0oAZzF0WL7qa+Qe3GX6z3rxo2CVIkiaBq7IlSWqIS8ySNI1cf/0S5u04e9hl9LT5rNmcdsY5wy5j6AxmSZpG8p5lzW6ym3dom/v4TDZXZUuS1BCDWZKkhhjMkiQ1xG3MkqQmtLxjGkzezmnNBXNE7AZ8EpgBfCkzDxtySZKkSdDyjmkweTunNbUqOyJmAJ8Bng/sALw8InYYblWSJE2epoIZeAJwZWb+ITPvBL4B7DnkmiRJmjStBfMs4Jqu29fWNkmSpoXIzGHXcK+I2At4Xma+rt5+NfCEzHxz1zQHAQfVm9sDv5n0QnvbBLhx2EUMmGOcOqbDOB3j1DEVxzk7M3tuUG9t569rga26bm8JXNc9QWYeCRw5mUX1IyLmZ+a8YdcxSI5x6pgO43SMU8d0GWdHa6uyLwC2i4htIuL+wD7AqUOuSZKkSdPUEnNmLouINwHfp/xc6qjMXDDksiRJmjRNBTNAZp4OnD7sOlZDc6vXB8AxTh3TYZyOceqYLuMEGtv5S5Kk6a61bcySJE1rBnMPEbFbRPwmIq6MiHf2uP9pEXFRRCyLiJd2tc+NiJ9HxIKIuCQi9u66b5uIOC8ifhcRJ9Sd24ZmQGM8OiKuioiL62XuZI1nNOMY5+yIuLCOY0FEvKHrvsdFxKW1z09FREzWeHoZ0BjPrn125uWmkzWeXlZ3jF33bxQRiyLiiK62puZjrWkQ45wy8zIi7u4ax6ld7U19vo5bZnrpulB2Ovs98DDg/sCvgB1GTDMH2Ak4FnhpV/sjgO3q9S2AxcDG9faJwD71+ueBf5qCYzy6e9phX8Y5zvsD69TrGwALgS3q7fOBJwIBfA94/hQc49nAvGHPw/GOsev+TwLHA0d0tTUzHwc8zikzL4G/jtJvM5+vE3Fxifm+VnlY0MxcmJmXAPeMaP9tZv6uXr8OuAGYWb+JPxP4Zp30GOBFgx3GSk34GCen7DEbzzjvzMw76s11qGuXImJzYKPM/HmWT4FjWXPnZc8xNmi1xwhlyRjYDDijq621+QgDGGeDxjXGXhr8fB23Vv8Rh2lCDgsaEU+gfCP8PfAQ4C+ZuWw8fU6gQYyx45C6ivvjEbHO+Moct3GNMyK2iohLah8frl9EZtV+VqvPARjEGDu+UlcZ/seQV/Ou9hgj4n7A4cA7evTZ0nyEwYyzY42fl9W6ETE/In4REZ3wbe3zddwM5vvq9aYd067r9dv4ccD+mXnPRPQ5wQYxRoCDgUcCjwceDPz7eIqcAOMaZ2Zek5k7AQ8H9ouIzcbb5wAMYowAr8zMfwCeWi+vHnelq288Y/xn4PTMvGZEe2vzEQYzTpg68xJg6yxHAHsF8ImI2HYC+myOwXxfqzws6MpExEbAd4H3ZOYvavONwMYR0fnd+Jj6HIBBjJHMXJzFHcBXKKuthmlc4+yoS5ELKB9q19Z+xtXnBBrEGMnMRfXvrZRtlsOcl+MZ4xOBN0XEQuCjwL4RcRjtzUcYzDin0rzsvE/JzD9Qtp0/lvY+X8fNYL6v1T4saJ3+W8CxmXlSp71uw/oR0NnDcD/g2xNa9dhM+BjrfZvXv0HZxnPZhFY9duMZ55YRsV69/iDgycBvMnMxcECGnQQAAAiJSURBVGtE7FLHuS9r7rzsOcaIWCsiNqntawMvYLjzcrXHmJmvzMytM3MO8K+U9+07G5yPMIBxTqV5GREP6mweq2N6MvDrBj9fx2/Ye5+1eAF2B35L2Xb67tr2AeCF9frjKd/8bgNuAhbU9lcBdwEXd13m1vseRtkL9ErgJOresFNsjD8ELqX8438V2GANnpfPAS6h7DV6CXBQV5/z6hh/DxxBPVDPVBkjsD5wYW1bQNnTd8aaOMYRfbyGFfdWbmo+DmKcU2leAk+qny+/qn8P6Oqzqc/X8V488pckSQ1xVbYkSQ0xmCVJaojBLElSQwxmSZIaYjBLktQQg1lD1XW2mF/VM8o8aYL7f9dE9jei79dExNJa/xUR8f/G0deuEfGdev2Fvc6600cfcyNi967bq9XPeEXEU6Ocrerizu+kR9z/4ojIiHjkZNc2VnW+3BwRv4yIyyPivQN+vtMjYuN6+edBPpfaZTBr2G7PzLmZ+RjKIT0PneD+ewZzFBPx/j8hM+dSDnbw7ojYalUPWJXMPDUzD1uNh86l/EZ0vP2M1yuBj9b5enuP+18O/JRycIlxi4gZE9HPSvwkMx9L+d3zq+rJIlap60hUfcvM3TPzL8DGlMNsahoymNWSjYA/w73B+d8RcVmUc+buvYr2zSPinLqUdlldajsMWK+2fS0i5tSlns8CFwFbRcTnohwUf0FEvL9TSEQsjIgPR8T59fLwlRWemTdRDm7QOfrZzIg4OSIuqJcn1/YnRMS5dQns3IjYfmRfdUn8iHr94q7L7RHx9F591KMofQDYu06794h+ZkfEWVFOMHJWRGxd24+Oci7icyPiD1HPf9vr9exR57NqDZdGxFERsU5EvA54GfCfEfG1Ho/ZgPIl5gC6gjnKOXS7l/aPjoiXRMSMOr8vqLW/vt6/a0T8KCKOpxxsgog4Jcr5pRdExEFdfR0QEb+Ncl7iL3a9Jj3n0Urm8W2Ug3VsGxHrRsRX6th/GRHP6Jp3J0XEacAZY3m/1vaFUY5qdVh9novr44+LiHvPwlTfzy9cWb1agw37CCdepvcFuJty9LArgJuBx9X2lwBnUs7fuhnwR0rojdb+dpYfRWgGsGG9/teu55pDOZXcLl1tD+56zNnATvX2wq7+9gW+06P217D8CEtb13GsW28fDzyl677L6/WNgLXq9WcDJ9fru3aegxFHqKptewA/AdZeSR8rPG5EfacB+9XrrwVOqdePphwp6X7ADpRT8jHa69nV97qUswQ9ot4+FnhbV589z8tNOXLcl+v1c4Gd6/UXA8fU6/evfa8HHEQ5JjuUU1POB7apr9dtwDY95uV6lCN6PYRyzvCFlJOqrF1fw85r0nMejai3e748pPa1Y319vlLbH0l5H65bX/Nru2oZ6/t1IbAJ5b16WVcdT++aZw8Eruq8B7xMvcuYV7VIE+z2LKuCiYgnAsdGxKOBpwBfz8y7gSUR8WPKofpGa78AOCrK8YBPycyLR3m+q7PrxBvAy+rS1VqUD8wdKIcvBPh619+Pj9Lf3nVpaXvgwMz8e21/NrBDLD/D3kYRsSHlQ/WYiNiOcgactVf1AtVp/xt4ZmbeFREPHWsflJMc/GO9fhzwka77TslyhrBfx/KzS63q9dweuCozf1tvHwO8EfjEKup4edc036i3LwK+B3wqyrGQdwPOyczbI+K5wE6dJXnK67cdcCdwfmZe1dX3WyLixfX6VnW6hwI/zsw/AUTEScAj6jQ951GWkz10e2pE/JLype6wzFwQER8CPg2QmVdExNVd/Z7ZeT7G/36lPsePI+IzEbEpZT6enMtPc6gpxlXZakZm/pyytDCT3qdyY7T2zDwHeBqwCDguIvYd5fG33dtRxDaUA/4/K8upD79LWeq5t9tRrnc7ITN3pJyV6fAamlD+t56YZTvr3MycVT/wPwj8KDMfTVkKXrd3t/fWuD5wIiX0O2fMGVMfo+gezx3dTwl9vZ5jPqdvRDyEckL7L0U5C9I7KF9son6hORt4HrA3JbQ7z/Pmrtdxm8w8o97XPS93pQTtE7Psr/BLyuuysjpHm0cj/SQzH5uZj8vMz/cx/tu6ro/3/drtOMr2+/0pZ2/TFGUwqxlR9tKdQTlw/TmUD+0ZETGT8iF2/mjtETEbuCEzvwh8Gdi5dntXXSrpZSPKh+jNdUnx+SPu37vr789XVnv9UnEc8NbadAbwpq6xza1XH0j5MIay2nNVvkJZZfqTrrbR+rgV2HCUfs5l+TbdV1J2vhrVSl7PjiuAObF82/urgR+vrE/K2X+OzczZmTknM7eirJJ9Sr3/G5TQeSrw/dr2feCfOvMwIh5Rv6yM9EDgz5n5t/o+2qW2nw88PcqZidairFruGG0e9eMcyutIRDyCsir8N6NMN5b3a0eveXk08DaAzFwwhlq1hjGYNWydnbMuBk6gbAe9m3Jqyc6Zj34I/FtmXr+S9l2Bi+sqx5dQzqIDcCRwSa8dkTLzV5QlqwXAUcDPRkyyTkScRwnbfn4K9WFg/7rK+i3AvLrD0q+BN9RpPgIcGhE/o3wJGVX98H4p8NpYvgPYvJX08SPKqtmLOzsZdXlLre0SSoi+lZXbld6vJwB1CXd/4KSIuJSymvfzIzsZ4eWU+dftZMpJ76EE5dOAH2TmnbXtS8CvgYsi4jLgC9BzE9z/AmvV8X0Q+EWtcxHwX8B5wA9qXzfXx4w2j/rxWWBGHfsJwGuynId8pLG+X6l13wT8rO4Y9t+1bQlwOS4tT3meXUrqoa5qnZeZNw67Fo1PRGyQmX+tS8zfAo7KzJFfEJoXEQ+g7IG+c2bevKrpteZyiVnSVPe+ukbmMsqq81OGXM+YRcSzKZsPPm0oT30uMUuS1BCXmCVJaojBLElSQwxmSZIaYjBLktQQg1mSpIYYzJIkNeT/AxUfqFRPHCLJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import random                          # import random package\n",
    "L = 1000                               # set the number of bootstrap realizations   \n",
    "por_avg_real = []                      # declare an empty list to store the bootstrap realizations of the statistic \n",
    "for k in range(0,L):                   # loop over the L bootstrap realizations\n",
    "    samples = random.choices(df['Porosity'].values, k=len(df)) # n Monte Carlo simulations\n",
    "    por_avg_real.append(np.average(samples)) # calculate the statistic of interest from the new bootstrap dataset\n",
    "plt.hist(por_avg_real,color = 'darkorange',alpha = 0.8,edgecolor = 'black') # plot the distribution, could also calculate any summary statistics\n",
    "plt.xlabel('Boostrap Realizations of Average Porosity'); plt.ylabel('Frequency'); plt.title('Uncertainty Distribution for Average Porosity')\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.1, wspace=0.2, hspace=0.2); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we proceed with a more complicated, robust work by first quantifying the spatial bias in the samples and assigned data weights to mitigate this bias.\n",
    "\n",
    "#### Declustering\n",
    "\n",
    "Let's calculate some declustering weights. There is a demonstration on declustering here https://git.io/fhgJl if you need more information. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 58 data with:\n",
      "   mean of      0.133789477 \n",
      "   min and max  0.065626061 and 0.22140153\n",
      "   standard dev 0.039179360550857516 \n",
      "Declustered mean = 0.121 and declustered standard deviation = 0.033\n"
     ]
    }
   ],
   "source": [
    "wts, cell_sizes, dmeans = geostats.declus(df,'X','Y','Porosity',iminmax = 1, noff= 10, ncell=100,cmin=10,cmax=2000)\n",
    "df['Wts'] = wts                            # add weights to the sample data DataFrame\n",
    "df.head()                                  # preview to check the sample data DataFrame\n",
    "\n",
    "def weighted_avg_and_std(values, weights): # function to calculate weighted mean and st. dev., from Eric O Lebigot, stack overflow,\n",
    "    average = np.average(values, weights=weights)\n",
    "    variance = np.average((values-average)**2, weights=weights)\n",
    "    return (average, math.sqrt(variance))\n",
    "\n",
    "sample_avg, sample_stdev = weighted_avg_and_std(df['Porosity'],df['Wts'])\n",
    "print('Declustered mean = ' + str(round(sample_avg,3)) + ' and declustered standard deviation = ' + str(round(sample_stdev,3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### A Couple of Bootstrap Realizations\n",
    "\n",
    "We will attempt boostrap by-hand and manually loop over $L$ realizations and draw $n$ samples to calculate the summary statistics of interest, mean and variance. The choice function from the random package simplifies sampling with replacement from a set of samples with weights.\n",
    "\n",
    "This command returns a ndarray with k samples with replacment from the 'Porosity' column of our DataFrame (df) accounting for the data weights in column 'Wts'.\n",
    "```p\n",
    "samples1 = random.choices(df['Porosity'].values, weights=df['Wts'].values, cum_weights=None, k=len(df))\n",
    "```\n",
    "\n",
    "It is instructive to look at a couple of these realizations from the original declustered data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bootstrap means, realization 1 = 0.12224859570689653 and realization 2 = 0.1190454733448276\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABU4AAAF6CAYAAADGYcJQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde7xld1kf/s9DBoSEEEAGZCYMAYRoRBvomB+KF+RSQa6+igoFRMXG6k8gSotgTcHGW1vKpfXXYuR+EdSIiIqVCEakBWIIQSDBBiGSORPIIMQwIVxCnt8faw2c7Jwzc2bOPmfvfc77/Xrt19nrsr/rWd+9Zj+znnWr7g4AAAAAAF9zi1kHAAAAAAAwbxROAQAAAAAmKJwCAAAAAExQOAUAAAAAmKBwCgAAAAAwQeEUAAAAAGCCwilTUVUfrqoHzTqORVVVD6qqfbOO42hV1SlV1VW1Y9axrMa2CcwDv0VstKr6s6p66qzjADgW8iQbTZ7kWCmcckRVdUVVPXRi3I9V1bsODXf3t3T3BUdoZ+6LbIczrvNXqurg+Pp4Vb2yqu4z69gOWem7mqUxnuvH/vrU2F+33cwYlm+bVfX8qnrdZi4f2PrkycEKefJjVfXTU2r7VVX1K2ucd9N/66vqgqr6wrjen66qN1XVXTczhu5+RHe/eoznJtvf0aqq76uqv6yqf6qqK6YWJLAtyZMDeXJL5cl/V1UfqqrPjXWBfze9SJk3CqdsGZuUQN/d3bdNclKShya5Psn7quq+m7DsDVWDjfhNePTYZ/dP8u1JfuloG1jU/xwBzJPNzJPj7/7jk/znqrrfJix3zTYw3/3suN73SXL7JC862gbmKN9dl+QVSewIAtuGPDmQJ9ekkvxokjskeXiSn62qJ8w2JDaKwilTsfwoYlWdUVUXVdW141mGLxxne+f495rxSNN3VNUtquqXquofqurqqnpNVZ20rN0fHaf9Y1WdPbGc51fVeVX1uqq6NsmPjct+d1VdU1VXVdVvVtWtlrXXVfUzVXX5eHTonKq61/iZa6vq95bPv5ru/kp3/313/0ySv0ry/GXLeEBV/Z8xhg/UsktOquqONZx1ub+qPltVb16lP7uqvnHZ8FePIFbVnarqT8b2P1NVfz3242uT7Enyx2P/PnsN8VxQVb9aVf87yeeT3LOqTqqql4/9t1RVv1JVx43zH1dVLxiPEn4sySOP1FfL+mwpyZ8lue/Y1q6qesu4Dh+tqn+9LK6Vvtuvq6oXj323f3z/dYfrk3HaFVX10Kp6eJJfTPIjY/98oKp+qKreN9H3z1rtewE4VtstTyZJd1+c5LIk37ys/cfUcDnmNWMOWj7tm8dx14zzPGYcf2aSJyV59tgvfzyO/4UxT32uqv6uqh6y0m/9OO9K+e7Hq+qy8fMfq6qfWhbLg6pqX1X94pjzrqiqJ61xvT+T5A/ytXx30vi9HRi/q19alqN+rKr+d1W9qKo+k+T5h/vOq+rW4/f5j2M//U1V3WXZOv7k2KcvTfIdYx9cU1XfPm5rX93hrKp/WVWXrLIOF3b3a5N8bC3rDLBe8uRX25cnFyNP/ufuvri7b+juv0vyR0keuJb1ZwF1t5fXYV9Jrkjy0IlxP5bkXSvNk+TdSZ4yvr9tkgeM709J0kl2LPvcTyT5aJJ7jvO+Kclrx2mnJTmY5LuS3CrJC5J8edlynj8OPy7DQYDbJPnnSR6QZMe4vMuSnLVseZ3kLUlul+RbknwxydvH5Z+U5NIkT12lH26yzhPr8Knx/e4k/5jkB8aYHjYO7xyn/2mS381wZOqWSb53HP+gJPsm4vzGZcOvSvIr4/tfz/BDf8vx9d1JaqXvag3xXJDkE2Nf7Bjbe3OS30pyQpI7J7kwyU+N8/+bJB9Jcrckd0zyl5Pf6WrbzviZDyc5Zxz+qyT/I8mtk5ye5ECShxzmu/2PSd4zxrQzyf9Z1taa+mRs93XL4vu6JJ9J8s3Lxr0/yb+c9b87Ly+vxXlN/vaO426SM7IN82SGqwyuSXKfcfg+Gc5kfNj4W/3scd1uNQ5/NMPO3K2SPDjJ55KcOn72VRnz4Dh8apIrk+xa1nf3Wrber5uI7YLcPN89Msm9Mpw18r0ZdhTvP87/oCQ3JHnhmCu+d4z91FXW/YIkPzm+v1OSdyz7nl6TYYfqxDHO/5vkacv67IYkTx/jus0RvvOfSvLHSY5Pctz4fd5uhRhu8l2M4y5N8ohlw3+Y5FlH2LYfmuSKWf8b8/LyWuxX5MnV1lme7MXOk+N8lWEf8t/M+t+a18a8nHHKWr15PBJzTVVdk6HgtZovJ/nGqrpTdx/s7vccZt4nJXlhd3+suw8meW6SJ4xHeh6f5I+7+13d/aUk/yFDolru3d395u6+sbuv7+73dfd7ejjyc0WGAuD3TnzmP3X3td394SQfSvK2cfn/lOGMyKO9VGJ/hiJikjw5yVu7+61jTOcnuSjJD9RwD5dHZPhB/Wx3f7m7/+ool5UM/XvXJHcf2/jr7p7sl0NWjWfZPK/q7g939w3jejwiw38OruvuqzNcQnHosoMfTvLi7r6yhyOFv76GeN88bjPvylAs/bWquluG/8D8Qnd/obsvSfKyJE9Z9rmbfLcZtpX/2N1Xd/eBJL+8bP6j6ZOv6u4vZihkPzlJqupbMiTrP1nDegEsJ08OHjD2wcEMB95em+TycdqPJPnT7j6/u7+cYQf2Nkm+M8NO6m2T/EZ3f6m735Hht/iJqyznKxl21E6rqlt29xXd/feHiStZlu/GXPGnPVw90mM+fluGA2/Lnd3dXxyn/2mGPLia/zZ+9x9IclWSn6/hio0fSfLc7v7c2Of/NTfNd/u7+7+PcR3Kd6t9519O8vUZDq5+Zfw+rz3Ceh/y6nwt390xyfcn+Z01fhZgveTJgTy59fLk8zMU3l+5xuWwYBROWavHdfftD72S/Mxh5n1ahqNlHxlPjX/UYebdleQflg3/Q4YjSXcZp115aEJ3fz7D2ZLLXbl8oKruU8Ml258cL7f4tQxHtJb71LL3168wfLQPL9qd4azFJLl7kh+a+E/Bd2Uo6t0tyWe6+7NH2f6k/5LhCNvbxksmnnOYeQ8XzyFXTsx/yyRXLZv/tzKc5ZlMfCe56Xe3mkPbzt27+2fGZLcrQ198bqKt3avEdWjZk9vKrvH90fTJpFcn+VdVVRkS9O+NBVWAoyFPDt4z9sFtk3xDhjNXfm2ldenuG8f4dh9al3Hc8nVdnhey7LMfTXJWhp2Vq6vqjVW1a6V5l5nsi0dU1XtquMXLNRkOKi7vi89293UT8RxuGc8Y1313dz9pPMh3pwxnBk1+h0eb7w59569N8udJ3ljDbWv+c1Xd8jAxLfe6JI+u4SGNP5zkr7v7qjV+FmC95MmBPLmF8mRV/WyGe50+0j7k1qVwytR19+Xd/cQMxbb/lOS8qjohNz+6lwxna9592fCeDKfifyrDUaiTD02oqttkOHp0k8VNDP/PDJeS37u7b5fhUoY69rVZkx9M8tfj+yszXCZw+2WvE7r7N8Zpd6yq26+hzc9nuLzgkG849GY8Eves7r5nkkdnOFL3kEOTJ9o5XDxZ4TNXZrjc5E7L5r9dd3/LOP2qDAXgQ/asYV1Wsj9DX5w40dbSKnEd+szktrI/OWKfLHezbXA8gv2lDEdP/1WGZAuwYbZLnuzuT2W4h9mjx1E3WZfxgNXdMvz2709yt7rpwyiW54WVfr9/p7u/a2yzM/TlivNOjq/hHtl/kOFsnruMO/FvzU374g7j97I8nv2rre8qPp3h7JfJ7/Bo890NGW4L9OXu/uXuPi3DGUiPyrDDNmml/lrKcPnrD2Y4UCjfAXNJnvxqvPLkyvHORZ6sqp9I8pwMt5vbd7h5WWwKp0xdVT25qnaOR8OuGUd/JcM9LG/McC+SQ96Q5Oeq6h7jkZ1fS/K7PVw2fl6GIz7fWcMNtn85R05aJya5NsnBqvqmJD89tRVbpoaHJN2jqv57hvu7/PI46dBRqu8f57l1DTfOPnk8WvVnSf5HVd2hqm5ZVd+zyiIuyXAW5HE13MD7q5eHVNWjquobx0R6bYa+/co4+VO5af+uGs9KCx1jfFuS/1pVt6vhxtv3qqpDy/+9JM+oqpOr6g4ZEsVR6+4rM9yj9NfHmL4tw5Hl1x/mY29I8ktVtbOq7pThUpvXJUfsk+U+leSUuvlTIl+T5DeT3NDd7zqWdQJYq+2QJ5Okqr4+w87Hh8dRv5fkkTU8nOKWSZ6V4WDd/0ny3gz3Rnv2mB8flGFH8o3jZ2+S36rq1Kp68Lhj94UMZ/gsz4Ur/dYvd6sMlzAeSHJDVT0iyb9YYb5frqpbVdV3Z9j5+v2j6YPu/sq43r9aVSdW1d2T/HzG/LWKVb/zqvq+qvrWGi5tvDbDzuZq+e7kuvkDSl6T4Z5535rh3m0rGvP/rTNchVJjrl7Tw04A1kuelCcz/3nySeNyH9bdHqS4xSmcshEenuTDNdy35SVJntDDfSw/n+RXk/zvGi4Df0CSV2Q4kvPOJB/P8KP+9CTp4Z4xT8+QDK7KcPPrqzMkj9X82wxnDX4uyW9nuH/lNH3HuF7XZri59O2SfHt3f3CM+cokj81wZPJAhjM4/12+9m/tKRl+vD8yrstZqyznmRkS4TUZ7uGy/Cnv907yFxludP7uJP+juy8Yp/16huLiNVX1b9cQz0p+NEOivDTJZzP8h+PQpf2/neHShw8kuTjDjbiP1RMz3E90f4ak9Lwe7sG6ml/JcH/Wv03ywXH5vzJOO1yfLHcokf9jVV28bPxrMzzV0dk3wGbY8nlyXLfLMuSeQ/H+XYZ7h/33DGeYPDrJo3u4V9uXkjwmw322P53h3nc/2t0fGdt9eYb7tF1TVW/OsDP3G+O8n8xwVtIvjvOu9lv/VeOtYp6RYWftsxn65C0Ts31ynLY/w4G9f7MsnqPx9Aw7ux/LcL/v38nwva5m1e88wxUo52X4f8hlGe4dvtLO5Tsy7Ih/sqo+vWz8H2Y4S+cPJy6vnPQ9GXay35rhTJ7rMxxYBdgM8qQ8Oe958lcynL38N4e+z6p66eFWksV16KnTMPfGo0nXZLhs4uOzjoeto4bLdq7O8JTIy480P8A8kienZzyT53XdveIVGousqv4+yU9191/MOhaAzSRPTo88yXbijFPmWlU9uqqOr+HeKS/IcKbhFbONii3op5P8jaIpsGjkSY5GVf3LDPd1e8esYwHYDPIkR0OeZCU7Zh0AHMFjM5yGXxku035CO02aKaqqKzJsX4+bcSgAx0KeZE2q6oIkpyV5ysRTmQG2MnmSNZEnWY1L9QEAAAAAJrhUHwAAAABggsIpAAAAAMCEhbjH6Z3udKc+5ZRTZh0GAJvgfe9736e7e+es41gUciTA9iFHHj15EmD72Ig8uRCF01NOOSUXXXTRrMMAYBNU1T/MOoZFIkcCbB9y5NGTJwG2j43Iky7VBwAAAACYoHAKAAAAADBB4RQAAAAAYILCKQAAAADABIVTAAAAAIAJCqcAAAAAABMUTgEAAAAAJiicAgAAAABMUDgFAAAAAJigcAoAAMDcqKpXVNXVVfWhZePuWFXnV9Xl4987zDJGALYHhVMAAADmyauSPHxi3HOSvL27753k7eMwAGwohVMAAADmRne/M8lnJkY/Nsmrx/evTvK4TQ0KgG1J4RQAAIB5d5fuvipJxr93nnE8AGwDO2YdAPPtjNNPy76lpam0dfLu3bnwkkun0hYAzAN5EmD+VNWZSc5Mkj179sw4msUyzbz25S99Ibe81a2n0pYcCcyKwimHtW9pKfvPOWEqbe06ezoJGADmhTwJsGk+VVV37e6rququSa5ebcbuPjfJuUmyd+/e3qwAt4Jp5rUTz7o2B/7T10+lLTkSmBWX6gMAADDv3pLkqeP7pyb5oxnGAsA2oXAKAADA3KiqNyR5d5JTq2pfVT0tyW8keVhVXZ7kYeMwAGwol+oDAAAwN7r7iatMesimBgLAtueMUwAAAACACQqnAAAAAAATFE4BAAAAACYonAIAAAAATFA4BQAAAACYoHAKAAAAADBB4RQAAAAAYILCKQAAAADABIVTAAAAAIAJCqcAAAAAABMUTgEAAAAAJiicAsAMVNUrqurqqvrQCtP+bVV1Vd1pFrEBAACgcAoAs/KqJA+fHFlVd0vysCSf2OyAAAAA+BqFUwCYge5+Z5LPrDDpRUmenaQ3NyIAAACWUzgFgDlRVY9JstTdH5h1LAAAANvdjlkHAAAkVXV8kn+f5F+sYd4zk5yZJHv27NngyAAAALYnZ5wCwHy4V5J7JPlAVV2R5OQkF1fVN0zO2N3ndvfe7t67c+fOTQ4TAABge3DGKQDMge7+YJI7Hxoei6d7u/vTMwsKAABgG3PGKQDMQFW9Icm7k5xaVfuq6mmzjgkAAICvccYpAMxAdz/xCNNP2aRQAAAAWIEzTgEAAAAAJiicAgAAAABMUDgFAAAAAJigcAoAAAAAMEHhFAAAAABgwo5ZBwAAAACwnZ1x+mnZt7Q0lbZO3r07F15y6VTagu1O4RQAAABghvYtLWX/OSdMpa1dZ0+nAAu4VB8AAAAA4GY2rHBaVa+oqqur6kPLxt2xqs6vqsvHv3fYqOUDAAAAAByrjTzj9FVJHj4x7jlJ3t7d907y9nEYAAAAAGCubFjhtLvfmeQzE6Mfm+TV4/tXJ3ncRi0fAAAAAOBYbfY9Tu/S3Vclyfj3zpu8fAAAAACAI5rbh0NV1ZlVdVFVXXTgwIFZhwMAAAAAbCObXTj9VFXdNUnGv1evNmN3n9vde7t7786dOzctQAAAAACAzS6cviXJU8f3T03yR5u8fAAAAACAI9qwwmlVvSHJu5OcWlX7quppSX4jycOq6vIkDxuHAQAAAADmyo6Nari7n7jKpIds1DIBAAAAAKZhbh8OBQAAAAAwKwqnAAAAAAATFE4BAAAAACYonAIAAAAATFA4BQAAAACYoHAKAAAAADBB4RQAAAAAYILCKQAAAADABIVTAAAAAIAJCqcAAAAAABMUTgEAAAAAJiicAgAAAABMUDgFAAAAAJigcAoAAAAAMEHhFAAAAABggsIpAMxAVb2iqq6uqg8tG/dfquojVfW3VfWHVXX7WcYIAPOmqn6uqj5cVR+qqjdU1a1nHRMAW5fCKQDMxquSPHxi3PlJ7tvd35bk/yZ57mYHBQDzqqp2J3lGkr3dfd8kxyV5wmyjAmArUzgFgBno7ncm+czEuLd19w3j4HuSnLzpgQHAfNuR5DZVtSPJ8Un2zzgeALawHbMOAABY0U8k+d2VJlTVmUnOTJI9e/ZsZkwAMDPdvVRVL0jyiSTXJ3lbd79tcj55cuu5/rqD2bXzpKm0dfLu3bnwkkun0haw9SmcAsCcqap/n+SGJK9faXp3n5vk3CTZu3dvb2JoADAzVXWHJI9Nco8k1yT5/ap6cne/bvl88uTWc8ONN2b/OSdMpa1dZy9NpR1ge3CpPgDMkap6apJHJXlSd9vZA4CveWiSj3f3ge7+cpI3JfnOGccEwBamcAoAc6KqHp7kF5I8prs/P+t4AGDOfCLJA6rq+KqqJA9JctmMYwJgC1M4BYAZqKo3JHl3klOral9VPS3JbyY5Mcn5VXVJVb10pkECwBzp7vcmOS/JxUk+mGF/9tyZBgXAluYepwAwA939xBVGv3zTAwGABdLdz0vyvFnHAcD24IxTAAAAAIAJCqcAAAAAABMUTgEAAAAAJiicAgAAAABMUDgFAAAAAJigcAoAAAAAMEHhFAAAAABggsIpAAAAAMAEhVMAAAAAgAkKpwAAAAAAExROAQAAAAAmKJwCAAAAAExQOAUAAAAAmKBwCgAAAAAwQeEUAAAAAGCCwikAAAAAwASFUwAAAACACQqnAAAAAAATFE4BAAAAACYonAIAAAAATJhJ4bSqfq6qPlxVH6qqN1TVrWcRBwAAAADASja9cFpVu5M8I8ne7r5vkuOSPGGz4wAAAAAAWM2sLtXfkeQ2VbUjyfFJ9s8oDgAAAACAm9n0wml3LyV5QZJPJLkqyT9199s2Ow4AAAAAgNXM4lL9OyR5bJJ7JNmV5ISqevIK851ZVRdV1UUHDhzY7DABAAAAgG1sFpfqPzTJx7v7QHd/Ocmbknzn5EzdfW537+3uvTt37tz0IAEAAACA7WsWhdNPJHlAVR1fVZXkIUkum0EcAAAAAAArmsU9Tt+b5LwkFyf54BjDuZsdBwAAAADAanbMYqHd/bwkz5vFsgEAAAAAjmQWl+oDAAAAAMw1hVMAAAAAgAkKpwAwA1X1iqq6uqo+tGzcHavq/Kq6fPx7h1nGCAAAsJ0pnALAbLwqycMnxj0nydu7+95J3j4OAwAAMAMKpwAwA939ziSfmRj92CSvHt+/OsnjNjUoAAAAvkrhFADmx126+6okGf/eecbxAAAAbFsKpwCwYKrqzKq6qKouOnDgwKzDAQAA2JIUTgFgfnyqqu6aJOPfq1eaqbvP7e693b13586dmxogAADAdqFwCgDz4y1Jnjq+f2qSP5phLAAAANuawikAzEBVvSHJu5OcWlX7quppSX4jycOq6vIkDxuHAQAAmIEdsw4AALaj7n7iKpMesqmBAAAAsCJnnAIAAAAATHDGKdveGaefln1LS1Np6+Tdu3PhJZdOpS0AAAAAZkfhlIU0zWLndQcP5p9edJeptLXr7OnEBAAAAMBsKZxuQdMuKiYnTKWtadq3tJT950wnrhPPunYq7QAAAACwdSicbkGKigAAAACwPh4OBQAAAAAwQeEUAAAAAGCCwikAAAAAwASFUwAAAACACQqnAAAAAAATFE4BAAAAACYonAIAALAQqur2VXVeVX2kqi6rqu+YdUwAbF07Zh0AAAAArNFLkvyv7n58Vd0qyfGzDgiArWtNZ5xW1X03OhAAWERyJACsbpp5sqpul+R7krw8Sbr7S919zbTaB4BJa71U/6VVdWFV/UxV3X5DIwKAxSJHAsDqppkn75nkQJJXVtX7q+plVXXCFGIEgBWt6VL97v6uqrp3kp9IclFVXZjkld19/oZGBwBzTo5cLGecflr2LS1Nrb3rDh5MMn/77NNcz5N3786Fl1w6lbaA7WfKeXJHkvsneXp3v7eqXpLkOUnOXj5TVZ2Z5Mwk2bNnz7rih8OZZr6d1/9TwHa35nucdvflVfVLSS5K8t+S3K+qKskvdvebNipAAJh3cuTi2Le0lP3nTG+n5MSzrp1aW9M0zfXcdfb0Cs3A9jTFPLkvyb7ufu84fF6Gwunk8s5Ncm6S7N27t9cVPBzGNPPtvP6fAra7td7j9Nuq6kVJLkvy4CSP7u5vHt+/aAPjA4C5JkcCwOqmmSe7+5NJrqyqU8dRD0nilHgANsxazzj9zSS/neGI4PWHRnb3/vHIIQBsV3IkAKxu2nny6UleX1W3SvKxJD8+nTAB4ObWWjj9gSTXd/dXkqSqbpHk1t39+e5+7YZFBwDzT44EgNVNNU929yVJ9k45RgBY0Zou1U/yF0lus2z4+HEcAGx3ciQArE6eBGBhrbVweuvuPnhoYHx//MaEBAALRY4EgNXJkwAsrLUWTq+rqvsfGqiqf57k+sPMDwDbhRwJAKuTJwFYWGu9x+lZSX6/qvaPw3dN8iMbExIALBQ5EgBWJ08CsLDWVDjt7r+pqm9KcmqSSvKR7v7yhkYGAAtAjgSA1cmTACyytZ5xmiTfnuSU8TP3q6p092s2JCoAWCxyJACsTp4EYCGtqXBaVa9Ncq8klyT5yji6k0h2AGxrciQArE6eBGCRrfWM071JTuvu3shgAGAByZEAsDp5EoCFdYs1zvehJN+wkYEAwIKSIwFgdfIkAAtrrWec3inJpVV1YZIvHhrZ3Y/ZkKgAYHHIkQCwOnkSgIW11sLp8zcyCABYYM+fdQAAMMeeP+sAAOBYralw2t1/VVV3T3Lv7v6Lqjo+yXEbGxoAzL+NyJFV9XNJfjLDwzM+mOTHu/sL648WADaXfUkAFtma7nFaVf86yXlJfmsctTvJmzcqKABYFNPOkVW1O8kzkuzt7vtm2Ll8wnrjBIBZsC8JwCJb68Oh/t8kD0xybZJ09+VJ7nysC62q21fVeVX1kaq6rKq+41jbAoAZm2qOHO1Icpuq2pHk+CT719keAMzKRuRJANgUay2cfrG7v3RoYNyR63Us9yVJ/ld3f1OSf5bksnW0BQCzNNUc2d1LSV6Q5BNJrkryT939tnVHCQCzMe19SQDYNGstnP5VVf1ihrNfHpbk95P88bEssKpul+R7krw8Sbr7S919zbG0BQBzYGo5Mkmq6g5JHpvkHkl2JTmhqp48Mc+ZVXVRVV104MCBdYQOABtuqnkSADbTWgunz0lyIMMDKn4qyVuT/NIxLvOeY1uvrKr3V9XLquqEY2wLAGZtmjkySR6a5OPdfaC7v5zkTUm+c/kM3X1ud+/t7r07d+5cx6IAYMNNO08CwKbZsZaZuvvGJL89vqaxzPsneXp3v7eqXpIhmZ69fKaqOjPJmUmyZ8+eKSwWAKZvyjkyGS7Rf8D41OHrkzwkyUVTahsANtUG5EkA2DRrKpxW1cezwn1ouvuex7DMfUn2dfd7x+HzMhROJ9s+N8m5SbJ37173wAFgLk05R2Y8qHhekouT3JDk/RnzIQAsmmnnSQDYTGsqnCbZu+z9rZP8UJI7HssCu/uTVXVlVZ3a3X+X4UyaS4+lLQCYA1PLkYd09/OSPG89bQDAnJh6ngSAzbKme5x29z8uey1194uTPHgdy316ktdX1d8mOT3Jr62jLQCYmQ3IkQCwZciTACyytV6qf/9lg7fIcNTwxGNdaHdfkpseeQSAhTTtHAkAW4k8CcAiW+ul+v912fsbklyR5IenHg0ALB45EgBWJ08CsLDWVDjt7u/b6EAAYBHJkQCwOnkSgEW21kv1f/5w07v7hdMJBwAWixwJAKuTJwFYZGu9VH9vkm9P8pZx+NFJ3pnkyo0ICgAWiBwJAKuTJzfYGaefln1LS1Np67qDB5OcMJW2ALaCtRZO75Tk/t39uSSpqucn+f3u/smNCgwAFoQcCQCrkyc32L6lpew/ZzrFzhPPunYq7Wq8AEAAACAASURBVABsFbdY43x7knxp2fCXkpwy9WgAYPHIkQCwOnkSgIW11jNOX5vkwqr6wySd5AeTvGbDogKAxSFHAsDq5EkAFtaaCqfd/atV9WdJvnsc9ePd/f6NCwsAFoMcCQCrkycBWGRrvVQ/SY5Pcm13vyTJvqq6xwbFBACLRo4EgNXJkwAspDWdcVpVz8vwNMRTk7wyyS2TvC7JAzcuNLaa6687mF07T5pKW572CMwLORIAVidPArDI1nqP0x9Mcr8kFydJd++vqhM3LCq2pBtuvNHTHoGtSI4EgNXJkwAsrLVeqv+l7u4MN/NOVTnVDwAGciQArE6eBGBhrbVw+ntV9VtJbl9V/zrJXyT57Y0LCwAWhhwJAKuTJwFYWGu6VL+7X1BVD0tybYZ70/yH7j5/QyMDgAUgRwLA6uRJABbZEQunVXVckj/v7ocmkeAAYCRHAsDq5EkAFt0RL9Xv7q8k+XxVTedx6ACwRciRALA6eRKARbemS/WTfCHJB6vq/CTXHRrZ3c/YkKgAYHHIkQCwOnkSgIW11sLpn44vAOCm5EgAWJ08CcDCOmzhtKr2dPcnuvvVmxUQACwCORIAVidPArAVHOkep28+9Kaq/mCDYwGARSJHAsDq5EkAFt6RCqe17P09NzIQAFgwciQArE6eBGDhHalw2qu8B4DtTo4EgNXJkwAsvCM9HOqfVdW1GY4W3mZ8n3G4u/t2GxodAMwvORIAVidPArDwDls47e7jNisQAFgkciQArE6eBGArONIZp8BRuP66g9m186SptHXy7t258JJLp9IWAAAAAEdH4RSm6IYbb8z+c06YSlu7zl6aSjsAAAAAHL0jPRwKAAAAAGDbUTgFAAAAAJigcAoAAAAAMME9TufEGaefln1L07mn5XUHDyaZzn02Adh8VXX7JC9Lct8kneQnuvvds40KAOZDVR2X5KIkS939qFnHA8DWpXA6J/YtLU3toUInnnXtVNoBYGZekuR/dffjq+pWSY6fdUAAMEeemeSyJLebdSAAbG0u1QeAOVJVt0vyPUleniTd/aXuvma2UQHAfKiqk5M8MsOVGQCwoRROAWC+3DPJgSSvrKr3V9XLqsr9VwBg8OIkz05y46wDAWDrUzgFgPmyI8n9k/zP7r5fkuuSPGf5DFV1ZlVdVFUXHThwYBYxAsCmq6pHJbm6u993hPnkSQCmQuEUAObLviT7uvu94/B5GQqpX9Xd53b33u7eu3Pnzk0PEABm5IFJHlNVVyR5Y5IHV9XrJmeSJwGYFoVTAJgj3f3JJFdW1anjqIckuXSGIQHAXOju53b3yd19SpInJHlHdz95xmEBsIXtmHUAAMDNPD3J66vqVkk+luTHZxwPAADAtqNwCgBzprsvSbJ31nEAwLzq7guSXDDjMADY4lyqDwAAAAAwQeEUAAAAAGCCwikAAAAAwASFUwAAAACACQqnAAAAAAATFE4BAAAAACbMrHBaVcdV1fur6k9mFQMAAAAAwEpmecbpM5NcNsPlAwAAAACsaCaF06o6Ockjk7xsFssHAAAAADicWZ1x+uIkz05y44yWDwAAAACwqk0vnFbVo5Jc3d3vO8J8Z1bVRVV10YEDBzYpOgAAAACA2Zxx+sAkj6mqK5K8McmDq+p1kzN197ndvbe79+7cuXOzYwQAAAAAtrFNL5x293O7++TuPiXJE5K8o7ufvNlxAAAAAACsZlb3OAUAAAAAmFs7Zrnw7r4gyQWzjAEAAAAAYNJMC6cAADBrZ5x+WvYtLU2lrZN3786Fl1w6lbYAAJgthVMAALa1fUtL2X/OCVNpa9fZ0ynAAgAwe+5xCgAAAAAwQeEUAAAAAGCCwikAAAAAwASFUwAAAACACQqnAAAAAAATFE4BAAAAACYonAIAAAAATFA4BQAAAACYsGPWAQAAAABshuuvO5hdO0+aSlvXHTyY5ISptDWvzjj9tOxbWppKWyfv3p0LL7l0Km3BZlE4BQAAALaFG268MfvPmU6x88Szrp1KO/Ns39LS1Ppr19nTKcDCZnKpPgAAAADABIVTAAAAAIAJCqcAAAAAABMUTgEAAAAAJiicAgAAAABMUDgFgDlTVcdV1fur6k9mHQsAAMB2pXAKAPPnmUkum3UQAAAA25nCKQDMkao6Ockjk7xs1rEAAABsZztmHQAAcBMvTvLsJCeuNkNVnZnkzCTZs2fPJoUF63fG6adl39LSVNo6effuXHjJpVNpCwAAVqJwCgBzoqoeleTq7n5fVT1otfm6+9wk5ybJ3r17e5PCg3Xbt7SU/eecMJW2dp09nQIsAACsxqX6ADA/HpjkMVV1RZI3JnlwVb1utiEBAABsTwqnADAnuvu53X1yd5+S5AlJ3tHdT55xWAAAANuSwikAAAAAwAT3OAWAOdTdFyS5YMZhAAAAbFvOOAUAAAAAmKBwCgAAAAAwQeEUAAAAAGCCwikAAAAAwASFUwAAAACACQqnAAAAAAATFE4BAAAAACYonAIAAAAATFA4BQAAAACYoHAKAAAAADBB4RQAAAAAYILCKQAAAHOvqu5WVX9ZVZdV1Yer6pmzjgmArW3HrAMAAACANbghybO6++KqOjHJ+6rq/O6+dNaBAbA1OeMUAACAudfdV3X3xeP7zyW5LMnu2UYFwFbmjFMAYK6dcfpp2be0NJW2rjt4MMkJU2lr2q6/7mB27TxpKm3N83oCTENVnZLkfkneu8K0M5OcmSR79uyZ2jKnmY9O3r07F17iRFng8PzuzJ7CKQAw1/YtLWX/OdMpAp541rVTaWcj3HDjjdtiPQHWq6pum+QPkpzV3Tf7wevuc5OcmyR79+7taS13mvlo19nTKYQAW5vfndlTOAXYYI4SAgBMR1XdMkPR9PXd/aZZxwPA1qZwCrDBHCUEAFi/qqokL09yWXe/cNbxALD1bfrDoarqblX1l1V1WVV9uKqeudkxAAAAsHAemOQpSR5cVZeMrx+YdVAAbF2zOOP0hiTP6u6Lq+rEJO+rqvO727WnAAAArKi735WkZh0HANvHpp9x2t1XdffF4/vPJbksye7NjgMAAAAAYDWbXjhdrqpOSXK/JO+dZRwAAAAAAMvNrHBaVbfN8DTEs7r72hWmn1lVF1XVRQcOHNj8AAEAAACAbWsmhdOqumWGounru/tNK83T3ed2997u3rtz587NDRAAAAAA2NY2vXBaVZXk5Uku6+4XbvbyAQAAAACOZBZnnD4wyVOSPLiqLhlfPzCDOAAAAAAAVrRjsxfY3e9KUpu93EPOOP207FtamkpbJ+/enQsvuXQqbQEAAAAA82PTC6eztm9pKfvPOWEqbe06ezoFWAAAAABgvszk4VAAAAAAAPNM4RQAAAAAYILCKQAAAADABIVTAJgjVXW3qvrLqrqsqj5cVc+cdUwAAADb0bZ7OBQAzLkbkjyruy+uqhOTvK+qzu/uS2cdGAAAwHbijFMAmCPdfVV3Xzy+/1ySy5Lsnm1UAAAA24/CKQDMqao6Jcn9krx3tpEAAABsPy7VB4A5VFW3TfIHSc7q7msnpp2Z5Mwk2bNnz9SWecbpp2Xf0tJU2jp59+5ceIm7C2wF1193MLt2njSVtq47eDDJCVNpi6OzHf59b4d1BAA2l8IpsGXYYWKrqKpbZiiavr673zQ5vbvPTXJukuzdu7entdx9S0vZf850ilq7zp7Ov0Vm74Ybb5zadnHiWdceeSY2xHb4970d1hEA2FwKp7ANTLOgmMxvUdEOE1tBVVWSlye5rLtfOOt4AAAAtiuFU9gGpllQTBQVYYM9MMlTknywqi4Zx/1id791hjEBAABsOwqnADBHuvtdSWrWcQAAAGx3t5h1AAAAAAAA80bhFAAAAABggsIpAAAAAMAEhVMAAAAAgAkKpwAAAAAAExROAQAAAAAmKJwCAAAAAExQOAUAAAAAmKBwCgAAAAAwQeEUAAAAAGCCwikAAAAAwIQdsw5gkV1/3cHs2nnSVNq67uDBJCdMpS0AAAAAYH0UTtfhhhtvzP5zplPsPPGsa6fSDgAAAACwfi7VBwAAAACYoHAKAAAAADDBpfoAwNS5DziwFmecflr2LS1NpS2/FQDAtCmcAjNlhwm2JvcBB9Zi39KS3woAYG4pnAIzZYcJAAAAmEfucQoAAAAAMEHhFAAAAABggsIpAAAAAMAEhVMAAAAAgAkeDgVz6vrrDmbXzpOm0panzQMAAAAcHYVTmFM33Hijp81zM9MsqJ+8e3cuvOTSqbQFAAAAW43CKcACmWZBfdfZS1NpBwAAALYi9zgFAAAAAJigcAoAAAAAMEHhFAAAAABggsIpAAAAAMAEhVMAAAAAgAkKpwAAAAAAE2ZSOK2qh1fV31XVR6vqObOIAQDmlTwJACuTIwHYTJteOK2q45L8f0kekeS0JE+sqtM2Ow4AmEfyJACsTI4EYLPN4ozTM5J8tLs/1t1fSvLGJI+dQRwAMI/kSQBYmRwJwKaaReF0d5Irlw3vG8cBAPIkAKxGjgRgU1V3b+4Cq34oyfd390+Ow09JckZ3P31ivjOTnDkO3jfJhzY10Om5U5JPzzqIdRD/7Cxy7In4Z2mRY0+SU7v7xFkHMStryZNbKEcmi7+9LnL8ixx7stjxL3LsifhnSY60L7lIxD87ixx7Iv5ZWuTYkw3Ikzum2dga7Utyt2XDJyfZPzlTd5+b5NwkqaqLunvv5oQ3XYsceyL+WVrk2BPxz9Iix54M8c86hhk7Yp7cKjkyEf8sLXLsyWLHv8ixJ+KfJTnSvuQiEf/sLHLsifhnaZFjTzYmT87iUv2/SXLvqrpHVd0qyROSvGUGcQDAPJInAWBlciQAm2rTzzjt7huq6meT/HmS45K8ors/vNlxAMA8kicBYGVyJACbbRaX6qe735rkrUfxkXM3KpZNsMixJ+KfpUWOPRH/LC1y7Mnix79uR5knF72/xD87ixx7stjxL3LsifhnaZFjnwr7kgtF/LOzyLEn4p+lRY492YD4N/3hUAAAAAAA824W9zgFAAAAAJhrm144raqHV9XfVdVHq+o5K0z/uqr63XH6e6vqlHH8KVV1fVVdMr5euuwz/7yqPjh+5r9VVc1h/E9aFvslVXVjVZ0+TrtgbPPQtDvPMP7vqaqLq+qGqnr8xLSnVtXl4+upy8ZvSv8fa+xVdXpVvbuqPlxVf1tVP7Js2quq6uPL+v70jYh9PfGP076yLMa3LBt/j3E7u3zc7m41T7FX1fdNbPdfqKrHjdPmqe9/vqouHbePt1fV3ZdNm+l2v57452HbX2ffz3S7n5U19Nnc5sl1xC5Hzij2efidWE/847SZ/1aso/+/r+TJmcS+QNu+PLnMGvprbnPkOuOfeZ481t+5cdoi/H9anrQveazxy5EbZJ19P73tvrs37ZXhBt5/n+SeSW6V5ANJTpuY52eSvHR8/4Qkvzu+PyXJh1Zp98Ik35GkkvxZkkfMW/wT83xrko8tG74gyd456f9Tknxbktckefyy8XdM8rHx7x3G93fYrP5fZ+z3SXLv8f2uJFcluf04/Krl885j34/TDq7S7u8lecL4/qVJfnreYp/Yhj6T5Pg57PvvWxbXT+drvzsz3e6nEP9Mt/31xD7r7X5WrzX22VzmyfXEPjGPHLm5scuRM45/YjuSJzcv9kXZ9uXJo+uvucyR641/Yp5Nz5NrjP2UzGGOnEL8i/JbsWL84zT7khsbvxw5Z30/7e1+s884PSPJR7v7Y939pSRvTPLYiXkem+TV4/vzkjzkcNX3qrprktt197t7WPPXJHnc9ENPMr34n5jkDRsU4+EcMf7uvqK7/zbJjROf/f4k53f3Z7r7s0nOT/LwTez/Y469u/9vd18+vt+f5OokOzcgxsNZT9+vaNyuHpxhO0uG7W6u+n7C45P8WXd/fgNiPJy1xP+Xy+J6T5KTx/ez3u7XFf8cbPvr6fsVbeJ2PyuLnCflSDnyWC1yjkzkyXnf9uc1Ryby5NFa5ByZLHaeXOQcua74F+i3Yl7zpBwpRx6rucmRm1043Z3kymXD+8ZxK87T3Tck+ackXz9Ou0dVvb+q/qqqvnvZ/PuO0Oa0rDf+Q34kN092rxxPIT57o07TztriP9rPblb/ryf2r6qqMzIcrfj7ZaN/dTy1+0VV9XXrC3NV643/1lV1UVW959DlCRm2q2vG7exY2lyrqfR9hqPmk9v9PPb90zIc9TvcZ+ftd2e55fF/1Yy2/fXGPsvtflYWOU/KkTf/7Dz1/RHJkcdMnrz5Z+d125+nHJnIk0drkXPkTWI7zLLmNU8uco48XAxHZYF+KybZlzx2cuTibPcbliM3u3C60o94r3Geq5Ls6e77Jfn5JL9TVbdbY5vTsp74h4lV/0+Sz3f3h5ZNf1J3f2uS7x5fT1lvoKtYT1+t9tnN6v91L2c8svPaJD/e3YeOZj03yTcl+fYMp9D/wnqCPNziVxh3NPHv6e69Sf5VkhdX1b2m0OZaTavvvzXJny8bPXd9X1VPTrI3yX85wmfn7XdnmPHm8R8aP6ttf72xz3K7n5VFzpNy5M0/O099f/gG5Mj1kCdv/tm52/bnMEcm8uTRWuQcebjY1jzPDPPkIufIw8Ww9gYW5LdiFfYlj50cuQDb/UbnyM0unO5Lcrdlwycn2b/aPFW1I8lJST7T3V/s7n9Mku5+X4Zq933G+ZefjrtSm9NyzPEvm36zIyXdvTT+/VyS38lwSvJGWEv8R/vZzer/9cSe8T9Gf5rkl7r7PYfGd/dVPfhikldmPvv+0Onx6e6PZbiP0f2SfDrJ7cft7KjbPArrin30w0n+sLu/fGjEvPV9VT00yb9P8pgxpsN9dt5+d1aLf9bb/rpin/F2PyuLnCflyJt/dp76flVy5LrJkzf/7Fxt+/9/e/cfbEdd3nH8/QkJ5AeEGJJRrLTxWiQFKrFGBwRrcNBqKioDbUbQISC2FGma2NZpm9amUCtCFQqMw1AK1x/gEAtqSgZQSIIUCGBifiESIGoLUrQprSakSOHpH9/nJJuTc0/OPefmnnNzP6+Znbu7Z/e7z+7Zu8/Z7+5+t0dzZMfx98C+P9xGco7cLbYmy+rVPDmSc2SzGFoyUo4VA/G5ZEecI3t8vx+WHBn7uEHXageMpTSI+1p2Ne56TN00H2P3BrGXZv904IDs7wOeBqbm8MPA8exqWHdur8Wfw2Pyy++rK3Na9o+jtLVwfrfir0zbz56Nev+A0qjxK7J/2LZ/h7EfCNwNLGww7eH5V8AVwCU9uO1fARyU/dOAx8lGkYGvsnvDxhf0UuyV8auBk3t121MOok+SDWD3yn4/BPF3dd/vMPau7vfd6lrcZj2ZJzuJPYedI7sTu3NkF+OvjHeeHP7YR8S+3yT+ru/7w921uL16Mkd2Gn8Ody1PthJ7Zdp+eihHDkH8I+JY0SR+n0vu4/hxjuzFbT+k+/2Qr1wLKz8X2JwrtzjHXUSpHQYYnyvyBOVNY305/nTgkdxYa4FTK2XOBjZlmVcD6rX487M5wOq68iYBa4ANuX7/QCb1LsX/ZkpC3g5sBR6pzHturtcTlFu1h3X7txs78CHgRWBdpZuVn60ANmb8XwYO7rVtD7w1Y1yffz9SKbMv97Mncr87qJdiz89mUH6cjqkrs5e2/V3As5X9Y1mv7PedxN8L+34HsXd9v+9W18I269k82W7s+dkcnCOHPXZ64DjRYfw9cazocN+ZgfPksMc+gvZ958nBba+ezZGdxJ+fzaGLebKF2Hs2R3YSPyPnWNGzebLDfWcGzpHDHvsI2u+HJUcqZzQzMzMzMzMzMzOzNNxtnJqZmZmZmZmZmZn1PFecmpmZmZmZmZmZmdVxxamZmZmZmZmZmZlZHVecmpmZmZmZmZmZmdVxxamZmZmZmZmZmZlZHVec2qgi6SVJ6yRtkvRVSRP34bIuknRK9i8c7LJUrJA0OYcXSHpU0o1DENt8Sa+uDF8n6eg2y7pQ0jmdxmRmNhpV8tJ6SWslvbXNcmZIOnMv00yRdEF7kbYUw1GSVuX6PCrp2n21rFzeHEm3DVFZ8yX9NGP/vqRFQxGXpPdJ+rM2ypglaW67MTQo71OS/l3StqEq08xsX3OO7Gh5zpGtlTVR0vJcr0ckXTIU5dr+xRWnNtrsiIhZEXEs8Avg/FZmykrMQf2/RMQnI+KuHFwIDLaSdi6wPiJ+lsMXAHMj4qy62MYOslyA+cDOitOIOC8ivtdGOQDXAwvanNfMbLSr5aXjgD8HPt1mOTOApieFwBRKLtmDpAPaXG7VlcDluT6/Blw1BGUOp5sjYhZwIrBY0hGdFhgRyyKinZOwWZTfAUPlX4C3DGF5ZmbDwTmyd+zPOfLvI2Im8EbgREnvGcKybT/gilMbze4FfhVA0sfzLtRNkhbmuBl5NfDzwFrgCEkflLQxp/tMTneApP4ct7F2BS7HnSFpAaWScqWklZI+IunyWhCSPirpcw3iOwv4Rk5zDdAHLJO0SNISSddK+ibwxYz13rwSu9vVWEmfyLjWS7pE0hnAbODGvGo4Ia9+zs7p91jHHL8t71hZL2m1pFcCRMTzwA8l+YTMzKwzk4HnYOcFu8squWVes/HAJcDb8ri+SNIxkh7K4Q2SjsxpXpfjLsu7PlZKugnYmOV/XdKavOvi92qBZQ74bOaYuyVNbxD/4cBTtYGIqJXZMEfl8u+RtFTS5sxRZ2XcGyW9Lqfrl3RNlrFZ0nvrFyxpkqTrJT0s6buS3t/ulxARW4Encn2QNF3SLVn2w5JOzPFvkXR/Lu9+SUc1iGu+pKuzf12l2yHp7Y3KkHQgcBEwL6edJ2lqfjcbMge/Ictckuu9StIWld8cjdZpdUQ80+42MTPrAc6RzpFDniMj4vmIWJn9v6Cc97+m3e1j+6mIcOdu1HTAtvw7llIp+QfAmyjJcBJwMPAI5WrTDOBl4Pic59XAvwHTc/4VwAdy/m9VljEl//YDZ2T/D4Fp2T8JeBIYl8P3A7/eINYfAYdUhqtlLAHWABNyeCIwPvuPBL6T/e/J8ifm8NT8uwqYXSl7FaUyteE65jQBnJr9lwJ/WZl/MfDH3f5+3blz526kdcBLwDrg+8D/AG/K8acD3wIOAF6Zx+bDm4yfA9xWKfcq4KzsPxCYkHltU2WaOcB24LWVcbU8MQHYBByWw1Ep75PA1Q3W5Zxch9uBRZV8OFCOmgP8d8Z/EPA08Df52R8BV2R/P3AH5YL/kZQTz/HVdQb+DvhQ9k8BNgOTBvE9zK+tE/DL+Z3UYr4JOKny2aPZPxkYm/2nALdU1uu2+nIryzqVcvF2XJMydpsvv8+/zv53AOuyfwklzx8ETAO2kr8vBljPbd3e5925c+eu1Q7nSOfI4c2RU4AtQF+39313vdW184iv2Ug2QdK67L8X+CdK5enXImI7gKRbgbcBy4AfRcTqnP7NwKqI+GlOdyPwm8DFQJ+kq4DlwDebBRAR2yWtAN4r6VHKwXtjg0mnRsTPmxS1LCJ2ZP844GpJsyg/MF6f408BbohyVygR8V/NYmuyjl+nNG1QaydnDfDOynw/AWbupWwzM9vTjiiPviHpBMpTBMcCJwFfiYiXgGcl3UM5Rg80/md15T5AeZTuNcCtEfG4pEbLfygiflAZXiDptOw/gnIStpVyIfHmHP9l4Nb6giLiBkl3Au8G3g/8vqTjGDhHATwceSekpCfZlUM3AidXplsaES8Dj0vawp45513A+yT9SQ6PJ0/gGq30AOZJOhk4CvhoRPxvjj8FOLqy/SZLOgQ4FPhC3qkUuZ5N5bSXAe+IiBclvarFMk6iVAgQESskHSbp0PxseUS8ALwg6SeUyoKnBijHzGwkcY50jhyWHKnS/N1XgCsjYsveYrXRxRWnNtrsTL41GiBLpu3VSRtNEBHPZdL7LeBjwO8C5+4ljuuAv6BcPb1hgGn+T9KYTIJ7i20R8CxwHOVqYy2RiZJkWtVsW7wYEbWyXmL348d4YMees5iZWasi4gFJ0yh3/Q90PG52nK6WdZOkB4HfBu6UdB7lLop6O3OJpDmUE6ATIuJ5Sasox/eGixhguT+mtH19vaRNwLGUu0ca5SiAFyr9L1eGX2b3PFO/vPphAadHxGMDxIukT1G2B/W/BdLNEXFhnpwvl3R7RPxHxnxC5WJlrbyrgJURcZqkGZSnNwYkaRKwlHLC+eMcfXGLZTT63mvboLoN6/Ozmdl+wTnSObJJGUORI68FHo+IK5rFaaOT2zg1g28DH1B5o94k4DTK3aj1HgTeLmmaSgPhHwTuyQQ+JiJuAf4K+I0G8/4cOKQ2EBEPUq5Snkm5stXIY5R2TVtxKPBMVrJ+mPJ4CpSrkudKmgggaWqjePa2ji0s//WUx1XMzKxNkmZSjt9bKblpnko72tMpd/8/1GT8bsd1SX3Aloi4kvIExRvqp2ngUOC5PCGcCRxf+WwMcEb2nwn8a4P43y1pXPa/CjiM8mjhQDlqMH5H0hiVNt36KDmy6k7gD2sXQyW9sb6AiFgc5aUcjU4Iq9M9AHyJ8igklFx6Ye3zvCsIyno9nf3zW1iHGyhPgVR/YwxURv139W1K2+e1k/f/jF0vjzQz2+85RzblHNlBjpT0t7msha3OY6OLK05t1IuItZS2YR6iVBxeFxHfbTDdM5S3Oa4E1gNrI+IbwC8Bq1SaAOjPaepdC9wuaWVl3FLgvoh4boDQllPagGnF54GzJa2mVGJuz5jvoPwY+E7GV3s8ox+4RvlyqBbWcW9OBO5qMVYzM9tlQh6L11Ee8zs7HzH8GrCBcixeAXwi7+wYaPwGypMK61VeUjgP2JTlzgS+GOWFDvepvDTjsgax3AGMlbSBcpfH6spn24FjJK2htB92UYP535XLXE85SfvTjK1hjhqkxygX8m4Hzq88IlhzMeURvg15F8/FbSyj6jPAOfm44QJgtspLJ74HnJ/TXAp8WtJ97OVEaaGDOQAAARtJREFUV9KvUE6qz9Wul1/MblLGSsqjj+tUXm6ypBYD5QUmZw9mZSRdKukpYKKkpyQtGcz8ZmZd4hzZGufINnOkSnMNi4GjgbVZ5nmtzm+jg3Y9eWtmw0nSbcDlEXH3AJ8fTkni72z0ea/IK5Yfj4gPdzsWMzPbNyRti4iDu7TsfsqLJP65G8s3MzNrxjnSbP/mO07NhpmkKZI2U9pbbVhpCjvv/vxHSZOHL7q2TKM0UWBmZmZmZmZmtt/wHadmZmZmZmZmZmZmdXzHqZmZmZmZmZmZmVkdV5yamZmZmZmZmZmZ1XHFqZmZmZmZmZmZmVkdV5yamZmZmZmZmZmZ1XHFqZmZmZmZmZmZmVkdV5yamZmZmZmZmZmZ1fl/uAOe+a4D/EQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "samples1 = random.choices(df['Porosity'].values, weights=df['Wts'].values, cum_weights=None, k=len(df))\n",
    "samples2 = random.choices(df['Porosity'].values, weights=df['Wts'].values, cum_weights=None, k=len(df))\n",
    "\n",
    "print('Bootstrap means, realization 1 = ' + str(np.average(samples1)) + ' and realization 2 = ' + str(np.average(samples2)))\n",
    "\n",
    "plt.subplot(131)\n",
    "GSLIB.hist_st(df['Porosity'],pormin,pormax,False,False,20,df['Wts'],'Porosity (fraction)','Histogram Declustered Porosity')\n",
    "\n",
    "plt.subplot(132)\n",
    "GSLIB.hist_st(samples1,pormin,pormax,False,False,20,None,'Bootstrap Sample - Realizaton 1','Histogram Bootstrap Porosity 1')\n",
    "\n",
    "plt.subplot(133)\n",
    "GSLIB.hist_st(samples2,pormin,pormax,False,False,20,None,'Bootstrap Sample - Realizaton 2','Histogram Bootstrap Porosity 2')\n",
    "\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=3.0, top=1.1, wspace=0.2, hspace=0.2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the bootstrap distributions vary quite a bit from the original.\n",
    "\n",
    "#### Summarizations Over Bootstrap Realizations\n",
    "\n",
    "Let's make a loop to conduct $L$ resamples and calculate the average and standard deviation for each  ($m^\\ell$, $\\sigma^2_{\\ell}$,  for $\\ell = 0,\\dots,L-1$). We then summarization over these $L$ realizations.  \n",
    "\n",
    "I did not find any built-in, concise functions to accomplish this, i.e. with a single line of code, so we are going to do it by hand.  \n",
    "\n",
    "To understand this code there are just a couple of Python concepts that you need to add to your Python arsenal.\n",
    "\n",
    "1. declaring arrays - NumPy has a lot of great array (ndarray) functionality. There are build in functions to make a ndarray of any length (and dimension). This includes 'zeros', 'ones' and 'rand', so when we use this code:\n",
    "\n",
    "```p\n",
    "mean = np.zeros(L); stdev = np.zeros(L)\n",
    "```\n",
    "\n",
    "    we're making arrays of length $L$ pre-populated with zeros.\n",
    "\n",
    "2. For Loops - when we are using the command below, we are instructing the computer to loop over all the indented code below the command for $l = 0,1,2,\\ldots,L-1$ times. For each loop the $l$ variable increments, so we can use this to save each result to a different index in the arrays mean and stdev. Note, Python arrays index starting at 0 and stop at the length - 1.\n",
    "\n",
    "```p\n",
    "for l in range(0, L): \n",
    "```\n",
    "\n",
    "    we are running each bootstrap resampled realization, calculating the average and standard deviation and storing them in the arrays that we already declared."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABU4AAAGWCAYAAACn5xmpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdf5xeZ10n/M+XprU0hJYfoTYTQkH7IBGl8MSCi2il4AoKxZUi4I8IxeqqaFf2oZW1Up+ClteqyLorWkAIP4QWBIqorFipLq7bWCCKpvIgUEomoQ3Qmk4olDTX88c5oXfvziT3JHNmJpn3+/U6rzk/rnPO95w599zXfM+5rlOttQAAAAAAcLf7LHUAAAAAAADLjcQpAAAAAMAYiVMAAAAAgDESpwAAAAAAYyROAQAAAADGSJwCAAAAAIyROAVWpKr656o6exH2s6GqZqrquKH3NYSqelJVfWKp4wAAOFaoh07mWKqHVtWlVfXWBdzeT1bVhxdqe/Pc949W1V8cwfp/XlWbFzImGJLEKRyhqrqxqu7oKyW3VtWfVtVDF2C7p1dVq6pVE5ZvVfXNR7rfSc31Zd2fj6csVhz9Pud1rpKktfatrbVrJ9z+YR9Ta+2m1tr9Wmt3zXfdqjq7qvb319btVfWJqnrB4cRxuFpr/6u19siRmI7491tVq/tj+rMjjxAAVi710HvNVw+9537UQ0dU1flV9S/98dzcf17W9MveVFWvWIi4l9rINTnTDzdX1fur6qkLsf3W2ttaa983YSz3Shi31p7WWtuyELHAYpA4hYXxjNba/ZKcluTmJL+7xPHcy3wqc0eTY/W4Ruzsr637J7koyeuqauN8NrAMz9Gzk3w1yfdV1WlD7GAZHjMADEU9dIkcq8c14piph1bV9yT59STPa62tSfKoJFctbVSTOYInhk/pf3+PSfLBJO+pqp9csMBghZA4hQXUWvtKkncl+XqFoqpOrqo3V9XuqvpsVf1KVd2nX3affvqzVXVLX+7kftW/6X/e1t8p/M6q+uaq+uuq+req+kJVXdlv50DZf+jL/kh/l3hHVV1UVZ9P8saqekB/t3F3/1TC+6tq/Uis11bVb1TV1n4fV1fVAw/3fPR3bv9Hfzf39qq6rqq+aWT5t1bVB6vqS/2d0JeNnJeLq+pTVfXFqrrqQBwjd1DPr6qbkvzVHOfqm6rqr/r1v1BVb6uqU0b2/fU71v2d0Kv68397dc2nNvXL3pJkQ5I/6bf90v54Xjx2rP9YVc+a5Rzc4ymE/hxfVlV/2+/rL6rqwYc6l63z3iS3pr++quqZfay39dt91NjxXVRV/5hkb1WtqqpH9eVu69d75kj5p1fV9j6m6ar6z/38s6tqx0KcixGbk/x+kn9M8qMj611cVe8a29Zrquq/9eMnV9UbqmpXH+Mrqq9IVvfkyd9W1aur6ktJLp3gGnhcVX2sP+Z3VtWVNfKkQVX9YFVt68/X/66qbz/U7wkAlop66D2Veqh66D19R5K/a619rD+mL7XWtrTWbq+qC9LVSV/ab/tP+m0duA5u7+P7oZH9/GRVfbiqfrO/nj9TVU8bWf7w/vNye1V9MMk9znN1dc/P99f631TVt44se1NVvbaq/qyq9ib53qp6UFW9r6r2VNXWJN+UCbXWPt9ae02SS5O8qu7+G7Cuqv64/0x+pqp+YWT+HaOfv6p6bH8tH19jT3xXV1//XB/bR6rqSf3870/ysiQ/0p/Xf+jnX1tVL+rH5/w7NHL9bq6qm/r9/5dJjxsWTGvNYDAcwZDkxiRP6cdPSrIlyZtHlr85ydVJ1iQ5Pcn/l+T8ftkLk/xrkkckuV+Sdyd5S7/s9CQtyaqRbb09yX9Jd9PjxCTfNbKsJfnmkemzk+xL8qok35DkvkkelOSH+zjXJHlnkveOrHNtkukkj06yOskfJ3nrHMf9k0k+fIjz8aYkX0pyVpJVSd6W5B39sjVJdiV5SX8sa5I8vl92YZL/k2R9H/sfJHn72Hl5cx/jfec4V9+c5Kn9+mvTVWp/Z444L03ylSRPT3Jckt9I8n9mK9tPPyfJdSPTj0nyxSQnzHI+7hFbf44/leT/6mO/Nsnlc5zjs5Ps6Mfvk+SHknwtySP79ff2x3h8kpemu5ZOGIl5W5KH9vs5vl/+siQnJHlyktuTPLIvvyvJk/rxByR53HgMR3ou+uUbkuxPV+l+SZJ/HFn2sCRfTnL/fvq4Pq4n9NPv7a+F1UkekmRrkp8euR73JXlxumvtvge7Bvpz8Nkkv9ifm/+Q5M4kr+iXPy7JLUke38exuT/2b1jqvzkGg8FgMBwYoh56sPPxpqiH3iO2rOB6aJInJbkjya8leWLG6nTprpdXjM07L8m6/vh/pD/m00auwa8l+an+9/Yfk+xMUv3yv0vy2/018N398b51ZNsvTHfdfUOS30mybSyWf+vjPPB5e0e6J2RXp/uMTGeWz8Bcn99+/iP6+Y/qt/uRJL/a/04ekeTTSf59X/avkvzUyLr/Ncnvz/b5S/Jj6T7fq9J9pj6f5MSR6/utY3Fcm+RF8/g79Lp019Fj0rVae9RS/+01rKxhyQMwGI72of8Cn0lyW7oK4s4k39YvO67/475xpPxPJ7m2H78myc+OLHtk/wW8arYvvHSVtCuSrJ8ljtkqrHce+NKaI/Yzk9w6Mn1tRipP6ZJbdyY5bpZ17/GFOXY+Riusrx9Z9vQk/9KPPy/Jx+aI64Yk54xMnzbLeXnEyPJ7natZtvms0f3l3hXWvxw77jtmK9tPf0O6ivgZ/fRvJvm9OfZ7j9j6c/wrI8t/NskH5lj37HRJxtv6/W1L8tx+2SVJrhope590FaizR2J+4cjyJ6WrxNxnZN7bk1zaj9+U7tq8/ywxHKzCOvG56Jf/SvqKYbqK6F1JHjuy/MNJfqIff2qST/Xjp6b7LN13pOzzknxo5Hq86RCf1a9fA+kqsNPpK7cj+z6QOH1tksvG1v9Eku852D4MBoPBYFjMIeqhs50P9dA5Yot66NOS/El/TDPpEpvHjVwvr5hr3b7MtiTnjlyD/zqy7KT+XH9jugcF9iVZPbL8jzL3jYBT+nVPHoll9AbIcemuwW8ZmffrmX/i9MR+/hPTPRxw09jyX07yxn78RUn+qh+vJJ9L8t0H+/yNbOfWJI8Zub4Pljid5O/Q+pHlW9NfhwbDYg2a6sPCeFZr7ZR0X94/n+Svq+ob0zXJOPBk2wGfTTLVj6+bZdmqdEmi2bw03RfX1r6JywsPEdfu1jXbSpJU1UlV9Qd9U4g96e5+n1L37Dfnc2PxHJ+xpiW9ff2yccen+7I74PMj419Odycx6e5Af2qOuB+Wrg+e26rqtnQV2Ltyz/PyuVnX7FXVQ6rqHX1znz1J3jrHccwV54k1R59MrbWvprvj+2N9U5fnJXnLweI5xL7uN1fBdH1LndJae2Br7czW2jv6+fe4dlpr+9Odk6mRdUfP0bokn+vLHTB6Lf5wun8oPts3K/rOSQ7kMM7FT6R74iOttZ1J/jrd05wH/FG/jSR5fj+ddNfE8Ul2jVwXf5DuydMD7nFNHOIaWJdkurXW5lj/YUlecmBf/f4e2q8HAMuJeujd1EMPbcXWQ1trf95ae0aSByY5N10C8EVzla+qn6i7u226Ld2TnqO/x6+fy9bal/vR+/XHe2trbe/Y8R7Y7nFVdXnfDcCedAnhjG179PytTffZHP98zNeB8/2ldNf5urG67sty93X+riTfWVXr0j1w0JL8r9k2WlUvqaob+m4Hbktycg5+vY+a5O/QfK5ZWHASp7CAWmt3tdbena5y9V1JvpCu8vawkWIb0t2RTbqnAsaX7UvXsf9oQufA9j/fWvup1tq6dHdlf68O/gbT8W28JN1dvMe31u6f7ksw6SrBB4y+iXVDH/8XZtn2TUk2VNXX162qk9Ilsib5Iv9c5u6b53NJntZX1A4MJ7bWpkfKtDnGD/iNfv6398f6Y7nncc7HbNvfkq4vpHOSfLm19neHue3DdY9rp/89PDR3X1vJPePemeShB/o06n39Wmyt/X1r7dx0v7/3Zu7O8g/7XFTVv0tyRpJfrq5Pp8+nu9v9vJF/Dt6Z5Ozq+jz7odydOP1cuqdmHjxyTdy/tfatI7sYj+1g18CuJFOj12/uee1/Lskrx67Bk1prb5/9tADA0lIPVQ9dREddPfQeG2ltf2vtmnTN0R8927ar6mHpmoj/fJIH9Tcn/imT/R53JXlAVa0embdhZPz56RK3T0mXZDz9wG5HwxwZ353uszn++ZivH0rXFdUn0l3nnxm7zte01p6eJK2125L8RbruEJ6frruKe53/vj/Ti/pyD+jP07+NHMtsv7NRB/s7BMuCxCksoOqcm65vnhtaa3el++J/ZVWt6b+AfyndXeeka6Lyn6rrPPx+6ZpcXNla25fuC3J/uv5eDmz/vLq7E/1b030R3dVP3zxadg5r0vXtc1t1nX2/fJYyP1ZVG/vK5/+b5F39cYy7Ll1/TBdX1Yl9xeDyJNdnsgrr+5N8Y1VdWFXf0J+fx/fLfj/dOXtYf9xr+/M6l3udq/5YZ/pjnUry/0wQ01zudW77Stn+JL+V+d3lXyhXJfmBqjqnqo5P98/IV5P87znKX5euX6aXVtep+9lJnpHkHVV1QlX9aFWd3Fr7WpI9ufu6Gnck52Jzujd6bkzXPO/MdJXVk9I1nUprbXe65jtvTFeZu6Gfvytd5e23qur+1XUk/03VvSF1Lge7Bv6uP8afr+6FBeem6wPtgNcl+Zmqenz/uV5dVT9QVWsOsj8AWDLqoeqhi+ioq4dW1blV9dzqXlJWVXVWku9J15/tbNtene4a392v/4LcnWQ9qNbaZ9Ndi7/WH9939cd7wJp05+uL6erBv36I7d2Vru/PS6t7cntj7tli66Cq6tSq+vl0n7lf7p/83ZpkT3Uv8bpvdU/BPrqqvmNk1T9K11rsh3P3wwzj1qRLdO5OsqqqfjXJ/UeW35zk9LGk+aiD/R2CZUHiFBbGn1TVTLov+lcm2dxa++d+2YvTVRQ+na4PxT9K8of9sj9M9wX/N0k+k64C+OLk6809Xpnkb6trPvGEdG+DvK7f1/uS/GJr7TP9ti5NsqUv+5w54vyddB1rfyFdJeEDs5R5S7p+dT6frh+cX5htQ33TmB9I3/dQf3zrkjxntruRs6x/e7o+LJ/R7+uTSb63X/ya/vj+oqpu72N9/Gzb6bc127n6tXQv+Pm3JH+arrJxuH4jya/02/7PI/PfnOTbcvc/IIumtfaJdE8v/G663+czkjyjtXbnHOXvTPLMdAnKLyT5vXR9if5LX+THk9xYXXOhn+m3PZvDOhdVdWK6O9G/2z+xcmD4TLprbry5/lNy7wraT6Rrcrg93T9s70rX79hc5rwG+vPxH5Kcn66fqx9L90/UV/vl16fr7P+/9/v613TNuQBguVEPVQ9dVEdbPbR3a7q63SfTfVbemuS/ttbe1i9/Q5KN/bbf21rbni4Z+3fpkn/fluRvD7L9cc9Pd918KV3C8s1j8X423RO323N38vZgfj5dE/XPp/uMvHGCdW6rqr1JPp6uK4TzWmt/mHw9GfuMdA8yfCbd7+X16Z6APeB96VqL3dxa+4c59vE/k/x5uhfPfTbd35HRLgXe2f/8YlV9dJb15/w7BMtFTfC9AqwQVXVtus67X7/UsRwNquonklzQWvuupY5lqR0L56Kqrkv3ttBJKqIAwAJSD52fY6HutVCcC2BInjgFOAx9E7KfTfd22RXtaD0XVfU9VfWNfVP9zUm+PbM//QIAsGwcrXWvITgXwNAkTgHmqar+fbp+fG7O3P39rAhH+bl4ZJJ/SNeM7iVJnt33pQoAsCwd5XWvBeVcAItBU30AAAAAgDGeOAUAAAAAGLNqqQOYxIMf/OB2+umnL3UYAADL2kc+8pEvtNbWLnUcLCx1YQCAQxuiLnxUJE5PP/30XH/99UsdBgDAslZVn13qGFh46sIAAIc2RF1YU30AAAAAgDESpwAAAAAAYyROAQAAAADGSJwCAAAAAIyROAUAAAAAGCNxCgAAAAAwRuIUAAAAAGCMxCkAAAAAwBiJUwAAAACAMRKnAAAAAABjJE4BAAAAAMZInAIAAAAAjJE4BQAAAAAYI3EKAAAAADBG4hQAAAAAYMyqpQ4AloOzztyYHdPTE5dfPzWVrdu2DxgRAADAyuT/M2C5GDRxWlX/KcmLkrQkH0/ygiSnJXlHkgcm+WiSH2+t3TlkHHAoO6ans/Oy1ROXX3fJ5F/iAAAATM7/Z8ByMVhT/aqaSvILSTa11h6d5Lgkz03yqiSvbq2dkeTWJOcPFQMAAAAAwOEYuo/TVUnuW1WrkpyUZFeSJyd5V798S5JnDRwDAAAAAMC8DJY4ba1NJ/nNJDelS5j+W5KPJLmttbavL7YjydRs61fVBVV1fVVdv3v37qHCBAAAAAC4l8H6OK2qByQ5N8nDk9yW5J1JnjZL0Tbb+q21K5JckSSbNm2atQwAAACwst2xdybr1p48UVkvkgLmY8iXQz0lyWdaa7uTpKreneTfJTmlqlb1T52uT7JzwBgAAACAY9i+/fsnfpmUF0kB8zFkH6c3JXlCVZ1UVZXknCTbk3woybP7MpuTXD1gDAAAAAAA8zZkH6fXpXsJ1EeTfLzf1xVJLkryS1X1r0kelOQNQ8UAAAAAAHA4hmyqn9bay5O8fGz2p5OcNeR+AQAAAACOxJBN9QEAAAAAjkoSpwAAAAAAYyROAQAAAADGSJwCAAAAAIyROAUAAAAAGCNxCgAAAAAwZtVSBwCTOOvMjdkxPT1x+fVTU9m6bfuAEQEAAABwLJM45aiwY3o6Oy9bPXH5dZdMnmQFABhCVZ2S5PVJHp2kJXlhkk8kuTLJ6UluTPKc1tqtSxQiAAAHoak+AAAM4zVJPtBa+5Ykj0lyQ5KLk1zTWjsjyTX9NAAAy5DEKQAALLCqun+S707yhiRprd3ZWrstyblJtvTFtiR51tJECADAoUicAgDAwntEkt1J3lhVH6uq11fV6iSnttZ2JUn/8yGzrVxVF1TV9VV1/e7duxcvagAAvk7iFAAAFt6qJI9L8trW2mOT7M08muW31q5orW1qrW1au3btUDECAHAQEqcAALDwdiTZ0Vq7rp9+V7pE6s1VdVqS9D9vWaL4AAA4BIlTAABYYK21zyf5XFU9sp91TpLtSd6XZHM/b3OSq5cgPAAAJrBqqQMAAIBj1IuTvK2qTkjy6SQvSPfgwlVVdX6Sm5Kct4TxAQBwEBKnAAAwgNbatiSbZll0zmLHAgDA/GmqDwAAAAAwRuIUAAAAAGCMxCkAAAAAwBiJUwAAAACAMRKnAAAAAABjJE4BAAAAAMZInAIAAAAAjJE4BQAAAAAYI3EKAAAAADBG4hQAAAAAYMyqpQ4AODJnnbkxO6anJy6/fmoqW7dtHzAiAAAAgKOfxCkc5XZMT2fnZasnLr/uksmTrAAAAAArlab6AAAAAABjJE4BAAAAAMZInAIAAAAAjBmsj9OqemSSK0dmPSLJryZ5cz//9CQ3JnlOa+3WoeIAAAAAls58X2i7d2YmyeTvcQAYymCJ09baJ5KcmSRVdVyS6STvSXJxkmtaa5dX1cX99EVDxQEAAAAsnfm+0HbNhXsGjAZgcovVVP+cJJ9qrX02yblJtvTztyR51iLFAAAAAAAwkcVKnD43ydv78VNba7uSpP/5kNlWqKoLqur6qrp+9+7dixQmAAAAAMAiJE6r6oQkz0zyzvms11q7orW2qbW2ae3atcMEBwAAAAAwi8V44vRpST7aWru5n765qk5Lkv7nLYsQAwAAAADAxBYjcfq83N1MP0nel2RzP745ydWLEAMAAAAAwMRWDbnxqjopyVOT/PTI7MuTXFVV5ye5Kcl5Q8YAy8FZZ27Mjunpicqun5rK1m3bB44IAAAAgIMZNHHaWvtykgeNzftiknOG3C8sNzump7PzstUTlV13yWQJVgAAAACGsxhN9QEAAAAAjioSpwAAAAAAYyROAQAAAADGSJwCAAAAAIyROAUAAAAAGCNxCgAAAAAwRuIUAAAAAGCMxCkAAAAAwBiJUwAAAACAMRKnAAAAAABjJE4BAAAAAMasWuoAAAAAgKPHWWduzI7p6YnL752ZSbJ6uIDm4Y69M1m39uSJy6+fmsrWbdsHjAhYziROAQAAgIntmJ7OzssmT4SuuXDPgNHMz779++cV+7pLJk8QA8ceTfUBAAAAAMZInAIAAAAAjJE4BQAAAAAYI3EKAAAAADBG4hQAAAAAYIzEKQAAAADAmFVLHQAAAByLqurGJLcnuSvJvtbapqp6YJIrk5ye5MYkz2mt3bpUMQIAMDdPnAIAwHC+t7V2ZmttUz99cZJrWmtnJLmmnwYAYBmSOAUAgMVzbpIt/fiWJM9awlgAADgITfUBAGAYLclfVFVL8gettSuSnNpa25UkrbVdVfWQ2VasqguSXJAkGzZsWKx4AVhkZ525MTumpycu/7U7v5LjTzhxorLrp6ayddv2ww0NiMQpS2g+XxB7Z2aSrB42IACAhfXE1trOPjn6war6l0lX7JOsVyTJpk2b2lABArC0dkxPZ+dlk/+vu+bCPdn9qgdNVHbdJZMnZIHZSZyyZObzBbHmwj0DRwMAsLBaazv7n7dU1XuSnJXk5qo6rX/a9LQktyxpkAAAzEkfpwAAsMCqanVVrTkwnuT7kvxTkvcl2dwX25zk6qWJEACAQ/HEKQAALLxTk7ynqpKuzv1HrbUPVNXfJ7mqqs5PclOS85YwRgAADkLiFAAAFlhr7dNJHjPL/C8mOWfxIwIAYL401QcAAAAAGCNxCgAAAAAwRuIUAAAAAGCMPk45Jt2xdybr1p48cfm9MzNJVg8XEAAAAABHlUETp1V1SpLXJ3l0kpbkhUk+keTKJKcnuTHJc1prtw4ZByvPvv37s/OyyROhay7cM2A0AAAAHI3m+1DO+qmpbN22fcCIgMU09BOnr0nygdbas6vqhCQnJXlZkmtaa5dX1cVJLk5y0cBxAAAAAMzLfB/KWXfJ9IDRAIttsD5Oq+r+Sb47yRuSpLV2Z2vttiTnJtnSF9uS5FlDxQAAAAAAcDiGfDnUI5LsTvLGqvpYVb2+qlYnObW1titJ+p8PmW3lqrqgqq6vqut37949YJgAAAAAAPc0ZOJ0VZLHJXlta+2xSfama5Y/kdbaFa21Ta21TWvXrh0qRgAAAACAexmyj9MdSXa01q7rp9+VLnF6c1Wd1lrbVVWnJbllwBgAAAAAFsV8Xya1d2YmyeR9qAKLa7DEaWvt81X1uap6ZGvtE0nOSbK9HzYnubz/efVQMQAAAAAslvm+TGrNhXsGjAY4UkM+cZokL07ytqo6Icmnk7wgXfcAV1XV+UluSnLewDEAAAAAAMzLoInT1tq2JJtmWXTOkPsFAAAAADgSQ74cCgAAAADgqCRxCgAAAAAwRuIUAAAAAGCMxCkAAAAAwBiJUwAAAACAMRKnAAAAAABjVi11AHA0umPvTNatPXni8ntnZpKsXvJtH872109NZeu27ROXBwAAADgWSJzCYdi3f392XjZ5snLNhXuWxbYPZ/vrLpme1/YBAAAAjgWa6gMAAAAAjJE4BQAAAAAYI3EKAAAAADBG4hQAAAAAYIzEKQAAAADAGIlTAAAAAIAxEqcAAAAAAGMkTgEAAAAAxkicAgAAAACMkTgFAAAAABizaqkD4Nhx1pkbs2N6euLye2dmkqweLiAAAAAAOEwSpyyYHdPT2XnZ5InQNRfuGTAaAAAAADh8muoDAAAAAIyROAUAAAAAGKOpPgAAAMAx5o69M1m39uSJy6+fmsrWbdsHjAiOPhKnAAAAAMeYffv3z+s9JOsumfxlz7BSaKoPAAAAADBG4hQAAAAAYIzEKQAAAADAGIlTAAAAAIAxEqcAAAAAAGMkTgEAYCBVdVxVfayq3t9PP7yqrquqT1bVlVV1wlLHCADA7CROAQBgOL+Y5IaR6VcleXVr7YwktyY5f0miAgDgkCROAQBgAFW1PskPJHl9P11JnpzkXX2RLUmetTTRAQBwKKuG3HhV3Zjk9iR3JdnXWttUVQ9McmWS05PcmOQ5rbVbh4wDAACWwO8keWmSNf30g5Lc1lrb10/vSDI124pVdUGSC5Jkw4YNA4cJAMkde2eybu3JE5dfPzWVrdu2DxgRLL1BE6e9722tfWFk+uIk17TWLq+qi/vpixYhDgAAWBRV9YNJbmmtfaSqzj4we5aibbb1W2tXJLkiSTZt2jRrGQBYSPv278/Oy1ZPXH7dJdMDRgPLw2IkTsedm+TsfnxLkmsjcQoAwLHliUmeWVVPT3JikvunewL1lKpa1T91uj7JziWMEQCAgxi6j9OW5C+q6iN9c6MkObW1titJ+p8PmW3Fqrqgqq6vqut37949cJgAALBwWmu/3Fpb31o7Pclzk/xVa+1Hk3woybP7YpuTXL1EIQIAcAhDJ06f2Fp7XJKnJfm5qvruSVdsrV3RWtvUWtu0du3a4SIEAIDFc1GSX6qqf03X5+kbljgeAADmMGhT/dbazv7nLVX1niRnJbm5qk5rre2qqtOS3DJkDAAAsJRaa9em654qrbVPp6sTAwCwzA32xGlVra6qNQfGk3xfkn9K8r50zZISzZMAAAAAgGVoyCdOT03ynqo6sJ8/aq19oKr+PslVVXV+kpuSnDdgDByBs87cmB3Tk78lb+/MTJLJ38AHAAAAAMvVYInTvhnSY2aZ/8Uk5wy1XxbOjunp7Lxs8kTomgv3DBgNAAAAACyeiZrqV9Wjhw4EAACWI3VhAICVadI+Tn+/qrZW1c9W1SmDRgQAAMuLujAAwAo0UVP91tp3VdUZSV6Y5Pqq2prkja21Dw4aHbDk7tg7k3VrT56o7PqpqWzdtn3giABgcakLAwCsTBP3cdpa+2RV/UqS65P8tySPre7NTy9rrb17qACBpbVv//6J+7pdd8nkLxMDgKOJujAAwMozaR+n315Vr05yQ5InJ3lGa+1R/firB4wPAACWlLowAMDKNOkTp/89yevS3VG/48DM1trO/s47AAAcq9SFAQBWoEkTp1KZsdsAACAASURBVE9Pckdr7a4kqar7JDmxtfbl1tpbBosOAACWnrowAMAKNFFT/SR/meS+I9Mn9fMAAOBYpy4MALACTZo4PbG1NnNgoh8/aZiQAABgWVEXBgBYgSZNnO6tqscdmKiq/zvJHQcpDwAAxwp1YQCAFWjSPk4vTPLOqtrZT5+W5EeGCQkAAJYVdWEAgBVoosRpa+3vq+pbkjwySSX5l9ba1waNDAAAlgF1YQCAlWnSJ06T5DuSnN6v89iqSmvtzYNEBQAAy4u6MLCgzjpzY3ZMT09c/mt3fiXHn3DiIOXnu+29MzNJVk9cHuBoNVHitKrekuSbkmxLclc/uyVRWQQA4JimLgwMYcf0dHZeNnnycc2Fe7L7VQ8apPzhbBtgJZj0idNNSTa21tqQwQAAwDKkLgwAsALdZ8Jy/5TkG4cMBAAAlil1YQCAFWjSJ04fnGR7VW1N8tUDM1trzxwkKgAAWD7UhQEAVqBJE6eXDhkEAAAsY5cudQAAsNzcsXcm69aePHH59VNT2bpt+4ARwcKbKHHaWvvrqnpYkjNaa39ZVSclOW7Y0AAAYOmpCwPAve3bv39eLzhbd8n0gNHAMCbq47SqfirJu5L8QT9rKsl7hwoKAACWC3VhAICVadKXQ/1ckicm2ZMkrbVPJnnIUEEBAMAyoi4MALACTZo4/Wpr7c4DE1W1KkkbJiQAAFhW1IUBAFagSROnf11VL0ty36p6apJ3JvmT4cICAIBlQ10YAGAFmjRxenGS3Uk+nuSnk/xZkl8ZKigAAFhG1IUBAFagVZMUaq3tT/K6fgAAgBVDXRgAYGWaKHFaVZ/JLP04tdYeseARAQDAMqIuDACwMk2UOE2yaWT8xCTnJXngwocDAADLjrowAMAKNFEfp621L44M062130ny5IFjAwCAJacuDACwMk3aVP9xI5P3SXfXfc0gEQEAwDKiLgwAsDJN2lT/t0bG9yW5MclzFjwaAABYftSFAQBWoIkSp6217x06EAAAWI7UhQEAVqZJm+r/0sGWt9Z++yDrHpfk+iTTrbUfrKqHJ3lHug71P5rkx1trd04eMgAALJ4jqQsDAHD0mujlUOn6cfqPSab64WeSbEzXt9Oh+nf6xSQ3jEy/KsmrW2tnJLk1yfnzCRgAABbZkdSFAQA4Sk3ax+mDkzyutXZ7klTVpUne2Vp70cFWqqr1SX4gySuT/FJVVbo3kD6/L7IlyaVJXjvvyAEAYHEcVl0YAICj26RPnG5IMtqc/s4kp0+w3u8keWmS/f30g5Lc1lrb10/vSHfX/l6q6oKqur6qrt+9e/eEYQIAwII73LowAABHsUmfOH1Lkq1V9Z4kLckPJXnzwVaoqh9Mcktr7SNVdfaB2bMUbbOt31q7IskVSbJp06ZZywAAwCKYd10YAICj30SJ09baK6vqz5M8qZ/1gtbaxw6x2hOTPLOqnp7kxCT3T/cE6ilVtap/6nR9kp2HFzoAAAzvMOvCAAAc5SZtqp8kJyXZ01p7TZIdVfXwgxVurf1ya219a+30JM9N8lettR9N8qEkz+6LbU5y9fzDBgCARTWvujAAAEe/iRKnVfXyJBcl+eV+1vFJ3nqY+7wo3Yui/jVdn6dvOMztAADA4Ba4LgwAwFFi0j5OfyjJY5N8NElaazuras2kO2mtXZvk2n7800nOmleUAACwdI6oLgwAwNFp0qb6d7bWWvoXOVXV6uFCAgCAZWXedeGqOrGqtlbVP1TVP1fVr/XzH15V11XVJ6vqyqo6YeDYAQA4TJMmTq+qqj9I92Knn0ryl0leN1xYAACwbBxOXfirSZ7cWntMkjOTfH9VPSHJq5K8urV2RpJbk5w/YNwAAByBiZrqt9Z+s6qemmRPkkcm+dXW2gcHjQw46tyxdybr1p48cfn1U1PZum37gBEBwJE7nLpw/4TqTD95fD+0JE9O8vx+/pYklyZ57QBhAwBwhA6ZOK2q45L8z9baU5JIlgJz2rd/f3ZeNnlPHusumR4wGgA4ckdSF+7X/UiSb07yP5J8KsltrbV9fZEdSabmWPeCJBckyYYNGw4veAAAjsghm+q31u5K8uWqmvwxMgAAOAYcSV24tXZXa+3MJOvTvRz1UbMVm2PdK1prm1prm9auXTvfXQMAsAAmaqqf5CtJPl5VH0yy98DM1tovDBIVAAAsH0dUF26t3VZV1yZ5Qrp+Ulf1T52uT7JzgHgBAFgAkyZO/7QfAABgpZl3Xbiq1ib5Wp80vW+Sp6R7MdSHkjw7yTuSbE5y9QLHCgDAAjlo4rSqNrTWbmqtbVmsgAAAYDk4wrrwaUm29P2c3ifJVa2191fV9iTvqKpXJPlYkjcsYMgAACygQz1x+t4kj0uSqvrj1toPDx8SAAAsC4ddF26t/WOSx84y/9Pp+jsFAGCZO9TLoWpk/BFDBgIAAMuMujAAwAp2qMRpm2McAACOderCAAAr2KGa6j+mqvaku9t+3348/XRrrd1/0OgAAGDpqAsDAKxgB02cttaOW6xAAABgOVEXBgBY2Q71xCkAAABwCGeduTE7pqcnLr93ZibJ6uECgmXmjr0zWbf25InKrp+aytZt2weOCA5N4hQAAACO0I7p6ey8bPJE6JoL9xy6EBxD9u3fP/FnZN0lk9+EgCEd6uVQAAAAAAArjsQpAAAAAMAYiVMAAAAAgDESpwAAAAAAY7wcagXxlkcAAAAAmIzE6QriLY8AAAAAMBlN9QEAAAAAxkicAgAAAACMkTgFAAAAABgjcQoAAAAAMEbiFAAAAABgjMQpAAAAAMAYiVMAAAAAgDESpwAAAAAAYyROAQAAAADGSJwCAAAAAIyROAUAAAAAGDNY4rSqTqyqrVX1D1X1z1X1a/38h1fVdVX1yaq6sqpOGCoGAAAAAIDDsWrAbX81yZNbazNVdXySD1fVnyf5pSSvbq29o6p+P8n5SV47YBwAAAAwb2eduTE7pqcnKrt3ZibJ6mEDghXijr0zWbf25InLr5+aytZt2weMiJVqsMRpa60lmeknj++HluTJSZ7fz9+S5NJInAIAALDM7Jiezs7LJkuGrrlwz8DRwMqxb//+iT97SbLukslucMB8DdrHaVUdV1XbktyS5INJPpXkttbavr7IjiRTc6x7QVVdX1XX7969e8gwAQAAAADuYdDEaWvtrtbamUnWJzkryaNmKzbHule01ja11jatXbt2yDABAAAAAO5h0MTpAa2125Jcm+QJSU6pqgNdBKxPsnMxYgAAAAAAmNRgidOqWltVp/Tj903ylCQ3JPlQkmf3xTYnuXqoGAAAAAAADsdgL4dKclqSLVV1XLoE7VWttfdX1fYk76iqVyT5WJI3DBgDAAAAAMC8DZY4ba39Y5LHzjL/0+n6OwUAAAAAWJYWpY9TAAAAAICjicQpAAAAAMAYiVMAAAAAgDESpwAAAAAAYyROAQAAAADGSJwCAAAAAIyROAUAAAAAGCNxCgAAAAAwRuIUAAAAAGCMxCkAAAAAwJhVSx0AwKTOOnNjdkxPT1x+/dRUtm7bPmBEADC7qnpokjcn+cYk+5Nc0Vp7TVU9MMmVSU5PcmOS57TWbl2qOAEAmJvEKXDU2DE9nZ2XrZ64/LpLJk+yAsAC25fkJa21j1bVmiQfqaoPJvnJJNe01i6vqouTXJzkoiWMEwCAOWiqDwAAC6y1tqu19tF+/PYkNySZSnJuki19sS1JnrU0EQIAcCgSpwAAMKCqOj3JY5Ncl+TU1tqupEuuJnnIHOtcUFXXV9X1u3fvXqxQAQAYIXEKAAADqar7JfnjJBe21vZMul5r7YrW2qbW2qa1a9cOFyAAAHOSOAUAgAFU1fHpkqZva629u599c1Wd1i8/LcktSxUfAAAHJ3EKAAALrKoqyRuS3NBa++2RRe9Lsrkf35zk6sWODQCAyaxa6gAAAOAY9MQkP57k41W1rZ/3siSXJ7mqqs5PclOS85YoPgAADkHiFFgyd+ydybq1J09cfu/MTJLVwwUEAAuktfbhJDXH4nMWMxYAAA6PxCmwZPbt35+dl02eCF1z4cTv1AAAAAA4Ivo4BQAAAAAYI3EKAAAAADBG4hQAAAAAYIw+TgEAAAA4as33xcPrp6ayddv2ASPiWCFxCgAAAMBRa74vHl53yfSA0XAs0VQfAAAAAGCMxCkAAAAAwBiJUwAAAACAMfo4BQAAAGDF8DIpJiVxCgAAAMCK4WVSTEpTfQAAAACAMRKnAAAAAABjBkucVtVDq+pDVXVDVf1zVf1iP/+BVfXBqvpk//MBQ8UAAAAAAHA4hnzidF+Sl7TWHpXkCUl+rqo2Jrk4yTWttTOSXNNPAwAAAAAsG4O9HKq1tivJrn789qq6IclUknOTnN0X25Lk2iQXDRUHAAAAJMlZZ27MjunJX/Kyd2YmyeQvkAHg2DJY4nRUVZ2e5LFJrktyap9UTWttV1U9ZI51LkhyQZJs2LBhMcJcFub7Rf61O7+S4084caKyvvQBAICVbMf09LzepL3mwj0DRgPAcjd44rSq7pfkj5Nc2FrbU1UTrddauyLJFUmyadOmNlyEy8vhfJHvftWDJi4LAAAAABzakH2cpqqOT5c0fVtr7d397Jur6rR++WlJbhkyBgAAAACA+RoscVrdo6VvSHJDa+23Rxa9L8nmfnxzkquHigEAAAAA4HAM2VT/iUl+PMnHq2pbP+9lSS5PclVVnZ/kpiTnDRgDsILdsXcm69aePHH59VNT2bpt+4ARAQBwMPN954P6GwBDGixx2lr7cJK5OjQ9Z6j9Ahywb//+efUZvO6SySvpAAAsvPm+80H9DYAhDdrHKQAAAADA0UjiFAAAAABgjMQpAAAAAMCYIV8OBQAAAABHNS8eXrkkTgEAAABgDl48vHJpqg8AAAAAMEbiFAAAAABgjMQpAAAAAMAYfZwCAAAwiLPO3Jgd05P39bd3ZibJ5P0IzveFLfPdPgArm8QpAAAAg9gxPT2vF6qsuXDPvLY/3xe2zHf7AKxsmuoDAAAAAIyROAUAAAAAGCNxCgAAAAAwRuIUAAAAAGCMxCkAAAAAwBiJUwAAAACAMRKnAAAAAABjJE4BAAAAAMasWuoAAAAAAOBYccfemaxbe/JEZddPTWXrtu0DR8ThkjgFAAAAgAWyb//+7Lxs9URl110yPXA0HAlN9QEAYIFV1R9W1S1V9U8j8x5YVR+sqk/2Px+wlDECAHBwEqcAALDw3pTk+8fmXZzkmtbaGUmu6acBAFimJE4BAGCBtdb+JsmXxmafm2RLP74lybMWNSgAAOZFH6cAALA4Tm2t7UqS1tquqnrIXAWr6oIkFyTJhg0bFik8AGC5O+vMjdkxPXm/qF4+dWQkTgEAYJlprV2R5Iok2bRpU1vicACAZWLH9PTEL55KvHzqSGmqDwAAi+PmqjotSfqftyxxPAAAHIQnTgc230eo987MJJn8zgEAAEeN9yXZnOTy/ufVSxsOAAAHI3E6sPk+Qr3mwj0DRgMAwGKoqrcnOTvJg6tqR5KXp0uYXlVV5ye5Kcl5SxchAACHInEKAAALrLX2vDkWnbOogQAAcNj0cQoAAAAAMMYTpwC9O/bOZN3akycqu35qKlu3bR84IgAAAGCpDJY4rao/TPKDSW5prT26n/fAJFcmOT3JjUme01q7dagYAOZj3/79E/dJvO6SyV/6BgAAABx9hmyq/6Yk3z827+Ik17TWzkhyTT8NAAAAALCsDJY4ba39TZIvjc0+N8mWfnxLkmcNtX8AAAAAgMO12H2cntpa25UkrbVdVfWQuQpW1QVJLkiSDRs2LFJ4AAAAALA45vOujSTZOzOTZLIu5jhyy/blUK21K5JckSSbNm1qSxwOAAAAACyo+bxrI0nWXLhnwGgYN2Qfp7O5uapOS5L+5y2LvH8AAAAAgENa7MTp+5Js7sc3J7l6kfcPAAAAAHBIgzXVr6q3Jzk7yYOrakeSlye5PMlVVXV+kpuSnDfU/ody1pkbs2N6euLy+p6AY9N8+6FZPzWVrdu2DxgRAAAAsJAGS5y21p43x6JzhtrnYtgxPa3vCWDe/dCsu2TyGy4AAADA0lvspvoAAAAAAMuexCkAAAAAwJjBmuofLfRZCgAAAACMW/GJU32WAgAAAADjVnziFAAAYCWbbyu89VNT2bpt+4ARAcDyIHEKAACwgs23Fd66SyZPsgLA0czLoQAAAAAAxkicAgAA/P/t3Xe0ZVV9wPHvj0EY6UgLSFUpC5UFgiBSHAELxoKBpaDSokE0hGCCxpqAGheG2LFhoVlAJBhiA1FAbPQpDAiCjpGBKKggM4K0X/44+8J5l1duebfO97PWWe/cU/bdZ+997v29fc/ZR5IkqYkdp5IkSZIkSZLUxI5TSZIkSZIkSWpix6kkSZIkSZIkNbHjVJIkSZIkSZKa2HEqSZIkSZIkSU3sOJUkSZIkSZKkJisPOgOStCK4b/kyNtlg7Za33/TJT+bK+Tf0MEeSJEmSJGk6dpxKUh889Mgj3P6+1VvefpP3LO1hbiRJkiRJ0ky8VV+SJEmSJEmSmthxKkmSJEmSJElN7DiVJEmSJEmSpCZjOcbprjtuz21LWxsfcPmyZUDr4w5K0jBq53PPB09JkiRJkibj/5YTjWXH6W1Ll7b8EJY1j/tTj3MjSb3XzueeD56SJEmSJE3G/y0n8lZ9SZIkSZIkSWpix6kkSZIkSZIkNRnLW/UladTdt3wZm2ywdsvbtzNec7tprwjj1kiSJEmS1MyOU0kaQg898kjL48pAe+M1t5v2ijBujSRJkiRJzew4lSRJkqQh1s4TjgEefOB+nrDK3Ja3b+fOFWjv7pV205YkaZjYcSpJkiRJQ6ydJxxDdSfKnR9cr63t29HO3Svtpi1J0jCx41SSNC3HRJUkSZIkrYjsOJUkTcsxUSVJkiRJK6KVBp0BSZIkSZIkSRo2XnEqSZIkSV3q5QOcfMCSJKlT7Q695nfORCPRcXrjDYutZEkaEY6JKklaEfXyAU4+YEmS1Kl2h17zO2eigXScRsSLgY8Bc4DPZ+ZJ023/4IMPcPv71m05fStZkgbHMVElaXrtxsKSJEkajL6PcRoRc4BPAvsD2wOHRMT2/c6HJEmS1G/GwpIkSaNjEA+H2hW4JTN/mZkPAGcDrxhAPiRJkqR+MxaWJEkaEZGZ/X3DiIOAF2fmG8rrQ4HdMvOYpu2OAo4qL58BXN/XjI6e9YG7Bp2JEWA5tcZympll1BrLqTWW08wso9Zsm5lrDjoTmpqx8Ejzc2g4WA/DwXoYDtbD8LAuhsOsx8KDGOM0Jln2uN7bzDwVOBUgIq7OzF16nbFRZhm1xnJqjeU0M8uoNZZTayynmVlGrYmIqwedB83IWHhEWQ/DwXoYDtbDcLAehod1MRx6EQsP4lb924DNaq83BW4fQD4kSZKkfjMWliRJGhGD6Di9Ctg6IraKiFWAg4ELBpAPSZIkqd+MhSVJkkZE32/Vz8yHIuIY4EJgDvDFzFw8w26n9j5nI88yao3l1BrLaWaWUWssp9ZYTjOzjFpjOQ05Y+GRZj0MB+thOFgPw8F6GB7WxXCY9Xro+8OhJEmSJEmSJGnYDeJWfUmSJEmSJEkaanacSpIkSZIkSVKTvnecRsSLI+KmiLglIt4+yfq9I+LaiHgoIg5qWvfdiLg7Ir7ZtHyriLgiIn4REeeUgfZHWo/K6fSI+FVEzC/Tjr0+jl7qtIwiYseI+GlELI6IhRHx6to629Jjy6crp7FqS9BVOW0REdeUclgcEUfX1u0cEYtKmh+PiOjX8fRCj8ro0pJmoy1t2K/j6ZVuPr/L+rUiYmlEnFJbNlZtCXpWTmPVnrqMBR6ulcMFteVj9z03alqo11VL3dxS6mrLsny9iLgkIpbV231ZN1Ztvx+6qIcXlO+0ReXvPrV9xu6zutd6VA+eDx3ooi52rZX1goh4Zatp6vF6VA9LyrkyPyKu7t/RjK5O66G2fvPyfX18q2nq8XpUD+2fD5nZt4lqAPxbgacAqwALgO2bttkS2AE4Ezioad2+wMuAbzYt/xpwcJn/DPCmfh7XCJXT6c3bjurUTRkB2wBbl/lNgDuAdWxLbZXT2LSlWSinVYBVy/wawBJgk/L6SmB3IIDvAPsP+liHsIwuBXYZ9PENQznV1n8M+ApwSm3Z2LSlHpfT2LSnbssIWDZFumP1PTdqU4v1+mbgM2X+YOCcMr86sCdwdL3dl3Vj0/ZHoB52qn2HPQNYWttnrD6rR7gePB/6WxerASuX+Y2B31E9hHrGNJ16Xw/l9RJg/UEf36hM3dRDbf15wLnA8a2m6dT7eijL2j4f+n3F6a7ALZn5y8x8ADgbeEV9g8xckpkLgUead87M7wP31peVX3L3Ab5eFp0BHNCDvPfTrJfTGOq4jDLz5sz8RZm/nepLZQPbUmvl1J9s91035fRAZv6lvFyVciV/RGwMrJWZP83qE/pMRrs9zXoZjamuPr8jYmdgI+Ci2rJxa0vQg3IaQ12V0WTG9Htu1MxYr+X1GWX+68C+ERGZuTwzfwTc37/sjq1u6uG6EhcBLAbmlitexvGzutdmvR76kuvx1E1d/DkzHyrL5wKNp0+3kqYm6kU9qH0d1wNARBwA/JLqs6mdNDVRL+qhI/3+5/XJwG9qr28ry7qxHnB37UNiNtIctF6UU8O/R3Xb9UdGPLiYlTKKiF2pfr24FdvSlJrKqWFc2hJ0WU4RsVlELCxpfLAE8k8u6XSU5hDqRRk1nFZulXhP44tuhHVcThGxEvAh4K2TpDlObQl6U04N49Keuv38nhsRV0fEz0rgCOP5PTdqWqnXR7cpdXUPVd3NZFzafj/MVj0cCFxXfhwcx8/qXutFPTR4PrSnq7qIiN0iYjGwCDi6rO/l/7Tjqhf1AFUn6kVRDWtxVA/zPy46roeIWB34F+DEDtLURL2oB+jgfOh3x+lkX1rd/hLSizQHrVfH9A5gO+DZwJOoGtKo6rqMypUBZwFHZuYjs5HmEOpFOcF4tSXospwy8zeZuQPwNODwiNio2zSHUC/KCOC1mflMYK8yHdp1Tgerm3J6M/DtzPxN0/Jxa0vQm3KC8WpP3db75pm5C/Aa4KMR8dRZSFPda6UOOqmncWr7/dB1PUTE04EPAm9sI01N1It6AM+HTnRVF5l5RWY+nep/g3dExNwW09REvagHgD0y81nA/sDfR8Tes5XhMdVNPZwIfCQzl3WQpibqRT1AB+dDvztObwM2q73eFLh9im1bdRewTkSsPItpDlovyonMvCMrfwFOo7r0eVR1VUYRsRbwLeDdmfmzsti21GSKchq3tgSzdM6VqygXUwXpt5V0ukpziPSijMjMpeXvvVTjVa7IbWl34JiIWAL8J3BYRJzE+LUl6E05jVt76uqca1zVnZm/pBrvbyfG83tu1LRSr49uU+pqbeAP0yU6Zm2/H7qqh4jYFDgfOCwzb61tP26f1b3Wi3rwfOjMrHw2ZeaNwHKqcWd78j/tmOtFPdRjgt9RnTOeE9Prph52A/6jxKnHAe+MiGNaTFMT9aIeOjof+t1xehWwdVRPdF2FavDWC2bYZ1plDKFLgMbTZA8H/rurXA7erJcTPHrlYGOMswOA67tNc4A6LqOy/fnAmZl5bmO5bWmiqcqprBuntgTdldOmEfHEMr8usAdwU2beAdwbEc8p5XQYo92eZr2MImLliFi/LH8C8FJW4LaUma/NzM0zc0vgeKpz7+1j2JagB+U0hu2pm3Nu3cYQKqVM9gBuGNPvuVHTSr1eQFU3UNXVD0rdTWoM234/dFwPEbEO1Y/K78jMHzc2HtPP6l6b9XrwfOhYN3WxVeMHuYjYAtiW6uErPfmfdszNej1ExOoRsWZZvjrwQjwnZtJxPWTmXpm5ZYlTPwp8IDNPaTFNTTTr9dDx+ZD9fzLWS4CbqcZKfFdZ9l7g5WX+2VS9xsuB3wOLa/teDtwJ3Fe2eVFZ/hSqp1jeQvXErFX7fVwjUk4/oBrv5HrgS8Aagz7OQZQR8DrgQWB+bdrRttRWOY1VW+qynF4ALKR6yt9C4KhamruUMroVOAWIQR/nMJUR1ROirynLFlM9JX3OoI9zUOXUlMYRTHxa/Fi1pV6U0zi2py7OueeWz+gF5e/ra2mO3ffcqE0t1OvcUje3lLp6Sm3fJVRXUiwrdb/9OLb9Ya4H4N3lnKvHRxuWdWP3WT1q9eD5MJC6OLSU9XzgWuCA6dJ06m89UH3vLyjTYuuht/XQlMYJTHyau+fDgOuh0/Mhys6SJEmSJEmSpKLft+pLkiRJkiRJ0tCz41SSJEmSJEmSmthxKkmSJEmSJElN7DiVJEmSJEmSpCZ2nEqSJEmSJElSEztOJbUtIl4ZERkR2w06LzOJiHkRcU9EXBcRN0bEv/X4/b4dEeuU6c0d7L9xRHyz9vqrEbEwIt4yC3l7Z9Prn3SR1tkRsXW3eZIkSepWRLwrIhaXmGl+ROxWlh8XEavN4vssiYj1u9h/Xj3Oa1reiFdviogfRsRLu3ifoyPisBm2OSAitq+9fm9E7Nfpe9bSMfZ+LC1jb2kM2HEqqROHAD8CDp6NxCJizmykM43LM3MnYBfgdRGxcys7RcTK7b5RZr4kM+8G1gHaDt6AfwI+V97/r4DnZuYOmfmRbvMGTAjeMvO5HaTR8GngbV3sL0mS1LWI2B14KfCszNwB2A/4TVl9HDBrHacd5K2dGPfyzNwpM7cFjgVOiYh9O3nfzPxMZp45w2YHAI92nGbmv2bmxZ283ySMvSvG3tIYsONUUlsiYg1gD+D11DpOI+KciHhJ7fXpEXFgRMyJiJMj4qry6+0by/p5EXFJRHwFWFSWfSMirilXDBxVS+v1EXFzRFwaEZ+LiFPK8g0i4ryS9lURscd0ec/M5cA1wFMjYm5EnBYRi8ov4s8vaR4REedGxP8AF0Xl5Ii4vmz76rLdxuVqu6iKdwAAB81JREFUgPll3V5leeNKhJPK+8wv+58VEa+oHdOXI+Llk2TzQOC7Zf4iYMOSxl7l+D8QEZcB/xgRL4uIK0r+L46IjRp1VDu2haUeTgKeWNL6ctluWfk71THOK+/59Yj4eclzlLxdDuzXYRApSZI0WzYG7srMvwBk5l2ZeXtEHAtsAlwSEZcARMSnI+LqEmue2EigxG8nRsS1JRbarixfLyIuKrHWZ4Go7TNV3Losqqs3rwB2j4gXlzjqR8DftHJAmTkfeC9wTEnzcTFvRKxU8r1O7b1viYiNIuKEiDi+LPu7ss+CksZqEfFc4OXAySU2fGpUsftBZZ99yzEviogvRsSq05XTNMdh7G3sLY2+zHRycnJqeQJeB3yhzP+E6td9gFcCZ5T5Vah+6X8icBTw7rJ8VeBqYCtgHrAc2KqW9pPK3ycC1wPrUQW8S4AnAU+gChpOKdt9BdizzG8O3DhJfucB3yzz65W0ng78M3BaWb4d8L/AXOAI4LZaXg4EvgfMATYq221c9n9X2WYOsGaZXwKsD2wJXF/Lx/OAb5T5tYFfASs35XUr4Jra6+Y0LgU+VXu9LhBl/g3Ah8r8B4GP1rcrf5c1vd+yGY5xHnAPsCnVD20/bZR32e97wM6DbpNOTk5OTk5OK+4ErAHMB24GPgU8r7ZuCbB+7XUjvptT4qodatv9Q5l/M/D5Mv9x4F/L/F8D2UiPSeLW8jqBV5X5uVQx8dZUna5fo8SlTccwr3k5sCMltmWKmBf4GHBkmd8NuLjMnwAcX+bXq6X5/tpxng4cVFt3OnBQLc/blOVnAsdNV05THQvG3sbeTk5jMHnFqaR2HQKcXebPLq8BvgPsU36R3h/4YWbeB7wQOCwi5gNXUAVQjfF5rszMX9XSPjYiFgA/AzYr2+0KXJaZf8jMB4Fza9vvR3Ub03zgAmCtiFhzkjzvFRHXUf2KfFJmLgb2BM4CyMyfA78Gtinbfy8z/1Dm9wS+mpkPZ+ZvgcuAZwNXAUdGxAnAMzPz3ukKLTMvA54WERuWMjsvMx9q2mxj4M7p0gHOqc1vClwYEYuAt1IFpVCVyydr7/3HGdKc6hihqqPbMvMRqn9Ktqzt9zuqjm1JkqSByMxlwM5UP9bfCZwTEUdMsfmrIuJa4DqquGn72rr/Kn+v4bF4Z2/gS+V9vgXUY6rJ4laAh4Hzyvx2wK8y8xeZmY20WhS1+ali3nOAV5dtDmZinNjwjIi4vMSLr+WxeHEq25Y831xen0FVDg2TlVMzY+/pGXtLI8TLvCW1LCLWA/ahCsCS6lfSjIi3Zeb9EXEp8CKqAO6rjd2ofpm+sCmteVRXnNZf7wfsnpl/LmnNZWLQ2Gylsv19M2T98sxsHmB/unSX1+Yn3S4zfxgRe1NdfXBWRJycM48ldRZVwHow8LeTrL+P6pinU8/bJ4APZ+YFpfxOqOU5Z0inbrqy+Ett/mEmfm/MpcqzJEnSwGTmw1RXB15aOrUOp7qC8lERsRVwPPDszPxjRJzOxLirEfM0xzuPi6mmiVsB7i/5mXL/Fu0E3FjmJ415I+KnVJ2DG1CNWfr+SdI5HTggMxeUDuV5M7zvdHEhTF1Odcbe0zP2lkaIV5xKasdBwJmZuUVmbpmZm1Hd9rJnWX82cCSwF9DoKL0QeFNEPAEgIraJiNUnSXtt4I8l+NwOeE5ZfiXwvIhYt4zpc2Btn4soYz+VtHds41h+SBVIERHbUN32dNMU2706qrFaN6D6xf3KiNgC+F1mfg74AvCspv3uBZqvfj2d6iEFlF/em93M1L/cT2ZtYGmZP7y2vLlc1i2zDzbqocmkx9jC+28DTHYckiRJfRER28bEp43vSHU1I0yMx9ai6gS7p4xNuX8Lydfjxf2pbtWGqePWZj8HtoqIp5bXh0yx3QQRsQPwHh67inHSmLdcxXo+8GGq2/d/P0lyawJ3lBjwtbXlk8WqjTxvGRFPK68PpboislvG3o8x9pZGiB2nktpxCFVwVnce8JoyfxHVF//FmflAWfZ54Abg2oi4Hvgsk/86/V1g5YhYCLyP6rYnMnMp8AGq2/wvLmndU/Y5FtilDMJ+A3B0G8fyKWBOuSrhHOCILA8VaHI+sBBYAPwAeFtm/h/Vr/Xzy21IB1KNMfWoErj+uAz6fnJZ9luqKwdOmyxDWQ2gf2stUJ3JCcC5EXE5cFdt+fuBdct7LwCeX5afCixsDFDfwjFOqfzDcV9m3tFiXiVJknphDeCMiLihxJHb89iVgKcC34mISzJzAdUt+ouBLwI/biHtE4G9y+39L6QaixKmiFubZeb9VEMIfCuqh0P9erLtir3KQ4duouowPTYzv1/WTRfznkP1DILJbtOHqgP2CqrxMX9eW3428Nbyno2O3Uaej6SKMRcBjwCfmSbfrTL2nvkYp2TsLQ1OY2BjSRpaEbFGZi4rV5yeD3wxM5s7cIdeRKwGLKJ6oNY9U2zzSqpB39/d18y1KSLeAvwpM78w6LxIkiRJzYy9Jc0GrziVNApOKIPhX081NMA3BpyftkXEflS/8n9iqsANoHQIL+lXvrpwN9XDAiRJkqShYuwtabZ4xakkSZIkSZIkNfGKU0mSJEmSJElqYsepJEmSJEmSJDWx41SSJEmSJEmSmthxKkmSJEmSJElN7DiVJEmSJEmSpCb/D4qq9KdgAJExAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Summary Statistics for Bootstrap Porosity Mean Realizations:\n",
      "DescribeResult(nobs=1000, minmax=(0.10831378589655172, 0.1348213317413793), mean=0.12148225985881035, variance=1.856591604516782e-05, skewness=0.14690316583288418, kurtosis=-0.06206387559881854)\n",
      "P10: 0.116, P50: 0.121, P90: 0.127\n",
      "\n",
      "Summary Statistics for Bootstrap Porosity Standard Deviation Realizations:\n",
      "DescribeResult(nobs=1000, minmax=(0.02066834194552453, 0.042905367804465244), mean=0.03251750740895691, variance=1.6943900665148023e-05, skewness=-0.0069496954114074095, kurtosis=-0.37886276009681374)\n",
      "P10: 0.027, P50: 0.032, P90: 0.038\n"
     ]
    }
   ],
   "source": [
    "L = 1000                                   # set the number of realizations\n",
    "mean = np.zeros(L); stdev = np.zeros(L)    # declare arrays to hold the realizations of the statistics\n",
    "for l in range(0, L):                      # loop over realizations\n",
    "    samples = random.choices(df['Porosity'].values, weights=df['Wts'].values, cum_weights=None, k=len(df))\n",
    "    mean[l] = np.average(samples)\n",
    "    stdev[l] = np.std(samples)\n",
    "    \n",
    "plt.subplot(121)\n",
    "GSLIB.hist_st(mean,0.11,0.15,False,False,50,None,'Average Porosity (fraction)','Bootstrap Uncertainty in Porosity Average')\n",
    "\n",
    "plt.subplot(122)\n",
    "GSLIB.hist_st(stdev,0.015,0.045,False,False,50,None,'Standard Deviation Porosity (fraction)','Bootstrap Uncertainty in Porosity Standard Deviation')\n",
    "\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=3.0, top=1.2, wspace=0.2, hspace=0.2)\n",
    "plt.show()   \n",
    "    \n",
    "print('Summary Statistics for Bootstrap Porosity Mean Realizations:')\n",
    "print(stats.describe(mean))\n",
    "print('P10: ' + str(round(np.percentile(mean,10),3)) + ', P50: ' + str(round(np.percentile(mean,50),3)) + ', P90: ' + str(round(np.percentile(mean,90),3))) \n",
    "\n",
    "print('\\nSummary Statistics for Bootstrap Porosity Standard Deviation Realizations:')\n",
    "print(stats.describe(stdev))\n",
    "print('P10: ' + str(round(np.percentile(stdev,10),3)) + ', P50: ' + str(round(np.percentile(stdev,50),3)) + ', P90: ' + str(round(np.percentile(stdev,90),3))) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comments\n",
    "\n",
    "This was a basic demonstration of bootstrap. Much more could be done, you could replace the statistics, average and standard deviation with any other statistics, for example P90, kurtosis, P13 etc. I have other demonstrations on the basics of working with DataFrames, ndarrays, univariate statistics, plotting data, declustering, data transformations, trend modeling and many other workflows available at https://github.com/GeostatsGuy/PythonNumericalDemos and https://github.com/GeostatsGuy/GeostatsPy. \n",
    "  \n",
    "I hope this was helpful,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "Michael Pyrcz, Ph.D., P.Eng. Associate Professor The Hildebrand Department of Petroleum and Geosystems Engineering, Bureau of Economic Geology, The Jackson School of Geosciences, The University of Texas at Austin\n",
    "\n",
    "#### More Resources Available at: [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}