import pandas as  pd

df = pd.read_csv("data/telco_churn.csv")

df.shape

df["SeniorCitizen"] = df["SeniorCitizen"].astype("object")

pd.options.display.max_rows = None

df["TotalChargesConv"] = pd.to_numeric(df["TotalCharges"], errors='coerce')

df["TotalChargesConv"].isnull().sum()

df.head()

def missing_data(data):
    total = data.isnull().sum().sort_values(ascending = False)
    percent = (data.isnull().sum()/data.isnull().count()*100).sort_values(ascending = False)
    return pd.concat([total, percent], axis=1, keys=['Total', 'Percent'])

missing_data(df)

df["TotalChargesConv"] = df["TotalChargesConv"].fillna(df["TotalChargesConv"].mean())

del df["TotalCharges"]

import numpy as np
# Keeping 15% for the validation set, thats about 1000 instances to test on 
msk = np.random.rand(len(df)) < 0.15
df["is_test"] = msk

cols_drop = ["customerID","Churn", "is_test"]
id_target = df[cols_drop]

df = pd.get_dummies(df.drop(cols_drop, axis=1))
df = pd.concat((id_target, df),axis=1)

df["Churn"] = df["Churn"].apply(lambda x: 0 if x=="No" else 1)

## Immediately Remove High Correlates
corr_matrix = df.corr().abs()

# Select upper triangle of correlation matrix
upper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(np.bool))

# Find index of feature columns with correlation greater than 0.90
to_drop = [column for column in upper.columns if any(upper[column] > 0.999999)]

to_drop

df = df.drop(to_drop, axis=1)

train = df[df["is_test"]==False]

X_train = df[df["is_test"]==False].drop(cols_drop,axis=1)
y_train = df[df["is_test"]==False]["Churn"]

from sklearn.model_selection import KFold
from sklearn.metrics import roc_auc_score
from lightgbm import LGBMClassifier
import gc
import random

params = {
         'boosting_type': 'gbdt',
          'max_depth': 6,
          'objective': 'binary',
          'n_estimators': 5000, 
          'nthread': 12,
          'num_leaves': 32,
          'learning_rate': 0.02,
          'subsample': 0.85,
          'colsample_bytree': 0.95,
          'reg_alpha': 0.5,
          'reg_lambda': 0.8,
          'min_split_gain': 0.05,
          'min_child_weight': 40,
          'min_child_samples': 5,
          'metric': 'auc'}



def pars(n_estimator_bottom, n_estimator_top):
    return random.sample(range(n_estimator_bottom, n_estimator_bottom+ n_estimator_top), 8) + random.sample(range(n_estimator_bottom, int(n_estimator_bottom + n_estimator_top/np.log(n_estimator_top))), 6) + random.sample(range(n_estimator_bottom, int(n_estimator_bottom+(n_estimator_top/np.log(n_estimator_top))/np.log(n_estimator_top))), 4)

    #return random.sample(range(n_estimator_bottom, n_estimator_bottom+ n_estimator_top), 4) + random.sample(range(n_estimator_bottom, int(n_estimator_bottom + n_estimator_top/np.log(n_estimator_top))), 3) + random.sample(range(n_estimator_bottom, int(n_estimator_bottom+(n_estimator_top/np.log(n_estimator_top))/np.log(n_estimator_top))), 2)



mdl = LGBMClassifier(
    **params
)

gridParams = {
    'learning_rate': [0.02, 0.05],
    'n_estimators': pars(10,15000),
    'num_leaves': pars(20,100),
    'boosting_type' : ['gbdt'],
    'random_state' : [405], # Updated from 'seed'
    'colsample_bytree' : [round(random.uniform(0.7,1),2) for _ in range(18)] ,
    'subsample' : [round(random.uniform(.2,.95),2) for _ in range(18)] ,
    'reg_alpha' : [round(random.uniform(0.1,4),2) for _ in range(18)],
    'reg_lambda' : [round(random.uniform(0.1,4),2) for _ in range(18)],
    }

X_test = df.loc[df.is_test, :].drop(cols_drop, axis=1)
y_test = df.loc[df.is_test, :]["Churn"]

