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

# ***
# **Algorithmes d'optimisation -- L3 MINT et doubles licences -- Université Paris-Sud**
# ***
# 
# $\newcommand{\Rsp}{\mathbb{R}}
# \newcommand{\nr}[1]{\|#1\|}
# \newcommand{\abs}[1]{|#1|}
# \newcommand{\eps}{\varepsilon}
# \newcommand{\sca}[2]{\langle#1|#2\rangle}
# \newcommand{\D}{\mathrm{D}}
# \newcommand{\hdots}{\dots}
# \newcommand{\cond}{\mathrm{cond}}$
# 
# # TP 2: Régression logistique
# 
# ## Partie I: Descente de gradient avec rebroussement d'Armijo ##
# 
# Étant donné un point $x^{(k)}$ et une direction de descente $d^{(k)}$ le pas d'Armijo est défini de la manière suivante: 
# 
# $$ t^{(k)}_{\mathrm{arm}} = \max \left\{ t  \mid t = \beta^\ell, \ell\in \mathbb{N}, \hbox{ tq } f(x^{(k)} + t d^{(k)}) \leq f(x) + \alpha t\sca{\nabla f(x)}{d^{(k)}} \right\}, $$
# 
# où on prendra $\alpha = 0.3, \beta = 0.5$. Algorithmiquement, on procède de la manière suivante pour calculer $t^{(k)}_{\mathrm{arm}}$ en fonction de $x=x^{(k)}$ et $v = d^{(k)}$:
# 
# $$
# \left|
# \begin{aligned}
# &\mathrm{pas\_armijo}(f,x,v)\\
# &\qquad t \gets 1\\
# &\qquad m \gets \sca{\nabla f(x)}{v}\\
# &\qquad \textbf{while } f(x + t v) > f(x) + \alpha t m \\
# &\qquad\qquad t \gets \beta t\\
# &\qquad \mathrm{return}~t
# \end{aligned}\right.
# $$
# 
# L'algorithme de descente de gradient avec rebroussement d'Armijo s'écrit alors de la manière suivante:
# 
# $$
# \begin{cases}
# d^{(k)} = - \nabla f(x^{(k)}) \\
# t^{(k)} = \mathrm{pas\_armijo}(f,x^{(k)}, d^{(k)})\\
# x^{(k+1)} = x^{(k)} + t^{(k)} d^{(k)}
# \end{cases}$$
# 
# 
# **QI.1**)  Écrire les fonctions suivantes:
# - `pas_armijo(f,gradf,x)` prend en entrée $f:\Rsp^N\to\Rsp$, `gradf`$=\nabla f:\Rsp^N\to\Rsp^N$, $x,d\in\Rsp^N$. On choisira `alpha`$=0.3$, `beta`$=0.5$ et on retournera le pas $t$. 
# - `gradient_armijo(f,gradf,x0,err=1e-6,maxiter=1000)` prenant en argument la fonction `f`, `gradf`, `x0`=$x^{(0)}$ et un critère d'arrêt `err`$\in\Rsp$, et qui effectue l'algorithme de descente de gradient avec pas d'Armijo en arrêtant la boucle dès que $\nr{d^{(k)}} \leq$ `err` où dès que $k\geq$ `maxiter`. Cette fonction retournera le point $x$ trouvé, un vecteur contenant la suite  $f(x^{(k)})$ et un second vecteur contenant $\nr{\nabla f(x^{(k)})}$. Indications:
#     - **Pour la fonction gradient_armijo, ne pas hésiter à copier-coller gradient_optimal et à faire les modifications nécessaires**
#     - **Ne stocker que l'itération courante et pas la suite des itérées**

# In[13]:


# on importe les modules numpy et pyplot
import numpy as np
import matplotlib.pyplot as plt
# les deux commandes suivante paramètrent l'affichage des figures
get_ipython().run_line_magic('matplotlib', 'inline')
plt.rcParams['figure.figsize'] = [9.,6.]

# les paramètres alpha et beta comportent une valeur par défaut, c'est-à-dire que 
# pas_armijo(f,x,d,m) est équivalent à pas_armijo(f,x,d,m,0.3,0.5)
def pas_armijo(f,gradf,x,v):
    # <completer>
    
