<jupyter_start><jupyter_text>Convolutional Neural Network (CNN) Percolation Threshold Demonstration <NAME>, Senior Petroleum Engineering Undergraduate Student, The University of Texas at Austin [LinkedIn](https://www.linkedin.com/in/john-mccarthy2) <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) <NAME>, Graduate Candidate, The University of Texas at Austin [LinkedIn](https://www.linkedin.com/in/honggeun-jo/?originalSubdomain=kr) | [GitHub](https://github.com/whghdrms) | [Twitter](https://twitter.com/HonggeunJ) Workflow for training CNN to classify imagesThis workflow demonstrates design and training of a CNN model to classify for a variety of synthetic labeled images.* a very simple, fast run time, well-documented toy problem to support experiential learning * I demonstrate the impact of scale of the image features, the ability to work with multiple scales  ___**Can we train a simple CNN to detect the extent of a randomly located artifact in a set of uniform images**?___We design a simple network with convolution over multiple filters (variable size) to demonstrate this. Convolutional Neural NetworksExtension of the artifical neural network, based on the visual cortex:* extraction of features from overlapping receptive fields, over a hierarchy (not shown) and then recompose the whole image, our perception.* We don’t perceive all the ‘pixels’, our visual cortex interprets and summarizes patterns. Let’s make a machine to do this.**Regularization**: a constraint to reduce the sensitivity of the model to the data, to reduce model varianceReceptive Fields* the use of receptive fields is a form of regularization* massive reduction in connections, weights and model parameters* effectively shrinking these potential weights to zero* while integrating / focusing on pixel patterns!We have access to operators to move from layer to layer (feature maps to feature maps) in our convolutional neural networks.  The common operators include:* **Convolution** – a weighting window / filter designed to extract features* **Pooling** – reduction in dimensionality, increase local translation invariance * **Depth-wise Pooling, Down Sampling** – 1x1 filter that combine channels, feature maps* **Activation** – use of an activation function to apply a nonlinear transformation to impart nonlinearity to the system* **Full-connected, feed forward** – see previous lecture For a demonstration of all of these operators, check out this [Python Jupyter Notebook](https://github.com/GeostatsGuy/PythonNumericalDemos/blob/master/SubsurfaceDataAnalytics_Convolution_Operators.ipynb) Objective Michael teaches data analytics, geostatistics and machine learning. To demonstrate the basic construction of a convolutional neural networks, training and prediction.  * Michael uses these examples in his lecture notes, see the lecture posted on his [YouTube channel](https://www.youtube.com/watch?v=za2my_XDoOs&ab_channel=GeostatsGuyLectures).* for any to gain experiential learning with the nuts and bolts of convolutional neural networksNote, we only demonstrate the basics of construction, training and prediction. There is much more about the design and training that we do not cover. 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 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.  Load the required librariesThe following code loads the required libraries.We will need some standard packages. These should have been installed with Anaconda 3.* [numpy](https://numpy.org/): To generate arrays * [matplotlib](https://matplotlib.org/): Vilsualization purpose * sklearn: for model metrics, confusion matrices, one hot encoderWe also will need tensor flow, this will require an install as it is not available in Anaconda* [tensorflow > 2.0.0](https://www.tensorflow.org/learn): Design, compile and train neural network models<jupyter_code>import numpy as np                                      # ndarrays for gridded data, our images
import matplotlib.pyplot as plt                         # plotting                   
from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) # control of axes ticks
import pandas as pd
from sklearn.metrics import f1_score, accuracy_score, recall_score, confusion_matrix # model metrics
import scipy
from scipy.ndimage import gaussian_filter               # Gaussian filter for smoothing our images
from sklearn.preprocessing import OneHotEncoder         # one hot encoder for our response feature
import geostatspy.GSLIB as GSLIB                        # GSLIB utilies, visualization and wrapper
import geostatspy.geostats as geostats                  # GSLIB methods convert to Python 
import random                                           # pseudo-random values

seed = 73073                                            # set the random number seed for repeatability
np.random.seed(seed)                                    

cmap = plt.cm.inferno                                   # color map for plots

import warnings
warnings.filterwarnings("ignore")
import tensorflow as tf                                 # import tensor flow
from tensorflow import keras
from tensorflow.keras import layers
from tensorflow.keras.models import Model
import visualkeras                                      # visualize CNN architectures

assert tf.__version__.startswith('2.')                  # this will give you error if tensorflow < 2.0
tf.keras.backend.set_floatx('float32')                  # default float to be 'float32' (Tensorflow only works with float8, 16, 32)

