predicting_seismic_bumps.py

# -*- coding: utf-8 -*-
"""Predicting_Seismic_Bumps.ipynb
Automatically generated by Colaboratory.
Original file is located at
    https://colab.research.google.com/drive/1fIvMom1iQUPN7K_ODtnq9Kb41ZfKH_xK
Nabanita Roy
[Predicting Hazardous Seismic Bumps Part I : EDA, Feature Engineering & Train Test Split for Unbalanced Datase](https://towardsdatascience.com/predicting-hazardous-seismic-bumps-using-supervised-classification-algorithms-part-i-2c5d21f379bc)
[Predicting Hazardous Seismic Bumps Part II: Training & Tuning Supervised ML Classifiers and Model Performance Analysis](https://towardsdatascience.com/predicting-hazardous-seismic-bumps-part-ii-training-supervised-classifier-models-and-8b9104b611b0)
"""



"""## Imports"""

import numpy as np
import pandas as pd

import scipy
from scipy.io import arff
from scipy import stats
import joblib

# scipy.optimize

pd.set_option('display.max_rows', 500)
pd.set_option('display.max_columns', 500)

import matplotlib.pyplot as plt
import seaborn as sns

import time

from sklearn.preprocessing import OneHotEncoder, LabelEncoder
from sklearn.model_selection import KFold, cross_val_score, GridSearchCV, train_test_split, StratifiedKFold, cross_validate, StratifiedShuffleSplit
from sklearn.neighbors import KNeighborsClassifier
from sklearn.metrics import confusion_matrix, precision_score, recall_score, f1_score, accuracy_score
from sklearn.svm import SVC

from collections import Counter
from imblearn.over_sampling import SMOTE # doctest: +NORMALIZE_WHITESPACE

from xgboost import XGBClassifier

"""# Mount Google drive"""

from google.colab import drive
drive.mount('/content/drive')

"""## Load Data"""

!wget https://archive.ics.uci.edu/ml/machine-learning-databases/00266/seismic-bumps.arff

data = arff.loadarff('seismic-bumps.arff')

df = pd.DataFrame(data[0])

df.shape

df.head()

df.columns

"""## Data Preprocessing
### Categorizing columns in types of features and labels
"""

label = 'class'
col_list_categorical = ['seismic', 'seismoacoustic', 'shift', 'ghazard']
col_list_numerical = ['genergy', 'gpuls', 'gdenergy', 'gdpuls', 'energy', 'maxenergy']
col_list_discrete = ['nbumps', 'nbumps2', 'nbumps3', 'nbumps4', 'nbumps5', 'nbumps6', 'nbumps7', 'nbumps89']
len([label]) + len(col_list_categorical) + len(col_list_discrete) + len(col_list_numerical)

"""### Decoding bytes to utf-8"""

for col in df.columns:
    if col in col_list_categorical + [label]:
        df[col] = df[col].str.decode("utf-8")
        if col == label:
            print(df[col].value_counts())
            df[col] = df[col].astype(int)

df.dtypes

df.head()

df.info(verbose=True, null_counts=True)

df.describe()

df.columns

# Dataframe Checkpoint
df_raw = df.copy(deep=True)

"""### Target Distributions Check"""

df['class'].value_counts() * 100 /len(df['class'])

colors = sns.color_palette([sns.color_palette("gist_earth_r")[2], sns.color_palette("gist_earth_r")[5]])
# The colors are:
sns.palplot(sns.color_palette("gist_earth_r")[2]), sns.palplot(sns.color_palette("gist_earth_r")[5])

sns.countplot(x=label, data=df, palette=colors)
plt.xlabel('CLASS')
plt.ylabel('COUNT')

"""### Categorical Feature Encoding
#### Categorical Variables Distributions
[Attribute information](https://archive.ics.uci.edu/ml/datasets/seismic-bumps):
1. seismic: result of shift seismic hazard assessment in the mine working obtained by the seismic
method (a - lack of hazard, b - low hazard, c - high hazard, d - danger state);
2. seismoacoustic: result of shift seismic hazard assessment in the mine working obtained by the
seismoacoustic method;
3. shift: information about type of a shift (W - coal-getting, N -preparation shift);
"""

for col in col_list_categorical:
    print(df[col].value_counts())

"""#### Contingency tables for categorical variables"""

df.seismic.value_counts()

data_crosstab = pd.crosstab(df['seismic'], df[label], colnames=['class'])
data_crosstab

"""Larger fraction in (b | 1)"""

data_crosstab = pd.crosstab(df['seismoacoustic'], df[label], colnames=['class'])
data_crosstab

101/1479 , 66/890

"""Not much bias"""

data_crosstab = pd.crosstab(df['shift'], df[label])
data_crosstab

17/904 , 153/1510

"""10% of the (W | 0) is 1."""

data_crosstab = pd.crosstab(df['ghazard'], df[label])
data_crosstab

2186/156, 198/14

"""Similar ratios"""



"""#### Encode hazard assesment with numerical values: the assessment coding being graded from Low to High, coding numerically adds meaningful information.
[Attribute information](https://archive.ics.uci.edu/ml/datasets/seismic-bumps):
1. seismic: result of shift seismic hazard assessment in the mine working obtained by the seismic
method (a - lack of hazard, b - low hazard, c - high hazard, d - danger state);
2. seismoacoustic: result of shift seismic hazard assessment in the mine working obtained by the
seismoacoustic method;
3. shift: information about type of a shift (W - coal-getting, N -preparation shift);
"""

def numeric_encoder(col_name, df):
    try:
      df[col_name + "_enc_0"] = df[col_name].map({'a': 0, 'b': 1, 'c' : 2, 'W' : 0, 'N' : 1})
      df.drop(columns=col, inplace=True)
      return df
    except Exception as e:
        print(e)
        print("Error in numeric encoding feature : " + col_name)

for col in col_list_categorical:
    df = numeric_encoder(col, df)
df.shape

"""#### One-hot encoder: REPLACED With Numeric encoding above"""

def one_hot_encoder(col_name, df):
    label_encoder = LabelEncoder()
    onehot_encoder = OneHotEncoder(drop='first', sparse=False)
    try:
        encoded_array = label_encoder.fit_transform(df[col_name])
        print(label_encoder.classes_)
#         SKLEARN IMPLEMENTATION
        encoded_array_reshaped = encoded_array.reshape(len(encoded_array), 1)
        one_hot_encoded_array = onehot_encoder.fit_transform(encoded_array_reshaped)
#         one_hot_encoded_array = to_categorical(encoded_array)
        print(one_hot_encoded_array)
        num_features = one_hot_encoded_array.shape[1]
        print("Number of encoded columns to add: ", num_features)
        new_enc_col_names = [col + '_enc_' + str(num) for num in range(0, num_features)] 
        df_enc = pd.DataFrame(one_hot_encoded_array)
        print("Shape of encoded df: ", df_enc.shape)
        df_enc.columns = new_enc_col_names
        print("New column names: ", new_enc_col_names)
        df = pd.concat([df, df_enc], axis=1)
        df.drop(columns=col, inplace=True)
        return df
    except Exception as e:
        print(e)
        print("Error in encoding feature : " + col_name)

#for col in col_list_categorical:
#    df = one_hot_encoder(col, df)
df.shape

df_raw['shift'].value_counts()

df.shift_enc_0.value_counts()

df.head()

df.info()

"""Notice that after one hot encoding, thate object types changed to floats.
### Numeric Feature Scaling
"""

df[col_list_numerical].describe()

df_corr = df[col_list_numerical].corr()

df_corr

plt.figure(figsize=[8, 8])
sns.heatmap(data=df_corr, vmin=-1, vmax=1, cmap='gist_earth_r', annot=True, square=True, linewidths=1)
plt.xticks(rotation=90)
plt.yticks(rotation=0)

"""#### genergy and gpuls"""

sns.distplot(df['genergy'], hist=True)

sns.distplot(np.log(df['genergy']), hist=True)

sns.distplot(df['gpuls'], hist=True)

sns.distplot(np.log(df['gpuls']), hist=True)

plt.figure(figsize=[10, 8])
sns.scatterplot(x='genergy', y='gpuls', hue='class', data=df)

