<jupyter_start><jupyter_text>Facies classification using plurality voting (e.g. multiclass majority voting) Contest entry by: <NAME> and <NAME>  [Original contest notebook](../Facies_classification.ipynb) by <NAME>, [Enthought](https://www.enthought.com/) The code and ideas in this notebook, by <NAME> and <NAME>, are licensed under a Creative Commons Attribution 4.0 International License. In this notebook we will attempt to predict facies from well log data using machine learnig classifiers.  The dataset comes from a class exercise from The University of Kansas on [Neural Networks and Fuzzy Systems](http://www.people.ku.edu/~gbohling/EECS833/).  This exercise is based on a consortium project to use machine learning techniques to create a reservoir model of the largest gas fields in North America, the Hugoton and Panoma Fields. For more info on the origin of the data, see [Bohling and Dubois (2003)](http://www.kgs.ku.edu/PRS/publication/2003/ofr2003-50.pdf) and [Dubois et al. (2007)](http://dx.doi.org/10.1016/j.cageo.2006.08.011). The dataset consists of log data from nine wells that have been labeled with a facies type based on observation of core. We will use this log data to train a support vector machine to classify facies types.   The planWe will created three classifiers with pretuned parameters:- best SVM in the competition (by our team's SVM submission) - best Random Forest in the competition (form the leading submission, by gccrowther) - multilayer perceptron (from previous notebooks, not submitted)We will then try to predict the facies using a plurality voting approach (plurality voting = multi-class majority voting).From the [scikit-learn website](http://scikit-learn.org/stable/modules/ensemble.htmlvoting-classifier): "The idea behind the voting classifier implementation is to combine conceptually different machine learning classifiers and use a majority vote or the average predicted probabilities (soft vote) to predict the class labels. Such a classifier can be useful for a set of equally well performing model in order to balance out their individual weaknesses".  Exploring the datasetFirst, we will examine the data set we will use to train the classifier.<jupyter_code>%matplotlib inline
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.colors as colors
from mpl_toolkits.axes_grid1 import make_axes_locatable

import pandas as pd
from pandas import set_option
set_option("display.max_rows", 10)
pd.options.mode.chained_assignment = None

from sklearn import preprocessing
from sklearn.metrics import f1_score, accuracy_score, make_scorer
from sklearn.model_selection import LeaveOneGroupOut
filename = 'facies_vectors.csv'
training_data = pd.read_csv(filename)
training_data<jupyter_output><empty_output><jupyter_text>This data is from the Council Grove gas reservoir in Southwest Kansas.  The Panoma Council Grove Field is predominantly a carbonate gas reservoir encompassing 2700 square miles in Southwestern Kansas.  This dataset is from nine wells (with 4149 examples), consisting of a set of seven predictor variables and a rock facies (class) for each example vector and validation (test) data (830 examples from two wells) having the same seven predictor variables in the feature vector.  Facies are based on examination of cores from nine wells taken vertically at half-foot intervals. Predictor variables include five from wireline log measurements and two geologic constraining variables that are derived from geologic knowledge. These are essentially continuous variables sampled at a half-foot sample rate. The seven predictor variables are:* Five wire line log curves include [gamma ray](http://petrowiki.org/Gamma_ray_logs) (GR), [resistivity logging](http://petrowiki.org/Resistivity_and_spontaneous_%28SP%29_logging) (ILD_log10),[photoelectric effect](http://www.glossary.oilfield.slb.com/en/Terms/p/photoelectric_effect.aspx) (PE), [neutron-density porosity difference and average neutron-density porosity](http://petrowiki.org/Neutron_porosity_logs) (DeltaPHI and PHIND). Note, some wells do not have PE.* Two geologic constraining variables: nonmarine-marine indicator (NM_M) and relative position (RELPOS)The nine discrete facies (classes of rocks) are: 1. Nonmarine sandstone2. Nonmarine coarse siltstone 3. Nonmarine fine siltstone 4. Marine siltstone and shale 5. Mudstone (limestone)6. Wackestone (limestone)7. Dolomite8. Packstone-grainstone (limestone)9. Phylloid-algal bafflestone (limestone)These facies aren't discrete, and gradually blend into one another. Some have neighboring facies that are rather close.  Mislabeling within these neighboring facies can be expected to occur.  The following table lists the facies, their abbreviated labels and their approximate neighbors.Facies |Label| Adjacent Facies:---: | :---: |:--:1 |SS| 22 |CSiS| 1,33 |FSiS| 24 |SiSh| 55 |MS| 4,66 |WS| 5,77 |D| 6,88 |PS| 6,7,99 |BS| 7,8Let's clean up this dataset.  The 'Well Name' and 'Formation' columns can be turned into a categorical data type.<jupyter_code>training_data['Well Name'] = training_data['Well Name'].astype('category')
training_data['Formation'] = training_data['Formation'].astype('category')
training_data['Well Name'].unique()<jupyter_output><empty_output><jupyter_text>These are the names of the 10 training wells in the Council Grove reservoir.  Data has been recruited into pseudo-well 'Recruit F9' to better represent facies 9, the Phylloid-algal bafflestone. Before we plot the well data, let's define a color map so the facies are represented by consistent color in all the plots in this tutorial.  We also create the abbreviated facies labels, and add those to the `facies_vectors` dataframe.<jupyter_code># 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','#A569BD',
       '#000000', '#000080', '#2E86C1', '#AED6F1', '#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]

def label_facies(row, labels):
    return labels[ row['Facies'] -1]
    
training_data.loc[:,'FaciesLabels'] = training_data.apply(lambda row: label_facies(row, facies_labels), axis=1)
training_data.describe()<jupyter_output><empty_output><jupyter_text>This is a quick view of the statistical distribution of the input variables.  Looking at the `count` values, most values have 4149 valid values except for `PE`, which has 3232.  We will drop the feature vectors that don't have a valid `PE` entry.<jupyter_code>PE_mask = training_data['PE'].notnull().values
training_data = training_data[PE_mask]
training_data.describe()<jupyter_output><empty_output><jupyter_text>Now we extract just the feature variables we need to perform the classification.  The predictor variables are the five log values and two geologic constraining variables, **and we are also using depth**. We also get a vector of the facies labels that correspond to each feature vector.<jupyter_code>y = training_data['Facies'].values
print y[25:40]
print np.shape(y)
X = training_data.drop(['Formation', 'Well Name','Facies','FaciesLabels'], axis=1)
print np.shape(X)
X.describe(percentiles=[.05, .25, .50, .75, .95])
scaler = preprocessing.StandardScaler().fit(X)
X = scaler.transform(X)<jupyter_output><empty_output><jupyter_text>Make performance scorersUsed to evaluate performance.<jupyter_code>Fscorer = make_scorer(f1_score, average = 'micro')<jupyter_output><empty_output><jupyter_text>Pre-tuned SVM classifier classifier and leave one well out average F1 scoreThis is the Support Vector Machine classifier from our [first submission](https://github.com/mycarta/2016-ml-contest/blob/master/MandMs/Facies_classification-M%26Ms_SVM_rbf_kernel_optimal.ipynb).<jupyter_code>from sklearn import svm
SVC_classifier = svm.SVC(C = 100, cache_size=2400, class_weight=None, coef0=0.0,
  decision_function_shape=None, degree=3, gamma=0.01, kernel='rbf',
  max_iter=-1, probability=True, random_state=49, shrinking=True,
  tol=0.001, verbose=False)
f1_svc = []

wells = training_data["Well Name"].values
logo = LeaveOneGroupOut()

for train, test in logo.split(X, y, groups=wells):
    well_name = wells[test[0]]
    SVC_classifier.fit(X[train], y[train])
    pred_svc = SVC_classifier.predict(X[test])
    sc = f1_score(y[test], pred_svc, labels = np.arange(10), average = 'micro')
    print("{:>20s}  {:.3f}".format(well_name, sc))
    f1_svc.append(sc)
    
print "-Average leave-one-well-out F1 Score: %6f" % (sum(f1_svc)/(1.0*(len(f1_svc))))<jupyter_output>CHURCHMAN BIBLE  0.542
      CROSS H CATTLE  0.347
            LUKE G U  0.440
               NEWBY  0.400
               NOLAN  0.494
          Recruit F9  0.721
             SHANKLE  0.483
           SHRIMPLIN  0.590
-Average leave-one-well-out F1 Score: 0.502174<jupyter_text>Pre-tuned multi-layer perceptron classifier and average F1 score<jupyter_code>from sklearn.neural_network import MLPClassifier
mlp_classifier = MLPClassifier(activation='logistic', alpha=0.01, batch_size='auto',
       beta_1=0.9, beta_2=0.999, early_stopping=False, epsilon=1e-08,
       hidden_layer_sizes=(100,), learning_rate='adaptive',
       learning_rate_init=0.001, max_iter=1000, momentum=0.9,
       nesterovs_momentum=True, power_t=0.5, random_state=49, shuffle=True,
       solver='adam', tol=0.0001, validation_fraction=0.1, verbose=False,
       warm_start=False)
f1_mlp = []

wells = training_data["Well Name"].values
logo = LeaveOneGroupOut()

for train, test in logo.split(X, y, groups=wells):
    well_name = wells[test[0]]
    mlp_classifier.fit(X[train], y[train])
    pred_mlp = mlp_classifier.predict(X[test])
    sc = f1_score(y[test], pred_mlp, labels = np.arange(10), average = 'micro')
    print("{:>20s}  {:.3f}".format(well_name, sc))
    f1_mlp.append(sc)
    
print "-Average leave-one-well-out F1 Score: %6f" % (sum(f1_mlp)/(1.0*(len(f1_mlp))))<jupyter_output>CHURCHMAN BIBLE  0.525
      CROSS H CATTLE  0.341
            LUKE G U  0.419
               NEWBY  0.415
               NOLAN  0.482
          Recruit F9  0.779
             SHANKLE  0.541
           SHRIMPLIN  0.575
-Average leave-one-well-out F1 Score: 0.509666<jupyter_text>Pre-tuned extra treesThis is the RF classifier with parameters tuned in the leading submission, by George Crowther, but without his engineered features.<jupyter_code>from sklearn.pipeline import make_pipeline
from sklearn.feature_selection import VarianceThreshold
from sklearn.ensemble import ExtraTreesClassifier

ET_classifier = make_pipeline(
    VarianceThreshold(threshold=0.49),
    ExtraTreesClassifier(criterion="entropy", max_features=0.71,
                         n_estimators=500, random_state=49))
f1_ET = []

wells = training_data["Well Name"].values
logo = LeaveOneGroupOut()

for train, test in logo.split(X, y, groups=wells):
    well_name = wells[test[0]]
    ET_classifier.fit(X[train], y[train])
    pred_cv = ET_classifier.predict(X[test])
    sc = f1_score(y[test], pred_cv, labels = np.arange(10), average = 'micro')
    print("{:>20s}  {:.3f}".format(well_name, sc))
    f1_ET.append(sc)
    
print "-Average leave-one-well-out F1 Score: %6f" % (sum(f1_ET)/(1.0*(len(f1_ET))))<jupyter_output>CH<NAME>  0.498
      CROSS H CATTLE  0.337
            LUKE G U  0.434
               NEWBY  0.408
               NOLAN  0.494
          Recruit F9  0.912
             SHANKLE  0.486
           SHRIMPLIN  0.614
-Average leave-one-well-out F1 Score: 0.522719<jupyter_text>Plurality voting classifier (multi-class majority voting)We will use a weighted approach, where the weights are somewhat arbitrary, but their proportion is based on the average f1 score of the individual classifiers.<jupyter_code>from sklearn.ensemble import VotingClassifier
eclf_cv = VotingClassifier(estimators=[
        ('SVC', SVC_classifier), ('MLP', mlp_classifier), ('ET', ET_classifier)], 
                        voting='soft', weights=[0.3,0.33,0.37])<jupyter_output><empty_output><jupyter_text>Leave one-well-out F1 scores<jupyter_code>f1_ens = []

wells = training_data["Well Name"].values
logo = LeaveOneGroupOut()

for train, test in logo.split(X, y, groups=wells):
    well_name = wells[test[0]]
    eclf_cv.fit(X[train], y[train])
    pred_cv = eclf_cv.predict(X[test])
    sc = f1_score(y[test], pred_cv, labels = np.arange(10), average = 'micro')
    print("{:>20s}  {:.3f}".format(well_name, sc))
    f1_ens.append(sc)
    
print "-Average leave-one-well-out F1 Score: %6f" % (sum(f1_ens)/(1.0*(len(f1_ens))))<jupyter_output><NAME>  0.554
      CROSS H CATTLE  0.351
            LUKE G U  0.451
               NEWBY  0.400
               NOLAN  0.501
          Recruit F9  0.912
             SHANKLE  0.519
           SHRIMPLIN  0.603
-Average leave-one-well-out F1 Score: 0.536423<jupyter_text>Comments Using the average F1 score from the leave-one-well out cross validation as a metric, the majority voting is superior to the individual classifiers, including the pre-tuned Random Forest from the leading submission. However, the Random Forest in the official leading submission was trained using additional new features engineered by George, and outperforms our majority voting classifier, with an F1 score of 0.580 against our 0.579. A clear indication, in our view, that the feature engineering is a key element to achieve the best possible prediction.  Predicting, displaying, and saving facies for blind wells<jupyter_code>blind = pd.read_csv('validation_data_nofacies.csv') 
X_blind = np.array(blind.drop(['Formation', 'Well Name'], axis=1)) 
X_blind = scaler.transform(X_blind) 
y_pred = eclf_cv.fit(X, y).predict(X_blind) 
blind['Facies'] = y_pred
def make_facies_log_plot(logs, facies_colors):
    #make sure logs are sorted by depth
    logs = logs.sort_values(by='Depth')
    cmap_facies = colors.ListedColormap(
            facies_colors[0:len(facies_colors)], 'indexed')
    
    ztop=logs.Depth.min(); zbot=logs.Depth.max()
    
    cluster=np.repeat(np.expand_dims(logs['Facies'].values,1), 100, 1)
    
    f, ax = plt.subplots(nrows=1, ncols=6, figsize=(8, 12))
    ax[0].plot(logs.GR, logs.Depth, '-g')
    ax[1].plot(logs.ILD_log10, logs.Depth, '-')
    ax[2].plot(logs.DeltaPHI, logs.Depth, '-', color='0.5')
    ax[3].plot(logs.PHIND, logs.Depth, '-', color='r')
    ax[4].plot(logs.PE, logs.Depth, '-', color='black')
    im=ax[5].imshow(cluster, interpolation='none', aspect='auto',
                    cmap=cmap_facies,vmin=1,vmax=9)
    
    divider = make_axes_locatable(ax[5])
    cax = divider.append_axes("right", size="20%", pad=0.05)
    cbar=plt.colorbar(im, cax=cax)
    cbar.set_label((17*' ').join([' SS ', 'CSiS', 'FSiS', 
                                'SiSh', ' MS ', ' WS ', ' D  ', 
                                ' PS ', ' BS ']))
    cbar.set_ticks(range(0,1)); cbar.set_ticklabels('')
    
    for i in range(len(ax)-1):
        ax[i].set_ylim(ztop,zbot)
        ax[i].invert_yaxis()
        ax[i].grid()
        ax[i].locator_params(axis='x', nbins=3)
    
    ax[0].set_xlabel("GR")
    ax[0].set_xlim(logs.GR.min(),logs.GR.max())
    ax[1].set_xlabel("ILD_log10")
    ax[1].set_xlim(logs.ILD_log10.min(),logs.ILD_log10.max())
    ax[2].set_xlabel("DeltaPHI")
    ax[2].set_xlim(logs.DeltaPHI.min(),logs.DeltaPHI.max())
    ax[3].set_xlabel("PHIND")
    ax[3].set_xlim(logs.PHIND.min(),logs.PHIND.max())
    ax[4].set_xlabel("PE")
    ax[4].set_xlim(logs.PE.min(),logs.PE.max())
    ax[5].set_xlabel('Facies')
    
    ax[1].set_yticklabels([]); ax[2].set_yticklabels([]); ax[3].set_yticklabels([])
    ax[4].set_yticklabels([]); ax[5].set_yticklabels([])
    ax[5].set_xticklabels([])
    f.suptitle('Well: %s'%logs.iloc[0]['Well Name'], fontsize=14,y=0.94)
make_facies_log_plot(blind[blind['Well Name'] == 'STUART'], facies_colors)
make_facies_log_plot(blind[blind['Well Name'] == 'CRAWFORD'], facies_colors)
np.save('ypred.npy', y_pred)<jupyter_output><empty_output><jupyter_text>Displaying predicted versus original facies in the training dataThis is a nice display to finish up with, as it gives us a visual idea of the predicted faces where we have facies from the core observations.The plot we will use a function from the original notebook. Let's look at the well with the lowest F1 from the previous code block, CROSS H CATTLE, and the one with the highest F1 (excluding Recruit F9), which is SHRIMPLIN.<jupyter_code>def compare_facies_plot(logs, compadre, facies_colors):
    #make sure logs are sorted by depth
    logs = logs.sort_values(by='Depth')
    cmap_facies = colors.ListedColormap(
            facies_colors[0:len(facies_colors)], 'indexed')
    
    ztop=logs.Depth.min(); zbot=logs.Depth.max()
    
    cluster1 = np.repeat(np.expand_dims(logs['Facies'].values,1), 100, 1)
    cluster2 = np.repeat(np.expand_dims(logs[compadre].values,1), 100, 1)
    
    f, ax = plt.subplots(nrows=1, ncols=7, figsize=(9, 12))
    ax[0].plot(logs.GR, logs.Depth, '-g')
    ax[1].plot(logs.ILD_log10, logs.Depth, '-')
    ax[2].plot(logs.DeltaPHI, logs.Depth, '-', color='0.5')
    ax[3].plot(logs.PHIND, logs.Depth, '-', color='r')
    ax[4].plot(logs.PE, logs.Depth, '-', color='black')
    im1 = ax[5].imshow(cluster1, interpolation='none', aspect='auto',
                    cmap=cmap_facies,vmin=1,vmax=9)
    im2 = ax[6].imshow(cluster2, interpolation='none', aspect='auto',
                    cmap=cmap_facies,vmin=1,vmax=9)
    
    divider = make_axes_locatable(ax[6])
    cax = divider.append_axes("right", size="20%", pad=0.05)
    cbar=plt.colorbar(im2, cax=cax)
    cbar.set_label((17*' ').join([' SS ', 'CSiS', 'FSiS', 
                                'SiSh', ' MS ', ' WS ', ' D  ', 
                                ' PS ', ' BS ']))
    cbar.set_ticks(range(0,1)); cbar.set_ticklabels('')
    
    for i in range(len(ax)-2):
        ax[i].set_ylim(ztop,zbot)
        ax[i].invert_yaxis()
        ax[i].grid()
        ax[i].locator_params(axis='x', nbins=3)
    
    ax[0].set_xlabel("GR")
    ax[0].set_xlim(logs.GR.min(),logs.GR.max())
    ax[1].set_xlabel("ILD_log10")
    ax[1].set_xlim(logs.ILD_log10.min(),logs.ILD_log10.max())
    ax[2].set_xlabel("DeltaPHI")
    ax[2].set_xlim(logs.DeltaPHI.min(),logs.DeltaPHI.max())
    ax[3].set_xlabel("PHIND")
    ax[3].set_xlim(logs.PHIND.min(),logs.PHIND.max())
    ax[4].set_xlabel("PE")
    ax[4].set_xlim(logs.PE.min(),logs.PE.max())
    ax[5].set_xlabel('Facies')
    ax[6].set_xlabel(compadre)
    
    ax[1].set_yticklabels([]); ax[2].set_yticklabels([]); ax[3].set_yticklabels([])
    ax[4].set_yticklabels([]); ax[5].set_yticklabels([])
    ax[5].set_xticklabels([])
    ax[6].set_xticklabels([])
    f.suptitle('Well: %s'%logs.iloc[0]['Well Name'], fontsize=14,y=0.94)
eclf_cv.fit(X,y)
pred = eclf_cv.predict(X)

X = training_data
X['Prediction'] = pred
compare_facies_plot(X[X['Well Name'] == 'CROSS H CATTLE'], 'Prediction', facies_colors)
compare_facies_plot(X[X['Well Name'] == 'SHRIMPLIN'], 'Prediction', facies_colors)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1254739028823682, 'fraction_numerical': 0.03318618540638795, 'mean_word_length': 4.262640065591691, 'pattern_counts': {'":': 0, '<': 56, '<?xml version=': 0, '>': 60, 'https://': 2, 'lorem ipsum': 0, 'www.': 4, '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': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '350442', 'license': 'Apache-2.0', 'n_tokens_mistral': 6961, 'n_tokens_neox': 6492, 'n_words': 1703, 'path': 'MandMs/02_Facies_classification-MandMs_plurality_voting_classifier.ipynb', 'repo_name': 'Arnoldsynchron/2016-ml-contest', 'size': '19299'}