Title, authors, credit...
In [1]:
%matplotlib notebook
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import tensorflow as tf
import random
from sklearn.kernel_ridge import KernelRidge
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import f1_score, confusion_matrix
import classification_utilities as utils
from ConvNet import inception_net
In [2]:
FROM_SCRATCH = False
In [3]:
# Some parameters
#feature_names = ['GR', 'ILD_log10', 'DeltaPHI', 'PHIND', 'PE', 'NM_M', 'RELPOS', 'Depth']
feature_names = ['GR', 'ILD_log10', 'DeltaPHI', 'PHIND', 'PE', 'NM_M', 'RELPOS']
facies_names = ['SS', 'CSiS', 'FSiS', 'SiSh', 'MS', 'WS', 'D', 'PS', 'BS']
facies_colors = ['#F4D03F', '#F5B041','#DC7633','#6E2C00', '#1B4F72','#2E86C1', '#AED6F1', '#A569BD', '#196F3D']
approximate_neighbors = [[2], [1, 3], [2], [5], [4, 6], [5, 7, 8], [6, 8, 9], [6, 7, 9], [7, 8]]
# Load data from file
training_data = pd.read_csv('../facies_vectors.csv').drop('Formation', axis=1)
pd.set_option("display.max_rows", 4, "float_format", lambda v: "%.2f" % v )
training_data
Out[3]:
In [4]:
# Well labels
well_name = training_data['Well Name']
# Features and labels
X = training_data[feature_names].values
y = training_data['Facies'].values
print(X.shape) # features
print(y.shape) # labels
In [5]:
# Not sure if good idea to use F9 well?
well_to_be_used = []
for well in well_name.unique():
if well != 'Recruit F9':
well_to_be_used.append(well)
random.shuffle(well_to_be_used)
#well_to_be_used.append('Recruit F9')
In [6]:
if FROM_SCRATCH is True:
# Fit KernelRidge with parameter selection based on 5-fold cross validation
krr_param_grid = {"alpha": [1e0, 1e-1, 1e-2, 1e-3],
"kernel": ['laplacian'],
"gamma": np.linspace(0.01,1,10)}
krr = GridSearchCV(KernelRidge(), cv=5, param_grid=krr_param_grid)
data_missing_vals = training_data[feature_names].copy()
data_no_missing_vals = data_missing_vals.dropna(axis = 0, inplace=False)
f_ = data_no_missing_vals.loc[:, data_no_missing_vals.columns != 'PE']
l_ = data_no_missing_vals.loc[:, 'PE']
krr.fit(f_, l_)
X[np.array(data_missing_vals.PE.isnull()),
4] = krr.predict(data_missing_vals.loc[data_missing_vals.PE.isnull(),:].drop('PE',axis=1,inplace=False))
else:
del X
X = np.load('X_after_krr.npy')
In [ ]:
In [7]:
training_data_imp = pd.DataFrame(X,columns=feature_names, index=well_name)
training_data_imp['Facies'] = training_data['Facies'].values
training_data_imp['Depth'] = training_data['Depth'].values
In [8]:
X.shape
Out[8]:
In [9]:
display_step = 10
batch_size = 32
training_epochs = 71#71
sequence_length = 31
perc_overlap = 30
perturbation_ratio = 0.08
initial_learning_rate = 0.001
epoch_decrease_every = 25
decrease_factor = 0.96
beta = .02#0.02
clip_norm=1.5
dropout_fc = 0.5 # probability to keep units in fully connected layers
dropout_conv = 0.5 # probability to keep units in fully conv layers
In [10]:
n_runs = 3
f1_blind = {}
# loop over blind well
for blind_well in ['NEWBY']:#well_to_be_used:
print(blind_well)
f1_blind[blind_well] = []
for i in range(n_runs):
tf.reset_default_graph()
print('random realisation #%d' % (i+1))
# training wells
wells_for_training = [well for well in well_to_be_used if well != blind_well]
# Generate sequences
TST_X, TST_Y, scaler = utils.generate_sequences(training_data_imp, wells_for_training,feature_names,
sequence_length=sequence_length,perc_overlap=perc_overlap,
flip_ud=False)
Y_blind = utils.generate_sequences(training_data_imp, [blind_well],feature_names,
sequence_length=sequence_length,perc_overlap=perc_overlap,
flip_ud=False)[1]
X_blind, Y_blind_absolute_idx_pos = utils.generate_sequences(training_data_imp, [blind_well],feature_names,
sequence_length=sequence_length,
with_labels=False,perc_overlap=perc_overlap,
flip_ud=False, scaler=scaler)
# Take some data from training set for intermediate testing
X_train, Y_train, X_test, Y_test = utils.split_training_testing(TST_X, TST_Y, percentage_for_testing=10)
total_batch = X_train.shape[0]//batch_size
# Parameters
# learning rate
global_step = tf.Variable(0, trainable=False)
initial_learning_rate = initial_learning_rate
lr = tf.train.exponential_decay(initial_learning_rate, global_step, epoch_decrease_every*total_batch,
decrease_factor, staircase=True)
# Network i/o
n_input = sequence_length
n_features = len(feature_names)
n_classes = len(facies_names)
x_net = tf.placeholder(tf.float32, [None, 1, n_input, n_features])
y_net = tf.placeholder(tf.float32, [None, n_classes])
keep_prob_fc = tf.placeholder(tf.float32) #dropout
keep_prob_conv = tf.placeholder(tf.float32)
is_training = tf.placeholder(tf.bool)
# Construct model
pred = inception_net(x_net, keep_prob_fc, is_training, dropout_conv=keep_prob_conv, clip_norm=clip_norm)
# Loss and optimizer
conv_weights = [var for var in tf.global_variables() if 'kernel' in var.name]
weights_constraint = tf.reduce_sum(beta*tf.pack([tf.nn.l2_loss(w) for w in conv_weights]))
loss = tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits(logits=pred, labels=y_net))
cost = tf.add(loss,weights_constraint)
optimizer = tf.train.AdamOptimizer(learning_rate=lr).minimize(cost, global_step=global_step)
# Evaluate model
correct_pred = tf.equal(tf.argmax(pred, 1), tf.argmax(y_net, 1))
accuracy = tf.reduce_mean(tf.cast(correct_pred, tf.float32))
# Initializing the variables
init = tf.global_variables_initializer()
# Launch the graph
with tf.Session() as sess:
sess.run(init)
# Training cycle
for epoch in range(training_epochs):
X_train, Y_train, X_test, Y_test = utils.split_training_testing(TST_X, TST_Y,
percentage_for_testing=10)
# shuffle training and testing set
idx_train = np.random.permutation(np.arange(X_train.shape[0]))
idx_test = np.random.permutation(np.arange(X_test.shape[0]))
x_train = X_train[idx_train,...]
y_train = Y_train[idx_train]
x_test = X_test[idx_test,...]
y_test = Y_test[idx_test]
# perturbe training labels
y_train = utils.perturbe_label(y_train, approximate_neighbors, perturbation_ratio=perturbation_ratio)
# Loop over all batches
for i in range(total_batch):
batch_idx = range(i*batch_size,(i+1)*batch_size)
batch_x = x_train[batch_idx,...]
batch_y = utils.dense_to_one_hot(y_train[batch_idx],n_classes)
# Run optimization op (backprop)
sess.run(optimizer, feed_dict={x_net: batch_x,
y_net: batch_y,
keep_prob_fc: dropout_fc,
keep_prob_conv:dropout_conv,
is_training:True})
# Calculate loss and accuracy
if epoch % display_step == 0:
#rand_idx = np.random.permutation(np.arange(x_test.shape[0]))[:batch_size]
#batch_x_test = x_test[rand_idx,...]
#batch_y_test = utils.dense_to_one_hot(y_test[rand_idx],n_classes)
loss, acc = sess.run([cost, accuracy], feed_dict={x_net: x_test,
y_net: utils.dense_to_one_hot(y_test,n_classes),
keep_prob_fc: 1.,
keep_prob_conv:1.,
is_training:False})
print("Epoch:", '%04d' % (epoch+1), ", Minibatch Loss= " , \
"{:.6f}".format(loss) , ", Training Accuracy= " , \
"{:.5f}".format(acc))
# Calculate accuracy for blind well
acc , y_blind_pred = sess.run([accuracy, pred], feed_dict={x_net: X_blind,
y_net: utils.dense_to_one_hot(Y_blind,n_classes),
keep_prob_fc: 1.,
keep_prob_conv:1.,
is_training:False})
y_blind_pred_sm = sess.run(tf.nn.softmax(y_blind_pred))
y_blind_pred_final = utils.redundant_to_final(training_data_imp.ix[blind_well]['Facies'].shape[0],
utils.one_hot_to_dense(y_blind_pred_sm),
Y_blind_absolute_idx_pos)
y_blind_truth = training_data_imp.ix[blind_well]['Facies'].values
f1_blind[blind_well].append(f1_score(y_blind_truth, y_blind_pred_final, average='micro'))
In [11]:
y_blind_pred.shape
Out[11]:
In [12]:
print('Number of runs per well: %d' % n_runs)
print('\n')
for well in f1_blind:
mean = np.mean(np.array(f1_blind[well]))
std = np.std(np.array(f1_blind[well]))
median = np.median(np.array(f1_blind[well]))
minimum = np.min(np.array(f1_blind[well]))
maximum = np.max(np.array(f1_blind[well]))
print('%s f1 scores: ' % well)
print('mean: %f, std: %f, median: %f, min: %f, max:%f' % (mean, std, median, minimum, maximum))
print('\n')
In [13]:
#fig, ax = plt.subplots()
#ax.imshow(X_train[np.random.randint(X_train.shape[0]),0,:,:], aspect='auto')