# Applying log transform

df['log_t_genergy'] = np.log(df['genergy'])
df['log_t_gpuls'] = np.log(df['gpuls'])

plt.figure(figsize=[10, 8])
sns.scatterplot(x='log_t_genergy', y='log_t_gpuls', hue='class', data=df)

df.drop(columns=['log_t_genergy', 'log_t_gpuls'], inplace=True)

"""#### gdenergy and gdpuls"""

plt.figure(figsize=[10, 8])
sns.scatterplot(x='gdenergy', y='gdpuls', hue='class', data=df)

"""There are negative shifts in gdenergy and gdpuls because of which log transformation is not going to be the best choice. So, trying out square root transformation."""

sns.distplot(df['gdenergy'], hist=True)

"""There are some negative shifts here."""

sns.distplot(np.sqrt(df['gdenergy']), hist=True)
plt.show()
sns.distplot(np.sqrt(df['gdpuls']), hist=True)

""".. and trying out standard transformation"""

sns.distplot(stats.zscore(df['gdenergy']), hist=True)
plt.show()
sns.distplot(stats.zscore(df['gdpuls']), hist=True)

"""#### energy and maxenergy"""

plt.figure(figsize=[10, 8])
sns.scatterplot(x='energy', y='maxenergy', hue='class', data=df)

"""A strong linear relationship is here betwen the two features."""

sns.distplot(df['energy'], hist=True)
plt.show()
sns.distplot(df['maxenergy'], hist=True)
plt.show()

"""There is a lot of zeroes in this feature. Therefore, using zscore again to check the distribution."""

sns.distplot(stats.zscore(df['energy']), hist=True)
plt.show()
sns.distplot(stats.zscore(df['maxenergy']), hist=True)
plt.show()

stats.zscore(df['energy']).mean(), stats.zscore(df['maxenergy']).mean()

"""Trying log transformation after adding a constant 1"""

sns.distplot(np.log(df['energy']+1), hist=True)
plt.show()
sns.distplot(np.log(df['maxenergy']+1), hist=True)
plt.show()

"""For energy and maxenergy, there is little differences in the distributions as in the plots and there is almost 100% correlation in the data - so I would drop maxenergy and keep energy only. I would translate using a constant 1 and then use log transformation on energy.
**Summarizing the chosen transformation in a dictionary**
"""

def shifted_log_func(df_col):
    return np.log(1 + df_col)

dict_num_cols_trnsfm = {'genergy': np.log,
                        'gpuls' : np.log,
                        'gdenergy': stats.zscore,
                        'gdenergy': stats.zscore, 
                        'energy': shifted_log_func}

for col_names, transfm_func in dict_num_cols_trnsfm.items():
    df['scaled_' + col_names] = transfm_func(df[col_names])
df.drop(columns=col_list_numerical, inplace=True)
df.head()

df[[col for col in df.columns if 'scaled_' in col]].describe()

"""### Discrete Features"""

for col in col_list_discrete:
    sns.countplot(df[col])
    plt.show()

"""Since nbumps 6 , 7 & 89 contain zero, I will drop them as they have no information."""

df.drop(columns=['nbumps6', 'nbumps7', 'nbumps89'], inplace=True)

df.head()

col_list_discrete = list(set(col_list_discrete) - set(['nbumps6', 'nbumps7', 'nbumps89']))
col_list_discrete

for each_col in col_list_discrete:
    data_crosstab = pd.crosstab(df[each_col], df[label], colnames=['class'])
    print(data_crosstab)
    print('-----')

"""## Correlations between features and the target"""

df.shape

df.describe()

df_corr = df.corr()

# Generate a mask for the upper triangle
mask = np.triu(np.ones_like(df_corr, dtype=bool))

plt.figure(figsize=[12, 12])
sns.heatmap(data=df_corr, mask=mask, vmin=-1, vmax=1, cmap='gist_earth_r', annot=True, square=True, linewidths=1)
plt.xticks(rotation=90)
plt.yticks(rotation=0)

"""## Splitting Traning and Test data using Stratified Shuffle Split"""

df.shape

X_cols = list(set(df.columns) - set([label]))
X_cols

# OLD CODE FOR DEMO
X_demo = df[list(set(df.columns) - set([label]))]
y_demo1 = df[label]
y_demo2 = df[[label]]
X_demo.shape, y_demo1.shape, y_demo2.shape

type(X_demo), type(y_demo1), type(y_demo2)

X = df[list(set(df.columns) - set([label]))].values
y = df[label]
X.shape, y.shape

type(X), type(y)

stratified_split = StratifiedShuffleSplit(n_splits=1, test_size=0.20)

for train_idx, test_idx in stratified_split.split(X, y):
    print(len(train_idx))
    print(len(test_idx))
    y_train= y[train_idx]
    X_train = X[train_idx]
    X_test, y_test = X[test_idx], y[test_idx]

len(X_train), len(y_train), len(X_test), len(y_test)

print("Training Set Target Class Distribution:")
print(y_train.value_counts()/len(y_train))
print("Test Set Target Class Distribution:")
print(y_test.value_counts()/len(y_test))

X_train = pd.DataFrame(X_train, columns=X_cols)
X_test = pd.DataFrame(X_test, columns=X_cols)

X_train_demo, X_test_demo, y_train_demo, y_test_demo = train_test_split(X, y, test_size=0.33, random_state=42, shuffle=False)
X_train_demo.shape, X_test_demo.shape, y_train_demo.shape, y_test_demo.shape, y_train_demo.value_counts()/len(y_train_demo), y_test_demo.value_counts()/len(y_test_demo)

print("Training Set Target Class Distribution:")
print(y_train_demo.value_counts()/len(y_train_demo))
print("Test Set Target Class Distribution:")
print(y_test_demo.value_counts()/len(y_test_demo))

X_train_demo, X_test_demo, y_train_demo, y_test_demo = train_test_split(X, y, test_size=0.33, 
                                                                        random_state=42, shuffle=True)
X_train_demo.shape, X_test_demo.shape, y_train_demo.shape, y_test_demo.shape, y_train_demo.value_counts()/len(y_train_demo), y_test_demo.value_counts()/len(y_test_demo)

print("Training Set Target Class Distribution:")
print(y_train_demo.value_counts()/len(y_train_demo))
print("Test Set Target Class Distribution:")
print(y_test_demo.value_counts()/len(y_test_demo))

"""## NEW Generate Training Set with high Class 1 with SMOTE"""

counter = Counter(y_train)
print(counter)

#sns.scatterplot(x='scaled_genergy', y='scaled_gpuls', hue=class, data=X_train)

# View current data
plt.figure(figsize=[10, 8])


for lb, _ in counter.items():
	row_ix = np.where(y_train == lb)[0]
	plt.scatter(X_train.loc[row_ix, ['scaled_genergy']], 
	            X_train.loc[row_ix, ['scaled_gpuls']], label=str(lb))

	#plt.scatter(X[y == label].scaled_genergy, X[y == label].scaled_gpuls, label=str(label))
plt.title("Original data: scaled_genergy vs. scaled_gpuls")
plt.legend()
plt.show()

plt.figure(figsize=[10, 8])

sns.distplot(X_train['nbumps'])
plt.title("Original data: nbumps distribution")

plt.show()

print('Original dataset shape %s' % Counter(y))

sm = SMOTE(random_state=42)
X_res, y_res = sm.fit_resample(X_train, y_train)

print('Resampled dataset shape %s' % Counter(y_res))

X_res = pd.DataFrame(X_res, columns=X_cols)

plt.figure(figsize=[10, 8])
plt.title("Synthetic Minority Oversampling Technique data: scaled_genergy vs. scaled_gpuls")

for lb, _ in counter.items():
	row_ix = np.where(y_res == lb)[0]
	plt.scatter(X_res.loc[y_res == lb, ['scaled_genergy']], 
	            X_res.loc[y_res == lb, ['scaled_gpuls']], label=str(lb))



plt.legend()
plt.show()

"""*Some non-natural number of Bumps! (float values!)*"""

plt.figure(figsize=[10, 8])

sns.distplot(X_res['nbumps'])
plt.title("SMOTE data: nbumps distribution")

