<jupyter_start><jupyter_text>**Deeplearning in Prediction of Geological Facies(Convolution Neural Nets)** Prediction of geological surface using available data found in the github [seg](www.github.c0m/seg) repository<jupyter_code>from google.colab import drive
drive.mount('/content/gdrive')
# ls '/content/gdrive/My Drive/Data'
%matplotlib inline
import pandas as pd
import numpy as np
import time
import warnings
warnings.filterwarnings('ignore')
from pandas import set_option

from sklearn import preprocessing
from sklearn.metrics import confusion_matrix, f1_score,accuracy_score

from collections import Counter
import operator

from keras.models import Model, Sequential
from keras.layers import Convolution1D, Convolution2D, Dense, Input, Dropout, Flatten, MaxPooling2D, Activation, normalization
from keras.optimizers import Nadam
from keras.utils import np_utils
from keras.utils.np_utils import to_categorical
from keras.callbacks import History

from keras.wrappers.scikit_learn import KerasClassifier
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import KFold, StratifiedKFold
np.random.seed(123)
def accuracy(conf):
    total_correct = 0.
    nb_classes = conf.shape[0]
    for i in np.arange(0,nb_classes):
        total_correct += conf[i][i]
    acc = total_correct/sum(sum(conf))
    return acc

def label_facies(row, labels):
    return labels[ row['Facies'] -1]
  
def accuracy_adj(conf, adjacent_facies):
  total_correct = 0
  nb_classes = conf.shape[0]
  for i in np.arange(0, nb_classes):
    total_correct +=conf[i][i]
    for j in adjacent_facies[i]:
      total_correct += conf[i][j]
  return total_correct/sum(sum(conf))

file = pd.read_csv('/content/gdrive/My Drive/Data/train_test_data.csv')
print(file.isna().mean().sort_values(ascending=False))
file.head()<jupyter_output>PE           0.184174
RELPOS       0.000000
NM_M         0.000000
PHIND        0.000000
DeltaPHI     0.000000
               ...   
GR           0.000000
Depth        0.000000
Well Name    0.000000
Formation    0.000000
Facies       0.000000
Length: 11, dtype: float64<jupyter_text>Prediction of PE values<jupyter_code>file[file['PE'].notnull()]
def prepare_feature_vectors(data, features, window_width):

    raw_feature_vectors = data[features]
    well_labels = data['Well Name']
    num_features = np.shape(raw_feature_vectors)[1]

    output = np.zeros((1, window_width, num_features))
    for x in well_labels.unique():
        well = raw_feature_vectors[well_labels == x].values
        well = np.concatenate((np.repeat(well[0:1], np.floor((window_width-1)/2.0), axis=0), well,
                              np.repeat(well[-1:None], np.floor(window_width/2.0), axis=0)), axis=0)

        tmp = np.zeros((np.size(well, axis=0) - window_width + 1, window_width, num_features))
        for i in np.arange(np.size(well, axis=0) - window_width + 1):
            tmp[i] = np.reshape(well[i: i + window_width], (window_width, num_features))

        output = np.append(output, tmp, axis=0)

    return output[1:]<jupyter_output><empty_output><jupyter_text>Remove the P.E values of Stuart and Crawford Wells<jupyter_code>def file_copy(file):
  data = file.copy()
  data = data[data['Well Name'] != 'STUART']
  data = data[data['Well Name'] != 'CRAWFORD']
  return data<jupyter_output><empty_output><jupyter_text>Split Data into Training and Testing<jupyter_code>data['PE'].isnull().sum()
window_width = 15 #Window around the central value
data = file_copy(file)
feature_list = ['GR', 'ILD_log10', 'DeltaPHI', 'PHIND', 'NM_M', 'RELPOS']
X_train = prepare_feature_vectors(data[data['PE'].notnull()], feature_list, window_width)
num_train_samples = np.asarray(np.shape(X_train))[0]

X_test = prepare_feature_vectors(data[data['PE'].isnull()], feature_list, window_width)
num_test_samples = np.asarray(np.shape(X_test))[0]

print('Training Samples=', num_train_samples, '   Test Samples=', num_test_samples)<jupyter_output>Training Samples= 3232    Test Samples= 917<jupyter_text>Neural Network to Predict the Missing PE Values<jupyter_code># define neural network to perform regression on PE
num_filters = 12
dropout_prob = 0.6666