from sklearn.model_selection import KFold, StratifiedKFold
from sklearn.metrics import roc_auc_score

n_splits = 5
cvv = StratifiedKFold(n_splits=n_splits, random_state=42)

oof_preds = np.zeros(X_train.shape[0])

sub = df.loc[df.is_test, :][cols_drop]
sub["Churn_Predict"] = 0
feature_importances = pd.DataFrame()
avg_iter = 0
for i, (fit_idx, val_idx) in enumerate(cvv.split(X_train, y_train)):
    
    X_fit = X_train.iloc[fit_idx]
    y_fit = y_train.iloc[fit_idx]
    X_val = X_train.iloc[val_idx]
    y_val = y_train.iloc[val_idx]
    
    model = LGBMClassifier(
**params
    )

    model.fit(
        X_fit,
        y_fit,
        eval_set=[(X_fit, y_fit), (X_val, y_val)],
        eval_names=('fit', 'val'),
        eval_metric='auc',
        early_stopping_rounds=150,
        verbose=False
    )
    
    print("itteration: ", model.best_iteration_)
    avg_iter += model.best_iteration_
    oof_preds[val_idx] = model.predict_proba(X_val, num_iteration=model.best_iteration_)[:, 1]
    sub['Churn_Predict'] += model.predict_proba(X_test, num_iteration=model.best_iteration_)[:,1]
    
    fi = pd.DataFrame()
    fi["feature"] = X_train.columns
    fi["importance"] = model.feature_importances_
    fi["fold"] = (i+1)
    
    feature_importances = pd.concat([feature_importances, fi], axis=0)
    
    print("Fold {} AUC: {:.8f}".format(i+1, roc_auc_score(y_val, oof_preds[val_idx])))
    
print('Full AUC score %.8f' % roc_auc_score(y_train, oof_preds)) 
print('Out of sample AUC: ', roc_auc_score(y_test, sub['Churn_Predict']/n_splits))
print("avg iteration: ",(avg_iter/n_splits))

## The splits seems to work well.
## You can take the average iteration above or use the process below.
## You can retrain below on an ensemble of different seeds to produce the final model

## This will work if ou change parameters to get the new iteration score. 

from lightgbm import cv
from lightgbm import Dataset


def get_score(X, y, usecols, params,seeder,depther):  
     dtrain = Dataset(X[usecols], y)
     params["max_depth"] = depther
     ## These are forced for variable importance
     params["feature_fraction"] = 0.9
     params["bagging_fraction"] = 0.9

     eval =  cv(params,
             dtrain,
             nfold=5,
             stratified=True,
             num_boost_round=20000,
             early_stopping_rounds=160, ## After it stopped how long should go on. 
             verbose_eval=20,
             #verbose_eval=-1,
             seed = seeder,
             show_stdv=True)
     return max(eval['auc-mean']), eval["auc-mean"].index(max(eval["auc-mean"]))
  
seed_dict = {}
for seeder, depther in zip([5, 1, 9 ,12, 20], [6,6,6,6,6]):  ## Here is where you add parameter adaptions
    score, seed_best = get_score(X_train,y_train ,list(X_train.columns), params, seeder,depther)
    seed_dict[seeder] = seed_best

## Lets see if a full dataset prediction validation outperforms.

## This process gives you 15% extra data, but you have no validation set. 
# This can of course be improved by saving the best iterations, which is what I would do
dip  = {}
dip[1] = [6, 5]
dip[2] = [6, 1]
dip[3] = [6, 9]
dip[4] = [6, 12]

dip[5] = [6, 20]

n_splits = 5
cvv = StratifiedKFold(n_splits=n_splits, random_state=42)

sub = df.loc[df.is_test, :][cols_drop]
sub["Churn_Predict"] = 0

feature_importances = pd.DataFrame()