plt.show()



"""## Save models to a particular folder"""

import datetime

def model_store_location():
    return "model-store-{}".format(datetime.date.today())

model_store_location = model_store_location()
print(model_store_location)
!mkdir {model_store_location}
model_store_location += '/'
model_store_location

"""## Applying ML model to Predict Hazardrous Seismic Bump
### K-Neighbors Classifier using Stratified KFold Cross Validation
#### Establishing Baseline
"""

model = KNeighborsClassifier()
skf = StratifiedKFold(n_splits=5)

'''
Demo of split using original features and label columns
split is 90-10
'''

for train, test in skf.split(X, y):
    print(len(train))
    print(y[train].value_counts()/len(train))
    print('----')
    print(len(test))
    print(y[test].value_counts()/len(test))
    print('----')
    break

"""#### Using cross_val_score for performance metrics"""

scores = cross_val_score(model, X_train, y_train, scoring='accuracy', cv=skf, n_jobs=-1, error_score='raise')
scores

scores.mean(), scores.std()

scores = cross_val_score(model, X_train, y_train, scoring='f1', cv=skf, n_jobs=-1, error_score='raise')
scores.mean(), scores.std()

"""I was getting an error "pos_label=1 is not a valid label: array(['0', '1'], dtype='<U1')" which I resolved by typecasting the labels into integers. This blocked me for 1 whole day sicne I had run the code before which successfully ran. This is when i noticed one-hot encoders transformed all the features I fed into it were transformed to float including my label column. So when I ran the code, it worked. But when I  excluded it, the object type of the label column caused this trouble."""

model.n_neighbors, model.weights

t0= time.perf_counter()

model.fit(X_train, y_train)

t1 = time.perf_counter() - t0

y_pred = model.predict(X_test)

confusion_matrix(y_test, y_pred)

y_test.value_counts()

pd.Series( y_pred).value_counts()

"""#### Add model results to a DataFrame"""

df_results = pd.DataFrame(index = ['scoring_technique', 'algorithm', 'n_neighbors', 'weights', 'leaf_size', 'accuracy', 'precision', 'recall', 'f1-score', 'Time_to_run_seconds'])

df_results['KNC Baseline'] = ['None', model.algorithm,
                          model.n_neighbors, 
                          model.weights, 
                          model.leaf_size,
                          accuracy_score(y_test, y_pred), 
                          precision_score(y_test, y_pred), 
                          recall_score(y_test, y_pred), 
                          f1_score(y_test, y_pred), t1]

df_results

import sklearn
sorted(sklearn.metrics.SCORERS.keys())

"""#### Using cross_validate for performance metrics"""

scores = cross_validate(model, X_train, y_train, scoring=['precision', 'recall', 'f1'], cv=skf, n_jobs=-1, error_score='raise')
scores.keys()

scores['test_precision'].mean(), scores['test_recall'].mean(), scores['test_f1'].mean()

scoring = {'precision': 'precision', 
           'recall': 'recall',
           'roc_auc': 'roc_auc'}

scores = cross_validate(model, X_train, y_train, scoring=scoring, cv=10, return_train_score=True)
sorted(scores.keys())

scores['test_precision'].mean(), scores['test_recall'].mean(), scores['test_roc_auc'].mean()

"""#### Tuning model using Grid Search CV using 'Precision'"""

n_neighbors = [1, 2, 3, 4, 5]
weights = ['uniform', 'distance']
algorithm = ['ball_tree', 'kd_tree', 'brute']
leaf_size = [10, 20, 30 , 40, 50]

param_grid = dict(n_neighbors=n_neighbors, weights=weights, algorithm=algorithm, leaf_size=leaf_size)

grid = GridSearchCV(estimator=model, 
                    param_grid=param_grid, 
                    cv=StratifiedKFold(shuffle=True), 
                    scoring=['precision'],
                    refit='precision',
                    verbose=1,
                    n_jobs=-1)

t0= time.perf_counter()
grid_result = grid.fit(X_train, y_train)

t1 = time.perf_counter() - t0

grid_result.best_estimator_

print("Time elapsed: ", t1) # CPU seconds elapsed (floating point)

grid_result.best_estimator_

grid_result.best_score_

n_neighbors, weights = grid_result.best_estimator_.n_neighbors , grid_result.best_estimator_.weights

n_neighbors, weights

grid_result.scorer_

grid_result.n_splits_

model = grid_result.best_estimator_

model

file_name = 'seismic_bump_k_neighbors_'+ str(n_neighbors) + '_' + weights + '_' + list(grid_result.scorer_.keys())[0]+'.sav'
joblib.dump(model, model_store_location + file_name)

y_pred = model.predict(X_test)

y_pred.shape

confusion_matrix(y_test, y_pred)

tn, fp, fn, tp = confusion_matrix(y_test, y_pred).ravel()
(tn, fp, fn, tp)

accuracy_score(y_test, y_pred), recall_score(y_test, y_pred), precision_score(y_test, y_pred), f1_score(y_test, y_pred)

df_results['KNC Model 1'] = [grid.get_params(deep=True)['scoring'], grid_result.best_estimator_.algorithm,
                         n_neighbors, weights, grid_result.best_estimator_.leaf_size,
                         accuracy_score(y_test, y_pred), 
                         precision_score(y_test, y_pred), recall_score(y_test, y_pred), 
                         f1_score(y_test, y_pred), t1]

df_results

"""#### Tuning model using Grid Search CV using 'Precision'  again"""

t0= time.perf_counter()

grid = GridSearchCV(estimator=model, 
                    param_grid=param_grid, 
                    cv=StratifiedKFold(shuffle=True), 
                    scoring=['precision'],
                    refit='precision',
                    verbose=10, n_jobs=-1)
grid_result = grid.fit(X_train, y_train)

t1 = time.perf_counter() - t0

n_neighbors, weights = grid_result.best_estimator_.n_neighbors , grid_result.best_estimator_.weights
file_name = 'seismic_bump_k_neighbors_'+ str(n_neighbors) + '_' + weights + '_' + list(grid_result.scorer_.keys())[0]+'.sav'
joblib.dump(model, model_store_location + file_name)
print(grid_result.best_score_)
print(grid_result.best_estimator_)
y_pred = grid_result.best_estimator_.predict(X_test)

accuracy_score(y_test, y_pred), recall_score(y_test, y_pred), precision_score(y_test, y_pred), f1_score(y_test, y_pred)

df_results['KNC Model 2'] = [grid.get_params(deep=True)['scoring'], grid_result.best_estimator_.algorithm,
                         n_neighbors, weights, grid_result.best_estimator_.leaf_size,
                         accuracy_score(y_test, y_pred), precision_score(y_test, y_pred), 
                         recall_score(y_test, y_pred), f1_score(y_test, y_pred), t1]

"""#### Tuning model using Grid Search CV using 'Recall'"""

t0= time.perf_counter()

grid = GridSearchCV(estimator=model, 
                    param_grid=param_grid, 
                    cv=StratifiedKFold(shuffle=True), 
                    scoring=['recall'],
                    refit='recall',
                    verbose=1, n_jobs=-1)
grid_result = grid.fit(X_train, y_train)

t1 = time.perf_counter() - t0

n_neighbors, weights = grid_result.best_estimator_.n_neighbors , grid_result.best_estimator_.weights
file_name = 'seismic_bump_k_neighbors_'+ str(n_neighbors) + '_' + weights + '_' + list(grid_result.scorer_.keys())[0]+'.sav'
joblib.dump(model, model_store_location + file_name)
print(grid_result.best_score_)
print(grid_result.best_estimator_)
y_pred = grid_result.best_estimator_.predict(X_test)

accuracy_score(y_test, y_pred), recall_score(y_test, y_pred), precision_score(y_test, y_pred), f1_score(y_test, y_pred)

df_results['KNC Model 3'] = [grid.get_params(deep=True)['scoring'], grid_result.best_estimator_.algorithm,
                         n_neighbors, weights, grid_result.best_estimator_.leaf_size,
                         accuracy_score(y_test, y_pred), 
                         precision_score(y_test, y_pred), recall_score(y_test, y_pred), f1_score(y_test, y_pred), t1]

