PythonDataBasic_Univariate_Statistics.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",
    "## Data Analytics \n",
    "\n",
    "### Basic Univariate Statistics 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"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data Analytics: Basic Univariate Statistics\n",
    "\n",
    "Here's a demonstration of calculation of univariate statistics in Python. This demonstration is part of the resources that I include for my courses in Spatial / Subsurface Data Analytics and Geostatistics at the Cockrell School of Engineering and Jackson School of Goesciences at the University of Texas at Austin.  \n",
    "\n",
    "We will cover the following statistics:\n",
    "\n",
    "#### Measures of Centrality\n",
    "* Arithmetic Average / Mean\n",
    "* Median\n",
    "* Mode (most frequent binned)\n",
    "* Geometric Mean\n",
    "* Harmonic Mean\n",
    "* Power Law Average\n",
    "\n",
    "#### Measures of Dispersion\n",
    "* Population Variance\n",
    "* Sample Variance\n",
    "* Population Standard Deviation\n",
    "* Sample Standard Deviation\n",
    "* Range\n",
    "* Percentile w. Tail Assumptions\n",
    "* Interquartile Range\n",
    "\n",
    "#### Tukey Outlier Test\n",
    "* Lower Quartile/P25\n",
    "* Upper Quartile/P75\n",
    "* Interquartile Range\n",
    "* Lower Fence\n",
    "* Upper Fence\n",
    "* Calculating Outliers\n",
    "\n",
    "#### Measures of Shape\n",
    "* Skew\n",
    "* Excess Kurtosis\n",
    "* Pearson' Mode Skewness\n",
    "* Quartile Skew Coefficient\n",
    "\n",
    "#### Nonparmetric Cumulative Distribution Functions (CDFs)\n",
    "* plotting a nonparametric CDF\n",
    "* fitting a parametric distribution and plotting\n",
    "\n",
    "I have a lecture on these univariate statistics available on [YouTube](https://www.youtube.com/watch?v=wAcbA2cIqec&list=PLG19vXLQHvSB-D4XKYieEku9GQMQyAzjJ&index=11&t=0s).   \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.  The dataset is available on my GitHub account in my GeoDataSets repository at:\n",
    "\n",
    "* Tabular data - [2D_MV_200wells.csv](https://github.com/GeostatsGuy/GeoDataSets/blob/master/2D_MV_200wells.csv)\n",
    "\n",
    "#### Importing Packages\n",
    "\n",
    "We will need some standard packages. These should have been installed with Anaconda 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np                        # ndarrys for gridded data\n",
    "import pandas as pd                       # DataFrames for tabular data\n",
    "import os                                 # set working directory, run executables\n",
    "import matplotlib.pyplot as plt           # plotting\n",
    "from scipy import stats                   # summary statistics\n",
    "import scipy                              # statistics\n",
    "import statistics as stats                # statistics like the mode\n",
    "from scipy.stats import norm              # fitting a Gaussian distribution"
   ]
  },
  {
   "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). Set this to your working directory, with the above mentioned data file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(\"d:/PGE383\")                     # set the working directory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Loading Data \n",
    "\n",
    "Let's load the provided multivariate, spatial dataset.  '2D_MV_200wells.csv' is available at https://github.com/GeostatsGuy/GeoDataSets.  It is a comma delimited file with X and Y coordinates,facies 1 and 2 (1 is sandstone and 2 interbedded sand and mudstone), porosity (fraction), permeability (mDarcy) and acoustic impedance (kg/m2s*10^6). 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)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>X</th>\n",
       "      <th>Y</th>\n",
       "      <th>facies_threshold_0.3</th>\n",
       "      <th>porosity</th>\n",
       "      <th>permeability</th>\n",
       "      <th>acoustic_impedance</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>565</td>\n",
       "      <td>1485</td>\n",
       "      <td>1</td>\n",
       "      <td>0.1184</td>\n",
       "      <td>6.170</td>\n",
       "      <td>2.009</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2585</td>\n",
       "      <td>1185</td>\n",
       "      <td>1</td>\n",
       "      <td>0.1566</td>\n",
       "      <td>6.275</td>\n",
       "      <td>2.864</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2065</td>\n",
       "      <td>2865</td>\n",
       "      <td>2</td>\n",
       "      <td>0.1920</td>\n",
       "      <td>92.297</td>\n",
       "      <td>3.524</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3575</td>\n",
       "      <td>2655</td>\n",
       "      <td>1</td>\n",
       "      <td>0.1621</td>\n",
       "      <td>9.048</td>\n",
       "      <td>2.157</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1835</td>\n",
       "      <td>35</td>\n",
       "      <td>1</td>\n",
       "      <td>0.1766</td>\n",
       "      <td>7.123</td>\n",
       "      <td>3.979</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      X     Y  facies_threshold_0.3  porosity  permeability  \\\n",
       "0   565  1485                     1    0.1184         6.170   \n",
       "1  2585  1185                     1    0.1566         6.275   \n",
       "2  2065  2865                     2    0.1920        92.297   \n",
       "3  3575  2655                     1    0.1621         9.048   \n",
       "4  1835    35                     1    0.1766         7.123   \n",
       "\n",
       "   acoustic_impedance  \n",
       "0               2.009  \n",
       "1               2.864  \n",
       "2               3.524  \n",
       "3               2.157  \n",
       "4               3.979  "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv(\"2D_MV_200wells.csv\")                      # read a .csv file in as a DataFrame\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": [
    "Let's extract one of the features, porosity, into a 1D ndarray and do our statistics on porosity.\n",
    "\n",
    "* then we can use NumPy's statistics methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "por = df['porosity'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's go through all the univariate statistics listed above one-by-one.\n",
    "\n",
    "#### Measures of Central Tendency\n",
    "\n",
    "Let's start with measures of central tendency.\n",
    "\n",
    "##### The Arithmetic Average / Mean\n",
    "\n",
    "\\begin{equation}\n",
    "\\overline{x} = \\frac{1}{n}\\sum^n_{i=1} x_i\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity average is 0.15.\n"
     ]
    }
   ],
   "source": [
    "por_average = np.average(por)\n",
    "print('Porosity average is ' + str(round(por_average,2)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Median\n",
    "\n",
    "\\begin{equation}\n",
    "P50_x = F^{-1}_{x}(0.50)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity median is 0.15.\n"
     ]
    }
   ],
   "source": [
    "por_median = np.median(por)\n",
    "print('Porosity median is ' + str(round(por_median,2)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Mode\n",
    "\n",
    "The most common value. To do this we should bin the data, like into histogram bins/bars.  To do this we will round the data to the 2nd decimal place.  We are assume bin boundaries, $0.01, 0.02,\\ldots, 0.30$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity mode is 0.14.\n"
     ]
    }
   ],
   "source": [
    "por_mode = stats.mode(np.round(por,2))\n",
    "print('Porosity mode is ' + str(round(por_mode,2)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Geometric Mean\n",
    "\n",
    "\\begin{equation}\n",
    "\\overline{x}_G = ( \\prod^n_{i=1} x_i )^{\\frac{1}{n}}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity geometric mean is 0.15.\n"
     ]
    }
   ],
   "source": [
    "por_geometric = scipy.stats.mstats.gmean(por)\n",
    "print('Porosity geometric mean is ' + str(round(por_geometric,2)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Harmonic Mean\n",
    "\n",
    "\\begin{equation}\n",
    "\\overline{x}_H = \\frac{n}{\\sum^n_{i=1} \\frac{1}{x_i}}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity harmonic mean is 0.14.\n"
     ]
    }
   ],
   "source": [
    "por_hmean = scipy.stats.mstats.hmean(por)\n",
    "print('Porosity harmonic mean is ' + str(round(por_hmean,2)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Power Law Average\n",
    "\n",
    "\\begin{equation}\n",
    "\\overline{x}_p = (\\frac{1}{n}\\sum^n_{i=1}{x_i^{p}})^\\frac{1}{p}\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity law mean for p = 1.0 is 0.15.\n"
     ]
    }
   ],
   "source": [
    "power = 1.0\n",
    "por_power = np.average(np.power(por,power))**(1/power)\n",
    "print('Porosity law mean for p = ' + str(power) + ' is ' + str(round(por_power,2)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Measures of Dispersion\n",
    "\n",
    "##### Population Variance\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma^2_{x} = \\frac{1}{n}\\sum^n_{i=1}(x_i - \\mu)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity population variance is 0.0011.\n"
     ]
    }
   ],
   "source": [
    "por_varp = stats.pvariance(por)\n",
    "print('Porosity population variance is ' + str(round(por_varp,4)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Sample Variance\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma^2_{x} = \\frac{1}{n-1}\\sum^n_{i=1}(x_i - \\overline{x})^2\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity sample variance is 0.0011.\n"
     ]
    }
   ],
   "source": [
    "por_var = stats.variance(por)\n",
    "print('Porosity sample variance is ' + str(round(por_var,4)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Population Standard Deviation\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_{x} = \\sqrt{ \\frac{1}{n}\\sum^n_{i=1}(x_i - \\mu)^2 }\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity sample variance is 0.0329.\n"
     ]
    }
   ],
   "source": [
    "por_stdp = stats.pstdev(por)\n",
    "print('Porosity sample variance is ' + str(round(por_stdp,4)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Sample Standard Deviation\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_{x} = \\sqrt{ \\frac{1}{n-1}\\sum^n_{i=1}(x_i - \\mu)^2 }\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity sample variance is 0.0329.\n"
     ]
    }
   ],
   "source": [
    "por_std = stats.stdev(por)\n",
    "print('Porosity sample variance is ' + str(round(por_std,4)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Range\n",
    "\n",
    "\\begin{equation}\n",
    "range_x = P100_x - P00_x\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity range is 0.17.\n"
     ]
    }
   ],
   "source": [
    "por_range = np.max(por) - np.min(por)\n",
    "print('Porosity range is ' + str(round(por_range,2)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Percentile\n",
    "\n",
    "\\begin{equation}\n",
    "P(p)_x = F^{-1}_{x}(p)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity 13th percentile is 0.11.\n"
     ]
    }
   ],
   "source": [
    "p_value = 13\n",
    "por_percentile = np.percentile(por,p_value)\n",
    "print('Porosity ' + str(int(p_value)) + 'th percentile is ' + str(round(por_percentile,2)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Inter Quartile Range\n",
    "\n",
    "\\begin{equation}\n",
    "IQR = P(0.75)_x - P(0.25)_x\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity interquartile range is 0.04.\n"
     ]
    }
   ],
   "source": [
    "por_iqr = scipy.stats.iqr(por)\n",
    "print('Porosity interquartile range is ' + str(round(por_iqr,2)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Tukey Test for Outliers\n",
    "\n",
    "Let's demonstrate the Tukey test for outliers based on the lower and upper fences.\n",
    "\n",
    "\\begin{equation}\n",
    "fence_{lower} = P_x(0.25) - 1.5 \\times [P_x(0.75) - P_x(0.25)]\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "fence_{upper} = P_x(0.75) + 1.5 \\times [P_x(0.75) - P_x(0.25)]\n",
    "\\end{equation}\n",
    "\n",
    "Then we declare samples values above the upper fence or below the lower fence as outliers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity outliers by Tukey test include [0.06726 0.05    0.06092].\n",
      "Porosity outlier indices by Tukey test are [110 152 198].\n"
     ]
    }
   ],
   "source": [
    "p25, p75 = np.percentile(por, [25, 75])\n",
    "lower_fence = p25 - por_iqr * 1.5\n",
    "upper_fence = p75 + por_iqr * 1.5\n",
    "por_outliers = por[np.where((por > upper_fence) | (por < lower_fence))[0]]\n",
    "print('Porosity outliers by Tukey test include ' + str(por_outliers) + '.')\n",
    "por_outliers_indices = np.where((por > upper_fence) | (por < lower_fence))[0]\n",
    "print('Porosity outlier indices by Tukey test are ' + str(por_outliers_indices) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Measures of Shape\n",
    "\n",
    "##### Pearson's Mode Skewness\n",
    "\n",
    "\\begin{equation}\n",
    "skew = \\frac{3 (\\overline{x} - P50_x)}{\\sigma_x}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity skew is -0.03.\n"
     ]
    }
   ],
   "source": [
    "por_skew = (por_average - por_median)/por_std\n",
    "print('Porosity skew is ' + str(round(por_skew,2)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Population Skew, 3rd Central Moment\n",
    "\n",
    "\\begin{equation}\n",
    "\\gamma_{x} = \\frac{1}{n}\\sum^n_{i=1}(x_i - \\mu)^3\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity 3rd cenral moment is -1.22e-05.\n"
     ]
    }
   ],
   "source": [
    "por_cm = scipy.stats.moment(por,moment=3)\n",
    "print('Porosity 3rd cenral moment is ' + str(round(por_cm,7)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Quartile Skew Coefficient\n",
    "\n",
    "\\begin{equation}\n",
    "QS = \\frac{(P75_x - P50_x) - (P50_x - P25_x)}{(P75_x - P25_x)}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity quartile skew coefficient is 0.14.\n"
     ]
    }
   ],
   "source": [
    "por_qs = ((np.percentile(por,75)-np.percentile(por,50))\n",
    "          -(np.percentile(por,50)-np.percentile(por,25))) /((np.percentile(por,75))-np.percentile(por,25))\n",
    "print('Porosity quartile skew coefficient is ' + str(round(por_qs,2)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot the Nonparametric CDF\n",
    "\n",
    "Let's demonstrate plotting a nonparametric cumulative distribution function (CDF) in Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAccAAAGWCAYAAAANLkjgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de3yU9Zn//9dFyIAmAkqCiBCLoOlhFUGFVouida0upbQ/2J5ExdV6KLZLa7d1t6672D19u+t3tV1atcYi9uCXNrXiYas9GKW13lUGNdU2ikEJp8ZEBCYCOV2/P+YeHMYcJgmTOb2fj0ceM3PPfd+55nbMxfW5Pwdzd0RERORtI7IdgIiISK5RchQREUmh5CgiIpJCyVFERCSFkqOIiEgKJUcREZEUSo4iecTMbjOzf8x2HINhZi+Y2bxsxyGSDiVHGVZm9qqZ/dnMypK2XWFmdVkMK+vMbKmZ/aa//dz9anf/+iDO/6qZ7TWzWHj9v2dm5YOLdnDc/X3uXhfG889m9v2hnM/MPmxmT5jZHjN73cweN7OPhu8tNbOu8PPGzGxT+JlPTDr+XWbmSfvEzOy5IX1IKRhKjpINI4G/zXYQfTGzkdmOIZWZlQzxFAvcvRyYBZwO3DCIGHLiupjZYuDHwGpgMnA0cCOwIGm334WfdyxwHrAXWG9mf5FyunHuXh7+zMh89JIPlBwlG/4T+LKZjevpTTM7w8yeNrNd4eMZSe/VmdnXzey3YcXwqJlVhO8lKoErzWybmW03s+uSjp1tZr8zszfD9/7HzCJJ77uZLTOzl4GXw223mlmTme02s/VmNjdp/382sx+b2ffDWOrN7EQz+3szaw6POz9p/7FmVhP+7q1m9i9mVmJm7wFuAz4QVi9vhvuvMrPvmNnDZtYGnBNu+5ekcy40s2fD+F4xswv6u/juvhX4X+AvwnNMMrO1ZvaGmW00s8+mfMafhJ9xN7DUzEaZ2S3hNd4WPh8V7l9hZg+G1/gNM1tnZiPC9141s/PCGP8B+GSiWjOzvzaz9Snfg+vM7Gc9fD8M+L/A1939Tnff5e7d7v64u382dX9373L3V9z9c8DjwD/3d41ElBwlG54B6oAvp75hZkcBDwHfBMYT/yP4kJmNT9rtM8BlwAQg0sN5zgFOAM4Hrjez88LtXcAXgQrgA8CHgM+lHPsxYA7w3vD108ApwFHAD4Efm9nopP0XAPcARwIbgEeI/391LHATcHvSvncDncB0YGYY3xXu/kfgasJKx92T/9HwGeBfgSOAg5pdzWw28crp74BxwFnAq/TDzKYAfxXGC/AjYAswCVgM/JuZfSjpkIXAT8Lf8QPga8D7w+syA5jN21XodeG5KolXc/8AHDRHpbv/HPg34P8lVWtrganhPxQSlhC/tqmqgSlhTAP1U2Buv3tJ0VNylGy5Efi8mVWmbJ8PvOzu97h7p7v/CPgTBzeXfc/dX3L3vcAa4n+kk61w9zZ3rwe+B3wawN3Xu/tT4XlfJZ64zk459t/d/Y3w3Lj79929NTzmZmAU8T/OCevc/RF37yTezFcJ/Ie7dwD3Au8ys3FmdjRwIbA8jK0Z+G/gU/1cp/vd/bdhZbQv5b3Lgbvc/Rfh+1vd/U99nOtnYVX6G+IV1L+FifKDwFfdfZ+7PwvcCVycdNzv3P1n4e/YC1wE3OTuze7+OrAiaf8O4BjgOHfvcPd1nsYEzu6+H/h/xBMiZvY+4F3Agz3snviH0vb+ztuDbcT/oZOsJax03zSzd/yDTYqTkqNkhbv/gfgfvutT3poEvJay7TXilVjCjqTnbwGpHUuaUo6dBBA2eT5oZjvCJsJ/I15F9nZsomnvj2ET75vE718lH/PnpOd7gRZ370p6TRjfcUApsD3xh5h4cp5A35r6eG8K8Eo/xyf7mLuPc/fj3P1zYaKbBLzh7nuS9ku93qkxpP43OnCNiTeZbwQeNbNGM0v979uXu4HPhM2mFwNrwqSZqjV8PGYA5044FngjZVtFeF3Guft/DeKcUoCUHCWb/gn4LAf/Id5GPJEkqwK2DuC8U1KO3RY+/w7xKvQEdx9DvMnPUo49UOWE9xe/CnwCODJs7tzVwzHpaAL2c/Af4jHu/r7U39tbPL2cc9ogYkm2DTjKzI5I2pZ6vVNjSP1vdOAau/sed7/O3Y8nXu1/KaWJtrdz4u5PAe3Emz0/Q89NqgANxD/7ot4+VB8+DqwbxHFSZJQcJWvcfSPxprQvJG1+GDjRzD5jZiPN7JPE7//11LzWm380s8PDprnLwt8B8ft2u4GYmb0buKaf8xxB/B7h68BIM7sRGDOAOA5w9+3Ao8DNZjbGzEaY2TQzSzTr/hmYnNxBKA01wGVm9qHwfMeGn2sgcTUBTwL/bmajzexk4s21P+jjsB8BN5hZpcU7Q90IfB/AzD5iZtPD6m838fu8XT2c48/Em5xT/watBv4H6HT3Hoe2hM20XyL+3/mypOv5QTO7I3X/sNPTVDP7FjCPeDOwSJ+UHCXbbgIOjHl091bgI8Q7drQCXwE+4u4tAzjn48Sb9n4F/Je7Pxpu/zLximQP8F3eTpq9eYR4r86XiDcd7qPvZs7+XEK8A9GLwE7iHUoSTYO/Bl4AdphZWp/V3X9PPPn/N/GK9nHeWXWn49PE7+9tA+4D/sndf9HH/v9CvFPV80A9EA23Qbwj1C+BGPA74NuJsY0pfhw+tppZNGn7PcR70fZWNQLg7j8BPgn8TRj3n8MY7k/a7QNmFiOepOuI/8Pm9PBetEifTIsdS6Ews3cBm4DSsIOM5BkzOwxoBma5+8vZjkeKlypHEckl1wBPKzFKtmUsOZrZXRYfCP2HXt43M/tmOOj4eTOblalYRCT3mdmrxGdOuq6fXUUyLpOV4yqgr9k6LiR+f+IE4EriPQlFBs3dX3V3U5NqfnL3d4XDTDb0v7dIZmUsObr7E7xzPFGyhcBqj3sKGGdmgxm3JCIickhlcxLhYzm459+WcNs7Zr0wsyuJV5eMHj361KqqqmEJsNB0d3czYoRuMw+Urtvg6doNTqavW3d394Hf0d/v6ezsxDo7KTGjyx0fOZKRI3Ni/vkevfTSSy3unjrz1oBl8xP2NJC6x66z7n4HcAdAdXW1NzQ0ZDKuglVXV8e8efOyHUbe0XUbPF27wcnkdQuCgNuXLWM+8UmMr1q5kjlz5hyy/bPNzFJn2BqUbCbHLRw8k8lk3p7JREREMqA+GuXsWIzTS0uJdXRQH432mezmzJkDK1dSH41y1axZOZ0YD6VsJse1wLVmdi/xVRB2hbOIiIhIhlgkwl2vvcZbXV3cW1LCJZH+J2WaM2dO0STFhIwlRzP7EfGpmirMbAvxeTRLAdz9NuLThP0V8ZlM3iI+04eIiGRQY0MDM9xpKylhhjuNuk3Vo4wlR3f/dD/vO7AsU79fRETeyYkv+vkB4AnCFa/lHdSNTESkiEyrrmYr8GhXF1vD1/JOudsfV0REDrnGhgbOGzuWSeXlnBiJ4O3t2Q4pJ6lyFBEpEkEQ8PwDD/BSLEbHjh1sAE6apZk7e6LKUUQkDwVBQH00ykkDGF5RH42ytKyMyTNncvfmzUxfsKDoeqGmS5WjiEieCYKAm5cupfPf/52bly4lCIK0jjtp1iweArbs28e+CRNYuGhRZgPNY6ocRUTyzNraWmZv3szC0lJira2sra1NqwIs1gH9g6HkKCKSo3prOnXgt8BUd34LvGcA5yzGAf2DoWZVEZEclJjT9MiaGm5ftuygptOFixbRWVXFr8aOpbOqSs2jGaDKUUQkByRXiRDvPDMfWHTMMbB9+0FzoM6ZM4cbVq2iPhrlUjWPZoQqRxGRLEutEtva2g50nqndvp2HeOeQizlz5nDFNdcoMWaIKkcRkSxbW1vL7OZmzq2qgn372PvWW+o8k2VKjiIiWVRTU8Pjq1dz1JtvEmtt5fdVVVx6+OGAOs9kk5KjiEiWBEHA3StWsHTPHkaNHEldWRnTFyygrKws26EVPd1zFBHJkvpolHlmPNbdzVudnTSPHq2epzlCyVFEJEssEmFdSwvHuXNrdzcnL1miZtQcoWZVEZEsaWxoYHZ5OUeXl3NxJELl5MnZDklCqhxFRLJAK2TkNlWOIiJZsLa2lvltbUybPp37du7UChk5RslRRCSDepofNVE1vtjaypmtrbRUVXGpOuLkFDWriohkSG/zoybWVbx+5kwax49X1ZiDlBxFRDKkPhrl7FiM0/fv5+xYjPpoFNC6ivlAzaoiIhlikQirmproAH4ALIlEAK2rmA+UHEVEMsTb25lbUUFrRwdzS0vx9vYD72lquNymZlURkQzZtGULj+3Ywbg33mBdSwsWVo6S+5QcRUQyIAgC1t1zD3PN2GbG3IqKgypHyW1KjiIiGZCYN/U1Mya5a5B/nlFyFBHJAM2bmt/UIUdE5BALgoDf/PKXzK2ooGrMGC7u7NS8qXlGlaOIyCGUGPh/4rPPsq6lhdKODl4uK1OTap5R5Sgicgitra1ldnMzn6yqAqDulFO4evlyNanmGSVHEZFDIAgC7q+tJVizhhdbW4m1trKhqorrlBjzkpKjiMgQJZpSRzc3c2FrKzO00kbe0z1HEZEhCIKA2265hbNjMS6tquK3wCs7d2rO1DynylFEZJAOdL5pa2NVUxNLgc6qKpoWLOCqRYtUNeYxJUcRkUFIVIzzYjEuPeEEIN755gbdYywISo4iIgOUWjECvFxerl6pBUTJUURkgDRco/ApOYqI9CMIAuqj0QMD+Z9/4AEN1yhwSo4iIn1INKHOB26vqeGYc89laVkZk2fO5O7NmzVco0BpKIeISB8STajnjh7NfMCBh4At+/ZpuEYBU+UoItKLIAgOakL9fVUV1y1aBIsWUR+NctWsWaoaC5SSo4hIL+qjURa7M27iRNbGYgc1oSopFjY1q4qI9MIiEVY1NfH6jh00xmJMq67OdkgyTFQ5ioj0wtvbmVtRQWtHB3NLS/H29myHJMNElaOISC82bdnCYzt2MO6NN1jX0oJFItkOSYaJkqOISA+CIGDdPfcw14xtZsytqFDlWESUHEVEelAfjTLPjNfMmOTOBjgwCYAUPiVHEZEeWCTCupYWjnPn1u5uTl6yRD1Ui4g65IiI9KCxoYHZ5eUcXV7OxZEIlZMnZzskGUaqHEVEUiQG/78Ui9GxY4eaVIuQKkcRkRRra2uZ39bGtOnTuW/nTs2fWoSUHEVEktTU1PD46tUc9eabnNXaSktVFZdq/tSio+QoIhIKgoC7V6xg6Z49jBo5krqyMlWNRUr3HEVEQonhG491d/NWZyfNo0dr1Y0ipeQoIhLS8A1JULOqiEgoMZdqWUcHi0tLqdLwjaKlylFEJKS5VCVByVFEBM2lKgdTchQRQXOpysGUHEVEUGccOVhGk6OZXWBmDWa20cyu7+H9sWb2gJk9Z2YvmNllmYxHRKQ3B+ZSPfZYLp42janqjFPUMpYczawEWAlcCLwX+LSZvTdlt2XAi+4+A5gH3GxmugMuIsNKc6lKqkwO5ZgNbHT3RgAzuxdYCLyYtI8DR5iZAeXAG0BnBmMSEXmH+miUxe6MmziRtbGYZsURzN0zc2KzxcAF7n5F+PpiYI67X5u0zxHAWuDdwBHAJ939oR7OdSVwJUBlZeWpa9asyUjMhS4Wi1FeXp7tMPKOrtvg5cu127ZtG3t27GA88IYZR02ZQkVFRdbiyZfrlovOOeec9e5+2lDPk8nK0XrYlpqJPww8C5wLTAN+YWbr3H33QQe53wHcAVBdXe3z5s079NEWgbq6OnTtBk7XbfDy4doFQcBt11/PB5ub2WaGT5jAiOuvZ/HixVmLKR+uW6HLZIecLcCUpNeTgW0p+1wG/NTjNgKbiFeRIiLDQkM4pCeZTI5PAyeY2dSwk82niDehJtsMfAjAzI4GqoHGDMYkInIQDeGQnmSsWdXdO83sWuARoAS4y91fMLOrw/dvA74OrDKzeuLNsF9195ZMxSQikurAEI7yci6ORKjUEA4hwxOPu/vDwMMp225Ler4NOD+TMYiI9CYxhGNkLEZlLMaGqiquU5OqoFU5RKSIaQiH9EbTx4lI0bJIhLtee41tTU1s3L2badXV2Q5JcoSSo4gUrcaGBma401ZSwgx3Ghsash2S5AglRxEpWg68ChwbPmZmShTJR0qOIlK0plVXsxV4tKuLreFrEVCHHBEpYo0NDZw3diyTyss5MRLR4sZygCpHESlKWolD+qLKUUSKkoZxSF+UHEWkqARBQH00yqYtW3iyqYmLiM9ZeYbuN0oSJUcRKRpBEHD7smXMB1Y3NzOvogLGjOH8zk7db5SD6J6jiBSN+miUs2MxTt+/n3lmbACOHDWKl8vKdL9RDqLKUUSKhkUirGpqogNYB5yxfDk7J0/mqlmzdL9RDqLkKCJFIQgCfvPLXzI3qSm1cvJkrrjmmmyHJjlIzaoiUvAS9xpPfPZZ1rW0UNrRoaZU6ZMqRxEpaEEQcNsttzAvFuPSE04AoO6UU7h6+XI1pUqvlBxFpGDV1NTws5tuYiawqiW+jvrL5eVKjNIvJUcRKUhBEHD3ihVc0trKaaWlUFFB3cyZSoySFiVHESlI9dEo88x4rLubrn372DB2LDcoMUqa1CFHRAqSRSKsa2nhOHdu7e7m5CVLlBglbaocRaQgJKaFOykcs+jt7ZxfVUXlyJFcHA7bEEmXkqOI5L0gCLh56VLObWvj5rIyrlu1ipNmzeL2sjJOBJ4cNYpzNWxDBkDJUUTy3traWmZv3szC0lJira2sra3lX7/xDVi5kvpoVDPgyIApOYpI3nPgia4uju3u5gl3/iLcPmfOHCVFGRR1yBGRvDetupqtwKNdXWwNX4sMhSpHEcl73t7O/KOPpqyjgymlpVp+SoZMlaOI5L3EsI3xu3axrqUFi0SyHZLkOVWOIpL3GhsamF1eTlt5OedHIqocZciUHEUkbwVBwP21tQRr1jAmFqMyFmNDVRXXadiGDJGSo4jkpcQyVKObm7mwtZUZ06dz386dTF+wQD1UZciUHEUkL9VHo5wdizGutJSari7Kd+5k34QJLFy0KNuhSQFQhxwRyUsWibCqqYnXd+xgKxCcdRZXrVypqlEOCVWOIpKXkjvhLI5EqDzzTCVGOWRUOYpI3gmCgOcfeICXYjE6duxgA3CSOuHIIaTKUUTyztraWua3tTFNnXAkQ1Q5ikheCYKAp9as4cEdO9jQ0EBLWZk64cghp8pRRPLK2tpaLnz9dd4TibC6vZ3DTz1VVaMccqocRSSvJFbg2N3Vxd4RI5g0aVK2Q5ICpOQoInlFK3DIcFCzqojklcaGBs4bO5ZJ5eWcqHlUJUNUOYpI3tAQDhkuqhxFJG/UR6MsdmfcxImsjcU0hEMyRpWjiOQNi0S467XX2NbUxMbdu3W/UTJGyVFE8kZjQwMz3GkrKWGGO40NDdkOSQqUkqOI5I2t27bR2N3NlJISmkpK8GwHJAWr3+RoZs+Y2TIzO3I4AhIR6UkQBLSuX8+eESP4aXs7uysrNTOOZEw6leOngEnA02Z2r5l92Mwsw3GJiBwkMZ/qDdXVjJ84kTmf+IQ640jG9Jsc3X2ju38NOBH4IXAXsNnMVpjZUZkOUEQkMYTjkdZWntu4UfOpSsaldc/RzE4Gbgb+E6gFFgO7gV9nLjQRkbhE1fi56dNpHD9eQzgk4/od52hm64E3gRrgenffH74VmNmZmQxORIpbEATcX1tLsGYNL7a2cmZrKy1VVVyqqlEyLJ1JAP7a3RuTN5jZVHff5O7/X4biEpEiFwQBty9bxujmZi5sbWWG1m6UYZROs+pP0twmInLI1EejnB2L8eHSUp7o6uKVnTvZN2GC7jXKsOi1cjSzdwPvA8aaWXKFOAYYnenARKS4WSTCqqYmLgK2AsFZZ3H15z+vqlGGRV/NqtXAR4BxwIKk7XuAz2YyKBGRxoYGZpeX01ZezuJIhMozz1RilGHTa3J09/uB+83sA+7+u2GMSUSKXE1NDY+vXs1Ru3ZRGYuxoaqK67T6hgyjvppVv+Lu3wA+Y2afTn3f3b+Q0chEpCgFQcDdK1awdM8eRo0cSV1ZmTrhyLDrq1n1j+HjM8MRiIgIxDvizDPjse5uznKnefRorlQnHBlmfTWrPhA+3j184YhIsbNIhHUtLXzAnVu7u/n4kiWqGmXY9dWs+gD0Pum9u380IxGJSFHz9nbmVlRQ1tHB4tJSqiZPznZIUoT6alb9r2GLQkQklKgcLwJ+DiyJRLIdkhShvppVHx/OQERE4OAhHOdHInh7e7ZDkiLU6ww5ZrYmfKw3s+eTfurN7Pl0Tm5mF5hZg5ltNLPre9lnnpk9a2YvmJkSskgRC4KAp9as4YWdO9m/dSsbgJM0hEOyoK9m1b8NHz8ymBObWQmwEvhLYAvx9SDXuvuLSfuMA74NXODum81swmB+l4gUhrW1tVz4+uu8JxJhdXs7h596qjrjSFb0Wjm6+/bw8TVgPzADOBnYH27rz2xgo7s3uns7cC+wMGWfzwA/dffN4e9qHvhHEJFCsXXbNuo6OogB7aWlTJo0KdshSZEy9147pMZ3MLsCuJH42o0GnA3c5O539XPcYuIV4RXh64uBOe5+bdI+twClxOdwPQK41d1X93CuK4ErASorK09ds2ZN2h9Q3haLxSgvL892GHlH123wBnLt2tra2P7qq3Tt308p0FlayrHHH09ZWVlmg8xB+s4N3jnnnLPe3U8b6nnSWbLq74CZ7t4KYGbjgSeBPpMj8USaKjUTjwROBT4EHAb8zsyecveXDjrI/Q7gDoDq6mqfN29eGmFLqrq6OnTtBk7XbfAGcu3u/M53KP3Wt6js6GBtLMbRl17KlZ//fGYDzFH6zmVfOslxC/HJxhP2AE1pHjcl6fVkYFsP+7S4exvQZmZPEG++fQkRKSrJq3A0AmdUV2c7JClifU0C8KXw6VYgMLP7iVd+C4Hfp3Hup4ETzGxqeI5PEb/HmOx+4H/MbCQQAeYA/z2gTyAiBUFDOCSX9FU5HhE+vhL+JNyfzondvdPMrgUeAUqAu9z9BTO7Onz/Nnf/o5n9HHge6AbudPc/DPRDiEh+SwzhOGznTsbv2sWGqVO1CodkVV+TAKwY6snd/WHg4ZRtt6W8/k/gP4f6u0Qkf2kIh+Safu85mlkl8BXiPUpHJ7a7+7kZjEtEiogDT3R1cawZe0eMYJqGcEiW9TrOMckPgD8BU4EVwKvE7yeKiBwS06qr2Qo82tXF1vC1SDal01t1vLvXmNnfhvOtPq5p3kTkUGpsaOC8sWOZVF7OieqMIzkgncqxI3zcbmbzzWwm8WEZIiJDFgQBzz/wAC/FYnTs2KH5VCUnpFM5/ouZjQWuA74FjAG+mNGoRKRorK2tZX5bG9OmT+e+nTuZvmCBOuNI1vWbHN39wfDpLuCczIYjIsUkMYTjuR07OKu5mZapU7l00aJshyXSf7OqmR1vZg+YWYuZNZvZ/WZ2/HAEJyKFLTGE45pIhPXd3RrCITkjnWbVHxJfeurj4etPAT8iPpuNiMiABUFAfTTK1m3beBGYWlJCe2mphnBIzkgnOZq735P0+vvhzDciIgMWBAG3L1vGfKC1rY23Kiv5VXc3nWVlLFSTquSIvuZWPSp8+piZXU98PUYHPgk8NAyxiUgBWltby+zmZs6tqgJg/YIFHD91KpfOmqUmVckZfVWO64knw8TSU1clvefA1zMVlIgUpsSwjRdbW4m1tvL7qiquW7RISVFyTl9zq04dzkBEpLAFQcBtt9zCYnfePXMmd2/erGEbkrPSmVu1FLgGOCvcVAfc7u4dvR4kIpKkpqaGn910EzOBVS0tLAX2TZige4ySs9LpkPMdoBT4dvj64nDbFZkKSkQKRxAE3L1iBZe0tnJaaSlUVFA3cyZXL1+uqlFyVjrJ8XR3n5H0+tdm9lymAhKRwlIfjTLPjMe6u+nat48NY8dygxKj5Lh05lbtMrNpiRfhBABdmQtJRAqJRSKsa2nhOHdu7e7m5CVLlBgl56VTOf4d8eEcjcR7rh4HXJbRqESkYHh7O3MrKijr6GBxaSlVk7VugeS+PpOjmY0A9gInANXEk+Of3H3/MMQmIgXAIhEeb27mU11dPFxSwiWRSLZDEulXn8nR3bvN7GZ3/wDw/DDFJCIFpLGhgRnutJWUMMOdxoaGbIck0q907jk+amaLzMz631VE5GBbt22jsbubKSUlNJWU4NkOSCQN6dxz/BJQBnSa2T7iTavu7mMyGpmI5L22tjZa169nz4gR/LS9nd2TJmlso+SFdNZzPGI4AhGRwrPrzTfjCxlXV3Pfzp1M+8Qn1FNV8kKvzapmNsHMbjGzB83s38xMlaKIpC0IAt56800eaW3luY0badGqG5JH+rrnuBpoA74FHAF8c1giEpGCUB+NciSwdOJE/lhernlUJa/01aw60d2/Fj5/xMyiwxGQiBQGi0Ro2bWL7qYmNpaUcEZ1dbZDEklbX5WjmdmRZnZUuLZjScprEZFeNTY0cDhoCIfkpb4qx7HE13RMHsKRqB4dOD5TQYlI/tu6bRtTJ05kSkkJv3XnPdkOSGQA+lrP8V3DGIeIFJAgCGhdv56uWbM0hEPyUjrjHEVEBqQ+GmWxO+MiEcqPPFJDOCTvpDNDjojIgFgkwqqmJjo6OmiMxZimzjiSZ5QcReSQS6zE0VlSwtyKCry9PdshiQxIWsnRzD5oZpeFzyvNbGpmwxKRfLZpyxYe27GDkZ2drGtpwbQSh+SZfpOjmf0T8FXg78NNpcD3MxmUiOSvIAhYd889zDWjHVQ5Sl5Kp3L8OPBR4rPl4O7biM+YIyLyDmtra/nA/v1sGTmSiBkbgJNmzcp2WCIDkk5v1XZ3dzNzADMry3BMIpKngiDg+QceYGQsxpSuLlpKSvjYjTeqp6rknXSS4xozux0YZ2afBf4G+G5mwxKRfHRgCMfEiayNxThi/HguuuiibIclMmDpLFn1X2b2l8BuoBq40d1/kfHIRCTvJIZwXAQ0AheOGpXtkEQGpd/kaGZfBH6shCgi/WlsaGB2eTlt5eWcH4mAe7ZDEhmUdDrkjCG+Ksc6M2xtLE0AABtMSURBVFtmZkdnOigRyT81NTU8vno1f9y1i44dO9gAHHb44dkOS2RQ+k2O7r7C3d8HLAMmAY+b2S8zHpmI5I0gCLh7xQqW7tnDX48ceWD9xrIy9d+T/DSQGXKagR1AKzAhM+GISD6qj0aZZ8Zj3d281dlJ8+jRmmhc8lo6kwBcY2Z1wK+ACuCz7n5ypgMTkfxhkQjrWlo4zp1bu7s5eckSDd+QvJbOUI7jgOXu/mymgxGR/JSYS7Wso4PFpaVUTZ6c7ZBEhqTXytHMxoRPvwFsNrOjkn+GJzwRyQeJuVTHvfGG5lKVgtBX5fhD4CPAesABS3rPgeMzGJeI5InkuVS3mWkuVSkIvSZHd/9I+KgVOESkV4nOOK+YcZY7DwLzNZeq5Ll0OuT8Kp1tIlKc1BlHClGvlaOZjQYOByrM7EjeblYdQ3y8o4gUuSAI+M0vf8ncigqqxozh4s5OKtUZRwpAX/ccrwKWE0+E63k7Oe4GVmY4LhHJcUEQcPuyZZzY1sajLS0sPewwXi4v51w1qUoB6Oue463ArWb2eXf/1jDGJCJ5YG1tLbObm/lkVRUAdaecwtXLl6tJVQpCOqtyfMvM/gJ4LzA6afvqTAYmIrkpCALur60lWLOGF1tbibW2sqGqiuuUGKWApLMqxz8B84gnx4eBC4HfAEqOIkUm0ZQ6urmZC1tbmTF9Ovft3Mn0BQuUGKWgpDO36mLgQ8AOd78MmAFokTaRIlQfjXJ2LMaHS0t5oquLV3buZN+ECZpHVQpOOslxr7t3A53hrDnNaAIAkaKUWMz49R072AoEZ53FVStXqmqUgpPO3KrPmNk44LvEe63GgN9nNCoRyUmJOVRbOzqYX1pK1ZlnKjFKQUqnQ87nwqe3mdnPgTHu/nxmwxKRXJQY8H8R8HNgieZQlQLV1yQAvQ5WMrNZ7h7NTEgikqsaGxqYXV5OW3k550cimkNVClZflePNfbznwLmHOBYRyWFBEPDUmjUctnMn43ftYsPUqVynAf9SoPqaBOCc4QxERHLb2tpaLnz9dd4TibC6vZ3DTz1V9xulYKUzzvGSnrZrEgCR4uLAE11dHGvG3hEjmDZJUyxL4Uqnt+rpSc9HEx/zGEWTAIgUlWnV1TwCPNrVxdYRI1hYXZ3tkEQyJp3eqp9Pfm1mY4F70jm5mV0A3AqUAHe6+3/0st/pwFPAJ939J+mcW0SGT2L1jflHH03VmDGc2NmpzjhS0NKZBCDVW8AJ/e1kZiXEV++4kPjUc582s/f2st//AR4ZRCwikmEHVt949lnWtbRQ2tHBy2VlnKTOOFLA0rnn+ADx2w0QT6bvBdakce7ZwEZ3bwzPcy+wEHgxZb/PA7Uc3HwrIjkiMWXcOaWl7K2ooG7mTK2+IQXP3L3vHczOTnrZCbzm7lv6PbHZYuACd78ifH0xMMfdr03a51jgh8SHhdQAD/bUrGpmVwJXAlRWVp66Zk06uVlSxWIxysvLsx1G3in269bS0kLr5s0cBbwBjK+qoqKiIq1ji/3aDZau2+Cdc8456939tKGeJ517jo8DhPOqjgyfH+Xub/RzqPWwLTUT3wJ81d27zHra/UAMdwB3AFRXV/u8efP6C1t6UFdXh67dwBX7dfvaV77Cnm9+kyozNrtzxBe+wL9+4xtpHVvs126wdN2yL51m1SuBrwN7gW7iSc/pf/LxLcCUpNeTgW0p+5wG3Bsmxgrgr8ys091/llb0IpJxW7dto6W7mzMjEX7rznuyHZDIMEhnKMffAe9z95YBnvtp4AQzmwpsBT4FfCZ5B3efmnhuZquIN6sqMYrkiCAIaF2/nj0jRvDT9nZ2T5qk5amkKKSTHF8h3kN1QNy908yuJd4LtQS4y91fMLOrw/dvG+g5RWR41UejLHZn3LHHsjYWY9onPqGOOFIU0kmOfw88aWYBsD+x0d2/0N+B7v4w8HDKth6TorsvTSMWERlGifUbLwIagTM08F+KRDrJ8Xbg10A98XuOIlIkktdvnFtaqoH/UjTSSY6d7v6ljEciIjkjCALqo1E2bdnCk1q/UYpQOsnxsbDH6gMc3Kza31AOEclDiRlx5gOrm5uZV1EBY8ZwvqaMkyKSTnJM9DD9+6Rt6QzlEJE8lJgR5/TSUuaZsQGYMWoUT44axbmaMk6KRDqTAEztbx8RKRyJTjgdwDrgjOXL2Tl5MlfNmqWeqlI0tJ6jiBwktRNO1eTJXHHNNdkOS2RYaT1HETnIpi1beGLHDpaY8fOSEnXCkaKU0fUcRSS/BEHAunvuYa4Z28yYW1GhTjhSlDK2nqOI5J/6aJR5ZrxmxiR3NoDWbZSilMn1HEUkz1gkwrqWFj7gzq3d3Xx8yRJ1wpGilM49x/9Kep72eo4ikn8aGxqYXV7O0eXlXByJUDl5crZDEsmKXpOjmU0Hjk6s55i0fa6ZjXL3VzIenYgMmyAIeP6BBxgZi1EZi7Ghqorr1KQqRaqve463AHt62L43fE9ECkhiBY6lEyfyx/Jypi9YoCZVKVp9Nau+y92fT93o7s+Y2bsyFpGIZMWmLVt4YtMmlpjRWFKiFTikqPVVOY7u473DDnUgIpI9GsIhcrC+kuPTZvbZ1I1mdjmwPnMhichw0xAOkYP11ay6HLjPzC7i7WR4GhABPp7pwERk+GgIh8jBek2O7v5n4AwzOwf4i3DzQ+7+62GJTESGjYZwiBwsnenjHgMeG4ZYRCQLNIRD5J3SmQRARApYYgjHuIkTWRuLaQiHCIObW1VECohFItz12mtsa2pi4+7dTNMQDhElR5Fi19jQwAx32kpKmOFOY0NDtkMSyTolR5Ei58CrwLHho/e1s0iRUHIUKXLTqqvZCjza1cXW8LVIsVOHHJEi5+3tzD/6aMo6OphSWqqZcURQ5ShS9BITAIzftYt1LS1YJJLtkESyTpWjSJHz9nbmVlTQ2tHBXFWOIoAqR5Git2nLFh7bsYNxb7yhylEkpOQoUsS0GodIz5QcRYrY2tpaPrB/P1tGjmTKiBFajUMkpOQoUqQSc6q+FIsxrrOT7x5xBB+78UZNHSeCOuSIFK3UOVU/dMklXH755dkOSyQnqHIUKVIWibCqqYnXd+ygMRbT4H+RJKocRYpUYg3HtvJyzo9E1BFHJIkqR5EiFAQBT61Zwws7d7J/61Z1xBFJocpRpAitra3lwtdf5z2RCKvb2zn81FPVEUckiZKjSJEJgoDnnn2WEV1dTB01ivbSUqZNmpTtsERyipKjSBEJgoDbly3jzLY2fgI8PHo0nUcdxcJFi7IdmkhOUXIUKSL10Shnx2KcU1rK3qOPpumMM7hh+XI1qYqkUHIUKSKbtmzhiU2b2GfGupISlpx3nhKjSA/UW1WkSGgeVZH0KTmKFAnNoyqSPiVHkSKgeVRFBkb3HEWKgOZRFRkYVY4iRUDzqIoMjCpHkSLg7e3MraigtaODuaWl6ogj0g9VjiJFwCIRHm9uZkxLC483N2ORSLZDEslpSo4iRaCxoYEZ7rSVlDDDncaGhmyHJJLTlBxFioADrwLHho+ezWBE8oCSo0gRmFZdzVbg0a4utoavRaR36pAjUgQaGxo4b+xYJpWXc6IWNhbplypHkQKXPAFAx44dmhlHJA1KjiIFLAgCbrvlFha7c/3MmTSOH8/0BQs0M45IP9SsKlKgampq+NlNNzETWNXSwlJg34QJWrtRJA1KjiIFKAgC7l6xgktaWzmttBQqKqibOZOrtXajSFqUHEUKUGIFjifMoKODDaBFjUUGQPccRQqMVuAQGTolR5ECktoBp3PCBK3AITIIalYVKRBBEHD7smWc2NbGqqYmdcARGQIlR5E8FwQB9dEomzZtYj6waPp0AOpOOUUdcEQGKaPJ0cwuAG4FSoA73f0/Ut6/CPhq+DIGXOPuz2UyJpFCkhiusfiww3jejJcBtm/n5bIyJUaRIchYcjSzEmAl8JfAFuBpM1vr7i8m7bYJONvdd5rZhcAdgP5vFklD8nCNGaWlMHEiDR/9KDunTuWqWbOUGEWGIJOV42xgo7s3ApjZvcBC4EBydPcnk/Z/CpicwXhECkai4808M54sLYWODu7fu5cbFi1SUhQ5BMw9M4vXmNli4AJ3vyJ8fTEwx92v7WX/LwPvTuyf8t6VwJUAlZWVp65ZsyYjMRe6WCxGeXl5tsPIO7l23dra2mjZvJlR3d3sbm+nfORI9gDjjzmGioqKbId3kFy7dvlC123wzjnnnPXuftpQz5PJytF62NZjJjazc4DLgQ/29L6730G8yZXq6mqfN2/eIQqxuNTV1aFrN3C5dt2+9pWvMOXee/loVRW3tbbyuxzueJNr1y5f6LplXyaT4xZgStLrycC21J3M7GTgTuBCd2/NYDwieS8xwP/F1lZira1sqKriuhxNjCL5LJOTADwNnGBmU80sAnwKWJu8g5lVAT8FLnb3lzIYi0hBWFtby/y2Nj43fbpW2BDJoIxVju7eaWbXAo8QH8pxl7u/YGZXh+/fBtwIjAe+bWYAnYeirVikENXU1PD46tUc9eabnNXaSktVFZdqgL9IRmR0nKO7Pww8nLLttqTnVwDv6IAjIgdLDNtYumcPo0aOpK6sTFWjSAZpblWRHJc6bGO/O82jR2taOJEM0vRxIjkseb7UR1tamFtRwWp3LtUqGyIZpeQokoN6my/1pVNO4T/VO1Uk45QcRXJMEATcvHQp57a18dSIEbx82GGaL1VkmCk5iuSYtbW1zN68mYWlpcQ6OvjDxz/OzjPP1HypIsNIyVEkx2zdto3nOjqYYsZvgfdMmsQV11yT7bBEiop6q4rkkCAIaF2/nj0jRvDT9nZ2V1aqV6pIFqhyFMkByR1wlpaVMfm007h782amfeITakoVyQIlR5EsSwzXmA/c39bGy8AngX0TJqhqFMkSJUeRLHnHcI1jjoHt21l/7rlasFgky5QcRbKgt+EaDwFXacFikaxTchTJAg3XEMltSo4iWaDhGiK5TUM5RIaZhmuI5D5VjiIZluh4c1LYZFofjWq4hkiOU+UokkGJYRpH1tRw+7JlBEHASbNm8RCwZd8+DdcQyVGqHEUyqD4a5exYjNPDjjf10Wj83uLKldRHo+qAI5KjlBxFMmjTli08sWkT+8y4t6SEJZEIAHPmzFFSFMlhalYVyZAgCFh3zz3MNWObGXMrKvD29myHJSJpUHIUyZC1tbV8YP9+towcyZQRI9gAnDRrVrbDEpE0KDmKZEAQBDz/wAO8FIsxrrOT7x5xBB+78UY1pYrkCSVHkUMsCAJuu+UWFrtz/cyZdE6YwIcuuYTLL78826GJSJrUIUfkEEoM3TixrY1VTU0sRatriOQjJUeRQyRRMc6Lxbj0hBMAqDvlFK5evlzNqSJ5RslR5BCoqanhZzfdxExgVUsLAC+XlysxiuQpJUeRIQqCgLtXrOCS1lZOKy2FigrqZs5UYhTJY0qOIkNUH40yz4zHurvp2rePDWPHcoMSo0heU29VkSGySIR1LS0c586t3d2cvGSJEqNInlPlKDJE3t7O3IoKyjo6WFxaStXkydkOSUSGSJWjyBAlKsfxu3axrqUFC+dPFZH8pcpRZIBS12f09nbOr6qCkSM5v7NT86eKFAAlR5EBSAzZWHzYYdxeXg4rV3LSrFncXlbGicCTo0ZxruZPFcl7So4iaUoesjGjtBQmTtT6jCIFSslRJE2pQzYe3LuXG8IqUeszihQWdcgRSZOGbIgUD1WOIr1IdLw5pqoK4EDHm8qRI7m4s5NKDdkQKVhKjiKh5F6oALcvW8Z8oOXyywmCQB1vRIqIkqMIby81NR+4vaaGY849l/nAomOO4WegjjciRUbJUYR48js7FuP00lJiHR00AA8BbN9ODA5Uk+p4I1IclBxFiHe2WdXURAfwA2BJdTULFy06cM9RCVGkuCg5StFKvsfY0yw3iSqxrq4u26GKyDDTUA4pSol7jEfW1HD7smVYJMLLZWUcOWoUL5eVHWhGFZHipMpRikZypVgfjR7ocMP27exsb+cqdbYRkZCSoxSU1EnBk7cn90Y985prDnS4eQgOJEQlRREBJUcpIKkJkJUrDyQ7VYoiMhBKjpL3EtXipk2bDkqA9dHogaR30qxZ8YSpSlFE0qDkKHktuVq8v62Nl+GgBJgwZ84cDeAXkbQpOUre6atjzfpzz2Xn1Kk9JkBViiKSLiVHySv9dqxZtEgJUESGTMlR8sra2lpmNzdzblUV7NunjjUikhFKjpI3giDgqTVreG7HDnY3N/PM1Klcp441IpIBSo6SN9bW1nLh66/znkiE1e3tHH7qqUqKIpIRSo6S8xIdcLZu28aLwNSSEtpLS5k2aVK2QxORAqXkKDktuQNOa1sbb1VW8qvubjrLyli4aFG2wxORAqXkKDmpt4H96xcs4PipU7lUnW9EJIOUHCXn9DmwX0M1RGQYKDlKVvQ2QTi8cx7Uvgb2i4hkgpKjHHJ9Jb7E+71NEA49zIOqalFEhpmSoxxS/SU+eGdlmDxBOGgeVBHJvhHZDkDyRxAE3Pmd7xAEQa/7JCe++eHrVCfNmsVDQG1YGZ6UNEF4wpw5c7jimmuUGEUkKzKaHM3sAjNrMLONZnZ9D++bmX0zfP95M3vnX0nJCYmK8MiaGm5ftqzXBJlu4rtq5Up2Xn45V/VQWYqIZFvGmlXNrARYCfwlsAV42szWuvuLSbtdCJwQ/swBvhM+yiEWBAEtr79OEASDSkb9NYUmpNskqinfRCSXZbJynA1sdPdGd28H7gUWpuyzEFjtcU8B48zsmAzGVJQSVd/IlpY+q76+pFMRJqhJVETyXSY75BwLNCW93sI7q8Ke9jkW2J68k5ldCVwZvtxvZn84tKEWNoPKiVDx0Pr1VgK+6v3vb3F4fRCnKlsFhzu89b33v7/tUMeZwyqAlmwHkad07QZH123wqg/FSTKZHK2HbT6IfXD3O4A7AMzsGXc/bejhFR9du8HRdRs8XbvB0XUbPDN75lCcJ5PNqluAKUmvJwPbBrGPiIjIsMpkcnwaOMHMpppZBPgUsDZln7XAJWGv1fcDu9x9e+qJREREhlPGmlXdvdPMrgUeAUqAu9z9BTO7Onz/NuBh4K+AjcBbwGVpnPqODIVcDHTtBkfXbfB07QZH123wDsm1M/d33OITEREpapohR0REJIWSo4iISIqcSo5DmW6uv2ML2RCv26tmVm9mzx6qLtD5JI1r924z+52Z7TezLw/k2EI2xOum71zf1+6i8P/T583sSTObke6xhWyI123g3zl3z4kf4p12XgGOByLAc8B7U/b5K+B/iY+PfD8QpHtsof4M5bqF770KVGT7c+TwtZsAnA78K/DlgRxbqD9DuW7he/rO9X3tzgCODJ9fqL9zQ7tu4esBf+dyqXIcynRz6RxbqDRN3+D1e+3cvdndnwY6BnpsARvKdSt26Vy7J919Z/jyKeLjv9M6toAN5boNSi4lx96mkktnn3SOLVRDuW4Qn5HoUTNbH07TV0yG8r3Rd+5tA/3s+s69rb9rdznxVp/BHFtIhnLdYBDfuVxa7Hgo082lNQ1dgRrqNH1nuvs2M5sA/MLM/uTuTxzSCHPXUL43+s4dbCCfXd+5g/V47czsHOJ/5D840GML0FCuGwziO5dLleNQppsr5mnohjRNn7snHpuB+4g3XxSLoXxv9J1724A+u75z/V87MzsZuBNY6O6tAzm2QA3lug3qO5dLyXEo082lc2yhGvR1M7MyMzsCwMzKgPOBYlrxZCjfG33nBvHZ9Z3r/9qZWRXwU+Bid39pIMcWsEFft8F+53KmWdWHMN1cb8dm4WMMu6FcN+Bo4D4zg/h34Yfu/vNh/ghZk861M7OJwDPAGKDbzJYT7yW3W9+5gV834ksx6TvX9/+vNwLjgW+H16nT3U/T37nBXTcG+XdO08eJiIikyKVmVRERkZyg5CgiIpJCyVFERCSFkqOIiEgKJUcREZEUSo4iKcysK5y9/w9m9mMzOzyDv+smMzsvfL58oL8rHLv6azMbE77+gpn90cx+cAhiW2pmk5Je32lm7x3kua41s8v631MkN2goh0gKM4u5e3n4/AfAenf/v2kcZ8T/n+oe5O99FTjN3VsGcMx84Dx3/2L4+k/Ahe6+KWW/ke7eOcB46oivqDHkZaXCpP9bd5851HOJDAdVjiJ9WwdMBzCzL4XV5B/CQe2Y2bvCSu3bQBSYYmaftvjacX8ws/8T7ldiZqvCbfVmlkhmq8xssZl9AZgEPGZmj5nZ5Wb234kgzOyzZtZTgr4IuD/c5zbiS/qsNbMvmtk/m9kdZvYosDqMdZ2ZRcOfM5LO/5UwrufM7D/MbDFwGvCDsIo+zMzqzOy0cP93fMZwe8zM/jU8z1NmdjSAu78FvGpmxTRVnOSzoa6zpR/9FNoPEAsfRxJPPNcApwL1QBlQDrwAzATeBXQD7w+PmQRsBirD438NfCw8/hdJv2Nc+LgKWBw+f5Vwzbnw97wClIavnwRO6iHW14Ajkl4nn+OfgfXAYeHrw4HR4fMTgGfC5xeG5z88fH1U+FhHvJIl+XVvnzHcx4EF4fNvADckHf814Lps//fVj37S+VHlKPJOh5nZs8SnP9sM1BCf4f8+d29z9xjxORznhvu/5vF1MiG+wG+du7/u8WbMHwBnAY3A8Wb2LTO7ANjdVwDu3kY86XzEzN5NPEnW97DrUe6+p49TrXX3veHzUuC7ZlYP/Jj4dG4A5wHf83h1h7u/0VdsfXxGgHbgwfD5euL/eEhoJp5YRXJezsytKpJD9rr7KckbwvuJvWlL3rWnHdx9p5nNAD4MLAM+AfxNP3HcCfwD8Cfge73s02lmI7z3+5zJsX0R+DMwg/gtlX1JMQ+k80Ff16LD3RPn6uLgvzGjgb3vPEQk96hyFEnPE8DHzOzwcGb/jxO/H5kqAM42swozKwE+DTxuZhXACHevBf4RmNXDsXuAIxIv3D0gvkzPZ4Af9RJXA/H7jOkYC2wPE+nFxCdwBngU+JtET1kzO6qnePr7jGn8/hMprhU4JI8pOYqkwd2jxO8P/p54crjT3Tf0sN924O+Bx4DngKi730981fK6sLl2VbhPqjuA/zWzx5K2rSHey3NnL6E9BMxL82N8G7jUzJ4inqjawph/Tnz5n2fC+L4c7r8KuC3RISeNz9ifM4FfphmrSFZpKIdIDjOzB4H/dvdf9fL+McBqd//L4Y1sYMxsJvAld78427GIpEOVo0gOMrNxZvYS8fufPSZGOFDFfTcxCUAOqyDenCySF1Q5ioiIpFDlKCIikkLJUUREJIWSo4iISAolRxERkRRKjiIiIin+f8IWYQToxtrKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# sort the data:\n",
    "por_sort = np.sort(por)\n",
    "\n",
    "# calculate the cumulative probabilities assuming known tails\n",
    "p = np.arange(len(por)) / (len(por) - 1)\n",
    "\n",
    "# plot the cumulative probabilities vs. the sorted porosity values\n",
    "plt.subplot(122)\n",
    "plt.scatter(por_sort, p, c = 'red', edgecolors = 'black', s = 10, alpha = 0.7)\n",
    "plt.xlabel('Porosity (fraction)'); plt.ylabel('Cumulative Probability'); plt.grid(); \n",
    "plt.title('Nonparametric Porosity CDF')\n",
    "plt.ylim([0,1]); plt.xlim([0,0.25])\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.2, hspace=0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fit a Gaussian Distribution\n",
    "\n",
    "Let's fit a Gaussian distribution\n",
    "\n",
    "* we get fancy with Maximuum Likelihood Estimation (MLE) for the Gaussian parametric distribution fit mean and standard deviation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAccAAAGWCAYAAAANLkjgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdeXhU5dnH8e+dZQgQwpawE3bijoCAGwWstW64FGqtG/iKK63SYqu2LkVxaalVa1EW2VsXJKIgVqyVAFYZwYCgUhDCvoWEJUyAbPO8f8xAQyRhWCaTSX6f65prZs45c3LPAXJzP+dZzDmHiIiI/E9MpAMQERGpapQcRUREylByFBERKUPJUUREpAwlRxERkTKUHEVERMpQchSJImY2xswei3QcJ8LMvjGzvpGOQyQUSo5SqcxsvZntMLO6pbYNMbOMCIYVcWY22Mw+PdZxzrl7nHNPncD515vZATPzBa//JDNLPLFoT4xz7kznXEYwnj+Y2d9P5nxm9mMzW2Bm+8xsp5nNN7NrgvsGm1lJ8Pv6zGxd8Dt3LvX5tmbmSh3jM7OvTupLSrWh5CiREAc8EOkgKmJmcZGOoSwziz3JU/R3ziUC3YAewKMnEEOVuC5mNhB4G5gKtAKaAo8D/Usd9nnw+9YHLgUOAF+a2VllTtfAOZcYfHQJf/QSDZQcJRJGAQ+aWYOj7TSzC81ssZntDT5fWGpfhpk9ZWb/CVYMH5lZcnDfoUrgLjPbambbzGx4qc/2NLPPzWxPcN/fzMxTar8zs6Fm9h3wXXDbS2a2yczyzOxLM+td6vg/mNnbZvb3YCwrzKyzmT1iZtnBz11W6vj6ZjYh+LO3mNlIM4s1s9OBMcAFweplT/D4yWb2qpl9YGb5QL/gtpGlznmtmS0LxrfWzC4/1sV3zm0B/gmcFTxHCzObZWa7zGyNmd1Z5jvOCH7HPGCwmdUysxeD13hr8HWt4PHJZvZ+8BrvMrOFZhYT3LfezC4Nxvg74GeHqjUz+6mZfVnm78FwM3v3KH8/DPgL8JRz7jXn3F7nnN85N985d2fZ451zJc65tc65+4D5wB+OdY1ElBwlEpYAGcCDZXeYWSNgDvBXoDGBX4JzzKxxqcNuAm4HmgCeo5ynH9AJuAx42MwuDW4vAX4FJAMXAD8E7ivz2euAXsAZwfeLgXOBRsDrwNtmllDq+P7ANKAhsBSYS+DfVUvgSWBsqWOnAMVAR6BrML4hzrmVwD0EKx3nXOn/NNwEPA3UA45odjWzngQqp98ADYAfAOs5BjNrDVwZjBfgDWAz0AIYCDxjZj8s9ZFrgRnBn/EP4PfA+cHr0gXoyf+q0OHBc6UQqOZ+BxwxR6Vz7kPgGeCtUtXaLKBd8D8Kh9xC4NqWlQa0DsZ0vN4Beh/zKKnxlBwlUh4HfmlmKWW2XwV855yb5pwrds69AfyXI5vLJjnnVjvnDgDTCfySLm2Ecy7fObcCmAT8HMA596VzblHwvOsJJK4+ZT77rHNuV/DcOOf+7pzLDX7meaAWgV/Ohyx0zs11zhUTaOZLAZ5zzhUBbwJtzayBmTUFrgCGBWPLBl4AbjzGdXrPOfefYGV0sMy+O4CJzrl/Bfdvcc79t4JzvRusSj8lUEE9E0yUFwMPOecOOueWAa8Bt5b63OfOuXeDP+MAcDPwpHMu2zm3ExhR6vgioDnQxjlX5Jxb6EKYwNk5VwC8RSAhYmZnAm2B949y+KH/KG071nmPYiuB/+iUlhOsdPeY2ff+wyY1k5KjRIRz7msCv/geLrOrBbChzLYNBCqxQ7aXer0fKNuxZFOZz7YACDZ5vm9m24NNhM8QqCLL++yhpr2VwSbePQTuX5X+zI5Srw8AOc65klLvCcbXBogHth36RUwgOTehYpsq2NcaWHuMz5d2nXOugXOujXPuvmCiawHscs7tK3Vc2etdNoayf0aHrzGBJvM1wEdmlmVmZf98KzIFuCnYbHorMD2YNMvKDT43P45zH9IS2FVmW3LwujRwzv35BM4p1ZCSo0TSE8CdHPmLeCuBRFJaKrDlOM7busxntwZfv0qgCu3knEsi0ORnZT57uMoJ3l98CLgBaBhs7tx7lM+EYhNQwJG/iJOcc2eW/bnlxVPOOTucQCylbQUamVm9UtvKXu+yMZT9Mzp8jZ1z+5xzw51z7QlU+78u00Rb3jlxzi0CCgk0e97E0ZtUAVYR+O4DyvtSFbgeWHgCn5MaRslRIsY5t4ZAU9r9pTZ/AHQ2s5vMLM7Mfkbg/t/RmtfK85iZ1Qk2zd0e/BkQuG+XB/jM7DTg3mOcpx6Be4Q7gTgzexxIOo44DnPObQM+Ap43syQzizGzDmZ2qFl3B9CqdAehEEwAbjezHwbP1zL4vY4nrk3AZ8CzZpZgZucQaK79RwUfewN41MxSLNAZ6nHg7wBmdrWZdQxWf3kE7vOWHOUcOwg0OZf9HTQV+BtQ7Jw76tCWYDPtrwn8Od9e6npebGbjyh4f7PTUzsxeBvoSaAYWqZCSo0Tak8DhMY/OuVzgagIdO3KB3wJXO+dyjuOc8wk07f0b+LNz7qPg9gcJVCT7gPH8L2mWZy6BXp2rCTQdHqTiZs5juY1AB6Jvgd0EOpQcahr8BPgG2G5mIX1X59wXBJL/CwQq2vl8v+oOxc8J3N/bCswEnnDO/auC40cS6FS1HFgBZAa3QaAj1MeAD/gceOXQ2MYy3g4+55pZZqnt0wj0oi2vagTAOTcD+Bnwf8G4dwRjeK/UYReYmY9Aks4g8B+bHsF70SIVMi12LNWFmbUF1gHxwQ4yEmXMrDaQDXRzzn0X6Xik5lLlKCJVyb3AYiVGibSwJUczm2iBgdBfl7PfzOyvwUHHy82sW7hiEZGqz8zWE5g5afgxDhUJu3BWjpOBimbruILA/YlOwF0EehKKnDDn3HrnnKlJNTo559oGh5ksPfbRIuEVtuTonFvA98cTlXYtMNUFLAIamNmJjFsSERE5pSI5iXBLjuz5tzm47XuzXpjZXQSqSxISErqnpqZWSoDVjd/vJyZGt5mPl67bidO1OzHhvm5+v//wzzjWzykuLsaKi4k1o8Q5XFwccXEnlzpKSkqOePj9/iNeH+3hnDvidQVynHNlZ946bpFMjkcbSH3Ub+ycGweMA0hLS3OrVq0KZ1zVVkZGBn379o10GFFH1+3E6dqdmHBeN6/Xy9ihQ7mKwCTGd48eTa9evU7J8SUlJWzevJn169ezYcMGNmzYwObNm9m6dSvbtm1j69atZGdnU1JytKGvYGYkJSVRv359kpKSSEpKIjExkbp16x5+1KlTh9q1ax9+JCQkUKtWrcOPG2+8sewMWyckkslxM0fOZNKK/81kIiIiYbAiM5M+Ph894uPxFRWxIjOzwuTYq1cvGD2aFZmZ3N2tG7169SI/P59vvvmGFStWsHLlSlavXs13333H2rVrKSoqOuLzKSkptGjRghYtWnDuuefStGlTUlJSDj8aN25Mw4YNadSoEUlJSSddMd9447GmKw5NJJPjLOAXZvYmgVUQ9gZnERERkTAxj4eJGzawv6SEN2Njuc1T8aRMPp+PgoIC9uTn8/zzz5OZmUlWVtbhps2EhAQ6duzI6aefzrXXXkuHDh1o27Ytbdq0ITU1lYSEhArPX1WFLTma2RsEpmpKNrPNBObRjAdwzo0hME3YlQRmMtlPYKYPEREJo6xVq+jiHPmxsXRxjqwyt6l2797NwoULmT9/PhkZGSxbtgy/3w9AmzZtOO+887jttts4++yzOfvss2nfvn21vK8ctuTonPv5MfY7YGi4fr6IiHyfI7Do5wXAAuBM58jMzGTOnDnMmTOHL774AucctWrV4vzzz+f3v/89vXr1okePHjRpcqxFZKqPSDariohIJeuQlsaHzjGtuJjFzrFg/Hie+/OfMTN69OjB448/Tr9+/ejVq1fUNomeCkqOIiI1xNdff824MWP4rriYpX4/cWac064d9w8bxhVXXFGjKsNjUXIUEanG9u/fz9tvv83YsWP5/PPPMaBpbCzXeDzsb9OGh8eMqbC3ak2l5CgiEoW8Xi8rMjM5Ozi8oqytW7fy0ksvMXbsWPbu3UtaWho/HTCAK9as4YzERKZs3Eib665TYixH9etiJCJSzXm9Xp4fPJjiZ5/l+cGD8Xq9h/f997//ZciQIbRr144///nP/PjHP2b+/PmsXLmS4b/5DQvj4th88CAHmzTh2gEDIvgtqjZVjiIiUWZWejo9N27k2vh4fLm5zEpPJyUlhccff5zXX3+dWrVqMWTIEIYPH0779u0Pf+5oA/rl6JQcRUSqqPKaTh3wH6Cdc/zbObI//pg/vfAC8fHx/Pa3v2X48OGkpBx9etFevXopKYZAzaoiIlXQoTlNG06YwNihQ49oOr12wAAKW7fmz2b86+BBvlq+nCFDhrBmzRqee+65chOjhE6Vo4hIFVC6SoTAHKhXAQOaN4dt246YA/XAgQN8W1DA+r176dOnD+PHj6dTp04RjL76UeUoIhJhZavE/Px8zu7WjTlA+rZtzAHO7taNnJwcbr31Vvr164eZ8f7775ORkaHEGAaqHEVEImxWejo9s7O5JDUVDh7kwP793+s8s2vXLq677jpyc3N59NFHeeSRR6hTp06kQ6+2lBxFRCJowoQJzJ86lUZ79uDLzeWL1FQGBZNer169OOuss3jwwQcZM2YMZ511Fh9++CFdunSJcNTVn5KjiEiEeL1epowYweB9+6gVF0dG3bp07N+funXrArBs2TJuuOEG1qxZw4MPPshTTz1Vo+c7rUy65ygiEiErMjPpa8Y8v5/9xcVkJyQcHpg/bdo0LrjgAvbv388nn3zCqFGjlBgrkSpHEZEIMY+HhTk5XOAcL/n9XH/LLXTt2pUbbriB9957j759+/LWW29pQvAIUHIUEYmQrFWr6JmYSNPERG71eKiVlES/fv347LPPePDBB3n22WeJi9Ov6UjQVRcRiQCv18vy2bOJ8/lI8flY0KwZK15+md27d/PEE0/whz/8IdIh1mhKjiIiETArPZ2r8vPp0LEjr27fzvxt20hMTGTevHns378/0uHVeEqOIiJhdLT5UQ9Vjd/m5pK0YwfvFRXRomVLMjIy6NChAxkZGZENWtRbVUQkXMqbH3VFZiaD69bl7JYtmVpYSLNmzcjMzKRDhw4RjlgOUXIUEQmTFZmZ9PH56FFQQB+fjxWZmUBgKrgXcnIY+d13tKhXj9ffeEOThVcxalYVEQkT83iYvGkTRcA/gFs8HgC++uor/rNhA2eecQajX3mFPn36RDRO+T4lRxGRMHGFhfROTia3qIje8fG4wkJeffVV7rvvPq688krS09M1sL+KUnIUEQmTdZs3s2D7dm4x48PYWNotXsykSZPo378/b7/9NrVq1Yp0iFIO3XMUEQkDr9fLwmnT6G3GVjOa1q3LlClTuPTSS5UYo4CSo4hIGByaN3WDGftLSng7J4e0tDTeeecdJcYooOQoIhIGh+ZNrVNSwl+KimjQsCHz58+nXr16kQ5NQqDkKCJyinm9Xj79+GO6NGjADOdIjIvjN7/6lYZrRBF1yBEROYUODfxvu28fz2/fTpFz/Pj00+n3ox9FOjQ5DqocRUROoVnp6fTYsYPlu3bh8/u56OKLeWTixMNTx0l0UOUoInIKeL1e3ktPxzt9Ott27GBlURFnp6Tw9KhRSoxRSMlRROQkHWpKTcjOptWOHXxSVETnOnW4atAgJcYopeQoInISvF4vY158kb4+H7WSk7ll0ybaejz0SkvjuoEDIx2enCAlRxGRE3SoYuycn8+EjRvJKikhNjaW/nfcwc2qGqOakqOIyAkoXTEO6tSJ9G3b2LJ3Ly+99BL3339/pMOTk6TkKCJynEpXjJM3beLTvXv5cu9ebr/9diXGakLJUUTkOM1KT6dndjY/S01lR0EBozdupFu3bowfPz7SockpouQoInIMXq+XFZmZnN2tGwDLZ8/m29xcdufkMKmkhHpJSbz//vvExsZGOFI5VZQcRUQqcKgJ9Spg7IQJNL/kEgbXrUurrl2586uvyDtwgLmzZ9O8efNIhyqnkJKjiEgFDjWhXpKaCgcP8iUwB0jNzmaFz8fAgQO57LLLIh2mnGJKjiIi5fB6vYebUH25uXyRmsrwAQPYd9ll/OyGG0hNTWXKlCmRDlPCQMlRRKQcKzIzGegcDZo1Y5bPR8f+/enVqxeDBg1ib14e//zwQ+rUqRPpMCUMNPG4iEg5zONh8qZN7Ny+nSyfjw7BxYqnTp3K7373O3r27BnpECVMlBxFRMrhCgvpnZxMbv369E5OJm/XLu6++266devGY489FunwJIzUrCoiUo51mzezYPt2bjHjw9hYCmbOZM+ePXzyySfEx8dHOjwJI1WOIiJH4fV6WThtGr3N2GpGy7p1+fzzz3nooYc4++yzIx2ehJmSo4jIUazIzKSvGRvMSPb7mbl7N61bt+bRRx+NdGhSCZQcRUSOwjweFubk0MY5RhQXs7+4mKlTp5KQkBDp0KQSKDmKiBxF1qpV9ExMpCg5md1+PxdfeCF9+/aNdFhSSdQhR0SkjEOD/2P37ePLggLiY2N57IknIh2WVCJVjiIiZcxKT+eq/HzaNWrEZr+fyy6/XFPE1TCqHEVESpkwYQLzp04lcfduPi0sJLl2bX6nTjg1jipHEZEgr9fLlBEjGLxvHyV+P/uBa264gfPPPz/SoUklU3IUEQk6NHzjneJi/l1cTGpiInfde2+kw5IIULOqiEiQeTws2LmTdUVFxADXDxlCr169Ih2WRIAqRxGRIFdYSOM6ddjoHD9u2JAzO3eOdEgSIUqOIiJBq9etY05uLq3NyD94EPN4Ih2SRIiaVUVECHTGeXvcOAqAfvHxtElJwRUWRjosiRBVjiIiwGcLFrAtL4+WMTFcaMZS4Oxu3SIdlkSIkqOICDBn7lwKnOPK2Fhe8vs555Zb1BmnBgtrcjSzy81slZmtMbOHj7K/vpnNNrOvzOwbM7s9nPGIiBzNmjVrmDdvHl0SEjizdWtu7dCBdq1aRTosiaCwJUcziwVGA1cAZwA/N7Mzyhw2FPjWOdcF6As8b2a6Ay4ileqee+7BnKO5GUXbt6tJVcJaOfYE1jjnspxzhcCbwLVljnFAPTMzIBHYBRSHMSYRkSN4vV7+/e9/c3WjRtzTvDkrExPp2L+/mlRruHD2Vm0JbCr1fjNQ9m/b34BZwFagHvAz55y/7InM7C7gLoCUlBQyMjLCEW+15/P5dO1OgK7biavq1845x7Bhw6hfvz79fv1rihIS6GVGo9atIxp3Vb9uNUE4k6MdZZsr8/7HwDLgEqAD8C8zW+icyzviQ86NA8YBpKWlOa2pdmIyMjK0Ht0J0HU7cVX92s2dO5fly5fToUEDdjz5JLvMcE2aEPPwwwwcODBicVX161YThLNZdTPQutT7VgQqxNJuB95xAWuAdcBpYYxJRAQIVI2PPfYYjRs14oZ69dhgRgvndL9RgPAmx8VAJzNrF+xkcyOBJtTSNgI/BDCzpkAakBXGmEREAJg9ezaLFy+m/zXX8HluLm2c0xAOOSxszarOuWIz+wUwF4gFJjrnvjGze4L7xwBPAZPNbAWBZtiHnHM54YpJRATA7/fz2GOP0bFjR5o1bkyTxESaJiZyq8dDioZwCGGePs459wHwQZltY0q93gpoeW0RqVTp6eksX76cP/zhDyx5803ifD5SfD6WpqYyXE2qguZWFZEapqSkhCeeeIIzzjiD5snJDHSOBs2aMcvn0xAOOUzJUURqlDfeeIOVK1cyY8YM9uzZw8QNG7ixpIQ1sbFcmJYW6fCkitDcqiJSYxQXFzNixAjOPfdcrr/+erJWraKLc+THxtLFObJWrYp0iFJFqHIUkRrjrbfeYs2aNcycOZOYmBgcsB64AFgAnBXR6KQqUeUoIjWC3+/n6aef5qyzzuKaa64BoENaGluAj0pK2BJ8LwKqHEWkhnjnnXdYuXIlb7zxBjExgboga9UqLq1fnxaJiXT2eLS4sRymylFEqj3nHCNHjqRz58789Kc/BQITji+fPZvVPp9W4pDvUeUoItXenDlz+Oqrr5g8eTKxsbEArMjM1DAOKZeSo4hUa845nnrqKdq2bctNN92E1+tlRWYm6zZv5rNNm7iZwJyVGsYhpSk5iki19vHHH/PFF18wduxYMjMzGTt0KFcBU7Oz6ZucDElJXFZcrPuNcgTdcxSRam3kyJG0bNmSQYMGsSIzkz4+Hz0KCuhrxlKgYa1afFe3ru43yhFUOYpItbVo0SIWLFjAX/7yF2rVqoV5PEzetIkiYCFw4bBh7G7Viru7ddP9RjmCkqOIVFujRo2iYcOG3HnnnXi9Xj79+GN6l2pKTWnViiH33hvpMKUKUrOqiFRLq1evZubMmdx333188803jB06lM7LlrEwJ4f4oiI1pUqFVDmKSLX0/PPP4/F4uPDCCxnz4ov09fkY1KkTABnnnss9w4apKVXKpeQoItXOjh07mDJlCueffz6v3nsvXYHJOYF11L9LTFRilGNSchSRaufll1+msLCQ/FWruCUvj/Pi4yE5mYyuXZUYJSS65ygi1YrP52P06NF07dKFKzwe5vn9eA8eZCkoMUrIlBxFpFp57bXX2LNnD5dfeSULc3Jo4xwv+f2cc8stSowSMjWriki14PV6WbZ4MX/84x/p3bs37Vq1IjE1lZS4OG4NDtsQCZWSo4hEPa/Xy/ODB1N/50625+by61//mrO7dWNs3bp0Bj6rVYtLNGxDjoOSo4hEvVnp6fTcuJG3CgpoZMauHTsCTaijR7MiM1Mz4MhxU3IUkajngPeKilhSUsKZsbFYcDHjXr16KSnKCVGHHBGJeh3S0vimpIR4IC42lg5afkpOkipHEYl6uTt2sMfv56K6denXuLGWn5KTpspRRKLevIULccBVJSUszMnBPJ5IhyRRTpWjiES1/Px8MjIy6FyrFrVbtOAyj0eVo5w0VY4iErW8Xi8DrruOgwcP0hQo2r6dpaDVNuSkKTmKSFTyer2Mue8+vBkZtIqJ4dHOnclq3JiO/furh6qcNCVHEYlKKzIzScnOZk9xMckxMWTt2cPBJk24dsCASIcm1YCSo4hEJfN4mLhtG/UAYmLw/uAH3D16tKpGOSXUIUdEolLmF1+wq6SESxs0oF+zZqRcdJESo5wyqhxFJOp4vV7mzpiBA3oGl6NSJxw5lVQ5ikjUSX/zTbbv2cPFSUnsql9fnXDklFPlKCJRxev18v60aeT7/Zxx4AA5deuqE46cckqOIhJV3psxg7xdu0g1Y7dz1OneXVWjnHJKjiISVbZs28YW57g0NpaDsbG0aNEi0iFJNaTkKCJRZdWaNcQABcAW0AocEhbqkCMiUWPnzp0sWbKEc2vXpkfz5pypeVQlTFQ5ikjUePLJJykpKaGBc5pHVcJKyVFEooLf7+ett94irXZt7m/RgpWJiRrCIWGj5CgiUeHjjz9m586dWHExWzdtYk1enu43StgoOYpIVHj11VepU6cOl5iRHxtLF+fIWrUq0mFJNaXkKCJV3ubNm5k9ezbt27Zlg3O0jo1lU2wsLtKBSbV1zORoZkvMbKiZNayMgEREynrttdfw+/2kHDzIvpgY3iksJC8lRTPjSNiEUjneCLQAFpvZm2b2YzOzMMclIgJAcXEx48ePp13bttxQVMSjaWk0btaMXjfcoM44EjbHTI7OuTXOud8DnYHXgYnARjMbYWaNwh2giNRss2fPZuvWrTQqLGRubi5frVmj+VQl7EK652hm5wDPA6OAdGAgkAd8Er7QRERg7Nix1EtMZLAZ93XsSFbjxhrCIWF3zBlyzOxLYA8wAXjYOVcQ3OU1s4vCGZyI1GwzZ85k7ty5tK1fn4937eLArl3kpKYySFWjhFko08f91DmXVXqDmbVzzq1zzv0kTHGJSA3n9XoZcd99ANxUUEDfzp2ZuXu3qkapFKE0q84IcZuIyCmzbPFiNuTk0K12bVY4x9rduznYpInuNUqlKLdyNLPTgDOB+mZWukJMAhLCHZiI1Gxfr1rFnuJizo+JYZEZ3h/8gHt++UtVjVIpKmpWTQOuBhoA/Utt3wfcGc6gREQ+/OAD6sbE0KFlS1rVqkXKRRcpMUqlKTc5OufeA94zswucc59XYkwiUsONGjWKNVlZdIyNxb9jB0tTUxmu1TekElXUrPpb59yfgJvM7Odl9zvn7g9rZCJSI3m9Xv42ciQA93g8fKvVNyQCKmpWXRl8XlIZgYiIAHy1ZAl78vNpFhNDot9PdkICd6kTjlSyippVZwefp1ReOCJS03373XfklZRwSVwcL/n9XH/LLaoapdJV1Kw6G8qf9N45d01YIhKRGm1BRgZ1YmI4LzmZsz0eUlu1inRIUgNV1Kz650qLQkQE2LFjB8u//prmMTE0zcvjH8AtHk+kw5IaqKJm1fmVGYiIyNSpUykpKeHHjRqR36ABl3k8uMLCSIclNVC5M+SY2fTg8wozW17qscLMlodycjO73MxWmdkaM3u4nGP6mtkyM/vGzJSQRWoo5xyjR4+mfq1abN+3j4ItW1gKnK0hHBIBFTWrPhB8vvpETmxmscBo4EfAZgLrQc5yzn1b6pgGwCvA5c65jWbW5ER+lohEv88++4wNGzbwU4+HQR4PUwsLqdO9uzrjSESUWzk657YFnzcABUAX4BygILjtWHoCa5xzWc65QuBN4Noyx9wEvOOc2xj8WdnH/xVEpDqYMGECcXFx5Pn9+IDC+HhatGgR6bCkhjLnyu2QGjjAbAjwOIG1Gw3oAzzpnJt4jM8NJFARDgm+vxXo5Zz7RaljXgTiCczhWg94yTk39Sjnugu4CyAlJaX79OnTQ/6C8j8+n4/ExMRIhxF1dN1OXKjXLj8/nwEDBtCjWzduGTiQeKA4Pp6W7dtTt27d8Adaxejv3Inr16/fl8658072PKEsWfUboKtzLhfAzBoDnwEVJkcCibSsspk4DugO/BCoDXxuZoucc6uP+JBz44BxAGlpaa5v374hhC1lZWRkoGt3/HTdTlyo1278+PEUFBTQd8sWmo0ezSyfj6aDBnHXL6YwNF4AACAASURBVH8Z/iCrIP2di7xQkuNmApONH7IP2BTi51qXet8K2HqUY3Kcc/lAvpktINB8uxoRqTEmTJhAixYtWLhzJ23MyAIuTEuLdFhSg1U0CcCvgy+3AF4ze49A5Xct8EUI514MdDKzdsFz3EjgHmNp7wF/M7M4wAP0Al44rm8gIlHtm2++wev18sO+fen+7bfkJyZqCIdEXEWVY73g89rg45D3Qjmxc67YzH4BzAVigYnOuW/M7J7g/jHOuZVm9iGwHPADrznnvj7eLyEi0etQR5yCNWv4ZvduGu/dy9J27bQKh0RURZMAjDjZkzvnPgA+KLNtTJn3o4BRJ/uzRCT6FBYWMm3aNDq0a8e1W7ZwuoZwSBVxzHuOZpYC/JZAj9KEQ9udc5eEMS4RqQHmzJlDTk4O/X7wAxZs3EhLMw7ExNBBQzgkwsod51jKP4D/Au2AEcB6AvcTRUROysSJE2nevDmX/vjHbAE+KilhC9BBnXEkwkLprdrYOTfBzB4Izrc6X9O8icjJ2rZtG//85z958MEH2bBmDZfWr0+LxEQ6qzOOVAGhVI5FwedtZnaVmXUlMCxDROSETZs2jZKSErp168by2bNZ7fNRtH275lOVKiGUynGkmdUHhgMvA0nAr8IalYhUa845Jk2axIUXXshXS5ZwVX4+HTp2ZObu3XTs31+dcSTijlk5Oufed87tdc597Zzr55zr7pybVRnBiUj1tGjRIv773//Sp08fFk2fzvvbt7N01Spy6tbl2gEDIh2eyLGTo5m1N7PZZpZjZtlm9p6Zta+M4ESkepo0aRJ16tShMD+fK3bu5F6Phy/9fg3hkCojlHuOrwPTgWZAC+Bt4I1wBiUi1Vd+fj7/+Mc/6HL22eTk5vIf4GBsrFbhkCollORozrlpzrni4OPvfH8CcRGRkIwaNYr9+/dz+e7d5H75JXkpKfy7fn2KU1PVpCpVRkVzqzYKvpxnZg8TWI/RAT8D5lRCbCJSDU2eNInk2FiGJieTUVDAl/37075dOwZ166YmVakyKuqt+iWBZHho6am7S+1zwFPhCkpEqqeZM2eyYeNGTouPZ9KyZXyRmsrwAQOUFKXKqWhu1XaVGYiIVG9er5cnH38cA5495xw+ys7WsA2pskKZWzUeuBf4QXBTBjDWOVdU7odEREqZMGECM0eMYP3WrdSPiWGvz8fBJk10j1GqrFA65LwKdAdeCT66B7eJiByT1+tlyogRnJ6dzZ6SEvo0bEhG167cPXq0qkapskKZIaeHc65LqfefmNlX4QpIRKqXFZmZ9DXj70VFJACFtWtzz7BhSoxSpYVSOZaYWYdDb4ITAJSELyQRqU7M42Hezp1s9vtJMKPrbbcpMUqVF0rl+BsCwzmyCPRcbQPcHtaoRKTacIWFNKpdm6IDB7iuSRPatdK6BVL1VZgczSwGOAB0AtIIJMf/OucKKiE2EakGzONh/u7dNAHW7tnDxR5PpEMSOaYKm1Wdc37geedcgXNuuXPuKyVGETkeXy5axF7n6B4Xx7lA1qpVkQ5J5JhCuef4kZkNMDM79qEiIkf6zOsF4Or4eDbFxmruSYkKodxz/DVQFyg2s4MEmladcy4prJGJSNTbt28fq1eupKEZ84uLyWvRQmMbJSocMzk65+pVRiAiUv0sWrSIA8XFDEtNZY9zdLjhBvVUlahQbrOqmTUxsxfN7H0ze8bMVCmKSMi8Xi/zPv6YOKDezp1ayFiiSkX3HKcC+cDLQD3gr5USkYhUC4s/+4zMZcu4pF491tSrp3lUJapU1KzazDn3++DruWaWWRkBiUj18OXy5RQUFnJOSQlfFBVxYVpapEMSCVlFydHMrCH/W7IqtvR759yucAcnItHr3x9/TLOUFJrk5dHFOQ3hkKhSUbNqfQJrOh56JAGZwddLwh+aiESrtWvXsmnzZi7o0YPUuDgN4ZCoU9F6jm0rMQ4RqUaeffZZAHp27847//qXhnBI1AllEgARkZD5/X7enTmTs+rU4fQmTUhs2JBeGsIhUUbJUUROqfnz55O7axclxcUUFRWR5fPRQZ1xJMqEMkOOiEjIJk+eTO2EBK5NTqY4Npbeycm4wsJIhyVyXEKqHM3sYjO7Pfg6xczahTcsEYlG+/btY8aMGaSddhqfZmcTV1zMwpwcTCtxSJQ5ZnI0syeAh4BHgpvigb+HMygRiU4zZsxg//79uK1b6W1GIahylKgUSuV4PXANgdlycM5tJTBjjojIESZPnkyjhg25zDk2x8XhMWMpcHa3bpEOTeS4hJIcC51zDgLDlMysbnhDEpFotHbtWhYsWECTuDi+y8+nQXExObGxXPf44+qpKlEnlA45081sLNDAzO4E/g8YH96wRCTaTJ06FTPjrnr1aF+vHrN8Puo1bszNN98c6dBEjlsoS1b92cx+BOQBacDjzrl/hT0yEYkafr+fqVOncsYZZzBr3TpuBrKAK2rVinRoIifkmMnRzH4FvK2EKCLlWbBgAevXr+faq68mbedO8hMTuczjAadJ4yQ6hXLPMYnAqhwLzWyomTUNd1AiEl0mT55M7dq12fHFF6zcu5ei7dtZCtSuUyfSoYmckGMmR+fcCOfcmcBQoAUw38w+DntkIhIVfD4f06dPp35MDHf4fPw0Lo6ViYl07N+funXVf0+i0/FMH5cNbAdygSbhCUdEos3bb7/NgQMHuCIxkXl+P/uLi8lOSNBE4xLVQpkE4F4zywD+DSQDdzrnzgl3YCISHSZPnkzTpk3JysujjXO85Pdzzi23aPiGRLVQhnK0AYY555aFOxgRiS6HxjZef801nLVsGXWLihgYH09qq1aRDk3kpJRbOZpZUvDln4CNZtao9KNywhORqmzy5MnExMTQqk0b5m3fToNduzSXqlQLFVWOrwNXA18SmB3HSu1zQPswxiUiVZzf72fKlCn07NmTZe++S28ztpppLlWpFsqtHJ1zVwef2znn2gefDz2UGEVquE8++YRNmzZxzhln0NeMDWa0cE5zqUq1EEqHnH+Hsk1EapZJkybRoEEDuvbowcKcHHXGkWql3GZVM0sA6gDJZtaQ/zWrJhEY7ygiNdSePXt45513uPLKK/HOn0/v5GRSk5K4tbiYFHXGkWqgonuOdwPDCCTCL/lfcswDRoc5LhGpwqZPn87Bgwcp/PZbOgMf5eQwuHZtvktM5BI1qUo1UG5ydM69BLxkZr90zr1ciTGJSBU3adIkkpOTudLn48Y2bQDIOPdc7hk2TE2qUi2EsirHy2Z2FnAGkFBq+9RwBiYiVdObb77JokWL6NigAR/t2kX+rl0sTU1luBKjVCOhrMrxBNCXQHL8ALgC+BRQchSpYbxeL08PG4YBPy8ooHenTszcvZuO/fsrMUq1EsrcqgOBHwLbnXO3A10ALdImUgMtW7yYTTk59KhTh2V+P2t37+ZgkyaaR1WqnVCS4wHnnB8oDs6ak40mABCpkb5ZvZq9JSX0KC5mC+D9wQ+4e/RoVY1S7YQyt+oSM2sAjCfQa9UHfBHWqESkSvp0wQLqxsTQolEjrvJ4SL3oIiVGqZZC6ZBzX/DlGDP7EEhyzi0Pb1giUtVkZ2ez/OuvaRETQ5O8PP4B3KI5VKWaqmgSgHIHK5lZN+dcZnhCEpGq6O9//zslJSX8qFEj8hs04DKPR3OoSrVVUeX4fAX7HHDJKY5FRKoo5xyjR48myeNhx759FOTns7RdO4ZrwL9UUxVNAtCvMgMRkapryZIlZGVlMcDj4XaPh6mFhdTp3l33G6XaCmWc421H265JAERqjokTJxIXF0c+kFdSwoGYGDq00BTLUn2F0lu1R6nXCQTGPGaiSQBEaoT9+/fz+uuv06NHD7IzM/mopIQtMTFcm5YW6dBEwiaU3qq/LP3ezOoD00I5uZldDrwExAKvOeeeK+e4HsAi4GfOuRmhnFtEKsfMmTPJy8sjKSGBnk2bkpqUROfiYnXGkWotlEkAytoPdDrWQWYWS2D1jisITD33czM7o5zj/gjMPYFYRCTM/vKXv1DP46H31q0szMkhvqiI7+rW1YLGUq2Fcs9xNoHeqRBIpmcA00M4d09gjXMuK3ieN4FrgW/LHPdLIJ0jm29FpApYs2YNmZmZDEhO5tb4eAqSk8no2lWrb0i1F8o9xz+Xel0MbHDObQ7hcy2BTaXebwaO+NdkZi2B6wkMCyk3OZrZXcBdACkpKWRkZITw46Usn8+na3cCavJ1Gz9+PGZGj/vuw1u/Pi2Bc1JTOXDgQEjXpCZfu5Oh6xZ5odxznA8QnFc1Lvi6kXNu1zE+akfZ5sq8fxF4yDlXYna0ww/HMA4YB5CWlub69u17rLDlKDIyMtC1O3419boVFxdz00030aF9e7b88Y/EmrHROerdfz9P/+lPIZ2jpl67k6XrFnnHvOdoZneZ2Q5gObCEwPyqS0I492agdan3rYCtZY45D3jTzNYTWP3jFTO7LoRzi0iY/fOf/2Tbtm20btmSLL+f1rGxbIqN/d7/cEWqo1CaVX8DnOmcyznOcy8GOplZO2ALcCNwU+kDnHPtDr02s8nA+865d4/z54hIGLz22ms0btyY2jt2sC8mhncKC8lr0ULLU0mNEEpv1bUEeqgeF+dcMfALAr1QVwLTnXPfmNk9ZnbP8Z5PRCrPtm3bmDNnDj26d+cG4NctW5LYsCG9brhBHXGkRgilcnwE+MzMvEDBoY3OufuP9UHn3AfAB2W2jSnn2MEhxCIilWDKlCmUlJTQu08fJj/9NDcDWcCFGvgvNUQoyXEs8AmwAvCHNxwRiTTnHK+99hp9+vShScOG9E5OJreoiN7x8Rr4LzVGKMmx2Dn367BHIiJVwvz581m7di19L76YdZs381lODjcDH6L1G6XmCCU5zguOM5zNkc2qxxrKISJR6Nlnn8UTG8sly5cz5pNP6JucDElJXKYp46QGCSU5Huph+kipbQ5of+rDEZFIys3N5ZNPPqFvvXpcXFTEf81YCnSpVYvPatXiEk0ZJzVEKJMAtDvWMSJSPUyZMoXi4mL2HDzIR1lZLAQuHDaM3a1acXe3buqpKjWG1nMUESDQEWfcuHF0aNeOK0pKDnfCSW3ViiH33hvp8EQqVSjjHHuUevQG/gBcE8aYRCQCFi5cyKpVqzj9zDOZt307DXbtYmFODqZOOFIDhXU9RxGJHuPGjSMxMZFdS5fSx4ytZvROTlYnHKmRwraeo4hEj9zcXGbMmEGP7t35YWwsG8xo4RxLQes2So0UzvUcRSRKTJs2jYKCAvr060fG4sVc4Bwv+f1cf8st6oQjNVI413MUkShwqCPO+eefT2F+Pj0TE2mamMitHg8prVpFOjyRiCi3WdXMOprZRc65+aUe/wHamVmHSoxRRMLo008/ZeXKlVx66aUsnz2b1T4fRdu3q0lVarSK7jm+COw7yvYDwX0iUg28+uqrJCUl0axxYwY6x+BmzViZmEjH/v3VpCo1VkXJsa1zbnnZjc65JUDbsEUkIpUmOzubGTNmMHjwYLbu3Mlr69axfcsWsnw+OmgFDqnBKkqOCRXsq32qAxGRyjdx4kSKioq44IILWDhtGr01hEMEqDg5LjazO8tuNLM7gC/DF5KIVIaSkhLGjh1L37598e3eTV8zDeEQCaqot+owYKaZ3cz/kuF5gAe4PtyBiUh4zZ07l/Xr1/PHP/6Rffv2sTAnR0M4RILKrRydczuccxcCI4D1wccI59wFzrntlROeiITLK6+8QtOmTbnuuuvIWrUqMISjZUtu7dCBdhrCITVcKNPHzQPmVUIsIlJJ1q9fzwcffMDvf/97li5dyvLZs4nz+Ujx+ViamspwNalKDXci08eJSJQbN24cZsZdd93FisxMDeEQKUPJUaSGKSwsZMKECVx99dW0bt0a83iYuGEDWzdtYk1enoZwiKDkKFLjpKenk52dzb3BNRqzVq2ii3Pkx8bSxTmyVq2KcIQikafkKFLD/PWvf6VTp05cdtllQGBVgfVAy+CzK/eTIjWHkqNIDbJkyRIWLVrE0KFDiYkJ/PPvkJbGFuCjkhK2BN+L1HShrMohItXEyy+/TGJiIoMHDz68zRUWclXTptQtKqJ1fLxmxhFBlaNIjZGdnc2bb77JoEGDqF+//uHt5vGwMCeHxnv3sjAnB/N4IhilSNWgylGkhhg/fjyFhYX84he/OGK7Kyykd3IyuUVF9FblKAKochSpEYqKinj11Vf50Y9+xGmnnXbEvnWbNzNv+3Ya7NqlylEkSJWjSA0wc+ZMtmzZwquvvnrEdq/Xq9U4RI5ClaNIDfDyyy/Tvn17rrzyyiO2z0pP54KCAjbHxdE6JkarcYgEKTmKVHPLli3j008/ZejQocTGxh7e7vV6WT57Nqt9PhoUFzO+Xj2ue/xxTR0ngppVRaq9F154gbp16/J///d/R2w/NKdqg2bNmOXz8cPbbuOOO+6IUJQiVYsqR5FqbOvWrbzxxhvccccdNGjQ4Ih95vEwedMmdm7fTpbPp8H/IqWochSpxkaPHk1xcTEPPPDA9/YdWsMxPzGRyzwedcQRKUWVo0g1tX//fsaMGcP1119P+/btj9jn9XpZNH063+zeTcGWLeqII1KGKkeRamrq1Kns2rWLX/3qV9/bNys9nSt27uR0j4ephYXU6d5dHXFESlFyFKmG/H4/L7zwAj169OCiiy46Yp/X6+WrZcuIKSmhXa1aFMbH06FFiwhFKlI1KTmKVEMffPABq1ev5o033sDMDm/3er2MHTqUi/LzmQF8kJBAcaNGXDtgQOSCFamClBxFqqG//OUvtG7dmgFlkt6KzEz6+Hz0i4/nQNOmbLrwQh4dNkxNqiJlKDmKVDNLly5l3rx5jBo1ivj4+CP2rdu8mQXr1nHQjIWxsdxy6aVKjCJHod6qItXMn/70J5KSkhgyZMgR2zWPqkjolBxFqpG1a9cyffp07rnnnu8N+tc8qiKhU3IUqUaef/554uLiGDZs2BHbNY+qyPHRPUeRamLHjh1MnDiRQYMG0bx58yP2aR5VkeOjylGkmvjrX/9KYWEhDz744Pf2aR5VkeOjylGkGsjLy+OVV17hJz/5CZ07d/7efldYSO/kZHKLiugdH6+OOCLHoMpRpBoYN24ce/bs4aGHHjrqfvN4mJ+dTVJODvOzszGPp5IjFIkuSo4iUa6goIAXXniBSy65hB49ehz1mKxVq+jiHPmxsXRxjqxVqyo5SpHoouQoEuWmTJnC1q1befjhh8s9xgHrgZbBZ1cpkYlELyVHkShWVFTEM888Q69evbj00kvLPa5DWhpbgI9KStgSfC8i5VOHHJEoNm3aNDZs2MArr7xyxATjZWWtWsWl9evTIjGRzlrYWOSYVDmKRKni4mKefvppunfvzhVXXFHucaUnACjavl0z44iEQJWjSJR6/fXXycrK4r333iu3avR6vYx58UUGOsdpXbsyZeNGOvbvr5lxRI5ByVEkCpWUlDBy5Ei6dOlC//79j3rMhAkTePfJJ+kKTM7JYTBwsEkTrd0oEgIlR5Eo9NZbb/Hdd9+Rnp5+1KrR6/UyZcQIbsvN5bz4eEhOJqNrV+7R2o0iIVFyFIkyh6rGs846i+uuu+6oxxxagWOBGRQVsRS0qLHIcVByFIkyb731FitXruStt94iJub7feoOdcCJ8/loXVLC+AYNuEsrcIgcF/VWFYkiRUVFPP7445xzzjkMHDjwe/tLd8B5uGtXips00QocIidAlaNIFJk0aRJr165l9uzZ36savV4vY4cOpXN+PpM3bVIHHJGToOQoEiUOHjzIk08+yQUXXMBVV111eLvX62VFZibr1q3jKmBAx44AZJx7rjrgiJygsCZHM7sceAmIBV5zzj1XZv/NwKFlBHzAvc65r8IZk0i0evXVV9myZQvTpk073EP10HCNgbVrs9yM7wC2beO7unWVGEVOQtiSo5nFAqOBHwGbgcVmNss5922pw9YBfZxzu83sCmAcoH/NImXs27ePZ555hksvvZR+/foBRw7X6BIfD82aseqaa9jdrh13d+umxChyEsJZOfYE1jjnsgDM7E3gWuBwcnTOfVbq+EVAqzDGIxK1XnzxRXJycnj66aeB/3W86WvGZ/HxUFTEewcO8OiAAUqKIqeAOReexWvMbCBwuXNuSPD9rUAv59wvyjn+QeC0Q8eX2XcXcBdASkpK9+nTp4cl5urO5/ORmJgY6TCiTqSv2969e7n55pvp2rUrTz31FPn5+eRs3Egtv5+8wkIS4+LYBzRu3pzk5OSIxXk0kb520UrX7cT169fvS+fceSd7nnBWjkeb7PGomdjM+gF3ABcfbb9zbhyBJlfS0tJc3759T1GINUtGRga6dscv0tftgQce4MCBA4wZM4YzzzyT3//2t7R+802uSU1lTG4un1fhjjeRvnbRStct8sKZHDcDrUu9bwVsLXuQmZ0DvAZc4ZzLDWM8IlFn9erVvPLKKwwZMoQzzzzz8AD/b3Nz8eXmsjQ1leFVNDGKRLNwTgKwGOhkZu3MzAPcCMwqfYCZpQLvALc651aHMRaRqPTQQw+RkJDAk08+CQSmhbsqP5/7OnYkq3FjrbAhEiZhqxydc8Vm9gtgLoGhHBOdc9+Y2T3B/WOAx4HGwCvBrunFp6KtWKQ6mD9/Pu+++y4jR46kadOmTJgwgflTp9Jozx5+kJtLTmoqgzTAXyQswjrO0Tn3AfBBmW1jSr0eAnyvA45ITef3+3nwwQdp2bIlv/rVrw4P2xi8bx+14uLIqFtXVaNIGGluVZEq6I033mDJkiU888wzrFix4ohhGwXOkZ2QoGnhRMJI08eJVDEHDhzgd7/7Hd26daNTp06H50v9KCeH3snJTHWOQVplQySslBxFqpjnnnuOjRs3csNPfsL7M2ceMV/q6nPPZZR6p4qEnZKjSBWydu1annvuOVonJdEhPZ23Y2L4rnZtzZcqUsmUHEWqCOcc999/PzjHoMJCrs3Lw1dUxNfXX8/uiy7SfKkilUjJUaSKmD17Nh988AHndevG0hUrWFBczH+A01u0YMi990Y6PJEaRb1VRaqAAwcO8MADD9CuXTua5OezLyaGdwoLyUtJUa9UkQhQchSpAn75y1+yfv16LuzZk/9LTORP551H42bN6HXDDWpKFYkAJUeRCHv77beZPHEiFzdsyN6vvuKt/Hw2HzzIwSZNVDWKRIjuOYpEiNfr5aslSxj53HMkxMQwvXt3Ptu9my8vuUQLFotEmJKjSAR4vV6eHzyYuB072LR7N+c1a8Znu3czB7hbCxaLRJySo0gEzEpPp/P69bx08CAdY2I445JL2H3xxaoWRaoIJUeRCNi8ZQsfFhRQBLTzeGjZsqWGa4hUIeqQI1LJvF4vK+bNI9s5TouJoaRpU3W8EaliVDmKhJnX62VFZiZnB5tMP503jzXZ2ZyRmMj5DRrQSMM1RKocJUeRMPJ6vYwdOpSrgLETJuD+9jfSZ81iv9/PnampLK9dW1WjSBWk5CgSRisyM+nj89EjPh5fUREvv/gin3/+OcOHDyexQwd1wBGpopQcRcJo3ebNLFi3joNmTDBjWVYWl19+OaNGjcLMIh2eiJRDHXJEwsTr9bJw2jR6m7EJ2FRSQq1atZg4caISo0gVp8pRJExmpadzQUEBm+PiyCkoYHtxMX8cOZLmzZtHOjQROQYlR5Ew8Hq9LJ89mzifD39REf8sLqZ379789re/jXRoIhICNauKnGJer5cxL77IQOcYfMYZfFRSQkpyMh9++GGkQxOREKlyFDmFDg3d6Jyfz8SNG9nsHCXA6FdeoU6dOpEOT0RCpMpR5BQ5VDH28fl4pGNHCmvVIuvgQUaMHMlPf/rTSIcnIsdBlaPIKTBhwgTeffJJugKTc3JYlJfHoj17uPHGG3nkkUciHZ6IHCclR5GT5PV6mTJiBLfl5nJefDzbk5J4bedOzjnnHKZOnRrp8ETkBCg5ipykFZmZ9DVjnt9P7oEDTPP5aBzsgBMfHx/p8ETkBOieo8hJMo+HhTk5NPf7+UNREf7YWObOnavxjCJRTMlR5CS5wkIubtyYD2NiKHCOu+64gy5dukQ6LBE5CWpWFTlJ5vEwZds2NhUX0yE+nnO7d490SCJykpQcRY5T2fUZP/rnP9lUXEyfRo24LCUFV1gY6RBF5CQpOYoch0NDNgbWrs3YxERmX34509PTadugAfe1bcuHZvywW7dIhykiJ0nJUSREpYds/H97dx8cVZXmcfz75IWBtHESNkpwgBoUR0oKGaMI4vqyAgFHkWxp4boju4MvoMXU1OiwO4usSm051qhV67CILxlAyjUjhZsF4ssos6hIrUsGiIMwI0IGHJGEFwFBmoROJ2f/uDexaTqh053Qnc7vU3Wr+3afc/vcw0k/nNvnnDsqN5eNgQBPPPEEkyZNYt68eXy6bZvuzyiSIRQcReIUOWVj64kTPHf0KCMvu4zKykoCgQDXXnttqosoIl1Eo1VF4tQ6ZSPc3MzCpibOHzCAdevWEQgEUl00Eeli6jmKtKN14M3AIUMAb8rGoMJCXq2ro/hb3+KfHnqIgoKCFJdSRLqDgqOIL3IUKsCLs2dzM/DlPfdQXV3NvkOHqNi7l2GBAKOHDWPc9dentsAi0m0UHEX45lZTNwMvLlnCwBtv5GbgtoEDWQU8/+yzvFxRwfcvv5y7p09n9LhxGngjksEUHEXwBttcf/w4o3NzOd7UxKfAm4Crq+OVN95g9VtvUVpaysqVK3VfRpFeQANyRPAG2yzbs4c1u3axbM8eLrrkEmY88wxPZ2ez+q23mDFjBq+//roCo0gvoZ6j9FqRvzG6UIjSIUMgJ4fScJhjhw/z8MMP8/uNG7n33nspLy/HzFJdZBE5SxQcpVeK/o3xmgceYGcgwPeAquZmNv/qVxw6dIjly5czYMAABUaRXkaXVaXXqK6u0roUOgAADOFJREFUZvHzz7f1GFsH3NyMN01j5rPPsvLii3m7tpasrCzWrVvHHXfckepii0gKKDhKRokMgNGvvzh7NoVLlvDi7NlYnz68CVTW1/MmcNHw4SxcuJCK5csZP348H330kUajivRiuqwqGSP6UimLFrUFuMieIvX1HAmFmLVoEVtrahiXm8v9999PbW0tjz/+OHPnziUrS/9vFOnN9A0gPV5rb7GqsvKUS6Vba2ra0owsKTmlpziypIRRo0axY/duZs2aRWNjI2vXrmXevHkKjCKinqP0bJG9xdXBIDsB/AA4K+LWUWPGjAG/p9j6eklJCZ988gkzZ87k6aef5txzz03FKYhIGlJwlB4ncgpG9OXSzTfeyJGhQ2PeOmrMmDEMHTqUuXPnsnTpUgYNGsQ777xDaWlpak5ERNKWgqP0KLGmYLwJ3/QWb7st5kCacDjMokWLeOyxxwgGg8yZM4dHHnlEvUURiUnBUXqUqspKrjpwgBuHDIHGxlMG1sTqLba0tPDaa68xf/58tm/fTmlpKQsWLGD48OEpOgMR6QkUHKXHqK6uZsOKFWzZt49jBw6waehQfuYHxOig6JyjqqqKRx99lI8//pgRI0awatUqbr31Vk3oF5Ez0rA86TGqKiu56eBBHujTh80tLeRdccVpQTEUCvHKK69QUlJCWVkZDQ0NVFRUsGXLFqZOnarAKCJxUc9R0l7rAJy9dXX8CRianU0oN5eLLrigLc3hw4cpLy9n4cKF1NXVcemll/LSSy9x1113kZOjZi4inaNvDUlrkQNwDgWDnDjvPNa2tBAOBLilrIw1a9awdOlSVq5cSSgUYuLEiSxZsoRJkyaplygiCVNwlLTU2lvcvXv3KVM1Nt1yC7m5uZyoq+POO+/k888/p7CwkFmzZnHfffcxcuTIVBddRDKAgqOkneiJ/Tuc49OdO3n1q6848uqr7N27l+zsbMaPH89TTz3F1KlT6du3b6qLLSIZRMFRUiJyIn/0oJqtNTVcFQpxLCuLw/v3U9PQwH81NpKdnc2ECROYP38+ZWVlFBUVpaj0IpLpFByly3UU+Frfj5zIH16wgEAgwIYNG1i/fj1r165l//79APTLzWXCpElMnz6dCRMm0L9//7N8NiLSGyk4Spfq6M4Y4I0qrVyxgpyDB1lrxvqDB/nPG24gHA4DUFxczHXXXcfgwYPJ79uXyVOmMHbs2BSdjYj0VgqOErcz9QgBtmzaxLimJgbn51NYX8+TTzxBYVERtbW17Nixg3379rWl7ZeVxbfz8pg2bRpTpkxh9OjRXHjhhRplKiIp163B0cwmAwuAbGCxc+6XUe+b//4PgBPAj5xzNacdSFLGOUdDQwPvvvsuS37+c0rCYX4TDjNm2jTy8/Opr6+nvr6effv2UVdXx549e9p6gQDs2kVxcTHDhg1j8uTJjBgxghEjRnDy5EkO1tVxWYyJ/CIiqdZtwdHMsoFFwETgC2CjmVU55/4Ukewm4GJ/GwM87z9mJOfcKVtLS0vM/ZaWlratubm57THWFg6H27ampqZTtlAoRCgU4uTJk2zfvp39+/ezatUqBgwYQENDA42NjQSDQU6cOEEwGCQYDPL111+3bceOHePo0aOEQqG2c1jlP773S+//OQUFBQwcOJDi4mKuvvpqpk2bRnNzM43HjzPmmmsoKyvjnHPOSUFti4gkrjt7jlcBtc65XQBmthyYCkQGx6nAy845B2wwswIzG+icq2/voDt37kzqy9b7qPjej07b3n6sx9Ytns9MlX79+pGXl0cgECAvL4+8vDzy8/MZNGgQ+fn55OfnU1hYSEFBAV999RUbVqzg+pwctuTkcM+TTzJx4kRNoRCRjNSdwfE7wJ6I/S84vVcYK813gFOCo5nNBGb6uyeDweC2ri1qZjM4rxiKmsGywe2DLx0cbGhooKGhgUOHDsV7qMAHkOfgxOpbbw12Z5nTTBHwZaoL0UOp7hKjekvcJV1xkO4MjrFGVUR3oeJJg3OuHCgHMLNNzrkrky9e76O6S4zqLXGqu8So3hJnZpu64jjdeVeOL4DBEfuDgLoE0oiIiJxV3RkcNwIXm9lQM+sD/B1QFZWmCvgH84wFjnb0e6OIiMjZ0G2XVZ1zYTP7MfAO3lSOpc65P5rZ/f77LwBv4U3jqMWbyjEjjkOXd1ORewPVXWJUb4lT3SVG9Za4Lqk7S9eRlCIiIqnSnZdVRUREeiQFRxERkShpFRzNbLKZfWpmtWb2LzHeNzP7D//9j82sJN68mSzJevvMzLaa2R+6agh0TxJH3Q03s/8zs5NmNqczeTNZkvWmNtdx3f3Q/zv92Mw+NLNR8ebNZEnWW+fbXPSSZqna8Abt/Bm4EOgDbAEujUrzA+C3ePMjxwLV8ebN1C2ZevPf+wwoSvV5pHHdnQ+MBn4BzOlM3kzdkqk3/z21uY7rbhxQ6D+/Sd9zydWbv9/pNpdOPce25eaccyGgdbm5SG3LzTnnNgAFZjYwzryZKpl66+3OWHfOuQPOuY1AU2fzZrBk6q23i6fuPnTOHfF3N+DN/44rbwZLpt4Skk7Bsb2l5OJJE0/eTJVMvYG3ItEaM9vsL9PXmyTTbtTmvtHZc1eb+8aZ6u4evKs+ieTNJMnUGyTQ5tLpfo7JLDcX1zJ0GSrZZfqucc7Vmdn5wO/MbLtz7oMuLWH6SqbdqM2dqjPnrjZ3qph1Z2Z/g/cl/9edzZuBkqk3SKDNpVPPMZnl5nrzMnRJLdPnnGt9PACsxLt80Vsk027U5r7RqXNXmztz3ZnZZcBiYKpz7lBn8maoZOotoTaXTsExmeXm4smbqRKuNzMLmFk+gJkFgFKgN93xJJl2ozaXwLmrzZ257sxsCPDfwHTn3I7O5M1gCddbom0ubS6ruiSWm2svbwpO46xLpt6AAcBKMwOvLfzGOff2WT6FlImn7sysGNgEnAu0mNlP8UbJHVOb63y94d2KSW2u47/XR4G/Ap7z6ynsnLtS33OJ1RsJfs9p+TgREZEo6XRZVUREJC0oOIqIiERRcBQREYmi4CgiIhJFwVFERCSKgqNIFDNr9lfv32Zmr5lZXjd+1r+Z2QT/+U87+1n+3NV3zexcf/8nZvaJmVV0Qdl+ZGYXROwvNrNLEzzWj81sxplTiqQHTeUQiWJmx51z5/jPK4DNzrl/jyOf4f1NtST4uZ8BVzrnvuxEnpuBCc65B/397cBNzrndUelynHPhTpbnfbw7aiR9Wyk/6P+vc+7yZI8lcjao5yjSsfXAMAAze8jvTW7zJ7VjZt/1e2rPATXAYDO707x7x20zsyf9dNlmtsx/bauZtQazZWZ2u5n9BLgAeM/M3jOze8zsmdZCmNl9ZhYrQP8QWO2neQHvlj5VZvagmc03s3IzWwO87Jd1vZnV+Nu4iOP/s1+uLWb2SzO7HbgSqPB70f3M7H0zu9JPf9o5+q8fN7Nf+MfZYGYDAJxzJ4DPzKw3LRUnPVmy99nSpi3TNuC4/5iDF3geAK4AtgIB4Bzgj8DlwHeBFmCsn+cC4HPgPD//u0CZn/93EZ9R4D8uA273n3+Gf885/3P+DOT6+x8CI2OU9S9AfsR+5DHmA5uBfv5+HtDXf34xsMl/fpN//Dx/v7//+D5eT5bI/fbO0U/jgCn+86eAf43IPw/4War/fbVpi2dTz1HkdP3M7A94y599DizBW+F/pXMu6Jw7jreG47V++r847z6Z4N3g933n3EHnXcasAK4DdgEXmtlCM5sMHOuoAM65IF7QucXMhuMFya0xkvZ3zn3dwaGqnHMN/vNc4NdmthV4DW85N4AJwEvO693hnDvcUdk6OEeAEPCG/3wz3n8eWh3AC6wiaS9t1lYVSSMNzrnvR77g/57YnmBk0lgJnHNHzGwUMAmYDUwD7j5DORYDDwPbgZfaSRM2syzX/u+ckWV7ENgPjML7SaUxosydGXzQUV00Oedaj9XMqd8xfYGG07OIpB/1HEXi8wFQZmZ5/sr+f4v3e2S0auB6Mysys2zgTmCdmRUBWc65SuARoCRG3q+B/NYd51w13m16/h54tZ1yfYr3O2M8vg3U+4F0Ot4CzgBrgLtbR8qaWf9Y5TnTOcbx+d+jd92BQ3owBUeRODjnavB+H/w9XnBY7Jz7KEa6emAu8B6wBahxzq3Gu2v5+/7l2mV+mmjlwG/N7L2I11bgjfI80k7R3gRuiPM0ngP+0cw24AWqoF/mt/Fu/7PJL98cP/0y4IXWATlxnOOZXAP8T5xlFUkpTeUQSWNm9gbwjHNubTvvDwReds5NPLsl6xwzuxx4yDk3PdVlEYmHeo4iacjMCsxsB97vnzEDI7T14n7dughAGivCu5ws0iOo5ygiIhJFPUcREZEoCo4iIiJRFBxFRESiKDiKiIhEUXAUERGJ8v/stfj1n/oTTQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "por_values = np.linspace(0.0,0.25,100) \n",
    "fit_mean, fit_stdev = norm.fit(por,loc = por_average, scale = por_std) # fit MLE of the distribution \n",
    "cumul_p = norm.cdf(por_values, loc = fit_mean, scale = fit_stdev)\n",
    "\n",
    "# plot the cumulative probabilities vs. the sorted porosity values\n",
    "plt.subplot(122)\n",
    "plt.scatter(por_sort, p, c = 'red', edgecolors = 'black', s = 10, alpha = 0.7)\n",
    "plt.plot(por_values,cumul_p, c = 'black')\n",
    "plt.xlabel('Porosity (fraction)'); plt.ylabel('Cumulative Probability'); plt.grid(); \n",
    "plt.title('Nonparametric Porosity CDF')\n",
    "plt.ylim([0,1]); plt.xlim([0,0.25])\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.2, hspace=0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comments\n",
    "\n",
    "This was a basic demonstration of univariate statistics in Python.\n",
    "\n",
    "I have other demonstrations on the basics of working with DataFrames, ndarrays, univariate statistics, plotting data, declustering, data transformations, trend modeling and many other workflows available at [Python Demos](https://github.com/GeostatsGuy/PythonNumericalDemos) and a Python package for data analytics and geostatistics at [GeostatsPy](https://github.com/GeostatsGuy/GeostatsPy). \n",
    "  \n",
    "I hope this was helpful,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "#### The Author:\n",
    "\n",
    "### Michael Pyrcz, Associate Professor, University of Texas at Austin \n",
    "*Novel Data Analytics, Geostatistics and Machine Learning Subsurface Solutions*\n",
    "\n",
    "With over 17 years of experience in subsurface consulting, research and development, Michael has returned to academia driven by his passion for teaching and enthusiasm for enhancing engineers' and geoscientists' impact in subsurface resource development. \n",
    "\n",
    "For more about Michael check out these links:\n",
    "\n",
    "#### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)\n",
    "\n",
    "#### Want to Work Together?\n",
    "\n",
    "I hope this content is helpful to those that want to learn more about subsurface modeling, data analytics and machine learning. Students and working professionals are welcome to participate.\n",
    "\n",
    "* Want to invite me to visit your company for training, mentoring, project review, workflow design and / or consulting? I'd be happy to drop by and work with you! \n",
    "\n",
    "* Interested in partnering, supporting my graduate student research or my Subsurface Data Analytics and Machine Learning consortium (co-PIs including Profs. Foster, Torres-Verdin and van Oort)? My research combines data analytics, stochastic modeling and machine learning theory with practice to develop novel methods and workflows to add value. We are solving challenging subsurface problems!\n",
    "\n",
    "* I can be reached at mpyrcz@austin.utexas.edu.\n",
    "\n",
    "I'm always happy to discuss,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "Michael Pyrcz, Ph.D., P.Eng. Associate Professor The Hildebrand Department of Petroleum and Geosystems Engineering, Bureau of Economic Geology, The Jackson School of Geosciences, The University of Texas at Austin\n",
    "\n",
    "#### More Resources Available at: [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}