<jupyter_start><jupyter_text>Subsurface Data Analytics  Isotonic Regression for Subsurface Data Analytics 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)  PGE 383 Exercise: Isotonic Regression for Subsurface Modeling in Python Here's a simple workflow, demonstration of isotonic regression for subsurface modeling workflows. This should help you get started with building subsurface models with data analytics and machine learning. Here's some basic details about isotonic regression.  Isotonic RegressionIsotonic regression for prediction.  Here are some key aspects of linear regression:**Stepwise Linear Parametric Model*** the fit model is a piece-wise linear model at specific predictor feature isotonic constraints / thresholds under the constraint that $x_1 \ge ... \ge x_k$ where the resulting model must be monotonically increasing, $\hat{f}(x_1) \ge ... \ge \hat{f}(x_k)$.* the loss function is: \begin{equation}min \sum_{i=1}^{n} w_i (y_i - \hat{f}(x)_i)^2\end{equation}Subject to the constraint that:\begin{equation}\hat{f}(x_1) \ge ... \ge \hat{f}(x_n)\end{equation}**Assumptions*** **Model Parameters** - parameterized by response at each isotonic constraint / threshold* **Monotonic Increasing** - may not have a negative slope * **Piece-wise Linear Model** - assumes linear interpolation between each point Workflow GoalsLearn the basics of isotonic regression in Python for analysis, modeling and prediction of porosity from density. This includes:* Basic Python workflows and data preparation* Training / fitting a ridge regression model* Checking the model and learning about the impact of hyperparameters Objective Content to support the PGE 383: Subsurface Machine Learning class.   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 numpy as np                                          # arrays and matrix math
import scipy.stats as st                                    # statistical methods
import pandas as pd                                         # DataFrames
import matplotlib.pyplot as plt                             # for plotting
from sklearn.isotonic import IsotonicRegression
from sklearn.model_selection import train_test_split        # train and test split
from sklearn.metrics import explained_variance_score<jupyter_output><empty_output><jupyter_text>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.    Set the working directoryI 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).  Also, in this case make sure to place the required (see below) data file in this working directory.<jupyter_code>os.chdir("c:/PGE383")                                       # set the working directory<jupyter_output><empty_output><jupyter_text>Loading DataLet's load the provided dataset. 'GrainSize_por.csv' is available at https://github.com/GeostatsGuy/GeoDataSets. It is a comma delimited file with 20 grain size ($mm$) and porosity (as a percentage) measures from the subsurface. We load the data file with the pandas 'read_csv' function into a data frame we called 'df'.  We preview each with the head function from Pandas DataFrames.<jupyter_code>df = pd.read_csv("GrainSize_Por.csv")                    # read a .csv file in as a DataFrame
df<jupyter_output><empty_output><jupyter_text>Test and Train SplitLet's separate the data into test and train and visualize it.<jupyter_code>X_train, X_test, y_train, y_test = train_test_split(df.iloc[:,0], df.iloc[:,1], test_size=0.25, random_state=73073)
plt.subplot(111)
plt.scatter(X_train, y_train,  color='black', s = 20, alpha = 0.5, label='train')
plt.scatter(X_test, y_test,  color='blue', s = 20, alpha = 0.5, label='test')
plt.title('Porosity from Grainsize with Training Data'); plt.xlabel('Grain Size (mm)'); plt.ylabel('Porosity (%)')
plt.xlim(0,100)#; plt.ylim(0,1500000)
plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.0, wspace=0.2, hspace=0.2)
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>The data is generally monotonic increasing with noise. The training and testing split looks reasonable.Let's assign the isotonic constraints / thresholds and visualize them.<jupyter_code>X_predict = np.linspace(5,95,10)
plt.subplot(111)
plt.scatter(X_train, y_train,  color='black', s = 20, alpha = 0.5, label='train')
plt.scatter(X_test, y_test,  color='blue', s = 20, alpha = 0.5, label='test')
for threshold in X_predict:
    plt.axvline(threshold, color='red', alpha = 0.3)
plt.title('Porosity from Grainsize'); plt.xlabel('Grain Size (mm)'); plt.ylabel('Porosity (%)')
plt.xlim(0,100)#; plt.ylim(0,1500000)
plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.0, wspace=0.2, hspace=0.2)
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Fit the ModelLet's fit the model with these isotonic constraints / thresholds.<jupyter_code>iso_reg = IsotonicRegression().fit(X_train, y_train)
y_predict = iso_reg.predict(X_predict)<jupyter_output><empty_output><jupyter_text>We can visualize the model parameters, the $f(x_1),...,f(x_k)$ estimates at the thresholds.<jupyter_code>y_predict<jupyter_output><empty_output><jupyter_text>Now we can visualize the model with a regular line plot over the model parameters.<jupyter_code># Plot model fit
plt.subplot(111)
plt.scatter(X_train, y_train,  color='black', s = 20, alpha = 0.5, label='train')
plt.scatter(X_test, y_test,  color='blue', s = 20, alpha = 0.5, label='test')
plt.plot(X_predict,y_predict, color='red', linewidth=1, label='model')
for threshold in X_predict:
    plt.axvline(threshold, color='red', alpha = 0.2)