# la descente de gradient prend en argument:
# f = la fonction à évaluer
# gradf = le gradient de f
# x0 = le point de départ
def gradient_armijo(f,gradf,x0,err=1e-6,maxiter=1000):
    # ne pas hésiter à copier-coller et adapter gradient_optimal...
    # <completer>


# **QI.2**) Tester la fonction `gradient_armijo` en dimension $N=1$, $f(x) = x^{2} + e^{x}$, $x_0 = 1$ et en prenant pour critère d'arrêt `err`$= 10^{-8}$. Vous devriez trouver $x^* = -0.3517...$. 

# In[14]:


# test de la fonction pas_armijo
# <completer>


# ## Partie II: Régression logistique avec rebroussement d'Armijo en dimension $N=2$
# 
# L'explication de la méthode est à lire dans la feuille de TD.
# On commence par construire les données $(x_a)_{a \in A} \subseteq \Rsp^2$ où :
# - $A = \{1,\hdots,n\} = A_0\sqcup A_1$, $A_0 = \{0, \hdots, m-1,\}$, $A_1 = \{m, 2m-1\}$, $n = 2m$
# - Les points $(x_a)_{a\in A_0}$ sont tirés selon une gaussienne centrée en $(0.5,0.5)$
# - Les points $(x_a)_{a\in A_1}$ sont tirés selon une gaussienne centrée en $(-0.5,-0.5)$.
# 
# On cherche alors à construire, par optimisation de la fonction $F$, un vecteur $w\in\Rsp^2$ tel
# que $\sigma(\sca{w}{x_a}) \simeq 0$ si $a\in A_0$ et  $\sigma(\sca{w}{x_a})\simeq 1$ si $a\in A_1$.

# In[18]:


m = 30
n = 2*m
X = np.vstack((np.hstack((-0.5+.2*np.random.randn(m,1), -0.5 + .2*np.random.randn(m,1))),
               np.hstack((0.5+.2*np.random.randn(m,1), 0.5 + .2*np.random.randn(m,1)))))
A0 = range(0,m)
A1 = range(m,2*m)
plt.plot(X[A0,0],X[A0,1],'.r',label='A0')
plt.plot(X[A1,0],X[A1,1],'.g',label='A1')
plt.legend()
plt.axis('equal');


# **QII.1**: Écrire deux fonctions `F(w,X,A0,A1,gamma)` et `gradF(w,X,A0,A1,\gamma)` calculant la valeur et le gradient de la fonctionnelle 
# 
# $$ F(w) = - \left(\sum_{a\in A_0} \log(1-\sigma(\sca{w}{x_a})) + \sum_{a\in A_1} \log(\sigma(\sca{w}{x_a})\right) + \frac{\gamma}{2} \nr{w}^2. $$
# 
# Tester que `gradF` correspond bien au gradient de `F` en utilisant la fonction `verifier_gradient` fournie ci-dessous.

# In[21]:


# on définit une fonction sigma(t) correspondant à la sigmoïde
def sigma(t):
    # <completer>

def F(w,X,A0,A1,gamma):
    # <completer>

def gradF(w,X,A0,A1,gamma):
    # <completer>

""" la fonction verifier_gradient calcule le gradient de f en un point x0 de deux manières:
    1) via la fonction gradf donnée en argument 
    2) via des différences finies: df/de_i (x_0) ~= (f(x_0+eps e_i) - f(x_0))/eps
    et compare ces deux résultats: si la fonction gradf est correctement programmée, la différence
    entre les deux méthode de calcul devrait être petite """
def verifier_gradient(f,gradf,x0):
    N = len(x0)
    gg = np.zeros(N)
    for i in range(N):
        eps = 1e-5
        e = np.zeros(N)
        e[i] = eps
        gg[i] = (f(x0+e) - f(x0-e))/(2*eps)
    print('erreur numérique dans le calcul du gradient: %g (doit être petit)' % np.linalg.norm(gradf(x0)-gg))
    
# on vérifie l'implémentation de gradF en utilisant la fonction verifier_gradient
gamma = .1
f = lambda w: F(w,X,A0,A1,gamma)
gradf = lambda w: gradF(w,X,A0,A1,gamma)
verifier_gradient(f, gradf, np.random.rand(2))


# **QI.2**: Calculer $w=\arg\min_{\Rsp^2} F$ en utilisant la fonction gradient_armijo. Tracer sur un graphe l'évolution de la norme du gradient au cours des itérations de la descente de gradient, soit $(k,\nr{\nabla F(w^{(k)})})$, où l'axe $y$ sera en échelle logarithmique.

