<p align="center">
    

</p>

## Subsurface Data Analytics 

### Multidimensional Scaling for Subsurface Data Analytics in Python 


#### Michael Pyrcz, 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)


### Subsurface Machine Learning: Multidimensional Scaling for Subsurface Data Analytics 

Here's a demonstration of mulidimensional scaling for subsurface modeling in Python. This is part of my Subsuface Machine Learning Course at the Cockrell School of Engineering at the University of Texas at Austin.  

#### Multidimensional Scaling

A powerful ordination method in inferential statistics / information visualization for exploring / visualizing the similarity (conversely the difference) between individual samples from a high dimensional dataset.  

* beyond 2 or 3 features it is difficult to visualize the relationship between samples

* for 2 features we can easily visualize the relationships between samples with a scatter plot

* for 3 features we can either visualize in 3D or include color or matrix scatter plots

Multidimensional scaling projects the $m$ dimensional data to $p$ dimensions such that $p << m$.

* ideally we are able to project to $p=2$ to easily explore the relationships between the samples

While principal component analysis (PCA) operates with the covariance matrix, multidimensional scaling operates with the distance / dissimilarity matrix.

* you don't need to know the actual feature values, just the distance or dissimilarity between the samples

* as with any distance in feature space, we consider feature standardization and weighting

* we may also work with a variety of dissimilarity measures

#### Classical Multidimensional Scaling

Based on Euclidian distance between samples.

The Steps:

1. calculate the square distance matrix, $D^{(2)} = \left[ \delta_{i,j}^2 \right]$.

2. apply double centering $B = - \frac{1}{2} J D^{(2)} J$

3. solve for the Eigen Values, $\lambda_1,\ldots,\lambda_p$

4. solve for the projected coordinates, $x^{'}_1, \dots ,x^{'}_p$

