#!/usr/bin/env python
# coding: utf-8

# Importation des packets.

# In[1]:


import numpy as np
from scipy.stats import norm
from sklearn import datasets
get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib
import matplotlib.pyplot as plt 
from sklearn import linear_model
import os

cwd = os.getcwd()

dirname = "../srcimages/"
if not os.path.exists(dirname):
    os.mkdir(dirname)
imageformat = '.pdf'

def my_saving_display(fig, dirname, filename, imageformat):
    """"saving faster"""
    dirname + filename + imageformat
    image_name = dirname + filename + imageformat
    fig.savefig(image_name)


# Importation de la base de données ''diabetes''.

# In[30]:


datasets.load_diabetes()

diabetes = datasets.load_diabetes()
diabetes_X = diabetes.data
diabetes_y = diabetes.target


# Les variables sont standardisées de la façon suivante:

# In[33]:


print(diabetes_X.sum(axis = 0))
print(np.dot(np.transpose(diabetes_X),diabetes_X))


# Regression avec sk learn

# In[3]:


regr = linear_model.LinearRegression()
regr.fit(diabetes_X , diabetes_y)

print(regr.coef_)
print(regr.intercept_)


# Premier graphe

# In[4]:


ypred_tmp = np.dot(diabetes_X , regr.coef_ )
ypred = ypred_tmp +  regr.intercept_ 
#ypred = regr.predict(diabetes_X)


plt.figure() 
plt.plot(ypred_tmp , ypred)
plt.scatter(ypred_tmp , diabetes_y , s =3 )
plt.title("y et y_pred vs xbeta")
plt.show()


# Méthode Forward : à chaque boucle sur k, on selectionne la variable avec la plus grande p-valeur

# In[5]:


diabetes_X_aug = np.column_stack( (np.ones( (diabetes_X.shape[0], 1 )), diabetes_X )) 

p = diabetes_X_aug.shape[1]
n = diabetes_X_aug.shape[0]

test = np.zeros((p,p))
pval_mem = np.zeros(p)
pval = np.zeros((p,p))
resids = diabetes_y

var_sel = []
var_remain = list(range(p))
in_test = []


regr = linear_model.LinearRegression(fit_intercept = False)


for k in range(p): 
    resids_mem = np.zeros((p,n))

    for i in var_remain:
        xtmp = diabetes_X_aug [:,[i]]
        regr.fit(xtmp , resids)
        
        #calcul de (x'x)
        xx = np.sum( diabetes_X_aug [:,[i]] ** 2 )    
        resids_mem [i,:] = regr.predict(xtmp) - resids
        sigma2_tmp = np.sum (resids_mem [i,:] ** 2) / xx
        test[k,i] = np.sqrt(n) * np.abs(regr.coef_) / (np.sqrt( sigma2_tmp )) 
        pval[k,i] =   2 * (1 - norm.cdf(test[k,i]))

    ####separe en deux vecteurs la listes des variables séléctionnées et les autres
    best_var = np.argmax(test[k,:])
    var_sel.append(best_var)
    resids = resids_mem[best_var,:]
    pval_mem [k] = pval[k,best_var]
    var_remain = np.setdiff1d(var_remain,var_sel)



# Pour les 3 premières étapes de la méthode (k=0,1,2), on trace les valeurs de la statistique de test calculée, enregistrée dans la matrice test.

# In[6]:


def foo(s1):
    return "step %s" % s1

fig = plt.figure()
for k in range(3): 
    lab_tmp = foo(k)
    plt.plot(np.arange(p),test[k,:], '-o',label = lab_tmp)
    plt.axis(xmin = -.5,xmax = 9.5,ymin = -1)
plt.legend(loc=1)
plt.title("values of the t-stat at each steps")
plt.xlabel("features")
plt.show()

filename = 'forw_first_steps'
my_saving_display(fig, dirname, filename, imageformat)


# On peut maintenant regarder, pour chaque étape, la pvaleurs associée à la variable séléctionnée. On voit que cette pvaleur grandie; et que si l'on se fixe un niveau de 0.1 pour le test, on rejette l'hypothèse nulle (H_0 : beta_k = 0) uniquement pour les 6 premières étapes.

# In[7]:


fig2 = plt.figure() 
for k in range(3): 
    plt.plot(np.arange(p),pval_mem, 'o')
    plt.plot([-0.5,9.5],[.1,.1],color = "blue" )
    plt.axis(xmin = -.5,xmax = 9.5,ymin = -.03)
plt.title("plot of the pvalues")
plt.xlabel("steps")
plt.show()

filename = 'pvalues'
my_saving_display(fig2, dirname, filename, imageformat)


# En pratique cela signifie que notre méthode sélectionne les 7 variables suivantes (0 correspondant à l'intercept):

# In[8]:


print("variables sélectionnées (dans l'ordre): " , 
      np.array(var_sel)[pval_mem<.1])


# In[9]:


diabetes_X_aug = np.column_stack( (np.ones( (xtmp.shape[0], 1 )), diabetes_X )) 


# In[ ]: