%matplotlib inline
# to install watermark magic command: pip install ipyext
%load_ext watermark 
%watermark -v -p numpy,scipy,pandas,matplotlib,seaborn,sklearn

"""
Explanation: Facies classification using Random Forest
Contest entry by geoLEARN:
Martin Blouin, Lorenzo Perozzi, Antoine Caté <br>
in collaboration with  Erwan Gloaguen
Original contest notebook by Brendon Hall, Enthought
In this notebook we will train a machine learning algorithm to predict facies from well log data.  The dataset comes from a class exercise from The University of Kansas on Neural Networks and Fuzzy Systems.  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) and Dubois et al. (2007). 
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 Random Forest model to classify facies types.
Original approach

We focus on feature engineering instead of focusing on algorithm optimization.
The prediction approach has been split in two steps. 1) try to predict non-marine facies; 2) try to predict marine facies. 
Validation scoring has been computed after the prediction results of the two approaches have been stacked together.

N.B We are using exhaustively an interpreted variable (NM_M). We are well aware that this may not be representative of a true machine learning experience where you don't have access to interpreted information. But for now, since we have access to it, why not use it. Moreover, we could easily adapt our approach in a 2 steps prediction where we would first discriminate marine from non-marine sediments and then classify each facies. 
Print the versions of  the used packages
End of explanation
"""


import pandas as pd
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
from sklearn import preprocessing

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

# turn off ipython warnings 
import warnings
warnings.filterwarnings('ignore')

filename = '../facies_vectors.csv'
training_data = pd.read_csv(filename)
training_data.describe()

"""
Explanation: Exploring the dataset
First, we will examine the data set we will use to train the classifier.
End of explanation
"""


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

"""
Explanation: 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. Validation (test) data (830 examples from two wells) have 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 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 (GR), resistivity logging (ILD_log10),
photoelectric effect (PE), neutron-density porosity difference and average neutron-density porosity (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 sandstone
2. Nonmarine coarse siltstone 
3. Nonmarine fine siltstone 
4. Marine siltstone and shale 
5. Mudstone (limestone)
6. Wackestone (limestone)
7. Dolomite
8. 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| 2
2 |CSiS| 1,3
3 |FSiS| 2
4 |SiSh| 5
5 |MS| 4,6
6 |WS| 5,7
7 |D| 6,8
8 |PS| 6,7,9
9 |BS| 7,8
Let's clean up this dataset.  The 'Well Name' and 'Formation' columns can be turned into a categorical data type.
End of explanation
"""


# 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]

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

"""
Explanation: 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.
End of explanation
"""


def make_dist_mar_vars(wells_df):
 
    grouped = wells_df.groupby(['Well Name'])
    new_df = pd.DataFrame()
    
    for key in grouped.groups.keys():
    
        NM_M = grouped.get_group(key)['NM_M'].values

        #We create a temporary dataframe that we reset for every well
        temp_df = pd.DataFrame()
        temp_df['Depth'] = grouped.get_group(key)['Depth']
        temp_df['Well Name'] = [key for _ in range(len(NM_M))]

        #We initialize new variables
        dist_mar_up = np.zeros(len(NM_M))
        dist_mar_down = np.zeros(len(NM_M))

        # A variable counting the interval from the upper marine deposit an one for bottom lower deposit
        # We initialize them to -99999 since we do not know what's above the first log
        count = -99999
        
        #we build them in two seperate loops
        for i in range(len(NM_M)):

            if ((NM_M[i] == 1) & (count>-99999)):

                count+=0.5
                dist_mar_up[i] += count

            elif NM_M[i] == 2:

                count=0

            else:
                dist_mar_up[i] = count

        #********************************************#
        #we reset count
        count = -99999
        
        for i in range(len(NM_M)-1,-1,-1):
               
            if ((NM_M[i] == 1) & (count>-99999)):

                count+=0.5
                dist_mar_down[i] += count

            elif NM_M[i] == 2:            
                count=0

            else:
                dist_mar_down[i] = count
                
        #********************************************#       

        temp_df['dist_mar_up'] = dist_mar_up
        temp_df['dist_mar_down'] = dist_mar_down
    
        # We append each well variable to a larger dataframe
        # We use a dataframe to preserve the index
        new_df = new_df.append(temp_df) 
        
    new_df = new_df.sort_index()
    new_df = new_df.drop(['Well Name','Depth'],axis=1)
    
    #We don't use merge as it creates duplicates for curious reasons that we later have to drop
    return pd.concat([wells_df,new_df],axis=1)

training_data = make_dist_mar_vars(training_data)
training_data.head()

"""
Explanation: We build a function to create new variables representing the absolute distance of a deposit to an upper marine context (younger) and one for the distance to a deeper marine context (older)
Sedimentary facies tend to be deposited in sequences that reflect the evolution of the paleo-environment (variations in water depth, water temperature, biological activity, currents strenght, detrital input, ...). Each facies represents a specific depositional environment and is in contact with facies that represent a progressive transition to an other environment. The distance to a marine context introduces information on the depositional environment of samples neighbouring the predicted sample and thus improves the classification.
This new variable illustrates our approach, which favorises the combination of machine learning algorithms with feature engineering adapted to the target variables.
End of explanation
"""


