SubsurfaceDataAnalytics_Spatial_Bootstrap.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<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",
    "### Spatial Bootstrap with the Number of Effective Data\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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: Spatial Bootstrap for Subsurface Data Analytics in Python \n",
    "\n",
    "Here's a simple workflow, demonstration of spatial bootstrap for subsurface modeling workflows. This should help you get started with building subsurface models that integrate uncertainty in the sample statistics.  \n",
    "\n",
    "#### 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.  \n",
    "\n",
    "#### Spatial Bootstrap\n",
    "\n",
    "Journel (1993) developed methods for spatial bootstrap based on bootstrap resampling accounting for the locations of the data and the spatial continuity model.  \n",
    "\n",
    "* One method to perform spatial bootstrap for uncertainty of a statistic is to adjust the number of data, $n$, to the **number of effective data**, then use the number of effective data instead of number of data of resamples with replacement for each bootstrap realization.\n",
    "\n",
    "This number of effectve data may be calculated by:\n",
    "\n",
    "* building multiple unconditional simulated realizations at the data locations only \n",
    "\n",
    "* calculating the variance of the average of each simulated realization\n",
    "\n",
    "* then calculate the number effective data by manipulating the standard error in the average equation:\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma^2_\\overline{x} = \\frac{\\sigma^2}{n}\n",
    "\\end{equation}\n",
    "\n",
    "to be expressed as:\n",
    "\n",
    "\\begin{equation}\n",
    "n^{'} = \\frac{\\sigma^2}{\\sigma^2_\\overline{x}}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "where $\\sigma^2_\\overline{x}$ is the variance of the average over bootstrap realizations $\\ell = 1,\\ldots,L$, $\\sigma^2$ is the variance / sill of the problem.\n",
    "\n",
    "#### Simulating Only at the Data Locations\n",
    "\n",
    "We ultilize LU Simulation to efficiently calculate Gaussian realizations only at the data locations, based on the Lower-Upper Decomposition of the left-hand / redudancy martix of the simple kriging.  \n",
    "\n",
    "* We multiple the Lower matrix ($n \\times n$) with a random Gaussian vector ($1 \\times n$) to calculate each realization at the data locations.\n",
    "\n",
    "* We then take the average of each realization and then calculate the variance of the average over enough realizations.\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": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geostatspy.geostats as geostats\n",
    "import geostatspy.GSLIB as GSLIB"
   ]
  },
  {
   "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": 21,
   "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",
    "import scipy                              # statistics\n",
    "import scipy.linalg                       # linear algebra library\n",
    "import math                               # trig, etc.\n",
    "import random                             # bootstrap"
   ]
  },
  {
   "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": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    " cmap = plt.cm.inferno"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Declare Functions\n",
    "\n",
    "I have included the n_effective function and dependencies here.  I will add this shortly to the GeostatsPy Package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numba import jit # for numerical speed up\n",
    "\n",
    "@jit(nopython=True)\n",
    "def cova2(x1, y1, x2, y2, nst, c0, pmx, cc, aa, it, ang, anis, rotmat, maxcov):\n",
    "    \"\"\"Calculate the covariance associated with a variogram model specified by a\n",
    "    nugget effect and nested variogram structures.\n",
    "    :param x1: x coordinate of first point\n",
    "    :param y1: y coordinate of first point\n",
    "    :param x2: x coordinate of second point\n",
    "    :param y2: y coordinate of second point\n",
    "    :param nst: number of nested structures (maximum of 4)\n",
    "    :param c0: isotropic nugget constant (TODO: not used)\n",
    "    :param pmx: TODO\n",
    "    :param cc: multiplicative factor of each nested structure\n",
    "    :param aa: parameter `a` of each nested structure\n",
    "    :param it: TODO\n",
    "    :param ang: TODO: not used\n",
    "    :param anis: TODO\n",
    "    :param rotmat: rotation matrices\n",
    "    :param maxcov: TODO\n",
    "    :return: TODO\n",
    "    \"\"\"\n",
    "    EPSLON = 0.000001\n",
    "\n",
    "    # Check for very small distance\n",
    "    dx = x2 - x1\n",
    "    dy = y2 - y1\n",
    "\n",
    "    if (dx * dx + dy * dy) < EPSLON:\n",
    "        cova2_ = maxcov\n",
    "        return cova2_\n",
    "\n",
    "    # Non-zero distance, loop over all the structures\n",
    "    cova2_ = 0.0\n",
    "    for js in range(0, nst):\n",
    "        # Compute the appropriate structural distance\n",
    "        dx1 = dx * rotmat[0, js] + dy * rotmat[1, js]\n",
    "        dy1 = (dx * rotmat[2, js] + dy * rotmat[3, js]) / anis[js]\n",
    "        h = math.sqrt(max((dx1 * dx1 + dy1 * dy1), 0.0))\n",
    "        if it[js] == 1:\n",
    "            # Spherical model\n",
    "            hr = h / aa[js]\n",
    "            if hr < 1.0:\n",
    "                cova2_ = cova2_ + cc[js] * (1.0 - hr * (1.5 - 0.5 * hr * hr))\n",
    "        elif it[js] == 2:\n",
    "            # Exponential model\n",
    "            cova2_ = cova2_ + cc[js] * np.exp(-3.0 * h / aa[js])\n",
    "        elif it[js] == 3:\n",
    "            # Gaussian model\n",
    "            hh = -3.0 * (h * h) / (aa[js] * aa[js])\n",
    "            cova2_ = cova2_ + cc[js] * np.exp(hh)\n",
    "        elif it[js] == 4:\n",
    "            # Power model\n",
    "            cov1 = pmx - cc[js] * (h ** aa[js])\n",
    "            cova2_ = cova2_ + cov1\n",
    "    return cova2_\n",
    "\n",
    "@jit(nopython=True)\n",
    "def setup_rotmat(c0, nst, it, cc, ang, pmx):\n",
    "    \"\"\"Setup rotation matrix.\n",
    "    :param c0: nugget constant (isotropic)\n",
    "    :param nst: number of nested structures (max. 4)\n",
    "    :param it: TODO\n",
    "    :param cc: multiplicative factor of each nested structure\n",
    "    :param ang: TODO\n",
    "    :param pmx: TODO\n",
    "    :return: TODO\n",
    "    \"\"\"\n",
    "    PI = 3.141_592_65\n",
    "    DTOR = PI / 180.0\n",
    "\n",
    "    # The first time around, re-initialize the cosine matrix for the variogram\n",
    "    # structures\n",
    "    rotmat = np.zeros((4, nst))\n",
    "    maxcov = c0\n",
    "    for js in range(0, nst):\n",
    "        azmuth = (90.0 - ang[js]) * DTOR\n",
    "        rotmat[0, js] = math.cos(azmuth)\n",
    "        rotmat[1, js] = math.sin(azmuth)\n",
    "        rotmat[2, js] = -1 * math.sin(azmuth)\n",
    "        rotmat[3, js] = math.cos(azmuth)\n",
    "        if it[js] == 4:\n",
    "            maxcov = maxcov + pmx\n",
    "        else:\n",
    "            maxcov = maxcov + cc[js]\n",
    "    return rotmat, maxcov\n",
    "\n",
    "def n_effective(df,xcol,ycol,seed,nreal,vario):\n",
    "    \"\"\"Calculate the number of effective data from spatial locations and spatial continuity model\n",
    "    Used in bootstrap to account for spatial continuity, use n effective instead of number of data\n",
    "    :param df: source DataFrame\n",
    "    :param xcol: column with the X locations\n",
    "    :param ycol: column with the Y locations\n",
    "    :param seed: random number seed for the random sampling\n",
    "    :param nreal: number of realizations to sample the variance of the average \n",
    "    :param vario: variogram model as a dictionary, see the GeostatsPy Package's GSLIB.make_variogram() function \n",
    "    :return: n_eff as effective number of data\n",
    "    \"\"\" \n",
    "\n",
    "# Set constants\n",
    "    np.random.seed(seed)\n",
    "    PMX = 9999.9\n",
    "    \n",
    "# load the variogram\n",
    "    nst = vario['nst']\n",
    "    cc = np.zeros(nst); aa = np.zeros(nst); it = np.zeros(nst)\n",
    "    ang = np.zeros(nst); anis = np.zeros(nst)\n",
    "    c0 = vario['nug']; \n",
    "    cc[0] = vario['cc1']; it[0] = vario['it1']; ang[0] = vario['azi1']; \n",
    "    aa[0] = vario['hmaj1']; anis[0] = vario['hmin1']/vario['hmaj1'];\n",
    "    if nst == 2:                                   # include 2nd structure if present (optional)\n",
    "        cc[1] = vario['cc2']; it[1] = vario['it2']; ang[1] = vario['azi2']; \n",
    "        aa[1] = vario['hmaj2']; anis[1] = vario['hmin2']/vario['hmaj2'];\n",
    "    \n",
    "# Set up the rotation matrix\n",
    "    rotmat, maxcov = setup_rotmat(c0,nst,it,cc,ang,PMX)\n",
    "    \n",
    "# Load the data\n",
    "    nd = len(df)\n",
    "    x = df[xcol].values\n",
    "    y = df[ycol].values\n",
    "    \n",
    "# Calculate Symmetric Covariance Array - assuming variogram with spherical structure with range specified\n",
    "    cov = np.zeros((nd,nd))\n",
    "    var_range = 100.0\n",
    "    for i in range(0, nd):\n",
    "        x1 = x[i]; y1 = y[i]\n",
    "        for j in range(0, nd):\n",
    "            x2 = x[j]; y2 = y[j]\n",
    "            cova = cova2(x1, y1, x2, y2, nst, c0, PMX, cc, aa, it, ang, anis, rotmat, maxcov)\n",
    "            cov[i,j] = cova\n",
    "            \n",
    "# Lower and upper deconvolution            \n",
    "    P, L, U = scipy.linalg.lu(cov) \n",
    "    \n",
    "# Build realization and calculate the average    \n",
    "    average_array = np.zeros(nreal)\n",
    "    rand = np.zeros((nd)) \n",
    "    for l in range(0, nreal):\n",
    "        rand = np.random.normal(loc = 0.0, scale = 1.0, size = nd)\n",
    "        realization = np.matmul(L,rand)\n",
    "        average_array[l] = np.average(realization)\n",
    "        \n",
    "# Back out the number of effecitve data useing the standard error in the average\n",
    "    var_average = np.var(average_array)\n",
    "    n_eff = max(min(1.0/var_average, nd),1.0)    # filter n effective less than 1.0 or greater than number of data\n",
    "    \n",
    "    return n_eff\n",
    "    "
   ]
  },
  {
   "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": 8,
   "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.  \n",
    "\n",
    "* Add a null feature for plotting the data locations, data values are not needed in this calculation / workflow for number of effective data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"12_sample_data.csv\")                      # read a .csv file in as a DataFrame\n",
    "df['null'] = np.zeros(len(df))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's drop some samples so that we increase the variations in bootstrap samples for our demonstration below.  \n",
    "\n",
    "* We will resample 25% of the original data file.  \n",
    "\n",
    "* I use a random number seed so that we all get the same sample to work with.  \n",
    "\n",
    "```python \n",
    "df.sample(frac = 0.1, random_state = 73073)\n",
    "```\n",
    "\n",
    "where *frac* is the fraction to select and *random_state* is the random number seed.\n",
    "\n",
    "To get different random 25% subset of the data, just drop the 'rand_state' parameter or change the value from '73073' to any other integer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using 48 number of samples\n"
     ]
    }
   ],
   "source": [
    "df = df.sample(frac = 0.1, random_state = 73073)                  # 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 plot the data locations below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfkAAAHtCAYAAAAa8QdpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xl4lOXd9vHvb7bsYQcBsWjFHRcMagVc0LpWVBQFV3Cjllofq9S2+lh9qrUudV+pilqtWqxb1VbRan3rVgMoKlZBUURB2ZPMJJntev/IgIkga5Jr5s75OY4cmbnmnnvODBPO3Ls55xAREZHgCfkOICIiIm1DJS8iIhJQKnkREZGAUsmLiIgElEpeREQkoFTyIiIiAaWSF5HAMLMTzex53zlE8oVKXmQTmNmnZnZgO7zOpWb2wHpkqTezWjNbbmavmdmPzWy9fs/NrL+ZOTOLbEJOZ2Zbb+zzN/C1VsvrnHvQOXdQe7y+SCFQyYsEyxHOuQrge8DvgQuBu/1GEhFfVPIircTMxprZv83sWjNbZmZzzezQZo+/bGZXmtl/zGyFmT1pZl1zj+1nZvO/Nb9PzexAMzsE+DVwvJnVmdk768rinFvhnHsKOB441cx2ys3zcDObYWY1Zva5mV3a7Gmv5L4vz73OD8zs+2b2TzNbYmaLzexBM+u8Ee9NyMwuNrPPzOxrM7vfzDo1e3xobs3D8lyusRuZd6yZ/bvZfPc2s7dy7/dbZrZ3s8deNrPfmtmrubUfz5tZ99xjxWb2QO7nXp57bq8N/blFfFPJi7SuPYEPge7A1cDdZmbNHj8FOA3oA6SBm9Y1Q+fcP4DfAY8458qdc7usbxjn3H+A+cCw3FA8l6EzcDhwtpkdlXtsn9z3zrnXeR0w4Mpc3u2BfsCl6/v6zYzNfe0PbAWUA7cAmNkWwN+Bm4EewK7A2xuZd5XcH1DP0PQedwOuA54xs27NJjsBGAf0BGLABbnxU4FOuZ+3G/BjoH4jfm4Rr1TyIq3rM+fcH51zGeA+oDfQfAnwT86595xzceB/gePMLNzGmb4EugI45152zr3rnMs652YCDwH7ftcTnXNznHNTnXONzrlFNBXld06/FicC1znnPnHO1QG/AkbntqefCLzgnHvIOZdyzi1xzr29MXm/5XBgtnPuT865tHPuIeC/wBHNppnsnPvIOVcP/IWmPzAAUjSV+9bOuYxzbppzrmYjfm4Rr1TyIq1r4cobzrlE7mZ5s8c/b3b7MyBK01J/W+oLLAUwsz3N7CUzW2RmK2haQv3O1zeznmb2sJl9YWY1wAMbmbcPTT/vSp8BEZr+AOoHfPwdr79Bedfxmitft2+z+wub3U7wzb/Vn4DngIfN7Eszu9rMouv5uiJ5QyUv0r76Nbu9BU1LjItpWi1duvKB3NJ9j2bTbtTlIs1sME2ltnI79Z+Bp4B+zrlOwB00rZL/rte4Mje+s3OuEjip2fQb4kuadgZcaQuaNld8RdMfPt//judtaN61vebK1/1iXWFzaxQuc87tAOwN/IimzQYiBUUlL9K+TjKzHcysFPg/4NHcqv2PgOLcjmZR4GKgqNnzvgL6b8DhcJVm9iPgYeAB59y7uYcqgKXOuQYz24OmbdIrLQKyNG0zp9n0dTTt3NYXmLgeLx/L7bi28itM02r288xsSzMr55t9DNLAg8CBZnacmUXMrJuZrVxtvqF5m3sW2MbMTsjN93hgB+Dpdf0AZra/mQ3MZa+h6Y+xzHr87CJ5RSUv0r7+BNxL02riYuBn0LQ3PPAT4C6aljTjNO0wt9KU3PclZjZ9LfP/m5nV0rR0fBFN29DHNXv8J8D/5aa5hKbt0OQyJIArgFdze5TvBVwGDAJW0LQT22Pr8TO+T9NOaiu/xgH35H72V4C5QANwTu515wGHAefTtFnhbWDlzoUbmpdmjy+haQn8fGAJ8AvgR865xevxM2wGPEpTwX8A/IumTRUiBcWc26i1gCKygczsZZqWqu/ynUVEOgYtyYuIiARUm5W8md2TO+nFe83GuprZVDObnfveJTduZnaTmc0xs5lmNqjZc07NTT/bzE5tq7wiIiJB02ar681sH5p22LnfObfybFtX07QTze/N7JdAF+fchWZ2GE3b5w6j6WQiNzrn9sydzKIaqKJpT9ppwO7OuWVtElpERCRA2mxJ3jn3Crljc5s5kqYThJD7flSz8ftdkzeAzmbWGzgYmOqcW5or9qnAIW2VWUREJEjae5t8L+fcAoDc95658b60PEnI/NzYd42LiIjIOmz0JSVb2ZpOruHWMr76DMzOAs4CKCsr23277bZrvXQiIiLtYNq0aYudcz3WPeX6ae+S/8rMejvnFuRWx3+dG59PyzOBbU7T2armA/t9a/zlNc3YOTcJmARQVVXlqqurWzd5AchmszzyyCNMnnwvjY0NHHbYYUyYMIHy8vJ1P1lERLwzs2+finmTtPfq+qdouroTue9PNhs/JbeX/V7Aitzq/OeAg8ysS25P/INyY7IGl112GZdd9lsWLVpMXV2Ce+65l+OPP550Ou07moiIeNCWh9A9BLwObGtm883sdOD3wA/NbDbww9x9aDr95CfAHOCPNJ3lCufcUuC3wFu5r//Ljcm3LF68mCeeeJKysnLC4QhmISoqKvn003m8/PLLvuOJiIgHbba63jk35jseOmAN0zpgwnfM5x6aTokpa/HJJ5+QzTpaXroc0ukMM2bM4MADD/SUTEREfNEZ7wKiX79+hMOr76cYDofQTogiIh2TSj4gevfuzdChQ6mrqyWbzeKco66ujp49e3DQQQf5jiciIh6o5APk+uuv58wzT6eoKEIo5DjooAOYMuUvFBUVrfvJIiISOIG8Cl1HPYROREQKm5lNc85Vtdb8tCQvIiISUCp5ERGRgFLJi4iIBJRKXkREJKBU8iIiIgGlkhcREQkolbyIiEhAqeRFREQCSiUvIiISUCp5kYD54IMPOProkey2227suede3HrrrWSzWd+xRMSDNrvUrIi0vy+//JITTzyRVCpLUVERqVSGW265lbq6Oi688ELf8USknWlJXiRAJk+eTCLRsOqiRKFQiNLScv7yl7/Q2NjoOZ2ItDeVvEiAfPTRR8RisRZjZkY261ixYoWnVCLii0peJED23ntvkslki7FsNkssFqVbt26eUomILyp5kQAZM2YMvXr1pKZmBdlslsbGRhoaEkycOJFwOOw7noi0M5W8SIBUVlby+OOPMXbsKXTpUsnOO+/I3XffxbHHHus7moh4YM453xlaXVVVlauurvYdQ0REZIOY2TTnXFVrzU9L8iIiIgGlkhdpBZlMhjvuuINhw4axxx57cO655/L111/7jiUiHZxKXqQVXHzxxdx4483E4w1kMvD88y9w7LGjSCQSvqOJSAemkhfZREuXLuXvf/8HpaVlhEIhzIyysnIWLVrMM8884zueiHRgKnmRTTR//nycazrpzLe99957HhKJiDTRuetFNtH3vvc9zMA516LozWDXXXf1mCyYnHO88cYbPPHEE5SXlzN69GgGDBjgO5ZIXtKSvMgm6tSpE8ccM5J4vI5MJoNzjtraGvr06c1hhx3mO16gOOe48MILOeOMM/nb357hgQf+zMiRI5kyZYrvaCJ5SSUv0gouuugifvOb/6V7964UF0cZPfo4pkyZsupCMdI63n33XZ599u8UF5dSWlpGRUUl0WgxV175e+3kKLIGWl0v0gpCoRBjxoxhzJgxvqME2j//+U+SyTTFxd9sFgmFQjQ2NvLee++xxx57eEwnkn+0JC8iBaN79+6EQqvv4BgKhejUqZOHRCL5TSUvIgXj8MMPp7S0hFQqtWoskUjQt28fttlmG4/JRPKTSl5ECkaXLl344x8nUVFRRjqdJJ1Osu22WzN58uQ1HsIo0tFpm7yIFJTdd9+dV175F5999hnFxcX07t3bdySRvKWSF5GCEwqF2HLLLdv9dRsaGnj11Vepr6/nBz/4Ad26dWv3DCIbQiUvIrIeZs6cyRlnnEFtbZx0Ok1paTETJ07kpJNO8h1N5Dtpm7yIyDpkMhl+/OMf09iYpqysnE6dOhOJxLjqqquZN2+e73gi30klLyKyDu+//z51dQmi0eiqMbMQjY1J/va3v3lMJrJ2KnkREZGA0jZ5Eckrc+fO5f7772fBggUccsghHH744S2WoH3YcccdKS8vJR5vWJXFuSxFRTGOOOIIr9lE1kZL8iKSN1566SWOPPIoHnroEf7979f41a8uYsyYMS1OfuNDOBzmjjvuIBaLEI/XUVOzgnQ6xYUX/oItttjCazaRtTHnnO8Mra6qqspVV1f7jiEiGyCbzTJ06DDq6xsJh8OrxuPxWq6++qq8WGJuaGjgtddeo76+nr322kuH0EmrM7Npzrmq1pqfVteLSF5YsGAB9fX1hMMtV82HwxGeffbZvCj54uJihg8f7juGyHpTyUvg1NfXM23aNGKxGIMGDSIS0ce8EFRUVLCmNYvpdJo+ffp4SCRS+PS/nwTK1KlT+eUvf0lDQxIzo6KinLvu+iM77rij72iyDpWVlfzgB3vx8suvUFZWDjQdnx6LRTn55JM9pxMpTNrxTgJj6dKlXHDBRLJZo6SklOLiEhKJes4880zS6bTveLIe/vCHPzB8+H5kMkkymRSlpUXcdNON9O/f33MykcKkJXkJjOeee476+gYqK7+5rng0GqO2to7p06ezxx57eEwn66O0tJRbb72Vmpoaamtr6d27N6GQlkVENpZKXgKjsbFxjdt0nXPeD8GSDVNZWUllZaXvGCIFT38iS2D88Ic/pLi4qEXRZzIZioqK2H333T0mExHxQyUvgdG3b1/OO+9/SCYbWLFiOTU1K3AuwzXXXE1xcbHveCIi7U4nw5HA+fLLL3nhhReIxWIcdNBBdO3a1XckEZH1opPhiKxDnz59OOWUU3zHEBHxTqvrRUREAkolLyIiElAqeRERkYBSyYuIiASUSt6TuXPncu655zJ8+HAmTJjA7NmzfUcKhGQyyV133cUhhxzCiBEj+Otf/0o2m/UdS0TECx1C58GcOXMYNWoUDQ1JSkvLqK9PEI1GePDBB9hpp518xytYzjlOOOEEpk9/m9LSMpzL0tDQwMiRR/G73/3OdzwRkXVq7UPotCTvwZVXXkkymaasrBwzo7S0jGzW8dvf/tZ3tIL21ltv8e6771NRUUk4HCYSiVJWVs4zzzzLwoULfccTEWl3KnkPZs2aRUlJaYuxoqJi5s791E+ggHjnnXdobGxsMWZmpFJpPvzwQ0+pRET8Ucl70Lt3b5LJZIuxdDpFt246M9um6N+/P7FYdLXxaDRCnz59PCQSEfFLJe/B+eefTzqdXHWN80wmQzLZyHnnnec5WWHbb7/96NWrJ/F4HOcczjnq6mrZYYftGTBggO94IiLtTiXvwZAhQ7jxxhvo1q0zmUySTp0quPrqqzjooIN8Ryto0WiUhx9+mH32GUI2m8K5DCNHHsXdd9/tO5qIiBfau15ERCRPaO96ERERWS8qeRERkYBSyYuIiASUSl5ERCSgVPIiIiIBpZIXEREJKJW8iIhIQKnkRUREAkolLyIiElAqeRERkYBSyYuIiASUSl5ERCSgVPIiIiIBpZIXEREJKC8lb2bnmdn7ZvaemT1kZsVmtqWZvWlms83sETOL5aYtyt2fk3u8v4/MIiIihabdS97M+gI/A6qcczsBYWA0cBVwvXNuALAMOD33lNOBZc65rYHrc9OJiIjIOvhaXR8BSswsApQCC4DhwKO5x+8DjsrdPjJ3n9zjB5iZtWNWERGRgtTuJe+c+wK4FphHU7mvAKYBy51z6dxk84G+udt9gc9zz03npu/27fma2VlmVm1m1YsWLWrbH0JERKQA+Fhd34WmpfMtgT5AGXDoGiZ1K5+ylse+GXBuknOuyjlX1aNHj9aKKyLSYSSTSZ5//nkeeOABPvroI99xpBVEPLzmgcBc59wiADN7DNgb6GxmkdzS+ubAl7np5wP9gPm51fudgKXtH1tEJLg+++wzTjzxRJYsWUYqlaKkpJhDDjmYq666ilBIB2IVKh//cvOAvcysNLdt/QBgFvAScGxumlOBJ3O3n8rdJ/f4P51zqy3Ji4jIxvvZz37G8uU1lJWV07lzF2KxYp5++ln++c9/+o4mm8DHNvk3adqBbjrwbi7DJOBC4OdmNoembe53555yN9AtN/5z4JftnVlEJMiWL1/OZ5/No7i4ZNWYmREOh3nkkUc8JpNN5WN1Pc653wC/+dbwJ8Aea5i2ARjVHrlECkVNTQ2pVIpu3VbbB1Vkg33X6njnIBqNtnMaaU3a0CJSQBYvXsyJJ57IPvvsywEHHMihhx7K7NmzfceSAldZWcn2229HIpFYNeacw7ksJ598ssdksqlU8iIFwjnH2LHjmDFjJpFIjEgkxvz5CzjppJNa/OcssjFuvPFG+vXrS2NjA/F4Lel0I6edNpYf/OAHvqPJJvCyul5ENtwHH3zAvHnzKC0tXTVWVFRETU0dzz//PEcdddRani2ydj179uTvf3+Wd955h8WLF7PLLrvQs2dP37FkE6nkRQrEkiVLyGSyq42n02kWLlzoIZEETSgUYrfddvMdQ1pRoFfXL1y4kF//+tfsv//+nHLKKUyfPt13JJGNNnDgQKLRMN8+grS4uIghQ4Z4SiUSbO+++y6nn346++23PxdccAHz58/3HWmDBLbkv/76a4488igee+wJamriTJ/+DieffApTp071HU1ko3Tu3Jnx48fT0JAgkUjQ0FBPIlHH8OH7s9NOO/mOJxI4r732GmPGnMAbb7xFbW2cp59+lqOPPrqgit6CeF6Zqqoqd+SRR/Lggw9RXl6xajydTtO5cyUvvfRPdI0bKVTTp0/n/vvvp76+nlGjRjF8+HCdkUykDRxyyKEsWPBVi8MI6+pqOeKIw7nmmmva5DXNbJpzrqq15hfYbfJvvvkmJSWlLcYikQgrViwnlUoRi8U8JRPZNIMGDWLQoEG+Y4gE3ldffbXaeQLKysqprq72lGjDBfbP/wEDBtDQ0NBiLJvNEosV6eQOIiKyTuXlZWQymRZjDQ0NfP/73/eUaMMFtuQnTJhANBommUwCkM1maGhIMH78WVpVLyIi6/Szn/2Mhob6VUWfSqUwc5x33nmek62/wJb81ltvzeTJ97DFFpuTyaQoLS3mV7/6JePGjfMdTURECsCoUaP4v/+7lIqKUjKZJH36bMYf/ziJHXfc0Xe09RbYHe8KaZuJiIgItP6Od4FdkhcREenoVPIiIiIBpZIXEREJKJW8iIhIQKnkRUREAkol3wrmzJnD+PHjGTp0KCeeeCIzZszwHWmVt99+m5NOOomhQ4dy5plnMnv2bN+RRESknegQuk308ccfc8wxx5JMpikpKSGZbMS5LJMm3cnee+/dLhm+yxtvvMEZZ5yJWYhYrIiGhnqi0TBTpkxhwIABXrOJiMjqdAhdnrn22mtJJtOUlpZiZhQVFRMOR7jiiit8R+Pyy68gHI5QVFSMmVFSUkoymeHaa6/1HU1ERNqBSn4Tvffee5SUlLQYi0ZjLFz4ladE31i4cAHRaMsL8ZSUlPD+++97SiQiIu1JJb+JttpqKxobW14IJ51O06VLF0+JvtG1a1fS6XSLscbGRrbccktPiURkfcyfP5/bbruNm266iU8++cR3HClgKvlNdMEFF2DGqgvhpNNpkskGLrjgfM/J4PzzzyeZbCCdTgErM2aZOHGi32Ai8p0effRRDj30MG644SZuueU2Row4kjvuuMN3LClQKvlNNHDgQO6552623HILzLL07NmN66+/jkMOOcR3NA4++GCuv/46evbsgVmW/v37cc89d7Pzzjv7jiYia1BTU8Pll19BNFpEZWUnKis7UVRUwq233saXX37pO54UoIjvAEEwePBgnnzySd8x1ujggw/m4IMP9h1DRNbDW2+9RWNjkrKyb/alMTPq6xv417/+xZgxYzymk0KkJXkRkTxRUlJCOBxebTwcDlNWVuYhkRQ6lbyISJ4YPHgwlZXlq/bxAUinU5SXlzJ8+HCPyaRQqeRFRPJENBpl8uTJdO5cQSrVSDLZSHFxjDvvvJPy8nLf8aQAaZu8iEge2XbbbXn55ZeZNWsWmUyGHXfckUhE/1XLxtEnR0Qkz4RCIXbaaSffMSQAtLpeREQkoFTyIiIiAaWSFxERCSiVvIiISECp5EU8aWho4De/+Q277747u+66G2PHjtWpS0WkVankRTz56U9/yl/+8igQJhKJ8Z//VHPcccdRX1/vO5qIBIRKXsSDL774grfeqqasrBwzA6C0tIwlS5bx97//3XM6EQkKlbyIB1988QXZrFttPJt1fPzxxx4SiUgQBb7k6+rqeOedd1iwYEGrznfevHnMnDmThoaGVp2vdAwDBgwgHA7hXMuij0TCDB482FMqEQmaQJ/x7vbbb2fSpElkMg4zGDy4iptvvpmSkpKNnueyZcs466yzmD17Ds45IpEwv/rVrzj22GNbMbkEXZcuXTj++OO4//4/EYsVEwqFSCTi7LDDdgwbNsx3PBEJCPv2kkQQVFVVuT/84Q+cffZPKCkpW7XNs66ulqOPPpIrr7xyo+d90kknMW3a26su+5jNZkkmG3j00Slst912rZJfOgbnHM888wx333038Xico48+mrFjx27SH6EiUtjMbJpzrqrV5hfUkt9+++15++2ZFBUVrxp3zpHNppg2bRqh0IZvqVi2bBn77z+cSCTWYryuro5jjjmKK664YpOzi4hIx9XaJR/YbfI1NbWEw+EWY2aGc5DJZDZqnt91aFM4HGbFihUbNU8REZG2EtiSP+qoI1cr5YaGer7//a2IRqMbNc/evXvTuXMnUqlUi3Hnshx55JEbnVVERKQtBLbkTzjhBAYO3JH6+gR1dXXU1dVSVlbKNddcs9HzNDOuu+46wmGora2hrq6WhoYEw4YN4YADDmjF9CIiIpsusNvkq6uryWQyvPzyy7z66qtstdVWjBgxgsrKyk2e/5IlS3jiiSdYsGABP/zhDxk8ePBGbeMXERFpTjverYeVJS8iIlJItOOdiIiIrBeVvIiISECp5EVERAJKJS8iIhJQKnkRaVNz5szhzTffpKamxncUkQ4n0BeoERF/li1bxumnn87HH39CNpslEgkzYcIEzjjjDN/RRDoMLcmLSJv4+c9/zgcffEQ0WkRRUQmhUJTrr7+B6dOn+44m0mGo5EWk1dXV1TFjxjdXa4SV144wJk+e7DGZSMei1fUi0uq+fX2HlcLhEPF4vM1f3znHrFmzmDNnDltuuSUDBw5cdclpkY5EJS8ira5Lly707r0ZX3yxgFisaNV4Op3mmGOOadPXbmxs5IwzzuCdd2aSSNRTUlLMdttty3333UdpaWmbvrZIvtHqehFpEzfeeCPFxUXU1dVSW1tDQ0M9++wzlEMPPbRNX/f222/nrbemEYsV07lzF4qKSpg58z2uu+66Nn1dkXykc9eLSJupq6vj2Wef5YsvvmDfffdlt912a/PV5vvssy91dYkWF41yzhGLhXn99dfb9LVFNlVrn7teq+tFpM2Ul5dz3HHH+Y4h0mGp5EUkUEaOPJo77phEeXnFqrG6ulpOPvlEj6k2zPTp03nwwQdJp9Mcd9xxbLXVVjzwwAPMmTOHYcOGMXLkSO1fIOtFq+tFJFCSySSnn376qh3vSktL2HbbbQpmx7ubb76ZO+64E+eaDjtsaKgnHo/TqVNnotEY6XSKfv368thjj1FRUbHuGUpB0ep6EZG1iMVi3H///cyaNYuPP/6Y/v37F8whdIsXL+auu+6muLh0Vd7587+gsbGeHj16UlRUDBQzb9587r33Xs455xy/gSXvae96EQkcM2PHHXdkxIgR7LzzzgVR8AAzZswgmUyuyutc0zkHnDPq6upWTVdaWsY//vGcr5hSQFTyIiJ5okuXLkQi36xgNWv6g8XMWoynUkl69erpI6IUGJW8iEieGDRoED179qC+PrFqrKKiAjO3akfCbDaLcxl+8pOf+IopBUQlLyKSJ0KhEH/605/YZputSaeTZDJJdtllJ047bRyQJZNJEYkYl112GVVVrbZvlgSYdrwTEckjffr04fHHH2fx4sVkMhl69eoFQDweZ9myZWy22WYtVt2LrI0+KSIieah79+4t7peVlbW4qp/I+gjk6vpZs2Zx+OE/0nWrRUSkQwtkyZsZn3/+Baeddhpz5871HUdERMSLQJY8NJ0Qo6EhyZ133uk7ioiIiBeBLXmA4uJi5syZ4zuGiIiIF4Eu+fr6evbee2/fMURERLwIZMk7B7W1NXTv3o2xY8f6jiMiIuJFIEs+Eglz4oljePLJJ+jatavvOCIiIl7oUrMiIiJ5orUvNetlSd7MOpvZo2b2XzP7wMx+YGZdzWyqmc3Ofe+Sm9bM7CYzm2NmM81skI/MrW3u3Lmce+65DB8+nAkTJjB79mzfkUREJGB8ra6/EfiHc247YBfgA+CXwIvOuQHAi7n7AIcCA3JfZwG3t3/c1jVnzhxGjhzJ88+/wIoVdbz00r849thRvPfee76jiYhIgLR7yZtZJbAPcDeAcy7pnFsOHAncl5vsPuCo3O0jgftdkzeAzmbWu51jt6orr7ySZDJNWVk5ZkZpaRnZrOO3v/2t72giIhIgPpbktwIWAZPNbIaZ3WVmZUAv59wCgNz3lRdL7gt83uz583NjBWvWrFmUlJS2GCsqKmbu3E/9BBIRkUDyUfIRYBBwu3NuNyDON6vm18TWMLba3oJmdpaZVZtZ9aJFi1onaRvp3bs3yWSyxVg6naJbNx0J4JwjmUySzWZ9RxERKXg+Sn4+MN8592bu/qM0lf5XK1fD575/3Wz6fs2evznw5bdn6pyb5Jyrcs5V9ejRo83Ct4bzzz+fdDpJOp0GIJPJkEw2ct5553lO5tcrr7zCAQccwB577Mmee+7JLbfcorIXEdkE7V7yzrmFwOdmtm1u6ABgFvAUcGpu7FTgydztp4BTcnvZ7wWsWLlav1ANGTKEG2+8gW7dOpPJJOnUqYKrr76Kgw46yHc0b2bNmsVcefWGAAAgAElEQVSECT9l6dIVRCIxslm49dbbue2223xHExEpWF6OkzezXYG7gBjwCTCOpj84/gJsAcwDRjnnlpqZAbcAhwAJYJxzbq0Hwes4+cIzYcIEXnrpFUpLv9lXIZvNEg4b//nPmzR9DEREgq21j5OPtNaMNoRz7m1gTT/EAWuY1gET2jyUeDVv3jyKiopajIVCIVKpJKlUilgs5imZiEjhCuRpbaXwDB06lEQi3mJs5c6IKngRkY2jkpe8cOaZZ9KjRzdqa2vIZNIkEnEymTSXXXaZ72giIgVLJS95oWvXrjz11FOcdtpY+vXrywEH7M+UKX9h6NChvqOJiBQsL9vkRdaka9euXHjhhb5jiIgEhko+jzU2NvLvf/+beDzOnnvuSa9evXxHEhGRAqKSz1OzZs1i3LjTqKmpJZ1OU1JSzIQJP2H8+PG+o4mISIHQNvk8lM1m+fGPf0x9fSNlZeV06tSZaLSIm2++hY8++sh3PBERKRAq+Tw0Z84campqWxw6ZmakUmkef/xxj8lERKSQqOTzlI8zEYqISLCo5PPQ1ltvTefOnVtcqc45RywWYeTIkR6TiYhIIVHJ56FQKMSkSXdSWlpMIhGnpmYF6XSS//mf/2HAgAG+44mISIHQ3vV5atttt+WVV/7F66+/Tl1dHXvuuSfdu3f3HUtERAqISj6PxWIx9t13X98xRNpUNpvl008/JRaLsfnmm/uOI7JOqVSKuXPn0rVr17xf+FLJS95YvHgxd955J6+++irf+973OOecc9hhhx18x5I2NHPmTH76059SU1MLQN++fbjzzjtV9pK3nnzySS6//HJSqTQAu+8+iJtuuomysjLPydbMy/Xk25quJ194lixZwogRI1i8eCmlpWUkk42YwW233co+++zjO560gbq6Ovbbb39SqQzRaBSAhoZ6evXqwdSpUwmFtMuQ5JcPPviAUaOOIxYrXvX5jMfjDB++L7feemurvEZrX09ev0WSFyZNmsTixUupqKgkHA5TUlJKJBLj0ksv9R1N2sjzzz9PXV18VcEDFBeXsGjRYmbOnOkxmciaTZ48mUwm2+IP0LKyMl599TXi8fhanumPSl7ywquvvkppacvVXZFIhKVLl7U4lFCCY+nSpWQymdXGM5kMNTU1HhKJrN3ixYuJRFbfym0GiUTCQ6J1U8lLXthqq61obGxoMZbNZonFoi2W9CQ49t13X0pKilqc+CmbzRKNRtl11109JhNZs8MPP5x0Ot1iLJVK0blz57zdAU8lL3nhJz/5CeFwiFQqBTT9Z19fn+C0007DzDynk7YwYMAAjjnmGBoaEsTjddTW1pBMNjBx4gVUVlb6jieymhEjRrDLLgOJx+tIJBLU1tYQCjmuueaavP1/SjveSd547bXXuPTSS/n6668pKipi3LhxjB8/Pm9/eWTTOeeYPn06jz/+OEVFRRx//PFss802vmOJfKdMJsOLL77I1KlT2XzzzRk9enSrXga8tXe8U8lL3slkMoRCIZW7iHQ4rV3yOk5e8k44HPYdQUQkELRNXkREJKBU8iIiIgGlkhcREQkolbyIiEhArXPHOzMrBn4EDAP6APXAe8Azzrn32zaeiIiIbKy1lryZXQocAbwMvAl8DRQD2wC/z/0BcL5zTieaFhERyTPrWpJ/yzl36Xc8dp2Z9QS2aN1IIiIi0hrWWvLOuWfW8fjXNC3di4iISJ5Zr5PhmFkVcBHwvdxzDHDOuZ3bMJuIiIhsgvU9492DwETgXSDbdnFERESktaxvyS9yzj3VpklERESkVa1vyf/GzO4CXgQaVw465x5rk1QiIiKyyda35McB2wFRvlld7wCVvIjIJkin07z44ou88MIL9O3bl9GjR7PZZpv5jtXuGhoaWLhwIT179qS0tNR3nMBY35LfxTk3sE2TiIh0MKlUipNPPpl33nkXsxCZTJr77ruPP/7xj1RVtdrVRvOac45bbrmFe++9l2w2CxhjxozmggsuIBTSSVk31fq+g2+Y2Q5tmkREpIN56qmneOeddykrK6e0tJSKikqyWWPixIk453zHaxdPPPEEt912OxAmHI4RCkWYPPleHnjgAd/RAmF9S34o8LaZfWhmM83sXTPTWe5ERDbBM888QyTScoVqNBpl+fLlLF682FOq9nXnnXdSVFSCmQFgZhQXl3LPPfd4ThYM67u6/pA2TSEi0gF1796ddDpNUVHLcefoMNula2pqCIfDLcZCoRCJRMJTomBZ65K8mZUDOOc+W9NX82lERGTDjBs3jnA4lNsW3SQejzNkyN6UlZV5TNZ+hg4dSjxe12IskYiz++67e0oULOtaXf+kmf3BzPYxs1WfODPbysxON7Pn0FK+iMhG2X777bniissJhyGdTpJOJ9l77z25+uqrfUdrNxdccAHdu3elrq6WhoZ66upq6dy5ExdffLHvaIFg69q5w8wOA04EhgBdgDTwIfAMcLdzbmFbh9xQVVVVrrq62ncMEZH1kkql+PTTT+nSpQvdu3f3Hafd1dTU8Ne//pVp06YxcOBAjj/+eDp37uw7lhdmNs0512qHVqyz5AuRSl5ERApRa5e8DkIUEREJKJW8iIhIQK1r7/pnzax/+0QRERGR1rSuJfl7gefN7CIzi7ZDHhEREWklaz0ZjnPuL2b2DHAJUG1mf6LZ9eSdc9e1cT4RERHZSOtzxrsUEAeKgAqalbyIiIjkr7WWvJkdAlwHPAUMcs7pPIMiIiIFYl1L8hcBo5xz77dHGBEREWk969omP6y9goiIiEjr0nHyIiIiAaWSFxERCSiVvIiISECp5EVERAJKJS8iIhJQKnkREZGAUsmLiIgElEpeREQkoFTyIiIiAaWSFxERCSiVvIiISECp5EVERAJKJS8iIhJQKnkREZGAUsmLiIgElEpeREQkoFTyIiIiAaWSFxERCSiVvIiISECp5EVERAJKJS8iIhJQKnkREZGAUsmLiIgElEpeREQkoFTyIiIiAeWt5M0sbGYzzOzp3P0tzexNM5ttZo+YWSw3XpS7Pyf3eH9fmUVEOrply5bx7LPPMnXqVBKJhO84sg4+l+TPBT5odv8q4Hrn3ABgGXB6bvx0YJlzbmvg+tx0IiLSzh577DH2229/zjvvfM4551yGDRvGf/7zH9+xZC28lLyZbQ4cDtyVu2/AcODR3CT3AUflbh+Zu0/u8QNy04uISDv56quvuPTSywiHo1RUVFJeXkEmAxMm/JRkMuk7nnwHX0vyNwC/ALK5+92A5c65dO7+fKBv7nZf4HOA3OMrctOLiEg7eeGFF2hoaCQU+qY2IpEI9fX1TJs2zWMyWZt2L3kz+xHwtXOu+adiTUvmbj0eaz7fs8ys2syqFy1a1ApJRURECpuPJfkhwAgz+xR4mKbV9DcAnc0skptmc+DL3O35QD+A3OOdgKXfnqlzbpJzrso5V9WjR4+2/QlERDqYAw88kOLiIrLZ7KqxdDpNSUkJu+++u8dksjbtXvLOuV855zZ3zvUHRgP/dM6dCLwEHJub7FTgydztp3L3yT3+T+fcakvyIiLSdnr16sWll/6GbDZFbW0N8Xgd4bBx6623EIvFfMeT7xBZ9yTt5kLgYTO7HJgB3J0bvxv4k5nNoWkJfrSnfCIiHdrIkSPZf//9ef3114lGowwZMoTS0lLfsWQtLIgLxVVVVa66utp3DBERkQ1iZtOcc1WtNb98WpIXEVkviUSCadOmUVxczKBBgwiHw74jieQllbyIFJSnn36aSy65hMbGJKFQmE6dKrjnnnvYZpttfEcTyTs6d72IFIyFCxfy619fBIQoKSmjqKiYmpo6zjjjjBZ7fYtIE5W8iBSMp59+OndClm9Wz8diRSxfvoL33nvPYzKR/KSSF5GC0djY+B2PGKlUql2ziBQClbyIFIxDDz2UWCxK86OCMpk0xcVF7Lzzzh6TieQnlbyIFIytttqKs8/+MclkAytWLKe2tgZwXHfdH4hGo77jieQdHScvIgVn3rx5vPTSS5SUlHDwwQfTqVMn35FEWoWOkxeRDm+LLbbg1FNPXfeEIh2cVteLiIgElEpeREQkoFTyIiIiAaWSFxERCSiVvKzT0qVLueKKKzjggAM44YQTeP31131HEhGR9aCSl7Wqqanh6KOP5r77/sSyZTW8++4sTj/9DKZMmeI7moiIrINKXtbqoYceYuHCr6ms7EQoFKKoqIji4lKuvfZaMpmM73giIrIWKnlZq9dee42ioqIWY6FQiGQyxZIlSzylEhGR9aGSl7XadtttSSaTLcacc4RCprOMiYjkOZW8rNXYsWMpLS1edfWvbDZLIlHHcccdt9oSvoiI5BeVvKxVnz59ePDBBxkw4PtkMklisQg//ekEJk6c6DuaiIisg85dL+u0/fbb89hjf/UdQ0RENpCW5EVERAJKJS8iIhJQKnkREZGAUsmLiIgElEpeREQkoFTyIiIiAaWSFxERCSiVvIiISECp5EVERAJKJR8AtbW1q11ERkRERCVfwN59910OPvhg9tlnH/bcc09+/vOfU19f7zuWiIjkCZ27vkAtWrSIU08dSzqdJRaL4Zzj739/jrq6OiZNmuQ7noiI5AEtyReoRx55hESinlgsBoCZUVZWzltvvcXixYs9pxMRkXygki9Qn3/+OaHQ6v98zsHSpUs9JJJ8kkql+Oijj/QHn0gHp5IvUPvtt99qY9lslkgkTP/+/ds9j+SPJ598kr333pvjjjueAw/8IaeffjrxeNx3LBHxQCVfoA488EB22GE76uqa9qxPJBI0Nib4xS9+sWoVvnQ8H3zwARdddDGZDEQiMSKRGK+99ia/+MUvfEcTEQ+0412Bikaj/PnPf+bRRx/lb3/7G926deOMM85g11139R1NPJo8eTKZTLbFppyysjJeffU14vE4ZWVlHtOJSHtTyRewWCzGCSecwAknnOA7iuSJxYsXE4ms/mttBolEQiUv0sFodb1IgBx++OGk0+kWY6lUis6dO9O9e3dPqUTEF5W8SICMGDGCXXYZSDxeRyKRoLa2hlDIcc0112BmvuOJSDvT6nqRAIlGozz44IO8+OKLTJ06lc0335zRo0fTq1cv39FExANzzvnO0OqqqqpcdXW17xgiIiIbxMymOeeqWmt+Wl0vspGccyxatIja2lrfUURE1kir60U2wsyZMznvvPNYunQpzsHuuw/i+uuvp7Ky0nc0EZFVtCQvsoGWLl3KuHGnsXjxUsLhlSeceYPx48f7jiYi0oJKXmQDPfHEE9TVxYlGvzmzYFlZObNmfcD8+fM9JhMRaUklL7KB5s+f/x2Ho5kuDiQieUUlL7KBhg8fTjjc8len6eJAIQYMGOAplYjI6lTyIhto7733ZvDgKuLxOhoaGojH4zQ21vOLX/yCkpIS3/FERFbR3vUiGygUCnH33Xfzj3/8gyeeeIJOnToxduxYdtppJ9/RRERa0MlwRERE8oROhiMiIiLrRSUvIiISUCp5ERGRgFLJi4iIBJRKXkREJKBU8iIiIgGlkhcREQkolbyIiEhAqeRFREQCSiUvIiIF5/333+fUU09lyJChnHbaafz3v//1HSkvqeRFRKSgvP3224wePYZp02bQ2JjizTerOe6445k1a5bvaHlHJS8iIgXliiuuAEIUFzdd9bGkpIRs1vG73/3Ob7A8pJIXEZGCMm/ePIqKilqMFReX8Mknn3hKlL90qVmR9ZBKpXjhhRd455132HHHHTn44IOJxWK+Y4l0SL169eKzzz4nGv3mdzCZbKR//609pspPKnmRdaipqWHUqOOYP/8LMpkMoZBx/fXXM2XKFLp16+Y7nkiHM3HiRH7847MBIxqNkkqlyGTSTJw40Xe0vKPV9SLrcMMNN/DZZ/MoKSmlvLyC0tJyvvxyIVdeeaXvaCId0rBhw7j11lvo02czzLJsvnlv7rzzDvbcc0/f0fKOOed8Z2h1VVVVrrq62ncMCYghQ4bQ0JDCzFaNOeeIRIw333zTYzIRCRozm+acq2qt+WlJXgTIZrPMmTOH+fPnr/ZYcXEx2Wy2xZhzjmg02l7xREQ2ikpeOry33nqLYcOGceyxo/jRj47giCNGsGDBglWPn3zyydTXJ2i+1iuRiHP88cf7iCsist5U8tKhLV26lLPOGk88Xk8kEiMSifHJJ3MZN+60VaV+yimnMHLkUWQySZLJRlKpRg499GDOPvtsz+lFRNZOe9dLh/bMM8+QSNRTUVG5aqykpJQFCxbw4Ycfst122xEKhfjd737Heeedx2effUa/fv3o1auXx9QiIutHJS8d2pIlS1jTzqfZrKOmpqbFWI8ePejRo0d7RRMR2WRaXS8d2oEHHkgsFm1R9Nlslmg0wsCBAz0mExHZdCp56dB23HFHRow4goaGBPF4HbW1NaRSjfzv/15MSUmJ73giIptEx8lLh+eco7q6mqeeeory8nJGjRrFVltt5TuWiHRArX2cvLbJS4dnZgwePJjBgwf7jiIi0qrafXW9mfUzs5fM7AMze9/Mzs2NdzWzqWY2O/e9S27czOwmM5tjZjPNbFB7Z5b8sWTJEu68804uueQS/vWvf612khoREfmGjyX5NHC+c266mVUA08xsKjAWeNE593sz+yXwS+BC4FBgQO5rT+D23HfpYGbOnMm4ceOoq0sQCoX4618fZ7fdduHee+8lEtFKKRGRb2v3JXnn3ALn3PTc7VrgA6AvcCRwX26y+4CjcrePBO53Td4AOptZ73aOLZ455/j5z88nnc5SUVFJWVk5JSWlVFdP5+mnn/YdT0QkL3ndu97M+gO7AW8CvZxzC6DpDwGgZ26yvsDnzZ42Pzf27XmdZWbVZla9aNGitowtHixbtozFixe3uH40QDQa48knn/SUSkQkv3kreTMrB/4K/I9zrmZtk65hbLVDApxzk5xzVc65Kp2wJHiKiopYwz87mUyaLl26tH8gEZEC4KXkzSxKU8E/6Jx7LDf81crV8LnvX+fG5wP9mj19c+DL9soq+aGsrIy99tqLeDy+aiybzWIGY8eO9RdMRCSP+di73oC7gQ+cc9c1e+gp4NTc7VOBJ5uNn5Lby34vYMXK1frSsVxzzTVUVe1GOt1IKpXELMsll/wvO++8s+9oIiJ5qd1PhmNmQ4H/B7wLrDz+6dc0bZf/C7AFMA8Y5Zxbmvuj4BbgECABjHPOrfVMNzoZTrAtWLCAJUuWsPXWW1NcXOw7johIqyn4k+E45/7NmrezAxywhukdMKFNQ0lB6d27N7176wALEZF10bnrRUREAkolLyIiElAqeRERkYBSyYuIiASUSl5ERCSgVPIiIiIBpZIXEREJKJW8iIhIQKnkRUREAkolLyIiElAqeRERkYBSyYuIiARUu1+gRkTa3ueff86MGTPo1q0be+21F+Fw2HckkcAopN8vlbxIgDjnuPzyy5ky5VESiXqKioro0aMbDz74IH379vUdT6SgFeLvl1bXiwTIv//9bx566GGi0SI6d+5CSUkpixYt5ZxzzvEdTaTgFeLvl0peJED+/Oc/YxbGzFaNlZSU8Mknc1m2bJnHZCKFrxB/v1TyIiIiAaWSFwmQE044AecyOOdWjdXX17PVVlvSpUsXj8lECl8h/n6p5EUCZOjQoYwZM5pUqpHly5fR0JCgR4+u3Hzzzb6jiRS8Qvz9suZ/kQRFVVWVq66u9h1DxJtCOsRHpNC05e+XmU1zzlW11vx0CJ1IAPXr149+/fr5jiESSIX0+6XV9SIiIgGlkhcREQkolbyIiEhAqeRFREQCSiUvIiISUCp5ERGRgFLJi4iIBJRKXkREJKBU8iIiIgGlkhcREQkolbyIiEhAqeRFREQCSiUvIiISUCp5ERGRgFLJi4iIBJRKXkREJKBU8iIiIgGlkhcREQkolbyIiEhAqeRFREQCSiUvIiISUCp5ERGRgIr4DtBRTJ8+nUmTJrFkyVIOP/wwRo8eTXFxse9YIiISYCr5dvDwww/z299eDoSIRCLMnHkVjz/+OFOmTCEWi/mOJyIiAaXV9W0smUxy7bXXUlRUQklJCdFolIqKSj76aA7PPfec73giIhJgKvk29tlnn5FOZwiFWr7VZiFefvllP6FERKRDUMm3sW7duq1xPJvNsOWWW7ZzGhER6UhU8m2sa9euDBmyN/F4Hc45oGkVfllZCccff7zndCIiEmQq+Xbwhz/8gaOOGkE2myKdTrL55r25//776dGjh+9oIiISYLZy6TJIqqqqXHV1te8Yq0mlUqRSKUpLS31HERGRPGRm05xzVa01Px1C146i0SjRaNR3DBER6SC0ul5ERCSgVPIiIiIBpZIXEREJKJW8R/F4nI8++oja2lrfUcSjL7/8ko8//phsNus7iogEjHa888A5x7XXXsuf//zQqmPnR448mosvvni1M+NJcH311VeMHz+euXM/xcwoKyvlhhtuYPDgwb6jiUhAqFE8mDJlCvfcM5lQKEIkEiMcjvLQQw8zefJk39GknTjnOO200/joo4+JRouIRGLE4/WMHz+epUuX+o4nIgGhkvfgrrvuori4BDMDwMwoKSnl/vvv95xM2svs2bOZP/+LFudMiESixOP1PPXUUx6TiUiQqOQ9qKuLEwqFW4yZhaivr/eUSNpbbW0t2eyaT0SlJXkRaS0qeQ/2229f4vG6FmOJRJw99tjDUyJpbzvssAOxWKTFznbOOaLRCAceeKDHZCISJCp5D84//3w226wn8Xgd9fUJ6upq6datCxdddJHvaNJOSkpKuPjii0mlGqmtrSUer6OhIcFhhx3KwIEDfccTkYDQues9SSQSPPHEE8yYMYOddtqJkSNHUlFR4TuWtLNPPvmERx55hJqaGo466ij22GOPVftqiEjH09rnrlfJi4iI5InWLnmtrhcREQkolbyIiEhAqeRFREQCSiUvIiISUCp5ERGRgFLJi4hIh7R8+XJ+//vfs//+wxk1ahQvv/yy70itTiUvIiIdTiKR4JhjjmHy5PtYsaKWDz+cw9ln/4R7773Xd7RWpZIXEZEO57HHHuOLLxZQUVFJKBQiFotRWlrGLbfcSjKZ9B2v1ajkRUSkw3n99deJRCItxsxCZDIZFi5c6ClV61PJi4hIh7PNNtuQTqdbjK08A2y3bt18RGoTKnkREelwTjzxRCoqymhoaLrEdzabJZGIM2LEEZSVlXlO13pU8iIi0uF0796dhx56iJ122oFMJkU0GuLss8dzySWX+I7WqiLrnkRERCR4tt56ax5++GHfMdpUwSzJm9khZvahmc0xs1/6ziMiIpLvCqLkzSwM3AocCuwAjDGzHfymEhERyW8FUfLAHsAc59wnzrkk8DBwpOdMIiIiea1QSr4v8Hmz+/NzYyIiIvIdCmXHO1vDmGsxgdlZwFm5u41m9l6bp+qYugOLfYcIIL2vbUfvbdvQ+9o2tm3NmRVKyc8H+jW7vznwZfMJnHOTgEkAZlbtnKtqv3gdh97btqH3te3ovW0bel/bhplVt+b8CmV1/VvAADPb0sxiwGjgKc+ZRERE8lpBLMk759Jm9lPgOSAM3OOce99zLBERkbxWECUP4Jx7Fnh2PSef1JZZOji9t21D72vb0XvbNvS+to1WfV9t5Qn5RUREJFgKZZu8iIiIbKDAlbxOf7vxzKyfmb1kZh+Y2ftmdm5uvKuZTTWz2bnvXXLjZmY35d7rmWY2yO9PkN/MLGxmM8zs6dz9Lc3szdz7+khup1LMrCh3f07u8f4+c+c7M+tsZo+a2X9zn90f6DO76czsvNz/A++Z2UNmVqzP7MYxs3vM7Ovmh3ZvzGfUzE7NTT/bzE5dn9cOVMnr9LebLA2c75zbHtgLmJB7/34JvOicGwC8mLsPTe/zgNzXWcDt7R+5oJwLfNDs/lXA9bn3dRlwem78dGCZc25r4PrcdPLdbgT+4ZzbDtiFpvdYn9lNYGZ9gZ8BVc65nWja4Xk0+sxurHuBQ741tkGfUTPrCvwG2JOms8D+ZuUfBmsTqJJHp7/dJM65Bc656bnbtTT9Z9mXpvfwvtxk9wFH5W4fCdzvmrwBdDaz3u0cuyCY2ebA4cBdufsGDAcezU3y7fd15fv9KHBAbnr5FjOrBPYB7gZwziWdc8vRZ7Y1RIASM4sApcAC9JndKM65V4Cl3xre0M/owcBU59xS5/5/e3cTGlcVhnH8/2KltRW0KgoiWKMUF4qNixKsYvGjiyK66S5UqaIbFaNgQcSFKwXr18YiVFyICH4UFTe2anHhorVFMaJSUxJiJZou2ggqtZLHxTnT3kwmyZ2PTOqd5weXmXvumcmdw0veOefeOUfHgb3M/uIwS9WSvKe/7ZA83NYP7AcukzQB6YsAcGmu5vYu7xVgOzCd9y8GTkj6N+8X2+50u+bjU7m+zdYHHAPezJdCdkXEKhyzbZH0K7ADGCcl9yngEI7ZTmo2RluK3aol+QWnv7WFRcT5wAfAkKQ/5qvaoMztXSci7gImJR0qFjeoqhLHbKZlwI3ATkn9wJ+cGfZsxG1bQh4Gvge4CrgcWEUaRq7nmO28udqypTauWpJfcPpbm19EnEtK8G9L2p2Lf68NaebHyVzu9i5nA3B3RIyRLiHdRurZX5iHQmFm251u13z8AmYP9VlyFDgqaX/ef5+U9B2z7bkDGJV0TNIpYDdwE47ZTmo2RluK3aoleU9/24Z8De0N4EdJLxUOfQzU7uS8D/ioUH5vvht0AJiqDT/ZGZKeknSFpDWkmPxC0iCwD9iSq9W3a629t+T67hU1IOk34JeIqC3qcTvwA47Zdo0DAxGxMv9fqLWrY7Zzmo3RT4FNEbE6j7RsymXzk1SpDdgMHAaOAE8v9fn8nzbgZtLwz3fAt3nbTLq29jnwc368KNcP0q8ZjgDDpDtxl/xznM0bsBH4JD/vAw4AI8B7wPJcviLvj+TjfUt93mfzBqwDDsiLTSgAAAIESURBVOa4/RBY7ZjtSLs+C/wEfA+8BSx3zLbclu+Q7m04ReqRP9BKjAL35zYeAbaV+due8c7MzKyiqjZcb2ZmZpmTvJmZWUU5yZuZmVWUk7yZmVlFOcmbmZlVlJO8mZlZRTnJm/WwSMsLj+YVrsgTbYxGxJUN6p4XEV/m1R7Lvv8jEbGtk+dsZuX5d/JmPS4itgPXSHooIl4HxiQ916Dew8AySa828d4rga+U5pU3sy5zT97MXiZNYTpEmvXwxTnqDZKn3oyIjblX/25EHI6I5yNiMCIORMRwRFwNIOkvYCwi1nfjg5jZTE7yZj1OaQGSJ0nJfkjSP/V18loQfZLGCsU3AI8B1wNbgbWS1gO7gEcL9Q4CtyzO2ZvZfJzkzQzSMqITwHVzHL8EOFFX9rWkCUknSfNs78nlw8CaQr1J0nKlZtZlTvJmPS4i1gF3AgPA47XlL+v8TVqEpOhk4fl0YX+atM57zYr8ejPrMid5sx6WlxHdSRqmHwdeAHbU15N0HDgnIuoTfRlrSSuZmVmXOcmb9bYHgXFJe/P+a8C1EXFrg7p7SDfmNWsD8FmL52dmbfBP6MyslIjoB56QtHUxX2NmneOevJmVIukbYF8zk+GQbth7ZpFOycwW4J68mZlZRbknb2ZmVlFO8mZmZhXlJG9mZlZRTvJmZmYV5SRvZmZWUf8Bw2Fiyhcw8TwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.subplot(111)   \n",
    "im = plt.scatter(df['X'],df['Y'],s=None,c=df['null'],marker=None,cmap=cmap,norm=None,vmin=0,vmax=1.0,alpha=0.8,\n",
    "        linewidths=0.8,verts=None,edgecolors=\"black\",)\n",
    "plt.title('Input Data Locations')\n",
    "plt.xlim(0, 1000); plt.ylim(0, 1000)\n",
    "plt.xlabel('X (m)'); plt.ylabel('Y (m)')\n",
    "\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.5, wspace=0.3, hspace=0.3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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": 19,
   "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",
       "      <th>AI</th>\n",
       "      <th>null</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>72</th>\n",
       "      <td>250</td>\n",
       "      <td>50</td>\n",
       "      <td>0</td>\n",
       "      <td>0.139637</td>\n",
       "      <td>0.347182</td>\n",
       "      <td>4747.274043</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>153</th>\n",
       "      <td>650</td>\n",
       "      <td>750</td>\n",
       "      <td>0</td>\n",
       "      <td>0.170732</td>\n",
       "      <td>10.720560</td>\n",
       "      <td>4535.625583</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>258</th>\n",
       "      <td>80</td>\n",
       "      <td>669</td>\n",
       "      <td>1</td>\n",
       "      <td>0.244345</td>\n",
       "      <td>3222.716042</td>\n",
       "      <td>2696.102930</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>56</th>\n",
       "      <td>200</td>\n",
       "      <td>150</td>\n",
       "      <td>0</td>\n",
       "      <td>0.167125</td>\n",
       "      <td>3.042590</td>\n",
       "      <td>5500.997419</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>303</th>\n",
       "      <td>60</td>\n",
       "      <td>929</td>\n",
       "      <td>1</td>\n",
       "      <td>0.216253</td>\n",
       "      <td>400.298484</td>\n",
       "      <td>3959.934912</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>374</th>\n",
       "      <td>750</td>\n",
       "      <td>699</td>\n",
       "      <td>0</td>\n",
       "      <td>0.159171</td>\n",
       "      <td>9.224292</td>\n",
       "      <td>5263.064063</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>333</th>\n",
       "      <td>340</td>\n",
       "      <td>549</td>\n",
       "      <td>0</td>\n",
       "      <td>0.170881</td>\n",
       "      <td>84.471492</td>\n",
       "      <td>2918.232227</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>50</td>\n",
       "      <td>850</td>\n",
       "      <td>1</td>\n",
       "      <td>0.237154</td>\n",
       "      <td>39.837129</td>\n",
       "      <td>3074.562617</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>54</th>\n",
       "      <td>200</td>\n",
       "      <td>250</td>\n",
       "      <td>1</td>\n",
       "      <td>0.188043</td>\n",
       "      <td>57.804784</td>\n",
       "      <td>4997.078597</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>311</th>\n",
       "      <td>250</td>\n",
       "      <td>579</td>\n",
       "      <td>1</td>\n",
       "      <td>0.215039</td>\n",
       "      <td>625.510541</td>\n",
       "      <td>2693.691341</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>50</td>\n",
       "      <td>450</td>\n",
       "      <td>1</td>\n",
       "      <td>0.200389</td>\n",
       "      <td>265.636019</td>\n",
       "      <td>3454.389302</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>269</th>\n",
       "      <td>540</td>\n",
       "      <td>859</td>\n",
       "      <td>0</td>\n",
       "      <td>0.173662</td>\n",
       "      <td>39.910675</td>\n",
       "      <td>3414.962471</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>279</th>\n",
       "      <td>580</td>\n",
       "      <td>779</td>\n",
       "      <td>0</td>\n",
       "      <td>0.168371</td>\n",
       "      <td>11.171392</td>\n",
       "      <td>3646.918994</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       X    Y  Facies  Porosity         Perm           AI  null\n",
       "72   250   50       0  0.139637     0.347182  4747.274043   0.0\n",
       "153  650  750       0  0.170732    10.720560  4535.625583   0.0\n",
       "258   80  669       1  0.244345  3222.716042  2696.102930   0.0\n",
       "56   200  150       0  0.167125     3.042590  5500.997419   0.0\n",
       "303   60  929       1  0.216253   400.298484  3959.934912   0.0\n",
       "374  750  699       0  0.159171     9.224292  5263.064063   0.0\n",
       "333  340  549       0  0.170881    84.471492  2918.232227   0.0\n",
       "1     50  850       1  0.237154    39.837129  3074.562617   0.0\n",
       "54   200  250       1  0.188043    57.804784  4997.078597   0.0\n",
       "311  250  579       1  0.215039   625.510541  2693.691341   0.0\n",
       "8     50  450       1  0.200389   265.636019  3454.389302   0.0\n",
       "269  540  859       0  0.173662    39.910675  3414.962471   0.0\n",
       "279  580  779       0  0.168371    11.171392  3646.918994   0.0"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head(n = 13)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Number of Effective Data\n",
    "\n",
    "Let's try out our program for a variety of spatial continuity models.\n",
    "\n",
    "* **100% nugget** - should have number of effective data equal to the number of data as they are all independent.\n",
    "\n",
    "* **Zonal Anisotropy** - should have number of effective data equal to 1.0 as they are all perfectly redundnat.\n",
    "\n",
    "* **Range of Fraction of Model Size** - should have a low number of effective data then increasing as the fraction decreases.\n",
    "\n",
    "* **Other** - testing the code.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100% Nugget:                  Number of effective data = 48 / 48.\n",
      "Zonal Anisotropy:             Number of effective data = 1.0 / 48.\n",
      "Range of Model Size:          Number of effective data = 2.14 / 48.\n",
      "Range of 1/2 Model Size:      Number of effective data = 5.69 / 48.\n",
      "Range of 1/4 Model Size:      Number of effective data = 15.89 / 48.\n",
      "Nested Variogram:             Number of effective data = 24.32 / 48.\n",
      "Nested Anisotropic Variogram: Number of effective data = 42.5 / 48.\n"
     ]
    }
   ],
   "source": [
    "vario = GSLIB.make_variogram(nug=1.0,nst=1,it1=1,cc1=0.0,azi1=0.0,hmaj1=100,hmin1=100)\n",
    "n_eff = n_effective(df,'X','Y',seed=73073,nreal=1000,vario=vario)\n",
    "print('100% Nugget:                  Number of effective data = ' + str(round(n_eff,2)) + ' / ' + str(len(df)) + '.')\n",
    "\n",
    "vario = GSLIB.make_variogram(nug=0.0,nst=1,it1=1,cc1=1.0,azi1=0.0,hmaj1=1000000,hmin1=1000000)\n",
    "n_eff = n_effective(df,'X','Y',seed=73073,nreal=1000,vario=vario)\n",
    "print('Zonal Anisotropy:             Number of effective data = ' + str(round(n_eff,2)) + ' / ' + str(len(df)) + '.')\n",
    "\n",
    "vario = GSLIB.make_variogram(nug=0.0,nst=1,it1=1,cc1=1.0,azi1=0.0,hmaj1=1000,hmin1=1000)\n",
    "n_eff = n_effective(df,'X','Y',seed=73073,nreal=1000,vario=vario)\n",
    "print('Range of Model Size:          Number of effective data = ' + str(round(n_eff,2)) + ' / ' + str(len(df)) + '.')\n",
    "\n",
    "vario = GSLIB.make_variogram(nug=0.0,nst=1,it1=1,cc1=1.0,azi1=0.0,hmaj1=500,hmin1=500)\n",
    "n_eff = n_effective(df,'X','Y',seed=73073,nreal=1000,vario=vario)\n",
    "print('Range of 1/2 Model Size:      Number of effective data = ' + str(round(n_eff,2)) + ' / ' + str(len(df)) + '.')\n",
    "\n",
    "vario = GSLIB.make_variogram(nug=0.0,nst=1,it1=1,cc1=1.0,azi1=0.0,hmaj1=250,hmin1=250)\n",
    "n_eff = n_effective(df,'X','Y',seed=73073,nreal=1000,vario=vario)\n",
    "print('Range of 1/4 Model Size:      Number of effective data = ' + str(round(n_eff,2)) + ' / ' + str(len(df)) + '.')\n",
    "\n",
    "vario = GSLIB.make_variogram(nug=0.0,nst=2,it1=1,cc1=0.5,azi1=0.0,hmaj1=50,hmin1=50,it2=1,cc2=0.5,azi2=0.0,hmaj2=250,hmin2=250)\n",
    "n_eff = n_effective(df,'X','Y',seed=73073,nreal=1000,vario=vario)\n",
    "print('Nested Variogram:             Number of effective data = ' + str(round(n_eff,2)) + ' / ' + str(len(df)) + '.')\n",
    "\n",
    "vario = GSLIB.make_variogram(nug=0.0,nst=2,it1=1,cc1=0.5,azi1=0.0,hmaj1=50,hmin1=5,it2=1,cc2=0.5,azi2=0.0,hmaj2=250,hmin2=25)\n",
    "n_eff = n_effective(df,'X','Y',seed=73073,nreal=1000,vario=vario)\n",
    "print('Nested Anisotropic Variogram: Number of effective data = ' + str(round(n_eff,2)) + ' / ' + str(len(df)) + '.')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6AAAAGWCAYAAABmee4oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XuYJGV58P/vDcuKIDAcVuW0rgckonhilU2MiKK+ggfwDSi+RgExG+M5mleIMa/ExATzSzzlgKIYWQ8gGg9oNEpQNFEWnVWEyGiWBWTXXWA4DIdFXNe9f3/UM2zvbHdvz0x3TU/393NdfU1XPdXV91NTXU/dVU9VRWYiSZIkSVKv7TTXAUiSJEmShoMJqCRJkiSpFiagkiRJkqRamIBKkiRJkmphAipJkiRJqoUJqCRJkiSpFiag0jwQET+JiKNr+J7FEXFPROzc6+/qhYh4ekT8bK7jkDR3IuLtEfHRDqf9eET8Va9jUn0i4uUR8Y2avutDEfHndXxXL0TE1yLilLmOoxsiIiPiUV2c3w0R8exuzW+a3z3jfb75sh9kAjpkyg/qlyXJuCMi/i0iDu7CfJeUH/+CDqfv6oaig+87NSL+q8n42jcw011WAJn52My8rMP5z7hOmXljZj4oM38z3c9GxNERsaWsW3dHxM8i4rSZxDFTmfmfmXloQ0yz/v9GxO6lTl+dfYTScIqI342I70XEnRFxe0R8NyKe0oX5Hh0R6xrHZeZfZ+aruzDvUyPiN+X3f09EXBcRfzTb+ZZ5d5z4RsRZEfHJbnxvp5ptO1u1ozXEMq2DBJn5qcx8bofznlWdMvM1mfmXM/lsRFwWEfeVdevWiPh8ROw/01hmIjOPzczzSzyzWhYRcVBE/Gupy50RcXVEnFrKpr3f08/KOrmp7OvcHRH/HRF/ExF7dWP+09zn22Z/eup+UL8yAR1OL8zMBwH7AzcD/zDH8WxnUDZSUw1qvRqsL+vWnsAZwEci4rDpzKAPl9GJwK+A5/Zq56AP6yx1TUTsCXyFqq3ZBzgQ+Auq31W/u7wclHsQ1bbgbyPiSXMdVKOoDOT+3HztjTMNry/r1qOBEeB9051BH7UfnwDWAg8D9gVeSbWP2ddmsfz+NjP3ABYBpwHLgO9GxO5dC26ADeQGS53JzPuAzwH3JwgRsVdErIiI8Yj4eUS8Y7Jhi4idyvDPI+KWMt3k0Z7vlL8T5Wjeb0fEoyLi2+VI2K0R8Zkyn8lpf1ymfenkUeyIOCMibgL+JSL2joivlFjuKO8Paoj1snLE6fvlO74UEfvMdHmUI1r/FNVZ4bsj4oqIeGRD+WMj4pJy9P7miHh7w3I5MyLWRMRtEXHRZBwNR/1Oj4gbgW+2WFaPjIhvls/fGhGfioiRhu++/2h0ORp+UVn+d0fVVWNpKfsEsBj4cpn320p93jClrldFxAlNlsE2RynLMv7LqM5W3B0R34iI/Xa0LLPyReAOyvoVES8qsU6U+T5mSv3OiIirgI0RsSAiHlOmmyife1HD9MdFxDUlpl9ExJ+U8fefDZntsmhwCvAh4Crg5Q2fOzMiPjdlXh+IiA+W93tFxHkRsaHE+FdRdqaiOtL83Yh4X0TcDpzVwTrw5Ij4UanzZyPiM9FwViAiXhARV5bl9b2IePyO/k9STR4NkJkXZOZvMvOXmfmNzLwKtvk9/EPZlv80Io6Z/HBEnBYRY2Xdvy4i/rCM3x34GnBAbD1LeUBMOWNYfi83lXl/JyIeO5NKZOYPgTGgcdvVbrvWdBsWEcuptiVvKzF/uYw/o2wrJnuQHBMRzwPeDry0TPvjMu1lEfHuiPgucC/wiFbLqUw/2ca+vWxfboiI+7dn09XQVpwSETeWef5ZQ/nO5bvWlHhWReltFRG/FVvb0p9FxEsaPvfxiDgnIr4aERuB01ssqzMb5n1NRLy4YR7bnMkrcb4mIlZHtS/xT1F5DNW2/bfLvCci4ilRte8LGj7/exFxZYvlcP/Z2YZl/Nao9pE2RIe9gDLzduBfgceVebXbF2vWfrTcP4uIXSPik1G1LRMR8YOIeEgpuywiXt2NZQE8Bfh4Zm7MzM2Z+aPM/Fopm+l+z59E1UbfGVWbt2tD+f8ty3h9RLxqyv/l+VG1l3dFxNqIOKuhrNl+GRHxirL8bmtclzv4392XmT8AXkSVeN//P4+IV0X1m7wjIr4eEQ8r4z8UEX83JeYvRcRbGuo+uc/31Ii4vPxPNkTEP0bEwlLWcn+6Yb7t9qXa7vf2VGb6GqIXcAPw7PJ+N+B8YEVD+QrgS8AewBLgf4DTS9mrgGuBRwAPAj4PfKKULQESWNAwrwuAP6M60LEr8LsNZQk8qmH4aGAz8B7gAcADqX7Iv1fi3AP4LPDFhs9cBvyCaoO9O9XG+5Mt6n0q8F87WB4fB24HngosAD4FXFjK9gA2AG8tddkDOLKUvRlYCRxUYv8wcMGU5bKixPjAFsvqUcBzyucXUW2s398izrOA+4DjgJ2BvwFWNpu2DL8EuKJh+AnAbcDCJstjm9jKMl5DtQP5wDJ8dotlfDSwrrzfCXgx8Gvg0PL5jaWOuwBvo1qXFjbEfCVwcPmeXUr524GFwLOAu4FDy/QbgKeX93sDT54aw2yXRSlfDGyhSqLfClzVUPYwqh2/PcvwziWuZWX4i2Vd2B14MPB94A8b1sfNwBuo1rUHtlsHyjL4OfCmsmz+N7AJ+KtS/mTgFuDIEscppe4PmOttji9fVD0ibqNqb44F9p5SPvl7+OOyfr8UuBPYp5Q/H3gkEMAzyu+u6W++jDuLhraAqu3ao/y23g9c2VD28cnfUZO4T6Wh3aDawZ4AHl2GW27X2PE2bJvvpdpOrgUOKMNLgEc2q08ZdxlwI/DYsg3ZpYPltBl4b1kOzyixH9qi7jfQsO2cujzY2lZ8hGr79QSqM9qPKeX/F7i61CtK+b5U28O1VDvpC6i2XbcCj21YLncCT2PrvsN2/yPgJOCAMs1LS132b/F/S6oz8CNU2/Rx4HnNpi3jrgGObRj+AvDWFsvp/tgalvG7yv/juPI/2LvFZy8DXl3e70eVCE3uU7XbFzuV7duPdvtnfwh8mWpfamfgCLa2W40xzHZZ/AfwXeBkYPGUsiXMbL/n++X/vA/VwZ/XlLLnUZ1dndz/+zQN+5Xlf3F4WT8eX6Y9YUosjftlhwH3AEeVeN5blvGzW9T1/v/7lPErgM+U9yeU/8ljyv/pHcD3StlRVL+DKMN7A79k6+//Brbu8x1BdXZ1QYl9DHjzlPV76v705L5YJ9uhpvu9vX7NecPkq95XWanvoWpENwPrgcNL2c5UDchhDdP/IXBZeX8p8NqGskOpEozJH8XUjcsK4FzgoCZxNPvBbAJ2bRP7E4E7GoYvoyEZKhuQTcDOTT57Kp0loB9tKDsO+Gl5/zLgRy3iGgOOaRjev8lyeURD+XbLqsk8T2j8PrZPQP9jSr1/2WzaMvyAsoE5pAz/HfDPLb53m9jKMn5HQ/lrgX9v8dmjqZK1ifJ9VwInl7I/By5qmHYnqoMHRzfE/KqG8qcDNwE7NYy7ADirvL+Rat3cs0kM7RLQjpdFKX8HZWeVqhH8DfCkhvL/Al5Z3j8HWFPeP4Tqt/TAhmlfBnyrYX28cQe/1fvXAarG6heUxqrhuyd3fM4B/nLK538GPKPdd/jyVdeLaifs48A6qrbnYuAhpexUqraocf3+PvCKFvP6IvCm8n6b33wZdxatD0aOlG3cXmX447RPQDeXbdo95XP/wNadxpbbtQ62Ydt8L9XO+C3As4FddlQfqm3zu3awzKcup83A7g3lFwF/3uKzN9BZAnpQQ/n32brN/xlwfJP5vhT4zynjPgy8s2G5rJhS3vJ/1DDNlZPfR/MEtPEA+EXAmc2mLePOAD5V3u9DlUTu3+J774+tLONfsu1+0C2Ug5JNPntZmfdEWW8+RZWI7Whf7FSmtB+03z97FfA94PEtYmiXgE5nWewNnA38hKqtvBJ4ypT1Zbr7Pb/fMPy3wIfK+4+x7f7fo5myXzll3u8H3jcllsb9sv9HQ+JFlZhuYvoJ6NnAJeX91ygHDcrwTmX5PYzqoMyNwFGl7A+Ab7b7/TWUvRn4wpT1u1UC2sl2qOl+b69fdsEdTidk5gjVzvjrgW9HxEOpjsBNnmmZ9HOq63Wg2gGfWraAame7mbdR/ci+X077v6rFdJPGs+oWDEBE7BYRHy5dIu6iOjo2EtteE7J2Sjy7lHpMtbmUTbUL1UZ60k0N7++lOpII1Zm5NS3ifhjwhdK9YYIqIf0N2y6XtU0/WUTEgyPiwqi6X90FfLJFPVrFuWu0uI4hM39F1eD+funC8zKqazU61WqZNLM+M0cyc5/MfGJmXljGb7PuZOYWqmVyYMNnG5fRAcDaMt2kxnXx96g2lD+Pqpv3b3dSkRksi1dS7RSQmeuBb1OdXZz06TIPgP9ThqFaJ3YBNjSsFx+mOhM6aZt1YgfrwAHAL7K0EE0+/zDgrZPfVb7v4PI5ac5l5lhmnpqZB1GdtTiAaqdw0tT1++dlGiLi2IhYGVWXzQmq3/4OLwUon905Is6OqrvmXVQ7dnT6eareJSNZXaf3UKozjn9dytpt13a0DdtGZl5LtWN5FnBL2Rbs6Pc7dRuyo+V0R2ZunBJPq+9o1mZObS9h+m3mw4Ajp2yrXk61bJvWq5mIeGVsveRggmqdmk6b2a4d+yTwwoh4EFWvmf/MzA07iqm4LTM3T+O73ljWrwMz8+WZOc6O98Vg+2XUbv/sE8DXgQtLd9W/jYhm+0PNdLwsMvOOzDwzMx9bvvdK4IsREc2m73C/p9X/7QC23/9rnPeREfGtqLow3wm8psm8t9vnaKjLRqpeG9N1INVBbqjW9Q80rKO3U+0TH1i2dRey7f7Dp5rNMCIeHdUlaDeV5fTXTerSSifboen8NrrGBHSIZXUtzuepkqXfpeoG82uqH82kxVRH5qA6Qj21bDNV14bGHYfJ+d+UmX+QmQdQHb3752h/59up83gr1VG8IzNzT6qzQFD9gCc13sF3cYn/1ibzvhFY3LghjIjdqBKCnzeZfqq1VF2bWpUdWxqRydeumfmLhmmyxftJf1PGP77U9ffZtp7T0Wz+51M18scA92bm5TOc90xts+6U/8PBbF23YNu41wMHx7Y31rh/XczMH2Tm8VT/vy9SJZXNzHhZRMTvAIcAf1o2/DdRdXF9WUOy/1ng6KiuTX4xWxPQtVRHsPdrWCf2LA1zq9jarQMbgAOnNOSN6/5a4N1T1sHdMvOC5otFmjuZ+VOqI++Paxg9df1eDKyPiAdQXV7xd1RnTEeAr7L1t9HsN97o/wDHU51Z3Ivq7AfMYPuamTeXWF5YRrXbrrXdhjWLOzM/nZm/W+aZVJekNJ126vgOlhPA3rHtDVIWlzibuZGty2rSw+msvYTWbeZa4NtTtlUPyszGuwtPre82w+U6uo9QHUDft9T1v5lZm9ns//AL4HKqbformN4B227Y0b4YbB93y/2zzPx1Zv5FZh4G/A7wAqqDq1N1bVlk5q1U6+Jk99lu7/dsYPv9v0afpuplcXBm7kV1fevUeTfGtM38yv7hvh3GMvmZB1FtZ/6zjFpLddlN47r+wMz8Xim/ADixrM9HUv1+mzkH+ClVz609qbrTdrqcdrQdmjMmoEMsKsdTdZsYy+rRGxcB746IPcqP4i1UR6Wg+rH8cUQ8vPzQ/pqqr/tmqmsqtlBdfzA5/5Ni602D7qD6sU8+3uPmxmlb2IOqO8tEVDf1eWeTaX4/Ig4rG4t3AZ/L5o8QuYLquskzo7ogf3eqrhKjdNagfgV4aES8OSIeUJbPkaXsQ1TL7GGl3ovKcm1lu2VV6npPqeuBVNfPzNR2y7YkWVuAv6f+xhSq9er5Ud1UYxeqgwu/ouoW1MwVVNf0vC0idonqeVgvpDqCuzCq57ztlZm/Bu5i63o11WyWxSnAJVRdnJ9YXo+juo7m2DKvcaouTP8CXJ+ZY2X8BuAbwN9HxJ5R3SDikRHxjDbf124duLzU8fVR3aDpeKprNiZ9BHhNOeobUT065vkRsUeb75NqEdVNZ9462R5EdTOal1FdOz/pwcAby+/9JKouu1+lOhP0AKrt5uaIOBZofMTGzcC+0frxB3tQbWtuo/rt/nWL6Tqpx75UO+I/KaPabddabsMa4m5sLw+NiGeVRPI+qravsb1cEu3vdLuj5TTpL8o29OlUichnW8zvM8Cby/8uorrR3asa4t+RjwJ/GRGHlM8/viy/rwCPjuqGL7uU11Oi4eZNTUzdju9OtT8xDtVNqtj2YMZ03AwcFOWmLg1WUPXiOpzqusfadLAv1kzL/bOIeGZEHB5V77G7qJLbZm3mrJZFRLwnIh5X2qg9gD8Crs3M2+j+fs9FwKkN+39T9w/3AG7PzPsi4qlUB6La+RzwgqgeF7WQan+yoxyp7BMeQXUw/A6q/QGo9g3/NMpNz6K6sdRJk5/LzB9RLZePAl/PzIkWX7EH1f/tnoj4Larl2qjd/vSOtkNzxgR0OH05Iu6hWqHfDZySmZMN6huoVtbrqK4x+zRVX3vK309QdYW9nqqRfANAZt5b5vXdqLobLKO6YcMV5bsuproW5foyr7OA88u0998Bb4r3U10cfivVjsq/N5nmE1RH0m+iulnBG5vNqHS9fD6lb3yp3wHAS6Z0+2oqM++musbvheW7VgPPLMUfKPX7RkTcXWI9stl8yryaLau/oLoZw53Av1HdQGCm/gZ4R5n3nzSMX0HVgNT6PDmAzPwZ1dHNf6D6f76Q6nFAm1pMv4nqjnLHlun/mepay5+WSV4B3BBVd5TXlHk3M6NlEdWd9l4C/EM5kz/5up5qnZvaDffZbD37OemVVDuF11A1Sp+juj64lZbrQFke/5vqjpATpb5foTzGIjNHqa4h+cfyXddSXc8j9YO7qbaJV0R1Z9OVVGes3towzRVUPQ5updo+npiZt5Vt7xupdjjvoNqRvHjyQ2WbcAFwXfmdT+1SuoLqIOMvqH6LK5meybuC3kN1ecU4W9u9ltu1DrZh5wGHlZi/SJU8nl2mvYkqIX97mXYySbwtIn7YLMgdLafiplK2nqq732sa4pnqI1Q70l+m2iatAP4sM5u1w828t8TyDap9jfOorom/myoxPrnEcRNbbz7YyjbLKjOvoTqAeDnVzvfhVDe/mYlvUh1QuCkiGntPfYFyec2Ubst1abcv1kzL/TOq7s2fo/o/jFFdStKs7ZvtstitTDtR4n4Y1W+g6/s9Wd1d9/0l5mvL30avBd5V9sn+H617SU3O7yfA66iW8waq38m6dp+hSurupupauwJYBfzO5DLKzC9QrdsXln2V/6YcvG5wAc33Hxr9CdXv+W6q3+VnppSfRYv96Q62Q3Nm8kJ6ad6JiMuobszw0bmOZT6IiFcCy0sXr6E2CMsiIq6guiHDv+xwYqmPRfWw+lfP599jvytnPj6Z1TW46kBErKHqQvkfcx3LXHNZqNs8AyoNgdJF5bVUdyUeavN1WUTEMyLioaV70ylUt5bv9GyEJKlDEfF7VN18p55ZGzouC/WCCag04CLif1F1G7uZ9t08Bt48XxaHAj+m6q70Vqouip3emVGS1IHSu+oc4HVT7h46dFwW6hW74EqSJEmSauEZUEmSJElSLZo+uH6+22+//XLJkiVzHYYkqQ+tWrXq1sxcNNdxzEe2r5KkZqbTtg5kArpkyRJGR0fnOgxJUh+KiE6e/asmbF8lSc1Mp221C64kSZIkqRYmoJIkSZKkWpiASpIkSZJqYQIqSZIkSaqFCagkSZIkqRYmoJIkSZKkWpiASpIkSZJqYQIqSZIkSaqFCagkSZIkqRYmoJIkSZKkWpiASpIkSZJqYQIqSZIkSaqFCagkSZIkqRYmoJIkSZKkWpiASpIkSZJqsWCuA5DUmatXrmTTxETL8oUjIxy+bFmNEUmSNLzatcu2yVJrJqDSPLFpYoIjFi1qWb5qfLzGaCRJGm7t2mXbZKk1u+BKkiRJkmphAipJkiRJqoUJqCRJkiSpFiagkiRJkqRamIBKkiRJkmphAipJkiRJqoUJqCRJkiSpFiagkiRJkqRaLJjrACRJkqR+dPXKlWyamGhadv3YGEcsWlRzRNL8ZwIqSZIkNbFpYqJlkrl6dLTmaKTBYAIqDYF2R3AXjoxw+LJlNUckSZKkYWQCKg2INWNjLcuuHxvjxKOOalq2any8VyFJklSLdgdaof6Dre3aZA/8atiZgEoDYsvGjXYTkiQNpXZdZaH+g63t2mQP/GrYeRdcSZIkSVItTEAlSZIkSbUwAZUkSZIk1aJn14BGxMeAFwC3ZObjyrh9gM8AS4AbgJdk5h0REcAHgOOAe4FTM/OH5TOnAO8os/2rzDy/VzFLkiRp8HhTIKl/9PImRB8H/hFY0TDuTODSzDw7Is4sw2cAxwKHlNeRwDnAkSVhfSewFEhgVURcnJl39DBuSZIkDRBvCiT1j54loJn5nYhYMmX08cDR5f35wGVUCejxwIrMTGBlRIxExP5l2ksy83aAiLgEeB5wQa/iliRJkuZCvz1ORuqFuh/D8pDM3ACQmRsi4sFl/IHA2obp1pVxrcZvJyKWA8sBFi9e3OWwJUmSpN7qt8fJSL3QLzchiibjss347UdmnpuZSzNz6aI2P1xJkiRJ0tyoOwG9uXStpfy9pYxfBxzcMN1BwPo24yVJkiRJ80zdCejFwCnl/SnAlxrGvzIqy4A7S1fdrwPPjYi9I2Jv4LllnCRJkiRpnunlY1guoLqJ0H4RsY7qbrZnAxdFxOnAjcBJZfKvUj2C5Vqqx7CcBpCZt0fEXwI/KNO9a/KGRJIkSZKk+aWXd8F9WYuiY5pMm8DrWsznY8DHuhiaJEmSJGkO9MtNiCRJkiRJA84EVJIkSZJUCxNQSZIkSVItTEAlSZIkSbUwAZUkSZIk1cIEVJIkSZJUCxNQSZIGRETcEBFXR8SVETFaxu0TEZdExOryd+8yPiLigxFxbURcFRFPntvoJUnDwARUkqTB8szMfGJmLi3DZwKXZuYhwKVlGOBY4JDyWg6cU3ukkqShYwIqSdJgOx44v7w/HzihYfyKrKwERiJi/7kIUJI0PExAJUkaHAl8IyJWRcTyMu4hmbkBoPx9cBl/ILC24bPryrhtRMTyiBiNiNHx8fEehi5JGgYL5joASXNrzdhYy7KFIyMcvmxZjdFImqWnZeb6iHgwcElE/LTNtNFkXG43IvNc4FyApUuXblcuSdJ0mIBKQ27Lxo0csWhR07JVnu2Q5pXMXF/+3hIRXwCeCtwcEftn5obSxfaWMvk64OCGjx8ErK81YKkPtDsQe/3YWMs2UtLMmIBKkjQAImJ3YKfMvLu8fy7wLuBi4BTg7PL3S+UjFwOvj4gLgSOBOye76krDpN2B2NWjozVHIw0+E1BJkgbDQ4AvRARU7funM/PfI+IHwEURcTpwI3BSmf6rwHHAtcC9wGn1hyxJGjYmoJIkDYDMvA54QpPxtwHHNBmfwOtqCE2SpPt5F1xJkiRJUi1MQCVJkiRJtTABlSRJkiTVwmtAJUmSpJr42BcNOxNQSZIkqSY+9kXDzi64kiRJkqRamIBKkiRJkmphAipJkiRJqoUJqCRJkiSpFt6ESOojV69cyaaJiaZl3hlPkiRJ850JqNRHNk1MeGc8SZKa8CCtNBhMQCVJktT3PEgrDQavAZUkSZIk1cIEVJIkSZJUCxNQSZIkSVItTEAlSZIkSbUwAZUkSZIk1cIEVJIkSZJUCxNQSZIkSVItTEAlSZIkSbUwAZUkSZIk1cIEVJIkSZJUCxNQSZIkSVItFsx1AJIkSRoeV69cyaaJiaZlC0dGOHzZspojklQnE1BJkiTVZtPEBEcsWtS0bNX4eM3RSKqbCaikltaMjbUs8yi1JEmSpssEVFJLWzZu9Ci1JEmSusabEEmSJEmSamECKkmSJEmqhQmoJEmSJKkWJqCSJEmSpFqYgEqSJEmSamECKkmSJEmqhQmoJEmSJKkWJqCSJEmSpFosmOsAJEmSJO3YmrGxlmULR0Y4fNmyGqORZsYEVJIkSZoHtmzcyBGLFjUtWzU+XnM00szMSQIaEX8MvBpI4GrgNGB/4EJgH+CHwCsyc1NEPABYARwB3Aa8NDNvmIu4JUmS1DvtzvBdPzbWMvmSNH/UnoBGxIHAG4HDMvOXEXERcDJwHPC+zLwwIj4EnA6cU/7ekZmPioiTgfcAL607bkmSJPVWuzN8q0dHa45GUi/M1U2IFgAPjIgFwG7ABuBZwOdK+fnACeX98WWYUn5MRESNsUqSJEmSuqD2BDQzfwH8HXAjVeJ5J7AKmMjMzWWydcCB5f2BwNry2c1l+n2nzjcilkfEaESMjtsHXpIkSZL6Tu0JaETsTXVW8+HAAcDuwLFNJs3Jj7Qp2zoi89zMXJqZSxd5fYAkSZIk9Z256IL7bOD6zBzPzF8Dnwd+BxgpXXIBDgLWl/frgIMBSvlewO31hixJkiRJmq25SEBvBJZFxG7lWs5jgGuAbwEnlmlOAb5U3l9chinl38zM7c6ASpIkSZL621xcA3oF1c2Efkj1CJadgHOBM4C3RMS1VNd4nlc+ch6wbxn/FuDMumOWJEmSJM3enDwHNDPfCbxzyujrgKc2mfY+4KQ64pIkSZIk9c5cPYZFkiRJkjRkTEAlSZIkSbUwAZUkSZIk1cIEVJIkSZJUizm5CZE0zK5euZJNExNNy64fG+OIRYtqjkjSoIiInYFR4BeZ+YKIeDhwIbAP1d3nX5GZmyLiAcAK4AjgNuClmXnDHIUtSRoiJqBSzTZNTLRMMlePjtYcjaQB8yZgDNizDL8HeF9mXhgRHwJOB84pf+/IzEdFxMllupfORcCSpOFiF1xJkgZARBwEPB/4aBkO4FlUz94GOB84obw/vgxTyo8p00uS1FMmoJIkDYb3A28DtpThfYGJzNxchtcBB5b3BwJrAUr5nWX67UTE8ogYjYjR8fHxXsUuSRoSdsGVJGmei4gXALdk5qqIOHpydJNJs4OybUdmngucC7B06dKm00iae2vGxlqWLRwZ4fBly2qMRmrNBFSSpPmDFgZgAAAgAElEQVTvacCLIuI4YFeqa0DfD4xExIJylvMgYH2Zfh1wMLAuIhYAewG31x+2pG7ZsnFjy3tMrLL3gvqIXXAlSZrnMvNPM/OgzFwCnAx8MzNfDnwLOLFMdgrwpfL+4jJMKf9mZnp2U5LUcyagkiQNrjOAt0TEtVTXeJ5Xxp8H7FvGvwU4c47ikyQNGbvgSpI0QDLzMuCy8v464KlNprkPOKnWwCRJwjOgkiRJkqSamIBKkiRJkmphAipJkiRJqoUJqCRJkiSpFiagkiRJkqRamIBKkiRJkmphAipJkiRJqoUJqCRJkiSpFiagkiRJkqRamIBKkiRJkmphAipJkiRJqoUJqCRJkiSpFiagkiRJkqRamIBKkiRJkmphAipJkiRJqoUJqCRJkiSpFiagkiRJkqRamIBKkiRJkmphAipJkiRJqoUJqCRJkiSpFiagkiRJkqRaLJjrACRJkiT1zpqxsZZlC0dGOHzZshqj0bAzAZUkSZIG2JaNGzli0aKmZavGx2uORsPOLriSJEmSpFqYgEqSJEmSamECKkmSJEmqhQmoJEmSJKkWJqCSJEmSpFqYgEqSJEmSamECKkmSJEmqhQmoJEmSJKkWHSWgEfG4XgciSZIqtruSpEHV6RnQD0XE9yPitREx0tOIJEmS7a4kaSB1lIBm5u8CLwcOBkYj4tMR8ZyeRiZJ0pCy3ZUkDaoFnU6Ymasj4h3AKPBB4EkREcDbM/PzvQpQkqRhZLur+ezqlSvZNDHRtOz6sTGOWLSo5ogk9YuOEtCIeDxwGvB84BLghZn5w4g4ALgcsCGUJKlLbHc1322amGiZZK4eHa05Gkn9pNMzoP8IfITqqOsvJ0dm5vpydFbSkFkzNtaybOHICIcvW1ZjNNLAsd2VJA2kThPQ44BfZuZvACJiJ2DXzLw3Mz/Rs+gk9a0tGze2PLq9any85mikgWO7K0kaSJ3eBfc/gAc2DO9WxkmSpO6z3ZUkDaROE9BdM/OeyYHyfrfehCRJ0tCz3ZUkDaROE9CNEfHkyYGIOAL4ZZvp24qIkYj4XET8NCLGIuK3I2KfiLgkIlaXv3uXaSMiPhgR10bEVY1xSJI0oLra7kqS1C86vQb0zcBnI2J9Gd4feOksvvcDwL9n5okRsZDqqO7bgUsz8+yIOBM4EzgDOBY4pLyOBM4pfyVJGlTdbnclSeoLHSWgmfmDiPgt4FAggJ9m5q9n8oURsSdwFHBqmfcmYFNEHA8cXSY7H7iMKgE9HliRmQmsLGdP98/MDTP5fkmS+l03211JkvpJp2dAAZ4CLCmfeVJEkJkrZvCdjwDGgX+JiCcAq4A3AQ+ZTCozc0NEPLhMfyCwtuHz68q4bRLQiFgOLAdYvHjxDMKSJKmvdKvdlSSpb3SUgEbEJ4BHAlcCvymjE5hJQ7gAeDLwhsy8IiI+QNXdtuXXNxmX243IPBc4F2Dp0qXblUuSNF90ud2VJKlvdHoGdClwWOkGO1vrgHWZeUUZ/hxVAnrzZNfaiNgfuKVh+oMbPn8QsB5JkgZXN9tdSZL6Rqd3wf1v4KHd+MLMvAlYGxGHllHHANcAFwOnlHGnAF8q7y8GXlnuhrsMuNPrPyVJA65r7a4kSf2k0zOg+wHXRMT3gV9NjszMF83we98AfKrcAfc64DSqZPiiiDgduBE4qUz7VeA44Frg3jKtJEmDrNvtriRJfaHTBPSsbn5pZl5J1b1oqmOaTJvA67r5/ZIk9bmz5joASZJ6oaMuuJn5beAGYJfy/gfAD3sYlyRJQ2sm7W5E7BoR34+IH0fETyLiL8r4h0fEFRGxOiI+U3ofEREPKMPXlvIlPa2UJEl0mIBGxB9Q3Szow2XUgcAXexWUJEnDbIbt7q+AZ2XmE4AnAs8r9054D/C+zDwEuAM4vUx/OnBHZj4KeF+ZTpKknur0JkSvA54G3AWQmauBB7f9hCRJmqlpt7tZuacM7lJeCTyLKpkFOB84obw/vgxTyo+JiGaPPpMkqWs6vQb0V5m5abJdiogFNHkWpyRJ6ooZtbsRsTOwCngU8E/AGmAiMzeXSdZRnU2l/F0LkJmbI+JOYF/g1inzXA4sB1i8ePHsaqWBcvXKlWyamGhadv3YGEcsWlRzRJLmg04T0G9HxNuBB0bEc4DXAl/uXViSJA21GbW7mfkb4IkRMQJ8AXhMs8nK32ZnO7dLcjPzXOBcgKVLl3rwWffbNDHRMslcPTpaczSS5otOu+CeCYwDVwN/SPVolHf0KihJkobcrNrdzJwALgOWASPlDCrAQcD68n4dcDDcf4Z1L+D2LsQuSVJLHZ0BzcwtwEfKS5Ik9dBM2t2IWAT8OjMnIuKBwLOpbiz0LeBE4ELgFOBL5SMXl+HLS/k3y6PPJEnqmY4S0Ii4nubdch7R9YikAeB1MZJmY4bt7v7A+eU60J2AizLzKxFxDXBhRPwV8CPgvDL9ecAnIuJaqjOfJ3ezDpIkNdPpNaBLG97vCpwE7NP9cKTB4HUxkmZp2u1uZl4FPKnJ+OuApzYZf1+ZryRJtenoGtDMvK3h9YvMfD/Vbd0lSVKX2e5KkgZVp11wn9wwuBPVkdk9ehKRJElDznZXkjSoOu2C+/cN7zcDNwAv6Xo0kiQJbHclSQOq07vgPrPXgUiSpIrtrqS6rBkba1m2cGSEw5ctqzEaDYNOu+C+pV15Zr63O+FIkiTbXUl12bJxY8sbJ64aH685Gg2D6dwF9ylUzwwDeCHwHWBtL4KSJGnI2e5KkgZSpwnofsCTM/NugIg4C/hsZr66V4FJkjTEbHclSQOpo8ewAIuBTQ3Dm4AlXY9GkiSB7a4kaUB1egb0E8D3I+ILQAIvBlb0LCpJkoab7a4kaSB1ehfcd0fE14Cnl1GnZeaPeheWJEnDy3ZXkjSoOu2CC7AbcFdmfgBYFxEP71FMkiTJdleSNIA6SkAj4p3AGcCfllG7AJ/sVVCSJA0z211J0qDq9Azoi4EXARsBMnM9sEevgpIkacjZ7kqSBlKnCeimzEyqGyEQEbv3LiRJkoae7a4kaSB1moBeFBEfBkYi4g+A/wA+0ruwJEkaara7kqSB1OldcP8uIp4D3AUcCvy/zLykp5FJkjSkbHclSYNqhwloROwMfD0znw3Y+EmS1EO2u5KkQbbDLriZ+Rvg3ojYq4Z4JEkaara7kqRB1lEXXOA+4OqIuIRyRz6AzHxjT6KSJGm42e5KkgZSpwnov5WXJEnqPdtdSdJAapuARsTizLwxM8+vKyBJkoaV7a4kadDt6BrQL06+iYh/7XEskiQNO9tdSdJA21ECGg3vH9HLQCRJku2uJGmw7SgBzRbvJUlS99nuSpIG2o5uQvSEiLiL6ojsA8t7ynBm5p49jU6SpOFiu6u+cvXKlWyamGhadv3YGEcsWlRzRJLmu7YJaGbuXFcgkiQNO9td9ZtNExMtk8zVo6M1RyNpEOyoC64kSZIkSV1hAipJkiRJqsWOrgGVpGlbMzbWsmzhyAiHL1tWYzSSJKnb2l0fbFuvdkxAJXXdlo0bW14ztGp8vOZoJEnSTLQ7oHz92BgnHnVU0zLberVjAipJkiRpO+0OKHsTKs2U14BKkiRJkmphAipJkiRJqoUJqCRJkiSpFiagkiRJkqRamIBKkiRJkmphAipJkiRJqoUJqCRJkiSpFj4HVJIkaUhdvXIlmyYmWpZfPzbW8jmQkjQTJqCSJElDatPERNsEc/XoaI3RSBoGdsGVJEmSJNXCBFSSJEmSVAsTUEmSJElSLUxAJUmSJEm1mLMENCJ2jogfRcRXyvDDI+KKiFgdEZ+JiIVl/APK8LWlfMlcxSxJkiRJmrm5PAP6JmCsYfg9wPsy8xDgDuD0Mv504I7MfBTwvjKdJEmSJGmemZMENCIOAp4PfLQMB/As4HNlkvOBE8r748swpfyYMr0kSZIkaR6Zq+eAvh94G7BHGd4XmMjMzWV4HXBgeX8gsBYgMzdHxJ1l+lsbZxgRy4HlAIsXL+5p8BK0f3i3D+6WJEmStld7AhoRLwBuycxVEXH05Ogmk2YHZVtHZJ4LnAuwdOnS7cqlbmv38G4f3C1JkiRtby7OgD4NeFFEHAfsCuxJdUZ0JCIWlLOgBwHry/TrgIOBdRGxANgLuL3+sCVJkiRJs1H7NaCZ+aeZeVBmLgFOBr6ZmS8HvgWcWCY7BfhSeX9xGaaUfzMzPcMpSVIREQdHxLciYiwifhIRbyrj94mIS8od5i+JiL3L+IiID5Y7zF8VEU+e2xpIkoZFPz0H9AzgLRFxLdU1nueV8ecB+5bxbwHOnKP4JEnqV5uBt2bmY4BlwOsi4jCqNvPScof5S9nahh4LHFJey4Fz6g9ZkjSM5uomRABk5mXAZeX9dcBTm0xzH3BSrYFJkjSPZOYGYEN5f3dEjFHdxO944Ogy2flUbe4ZZfyK0qNoZUSMRMT+ZT6SJPVMP50BlSRJsxQRS4AnAVcAD5lMKsvfB5fJ7r/DfNF49/mp81seEaMRMTo+Pt6rsCVJQ8IEVJKkARERDwL+FXhzZt7VbtIm45reXyEzz83MpZm5dJGPl5IkzZIJqCRJAyAidqFKPj+VmZ8vo2+OiP1L+f7ALWX85B3mJzXefV6SpJ4xAZUkaZ6LiKC6ad9YZr63oajxTvJT7zD/ynI33GXAnV7/KUmqw5zehEiSJHXF04BXAFdHxJVl3NuBs4GLIuJ04Ea23tTvq8BxwLXAvcBp9YYrSRpWJqCSJM1zmflfNL+uE+CYJtMn8LqeBiVJUhN2wZUkSZIk1cIEVJIkSZJUCxNQSZIkSVItTEAlSZIkSbUwAZUkSZIk1cIEVJIkSZJUCxNQSZIkSVItTEAlSZIkSbUwAZUkSZIk1cIEVJIkSZJUCxNQSZIkSVItTEAlSZIkSbUwAZUkSZIk1WLBXAcgSZIkaXCsGRtrW75wZITDly2rKRr1GxNQSZIkSV2zZeNGjli0qGX5qvHxGqNRv7ELriRJkiSpFiagkiRJkqRa2AVXUq28LkSSJGl4mYBKqpXXhUhSva5euZJNExNNy64fG2u7TZakbjMBlSRJGmCbJiZaJpmrR0drjkbSsPMaUEmSJElSLUxAJUmSJEm1MAGVJEmSJNXCBFSSJEmSVAsTUEmSJElSLUxAJUmSJEm1MAGVJEmSJNXCBFSSJEmSVAsTUEmSJElSLUxAJUmSJEm1WDDXAUj97OqVK9k0MdG07PqxMY5YtKjmiCRJ2p7tlaT5wgRUamPTxETLRnv16GjN0UiS1JztlaT5wi64kiRJkqRamIBKkiRJkmphAipJkiRJqoUJqCRJkiSpFiagkiRJkqRamIBKkiRJkmphAipJkiRJqoUJqCRJkiSpFiagkiRJkqRaLJjrACRJkiQNjzVjYy3LFo6McPiyZTVGo7qZgEqSJEmqzZaNGzli0aKmZavGx2uORnUzAZUkSZLUFzw7OvhMQCVJkiT1Bc+ODj5vQiRJkiRJqkXtZ0Aj4mBgBfBQYAtwbmZ+ICL2AT4DLAFuAF6SmXdERAAfAI4D7gVOzcwf1h23pHrY9UaSJGlwzUUX3M3AWzPzhxGxB7AqIi4BTgUuzcyzI+JM4EzgDOBY4JDyOhI4p/yVNIDseiNJkjS4ak9AM3MDsKG8vzsixoADgeOBo8tk5wOXUSWgxwMrMjOBlRExEhH7l/lIkiQNhatXrmTTxETTsuvHxloevJOkfjKnNyGKiCXAk4ArgIdMJpWZuSEiHlwmOxBY2/CxdWXcNgloRCwHlgMsXry4p3FLktRvIuJjwAuAWzLzcWWcl7cMkE0TEy2TzNWjozVHI0kzM2c3IYqIBwH/Crw5M+9qN2mTcbndiMxzM3NpZi5d5BFASdLw+TjwvCnjzqS6vOUQ4NIyDNte3rKc6vIWSZJ6bk4S0IjYhSr5/FRmfr6Mvjki9i/l+wO3lPHrgIMbPn4QsL6uWCVJmg8y8zvA7VNGH091WQvl7wkN41dkZSUwMtkGS5LUS7UnoKXbz3nAWGa+t6HoYuCU8v4U4EsN418ZlWXAnV7/KUlSR7a5vAXY0eUt24mI5RExGhGj494ITJI0S3NxBvRpwCuAZ0XEleV1HHA28JyIWA08pwwDfBW4DrgW+Ajw2jmIWZKkQdLR5S3gJS6SpO6ai7vg/hfNGz6AY5pMn8DrehqUJEmD6ebJO8d7eYskqR/M2U2IJElSz3l5iySpr8zpY1ikfuBz1SQNgoi4gOp52vtFxDrgnVSXs1wUEacDNwInlcm/SvUIlmupHsNyWu0BS5KGkgmohp7PVZM0CDLzZS2KvLxFktQ37IIrSZIkSaqFZ0AlSZIk9b01Y2MtyxaOjHD4smU1RqOZMgGVJEmS1Pe2bNzY8rKpVT6neN6wC64kSZIkqRYmoJIkSZKkWtgFV5IkqU/4aDBJg84EVJIkqU/4aDBJg84uuJIkSZKkWpiASpIkSZJqYQIqSZIkSaqFCagkSZIkqRYmoJIkSZKkWpiASpIkSZJqYQIqSZIkSaqFzwGVJEmq0dUrV7JpYqJp2fVjYy2fAypJg8AEVJIkqUabJiZaJpmrR0drjkaS6mUXXEmSJElSLUxAJUmSJEm1sAuupHljzdhYy7KFIyMcvmxZjdFIkiRpukxAJc0bWzZubHnd1Krx8ZqjkSRJ/cKD1POHCagkSZKkec2D1POH14BKkiRJkmrhGVANvHbPWwOfuSZJkiTVxQRUA6/d89bAZ65JkiRJdbELriRJkiSpFiagkiRJkqRamIBKkiRJkmphAipJkiRJqoUJqCRJkiSpFiagkiRJkqRa+BgWSZIkSQNrzdhYy7KFIyMcvmxZjdHIBFSSJKnLrl65kk0TE03Lrh8ba/t8akndtWXjxpa/uVXj4zVHIxNQSZKkLts0MdFyh3f16GjN0UhS//AaUEmSJElSLUxAJUmSJEm1sAuuJEmSJE3R7lpub140cyagkiRJkjRFu2u5vXnRzJmAShoI3mJdkiSp/5mAShoI3mJdUp3adc0DH7UizRftDmD7O+4NE1BJkqQmdvQszxOPOqrlZ33UijQ/tDuA7e+4N0xANRB84Lckqdt8lqckdZ8JqAaCOwmSJElS//M5oJIkSZKkWngGVJIkSZKmod3Ni8A78LdjAqp5w+s8NVM+okWSJHVTu5sXgXfgb8cEVPOG13lqpnxEizTc2h3A9CCUJNXLBFSSJA20dgcwPQglqRfsfdWaCaikoWYDIUmSus3eV62ZgEoaajYQ0nBrdxDK+wtIUveZgKqveKMhSVKd2h2E8v4Cknph2HtfzZsENCKeB3wA2Bn4aGaePcchqQe80ZAk1We+ta3eTEjSIBj23lfzIgGNiJ2BfwKeA6wDfhARF2fmNXMbmaRBNuxHKDXY5mPb2u4g5UXf+Y49aCSpiX47eDcvElDgqcC1mXkdQERcCBwP9G0jOQzarcw3rl3L4oMPnnaZOwnqJ+2OULbb2YWZb9D7rZHQQJuTtrXdOg4zX8/tSitpELQ7+N1uH7rdtrPf7gQemVn7l05XRJwIPC8zX12GXwEcmZmvb5hmObC8DB4K/KwLX70fcGsX5jNfDFN9h6muMFz1Haa6wnDVt1t1fVhmDv2Rrk7a1jK+2+3rMK2zMFz1Haa6wnDVd5jqCsNV39rb1vlyBjSajNsmc87Mc4Fzu/qlEaOZubSb8+xnw1TfYaorDFd9h6muMFz1Haa61mSHbSt0v30dtv/jMNV3mOoKw1XfYaorDFd956KuO9X5ZbOwDmg833wQsH6OYpEkaRDYtkqSajdfEtAfAIdExMMjYiFwMnDxHMckSdJ8ZtsqSardvOiCm5mbI+L1wNepbhX/scz8SQ1f3dUuvfPAMNV3mOoKw1XfYaorDFd9h6muPWfbWpthqu8w1RWGq77DVFcYrvrWXtd5cRMiSZIkSdL8N1+64EqSJEmS5jkTUEmSJElSLYY2AY2I50XEzyLi2og4s0n5URHxw4jYXJ6V1li2OCK+ERFjEXFNRCypK+6ZmGldI+KZEXFlw+u+iDih3uinb5b/27+NiJ+U/+0HI6LZYwr6xizr+p6I+O/yeml9Uc9cB/V9S/lNXhURl0bEwxrKTomI1eV1Sr2RT98s6/rvETEREV+pN+qZm2l9I+KJEXF5+d1eNV/W5UFl27pNuW3r1jLb1j5m27pNuW0rNbStmTl0L6qbLawBHgEsBH4MHDZlmiXA44EVwIlTyi4DnlPePwjYba7r1Ku6NkyzD3B7P9d1tvUFfgf4bpnHzsDlwNFzXace1fX5wCVUNyLbHRgF9pzrOnWhvs+cXEeBPwI+U97vA1xX/u5d3u8913XqRV3L8DHAC4GvzHVdavjfPho4pLw/ANgAjMx1nYbxNdv2BtvWOa9XL+qLbatta5+8ZlPXMmzb2qXYhvUM6FOBazPzuszcBFwIHN84QWbekJlXAVsax0fEYcCCzLykTHdPZt5bU9wzMeO6TnEi8LU+ryvMrr4J7Er1I30AsAtwc+9DnrHZ1PUw4NuZuTkzN1JtlJ5XR9Cz0El9v9Wwjq6keq4hwP8CLsnM2zPzDqodhH6u72zqSmZeCtxdV7BdMOP6Zub/ZObq8n49cAuwqLbI1ci2tYFt69YibFv7mW1rA9vWetrWYU1ADwTWNgyvK+M68WhgIiI+HxE/ioj/LyJ27nqE3TObujY6GbigKxH11ozrm5mXA9+iOsqzAfh6Zo51PcLumc3/9sfAsRGxW0TsR3UE7OAdfGauTbe+pwNfm+Fn59ps6jofdaW+EfFUqp3cNV2NTp2ybZ0+29b+Y9tq2zrJtpXetK3z4jmgPdDs2oNOn0ezAHg68CTgRuAzwKnAeV2JrPtmU9dqBhH7A4dTPSuu3824vhHxKOAxbD3adUlEHJWZ3+lWcF0247pm5jci4inA94Bxqi5Rm7sYWy90XN+I+H1gKfCM6X62T8ymrvPRrOtbtlOfAE7JzHZnnNQ7tq3TmYFtq21rf7BtbTahbevk+J60rcN6BnQd2x6ROghYP43P/qiczt4MfBF4cpfj66bZ1HXSS4AvZOavuxZV78ymvi8GVpauX/dQHQVa1uX4umlW/9vMfHdmPjEzn0O1kVrd5fi6raP6RsSzgT8DXpSZv5rOZ/vIbOo6H82qvhGxJ/BvwDsyc2WPY1Vrtq3TY9van2xbbVttW+lt2zqsCegPgEMi4uERsZCqC8zF0/js3hEx2Q/6WcA1PYixW2ZT10kvY350EYLZ1fdG4BkRsSAidqE6CtTP3YRmXNeI2Dki9i3vH091M4Vv9CzS7thhfSPiScCHqTaitzQUfR14bkTsHRF7A8+lv886zKau89GM61um/wKwIjM/W2PM2p5t6/TYtvYn21bbVtvWXret2Qd3aZqLF3Ac8D9U/Zn/rIx7V/kHADyF6sjBRuA24CcNn33O/9/e/cfWWdVxHH9/6IBON+f4XQOui7IsSgyZmmjYWDGLGEXNHMlAlzhA8UcUNRr+UBKXSMjMIkGHJIhjgwW2OZdNJArb1P1AdJuD0nU/JBqqsij4AwedFRC//vF8mz1ebm/b2/beUj6v5KTnOfec85zz3KbfPuc+z3OBLuAAsAY4pdnzGcO5tgNHgZOaPY+xni/F08JupwiMh4Cbmz2XMZxra87xEMVN5xc2ey6jNN/tFA+36Mx0X6nt1cDvMl3V7LmM8Vx3U1z+1Zfv/6XNns9YzRdYArxYKu98pfw+T8Q0wnjj2DqO0wjijWPrOE8jjDeOreM41Ttfxji2KndiZmZmZmZmNqZerZfgmpmZmZmZWYP5BNTMzMzMzMwawiegZmZmZmZm1hA+ATUzMzMzM7OG8AmomZmZmZmZNYRPQG3Ck7RQUkia3eyxDEZSh6Rjkh6VdFjS18d4fz+R9PpMn62jfZuk+0vb6yR1SfrSKIztqxXbD4+gr/WSzh/pmMzMrODYWnN/jq1mNfhrWGzCk/QDoA34WUQsG4X+WiLipREPrHrfHcBXIuIySa+l+N6lKyJi/xDaToqI/9S533bg/oi4YJjtVgAPRcSPJJ0D7ImIGaMxNkm9ETFlOG1q9DUfWBIRnxyN/szMXu0cW4e033YcW81exp+A2oQmaQpwEXANcEWpfIOk95e210haJKlF0gpJ+3K18VP5eoekX0i6l+JL0pG0RdJ+SQclXVvq6xpJj0vaIekOSbdm+ZmSNmXf+yRdVGvsEXEc2A+8SVKrpNWSDuQK7iXZ51JJGyX9GNiqwgpJ3Vl3cdZrk7RLUme+Ni/LeySdASzP/XRm+7WSPlya0z2SPlRlmIuABzK/FTgr+5iX879J0k7gC5I+KGlPjn+7pLP736PS3LryfVgOTM6+7sl6vflzoDl25D5/KOlIjlk5tt3AAkmTah1zMzMbnGOrY2uOzbHV6hMRTk4TNgFLgFWZfxiYk/mFwF2ZPwX4EzAZuBa4IctPBX4DzAQ6gOPAzFLfp+XPyUA3cDrwBqAHOA04meKP861Z715gbubfCByuMt4OitVSsr8e4K3Al4HVWT4b+CPQCiwFniyNZRGwDWgBzs56bdn+a1mnBZia+R7gDKAd6C6NYz6wJfPTgCeASRVjnQnsL21X9rEDuK20PZ0TV118AvhW5r8J3FKulz97K/bXO8gcO4BjwLkUi2u/6j/e2W4b8PZm/046OTk5vdITjq2OrSfaO7Y6DTt5xcImuiuBWzK/PrcfAX4KfEfSqcD7gF0R0SfpvcDbJF2ebaYB5wMvAHsj4olS39dJWpj587LeOcDOiPgHgKSNwKysswB4y4mFQ14naWpEPFcx5nmSHgX+CyyPiIOSbgRWAkTEEUl/KPW7rX9/wFxgXRSXMT2VK6TvBPYBd0o6mSL4ddY6aBGxU9J3JZ0FfATYFC+/zKcN+GutfoANpfy5wAZJbRT/mPQfywWUVtAj4plB+hxojs9SvEdPAkjqpAjcD2W7pyn+iRn0kpPGxKwAAALLSURBVCszM6vJsdWx1bHV6uYTUJuwJJ0OvAe4QFJQrOqFpOsj4t+SdgCXAouBdf3NgM9HxIMVfXVQrNKWtxcA746If2Vfrdl+ICdl/b5Bhr47Ii6rnE6N+sdL+ar1ImKXpIuBDwBrJa2IiLsHGcda4GMUAezqKq/3Ucy5lvLYVgI3R8R9efyWlcY8nJvRax2L50v5l/j/v3GtFGM2M7M6Obae4NgKOLZaHXwPqE1klwN3R8SMiGiPiPMoVgbn5uvrgauAeUB/UHwQ+EyuZiJplooHFlSaBjyTAXI28K4s3wvMlzQ974lYVGqzFfhc/4akC4cxl10UAQtJsyguM/rtAPUWq7jf5kzgYmCvpBnA0xFxB7AKmFPR7jlgakXZGuCLABFxsMq+HqdYBR2qacDRzH+8VF55XKZn9sX+96FC1TkOYf+zgGrzMDOzoXNsdWwtc2y1YfMJqE1kVwKbK8o2AR/N/FaKP7DbI+KFLPs+cAh4RFI3cDvVrxR4AJgkqQv4BvBrgIg4CtwE7AG2Z1/Hss11wDvyYQCHgE8PYy63AS2SDlBcerM0Ip6vUm8z0AU8BvwcuD4i/kJxD0dnXn60CPh2uVFE/B34ZT58YEWWPQUcBlZXG1AUD3L4vaQ3D3EOy4CNknYDfyuV3whMz30/BlyS5d8DuvoflDCEOQ4oH8rQFxF/HuJYzcysOsdWx1bAsdXq569hMRtlkqZERG+u0m4G7oyIymA97kl6DcVTCedExLEB6iykePjADQ0d3DCp+O60ZyNiVbPHYmZmw+fYOv44tlq9/Amo2ehbljfpd1NclrSlyeMZNkkLgCPAyoECJEAG/55GjWsE/gnc1exBmJlZ3Rxbxx/HVquLPwE1MzMzMzOzhvAnoGZmZmZmZtYQPgE1MzMzMzOzhvAJqJmZmZmZmTWET0DNzMzMzMysIXwCamZmZmZmZg3xP1CxNicp9BZWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "vario = GSLIB.make_variogram(nug=0.0,nst=1,it1=1,cc1=1.0,azi1=0.0,hmaj1=250,hmin1=250)\n",
    "n_eff = n_effective(df,'X','Y',seed=73073,nreal=1000,vario=vario)\n",
    "\n",
    "L = 10000                                      # set the number of realizations for uncertainty calculation\n",
    "mean = np.zeros(L); spatial_mean = np.zeros(L) # declare arrays to hold the realizations of the statistics\n",
    "P10 = np.zeros(L)                          \n",
    "P50 = np.zeros(L); P90 = np.zeros(L)                         \n",
    "for l in range(0, L):                          # loop over realizations\n",
    "    samples = random.choices(df['Porosity'].values, weights=None, cum_weights=None, k=len(df))\n",
    "    spatial_samples = random.choices(df['Porosity'].values, weights=None, cum_weights=None, k=int(n_eff))\n",
    "\n",
    "    mean[l] = np.average(samples)\n",
    "    spatial_mean[l] = np.average(spatial_samples)\n",
    "    \n",
    "plt.subplot(121)\n",
    "GSLIB.hist_st(mean,0.16,0.22,False,False,50,None,'Average Porosity (fraction)','Bootstrap Uncertainty in Porosity Average')\n",
    "\n",
    "plt.subplot(122)\n",
    "GSLIB.hist_st(spatial_mean,0.16,0.22,False,False,50,None,'Average Porosity (fraction)','Spatial Bootstrap Uncertainty in Porosity Standard Deviation')\n",
    "\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.2, hspace=0.2)\n",
    "plt.show()   \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the expansion of variance when we account for spatial continuity.  Experiment with the variogram model and rerun this result.  \n",
    "\n",
    "* What happens when the range is increased?\n",
    "* What happens when the nugget is added and increased?\n",
    "\n",
    "#### Comments\n",
    "\n",
    "This was a basic demonstration of spatial bootstrap with a comparison to traditional bootstrap with a developed function to calculate the number of effective data with LU-simulation.\n",
    "\n",
    "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",
    "#### The Author:\n",
    "\n",
    "### Michael Pyrcz, Associate Professor, University of Texas at Austin \n",
    "*Novel Data Analytics, Geostatistics and Machine Learning Subsurface Solutions*\n",
    "\n",
    "With over 17 years of experience in subsurface consulting, research and development, Michael has returned to academia driven by his passion for teaching and enthusiasm for enhancing engineers' and geoscientists' impact in subsurface resource development. \n",
    "\n",
    "For more about Michael check out these links:\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",
    "#### Want to Work Together?\n",
    "\n",
    "I hope this content is helpful to those that want to learn more about subsurface modeling, data analytics and machine learning. Students and working professionals are welcome to participate.\n",
    "\n",
    "* Want to invite me to visit your company for training, mentoring, project review, workflow design and / or consulting? I'd be happy to drop by and work with you! \n",
    "\n",
    "* Interested in partnering, supporting my graduate student research or my Subsurface Data Analytics and Machine Learning consortium (co-PIs including Profs. Foster, Torres-Verdin and van Oort)? My research combines data analytics, stochastic modeling and machine learning theory with practice to develop novel methods and workflows to add value. We are solving challenging subsurface problems!\n",
    "\n",
    "* I can be reached at mpyrcz@austin.utexas.edu.\n",
    "\n",
    "I'm always happy to discuss,\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}