ba= 0 
for i, r in enumerate(range(n_splits)):
    ba = ba + 1
    
    params["num_boost_round"] = seed_dict[dip[ba][1]]  
    params["seed"] = dip[ba][1]
    params["max_depth"] = dip[ba][0] # if you want to customise uncomment
    
    model = LGBMClassifier(
**params
    )

    model.fit(
        X_train,
        y_train,
        eval_metric='auc',
        verbose=False
    )
    
    preds = model.predict_proba(X_test, num_iteration=seed_dict[dip[ba][1]] )[:,1]
    sub['Churn_Predict'] += preds
    
    fi = pd.DataFrame()
    fi["feature"] = X_train.columns
    fi["importance"] = model.feature_importances_
    fi["fold"] = (i+1)
    
    feature_importances = pd.concat([feature_importances, fi], axis=0)
    print("Fold {} AUC: {:.8f}".format(i+1, roc_auc_score(y_test, preds)))
  

print('Out of sample AUC: ', roc_auc_score(y_test, sub['Churn_Predict']/n_splits))

feat_imp = feature_importances.groupby("feature").mean().sort_values("importance", ascending=False)

feat_imp.head(10)

# if is_test shows significance it would point to leakage issue 

## Often Times the cross-validation splits can just be completely ignored.

#180

import numpy as np
l = [x for x in seed_dict.values()]
best_iterations = np.mean(l)

best_iterations

### best_iterations  394.6 for inst
### 252 -  litle subsample

import lightgbm as lgb
d_train = lgb.Dataset(X_train, label=y_train)
d_test = lgb.Dataset(X_test, label=y_test)
# 1.1 for no leaks, 0.7 used because best model feature importance 
# this slightly outperforms enemble but it leaks
params["num_boost_round"] = int(best_iterations*1.1) 
params["num_boost_round"] = int(best_iterations*1.1)


## Here you have to force feature fractions for predictor analysis

import shap

params["num_boost_round"] = int(best_iterations*1.1)
splits = 5 
f = -1
for seeds in [15,1,6,7,2]:
    f = f +1
    print(seeds)
    params["feature_fraction_seed"] = seeds
    
    params["random_seed"] = seeds + 1
    
    model = lgb.train(params, d_train, verbose_eval=1000)

    shap_values = shap.TreeExplainer(model).shap_values(df.drop(cols_drop, axis=1))

    shap_fram = pd.DataFrame(shap_values[:,:-1], columns=list(X_train.columns))

    shap_new = shap_fram.sum().sort_values().to_frame()

    shap_new.columns = ["SHAP"]

    shap_new["Largest_Effect"] = shap_new["SHAP"].apply(lambda x: "-" if x<0else "+")

    shap_new["SHAP_abs"] = shap_new["SHAP"].abs()

    direction = df.corr()["Churn"].to_frame()
    direction.columns = ["corr"]
    shap_new = pd.merge(shap_new, direction, left_index=True, right_index=True, how="left")

    shap_new["direction"] = shap_new["corr"].apply(lambda x: "-" if x<0else "+")
    shap_new["abs_corr"] = shap_new["corr"].abs()
    
    if f==0:
        shap_fin = shap_new
    else:
        shap_fin = shap_fin + shap_new

from lightgbm import plotting

plotting.create_tree_digraph(model, tree_index=1)

plotting.plot_tree(model)

model.plot_metric()

shap_fin.loc[:,["SHAP","SHAP_abs","corr","abs_corr"]] = shap_fin.loc[:,["SHAP","SHAP_abs","corr","abs_corr"]]/5

shap_fin["SHAP_corr_corrrected"] = 0

shap_fin.head()

b = "Contract_Month-to-month"
das  = df.corr()[df.corr().index.isin(list(shap_fin.index))].abs()

shap_fin.head()

for b in das.iloc[:,2:].columns:
    cora = df.corr()[df.corr().index.isin(list(shap_fin.index))][b].abs()
    for c, v in zip(cora.index, cora.values):
        shap_fin["SHAP_corr_corrrected"].loc[c] = shap_fin["SHAP_corr_corrrected"].loc[c] + (shap_fin["SHAP_abs"].loc[b]/shap_fin["SHAP_abs"].max())*v

shap_fin["SHAP_corr"] = shap_fin["SHAP_corr_corrrected"].copy()