LE = preprocessing.LabelEncoder()
training_data['Formation_category'] = LE.fit_transform(training_data.Formation)

"""
Explanation: Encoding Formation labels
The nature of the 'Formation' variable is not entirely clear, and we do not know how it has been obtained. In this Notebook, we include it in the predictor variables. However, if 'Formation' appears to be in part derived from 'Facies', it should be removed from the predictor variables.
Formation labels are encoded as integer values.
End of explanation
"""


import seaborn as sns
training_data['NM_M_label'] = training_data['NM_M'].apply(lambda x: 'non_marine' if x == 1 else 'marine')

plt.figure(figsize=(8,6))

sns.countplot(x="FaciesLabels", hue="NM_M_label", 
              order=['SS', 'CSiS', 'FSiS', 'SiSh','MS','WS','D','PS','BS'], 
              data=training_data, palette = "RdBu");

"""
Explanation: Let's see how well NM_M variable discriminates facies
End of explanation
"""


training_data.head()

training_data = training_data.drop(['NM_M_label','FaciesLabels'], axis=1)
NM_training_data = training_data[training_data['NM_M'] == 1]
M_training_data = training_data[training_data['NM_M'] == 2]

"""
Explanation: Now, as the NM_M variable discriminate fairly well facies, we assume that facies 1-3 (SS to FSiS) are non marine and facies 4-9 (SiSh to BS) are marine.
End of explanation
"""


NM_training_data.replace(to_replace=np.nan,value=-99999,inplace=True)
M_training_data.replace(to_replace=np.nan,value=-99999,inplace=True)

"""
Explanation: Instead of dropping samples with the NaN value for the PE variable, we decided to replace NaN with an 'out of the range' value (i.e. -99999). This allow the classifier to take into account all the samples despite missing values for the PE variable.
End of explanation
"""


nm_drop_list = ['Formation', 'Well Name', 'Depth', 
                'Facies','DeltaPHI','PE','NM_M','RELPOS']

m_drop_list = ['Formation', 'Well Name', 'Depth', 
               'Facies','dist_mar_down','dist_mar_up','PE','NM_M','PHIND']

from sklearn import ensemble
from sklearn import metrics
clf = ensemble.RandomForestClassifier(n_estimators=500,n_jobs=-1)

names = list(np.unique(training_data['Well Name']))

nm_grouped = NM_training_data.groupby(['Well Name'])
m_grouped = M_training_data.groupby(['Well Name'])

new_df = pd.DataFrame()

scores = []

for name in names:

    temp_df = pd.DataFrame()

    # We need to isolate Recruit F9 since it has only marine facies
    if name == 'Recruit F9':

        #Build list of well names and remove blind well (and remove recruit F9 for NM)
        m_train_names = names.copy()
        m_train_names.remove(name)

        # Do it for marine sediments
        m_test = m_grouped.get_group(name)
        m_test_depth = m_test.Depth
        m_X_test = m_test.drop(m_drop_list, axis=1).values
        y_test = m_test['Facies'].values

        m_train = pd.DataFrame()

        for train_name in m_train_names:
            m_train = m_train.append(m_grouped.get_group(train_name))

        id_train_M = m_train.Facies >= 4
        m_train = m_train [id_train_M]
        m_X_train = m_train.drop(m_drop_list, axis=1).values
        m_y_train = m_train['Facies'].values

        #Then we do random forest classification and prediction
        clf.fit(m_X_train, m_y_train)
        y_pred = clf.predict(m_X_test)

    else:

        #Build list of well names and remove blind well (and remove recruit F9 for NM)
        m_train_names = names.copy()
        m_train_names.remove(name)   

        nm_train_names = m_train_names.copy()
        nm_train_names.remove('Recruit F9')


        # Do it for non-marine sediments
        nm_test = nm_grouped.get_group(name)
        nm_test_depth = nm_test.Depth
        nm_X_test = nm_test.drop(nm_drop_list, axis=1).values
        nm_y_test = nm_test['Facies'].values

        nm_train = pd.DataFrame()

        for train_name in nm_train_names:
            nm_train = nm_train.append(nm_grouped.get_group(train_name))

        id_train_NM = nm_train.Facies <= 3
        nm_train = nm_train [id_train_NM]
        nm_X_train = nm_train.drop(nm_drop_list, axis=1).values
        nm_y_train = nm_train['Facies'].values

        #Then we do random forest classification and prediction
        clf.fit(nm_X_train, nm_y_train)
        nm_y_pred = clf.predict(nm_X_test)

        print(clf.feature_importances_)

        #*********************************************************************#

        # Do it for marine sediments
        m_test = m_grouped.get_group(name)
        m_test_depth = m_test.Depth
        m_X_test = m_test.drop(m_drop_list, axis=1).values
        m_y_test = m_test['Facies'].values

        m_train = pd.DataFrame()

        for train_name in m_train_names:
            m_train = m_train.append(m_grouped.get_group(train_name))

        id_train_M = m_train.Facies >= 4
        m_train = m_train [id_train_M]
        m_X_train = m_train.drop(m_drop_list, axis=1).values
        m_y_train = m_train['Facies'].values

        #Then we do random forest classification and prediction
        clf.fit(m_X_train, m_y_train)
        m_y_pred = clf.predict(m_X_test)

        print(clf.feature_importances_)

        #================================================================#
        #                         combine results
        #================================================================#

        y_test = np.hstack((nm_y_test,m_y_test))
        y_pred = np.hstack((nm_y_pred,m_y_pred))

    #Scoring
    conf_mat = metrics.confusion_matrix(y_test,y_pred)
    print(conf_mat)
    try:
        score = metrics.f1_score(y_test, y_pred,average='weighted')
    except:
        #this exception is for Recruit F9
        score = metrics.f1_score(y_test, y_pred)
    scores.append(score)
    print('********')
    print('Blind well is {0}, F1 score : {1:.4%}\n'.format(name,score))

    if name == 'Recruit F9':
        depth = m_test_depth
    else:
        depth = np.hstack((nm_test_depth,m_test_depth))

    idx = np.argsort(depth)
    temp_df['Depth'] = depth[idx]
    temp_df['True Facies'] = y_test[idx]
    temp_df['Predicted Facies'] = y_pred[idx]

    temp_df['Well Name'] = [name for _ in range(len(depth))]
    new_df = new_df.append(temp_df)

