<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Facies classification using Machine Learning
LA Team Submission 2 ##
Lukas Mosser, Alfredo De la Fuente
In this python notebook we explore a facies classification model using a Deep Neural Network arquitecture that takes into account spatial dependencies to outperform the prediction model proposed in the prediction facies from wel logs challenge. 
Problem Modeling
The dataset we will use comes from a class excercise 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 we will use is log data from nine wells that have been labeled with a facies type based on oberservation of core.  We will use this log data to train a classifier to predict facies types. 
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
Step1:  We load the training and testing data to preprocess it for further analysis.
Step2:  We fill the missing data values in the PE field with zero and proceed to normalize the data that will be fed into our model.
Step3:  In order to start training stage, it is required to format the data by considering a sliding window over the depth component in order to classify a given set of features at some specific depth for each well in the training set.
Step4:  Data Analysis
Our model will be composed by a Deep Neural Network of an input layer, two hidden layers and an output layer.
As input we will consider a set of features extracted by padding a depth interval segment, that way we take into account spatial dependencies.
Step5:  In order to evaluate our classification model accurary we will use the our following defined metrics, based on the confusion matrix once the classification is performed. The first metric only considers misclassification error and the second one takes into account the fact that facies could be misclassified if they belong to a same group with similar geological characteristics.
Step6:  Once the set of parameters are fixed, the training stage of our model begins. We perform a Cross Validation routine to evaluate the performance of the model.
Step7:  Prediction
We obtain the predictions for test data.
<ASSISTANT_TASK:>
Python Code:
%%sh
pip install pandas
pip install scikit-learn
pip install keras
from __future__ import print_function
import numpy as np
%matplotlib inline
import pandas as pd
import matplotlib.pyplot as plt
from keras.preprocessing import sequence
from keras.models import Sequential
from keras.layers import Dense, Dropout, Activation
from keras.wrappers.scikit_learn import KerasClassifier
from keras.utils import np_utils
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import KFold
from classification_utilities import display_cm, display_adj_cm
from sklearn.metrics import confusion_matrix, f1_score
from sklearn import preprocessing
Explanation: Facies classification using Machine Learning
LA Team Submission 2 ##
Lukas Mosser, Alfredo De la Fuente
In this python notebook we explore a facies classification model using a Deep Neural Network arquitecture that takes into account spatial dependencies to outperform the prediction model proposed in the prediction facies from wel logs challenge. 
Problem Modeling
The dataset we will use comes from a class excercise 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 we will use is log data from nine wells that have been labeled with a facies type based on oberservation of core.  We will use this log data to train a classifier to predict facies types. 
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 (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
Data Preprocessing
Let's import all the libraries that will be particularly needed for the analysis.
End of explanation
filename = 'train_test_data.csv'
training_data = pd.read_csv(filename)
training_data.head(10)
Explanation: We load the training and testing data to preprocess it for further analysis.
End of explanation
# Set 'Well Name' and 'Formation' fields as categories
training_data['Well Name'] = training_data['Well Name'].astype('category')
training_data['Formation'] = training_data['Formation'].astype('category')
# Fill missing values and normalize for 'PE' field
training_data['PE'] = training_data['PE'].fillna(value=0)
mean_pe = training_data['PE'].mean()
std_pe = training_data['PE'].std()
training_data['PE'] = (training_data['PE']-mean_pe)/std_pe
# Normalize the rest of fields (GR, ILD_log10, DelthaPHI, PHIND,NM_M,RELPOS)
correct_facies_labels = training_data['Facies'].values
feature_vectors = training_data.drop(['Formation', 'Depth'], axis=1)
well_labels = training_data[['Well Name', 'Facies']].values
data_vectors = feature_vectors.drop(['Well Name', 'Facies'], axis=1).values
scaler = preprocessing.StandardScaler().fit(data_vectors)
scaled_features = scaler.transform(data_vectors)
data_out = np.hstack([well_labels, scaled_features])
Explanation: We fill the missing data values in the PE field with zero and proceed to normalize the data that will be fed into our model.
End of explanation
def preprocess(data_out):
    data = data_out
    well_data = {}
    well_names = list(set(data[:, 0]))
    for name in well_names:
        well_data[name] = [[], []]
    for row in data:
        well_data[row[0]][1].append(row[1])
        well_data[row[0]][0].append(list(row[2::]))
    # Sliding window
    positive_lag = 5
    negative_lag = 5
    chunks = []
    chunks_test = []
    chunk_length = positive_lag+negative_lag+1 
    chunks_facies = []
    for name in well_names:
        if name not in ['STUART', 'CRAWFORD']:
            test_well_data = well_data[name]
            log_values = np.array(test_well_data[0])
            log_values_padded = np.lib.pad(log_values, (negative_lag,positive_lag), 'edge')[:, negative_lag:-positive_lag]
            facies_values =  np.array(test_well_data[1])
            for i in range(log_values.shape[0]):
                chunks.append(log_values_padded[i:i+chunk_length, :])
                chunks_facies.append(facies_values[i])
        else:
            test_well_data = well_data[name]
            log_values = np.array(test_well_data[0])
            log_values_padded = np.lib.pad(log_values, (negative_lag,positive_lag), 'edge')[:, negative_lag:-positive_lag]
            facies_values =  np.array(test_well_data[1])
            for i in range(log_values.shape[0]):
                chunks_test.append(log_values_padded[i:i+chunk_length, :])
    
    chunks_facies = np.array(chunks_facies, dtype=np.int32)-1
    X_ = np.array(chunks)
    X = np.zeros((len(X_),len(X_[0][0]) * len(X_[0])))
    for i in range(len(X_)):
        X[i,:] = X_[i].flatten()
        
    X_test = np.array(chunks_test)
    X_test_out = np.zeros((len(X_test),len(X_test[0][0]) * len(X_test[0])))
    for i in range(len(X_test)):
        X_test_out[i,:] = X_test[i].flatten()
    y = np_utils.to_categorical(chunks_facies)
    return X, y, X_test_out
Explanation: In order to start training stage, it is required to format the data by considering a sliding window over the depth component in order to classify a given set of features at some specific depth for each well in the training set.
End of explanation
np.random.seed(1337) # For reproducibility
X_train, y_train, X_test = preprocess(data_out)
# Set parameters
input_dim = 77
hidden_dim_1 = 128
hidden_dim_2 = 32
output_dim = 9 
batch_size = 32
nb_epoch = 10
def proposed_model():
    # Define the model
    model = Sequential()
    model.add(Dense(128, input_dim=77, init='normal', activation='relu'))
    model.add(Dense(32, input_dim=128, init='normal', activation='relu'))
    model.add(Dropout(0.3))
    model.add(Dense(9, init='normal', activation='softmax'))
    # Compile model
    model.compile(loss='categorical_crossentropy', optimizer='adam', metrics=['accuracy'])
    return model
Explanation: Data Analysis
Our model will be composed by a Deep Neural Network of an input layer, two hidden layers and an output layer.
As input we will consider a set of features extracted by padding a depth interval segment, that way we take into account spatial dependencies.
End of explanation
def accuracy(conf):
    total_correct = 0.
    nb_classes = conf.shape[0]
    for i in np.arange(0,nb_classes):
        total_correct += conf[i][i]
    acc = total_correct/sum(sum(conf))
    return acc
adjacent_facies = np.array([[1], [0,2], [1], [4], [3,5], [4,6,7], [5,7], [5,6,8], [6,7]])
facies_labels = ['SS', 'CSiS', 'FSiS', 'SiSh', 'MS','WS', 'D','PS', 'BS']
def accuracy_adjacent(conf, adjacent_facies):
    nb_classes = conf.shape[0]
    total_correct = 0.
    for i in np.arange(0,nb_classes):
        total_correct += conf[i][i]
        for j in adjacent_facies[i]:
            total_correct += conf[i][j]
    return total_correct / sum(sum(conf))
Explanation: In order to evaluate our classification model accurary we will use the our following defined metrics, based on the confusion matrix once the classification is performed. The first metric only considers misclassification error and the second one takes into account the fact that facies could be misclassified if they belong to a same group with similar geological characteristics.
End of explanation
estimator = KerasClassifier(build_fn=proposed_model, nb_epoch=10, batch_size=50, verbose=0)
# Cross Validation
kfold = KFold(n_splits=5, shuffle=True)
results = cross_val_score(estimator, X_train, y_train, cv=kfold)
print(' Cross Validation Results')
print( results )
plt.hist(results, bins= 5) 
plt.title("Model Accuracy")
plt.show()
print ( "Cross Validation Average Accuracy %.3f" %  results.mean() )
print ( "Cross Validation Standard Deviation %.3f" %  results.std() )
# Load the model
model = proposed_model()
#Train model
model.fit(X_train, y_train, nb_epoch=10, verbose=2, shuffle = True)
# Predict Values on Training set
y_predicted = model.predict( X_train , batch_size=32, verbose=0)
# Print Report
# Format output [0 - 8 ]
y_ = np.zeros((len(y_train),1))
for i in range(len(y_train)):
    y_[i] = np.argmax(y_train[i])
y_predicted_ = np.zeros((len(y_predicted), 1))
for i in range(len(y_predicted)):
    y_predicted_[i] = np.argmax( y_predicted[i])
    
# Confusion Matrix
conf = confusion_matrix(y_, y_predicted_)
# Print Results
print ("\nModel Report")
print ("-Accuracy: %.6f" % ( accuracy(conf) ))
print ("-Adjacent Accuracy: %.6f" % ( accuracy_adjacent(conf, adjacent_facies) ))
print ("\nConfusion Matrix")
display_cm(conf, facies_labels, display_metrics=True, hide_zeros=True)
Explanation: Once the set of parameters are fixed, the training stage of our model begins. We perform a Cross Validation routine to evaluate the performance of the model.
End of explanation
y_test = model.predict( X_test , batch_size=32, verbose=0)
final_predictions = np.zeros((len(y_test),1))
for i in range(len(y_test)):
    final_predictions[i] = np.argmax(y_test[i]) + 1 
last_submission = pd.read_csv('Prediction.csv')
ls = last_submission['Facies'].as_matrix()
conf = confusion_matrix(ls, final_predictions)
# Print Results
print ("\nModel Report")
print ("-Accuracy: %.6f" % ( accuracy(conf) ))
print ("-Adjacent Accuracy: %.6f" % ( accuracy_adjacent(conf, adjacent_facies) ))
print ("\nConfusion Matrix")
display_cm(conf, facies_labels, display_metrics=True, hide_zeros=True)
# Store results
test_data = pd.read_csv('../validation_data_nofacies.csv')
test_data['Facies'] = final_predictions
test_data.to_csv('Prediction2.csv')
Explanation: Prediction
We obtain the predictions for test data.
End of explanation
<END_TASK>