physical_devices = tf.config.list_physical_devices('GPU') # constrain memory use to avoid a CUDNN_STATUS_ALLOC_FAILED error
#tf.config.experimental.set_memory_growth(physical_devices[0], True)<jupyter_output><empty_output><jupyter_text>Define FunctionsLet's define a couple of convenience functions for workflow brevity and read-ability* train and test accuracy over epochs performance plots* train and test prediction acurracy, confusion matrix plots* visualize the trained CNN filters* visualize the feature maps for a specific image classification case* visualize the last layer, the probability output from the softmax activation function  Code below taken from the Simple Convolutional Neural Network Scale Demonstration [Workflow](https://github.com/GeostatsGuy/PythonNumericalDemos/blob/master/SubsurfaceDataAnalytics_ConvolutionalNeuralNetworks_Classifier_Scale.ipynb)<jupyter_code>def model_performance(history,nepoch):                  # visualize error and loss, training and testing over Epochs     
    plt.subplot(1,2,1)
    plt.plot(history.history['loss'],c='red')
    plt.plot(history.history['val_loss'],c='black')
    plt.title('CNN Model Loss')
    plt.ylabel('Testing Categorical Entropy'); plt.xlabel('Epoch'); plt.grid()
    plt.legend(['train', 'test'], loc='upper right'); plt.xlim([0,nepoch]); plt.ylim([0,1])
    
    plt.subplot(1,2,2)
    plt.plot(history.history['sparse_categorical_accuracy'],c='red')
    plt.plot(history.history['val_sparse_categorical_accuracy'],c='black')
    plt.title('CNN Model Accuracy')
    plt.ylabel('Testing Proportion Correct Classification'); plt.xlabel('Epoch'); plt.grid()
    plt.legend(['train', 'test'], loc='upper right'); plt.xlim([0,nepoch]); plt.ylim([0,1])
    plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.1, wspace=0.3, hspace=0.3)
    plt.show()
    
def model_cross_validation(model,X_train,X_test,y_train,y_test): # visualize misclassification error for training and testing data 
    y_train_predict = model.predict(X_train[:,:,:,0].reshape([ntrain,nx,ny,-1])) # predict over the training images
    y_train_predict = np.argmax(y_train_predict,axis = 1) # assign the maximum probability category
    
    y_test_predict = model.predict(X_test[:,:,:,0].reshape([ntest,nx,ny,-1])) # predict over the testing images 
    y_test_predict = np.argmax(y_test_predict,axis = 1) # assign the maximum probability category
    
    plt.subplot(121)                                                
    train_confusion_matrix = confusion_matrix(y_train,y_train_predict) # calculate and plot train confusion matrix
    train_confusion_matrix = train_confusion_matrix/np.sum(train_confusion_matrix) 
    im = plt.imshow(train_confusion_matrix,cmap=cmap,vmin=0.0,vmax=1.0); plt.xlabel('Prediction'); plt.ylabel('Truth'); plt.title('Training Confusion Matrix')
    ax = plt.gca(); ax.set_xticks([0,1]); ax.set_xticklabels(name)
    ax.set_yticks([0,1]); ax.set_yticklabels(name)
    cbar = plt.colorbar(im, orientation="vertical", ticks=np.linspace(0, 1.0, 5))
    cbar.set_label('Proportion', rotation=270, labelpad=20)
    
    plt.subplot(122)
    test_confusion_matrix = confusion_matrix(y_test,y_test_predict) # calculate and plot train confusion matrix
    test_confusion_matrix = test_confusion_matrix/np.sum(test_confusion_matrix) 
    im = plt.imshow(test_confusion_matrix,cmap=cmap,vmin = 0,vmax = 1.0)
    plt.xlabel('Prediction'); plt.ylabel('Truth'); plt.title('Testing Confusion Matrix')
    ax = plt.gca(); ax.set_xticks([0,1]); ax.set_xticklabels(name)
    ax.set_yticks([0,1]); ax.set_yticklabels(name)
    cbar = plt.colorbar(im, orientation="vertical", ticks=np.linspace(0, 1.0, 5))
    cbar.set_label('Proportion', rotation=270, labelpad=20)
    
    plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.1, wspace=0.3, hspace=0.3)
    plt.show()
    
    print(f'            Train   Test')                  # print a table with train and test accuracy
    print(f'precision: {accuracy_score(y_train.flatten(), y_train_predict):.4f}, {accuracy_score(y_test.flatten(), y_test_predict):.4f}')
    print(f'recall:    {recall_score(y_train.flatten(), y_train_predict, average = "weighted"):.4f}, {recall_score(y_test.flatten(), y_test_predict, average = "weighted"):.4f}')
    print(f'f1 score:  {f1_score(y_train.flatten(), y_train_predict, average = "weighted"):.4f}, {f1_score(y_test.flatten(), y_test_predict, average = "weighted"):.4f}')

def visualize_filters(model):                           # visualize CNN filters
    for layer in model.layers:
        if 'conv' not in layer.name:                    # check if layers is a convolution
            continue
        filters, biases = layer.get_weights()
        fmin, fmax = filters.min(), filters.max()       # normalize the filter weights for visualization
        filters = (filters - fmin) / (fmax - fmin)
        print('Convolution Layer Name: ' + layer.name + ', Shape: ' + str(filters.shape[0]) + ' x ' + str(filters.shape[1]) + ', Number of Channels: ' + str(filters.shape[3]) + '.')
        nch = filters.shape[3]
        for ich in range(0,nch):
            plt.subplot(1,nch,ich+1)
            im = plt.imshow(filters[:,:,0,ich],cmap=cmap,vmin=0.0,vmax=1.0)
            ax = plt.gca(); ax.set_xticks([0,1,2]); ax.set_xticklabels([-1,0,1])
            ax.set_yticks([0,1,2]); ax.set_yticklabels([-1,0,1])
            cbar = plt.colorbar(im, orientation="vertical", ticks=np.linspace(0.0, 1.0, 5))
            cbar.set_label('Weights', rotation=270, labelpad=20)
            plt.title('Filter ' + str(ich+1))
        plt.subplots_adjust(left=0.0, bottom=0.0, right=0.5*nch, top=0.5, wspace=0.4, hspace=0.3)
        plt.show()

def visualize_convolution_feature_maps(model,image):    # visualize CNN feature maps for a prediction case
    ilayer = -1
    for layer in model.layers:
        ilayer = ilayer + 1
        # check for convolutional layer
        if 'conv' not in layer.name:
            continue
        model_trunc = Model(inputs=model.inputs, outputs=model.layers[ilayer].output) # truncate model to output feature maps
        feature_maps = model_trunc.predict(X[image,:,:,0].reshape([1,nx,ny,-1]))        
        print('Convolution Layer Name: ' + layer.name + ', Shape: ' + str(feature_maps.shape[1]) + ' x ' + str(feature_maps.shape[2]) + ', Number of Channels: ' + str(feature_maps.shape[3]) + '.')
    
        nch = feature_maps.shape[3]
        plt.subplot(1,nch+1,1)
        im = fig = plt.imshow(X[image,:,:,0], cmap=cmap)
        fig.axes.get_xaxis().set_visible(False); fig.axes.get_yaxis().set_visible(False)
        plt.title('Image ' + str(image) + ': ' + names[y[image]])
    
        for ich in range(0,nch):
            plt.subplot(1,nch+1,ich+2)
            im = plt.imshow(feature_maps[0, :, :, ich], vmin = -1,vmax = 1.0,cmap=cmap)
            ax = plt.gca(); ax.set_xticks([]); ax.set_yticks([])
            plt.title('Convolution Feature Map ' + str(ich+1))
            cbar = plt.colorbar(im, orientation="vertical", ticks=np.linspace(-1.0, 1.0, 5))
            cbar.set_label('Feature Map Values', rotation=270, labelpad=20)
        
        plt.subplots_adjust(left=0.0, bottom=0.0, right=0.5*nch, top=0.35, wspace=0.4, hspace=0.3)
        plt.show()
        
def visualize_max_pooling_feature_maps(model,image):    # visualize CNN feature maps for a prediction case
    ilayer = -1
    for layer in model.layers:
        ilayer = ilayer + 1
        # check for convolutional layer
        if 'max_pooling' not in layer.name:
            continue
        model_trunc = Model(inputs=model.inputs, outputs=model.layers[ilayer].output) # truncate model to output feature maps
        feature_maps = model_trunc.predict(X[image,:,:,0].reshape([1,nx,ny,-1]))        
        print('Max Pooling Layer Name: ' + layer.name + ', Shape: ' + str(feature_maps.shape[1]) + ' x ' + str(feature_maps.shape[2]) + ', Number of Channels: ' + str(feature_maps.shape[3]) + '.')
    
        nch = feature_maps.shape[3]
        plt.subplot(1,nch+1,1)
        fig = plt.imshow(X[image,:,:,0], cmap=cmap)
        fig.axes.get_xaxis().set_visible(False); fig.axes.get_yaxis().set_visible(False)
        plt.title('Image ' + str(image) + ': ' + names[y[image]])
    
        for ich in range(0,nch):
            plt.subplot(1,nch+1,ich+2)
            im = plt.imshow(feature_maps[0, :, :, ich], cmap=cmap)
            ax = plt.gca(); ax.set_xticks([]); ax.set_yticks([])
            plt.title('Max Pooling Feature Map ' + str(ich+1))
        
        plt.subplots_adjust(left=0.0, bottom=0.0, right=0.5*nch, top=0.5, wspace=0.4, hspace=0.3)
        plt.show()

def visualize_last_layers(model, image):                # visualize last layer, category probabilities
    ilayer = -1; llayer = -1
    for layer in model.layers:                          # get the last dense layer
        ilayer = ilayer + 1
        # check for convolutional layer
        if 'dense' not in layer.name:
            continue
        else:
            llayer = ilayer
    if llayer == -1:
        print('No dense layer found'); return -1
    model_trunc = Model(inputs=model.inputs, outputs=model.layers[ilayer].output) # truncate model to output feature maps
    output_layer = model_trunc.predict(X[image,:,:,0].reshape([1,ny,nx,-1]))        
    output_layer = output_layer.flatten()*100
    
    
    plt.subplot(1,2,1)
    fig = plt.imshow(X[image,:,:,0], cmap=cmap)
    fig.axes.get_xaxis().set_visible(False); fig.axes.get_yaxis().set_visible(False)
    plt.title('Image ' + str(image) + ': ' + names[y[image]])

    plt.subplot(1,2,2)
    plt.title('Predicted Category Probabilities')
    plt.bar(x = ['Connected', 'Not Connected'], height = output_layer,color='red',edgecolor='black',alpha=0.2)
    plt.xlabel('Image Categories'); plt.ylabel('Probability (\%)')
    plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.1, wspace=0.4, hspace=0.3)
    plt.show()              
    
def visualize_testing_predictions(model):               # visualize testing images and predicted probabilities 
    y_test_predict = model.predict(X_test[:,:,:,0].reshape([ntest,nx,ny,-1])) # predict over the testing images 
    y_test_predict_cat = np.argmax(y_test_predict,axis = 1) # assign the maximum probability category
    
    c, r = 4, 5                                         
    plt.figure(figsize = (15,15))
    for i in range(c*r):
        plt.subplot(r,c,i+1)
        fig = plt.imshow(X_test[i,:,:,0],cmap=cmap)
        fig.axes.get_xaxis().set_visible(False); fig.axes.get_yaxis().set_visible(False)
        icat = y_test_predict_cat[i]
        plt.title('Image ' + str(i) + ': ' + names[y_test[i]] + ', ' + str(np.round(y_test_predict[i,icat]*100,2)) + "% " + names[icat] )

            
    plt.subplots_adjust(left=0.0, bottom=0.0, right=1.3, top=1.6, wspace=0.01, hspace=0.3)
    plt.show()
    
def prediction_boxplot(model):                             # visualize the box plots for predicted probabilities by true category
    y_train_predict = model.predict(X_train[:,:,:,0].reshape([ntrain,ny,nx,-1])) # predict over the testing images 
    y_test_predict = model.predict(X_test[:,:,:,0].reshape([ntest,ny,nx,-1])) # predict over the testing images 
    c = 'red'
    
    plt.subplot(121)
    box_train = plt.boxplot([y_train_predict[:,0][y_train[:] == 0],y_train_predict[:,1][y_train[:] == 1]],labels=['Connected','Not Connected'])
    for i in range(1,4): plt.plot([i+.5,i+.5],[0,1],color='k',alpha=0.3,linestyle='--',linewidth=1)
    plt.gca().yaxis.grid(True, which='major',linewidth = 1.0); plt.gca().yaxis.grid(True, which='minor',linewidth = 0.2) # add y grids
    plt.gca().tick_params(which='major',length=7); plt.gca().tick_params(which='minor', length=4)
    plt.gca().yaxis.set_minor_locator(AutoMinorLocator())   
    plt.ylim([0,1]); plt.xlim([.5,2.5]); plt.ylabel('Prediction Probability of Correct Category'); plt.xlabel('Correct Image Label')
    plt.title('Training Image Correct Category Probabilities')

    
    plt.subplot(122)
    box_train = plt.boxplot([y_test_predict[:,0][y_test[:] == 0],y_test_predict[:,1][y_test[:] == 1]],labels=['Connected','Not Connected'])
    for i in range(1,4): plt.plot([i+.5,i+.5],[0,1],color='k',alpha=0.3,linestyle='--',linewidth=1)
    plt.gca().yaxis.grid(True, which='major',linewidth = 1.0); plt.gca().yaxis.grid(True, which='minor',linewidth = 0.2) # add y grids
    plt.gca().tick_params(which='major',length=7); plt.gca().tick_params(which='minor', length=4)
    plt.gca().yaxis.set_minor_locator(AutoMinorLocator())   
    plt.ylim([0,1]); plt.xlim([.5,2.5]); plt.ylabel('Prediction Probability of Correct Category'); plt.xlabel('Correct Image Label')
    plt.title('Testing Image Correct Category Probabilities')
    
    plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.0, wspace=0.3, hspace=0.3)
    plt.show()<jupyter_output><empty_output><jupyter_text>Depth-first search (DFS) algorithmode used to identify images exceeding the percolation threshold:<jupyter_code>def maxAreaOfIsland(grid):
        max_area = 0
        for r in range(len(grid)):
            for c in range(len(grid[r])):
                if grid[r][c] == 1:
                    max_area = max(max_area, explore(grid, r, c))
        return max_area
    
def explore(grid, r, c):
    if (r < 0 or r >= len(grid)) or (c < 0 or c >= len(grid[r])):
        return 0
    if grid[r][c] == 0:
        return 0
    # mark as visited
    grid[r][c] = 0
    return 1 + explore(grid, r - 1, c) + explore(grid, r, c - 1) \
        + explore(grid, r + 1, c) + explore(grid, r, c + 1)<jupyter_output><empty_output><jupyter_text>Convolutional Neural Network ExampleLet's try a more complicated problem that will require our model to be location invariant and multiscale. * we will create synthetic images that imitate seismic tomography* we will identify geobodies exceeding the **percolation threshold** using a CNN* the percolation threshold is a mathematical concept in percolation theory that describes the formation of long-range connectivity in random systems* The percolation threshold of a given system is dependent upon the shape of its lattice. We will be generating 2d square lattices that have a threshold of 1/2 Create ImagesWe make 100 simple images, 80 train and 20 test 12x12 binary, 0 and 1 images.* we do this for a very rapid, fast toy problemWe make the random training and testing images and visualize some of the testing models with their labels.<jupyter_code>np.random.seed(seed)                                    # make synthetic images
random.seed(seed)
nmodel = 100; ptest = 0.2; nx = 12; ny = 12
ntest = int(nmodel*ptest); ntrain = nmodel - ntest
X = np.random.rand(nmodel,nx,ny,1)
y = np.zeros(nmodel,dtype = int); name = ['anistropic goebody ']
X_train = np.zeros([ntrain,nx,ny,1]); y_train = np.zeros(ntrain,dtype = int)
X_test = np.zeros([ntest,nx,ny,1]); y_test = np.zeros(ntest,dtype = int)
filter_kernel = [[1/144, 1/144, 1/144, 1/144, 1/144, 1/144, 1/144, 1/144]]

    
for i in range(0,nmodel):
    image = X[i,:,:,0]
    X[i,:,:,0] = scipy.signal.convolve2d(image, filter_kernel, mode='same', boundary='fill', fillvalue=0)

for i in range(0,nmodel):
    data = np.concatenate(X[i,:,:,0], axis=None)
    mean = sum(data)/ len(data)
    target_std = .03
    target_value = .15
    std = np.std(data)
    X[i,:,:,0] = (target_std / std) * (X[i,:,:,0] - mean) + target_value
    
