# Load libraries
import sys
print("Python version:\t\t{}".format(sys.version))
import pandas as pd
print("pandas version:\t\t{}".format(pd.__version__))
import matplotlib
print("matplotlib version:\t{}".format(matplotlib.__version__))
import numpy as np
print("numpy version:\t\t{}".format(np.__version__))
import scipy as sp
print("scipy version:\t\t{}".format(sp.__version__))
import sklearn
print("sklearn version:\t{}".format(sklearn.__version__))

# Visualization
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(style='whitegrid', context='notebook', palette='deep')
%matplotlib inline

# Models
from sklearn.model_selection import cross_val_score,GridSearchCV, StratifiedKFold, train_test_split
from sklearn.preprocessing import StandardScaler, LabelEncoder, OneHotEncoder

from sklearn.tree import DecisionTreeClassifier, DecisionTreeRegressor, ExtraTreeClassifier
from sklearn.naive_bayes import GaussianNB, BernoulliNB
from sklearn.tree import DecisionTreeClassifier
from sklearn.svm import SVC, NuSVC, LinearSVC
from sklearn.linear_model import LogisticRegression, PassiveAggressiveClassifier, RidgeClassifier, SGDClassifier, Perceptron
from sklearn.neighbors import KNeighborsClassifier
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis, QuadraticDiscriminantAnalysis
from sklearn.neural_network import MLPClassifier
from sklearn.ensemble import AdaBoostClassifier, VotingClassifier, ExtraTreesClassifier, RandomForestClassifier, GradientBoostingClassifier, BaggingClassifier
from sklearn.gaussian_process import GaussianProcessClassifier
from sklearn.decomposition import PCA

from xgboost import XGBClassifier

# Tools
from sklearn.metrics import mean_squared_error, roc_auc_score, roc_curve, auc
from sklearn.metrics import accuracy_score,classification_report, precision_recall_curve, confusion_matrix
from sklearn.model_selection import cross_val_score, GridSearchCV, learning_curve

from collections import Counter
from tqdm import tqdm
import string


# Do not show warnings (added after Notebook was finalized)
import warnings
warnings.filterwarnings('ignore')

# Import training and testing csv datasets
train = pd.read_csv('./Data/train.csv')
test = pd.read_csv('./Data/test.csv')

# Inspect data
train.head()

# Inspect data
test.head()

# tools functions: concat and divide
def concat_df(train, test):
    return pd.concat([train, test], sort=True).reset_index(drop=True)


def divide_df(full_df):
    return full_df.loc[:890], full_df.loc[891:].drop(['Survived'], axis=1)


df_all = concat_df(train, test)
# create a clone to test different transformations
df_exp = df_all.copy()

# Remove the Passengerid from the set as it does not need to be included in the models
passengerId = test.PassengerId

train = train.drop(['PassengerId'],axis=1)
test = test.drop(['PassengerId'],axis=1)

# list information for each feature (type, number of nun-null records)
train.info()

# list information for each feature (type, number of nun-null records)
test.info()

# compute percentage of missing values

# compute number of missing records
missing_total = train.isnull().sum().sort_values(ascending=False)

# convert to percentages
missing_percentage = missing_total/train.shape[0]*100

# display missing record %
print('Missing values in training set:')
pd.concat([missing_total,missing_percentage],keys=['Count','Percentage'],axis=1)

# compute percentage of missing values

# compute number of missing records
missing_total = test.isnull().sum().sort_values(ascending=False)

# convert to percentages
missing_percentage = missing_total/train.shape[0]*100

# display missing record %
print('Missing values in test set:')
pd.concat([missing_total,missing_percentage],keys=['Count','Percentage'],axis=1)

train[['Sex','Survived']].groupby(['Sex'],as_index=False).agg(['mean','count'])

train[['Pclass','Survived']].groupby(['Pclass'],as_index=False).agg(['mean','count'])

train[['SibSp','Survived']].groupby(['SibSp'],as_index=False).agg(['mean','count'])

train[['Parch','Survived']].groupby(['Parch'],as_index=False).agg(['mean','count'])

# compute all correlations
df_all_corr = df_all.corr().abs()

# display correlations with age
print("Correlation with Age feature")
df_all_corr['Age'].sort_values(ascending=False)[1:]