df_results

"""#### Tuning model using Grid Search CV using 'f1'"""

t0= time.perf_counter()
grid = GridSearchCV(estimator=model, 
                    param_grid=param_grid, 
                    cv=StratifiedKFold(shuffle=True), 
                    scoring=['f1'],
                    refit='f1',
                    verbose=1)
grid_result = grid.fit(X_train, y_train)

t1 = time.perf_counter() - t0

n_neighbors, weights = grid_result.best_estimator_.n_neighbors , grid_result.best_estimator_.weights
file_name = 'seismic_bump_k_neighbors_'+ str(n_neighbors) + '_' + weights + '_' + list(grid_result.scorer_.keys())[0]+'.sav'
joblib.dump(model, model_store_location + file_name)
print(grid_result.best_score_)
print(grid_result.best_estimator_)
y_pred = grid_result.best_estimator_.predict(X_test)

accuracy_score(y_test, y_pred), recall_score(y_test, y_pred), precision_score(y_test, y_pred), f1_score(y_test, y_pred)

df_results['KNC Model 4'] = [grid.get_params(deep=True)['scoring'], grid_result.best_estimator_.algorithm,
                         n_neighbors, weights, grid_result.best_estimator_.leaf_size, 
                         accuracy_score(y_test, y_pred), precision_score(y_test, y_pred), 
                         recall_score(y_test, y_pred), f1_score(y_test, y_pred), t1]

"""#### Tuning model using Grid Search CV using 'Precision' again"""

t0= time.perf_counter()

grid = GridSearchCV(estimator=model, 
                    param_grid=param_grid, 
                    cv=StratifiedKFold(shuffle=True), 
                    scoring=['precision'],
                    refit='precision',
                    verbose=10, n_jobs=-1)
grid_result = grid.fit(X_train, y_train)

t1 = time.perf_counter() - t0

n_neighbors, weights = grid_result.best_estimator_.n_neighbors , grid_result.best_estimator_.weights
file_name = 'seismic_bump_k_neighbors_'+ str(n_neighbors) + '_' + weights + '_' + list(grid_result.scorer_.keys())[0]+'.sav'
print(file_name)
joblib.dump(model, model_store_location + file_name)
print(grid_result.best_score_)
print(grid_result.best_estimator_)
y_pred = grid_result.best_estimator_.predict(X_test)

accuracy_score(y_test, y_pred), recall_score(y_test, y_pred), precision_score(y_test, y_pred), f1_score(y_test, y_pred)

df_results['KNC Model 5'] = [grid.get_params(deep=True)['scoring'], grid_result.best_estimator_.algorithm,
                         n_neighbors, weights, grid_result.best_estimator_.leaf_size,
                         accuracy_score(y_test, y_pred), 
                         precision_score(y_test, y_pred), recall_score(y_test, y_pred), f1_score(y_test, y_pred), t1]

"""#### Run model with SMOTE augmented dataset"""

#model = KNeighborsClassifier(n_neighbors=1, leaf_size=10)
t0= time.perf_counter()

model = KNeighborsClassifier(algorithm='ball_tree', leaf_size=10, n_neighbors=4, weights='distance')
model.fit(X_res, y_res)

t1 = time.perf_counter() - t0

y_pred = model.predict(X_test)

y_ptraining = model.predict(X_res)
df_results['KNC SMOTE 1'] = ['None', model.algorithm,
                          model.n_neighbors, 
                          model.weights, 
                          model.leaf_size, 
                          accuracy_score(y_test, y_pred),
                          precision_score(y_test, y_pred), 
                          recall_score(y_test, y_pred), 
                          f1_score(y_test, y_pred), t1]

df_results

confusion_matrix(y_test, y_pred)

"""#### Tuning model using SMOTE dataset and Grid Search CV using 'Recall'"""

t0= time.perf_counter()

grid = GridSearchCV(estimator=model, 
                    param_grid=param_grid, 
                    cv=StratifiedKFold(shuffle=True), 
                    scoring=['recall'],
                    refit='recall',
                    verbose=1, n_jobs=-1)
grid_result = grid.fit(X_res, y_res)

t1 = time.perf_counter() - t0

n_neighbors, weights = grid_result.best_estimator_.n_neighbors , grid_result.best_estimator_.weights
file_name = 'seismic_bump_k_neighbors_'+ str(n_neighbors) + '_' + weights + '_' + list(grid_result.scorer_.keys())[0]+'.sav'
joblib.dump(model, model_store_location + file_name)
print(grid_result.best_score_)
print(grid_result.best_estimator_)
y_pred = grid_result.best_estimator_.predict(X_test)

confusion_matrix(y_test, y_pred)

df_results['KNC Optimised SMOTE'] = ['[recall]', model.algorithm,
                          model.n_neighbors, 
                          model.weights, 
                          model.leaf_size, 
                          accuracy_score(y_test, y_pred),
                          precision_score(y_test, y_pred), 
                          recall_score(y_test, y_pred), 
                          f1_score(y_test, y_pred), t1]

df_results

"""## Random Forest
### Optimise for Recall with raw Dataset
"""

from sklearn.ensemble import RandomForestClassifier
model = RandomForestClassifier()

param_grid = {
    'bootstrap': [True],
    'max_depth': [3, 5, 10],
    'max_features': ['auto', 'sqrt'],
    'min_samples_leaf': [3, 4, 5, 6],
    'min_samples_split': [4, 8, 10, 12],
    'n_estimators': [100, 200, 300, 500]
}

# WAS     'max_features': [2, 3, 4],
#     'min_samples_split': [8, 10, 12],
small_param_grid = {
    'bootstrap': [True],
    'max_depth': [3, 5, 10]
}

small_param_grid

X_train.shape

# was GridSearchCV

from sklearn.model_selection import RandomizedSearchCV

t0= time.perf_counter()

grid = RandomizedSearchCV(estimator=model, 
                    param_distributions=param_grid, 
                    cv=StratifiedKFold(shuffle=True), 
                    scoring=['recall'],
                    refit='recall',
                    n_iter = 100, verbose=1, random_state=42, n_jobs = -1)

grid_result = grid.fit(X_train, y_train)

t1 = time.perf_counter() - t0

# summarize the best score and configuration
print("Best: %f using %s" % (grid_result.best_score_, grid_result.best_params_))
# summarize all scores that were evaluated
means = grid_result.cv_results_['mean_test_recall']
stds = grid_result.cv_results_['std_test_recall']
params = grid_result.cv_results_['params']
#for mean, stdev, param in zip(means, stds, params):
    #print("%f (%f) with: %r" % (mean, stdev, param))

grid_result.best_estimator_

file_name = 'seismic_bump_rf_model.sav'
joblib.dump(model, model_store_location + file_name)

grid_result.best_estimator_

"""```
RandomForestClassifier(bootstrap=True, ccp_alpha=0.0, class_weight=None,
                       criterion='gini', max_depth=10, max_features='sqrt',
                       max_leaf_nodes=None, max_samples=None,
                       min_impurity_decrease=0.0, min_impurity_split=None,
                       min_samples_leaf=3, min_samples_split=4,
                       min_weight_fraction_leaf=0.0, n_estimators=300,
                       n_jobs=None, oob_score=False, random_state=None,
                       verbose=0, warm_start=False)
```
"""

y_pred = grid_result.best_estimator_.predict(X_test)
accuracy_score(y_test, y_pred), recall_score(y_test, y_pred), precision_score(y_test, y_pred), f1_score(y_test, y_pred)

confusion_matrix(y_test, y_pred)

"""### Try a few params"""

def rf(xs, y, xt, yt, n_estimators=100, max_depth=80,
       max_features='sqrt', min_samples_leaf=3, min_samples_split=8, **kwargs):
    model = RandomForestClassifier(n_jobs=-1, n_estimators=n_estimators, max_depth=max_depth,
        max_features=max_features, min_samples_split=min_samples_split,
        min_samples_leaf=min_samples_leaf, oob_score=True)
    model.fit(xs, y)
    yp = model.predict(xt)
    print(confusion_matrix(yt, yp))
    return recall_score(yt, yp)

