<jupyter_start><jupyter_text>Facies classification using machine learning techniquesCopy of <NAME>'s "Try 2", augmented, by <NAME> (Ausar Geophysical), with an ML estimator for PE in the wells where it is missing (rather than just using the mean).In the following, we provide a possible solution to the facies classification problem described at https://github.com/seg/2016-ml-contest.The proposed algorithm is based on the use of random forests combined in one-vs-one multiclass strategy. In particular, we would like to study the effect of:- Robust feature normalization.- Feature imputation for missing feature values.- Well-based cross-validation routines.- Feature augmentation strategies. Script initializationLet us import the used packages and define some parameters (e.g., colors, labels, etc.).<jupyter_code># Import
from __future__ import division
%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rcParams['figure.figsize'] = (20.0, 10.0)
inline_rc = dict(mpl.rcParams)
from classification_utilities import make_facies_log_plot

import pandas as pd
import numpy as np
#import seaborn as sns

from sklearn import preprocessing
from sklearn.model_selection import LeavePGroupsOut
from sklearn.metrics import f1_score
from sklearn.multiclass import OneVsOneClassifier
from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor

from scipy.signal import medfilt
import sys, scipy, sklearn
print('Python:  ' + sys.version.split('\n')[0])
print('         ' + sys.version.split('\n')[1])
print('Pandas:  ' + pd.__version__)
print('Numpy:   ' + np.__version__)
print('Scipy:   ' + scipy.__version__)
print('Sklearn: ' + sklearn.__version__)
# Parameters
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']<jupyter_output><empty_output><jupyter_text>Load dataLet us load training data and store features, labels and other data into numpy arrays.<jupyter_code># Load data from file
data = pd.read_csv('../facies_vectors.csv')
# Store features and labels
X = data[feature_names].values  # features
y = data['Facies'].values  # labels
# Store well labels and depths
well = data['Well Name'].values
depth = data['Depth'].values<jupyter_output><empty_output><jupyter_text>Data inspectionLet us inspect the features we are working with. This step is useful to understand how to normalize them and how to devise a correct cross-validation strategy. Specifically, it is possible to observe that:- Some features seem to be affected by a few outlier measurements.- Only a few wells contain samples from all classes.- PE measurements are available only for some wells.<jupyter_code># Define function for plotting feature statistics
def plot_feature_stats(X, y, feature_names, facies_colors, facies_names):
    
    # Remove NaN
    nan_idx = np.any(np.isnan(X), axis=1)
    X = X[np.logical_not(nan_idx), :]
    y = y[np.logical_not(nan_idx)]
    
    # Merge features and labels into a single DataFrame
    features = pd.DataFrame(X, columns=feature_names)
    labels = pd.DataFrame(y, columns=['Facies'])
    for f_idx, facies in enumerate(facies_names):
        labels[labels[:] == f_idx] = facies
    data = pd.concat((labels, features), axis=1)

    # Plot features statistics
    facies_color_map = {}
    for ind, label in enumerate(facies_names):
        facies_color_map[label] = facies_colors[ind]

    sns.pairplot(data, hue='Facies', palette=facies_color_map, hue_order=list(reversed(facies_names)))
# Feature distribution
# plot_feature_stats(X, y, feature_names, facies_colors, facies_names)
# mpl.rcParams.update(inline_rc)
# Facies per well
for w_idx, w in enumerate(np.unique(well)):
    ax = plt.subplot(3, 4, w_idx+1)
    hist = np.histogram(y[well == w], bins=np.arange(len(facies_names)+1)+.5)
    plt.bar(np.arange(len(hist[0])), hist[0], color=facies_colors, align='center')
    ax.set_xticks(np.arange(len(hist[0])))
    ax.set_xticklabels(facies_names)
    ax.set_title(w)
# Features per well
for w_idx, w in enumerate(np.unique(well)):
    ax = plt.subplot(3, 4, w_idx+1)
    hist = np.logical_not(np.any(np.isnan(X[well == w, :]), axis=0))
    plt.bar(np.arange(len(hist)), hist, color=facies_colors, align='center')
    ax.set_xticks(np.arange(len(hist)))
    ax.set_xticklabels(feature_names)
    ax.set_yticks([0, 1])
    ax.set_yticklabels(['miss', 'hit'])
    ax.set_title(w)<jupyter_output><empty_output><jupyter_text>Feature imputationLet us fill missing PE values. This is the only cell that differs from the approach of Paolo Bestagini. Currently no feature engineering is used, but this should be explored in the future.<jupyter_code>def make_pe(X, seed):
    reg = RandomForestRegressor(max_features='sqrt', n_estimators=50, random_state=seed)
    DataImpAll = data[feature_names].copy()
    DataImp = DataImpAll.dropna(axis = 0, inplace=False)
    Ximp=DataImp.loc[:, DataImp.columns != 'PE']
    Yimp=DataImp.loc[:, 'PE']
    reg.fit(Ximp, Yimp)
    X[np.array(DataImpAll.PE.isnull()),4] = reg.predict(DataImpAll.loc[DataImpAll.PE.isnull(),:].drop('PE',axis=1,inplace=False))
    return X<jupyter_output><empty_output><jupyter_text>Feature augmentationOur guess is that facies do not abrutly change from a given depth layer to the next one. Therefore, we consider features at neighboring layers to be somehow correlated. To possibly exploit this fact, let us perform feature augmentation by:- Aggregating features at neighboring depths.- Computing feature spatial gradient.<jupyter_code># Feature windows concatenation function
