PythonDataBasics_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"
   ]
  },
  {
   "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": 15,
   "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",
    "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": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(\"c:/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": 3,
   "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": 3,
     "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": 4,
   "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": 5,
   "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": 6,
   "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": 7,
   "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": 8,
   "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": 9,
   "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": 10,
   "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": 11,
   "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": 14,
   "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": 15,
   "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": 16,
   "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": 17,
   "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": 18,
   "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": 19,
   "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": 20,
   "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": 21,
   "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": 22,
   "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": 23,
   "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": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAccAAAGWCAYAAAANLkjgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dfXycVZ3//9enaaaFhFIgaUtvAqWFrAqWFmgXaqHcLMJWrPzaRdQi1SKIRbeK6xdXlhV01/26y3dh3a4ULRbQla1GJQgriBKoIpfQKRJBQ0sKvUljmlBLJ7RN0nx+f8w1ZTrkZpJ0Mnfv5+ORR2auuebKZw6hn5xznfM55u6IiIjIW0ZkOwAREZFco+QoIiKSQslRREQkhZKjiIhICiVHERGRFEqOIiIiKZQcRfKImd1lZv+Q7TgGw8xeNLP52Y5DJB1KjjKszOxVM/uTmZUlHbvGzOqyGFbWmdlSM/tVf+e5+yfd/SuDuP6rZrbXzGJh+3/HzMoHF+3guPu73L0ujOfLZvbdoVzPzN5rZk+Z2R4z22lmT5rZ+8PXlprZgfDzxsxsc/iZT0l6/4lm5knnxMzsd0P6kFIwlBwlG0YCf5vtIPpiZiOzHUMqMysZ4iUuc/dyYBZwFnDzIGLIiXYxs8XAD4D7gMnAeOAW4LKk034Tft6jgYuAvcB6Mzs15XJj3b08/JqR+eglHyg5Sjb8K/B5Mxvb04tmdo6ZPWtmu8Pv5yS9VmdmXzGzX4c9hsfMrCJ8LdETuNbMmsxsh5ndmPTe2Wb2GzP7c/jaf5pZJOl1N7PlZrYR2Bgeu9PMtprZG2a23szmJZ3/ZTP7gZl9N4yl3sxOMbMvmllL+L6Lk84/2sxWhz97u5l91cxKzOwdwF3A2WHv5c/h+WvM7Jtm9oiZtQPnh8e+mnTNhWb2fBjfK2Z2SX+N7+7bgf8FTg2vMdHMas3sdTPbZGafSPmMPww/4xvAUjMbZWZ3hG3cFD4eFZ5fYWY/Ddv4dTNbZ2YjwtdeNbOLwhj/HvhgordmZn9jZutTfg9uNLOf9PD7YcD/A77i7t92993u3u3uT7r7J1LPd/cD7v6Ku38KeBL4cn9tJKLkKNnwHFAHfD71BTM7FngY+A/gOOL/CD5sZsclnfZh4GPAOCDSw3XOB04GLgZuMrOLwuMHgM8CFcDZwIXAp1Le+wFgDvDO8PmzwOnAscB/Az8ws9FJ518G3A8cA2wAHiX+/9Uk4DZgVdK59wJdwHRgZhjfNe7+B+CThD0dd0/+o+HDwD8BRwGHDLua2WziPae/A8YC5wKv0g8zmwL8dRgvwPeBbcBEYDHwz2Z2YdJbFgI/DH/G94AvAX8ZtssMYDZv9UJvDK9VSbw39/fAITUq3f1nwD8D/5PUW6sFpoZ/KCQsId62qaqBKWFMA/UjYF6/Z0nRU3KUbLkF+LSZVaYcXwBsdPf73b3L3b8P/JFDh8u+4+4vu/teYC3xf6ST3eru7e5eD3wH+BCAu69392fC675KPHGdl/Ler7n76+G1cffvuntb+J7bgVHE/3FOWOfuj7p7F/FhvkrgX9y9E3gAONHMxprZeOBSYEUYWwvw78CV/bTTg+7+67BntC/ltWXAPe7+8/D17e7+xz6u9ZOwV/or4j2ofw4T5XuA/+Pu+9z9eeDbwFVJ7/uNu/8k/Bl7gY8At7l7i7vvBG5NOr8TOB44wd073X2dp1HA2d33A/9DPCFiZu8CTgR+2sPpiT+UdvR33R40Ef9DJ1lr2NP9s5m97Q82KU5KjpIV7v574v/w3ZTy0kTgtZRjrxHviSU0Jz1+E0idWLI15b0TAcIhz5+aWXM4RPjPxHuRvb03MbT3h3CI98/E718lv+dPSY/3Aq3ufiDpOWF8JwClwI7EP8TEk/M4+ra1j9emAK/08/5kH3D3se5+grt/Kkx0E4HX3X1P0nmp7Z0aQ+p/o4NtTHzIfBPwmJk1mlnqf9++3At8OBw2vQpYGybNVG3h9+MHcO2EScDrKccqwnYZ6+7/NohrSgFScpRs+kfgExz6D3ET8USSrArYPoDrTkl5b1P4+JvEe6Enu/sY4kN+lvLeg72c8P7i/wGuAI4Jhzt39/CedGwF9nPoP8Rj3P1dqT+3t3h6uea0QcSSrAk41syOSjqW2t6pMaT+NzrYxu6+x91vdPeTiPf2P5cyRNvbNXH3Z4AO4sOeH6bnIVWABuKffVFvH6oPlwPrBvE+KTJKjpI17r6J+FDaZ5IOPwKcYmYfNrORZvZB4vf/ehpe680/mNmR4dDcx8KfAfH7dm8AMTP7C+D6fq5zFPF7hDuBkWZ2CzBmAHEc5O47gMeA281sjJmNMLNpZpYY1v0TMDl5glAaVgMfM7MLw+tNCj/XQOLaCjwNfM3MRpvZu4kP136vj7d9H7jZzCotPhnqFuC7AGb2PjObHvb+3iB+n/dAD9f4E/Eh59R/g+4D/hPocvcel7aEw7SfI/7f+WNJ7fkeM7s79fxw0tNUM/sGMJ/4MLBIn5QcJdtuAw6ueXT3NuB9xCd2tAFfAN7n7q0DuOaTxIf2fgH8m7s/Fh7/PPEeyR7gW7yVNHvzKPFZnS8THzrcR9/DnP35KPEJRC8Bu4hPKEkMDf4SeBFoNrO0Pqu7/5Z48v934j3aJ3l7rzsdHyJ+f68J+DHwj+7+8z7O/yrxSVUvAPVANDwG8YlQjwMx4DfAfyXWNqb4Qfi9zcyiScfvJz6LtrdeIwDu/kPgg8DHw7j/FMbwYNJpZ5tZjHiSriP+h81Z4b1okT6ZNjuWQmFmJwKbgdJwgozkGTM7AmgBZrn7xmzHI8VLPUcRySXXA88qMUq2ZSw5mtk9Fl8I/fteXjcz+49w0fELZjYrU7GISO4zs1eJV066sZ9TRTIukz3HNUBf1TouJX5/4mTgWuIzCUUGzd1fdXfTkGp+cvcTw2UmG/o/WySzMpYc3f0p3r6eKNlC4D6PewYYa2aDWbckIiJyWGWziPAkDp35ty089raqF2Z2LfHeJaNHjz6jqqpqWAIsNN3d3YwYodvMA6V2Gzy13eBkut26u7sP/oz+fk5XVxfW1UWJGQfc8ZEjGTkyJ+rP9+jll19udffUylsDls1P2NNC6h6nzrr73cDdANXV1d7Q0JDJuApWXV0d8+fPz3YYeUftNnhqu8HJZLsFQcCq5ctZQLyI8XUrVzJnzpzDdn62mVlqha1ByWZy3MahlUwm81YlExERyYD6aJTzYjHOKi0l1tlJfTTaZ7KbM2cOrFxJfTTKdbNm5XRiPJyymRxrgRvM7AHiuyDsDquIiIhIhlgkwj2vvcabBw7wQEkJH430X5Rpzpw5RZMUEzKWHM3s+8RLNVWY2TbidTRLAdz9LuJlwv6aeCWTN4lX+hARkQxqbGhghjvtJSXMcKdRt6l6lLHk6O4f6ud1B5Zn6ueLiMjbbW9qorW7m7mRCL925x39v6UoaRqZiEiRCIKAtvXr2TNiBD/q6OCNykoWLhrM5iaFL3fn44qIyGFVH42y2J2xkyZRG4sx7Yoriu5eYrrUcxQRKRIWibBm61Z2NjfTGIsxrbo62yHlLPUcRUTyUBAE1EejnDaA5RXe0cHFVVUwciQXd3XhHR0ZjjJ/qecoIpJngiDg9qVL6fra17h96VKCIEjrfafNmsXGsjKOGTWKjWVlnDZL+z30Rj1HEZE8U1tTw+wtW1hYWkqsrY3ampq0eo/FuqB/MJQcRURyVG9Dpw78Gpjqzq9hQMsxinFB/2BoWFVEJAclapoes3o1q5YvP2TodOGiRXRVVfGLo4+mq6pKyzEyQD1HEZEs66mHWB+NsgBYdPzxsGPHITVQ58yZw81r1lAfjXK1hkczQslRRCSLkne9WLV6NaxcCcQnz6xavRp27IjvhpEyeUbDo5ml5CgikkU97ZIx/R3v0OSZLFNyFBHJos3btvHU5s3sM+OBkhKWJO2Sod5h9mhCjohIlgRBwLr772eeGU1mzKuo0ML8HKHkKCKSJfXRKPPNeM2Mie5sAC3MzxFKjiIiWWKRCOtaWznBnTu7u3n3kiUaRs0RuucoIpIljQ0NzC4vZ3x5OVdFIlROnpztkCSknqOISBYEQcALDz3Ey7EYnc3NGlLNMeo5iohkQW1NDQva25k2fTo/3rWL6ZddpiHVHKLkKCKSQT1Vv0n0Gl9qa2NuWxutVVVcrRJwOUXDqiIiGdJbfdT6aJSlZWXcNHMmjccdp15jDlJyFBHJkIPVb/bv57xYjPpoFIjfW3wY2LZvH/vGjVPh8BykYVURkQyxSIQ1W7fSCXwPDla/UWm43KfkKCKSId7RwbyKCto6O5lXWnpI9RuVhsttGlYVEcmQzdu28URzM2Nff511ra1YUt1UyW1KjiIiGaC6qflNyVFEJANUNzW/KTmKiGSA6qbmN03IERE5zIIg4FePP868igqqxozhqq4u1U3NM+o5iogcRomF/6c8/zzrWlsp7exkY1mZhlTzjHqOIiKHUW1NDbNbWvhgVRUAdaefzidXrNCQap5RchQRGaJE/VSLRA7WTI21tbGhqooblRjzkpKjiMgQJIZRFwB3trSw7Mgj+YuZM7l3yxbVTM1juucoIjIEiWHUC0aP5tKyMn64d69qphYA9RxFRAYpeeupxDDqB265hV0dHaqZmueUHEVEBiEIAu664w4Wux8yjLps2bJshyaHgZKjiMgAHVyu0d7Omq1bWQoaRi0wSo4iIgOk5RqFT8lRRKQfiaUaiYX8Wq5R+JQcRUT6kLxUY9Xq1Rx/wQUsLStjspZrFDQt5RAR6UPyUo0FgAMPg5ZrFDj1HEVEepG6VOO3VVXcuGgRLFpEfTSq5RoFTMlRRKQX9dEoi90ZO2ECtbHYIUOoSoqFTcOqIiK9sEiENVu3srO5mcZYjGnV1dkOSYaJeo4iIr3wjg7mVVTQ1tnJvNJSvKMj2yHJMFHPUUSkF5u3beOJ5mbGvv4661pbsUgk2yHJMFFyFBHpQRAErLv/fuaZ0WTGvIoK9RyLiJKjiEgPamtqOHv/fraNHMmUESPYAAeLAEjhU3IUEUmRWMLxcizG2K4uvnXUUXzglls0Q7WIaEKOiEiK1CUcF370o9pto8io5ygikkJLOEQ9RxGRFI0NDcwuL6e9vJyLIxFNxClC6jmKiCRZvXo1T953H3/YvZvO5mZNxClS6jmKiISCIODeW29l6Z49jBo5krqyMu26UaTUcxQRCdVHo8w344nubt7s6qJl9GjtulGklBxFREIWibCutZUT3Lmzu5t3L1miXmOR0rCqiEgoUUu1rLOTxaWlVE2enO2QJEvUcxQRCamWqiQoOYqIoFqqciglRxERVEtVDqXkKCJFT7VUJVVGk6OZXWJmDWa2ycxu6uH1KjN7wsw2mNkLZvbXmYxHRKQniVqqSydMoP3oo1VLVTKXHM2sBFgJXAq8E/iQmb0z5bSbgbXuPhO4EvivTMUjItIb1VKVVJlcyjEb2OTujQBm9gCwEHgp6RwHxoSPjwaaMhiPiEiPEks42jo7mVdaqok4grl7Zi5sthi4xN2vCZ9fBcxx9xuSzjkeeAw4BigDLnL39T1c61rgWoDKysoz1q5dm5GYC10sFqO8vDzbYeQdtdvg5UvbNTU1sae5meOA1804dsoUKioqshZPvrRbLjr//PPXu/uZQ71OJnuO1sOx1Ez8IWCNu99uZmcD95vZqe7efcib3O8G7gaorq72+fPnZyLegldXV4fabuDUboOXD20XBAF33XQT72lpockMHzeOETfdxOLFi7MWUz60W6HL5IScbcCUpOeTefuw6TJgLYC7/wYYDWTvzzURKTpawiE9yWRyfBY42cymmlmE+ISb2pRztgAXApjZO4gnx50ZjElE5CAt4ZDeZGxY1d27zOwG4FGgBLjH3V80s9uA59y9FrgR+JaZfZb4kOtSz9RNUBGRFIklHGMnTKA2FtMSDjkoo4XH3f0R4JGUY7ckPX4JmJvJGEREepNYwvERoBE4R0s4JKRdOUSkaGkJh/RG5eNEpGhZJMKTLS2MaW3lyZYW7cIhByk5ikjRamxoYIY77SUlzHCnsaEh2yFJjlByFJGitb2picbubqaUlLC1pORtC7GleCk5ikhRCoKAtvXr2TNiBD/q6OCNykoWLlqU7bAkR2hCjogUpYPLOCZNojYWY9oVV2h9oxyknqOIFCXtxCF9Uc9RRIqSlnFIX5QcRaSoBEFAfTTK5m3beLq1lY8APwOWaBmHJFFyFJGiEQQBq5YvZwFwX0sL8ysqYMwYLu7qUs9RDqF7jiJSNOqjUc6LxThr/37mm7EBOGbUKDaWlWknDjmEeo4iUjQSk3A6gXXAOStWsGvyZK6bNUszVeUQSo4iUhSCIOBXjz/OvKSh1MrJk7nm+uuzHZrkIA2rikjBS9xrPOX551nX2kppZ6eGUqVP6jmKSEELgoC77riD+bEYV598MgB1p5/OJ1es0FCq9ErJUUQK1urVq/nJbbcxE1jT2grAxvJyJUbpl5KjiBSkIAi499Zb+WhbG2eWlkJFBXUzZyoxSlqUHEWkINVHo8w344nubg7s28eGo4/mZiVGSZMm5IhIQbJIhHWtrZzgzp3d3bx7yRIlRkmbeo4iUhASZeFOC9csekcHF1dVUTlyJFeFyzZE0qXkKCJ5LwgCbl+6lAva27m9rIwb16zhtFmzWFVWxinA06NGcYGWbcgAKDmKSN6rralh9pYtLCwtJdbWRm1NDf/09a/DypXUR6OqgCMDpuQoInnPgacOHGBSdzdPuXNqeHzOnDlKijIompAjInlvWnU124HHDhxge/hcZCjUcxSRvOcdHSwYP56yzk6maONiOQzUcxSRvJdYtnHc7t2sa23FtHGxDJF6jiKS9xobGphdXk57eTkXRyLqOcqQKTmKSN4KgoAHa2oI1q5lTCxGZSzGhqoqbtSyDRkiJUcRyUuJbahGt7RwaVsbM6ZP58e7djH9sss0Q1WGTMlRRPJSfTTKebEYY0tLWX3gAOW7drFv3DgWLlqU7dCkAGhCjojkJYtEWLN1Kzubm9kOBOeey3UrV6rXKIeFeo4ikpeSJ+EsjkSonDtXiVEOG/UcRSTvBEHACw89xMuxGJ3NzWwATtMkHDmM1HMUkbxTW1PDgvZ2pmkSjmSIeo4ikleCIOCZtWv5aXMzGxoaaC0r0yQcOezUcxSRvFJbU8OlO3fyjkiE+zo6OPKMM9RrlMNOPUcRySvbm5qo6+wkBnSUljJx4sRshyQFSMlRRPJGEAS0rV/PnhEj+FFHB29UVmpIVTJCw6oikjfqo1EWuzN20iRqYzGmXXGFhlQlI9RzFJG8kbzwvzEW076NkjHqOYpI3vCODuZVVNDW2ck87dsoGaSeo4jkDYtEeLKlhTGtrTzZ0qJ9GyVjlBxFJG80NjQww532khJmuNPY0JDtkKRAKTmKSN7Y3tREY3c3U0pK2FpSgmc7IClY/SZHM3vOzJab2THDEZCISE+0jEOGUzo9xyuBicCzZvaAmb3XzCzDcYmIHCJRT/Xm6mqOmzCBOVrGIRnUb3J0903u/iXgFOC/gXuALWZ2q5kdm+kARUQSu3A82tbG7zZtUj1Vybi07jma2buB24F/BWqAxcAbwC8zF5qISFyi1/ip6dNpPO447cIhGdfvOkczWw/8GVgN3OTu+8OXAjObm8ngRKS4BUHAgzU1BGvX8lJbG3Pb2mitquJq9Rolw9IpAvA37t6YfMDMprr7Znf//zIUl4gUuSAIWLV8OaNbWri0rY0Z2rtRhlE6w6o/TPOYiMhhUx+Ncl4sxntLS3nqwAFe2bWLfePG6V6jDItee45m9hfAu4CjzSy5hzgGGJ3pwESkuCXqqH4E2A4E557LJz/9afUaZVj0NaxaDbwPGAtclnR8D/CJTAYlItLY0MDs8nLay8tZHIlQOXeuEqMMm16To7s/CDxoZme7+2+GMSYRKXKrV6/myfvu49jdu6mMxdhQVcWNs2ZlOywpIn0Nq37B3b8OfNjMPpT6urt/JqORiUhRCoKAe2+9laV79jBq5Ejqyso0CUeGXV/Dqn8Ivz83HIGIiEB8Is58M57o7uZcd1pGj+ZaTcKRYdbXsOpD4fd7hy8cESl2FomwrrWVs925s7uby5csUa9Rhl1fw6oPQe9F7939/RmJSESKWmIizvjycq6KRKicPDnbIUkR6mtY9d+GLQoREd6qoToyFtNEHMmqvoZVnxzOQERE6qNRFrszdsIEamMxTcSRrOm1Qo6ZrQ2/15vZC0lf9Wb2QjoXN7NLzKzBzDaZ2U29nHOFmb1kZi+a2X8P7mOISCGwSIR7XnuNpq1b2fTGG0yrrs52SFKk+hpW/dvw+/sGc2EzKwFWAn8FbCO+H2Stu7+UdM7JwBeBue6+y8zGDeZniUhhaGxoYIY77SUlzHCnsaEh2yFJkeq15+juO8LvrwH7gRnAu4H94bH+zAY2uXuju3cADwALU875BLDS3XeFP6tl4B9BRArF9qYmGru7mVJSwtaSkt5nBIpkmLn3/etnZtcAtxDfu9GA84Db3P2eft63GLjE3a8Jn18FzHH3G5LO+QnwMjAXKAG+7O4/6+Fa1wLXAlRWVp6xdu3atD+gvCUWi1FeXp7tMPKO2m3wBtJ27e3t7Hj1VQ7s308p0FVayqSTTqKsrCyzQeYg/c4N3vnnn7/e3c8c6nXS2bLq74CZ7t4GYGbHAU8DfSZH4ok0VWomHgmcDMwHJgPrzOxUd//zIW9yvxu4G6C6utrnz5+fRtiSqq6uDrXdwKndBm8gbfftb36T0m98g8rOTmpjMcZffTXXfvrTmQ0wR+l3LvvS2bJqG/Fi4wl7gK1pvm9K0vPJQFMP5zzo7p3uvhloIJ4sRaTIJHbh2NncTGMspsk4klV9FQH4XPhwOxCY2YPEe34Lgd+mce1ngZPNbGp4jSuBD6ec8xPgQ8AaM6sATgEaEZGik7wLx8WRCN7Rke2QpIj11XM8Kvx6hXgSSwyJPgjs6O/C7t4F3AA8SrxO61p3f9HMbjOzRHWdR4E2M3sJeAL4u8TwrYgUjyAIeGbtWl7ctYv927ezAThNi/8li/oqAnDrUC/u7o8Aj6QcuyXpsQOfC79EpEjV1tRw6c6dvCMS4b6ODo484wwt/pes6ndCjplVAl8A3gWMThx39wsyGJeIFBEHnjpwgElm7B0xgmkTJ2Y7JCly6UzI+R7wR2AqcCvwKvH7iSIih8W06mq2A48dOMD28LlINqWzlOM4d19tZn8b1lt90sxUd1VEDpvGhgYuOvpoJpaXc4om40gOSKfn2Bl+32FmC8xsJvFlGSIiQ5bYiePlWIzO5mZNxpGckE7P8atmdjRwI/ANYAzw2YxGJSJFo7amhgXt7UybPp0f79qlnTgkJ/SbHN39p+HD3cD5mQ1HRIpJYgnH75qbObelhdapU7l60aJshyXS/7CqmZ1kZg+ZWauZtZjZg2Z20nAEJyKFLbGE4/pIhPXd3VrCITkjnWHV/ya+9dTl4fMrge8D+g0WkUEJgoD6aJTtTU28BEwtKaGjtFRLOCRnpJMczd3vT3r+XTO7odezRUT6EAQBq5YvZwHQ1t7Om5WV/KK7m66yMhZqSFVyRF+1VY8NHz5hZjcR34/RgQ8CDw9DbCJSgGprapjd0sIFVVUArL/sMk6aOpWrZ83SkKrkjL56juuJJ8PE1lPXJb3mwFcyFZSIFKbEso2X2tqItbXx26oqbly0SElRck5ftVWnDmcgIlK4EvcYN2/ezNKyMibPnMm9W7Zo2YbkrHRqq5YC1wPnhofqgFXu3tnrm0REQkEQcPvSpVzQ3s4zI0aw8Ygj+CCwb9w43WOUnJXOhJxvAqXAf4XPrwqPXZOpoESkcNTW1DB7yxYWlpYS6+zk95dfzq65c7lO9xglh6WTHM9y9xlJz39pZr/LVEAiUlgO7rjR3c1T7pw6cSLXXH99tsMS6VM6tVUPmNm0xJOwAMCBzIUkIoVEO25IPkqn5/h3xJdzNBKfuXoC8LGMRiUiBcM7OlgwfjxlnZ1MKS3VjhuSF/pMjmY2AtgLnAxUE0+Of3T3/cMQm4gUgM3btvFUczNLzPhZSQlLIpFshyTSrz6To7t3m9nt7n428MIwxSQiBaK9vZ1199/PPDOazJhXUaGeo+SFdO45PmZmi8zM+j9VROQte998k/lmvGbGRHft1Sh5I517jp8DyoAuM9tHfGjV3X1MRiMTkfxnxrrWVs52587ubi5fskTLNyQvpLOf41HDEYiIFJYgCIjt2cO8igqqxozhqq4uKidPznZYImnpdVjVzMaZ2R1m9lMz+2czU09RRNKS2Hlj9Jtvsq61ldLOTjaWlWlIVfJGX/cc7wPagW8ARwH/MSwRiUjeq49GOS8W49hwEk7dzJlct3KlhlQlb/SVHCe4+5fc/VF3/zTw7uEKSkTym0Ui3PPaa+zet48nW1p4z0UXKTFKXunrnqOZ2TG8tWVVSfJzd38908GJSH5qbGhghjvdwAx3Ghsash2SyID0lRyPJr6nY/ISjmj43YGTMhWUiOQ3B14F5offT81iLCKD0dd+jicOYxwiUkCmVVfzKPAGsB1YqHqqkmfSKQIgIjIgiXqqo0eOZMH48aqKI3lHyVFEDjuLRFjX2srIAwdY19qKqZ6q5BklRxE57BobGphdXk53aSkXV1Wp5yh5J63kaGbvMbOPhY8rzWxqZsMSkXwVBAHPrF3Li7t24R0dqqcqeanf8nFm9o/AmcS3rPoOUAp8F5ib2dBEJB/V1tRw6c6dvCMS4U/AkWecoTWOknfS6TleDryfeLUc3L2JeMUcEZG32d7URF1nJzHAzZg4cWK2QxIZsHSSY4e7O/GlS5hZWWZDEpF8FQQBbevXs2fECH7U0cGBkSNZuGhRtsMSGbB0tqxaa2argLFm9gng48C3MhuWiOSj+miUxe6MnTSJ2liMsmOP1ZCq5KV0tqz6NzP7K+LreauBW9z95xmPTETyjkUirNm6lY8AjcClo0ZlOySRQUlnQs5ngR8oIYpIfxJLONrLy7k4EgH3bIckMijp3HMcAzxqZuvMbLmZjc90UCKSf1avXs2T993HH8bdkjQAABxhSURBVHbvprO5mQ3AEUceme2wRAal3+To7re6+7uA5cBE4EkzezzjkYlI3giCgHtvvZWle/bwNyNH8ofycqZfdhllZZq/J/lpIBVyWoBmoA0Yl5lwRCQf1UejzDfjie5u3uzqomX0aM1SlbzWb3I0s+vNrA74BVABfMLdtfGxiByUqKV6gjt3dnfz7iVLNEtV8lo6SzlOAFa4+/OZDkZE8lNiIs748nKuikSonDw52yGJDEmvPUczGxM+/DqwxcyOTf4anvBEJNcl11Ldv327aqlKQeir5/jfwPuA9cSr41jSaw6clMG4RCRPJNdSva+jQ7VUpSD0mhzd/X3hd+3AISK9cuCpAweYZMbeESOYplqqUgDSmZDzi3SOiUhxmlZdzXbgsQMH2B4+F8l3vfYczWw0cCRQYWbH8Naw6hji6x1FpMgFQcCvHn+cBePHUzVmDKd0dWljYykIfd1zvA5YQTwRruet5PgGsDLDcYlIjguCgFXLl3NKezuPtbay9Igj2FhezgWajCMFoK97jncCd5rZp939G8MYk4jkgdqaGma3tPDBqioA6k4/nU+uWKHJOFIQ0tmV4xtmdirwTmB00vH7MhmYiOSmIAh4sKaGYO1aXmprI9bWxoaqKm5UYpQCks6uHP8IzCeeHB8BLgV+BSg5ihSZxFDq6JYWLm1rY8b06fx41y6mX3aZEqMUlHRqqy4GLgSa3f1jwAxAm7SJFKH6aJTzYjHeW1rKUwcO8MquXewbN051VKXgpJMc97p7N9AVVs1pQQUARIpSYjPjnc3NbAeCc8/lupUr1WuUgpNObdXnzGws8C3is1ZjwG8zGpWI5CTv6GBeRQVtnZ0sKC2lau5cJUYpSOlMyPlU+PAuM/sZMMbdX8hsWCKSixK7b3wE+BmwJBLJdkgiGdFXEYBeFyuZ2Sx3j2YmJBHJRYkF//MqKmDMGC7Wgn8pYH31HG/v4zUHLjjMsYhIjgqCgNuXLmXq66/z5Btv8HEt+JcC11cRgPOHMxARyV21NTXM3rKFj5SW8jV3fjhlCjd/9au63ygFK511jh/t6biKAIgUj4M7b3R38ypw2syZSoxS0NKZrXpW0uPRxNc8RlERAJGiMa26mkcJd94YMYKF2nlDClw6s1U/nfzczI4G7k/n4mZ2CXAnUAJ8293/pZfzFgM/AM5y9+fSubaIDJ/GhgYuOvpoJpaXc0okook4UvDSKQKQ6k3g5P5OMrMS4rt3XEq89NyHzOydPZx3FPAZIBhELCKSYUEQ8MJDD/FyLEZnczMbgNM0EUcKXDr3HB8ifssB4sn0ncDaNK49G9jk7o3hdR4AFgIvpZz3FeDrwOfTjFlEhlF9NMpid8ZOmEBtLKY6qlIU0rnn+G9Jj7uA19x9WxrvmwRsTXq+DTjk/ygzmwlMcfefmlmvydHMrgWuBaisrKSuri6NHy+pYrGY2m4Qir3dxlZW0nbttXQQ/4v3uKqqtNuj2NtusNRu2ZfOPccnAcK6qiPDx8e6++v9vNV6OOYHXzQbAfw7sDSNGO4G7gaorq72+fPn9/cW6UFdXR1qu4Er9nb70he+wJ7/+A+qzNjizlGf+Qz/9PWvp/XeYm+7wVK7ZV+/9xzN7Foz+xPwAvAc8fqq6Uya2QZMSXo+GWhKen4UcCpQZ2avAn8J1JrZmemFLiLDYXtTE43d3UwpKWFrSclbf+GKFLB0hlX/DniXu7cO8NrPAieb2VRgO3Al8OHEi+6+G6hIPDezOuDzmq0qkjuCIKBt/Xr2jBjBjzo6eGPiRG1PJUUhneT4CvEZqgPi7l1mdgPwKPGlHPe4+4tmdhvwnLvXDvSaIjK8Dk7GmTSJ2liMaVdcock4UhTSSY5fBJ42swDYnzjo7p/p743u/gjwSMqxW3o5d34asYjIMErs3/gRoBE4R4v/pUikkxxXAb8E6oHuzIYjIrkkef/GeaWlWvwvRSOd5Njl7p/LeCQikjOCIKA+GmXztm08rf0bpQilkxyfCNcZPsShw6r9LeUQkTwUBAGrli9nAXBfSwvztX+jFKF0kmNihukXk445cNLhD0dEsq0+GuW8WIyzSkuZb8YGYMaoUTw9apT2b5SikU4RgKnDEYiI5IbEJJxOYB1wzooV7Jo8metmzdJMVSka2s9RRA6ROgmnavJkrrn++myHJTKstJ+jiBxi87ZtPNXczBIzflZSokk4UpQyup+jiOSXIAhYd//9zDOjyYx5FRWahCNFKWP7OYpI/qmtqeHs/fvZNnIkU0aM0N6NUrQyuZ+jiOSRxKbGI2Mxphw4wLfGjuXaW27RJBwpSpncz1FE8khtTQ0L2tuZNn06P961iwuvvJJly5ZlOyyRrOg1OZrZdGB8Yj/HpOPzzGyUu7+S8ehEZFgkeo0vtbUxt62N1qoqrtbuG1LE+rrneAewp4fje8PXRKRAJHbfWDphAn8oL2f6ZZdpOFWKWl/J8UR3fyH1YLjf4okZi0hEhp1FItzz2ms0bd3KpjfeYJp235Ai11dyHN3Ha0cc7kBEJHsaGxqY4U57SQkz3GlsaMh2SCJZ1VdyfNbMPpF60MyWAeszF5KIDDcHXgUmhd+9r5NFikBfs1VXAD82s4/wVjI8E4gAl2c6MBEZPtOqq3kUeOzAAbaPGMFCDatKkes1Obr7n4BzzOx84NTw8MPu/sthiUxEhk1jQwMXHX00E8vLOSUSUVUcKXrplI97AnhiGGIRkSxIXvxfGYuxoaqKG1UVR4pcOkUARKSAJZZxjJ0wgdpYTMs4RBhcbVURKSBaxiHydkqOIkVOyzhE3k7JUaTIaRmHyNspOYoUuWnV1WwnXMYRPhcpdpqQI1LkvKODBePHU9bZyZTSUi3jEEE9R5GiZ5EI61pbOW73bta1tmKRSLZDEsk69RxFipx3dDCvooK2zk7mqecoAqjnKFL0Nm/bxhPNzYx9/XX1HEVCSo4iRSwIAtbdfz/zzGgyY15FhXqOIig5ihS12poazt6/n20jRzJlxAg2AKepdJyIkqNIsUrUVH05FmNsVxffOuooPnDLLSodJ4Im5IgUrdqaGha0tzNt+nR+vGsXF155JcuWLct2WCI5QT1HkSKU6DU+2tbG7zZtorWsjIWLFmU7LJGcoZ6jSBHSThwifVPPUaQIaScOkb4pOYoUIe3EIdI3JUeRIhMEAb97/nleBaZGImwtKdFOHCIpdM9RpIgEQcCq5cuZ297OD4FHRo+m69hjNRlHJIWSo0gRqY9GOS8W4/zSUvaOH8/Wc87h5hUrNBlHJIWSo0gR2bxtG09t3sw+M9aVlLDkoouUGEV6oHuOIkVCdVRF0qfkKFIkVEdVJH1KjiJFQHVURQZG9xxFikBqRZwLP/pR1VEV6YN6jiJFwCIR1mzdys7mZhpjMVXEEemHeo4iRcA7OphXUUFbZyfzSks1EUekH+o5ihSBzdu28URzM2Nff511ra1YJJLtkERympKjSIHTEg6RgVNyFClw9dEo8814zYyJ7lrCIZIGJUeRAmeRCOtaWznBnTu7u3n3kiVawiHSD03IESlwick4ZZ2dLC4tpWry5GyHJJLz1HMUKXCJnuNxu3drMo5ImtRzFClgQRDwq8cfZ15FBYwZw8VdXZqMI5IG9RxFCtTq1av56hVXMOXpp1nX2kppZycby8o0GUckDeo5ihSgIAi499Zb+WhbG2eWlkJFBXUzZ/JJ7d0okhYlR5EClNiB4ykz6OxkA2hTY5EB0LCqSIHRDhwiQ6fkKFJAgiDgrjvuYLE7N82cSde4cdqBQ2QQNKwqUiCCIGDV8uWc0t7Omq1bWQrsGzeOhYsWZTs0kbyj5CiS54IgoD4aZfPmzSwAFk2fDkDd6adrAo7IIGU0OZrZJcCdQAnwbXf/l5TXPwdcA3QBO4GPu/trmYxJpJCsXr2an9x2G4uPOIIXzNgIsGMHG8vKlBhFhiBjydHMSoCVwF8B24BnzazW3V9KOm0DcKa7v2lm1wNfBz6YqZhECknyco0ZpaUwYQIN738/u6ZO5bpZs5QYRYYgkz3H2cAmd28EMLMHgIXAweTo7k8knf8MsCSD8YgUjMTEm/lmPF1aCp2dPLh3LzcvWqSkKHIYmLtn5sJmi4FL3P2a8PlVwBx3v6GX8/8TaHb3r/bw2rXAtQCVlZVnrF27NiMxF7pYLEZ5eXm2w8g7udZu7e3ttG7Zwqjubt7o6KB85Ej2AMcdfzwVFRXZDu8QudZ2+ULtNnjnn3/+enc/c6jXyWTP0Xo41mMmNrMlwJnAeT297u53A3cDVFdX+/z58w9TiMWlrq4Otd3A5VK7BUHAvatWMX/DBq48+WS+tmkTv8nhiTe51Hb5RO2WfZlMjtuAKUnPJwNNqSeZ2UXAl4Dz3H1/BuMRyWupSzUANpaX52xiFMlnmUyOzwInm9lUYDtwJfDh5BPMbCawivjwa0sGYxHJe/XRKOfFYpxfWspe1UoVyaiMJUd37zKzG4BHiS/luMfdXzSz24Dn3L0W+FegHPiBmQFscff3ZyomkXy2eds2ntq8mX1mrCspYclFFykximRIRtc5uvsjwCMpx25JenxRJn++SKEIgoB199/PPDOazJhXUaF9GUUySLVVRXJc8rKN5tGjmTJiBBtA+zKKZJDKx4nksORJOI+1tjKvooL73Llau2yIZJSSo0gO6q1e6sunn86/ahKOSMYpOYrkmCAIuH3pUi5ob+eZESPYeMQRqpcqMsyUHEVyTG1NDbO3bGFhaSmxzk5+f/nl7Jo7V/VSRYaRkqNIjtne1MTvOjuZYsavgXdMnMg111+f7bBEiopmq4rkkCAIaFu/nj0jRvCjjg7eqKzUZsUiWaCeo0gOSJ6As7SsjMlnnsm9W7Yw7YorNJQqkgVKjiJZlliusQB4sL2djcQ3Nd03bpx6jSJZouQokiVvW65x/PGwYwfrL7hAGxaLZJmSo0gW9LZc42HgOm1YLJJ1So4iWaDlGiK5TclRJAu0XEMkt2kph8gw03INkdynnqNIhiUm3pwWDpnWR6NariGS49RzFMmgxDKNY1avZtXy5QRBwGmzZvEwsG3fPi3XEMlR6jmKZFB9NMp5sRhnhRNv6qPR+L3FlSupj0Y1AUckRyk5imRIEASs+/WvefW11+gcMYLvAUsiEQDmzJmjpCiSw5QcRTIgMZxa1tLCDHfaJ0zg4kgE7+jIdmgikgbdcxTJgMRw6ntLS3kVGNXZycayMk6bNSvboYlIGpQcRTLAIhHWbN3KzuZmtgPBuedy3cqVGkoVyRMaVhU5zIIg4FePP868igoYM4bFXV1Uzp2rxCiSR9RzFDmMEvcaT3n+eda1tlKq4VSRvKSeo8hh8LYdNqZPB6Du9NP55IoV6jWK5BklR5Eh6m2HjY1lZUqMInlKyVFkiLTDhkjhUXIUGSIHnjpwgEnd3TzlzqnaYUMk72lCjsgQTauuZjvw2IEDbA+fi0h+U89RZIi8o4MF48dT1tnJlNJSVcERKQDqOYoMkUUirGtt5bjdu1nX2oqF9VNFJH+p5ygyQKn7M3pHBxdXVcHIkVzc1aWeo0gBUHIUGYDVq1fzk9tuY/ERR7CqvBxWruS0WbNYVVbGKcDTo0ZxgRb8i+Q9JUeRNAVBwL233spH29qYUVoKEyZof0aRAqXkKJKm2poazt6/n6fMoLOTB/fu5eawl6j9GUUKiybkiKQhCAJeeOghXo7FGNvVxbeOOooP3HKLEqJIgVLPUaQXiYk3x1dVsekPf2BpWRmTZ87k3i1buPDKK1m2bFm2QxSRDFFyFAklz0IFWLV8OQuA1mXLsEiEh4EF+/axb9w4Fi5alNVYRSSzlBxFeGurqQXAqtWrOf6CC+K7axx/PD8BWjs6uE6TbkSKhpKjCFAfjXJeLMZZYfHwBuBhgB07iMHBNY1KiiLFQclRhHiVmzVbt9IJfA9YUl3NwkWLDt5zVFIUKS5KjlKU0qlyk+gp1tXVZTtcERlmWsohRSdxf/GY1atZtXw5QRBw2qxZbCwr45hRo9hYVnZwUo6IFCf1HKVoJHqLmzdvPjjZhh07VOVGRN5GyVEKRupQaepridmoD7a3sxFgxw4eBq5TlRsRSaHkKAUhdSkGK1cekujqo9FDeovrL7iAXVOnqpcoIj3SPUfJe0EQcNcdd3BeLMai449nAfFkmOy0WbN4GKgJe4sLFy3imuuvV2IUkR6p5yh5LdFjPKW9nTVbtwLwZHn5waHShDlz5uieooikTclR8k7yvcWDw6XTpwNQd/rpfHLFih6Tn+4piki6NKwqeSV1GUai5mnNjh1sLCvrNTGKiAyEeo6SV2prapjd0sIFVVWwbx+7VPNURDJAyVHyRhAEPLN2Lb9rbuaNlhaemzqVG1XzVEQyQMlR8kZtTQ2X7tzJOyIR7uvo4MgzzlBSFJGMUHKUnJeYgLO9qYmXgKklJXSUljJt4sRshyYiBUrJUXJa8uL+tvZ23qys5Bfd3XSVlWnDYRHJGCVHyUm91UFdf9llnDR1Kldr8o2IZJCSo+ScPuugLlqkpCgiGafkKFnRV5Fw1UEVkWxTcpTDrq/El3i9ryLhp82aFT+u3qKIZImSoxxW/SU+eHvPsD4aPeQc1UEVkWxT+ThJWxAEfPub3yQIgl7PSU58Pe2OAW/fIeO0lCLhEE+Q2jVDRLIlo8nRzC4xswYz22RmN/Xw+igz+5/w9cDMTsxkPDJ4qTVNe0uQ6Sa+61auZNeyZVzXQ89SRCTbMjasamYlwErgr4BtwLNmVuvuLyWdtgzY5e7TzexK4P8CH8xUTMUsCAJad+4kCIJBJaP+hkIT0h0SVck3Ecllmew5zgY2uXuju3cADwALU85ZCNwbPv4hcKGZWQZjKkqJXt/I1tY+e319SadHmKAhURHJd+bumbmw2WLgEne/Jnx+FTDH3W9IOuf34TnbwuevhOe0plzrWuDa8OmpwO8zEnSBMqicABUHwErAm6HVYecgLlVmcKTDm0D74Y4zh1UArf2eJT1R2w2O2m3wqt39qKFeJJOzVXvqAaZm4nTOwd3vBu4GMLPn3P3MoYdXfNR2g6N2Gzy13eCo3QbPzJ47HNfJ5LDqNmBK0vPJQFNv55jZSOBo4PUMxiQiItKvTCbHZ4GTzWyqmUWAK4HalHNqgavDx4uBX3qmxnlFRETSlLFhVXfvMrMbgEeBEuAed3/RzG4DnnP3WmA1cL+ZbSLeY7wyjUvfnamYi4DabnDUboOnthsctdvgHZa2y9iEHBERkXylCjkiIiIplBxFRERS5FRyHEq5OTP7Yni8wczeO5xxZ9tg283MTjSzvWb2fPh113DHnm1ptN25ZhY1s65w7W7ya1eb2cbw6+rU9xayIbbbgaTfudRJegUvjbb7nJm9ZGYvmNkvzOyEpNf0Oze4dhv475y758QX8Uk7rwAnARHgd8A7U875FHBX+PhK4H/Cx+8Mzx8FTA2vU5Ltz5QH7XYi8Ptsf4Ycb7sTgXcD9wGLk44fCzSG348JHx+T7c+U6+0WvhbL9mfI8bY7HzgyfHx90v+v+p0bRLuFzwf8O5dLPcehlJtbCDzg7vvdfTOwKbxeMVCZvsHrt+3c/VV3fwHoTnnve4Gfu/vr7r4L+DlwyXAEnQOG0m7FLp22e8Ld3wyfPkN8jTjod26w7TYouZQcJwFbk55vC4/1eI67dwG7gePSfG+hGkq7AUw1sw1m9qSZzct0sDlmKL83+p17y0A/+2gze87MnjGzDxze0HLeQNtuGfC/g3xvIRlKu8EgfudyabPjoZSbS6sMXYEaSrvtAKrcvc3MzgB+Ymbvcvc3DneQOWoovzf6nTvUQD57lbs3mdlJwC/NrN7dXzlMseW6tNvOzJYAZwLnDfS9BWgo7QaD+J3LpZ7jUMrNpfPeQjXodguHodsA3H098TH9UzIece4Yyu+NfufeMqDP7u5N4fdGoA6YeTiDy3FptZ2ZXQR8CXi/u+8fyHsL1FDabXC/c9m+0Zp0w3Qk8RvMU3nrhuu7Us5ZzqETS9aGj9/FoRNyGimeCTlDabfKRDsRv9G9HTg2258pl9ou6dw1vH1CzmbiEyOOCR8XRdsNsd2OAUaFjyuAjaRMrCjkrzT/f51J/A/Vk1OO63ducO02qN+5rH/olA/x18DL4Qf8UnjsNuJ/BQCMBn5AfMLNb4GTkt77pfB9DcCl2f4s+dBuwCLgxfAXLQpclu3PkoNtdxbxv1rbgTbgxaT3fjxs003Ax7L9WfKh3YBzgPrwd64eWJbtz5KDbfc48Cfg+fCrVr9zg2+3wf7OqXyciIhIily65ygiIpITlBxFRERSKDmKiIikUHIUERFJoeQoIiKSQslRJEVSBf/fm9kPzOzIDP6s28KFy5jZioH+LIv7pZmNCZ9/xsz+YGbfOwyxLTWziUnPv21m7xzktW4ws48NNSaR4aKlHCIpzCzm7uXh4+8B6939/6XxPiP+/9Sgim2b2avAme7eOoD3LAAucvfPhs//SHyd7+aU80Z6vK7uQOKpAz7v7s8N5H29XOtI4NfuXkzVcCSPqeco0rd1wHQ4uF/c78OvFeGxE8Oe2n8RL6Qwxcw+ZGb14Xn/NzyvxMzWhMfqzSyRzNaY2WIz+wwwEXjCzJ4ws2Vm9u+JIMzsE2bWU4L+CPBgeM5dxCsd1ZrZZ83sy2Z2t5k9BtwXxrou3GcxambnJF3/C2FcvzOzf7H4HoxnAt8Le9FHmFmdmZ0Znv+2zxgej5nZP4XXecbMxgN4fLeEV82sWHbLkXyX7aoH+tJXrn0R7v1GvGTVg8T3hjuDeHWNMqCceGWhmcT3LewG/jJ8z0RgC/HSfCOBXwIfCN//86SfMTb8voawvBrwKlARPi4jXgmkNHz+NHBaD7G+BhyV9Dz5Gl8G1gNHhM+PBEaHj08GngsfXxpeP7EX3rHh9zriPVmSn/f2GcNznLDSEvB14Oak938JuDHb/331pa90vtRzFHm7I8zseeA54klgNfAe4Mfu3u7uMeBHQGKLr9fc/Znw8VlAnbvv9Pgw5veAc4nXhTzJzL5hZpcAfe584u7txJPO+8zsL4gnyfoeTj3W3ff0calad98bPi4FvmVm9cTLCSbuH14EfMfDvfDc/fW+YuvjMwJ0AD8NH68n/sdDQgvxxCqS83JpyyqRXLHX3U9PPtDP5tDtyaf2dIK77zKzGcQ3rF0OXEG8TmZfvg38PfBH4Du9nNNlZiO89/ucybF9lnjtyRnEb6nsS4p5IJMP+mqLTndPXOsAh/4bMxrY+/a3iOQe9RxF0vMU8AEzO9LMyoDLid+PTBUA55lZhZmVAB8CnjSzCmCEu9cA/wDM6uG9e4CjEk/cPSC+Tc+Hge/3ElcD8fuM6Tga2BEm0quAkvD4Y8DHEzNlzezYnuLp7zOm8fNPAX6fZqwiWaXkKJIGd48Svz/4W+LJ4dvuvqGH83YAXwSeINztxN0fJL5reV04XLsmPCfV3cD/mtkTScfWEp/luauX0B4G5qf5Mf4LuNrMniGeqNrDmH8G1ALPhfF9Pjx/DXBXYkJOGp+xP3OJ75wgkvO0lEMkh5nZT4F/d/df9PL68cB97v5XwxvZwJjZTOBz7n5VtmMRSYd6jiI5yMzGmtnLxO9/9pgY4WAv7luJIgA5rIL4cLJIXlDPUUREJIV6jiIiIimUHEVERFIoOYqIiKRQchQREUmh5CgiIpLi/wfQHMO0vwWPMQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAccAAAGWCAYAAAANLkjgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXiU5dXH8e/JJIGQELaEHURZUizIIouKWFyrVaQKRVSsWK1IqRa3vrZafbXWBaVgFRVoFHBHsArqW6nVAG4jEpS4IcgaQggJYZmQkO1+/5gJhhjCsEwmk/w+15UrM882Jw+ak3M/92LOOUREROQHUeEOQEREpK5RchQREalCyVFERKQKJUcREZEqlBxFRESqUHIUERGpQslRJIKY2dNm9pdwx3EkzOwrMxsW7jhEgqHkKLXKzDaY2TYzi6+07TozSwtjWGFnZuPM7INDHeecu8E599cjuP4GMys0M1/g/j9rZglHFu2Rcc791DmXFojnf83s+aO5npn93MyWmtkeM9tuZkvM7OLAvnFmVhb4eX1mtj7wM/eodH4XM3OVjvGZ2RdH9UNKvaHkKOEQDfwh3EHUxMyiwx1DVWbmOcpLDHfOJQD9gYHAXUcQQ524L2Y2CngVmAt0BNoAdwPDKx32ceDnbQacAxQCK8ysV5XLNXfOJQS++oQ+eokESo4SDo8At5lZ8+p2mtlpZrbczHYFvp9WaV+amf3VzD4MVAyLzSwpsK+iErjezLLMbKuZ3Vrp3EFm9rGZ7Qzse8LMYivtd2Y20czWAGsC2x4zs81mttvMVpjZ0ErH/6+ZvWpmzwdiyTCzHmb2JzPLCZx3XqXjm5lZauCzt5jZ/WbmMbOewNPAqYHqZWfg+Nlm9pSZvW1mBcCZgW33V7rmCDP7PBDf92Z2/qFuvnNuC/B/QK/ANdqb2UIz22Fma83st1V+xvmBn3E3MM7MGpnZtMA9zgq8bhQ4PsnM3gzc4x1mtszMogL7NpjZOYEY/wxcVlGtmdmvzGxFlf8ObjWz16v578OAvwN/dc790zm3yzlX7pxb4pz7bdXjnXNlzrnvnXO/A5YA/3uoeySi5Cjh8BmQBtxWdYeZtQTeAv4BtML/S/AtM2tV6bArgGuA1kBsNdc5E+gOnAfcYWbnBLaXATcDScCpwNnA76qc+0tgMHBi4P1yoC/QEngReNXMGlc6fjjwHNACWAm8g///qw7AfcCMSsfOAUqBbkC/QHzXOee+AW4gUOk45yr/0XAF8DegKXBAs6uZDcJfOd0ONAfOADZwCGbWCfhFIF6Al4BMoD0wCnjAzM6udMoIYH7gM14A7gROCdyXPsAgfqhCbw1cKxl/Nfdn4IA5Kp1z/wYeAF6pVK0tBI4P/KFQYSz+e1tVCtApENPheg0YesijpMFTcpRwuRu40cySq2y/EFjjnHvOOVfqnHsJ+JYDm8uedc5955wrBObh/yVd2b3OuQLnXAbwLHA5gHNuhXPuk8B1N+BPXD+rcu6DzrkdgWvjnHveOZcXOGcK0Aj/L+cKy5xz7zjnSvE38yUDDznnSoCXgS5m1tzM2gAXAJMCseUAU4Exh7hPbzjnPgxURkVV9l0LPOOc+09g/xbn3Lc1XOv1QFX6Af4K6oFAojwd+B/nXJFz7nPgn8BVlc772Dn3euAzCoErgfuccznOue3AvZWOLwHaAcc550qcc8tcEBM4O+f2Aa/gT4iY2U+BLsCb1Rxe8YfS1kNdtxpZ+P/QqSw3UOnuNLMf/cEmDZOSo4SFc+5L/L/47qiyqz2wscq2jfgrsQrZlV7vBap2LNlc5dz2AIEmzzfNLDvQRPgA/iryYOdWNO19E2ji3Yn/+VXlc7ZVel0I5Drnyiq9JxDfcUAMsLXiFzH+5Nyamm2uYV8n4PtDnF/ZL51zzZ1zxznnfhdIdO2BHc65PZWOq3q/q8ZQ9d9o/z3G32S+FlhsZuvMrOq/b03mAFcEmk2vAuYFkmZVeYHv7Q7j2hU6ADuqbEsK3JfmzrlHj+CaUg8pOUo43QP8lgN/EWfhTySVdQa2HMZ1O1U5Nyvw+in8VWh351wi/iY/q3Lu/ion8Hzxf4DRQItAc+euas4JxmZgHwf+Ik50zv206uceLJ6DXLPrEcRSWRbQ0syaVtpW9X5XjaHqv9H+e+yc2+Ocu9U5dwL+av+WKk20B7smzrlPgGL8zZ5XUH2TKsBq/D/7yIP9UDW4BFh2BOdJA6PkKGHjnFuLvyntpkqb3wZ6mNkVZhZtZpfhf/5XXfPawfzFzJoEmuauCXwG+J/b7QZ8ZvYTYMIhrtMU/zPC7UC0md0NJB5GHPs557YCi4EpZpZoZlFm1tXMKpp1twEdK3cQCkIqcI2ZnR24XofAz3U4cW0GPgIeNLPGZnYS/ubaF2o47SXgLjNLNn9nqLuB5wHM7CIz6xao/nbjf85bVs01tuFvcq76O2gu8ARQ6pyrdmhLoJn2Fvz/ztdUup+nm9nMqscHOj0db2aPA8PwNwOL1EjJUcLtPmD/mEfnXB5wEf6OHXnAH4GLnHO5h3HNJfib9v4LPOqcWxzYfhv+imQPMIsfkubBvIO/V+d3+JsOi6i5mfNQfo2/A9HXQD7+DiUVTYPvAV8B2WYW1M/qnPsUf/Kfir+iXcKPq+5gXI7/+V4W8C/gHufcf2o4/n78napWARlAemAb+DtCvQv4gI+BJyvGNlbxauB7npmlV9r+HP5etAerGgFwzs0HLgN+E4h7WyCGNyoddqqZ+fAn6TT8f9gMDDyLFqmRabFjqS/MrAuwHogJdJCRCGNmcUAO0N85tybc8UjDpcpRROqSCcByJUYJt5AlRzN7xvwDob88yH4zs38EBh2vMrP+oYpFROo+M9uAf+akWw9xqEjIhbJynA3UNFvHBfifT3QHrsffk1DkiDnnNjjnTE2qkck51yUwzGTloY8WCa2QJUfn3FJ+PJ6oshHAXOf3CdDczI5k3JKIiMgxFc5JhDtwYM+/zMC2H816YWbX468uady48cmdO3eulQDrm/LycqKi9Jj5cOm+HTnduyMT6vtWXl6+/zMO9TmlpaVYaSkeM8qcw0VHEx19dKmjrKzsR1/l5eX7v1d8lZWV4Zzb/77i9SE6kuY656rOvHXYwpkcqxtIXe1P7JybCcwESElJcatXrw5lXPVWWloaw4YNC3cYEUf37cjp3h2ZUN43r9fLjIkTuRD/JMbjp09n8ODBx+T4kpISNm3axIYNG9i4cSMbN24kMzOTrVu3snXrVrKyssjNzaW8vLza86OiomjWrBktWrQgMTGRpk2bkpCQQHx8/P7vcXFxB3w1btyYRo0a7f+67LLLqs6wdUTCmRwzOXAmk478MJOJiIiEQEZ6Oj/z+RgYE4OvpISM9PQak+PgwYNh+nQy0tMZ378/gwcPZteuXXz55ZdkZGTwzTff8N1337FmzRo2bNhAWdkPcz6YGW3atKFDhw507NiRQYMG0aZNG5KTk0lOTiYpKYlWrVrRokULWrZsSdOmTfHPH3HkLrvssqM6v0I4k+NC4Pdm9jL+VRB2BWYRERGRELHYWJ7ZuJG9ZWW87PHw69iaJ2XKz89n9+7dbNuxgzcffJD09HQ2b/7hiVhCQgLdu3fn5JNPZsyYMXTt2pUuXbpw3HHH0bFjR2IPcf26KmTJ0cxewj9VU5KZZeKfRzMGwDn3NP5pwn6BfyaTvfhn+hARkRBat3o1fZyjwOOhj3Osq/KYatu2bSxZsoQlS5aQlpbG119/vX9fjx49OP300+nduzcnnXQSvXv3plOnTkdd7dVFIUuOzrnLD7HfARND9fkiIvJjW7KyyC0vZ0hsLB86R0p5OR988AFvvfUWb7/9NqtWrQL8FeGQIUO44oorGDx4MAMGDKB582rXJ6+XwtmsKiIitcjr9ZK3YgW7zJheVMSaxo15b9YsHp4yBY/Hw+mnn85DDz3EsGHD6N+/PzExMeEOOWyUHEVEGgDnHK/Pn49v61a+KC1lT3k5jUpK+NUllzB8+HDOO++8BlUZHoqSo4hIPbZz506ef/55Zs6cSUZGBlHAyR4PhbGxTHzsMW644YZwh1gnKTmKiEQgr9dLRno6vQPDK6pau3YtU6ZMYc6cORQWFjJgwACuuuIKTli+nI6NGrG9tJRorcp0UJq6QkQkwni9XqaMG0fpgw8yZdw4vF7v/n3Lly/nV7/6FT169ODZZ5/lyiuvZMWKFSxfvpyJN93EpsREWjRqxJr4eHr313oPB6PKUUQkwixcsIBBmzYxIiYGX14eCxcsoHHjxvz5z3/m7bffplmzZtxxxx3cdNNNtG3bdv951Q3ol+opOYqI1FEHazp1wIfA8c6xuLycLW++yQOPPELz5s158MEH+d3vfkdiYmK11xw8eLCSYhDUrCoiUgdVzGnaIjWVGRMnHtB0OmLkSIo6duSvzvFuURHfr1vHn/70J9atW8cdd9xx0MQowVPlKCISZtVViBnp6VwIjGzXDrZu3T8HqnOOzMxM0nfuJGfPHoYPH86MGTNo104r/h1LqhxFRMLoYBVi7/79eQtYsHUrbwXeb9y4kV/84heMGjWKdu3a8dFHH7Fw4UIlxhBQ5SgiEkbVrZLRrWfPAzrPXN+vH99++y3nnnsuzjmmTZvGxIkTj3pdRTk43VkRkTBan5nJ0vXrKTLjZY+HsZVWsRg8eDBdu3Zl/PjxvPbaawwdOpS5c+fSpUuX8AXcQKhZVUQkTLxeL8uee46hZmSZMTQpCVdcvH//e++9R69evXjzzTeZPHky77//vhJjLVFyFBEJk4z0dIaZsdGM9s6xEv+zRecckydP5txzz6Vly5YsX76c22+/HY/HE+6QGww1q4qIhInFxrIsN5dTneOx8nIuGTuWnj17ctFFF7Fs2TJGjx5NamoqCQkJ4Q61wVFyFBEJk3WrVzMoIYE2CQlcFRtLeWwsgwcP5rvvvmPKlCncfPPN9XIh4Uig5CgiEgZer5dVixYR7fOR7PPxbnIy6dOmER0dzZQpU5g0aVK4Q2zQ9MxRRCQMFi5YwIUFBfyuWzf+26QJS7dsoVWrVnz88cf07ds33OE1eKocRURCqLrZbyqqxq/z8iA7m3+XlHDiiSeSlpZGcnIymZmZYY5aVDmKiITIwWa/yUhPZ1x8PK2Tkni9pIRuXbuyfPlykpOTwxyxVFByFBEJkf2z3+zbx898PjLS0wHo1a8f92dlMXPTJrq2bMnsOXNo0qRJmKOVytSsKiISIhYby+zNmykBXgDGxsbinOOtt97i861bGXLqqUx+9FFOO+20cIcqVSg5ioiEiCsuZmhSEnklJQyNiaF83z7+/Oc/89BDD3HdddcxY8YMoqLUgFcXKTmKiITI+sxMlmZnM9aMf3s8NFu8mDfeeIMbbriB6dOnKzHWYfqXEREJgarzpjZp1Ig33niDq6++WokxAqhyFBEJgYp5U783I760lLfz8xkyZAizZs1SYowA+hcSEQmBinlTXVkZ/ywtpWPHjixevJiYmJhwhyZBUHIUETnGvF4vH7z7Lj0SE/lXWRmtGzXitkmTNFwjgqhZVUTkGKoY+N9+1y6mbNtGnMfDz3r04JTTTw93aHIYVDmKiBxDCxcsoP+2bbyfk0MZcPpZZ3HLrFn7p46TyKDKUUTkKFXMn2qxsXyxcCEvbN3KxrIyTmnfnrvuu0+JMQIpOYqIHIWKZtQLgcdycui0bx8by8o4OTGRM6+8UokxQqlZVUTkKCxcsIBBOTmc1bgxKcCLOTkMat6c3t26MWLkyHCHJ0dIlaOIyBGqvPTU1txcni8upm3btlz1xz8y8LTTVDVGMCVHEZEj4PV6eXraNEY5R/e+fRm9YgWlZixevJjevXuHOzw5SkqOIiKHqeI5Y4+CAmZv3kyL3Fy2BCYVV2KsH/TMUUTkMFU8Z7yhVSu6tmrFv/LyuOCCC7j//vvDHZocI6ocRUQOoWKoRu/+/QH2P2fMDjxn7NKlC6+++ipmFuZI5VhRchQRqUHloRozUlNpd9ZZjIuPp0Pfvly+ciWlwMKFC4mPjw93qHIMKTmKiNRg/1CNzp2hqIgVwFtAky1b2FBYyI033qjnjPWQkqOIyEFUHqrhy8vj086duXXkSLJOPZUrrriCfv36MW3atHCHKSGg5CgichAZ6emMco7mbduy0Oej2/DhDBw4kLPOOovo6Ghee+01rc1YTyk5iogchMXGMnvzZq4E1gGnpaQwbdo0lixZQmpqKl26dAlzhBIqSo4iIgfhiosZmpREXkkJQ2Ni2Lp5M/c/9BAXX3wx11xzTbjDkxBSchQROYj1mZkszc5mrBn/FxXF1pdeIi4ujhkzZmjYRj2n5CgiUg2v18uy555jqBlZZjSNi2PJd98xa9Ys2rZtG+7wJMT0JFlEpBoLFyzg1H37yIyOJtGMd/Lz6devH9dee224Q5NaoMpRRKSKiiEc0T4fncrKeABwHg8vvfSSmlMbCFWOIiJVVAzhGNe2Ld80bsyOffu45557SElJCXdoUktUOYqIVFExhGOkcywrKqJ9+/b88Y9/DHdYUotUOYqIVLFu9WoGJSSwsFEj9jnH2DFjiI2NDXdYUotUOYqIVJKamsqSuXNptHMnS4qL6ZyYyKWjR4c7LKllqhxFRAK8Xi9z7r2XcXv2sKOsDDNjxJVXMnjw4HCHJrVMyVFEJCAjPZ1hZjxfUsLnZWV0b9aMK6++OtxhSRioWVVEJMBiY1myfTsZpaXEAhePH6+qsYFS5SgiEuCKi4lu1Ih85xiZnEy3444Ld0gSJkqOIiIBGd9+S9rOnZwUFcXWggJMPVQbLDWriojg74wzLzUVgKHR0bRMSsIVF4c5KgkXVY4iIsAzM2eyraCA7tHR9I2OZiXQu3//cIclYaLkKCINntfr5Y158/AAQ6OimNW0Kb+8+251xmnAQpoczex8M1ttZmvN7I5q9nc2s/fNbKWZrTKzX4QyHhGR6ix45RW2+XyMadGC8ubNOfvXv9bqGw1cyJKjmXmA6cAFwInA5WZ2YpXD7gLmOef6AWOAJ0MVj4hIdZxzzH/9dWLNOKWoiHU+H101wXiDF8oOOYOAtc65dQBm9jIwAvi60jEOSAy8bgZkhTAeEZEfefXVV1m/fj0Xt2yJr1EjhsbEqCOOYM650FzYbBRwvnPuusD7q4DBzrnfVzqmHbAYaAHEA+c451ZUc63rgesBkpOTT543b15IYq7vfD4fCQkJ4Q4j4ui+Hbm6fu9KSkoYN24cHo+HP910E8lRUewwo2WnTiQlJYUtrrp+3+qyM888c4VzbsDRXieUlWN1K4JWzcSXA7Odc1PM7FTgOTPr5ZwrP+Ak52YCMwFSUlLcsGHDQhFvvZeWlobu3eHTfTtydf3ezZgxg6ysLH6alMTmu+4i2wzXujVRd9zBqFGjwhZXXb9vDUEoO+RkAp0qve/Ij5tNrwXmATjnPgYaA+H7c01EGoyioiLuv/9+OrRvzy/MyIyOplNUlIZwCBDa5Lgc6G5mx5tZLP4ONwurHLMJOBvAzHriT47bQxiTiAgAM2fOJDMzk04eD2sKCmheWqohHLJfyJpVnXOlZvZ74B3AAzzjnPvKzO4DPnPOLQRuBWaZ2c34m1zHuVA9BBURCdi7dy8PPPAAKd27c0NUFM1jYljo82kIh+wX0unjnHNvA29X2XZ3pddfA0NCGYOISFXTp09n27ZtjBs3jtmPP86VwDrgNA3hkADNrSoiDcqePXt4+OGHOf/88+l23HHEJiWRV1KiIRxyAE0fJyINymOPPUZeXh5//etf/es35uSQmJvLkpwcrcIh+yk5ikiDsWvXLh599FFGjBjBgAEDWLd6NX2co8DjoY9zrFu9OtwhSh2h5CgiDcYTTzzBrl27uOeeewDYkpXFuvJyOnk8bPZ4fjQQWxouJUcRaRB8Ph9Tp07lwgsvpF+/fni9XvJWrGBPVBSvFRezOzmZESNHhjtMqSPUIUdEGoSnn36avLw87rrrLgAy0tMZ5RzNO3Rgoc9H19GjNb5R9lPlKCL1XmFhIY8++ijnnHMOp5xyCgAWG8vszZvZnp2tlTjkR1Q5iki9l5qayrZt23jllVf2b3PFxQzVMA45CCVHEanX9u3bx8MPP8zpp5/OGWecgdfrJSM9nfWZmXyUm8uVwL+BsRrGIZUoOYpIvTZ37lwyMzNJTU3l008/ZcbEiVwIzM3JYVhSEiQmcl5pqSpHOYCeOYpIvVVaWspDDz3EwIEDOffcc8lIT+dnPh8D9+1jmBkrgRaNGrEmPl4rccgBVDmKSL21YMEC1q1bx6OPPoqZ7e+EUwIsA06bNIn8jh0Z37+/eqrKAZQcRaRecs4xefJkevTowYgRI/B6vXzw7rsMrdSUmtyxI9dNmBDuUKUOUrOqiNRL7733Hunp6dx+++0sX76cGRMn0uPzz1mWm0tMSYmaUqVGqhxFpF6aPHkybdu2pUePHjw9bRrDfD6u7t4dgLS+fblh0iQ1pcpBKTmKSL3z+eefs3jxYi699FIeueoq+gGzc3MBWJOQoMQoh6TkKCL1zqOPPkqTJk3Y6vXym/x8BsTEQFISaf36KTFKUPTMUUTqlY0bN/Lyyy8z5NRTOcfj4f3ycrxFRawEJUYJmpKjiNQrU6dOxcw49+c/Z1luLsc5x2Pl5Zw0dqwSowRNzaoiUi94vV68H3zAzJkzufzyy2mRkMB5nTuTHB3NVYFhGyLBUnIUkYjn9XqZMm4cpVu3UlhYyLnnnkuPHj2YER9PD+CjRo04S8M25DAoOYpIxFu4YAEnb9zIP4qK6BoVxbcZGVx11VUwfToZ6emaAUcOm5KjiEQ8B7xaUkKWcwz0eHCB7YMHD1ZSlCOiDjkiEvG6pqTwbVkZTYHSqCgtXCxHTZWjiES8tatXU+AcwxMT6du8uZafkqOmylFEIt5/09LwAOeVlLAsNxfTwsVylFQ5ikhE27x5M5999hkD4uIoa9eO82JjVTnKUVPlKCIRy+v1MnrUKJxztHCOkuxsVoJW25CjpuQoIhHJ6/UyfcIEVn72Gb09Hm7r0YN1rVrRbfhw9VCVo6bkKCIRKSM9neisLPaVl9M0Korv8/Mpat2aESNHhjs0qQeUHEUkMsXEMC8nh45mFJrhPeMMxk+frqpRjgl1yBGRiLT0vfcocI6LkpLok5RE8pAhSoxyzKhyFJGI4/V6+e+bbxIDnFRQoE44csypchSRiPPCnDlk79nDJUlJZMbFqROOHHOqHEUkoni9Xha9+CLlQNfdu8mNj1cnHDnmlBxFJKL8a9488nbtYqDHwwbnaHLyyaoa5ZhTchSRiPKx18se4OzoaIpjYmjfvn24Q5J6SMlRRCKG1+slIz2dxmZ8X1bG7uRkNalKSKhDjohEjP9btIj8wkKuadkSi46m2+jRalKVkFDlKCIRY+lHH2FAn8JC1vl8WrdRQkaVo4hEBJ/Px8cff0zv+HiKEhMZGhOj1TckZFQ5ikhEePHFFykqKiKqpITE3FyW5ORo3UYJGSVHEanznHM89dRTtE5O5nSgwOOhj3OsW7063KFJPaXkKCJ13qeffsrnn3/OCV26sN45Onk8bPZ4cOEOTOqtQyZHM/vMzCaaWYvaCEhEpKqnnnqKJk2a0HznTvZERfFacbGGcUhIBVM5jgHaA8vN7GUz+7mZWYjjEhEBYMeOHbzyyiukdO/OiKIi7kpJoVXbtgzWMA4JoUMmR+fcWufcnUAP4EXgGWCTmd1rZi1DHaCINGxz5syhqKiIprt28U5eHl+sXav5VCXkgnrmaGYnAVOAR4AFwChgN/Be6EITkYbOOceMGTPo0K4dl5eV8btu3VjXqpVW4ZCQO+Q4RzNbAewEUoE7nHP7Aru8ZjYklMGJSMM2Y8YMVq9eTc9WrXhn1y6G5OWR27kzV6tqlBALZhKAXznn1lXeYGbHO+fWO+cuDVFcItLAeb1eHr3zTmLNGLt3LwO7d+df+fmqGqVWBNOsOj/IbSIix8wny5axMT+fs5o25ZPycr7Pz6eodWs9a5RacdDK0cx+AvwUaGZmlSvERKBxqAMTkYbtkxUrKHWOgcXFvAV4zziDG268UVWj1IqamlVTgIuA5sDwStv3AL8NZVAi0rA553j33XdpGx1Nq/btGRUbS/KQIUqMUmsOmhydc28Ab5jZqc65j2sxJhFp4O68805yc3M5KTqakuxsVnbuzK39+4c7LGlAampW/aNzbjJwhZldXnW/c+6mkEYmIg2S1+vlmX/8g0bAxNhYvAkJ6oQjta6mZtVvAt8/q41AREQAvB9+SN7evRzn8eDKyshp3Jjr1QlHallNzaqLAt/n1F44ItLQfZqeTqlznB4VxWPl5VwydqyqRql1NTWrLoKDT3rvnLs4JBGJSIP27n/+Q+voaPp27EhKo0Ykd+wY7pCkAaqpWfXRWotCRAT/PKrbcnI4KSaG0m3b1BFHwqamZtUltRmIiMjsZ54h2ow/tm/Pe4WF6ogjYXPQGXLMbF7ge4aZrar0lWFmq4K5uJmdb2arzWytmd1xkGNGm9nXZvaVmb14ZD+GiES6wsJCvMuX08KM3VlZrN29m64pKeEOSxqomppV/xD4ftGRXNjMPMB04FwgE/96kAudc19XOqY78CdgiHMu38xaH8lniUjke+211ygsLOTnMTEUeDz0cY51q1eHOyxpoA5aOTrntga+bwT2AX2Ak4B9gW2HMghY65xb55wrBl4GRlQ55rfAdOdcfuCzcg7/RxCR+iA1NZWEhASKnaOTx8Nmj+fgPQJFQiyYJauuA+7Gv3ajAY+b2X3OuWcOcWoHYHOl95lA1YcHPQKf8SHgAf7XOffvamK4HrgeIDk5mbS0tEOFLdXw+Xy6d0dA9+3IBXvvtmzZwvvvv8+Iiy7igjPPpAC4ICaGDiec0CDvvf6bC79glqy6HejnnMsDMLNWwEfAoZKjVbOt6h+C0UB3YBjQEVhmZr2cczsPOMm5mcBMgJSUFDds2LAgwpaq0tLS0L07fLpvRy7Ye3fnnXdiZvzs669p+/XXLPT5aHP11Vx/442hD7IO0n9z4RfMklWZ+Ccbr7CHAyvCms7rVOl9RyCrmmPecM6VOOfWA8pL8eAAACAASURBVKvxJ0sRaSBKS0uZPXs2vXr1YmF2Ntuzs1nn86kzjoRVTZMA3BJ4uQXwmtkb+Cu/EcCnQVx7OdDdzI4PXGMMcEWVY14HLgdmm1kS/mbWdYhIg/HOO++QlZXFKQMH0m3bNgoSEjgvNhZXXBzu0KQBq6lybBr4+h5/EqtoEn0D2HqoCzvnSoHfA+/gn6d1nnPuKzO7z8wqZtd5B8gzs6+B94HbK5pvRaRhSE1NpUWLFuxYuZKv8vPZt2ULK4HeGvwvYVTTJAD3Hu3FnXNvA29X2XZ3pdcOuCXwJSINTE5ODosWLaJ/375c+PXX9IyNZW5xMU1OPlmD/yWsgumtmgz8Efgp0Lhiu3PurBDGJSINwAsvvEBpaSkn9e7N0owMOphRGBVF1/btwx2aNHDBdMh5AfgWOB64F9iA/3miiMgRc87xzDPPMHDgQE4ZMoQtwOKyMraAOuNI2AUzlKOVcy7VzP4QmG91iZlp3lUROSorVqzgyy+/5KmnnmLd6tWc06wZ7RMS6KHOOFIHBFM5lgS+bzWzC82sH/5hGSIiR+yZZ56hcePGdOvWjVWLFvGdz0dJdrY640idEEzleL+ZNQNuBR4HEoGbQxqViNRrhYWFvPTSS1x66aW8v3gxFxYU0LVbN/6Vn6+VOKROOGRydM69GXi5CzgztOGISEPw+uuvs3PnTk455RRenzKFL7KzOSMnh9zjj+fqkSPDHZ7IoZtVzewEM1tkZrlmlmNmb5jZCbURnIjUT88++yydO3dm66ZNXLB9OxNiY1lRXq4hHFJnBPPM8UVgHtAWaA+8CrwUyqBEpP7atGkT7777Lif16kXW1q18CBR5PBTHxNBeQzikjggmOZpz7jnnXGng63l+PIG4iEhQ/va3v+Gc44JNm8hbsYLdycn8t1kzSjt3ZoSaVKWOqGlu1ZaBl++b2R3412N0wGXAW7UQm4jUM+Xl5cx75RV6NGrE5c2a0aaoiBXDh3PC8cdzdf/+alKVOqOmDjkr8CfDiqWnxlfa54C/hiooEamfZs6cyc5duzghNpbUlSv5tHNnbh05UklR6pya5lY9vjYDEZH6y+v1kpGezlNPPklcVBRT+vdn3pYtGrYhdVYwc6vGABOAMwKb0oAZzrmSg54kIhLg9XqZMm4cp+3Zw1dZWXROTCSvpISi1q31jFHqrGAmAXgKiAGeDLy/KrDtulAFJSL1x8IFCxi0aRPmHGXOkTJoEPmXXMJ4PWOUOiyY5DjQOden0vv3zOyLUAUkIvWLA5aWlbG2uJh4M/r26cN1EyaEOyyRGgUzlKPMzLpWvAlMAFAWupBEpD7pmpLCOuf4xjkSo6Lo9pOfhDskkUMKpnK8Hf9wjnX4e64eB1wT0qhEpN5wxcW0atwYKy7mV23aaMUNiQg1JkcziwIKge5ACv7k+K1zbl8txCYi9cC6zZvx7t7NT8xYtXMnJ8XGhjskkUOqMTk658rNbIpz7lRgVS3FJCL1REFBAW/+85/sA3pGR/PTpCRVjhIRgnnmuNjMRpqZHfpQEZEfFO7di+3dSwwwzExrNUrECOaZ4y1APFBqZkX4m1adcy4xpJGJSMQrLCriy4IC+kRF8ZRzXDJ2rIZvSEQIZj3HprURiIjUL16vl7T336cc+HmXLiTGxpLcsWO4wxIJSk0Tj7cG/gx0w/+88SHn3O7aCkxEIpfX62XGxIl8vGsXcWakREezND6es9SkKhGipmeOc4EC4HGgKfCPWolIRCJeRno6PXfs4Nu1azklMZEl/fszfvp0NalKxKipWbWtc+7OwOt3zCy9NgISkchnsbE8sWkTZsbeoiJOP+ccJUaJKDVVjmZmLcysZWBtR0+V9yIi1fr+22/JLyujd/fuDALWrV4d7pBEDktNlWMz/Gs6Vh7CUVE9OuCEUAUlIpFt4+bN7AHOGjiQ7zdupFe4AxI5TDWt59ilFuMQkXpk/aZNRAE9e/dm6fz5jEhJCXdIIoclmEkARESCtmfPHlasWEHfhAQS4+K4UPOpSgRSchSRY2r+/PkUFxfjSkqILitjWW4upvlUJcIEM0OOiEjQZs+eTcsWLTgnOprymBjO69xZlaNEnKAqRzM73cyuCbxONrPjQxuWiESi77//nqVLl9LcOb7euRNXXKz5VCUiHbJyNLN7gAH4l6x6FogBngeGhDY0EYk0c+fOBeDyoiJOjY1lG9Dk5JM1xlEiTjCV4yXAxfhny8E5l4V/xhwRkf3Ky8uZM2cO7du14/OyMnyAM6N9+/bhDk3ksAWTHIudcw7/2EbMLD60IYlIJFqyZAkbN26kjRl7oqJ4rbiYsuhoRowcGe7QRA5bMB1y5pnZDKC5mf0W+A0wK7RhiUikmT17NnGNGzMhPp7WjRuz0OcjvmVLNalKRApmyapHzexcYDf+5453O+f+E/LIRCRi+Hw+FixYwMBBg3jxs8+4ElgHXNCoUbhDEzkiwXTIuRl4VQlRRA5m/vz5FBQU0KVjR9p++y0FCQmcFxsLzoU7NJEjEswzx0T8q3IsM7OJZtYm1EGJSGSZPXs2bdq0Ye277/LNrl2UZGezEohr0iTcoYkckUMmR+fcvc65nwITgfbAEjN7N+SRiUhEWLduHUuWLCGuqIhrfD5+FR3NNwkJdBs+nPh49d+TyHQ408flANlAHtA6NOGISKSZM2cOZsaI+HjeLy9nb2kpOY0bq5eqRLRDJkczm2BmacB/gSTgt865k0IdmIjUfRVjG3v27MkXO3dynHM8Vl7OSWPHqpeqRLRghnIcB0xyzn0e6mBEJLKkpaWxceNGRlx0ESm5ubRJSOCq2FiSO3YMd2giR+WglaOZJQZeTgY2mVnLyl+1E56I1GXPPvssCQkJ7PziC77Kz2ffli2aS1XqhZoqxxeBi4AV+GfHsUr7HHBCCOMSkTpu9+7dLFiwgJ/06MFFa9bQMzaWucXFmktV6oWDJkfn3EWB71qBQ0R+ZN68eRQWFtK7Vy+WfvstHcwojIqiq+ZSlXogmA45/w1mm4g0LM8++yw9e/Zk6LBhbAEWl5WxBeiakhLu0ESO2kErRzNrDDQBksysBT80qybiH+8oIg3U6tWr+eijj5g4cSIf/ve/XNimDZ0TE+lRWqqFjaVeqOmZ43hgEv5EuIIfkuNuYHqI4xKROmzOnDl4PB5yly2jT3Exi3NzGRcXx5qEBM5SZxypB2p65vgY8JiZ3eice7wWYxKROqysrIy5c+dy/HHHMSw/n8s6dwYgrW9fbpg0SZ1xpF4IZlWOx82sF3Ai0LjS9rmhDExE6qbHH3+cLVu20CspiXf27sWXl8fKzp25VYlR6pFgVuW4BxiGPzm+DVwAfAAoOYo0MF6vl8fuu49GUVFcWVDAyd2786/8fLoNH67EKPVKMHOrjgLOBrKdc9cAfQAt0ibSAH2ybBmZO3dybkICH5WX831+PkWtW2seVal3gkmOhc65cqA0MGtODpoAQKRB+nTlSkqdo39xMVsA7xlnMH76dFWNUu8EM7fqZ2bWHJiFv9eqD/g0pFGJSJ30wbJltI2JoUmLFlwYE0PnIUOUGKVeCqZDzu8CL582s38Dic65VaENS0TqmpUrV7Jp82a6x8TQatcu/g2MjY0Nd1giIVHTJAAHHaxkZv2dc+mhCUlE6qIHH3yQqKgoRrRuDS1acJ4G/Es9VlPlOKWGfQ446xjHIiJ11LJly3hjwQJ6xMTwaV4evRISNOBf6rWaJgE4szYDEZG6a+qjj1JcXs7dHg8fl5Uxv1Mn7rr/fj1vlHormHGOv65uuyYBEGk4vsjIoDFQXlbGBqB3v35KjFKvBdNbdWCl143xj3lMR5MAiDQImzdvZv2GDbTxeHi3vJwtUVGM0MobUs8F01v1xsrvzawZ8FwwFzez84HHAA/wT+fcQwc5bhTwKjDQOfdZMNcWkdoxd+5cnHNc3Lw5JzZrRo/YWHXEkXovmEkAqtoLdD/UQWbmwb96xwX4p5673MxOrOa4psBNgPcIYhGRECovL+fJJ5+kVZMm5BQWUpKdzUqgtzriSD0XzDPHRfh7p4I/mZ4IzAvi2oOAtc65dYHrvAyMAL6uctxfgcnAbUHGLCK1JC0tjaysLMa3a8cFcXEs9Pk0j6o0CME8c3y00utSYKNzLjOI8zoAmyu9zwQO+D/KzPoBnZxzb5rZQZOjmV0PXA+QnJxMWlpaEB8vVfl8Pt27I9CQ79vf/vY34uPj6fWHP1AcE8MgoFXnzkHfj4Z8746G7lv4BfPMcQlAYF7V6MDrls65HYc41arZ5vbvNIsCpgLjgohhJjATICUlxQ0bNuxQp0g10tLS0L07fA31vu3YsYMPP/yQXieeyHf33EORGZuco+lNN/G3yZODukZDvXdHS/ct/A75zNHMrjezbcAq4DP886sG02kmE+hU6X1HIKvS+6ZALyDNzDYApwALzWxAcKGLSCi98MIL7Nu3j/bt2rGuvJxOHg+bPZ4f/sIVqceCaVa9Hfipcy73MK+9HOhuZscDW4AxwBUVO51zu4Ckivdmlgbcpt6qIuHnnGPWrFn07NkTt24de6KieK24mN3t22t5KmkQgumt+j3+HqqHxTlXCvweeAf4BpjnnPvKzO4zs4sP93oiUns+++wzMjIy6N+nD6Oc45YOHUho0YLBo0erM440CMFUjn8CPjIzL7CvYqNz7qZDneicext4u8q2uw9y7LAgYhGRWpCamkpcXByDTzuN2QsXciWwDjhNg/+lgQgmOc4A3gMygPLQhiMi4VZQUMCLL77I6NGjiYuOZmhSEnklJQyNidHgf2kwgkmOpc65W0IeiYjUCa+++ip79uyhTatWrM/M5KPcXK4Erd8oDUowyfH9wDjDRRzYrHqooRwiEoH+/ve/06xRIwampfGP7dsZlpQEiYlav1EalGCSY0UP0z9V2uaAE459OCISTl9++SUZGRmMSU5mUHExw8xYCfRp1IiPGjXS+o3SYAQzCcDxtRGIiITfzJkziY6OZpPPx+KCApYBp02aRH7Hjozv3189VaXB0HqOIgLA3r17ee655+jfty9n5+Ts74TTuWNHrpswIdzhidSqYMY5Dqz0NRT4X0DjFEXqmfnz57Nz5066p6TwfnY2zXfsYFluLqZOONIAhXQ9RxGJHDNnzqRTp05sXLKEoWZkmTE0KUmdcKRBCtl6jiISOb766is+/PBDup9wAqcVF5MZHU2nqCit3SgNVijXcxSRCDFr1ixiYmKI3rKF73w+OpWVMat5c66/+251wpEGKZTrOYpIBCgsLGTu3Ll079qVSwoK6NqtG//Kz+fsMWO49tprwx2eSFgcNDmaWTegTcV6jpW2DzWzRs6570MenYiE3Pz588nPzyclIYF38vIYkpdHbufOXK3VN6QBq+mZ4zRgTzXbCwP7RKQeeOqpp2jdujXj4+IY17Yt3yQk0G34cDWnSoNWU3Ls4pxbVXVjYL3FLiGLSERqzRdffMHHH3/MsGHDeHbTJrI2b2bt7t101eob0sDV9MyxcQ374o51ICJS+5566ikaN25Mh9ataeMcBR4PfZxj3erV4Q5NJKxqqhyXm9lvq240s2uBFaELSURqw+7du3n++ecZM2YMjePi2AB0ADbwQ/d0kYaqpspxEvAvM7uSH5LhACAWuCTUgYlIaD3//PMUFBQwYcIEMjIyeAdYXFbGlqgoRqhZVRq4gyZH59w24DQzOxPoFdj8lnPuvVqJTERCxjnHk08+Sf/+/Rk4cCBvzJ/POc2a0T4hgR6xsZoVRxq8YKaPex94vxZiEZFa8sEHH/DVV1/xz3/+k08//ZRVixYR7fOR7POxsnNnbtWsONLAHcn0cSIS4Z566imaNWvGmDFjyEhPZ5RzGsYhUomSo0gDk5OTw/z587n66quJj4/HYmN5ZuNGDeMQqUTJUaSBSU1NpaSkhBtuuAGAdatX00fDOEQOoOQo0oCUlpby5JNPcvbZZ9OzZ0/AP2xjAxrGIVKZkqNIA/L666+TmZnJjTf+sExr15QUthAYxhF4L9LQBbMqh4jUE48//jhdunThoosu2r/NFRdzYZs2xJeU0CkmRsM4RFDlKNJgrFq1iqVLlzJx4kQ8Hs/+7RYby7LcXFrt2sWy3FwsNjaMUYrUDaocRRqIxx9/nLi4OH7zm98csN0VFzM0KYm8khKGqnIUAVQ5ijQIeXl5vPDCC4wdO5aWLVsesG99ZibvZ2fTfMcOVY4iAUqOIg1AamoqhYWFB3TEAfB6vSx77jmGmpFlxtCkJFWOIig5itR7ZWVlPPnkkwwbNozevXsfsG/hggWcum8fmdHRdIqKYiXQW1PHiSg5itR3ixYtYuPGjdVWjasWLeI7n4/mpaXMatqUX959t6aOE0EdckTqvalTp3Lcccdx8cUXH7B94YIFXFhQQNdu3fhXfj5njxnDtddeG6YoReoWVY4i9dhnn33G0qVL+cMf/kB09A9/C1dUje/k5fHF2rXkxsczYuTIMEYqUreochSpx6ZOnUrTpk1/VBFWrMTRvG1bFvp8WolDpApVjiL1VGZmJvPmzeO3v/0tiYmJB+zTShwiNVNyFKmnnnjiCcrLy3/UEQe0EofIoSg5itRDPp+PGTNmMHLkSLp06XLAPq/Xyxeff84G4PjYWDZ7PFqJQ6QKPXMUqYdmz57Nzp07ueWWWw7Y7vV6mTFxIkMKCpgPvN24MaUtW6ozjkgVSo4i9UxZWRnTpk3j1FNP5ZRTTjlgX0Z6Oj/z+TgzJobCNm3YfNpp3DVpkjrjiFSh5ChSzyxcuJDvv/+ehx566Ef71mdmsnT9eorMWObxMPacc5QYRaqhZ44i9YhzjocffpgTTjiBX/7ylwfs0zyqIsFTchSpR5YuXYrX6+W22247YNA/aB5VkcOh5ChSjzz88MO0bt2acePGHbBd86iKHB4lR5F64osvvuD//u//+MMf/kBcXNwB+ypmxBnXti0FzZpx9q9/rXlURWqg5ChST0yePJmEhAQmTJjwo30WG8vszZvZnp3NOp9PM+KIHIJ6q4rUA+vXr+eVV15h0qRJtGjR4kf7XXExQ5OSyCspYWhMjDriiByCKkeRemDKlClERUVx8803V7t/fWYm72dn03zHDpbl5mKxsbUcoUhkUXIUiXDbt2/nmWee4aqrrqJDhw4/2q8hHCKHT8lRJMJNnTqVoqIibr/99mr3Z6SnM8yMjWa0d05DOESCoOQoEsF27NjB448/zujRo/nJT35S7TEWG8uy3FyOc47Hyss5aexYDeEQOQR1yBGJYNOmTcPn83HXXXcd9JiKzjjxJSWMiomhc8eOtRihSGRS5SgSoXbu3Mljjz3GyJEj6dWr10GPq6gcW+3apc44IkFS5SgSof7xj3+we/fuGqtGr9fLB+++y9CkJEhM5LzSUnXGEQmCKkeRCLR7926mTp3KiBEj6Nu3b7XHpKamcv/o0XT66COW5eYSU1LCmvh4dcYRCYIqR5EI9MQTT7Bz507+8pe/VLvf6/Uy5957+XVeHgNiYiApibR+/bhBazeKBEXJUSTC7NmzhylTpnDhhRdy8sknV3tMxQocS82gpISVoEWNRQ6DmlVFIszjjz/Ojh07aqwatQKHyNFRchSJIPn5+UyePJnhw4dXm+y8Xi9PT5vGKOe4o18/Slu31gocIkdAzaoiEeSRRx5h165d3H///T/a5/V6mTFxIj0KCpi9eTPjgKLWrRkxcmStxykS6ZQcRSJEdnY2jz32GJdffjknnXTS/u1er5eM9HTWr1/PhcDIbt0ASOvbVx1wRI5QSJOjmZ0PPAZ4gH865x6qsv8W4DqgFNgO/MY5tzGUMYlEqgcffJB9+/Zx77337t+WmprK6/fdx6i4OFaZsQZg61bWxMcrMYochZAlRzPzANOBc4FMYLmZLXTOfV3psJXAAOfcXjObAEwGLgtVTCKRauPGjTz99NNcc801dO/eHThwuEafmBho25bVF19M/vHHM75/fyVGkaMQyspxELDWObcOwMxeBkYA+5Ojc+79Ssd/AowNYTwiEeu+++4D4O677wZ+6HgzzIyPYmKgpIQ3Cgu5a+RIJUWRY8Ccc6G5sNko4Hzn3HWB91cBg51zvz/I8U8A2c65H/U0MLPrgesBkpOTT543b15IYq7vfD4fCQkJ4Q4j4oT7vm3atIlrrrmGSy+9lIkTJ1JQUEDupk00Ki9nd3ExCdHR7AFatWtHUlJS2OKsTrjvXaTSfTtyZ5555grn3ICjvU4oK0erZlu1mdjMxgIDgJ9Vt985NxOYCZCSkuKGDRt2jEJsWNLS0tC9O3zhvm8jRowgPj6e6dOns379eubMmMGwlSsZ0707D65dy8d1uONNuO9dpNJ9C79QJsdMoFOl9x2BrKoHmdk5wJ3Az5xz+0IYj0jESUtLY+HChTzwwAOsX7/+gKEaAGsSEupsYhSJZKGcBGA50N3MjjezWGAMsLDyAWbWD5gBXOycywlhLCIRp7y8nFtvvZXOnTszadIkMtLT+ZnPx5XR0QwNzJU6fvp0JUaREAhZ5eicKzWz3wPv4B/K8Yxz7iszuw/4zDm3EHgESABeNTOATc65i0MVk0gkef7550lPT+f5558nLi6O9ZmZLF2/niIzlnk8jD3nHCVGkRAJ6ThH59zbwNtVtt1d6fU5ofx8kUi1d+9e7rzzTgYMGMDll1+O1+tl2XPPMdSMLDOGJiVpXUaRENLcqiJ10NSpU8nMzGTKlCksX758/7CN7MaN6RQVxUrQuowiIaTp40TqmOzsbB566CF++ctf0qhRo/2dcBbn5jI0KYm5znG1VtkQCSklR5E65n/+538oKipiQN++LFyw4ID5Ur/r25dH1DtVJOSUHEXqkA8++IC5c+fyk1ataJWayntRUayJi9N8qSK1TMlRpI4oLS1l4sSJJDZtylV79zKitBRfSQlfXnIJ+UOGaL5UkVqk5ChSRzz55JOsWrWKYUOH8tEnn9A1KooPgZ7t23PdhAnhDk+kQVFvVZE6YNu2bfzlL39h8ODBxOfksCcqiteKi9mdnKzFikXCQMlRpA74zW9+Q0FBAX179eKahAQmDxhAq7ZtGTx6tJpSRcJAyVEkzJ5++mnefvttLk5OZsuHH/JKQQGZRUUUtW6tqlEkTPTMUSRMvF4vKz75hL/cey9tYmN5rm9f/r19OyvOOksLFouEmZKjSBh4vV6mjBvHri1b2LFnD2d06sS/t2/nLWC8FiwWCTslR5EwWLhgAe03bOC1oiIGeTwcf8YZGq4hUocoOYqEwebMTN4oKiLRjKTYWNpruIZInaIOOSK1zOv18ul//sNuoGdUlDreiNRBqhxFQszr9ZKRnk7vQJPpv998k+/z8jirVStSmjShuYZriNQ5So4iIeT1epkxcSIXAjNSUymeOpUXXnmFaI+HKzp04MOYGFWNInWQkqNICGWkp/Mzn4+BMTH4Skr42333sWbNGh555BFcfLw64IjUUUqOIiHi9XpZ9uGHbNi4kZKoKJ4oK2PVd98xfvx4brvttnCHJyI1UHIUCYGK5tT4nBz6OMf25GQ2bNlC6+RkpkyZEu7wROQQlBxFQqCiObV5TAz/dI6Pt29nd2kp0x5+mPj4+HCHJyKHoKEcIiFgsbHM3ryZ7dnZrCov57O9e7l+/HjGjRsX7tBEJAhKjiLHmNfr5YN332VoUhIb27Ujq6yME3v2ZPr06eEOTUSCpGZVkWOo4lljj4IC3tq+ne9KS4mNjmbK3/+Ox+MJd3giEiQlR5FjoGKg//r167kQuKRrV+ZmZpJbVMTMmTM5//zzwx2iiBwGJUeRo1SxwsZZBQV8EhXFmrg4/vX553zr83HzzTdz3XXXhTtEETlMeuYocpQWLljAoE2bGLF7Nxds305+u3a8mJXF2WefrWEbIhFKlaPIUXLA0rIyOpSX82ZZGR9/+CF9+vTh9ddfx8zCHZ6IHAFVjiJHqWtKCluAN0pL+bikhLi4OBYtWkRCQkK4QxORI6TkKHKUXHEx5yUn87EZ5WbcOGECHTt2DHdYInIU1KwqcpRcdDTPbN1KTlkZP42NpUu3buEOSUSOkpKjyGGqvD7joEGDeOXFF8kpK+PSNm0Y0KIFrrg43CGKyFFSchQ5DKmpqbx+332MiotjRkICs/r25b9pafy0dWsu79CBt804q3//cIcpIkdJyVEkSF6vlzn33vv/7d19dFTlncDx78+EBDJEDKK8mg0UdSFaJaLBl1WgikqQN9HKrt1Ystbt4vGI7Vnpii9Ie7S4WtdtaKQgyB6wxSKSVtkqGtSjJjTEF1CkpERX80ICRCETSDLMb/+4d+IwJGEySZhJ8vucc8/cO/Pcm+f5ccMvz713nod/PnCAi/r04c8JCbywfTvz5s3jzjvvZOcHH9j8jMb0EJYcjQlT/oYNXN7QwNsivHvkCC988w2TJk1i+fLlxMXFMWHChGhX0RjTSSw5GhOGoqIiPv7jH4mvq+NwYyNv+Hykp6ezefNmGzPVmB7IkqMxrQg8eDM0NZXSXbu4w+OhcMgQ/nPvXkampbFt2zYSExOjXU1jTBew5GiMK/gpVIBn588nC9ifk4MkJPB4VRXF5eWcM2AAz69ZQ1JSUnQrbIzpMpYcjeHbqaaygGdXrmTo5MlkATcPHcpGVVZt3EhxeTnjMzJ4+plnuPLKK6NdZWNMF7LkaAywo6SEa+rquLRPH+qamtgNvAL4KipYtnYtb7/3HtnZ2axcudLuMRrTC1hyNAaQhARWf/klTcBa4Pbzz2fylCncc/fd7CopYdGiRSxevJjTTrMRF43pDSw5ml4p+P5iZmamMz5qairExzPF56O6ooInn3yS0r17uf/++1myZEm0q2yMOYUsOZpeJ/T+Irm5XJiRwbMeD+cBLzY0sO2ppxARXnvttWhX1xgTBXaNyPQaRUVFrPjNb8jfsKH5YZssgviyjwAADchJREFUnPuNmZmZ/Mszz/D84MG8XlrK8OHDKSwsZOLEiVGutTEmGqznaHqM0EuloZ8FeoubvF72AFRW8gpwV0YG1dXVPPzww2zZsoXs7GyWLVtmX9Uwphez5Gh6hJYulQYnyB0lJc29RSor2T55MrUjR3JXRgY1NTXMmjWL2tpaVqxYwbx58xCRqLXFGBN9dlnVdHtFRUXkPf0019TVHXepNNiFGRm8Amxwe4szbr6ZOXPnkpeXx0033cSgQYMoLCwkJyfHEqMxxnqOpnsL9BjP83pZ/eWXALzVvz93hUwblZmZCbm57Cgp4a6MDA4ePMgFF1xAVVUVixYt4sEHHyQhISEaTTDGxCBLjqbbCb632Hy5dPRoALZefDH/eu+9LU4blZmZydlnn82CBQvYtGkTY8eO5eWXX2b8+PGnuAXGmFhnl1VNtxLoKaasXMmz8+cjCQnNl0v3eDytJsb6+noeeughxowZw5YtW3jssccoKSmxxGiMaZH1HE23kr9hA5dVVzM5NRWOHqW2sZG7gi6XhibGxsZGVq1axZIlSygvL2fu3LksXbqUESNGRKkFxpjuwJKj6TaKioooXL+ej6qqOFRdTfHIkfzETYihSdHn87F27VoWL15MWVkZl19+OevWrePqq6+OUu2NMd2JXVY13Ub+hg3cWFPDjxMS2O73k3TJJSckRa/XS25uLmPGjOGOO+4gJSWFV199lXfffdcSozEmbNZzNDEv8ABOeUUFnwIj4+Jo7NOH7wwb1lymvLyc3Nxc8vLyqK2tJTMzkyeeeIIZM2bYVzOMMe1mydHEtOAv9x/weqk/6yze8PvxeTxMnT6dl156ieeee47NmzcDMGvWLO677z6uuOKK6FbcGNOtWXI0MSnQWywrKztuZJviadPw+/0kffUVs2fPpqamhmHDhrFw4UJycnIYNWpUtKtujOkBLDmamBM6Dupnfj8fHzrEi19/Tc3q1ezfv5/ExESysrLIyclhypQpxMfbqWyM6Tz2P4qJirYGCf94+3bGHT3KPhGqKiv5qL6el5qaSExMZOrUqdxyyy1kZWVx+umnR6n2xpiezpKj6XRtJb7A58GDhDc89RR9+/bl/fff55133qGgoICDBw8C0D8hgenTp5Odnc2kSZPo37//KW6NMaY3suRoOlVbs2OoKvv27WP9unVIdTWvAgX797Nm0iSO+f0ApKWlkZWVxfDhw+nfty/fu/56JkyYEL0GGWN6JUuOJmwn6xECfFRczITGRoYmJzOgspLHfv5zBgwcSGlpKbt37+bAgQPNZZPj4kj2eLj91luZPn06l112mY1cY4yJCV2aHEXkBuC/gDhghao+HvJ5IrAGuAQ4AHxfVT/vyjqZ9lFVvF4vBQUFrFy4kHE+H+t8Pi6dMwePx0NVVRWVlZVUVVVRXl5ORUUFx44d+/YAZWWMGDGC0aNHM3v2bNLT00lPT6ehoYGKL77guy18kd8YY6Kty5KjiMQBucB1wFfAX0QkX1U/DSqWA9Sq6mgRuQ34JfD9rqpTtKnqcYvf729x2+/34/f7OXbs2AmvgcXn8zW/Bi9NTU00NjbS1NREU1MTDQ0NNDY2smvXLqr37WPjxo0MHjyYI0eOcPToUbxeL16vl/r6erxeL4cOHeLw4cPNy9dff43P52tuwyb3tWDpUgDOPPNMhg4dypAhQ5g4cSKpqan4fD6OHD7MhKuuYubMmfTr1y8K0TbGmMh1Zc/xMqBUVfcCiMjvgBlAcHKcATzirv8B+LWIiKpqawfds2dPxA9ltHHYFj8Pd7ul18ASzs+NBhGhX79+eDwekpKSSEpKwuPxkJycTFpaGsnJySQnJ5OSksIZZ5xBbW0thevXMzE+ng/i47lz6VKuu+46mwPRGNMjdWVyHA58GbT9FRB6/ay5jKr6ROQb4Exgf3AhEfkR8CN3s8Hr9e7skhr3UAJnDYFBx0DiQKtgv6rW1NfXU19f355Ded6CJIX6TdOmebuqvjFoECHnpAmbxS4yFrfInd8ZB+nK5NjSgJahXahwyqCqy4HlACJSrKo2CV8ELHaRsbhFzmIXGYtb5ESkuDOO05WzcnwFnBO0PQKoaK2MiMQDA4CDXVgnY4wx5qS6Mjn+BThXREaKSAJwG5AfUiYfyHbX5wBvtnW/0RhjjDkVuuyyqnsP8W7gzzhf5XhOVT8RkUeBYlXNB1YC/yMipTg9xtvCOPTyrqpzL2Cxi4zFLXIWu8hY3CLXKbET66gZY4wxx+vKy6rGGGNMt2TJ0RhjjAkRU8lRRG4Qkd0iUioiC1v4PFFEfu9+XiQiaUGf/cx9f7eIXH8q6x1tkcZNRNJE5IiIfOgueae67tEWRuyuFpESEfGJyJyQz7JFZI+7ZIfu25N1MG7Hgs650If0erwwYnefiHwqIh+LyBsi8ndBn9k5F1nc2n/OhQ5pFq0F56GdvwGjgATgI2BsSJl/A/Lc9duA37vrY93yicBI9zhx0W5TN4hbGrAz2m2I8dilAd/FGQN4TtD7A4G97muKu54S7TbFetzcz+qi3YYYj90kIMld/3HQ76udcxHEzd1u9zkXSz3H5uHmVLURCAw3F2wG8Ly7/gfgeyIi7vu/U9UGVS0DSt3j9QYdiVtvd9LYqernqvox4A/Z93rgdVU9qKq1wOvADaei0jGgI3Hr7cKJXYGqBoauKsT5jjjYORdp3CISS8mxpeHmhrdWRlV9QGC4uXD27ak6EjeAkSLygYi8JSL/0NWVjTEdOW/snPtWe9veV0SKRaRQRGZ2btViXntjlwNsjnDfnqQjcYMIzrlYms+xI8PNhTUMXQ/VkbhVAqmqekBELgFeFpF0VT3U2ZWMUR05b+ycO1572p6qqhUiMgp4U0R2qOrfOqlusS7s2InI7cB44Jr27tsDdSRuEME5F0s9x44MNxfOvj1VxHFzL0MfAFDV7TjX9M/r8hrHjo6cN3bOfatdbVfVCvd1L7AVGNeZlYtxYcVORK4FHgCmq2pDe/btoToSt8jOuWjfaA26YRqPc4N5JN/ecE0PKTOf4x8sWe+up3P8Azl76T0P5HQkbmcF4oRzo7scGBjtNsVS7ILKrubEB3LKcB6MSHHXe0XsOhi3FCDRXR8E7CHkwYqevIT5+zoO5w/Vc0Pet3MusrhFdM5FvdEhjZgK/NVt4APue4/i/BUA0Bd4EeeBm23AqKB9H3D32w3cGO22dIe4ATcDn7gnWglwU7TbEoOxuxTnr1YvcAD4JGjfeW5MS4EfRrst3SFuwBXADvec2wHkRLstMRi7LcA+4EN3ybdzLvK4RXrO2fBxxhhjTIhYuudojDHGxARLjsYYY0wIS47GGGNMCEuOxhhjTAhLjsYYY0wIS47GhAgawX+niLwoIkld+LMedb+4jIjc296fJY43ReR0d/seEdklIms7oW53iMiwoO0VIjI2wmPdLSI/7GidjDlV7KscxoQQkTpV7e+urwW2q+pTYewnOL9TEQ22LSKfA+NVdX879skCrlXVBe72Zzjf8y0LKRevzri67anPVuCnqlrcnv1aOVYS8K6q9qbRcEw3Zj1HY9r2DjAamueL2+ku97rvpbk9tWU4AymcIyJzRWSHW+6Xbrk4EVntvrdDRALJbLWIzBGRe4BhQIGIFIhIjoj8KlAJEblTRFpK0P8EbHLL5OGMdJQvIgtE5BERWS4irwFr3Lq+486zWCIiVwQd/9/den0kIo+LMwfjeGCt24vuJyJbRWS8W/6ENrrv14nIL9zjFIrIYAB1Zkv4XER6y2w5pruL9qgHttgSawvu3G84Q1Ztwpkb7hKc0TU8QH+ckYXG4cxb6AcmuPsMA/4PZ2i+eOBNYKa7/+tBP+MM93U17vBqwOfAIHfdgzMSSB93+z3gwhbq+gWQHLQdfIxHgO1AP3c7Cejrrp8LFLvrN7rHD8yFN9B93YrTkyV4u7U2umUUd6QlYCmwKGj/B4CfRPvf1xZbwlms52jMifqJyIdAMU4SWAlcBWxUVa+q1gEvAYEpvr5Q1UJ3/VJgq6rWqHMZcy1wNc64kKNE5L9F5AagzZlPVNWLk3Smicjf4yTJHS0UHaiqh9s4VL6qHnHX+wC/FZEdOMMJBu4fXgusUncuPFU92Fbd2mgjQCPwJ3d9O84fDwHVOInVmJgXS1NWGRMrjqjqxcFvnGRyaG9w0ZYKqGqtiFyEM2HtfOBWnHEy27IC+A/gM2BVK2V8InKatn6fM7huC3DGnrwI55bK0aA6t+fhg7Zi0aSqgWMd4/j/Y/oCR07cxZjYYz1HY8LzNjBTRJJExAPMwrkfGaoIuEZEBolIHDAXeEtEBgGnqeoG4EEgo4V9DwPJgQ1VLcKZpucfgRdaqddunPuM4RgAVLqJ9AdAnPv+a8C8wJOyIjKwpfqcrI1h/PzzgJ1h1tWYqLLkaEwYVLUE5/7gNpzksEJVP2ihXCXwM6AAd7YTVd2EM2v5Vvdy7Wq3TKjlwGYRKQh6bz3OU561rVTtFWBimM1YBmSLSCFOovK6df5fIB8oduv3U7f8aiAv8EBOGG08mStxZk4wJubZVzmMiWEi8ifgV6r6RiufDwXWqOp1p7Zm7SMi44D7VPUH0a6LMeGwnqMxMUhEzhCRv+Lc/2wxMUJzL+63gUEAYtggnMvJxnQL1nM0xhhjQljP0RhjjAlhydEYY4wJYcnRGGOMCWHJ0RhjjAlhydEYY4wJ8f+tAMOtcBr8kQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}