plt.plot([rf(X_train, y_train, X_test, y_test, max_features='auto', max_depth=i) for i in [80, 100, 120, 200]])

"""### Optimise for Recall with SMOTE Dataset - VERY LONG , GRID SEARCH COMMENTED OUT"""

t0= time.perf_counter()

model = RandomForestClassifier()

grid = GridSearchCV(estimator=model, 
                    param_grid=param_grid, 
                    cv=StratifiedKFold(shuffle=True), 
                    scoring=['recall'],
                    refit='recall',
                    verbose=2,
                    n_jobs=-1)
#grid_result = grid.fit(X_res, y_res)

t1 = time.perf_counter() - t0


# summarize the best score and configuration
print("Best: %f using %s" % (grid_result.best_score_, grid_result.best_params_))
# summarize all scores that were evaluated
#means = grid_result.cv_results_['mean_test_recall']
#stds = grid_result.cv_results_['std_test_recall']
#params = grid_result.cv_results_['params']
#for mean, stdev, param in zip(means, stds, params):
#    print("%f (%f) with: %r" % (mean, stdev, param))

#grid_result.best_estimator_

"""Result:
```
Fitting 5 folds for each of 384 candidates, totalling 1920 fits
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 2 concurrent workers.
[Parallel(n_jobs=-1)]: Done  37 tasks      | elapsed:   19.9s
[Parallel(n_jobs=-1)]: Done 158 tasks      | elapsed:  1.4min
[Parallel(n_jobs=-1)]: Done 361 tasks      | elapsed:  3.2min
[Parallel(n_jobs=-1)]: Done 644 tasks      | elapsed:  5.7min
[Parallel(n_jobs=-1)]: Done 1009 tasks      | elapsed:  9.5min
[Parallel(n_jobs=-1)]: Done 1454 tasks      | elapsed: 14.7min
[Parallel(n_jobs=-1)]: Done 1920 out of 1920 | elapsed: 20.8min finished
RandomForestClassifier(bootstrap=True, ccp_alpha=0.0, class_weight=None,
                       criterion='gini', max_depth=10, max_features='sqrt',
                       max_leaf_nodes=None, max_samples=None,
                       min_impurity_decrease=0.0, min_impurity_split=None,
                       min_samples_leaf=4, min_samples_split=10,
                       min_weight_fraction_leaf=0.0, n_estimators=100,
                       n_jobs=None, oob_score=False, random_state=None,
                       verbose=0, warm_start=False)
Best: 0.917663 using {'bootstrap': True, 'max_depth': 10, 'max_features': 'auto', 'min_samples_leaf': 3, 'min_samples_split': 8, 'n_estimators': 100}
```
"""

model = RandomForestClassifier(bootstrap=True, ccp_alpha=0.0, class_weight=None,
                       criterion='gini', max_depth=10, max_features='auto',
                       max_leaf_nodes=None, max_samples=None,
                       min_impurity_decrease=0.0, min_impurity_split=None,
                       min_samples_leaf=3, min_samples_split=8,
                       min_weight_fraction_leaf=0.0, n_estimators=100,
                       n_jobs=None, oob_score=False, random_state=None,
                       verbose=0, warm_start=False)

model.fit(X_res, y_res)

y_pres = model.predict(X_res)
print( confusion_matrix(y_res, y_pres))
print( recall_score(y_res, y_pres))

y_pred = model.predict(X_test)

print( confusion_matrix(y_test, y_pred))

df_results['RFC SMOTE'] = ['Recall', '',
                          'max_features=' + str(model.max_features), 
                          'min_samples_split=' + str(model.min_samples_split), 
                          'min_samples_leaf=' + str(model.min_samples_leaf),
                          accuracy_score(y_test, y_pred), 
                          precision_score(y_test, y_pred), 
                          recall_score(y_test, y_pred), 
                          f1_score(y_test, y_pred), 'N\A']

model = RandomForestClassifier(bootstrap=True, ccp_alpha=0.0, class_weight=None,
                       criterion='gini', max_depth=10, max_features='sqrt',
                       max_leaf_nodes=None, max_samples=None,
                       min_impurity_decrease=0.0, min_impurity_split=None,
                       min_samples_leaf=4, min_samples_split=10,
                       min_weight_fraction_leaf=0.0, n_estimators=100,
                       n_jobs=None, oob_score=False, random_state=None,
                       verbose=0, warm_start=False)

model.fit(X_res, y_res)

y_pres = model.predict(X_res)
y_pred = model.predict(X_test)
print("Training set Recall score: {}, Test set Recall score: {}".format( recall_score(y_res, y_pres), recall_score(y_test, y_pred)))


print( confusion_matrix(y_test, y_pred))

df_results['RFC SMOTE'] = ['Recall', '',
                          'max_features=' + str(model.max_features), 
                          'min_samples_split=' + str(model.min_samples_split), 
                          'min_samples_leaf=' + str(model.min_samples_leaf),
                          accuracy_score(y_test, y_pred), 
                          precision_score(y_test, y_pred), 
                          recall_score(y_test, y_pred), 
                          f1_score(y_test, y_pred), 'N\A']

"""### Looks like optimised RFC overfitting, try out a smaller max_leaf_nodes"""

model = RandomForestClassifier(max_leaf_nodes=10)
model.fit(X_res, y_res)

y_pres = model.predict(X_res)
y_pred = model.predict(X_test)
print("Training set Recall score: {}, Test set Recall score: {}".format( recall_score(y_res, y_pres), recall_score(y_test, y_pred)))

print( confusion_matrix(y_test, y_pred))

df_results['RFC SMOTE max_leaf_nodes=10'] = ['max_depth' + str(model.max_depth),                                   
                          'max_leaf_nodes=' + str(model.max_leaf_nodes),
                          'max_features=' + str(model.max_features), 
                          'min_samples_split=' + str(model.min_samples_split), 
                          'min_samples_leaf=' + str(model.min_samples_leaf),
                          accuracy_score(y_test, y_pred), 
                          precision_score(y_test, y_pred), 
                          recall_score(y_test, y_pred), 
                          f1_score(y_test, y_pred), 'N\A']

20 / (20 + 14)

model = RandomForestClassifier(max_depth=3)
model.fit(X_res, y_res)

y_pred = model.predict(X_test)

print( confusion_matrix(y_test, y_pred))

df_results['RFC SMOTE max_depth=3'] = ['max_depth' + str(model.max_depth),                                   
                          'max_leaf_nodes=' + str(model.max_leaf_nodes),
                          'max_features=' + str(model.max_features), 
                          'min_samples_split=' + str(model.min_samples_split), 
                          'min_samples_leaf=' + str(model.min_samples_leaf),
                          accuracy_score(y_test, y_pred), 
                          precision_score(y_test, y_pred), 
                          recall_score(y_test, y_pred), 
                          f1_score(y_test, y_pred), 'N\A']

df_results.T.to_csv(model_store_location + "FinalResults.csv")
df_results

"""### Feature importances"""

X_cols

rfc_fi = pd.DataFrame({'cols':X_cols, 'imp':model.feature_importances_}
                       ).sort_values('imp', ascending=False)
rfc_fi

"""## SVM"""

from sklearn.svm import SVC 
model = SVC() 

param_grid = {'C': [0.1, 1, 10, 100, 1000],  
              'gamma': [1, 0.1, 0.01, 0.001, 0.0001], 
              'kernel': [ 'linear', 'poly', 'rbf', 'sigmoid']}

"""### SVM Optimisation with Raw dataset
Long optimisation, comment out here
"""

t0= time.perf_counter()

grid = GridSearchCV(estimator=model, 
                    param_grid=param_grid, 
                    cv=StratifiedKFold(shuffle=True), 
                    scoring=['f1'],
                    refit='f1',
                    verbose=1, n_jobs=-1)
#grid_result = grid.fit(X_train, y_train)

t1 = time.perf_counter() - t0

#grid_result.best_estimator_

# defining parameter range 
small_param_grid = {'C': [1],  
              'gamma': [1], 
              'kernel': ['poly', 'rbf', 'sigmoid']}

t0= time.perf_counter()
grid = GridSearchCV(estimator=model, 
                    param_grid=small_param_grid, 
                    cv=StratifiedKFold(shuffle=True), 
                    scoring=['f1'],
                    refit='f1',
                    verbose=3, n_jobs=-1)
