<jupyter_start><jupyter_text>scikit learn for Prediction  scikit learn for Subsurface Modeling in Python  <NAME>, Associate Professor, University of Texas at Austin  [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) PGE 383 Exercise: scikit learn for Subsurface Modeling in Python Here's a simple workflow, demonstration of scikit learn for subsurface modeling workflows. This should help you get started with building subsurface models that data analytics and machine learning. Here's some basic details about scikit learn.   scikit learnDemonstration of scikit learn for machine learning.In this workflow we demonstrate the plug and play nature of scikit learn machine learning models.  For an unconventional dataset we demonstrate the following steps:        1. instantiation    2. fitting    3. prediction    4. cross validationWe will work with the following regression methods:* linear regression* multilinear regression* decision tree regression* support vector machine regression Getting StartedHere's the steps to get setup in Python with the GeostatsPy package:1. Install Anaconda 3 on your machine (https://www.anaconda.com/download/). 2. From Anaconda Navigator (within Anaconda3 group), go to the environment tab, click on base (root) green arrow and open a terminal. 3. In the terminal type: pip install geostatspy. 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. 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.  Import Required PackagesLet's import the GeostatsPy package.<jupyter_code>import os                                                 # to set current working directory 
import math                                               # basic calculations like square root
from sklearn.model_selection import train_test_split      # train and test split
from sklearn import svm                                   # support vector machine methods
from sklearn import tree                                  # tree program from scikit learn (package for machine learning)
from sklearn.metrics import mean_squared_error, r2_score  # specific measures to check our models
import pandas as pd                                       # DataFrames and plotting
import pandas.plotting as pd_plot
import numpy as np                                        # arrays and matrix math
import matplotlib.pyplot as plt                           # plotting
from intake import cat                                    # data catalogue<jupyter_output><empty_output><jupyter_text>Declare functionsLet's define a couple of functions to streamline plotting correlation matrices and visualization of a machine learning regression model responce over the 2 predictor features.<jupyter_code>def plot_corr(dataframe,size=10):                         # plots a graphical correlation matrix 
    corr = dataframe.corr()
    fig, ax = plt.subplots(figsize=(size, size))
    im = ax.matshow(corr,vmin = -1.0, vmax = 1.0)
    plt.xticks(range(len(corr.columns)), corr.columns);
    plt.yticks(range(len(corr.columns)), corr.columns);
    plt.colorbar(im, orientation = 'vertical')
    plt.title('Correlation Matrix')
    
def visualize_model(model,xfeature,yfeature,response,title,):# plots the data points and the decision tree prediction 
    n_classes = 10
    cmap = plt.cm.RdYlBu
    plot_step = 0.02
    x_min, x_max = min(xfeature) - 1, max(xfeature) + 1
    y_min, y_max = min(yfeature) - 1, max(yfeature) + 1
    resp_min = round(min(response)); resp_max = round(max(response));
    xx, yy = np.meshgrid(np.arange(x_min, x_max, plot_step),
                     np.arange(y_min, y_max, plot_step))
    z_min = round(min(response)); z_max = round(max(response))
    Z = model.predict(np.c_[xx.ravel(), yy.ravel()])
    Z = Z.reshape(xx.shape)
    cs = plt.contourf(xx, yy, Z, cmap=cmap,vmin=z_min, vmax=z_max)

    im = plt.scatter(xfeature,yfeature,s=None, c=response, marker=None, cmap=cmap, norm=None, vmin=z_min, vmax=z_max, alpha=0.8, linewidths=0.3, verts=None, edgecolors="black")
    plt.title(title)
    plt.xlabel(xfeature.name)
    plt.ylabel(yfeature.name)
    cbar = plt.colorbar(im, orientation = 'vertical')
    cbar.set_label(response.name, rotation=270, labelpad=20)
    return(plt)<jupyter_output><empty_output><jupyter_text>Read the data table<jupyter_code>my_data = cat.unconv_MV().read()                    # load the comma delimited data file<jupyter_output><empty_output><jupyter_text>Let's visualize the first several rows of our data stored in a DataFrame so we can make sure we successfully loaded the data file.<jupyter_code>my_data.head()                                            # preview the first 5 rows of the dataframe<jupyter_output><empty_output><jupyter_text>Let's remove the well index and check the summary summary statistics.<jupyter_code>my_data = my_data.iloc[:,1:8]                             # copy all rows and columns 1 through 8, note 0 column is removed