df_exp['Cat'] = df_exp['Sex'].astype(str) + df_exp['Pclass'].astype(str)
df_exp.groupby(['Cat']).median()['Age']

# fill missing age with groups Sex + PClass
df_all['Age'] = df_all.groupby(['Sex', 'Pclass'
                               ])['Age'].apply(lambda x: x.fillna(x.median()))

assert df_all['Age'].isnull().sum()==0

# Distribution of the data
train['Embarked'].value_counts(dropna=False)

train[train['Embarked'].isnull()]

df_all[(df_all['Pclass'] == 1) & (df_all['Sex'] == 'female') &
       (df_all['Fare'] >= 50) & (df_all['Fare'] <= 100)].groupby(
           ['Embarked']).count()['Age']

# Fill missing values
df_all['Embarked'].fillna('S',inplace=True)

# Display record corresponding to missing value
df_all[df_all.Fare.isnull()]

# Extract median
subset_med = df_all[(df_all['Pclass'] == 3) & (df_all['SibSp'] == 0) &
                  (df_all['Parch'] == 0)]['Fare'].median()

# Replace missing value
df_all['Fare'] = df_all['Fare'].fillna(subset_med)

# Feature inspection
df_all['Cabin'].sort_values().head(5)

# Extract deck name from cabine feature (replace with "M" (missing) if Cabin is null)
df_all['Deck'] = df_all['Cabin'].apply(lambda x: x[0] if pd.notnull(x) else "M")

df_plot = df_all.groupby(['Deck', 'Pclass']).size().reset_index().pivot(columns='Pclass', index='Deck', values=0)
df_plot.fillna(0)

df_plot = df_plot.div(df_plot.sum(axis=1), axis=0)

fig, ax = plt.subplots(figsize=(16,6))
#ax = sns.countplot(x='Deck', data=df_all, hue='Pclass')
df_plot.plot(kind='bar', stacked=True, ax = ax)
ax.set_title('Class Distribution per Decks', fontsize=15)
ax.set_ylabel("Fraction of deck")
ax.legend()
ax.set_ylim(0,1);

# move T Deck passenger to Deck A
idx = df_all[df_all['Deck'] == 'T'].index
df_all.loc[idx, 'Deck'] = 'A'

df_plot, _ = divide_df(df_all)
df_plot = df_plot.groupby(['Deck', 'Survived']).size().reset_index().pivot(
    columns='Survived', index='Deck').fillna(0).astype(int)

df_plot = df_plot.div(df_plot.sum(axis=1), axis=0)

fig, ax = plt.subplots(figsize=(16,6))
#ax = sns.countplot(x='Deck', data=df_all, hue='Pclass')
df_plot.plot(kind='bar', stacked=True, ax = ax, color=['r', 'b'])
ax.set_title('Survival Rate per Decks', fontsize=15)
ax.set_ylabel("Survival Rate")
ax.legend()
ax.set_ylim(0,1);

df_all['Deck'] = df_all['Deck'].replace(['A', 'B', 'C'], 'ABC')
df_all['Deck'] = df_all['Deck'].replace(['D', 'E'], 'DE')
df_all['Deck'] = df_all['Deck'].replace(['F', 'G'], 'FG')

# drop cabin
df_all = df_all.drop(['Cabin'], axis=1)

train, test = divide_df(df_all)

# Plot the box plot
pal = {'female':"salmon",'male':"skyblue"}
plt.subplots(figsize = (6,6))
ax = sns.barplot(x = "Sex", y = "Survived", data=train, palette = pal)
plt.title("Impact of gender on the survival rate",fontsize=15)
plt.ylabel("Survival rate",fontsize=15)
plt.xlabel("Sex",fontsize=15);

print("Male mean survival rate = \t {:.2f}%".format(train[train['Sex'] =='male']['Survived'].mean()*100))
print("Female mean survival rate = \t {:.2f}%".format(train[train['Sex'] != 'male']['Survived'].mean()*100))

# Plot kernel density plot
fig, ax = plt.subplots(figsize=(16,6))
ax = sns.kdeplot(train.loc[(train['Survived']==1),'Age'],shade=True,color='blue',label='Survived');
ax = sns.kdeplot(train.loc[(train['Survived']==0),'Age'],shade=True,color='red',label='Not Survived')
ax.set_ylabel('Frequency',fontsize=15)
ax.set_xlabel('Age',fontsize=15);