# In[22]:


# <completer>


# **QI.3**: Soit $w$ le minimiseur approché de $F$ calculé dans la question précédente. Tracer sur un graphe:
# - les points $(x_a)_{a\in A_0}$ en rouge
# - les points $(x_a)_{a\in A_1}$ en bleu
# - quelques lignes de niveau de la fonction $u:x\mapsto \sigma(\sca{w}{x})$ en utilisant la fonction fournie `afficher_lignes_niveau(w)`.
# 
# Faire varier le paramètre de régularisation $\gamma \in \{0.01,0.1,1,10\}$ et observer:
# - comment varie le nombre d'itérations.
# - comment varie $\nr{w}$ (où $\nr{w}$ est le minimiseur).
# - quel est l'effet sur la fonction $u:x\mapsto \sigma(\sca{w}{x_i})$.
# 
# <!--
# %On peut expliquer le premier point en se rappelant que
# %Si on note $F_\gamma$ la fonction avec le paramètre $\gamma$, montrer que pour $\gamma' \geq \gamma$ on a 
# %$$ \forall w\in\Rsp^d,\forall v\in \Rsp^d, \sca{D^2 F_\gamma(w)v}{v} \leq \sca{D^2 F_{\gamma'}(w)v}{v} $$
# %En déduire que si $m \leq \D^2 F_\gamma \leq M$ sur $S\subseteq \Rsp^d$, alors 
# %$$ m+(\gamma'-\gamma) \leq \D^2 F_\gamma \leq M + \gamma'-\gamma $$
# %, en notant $\lambda_1(A) < \hdots \leq \lambda_d(A)$ les valeurs propres de $A$,
# %$$ \frac{\lambda_d(A')}{\lambda_1(A')} \leq \frac{\lambda_d(A)}{\lambda_1(A)} $$
# $$
# -->

# In[23]:


# on définit une fonction permettant d'afficher les lignes de niveaux de u_w: R^2 -> R
# definie par u_w(x) = \sigma(<x|w>)
# NB: u peut être définie en utilisant lambda, par exemple:
#    u = lambda (x,y): ... où x,y sont les coordonnées du point
def afficher_lignes_niveau(w):
    Xcontour,Ycontour = np.meshgrid(np.linspace(-1., 1., 100),
                                    np.linspace(-1., 1., 100))
    u = lambda x,y: sigma(x*w[0] + y*w[1])
    Zcontour = u(Xcontour,Ycontour)
    plt.axes([-1.5, -1.5, 1.5, 1.5])
    plt.axis('equal')
    p = plt.contour(Xcontour, Ycontour, Zcontour, cmap='RdBu')
    plt.clabel(p, inline=1, fontsize=10)
    
# <completer>



# ## Partie III: Régression logistique avec méthode de Newton en dimension $N=2$
# 
# **QIII.1**: Écrire une fonction `hessF` calculant la hessienne de $F$, en utilisant la formule obtenue en TD. Vérifier le bon fonctionnement de `hessF` via la fonction `verifier_hessienne` fournie ci-dessous.

# In[7]:


""" cette fonction vérifie numériquement que hessf (=fonction évaluant la hessienne de f) est bien
    la dérivée de grad (= fonction évaluant le gradient de f), en un point x0 """
def verifier_hessienne(gradf,hessf,x0):
    N = len(x0)
    H = np.zeros((N,N))
    for i in range(N):
        eps = 1e-5
        e = np.zeros(N)
        e[i] = eps
        H[i,:] = (gradf(x0+e) - gradf(x0-e))/(2*eps)
    print('erreur numerique dans le calcul de la hessienne: %g (doit etre petit)' % np.sum((H-hessf(x0)))**2)
    
def hessF(w,X,A0,A1,K):
# <completer>

# on vérifie que la hessienne est correcte en utilisant la fonction verifier_hessienne
# <completer>