#grid_result = grid.fit(X_train, y_train)

t1 = time.perf_counter() - t0

#grid_result.best_estimator_

"""Result:
```
Fitting 5 folds for each of 3 candidates, totalling 15 fits
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 2 concurrent workers.
[Parallel(n_jobs=-1)]: Done  15 out of  15 | elapsed: 83.6min finished
SVC(C=1, break_ties=False, cache_size=200, class_weight=None, coef0=0.0,
    decision_function_shape='ovr', degree=3, gamma=1, kernel='poly',
    max_iter=-1, probability=False, random_state=None, shrinking=True,
    tol=0.001, verbose=False)
```
"""

model = grid_result.best_estimator_

file_name = 'seismic_bump_svm_model.sav'
joblib.dump(model, model_store_location + file_name)

#grid_result.best_estimator_

model = SVC(C=1, break_ties=False, cache_size=200, class_weight=None, coef0=0.0,
    decision_function_shape='ovr', degree=3, gamma=1, kernel='poly',
    max_iter=-1, probability=False, random_state=None, shrinking=True,
    tol=0.001, verbose=False)

model.fit(X_train, y_train)

y_pred = model.predict(X_test)
accuracy_score(y_test, y_pred), recall_score(y_test, y_pred), precision_score(y_test, y_pred), f1_score(y_test, y_pred)

confusion_matrix(y_test, y_pred)

df_results['SVC Original Dataset'] = ['f1', 'SVM', 'C: {}'.format(model.get_params()['C']), 
                                      'gamma: {}'.format(model.get_params()['gamma']),
                                      'kernel: {}'.format(model.get_params()['kernel']),
                                      accuracy_score(y_test, y_pred), 
                                      precision_score(y_test, y_pred), 
                                      recall_score(y_test, y_pred), 
                                      f1_score(y_test, y_pred), (83*60)]

df_results

"""### SVM Optimisation with SMOTE Augmented dataset - VERY LONG GRIDSEARCH -> COMMENTED OUT , results in"""

model = SVC() 

param_grid = {'C': [0.1, 1, 10],  
              'gamma': [1, 0.1, 0.01], 
              'kernel': [ 'linear', 'poly', 'rbf', 'sigmoid']}

# defining parameter range 
t0= time.perf_counter()

grid = GridSearchCV(estimator=model, 
                    param_grid=param_grid, 
                    cv=StratifiedKFold(shuffle=True), 
                    scoring=['recall'],
                    refit='recall',
                    verbose=2, n_jobs=-1)
#grid_result = grid.fit(X_res, y_res)

t1 = time.perf_counter() - t0


# summarize the best score and configuration
print("Best: %f using %s" % (grid_result.best_score_, grid_result.best_params_))
# summarize all scores that were evaluated
#means = grid_result.cv_results_['mean_test_score']
#stds = grid_result.cv_results_['std_test_score']
#params = grid_result.cv_results_['params']
#for mean, stdev, param in zip(means, stds, params):
#    print("%f (%f) with: %r" % (mean, stdev, param))

#grid_result.best_estimator_

"""Results:
```
Fitting 5 folds for each of 36 candidates, totalling 180 fits
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 2 concurrent workers.
[Parallel(n_jobs=-1)]: Done  37 tasks      | elapsed: 104.3min
[Parallel(n_jobs=-1)]: Done 158 tasks      | elapsed: 614.2min
[Parallel(n_jobs=-1)]: Done 180 out of 180 | elapsed: 615.1min finished
Best: 0.973589 using {'C': 10, 'gamma': 1, 'kernel': 'rbf'}
```
"""

model = SVC(C=10, break_ties=False, cache_size=200, class_weight=None, coef0=0.0,
    decision_function_shape='ovr', degree=3, gamma=1, kernel='rbf', max_iter=-1,
    probability=False, random_state=None, shrinking=True, tol=0.001,
    verbose=True)
model.fit(X_res, y_res)

y_pred = model.predict(X_test)
print(confusion_matrix(y_test, y_pred))

y_pres = model.predict(X_res)
print(confusion_matrix(y_res, y_pres))

df_results['SVC SMOTE'] = ['[recall]', 'SVM', 'C: {}'.format(model.get_params()['C']), 
                           'gamma: {}'.format(model.get_params()['gamma']),
                           'kernel: {}'.format(model.get_params()['kernel']),
                          accuracy_score(y_test, y_pred), 
                          precision_score(y_test, y_pred), 
                          recall_score(y_test, y_pred), 
                          f1_score(y_test, y_pred), (530*60)]

df_results

"""# Decision Tree
Optimisation with Raw dataset
"""

from sklearn.tree import DecisionTreeClassifier, export_graphviz

t0= time.perf_counter()

DT_classifier = DecisionTreeClassifier()
DTree = DT_classifier.fit(X_train, y_train)

t1 = time.perf_counter() - t0

predictedDT = DT_classifier.predict(X_test)

from sklearn.metrics import accuracy_score, classification_report, confusion_matrix

DT_result = accuracy_score(y_test, predictedDT) * 100
print("Accuracy score using DecisionTreeClassifier:", DT_result)
print(confusion_matrix(y_test, predictedDT))
print(classification_report(y_test, predictedDT))

df_results['Decision Tree Raw Data'] = ['None', '',
                          'max_features=' + str(DT_classifier.max_features), 
                          'min_samples_split=' + str(DT_classifier.min_samples_split), 
                          'min_samples_leaf=' + str(DT_classifier.min_samples_leaf),
                          accuracy_score(y_test, predictedDT),
                          precision_score(y_test, predictedDT), 
                          recall_score(y_test, predictedDT), 
                          f1_score(y_test, predictedDT), t1]

df_results

from sklearn.externals.six import StringIO
from IPython.display import Image 
import pydotplus

# visualizing data 
dot_data = StringIO()

export_graphviz(
    DTree,
    out_file = dot_data,
    rounded = True, 
    special_characters = True,
    filled = True)

graph = pydotplus.graph_from_dot_data(dot_data.getvalue())
Image(graph.create_png())

"""Optimisation with SMOTE Augmented dataset"""

t0= time.perf_counter()

DTreeSM = DT_classifier.fit(X_res, y_res)

t1 = time.perf_counter() - t0

predictedDTSM = DT_classifier.predict(X_test)

SmoteDT = accuracy_score(y_test, predictedDTSM) * 100
print("Accuracy score using DecisionTreeClassifier:", SmoteDT)
print(confusion_matrix(y_test, predictedDTSM))
print(classification_report(y_test, predictedDTSM))

df_results['Decision Tree SMOTE Data'] = ['None', '',
                          'max_features=' + str(DTreeSM.max_features), 
                          'min_samples_split=' + str(DTreeSM.min_samples_split), 
                          'min_samples_leaf=' + str(DTreeSM.min_samples_leaf),
                          accuracy_score(y_test, predictedDTSM),
                          precision_score(y_test, predictedDTSM), 
                          recall_score(y_test, predictedDTSM), 
                          f1_score(y_test, predictedDTSM), t1]

df_results

# visualizing data 
dot_data = StringIO()

export_graphviz(
    DTreeSM,
    out_file = dot_data,
    rounded = True, 
    special_characters = True,
    filled = True)

graph = pydotplus.graph_from_dot_data(dot_data.getvalue())
Image(graph.create_png())

"""### max_leaf_nodes=10 for Visualisation and to avoid overfitting"""

t0= time.perf_counter()

DTM_classifier = DecisionTreeClassifier(max_leaf_nodes=10)
DTreeSM = DTM_classifier.fit(X_res, y_res)

t1 = time.perf_counter() - t0

predictedDTSM = DTM_classifier.predict(X_test)

df_results['DecisionTree SMOTE max_leaf_nodes=10'] = ['None', '',
                          'max_features=' + str(DTreeSM.max_features), 
                          'min_samples_split=' + str(DTreeSM.min_samples_split), 
                          'min_samples_leaf=' + str(DTreeSM.min_samples_leaf),
                          accuracy_score(y_test, predictedDTSM),
                          precision_score(y_test, predictedDTSM), 
                          recall_score(y_test, predictedDTSM), 
                          f1_score(y_test, predictedDTSM), t1]

