# import libraries
import pandas as pd
import numpy as np
import sys
np.set_printoptions(threshold=sys.maxsize)

import matplotlib.pyplot as plt
import seaborn as sns
# sns.set_style("whitegrid")
sns.set()
mycols = ["#66c2ff", "#00cc99", "#85e085", "#ffd966", "#ffb3b3", "#dab3ff", "#c2c2d6"]
sns.set_palette(palette = mycols, n_colors = 7)

from scipy import stats
from scipy.stats import skew, norm
from scipy.special import boxcox1p
from scipy.stats.stats import pearsonr
from scipy.stats import ttest_ind

from sklearn.model_selection import GridSearchCV, ShuffleSplit

import xgboost as xgb

import math

import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
warnings.simplefilter(action='ignore', category=RuntimeWarning)
warnings.filterwarnings("ignore", category=DeprecationWarning)

from sklearn.base import TransformerMixin, BaseEstimator

decimals = 2
%matplotlib inline

# load data
train = pd.read_csv('data/train.csv')
test = pd.read_csv('data/test.csv')

# save Id
train_ID = train['Id']
test_ID = test['Id']

# drop iD
del train['Id']
del test['Id']

# investigate shapes
print('Train shape:', train.shape)
print('Test shape:', test.shape)

# list features
print('List of features:')
print(train.columns)

# explore data
print("Here are a few observations: ")
train.head()

# explore data
train.info()

# explore data
test.info()

# explore data
train.describe()

def t_test_feature(df, feature_name, class_feature_name, same_variance=False, sign_level=0.05, ax=None):
    '''
    Perform t-test hypothesis testing and return a heatmap of the results:
    
    Significance level (alpha) : 0.05
    Null hypothesis (Ho): The two classes of the target feature have the same mean for the selected feature.
    Alternative hypothesis (Ha): The two classes of the target feature do not have the same mean for the selected feature.
    
    Inputs:
    df--input dataframe (pandas.DataFrame)
    feature_name--name of the feature to be considered for the hypothesis testing (string)
    class_feature_name--name of the target feature (string). The records are divided into classes corresponding to the classes of the class_feature_name.
    same_variance--should the variances of the two populations be considered equal (boolean)
    sign_level--threshold used to reject the null hypothesis.
    ax--axes to be used for plot (matplotlib.axes)
    
    Outputs:
    pvalues--pvalues of the hypothesis testing.
    plot--seaborn heatmap of the results of the hypothesis testing.
    '''
    
    # isolate classes of Cover_Type:
    classes = df[class_feature_name].unique().tolist()
    classes.sort()
    
    # iterate over classes and perform T-test hypothesis testing for identical mean
    pvalues =  np.full([len(classes),len(classes)], np.nan)
    
    for class_i in classes:
        for class_j in classes:
            if class_i>class_j:
                
                # isolate records for class i
                records_i = df.loc[df[class_feature_name]==class_i,feature_name]
                
                # isolate records for class j
                records_j = df.loc[df[class_feature_name]==class_j,feature_name]
                
                # compute pvalue and store in array
                _, pvalue = ttest_ind(records_i, records_j, axis=0, equal_var=True, nan_policy='omit')
                pvalues[class_i-1,class_j-1] = pvalue
                
    # filter p-value (pvalue<sign_level)
    pvalues = pvalues>sign_level
    
    # generate a mask for the upper triangle
    mask = np.zeros_like(pvalues, dtype=np.bool)
    mask[np.triu_indices_from(mask)] = True
    
    # plot heatmat
    cmap = sns.diverging_palette(10, 220, n=2)
    
    # create plot if no plot exists
    if ax is None:
        fig, ax = plt.subplots(figsize=(9,8))
    
    # plot heatmap
    ax = sns.heatmap(pvalues, mask=mask, cmap=cmap, vmax=1.0, center=0.5,
                     square=True, linewidths=.5, cbar_kws={"shrink": .3})
    
    # set colorbar properties
    ax.collections[0].colorbar.set_ticks([0,0.25,0.5,0.75 ,1])
    ax.collections[0].colorbar.ax.set_yticklabels(['','p<α','','p>α',''])
    ax.collections[0].colorbar.ax.tick_params(labelsize=15)
    
    # set plot title and tick properties
    ax.set_title('T-test: Identical Mean (α=0.05)')
    ax.set_yticklabels(classes)
    ax.set_xticklabels(classes)
    
    return pvalues

