%reload_ext autoreload
%autoreload 2

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

import matplotlib.pylab as pylab
import pandas as pd
from operator import itemgetter

from warnings import filterwarnings

pylab.rcParams['figure.figsize'] = 5, 5
plt.rc("font", size=10)


plt.rc('xtick', color='k', labelsize='medium', direction='in')
plt.rc('xtick.major', size=8, pad=12)
plt.rc('xtick.minor', size=8, pad=12)

plt.rc('ytick', color='k', labelsize='medium', direction='in')
plt.rc('ytick.major', size=8, pad=12)
plt.rc('ytick.minor', size=8, pad=12)

df_h_n=pd.read_csv("data/Human_nuisance.csv", index_col=0)
df_h_n.rename(columns={ "Breeding density(individuals per ha)":"Breeding",
                       "Number of pedestrians per ha per min":"N"},inplace=True)

fig, ax = plt.subplots(figsize=(8,6))
ax.plot(df_h_n.N,df_h_n.Breeding,'ro')
ax.set_ylabel("Breeding density_h_n(individuals per ha)")
ax.set_xlabel("Number of pedestrians per ha per min")
plt.show()

fig, ax = plt.subplots(figsize=(8,6))
ax.plot(df_h_n.N[:-2],df_h_n.Breeding[:-2],'ro')
ax.set_ylabel("Breeding density(individuals per ha)")
ax.set_xlabel("Number of pedestrians per ha per min")
plt.show()

import statsmodels
import statsmodels.api as sm
import statsmodels.formula.api as smf

list_co=[]#list of covariable


df_nuisance=df_h_n.copy()

models = []
logLikelihoods = []

for i in range(1,5):

    df_nuisance["N"+str(i)]=df_nuisance.N**i
    
    list_co.append( "N"+str(i) )

    ## create the model, without the last 2 points
    model = smf.ols("Breeding ~ " + "+".join(list_co) , data = df_nuisance.iloc[:-2, :] )
    results = model.fit()#we do the actual fit
    
    models.append("+".join(list_co))
    logLikelihoods.append( results.llf )
    
logLikelihoods

import scipy.stats as stats 

for i in range(1,len(logLikelihoods)):
    #calculating the pvalue for the LRT between the models
    pval=1-stats.chi2.cdf(2*(logLikelihoods[i]-logLikelihoods[i-1]),1)
        
    print("model {:<10} - model {:<12} : LRT p-value = {:.2e}".format(models[i-1],
                                                                     models[i],
                                                                     pval))

from sklearn.metrics import r2_score
from sklearn.metrics import mean_squared_error
fig, ax = plt.subplots(1, 2,figsize=(14,7))
## Quadratic model

for i, variables in enumerate( ["N1+N2","N1+N2+N3"] ):

    model = smf.ols("Breeding ~ "+variables , data = df_nuisance.iloc[:-2] )
    results = model.fit()#we do the actual fit
    
    predicted = results.predict( df_nuisance )
    
    ## computing goodness of fit metrics 
    R2_2_h_n = r2_score(df_nuisance.Breeding[:-2], predicted[:-2])
    MSE_2_h_n = mean_squared_error(df_nuisance.Breeding[:-2],predicted[:-2])

    R2_22_h_n = r2_score(df_nuisance.Breeding[-2:], predicted[-2:])
    MSE_22_h_n = mean_squared_error(df_nuisance.Breeding[-2:], predicted[-2:])


    ax[i].scatter( x = df_nuisance.N , y = df_nuisance.Breeding , c = ['teal']*(df_h_n.shape[0]-2) + ['orange']*2 )
    ax[i].plot(df_nuisance.N[:-2],predicted[:-2] ,'r-')
    ax[i].plot(df_nuisance.N[-3:],predicted[-3:] ,'r--')

    ax[i].set_title(variables+'\nall but 2 : R2={0:.2f}, MSE={1:.2f}\n last 2 : R2={2:.2f}, MSE={3:.2f}'.format(R2_2_h_n,MSE_2_h_n,R2_22_h_n,MSE_22_h_n),fontsize=13)

from utils import poly_fit_train_test

cubic=[]
quadratic=[]

## we reshape N to ake it compatible with the sklearn functions we use there
Nreshaped = np.array( df_nuisance.N ).reshape(-1,1)

for i in range(3):# here we are fitting our model and checking it on different random subsample of the data 
    seed = np.random.randint( 10**4 )
    fig, ax = plt.subplots( 1 , 2 , figsize=(10,5) )
    
    cubic_metrics = poly_fit_train_test( Nreshaped,
                                        df_nuisance.Breeding, seed = seed, deg = 3 , ax = ax[0])#this contain the fit and some scoring metric
    quad_metrics  = poly_fit_train_test( Nreshaped,
                                        df_nuisance.Breeding, seed = seed, deg = 2 , ax = ax[1])
    
    cubic.append(cubic_metrics)
    quadratic.append(quad_metrics)
   

cubic=[]
quadratic=[]
for i in range(500):# same as before but on way more different split 
    temp3=poly_fit_train_test(Nreshaped,df_nuisance.Breeding,deg=3,ax=None)
    temp2=poly_fit_train_test(Nreshaped,df_nuisance.Breeding,deg=2,ax=None)
    
    if min(temp3)>0 and min(temp2)>0:
        cubic.append(temp3)
        quadratic.append(temp2)