# **QIII.2**: Écrire une fonction `newton_armijo(f,gradf,hessf,x0,err=1e-6,maxiter=1000)` implémentant la méthode de Newton avec rebroussement d'Armijo, c'est-à-dire l'algorithme suivant:
# 
# $$
# \begin{cases}
# g^{(k)} = \nabla F(x^{(k)}) \\
# d^{(k)} = -[D^2 F(x^{(k)})]^{-1} \nabla F(x^{(k)}) \\
# t^{(k)} = \mathrm{pas\_armijo}(f,x^{(k)},d^{(k)})\\
# x^{(k+1)} = x^{(k)} + t^{(k)} d^{(k)}
# \end{cases}
# $$
# 
# Il ne faut pas hésiter à copier-coller la fonction gradient_armijo et à changer ce qui est nécessaire. Pour la résolution du système linéaire, on utilisera `np.linalg.solve`. 
# 
# Utiliser cette fonction pour calculer le minimum de $F$. Comme dans la question Q1.2, tracer la norme du gradient en fonction de l'itération. Vérifier également que le nombre d'itérations dépend peu du choix du paramètre $\gamma$.

# In[28]:


def newton_armijo(f,gradf,hessf,x0,err=1e-6,maxiter=500):
    # la boucle est tres similaire à celle de gradient_armijo: la seule chose qui
    # change est le choix de la direction de descente. 
    # NB: on utilise np.linalg.solve pour résoudre le système linéaire apparaissant dans la méthode
    # <completer>


# comparer le nombre d'iterations dans la methode de Newton pour gamma dans [0.01,0.1,1,5]
# <completer>


# ## Partie IV: Application à la classification de caractères manuscrits
# 
# L'objectif de cette partie est d'appliquer la méthode décrite dans les parties précédentes (où l'on traitait des données en dimension $N=2$) à un jeu de données plus réaliste et de dimension beaucoup plus grande. L'objectif est de distinguer des images de caractères manuscrits. Pour cela, on a demandé à 6000 personnes d'écrire les chiffres de zéro à neuf. Ici, l'on cherchera uniqument à distinguer les chiffres zéro manuscrits des chiffres un. Chaque point de donnée $x_a$ est construit de la manière suivante: 
# 
# - On demande à une personne d'écrire un chiffre (zéro, un, ..., neuf), que l'on stocke sous la forme d'une image en niveaux de gris de taille 28x28 pixels. On représente cette image par une matrice carrée $p_a \in \mathcal{M}_{28,28}(\Rsp)$ (où $(p_{a})_{ij} = 0$ si le pixel $(i,j)$ est noir et $1$ s'il est blanc). 
# - La matrice carrée $p_a \in \mathcal{M}_{28,28}(\Rsp)$ est convertie en un vecteur $x_a\in\Rsp^{N}$ de taille $N = 28\times 28 = 784$. Pour cela, on met bout-à-bout les lignes de la matrice $p$, via la fonction `np.reshape`.
# - Enfin, on stocke séparément un entier correspondant au chiffre manuscrit (on appelle ça le "label", ou étiquette, de l'image).
# 
# Ce faisant, on obtient donc deux tableaux: 
# - `images` possède $60000$ lignes (soit $10$ chiffres zéro, un, $\hdots$ dessinés par $6000$ personnes) et $784=28^2$  colonnes (images de $28\times 28$ pixels stockées sous formes de vecteurs). `images[i]` correspond alors à la $i$ème image, représentée par un vecteur de $\Rsp^{768}$.
# - `labels` possède $60000$ lignes et une seule colonne: `label[i]` contient le numéro du chiffre manuscrit représenté par l'image `images[i]`
# 
# La fonction `read_mnist()` ci-dessous télécharge ces données et retourne les tableaux `images` et `labels`. On fournit une fonction d'affichage: `show_mnist(images[i])` affiche la $i$ème image.

# In[14]:


# import du jeu de données MNIST de Yan Lecun: 60000 chiffres manuscrits avec labels
import os
import struct
import urllib
import gzip
import shutil
try:
    from urllib.request import urlretrieve
except ImportError:
    from urllib import urlretrieve

def gunzip(iname,oname):
    with gzip.open(iname, 'rb') as f_in:
        with open(oname, 'wb') as f_out:
            shutil.copyfileobj(f_in, f_out)
def download_mnist(dest):
    src = 'http://yann.lecun.com/exdb/mnist'
    for f in ['train-images-idx3-ubyte', 'train-labels-idx1-ubyte']:
        urlretrieve('%s/%s.gz' % (src,f), '%s/%s.gz' % (dest,f))
        gunzip('%s/%s.gz' % (dest,f),'%s/%s' % (dest,f))
def read_mnist():
    d = os.path.expanduser('~/.m315-mnist-data')
    if not os.path.isdir(d):
        os.makedirs(d)
        download_mnist(d)
    with open('%s/train-labels-idx1-ubyte' % d, 'rb') as flbl:
        magic, num = struct.unpack(">II", flbl.read(8))
        labels = np.fromfile(flbl, dtype=np.int8)
    with open('%s/train-images-idx3-ubyte' % d, 'rb') as fimg:
        magic, num, rows, cols = struct.unpack(">IIII", fimg.read(16))
        images = np.fromfile(fimg, dtype=np.uint8).reshape(len(labels), rows*cols)/255.
    return images,labels
def show_mnist(image, cmap='gray'):
    plt.figure()
    imgplot = plt.imshow(np.reshape(image,(28,-1)), interpolation='nearest', cmap=cmap)
   
images,labels = read_mnist()


# **QIV.0** Passer quelques minutes à jouer avec les fonctions show_mnist() et print() (afficher les labels, et aussi les vecteurs images[i]) pour comprendre le jeu de données.

# In[15]:


# exemple d'image: faire varier un peu i pour comprendre le jeu de données
i = 16
x = images[i]
p = np.reshape(images[i],(28,28))
show_mnist(p)
print(x.shape)
print(p.shape)
print(labels[i])


# **QIV.1** Construire un tableau `X0` contenant $m=100$ images du chiffre zéro extraite du tableau `images`, et un tableau `X1` contenant également $m$ représentants du chiffre un également extraits de `images` (`X0` et `X1` possèderont donc $m$ lignes et $784$ colonnes). Construire ensuite les données du problème d'optimisation:
# - Un tableau `X` de $n=2m$ lignes et $784$ colonnes obtenu en mettant à bout `X0` et `X1` (on pourra utiliser np.vstack). 
# - Un tableau `A0` (resp. `A1`) contenant la liste des indices correspondant aux chiffres zéro (resp un).

# In[16]:


# <completer>


# **QIV.2** Résoudre le problème d'optimisation $\min_{w\in\Rsp^d} F(w)$ pour $\gamma=0.1$ en utilisant la méthode de Newton.

# In[12]:


# <completer>


# Étant donné $w^* \in \Rsp^{768}$ l'unique solution du problème d'optimisation, on notera $u^*: x\mapsto\sigma(\sca{x}{w^*})$. Le principe de la régression logistique (cf intro du TD) est que  $u^*(x_a) \simeq 1$ si $a \in A_0$ (c'est-à-dire si l'image $x_a$ représente le chiffre zéro) et $u^*(x_a) \simeq 0$ si $a\in A_1$ (si l'image $x_a$ représente le chiffre un). Comme on a utilisé seulement une petite fraction des images pour construire $w$ (on a choisi $n=100$ images du chiffre zéro alors que la base de données en contient $6000$!), on va maintenant chercher à valider la méthode en comparant la catégorie ($0$ où $1$) devinée par la fonction $u^*$ à celle enregistrée dans le tableau `labels` sur l'ensemble des zéros et des uns manuscrits contenus dans la base de données.
# 
# **QIV.3** Vérifier la pertinence de la méthode de la manière suivante :
# - Pour chaque image MNIST $x = images[i]$ représentant un $0$ où un $1$ (pour tout $i$ tel que labels[i] $\in \{0,1\}$), calculer $u^*(x)$, et définir `label_devine`=0 si $u(x)< 1/2$ et `label_devine`=$1$ sinon.
# - Comparer  label_devine au véritables label (i.e. `labels[i]`).
# - Calculer enfin la proportion de cas où `label_devine=labels[i]`: ce nombre mesure la proportion de cas où la régression logistique permet de reconnaître correctement le chiffre représenté par l'image.

# In[108]:


# <completer>


# **QIV.4:** Afficher le vecteur $w$ sous la forme d'une image en utilisant show_mnist(). Interpréter.

# In[104]:


# <completer>


# **QIV.5** Adapter l'approche à la séparation des chiffres 0 de tout les autres chiffres. (on pourra prendre m=100 chiffres $0$ et $m$ chiffres de chaque autre classe, soit 1000 points de donnée au total).

# In[107]:


# <completer>