def augment_features_window(X, N_neig):
    
    # Parameters
    N_row = X.shape[0]
    N_feat = X.shape[1]

    # Zero padding
    X = np.vstack((np.zeros((N_neig, N_feat)), X, (np.zeros((N_neig, N_feat)))))

    # Loop over windows
    X_aug = np.zeros((N_row, N_feat*(2*N_neig+1)))
    for r in np.arange(N_row)+N_neig:
        this_row = []
        for c in np.arange(-N_neig,N_neig+1):
            this_row = np.hstack((this_row, X[r+c]))
        X_aug[r-N_neig] = this_row

    return X_aug
# Feature gradient computation function
def augment_features_gradient(X, depth):
    
    # Compute features gradient
    d_diff = np.diff(depth).reshape((-1, 1))
    d_diff[d_diff==0] = 0.001
    X_diff = np.diff(X, axis=0)
    X_grad = X_diff / d_diff
        
    # Compensate for last missing value
    X_grad = np.concatenate((X_grad, np.zeros((1, X_grad.shape[1]))))
    
    return X_grad
# Feature augmentation function
def augment_features(X, well, depth, seed=None, pe=True, N_neig=1):
    seed = seed or None
    
    if pe:
        X = make_pe(X, seed)
    # Augment features
    X_aug = np.zeros((X.shape[0], X.shape[1]*(N_neig*2+2)))
    for w in np.unique(well):
        w_idx = np.where(well == w)[0]
        X_aug_win = augment_features_window(X[w_idx, :], N_neig)
        X_aug_grad = augment_features_gradient(X[w_idx, :], depth[w_idx])
        X_aug[w_idx, :] = np.concatenate((X_aug_win, X_aug_grad), axis=1)
    
    # Find padded rows
    padded_rows = np.unique(np.where(X_aug[:, 0:7] == np.zeros((1, 7)))[0])
    
    return X_aug, padded_rows
# Augment features
X_aug, padded_rows = augment_features(X, well, depth)<jupyter_output><empty_output><jupyter_text>Generate training, validation and test data splitsThe choice of training and validation data is paramount in order to avoid overfitting and find a solution that generalizes well on new data. For this reason, we generate a set of training-validation splits so that:- Features from each well belongs to training or validation set.- Training and validation sets contain at least one sample for each class.<jupyter_code># Initialize model selection methods
lpgo = LeavePGroupsOut(2)

# Generate splits
split_list = []
for train, val in lpgo.split(X, y, groups=data['Well Name']):
    hist_tr = np.histogram(y[train], bins=np.arange(len(facies_names)+1)+.5)
    hist_val = np.histogram(y[val], bins=np.arange(len(facies_names)+1)+.5)
    if np.all(hist_tr[0] != 0) & np.all(hist_val[0] != 0):
        split_list.append({'train':train, 'val':val})
            
# Print splits
for s, split in enumerate(split_list):
    print('Split %d' % s)
    print('    training:   %s' % (data['Well Name'][split['train']].unique()))
    print('    validation: %s' % (data['Well Name'][split['val']].unique()))<jupyter_output>Split 0
    training:   ['SHRIMPLIN' 'SHANKLE' 'LUKE G U' 'KIMZEY A' 'CROSS H CATTLE' 'NOLAN'
 'Recruit F9' 'NEWBY']
    validation: ['ALEXANDER D' 'CHURCHMAN BIBLE']
Split 1
    training:   ['SHRIMPLIN' 'SHANKLE' 'LUKE G U' 'KIMZEY A' 'NOLAN' 'Recruit F9' 'NEWBY'
 'CHURCHMAN BIBLE']
    validation: ['ALEXANDER D' 'CROSS H CATTLE']
Split 2
    training:   ['SHRIMPLIN' 'SHANKLE' 'LUKE G U' 'CROSS H CATTLE' 'NOLAN' 'Recruit F9'
 'NEWBY' 'CHURCHMAN BIBLE']
    validation: ['ALEXANDER D' 'KIMZEY A']
Split 3
    training:   ['SHRIMPLIN' 'SHANKLE' 'LUKE G U' 'KIMZEY A' 'CROSS H CATTLE' 'Recruit F9'
 'NEWBY' 'CHURCHMAN BIBLE']
    validation: ['ALEXANDER D' 'NOLAN']
Split 4
    training:   ['SHRIMPLIN' 'LU<NAME>' 'KIMZEY A' 'CROSS H CATTLE' 'NOLAN' 'Recruit F9'
 'NEWBY' 'CHURCHMAN BIBLE']
    validation: ['ALEXANDER D' 'SHANKLE']