cubic_known , cubic_new =  list( zip( *cubic ) )
quadratic_known , quadratic_new =  list( zip( *quadratic ) )

fig, ax = plt.subplots(1, 2,figsize=(14,7))
ax[0].hist([cubic_known,quadratic_known],label=['cubic','quadratic'])
ax[0].set_title('Known')
ax[0].legend(loc='best')
ax[0].set_xlabel('R2')

ax[1].hist([cubic_new,quadratic_new],label=['cubic','quadratic'])
ax[1].set_title('New')
ax[1].legend(loc='best')
ax[1].set_xlabel('R2')
plt.show()

from sklearn.datasets import load_diabetes
diabetes = load_diabetes()

df_diabetes=pd.DataFrame(diabetes['data'],
                         columns=  diabetes['feature_names'])

df_diabetes.head()

%%time
from sklearn.linear_model import SGDRegressor
from sklearn.metrics import r2_score

# the SGDRegressor from sklearn
# uses gradient descent to find the best fit with the 
# objective function modified with the regularization term

logalphas = []

coef_dict = {'name' : [],
             'val' : [],
             'log-alpha' : []}
r2 = []

for alpha in np.logspace(-1,1,50):

    reg_diabetes = SGDRegressor( penalty='l1' , alpha = alpha , max_iter=10000 )
    reg_diabetes.fit( df_diabetes , diabetes['target'] )
    
    logalphas.append(np.log10(alpha))
    r2.append( r2_score( diabetes['target'] , reg_diabetes.predict(df_diabetes) ) )
    
    coef_dict['name'] += list( df_diabetes.columns )
    coef_dict['val'] += list( reg_diabetes.coef_ )
    coef_dict['log-alpha'] += [np.log10(alpha)]* len(df_diabetes.columns )

coef_df = pd.DataFrame(coef_dict)

import seaborn as sns

fig,ax = plt.subplots(1,2,figsize = (14,7))

ax[0].plot(logalphas , r2)
ax[0].set_xlabel("log10( alpha )")
ax[0].set_ylabel("R2")

sns.lineplot( x = 'log-alpha' , y='val' , hue = 'name' , data= coef_df , ax = ax[1])

fig.suptitle("regression of diabetes data with an L1 regularization.")

df = pd.read_csv("data/potato_data.phenotypic.csv" , index_col=0)
y = df["Flesh Colour"]
y.describe()

dfTT = pd.read_csv("data/potato_data.transcriptomic.top200norm.csv" , index_col = 0)
dfTT.head()

%%time
X = np.array( dfTT )

I = list( range(X.shape[0]))
np.random.shuffle( I ) 

# we will use the first 30 points to evaluate the model
X_valid = X[ I[:30] , :] 
y_valid = y.iloc[I[:30]]

# we will use the rest to train the model
X_train = X[ I[30:] , :] 
y_train = y.iloc[I[30:]]


logalphas = []

r2_train = []
r2_valid = []

for alpha in np.logspace(-3,2,200):

    reg_diabetes = SGDRegressor( penalty='l1' , alpha = alpha , max_iter=2000 )
    reg_diabetes.fit( X_train , y_train )
    
    logalphas.append(np.log10(alpha))
    r2_train.append( r2_score( y_train , reg_diabetes.predict(X_train) ) )
    r2_valid.append( r2_score( y_valid , reg_diabetes.predict(X_valid) ) )
    
## plotting and reporting 
bestI = np.argmax(r2_valid)
bestLogAlpha = logalphas[bestI]
bestR2_valid = r2_valid[bestI]

fig,ax = plt.subplots(figsize = (10,5))
fig.suptitle("best alpha : {:.2f} - validation R2 : {:.2f}".format(10**bestLogAlpha , bestR2_valid))
ax.plot( logalphas, r2_train , label='train set' )
ax.plot( logalphas, r2_valid , label='validation set' )
ax.scatter( [bestLogAlpha] , [bestR2_valid]  , c='red')
ax.set_xlabel("log10( alpha )")
ax.set_ylabel("R2")
ax.set_ylim(-0.1,1.1)
ax.legend()

## I make a helper function with tests a bunch for a given validation set of alpha values 
def testAlphas( X_train , y_train , X_valid , y_valid , alphas):
    r2_valid = []

    for alpha in alphas:
        reg_diabetes = SGDRegressor( penalty='l1' , alpha = alpha , max_iter=2000 )
        reg_diabetes.fit( X_train , y_train )
        r2_valid.append( r2_score( y_valid , reg_diabetes.predict(X_valid) ) )
        
    return r2_valid

k = 4 ## 4-fold