print("="*30)
print('*********** RESULT ***********')
print("="*30)
print('\nAverage  F1-score is {:.4%}'.format(np.mean(scores)))

"""
Explanation: We drop unused variables for marine and non-marine dataset (we prepare the list here)
Non-marine facies are dominated by clastic sedimentary rocks, and marine facies by carbonate sedimentary rocks. The best predictive variables for each group differ. Thus, marine and non-marine rocks are classified separately here, and a different set of predictive variables is chosen for each group.
End of explanation
"""


filename = '../validation_data_nofacies.csv'
test_data = pd.read_csv(filename)
test_data

"""
Explanation: Applying the classification model to new data
End of explanation
"""


# absoulte distance
test_data = make_dist_mar_vars(test_data)
# formation category encoding
test_data['Formation_category'] = LE.fit_transform(test_data.Formation)

"""
Explanation: Create new variables for the test dataset

we compute the absolute distance of a deposit to an upper marine context (younger) and one for the distance to a deeper marine context (older)
we encode the formation category
End of explanation
"""


#  nm = non-marine and m = marine
nm_drop_list = ['Formation', 'Well Name', 'Depth','DeltaPHI','PE',
                'NM_M','RELPOS']
m_drop_list = ['Formation', 'Well Name', 'Depth',
               'dist_mar_down','dist_mar_up','PE','NM_M','PHIND']

"""
Explanation: We train our classifier to the entire available dataset (10 wells), keeping the separation between marine and non-marine sediments
End of explanation
"""


grouped = test_data.groupby('Well Name')
new_df = pd.DataFrame()
    
for name in grouped.groups.keys():
    
    temp_df = pd.DataFrame()
    temp_df['Depth'] = grouped.get_group(name)['Depth']

    nm_test = grouped.get_group(name)[test_data['NM_M'] == 1]
    m_test = grouped.get_group(name)[test_data['NM_M'] == 2]

    # Do it for non-marine sediments
    nm_test_depth = nm_test.Depth
    nm_X_test = nm_test.drop(nm_drop_list, axis=1).values

    id_train_NM = NM_training_data.Facies <= 3


    nm_train = NM_training_data[id_train_NM]

    nm_X_train = nm_train.drop(nm_drop_list, axis=1)
    nm_X_train = nm_X_train.drop('Facies', axis=1).values
    nm_y_train = nm_train['Facies'].values

    #Then we do random forest classification
    clf.fit(nm_X_train, nm_y_train)
    nm_y_pred = clf.predict(nm_X_test)

    #*********************************************************************#

    # Do it for marine sediments
    m_test_depth = m_test.Depth
    m_X_test = m_test.drop(m_drop_list, axis=1).values

    id_train_M = M_training_data.Facies >= 4
    m_train = M_training_data[id_train_M]

    m_X_train = m_train.drop('Facies', axis=1)
    m_X_train = m_X_train.drop(m_drop_list,axis=1).values
    m_y_train = m_train['Facies'].values

    #Then we do random forest classification
    clf.fit(m_X_train, m_y_train)
    m_y_pred = clf.predict(m_X_test)

    #================================================================#
    #                         combine results
    #================================================================#

    y_pred = np.hstack((nm_y_pred,m_y_pred))


    depth = np.hstack((nm_test_depth,m_test_depth))
    idx = np.argsort(depth)
    
    temp_df['Predicted Facies'] = y_pred[idx]
    temp_df['Well Name'] = [name for _ in range(len(depth))]
    
    new_df = new_df.append(temp_df)
new_df = new_df.sort_index()

new_df.to_pickle('Prediction_blind_wells.pkl')
new_df

"""
Explanation: Train and making prediction for the the well test
End of explanation
"""