my_data.describe().transpose()                            # calculate summary statistics for the data<jupyter_output><empty_output><jupyter_text>It is good that we checked the summary statistics, because we have some negative values for brittleness and total organic carbon. The is physically imposible.  The values must be in error. We know the lowest possible values are 0.0, so we will truncate on 0.0.  We use the *get_numerical_data()* DataFrame member function to get a shallow copy of the data from the DataFrame.  Since it is a shallow copy, any changes we make to the copy are made to the data in the original DataFrame.  This allows us to apply this simple conditional statement to all the data values in the DataFrame all at once.<jupyter_code>num = my_data._get_numeric_data()                         # get the numerical values
num[num < 0] = 0                                          # truncate negative values to 0.0
my_data.describe().transpose()                            # calculate summary statistics for the data<jupyter_output><empty_output><jupyter_text>This dataset has variables from 1,000 unconventional wells including well average porosity, log transform of permeability (to linearize the relationships with other variables), accoustic impedance (kg/m2s*10^6), brittness ratio (%), total organic carbon (%), vitrinite reflectance (%), and initial production 90 day average (MCFPD).  Note, the dataset is synthetic. Calculate the correlation matrix For multivariate analysis it is a good idea to check the correlation matrix.  We can calculate it and view it in the console with these commands.<jupyter_code>corr_matrix = np.corrcoef(my_data, rowvar = False)
print(np.around(corr_matrix,2))                           # print the correlation matrix to 2 decimals<jupyter_output>[[ 1.    0.81 -0.51 -0.25  0.71  0.08  0.69]
 [ 0.81  1.   -0.32 -0.15  0.51  0.05  0.57]
 [-0.51 -0.32  1.    0.17 -0.55  0.49 -0.33]
 [-0.25 -0.15  0.17  1.   -0.24  0.3  -0.07]
 [ 0.71  0.51 -0.55 -0.24  1.    0.31  0.5 ]
 [ 0.08  0.05  0.49  0.3   0.31  1.    0.14]
 [ 0.69  0.57 -0.33 -0.07  0.5   0.14  1.  ]]<jupyter_text>Note the 1.0 diagonal resulting from the correlation of each variable with themselves.  Let's use our function declared above to make a graphical correlation matrix visualization.  This may inprove our ability to spot features.  It relies on the built in correlation matrix method with Numpy DataFrames and MatPlotLib for plotting.<jupyter_code>plot_corr(my_data,10)                                      # using our correlation matrix visualization function
plt.show()<jupyter_output><empty_output><jupyter_text>This looks good.  There is a mix of correlation magnitudes. Of course, correlation coeffficients are limited to degree of linear correlations.  For more complete information, let's look at the matrix scatter plot from the Pandas package.<jupyter_code>pd_plot.scatter_matrix(my_data, alpha = 0.1,              # pandas matrix scatter plot
    figsize=(10, 10),color = 'black', hist_kwds={'color':['grey']})