## ugly initialization , we have to handle the fact that the number of sample is not perfectly divisible by 4
I = np.array( ( [0,1,2,3] * ((X.shape[0]//k)+1) )[:X.shape[0]] )

## shuffle randomly
np.random.shuffle( I ) 
print( I )

alphas = np.logspace(-1,1,200)

cumulated_r2 = np.zeros( alphas.shape )


fig,ax = plt.subplots(figsize = (10,5))


for fold in range(k):

    X_valid = X[ I==fold , :] 
    y_valid = y[I==fold]

    # we will use the rest to train the model
    X_train = X[ I!=fold , :] 
    y_train = y[I!=fold]
    
    r2 = testAlphas( X_train , y_train , X_valid , y_valid , alphas)
    cumulated_r2 += r2
    
    
    ax.plot( np.log10(alphas) , r2 , label = 'fold '+str(fold) , linestyle='dotted' )

## cross-validated r2 is the average across each fold:
cross_valid_r2 = cumulated_r2 / k

bestI = np.argmax(cross_valid_r2)
bestLogAlpha = np.log10( alphas[bestI] )
bestR2_valid = cross_valid_r2[bestI]

ax.plot( np.log10(alphas), cross_valid_r2 , label='cross-validated r2' )
ax.scatter( [bestLogAlpha] , [bestR2_valid]  , c='red')
ax.set_xlabel("log10( alpha )")
ax.set_ylabel("R2")
ax.set_ylim(-0.1,1.1)
ax.legend()

fig.suptitle("best alpha : {:.2f} - cross-validated R2 : {:.2f}".format(10**bestLogAlpha , bestR2_valid))

df_diabetes.head()

diabetes['target'][:5]

from sklearn.model_selection import train_test_split

#always split your dataset and do the fitting on the training
# the train_test_split will split the data into train and test set
# let's keep 25% of the data aside for the final test
X_diabetes_train, X_diabetes_test, y_diabetes_train, y_diabetes_test = train_test_split(df_diabetes, 
                                                                                        diabetes['target'],
                                                                                        test_size = 0.25,
                                                                                        random_state=1425)
print('train set size', len(y_diabetes_train))
print('test  set size', len(y_diabetes_test ))

%%time
from sklearn.linear_model import SGDRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline
from sklearn.model_selection import GridSearchCV

# while not strictly necessary for this dataset, 
# regression model normally need the X values to be properly scaled.
# this is something that was done under the hood by statsmodels, but with sklearn
# we have to handle it explicitely.
#
# NB : in that particular dataset the data is already scaled, 
#      but we wanted to show this important step anyway

# in our model we will also consider polynomials for our features,
# the PolynomialFeatures object generates all polynomials of a given degree for us


pipeline_reg_diabetes=Pipeline([('scalar',StandardScaler()),
                                ('poly',PolynomialFeatures(include_bias=False)),
                                ('model',SGDRegressor())])



# define the hyperparameters you want to test
# with the range over which you want it to be tested.
# Here the hyperparamters 
#   * the degree of the polynomials
#   * the form of regularization : l1 or l2
#   * the strength of regularization : alpha
grid_values = {'poly__degree': np.arange(1,4,1),
               'model__penalty':['l1','l2'],
               'model__alpha':np.logspace(-1,1,100)}


#Feed the pipeline and set of values to the GridSearchCV with the 
# score over which the decision should be taken (ir R^2).
# and the cross-validation scheme, here the number of fold in a stratified k-fold strategy
grid_reg_diabetes_r2 = GridSearchCV(pipeline_reg_diabetes, 
                                    param_grid = grid_values, 
                                    scoring='r2', 
                                    cv = 5,
                                    n_jobs=-1)

# Where the actual fit happens
#  the gridSearchCV object will go through each hyperparameter value combination
#  and fit + evaluate each fold, and averages the score across each fold
#
#  it then finds the combination that gave the best score and
#  use it to re-train a model with the whole train data
grid_reg_diabetes_r2.fit(X_diabetes_train, y_diabetes_train)

#get the best cross-validated score 
print('Grid best score ('+grid_reg_diabetes_r2.scoring+'): ',
      grid_reg_diabetes_r2.best_score_)
# print the best parameters
print('Grid best parameter :')
for k,v in grid_reg_diabetes_r2.best_params_.items():
    print(' {:>20} : {}'.format(k,v))

# get the equivalent score on the test set
y_diabetes_decision_fn_scores_r2=grid_reg_diabetes_r2.score(X_diabetes_test,
                                                            y_diabetes_test)

print('Grid best parameter (max. '+grid_reg_diabetes_r2.scoring+') model on test: ',
      y_diabetes_decision_fn_scores_r2)

# here is a recipe to get the name of the features, with their power 

best_reg = grid_reg_diabetes_r2.best_estimator_['model']
poly     = grid_reg_diabetes_r2.best_estimator_['poly']

coef_names = ['_'.join([df_diabetes.columns[j]+'^'+str(poly.powers_[i][j])for j in range(len(df_diabetes.columns)) if poly.powers_[i][j]>0]) for i in range(len(poly.powers_)) ] 

sorted_features=sorted( [(coef_names[i],abs(best_reg.coef_[i])) for i in range(len(poly.powers_))] ,
                       key=itemgetter(1),reverse=True)

print('Important features')

for feature, weight in sorted_features:
    print('\t{:>10}\t{:.3f}'.format(feature,weight) )

## train prediction
y_diabetes_train_pred = grid_reg_diabetes_r2.predict(X_diabetes_train)
## test prediction
y_diabetes_test_pred  = grid_reg_diabetes_r2.predict(X_diabetes_test)


plt.scatter( y_diabetes_train_pred , y_diabetes_train , c="blue" , label='train set'  )
plt.scatter( y_diabetes_test_pred , y_diabetes_test , c="red" , label='test set'  )
m,M = min(y_diabetes_train_pred) , max(y_diabetes_train_pred)
plt.plot( [m,M] , [m,M] , 'k--' )
plt.xlabel( 'predicted values' ) 
plt.ylabel( 'real values' ) 
plt.legend()

grid_reg_diabetes_r2.best_estimator_['scalar']

X1 = np.concatenate( [ np.random.randn(300) , np.random.randn(300)+4 ])
y = np.array( [0]*300 + [1]*300 )

sns.histplot( x=X1,hue = y )

from sklearn.linear_model import LogisticRegression

# do not forget to scale the data
X1_norm = StandardScaler().fit_transform(X1.reshape( X1.shape[0] , 1 ))

fig,ax = plt.subplots( figsize = (10,5) )

ax.scatter( X1_norm , y , c = y )

for alpha in [0.01,0.1,1,10]:
    
    # this implementation does not take alpha but rather C = 1/alpha
    C = 1/alpha
    lr = LogisticRegression( penalty = 'l2' , C = C )
    lr.fit(X1_norm , y)
    
    proba = lr.predict_proba(np.linspace(-2,2,100).reshape(-1, 1))
    ax.plot( np.linspace(-2,2,100) , proba[:,1] , label = 'alpha = {}'.format(alpha) )
ax.legend()

from sklearn.metrics import accuracy_score
y_predicted = lr.predict(X1_norm)

print( f"Accuracy with a threshold of 0.5 : {accuracy_score(y,y_predicted):.2f}"  )

pd.crosstab( y , y_predicted )

threshold = 0.2
y_predicted = lr.predict_proba(X1_norm)[:,1] > threshold
print( f"Accuracy with a threshold of {threshold} : {accuracy_score(y,y_predicted):.2f}"  )
pd.crosstab( y , y_predicted )

from sklearn.datasets import load_breast_cancer
#loading the dataset which is comprised in scikit learn already
data = load_breast_cancer()

## we reduce the features because otherwise this problem is a bit too easy ;-)
m = list( map( lambda x : x.startswith("mean ") , data["feature_names"] ) )


X_cancer=data['data'][:,m]
y_cancer= 1-data['target']

#making it into a dataframe
breast_cancer_df=pd.DataFrame(X_cancer,
    columns=data["feature_names"][m])

breast_cancer_df["target"]=y_cancer

breast_cancer_df.head()

breast_cancer_df.target.value_counts()

#split your data

# stratify is here to make sure that you split keeping the repartition of labels unaffected
X_train_cancer, X_test_cancer, y_train_cancer, y_test_cancer = train_test_split(X_cancer,
                                                                                y_cancer,
                                                                                random_state=0,
                                                                                stratify=y_cancer)


print("fraction of class malignant in train",sum(y_train_cancer)/len(y_train_cancer))
print("fraction of class malignant in test ",sum(y_test_cancer)/len(y_test_cancer) )
print("fraction of class malignant in full ",sum(y_cancer)/len(y_cancer))

%%time
## create the pipeline of data handling :
## scaling, then logistic regression

## in order for the LogisticRegression object takes solver="liblinear" argument
## in order to be able to test the l1 and l2 norm
## other solver exists, and are more or less performant / have different capabilities
pipeline_lr_cancer=Pipeline([('scaler',StandardScaler()),
                             ('model',LogisticRegression(solver = 'liblinear'))])

from sklearn.model_selection import GridSearchCV
from sklearn.metrics import classification_report

grid_values = {'model__C': np.logspace(-5,2,100),
               'model__penalty':['l1','l2']}

grid_lr_cancer_acc = GridSearchCV(pipeline_lr_cancer, 
                                  param_grid = grid_values,
                                  scoring='accuracy',
                                  cv=5,
                                  n_jobs=-1)

grid_lr_cancer_acc.fit(X_train_cancer, y_train_cancer)#train your pipeline

#get the best cross-validated score 
print('Grid best score ('+grid_lr_cancer_acc.scoring+'): ',
      grid_lr_cancer_acc.best_score_)
# print the best parameters
print('Grid best parameter :')
for k,v in grid_lr_cancer_acc.best_params_.items():
    print(' {:>20} : {}'.format(k,v))

1/0.02

# calculate the score of your trained pipeline on the test
cancer_acc = grid_lr_cancer_acc.score(X_test_cancer,y_test_cancer)
print('Grid best parameter (max.'+grid_lr_cancer_acc.scoring+') model on test: ', 
      cancer_acc) 

y_test_cancer.shape

# we can then access the coefficient of the model, to assess the importance of the different parameters:

w_lr_cancer=grid_lr_cancer_acc.best_estimator_['model'].coef_[0]

sorted_features=sorted( zip( breast_cancer_df.columns , np.abs( w_lr_cancer ) ),
                       key=itemgetter(1),
                       reverse=True)

print('Features sorted per importance in discriminative process')
for f,ww in sorted_features:
    print('{:>25}\t{:.3f}'.format(f,ww))

from sklearn.metrics import accuracy_score, confusion_matrix
import seaborn as sns


y_cancer_test_pred = grid_lr_cancer_acc.predict(X_test_cancer)

# get the confusion matrix:
confusion_m_cancer = confusion_matrix(y_test_cancer, y_cancer_test_pred)

## recipe to plot the confusion matrix : 
plt.figure(figsize=(5,4))
sns.heatmap(confusion_m_cancer, annot=True)
plt.title('Accuracy:{0:.3f}'.format(accuracy_score(y_test_cancer,y_cancer_test_pred)))
plt.ylabel('True label')
plt.xlabel('Predicted label')

from scipy.special import expit
from sklearn.metrics import roc_curve
from sklearn.metrics import auc

# 1. the decision_function gives you the "score" for a point to be in a class
y_score_lr_cancer = grid_lr_cancer_acc.decision_function(X_test_cancer)

#   for the logistic regression, that score is the logit, so can convert it back to 
#   probabilities with the expit function (which is the inverse of the logit)
y_proba_lr_cancer = expit(y_score_lr_cancer )

# 2. this calculates the ROC curve : TPR and FPR for each threshold of score
fpr_lr_cancer, tpr_lr_cancer, threshold_cancer = roc_curve(y_test_cancer, y_proba_lr_cancer)

# we find the point corresponding to a 0.5 theshold
keep = np.argmin( np.abs( threshold_cancer - 0.5 ) )

# we compute the area under the ROC curve
roc_auc_lr_cancer = auc( fpr_lr_cancer, tpr_lr_cancer )

# 3. plotting 
plt.figure()
plt.xlim([-0.01, 1.01])
plt.ylim([-0.01, 1.01])
plt.plot(fpr_lr_cancer, tpr_lr_cancer, lw=3, label='LogRegr ROC curve\n (area = {:0.2f})'.format(roc_auc_lr_cancer))
plt.plot(fpr_lr_cancer[keep], tpr_lr_cancer[keep],'ro',label='threshold=0.5')
plt.xlabel('False Positive Rate', fontsize=16)
plt.ylabel('True Positive Rate', fontsize=16)
plt.title('ROC curve (logistic classifier)', fontsize=16)
plt.legend(loc='lower right', fontsize=13)
plt.plot([0, 1], [0, 1], color='navy', lw=3, linestyle='--')
plt.show()

X1 = np.concatenate( [ np.random.randn(300) , np.random.randn(30)+2 ])
y = np.array( [0]*300 + [1]*30 )

# do not forget to scale the data
X1_norm = StandardScaler().fit_transform(X1.reshape( X1.shape[0] , 1 ))

fig,ax = plt.subplots(1,2, figsize = (14,5) )

sns.histplot( x=X1,hue = y , ax =ax [0])


ax[1].scatter( X1_norm , y , c = y )

for alpha in [0.01,0.1,1,10]:
    
    # this implementation does not take alpham but rather C = 1/alpha
    C = 1/alpha
    lr = LogisticRegression( penalty = 'l2' , C = C )
    lr.fit(X1_norm , y)
    
    proba = lr.predict_proba(np.linspace(-2,3,100).reshape(-1, 1))
    ax[1].plot( np.linspace(-2,3,100) , proba[:,1] , label = 'alpha = {}'.format(alpha) )
ax[1].legend()

y_predicted = lr.predict(X1_norm)
print( f"Accuracy with a threshold of 0.5 : {accuracy_score(y,y_predicted):.2f}"  )
pd.crosstab( y , y_predicted )

from sklearn.metrics import recall_score

## RECALL = TP / (TP + FN)

recall_list = []
acc_list = []

imbalance_list = np.linspace(0,1,50)

alpha = 1
for imbalance in imbalance_list:

    n0 = 300
    n1 = int( n0 * (1 - imbalance) )
    if n1 == 0:
        n1 = 1
    
    X1 = np.concatenate( [ np.random.randn(n0) , np.random.randn(n1)+2 ])
    y = np.array( [0]*n0 + [1]*n1 )

    X1_norm = StandardScaler().fit_transform(X1.reshape( X1.shape[0] , 1 ))
    
    C = 1/alpha
    lr = LogisticRegression( penalty = 'l2' , C = C )
    lr.fit(X1_norm , y)
    
    y_predicted = lr.predict(X1_norm)
    
    recall_list.append( recall_score( y , y_predicted ) )
    acc_list.append( accuracy_score(y,y_predicted) )

        
fig,ax=plt.subplots(figsize = (10,4))
ax.plot( imbalance_list , acc_list , label='accuracy' )
ax.plot( imbalance_list , recall_list , label='recall' )
ax.set_xlabel("imbalance")
ax.legend()

from sklearn.metrics import balanced_accuracy_score


def check_imbalance_effect( imbalance_list , class_weight = None):
    
    recall_list = []
    balanced_acc_list = []
    acc_list = []
    for imbalance in imbalance_list:

        n0 = 300
        n1 = int( n0 * (1 - imbalance) )
        if n1 == 0:
            n1 = 1

        X1 = np.concatenate( [ np.random.randn(n0) , np.random.randn(n1)+2 ])
        y = np.array( [0]*n0 + [1]*n1 )

        X1_norm = StandardScaler().fit_transform(X1.reshape( X1.shape[0] , 1 ))

        # LR
        lr = LogisticRegression( penalty = 'l2' , C = 1 , class_weight=class_weight)
        lr.fit(X1_norm , y)

        y_predicted = lr.predict(X1_norm)

        recall_list.append( recall_score( y , y_predicted )  )
        acc_list.append( accuracy_score(y,y_predicted) )
        balanced_acc_list.append( balanced_accuracy_score(y,y_predicted) )

    return recall_list , acc_list , balanced_acc_list

imbalance_list = np.linspace(0,1,50)

### first, we see what happens without class_weight=None

recall_list , acc_list , balanced_acc_list = check_imbalance_effect( imbalance_list , 
                                                               class_weight = None)
    
        
fig,ax=plt.subplots(1,2,figsize = (12,4))
ax[0].plot( imbalance_list , acc_list , label='accuracy - class_weight=None' )
ax[0].plot( imbalance_list , recall_list , label='recall - class_weight=None' )
ax[0].plot( imbalance_list , balanced_acc_list , label='balanced_accuracy - class_weight=None' )
ax[0].set_xlabel("imbalance")
ax[0].set_ylim(0,1)
ax[0].legend()
## now, with class weight 

recall_list , acc_list , balanced_acc_list = check_imbalance_effect( imbalance_list , 
                                                               class_weight = 'balanced')
            
ax[1].plot( imbalance_list , acc_list , label='accuracy - class_weight=balanced' )
ax[1].plot( imbalance_list , recall_list , label='recall - class_weight=balanced' )
ax[1].plot( imbalance_list , balanced_acc_list , label='balanced_accuracy - class_weight=balanced' )
ax[1].set_xlabel("imbalance")
ax[1].set_ylim(0,1)
ax[1].legend()

%%time
## create the pipeline of data handling :
## scaling, then logistic regression

## in order for the LogisticRegression object takes solver="liblinear" argument
## in order to be able to test the l1 and l2 norm
## other solver exists, and are more or less performant / have different capabilities
pipeline_lr_cancer=Pipeline([('scaler',StandardScaler()),
                             ('model',LogisticRegression(solver = 'liblinear' ))])

from sklearn.model_selection import GridSearchCV
from sklearn.metrics import classification_report

grid_values = {'model__C': np.logspace(-5,2,100),
               'model__penalty':['l1','l2'],
               'model__class_weight':[None,'balanced']}

grid_lr_cancer_acc = GridSearchCV(pipeline_lr_cancer, 
                                  param_grid = grid_values,
                                  scoring='balanced_accuracy',
                                  cv=5,
                                  n_jobs=-1)

grid_lr_cancer_acc.fit(X_train_cancer, y_train_cancer)#train your pipeline

#get the best cross-validated score 
print('Grid best score ('+grid_lr_cancer_acc.scoring+'): ',
      grid_lr_cancer_acc.best_score_)
# print the best parameters
print('Grid best parameter :')
for k,v in grid_lr_cancer_acc.best_params_.items():
    print(' {:>20} : {}'.format(k,v))

# calculate the score of your trained pipeline on the test
cancer_acc = grid_lr_cancer_acc.score(X_test_cancer,y_test_cancer)
print('Grid best parameter (max.'+grid_lr_cancer_acc.scoring+') model on test: ', 
      cancer_acc) 

# we can then access the coefficient of the model, to assess the importance of the different parameters:

w_lr_cancer=grid_lr_cancer_acc.best_estimator_['model'].coef_[0]

sorted_features=sorted( zip( breast_cancer_df.columns , np.abs( w_lr_cancer ) ),
                       key=itemgetter(1),
                       reverse=True)

print('Features sorted per importance in discriminative process')
for f,ww in sorted_features:
    print('{:>25}\t{:.3f}'.format(f,ww))

from sklearn.metrics import accuracy_score, confusion_matrix
import seaborn as sns


y_cancer_test_pred = grid_lr_cancer_acc.predict(X_test_cancer)

# get the confusion matrix:
confusion_m_cancer = confusion_matrix(y_test_cancer, y_cancer_test_pred)

## recipe to plot the confusion matrix : 
plt.figure(figsize=(5,4))
sns.heatmap(confusion_m_cancer, annot=True)
plt.title('Accuracy:{0:.3f}'.format(accuracy_score(y_test_cancer,y_cancer_test_pred)))
plt.ylabel('True label')
plt.xlabel('Predicted label')

df_heart=pd.read_csv('data/framingham.csv')

df_heart.dropna(axis=0,inplace=True) # removing rows with NA values.

print(df_heart.shape)
df_heart.head()

##separation in X and y
X_heart = df_heart.drop( columns = "TenYearCHD" )
y_heart = df_heart[ "TenYearCHD" ]

# %load -r -20 solutions/solution_03_01.py

# %load -r 22-30 solutions/solution_03_01.py

# %load -r 32-50 solutions/solution_03_01.py

# %load -r 51-71 solutions/solution_03_01.py

# %load -r 72-96 solutions/solution_03_01.py

# %load -r 97- solutions/solution_03_01.py

np.random.seed(10815657)


fig,ax = plt.subplots(1,2,figsize=(15,7))

# case 1 
norm1=0.2*np.random.randn(100)-2
norm2=0.8*np.random.randn(100)+2.5
ax[0].plot(norm1,[1]*100,'ro',markersize=10)
ax[0].plot(norm2,[1]*100,'bo')

min2 = min( norm2 )
max1 = max( norm1 )

ax[0].axvline(min2, color='k', linestyle='--', label='defined by the most extreme blue point')
ax[0].axvline( (min2+max1)/2,color='k',linestyle='-.',label='middle of extreme of two classes')
ax[0].legend(loc='best')


# case 2
cauch=0.8*np.random.standard_cauchy(10)-2
norm=1*np.random.randn(100)+2.5

ax[1].plot(cauch,[1]*10,'ro',markersize=10)
ax[1].plot(norm,[1]*100,'bo')

min2 = min( norm )
max1 = max( cauch )

ax[1].axvline(min2, color='k', linestyle='--', label='defined by the most extreme blue point')
ax[1].axvline( (min2+max1)/2,color='k',linestyle='-.',label='middle of extreme of two classes')
ax[1].legend(loc='best')

from sklearn.datasets import make_blobs
X2, y2 = make_blobs(n_samples=(250,250), centers=[[-1,-1],[1,1]], random_state=6)
plt.scatter(X2[:,0],X2[:,1],c=y2)
plt.xlim(min(X2[:,1])-0.5,max(X2[:,1])+0.5)
plt.ylim(min(X2[:,1])-0.5,max(X2[:,1])+0.5)
plt.show()

from utils import contour_SVM

contour_SVM(X2,y2,c=100,ker='linear')

contour_SVM(X2,y2,c=0.01,ker='linear')

contour_SVM(X2,y2,c=1,ker='rbf',gam=10)

contour_SVM(X2,y2,c=1,ker='rbf',gam=1)

contour_SVM(X2,y2,c=1,ker='poly',deg=5)

grid_values = {"model__kernel": ['linear', 'rbf', 'poly'],
                 "model__C":np.logspace(-2, 2, 10),
                 "model__degree":[2,3,4,5],
                 "model__gamma": np.logspace(-2,1,10)}

%%time
from sklearn.svm import SVC

## set up the pipeline as usual
pipe = Pipeline([('scalar',StandardScaler()),
                 ("classifier", SVC())])

# the grid of parameter values is not a dictionnary, but now a list 
#  of dictionnaries : smaller grid to explore independently

# for each of these grid, we can re-define locally the classifier
grid_param = [  {"classifier": [SVC(class_weight='balanced', probability=True, kernel='linear')],
                 "classifier__C":np.logspace(-2, 2, 10)},
                {"classifier": [SVC(class_weight='balanced', probability=True, kernel='rbf')],
                 "classifier__C":np.logspace(-2, 2, 10),
                 "classifier__gamma": np.logspace(-2,1,10)},
                {"classifier": [SVC(class_weight='balanced', probability=True, kernel='poly')],
                 "classifier__C":np.logspace(-2, 2, 10),
                 "classifier__degree":[2,3,4,5]}]
## NB : we use probability = True in order to make ROC auc computation possible later on

grid_svm = GridSearchCV(pipe, 
                        grid_param, 
                        cv=5, 
                        verbose=0,
                        n_jobs=-1,
                        scoring='accuracy') 

grid_svm.fit(X_train_cancer, y_train_cancer)

#get the best cross-validated score 
print('Grid best score ('+grid_svm.scoring+'): ',
      grid_svm.best_score_)
# print the best parameters
print('Grid best kernel    :' , grid_svm.best_params_["classifier"].kernel )
print('Grid best parameter :')
for k,v in grid_svm.best_params_.items():
    print(' {:>20} : {}'.format(k,v))

%%time
from sklearn.svm import SVC

scoring_metric = "accuracy"

# linear kernel
pipe1 = Pipeline([('scalar',StandardScaler()),
                 ("classifier", SVC(class_weight='balanced', probability=True, kernel='linear'))])
grid_param1 = { "classifier__C":np.logspace(-2, 2, 10)}

grid_svm_linear = GridSearchCV(pipe1, grid_param1, 
                        cv=5, verbose=0, n_jobs=-1, scoring=scoring_metric) 

# rbf kernel
pipe2 = Pipeline([('scalar',StandardScaler()),
                 ("classifier", SVC(class_weight='balanced', probability=True, kernel='rbf'))])
grid_param2 = { "classifier__C":np.logspace(-2, 2, 10) , 
               "classifier__gamma": np.logspace(-2,1,10)}

grid_svm_rbf = GridSearchCV(pipe2, grid_param2, 
                        cv=5, verbose=0, n_jobs=-1, scoring=scoring_metric) 

# poly kernel
pipe3 = Pipeline([('scalar',StandardScaler()),
                 ("classifier", SVC(class_weight='balanced', probability=True, kernel='poly'))])
grid_param3 = { "classifier__C":np.logspace(-2, 2, 10) , 
               "classifier__degree": [2,3,4,5]}

grid_svm_poly = GridSearchCV(pipe3, grid_param3, 
                        cv=5, verbose=0, n_jobs=-1, scoring=scoring_metric) 



grid_svm_linear.fit(X_train_cancer, y_train_cancer)
grid_svm_rbf.fit(X_train_cancer, y_train_cancer)
grid_svm_poly.fit(X_train_cancer, y_train_cancer)


print('linear Grid best score ('+grid_svm_linear.scoring+'): ',grid_svm_linear.best_score_)
print('   rbf Grid best score ('+grid_svm_rbf.scoring+'): '   ,grid_svm_rbf.best_score_)
print('  poly Grid best score ('+grid_svm_poly.scoring+'): '  ,grid_svm_poly.best_score_)

y_cancer_test_score=grid_svm.score(X_test_cancer,y_test_cancer)

print('Grid best parameter (max.'+grid_svm.scoring+') model on test: ', y_cancer_test_score)

y_cancer_pred_test=grid_svm.predict(X_test_cancer)

confusion_m_cancer_SVMC = confusion_matrix(y_test_cancer, y_cancer_pred_test)

plt.figure(figsize=(5,4))
sns.heatmap(confusion_m_cancer_SVMC, annot=True)

plt.ylabel('True label')
plt.xlabel('Predicted label')
plt.title("test {} : {:.3f}".format(grid_svm.scoring , y_cancer_test_score))

y_cancer_score_SVMC = grid_svm.decision_function(X_test_cancer)
fpr_SVMC_cancer, tpr_SVMC_cancer, thre_SVMC_cancer = roc_curve(y_test_cancer, y_cancer_score_SVMC)
roc_auc_SVMC_cancer = auc(fpr_SVMC_cancer, tpr_SVMC_cancer)

proba=expit(thre_SVMC_cancer)
for i in range(len(proba)):
    if abs(proba[i]-0.5)<0.1:
        keep=i
        break
        

plt.figure()
plt.xlim([-0.01, 1.00])
plt.ylim([-0.01, 1.01])
plt.plot(fpr_SVMC_cancer, tpr_SVMC_cancer, lw=3, label='SVC ROC curve\n (area = {:0.2f})'.format(roc_auc_SVMC_cancer))
plt.plot(fpr_SVMC_cancer[keep], tpr_SVMC_cancer[keep],'ro',label='threshold=0.5')
plt.xlabel('False Positive Rate', fontsize=16)
plt.ylabel('True Positive Rate', fontsize=16)
plt.title('ROC curve (SVM classifier)', fontsize=16)
plt.legend(loc='lower right', fontsize=13)
plt.plot([0, 1], [0, 1], color='navy', lw=3, linestyle='--')
#plt.axes().set_aspect('equal')
plt.show()

%%time
from sklearn.inspection import permutation_importance

r = permutation_importance(grid_svm.best_estimator_, X_test_cancer, y_test_cancer,
                            n_repeats=1000,
                            random_state=132987)

for i in r.importances_mean.argsort()[::-1]:
    if r.importances_mean[i] - 2 * r.importances_std[i] > 0:
        print(f"{breast_cancer_df.columns[i]:<8}"
               f"{r.importances_mean[i]:.3f}"
               f" +/- {r.importances_std[i]:.3f}")

%%time
from sklearn.svm import SVR
## SVR for regression

pipeline_SVR=Pipeline([('scalar',StandardScaler()),
                       ('model',SVR(kernel='rbf'))])

grid_values = {"model__C":np.logspace(-2, 2, 20),
               "model__gamma": np.logspace(-2,2,20)}

## don't forget to change the metric to one adapted for regression:
grid_SVR_diabetes = GridSearchCV(pipeline_SVR, 
                                 param_grid = grid_values, 
                                 scoring='r2', n_jobs=-1)

grid_SVR_diabetes.fit(X_diabetes_train, y_diabetes_train)

print('Grid best score ('+grid_SVR_diabetes.scoring+'): ', grid_SVR_diabetes.best_score_)
print('Grid best parameter (max.'+grid_SVR_diabetes.scoring+'): ', grid_SVR_diabetes.best_params_)

y_diabetes_test_score=grid_SVR_diabetes.score(X_diabetes_test,y_diabetes_test)
print('Grid best parameter (max.'+grid_SVR_diabetes.scoring+') model on test: ', y_diabetes_test_score)

## train prediction
y_diabetes_train_pred = grid_SVR_diabetes.predict(X_diabetes_train)
## test prediction
y_diabetes_test_pred  = grid_SVR_diabetes.predict(X_diabetes_test)


plt.scatter( y_diabetes_train_pred , y_diabetes_train , c="blue" , label='train set'  )
plt.scatter( y_diabetes_test_pred , y_diabetes_test , c="red" , label='test set'  )
m,M = min(y_diabetes_train_pred) , max(y_diabetes_train_pred)
plt.plot( [m,M] , [m,M] , 'k--' )
plt.xlabel( 'predicted values' ) 
plt.ylabel( 'real values' ) 
plt.legend()