print("Toddler survival rate = \t{:.2f}%".format(train[train['Age'] <= 2]['Survived'].mean()*100,2))
print("Children survival rate = \t{:.2f}%".format(train[(train['Age'] > 2) & (train['Age'] <= 12)]['Survived'].mean()*100,2))
print("Teenager survival rate = \t{:.2f}%".format(train[(train['Age'] > 12) & (train['Age'] <= 18)]['Survived'].mean()*100,2))
print("Young Adult survival rate = \t{:.2f}%".format(train[(train['Age'] > 18) & (train['Age'] <= 34)]['Survived'].mean()*100,2))
print("Adult survival rate = \t\t{:.2f}%".format(train[(train['Age'] > 34) & (train['Age'] <= 50)]['Survived'].mean()*100,2))
print("pre-senior survival rate = \t{:.2f}%".format(train[(train['Age'] > 50) & (train['Age'] <= 70)]['Survived'].mean()*100,2))
print("senior survival rate = \t\t{:.2f}%".format(train[train['Age'] > 70]['Survived'].mean()*100,2))

# Plot the survival rate per class
pal = {1:"gold",2:"silver",3:'sandybrown'}
plt.subplots(figsize = (6,6))
ax = sns.barplot(x = "Pclass", y = "Survived", data=train, palette = pal)
plt.title("Impact of class on the survival rate",fontsize=15)
plt.ylabel("Survival rate",fontsize=15)
plt.xlabel("Class",fontsize=15);

print("Upper class survival rate =\t{:.2f}%".format(train[train['Pclass'] == 1]['Survived'].mean()*100))
print("Middle class survival rate =\t{:.2f}%".format(train[train['Pclass'] == 2]['Survived'].mean()*100))
print("Lower class survival rate =\t{:.2f}%".format(train[train['Pclass'] == 3]['Survived'].mean()*100))

# Plot the survival rate per class
pal = {1:"gold",2:"silver",3:'sandybrown'}
plt.subplots(figsize = (6,6))
ax = sns.barplot(x = "Pclass", y = "Fare", data=train, palette = pal)
plt.title("Impact of the passenger class on the fare",fontsize=15)
plt.ylabel("Fare",fontsize=15)
plt.xlabel("Pclass",fontsize=15);

print('Median fare:')
print("1st class =\t$${:.2f}".format(train[train['Pclass'] == 1]['Fare'].median()))
print("2nd class =\t$${:.2f}".format(train[train['Pclass'] == 2]['Fare'].median()))
print("3nd class =\t$${:.2f}".format(train[train['Pclass'] == 3]['Fare'].median()))

# Box plot
fig, ax = plt.subplots(figsize=(16,6))
sns.boxplot(x = 'Pclass', y = 'Fare', hue = 'Survived', data = train, ax = ax)
ax.set_title('Pclass vs Fare Survival Comparison',fontsize=15)
ax.set_ylim(0, 300);

fig = plt.figure(figsize=(16,6),)
ax=sns.distplot(train['Fare'] , color='blue',kde=False,bins=20)
ax.set_ylabel("Count",fontsize=15)
ax.set_title('Fare distribution',fontsize=15)
ax.set_xticks(range(0,600,25));

# Kernel Density Plot
fig = plt.figure(figsize=(16,6))
ax=sns.distplot(train.loc[(train['Survived'] == 0),'Fare'] , color='red',label='Not Survived',kde=False,bins=21)
ax=sns.distplot(train.loc[(train['Survived'] == 1),'Fare'] , color='blue', label='Survived',kde=False,bins=41)
ax.legend(['Not Survived','Survived'])
plt.xlabel("Fare",fontsize=15)
plt.ylabel("Frequency of Passenger Survived",fontsize=15)
plt.title('Fare Distribution Survived vs Non Survived',fontsize=15);

train[['SibSp','Survived']].groupby(['SibSp'],as_index=False).mean().sort_values(by='Survived',ascending =False)

fig, ax = plt.subplots(figsize=(16,6))
sns.barplot(x="SibSp", y="Survived",color='r',data=train);

train[['Parch','Survived']].groupby(['Parch'],as_index=False).mean().sort_values(by='Survived',ascending =False)