plt.show()<jupyter_output><empty_output><jupyter_text>Working with Only Two Predictor FeaturesLet's simplify the problem to 2 features, Porosity and Brittleness, to predict one response feature, Production rate.  We will also reduce the number of wells from 1,000 to 500.  By working with only 2 predictor features, it is very easy to visualize the segmentation of the feature space (it is only 2D and the model can be interogated exhaustively on a single plot).<jupyter_code>X = my_data.iloc[0:500,[0,3]]              # extract porosity and brittleness 500 samples to a predictor array
y = my_data.iloc[0:500,[6]]                # extract production 500 samples to a response array<jupyter_output><empty_output><jupyter_text>Let's check the summary statistics of the predictor and response features.<jupyter_code>X.describe().transpose()                   # calculate summary statistics for the data
y.describe().transpose()                   # calculate summary statistics for the data<jupyter_output><empty_output><jupyter_text>Now let's withhold 100 samples as testing data and retain the remaining 400 as training data.<jupyter_code>X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=73073)
n_train = len(X_train)
n_test = len(X_test)
print('Number of training ' + str(n_train) + ', number of test ' + str(n_test))<jupyter_output>Number of training 400, number of test 100<jupyter_text>Let's compare the univariate statistics of Porosity, Brittleness and Producton training and testing datasets.    * let's check for bias and extrapolation.<jupyter_code>X_train.describe().transpose()                   # calculate summary statistics for the data
X_test.describe().transpose()                   # calculate summary statistics for the data
y_train.describe().transpose()                   # calculate summary statistics for the data
y_test.describe().transpose()                   # calculate summary statistics for the data<jupyter_output><empty_output><jupyter_text>Now let's plot the training and testing dataset distributions to check coverage and extrapolation.<jupyter_code>plt.subplot(231)
plt.hist(X_train["Por"], alpha = 0.2, color = 'red', edgecolor = 'black', bins=20)
plt.title('Porosity Train Dataset (%)')

plt.subplot(232)
plt.hist(X_train["Brittle"], alpha = 0.2, color = 'red', edgecolor = 'black', bins=20)
plt.title('Britteness Train Dataset (%)')

plt.subplot(233)
plt.hist(y_train["Production"], alpha = 0.2, color = 'red', edgecolor = 'black', bins=20)
plt.title('Production Train Dataset (MCFPD)')

plt.subplot(234)
plt.hist(X_test["Por"], alpha = 0.2, color = 'red', edgecolor = 'black', bins=20)
plt.title('Porosity Test Dataset (%)')

plt.subplot(235)
plt.hist(X_test["Brittle"], alpha = 0.2, color = 'red', edgecolor = 'black', bins=20)
plt.title('Britteness Test Dataset (%)')

plt.subplot(236)
plt.hist(y_test["Production"], alpha = 0.2, color = 'red', edgecolor = 'black', bins=20)
plt.title('Production Test Dataset (MCFPD)')

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.2, top=1.2, wspace=0.3, hspace=0.3)
plt.show()<jupyter_output><empty_output><jupyter_text>The distributions are well behaved, we cannot observe obvious gaps nor truncations.  Let's look at a scatter plot of Porosity vs. Brittleness with points colored by Production.  * Let's plot the training and testing datasets to check coverage and extrapolation in the features space.<jupyter_code>plt.subplot(121)
im = plt.scatter(X_train["Por"],X_train["Brittle"],s=None, c=y_train["Production"], marker=None, cmap=None, norm=None, vmin=None, vmax=None, alpha=0.8, linewidths=0.3, verts=None, edgecolors="black")
plt.title('Training Data: Production vs. Brittleness and Porosity'); plt.xlabel('Porosity (%)'); plt.ylabel('Brittleness (%)')
cbar = plt.colorbar(im, orientation = 'vertical'); plt.xlim(5,25); plt.ylim(0,100)
cbar.set_label("Production", rotation=270, labelpad=20)

plt.subplot(122)
im = plt.scatter(X_test["Por"],X_test["Brittle"],s=None, c=y_test["Production"], marker=None, cmap=None, norm=None, vmin=None, vmax=None, alpha=0.8, linewidths=0.3, verts=None, edgecolors="black")
plt.title('Testing Data: Production vs. Brittleness and Porosity'); plt.xlabel('Porosity (%)'); plt.ylabel('Brittleness (%)')
cbar = plt.colorbar(im, orientation = 'vertical'); plt.xlim(5,25); plt.ylim(0,100)
cbar.set_label("Production", rotation=270, labelpad=20)

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.2, top=1.2, wspace=0.3, hspace=0.2)
plt.show()<jupyter_output><empty_output><jupyter_text>Building a Linear Regression ModelLet's build our first machine learning model with scikit learn.  We will start with linear regression.  For this model we will pick one predictor feature and one response feature.<jupyter_code># Linear Regression Model with scikit learn
from sklearn import linear_model