def plot_3_plots(feature):
    '''
    Plot:
    1--boxplot of feature data per Cover Type
    2--histogram of feature data per Cover Type
    3--t-test hypothesis testing resuts   
    '''
    sns.set()
    sns.set_palette(palette = mycols, n_colors = 7)
    
    fig, axarr = plt.subplots(1,3,figsize =(16, 6))
    fig.suptitle(feature.upper(),fontsize=15)
    
    # boxplot
    sns.boxplot(x='Cover_Type', y=feature,data=train, ax=axarr[0])
    axarr[0].set_ylim(0,train[feature].max()*1.10)
    
    # barplot   
    sns.barplot(x='Cover_Type', y=feature,data=train, ax=axarr[1], ci='sd')
    axarr[1].set_ylim(axarr[0].get_ylim())
    axarr[1].set_yticklabels([])
    axarr[1].set_ylabel('')
    
    # ttest   
    t_test_feature(train, feature, 'Cover_Type', same_variance=False, sign_level=0.05, ax=axarr[2])
    
    return None;

def plot_continuous_features(df, feature_name, class_feature_name):
    '''
    Return multiple plots used to inspect a continuous feature. One stacked histogram and one kde plot are produced.
    
    Inputs:
    df--input dataframe (pandas.DataFrame)
    feature_name--name of the continuous feature to be inspected (string)
    class_feature_name--name of the class feature (string)
    
    Output:
    None
    '''
    sns.set()
    # isolate feature records for each class
    values_per_class = [df.loc[df[class_feature_name]==i,feature_name] for i in np.sort(df[class_feature_name].unique())]

    # create figure
    fig, axes = plt.subplots(2, 1, figsize=(16,12), sharex=True)
    axes[0].set_title(feature_name.upper(),fontsize=15)
    
    # plot histogram, assign legend and scale y axis
    (n, bins, patches) = axes[0].hist(values_per_class,stacked=True, bins=50, color=mycols[0:7])
    
    # add legend to plot
    axes[0].legend(['1 - Spruce/Fir','2 - Lodgepole Pine','3 - Ponderosa Pine',
                    '4 - Cottonwood/Willow','5 - Aspen','6 - Douglas-fir',
                    '7 - Krummholz']);

    # histogram upper value and round to upper 100
    y_max = n[-1].max() * 1.05
    y_max = int(math.ceil(y_max / 100.0)) * 100   
    axes[0].set_ylim((0,y_max))
    
    for ix, el in enumerate(values_per_class):
        sns.kdeplot(el,shade=True,
                    color=mycols[ix],
                    ax=axes[1])

    # add legend to plot
    labels = ['1 - Spruce/Fir','2 - Lodgepole Pine','3 - Ponderosa Pine',
                    '4 - Cottonwood/Willow','5 - Aspen','6 - Douglas-fir',
                    '7 - Krummholz']
    axes[1].legend(labels);
        
    # Initialize the FacetGrid object
    df = train[[feature_name,class_feature_name]]
    
    
    #sns.p(palette = mycols, n_colors = 7)
    sns.set(style="white", rc={"axes.facecolor": (0, 0, 0, 0)})
    #pal = sns.cubehelix_palette(10, rot=-.25, light=.7)
    g = sns.FacetGrid(df, row=class_feature_name, hue=class_feature_name, aspect=15, height=1, palette=mycols)
    
    # Draw the densities in a few steps
    g.map(sns.kdeplot, feature_name, clip_on=False, shade=True, alpha=1, lw=1.5, bw=.2)
    g.map(sns.kdeplot, feature_name, clip_on=False, color="w", lw=2, bw=.2)
    g.map(plt.axhline, y=0, lw=2, clip_on=False)

    g.map(label, feature_name)
    
    # Set the subplots to overlap
    g.fig.subplots_adjust(hspace=-0.50)
    # Remove axes details that don't play well with overlap
    g.set_titles("")
    g.set(yticks=[])
    g.despine(bottom=True, left=True)
    plt.show();
    return None;