fig, ax = plt.subplots(figsize=(16,6))
sns.barplot(x="Parch", y="Survived",color='b', data=train);

train[['Embarked','Survived']].groupby(['Embarked'],as_index=False).mean().sort_values(by='Survived',ascending =False)

# Turning the Sex feature into a boolean classifier
train['Sex'] = train['Sex'].apply(lambda x: 0 if x == "female" else 1)
test['Sex'] = test['Sex'].apply(lambda x: 0 if x == "female" else 1)

train.describe()

train[['Sex', 'Survived']].groupby("Sex").mean()

train[['Pclass', 'Survived']].groupby("Pclass").mean()

# Feature correlation
train.corr()['Survived'].sort_values()

# Compute the correlation matrix
corr = train.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=(10, 12))

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

# Draw the heatmap with the mask and correct aspect ratio
ax = sns.heatmap(corr, mask=mask, cmap=cmap, vmax=1.0,vmin=-1.0, center=0,annot=True,
            square=True, linewidths=.5, cbar_kws={"shrink": .5})

ylabels = [x.get_text() for x in ax.get_yticklabels()]
plt.yticks(np.arange(len(ylabels))+0.5,ylabels, rotation=0, fontsize="10", va="center")
ax.set_ylim(8,0.5);

# extrace new feature using regular expression
df_all['Title'] = df_all['Name'].str.extract(r' ([A-Za-z]+)\.',expand=False)
df_all['Is_Married'] = 0
df_all['Is_Married'].loc[df_all['Title']=='Mrs'] = 1

fig, ax = plt.subplots(figsize=(16,6))
sns.barplot(x=df_all['Title'].value_counts().index, y=df_all['Title'].value_counts().values, ax=ax)
ax.set_title("Title Distribution");

df_plot, _ = divide_df(df_all)
df_plot = df_plot.groupby(['Title', 'Survived']).size().reset_index().pivot(columns='Survived', index='Title').fillna(0).astype(int)

df_plot = df_plot.div(df_plot.sum(axis=1), axis=0)

fig, ax = plt.subplots(figsize=(16,6))
#ax = sns.countplot(x='Deck', data=df_all, hue='Pclass')
df_plot.plot(kind='bar', stacked=True, ax = ax, color=['r', 'b'])
ax.set_title('Survival Rate per Decks', fontsize=15)
ax.set_ylabel("Survival Rate")
ax.legend(['Not Survived', 'Survived'], loc="upper right")
ax.set_ylim(0,1);

df_all['Title'] = df_all['Title'].replace(
    ['Miss', 'Mrs', 'Ms', 'Mlle', 'Lady', 'Mme', 'Countess', 'Dona'],'Miss/Mrs/Ms')
df_all['Title'] = df_all['Title'].replace(
    ['Dr', 'Col', 'Major', 'Jonkheer', 'Capt', 'Sir', 'Don', 'Rev'], 'Rare')

df_plot, _ = divide_df(df_all)
df_plot = df_plot.groupby(['Title', 'Survived']).size().reset_index().pivot(
    columns='Survived', index='Title').fillna(0).astype(int)

df_plot = df_plot.div(df_plot.sum(axis=1), axis=0)

fig, axes = plt.subplots(1, 2, figsize=(16,6))
df_plot.plot(kind='bar', stacked=True, ax = axes[0], color=['r', 'b'])
axes[0].set_title('Survival Rate per Decks', fontsize=15)
axes[0].set_ylabel("Survival Rate")
axes[0].legend(['Not Survived', 'Survived'], loc="upper right")
axes[0].set_ylim(0,1)

sns.barplot(x=df_all['Title'].value_counts().index, y=df_all['Title'].value_counts().values, ax=axes[1])
axes[1].set_title("Title Distribution");

df_all['Family_Size'] = df_all['Parch'] + df_all['SibSp'] + 1

fig, axes = plt.subplots(1, 2, figsize=(16, 6))

sns.barplot(x=df_all['Family_Size'].value_counts().index,
            y=df_all['Family_Size'].value_counts(),
            ax=axes[0])
axes[0].set_title("Family Size Distribution", fontsize=15)
sns.countplot(x='Family_Size',
              hue='Survived',
              data=df_all,
              ax=axes[1],
              palette=sns.diverging_palette(10, 255, sep=80, n=2))
