GeostatsPy_Monte_Carlo_simulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "\n",
    "<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",
    "### Monte Carlo Simulation for Subsurface Data Analytics in Python \n",
    "\n",
    "\n",
    "#### Michael Pyrcz, Associate Professor, The 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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Monte Carlo Simulation\n",
    "\n",
    "Definition: random sampling from a distribution\n",
    "\n",
    "Procedure: \n",
    "\n",
    "1. Model the representative distribution (CDF)\n",
    "2. Draw a random value from a uniform [0,1] distribution (p-value)\n",
    "3. Apply the inverse of the CDF to calculate the associated realization\n",
    "\n",
    "In practice, Monte Carlo simulation refers to the workflow with multiple realizations drawn to buld an uncertainty model. \n",
    "\n",
    "\\begin{equation}\n",
    "X^\\ell = F_x(p^\\ell),  \\, \\forall \\, \\ell = 1,\\ldots, L\n",
    "\\end{equation}\n",
    "\n",
    "where $X^\\ell$ is the realization of the variable $X$ drawn from its CDF, $F_x$, with cumulative probability, p-value, $p^\\ell$.  \n",
    "\n",
    "It would be trivial to apply Monte Carlo simulation to a single variable, after many realizations one would get back the original distribution. The general approach is to:\n",
    "\n",
    "1. Model all distributions for the input, variables of interest $F_{x_1},\\ldots,F_{x_m}$.\n",
    "2. For each realization draw $p^\\ell_{1},\\ldots,p^\\ell_{m}$, p-values\n",
    "3. Apply the inverse of each distribution to calculate a realization of each variable, $X^\\ell_j = F_{x^\\ell_j}^{-1}(p^\\ell_j),  \\, \\forall \\, j = 1,\\ldots$, $m$ variables.\n",
    "4. Apply each set of variables for a $\\ell$ realization to the transfer function to calculate the ouptput realization, $Y^\\ell = F(X_1^\\ell,\\ldots,X_m^\\ell)$.\n",
    "\n",
    "Monte Carlo Simulation (MCS) is extremely powerful\n",
    "\n",
    "* Possible to easily simulate uncertainty models for complicated systems \n",
    "* Simulations are conducted by drawing values at random from specified uncertainty distributions for each variable\n",
    "* A single realization of each variable, $X_1^\\ell, X_2^\\ell,\\ldots,X_m^\\ell$ is applied to the transfer function to calculate the realization of the variable of interest (output, decision criteria):\n",
    "\n",
    "\\begin{equation}\n",
    "Y^\\ell = F(X_1^\\ell,\\ldots,X_m^\\ell), \\, \\forall \\, \\ell = 1,\\ldots, L\n",
    "\\end{equation}\n",
    "\n",
    "* The MCS method builds empirical uncertainty models by random sampling\n",
    "\n",
    "Let’s take a simple example, $OIP$ is oil-in-place calculated as the product of reservoir volume, $V$, average porosity, $\\overline{\\phi}$, and oil saturaton, $\\overline{s_o}$:\n",
    "\n",
    "\\begin{equation}\n",
    "OIP^\\ell = V^\\ell \\cdot \\overline{\\phi}^\\ell \\cdot \\overline{s_o}^\\ell, \\, \\forall \\, \\ell = 1,\\ldots, L\n",
    "\\end{equation}\n",
    "\n",
    "It would be difficult to directly calculate the OIP distribution as a combination of all these different distributions. \n",
    "* The distributions could all have different forms (parametric or non-parametric)\n",
    "* We use MCS to empirically work this out by sampling\n",
    "* Repeat to calculate enough realizations for analysis.\n",
    "\n",
    "How many realizations, $L$?\n",
    "\n",
    "The answer is enough! If the MCS computational cost is low then **many** is the right answer. If too few realizations are calculated then the summary statistics and the entire CDF of the output, decision criteria may be incorrect. This is caused by fluctuations due to not enough samples (see the 'Law of Small Numbers').\n",
    "\n",
    "The MCS method is very powerful. You can simulate output distributions that could not be calculated analytically.  \n",
    "\n",
    "#### Limitations\n",
    "\n",
    "The MCS method above assumes:\n",
    "1. **representativity** - the distribution is representative\n",
    "2. **independence** - the variables are independent of eachother\n",
    "3. **stationarity** - all realizations for each variable are from the same distribution \n",
    "\n",
    "#### Objective \n",
    "\n",
    "I want to provide hands-on experience with building data analytics and machine learning workflows. Python provides an excellent vehicle to accomplish this. Also, I have coded a package called GeostatsPy with GSLIB: Geostatistical Library (Deutsch and Journel, 1998) functionality that provides basic building blocks for building subsurface modeling workflows. \n",
    "\n",
    "The objective is to remove the hurdles of subsurface modeling workflow construction by providing building blocks and sufficient examples. This is not a coding class per se, but we need the ability to 'script' workflows working with numerical methods.    \n",
    "\n",
    "#### Getting Started\n",
    "\n",
    "Here's the steps to get setup in Python with the GeostatsPy package:\n",
    "\n",
    "1. Install Anaconda 3 on your machine (https://www.anaconda.com/download/). \n",
    "2. From Anaconda Navigator (within Anaconda3 group), go to the environment tab, click on base (root) green arrow and open a terminal. \n",
    "3. In the terminal type: pip install geostatspy. \n",
    "4. Open Jupyter and in the top block get started by copy and pasting the code block below from this Jupyter Notebook to start using the geostatspy functionality. \n",
    "\n",
    "There are examples below with these functions. You can go here to see a list of the available functions, https://git.io/fh4eX, other example workflows and source code. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geostatspy.GSLIB as GSLIB          # GSLIB utilies, visualization and wrapper\n",
    "import geostatspy.geostats as geostats    # GSLIB methods convert to Python        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also need some standard packages. These should have been installed with Anaconda 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np                        # ndarrys for gridded data\n",
    "import pandas as pd                       # DataFrames for tabular data\n",
    "import os                                 # set working directory, run executables\n",
    "import matplotlib.pyplot as plt           # for plotting\n",
    "from scipy import stats                   # summary statistics\n",
    "import math                               # trig etc.\n",
    "import random"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set the working directory\n",
    "\n",
    "I always like to do this so I don't lose files and to simplify subsequent read and writes (avoid including the full address each time). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#os.chdir(\"c:/PGE383\")                     # set the working directory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simple Monte Carlo Simulation Example\n",
    "\n",
    "Let's work out the OIP uncertainty example above. Let's work with parametric distributions for average porosity, $apor$, volume, $vol$ and oil saturation, $so$. Here's some parameters for these distributions. Note, these parameters could originate from a bootstrap workflow, analogs, etc. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "apor_mean = 0.15; apor_stdev = 0.02       # Gaussian mean and standard deviation porosity (fraction)\n",
    "vol_mu = 13.0; vol_sigma = 0.5            # LogNormal mu and sigma volume (m^3)\n",
    "so_min = 0.2; so_max = 0.4                # Uniform minimum and maximum oil saturation (fraction)\n",
    "L = 1000                                  # Number of MCS realizations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's set the minimum and maximum values for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "apor_min = 0.1; apor_max = 0.2            # average porosity min and max\n",
    "vol_min = 0.0; vol_max = 2.0              # vol min and max"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the *NumPy* package we have handy methods for Monte Carlo simulation from parametric distributions. We can actually draw all $L$ realizations at once for each variable and store them in ndarrays (each ndarray with realizations $\\ell = 1,\\ldots,L$).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "apor = np.random.normal(apor_mean, apor_stdev, size=L) # average porosity MCS simulation L times and store in array \n",
    "vol = np.random.lognormal(vol_mu, vol_sigma, size=L)/1.0e6   # volume ...\n",
    "so = np.random.uniform(so_min, so_max, size=L)         # saturation oil"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the distributions of the realizations of each variable to make sure the match the form of the parametric distributions that we selected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA54AAAGbCAYAAACs3U99AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdf7xldV3v8debAQRx+D0aMyOiRepkOtI4UqYRaiGVYOlV8yr+yNFSk/Lmr6tJVy3tZv64drUxDTAVyR9pliWpxLUUAh2RHxaKCDODMIqIM4Iwzuf+sdaBPWf2ObPPj3X2OXu/no/Hfpy9fuy1PmvtvT9nffb6ru9KVSFJkiRJUlf2GXYAkiRJkqTRZuEpSZIkSeqUhackSZIkqVMWnpIkSZKkTll4SpIkSZI6ZeEpSZIkSeqUhacWpSSfTHLasOPYmyRnJPmb9vnRSbYnWTbP63hkkv+cz2VKml6SE5JsHnYcXUpyTJJKsm873EneTXJ5khPme7mS5q6rY5cZrP+dSV7dPh9a3k1yTZLHtM9fmeSvOljHnds6riw850GS85N8N8ndhh3LfGi357Y2EX07yUeSHLWQMVTV46rqrDaeZyb53GyXleTMJLe323NTkvOSPGD+om1U1bVVdY+q+tFcltMeCP5Ez3L/X1Xdf+4RStMbpVyW5C+TnN1n/IOT/DDJ4cOIa6baA7Fdbf76fpL/TPKsLtbVm3dnq823r5u03J+qqvPnFJxG1ijlHYAkP5XkU+023ZzkkiQnD/ja85P8Vsfx3Vlgwfwdu0yxriT5gyRXJbk1ybVJ3tD7XlfV86vqtQMur5LsaPPhliR/3kXBXFV/XFVzeh/6HbvOZFtHlYXnHCU5BngkUMDjO1rHvl0sdy9eWFX3AH4SOBR480wXMKS4p/Kn7fasArYA7x5yPNKiMoK57Ezg15McNGn8M4BPVNVNCxjLXG1t89fBwO8B70rij1Fa8kYw7wD8PXAecC/gnsDvArcsxIoX2XEXwNuADTR5dznwOOBE4Nw5LPMhbT78BeDJwLPnGqQWjoXn3D0D+ALNQc6dTZSSHJ/kW72/xCR5QpJL2+f7JHl5kq8n+U6Scyd+ge9p/vScJNcCn2nH/227zO8luSDJT/Us+4gkf5/kliT/keR1vb+0JHlAe6bvpvYX8/82yMa1B2cfBh7ULueQJGcn2Zbkm0lelWSfdtozk/xbkjcnuQk4o93OV7Xz3ti+9pB2/gOS/E27/Te3cd+rnXZ+kt9K8kDgncDPtr9w3ZzkYUlu6E2wSX4jyaYBtudWmoS3tnd8kmcnubL9hfKfk9ynZ9pbk1zX7ttLkjyy37J73rd9k0zEO/G4Lck17Xzrk3y+3Zbrk7w9yf7ttAvaxX25fd2TM6npSZIHtvvn5jRN2B7fM+3MJH+R5B/SnB25MMmPt9PSvjc3tp+hS5M8aG/7TGNjpHJZVX2e5kem3+h57TLgN4GJ1hR3S/KWJFvbx1syxVmXTGqJkJ4zexPf0SQvbb9f1yc5NcnJSf6rjfWVPa+dcp9Npxr/CNwEPHiQfZLkV5J8qd2f1yU5Y6rlp+dsS5KJHDTxqLTNZad6/5JsAJ4GvLR9zd+343ubsE25z3v240t69uOzeuI7OckVbW7bkuR/7G2fadEbqbyT5EjgvsC7qur29vFvVfW5dvphST6R5hjqu+3z1e2019MU4W9vvz9vz6Tm8O18vd/TfsddP57kM+1++XaS9yU5tJ3/vcDRwN+363jp5HUkWZnk4+22fi3Jc3vWfUa7r89uv4eXJ1k3xb44Fvgd4GlV9fmq2llVl9Pk5JOSnNjOt0criUFU1deAf6PneC7NMeq729yxpX0fl7XTptwvfWI/I3ddRjXxfkw8dqbNoz2fwe+3uekJ7fg9jl37bWuS57b7+KZ2n6/smVZJnp/mbPF30xzbpZ32E0n+tf0sfzvJB2e6/4amqnzM4QF8jeaL9TPAHcC9eqZ9HXhsz/DfAi9vn59Ok2xXA3cD/hL4QDvtGJpf/84GDgIObMc/m+YXo7sBbwE29Sz7nPZxd2ANcB3wuXbaQe3ws4B9geOAbwM/NcU2nQ/8Vvv8SJqk/d52+GzgY20cxwD/BTynnfZMYCfwonY9B7Yxfw24H3AP4CM9y3oezS+DdweWtfvw4D4xPHNiW3pivAJ4XM/wR4GXTLE9ZwKv69kX7wW+3DP91DbGB7Zxvwr4957p/x04op32EuBbwAHttDOAv5n0vu07af37tdvzJ+3wzwDHt8s7BrgSOL1n/gJ+omf4BGBzz7K+BrwS2J/ml8PvA/fv2dabgPXt8t8HnNNO+2XgEpoz2Gm396hhf4d8LI4Ho5nL/ifwLz3DvwxsA/Zrh/9XG/s9gRXAvwOvbafd+b1rhyd/L8/krrxyAk3u+8P2O/rcdj3vb7fzp4DbgPvtbZ/12Ybe7/8+NGeFdgEPHWSftK//6fa1DwZuAE6d9P7s2w6fT5t3J8WwAfgqd+Xn6d6/O/dLz7hrgMcMuM93tvPsB5wM/AA4rJ1+PfDI9vlhwHHD/t74mNuDEcs7NP9brwI+QXNsca9J04+gKbzu3sbyt8Df9Uzf7Ts4+Ts6eR76H3f9BPDYdjtXABcAb+l5/Z3fx37rAP4V+L/AATRF3Tbg0e20M2hy2ck0x21/Anxhivf2+cA3p5j2r9x1THQmu+fSzf1e006/Mw8DD6DJCb/XM/3v2s/CQTQ55iLgee20gfcLPcd2k9Y/sT8m8u+TgJU0+fXJwA7a4yr6H7v2buuJNJ+j49qY/g9wwaRt/QTNMdvR7XpPaqd9gOb/2z7t+/Tzw/4uD/ydH3YAS/kB/DxNojyyHf7qpC/A64D3tM+Xtx/I+7TDV058kdvho9plTRQjRXuQMsW6D23nOaT98t9BW3z0rHsiaT4Z+H+TXv+XwGumWPb5NP/sb6Y5Y/C+9ku6DPghsKZn3ucB57fPnwlcO2lZnwZ+p2f4/j3b+Wyag44HTxHDdIXny4D3tc8Pb+PtW0S1X/Tb2u3ZBXyjd53AJ2mL53Z4n3Z595lied+laeoBgxWe7wD+AdhniuWdDny0Z3i6wvORNIXvPj3TPwCc0bOtf9Uz7WTgq+3zE2l+KDh+qlh8jOeD0c1lR7fLW90Ovw94a8/0rwMn9wz/MnBN+/zO7107vLfC81ZgWc8+KuDhPfNfwl0F35T7rM82nECTt26myb8/Yvcfqma6T94CvLl9PvH+TFl4tp+NG4Gf3Nv7N3m/9MxzDXcd0O1tn9/K7gfZNwLHt8+vpfmfc/CwvzM+5v5gdPPOauDt7Wd9F02Bc+wU864FvtszvNt3cPJ3dPI89Dnu6rOOU4Ev9Qzf+X2cvA7g3jQ5ZnnP9D8Bzmyfn8HuP+atAW6dYr2vYuqi9Byas8K75QwGKzxvaT8LRXP8c7d22r1ocuSBPfM/FfjsTPcLfQpPmuPga4CnTBPfJuCUnvdmusLz3TSXgU1Mu0f7OTymZ1t/vmf6udz1w8vZwEba/21L6WFT27k5DfhUVX27HX4/PU1F2uFfT9OM6NeBL1bVN9tp9wE+mqa55M00SfRHNF+cCddNPEmyLM0F2V9PcgvNhx+aM5IraBLGdf1e267r4RPratf3NODHptm2362qQ6tqVVU9raq2tevaH/hmz3zfpLlust96ofklaPL8+7bb+V7gn4Fz0jS5+tMk+00TU6+/AX4tyT2A/0bzT+H6aeb/s6o6lCbB3kpTAE+4D/DWnn1zE82vlqsA0jT7urJt0nAzzT+qIwcJMsnzaBLpb1bVrnbcT7bNa77Vvpd/POjyaPbndRPLak1+D77V8/wHNMmMqvoMzT/DvwBuSLIxycEDrlejbSRzWVVdS3PQ99/bXHEqbTPbVr/8tJLZ+U7d1TnHre3fG3qm30r7XWSwfdZra5u/Dqa5ZurEnmnT7pMkD0/y2TRN+75HcxZi0Px1b5qDndOq6r/acdO9f4PY2z7/TlXt7Bm+M4fRnCk6Gfhm28zsZwdcpxanUc07m6vqhVX14+1rd9AUCiS5e5qOz77ZxnEBcGjm1kHObsddSe6Z5Jy2qektNMdLM/l+3lRV3+8Zt7djjAPS/9rSb9P8INDPUe302TiOJic8GXg4zdlNaPb1fsD1Pe/TX9Kc+ZzTfmmPTT8EvL+qzukZ/4wkm3rW96BBl8mkXFhV24HvMMDxHPBSmuPUi9I0d14y17laeM5SkgNpCp5faAuIb9F0+vCQJA8BqKoraD5Uj6O5ruj9PYu4jqap6KE9jwOqakvPPNXz/DeBU4DH0BQ+x0yEQnP6fSfNr2wT7j1pXf86aV33qKrfnuFmf5vm15j79Iw7muasaL+YAbb2mX8ncENV3VFVf1RVa4CfA36V5nqPySYvk3Y/fR54AvB0miJ2r9oD0RfTFJoHtqOvo2mK0bt/Dqyqf09zPefLaN7rw9qDv+/R7Pdpta99Lc2vX9/rmfQOml92j62qg2maze51ea2twL3TXlfbmvweTKmq3lZVP0PT9O8ngT8YcL0aUWOQy86iySu/AXyjqr7YM61ffto6xXJ+QNM8bsJ0P9ztzSD7bA9V9UOafPTTSU7tWdZ0++T9wMeBe1fVITTXHQ2Svw6kabb2lqr6ZM+k6d4/6JOvJ5nJPt9NVf1HVZ1CcyD5d8ytgxIN0RjkHdptuI7mx96J/hReQvPD98Pb//+P6oljcszQFK0wfe6Z/Jo/acc9uF3Hf2f37/x039GtwOFJlveMG/gYY5LP0ByvrO8d2f6gdTxNi7hZqca5NMeBf9iOvo7mjOeRPe/TwVU1cS3v3vbLdP4PzWVNr+rZjvsA7wJeCBzRHh9exixzYZqO8I5ggH1dVd+qqudW1UqaViD/Nz19ECxmFp6zdyrNr2traJpKrKW5Zu7/sXvx9H6aHs0eRdOWf8I7gde3H1ySrEhyyjTrW07zhfoOTQL644kJ7a/sH6G5qPzuaW4V0hvDJ4CfTPL0JPu1j4elufh5YO16zm3jXt7G/vs0vxpN5QPA7yW5b3vG4Y+BD1bVziS/mOSn21/6bqEpavt1530DsDptBzw9zqb51eenaa7xHHQ7zqP5wm9oR70TeEXu6iDjkCRPaqctp/mHtA3YN8kf0px1mFabWD8IPGPiTEGP5TTbu719ryb/87qB5prYfi6k+Uf00vZ9PAH4NZpmK3uL6WHt2Y/92mXcRv/9rfEy6rnswzQHkX/E7mc7oclPr2pjPpLmAGaqfLYJ+M32zMlJND0qztZM99mdqup24E3cdbC1t32ynOYMxm3tAeBvDhjje2ia6f/ppPFTvn+t6fIXzGyf3ynJ/kmeluSQqrqDJoeav5aukcw7aToP+qM0nb/s037Gn01zPepEHLcCN6fpDOk1kxax2/enbW22habVxrL2zNaPT7OdE+vY3q5jFXv+wDzld7QtlP8d+JM0HUA+GHgOzWUKM9Ie+7wTeF+azqKWtcdZH6ZprvsvM11mH28ANiT5sbbV26eANyU5uN3/P55kIlfvbb/0labl2i/Q03KtdRBNcbmtne9Z3PUDA0x97Drh/cCzkqxNc1b/j4ELq+qaAWJ6UtpOqWgu/yqWSD608Jy904C/rub+R9+aeNA0ZXxa7mp28AGappaf6WlOAvBWml+hP5Xk+zRJ6eHTrO9sml/+ttB0rPOFSdNfSPMr3rdozv59gCbJ0jaZ+CXgKTQF17eAN9JczDxTL6IpWq4GPkfzxXnPNPO/p43nApprK29rlwHNr3YfojmAuJLmYvN+ByCfAS4HvpWkdx9+lLa5TVXt6PO66fxvmuLtblX1UZr9cU6a5heX0fzCCk1T4E/SXBv5zTb+yc2J+3k07fblrp7QLm+n/Q+ag7/v0/xaNrk3sjOAs9I03dit57z2oPPxbXzfpukA4BlV9dUBYjq4Xd932235DvBnA7xOo22kc1mbGyaKz8kHT68DLgYuBb4CfLEd18+LaX7kmWhm93fTbOPezHSfTfYe4OgkvzbAPvkd4H+16/lDBj9L+BTgCdm9N8dHsvf3793AmjZ/9dtHM9nnkz0duKbN08+nOWOhpWlU887tNGdT/4Xm2OaydjnPbKe/haYDoG+3MfzTpNe/FXhiml5M39aOey5NkfQdmtZK/z7NdkLzI9txNK2z/oGmqO71JzQ//tyc/j1DP7Xdhq00x1mvaX+wn40XAn9Fc2y3nWZ7z6ent/G5qKqv0Bw7ThSRz6C5JOwKmmOdD3FXc9+97ZepPJWmUN/akwtf2Z6RfxPNWdcbaE6C/FvP66Y6dp2I/dPAq2n+P11P84PCUwaM6WHAhUm203wPXlxV3xjwtUOVqr2dCdZSlOSNwI9V1Wl7nXkJS/J1mmay8/HLmaRFZlxymaTFw7wjdcMzniMizT2mHpzGepqmEQM3P12KkvwGTfOCzww7FknzYxxzmaThMu9IC6PTwjPJoUk+lOSraXoF/dkkh6e5Ce9V7d/DuoxhjCynaTawg6Yp1Zto7rc5kpKcT9NJzwsmtbmXZiTJe9LcrP6ynnH/u81blyb5aHpuMp3kFWlu+PyfSX55OFGPtLHKZVLX2mvrvpTkE+3wfZNc2B6HfXCaa9DGiXlHWgCdNrVNchbNbS7+qk1sd6fpwfOmqnpDkpfT9BT6ss6CkKRpJHkUzbUnZ1fVg9pxv0RzTdHOtskVVfWyJGtorv1ZT9MV+r/Q3N9wSVzUL2n8JPl9YB3N/U9/Ncm5wEeq6pwk7wS+XFXvGG6UksZBZ2c809wf8FE0HQ1QVbdX1c003VlP9Cx4Fk3PZpI0FFV1Ac29W3vHfaruupfgF7irm/1TgHOq6ofthfxfoylCJWnRaXu+/BWaDl5IEpr7wH6oncXjMEkLpt8NX+fL/Wi6GP7rNPdkuoSmV8B7tV0eU1XXJ7lnvxcn2UB7u4uDDjroZx7wgAd0GKqkxeaSSy75dlWtGHYcNF3hT/Q8vIrde0PczO43e76TOUwab4skh72F5rZjE/dlPAK4ueeHNXOYpL66yGFdFp770nRb/KKqujDJW4GXD/riqtoIbARYt25dXXzxxd1EKWlRSvLNRRDD/6S5j+vEbTj63Wy67/UK5jBpvA07hyX5VeDGqrokzT2fwRwmaUBd5LAuC8/NwOaqurAd/hBN4XlDkqPas51HATd2GIMkzUqS04BfBR5dd10Mv5nmfpATVtPc60ySFptHAI9PcjJwAM29nN8CHJpk3/aspzlM0oLp7BrP9kbA1yW5fzvq0TQ3dP04zY2Daf/aa5ikRSXJScDLgMdX1Q96Jn0ceEqSuyW5L3AscNEwYpSk6VTVK6pqdVUdQ3Nj+s9U1dOAzwJPbGfzOEzSgunyjCfAi4D3tT3aXg08i6bYPTfJc4BrgSd1HIMkTSnJB4ATgCOTbAZeA7wCuBtwXtMXB1+oqudX1eVtj5BX0DTBfYE92kpaYl4GnJPkdcCXaDuBlKSudVp4VtUmmi68J3t0l+uVpEFV1VP7jJ7yQKyqXg+8vruIJGl+VdX5wPnt86uxN25JQ9BZU1tJkiRJksDCU5IkSZLUMQtPSZIkSVKnLDwlSZIkSZ2y8JQkSZIkdcrCU5IkSZLUKQtPSZIkSVKnLDwlSZIkSZ2y8JQkSZIkdcrCU5IkSZLUKQtPSZIkSVKnLDwlSZIkSZ2y8JQkSZIkdcrCU5IkSZLUKQtPSZIkSVKnLDwlSZIkSZ2y8JQkSZIkdcrCU5IkSZLUKQtPSZIkSVKnLDwlSZIkSZ2y8JQkSZIkdcrCU5IkSZLUKQtPSZIkSVKnLDwlSZIkSZ2y8JQkSZIkdcrCU5IkSZLUKQtPSZIkSVKnLDwlSZIkSZ2y8JQkSZIkdcrCU5IkSZLUKQtPSZIkSVKnLDwlSZIkSZ2y8JQkSZIkdcrCU5IkSZLUKQtPSZKkEZPkgCQXJflyksuT/FE7/swk30iyqX2sHXasksbDvsMOQJIkSfPuh8CJVbU9yX7A55J8sp32B1X1oSHGJmkMWXhKkiSNmKoqYHs7uF/7qOFFJGnc2dRWkiRpBCVZlmQTcCNwXlVd2E56fZJLk7w5yd2GGKKkMWLhKUmSNIKq6kdVtRZYDaxP8iDgFcADgIcBhwMv6/faJBuSXJzk4m3bti1YzJJGl4WnJEnSCKuqm4HzgZOq6vpq/BD4a2D9FK/ZWFXrqmrdihUrFjBaSaPKwlOSJGnEJFmR5ND2+YHAY4CvJjmqHRfgVOCy4UUpaZzYuZAkSdLoOQo4K8kymhMN51bVJ5J8JskKIMAm4PnDDFLS+LDwlCRJGjFVdSnw0D7jTxxCOJJkU1tJkiRJUrcsPCVJkiRJnbLwlCRJkiR1ysJTkiRJktQpC09JkiRJUqcsPCVJkiRJneq08ExyTZKvJNmU5OJ23OFJzktyVfv3sC5jkKTpJHlPkhuTXNYzrm+eSuNtSb6W5NIkxw0vckmSpKVjIc54/mJVra2qde3wy4FPV9WxwKfbYUkaljOBkyaNmypPPQ44tn1sAN6xQDFKkiQtacNoansKcFb7/Czg1CHEIEkAVNUFwE2TRk+Vp04Bzq7GF4BDkxy1MJFKkiQtXV0XngV8KsklSTa04+5VVdcDtH/v2e+FSTYkuTjJxdu2bes4TEnazVR5ahVwXc98m9txezCHSZIk3aXrwvMRVXUcTfO0FyR51KAvrKqNVbWuqtatWLGiuwglaXDpM676zWgOkyRJukunhWdVbW3/3gh8FFgP3DDRNK39e2OXMUjSLEyVpzYD9+6ZbzWwdYFjkyRJWnI6KzyTHJRk+cRz4JeAy4CPA6e1s50GfKyrGCRplqbKUx8HntH2bns88L2JJrmSJEma2r4dLvtewEeTTKzn/VX1T0n+Azg3yXOAa4EndRiDJE0ryQeAE4Ajk2wGXgO8gf556h+Bk4GvAT8AnrXgAUuSJC1BnRWeVXU18JA+478DPLqr9UrSTFTVU6eYtEeeqqoCXtBtRJIkSaNnGLdTkSRJkiSNEQtPSZIkSVKnLDwlSZIkSZ2y8JQkSZIkdcrCU5IkSZLUKQtPSZIkSVKnLDwlSZIkSZ2y8JQkSZIkdcrCU5IkSZLUKQtPSZIkSVKnLDwlSZIkSZ3ad9gBSDO1fu0aNm/Zstu41atWcdGmK4YUkSRJkqTpWHhqydm8ZQtbX3vQbuNWvnrLFHNLkiRJGjab2kqSJEmSOuUZT0kaIzZVlyRJw2DhKUljxKbqkiRpGGxqK0mSJEnqlIWnJEmSJKlTFp6SJEmSpE5ZeEqSJEmSOmXhKUmSNGKSHJDkoiRfTnJ5kj9qx983yYVJrkrywST7DztWSePBwlOSJGn0/BA4saoeAqwFTkpyPPBG4M1VdSzwXeA5Q4xR0hix8JQkSRox1djeDu7XPgo4EfhQO/4s4NQhhCdpDFl4SpIkjaAky5JsAm4EzgO+DtxcVTvbWTYDq6Z47YYkFye5eNu2bQsTsKSRZuEpSZI0gqrqR1W1FlgNrAce2G+2KV67sarWVdW6FStWdBmmpDFh4SlJkjTCqupm4HzgeODQJPu2k1YDW4cVl6TxYuEpSZI0YpKsSHJo+/xA4DHAlcBngSe2s50GfGw4EUoaN/vufRZJkiQtMUcBZyVZRnOi4dyq+kSSK4BzkrwO+BLw7mEGKWl8WHhKkiSNmKq6FHhon/FX01zvKUkLysJTI+HWHdtZueKQ3catXrWKizZdMaSIJEmSJE2w8NRI2LlrF1tfe9Bu41a+esuQopEkSZLUy8JTkiRJWgTWr13D5i27/3BuCy6NCgtPSZIkaRHYvGWLLbg0srydiiRJkiSpU57x1LwaxyYi47jNkiRJ0kxYeGpejWMTkXHcZkmSJGkmbGorSZIkSeqUhackSZIkqVMWnpIkSZKkTll4SpIkSZI6ZedCkiRJksRo3a1gsW2LhackSZIkMVp3K1hs22JTW0mSJElSpyw8JUmSJEmdsvCUJEmSJHXKwlOSJEmS1CkLT0mSJElSpyw8JUmSJEmd6rzwTLIsyZeSfKIdvm+SC5NcleSDSfbvOgZJmo0kv5fk8iSXJflAkgPMYZIkSTO3EGc8Xwxc2TP8RuDNVXUs8F3gOQsQgyTNSJJVwO8C66rqQcAy4CmYwyRJ01i/dg0rVxyy22P92jXDDkuLyLh+RvbtcuFJVgO/Arwe+P0kAU4EfrOd5SzgDOAdXcYhSbO0L3BgkjuAuwPXYw6TJE1j85YtbH3tQbuNW/nqLUOKRovRuH5Guj7j+RbgpcCudvgI4Oaq2tkObwZWdRyDJM1YVW0B/gy4lqbg/B5wCeYwSZKkGeus8Ezyq8CNVXVJ7+g+s9YUr9+Q5OIkF2/btq2TGLW4jGuzAy1OSQ4DTgHuC6wEDgIe12dWc5gkSdJedNnU9hHA45OcDBwAHExzBvTQJPu2ZwxWA1v7vbiqNgIbAdatW9f3wE6jZVybHWjRegzwjaraBpDkI8DPYQ6TJEmasc7OeFbVK6pqdVUdQ9Mhx2eq6mnAZ4EntrOdBnysqxgkaQ6uBY5Pcvf2+vRHA1cw5jnMlgmSJGk2Ou1caAovA85J8jrgS8C7hxCDJE2rqi5M8iHgi8BOmny1EfgHxjiH2TJBkiTNxoIUnlV1PnB++/xqYP1CrFeS5qKqXgO8ZtJoc5gkSdIMDeOMp0bA+rVr2Lxlz7McO7Zvp+mDZeavH/S1861fLKtXreKiTVcseCySJGk8LKZjoVHjsd3iZOGpWenX3A5g+em3zPr1g752vtl0UJIkLbTFdCw0ajy2W5y6vo+nJEmSJGnMWXhKkiRJkjpl4SlJkiRJ6pSFpyRJkiSpU3YupJF1647trFxxyG7j7NFMGpw9LkpLV5J7A2cDPwbsAjZW1VuTnAE8F9jWzvrKqvrH4UQpaZxYeGpk7dy1yx7NpDmwx0VpSdsJvKSqvphkOXBJkvPaaW+uqj8bYmySxpCFpyRJ0oipquuB69vn309yJbBquFFJGmcWnpIkSSMsyTHAQ4ELgUcAL0zyDOBimrOi3+3zmg3ABoCjjz56wWLVwuh3KYWXI42vfp+HLlh4SpIkjagk9wA+DJxeVbckeQfwWqDav28Cnj35dVW1EdgIsG7dulq4iLUQ+l1K4eVI46vf59l8pRsAACAASURBVCG/Pf+X1tirrSRJ0ghKsh9N0fm+qvoIQFXdUFU/qqpdwLuA9cOMUdL4sPCUJEkaMUkCvBu4sqr+vGf8UT2zPQG4bKFjkzSebGorSZI0eh4BPB34SpJN7bhXAk9Nspamqe01wPOGE56kcWPhKUmSNGKq6nNA+kzynp2ShsKmtpIkSZKkTll4SpIkSZI6ZeEpSZIkSeqU13iOCW8ULEmSxt04Hg/12+Yd27cDB/V/QYfrHfV9relZeI4JbxQsSZLG3TgeD/Xb5uWn3zKU9Y76vtb0bGorSZIkSeqUhackSZIkqVMWnpIkSZKkTll4SpIkSZI6ZeEpSZIkSeqUhackSZIkqVPeTkWSxtytO7azcsUhe4xfiPu8SdJMeG9Ijbul/B2w8JSkMbdz16497rUGC3OfN0maCe8NqXG3lL8DAzW1TfKgrgORpLkyV0kaReY2SaNg0DOe70yyP3Am8P6qurm7kKTu9GtSaHPCkWKukjSKzG2SlryBCs+q+vkkxwLPBi5OchHw11V1XqfRSfOsX5NCmxOODnOVpFFkbpM0Cgbu1baqrgJeBbwM+AXgbUm+muTXuwpOkmbKXCVpFJnbJC11g17j+eAkbwauBE4Efq2qHtg+f3OH8UnSwMxVkkaRuU3SKBj0Gs+3A+8CXllVt06MrKqtSV7VSWSSNHPmKkmjyNwmackbtPA8Gbi1qn4EkGQf4ICq+kFVvbez6CRpZsxVkkaRuU3SkjfoNZ7/AhzYM3z3dpykOVi/dg0rVxyy22P92jXDDmspM1dJGkVLIrf5P03SdAY943lAVW2fGKiq7Unu3lFM0thYyjcBXqTMVZJG0ZLIbf5PkzSdQc947khy3MRAkp8Bbp1mfkkaBnOVpFFkbpO05A16xvN04G+TbG2HjwKe3E1IkjRr5qohuHXHdlauOGS3catXreKiTVcMKSJp5JjbJC15AxWeVfUfSR4A3B8I8NWquqPTyCRphsxVw7Fz1y6b10kdMrdJGgWDnvEEeBhwTPuahyahqs7uJCpJmj1zlaRRZG6TtKQNVHgmeS/w48Am4Eft6AJMeNqrfs3wdmzfDhzU/wXSLJmrJI0ic5ukUTDoGc91wJqqqi6D0Wjq1wxv+em3DCkajThzlaRRZG6TtOQN2qvtZcCPdRmIJM0Dc5WkUWRuk7TkDXrG80jgiiQXAT+cGFlVj+8kKkmaHXNVj/Vr17B5y+6d/NjMXVqShprbrrzi8s57ru6Xr4bZO/Z8xzPflx0ttv016hai9/Zx6CF+0MLzjC6DkKR5csawA1hM+t3M3Wbu0pJ0xjBXfscdt7P1tYftNm6+e67ul6+G2Tv2fMcz35cdLbb9NeoWovf2ceghftDbqfxrkvsAx1bVvyS5O7Cs29CkxaXfr4vgGaTFxFwlaRSZ2ySNgkF7tX0usAE4nKZXtVXAO4FHdxeatLj0+3URPIO0mJirJI0ic5ukUTBo50IvAB4B3AJQVVcB9+wqKEmaJXOVpFFkbpO05A1aeP6wqm6fGEiyL839o6aU5IAkFyX5cpLLk/xRO/6+SS5MclWSDybZf/bhS9JuZpyrppPk0CQfSvLVJFcm+dkkhyc5r81h5yU5bO9LEjTN1VeuOGS3x/q1a4YdlrQUzOY47N5JPtvmrsuTvLgdbw6TNBSDFp7/muSVwIFJHgv8LfD3e3nND4ETq+ohwFrgpCTHA28E3lxVxwLfBZ4zu9AlaQ+zyVXTeSvwT1X1AOAhwJXAy4FPtzns0+2wBjDRXL330e+6aUl7mE1u2wm8pKoeCBwPvCDJGsxhkoZk0MLz5cA24CvA84B/BF413Quqsb0d3K99FHAi8KF2/FnAqTOMWZKmMuNcNZUkBwOPAt4NUFW3V9XNwCk0uQvMYZIWxmyOw66vqi+2z79P88PZKsxhkoZk0F5tdwHvah8DS7IMuAT4CeAvgK8DN1fVznaWzTRJsN9rN9BcSM/RRx89k9VqDvr13HrH7bex3/4H7DZuoXpyne/7Xmm0zTZXTeF+NAd6f53kITS57MXAvarq+nZ91yfpe53VQucw79kpja655rYkxwAPBS5kFjls2T6ZzWolaTeD9mr7DfpcS1BV95vudVX1I2BtkkOBjwIP7DfbFK/dCGwEWLdu3ayv0dLMTHXfv21vPGKPcQthvu97pdE221w1hX2B44AXVdWFSd7KDJqkLXQO856d0uiaS25Lcg/gw8DpVXVLMlgR2ZvD9t9vmcdhkuZsoMITWNfz/ADgSTRdeg+kqm5Ocj7NNQaHJtm3Peu5Gtg66HIkaS/mlKsm2QxsrqoL2+EP0RSeNyQ5qj1TcBRw46yjlaTBzCq3JdmPpuh8X1V9pB1tDpM0FANd41lV3+l5bKmqt9BcqzmlJCvaM50kORB4DM31BZ8FntjOdhrwsVlHL0k9ZpOrplnWt4Drkty/HfVo4Arg4zS5C8xhkhbALI/DQnON+pVV9ec9k8xhkoZi0Ka2x/UM7kPzy9vyvbzsKOCs9jrPfYBzq+oTSa4AzknyOuBLtB13SNJczTJXTedFwPva2z5dDTyrXe65SZ4DXEtz5kGSOjPL3PYI4OnAV5Jsase9EngD5jBJQzBoU9s39TzfCVwD/LfpXlBVl9JcyD55/NXA+gHXK0kzMeNcNZ2q2sTuTdwmPHq2y5SkWZjNcdjngKku6DSHSVpwg/Zq+4tdByJJc2WuWjzskVqaP+Y2SaNg0Ka2vz/d9EnXDkjSUJirFg97pJbmj7lN0iiYSa+2D6O5IB3g14ALgOu6CEqSZslcJWkUmdskLXmDFp5HAsdV1fcBkpwB/G1V/VZXgUnSLJirJI0ic9si0O8SgtWrVnHRpisWfL3DvHRh/do1bN6yZbdx870f+q3jjttvY7/9D5jX9fZbz6D7dtDPw1zWMWoGLTyPBm7vGb4dOGbeo5GkuTFXSRpF5rZFoN8lBCtfvWWKubtd7zAvXdi8ZUvn+6HfOpaffgvb3njEvK53qvUMYtDPw1zWMWoGLTzfC1yU5KNAAU8Azu4sKkmaHXOVpFFkbpO05A3aq+3rk3wSeGQ76llV9aXuwpKkmTNXSRpF5jZJo2CfGcx7d+CWqnorsDnJfTuKSZLmwlwlaRSZ2yQtaQMVnkleA7wMeEU7aj/gb7oKSpJmw1wlaRSZ2ySNgkHPeD4BeDywA6CqtgLLuwpKkmbJXCVpFJnbJC15gxaet1dV0VzQTpLx6/9X0lJgrpI0isxtkpa8QXu1PTfJXwKHJnku8GzgXd2FpUEtxL2UxtFc75k1rHt9yVwlaSQt2dy2EP8P+62j3z0fYXzvnzhb3oNS82nQXm3/LMljgVuA+wN/WFXndRqZBrIQ91IaR3O9Z9aw7vU17sxVkkbRUs5tC/H/cKr/2ZPv+TgxXoPzHpSaT3stPJMsA/65qh4DLIkkJ2n8mKskjSJzm6RRsddrPKvqR8APkhyyt3klaVjMVZJGkblN0qgY9BrP24CvJDmPtkc1gKr63U6ikqTZMVdJGkXmNklL3qCF5z+0D0lazMxVkkaRuU3Skjdt4Znk6Kq6tqrOWqiAJGmmzFWSRpG5TdIo2ds1nn838STJhzuORZJmy1wlaRSZ2ySNjL0Vnul5fr8uA5GkOTBXSRpF5jZJI2Nv13jWFM8laTExV0kaRYs2t926YzsrV+ze0e6O7duBg/q/YBEbpW0Zlrnsw/Vr17B5y+73du1i/y/F93kpxjydvRWeD0lyC80vbge2z2mHq6oO7jQ6SRqMuUrSKFq0uW3nrl1sfe3uB7/LT79lirkXt1HalmGZyz7cvGXLguz/pfg+L8WYpzNt4VlVyxYqEEmaLXOVpFFkbpM0SvZ2jackSZIkSXNi4SlJkiRJ6pSFpyRJkiSpUxaekiRJkqROWXhKkiRJkjpl4SlJkiRJ6tTe7uOpETZqN6WVJEmLn8cfGlVL4bM9zBgtPMfYqN2UVpIkLX4ef2hULYXP9jBjtKmtJEmSJKlTFp6SJEmSpE5ZeEqSJEmSOmXhKUmSNGKSvCfJjUku6xl3RpItSTa1j5OHGaOk8WLnQiNoKfSoJUmSOnUm8Hbg7Enj31xVf7bw4UgadxaeI2gp9KglSZK6U1UXJDlm2HFI0gQLT0mSpPHxwiTPAC4GXlJV3+03U5INwAaAZftkAcNbWhZTK7N+sQwzHi2cpfLeW3hKkiSNh3cArwWq/fsm4Nn9ZqyqjcBGgP33W1YLFeBSs5hamfWLBWz1Ng6Wyntv50KSJEljoKpuqKofVdUu4F3A+mHHJGl8eMZzEVi/dg2bt2zZbdzqVau4aNMVQ4pIS4GfG0nSTCQ5qqqubwefAFw23fySNJ8sPBeBzVu27HF6fOWrt0wxt9TwcyNJmkqSDwAnAEcm2Qy8BjghyVqaprbXAM8bWoCSxo6FpyRJ0oipqqf2Gf3uBQ9EkloWnpKkoenXE1+/JuM2LZckaWmz8JQkDU2/nvj6NRm3abkkSUubvdpKkiRJkjrlGU9JkiQtOv2a4u/Yvh3Y836FmruF2N++p+PNwlOSJEmLTr+m+MtPv2VI0Yy+hdjfvqfjzaa2kiRJkqROdVZ4Jrl3ks8muTLJ5Ule3I4/PMl5Sa5q/x7WVQySNFdJliX5UpJPtMP3TXJhm8M+mGT/YccoSZK02HV5xnMn8JKqeiBwPPCCJGuAlwOfrqpjgU+3w5K0WL0YuLJn+I3Am9sc9l3gOUOJSpIkaQnprPCsquur6ovt8+/THLitAk4BzmpnOws4tasYJGkukqwGfgX4q3Y4wInAh9pZzGGSJEkDWJDOhZIcAzwUuBC4V1VdD01xmuSeU7xmA7AB4Oijj16IMKVO2ZPbkvQW4KXA8nb4CODmqtrZDm+m+UFtD+YwSZKku3ReeCa5B/Bh4PSquqU5YbB3VbUR2Aiwbt266i5CaWHYk9vSkuRXgRur6pIkJ0yM7jNr3/xkDpMkSbpLp4Vnkv1ois73VdVH2tE3JDmqPdt5FHBjlzFI0iw9Anh8kpOBA4CDac6AHppk3/as52pg6xBjlCRJWhK67NU2wLuBK6vqz3smfRw4rX1+GvCxrmKQpNmqqldU1eqqOgZ4CvCZqnoa8Fngie1s5rAOTDRL7300zdIlSdJS1eUZz0cATwe+kmRTO+6VwBuAc5M8B7gWeFKHMUjSfHsZcE6S1wFfovmBTfPIZumSJI2ezgrPqvoc/a+HAnh0V+uVpPlWVecD57fPrwbWDzMeSZKkpabL+3hKkiRJkmThKUmSJEnqloWnJEmSJKlTFp6SJEmSpE5ZeEqSJEmSOmXhKUmSJEnqlIWnJEmSJKlTFp6SJEmSpE5ZeEqSJEmSOmXhKUmSJEnqlIWnJEmSJKlT+w47AEm7u3XHdlauOGS3catXreKiTVcMKSJJkiRpbiw8pUVm565dbH3tQbuNW/nqLUOKRpIkSZo7m9pKkiRJkjrlGU9JWmLWr13D5i27nwXfsX07cFD/F0iSNIL6XZ4E/k9crCw8JWmJ2bxlyx7NsZeffsuQopEkaTj6XZ4E/k9crGxqK0mSJEnqlGc8lxCb10mSJElaiiw8lxCb10mSJElaimxqK0mSNGKSvCfJjUku6xl3eJLzklzV/j1smDFKGi8WnpIkSaPnTOCkSeNeDny6qo4FPt0OS9KCsPCUJEkaMVV1AXDTpNGnAGe1z88CTl3QoCSNNQtPSZKk8XCvqroeoP17zyHHI2mM2LnQItXvhrj2YCtJkhZCkg3ABoBl+2TI0UgaBRaei1S/G+Lag60kSZqDG5IcVVXXJzkKuHGqGatqI7ARYP/9ltVCBShpdNnUVpIkaTx8HDitfX4a8LEhxiJpzFh4SpIkjZgkHwA+D9w/yeYkzwHeADw2yVXAY9thSVoQNrWVJEkaMVX11CkmPXpBA5Gklmc8JUmSJEmdsvCUJEmSJHXKwlOSJEmS1CkLT0mSJElSpyw8JUmSJEmdsvCUJEmSJHXK26lII+TWHdtZueKQ3catXrWKizZdMaSIJEmSJAtPaaTs3LWLra89aLdxK1+9ZUjRSJIkSQ2b2kqSJEmSOmXhKUmSJEnqlIWnJEmSJKlTFp6SJEmSpE7ZuZC0BPTrrXbH9u3AQf1fIEmSJC0iFp7SEtCvt9rlp98ypGgkSZKkmbGprSRJkiSpUxaekiRJkqROWXhKkiRJkjpl4SlJkiRJ6pSdC0mSRsb6tWvYvGXLbuNWr1rFRZuuGFJEkiQJLDwlSSNk85Yte/QAvfLVW6aYW5IkLZTOmtomeU+SG5Nc1jPu8CTnJbmq/XtYV+uXpLlIcu8kn01yZZLLk7y4HW8ekyRJmqEur/E8Ezhp0riXA5+uqmOBT7fDkrQY7QReUlUPBI4HXpBkDeaxRePWHdtZueKQ3R47tm8fdliSJKmPzpraVtUFSY6ZNPoU4IT2+VnA+cDLuopBkmarqq4Hrm+ffz/JlcAqzGOLxs5du/ZoVrv89FuGFI0kSZrOQl/jea/2YI6quj7JPaeaMckGYAPA0UcfvUDhSdKe2h/RHgpcyIB5zBy2uNkJkSRJC2vRdi5UVRuBjQDr1q2rIYcjaUwluQfwYeD0qrolyUCvM4ctbnZCJEnSwlro+3jekOQogPbvjQu8fkkaWJL9aIrO91XVR9rR5jFJkqQZWujC8+PAae3z04CPLfD6JWkgaU5tvhu4sqr+vGeSeUySJGmGOmtqm+QDNB1wHJlkM/Aa4A3AuUmeA1wLPKmr9UvSHD0CeDrwlSSb2nGvxDwmSZI0Y132avvUKSY9uqt1StJ8qarPAVNd0GkekyRJmoFF27nQKOjXa+Idt9/GfvsfsNu45r5zu3dyIUlaWBP3BZ2sX97u1wOuPeVKkjQ1C88O9es1cfnpt7DtjUfsMU6SNFz97gsK/fN2vx5w7SlXkqSpLXTnQpIkSZKkMeMZT0nSSOvXhNZLHCRJWlgWnpKkkdavCa2XOEiStLBsaitJkiRJ6pRnPGfIXgslSdJSluQa4PvAj4CdVbVuuBFJGgcWnjNkr4WSJGkE/GJVfXvYQUgaHza1lSRJkiR1ysJTkiRpvBTwqSSXJNnQb4YkG5JcnOTiXbtqgcOTNIpsaitJkjReHlFVW5PcEzgvyVer6oLeGapqI7ARYP/9lll5Spozz3hKkiSNkara2v69EfgosH64EUkaB57xlEbcrTu2s3LFIbuNG7Qn5n69OM/k9ZKkxSXJQcA+VfX99vkvAf9ryGFJGgMWntKI27lr16x7Yu7Xi/NMXi9JWnTuBXw0CTTHge+vqn8abkiSxoGFpyRJ0pioqquBhww7Dknjx8JTkhaxfs2dd2zfDux5JloLp18Tdt8XSZKmZuEpSYtYv+bOy0+/ZUjRaEK/Juy+L5IkTc3CUxIwszNrc+mwSJIkSePHwlMSMLMza3PpsEiSJEnjx/t4SpIkSZI6ZeEpSZIkSeqUhackSZIkqVMWnpIkSZKkTtm5UKtfj5533H4b++1/wG7j+vXy2a+Hz6nmlSRJkqRxY+HZmqpHz21vPGKPcZP16+FzqnklSZIkadzY1FaSJEmS1CnPeEqSNGT9LvdYvWoVF226YkgRSZI0vyw8JUkasn6Xe6x89ZYp5pYkaemxqa0kSZIkqVOe8ZTUGZsPSpIkCSw8JXXI5oOSJEkCm9pKkiRJkjrmGU9JWiT6NU3esX07sOd9grU03LpjOytXHLLbOJubS5LGkYWnJC0S/ZomLz/9liFFo/mwc9cum5tLkoRNbSVJkiRJHVuyZzwH7S2z33x33H4b++1/wG7jbM6mcdKv+d8wvwP9vqeSJEkaHUu28By0t8ypmq5te+MRe4yTxkW/5n/D/A70+57mt/1OSpIkjQqb2kqSJEmSOrVkz3hKkjTK7BFXkjRKLDwlSVqE7BFXkjRKlkTheeUVlw/UEcpi6zBFGieDfv/8nmrcLcR3YNAO+CRJWihLovC8447b2fraw3Yb168jlMXWYYo0Tgb9/vk91bhbiO/AoB3wSZK0UOxcSJIkSZLUqSVxxlOSJA23wyGb70qS5sLCU5KkJWKYHQ7ZfFeSNBc2tZUkSZIkdcoznpIkLWFz6VH6jttvY7/9D9hjmYu9Ca3NfiVp6bHwlCRpCZtrj9Lb3njEHvMu9ia0NvuVpKXHpraSJEmSpE4N5YxnkpOAtwLLgL+qqjcMIw5Jmg1zmEbdfDff7Tdu0Kax/ZrV9oulC6PapNccJmkYFrzwTLIM+AvgscBm4D+SfLyqlnYWlzQWzGEaB/PdfLffuEGbxvZrVtsvli6MYpNec5ikYRlGU9v1wNeq6uqquh04BzhlCHFI0myYwyQtZeYwSUORqlrYFSZPBE6qqt9qh58OPLyqXjhpvg3AhnbwQcBlCxpod44Evj3sIObBqGwHuC2L1f2ravmwg5hsieawYX4uhv2ZdP2+98NiDpsfw34fXf/4rn+ctx06yGHDuMYzfcbtUf1W1UZgI0CSi6tqXdeBLYRR2ZZR2Q5wWxarJBcPO4YpLLkcNsz1j/O2j/v6x3nbJ9Y/rHXvxZLKYYvhfXT947n+cd72ifXP9zKH0dR2M3DvnuHVwNYhxCFJs2EOk7SUmcMkDcUwCs//AI5Nct8k+wNPAT4+hDgkaTbMYZKWMnOYpKFY8Ka2VbUzyQuBf6bpxvs9VXX5Xl62sfvIFsyobMuobAe4LYvVotyWJZrDhrn+cd72cV//OG/7Ylh/X0swhw17P7r+8V3/OG97J+tf8M6FJEmSJEnjZRhNbSVJkiRJY8TCU5IkSZLUqaEWnklOSvKfSb6W5OV9pj8qyReT7GzvO9U77bQkV7WP0xYu6v5muy1J1ib5fJLLk1ya5MkLG/me5vK+tNMPTrIlydsXJuKpzfEzdnSSTyW5MskVSY5ZqLj7meO2/Gn7GbsyyduS9OtOf8EMsC2/3+7zS5N8Osl9eqYtmu/+ANtxtyQfbKdf2PsZSvKKdvx/JvnljtY/3X78UZJN7WNWHYsMsP5nJtnWs57f6pk25/dxgPW/uWfd/5Xk5p5pc9r+JO9JcmOSvvc2TONtbWyXJjmuZ9p8bPve1v+0dr2XJvn3JA/pmXZNkq+02z7j7vIHWPcJSb7Xs3//sGfatO/ZPK3/D3rWfVn7Xh/eTpvTtrfLuHeSz7b59PIkL+4zT6fv/1zMMW/0jT3Jz7T79WvZy/+Y2a4/0xwvJTkzyTd63ve1HW1/37yRpqOmC9v98sE0nTbN57b/Ys96NyW5LcmpHWz783u+H59LsqZnWt//WTP5Ts92/Ukem+SSdtolSU7sec357TIntv+eHaz/mCS39qzjnT2vGeizP4d1P23Se79r4j2ez23vme+JSSrJup5xc37v71RVQ3nQXND+deB+wP7Al4E1k+Y5BngwcDbwxJ7xhwNXt38Pa58ftkS35SeBY9vnK4HrgUOX4rb0TH8r8H7g7cPajvnYFuB84LHt83sAd1+K2wL8HPBv7TKWAZ8HTljk2/KLE/sb+G3gg+3zRfPdH3A7fgd4Z/v8KT3bsaad/27AfdvlLFuo/dgOb1+A7X9mvzwwH+/jIOufNP+LaDpRma/tfxRwHHDZFNNPBj5Jc8/E44EL5/MzPMD6f25iucDjJtbfDl8DHNnhtp8AfGKu79ls1z9p3l8DPjNf294u4yjguPb5cuC/+nz2O33/5xB7J/kXuAj42XZ7Pwk8roP1T3m8BJxJn+OR+Vx/O9w3bwDnAk9pn78T+O35XnfPPIcDN/XMN5/bfnDP88cD/9Q+7/s/a5BlztP6HwqsbJ8/CNjSM9/5wLqOt/8Yps53e/3sz2Xdk+b5aeDqLra9nW85cAHwhYnlzsd73/sY5hnP9cDXqurqqrodOAc4pXeGqrqmqi4Fdk167S8D51XVTVX1XeA84KSFCHoKs96Wqvqvqrqqfb4VuBFYsTBh///2zjzejqLK498fCZCwBUgQIyBhR5BACKJIomGZwWFEZMIMIIhRRlBBREQdBCUjLiijgiyKgAQQCcMSBhElJBIWIbJlZREhRA1mQFCWACIwZ/445yb9bu76bt/3Hh/O9/Ppz62uru10ddc9VXWquiad1AuSxgIbAtP7orBN6LUsMdI02MxuinDLzOzFPip3LTqpFwOG4A3D6sCqwBPdL3JdWpHl5sL9no1/Zw4G1rvfVI44vzjcVwF7xWjo/sBUM3vZzB4DHon0Ss2/wX0sg1bkr0cZ9dhu/ocAl7eZR13M7FZc+avH/sAl5swG1pU0kpKe4Wb5m9kdkT6UXPctyF6PTp6Z3uZfar1H/kvN7L5wPw88CGxUFayr9d8Bpbe/Idc6ZnanuaZ6CfDBsvMvSV8qvd2MNn1PvI0Hb/NryV9W3gcCv+iFTtJK/s8VTtfE9Qeo/5/Vzjvd6/zNbE7UOcD9wBBJq7che0f516ONZ7+svHvbnrVaT6cC3wb+VvAro+6X058dz42APxbOl7Byw92NuN2glPJI2hXvHDxaUrl6Q69lkbQK8B3g810oV2/opF62Bp6RdI2kOZJOlzSo9BK2Tq9lMbM7gZvx0eGlwI1m9mDpJWyddmU5Ah9F7E3cbtJKWZaHMbNXgWeB4S3GLSP/IsX7CP7HfY+k2RWTrS7lPzHMxq6SVPlofZ/KH+ZqmwG/Knh3Kn9vy9cfz3B13RswPUzWjuxSnrtJmifpF5K2D78+lV3SGnin7uqCd6myy83nxwC/qbo0kOq/lXLVo5X2d6Nwt5Jmp+0WUFdf+nq0Nd9r0CnpRrs5HHgm2vhGaZYiO249U935KE12SUdLehTvgBzbJG47MnWSf5GJwBwze7ngd1GYmn65nqlrCflvFvrgLZLGF9Js5dkvS/aDWLnuS5Fd0hhgEzO7vsW4vWrL+rPjWevmtPptl07idoOOyxOjJpcCHzWzlWYS+5BOZPkUcIOZ/bFpyL6hE1kGA+OBE4B34KYEk8opVq/otSySrsGLiAAAEsJJREFUtgTeho+cbgTsKek9JZatXVqWRdJhwC7A6e3G7QNaKUu9MGXI0cl9BHirme0CfAg4Q9IWXcj/Z8AoMxsNzGDF7G+fyo8raleZ2WsFv07l7235+vQZlrQHrsB+seC9u5ntjJvgHt2F9uA+YFMz2xE4C7i2UpwaYbv5/u4H/NrMirOjpckuaS28U3tc1WwFDJD6r0E32t92ZOq03aqnL50IbIv/X69Pz+e9zPxrtRutplmW7Dvg32CtUKrsZnaOmW0R6ZzcJG7pdV8nf0/AB7G+BRxV8D7UzHbAdbbxwIe7kP9SvO7HAMcDP5W0Tqtpdpi3JyC9E3jRzIpr20uRPSaOvgd8ro24vWrL+rPjuQTYpHC+MfCnOmHLjNsNOipPPLw/B04Ok5z+pBNZdgOOkbQY+C/gcEmnlVu8tuj0GZsTJgSv4orTzk3idJNOZDkAmG1uLrwMH0F9V8nla4eWZJG0N3AS8IHCyOZAevdbKcvyMJIGA8NwE8Ey5OjkPlZM1TCzRfg6kTFl529mTxfyPB8Y207ZO82/wEozBCXI39vy9dkzLGk0cAGwv5k9XfEvyP4kMI32zbwbYmbPRVuDmd0ArCppBH3//jaq945kl7Qq3um8zMyuqRGk3+u/Dt1of5fQ0yS0kUwdtVv19KUwf7YIexH167Ub7eZTuCn14EZpdpp38G/ANDN7pVCmUmUvMJUVZqNlPM+d5I+kjfF39nAzWz7TbWaPx+/z+P4ipcsfZqZPh/tefKZ9a1p/9juSPajVnpUl+9r42tlZocO/C7guNhgqty2zJotAu3XgM0qLcPOnyqLU7euEncLKmws9hi9uXy/c679OZVkNmImPmPZbfZQhS9W1SfT/5kKd1MugCL9BnF8EHP06leUgfLZpML6+cyaw30CWBf8zf5TYSKLgP2De/RblOJqemwv9d7i3p+di/UW0v7lQJ/dxPWD1cI8Afkebm7y0mP/IgrsyAFJKPbb6TgDb4BvKqEz5I+4o6m848c/03FzmrrKf4Sb5vxVfi/PuKv81gbUL7juA95Wc95sr9xtXhP4Q96HldqyT/ON6ZZBnzS7ILnwt1xkNwnS9/nt537rS/gJ3h5yVDVb27UL+dfUloq2J/M8ATutC/nXbDeBKem4u9Kky8y5cnw3s0UXZtyq49wPuCXfN/6xW0iwp/3Uj/MQaaY4I96r4OttPdCH/DYj/aNwC7nHaePY7yTvOV8E7ept3S/aq8LNYsblQx3XfI+1mAbp54Lu+PYy/ZCeF31fxUR5ws4ElwAvA08D9hbgfw/9UH8HNLV6XsgCHAa8AcwvHTq9HWarSmEQ/dzxLeMb+AZgPLMA7c6u9HmWJBuI8fAOMB4Dvvg7qZQa+AVLlnbiuEHfAvPstyDEEV0gewXe+K/5pnBTxfkudHSC7dR/xHU8X4H8UC4AjupT/N/GNIObh64y3LbMem+Uf55OpUsTKkB8feV6Kt99LcHPWTxB//LgSck6UbQGFnQdLkr1Z/hcAfy3UfUWB2jzknhd1c1IX8j6mUO+zKXR+a9VZ2flHmEn4hhjFeB3LHumMw03K5hfu7759Wf+dHM3eG3rR/uJmoQsjzbMpDPSUlT8N9CV8/faCKMNPgLW6kH/ddiOerbvivlxJdFBLvvej8A7PKlVplin7mfi7MRdvs7cvxK35n1UrzbLzx81OX6iq+zfhA0j34u/i/RG/7iBuB/lPZEWbdh+FwXtafPY7vPcTiIHbgl+psleFnUXPNqvjuq8clRHJJEmSJEmSJEmSJOkK/bnGM0mSJEmSJEmSJHkDkB3PJEmSJEmSJEmSpKtkxzNJkiRJkiRJkiTpKtnxTJIkSZIkSZIkSbpKdjyTJEmSJEmSJEmSrpIdzwGMpAMkmaRt+7sszZA0QdKzkuZIelDSKV3O7wZJ68bxqV7EHynp+sL55ZLmS/psCWX7UtX5HR2kNVXSVp2WKUmSJEmS8pB0kqT7Q3eYK+mdTcJPkvSWEvMfJelDhfNdJH2/pLSHSbpE0qNxXCJpWFx7i6Srwj2hqEtVpTFG0gXhXl3SjLhPB3VYth56X7E8vUxvhqT1OilT0jrZ8RzYHALcjn94vmMkDSojnQbcZmZj8G8aHSZpbCuRJA1uNyMz29fMnsE/Ktx2xxM4Hjg/8n8z/p250Wb2vU7LBvToeJrZu3uRRoUfAF/oIH6SJEmSJCUiaTfg/cDOZjYa2Bv4Y5Nok4C2Op5NdJBRwPKOp5ndY2bHtpN+Ay4EFpnZFma2BfAY/m1gzOxPZnZgC2l8CTgr3GOAVc1sJzO7ohioF7ppD72vjfLU41J6p0cmvSA7ngMUSWsBu+Mfxj644H+FpH0L51MkTZQ0SNLpku6O0bej4voESTdL+in+gWEkXSvp3hipO7KQ1hGSHpY0S9L5ks4O/w0kXR1p3y1p90ZlN7MX8I/abiFpiKSLJC2I2dA9Is1Jkq6U9DNgupzTJS2MsAdFuJGSbo1RsoWSxof/YkkjgNMin7kR/1JJ+xdkukzSB2oUcyLwy3BPB94UaYwP+b8h6RbgM5L2k/SbKP8MSRtW6qgg2/yoh9OAoZHWZRFuWfzWk3FC5HmVpIeizIqy3Qbs3csOcJIkr2MkrSLp65LOkvSRDtN6m6QfRjvzybLKmCRvUEYCT5nZywBm9pSZ/QlA0ldCV1oo6Ufx338gPih/WegHQwt6TGW2cla4J0e86cAl8pnN2yTdF0dlMPs0YHyk91kVZh8lrR+63nxJsyWNLqT949A5FklaqaMqaUtgLHBqwfurwC6StojyLGx0cyStDYw2s3mS3gT8BNgpyrpFyP4VSbcD/yrp43HP5oW+uUaks6GkaeE/L2Sv1vuWl0eNdc5rJP1S0u8kfbtQ3OvwiZ6kLzCzPAbgARwGXBjuO/BRNYADgIvDvRo+wjYUOBI4OfxXB+4BNgMmAC8AmxXSXj9+hwILgeH4KNxiYH1gVbzDc3aE+ykwLtxvBR6sUd4JwPXhHh5pbQ98Drgo/LcF/gAMwUf+lhTKMhG4CRgEbBjhRkb8kyLMIGDtcC8GRuAjfgsL5XgvcG24h+GjdIOryroZcG/hvDqNWcC5hfP1AIX734HvhPtbwBnFcPG7rCq/ZU1knAA8C2yMDwbdWbnfEe8mYGx/P5N55DGQjnhP96nyO6747taIs6ybZaqT51DgFmBQnP8YeLLY5oT/+4DfAo8A/xF+BwBTgO8Ce5VUnlXw2YzVgFur28c88sij+QGsBcwFHgbOBd5buLZ+wX0psF+4ZwG7FK4tBkaEexdgVrgn44P3Q+N8DWBIuLcC7gn3BELvqj7HZxpPCfeewNxC2nfgeuII4Gl8JrIo2weAaTVknhbXRlXar+oyFMLuAVxdq2wF2b9QOB9ecH8N+HS4rwCOC/cgXK9bnn/4F8vTSOdcFPGHAL8HNimk8btiGfLo3pEzngOXQ4Cp4Z7KitGYXwB7Slod+CfgVjN7CfhH4HBJc4Hf4J2/ytrAu8zssULax0qaB8wGNolwuwK3mNlfzOwV4MpC+L2BsyPt64B1YjSrmvGS5uAziKeZ2f3AOLzhxcwewl/2rSP8TWb2l3CPAy43s9fM7AlcUXsHcDfwUUmTgR3M7PlGN83MbgG2jBG2Q/CG79WqYCOBPzdKB2/sKmwM3ChpAfB5vEMNfl/OKeT91yZp1pMRvI6WmNn/4X9mowrxnqRN85wkeQNwOSsvQzg4/AcSHwOuMbPX4nwK3slcjtzU7By8Td8OOETSdsA2wJ1mdjzQ8SxlWH/cDsw0s78DM4GO1lslyRsRM1uGzwoeiesTV0iaFJf3CCupBXinb/vaqTTkutDtwCcDzo/0rsTbiGYUda9fAcMVazSBn5vZy2b2FK5fbFgVV4DVSLOefy3a1bPeHrO6C4BDWXHP9sSXHBG607NN0mykc840s2fN7G/AA8CmhXipZ/UR2fEcgEgajr9sF0hajHd2DpKkeGFmAfvgCkOlcyp8hGinODYzs+lx7YVC2hPwDtNuZrYjMAcf/amYdtZilQhfSXujOh3A28xsjJmNNbMfFspVjxcK7prhzOxW4D3A48Clkg5vkF6FS/GG66PARTWuv4TL3Ihi2c7CZ393AI4qxG2nEa6Er8fLBfdrQNG0dghe5iRJVnAV8P4YhEPSKFxxuF3S8WHmtlDScdURq03FJJ0Qg1uVaw9JuiDiXyZpb0m/DhOtXSPcYZLuCnOv81R/ndKhwP9UTqJN+0tVmF2BR8xsUXQIpwL741YhlQGt16iBfMnC2ZJul/R7SePkG4E8LOnCYlgzu858zfmh4XVtwZ0kSRtER2iWmZ0CHANMlDQEnwE9MHSG86mvb7zKCj28OkxRB/ks8ASwIz4zuloLxaulb1T0lUb6BsD9wBhJy/sI4d4ReLCFvKF9PWsKcEzcs/9sIW49Us8a4GTHc2ByIHCJmW1qZqPMbBPcZHRcXJ+Kd6rGAzeG343AJyWtCiBpa0lr1kh7GPBXM3tRvlvuu8L/LuC9ktaTryecWIgzHW9UibR3akOWWwnFRtLWuKnub+uEO0i+VnUDvLN5l6RNgSfN7HzcPGznqnjPA9Wzr1Nwkzti1rWah+k5o9iMYXjHF6C4zqr6vlR2RXulUg9V1JSxhfy3xv8IkiQJzOxp/P2pzB4ejI+g74y3j+/E27ePSxrTZvJbAmcCo3FzrQ/h7e8JwJckvQ0f+NvdzHbClZiVOnCSVgM2N7PFTfLbiJ4bkywJv2uAfSSdhbcftdgB3wRkHHAx3k5+EXg78C+FjvkESd+XdB5wQ8RdyAqriyRJWkTSNuq54/xO+OxapcP0lHyvjuKmN9X6ymJ81hR66lzVDAOWhkXUh3GT01rpFSnqXhPw9ajPNchjOWb2CD4pcXLB+2TgvrjWCg/i7WirrA0sDd2p2JbOJKw9Qndah9blbqRzLkeSgDfj9ZF0mex4DkwOwW3pi1zNit3LpuOdlhkxOg6+29gDwH0xkn8eK49igW+oM1jSfHzh+GwAM3sc+AZupjsj0qqYNByLLyqfL+kB4BNtyHIuMCjMJ64AJlksxq9iGjAfmAf8Crf9/198XcDcMOGdiCuDywnl89cxM3F6+D2BN3q1Zjsx3/zoUfkC+laYDFwp6TbgqYL/14D1Iu95+JoGgB8B8xWbC7UgY13kGxm9ZGZLWyxrkryRKJrbVsxsx+Hrk14Ic7hr8EG6dnjMzBaEonc/bqJl+AZto4C9cIXx7liCsBeweY10RgDPtJBfzdkJM3vRzI4ws0+b2TkrRfLZlXWBM8LrJXxvgKXx3/Ai8PdIbJaZHWtmR1XSCvPfv9dZOpEkSX3WAi6W9EDoU9sBk8132z8fbyuuxZcLVZgC/DCsJIbiM3tnhm5R06IhOBf4iKTZ+EB0ZaZwPvCqfNOd6k/BTSb0NnwznnY3JzsC2FrSI5IejXyPaDVymLkOa6Nt+TKuf94EPFTw/wxuurwAX/e6fS29r0CrOmeRscDsGsuyki5Q2TAlSZC0lpktixnPacCPzay6Azzgke+GtgDfkKnmegBJB+Ab9pxc6/pAIf5MnjOzC5sGTpI3GDGjsAif9bzczLYJ09r1zewrEeZU4M9m9n1Jy8xsLUkbA9PNbLsIczK+yc7kMNm93szeHtemxPlVlWv4wN5bzOzEJuVbD5hjZqOq/Kvz2A1XWveJ8xMBzOybTdIfC5xqZvvG+VTgB2Z2S8g4zcwazmhKegoYGWv7kyRJSiH0l+fN7IL+LksjJJ2Jr6md2d9leSOQM55Jkckxer8QN+29tp/L0zaS9sZHy85qtAg9OtSL+6pcHfAMbj6XJEkVMaM5C98ptrKp0K3AByWtEcsNDsB36S7yBP4JpeFhivr+NrOeCRwo38Ss8umCTasDxYZjg2JmshF3A1tJ2izMcw/GN3Jrxg64BUWF0fgsCPh6rPkrxSgg30/gz9npTJKkC/yAnusqByoLs9PZd+S3AZPlmNkJ/V2GTjGzGbhNfythB/QoHICZ1TQXTpJkOZfj5rQHA5jZfTFLWVk/fYGZzSlGMLNXJH0VN+16jJ6mXU0xswdilnR6bLrxCnA0vsarmum4+e8MAEmX40sIRkhagn/y4EJJx+Br9Qfh1iatrOveoSJndG6H2ordtYud0HrswYr1nkmSJKURm2Fe2t/laEbsIZL0EWlqmyRJkiRdIjY2Ot7MPtzfZalG0jXAiWbWcPONJEmSJCmDNLVNkiRJki4Rs603q/7nVvqFMOm9NjudSZIkSV+RM55JkiRJkiRJkiRJV8kZzyRJkiRJkiRJkqSrZMczSZIkSZIkSZIk6SrZ8UySJEmSJEmSJEm6SnY8kyRJkiRJkiRJkq6SHc8kSZIkSZIkSZKkq2THM0mSJEmSJEmSJOkq2fFMkiRJkiRJkiRJusr/A7qfRu/R3SZNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.subplot(131)\n",
    "GSLIB.hist_st(apor,apor_min,apor_max,log=False,cumul=False,bins=50,weights=None,xlabel=\"Average Porosity (fraction)\",title=\"Average Porosity Realizations\")\n",
    "plt.ylim(0.0,60)\n",
    "\n",
    "plt.subplot(132)\n",
    "GSLIB.hist_st(vol,vol_min,vol_max,log=False,cumul=False,bins=50,weights=None,xlabel=\"Volume ($10^6$ $m^3$)\",title=\"Average Volume Realizations\")\n",
    "plt.ylim(0.0,120)\n",
    "\n",
    "plt.subplot(133)\n",
    "GSLIB.hist_st(so,so_min,so_max,log=False,cumul=False,bins=50,weights=None,xlabel=\"Saturation Oil (fraction)\",title=\"Average Saturation Oil Realizations\")\n",
    "plt.ylim(0.0,40)\n",
    "\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.2, hspace=0.2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks good, the shapes are Gaussian, lognormal and uniform and the central tendency and dispersion make sense given the parameters that we selected.\n",
    "\n",
    "Now we can use broadcast methods to calculate the output realizations of $OIP$, based on this equation.\n",
    "\n",
    "\\begin{equation}\n",
    "OIP^\\ell = V^\\ell \\cdot \\overline{\\phi}^\\ell \\cdot \\overline{s_o}^\\ell \\cdot 6.29 \\quad \\forall \\, \\ell = 1,\\ldots, L\n",
    "\\end{equation}\n",
    "\n",
    "where 6.26 $bbls/m^3$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe4AAAGbCAYAAAAC+HMkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de7xcdX3v/9dbIiLIRSUoScCgRSx6FGik9vpT0RY9VfTUC9ZWtFRq6+VQPafej/yKbbW1WlsvLSoFtIJILdDWWpGqtB4Bw0W5VUFASIIQFcQAAiGf88daW4bNzt6TZM/M/iav5+OxH5lZ6ztrfb4zSd77+11r1kpVIUmS2vCASRcgSZKGZ3BLktQQg1uSpIYY3JIkNcTgliSpIQa3JEkNMbilBS7J3knWJdlu0rWMQ5Jrkzyjf/yWJB8dwT7+Jsnb53u70jgY3Nrq9UFwRx9+301yQpKHTLquYVXVdVX1kKq6Z1Nel+TlSS5Jcnvf7w8n2W1g/TFJPjHwvJLc1r9Pq5O8d2O/LGxK2y1RVX9SVb+zJdvo34f/nLbdV1XVsVtWnTQZBre2Fc+pqocABwAHAm8e586TLBrRdpPkfv+Ok7wBeDfwv4FdgacAjwLOSrL9LJt8Uv8+HQL8BvDKIdr+f8CLgd/evF5I2hQGt7YpVfVd4N/oAhyAJA9K8p4k1yW5sZ9GfXC/bvck/5zkliQ/SPIfU0GZZEmSf0iyNsk1SV43sM1jkpyW5BNJbgXe0o/6HzbQ5sAk30vywCQPSPK2JN9JclOSk5Ls2rdb3o9wF/XPv5Tkj5N8BbgdePRgH5PsAvz/wGur6nNVdXdVXQu8iC68f3OI9+m/gP8AnjBE26uAr0x7T3dN8rEkN/Qj8ndOjciTPCbJvyf5ft//vx+cCZjWl5/MCiT5QD/Cn/pZn+SYft2bknw7yY+SXJ7k+f3ynwb+Bvi5/jW39MtPSPLOgf28MslV/Wd8ZpIlA+sqyauSXJnk5iQfTJJ+3U8l+XKSH/Z9+dRc75e0pQxubVOSLAOeBVw1sPjdwGPpguengKXA/+nXvQFYBSwGHgG8Bag+vP8J+Hrf/hDg6CS/OrDdw4DTgN2APwe+Cvz6wPrfAE6rqruBl/c/T6ML4ocAH5ilK78FHAXsDHxn2rqfB3YAPjO4sKrWAf8KPHOW7QKQZH/gl4CLhmj7uL7t4Ht6IrCe7v08EPgVYGrKO8CfAkuAnwb2Ao6Zaz9V9Zr+kMFDgF8EbgbO6Fd/u69hV7pfWj6RZM+qugJ4FfDV/rX3+wUhydP7el4E7En3fp4yrdmvAU8GntS3m/qcjwU+DzwUWAb89Vz9kLaUwa1txelJfgRcD9wEvAO6qWa66eA/qKofVNWPgD8BDu9fdzfdf+aP6keu/1HdBf6fDCyuqj+qqruq6mrgIwOvgy4sTq+qDVV1B/BJ4CUD+z28XwbwUuC9VXV1H7BvBg6fZYr9hKq6rKrW98E/aHfge1W1fobX3dCv35gLk9xM90vJR4G/m6PtbcAVwJeAD/V9ewTdL0dHV9VtVXUT8L6+v1TVVVV1VlXdWVVrgffSTbcPJcli4HS6GYWL+m1+uqrW9O/1p4ArgYOH3ORLgeOr6sKqupPuvf+5JMsH2ryrqm6pquuAL3Lv7MLddLMYS6rqx1V1n2Pp0igY3NpWPK+qdgaeCjyOe8NrMbAjcEE/HX4L8Ll+OXQj5auAzye5Osmb+uWPApZMvaZ/3VvoRuVTrp9Ww2l0gbAE+GWg6KajoRt9Do6cvwMsmra9QdO3Peh7wO4bCf09+/Ubc1BVPbSqHlNVb6uqDbO1pZsZeDHws8BO/fJHAQ8Ebhh4b/4W2AMgyR5JTumn0G8FPsHsv0z8RJIH0r2Pn6yqUwaWvyzJxQP7e8Kw22Tae9//4vR9upmUKd8deHx732+AP6SbQTg/yWVJPM6vkTO4tU2pqi8DJwDv6Rd9D7gDeHxV7db/7NpPx1JVP6qqN1TVo4HnAK9PcghdcF4z8Jrdqmrnqnr24O6m7fsWumnVF9FNk59c996ebw1d4E3Zm26q+caNdWWWbn4VuBP4H4MLk+xENxI+e5bXbpLqnNrvc+rwwvX9/ncfeG92qarH9+v/tK//iVW1C90x9wy5y78GfgS8bWpBkkfRzXa8Bnh4Px1+6cA257oF4n3e+/59ejiweq5iquq7VfXKqloC/C7woSQ/NWRfpM1icGtb9JfAM5Mc0I8oPwK8L8nUiHDp1LHqJL/Wn4AU4Fbgnv7nfODWJG9M8uAk2yV5QpInz7HvTwIvozvW/cmB5ScDf5Bkn3RfVfsT4FMbme6eVVX9kO44718nObQ/+W058Gm64/Uf39RtDuFdwFFJHllVN9D9gvIXSXbpT7x7TJKp6fCdgXXALUmW0p35Pqckv0s3pf4b02YCdqIL57V9u1dw35PqbgSWZeNn038SeEWSA5I8iO69P68/oW+uml7YnzcB3TH3ovv7IY2Mwa1tTn9c9SRg6gIcb6SbDj+3n7r9ArBfv27f/vk6ulHlh6rqS/13qp9Dd6zzGrqR+0fpTo6azZn9Nm+sqq8PLD+eLlDP6bf3Y+C1W9DHP6Obun8P3S8c59GNhA/pj+POq6q6BPgy94bwy4DtgcvpAu00uml66H6pOAj4IfAvTDuJbhYvoTtxb83AmeVvqarLgb+g+3xuBP4b3VnuU/4duAz4bpL7HSaoqrPp/i78A905AI/hvucqzObJwHlJ1tF9tv+zqq4Z8rXSZsm9M3WSJGmhc8QtSVJDDG5JkhpicEuS1BCDW5Kkhozkxgfjsvvuu9fy5csnXYYkSZvkggsu+F5VLZ675f01HdzLly9n5cqVky5DkqRNkmT6PQaG5lS5JEkNMbglSWqIwS1JUkMMbkmSGmJwS5LUEINbkqSGGNySJDXE4JYkqSEGtyRJDTG4JUlqiMEtSVJDDG5JkhpicEuS1BCDW5KkhhjckiQ1xOCWJKkhBrckSQ0xuCVJaojBLUlSQwxuSZIaYnBLktQQg1uSpIYY3JIkNcTgliSpIQa3JEkNGVlwJzk+yU1JLp22/LVJvpnksiR/NrD8zUmu6tf96qjqkiSpZYtGuO0TgA8AJ00tSPI04DDgiVV1Z5I9+uX7A4cDjweWAF9I8tiqumeE9UmS1JyRjbir6hzgB9MW/x7wrqq6s29zU7/8MOCUqrqzqq4BrgIOHlVtkiS1atzHuB8L/FKS85J8OcmT++VLgesH2q3ql91PkqOSrEyycu3atSMuV5KkhWXcwb0IeCjwFOB/A6cmCZAZ2tZMG6iq46pqRVWtWLx48egqlSRpARp3cK8CPlOd84ENwO798r0G2i0D1oy5NkmSFrxxB/fpwNMBkjwW2B74HnAmcHiSByXZB9gXOH/MtUmStOCN7KzyJCcDTwV2T7IKeAdwPHB8/xWxu4AjqqqAy5KcClwOrAde7RnlkiTdX7rcbNOKFStq5cqVky5DkqRNkuSCqlqxOa/1ymmSJDXE4JYkqSEGtyRJDTG4JUlqiMEtSVJDDG5JkhpicEuS1BCDW5KkhhjckiQ1xOCWJKkhBrckSQ0xuCVJaojBLUlSQwxuSZIaYnBLktQQg1uSpIYY3JIkNcTgliSpIQa3JEkNMbglSWqIwS1JUkMMbkmSGmJwS5LUEINbkqSGGNySJDXE4JYkqSEGtyRJDTG4JUlqiMEtSVJDDG5JkhoysuBOcnySm5JcOsO6/5WkkuzeP0+Sv0pyVZJvJDloVHVJktSyUY64TwAOnb4wyV7AM4HrBhY/C9i3/zkK+PAI65IkqVkjC+6qOgf4wQyr3gf8IVADyw4DTqrOucBuSfYcVW2SJLVqrMe4kzwXWF1VX5+2ailw/cDzVf0ySZI0YNG4dpRkR+CtwK/MtHqGZTXDMpIcRTedzt577z1v9UmS1IJxjrgfA+wDfD3JtcAy4MIkj6QbYe810HYZsGamjVTVcVW1oqpWLF68eMQlS5K0sIwtuKvqkqrao6qWV9VyurA+qKq+C5wJvKw/u/wpwA+r6oZx1SZJUitG+XWwk4GvAvslWZXkyFmafxa4GrgK+Ajw+6OqS5Kklo3sGHdVvWSO9csHHhfw6lHVIknS1sIrp0mS1BCDW5KkhhjckiQ1xOCWJKkhBrckSQ0xuCVJaojBLUlSQwxuSZIaYnBLktQQg1uSpIYY3JIkNcTgliSpIQa3JEkNMbglSWqIwS1JUkMMbkmSGmJwS5LUEINbkqSGGNySJDXE4JYkqSEGtyRJDTG4JUlqiMEtSVJDDG5JkhpicEuS1BCDW5KkhhjckiQ1xOCWJKkhBrckSQ0xuCVJaojBLUlSQ0YW3EmOT3JTkksHlv15kv9K8o0k/5hkt4F1b05yVZJvJvnVUdUlSVLLRjniPgE4dNqys4AnVNUTgW8BbwZIsj9wOPD4/jUfSrLdCGuTJKlJIwvuqjoH+MG0ZZ+vqvX903OBZf3jw4BTqurOqroGuAo4eFS1SZLUqkke4/5t4F/7x0uB6wfWreqX3U+So5KsTLJy7dq1Iy5RkqSFZSLBneStwHrg76cWzdCsZnptVR1XVSuqasXixYtHVaIkSQvSonHvMMkRwK8Bh1TVVDivAvYaaLYMWDPu2iRJWujGOuJOcijwRuC5VXX7wKozgcOTPCjJPsC+wPnjrE2SpBaMbMSd5GTgqcDuSVYB76A7i/xBwFlJAM6tqldV1WVJTgUup5tCf3VV3TOq2iRJalXuna1uz4oVK2rlypWTLkOSpE2S5IKqWrE5r/XKaZIkNcTgliSpIQa3JEkNMbglSWqIwS1JUkMMbkmSGmJwS5LUEINbkqSGjP1a5duqgw/Yn1WrV8/aZtnSpZx/8eVjqkiS1CKDe0xWrV7NmmN3mrXNkrfPHuySJDlVLklSQwxuSZIaYnBLktQQg1uSpIYY3JIkNcTgliSpIQa3JEkNMbglSWqIwS1JUkMMbkmSGmJwS5LUEINbkqSGGNySJDXE4JYkqSEGtyRJDTG4JUlqiMEtSVJDDG5JkhpicEuS1BCDW5KkhowsuJMcn+SmJJcOLHtYkrOSXNn/+dB+eZL8VZKrknwjyUGjqkuSpJaNcsR9AnDotGVvAs6uqn2Bs/vnAM8C9u1/jgI+PMK6JElq1siCu6rOAX4wbfFhwIn94xOB5w0sP6k65wK7JdlzVLVJktSqcR/jfkRV3QDQ/7lHv3wpcP1Au1X9MkmSNGDRMI2SPKGqLp275WbLDMtqI7UcRTedzt577z3CksbvjtvWsWTxrrO2WbZ0KedffPmYKpIkLTRDBTfwN0m2pztu/cmqumUz93djkj2r6oZ+KvymfvkqYK+BdsuANTNtoKqOA44DWLFixYzh3qr1Gzaw5tidZm2z5O2rx1SNJGkhGmqqvKp+EXgpXbiuTPLJJM/cjP2dCRzRPz4COGNg+cv6s8ufAvxwakpdkiTda9gRN1V1ZZK3ASuBvwIOTBLgLVX1mentk5wMPBXYPckq4B3Au4BTkxwJXAe8sG/+WeDZwFXA7cArNrtHW7m5ptOdSpekrduwx7ifSBem/x04C3hOVV2YZAnwVeB+wV1VL9nI5g6ZoW0Brx626G3ZXNPpTqVL0tZt2BH3B4CP0I2u75haWFVr+lG4JEkag2GD+9nAHVV1D0CSBwA7VNXtVfXxkVUnSZLuY9jvcX8BePDA8x37ZZIkaYyGDe4dqmrd1JP+8Y6jKUmSJG3MsMF92+CNP5L8DHDHLO0lSdIIDHuM+2jg00mmLoqyJ/Di0ZQkSZI2ZqjgrqqvJXkcsB/d5Un/q6ruHmllkiTpfoa+AAvwZGB5/5oDk1BVJ42kKkmSNKNhL8DyceAxwMXAPf3iAgxuSZLGaNgR9wpg//4KZ5IkaUKGDe5LgUcC3vhjIw4+YH9Wrd745UZvW7cOmP3OX5IkzWXY4N4duDzJ+cCdUwur6rkjqapBq1avnvUa4jsffesYq5Ekba2GDe5jRlmEJEkazrBfB/tykkcB+1bVF5LsCGw32tIkSdJ0Q105LckrgdOAv+0XLQVOH1VRkiRpZsNe8vTVwC8AtwJU1ZXAHqMqSpIkzWzY4L6zqu6aepJkEd33uCVJ0hgNG9xfTvIW4MFJngl8Gvin0ZUlSZJmMmxwvwlYC1wC/C7wWeBtoypKkiTNbNizyjcAH+l/JEnShAx7rfJrmOGYdlU9et4rkiRJG7Up1yqfsgPwQuBh81+OJEmazVDHuKvq+wM/q6vqL4Gnj7g2SZI0zbBT5QcNPH0A3Qh855FUJEmSNmrYqfK/GHi8HrgWeNG8VyNJkmY17FnlTxt1IZIkaW7DTpW/frb1VfXe+SlHW+qO29axZPGus7ZZtnQp5198+ZgqkiTNp005q/zJwJn98+cA5wDXj6Iobb71GzbMel9wgCVvXz2maiRJ823Y4N4dOKiqfgSQ5Bjg01X1O6MqTJIk3d+wlzzdG7hr4PldwPJ5r0aSJM1q2BH3x4Hzk/wj3RXUng+cNLKqJEnSjIa9AMsfA68AbgZuAV5RVX+yuTtN8gdJLktyaZKTk+yQZJ8k5yW5Msmnkmy/uduXJGlrNexUOcCOwK1V9X5gVZJ9NmeHSZYCrwNWVNUTgO2Aw4F3A++rqn3pfkE4cnO2L0nS1myo4E7yDuCNwJv7RQ8EPrEF+11Ed2/vRXS/ENxAdwnV0/r1JwLP24LtS5K0VRp2xP184LnAbQBVtYbNvORpVa0G3gNcRxfYPwQuAG6pqvV9s1XA0plen+SoJCuTrFy7du3mlCBJUrOGDe67qqrob+2ZZPYvCs8iyUOBw4B9gCXATsCzZmh6v9uIAlTVcVW1oqpWLF68eHPLkCSpScMG96lJ/hbYLckrgS8AH9nMfT4DuKaq1lbV3cBngJ/vtz11lvsyYM1mbl+SpK3WsNcqf0+SZwK3AvsB/6eqztrMfV4HPCXJjsAdwCHASuCLwAuAU4AjgDM2c/uSJG215gzuJNsB/1ZVzwA2N6x/oqrOS3IacCHdncYuAo4D/gU4Jck7+2Uf29J9SZK0tZkzuKvqniS3J9m1qn44HzutqncA75i2+Grg4PnYviRJW6thr5z2Y+CSJGfRn1kOUFWvG0lVkiRpRsMG97/0P5IkaYJmDe4ke1fVdVV14rgKkiRJGzfX18FOn3qQ5B9GXIskSZrDXMGdgcePHmUhkiRpbnMFd23ksSRJmoC5Tk57UpJb6UbeD+4f0z+vqtplpNVJkqT7mDW4q2q7cRUiSZLmtin345YkSRNmcEuS1JBhL8Circgdt61jyeJdZ22zbOlSzr/48jFVJEkalsG9DVq/YQNrjp39lupL3r56TNVIkjaFwT2Egw/Yn1WrZw+y29atA2YPQ0mStpTBPYRVq1fPOULd+ehbZ10vSdJ88OQ0SZIaYnBLktQQg1uSpIYY3JIkNcTgliSpIQa3JEkNMbglSWqIwS1JUkMMbkmSGmJwS5LUEINbkqSGGNySJDXE4JYkqSEGtyRJDTG4JUlqiMEtSVJDJhLcSXZLclqS/0pyRZKfS/KwJGclubL/86GTqE2SpIVsUiPu9wOfq6rHAU8CrgDeBJxdVfsCZ/fPJUnSgLEHd5JdgF8GPgZQVXdV1S3AYcCJfbMTgeeNuzZJkha6SYy4Hw2sBf4uyUVJPppkJ+ARVXUDQP/nHhOoTZKkBW0Swb0IOAj4cFUdCNzGJkyLJzkqycokK9euXTuqGiVJWpAmEdyrgFVVdV7//DS6IL8xyZ4A/Z83zfTiqjquqlZU1YrFixePpWBJkhaKsQd3VX0XuD7Jfv2iQ4DLgTOBI/plRwBnjLs2SZIWukUT2u9rgb9Psj1wNfAKul8iTk1yJHAd8MIJ1SZJ0oI1keCuqouBFTOsOmTctUiS1BKvnCZJUkMMbkmSGmJwS5LUEINbkqSGGNySJDXE4JYkqSEGtyRJDTG4JUlqiMEtSVJDDG5JkhpicEuS1BCDW5KkhhjckiQ1xOCWJKkhBrckSQ0xuCVJasiiSReghemO29axZPGuG11/910/5oHb7zDrNpYtXcr5F18+36VJ0jbN4NaM1m/YwJpjd9ro+p2PvpW17374rNtY8vbV812WJG3znCqXJKkhBrckSQ0xuCVJaojBLUlSQwxuSZIaYnBLktQQg1uSpIYY3JIkNcTgliSpIQa3JEkNMbglSWqIwS1JUkMMbkmSGjKx4E6yXZKLkvxz/3yfJOcluTLJp5JsP6naJElaqCY54v6fwBUDz98NvK+q9gVuBo6cSFWSJC1gEwnuJMuA/w58tH8e4OnAaX2TE4HnTaI2SZIWskmNuP8S+ENgQ//84cAtVbW+f74KWDrTC5MclWRlkpVr164dfaWSJC0gYw/uJL8G3FRVFwwunqFpzfT6qjquqlZU1YrFixePpEZJkhaqRRPY5y8Az03ybGAHYBe6EfhuSRb1o+5lwJoJ1CZJ0oI29hF3Vb25qpZV1XLgcODfq+qlwBeBF/TNjgDOGHdtkiQtdJMYcW/MG4FTkrwTuAj42Dh2evAB+7Nq9epZ29y2bh2w0zjK2arccds6lizeddY2y5Yu5fyLLx9TRZLUvokGd1V9CfhS//hq4OBx17Bq9WrWHDt7KO989K1jqmbrsn7Dhjnf2yVvn/2XJknSfXnlNEmSGmJwS5LUEINbkqSGGNySJDXE4JYkqSEGtyRJDTG4JUlqiMEtSVJDDG5JkhpicEuS1BCDW5KkhhjckiQ1xOCWJKkhBrckSQ0xuCVJaojBLUlSQwxuSZIaYnBLktQQg1uSpIYY3JIkNWTRpAuQZnPwAfuzavXqWdssW7qU8y++fEwVSdJkGdxa0FatXs2aY3eatc2St88e7JK0NXGqXJKkhhjckiQ1xOCWJKkhBrckSQ0xuCVJaojBLUlSQwxuSZIaYnBLktSQsV+AJclewEnAI4ENwHFV9f4kDwM+BSwHrgVeVFU3j7s+bZu8QpukVkziymnrgTdU1YVJdgYuSHIW8HLg7Kp6V5I3AW8C3jiB+rQN8gptklox9qnyqrqhqi7sH/8IuAJYChwGnNg3OxF43rhrkyRpoZvoMe4ky4EDgfOAR1TVDdCFO7DH5CqTJGlhmthNRpI8BPgH4OiqujXJsK87CjgKYO+99x5dgRqLO25bx5LFu250/W3r1gGzT2HPtQ2PTUvamkwkuJM8kC60/76qPtMvvjHJnlV1Q5I9gZtmem1VHQccB7BixYoaS8EamfUbNsx6bHnno2/d4m14bFrS1mTsU+XphtYfA66oqvcOrDoTOKJ/fARwxrhrkyRpoZvEiPsXgN8CLklycb/sLcC7gFOTHAlcB7xwArVJkrSgjT24q+o/gY0d0D5knLVIktSaiZ2cJo3LXCevwXAnwUnSQmBwa6s318lrMNxJcJK0EBjc0pD82pmkhcDglobk184kLQTeHUySpIYY3JIkNcTgliSpIQa3JEkN2epPTjv4gP1ZtXrjJw35/V3Nl2G+L+6Z55K21FYf3KtWr97im1hIwxjm++KeeS5pSzlVLklSQwxuSZIaYnBLktQQg1uSpIZs9SenSa2Z65sQnpkubdsMbmmBmeubEJ6ZLm3bnCqXJKkhBrckSQ0xuCVJaojBLUlSQwxuSZIaYnBLktQQvw4mjdEwdxDzjnWSZmNwS2M0zB3EvGOdpNkY3FJj5uu+33Ndoe3uu37MA7ffYYv3I2l+GdxSY+brvt/D3Kt+7bsfvsX7kTS/DG5pK7S1HUv3+u3SvQxuaSu0tR1L9/rt0r0MbknqObJXCwxuSeo5slcLDG5Jm22uY+kLaYQ612ga2jrur23XggvuJIcC7we2Az5aVe+acEmSNmKuY+kLaYQ612ga2jrur23XggruJNsBHwSeCawCvpbkzKpaGL+yS5p3joQXNo/7LzwLKriBg4GrqupqgCSnAIcB/q2QtlKOhBc2j/svPKmqSdfwE0leABxaVb/TP/8t4Ger6jUDbY4CjuqfPgG4dOyFjt7uwPcmXcQI2K+22K+22K+27FdVO2/OCxfaiDszLLvPbxZVdRxwHECSlVW1YhyFjZP9aov9aov9asvW3K/Nfe1Cu63nKmCvgefLgDUTqkWSpAVnoQX314B9k+yTZHvgcODMCdckSdKCsaCmyqtqfZLXAP9G93Ww46vqsllectx4Khs7+9UW+9UW+9UW+zXNgjo5TZIkzW6hTZVLkqRZGNySJDWkieBOcmiSbya5KsmbZlj/oCSf6tefl2T5+KvcdEP065eTXJhkff8d9yYM0a/XJ7k8yTeSnJ3kUZOoc1MN0a9XJbkkycVJ/jPJ/pOoc1PN1a+Bdi9IUkma+GrOEJ/Xy5Os7T+vi5P8ziTq3FTDfF5JXtT/G7ssySfHXePmGOLzet/AZ/WtJLdMos5NNUS/9k7yxSQX9f8nPnvOjVbVgv6hO0nt28Cjge2BrwP7T2vz+8Df9I8PBz416brnqV/LgScCJwEvmHTN89ivpwE79o9/byv6vHYZePxc4HOTrns++tW32xk4BzgXWDHpuufp83o58IFJ1zqCfu0LXAQ8tH++x6Trno9+TWv/WrqTlyde+zx8XscBv9c/3h+4dq7ttjDi/sllUKvqLmDqMqiDDgNO7B+fBhySZKaLuSwkc/arqq6tqm8AGyZR4GYapl9frKrb+6fn0n1ff6Ebpl+D1+XciWkXD1qghvn3BXAs8GfAj8dZ3BYYtl+tGaZfrwQ+WFU3A1TVTWOucXNs6uf1EuDksVS2ZYbpVwG79I93ZYhrl7QQ3EuB6weer+qXzdimqtYDPwQePpbqNt8w/WrRpvbrSOBfR1rR/BiqX0leneTbdCH3ujHVtiXm7FeSA4G9quqfx1nYFhr27+Gv99OTpyXZa4b1C80w/Xos8NgkX0lybn/HxYVu6P83+kNr+wD/Poa6ttQw/ToG+M0kq4DP0s0mzKqF4J7zMqhDtlloWqx5GEP3K8lvAiuAPx9pRfNjqH5V1Qer6jHAG4G3jbyqLTdrv5I8AHgf8IaxVTQ/hvm8/glYXlVPBL7AvbN2C9kw/VpEN13+VLqR6UeT7DbiurbUpvx/eDhwWlXdM8J65ssw/XoJcEJVLQOeDXy8/3e3US0E968QeW0AAAZzSURBVDCXQf1JmySL6KYbfjCW6jbf1np516H6leQZwFuB51bVnWOqbUts6ud1CvC8kVY0P+bq1850N/P5UpJrgacAZzZwgtqcn1dVfX/g795HgJ8ZU21bYtj/D8+oqrur6hrgm3RBvpBtyr+vw2ljmhyG69eRwKkAVfVVYAe6G6ts3KQP3g9xcH8RcDXd1MjUwf3HT2vzau57ctqpk657Pvo10PYE2jk5bZjP60C6Ezb2nXS989yvfQcePwdYOem656Nf09p/iTZOThvm89pz4PHzgXMnXfc89etQ4MT+8e50U7UPn3TtW9qvvt1+wLX0Fw9b6D9Dfl7/Cry8f/zTdME+a/8m3rEhO/9s4Fv9f/Zv7Zf9Ed1oDbrfUD4NXAWcDzx60jXPU7+eTPcb223A94HLJl3zPPXrC8CNwMX9z5mTrnme+vV+4LK+T1+cLQAX0s9c/ZrWtongHvLz+tP+8/p6/3k9btI1z1O/ArwXuBy4BDh80jXPR7/658cA75p0rfP8ee0PfKX/e3gx8CtzbdNLnkqS1JAWjnFLkqSewS1JUkMMbkmSGmJwS5LUEINbkqSGGNySJDXE4JYkqSGLJl2AJMFProt+LN2dklZWVQvXDpfGzhG3tImS3JPk4iSXJvmnhXYDhyT/d8h2y5KckeTKJN9O8v4k2w+sXzfweLDPn06y4wzbe3CSLyfZrn9+fJKbklw6rd2hSb6Z5KokbxpYdRjdnZPuprti4PTtL5++rdmWz9SPGdZtn+Sc/h4HUhMMbmnT3VFVB1TVE+huZvPqUewknU3+N1pVPz/Xtvr71X8GOL2q9qW7FeRDgD/eyGYH+3wX8KoZ2vw28Jm6965NJ9BdN3uwju2ADwLPorvU40uS7N+v3g/4alW9Hvi9Yfq6paq7R/LZwIvHsT9pPhjc0pb5Kv39dZP8ZpLz+5Hp3w6MPHdK8i9Jvt6PWF+8sfb96PGKJB8CLgQ+luT3p3aW5Jgkb+gfv77f3qVJjh5os67/c/q2Bu9S9HTgx1X1dwB92P4B8Nszjaan+Q/gp2ZY/lLgjKknVXUO979L38HAVVV1dR+ap9CNtKEbZd/cP97YLRsXJTlx4B7aO86xHNj4Z9A7va9daoLBLW2mPpgPobvN5U/Tjdp+oaoOoAueqTA4FFhTVU/qR6yfm6P9fsBJVXUg8FfcdzT4IuDTSX4GeAXws3S32nxlkgNnKPMn26qq7wwsfzxwwWDDqroVuI6ZQ3mqz4voRsuXTFu+Pd3Nfa7d2Gt7S+nuVjVlVb8MuhmAX03y18A5G3n9fsBx1d1D+1bg9+dYPuV+n8HAukvpbugjNcHgljbdg5NcTHfHtocBZ9EF+M8AX+vXHQI8um9/CfCMJO9O8ktV9cM52n+nqs4FqKqLgD2SLEnyJODmqroO+EXgH6vqtqpaRxd6vzRDrT/Z1jQBZrrD0MaWT/V5JV24f2za+t2BW2Z43Uzbn64Aqur2qjqyql5bVR/cyOuvr6qv9I8/Qfc+zLZ8ykyfAf1+7wHuSrLzEPVLE+cJGdKmu6OqDkiyK/DPdMe4i+4eyG+e3riqvtWPkJ8N/GmSz9NNCd+vfZLldLdxHXQa8ALgkXRTyzBzAM5k+ramXAb8+rR970I3nf7tGdrf0c8MbMwddLfXncsq7jtlv4zu/sPDmv5LRc2xvHsyw2dQVX800ORBwI83oQ5pYhxxS5upH7W9DvhfdFO7L0iyB0CShyV5VP94CXB7VX0CeA9wEN0JUTO2n8EpwOF04X1av+wc4HlJdkyyE/B8umPPwzob2DHJy/r9bwf8BXBCVd2+CdsBoKpuBrZLMld4fw3YN8k+/fT64cCZm7CrvZP8XP/4JcB/zrEc2OhnMLXu4cDaqrp7E+qQJsbglrZAP5X9deCJwNuAzyf5Bt30+Z59s/8GnN9PNb8VeGdVXT5L++n7uAzYGVhdVTf0yy6kO2v7fOA84KN9LcPWXXRh/8IkVwLfohtxvmX43t/P5xmYok5yMt3Je/slWZXkyKpaD7wG+DfgCuDUvn/DugI4on/PHgZ8eI7lU+73GQysexrw2U2oQZqodP9+JWnL9CfHvb6qfmvStWyKJJ8B3lxV35x0LdIwHHFLmhf9iP+LU1+Da0E/XX+6oa2WOOKWJKkhjrglSWqIwS1JUkMMbkmSGmJwS5LUEINbkqSGGNySJDXE4JYkqSH/D1NwItkfSFWMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "OIP = apor * vol * so * 6.29\n",
    "\n",
    "plt.subplot(111)\n",
    "GSLIB.hist_st(OIP,0.0,0.8,log=False,cumul=False,bins=50,weights=None,xlabel=\"Reservoir OIP ($10^6$ bbls)\",title=\"Reservoir OIP Realizations\")\n",
    "plt.ylim(0.0,175)\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.2, wspace=0.2, hspace=0.2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's not a bad idea to check your calculation. Let's compare the first and last realizations of $OIP$ to the results from performing the calculation by hand. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Check realization 0, broadcast method = 0.066 compare to by-hand 0.066\n",
      "Check realization L-1, broadcast method = 0.053 compare to by-hand 0.053\n"
     ]
    }
   ],
   "source": [
    "print('Check realization 0, broadcast method = ' + str(round(OIP[0],3)) + ' compare to by-hand ' + str(round(apor[0]*vol[0]*so[0]*6.29,3)))\n",
    "print('Check realization L-1, broadcast method = ' + str(round(OIP[L-1],3)) + ' compare to by-hand ' + str(round(apor[L-1]*vol[L-1]*so[L-1]*6.29,3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results check out.  \n",
    "\n",
    "Now let's look at what happens if too few realizations are used!  Set the $L$ smaller, calculate 3 sets of $L$ realizations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = 100                                                # set L small \n",
    "# First set of realizations\n",
    "apor = np.random.normal(apor_mean, apor_stdev, size=L) # average porosity MCS simulation L times and store in array \n",
    "vol = np.random.lognormal(vol_mu, vol_sigma, size=L)/1.0e6 # volume ...\n",
    "so = np.random.uniform(so_min, so_max, size=L)         # saturation oil\n",
    "OIP1 = apor * vol * so * 6.29                          # calculate OIP\n",
    "\n",
    "# Second set of realizations\n",
    "apor = np.random.normal(apor_mean, apor_stdev, size=L) # average porosity MCS simulation L times and store in array \n",
    "vol = np.random.lognormal(vol_mu, vol_sigma, size=L)/1.0e6 # volume ...\n",
    "so = np.random.uniform(so_min, so_max, size=L)         # saturation oil\n",
    "OIP2 = apor * vol * so * 6.29                          # calculate OIP\n",
    "\n",
    "# Third set of realizations\n",
    "apor = np.random.normal(apor_mean, apor_stdev, size=L) # average porosity MCS simulation L times and store in array \n",
    "vol = np.random.lognormal(vol_mu, vol_sigma, size=L)/1.0e6 # volume ...\n",
    "so = np.random.uniform(so_min, so_max, size=L)         # saturation oil\n",
    "OIP3 = apor * vol * so * 6.29                          # calculate OIP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the distributions, and calculate and compare the average for each of the 3 sets.\n",
    "\n",
    "* note that I use the subplot method below that allows me to share the y axis of all the plots to improve our ability to compare the distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5AAAAGWCAYAAAAOimwQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3debglZ1Uv/u8iCWOAENJg6JCEITIEJUATRZAZjFwDeAUERUGQ4ADKT39eCJer4BRREAdQCcJNBJnVC8EpYQgRZbCBMOMFgUCaTicMAQIRSLLuH1VNdk76dNdJn33O2d2fz/Psp2tXvbtqvXtYfVZNb3V3AAAAYE+utd4BAAAAsBgUkAAAAEyigAQAAGASBSQAAACTKCABAACYRAEJAADAJArIdVJVR1bVJVV1wHrHshaq6jNV9cBx+plV9Zdz2MZfVNX/Wu31LqKq6qq67XrHwb5LDpPD5kkOY97kMDlsnvb1HLawBeT4Q7h0/PFfUFWnVdXB6x3XVN392e4+uLsvX8nrqurxVfWhqvrG2O8/r6pDZpY/u6peMfO8q+rr4/u0rar+cLlkuZK2e6O7f7e7f3Zv1jG+D+9Yst6f6+7f2rvo9iqma1fV86vq/PE9/HRVvWC94rkm9qYPS797y7R5SlVtrapvVtVpqxL0gpLD5DA5bPXNM4dV1XWq6qVVdV5Vfa2q3l9VP7x60S8WOUwOk8NW3xr8HfaKqtpeVV+tqv9bVdfoe7CwBeToxO4+OMlxSe6S5OS13HhVHTin9VZVXe2zqapfTfLcJL+W5MZJvj/JUUnOqqpr72aVdx7fpwck+YkkT5rQ9j5JfjzJE65ZL/ZLJyfZkuT4JDdMcr8k71/XiFZu3n34fJLfTvKyVVznIpPD5LCNRA7bvQOTfC7Dd+vGSf5XktdW1dGrtP5FJIfJYRuJHLZnpyQ5urtvlOShSX67qu624rV090I+knwmyQNnnv9+kr+feX6dJM9L8tkkO5L8RZLrjcsOS/KmJBcn+VKSf0lyrXHZLZL8TZKLknw6yS/NrPPZSV6f5BVJvprk15NcmuTQmTZ3SfKFJAdlKNCfleS8JBcm+askNx7bHZ2kkxw4Pj87ye8k+ddxnbdd0t8bJbkkyaOWzD94XPcTZmJ8xczynl1XktcleeEy7+nStq9N8qKZ5zdO8tIk25Nsy1AIHDAuu02Styb54tj/v05yyK4+r9kYk7xw7NfOx2VJnj0ue0aS/0zytSQfTfKj4/w7JPmvJJePr7l4nH9akt+e2eaTknxy/IzfmOQWS/r6c0k+keTLSV6UpMZlt03y9iRfGfvymonfyTcledpulu+yP+Oyx4+f/QsyfC8/leQHxvmfGz/jx820Py3Dd/qscX1vT3LUrj7L7Oa3cA36sMvfR5ITknwrybfHz+QDe3ivfjvJaeudR9bzETlMDpPDFjaHzazvg0l+bL3zyXo8IofJYXLYQuewJLcbv0uP2lPbpY9FPwKZJKmqI5L8cIYv6U7PTfLdGfaK3TbJ5gyJJkl+Ncn5STYluXmSZybpcW/TGUk+MLZ/QJKnVdUPzaz3YRmS1yFJ/iDJO5P82Mzyn0jy+u7+doYv3eMz7D24dYYk88LddOWnkpyUYY/DeUuW/UCS6yb529mZ3X1Jkn9M8qDdrDdJUlV3TPKDmbAno6puP7adfU9Pz5BYbpshQT84yc5D35Vhr8YtMiSWW2ZIULvV3U/p4RSSg5PcK0MSecO4+D/HGG6c5DlJXlFVh3f3xzIknXeOrz1k6Xqr6v5jPI9KcniG9/PVS5r9SJK7J7nz2G7n5/xbSc5McpMkRyT50z31Y/SuJL9SVb9QVd9TVbVk+S77M7P8+zL8MXLTJK8c4717hvf7sUleuOT0oJ8cYz0sybkZ/rPYld39Fib3YXe/j+7+pyS/myHJH9zdd15m/eyCHCaH7SJ+OexKGzKHVdXNx7g+sqe2+zo5TA7bRfxy2JU2VA6rqj+rqm8k+XiGAvIflmu7rJVWnBvlkWFPyiUZqv5O8paMe1oy/Ii+nuQ2M+3vkeTT4/RvZvhxLN279H1JPrtk3slJ/ndfucfmnCXLfzbJW2e2+7kk9x6fvyXJLyyp9L+d4TSYo3P1PV+/uZv+PjbJBcss+70kZ83EuHTP11czJIT/zLC36lrLrGdn26+P069Kcp1x2c2TfDMze0ySPCbJ25ZZ18OTvH/J53W1PV8zyzeNbR69m/fg3CQPG6cfn+QdS5aflnHPV4Y9dL8/s+zg8b0/eqav95pZ/tokzxin/yrJqUmOWOF38oAkv5hhD9Y3M5yu+bgV9OcTM8u+Z4zx5jPzvpjkuJm+vnpJ/y5PcsuZ/t02e/gtrKQPmfb7eMVy/V3yOkcg5bDZZXJYy2FZrBx2UJI3J3nxSt7jfekROWx2mRzWclgWK4cdkGGHwbOSHLSS97l78Y9APry7b5jkvklun2EPQDL8CK6f5L1VdXFVXZzkn8b5ybDH6pNJzqyqT1XVM8b5RyW5xc7XjK97ZoYf7U6fWxLD65Pco6pukeTeGb4w/zIuu0WuugfrvAxJ6+bZtaXrnvWFJIctc77/4ePy5dy1u2/S3bfp7md19xW7a5vhR/DjGb6oNxjnH5XhP8ztM+/Ni5PcLEmq6mZV9erxou+vZji95LCrrX0XquqgDO/jK7v71TPzf7qqzp3Z3p2mrjNL3vse9hB+McMem50umJn+xtjvJPkfGX7w76mqj1TVpOsPuvvy7n5Rd98zw57R30nysqq6w8T+7JiZvnRc59J5s3u+vvN9Gfv3pbHfs/b0W1hJH6b8PlgZOWwgh12dHDbYcDlsPArw8gyniz1lJa/dB8lhAzns6uSwwYbLYTPbeUeGI7w/v9LXL3oBmSTp7rdn2BPwvHHWFzJ8yMd29yHj48Y9HJ5Pd3+tu3+1u2+d5MQMh4ofkOGL8OmZ1xzS3Tfs7ofMbm7Jti/OcJj9URlOm3hVj6V9hr0GR800PzLDqQezX8arrG433Xxnhj0R/312ZlXdIMNpI2/ZzWtXpAevHbe58xD758btHzbz3tyou48dl58yxv+9PVyY+9gMP/4p/jTDHsxn7ZxRVUcleUmG/5xv2sPpER+eWefu3qtkyXs/vk83zXDNwG519wXd/aTuvkWSJyf5s1rhrZi7+9LuflGGPY53nNCfa+KWOyfGUyoOzdDvWbv9LaykD9nz72NPnwnLkMPksF2QwwYbKoeNp5O9NMMfbD/Ww2mS+z05TA7bBTlssKFy2C4cmOH62RXZJwrI0R8leVBVHdfDnp2XJHlBVe3cM7O5xnPoq+pHquq2438EX81wyPnyJO9J8tWqenpVXa+qDqiqO1XV3few7Vcm+ekM5+C/cmb+q5L8f1V1q/GLtfPc5MtW2rnu/kqG87X/tKpOqKqDarjz2+syXEfw8pWuc4LfS3JSVX1Xd2/PkKCfX1U3qqprVdVtquo+Y9sbZryQuqo2Z7hD2R5V1ZMz3GnsJ/qqe+RukOGHcNHY7mcy7CnaaUeSI2r5u569MsnPVNVxVXWdDO/9u7v7MxNiemQN13Mkw4+2M3w/UlVnV9Wzl3nd06rqvuN358CqelyG9+X9E/pzTTykqu41vge/laF/V9l7uqffwgr7sKffx44kR9cu7lw3s/4Dq+q6GU6dOKCqrrvM3tz9kRy2+uQwOWxVc1iSP89wfdmJ3X3pXvR9XySHrT45TA5btRxWw1HqR1fVweNrfyjDadBvXWnH95kCsrsvynDO9M4BTJ+e4fSId9VwKP/NGc59T5JjxueXZNi782fdfXYPYwGdmOEi109n2Gvwlxkutt2dN47r3NHdH5iZ/7IMCeWccX3/leSpe9HH389wqPp5GRLuuzPsjXhAd3/zmq53N9v7UIa7Su1MQj+d5NoZ7lz15QynO+y8+Pg5GU67+EqSv8+Si8x34zEZLmz/fA3j3VxSVc/s7o8meX6Gz2dHhnPR/3XmdW/NcOOCC6rqaqeNdPdbMnwX/ibDBcK3SfLoiTHdPcm7q+qSDJ/tL3f3p8dlt1wSx6xLx5gvyPDd+cUMe6g/NaE/18Qrk/xGhlMm7pbhYu5d2d1vYSV92NPv43Xjv1+sqvcts/5njdt4Roa9o5dmZo/n/kwOk8OWxC6HXWlD5LAajmA8eXz9BTOf93Jx71fkMDlsSexy2JU2RA7LUED/fIYdHl/O8D1+Wne/YRdtd2vn7XKBPRj3hr2uu++xAWI5Lcn53a34AiaRw4BFJodtHE4dg4m6+/wMd84CWDhyGLDI5LCNY585hRUAAID5cgorAAAAkzgCCQAAwCQLcQ3kYYcd1kcfffR6hwGsofe+971f6O5dDrS7aOQw2P/IYcAi210OW4gC8uijj87WrVvXOwxgDVXVeesdw2qRw2D/I4cBi2x3OcwprAAAAEyigAQAAGASBSQAAACTKCABAACYRAEJAADAJApIAAAAJlFAAgAAMIkCEgAAgEkUkAAAAEyigAQAAGASBSQAAACTKCABAACYRAEJAADAJApIAAAAJlFAAgAAMMmB6x0AAGvrxAffO9u3nTe5/eGbj8oZZ54zx4gAgEWhgATYz2zfdl62nrxpcvstp0wvNgGAfZtTWAEAAJhEAQkAAMAkCkgAAAAmmXsBWVUHVNX7q+pN4/NbVdW7q+oTVfWaqrr2vGMAAABg763FEchfTvKxmefPTfKC7j4myZeTPHENYgAAAGAvzbWArKojkvy3JH85Pq8k90/y+rHJ6UkePs8YAAAAWB3zPgL5R0n+R5Irxuc3TXJxd182Pj8/yeZdvbCqTqqqrVW19aKLLppzmACrSw4DFpkcBixnbgVkVf1Ikgu7+72zs3fRtHf1+u4+tbu3dPeWTZumj1cGsBHIYcAik8OA5Rw4x3XfM8lDq+ohSa6b5EYZjkgeUlUHjkchj0jy+TnGAAAAwCqZ2xHI7j65u4/o7qOTPDrJW7v7J5O8LckjxmaPS/KGecUAAADA6lmPcSCfnuRXquqTGa6JfOk6xAAAAMAKzfMU1u/o7rOTnD1OfyrJ8WuxXQAAAFbPehyBBAAAYAEpIAEAAJhEAQkAAMAkCkgAAAAmUUACAAAwiQISAACASRSQAAAATKKABAAAYBIFJAAAAJMoIAEAAJhEAQkAAMAkCkgAAAAmUUACAAAwiQISAACASRSQAAAATKKABAAAYBIFJAAAAJMoIAEAAJhEAQkAAMAkCkgAAAAmUUACAAAwiQISAACASRSQAAAATKKABAAAYBIFJAAAAJMoIAEAAJhEAQkAAMAkCkgAAAAmUUACAAAwiQISAACASRSQAAAATDK3ArKqrltV76mqD1TVR6rqOeP806rq01V17vg4bl4xAAAAsHoOnOO6v5nk/t19SVUdlOQdVfWP47Jf6+7Xz3HbAAAArLK5FZDd3UkuGZ8eND56XtsDAABgvuZ6DWRVHVBV5ya5MMlZ3f3ucdHvVNUHq+oFVXWdZV57UlVtraqtF1100TzDBFh1chiwyOQwYDlzLSC7+/LuPi7JEUmOr6o7JTk5ye2T3D3JoUmevsxrT+3uLd29ZdOmTfMME2DVyWHAIpPDgOWsyV1Yu/viJGcnOaG7t/fgm0n+d5Lj1yIGAAAA9s4878K6qaoOGaevl+SBST5eVYeP8yrJw5N8eF4xAAAAsHrmeRfWw5OcXlUHZChUX9vdb6qqt1bVpiSV5NwkPzfHGAAAAFgl87wL6weT3GUX8+8/r20CAAAwP2tyDSQAAACLTwEJAADAJApIAAAAJlFAAgAAMIkCEgAAgEkUkAAAAEyigAQAAGASBSQAAACTKCABAACYRAEJAADAJApIAAAAJlFAAgAAMIkCEgAAgEkUkAAAAEyigAQAAGASBSQAAACTKCABAACYRAEJAADAJApIAAAAJlFAAgAAMIkCEgAAgEkUkAAAAEyigAQAAGASBSQAAACTKCABAACYRAEJAADAJApIAAAAJlFAAgAAMIkCEgAAgEkUkAAAAEyigAQAAGCSuRWQVXXdqnpPVX2gqj5SVc8Z59+qqt5dVZ+oqtdU1bXnFQMAAACrZ55HIL+Z5P7dfeckxyU5oaq+P8lzk7ygu49J8uUkT5xjDAAAAKySuRWQPbhkfHrQ+Ogk90/y+nH+6UkePq8YAAAAWD1zvQayqg6oqnOTXJjkrCT/meTi7r5sbHJ+ks3LvPakqtpaVVsvuuiieYYJsOrkMGCRyWHAcuZaQHb35d19XJIjkhyf5A67arbMa0/t7i3dvWXTpk3zDBNg1clhwCKTw4DlrMldWLv74iRnJ/n+JIdU1YHjoiOSfH4tYgAAAGDvzPMurJuq6pBx+npJHpjkY0neluQRY7PHJXnDvGIAAABg9Ry45ybX2OFJTq+qAzIUqq/t7jdV1UeTvLqqfjvJ+5O8dI4xAAAAsErmVkB29weT3GUX8z+V4XpIAAAAFsiaXAMJAADA4lNAAgAAMMk8r4EEAIDdOvHB9872bedNbn/45qNyxpnnzDEiYHcUkAAArJvt287L1pOnjzW55ZTpxSaw+pzCCgAAwCQKSAAAACZRQAIAADCJAhIAAIBJFJAAAABMooAEAABgEgUkAAAAkyggAQAAmEQBCQAAwCQKSAAAACZRQAIAADCJAhIAAIBJFJAAAABMooAEAABgEgUkAAAAkyggAQAAmEQBCQAAwCQKSAAAACZRQAIAADCJAhIAAIBJFJAAAABMooAEAABgkgPXO4BFdOKD753t286b3P7wzUfljDPPmWNEAAAA86eAvAa2bzsvW0/eNLn9llOmF5sAAAAblVNYAQAAmEQBCQAAwCQKSAAAACaZWwFZVbesqrdV1ceq6iNV9cvj/GdX1baqOnd8PGReMQAAALB65nkTncuS/Gp3v6+qbpjkvVV11rjsBd39vDluGwAAgFU2twKyu7cn2T5Of62qPpZk87y2BwAAwHytyTAeVXV0krskeXeSeyZ5SlX9dJKtGY5SfnkXrzkpyUlJcuSRR65FmACrZi1z2ErHpr1wx44k04ciAvY//g4DljP3ArKqDk7yN0me1t1frao/T/JbSXr89/lJnrD0dd19apJTk2TLli097zgBVtNa5rCVjk27+anb5hgNsC/wdxiwnLnehbWqDspQPP51d/9tknT3ju6+vLuvSPKSJMfPMwYAAABWxzzvwlpJXprkY939hzPzD59p9qNJPjyvGAAAAFg9k05hrao7dfdKC717JvmpJB+qqnPHec9M8piqOi7DKayfSfLkFa4XAACAdTD1Gsi/qKprJzktySu7++I9vaC735GkdrHoH6aHBwAAwEYx6RTW7r5Xkp9McsskW6vqlVX1oLlGBgAAwIYy+RrI7v5EkmcleXqS+yT5k6r6eFX993kFBwAAwMYxqYCsqu+tqhck+ViS+yc5sbvvME6/YI7xAQAAsEFMvQbyhRmG3Hhmd1+6c2Z3f76qnjWXyAAAANhQphaQD0lyaXdfniRVda0k1+3ub3T3y+cWHQAAABvG1Gsg35zkejPPrz/OAwAAYD8xtYC8bndfsvPJOH39+YQEAADARjS1gPx6Vd1155OquluSS3fTHgAAgH3M1Gsgn5bkdVX1+fH54Ul+fD4hAQAAsBFNKiC7+9+r6vZJbpekkny8u78918gAAADYUKYegUySuyc5enzNXaoq3f1Xc4kKAACADWdSAVlVL09ymyTnJrl8nN1JFJAAAAD7ialHILckuWN39zyDAQAAYOOaehfWDyf5rnkGAgAAwMY29QjkYUk+WlXvSfLNnTO7+6FziQoAAIANZ2oB+ex5BgEAAMDGN3UYj7dX1VFJjunuN1fV9ZMcMN/QAAAA2EgmXQNZVU9K8vokLx5nbU7yf+YVFAAAABvP1FNYfzHJ8UnenSTd/YmqutncolpjJz743tm+7bzJ7S/csSPJpvkFBAAAsAFNLSC/2d3fqqokSVUdmGEcyH3C9m3nZevJ0wvCzU/dNsdoAAAANqapw3i8vaqemeR6VfWgJK9Lcsb8wgIAAGCjmVpAPiPJRUk+lOTJSf4hybPmFRQAAAAbz9S7sF6R5CXjAwAAgP3QpAKyqj6dXVzz2N23XvWIAAAA2JCm3kRny8z0dZM8Msmhqx8OAAAAG9WkayC7+4szj23d/UdJ7j/n2AAAANhApp7CeteZp9fKcETyhnOJCAAAgA1p6imsz5+ZvizJZ5I8atWjAQAAYMOaehfW+807EAAAADa2qaew/srulnf3H65OOAAAAGxUK7kL692TvHF8fmKSc5J8bh5BAQAAsPFMLSAPS3LX7v5aklTVs5O8rrt/drkXVNUtk/xVku9KckWSU7v7j6vq0CSvSXJ0xmspu/vL17QDAAAArI1Jw3gkOTLJt2aefytDAbg7lyX51e6+Q5LvT/KLVXXHJM9I8pbuPibJW8bnAAAAbHBTj0C+PMl7qurvknSSH81wdHFZ3b09yfZx+mtV9bEkm5M8LMl9x2anJzk7ydNXGjgAAABra+pdWH+nqv4xyQ+Os36mu98/dSNVdXSSuyR5d5Kbj8Vlunt7Vd1smdeclOSkJDnyyCOnbgpgQ5DDgEUmhwHLmXoKa5JcP8lXu/uPk5xfVbea8qKqOjjJ3yR5Wnd/derGuvvU7t7S3Vs2bdq0gjAB1p8cBiwyOQxYzqQCsqp+I8NppiePsw5K8ooJrzsoQ/H41939t+PsHVV1+Lj88CQXrjRoAAAA1t7UI5A/muShSb6eJN39+SQ33N0LqqqSvDTJx5aME/nGJI8bpx+X5A0rCRgAAID1MfUmOt/q7q6qTpKqusGE19wzyU8l+VBVnTvOe2aS30vy2qp6YpLPJnnkCmMGAABgHUwtIF9bVS9OckhVPSnJE5K8ZHcv6O53JKllFj9geogAAABsBFPvwvq8qnpQkq8muV2SX+/us+YaGQAAABvKHgvIqjogyT939wOTKBoBAAD2U3ssILv78qr6RlXduLu/shZB7WsuuGBHthx71OT2h28+Kmecec4cIwIAAFi5qddA/leGm+GclfFOrEnS3b80l6j2MX3FZdl68vQxlLacct4cowEAALhmphaQfz8+AAAA2E/ttoCsqiO7+7PdffpaBQQAAMDGdK09LP8/Oyeq6m/mHAsAAAAb2J4KyNlxHG89z0AAAADY2PZUQPYy0wAAAOxn9nQTnTtX1VczHIm83jid8Xl3943mGh0AAAAbxm4LyO4+YK0CAQCA9Xbig++d7dumD6lm/G72N1OH8QAAgH3e9m3nGb8bdmNP10ACAABAEgUkAAAAEykgAQAAmEQBCQAAwCQKSAAAACZRQAIAADCJYTwAAFgYF1ywI1uOPWpy+y996eIceughk9tfuGNHkunDeMD+RgEJAMDC6CsuW9E4jZufui1bTz5mRe2B5TmFFQAAgEkUkAAAAEyigAQAAGASBSQAAACTKCABAACYRAEJAADAJApIAAAAJlFAAgAAMIkCEgAAgEkUkAAAAEyigAQAAGCSuRWQVfWyqrqwqj48M+/ZVbWtqs4dHw+Z1/YBAABYXfM8AnlakhN2Mf8F3X3c+PiHOW4fAACAVTS3ArK7z0nypXmtHwAAgLW1HtdAPqWqPjie4nqT5RpV1UlVtbWqtl500UVrGR/AXpPDgEUmhwHLWesC8s+T3CbJcUm2J3n+cg27+9Tu3tLdWzZt2rRW8QGsCjkMWGRyGLCcNS0gu3tHd1/e3VckeUmS49dy+wAAAFxza1pAVtXhM09/NMmHl2sLAADAxnLgvFZcVa9Kct8kh1XV+Ul+I8l9q+q4JJ3kM0mePK/tAwAAsLrmVkB292N2Mful89oeAAAA87Ued2EFAABgASkgAQAAmGRup7ByzV1wwY5sOfaoye0P33xUzjjznDlGBAAAoIDckPqKy7L15OljLm055bw5RgMAADBwCisAAACTKCABAACYRAEJAADAJApIAAAAJlFAAgAAMIkCEgAAgEkM4wHAbhmbFgDYSQEJwG4ZmxYA2MkprAAAAEyigAQAAGASBSQAAACTKCABAACYRAEJAADAJApIAAAAJlFAAgAAMIkCEgAAgEkUkAAAAEyigAQAAGASBSQAAACTKCABAACYRAEJAADAJApIAAAAJlFAAgAAMIkCEgAAgEkUkAAAAEyigAQAAGASBSQAAACTKCABAACYZG4FZFW9rKourKoPz8w7tKrOqqpPjP/eZF7bBwAAYHXN8wjkaUlOWDLvGUne0t3HJHnL+BwAAIAFMLcCsrvPSfKlJbMfluT0cfr0JA+f1/YBAABYXWt9DeTNu3t7koz/3my5hlV1UlVtraqtF1100ZoFCLAa5DBgkclhwHI27E10uvvU7t7S3Vs2bdq03uEArIgcBiwyOQxYzloXkDuq6vAkGf+9cI23DwAAwDW01gXkG5M8bpx+XJI3rPH2AQAAuIbmOYzHq5K8M8ntqur8qnpikt9L8qCq+kSSB43PAQAAWAAHzmvF3f2YZRY9YF7bBAAAYH427E10AAAA2FgUkAAAAEwyt1NYAdg/XXDBjmw59qjJ7Q/ffFTOOPOcOUYEAKwWBSQAq6qvuCxbT54+btyWU86bYzQAwGpyCisAAACTKCABAACYRAEJAADAJApIAAAAJlFAAgAAMIkCEgAAgEkUkAAAAEyigAQAAGASBSQAAACTKCABAACYRAEJAADAJApIAAAAJlFAAgAAMIkCEgAAgEkOXO8AAABgUV1wwY5sOfaoye0P33xUzjjznDlGBPOlgAQAgGuor7gsW0/eNLn9llPOm2M0MH9OYQUAAGASBSQAAACTKCABAACYRAEJAADAJApIAAAAJlFAAgAAMIlhPPYBxh8CAADWggJyH2D8IQAAYC04hRUAAIBJFJAAAABMooAEAABgknW5BrKqPpPka0kuT3JZd29ZjzgAAACYbj1vonO/7v7COm4fAACAFXAKKwAAAJOs1xHITnJmVXWSF3f3qUsbVNVJSU5KkiOPPHJFKz/xwffO9m3Th6q4cMeOJNOHwQDYk73JYQDrTQ4DlrNeBeQ9u/vzVXWzJGdV1ce7+yoj249F5alJsmXLll7JyrdvO29F4yJufuq2laweYI/2JocBrDc5DFjOukwkjKoAAAxzSURBVJzC2t2fH/+9MMnfJTl+PeIAAABgujUvIKvqBlV1w53TSR6c5MNrHQcAAAArsx6nsN48yd9V1c7tv7K7/2kd4gAAAGAF1ryA7O5PJbnzWm8XAACAvWMYDwAAACZRQAIAADDJeg3jwTq64IId2XLsUSt6zeGbj8oZZ56z54bAXlnpOLaJsWwBgLWjgNwP9RWXrWiczCTZcsrK/qAFrpmVjmObGMsWAFg7TmEFAABgEgUkAAAAkyggAQAAmEQBCQAAwCQKSAAAACZRQAIAADCJYTwAAFgVxrLds5WOx/2lL12cQw89ZEXbMH4386SABABgVRjLds9WOh735qduy9aTj1nRNozfzTw5hRUAAIBJFJAAAABMooAEAABgEgUkAAAAkyggAQAAmEQBCQAAwCSG8WCSlY5ZtNLxh1Y6bpTxjWDfYUw0gH2bv/P2LQpIJlnpmEUrHX9opeNGGd8I9h3GRAPYt/k7b9/iFFYAAAAmUUACAAAwiQISAACASRSQAAAATKKABAAAYBIFJAAAAJMYxgOAfd68x5qcd3tjogEbyUrHdbxwx44k04fxYGNTQAKwz5v3WJPzbm9MNGAjWem4jpufum2O0bDWnMIKAADAJApIAAAAJlFAAgAAMMm6FJBVdUJV/UdVfbKqnrEeMQAAALAya15AVtUBSV6U5IeT3DHJY6rqjmsdBwAAACuzHkcgj0/yye7+VHd/K8mrkzxsHeIAAABgBaq713aDVY9IckJ3/+z4/KeSfF93P2VJu5OSnDQ+vV2S/1jBZg5L8oVVCHcj2x/6mOjnvmYl/Tyquxd20Ki9yGG+C/sW/dx3rLSPcti+TT/3LftDP1cth61HAfnIJD+0pIA8vrufuorb2NrdW1ZrfRvR/tDHRD/3NftLP/fG/vIe6ee+ZX/o5/7Qx9Wwv7xP+rlv2R/6uZp9XI9TWM9PcsuZ50ck+fw6xAEAAMAKrEcB+e9JjqmqW1XVtZM8Oskb1yEOAAAAVuDAtd5gd19WVU9J8s9JDkjysu7+yCpv5tRVXt9GtD/0MdHPfc3+0s+9sb+8R/q5b9kf+rk/9HE17C/vk37uW/aHfq5aH9f8GkgAAAAW03qcwgoAAMACUkACAAAwycIWkFV1QlX9R1V9sqqesYvl16mq14zL311VR699lHtvQj9/pao+WlUfrKq3VNVR6xHn3tpTP2faPaKquqoW7lbLU/pYVY8aP8+PVNUr1zrG1TDhO3tkVb2tqt4/fm8fsh5xrjc57DvL5bAFIYd9Z7kcFjlsZrkctiDksO8s3/sc1t0L98hw853/THLrJNdO8oEkd1zS5heS/MU4/egkr1nvuOfUz/sluf44/fP7aj/HdjdMck6SdyXZst5xz+GzPCbJ+5PcZHx+s/WOe079PDXJz4/Td0zymfWOe4O+T3LYgjzksO+0kcP2k4ccdpU2ctgCPOSwq7TZ6xy2qEcgj0/yye7+VHd/K8mrkzxsSZuHJTl9nH59kgdUVa1hjKthj/3s7rd19zfGp+/KMK7mopnyeSbJbyX5/ST/tZbBrZIpfXxSkhd195eTpLsvXOMYV8OUfnaSG43TN87+OQ6sHDaSwxaGHHYlOUwO+w45bGHIYVfa6xy2qAXk5iSfm3l+/jhvl226+7IkX0ly0zWJbvVM6eesJyb5x7lGNB977GdV3SXJLbv7TWsZ2Cqa8ll+d5Lvrqp/rap3VdUJaxbd6pnSz2cneWxVnZ/kH5I8dW1C21DksF2TwzYuOexKz44cJoftmhy2cclhV3p29jKHrfk4kKtkV3uwlo5HMqXNRje5D1X12CRbktxnrhHNx277WVXXSvKCJI9fq4DmYMpneWCG0yfum2EP5r9U1Z26++I5x7aapvTzMUlO6+7nV9U9krx87OcV8w9vw5DDljaUwzY6OexKcpgcdvWGcthGJ4ddaa9z2KIegTw/yS1nnh+Rqx9+/U6bqjowwyHaL61JdKtnSj9TVQ9M8j+TPLS7v7lGsa2mPfXzhknulOTsqvpMku9P8sYFu4B76nf2Dd397e7+dJL/yJDIFsmUfj4xyWuTpLvfmeS6SQ5bk+g2Djlshhy2EOSwK8lhcthVyGELQQ670t7nsPW+2POaPDLsIfhUklvlygtEj13S5hdz1Yu3X7vecc+pn3fJcLHsMesd7zz7uaT92Vm8i7enfJYnJDl9nD4swykIN13v2OfQz39M8vhx+g5jYqv1jn0Dvk9y2II85LDvtJHD9pOHHHaVNnLYAjzksKu02esctu4d3Ys36CFJ/u/4o/2f47zfzLD3Jxmq6dcl+WSS9yS59XrHPKd+vjnJjiTnjo83rnfM8+jnkrYLl7gmfpaV5A+TfDTJh5I8er1jnlM/75jkX8ekdm6SB693zBv0fZLDFughh8lh+9tDDpPD1jvmOXyWctjER40rAgAAgN1a1GsgAQAAWGMKSAAAACZRQAIAADCJAhIAAIBJFJAAAABMooDcx1XV5VV1blV9uKrOqKpD1jumWVX1m+MAvFPanlRVHx8f76mqe80sO3vngLZV9Zmq+lBVfaCqzqyq71pmfa+vqluP05cs0+bsXQ2UW1WPr6oX7ibW76mq06b0C1ieHCaHwSKTw+SwfZECct93aXcf1913SvKlDAP7zkVVHbjS13T3r3f3m3exrgOWPP+RJE9Ocq/uvn2Sn0vyyuWSUpL7dfedk2xN8sxdrP/YJAd096dWGvMU3f2hJEdU1ZHzWD/sR+QwOQwWmRwmh+1zFJD7l3cm2bzzSVX9WlX9e1V9sKqeM867QVX9/bjX6MNV9ePj/LtV1dur6r1V9c9Vdfg4/+yq+t2qenuS/znudbrWuOz6VfW5qjqoqo6rqneN2/q7qrrJ2Oa0qnrEOP2Zqvr1qnpHkkcuif3pSX6tu7+QJN39viSnZ8+J+Jwkt93F/J9M8obZGVX1/Kp6X1W9pao2zSx6bFX92/h+HL90RVX1yHHZB6rqnJlFZyR59B7iA6aTw64kh8HikcOuJIctMAXkfqKGPUkPSPLG8fmDkxyT5PgkxyW5W1XdO8kJST7f3Xce95b9U1UdlORPkzyiu++W5GVJfmdm9Yd09326+zlJPpDkPuP8E5P8c3d/O8lfJXl6d39vkg8l+Y1lQv2v7r5Xd796yfxjk7x3ybyt4/zd+ZFxe0vdc8n6bpDkfd191yRvXxLfDbr7B5L8Qoa+L/XrSX5o3NP20CXx/eAe4gMmkMOuRg6DBSKHXY0ctsAUkPu+61XVuUm+mOTQJGeN8x88Pt6f5H1Jbp8hkX0oyQOr6rlV9YPd/ZUkt0typyRnjet6VpIjZrbxmiXTPz5OPzrJa6rqxhmS29vH+acnufcy8b5mmfm7Ukl6mWVvG2O9UZJTdrH88CQXzTy/Ymbbr0hyr5llr0qS7j4nyY3q6tcv/GuS06rqSUlmT/m4MMktJvQDWJ4cJofBIpPD5LB9zorPlWbhXNrdx43J400ZTjX4kww/+lO6+8VLX1BVd0vykCSnVNWZSf4uyUe6+x7LbOPrM9NvHF93aJK7JXlrkoNXEO/Xl5n/0Zn17XTXcf6u3G/naRbLuDTJdXezvJeZvtrz7v65qvq+JP8tyblVdVx3f3Fc/6W72QawZ3LYrslhsBjksF2TwxaYI5D7iXEP1i8l+f/HUyH+OckTqurgJKmqzVV1s6q6RZJvdPcrkjwvQ3L4jySbquoeY9uDarj4eVfbuSTJe5L8cZI3dffl47a/XFU7TyP4qQynJ6zE7yd5blXddIzhuCSPT/JnK1zPTh/LVc/Jv1aSR4zTP5HkHTPLdl5/cK8kXxn78x1VdZvufnd3/3qSLyS55bjou5N8+BrGB8yQw65GDoMFIoddjRy2wByB3I909/ur6gNJHt3dL6+qOyR5Z1UlySVJHpvhx/wHVXVFkm8n+fnu/tZ4gfWfjHvQDkzyR0k+ssymXpPkdUnuOzPvcUn+oqqun+RTSX5mhbG/sao2J/m3quokX0vy2O7evpL1zPj7Mb6ddx77epJjq+q9Sb6SK0//SIak+28ZTsN4wi7W9QdVdUyGvYlvyXD9QZLcb9wOsArksKuQw2DByGFXIYctsOpe7tRl2HdV1fWSvC3JPbv78jms/zoZ9u7dq7svW+31A/s3OQxYZHLYYlNAst+qqh9K8rHu/uwc1n1Mks3dffZqrxsgkcOAxSaHLS4FJAAAAJO4iQ4AAACTKCABAACYRAEJAADAJApIAAAAJlFAAgAAMMn/A9SEFDVbDYdtAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The average OIP for large number of realizations = 0.146\n",
      "The average OIP for small number of realizations = 0.135, 0.121,0.142\n"
     ]
    }
   ],
   "source": [
    "f, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True)\n",
    "ax1.hist(OIP1,bins=np.linspace(0.0,0.8,20),color=\"darkorange\",edgecolor='black',alpha=0.8); ax1.set_xlabel(\"Reservoir OIP (bbls)\"); ax1.set_ylabel(\"Frequency\"); ax1.set_title(\"Reservoir OIP Realizations, Sample Set 1\")\n",
    "    \n",
    "ax2.hist(OIP2,bins=np.linspace(0.0,0.8,20),color=\"darkorange\",edgecolor='black',alpha=0.8); ax2.set_xlabel(\"Reservoir OIP (bbls)\"); ax2.set_title(\"Reservoir OIP Realizations, Sample Set 2\")\n",
    "\n",
    "ax3.hist(OIP3,bins=np.linspace(0.0,0.8,20),color=\"darkorange\",edgecolor='black',alpha=0.8); ax3.set_xlabel(\"Reservoir OIP (bbls)\"); ax3.set_title(\"Reservoir OIP Realizations, Sample Set 3\")\n",
    "\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.2, hspace=0.2)\n",
    "plt.show()\n",
    "\n",
    "print('The average OIP for large number of realizations = ' + str(round(np.average(OIP),3)))\n",
    "print('The average OIP for small number of realizations = ' + str(round(np.average(OIP1),3)) + ', ' + str(round(np.average(OIP2),3)) + ',' + str(round(np.average(OIP3),3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try setting the number of realizations, $L$, larger and you will see that the 3 sets of realizations converge, they become stable. This is why we set $L$ large. \n",
    "\n",
    "#### Non-Parametric Monte Carlos Simulation\n",
    "\n",
    "I have demonstrated the case of MCS with parametric distributions. Of course, we may have non-parametric distributions. In this case we are working with a list of values (typically in an 1D ndarray). For this case we simply substitute in the command:\n",
    "\n",
    "```p\n",
    "apor = np.random.choice(array_average_porosity,size=L)\n",
    "```\n",
    "\n",
    "The result is a ndarray with $L$ Monte Carlo simulations from the array 'array_average_porosity'.  This is simply $L$ random draws with replacement (so the same value may be selected again), from the list of values in the ndarray."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comments\n",
    "\n",
    "This was a basic demonstration of Monte Carlo simulation for uncertainty analysis. A lot more could be done, for example, more complicated transfer functions and a combination of non-parametric and parametric distributions. Also, one could integrate relationships between the variables (we assumed independent here).\n",
    "\n",
    "I have other demonstrations on the basics of working with DataFrames, ndarrays, univariate statistics, plotting data, declustering, data transformations, trend modeling, multivariate analysis and many other workflows available at https://github.com/GeostatsGuy/PythonNumericalDemos and https://github.com/GeostatsGuy/GeostatsPy. \n",
    "  \n",
    "I hope this was helpful,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "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"
   ]
  }
 ],
 "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
}