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

# ***
# ** Algorithmes d'optimisation -- L3 MINT et doubles licences 2017/2018 -- 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 5: Projection, débruitage et algorithme d'Uzawa
# 
# ## Partie I: Construction de l'algorithme
# 
# Dans ce TP, on cherche à appliquer l'algorithme d'Uzawa au calcul de la projection d'un point $x_0\in \Rsp^d$ sur un convexe $K$ de la forme  
# 
# $$ K = \{ x \in \Rsp^d\mid Ax \leq b\} $$
# 
# où $A$ est une matrice possédant $k$ lignes notées $a_1,\hdots,a_k \in \Rsp^d$ et $b\in\Rsp^k$. L'inégalité $A x \leq b$ doit être comprise composante à composante, i.e. $x \in K$ ssi $\forall i, c_i(x) := \sca{a_i}{x} - b_i \leq 0$. On regarde donc le problème de minimisation 
# 
# $$ (P) := \min_{x\in K} \frac{1}{2} \nr{x - x_0}^2 $$
# 
# On rappelle que le lagrangien de ce problème est donné par
# 
# $$ \ell: (x,\lambda)\in\Rsp^d\times \Rsp^k_+ \mapsto f(x) + \sum_{1\leq i\leq k} \lambda_i c_i(x) $$
# 
# Dans la partie suivante, on verra des applications à des modèles de régression. 
# 
# **Q1)** Montrer que le lagrangien associé au problème (P) peut s'écrire $\ell(x,\lambda) = \frac{1}{2}\nr{x - x_0}^2  + \sca{\lambda}{A x - b}$. En déduire que le problème dual est  donné par
# 
# $$(D) := \max_{\lambda \in\Rsp_+^k} - h(\lambda) \hbox{ où } h(\lambda) =  \frac{1}{2} \nr{A^T \lambda - x_0}^2 - \frac{1}{2}\nr{x_0}^2 + \sca{\lambda}{b} $$
# 
# 
# 
# **Q2)** Montrer que $\nabla h(\lambda) = A(A^T \lambda - x_0) + b$. Étant donné $\lambda \in \Rsp^k$, donner l'expression de l'unique solution $x_\lambda$ du problème de minimisation $\min_{x \in \Rsp^d} \ell(x,\lambda).$ Montrer que $\nabla h(\lambda) = b - A x_\lambda$.
# 
# 
# 
# L'algorithme du gradient projeté pour le problème dual (D) (aussi appelé algorithme d'Uzawa) est donné par 
# 
# $$ \begin{cases}
# \lambda^{(0)} = 0 \in \Rsp^k \\
# g^{(k)} = \nabla h(\lambda^{(k)})\\
# \lambda^{(k+1)} = p_{\Rsp_+^k}(\lambda^{(k)} - \tau g^{(k)})\\
# x^{(k+1)} = x_0 - A^T \lambda^{(k+1)} \quad (\in \arg\min_{x\in\Rsp^d} \ell(x,\lambda))
# \end{cases}
# $$
# 
# où l'on rappelle que $p_{\Rsp_+^k}(v) = (\max(v_1,0),\hdots,\max(v_k,0))$. L'algorithme est arrêté lorsque $\nr{x^{(k)} - x^{(k+1)}}\leq \eps$.
# 
# 
# **Q3)** Montrer que si $S(\lambda) := p_{\Rsp_+^k}(\lambda - \tau \nabla G(\lambda))$ est $L$-contractante pour $L<1$, alors la suite $(\lambda^{(k)})$ converge vers l'unique maximiseur $\lambda^*$ de (D).
# 
# 
# On rappelle que si $\lambda^*$ est un maximiseur de (D), alors $x^{*}$ (minimiseur de (P)) est l'unique point dans $\arg\min_{x\in\Rsp^d} \ell(x,\lambda^*)$.
# 
# **Q4)** Écrire l'algorithme sous la forme d'une fonction projection_convexe(A,b,x0,tau,err=1e-6) calculant les itérées de l'algorithme du gradient projeté et arrêtant la boucle dès que $\nr{x^{(k)}- x^{(k+1)}} \leq$ err. Cette fonction retournera $\lambda^{(k)}, x^{(k)}$. Tester cette fonction sur les deux convexes suivants :
# 
# $$ K_1 = \{ x \in\Rsp^2 \mid \nr{x}_\infty \leq 1 \} $$
# 
# $$ K_2 = \{ x \in \Rsp^2\mid \nr{x}_1 \leq 1 \} $$
# 
# On commencera par décrire ces convexes sous la forme $K = \{x \mid Ax \leq b \}$, où $A$ et $b$ sont à déterminer, et on prendra $\tau = 0.1$. On vérifiera la validité du calcul de deux manières:
# 
# - visuellement, en affichant le segment reliant x0 à son projeté, pour un assez grand nombre (100) de points x0 choisis aléatoirement dans $[-4,4]^2$.
# 
# - en vérifiant que la solution $x,\lambda$ = projection_convexe(A,b,x0,tau) satisfait (à erreur numérique près) les quatre conditions du théorème de Karush-Kuhn-Tucker (à savoir: $x \in K$ (soit $Ax \leq b$), $\lambda \geq 0$, $\forall i, (A x - b)_i \lambda_i = 0$ et $\nabla_x \ell(x,\lambda) = 0$).
# 

# In[1]:


import numpy as np
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')
# la commande suivante agrandit les figures
plt.rcParams['figure.figsize'] = [9.,6.]

# <completer>


# In[2]:


# <completer>


# ## Partie II: Régression isotone
# 
# Nous allons considérer deux problèmes de débruitage consistant simplement à projeter sur un convexe. Les données sont par exemple des séries temporelles $y = (y_1,\hdots,y_n)\in \Rsp^n$, mesurées avec un bruit. On sait que les données réelles appartiennent à un certain ensemble convexe $K$ de $\Rsp^n$. Deux exemples
# 
# - régression isotone: $K = \{ x\in \Rsp^n \mid \forall 1\leq i < n,~x_{i+1}\geq x_i \}$
# - régression convexe: $K = \{ x\in \Rsp^n \mid \forall 1 < i < n,~x_{i} \leq \frac{1}{2} (x_{i-1} + x_{i+1}) \}$.
# 
# À cause du bruit, le vecteur $y$ mesuré n'appartient pas à $K$, et on va donc le débruité en le projetant sur $K$, soit:
# 
# $$ (P) := \min_{x\in K} \frac{1}{2} \nr{x - y}^2 $$
# 
# **Q1)** Implémenter la régression isotone en utilisant projection_convexe() avec tau=0.1 et avec le vecteur $y$ donné ci-dessous. Que peut-on dire expérimentalement de la solution $x^*$ aux points $\{i,i+1\}$ où $y_{i}\geq y_{i+1}$ ? Démontrer ce résultat en utilisant les conditions du théorème KKT.

# In[4]:


t = np.linspace(0,1,n)
y = t**2 + .3*np.random.rand(n)

# <completer>


# **Q2)** Implémenter la régression convexe en utilisant projection_convexe() avec tau=0.1 et avec le vecteur $y$ donné ci-dessous. 

# In[5]:


# <completer>