# Define and use a simple function to label the plot in axes coordinates
def label(x, color, label):
    ax = plt.gca()
    ax.text(0, .2, label, fontweight="bold", color=color,
            ha="left", va="center", transform=ax.transAxes, fontsize=20)

def plot_multiple_boolean_subclasses(df, feature_root_name, class_feature_name):
    '''
    Plot the histogram of any boolean feature divided into subclasses.
    
    Inputs:
    df--input dataframe (pandas.DataFrame)
    feature_root_name--root of the feature divided into subclasses (string)
    class_feature_name--feature to predict (string)
    '''
    
    # create list of features based on root name
    features = [col for col in df.columns if feature_root_name in col]
    
    # group by and count number of positive record per subclass
    grouped_df = df.groupby(class_feature_name)[features].sum().T

    # define width of plot based on number of subclasses
    subclass_count = len(features)
    plot_width = min(max(8,subclass_count+4), min(16,subclass_count+4))

    # create figure and plot results
    fig,ax = plt.subplots(figsize=(plot_width, 6))
    grouped_df.plot.bar(stacked=True,color=mycols,ax=ax);
    
    return None;

# Plot sample distribution
fig = plt.figure(figsize=(8, 6))
train['Cover_Type'].value_counts().sort_index().plot(kind='bar',color=mycols)
plt.title('Sample distribution',fontsize=20)
plt.xlabel('Cover_Type',fontsize=15)
plt.ylabel('Count',fontsize=15);
print(train['Cover_Type'].value_counts().sort_index());

plot_continuous_features(train,'Elevation','Cover_Type')

plot_3_plots('Elevation')

plot_continuous_features(train,'Aspect','Cover_Type')

plot_3_plots('Aspect')

plot_continuous_features(train,'Slope','Cover_Type')

plot_3_plots('Slope')

plot_continuous_features(train,'Horizontal_Distance_To_Hydrology','Cover_Type')

plot_3_plots('Horizontal_Distance_To_Hydrology')

train['log_Horizontal_Distance_To_Hydrology'] = train['Horizontal_Distance_To_Hydrology'].apply(lambda x: np.log(x) if x>0 else 0)
plot_continuous_features(train,'log_Horizontal_Distance_To_Hydrology','Cover_Type')

plot_continuous_features(train,'Vertical_Distance_To_Hydrology','Cover_Type')

plot_3_plots('Vertical_Distance_To_Hydrology')

train['log_Vertical_Distance_To_Hydrology'] = train['Vertical_Distance_To_Hydrology'].apply(lambda x: np.log(x) if x>0 else 0)
plot_continuous_features(train,'log_Vertical_Distance_To_Hydrology','Cover_Type')

plot_continuous_features(train,'Horizontal_Distance_To_Roadways','Cover_Type')

plot_3_plots('Horizontal_Distance_To_Roadways')

plot_continuous_features(train,'Hillshade_9am','Cover_Type')

plot_3_plots('Hillshade_9am')

plot_continuous_features(train,'Hillshade_Noon','Cover_Type')

plot_3_plots('Hillshade_Noon')

plot_continuous_features(train,'Hillshade_3pm','Cover_Type')

plot_3_plots('Hillshade_3pm')

plot_continuous_features(train,'Horizontal_Distance_To_Fire_Points','Cover_Type')

plot_3_plots('Horizontal_Distance_To_Fire_Points')