c, r = 7, 3                                             # visualize labelled synthetic images
plt.figure(figsize = (15,15))
for i in range(c*r):
    plt.subplot(r,c,i+1)
    fig = plt.imshow(X[i,:,:,0],cmap=cmap)
    fig.axes.get_xaxis().set_visible(False); fig.axes.get_yaxis().set_visible(False)
    plt.title(name[y[i]] + str(i))

plt.subplots_adjust(left=0.0, bottom=0.0, right=1.3, top=0.5, wspace=.5, hspace=0.3)
plt.show()
hist_data2 = np.concatenate(X[1,:,:,0], axis=None)
GSLIB.hist_st(hist_data2,np.min(hist_data2),np.max(hist_data2),log=False,cumul = False,bins=10,weights = None,xlabel='porosity (%)',title='porosity')<jupyter_output><empty_output><jupyter_text>Number of Porous and Non-Porous Cells in Samples:<jupyter_code>plt.figure(figsize = (15,15))
for i in range(c*r):
    plt.subplot(r,c,i+1)
    non_porous = [x for x in np.concatenate(X[i,:,:,0], axis=None) if x <= .15]
    porous = [x for x in np.concatenate(X[i,:,:,0], axis=None) if x >= .15]
    plt.hist(non_porous, bins=5,label='non_porous',alpha=.5)
    plt.hist(porous, bins=5,label='porous',alpha=.5)
    plt.title(name[y[i]] + str(i))
    plt.legend()
plt.subplots_adjust(left=0.0, bottom=0.0, right=1.3, top=0.5, wspace=.5, hspace=0.3)
plt.show()<jupyter_output><empty_output><jupyter_text>Number of Porous and Non-Porous Cells in All Samples:<jupyter_code>plots = []
for i in range(0,nmodel):
    plots.append(X[i,:,:,0])
non_porous = [x for x in np.concatenate(plots, axis=None) if x <= .15]
porous = [x for x in np.concatenate(plots, axis=None) if x >= .15]
plt.hist(non_porous, bins=5,label='Non Porous',alpha=.5)
plt.hist(porous, bins=5,label='Porous',alpha=.5)
plt.title("Number of Porous and Non-Porous Cells in All Samples")
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Proportion of Sand in All Samples:<jupyter_code>z = np.zeros(nmodel,dtype = int)
for i in range(0,nmodel):
    occurrences_more_than_15 = X[i,:,:,0] > .15
    z[i] = occurrences_more_than_15.sum() 
z = z /144 * 100
plt.hist(z[:], bins=5,alpha=1)
plt.title("Proportion of Sand in All Samples")
plt.xlabel('Proportion of Sand in Samples (%)')
plt.ylabel('Frequency')
plt.show()<jupyter_output><empty_output><jupyter_text>Create image labels:<jupyter_code>names = ['connected', 'not connected']
for i in range(0,nmodel):
    array = X[i,:,:,0]
    new_array = np.where(array < .15, 0, 1)
    new_array = new_array.tolist()
    if maxAreaOfIsland(new_array) > 72:
        y[i] = 0
    else:
        y[i] = 1
not_connected_occurrences = np.count_nonzero(y == 1)
print(str(not_connected_occurrences) + ' images are not connected')
connected_occurrences = np.count_nonzero(y == 0)
print(str(connected_occurrences) + ' images are connected')
c, r = 7, 3                                             # visualize labelled synthetic images
plt.figure(figsize = (15,15))
for i in range(c*r):
    plt.subplot(r,c,i+1)
    fig = plt.imshow(X[i,:,:,0],cmap=cmap)
    fig.axes.get_xaxis().set_visible(False); fig.axes.get_yaxis().set_visible(False)
    plt.title('Image '+ str(i) + str(': ') + names[y[i]])

plt.subplots_adjust(left=0.0, bottom=0.0, right=1.3, top=0.5, wspace=.5, hspace=0.3)
plt.show()<jupyter_output><empty_output><jupyter_text>Create binary images to increase speed of CNN model* the yellow is representative of porous media exceeding 15% porosity * the purple represents porous media below 15% porostiy* this is done so that the ideal porous media for connectivity is highlighted<jupyter_code>c, r = 7, 3                                             # visualize labelled synthetic images
plt.figure(figsize = (15,15))
for i in range(0,nmodel):
    array = X[i,:,:,0]
    X[i,:,:,0] = np.where(array < .15, 0, 1)
    for i in range(c*r):
        plt.subplot(r,c,i+1)
        fig = plt.imshow(X[i,:,:,0],cmap=cmap)
        fig.axes.get_xaxis().set_visible(False); fig.axes.get_yaxis().set_visible(False)
        plt.title('Image '+ str(i) + str(': ') + names[y[i]])