df_results

# visualizing data 
dot_data = StringIO()

export_graphviz(
    DTreeSM,
    out_file = dot_data,
    feature_names=X_cols,
    rounded = True, 
    special_characters = True,
    filled = True)

graph = pydotplus.graph_from_dot_data(dot_data.getvalue())
Image(graph.create_png())

"""Diagram shows 0 and 1 values"""

Counter( y_res[(X_res.nbumps > 1.0) & (X_res.seismic_enc_0 <= 0.001)])

Counter(y_res)

t0= time.perf_counter()

DTM_classifier = DecisionTreeClassifier(max_depth=3)
DTreeSM = DTM_classifier.fit(X_res, y_res)

t1 = time.perf_counter() - t0

predictedDTSM = DTM_classifier.predict(X_test)

df_results['DecisionTree SMOTE max_depth=3'] = ['None', '',
                          'max_features=' + str(DTreeSM.max_features), 
                          'min_samples_split=' + str(DTreeSM.min_samples_split), 
                          'min_samples_leaf=' + str(DTreeSM.min_samples_leaf),
                          accuracy_score(y_test, predictedDTSM),
                          precision_score(y_test, predictedDTSM), 
                          recall_score(y_test, predictedDTSM), 
                          f1_score(y_test, predictedDTSM), t1]

df_results

# visualizing data 
dot_data = StringIO()

export_graphviz(
    DTreeSM,
    out_file = dot_data,
    feature_names=X_cols,
    rounded = True, 
    special_characters = True,
    filled = True)

graph = pydotplus.graph_from_dot_data(dot_data.getvalue())
Image(graph.create_png())

"""### Bit strange that Shift type is so predictive, could be because of correlations with nbump, nbump2.. or Leakage? Try and remove it from features."""

X_reduced_res = X_res.drop('shift_enc_0', axis=1)
X_reduced_test = X_test.drop('shift_enc_0', axis=1)

t0= time.perf_counter()

DTM_classifier = DecisionTreeClassifier(max_leaf_nodes=10)
DTreeSM = DTM_classifier.fit(X_reduced_res, y_res)

t1 = time.perf_counter() - t0

predictedDTSM = DTM_classifier.predict(X_reduced_test)

df_results['DecisionTree SMOTE-shift max_leaf_nodes=10'] = ['None', '',
                          'max_features=' + str(DTreeSM.max_features), 
                          'min_samples_split=' + str(DTreeSM.min_samples_split), 
                          'min_samples_leaf=' + str(DTreeSM.min_samples_leaf),
                          accuracy_score(y_test, predictedDTSM),
                          precision_score(y_test, predictedDTSM), 
                          recall_score(y_test, predictedDTSM), 
                          f1_score(y_test, predictedDTSM), t1]

df_results

X_red_cols = list(set(X_cols) - set(['shift_enc_0']))
# visualizing data 
dot_data = StringIO()

export_graphviz(
    DTreeSM,
    out_file = dot_data,
    feature_names=X_red_cols,
    rounded = True, 
    special_characters = True,
    filled = True)

graph = pydotplus.graph_from_dot_data(dot_data.getvalue())
Image(graph.create_png())

cm = confusion_matrix(y_test, predictedDTSM)
cm

precision_class1 = cm[1, 1] / (cm[0, 1] + cm[1, 1])
recall_class1 = cm[1, 1] / (cm[1, 0] + cm[1, 1])
precision_class1, recall_class1

precision_class2 = cm[0, 0] / (cm[1, 0] + cm[0, 0])
recall_class2 = cm[0, 0] / (cm[0, 1] + cm[0, 0])
precision_class2, recall_class2

df_results.T.to_csv(model_store_location + "FinalResults.csv")

"""# AdaBoost"""

from sklearn.model_selection import cross_val_score
from sklearn.ensemble import AdaBoostClassifier


clf = AdaBoostClassifier(n_estimators=100, learning_rate=1.0, random_state=0)
scores = cross_val_score(clf, X_train, y_train, scoring='f1', cv=skf, n_jobs=-1, error_score='raise')
scores.mean(), scores.std()

"""### AdaBoost Baseline"""

t0= time.perf_counter()

clf.fit(X_train, y_train)

t1 = time.perf_counter() - t0

y_pred = clf.predict(X_test)

print(confusion_matrix(y_test, y_pred))

df_results['AdaBoost Raw data'] = ['None', str(clf.algorithm),
                          'learning_rate=' + str(clf.learning_rate), 
                          'n_estimators=' + str(clf.n_estimators), 
                          '',
                          accuracy_score(y_test, y_pred),
                          precision_score(y_test, y_pred), 
                          recall_score(y_test, y_pred), 
                          f1_score(y_test, y_pred), t1]

df_results

t0= time.perf_counter()

clf.fit(X_res, y_res)

t1 = time.perf_counter() - t0

y_pred = clf.predict(X_test)

print(confusion_matrix(y_test, y_pred))

df_results['AdaBoost SMOTE data'] = ['None', str(clf.algorithm),
                          'learning_rate=' + str(clf.learning_rate), 
                          'n_estimators=' + str(clf.n_estimators), 
                          '',
                          accuracy_score(y_test, y_pred),
                          precision_score(y_test, y_pred), 
                          recall_score(y_test, y_pred), 
                          f1_score(y_test, y_pred), t1]

df_results

"""### AdaBoost SMOTE  DTD_Max_Depth=3"""

clf = AdaBoostClassifier(base_estimator=DecisionTreeClassifier(max_depth=3) , 
                         n_estimators=100, learning_rate=0.1, random_state=0)


t0= time.perf_counter()

clf.fit(X_res, y_res)

t1 = time.perf_counter() - t0

y_pred = clf.predict(X_test)

print(confusion_matrix(y_test, y_pred))

df_results['AdaBoost SMOTE DTD_Max_Depth=3'] = ['None', str(clf.algorithm),
                          'learning_rate=' + str(clf.learning_rate), 
                          'n_estimators=' + str(clf.n_estimators), 
                          '',
                          accuracy_score(y_test, y_pred),
                          precision_score(y_test, y_pred), 
                          recall_score(y_test, y_pred), 
                          f1_score(y_test, y_pred), t1]

df_results

"""### AdaBoost SMOTE Finetune"""

# define the model with default hyperparameters
model = AdaBoostClassifier()

# define the grid of values to search
grid = dict()
grid['n_estimators'] = [10, 50, 100, 500]
grid['learning_rate'] = [0.0001, 0.001, 0.01, 0.1, 1.0]
skf = StratifiedKFold(n_splits=5)

t0= time.perf_counter()

# define the grid search procedure
grid_search = GridSearchCV(estimator=model, param_grid=grid, n_jobs=-1, cv=skf, scoring='recall')
# execute the grid search
grid_result = grid_search.fit(X_res, y_res)

t1 = time.perf_counter() - t0


# summarize the best score and configuration
print("Best: %f using %s" % (grid_result.best_score_, grid_result.best_params_))
# summarize all scores that were evaluated
means = grid_result.cv_results_['mean_test_score']
stds = grid_result.cv_results_['std_test_score']
params = grid_result.cv_results_['params']
for mean, stdev, param in zip(means, stds, params):
    print("%f (%f) with: %r" % (mean, stdev, param))

y_pred = grid_result.best_estimator_.predict(X_res)
print(confusion_matrix(y_res, y_pred))
print(recall_score(y_res, y_pred, pos_label=1))

clf = AdaBoostClassifier(n_estimators=500, learning_rate=1.0)

clf.fit(X_res, y_res)

y_pred = clf.predict(X_test)

print(confusion_matrix(y_test, y_pred))

df_results['AdaBoost SMOTE finetuned'] = ['f1', str(clf.algorithm),
                          'learning_rate=' + str(clf.learning_rate), 
                          'n_estimators=' + str(clf.n_estimators), 
                          '',
                          accuracy_score(y_test, y_pred),
                          precision_score(y_test, y_pred), 
                          recall_score(y_test, y_pred), 
                          f1_score(y_test, y_pred), t1]

