PythonDataBasics_Bivariate_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 Bivariate 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 Bivariate Statistics\n",
    "\n",
    "Here's a demonstration of calculation of bivariate 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",
    "#### Bivariate Statistics\n",
    "* Covariances\n",
    "* Pearson Product Momment Correlation Coefficient\n",
    "* Spearman Rank Correlation Coefficient\n",
    "\n",
    "I have a lecture on these bivariate statistics available on [YouTube](https://www.youtube.com/watch?v=wZwYEDqB4A4&list=PLG19vXLQHvSB-D4XKYieEku9GQMQyAzjJ&index=21&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": 11,
   "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.stats import pearsonr          # Pearson product moment correlation\n",
    "from scipy.stats import spearmanr         # spearman rank correlation"
   ]
  },
  {
   "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": 12,
   "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": 13,
   "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": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#df = pd.read_csv(\"2D_MV_200wells.csv\")                      # read a .csv file in as a DataFrame\n",
    "df = pd.read_csv('https://raw.githubusercontent.com/GeostatsGuy/GeoDataSets/master/2D_MV_200wells.csv')\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": [
    "#### Calculate the Rank Transform\n",
    "\n",
    "Let's demonstrate the rank transform, used to calculate the rank correlation coefficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5YAAAGHCAYAAADLHs1pAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABy+0lEQVR4nO39e3xU133v/78+I6HLoBu6IEDiagtjIxtfkLHj2MYhaRwnDXWapE56GjdN4/b80qY5ab5N0p7vaU7O8feX06ZpT9vz7Tm4SeO0dhInLsFxcOJLjC+xwcIYDBiwwIAZIXSDQQhdQJr1/UMjPIiRNKO57Zl5Px8PHkhbe2bWHo323p/1+ay1zDmHiIiIiIiIyEz5Mt0AERERERERyW4KLEVERERERCQhCixFREREREQkIQosRUREREREJCEKLEVERERERCQhCixFREREREQkIYWZbkAiamtr3ZIlSzLdDBER8aBXX321xzlXl+l2ZCNdX0VEJJqprq1ZHVguWbKE7du3Z7oZIiLiQWZ2NNNtyFa6voqISDRTXVtVCisiIiIiIiIJUWApIiIiIiIiCVFgKSIiIiIiIglRYCkiIiIiIiIJUWApIiIiIiIiCVFgKSIiIiIiIglRYCkiIiIiIiIJUWApIiIiIiIiCVFgKSIiIiIiIglRYCkiIiIiIiIJUWApIiKSJcxsoZk9a2b7zGyvmf1JeHu1mT1lZm3h/+dEPOarZnbQzA6Y2fsz13oREcllhZlugIiI5Ib2QIC21lb6u7spq6ujqaWFhsbGTDcr14wAf+qc22Fm5cCrZvYU8LvAM865b5jZV4CvAF82s6uAe4CVwALgaTNb7pwbzVD7RUQkybxy/VVgKSIiCWsPBNi/aROrqqqorq/nZH8/uzZtgvXrFVwmkXOuA+gIf33GzPYBDcB6YG14tweBLcCXw9t/4JwbBg6b2UHgRuDl9LZcRERmYrqg0UvXX5XCiohIwtpaW1lVVUVtRQU+n4/aigpWVVXR1tqa6ablLDNbAlwHbAPqw0HnePA5N7xbA3As4mGB8DYREfG48aCxeWCAu+rraR4YYP+mTbQHAhf28dL1VxlLERFJWH93N9X19Rdtqy4ro7+zM0Mtym1mVgY8CnzBOddnZpPuGmWbm+Q57wPuA1i0aFEymikiIgmIDBqBsaAR2NPaeiEbOfH6297Tw4EDB9gaDj7TWRarjKWIiCSsrK6Ok/39F2072d9PWV1dhlqUu8xsFmNB5UPOuX8Pb+40s/nhn88HusLbA8DCiIc3AsejPa9zboNzbrVzbnWdfm8iIhnX391NdVnZRduqy8ro7+6+8H3k9be9p4f927axtK+PDy5cGDXDmUoKLEVEJGFNLS3sCgbp6esjFArR09fHrmCQppaWTDctp9hYavLbwD7n3LcifvQYcG/463uBTRHb7zGzYjNbCjQBr6SrvSIiMnOxdNpGXn8PHDjAMp+PvlCIxuXLsXPnmNXWxsP338+WjRtTHmAqsBQRkYQ1NDayYv169vj9bO7sZI/fzwpN3JMKtwC/A7zHzHaG/90FfAN4n5m1Ae8Lf49zbi/wCPAG8HPgc5oRVkQkO0QGjce6utj0zDP8y+OP09vbeyFIjLz+bg0E6KuspHHNGgAC27Zxjc/HFWZpyV5qjKWIiCRFQ2OjAskUc869SPRxkwDrJnnM/cD9KWuUiIikRENjI6xfz7NPPEH7s89y/dy5/O5734uVlFw082vk9bdhYICaigp2v/wyS2fP5hxQER6nOXF8ZrIpsBQREREREfGghsZG6mprueODH7wwiQ9wIUgELixHMmjGlq4u1i5ZwtlgkKHiYnYPDrKiuRlI/aR6CixFREREREQyZLq1Kiebef3Y/v2MHj9+0RqWzzjHc4ODHAmF6HSO1WvW0FBbC6R+Uj2NsRQREREREcmAWNaqnGwSn+DJk5esYblu6VJqamq45z//Z8qamiguKkrbpHoKLEVERERERDIgcq3K8eBwVVUVbeEyV5h85vWKOXMmXY5kfFKf5wYH+dsnnuBfXnyRk0VFKT0WlcKKiIiIiIhkwGRlrpFjIccn8dnT2kp/ZydldXWsuP12ClpbOdnff9HYy4nlrlXDw9x7661Ul5Vxsr//okl/kk2BpYiIiIiISAaMl7lGBod7jx7lcHs7j2/YcNGYy2jB4K5Nm1gF7wSOwSArbr8duDgbCqR8ZlgFliIiIiIikjbTTVaTT5paWi4KDvcePcrzv/oV77/lFi4PT8gzWZZxskzm+H6xZEOTSYGliIiIiIikxfhkNZEzmaayPNPrJgaHh9vbef8tt7B86VJg+izjVGtIR8uGpnJmWAWWIiIiIiKSFukuz0ynmWZiI4PDxzds4PIkZRknZkMnlsommwJLEREREZEUUMnnpdJdnpkuycrEJjPLOF2pbLIpsBQRERERSTKVfEaX7vLMdElWJjbZWcapSmWTTetYioiIiIgkWSzrE+ajydZkbGppyXTTEtLf3T3pmpLxGF9/co/fz0P79/P9118ncOYMba2ttAcCyWxy0iljKSIiIiKSZLla8pmodJdnpsvETGxvTw97d+1i3/AwZRs3xlUGPb7f6PHjfKCxMS1rUCaDAksRERERkSTL1ZLPZEhneWa6RJawuqEh9j7/PKfNuOfWWykdGIg7KMzGSY4UWIqIiIiIJFm6Z+SUzIrMxL76wgtcPmcOq1etoqG2FiCmoDBysqd9r73GZatXQ0THhNcz3gosRURERESSLFdLPmVy45nY/u5u7qqvx+d7Zzqb6YLCiZM9VRQXs+OFF/DdfvuF4NTrGW8FliIiIiIiKZCLJZ8yvanKoCdbgmZi6evKVavYvWUL23ftYv4dd2RFxluzwoqIiIiIiCTJZDPfljU28uJ3v8vp554jtHMnp597jhe/+13aA4FLZpWtqa1l5W23cXBoiM2dnezx+1nh4Yl7QBlLERERERGRpJmsDPqlJ56g5OBBbqmpobq6mpODg/zq4EFeeuIJ6qJkOa2khBve/37W3n13Bo8mdspYioiIiIiIJFFDYyNr776b6+66C4DXNm/mxUcf5YrCQmr9fnxm1Pr93DxnDkd27MiJ9T2VsRQRERERkZw22djGVL9m5IQ8b54/z8m33qJ39mxqwpnJIjN8zuXEZE8KLEVEREREJGdNDPBO9vfHva7kTEyckOfyK69kZOdO2g4fZs7VV9M/NMTO3l6W3HorkP2TPSmwFBERERGRnDUxwKutqIhpXclE9Xd3U11ff+H7G264gafb2+k9eJAjp08zWl7O8MqVvP/OO1PWhnRSYCkiIiIiIjmpPRCg9ckn6TOjoqqKpqYmGmprp11XMhmiLTsyWlhI8eLF+BctwhUVUVhentI2pJMCSxERERERyTnjJbBriou5xudjZHiYXdu2wZo1FBcVUVZXl9LXb2ppYdemTawCqsvK2L5rF8v8fq6+6y5qamsB6OnrS3nmNF0UWIqIiIiIRJGJCV8kecZLYG3VKo5u28bS2bO5urSUl3ftoqypiRW3357S1584Ic/BoSE+ddttF4JKIC2Z03RRYCkiIiIiMkGmJnyR5Bkf4+jz+WDNGo62tXE2GORgKMQ9afo9Rk7IU7ZxIzYwcNHPT/b3pzxzmi4KLEVEREREJsjUhC+5JpNZ38gxjjW1tdTU1tLT18eQ35+R3+HE0tiT/f3sCgZTnjlNFwWWIiJZTGVaIiKpMXFGT8itssV0yHTWN9WBXLzX4FxYq3IqCixFRLJUpi/YIiK5LNqMnrlUtpgOmc76pjKQm+k1ONvXqpyKAksRkSyV6Qu2ZIaZfQf4ENDlnGsOb/shcEV4lyog6Jy71syWAPuAA+GfbXXO/WF6WyySnXK9bDEdvJD1TTSQmywrqWvwpRRYiohkKS9csCUjvgv8I/C98Q3Oud8a/9rM/gY4HbH/IefctelqnEiuyPWyxXTIZNY3GUNFpspK6hp8KQWWIiJZSmVa+ck593w4E3kJMzPg48B70tookRyVy2WL6ZCprG+yhopMlZXUNfhSvkw3QEREZqappYVdwSA9fX2EQiF6+vrYFQzS1NKS6aZJ5twKdDrn2iK2LTWz18zsOTO7NVMNE5H809DYyIr169nj97O5s5M9fj8r0jAPQGRA6PP5xgLCqiraWlvjep7+7m6qy8ou2lZdVkZ/d/dF1+Duri6ef+YZ/u1nP6O7p4f2QCCZh5M1Up6xNLMCYDvQ7pz7kJlVAz8ElgBHgI87506F9/0q8BlgFPi8c+4XqW6fiEi2UpmWRPEJ4PsR33cAi5xzvWZ2A/ATM1vpnOub+EAzuw+4D2DRokVpaayI5L5MZH2TVaY6VVZy/Br83M9/zlvPPMPy+np+c906SktK8nYivXSUwv4JYxMHjP9GvgI845z7hpl9Jfz9l83sKuAeYCWwAHjazJY750bT0EYRkaykMi0ZZ2aFwEeAG8a3OeeGgeHw16+a2SFgOWMdvhdxzm0ANgCsXr3apaPNItPRkkoyE8kqU52ulLehsZGamhpu/9CHLnqtfJ3EJ6WlsGbWCHwQ+OeIzeuBB8NfPwj8RsT2Hzjnhp1zh4GDwI2pbJ+IiEgOeS+w3zl3oQbLzOrClUOY2TKgCXgrQ+0Ticv4OLnmgQHuqq+neWCA/Zs25W2ZocQuWUNFYinlnapcNpXaAwG2bNzI4xs2sGXjRk/8XaQ6Y/l3wJ8B5RHb6p1zHQDOuQ4zmxve3gBsjdgvEN4mIiIiYWb2fWAtUGtmAeAvnXPfZqzq5/sTdr8N+LqZjTA2zOQPnXMn09lekZnScg4yU8kcKjJdZVAmJvHx6jrWKQsszWx8ja1XzWxtLA+Jsu2SUhyNARERkXzmnPvEJNt/N8q2R4FHU90mkVTI1+UcVP6bHOkaKpKJmW+92umSyozlLcCHzewuoASoMLN/AzrNbH44Wzkf6ArvHwAWRjy+ETg+8Uk1BkREREQk9+Xjcg5ezUR5jZeC70xMpOfVTpeUBZbOua8CXwUIZyy/5Jz7D2b218C9wDfC/28KP+Qx4GEz+xZjk/c0Aa+kqn0iIiIi4l2ZWgMxk7yaifISLwbf6Z5Iz6udLumYFXaibwCPmNlngLeBjwE45/aa2SPAG8AI8DnNCCsiIiLiTanOGuXjkkpezUR5iYJv73a6pCWwdM5tAbaEv+4F1k2y3/3A/elok4iIiIjMTLqyRvm2pJJXM1FeMN6R8dKjjxJqbOSKK66gobYWyL/g26udLpnIWIqIiIhIFlPWKDW8monKtMiOjDkLF1Jx+jT7t22DNWtoqK3Ny+Dbi50uCixFREREJC4q2UwNr2aiMi2yI8OWLyewbRvLfD4OHDhAcVFRWoJvL00Y5FUKLEVEREQkLirZTB0vZqIyLbIjo6a2FtasIfDmm2w9dgzfmjUpD769OGGQF/ky3QARERERyS5NLS3sCgbp6esjFArR09fHrmCQppaWTDdNctB4R8a4mtpaGpqbuWLtWgBe27yZLRs30h4IpOT1IzOmPp9vrPS7qoq21taUvF62UsZSREREROKiks3ckC3lndHGnm45coRR52guLU15FlGl37FRYCkiIiIicVPJZnbLpvLOaB0Zbu5c1vn9aZlASqXfsVFgKSIiIiKSZ7JtZt+JHRmPb9hAdVnZRfskM4sYmc0dNGNLVxdrlyzRbL1TUGApIiIiIpIHIoOlAzt2sLKlBSKycNlU3pnKLGK0bO4zzvHc4CDFZ8+q9HsSCixFRERERHLcxGBptKSEvc8/z9Vr147NtEp2lXemcs3PaNncdUuXssfvZ+3ddyf8/LlKgaWIiIiISI6bGCytXrWK7c89x95du3j3HXdkXXlnKieQmm6ynmyZ9CjdFFiKiIiIiCfohj11JgZLDbW1hG69lR9s305fls7sm6oJpKYqs82mSY/STetYioiIiEjGjd+wNw8McFd9Pc0DA+zftCllaxPmm4lrQQKUlpTQ8mu/xofuu4+1d9+d94HRuKnWadWalpNTYCkiIiIiGacb9tSaKljKF+2BAFs2buTxDRvYsnHjpJ0WDY2NrFi/nj1+P5s7O9nj97MinJHs7+6OPhttd3c6DsHTVAorIiIiIhmnRehTK5VjErNBvCWsk5XZak3LySmwFBEREZGM0w176qVqTGI2SNa6namcjTYWXh6HrMBSREREJEd4+aZzOpm+YZfclqyMeCYzv16fOEiBpYiIiEgO8PpN53TyvVRTUiuZGfFMZX6TlXVNFQWWIiIiIjnA6zedscjnUk2ZWqLZ+FzIiHt9HLICSxERSZpsLsMTyXZev+lMFZ13cl8ysvG5kBH3+jhkBZYiIpIU2V6GJ5LtvH7TmQo67+SHZGXjsz0j7vWsqwJLERFJilwowxPJZl6/6UwFnXeyQ6JZ5XzNxk/k9ayrAksREUkKXfhFMsvrN52poPOO9yUjq5yP2fjJeDnrqsBSRESSQhd+kczz8k1nKui8kx6JZByTkVXOx2x8NlJgKSIiSaELv4ikm847yTcxiCxrbOT0K6/MOOOYjKxyPmbjs5ECSxERSQpd+EUk3XTeSa5oZas/euAB3t3cPOOMY7KyyvmWjc9GCixFRCRpdOEXkXTTeWfmJmYnu3t6uGNC2eplIyN0Hz8OS5deeFw8GUdllfOHAksRERERkSyRrHU7o2Un/+3ZZxlctw4isov19fW80dV10WPjyTgqq5w/FFiKiIhkETP7DvAhoMs51xze9jXgs0B3eLc/d85tDv/sq8BngFHg8865X6S90SKSFMlctzPapDrXz53Ljt27Wbhu3YX9SufPp6u3l56+vhlnHJVVzg8KLEVERLLLd4F/BL43YfvfOue+GbnBzK4C7gFWAguAp81suXNuNB0NFZHkSua6nRMn1ent6YHz53m2tZU5wJVXX42VlHCsoIBbPvtZ9gQCyjjKlBRYioiIZBHn3PNmtiTG3dcDP3DODQOHzewgcCPwcqraJyKpk8x1OyMn1ent6SGwbRsLfT5uaWnhFPDdp59m2bp13DSeDW1pSdJReFuySo3zkS/TDRAREZGk+CMze93MvmNmc8LbGoBjEfsEwttEJAuNB4ORZrpuZ1NLC7uCQXr6+gi8+SYVPh9vhUK8a/Vq1q9bx6c/9CFqamryKqgaLzVuHhjgrvp6mgcG2L9pE+2BQKablhUUWIqIiGS/fwIuA64FOoC/CW+3KPu6aE9gZveZ2XYz297d3R1tFxHJsMhgMBQK0dPXx65gkKYZZBMbGhtZsX49e/x+fnbsGIcrKlixZg0NtbVAOBOaZ+eCyFJjn883VmpcVUVba+uFfdoDAbZs3MjjGzawZeNGBZ0RVAorIiKS5ZxzF+rgzOwB4PHwtwFgYcSujcDxSZ5jA7ABYPXq1VGDTxHJrGTPsBo5qU7zwEDCa01mu+lKjZM5eVIuUmApIiKS5cxsvnOuI/zt3cCe8NePAQ+b2bcYm7ynCXglA00USQqNf0vNDKtaa3JM5LjTcZEBdjInT8pFCixFRESyiJl9H1gL1JpZAPhLYK2ZXctYmesR4A8AnHN7zewR4A1gBPicZoSVbKVsUWxmEnxrrckx0wXYyZw8KRcpsBQREckizrlPRNn87Sn2vx+4P3UtEkkPZYuml0jwrbUmpw+wp8to5jsFliIiIiLiecoWTU/Bd+KmCrCTXTKca6XdCixFRERExPNSlS3KpZt7Bd+plcyS4Vws7VZgKSIiIiKel4oJZnLt5l6lmqmXrJLhXMwuK7AUERERySPZmqFLxQQzuXZzr9lds0cuZpcVWIqIiIjkiWzP0CV7gplcu7n3wuyu2dpxkW65mF1WYCkiIiKSJ3ItQ5eoXLy5z+TsrtnecZFOuZhdVmApIiIikidyLUOXqFy8uc+kbO+4SGe21QvZ5WRTYCkiIiKSJ3IxQ5eIXLy5z6Rs7rjIRLY119YOVWApIiIikieyLUOXjgxSrt3cZ1I2d1xke7bVCxRYioiIiOSJbMrQabxe9kmk4yLTk/5kc7bVKxRYioiIiOSRbMnQKYOUfWbaceGFToRszrZ6hQJLEREREfGcXM0gZTozl2oz6bjwQidCtpWJe5ECSxERERHxnFzMIHkhM5cKiQbLXuhEyKYyca9SYCkiIiIinpOLGSQvZOaSLRnBslc6EbKlTNyrFFiKiIiIiOfkYgYp1sxcNpXLJiNYzsVOhHykwFJEREREPCnXMkixZOayrVw2GWWsudiJkI8UWIqIiIiIpEEsmblsK5dNVhlrrnUiZEoms90KLEVEREQkZtlUpuk1sWTmvDCRTTxUxuodmc52K7AUERERkZhk+sY1F0yXmfPKRDaxUhmrd2Q6263AUkRERERiykRm+sY1H2RjBlBlrN6Q6Wy3AksRERGRPBdrJjLTN675IB0ZQJUz56ZMZ7sVWIqIiIjkuVgzkZm+cc0XqcwAqpw5d2U6263AUkRERCTPxZqJzPSNqyRO5czvyLXMbabHu6YssDSzEuB5oDj8Oj92zv2lmVUDPwSWAEeAjzvnToUf81XgM8Ao8Hnn3C9S1T4RERERGRNrJjLTN66SuFwvZ441WMzVzG0mx7umMmM5DLzHOddvZrOAF83sCeAjwDPOuW+Y2VeArwBfNrOrgHuAlcAC4GkzW+6cG01hG0VERETyXjyZSE3Ukt1yuZw5nmBRmdvkS1lg6ZxzQH/421nhfw5YD6wNb38Q2AJ8Obz9B865YeCwmR0EbgReTlUbRURERESZyHySSDmz10tH4wkWcz1zmwkpHWNpZgXAq8DlwP9yzm0zs3rnXAeAc67DzOaGd28AtkY8PBDeJiIiIiIppkxkfphpJ0I2lI7GEyzmcuY2U1IaWIbLWK81sypgo5k1T7G7RXuKS3Yyuw+4D2DRokXJaKaIiEjahOcg+BBwK2NDPwaBPcDPnHN7M9k2ySyvZ4Mkd8ykEyEbSkfjCRY1EVXypWVWWOdc0My2AHcCnWY2P5ytnA90hXcLAAsjHtYIHI/yXBuADQCrV6++JPAUERHxKjP7GvDrjA0D2cbYNbAEWA58Ixx0/qlz7vVMtVEyIxuyQZK9ktFpkQ2lo/GOFVb5d3JNGVgm0qtqZnXA+XBQWQq8F/gfwGPAvcA3wv9vCj/kMeBhM/tW+LWagFdmeFwiIiJe1Oqc+9okP/tWeHiIynHyUDZkgyQ7JavTIhtKR+MNFlX+nVyTBpZJ6FWdDzwYHmfpAx5xzj1uZi8Dj5jZZ4C3gY8BOOf2mtkjwBvACPA5zQgrIiK5xDn3MwAzW+KcOxL5MzNrcc618k4lT1Rm9h3GOn27nHPN4W1/zdg1+xxwCPh0uGN3CbAPOBB++Fbn3B8m74gkUeOZpJcefZRQYyNXXHEFDbW1gPeyQZKdktVpkS2lowoWM2eqjGVCvarhgPO6KNt7gXWTPOZ+4P4p2iQiIpIL/t3Mft051w5gZrcD/whcHcNjvxve93sR254CvuqcGzGz/wF8lbEZ1wEOOeeuTVbDJXkiM0lzFi6k4vRp9m/bBmvW0FBb67lskGSnZJWwxpsN1Jjh/DNpYDneqzrFz7uYpldVREREovoD4Cdm9uvA9cD/A9wVywOdc8+HM5GR256M+HYr8NEktVNSKDKTZMuXE9i2jWU+HwcOHKC4qMiT2SDJPsksYY01G6gxw/nJN9UPzexeM9thZmfD/7ab2afS1TgREZFcFC55/TzwJPA14H3OuWNJevrfA56I+H6pmb1mZs+Z2a1Jeg1Jgv7ubqrLygCoqa2lcc0a+ior2RoIsMfvZ4VuwiUJmlpa2BUM0tPXRygUoqevj13BIE0tLSl7zchOE5/PN1Z+W1VFW2vrJfu2BwJs2biRxzdsYMvGjbQHAilrl6TWVGMsPwV8AfgisIOx5UCuB/7azHDOfW+yx4qIiMilzOynXLyUlh84DXw7fG39cILP/xeMzVPwUHhTB7DIOddrZjcwliVd6Zzri/JYLeeVZhMzSTW1tbiiIt51442svfvuDLdOvCAZ5aSZmP001vJbZTZzy1RjLP9/wN0TJhf4pZn9JvADLh7bISIiItP7Zqqe2MzuZWxSn3XOOQfgnBsGhsNfv2pmhxibhG/7xMdrOa/0y5bJUCQzkhl0pXtCm1jLbzUbcm6ZKrCsmDhjHYBz7oiZVUTZX0RERKbgnHsuFc9rZncyNlnP7c65gYjtdcBJ59yomS1jbCmvt1LRBomf1tGTqWRz0BVrp8lUmU1N/pN9pgosB2f4MxEREZmCmX2EsbWd5zI21MQA55ybtuPWzL4PrAVqzSwA/CVjs8AWA0+ZGbyzrMhtwNfNbAQYBf7QOXcy+UckM6WlEWQyyZrNNRNi7TSZLLM57POpRDYLTRVYXmlm0daoNGBZitojIiKSD/4K+HXn3L54H+ic+0SUzd+eZN9HgUfjfQ0RybxkzOaayaxfLJ0mk2U2Q0VFWZutzWdTBpZpa4WIiEh+6ZxJUCniJSpVTK1Ex+Bmw8Q4k2U2X9u8+cKMyeOyJVubz6Zax/Jo5PfhcZVTBaIiIiISm+1m9kPgJ4Qn1wFwzv17xlokEodsCFqyXaJjcLNljGa0zGZbEtfelPSZNlA0sz8Avs7YuMrxWeIcKocVEZEJlMGIWQUwAPxaxDYHKLCUrJAtQYvXxHuOTGQMbjaP0dSMydkplgzkl4CVzrmeVDdGRESylzIYsXPOfTrTbRBJRDYHLckSb5CY7nNkMsZoZopmTM5OsQSWhxjrVRUREZmUMhixM7MS4DPASqBkfLtz7vcy1iiROGRz0JIMMwkS032OzPasn2ZMzj6xBJZfBV4ys21cPA7k8ylrlYiIZB1lMOLyr8B+4P2MDTf5bUCT+UjWyPagJVEzCRLTfY5MRtZPwxskHrEElv8H+CWwGwiltjkiIpKt8j2DEafLnXMfM7P1zrkHzexh4BeZbpRIrPK9VHEmQWImzpGJZP00vEHiFUtgOeKc+2LKWyIiIlkt3zMYcTof/j9oZs3ACWBJ5pojEr98LlWcSZA48Rx58OhRntmzB/+yZWzZuNFz2UANb5B4xRJYPmtm9wE/5eJS2JMpa5WIiGSdfM9gxGmDmc0B/jPwGFAG/N+ZbZLEQyWC+W0mHWmR58hj+/dz8tAh3tPczMrFiz2ZDdTwBolXLIHlJ8P/fzVim5YbERGRS+RzBiNWZuYD+pxzp4Dn0fU066hEUGbakTZ+jtyycSMfaGz0dDZQwxskXtMGls65peloiIiISD5wzoXM7I+ARzLdFpkZlQgKeGeNyVRlzzW8QeI1aWBpZu92zr04xc8rgEXOuT0paZmIiEjuesrMvgT8EDg7vlHDTLKDSgQlmngCvGRlA1OZPdfwBonXVBnL3zSzvwJ+DrwKdDO21tblwB3AYuBPU95CERGR3DO+XuXnIrZpmEmWUImgTBRvgJesbGCqs+ca3iDxmDSwdM79p/DEAh8FPgbMBwYZW2fr/0yVzRQRkczT5CLepWEm2S3fSgR1LplevAFesrKByp6Ll0w5xjI8scAD4X8iIpIlNLmIN2mYSW7IpxJBnUtiM5MALxnZwGzLnquTIrfFMiusiIhkGU0u4lkaZpIj8qVEUOeS2GQqwMum7Lk6KXKfAksRkRyk8ihv0jATyTY6l8QmUwHeVNlzr2UH1UmR+xRYiojkoGwrj8onGmYi2UTnkthksjw6Wvbci9lBdVLkvmkDSzPbDvwL8HD4YigiIh6XTeVRIuJdOpfELl3l0bFkIr2YHVQnRe6LJWN5D/BpoDUiyHzSOedS2jIREZmxfJpcRERSR+eS1IunZDXWTKQXs4PqpMh90waWzrmDwF+Y2f8NfAj4DhAys+8A/1OLOYuIeFO+TC4iIqmlc0nqxFuyGmsm0ovZQXVS5L6Yxlia2TWMZS3vAh4FHgLeDfwSuDZVjRMREclFGmYiIhB/yepUmcjIzOegGVu6uli7ZImnsoPqpMhtsYyxfBUIAt8GvuKcGw7/aJuZ3ZLCtomIiOQqDTMRkbhLVifLRA77fJdkPp9xjucGByk+ezZl2UGvzTwrmRVLxvJjzrm3IjeY2VLn3GHn3EdS1C4REZGcpWEmIgLxl6xONk4xVFR0SeZz3dKl7PH7WXv33SlpuxdnnpXM8sWwz49j3CYiIiIxCg8z+RvgrxkbZvJRoI+xYSYikgeaWlrYFQzS09dHKBSip6+PXcEgTS0tUfdvaGxkxfr17PH72dzZyR6/nxXr11PqHNVlZRftW11WRn93d8raHlnG6/P5xsp4q6poa21N2Wtmi/ZAgC0bN/L4hg1s2biR9kAg001Ki0kzlma2AlgJVJpZZGayAihJdcNERERylYaZiBeprDH9ZjKhTbRxim0ZmKzHizPPekE+Z3KnKoW9grHynCrg1yO2nwE+m8I2iYiI5DoNMxFPyeeb4UxLxoQ2mVjKw4szz3qBF9cQTZdJA0vn3CZgk5nd7Jx7OY1tEhERyXU/Bq6Psu2GDLRFJK9vhnNBJpby8OK6lF7IuudzJneqUtg/c879FfBJM/vExJ875z6f0paJiIjkGA0zEa/K55vhXJHupTy8ti6lV7Lu+ZzJnaoUdl/4/+3paIiIiEge0DAT8aR8vhmWmfPSupSJZN2Tmen0YiY3XaYqhf1p+P8H09ccERGR3KVhJuJV+XwzLN4Wa9A306x7sjOdXsvkptNUpbA/BSZdqNk59+GUtEhERCRHJWOYSXityw8BXc655vC2auCHwBLgCPBx59yp8M++CnwGGAU+75z7RXKORnJJPt8Mi3fFE/TNNOueivHFXsrkptNUpbDfTFsrRERE8kMyhpl8F/hH4HsR274CPOOc+4aZfSX8/ZfN7CrgHsbGdS4Anjaz5c650QReX6bhhQlEZiJfb4bFu+IJ+maaddf44uSZqhT2uXQ2REREJNclY5iJc+55M1syYfN6YG346weBLcCXw9t/EF4n87CZHQRuBFSGmyJemUBEJJXS1XkST9A306y7xhcnz1SlsI845z5uZru5uCTWAOecuyblrRMREckhKRxmUu+c6wg/R4eZzQ1vbwC2RuwXCG+TFNGyHZLr0tl5Em/QN5Osu8YXJ89UpbB/Ev7/Q+loiIiISB5I9zATi7ItamBrZvcB9wEsWrQolW3KaSqrk1yXzs6TVAV9EzOulTfeyJ5AQOOLEzRVKex4z+dRM5vHWOmMA1qdcyfS1D4REZGckcJhJp1mNj+crZwPdIW3B4CFEfs1AscnadsGYAPA6tWrJ82qytRUVie5Lp2dJ6mYVCpqxvWVV1ihcvWETZWxBMDMfh/4L8AvGev5/Acz+7pz7jupbpyIiEguSeEwk8eAe4FvhP/fFLH9YTP7FmOT9zQBr8zwNSQGKqvLXdk6KVOypbvzJNmTSqlcPXWmDSyB/wu4zjnXC2BmNcBLgAJLERGR+CQ8zMTMvs/YRD21ZhYA/pKxgPIRM/sM8DbwMQDn3F4zewR4AxgBPqcZYVNLy3bkpnjGFeZ6AJrtnScqV0+dWALLAHAm4vszwLHUNEdERCR3JWOYiXPukvUvw9ZNsv/9wP0zaK7MkJbtyD2xZrnyYVbgbO88Ubl66kw1K+wXw1+2A9vMbBNjF7/1qIxGRERkxjTMRCS7xJrlamttZe7oKHv27qX/9GnKKiuZO28ebTlWZpnNnSfZnnH1sqkyluXh/w+F/43bFGVfERERiZ2GmUjOy6WS0FizXEcPHOD84cNcV1ZG9Zw5nBwc5LW9ezk+MJDuJssksj3j6mVTzQr7X9PZEBERkTyiYSaS03KtJDTWLFffqVPcVlBArd8PQK3fz+WDg+w/dSoDrZbJZHPG1ctimRW2DvgzYCVQMr7dOfeeFLZLREQk52iYiWSTRDKOuTbzZqxZrqrqaoLBIH0DA5SVlNA/NERwdJSq6uoMtVwkfWKZvOch4IeMzWD3h4xNY96dykaJiIjkKA0zkayQaMYxF2fejCXLtXD5cmaXlnL0xAkGg0FKKyuZvXgxCxcunPJxXpRLpcySHrEEljXOuW+b2Z+EF3Z+zsxStcCziIhIztIwE8kWiWYc83XmzaaWFvYfP86qlSsvLpltacl00+KSa6XMkh6xBJbnw/93mNkHgeOAPlEiIiIzpGEm4nWJZhzzdebNXJkYJp2lzMqM5o5YAsv/bmaVwJ8C/wBUAP8ppa0SERHJbRpmIp6WaMYxVwKsmciFiWHSVcqszGhumTawdM49Hv7yNHBHapsjIiKSFzTMRDwtGRnHXAiw8lW8HQszzTrm2iRP+c433Q5mttzMnjGzPeHvrzGz/5z6pomIiOSsi4aZmNl1aJiJeEhDYyMr1q9nj9/P5s5O9vj9rMhAFqk9EGDLxo08vmEDWzZupD0QSOvr56umlhZ2BYP09PURCoXo6etjVzBIU5SxouNZx+aBAe6qr6d5YID9mzbF9Lvq7+6muqzsom3VZWX0d6uAIxvFUgr7AGMLOf8fAOfc62b2MPDfU9kwERGRHKZhJuJ5mc44qkwyc+IpZU4k65ivkzzlqlgCS79z7hUzi9w2kqL2iIiI5DwNMxFNWDI9lUlmVqwdC4mMx8zXSZ5yVSyBZY+ZXcbYAs6Y2UeBjukeZGYLge8B84AQsME59z/NrJqxCQuWAEeAjzvnToUf81XgM8Ao8Hnn3C/iPSAREUk/3STHx8yWA/8E1Dvnms3sGuDDzjlVA+UBZeJiE0/AonNQ5iSSdcznSZ5yUSyB5eeADcAKM2sHDgP/IYbHjQB/6pzbYWblwKtm9hTwu8AzzrlvmNlXgK8AXzazq4B7GJt6fQHwtJktd86Nxn1UIiI5yos3T7pJnhENM8ljysTFJtaAReegzEo065jpkmtJnlhmhX0LeK+ZzQZ8zrkzsTyxc66DcGbTOXfGzPYBDcB6YG14tweBLcCXw9t/4JwbBg6b2UHgRuDleA5IRCRXefXmSTfJM6JhJnksXUs5ZLtYA5ZYz0Fe7JjLBco6yrhpA0szqwI+xVjpauH4RdA59/lYX8TMlgDXAdsYK/sZDzg7zGxueLcGYGvEwwLhbSIigncDON0kz8iMhplIbtCEJbGJNWCJ5Rzk1Y65XJHprKM6DbwhllLYzYwFfLsZGysZFzMrAx4FvuCc65vQO3vRrlG2uSjPdx9wH8CiRYvibY6ISNbyagCnm+QZmekwE8kBmrAkdrEELLGcg7zaMSeJU6eBd8QSWJY45744kyc3s1mMBZUPOef+Pby508zmh7OV84Gu8PYAsDDi4Y3A8YnP6ZzbwNjFmNWrV18SeIqI5CqvBnC6SY7fTIeZSG5Q6WByxXIOmq5jThmv7KVOA++IJbD8VzP7LPA4MDy+0Tl3cqoH2Vhq8tvAPufctyJ+9BhwL/CN8P+bIrY/bGbfYmzynibglRiPQ0Qk53k1gNNNcvySMcxEslumSwdzSSznoKk65pTxym5erebJR7EElueAvwb+gndKUx2wbJrH3QL8DrDbzHaGt/05YwHlI2b2GeBt4GMAzrm9ZvYI8AZjExh8TjPCioi8w8sBnG6S45bQMBORmcjlrNx056CpOuaU8cpuXq3myUexBJZfBC53zvXE88TOuReJPm4SYN0kj7kfuD+e1xERyScK4HLGjIeZiMyE17NyqQ56p+qYe23zZmW8sphXq3nyUSyB5V5gINUNERERySMzGmYiMlNezsqlK+idrGNOGa/s5uVqnnwTS2A5Cuw0s2e5+OKncSAiIiIzM9NhJiIz4uVxaJkOer2S8crlUuVUUzWPN8QSWP4k/E9ERESSY0bDTERmystZuUwHvV7IeHm9VDmdFGBnrykDSzMrAH7HOffeNLVHREQkH2iYiaSVV7Jy0Xgh6M10xivTWVuvUICd3aYMLJ1zo2Y2YGaVzrnT6WqUiIhIjtMwE0mrZGTlUpVJ8nLQmy6Zztp6hQLs7BZLKewQY0uGPAWcHd+oi5+IiMiM/QQNM5E0SyQrl8pMkhdKUTPNC1lbL1CAnd1iCSx/Fv4nIiIiCdIwE8lGycgkTZXxzHQpaqYpaztGAXZ2mzawdM49aGalwCLn3IE0tElERCRnaZiJZKNEM0kaOzc1ZW3HKMDObtMGlmb268A3gSJgqZldC3zdOffhFLdNREQkV2mYiWSVmWSSIjOUhw8dYm1DQ06OnUvW2NN8z9qCAuxsF0sp7NeAG4EtAM65nWa2NIVtEhGRPJRnU8wnfZiJmV0B/DBi0zLgvwBVwGeB7vD2P3fObU7ma0vuizeTNDFD+fTWrZwIBqkuL6ehthbIjbFzysQmnwLs7BVLYDninDttZpHb3GQ7i4iIxCvfbs5SMcwk/DzXwoVxnO3ARuDTwN86576ZjNeR3DXdGMh4MkkTx2TOnzePsmCQtra2C4FlLoyd0yymIu+IJbDcY2afBArMrAn4PPBSapslIiL5JN9uztIwzGQdcMg5d3RCx7BIVLF07sSTSZo4JnNBUxNvb91KR2cnoVAoZ8bOaRZTkXf4Ytjnj4GVjK2z9TBwGvhCCtskIiJ5pr+7m+qysou2VZeV0d/dPckjst7XGBtmEoSxYSZAMoeZ3AN8P+L7PzKz183sO2Y2J4mvIzkisnPH5/ONde5UVdHW2jqj5xsfkzmupraW2StX0lNby+bOTvb4/azIcEVCeyDAlo0beXzDBrZs3Eh7IBD3c0w8TsiNTKzITEyasTSzEuAPgcuB3cDNzrmRdDVMRETyRx5OMZ+yYSZmVgR8GPhqeNM/Af8t/Pz/Dfgb4PeiPO4+4D6ARYsWJaMpkkWSnXmLNibzWEEBH/njP/ZEFUKyyu/TNYtpno1Blyw1VcbyQWA1Y0HlBxgr2REREUm6ppYWdgWD9PT1EQqF6OnrY1cwSFNLS6ablioXDTMxs38gecNMPgDscM51AjjnOp1zo865EPAAY5nSSzjnNjjnVjvnVtflbkAvk0h25q2hsZEV69ezx+/3TIYyUrIytOk4zvEguHlggLvq62keGGD/pk0zyrCKpNJUYyyvcs5dDWBm3wZeSU+TREQk3+ThFPN/DPwF7wwz+QXw35P03J8gogzWzOY75zrC394N7EnS60gOSUXmzcuze84kQztZ1jDVx5lvY9Ale00VWJ4f/8I5N6LB/yIikkpevglNllQPMzEzP/A+4A8iNv9VeHIgBxyZ8DMRIP86d+Itv8/kzNWaIEiyxVSB5Soz6wt/bUBp+HsDnHOuYvKHioiISBQPMtZx+wJjJatXksQJ8ZxzA0DNhG2/k6znl8xJxxi7fOjcGRdvhjaTWcM8HIMuWWrSwNI5V5DOhoiIiOQBDTORuOXbOq+pEC0wXxFHhjaTWcN0TRAkkqhY1rEUERGP0QyBWUvDTCRuGmOXmMkC8xXr17P27rtjeo5MZg3zrUxZspcCSxGRLKPsRVbTMBOJm8bYJWYmgfnEzruyxkZ2vfJKxrKG+VSmLNlLgaWISJZR9iJ7aZiJzERktqy9p4e2tjY6Ojvpqa2lPRDQ3/004g3Mo3bevfIKlTfeyJ5AQFlDkUkosBQRyTLKXojkl/ExdnN7ezmxdy+XFxRQWVjI7IYG9qtaYVrxlrFO2nkXCMRcOiuSj3yZboCIiMQn2QuZi4i3NTQ2smL9era0t+POn6e/qopFN93E8qVLWVVVRVtra6ab6GlNLS3sCgbp6esjFArR09fHrmCQppaWqPv3d3dTXVZ20bbqsjL6u7tn3Ib2QIAtGzfy+IYNbNm4kfZAYMbPJeJVyliKiGQZzRAokn8aGhtZetllvPdd78LneycvoGqF6cUz+U17IMDhQ4d4eutW5s+bx4KmJmpqaxPqvNO4eMkXCixFRLKMZggUyU9az3DmYpn8ZjwAXNvQwIlgkLJgkLe3bqV35UqOFRTMuPNO4+IlXyiwFBHJQpohUCT/qFrhYsledikyAKwuL39nkqT2dj7yx3884+fWuHjJFwosRURERLKAqhXekYry0sgAsKG2lobaWkKhEJs7OxN6j5VplnyhwFJE8kqye7hFRFJlsvOVzlmpKS9NVQCoTLPkC80KKyJ5Y7yHu3lggLvq62keGGD/pk2anU9E0ibW2UF1vppaKmZujXf22FiNz+q7x+9nc2cne/x+VmjiHslByliKSN7QBAoikknxlG/qfDW1VGQXU1lqrEyz5AMFliKSNzSBgoik08RS1t7eXm6PMVjU+WpqqSovVQAoMnMqhRWRvDHewx1JEyiISCpElrKuKSjA99xzvLBhA2+89BK9PT0X9pusfFPnq6mpvFTEe5SxFJG8oQkURCRdxktZ7dw52ltbuXb2bE7W1zNw8CCB0VFYs4aa2tpJg0Wdr6an7KKItyiwFJG84aWp+jU7rUhuGy9l3bttG0tnz6bC7+eGpUt5dOdOrvD5CLz5Jq6oaNJgcabnK51bRCRTFFiKSF7xQg93KtZfExFvGS9lHTx9mrKqKgBKZ82i4ZprOFxWxtZjx3jXjTdOGSzGe77SuUVEMkmBpYhImmm2R5HcN17KOmvWLPoGBxkxY9fZs7xrzRqKi4rwrVnD2rvvTuprZuO5RRlWkdyhyXtERNIsFeuvieS6WNd/9IrxyWW6m5t5oKODX42OsrylheJw+WuiayNGk23nFq3VKZJblLEUkYSpxzk+qVh/TSSXZWuJZ0NjI7/5+79P+5130tbayq7ubsr8/pSN7c62c0s2ZlhFZHIKLEUkIdl6w5dJmu1RJD7ZHoCkY2x3eyBAb28v//LMMyyvr+f6q6+mtKTE0+cWrdUpklsUWIpIQrL9hi8TvDQ7rUg2mC4AyfeqifEOvturqrjtve9l3+7dPPrMMzTccQfv8nAnX7ZlWEVkagosRSQh6nGeGS/MTiuSLaYKQFQ1MaGDr6KCunXruKqvjz1+v6ffA1VviOQWBZYikhD1OItIqk0VgKhqIns7+GZSvZHv2WkRL1NgKSIJUY9z5ugGS/LFVAHIa5s3Z2VQlUxe6uCL97wUT/WGstMi3qbAUkQSovGCmaEbLMk3kwUgXgqqMsUrHXypPi8pOy3ibQosRSRhGi+YfrrBEhnjlaBqOqmsMPBKB1+qz0vZWvIrki8UWIqIZCHdYImM8UpQNZV0VBh4oYMv1eclZadFvE2BpYhIFtINlsg70rVO5EwzjvlSYZDq81K2ZKdF8pUCSxERD4j3plU3WCLpk2jGMVsqDBIt1031eSkbstMi+UyBpYhIhs3kplU3WBKNmR0BzgCjwIhzbrWZVQM/BJYAR4CPO+dOZaqN2SjRjGM2VBgko1w3HeclL5T8ikh0CixFRDJsJjetWmpEpnCHc64n4vuvAM84575hZl8Jf//lzDQtOyWaccyGCoNkletGC/x0vhLJDwosRUQyLN6bVi01InFaD6wNf/0gsAUFlnFJNOMYbyYvE4FYqsp1db4SyR8KLEVEMizem9Z8mQhEZsQBT5qZA/6Pc24DUO+c6wBwznWY2dyMtjALJSPjGGsJZ6YCsVSV6+p8JZI/fJlugIhIvmtqaWFXMEhPXx+hUIievj52BYM0tbRE3b+/u5vqsrKLtlWXldHf3Z2O5oq33eKcux74APA5M7st1gea2X1mtt3Mtnfrs3TBePbweH8/33/9df513z72+P2sSFGgFxmI+Xy+sUCsqoq21takv1akeM9DsdL5SiR/KGMpIpJh8ZbJZcNEIJIZzrnj4f+7zGwjcCPQaWbzw9nK+UDXJI/dAGwAWL16tUtXm70sMnt424oVFzKVqSxNzdQMsqmaeGeq85XGXorkFgWWIiIeEM9Mh9kwEYikn5nNBnzOuTPhr38N+DrwGHAv8I3w/5sy18rskokyzkx2HKVixtXJzleVy5dr7KVIjlFgKSKSZbTUiEyiHthoZjB2fX/YOfdzM2sFHjGzzwBvAx/LYBuzSiayh+nqOEpXtnCy85XGXorkHgWWIiJZSGu5yUTOubeAVVG29wLr0t+i7JeJ7GE6Oo7SPUFQtPPVa5s3Z6TkV0RSJ2WBpZl9B/gQ0OWcaw5vm3SRZjP7KvAZxhZ1/rxz7hepapuIiIhkRjaNq8tU2XmqO468kC3UWHGR3JPKWWG/C9w5Ydv4Is1NwDPh7zGzq4B7gJXhx/y/ZlaQwraJiIhImo1nypoHBrirvp7mgQH2b9pEeyCQ6aZF1dDYyIr169nj97O5szOls8Gmkxdmak3VLLSp1B4IsGXjRh7fsIEtGzd69nMrkikpy1g65543syUTNk+2SPN64AfOuWHgsJkdZGwmu5dT1T4RERFJLy9kyuKVi2XnXsgWZttY8UytLyqSTdI9xnKyRZobgK0R+wXC20RERCRHZGopDblYpmeWnlgOfd1dd3k+OMvGThGRdPPK5D0WZVvUNbTM7D7gPoBFixalsk0iIiKSRF7IlElms4XZmvnzYqdINo1XlvyQ7sByskWaA8DCiP0agePRnkALOIuIiGSnTGfK5B2ZKvHN1syf1zpFsjVAl9yWysl7ohlfpBkuXqT5MeAeMys2s6VAE/BKmtsmIiIiKZSrk+FI7NI1cVCyJ9rx2mRDkQG6z+cbC9Crqmhrbc1Ie0QgtcuNfJ+xiXpqzSwA/CXwDaIs0uyc22tmjwBvACPA55xzo6lqm4iIF6msSfJBLk6GI7FLR+YvFdk8r0025MXSXJFUzgr7iUl+FHWRZufc/cD9qWqPiIiXqaxJJL/lS8dSOsqhU1Vu66VOEa+V5oqAdybvERHJa22trcwdHWXP3r30nz5NWWUlc+fNo83j445EJHH51LGU6sxfeyBA65NP0mdGRVUVTU1NNNTW5lw2T+OVxYsUWIqIeMDRAwc4f/gw15WVUT1nDicHB3lt716ODwxkumkikmLZOqHNTKUq8zceoK8pLuYan4+R4WF2bdsGa9ZQXFSUU9k8r5XmioACSxERT+g7dYrbCgqo9fsBqPX7uXxwkP2nTmW4ZSK5L9NlqBovlxzjAbqtWsXRbdtYOns2V5eW8vKuXZQ1NeVcNs9LpbkioMBSRMQTqqqrCQaD9A0MUFZSQv/QEMHRUaqqqzPdNJGclmgZajKCUo2XS47xAN3n88GaNRxta+NsMMjBUIh7crCsWMRrFFiKiHjAwuXLmV1aytETJxgMBimtrGT24sUsXLhw+geLyIwlUoaarLGRGi+XHJEBek1tLTW1tfT09THk9yuoFEkDBZYiIh7Q1NLC/uPHWbVy5cU3lhlaI00knTJZippIGWqyxkZqvFxyKEAXySwFliIiHqAbS8lXmZ4RNZEy1GSOjdR4ucRFO49WLl9OW2srr23enNPLuIh4gQJLEREPyPTkISKZkukZURPJcmlspPdEBuiZ7rQQyTcKLEVEppHqoE83P5LPMj0jaiLVAuksvVTnU/wy3Wkhkm8UWIqITCEdQZ9ufiSfeSHrN9My1HSVsKvzaWYy3Wkhkm8UWIqITCEdQZ9ufiSfZfuEK+kYG5nozLX5mun0QqeFSD7xZboBIiJe1t/dTXVZ2UXbqsvK6O/uTtprjN/8RNLNj+SLhsZGVqxfzx6/n82dnezx+1mxfj0AWzZu5PENG9iycSPtgUCGW5o5Mz0PjWc6mwcGuKu+nuaBAfZv2pQ372VTSwu7gkF6+voIhUL09PWxKxikSbNti6SEMpYiIlNIR493tmdsRBI1Meun0s+LzfQ8lO9l9pptWyS9FFiKiEwhHUGfbn5ELpbvAdFEU52Hpip1VZm9lnERSScFliIiU0hX0KebH5F3pCMgyqaxh5Odh4ApM7saYygi6aTAUkRkGgr6RNIr1QFRNpbaRjsPbdm4ccrMrsrsRSSdFFiKiIiIp6Q6IJqs1PbZJ56grrY2K7KYMH1mV2X2IpJOCixFRETEU6YLiBItY40WkA0ODdH+7LPc8cEPZk0WM5bMriouRCRdFFiKiHhENo35Ekm1yQKiZJSxRgvIduzezfVz52bVhEEqdRURL9E6liIiHpDv681JfmoPBOJeqzKyjNXn840FgFVVtLW2xvy60dY3fLOzkyuvvvqi/ZK9Zm2yTbYGqFcDYRHJbcpYioh4gJZXkHwz08xjMmaMjVZqu2zdOqyk5KL9smEGVZW6iohXKLAUEfEArTcniTKzhcD3gHlACNjgnPufZvY14LPAeOrtz51zmzPTynfMtDMlWTPGTgzI2gOBrCorVem8iHiNAksRiVu+3dCk43i13pwkwQjwp865HWZWDrxqZk+Ff/a3zrlvZrBtl5hpZ0qqxhUmYwbVdJ0bs3G5FBHJfQosRSQu2XZDk+iNXrqOV5NwSKKccx1AR/jrM2a2D2jIbKsmN9POlFQuoZFIWWk6z40qnRcRL1JgKSJxyaYbmmTc6CXjeGMJbrXenCSTmS0BrgO2AbcAf2RmnwK2M5bVPJXB5gGJdaYkGgCmIquYznOjSudFxIsUWIpIXLLphiYZN3qJHm88wa0m4ZBkMLMy4FHgC865PjP7J+C/AS78/98AvxflcfcB9wEsWrQo5e1MZWfKZMFjKrOK6Tw3qnReRLxIgaWIxCWbbmiScaOX6PFmU4ZXsp+ZzWIsqHzIOffvAM65zoifPwA8Hu2xzrkNwAaA1atXu9S3NjWdKVMFj6n8e0znuVGl8yLiRVrHUkTiEm39t13BIE0tLZlu2iXGb/QixXujl+jx9nd3U11WdtE2r6+NJ9nJzAz4NrDPOfetiO3zI3a7G9iT7ral01TrXKbq77E9EKC7p4d/+9nPeP6ZZ+ju6krpuVHrV4qIFyljKSJxyaaxgMno1U/0eLMpwytZ7xbgd4DdZrYzvO3PgU+Y2bWMlcIeAf4gE41Ll6kqFVLx9zieIb2jqorBdevYsXs33336aZatW8dNKQz2VDovIl6jwFJE4pYtNzTJCoITOV6VrEm6OOdeBCzKjzK+ZmU6TRU8Rvt73HLkCG7uXB7fsGFGk/lcVF5bUcHCdevo6etjj9+fFedJEZFkUWApIjkt00FwNmV4RXLBVJ05E/8eh30+Rp1jnd//zr6bNtFx4430BwIxzRybTROaiYikkgJLEZEZiGfJgkwHtyL5ZLrOnMi/xy0bN9JcWnrRZD5ze3t58YEH+Njtt8c0c6zK3UVExiiwFBGJUzoXQheR+MXamRMt29h9/DiXjYzEPHOs18vdU7Vup4jIRAosRUTipCVERHJDtGxjZ3c3V01R2hotUFvh0XJ3dYKJSDopsBQRiZPGVInkhmjZxq6CAm6YP/+i/cZLWycL1FasX8/au+/OzEFMQZ1gIpJOCixFROKkMVUiuSHaeMxbPvtZjr3yCtV9fZeUtmZboJbOTjCV3IqIAksRkTh5fUyViIyJJdiJNh6zff78qKWtr23enFXVCunqBFPJrYiAAksRkbhpCRER70sk2Jls8p9sq1ZIVydYtmVyRSQ1FFiKiExhsoyHlhAR8bZUBDvZVq2Qrk6wbBt3rrJdkdRQYCmSh3RRjU2yy7v0voukTyqCnWysVkhHJ1g2ZXJVtiuSOgosRfJMpi+q2RRcJTPjken3fap2ZcvvQ7JTpj5jqQp2VK1wqWzK5KpsVyR1fJlugIikV+RF1efzjV1Uq6poa21N+WuPB1fNAwPcVV9P88AA+zdtoj0QSPlrx6o9EGDLxo08vmEDrU8+yeDQ0EU/ry4ro7+7O+7nzeT7Ppls+H1IdsvkZ6yppYVdwSA9fX2EQiF6+vrYFQzS1NIS83NEng+2bNyov41JNDQ2jq3l6fezubOTPX4/KzyaAezv7qa6rOyibTM9r4vIxZSxFMkzmRwL4/We4olZxYriYna88AK+22+nobYWmHnGw4tjkLz++5Dsl8nPWKJlq16tMvCqbMnkZlPZrki2UWApkmcyeVH1YnAVaeJN8MpVq9i9ZQvbd+1i/h13JFTe5cWbGa//PiT7peozFmt5bSLBjjpeclM2le2KZBsFliI5JJabrUxeVL0YXEXq7+5msKCALXv30n/6NGWVlVQ3N/PqoUNsTnCiDi/ezHj99yHZLxmfsYnntbLGRk6/8krKM4nqeMlN2TgBk0i2UGApkiNiLdvK5EXVi8FVpEEzdrzwAtcUFeFOnqTrzTd5aWSE8g98gA/dd98l+8czKYkXb2a8/vuQ7JfoZyzaee1HDzzAu5ubU55JVMdL7sqWsl2RbKPAUiRHxFO2lamL6kyDq3TNKukzo/DsWfrffpul5eVU+/0ET55k+86dtAcCF73m9tZWWh94gMtGRlhYX0/p0BD7jx+fMmvitZsZLwa7klsS/YxFO69dNjJC9/HjsHQpvT09HG9r42wwyKuhUFLPDep4ERGJjwJLkRwRS9mWF5aWiCe4ag8EeOmJJ2h/9lmunzuXNVdfjQ0MpGwCjeJQiPK5czk7PMwbIyOU+v3ccN11nDp7lraIAL09EOBXDzzAhwsLWVxXR//QEIffeIOFV111yX6Zfr+n47VgV3JPIp+xaOe1+vp63ujqorenh8C2bSydPZuh4mI6nePlBx/EzZ1LqXMJ/82p40VEJD4KLEVyxHRlW9k2w+F4e4va2vhkfT0lPh+HW1tpXLOGVVVVKZlAo6yujqHBQd599dX4fGOrMfUMDFBYUkLrk09eCBC7e3qYOzrK4rlz8ZlR4fezFDjc0UF/SclF7ff6+50Nwa/kr2jntdL58+nq7WXvrl1c4/dzDtg9OMhlTU30vP46p7u7uWt8sq0E/+bU8SIiEjsFliJZZKogYLqyrVhKZb0UZIy3d+v589TOmYPPjKXA0bY2Vq5Zk5IJNJpaWnjoRz8itH07VQUFnCss5GBhIRQUsKahgXeHA8R/e/ZZGubM4eTgILV+PwBlJSV0dnRQduONF7XfyzNKZkvwK/kr2nntWEEBt3z2szz7ne8QNKOiqooVzc2cbGvj2poaXhwefmetWN75m/PS+U1EJBcpsBTJEtMFAdOVbU1XKpuJIGOqG73x9pZVVrK3q4sdBw+y/8gRus+epfKxxxhdvjymm8N4byYrGxtpO36c+sJCigoKONrRQXNDAytXrbpws3r93Lns6uhgczDI8lCIuZWVnPX7OVRSwofDi69nw4yS2RD8Sn6b6rzWHwjQPDBw4fPbvn0754qKKKusvPD48b85daKIiKSeAkuRLBFLEDBV2Va0krK9R49yuL2dxzds4PChQ6xtaEhbkDHxRu/g0aM89sQT+JctY/EVVzBoNlbKW1PDo5s3c8WJE9xXUEDnyAgvd3fT5fMx8vrrU06YE+/NZFtrK7+5ahV25ZUcb2tj8PRpTp89S1F9PTW1tRf2q29s5MiLL/LRK6/k7MAAz3V3s29khLVf/eqFzMjhQ4f44bZtzK+vp6mpiYba2oRnlEx2xiUbgl/JTzNZOmlo1ixePnWK1ddcc2Gf8b85daKIiKSeAkuRLBEtCBgcGqL1xRdjCjQib8IGh4Z46qWXeHPnTm677jpali+nqLubE8Eg1eXlNISDqIkZzWQGNZE3er09PQy+8QZ3FhZyOBikeWCALV1dPOMcJT09NIRCLCspob2/nzkVFfxGVRXbgH2vv84nPvIRnn3iCepqay9pW7w3k+Pvsc/nuxBIFldU8EZX14V92nt6eHrLFoqKi9k+MEDF7NksveIKVs+bR+/58xeC2fc3NHA2GKQqGOSNrVs5uXIlXQUFM55RMhUZFy2nIJky1flkpksnDTc3M3LiBMVFRYRCoYuGA7y2ebM6UUREUkyBpUiWmBgEtPf0sOOFF1hTVXVh7N9Ugcb4Tdiz4VlWZ/X1cd9111FbWcnh1lbKZs9m/ugobW1tFwLL8SAj0aAm2gLnr/7iF5T4fHRUVXHmzBmaZ8+mrLSUvadOUVtRwdolS3hucJCDgQAlg4MsWrCAof5+FoeP//WBAc709jI4NETb449Tv3QptefOca6oiJffeIOb7703rozceJbx6a1bmT9vHguamqiprb0wUUhPXx+DQ0PseOEFKk6e5AvXX09BcTG7zp6lqamJ+dXVbO7svDhgLi/neFsb7sQJtrS385E//uMZB4GpyLhoOQXJhOnOJ4ksndQeCEQtm21TJ4qISMopsBTJkGjBVn8gMGlGcGIQsH3XLiqdY+WqVZw6eZKOtjZOHzrE3z39NNfcdBMLly+/5DkaGhupq63ljg9+kK0vvcSSiElxXh8eJjg6Skdn5yW9/YkENRNvIvcePcrzP/gB88rKWFZRQcnwME/v2sXKVas4CZRVVtLb00P7m29y4NgxRsvLOVlayi/7+vCfP487e5ZCoD0YpDMU4jsPP8ycgQEW1ddz+uRJhk+fZmhkhCdmzeLypqaYbibH27i2oYETwSBlwSBvb93KGwsWsKO9nYKaGr7/+ut0Hz/OqqIizlRW8vyBA9RVVTF3zhza2tooXrmSsrq6i4LZmtpaamprWRkKca6z03Nlq6lYTkETpMh0pjufJPJZjxZobtm4kaMHDrDvrbdY19zM5YsXqxNFRCQFFFiKZEC08YVP/eAH3HbLLdw2ftMzISM4MQg4ODTEp267DYDAtm1UhkIsDgY5PTRE81tvMbu0NOr4w8hJccZnNS0rKaFoaIjZK1fS097O5glBRiJlZBNvIntPnOCumhpeHxjg52+8wfJQiFkDA7Tu3o2tWEH9ggUEtm2jwufjsooKhru7GRgdpW5ggIbiYrYcP07H+fMMV1TwG1deyeuvvkp5eTlvHjjATXPmUDZnDifOnOFbP/0p137rW+x65ZWoGbnIAGh8fOmC8nI6Z83iqYMHaT9zhsF9+/j9j370wo3o3+3bR0F5OR9YuJDQ0aPM6u9nf18fbadPM9rQcCEIT0VmJFVlq8lcTkETpEgspgockzk+OfLzeNuVV7LX7+eJPXuoHhxk4fLlWpNSRCTJFFiKxCFZ2ZiJwdZwONg6euIEvqVLJ80IRgYBZRs3YgMDHN+7l6WzZ7Nv/376AwHeHhzk3KlTuGPHuOXOO2mb8BzjAUpTUxO7tm1jFVDoHEOzZnGsoCBqueZkQc2gGY/+8z9zZPt2QmYsuf563vWBD1w0tf9Ljz5KqLGRuXV1uN5edr/0EtV+PwO9vdy4ahUHAgF2nz3L60ePcn1hIcNbt1IwMEC/309dUxP3LlpE/5w5PNfWxsGuLo4ODbHb5+P6+npGa2qobWhgYTBI/+AgFfPmAVBYVETjrFn0BwKsCJf/HnnxRXzOsWT1ajo6Ojj9yisXAqCnt27lUCBAhxk31dfznhtuYPPu3fR0d1NTXn5hNtjyUIianh6WLltGX2kpJ48fp7CzkyN+P78RETylorw0G8pWNUGKxGKy88mwz5fU8ckTP49XL13K/Joa9vj9rL377qQek4iIKLAUiVks2ZhYA8+JPfaDp09zZVUVe0+fvrBtuozgeKAxeuIElc7Run07Z4aGWFdVxUqfj5cPHeLQtm2cHhiI+rhVVVUsb2nhV7t382ZnJ8vWreOmO++M2t6JQc3Bo0d5dOtWzvb2suLcOX5t3jzO9ffz0s6dfO/ZZ1n1yU9SfPQoq6qqmLNwIb72dl598UXWXHUVi/x+3tyzh9DoKJ2Fhfh6erh5YICuUIgjra0sLyzkhsWLKa+s5KHXXqO3rIymBQtoHB6m9oorWNPWxruKirhm4UJeaGujp6+P7iNHOFdczHWFhZwvL+f1kRGuuvJK+ru7x97Lc+e4493vvhCQ/eiBB3h3c/OFG8758+YRfPVVzhQXc6a6mh1vv80rbW0sKCzkh488wqJlyxg0o7e3l9c7Oug7eZLasjJml5cza9Eirrn66kkzy1OVl8bTUZGKstVk0yyzMpXxz/tkZamhoqKkjk/W51FEJL0UWIrEaLpszHSB58TSy4NDQyxfuhSA0spKjgWDlFVVXXi9k/39tJ86xT98+cucaW+nvKGBmz/6UVaH10kEOFlUxDOHDjG8ezfB/n5WO8eJ8+dxp09zfs4cDvzqV+zZuZO+U6e47NZbOdvRwZEdO+ju6uLfz56lrLSUxuXLueOLX2R1S8uF8UiR4z4P79zJkR076Onq4qGeHk51dzPQ24vfjLUFBZQWF/PKnj3UFhezoLCQXadP88ODB/mta67hreJiRnw+2o8eZXVJCQfa2znR2cmrp0+zprCQ1u3bqXWO4+XlNDrHSqCutJThUIjr587lzWCQ7Xv2ML+6moH+fmqLi9nb18fh/n5mnzjBFaWluMFBzvh8tPb30/3WW5QWFbH4qqu49YorOD7JMgOXjYzw1oEDcOIEg6dPM+LzcfTUKbpCIbrffpuFIyMMnzrF4OgoC7u7KR0a4uz589R0d1NXUMC8kRECXV34CwsZKirCX1Nz0WcllvLSmZSNJrNsNRU0y6xMJpay1Nc2b6a6rAxIzvhkfR5FRNLLc4Glmd0J/E+gAPhn59w3MtyknJDJCTWS9drpOIb2QICtP/951NLO6Xq/xwMYO3eOrU8/TWdbG11nzvDjn/6UK669lsG332ZRaSnlxcWUDw/zozfe4GPvex+XL15M8bx5bD50iNtWrOBYVxc7du/mmT178HV2ckd1NRVmdOzZw3c3beKHLS00NDRQfPIkzQ0N1Jw7x9tnzrB0ZITrgUVmPDk8TLCvj+bSUipHRyl/8kn+/eGHGa2spKK/n4HeXgpmzcLX1ERxUREPf/GLPFxXx4KSEj58001cvngxv9y5k+/dfz+LneOm6mo4dYo3u7rYfuYM1efPc/rcOVwoxEHnmF1YyPzKSq4uLeVYby9HOzqoqKhg9Zo19A8NERgaosfv5/G2Nq4sLua6ujr6Ojs5PTrKucJC5g0NcWpkhJKiIvzAkd5eBoaHuWn+fP5u/34aenpwwPGDB+krLubygQGaQiGOHD/O22YUVlTwwcFBSktKqJ8zh1dOn+a1vj7e3dISdXyo3+9n96uvcse73kVZVRX9Q0PsHhnhVx0dfKqmhqKSEmYDC8+d46qqKl49doyVRUW4UIjuoiLes3gxlznH684xUFZGyLm4P2u5WDaaDeW6knrRztWxlKUme+ZWfR5FRNLLU4GlmRUA/wt4HxAAWs3sMefcG6l4vVQGKjN57ngfE+v+0Wbl/PETT1B92WUsXL6cgVmzeO1nP+PE3r0UlJZy1bp1fPCTnwSIa9bSaO0a9vkYOXGCdUuXXpKV6ejo4OUf//iibNz8+fOjHtP4MSwcHeXEgQO8tX8/vxgeZumHP8wHP/lJGhob2fz44/zsH/+R4KFD+GbPZuG7380N115LcSjEiWCQjsOH6XjrLQaCQarnzGHpqlUXZQC3t7by5N/+LcWHD3NFXR0L587l1Sef5F+efZbFq1YRPHmSvUNDVJeX89Krr3Lk4EGGz5/Hli9nXmsrrU8+ScfZswSPHqXp/HneXVZGb18fx7u66Ojt5XIzVpSVcVlzMyOlpWzu6+OH+/ZxXUkJw34/s264gUdbW+k7coRbli+n4vx5rhsYoPrECXwDAzQ4R7Nz7BweJlBXR+HJkxw8fZrKs2f5uHOMmvGmGQOhEL5QiGKgeGSE5T4fC4JBes6cYXtnJ79VXs61ZWUcGhriO6+9xvbjx/nyypU89fbbXFVSQndBAefOn+dXzzxDwfHjdA0P89O2NurOn2fOrFmUDwywYNYs3jNrFvOGhmgPhTgeCtExNMSyqiqKgGtmzaK9qwufz0eF389V9fUEhofx+/38h8WLOdjezkhvLwOFhVSZcSQUosnnw5wjODiIq6uja9Ys3jp5kqGmJnpvvZVf/ehH3FhUxK3NzfTs2UNoaIhgVxeHzPjzq6+mvqiIXx4/TkF9PVXDwxTOmzfpMgPBwUFGSko4B2DGOWCopIQ5paWsu+wyDnV301RUxFUFBZwoLOTQwAA3zJlDvc/H7Lo69syaRd/ZsxwIhbjn1lvZNTo63SnhErlYppcN5bqSWpHXm8GCAnY89xwPP/IIZwoLuWztWoj4O5z4eU92IKjPo4hIenkqsARuBA46594CMLMfAOuBpAeWqZy9cCbPHe9j4tk/sqe4vaeHrjfe4MOFhfQFg/S9/jrP/ehHtPj9/Md58+g7d44nHn2UHxw+TGNDw4WAcHyJiPfdcss7Y2KivN7Edr347LOcCgYZbmi4MAHKKuDRhx5iZMcOPlRTw8JFizgWDPLjr3+dshUr+Ni1115yTG2trSwcHaXzlVco6ujgN8rKCM6axU83b+bl8+cZaGhg7ze/yUfOnWN1VRVtZ87w0EMPceyNN3jXLbfQ99RTFA4OsuTsWd5XXMxwXx/Mns22b34TvvQl5s+fT+sDD3Brby/XzZ3LkBlbDh6E8+dZXVbG/GCQ0oYGfvLkkxQ4x7KBAT5WXMzpUIidb7/Ns1//Oivr62k8dYrRvj5OjI7SPjxMXVER15eW8lxXF011dVxZUUFXRwdLVqzgzsZGHhwa4rq77mL/pk3c3tBA+6lTVJSX81YoRODYMX5zYIDSvj6qCwoo8fmoO3eOnwUCLAoGua6oiDsrKng+GKTSOSpmzcI3OkrQOeqBo8CsggIWOEfhyAjnz56lxedjbihESWEhl/t8/LoZ/3L6NEsrKjjf3s71DQ0cOX2ap55+mqFAgMUDA7xvZIRiMw6GQrwyMMDp0VGuKSxkRVERO4eHKQWuGh3lhfPnKTh3jkKfj6tLS/nV6Cg9AwNUl5binzuXLa+9RvHs2RQ6R2hkhPKyMoLnz3Pq3Dm6Cgu5vaqKzadOsXh0lJLSUlxdHft8Pn7/S19idUsL/T091Bw+TPGsWZRXV1Nw9ixD9fXU9fdzeWUl/efOcfmyZSxetIiSoiKOhUJA9JvVg4ODvOeuu9hz6hT9p05RVlnJ0ssu40xREW86R3tfHyWlpVSUluIDlsydS+3cuRw/coS68nJuWrGCnoEBKoqLKS0poczvj/s8katlel4v181X6aoGGr/eDJ87x5utrdwyeza3z5/Pz958kx0vvIDv9tsvWSd3XCoCQX0eRUTSx5fpBkzQAByL+D4Q3pZ0kcHWhYCnqoq21taMPHe8j4ln//7u7gvjVtra2lg1ezaL58xhuK+PQ6+/zq+Z0QhUl5aypLKSD1ZX0//KK5T09Fx4/vElIoZPnJjy9Sa2q+T8eW4Or/M3rrqsjDeefZa7ampYWlNDYUEBS2tquOncOfr27Yt6TP3d3Qx2dDCrr4+rysupKilhUUUF82bNYm53N7/8p3/i18y4ubqaypISas34dHk5J99+mz1PPcV7qqvxnT3Lzc5xXXU1V/j9+Do6uKumhpd//GPaWlu5bGSEuoIC/MXF1BYXM3dggOrz55lTUMBwXx/Lly6led48Ok+cYH5hIZSXs/Saa6grKOCmc+eoKCnhcF8fpaEQ15eUsLO3lzec4/q5c/GNjuLOnaOkqIhz4cl0iszwOXfRezbc18fiOXNYNXs2p8+epX9ggLpwUFkIzCosxIVCXB0K4UZGmFdURFVBAeU+H4OhEObzsY+xnpgOYHlFBUtLSwmFQpwdHaWxoACA0VCIkdFR5vp8FIRC9J87R3lZGQNA4cgIx99+mwU+H9eEQtTMmkVNURGrCgtZEAoxFH6MH6gtKuLtwkJed459ZhTV1LCwpobB4mKWrF7NnuJiNp86xe6yMpb91m9Rec01/Kynh9OzZlG1bBlzysooKC5msKyM1woKOFxezs7aWjYDj5eXsyYcVAIsvuIK5q1cyZ7iYl6ZPZtnQiFqLruM4rIy9p08yaGhIcrnzOHw2bOUzp9/4Ya1obGRFevXs8fvZ3NnJ3v8fhruuIOlCxaw9uab+dCdd7L25pupqqmh8bLLOLVgAZXLlnHFokW8fP48+0ZHuX7FCp49dYpjs2czXF7O4d5eXuvvp2bePHYFgzRFjHuNVVNLC7uCQXr6+giFQvT09c34uUSmElEN9AHgKuATZnZVKl5r/Hozfq2p9fupKC1lQU0Nlc6xfdeuKT/vDY2NrL37bj50332svftuBYUiIlnEaxlLi7LtosFLZnYfcB/AokWLZvxCqSxDm8lzx/uYePaPzIz0nz5N9Zw59A8OUlpZyZm9e5k/axZnIkr5FpSVMXroEHbu3Duvd/o0C6uq2DHNrKUT21VaWUnR0BD9EY872d+PGxhgYcRENQA1Ph8jZ85EPaayujo6t26lcmSEsnCQfPLcOWoqKyk6d47hri7mz5lDSeHYR/rc+fOsKC1l6MwZhk+epG7hQkbOnaMm/POKkhKGw8d05u236e/uZmF9PSM9PQydO4e/uJii0VHOjYzgCgvxV1aOtaeoiMqyMm5cswafjX1cR0ZGqPH56HeOmquv5tC2bZQPD3PY5+PdtbWUFhZSXl1NeyhEV38/s/x++gYG2Nnby5Jbb73oPSutrKR/aIjq0lJm+f0809dHrc/HPOc4ATw1OkqxGdVmnCgs5O3BQRrKyzlx5gwnz5+nqLiYRQUFnAROhUL0FRbiM6NzdJQDhYUsASgqInj+PAOjo+wDSoqLOTw0xOorr+SVo0epLCyksKiI4vPnAej1+VhQWMh5YAQ4Ey65nT0ygq+wkGKfj8Ds2SxftIhZc+ZweGCArtmz+Y0bbmBlxIyPq9evB+DlBx/EvfUWO9vaGKmpodfvZ0FREV0FBXxg9WpW3nADxxcuvGQ5gKaWFvYfP86qlSu5bc0aDh49yjN79jBr3jy+/+abXF9Zia+mhtL58zlWUMCKiBvWaIumT8xidtXVUVpby9yKCvYXFPD8jh248nIar7iCvT4fbzY20nDttbztHK+fPEnFnDnMWriQFTMsoVeZnqRR2qqBxq8349cagP6hIWobGph32WV8r7X1knVyRUQkN3gtsAwACyO+bwSOR+7gnNsAbABYvXp1/DNmhKWyDG0mzx3vY+LZP7IU0F9eztGTJ+kLhWhsbqb8wAE6enuZXVp6Yf/j/f0UlJfjioreeb3wrKWlE2Ytnfh6E9u1oKmJnVu2MDpnDqFQ6EKQMb+5mWPBIEsjZtPsDYUoLC+/6PnGX6OppYXHnniCpaOj9A0PM2LGruFhGubOpauoiOK5c+kYGGDJyAj+WbMomjWL/WfPUuL3U1xdTffZsxQWFdEbnmSlb2iI4ooKjgWDlDc0UFZXR+nQEN2BAIfa21laXk4wFCJw7hzzKypY0NQEwLmiIgrLyzk5OEhtuPyxsLCQ3uFhqioruaypif2nTzN6+DDNc+dyamSEV0+epHbpUk6fP8/TAwPMr63l2OgoPcuXc3N4ncnx92xBUxOHt22jYnCQ5uXLOX7mDI+cPUuRGSWFhQybMc/n40xBAauWLOGJ48e5sqqKAefYMzJCuxllBQUsnD2bmxYu5PuBAGXBIEMlJSxbu5bnW1spARqKi+nw+Xisv5/Z8+ZR3NjIvIoK9vX18XpREedmz+bEm2+yfN48Rs6fZ/fAAF1FRRyfPZuzw8M8DZwsK6O2vJy+2bPpra2laO5cdsyezZLVq7lh1Sp6A4GoN5A333svba2tHHvzTQJHj3I2EKBl2TLWXn01pSUlY2OqomTtLgnEFi7kwx/5yEUz7h4Lj82dLtiLFtTdfO+9wFjWvbSkhCuuvRafGcWhEGV1ddyWgsmiVKYnaRKtGmhNKl7owjJIs2bRc/YsJT4fh8+epbG5GVdUxA3vf7/WkBQRyVFeCyxbgSYzWwq0A/cAn0zFC6VytriZPHe8j4ln/8ib6GNVVRw4eZJ1zc3Mqa7msmuu4fEDB2gBTg4Ojo2xPHWKshtvZKi2lp6+PqrLyqgJz1r6vhUrLgoQJ77exHa5oiJ6li/HzZ17UZBReeONbP7mN7kLWFhVxbFgkK1FRVSsWHHhNSNfo6GxkZbPfpaf/+3fsv3wYZpra1ne2MiJ4WF6Fi/mPf/xP/LkN78JJ0+yurKSHud46MwZFl93Hc233MIvn3qK0OzZvHz2LGUnTzIcCsGSJWzu7eXmT3+a+fPns//4cRa2tNBVWcmjb7zBjlmzKFu0iPobb2ROdTU9fX101dVRMWsWv+ru5uZQiCIzzhYVsdc51s+bx5zqaupaWnhsdJQ+n4+9ZhSUlbFo8WL8NTX4zBgMhSioq+PmiEDlwntWXU3vVVfx8z17qLjuOnpDIU63tbGiqIhZoRCdg4N0OMcrfj/Lq6pYu2IFuzs7+VV7O7OvvZZrw7PFXt/YSMngIJWHDvH6qVM0rl5N84038vYnP8n2TZt47vhxZi9YwPvWr6f0/Hl+umMHvoEBltx5J/feeScAj//937PrxRe5vKeHufPnU+r3c3lDA/OXLKGgro7Thw/T7RxLVq/m7mjrX05S0hkte9jW2squ7m7K/P4psxiTBWIzCdCmei6RHDNtNRAkpyJo/Hrz0hNP8PCzz3L93Llc2dKCKyrSjKwiIjnO3AymyU8lM7sL+DvGJhj4jnPu/sn2Xb16tdu+ffuMXytfZoWd7nGpnBV2qv22t7bGPCts5HO/9MQTHNmxA184qLkpHNQkOitstHZPfA/Gt01sw+JVq2J6b+L9XbYHAvzs4Yd545lnGB0cZN7Klbz/934P4JL3bqrjmOnyLi898QR7n3+evq4uyubO5erbbruw/IqI15nZq8651ZluR6aZ2c3A15xz7w9//1UA59z/f7LHJHp9hcwucyUiIqkx1bXVc4FlPJJx4RMRkdykwHKMmRUCbwLrGKsGagU+6ZzbO9ljdH0VEZFoprq2eq0UVkRERJLIOTdiZn8E/IJ3qoEmDSpFRERmQoGliIhIjnPObQY2Z7odIiKSu7y2jqWIiIiIiIhkGQWWIiIiIiIikhAFliIiIiIiIpIQBZYiIiIiIiKSEAWWIiIiIiIikhAFliIiIiIiIpIQBZYiIiIiIiKSEAWWIiIiIiIikhAFliIiIiIiIpIQc85lug0zZmbdwNEkPFUt0JOE5/GKXDqeXDoWyK3jyaVjgdw6nlw6Fpj58Sx2ztUluzH5IAnX12z8DGZbm7OtvaA2p4vanB7Z1uZktXfSa2tWB5bJYmbbnXOrM92OZMml48mlY4HcOp5cOhbIrePJpWOB3DuefJCNv7Nsa3O2tRfU5nRRm9Mj29qcjvaqFFZEREREREQSosBSREREREREEqLAcsyGTDcgyXLpeHLpWCC3jieXjgVy63hy6Vgg944nH2Tj7yzb2pxt7QW1OV3U5vTItjanvL0aYykiIiIiIiIJUcZSREREREREEpJzgaWZ3WlmB8zsoJl9JcrPzcz+Pvzz183s+oifHTGz3Wa208y2R2yvNrOnzKwt/P8crx+PmV0RPo7xf31m9oXwz75mZu0RP7vLI8eywsxeNrNhM/tSLI/1+O8m6vGY2UIze9bM9pnZXjP7k4ifZeR3k8jxhH/mqb+dBH43nvu7ifF4fjv89/+6mb1kZqume6yHfzdRj8Wrfzdysel+v14w2WfJ65+jaOfZTF4DpzPZ+dRr77OZfcfMusxsT8S2Sd9XM/tq+PN9wMze75H2/rWZ7Q+fNzeaWVV4+xIzG4x4r/93uts7RZsn/Rxk+j2eos0/jGjvETPbGd7ulfd5snNb+j7Pzrmc+QcUAIeAZUARsAu4asI+dwFPAAbcBGyL+NkRoDbK8/4V8JXw118B/kc2HM+E5znB2LozAF8DvuTB381coAW4P7J9Uz3W47+byY5nPnB9+Oty4M2I40n77ybR4wn/zDN/O4key4TnyejfTRzH8y5gTvjrD4yfB7z2t5PgsXju70b/4v/9euHfZJ8lr3+Oop1nM/F3nMBn4wSw2GvvM3AbcD2wZ7r3Nfw52QUUA0vDn/cCD7T314DC8Nf/I6K9SyL389h7HPVz4IX3eLI2T/j53wD/xWPv82TntrR9nnMtY3kjcNA595Zz7hzwA2D9hH3WA99zY7YCVWY2f5rnXQ88GP76QeA3ktjmqSTreNYBh5xziSx2nahpj8U51+WcawXOx/FYz/5uJjse51yHc25H+OszwD6gIT3NnlQiv5+pZOL3k6xj8cLfDcR2PC85506Fv90KNMbwWK/+bqIei0f/buRisVyzMi7HPkuZugbGyyvn00s4554HTk7YPNn7uh74gXNu2Dl3GDjI2Oc+baK11zn3pHNuJPxt5DXAEyZ5jyeT8fcYpm6zmRnwceD7aW3UNKY4t6Xt85xrgWUDcCzi+wCXXiym2scBT5rZq2Z2X8Q+9c65Dhj7pTGW7UiHRI9n3D1c+uH/o3DJxHfSVDoTSztn8lgv/26mZWZLgOuAbRGb0/27gcSPx0t/O0n53eCNvxuI/3g+w1gVw3SPzYbfTeSxXOChvxu5WLL+9tImymfJy5+jaOfZTF0D4zXxfOrl9xkmf1+z4TP+e1x83lxqZq+Z2XNmdmumGjWJaJ+DbHiPbwU6nXNtEds89T5POLel7fOca4GlRdk2cdrbqfa5xTl3PWPlV58zs9uS2bgZSPR4MLMi4MPAjyJ+/k/AZcC1QAdj6fxUi+VYUvHYVEm4TWZWBjwKfME51xfenInfDSR+PF7620nG78YrfzcQx/GY2R2MBWNfjvexaZLIsYxv99LfjVzMa5+3KUX5LHn9c+Sl82zMopxPvf4+T8XTn3Ez+wtgBHgovKkDWOScuw74IvCwmVVkqn0TTPY58PR7HPYJLu4o8dT7PMl1MuquUbYl9F7nWmAZABZGfN8IHI91H+fc+P9dwEbeSQd3jpeXhv/vSnrLo0voeMI+AOxwznWOb3DOdTrnRp1zIeAB0lNiEMuxzOSxXv7dTMrMZjH2R/+Qc+7fx7dn6HcDCR6Px/52EjqWMK/83UCMx2Nm1wD/DKx3zvXG8FjP/m4mORYv/t3IxZLxt5cW0T5LXv8cTXKezdQ1MB4XnU+9/j6HTfa+evYzbmb3Ah8Cftu5sQF04RLH3vDXrzI2hm555lr5jik+B559jwHMrBD4CPDD8W1eep8nuU6m7fOca4FlK9BkZkvDPWT3AI9N2Ocx4FM25ibgtHOuw8xmm1k5gJnNZmwg9J6Ix9wb/vpeYFOqDyRsxscT8fOJvSrjH6pxd/POcaZSLMcyk8d6+XcTVbg2/9vAPufctyb8LBO/G0jseLz2t5PIZ22cV/5uIIbjMbNFwL8Dv+OcezPGx3rydzPZsXj070Yuloy/vZSb7LPk5c/RFOfZTF0D43HR+dTL73OEyd7Xx4B7zKzYzJYCTcArGWjfRczsTsaqOz7snBuI2F5nZgXhr5cx1t63MtPKi03xOfDkexzhvcB+51xgfINX3ucprpPp+zy7DM9glOx/jM2S+iZjvQV/Ed72h8Afhr824H+Ff74bWB3evoyxmZF2AXvHHxv+WQ3wDNAW/r/a68cT/pkf6AUqJzznv4b3fT38oZrvkWOZx1jvSR8QDH9dMdljs+B3E/V4gHczVmrwOrAz/O+uTP5uEjwez/3tJPhZ89TfTYzH88/AqYjP0/apHuvx303UY/Hq343+Tf/79dq/yT5LXv4cTXaezdTfcRztvuR86rX3mbGgt4OxydwCjJXgT/q+An8R/nwfAD7gkfYeZGys3Pjn+X+H9/3N8OdlF7AD+HUPvceTfg4y/R5P1ubw9u+OX68i9vXK+zzZuS1tn2cLP6mIiIiIiIjIjORaKayIiIiIiIikmQJLERERERERSYgCSxEREREREUmIAksRERERERFJiAJLERERERERSYgCS5EpmNmome00sz1m9iMz86fodVab2d+Hv15rZu+awXN8wcw+Ff56Rbjdr5nZZQm27Vozuyvi+w+b2Vdm+Fx1ZvbzRNojIiK5L1uvv0lozxIzu2RtTzO72sy+m4zXEEkVBZYiUxt0zl3rnGsGzjG23t60zKwwnhdxzm13zn0+/O1aIK4LW/j1fg94OLzpN4BNzrnrnHOHIvYzM4v37/5axtZBGm/rY865b8T5HOOP7QY6zOyWmTxeRETyRrZef2fclqk453YDjWa2KFnPKZJsCixFYvcCcLmZVZvZT8zsdTPbambXAJjZ18xsg5k9CXzPzBab2TPh/Z4ZvxiY2cfCPbC7zOz58La1Zva4mS1h7OL5n8I9tbea2WEzmxXer8LMjox/H+E9wA7n3Eg4u/gF4PfN7Nlw7+c+M/t/GVu4d6GZ/ZOZbTezvWb2X8efxMxazOylcNteMbNK4OvAb4Xb81tm9rtm9o/h/Sc7xu+a2d+Hn+stM/toRFt/Avx2En8vIiKS27Li+hveb4uZ/T9m9hzwJ2b262a2LVxB9LSZ1Ue0+Tvh/d8ys89PeF7MbFn4cS3hTT8F7knuWyuSPAosRWIQ7nX8ALAb+K/Aa865a4A/B74XsesNwHrn3CeBfwS+F97vIeDvw/v8F+D9zrlVwIcjX8c5dwT438DfhntqXwC2AB8M73IP8Khz7vyEJt4CvBp+js0Rz3FH+OdXhNtynXPuKPAXzrnVwDXA7WZ2jZkVAT8E/iTctvcCZ8Pt/WG4PT+c8LqTHSPAfODdwIeAyAznduBWREREppFN198IVc65251zfwO8CNzknLsO+AHwZxH7rQDeD9wI/GVk0GpmVwCPAp92zrWGN+v6KZ6mwFJkaqVmtpOxk/nbwLcZC5b+FcA590ugJpzZA3jMOTcY/vpm3imN+dfw4wB+BXzXzD4LFMTQhn8GPh3++tPAv0TZZz7QPcVzHHXObY34/uNmtgN4DVgJXMVY8NkxfgFzzvWN98BOYbJjBPiJcy7knHsDqI/Y3gUsmOZ5RUQkv2Xz9TeyE7YR+IWZ7Qb+L8auueN+5pwbds71MHZtHL9W1gGbgP/gnNsZsb+un+JpSav9FslRg865ayM3mJlF2c+F/z87xXM5AOfcH5rZGsZ6QXea2bVTPAbn3K/C5ay3AwXOuUsG9QODQMkUT3OhXWa2FPgS0OKcO2VjkwGUABZxHDMV+fjhiK8j37OScHtFREQmk83X38i2/APwLefcY2a2FvhaxM8ir5OjvHNffho4xlg2dG/EPrp+iqcpYykSv+cJjxEMXyR6nHN9UfZ7iXfGQvw2Y+UwmNllzrltzrn/AvQACyc87gxQPmHb94DvE723FGAfcHmM7a9g7KJ3OjzW4wPh7fuBBeNjOcysPFyCFK0946Ie4zSWA9EuziIiIlPJxutvJdAe/vreKfaLdI6xSfg+ZWafjNiu66d4mgJLkfh9DVhtZq8zNnZwsgvF54FPh/f7HeBPwtv/2sx229h04s8DuyY87qfA3eOTB4S3PQTMYeziFs0TwG2xNN45t4uxEti9wHcYKw3COXcO+C3gH8xsF/AUY72jzwJXjU/eE+MxTuUO4GextFVERCTC18i+6+/XgB+Z2QuMBbMxcc6dZWyOgv9kZuvDm3X9FE8z5xKtfBORVAvPqrreOfc7U+yzEfgz51xb+loWv/BMfOudc6cy3RYREZGpeOX6a2bFwHPAu2OY/0AkIxRYinicmf0DY+Wqdznn3pxivyuAeufc82lrXJzMrA64xTn3k0y3RUREZCpeuv6aWRPQ4JzbkqrXEEmUAksRERERERFJiMZYioiIiIiISEIUWIqIiIiIiEhCFFiKiIiIiIhIQhRYioiIiIiISEIUWIqIiIiIiEhCFFiKiIiIiIhIQv4/IZtQzAUy51YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df['Rpor'] = df['porosity'].rank()\n",
    "df['Rperm'] = df['permeability'].rank()\n",
    "# df.head()\n",
    "\n",
    "plt.subplot(121)\n",
    "plt.scatter(df['porosity'],df['permeability'],c='red',alpha=0.3,edgecolor='black')\n",
    "plt.xlabel('Porosity (fraction)'); plt.ylabel('Permeability (mD)')\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.scatter(df['Rpor'],df['Rperm'],c='red',alpha=0.3,edgecolor='black')\n",
    "plt.xlabel('Porosity (rank)'); plt.ylabel('Permeability (rank)')\n",
    "\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.2, hspace=0.3); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's extract two 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": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "por = df['porosity'].values\n",
    "perm = df['permeability'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a quick look at the data. It is always a good idea to visualize data before we calculate statistics.\n",
    "\n",
    "* check for linearity with respect to the related assumptions for each measure\n",
    "\n",
    "* chech for outliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEWCAYAAACNJFuYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA730lEQVR4nO3deZwcZ3no+9/T+0xPzyJptFqbbclajGzLkneQhA02xMSQQHAg4CQQX07IhQRIAsm558I98T0knHBOgIQbJyxmFY4h2IDt2BGSvEhepJEsW/s+WkbSSLP19HT3dHU/94+ukVuz9PQsvczo+X4+/enu6lqeqqmpp+p933pLVBVjjDEml6fcARhjjKk8lhyMMcYMYMnBGGPMAJYcjDHGDGDJwRhjzACWHIwxxgxgycFcNkTk90XkhTy/PyUiDww2roh0i8iVpYhzIhKRD4vIM+WOw4wfSw5mVETkmIjE3YPmWRH5jojUlDuusVDVd6nqI0P8VqOqRwBE5Lsi8jeljS4/EdkkIgn373FeRH4mIrNKtXxV/aGqvjMnHhWRq0u1fDP+LDmYsXiPqtYAK4HVwH8dycSSZfvg+PkT9++xGKgH/tdIJra/h8llO4IZM1U9BTwFXAsgIreIyBYR6RCR10Rkbd+47hnuQyLyItADXOmeZf6xiBwUkaiI/HcRuUpEtopIl4g8KiKBnHncKyI73flvEZEVOb99XkQOu/PZIyLv6xeuiMjXRaRTRPaJyJ39Yvv4YOvYdyYsIg8CHwb+wj1L/4WI/LmI/LTf+F8Xkf89yHw+LyKP9Rv2DyLyNffz74vIETf+oyLy4XzbfjCq2gb8lDf/HreJyKvuOr8qIrf1W+f+f4984w8aX24xnIg8547+mruNPigib4jIe3Lm43evcK4f6fqZElFVe9lrxC/gGHCX+3kusBv478Ac4ALwbrInH+9wvze6424CmoHlgA/wAwo8AdS6w5PABuBKoA7YAzzgTr8SOAfcDHiBB9xYgu7vHwBmu8v+IBADZrm//T7gAH/mLveDQCcwJSe2j+eM+0LO+ipwtfv5u8Df5Pw2y11Ovfvd58Z44yDbbT7Zg3Ct+90LtAC3AGGgC7gmZ77LC/x75MY+Dfg18H1gCtAOfMSN63fd71OH+HvMGGr8fPHl217u978AfpLz/T7g9XLvx/Ya+mVXDmYsfi4iHcALwGbg/wV+D3hSVZ9U1YyqPgtsI5ss+nxXVXerqqOqKXfY36pql6ruBt4AnlHVI6raSfaq5AZ3vD8C/llVX1bVtGbrCJJkD66o6r+p6ml32T8BDgI35Sz7HPC/VTXl/r4f+I2xbARVbQGeI5uYAO4Bzqvq9kHGPQ40Ae91B70d6FHVl9zvGeBaEalS1RZ3exTqa+7f4zWyCeczZNftoKp+393ePwb2Ae/Jme7i3wN45zDjjza+HwDvFpFa9/tHyCYvU6EsOZixeK+q1qvqfFX9Y1WNkz0z/oBb5NPhHqzuIHuW2efEIPM6m/M5Psj3vsru+cBn+81/LtmrBUTkozlFTh1ki1am5czrlKrm9jZ5vG/aMXqEbGLEfc934PsR2TNygA+531HVGNmrmU8ALSLyKxFZMoIYPuX+Peao6odVtZXsuh3vN95xsld4fXL/HkOOP5b4VPU08CLw2yJSD7wL+GGB62XKwJKDGW8ngO+7B6m+V1hVv5wzzli6Aj4BPNRv/tWq+mMRmQ/8C/AnZItN6slehUjO9HNEJPf7POD0CGMYLP6fAytE5FrgXvIf+P4NWCsiVwDvw00OAKr6H6r6DrLJdJ+7PmNxmmxCzTUPOJXzPXd98o4/xvj6EugHgK2arasyFcqSgxlvPwDeIyJ3i4hXREIi0ncgHA//AnxCRG52W9eEReQ3RCRCtkxcgVYAEfkD3ErZHNOBT7kVoh8AlgJPjjCGs2TrQy5S1QTwGNkD/Suq2jzUxO4Z/SbgO8BRVd3rxjtDRH5TRMJki8q6gfQIY+vvSWCxiHxIRHwi8kFgGfDLkY4/wvgGbCOyCXQl8Gnge2NaK1N0lhzMuFLVE2QrG/+K7EH6BPDnjNO+pqrbyNY7fINsRekhspWhqOoe4O+BrWQPTm8hW5SR62VgEXAeeAh4v6peGGEY3wKWuUVXP88Z/oi7zELK0n8E3EXOVQPZbfRZsmfvbcAa4I8BROStItI9wjhx1+1ed74XyFYM36uq50cx/pDxDeKLwCPuNvodd95xsq2oFgI/G+m6mNKSS4tfjTGjJSLzyBa1zFTVrnLHU4lE5L8Bi1X194Yd2ZSVr9wBGDMZSPbmsc8A6y0xDE5EpgAfI9tSyVQ4K1YyZozcMvgusvd0/N9lDqciicgfkS1ifEpVnxtufFN+VqxkjDFmALtyMMYYM8CErnOYNm2aLliwYNTTx2IxwuHw+AVURBMpVrB4i83iLa7JHu/27dvPq2pj3pHK3X/HWF433nijjsXGjRvHNH0pTaRYVS3eYrN4i2uyxwtsU+tbyRhjzEhZcjDGGDOAJQdjjDEDWHIwxhgzgCUHY4wxA1hyMMaYMnEch1gshuM45Q5lgAl9n4MxxkxU7W1tNO/YQSCVotfvZ97KlTQ0NJQ7rIvsysEYY0rMcRyad+xgcSjE0sZGFodCNDc1VdQVhCUHY4wpsWQySSCVoioUAqAqFCKQSpFMJssc2ZssORhjTIkFg0F6/X7iiQQA8USCXr+fYDBY5sjeZHUOxhhTYj6fj3krV3KgqYlANHqxzsHnq5xDcuVEYowxl5GGhgYia9aQTCYJBoMVlRjAkoMxxpSNz+eruKTQx+ocjDHGDGDJwRhjzACWHIwxxgxgycEYY8wAlhyMMcYMYMnBGGPMAJYcjDHGDGDJwRhjzACWHIwxxgxgycEYY8wAlhyMMcYMYMnBGGPMAJYcjDHGDGDJwRhjzACWHIwxxgxgycEYc9lzHIdYLIbjOOUOpWJU5lMmjDGmRNrb2mjesYNAKnXxcZ3GkoMx5jLmOA7NO3awOBSiqr6eeCLBgaYm8HrLHVrZWbGSMeaylUwmCaRSVIVCAFSFQgRSKTKZTJkjK7+iJwcR8YrIDhH5pft9iog8KyIH3feGnHG/ICKHRGS/iNxd7NiMMZe3YDBIr99PPJEAIJ5I0Ov34/HYeXMptsCngb053z8PbFDVRcAG9zsisgy4H1gO3AP8k4jYtZ0xpmh8Ph/zVq7kQCLB3tZWDiQSVufgKmqdg4hcAfwG8BDwGXfwfcBa9/MjwCbgL93h61U1CRwVkUPATcDWYsZojLm8NTQ0EFmzhmQySTAYxOezqlgAUdXizVzkMeB/ABHgc6p6r4h0qGp9zjjtqtogIt8AXlLVH7jDvwU8paqP9Zvng8CDADNmzLhx/fr1o46vu7ubmpqaUU9fShMpVrB4i83iLa7JHu+6deu2q+qqfOMULUWKyL3AOVXdLiJrC5lkkGEDMpeqPgw8DLBq1Spdu7aQWQ9u06ZNjGX6UppIsYLFW2wWb3FZvMUtVrod+E0ReTcQAmpF5AfAWRGZpaotIjILOOeOfxKYmzP9FcDpIsZnjDFmCEWrkFbVL6jqFaq6gGxF869V9feAJ4AH3NEeAB53Pz8B3C8iQRFZCCwCXilWfMYYY4ZWjpqXLwOPisjHgGbgAwCqultEHgX2AA7wSVVNlyE+Y4y57JUkOajqJrKtklDVC8CdQ4z3ENmWTcYYY8rI7vQwxhgzgCUHY4wxA1hyMMYYM4AlB2OMMQNYcjDGGDOAJQdjjDEDWHIwxhgzgCUHY4wxA1hyMMaYMnIch1gshuM45Q7lEtZxuTHGlEl7WxvNO3YQSKXo9fuZt3IlDQ0Nw09YAnblYIwxZeA4Ds07drA4FGJpYyOLQyGam5oq5grCkoMxxpRBMpkkkEpRFQoBUBUKEUilSCaTZY4sy5KDMcaUQTAYpNfvJ55IABBPJOj1+wkGg2WOLMvqHIwxpgx8Ph/zVq7kQFMTgWj0Yp1DpTzDujKiMMaYy1BDQwORNWtIJpMEg8GKSQxgycEYY8rK5/NVVFLoY3UOxhhjBrDkYIwxZgBLDsYYYwaw5GCMMWYASw7GGGMGsORgjDFmAEsOxhhjBsjbuFZEQsC9wFuB2UAceAP4laruLn54xhhjymHI5CAiXwTeA2wCXgbOASFgMfBlN3F8VlV3FT9MY4wxpZTvyuFVVf3iEL99VUSmA/PGPyRjjDHlNmRyUNVf5ZtQVc+RvZowxhgzyeStkBaRB0SkSURi7mubiHy0VMEZY4wpj3x1Dh8F/hT4DNAECLAS+IqIoKrfK0mExhhjSi7flcMfA+9T1Y2q2qmqHar6a+C33d+MMcZMUvmSQ62qHus/0B1WW6yAjDHGlF++5BAf5W/GGGMmuHxNWZeKyGD3MAhwZZHiMcYYUwHyJoeSRWGMMaaiDFmspKrHc19AOxDNeeUlIiEReUVEXhOR3SLyJXf4FBF5VkQOuu8NOdN8QUQOich+Ebl77KtnjDFmNIbteE9E/g8ROQvsAra7r20FzDsJvF1VrwOuB+4RkVuAzwMbVHURsMH9jogsA+4HlgP3AP8kIt4Rr5ExxpgxK6RX1s8By1V1gaoudF/D1jloVrf71e++FLgPeMQd/gjwXvfzfcB6VU2q6lHgEHBT4atijDFmvIiq5h9B5Gngt1S1Z8Qzz575bweuBv5RVf9SRDpUtT5nnHZVbRCRbwAvqeoP3OHfAp5S1cf6zfNB4EGAGTNm3Lh+/fqRhnVRd3c3NTU1o56+lCZSrGDxFpvFW1yTPd5169ZtV9VV+cbJ22W36wvAFhF5mWxREQCq+qnhJlTVNHC9iNQD/y4i1+YZXQabxSDzfBh4GGDVqlW6du3a4cIY0qZNmxjL9KU0kWIFi7fYLN7isngLSw7/DPwaeB3IjGYhqtohIpvI1iWcFZFZqtoiIrN4s/O+k8DcnMmuAE6PZnnGGGPGppDk4KjqZ0Y6YxFpBFJuYqgC7gL+FngCeAD4svv+uDvJE8CPROSrZB8stAh4ZaTLNcYYM3aFJIeNbjn/L7i0WKltmOlmAY+49Q4e4FFV/aWIbAUeFZGPAc3AB9z57RaRR4E9gAN80i2WMsYYU2KFJIcPue9fyBmmDHOXtPuEuBsGGX4BuHOIaR4CHiogJmOMMUU0bHJQ1YWlCMQYY0zlGPI+BxG5I9+EIlI7TOsjY4wxE1S+K4ffFpG/A54me69CKxAie8/COmA+8NmiR2iMMZOY4zgkk0mCwSA+XyEl/aWR7xnSf+b2e/R+spXGs8h21b0X+GdVfaE0IRpjzOTU3tZG844dBFIpev1+5q1cSUNDw/ATlkDeNKWq7cC/uC9jjDHjxHEcmnfsYHEoRFV9PfFEggNNTUTWrKmIK4hC+lYyxhgzzpLJJIFUiqpQCICqUIhAKkUymRxmytKw5GCMMWUQDAbp9fuJJxIAxBMJev1+gsFgmSPLKv+1izHGXIZ8Ph/zVq7kQFMTgWj0Yp1DJRQpQQHJQUS2Ad8BfuTWQRhjjBkHDQ0NRNasqcjWSoUUK91Ptq+jV0VkvYjcLSKD9aBqjDGTiuM4xGIxHMcp2jJ8Ph/hcLiiEgMUdof0IeCvReT/Au4Fvg1kROTbwD8U0MeSMcZMOJXczLQUCqqQFpEVwN8DXwF+Svbehy6yXXkbY8yk09fMdGljI4tDIZqbmop6BVFpCqlz2A50AN8CPq+qfe2sXhaR24sYmzHGlEUmk8k2M62vB9xmptEoyWSy4op/iqWQtfyAqh7JHSAiC1X1qKr+VpHiMsaYsvF4PBebmVaFQhXXzLQUCilWeqzAYcYYM2nMW7mSA4kEe1tbOZBIVFQz01IYck1FZAmwHKgTkdwrhFqyHfAZY8ykVcnNTEsh39peQ7Z1Uj3wnpzhUeCPihiTMcZUBJ/Pd9klhT75emV9HHhcRG5V1a0ljMkYY0yZ5StW+gtV/TvgQyLyu/1/V9VPFTUyY4wxZZPvemmv+76tFIEYY4ypHPmKlX7hvj9SunCMMcZUgnzFSr8AdKjfVfU3ixKRMcaYsstXrPQ/SxaFMcaYipKvWGlzKQMxxhhTOfIVKz2qqr8jIq9zafGSAKqqK4oenTHGmLLIV6z0aff93lIEYowxpnIM2beSqra478eBJHAdsAJIusOMMcZMUsN2vCciHwdeAX6L7HMcXhKRPyx2YMYYY8qnkE5D/hy4QVUvAIjIVGAL2SfCGWOMmYQK6bL7JNnO9vpEgRPFCccYY0wlyNda6TPux1Nkn/r2ONlWS/eRLWYyxhgzSeUrVoq474fdV5/HixeOMcaYSpDvJrgvlTIQY4wxlaOQ1kqNIvIVEXlSRH7d9ypgurkislFE9orIbhH5tDt8iog8KyIH3feGnGm+ICKHRGS/iNw9tlUzxpjK5TgOsVgMx3HKHcqgCqmQ/iGwD1gIfAk4BrxawHQO8FlVXQrcAnxSRJYBnwc2qOoiYIP7Hfe3+8k+mvQe4J9ExDuitTHGmAKV8+Dc3tbG7s2baX7+eXZv3kx7e3vJYxhOIclhqqp+C0ip6mZV/UOyB/u8VLVFVZvcz1Gyz4eYQ7ZCu68b8EeA97qf7wPWq2pSVY8Ch4CbRrIyxhhTiHIenB3HoXnHDhaHQixtbGRxKERzU1PFXUEUkhxS7nuLiPyGiNwAXDGShYjIAuAG4GVgRs7d1y3AdHe0OVzaRPakO8wYY8ZNuQ/OyWSSQCpFVSgEQFUoRCCVIplMlmT5hRLVIR/ZkB1B5F7geWAu8HWgFviSqj5R0AJEaoDNwEOq+jMR6VDV+pzf21W1QUT+Ediqqj9wh38LeFJVf9pvfg8CDwLMmDHjxvXr1xe2poPo7u6mpqZm1NOX0kSKFSzeYrN4Ry+TydDb3U3I92Z7nITjEKipwePJni8XO954NErI40FEUFUSmQxVkcjwEw5hpPGuW7duu6quyjfOsHdIq+ov3Y+dwLqClw6IiB/4KfBDVf2ZO/isiMxS1RYRmQWcc4efJJuA+lwBnB4knoeBhwFWrVqla9euHUlIl9i0aRNjmb6UJlKsYPEWm8U7eo7jsHvzZhaHQlSFQsQTCQ4kEix/29vwuQmj2PG2t7fT3NREIJWi1+9n3sqVNDQ0DD/hEIoRbyGtlRaLyAYRecP9vkJE/msB0wnwLWCvqn4156cngAfczw/w5n0TTwD3i0hQRBYCi7Cb7Ywx48zn8zFv5UoOJBLsbW3lQCLBvJUrLyaGUmhoaGD5mjXMe+tbWb5mzZgSQ7EUsjX+hWz/Sv8MoKq7RORHwN8MM93twEeA10Vkpzvsr4AvA4+KyMeAZuAD7nx3i8ijwB6yLZ0+qarpka2OMcYMr6GhgciaNSSTSYLBYEkTQx+fz1eW5RaqkMiqVfWV7IXARcPW3KjqC2QfDDSYO4eY5iHgoQJiMsaYMan0g3O5FdJa6byIXIX7NDgReT/QUtSojDHGlFUhafOTZCuAl4jIKeAo8HtFjcoYY0xZFdJa6Qhwl4iEAY97Q5sxxpgycRyn6PUlw85VROqBjwILAF9f3YOqfqooERljjBlSe1sbzTt2XNIMthgKSTlPAi8BrwOZokRhjDFmWLl3d1fV12fv0WhqAu/4d0NXSHIIqepnhh/NGGNMMV3seqO+HnC73ohGychQDUNHr5DWSt8XkT8SkVlud9tTRGTKuEdijDEmr2AwSK/fTzyRACCeSNDr91/s9mM8FTLHXuArwFZgu/vaNu6RGGOMyWuou7uLsqwCxvkMcLWqni9KBMYYYwpWqru7C5nrbqCnKEs3xpgyK0Wz0PFWiru7C5l7GtgpIhuBix2OW1NWY8xEN1iz0ErsBK8cCkkOP3dfxhgzaQzVLDSyZk25Q6sIeZOD+wznj6jqXSWKxxhjSmKoZqGV9kS2csnbWsntMrtHROpKFI8xxpTEUM1Cg8FgmSOrDIUUKyXIPpPhWSDWN9DqHIwxE9nFZqFNTQSi0Yt1DhOlUrrYCtkKv3JfxhhTkInSAqgSHvpTqQrplfUREakC5qnq/hLEZIyZwErVMdx4sYf+DK6QZ0i/B9gJPO1+v15EnihyXMaYCSi3BdDSxkYWh0I0NzWVOywzCoV0n/FF4CagA0BVdwILixaRMWbCutgCKBQC3BZAqRSZTPE7dHYch1gshuMM+xRjU4BCrqUcVe3s9wxpLVI8xpgJLLcFUFUoVNSO4XLZzWzjr5C/2Bsi8iHAKyKLROTrwJYix2WMmYBK2TFcn6GKsuwKYmwKuXL4P4G/Jtt1xo+A/wD+pphBGWMKV2ktg0rdAijfzWyVsD0mqiG3nIiEgE8AV5N9Ctytqmqp2JgKUqnFKaVsATRUUZbdzDY2+YqVHgFWkU0M7wL+Z0kiMsYUxIpTsoYqyrKrhrHJt/WWqepbAETkW8ArpQnJGFMIK05501iLsiqtaK4S5NsKqb4PqupIEZ5RaowZPStOudRoi7Im2k17pZJvS14nIl3uZwGq3O8CqKrWFj06Y8yQrG+gsRuq22683nKHVnZD7kWqalvHmApnfQMNrtBioqGK5jJWUlJQU1ZjTAWzvoEuNZIWXF6vl07HIdrdTaSmpmQ37U0EtgWMMZPGSFpwtbe1sf/FF/EkEjzzyiu8dORISW7amyjsdMMYM2kU2oLrkrqGhQtZ2tjI7u5urrn9dkJuv1CXO7tyMMZMGoU+3a1/B4GRmhrqvF7S6XTJY65UlhyMMZNGoTfE2SNCh2fFSsaYSaWQFlzWDHh4RdsSIvJt4F7gnKpe6w6bAvwEWAAcA35HVdvd374AfAxIA59S1f8oVmzGmMmtkBZc1gw4v2IWK30XuKffsM8DG1R1EbDB/Y6ILAPuB5a70/yTiNh9FsaYovL5fITDYUsMgyhaclDV54C2foPvI9uhH+77e3OGr1fVpKoeBQ6RffqcMcaYMhDV4j3UTUQWAL/MKVbqUNX6nN/bVbVBRL4BvKSqP3CHfwt4SlUfG2SeDwIPAsyYMePG9evXjzq+7u5uampqRj19KU2kWMHiLTaLt7gme7zr1q3brqqr8o1TKddSg92rPmjWUtWHgYcBVq1apWvXrh31Qjdt2sRYpi+liRQrWLzFZvEObyw9rdr2LX1yOCsis1S1RURmAefc4SeBuTnjXQGcLnFsxphJor2tjSPbtkFPD1RXc+Xq1RXxEKSJpNT3OTwBPOB+fgB4PGf4/SISFJGFwCLs+RHGmFFwHIc3Nm/G2buX6hMncPbuzX6/zB6CNFZFSw4i8mNgK3CNiJwUkY8BXwbeISIHgXe431HV3cCjwB7gaeCTqmq3KhpjRiwWi9G+bx8rIhGWNjSwIhKhfe9eYrFYuUObUIpWrKSqvzvET3cOMf5DwEPFiscYc/nwiCBuYxtRxWNdcI9YpVRIG2PMuAiHw9QsWsTrhw9THwwSB+qWLCEcDpc7tAnF+lYyxpSV4zjEYrFxqxOIdnXh9Xo5mUqxq7OTtiuu4No1a+xGtxGyrWWMKZuRPJinEH1dca9qbMR/1110RKMcz2SIRCLjGPXlwa4cjDFlMZIH8xQqtytun8/HtIYGqjMZksnkOEZ+ebDkYIwpi/7PVKgKhQikUmM6kFtX3OPHkoMxpiyKcSAv9HkOZni2xYwxZVGsZypYV9zjw7aaMaZsinUgL+R5DiY/23rGmLIq1oF8LB3vGUsOxphJaLybyF6OrELamAlqPG8cK6fxvgmuGE1kL0d25WDMBNPe1kY8GqX5+ecr6qx4NMU4oz3Dz7esi01k6+sBt4lsNEoymbTipRGwLWXMBNJ3VhzyeFja2Eg8keBAUxORMncPMZqDfO4ZflV9fcHrMtyycpvIVoVCdq/DKFmxkjETSN9Zsbi9jI7HjWNjNdpinNHcBFfIsuxeh/FhW8uYCaTvrFgr6A7g0RbjjOYMv9Bl2b0OY2dXDsZMIH1nxYlMpmLOikd7p3O+M/yhKqlHsiyfz0c4HLbEMEq21YyZYBoaGqiKRJi3enVFnBWP5U7nwc7w89UpFOuuajOQbVFjJqhKenjNWIpxcm+CK6SS2oqMSsO2qjFmXIzHnc6F1ilY9xjFZ3UOxpiKYV1uVw5LDsaYsulf8WzNUCuHbXFjyqxUHcRVWkd0Q1U8j6ROodzrVO7lF9PkWhtjJphSdRBXSR3R9V0tHNm2jWXh8KAVz4XUKZR7ncq9/GKzYiVjyqRUHcRVUkd07W1t7N68mcMbN3Jq+3YSvb3AyO/0Lvc6lXv5pWDJwZgiytfjaDGeoTyYUi1nOH0H1HnA7FCIWX4/B3btwnGcARXPw/XUWu51KvfyS8GKlYwpkkrpIG40yylGWXoymaT1+HFOvPYa1b29tGUydNXW4mtpwReJXKx4bj13jsOvvkoYyIRCgxbXlLtzvXIvvxQsORhTBIX2ODr1mmvYs2cP1UW823ckdxU7jsOZlhbO7t1LdSYzrmXpqsrel17id2prmTFtGmc7O1nf0cHM1aupqakhHA7Teu4czz3yCMs8HrzV1dQtXEhzv+3Wl7hmr1jBgV27ynKndN823fPqq0hrK1pdzZWrV0+qSunJsybGlMFQZ9jD3cyVe1WBx0Pk2muZOXPmxX6FxvusvZAWQO1tbRzZto1T27ezoLqaOStWIB4Pu7ds4fo77yTkFqGMdFv06enpYXZjI22pFG0XLkAgwJSaGg689BIzw2F6PB4utLVxldfL0unTSSSTNB89imfevEG3W6/fz+wVK6iuri5PayHV7Ftpl1oylhyMGaV8xUb5ih0Gu6rYs2cPkUiE3mSS06+/XpQWMH0tgPrK8/uKQJLJJF6vl+YdO5jv8VBdU8PcYJBXNm/G5/PRmUqxU5Vrbr99yFgKablTW1tLTyBAtKuLiM9HRyxGSzTKvfX1NNTX09bRQfOhQ/RGIrRHo0Sqq0l3dBBzt2cikWD/li0sj0SI9F2N7drF8lE8y2KsCbjvb7gsHKZq6tSKea7GeJoca2FMiQ1VbFR1++2k02mCweCQRTmxWOySq4ruWIzDW7eS6ujg1PHj3HrttcyePn3QA07fQW20cg/i5+NxEGFaKERnOg09PdTPm8feZJIDTU3seeklIpEIscZG3r106cXiHeCSA+tIHtoTnjOH5jNnmOI4tKbTTJ0xg0hNDQC1NTUkEgnO9PTQ2tEBmQxts2Zxz003Ee3qYv/WrXTs3MmRqVOZt2QJDbW1o3rC23g0Qb0cnjY3OdbCmHGWe2YJDDjLHOzgED91ip0bNlDn9dLj8TDtmmu44vrrB3QdnXtV0d3Tw4bHH2daJkMwEmF+MknrsWNMnzIlb1FUPJWivb29oKet9cXuOA77t25leU0NVTU1pPbvJw4suuUW4okET73xBvNqaoj19HBy925Web1cMWUKx/x+Xtu0iWvuvptTJ09yctcugqkUEolw5erVBAKBYQ+U6XSanf/5n1SfOcOcpUsJX3EFN9bX8+wLL3C+o4OZ06bR3NLCkcOHuTocJhAI4Js5kyuuvZZIJML+F19keU0NLVOn0gg079uHb9myEVcCj/bpc/2NZ4V0pd5IVzmRGDOOxvIPN9TZde5ZZv+DQ7S7m1PHj/POm25C02le27KFn//gB8xZuJDapUu54c47aWxsBC6tzDy8dSuRdJpbbrgBEeGl48epCwRIplKkHGfIoqgzZ87Q7F6p9F1J5Cag3IrlQCpFSzRKOp1GDx/mWF0d9QsW4HccxOslFo/TE4uRjkZ5/JlnOLtzJzNTKZJTp+IPh/G2t9PV3c3W6mq8W7awKBol4/fjmzOHN3p6WH3PPfT6/US7u/F4vWTS6UvijsViJGMxbqqr40B9PY3pNMcOHeKsz0fGcXhu504iU6awd+tWlgcCzJ0yhYYZM2gPBgmEQnR1dRFIpYjU16NLlnBm3z5aLlwgHo1yze23j+jvO15n/OPVdXgl30hnycFMOvn+4fIljf537vr7nV2nHOeSs8zcg0On4zB34UKqQiF2bNpEfPt25kSjpNvbiTY3s/H4cd760Y9SW1tLMBikoaGBq2++mXQ0CvX19GYy1FZXU9vYyL5kkvTZs/jcM3OAtrY2PInExYOaAO2nT7P5sceIHT1KTzLJjBUrWPWOd4AqR7Zto6WpiXAmQ8jnI3boED09PfhCIRKqbHjmGcK1tfSEQpyLxYgeOYLj9XL1/PksXbSI6L59LIpE+PWOHXiCQc7U13NNby9nDx5k6bJl1IdCnDl/nr27d5Nct47ahQt55qc/JZxMEgsGWfn+9xPt6uLItm30trXRC7SpkkkmeeXgQXYdOsStt97KW2+9FZ/Pxy+ef55b581jWixGbTrNybNn0cZGYqrU1tbS4ibi2tpaZNky2rq7B60kH+6kYDzP+Mfadfh4XcUUS/kjMKNW6svR0S6vWHEmEgm6urqora29eJDI9w8X7eqieccOPIkEXek0M5cvZ9q0aYTDYRzHYdvTTyPJJK0HD7Lo5ptBlYjHgx9IplIE/X6caJTOzk4AYrEY05YsQUSYW1vL/i1baD51igOvvkrgwAFOdnezuqGBjOOQmjmTX3396yxasoTglCnMestb8Pl8nGpvx9/RQfPp0/QkEnQ1NrJk9Wr8Xi8KdLS3c+TVV9FolEP79xNMJPCr0tXby74nn2R6Os2UWIy5fj+7Xn+dzs5OZs2YwaJQiNqqKjInTnCsq4vA+fM4LS30BIO8FI3ii8U4XlVFuqaG9IkT1AWDzJk7l/PHjvGO66/n6fZ2ntq1iwtuQrph3jyOvfYax48d45XOThrq6gjOmUMyEsFxHLqOHuWdN92Ez+vFSac5cPAgx7q6mHrqFLXAuYYGNj3zDLctXky710swnebwzp0EfT66PB5ihw7RHgzS0dFBRISjySTJG27gvptvJuTe65B7ln7NbbcNSAyFnIWP98OCxtJ1eKXXW5Q/gn5E5B7gHwAv8K+q+uViLMdxHDKZDI7jjPkPMZKD32jH7W+4Jn2JRIK2tjaCwSB1dXUDmkjCwHL0fHEO9o9XVVU16ME5t8K0b7q+A/IV113HnDlzLi6zu7ubM2fOEIlEqKuru1iZ6zhOtjghEKC3t/eSZQAc3L+fVx57jHrHIVFdzYr3vpfZs2fjOM7Ff7hEIkG0p4dMPE5nZycHt25lLnBu3z7ad+/m19/8JgtuvJG6q6/GCQaZ1tJCOhAgqMqxXbu4ZvVqopkMsXSaqrY2Th06xImODl564QWiR47Qc+YMPX4/V69axYxlyzh/4QLt27dzePNmvMkkV/l8tMXjHD5zhmhbG9PTaU6/+ioX0mlO+v0EUim0txf1eKj1+Wjwejl58iTvXLWKxXPn0t7RwaMPP8w14XC2nmHvXr72+OMkzp5lzec+x7mmJjLd3Vw1dSqL5s0jpcqzP/4xN91yC1PmzKGjp4dpqRSZ9na6k0m6u7oIdneTjsepmjGDG6dP53RnJ3Pq6ph91VXMCYV4dN8+OmMxIlOnMuu22zhz4QJL/H4ywNELF7gpEqE6nSbS28t/7NvH6rvuwufzXSz26aPnz3Nhzx5unTOHUDDIEa+X/a+9RsepU8yvrqYmEGBZRwfs2EF3JoMmkzReeSXVySQnkknSs2ez4LrrLh7chztLH8lZeKU8LKjSb6SrqOQgIl7gH4F3ACeBV0XkCVXdM57L6Ttg9cZi7N68eUzlfCMpMxzLuOmc3/r/I5w+d47nvvtdFi9bRiYUQhsaeP0XvyB17Bi9Xi/z1q5l6W230XX0aN5y9KGWPXvFCk7v2nXJP97GX/yC7pYWalOpi8UI9XV1l1SYtra2cnrXLmb29nJ6715ihw7x1BNPcM273sX1b387x48c4dl/+Ad8LS20qVJ/7bW89Y47aIvHOXvkCJmODo4fOMA1S5YQmDmTle9/PwsXLqTl9Gme/vu/5+ZAgGBNDSrC01/5Cm+7916cUIieRII9hw6xd8sWAvE4exIJrj5zhqozZzh3/jwzHIcVjsOMSITOzk5OPP88U9/2Nq6dOpVUOs22ZJLD0SjOuXOcmzqV0wcPsvff/51IKES4sZHozp0sPniQht5ejomw78IFDm7cyEyfjwX19ajHwyy/n/PpNJ5MhiPRKPXHjnHVnDnMcRwOnztH67lzvH3RIrpiMY709jK1sZF3rlrFK/v2ceS556itq+O1bds4+uSTJAMBAj09hFXxdXVxXVUVNSLMb28nGosRq64m09vL6dZWFk2bxpRIhJkeDzGvlzdiMU53d+P3+Yim0yzyeOjyeJBYjLaWFqK9vSSmTGHGggW0njpFIhzm19EodQ0NJHt6uHbBApK7d3MqGqWnqopF8+dzqrWV1kiEupkzWbhqFeFweMBBLun34/V4UBGcdBoFZk6bxrK5c1leV8cTL75IOhymN5MhJMKMqVM5mU5z9ezZRLu7WXfvvfT4/ZecRec7Sx/pWXglPCyo0h95WhlRvOkm4JCqHgEQkfXAfcC4JYfcA2t7MsniUGjU5XwjOVsZ67in29ouXuXk/iM4jkPrsWMs83hYUFdHIpnke9/8JjeHQlw/dy6J3l7+88UXeam5mXtvv52qmhq69u+nisHL0Qdb9muvvEIYLinvPr55M+9bsYKZc+bQHo3y7KOPMn/FClbU1l6sMD388stUZTK0nTiB78wZ3jZ9OjO6uogfPszrXi/bn3iCt8fjLJs3j73Hj7Pl1VeR2bNxjh5ljsdDW2srH6+t5fTp0yxYsIDnH3uMxv/yX9i/dSvXqnLTjBlEk0me2bOHpeEwc6uqqKquZmtHB888+yx3V1URaGgg1NFBcscOZl55JZ5UijPnzjHN5yNQXU1EhDpVRIRYLEZdJEJEBFmxgvm3345u387K6dO5UFtLQzDI06+8gnPqFEsch5nV1UzNZNhx8iRLGxtZNns2c8NhzmYyNLqVqcdiMSLBID0+Hwvq6oh2ddETjzPT4+EttbU0x2I0x+PUAj2pFITDeBIJ9u3cif/ECWZlMrxDlVbHwZ9Os6O3l4wq1cCaSISNqRQvdHbSlkjgeL1Uz5vHkuuvp/nYMc739hJcs4b0668z68wZFixejJ4/T+3p0yQ9HqSqit5gkP3pNAtU6Zk9m2WrV3PT3Xez74UXmA20HjvGqfp6WsNhFsyaxRW1tQTmz6d+9my6MxnmzZs36EFu8a23kslk2HX8OEHHIdnYyNW33krswgXOxeNkgCsWL6YrFCKWydAdjzN10SIa/H6u9PmYUltLh+MUfBZd6WfhQ6mUq5jBiGrl3N8nIu8H7lHVj7vfPwLcrKp/kjPOg8CDADNmzLhx/fr1I1pGJpOht7ubkM9Ht+NQ4/ORcBwCNTV4PCPrhzB3Xn2GmtdYx+1KpaiJRC6OG49GCXk8qCrJWAyPx0Owuhonk6Hz/HnCPt/F6bt6e0n7fNTV1QHQ29ODB/BXVyMezyVxDLbseCpFBqj2+RARelMpou3tTI1EQASAzkQCb1UVNe4/Y7fj4FUlrYo3kcDrOPg9HhKqePx+Eh4PvV1d1Ivg83hwUim602m84TCSTiMipFMp6vx+4o5DIBIh5jiEpk4lHY/jdHUR9ngQETp7epBAgNqpU7MH+WSSVDxOXSBABkglk6TTafyRCMnubnoTCbweD9WBACkRkpkM3kiEoONAJkMCqJk6FY/Hk90WXi/Jnh78IrR3d5Pq7qY2k8Hv8ZBUpT2dptrvx+fzURMM0t7dDaqIx4MAjghpj4eIx0Ov45BMp/ED9aEQaVXakkm8fn+2vb/XS1cySdDvx+M4pHt7CTkOve7ZdyyTwS9CaOZMfGfPEhchIUIwFKI3nWbKtGmEgkEymQzxdJrq2locxyF6/jy+3l58gJNMkspkyHi9hMJhUkCgqgqP10swHMbr9ZJOp+nt6YFMhlQmg4hAOk28u5uA3w9eL9UNDQQCgQH7bt8+mk6nScZiaCZDSpWw14uIkIjFiMXj+L1eQoEAaSClSjAUItXbiz8QwOMmb6/XW/D/Y1/MooqKjHj6XN3d3dS4919MBCONd926ddtVdVW+cSonTWXJIMMuyV6q+jDwMMCqVat07dq1I1qA4zjs3ryZxaEQL3d2cnNdHQcSCZa/7W2junLom1ff2cpQ8xrruJva2nhbzrjt7e00NzXhSSQ40NzMHW95CzMaG2nv6OB7W7ZkrxwaGrJXDp2dpObO5a7bb6cqFGL7Sy9RBSzru3LIiWOoOPuKlgKpFF2BAGePH+ft06fT4N7N+mxX15tXDqEQG8+cYUooxPRly3j+e9/Dd/Qoc2pqWDJvHq2qdMyfz/YXXuDKaJRltbXsPX6c/fE4N993H2ePHiXkXjlc7fVyGlhw2208n0jw9ve8h8Mvv0y4t5fjL7xAprubnWfOcP9HPsKSmTOJJxK83tnJsf37ub6ujkh1NVubmog6Dms+/GEuBAI8uXUroUyGM+fOUTt7NvWLF5MRYU51NTHgqptuorGx8ZJtkaqu5tiuXTS3tLB9yxamtrQQSaXoVmVfby/3LVvGtKuu4vy5czS1teGZNo0aVRI9PaxcsQKP38+xnTs5F43S6fFwVV0drefP4wuHORIO07hwIUsaGohXVbHi3ntpO3iQwIEDRI8cwbN7Nz29vZzPZDjh83HqzBne9dnP0vOzn+GfORNnwQJmXHcdc1asoPf06WxxoHsTXl9xYWtrKzs3bKD74EHOHD3KFXPnct3q1TREIhxIJrn65psvaQrbtx/2r6NSVXp6egbUAeX7H0kmk2zdupVGEQKpFD1TphCaM4fu5mZ8ySRUVzPvhhuodg/mffVOo71zeTzOwjdt2sRIjy3lVIx4K+3K4Vbgi6p6t/v9CwCq+j8GG3/VqlW6bdu2ES+n78DaEosxKxweW52DO6+C6hHGMO554M4777xknL5/hJ6enosH7l6/H506ldcff/zSOofbb6fryJFsnYP7fN4h6xyGiDP3H+/EiRM0PfbYJU0X6+vrL053LJXiljvuoKGhgdbWVnZs2EDs0CH8Ph91S5Zw7Zo1HD9+nGe/+tVsnQNQv3x5ts4hkeDs4cNkOjs5vn//gDqHvviIxWhPpZixbNmbB0U33o6OjovxnU2nmXHVVcyfNu1iHUogELhYTBcOh3nhhRdYvXr1gINK7rbo8XiYsWwZ0WiUjd/5DvEjR3ACAaqXL6e2t5dGr5c24C333cdVV11FV1cX8XicjsOHCaRSXOjtZfo11xCJRDi+YwfEYiS9Xpa99a00NDRcUrnf3t7OG5s3c2rHDk4dPgzpNNMWLGD6W97ClKuu4uypU0x3HDw1NcxesaKgfpn6muq2t7dz4cCBce9YL59NmzZxxx13DLizuhKLU2DyJwcRGfbKAVWtmBfZK5kjwEIgALwGLB9q/BtvvFFHK5VK6YYNGzSVSo16Hrnz6u7uLmheox1348aNI5pvPB7XU6dO6fnz5y8Oyx1nuDgKiTMej+vZs2c1Ho8PmK5/vKlUSjs6OrSjo+OSeUajUT148KCeOXNG4/H4xWX2zbu9vX3AMgaLb7B4c+Mbbn3ybd+h5p27fQfbFvmmL2T79m2z8+fP6/nz5y/Zdhs3bix4Pyp0nYppuP230kz2eIFtOszxuKLStao6IvInwH+Qbcr6bVXdXYxl+Xw+PB7PuJyxjKTlQ6nGDYVCzJ49O+84+eZdyLJDodCAooWhpvPl1Hnkqqmp4eqrrx50/HzFFoOtS//l9o9vPNuy99+++eIdbPpCtu9Q26xPOBweLvS88660s3VTWSpu71DVJ4Enyx2HMcZczuwxocYYYwaw5GCMMWYASw7GGGMGsORgjDFmgIq6z2GkRKQVOD6GWUwDzo9TOMU2kWIFi7fYLN7imuzxzlfVxnwjTOjkMFYisk2HuxGkQkykWMHiLTaLt7gsXitWMsYYMwhLDsYYYwa43JPDw+UOYAQmUqxg8RabxVtcl328l3WdgzHGmMFd7lcOxhhjBmHJwRhjzACTJjmIyD0isl9EDonI5wf5XUTka+7vu0RkZc5vx0TkdRHZKSLbcoZPEZFnReSg+z5und6PNl4RucaNs+/VJSJ/6v72RRE5lfPbu0sY7xIR2SoiSRH5XCHTlnn7DhqviMwVkY0isldEdovIp3N+q9TtW9L9dwzbtlL33Q+7/2O7RGSLiFw33LRl3ncHjXfc993h+vSeCC+y3XsfBq7kzedALOs3zruBp8g+be4W4OWc344B0waZ798Bn3c/fx7420qIt998zpC9oQXgi8DnyrR9pwOrgYdyY8g3bZm371DxzgJWup8jwIGceCtu+5Z6/x1rrBW6794GNLif39X3v1bB++5Q8Y7rvjtZrhxuAg6p6hFV7QXWA/f1G+c+4Hua9RJQLyKzhpnvfcAj7udHgPdWWLx3AodVdSx3iY9LvKp6TlVfBVIjmLZs23eoeFW1RVWb3M9RYC8wZ5ziGvd4h1GM7TtesVbSvrtFVdvdry8BVxQwbTn33UHjHe99d7IkhznAiZzvJxm4UfKNo8AzIrJdRB7MGWeGqrZAdsOTPSOqhHj73A/8uN+wP3EvN789jpe6hcQymmnLuX2HJSILgBuAl3MGV9r2hdLuv+OybancffdjZK/Yh5u2Uvbd3HgvGo99d7IkBxlkWP82uvnGuV1VV5K9RPukiLxtPIMbxFjjRUQCwG8C/5bz+zeBq4DrgRbg78cUZYGxFHHa0RrzMkWkBvgp8Keq2uUOrsTtC6Xdf8dj21bkvisi68gebP9ypNOOo7HE2zd8XPbdyZIcTgJzc75fAZwudBxV7Xs/B/w72Us7gLN9RTnu+7lKiNf1LqBJVc/2DVDVs6qaVtUM8C+8uR6liHc005Zz+w5JRPxk/7l+qKo/6xteodu31PvvmGJ1Vdy+KyIrgH8F7lPVCwVMW9Z9d4h4x3XfnSzJ4VVgkYgsdM9K7gee6DfOE8BHJesWoFNVW0QkLCIRABEJA+8E3siZ5gH38wPA4+WON+f336XfZXm/Oon38eZ6lCLe0Uxbzu07KBER4FvAXlX9ar/fKm77lmH/Hcu+0Kei9l0RmQf8DPiIqh4ocNqy7btDxTvu++541LBXwots654DZGv6/9od9gngE+5nAf7R/f11YJU7/EqyLQJeA3b3Tev+NhXYABx036eUO173t2rgAlDXb57fd8fd5e5Qs0oY70yyZz1dQIf7uXaoaStg+w4aL3AH2cv4XcBO9/XuSt2+5dh/x7gvVOK++69Ae87fe1u+aStg3x003vHed637DGOMMQNMlmIlY4wx48iSgzHGmAEsORhjjBnAkoMxxpgBLDkYY4wZwJKDmTBEJO32KPmGiPybiFQXaTmrRORr7ue1InLbKObxpyLyUffzEjfuHSJy1Rhjuz63R00R+U0ZpOfOAufVKCJPjyUeM3lZcjATSVxVr1fVa4Fesm2/hyUivpEsRFW3qeqn3K9ryfaCWTB3eX8I/Mgd9F7gcVW9QVUP54wnIjLS/8HrybaD74v1CVX98gjn0TdtK9AiIrePZnozuVlyMBPV88DVku1X/+duh2Ivud0K9PVf/7CIPAN8T0Tmi8gGd7wN7l2miMgH3CuR10TkOXfYWhH5pWQ7L/sE8Gfumf9bReSo20UBIlIr2Wcp+PvF9nay3UM47ln+nwIfl2xf+wsk29/+PwFNwFwR+aaIbJNsH/xf6puJiKyWbH/9r4nIKyJSB/w/wAfdeD4oIr8vIt9wxx9qHb8r2WeDbBGRIyLy/pxYfw58eBz/LmaSsORgJhz3zPxdZO/4/BKwQ1VXAH8FfC9n1BvJ9j3zIeAbZLtAXwH8EPiaO85/A+5W1evIdgZ3kaoeA/4/4H+5VyzPA5uA33BHuR/4qar275r6dmC7O48nc+axzv39GjeWGzTbZfVfq+oqYAWwRkRWuF0n/AT4tBvbXUDMjfcnbjw/6bfcodYRsn393wHcC+ReaWwD3oox/VhyMBNJlYjsJHtAaybbj8wdZLsGQFV/DUx1z7ABnlDVuPv5Vt4s5vm+Ox3Ai8B3ReSPyD5oZTj/CvyB+/kPgO8MMs4soDXPPI5r9hkdfX5HRJqAHcByYBnZBNKi2ecioKpdquoME9tQ6wjwc1XNqOoeYEbO8HPA7GHmay5DIyqLNabM4qp6fe4At7Ox/vr6hInlmZcCqOonRORmslcDO0Xk+jzToKovukVDawCvqg7WgVkcCOWZzcW4RGQh8Dlgtaq2i8h33WmFsXcPnTt9Mudz7jYLufEacwm7cjAT3XO4ZeYishY4r2/2YZ9rC9liINzxX3CnuUpVX1bV/wac59LukgGiZB+5mOt7ZHsVHeyqAbJP4Lq6wPhrySaLThGZQba4DGAfMFtEVrtxRtzitMHi6TPoOg5jMePXA6qZRCw5mInui8AqEdlFtiz9gSHG+xTwB+54HwE+7Q7/ioi8LiJvkE00r/Wb7hfA+/oqpN1hPwQaGPgksz5PAQU9cEdVXyNbnLQb+DbZYi40+4jIDwJfF5HXgGfJnuVvBJb1VUgXuI75rAN+VUis5vJivbIaM0Jua5/7VPUjecb5d+AvVPVg6SIbObeF1n365jOJjQEsORgzIiLydbJFP+/WSx8M03+8a8g+Z/i5kgU3QiLSSPYRoz8vdyym8lhyMMYYM4DVORhjjBnAkoMxxpgBLDkYY4wZwJKDMcaYASw5GGOMGeD/B2sTUVf5qTZDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(por,perm,c = 'red',alpha=0.2, s = 20,edgecolors = 'black')\n",
    "plt.xlabel('Porosity (fraction)'); plt.ylabel('Permeability (mD)')\n",
    "plt.title('Permeability vs. Porosity'); plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's go through all the bivariate statistics listed above one-by-one.\n",
    "\n",
    "#### Bivariate Statistics\n",
    "\n",
    "We will cover a variety of measures of correlation.\n",
    "\n",
    "##### Covariance\n",
    "\n",
    "\\begin{equation}\n",
    "C_{X,Y} = \\frac{\\sum^n_{i=1}(x_i - \\overline{x})(y_i - \\overline{y})}{(n-1)}\n",
    "\\end{equation}\n",
    "\n",
    "We can use Numpy ot calculate the covariance matrix including:\n",
    "\n",
    "* sample variances on the diagonal\n",
    "* sample covariance on the off-diagonal \n",
    "\n",
    "We must put the two features into an $2 \\times n$ array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity sample variance is 0.00108557\n",
      "Permeability sample variance is 4156.4\n",
      "Sample covariance of permeability and porosity is 1.07\n"
     ]
    }
   ],
   "source": [
    "cov_matrix = np.cov(np.array([por,perm]))\n",
    "\n",
    "print('Porosity sample variance is ' + str(round(cov_matrix[0,0],8)))\n",
    "print('Permeability sample variance is ' + str(round(cov_matrix[1,1],2)))\n",
    "print('Sample covariance of permeability and porosity is ' + str(round(cov_matrix[0,1],2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Pearson Product Moment Correlation Coefficient\n",
    "\n",
    "\\begin{equation}\n",
    "C_{X,Y} = \\frac{\\sum^n_{i=1}(x_i - \\overline{x})(y_i - \\overline{y})}{(n-1)\\sigma_X \\sigma_Y}\n",
    "\\end{equation}\n",
    "\n",
    "We can use Numpy ot calculate the correlation matrix including:\n",
    "\n",
    "* sample variances on the diagonal\n",
    "* sample correlation on the off-diagonal \n",
    "\n",
    "We must put the two features into an $2 \\times n$ array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Standardized porosity sample variance is 1.0\n",
      "Standardized permeability sample variance is 1.0\n",
      "Sample Pearson product moment correlation coefficient of permeability and porosity is 0.51\n"
     ]
    }
   ],
   "source": [
    "corr_matrix = np.corrcoef(np.array([por,perm]))\n",
    "\n",
    "print('Standardized porosity sample variance is ' + str(round(corr_matrix[0,0],8)))\n",
    "print('Standardized permeability sample variance is ' + str(round(corr_matrix[1,1],2)))\n",
    "print('Sample Pearson product moment correlation coefficient of permeability and porosity is ' + str(round(corr_matrix[0,1],2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use Scipy to calculate the product moment correlation between 2 1D, paired arrays. \n",
    "\n",
    "* we also get the p-value for the significance of the measure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample Pearson product moment correlation coefficient of permeability and porosity is 0.51\n",
      "Sample Pearson product moment correlation coefficient of permeability and porosity p-value is 0.0\n"
     ]
    }
   ],
   "source": [
    "corr, corr_p_value = pearsonr(por,perm)\n",
    "\n",
    "print('Sample Pearson product moment correlation coefficient of permeability and porosity is ' + str(round(corr,2)))\n",
    "print('Sample Pearson product moment correlation coefficient of permeability and porosity p-value is ' + str(round(corr_p_value,10)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Spearman Rank Correlation Coefficient\n",
    "\n",
    "\\begin{equation}\n",
    "\\rho_{R_X,R_Y} = \\frac{\\sum^n_{i=1}(R_{x_i} - \\overline{R_x})(R_{y_i} - \\overline{R_y})}{(n-1)\\sigma_{R_X} \\sigma_{R_Y}}\n",
    "\\end{equation}\n",
    "\n",
    "We can use Scipy to calculate the rank correlation between 2 1D, paired arrays. \n",
    "\n",
    "* we also get the p-value for the significance of the measure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample Spearman rank correlation coefficient of permeability and porosity is 0.82\n",
      "Sample Spearman rank correlation coefficient of permeability and porosity p-value is 0.0\n"
     ]
    }
   ],
   "source": [
    "rank_corr = spearmanr(por,perm)[0]\n",
    "rank_corr_p_value = spearmanr(por,perm)[1]\n",
    "\n",
    "print('Sample Spearman rank correlation coefficient of permeability and porosity is ' + str(round(rank_corr,2)))\n",
    "print('Sample Spearman rank correlation coefficient of permeability and porosity p-value is ' + str(round(rank_corr_p_value,4)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comments\n",
    "\n",
    "This was a basic demonstration of bivariate 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 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}