Interactive_Declustering.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
    "## Interactive Spatial Data Declustering Demonstration\n",
    "\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)\n",
    "\n",
    "\n",
    "### The Interactive Workflow\n",
    "\n",
    "Here's a simple workflow for calculating and evaluating spatial data declustering weights smapling bias due to clustered sampled mitigation.\n",
    "\n",
    "\n",
    "### Basic Univariate Summary Statistics and Data Distribution Representativity Plotting in Python with GeostatsPy\n",
    "\n",
    "Here's a simple workflow with some basic univariate statistics and distribution representativity. This should help you get started data declustering to address spatial sampling bias.\n",
    "\n",
    "#### Geostatistical Sampling Representativity\n",
    "\n",
    "In general, we should assume that all spatial data that we work with is biased.\n",
    "\n",
    "##### Source of Spatial Sampling Bias\n",
    "\n",
    "Data is collected to answer questions:\n",
    "* how far does the contaminant plume extend? – sample peripheries\n",
    "* where is the fault? – drill based on seismic interpretation\n",
    "* what is the highest mineral grade? – sample the best part\n",
    "* who far does the reservoir extend? – offset drilling\n",
    "and to maximize NPV directly:\n",
    "* maximize production rates\n",
    "\n",
    "**Random Sampling**: when every item in the population has a equal chance of being chosen. Selection of every item is independent of every other selection.\n",
    "Is random sampling sufficient for subsurface?  Is it available?\n",
    "* it is not usually available, would not be economic\n",
    "* data is collected answer questions\n",
    "    * how large is the reservoir, what is the thickest part of the reservoir \n",
    "* and wells are located to maximize future production\n",
    "    * dual purpose appraisal and injection / production wells!\n",
    "\n",
    "**Regular Sampling**: when samples are taken at regular intervals (equally spaced).  \n",
    "* less reliable than random sampling.\n",
    "* Warning: may resonate with some unsuspected environmental variable.\n",
    "\n",
    "What do we have?\n",
    "* we usually have biased, opportunity sampling \n",
    "* we must account for bias (debiasing will be discussed later)\n",
    "\n",
    "So if we were designing sampling for representativity of the sample set and resulting sample statistics, by theory we have 2 options, random sampling and regular sampling.\n",
    "\n",
    "* What would happen if you proposed random sampling in the Gulf of Mexico at $150M per well?\n",
    "\n",
    "We should not change current sampling methods as they result in best economics, we should address sampling bias in the data.\n",
    "\n",
    "Never use raw spatial data without access sampling bias / correcting.\n",
    "\n",
    "##### Mitigating Sampling Bias\n",
    "\n",
    "In this demonstration we will take a biased spatial sample data set and apply declustering using **GeostatsPy** functionality.\n",
    "\n",
    "#### Objective \n",
    "\n",
    "In the PGE 383: Stochastic Subsurface Modeling class 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": 105,
   "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",
    "from ipywidgets import interactive        # widgets and interactivity\n",
    "from ipywidgets import widgets                            \n",
    "from ipywidgets import Layout\n",
    "from ipywidgets import Label\n",
    "from ipywidgets import VBox, HBox\n",
    "from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) # control of axes ticks\n",
    "plt.rc('axes', axisbelow=True)            # set axes and grids in the background for all plots\n",
    "from matplotlib.patches import Rectangle  # drawing shapes on plots\n",
    "from statsmodels.stats.weightstats import DescrStatsW"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you get a package import error, you may have to first install some of these packages. This can usually be accomplished by opening up a command window on Windows and then typing 'python -m pip install [package-name]'. More assistance is available with the respective package docs.  \n",
    "\n",
    "#### Declare Functions\n",
    "\n",
    "These functions read in the multiple realizations and produce local statistical summaries that we will cover below.  They will shortly be added to GeostatsPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [],
   "source": [
    "def declus(df, xcol, ycol, noff, csize):\n",
    "    # modified from GeostatsPy Python Package to run 1 cell size and report all cell offsets\n",
    "    # Load data and set up arrays\n",
    "    nd = len(df)\n",
    "    x = df[xcol].values\n",
    "    y = df[ycol].values\n",
    "    wt = np.zeros([nd,noff])\n",
    "    wtopt = np.ones(nd)\n",
    "    index = np.zeros(nd, np.int32)\n",
    "    roff = float(noff)\n",
    "\n",
    "    # Calculate data extents\n",
    "    xmin = np.min(x); xmax = np.max(x); ymin = np.min(y); ymax = np.max(y)\n",
    "    xmin = 0.0; ymin = 0.0; xmax = 1000.0; ymax = 1000.0\n",
    "    # Define a \"lower\" origin to use for the cell sizes\n",
    "    xo1 = xmin\n",
    "    yo1 = ymin\n",
    "\n",
    "    # Only 1 cell size\n",
    "    ncellx = int((xmax - (xo1 - csize)) / csize) + 1\n",
    "    ncelly = int((ymax - (yo1 - csize)) / csize) + 1\n",
    "    ncellt = ncellx * ncelly\n",
    "    cellwt = np.zeros(ncellt)\n",
    "    \n",
    "    # Loop over all the origin offsets selected\n",
    "    xfac = min((csize / roff), (0.5 * (xmax - xmin)))\n",
    "    yfac = min((csize / roff), (0.5 * (ymax - ymin)))\n",
    "    for kp in range(0, noff ):\n",
    "        xo = xo1 - (float(kp)) * xfac\n",
    "        yo = yo1 - (float(kp)) * yfac\n",
    "\n",
    "        # Initialize the cumulative weight indicators\n",
    "        cellwt.fill(0.0)\n",
    "\n",
    "        # Determine which cell each datum is in\n",
    "        for i in range(0, nd):\n",
    "            icellx = int((x[i] - xo) / csize) + 1\n",
    "            icelly = int((y[i] - yo) / csize) + 1\n",
    "            icell = icellx + (icelly - 1) * ncellx - 1\n",
    "            index[i] = icell\n",
    "            cellwt[icell] = cellwt[icell] + 1.0\n",
    "\n",
    "        # The weight assigned to each datum is inversely proportional to the\n",
    "        # number of data in the cell. We first need to get the sum of\n",
    "        # weights so that we can normalize the weights to sum to one\n",
    "        sumw = 0.0\n",
    "        for i in range(0, nd):\n",
    "            ipoint = index[i]\n",
    "            sumw = sumw + (1.0 / cellwt[ipoint])\n",
    "        sumw = 1.0 / sumw\n",
    "\n",
    "        # Accumulate the array of weights (that now sum to one)\n",
    "        for i in range(0, nd):\n",
    "            ipoint = index[i]\n",
    "            wt[i,kp] = (1.0 / cellwt[ipoint]) * sumw\n",
    "\n",
    "    # End loop over all offsets\n",
    "\n",
    "    # Correct weights over all offsets\n",
    "    for ioff in range(0, noff):\n",
    "        sumw = np.sum(wt[:,ioff])\n",
    "        wt[:,ioff] = wt[:,ioff]/sumw\n",
    "                        \n",
    "    return wt\n",
    "\n",
    "def weighted_percentile(data,perc,weights=None): # from Luca Jokull (https://stackoverflow.com/questions/21844024/weighted-percentile-using-numpy)\n",
    "    if weights is None:\n",
    "        return nanpercentile(data,perc)\n",
    "    else:\n",
    "        d=data[(~np.isnan(data))&(~np.isnan(weights))]\n",
    "        ix=np.argsort(d)\n",
    "        d=d[ix]\n",
    "        wei=weights[ix]\n",
    "        wei_cum=100.*np.cumsum(wei*1./sum(wei))\n",
    "        return np.interp(perc,wei_cum,d)"
   ]
  },
  {
   "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": 52,
   "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": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('sample_data_biased.csv')     # load our data\n",
    "df['Porosity'] = df['Porosity'] * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We loaded our file into our DataFrame called 'df'. But how do you really know that it worked? Visualizing the DataFrame would be useful and we already leard about these methods in this demo (https://git.io/fNgRW). \n",
    "\n",
    "#### Resample from the Data\n",
    "\n",
    "Option to reduce the number of spatial sample data to explore the impact on spatial uncertainty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "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>index</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>0</th>\n",
       "      <td>14</td>\n",
       "      <td>400</td>\n",
       "      <td>900</td>\n",
       "      <td>1</td>\n",
       "      <td>22.366071</td>\n",
       "      <td>2372.383732</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>246</td>\n",
       "      <td>730</td>\n",
       "      <td>169</td>\n",
       "      <td>1</td>\n",
       "      <td>13.372385</td>\n",
       "      <td>7.350097</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>221</td>\n",
       "      <td>490</td>\n",
       "      <td>959</td>\n",
       "      <td>1</td>\n",
       "      <td>19.656596</td>\n",
       "      <td>403.008633</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>136</td>\n",
       "      <td>0</td>\n",
       "      <td>309</td>\n",
       "      <td>0</td>\n",
       "      <td>8.777591</td>\n",
       "      <td>0.665880</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>31</td>\n",
       "      <td>800</td>\n",
       "      <td>300</td>\n",
       "      <td>1</td>\n",
       "      <td>14.708393</td>\n",
       "      <td>25.729277</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   index    X    Y  Facies   Porosity         Perm\n",
       "0     14  400  900       1  22.366071  2372.383732\n",
       "1    246  730  169       1  13.372385     7.350097\n",
       "2    221  490  959       1  19.656596   403.008633\n",
       "3    136    0  309       0   8.777591     0.665880\n",
       "4     31  800  300       1  14.708393    25.729277"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seed = 73061\n",
    "ndata = 50\n",
    "df = df.sample(ndata,replace = False,random_state = seed) # extract 50 samples\n",
    "df = df.reset_index()                                     # reset the record index\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Summary Statistics for Tabular Data\n",
    "\n",
    "The table includes X and Y coordinates (meters), Facies 1 and 2 (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>index</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>count</th>\n",
       "      <td>50.000000</td>\n",
       "      <td>50.000000</td>\n",
       "      <td>50.000000</td>\n",
       "      <td>50.00000</td>\n",
       "      <td>50.000000</td>\n",
       "      <td>50.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>162.180000</td>\n",
       "      <td>481.000000</td>\n",
       "      <td>562.980000</td>\n",
       "      <td>0.86000</td>\n",
       "      <td>14.050819</td>\n",
       "      <td>223.418381</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>84.369133</td>\n",
       "      <td>242.219752</td>\n",
       "      <td>308.571659</td>\n",
       "      <td>0.35051</td>\n",
       "      <td>3.873516</td>\n",
       "      <td>487.955880</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>14.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>9.000000</td>\n",
       "      <td>0.00000</td>\n",
       "      <td>5.854787</td>\n",
       "      <td>0.075819</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>97.000000</td>\n",
       "      <td>300.000000</td>\n",
       "      <td>302.250000</td>\n",
       "      <td>1.00000</td>\n",
       "      <td>11.088426</td>\n",
       "      <td>4.204787</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>163.500000</td>\n",
       "      <td>455.000000</td>\n",
       "      <td>609.500000</td>\n",
       "      <td>1.00000</td>\n",
       "      <td>13.378642</td>\n",
       "      <td>25.234794</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>239.250000</td>\n",
       "      <td>662.500000</td>\n",
       "      <td>877.750000</td>\n",
       "      <td>1.00000</td>\n",
       "      <td>15.880765</td>\n",
       "      <td>82.676408</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>284.000000</td>\n",
       "      <td>920.000000</td>\n",
       "      <td>999.000000</td>\n",
       "      <td>1.00000</td>\n",
       "      <td>22.366071</td>\n",
       "      <td>2372.383732</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            index           X           Y    Facies   Porosity         Perm\n",
       "count   50.000000   50.000000   50.000000  50.00000  50.000000    50.000000\n",
       "mean   162.180000  481.000000  562.980000   0.86000  14.050819   223.418381\n",
       "std     84.369133  242.219752  308.571659   0.35051   3.873516   487.955880\n",
       "min     14.000000    0.000000    9.000000   0.00000   5.854787     0.075819\n",
       "25%     97.000000  300.000000  302.250000   1.00000  11.088426     4.204787\n",
       "50%    163.500000  455.000000  609.500000   1.00000  13.378642    25.234794\n",
       "75%    239.250000  662.500000  877.750000   1.00000  15.880765    82.676408\n",
       "max    284.000000  920.000000  999.000000   1.00000  22.366071  2372.383732"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Specify the Area of Interest\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",
    "vmin = 0.05; vmax = 25.0;               # range of porosity values\n",
    "cmap = plt.cm.inferno                   # color map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualizing Tabular Data with Location Maps¶ \n",
    "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": "iVBORw0KGgoAAAANSUhEUgAAAf8AAAGDCAYAAAAlERkEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3hUVfrA8e87SSaNhB460kEEAQmgImJDECzrqj9FVgULNuysgt3Fhg17QcSKqKuIDQtWdFdcitgAAelFIZRACKnz/v64FxhCCplMcjOZ9/M89yH3zDl33onIO+eec88RVcUYY4wx0cPndQDGGGOMqVqW/I0xxpgoY8nfGGOMiTKW/I0xxpgoY8nfGGOMiTKW/I0xxpgoY8nfRDUReUlE7nZ/PkZE1nodU6QTkd9E5Biv4zDGlMySv4koIjJWRGYUKVtaQtk5YX5vFZGdIpIlIptF5AsRObsc7Sv9y4WIrBSRXW6Mf4nIiyJSqzLfsyhVPURVv3bjuVNEXqvK9zfGlM2Sv4k0s4C+IhIDICKNgTjgsCJl7dy64dZNVWsBHYGXgCdF5I5KeJ+KOMWN8TCgF3BreS8gIrFhj8oYU21Y8jeRZg5Osu/unh8NfAX8XqTsD1VdDyAinURkpohsEZHfReT/KhqEqmao6qvA5cBYEanvvtcIEVkkIjtEZLmIXOqWJwMfA03dXnmWiDQVkd4i8r2IbBORDSLypIj4KxqfG+M69z27uDE0FZH33d/DMhG5ZHddt4f+toi8JiLbgeEiEi8ij4rIevd4VETi3foNRORDN+4tIvKtiPjc11aKyAkiMgi4GTjb/bw/ichZIjIvOE4RuUFEpofjMxtjDowlfxNRVDUP+AEnweP++S3wXZGyWbAn6c4EXgfSgKHA0yJySJhCeg+IBXq75xuBk4FUYAQwQUQOU9WdwEnAelWt5R7rgULgOqABcARwPHBFOAITkRbAYOBHt2gqsBZoCpwJ3Csixwc1OQ14G6gDTAFuAQ7H+VLVzf2Mu+8i3OBeqyHQCCfJ77NWuKp+AtwLvOl+3m7A+0BrETk4qOo/gFfD8JGNMQfIkr+JRN+wN9H3w0n+3xYp+8b9+WRgpaq+qKoFqjofeAcn+VWYquYDGUA99/wjVf1DHd8An7nxlNR+nqrOdmNbCTwH9K9gWNNFZBvOF6JvcJJ8C+Ao4CZVzVHVBcAk4Lygdt+r6nRVDajqLmAY8C9V3aiqm4C7gurnA02Ag1Q1X1W/1QPYKERVc4E3cRI+7pewVsCHFfzMxphysORvItEs4CgRqQs0VNWlwH+BI92yLuwd7z8I6OPent7mJsVhQONwBCIicTi93y3u+UkiMtu9Fb4Np+fdoJT2Hdzb53+6t9vvLam+iDwbNGRwcylh/U1V66jqQap6hZvImwJbVHVHUL1VQLOg8zVFrtPUrRNcv6n784PAMuAzd3hjTCnxFPUycK6ICM6XibfcLwXGmCpiyd9Eou+B2sBI4D8AqrodWO+WrVfVFW7dNcA3bjLcfdRS1cvDFMtpQAHwP3c8/B3gIaCRqtYBZgDi1i2uZ/wMsBhor6qpOLfPpZh6qOplQUMG95YzzvVAPRFJCSprCawLfoti2hxUpP56N5YdqnqDqrYBTgGuLzKEUNI1UdXZQB7OHZFzsVv+xlQ5S/4m4rg92bnA9Ti3+3f7zi0LnuX/IdBBRM4TkTj36FVkzLncRKSeiAwDngLGq+pmwA/EA5uAAhE5CTgxqNlfQH0RqR1UlgJsB7JEpBPOBMKwU9U1OHdH7hORBBE5FLgIZ2y/JFOBW0WkoYg0AG4HXgMQkZNFpJ3be9+OM3ehsJhr/AW02j0ZMMgrwJNAgap+V5HPZowpP0v+JlJ9gzOBLzhxfOuW7Un+7m3uE4FzcHqtfwLjcZJ0KH4SkSycW94XA9ep6u1B73U18BawFadX+35QLItxEupydwiiKTDarbcDeB5nPLyyDMUZX18PvAvcoaozS6l/N86XrJ+BX4D5bhlAe+BzIAvnTszTu5/tL+Lf7p+bRWR+UPmrOMMz1us3xgNyAHN0jDEmrEQkEefJiMPcORvGmCpkPX9jjBcuB+ZY4jfGG5WW/EVksohsFJFfg8rquYutLHX/rOuWi4g87i488rOIHBbU5gK3/lIRuaCy4jXGVA0RWQlcg7NWgDHVmoi0EJGv3MW7fhORa9zyO0VknYgscI/BJbQfJM7iYsuCn4oRkdYi8oOb296UMC3udaAq7ba/iByNMx74iqruXmHsAZzHje53fwl1VfUm95d2Fc5jUX2Ax1S1j4jUwxlzTMeZNTwP6KmqWyslaGOMMSaIiDQBmqjqfPdpmXnA34D/A7JU9aFS2sYAS4ABOItizQGGqupCEXkLmKaqb4jIs8BPqvpMZX+e3Sqt56+qs3CffQ5yGs4zvrh//i2o/BV3YZTZQB33Fz4QmKmqW9yEPxMYVFkxG2OMMcFUdYO7ONjuSb2L2Hd9jNL0Bpap6nJ3ddI3gNPcp2SOw1lRE/bNh1Wiqsf8G6nqBnB+oTgzs8H5RQYvMLLWLSup3BhjjKlSItIK6IGzxDjAKHeoevLuYewiSsph9YFtqlpQpLzKVJedu4pb1ERLKd//AiIjcRZ4ITk5uWenTp3CF50Jq+VLfqdFrbh9/uMGVFm3q5DW7Tp4FpcxpnLMmzcvQ1UbVsa1Bw7qphkZO8quWIL581b8BuQEFU1U1YlF64mzNfY7wLWqul1EngHG4eSkccDDwIVFmxXzluXKbZWlqpP/XyLSRFU3uLf1N7rla4EWQfWa4zyLvBY4pkj518Vd2P2PNREgPT1d586dG97ITdg8Ov5e6vz0Gad03LvC7pRf1pE66DyGXzzSw8iMMZVBRFaVXSs0GRnZ/PC/J0JuHxdzco6qppdWx13G+x1giqpOA1DVv4Jef57i96coKbdl4Axvx7q9/93lVaaqb/u/D+yesX8Bzo5ou8vPd2f9Hw5kusMCnwInikhd95bKiW6ZiWCXXn0d8+JbcPfsNbz12zru+H4165p0YdgFI7wOzRgTcQpR3RbyURZ3fP4FYJGqPhJU3iSo2unAr0Xb4kzwa+/O7PfjLDb2vrsJ1lfs3WAsOB9WiUrr+YvIVJxeewMRWQvcAdwPvCUiFwGrgbPc6jNwZvovA7JxtkJFVbeIyDicXyA4O4wVnURoqqmPPviAN59/nuwd22nVsSNX3HgTbdq0ITExkSdefJVFixaxatUqTmzfnrZt23odrjEmEqmyd+i8UvTF2YDqFxFZ4JbdDAwVke44t+tXApcCuCt3TlLVwapaICKjcDqtMcBkVf3NvcZNwBsicjfOttsvVOaHKKpGrvBnt/299+Zrr/HtpIlc2KEdtRPi+T1jMy+s+5Mn33iTJk2alH0BY0yNISLzyrq1HqqePVvpf364PeT2iXEXVVps1Zmt8GfCLhAI8Oak57nikE7UTnCW0O/YoD6n16/D1Jde9Dg6Y0zNEgOkVuCITtVltr+pQbKyskgRIT4mZp/yDvXr8f1vCz2KyhhTEykF2Ghw+VnyN2FXq1YtdqiSW1BAfOzev2K/Z2yhXZ8jPIzMGFMTVfKYf41kt/1N2Pl8Ps699DKe+G0xW3btQlX5deMm3tuaybkjij4Ga4wxFaFooCDkI1pZz99UijOHDqV2vXo8/dyzZG1bTbsuh/Dw+Idp1KiR16FVS4FAgDlz5rBp0ya6detGixYtym5kjHFn+xd6HUXEseRvKs2AgQMZMHCg12FUexs3buSfVw6jU92NtKidx70vxtOm5ymMHvsvnEeMjTElkhh8UsfrKCKOJX9jPHbfHdczqu9GDmubDCRzTj/l7mnv8uWX/Tn++BO8Ds+Y6k0LCAQyvI4i4tiYvzEeys7OZvPa39zE7xARzj3Sz2fvT/EwMmMih2pByEe0sp6/MR4r7s6+CNTEBbiMCT+FQL7XQUQc6/kb46GkpCTqNOnEzyt27ilTVd6cnccJJw/1MDJjIoO6y/taz798rOdvjMfG3PkIN145jEMXZtCyTh7f/pFA00OGMGDAiV6HZky1JxKLT+p7HUbEseRvjMeaNGnCS/+eyffff8+mTZu4fmQ32rRp43VYxkQE1XwChX+VXdHsw5K/MdVATEwMRx11lNdhGBN5VCFgz/mXlyV/Y4wxkS2Kx+5DZcnfGGNMxBKJJcbX0OswIo4lf2OMMZFL8wkUbPA6iohjyd8YY0zkUoUo3qAnVJb8jYkweXl5zJkzh0AgQK9evUhISPA6JFOCwsJC5s+fT1ZWFj179iQ1NdXrkGogS/6hsORvTASZN28uD9wxisPb5OIT5Yn74rjm5ofpe9TRXodmili+fDnXXXwlDXMTSCSWBwsyuPD6Kzjz7LO8Dq3GEdvVr9ws+RsTIXJychh/+xU8OSKW+m4P8vydBVx233V0nfqV9SqrEVXln5dfw0nxnWlUrx4A/QIFvPLQs/TsnU7r1q09jrAmicUX09jrICKOLe9rPBMIBNi+fTuBQMDrUCLC7NmzObpDLvVT4/aUpSbHMvCQPGbNmuVhZKaoP/74g8QsaJRcb09ZnC+WwxJb8tH0DzyMrAbSfAL560I+opX1/E2VU1Wmvvoy016aRGqMskNjGXb5KP52ht0OLU1BQQFxvv03+/HHKgUFNuZZneTn5xNTTN8qVmLIz7NNaMJLERvzLzfr+Zsq9/60d/hxyrM8dkQzxh9xEA/3asRnzzzEN19/5XVoxcrLy+Pzzz/n9SmvsmjRIs/i6NOnD18tjiM7Z+/4Zl5+gE9+iaNfv36exWX217FjR7bG5ZCZu3fDpoAGWLBrNQNPOcnDyGqoQGHoR5Synr+pcm+/OIm7urfAHxMDQGJcLFd0bcoTzz1F/2OO9SyuzMxMsrOzady4MeLus7t27Vr+ec15HNkhh+YNAjz/INRqfAR33v0IPl/VfndOSUnhoqvv5tInbubkbvnE+ODDBbGcdeEt1K9vG5tUJz6fj7sm3M8tV9xAB2lIosSxKP9PBgw9lc6dO3sdXs2iaiv8hcCSv6lyudlZpMTX26csrVYCWzd7s1DH9u3bGXfnaDZtWEhqrVg2ZPho3aYTG9YuZfmypTw8qjHpXZ0JRaf0h4en/JdPPpnB4MEnV3msAwaexGHpvfnqyy8IBAp5+KrjSUtLq/I4TNl69OjBWzPf58svvmDH9h1cdnQ/DjroIK/DqnFE4oiJbeZ1GBHHkr+pck1at2VZxmbaNdg7O33++s107t7Tk3huu/kqhhyxlf59DmFHVh6X3TSTjonLueGyzlx3/04ap2xlc0YMtVJS8cf5OfO4Ojz76dueJH+A+vXrc+ZZ/+fJe5vyqVWrFqeedprXYdRoqnkE8lZ7HUbEseRvqtyVN93CrSMvYNhBOXRokMovf2Xy9oZ8Jvzr+iqPZf369RRkr6R/n04AfDBzBX/vl8CA3slk7dyOCAiFbN64jpydm8nNh8ycVGLjrAdnTPVgu/qFwpK/qXIdOnRgwpS3mfriC3y+ZDHtu/XhyUdGeHL7etu2baTV3/vo3O9/ZDB8QBKxsT7iCpVN2wr5bcUujumeSGxcDKpw5SPLOeK0kVUeq1fy8/P58vPP+W3+PJq3bcth6b345vPP2bFtK0efOJAePXrsmSNhTJVTkEp8XFhEWgCvAI2BADBRVR8TkQeBU4A84A9ghKpuO5C27mvjgNPc8o3AcFVdX2kfpAhL/sYTzZo1Y/Stt3sdBu3atWPhslxycgtIiI+ledNUlq7eTJI/nvikJGol+Pj3VzuZ9VMObZsl8v3CPNIa1Gb18p+9Dr1KZGdnc8X5w2iXk0m3Osks+noG91y7hPO7H0Kr2qm89PnHfNTnKG69+x77AmA8IeJH4ppX5lsUADeo6nwRSQHmichMYCYwVlULRGQ8MBa46UDaqupC4EFVvc35DHI1cDtwWWV+kGCW/E1U8/v9nH/RaG64bzznn5ZGl04Nueuh37l1eBN6dE2ibmoMd42oz6K1fgp9KQw+LpXCgDLhwy1eh14lXpv8AocXZHHywa1AIXX7Fjp3a8MHf23irM7t6dWsEQ/N/pZ58+aRnp7udbgmCqnmormrKvH6ugHY4P68Q0QWAc1U9bOgarOBMw+0LbBQVbcHVU0G9l/EoxJZ8jdRb8jJp9G2XUfef/d1duzYynmXjee9OV8x6bM/WPJHgM3ZqRx7eDNwO7Yvffgnh/c774CunZGRwZIlS2jWrFlEzvT+7rNPua1NE8BZXjhehBZ1U3hu6d67k8c0rMO3Mz+z5G88IQpSRWP+ItIK6AH8UOSlC4E3y9tWRO4BzgcygSp9ztmSvzFAp06d6DT2X3sLLrwYgJ9++on7bx/JucdtokWjeL77OZufNzTj8TH7fcnfh6ry2EPj+Onbd+nRMsCbm3zEpXVn3ANPk5iYWJkfJaySatViR14eiXGx+Hw+AiiFAd2ni7Ijv4DkFNtXwHilwhP+GojI3KDziao6sWglEakFvANcG9xrF5FbcG7vTynpDUpqq6q3ALeIyFhgFHBHRT5IedgKf8aUolu3bvzzjmd46791uGXiZtYXHsWDj04ucxvdGR++T9aiN5l4YQJXDEjm4XMT6ZMyh6cevbeKIg+P08+7gClL1hJQxR/vJxATyzsrN5DewnmuemdePjM2bmPI6ad7HKmJZhIoDPkAMlQ1PegoLvHH4STvKao6Laj8AuBkYJiqFnvbvqS2RbwOnFGx30L5WM/fmFLMmfM/Hr/nckYcK7RoWIvvFn7JtZcv4qlJ/yYpKanEdh++8yK3HZu0zyS403qlMOy5j2HsuKoIPSxOHDSIlUuXMPrtN+mYksgf23NYviNA95R4nl60gqW78hl15900a2aLrBiPiB+Jb1l5l3f+J34BWKSqjwSVD8KZ4NdfVbPL09Z9rb2qLnVPTwUWV0b8JbHkb0wJVJUnH7yFB4cn07COH4A2TROJ+2oD777zJsPOG1Fi25xdu0hO2PfGms8nCJH1PLKIcOnV1zB0+AhWrlxJ48aNSUtLY8WKFWRlZdGpUyfi4uLKvpAxlSWQi+5aUZnv0Bc4D/hFRBa4ZTcDjwPxwEz3S/5sVb1MRJoCk1R1cEltVXUGcL+IdMR51G8VVTjTHyz5G1OiHTt2EM82GtapvU/5sd2SeWjmZ6Um//4nns70OU8wrN/etvOWZdGmU+9Ki7cypaamcuihh+45t/3oTXVSmc/5q+p37Jnuu48ZJdRfDwwuoy2qWqW3+Yuy5G9MCRITE9meDYGA4vPt/f93w+Y86jVsUmrbc4YN57orvmDdB7/Tu1Uhyzb6mLW6Ho88e3dlh21MdFFb4S8UlvyNKUFcXBzpR53E619+yLDj6yIiZO0q4LlP87juXxeX2jYhIYGnJr3J999/z8Kf59KyR1tePGEA8fHxVRS9MVFEK6/nX1NZ8jemFFdffyuPjC/gggmfkVbbx4btcVw86r4D2pbV5/PRt29f+vbtWwWRGhOlfH58Ca0qcIGij+xHB0v+xpTC7/cz5rZ7yc6+lczMTBo1aoTPZ0/IGlNtBHLR7OVeRxFxLPkbcwCSkpJKfbRvt/Xr1zPzk49QDXDcgJNo2bLyHkGKBEuWLOGbL74gMTGRgUOG0LBhQ69DqtFUlR9//JEfvvsPDRs3ZuBJg0hJSfE6rMqnNuZfXtaFMSZM3p/+NmMuGUTtPx6h/srHuHPUybwx5UWvw/LM0xMeZdyISwh8MJOMqe8y8rTT+fqrr7wOq8YKBALcdM11PH31zeR9MJtfn3mTYSedyuLFVfr4eNVThUAg9CNKWc/fmDDIzMzkjYn38PyIROLjkgEY0F0Z+eKjHDdgiCfbFXtpyZIlzH7nXUZ37oHPXejo8Lzm3H/7nRz5RV/8fr/HEdY8X3z+OdnzlzC8dc89ZYfu3M5do8cw9cPpHkZWyXzxSEJFHj1dUHaVGsiSvzFh8MMPP3Bs+zzi4/Yu+xsbIwzslMu3387ijDNK3wugpvl65kyOSq2/J/EDJPv9tI9P4pdffqFnz56ltDah+Pz9jziy/r7DTGnJqfBXFps2baq5Qy6BXDT7D6+jiDiW/I0Jg4SEBHbm7T+Ktis/hrqJZc8VqGkSkpLYVswt1dxAoMx9EUxoEpOSyCnI2q88XwM1+xFT1Srb1a8msTF/Y8Lg8MMP5/vVyfy5NW9P2ZYdBXy2JIH+/ft7GJk3Bg0ZwteZGWTn7/19rN2eyUZ/DAcffLCHkdVcfzv3bD7fvIKCoC9dv2Wsp1GHNqSm1vBdF23Mv9ys529MGPj9fm4bP4nRYy+lS9pOYnywYH0CN457muTkZK/Dq3JpaWmMGncn9901jg7xSeQGAvzp93H/M0/bo5KVpHv37gy+7AImPDeZdvF12VKYgzaqzSOPPON1aJVMozqJh0pK2IUwoqWnp+vcuXPLrmhMmBUWFvLLL78QCAQ49NBDiY2N7u/Xubm5/PzzzyQkJHDIIYdY4q8CO3bsYOHChdSrV4/27dt7HQ4AIjJPVdMr49o9D66tP0w+KuT2cUfOqLTYqrPo/pfJmDCLiYmhe/fuXodRbcTHx9OrVy+vw4gqKSkp9OnTx+swqk5hLmQtLbue2Yclf2OMMRHMNvYJhSV/Y4wxEUsUJFDzhq8rmyV/Y4wxkc0m/JWbJX9jjDGRyxcPye0qcIFVYQslkljyN8YYE7kCObBjiddRRBxL/sYYYyKXYrf9Q2DJ3xhjTGSrgevVVDZPkr+IXAdcjPOd7RdgBNAEeAOoB8wHzlPVPBGJB14BegKbgbNVdaUXcRtjjKlmfPFQqyKLGa0PWyiRpMqTv4g0A64GOqvqLhF5CzgHGAxMUNU3RORZ4CLgGffPraraTkTOAcYDZ1d13MYYY6qhQA5s/93rKCKOV2ttxgKJIhILJAEbgOOAt93XXwb+5v58mnuO+/rxIkH7hBpjjIleu8f8bWOfcqny5K+q64CHgNU4ST8TmAdsU9UCt9paoJn7czNgjdu2wK1fvypjNsYYU40FNPQjSlV58heRuji9+dZAUyAZOKmYqrv/qxTXy9/vv5iIjBSRuSIyd9OmTeEK1xhjTLWmoIHQjyjlxYS/E4AVqroJQESmAUcCdUQk1u3dN2fvLIy1QAtgrTtMUBvYUvSiqjoRmAjOrn6V/imMMcZ4z5cAKR0rcIH/hi2USOJF8l8NHC4iScAu4HhgLvAVcCbOjP8LgPfc+u+759+7r3+pNXEfYmOMMeVXmAOZi72OIuJUefJX1R9E5G2cx/kKgB9xeuwfAW+IyN1u2QtukxeAV0VkGU6P/5yqjtkYU7zc3Fzee3c6334+i4aNGjJ0+LBqs4e8iSJRPHYfKk+e81fVO4A7ihQvB3oXUzcHOKsq4jLGHLjc3FwuHHoByet9dKjVkszfN3H1p5dx4wO3cuxxx3odnokSqpU7dC8ik4GTgY2q2sUt6wY8C9QCVgLDVHV7Ce1jcO5ur1PVk90yAe7GyW2FwDOq+njlfYr9efWonzEmwr0//T2S1glHpnWjQVJd2tZtyakNj+ahu8YTiOJHqIwHKne2/0vAoCJlk4AxqtoVeBf4ZyntrwEWFSkbjjOXrZOqHowz3F2lbHnfamL16tXMnTuXevXq0bdvX+Li4rwOyZhS/efL72hfq8U+ZYlxCSTkxfLXX3/RpEkTjyIzUSUmAWp3qsAFZpf6qqrOEpFWRYo7ArPcn2cCnwK3FW0rIs2BIcA9wPVBL10OnKvq3LNQ1Y0hBF4hlvw9pqpMuH88cz74hG4JqWQS4HHu4ZHJz9OqVSuvwzOmRI2aNmbbb2tIS9677Iaqkh3IITU11cPITFQpyIGtFZrw10BE5gadT3SfHivNr8CpOBPTz8LpxRfnUeBGIKVIeVvgbBE5HdgEXK2qS8sdeQXYbX+P/fDDDyz+6HNu6JjOgFYdObPVwZxf7yBuv+4Gr0MzplTnnD+U+blL2JmXDTiJf/7mRfTs14vk5GSPozNRJVCBAzJUNT3oKCvxA1wIXCki83ASe17RCiKye57AvGLaxwM5qpoOPA9MPvAPGx6W/D328dvTOK5Bc4JXLG6RWpdAxlY2bqzyO0HGHLDWrVtz6yN38WnuHN7b9A1v/jWTOv2ac9s9d3odmok2WoEjlLdTXayqJ6pqT2Aq8Ecx1foCp4rISpwx/eNE5DX3tbXAO+7P7wKHhhZJ6Oy2v8d8Ph9azN9ARbAtDEx1d1S/o+j75UdkZGSQnJxMUlKS1yGZaBOTAHUqb8y/OCKSpqobRcQH3Ioz838fqjoWGOvWPwYYrar/cF+ejrOfzWSgP7AkpNArwJK/xwafdQbPXjeGdnUb4nOT/cptm4lrVJ+GDRt6HJ0xZRMR+7tqvFOQA1sqb5EfEZkKHIMzN2AtzmPqtUTkSrfKNOBFt25TYJKqDi7jsvcDU9zt7bNwtrivUpb8PdarVy/+d8YpPPjvd+mamMo2LWR1rDJh8vNeh2aMMRFBtfLukqrq0BJeeqyYuutxtqcvWv418HXQ+TacpwA8Y8m/Grjyums5Y+g5zJ8/n7p169K7d29iYmK8DssYYyKDLStRbpb8q4nGjRszeHBZd4qMMcbsQ4GAzY8qL0v+xhhjIldsRSf82a5+xhhjTGQpyEG3/O51FBHHkr8xxpgIJnbbPwSW/I0xxkQuBSpxtn9NZcnfGGNM5IpNgLoVGfPfFLZQIoklf2OMMZGrIAfdbGP+5WXJ3xhjTGSzMf9ys+RvKiw7O5ucnBzq1q1r+xEYY6qWPecfEkv+JmRZWVnce+utLJ07l6TYWApr1eLGe+/l0EOrfIMqY0zUkkpd3remsuRvQnbbddfRddMmhnfrBsDGnTu568ormTh9OvXr1/c4OmNMVIhNQOpXZMLfhrCFEkks+ZuQ/Pnnn2xZsoSjunTZU5aWnMzxderwwbRpDL/kEg+jM8ZEjYIcdFOV74gb8Sz5m5Bs2bKFhvHx+5U3Skjgj/XrPYjIGBONVEFtzL/cfF4HYCJT27Zt+WPXLvIKC/cpn7t1K+lHHeVRVMaYqKQS+hGlLPmbkMTHx3P+1VfzwK+/snDTJtZu387UpUvJbNaM/scc43V4xpioIWgg9CNa2W1/E7LTz848q4YAACAASURBVDqLth07Mu3VV9m+bRtHX3YZQ045hZiYGK9DM8ZECYmNx9egYwWusDpssUQSS/6mQg499FAOffBBr8MwxkQpzc8lsGmp12FEHEv+xhhjIlo0374PlSV/Y8IoPz+f+fPnU1hYSHp6On6/3+uQjKnxbJGf8rPkb0yY/Dh/PveMvpZDk334BB7ZUcD1d9/PkX3t6QdjKk1sQgXH/JeHLZRIYsnfmDDIzc3lnhuu4c5uTWiQlABAZk4eN9/8T7p++BkpKSkeR2hMDZWfS2CjjfmXlz3qZzwTCATIzs5GVb0OpcJmz55NekrMnsQPUDvBz9F1/cyaNcvDyEwwVSU7O5tAIOB1KCaM7FG/8rOev6lyqspbU6bw5qRJJGmAnNhYhl99DSefdprXoYUsLy+PONn/S0y8CHl5eR5EZIqa+elMJtw3gYKdhRTGFDJ0xDlcNPIi24kywik25h8KS/6myn0wfTrfTXqecZ07ER8Tw678fB576EFS69Th6P79vQ4vJIcffjgT787j7/kFJMY5/1vlFRby1eYcnurXz+PozNy5cxl/04McVe9oEuomUBAo4INnPyI+IZ7zLjjP6/BMRahAwG5il5clf1Pl3pz0PKM7tifeXQwoMS6OEe3a8vIzz0Rs8k9JSeHim25lzP3jOKFhIjHAzE27+L8rr6VBgwZehxf1nn/ieXrUOoyEWGdYJtYXS3q93rw2aYol/wgncfH40tpX4AqLwxZLJLHkb6pcTlYWKUUegUtLTmLrmrUeRRQeAwcP4bBevfny85kUFhQyYcAAGjdu7HVYBli/dj0tE9rsUxYXE0feLhuSiXSan0vhn8u8DiPi2L0SU+Watm7Diq3b9in7eeNGOnXr7lFE4dOwYUPOHnou5553niX+aqRX396syly5T1lmzjYaNW/kTUAmrFQl5KMsIjJZRDaKyK9Fyq8Skd9F5DcReaCEtnVE5G0RWSwii0TkCLe8u4jMFpEFIjJXRHqH5RdRDpb8TZUbNXYsz65azbz1G8jMyWX22nVM3ZjByGuv9To0U0NdNupSVsUvZ+nmJWTnZ7M6czX/y57NmDtv8jo0U1GKM+Yf6lG2l4BBwQUicixwGnCoqh4CPFRC28eAT1S1E9ANWOSWPwDcpardgdvd8yplyd9UuY4dOzJhyuss79KViTt2srFXH558401atGjhdWimhkpLS2PK9Cl0Pqsjf9RbQv1j6vD8WxPp3j3y7zaZyt3VT1VnAVuKFF8O3K+quW6djftFJZIKHA284NbJU9XdtzwVSHV/rg2sD+WTV4TUhGesi0pPT9e5c+d6HYYxxhhAROapanplXLtH87r69ahjQ25fZ+y7q4CMoKKJqjoxuI6ItAI+VNUu7vkC4D2cOwI5wGhVnVOkTXdgIrAQp9c/D7hGVXeKyMHAp4DgdMKPVNVVIX+IENiEP2OMMRFL83Mp3FChCX8ZIXwxiQXqAocDvYC3RKSN7tubjgUOA65S1R9E5DFgDHAbzp2D61T1HRH5P5y7AydU5EOUl932N8aYKJGXl8e8efNYuHBhjVhZEwAFDfhCPkK0Fpimjv8BAaDoM71rgbWq+oN7/jbOlwGAC4Bp7s//Bqp8wp/1/I0xJgp8+cWX3DX2HpIL6lJIAbH1Ajwz+UlatmzpdWgV5sEKf9OB44CvRaQD4GffoQNU9U8RWSMiHVX1d+B4nCEAcMb4+wNfu9ep8s0JLPkbY0wNt3HjRu4YfTfdk4/FnxgPwNbtGVx+4Sg+nPleRC9xLHEJxDZuV4ErLCj9+iJTgWOABiKyFrgDmAxMdh//ywMuUFUVkabAJFUd7Da/CpgiIn6c7QNHuOWXAI+JSCzOnIGRFfgAIbHkb4wxNdyH739Iw8JW+GPi95TVTWzAuq0+Fi1aROfOnT2MrmI0P4eCDX9U3vVVh5bw0j+KqbseGBx0vgDYbz6Bqn4H9AxXjKGw5G+MMTXczqxsYor55z6GOHJycjyIKJyie3e+UNmEP2OMqeEGDDqBTbJ6n0l+eQW5ZMdupUuXLh5GFgZauSv81VTW8zfGmBquU6dODBk6gI/emEmDQAsCFLAxZiXjHr4df5F9NiKNs6Wv9WPLy5K/McZEgX+OuYFT/jaEmZ98Qa2UZIacMp60tDSvw6owiYsnpnHbClxhTtlVaiBL/sYYEyU6depEp06dvA4jrDQ/l4L1y70OI+JY8jfGGBPBonvsPlSW/IuRmZnJpx9/zI7M7fQ79hg6dOjgdUjGGGOKo54s8hPxLPkXMW/uXMZdcz1HJNalVkwMd7/yOl1PGsDoW26O6IUwjDGmpor2R/1EZDvOJkEKJAG73J8BklU1pmgbS/5BAoEA9944lmtbdaFOQiIAR2ornvrkCxYMGUyPHj08jtAYY8w+4uKJbVqRCX/fhy0Ur6jq7u2BEZH5qnpY8HlxbSz5B1m2bBmNNGZP4gcQEfrVbcQXH82w5G+MMdWM5ueSv26F12FUKyISq6oF7mlccXUs+QeJj48nL1C4X3luYSHxCQkeRGSMMaZUKhDlt/2LmAW8LSIfA/3Yu5nQPmxlhCAtW7Ykt04tVmzbsqcsr7CAz7duYMjfT/cwMmOMMSWxFf72MRr4GOgKzMXZPng/1vMPIiLc99QTjB55GQ22rCc1Jo6Fu7YzfPR1tGnTxuvwTCVSVRYuXEhWVhZdu3YlKSnJ65CMMQdA/PHENmtdgSt8G7ZYvCQiPYCjcSb6faeqz5VW35PkLyJ1gElAF5xALwR+B94EWgErgf9T1a3iTLF/DGenpGxguKoWO4EhHJo3b87Ujz7g119/ZceOHdzSrRvJycmV9XamGli3bh1jrriIZppFaqzw0NZ8hl87hiGn/c3r0IwxZdC8XPLXrvQ6DE+JyLXAcGC6W/SiiLysqo+U1Marnv9jwCeqeqa7z3EScDPwhareLyJjgDHATcBJQHv36AM84/5ZaUSErl27VuZbmGpCVbnlqksZdRB0aNgUgJz8Qm587F66dO/BQQcd5HGExpiyRPujfsBFQC9VzQEQkftx1i0uMflX+Zi/iKTi3Jp4AUBV81R1G3Aa8LJb7WVgd7frNOAVdcwG6ohIkyoO29RQK1eupE7uVjo03POkDAlxMZzeMpEZ703zMDJjzAGxXf3AuYMe/Cx/DHuf8y+WFz3/NsAmnNsS3YB5wDVAI1XdAKCqG0Rk944TzYA1Qe3XumUbgi8qIiOBkeBM3DPmQOTm5pIQs/8/AElxMeTuyvYgImNMeShiu/rB88BsEXnXPf+7W1YiL5J/LHAYcJWq/iAij+Hc4i9JcV/N9vtGo6oTgYkA6enppX7jMWa3Dh06sDwnhq3ZedRNcrY2VVVmrNnJhVcO8Tg6Y0xZxO8ntnkrr8PwlKo+ISKzcB7tAximqj+V1saL5L8WWKuqP7jnb+Mk/79EpInb628CbAyq3yKofXNgfZVFa2o0n8/H6HEPMGbMNQxsHEsdv48v/syjzdFD6Natm9fhGWPKoHl55K1Z5XUYnhKRg4BtwAfBZapa4i+mypO/qv4pImtEpKOq/g4cj7MIwUKc5xHvd/98z23yPjBKRN7AmeiXuXt4wJhw6N2nD89N+5jPPv6Yrdu3cfUxx9W4bU+Nqclq0Nh9qD5g79r+ycBBwBKgc0kNvJrtfxUwxZ3pvxwYgTP58C0RuQhYDZzl1p2B85jfMpxH/UZUfbimpqtbty5nn3uu12EYY0IQ7clfVQ8NPheRXjiP/pXIk+SvqguA9GJeOr6YugpcWelBGWOMiTiqzmH2UtU5IvJsaXVshT9jjDERy+ePxx/lE/4ARKQxztC44jzj/zcREbcDvR9L/sYYYyJWIC+P3DWrvQ7DUyJyDnAfzqY+CjwBjFHVqSW1seRvjDEmogWifMwfZ4Xcnqq6BUBE6gFfA5b8jTHG1EBqy/u6MoN+3lZWZUv+xhhjIpb44/G3aOV1GF77CPhERHb39IfhPClXIkv+xhhjIlYgL5fc1ZW3yI+ITAZOBjaqahe3bBzOvjMBnAXphqvqfovPich4YPdSoeNU9U23fBRwLdAWaKiqGRWJUVXHisgpQH+36DFVfb+0NlG/ILIxxphIFvqmPge4PsBLwKAiZQ+q6qGq2h34ELh9v6hEhuAsZd8dZxb+P92N7QD+A5wAhOVbi4gcCvxXVUcDdwNrRaTU/G7J3xhjTEQLqIR8lEVVZwFbipRtDzpNpvgd9DoD36hqgaruBH7C/RKhqj+q6soQP25xXgTyRCQJ5zG/+3F3zi2JJX9jjDERzYstfUXkHhFZgzO+vl/PHyfZnyQiSSLSADiWffepCSefqu4ABgIzVPVEoGdpDWzM3xhjTMQSv5/4lgdV5BINRGRu0PlEd5fYUqnqLcAtIjIWGAXcUeT1z9xldv+Ls43990BBRQItLRwR6Yaz/P1Tu8tKa2DJ3xhjTMQK5OaRs6pCi/xkqGpxy80fqNdxZtvfUfQFVb0HuAdARF4HllbgfUozBufW/3zgM3duweTSGljyN8YYE9GqepEfEWmvqrsT+anA4mLqxAB1VHWzOyHvUOCzyohHVT8rcu3twGOltbHkb4wxJoJVbOy+zKs7z84fgzM8sBanhz9YRDriPOq3CrjMrZsOXKaqFwNxwLciAk4y/oeqFrj1rgZuBBoDP4vIDLdNqDFOxtnSdx+qOkJE7lLV/e5KWPI3xhgT0Soz+avq0GKKi51Jr6pzgYvdn3NwZvwXV+9x4PFwxYjzuGFJvimu0JK/McaYiCV+PwktW3odhqdUdVrRMhG5yn3ty+LaWPI3xhgTsQJ5eexatcbrMDzlrhh4CZAaVNxURK4HHlXV/cb/LfkbY4yJXFq5t/0jxJXASThzC8B5zO9rnLUFsotrcEDJ353E0A9oCuwCfgU+3719oDHGGOONA1upr4ZbX3TFQBHJKC1Hl5r8RWQ4cDWwApgH/A4kAEcBN4nIr8BtqlqhhyyNMcaYUEh8HAkHVdbCedWbiNRW1UxVPb5IeRIwpbS2ZfX8k4G+qrqrhDfuDrQHLPkbY4ypcpqbz66Va70Owyv/AzruPhGR3sBFwHHAB6U1LDX5q+pTZby+4MBjNMYYY8JLieox/x9F5APgO+D/gDU4K/1doaqFpTU80DH/1sBVQKvgNqp6aogBG2OMMWERrWP+qnqOiByLM9O/ATATWFxW4ocDn+0/HWdRgw9wVjQyxhhjqgXdf3G7qKGqXwFfiUhtYCjwiogUAC+q6qSS2h1o8s9xVySKaEuWLGHlypW0a9eONm3aeB2OMcUKBAJ8+smnfPbeRyQlJ3HGP87hsMMO8zosY6oln99PYqvonPAXTFUzgWeBZ0XkEJwd/kp0oMn/MRG5A2fjgNygN5sfaqBVKScnhxuvvIKClStok+Dn37tyqdu5C/dMmEBcXJzX4Rmzh6oyetQ1ZP24hh512pBbkM29l97EqZcP4/wLh3sdnjHVTiA3n+wVUTvhr1iq+hswurQ6B5r8uwLn4cwg3H3bX93zau+5xx+nw8YNnHRIpz1l/166hCkvvcTwSy7xMDJj9jV//nwyflzB35sfvqesRWpDXnjuZU4/6wxSUlI8jM6Y6keJ3jH/ivAdYL3TgTaq2l9Vj3WPiEj8AN9+8jEntNz3ttDJrVry6TvveBSRMcWb/e1/ae9P26csxhfDQf76LFy40KOojKneVCXkI1odaM//J6AOsLESY6k0GlBifPt+z4nx+QiozV001Utak0asKNx/WY1thdk0aNDAg4iMqf6s519+B5r8GwGLRWQO+475R8Sjfr369+fbn+bTv2XzPWVfrFpD/8GDPYzKmP0NGnwSLz7+HB2zm9EwqQ4Ai7esJr55Xdq2betxdMZUPz6/n+RWzcuuWJJvwxeLF0SkMc7CPkNwFvyJBbYC/wGmqOqM4todaPK/IxxBeuWKG27g6hHDWbp4Ka3j4/h9Vy470hrz2MhLvQ7NmH2kpKQwYfLT3HHDWPLX7SQ/UMBBXTow4cEnvQ7NmGqpMC+PrJXrvA7DEyJyJs4z/m8B/1DV5W55MtATOE1Ebii6/C+Uvba/qOObsupU6BNUstTUVCb/+21++OEHVvzxB0M7daJnz56I2K0iU/107NiRNz6cxpYtW/D7/dSqVcvrkGqEVatWMfGxp1j40y80bdGcS669ku7du3sdlgmDKB7A/URV3y5aqKo7gVnALBEp9h+Qsib8fSUiV4lIy+BCEfGLyHEi8jJwQahRVyWfz8cRRxzBuf/4B+np6Zb4TbVXr149S/xhsnr1aq4cOoLGv27j0gZ9SN+UyF0jr+O7WRF+z9cAoU/2i/QJf6qaVdJrIvJMaXXKuu0/CLgQmOou8bsNZ1e/GJxn/ifY+v7GmOpu0pPPMCC5LR3qNQOgaUp9zonrxZPjH+aoo/t5HJ2pCJ8/rmJj/t+FLxaviMhy2GeZQwWaicgxwDPFLdJXVvKfBlypqk+LSBzO2sG7VHVbmGI2xphKt3DBL1xYt+c+ZXUSksn6M9OjiEy4FObms2NFdI75B0kvcq7A18CpwMfAfsm/rNv+LwGfisjNAKq6wRK/MSbSHNS2NWt3ZOxTlp2fi79WokcRmXCK1tv+u6nqliLHVmCyqq4C5hXXpqwtfd8SkY+A24G5IvIqQXMrVPWRMMZvjDGV4qKrLmfsiCs5259A/cRUdubn8M7aeQwfc4XXoZkwiOIJfwC4y+8XS1XPK678QB71ywd2AvFACvZ7rpFycnJYsGAB8fHxdOvWDZ/vQBd/NKb669y5M7c9+QCP3/cQW1dvIqFWEuffPIpTTouIpUpMGWpKD74CdgT9nIxzu7/UvXfKetRvEPAI8D5wmKpmVzRCU/3M/HQmE+68j1YxdcmjkAx/Lg9NfIJ27dp5HZoxYZPeqxevTHvT6zBMmPni/dSqyIS/78MXi1eK3oUXkfHAB6W1KavnfwtwlrtDkKmBNm7cyITb7uGCZkeTEOsHYFP2NkZfehXTZn5kdwCMMdVaYW6eTfgrQlXzRGSTiMSoamFxdUr9l11V+1nir9k+nfEx3f3N9iR+gIZJdaib4+fXX3/1MDJjwk9VycvLo5qvS2bKRdAKHDWBiLQTkUdE5A4RSRWRBOCWkhI/HPiufqaG2pW9C79v/xtAfl8Mubm5xbQwJjJ9+fnnnD3wJM47/kROP+Y4XnvxJfsSUEMEVEI+aoh3gDU44/1P4szNe7W0Bge6tr+poY4fOIB/vvQ23bQtPnG+C2bn57CmMJNu3bp5HJ0x4TF//nwm3XYX13ToQmp8AnmFBbz84qvExcVx9j+GeR2eqQCfP45arZuGfoEfSn9ZRCYDJwMbVbWLW3YWcCdwMNBbVeceaNsir48GHgQaqmpG0dfLYaeqTnCv+bN72z+ptAaW/KNc27ZtGXj+Gbzyyjt0jW9KnhbyU/46xoy/E7/fX/YFjIkArzz9DOe2bEtqfAIA/phY/tHuYMZPftGSf4QrzM1n+/L1lfkWL+H0pl8JKvsV+DvwXAhtARCRFsAAYHUYYvxKREa471MoImXO1rbkb7h01OWcOGQQX33+JUnJidw4aBD16tXzOixjwmbD2rU0a95hn7LE2DgCOTa0VRNU5ti9qs4SkVZFyhYBZe4RU1zbIBOAG4H3KhojcBXOLf/ngDxgKjCqtAaW/A0ArVu3pvUlF3kdhjGVonuf3iyYv5D0JnsfCdu4M4vajdM8jMqES6SN3YvIqcA6Vf0pHJvMqWpqedtY8jfG1HgjLr+cK845l4JAgC4NG7Eqcxv//ms1tzzxmNehmTCo4LzNBiISPGY/UVUnViyikrlj8bcAJ4b5ugPdayrwhap+XFp9S/7GmBqvcePGPPvWG7z2wgs8P3c+Ldq1ZvyD/7KFrGoAX3wcKW0qMOFvLhmqWnRjnMrUFmgN7O71Nwfmi0hvVf0zlAuKyI04q/pNxtnd7xYR6aqqD5TUxpK/MSYqpKWlcf3YsV6HYcKsMDefzOUbvA7jgKnqL8Ce8SYRWQmkV3C2/wXuNXa513wNmAOUmPztOX9jjDERLPQd/Q5kTwARmYqzCHBHEVkrIheJyOkishY4AvhIRD516zYVkRmlta2UXwHk7k78AKqaC5S4wA9Yz9+YGi07O5vFixdTr149WrVq5XU41YqqsmzZMnbu3Ennzp3t0dYIVpkT/lR1aAkvvVtM3fXA4ANoG9ymVcjB7TVDROq6W/kiInUAG/M3Jhq9PfUNXn3iadonpbI5L4fYZo144OknqVu3rteheW7Dhg3889LLSd6xkxRfLMvyshl12y2ccGJY52CZKqDuEc1U9dYi59tEZF5pbSz5G1MD/frrr7z7xDOMbZ9OXEwMAD9tWscdo2/k8Ree9zg6b6kqN11xJafHpdK+XRsAsvPzefDOcRx8yCE0a9bM4whNecTEx5FakQl/pW58GxlE5O844/67H/lToJeIjAJeUtWXi7ax5G9MDTT99TcYWL/FnsQP0K1hMz75fQ6ZmZnUrl3bw+i8tXLlSuK3bqd9+1Z7ypLi4jihdkM+fPddLh1V6tooppopzMln2x+RM+GvktwLXAZsd88VeB0YDRS75aElf2NqoJ07dpDsj9+vPMEXQ05OTlQn/127dpEUs/8/fclxfjK27/AgIlMRirOLTZTLVtWvgwtEZJeqlnjr32b7G1MDHXvyYP6zed/1zjOyd5KbnEBaWnSvatehQwdWFeSSlbd3aV9V5bttGzlm0EAPIzOhqszZ/hHiyAMs28N6/sbUQCcMGMCn09/nxYU/0T25Plvyc/hv9hb+9cwTZa5HXtPFxsbyz3H/4sExYzk2tSG1YuP4PjOD1sf2o0ePHl6HZ8opJj6O2m2ahH6BBeGLxUMzpPj/sY8VkYmqOrLoC54lfxGJAebirG98soi0Bt4A6uFMwTjP3ZYwHmenop7AZuBsVV3pUdjGRASfz8dDTz/JnDlzmP3NLJqlpfHiqafYTH9X36P70X76NGa89z47MjO56sQBdO3aNeq/GEWiwtx8tv0R0sJ4NcnoUl57pLhCL3v+1wCL2Ds7cTwwQVXfEJFngYuAZ9w/t6pqOxE5x613thcBGxNJRITevXvTu3dvr0OpltLS0hh+ycVeh2HCINrH/FV1vog0BvrgTIOYo6ob3NcWF9fGkzF/EWkODAEmuecCHAe87VZ5Gfib+/Np7jnu68eXcHvDGGNMlFGEgIZ+1ARux/h74O/uMVtESl1gyKue/6M4+xinuOf1gW2qWuCerwV2P2zbDFgDoKoFIpLp1t9nHWQRGQmMBGjZsmWlBm+MMab6iPZFfoCbgZ6qugVAROoBXwNTS2pQ5clfRE4GNqrqPBE5ZndxMVX1AF7bW+BswTgRID093f4uGGNMFIiNj6Nu2wpM+Ps5fLF4LDPo521lVfai598XOFVEBgMJOGP+jwJ1RCTW7f03B3Y/p7QWaAGsFZFYoDawperDNsYYU90U5OazZVnUT/ibAXzibiQEMMwtK1GVj/mr6lhVbe5uZnAO8KWqDgO+As50q10AvOf+/L57jvv6l6pqPXtjjDF7FvkJ9agJVHUM8BjQ2T0ec8tKVJ2e878JeENE7gZ+BF5wy18AXhWRZTg9/nM8is8YY0x1o9SkxXpCpqofAh8eaH1Pk7+7HOHX7s/Lgf2eSVLVHOCsKg3MGGNMxAhE+b1gEdnO3vlxcYAf2KmqKSW1qU49f2OMMaZcYuPjqNeucegX+C18sXhFVVODz905dba8rzHGmJqpIDefzTbhbx+qOkNE7gVuLamOJX9jjDERS6HGLNYTKhE5I+g0Bmc5/OzS2ljyN8YYE9Fqyqz9ChgS9HMBsBJnddwSWfI3xhgTsWLj46jftgJj/gvDF4tXVPXC8rax5G+MMSZiFeTmkxGlY/4icjywWVWL3ZhYRBKBC1X1qaKvWfI3xhgTsXYv8hOllgLjRKQL8B9gCZADpOGM+7cGHi+uoSV/Y4wxEU2L3QKm5lPV1cAF7kY+JwBdgESc5fHvU9W5JbW15G+MMSaiRfsiP+5ufm+5xwGx5G+MMSZixcTHUb8ii/wsDl8skcSSvzEmam3evJmdO3fSvHlzfL4q3+fMhEFBTj6blkbnhL+KsORvjIk6W7du5fbrb2DzsuXUio1jS4xw071306v3ftuLmGpPonbMvyIs+Rtjos6YUVfTN6uA7p0OAyAzJ4d7r7uB5959h7S0NI+jM+Vlm7yXn93nMsZElTVr1pC3dj3d05rsKaudkMCxqQ34cPp0DyMzodj9qF+ox4EQketE5DcR+VVEpopIQpHXLxORX0RkgYh8JyKd3fJhbtnuIyAi3Sv+qSvOev7GmKiSmZlJnVj/fuX1/An8tXGTBxGZioiNj6VBRSb8LSn9ZRFpBlwNdFbVXSLyFnAO8FJQtddV9Vm3/qnAI8AgVZ0CTHHLuwLvlbQgT1Wz5G+MiSodOnRgRV42OQUFJMTu/Sfwf5kZDD3+OA8jM6EoyC2oigl/sUCiiOQDSTjP0e+hqtuDTpNxbkgUNRSYWmkRlpMlf2NMVPH7/Yy8cTQP3/cAJ9VvSoo/nv9u/pP4QzrSp08fr8MzIajgCn8NRCR4MZyJqjpx94mqrhORh4DVwC7gM1X9rOhFRORK4HrADxT3LfJsythspypZ8jfGlGrnzp1Mf2ca/9/encdZVdd/HH99hrmzMCzDNuwwssmoFSghbmQaBphppbmDS5oLqb8styyXrDQttRQNwXJLK1HDHTSsSEVBEReQXUQQkHUQZr2f3x/3gDPsc+fee+be+37yOI+593vOufdzv5yZzznn+73f7zvTZ1DarzcnnnpKwjvFff755/xz4kTem/Em3fv04XunnJrUjncjjzuOvv3789Tf/sbi9RsYeeEovn7UUfq6X5pqZG//z9x90K5WmlkbYkl7H2A98A8zO8PdH64XQ2z8/LvN7DTgWmB0ndc4Q9JX8gAAHm1JREFUGNjs7u81JtBEUvIXkV1av349Pzj5DMoqmvOlViV8MvcNfvD3p7jtz/fSr1+/hL3HBaeewpBmcHTbtix5eTEXTpzIb+4bn7D32Jm+ffvy02uvTdrrS2rk5kca1+Y/f49bfANY7O6rAczsCeBQ4OFdbP8YcM92ZafQhG75g5K/iOzGn/80noOqW3Nw194AdG3Zlq7lxdx2/U2M++uDCXmPB8aNY1heLkf27A5AaXFrSluv544bb2Dsw48k5D0kc1VXVrMquW3+S4EhZtac2G3/o4F6Y+abWV9333oacSx1TinMLAc4CRiazCAbKiPvcW3YsIGlS5eGHYZI2nt96n84sENpvbKuLduyYnHifr+mvzKVQ7t1qVdWWlzMyiUfJew9JLO5x7/s+bV9OvA48BbwLrG8Oc7Mbgx69gOMCb4KOItYu//oOi8xFFjm7osS+ZkbKyOv/NcvX8tF3zufg4cfwrU3/gIzjf4kEo827dqydtMmOha13lZWXVuLRZol8D3a89nmzXRp2XJbWVVtLZ6buPeQzJbsKX3d/Trguu2Kf1Fn/aW72fcVYEhyIotfRl75N48U8q2Sr/Pu828zZcqUsMMRSVun//Acnln5PtW1tQC4O5OXf8CxJ303Ye9xynnn8dDCj6gK3iPqzj8WLGLE97+fsPeQzNWYq/5sHhkwI6/8AcyMr7Tqz5OPPsExxxwTdjgiaemIoUNZccm5/PHucbRvVsia6s0MHXkM5118YWLf44cX8LNxf6JzXoSVFZUcOmIE516QuPdIprVr1/LK1Km4O1878kjat28fdkhZJVIQoaRvIzr8Namb8amTsckfwDA8muwbQiKZ7funncoJJ36PVatW0bZtW5o3b5749zj9dE446aSkvkcyTHnhBe696Qa+VlxEsxzjwjt+x+ifXsW3jm8yX+fOeNUV1XyqWf0aLGOTv7sze+NcRp98ftihiKS9vLw8unXrlvbvkUjl5eX86de/5FcD+1MYif0pHVZTyzW33syhRxxB27ZtQ44we0Sz+PZ9vDKyzX9z9RaeXfUKfY7aj28O/2bY4YhIBnrttdcY3LJgW+IHyM9txuGti5g2bVqIkWUXb+SSrTLyyr915zbc8dhd9O7dO+xQRCRD5ebmUrOT7FGLk5ubkX9amyxd+TdcRh6hxcXFSvwiklSHHHIId22q5NiKCooLYjO8lldW8b+NWzj/iCNCji57RPIjdGpMh78lCQslrWRk8hcRSbbCwkKuuvV33HDF5QwsKsCAtz6v4Me/+g0t64xZIMlVXVnNinnq8NdQSv4iTcCqVav45+MTWblsGYMOP4yjhw0jEomEHZbsweCDD+bRyS8zffp03J3LBw+msLAw7LCyipP8QX4ykZK/SMhmz57NdReO4Rut29G/sIj/Tn+Lxx98iLsf+Av5+flhhyd7kJ+fz9ChTWrY9qyTzYP1xEvJXyRkv/3ZzxlTui8lRS0A2LddB/65aB5PTXyCk087NeToRJq2SH6ETv0a0eafpdPAKPmLhKi8vJzo+g2UdO5Vr/yQjl148oUXlfxF9qC6sprlavNvMCV/kRDl5+ezpbYWd683AdX6ii0U9ygJMTKRNJHlY/THKyMH+RFJF3l5eQw44jBe+eSL6Wura2t5asVSThw9KsTIRNKDo4l94qErf5GQXXH9dVx/xZX8etbbdMgrYGnVFs768SUMGDAg7NBE0oJ6+zeckr9IyAoLC7nlj39gzZo1rFmzhtLSUvLy8sIOSyQtRAoaOcjPssTFkk6U/EWaiHbt2tGuXbuwwxBJK9UV6vAXDyV/ERFJW9k+QU+8lPxFRCStaWKfhlPyFxGRtJbNvfbjpeQvIiJpK5Ifocu+jejwtyJxsaQTJX8REUlbVZXVLPtQHf4aSslfRETSl0NU9/0bTMlfRETSmlJ/wyn5B9asWcO8efPo2rUrPXr0CDscERHZC5GCXLo2Zla/LG0xyPrk7+78/uZbeeWpF+mR15bVNeV0KOvJbXffQWFhYdjhiYjIblRV1PCx2vwbLOsn9nn+2eeY/dR/OKf7kQzr9BVO63Y4rRZs5ne//m3YoYmIyB55o/5lq6xP/o8/9BhHluxfbzrVr5b0ZdqUqSFGJSIieyvq8S97w8yWmNm7ZjbLzGbsZP2RZrYhWD/LzH5RZ92lZvaemb1vZpcl7lM3Ttbf9q/YvIX8ZpF6ZTmWA1F2mGNdRESanhRdv3/d3T/bzfr/uvu36haY2QHAecBgoAp4wcyedff5SYxzr2R98v/Gt0cw4y+TGdrlS9vKFq1bTt8vlynxi4g0cZGCCN327Rj/C6xKXCw7UQa87u6bAczs38B3gNDblbM++Z8+6gwuePkVnl46g30i7VhVW86SyEbuvfH+sEMTEZE9qKqoZuncRnX4a7/drfxx7j5uu20cmGxmDvxpJ+sBDjGzd4DlwE/c/X3gPeBXZtYO2AKMBHZoNghD1if//Px8Jvz1AV599VVmz5zFgF49+cawYRQUFIQdmjQRs2fP5tH7xrJqxXIOPOQITj37XNq2bRt2WCISaORt/8/cfdAetjnM3ZebWQkwxczmuvt/6qx/C+jp7pvMbCTwFNDX3eeY2S3AFGAT8A5Q07hwEyPrO/wB5OTkcPjhh3PRpWP41nHHKfHLNi9Pmcztl53LcdEl/KJfM9rPeoaLTj+RDRs2hB2aiATc41/27vV9efBzFfAksTb8uus3uvum4PFzQMTM2gfPJ7j7ge4+FFgLhN7eD7ryF9kld2fcbb/hloO70aog1il0WJ+OVEU/5e9/fYjzLhwTcoQiEimI0K1/Iwb5mbb71WZWBOS4e3nw+Bjgxu226QSsdHc3s8HELqzXBOtK3H2VmfUAvgscEn+wiZPyK38z625mU81sTvDVh0uD8rZmNsXM5gc/2wTlZmZ/MLMFZjbbzA5Mdczpat26dTz2178ybuy9zJkzJ+xw0k55eTktvWpb4t9qcJdi3p3+akhRiUg9wdj+8S57oSMwLWjPfwN41t1fMLMLzOyCYJsTgfeCbf4AnOK+7cUnmtkHwNPAxe6+LrEVEJ8wrvxrgMvd/S0zawnMNLMpwFnAy+5+s5ldBVwFXAmMAPoGy8HAPcFP2Y0333iD6y65ggGRzhRahBsefJIBI4/gyp//TN9i2EvNmzdnQ7VTG3Wa5XxRZ0vWfU6Xnl8OMTIR2SrW4S958/K6+yLgKzspv7fO47uAu3ax/xFJC64RUn7l7+4r3P2t4HE5MAfoChwPPBBs9gBwQvD4eOBBj3kdKDazzikOO61Eo1F+ecW1jOo0hMM6l3Fgpz6c0f1Q3n1uGrNmzQo7vLSRm5vL0SecxPh3PqGqJgrAyvIKHlhYzvdHnxNydCKylTdiyVahtvmbWSkwEJgOdHT3FRA7QQh6VULsxODjOrstC8qSd6qX5ubPn0+72gJa5jXfVmZmDGzejSnPPM/AgQNDjC69nP+jS7k/EmHM44+S57XktWnP5beOpVevXmGHJiLEEng0q9N4fEJL/mbWApgIXObuG3dzK3pnK3b4nzaz84HzgayflS8/P5+q6I7fJqmK1tCyqPlO9pBdycnJ4QcXjeHcCy8mGo3SrFmzsEMSkTryCnLp3pgOf/9LXCzpJJTkb2YRYon/EXd/IiheaWadg6v+znwx7tIyoHud3bsRG0ShnmDQhXEAgwYNyurTwJ49e+LtmvPxxtV0b9UBgKraal7f/BF/+O714QaXpsxMiV+kCaqqqOajJLb5Z6owevsbMAGY4+6/r7NqEjA6eDwa+Ged8lFBr/8hwIatzQOyc2bGrffcyb9ylvL48jd5dsUsxi+fxnnXXEZpaWnY4YmIJIyT/Il9MlEYg/wcBpwJHFVnBqSRwM3AMDObDwwLngM8BywCFgD3AReFEHPaKSkp4fCjjmTh5tXMWvcRPfr1YsBB+pakiGQeTenbcCm/7e/u09h5Oz7A0TvZ3oGLkxpUBrrxml+w+bVF/LjXcHJzmrF4xadcfMa5PPLMRFq2bBl2eCIiCaMOfw2nEf4y0Nq1a5n93zc4r8eR28r2Ke5E2Yq1PP3UPzntzDPCC05EJIHyCiL0bEyHvywdr0vJPwOtWLGCjnmtdijvXFDMknkLQ4hIJPXcnZdfeoln//YP3KMM/953OWb4cHJyNKVJJqmsqGaxOvw1mH4LMlDPnj35pGo9UY/WK1+85TO+NEjf8ZfscOsvb2LSTbdwzKYaRmx2Xrr1Dm68+pqww5IkUJt/wyn5Z6AWLVpw3Okn8vjS11lfsYnqaA3TP/2QT4tr+OaI4WGHJ5J0y5Yt453JL/GDfgfQrVVrurRsxag++7Hk1eksWLAg7PAkoZyoxb9kK932z1A/HHMRpX168diEhyjfuJEjTxjGFeedQ15eXtihiSTd7Nmz+VJhy3rzWJgZXy4oYtbbb9OnT58Qo5NEyiuIUNqINv//apAfySRmxvARIxg+YkTYoYikXElJCatqq3coXxWtYVDHjiFEJMlSWVHNork7jPsme6Db/iKScQ488EBWFhXw3upPt5V9uGY1iyPGoYceGmJkkgxq8284JX9JmWg0yqJFi1i+XGfpklw5OTnccf94Xm/Xghven8mN789kalEud/z5fnJzdcMzszhONO4lW+m3QFLirZkzueWan9A1t5bNNbXUFHfipjvvplOnRnw/V2Q3OnTowJ3j76OiogJ3p7CwMOyQJBmMrO64Fy8lf0m69evX85vLx3DToC60a54PwAcr13PlhefxlycmsZsZHUUaraCgIOwQJIliHf46x73/9GkJDCaNKPlL0r343HMcU5K3LfED7NexmA7LPmbu3LmUlZWFGJ2IpLNYh79Pwg4j7Sj5S9JtWLeG9vk7TofbJmKUl5eHEJGIZA4nmsVt9/FShz9JukO/9nWmrtxCbI6mmMqaWt7ZUM0BBxwQYmQiku4ciDbiX7bSlb8k3f7770/nIUdxyxtTGd6tJZura5m4tJzTx/yE5s2bhx2eiKQ5t+xN4vFS8pekMzOuvuEm3nzzTaY+O4mCoiJ+/vOTNMpahotGo0x+cTIvTXqW5kVFfOf0kxk4UHNLSGLlF0TYp3/XuPeftTqBwaQRJX9JCTNj8ODBDB48OOxQJAXcnSt+dBmbZy1mcJueVNSs4+YLfsK3LzqL00efGXZ4kkEqK6pYOPfjsMNIO0r+IpJwM2fOZO2sBZza46vbykpbl3DXPRM4/nvfoUWLFiFGJ5lka5u/NIw6/IlIwk2f9ir75XeoV9YsJ4d98trwwQcfhBSVZCbHqY17yVa68heRhOvQqYR3arfsUL4hWkG7du1CiEgyVX5BHr36d4t7//fV5i8ikhjDR47kwT+OY//NG+nQvBUAH3y2jNwubendu3fI0UkmqayoYoHa/BtMt/1FJOFatWrFbRPG8nTNEu776FXGLpnGgh65/O5Pd4UdmmSY2LQ+tXEve2Jm3c1sqpnNMbP3zezSnWzT2syeNrN3gm3ODsoHmNlrQdlsMzs5CVUQF135i0hS9O/fn0effZI1a9YQiURo1apV2CFJhkry7Hw1wOXu/paZtQRmmtkUd6/beeVi4AN3P87MOgAfmtkjwGZglLvPN7Muwb4vuvv6ZAa8N5T8RSSp1MYvyeVJTf7uvgJYETwuN7M5QFegbvJ3oKXFZilrAawFatx9Xp3XWW5mq4AOgJK/iIhIvPIL8ujTv3vc+89vQIc/MysFBgLTt1t1FzAJWA60BE529+h2+w4G8oCFcQebQEr+0uTMmzePe2+9laULFtCmQwdGjxnD4UOHhhbPC8+/wMP3jmfjug3s++X9+dEVP6ZHjx6hxSMiX6isqGLe3I8a8xLtzWxGnefj3H3c9huZWQtgInCZu2/cbvU3gVnAUUBvYIqZ/XfrdmbWGXgIGL39SUFY1OFPmpRFixZx9TlnM+zzcm7ev4yzigoZf83VvDR5cijx/P3Rx3jw+ts51nrxw06H0X3eFi4+7WxWrlwZSjwiUp83/nv+n7n7oDrLzhJ/hFjif8Tdn9hJGGcDT3jMAmAx0D/YtxXwLHCtu7+erHpoKCV/aVL+MnYsZ3bpQt+2bQHoWFTEJWX9mXD77SmPJRqN8uDY+zi5+2Ba5xdhZvRp04VDIt155P4HUh6PiOycE4172ZOgHX8CMMfdf7+LzZYCRwfbdwT2BRaZWR7wJPCgu/8jIR82QXTbX5qUBR+8z6k9699Sb5WfT1V5Oe5O7PcwNTZt2kRhNJdIs/q/JqWtOvDvdzVKnUjT4EQ9qSP1HQacCbxrZrOCsmuAHgDufi/wS+AvZvYuYMCV7v6ZmZ0BDAXamdlZwb5nufssQqbkL01K77L9mLdsKft3+GJo2PLKSiItWqQ08QO0aNGCLVZDdW1NvROAj8o/o9/BZSmNRUR2Lr8gn779e8a9/0fTdr/e3acRS+i722Y5cMxOyh8GHo47uCRS8pcm5ayLLuLK0aO4IBKhtLiYNZs3M27+As666uqUx5KTk8Np55/NP+5+iOM6D6BFXiFLNqzkf5Ufcd/Zv055PCKyo8qKSubNXRx2GGlHyV+alN69e/PL+8Zzz62/Zfl7H9CybVtG33AjRx51VCjxnDbqDFq1Kebheyfw+Zpy+pTtyx+vnkDnzp1DiUdE6tva4U8axtw97BgSbtCgQT5jxow9bygiIklnZjPdfVAyXjsvt4WXFH857v0/WfNa0mJrynTlLyIiaaugIJ9+/Uvj3v+T/72WuGDSiJK/iIikrYqKSj6cuyjsMNKOkr+IiKQxp4kMmpdWlPxFRCRtORBN7qx+GUnJX0RE0lhyZ/XLVEr+IiKStgoK8tm3ER3+Vr36duKCSSNK/iIikrYqKiqZO6dJzJKbVpT8RUQkjem2fzyU/EVEUiAajfLvV17h1anTaN+5hO+c+F1KSkrCDivtOai3fxyU/EVEkqympoYx515A5Ydr6FfYlUVVCxn1wN/41d23ctCgrBtcLsEcT+6sfhlJyV9EJMleeP55aj9cz/Aug7eV9a7sxk1X38ATkyelfMbKTFJQkE///r3j3v/V6R8mMJr0oeQvIpJkU597iQNaldYra53fgmbrali9erVu/zdCRUUlc+bODzuMtKPkLyKSZK3bFLOpajXQrl55RW01hYWF4QSVKVxt/vHICTsAEZFMd9KZp/D6pg+pqq3eVvbeZ4voe+B+tGzZMsTIMkGst3+8S7bSlb+ISJKVlZVx/rWXMPaWO+mQ05ry2i106t+D3/72V2GHlvZibf594t5/+huLExhN+lDyFxFJgeNOOJ5vjhzBwoULadOmDZ06dQo7pIywpaKSOXPnhR1G2lHyFxFJkby8PMrKysIOI8M4qM2/wZT8RUQkzSn5N5SSv4iIpDX19m84JX8REUlbBQX59C/rG/f+b765IoHRpA8lfxERSVsVFRXMmZOdo/Q1hpK/iEgINm7cyLNPP8Onn6zg4MMPYciQIeTkaOiV+Oi2f0Mp+YuIpNi8efO49KyL6OOdad2sBWMnTuPBsg7cNX4subn6s9ww6u0fD51mioik2PU/vZZhRYP4askB9GtXyjGdhlA1Zx3PTHom7NDSkkb4a7i0OcU0s+HAnUAzYLy73xxySCIiDbZx40bKV66nfcc29cq/1KoPkyc9zwnfPSGkyNJTQUEBZWX7xr3/jBmv73GbPeUfM8sHHgQOAtYAJ7v7EjMrBeYAWzslvO7uF8QdbAKlRfI3s2bA3cAwYBnwpplNcvcPwo1MRKRhIpEI1dEd55/fUlNJUYuiECJKb7EOf3OS9vp7mX/OBda5ex8zOwW4BTg5WLfQ3QckLcA4pctt/8HAAndf5O5VwGPA8SHHJCLSYIWFhez/1S8xZ+3CbWW10VreKH+fk886LcTI0pWDN2LZs73JP8cDDwSPHweONjNL2EdMgrS48ge6Ah/Xeb4MODikWEREGuW639zI5RddxlML/k3rZkUsr17DaReOYtCgQWGHlnYKCgroX9Y/7v1nznxjT5vsTf7Zto2715jZBr6Yv3kfM3sb2Ahc6+7/jTvYBEqX5L+zM6h6p2xmdj5wfvC00szeS3pU6as98FnYQTRRqpvdU/3sWqPqZvKbUznrB2cnMJwmJf5G+T3YvHnzizNnvtG+ES9RYGYz6jwf5+7j6jzfY/7ZzTYrgB7uvsbMDgKeMrP93X1jI+JNiHRJ/suA7nWedwOW190g+M8aB2BmM9xdp9C7oPrZNdXN7ql+dk11s2vbJdeEcvfhyXrtwB7zT51tlplZLtAaWOvuDlQGcc40s4VAPyBp9bG30qXN/02gr5ntY2Z5wCnApJBjEhGRzLc3+WcSMDp4fCLwL3d3M+sQdBjEzHoBfYFFKYp7t9Liyj9oQxkDvEjsqxb3u/v7IYclIiIZblf5x8xuBGa4+yRgAvCQmS0A1hI7QQAYCtxoZjVALXCBu69N/afYkfne9XZMK2Z2/nZtNlKH6mfXVDe7p/rZNdXNrqlump6MTP4iIiKya+nS5i8iIiIJknHJ38yGm9mHZrbAzK4KO55UM7PuZjbVzOaY2ftmdmlQ3tbMppjZ/OBnm6DczOwPQX3NNrMDw/0EyWdmzczsbTN7Jni+j5lND+rmb0GnHswsP3i+IFhfGmbcqWBmxWb2uJnNDY6hQ3TsxJjZ/wW/U++Z2aNmVpDNx46Z3W9mq+p+rTqeY8XMRgfbzzez0Tt7L0m8jEr+9sUwjCOA/YBTzWy/cKNKuRrgcncvA4YAFwd1cBXwsrv3BV4OnkOsrvoGy/nAPakPOeUuJTbe9la3ALcHdbOO2FCdUGfITuD2YLtMdyfwgrv3B75CrJ6y/tgxs67AJcAgdz+AWMevrcO4Zuux8xdg+6/ZNehYMbO2wHXEBs0ZDFy39YRBkiujkj8aBhh3X+HubwWPy4n98e5K/eEnHwC2zh5yPPCgx7wOFJtZ5xSHnTJm1g04FhgfPDfgKGJDcsKOdZNWQ3Y2hpm1ItY7eQKAu1e5+3p07GyVCxQG3+NuTmwAl6w9dtz9P8R6ttfV0GPlm8AUd1/r7uuAKex4QiFJkGnJf2fDMHYNKZbQBbcaBwLTgY7uvgJiJwhASbBZttXZHcAVsG0uz3bAenevCZ7X/fz1huwE6g7ZmYl6AauBPwfNIuPNrAgdO7j7J8BtwFJiSX8DMBMdO9tr6LGSNcdQU5NpyX9vhmHMCmbWApgIXLaHoSSzps7M7FvAKnefWbd4J5v6XqzLRLnAgcA97j4Q+JwvbtvuTNbUT3Ar+nhgH6ALUETsVvb2svXY2ZNd1YfqKSSZlvz3ZhjGjGdmEWKJ/xF3fyIoXrn1lmzwc1VQnk11dhjwbTNbQqxJ6ChidwKKg1u5UP/zb6sbqzNkZyoDTrFlwDJ3nx48f5zYyYCOHfgGsNjdV7t7NfAEcCg6drbX0GMlm46hJiXTkn/WDwMctCtOAOa4++/rrKo7/ORo4J91ykcFvXGHABu23rbLNO5+tbt3c/dSYsfGv9z9dGAqsSE5Yce62WHIzhSGnFLu/inwsZltnYTlaOADdOxA7Hb/EDNrHvyOba0bHTv1NfRYeRE4xszaBHdXjgnKJNncPaMWYCQwD1gI/CzseEL4/IcTu202G5gVLCOJtTe+DMwPfrYNtjdi35BYCLxLrDdz6J8jBfV0JPBM8LgX8AawAPgHkB+UFwTPFwTre4UddwrqZQCxSUdmA08BbXTsbKubG4C5wHvAQ0B+Nh87wKPE+j9UE7uCPzeeYwU4J6inBcDZYX+ubFk0wp+IiEiWybTb/iIiIrIHSv4iIiJZRslfREQkyyj5i4iIZBklfxERkSyj5C8iIpJllPxFUshiUy4vDmYzIxjcZLGZ9TSzzhZMM9yA17vNzI5KTrQikqmU/EVSyN0/Jjad6c1B0c3AOHf/CPgxcF8DX/KP7H78fRGRHWiQH5EUC+ZemAncD5wHDHT3KjNbBJS5e6WZnUVsOtRmwAHA74A84EygEhjp7muD15sJHOux4XlFRPZIV/4iKeaxiWF+CtxObNbFKjPbB1jn7pV1Nj0AOA0YDPwK2Oyx2fZeA0bV2e4tYpMWiYjsFSV/kXCMIDYu+gHB887A6u22meru5e6+mth88E8H5e8CpXW2W0VsmlkRkb2i5C+SYmY2ABgGDAH+L5j6dAuxyWDqqnsXIFrneRTIrbOuINhfRGSvKPmLpFAwHew9xG73LwVuBW4jNhNlaZwv24/YTHMiIntFyV8ktc4Dlrr7lOD5WKA/MAhYaGZ9GvJiQefBPsSm4RUR2Svq7S/SRJjZd4CD3P3aBu5zoLv/PHmRiUimyd3zJiKSCu7+pJm1a+BuucS+Bigistd05S8iIpJl1OYvIiKSZZT8RUREsoySv4iISJZR8hcREckySv4iIiJZ5v8B9Y2O+TKe0JkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x23311501588>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "GSLIB.locmap(df,'X','Y','Porosity',xmin,xmax,ymin,ymax,vmin,vmax,'Well Data - Porosity','X(m)','Y(m)','Porosity (fraction)',cmap,'locmap_Porosity')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look carefully, and you'll notice the the spatial samples are more dense in the high porosity regions and lower in the low porosity regions.  There is preferential sampling.  We cannot use the naive statistics to represent this region.  We have to correct for the clustering of the samples in the high porosity regions. \n",
    "\n",
    "Let's try cell declustering. We can interpret that we will want to minimize the declustering mean and that a cell size of between 100 - 200m is likely a good cell size, this is 'an ocular' estimate of the largest average spacing in the sparsely sampled regions.   \n",
    "\n",
    "Let's check out the declus program reimplimented from GSLIB."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {},
   "outputs": [],
   "source": [
    "wmin = 0.0; wmax = 3.0\n",
    "noff = 3\n",
    "csize = 50.0\n",
    "\n",
    "# interactive calculation of the sample set (control of source parametric distribution and number of samples)\n",
    "style = {'description_width': 'initial'}\n",
    "l = widgets.Text(value='                                              Cell-based Declustering, Michael Pyrcz, Associate Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n",
    "csize = widgets.FloatSlider(min = 5.0, max = 10000, value = 100.0, step = 20.0, description = 'Cell Size',style=style,orientation='horizontal',\n",
    "                          layout=Layout(width='250px', height='50px'),continuous_update = False)\n",
    "csize.style.handle_color = 'gray'\n",
    "\n",
    "noff = widgets.IntSlider(min=1, max = 30, value = 5, step = 5, description = 'Number of Offsets',style=style,\n",
    "                        orientation='horizontal',layout=Layout(width='250px', height='50px'),continuous_update=False)\n",
    "noff.style.handle_color = 'gray'\n",
    "\n",
    "uipars = widgets.HBox([csize,noff],) \n",
    "ui = widgets.VBox([l,uipars],)\n",
    "\n",
    "def f_make_declus(csize,noff,): # function to take parameters, make sample and plot\n",
    "    \n",
    "    wts =  declus(df,'X','Y',noff=noff,csize=csize)*len(df)\n",
    "    avg_wts = np.average(wts,axis=1)\n",
    "    \n",
    "    plt.subplot(131)\n",
    "    im = plt.scatter(df['X'],df['Y'],c=wts[:,0],cmap=cmap,vmin=wmin,vmax=wmax,alpha=0.8,linewidths=0.8,\n",
    "            verts=None,edgecolors=\"black\",)\n",
    "    plt.title('Spatial Sample Data with Weights')\n",
    "    plt.xlim(xmin, xmax); plt.ylim(ymin, ymax); plt.xlabel('X (m)'); plt.ylabel('Y (m)')\n",
    "    cbar = plt.colorbar(im, orientation=\"vertical\", ticks=np.linspace(wmin, wmax, 10))\n",
    "    cbar.set_label('Weights',rotation=270, labelpad=20)\n",
    "     \n",
    "    if csize < 1000:\n",
    "        major_ticks = np.arange(0, 1001, csize)\n",
    "        plt.gca().set_xticks(major_ticks); plt.gca().set_yticks(major_ticks)\n",
    "        plt.grid()   \n",
    "\n",
    "    plt.subplot(132)\n",
    "    plt.hist(df['Porosity'],color='red',alpha=0.4,bins = np.linspace(vmin,vmax,20),edgecolor = 'black',label = 'Naive')\n",
    "    plt.hist(df['Porosity'],weights = avg_wts,color='gray',alpha=0.4,bins = np.linspace(vmin,vmax,20),edgecolor = 'black',label= 'Weighted')\n",
    "    plt.ylabel('Frequency'); plt.xlabel('Porosity (%)'); plt.title('Naive and Declustered Distributions')\n",
    "    plt.legend()\n",
    "    \n",
    "    plt.xlim(0,30); plt.ylim(0.0,int(len(df)*0.5))\n",
    "    \n",
    "    avg = np.average(df['Porosity'])\n",
    "    declus_avg = np.average(df['Porosity'],weights = avg_wts)\n",
    "    P10 = weighted_percentile(df['Porosity'],10,weights=np.ones(len(df)))\n",
    "    P90 = weighted_percentile(df['Porosity'],90,weights=np.ones(len(df)))\n",
    "    \n",
    "    plt.gca().add_patch(Rectangle((21.7,len(df)*0.345),7.8,len(df)*0.07,facecolor='white',edgecolor='black',linewidth=0.5)); \n",
    "    plt.text(22,len(df)*0.4,'Naive Mean   ' + str(round(avg,1))); \n",
    "    plt.text(22,len(df)*0.38,'Declus. Mean ' + str(round(declus_avg,1))); \n",
    "    plt.text(22,len(df)*0.36,'Correction    ' + str(round((declus_avg-avg)/avg*100,0)) + '%'); \n",
    "    \n",
    "    plt.gca().add_patch(Rectangle((21.7,len(df)*0.235),4.5,len(df)*0.08,facecolor='white',edgecolor='black',linewidth=0.5)); \n",
    "    plt.text(22,len(df)*0.3,'Mean ' + str(round(np.average(df['Porosity'],weights = avg_wts),1))); \n",
    "    plt.text(22,len(df)*0.28,'P90    ' + str(round(P90,1))); \n",
    "    plt.text(22,len(df)*0.26,'P10    ' + str(round(P10,1))); \n",
    "    plt.text(22,len(df)*0.24,'n        ' + str(len(df)))\n",
    "\n",
    "    plt.plot([avg,avg],[0.0,45],color = 'red')    \n",
    "    plt.plot([P10,P10],[0.0,45],color = 'black',linestyle='dashed')\n",
    "    plt.plot([declus_avg,declus_avg],[0.0,45],color = 'black')\n",
    "    plt.plot([P90,P90],[0.0,45],color = 'black',linestyle='dashed')\n",
    "    \n",
    "    plt.subplot(133)\n",
    "    plt.boxplot(wts.T)\n",
    "    plt.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.2, wspace=0.2, hspace=0.2)\n",
    "    for label in (plt.gca().get_xticklabels() + plt.gca().get_yticklabels()):\n",
    "        label.set_fontname('Arial')\n",
    "        label.set_fontsize(8)\n",
    "    plt.scatter(np.linspace(1,len(df),len(df)),avg_wts,color='blue',alpha=0.8,marker='x',s=50,linewidth=1,label='Average')\n",
    "    plt.xlabel('Data Index'); plt.ylabel('Data Weight Over Cell Mesh Offsets'); plt.title('Weights Variation Over Random Cell Mesh Offsets')\n",
    "    plt.xlim(0,len(df)+1); #plt.ylim(0,1.5)\n",
    "    plt.legend(loc = 'lower right')\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",
    "# connect the function to make the samples and plot to the widgets    \n",
    "interactive_plot = widgets.interactive_output(f_make_declus, {'csize':csize, 'noff':noff})\n",
    "interactive_plot.clear_output(wait = True)               # reduce flickering by delaying plot updating\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interactive Cell-based Declustering Demostration\n",
    "\n",
    "* select the cell size and number of cell mesh offsets and visualize the declustering method \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) | [GeostatsPy](https://github.com/GeostatsGuy/GeostatsPy)\n",
    "\n",
    "### The Inputs\n",
    "\n",
    "* **Cell Size**: the size the of the cells in the mesh, **Number of Offsets**: number of cell mesh offsets to average to calculate the data weights "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9071467a33bd42d385d4544edf62e087",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='                                              Cell-based Declustering, Michael Pyrc…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "788677b47308424a8848f4ecced3072c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output(outputs=({'output_type': 'display_data', 'data': {'text/plain': '<Figure size 432x288 with 4 Axes>', 'i…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(ui, interactive_plot)                            # display the interactive plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([   0,   10,   20,   30,   40,   50,   60,   70,   80,   90,  100,\n",
       "        110,  120,  130,  140,  150,  160,  170,  180,  190,  200,  210,\n",
       "        220,  230,  240,  250,  260,  270,  280,  290,  300,  310,  320,\n",
       "        330,  340,  350,  360,  370,  380,  390,  400,  410,  420,  430,\n",
       "        440,  450,  460,  470,  480,  490,  500,  510,  520,  530,  540,\n",
       "        550,  560,  570,  580,  590,  600,  610,  620,  630,  640,  650,\n",
       "        660,  670,  680,  690,  700,  710,  720,  730,  740,  750,  760,\n",
       "        770,  780,  790,  800,  810,  820,  830,  840,  850,  860,  870,\n",
       "        880,  890,  900,  910,  920,  930,  940,  950,  960,  970,  980,\n",
       "        990, 1000])"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.arange(0, 1001, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comments\n",
    "\n",
    "This is an interactive demonstration of cell-based declustering to correct for sampling bias. Much more could be done, I have other demonstrations on the basics of working with DataFrames, ndarrays and many other workflows availble 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
}