plot_multiple_boolean_subclasses(train, 'Wilderness_Area', 'Cover_Type')

plot_multiple_boolean_subclasses(train, 'Soil_Type', 'Cover_Type')

# Data distribution
train.sum(axis=0).sort_values(ascending=True)[0:10]

# drop empty soil types
def drop_empty_soil(train, test):
    '''drop the Soil_Type7 and Soil_Type15 features of train and test DataFrame.'''
    train = train.drop(['Soil_Type7', 'Soil_Type15'], axis=1)
    test = test.drop(['Soil_Type7', 'Soil_Type15'], axis=1)
    
    return train, test

# drop empty features
train, test = drop_empty_soil(train, test)

# Compute the correlation matrix
corr = train[['Elevation', 'Aspect', 'Slope',
       'Horizontal_Distance_To_Hydrology', 'Vertical_Distance_To_Hydrology',
       'Horizontal_Distance_To_Roadways', 'Hillshade_9am', 'Hillshade_Noon',
       'Hillshade_3pm', 'Horizontal_Distance_To_Fire_Points']].corr()

# Generate a mask for the upper triangle
mask = np.zeros_like(corr, dtype=np.bool)
mask[np.triu_indices_from(mask)] = True

# Set up the matplotlib figure
f, ax = plt.subplots(figsize=(15, 13))

# Generate a custom diverging colormap
cmap = sns.diverging_palette(220, 10, as_cmap=True)

# Draw the heatmap with the mask and correct aspect ratio
g = sns.heatmap(corr, mask=mask, cmap=cmap, vmax=1.0,vmin=-1.0, center=0,annot=True,
            square=True, linewidths=.5, cbar_kws={"shrink": .3})
g.set_xticklabels(g.get_xticklabels(), rotation = 45, ha="right");

# load data
train = pd.read_csv('data/train.csv')
test = pd.read_csv('data/test.csv')

# save Id
train_ID = train['Id']
test_ID = test['Id']

# create X and y dataframe
X_train = train.drop(['Cover_Type'],axis=1).copy()
y_train = train['Cover_Type'].copy()
X_test = test.copy()

# drop iD
del X_train['Id']
del X_test['Id']

# isolate categorical and numerical feature
# different pipelines will be used
numerical_features = [
    'Elevation','Aspect','Slope','Horizontal_Distance_To_Hydrology',
    'Vertical_Distance_To_Hydrology','Horizontal_Distance_To_Roadways',
    'Hillshade_9am','Hillshade_Noon','Hillshade_3pm','Horizontal_Distance_To_Fire_Points']

catergorical_features = list(X_train.columns.difference(numerical_features))

new_numerical_features = [
    'log_Horizontal_Distance_To_Hydrology', 'log_Vertical_Distance_To_Hydrology',
    'Slope_sq', 'Elevation_sq', 'Aspect_sq', 'Horizontal_Distance_To_Hydrology_sq',
    'Vertical_Distance_To_Hydrology_sq','Horizontal_Distance_To_Roadways_sq',
    'Horizontal_Distance_To_Fire_Points_sq','Hillshade_mean','Euclidean_Distance_to_Hydrology',
    'Cosine_Slope','Sine_Slope','H_to_F_sum','H_to_F_sub','H_to_R_sum','H_to_R_sub','R_to_F_sum',
    'R_to_F_sub','Slope_Hydrology']

new_categorical_features = ['Has_Soil_Type']

full_feature_list_with = numerical_features+new_numerical_features+catergorical_features+new_categorical_features
full_feature_list_without = numerical_features+catergorical_features+new_categorical_features

# import scalers
from sklearn.preprocessing import MinMaxScaler, StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.pipeline import FeatureUnion

class DataFrameSelector(BaseException, TransformerMixin):
    
    def __init__(self, attributes_names):
        self.attributes_names = attributes_names
        
    def fit(self, X, y=None):
        return self
    
    def transform(self, X, y=None):
        return X[self.attributes_names]

