Interactive_Simulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
    "## Interactive Sequential Gaussian Simulation Demonstration\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",
    "\n",
    "\n",
    "### The Interactive Workflow\n",
    "\n",
    "Here's a simple workflow for calculating the sequential Gaussian simulation with sequential use of simple kriging estimate and the estimation variance for a local uncertainty model and Monte Carlo simulation.\n",
    "\n",
    "* we use a 'toy problem' with only 3 data for speed and interpretability of the results\n",
    "\n",
    "#### Spatial Estimation\n",
    "\n",
    "Consider the case of making an estimate at some unsampled location, $𝑧(\\bf{u}_0)$, where $z$ is the property of interest (e.g. porosity etc.) and $𝐮_0$ is a location vector describing the unsampled location.\n",
    "\n",
    "How would you do this given data, $𝑧(\\bf{𝐮}_1)$, $𝑧(\\bf{𝐮}_2)$, and $𝑧(\\bf{𝐮}_3)$?\n",
    "\n",
    "It would be natural to use a set of linear weights to formulate the estimator given the available data.\n",
    "\n",
    "\\begin{equation}\n",
    "z^{*}(\\bf{u}) = \\sum^{n}_{\\alpha = 1} \\lambda_{\\alpha} z(\\bf{u}_{\\alpha})\n",
    "\\end{equation}\n",
    "\n",
    "We could add an unbiasedness constraint to impose the sum of the weights equal to one.  What we will do is assign the remainder of the weight (one minus the sum of weights) to the global average; therefore, if we have no informative data we will estimate with the global average of the property of interest.\n",
    "\n",
    "\\begin{equation}\n",
    "z^{*}(\\bf{u}) = \\sum^{n}_{\\alpha = 1} \\lambda_{\\alpha} z(\\bf{u}_{\\alpha}) + \\left(1-\\sum^{n}_{\\alpha = 1} \\lambda_{\\alpha} \\right) \\overline{z}\n",
    "\\end{equation}\n",
    "\n",
    "We will make a stationarity assumption, so let's assume that we are working with residuals, $y$. \n",
    "\n",
    "\\begin{equation}\n",
    "y^{*}(\\bf{u}) = z^{*}(\\bf{u}) - \\overline{z}(\\bf{u})\n",
    "\\end{equation}\n",
    "\n",
    "If we substitute this form into our estimator the estimator simplifies, since the mean of the residual is zero.\n",
    "\n",
    "\\begin{equation}\n",
    "y^{*}(\\bf{u}) = \\sum^{n}_{\\alpha = 1} \\lambda_{\\alpha} y(\\bf{u}_{\\alpha})\n",
    "\\end{equation}\n",
    "\n",
    "while satisfying the unbaisedness constraint.  \n",
    "\n",
    "#### Kriging\n",
    "\n",
    "Now the next question is what weights should we use?  \n",
    "\n",
    "We could use equal weighting, $\\lambda = \\frac{1}{n}$, and the estimator would be the average of the local data applied for the spatial estimate. This would not be very informative.\n",
    "\n",
    "We could assign weights considering the spatial context of the data and the estimate:\n",
    "\n",
    "* **spatial continuity** as quantified by the variogram (and covariance function)\n",
    "* **redundancy** the degree of spatial continuity between all of the available data with themselves \n",
    "* **closeness** the degree of spatial continuity between the avaiable data and the estimation location\n",
    "\n",
    "The kriging approach accomplishes this, calculating the best linear unbiased weights for the local data to estimate at the unknown location.  The derivation of the kriging system and the resulting linear set of equations is available in the lecture notes.  Furthermore kriging provides a measure of the accuracy of the estimate!  This is the kriging estimation variance (sometimes just called the kriging variance).\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma^{2}_{E}(\\bf{u}) = C(0) - \\sum^{n}_{\\alpha = 1} \\lambda_{\\alpha} C(\\bf{u}_0 - \\bf{u}_{\\alpha})\n",
    "\\end{equation}\n",
    "\n",
    "What is 'best' about this estimate? Kriging estimates are best in that they minimize the above estimation variance. \n",
    "\n",
    "#### Properties of Kriging\n",
    "\n",
    "Here are some important properties of kriging:\n",
    "\n",
    "* **Exact interpolator** - kriging estimates with the data values at the data locations\n",
    "* **Kriging variance** can be calculated before getting the sample information, as the kriging estimation variance is not dependent on the values of the data nor the kriging estimate, i.e. the kriging estimator is homoscedastic. \n",
    "* **Spatial context** - kriging takes into account, furthermore to the statements on spatial continuity, closeness and redundancy we can state that kriging accounts for the configuration of the data and structural continuity of the variable being estimated.\n",
    "* **Scale** - kriging may be generalized to account for the support volume of the data and estimate. We will cover this later.\n",
    "* **Multivariate** - kriging may be generalized to account for multiple secondary data in the spatial estimate with the cokriging system. We will cover this later.\n",
    "* **Smoothing effect** of kriging can be forecast. We will use this to build stochastic simulations later.\n",
    "\n",
    "#### Spatial Continuity \n",
    "\n",
    "**Spatial Continuity** is the correlation between values over distance.\n",
    "\n",
    "* No spatial continuity – no correlation between values over distance, random values at each location in space regardless of separation distance.\n",
    "\n",
    "* Homogenous phenomenon have perfect spatial continuity, since all values as the same (or very similar) they are correlated. \n",
    "\n",
    "We need a statistic to quantify spatial continuity! A convenient method is the Semivariogram.\n",
    "\n",
    "#### The Semivariogram\n",
    "\n",
    "Function of difference over distance.\n",
    "\n",
    "* The expected (average) squared difference between values separated by a lag distance vector (distance and direction), $h$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\gamma(\\bf{h}) = \\frac{1}{2 N(\\bf{h})} \\sum^{N(\\bf{h})}_{\\alpha=1} (z(\\bf{u}_\\alpha) - z(\\bf{u}_\\alpha + \\bf{h}))^2  \n",
    "\\end{equation}\n",
    "\n",
    "where $z(\\bf{u}_\\alpha)$ and $z(\\bf{u}_\\alpha + \\bf{h})$ are the spatial sample values at tail and head locations of the lag vector respectively.\n",
    "\n",
    "* Calculated over a suite of lag distances to obtain a continuous function.\n",
    "\n",
    "* the $\\frac{1}{2}$ term converts a variogram into a semivariogram, but in practice the term variogram is used instead of semivariogram.\n",
    "* We prefer the semivariogram because it relates directly to the covariance function, $C_x(\\bf{h})$ and univariate variance, $\\sigma^2_x$:\n",
    "\n",
    "\\begin{equation}\n",
    "C_x(\\bf{h}) = \\sigma^2_x - \\gamma(\\bf{h})\n",
    "\\end{equation}\n",
    "\n",
    "Note the correlogram is related to the covariance function as:\n",
    "\n",
    "\\begin{equation}\n",
    "\\rho_x(\\bf{h}) = \\frac{C_x(\\bf{h})}{\\sigma^2_x}\n",
    "\\end{equation}\n",
    "\n",
    "The correlogram provides of function of the $\\bf{h}-\\bf{h}$ scatter plot correlation vs. lag offset $\\bf{h}$.  \n",
    "\n",
    "\\begin{equation}\n",
    "-1.0 \\le \\rho_x(\\bf{h}) \\le 1.0\n",
    "\\end{equation}\n",
    "\n",
    "#### Sequential Gaussian Simulation\n",
    "\n",
    "With sequential Gaussian simulation we build on kriging by:\n",
    "\n",
    "* adding a random residual with the missing variance\n",
    "\n",
    "* sequentially adding the simulated values as data to correct the covariance between the simulated values\n",
    "\n",
    "I have more on this topic at [Simulation YouTube Lecture](https://www.youtube.com/watch?v=3cLqK3lR56Y&list=PLG19vXLQHvSB-D4XKYieEku9GQMQyAzjJ&index=45&t=813s).\n",
    "\n",
    "#### Objective \n",
    "\n",
    "In the PGE 383: Stochastic Subsurface Modeling class I want to provide hands-on experience with building subsurface modeling workflows. Python provides an excellent vehicle to accomplish this. I have coded a package called GeostatsPy with GSLIB: Geostatistical Library (Deutsch and Journel, 1998) functionality that provides basic building blocks for building subsurface modeling workflows. \n",
    "\n",
    "The objective is to remove the hurdles of subsurface modeling workflow construction by providing building blocks and sufficient examples. This is not a coding class per se, but we need the ability to 'script' workflows working with numerical methods.    \n",
    "\n",
    "#### Getting Started\n",
    "\n",
    "Here's the steps to get setup in Python with the GeostatsPy package:\n",
    "\n",
    "1. Install Anaconda 3 on your machine (https://www.anaconda.com/download/). \n",
    "2. From Anaconda Navigator (within Anaconda3 group), go to the environment tab, click on base (root) green arrow and open a terminal. \n",
    "3. In the terminal type: pip install geostatspy. \n",
    "4. Open Jupyter and in the top block get started by copy and pasting the code block below from this Jupyter Notebook to start using the geostatspy functionality. \n",
    "\n",
    "You will need to copy the data file to your working directory.  They are available here:\n",
    "\n",
    "* Tabular data - sample_data.csv at https://git.io/fh4gm.\n",
    "\n",
    "There are exampled 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. \n",
    "\n",
    "#### Load the required libraries\n",
    "\n",
    "The following code loads the required libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "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": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import os                                               # to set current working directory \n",
    "import sys                                              # supress output to screen for interactive variogram modeling\n",
    "import io\n",
    "import numpy as np                                      # arrays and matrix math\n",
    "import pandas as pd                                     # DataFrames\n",
    "import matplotlib.pyplot as plt                         # plotting\n",
    "from matplotlib.pyplot import cm                        # color maps\n",
    "from matplotlib.patches import Ellipse                  # plot an ellipse\n",
    "import math                                             # sqrt operator\n",
    "import random                                           # random simulation locations\n",
    "from copy import copy                                   # copy a colormap\n",
    "from scipy.stats import norm\n",
    "from ipywidgets import interactive                      # widgets and interactivity\n",
    "from ipywidgets import widgets                            \n",
    "from ipywidgets import Layout\n",
    "from ipywidgets import Label\n",
    "from ipywidgets import VBox, HBox\n",
    "from scipy.stats import norm                            # Gaussian distribution\n",
    "import scipy.stats as stats                             # trimmed statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you get a package import error, you may have to first install some of these packages. This can usually be accomplished by opening up a command window on Windows and then typing 'python -m pip install [package-name]'. More assistance is available with the respective package docs.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simple, Simple Kriging Function\n",
    "\n",
    "Let's write a fast Python function to take data points and unknown location and provide the:\n",
    "\n",
    "* **simple kriging estimate**\n",
    "\n",
    "* **simple kriging variance / estimation variance**\n",
    "\n",
    "* **simple kriging weights**\n",
    "\n",
    "This provides a fast method for small datasets, with less parameters (no search parameters) and the ability to see the simple kriging weights.\n",
    "\n",
    "* we use it here for fast, flexible application of sequential simulation\n",
    "\n",
    "* the method will not work with only ones simulation location so we send 2 and only use the first result (the 2nd is always a dummy location in the workflow below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simple_simple_krige(df,xcol,ycol,vcol,dfl,xlcol,ylcol,vario,skmean):\n",
    "# load the variogram\n",
    "    nst = vario['nst']; pmx = 9999.9\n",
    "    cc = np.zeros(nst); aa = np.zeros(nst); it = np.zeros(nst)\n",
    "    ang = np.zeros(nst); anis = np.zeros(nst)\n",
    "    nug = vario['nug']; sill = nug \n",
    "    cc[0] = vario['cc1']; sill = sill + cc[0]\n",
    "    it[0] = vario['it1']; ang[0] = vario['azi1']; \n",
    "    aa[0] = vario['hmaj1']; anis[0] = vario['hmin1']/vario['hmaj1'];\n",
    "    if nst == 2:\n",
    "        cc[1] = vario['cc2']; sill = sill + cc[1]\n",
    "        it[1] = vario['it2']; ang[1] = vario['azi2']; \n",
    "        aa[1] = vario['hmaj2']; anis[1] = vario['hmin2']/vario['hmaj2'];    \n",
    "\n",
    "# set up the required matrices\n",
    "    rotmat, maxcov = geostats.setup_rotmat(nug,nst,it,cc,ang,pmx)    \n",
    "    ndata = len(df); a = np.zeros([ndata,ndata]); r = np.zeros(ndata); s = np.zeros(ndata); rr = np.zeros(ndata)\n",
    "    nest = len(dfl)\n",
    "\n",
    "    est = np.zeros(nest); var = np.full(nest,sill); weights = np.zeros([nest,ndata])\n",
    "\n",
    "# Make and solve the kriging matrix, calculate the kriging estimate and variance \n",
    "    for iest in range(0,nest):\n",
    "        for idata in range(0,ndata):\n",
    "            for jdata in range(0,ndata):\n",
    "                a[idata,jdata] = geostats.cova2(df[xcol].values[idata],df[ycol].values[idata],df[xcol].values[jdata],df[ycol].values[jdata],\n",
    "                                        nst,nug,pmx,cc,aa,it,ang,anis,rotmat,maxcov)\n",
    "            r[idata] = geostats.cova2(df[xcol].values[idata],df[ycol].values[idata],dfl[xlcol].values[iest],dfl[ylcol].values[iest],\n",
    "                                        nst,nug,pmx,cc,aa,it,ang,anis,rotmat,maxcov)\n",
    "            rr[idata] = r[idata]\n",
    "        \n",
    "        s = geostats.ksol_numpy(ndata,a,r)    \n",
    "        sumw = 0.0\n",
    "        for idata in range(0,ndata):                          \n",
    "            sumw = sumw + s[idata]\n",
    "            weights[iest,idata] = s[idata]\n",
    "            est[iest] = est[iest] + s[idata]*df[vcol].values[idata]\n",
    "            var[iest] = var[iest] - s[idata]*rr[idata]\n",
    "        est[iest] = est[iest] + (1.0-sumw)*skmean\n",
    "    return est,var,weights "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Interactive Sequential Simulation to Random Points Method\n",
    "\n",
    "For this first interactive method we will perform sequential simulation:\n",
    "\n",
    "* at **nsim** random point locations in the area of interest.\n",
    "\n",
    "The following code includes:\n",
    "\n",
    "* dashboard with number of simulation locations, variogram model and data locations \n",
    "\n",
    "* plots of variogram model, data locations with point scaled by weights and uncertainty distribution at the unknown location\n",
    "\n",
    "Let's first set up the model area of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X extents [0.0,1000.0] and Y entents [0.0,1000.0]\n"
     ]
    }
   ],
   "source": [
    "csiz = 100; xmn = csiz * 0.5; nx = 10; ymn = csiz * 0.5; ny = 10 \n",
    "xmin = xmn - csiz * 0.5; xmax = xmin + nx * csiz\n",
    "ymin = ymn - csiz * 0.5; ymax = ymin + ny * csiz\n",
    "print('X extents [' + str(xmin) + ',' + str(xmax) + '] and Y entents [' + str(ymin) + ',' + str(ymax) + ']')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's set up our dash board."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings; warnings.simplefilter('ignore')\n",
    "\n",
    "# dashboard: number of simulation locations and variogram parameters\n",
    "style = {'description_width': 'initial'}\n",
    "l = widgets.Text(value='                                              Sequential Simulation, Michael Pyrcz, Associate Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n",
    "nsim = widgets.IntSlider(min = 0, max = 99, value = 5, step = 1, description = 'nsim',orientation='vertical',\n",
    "                          layout=Layout(width='25px', height='200px'),continuous_update = False)\n",
    "nsim.style.handle_color = 'gray'\n",
    "nug = widgets.FloatSlider(min = 0, max = 1.0, value = 0.0, step = 0.1, description = 'nug',orientation='vertical',\n",
    "                          layout=Layout(width='25px', height='200px'),continuous_update = False)\n",
    "nug.style.handle_color = 'gray'\n",
    "it1 = widgets.Dropdown(options=['Spherical', 'Exponential', 'Gaussian'],value='Spherical',\n",
    "    description='Type1:',disabled=False,layout=Layout(width='180px', height='30px'), style=style,continuous_update = False)\n",
    "\n",
    "azi = widgets.FloatSlider(min=0, max = 360, value = 0, step = 22.5, description = 'azi',\n",
    "                        orientation='vertical',layout=Layout(width='40px', height='200px'),continuous_update = False)\n",
    "azi.style.handle_color = 'gray'\n",
    "hmaj1 = widgets.FloatSlider(min=0.01, max = 10000.0, value = 100.0, step = 25.0, description = 'hmaj1',\n",
    "                        orientation='vertical',layout=Layout(width='40px', height='200px'),continuous_update = False)\n",
    "hmaj1.style.handle_color = 'gray'\n",
    "hmin1 = widgets.FloatSlider(min = 0.01, max = 10000.0, value = 100.0, step = 25.0, description = 'hmin1',\n",
    "                        orientation='vertical',layout=Layout(width='40px', height='200px'),continuous_update = False)\n",
    "hmin1.style.handle_color = 'gray'\n",
    "uikvar = widgets.HBox([nsim,nug,it1,azi,hmaj1,hmin1],)     \n",
    "\n",
    "# dashboard: data locations \n",
    "x1 = widgets.FloatSlider(min=0.0, max = 1000.0, value = 100.0, step = 1.0, description = 'x1',orientation='horizontal',\n",
    "                         layout=Layout(width='180px', height='30px'),readout_format = '.0f',style=style,continuous_update = False)\n",
    "x1.style.handle_color = 'darkblue'\n",
    "y1 = widgets.FloatSlider(min=0.0, max = 1000.0, value = 100.0, step = 1.0, description = 'y1',orientation='vertical',\n",
    "                         layout=Layout(width='90px', height='180px'),readout_format = '.0f',style=style,continuous_update = False)\n",
    "y1.style.handle_color = 'darkblue'\n",
    "uik1 = widgets.VBox([x1,y1],)\n",
    "\n",
    "x2 = widgets.FloatSlider(min=0.0, max = 1000.0, value = 500.0, step = 1.0, description = 'x2',orientation='horizontal',\n",
    "                         layout=Layout(width='180px', height='30px'),readout_format = '.0f',style=style,continuous_update = False)\n",
    "x2.style.handle_color = 'darkred'\n",
    "y2 = widgets.FloatSlider(min=0.0, max = 1000.0, value = 800.0, step = 1.0, description = 'y2',orientation='vertical',\n",
    "                         layout=Layout(width='90px', height='180px'),readout_format = '.0f',style=style,continuous_update = False)\n",
    "y2.style.handle_color = 'darkred'\n",
    "uik2 = widgets.VBox([x2,y2],)\n",
    "\n",
    "x3 = widgets.FloatSlider(min=0.0, max = 1000.0, value = 900.0, step = 1.0, description = 'x3',orientation='horizontal',\n",
    "                         layout=Layout(width='180px', height='30px'),readout_format = '.0f',style=style,continuous_update = False)\n",
    "x3.style.handle_color = 'darkorange'\n",
    "y3 = widgets.FloatSlider(min=0.0, max = 1000.0, value = 200.0, step = 1.0, description = 'y3',orientation='vertical',\n",
    "                         layout=Layout(width='90px', height='180px'),readout_format = '.0f',style=style,continuous_update = False)\n",
    "y3.style.handle_color = 'darkorange'\n",
    "uik3 = widgets.VBox([x3,y3],)\n",
    "\n",
    "uipars = widgets.HBox([uikvar,uik1,uik2,uik3],) \n",
    "uik = widgets.VBox([l,uipars],)\n",
    "\n",
    "def convert_type(it):\n",
    "    if it == 'Spherical': \n",
    "        return 1\n",
    "    elif it == 'Exponential':\n",
    "        return 2\n",
    "    else: \n",
    "        return 3\n",
    "\n",
    "def f_make_krige(nsim,nug,it1,azi,hmaj1,hmin1,x1,y1,x2,y2,x3,y3): # function to take parameters, make sample and plot\n",
    "    text_trap = io.StringIO()                           # suppress all text function output to dashboard to avoid clutter \n",
    "    sys.stdout = text_trap\n",
    "    cmap = cm.inferno\n",
    "    np.random.seed(seed = 73073)                        # ensure same results for all runs\n",
    "    it1 = convert_type(it1)\n",
    "    nst = 1; xlag = 10; nlag = int(hmaj1/xlag); c1 = 1.0-nug\n",
    "    vario = GSLIB.make_variogram(nug,nst,it1,c1,azi,hmaj1,hmin1) # make model object\n",
    "    index_maj,h_maj,gam_maj,cov_maj,ro_maj = geostats.vmodel(nlag,xlag,azi,vario)   # project the model in the major azimuth                                                  # project the model in the 135 azimuth\n",
    "    index_min,h_min,gam_min,cov_min,ro_min = geostats.vmodel(nlag,xlag,azi+90.0,vario) # project the model in the minor azimuth\n",
    "    \n",
    "    seed = 73073\n",
    "\n",
    "# make hard data dataframe and hard code the data values\n",
    "    x = [x1,x2,x3]; y = [y1,y2,y3]; value = [-2.0,0.0,2.0] \n",
    "    df = pd.DataFrame({'X':x,'Y':y,'Value':value})\n",
    "    ndata = len(df); skmean = np.average(df['Value'].values)\n",
    "\n",
    "# make simulation locations dataframe\n",
    "    random.seed(a = seed)\n",
    "    xl = random.sample(range(0, 1000), nsim); \n",
    "    random.seed(a = seed+1)\n",
    "    yl = random.sample(range(0, 1000), nsim); valuel = np.full(nsim,-9999)\n",
    "    dfl = pd.DataFrame({'X':xl,'Y':yl, 'Value':valuel},dtype=np.single)\n",
    "    dfl_temp = pd.DataFrame({'X':[-9999,9999],'Y':[-9999,9999], 'Value':[-9999,-9999]},dtype=np.single)\n",
    "  \n",
    "    sim = np.zeros(len(dfl)); sk_est = np.zeros(len(dfl)); sk_var = np.zeros(len(dfl)); sk_std = np.zeros(len(dfl))\n",
    "    sk_weights = np.zeros([ndata,len(dfl)])\n",
    "     \n",
    "# perform sequential simulation\n",
    "    for isim in range(0,len(dfl)):\n",
    "        dfl_temp.loc[0,'X'] = dfl.loc[isim,'X']; dfl_temp.loc[0,'Y'] = dfl.loc[isim,'Y']; # copy current data to first data / method needs atleast 2 data\n",
    "        sk_est_temp, sk_var_temp, sk_weights_temp = simple_simple_krige(df,'X','Y','Value',dfl_temp,'X','Y',vario,skmean=skmean)\n",
    "        sk_est[isim] = sk_est_temp[0]; \n",
    "        sk_var[isim] = sk_var_temp[0]; \n",
    "        sk_weights[:,isim] = sk_weights_temp[0,:ndata]\n",
    "        if sk_var[isim] == 0: \n",
    "            sk_std[isim] = 0.0\n",
    "        else:\n",
    "            sk_std[isim] = math.sqrt(sk_var[isim])\n",
    "        sim[isim] = norm.rvs(loc=sk_est[isim], scale=sk_std[isim], size=1)[0]  # random seedset at the start \n",
    "        #df = df.append({'X': dfl.loc[isim,'X'],'Y': dfl.loc[isim,'Y'],'Value': sim[isim]}, ignore_index=True) # append is removed\n",
    "        df = pd.concat([df, pd.DataFrame.from_records([{'X': dfl.loc[isim,'X'],'Y': dfl.loc[isim,'Y'],'Value': sim[isim]}])]\n",
    "            ,ignore_index=True)\n",
    "        dfl.at[isim,'Value'] = float(sim[isim])\n",
    "    \n",
    "# plot the variogram model\n",
    "    xlag = 10.0; nlag = int(hmaj1/xlag)\n",
    "    plt.subplot(1,3,1)\n",
    "    plt.plot([0,hmaj1*1.5],[1.0,1.0],color = 'black')\n",
    "    plt.plot(h_maj,gam_maj,color = 'black',label = 'Major ' + str(azi))    \n",
    "    plt.plot(h_min,gam_min,color = 'black',label = 'Minor ' + str(azi+90.0))\n",
    "    deltas = [22.5, 45, 67.5]; \n",
    "    ndelta = len(deltas); hd = np.zeros(ndelta); gamd = np.zeros(ndelta);\n",
    "    color=iter(cm.plasma(np.linspace(0,1,ndelta)))\n",
    "    for delta in deltas:\n",
    "        index,hd,gamd,cov,ro = geostats.vmodel(nlag,xlag,azi+delta,vario);\n",
    "        c=next(color)\n",
    "        plt.plot(hd,gamd,color = c,label = 'Azimuth ' + str(azi+delta))\n",
    "    plt.xlabel(r'Lag Distance $\\bf(h)$, (m)')\n",
    "    plt.ylabel(r'$\\gamma \\bf(h)$')\n",
    "    plt.title('Interpolated NSCORE Porosity Variogram Models')\n",
    "    plt.xlim([0,hmaj1*1.5])\n",
    "    plt.ylim([0,1.4])\n",
    "    plt.legend(loc='upper left')\n",
    "\n",
    "# plot the data and simulated values on a scatter plot \n",
    "    sk_weights_avg = np.mean(sk_weights,axis = 1)\n",
    "    plt.subplot(1,3,2)\n",
    "    for idata in range(0,len(df)):\n",
    "        if idata < ndata:\n",
    "            plt.scatter([df.loc[idata,'X']],[df.loc[idata,'Y']],marker='^',\n",
    "                        c = [df.loc[idata,'Value']], cmap = cmap, vmin = -2.0, vmax = 2.0, edgecolors = 'black',\n",
    "                        s = 100,label = 'Original Data')\n",
    "        else: \n",
    "            plt.scatter([df.loc[idata,'X']],[df.loc[idata,'Y']],\n",
    "                        c = [df.loc[idata,'Value']], cmap = cmap, vmin = -2.0, vmax = 2.0, edgecolors = 'black',\n",
    "                        label = 'Simulated Values')\n",
    "    \n",
    "    ax = plt.gca()\n",
    "    plt.xlabel('X(m)'); plt.ylabel('Y(m)')\n",
    "    plt.title('Sequential Simulation - Data and Unknown Locations')\n",
    "    plt.xlim([0,1000])\n",
    "    plt.ylim([0,1000])\n",
    "    plt.colorbar()\n",
    "\n",
    "    if nsim < 10:\n",
    "        for i, txt in enumerate(np.round(dfl['Value'].values,2)):\n",
    "            plt.annotate(txt, (dfl.loc[i,'X']-40, dfl.loc[i,'Y']-40))        \n",
    "\n",
    "    ellipse = Ellipse((500, 500),width=hmin1*2.0,height=hmaj1*2.0,angle = 360-azi,facecolor='none',edgecolor='black',alpha = 1.0,linestyle='--')\n",
    "    ax = plt.gca()\n",
    "    ax.add_patch(ellipse)\n",
    "\n",
    "# plot the distribution of the simulated values   \n",
    "    plt.subplot(1,3,3)\n",
    "    plt.hist(sim,bins = np.linspace(-3.0,3.0,20),alpha=0.2,color=\"red\",edgecolor=\"black\")\n",
    "    plt.xlim([-3.0,3.0]); plt.ylim([0,nsim/2])\n",
    "    plt.title('Uncertainty Model at Unknown Location')\n",
    "    plt.xlabel('Value'); plt.ylabel('Frequency')\n",
    "    \n",
    "    ax = plt.gca()\n",
    "    ax.annotate('Simulations: Mean = ' + str(np.round(np.average(sim),2)), (-2.8, nsim*0.05))\n",
    "    ax.annotate('Simulations: Standard Deviation = ' + str(np.round(np.std(sim),2)), (-2.8, nsim *0.02))\n",
    "    plt.subplots_adjust(left=0.0, bottom=0.0, right=2.2, top=0.9, wspace=0.3, hspace=0.3)\n",
    "    plt.show()\n",
    "    \n",
    "# connect the function to make the samples and plot to the widgets    \n",
    "interactive_plot = widgets.interactive_output(f_make_krige, {'nsim':nsim,'nug':nug, 'it1':it1, 'azi':azi, 'hmaj1':hmaj1, 'hmin1':hmin1, \n",
    "                                                      'x1':x1, 'y1':y1, 'x2':x2, 'y2':y2, 'x3':x3, 'y3':y3,})\n",
    "interactive_plot.clear_output(wait = True)               # reduce flickering by delaying plot updating"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interactive Sequential Simulation to Random Points Demostration\n",
    "\n",
    "* select the variogram model and the data locations and observe the outputs from sequential simulation \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) | [GeostatsPy](https://github.com/GeostatsGuy/GeostatsPy)\n",
    "\n",
    "### The Inputs\n",
    "\n",
    "Select the variogram model and the data locations:\n",
    "\n",
    "* **nug**: nugget effect\n",
    "\n",
    "* **c1 **: contributions of the sill\n",
    "\n",
    "* **hmaj1 / hmin1 **: range in the major and minor direction\n",
    "\n",
    "* **(x1, y1),...(x3,y3) **: spatial data locations  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "23d6a33e436f49de9f3a5cdb56d7c793",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='                                              Sequential Simulation, Michael Pyrcz,…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5f131774884e4df68a4b2834f0801aa6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output(outputs=({'output_type': 'display_data', 'data': {'text/plain': '<Figure size 640x480 with 4 Axes>', 'i…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(uik, interactive_plot)                            # display the interactive plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Interactive Sequential Simulation to a Regular Grid Method\n",
    "\n",
    "Let's repeat the previous interactive demonstration, but this time we will simulate on a random set of nodes on a regular grid with our point data.\n",
    "\n",
    "* this is more similar to current practice with most spatial modeling software\n",
    "\n",
    "The following code includes:\n",
    "\n",
    "* dashboard with number of simulated nodes, variogram model and data locations \n",
    "\n",
    "* plot of the point data and the simulated model on a regular grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings; warnings.simplefilter('ignore')\n",
    "\n",
    "# dashboard: number of simulation grid nodes and the variogram model\n",
    "style = {'description_width': 'initial'}\n",
    "l = widgets.Text(value='                                              Sequential Simulation, Michael Pyrcz, Associate Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n",
    "nsim = widgets.IntSlider(min = 0, max = 100, value = 5, step = 1, description = 'nsim',orientation='vertical',\n",
    "                          layout=Layout(width='40px', height='200px'),continuous_update=False)\n",
    "nsim.style.handle_color = 'gray'\n",
    "nug = widgets.FloatSlider(min = 0, max = 1.0, value = 0.0, step = 0.1, description = 'nug',orientation='vertical',\n",
    "                          layout=Layout(width='25px', height='200px'),continuous_update=False)\n",
    "nug.style.handle_color = 'gray'\n",
    "it1 = widgets.Dropdown(options=['Spherical', 'Exponential', 'Gaussian'],value='Spherical',\n",
    "    description='Type1:',disabled=False,layout=Layout(width='180px', height='30px'), style=style)\n",
    "\n",
    "seed = widgets.IntText(value=73074,description='Seed:',disabled=False,layout=Layout(width='180px', height='30px'),continuous_update=False)\n",
    "\n",
    "azi = widgets.FloatSlider(min=0, max = 360, value = 0, step = 22.5, description = 'azi',\n",
    "                        orientation='vertical',layout=Layout(width='40px', height='200px'),continuous_update=False)\n",
    "azi.style.handle_color = 'gray'\n",
    "hmaj1 = widgets.FloatSlider(min=0.01, max = 10000.0, value = 100.0, step = 25.0, description = 'hmaj1',\n",
    "                        orientation='vertical',layout=Layout(width='40px', height='200px'),continuous_update=False)\n",
    "hmaj1.style.handle_color = 'gray'\n",
    "hmin1 = widgets.FloatSlider(min = 0.01, max = 10000.0, value = 100.0, step = 25.0, description = 'hmin1',\n",
    "                        orientation='vertical',layout=Layout(width='40px', height='200px'),continuous_update=False)\n",
    "hmin1.style.handle_color = 'gray'\n",
    "\n",
    "uikvarb = widgets.VBox([it1,seed],)\n",
    "\n",
    "uikvar = widgets.HBox([nsim,nug,uikvarb,azi,hmaj1,hmin1],)                   # basic widget formatting   \n",
    "\n",
    "# dashboard: data locations\n",
    "x1 = widgets.FloatSlider(min=0.0, max = 1000.0, value = 100.0, step = 1.0, description = 'x1',orientation='horizontal',\n",
    "                         layout=Layout(width='180px', height='30px'),readout_format = '.0f',style=style,continuous_update=False)\n",
    "x1.style.handle_color = 'blue'\n",
    "y1 = widgets.FloatSlider(min=0.0, max = 1000.0, value = 100.0, step = 1.0, description = 'y1',orientation='vertical',\n",
    "                         layout=Layout(width='90px', height='180px'),readout_format = '.0f',style=style,continuous_update=False)\n",
    "y1.style.handle_color = 'blue'\n",
    "uik1 = widgets.VBox([x1,y1],)\n",
    "\n",
    "x2 = widgets.FloatSlider(min=0.0, max = 1000.0, value = 500.0, step = 1.0, description = 'x2',orientation='horizontal',\n",
    "                         layout=Layout(width='180px', height='30px'),readout_format = '.0f',style=style,continuous_update=False)\n",
    "x2.style.handle_color = 'red'\n",
    "y2 = widgets.FloatSlider(min=0.0, max = 1000.0, value = 800.0, step = 1.0, description = 'y2',orientation='vertical',\n",
    "                         layout=Layout(width='90px', height='180px'),readout_format = '.0f',style=style,continuous_update=False)\n",
    "y2.style.handle_color = 'red'\n",
    "uik2 = widgets.VBox([x2,y2],)\n",
    "\n",
    "x3 = widgets.FloatSlider(min=0.0, max = 1000.0, value = 900.0, step = 1.0, description = 'x3',orientation='horizontal',\n",
    "                         layout=Layout(width='180px', height='30px'),readout_format = '.0f',style=style,continuous_update=False)\n",
    "x3.style.handle_color = 'darkorange'\n",
    "y3 = widgets.FloatSlider(min=0.0, max = 1000.0, value = 200.0, step = 1.0, description = 'y3',orientation='vertical',\n",
    "                         layout=Layout(width='90px', height='180px'),readout_format = '.0f',style=style,continuous_update=False)\n",
    "y3.style.handle_color = 'darkorange'\n",
    "uik3 = widgets.VBox([x3,y3],)\n",
    "\n",
    "uipars = widgets.HBox([uikvar,uik1,uik2,uik3],) \n",
    "uik = widgets.VBox([l,uipars],)\n",
    "\n",
    "def convert_type(it):\n",
    "    if it == 'Spherical': \n",
    "        return 1\n",
    "    elif it == 'Exponential':\n",
    "        return 2\n",
    "    else: \n",
    "        return 3\n",
    "\n",
    "def f_make_krige2(nsim,nug,it1,seed,azi,hmaj1,hmin1,x1,y1,x2,y2,x3,y3): # function to take parameters, make sample and plot\n",
    "    text_trap = io.StringIO()                           # suppress all text function output to dashboard to avoid clutter \n",
    "    sys.stdout = text_trap\n",
    "    cmap = cm.inferno                                                   \n",
    "    it1 = convert_type(it1)\n",
    "    nst = 1; xlag = 10; nlag = int(hmaj1/xlag); c1 = 1.0-nug\n",
    "    vario = GSLIB.make_variogram(nug,nst,it1,c1,azi,hmaj1,hmin1) # make model object\n",
    "    index_maj,h_maj,gam_maj,cov_maj,ro_maj = geostats.vmodel(nlag,xlag,azi,vario)   # project the model in the major azimuth                                                  # project the model in the 135 azimuth\n",
    "    index_min,h_min,gam_min,cov_min,ro_min = geostats.vmodel(nlag,xlag,azi+90.0,vario) # project the model in the minor azimuth\n",
    "    \n",
    "# make data dataframe\n",
    "    x = [x1,x2,x3]; y = [y1,y2,y3]; value = [-1.5,0.0,1.5] \n",
    "    df = pd.DataFrame({'X':x,'Y':y,'Value':value})\n",
    "    ndata = len(df); skmean = np.average(df['Value'].values)\n",
    "\n",
    "# make simulation nodes dataframe   \n",
    "    random.seed(a = seed)                               # ensure same results for all runs, you can sequentially add / remove nodes\n",
    "    if nsim == 100:\n",
    "        icelll = np.linspace(0, nx*ny-1, 100) \n",
    "        random.shuffle(icelll)\n",
    "    else:\n",
    "        random.seed(seed)\n",
    "        icelll = np.asarray(random.sample(range(0, nx*ny-1), nsim),dtype = np.int32)\n",
    "    iyl = np.around(icelll / nx-0.49,0); yl = iyl * csiz + ymn\n",
    "    ixl = np.around(icelll - iyl * nx , 0); xl = ixl * csiz + xmn\n",
    "    valuel = np.full(nsim,-9999)\n",
    "    dfl = pd.DataFrame({'X':xl,'Y':yl, 'Value':valuel},dtype=np.single)\n",
    "    dfl_temp = pd.DataFrame({'X':[-9999,9999],'Y':[-9999,9999], 'Value':[-9999,-9999]},dtype=np.single)\n",
    "  \n",
    "    np.random.seed(seed = seed)\n",
    "    sim = np.zeros(len(dfl)); sk_est = np.zeros(len(dfl)); sk_var = np.zeros(len(dfl)); sk_std = np.zeros(len(dfl))\n",
    "    sk_weights = np.zeros([ndata,len(dfl)])\n",
    "\n",
    "# perform sequential simulation\n",
    "    for isim in range(0,len(dfl)):\n",
    "        dfl_temp.loc[0,'X'] = dfl.loc[isim,'X']; dfl_temp.loc[0,'Y'] = dfl.loc[isim,'Y']; # copy current data to first data / method needs atleast 2 data\n",
    "        sk_est_temp, sk_var_temp, sk_weights_temp = simple_simple_krige(df,'X','Y','Value',dfl_temp,'X','Y',vario,skmean=skmean)\n",
    "        sk_est[isim] = sk_est_temp[0]; \n",
    "        sk_var[isim] = sk_var_temp[0]; \n",
    "        sk_weights[:,isim] = sk_weights_temp[0,:ndata]\n",
    "        if sk_var[isim] == 0: \n",
    "            sk_std[isim] <= 0.0\n",
    "        else:\n",
    "            sk_std[isim] = math.sqrt(sk_var[isim])\n",
    "        sim[isim] = norm.rvs(loc=sk_est[isim], scale=sk_std[isim], size=1)[0]  # random seedset at the start \n",
    "  \n",
    "        #df = df.append({'X': dfl.loc[isim,'X'],'Y': dfl.loc[isim,'Y'],'Value': sim[isim]}, ignore_index=True) # append has been removed\n",
    "        df = pd.concat([df, pd.DataFrame.from_records([{'X': dfl.loc[isim,'X'],'Y': dfl.loc[isim,'Y'],'Value': sim[isim]}])]\n",
    "                      ,ignore_index=True)\n",
    "    \n",
    "        dfl.at[isim,'Value'] = float(sim[isim])\n",
    "\n",
    "# make the 2D simulated model on a regular grid        \n",
    "    plt.subplot(121)\n",
    "    model = np.full([ny,nx],-999.9)            \n",
    "    for idata in range(len(df)-1,-1,-1):\n",
    "        ix = int(df.loc[idata,'X']/csiz); iy = int(df.loc[idata,'Y']/csiz);\n",
    "        model[ny - iy - 1, ix] = df.loc[idata,'Value']\n",
    "        \n",
    "    ax = plt.gca()\n",
    "    plt.xlabel('X(m)'); plt.ylabel('Y(m)')\n",
    "    plt.title('Sequential Simulation - Data, Simulated Values and Random Path')  \n",
    "    palette = copy(plt.cm.inferno)\n",
    "    palette.set_under('r', 0.0)\n",
    "    palette.set_over('r', 0.0)\n",
    "    im = plt.imshow(model,interpolation = None,extent = [0,1000,0,1000], vmin = -3.0, vmax = 3.0,cmap = palette)\n",
    "    plt.scatter(df['X'].values[:ndata],df['Y'].values[:ndata],marker='^',c=df['Value'].values[:ndata], vmin = -2.0, vmax = 2.0, cmap = cmap, edgecolors = 'black',s = 300,label = 'Original Data')\n",
    "    plt.xlim([0,1000]); plt.ylim([0,1000])\n",
    "    for idata in range(len(df)-1,-1,-1):\n",
    "        x = df.loc[idata,'X'];y = df.loc[idata,'Y']\n",
    "        ix = int(x/csiz); iy = int(y/csiz)\n",
    "        xc = csiz*ix + csiz*0.45; yc = csiz*iy + csiz*0.5; \n",
    "#         if idata < 3:\n",
    "#             #plt.annotate('D'+str(idata+1),[xc-15,yc],color='white',weight='bold')\n",
    "#         else:\n",
    "        if idata > 2:\n",
    "            plt.annotate(idata-2,[xc-10,yc],color='white')\n",
    "\n",
    "    cbar = plt.colorbar(im,ax = plt.gca()) # Similar to fig.colorbar(im, cax = cax)\n",
    "    plt.gca().set_aspect('auto')\n",
    "    cbar.set_label('Simulated Values')\n",
    "    \n",
    "# plot the variogram modle for visualization    \n",
    "    ellipse1 = Ellipse((x1, y1),width=hmin1*2.0,height=hmaj1*2.0,angle = 360-azi,facecolor='none',edgecolor='darkblue',lw=3,alpha = 0.6)\n",
    "    ellipse2 = Ellipse((x2, y2),width=hmin1*2.0,height=hmaj1*2.0,angle = 360-azi,facecolor='none',edgecolor='darkred',lw=3,alpha = 0.6)\n",
    "    ellipse3 = Ellipse((x3, y3),width=hmin1*2.0,height=hmaj1*2.0,angle = 360-azi,facecolor='none',edgecolor='darkorange',lw=3,alpha = 0.6)\n",
    "    ax = plt.gca()\n",
    "    ax.add_patch(ellipse1); ax.add_patch(ellipse2); ax.add_patch(ellipse3);\n",
    "    \n",
    "    x_values = np.linspace(-3.0,3.0,100)                                          # get an array of x values\n",
    "    p_values = norm.pdf(x_values, loc = 0.0, scale = 1.0)    \n",
    "    \n",
    "    plt.subplot(122)\n",
    "    plt.hist(model.flatten(),color='red',alpha=0.8,edgecolor='black',bins=np.linspace(-3,3,10),density =True)\n",
    "    plt.plot(x_values,p_values,color='red')\n",
    "    plt.xlim(-3,3); plt.ylim(0,0.6)\n",
    "    plt.title('Distribution of Seuqential Gaussian Simulated Values')\n",
    "    plt.xlabel('Simulated Gaussian Values'); plt.ylabel('Normalized Frequence')\n",
    "    \n",
    "    plt.gca().annotate('Simulation Mean = ' + str(np.round(stats.tmean(model.flatten(),limits=(-5,5)),2)), (0.9, 0.55))\n",
    "    plt.gca().annotate('Simulation StDev. = ' + str(np.round(stats.tstd(model.flatten(),limits=(-3,3)),2)), (0.9, 0.52))\n",
    "    \n",
    "    plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.3, hspace=0.3)\n",
    "    plt.show()\n",
    "    \n",
    "# connect the function to make the samples and plot to the widgets    \n",
    "interactive_plot = widgets.interactive_output(f_make_krige2, {'nsim':nsim,'nug':nug, 'it1':it1,'seed':seed,'azi':azi, 'hmaj1':hmaj1, 'hmin1':hmin1, \n",
    "                                                      'x1':x1, 'y1':y1, 'x2':x2, 'y2':y2, 'x3':x3, 'y3':y3,})\n",
    "interactive_plot.clear_output(wait = True)               # reduce flickering by delaying plot updating"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interactive Sequential Simulation to Model Regular Grid Demonstration\n",
    "\n",
    "* select the variogram model and the data locations and observe the outputs from sequential simulation \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) | [GeostatsPy](https://github.com/GeostatsGuy/GeostatsPy)\n",
    "\n",
    "### The Inputs\n",
    "\n",
    "Select the simulation nodes, variogram model and the data locations:\n",
    "\n",
    "* **nsim**: number of simulated nodes, for computational speed up use less nodes\n",
    "\n",
    "* **nug**: nugget effect\n",
    "\n",
    "* **c1 **: contributions of the sill\n",
    "\n",
    "* **hmaj1 / hmin1 **: range in the major and minor direction\n",
    "\n",
    "* **(x1, y1),...(x3,y3) **: spatial data locations  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ecf1423970ad4fe3b0eb914a3c38127b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='                                              Sequential Simulation, Michael Pyrcz,…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "70416ed18def48d1a00a1af209b2dbe6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output(outputs=({'output_type': 'display_data', 'data': {'text/plain': '<Figure size 640x480 with 3 Axes>', 'i…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(uik, interactive_plot)                            # display the interactive plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comments\n",
    "\n",
    "This was an interactive demonstration of sequential Gaussian simulation for spatial data analytics. Much more could be done, I have other demonstrations on the basics of working with DataFrames, ndarrays, univariate statistics, plotting data, declustering, data transformations and many other workflows available at https://github.com/GeostatsGuy/PythonNumericalDemos and https://github.com/GeostatsGuy/GeostatsPy. \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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}