Split 5
    training:   ['SHRIMPLIN' '<NAME>' 'SHANKLE' 'LUKE G U' 'KIMZEY A' 'NOLAN'
 'Recruit F9' 'NEWBY']
    validation: ['CROSS H CATTLE' 'CHURCHMAN BIBLE'[...]<jupyter_text>Classification parameters optimizationLet us perform the following steps for each set of parameters:- Select a data split.- Normalize features using a robust scaler.- Train the classifier on training data.- Test the trained classifier on validation data.- Repeat for all splits and average the F1 scores.At the end of the loop, we select the classifier that maximizes the average F1 score on the validation set. Hopefully, this classifier should be able to generalize well on new data.<jupyter_code># Parameters search grid (uncomment parameters for full grid search... may take a lot of time)
N_grid = [100]  # [50, 100, 150]
M_grid = [10]  # [5, 10, 15]
S_grid = [25]  # [10, 25, 50, 75]
L_grid = [5] # [2, 3, 4, 5, 10, 25]
param_grid = []
for N in N_grid:
    for M in M_grid:
        for S in S_grid:
            for L in L_grid:
                param_grid.append({'N':N, 'M':M, 'S':S, 'L':L})
# Train and test a classifier
def train_and_test(X_tr, y_tr, X_v, well_v, clf):
    
    # Feature normalization
    scaler = preprocessing.RobustScaler(quantile_range=(25.0, 75.0)).fit(X_tr)
    X_tr = scaler.transform(X_tr)
    X_v = scaler.transform(X_v)
    
    # Train classifier
    clf.fit(X_tr, y_tr)
    
    # Test classifier
    y_v_hat = clf.predict(X_v)
    
    # Clean isolated facies for each well
    for w in np.unique(well_v):
        y_v_hat[well_v==w] = medfilt(y_v_hat[well_v==w], kernel_size=5)
    
    return y_v_hat
# For each set of parameters
# score_param = []
# for param in param_grid:

#     # For each data split
#     score_split = []
#     for split in split_list:

#         # Remove padded rows
#         split_train_no_pad = np.setdiff1d(split['train'], padded_rows)

#         # Select training and validation data from current split
#         X_tr = X_aug[split_train_no_pad, :]
#         X_v = X_aug[split['val'], :]
#         y_tr = y[split_train_no_pad]
#         y_v = y[split['val']]

#         # Select well labels for validation data
#         well_v = well[split['val']]

#         # Train and test
#         y_v_hat = train_and_test(X_tr, y_tr, X_v, well_v, param)

#         # Score
#         score = f1_score(y_v, y_v_hat, average='micro')
#         score_split.append(score)

#     # Average score for this param
#     score_param.append(np.mean(score_split))
#     print('F1 score = %.3f %s' % (score_param[-1], param))

# # Best set of parameters
# best_idx = np.argmax(score_param)
# param_best = param_grid[best_idx]
# score_best = score_param[best_idx]
# print('\nBest F1 score = %.3f %s' % (score_best, param_best))<jupyter_output>F1 score = 0.565 {'S': 25, 'M': 10, 'L': 5, 'N': 100}

Best F1 score = 0.565 {'S': 25, 'M': 10, 'L': 5, 'N': 100}<jupyter_text>Predict labels on test dataLet us now apply the selected classification technique to test data.<jupyter_code>param_best = {'S': 25, 'M': 10, 'L': 5, 'N': 100}
# Load data from file
test_data = pd.read_csv('../validation_data_nofacies.csv')
# Prepare test data
well_ts = test_data['Well Name'].values
depth_ts = test_data['Depth'].values
X_ts = test_data[feature_names].values
y_pred = []
print('o' * 100)
for seed in range(100):
    np.random.seed(seed)

    # Make training data.
    X_train, padded_rows = augment_features(X, well, depth, seed=seed)
    y_train = y
    X_train = np.delete(X_train, padded_rows, axis=0)
    y_train = np.delete(y_train, padded_rows, axis=0) 
    param = param_best
    clf = OneVsOneClassifier(RandomForestClassifier(n_estimators=param['N'], criterion='entropy',
                             max_features=param['M'], min_samples_split=param['S'], min_samples_leaf=param['L'],
                             class_weight='balanced', random_state=seed), n_jobs=-1)
    
    # Make blind data.
    X_test, _ = augment_features(X_ts, well_ts, depth_ts, seed=seed, pe=False)

    # Train and test.
    y_ts_hat = train_and_test(X_train, y_train, X_test, well_ts, clf)
    
    # Collect result.
    y_pred.append(y_ts_hat)
    print('.', end='')
    
np.save('100_realizations.npy', y_pred)<jupyter_output>....................................................................................................
{'fraction_non_alphanumeric': 0.13382084095063984, 'fraction_numerical': 0.01579524680073126, 'mean_word_length': 3.9318427695636493, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 34, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '40147', 'license': 'Apache-2.0', 'n_tokens_mistral': 4944, 'n_tokens_neox': 4739, 'n_words': 1303, 'path': 'ar4/ar4_submission2_VALIDATION.ipynb', 'repo_name': 'Arnoldsynchron/2016-ml-contest', 'size': '13698'}