cnn = Sequential()
cnn.add(Convolution1D(num_filters, 1, border_mode='valid', input_shape=(window_width, len(feature_list))))
cnn.add(normalization.BatchNormalization())
cnn.add(Activation('tanh'))
cnn.add(Convolution1D(num_filters, 3, border_mode='valid'))
cnn.add(normalization.BatchNormalization())
cnn.add(Activation('tanh'))
cnn.add(Dropout(dropout_prob / 2))

cnn.add(Flatten())
cnn.add(Dense(4*num_filters))
cnn.add(normalization.BatchNormalization())
cnn.add(Activation('tanh'))
cnn.add(Dropout(dropout_prob))

cnn.add(Dense(1))
cnn.compile(loss='mse', optimizer='rmsprop', metrics=['acc'])
cnn.summary()

# save initial weights, which are random
initial_weights = cnn.get_weights()
Y_train = data[data['PE'].notnull()]['PE'].values
# define training parameters and prepare arrays to store training metrics
epochs_per_fold = 1000
num_fold = 5
roll_stride = np.ceil(num_train_samples/num_fold).astype(int)

cnn_hist = History()
hist = np.zeros((4, num_fold, epochs_per_fold))
f1scores = np.zeros(num_fold)
Y_test = np.zeros((num_test_samples, num_fold))


# shuffle input data
rand_perm = np.random.permutation(num_train_samples)
X_train = X_train[rand_perm]
Y_train = Y_train[rand_perm]
# use 5-fold cross validation and train 5 neural networks, ending up with 5 sets of predictions
for i in np.arange(num_fold):
    cnn.set_weights(initial_weights)
    X_train = np.roll(X_train, i*roll_stride, axis=0)
    Y_train = np.roll(Y_train, i*roll_stride, axis=0)

    cnn.fit(X_train, Y_train, batch_size=150, nb_epoch=epochs_per_fold, verbose=0,
                validation_split=1.0/num_fold, callbacks=[cnn_hist])

    # make predictions, i.e. impute PE
    Y_test[:, i] = cnn.predict(X_test)[:, 0]

    hist[:, i, :] = [cnn_hist.history['acc'], cnn_hist.history['val_acc'],
                     cnn_hist.history['loss'], cnn_hist.history['val_loss']]
    print("Accuracy  =", np.mean(hist[1, i, -100:]))
import matplotlib.pyplot as plt
# plot callbacks to evaluate quality of training
drop_values = 100
drop_hist = np.reshape(hist[:, :, drop_values:], (4, num_fold * (epochs_per_fold - drop_values)))
print("Mean Validation Accuracy  =", np.mean(hist[1, :, -drop_values:]))

plt.plot(drop_hist[0]); plt.plot(drop_hist[1])
plt.legend(['train', 'val'], loc='lower left')
plt.plot(drop_hist[2]); plt.plot(drop_hist[3])
plt.legend(['train', 'val'], loc='upper left')
# Update dataframe with imputed values by averaging results of 5 neural networks 
data['PE'][np.array(data['PE'].isnull())] = np.mean(Y_test, axis=1)

# Write intermediate data to file
data.to_csv('/content/gdrive/My Drive/Data/train_test_data_PEfilled.csv', index=False)
# Plot PE of all wells (original plus imputed) 
plt.plot(data['PE'])
data['PE'].isnull().sum() #No missing PE Values<jupyter_output><empty_output><jupyter_text>Facies Prediction<jupyter_code>file_2 = file.copy()
names = file['Well Name'].unique()
for name in names:
  if name not in ['STUART','CRAWFORD']:
    file_2 = file_2[file_2['Well Name'] != name]
file_2.shape
file_ii = pd.read_csv('/content/gdrive/My Drive/Data/train_test_data_PEfilled.csv')
training_data = pd.concat([file_ii, file_2], axis = 0)
print(training_data.shape)
training_data.head()<jupyter_output>(4979, 11)<jupyter_text>Data Wrangling<jupyter_code>set_option("display.max_rows", 10)
pd.options.mode.chained_assignment = None

training_data['Well Name'] = training_data['Well Name'].astype('category')
training_data['Formation'] = training_data['Formation'].astype('category')
training_data['Well Name'].unique()

# 1=sandstone  2=c_siltstone   3=f_siltstone 
# 4=marine_silt_shale 5=mudstone 6=wackestone 7=dolomite
# 8=packstone 9=bafflestone

facies_colors = ['#F4D03F', '#F5B041','#DC7633','#6E2C00',
       '#1B4F72','#2E86C1', '#AED6F1', '#A569BD', '#196F3D']