axes[1].legend(['Not Survived', 'Survived'], loc='upper right')
axes[1].set_title("Survival Rate based on Family Size", fontsize=15)

# group family sizes
family_map = {
    1: 'Alone',
    2: 'Small',
    3: 'Small',
    4: 'Small',
    5: 'Medium',
    6: 'Medium',
    7: 'Large',
    8: 'Large',
    11: 'Large'
}

df_all['Family_Size_Grp'] = df_all['Family_Size'].map(family_map)

fig, axes = plt.subplots(1, 2, figsize=(16, 6))

sns.barplot(x=df_all['Family_Size_Grp'].value_counts().index,
            y=df_all['Family_Size_Grp'].value_counts(),
            ax=axes[0])
axes[0].set_title("Family Size Distribution", fontsize=15)
sns.countplot(x='Family_Size_Grp', hue='Survived', data=df_all, ax=axes[1], palette=sns.diverging_palette(10, 255, sep=80, n=2))
axes[1].legend(['Not Survived', 'Survived'], loc='upper right')
axes[1].set_title("Survival Rate based on Family Size", fontsize=15);

df_all.groupby('Ticket').size().sort_values(ascending=False)

df_all['Ticket_Frequency'] = df_all.groupby('Ticket')['Ticket'].transform('count')

fig, axes = plt.subplots(1, 2, figsize=(16, 6))

sns.barplot(x=df_all['Ticket_Frequency'].value_counts().index,
            y=df_all['Ticket_Frequency'].value_counts(),
            ax=axes[0])
axes[0].set_title("Ticket Frequency Distribution", fontsize=15)
sns.countplot(x='Ticket_Frequency', hue='Survived', data=df_all, ax=axes[1], palette=sns.diverging_palette(10, 255, sep=80, n=2))
axes[1].legend(['Not Survived', 'Survived'], loc='upper right')
axes[1].set_title("Survival Rate based on  Ticket Frequency", fontsize=15);

df_all['Last_Name'] = df_all['Name'].str.extract(r'^([a-zA-Z\s\-\']+)\,*')

df_train, df_test = divide_df(df_all)

# 1. families and tickets occuring in both sets
non_unique_families = df_train.loc[df_train['Last_Name'].isin(df_test['Last_Name'].unique()), 'Last_Name']
non_unique_tickets = df_train.loc[df_train['Ticket'].isin(df_test['Ticket'].unique()), 'Ticket']

# 2. compute median survival rate and size for each ticket and family
df_family = df_train.groupby('Last_Name').median()[['Survived', 'Family_Size']]
df_ticket = df_train.groupby('Ticket').median()[['Survived', 'Ticket_Frequency']]

# 3. filter families with more that one member and present in both the train and test sets
df_family = df_family.loc[(df_family.index.isin(non_unique_families)) & (df_family['Family_Size']>1)][['Survived']]
df_ticket = df_ticket.loc[(df_ticket.index.isin(non_unique_tickets)) & (df_ticket['Ticket_Frequency']>1)][['Survived']]

# 3. compute mean survival rate
mean_survival_rate = np.mean(df_train['Survived'])

# 3.
# assign family survival rate to each passenger, same for tickets
df_train = df_train.merge(right=df_family, how='left', left_on='Last_Name',right_index=True, suffixes=('', '_y'))
df_train = df_train.rename(columns={"Survived_y": "Family_Survival"})
df_train = df_train.merge(right=df_ticket, how='left', left_on='Ticket',right_index=True, suffixes=('', '_y'))
df_train = df_train.rename(columns={"Survived_y": "Ticket_Survival"})

df_test = df_test.merge(right=df_family, how='left', left_on='Last_Name', right_index=True, suffixes=('', '_y'))
df_test = df_test.rename(columns={"Survived": "Family_Survival"})
df_test = df_test.merge(right=df_ticket, how='left', left_on='Ticket', right_index=True, suffixes=('', '_y'))
df_test = df_test.rename(columns={"Survived": "Ticket_Survival"})

# new feature to determine if a family has a family-based survival rate, same for tickets
df_train['Has_Family_Survival'] = (~df_train['Family_Survival'].isnull()).astype(int)
df_test['Has_Family_Survival'] = (~df_test['Family_Survival'].isnull()).astype(int)
df_train['Has_Ticket_Survival'] = (~df_train['Ticket_Survival'].isnull()).astype(int)
df_test['Has_Ticket_Survival'] = (~df_test['Ticket_Survival'].isnull()).astype(int)
     
