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

Data


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]:
Facies Well Name Depth GR ILD_log10 DeltaPHI PHIND PE NM_M RELPOS
0 3 SHRIMPLIN 2793.00 77.45 0.66 9.90 11.91 4.60 1 1.00
1 3 SHRIMPLIN 2793.50 78.26 0.66 14.20 12.56 4.10 1 0.98
... ... ... ... ... ... ... ... ... ... ...
4147 5 CHURCHMAN BIBLE 3122.00 51.47 0.96 3.08 7.71 3.15 2 0.66
4148 5 CHURCHMAN BIBLE 3122.50 50.03 0.97 2.61 6.67 3.29 2 0.65

4149 rows × 10 columns


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


(4149, 7)
(4149,)

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]:
(4149, 7)

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


NEWBY
random realisation #1
Epoch: 0001 , Minibatch Loss=  1.506495 , Training Accuracy=  0.42085
Epoch: 0011 , Minibatch Loss=  1.048651 , Training Accuracy=  0.62162
Epoch: 0021 , Minibatch Loss=  0.928421 , Training Accuracy=  0.68919
Epoch: 0031 , Minibatch Loss=  0.887138 , Training Accuracy=  0.72780
Epoch: 0041 , Minibatch Loss=  0.836614 , Training Accuracy=  0.73745
Epoch: 0051 , Minibatch Loss=  0.843847 , Training Accuracy=  0.72973
Epoch: 0061 , Minibatch Loss=  0.807839 , Training Accuracy=  0.74517
Epoch: 0071 , Minibatch Loss=  0.817065 , Training Accuracy=  0.74710
random realisation #2
Epoch: 0001 , Minibatch Loss=  1.438035 , Training Accuracy=  0.45174
Epoch: 0011 , Minibatch Loss=  1.051913 , Training Accuracy=  0.65058
Epoch: 0021 , Minibatch Loss=  0.939903 , Training Accuracy=  0.68533
Epoch: 0031 , Minibatch Loss=  0.917730 , Training Accuracy=  0.69305
Epoch: 0041 , Minibatch Loss=  0.905735 , Training Accuracy=  0.67375
Epoch: 0051 , Minibatch Loss=  0.840623 , Training Accuracy=  0.70270
Epoch: 0061 , Minibatch Loss=  0.834898 , Training Accuracy=  0.72201
Epoch: 0071 , Minibatch Loss=  0.798336 , Training Accuracy=  0.71429
random realisation #3
Epoch: 0001 , Minibatch Loss=  1.450656 , Training Accuracy=  0.48069
Epoch: 0011 , Minibatch Loss=  1.019614 , Training Accuracy=  0.62355
Epoch: 0021 , Minibatch Loss=  0.979476 , Training Accuracy=  0.65830
Epoch: 0031 , Minibatch Loss=  0.916164 , Training Accuracy=  0.66409
Epoch: 0041 , Minibatch Loss=  0.879585 , Training Accuracy=  0.70077
Epoch: 0051 , Minibatch Loss=  0.893836 , Training Accuracy=  0.67954
Epoch: 0061 , Minibatch Loss=  0.803441 , Training Accuracy=  0.73166
Epoch: 0071 , Minibatch Loss=  0.834823 , Training Accuracy=  0.72008

In [11]:
y_blind_pred.shape


Out[11]:
(651, 9)

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


Number of runs per well: 3


NEWBY f1 scores: 
mean: 0.570194, std: 0.013773, median: 0.578834, min: 0.550756, max:0.580994



In [13]:
#fig, ax = plt.subplots()
#ax.imshow(X_train[np.random.randint(X_train.shape[0]),0,:,:], aspect='auto')