shap_fin["SHAP_corr"] = shap_fin["SHAP_corr"]/(shap_fin["SHAP_corr"].max()/shap_fin["SHAP_abs"].max())

shap_fin["SHAP_corr_corrrected"] = shap_fin["SHAP_corr_corrrected"]/(shap_fin["SHAP_corr_corrrected"].max()/shap_fin["SHAP_abs"].max())

shap_fin["SHAP_corr_corrrected"] = shap_fin["SHAP_corr_corrrected"] + shap_fin["SHAP_abs"]

shap_fin["SHAP_corr_corrrected"] = shap_fin["SHAP_corr_corrrected"]/(shap_fin["SHAP_corr_corrrected"].max()/shap_fin["SHAP_abs"].max())

shap_fin.sort_values('SHAP_corr', ascending=False).head(10)

dans = pd.read_csv("amalgamated.csv")

dans.sort_values('SHAP_corr', ascending=False).head(10)

### There is a way to prove that it works - using all. 
#shap_fin.sort_values('SHAP_corr', ascending=False).to_csv("amalgamated.csv")

"SHAP values are sensitive to high correlations among different features. When features are correlated, their impact on the model score can be split among them in an infinite number of ways. This means that the SHAP values will be lower than if all but one of the correlated feature(s) had been removed from the model. "

shap_fin.sort_values("SHAP_corr_corrrected",ascending=False).head(20)

## I would argue that SHAP_corr are more causal in manner.
## Other shit is correlated with SHAP_corr, so if it falls away not too much issues.
shap_fin.sort_values("SHAP_corr",ascending=False).head(20)

## New Number I want to TRY
shap_fin = shap_fin.sort_values("SHAP_abs",ascending=False)

shap_new.to_csv("imp1.csv")

preds = model.predict(X_test, num_iteration=model.best_iteration)

print('Out of sample AUC: ', roc_auc_score(y_test, preds))

# To be truthfull the prediction excerse is not that interesting, what is more interesting
# fitting the data to the response, i.e. training the model and looking at the feature
# interactions it proposes. Feature importance is only interesting to the extent it is
# different to correlation plot, in that regard, higher dimensional interactions become
# important

## To get direction use general df
## df.corr()["Churn"]
shap.summary_plot(shap_values, df.drop(cols_drop, axis=1))

shap.TreeExplainer(xgb).shap_interaction_values(df.drop(cols_drop, axis=1))

shap_fram = pd.DataFrame(shap_values[:,:-1], columns=list(X_train.columns))

shap_new = shap_fram.sum().sort_values().to_frame()

shap_new.columns = ["SHAP"]

shap_new["Largest_Effect"] = shap_new["SHAP"].apply(lambda x: "-" if x<0else "+")

shap_new["SHAP_abs"] = shap_new["SHAP"].abs()

df.drop(cols_drop, axis=1).shape

print('Out of sample AUC: ', roc_auc_score(y_test, preds))

## Here at this point, you migth want to investigate individual mistakes from the validation
## set
## You can simply focus on train on validation instead of train on test.
## However I doubt any leakage would occur here as there is not much to be gained.

sub["Churn_Predict"] = sub['Churn_Predict']/n_splits
sub["Churn"] = y_test

sub["Error"] = sub["Churn"] - sub["Churn_Predict"]

## Worse false positive
sub.sort_values("Error").head(2)

## All worse false positives:
## Contract Month to Months seems to be the culprit, it would be good idea to look at its distribution 
## before and after the 

sub
shap.force_plot(shap_values[list(sub.sort_values("Error").head(100).index),:], df.drop(cols_drop, axis=1).iloc[list(sub.sort_values("Error").head(100).index),:], link="logit")

## The distributions are not different, so this is not a concern 

X_train["Contract_Month-to-month"].sum()/len(X_train)

X_test["Contract_Month-to-month"].sum()/len(X_test)

## Different model giving different amounts, run initial model for shap values.
shap.initjs()

shap.force_plot(shap_values[5354,:], df.drop(cols_drop, axis=1).iloc[5354,:], link="logit")