# fill null with mean survival rate, same for tickets
df_train['Family_Survival'] = df_train['Family_Survival'].fillna(mean_survival_rate)
df_test['Family_Survival'] = df_test['Family_Survival'].fillna(mean_survival_rate)
df_train['Ticket_Survival'] = df_train['Ticket_Survival'].fillna(mean_survival_rate)
df_test['Ticket_Survival'] = df_test['Ticket_Survival'].fillna(mean_survival_rate)

# 4. Combine survival rates (Family and Ticket)
for df in [df_train, df_test]:
    df['Survival_Rate'] = (df['Ticket_Survival'] + df['Family_Survival']) / 2
    df['Has_Survival'] = (df['Has_Ticket_Survival'] + df['Has_Family_Survival']) / 2  

df_all = concat_df(df_train, df_test)

df_all['Fare_Bin'] = pd.qcut(df_all['Fare'], 13)

fig, axs = plt.subplots(figsize=(16, 6))
sns.countplot(x='Fare_Bin', hue='Survived', data=df_all, palette=sns.diverging_palette(10, 255, sep=80, n=2))

plt.xlabel('Fare Bins')
plt.ylabel('Passenger Count')
plt.legend(['Not Survived', 'Survived'], loc='upper right')
plt.title('Count of Survival in {} Feature'.format('Fare Bin'), size=15)
plt.tight_layout()

df_all['Age_Bin'] = pd.qcut(df_all['Age'], 10)

fig, axs = plt.subplots(figsize=(16, 6))
sns.countplot(x='Age_Bin', hue='Survived', data=df_all, palette=sns.diverging_palette(10, 255, sep=80, n=2))

plt.xlabel('Age Bins')
plt.ylabel('Passenger Count')
plt.legend(['Not Survived', 'Survived'], loc='upper right')
plt.title('Count of Survival in {} Feature'.format('Age Bin'), size=15)
plt.tight_layout()

df_train, df_test = divide_df(df_all)

# features to be encoded
encoding = ['Embarked', 'Sex', 'Deck', 'Family_Size_Grp', 'Age_Bin', 'Fare_Bin', 'Title']

for feature in encoding:
    encoder = LabelEncoder()
    df_train[feature] = encoder.fit_transform(df_train[feature])
    df_test[feature] = encoder.transform(df_test[feature])

dummies = ['Pclass', 'Sex', 'Deck', 'Embarked', 'Title', 'Family_Size_Grp']

df_train = pd.get_dummies(df_train, prefix_sep='_', columns=dummies, drop_first=False)
df_test = pd.get_dummies(df_test, prefix_sep='_', columns=dummies, drop_first=False)

df_all = concat_df(df_train, df_test)

drop_columns = ['Family_Size', 'Survived', 'Name', 'Parch', 'Last_Name', 'PassengerId', 'SibSp',
                'Ticket', 'Ticket_Survival', 'Family_Survival', 'Has_Family_Survival', 'Has_Ticket_Survival',
                'Age', 'Fare']

drop_columns_test = ['Family_Size', 'Name', 'Parch', 'Last_Name', 'PassengerId', 'SibSp',
                     'Ticket', 'Ticket_Survival', 'Family_Survival', 'Has_Family_Survival', 'Has_Ticket_Survival',
                     'Age', 'Fare']

df_all = df_all.drop(columns=drop_columns)
df_all.head()

df_all.columns

X_train = StandardScaler().fit_transform(df_train.drop(columns=drop_columns))
y_train = df_train['Survived'].values
X_test = StandardScaler().fit_transform(df_test.drop(columns=drop_columns_test))

print('X_train shape: {}'.format(X_train.shape))
print('y_train shape: {}'.format(y_train.shape))
print('X_test shape: {}'.format(X_test.shape))

# PCA
pca = PCA(random_state=42)
X_pca = pca.fit_transform(X_train)

fig, ax = plt.subplots(figsize=(12,6))
ax.plot(np.cumsum(pca.explained_variance_ratio_))
ax.set_xlabel("Component")
ax.set_ylabel("Percent of Explained Variance")
ax.set_title("Cumulative Variance", fontsize=15)
ax.set_ylim(0,1)
ax.set_xlim(0,pca.n_components_);