facies_labels = ['SS', 'CSiS', 'FSiS', 'SiSh', 'MS','WS', 'D','PS', 'BS']
#facies_color_map is a dictionary that maps facies labels
#to their respective colors
facies_color_map = {}
for ind, label in enumerate(facies_labels):
    facies_color_map[label] = facies_colors[ind]
  #mapping facies labels to the dataframe
training_data.loc[:,'FaciesLabels'] = training_data.apply(lambda row: label_facies(row, facies_labels), axis=1)

correct_facies_labels = training_data['Facies'].values
feature_vectors = training_data.drop(['Formation', 'FaciesLabels', 'Facies','Well Name'], axis=1) #, 'RELPOS', 'NM_M', 'Depth', 'ILD_log10',  'DeltaPHI',   'PHIND'], axis=1)

well_labels = training_data[['Well Name', 'Facies']].values
well_name = training_data['Well Name'].values
depth = training_data['Depth'].values

training_data[:3]

scaler = preprocessing.StandardScaler().fit(feature_vectors)
scaled_features = scaler.transform(feature_vectors)

data_out = np.hstack([well_labels, scaled_features])<jupyter_output><empty_output><jupyter_text>Data Transform<jupyter_code>data = data_out
well_data = {}
well_names = list(set(data[:, 0])) #set of well names
for name in well_names:
    well_data[name] = [[], []] #creates a dictionary containing list of list for the well names
    
for row in data:
    well_data[row[0]][1].append(row[1]) #Append facies to the second list[]
    well_data[row[0]][0].append(list(row[2::])) #Append features to the first list[]

#For data padding
positive_lag = 10
negative_lag = 11

chunks_cnn = []
chunks_cnn_test = []
chunk_length = positive_lag+negative_lag+1 #were gonna predict middle facies
chunks_facies_cnn = []

for name in well_names:
    if name not in ['STUART', 'CRAWFORD']:  #Training set augumentation
        test_well_data = well_data[name]    #Returns a list of list
        log_values = np.array(test_well_data[0]) #values
        log_values_padded = np.lib.pad(log_values, (negative_lag,positive_lag), 'edge')[:, negative_lag:-positive_lag] #Pad values(11,10) and 
                                                                                        #returns index from 11:-10
        facies_values =  np.array(test_well_data[1])
        for i in range(log_values.shape[0]):
            chunk = log_values_padded[i:i+chunk_length, :]
            chunk_trans = chunk.T
            chunks_cnn.append(chunk_trans)
            chunks_facies_cnn.append(facies_values[i])
    else:
        test_well_data = well_data[name]
        log_values = np.array(test_well_data[0])
        log_values_padded = np.lib.pad(log_values, (negative_lag,positive_lag), 'edge')[:, negative_lag:-positive_lag]
        facies_values =  np.array(test_well_data[1])
        for i in range(log_values.shape[0]):
            chunk = log_values_padded[i:i+chunk_length, :]
            chunk_trans = chunk.T
            chunks_cnn_test.append(chunk_trans)

chunks_cnn = np.array(chunks_cnn)
chunks_cnn_test = np.array(chunks_cnn_test)

chunks_facies_cnn = np.array(chunks_facies_cnn, dtype=np.int32)-1

unique_facies = len(set(chunks_facies_cnn))
print (unique_facies, set(chunks_facies_cnn))
print (chunks_cnn.shape, chunks_cnn_test.shape)
print (chunks_facies_cnn.shape, chunks_cnn_test.shape)
X = chunks_cnn #Save data to X and Y for training
y = chunks_facies_cnn

X = X.reshape((chunks_cnn.shape[0], chunks_cnn.shape[1], chunks_cnn.shape[2], 1))