shap.force_plot(shap_values[345,:], df.drop(cols_drop, axis=1).iloc[345,:], link="logit")

sub["Error_Abs"] = sub["Error"].abs()

## Best True Negative
sub[sub["Churn"]==0].sort_values("Error_Abs").head(2)

shap.force_plot(shap_values[5804,:], df.drop(cols_drop, axis=1).iloc[5804,:], link="logit")

shap_ordered = shap_new.sort_values("SHAP_abs",ascending=False)

# Remember this direction says the average effect on the overall outcome.
# The direction is generally not too informative as it is not the same
# As the direction of the variable relationship with response, for that
# you would have to look at the chart. 
shap_ordered[shap_ordered["Direction"]=="+"].head()

## All this says is that low tenure is more informative than high tenure. 
shap_ordered[shap_ordered["Direction"]=="-"].head()

df["CurrentRateTotal"] = df["tenure"] * df["TotalChargesConv"]
df["RateDiff"] = df["CurrentRateTotal"] - df["TotalChargesConv"] 

df[df.index.isin([4484, 5804])]

## Worse false negative
sub.sort_values("Error").tail(2)

shap.force_plot(shap_values[4272,:], df.drop(cols_drop, axis=1).iloc[4272,:], link="logit")

## Best True Positive
sub[sub["Churn"]==1].sort_values("Error_Abs").head(2)

sub_full = pd.merge(sub,df.drop(cols_drop, axis=1),left_index=True, right_index=True, how="left" )

sub_full.shape

sub_full.head()

# Churn - Predict
sub_full[(sub_full["tenure"]<20) &(sub_full["Contract_Month-to-month"]==1)&(sub_full["MonthlyCharges"]>50) ]["Churn_Predict"].mean()

# Churn - Predict
sub_full[(sub_full["tenure"]<20) &(sub_full["Contract_Month-to-month"]==1)&(sub_full["MonthlyCharges"]>50) ]["Churn"].mean()

# No Chrun - Predict
sub_full[(sub_full["tenure"]>20) &(sub_full["Contract_Month-to-month"]==0)&(sub_full["MonthlyCharges"]<50) ]["Churn_Predict"].mean()

# No Churn - Actual
sub_full[(sub_full["tenure"]>20) &(sub_full["Contract_Month-to-month"]==0)&(sub_full["MonthlyCharges"]<50) ]["Churn"].mean()

df.head()

## New Features - Lets See Effect



X_train = df[df["is_test"]==False].drop(cols_drop,axis=1)
y_train = df[df["is_test"]==False]["Churn"]

X_test = df.loc[df.is_test, :].drop(cols_drop, axis=1)
y_test = df.loc[df.is_test, :]["Churn"]

## This process gives you 15% extra data, but you have no validation set. 
# This can of course be improved by saving the best iterations, which is what I would do
dip  = {}
dip[1] = [6, 5]
dip[2] = [6, 1]
dip[3] = [6, 9]
dip[4] = [6, 12]
dip[5] = [6, 20]

n_splits = 5
cvv = StratifiedKFold(n_splits=n_splits, random_state=42)

sub = df.loc[df.is_test, :][cols_drop]
sub["Churn_Predict"] = 0

feature_importances = pd.DataFrame()

ba= 0 
for i, r in enumerate(range(n_splits)):
    ba = ba + 1
    
    params["num_boost_round"] = seed_dict[dip[ba][1]]  
    params["seed"] = dip[ba][1]
    params["max_depth"] = dip[ba][0] # if you want to customise uncomment
    
    model = LGBMClassifier(
**params
    )

    model.fit(
        X_train,
        y_train,
        eval_metric='auc',
        verbose=False
    )
    
    preds = model.predict_proba(X_test, num_iteration=seed_dict[dip[ba][1]] )[:,1]
    sub['Churn_Predict'] += preds
    
    fi = pd.DataFrame()
    fi["feature"] = X_train.columns
    fi["importance"] = model.feature_importances_
    fi["fold"] = (i+1)
    
    feature_importances = pd.concat([feature_importances, fi], axis=0)
    print("Fold {} AUC: {:.8f}".format(i+1, roc_auc_score(y_test, preds)))
  