fig, ax = plt.subplots(figsize=(16,16))
cmap = sns.diverging_palette(220, 10, as_cmap=True)
cset = ax.imshow(
    pca.components_.T,
    cmap = cmap,
    vmin=-1,
    vmax=1)

ylabels = df_train.drop(columns=drop_columns).columns.values.tolist()
ax.set_yticks(np.arange(len(ylabels)))
ax.set_yticklabels(ylabels, rotation=0)
ax.set_ylim(-0.5,len(ylabels)-0.5)
ax.set_xticks(range(pca.n_components_))
ax.grid(False)
ax.set_xlabel("Principal Component", fontsize=15)
ax.set_ylabel("Contribution", fontsize=15)
ax.set_title("Contribution of Features to Components", fontsize=15)
fig.colorbar(cset, ax=ax, shrink=0.5)
plt.tight_layout();

# Cross validate model with Kfold stratified cross val
kfold = StratifiedKFold(n_splits=5, random_state=5, shuffle=True)

random_state = 42

# classifiers
classifiers_list = [
    #Ensemble Methods
    AdaBoostClassifier(AdaBoostClassifier(DecisionTreeClassifier(random_state=random_state),random_state=random_state)),
    BaggingClassifier(random_state=random_state),
    ExtraTreesClassifier(random_state=random_state),
    GradientBoostingClassifier(random_state=random_state),
    RandomForestClassifier(random_state=random_state),

    #Gaussian Processes
    GaussianProcessClassifier(random_state=random_state),
    
    #GLM
    LogisticRegression(random_state=random_state),
    PassiveAggressiveClassifier(random_state=random_state),
    RidgeClassifier(),
    SGDClassifier(random_state=random_state),
    Perceptron(random_state=random_state),
    MLPClassifier(random_state=random_state),
    
    #Navies Bayes
    BernoulliNB(),
    GaussianNB(),
    
    #Nearest Neighbor
    KNeighborsClassifier(),
    
    #SVM
    SVC(probability=True, random_state=random_state),
    NuSVC(probability=True, random_state=random_state),
    LinearSVC(random_state=random_state),
    
    #Trees    
    DecisionTreeClassifier(random_state=random_state),
    ExtraTreeClassifier(random_state=random_state),
    
    #Discriminant Analysis
    LinearDiscriminantAnalysis(),
    QuadraticDiscriminantAnalysis(),

    #xgboost: http://xgboost.readthedocs.io/en/latest/model.html
    XGBClassifier()    
    ]

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

# perform cross-validation
for clf in tqdm(classifiers_list):
    cv_results_list.append(cross_val_score(clf,
                                           X_train,
                                           y_train,
                                           scoring = "accuracy",
                                           cv = kfold,
                                           n_jobs=-1))
    clf.fit(X_train, y_train)
    pred = clf.predict(X_test)

# 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":[clf.__class__.__name__ for clf in classifiers_list]})                    

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

cv_res_df.set_index('Algorithm')

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

# store best models
best_models = []

def optimize_model(model, X_train, y_train, paramgrid, random_state=10, suffix='_best',
                  metric = 'accuracy', kfold=StratifiedKFold(n_splits=5, random_state=5, shuffle=True), verbose=1, other_args={}, print_best=True):
    
    # adjust base estimator parameters
    estimator = model(**other_args)
    
    # create k-fold-based grid-search
    gridsearch = GridSearchCV(estimator, param_grid=paramgrid, cv=kfold, scoring=metric, n_jobs=-1, verbose=verbose)
    
    # fit grid search
    gridsearch.fit(X_train, y_train)
    
    # store best model
    name = model.__name__ + suffix
    best_models.append((name, gridsearch.best_estimator_))
    
    # print (optional)
    if print_best:
        print(gridsearch)
    
    # display accuracy
    print('Best {} model archieves {:.2f}% {}'.format(model.__name__, 100 * gridsearch.best_score_, metric))

param_grid = {
    'max_depth': [3, 5, 6, 7],
    'learning_rate': [0.01, 0.025, 0.05, 0.075, 0,1, 0,15, 0.2],
    'n_estimators': [500, 1000, 1500],
}