plt.subplots_adjust(left=0.0, bottom=0.0, right=1.3, top=0.5, wspace=.5, hspace=0.3)
plt.show()
for i in range(0,nmodel):
    if i < ntest:
        X_test[i,:,:,0] = X[i,:,:,0]; 
        y_test[i] = y[i]
    else:
        X_train[i-ntest,:,:,0] = X[i,:,:,0]; 
        y_train[i-ntest] = y[i]<jupyter_output><empty_output><jupyter_text>Step 1. Define Classification Convolutional Neural Network Now we will define our convolutional neural network* with keras frontend for tensor flow it is not too dificult to design our network* the overall architure looks like this:![image info](https://miro.medium.com/max/2000/1*vkQ0hXDaQv57sALXAJquxA.jpeg)Image taken from [blog post](https://towardsdatascience.com/a-comprehensive-guide-to-convolutional-neural-networks-the-eli5-way-3bd2b1164a53)) by <NAME>.We cycle multiple layers of:* convolution and activation with ReLU* max poolingAfter multiple cycles we have learning the features in the images we then finish with:* flattening the feature into a 1D vector* feed-forward, fully-connected artificial neural network with 2 output nodes for the probability of each categoryAs specified below the model includes:* cycles of convolution with 3x3 kernels, stride = 1 for feature maps extent reduction (one half)* feature maps transition from 12x12x1 [nx,ny,nchannel] over multiple steps to reduce featured map extent while increasing the depth and then flattens to a vector fully connected with 2 output nodes We begin by tuning the hyperparameters of multiple models at once to see which model performs the best:<jupyter_code>def CNN_model(Num_channel):                                        # CNN design function
    model = tf.keras.Sequential()                       # define neural network model sequentially
    # Convolution Feature Map 1: (12x12x1) --> (12x12x6)
    model.add(layers.Conv2D(int(Num_channel), kernel_size=(3,3),strides=1,padding='same',input_shape=[ny,nx,1])) # -2 feature map extent given filter and no padding
    model.add(layers.ReLU())
    # Max Pooling Feature Map 1: (12x12x6) --> (6x6x6)    
    model.add(layers.MaxPooling2D((2, 2)))              # 1/2 extent
    # Convolution Feature Map 2: (6x6x6) --> (6x6x12)
    model.add(layers.Conv2D(int(Num_channel*2), kernel_size=(3,3),strides=1,padding='same',input_shape=[ny,nx,1])) # -2 feature map extent given filter and no padding
    model.add(layers.ReLU())
    # Max Pooling Feature Map 2: (6x6x12) --> (3x3x12)    
    model.add(layers.MaxPooling2D((2, 2)))              # 1/2 extent
    # Convolution Feature Map 3: (3x3x12) --> (3x3x24)
    model.add(layers.Conv2D(int(Num_channel*4), kernel_size=(3,3),strides=1,padding='same',input_shape=[ny,nx,1])) # -2 feature map extent given filter and no padding
    model.add(layers.ReLU())
    # Max Pooling Feature Map 3: (3x3x24) --> (1x1x12)    
    model.add(layers.MaxPooling2D((3, 3)))              # 1/2 extent
    model.add(layers.Flatten())
    # Output layer: (12) --> 2 (i.e., each node corresponds to the probability to be each class)
    model.add(layers.Dense(2, activation = 'softmax'))  # softmax activation function for classfier probabilities
    # Compile the Neural Network - define Loss and optimizer to tune the associated weights
    loss = tf.keras.losses.SparseCategoricalCrossentropy()
    metrics = tf.keras.metrics.SparseCategoricalAccuracy()
    model.compile(loss=loss, metrics=[metrics], optimizer='adam')
    return model

nepoch = 100; batch_size = 5                             # training parameters


for i in [4,5,6,7]:
    tf.random.set_seed(seed)
    model = CNN_model(i)
    history = model.fit(X_train[:,:,:,0].reshape([ntrain,ny,nx,-1]), y_train,
                          batch_size=batch_size, epochs=nepoch, verbose=0,
                          validation_data=(X_test[:,:,:,0].reshape([ntest,ny,nx,-1]), y_test))
    print(f'\n\n Result when {i} initial channels are used')
    model_performance(history,nepoch)                       # plot loss and accuracy over training Epochs<jupyter_output>Result when 4 initial channels are used<jupyter_text>From the results we can conclude that the models starting with 4 and 6 channels performed the best. We will go forth with the model starting with 6 channels* Starting with 6 channels allows for a diversity of filters that aid in identification of artifacts* Using (3,3) kernel filters reduces the amount of information coming from padding * When using larger kernel filters, a large amount of information would come from "padding" (which are simply zeros)* We will use feature maps that have (12x12), (6x6), and (3x3) sizes * Max pooling is used after each convolution to help find the feature map containing the most prominent features of the previous feature map* Optimize using sparse categorical crossentropy to compute the crossentropy loss between the labels and predictions* ReLu and Softmax are used to restrict values between [0,1], make probability predictions, and prevent the “vanishing gradient” problem caused by other activation functions* Epochs are reduced from 100 to prevent over training<jupyter_code>def CNN_model():                                        # CNN design function
    model = tf.keras.Sequential()                       # define neural network model sequentially
    # Convolution Feature Map 1: (12x12x1) --> (12x12x6)
    model.add(layers.Conv2D(6, kernel_size=(3,3),strides=1,padding='same',input_shape=[ny,nx,1])) # -2 feature map extent given filter and no padding
    model.add(layers.ReLU())
    # Max Pooling Feature Map 1: (12x12x6) --> (6x6x6)    
    model.add(layers.MaxPooling2D((2, 2)))              # 1/2 extent
    # Convolution Feature Map 2: (6x6x6) --> (6x6x12)
    model.add(layers.Conv2D(12, kernel_size=(3,3),strides=1,padding='same',input_shape=[ny,nx,1])) # -2 feature map extent given filter and no padding
    model.add(layers.ReLU())
    # Max Pooling Feature Map 2: (6x6x12) --> (3x3x12)    
    model.add(layers.MaxPooling2D((2, 2)))              # 1/2 extent
    # Convolution Feature Map 3: (3x3x12) --> (3x3x24)
    model.add(layers.Conv2D(24, kernel_size=(3,3),strides=1,padding='same',input_shape=[ny,nx,1])) # -2 feature map extent given filter and no padding
    model.add(layers.ReLU())
    # Max Pooling Feature Map 3: (3x3x24) --> (1x1x12)    
    model.add(layers.MaxPooling2D((3, 3)))              # 1/2 extent
    model.add(layers.Flatten()) 
    # Output layer: (12) --> 2 (i.e., each node corresponds to the probability to be each class)
    model.add(layers.Dense(2, activation = 'softmax'))  # softmax activation function for classfier probabilities
    # Compile the Neural Network - define Loss and optimizer to tune the associated weights
    loss = tf.keras.losses.SparseCategoricalCrossentropy()
    metrics = tf.keras.metrics.SparseCategoricalAccuracy()
    model.compile(loss=loss, metrics=[metrics], optimizer='adam')
    return model<jupyter_output><empty_output><jupyter_text>Instantiate, Train the Convolutional Neural Network and Visualize the Model Performance in Training and TestingThis includes the following steps:1. instantiate the CNN specified above 2. train it with the 80 images in the training set3. visualize the training and testing accuracy over the Epochs of training4. write out the summary of the convolutional neural network Below is the actual CNN model that we just created. If you go back and adjust the model, that change will be represented below<jupyter_code>model = CNN_model() 
visualkeras.layered_view(model, legend=True, scale_xy=20, scale_z=5, spacing=30)
nepoch = 94; batch_size = 5                             # training parameters

tf.random.set_seed(seed) 
model = CNN_model() 
history = model.fit(X_train[:,:,:,0].reshape([ntrain,ny,nx,-1]), y_train, 
                      batch_size=batch_size, epochs=nepoch, verbose=0, 
                      validation_data=(X_test[:,:,:,0].reshape([ntest,ny,nx,-1]), y_test))
    
model_performance(history,nepoch)                       # plot loss and accuracy over training Epochs
model.summary()                                         # write out the model summary<jupyter_output><empty_output><jupyter_text>Predict with the Trained Convolutional Neural NetworkNow we load the trained deeper CNN classifier and visualize its prediction performance* we predict over the training and testing image datasets* we predict with the category assigned the maximum probability by our model<jupyter_code>model_cross_validation(model,X_train,X_test,y_train,y_test) # predict and cross validate CNN classifier<jupyter_output><empty_output><jupyter_text>Interogate the Trained Convolutional Neural Network Visualize the Convolution FiltersWe can access and visualize the convolution weights. The follow code:* loops over the layers and finds the convolution layers* loops over the channels and plots all of the filters<jupyter_code>visualize_filters(model)                                # visualize all of the CNN filter weights<jupyter_output>Convolution Layer Name: conv2d_18, Shape: 3 x 3, Number of Channels: 6.<jupyter_text>These are small 3 x 3 filters, but we can recognize some common structures:* averaging / smoothing* gradient* edge / sharpen (sobel)and combinations of these structures.* it makes sense that averaging, edge and gradient filters would differentiate the random and smooth images. Visualize the Convolution Feature MapsLet's visualize the feature maps after each convolution for the classification prediction problem * you can specify the image below from 0 - 99<jupyter_code>image = 55                                              # selected prediction case image from 0 - 99
visualize_convolution_feature_maps(model, image)        # function to visualize the convolution feature maps<jupyter_output>Convolution Layer Name: conv2d_18, Shape: 12 x 12, Number of Channels: 6.<jupyter_text>These feature maps in each layer demonstrate the extraction of information from the original images:* note the identification of the artifact over each convolution filter. Visualize the Max Pooling Feature MapsLet's visualize the feature maps after max pooling operations for our classification prediction problem.* you can specify the image below from 0 - 99<jupyter_code>image = 55                                              # selected prediction case image from 0 - 99
visualize_max_pooling_feature_maps(model, image)        # function to visualize the max pooling feature maps<jupyter_output>Max Pooling Layer Name: max_pooling2d_18, Shape: 6 x 6, Number of Channels: 6.<jupyter_text>Visualize the Output LayerLet's visualize the output layer which contains the probability to be each categories (e.g., random or smooth)* you can specify the image below from 0 - 99<jupyter_code>image = 55                                              # selected prediction case image from 0 - 99
visualize_last_layers(model, image)                     # function to visualize the output layer<jupyter_output><empty_output><jupyter_text>Other Convolutional Neural Network Architecture ExamplesSome other examples of layers that could be added to the above CNN1. More dense layers, ANN feed forward```python    model.add(layers.Dense(512))    model.add(layers.BatchNormalization(momentum=0.8))    model.add(layers.ReLU())    model.add(layers.Dropout(0.25))```2. Max pooling 2D layers```python    model.add(layers.MaxPooling2D((2, 2)))``` Visualize the Testing Images and Predicted Classification ProbabilitiesLet's look at all the testing images and the predicted classification probabilities from our our trained CNN.* we report the true label, and then the greatest category probability, > 50% (argmax), and the associated category from our CNN * for example, 'connected, 64% connected' indicates a true label of connected (first connected) and the connected (second connected) category got the largest predicted probability and that probability was 64%. Given an argmax assignment, this is a correct prediction.<jupyter_code>visualize_testing_predictions(model)                    # visualize testing images and probability for correct category<jupyter_output><empty_output><jupyter_text>Summarize the Performance of our Image Classification CNNLet's summarize the performance of our CNN over all training and testing images * we use boxplots to report the training and testing images prediction probabilities * we show the predicted probability of connected or not connected geobody images<jupyter_code>prediction_boxplot(model)                                  # visualize the prediction probability box plots by correct category<jupyter_output><empty_output><jupyter_text>To better understand our model, we are going to test it's abilities by running two test cases on it Test Case 1: Inverted Images<jupyter_code>names = ['connected', 'not connected']
for i in range(0,len(X_test)):
    array = X_test[i,:,:,0]
    X_test[i,:,:,0] = np.where(array < .15, 1, 0)  

for i in range(0,len(X_test)):
    array = X_test[i,:,:,0]
    new_array = np.where(array < .15, 0, 1)
    new_array = new_array.tolist()
    if maxAreaOfIsland(new_array) > 72:
        y_test[i] = 0
    else:
        y_test[i] = 1

visualize_testing_predictions(model)<jupyter_output><empty_output><jupyter_text>Test Case 2: Images that lie outside of our original data's sand proportions<jupyter_code>XX = np.zeros([4,nx,ny,1])
yy_test = np.zeros(ntest,dtype = int)
c, r = 4, 1                                             # visualize labelled synthetic images
XX[0,0:12,4:6,0] = 1
XX[1,4:6,0:12,0] = 1
XX[2,2:10,0:12,0] = 1
XX[3,0:12,2:10,0] = 1
for i in range(0,len(XX)):
    array = XX[i,:,:,0]
    new_array = np.where(array < .15, 0, 1)
    new_array = new_array.tolist()
    if maxAreaOfIsland(new_array) > 72:
        yy_test[i] = 0
    else:
        yy_test[i] = 1
for i in range(c*r):
    plt.subplot(r,c,i+1)
    fig = plt.imshow(XX[i,:,:,0],cmap=cmap)
    fig.axes.get_xaxis().set_visible(False); fig.axes.get_yaxis().set_visible(False)
    plt.title(names[yy_test[i]])

plt.subplots_adjust(left=0.0, bottom=0.0, right=1.3, top=0.5, wspace=.5, hspace=0.3)
plt.show()
def visualize_testing_predictions2(model):               # visualize testing images and predicted probabilities 
    y_test_predict = model.predict(XX[:,:,:,0].reshape([4,nx,ny,-1])) # predict over the testing images 
    y_test_predict_cat = np.argmax(y_test_predict,axis = 1) # assign the maximum probability category
    
    c, r = 4, 1                                         
    plt.figure(figsize = (15,15))
    for i in range(c*r):
        plt.subplot(r,c,i+1)
        fig = plt.imshow(XX[i,:,:,0],cmap=cmap)
        fig.axes.get_xaxis().set_visible(False); fig.axes.get_yaxis().set_visible(False)
        icat = y_test_predict_cat[i]
        plt.title('Image ' + str(i) + ': ' + names[yy_test[i]] + ', ' + str(np.round(y_test_predict[i,icat]*100,2)) + "% " + names[icat] )

            
    plt.subplots_adjust(left=0.0, bottom=0.0, right=1.3, top=1.6, wspace=0.01, hspace=0.3)
    plt.show()
visualize_testing_predictions2(model)<jupyter_output><empty_output><jupyter_text>Conclusion:Although the model is able to distinguish porous media from non-porous in the original images, it struggles with edge detection around the sides of the synthetic images. Specifically, in Test Case 1, the model failed to correctly identify all of the  connected images as all the porous (yellow) cells lied on the exterior of the images. This highlights the importance of using padding when creating a model. The model performed well with images *exceeding* the train and test data sand proportions. This proves that the proportions of porous media have little to no effect on the model's effectivness. Convolutional neural networks such as this one are applicable to the oil and gas industry as they can save time and money when searching for profitable reservoirs.  CommentsThis was a very simple convolutional neural network workflow to support experiential learning with fast experientation.  The Texas Center for Data Analytics and Geostatistics has many other demonstrations on the basics of working with DataFrames, ndarrays, univariate statistics, plotting data, declustering, data transformations, trend modeling and many other workflows available [here](https://github.com/GeostatsGuy/PythonNumericalDemos), along with a package for geostatistics in Python called [GeostatsPy](https://github.com/GeostatsGuy/GeostatsPy).   We hope this was helpful,*<NAME>,* and *Michael**** More on <NAME> and the Texas Center for Data Analytics and Geostatistics: <NAME>, 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 <EMAIL>.I'm always happy to discuss,*Michael*<NAME>, 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)<jupyter_code>print(list(history.history.keys())[1])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11698650241266594, 'fraction_numerical': 0.024860465633408196, 'mean_word_length': 3.850825153705102, 'pattern_counts': {'":': 0, '<': 105, '<?xml version=': 0, '>': 118, 'https://': 35, 'lorem ipsum': 0, 'www.': 14, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 44, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '225483', 'license': 'MIT', 'n_tokens_mistral': 14806, 'n_tokens_neox': 13905, 'n_words': 4023, 'path': 'SubsurfaceDataAnalytics_ConvolutionalNeuralNetworks_Percolation_Threshold.ipynb', 'repo_name': 'GeostatsGuy/PythonNumericalDemos', 'size': '45053'}