SubsurfaceDataAnalytics_Confidence_Hypothesis.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",
    "### Confidence Intervals and Hypothesis Testing for Subsurface Data Analytics in Python \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": [
    "##### Reporting Uncertainty and Significance\n",
    "\n",
    "With confidence intervals and hypothesis testing we have the opportunity to report uncertainty and to report significance in our statistics. \n",
    "\n",
    "* report **uncertainty and significance** with our results\n",
    "\n",
    "This is a tutorial / demonstration of **Confidence Intervals and Hypothesis Testing in Python** for Subsurface Modeling.  In Python, the SciPy package, specifically the Stats functions (https://docs.scipy.org/doc/scipy/reference/stats.html) provide excellent tools for efficient use of statistics. \n",
    "\n",
    "This tutorial includes basic, typical confidence interval and hypothesis testing methods that would commonly be required for Engineers and Geoscientists including:\n",
    "\n",
    "1. Student-t confidence interval for the mean\n",
    "2. Student-t hypothesis test for difference in means (pooled variance)\n",
    "3. Student-t hypothesis test for difference in means (difference variances), Welch's t Test\n",
    "3. F-distribution hypothesis test for difference in variances \n",
    "\n",
    "#### Additional Resources\n",
    "\n",
    "These workflows are based on standard methods with their associated limitations and assumptions. For more information see:\n",
    "\n",
    "* [Confidence Intervals Lecture](https://www.youtube.com/watch?v=oaXCcTWcU04)\n",
    "\n",
    "* [Hypothesis Testing_Lecture](https://www.youtube.com/watch?v=rvt9UM148tQ) \n",
    "\n",
    "* Also, see the lecture and workflow on Bootstrap https://git.io/fhgUW for a general, empirical approach to assess uncertianty in statistics.\n",
    "\n",
    "I have provided various workflows for subsurface data analytics, geostatistics and machine learning:\n",
    "\n",
    "* [Python](https://git.io/fh4eX)\n",
    "\n",
    "* [Excel](https://github.com/GeostatsGuy/LectureExercises/blob/master/Lecture7_CI_Hypoth_eg_R.xlsx) \n",
    "* [R](https://github.com/GeostatsGuy/LectureExercises/blob/master/Lecture7_CI_Hypoth_eg.R)  \n",
    "\n",
    "and all of my University of Texas at Austin \n",
    "\n",
    "* [Lectures](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig/featured?view_as=subscriber)\n",
    "\n",
    "##### Caveats\n",
    "\n",
    "I have not included all the details, specifically the test assumptions in this document.  These are included in the accompanying course notes, recorded lectures listed above.\n",
    "\n",
    "#### Workflow Goal\n",
    "\n",
    "0. Introduction to Python in Jupyter including setting a working directory, loading data into a Pandas DataFrame.\n",
    "1. Learn the basics for working with confidence intervals and hypothesis testing in Python.  \n",
    "2. Demonstrate the efficiency of using Python and SciPy package for statistical analysis.\n",
    "3. Learn how to quantify uncertainty and significance in samples.\n",
    "\n",
    "#### Objective \n",
    "\n",
    "I want to provide hands-on experience with building subsurface modeling workflows. Python provides an excellent vehicle to accomplish this. I have coded a package called GeostatsPy with GSLIB: Geostatistical Library (Deutsch and Journel, 1998) functionality that provides basic building blocks for building subsurface modeling workflows. \n",
    "\n",
    "The objective is to remove the hurdles of subsurface modeling workflow construction by providing building blocks and sufficient examples. This is not a coding class per se, but we need the ability to 'script' workflows working with numerical methods.    \n",
    "\n",
    "#### Getting Started\n",
    "\n",
    "Here's the steps to get setup in Python with the GeostatsPy package:\n",
    "\n",
    "1. Install Anaconda 3 on your machine (https://www.anaconda.com/download/). \n",
    "2. From Anaconda Navigator (within Anaconda3 group), go to the environment tab, click on base (root) green arrow and open a terminal. \n",
    "3. In the terminal type: pip install geostatspy. \n",
    "4. Open Jupyter and in the top block get started by copy and pasting the code block below from this Jupyter Notebook to start using the geostatspy functionality. \n",
    "\n",
    "You will need to copy the data file to your working directory.  They are available here:\n",
    "\n",
    "* Tabular data - PorositySample2Units.csv at https://git.io/fhrM8\n",
    "\n",
    "#### Load the required libraries\n",
    "\n",
    "The following code loads the required libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "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": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os                                              # to set current working directory \n",
    "import numpy as np                                     # arrays and matrix math\n",
    "import scipy.stats as st                               # statistical methods\n",
    "import pandas as pd                                    # DataFrames\n",
    "import matplotlib.pyplot as plt                        # plotting"
   ]
  },
  {
   "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"
   ]
  },
  {
   "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).  Also, in this case make sure to place the required (see below) data file in this directory.  When we are done with this tutorial we will write our new dataset back to this directory.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(\"C:\\PGE337\")                                  # set the working directory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Loading Data\n",
    "\n",
    "Let's load the provided dataset. 'PorositySamples2Units.csv' is available at https://github.com/GeostatsGuy/GeoDataSets. It is a comma delimited file with 20 porosity measures from 2 rock units from the subsurface:\n",
    "\n",
    "* porosity (fraction) \n",
    "\n",
    "We load it with the pandas 'read_csv' function into a data frame we called 'df' and then preview it by printing a slice and by utilizing the 'head' DataFrame member function (with a nice and clean format, see below).\n",
    "\n",
    "We load it with the pandas 'read_csv' function into a data frame we called 'df' and then preview it by printing a slice and by utilizing the 'head' DataFrame member function (with a nice and clean format, see below).\n",
    "\n",
    "**Python Tip: using functions from a package** just type the label for the package that we declared at the beginning:\n",
    "\n",
    "```python\n",
    "import pandas as pd\n",
    "```\n",
    "\n",
    "so we can access the pandas function 'read_csv' with the command: \n",
    "\n",
    "```python\n",
    "pd.read_csv()\n",
    "```\n",
    "\n",
    "but read csv has required input parameters. The essential one is the name of the file. For our circumstance all the other default parameters are fine. If you want to see all the possible parameters for this function, just go to the docs [here](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.read_csv.html).  \n",
    "\n",
    "* The docs are always helpful\n",
    "* There is often a lot of flexibility for Python functions, possible through using various inputs parameters\n",
    "\n",
    "also, the program has an output, a pandas DataFrame loaded from the data.  So we have to specficy the name / variable representing that new object.\n",
    "\n",
    "```python\n",
    "df = pd.read_csv(\"PorositySample2Units.csv\")  \n",
    "```\n",
    "\n",
    "Let's run this command to load the data and then look at the resulting DataFrame to ensure that we loaded it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "scrolled": false
   },
   "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>Unit1_Por</th>\n",
       "      <th>Unit2_Por</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.21</td>\n",
       "      <td>0.20</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.17</td>\n",
       "      <td>0.26</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.15</td>\n",
       "      <td>0.20</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.20</td>\n",
       "      <td>0.19</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.19</td>\n",
       "      <td>0.13</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Unit1_Por  Unit2_Por\n",
       "0       0.21       0.20\n",
       "1       0.17       0.26\n",
       "2       0.15       0.20\n",
       "3       0.20       0.19\n",
       "4       0.19       0.13"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv(\"PorositySample2Units.csv\")           # read a .csv file in as a DataFrame\n",
    "df = df.rename(columns={'X1': 'Unit1_Por','X2':'Unit2_Por'})   # rename variables for clarity\n",
    "#print(df.iloc[0:5,:])                                 # display first 4 samples in the table as a preview\n",
    "df.head()                                              # we could also use this command for a table preview "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Summary Statistics\n",
    "\n",
    "It is useful to review the summary statistics of our loaded DataFrame.  That can be accomplished with the 'describe' DataFrame member function.  We transpose to switch the axes for ease of visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>count</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>min</th>\n",
       "      <th>25%</th>\n",
       "      <th>50%</th>\n",
       "      <th>75%</th>\n",
       "      <th>max</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Unit1_Por</th>\n",
       "      <td>20.0</td>\n",
       "      <td>0.1645</td>\n",
       "      <td>0.027810</td>\n",
       "      <td>0.11</td>\n",
       "      <td>0.1500</td>\n",
       "      <td>0.17</td>\n",
       "      <td>0.19</td>\n",
       "      <td>0.21</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Unit2_Por</th>\n",
       "      <td>20.0</td>\n",
       "      <td>0.2000</td>\n",
       "      <td>0.045422</td>\n",
       "      <td>0.11</td>\n",
       "      <td>0.1675</td>\n",
       "      <td>0.20</td>\n",
       "      <td>0.23</td>\n",
       "      <td>0.30</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           count    mean       std   min     25%   50%   75%   max\n",
       "Unit1_Por   20.0  0.1645  0.027810  0.11  0.1500  0.17  0.19  0.21\n",
       "Unit2_Por   20.0  0.2000  0.045422  0.11  0.1675  0.20  0.23  0.30"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.describe().transpose()                              # visualize summary statistics "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Extract Porosity from Well1 and Well2\n",
    "\n",
    "Here we extract the X1 and X2 unit porosity samples from the DataFrame into separate arrays called 'Por1' and 'Por2' for convenience.\n",
    "\n",
    "**Python Tip: extracting data from DataFrames**\n",
    "\n",
    "We can extract the data for a single feature with this command:\n",
    "\n",
    "```python\n",
    "1D_series = df['feature_name']\n",
    "```\n",
    "\n",
    "The 'series' retains information about the feature that was included in the DataFrame including the name and the indexing. This is fine, but some methods don't work with series so we can also extract the data as a 1D ndarray.  By adding the '.values' the series is converted to a 1D array.\n",
    "\n",
    "```python\n",
    "1D_ndarray = df['feature_name'].values\n",
    "```\n",
    "\n",
    "So "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "por1 = df['Unit1_Por'].values                                 # extract well 1 porosity to a ndarray\n",
    "por2 = df['Unit2_Por'].values                                 # extract well 2 porosity to a ndarray          "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could also work with the DataFrame directly and just use this command everytime we use a function that requires a ndarray.  For example:\n",
    "\n",
    "```python \n",
    "np.average(df['Unit1_Por'].values)\n",
    "```\n",
    "\n",
    "would return the average over unit 1.  I like to extract the ndarray so we can use more compact code:\n",
    "\n",
    "```python\n",
    "np.average(por)\n",
    "```\n",
    "\n",
    "to get the same result.\n",
    "\n",
    "**Coding Tip** the ndarrays are **shallow copies**, that means that if you change the array the DataFrame will also be changed. In other words that DataFrame and ndarray remained linked together.  \n",
    "\n",
    "* you can modify either the ndarray or DataFrame and they remain consistent\n",
    "\n",
    "#### Well Distributions\n",
    "\n",
    "We should first visualize the distributions.  It is convenient to plot the distributions on top of each other for comparision.\n",
    "\n",
    "**Coding Tip** this is our first time to use MatplotLib.  This exaple below makes the binned histogram and cumulative distribution function for the unit 1 and 2 porosity measures.\n",
    "\n",
    "* we use $plt$ to access the matplotlib.plotly functionality as imported above\n",
    "\n",
    "* after we declare a subplot all commands append and modify to the current plot\n",
    "\n",
    "* subplot(number of rows, number of columns, plot index) is the format\n",
    "\n",
    "```python \n",
    "plt.subplot(121)                # the first plot for a 1 x 2 matrix of plots\n",
    "```\n",
    "\n",
    "* we superimpose two histograms on top of eachother with different colours and transparency \n",
    "\n",
    "```python \n",
    "plt.hist(por1, facecolor='red',bins=np.linspace(0.0,0.4,20),alpha=0.2,density=False,edgecolor='black',label='Unit 1')\n",
    "plt.hist(por2, facecolor='blue',bins=np.linspace(0.0,0.4,20),alpha=0.2,density=False,edgecolor='black',label = 'Unit 2')\n",
    "```\n",
    "\n",
    "* then we add labels, set the axes, and include a legend\n",
    "\n",
    "Later on will will use the histogram function from reimplimented in geostatspy from GSLIB for simple, easy histograms.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5QAAAGWCAYAAAAHYcxqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XuYXWV58P/vnYQQQwIJZAwqSIKFWKQiJlJEsaN4AAtoKVqxJ1FMrXiq2Fc8vEKp7c+29vXAa4sRPIAHBFolWBgPwEZ9MRwSNSCIEEw0oGMERhgCTCa5f3/snbAzzGR21uy1D5Pv57rmyjo8a637njXsh3uvtZ4VmYkkSZIkSTtrSrsDkCRJkiR1JwtKSZIkSVIhFpSSJEmSpEIsKCVJkiRJhVhQSpIkSZIKsaCUJEmSJBViQSk1ICJ+EhG97Y4DICLeHxHntzuORkTE5yPiw+2OQ5LUeexbi7FvVaexoFTHi4i1EfFIRAxGRH9EfC4iZrUyhsx8VmZWavGcHRFfLLqvWj4vHbHsDRHx/QZj+efMPK223YKIyIiYtoPjHRoR34yI30ZEx7x4NiL+OCK+HxEDEfHriPhMRMxud1yStCuwb31CLPatUkEWlOoWJ2TmLOC5wPOAD+7sDnbUMUxym4BLgDe1O5AR9gI+DDwV+H1gP+Df2hqRJO1a7FuLs2+Vaiwo1VUy8x7gKuBQgIh4akQsj4j7I+KuiHjz1ra1bzsvi4gvRsSDwBsiYveI+HhE3Fv7+XhE7F5rPy8ivlH7Vu/+iPheREyprVsbES+NiGOB9wN/VvtW98cR8ZqIWFkfZ0ScERFfL5Jj3Tejfx0Rv6h9+/mBEXlt/Rb3u7V/B2rxPH+U39kdmXkB8JMGj/+JiPhlRDwYESsj4ugRx74kIi6MiIdqtystqVt/eESsqq37KjBjrONk5pczsy8zN2bmA8BngBc0EqMkqXnsW+1bpYmwoFRXiYj9gVcCP6wt+gqwnuo3cScD/xwRx9Rt8irgMmAO8CXgA8CRwHOAw4AjePwb2TNq++oB5lPt3La7jSUz+4B/Br6ambMy8zBgObAwIn6/rulfABdNMN0XAouAY4APjdj/Vi+q/TunFs8PJnhMgJuo/n72Br4MXBoR9Z3XicDFVH+ny4H/CxAR04GvU817b+BS4E934rgvosGOWZLUPPatT2DfKu0EC0p1i69HxADwfeA6qp3b/lQ7hvdm5qOZ+SPgfOAv67b7QWZ+PTO3ZOYjwJ8D52TmbzJzA/APde03AU8BDsjMTZn5vcwc97mIzHwM+CrVjo6IeBawAPjGBHP+h8x8JDN/DPyYaiddusz8Ymbel5nDmfnvwO5UO9+tvp+ZV2bmZqod3Na4jgR2Az5e+/1dRrUDHVdEvAz4a+BDTUtEkjQe+1b7VmnCLCjVLV6dmXMy84DMfGutA3sqcH9mPlTXbh3wtLr5X47Yz1NrberbP7U2/W/AXcC3IuLuiDhzJ+L7AvD6iAiqnegltc5wNMNUO4d6u1HtdOv9um56I9CSwRJqtxTdHhG/q/2Pxl7AvB3ENSOqz9A8FbhnxP8o1P+uxzrekVS/rT05M3828QwkSQ2yb7VvlSbMglLd7F5g79h+9LKnA/fUzY/8FvRe4IAR7e8FyMyHMvOMzDwQOAF494hbfMbaJ5m5AhgCjgZez45vyfkF1W9Z6y2kgQ6ikVgmovZMx3uB1wJzM3MO8DsgGtj8V8DTah3/Vk8f53iHU721542ZeXWxqCVJTWTfat8q7RQLSnWtzPwlcD3w/0XEjIh4NtXR1r60g82+AnwwInoiYh7V20C+CBARx0fE79U+tB8ENtd+RuoHFmwdVKDOhVSfeRjOzB0NU/5V4F0R8cyoWgK8keqzEztrA7AFOHCsBrVjzACm1+ZnbB0sYRSzqX7LuwGYFhEfAvZsMJYf1LZ9R0RMi4iTqD5HM1ZchwJ9wNsz84oGjyFJKpF9K2DfKu0UC0p1u1OofiN5L/A14KzM/PYO2n8YuBlYDdwCrKotAzgI+A4wSPUD/D+2vh9rhEtr/94XEavqll9EdYS88QYM+AzwOeAKqt9QXgh8oDYowU7JzI3APwH/L6oj6B05SrMDgEd4/KH8R4A7xtjlN6mO9Pczqt/qPsoTb20aK5Yh4CTgDcADwJ8B/72DTc6gOkjDBbVR9AYjwoEDJKn97FvtW6WGRQPPRUtqQEQ8CfgN8NzMvLPd8UiS1O3sW6XO5xVKqXn+FrjJDk+SpKaxb5U63LQydx4RfwecRvXh5luAUzPz0TKPKbVDRKyl+nD9q9sciqQuFRGfBY4HfpOZh46yPoBPUH1f4EbgDZm5amQ7abKwb5W6Q2lXKCPiacA7gCW1jnEq8Lqyjie1U2YuqA27/sPxW0vSqD4PHLuD9cdRfR7tIGAp8J8tiElqG/tWqTuUfcvrNOBJtffozKQ2hLQkSdpeZn4XuH8HTV4FXJhVK4A5EfGU1kQnSdLoSrvlNTPviYiPUn0v0CPAtzLzWyPbRcRSqt+0MmPGjMVPf/oOX63T8bZs2cKUKd39aKo5dAZz6Azm0Bl+9rOf/TYze9odR5s9je1HhlxfW/arkQ3tWzuPOXSGkTls2bIFumyAyqSxF1h2gi2Zo77VM6YEuaW7fu8jTYYcAO7++d0T7l9LKygjYi7Vb1MXAgPApRHxF5n5xfp2mbkMWAawaNGivOOOsUZc7g6VSoXe3t52hzEh5tAZzKEzmENniIgiLyefbEb7f8hR/2/GvrXzmENnGJnDyr4+Fvd013dVlf5+eufPb3cYDem7/k565h72hOX9M9cxf+MBbYioeSZDDgBL/vKQCfevZX7N9FLg55m5ITM3UX1nzlElHk+SpMlsPbB/3fx++CiJJKnNyiwofwEcGREzayPTHQPcXuLxJEmazJYDfxVVRwK/y8wn3O4qSVIrlfkM5Q0RcRmwChgGfkjt9htJkrS9iPgK0AvMi4j1wFnAbgCZeR5wJdVXhtxF9bUhp7YnUkmSHlfqeygz8yyqHWJhmzZtYv369Tz6aHe8vnKvvfbi9tubeyF2xowZ7Lfffuy2225N3a8kqXNk5injrE/g9GYcy77VvlWSmqXUgrIZ1q9fz+zZs1mwYAHVO2c720MPPcTs2bObtr/M5L777mP9+vUsXLiwafuVJO267FvtWyWpWTp+7OdHH32UffbZpys6vDJEBPvss0/XfIssSep89q32rZLULB1fUAK7bIe31a6evySp+Xb1vmVXz1+SmqUrCkpJkiRJUufp+GcoR7plxQqGBgaatr/pc+bwB0ceOeb6tWvXcvzxx3PrrbduW3b22Wcza9Ys3vOe94y53c0338yFF17IJz/5SSqVCtOnT+eoo574Gs6f/vSnnHrqqaxatYp/+qd/2uE+JUkqg32rJKmorisohwYGWNzT07T9rdywoWn7qrdkyRKWLFkCQKVSYdasWaN2envvvTef/OQn+frXv15KHJIkjce+VZJUlLe8TlBvby/vfe97OeKIIzj44IO5/vrrgWpHd/zxx7N27VrOO+88Pvaxj/Gc5zyH733ve9tt/+QnP5nnPe95DlsuSVKNfaskdY+uu0LZiYaHh7nxxhu58sor+chHPsIrXvGKbesWLFjAW97ylnFv45EkSY+zb5Wk7uAVynGMNQpc/fKTTjoJgMWLF7Nu3bqWxCVJUreyb5WkycOCchz77LMPDzzwwHbL7r//fubNm7dtfvfddwdg6tSpbN68uaXxSZLUbexbJWnysKAcx6xZs3jKU57C1VdfDVQ7vL6+Pl74whc2vI/Zs2fz0EMPlRWiJEldxb5VkiaPrnuGcvqcOU0dPW76nDnjtrnwwgs5/fTTOeOMMwA466yzeMYzntHwMU444QROPvlkLr/8cs4991yOPvrobet+/etfs2TJEh588EGmTJnCxz/+cW677Tb23HPPnU9GkqQC7Fs1WY18Jc7GoSFW9vVtm5++aVM7wmqaFavXMDC4pd1hjGnT8PR2h6AW6LqCckfvtSrLIYccwrXXXjvqukqlsm163rx5296p1dvbS29vLwAHH3wwq1evHnX7fffdl/Xr1zc1XkmSdoZ9qyarka/EqfT3N/UVOe02MLiFnrmHtTsM7eK85VWSJEmSVIgFpSRJkiSpEAtKSZIkSVIhFpSSJEmSpEIsKCVJkiRJhVhQSpIkSZIK6brXhqxYcQsDA0NN29+cOdM58sg/GHP92rVrOf7447cNWQ5w9tlnM2vWLN7znveMud3NN9/MhRdeyCc/+UkqlQrTp0/nqKOOekK7L33pS/zLv/wLUH3R83/+539y2GEO/yxJah37VklSUV1XUA4MDNHTs7hp+9uwYWXT9lVvyZIlLFmyBKi+T2vWrFmjdnoLFy7kuuuuY+7cuVx11VUsXbqUG264oZSYJEkajX2rJKkob3mdoN7eXt773vdyxBFHcPDBB3P99dcD1Y7u+OOPZ+3atZx33nl87GMf4znPeQ7f+973ttv+qKOOYu7cuQAceeSRvohZkrTLs2+VpO7RdVcoO9Hw8DA33ngjV155JR/5yEd4xStesW3dggULeMtb3jLubTwAF1xwAccdd1zZ4UqS1PHsWyWpO1hQjiMixl1+0kknAbB48WLWrVtX6DjXXnstF1xwAd///vcLbS9JUrewb5WkycNbXsexzz778MADD2y37P7772fevHnb5nfffXcApk6dyubNm3f6GKtXr+a0007j8ssvZ5999plYwJIkdTj7VkmaPCwoxzFr1iye8pSncPXVVwPVDq+vr48XvvCFDe9j9uzZPPTQQ6Ou+8UvfsFJJ53ERRddxMEHH9yUmCVJ6mT2rZI0eXTdLa9z5kxv6uhxc+ZMH7fNhRdeyOmnn84ZZ5wBwFlnncUznvGMho9xwgkncPLJJ3P55Zdz7rnncvTRR29bd84553Dffffx1re+FYBp06Zx880372QWkiQVZ98qSSqq6wrKHb3XqiyHHHII11577ajrKpXKtul58+Zte6dWb28vvb29ABx88MGsXr161O3PP/98zj///KbGK0nSzrBvlSQV5S2vkiRJkqRCuu4KpSRJkgRwy4oVDA0MjLl++qZNLYym+VasXsPA4Jbtlg3N352+NXcCsGl4/NvLpbJ1RUGZmWMOMb4ryMx2hyBJmmTsW+1bJ4OhgQEW9/S0O4zSDAxuoWfuYdst65+2jp65z2xTRNITdfwtrzNmzOC+++7bZT/4M5P77ruPGTNmtDsUSdIkYd9q3ypJzdLxVyj3228/1q9fz4YNG9odSkMeffTRpndQM2bMYL/99mvqPiVJuy77VvtWSWqWji8od9ttNxYuXNjuMBpWqVQ4/PDD2x2GJEljsm+VJDVLabe8RsSiiPhR3c+DEfGuso4nSZIkSWqt0q5QZuYdwHMAImIqcA/wtbKOJ0mSJElqrVYNynMMsCYz17XoeJIkSZKkkrXqGcrXAV8ZbUVELAWWAvT09FCpVFoUUjkGBwfNoQN0cg4PP/wImzePP7Ji5jDLl185brupU4M99nhSM0Jruk4+D40yB0mSpLGVXlBGxHTgROB9o63PzGXAMoBFixZlb29v2SGVqlKpYA7t18k59PWtZN99F4/brr+/wvz5veO227BhJb294++vHTr5PDTKHCRJksbWiltejwNWZWZ/C44lSZIkSWqRVhSUpzDG7a6SJEmSpO5VakEZETOBlwH/XeZxJEmSJEmtV+ozlJm5EdinzGNIkiRJktqjVa8NkSRJkiRNMhaUkiRJkqRCLCglSZIkSYWU/h5KSZIkaVe0YvUaBga3FN5+0/D0JkYjlcOCUpIkSSrBwOAWeuYe1u4wpFJ5y6skSZIkqRALSkmSJElSIRaUkiRJkqRCLCglSZIkSYVYUEqSJEmSCrGglCRJkiQVYkEpSZIkSSrEglKSJEmSVIgFpSRJkiSpEAtKSZIkSVIhFpSSJEmSpEIsKCVJkiRJhVhQSpLUISLi2Ii4IyLuiogzR1n/9Ii4NiJ+GBGrI+KV7YhTkqStLCglSeoAETEV+BRwHHAIcEpEHDKi2QeBSzLzcOB1wH+0NkpJkrZnQSlJUmc4ArgrM+/OzCHgYuBVI9oksGdtei/g3hbGJ0nSE0xrdwCSJAmApwG/rJtfD/zhiDZnA9+KiLcDewAvHW1HEbEUWArQ09NDpVJpdqwtNTg4aA4doBNz2Dg0RKW/v+H2g8PDO9V+oobm707/tHVN3efwlCH6ZzZ3n61mDpOLBaUkSZ0hRlmWI+ZPAT6fmf8eEc8HLoqIQzNzy3YbZS4DlgEsWrQoe3t7y4i3ZSqVCubQfp2Yw8q+Phb39DTcvtLfT+/8+SVGtL2+NXfSM/eZTd1n/8x1zN94QFP32WrmMLl4y6skSZ1hPbB/3fx+PPGW1jcBlwBk5g+AGcC8lkQnSdIoLCglSeoMNwEHRcTCiJhOddCd5SPa/AI4BiAifp9qQbmhpVFKklTHglKSpA6QmcPA24BvArdTHc31JxFxTkScWGt2BvDmiPgx8BXgDZk58rZYSZJaxmcoJUnqEJl5JXDliGUfqpu+DXhBq+OSJGksXqGUJEmSJBViQSlJkiRJKsSCUpIkSZJUiAWlJEmSJKkQC0pJkiRJUiEWlJIkSZKkQiwoJUmSJEmFWFBKkiRJkgqxoJQkSZIkFVJqQRkRcyLisoj4aUTcHhHPL/N4kiRJkqTWmVby/j8B9GXmyRExHZhZ8vEkSZIkSS1SWkEZEXsCLwLeAJCZQ8BQWceTJEmSJLVWmVcoDwQ2AJ+LiMOAlcA7M/Ph+kYRsRRYCtDT00OlUikxpPINDg6aQwfo5BweevB+Nj7463HbTdltmHVr/mfcdpt5jErloWaE1nSdfB4aZQ6SJEljK7OgnAY8F3h7Zt4QEZ8AzgT+d32jzFwGLANYtGhR9vb2lhhS+SqVCubQfp2cw7kf/jRHLXjRuO36p61j/tAB47a7fu13ec1rT2pGaE3XyeehUeYgSZI0tjIH5VkPrM/MG2rzl1EtMCVJkiRJk0BpBWVm/hr4ZUQsqi06BritrONJkiRJklqr7FFe3w58qTbC693AqSUfT5IkSZLUIqUWlJn5I2BJmceQJEmSJLVHmc9QSpIkSZImMQtKSZIkSVIhFpSSJEmSpEIsKCVJkiRJhVhQSpIkSZIKsaCUJEmSJBViQSlJkiRJKsSCUpIkSZJUiAWlJEmSJKkQC0pJkiRJUiEWlJIkSZKkQqa1OwBJkiRpNLesWMHQwMCY66dv2lTq8VesXsPA4JbC228ant7EaKTOZEEpSZKkjjQ0MMDinp62HX9gcAs9cw9r2/GlbuAtr5IkSZKkQiwoJUmSJEmFWFBKkiRJkgqxoJQkSZIkFWJBKUmSJEkqxIJSkiRJklSIBaUkSZIkqRALSkmSJElSIRaUkiRJkqRCLCglSZIkSYVYUEqSJEmSCrGglCRJkiQVYkEpSZIkSSrEglKSpCaKiJsj4vSImNvuWCRJKpsFpSRJzfU64KnATRFxcUS8IiKi3UFJklSGae0OQJKkySQz7wI+EBH/Gzge+CywJSI+C3wiM+9va4BSC92yYgVDAwOFt5++aVMTo3miFavXMDC4Zcz1m4anl3p8aTKwoJQkqcki4tnAqcArgf8CvgS8ELgGeE4bQ5NaamhggMU9Pe0OY0wDg1vomXtYu8OQupoFpSRJTRQRK4EB4ALgzMx8rLbqhoh4QfsikySp+SwoJUlqrtdk5t31CyJiYWb+PDNPaldQkiSVwUF5JElqrssaXCZJUtfzCqUkSU0QEc8EngXsFRH1VyL3BGa0JypJkspVakEZEWuBh4DNwHBmLinzeJIktdEiqqO6zgFOqFv+EPDmtkQkSVLJWnGF8sWZ+dsWHEeSpLbJzMuByyPi+Zn5g3bHI0lSK3jLqyRJTRAR/ysz/xV4fUScMnJ9Zr6jDWFJklSqsgvKBL4VEQl8OjOXjWwQEUuBpQA9PT1UKpWSQyrX4OCgObTJIw8/TG7eDMBwJlcuXz7hfQ4NDTF9enNfarzP/rPpn7lu3HbDU4Yaajf/Gft07Pnq1r+leuagnXB77d+b2xqFJEktVHZB+YLMvDcingx8OyJ+mpnfrW9QKzKXASxatCh7e3tLDqlclUoFc2iPlX19LN53XwAq/f30zp8/4X1efNVVvPq44ya8n3pnXvZlXvPy08Zt1z9zHfM3HjBuuzVrv8trTzm5GaE1Xbf+LdUzBzUqM6+o/fuFdsciSVKrlFpQZua9tX9/ExFfA44AvrvjrSRJ6j4RcQXVO3NGlZkntjAcSZJaorSCMiL2AKZk5kO16ZcD55R1PEmS2uyj7Q5AkqRWK/MK5XzgaxGx9Thfzsy+Eo8nSVLbZOZ17Y5BkqRWK62gzMy7gcPK2r8kSZ0kIi7JzNdGxC1sf+trAJmZz25gH8cCnwCmAudn5kdGafNa4OzaMX6cma9vRvySJBXha0MkSWqOd9b+Pb7IxhExFfgU8DJgPXBTRCzPzNvq2hwEvI/qoHcP1Aa9kySpbaa0OwBJkiaDzPxV7d91wGNU79J5NvBYbdl4jgDuysy7M3MIuBh41Yg2bwY+lZkP1I71m2bFL0lSEV6hlCSpiSLiNOBDwDVUb3c9NyLOyczPjrPp04Bf1s2vB/5wRJuDa8f4f1Rviz17tPEJfMdz59lVc9g4NESlv7+cgAoYHB7eLp6h+bvTP62R73s6R6Pvqe5k5jC5WFBKktRcfw8cnpn3AUTEPsD1wHgFZYyybORrSKYBBwG9wH7A9yLi0Mwc2G4j3/HccXbVHFb29bG4p6ecgAoY+Z7qvjV30jP3mW2MaOc1+p7qTmYOk4u3vEqS1FzrgYfq5h9i+yuPO9pu/7r5/YB7R2lzeWZuysyfA3dQLTAlSWoLr1BKktQEEfHu2uQ9wA0RcTnVK4yvAm5sYBc3AQdFxMLaPl4HjBzB9evAKcDnI2Ie1Vtg725C+JIkFWJBKUlSc8yu/bum9rPV5Y1snJnDEfE24JtUn4/8bGb+JCLOAW7OzOW1dS+PiNuAzcDfb721VpKkdrCglCSpCTLzH5qwjyuBK0cs+1DddALvrv1IktR2FpSSJDVRRPQA/wt4FjBj6/LMfEnbgpIkqSQOyiNJUnN9CfgpsBD4B2At1ecjJUmadCwoJUlqrn0y8wJgU2Zel5lvBI5sd1CSJJXBW14lSWquTbV/fxURf0z11R/7tTEeSZJKY0EpSVJzfTgi9gLOAM4F9gT+rr0hSZJUjoYKyog4NDNvLTsYSZK6XWZ+ozb5O+DF7YxFkqSyNfoM5XkRcWNEvDUi5pQakSRJXSwiDoyIKyLitxHxm4i4PCIObHdckiSVoaGCMjNfCPw5sD9wc0R8OSJeVmpkkiR1py8DlwD7Ak8FLgW+0taIJEkqScOjvGbmncAHgfcCfwR8MiJ+GhEnlRWcJEldKDLzoswcrv18Ech2ByVJUhkafYby2cCpwB8D3wZOyMxVEfFU4AfAf5cXoiRJnS8i9q5NXhsRZwIXUy0k/wz4n7YFJklSiRod5fX/Ap8B3p+Zj2xdmJn3RsQHS4lMkqTuspJqARm1+b+pW5fAP7Y8IkmSStZoQflK4JHM3AwQEVOAGZm5MTMvKi06SZK6RGYubHcMkiS1WqMF5XeAlwKDtfmZwLeAo8oISpKkbhURuwF/C7yotqgCfDozN7UtKGkXtWL1GgYGt2ybH5q/O31r7tw2v2l4ejvCkiaVRgvKGZm5tZgkMwcjYmZJMUmS1M3+E9gN+I/a/F/Wlp3WtoikXdTA4BZ65h62bb5/2jp65j6zjRFJk0+jBeXDEfHczFwFEBGLgUfG2UaSpF3R8zLzsLr5ayLix22LRpKkEjVaUL4LuDQi7q3NP4XqqHWSJGl7myPiGZm5BiAiDgQ2tzkmSZJK0VBBmZk3RcQzgUVUR6/7qc+CSJI0qr+n+uqQu6n2mQdQffWWJEmTTqNXKAGeByyobXN4RJCZF5YSlSRJXag2CvojwEFs/yXsY20NTJKkkjRUUEbERcAzgB/x+G07CVhQSpJUk5lbIuLfM/P5wOp2xyNJUtkavUK5BDgkM7PMYCRJmgS+FRF/Cvy3/aYkabJrtKC8FdgX+FWJsUiSNBm8G9gDGI6IR6ne9pqZuWd7w5IkqfkaLSjnAbdFxI3AtudAMvPEUqKSJKlLZebsdscgSVKrNFpQnl1mEJIkdbuIeDLwfuD3qD4/+ZHMfLC9UUmSVK4pjTTKzOuAtcButembgFUlxiVJUre5EHgYOBeYDXyyveFIklS+Rkd5fTOwFNib6mivTwPOA44pLzRJkrrKvpn5gdr0NyPCL14lSZNeo7e8ng4cAdwAkJl31m7tkSRJVRERc6kOwgMwtX4+M+9vW2SSJJWk0YLyscwciqj2kRExjep7KMcVEVOBm4F7MvP4QlFKktT59gJW8nhBCY8/HpLAgS2PSJKkkjVaUF4XEe8HnhQRLwPeClzR4LbvBG4HHC5dkjRpZeaCdscgSVKrNTQoD3AmsAG4Bfgb4Ergg+NtFBH7AX8MnF80QEmSJElSZ2roCmVmbgE+U/vZGR8H/hfV0e5GFRFLqQ74Q09PD5VKZScP0VkGBwfNoU02Dg1R6e8HYHB4eNv0ROz+jGc0ZT/1/uAlh9M/c9247YanDDXUbu7T92T58iubEdo2U6cGe+zxpAnvp1v/luqZgyRJ0tgaHeX154zyzGRmjvk8SEQcD/wmM1dGRO9Y7TJzGbAMYNGiRdnbO2bTrlCpVDCH9ljZ18finh4AKv399M6fP+F9XrxqFb3HHTfh/dQ78+Krec3LTxu3Xf/MdczfeMC47a74zo2c9s4PNSO0bTZsWElv7+IJ76db/5bqmYMkSdLYGn2Gcknd9AzgNVRfIbIjLwBOjIhX1rbZMyK+mJl/sfNhSpLUPSLihcBBmfm5iOgBZmXmz9sdlyRJzdbQM5SZeV/dzz2Z+XHgJeNs877M3K82SMHrgGssJiVJk11EnAW8F3hfbdFuwBfbF5EkSeVp9JbX59bNTqF6xXLM5yIlSdqF/QlwOLVXhmTmvRFhnylJmpQaveX13+umh4G1wGsbPUhmVoBKo+14fgcSAAAckElEQVQlSepiQ5mZEZEAEbFHuwOSJKksjY7y+uKyA5EkaZK4JCI+DcyJiDcDb2TnR0mXJKkrNHrL67t3tD4z/09zwpEkqbtl5kcj4mXAg8Ai4EOZ+e02hyVJUil2ZpTX5wHLa/MnAN8FfllGUJIkdauI+DvgUotISdKuoNGCch7w3Mx8CCAizqbaWY7/Mj1JknYtewLfjIj7gYuByzKzv80xSZJUioZeGwI8HRiqmx8CFjQ9GkmSulxm/kNmPgs4HXgqcF1EfKfNYUmSVIpGr1BeBNwYEV8DkuqQ6BeWFpUkSd3vN8CvgfuAJ7c5FkmSStHoKK//FBFXAUfXFp2amT8sLyxJkrpTRPwt8GdAD3AZ8ObMvK29UUmSVI5Gr1ACzAQezMzPRURPRCzMzJ+XFZgkSV3qAOBdmfmjdgciSVLZGn1tyFlUR3pdBHwO2A34IvCC8kKTJKl7RMSemfkg8K+1+b3r12fm/W0JTJKkEjV6hfJPgMOBVQCZeW9EzC4tKkmSus+XgeOBlVTHG4i6dQkc2I6gJEkqU6MF5VBmZkQkQETsUWJMkiR1ncw8vvbvwnbHIklSqzRaUF4SEZ8G5kTEm4E3Ap8pLyxJkrpTRFydmceMt0yaDG5ZsYKhgYEx10/ftKmF0Uhqh0ZHef1oRLwMeJDqc5QfysxvlxqZJEldJCJmUB3Abl5EzOXxW173pPo+SmnSGRoYYHFPT7vDkNRG4xaUETEV+GZmvhSwiJQkaXR/A7yLavG4kscLygeBT7UrKEmSyjRuQZmZmyNiY0TslZm/a0VQkiR1m8z8BPCJiHh7Zp7b7ngkSWqFRp+hfBS4JSK+DTy8dWFmvqOUqCRJ6lKZeW5EHAocAsyoW35h+6KSJKkcjRaU/1P7kSRJO1B7d3Mv1YLySuA44PuABaUkadLZYUEZEU/PzF9k5hdaFZAkSV3uZOAw4IeZeWpEzAfOb3NMkiSVYso467++dSIi/qvkWCRJmgweycwtwHBE7An8BjiwzTFJklSK8W55jbppO0NJksZ3c0TMofq+5pXAIHBje0OSJKkc4xWUOca0JEkaRWa+tTZ5XkT0AXtm5up2xiRJUlnGKygPi4gHqV6pfFJtmtp8ZuaepUYnSVKXiIjn7mhdZq5qZTySJLXCDgvKzJzaqkAkSepy/76DdQm8pFWBSJLUKo2+NkSSJO1AZr643TFIktRqFpSSJDVRRPzVaMsz0/dQSpImHQtKSZKa63l10zOAY4BVgAWlJGnSsaCUJKmJMvPt9fMRsRdwUSPbRsSxwCeAqcD5mfmRMdqdDFwKPC8zb55YxJIkFTel3QFIkjTJbQQOGq9RREwFPgUcBxwCnBIRh4zSbjbwDuCGJscpSdJO8wqlJElNFBFX8Pi7m6dQLQ4vaWDTI4C7MvPu2n4uBl4F3Dai3T8C/wq8pykBS5I0ARaUkiQ110frpoeBdZm5voHtngb8sm5+PfCH9Q0i4nBg/8z8RkSMWVBGxFJgKUBPTw+VSqXB0DvT4OCgOXSA0XLYODREpb+/PQEBDz/yGJs3j70+58+gf9q6bfPDU4bon7lu7A26gDl0hsmQQ7NYUEqS1ESZeR1AROxJrZ+NiL0z8/5xNo3RdrdtZcQU4GPAGxqIYRmwDGDRokXZ29vbSOgdq1KpYA7tN1oOK/v6WNzT056AgL7r72TfuYftuNHQ45P9M9cxf+MB5QZVMnPoDJMhh2axoJQkqYlqVwf/EXgE2EK1UEzgwHE2XQ/sXze/H3Bv3fxs4FCgEhEA+wLLI+JEB+aRJLWLBaUkSc3198CzMvO3O7ndTcBBEbEQuAd4HfD6rSsz83fAvK3zEVEB3mMxKUlqJ0d5lSSpudZQHdl1p2TmMPA24JvA7cAlmfmTiDgnIk5scoySJDVFaVcoI2IG8F1g99pxLsvMs8o6niRJHeJ9wPURcQPw2NaFmfmO8TbMzCuBK0cs+9AYbXsnFqYkSRNX5i2vjwEvyczBiNgN+H5EXJWZK0o8piRJ7fZp4BrgFqrPUEqSNGmVVlBmZgKDtdndaj859haSJE0Kw5n57nYHIUlSK5Q6KE9ETAVWAr8HfCozbxilje/K6jDdmsN9Dw5yz4OPAJC7BcvX/GLC+8ynL2j6+7X+4CWHN/Teokbfb/T8lxxEf3+lCZHVHXt4Y1P+Brr1b6meOaiAa2t92xVsf8vreK8NkSSp65RaUGbmZuA5ETEH+FpEHJqZt45o47uyOky35nDuhz/NUQteBED/tHXMH5r4u4Eu/c75nPLO14/fcCecefHVvOblp43brtH3G11xzY2c9s5TmhHaNhs2rKS3d/GE99Otf0v1zEEFbP3QeF/dskZeGyJJUtdpyWtDMnOgNrz5scCt4zSXJKlrZebCdscgSVKrlDnKaw+wqVZMPgl4KfAvZR1PkqROEBF/NdryzLyw1bFIklS2Mq9QPgX4Qu05yilU36f1jRKPJ0lSJ3he3fQM4BhgFWBBKUmadMoc5XU1cHhZ+5ckqRNl5tvr5yNiL+CiNoUjSVKpprQ7AEmSJrmNwEHtDkKSpDK0ZFAeSZJ2FRFxBY+/d3kKcAhwSfsikiSpPBaUkiQ110frpoeBdZm5vl3BSJJUJgtKSZKaICJ+D5ifmdeNWH50ROyemWvaFJpU2C0rVjA0MADAxqEhVvb1bbd++qZN7QhLUgexoJQkqTk+Drx/lOWP1Nad0NpwpIkbGhhgcU8PAJX+/m3TkrSVg/JIktQcC2ojnG8nM28GFrQ+HEmSymdBKUlSc8zYwbontSwKSZJayIJSkqTmuCki3jxyYUS8CVjZhngkSSqdz1BKktQc7wK+FhF/zuMF5BJgOvAnbYtKkqQSWVBKktQEmdkPHBURLwYOrS3+n8y8po1hSZJUKgtKSZKaKDOvBa5tdxySJLWCz1BKkiRJkgqxoJQkSZIkFWJBKUmSJEkqxIJSkiRJklSIBaUkSZIkqRALSkmSJElSIRaUkiRJkqRCLCglSZIkSYVYUEqSJEmSCpnW7gAkSZI0Oa1YvYaBwS2l7X/T8PTS9i2pMRaUkiRJKsXA4BZ65h7W7jAklchbXiVJkiRJhVhQSpIkSZIKsaCUJEmSJBViQSlJkiRJKsSCUpIkSZJUiAWlJEmSJKkQC0pJkiRJUiEWlJIkSZKkQiwoJUmSJEmFWFBKkiRJkgqxoJQkSZIkFWJBKUmSJEkqpLSCMiL2j4hrI+L2iPhJRLyzrGNJkiRJklpvWon7HgbOyMxVETEbWBkR387M20o8piRJkiSpRUq7QpmZv8rMVbXph4DbgaeVdTxJkiRJUmuVeYVym4hYABwO3DDKuqXAUoCenh4qlUorQirN4ODgpMzhkYcfJjdvbupxhoaGmD59etP2t8/+s+mfuQ6A4SlD26Yn4g9ecjiV/v4J72fkPhuJrdEcnv+Sg+jvrzQhsrpjD29syt/xZP3vodtMhhwkdaYVq9cwMLhlzPWbhpvXz0vqTKUXlBExC/gv4F2Z+eDI9Zm5DFgGsGjRouzt7S07pFJVKhUmYw4r+/pYvO++TT3OxVddxauPO65p+zvzsi/zmpefBkD/zHXM33jAhPf53Wu+zZ+/8/UT3k+9My++elucO9JoDldccyOnvfOUZoS2zYYNK+ntXTzh/UzW/x66zWTIQVJnGhjcQs/cw9odhqQ2KnWU14jYjWox+aXM/O8yjyVJkiRJaq0yR3kN4ALg9sz8P2UdR5IkSZLUHmVeoXwB8JfASyLiR7WfV5Z4PEmSJElSC5X2DGVmfh+IsvYvSZIkSWqvUp+hlCRJkiRNXhaUkiRJkqRCLCglSZIkSYVYUEqSJEmSCrGglCRJkiQVYkEpSZIkSSrEglKSJEmSVIgFpSRJkiSpEAtKSZIkSVIhFpSSJHWIiDg2Iu6IiLsi4sxR1r87Im6LiNURcXVEHNCOOCVJ2sqCUpKkDhARU4FPAccBhwCnRMQhI5r9EFiSmc8GLgP+tbVRSpK0PQtKSZI6wxHAXZl5d2YOARcDr6pvkJnXZubG2uwKYL8WxyhJ0namtTsASZIEwNOAX9bNrwf+cAft3wRcNdqKiFgKLAXo6emhUqk0KcT2GBwcNIeSPPLww+TmzWOuj0wq/f0ADA4Pb5veamj+7vRPW1dqjM00PGWI/pndE+9ozKEzTIYcmsWCUpKkzhCjLMtRG0b8BbAE+KPR1mfmMmAZwKJFi7K3t7dJIbZHpVLBHMqxsq+Pxfvu21DbSn8/vfPnb7esb82d9Mx9ZhmhlaJ/5jrmb+zuR4/NoTNMhhyaxYJSkqTOsB7Yv25+P+DekY0i4qXAB4A/yszHWhSbJEmj8hlKSZI6w03AQRGxMCKmA68Dltc3iIjDgU8DJ2bmb9oQoyRJ27GglCSpA2TmMPA24JvA7cAlmfmTiDgnIk6sNfs3YBZwaUT8KCKWj7E7SZJawlteJUnqEJl5JXDliGUfqpt+acuDkiRpB7xCKUmSJEkqxIJSkiRJklSIBaUkSZIkqRALSkmSJElSIRaUkiRJkqRCLCglSZIkSYVYUEqSJEmSCrGglCRJkiQVYkEpSZIkSSrEglKSJEmSVIgFpSRJkiSpkGntDkCSJEnF3LJiBUMDA4W3n75p0w7Xr1i9hoHBLQAMzd+dvjV3brd+0/D0wseWNDlYUEqSJHWpoYEBFvf0lLb/gcEt9Mw9DID+aevomfvM0o4lqTt5y6skSZIkqRALSkmSJElSIRaUkiRJkqRCLCglSZIkSYWUVlBGxGcj4jcRcWtZx5AkSZIktU+ZVyg/Dxxb4v4lSZIkSW1UWkGZmd8F7i9r/5IkSZKk9mr7eygjYimwFKCnp4dKpdLegCZocHBwUuZw34OD3PPgI009zmNP2Z/la37RtP0tOvrZ9M9cB8DwlKFt0xPdZzNj3LrPRmJrNIfnv+Qg+vsrTYjscYODD7B8ef+E95M5zPLlVwIwNDTE9OnNewH21KnBHns8qWn7A3j44UfYvDm3W1afQxFlxLmzJsPnkiRJ6kxtLygzcxmwDGDRokXZ29vb3oAmqFKpMBlzOPfDn+aoBS9q6nEu+vLHOfPUdzVvfxc/vr/+meuYv/GACe/zcxd/rakxwvZx7kijOVxxzY2c9s5TmhHaNqtWXcxxx/3JhPfT319h/vxeAK666mKOO+7VE97nVhs2rKS3d3HT9gfQ17eSfffdfp/1ORRRRpw7azJ8LkmSpM7kKK+SJEmSpEIsKCVJkiRJhZT52pCvAD8AFkXE+oh4U1nHkiRJkiS1XmnPUGZmcx/qkiRJkiR1FG95lSRJkiQVYkEpSZIkSSqk7a8NkSRJ2lXdsmIFQwMDhbe//Za72TCj+Pbj2TTcvPcHS5qcLCglSZLaZGhggMU9PYW33zBjgJ65hzUxIknaOd7yKkmSJEkqxIJSkiRJklSIBaUkSZIkqRALSkmSJElSIRaUkiRJkqRCLCglSZIkSYVYUEqSJEmSCrGglCRJkiQVYkEpSZIkSSrEglKSJEmSVIgFpSRJkiSpEAtKSZIkSVIhFpSSJEmSpEIsKCVJkiRJhVhQSpIkSZIKsaCUJEmSJBUyrd0BSJIk7apW3/4LNtw5UHj7TcPTmxiNJO08C0pJkqQ2GXw46VlwWLvDkKTCvOVVkiRJklSIBaUkSZIkqRALSkmSJElSIRaUkiRJkqRCLCglSZIkSYVYUEqSJEmSCrGglCRJkiQVYkEpSZIkSSrEglKSJEmSVIgFpSRJkiSpEAtKSZIkSVIhFpSSJEmSpEJKLSgj4tiIuCMi7oqIM8s8liRJ3W68fjMido+Ir9bW3xARC1ofpSRJjyutoIyIqcCngOOAQ4BTIuKQso4nSVI3a7DffBPwQGb+HvAx4F9aG6UkSdsr8wrlEcBdmXl3Zg4BFwOvKvF4kiR1s0b6zVcBX6hNXwYcExHRwhglSdpOZGY5O444GTg2M0+rzf8l8IeZ+bYR7ZYCS2uzhwK3lhJQ68wDftvuICbIHDqDOXQGc+gMizJzdruDKFMj/WZE3Fprs742v6bW5rcj9mXf2nnMoTOYQ2cwh84x4f51WrMiGcVo35g+oXrNzGXAMoCIuDkzl5QYU+nMoTOYQ2cwh84wWXJodwwt0Ei/ad/apcyhM5hDZzCHztGM/rXMW17XA/vXze8H3Fvi8SRJ6maN9Jvb2kTENGAv4P6WRCdJ0ijKLChvAg6KiIURMR14HbC8xONJktTNGuk3lwN/XZs+Gbgmy3p2RZKkBpR2y2tmDkfE24BvAlOBz2bmT8bZbFlZ8bSQOXQGc+gM5tAZzKELjNVvRsQ5wM2ZuRy4ALgoIu6iemXydQ3sejL87syhM5hDZzCHzjAZcoAm5FHaoDySJEmSpMmtzFteJUmSJEmTmAWlJEmSJKmQlhSUEXFsRNwREXdFxJmjrN89Ir5aW39DRCyoW/e+2vI7IuIVrYh3NEVziIgFEfFIRPyo9nNeq2Ovi3G8HF4UEasiYrj2PrT6dX8dEXfWfv565LatMsEcNtedh7YOENVAHu+OiNsiYnVEXB0RB9St65ZzsaMcOuJcNJDDWyLillqc34+IQ+rWdctn06g5dNNnU127kyMiI2JJ3bKOOA/tMBn61los9q/d8Zne8f2rfWvXnAf71hZoad+amaX+UB1YYA1wIDAd+DFwyIg2bwXOq02/DvhqbfqQWvvdgYW1/UwtO+Ym57AAuLXVMRfMYQHwbOBC4OS65XsDd9f+nVubnttNOdTWDbb7POxEHi8GZtam/7bu76mbzsWoOXTKuWgwhz3rpk8E+mrT3fTZNFYOXfPZVGs3G/gusAJY0knnoVN/b3R439qEPLrmbxj7107Jwb61M3Kwb+2AHGrtmtK3tuIK5RHAXZl5d2YOARcDrxrR5lXAF2rTlwHHRETUll+cmY9l5s+Bu2r7a7WJ5NApxs0hM9dm5mpgy4htXwF8OzPvz8wHgG8Dx7Yi6BEmkkMnaSSPazNzY212BdX30UF3nYuxcugUjeTwYN3sHjz+Avmu+WzaQQ6dopHPV4B/BP4VeLRuWaech3aYDH0r2L9202d6p/ev9q2dwb61M7S0b21FQfk04Jd18+try0Ztk5nDwO+AfRrcthUmkgPAwoj4YURcFxFHlx3sGCbyu+ym87AjMyLi5ohYERGvbm5oO2Vn83gTcFXBbcsykRygM85FQzlExOkRsYbqB+47dmbbFphIDtAln00RcTiwf2Z+Y2e3ncQmQ99Kg7HYv5ZvMvSv9q1ddB7sW0vX0r61tPdQ1hntW8SRVfxYbRrZthUmksOvgKdn5n0RsRj4ekQ8a8Q3G60wkd9lN52HHXl6Zt4bEQcC10TELZm5pkmx7YyG84iIvwCWAH+0s9uWbCI5QGeci4ZyyMxPAZ+KiNcDH6T6UvmuOg9j5NAVn00RMQX4GPCGnd12kpsMfSvYv3bKuZgM/at9axedB/vW0rW0b23FFcr1wP518/sB947VJiKmAXtRfWFzI9u2QuEcapeL7wPIzJVU70M+uPSIn2giv8tuOg9jysx7a//eDVSAw5sZ3E5oKI+IeCnwAeDEzHxsZ7ZtgYnk0CnnYmd/lxcDW7/x7arzUGdbDl302TQbOBSoRMRa4EhgeW3wgE45D+0wGfpWGozF/rV8k6F/tW/tovNQx761HK3tW7P8h0KnUX24eSGPPxT6rBFtTmf7B+4vqU0/i+0fCr2b9jycO5EcerbGTPXB2HuAvTsxh7q2n+eJgwb8nOqD6nNr092Ww1xg99r0POBORnk4uVPyoNoJrAEOGrG8a87FDnLoiHPRYA4H1U2fANxcm+6mz6axcui6z6Za+wqPDxzQEeehHT8NnvuO7lubkEfX/Q1j/9ruvyX71s7Iwb61A3IY0b7CBPrWViX1SuBntf8APlBbdg7Vb1YAZgCXUn3o80bgwLptP1Db7g7guFafkInmAPwp8JPaiVkFnNDBOTyP6rcSDwP3AT+p2/aNtdzuAk7tthyAo4BbaufhFuBN7cqhwTy+A/QDP6r9LO/CczFqDp10LhrI4RO1/35/BFxL3YdxF302jZpDN302jWhbodbpddJ56MTfG13Qt04kj276G8b+tVNysG/tjBzsWzsghxFtK0ygb43aRpIkSZIk7ZRWPEMpSZIkSZqELCglSZIkSYVYUEqSJEmSCrGglCRJkiQVYkEpSZIkSSrEglK7jIjYHBE/iohbI+LSiJhZ4rHOqb14mIh4184eK6quiYg9a/PviIjbI+JLTYjtDRHx1Lr58yPikIL7eltEnDrRmCRJ3cv+ddu+7V+1S/K1IdplRMRgZv7/7d1baBxVHMfx769JoUlaLCEiVgqlqPgiNVpFvD9UVFSoEIQqRaj4IIKYWgRv0BdBfbAPQiltpUshCIpKQ0Wpl6wWJZUkWlKhCmrqS/GCBZsY0dq/D3O2jOnuZrNK1nR+Hzjsmck5Z/6TJfvPXM7O0lQfAEYj4qUG+onsb+V0k9udIHu2z89z6HMnsC4i+tPyUbLnAH03o117RJyaYzxlYEtEjMylX42xOoFPIqL3345lZmYLk/PrmT5lnF+tgHyF0orqIHAxgKTN6azqEUmPpXWr0hnL7WQPpl0paYOk8dTuhdSuTVIprRuXVElQJUl9kh4FVgBDkoYkPShpWyUISQ9JqpZ07wf2pTY7gNXAoKR+SVsl7ZR0ANibYj0oaSyV63LjP5HiOizpeUl9wFpgIJ1N7pBUlrQ2tT9rH9P6SUnPpXGGJV0AEBG/AROSrvlP3hUzM1vonF+dX61oIsLFpRAFmEyv7WTJ5GHgKmAc6AKWAl8CvcAq4DRwbeqzAvgeOD/1/xBYn/q/l9vG8vRaAvpSfQLoSfUu4BtgcVr+FLi8SqzHgGW55fwYW4FRoCMtdwJLUv0SYCTV70jjd6bl7vRaJjujS3651j6mNgHcneovAs/k+j8NPN7q99fFxcXFpTXF+dX51aXYxVcorUg6JH0BjJB9sL8C3AC8FRFTETEJvAncmNofi4jhVL8aKEfET5HdAjMA3AR8C6yW9LKk24Ff6wUQEVNkieQuSZeRJb7xKk27I+JknaEGI2I61RcDuySNA68Dlfka64A9kZ3lJCJ+qRdbnX0E+APYn+qjZP8QVPxIlizNzKyYnF/rc361c1p7qwMwm0fTEXFFfkWav1HLVL5ptQYRcULSGuA24BHgXmDTLHHsBp4CjgJ7arQ5JWlR1J5Xko+tH/gBWEN2G/vvuZjnMkm63u/iz4iojPUX//zsWAJMn93FzMwKwvm1PudXO6f5CqUV3cfAekmdkrqAe8jmf8x0CLhZUo+kNmAD8JGkHmBRRLwBPAtcWaXvSWBZZSEiDgErgfuAV2vE9RXZvI5GnAccT8lxI9CW1h8ANil9A56k7mrxzLaPDWz/UuBIg7GamVkxOL/Oso8NbN/51RYEH1BaoUXEGNl8jM/IPvB3R8TnVdodB54EhoDDwFhE7AMuAsrpVp9SajPTTuAdSUO5da+RfXvbiRqhvQ3c0uBubAcekDRMlnymUszvAoPASIpvS2pfAnZUvjSggX2czfXA+w3GamZmBeD86vxqxeHHhpi1gKT9wLaI+KDGzy8E9kbErfMb2dxI6gU2R8TGVsdiZmbm/Go2/3yF0mweSVou6Wuy+SZVkx2cOZu5S+nBy/9jPWS3IpmZmbWM86tZ6/gKpZmZmZmZmTXFVyjNzMzMzMysKT6gNDMzMzMzs6b4gNLMzMzMzMya4gNKMzMzMzMza4oPKM3MzMzMzKwpfwPoGYZhdfmYVQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.subplot(121)\n",
    "plt.hist(por1, facecolor='red',bins=np.linspace(0.0,0.4,20),alpha=0.2,density=False,edgecolor='black',label='Unit 1')\n",
    "plt.hist(por2, facecolor='blue',bins=np.linspace(0.0,0.4,20),alpha=0.2,density=False,edgecolor='black',label = 'Unit 2')\n",
    "plt.xlim([0.0,0.4]); plt.ylim([0,8.0])\n",
    "plt.xlabel('Porosity (fraction)'); plt.ylabel('Frequency'); plt.title('Porosity Unit 1 and 2')\n",
    "plt.legend(loc='upper left')\n",
    "plt.grid(True)\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.hist(por1, facecolor='red',bins=np.linspace(0.0,0.4,50),histtype=\"stepfilled\",alpha=0.2,density=True,cumulative=True,edgecolor='black',label='Unit 1')\n",
    "plt.hist(por2, facecolor='blue',bins=np.linspace(0.0,0.4,50),histtype=\"stepfilled\",alpha=0.2,density=True,cumulative=True,edgecolor='black',label='Unit 2')\n",
    "plt.xlim([0.0,0.4]); plt.ylim([0,1.0])\n",
    "plt.xlabel('Porosity (fraction)'); plt.ylabel('Cumulative Probability'); plt.title('Porosity Unit 1 and 2')\n",
    "plt.legend(loc='upper left')\n",
    "plt.grid(True)\n",
    "\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.2, hspace=0.3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The histogram and cumulative distribution functions indicate the porosity distributions in both wells 1 and 2. Ocular inspection does indicate: \n",
    "\n",
    "1. not a lot of samples are available for each (n=20) \n",
    "2. the distributions look different\n",
    "\n",
    "Let's use statistics to learn more about the subsurface from these two wells."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Confidence Intervals\n",
    "\n",
    "Let's first demonstrate the calculation of the confidence interval for the sample mean at a 95% confidence level.  This could be interpreted as the interval over which there is a 95% confidence that it contains the true population mean.  We use the student's t distribution as we assume we do not know the variance and the sample size is small. \n",
    "\n",
    "\\begin{equation}\n",
    "x̅ \\pm t_{\\frac{\\alpha}{2},n-1} \\times \\frac {s}{\\sqrt{n}} \n",
    "\\end{equation}\n",
    "\n",
    "**Coding Tip** it is easy to combine text, variables and equations in printed output.\n",
    "\n",
    "* we just need to use the '+' to concatenate (add) the values and the 'str()' command to convert the values to strings.\n",
    "\n",
    "For example:\n",
    "\n",
    "```python \n",
    "print('The average porosity is ' + str(np.average(por1)) + '.')\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The confidence interval for the Well 1 mean porosity is 0.164 +/- 0.011, with a range of 0.154, 0.175\n",
      "The confidence interval for the Well 2 mean porosity is 0.2 +/- 0.018, with a range of 0.182, 0.218\n"
     ]
    }
   ],
   "source": [
    "ci_95_por1 = st.t.interval(alpha = 0.9, df = len(df)-1, loc=np.mean(por1), scale=st.sem(por1))\n",
    "print('The confidence interval for the Well 1 mean porosity is ' + str(round(np.mean(por1),3)) + ' +/- ' + str(round(ci_95_por1[1]-np.mean(por1),3))\n",
    "      + ', with a range of ' + str(round(ci_95_por1[0],3)) + ', ' + str(round(ci_95_por1[1],3)))\n",
    "\n",
    "ci_95_por2 = st.t.interval(alpha = 0.9, df = len(df)-1, loc=np.mean(por2), scale=st.sem(por2))\n",
    "print('The confidence interval for the Well 2 mean porosity is ' + str(round(np.mean(por2),3)) + ' +/- ' + str(round(ci_95_por2[1]-np.mean(por2),3))\n",
    "      + ', with a range of ' + str(round(ci_95_por2[0],3)) + ', ' + str(round(ci_95_por2[1],3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So with one line of code we got the confidence interval in the mean with the Student's distribution!\n",
    "\n",
    "```python\n",
    "ci_95_por1 = st.t.interval(alpha = 0.9, df = len(df)-1, loc=np.mean(por1), scale=st.sem(por1))\n",
    "```\n",
    "where:\n",
    "\n",
    "* alpha is the significance level (alpha should be 1 - significance)\n",
    "\n",
    "* df is the degrees of freedom\n",
    "\n",
    "* loc is the sample mean calculated with the numpy.mean of an array function\n",
    "\n",
    "* scale is the stand error of the mean calculated with the scipy.stats.sem function\n",
    "\n",
    "The output is an array with:\n",
    "\n",
    "```python\n",
    "ci_95_por2[0]                            # lower bound\n",
    "ci_95_por2[1]                            # upper bound\n",
    "```\n",
    "\n",
    "\n",
    "#### Using Confidence Intervals to Build an Uncertainty Model for by-Well Distributions\n",
    "\n",
    "The confidence interval for the mean could be applied directly as an uncertainty model! The uncertainty in mean porosity at each well could be applied to calculate the uncertainty in the mean porosity applied to model porosity around the two wells. This method would account for the limited number of porosity samples available. If more porosity samples we available then this uncertainty would decrease (note the $\\sqrt{n}$ in the denominator of the confidence interval in the mean calculation above.\n",
    "\n",
    "One possible work flow would be to apply the affine correction from the *GeostatsPy* package to formulate a the P5 and P95 porosity distributions for porosity in each well. Note, these distribution scenarios are based only on uncertainty in the mean porosity and assumes the distribution tails are not anchored (are free to shift)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA54AAAGWCAYAAAAQQ1yjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3X+cHXV56PHPQwIETFgUIgbDj1C4uQECCIs/qoRYtEY0WFus4I9GidAiAiqrRen1FxRbY6+iUFGKV9QigrcK9KJiiyt6bYUNRigGLgENJGwACUiiATbkuX+cSXKy7Nlzsruze3583q/XvvbMzPc783zP2Z3nPHNm5kRmIkmSJElSWXaY6AAkSZIkSe3NwlOSJEmSVCoLT0mSJElSqSw8JUmSJEmlsvCUJEmSJJXKwlOSJEmSVCoLT5UmIu6MiPnjuL0PR8Q/jeH61kfEAcXjr0TEBWO47ksj4n+M1fraQUTsWzznkyY6lrEWEb+OiFcVjz8WEV+f6JgktS7z67DrNr8OYn5Vs7DwbFPFP+KGYkfzUET8r4iYOp4xZOYhmdlbxDOqnUFE9EbEkxGxLiKeiIilEXFuROxctb0LM/NdDa6rbrvMnJqZ94005qrtvSMifjJo3X+VmeePdt1DbOtZCTwi9o+IjIjJY729OrE8a9zDycz7i+f8mQbWXeqYIuL7EfHBqukXFtsbat4LxnC7O0XEt4r/3xzPN5aSGmN+rbsu82vJzK8j2u5LI+IHEbE2Ih6JiGsiYsZYrV+NsfBsbwszcypwJHA08Dfbu4Lx3pnW8Z7MnAbMAM4BTgJuiIgYy4002ZhbUhs8hzcDx1ZNzwPuGmLePZm5Zoy3/RPgbcBYr1fS2DG/jkCTjbkltcFzOFH59bnAl4D9gf2AdcD/GsP1qwEWnh0gM1cD3wUOBYiIvSPiuuKoz4qIOHVz2+LI6bci4usR8QTwjojYOSI+GxEPFj+f3XwkNCL2jIh/jYjHi/X9OCJ2KJb9OiJeFRELgA8Dby6OEP8iIt4UEUur44yIcyLiOw2M53fFkd4TgJcBr6uK/evF4ynFGB4tYrs1IvaKiL8FjgEuLmK5uGifEXFGRNwD3FM178CqTe9ZHC1bFxE/ioj9inbPOjq4+ahvRMwBLgVeVmzv8WL5NkdOI+LU4rVYW7w2e1cty4j4q4i4JyIei4hLRvNmoHhdeiLi9oj4bUR8MyKmVC1/Q0QsK45831u8fkREV0RcHhH9EbE6Ii6I4rSd4ujr/42Iz0TEWuCbNcb9uoj4ebHuByLiY1Xb3eZ5LJ7D84v1rouIGyNiz6L5zcXvx4v1H1s8d3Or1vf8qHwqMX0ET9PNwMs3/y1T+Zv5LNA9aN7mOIiI1xfP2+MR8dOIOGx7N5qZT2fmZzPzJ0DdI9OSJpb51fxazfzakInKr9/NzGsy84nM/D1wMfDyEcSvUbDw7AARsQ9wPPDzYtY3gFXA3sCJwIURcVxVlzcA3wJ2B/4ZOA94KXAEcDjwYrYe3T2nWNd0YC8qCTCrt5+Z3wMuBL5ZnOpxOHAdMKtIHJu9Dfhao+PKzPuBPio7qMEWAV3APsAewF8BGzLzPODHVI7uTs3M91T1+RPgJcDBNTb5VuB8YE9gGZXnpl6My4tt/0exvd0Ht4mIPwI+Cfw5laPNK4GrBjV7PZWj6ocX7V5Tb9t1/DmwAJgFHAa8o4jlxcBXgQ9Qef3nAb8u+lwBbAQOBF4E/DFQfUrVS4D7gOdTeS2HGvfvgL8o1v064PSI+JNh4nwL8M5inTsBPcX8ecXv3Yv1/4jKc/a2qr4nA/+WmY/UfTae7RZgZyrP9+bt/QBYMWjezQARcSTwZeAvqfy9fRG4LqpOVZPUfsyv5tchmF+H1yz5dR5w5yjXoe1k4dnevlMcCfsJ8CMqCXAf4BXAX2fmk5m5DPgn4O1V/f4jM7+TmZsycwOVhPCJzHy42Ml8vKr9AJWd+X6ZOZCZP87MbRLjUDLzKSpH7d4GEBGHUDn94V+3c4wPAs8bYv4AlR3UgZn5TGYuzcwn6qzrk5m5thjzUP5PZt5cxH4elaON+2xnvEN5K/DlzLytWPeHinXvX9Xm7zLz8eLNwA+pvEkZjc9l5oOZuRa4vmp9i4tYflC8/qsz866I2At4LfDe4oj4w8BnqJyOtdmDmfn5zNxY6znMzN7MvKNY9+1U3qQdO1Tbwv/KzP9XrO9qhh/3FcBbqo6Yvp3teKM1KM6ngJ8B8yLieVQS8H1U3lRtnncwlf8rgFOBL2bmz4q/tyuAp6i8oZTUfsyv5tdazK/DaIb8Wnxi+hEqBwE0jiw829ufZObumblfZr672LnsDazNzHVV7VYCL6yafmDQevYu2lS333yqyhIqR6lujIj7IuLc7Yhv844sqOzEri52SNvjhcDaIeZ/Dfg+cFVUTl/6VETsWGddg8ddc3lmri+2u3ft5g3b5vkt1v0o274m1dc5/B6odSOLjcDgce4IbCp+6q1vH+DeIda7X7Ge/uJUl8epHHV8flWbes8fEfGSiPhhVC7s/y2Vo7Z7DtOl0XGTmT+jcsT32Ij471SOHF9XI447i1OI1kfEUEf0oXK0dR6VI/6bb+Lwk6p5D2Tm5tdtP+Cczc9N8fzsw9j8fUhqPuZX8yuYX4eKo6nza1RO8f4ucHZm/ngk69DIWXh2ngeB50XEtKp5+wKrq6YHH1F9kMo/fnX7BwEyc11mnpOZBwALgfcPOq2o1jrJzP8Enqayk3kL23n0rDgaehSVo2SD1z2QmR/PzIOBP6RyKs1f1IqlzvzNthx9jcodDJ9H5Xn4XTF716q21Xdiq7febZ7fiHgOlaPJq2v2qO1+Kke2q82ishPf9Ozmz/IA8Ac15j8F7Fm82do9M3fLzEOq2gwe51DjvpJKstonM7uoXKcykutpaj2nV1A5yv924FuZ+eSQnSt3hJxa/NRKPDdT+ducx9a/sf9L5ZqQLacBFR4A/rbqudk9M3fNzG9s37AktTDzq/l1OObXrSYkv0bl2uF/A87PzBF9YqvRsfDsMJn5APBT4JNRuUHAYVRO/xjueopvAH8TEdOLi88/Amy+ycDrI+LA4qjqE1RuiDLUTVEeAvavOk1js69SucB7Y1ZuqFJXROwaEccC11K5VuCGIdq8MiLmRuXi/CeonBq0Oa6HgAMa2dYgx0fEKyJiJyrXovwsMx8oTo9aDbwtIiZFxClsm1weAmYW/YZyJfDOiDiiuGbhwmLdvx5BjP8beF1E/HERy95UrhcafE1LLZcXsRwXETtE5Zbm/z0z+4EbgX+IiN2KZX9QvA61DDXuaVQ+EXiyuN7lLds/RAAeoXKEefDr+DXgjVSS41dHuO7NfkrlWpm3USTGzHys2Pbb2DYxXgb8VXHEOSLiOVG50cO0wSutJyo3G9l8M4qdiv/TMb2zpKSxZ37dEov5dWjm163GPb9GxAuBm4BLMvPSUcavEbLw7EwnUzlq9yDwbeCjmfmDYdpfQOUmA7cDdwC3FfMADqJy9Gg98B/AP2bx3WKDXFP8fjQibqua/zUqdwNs5MjTxRGxjsoO97NUksCCGkcaX0DlBg5PAMupXCuw+XvOLgJOjMod7D7XwHY3uxL4KJVTgI6icu3IZqdSuVbgUeAQKjvVzW6icgH7moj4zeCVZua/A/+jGE8/laR60uB2jcjMO6m8vp8s4vwPKtdSfLzB/rdQudnAZ4DfUnneNh8t/gsqNyD4JfAYled3uO/AGmrc7wY+UbyOH6FyXcl2y8od6f4W+L/FqTcvLeavovL3mQxxpH4E21hK5SYI/1W16MdUToG6uaptH5W/gYupPDcrKG4oMQJ3AxuonAr2/eLxfsP2kNQszK/m11r9za/bbmO88+u7qBTTH42tpwKvH9EANGKR9a9Tl0oTEbsADwNHZuY9Ex2PWl9EfJnKjRi2+3v1JKldmF811syvGq1W/xJatb7TgVtNihoLUblT4Z9SuR29JHUy86vGjPlVY6G0U20j4ssR8XBE/FeN5RERn4vKl/reHpXv6VEHiYhfA2dT+a4yaVQi4nwqp+wsycxfTXQ8UlnMr6rH/KqxZH7VWCntVNuImEfluoSvZuahQyw/HjiTyhcvvwS4KDNfUkowkiS1CfOrJKkVlfaJZ2bezNDf/7TZG6gkzSxu+717RAx3IbUkSR3P/CpJakUTeY3nC9n2C3FXFfP6BzeMiNOA0wCmTJly1L777jsuAZZl06ZN7LBDa99Q2DFMjKeeepLqsxQmT57Mxo0bt0yP+gyGCbrX2OQdd2TjwMDEbHyMOIbm8MCqVb/JzOkTHccEayi/mlubUzuMwzE0h/Ecw+D3J9urVt/Jk3dk48YG8lIT3yvV3LrVRBaeQ30v3ZB/Npn5JeBLALNnz8677767zLhK19vby/z58yc6jFFxDBOju3s/+vqev2W6t3cx8+dfvnX53DvoWzLyf+vuRRvo6xr/r4zsPetTzP/cB8d9u2PJMTSHgJUTHUMTaCi/mlubUzuMwzE0h/Ecw+D3J9vdv8b7l96nL2D+TvVvojtR718aYW7daiIP5awC9qmanknle68kSdLImV8lSU1nIgvP64C/KO6+91Lgt5n5rNNsJUnSdjG/SpKaTmmn2kbEN4D5wJ4RsQr4KLAjQGZeCtxA5Y57K4DfA+8sKxZJktqF+VWS1IpKKzwz8+Q6yxM4Yyy2NTAwwKpVq3jyySfHYnWl6+rqYvny5RMdxjamTJnCzJkz2XHHHSc6FEnSMMYrv5pbx4b5VZIqJvLmQmNm1apVTJs2jf3335+I5rywuNq6deuYNm3aRIexRWby6KOPsmrVKmbNmjXR4UiSmoC5dfTMr5K0VWvfJ7rw5JNPsscee7REYmxGEcEee+zRMke1JUnlM7eOnvlVkrZqi8ITMDGOks+fJGkwc8Po+RxKUkXbFJ6SJEmSpObUFtd4DrZw4Tz6+8fuO8RnzNiP66+/edg2kyZNYu7cuWzcuJE5c+ZwxRVXsOuuu/K9732Ps88+m2eeeYZ3vetdnHvuuQC84x3v4Ec/+hFdXV0AfOUrX+GII44Ys5glSRpL5lZJ0mi0ZeHZ37+Svr7pY7a+7u76iXaXXXZh2bJlALz1rW/l0ksv5eyzz+aMM87gBz/4ATNnzuToo4/mhBNOYJ99Kt/rvWTJEk488cQxi1OSpLKYWyVJo+GptiU45phjWLFiBbfccgsHHnggBxxwADvttBMnnXQS11577USHJ0lSyzG3SlJrs/AcYxs3buS73/0uc+fOZfXq1VuOwALMnDmT1atXb5k+77zzOOyww3jf+97HU089NRHhSpLU9MytktT6LDzHyIYNGzjiiCPo7u5m3333ZfHixVS+w3tbm+9u98lPfpK77rqLW2+9lbVr1/L3f//34x2yJElNzdwqSe2jLa/xnAjV16FsNnPmTB544IEt06tWrWLvvfcGYMaMGQDsvPPOvPOd7+TTn/70+AUrSVILMLdKUvvwE88SHX300dxzzz386le/4umnn+aqq67ihBNOAKC/vx+AzOQ73/kOhx566ESGKklSSzC3SlJrastPPGfM2K+hu+Vtz/pGYvLkyVx88cW85jWv4ZlnnuGUU07hkEMOYd26dbz1rW/lkUceITM54ogjuPTSS8csXrWuel9XsHbNQ3TPXbNlevHpG+g5846ty+8foHvRwIi3PyMT8MvOJT2buVWaOCP5OqPFi3vo6VkEwJr7HyI3bSwjNAA2PvEML9zp/hH3n7QDQ75/WfzhTfRcuKFuf9+/tIa2LDzrfS9YGdavXz/k/OOPP57jjz/+WfNvuummskNSC6r3dQXdc9fQt2Trv23v07HNdPeiAfq6RrPjdactaWjmVmnijOTrjHp7J2/pU3n/sHMZoVXWv2jDKN9/DK13Eg2u1/cvrcBTbSVJkiRJpbLwlCRJkiSVysJTkiRJklQqC09JkiRJUqksPCVJkiRJpbLwlCRJkiSVqi2/TmXhgnn0rx7D7xp74X5c/73hbyM/adIk5s6du2X6pJNO4txzzx2zGIby+OOPc+WVV/Lud797u/p97GMfY+rUqfT09JQUmSSp3Zhbh2dulaThtWXh2b96JX1Ltu+7jobT/YH6iXaXXXZh2bJlY7bNRjz++OP84z/+43YnR0mStpe5VZI0Gp5qW6Lf/va3zJ49m7vvvhuAk08+mcsuuwyAqVOncs4553DkkUdy3HHH8cgjjwBw7733smDBAo466iiOOeYY7rrrLgAeeugh3vjGN3L44Ydz+OGH89Of/pRzzz2Xe++9lyOOOIIPfOADACxZsoSjjz6aww47jI9+9KNbYvnbv/1bZs+ezate9aot8UiS1GrMrZLUmiw8x8iGDRs44ogjtvx885vfpKuri4svvph3vOMdXHXVVTz22GOceuqpAPzud7/jyCOP5LbbbuPYY4/l4x//OACnnXYan//851m6dCmf/vSntxxxPeusszj22GP5xS9+wW233cYhhxzC3/3d3/EHf/AHLFu2jCVLlnDjjTdyzz33cMstt7Bs2TKWLl3KzTffzNKlS7nqqqv4+c9/zr/8y79w6623TtjzJElSo8ytktQ+2vJU24lQ63SgV7/61VxzzTWcccYZ/OIXv9gyf4cdduDNb34zAG9729v40z/9U9avX89Pf/pT3vSmN21p99RTTwFw00038dWvfhWoXPPS1dXFY489ts22brzxRm688UZe9KIXAbB+/Xruuece1q1bxxvf+EZ23XVXAE444YQxHLkkSeUwt0pS+7DwLNmmTZtYvnw5u+yyC2vXrmXmzJlDtosINm3axO677z7i61kykw996EP85V/+5TbzP/vZzxIRI1qnJEnNxtwqSa3HU21L9pnPfIY5c+bwjW98g1NOOYWBgQGgkjS/9a1vAXDllVfyile8gt12241Zs2ZxzTXXAJVkt/lI7nHHHccXvvAFAJ555hmeeOIJpk2bxrp167Zs6zWveQ1f/vKXWb9+PQCrV6/m4YcfZt68eXz7299mw4YNrFu3juuvv37cxi9J0lgzt0pS62nLTzxnvHC/hu6Wtz3rq2fzdSibLViwgFNOOYV/+qd/4pZbbmHatGnMmzePCy64gJ6eHp7znOdw5513ctRRR9HV1cU3v/lNAP75n/+Z008/nQsuuICBgQFOOukkDj/8cC666CJOO+00Lr/8ciZNmsQXvvAFXvayl/Hyl7+cQw89lNe+9rUsWbKE5cuX87KXvQyo3GTh61//OkceeSRvfvObOeKII9hvv/045phjxuy5kSR1BnOruVWSRqMtC8963wtWhmeeeWbI+cuXL9/y+H/+z/8JsOVI6vnnn8/555+/TftZs2bxve9971nr2Wuvvbj22mufNf/KK6/cZvrss8/m7LPPfla78847j/POO6/OKCRJGpq51dwqSaPhqbaSJEmSpFK15SeerWDztSJqLwsXzqO/f+Snoq1d8xDdc9fUXn7/AN2LBrZML/7wJnou3LBlekYm4M0uJHUmc6taVb33DzNm1Om/4C76Vw9sM2/x6RvoOfOOSv+uHL5/z5P0PzJ8m+H4/kONsPCUxlB//0r6+qaPuH/33DX0Lan9b9m9aIC+rq079t5JbDPtTl+SpNYz2vcP/asHnvX+offpGPY9xTb9H8lB7ye2l+8/VJ+n2kqSJEmSSmXhKUmSJEkqlYWnJEmSJKlUbXmN58J58+hfOYbfNbbfflx/8/C3kZ80aRJz585l48aNzJkzhyuuuIJdd92Viy66iMsuu4zM5NRTT+W9730vAB/72Me47LLLmD69cj7/hRdeyPHHHz9mMUuSNJbMrZKk0WjLwrN/5Ur6po/8Au3BuhtItLvssgvLli0D4K1vfSuXXnopf/zHf8xll13GLbfcwk477cSCBQt43etexwte8AIA3ve+99HT0zNmcUqSVBZzqyRpNDzVtgTHHHMMK1asYPny5bz0pS9l1113ZfLkyRx77LF8+9vfnujwJElqOeZWSWptFp5jbOPGjXz3u99l7ty5HHroodx88808+uij/P73v+eGG27ggQce2NL24osv5rDDDuOUU07hsccem8CoJUlqXuZWSWp9Fp5jZMOGDRxxxBF0d3ez7777snjxYubMmcNf//Vf8+pXv5oFCxZw+OGHM3ly5ezm008/nXvvvZdly5YxY8YMzjnnnAkegSRJzcXcKkntoy2v8ZwI1dehVFu8eDGLFy8G4MMf/jAzZ84EYK+99trS5tRTT+X1r3/9+AQqSVKLMLdKUvvwE8+SPfzwwwDcf//9/Mu//Asnn3wyAP39/VvafPvb3+bQQw+dkPgkSWo15lZJaj1t+YnnjP32a+hueduzvpH6sz/7Mx599FF23HFHLrnkEp773Oeybt06PvjBD7Js2TIigv33358vfvGLYxavJEljzdwqSRqNtiw8630vWBnWr18/5Pwf//jHQ87/2te+VmY4kiSNKXOrJGk0PNVWkiRJklQqC09JkiRJUqnapvDMzIkOoaX5/EmSBjM3jJ7PoSRVtMU1nlOmTOHRRx9ljz32ICImOpyWk5k8+uijTJkyZaJDkSQ1CXPr6JlfO8fChfPo7x/5zbfWrnmI7rlrRt7//gG6Fw1sM2/xhzfRc+GGhvrPyAT8P1e52qLwnDlzJqtWreKRRx6Z6FAa8uSTTzZdEpoyZcqW70GTJMncOjbMr52hv38lfX3TR9y/e+4a+paM/G1596IB+rq2LRx7J/GsebVZdKp8bVF47rjjjsyaNWuiw2hYb28vL3rRiyY6DEmSajK3SpLGUttc4ylJkiRJak4WnpIkSZKkUll4SpIkSZJKZeEpSZIkSSqVhackSZIkqVQWnpIkSZKkUll4SpIkSZJKZeEpSZIkSSqVhackSZIkqVQWnpIkSZKkUll4SpIkSZJKVWrhGRELIuLuiFgREecOsXzfiPhhRPw8Im6PiOPLjEeSpFZnbpUktaLSCs+ImARcArwWOBg4OSIOHtTsb4CrM/NFwEnAP5YVjyRJrc7cKklqVWV+4vliYEVm3peZTwNXAW8Y1CaB3YrHXcCDJcYjSVKrM7dKklpSZGY5K444EViQme8qpt8OvCQz31PVZgZwI/Bc4DnAqzJz6RDrOg04DWD69OlHXX311aXEPF7Wr1/P1KlTJzqMUXEMQ1u+/A7mzJk88v6/3MCcmVF7+a82MWfS1un1e81k6kOrRry9ZtEO43AMzeGVZ/YszczuiY6jLObW2tohL0F7jKMTx1B2/q/bf9D7A2iPfbpjaA5jlVvLLDzfBLxmUHJ8cWaeWdXm/UUM/xARLwMuBw7NzE211jt79uy8++67S4l5vPT29jJ//vyJDmNUHMPQurv3o69v+sj7z72DviW1E1f3og30dW1NTL1nfYr5n/vgiLfXLNphHI6hOcQ9m9q98DS31tAOeQnaYxydOIay83/d/oPeH0B77NMdQ3MYq9xa5qm2q4B9qqZn8uzTfRYDVwNk5n8AU4A9S4xJkqRWZm6VJLWkMgvPW4GDImJWROxE5QYH1w1qcz9wHEBEzKGSHB8pMSZJklqZuVWS1JJKKzwzcyPwHuD7wHIqd9i7MyI+EREnFM3OAU6NiF8A3wDekWWd+ytJUoszt0qSWtXITyZvQGbeANwwaN5Hqh7/Enh5mTFIktROzK2SpFZU5qm2kiRJkiRZeEqSJEmSymXhKUmSJEkqlYWnJEmSJKlUFp6SJEmSpFJZeEqSJEmSSmXhKUmSJEkqlYWnJEmSJKlUFp6SJEmSpFJZeEqSJEmSSmXhKUmSJEkqlYWnJEmSJKlUFp6SJEmSpFJZeEqSJEmSSmXhKUmSJEkqlYWnJEmSJKlUFp6SJEmSpFJZeEqSJEmSSmXhKUmSJEkqlYWnJEmSJKlUFp6SJEmSpFJZeEqSJEmSSmXhKUmSJEkqlYWnJEmSJKlUFp6SJEmSpFJZeEqSJEmSSmXhKUmSJEkqlYWnJEmSJKlUFp6SJEmSpFJNnugApFaycOE8+vtX1lw+Y0ad/gvuon/1QM3la+8foHtR7eUzMoGoF6YkSdoOg/P74sU99PQsarj/2jUP0T13zYi3P6Mrh4+v50n6H6ndxvcHagUWntJ26O9fSV/f9JH3Xz1A35La/3bdiwbo6xoucZhUJEkaa4Pze2/v5O3K991z1wyb30er/5H0/YFanqfaSpIkSZJKZeEpSZIkSSqVhackSZIkqVQWnpIkSZKkUll4SpIkSZJKZeEpSZIkSSqVhackSZIkqVQWnpIkSZKkUll4SpIkSZJKZeEpSZIkSSqVhackSZIkqVQWnpIkSZKkUll4SpIkSZJKVbfwjIi+iDgjIp47HgFJktTuzK2SpE7TyCeeJwF7A7dGxFUR8ZqIiJLjkiSpnZlbJUkdZXK9Bpm5AjgvIv4H8Hrgy8CmiPgycFFmri05RmncLFw4j/7+lTWXr13zEN1z14x4/WvvH6B70UDN5TMyAd97Su3O3CqNr3r5fcaMOv0X3EX/6mHyd1cO37/nSfofGb7NcHx/oHZQt/AEiIjDgHcCxwP/G/hn4BXATcARpUUnjbP+/pX09U2vubx77hr6ljT0bzN0/0UD9HUNlzhMKlKnMLdK46defq/bf/XAqPJ//yNZJ//X4/sDtb66/0ERsRR4HLgcODcznyoW/SwiXl5mcJIktSNzqySp0zRy6OZNmXlf9YyImJWZv8rMPy0pLkmS2pm5VZLUURq5udC3GpwnSZIaY26VJHWUmp94RsR/Bw4BuiKi+ujrbsCUsgOTJKndmFslSZ1quFNtZ1O5097uwMKq+euAU8sMSpKkNmVulSR1pJqFZ2ZeC1wbES/LzP8Yx5gkSWpL5lZJUqca7lTbD2bmp4C3RMTJg5dn5lmlRiZJUpsxt0qSOtVwp9ouL373jUcgkiR1AHOrJKkjDXeq7fXF7yvGLxxJktqXuVWS1KmGO9X2eiBrLc/ME0qJSJKkNmVulSR1quFOtf30uEUhSVJnMLdKkjrScKfa/mg8A5Ekqd2ZWyVJnWqHWgsi4uri9x0RcXvVzx0RcXsjK4+IBRFxd0SsiIhza7T584j4ZUTcGRFXjmwYkiQ1P3OrJKlTDXeq7dkZoUAvAAAgAElEQVTF79ePZMURMQm4BHg1sAq4NSKuy8xfVrU5CPgQ8PLMfCwinj+SbUmS1CLMrZKkjlTzE8/M7C9+rwSeAg4HDgOeKubV82JgRWbel5lPA1cBbxjU5lTgksx8rNjWw9s/BEmSWoO5VZLUqSKz5s31Kg0i3gV8BLgJCOBY4BOZ+eU6/U4EFmTmu4rptwMvycz3VLX5DvD/gJcDk4CPZeb3hljXacBpANOnTz/q6quvbniAzWj9+vVMnTp1osMYlXYdw/LldzBnTu0TAZb/cgNzZsaIt7n8V5uYM2nE3Z9l/V4zmfrQqrFb4QRph3E4hubwyjN7lmZm90THUY+5dey1Q16C9hhHM46hXn4fbP36PZk69Tdb+zdZ/m9EO+zTHUNzGKvc2kjheTfwh5n5aDG9B/DTzJxdp9+bgNcMSo4vzswzq9r8KzAA/DkwE/gxcGhmPl5rvbNnz8677767kbE1rd7eXubPnz/RYYxKu46hu3s/+vqm1+zTPfcO+pY0nrie1X/RBvq6Rp64Bus961PM/9wHx2x9E6UdxuEYmkPcs6lVCk9z6xhrh7wE7TGOZhxDvfw+WG/vYubPv3xr/ybL/41oh326Y2gOY5Vba55qW2UVsK5qeh3wQIP99qmangk8OESbazNzIDN/BdwNHNTAuiVJamXmVklSR6l56CYi3l88XA38LCKupfKl128Abmlg3bcCB0XErGIdJwFvGdTmO8DJwFciYk/gvwH3bdcIJElqEeZWSVKnGu6cgWnF73uLn82ubWTFmbkxIt4DfJ/KNSZfzsw7I+ITQF9mXlcs++OI+CXwDPCBzacdSZLUhsytkqSOVLPwzMyPj3blmXkDcMOgeR+pepzA+4sfSZLamrlVktSp6l4lHRHTgQ8ChwBTNs/PzD8qMS5JktqWuVWS1GkaubnQPwN3AbOAjwO/pnKNiSRJGhlzqySpozRSeO6RmZcDA5n5o8w8BXhpyXFJktTOzK2SpI7SyBcSDRS/+yPidVRu2z6zvJAkSWp75lZJUkdppPC8ICK6gHOAzwO7Ae8rNSpJktqbuVWS1FHqFp6Z+a/Fw98Cryw3HEmS2p+5VZLUaepe4xkRB0TE9RHxm4h4OCKujYgDxiM4SZLakblVktRpGrm50JXA1cALgL2Ba4BvlBmUJEltztwqSeoojRSekZlfy8yNxc/XgSw7MEmS2pi5VZLUUWpe4xkRzyse/jAizgWuopIU3wz8n3GITZKktmJulSR1quFuLrSUSjKMYvovq5YlcH5ZQUmS1KbMrZKkjlSz8MzMWeMZiCRJ7c7cKknqVHW/TiUidgROB+YVs3qBL2bmQM1OkiSpJnOrNLYWLpxHf//KmstnzKjTf8Fd9K/e+u+3+PQN9Jx5x9b+XcNfgr2w50n6H6ndZkZWn+ggdaa6hSfwBWBH4B+L6bcX895VVlCSJLU5c6s0hvr7V9LXN33k/VcP0Ldk69vi3qdjm+m6/R9J+rqGKywtOqVG/qOOzszDq6ZviohflBWQJEkdwNwqSeoojXydyjMR8QebJ4ovuH6mvJAkSWp75lZJUkdp5BPPD1C57ft9VM4T2A94Z6lRSZLU3sytkqSOMmzhGRE7ABuAg4DZVJLjXZn51DjEJklS2zG3SpI60bCFZ2Zuioh/yMyXAbePU0ySJLUtc6skqRM1co3njRHxZxHh7bgkSRob5lZJUkdp5BrP9wPPATZGxJNUTgnKzNyt1MgkSWpf5lZJUkepW3hm5rTxCESSpE5hbpUkdZqap9pGxPMj4rMR8a8RcWFEeBRWkqRRMLdKkjrVcNd4fhX4HfB5YBrwuXGJSJKk9mVulSR1pOFOtX1BZp5XPP5+RNw2HgFJktTGzK2SpI40XOEZEfFcKjc8AJhUPZ2Za8sOTpKkNmNulSR1pOEKzy5gKVuTI8DmI7MJHFBWUJIktSlzqySpI9UsPDNz/3GMQ5KktmdulSR1quFuLiRJkiRJ0qhZeEqSJEmSSmXhKUmSJEkqVUOFZ0S8IiLeWTyeHhGzyg1LkqT2Zm6VJHWSuoVnRHwU+GvgQ8WsHYGvlxmUJEntzNwqSeo0jXzi+UbgBOB3AJn5IDCtzKAkSWpz5lZJUkdppPB8OjOTyveLERHPKTckSZLanrlVktRRGik8r46ILwK7R8SpwL8Bl5UbliRJbc3cKknqKJPrNcjMT0fEq4EngNnARzLzB6VHJklSmzK3SpI6Td3CMyLeB1xjQpQkaWyYWyVJnaaRU213A74fET+OiDMiYq+yg5Ikqc2ZWyVJHaVu4ZmZH8/MQ4AzgL2BH0XEv5UemSRJbcrcKknqNI184rnZw8Aa4FHg+eWEI0lSRzG3SpI6Qt3CMyJOj4he4N+BPYFTM/OwsgOTJKldmVslSZ2m7s2FgP2A92bmsrKDkSSpQ5hbJUkdpWbhGRG7ZeYTwKeK6edVL8/MtSXHJklSWzG3SpI61XCfeF4JvB5YCiQQVcsSOKDEuCRJakfmVklSR6pZeGbm64vfs8YvHEmS2pe5VZLUqepe4xkR/56Zx9WbJ7WChQvn0d+/EoDFi3vo6Vm0zfK1ax6ie+6amv3X3j9A96KBEW9/Rg7+gENSJzK3SpI6zXDXeE4BdgX2jIjnsvXd8m5UvnNMajn9/Svp65sOQG/v5C2PN+ueu4a+JbWPx3QvGqCvazSFo0Wn1MnMrZKkTjXcJ55/CbyXSiJcytbk+ARwSclxSZLUjsytkqSONNw1nhcBF0XEmZn5+XGMSZKktmRulSR1qrrXeGbm5yPiUOBgYErV/K+WGZgkSe3K3CpJ6jSN3Fzoo8B8KsnxBuC1wE8Ak6MkSSNgbpUkdZodGmhzInAcsCYz3wkcDuxcalSSJLU3c6skqaM0UnhuyMxNwMaI2A14GL/gWpKk0TC3SpI6St1TbYG+iNgduIzKHfjWA7eUGpUkSe3N3CpJ6iiN3Fzo3cXDSyPie8BumXl7uWFJktS+zK2SpE5Ts/CMiCOHW5aZt5UTkiRJ7cncKknqVMN94vkPwyxL4I/GOBZJktqduVWS1JFqFp6Z+crxDESSpHZnbpUkdapGvsfzL4aa75dcS5I0MuZWSVKnaeSutkdXPZ5C5XvHbsMvuZYkaaTMrZKkjtLIXW3PrJ6OiC7ga42sPCIWABcBk4B/ysy/q9HuROAa4OjM7Gtk3ZIktSpzqySp0+wwgj6/Bw6q1ygiJgGXAK8FDgZOjoiDh2g3DTgL+NkIYpEkqR2YWyVJba2Razyvp3KnPagUqgcDVzew7hcDKzLzvmI9VwFvAH45qN35wKeAngZjliSppZlbJUmdJjJz+AYRx1ZNbgRWZuaquiuunOKzIDPfVUy/HXhJZr6nqs2LgL/JzD+LiF6gZ6jTgSLiNOA0gOnTpx919dWN5ObmtX79eqZOnTrRYYxKq45h+fI7mDOncrxl/fo9mTr1N9su/+UG5syM2v1/tYk5k0oNcbus32smUx+q++/Y9NphHI6hObzyzJ6lmdk90XHUY24de62alwZrh3GUMYYVK+5mYODpmst33BEOPLD25ykr7nmSgYHa73l3nAQHztia/9fnC5kaq7f2X7WJgYHa8e0IHNhE7w+gPfbpjqE5jFVubeQazx8BRMRum9tHxPMyc22drkO9e9/yHx8ROwCfAd7RQAxfAr4EMHv27Jw/f369Lk2tt7cXxzAxenoW0dc3HYDe3sXMn3/5tsvPvIO+JbX/LXou3EBfV+3CdLz1nvUp5n/ugxMdxqi1wzgcg7aHuXXstWpeGqwdxlHGGKrz94j618nvg/U+fQHzd/qbrf0/3lz5vxHtsE93DO2lkVNtT6Nyys4GYBOVpJfAAXW6rgL2qZqeCTxYNT0NOBTojQiAFwDXRcQJ3gRBktTOzK2SpE7TyKGfDwCHZOZv6rbc1q3AQRExC1gNnAS8ZfPCzPwtsOfm6eFOB5Ikqc2YWyVJHaWRu9reS+Vue9slMzcC7wG+DywHrs7MOyPiExFxwvauT5KkNmJulSR1lEY+8fwQ8NOI+Bnw1OaZmXlWvY6ZeQNww6B5H6nRdn4DsUiS1A7MrZKkjtJI4flF4CbgDirXoUiSpNExt0qSOkojhefGzHx/6ZFIktQ5zK2SpI7SyDWeP4yI0yJiRkQ8b/NP6ZFJktS+zK2SpI7SyCeem++W96GqeY3c8l2SJA3N3CpJ6ih1C8/MnDUegUiS1CnMrZKkTlO38IyIvxhqfmZ+dezDkSSp/ZlbJUmdppFTbY+uejwFOA64DTA5SpI0MuZWSVJHaeRU2zOrpyOiC/haaRFJktTmzK2SpE7TyF1tB/s9cNBYByJJUgczt0qS2loj13heT+VOe1ApVA8Gri4zKEmS2pm5VZLUaRq5xvPTVY83Aiszc1VJ8UiS1AnMrZKkjlKz8IyIA4G9MvNHg+YfExE7Z+a9pUcnbaeFC+fR37+y5vK1ax6ie+4aABafvoGeM+/Ydvn9A3QvGqjZf0YmEGMSq6TOY26VJHWq4T7x/Czw4SHmbyiWLSwlImkU+vtX0tc3veby7rlr6FtS+bPvfTq2PN6yfNEAfV3DFZYWnZJGxdwqSepIw91caP/MvH3wzMzsA/YvLSJJktqXuVWS1JGGKzynDLNsl7EORJKkDmBulSR1pOEKz1sj4tTBMyNiMbC0vJAkSWpb5lZJUkca7hrP9wLfjoi3sjUZdgM7AW8sOzBJktqQuVWS1JFqFp6Z+RDwhxHxSuDQYvb/ycybxiUySZLajLlVktSp6n6PZ2b+EPjhOMQiSVJHMLdKkjrNcNd4SpIkSZI0ahaekiRJkqRSWXhKkiRJkkpl4SlJkiRJKpWFpyRJkiSpVBaekiRJkqRSWXhKkiRJkkpl4SlJkiRJKpWFpyRJkiSpVJMnOgBJkiS1toUL59HfvxKAxYt76OlZtF3919z/ELlpY83lkyZB99w1tfvfN0AOZO3+O0D3ooGG41n84U30XLhhy/SMTCAa7i/p2Sw8JUmSNCr9/Svp65sOQG/v5C2PG9U9dw19S3Ye8fa7Fw3Q1zV2J/L1ToK+rupC06JTGi1PtZUkSZIklcrCU5IkSZJUKgtPSZIkSVKpLDwlSZIkSaWy8JQkSZIklcrCU5IkSZJUKgtPSZIkSVKpLDwlSZIkSaWy8JQkSZIklcrCU5IkSZJUKgtPSZIkSVKpLDwlSZIkSaWy8JQkSZIklcrCU5IkSZJUKgtPSZIkSVKpJk90AJIkSWpuCxfOo79/Zc3lM2bU6b/gLvpXD9Tu35XD9+95kv5HareZkQnE8EFImlAWnpIkSRpWf/9K+vqmj7z/6gH6loz8bWf/I0lf13CFpUWn1Ow81VaSJEmSVCoLT0mSJElSqSw8JUmSJEmlsvCUJEmSJJXKwlOSJEmSVCoLT0mSJElSqSw8JUmSJEmlsvCUJEmSJJXKwlOSJEmSVCoLT0mSJElSqSw8JUmSJEmlsvCUJEmSJJXKwlOSJEmSVKpSC8+IWBARd0fEiog4d4jl74+IX0bE7RHx7xGxX5nxSJLU6sytkqRWVFrhGRGTgEuA1wIHAydHxMGDmv0c6M7Mw4BvAZ8qKx5JklqduVWS1KrK/MTzxcCKzLwvM58GrgLeUN0gM3+Ymb8vJv8TmFliPJIktTpzqySpJUVmlrPiiBOBBZn5rmL67cBLMvM9NdpfDKzJzAuGWHYacBrA9OnTj7r66qtLiXm8rF+/nqlTp050GKNS1hhWrLibgYGnR9z/mYEBJk0eZvnTyaSoPN5zxkx+079qm+U7AgdOGvHmx936vWYy9aFV9Rs2uXYYh2NoDq88s2dpZnZPdBxlMbfW1g65FZp3HMuX38GcOcMk2Crr1+/J1Km/2bb/LzcwZ2aMfPu/2sSccczP7bA/dAzNoR3GMFa5tbE9yMgMtXcZssqNiLcB3cCxQy3PzC8BXwKYPXt2zp8/f4xCnBi9vb04hqH19Cyir2/6iPt3z72DviW1/6y7F22gr6vyp9l71qd40yUfHPG2mkHvWZ9i/udaewzQHuNwDBon5tYa2iG3QvOOY3vyc2/vYubPv3zb/mcOn5/rbv/Crfl7PLTD/tAxNId2GMNYKbPwXAXsUzU9E3hwcKOIeBVwHnBsZj5VYjySJLU6c6skqSWVeY3nrcBBETErInYCTgKuq24QES8CvgickJkPlxiLJEntwNwqSWpJpRWembkReA/wfWA5cHVm3hkRn4iIE4pmS4CpwDURsSwirquxOkmSOp65VZLUqso81ZbMvAG4YdC8j1Q9flWZ25ckqd2YWyVJrajMU20lSZIkSbLwlCRJkiSVy8JTkiRJklQqC09JkiRJUqksPCVJkiRJpbLwlCRJkiSVysJTkiRJklQqC09JkiRJUqksPCVJkiRJpbLwlCRJkiSVysJTkiRJklSqyRMdgDrLwoXz6O9fWXP52jUP0T13zYjXv/b+AboXDdRcPiMTiBGvX5KkdlQvP8+YUaf/grvoX13Jv4tP30DPmXds278rh+/f8yT9j9RuY/6WWp+Fp8ZVf/9K+vqm11zePXcNfUtG/mfZvWiAvq7hEpNJS5Kkwerl57r9Vw9syd+9T8d25/L+R9L8LbU5T7WVJEmSJJXKwlOSJEmSVCoLT0mSJElSqSw8JUmSJEmlsvCUJEmSJJXKwlOSJEmSVCoLT0mSJElSqSw8JUmSJEmlsvCUJEmSJJXKwlOSJEmSVCoLT0mSJElSqSw8JUmSJEmlsvCUJEmSJJXKwlOSJEmSVCoLT0mSJElSqSZPdABqLQsXzqO/fyWLF/fQ07PoWcvX3P8QuWljzf4bn3iGF+50f83lk3aA7kUDI45vRiYQI+4vSVIz2px/R6pefp40Cbrnrhnx+tfeP7Alfy/+8CZ6LtywXf3N31L7s/DUdunvX0lf33R6eyfT1zf9Wcu7566hb8nONft3L9pAX1eZicWkJUlqP5vz70jVy8+j1b1oYEt+753ECHK9+Vtqd55qK0mSJEkqlYWnJEmSJKlUFp6SJEmSpFJZeEqSJEmSSmXhKUmSJEkqlYWnJEmSJKlUFp6SJEmSpFJZeEqSJEmSSmXhKUmSJEkqlYWnJEmSJKlUFp6SJEmSpFJZeEqSJEmSSmXhKUmSJEkqlYWnJEmSJKlUFp6SJEmSpFJZeEqSJEmSSjV5ogOQJElqdwsXzqO/f2XddosX99DTs+hZ89eueYjuuWtGvP219w/QvWhgxP3rmZEJRGnrl9T6LDwlSZJK1t+/kr6+6XXb9fZOHrJd99w19C0Z+du27kUD9HWVWRhadEoanqfaSpIkSZJKZeEpSZIkSSqVhackSZIkqVQWnpIkSZKkUll4SpIkSZJKZeEpSZIkSSqVhackSZIkqVQWnpIkSZKkUll4SpIkSZJKZeEpSZIkSSqVhackSZIkqVQWnpIkSZKkUpVaeEbEgoi4OyJWRMS5QyzfOSK+WSz/WUTsX2Y8kiS1OnOrJKkVlVZ4RsQk4BLgtcDBwMkRcfCgZouBxzLzQOAzwN+XFY8kSa3O3CpJalVlfuL5YmBFZt6XmU8DVwFvGNTmDcAVxeNvAcdFRJQYkyRJrczcKklqSZGZ5aw44kRgQWa+q5h+O/CSzHxPVZv/KtqsKqbvLdr8ZtC6TgNOKyYPBf6rlKDHz57Ab+q2am6OoTm0wxigPcbhGJrD7MycNtFBlMXcOqx2+PuF9hiHY2gOjqE5tMMYxiS3Th6LSGoY6ujq4Cq3kTZk5peALwFERF9mdo8+vInjGJqDY2ge7TAOx9AcIqJvomMombm1hnYYA7THOBxDc3AMzaFdxjAW6ynzVNtVwD5V0zOBB2u1iYjJQBewtsSYJElqZeZWSVJLKrPwvBU4KCJmRcROwEnAdYPaXAcsKh6fCNyUZZ37K0lS6zO3SpJaUmmn2mbmxoh4D/B9YBLw5cy8MyI+AfRl5nXA5cDXImIFlaOxJzWw6i+VFfM4cgzNwTE0j3YYh2NoDu0whprMrcNqhzFAe4zDMTQHx9AcHEOhtJsLSZIkSZIE5Z5qK0mSJEmShackSZIkqVxNVXhGxIKIuDsiVkTEuUMs3zkivlks/1lE7F+17EPF/Lsj4jXjGfegGEc0hojYPyI2RMSy4ufS8Y69KsZ6Y5gXEbdFxMbiO+Wqly2KiHuKn0WD+46XUY7hmarXYfBNO8ZNA2N4f0T8MiJuj4h/j4j9qpa1yusw3Bha5XX4q4i4o4jzJxFxcNWyVtkvDTmGVtovVbU7MSIyIrqr5jXF6zBRzK2t8Tdsbh0f5tbmeB2KWMyvLbBvqmo3+vyamU3xQ+UmCfcCBwA7Ab8ADh7U5t3ApcXjk4BvFo8PLtrvDMwq1jOpxcawP/BfLfI67A8cBnwVOLFq/vOA+4rfzy0eP7eVxlAsW98ir8MrgV2Lx6dX/S210usw5Bha7HXYrerxCcD3isettF+qNYaW2S8V7aYBNwP/CXQ30+vQzM8d5tZmGcP+mFubYQzm1uYZh/m1CcZQtBuT/NpMn3i+GFiRmfdl5tPAVcAbBrV5A3BF8fhbwHEREcX8qzLzqcz8FbCiWN94G80YmkXdMWTmrzPzdmDToL6vAX6QmWsz8zHgB8CC8Qh6kNGMoVk0MoYfZubvi8n/pPJ9ftBar0OtMTSLRsbwRNXkc4DNd2xrmf3SMGNoFo3sWwHOBz4FPFk1r1leh4libm0O5tbmYG5tHubX5jCu+bWZCs8XAg9UTa8q5g3ZJjM3Ar8F9miw73gYzRgAZkXEzyPiRxFxTNnB1jCa57KVXofhTImIvoj4z4j4k7ENrWHbO4bFwHdH2LcsoxkDtNDrEBFnRMS9VHbKZ21P33EwmjFAi+yXIuJFwD6Z+a/b27fNmVtb5G+4pL5jydzamq9DM+ZWML9Ci+ybxjK/lvY9niMw1JHJwUcFarVppO94GM0Y+oF9M/PRiDgK+E5EHDLoSMl4GM1z2Uqvw3D2zcwHI+IA4KaIuCMz7x2j2BrV8Bgi4m1AN3Ds9vYt2WjGAC30OmTmJcAlEfEW4G+ARY32HQejGUNL7JciYgfgM8A7trdvBzC3tsDfcIl9x5K5tcG+JWuH3Arm15bYN411fm2mTzxXAftUTc8EHqzVJiImA11Uvhy7kb7jYcRjKD6mfhQgM5dSOU/6v5Ue8bON5rlspdehpsx8sPh9H9ALvGgsg2tQQ2OIiFcB5wEnZOZT29N3HIxmDC31OlS5Cth8BLmlXocqW8bQQvulacChQG9E/Bp4KXBdcQOEZnkdJoq5tTX+hsvqO5bMrS30OjR5bgXza6vsm8Y2v2YTXGCclQtUJ1O5UHsWWy9uPWRQmzPY9uYBVxePD2Hbi1vvY2IuMh7NGKZvjpnKBb6rgec14xiq2n6FZ98A4VdULrp/bvG41cbwXGDn4vGewD0McZF1M4yBSrK4Fzho0PyWeR2GGUMrvQ4HVT1eCPQVj1tpv1RrDC23Xyra97L15gdN8TpM1E+Dr7+5tQnGUNX2K5hbJ/JvydzaPOMwvzbBGAa172UU+XXc/9DqDP544P8V/yznFfM+QeVoDcAU4BoqF6/eAhxQ1fe8ot/dwGtbbQzAnwF3Fi/gbcDCJh7D0VSOcvwOeBS4s6rvKcXYVgDvbLUxAH8I3FG8DncAi5t4DP8GPAQsK36ua8HXYcgxtNjrcFHxv7sM+CFVO+wW2i8NOYZW2i8NattLkRib6XVo1ucOc2uzjMHc2hxjMLc2zzjMr00whkFtexlFfo2ikyRJkiRJpWimazwlSZIkSW3IwlOSJEmSVCoLT0mSJElSqSw8Jen/t3f/oVWVcRzH3x+noJuWjEW0EEQqIghbWUS//zAqKjAYgYVURn9EIWoS9AskCKo/8o9ARBcOYQRJhcMorNxKihnbUrbIgkr7R/pBQm4ZpX774zx3nLZ77+7U637czwse7vOcPc9zvmdj++6cc597zMzMzKyqfOJpZmZmZmZmVeUTT7NRJJ2SdEDSoKSdkuqruK+X00OekbR2ovtSZq+kC1J7jaRvJXWcg9geldSca7dJuuoM53pa0mNnG5OZmU1Pzq0jczu3Ws3y41TMRpE0FBHzU70D6IuINyoYJ7LfqdNnuN/DZM9G+n0CY+4FlkfEutQ+RPYcpZ9G9ZsdEScnGE83sCEieicyrsRc9cAXEdFytnOZmdn049w6MqYb51arUb7jaVbePuAyAEnr05XaQUlr07bF6SroZrIHAC+StFLSQOr3WupXJ6k9bRuQVEhm7ZJaJa0BmoEuSV2SHpe0qRCEpCckFUvQDwO7Up8twBKgU9I6SRslbZW0B9iRYt0nqT+Vm3LzP5viOijpVUmtwDKgI12hniepW9Ky1H/MMabtQ5JeSfP0SLoYICL+Ag5LuuGc/FTMzGw6c251brVaFBEuLi65Agyl19lkiedJ4DpgAGgA5gPfAC3AYuA0cGMa0wz8DFyUxu8FVqTxH+f2sTC9tgOtqX4YaEr1BuAHYE5qfwlcXSTWI8CCXDs/x0agD5iX2vXA3FS/HOhN9XvS/PWp3Zheu8muEpNvlzrG1CeA+1P9deDF3PgXgGcm++fr4uLi4nL+i3Orc6uLi+94mo01T9IBoJcsCbwF3AK8HxHDETEEvAfcmvofiYieVL8e6I6I3yJ7+00HcBvwI7BE0puS7gb+LBdARAyTJZ37JF1JliQHinRtjIjjZabqjIgTqT4H2CZpANgJFNaULAe2R3bllIj4o1xsZY4R4B9gd6r3kf3zUPArWWI1M7Pa49xannOrzXizJzsAsynoRERck9+Q1piUMpzvWqxDRByTtBS4C3gKeBBYPU4cbcDzwCFge4k+JyXNitJrX/KxrQN+AZaSvc3+71zME1nsXe578W9EFOY6xf//xswFTowdYmZmNcC5tTznVpvxfMfTrDKfAysk1UtqAB4gW6My2n7gdklNkuqAlcBnkpqAWRHxLt7YQE4AAAFQSURBVPAScG2RsceBBYVGROwHFgEPAW+XiOs7srUnlbgQOJoS6SqgLm3fA6xW+tQ/SY3F4hnvGCvY/xXAYIWxmpnZzOfcOs4xVrB/51abNnziaVaBiOgnWzPyFVlyaIuIr4v0Owo8B3QBB4H+iNgFXAp0p7cZtac+o20FPpTUldv2Dtkn1h0rEdoHwB0VHsZm4BFJPWSJajjF/BHQCfSm+Dak/u3AlsIHIFRwjOO5GfikwljNzGyGc251brXa4sepmE1hknYDmyLi0xJfvwTYERF3nt/IJkZSC7A+IlZNdixmZlbbnFvNJofveJpNQZIWSvqebE1M0cQII1dItyk95HoKayJ7G5SZmdmkcG41m1y+42lmZmZmZmZV5TueZmZmZmZmVlU+8TQzMzMzM7Oq8omnmZmZmZmZVZVPPM3MzMzMzKyqfOJpZmZmZmZmVfUfyKfew6xVTi0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "stdev1 = np.std(por1); stdev2 = np.std(por2)\n",
    "por1r05 = GSLIB.affine(por1,tmean=ci_95_por1[0],tstdev=stdev1)\n",
    "por1r95 = GSLIB.affine(por1,tmean=ci_95_por1[1],tstdev=stdev1)\n",
    "por2r05 = GSLIB.affine(por2,tmean=ci_95_por2[0],tstdev=stdev2)\n",
    "por2r95 = GSLIB.affine(por2,tmean=ci_95_por2[1],tstdev=stdev2)\n",
    "\n",
    "plt.subplot(121)\n",
    "plt.hist(por1r05, facecolor='yellow',bins=np.linspace(0.0,0.4,50),histtype=\"stepfilled\",alpha=0.8,density=True,cumulative=True,edgecolor='black',label='P05')\n",
    "plt.hist(por1, facecolor='orange',bins=np.linspace(0.0,0.4,50),histtype=\"stepfilled\",alpha=0.8,density=True,cumulative=True,edgecolor='black',label='Expected')\n",
    "plt.hist(por1r95, facecolor='red',bins=np.linspace(0.0,0.4,50),histtype=\"stepfilled\",alpha=0.8,density=True,cumulative=True,edgecolor='black',label='P95')\n",
    "\n",
    "plt.xlim([0.0,0.4]); plt.ylim([0,1.0])\n",
    "plt.xlabel('Porosity (fraction)'); plt.ylabel('Cumulative Probability'); plt.title('Porosity Distribution Uncertainty - Well 1')\n",
    "plt.legend(loc='upper left')\n",
    "plt.grid(True)\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.hist(por2r05, facecolor='yellow',bins=np.linspace(0.0,0.4,50),histtype=\"stepfilled\",alpha=0.8,density=True,cumulative=True,edgecolor='black',label='P05')\n",
    "plt.hist(por2, facecolor='orange',bins=np.linspace(0.0,0.4,50),histtype=\"stepfilled\",alpha=0.8,density=True,cumulative=True,edgecolor='black',label='Expected')\n",
    "plt.hist(por2r95, facecolor='red',bins=np.linspace(0.0,0.4,50),histtype=\"stepfilled\",alpha=0.8,density=True,cumulative=True,edgecolor='black',label='P95')\n",
    "\n",
    "plt.xlim([0.0,0.4]); plt.ylim([0,1.0])\n",
    "plt.xlabel('Porosity (fraction)'); plt.ylabel('Cumulative Probability'); plt.title('Porosity Distribution Uncertainty - Well 2')\n",
    "plt.legend(loc='upper left')\n",
    "plt.grid(True)\n",
    "\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.2, hspace=0.3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We just calculated a scenario-based uncertainty model for the porosity distribution around wells 1 and 2. Of course, we could actually sample porosity means continuously from our confidence calculation as we have access to the complete distribution for uncertainty in the mean porosity of both wells.\n",
    "\n",
    "Communicating uncertainty is powerful, but always remember to state the assumptions. For example, here we assumed:\n",
    "\n",
    "1. The population distribution of porosity for each well is Gaussian distributed\n",
    "2. That the samples are independent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can check the Excel file linked above with the confidence interval calculated by hand and confirm that this result is correct.\n",
    "\n",
    "#### Are the Porosity Values from Unit 1 and Unit 2 Significantly Different?\n",
    "\n",
    "Let's first visually compare the porosity from the two units.  We will use a box and wisker plot.  \n",
    "\n",
    "* Box and Whiskers plot is a built in function with DataFrames\n",
    "\n",
    "```python\n",
    "df.boxplot()\n",
    "```\n",
    "\n",
    "* We return the 'plt' label so that we can further improve the plot\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'Porosity (fraction)')"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEXCAYAAACH/8KRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHzJJREFUeJzt3Xu0XlV97vHvY+SmXASTDoEQAppRRdHQboKjVlCMGuwxcFoUsCgoNeIoLT0clXCkWFHaiNZrUclRMIqIQKtECSdaDFqrQMI1JEgNMcI2IMFyvwQSnvPHWhsXL+/ee+299tqX5PmM8Y79rrnmnGu+m5f9y5xzrTllm4iIiOF6zlg3ICIiJrYEkoiIaCSBJCIiGkkgiYiIRhJIIiKikQSSiIhoJIEktgiSvizp78e6HXVIukrSX411O0aKpIcl7TvW7Yixk0AS44IkS3pJR9o/SLqgTnnbJ9r+WFnudZJ6B7ne6yUtk/SApHXDbvgIk3ScpOskPSipV9LZkp47QH5LeqT8Y/4bSZ+WNGk022x7R9try/Z8TdLHR/P6MfYSSGJr9QhwHvDBsW5Ih+cBfwdMBg4C3gB8YJAyr7K9Y5n3HcB7h3rR0Q4+sWVJIIkJoa+XIel/S7pH0l2S3l05/zVJH5f0fOAKYI/yX+kPS9qjsz7b19r+BrC25vUvkXR32YP5iaSXd1z7HEmXS3pI0jWSXlw5/0ZJvyjL/gug/q5j+0u2/8P2E7Z/A3wTeE2dNtr+BfAfwCvK676sHEa7X9IqSXM72vwlSUskPQK8XtIukr4uaYOkX0s6XdJzyvwvkfTj8jPcK+nblbpcnp8H/CXwofL3/j1JH5T0rx2/yy9I+mydzxQTQwJJTCQvAnYB9gROAM6RtGs1g+1HgMOA9eWQy46214/Ata8AZgB/AFxP8Qe+6hjgo8CuwBrgLABJk4F/BU6n6GXcTs3AUDoYWFUno6T9gNcCN0jaBvge8IOyzX8DfFPSH1aKvKNs507AT4EvUPx+9wUOAd4F9AXrj5V17QpMLfM+g+2FFL+Xs8vf+1uBC4A5kl5QtvG5wFHAN2p+/pgAEkhiInkSONP2k7aXAA8DfzhImRFh+zzbD9neCPwD8CpJu1Sy/FvZy9lE8cd0Zpn+FmC17UttPwl8Fri7zjXLHlcP8KlBsl4v6T6KwPEV4Hzg1cCOwIKyd/Mj4PsUAa/PZbb/0/ZTFL/bo4DTys+5Dvhn4J1l3ieBvYE9bD9u+6d1PoPtu4CfAG8rk+YA99q+rk75mBgSSGK82Axs05G2DcUfsD6/K/9Q93mU4o9lqyRNkrRA0u2SHgTWlacmV7JVg0O1XXsAd/adcLFK6p0MQtIRwALgMNv3DpL9j2zvavvFtk8vA8MewJ3l+z6/pujN9am2YzKwbZmnW/4PUQzJXVsOk71nsM9QsQg4tnx/LOmNbHESSGK8uAOY3pG2D8/8w1bXSC9p/Q7gcGA2xdDP9DK937mOiruAvfoOJKl63I2kOcD/Bd5qe+Uw2guwHtirb46jNA34TeW4+nu6l9/3Op6V3/bdtt9rew/gfcAXO++y61Jnn+8Cr5T0CuB/8OxhwZjgEkhivPg2cLqkqZKeI2k28Fbg0mHU9VvghR1DT89QXmN7il6PJG0vadt+su8EbAR+R3FX1T8OoS2XAy+X9Ofl/MDfUsz19NeuQyn+0P6F7WuHcJ1O11DcmfYhSdtIeh3F7/OibpltbwYuBs6StJOkvYFTKOY4kPQ2SVPL7PdRBIzNXar6LcUcS7Xuxyn+O14IXGv7jgafK8ahBJIYL84EfkYx6XsfcDbwl7ZvGWpF5d1L3wLWlncsPeuuLYpJ7MeAJRT/8n6MYjK5m69T9Ix+A6wGrh5CW+6lmB9YQBGIZgD/OUCRv6fo9Syp3HV2Rd3rVa77BDCX4saDe4EvAu8qfzf9+RuK4LOW4r/DhRS3SAMcCFwj6WFgMXCy7V91qeOrwH7l7/27lfRFwP5kWGuLpGxsFRFtkzQN+AXwItsPjnV7YmSlRxIRrSrnaU4BLkoQ2TL1u/RCRERT5QOiv6UYGpwzxs2JlmRoKyIiGsnQVkRENLJVDG1NnjzZ06dPH+tmRERMKNddd929tqcMlm+rCCTTp09nxYoVY92MiIgJRVKtB4IztBUREY0kkERERCMJJBER0UgCSURENJJAEhERjbQaSCTNkXSbpDWS5nc5f6KklZJulPTTcoe3vnOnleVuk/TmunVGRMToai2QSJoEnEOx+uh+wDHVQFG60Pb+tmdSrPb66bLsfsDRwMspllX4Yrm5UJ06IyJiFLXZI5kFrLG9tlzS+iKKzYGe1rGA2/P5/aY4h1Ms8LaxXKp6TVnfoHVGRMToavOBxD155laevcBBnZkk/TXFyqDbAodWylb3fOjl91t+DlpnWe88YB7AtGnTht76iBi3io0mhybrCranzR5Jt//Sz/ovafsc2y8GTgVOH6RsrTrLehfa7rHdM2XKoE/4R8QEYrvra+9Tv9/vuWhPmz2SXp65N/VUin2k+3MR8KUaZYdSZ0REtKzNHslyYIakfcq9sI+m2KLzaZJmVA7/DPhl+X4xcLSk7STtQ7E96bV16oyIiNHVWo/E9iZJJwFLgUnAebZXSToTWGF7MXCSpNnAkxT7dB9Xll0l6WKK/bE3AX9tezNAtzrb+gwRETG4Vlf/tb0EWNKRdkbl/ckDlD0LOKtOnRERMXbyZHtERDSSQBIREY0kkERERCMJJBER0UgCSURENJJAEhERjSSQREREIwkkERHRSAJJREQ0kkASERGNJJBEREQjCSQREdFIAklERDSSQBIREY0kkERERCMJJBER0UgCSURENNJqIJE0R9JtktZImt/l/CmSVku6WdKVkvYu018v6cbK63FJR5TnvibpV5VzM9v8DBERMbDWttqVNAk4B3gj0Assl7TY9upKthuAHtuPSno/cDZwlO1lwMyynt2ANcAPKuU+aPvSttoeERH1tdkjmQWssb3W9hPARcDh1Qy2l9l+tDy8GpjapZ4jgSsq+SIiYhxpM5DsCdxZOe4t0/pzAnBFl/SjgW91pJ1VDod9RtJ23SqTNE/SCkkrNmzYMJR2R0TEELQZSNQlzV0zSscCPcAnO9J3B/YHllaSTwNeChwI7Aac2q1O2wtt99jumTJlytBbHxERtbQZSHqBvSrHU4H1nZkkzQY+DMy1vbHj9NuB79h+si/B9l0ubATOpxhCi4iIMdJmIFkOzJC0j6RtKYaoFlczSDoAOJciiNzTpY5j6BjWKnspSBJwBHBLC22PiIiaWrtry/YmSSdRDEtNAs6zvUrSmcAK24sphrJ2BC4p4gJ32J4LIGk6RY/mxx1Vf1PSFIqhsxuBE9v6DBERMbjWAgmA7SXAko60MyrvZw9Qdh1dJudtHzqCTYyIiIbyZHtERDSSQBIREY0kkERERCMJJBER0UgCSURENJJAEhERjSSQREREIwkkERHRSAJJREQ0kkASERGNJJBEREQjCSQREdFIAklERDSSQBIREY0kkERERCMJJBER0UgCSURENNJqIJE0R9JtktZImt/l/CmSVku6WdKVkvaunNss6cbytbiSvo+kayT9UtK3y/3gIyJijLQWSCRNAs4BDgP2A46RtF9HthuAHtuvBC4Fzq6ce8z2zPI1t5L+CeAztmcA9wEntPUZIiJicG32SGYBa2yvtf0EcBFweDWD7WW2Hy0PrwamDlShJAGHUgQdgEXAESPa6oiIGJI2A8mewJ2V494yrT8nAFdUjreXtELS1ZL6gsULgfttbxqsTknzyvIrNmzYMLxPEBERg3pui3WrS5q7ZpSOBXqAQyrJ02yvl7Qv8CNJK4EH69ZpeyGwEKCnp6drnoiIaK7NHkkvsFfleCqwvjOTpNnAh4G5tjf2pdteX/5cC1wFHADcC7xAUl8A7FpnRESMnjYDyXJgRnmX1bbA0cDiagZJBwDnUgSReyrpu0rarnw/GXgNsNq2gWXAkWXW44DLWvwMERExiNYCSTmPcRKwFLgVuNj2KklnSuq7C+uTwI7AJR23+b4MWCHpJorAscD26vLcqcApktZQzJl8ta3PEBERgxt0jkRSD/BaYA/gMeAW4N9t//dgZW0vAZZ0pJ1ReT+7n3I/A/bv59xaijvCIiJiHOi3RyLpeEnXA6cBOwC3AfcAfwr8UNIiSdNGp5kRETFeDdQjeT7wGtuPdTspaSYwA7ijjYZFRMTE0G8gsX3OQAVt3zjyzYmIiImmzhzJFOC9wPRqftvvaa9ZERExUdR5IPEy4D+Afwc2t9uciIiYaOoEkufZPrX1lkRExIRU5zmS70t6S+stiYiICalOIDmZIpg8Lumh8tVtzauIiNgKDTq0ZXun0WhIRERMTLVW/y2XNDm4PLzK9vfba1JEREwkgw5tSVpAMby1unydXKZFRETU6pG8BZhp+ykASYsotsh91h7sERGx9am7+u8LKu93aaMhERExMdXpkfwTcIOkZRS7Hh5MsZBjRERErbu2viXpKuBAikByqu27225YRERMDP0GEkkvtf0LSX9UJvWWP/eQtIft69tvXkRsrV710R/wwGNPDrnc9PmX1867yw7bcNNH3jTka8QzDdQjOQWYB/xzl3MGDm2lRRERwAOPPcm6BX/W6jWGEnSifwMtIz+vfHuY7cer5yRtX6dySXOAzwGTgK/YXtBx/hTgr4BNwAbgPbZ/Xe518iVgZ4qFIs+y/e2yzNeAQ4AHymqOz5L2ERFjp85dWz+rmfYMkiYB5wCHAfsBx0jaryPbDUCP7VcClwJnl+mPAu+y/XJgDvBZSdU7xz5oe2b5ShCJiBhDA82RvAjYE9hB0gEUE+1Q9BKeV6PuWcCaco91JF0EHE7xUCMAtpdV8l8NHFum/1clz3pJ9wBTgPtrXDciIkbRQHMkbwaOB6ZSzJP0BZIHgf9To+49gTsrx73AQQPkPwG4ojNR0ixgW+D2SvJZks4ArgTm297Ypdw8ijkepk3L1vIREW0ZaI5kEbBI0l/Y/tdh1K0uae6aUToW6KGY+6im7w58Aziu78l6imdY7qYILguBU4Ezu7R/YXmenp6erteNiIjm6syR/HF1fkLSrpI+XqNcL7BX5XgqsL4zk6TZwIeBudWehaSdgcuB021f3Zdu+y4XNgLnUwyhRUTEGKkTSA6z/fTchO37KNbfGsxyYIakfSRtCxwNLK5mKOdezqUIIvdU0rcFvgN83fYlHWV2L38KOAK4pUZbIiKiJXWWSJkkabu+3oKkHYDtBitke5Okk4ClFLf/nmd7laQzgRW2FwOfBHYELiniAnfYngu8nWIplhdKOr6ssu82329KmkIxdHYjcGL9jxsRESOtTiC5ALhS0vkUcxzvARbVqdz2EmBJR9oZlfez+yl3QXndbufyIGRExDhSZ62tsyWtBN5A0Qv4mO2lrbcsIiImhFo7JNq+gi635kZERNTZIfHVkpZLeljSE5I2S3pwNBoXERHjX527tv4FOAb4JbADxdpYX2izURERMXHUHdpaI2mS7c3A+ZIGXWsrIiK2DnUCyaPlcx03SjobuAt4frvNioiIiaLO0NY7y3wnAY9QPK3+F202KiIiJo4BeyTlUvBn2T4WeBz46Ki0KiIiJowBeyTlnMiUcmgrIiLiWerMkawD/lPSYoqhLQBsf7qtRkVExMRRJ5CsL1/PAXZqtzkxXpRrnw2ZnRX7I7Y2A+2Q+A3b7wTut/25UWxTjAMDBYTp8y9n3YI/G8XWRMR4NtAcyR9L2ht4T7kHyW7V12g1MCIixreBhra+DPw/YF/gOp6546HL9IiI2Mr12yOx/XnbL6PYR2Rf2/tUXgkiEREBDBBIJO0IYPv9g+WJiIit10BzJJdJ+mdJB0t6ekkUSftKOkHSUmBO+02MiIjxbKChrTcAVwLvA1ZJekDS7yh2LnwRcJztSweqXNIcSbdJWiNpfpfzp0haLelmSVeWk/t9546T9MvydVwl/Y8lrSzr/LyGe59qRESMiAGfI+m2VW5d5fIq5wBvBHqB5ZIW215dyXYD0GP7UUnvB84GjirvCvsI0EMxsX9dWfY+4EvAPODqsm1zyKZbERFjps6ijcM1C1hje63tJ4CLgMOrGWwvs/1oeXg1MLV8/2bgh7b/uwwePwTmSNod2Nn2z1086PB14IgWP0NERAyi1n4kw7QncGfluBc4aID8J/D7nkW3snuWr94u6c8iaR5Fz4Vp06YNpd0RMQ7s9LL57L/oWSPiI3wNgDxc21SbgaTb3EXXx6UlHUsxjHXIIGVr12l7IbAQoKenJ+t2REwwD926oPUVFKbPv7zV+rcWdfZs/5Sklw+j7l6KvUv6TKVYs6uz/tnAh4G5tjcOUraX3w9/9VtnRESMnjpzJL8AFkq6RtKJknapWfdyYIakfcpl6I8GFlczSDoAOJciiNxTObUUeFO5NMuuwJuApbbvAh6S9Orybq13AZfVbE9ERLRg0EBi+yu2X0PxR3s6cLOkCyW9fpBymyh2VVwK3ApcbHuVpDMlzS2zfRLYEbhE0o3lUvXY/m/gYxTBaDlwZpkG8H7gK8Aa4HZyx1ZExJiqNUdS3sr70vJ1L3ATcIqk99k+ur9y3W4ftn1G5f3sAcqeB5zXJX0F8Io67Y6IiPYNGkgkfRp4K/Aj4B9tX1ue+oSk29psXEREjH91eiS3AKdXnveomjXC7YmIiAmmzmT7X3YGEUlXAth+oJVWRUTEhDHQDonbA88DJpd3TvU9w7EzsMcotC0iIiaAgYa23gf8HUXQuL6S/iDFGloRERH9B5Jyn/bPSfob218YxTbFKHrVR3/AA489OeRyQ3kieJcdtuGmj7xpyNeIiIlhoKGtQ23/CPiNpD/vPG/731ptWYyKBx57MstQREQjAw1tHUJxy+9bu5wzkEASEREDDm19pPz57tFrTkRETDR1Fm08WdLOKnxF0vWSMuAdERFAvedI3mP7QYqFE/8AeDewoNVWRUTEhFEnkPQ9P/IW4HzbN9F9X5CIiNgK1Qkk10n6AUUgWSppJ+CpdpsVERETRZ21tk4AZgJrbT8q6YUUw1sRERGDBxLbT0maCryj2EuKH9v+Xusti4iICaHOXVsLgJOB1eXrbyX9U9sNi4iIiaHO0NZbgJm2nwKQtAi4ATitzYZFRMTEUGeyHeAFlfd192xH0hxJt0laI2l+l/MHl8+lbJJ0ZCX99eXWu32vxyUdUZ77mqRfVc7NrNueiIgYeXV6JP8E3CBpGcVtvwdTozdSbs97DvBGoBdYLmmx7dWVbHcAxwMfqJa1vYxigh9Ju1Hsz/6DSpYP2r60RtsjIqJlAwYSFbPrPwVeDRxIEUhOtX13jbpnAWtsry3rugg4nGKeBQDb68pzA91OfCRwRT87NEZExBgbcGjLtoHv2r7L9mLbl9UMIgB7AndWjnvLtKE6GvhWR9pZkm6W9BlJ23UrJGmepBWSVmzYsGEYl42IiDrqzJFcLenAYdTd7el3D6kCaXdgf2BpJfk04KUUPaTdgFO7lbW90HaP7Z4pU6YM5bIRETEEdQLJ64GfS7q97AWslHRzjXK9wF6V46nA+iG27+3Ad2w/vfNS2Tuy7Y3A+RRDaBERMUbqTLYfNsy6lwMzJO0D/IZiiOodQ6zjGDom9iXtbvuucv7mCOCWYbYvIiJGwKA9Etu/prj9963l6wVl2mDlNgEnUQxL3QpcbHuVpDMlzQWQdKCkXuBtwLmSVvWVlzSdokfz446qvylpJbASmAx8fLC2REREewbtkUg6GXgvv98R8QJJC+vs4257CbCkI+2MyvvlFENe3cquo8vkvO1DB7tuRESMnrqLNh5k+xEASZ8Afg4MGkhi/NvpZfPZf9GznhUd4WsAtLsvfESMnTqBRMDmyvFmsh/JFuOhWxewbkG7f+Snz7+81fojYmzVCSTnA9dI+k55fATw1faaFBFRaPsfIbvssE2r9W8t6iwj/2lJVwF/StETebftG9puWERs3YbTU54+//LWe9jxbP0GEknbAycCL6G4Q+qL5Z1YERERTxvo9t9FQA9FEDkM+NSotCgiIiaUgYa29rO9P4CkrwLXjk6TIiJiIhmoR1JdliRDWhER0dVAPZJXSXqwfC9gh/JYFAsD79x66yIiYtzrN5DYnjSaDYmIiImp7la7ERERXSWQREREIwkkERHRSAJJREQ0kkASERGNJJBEREQjCSQREdFIq4FE0hxJt0laI+lZuydJOljS9ZI2STqy49xmSTeWr8WV9H0kXSPpl5K+LWnbNj9DREQMrLVAImkScA7Fgo/7AcdI2q8j2x3A8cCFXap4zPbM8jW3kv4J4DO2ZwD3UezgGBERY6TNHsksYI3ttbafAC4CDq9msL3O9s3AU3UqlCTgUODSMmkRxUZbERExRurskDhcewJ3Vo57gYOGUH57SSuATcAC298FXgjcX1lEsre8zrNImgfMA5g2bdoQm751yS50EdFEm4Gk277uHkL5abbXS9oX+JGklcCDXfJ1rdP2QmAhQE9Pz1Cuu1XJLnQR0VSbQ1u9wF6V46nA+rqFba8vf64FrgIOAO4FXiCpLwAOqc6IiBh5bQaS5cCM8i6rbYGjgcWDlAFA0q6StivfTwZeA6y2bWAZ0HeH13HAZSPe8oiIqK21QFLOY5wELAVuBS62vUrSmZLmAkg6UFIv8DbgXEmryuIvA1ZIuokicCywvbo8dypwiqQ1FHMmX23rM0RExODanCPB9hJgSUfaGZX3yymGpzrL/QzYv58611LcERYREeNAnmyPiIhGEkgiIqKRBJKIiGgkgSQiIhpJIImIiEYSSCIiopEEkoiIaCSBJCIiGkkgiYiIRhJIIiKikQSSiIhoJIEkIiIaSSCJiIhGEkgiIqKRBJKIiGgkgSQiIhpJIImIiEZaDSSS5ki6TdIaSfO7nD9Y0vWSNkk6spI+U9LPJa2SdLOkoyrnvibpV5JuLF8z2/wMERExsNa22pU0CTgHeCPQCyyXtLiy9zrAHcDxwAc6ij8KvMv2LyXtAVwnaant+8vzH7R9aVttj4iI+trcs30WsKbcYx1JFwGHA08HEtvrynNPVQva/q/K+/WS7gGmAPcTERHjSptDW3sCd1aOe8u0IZE0C9gWuL2SfFY55PUZSdv1U26epBWSVmzYsGGol42IiJraDCTqkuYhVSDtDnwDeLftvl7LacBLgQOB3YBTu5W1vdB2j+2eKVOmDOWyERExBG0Gkl5gr8rxVGB93cKSdgYuB063fXVfuu27XNgInE8xhBYREWOkzUCyHJghaR9J2wJHA4vrFCzzfwf4uu1LOs7tXv4UcARwy4i2OiIihqS1QGJ7E3ASsBS4FbjY9ipJZ0qaCyDpQEm9wNuAcyWtKou/HTgYOL7Lbb7flLQSWAlMBj7e1meIiIjBtXnXFraXAEs60s6ovF9OMeTVWe4C4IJ+6jx0hJsZEREN5Mn2iIhoJIEkIiIaSSCJiIhGEkgiIqKRBJKIiGgkgSQiIhpJIImIiEYSSCIiopEEkoiIaCSBJCIiGml1iZSYuIo1MQc4/4nu6faQdgqIGJaBvp/5bo6+BJLoKv/TxXiW7+f4kqGtiIhoJIEkIiIaSSCJiIhGEkgiIqKRBJKIiGgkgSQiIhpJIImIiEYSSCIiohFtDQ/2SNoA/Hqs27EFmQzcO9aNiOgi382RtbftKYNl2ioCSYwsSSts94x1OyI65bs5NjK0FRERjSSQREREIwkkMRwLx7oBEf3Id3MMZI4kIiIaSY8kIiIaSSCJiIhGEkgiIqKRBJItkKTpkm7pSPsHSR8YoEyPpM+X718n6U8q5w6WdL2kTZKOrHHtxyTdKGm1pC9LyvcsntbC9/OU8rt2s6QrJe09yLXz/Rxh+QUGALZX2P7b8vB1wJ9UTt8BHA9cWLO6223PBF4J7AccUaeQpEk164+tzCDfzxuAHtuvBC4Fzh6kunw/R1gCyVZG0lWSPiHpWkn/Jem1ZfrrJH1f0nTgROB/lf9qe63tdbZvBp4ayrVsbwJ+BrxEhU9KukXSSklHVa67TNKFwMqR/Kwx8Qzz+7nM9qNlFVcDU+tcK9/PkfPcsW5AjInn2p4l6S3AR4DZfSdsr5P0ZeBh259qchFJzwPeAJwB/DkwE3gVxXpIyyX9pMw6C3iF7V81uV5sMZp8P08ArqhzkXw/R04CyZapv4eD+tL/rfx5HTC9heu/WNKN5fUus32FpM8A37K9GfitpB8DBwIPAtfmf9KtSivfT0nHAj3AIYNkzfdzhCWQbJl+B+zakbYb0Pc/w8by52ba+Q70jUFXaYD8j7TQhhi/Rvz7KWk28GHgENsbB8me7+cIyxzJFsj2w8Bdkt4AIGk3YA7w05pVPATsNMLN+glwlKRJkqYABwPXjvA1YgIY6e+npAOAc4G5tu8ZZrPy/WwggWTL9S7g9LIL/yPgo7Zvr1n2e8D/7JvMlHSgpF7gbcC5klYNoz3fAW4Gbirb8yHbdw+jntgyjNj3E/gksCNwSZm2eBjtyfezgay1FRERjaRHEhERjWSyPYZF0v7ANzqSN9o+aCzaE1GV7+foytBWREQ0kqGtiIhoJIEkIiIaSSCJaMlIr3IbMV5lsj1iHLG9AlhRHr4OeJhiYcGIcSs9kogxMJxVbseyvREDSY8kYuyMyirMEW1LjySiPWO9CnPEqEggiWhPf6vc3lu+b3sV5ohRkUAS0ZJxugpzxIhLIIlo10iuchsxLmWJlIiIaCQ9koiIaCSBJCIiGkkgiYiIRhJIIiKikQSSiIhoJIEkIiIaSSCJiIhG/j9XwlPKyVd54gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt = df.boxplot(grid = False,fontsize = 10)\n",
    "plt.set_title('Unit 1 and 2 Porosity'); plt.set_xlabel(\"Unit\"); plt.set_ylabel('Porosity (fraction)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the porosity distributions seem quite different between the two units.  \n",
    "\n",
    "* Let's check and see if they are **statistically significantly different**\n",
    "\n",
    "#### Hypothesis Testing of Difference Between Unit 1 and 2\n",
    "\n",
    "The confidence intervals help with uncertainty in the distributions of porosity.  Now let's try to figure out if:\n",
    "\n",
    "1. could we combine units 1 and 2? \n",
    "2. did something change between the 2 units? \n",
    "\n",
    "Now, let's try the t test, hypothesis test for difference in means.  This test assumes that the variances are similar along with the data being Gaussian distributed (see the course notes for more on this).  This is our test:\n",
    "\n",
    "\\begin{equation}\n",
    "H_0: \\mu_{X1} = \\mu_{X2}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "H_1: \\mu_{X1} \\ne \\mu_{X2}\n",
    "\\end{equation}\n",
    "\n",
    "For the resulting t-statistic and p-value we run this command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The t statistic is -2.981 and the p-value is 0.005\n"
     ]
    }
   ],
   "source": [
    "t_pooled, p_pooled = st.ttest_ind(por1,por2)\n",
    "print('The t statistic is ' + str(round(t_pooled,3)) + ' and the p-value is ' + str(round(p_pooled,3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The p-value, $p$, is the symmetric interval probaiblity our outside.  In other words the $p$ reported is 2 x cumulative probaiblity of the t statistic applied to the sampling t distribution.  Another way to look at it, if one used the $\\pm t_{t_{statistic},.d.f}$ statistic as thresholds, $p$ is the probability being outside this symmetric interval. So we will reject the null hypothesis if $p \\lt \\alpha$.  From the p-value alone it is clear that we would reject the null hypothesis and accept the alternative hypothesis that the means are not equal.  \n",
    "\n",
    "In case you want to compare the t-statistic to t-critical, we can apply the inverse of the student's t distribution at $\\frac{\\alpha}{2}$ and $1-\\frac{\\alpha}{2}$ to get the upper and lower critcal values.       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The t crical lower and upper values are [-2.02  2.02]\n"
     ]
    }
   ],
   "source": [
    "t_critical = st.t.ppf([0.025,0.975], df=len(por1)+len(por2)-2)\n",
    "print('The t crical lower and upper values are ' + str(np.round(t_critical,2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can observe that, as expected, the t-statistic is outside the t-critcal interval.  These results are exactly what we got when we worked out the problem by hand in Excel, but so much more efficient!  \n",
    "\n",
    "**The porosity from the 2 units are significantly different.**\n",
    "\n",
    "Now let's try the t-test, hypothesis test for difference in means allowing for unequal variances, this is also known as the Welch's t test.  All we have to do is set the parameter 'equal_var' to false, note it defaults to true (e.g. the command above). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ttest_indResult(statistic=-2.9808897468855644, pvalue=0.005502572350112333)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "st.ttest_ind(por1, por2, equal_var = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once again we can see by $p$ that we will clearly reject the null hypothesis.  \n",
    "\n",
    "Let's now compare the variances with the F-test for difference in variances.  \n",
    "\n",
    "\\begin{equation}\n",
    "H_0: \\sigma^{2}_{X2} \\le \\sigma^{2}_{X1}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "H_1: \\sigma^{2}_{X2} \\gt \\sigma^{2}_{X1}\n",
    "\\end{equation}\n",
    "\n",
    "Note, by ordering the variances we eliminate the case of $\\sigma^{2}_{X2} \\lt \\sigma^{2}_{X1}$.\n",
    "\n",
    "Details about the test are available in the course notes (along with assumptions such as Gaussian distributed) and this example is also worked out by hand in the linked Excel workbook.  We can accomplish the F-test in with SciPy.Stats the function with one line of code if we calculate the ratio of the sample variances ensuring that the larger variance is in the numerator and get the degrees of freedom using the len() command, ensuring that we are consistent with the numerator degrees of freedom set as 'dfn' and the denominator degrees of freedom set as 'dfd'.  We take a p-value of $1-p$ since the test is configured to be a single, right tailed test.    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The p-value for the F-test for difference is variances is 0.019\n"
     ]
    }
   ],
   "source": [
    "p_value = 1 - st.f.cdf(np.var(por2)/np.var(por1), dfn=len(por2)-1, dfd=len(por1)-1)\n",
    "print('The p-value for the F-test for difference is variances is ' + str(round(p_value,3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once again we would clearly reject the null hypothesis since $p \\lt alpha$ and assume that the variances are not equal.  Note this is a single tail test so our p-critical is 0.05 for a 95% confidence level.\n",
    "\n",
    "#### Reporting the Hypothesis Test Results\n",
    "\n",
    "The difference in the means and variances were tested for statistical significance for porosity measured in wells 1 and 2.  In all cases we rejected the null hypotheses and adopted the alternative hypotheses that the porosity distributions' means and variances for well 1 and 2 are statistically significant in their differences. It is possible that Well 1 and 2 have drilled into different rock units or have encountered nonstationary behavoirs in a single rock unit. This result may be applied to update subsurface maps and consider a difference in the subsurface at each well for development decision making."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comments\n",
    "\n",
    "We are just scratching the surface for the application of confidence intervals and hypothesis tests to subsurface modeling.  Once again, there are a lot of details left out of the problem formulation and assumptions, see the course notes for more coverage.  By running the same confidence interval and hypothesis tests: 1) by hand in Excel and with 2) R and now here in 3) Python code, I hope this demonstration will enable and encourage more engineers and scientists to make these R and Python tools part of their common practice. I'm always happy to discuss,\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"
   ]
  }
 ],
 "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
}