class NewNumericalFeatureGenerator(BaseEstimator, TransformerMixin):
    '''
    Generate new features on input dataframe.
    '''
    def __init__(self, add_new_features=True):
        self.add_new_features = add_new_features 
    
    def fit(self, df, y=None):
        return self
    
    def transform(self, df, y=None):
        
        df_copy = df.copy()
        
        df_copy = df_copy.astype(np.float64)
        
        if self.add_new_features:
        
            # log features
            log_features = ['Horizontal_Distance_To_Hydrology', 'Vertical_Distance_To_Hydrology']
            
            for feature in log_features:
                df_copy['log_'+feature] = df_copy[feature].apply(lambda x: np.log(x) if x>0 else 0)
                
            # squared features
            squared_features = ['Slope', 'Elevation', 'Aspect', 'Horizontal_Distance_To_Hydrology',
                               'Vertical_Distance_To_Hydrology','Horizontal_Distance_To_Roadways',
                               'Horizontal_Distance_To_Fire_Points']
            
            for feature in squared_features:
                df_copy[feature+'_sq'] = df_copy[feature]**2
                
            # Hillshade_mean
            df_copy['Hillshade_mean'] = (df_copy['Hillshade_9am']+df_copy['Hillshade_Noon']+df_copy['Hillshade_3pm'])/3.
              
            # Euclidean_Distance_to_Hydrology
            df_copy['Euclidean_Distance_to_Hydrology'] = np.sqrt(df_copy['Horizontal_Distance_To_Hydrology_sq']+df_copy['Vertical_Distance_To_Hydrology_sq'])
                   
            # Cosine_Slope
            df_copy['Cosine_Slope'] = np.cos(df_copy['Slope'] * np.pi / 180)
            
            # Sine_Slope
            df_copy['Sine_Slope'] = np.sin(df_copy['Slope'] * np.pi / 180)
            
            # H_to_F_sum
            df_copy['H_to_F_sum'] = np.abs(df_copy['Horizontal_Distance_To_Hydrology'] + df_copy['Horizontal_Distance_To_Fire_Points'])
            
            # H_to_F_sub
            df_copy['H_to_F_sub'] = np.abs(df_copy['Horizontal_Distance_To_Hydrology'] - df_copy['Horizontal_Distance_To_Fire_Points'])
            
            # H_to_R_sum
            df_copy['H_to_R_sum'] = np.abs(df_copy['Horizontal_Distance_To_Hydrology'] + df_copy['Horizontal_Distance_To_Roadways'])
            
            # H_to_R_sub
            df_copy['H_to_R_sub'] = np.abs(df_copy['Horizontal_Distance_To_Hydrology'] - df_copy['Horizontal_Distance_To_Roadways'])
            
            # R_to_F_sum
            df_copy['R_to_F_sum'] = np.abs(df_copy['Horizontal_Distance_To_Roadways'] + df_copy['Horizontal_Distance_To_Fire_Points'])
            
            # R_to_F_sub
            df_copy['R_to_F_sub'] = np.abs(df_copy['Horizontal_Distance_To_Roadways'] - df_copy['Horizontal_Distance_To_Fire_Points'])
            
            # Slope_Hydrology
            df_copy['Slope_Hydrology'] = df_copy['Vertical_Distance_To_Hydrology'] / df_copy['Horizontal_Distance_To_Hydrology']
            df_copy['Slope_Hydrology'] = df_copy['Slope_Hydrology'].map(lambda x: 0 if np.isnan(x) else x)
        
        return df_copy

class NewCategoricalFeatureGenerator(BaseEstimator, TransformerMixin):
    '''
    Generate new features on input dataframe.
    '''
    def __init__(self):
        pass
    
    def fit(self, df, y=None):
        return self
    
    def transform(self, df, y=None):
        
        df_copy = df.copy()
            
        # Has_Soil_Type
        soil_types = [x for x in df_copy.columns if 'Soil_Type' in x]
        df_copy['Has_Soil_Type'] = df_copy[soil_types].sum(axis=1)
        
        # return numpy array
        df_copy_values = df_copy.values
        
        return df_copy_values