# Step 1. Instantiate the Model 
linear_reg = linear_model.LinearRegression()

# Step 2: Fit the Data on Training Data
linear_reg.fit(X_train["Por"].values.reshape(n_train,1), y_train["Production"]) # fit model
porosity_model = np.linspace(0.0,30.0,10)

# Print the model parameters
production_model = linear_reg.predict(porosity_model.reshape(10,1)) # predict with the fit model
print('Coefficients: ', str(round(linear_reg.coef_[0],3)) + ', Intercept: ', str(round(linear_reg.intercept_,3))) 

# Plot model fit
plt.figure(figsize=(8,6))
plt.scatter(X_train["Por"].values, y_train["Production"],  color='black', s = 20, alpha = 0.3)
plt.plot(porosity_model,production_model, color='red', linewidth=1)
plt.title('Linear Regression Production from Porosity on Training'); plt.xlabel('Porosity (%)'); plt.ylabel('Production (MCFPD)')
plt.xlim(5,25)#; plt.ylim(0,1500000)
# Step 3: - Make predictions using the testing dataset
y_pred = linear_reg.predict(X_test['Por'].values.reshape(n_test,1))

# Report the goodness of fit
print('Variance explained: %.2f' % r2_score(y_test, y_pred))

# Plot testing diagnostics 
plt.subplot(121)
plt.scatter(X_test['Por'].values, y_test['Production'].values,  color='black', s = 20, alpha = 0.3)
plt.scatter(X_test['Por'], y_pred, color='blue', s = 20, alpha = 0.3)
plt.title('Linear Regression Model Testing - Production from Porosity'); plt.xlabel('Porosity (%)'); plt.ylabel('Production (MCFPD)')
plt.xlim(5,25)#; plt.ylim(0,1500000)

y_res = y_pred - y_test['Production'].values

plt.subplot(122)
plt.hist(y_res, alpha = 0.2, color = 'red', edgecolor = 'black', bins=20)
plt.title('Linear Regression Model Prediction Error - Production from Porosity'); plt.xlabel('Production Estimation Error (MCFPD)'); plt.ylabel('Frequency')
plt.xlim(-4000,4000)#; plt.ylim(0,1500000)

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.2, top=1.2, wspace=0.3, hspace=0.2)
plt.show()<jupyter_output>Variance explained: 0.46<jupyter_text>Building a Multilinear Regression ModelLet's build our second machine learning model with scikit learn.  We will work with multilinear regression!  We will use both predictor features, porosity and brittleness, and one response feature, production.<jupyter_code># Linear Regression Model with scikit learn
from sklearn import linear_model

# Step 1. Instantiate the Model 
multilinear_reg = linear_model.LinearRegression()

# Step 2: Fit the Data on Training Data
multilinear_reg.fit(X_train.values.reshape(n_train,2), y_train["Production"]) # fit model

# Print the model parameters
print('Porosity Coef: ', str(round(multilinear_reg.coef_[0],3)) + ', Brittleness Coef: ', str(round(multilinear_reg.coef_[1],3)) + ', Intercept: ', str(round(multilinear_reg.intercept_,3))) 

# Plot model fit
plt.subplot(111)
plt = visualize_model(multilinear_reg,X_test["Por"],X_test["Brittle"],y_test["Production"],'Training Data and Multilinear Regression Model')
plt.subplots_adjust(left=0.0, bottom=0.0, right=1.2, top=1.2, wspace=0.3, hspace=0.2)
plt.show()

# Step 3: - Make predictions using the testing dataset
y_pred = multilinear_reg.predict(X_test.values.reshape(n_test,2))

# Report the goodness of fit
print('Variance explained: %.2f' % r2_score(y_test, y_pred))

y_res = y_pred - y_test['Production'].values