optimize_model(XGBClassifier, X_train, y_train, paramgrid=param_grid)

param_grid = {
    "max_depth": [None],
    "max_features": [5, 7, 'auto'],
    "min_samples_split": [4, 5, 6],
    "min_samples_leaf": [4, 5, 6],
    "bootstrap": [False],
    "n_estimators": [1500, 2000],
    "criterion": ["gini"]
}

optimize_model(ExtraTreesClassifier, X_train, y_train, paramgrid=param_grid)

param_grid = {
    "max_depth": [5, 6, 7],
    "min_samples_split": [4, 5, 6],
    "min_samples_leaf": [4, 5, 6],
    "max_features": [5, 7, 'auto'],
    "n_estimators": [1750],
    "criterion": ["gini"]
}

optimize_model(RandomForestClassifier, X_train, y_train, paramgrid=param_grid)

param_grid = {
    'kernel': ['rbf'],
    'gamma': [0.0005, 0.0008, 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1],
    'C': [1, 10, 50, 100, 200, 250, 500]
}

optimize_model(SVC, X_train, y_train, paramgrid=param_grid, other_args={'probability':True})

param_grid = {"C": np.logspace(-3, 3, 40), "penalty": ["l1", "l2"]}

optimize_model(LogisticRegression, X_train, y_train, paramgrid=param_grid)

param_grid = {
    "learning_rate": [0.01, 0.025, 0.05, 0.075, 0.1, 0.15, 0.2],
    "subsample": [0.5, 0.62, 0.8, 0.85, 0.87, 0.9, 0.95, 1.0]
}

optimize_model(GradientBoostingClassifier,
               X_train,
               y_train,
               paramgrid=param_grid)

param_grid = {
    "n_neighbors": [3, 5, 7, 11, 13, 15],
    "weights": ['uniform', 'distance'],
    "metric": ['manhattan'],
    "algorithm":['auto']
}

optimize_model(KNeighborsClassifier,
               X_train,
               y_train,
               paramgrid=param_grid)

LDA = LinearDiscriminantAnalysis()
LDA.fit(X_train, y_train)
LDA_best = LDA
best_models.append(("LDA_best", LDA_best))

predictions = {}
N_fold = 5


def generalize_predictions(model,
                           X_train,
                           y_train,
                           X_test,
                           skf=StratifiedKFold(n_splits=N_fold,
                                               random_state=N_fold,
                                               shuffle=True)):

    # store predicted probas for each folds
    probas = pd.DataFrame(np.zeros((len(X_test), N_fold * 2)),
                          columns=[
                              'Fold_{}_Class_{}'.format(i, j)
                              for i in range(1, N_fold + 1) for j in range(2)
                          ],
                          index=passengerId)

    # for each fold, train then predict
    for fold_id, (train_idx,
                  val_idx) in tqdm(enumerate(skf.split(X_train, y_train), 1)):

        # Fitting the model
        model.fit(X_train[train_idx], y_train[train_idx])

        # X_test probabilities
        probas.loc[:, 'Fold_{}_Class_0'.format(fold_id)] = model.predict_proba(
            X_test)[:, 0]
        probas.loc[:, 'Fold_{}_Class_1'.format(fold_id)] = model.predict_proba(
            X_test)[:, 1]

    # save results
    return (probas)

predictions = {name: generalize_predictions(model, X_train, y_train, X_test) for name, model in best_models}

def make_predictions(model_df, mode="hard", threshold=0.5):
    
    predictions = pd.DataFrame(np.zeros((len(X_test), 3)),
                     columns=['1', '0', 'pred'])
    
    # isolate probabilities of class 1
    class_one = [col for col in model_df.columns if col.endswith('Class_1')]
    
    # compute average of class 1 probabilities
    predictions['1'] = model_df[class_one].sum(axis=1) / N_fold
    predictions['0'] = 1. - predictions['1']
    
    if mode=="hard":
        predictions['pred'] = (predictions['1'] >= threshold).astype(int)
    else:
        predictions['pred'] = predictions['1']
    
    model_results = pd.DataFrame(columns=['PassengerId', 'Survived'])
    model_results['PassengerId'] = df_test['PassengerId']
    model_results['Survived'] = predictions['pred']
    return model_results

results = [(key, make_predictions(model)) for key, model in predictions.items()]