\begin{equation}
X^{'} = E_m \land^{\frac{1}{2}}_m
\end{equation}

where $E_m$ is the matrix of Eigen vectors and $\land_m$ is the diagonal matrix of Eigen values. 

General comments about classical multidimensional scaling:

* nonlinear dimensionality reduction

* no distribution assumption

* the transform may not be unique, and may be arbitrarily be translated, rotated and transformed (these do not change the pairwise distances $||x_i - x_j||$.

#### Metric Multidimensional Scaling

A generalization of classical multidimensional scaling with a variety of metrics and a loss function optimization.

* formulated as an optimization problem to minimize the squate difference between the original and projected pairwise distances

\begin{equation}
min_{x_1,\ldots,x_m} \sum_{i<j} \left( ||x_i - x_j|| - \delta_{i,j} \right)^2
\end{equation}

where $||x_i - x_j||$ are the pairwise distances in the projected space ($p$ dimensional) and $\delta_{i,j}$ are the pairwise distances in the original feature space.

General comments about metric multidimensional scaling:

* dissimilarity measure must be meaningful

* dimensionality reduction is performed such that the error in the sample pairwise distance is minimized

* there is a variant known as Nonmetric Multidimensional Scaling for ordinal features (categorical with ordering).

#### Checking Multidimensional Scaling Results

The multidimensional scaling approach minimizes the square difference of the pairwise distances between all of the data samples and eachother between the projected, lower dimensional, and original feature space.

* **stress** is defined as:

\begin{equation}
Stress_P(x_1,\ldots,x_n) = \left( \sum_{i \ne j = 1,\ldots,n} \left( ||x_i - x_j|| - \delta_{i,j} \right)^2 \right)^{\frac{1}{2}}
\end{equation}

* it is also useful to visualize the scatterplot of projected vs. original pairwise distances

#### Comparison with Principal Component Analysis

Principal component analysis takes the covariance matrix ($m \times m$) between all the features and finds the orthogonal rotation such that the **variance is maximized** over the ordered principle components.  

Multidimensional scaling takes the matrix of the pairwise distances ($n \times n$) between all the samples in feature space and finds the nonlinear projects such that the **error in the pairwise distances is minimized**.

#### Random Projection

An alternative multidimensional scaling that relies on an random ($p \times n$) projection matrix, with columns with unit length.

* first row is a random unit vector 

* second row is orthogonal to the first row

* continue 

This random approach works practically because of the The Johnson-Lindenstrauss lemma:

*In mathematics, the Johnson-Lindenstrauss lemma is a result concerning low-distortion embeddings of points from high-dimensional into low-dimensional Euclidean space. The lemma states that a small set of points in a high-dimensional space can be embedded into a space of much lower dimension in such a way that distances between the points are nearly preserved. The map used for the embedding is at least Lipschitz, and can even be taken to be an orthogonal projection.* - Wikipedia article.

#### Objective 

In the Stochastic Machine Learning class, I want to provide hands-on experience with solving complicated subsurface modeling problems with data analytics, machine learning. 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. 

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.    

#### Getting Started

Here'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. 

You will need to copy the data file to your working directory.  They are available here:

* Tabular data - unconv_MV_v4.csv at https://git.io/fhHLT.

There are examples below with GeostatsPy functions. You can go here to see a list of the available functions, https://git.io/fh4eX, and for other example workflows and source code. 


```python
import geostatspy.GSLIB as GSLIB          # GSLIB utilies, visualization and wrapper
import geostatspy.geostats as geostats    # GSLIB methods convert to Python        
```

We will also need some standard packages. These should have been installed with Anaconda 3.


```python
import numpy as np                        # ndarrys for gridded data
import pandas as pd                       # DataFrames for tabular data
import os                                 # set working directory, run executables
import matplotlib.pyplot as plt           # for plotting
from scipy import stats                   # summary statistics
import math                               # trigonometry etc.
import scipy.signal as signal             # kernel for moving window calculation
import random                             # for randon numbers
import seaborn as sns                     # for matrix scatter plots
from scipy import linalg                  # for linear regression


from sklearn.random_projection import johnson_lindenstrauss_min_dim
from sklearn.random_projection import SparseRandomProjection
from sklearn.metrics.pairwise import euclidean_distances
```

Here's a couple of functions that I found online for metrics that we will calculate below.

#### Set the working directory

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). 


```python
os.chdir("c:/PGE383")                     # set the working directory
```

#### Loading Tabular Data

Here's the command to load our comma delimited data file in to a Pandas' DataFrame object.  


```python
df = pd.read_csv('unconv_MV_v4.csv')      # load our data table 
df['TOC'] = np.where(df['TOC']<0.0, 0.0, df['TOC']) # set TOC < 0.0 as 0.0, otherwise leave the same
```

Visualizing the DataFrame would be useful and we already learned about these methods in this demo (https://git.io/fNgRW). 

We can preview the DataFrame by utilizing the 'head' DataFrame member function (with a nice and clean format, see below). With the head command, add parameter 'n=13' to see the first 13 rows of the dataset.  


```python
df.head(n=13)                             # we could also use this command for a table preview
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Well</th>
      <th>Por</th>
      <th>Perm</th>
      <th>AI</th>
      <th>Brittle</th>
      <th>TOC</th>
      <th>VR</th>
      <th>Prod</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>1</td>
      <td>12.08</td>
      <td>2.92</td>
      <td>2.80</td>
      <td>81.40</td>
      <td>1.16</td>
      <td>2.31</td>
      <td>1695.360819</td>
    </tr>
    <tr>
      <th>1</th>
      <td>2</td>
      <td>12.38</td>
      <td>3.53</td>
      <td>3.22</td>
      <td>46.17</td>
      <td>0.89</td>
      <td>1.88</td>
      <td>3007.096063</td>
    </tr>
    <tr>
      <th>2</th>
      <td>3</td>
      <td>14.02</td>
      <td>2.59</td>
      <td>4.01</td>
      <td>72.80</td>
      <td>0.89</td>
      <td>2.72</td>
      <td>2531.938259</td>
    </tr>
    <tr>
      <th>3</th>
      <td>4</td>
      <td>17.67</td>
      <td>6.75</td>
      <td>2.63</td>
      <td>39.81</td>
      <td>1.08</td>
      <td>1.88</td>
      <td>5288.514854</td>
    </tr>
    <tr>
      <th>4</th>
      <td>5</td>
      <td>17.52</td>
      <td>4.57</td>
      <td>3.18</td>
      <td>10.94</td>
      <td>1.51</td>
      <td>1.90</td>
      <td>2859.469624</td>
    </tr>
    <tr>
      <th>5</th>
      <td>6</td>
      <td>14.53</td>
      <td>4.81</td>
      <td>2.69</td>
      <td>53.60</td>
      <td>0.94</td>
      <td>1.67</td>
      <td>4017.374438</td>
    </tr>
    <tr>
      <th>6</th>
      <td>7</td>
      <td>13.49</td>
      <td>3.60</td>
      <td>2.93</td>
      <td>63.71</td>
      <td>0.80</td>
      <td>1.85</td>
      <td>2952.812773</td>
    </tr>
    <tr>
      <th>7</th>
      <td>8</td>
      <td>11.58</td>
      <td>3.03</td>
      <td>3.25</td>
      <td>53.00</td>
      <td>0.69</td>
      <td>1.93</td>
      <td>2670.933846</td>
    </tr>
    <tr>
      <th>8</th>
      <td>9</td>
      <td>12.52</td>
      <td>2.72</td>
      <td>2.43</td>
      <td>65.77</td>
      <td>0.95</td>
      <td>1.98</td>
      <td>2474.048178</td>
    </tr>
    <tr>
      <th>9</th>
      <td>10</td>
      <td>13.25</td>
      <td>3.94</td>
      <td>3.71</td>
      <td>66.20</td>
      <td>1.14</td>
      <td>2.65</td>
      <td>2722.893266</td>
    </tr>
    <tr>
      <th>10</th>
      <td>11</td>
      <td>15.04</td>
      <td>4.39</td>
      <td>2.22</td>
      <td>61.11</td>
      <td>1.08</td>
      <td>1.77</td>
      <td>3828.247174</td>
    </tr>
    <tr>
      <th>11</th>
      <td>12</td>
      <td>16.19</td>
      <td>6.30</td>
      <td>2.29</td>
      <td>49.10</td>
      <td>1.53</td>
      <td>1.86</td>
      <td>5095.810104</td>
    </tr>
    <tr>
      <th>12</th>
      <td>13</td>
      <td>16.82</td>
      <td>5.42</td>
      <td>2.80</td>
      <td>66.65</td>
      <td>1.17</td>
      <td>1.98</td>
      <td>4091.637316</td>
    </tr>
  </tbody>
</table>
</div>



This dataset has features from 200 unconventional wells including: 

0. well index
1. well average porosity (%) 
2. permeability (mD)
3. accoustic impedance (kg/m2s*10^6)
4. brittness ratio (%) 
5. total organic carbon (%) 
6. vitrinite reflectance (%)
8. normalized initial production 90 day average (MCFPD). 

Note, the dataset is synthetic, but has realistic ranges and general multivariate relationships.

Ranking features is really an effort to understand the features and their relationships with eachother.  We will start with basic data visualization and move to more complicated methods such are partial correlation and recursive feature elimination.

#### Summary Statistics

Let's check the summary statistics of our data.


```python
df.describe().transpose()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>count</th>
      <th>mean</th>
      <th>std</th>
      <th>min</th>
      <th>25%</th>
      <th>50%</th>
      <th>75%</th>
      <th>max</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Well</th>
      <td>200.0</td>
      <td>100.500000</td>
      <td>57.879185</td>
      <td>1.000000</td>
      <td>50.750000</td>
      <td>100.500000</td>
      <td>150.250000</td>
      <td>200.000000</td>
    </tr>
    <tr>
      <th>Por</th>
      <td>200.0</td>
      <td>14.991150</td>
      <td>2.971176</td>
      <td>6.550000</td>
      <td>12.912500</td>
      <td>15.070000</td>
      <td>17.402500</td>
      <td>23.550000</td>
    </tr>
    <tr>
      <th>Perm</th>
      <td>200.0</td>
      <td>4.330750</td>
      <td>1.731014</td>
      <td>1.130000</td>
      <td>3.122500</td>
      <td>4.035000</td>
      <td>5.287500</td>
      <td>9.870000</td>
    </tr>
    <tr>
      <th>AI</th>
      <td>200.0</td>
      <td>2.968850</td>
      <td>0.566885</td>
      <td>1.280000</td>
      <td>2.547500</td>
      <td>2.955000</td>
      <td>3.345000</td>
      <td>4.630000</td>
    </tr>
    <tr>
      <th>Brittle</th>
      <td>200.0</td>
      <td>48.161950</td>
      <td>14.129455</td>
      <td>10.940000</td>
      <td>37.755000</td>
      <td>49.510000</td>
      <td>58.262500</td>
      <td>84.330000</td>
    </tr>
    <tr>
      <th>TOC</th>
      <td>200.0</td>
      <td>0.991950</td>
      <td>0.478264</td>
      <td>0.000000</td>
      <td>0.617500</td>
      <td>1.030000</td>
      <td>1.350000</td>
      <td>2.180000</td>
    </tr>
    <tr>
      <th>VR</th>
      <td>200.0</td>
      <td>1.964300</td>
      <td>0.300827</td>
      <td>0.930000</td>
      <td>1.770000</td>
      <td>1.960000</td>
      <td>2.142500</td>
      <td>2.870000</td>
    </tr>
    <tr>
      <th>Prod</th>
      <td>200.0</td>
      <td>3864.407081</td>
      <td>1553.277558</td>
      <td>839.822063</td>
      <td>2686.227611</td>
      <td>3604.303507</td>
      <td>4752.637556</td>
      <td>8590.384044</td>
    </tr>
  </tbody>
</table>
</div>



Summary statistics are a critical first step in data checking. 

* this includes the number of valid (non-null) values for each feature (count removes all np.NaN from the totals for each variable).

* we can see the general behavoirs such as central tendency, mean, and dispersion, variance. 

* we can identify issue with negative values, extreme values, and values that are outside the range of plausible values for each property. 

* We can also establish the feature ranges for plotting.  We could calculate the feature range directly from the data with code like this:

```p
Pormin = np.min(df['Por'].values)          # extract ndarray of data table column
Pormax = np.max(df['Por'].values)          # and calculate min and max
```

but, this would not result in easy to understand color bars and axis scales, let's pick convenient round numbers. We will also declare feature labels for ease of plotting. 


```python
pormin = 6.0; pormax = 24.0; porname = 'Porosity (%)'; portitle = 'Porosity' # user specified min and max values, and labels for plotting
permmin = 0.0; permmax = 10; permname = 'Permeability (mD)'; permtitle = 'Permeability'                
AImin = 1.0; AImax = 5.0; AIname = 'Acoustic Impedance (kg/m2s*10^6)'; AItitle = 'Acoustic Impedance'
brmin = 10.0; brmax = 85.0; brname = 'Brittleness Ratio (%)'; brtitle = 'Brittleness'
TOCmin = 0.0; TOCmax = 2.2; TOCname = 'Total Organic Carbon (%)'; TOCtitle = 'Total Organic Carbon' 
VRmin = 0.9; VRmax = 2.9; VRname = 'Vitrinite Reflectance (%)'; VRtitle = 'Vitrinite Reflectance'
prodmin = 500.0; prodmax = 9000.0; prodname = 'Normalized Initial Production (MCFPD)'; prodtitle = 'Normalized Initial Production'
```

The data looks to be in pretty good shape and for brevity we skip outlier detection. Let's look at the distributions with a martix scatter plot from the Seaborn package.


```python
sns.pairplot(df,vars=['Por','Perm','AI','Brittle','TOC','VR','Prod'],markers='o',plot_kws={'alpha': 0.2})
plt.subplots_adjust(left=0.0, bottom=0.0, right=0.6, top=0.6, wspace=0.3, hspace=0.2)
plt.show()
```

#### Data Preparation 

Let's make an ordinal feature from the continuous production:

1. low
2. medium
3. high
4. very high 

production rates.  This will help us visualize the results as we proceed, we can look at wells with different levels of production projected into a variety of lower dimensional spaces with multidimensional scaling.


```python
bins = [0,2500,5000,7500,10000]                # assign the production bins (these are the fence posts)
labels = ['low', 'med', 'high', 'vhigh']       # assign the labels
category = pd.cut(df['Prod'],bins,labels=labels)     # make the 1D array with the labels for our data
df['tProd'] = category                                # add the new ordinal production feature to our DataFrames    
df.head()

```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Well</th>
      <th>Por</th>
      <th>Perm</th>
      <th>AI</th>
      <th>Brittle</th>
      <th>TOC</th>
      <th>VR</th>
      <th>Prod</th>
      <th>tProd</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>1</td>
      <td>12.08</td>
      <td>2.92</td>
      <td>2.80</td>
      <td>81.40</td>
      <td>1.16</td>
      <td>2.31</td>
      <td>1695.360819</td>
      <td>low</td>
    </tr>
    <tr>
      <th>1</th>
      <td>2</td>
      <td>12.38</td>
      <td>3.53</td>
      <td>3.22</td>
      <td>46.17</td>
      <td>0.89</td>
      <td>1.88</td>
      <td>3007.096063</td>
      <td>med</td>
    </tr>
    <tr>
      <th>2</th>
      <td>3</td>
      <td>14.02</td>
      <td>2.59</td>
      <td>4.01</td>
      <td>72.80</td>
      <td>0.89</td>
      <td>2.72</td>
      <td>2531.938259</td>
      <td>med</td>
    </tr>
    <tr>
      <th>3</th>
      <td>4</td>
      <td>17.67</td>
      <td>6.75</td>
      <td>2.63</td>
      <td>39.81</td>
      <td>1.08</td>
      <td>1.88</td>
      <td>5288.514854</td>
      <td>high</td>
    </tr>
    <tr>
      <th>4</th>
      <td>5</td>
      <td>17.52</td>
      <td>4.57</td>
      <td>3.18</td>
      <td>10.94</td>
      <td>1.51</td>
      <td>1.90</td>
      <td>2859.469624</td>
      <td>med</td>
    </tr>
  </tbody>
</table>
</div>



Let's take a look at the matrix scatter plot of our 3 features and the production levels.


```python
sns.pairplot(df[['Por','AI','TOC','tProd']],markers='o',hue = 'tProd')
plt.subplots_adjust(left=0.0, bottom=0.0, right=1.5, top=1.5, wspace=0.3, hspace=0.2)
plt.show()
```

#### Multidimensional Scaling

The multidimensional scaling method follows the sample pattern as other scikit-learn methods, we instantiate, fit and then apply or transform.

Let's run multidimensional scaling on our subset of features ($m = 3$) and project to only 2 features ($p = 2$).


```python
from sklearn.manifold import MDS
n_components = 2
embedding = MDS(n_components=2)
MDS_transformed = embedding.fit_transform(df[['Por','AI','TOC']])
MDS_transformed.shape
```




    (200, 2)



The output is 2 multidimensional scaling components.  We have projected our 3 features to 2 features to minimize the error in pairwise distance between the samples.  Let's add the 2 components to our DataFrame.


```python
df['MDS1'] = MDS_transformed[:,0]
df['MDS2'] = MDS_transformed[:,1]
df.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Well</th>
      <th>Por</th>
      <th>Perm</th>
      <th>AI</th>
      <th>Brittle</th>
      <th>TOC</th>
      <th>VR</th>
      <th>Prod</th>
      <th>tProd</th>
      <th>MDS1</th>
      <th>MDS2</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>1</td>
      <td>12.08</td>
      <td>2.92</td>
      <td>2.80</td>
      <td>81.40</td>
      <td>1.16</td>
      <td>2.31</td>
      <td>1695.360819</td>
      <td>low</td>
      <td>-1.463746</td>
      <td>2.531615</td>
    </tr>
    <tr>
      <th>1</th>
      <td>2</td>
      <td>12.38</td>
      <td>3.53</td>
      <td>3.22</td>
      <td>46.17</td>
      <td>0.89</td>
      <td>1.88</td>
      <td>3007.096063</td>
      <td>med</td>
      <td>-0.877417</td>
      <td>2.474143</td>
    </tr>
    <tr>
      <th>2</th>
      <td>3</td>
      <td>14.02</td>
      <td>2.59</td>
      <td>4.01</td>
      <td>72.80</td>
      <td>0.89</td>
      <td>2.72</td>
      <td>2531.938259</td>
      <td>med</td>
      <td>0.576820</td>
      <td>1.295791</td>
    </tr>
    <tr>
      <th>3</th>
      <td>4</td>
      <td>17.67</td>
      <td>6.75</td>
      <td>2.63</td>
      <td>39.81</td>
      <td>1.08</td>
      <td>1.88</td>
      <td>5288.514854</td>
      <td>high</td>
      <td>0.820377</td>
      <td>-2.573599</td>
    </tr>
    <tr>
      <th>4</th>
      <td>5</td>
      <td>17.52</td>
      <td>4.57</td>
      <td>3.18</td>
      <td>10.94</td>
      <td>1.51</td>
      <td>1.90</td>
      <td>2859.469624</td>
      <td>med</td>
      <td>1.061472</td>
      <td>-2.338752</td>
    </tr>
  </tbody>
</table>
</div>



Let's take a look at the samples projected into the new 2 dimensional feature space.  

* note the rotation, translation is arbitrary in this space, only the sample pairwise distances are relevant. 


```python
pairplot = sns.scatterplot(x = df['MDS1'],y = df['MDS2'],hue = df['tProd'],markers='o')
plt.subplots_adjust(left=0.0, bottom=0.0, right=1., top=1.3, wspace=0.3, hspace=0.2)
plt.show()
```

Some comments:

* the general relationship between samples is preserved.

* the general production transition from low to very high production is preserved.

Let's check our model:

* we will calculate the original and projected pairwise distances between all the samples

* we will cross plot the original vs the projects pairwise distances

* we will plot the distribution of the ratio between projects / original pariwise distances



```python
dists = euclidean_distances(df[['Por','AI','TOC']], squared=False).ravel()
nonzero = dists != 0   # select only non-identical samples pairs
dists = dists[nonzero]
projected_dists = euclidean_distances(MDS_transformed, squared=False).ravel()[nonzero]

plt.subplot(221)
plt.scatter(dists,projected_dists,c='red',alpha=0.2,edgecolor = 'black')
plt.arrow(0,0,200,200,width=0.02,color='black',head_length=0.0,head_width=0.0)
plt.xlim(0,15); plt.ylim(0,15)
plt.xlabel("Pairwise Distance: original space")
plt.ylabel("Pairwise Distance: projected space")
plt.title("Pairwise Distance: Projected to %d Components" % n_components)

rates = projected_dists / dists
print("Distance Ratio, mean: %0.2f, standard deviation %0.2f." % (np.mean(rates), np.std(rates)))

plt.subplot(222)
plt.hist(rates, bins=50, range=(0.5, 1.5),color = 'red', alpha = 0.2, edgecolor='k')
plt.xlabel("Distance Ratio: projected / original")
plt.ylabel("Frequency")
plt.title("Pairwise Distance: Projected to %d Components" % n_components)

plt.subplot(223)
plt.hist(dists, bins=50, range=(0., 15.),color = 'red', alpha = 0.2, edgecolor='k')
plt.xlabel("Pairwise Distance")
plt.ylabel("Frequency")
plt.title("Pairwise Distance: Original Data")

plt.subplot(224)
plt.hist(projected_dists, bins=50, range=(0., 15.),color = 'red', alpha = 0.2, edgecolor='k')
plt.xlabel("Pairwise Distance")
plt.ylabel("Frequency")
plt.title("Pairwise Distance: Projected to %d Components" % n_components)

plt.subplots_adjust(left=0.0, bottom=0.0, right=1.7, top=2.3, wspace=0.2, hspace=0.3)
plt.show()
```

We projected to a 2 dimensional feature space and did a pretty good job preserving the pairwise distances between the samples.

#### Observing Specific Samples

Let's reduce the number of wells and actually post the well numbers and observe their spacings in the original and projected features spaces.


```python
df_subset = df.iloc[:30,:]
```

Let's look at out dataset, we will visualize the 3 possible scatter plots between the three features with the sample points labeled by well number (1 through 30).


```python
plt.subplot(131)
pairplot = sns.scatterplot(x = df_subset['Por'],y = df_subset['AI'],hue = df_subset['tProd'],markers='o')
for i, txt in enumerate(df_subset['Well']):
    pairplot.annotate(txt, (df_subset['Por'][i], df_subset['AI'][i]))
    
plt.subplot(132)
pairplot = sns.scatterplot(x = df_subset['AI'],y = df_subset['TOC'],hue = df_subset['tProd'],markers='o')
for i, txt in enumerate(df_subset['Well']):
    pairplot.annotate(txt, (df_subset['AI'][i], df_subset['TOC'][i]))
    
plt.subplot(133)
pairplot = sns.scatterplot(x = df_subset['Por'],y = df_subset['TOC'],hue = df_subset['tProd'],markers='o')
for i, txt in enumerate(df_subset['Well']):
    pairplot.annotate(txt, (df_subset['Por'][i], df_subset['TOC'][i]))
    
plt.subplots_adjust(left=0.0, bottom=0.0, right=2., top=1., wspace=0.3, hspace=0.2)
plt.show()
```

We can now see our 30 wells with their indices plots on each of the 3 possible bivariate plots.  Let's perform multidimensional scaling down to 2 components.


```python
from sklearn.manifold import MDS

embedding = MDS(n_components=2)
MDS_transformed = embedding.fit_transform(df_subset[['Por','AI','TOC']])
MDS_transformed.shape
```




    (30, 2)



Now let's visualize the wells in the projected feature space.


```python
plt.subplot(221)
pairplot = sns.scatterplot(x = df_subset['Por'],y = df_subset['AI'],hue = df_subset['tProd'],markers='o')
for i, txt in enumerate(df_subset['Well']):
    pairplot.annotate(txt, (df_subset['Por'][i], df_subset['AI'][i]))
pairplot.legend(loc='upper left')  

plt.subplot(222)
pairplot = sns.scatterplot(x = df_subset['AI'],y = df_subset['TOC'],hue = df_subset['tProd'],markers='o')
for i, txt in enumerate(df_subset['Well']):
    pairplot.annotate(txt, (df_subset['AI'][i], df_subset['TOC'][i]))
pairplot.legend(loc='lower right')     
    
plt.subplot(223)
pairplot = sns.scatterplot(x = df_subset['Por'],y = df_subset['TOC'],hue = df_subset['tProd'],markers='o')
for i, txt in enumerate(df_subset['Well']):
    pairplot.annotate(txt, (df_subset['Por'][i], df_subset['TOC'][i]))
pairplot.legend(loc='upper left')     
    
plt.subplot(224)
pairplot = sns.scatterplot(x = MDS_transformed[:,0],y = MDS_transformed[:,1],hue = df_subset['tProd'],markers='o')
for i, txt in enumerate(df_subset['Well']):
    pairplot.annotate(txt, (MDS_transformed[i,0], MDS_transformed[i,1]))
pairplot.set_xlabel('MDS1'); pairplot.set_ylabel('MDS2')  
pairplot.legend(loc='lower right') 

plt.subplots_adjust(left=0.0, bottom=0.0, right=1.7, top=1.7, wspace=0.3, hspace=0.2)
plt.show()
```

Some observations:

* the transformation is not a rotation (like PCA)

* the intersample distance is approximately preserved 

\begin{equation}
\delta_{i,j} = \left( (\delta Por_{i,j}^{2}) + (\delta AI_{i,j}^{2}) + (\delta TOC_{i,j}^{2}) \right)^{\frac{1}{2}}
\end{equation}

#### Applications

The main benefit with multiple dimensional scaling is the ability to inspect high dimensional feature spaces for relationships between samples.

* we may observe that specific wells cluster together
* we may observe systematic transitions

#### Random Projection

Let's demonstrate the use of random projection.

* given the low dimensionalilty you will see quite poor performance for some randon seeds


```python
from sklearn.random_projection import GaussianRandomProjection
n_components = 2; seed = 73058
rp = GaussianRandomProjection(n_components=n_components,random_state = seed)

dists = euclidean_distances(df[['Por','AI','TOC']], squared=False).ravel()
nonzero = dists != 0   # select only non-identical samples pairs
dists = dists[nonzero]

projected_data = rp.fit_transform(df[['Por','AI','TOC']])
projected_dists = euclidean_distances(projected_data, squared=False).ravel()[nonzero]

plt.subplot(221)
plt.scatter(dists,projected_dists,c='red',alpha=0.2,edgecolor = 'black')
plt.arrow(0,0,200,200,width=0.02,color='black',head_length=0.0,head_width=0.0)
plt.xlim(0,15); plt.ylim(0,15)
plt.xlabel("Pairwise Distance: original space")
plt.ylabel("Pairwise Distance: projected space")
plt.title("Pairwise Distance: Projected to %d Components" % n_components)

rates = projected_dists / dists
print("Distance Ratio, mean: %0.2f, standard deviation %0.2f." % (np.mean(rates), np.std(rates)))

plt.subplot(222)
plt.hist(rates, bins=50, range=(0., 2.),color = 'red', alpha = 0.2, edgecolor='k')
plt.xlabel("Distance Ratio: projected / original")
plt.ylabel("Frequency")
plt.title("Pairwise Distance: Projected to %d Components" % n_components)

plt.subplot(223)
plt.hist(dists, bins=50, range=(0., 15.),color = 'red', alpha = 0.2, edgecolor='k')
plt.xlabel("Pairwise Distance")
plt.ylabel("Frequency")
plt.title("Pairwise Distance: Original Data")

plt.subplot(224)
plt.hist(projected_dists, bins=50, range=(0., 15.),color = 'red', alpha = 0.2, edgecolor='k')
plt.xlabel("Pairwise Distance")
plt.ylabel("Frequency")
plt.title("Pairwise Distance: Projected to %d Components" % n_components)

plt.subplots_adjust(left=0.0, bottom=0.0, right=1.7, top=2.3, wspace=0.2, hspace=0.3)
plt.show()
```

We should see improved performance (lower stress) for high dimensional problems (large $m$).

#### Comments

This was a basic demonstration of multidimensional scaling. A lot more could be done:  

* comparison to principal components analysis
* use of a dataset with larger number of features
* use of a dataset with natural clusters in multivariate space

I have other demonstrations on the basics of working with DataFrames, ndarrays, univariate statistics, plotting data, declustering, data transformations, trend modeling and many other workflows available at https://github.com/GeostatsGuy/PythonNumericalDemos and https://github.com/GeostatsGuy/GeostatsPy. 
  
I hope this was helpful,

*Michael*

#### The Author:

### Michael Pyrcz, Associate Professor, University of Texas at Austin 
*Novel Data Analytics, Geostatistics and Machine Learning Subsurface Solutions*

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. 

For more about Michael check out these links:

#### [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)

#### Want to Work Together?

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.

* 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! 

* 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!

* I can be reached at mpyrcz@austin.utexas.edu.

I'm always happy to discuss,

*Michael*

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

#### 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)



```python

```

9caeb8170e78ea3427c4882adff0cc629a4b9493
ipynb
Jupyter Notebook
SubsurfaceDataAnalytics_Multidimensional_Scaling.ipynb
ushareng/PythonNumericalDemos
a78fba8234760c6449de4a2e5c66348ad9f6e61c
MIT
2020-09-29T19:33:08.000Z
2021-08-08T15:22:21.000Z
SubsurfaceDataAnalytics_Multidimensional_Scaling.ipynb
victsnet/PythonNumericalDemos
cce91eb103c3f000e4fe5a6c51ff2a3d77c6b283
MIT
SubsurfaceDataAnalytics_Multidimensional_Scaling.ipynb
victsnet/PythonNumericalDemos
cce91eb103c3f000e4fe5a6c51ff2a3d77c6b283
MIT
2020-09-26T09:27:16.000Z
2020-11-10T16:40:33.000Z
Qwen/Qwen-72B
1. YES
2. YES