df_results

"""# XGBoost
## XGBoost model with raw data
"""

xgb_model = XGBClassifier(n_estimators=1000, learning_rate=0.05)
t0= time.perf_counter()

xgb_model.fit(X_train, y_train, 
             early_stopping_rounds=5, # stop overfitting
             eval_set=[(X_test, y_test)], 
             verbose=False)

t1 = time.perf_counter() - t0

y_pred = xgb_model.predict(X_test)

print(confusion_matrix(y_test, y_pred))

df_results['XGBoost'] = ['[Baseline]', 'XGBoost', '', '', '',
                          accuracy_score(y_test, y_pred), 
                          precision_score(y_test, y_pred), 
                          recall_score(y_test, y_pred), 
                          f1_score(y_test, y_pred), t1]

df_results

"""## XGBoost model with SMOTE data"""

xgb_model = XGBClassifier(n_estimators=1000, learning_rate=0.05)

# split SMOTE Data for Training and Eval set
train_X_res, val_X_res, train_y_res, val_y_res = train_test_split(X_res, y_res, random_state = 0)
t0= time.perf_counter()

xgb_model.fit(train_X_res, train_y_res, 
             early_stopping_rounds=5, # stop overfitting
             eval_set=[(val_X_res, val_y_res)], 
             verbose=False)

t1 = time.perf_counter() - t0


y_pred = xgb_model.predict(X_test)

print(confusion_matrix(y_test, y_pred))

df_results['XGBoost SMOTE'] = ['[Baseline]', 'XGBoost SMOTE', '', '', '',
                          accuracy_score(y_test, y_pred), 
                          precision_score(y_test, y_pred), 
                          recall_score(y_test, y_pred), 
                          f1_score(y_test, y_pred), t1]

df_results

16 / (16 +18)

from sklearn.feature_selection import SelectKBest, chi2

# Define our search space for grid search
search_space = [
  {
    'clf__n_estimators': [50, 100, 150, 200],
    'clf__learning_rate': [0.01, 0.1, 0.2, 0.3],
    'clf__max_depth': range(3, 10),
    'clf__colsample_bytree': [i/10.0 for i in range(1, 3)],
    'clf__gamma': [i/10.0 for i in range(3)],
    'fs__score_func': [chi2],
    'fs__k': [10],
  }
]

small_earch_space = [
  {
    'clf__n_estimators': [50, 100],
    'clf__learning_rate': [0.01, 0.1],
    'clf__max_depth': range(3, 10),
    'clf__colsample_bytree': [i/10.0 for i in range(1, 2)],
    'clf__gamma': [i/10.0 for i in range(2)],
    'fs__score_func': [chi2],
    'fs__k': [10],
  }
]
# Define cross validation
kfold = KFold(n_splits=10, random_state=42)

# Define grid search
grid = GridSearchCV(
  xgb_model,
  param_grid=search_space,
  cv=kfold,
  scoring=['recall'],
  refit='recall',
  verbose=1,
  n_jobs=-1
)
t0= time.perf_counter()

# Fit grid search
#grid_result = grid.fit(X_res, y_res)

t1 = time.perf_counter() - t0

#grid_result.best_estimator_

"""### LONG Optimisation, result here:
```
Fitting 10 folds for each of 672 candidates, totalling 6720 fits
/usr/local/lib/python3.6/dist-packages/sklearn/model_selection/_split.py:296: FutureWarning: Setting a random_state has no effect since shuffle is False. This will raise an error in 0.24. You should leave random_state to its default (None), or set shuffle=True.
  FutureWarning
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 2 concurrent workers.
[Parallel(n_jobs=-1)]: Done  46 tasks      | elapsed:  1.6min
[Parallel(n_jobs=-1)]: Done 196 tasks      | elapsed:  6.6min
[Parallel(n_jobs=-1)]: Done 446 tasks      | elapsed: 14.9min
[Parallel(n_jobs=-1)]: Done 796 tasks      | elapsed: 26.6min
[Parallel(n_jobs=-1)]: Done 1246 tasks      | elapsed: 41.8min
[Parallel(n_jobs=-1)]: Done 1796 tasks      | elapsed: 60.3min
[Parallel(n_jobs=-1)]: Done 2446 tasks      | elapsed: 82.0min
[Parallel(n_jobs=-1)]: Done 3196 tasks      | elapsed: 107.3min
[Parallel(n_jobs=-1)]: Done 4046 tasks      | elapsed: 135.7min
[Parallel(n_jobs=-1)]: Done 4996 tasks      | elapsed: 167.6min
[Parallel(n_jobs=-1)]: Done 6046 tasks      | elapsed: 202.9min
[Parallel(n_jobs=-1)]: Done 6720 out of 6720 | elapsed: 225.8min finished
XGBClassifier(base_score=0.5, booster='gbtree', clf__colsample_bytree=0.1,
              clf__gamma=0.0, clf__learning_rate=0.01, clf__max_depth=3,
              clf__n_estimators=50, colsample_bylevel=1, colsample_bynode=1,
              colsample_bytree=1, fs__k=10,
              fs__score_func=<function chi2 at 0x7faea0c6eb70>, gamma=0,
              learning_rate=0.05, max_delta_step=0, max_depth=3,
              min_child_weight=1, missing=None, n_estimators=1000, n_jobs=1,
              nthread=None, objective='binary:logistic', random_state=0,
              reg_alpha=0, reg_lambda=1, scale_pos_weight=1, seed=None,
              silent=None, subsample=1, verbosity=1)```
"""

xgb_model = XGBClassifier(base_score=0.5, booster='gbtree', clf__colsample_bytree=0.1,
              clf__gamma=0.0, clf__learning_rate=0.01, clf__max_depth=3,
              clf__n_estimators=50, colsample_bylevel=1, colsample_bynode=1,
              colsample_bytree=1, fs__k=10,
              fs__score_func=chi2, gamma=0,
              learning_rate=0.05, max_delta_step=0, max_depth=3,
              min_child_weight=1, missing=None, n_estimators=1000, n_jobs=1,
              nthread=None, objective='binary:logistic', random_state=0,
              reg_alpha=0, reg_lambda=1, scale_pos_weight=1, seed=None,
              silent=None, subsample=1, verbosity=1)


xgb_model.fit(train_X_res, train_y_res, 
             early_stopping_rounds=5, # stop overfitting
             eval_set=[(val_X_res, val_y_res)], 
             verbose=False)

y_pred = xgb_model.predict(X_test)

print(confusion_matrix(y_test, y_pred))

df_results['XGBoost SMOTE recall'] = ['[recall]', 'XGBoost SMOTE', '', '', '',
                          accuracy_score(y_test, y_pred), 
                          precision_score(y_test, y_pred), 
                          recall_score(y_test, y_pred), 
                          f1_score(y_test, y_pred), (225*60)]

df_results

"""## Features importance: compare XGBoost and RFC"""

xgb_model.feature_importances_

xgb_fi = pd.DataFrame({'cols':X_cols, 'imp':xgb_model.feature_importances_}
                       ).sort_values('imp', ascending=False)

xgb_fi["model"]="XGBoost"                       
xgb_fi

rfc_fi["model"]="Random Forest"

fi = rfc_fi.append(xgb_fi)
fi

f, ax = plt.subplots(figsize=(18, 5))

chart = sns.barplot(x="cols", y="imp", data=(fi.loc[fi.imp > 0.01]), hue="model", palette="Purples")

# Add a legend and informative axis label
#ax.legend(ncol=2, loc="best", frameon=True)
ax.set(xlim=(-1, 12), ylabel="Importance",
       xlabel="Features")
sns.despine(left=True, bottom=True)
chart.set_xticklabels(chart.get_xticklabels(), rotation=45, horizontalalignment='right')

# You can specify a rotation for the tick labels in degrees or with keywords.
#plt.xticks(fi["cols"], rotation='vertical')
# Pad margins so that markers don't get clipped by the axes
plt.margins(0.2)
# Tweak spacing to prevent clipping of tick-labels
plt.subplots_adjust(bottom=0.15)
plt.show()

"""# Final results"""

df_results.T.to_csv(model_store_location + "FinalResults.csv")

df_results