class SpecialScaler(BaseEstimator, TransformerMixin):
    '''
    Standardize features
    '''
    def __init__(self, scaling_type='MinMax'):
        
        if scaling_type=='MinMax':
            self.scaler = MinMaxScaler()
        
        elif scaling_type=='StandardScaler':
            self.scaler = StandardScaler()
            
        else:
            self.scaler = None
    
    def fit(self, df, y=None):
        if self.scaler is None:
            pass
        else:
            self.scaler.fit(df)
        return self
        
    def transform(self, df):
        if self.scaler is None:
            return df.values
        else:
            return self.scaler.transform(df)

num_pipeline_with_new_features = Pipeline([
    ('selector',DataFrameSelector(numerical_features)),
    ('feature_adder',NewNumericalFeatureGenerator(add_new_features=True)),
    ('scaler',SpecialScaler(scaling_type='StandardScaler'))
])

num_pipeline_without_new_features = Pipeline([
    ('selector',DataFrameSelector(numerical_features)),
    ('feature_adder',NewNumericalFeatureGenerator(add_new_features=False)),
    ('scaler',SpecialScaler(scaling_type='StandardScaler'))
])


cat_pipeline = Pipeline([
    ('selector',DataFrameSelector(catergorical_features)),
    ('feature_adder',NewCategoricalFeatureGenerator()),
])

full_pipeline_with_new_features = FeatureUnion(transformer_list=[('num_pipeline',num_pipeline_with_new_features),
                                                                 ('cat_pipeline',cat_pipeline)])

full_pipeline_without_new_features = FeatureUnion(transformer_list=[('num_pipeline',num_pipeline_without_new_features),
                                                                    ('cat_pipeline',cat_pipeline)])

X_train_prepared_with = full_pipeline_with_new_features.fit_transform(X_train)
X_train_prepared_without = full_pipeline_without_new_features.fit_transform(X_train)

from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import GridSearchCV

# create model
forest_clf = RandomForestClassifier(random_state=42)

# create gridsearch
param_grid = {'n_estimators':[3,10,30,50,100],'max_features':[2,4,6,8],'max_depth':[3,5,10,15]}

# train model
grid_search_with = GridSearchCV(forest_clf,param_grid,cv=5,scoring='accuracy',n_jobs=4,verbose=2)
grid_search_with.fit(X_train_prepared_with, y_train)

grid_search_without = GridSearchCV(forest_clf,param_grid,cv=5,scoring='accuracy',n_jobs=4,verbose=2)
grid_search_without.fit(X_train_prepared_without, y_train)

print("Best accuracy with new features",grid_search_with.best_score_)
print("Best accuracy without new features",grid_search_without.best_score_)

feature_importances_with = grid_search_with.best_estimator_.feature_importances_

for ix,item in enumerate(sorted(zip(feature_importances_with,full_feature_list_with),reverse=True)):
    print(ix+1, item)

feature_importances_without = grid_search_without.best_estimator_.feature_importances_

for ix,item in enumerate(sorted(zip(feature_importances_without,full_feature_list_without),reverse=True)):
    print(ix+1, item)

X_train_new_features = NewNumericalFeatureGenerator(add_new_features=True).fit_transform(X_train)

# Compute the correlation matrix
corr = X_train_new_features[new_numerical_features].corr()
corr = np.around(corr, decimals=2)

# Generate a mask for the upper triangle
mask = np.zeros_like(corr, dtype=np.bool)
mask[np.triu_indices_from(mask)] = True

# Set up the matplotlib figure
f, ax = plt.subplots(figsize=(16, 16))

# Generate a custom diverging colormap
cmap = sns.diverging_palette(220, 10, as_cmap=True)