plt.subplot(111)
plt.hist(y_res, alpha = 0.2, color = 'red', edgecolor = 'black', bins=20)
plt.title('Linear Regression Model Prediction Error - Production from Porosity and Brittleness'); plt.xlabel('Production Estimation Error (MCFPD)'); plt.ylabel('Frequency')
plt.xlim(-4000,4000)#; plt.ylim(0,1500000)

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.2, top=1.2, wspace=0.3, hspace=0.2)
plt.show()<jupyter_output>Variance explained: 0.47<jupyter_text>Including brittleness only resulted in a slight improvement.  * due to the nonlinear nature of brittleness Building a Decision Tree Regression ModelLet's build our third machine learning model with scikit learn.  We will work with a decision tree.  We will use both predictor features, porosity and brittleness, and one response feature, production.<jupyter_code># Decison Tree Model with scikit learn
from sklearn.tree import tree                                     # for accessing tree information

# Step 1. Instantiate the Model 
decision_tree_reg = tree.DecisionTreeRegressor(min_samples_leaf=5, max_depth = 5)

# Step 2: Fit the Data on Training Data
decision_tree_reg.fit(X_train.values.reshape(n_train,2), y_train["Production"]) # fit model

# Plot model fit
plt.subplot(111)
plt = visualize_model(decision_tree_reg,X_test["Por"],X_test["Brittle"],y_test["Production"],'Training Data and Decision Tree Model')
plt.subplots_adjust(left=0.0, bottom=0.0, right=1.2, top=1.2, wspace=0.3, hspace=0.2)
plt.show()
# Step 3: - Make predictions using the testing dataset
y_pred = decision_tree_reg.predict(X_test.values.reshape(n_test,2))

# Report the goodness of fit
print('Variance explained: %.2f' % r2_score(y_test, y_pred))

y_res = y_pred - y_test['Production'].values

plt.subplot(111)
plt.hist(y_res, alpha = 0.2, color = 'red', edgecolor = 'black', bins=20)
plt.title('Decision Tree Model Prediction Error - Production from Porosity and Brittleness'); plt.xlabel('Production Estimation Error (MCFPD)'); plt.ylabel('Frequency')
plt.xlim(-4000,4000)#; plt.ylim(0,1500000)

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.2, top=1.2, wspace=0.3, hspace=0.2)
plt.show()<jupyter_output>Variance explained: 0.88<jupyter_text>Building a Support Vector Machine Regression ModelLet's build our fourth machine learning model with scikit learn.  We will work with a support vector machine!  We will use both predictor features, porosity and brittleness, and one response feature, production.<jupyter_code># Support Vector Regression Model with scikit learn
from sklearn import linear_model

# Step 1. Instantiate the Model 
support_vector_reg = svm.SVR(kernel='poly', C=0.01, gamma='auto', degree=2, epsilon=.01,coef0=1,max_iter=1000)

# Step 2: Fit the Data on Training Data
support_vector_reg.fit(X_train.values.reshape(n_train,2), y_train["Production"]) # fit model

# Plot model fit
plt.subplot(111)
plt = visualize_model(support_vector_reg,X_test["Por"],X_test["Brittle"],y_test["Production"],'Training Data and Support Vector Model')
plt.subplots_adjust(left=0.0, bottom=0.0, right=1.2, top=1.2, wspace=0.3, hspace=0.2)
plt.show()

# Step 3: - Make predictions using the testing dataset
y_pred = support_vector_reg.predict(X_test.values.reshape(n_test,2))

# Report the goodness of fit
print('Variance explained: %.2f' % r2_score(y_test, y_pred))

y_res = y_pred - y_test['Production'].values

plt.subplot(111)
plt.hist(y_res, alpha = 0.2, color = 'red', edgecolor = 'black', bins=20)
plt.title('Support Vector Machine Model Prediction Error - Production from Porosity and Brittleness'); plt.xlabel('Production Estimation Error (MCFPD)'); plt.ylabel('Frequency')
plt.xlim(-4000,4000)#; plt.ylim(0,1500000)

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.2, top=1.2, wspace=0.3, hspace=0.2)
plt.show()<jupyter_output>Variance explained: 0.65
StarcoderdataJupyterStructuredCleanDedup