plt.title('Isotonic Regression Porosity from Grainsize with Training Data'); plt.xlabel('Grain Size (mm)'); plt.ylabel('Porosity (%)')
plt.xlim(0,100)#; plt.ylim(0,1500000)
plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.0, wspace=0.2, hspace=0.2)
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Make Predictions and Check the Model GoodnessNow let's make predictions at the testing data locations and check the model goodness* let's use the variance explainedWe will plot the estimates for the training data as 'predicted test'* note the iso_reg.transform() function does not seem to work properly* we will just use a linear interpolation function from numpy to make estimates of the response feature with the model at new predictor feature values<jupyter_code># The transform function should calculate the result for new observations, but doesn't seem to work
# let's just use piece-wise linear interpolation from NumPy.
#y_predict_test = iso_reg.transform(X_test)
y_predict_test = np.interp(X_test, X_predict, y_predict, left=None, right=None, period=None)
r2 = explained_variance_score(y_test,y_predict_test)
print('Variance Explained = ' + str(round(r2,3)))

# Plot test predictions
plt.subplot(111)
plt.scatter(X_train, y_train,  color='black', s = 20, alpha = 0.5, label='train')
plt.scatter(X_test, y_test,  color='blue', s = 20, alpha = 0.5, label='test')
plt.scatter(X_test, y_predict_test,  color='green', s = 20, alpha = 0.5, label='predicted test')
plt.plot(X_predict,y_predict, color='red', linewidth=1, label='model')
for threshold in X_predict:
    plt.axvline(threshold, color='red', alpha = 0.2)
plt.title('Isotonic Regression Porosity from Grainsize with Training Data'); plt.xlabel('Grain Size (mm)'); plt.ylabel('Porosity (%)')
plt.xlim(0,100)#; plt.ylim(0,1500000)
plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.0, wspace=0.2, hspace=0.2)
plt.legend()
plt.show()<jupyter_output>Variance Explained = 0.337<jupyter_text>Here's another model with very few isotonic constraints.<jupyter_code># Let's try varying the number of isotonic constraints, thresholds
number_isotonic_constraints = 2
X_predict = np.linspace(10,90,number_isotonic_constraints)
iso_reg = IsotonicRegression().fit(X_train, y_train)
y_predict = iso_reg.predict(X_predict)

# Make predictions at testing data locations and calculate the model goodness
y_predict_test = np.interp(X_test, X_predict, y_predict, left=None, right=None, period=None)
r2 = explained_variance_score(y_test,y_predict_test)
print('Variance Explained = ' + str(round(r2,3)))

# Plot test predictions
plt.subplot(111)
plt.scatter(X_train, y_train,  color='black', s = 20, alpha = 0.5, label='train')
plt.scatter(X_test, y_test,  color='blue', s = 20, alpha = 0.5, label='test')
plt.scatter(X_test, y_predict_test,  color='green', s = 20, alpha = 0.5, label='predicted_test')
plt.plot(X_predict,y_predict, color='red', linewidth=1, label='model')
for threshold in X_predict:
    plt.axvline(threshold, color='red', alpha = 0.2)
plt.title('Isotonic Regression Porosity from Grainsize with Training Data'); plt.xlabel('Grain Size (mm)'); plt.ylabel('Porosity (%)')
plt.xlim(0,100)#; plt.ylim(0,1500000)
plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.0, wspace=0.2, hspace=0.2)
plt.legend()
plt.show()<jupyter_output>Variance Explained = 0.27<jupyter_text>Now let's run multiple models with variable number of isotonic constraints to perform a rough hyperparameter tuning.* simple test and train split. cross validation* model goodness by variance explained<jupyter_code>isotonic_constraint_list = np.linspace(2,35,100)
r2 = []

for k in isotonic_constraint_list:
    X_predict = np.linspace(3,97,k)
    iso_reg = IsotonicRegression().fit(X_train, y_train)
    y_predict = iso_reg.predict(X_predict)
    y_predict_test = np.interp(X_test, X_predict, y_predict, left=None, right=None, period=None)
    r2.append(explained_variance_score(y_test,y_predict_test))

plt.subplot(111)
plt.plot(isotonic_constraint_list, r2, color='red', linewidth=1, label='model')
plt.title('Isotonic Regression Tuning Number of Isotonic Constraints'); plt.xlabel('Number of Isotonic Constraints'); plt.ylabel('Variance Explained in Testing')
plt.xlim(0,40); plt.ylim(0,0.6)
plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.0, wspace=0.2, hspace=0.2)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11243016759776536, 'fraction_numerical': 0.02653631284916201, 'mean_word_length': 5.25382096069869, 'pattern_counts': {'":': 0, '<': 44, '<?xml version=': 0, '>': 44, 'https://': 9, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '188806', 'license': 'MIT', 'n_tokens_mistral': 3931, 'n_tokens_neox': 3635, 'n_words': 1116, 'path': 'SubsurfaceDataAnalytics_IsotonicRegression.ipynb', 'repo_name': 'AndrewAnnex/PythonNumericalDemos', 'size': '11463'}