# Draw the heatmap with the mask and correct aspect ratio
g = sns.heatmap(corr, mask=mask, cmap=cmap, vmax=1.0,vmin=-1.0, center=0,annot=True,
            square=True, linewidths=.5, cbar_kws={"shrink": .2})
g.set_xticklabels(g.get_xticklabels(), rotation = 40, fontsize = 9, ha='right')
g.set_yticklabels(g.get_yticklabels(), fontsize = 9);

# add feature selection to grid search
def indices_of_top_k(arr, k):
    return np.sort(np.argpartition(np.array(arr),-k)[-k:])

class TopFeatureSelector(BaseEstimator, TransformerMixin):
    '''Select the k best feature from dataset'''
    
    def __init__(self, feature_importances, k):
        self.feature_importances = feature_importances
        self.k = k
        
    def fit(self, X, y=None):
        self.feature_indices_ = indices_of_top_k(self.feature_importances, self.k)
        return self
    
    def transform(self, X):
        return X[:, self.feature_indices_]

# first example
k = 45

# pipeline
preparation_and_feature_selection_pipeline_with = Pipeline([
    ('preparation',full_pipeline_with_new_features),
    ('feature_selection',TopFeatureSelector(feature_importances_with,k))
])

preparation_and_feature_selection_pipeline_without = Pipeline([
    ('preparation',full_pipeline_without_new_features),
    ('feature_selection',TopFeatureSelector(feature_importances_without,k))
])

preparation_and_feature_selection_pipeline_with.fit(X_train)
preparation_and_feature_selection_pipeline_without.fit(X_train)

X_train_prepared_with = preparation_and_feature_selection_pipeline_with.transform(X_train)
X_train_prepared_without = preparation_and_feature_selection_pipeline_without.transform(X_train)

# import modules
# Models
from sklearn.tree import DecisionTreeClassifier, DecisionTreeRegressor
from sklearn.model_selection import cross_val_score,GridSearchCV, StratifiedKFold, train_test_split
from sklearn.naive_bayes import GaussianNB
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier, BaggingClassifier
from sklearn.svm import SVC
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression, SGDClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.neural_network import MLPClassifier
from sklearn.ensemble import AdaBoostClassifier, VotingClassifier, ExtraTreesClassifier

from xgboost import XGBClassifier

# cross validation
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import accuracy_score

# cross validate model with Kfold stratified cross val
kfold = StratifiedKFold(n_splits=10)

# random state used for all models 
random_state = 42

# store models in list
classifiers_list = []

classifiers_list.append(LogisticRegression(random_state=random_state))
classifiers_list.append(KNeighborsClassifier())
classifiers_list.append(SVC(random_state=random_state))
classifiers_list.append(GaussianNB())
classifiers_list.append(DecisionTreeClassifier(random_state=random_state))
classifiers_list.append(BaggingClassifier(random_state=random_state))
classifiers_list.append(RandomForestClassifier(random_state=random_state))
classifiers_list.append(GradientBoostingClassifier(random_state=random_state))
classifiers_list.append(LinearDiscriminantAnalysis())
classifiers_list.append(MLPClassifier(random_state=random_state))
classifiers_list.append(AdaBoostClassifier(DecisionTreeClassifier(random_state=random_state),random_state=random_state))
classifiers_list.append(ExtraTreesClassifier(random_state=random_state))
classifiers_list.append(SGDClassifier(random_state=random_state))
classifiers_list.append(XGBClassifier(random_state=random_state))

# store cv results in list
cv_results_list = []
cv_means_list = []
cv_std_list = []

# perform cross-validation
for clf in classifiers_list:
    cv_results_list.append(cross_val_score(clf,
                                           X_train_prepared_with,
                                           y_train,
                                           scoring = "accuracy",
                                           cv = kfold,
                                           n_jobs=4))
    
    cv_results_list.append(cross_val_score(clf,
                                           X_train_prepared_without,
                                           y_train,
                                           scoring = "accuracy",
                                           cv = kfold,
                                           n_jobs=4))
        
    print("Training {}: Done".format(type(clf).__name__))