print('Out of sample AUC: ', roc_auc_score(y_test, sub['Churn_Predict']/n_splits))

## It has led to worse performance 
Out of sample AUC:  0.8648277170687089

#Start here for XGBoost two interactions:
import xgboost

d_train = xgboost.DMatrix(X_train, label=y_train)
d_test = xgboost.DMatrix(X_test, label=y_test)
    
params = {
    "eta": 0.01,
    "objective": "binary:logistic",
    "colsample_bytree":0.2,
    "colsample_bylevel":0.2,
    "subsample": 0.5,
    "base_score": np.mean(y_train),
    "eval_metric": "logloss"
}
xgb = xgboost.train(params, d_train, 15000, evals = [(d_test, "test")], verbose_eval=100, early_stopping_rounds=200)

import _pickle as pickle
#pickle.dump(shap_inter, open("shap_inter.p", "wb"))  # save it into a file 
shap_inter = pickle.load(open("shap_inter.p", "rb")) # open pickle

shap_inter[1,1,:]

import shap

shap.summary_plot(shap_inter, df.drop(cols_drop, axis=1))

## This is a very important interaction that can't be over-emphesised. 

df.head()

shap.dependence_plot(
    ("tenure", "gender_Female"),
    shap_inter, df.drop(cols_drop, axis=1),
    display_features=df.drop(cols_drop, axis=1)
)

shap.dependence_plot(
    ("tenure", "SeniorCitizen_0"),
    shap_inter, df.drop(cols_drop, axis=1),
    display_features=df.drop(cols_drop, axis=1)
)

shap.dependence_plot(
    ("tenure", "TotalChargesConv"),
    shap_inter, df.drop(cols_drop, axis=1),
    display_features=df.drop(cols_drop, axis=1)
)

df.head()

shap.dependence_plot(
    ("TotalChargesConv","Contract_Month-to-month"),
    shap_inter, df.drop(cols_drop, axis=1),
    display_features=df.drop(cols_drop, axis=1)
)

shap.dependence_plot(
    ("MonthlyCharges","Contract_Two year"),
    shap_inter, df.drop(cols_drop, axis=1),
    display_features=df.drop(cols_drop, axis=1)
)

shap.dependence_plot(
    ("TotalChargesConv","MonthlyCharges"),
    shap_inter, df.drop(cols_drop, axis=1),
    display_features=df.drop(cols_drop, axis=1)
)

cols_drop

kt

df[kt].head()

df["Churn"].value_counts()

import seaborn as sns

kt = list(shap_fin.sort_values("SHAP_corr",ascending=False).index[:4])
kt = []
kt.append("Churn")
kt.append("MonthlyCharges")
kt.append("TotalChargesConv")
g = sns.pairplot(df[kt], hue='Churn', palette = 'seismic',size=3,diag_kind = 'kde',diag_kws=dict(shade=True),plot_kws=dict(s=10) )
g.set(xticklabels=[])

## You are more likely to churn with high monthly charges and low total cost.
## You are not as likley to churn with high monthly charges and high monthly costs.
## See it shows the same relationship
## What it does unfortunately not show is the probability of being this and that as it is single points. 
## For that you might have to do a nearest neighbour transofmration or a strategised groupby.
## and that is simply what the tree provides, see it as a groupby function, in out example it is just not 
## entirely overfitted, so some information might be lost, but not necessarily if generalisability is important. 
## With scatter plots you have to do it with actual or expected values. 
## I have neverv seen expected value plots because to create it you would have to create some groupbys or knn's
## so inface your LightGBM is the expected value amounts - nothing can actually get as close as it, hence why it
## is used. And the expected value plot gives you a slightly better picture than colored scatter plots for categorical.
## However, for regression plots, it is going to show a very similar plot. - You can at a a later stage try that out 
## for regression plots. So I would recommend using the shap dependence plots when it comes to classification
## and maybe switching  to scater plots when it comes to regression. There is also an argument to just use the
## the scatter plots as the relationship might become moot as a result of another correlated vairable.