y = np_utils.to_categorical(y)
N = 128
cnnet = Sequential()
cnnet.add(Convolution2D(N, 1, 5, border_mode="same",activation="relu",input_shape=(chunks_cnn.shape[1], chunks_cnn.shape[2], 1)))
cnnet.add(MaxPooling2D(pool_size=(1, 2)))
cnnet.add(Dropout(0.25))
cnnet.add(Convolution2D(N, 1, 3, border_mode="same",activation="relu",input_shape=(chunks_cnn.shape[1], chunks_cnn.shape[2], 1)))
cnnet.add(MaxPooling2D(pool_size=(1, 2)))
#cnn.add(Dropout(0.5))
cnnet.add(Convolution2D(N, 2, 2, border_mode="same", activation="relu"))
#cnn.add(Convolution2D(N, 3, 1, border_mode="same", activation="relu"))
cnnet.add(MaxPooling2D(pool_size=(2, 2)))
cnnet.add(Dropout(0.8))
cnnet.add(Flatten())
cnnet.add(Dense(128, activation="relu"))
cnnet.add(Dropout(0.5))
cnnet.add(Dense(64, activation="relu"))
cnnet.add(Dropout(0.5))
cnnet.add(Dense(9, activation="softmax"))
cnnet.compile(loss="categorical_crossentropy", optimizer="Adamax", metrics=['acc'])
cnnet.summary()
intial_weights_ii = cnnet.get_weights()
cnnet_hist = History()
epoch_fold = 200
hist = np.zeros((4, 1, epoch_fold))
start = time.time()
cnnet.fit(X, y, nb_epoch=epoch_fold, validation_split=0.25, batch_size=64, verbose=1, shuffle=True, callbacks=[cnnet_hist])
hist[:,0,:] = [cnnet_hist.history['acc'], cnnet_hist.history['val_acc'], cnnet_hist.history['loss'], cnnet_hist.history['val_loss']]
print(f'Time to train Network {time.time() - start}')
hist = np.reshape(hist,(4,i*epoch_fold))
plt.plot(hist[0]); plt.plot(hist[1])
plt.legend(['train','val'], loc= 'lower left')
cnn.get_weights()
plt.plot(hist[2]); plt.plot(hist[3])
plt.legend(['train', 'val'], loc='upper left')<jupyter_output><empty_output><jupyter_text>Model Predictions<jupyter_code>y_predicted = cnnet.predict(X, batch_size=32, verbose=1)

y_preds = []
for row in y_predicted:
    index, value = max(enumerate(row), key=operator.itemgetter(1))
    y_preds.append(index)

print ("" )   
print (confusion_matrix(chunks_facies_cnn, y_preds))
print ('Model F1 score',f1_score(chunks_facies_cnn, y_preds, average='weighted'))
print('Model Accuracy Score',accuracy_score(chunks_facies_cnn,y_preds))
conf  = confusion_matrix(chunks_facies_cnn,y_preds)

adjacent_facies = np.array([[1], [0,2], [1], [4], [3,5], [4,6,7], [5,7], [5,6,8], [6,7]])
accuracy_adj(conf,adjacent_facies)
X_test = chunks_cnn_test

X_test = X_test.reshape((chunks_cnn_test.shape[0], chunks_cnn_test.shape[1], chunks_cnn_test.shape[2], 1))


y_predicted = cnnet.predict(X_test, batch_size=32, verbose=1)

y_preds = []
for row in y_predicted:
    index, value = max(enumerate(row), key=operator.itemgetter(1))
    y_preds.append(index)
y_preds = np.array(y_preds)+1   
print(y_preds)
test_data = pd.read_csv("/content/gdrive/My Drive/Data/validation_data_nofacies_online.csv")
test_data['Facies'] = pd.Series(y_preds)
test_data.to_csv("/content/gdrive/My Drive/Data/validation_data_with_facies_new.csv")
print (test_data.head())
set(y_preds)<jupyter_output><empty_output><jupyter_text>Prediction on Blind Data<jupyter_code>blind_data = pd.read_csv('/content/gdrive/My Drive/Data/blind_stuart_crawford_core_facies.csv')
blind_data.head()
def merge_test(blind_data, test_data):
  columns = blind_data.WellName.unique()
  merge_col = 0
  for column in columns:
    blind = blind_data[blind_data['WellName'] == column]
    test = test_data[test_data['Well Name']==column]
    if type(merge_col) is int:
      merge_col = pd.merge(blind, test, how = 'inner', left_on = 'Depth.ft', right_on= 'Depth')
    else:
      merge_col_2 = pd.merge(blind, test, how = 'inner', left_on = 'Depth.ft', right_on= 'Depth')
  y_true = np.concatenate((merge_col.LithCode.values, merge_col_2.LithCode.values), axis = 0)
  y_pred = np.concatenate((merge_col.Facies.values, merge_col_2.Facies.values), axis = 0)
  return accuracy_score(y_true, y_pred)*100
merge_test(blind_data, test_data)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14312989663572323, 'fraction_numerical': 0.024293896899071698, 'mean_word_length': 5.55023717119448, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '280753', 'license': 'Apache-2.0', 'n_tokens_mistral': 5780, 'n_tokens_neox': 5539, 'n_words': 1054, 'path': 'Deeplearning_CNN.ipynb', 'repo_name': 'Arnoldsynchron/2016-ml-contest', 'size': '15189'}