# store mean and std accuracy
for cv_result in cv_results_list:
    cv_means_list.append(cv_result.mean())
    cv_std_list.append(cv_result.std())
                      
cv_res_df = pd.DataFrame({"CrossValMeans":cv_means_list,
                          "CrossValerrors": cv_std_list,
                          "Algorithm":["LogReg_with","LogReg","KNN_with","KNN","SVC_with","SVC",
                                       "GaussianNB_with","GaussianNB","DecisionTree_with","DecisionTree",
                                       "Bagging_with","Bagging","RdmForest_with","RdmForest",
                                       "GradientBoost_with","GradientBoost","LDA_with","LDA",
                                       "NN_with","NN","AdaBoost_with","AdaBoost",
                                       "ExtraTrees","ExtraTrees_with",'SGDC','SGDC_with','XGBC_with','XGBC']})                    

cv_res_df = cv_res_df.sort_values(by='CrossValMeans',ascending=False)

# plot results
fig, ax = plt.subplots(figsize=(16,16))
g = sns.barplot("CrossValMeans",
                "Algorithm",
                data = cv_res_df,
                palette="Set3",
                orient = "h", ax=ax,
                **{'xerr':cv_std_list})
g.set_xlabel("Mean Accuracy")
g = g.set_title("Cross validation scores")

from sklearn.model_selection import GridSearchCV

prepare_select_and_predict_pipeline = Pipeline([
    ('preparation',full_pipeline_with_new_features),
    ('feature_selection',TopFeatureSelector(feature_importances_with,k)),
    ('extra_tree',ExtraTreesClassifier(random_state=42))
])

extra_tree_clf = ExtraTreesClassifier(random_state=42)

param_grid = [{'extra_tree__n_estimators':[50,100,150],
               'extra_tree__criterion':['gini','entropy'],
               'extra_tree__max_depth':[3,5,10,15,20],
               'extra_tree__min_samples_split':[2,3,5],
               'extra_tree__min_samples_leaf':[2,3,5],
               'extra_tree__max_features':['auto',None],
               'feature_selection__feature_importances':[feature_importances_with],
               'feature_selection__k':[10, 20, 30, 45],
               'preparation__num_pipeline__feature_adder__add_new_features':[True]},
              {'extra_tree__n_estimators':[50,100,150],
               'extra_tree__criterion':['gini','entropy'],
               'extra_tree__max_depth':[3,5,10,15,20],
               'extra_tree__min_samples_split':[2,3,5],
               'extra_tree__min_samples_leaf':[2,3,5],
               'extra_tree__max_features':['auto',None],
               'feature_selection__feature_importances':[feature_importances_without],
               'feature_selection__k':[10, 20, 30, 45],
               'preparation__num_pipeline__feature_adder__add_new_features':[False]}
             ]

grid_search_pipeline = GridSearchCV(prepare_select_and_predict_pipeline,param_grid,cv=5,scoring='accuracy',verbose=1,n_jobs=4)
grid_search_pipeline.fit(X_train, y_train)

grid_search_pipeline.best_score_

grid_search_pipeline.best_params_

submission_pipeline = Pipeline([
    ('preparation',full_pipeline_with_new_features),
    ('feature_selection',TopFeatureSelector(feature_importances_with,30))
])

X_train_prep = submission_pipeline.fit_transform(X_train)

final_model = ExtraTreesClassifier(random_state=42, criterion='entropy',
max_depth=20,
max_features=None,
min_samples_leaf=2,
min_samples_split= 2,
n_estimators=100)

final_model.fit(X_train_prep,y_train)

X_test_prep = submission_pipeline.transform(X_test)

# make prediction on test set
y_test_pred = final_model.predict(X_test_prep).astype(int)

# generate submission file 
submission_df = pd.DataFrame({'Id': test_ID,
                            'Cover_Type': y_test_pred})
submission_df.to_csv("extra_tree_full_pipe.csv", index=False)