{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Expériences numériques du chapitre 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Récupération du schéma via l'opérateur proximal "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Opérateur proximal de la norme de la norme $\\ell_1$ ordonnée :** Cet opérateur est la solution du problème d'optimisation \n",
    "$${\\rm prox}_{\\tau \\lambda}(y)=\\arg \\min_{b\\in \\mathbb{R}^p} \\left\\{ \\frac{1}{2}\\|y-b\\|_2^2+\\tau J_\\lambda(b) \\right\\}.$$\n",
    "Une implémentation naïve de cet opérateur est donnée ci-dessous"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# La suite de Cesàro permet de déterminer les premiers termes de l'opérateur. \n",
    "\n",
    "def Cesaro(y,Lambda):\n",
    "    p = len(y)\n",
    "    C = np.cumsum(y-Lambda)/np.linspace(1,p,p)\n",
    "    k = np.argmax(C)+1\n",
    "    if C[k] <= 0:\n",
    "        k = p\n",
    "        C[-1] = 0\n",
    "    return C[k-1] , k\n",
    "\n",
    "# Opérateur proximal lorsque y1 >= ... >= yp >= 0. \n",
    "\n",
    "def prox_L1_ordonnee_decroissant(y,Lambda):\n",
    "    p = len(y)\n",
    "    prox = np.zeros(p)\n",
    "    i = 0\n",
    "    while i<p :\n",
    "        C, k = Cesaro(y,Lambda)\n",
    "        prox[i:i+k] = C\n",
    "        y = y[k:]\n",
    "        Lambda=Lambda[k:]\n",
    "        i += k\n",
    "    return prox\n",
    "    \n",
    "# Opérateur proximal pour y quelconque.\n",
    "\n",
    "def prox_L1_ordonnee(y,Lambda):\n",
    "    p = len(y)\n",
    "    prox = np.zeros(p)\n",
    "    sign = np.sign(y)\n",
    "    perm = np.argsort(-np.abs(y))\n",
    "    y = np.abs(y)[perm]\n",
    "    prox[perm] = prox_L1_ordonnee_decroissant(y,Lambda)\n",
    "    prox *= sign\n",
    "    return(prox)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Implémentation rapide de l'opérateur proximal :** Cet algorithme est basé sur les maximum locaux de la suite de Cesàro et sur le depistage des composantes nulles de l'opérateur proximal. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prox_L1_ordonnee_rapide_depistage(y,Lambda):\n",
    "    Rep=list()\n",
    "    V=list()\n",
    "    p=len(y)\n",
    "    \n",
    "    y_sign = np.sign(y)\n",
    "    y = np.abs(y)\n",
    "    ord = np.argsort(-y)\n",
    "    y = y[ord]\n",
    "    \n",
    "    l = 0 \n",
    "       \n",
    "    positif = np.sum(1-np.cumprod(np.flip(y<Lambda)))  # Borne supérieure sur le nombre de composantes non-nulles\n",
    "    \n",
    "    if positif <=0:\n",
    "        return np.zeros(p)\n",
    "    \n",
    "    Lambda = y[:positif] - Lambda[:positif] \n",
    "    \n",
    "    while positif>0:\n",
    "        S = Lambda[0]\n",
    "        i = 1\n",
    "        while i<positif and   Lambda[i]  >= S:\n",
    "            S += Lambda[i]\n",
    "            i += 1\n",
    "        Rep.append(i)\n",
    "        V.append(S/i)\n",
    "        l += 1 \n",
    "        Lambda=Lambda[i:]\n",
    "        positif += -i\n",
    "        \n",
    "        while l>1 and V[-1]>V[-2]:\n",
    "            V[-2]=(Rep[-1]*V[-1]+Rep[-2]*V[-2])/(Rep[-2]+Rep[-1])   \n",
    "            Rep[-2]+=Rep[-1]\n",
    "            del V[-1]\n",
    "            del Rep[-1]\n",
    "            l += -1\n",
    "    k=0\n",
    "    for i in range(l):\n",
    "        y[k:k+Rep[i]] = max(V[i],0)\n",
    "        k=k+Rep[i]\n",
    "    \n",
    "    y[ord[:k]] =  y[:k].copy()\n",
    "    y[ord[:k]] *= y_sign[ord[:k]]\n",
    "    y[ord[k:]] = 0\n",
    "\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Image du \"six\" tiré de la base de données MNIST**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAALCUlEQVR4nO3dT4ic9R3H8c+n/rmoh6QZQ4ihayWHSqFRhlBIEYtUYi7Rg8UcJAVhPSgoeKjYgx5DqUoPRVhrMC1WEVTMIbSGIIgXcZQ0fxraWNlqzJKdkIPxZKPfHvaxbOLMzjjP88zzZL/vFwwz++ys882Yd56deWbm54gQgNXve00PAGA6iB1IgtiBJIgdSILYgSSunOaNrVu3LmZmZqZ5k0Aq8/PzOnv2rAd9r1TstrdL+r2kKyT9MSL2rHT9mZkZ9Xq9MjcJYAXdbnfo9yb+Nd72FZL+IOkuSTdL2mX75kn/ewDqVeYx+1ZJH0XExxHxpaRXJO2sZiwAVSsT+0ZJny77+lSx7SK2Z233bPf6/X6JmwNQRpnYBz0J8K3X3kbEXER0I6Lb6XRK3ByAMsrEfkrSpmVf3yDpdLlxANSlTOzvS9ps+0bbV0u6T9L+asYCULWJD71FxAXbD0v6m5YOve2NiOOVTQagUqWOs0fEAUkHKpoFQI14uSyQBLEDSRA7kASxA0kQO5AEsQNJEDuQBLEDSRA7kASxA0kQO5AEsQNJEDuQxFQ/Shqrjz3wU4tbgUVLL8aeHUiC2IEkiB1IgtiBJIgdSILYgSSIHUiC4+zJtfk4eVll/myr8Rg9e3YgCWIHkiB2IAliB5IgdiAJYgeSIHYgCY6zr3JtPo4+6lh2m2e/HJWK3fa8pPOSvpJ0ISK6VQwFoHpV7Nl/HhFnK/jvAKgRj9mBJMrGHpLesv2B7dlBV7A9a7tnu9fv90veHIBJlY19W0TcKukuSQ/Zvu3SK0TEXER0I6Lb6XRK3hyASZWKPSJOF+eLkt6QtLWKoQBUb+LYbV9j+7pvLku6U9KxqgYDUK0yz8avl/RGcSz0Skl/iYi/VjIVVo0y7wvnOHy1Jo49Ij6W9JMKZwFQIw69AUkQO5AEsQNJEDuQBLEDSfAW11WgyUNQTX7k8mr8uOc6sWcHkiB2IAliB5IgdiAJYgeSIHYgCWIHkuA4+2Ug63F0VIs9O5AEsQNJEDuQBLEDSRA7kASxA0kQO5AEx9lbgOPomAb27EASxA4kQexAEsQOJEHsQBLEDiRB7EASxA4kMTJ223ttL9o+tmzbWtsHbZ8sztfUOyaAssbZs78oafsl2x6XdCgiNks6VHwNoMVGxh4R70g6d8nmnZL2FZf3Sbq72rEAVG3Sx+zrI2JBkorz64dd0fas7Z7tXr/fn/DmAJRV+xN0ETEXEd2I6HY6nbpvDsAQk8Z+xvYGSSrOF6sbCUAdJo19v6TdxeXdkt6sZhwAdRn5fnbbL0u6XdI626ckPSlpj6RXbT8g6RNJ99Y55OWO96ujDUbGHhG7hnzrjopnAVAjXkEHJEHsQBLEDiRB7EASxA4kQexAEsQOJEHsQBLEDiRB7EASxA4kQexAEsQOJMGSzasAb2PFONizA0kQO5AEsQNJEDuQBLEDSRA7kASxA0kQO5AEsQNJEDuQBLEDSRA7kASxA0kQO5AEsQNJEDuQxMjYbe+1vWj72LJtT9n+zPbh4rSj3jEBlDXOnv1FSdsHbH82IrYUpwPVjgWgaiNjj4h3JJ2bwiwAalTmMfvDto8Uv+avGXYl27O2e7Z7/X6/xM0BKGPS2J+TdJOkLZIWJD097IoRMRcR3YjodjqdCW8OQFkTxR4RZyLiq4j4WtLzkrZWOxaAqk0Uu+0Ny768R9KxYdcF0A4jPzfe9suSbpe0zvYpSU9Kut32FkkhaV7Sg/WN2H62mx6hMW3+s/N5+hcbGXtE7Bqw+YUaZgFQI15BByRB7EASxA4kQexAEsQOJMGSzcm1+dBZWSv92TIelmPPDiRB7EASxA4kQexAEsQOJEHsQBLEDiTBcfYKjDpmW/ex7NV8rLwuo+6z1Xgcnj07kASxA0kQO5AEsQNJEDuQBLEDSRA7kASxA0kQO5AEsQNJEDuQBLEDSRA7kASxA0kQO5AE72fHisq+r5v32rfHyD277U2237Z9wvZx248U29faPmj7ZHG+pv5xAUxqnF/jL0h6LCJ+JOmnkh6yfbOkxyUdiojNkg4VXwNoqZGxR8RCRHxYXD4v6YSkjZJ2StpXXG2fpLtrmhFABb7TE3S2ZyTdIuk9SesjYkFa+gdB0vVDfmbWds92r9/vlxwXwKTGjt32tZJek/RoRHw+7s9FxFxEdCOi2+l0JpkRQAXGit32VVoK/aWIeL3YfMb2huL7GyQt1jMigCqM82y8Jb0g6UREPLPsW/sl7S4u75b0ZvXjrQ4RseKpzWyXOrXV5fz/ZFLjHGffJul+SUdtHy62PSFpj6RXbT8g6RNJ99YyIYBKjIw9It6VNOyf6DuqHQdAXXi5LJAEsQNJEDuQBLEDSRA7kARvcW2Bppd8vlyt1uPhdWHPDiRB7EASxA4kQexAEsQOJEHsQBLEDiTBcfbLAMeTUQX27EASxA4kQexAEsQOJEHsQBLEDiRB7EASxA4kQexAEsQOJEHsQBLEDiRB7EASxA4kQexAEuOsz77J9tu2T9g+bvuRYvtTtj+zfbg47ah/XACTGufDKy5IeiwiPrR9naQPbB8svvdsRPyuvvEAVGWc9dkXJC0Ul8/bPiFpY92DAajWd3rMbntG0i2S3is2PWz7iO29ttcM+ZlZ2z3bvX6/X25aABMbO3bb10p6TdKjEfG5pOck3SRpi5b2/E8P+rmImIuIbkR0O51O+YkBTGSs2G1fpaXQX4qI1yUpIs5ExFcR8bWk5yVtrW9MAGWN82y8Jb0g6UREPLNs+4ZlV7tH0rHqxwNQlXGejd8m6X5JR20fLrY9IWmX7S2SQtK8pAdrmA9ARcZ5Nv5dSYMWCD9Q/TgA6sIr6IAkiB1IgtiBJIgdSILYgSSIHUiC2IEkiB1IgtiBJIgdSILYgSSIHUiC2IEkiB1IwhExvRuz+5L+s2zTOklnpzbAd9PW2do6l8Rsk6pyth9ExMDPf5tq7N+6cbsXEd3GBlhBW2dr61wSs01qWrPxazyQBLEDSTQd+1zDt7+Sts7W1rkkZpvUVGZr9DE7gOlpes8OYEqIHUiikdhtb7f9T9sf2X68iRmGsT1v+2ixDHWv4Vn22l60fWzZtrW2D9o+WZwPXGOvodlasYz3CsuMN3rfNb38+dQfs9u+QtK/JP1C0ilJ70vaFRH/mOogQ9iel9SNiMZfgGH7NklfSPpTRPy42PZbSeciYk/xD+WaiPh1S2Z7StIXTS/jXaxWtGH5MuOS7pb0KzV4360w1y81hfutiT37VkkfRcTHEfGlpFck7WxgjtaLiHcknbtk805J+4rL+7T0l2XqhszWChGxEBEfFpfPS/pmmfFG77sV5pqKJmLfKOnTZV+fUrvWew9Jb9n+wPZs08MMsD4iFqSlvzySrm94nkuNXMZ7mi5ZZrw1990ky5+X1UTsg5aSatPxv20RcaukuyQ9VPy6ivGMtYz3tAxYZrwVJl3+vKwmYj8ladOyr2+QdLqBOQaKiNPF+aKkN9S+pajPfLOCbnG+2PA8/9emZbwHLTOuFtx3TS5/3kTs70vabPtG21dLuk/S/gbm+Bbb1xRPnMj2NZLuVPuWot4vaXdxebekNxuc5SJtWcZ72DLjavi+a3z584iY+knSDi09I/9vSb9pYoYhc/1Q0t+L0/GmZ5P0spZ+rfuvln4jekDS9yUdknSyOF/botn+LOmopCNaCmtDQ7P9TEsPDY9IOlycdjR9360w11TuN14uCyTBK+iAJIgdSILYgSSIHUiC2IEkiB1IgtiBJP4H6hmsJi30pJYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Position des pixels noirs\n",
    "\n",
    "ind=np.array([99, 100, 101, 102, 103, 125, 126, 127, 128, 129, 130, 131, 132, 133, 152, 153, 154, 155, \n",
    "156, 158, 159, 160, 161, 180, 181, 182, 183, 187, 188, 189, 207, 208, 209, 210, 234, 235, 236, 237, \n",
    "261, 262, 263, 264, 289, 290, 291, 292, 316, 317, 318, 319, 344, 345, 346, 347, 372, 373, 374, 400, \n",
    "401, 402, 428, 429, 430, 455, 456, 457, 458, 463, 464, 465, 466, 467, 483, 484, 485, 486, 489, 490, \n",
    "491, 492, 493, 494, 495, 496, 512, 513, 514, 515, 516, 517, 518, 519, 520, 521, 522, 523, 524, 525,\n",
    "540, 541, 542, 543, 544, 545, 546, 547, 548, 549, 550, 551, 552, 553, 568, 569, 570, 571, 572, 573, \n",
    "576, 577, 578, 579, 580, 581, 597, 598, 599, 600, 601, 602, 603, 604, 605, 606, 607, 608, 626, 627, \n",
    "628, 629, 630, 631, 632, 633, 634])\n",
    "\n",
    "six=np.zeros(784)\n",
    "six[ind]=1\n",
    "image_six=six.reshape(28, 28)\n",
    "plt.imshow(image_six, cmap=\"Greys\")\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Récupération du schéma lorsque l'image précédente est bruitée**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVP0lEQVR4nO3dbWyVZZoH8P8FtIhQ3hGRl50yYARNFvAE1lRGN+OOLx/USdQMH0Y2UZkPmswkfljjfhg/ms3OTCZxMwmsZnAz62SSGaMfzC5YTIgxGgoiwhJeFARKoa0F2vIOvfZDH5MO9rn+x/OcnnN27/8vadqeq/d57vOc5+ppz3W/mLtDRP7/G1fvDohIbSjZRRKhZBdJhJJdJBFKdpFETKjlwWbNmuWLFi3KjY8bV/nvnqJVhevXr4dxM8uNjR8/vuK25WB9Gxoayo01NTWFbdl5K3peo+c06jdrC4x93yPsOWXHjtqz+47O27Fjx9Db2zvqHRRKdjN7CMBvAYwH8O/u/mr084sWLcK2bdty45MnT664L1evXg3j7MLp6+sL45MmTcqN3XzzzWFblnDM4OBgGL948WJubM6cOWFbdlFevnw5jDPReWP3HbUFgCtXroTx6Jpgj5v9Ap8wIU6da9euhfHoemxubg7bXrp0KTfW1taWf8zwXgNmNh7AvwF4GMByAOvMbHml9yciY6vI/+yrARx29y/d/QqAPwJ4rDrdEpFqK5Ls8wEcH/H9iey2v2JmG8ysw8w6ent7CxxORIookuyjvQnwrX+E3H2ju5fcvTR79uwChxORIook+wkAC0d8vwDAyWLdEZGxUiTZdwBYamatZtYM4CcA3q1Ot0Sk2iouvbn7NTN7AcB/Y7j09oa774vajB8/HlOnTs2Ns3pyVC6ZOHFi2Jbd96xZs8J4VPtkxz5//nwYZyUmVo+ePn16bozVbFn5ipWgLly4EMajsiQrh7LHzcqlc+fOzY2xUi173Kx0x/oelWOLlAWj57tQnd3d3wPwXpH7EJHa0HBZkUQo2UUSoWQXSYSSXSQRSnaRRCjZRRJR0/ns7l5o2mE09Y/Vi1m9mdVdb7rpptwYq+GzKYusb6x9NFX0zJkzYVs2rZg9Jy0tLRW3Z7XoInPC2f2zY7MprOx6YWMniqh03Qe9soskQskukgglu0gilOwiiVCyiyRCyS6SiJqW3swsLCMVmY7JVnBlUxZZ+Swqb7HVZYtOM2VTaKPjs9IaKzGxFWCXLVsWxhcvXpwb27NnT9iWrYzLztuuXbtyY0VXri1aNoyux/7+/rBt9JxG15pe2UUSoWQXSYSSXSQRSnaRRCjZRRKhZBdJhJJdJBE1n+Ia1S9ZbTKqpUc7mQJ8Kiabshgdu+h0SbbjJ3tsUZzV+O+6664wzpaKbm1tDePRLkDRtGEAmDFjRhifMmVKGH/uuedyY59++mnYlsXZuA12TUTPGZvSzI6dR6/sIolQsoskQskukgglu0gilOwiiVCyiyRCyS6SiP9T89kvXbqUG2Pz2dn85GnTpoXxCKvRs5ore9xsvvzZs2dzY/fcc0/Ytru7O4wvX748jEdbcAPAsWPHcmNHjhwJ27K+s2NHWzqz8QXsemLYcx49p9F1Xs595ymU7GZ2FMAAgOsArrl7qcj9icjYqcYr+9+7e28V7kdExpD+ZxdJRNFkdwBbzGynmW0Y7QfMbIOZdZhZR09PT8HDiUiliiZ7m7uvAvAwgOfN7Ac3/oC7b3T3kruX2AKCIjJ2CiW7u5/MPncDeBvA6mp0SkSqr+JkN7PJZtbyzdcAfgRgb7U6JiLVVeTd+LkA3s5qxBMA/Ke7/1fU4OrVqzh16lRunNU2o7oqqz12dXWF8QULFoTxaA4xm69+/vz5MM7q6Gy++2233ZYbu++++8K2+/btC+NsfMKBAwfCeFTHZ+MTtm/fHsbXrFkTxiMDAwNhnK0hwObiR2MfgPg5Z9fTuXPncmPR/gcVJ7u7fwngbyttLyK1pdKbSCKU7CKJULKLJELJLpIIJbtIImo6xXXChAnh0sLjxsW/e6ISFFuGev78+WGclYGikgYrs7DtgaNSCsCnwD777LO5MVbOZKW3tWvXhvGPP/44jEdLURedlvzRRx+F8eh6YtcaK72x5ZzZY4umerNS7fTp0yvql17ZRRKhZBdJhJJdJBFKdpFEKNlFEqFkF0mEkl0kEQ21ZfOZM2fC9nPnzs2NsW1uWR2dtY/q+NHYgXKOzerJjz76aBhftGhRboyNP3jwwQfD+JtvvhnGi2wvzGrdbEvm/v7+MF5kOWg2hZVNO2bxaLlo9rjZc5pHr+wiiVCyiyRCyS6SCCW7SCKU7CKJULKLJELJLpKImtbZr169Gi4t3NLSEraP6rJsm1tWk505c2YYj+aUs7onq0W/+OKLYfzkyZNhfOLEibkx9rg/+OCDMM7m4hfBxh+M5VbYlW57XC5W44/6xtYvGBwczI1F6y7olV0kEUp2kUQo2UUSoWQXSYSSXSQRSnaRRCjZRRJR0zp7c3NzuL0wW4s7irPaZDQXHuC18qjGz2q20drpALB///4wzraTjsYYsDXIWT2YzTlnj73Sudfl3Dfb6rqvry83xsYPRGMXAP64iszFZ22jvoVr5Yf3CsDM3jCzbjPbO+K2mWa21cwOZZ9nsPsRkfoq58/43wN46IbbXgLQ7u5LAbRn34tIA6PJ7u7bAdz499BjADZnX28G8Hh1uyUi1VbpG3Rz3b0LALLPt+T9oJltMLMOM+vo6emp8HAiUtSYvxvv7hvdveTupTlz5oz14UQkR6XJftrM5gFA9jl/KpuINIRKk/1dAOuzr9cDeKc63RGRsULr7Gb2FoD7Acw2sxMAfgngVQB/MrNnABwD8GS5BywyLzyqV7N1vgcGBsI4q6t+9dVXubFo3XYAWLlyZRhn/97MmBFXNidMyH8aP/zww7Dt5cuXwzirdbM550XmlLNxF2y9/VtuyX0riY59KHpsdl6jPdbPnTsXtq10jQGa7O6+Lif0w4qOKCJ1oeGyIolQsoskQskukgglu0gilOwiiajpFNdr167h7NmzuXE2nfLixYu5MbZcM9tWObpvAGhtbc2NsX5HJSCAb+/LtrIuct8sXmS5ZiB+XopMf2X3zbDpsdGSzEBc7gT4eYmmJbOy3uTJk8N4Hr2yiyRCyS6SCCW7SCKU7CKJULKLJELJLpIIJbtIImpaZ29qagqnc7L6YjQt8Pjx42HbossSR1Me2TTPzs7OML5kyZIwzqapRlMiv/7667BtdE4BPhWUnbfoOWXTkqMxGQCvN0djCFidnMXZGAE2xTW6nth06+j51pbNIqJkF0mFkl0kEUp2kUQo2UUSoWQXSYSSXSQRNa2zu3tYB2TL90b17GgraIAvz8u2yY2W72W17La2tjB+6NChML506dIwHtWy2Tx+ds7ZXPoi4xPYdtItLS1hnG0nFs13Z8tznzhxIoyzx83GAETXMmsbPe6o/q9XdpFEKNlFEqFkF0mEkl0kEUp2kUQo2UUSoWQXSURN6+xmFtZd+/v7w/ZRLZytf87WAWdzq6O6KKtlHz58OIxPnTo1jEfbRQPxGuXsvLA6O6t1s7nX0fHZ2upsnYCmpqYwHvWdjR9gdXi2PgLrWzRfnj3u6HqJnk/6ym5mb5hZt5ntHXHbK2bWaWa7s49H2P2ISH2V82f87wE8NMrtv3H3FdnHe9XtlohUG012d98OoK8GfRGRMVTkDboXzGxP9md+7j84ZrbBzDrMrIONZRaRsVNpsv8OwPcBrADQBeBXeT/o7hvdveTupWixSREZWxUlu7ufdvfr7j4EYBOA1dXtlohUW0XJbmbzRnz7YwB7835WRBoDrbOb2VsA7gcw28xOAPglgPvNbAUAB3AUwM/KOZi7h/VJVneN2rI6OqubsnW+o/t/+umnw7bz5s0L43v27AnjpVIpjEc13SJ73pfT/o477gjj0RgE9nyzWviRI0fCeLTOABs/wK4nNr6ArSvf15f/nje7VqP7js4pTXZ3XzfKza+zdiLSWDRcViQRSnaRRCjZRRKhZBdJhJJdJBE1neJ6/fr1sOTAppmGZQWyxS4rIbHtoqOpg5999lnYNnrMAMJtrAGgo6MjjA8ODubG2DllZcNTp06FcVZiipZcPnnyZNh2zZo1Yfz2228P46tWrcqNtbe3h22nTZsWxouWDaOSJCsDR6XWqF96ZRdJhJJdJBFKdpFEKNlFEqFkF0mEkl0kEUp2kUQYq5NWU6lU8k8++SS/M6R2GS1LzJbuZcvzMlEte/r06WHb1avjtT3Ylsw7duwI43feeWdujJ3TAwcOhHG2fTCrlV+6dCk3xurovb29YXzfvn1hPKpls+tl9+7dYbzoFNeols6m10bbh7e1tWHnzp2jPul6ZRdJhJJdJBFKdpFEKNlFEqFkF0mEkl0kEUp2kUTUdD47w7bBjeaUs9okq3uyrYuj9qzfbGngCxcuhPG1a9eG8a1bt+bGnnjiibAt2276iy++CONFnjM2PuH8+fNh/Kmnngrj77//fm7s3nvvDduy8QXR+AGAryMQtWdt2doNefTKLpIIJbtIIpTsIolQsoskQskukgglu0gilOwiiahpnd3dw3nlbI7wwMBAbozVJlkt++DBg2F8xYoVuTFWo2dbNrM55Wxb5Wht9m3btoVtDx06FMbZeYvWGACAu+++OzfW2dkZtmW2bNkSxqNrItrOGSh2zgE+RiBal57tM8DWKMhDX9nNbKGZfWBm+81sn5n9PLt9ppltNbND2ed45IiI1FU5f8ZfA/Ciuy8D8HcAnjez5QBeAtDu7ksBtGffi0iDosnu7l3uviv7egDAfgDzATwGYHP2Y5sBPD5GfRSRKvhOb9CZ2fcArATwCYC57t4FDP9CAHBLTpsNZtZhZh09PT0FuysilSo72c1sCoA/A/iFu/eX287dN7p7yd1LbANDERk7ZSW7mTVhONH/4O5/yW4+bWbzsvg8AN1j00URqQZaerPh9/lfB7Df3X89IvQugPUAXs0+v1PGfYXltePHj4ftb7311twYK39NmTIljLPlnKN+nzt3Lmy7adOmQsdm2wdHZaQHHnggbMtKZ08++WQYj5bYBoCVK1fmxk6fPh22ZctUt7a2hvHFixfnxl577bWwLSsDs9Icm4ZaZCnpaLp1FCunzt4G4KcAPjez3dltL2M4yf9kZs8AOAYgvipEpK5osrv7hwDyqvg/rG53RGSsaLisSCKU7CKJULKLJELJLpIIJbtIImo6xXVoaCic+jdr1izavlLjxsW/19h9R/VkVgdnNdmjR4+G8SVLloTxZcuW5cY6OjrCtmyZatZ3NkU2GoPQ3NwctmVLTbe3t4fxqObMpomyLb5Z369cuRLGo76x6bNR36OYXtlFEqFkF0mEkl0kEUp2kUQo2UUSoWQXSYSSXSQRNd+yOap3s1p4tJQ0q02yuujMmTPDeFRnZ7XoaO4ywOcvs1p2hG0tzLYmLjL+AABaWlpyY0eOHAnbLly4MIyzWni0xgGrs7PrhS0VzdZXiJa5ZktJR9tsq84uIkp2kVQo2UUSoWQXSYSSXSQRSnaRRCjZRRJR0zq7mYX1S1Y3jdZ+Z3VRhtX4o2OzOnlTU1MYZ3OfWR0/qstOmjQpbMv6ztqz9dWjMQZs3Xc2RoBt0x09NrauOxtfEI0fAPhzFsXZOY/W+o/myeuVXSQRSnaRRCjZRRKhZBdJhJJdJBFKdpFEKNlFElHO/uwLAbwJ4FYAQwA2uvtvzewVAM8B6Ml+9GV3f4/dX1QjZDXfaO51f39/2JbNd2fzj6N522xuNKubsvEFrH3U9zlz5oRto7oswOvNFy5cCOPR+AQ2z7/o2Ilo33q21j+rw7PrhV3L0fXIxnxUum58OYNqrgF40d13mVkLgJ1mtjWL/cbd/7WM+xCROitnf/YuAF3Z1wNmth/A/LHumIhU13f6n93MvgdgJYBPspteMLM9ZvaGmc3IabPBzDrMrKOnp2e0HxGRGig72c1sCoA/A/iFu/cD+B2A7wNYgeFX/l+N1s7dN7p7yd1L7P9HERk7ZSW7mTVhONH/4O5/AQB3P+3u1919CMAmAKvHrpsiUhRNdht+e+91APvd/dcjbp834sd+DGBv9bsnItVSzrvxbQB+CuBzM9ud3fYygHVmtgKAAzgK4GfsjswsnBLJyhXRtEA21fLgwYNhfOnSpWE8KpUU2Z4X4FM1Ozs7w3j07xE7NpuKyR4bK2myMlKEXQ+sNBedF1buZEtFs+2kmaikyR539JxGsXLejf8QwGjFO1pTF5HGoRF0IolQsoskQskukgglu0gilOwiiVCyiySipktJDw0NhcsDs2mFUS2dTcWcPz+eu8PqydHUQTZVkz0uNj13wYIFYTyabnnq1KmwLRvCzOrobLnnIksmszEA7Lyy5zTCHnd3d3cYZ1Noo76xMSNs7EQevbKLJELJLpIIJbtIIpTsIolQsoskQskukgglu0girNKaXUUHM+sB8NWIm2YD6K1ZB76bRu1bo/YLUN8qVc2+/Y27jzp4oqbJ/q2Dm3W4e6luHQg0at8atV+A+lapWvVNf8aLJELJLpKIeif7xjofP9KofWvUfgHqW6Vq0re6/s8uIrVT71d2EakRJbtIIuqS7Gb2kJkdMLPDZvZSPfqQx8yOmtnnZrbbzDrq3Jc3zKzbzPaOuG2mmW01s0PZ51H32KtT314xs87s3O02s0fq1LeFZvaBme03s31m9vPs9rqeu6BfNTlvNf+f3czGAzgI4B8AnACwA8A6d/+fmnYkh5kdBVBy97oPwDCzHwAYBPCmu9+V3fYvAPrc/dXsF+UMd/+nBunbKwAG672Nd7Zb0byR24wDeBzAP6KO5y7o11OowXmrxyv7agCH3f1Ld78C4I8AHqtDPxqeu28H0HfDzY8B2Jx9vRnDF0vN5fStIbh7l7vvyr4eAPDNNuN1PXdBv2qiHsk+H8DxEd+fQGPt9+4AtpjZTjPbUO/OjGKuu3cBwxcPgFvq3J8b0W28a+mGbcYb5txVsv15UfVI9tEWc2uk+l+bu68C8DCA57M/V6U8ZW3jXSujbDPeECrd/ryoeiT7CQALR3y/AMDJOvRjVO5+MvvcDeBtNN5W1Ke/2UE3+xyvfFhDjbSN92jbjKMBzl09tz+vR7LvALDUzFrNrBnATwC8W4d+fIuZTc7eOIGZTQbwIzTeVtTvAliffb0ewDt17MtfaZRtvPO2GUedz13dtz9395p/AHgEw+/IfwHgn+vRh5x+LQbwWfaxr959A/AWhv+su4rhv4ieATALQDuAQ9nnmQ3Ut/8A8DmAPRhOrHl16tu9GP7XcA+A3dnHI/U+d0G/anLeNFxWJBEaQSeSCCW7SCKU7CKJULKLJELJLpIIJbtIIpTsIon4X6fjA4WFkGYMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAANoklEQVR4nO3dX4hc93nG8efRPwuUGKRqZYRjqjT4oqZQJQyi4BJcQoPtGzkXKdFFcMFYsbEhhkBr0ov4UrRNQy+KbKUSUUvqEEhMfGHaGBEwuQkeG9WWI1q7Rk0UC2lVXcgxyPJKby/22GzknXNG5//u+/3AMLPnzMx5NTuPzuy853d+jggBWP82DF0AgH4QdiAJwg4kQdiBJAg7kMSmPje2c+fO2LNnT5+bnFtVV8J2T5UA9Z05c0YXL15c9c3aKOy275X0j5I2SvrniDhUdv89e/ZoOp022WRnlpaWStdv3Lhx5rqh/yO4du3azHVldWP9mUwmM9fV/hhve6Okf5J0n6S7JB2wfVfd5wPQrSZ/s++T9FZEvB0RVyX9QNL+dsoC0LYmYb9d0q9X/Hy2WPY7bB+0PbU9XVxcbLA5AE00Cftqf6h+7FuuiDgSEZOImCwsLDTYHIAmmoT9rKQ7Vvz8KUnvNCsHQFeahP1lSXfa/rTtLZK+Iun5dsoC0LbarbeIWLL9uKT/0HLr7VhEvNFaZT3btKnXQw5aNWR7bcjjE8pajhJtxxs1eodHxAuSXmipFgAd4nBZIAnCDiRB2IEkCDuQBGEHkiDsQBJrt7mMuXTdix5yeO/QQ4vXGvbsQBKEHUiCsANJEHYgCcIOJEHYgSRova1zXQ/zPHz4cKfP38Sjjz46dAmjwp4dSIKwA0kQdiAJwg4kQdiBJAg7kARhB5Kgz74OXL9+vfZjn3nmmRYrGZcmxwCsxx49e3YgCcIOJEHYgSQIO5AEYQeSIOxAEoQdSII++zqwtLQ0c93Ro0d7rOTmVPWyxzxWfi1qFHbbZyS9K+mapKWImLRRFID2tbFn/7OIuNjC8wDoEH+zA0k0DXtI+qntV2wfXO0Otg/antqeLi4uNtwcgLqahv3uiPicpPskPWb78zfeISKORMQkIiYLCwsNNwegrkZhj4h3iusLkp6TtK+NogC0r3bYbW+z/ckPb0v6oqRTbRUGoF1Nvo2/TdJzxbS5myT9W0T8e9WDyqYQrpqCd8OG+h9Erl69Wrp+y5YttZ97aGOuvcm48LXch686x0CT93JdtcMeEW9L+uMWawHQIVpvQBKEHUiCsANJEHYgCcIOJNH7ENeupxCepcv2VNM2S0SUrq9qSQ7ZghrylMtjPt3zEK21KuOrCEAnCDuQBGEHkiDsQBKEHUiCsANJEHYgiVGdSnqtDkPtuqc6ZB/9kUceGWzba1nTYye6wJ4dSIKwA0kQdiAJwg4kQdiBJAg7kARhB5Lotc8eEbpy5crM9V320ctOYS0NN85ekp5++unBtj3mMeFr2RB99Crs2YEkCDuQBGEHkiDsQBKEHUiCsANJEHYgiV777La1devWPjf5kSH76BhG2bEVGd8PlXt228dsX7B9asWyHbZftP1mcb292zIBNDXPx/jvSbr3hmVPSjoREXdKOlH8DGDEKsMeES9JunTD4v2Sjhe3j0t6oN2yALSt7hd0t0XEOUkqrnfNuqPtg7antqeLi4s1Nwegqc6/jY+IIxExiYjJwsJC15sDMEPdsJ+3vVuSiusL7ZUEoAt1w/68pAeL2w9K+kk75QDoSmWf3fazku6RtNP2WUnfknRI0g9tPyTpV5K+3GWRY/Dee+/NXLdt27bSx2adP31oGXvpZSrDHhEHZqz6Qsu1AOgQh8sCSRB2IAnCDiRB2IEkCDuQRO9TNpcNO6w6/e7S0tLMdVVtlqZtmKr2GtCW69evl66vO0U4e3YgCcIOJEHYgSQIO5AEYQeSIOxAEoQdSKL3PnuTfnfZlM5lU0E33e7YlQ1jbTpVdUSUrh/j1MRrXVUfvep3MvN5az0KwJpD2IEkCDuQBGEHkiDsQBKEHUiCsANJ9N5n70rVVNBdjRGWqnvZQ2p6fEHV6zbk8QtVr3vZ77Tq+IAPPvigdP2mTeXR6fL4g6rfySzs2YEkCDuQBGEHkiDsQBKEHUiCsANJEHYgiXXTZ69SdwzwPNbzWPkx/9uqjo0o63VXvR82b95cq6Y+1P2dVO7ZbR+zfcH2qRXLnrL9G9sni8v9tbYOoDfzfIz/nqR7V1n+nYjYW1xeaLcsAG2rDHtEvCTpUg+1AOhQky/oHrf9WvExf/usO9k+aHtqe7q4uNhgcwCaqBv2w5I+I2mvpHOSvj3rjhFxJCImETFZWFiouTkATdUKe0Scj4hrEXFd0ncl7Wu3LABtqxV227tX/PglSadm3RfAOFT22W0/K+keSTttn5X0LUn32N4rKSSdkfS17kpsR5f94sOHD3f23GM35n972fn0hz7ffdmY9CbnVihTGfaIOLDK4qMd1AKgQxwuCyRB2IEkCDuQBGEHkiDsQBK9D3EdouWw3pUN16xqMY25ddZU2b/t4YcfLn1s1amiq1QNoR3ivU66gCQIO5AEYQeSIOxAEoQdSIKwA0kQdiAJd3mK5RtNJpOYTqe9ba9N9LJzKRseO2aTyUTT6XTVNyR7diAJwg4kQdiBJAg7kARhB5Ig7EAShB1IIs2UzVndcsstpevff//9nirB0NizA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EAS9NkLFy5cKF2/a9eunippV9M+etNx3Yzlv3lV55ioO9105Z7d9h22f2b7tO03bH+9WL7D9ou23yyut9eqAEAv5vkYvyTpGxHxh5L+RNJjtu+S9KSkExFxp6QTxc8ARqoy7BFxLiJeLW6/K+m0pNsl7Zd0vLjbcUkPdFQjgBbc1Bd0tvdI+qykX0i6LSLOScv/IUha9Y9a2wdtT21PFxcXG5YLoK65w277E5J+JOmJiLg87+Mi4khETCJisrCwUKdGAC2YK+y2N2s56N+PiB8Xi8/b3l2s3y2p/OtsAIOqPJW0l7/nPy7pUkQ8sWL530n6v4g4ZPtJSTsi4q/KnqvqVNJXr14trWXLli2l65voqt0xD9pT/RvzqaLLpjWXyqd7LjuV9Dx99rslfVXS67ZPFsu+KemQpB/afkjSryR9eY7nAjCQyrBHxM8lzdqtfaHdcgB0hcNlgSQIO5AEYQeSIOxAEoQdSGJUQ1y77KOPWVXPlz786sbcK2+irI/e6Hk7eVYAo0PYgSQIO5AEYQeSIOxAEoQdSIKwA0mMqs9epWycb9PeZJfj1Ztar/3kK1eulK7funVrT5WMS5Px7KWPq/UoAGsOYQeSIOxAEoQdSIKwA0kQdiAJwg4ksab67F2N813rys63v2lT+a94yNc0ax+9CuPZATRC2IEkCDuQBGEHkiDsQBKEHUiCsANJVIbd9h22f2b7tO03bH+9WP6U7d/YPllc7u++XKxmw4YNtS/r2bVr12ZeMprnoJolSd+IiFdtf1LSK7ZfLNZ9JyL+vrvyALRlnvnZz0k6V9x+1/ZpSbd3XRiAdt3U5zjbeyR9VtIvikWP237N9jHb22c85qDtqe3p4uJis2oB1DZ32G1/QtKPJD0REZclHZb0GUl7tbzn//Zqj4uIIxExiYjJwsJC84oB1DJX2G1v1nLQvx8RP5akiDgfEdci4rqk70ra112ZAJqa59t4Szoq6XRE/MOK5btX3O1Lkk61Xx6Atszzbfzdkr4q6XXbJ4tl35R0wPZeSSHpjKSvdVBfay5fvly6/tZbb+2pkpvHKZfr2bhx49AljMo838b/XNJqJ1V/of1yAHRlfR9VAeAjhB1IgrADSRB2IAnCDiRB2IEk1tSppJvoso8eEaXrq6aDrpqit0kfvWo4J73oerqaVrlL46sIQCcIO5AEYQeSIOxAEoQdSIKwA0kQdiAJV/WIW92YvSjpf1cs2inpYm8F3Jyx1jbWuiRqq6vN2n4/IlY9/1uvYf/Yxu1pREwGK6DEWGsba10StdXVV218jAeSIOxAEkOH/cjA2y8z1trGWpdEbXX1Utugf7MD6M/Qe3YAPSHsQBKDhN32vbb/y/Zbtp8cooZZbJ+x/XoxDfV04FqO2b5g+9SKZTtsv2j7zeJ61Tn2BqptFNN4l0wzPuhrN/T0573/zW57o6T/lvTnks5KelnSgYj4Za+FzGD7jKRJRAx+AIbtz0v6raR/iYg/Kpb9raRLEXGo+I9ye0T89Uhqe0rSb4eexruYrWj3ymnGJT0g6S814GtXUtdfqIfXbYg9+z5Jb0XE2xFxVdIPJO0foI7Ri4iXJF26YfF+SceL28e1/Gbp3YzaRiEizkXEq8XtdyV9OM34oK9dSV29GCLst0v69Yqfz2pc872HpJ/afsX2waGLWcVtEXFOWn7zSNo1cD03qpzGu083TDM+mteuzvTnTQ0R9tVOyDam/t/dEfE5SfdJeqz4uIr5zDWNd19WmWZ8FOpOf97UEGE/K+mOFT9/StI7A9Sxqoh4p7i+IOk5jW8q6vMfzqBbXF8YuJ6PjGka79WmGdcIXrshpz8fIuwvS7rT9qdtb5H0FUnPD1DHx9jeVnxxItvbJH1R45uK+nlJDxa3H5T0kwFr+R1jmcZ71jTjGvi1G3z684jo/SLpfi1/I/8/kv5miBpm1PUHkv6zuLwxdG2SntXyx7oPtPyJ6CFJvyfphKQ3i+sdI6rtXyW9Luk1LQdr90C1/amW/zR8TdLJ4nL/0K9dSV29vG4cLgskwRF0QBKEHUiCsANJEHYgCcIOJEHYgSQIO5DE/wOVclbnjREXFgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAALD0lEQVR4nO3dT4ic9R3H8c+n1h78c0iaMYQYulZyqBQaZQiFFLFIJeYSPVjMQVIQ1oOCgoeKPegxlKr0UIS1BtNiFUHFHEJrCIJ4EUdJ86ehjZWtrlmyE3IwkoONfnvYJ2WNOzvjPM8zz5P9vl+wzOwzs5lvhrzz7MxvZh5HhACsft9pegAAk0HsQBLEDiRB7EASxA4k8d1J3ti6detiampqkjcJpDI7O6szZ854uctKxW57u6TfS7pC0h8jYs9K15+amlKv1ytzkwBW0O12B1429q/xtq+Q9AdJd0q6SdIu2zeN++cBqFeZx+xbJX0YER9FxBeSXpa0s5qxAFStTOwbJX2y5Pu5YtvX2J623bPd6/f7JW4OQBllYl/uSYBvvPY2ImYiohsR3U6nU+LmAJRRJvY5SZuWfH+9pFPlxgFQlzKxvydps+0bbH9P0r2S9lczFoCqjb30FhEXbD8k6W9aXHrbGxHHK5sMQKVKrbNHxAFJByqaBUCNeLkskASxA0kQO5AEsQNJEDuQBLEDSRA7kASxA0kQO5AEsQNJEDuQBLEDSRA7kMREP0oaq8/58+ebHmGgq666qukRWoU9O5AEsQNJEDuQBLEDSRA7kASxA0kQO5AE6+zJtXmdvKwyf7fVuEbPnh1IgtiBJIgdSILYgSSIHUiC2IEkiB1IgnX2Va7N6+jD1rLbPPvlqFTstmclnZP0paQLEdGtYigA1atiz/7ziDhTwZ8DoEY8ZgeSKBt7SHrT9vu2p5e7gu1p2z3bvX6/X/LmAIyrbOzbIuIWSXdKetD2rZdeISJmIqIbEd1Op1Py5gCMq1TsEXGqOF2Q9LqkrVUMBaB6Y8du+2rb1148L+kOSceqGgxAtco8G79e0uu2L/45f4mIv1YyFVaNMu8LZx2+WmPHHhEfSfpJhbMAqBFLb0ASxA4kQexAEsQOJEHsQBK8xXUVaHIJqsmPXF6NH/dcJ/bsQBLEDiRB7EASxA4kQexAEsQOJEHsQBKss18Gsq6jo1rs2YEkiB1IgtiBJIgdSILYgSSIHUiC2IEkWGdvAdbRMQns2YEkiB1IgtiBJIgdSILYgSSIHUiC2IEkiB1IYmjstvfaXrB9bMm2tbYP2j5ZnK6pd0wAZY2yZ39B0vZLtj0m6VBEbJZ0qPgeQIsNjT0i3pZ09pLNOyXtK87vk3RXtWMBqNq4j9nXR8S8JBWn1w26ou1p2z3bvX6/P+bNASir9ifoImImIroR0e10OnXfHIABxo39tO0NklScLlQ3EoA6jBv7fkm7i/O7Jb1RzTgA6jL0/ey2X5J0m6R1tuckPSFpj6RXbN8v6WNJ99Q55OWO96ujDYbGHhG7Blx0e8WzAKgRr6ADkiB2IAliB5IgdiAJYgeSIHYgCWIHkiB2IAliB5IgdiAJYgeSIHYgCWIHkuCQzasAb2PFKNizA0kQO5AEsQNJEDuQBLEDSRA7kASxA0kQO5AEsQNJEDuQBLEDSRA7kASxA0kQO5AEsQNJEDuQxNDYbe+1vWD72JJtT9r+1Pbh4mtHvWMCKGuUPfsLkrYvs/2ZiNhSfB2odiwAVRsae0S8LensBGYBUKMyj9kfsn2k+DV/zaAr2Z623bPd6/f7JW4OQBnjxv6spBslbZE0L+mpQVeMiJmI6EZEt9PpjHlzAMoaK/aIOB0RX0bEV5Kek7S12rEAVG2s2G1vWPLt3ZKODbougHYY+rnxtl+SdJukdbbnJD0h6TbbWySFpFlJD9Q3YvudP3++6REa0+a/O5+n/3VDY4+IXctsfr6GWQDUiFfQAUkQO5AEsQNJEDuQBLEDSXDI5uTavHRW1kp/t4zLcuzZgSSIHUiC2IEkiB1IgtiBJIgdSILYgSRYZ6/AsDXbuteyV/NaeV2G3WercR2ePTuQBLEDSRA7kASxA0kQO5AEsQNJEDuQBLEDSRA7kASxA0kQO5AEsQNJEDuQBLEDSRA7kATvZ8eKyr6vm/fat8fQPbvtTbbfsn3C9nHbDxfb19o+aPtkcbqm/nEBjGuUX+MvSHo0In4k6aeSHrR9k6THJB2KiM2SDhXfA2ipobFHxHxEfFCcPyfphKSNknZK2ldcbZ+ku2qaEUAFvtUTdLanJN0s6V1J6yNiXlr8D0HSdQN+Ztp2z3av3++XHBfAuEaO3fY1kl6V9EhEfDbqz0XETER0I6Lb6XTGmRFABUaK3faVWgz9xYh4rdh82vaG4vINkhbqGRFAFYYuvdm2pOclnYiIp5dctF/Sbkl7itM3aplwFWj6o6bLaPNsZazGj4oeZpR19m2S7pN01PbhYtvjWoz8Fdv3S/pY0j21TAigEkNjj4h3JHnAxbdXOw6AuvByWSAJYgeSIHYgCWIHkiB2IAne4toCl/M6fJMyrpWXwZ4dSILYgSSIHUiC2IEkiB1IgtiBJIgdSIJ19ssA68moAnt2IAliB5IgdiAJYgeSIHYgCWIHkiB2IAliB5IgdiAJYgeSIHYgCWIHkiB2IAliB5IgdiCJobHb3mT7LdsnbB+3/XCx/Unbn9o+XHztqH9cAOMa5cMrLkh6NCI+sH2tpPdtHywueyYiflffeACqMsrx2eclzRfnz9k+IWlj3YMBqNa3esxue0rSzZLeLTY9ZPuI7b221wz4mWnbPdu9fr9fbloAYxs5dtvXSHpV0iMR8ZmkZyXdKGmLFvf8Ty33cxExExHdiOh2Op3yEwMYy0ix275Si6G/GBGvSVJEnI6ILyPiK0nPSdpa35gAyhrl2XhLel7SiYh4esn2DUuudrekY9WPB6Aqozwbv03SfZKO2j5cbHtc0i7bWySFpFlJD9QwH4CKjPJs/DuSvMxFB6ofB0BdeAUdkASxA0kQO5AEsQNJEDuQBLEDSRA7kASxA0kQO5AEsQNJEDuQBLEDSRA7kASxA0k4IiZ3Y3Zf0n+WbFon6czEBvh22jpbW+eSmG1cVc72g4hY9vPfJhr7N27c7kVEt7EBVtDW2do6l8Rs45rUbPwaDyRB7EASTcc+0/Dtr6Sts7V1LonZxjWR2Rp9zA5gcpreswOYEGIHkmgkdtvbbf/T9oe2H2tihkFsz9o+WhyGutfwLHttL9g+tmTbWtsHbZ8sTpc9xl5Ds7XiMN4rHGa80fuu6cOfT/wxu+0rJP1L0i8kzUl6T9KuiPjHRAcZwPaspG5ENP4CDNu3Svpc0p8i4sfFtt9KOhsRe4r/KNdExK9bMtuTkj5v+jDexdGKNiw9zLikuyT9Sg3edyvM9UtN4H5rYs++VdKHEfFRRHwh6WVJOxuYo/Ui4m1JZy/ZvFPSvuL8Pi3+Y5m4AbO1QkTMR8QHxflzki4eZrzR+26FuSaiidg3Svpkyfdzatfx3kPSm7bftz3d9DDLWB8R89LiPx5J1zU8z6WGHsZ7ki45zHhr7rtxDn9eVhOxL3coqTat/22LiFsk3SnpweLXVYxmpMN4T8oyhxlvhXEPf15WE7HPSdq05PvrJZ1qYI5lRcSp4nRB0utq36GoT188gm5xutDwPP/XpsN4L3eYcbXgvmvy8OdNxP6epM22b7D9PUn3StrfwBzfYPvq4okT2b5a0h1q36Go90vaXZzfLemNBmf5mrYcxnvQYcbV8H3X+OHPI2LiX5J2aPEZ+X9L+k0TMwyY64eS/l58HW96NkkvafHXuv9q8Tei+yV9X9IhSSeL07Utmu3Pko5KOqLFsDY0NNvPtPjQ8Iikw8XXjqbvuxXmmsj9xstlgSR4BR2QBLEDSRA7kASxA0kQO5AEsQNJEDuQxP8AIb2CsIWm+a8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(1) \n",
    "I = image_six.flatten() # on vectorialise l'image d chiffre 1\n",
    "y = I+0.10*np.random.randn(784)\n",
    "M=max(np.abs(y))\n",
    "\n",
    "# Représentation graphique de y\n",
    "\n",
    "plt.imshow(np.abs(y).reshape(28, 28), cmap=\"Greys\", vmin=0,vmax=M)\n",
    "plt.show()\n",
    "plt.close()\n",
    "\n",
    "# On éxécute SLOPE avec le l'hyper-paramètre quasi-sphérique \n",
    "\n",
    "Lambda=np.sqrt(range(1,785))-np.sqrt(range(0,784))\n",
    "\n",
    "# Représentation graphique de l'opérateur proximal lorsque tau = 4,77\n",
    "\n",
    "tau=4.77\n",
    "image_six_debruite = prox_L1_ordonnee(y,tau*Lambda)\n",
    "\n",
    "plt.imshow(np.abs(image_six_debruite).reshape(28, 28), cmap=\"Greys\",vmin=0,vmax=M)\n",
    "plt.show()\n",
    "plt.close()\n",
    "\n",
    "\n",
    "# Représentation graphique de l'opérateur proximal lorsque tau = 10,80\n",
    "\n",
    "tau=10.80\n",
    "image_six_debruite = prox_L1_ordonnee(y,tau*Lambda)\n",
    "\n",
    "plt.imshow(np.abs(image_six_debruite).reshape(28, 28), cmap=\"Greys\",vmin=0,vmax=M)\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparaison entre la condition d'irreprésentabilité et condition d'accessibilité"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Condition d'irreprésentabilité du SLOPE :** La condition d'irreprésentabilité est satisfaite pour le schéma $m\\in P^{\\rm slope}_p$ si\n",
    "$$X^T\\tilde X_m^+ \\tilde \\lambda_m \\in \\partial J_\\lambda(m) \\Leftrightarrow \\begin{cases} J_\\lambda^*(X^T\\tilde X_m^+ \\tilde \\lambda_m)\\le 1 \\\\\n",
    "m^TX^T\\tilde X_m^+ \\tilde \\lambda_m=J_\\lambda(m)\n",
    "\\end{cases}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Norme $\\ell_1$ ordonnée et norme $\\ell_1$ ordonnée duale :** C'est normes sont utiles pour tester la condition d'accessibilité \n",
    "$$J_\\lambda(b)=\\sum_{i=1}^{p}\\lambda_i|b|_{\\downarrow i}\\; \\text{ et }\\;J^*_\\lambda(b)=\\max\\left\\{ \\frac{|b|_{\\downarrow 1}}{\\lambda_1}, \\frac{|b|_{\\downarrow 1}+|b|_{\\downarrow 1}}{\\lambda_1+\\lambda_2},\\dots,\\frac{|b|_{\\downarrow 1}+\\dots+|b|_{\\downarrow p}}{\\lambda_1+\\dots+\\lambda_p}\n",
    "\\right\\}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def norme(b,Lambda):\n",
    "    b = - np.sort(-np.abs(b))  # b ordonné par valeur absolue décroissante \n",
    "    return np.sum(Lambda*b)\n",
    "\n",
    "def norme_duale(v,Lambda):\n",
    "    v = - np.sort(-np.abs(v)) # v ordonné par valeur absolue décroissante\n",
    "    v_cum = np.cumsum(v)  # somme cumulé de v\n",
    "    Lambda_cum = np.cumsum(Lambda) # somme cumulé de Lambda\n",
    "    return np.max(v_cum/Lambda_cum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxlUlEQVR4nO3deXwU9f3H8dcnF/d9ySE3iCigGEHRoq231qse9apXrdVqtfprq9YeWmtbbbXVelBrvYqKt0VrRaviLRCQG4FwhzPcEAJJNp/fHzPBNW6SDWSzm+z7+XjsY+f4zsxnJpv5zHxn5jvm7oiISPrKSHYAIiKSXEoEIiJpTolARCTNKRGIiKQ5JQIRkTSnRCAikuaUCCShzOxSM/soqn+7mfWtpvwcMzu6HuK6zczGJjOGeJjZEWY2xczaJzuWaGZ2tJkVVDN+jJn9KlbZeLevmWWY2Xgz+0FdxCxVy0p2AJIcZnYBcCMwCNgGTAfudPePqptub7l7y6gYngAK3P2XUeMPSOTyYzGz3sBEd++drBhiMbN9gd8Dp7j7xgQv6zagv7tfVBfzc/erqhm3e/vWsNw7gXfc/R91EZNUTYkgDZnZjcDNwFXABKAEOBE4HUhoImjozCzL3cui+g0wdy+v6/m7+wrgqLqYb0Pk7rckO4a04e76pNEHaANsB86ppkwT4K/AqvDzV6BJOO5ooAD4P2AdsBq4LGraDsB4YCswGbgD+ChqvAP9gSuBUoIktB14LRy/FDh2b+OIsU59gPcJzn7eBh4AxobjegNLo8pGx3Ab8CIwNlynK4CJBEerHwPF4foMCue7EZgPnBs1vyeAMeH4bWEcvSptk2uAhcCScNi3Cc7SNgOfAEOjyt8ErAznNR84JhyeQZDgFwEbgOeB9lHr6MAlwHJgPXBrOO7E8O9QGv4tZoTDLwPmhctZDPwwKoaK7f+LcF5LgQsrrfPvostW3r7VLLcN8M/wb7oS+B2Qmez/ncb8SXoA+tTzHzz45ysDsqop81vgM6Az0CncEd0Rjjs6nP63QDZwMrADaBeOHxfugFoAB4b/yF9LBGH37p1F1PilfLkT3uM4YqzTp8C9BMlldLhzG1tF2egYbgt3VGcQ7GibESSC5cABBGfVbYAV4Y4zCxge7hwPiFrPbeFymwD3xdgmbwPtw/kPJ0huI4FMgp330nDa/cJldQun7Q30C7t/Em6vHmHZvwPPRpVz4B/hMoYBu4D9o9ZzbKXtcArQDzCCM5MdwPBK279imx4FFAH7Vf7bUkUiqGa5r4axtyD4208mKgnpU/cfXSxOPx2A9R5VvRHDhcBv3X2duxcCtwPfixpfGo4vdfc3CI7m9jOzTOAs4NfuXuTus4En9yLWPYqj8kzMrCdwKPArd9/l7h8Ar9Uijk/d/VV3L3f34nDYE+4+J9yOJxKcUTzu7mXuPg14CTg7ah7/cfcP3H0XcCtweHgNoMIf3H1jOP8fAH9390nuHnH3Jwl22ocBEYId72Azy3b3pe6+KJzHDwmO8gvC5dwGnG1m0VXAt7t7sbvPAGYQJISY3P0/7r7IA+8DbwHfqFSsYpu+D/wHODeO7VklM+sCnAT8JPwNrQP+Apy3N/OV6ukaQfrZAHSsXNddSTdgWVT/snDY7nlUmnYH0JLgqD2L4Ig1eto9tadxxJrPJncvqjSvfWOUjWVFDcN6ASPNbHPUsCzgX7HKu/t2M9sYxrWi8vhwfpeY2Y+jhuUQnAW8b2Y/IdjJH2BmE4Ab3X1VON0rZhZ9vSICdInqXxPVXdX2AsDMTgJ+AwwkOBtqDsyKKhJrm0b/ffZEL4IzvNXB5RcIlx3rbyB1RGcE6edTYCdBVUdVKnYqFXqGw2pSSFBdEL2D7VlN+Zqavt3TOCpbDbQzsxZxxlVZrDijh60A3nf3tlGflu5+dVSZ3dvEzFoSVANFr0vl+d1ZaX7N3f1ZAHd/xt2PJNg2DtwVNd1JlaZr6u4ra7uOZtaE4Kzmz0AXd28LvEFQTVQh1jat7d+n8rZdQXD20zFqHVp7itzJ1VgpEaQZd98C/Bp40MzOMLPmZpZtZieZ2d1hsWeBX5pZJzPrGJaPec99pXlHgJeB28L5Diao367KWqDKZwr2NI4YcS0D8oDbzSzHzI4ETq3tfKrxOjDQzL4XbstsMzvUzPaPKnOymR1pZjkEF9AneXBXUCz/AK4ys5EWaGFmp5hZKzPbz8y+Fe6odxJcrI6E040B7jSzXgDhdjs9znVYC/Q2s4p9Qg5BFVQhUBaeHRwfY7qKbfoNggvcL8S5vJjLdffVBFVQ95hZ6/BZgn5mlrZ3T9UHJYI05O73EjxD8EuCf/QVwLUEF+kguEsjD5hJUBUwLRwWj2sJqhvWEFwwfLyasv8kqOvebGavxhi/N3FUdgHBxdeNBNUdT+3hfL7G3bcR7CTPIzgiXkNwlN4kqtgz4XI3AocQXP+oan55BNcJHgA2AfnApeHoJsAfCS5GryG4mPqLcNx9BHdsvWVm2wguHI+MczUqduAbzGxauE7XEVz430Sw/cZXmmZNOG4V8DRwlbt/EefyYi437L6YIBHNDef/ItC1lvOVWjB3vZhGJJFiPTgnkkp0RiAikuaUCERE0pyqhkRE0pzOCERE0lyDe6CsY8eO3rt372SHISLSoEydOnW9u3eKNa7BJYLevXuTl5eX7DBERBoUM6vyKX9VDYmIpLmEJQIze8zM1pnZ7CrGm5ndb2b5ZjbTzIYnKhYREalaIs8IniBolbEqJwEDws+VwMMJjEVERKqQsEQQNvVb3ev1TgeeCpu4/Qxoa2Z6jFxEpJ4l8xpBd77atGxBOExEROpRMhOBxRgW8+k2M7vSzPLMLK+wsDDBYYmIpJdkJoICvtpufQ+qaMvc3R9x91x3z+3UKeZtsCIisoeS+RzBeOBaMxtH0FTulrAtchGRRsPd2VJcysrNxazdupOSsnIi5VBWXk65O2URp9ydSDlEysuJlDsRr+iO+nYnt1c7Rg+s+4PhhCUCM3uW4KXVHc2sgKAt9mwAdx9D8LajkwnaWt9B8OJvEZEGpSxSzrptu1i5uZhVm4sp2BR8r9xczMqwu6gkUvOM4nDVUf0aViJw9/NrGO/ANYlavojI3thZGmFjUQkbi0rYtKNkd/f67btYtXknKzcFO/s1W3cSKf/q5c12zbPp1rYZfTq24Ij+HenRrhnd2jZjnzZNaZqVSVamkWFGZoaRlWFkVHyHw3Z/ovozDKLe41ynGlwTEyIie6poVxmLC4tYuXkHG4tK2Vi0i41FpV/Z0Vfs+HdUcRSfmWHs07op3ds1Y0Sf9nRvG+zku7drRve2TenWthnNcxrWrrVhRSsiUgN3Z0NRCfnrtpO/bjuLCsPvddtZtWXn18q3yMmkXYscOrTIoUPLHAZ0bkn7Fjm0a5FD+6hPu+ZBmdbNssnMSMyRebIoEYhIgxQpd1ZuKia/cFu4oy8iP9zpbyku3V2uWXYm/Tq3YESf9vTv3JJ+nVqyb/vmdGgZ7NybZmcmcS1SgxKBiDQIZZFyPl+xmQ8WFPLBwvV8sXoru8rKd4/v2DKHfp1acsrQrvTv1DLY6XduSdfWTcloZEfwdU2JQERS1oqNO/hgYSEfLCjkk/wNbNtVRobBwT3bcfHhvejfueXuo/y2zXOSHW6DpUQgIiljR0kZkxZv5P0FhXywsJDFhUUAdGvTlFOGdmX0wE4c0a8jbZpnJznSxkWJQESSxt2Zv3Yb788PdvxTlmyiJFJOk6wMDuvbgQtH9uKogR3p16llwm6dFCUCEalnJWXlfJRfyJuz1/D+gkLWbt0FwMAuLblkVC9GD+zEob3b6yJuPVIiEJGE21UW4cMF63lj9mrenruWbTvLaNU0i9EDO3HUgE58Y2BHurZpluww05YSgYgkxM7SCB8sKOSNWat5Z946tu0qo3XTLE44YB9OGdKVI/p3JCdLb8tNBUoEIlJndpZGmDi/Yue/lqKSCG2bZ3PSkH04eUhXRvXTzj8VKRGIyF4pLokwcf46/jNrNe9+sY4dJRHaNc/m1GHdOHlIVw7v14HsTO38U5kSgYjskanLNvHYx0t4d946iksjtG+Rw+kHdeeUIV05rG97srTzbzCUCESkVj5bvIG/vbuQj/M30LZ5Nt8ZHuz8R/TRzr+hUiIQkRq5Ox/nb+D+dxcyeclGOrZswq0n78+Fh/VscC1tytfpLygiVXJ3Ji4o5P53FvL58s3s07opt506mPNG9NR9/o2IEoGIfI278/bctTzwXj4zC7bQvW0zfnfGgZyT24MmWUoAjY0SgYjsVl7uvDlnDX97N595q7fSs31z7jprCGce3EO3fTZiSgQiQqTceX3mKh54N5+F67bTt1ML7j13GKcN66YLwGlAiUAkjZVFynl1+ioeei+fxeuLGNilJX87/2BOHtK10b2FS6qmRCCSprbsKOWqsVP5dPEGBndtzZiLhnP84H30Epc0pEQgkoaWri/i8iemULCpmLvPGso5uT3UzHMaUyIQSTOfLd7AVWOnYsDTPxjJob3bJzskSTIlApE08kLeCn7xyix6tm/OY5ceSq8OLZIdkqQAJQKRNFBe7vzprfk8PHERR/bvyIMXDqdNM73uUQJKBCKNXHFJhBuem86bc9Zwwcie3H7aAWoNVL5CiUCkEVu7dSdXPJnH7FVb+NW3B3P5Eb11UVi+RolApJGas2oLVzyZx5biUh69OJdj9u+S7JAkRSkRiDRCb89dy/XjPqdts2xevGoUg7u1TnZIksKUCEQaEXfn0Q+X8Pv/zmNo9zb84+JcOrdumuywJMUpEYg0EqWRcn797zk8O3k5Jw/Zh3vOOYhmOWopVGqmRCDSCGzZUcqPnpnKx/kbuPab/bnxuIFqKkLiltB7yMzsRDObb2b5ZnZzjPFtzOw1M5thZnPM7LJExiPSGC1dX8SZD3/M5CUbueecYfz0hP2UBKRWEnZGYGaZwIPAcUABMMXMxrv73Khi1wBz3f1UM+sEzDezp929JFFxiTQW7s4LeQXc/toccrIyePqKwxjRR81FSO0lsmpoBJDv7osBzGwccDoQnQgcaGXBjc0tgY1AWQJjEmkU1m3byS9ensX/5q3jsL7t+fM5w+jRrnmyw5IGKpGJoDuwIqq/ABhZqcwDwHhgFdAK+K67l1eekZldCVwJ0LNnz4QEK9JQ/HfWan7xyix2lET41bcHc9mo3qoKkr2SyEQQ65fplfpPAKYD3wL6AW+b2YfuvvUrE7k/AjwCkJubW3keImlhS3Ept4+fw8ufr2RI9zb85bvD6N+5VbLDkkYgkYmgANg3qr8HwZF/tMuAP7q7A/lmtgQYBExOYFwiDc5HC9fzsxdnsG7bLq4/ZgDXfqu/2guSOpPIRDAFGGBmfYCVwHnABZXKLAeOAT40sy7AfsDiBMYk0qAUl0S4680veOKTpfTt1IKXrx7FsH3bJjssaWQSlgjcvczMrgUmAJnAY+4+x8yuCsePAe4AnjCzWQRVSTe5+/pExSTSkExfsZkbn5vO4vVFXHZEb246cRBNs/WAmNS9hD5Q5u5vAG9UGjYmqnsVcHwiYxBpaEoj5fztnYU8OHERXVo14ZkrRjKqf8dkhyWNmJ4sFkkhC9Zu48bnpzN75VbOGt6D35w2mNZN9QIZSSwlApEUUF7uPPbxEu6eMJ+WTbIYc9EhnHjgPskOS9KEEoFIkq3btpPrn53Op4s3cOz+XfjDd4bQqVWTZIclaUSJQCSJ8pZu5EdPT2PrzlLuPmso5+T20BvEpN4pEYgkgbvzxCdLufM/8+jRrhlPXj6C/bvq5TGSHEoEIvVsR0kZN780i/EzVnHs/l2459xhtGmmC8KSPEoEIvVoceF2rh47jYXrtvGzE/bj6qP6qZ0gSTolApF6MmHOGn76/AyyMo0nLx/BNwZ0SnZIIoASgUjClUXKueftBTw8cRFDe7Th4YsOoXvbZskOS2Q3JQKRBNqwfRfXjfucj/M3cP6Invzm1MFqJkJSjhKBSIJMX7GZq8dOZUNRCXefPZRzc/eteSKRJFAiEKlj7s7Tk5bz29fm0rl1E16+ehQHdm+T7LBEqqREIFKHdpZGuPWV2bw0rYCjBnbivvMOom3znGSHJVItJQKROrJ8ww6uGjuVeWu2cv0xA7j+mAG6NVQaBCUCkTqwYO02zhnzKe7OY5ccyjcHdU52SCJxUyIQ2UtbdpRy5VN5ZGdm8NLVh9OrQ4tkhyRSK3rpqcheiJQ71437nJWbixlz0XAlAWmQdEYgshfueWs+7y8o5PdnDiG3d/tkhyOyR3RGILKHXp+5iocmLuL8ET25YGTPZIcjsseUCET2wLzVW/nZCzM5pFc7bjttcLLDEdkrSgQitbSpqIQr/5VH62ZZPHzhcJpkqckIadh0jUCkFsoi5fz42c9Zu2UXz/3wMDq3bprskET2mhKBSC3c9eYXfJS/nrvPHsrBPdslOxyROqGqIZE4vfr5Sv7x4RIuObyXGpCTRkWJQCQOs1du4aaXZjKiT3t++W1dHJbGRYlApAbrt+/iyqfy6NAih4cuHE52pv5tpHHRNQKRapRGyrnm6WlsKCrhxatG0bFlk2SHJFLnlAhEqnHnf+YxaclG7j13GEN66J0C0jjpHFekCi/kreCJT5by/SP78J3hPZIdjkjCKBGIxDB9xWZufXU2o/p14JaTBiU7HJGEUiIQqWTdtp1c9a+pdG7VhAcuGE6WLg5LI5fQX7iZnWhm880s38xurqLM0WY23czmmNn7iYxHpCYlZeX8aOw0NheX8Mj3cmnfQq+ZlMYv7ovFZnYkMMDdHzezTkBLd19STflM4EHgOKAAmGJm4919blSZtsBDwInuvtzM9FonSarbX5tD3rJN3H/+wQzu1jrZ4YjUiyrPCMzsgKju3wA3AbeEg7KBsTXMewSQ7+6L3b0EGAecXqnMBcDL7r4cwN3X1S58kbrz7+kreXrScn54VF9OG9Yt2eGI1JvqqoZ6mdkfw+4zgdOAIgB3XwW0qmHe3YEVUf0F4bBoA4F2ZjbRzKaa2cWxZmRmV5pZnpnlFRYW1rBYkdrbsqOUO16fy7B92/LzE3RxWNJLlVVD7v6GmUXC3hJ3dzNzADOL5318Fmu2MZZ/CHAM0Az41Mw+c/cFlWJ5BHgEIDc3t/I8RPbaH9/8gk07Snny8gPJzIj10xVpvKq9WOzuE8LO583s70BbM/sB8D/g0RrmXQBEt8zVA1gVo8yb7l7k7uuBD4Bh8QYvUhemLtvIs5OXc9mo3hzQTQ+NSfqJ664hd/8z8CLwErAf8Gt3v7+GyaYAA8ysj5nlAOcB4yuV+TfwDTPLMrPmwEhgXm1WQGRvlEbKufWV2XRt05QbjhuY7HBEkiKuu4bM7C53vwl4O8awmNy9zMyuBSYAmcBj7j7HzK4Kx49x93lm9iYwEygHHnX32XuxPiK18thHS/hizTb+/r1DaNFELa5IejL3mqvczWyauw+vNGymuw9NWGRVyM3N9by8vPperDRCBZt2cNy9H3BE/448eklussMRSSgzm+ruMX/o1R4CmdnVwI+AvmY2M2pUK+DjugtRpH65O7eNnwOgl89L2qvpXPgZ4L/AH4DoJ4O3ufvGhEUlkmAT5qzlf/PW8YuTB9GjXfNkhyOSVDUlAnf3pWZ2TeURZtZeyUAaou27yrj9tTkM2qcVlx3RJ9nhiCRdPGcE3wamEjwDEH2DtQN9ExSXSML85e0FrNm6kwcu0NvGRKCGRODu3w6/ddgkjcLslVt4/OMlnD+iJ4f0apfscERSQk0Xi4dXN97dp9VtOCKJEyl3bn1lFu1b5HCTmpEQ2a2mqqF7qhnnwLfqMBaRhHpm0jJmFGzhr989iDbNs5MdjkjKqKlq6Jv1FYhIIq3bupO735zPkf07cvpBallUJFpNVUPfcvd3zew7sca7+8uJCUukbt3xn3nsipRzxxkHYqZG5USi1VQ1dBTwLnBqjHEOKBFIyvtgQSGvzVjFT44dQJ+O8TScK5Jeaqoa+k34fVn9hCNSt3aWRvjlq7Pp27EFVx/dL9nhiKSkuG6iNrMOZna/mU0LXyBzn5l1SHRwInvrwffyWb5xB78740CaZGUmOxyRlBTv0zTjgELgLODssPu5RAUlUhfy121jzPuLOPPg7ozq3zHZ4YikrHjb3W3v7ndE9f/OzM5IQDwidcLdufWV2TTPyeLWU/ZPdjgiKS3eM4L3zOw8M8sIP+cC/0lkYCJ746VpK5m0ZCM3nzSIji2bJDsckZRW0+2j2/iyjaEbgbHhqAxgO/CbhEYnsgc2FZXw+zfmcUivdnw3d9+aJxBJczXdNdSqvgIRqSt/+O88thaXcueZB5KhF9GL1Cjud/OZWTtgANC0Ypi7f5CIoET21KTFG3g+r4AfHtWXQfu0TnY4Ig1CvO8svgK4HugBTAcOAz5FbQ1JClmxcQfXPvs5+7ZvxvXHDEh2OCINRrwXi68HDgWWhe0PHUxwC6lISthUVMIlj09mV2mExy45lOY5ehG9SLziTQQ73X0ngJk1cfcvgP0SF5ZI/HaWRvj+k1Mo2FTMo5ccyoAuurQlUhvxHjYVmFlb4FXgbTPbBKxKVFAi8YqUOz9+9nM+X7GZhy4Yzog+7ZMdkkiDE1cicPczw87bzOw9oA3wZsKiEomDu/Ob8bN5e+5abjt1MCcN6ZrskEQapNrcNTQcOJLguYKP3b0kYVGJxOGhiYsY+9lyfnhUXy7VS+hF9li8jc79GngS6AB0BB43s18mMjCR6rw4tYA/TZjPGQd102snRfZSvGcE5wMHR10w/iMwDfhdogITqcoHCwq5+aWZHNG/A3efPUwPjYnspXjvGlpK1INkQBNgUZ1HI1KD2Su3cPXYqQzo0ooxFx1CTla8P2ERqUpNbQ39jeCawC5gjpm9HfYfB3yU+PBEvrRi4w4ufXwKbZvn8MRlh9KqqV5AL1IXaqoaygu/pwKvRA2fmJBoRKqwsaiESx6bTGmknHFXjqRL66Y1TyQicamp0bknK7rNLAcYGPbOd/fSRAYmUqG4JMIVT06hYHMxz1wxkv6d9cCYSF2Kt62hownuGlpK0CT1vmZ2iRqdk0SLlDvXjQseGHv4wuHk9tYDYyJ1Ld4rbfcAx7v7Ue4+GjgB+EtNE5nZiWY238zyzezmasodamYRMzs7zngkDXz1gbEDOPFAPTAmkgjxJoJsd59f0ePuC4Bqr9SZWSbwIHASMBg438wGV1HuLmBCvEFLeqh4YOzqo/txyajeyQ5HpNGKNxFMNbN/mtnR4ecfBBeQqzMCyHf3xeFTyOOA02OU+zHwErAu7qil0at4YOzMg7vz8xPUvqFIIsWbCK4C5gDXETRJPTccVp3uwIqo/oJw2G5m1h04ExhT3YzM7EozyzOzvMJCtX7d2FU8MHZk/47cddZQzPTAmEgi1Xix2MwygKnufiBwby3mHeu/1yv1/xW4yd0j1f2zu/sjwCMAubm5lechjciOkjJufH46/Tu35OGLhuuBMZF6UGMicPdyM5thZj3dfXkt5l0ARL85vAdfb7o6FxgXJoGOwMlmVubur9ZiOdKIPP7xUtZvL+GRi3P1wJhIPYm3raGuBE8WTwaKKga6+2nVTDMFGGBmfYCVwHnABdEF3H13k5Fm9gTwupJA+tpSXMrf31/Esft3ZnjPdskORyRtxJsIbq/tjN29zMyuJbgbKBN4zN3nmNlV4fhqrwtI+vnHB4vZurOMG4/TxWGR+lRTW0NNCS4K9wdmAf9097J4Z+7ubwBvVBoWMwG4+6Xxzlcan/Xbd/HYx0s4dVg3BndrnexwRNJKTVfiniSox59F8DzAPQmPSNLSQ+8tYldZOTccOyDZoYiknZqqhga7+xAAM/snMDnxIUm6WbW5mLGfLePs4T3o26llssMRSTs1nRHsbliuNlVCIrXxt3cXAnCdzgZEkqKmM4JhZrY17DagWdhvgLu7KnNlryxdX8TzeQV877BedG/bLNnhiKSlmpqhzqyvQCQ9/eV/C8jJzOCab/ZPdigiaUuPbUrSfLFmK+NnrOKyI3rTqVWTZIcjkraUCCRp7nlrAS2bZPHD0f2SHYpIWlMikKSYvmIzb89dyw9H96VNczUlIZJMSgSSFH+eMJ8OLXK47Ig+NRcWkYRSIpB698mi9XyUv54ffbM/LZrE28qJiCSKEoHUK3fnzxPm07VNUy4c2TPZ4YgISgRSz96bv45pyzdz3TEDaJqtu5NFUoESgdSb8nLnTxMW0KtDc84+pEeywxGRkBKB1Js3Zq9m3uqt3HjcQLIz9dMTSRX6b5R6URYp5963FrBfl1acOrRbssMRkShKBFIvXv58JYvXF/F/xw8kI0MvoxdJJUoEknC7yiLc97+FDNu3LccN7pLscESkEiUCSbhxk1ewcnMxPz1+IGY6GxBJNUoEklA7Ssr427v5HNa3PUf275jscEQkBiUCSagnP1nG+u27+NkJ++lsQCRFKRFIwmwpLmXM+4v41qDOHNKrfbLDEZEqKBFIwvzzw8VsKS7l/44fmOxQRKQaSgSSEBu27+KfHy3hlKFdOaBbm2SHIyLVUCKQhHh44iKKSyPccKzOBkRSnRKB1Ln5a7bx1GfL+M7wHvTv3DLZ4YhIDZQIpE7tLI3w42en0bppNjedOCjZ4YhIHPRWEKlTd7w+lwVrt/PU5SP0QnqRBkJnBFJn3py9mqcnLeeHo/syemCnZIcjInFSIpA6sXJzMT9/cSZDe7Th/47fL9nhiEgtKBHIXiuLlHPDuOlEyp37zzuYnCz9rEQaEl0jkL32wHv5TF66kb98dxi9O7ZIdjgiUksJPXQzsxPNbL6Z5ZvZzTHGX2hmM8PPJ2Y2LJHxSN2bvGQj97+zkO8c3J0zD9brJ0UaooQlAjPLBB4ETgIGA+eb2eBKxZYAR7n7UOAO4JFExSN1b/OOEn4y7nN6tm/Ob884MNnhiMgeSuQZwQgg390Xu3sJMA44PbqAu3/i7pvC3s8AHVI2EO7OzS/NonD7Lu4//2BaNlEto0hDlchE0B1YEdVfEA6ryveB/yYwHqlDT09azptz1vDzEwYxtEfbZIcjInshkYdxsRqf95gFzb5JkAiOrGL8lcCVAD179qyr+GQPzV+zjTten8vogZ34/pF9kh2OiOylRJ4RFAD7RvX3AFZVLmRmQ4FHgdPdfUOsGbn7I+6e6+65nTrpQaVkqmhColXTLO45Z5heRC/SCCQyEUwBBphZHzPLAc4DxkcXMLOewMvA99x9QQJjkTryu/8ETUjcc+5BakJCpJFIWNWQu5eZ2bXABCATeMzd55jZVeH4McCvgQ7AQ+FrDMvcPTdRMcneeXP2GsZ+tpwrR/flKDUhIdJomHvMavuUlZub63l5eckOI+2s2lzMSfd9SM/2zXnp6lF6elikgTGzqVUdaOu/WWoUKXd+Mm46ZZFy7j9fTUiINDa6+Vtq9MC7QRMS9547jD5qQkKk0dGhnVRr8pKN3PfOAs48uDvfGa7n/UQaIyUCqVJFExL7tm/Ob08/INnhiEiCqGpIYqpoQmLdtl28dPUoWjXNTnZIIpIgOiOQmP76v4W8OWcNPz1hP4bt2zbZ4YhIAikRyNf844PF3PfOQs45pAdXfqNvssMRkQRTIpCveGbScu58Yx6nDOnKH88aqiYkRNKAEoHs9u/pK7n11VkcvV8n/vLdg8hUEhBJC0oEAsDbc9dy4/MzGNG7PWMuOkQPjYmkEf23Cx/nr+eaZ6ZxYLfWPHpJLk2zM5MdkojUIyWCNDd12SZ+8FQefTq04InLRug2UZE0pESQxuau2splj0+mc6sm/Ov7I2jXIifZIYlIEigRpKlFhdu5+LFJtGiSxdgrRtK5ddNkhyQiSaJEkIYKNu3gokcnAfD0FSPp0a55kiMSkWRSIkgz67bu5MJHJ1G0q4ynLh9J304tkx2SiCSZ2hpKI5uKSvjePydTuG0XY68YyeBurZMdkoikACWCNLFtZymXPj6ZJRuKeOLSQxnes12yQxKRFKGqoTRQXBLh+0/mMWfVVh66YDij+ndMdkgikkJ0RtDIlZSVc/XTU5mydCP3nXcwxw7ukuyQRCTF6IygEYuUOzc8N52J8wv5w5lDOG1Yt2SHJCIpSGcEjdCmohJenFrAM5OXs2R9Eb88ZX/OG9Ez2WGJSIpSImgk3J3PV2xm7GfLeH3makrKysnt1Y6fnbAfJw/pmuzwRCSFKRE0cEW7yvj39FWM/WwZc1dvpUVOJufm9uCiw3oxaB/dHioiNVMiaKDmr9nG05OW8fK0lWzfVcagfVrxuzMO5IyDu9Oyif6sIhI/7TEakF1lEd6cvYaxny1jytJN5GRl8O0hXbnwsF4M79kWM71IRkRqT4mgAVixcQdPT1rOC3kr2FBUQq8OzbnlpEGck7sv7dViqIjsJSWCFBUpd977Yh1jJy3j/QWFGHDM/l246LBefKN/R71LWETqjBJBilm3bSfPT1nBs5NXsHJzMZ1bNeHH3+zPeSN60q1ts2SHJyKNkBJBCnB3Plu8kbGTljFh9hrKyp0j+nfgl6fsz7GDu5Cdqef+RCRxlAiSaEtxKS9PK+DpScvJX7edNs2yuWRUby4Y2ZN+ah5aROqJEkESzCrYwtjPljF+xiqKSyMM27ctfzp7KKcO66YXx4tIvUtoIjCzE4H7gEzgUXf/Y6XxFo4/GdgBXOru0xIZU7IUl0R4beYqnv5sGTMKttAsO5PTD+rGhSN7MaRHm2SHJyJpLGGJwMwygQeB44ACYIqZjXf3uVHFTgIGhJ+RwMPhd0pwd8o9uIMnUu7sLI1QXPEpiXzZX1LpuzTCzqjubTvLeO+LdWzdWUb/zi257dTBnDm8B22aZSd7FUVEEnpGMALId/fFAGY2DjgdiE4EpwNPubsDn5lZWzPr6u6r6zqYifPXccfrc7+yY4+UO2XlTrn7V4ZFovr3lBk0y86kWXYmTbMzGT2wExcd1ouRfdrrwS8RSSmJTATdgRVR/QV8/Wg/VpnuwFcSgZldCVwJ0LPnnrWi2appNoP2aU1GhpGVYWSYkZkBmRkZZGZAVkbG14ZlZmSQGTWsaXZGsHPPyfzKd9Oo7uY5QX+TrAzt8EWkQUhkIoi1F6x8iB1PGdz9EeARgNzc3D06TD+kVzsO6aXXM4qIVJbIG9QLgH2j+nsAq/agjIiIJFAiE8EUYICZ9TGzHOA8YHylMuOBiy1wGLAlEdcHRESkagmrGnL3MjO7FphAcPvoY+4+x8yuCsePAd4guHU0n+D20csSFY+IiMSW0OcI3P0Ngp199LAxUd0OXJPIGEREpHpqxEZEJM0pEYiIpDklAhGRNKdEICKS5iy4XttwmFkhsCzZcVTSEVif7CBqoSHF25BihYYVb0OKFRpWvKkYay937xRrRINLBKnIzPLcPTfZccSrIcXbkGKFhhVvQ4oVGla8DSlWUNWQiEjaUyIQEUlzSgR145FkB1BLDSnehhQrNKx4G1Ks0LDibUix6hqBiEi60xmBiEiaUyIQEUlzSgTVMLNMM/vczF4P+9ub2dtmtjD8bhdV9hYzyzez+WZ2QtTwQ8xsVjjufkvQa8vMbGm4nOlmlpfK8YavJH3RzL4ws3lmdngKx7pfuE0rPlvN7CcpHO8NZjbHzGab2bNm1jSFY70+jHOOmf0kHJYysZrZY2a2zsxmRw2rs/jMrImZPRcOn2Rmvesi7j3i7vpU8QFuBJ4BXg/77wZuDrtvBu4KuwcDM4AmQB9gEZAZjpsMHE7wNrb/AiclKNalQMdKw1IyXuBJ4IqwOwdom6qxVoo7E1gD9ErFeAle87oEaBb2Pw9cmqKxHgjMBpoTtIL8P2BAKsUKjAaGA7MT8T8F/AgYE3afBzyXyN9vteuarAWn+ofgbWnvAN/iy0QwH+gadncF5ofdtwC3RE07IfzDdwW+iBp+PvD3BMW7lK8ngpSLF2gd7qws1WONEfvxwMepGi9fvgO8PcHO9fUw5lSM9Rzg0aj+XwE/T7VYgd58NRHUWXwVZcLuLIInka2uYq/NR1VDVfsrwQ+zPGpYFw/foBZ+dw6HV/wDVigIh3UPuysPTwQH3jKzqWZ2ZQrH2xcoBB63oNrtUTNrkaKxVnYe8GzYnXLxuvtK4M/AcmA1wRv/3krFWAnOBkabWQcza07wgqp9UzTWaHUZ3+5p3L0M2AJ0SFjk1VAiiMHMvg2sc/ep8U4SY5hXMzwRjnD34cBJwDVmNrqassmMN4vgdPthdz8YKCI4xa5KKmxbLHjd6mnACzUVjTGsXuIN66tPJ6ia6Aa0MLOLqpukipgSHqu7zwPuAt4G3iSoVimrZpKU+B1UY0/iS5XYlQiqcARwmpktBcYB3zKzscBaM+sKEH6vC8sXEBzNVOgBrAqH94gxvM65+6rwex3wCjAiReMtAArcfVLY/yJBYkjFWKOdBExz97VhfyrGeyywxN0L3b0UeBkYlaKx4u7/dPfh7j4a2AgsTNVYo9RlfLunMbMsoA3Bdqh3SgQxuPst7t7D3XsTVAe86+4XAeOBS8JilwD/DrvHA+eFdwH0IbjoNTk8ddxmZoeFdwpcHDVNnTGzFmbWqqKboF54dirG6+5rgBVmtl846BhgbirGWsn5fFktVBFXqsW7HDjMzJqHyzgGmJeisWJmncPvnsB3CLZvSsYapS7ji57X2QT7meQ84ZuMCxMN6QMczZcXizsQXEBeGH63jyp3K8GdAvOJumsByCXYKS8CHiABF4MI6t1nhJ85wK0pHu9BQB4wE3gVaJeqsYbLaQ5sANpEDUvJeIHbgS/C5fyL4C6WVI31Q4KDgBnAMam2XQkS02qglODo/ft1GR/QlKCqMZ/gzqK+ifj9xvNRExMiImlOVUMiImlOiUBEJM0pEYiIpDklAhGRNKdEICKS5pQIGjkzczO7J6r/p2Z2Wx3Ne3tdzKcu5m1B66sdw+5PEhNVzTFEf9dQdqKZ9Y6z7KVm9kAt4jjawhZz60O4vFHVjD/NzKp7ejyeZWSGzaeMjhr2lpmdszfzlYASQeO3C/hOxU6yvoVPTNYrd69ypyQJcTTBE8xfY2ZZ7j7e3f+4Nwtw9whBa50Pmlm2mZ0fDPaamvyQOCgRNH5lBO9PvaHyCDPrZWbvmNnM8LtnOPwJM3vYzN4zs8VmdpQFbbPPM7MnKs3jHjObFk7fKRw20cx+b2bvA9db0B77++ER3YSKR/QrzaePmX1qZlPM7I5K434WDp9pZrfXtMIVZxNmlmFmD1nQ3v3rZvaGmZ0djltqZreHsc8ys0Hh8Bbhuk6xoFG808PhB5jZZAveSTDTzAbEWHRh9Hd4pDzRvnz3wtPh06UQNCUQiZqm8jpcZmYLwm14RNTwTmb2UhjfFDM7Itb0UeVHmNkn4bp8Yl8+0R1d5ujw7/N8uMw/mtmF4frOMrN+YblTLWg3/3Mz+5+ZdbGgDf2rgBvCbfON8Pdzr5m9B9wVfUYTTvOKmc0IP6PC4RdFbd+/m1lm5Tg9aJbkE+A24PfANdWtu9RCsp5k06d+PsB2gqaflxK0ZfJT4LZw3GvAJWH35cCrYfcTBG0sGUEjZluBIQQHDlOBg8JyDlwYdv8aeCDsngg8FHZnE/zzdgr7vws8FiPO8cDFYfc1wPaw+3iCRGbh8l8HRseYfilhM9xR054NvBFOtw+wCTg7qvyPw+4fETaJTLCDuSjsbgssAFoAf4ta1xzCNv9r2PZHE7Qo2SOM4VPgyDim60rQXESncFkfR23bZyrmAfQE5lWx3Iqn4VsDWWH3scBLVZTfHC63CbASuD0cdz3w17C7HV8+FXsFcE/YfRvw06j5PRH+nSra4780Kv7ngJ+E3ZkEv8n9CX6L2eHwhyp+CzFibU/QUOGdyf7fakyfej9tl/rn7lvN7CngOqA4atThBG28QNAcwd1R415zdzezWcBad58FYGZzCNpon07QRPdzYfmxBI2cVagYvh/BS0jeDg+GMwke26/sCOCsqFjuCruPDz+fh/0tCdpx+aCG1QY4EnjB3cuBNeERarSKeKfy5XY4nqDBwZ+G/U0JdrifAreaWQ/gZXdfGMfyIWhvpgDAzKYTbLuPaphmJDDR3SvOLJ4DBobjjgUGf3liQWsza+Xu26qYVxvgyfAMxgkScyxTPGxe2cwWAW+Fw2cB3wy7ewDPhWd0OQTvlajKCx5U51T2LYL2dgjHbzGz7wGHAFPC9WrGl425VTaaILkeWM2ypZaUCNLHX4FpwOPVlIlub2RX+F0e1V3RX9XvJnr6ovDbgDnufngcMcZq78SAP7j73+OYPta01alYrwhfrpMBZ7n7/Epl55nZJOAUYIKZXeHu78YRQ/S2i15OTapq+yWD4GUmxVWMr+wO4D13PzOsxplYRbnKf+Pov39FzH8D7nX38WZ2NMGZQFWKqhlXmQFPuvst1RYKGlS8myCZPGZmJ7v7G7VYjlRB1wjShLtvJHh14fejBn9C0LoqwIXUfKRaWQZB9QvABVVMPx/oZGaHA1hwoe+AGOU+rhRLhQnA5WbWMpy+u4WtVsbhI+Cs8FpBF4IqkJpMAH5cUZdvZgeH332Bxe5+P0E11tA4Y9gTk4CjLXhpSzbB27wqvAVcW9FjZgfVMK82BFU9EFTR7I3oeV0SNXwb0CrOebwDXA277wRqHQ47275sjbS9mfWKMe2vgefd/QuC6ry/mFnT2q+GVKZEkF7uAaLvHroOuMzMZgLfI6gPro0i4AAzm0pwlPbbygXcvYQgWdxlZjMIqpRi3WFyPcELdaYQ7HAqpn+LoF7807Ca6kXi3+m8RNBq5Gzg7wQ72C01THMHQfXJTAteWl5x4fq7wOywemcQ8FScMdRaWEVzG0F11P8IzuQqXAfkhhes5xJcqK3O3cAfzOxjgmq5vXEb8IKZfUjwWsUKrwFnVlwsrmEe1wPfDP+WU4ED3H0u8EuCN+zNJHhZzVduKDCzwcCZwJ0A7j6dIGnftJfrJKDWR6VxM7OW7r7dzDoQNPV7hAfvRBCRkK4RSGP3upm1Jbi4eYeSgMjX6YxARCTN6RqBiEiaUyIQEUlzSgQiImlOiUBEJM0pEYiIpLn/B8s+ohNBCDv4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "NB=1000 # nombre de simulations pour calculer la probabilité que la condition d'irreprésentabilité soit satisfaite\n",
    "k=len(ind)\n",
    "Lambda=np.sqrt(range(1,785))-np.sqrt(range(0,784))\n",
    "Lambdatilde=np.sum(Lambda[0:k])\n",
    "N_ci = np.linspace(3500,10500,21,dtype=int)\n",
    "succes_ci=np.zeros(len(N_ci))\n",
    "tol = 1e-6\n",
    "\n",
    "for i in range(len(N_ci)):\n",
    "    n = N_ci[i]\n",
    "    for j in range(NB):\n",
    "        X = np.random.randn(n,784)    \n",
    "        Xtilde = X @ I\n",
    "        CI = Lambdatilde * (X.T @ Xtilde) / (Xtilde.T @ Xtilde)\n",
    "        CI.T @ I - norme(I, Lambda)\n",
    "        if norme_duale(CI, Lambda)<1+tol and np.abs(CI.T @ I - norme(I, Lambda))<tol:\n",
    "            succes_ci[i] = succes_ci[i]+1/NB\n",
    "\n",
    "plt.ylabel('Probabilité')  \n",
    "plt.xlabel('Nombre de lignes \"n\" de la matrice X')  # Titre de l'axe x\n",
    "plt.plot(N_ci, succes_ci)\n",
    "plt.title('Condition d\\'irreprésentabilité')\n",
    "plt.show()            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**L'algorithme de Douglas Rachford :** Cet algorithme permet de vérifier la condition d'accessibilité et de calculer l'estimateur par récupération $J_\\lambda$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def DR_SLOPE(A,b,Lambda,rtol=1e-6,it_max=1e2):\n",
    "    p = A.shape[1]\n",
    "    y = np.zeros(p)   \n",
    "    Apsi = np.linalg.pinv(A)\n",
    "    it = 0\n",
    "    while True:\n",
    "        if it >= it_max:\n",
    "            #print(f\"attention : nombre maximal d'itérations atteint ({it_max :.0})\")\n",
    "            break\n",
    "        if it>=1:\n",
    "            Jvieux = norme(x,Lambda)\n",
    "        x = prox_L1_ordonnee_rapide_depistage(y,Lambda)\n",
    "        #x = prox_SLOPE(y,Lambda)\n",
    "        z = 2*x-y + Apsi @ (b - A @ (2*x-y))\n",
    "        y = y + z - x\n",
    "        if it>=1:\n",
    "            Jnouv = norme(x,Lambda)\n",
    "        if it>=1 and np.abs(Jnouv - Jvieux) < rtol*Jnouv:\n",
    "            #print(f\"critère d'arrêt (tolérance relative) satisfait en {it} itérations\")\n",
    "            break\n",
    "        it = it+1    \n",
    "    return(x) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Condition d'accessibilité :** La condition d'accessibilité est satisfaite pour le schéma $m\\in P^{\\rm slope}_p$ si $Xm= Xb$ implique  $J_\\lambda(b)\\ge J_\\lambda(m)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvGUlEQVR4nO3dd3wc1bn/8c+jZsndWO4dF4wNxsamGDCYXkINhAAJLYCBYOCSSy6Q8KNcEgJJIBAglzjUQIwxwdSYmtBNce8FY1xk494tyyr7/P6YEV6EykrWarW73/frta+dnTkz+5xdaZ6dMzPnmLsjIiLpKyPRAYiISGIpEYiIpDklAhGRNKdEICKS5pQIRETSnBKBiEiaUyKQRsvMLjGzj6NebzezvaspP9fMRjZAXHeY2bPxfp9K3vdcM3vHzHL3YBsjzGxh1OulZnZcbcua2a/M7LGGilviS4lAas3MLjCzKeGO+Rsze8PMjoj3+7p7c3dfEsbwlJn9psLyge7+frzjiGZmPc1saQO8zxDgMuBMdy+q63bc/SN332dPy7r73e5+eRhbTzNzM8uKV9wSX0oEUitm9gvgAeBuoAPQHfgLcEYCw0p57j7d3U909x2JjqU2kjXudKNEIDEzs1bA/wLXuPsEd9/h7iXu/pq7/zIs08TMHjCzVeHjATNrEi4baWYFZvbfZrY2PJq4NGr7bc3sVTPbamZfAL0rvL+bWR8zGwX8BPif8KjktXB5dNNFneOopN69zOwDM9tmZu8A+dWUvdnMvgrLzjOzsyosv8LM5kctPzCc383MJpjZOjPbYGYPR63zs3CdTWb2lpn1COebmf0prMMWM5tlZvuFy04Jt7/NzFaa2Y3Rda8Q9kFh2U1m9mR5E04VZctjim4e+zB83hx+H8Ori1saHyUCqY3hQC7wUjVlfg0cCgwGDgAOBm6NWt4RaAV0IWgyeMTM2oTLHgGKgE7Az8LH97j7GOAfwO/D5qLT6jmOisYCUwkSwF3AxVGxLHX3nlFlvwJGhNu+E3jWzDoBmNmPgDuAi4CWwOnABjPLBF4HlgE9w5jGheucCfwK+CHQDvgIeC58rxOAI4F+QGvgx8CGcNnjwJXu3gLYD/hPFXWDIKmeSJB4+1X4nGJxZPjcOvw+Pq0hbmls3F0PPWJ6EOwwVtdQ5ivglKjXJwJLw+mRwE4gK2r5WoIddiZQAvSPWnY38HHUawf6hNNPAb+p8N5LgeP2JI5K6tMdKAWaRc0bCzwb42c2AzgjnH4LuL6SMsOBddHxRC17A7gs6nUGUAj0AI4BFoWfX0aF9ZYDVwItK8wfCRRU+Myuinp9CvBVNWXLP987yj8DguTlFT7PKuNO9N+xHt9/6IhAamMDkF/ZScEonQl+2ZZbFs77dhvuXhr1uhBoTvCrMQtYUWHduqprHJVtZ5N/t427yrjM7CIzm2Fmm81sM8Gv8fKmpG4ECaqibsCyCvGU6wE8GLW9jYABXdz9P8DDBEdSa8xsjJm1DNc7m2Cnvixs1hpeVcx8/zPvXFXBWqgy7nrYttQzJQKpjU8Jmm7OrKbMKoKdQLnu4byarCP45d2twrpVqanb3LrGUdE3QBsza1ZTXGEb+N+A0UBbd28NzCHYAUKww+1dyaorgO5VJNgVBE08raMeee4+CcDd/+zuQ4GBBM06vwznT3b3M4D2wMvA+GrqWPEzr+3nVNl3UW3c0rgoEUjM3H0LcBtBe/qZZtbUzLLN7GQz+31Y7DngVjNrZ2b5Yfkar7l39zJgAnBHuN0BRLXFV2INUOU9BXWNo5K4lgFTgDvNLMeCy2QrOycB0Ixgp7gOIDwBvV/U8seAG81saHiit0+YPL4gSDj3mFkzM8s1s8PDdR4FbjGzgeE2W4XnGjCzg8zsEDPLBnYQJOmyMM6fmFkrdy8BtgJl1VTzGjPramZ7EbTrP1+7T4l1QITvfh9Vxi2NjxKB1Iq73w/8guCE4jqCX36jCX51AvyGYMc5C5gNTAvnxWI0QfPMaoJzAE9WU/ZxYEDY9PByJcv3JI6KLgAOIWjeuB34e2WF3H0ecB/BkdMaYH/gk6jlLwC/JTjHsI3gM9srTIKnAX0I2vYLCE784u4vAfcC48xsK8ERxsnhJlsSHIFsImjS2QD8MVx2IbA0XOcq4KfV1G8s8DawJHzU6nNy98KwXp+E38ehNcQtjYy5a2AaEZF0piMCEZE0p0QgIpLmlAhERNKcEoGISJqr7sagRik/P9979uyZ6DBERJLK1KlT17t7u8qWJV0i6NmzJ1OmTEl0GCIiScXMqrwjXk1DIiJpTolARCTNKRGIiKQ5JQIRkTSnRCAikubilgjM7IlwCL05VSw3M/uzmS0Oh9g7MF6xiIhI1eJ5RPAUcFI1y08G+oaPUcD/xTEWERGpQtzuI3D3D82sZzVFzgD+7kH3p5+ZWWsz6+Tu38QrJhERd2dnSRmbC0vYVFjMlsIStu0qpSzilEacskiEsgiURSKURpzIt/M9qkzwaOjem4f13Isj+1V6T9geSeQNZV347hB5BeG87yUCMxtFcNRA9+7VDVolIqmuuDTCzuIyCktKKSwuC6aLyygsDl5vLixh887i4Lmw/Dl6XgnFZZF6icWs5jL16aqjeqdcIqjsI6w0vbr7GGAMwLBhwzSAgkiK2VlcxuK121m4ZhuL1mxj8drtbN1ZEuzoS3bv5HcWl1EaiW0XkJOVQZum2bRpmkOrvGx65TcLpptm0zovhzZNs2ndNJtWeTm0yM0iOzODzAzIzMggK8PIjHpkZRgZ4XNmhpFpwbM1dCaIk0QmggK+O1ZqV+o2pqyIJImSsghfr9/BwtXBDr/8ednGQspbWXKyMujdrvm3O+q8nCyaZmeSl5NJ0/CRl5NFXnb59O75TXOyaB3u6PNyMhNb2SSSyETwKjDazMYRDAO4RecHRFJHWcT5ZPF6ZhVsZuGa7SxavY0l67dTUhbs8TMzjF75zRjYuRVnDenKPh2b069DC3q0bUZmRmr80k4WcUsEZvYcMBLIN7MCgrFeswHc/VFgInAKsBgoBC6NVywi0nBWbyni+ckreH7yclZtKQKg21557NOhBcfu2559OragX4cW7N2uGU2y9Ku9MYjnVUPn17DcgWvi9f4i0nDKIs6Hi9bxj8+X858Fa4g4jOibz62nDuDIfu1o3iTpOjpOK/p2RKTO1mwtYvzkFYybvIKVm3eS3zyHK4/qzfkHdad726aJDk9ipEQgIrVSFnE++nIdYz9fzr8XrKUs4hzRJ59fnbIvxw/oQE6Weq5JNkoEIhKTtVuLGD9lBc99Efz6b9ssh8tH9OL8g7rTM79ZosOTPaBEICLVWr2liHvemM/rs76hNOIc1rstt5zSnxMGdNSv/xShRCAilSopi/DkJ1/z4LtfUhpxLjmsJz85tAe99Os/5SgRiMj3TPpqPbe/Mpcv127n2P7tuf20gTr5m8KUCETkW2u2FvHbf83n1Zmr6Nomj8cuGsZxAzokOiyJMyUCEaGkLMLTk5bywLtfUlwW4bpj+/Lzkb3JzdYNX+lAiUAkzX2+ZAO3vTKXhWu2MXKfdtxx2kBdBZRmlAhE0tTabUX8buICXpq+ki6t8xhz4VCOH9AhZXrUlNgpEYikmdKyCM98toz7317ErtIIo4/uwzVH91FvnWlMiUAkjcxYsZmbX5zFgtXbOLJfO+48faAuBxUlApF08dwXy7ntlTnkN2/Coz89kBMHdlQzkABKBCIpr7g0wh2vzWXs58sZ0Tefh84fQuumOYkOSxoRJQKRFLZ2WxE/f3YaU5Zt4qqjevPLE/fRoC/yPUoEIilq+vJNXPXsVLbuLOWh84dw2gGdEx2SNFJKBCIpaPzkFdz68hzat2zCi1cfxoDOLRMdkjRiSgQiKaSkLMJdr8/j758u4/A+bXn4/ANp00znA6R6SgQiKWL99l38/NlpfLF0I1eM6MVNJ/UnK1PdREvNlAhEUsCsgs1c+cxUNu4o5sHzBnPG4C6JDkmSiBKBSJJ7cWoBt7w0m3bNg/MB+3VpleiQJMkoEYgkqZKyCHdPnM+Tnyxl+N5tefiCIbRt3iTRYUkSUiIQSUKbdhRz9T+m8tmSjfzs8F786hSdD5C6UyIQSTKRiDP6uWlMW76Z+350AGcP7ZrokCTJ6SeESJL520dL+GTxBu48faCSgNQLJQKRJDKrYDN/eGshJw3syHkHdUt0OJIilAhEksSOXaVcP24G7Vo04Z6z91fPoVJvdI5AJEnc+dpclm7YwXNXHKreQ6Ve6YhAJAm8PmsV46cUcM3IPhy6d9tEhyMpRolApJEr2FTILRNmM7hba64/rm+iw5EUpEQg0oiVlkX4r3EzcIc/nzeEbN0rIHEQ178qMzvJzBaa2WIzu7mS5a3M7DUzm2lmc83s0njGI5JsHnnvK6Ys28RdZw6ke9umiQ5HUlTcEoGZZQKPACcDA4DzzWxAhWLXAPPc/QBgJHCfmeksmAgwZelGHvz3Is4a0oWzhuh+AYmfeB4RHAwsdvcl7l4MjAPOqFDGgRYWXAfXHNgIlMYxJpGksLWohOvHzaBLmzz+94yBiQ5HUlw8E0EXYEXU64JwXrSHgX2BVcBs4Hp3j1TckJmNMrMpZjZl3bp18YpXpFFwd3790hxWby3iwfOG0CI3O9EhSYqLZyKo7G4Xr/D6RGAG0BkYDDxsZt8bU8/dx7j7MHcf1q5du/qOU6RRmTBtJa/NXMUNx/XlwO5tEh2OpIF4JoICIPoe+K4Ev/yjXQpM8MBi4GugfxxjEmnUlq7fwW2vzOGQXntx9cg+iQ5H0kQ8E8FkoK+Z9QpPAJ8HvFqhzHLgWAAz6wDsAyyJY0wijVZxaYTrxk0nKzODP/14MJkZ6kJCGkbcuphw91IzGw28BWQCT7j7XDO7Klz+KHAX8JSZzSZoSrrJ3dfHKyaRxuxP7y5iVsEW/u8nB9K5dV6iw5E0Ete+htx9IjCxwrxHo6ZXASfEMwaRZDBp8Xoe/eArzj+4Gyfv3ynR4Uia0W2KIgm2aUcxN4yfQa/8Zvy/UyveaiMSf0oEIgnk7tz04iw27Sjhz+cNoWmOOgSWhqdEIJIg7s7dE+fz9rw1/M9J+7Bfl1aJDknSlH5+iCRAaVmEm16czYvTCrhoeA9+dnivRIckaUyJQKSBFZWUMXrsNN6dv5YbjuvHdcf20WhjklBKBCINaMvOEq54egqTl23krjP348JDeyQ6JBElApGGsnZbERc/MZnFa7fx0PlDOHVQ50SHJAIoEYg0iGUbdnDh41+wfvsunrjkIEb0VZ9Z0ngoEYjE2bxVW7noiS8oi0QYe8WhDO7WOtEhiXyHEoFIHH2+ZAOXPz2F5rlZjBs1nD7tWyQ6JJHvUSIQiZN35q1h9NhpdG2TxzOXHaL+g6TRUiIQiYMXpqzg5gmz2a9LK5685CD2aqYRWKXxUiIQqWd//eArfvfGAkb0zefRnw6lWRP9m0njpr9QkXri7tzzxgL++uESTh3UifvPHUxOlnpxkcZPiUCkHkQizs0TZjF+SgEXHtqDO04fqIFlJGkoEYjUg3fnr2H8lAKuObo3N56wj7qMkKSi41aRevDYR1/TtU0eNxzXT0lAko4SgcgemrFiM18s3cjPDu9FVqb+pST56K9WZA/97aMltMjN4tyDuiU6FJE6USIQ2QMrNhbyxuxvuOCQ7jTXZaKSpJQIRPbAk58sJcOMSw/TwDKSvJQIROpoy84Snp+8nNMP6EzHVrmJDkekzpQIROrouS+Ws6O4jMtH7J3oUET2iBKBSB0Ul0Z46pOlHN6nLQM6t0x0OCJ7RIlApA7+NXsVq7cWcYWOBiQFKBGI1JK7M+bDr+nXoTlH9dNIY5L8lAhEamnSVxuY/81WLj9ib91FLClBiUCklv720RLymzfhjCEafF5SgxKBSC0sWrON9xeu4+LhPWiSlZnocETqRcyJwMyOMLNLw+l2ZqY7aCTtPPbREnKzM/jpoT0SHYpIvakyEZjZwKjp24GbgFvCWdnAszVt3MxOMrOFZrbYzG6uosxIM5thZnPN7IPahS/ScNZuK+Ll6av40dButNHQk5JCqjsi6GFm94TTZwGnAzsA3H0V0KK6DZtZJvAIcDIwADjfzAZUKNMa+AtwursPBH5UhzqINIhnPl1GSSTCZUfoYFhSS5W9ZLn7RDMrC18Wu7ubmQOYWbMYtn0wsNjdl4TrjAPOAOZFlbkAmODuy8P3XFuHOojEXWFxKc98tozj9+1Az/xY/vxFkke15wjc/a1wcryZ/RVobWZXAO8Cj9Ww7S7AiqjXBeG8aP2ANmb2vplNNbOLYg9dpOG8OLWAzYUljDpSN5BJ6omp31x3/6OZHQ9sBfYBbnP3d2pYrbILrL2S9x8KHAvkAZ+a2Wfuvug7GzIbBYwC6N69eywhi9Sbsojz+MdfM7hba4b2aJPocETqXUxXDZnZve7+jrv/0t1vdPd3zOzeGlYrAKJH6ugKrKqkzJvuvsPd1wMfAgdU3JC7j3H3Ye4+rF073ckpDevd+WtYuqGQK0boBjJJTbFePnp8JfNOrmGdyUBfM+tlZjnAecCrFcq8AowwsywzawocAsyPMSaRBvG3D5fQtU0eJw7skOhQROKi2qYhM7sa+Dmwt5nNilrUAvikunXdvdTMRgNvAZnAE+4+18yuCpc/6u7zzexNYBYQAR5z9zl1r45I/Zq+fBNTlm3i9tMGaDxiSVk1nSMYC7wB/A6Ivg9gm7tvrGnj7j4RmFhh3qMVXv8B+ENM0Yo0sMc++pqWuVmcO0zjEUvqquknjrv7UuAaYFvUAzPbK76hiSTWio2FvDHnGy44pAfNNB6xpLBYjghOBaYSXPETfabMAV1LJynr8Y+/JsOMSw7rmehQROKq2kTg7qeGz7qVUtLKlsISxk9ZwemDNR6xpL6aThYfWN1yd59Wv+GINA5jv1hOYXEZlx+hg15JfTU1Dd1XzTIHjqnHWEQaheLSCE9N+poj+uRrPGJJCzU1DR3dUIGINBavzVzFmq27uPfsQYkORaRB1NQ0dIy7/8fMfljZcnefEJ+wRBKjLOL87aMlGo9Y0kpNTUNHAf8BTqtkmQNKBJJS/vbREhas3sZD5w9RdxKSNmpqGro9fL60YcIRSZwFq7dy/9uLOGlgR04d1CnR4Yg0mFg7nWtrZn82s2lhd9EPmlnbeAcn0lCKSyP84vmZtMzL4rdn7aejAUkrsXaeMg5YB5wNnBNOPx+voEQa2kP/+ZJ532zl7rP2p23zJokOR6RBxXrf/F7uflfU69+Y2ZlxiEekwU1fvom/vP8V5wztygkDOyY6HJEGF+sRwXtmdp6ZZYSPc4F/xTMwkYaws7iM/x4/k44tc7nttAE1ryCSgmq6fHQbu/sY+gXwbLgoA9gO3B7X6ETi7N43F7Bk/Q7GXn4ILXOzEx2OSELUdNVQi4YKRKShTVq8nqcmLeWSw3pyWJ/8RIcjkjAx961rZm2AvsC3PXC5+4fxCEok3rYWlfDLf85i7/xm3HRS/0SHI5JQMSUCM7scuJ5g3OEZwKHAp6ivIUlSd702j2+27OSfVx9GXk5mosMRSahYTxZfDxwELAv7HxpCcAmpSNJ5d94aXphawNUje3Ng9zaJDkck4WJNBEXuXgRgZk3cfQGwT/zCEomPjTuKuXnCbPbt1JLrj+2X6HBEGoVYzxEUmFlr4GXgHTPbBKyKV1Ai8eDu/Pql2WzdWcKzlx9MTpYGoxeBGBOBu58VTt5hZu8BrYA34xaVSBy8MmMVb8xZzU0n9ad/R40zIFKuNlcNHQgcQXBfwSfuXhy3qETq2eotRdz2yhyG9mjDqCM16phItFg7nbsNeBpoC+QDT5rZrfEMTKS+uDv/8+IsSsqc+350AJkZ6lBOJFqsRwTnA0OiThjfA0wDfhOvwETqyz8+X86Hi9Zx1xkD6ZnfLNHhiDQ6sZ4tW0rUjWRAE+Creo9GpJ4t27CDuyfOZ0TffH56aI9EhyPSKNXU19BDBOcEdgFzzeyd8PXxwMfxD0+k7soizn+Pn0lmhvH7cwZpjAGRKtTUNDQlfJ4KvBQ1//24RCNSjx77aAlTlm3i/nMPoFOrvESHI9Jo1dTp3NPl02aWA5TfgbPQ3UviGZjInliybjv3vbOIEwd24KwhXRIdjkijFmtfQyMJrhpaStAldTczu1idzkljFIk4t0yYTW5WBnedqWEnRWoS61VD9wEnuPtCADPrBzwHDI1XYCJ19fyUFXz+9UbuPXt/2rfIrXkFkTQX61VD2eVJAMDdFwEaxUManbVbi7h74nyG792Wc4d1S3Q4Ikkh1kQw1cweN7OR4eNvBCeQq2VmJ5nZQjNbbGY3V1PuIDMrM7NzYg1cpDK3vzqX4tIId/9wfzUJicQo1kRwFTAXuI6gS+p54bwqmVkm8AhwMjAAON/MvjcobFjuXuCt2MMW+b4356zmjTmruf64vvTSjWMiMavxHIGZZQBT3X0/4P5abPtgYLG7Lwm3Mw44gyCJRLsWeJFgvAOROtmys4TbXpnDvp1acsUI9SUkUhs1HhG4ewSYaWbda7ntLsCKqNcF4bxvmVkX4Czg0eo2ZGajzGyKmU1Zt07j4cj33fvmAtZv38W9Z+9Pdqa6lxapjVivGupEcGfxF8CO8pnufno161TWQOsVXj8A3OTuZdW157r7GGAMwLBhwypuQ9Lc50s2MPbz5VwxoheDurZOdDgiSSfWRHBnHbZdAERfttGV7w9mMwwYFyaBfOAUMyt195fr8H6ShopKyrhlwmy67ZXHDcdrxDGRuqipr6FcgpPCfYDZwOPuXhrjticDfc2sF7ASOA+4ILqAu/eKeq+ngNeVBKQ2Hv7PYpas38Ezlx1M05yYh9cQkSg1/ec8DZQAH7H76p/rY9mwu5ea2WiCq4EygSfcfa6ZXRUur/a8gEhNFqzeyqMffMXZB3ZlRN92iQ5HJGnVlAgGuPv+AGb2OPBFbTbu7hOBiRXmVZoA3P2S2mxb0ltZxLnpxdm0ysvm1h/sm+hwRJJaTZdXfNuxXC2ahETi7ulJS5m5YjO3nz6QNs1yEh2OSFKr6YjgADPbGk4bkBe+NsDdXSOAS4NbsbGQP769kGP6t+e0QZ0SHY5I0qupG+rMhgpEJBbuzq9fnoOBehYVqSe680aSyiszVvHhonX88sR96NJag82I1AclAkkaG7bv4s7X5jKke2suHN4z0eGIpAwlAkkav/nXfLbvKuXesweRmaEmIZH6okQgSeGDRet4afpKrh7Zh34dWiQ6HJGUokQgjd6OXaX8asJserdrxjVH9050OCIpR/fkS6P3h7cWsnLzTv551XCaZOlCNpH6piMCadTeW7iWpyYt5ZLDejKs516JDkckJSkRSKO1btsufvnCTPp3bMHNJ/dPdDgiKUtNQ9IoRSLOjS/MZFtRKWOvOJTcbDUJicSLjgikUXpy0lI+WLSOW08doKuEROJMiUAanbmrtnDvGws4bt8O/PSQ2o6QKiK1pUQgjcrO4jKue246bZpl8/tzBqkvIZEGoHME0qjc9a95LFm/g2cvO4S91L20SIPQEYE0Gm/O+Yaxny/nyiN7c3if/ESHI5I2lAikUfhmy05uenE2g7q24hcahF6kQSkRSMKVRZwbnp9BSVmEB88bQk6W/ixFGpLOEUjCPfrBV3y2ZCN/OGcQvfKbJTockbSjn16SUNOXb+L+dxZx6qBOnDO0a6LDEUlLSgSSMNt3lXL9uBl0bJnLb8/aX5eKiiSImoYkYW57eQ4Fmwp5/srhtMrLTnQ4ImlLRwSSEC9PX8mE6Su59pi+HKReRUUSSolAGtzyDYXc+vIchvVow7XH9El0OCJpT4lAGlRpWYTrn5+OGTxw3mCyMvUnKJJoOkcgDerP//6S6cs38+fzh9C1TdNEhyMi6IhAGtBnSzbw8HuLOWdoV04/oHOiwxGRkBKBNIhVm3cyeuw0erZtxh2nD0x0OCISRU1DEndFJWVc/exUikoijBs1lOZN9Gcn0pjoP1Liyt359UtzmFmwhTEXDqVPe402JtLYxLVpyMxOMrOFZrbYzG6uZPlPzGxW+JhkZgfEMx5peE9PWsqL0wq47ti+nDCwY6LDEZFKxC0RmFkm8AhwMjAAON/MBlQo9jVwlLsPAu4CxsQrHml4ny3ZwF3/ms9x+7bnv47tm+hwRKQK8TwiOBhY7O5L3L0YGAecEV3A3Se5+6bw5WeAeh1LESs37+Saf0yjR9um/OnHg8nIUD9CIo1VPBNBF2BF1OuCcF5VLgPeqGyBmY0ysylmNmXdunX1GKLEQ1FJGVc+M4VdpRHGXDiMFrnqR0ikMYtnIqjsJ6BXWtDsaIJEcFNly919jLsPc/dh7dq1q8cQpb65O7+aMJs5K7fywI8H06d980SHJCI1iOdVQwVAt6jXXYFVFQuZ2SDgMeBkd98Qx3ikATz5yVImTF/JDcf147gBHRIdjojEIJ5HBJOBvmbWy8xygPOAV6MLmFl3YAJwobsvimMs0gAmfbWe306czwkDOqgzOZEkErcjAncvNbPRwFtAJvCEu881s6vC5Y8CtwFtgb+Eg5KUuvuweMUk8VOwqZDRY6fTK78Z9517gE4OiySRuN5Q5u4TgYkV5j0aNX05cHk8Y5D421lcxpXPTKWkLMKYC4fq5LBIktGdxbJH3J1bJsxi3jdbefziYezdTieHRZKNOp2TPfL4x1/z8oxV/Pfx/Timv04OiyQjJQKps08Wr+fuifM5aWBHrjlaJ4dFkpUSgdTJio2FjB47jd7tmvPHcw8gPNkvIklIiUBqbceuUq58ZiqlEWfMRcPUrbRIktN/sNTK8g2FjHpmCovWbOPxSw6iV36zRIckIntIiUBi9uGidVz73HQAnv7ZwYzoq+4+RFKBEoHUyN0Z8+ES7n1zAf06tGDMhcPo3lYDz4ukCiUCqVZhcSk3vTib12au4gf7d+IPPxpE0xz92YikEv1HS5VWbCxk1DNTWbB6Kzed1J+rjtpbVweJpCAlAqnUJ4vXc83YaUQizpOXHMTIfdonOiQRiRMlAvkOd+fxj7/m7onz6dO+OWMuHEZPXRkkktKUCORbO4vLuHnCLF6ZsYqTBnbkj+ceoHsERNKA/ssFCLqRHvX3qcxfvZUbT+jHz0f2UVfSImlCiUCY9NV6Ro+dTklphMcvHqbO40TSjBJBGispi/D0pKX87o0F9GzblDEXDaO3upEWSTtKBGlo8drtjJ+yggnTCli/vZjjB3Tg/nMP0IAyImlKiSBN7NhVyr9mf8Pzk1cwddkmMjOMY/q358fDunFM//Y6HyCSxpQIUpi7M33FZsZPXsFrM1exo7iMvfObcfPJ/fnhgV1o3yI30SGKSCOgRJCCNmzfxUvTV/L85BV8uXY7edmZ/GBQJ358UDeG9Wiju4NF5DuUCFJEWcT58Mt1jJ+8gnfnr6GkzBncrTW/++H+nDqok9r/RaRKSgRJLhJxXpm5kvvfWcSKjTtp0zSbi4b35Nxh3dinY4tEhyciSUCJIIl99OU67nljAXNXbWVg55Y8csG+HDegPU2yMhMdmogkESWCJDR31RbueWMBH325nq5t8njwvMGcNqizrvwRkTpRIkgiBZsKue/tRbw8YyWt8rK59Qf7cuHwHjoCEJE9okSQBDYXFvPIe4t5etIyzODKI3tz9cjetMrTCWAR2XNKBI1YUUkZT09ayiPvLWbbrlLOObArNxzfj86t8xIdmoikECWCRqgs4rw8fSX3vb2QVVuKOHqfdtx0cn/6d2yZ6NBEJAUpESSAu1NUEmHbrhK2FZWyraiU7UWlbCsqYWNhMc98uowFq7cxqGsr/njuARzWOz/RIYtIClMiqAeRiLNhRzGrNu/kmy07WbW5iFWbd7JhRzHbinbv7LftKgl3+KWURrzK7XXbK4+Hzh/CD/bvpCuBRCTu4poIzOwk4EEgE3jM3e+psNzC5acAhcAl7j4tnjHVlruzbVdpsJPfXMSqLTsrTBexeksRxWWR76yXm51BuxZNaNEkmxa5WXRunUuL3BY0b5JFi9wsWuRm0zw3i5a5wevmYbnmTbLo1CqXrMyMBNVYRNJN3BKBmWUCjwDHAwXAZDN71d3nRRU7GegbPg4B/i98jqvi0gibCotZv30XG7YXs2FH8Lx+e/m8XWzYURzO28Wu0u/u5DMzjI4tc+ncOpfB3VrTaf9curTOo1OrPDq1CqZbN81Wnz4ikhTieURwMLDY3ZcAmNk44AwgOhGcAfzd3R34zMxam1knd/+mvoN5b+Fa7np9Hhu2F7NlZ0mlZXIyM2jbPCd4NGtCn/bNyW/ehPzmOXQOd/SdW+fSvkUumWqyEZEUEc9E0AVYEfW6gO//2q+sTBeg3hNB67xs9u3UkvxmObRt3uTbnX1+892vWzTJ0q94EUk78UwEle1RK54hjaUMZjYKGAXQvXv3OgUzpHsbHrmgTZ3WFRFJZfE8I1kAdIt63RVYVYcyuPsYdx/m7sPatWtX74GKiKSzeCaCyUBfM+tlZjnAecCrFcq8ClxkgUOBLfE4PyAiIlWLW9OQu5ea2WjgLYLLR59w97lmdlW4/FFgIsGlo4sJLh+9NF7xiIhI5eJ6H4G7TyTY2UfPezRq2oFr4hmDiIhUT3ctiYikOSUCEZE0p0QgIpLmlAhERNKcBedrk4eZrQOWJTqOGOUD6xMdRJykct0gteunuiWvPalfD3ev9EaspEsEycTMprj7sETHEQ+pXDdI7fqpbskrXvVT05CISJpTIhARSXNKBPE1JtEBxFEq1w1Su36qW/KKS/10jkBEJM3piEBEJM0pEYiIpDklgj1gZkvNbLaZzTCzKeG8vczsHTP7MnxuE1X+FjNbbGYLzezExEUemyrq9wczW2Bms8zsJTNrHVU+aepXWd2ilt1oZm5m+VHzkr5uZnZtGP9cM/t91PykqRtU+Xc52Mw+K59nZgdHlU+a+oXD9f4z/B+bb2bDG2Sf4u561PEBLAXyK8z7PXBzOH0zcG84PQCYCTQBegFfAZmJrkMd6ncCkBVO35us9ausbuH8bgRdpy8rX54KdQOOBt4FmoSv2ydj3aqp39vAyeH0KcD7yVg/4Gng8nA6B2jdEPsUHRHUvzMIvkzC5zOj5o9z913u/jXBGAwHf3/1xs3d33b30vDlZwSjykGK1A/4E/A/fHfI1FSo29XAPe6+C8Dd14bzU6FuEHxfLcPpVuwe6TBp6mdmLYEjgccB3L3Y3TfTAPsUJYI948DbZjY1HFcZoIOHo6yFz+3D+V2AFVHrFoTzGrPK6hftZ8Ab4XSy1e97dTOz04GV7j6zQtmkrxvQDxhhZp+b2QdmdlA4P9nqBpXX77+AP5jZCuCPwC3h/GSq397AOuBJM5tuZo+ZWTMaYJ8S14Fp0sDh7r7KzNoD75jZgmrKWiXzGvu1u9+rn7t/CGBmvwZKgX+EZZOtfpV9d78maPqqKBXqlgW0AQ4FDgLGm9neJF/doPL6nQPc4O4vmtm5BL+qjyO56pcFHAhc6+6fm9mDBE1BVam3uumIYA+4+6rweS3wEsFh2Roz6wQQPpcfghcQtD+X68ruw9dGqYr6YWYXA6cCP/GwsZIkq18ldTuKoJ11ppktJYh/mpl1JPnrdjBBHSZ44AsgQtCBWVLVDaqs38XAhLDIC+xuIkmm+hUABe7+efj6nwSJIe77FCWCOjKzZmbWonya4JfkHOBVgj9KwudXwulXgfPMrImZ9QL6Al80bNSxq6p+ZnYScBNwursXRq2SNPWrom6T3b29u/d0954E/2QHuvtqkr9uc4CXgWPC+f0ITkSuJ4nqBtXWbxVBMoegnl+G00lTv/BvbYWZ7RPOOhaYRwPsU9Q0VHcdgJfMDILPcay7v2lmkwkOuy8DlgM/AnD3uWY2nuCLLQWucfeyxIQek6rqt5jgKoV3wmWfuftVSVa/SutWVeFUqJuZ5QBPmNkcoBi4ODyaS6a6QdX12w48aGZZQBEwCpLuuwO4FvhH+H0tAS4l+MEe132KupgQEUlzahoSEUlzSgQiImlOiUBEJM0pEYiIpDklAhGRNKdEkOIs6EXzvqjXN5rZHfW07e31sZ362HbYI2V+OD0pPlHVHEP0cw1l3zeznjGWvcTMHq5FHCPN7PVYy++p8P0Oq2b56WZW3R2ysbxHZtilxJFR8942sx/tyXYloESQ+nYBP7SoLpUbUnhdd4Ny9yp3ShIXI4FKP3Mzy3L3V939nj15g/D6+J8Dj5hZtpmdH8z2F/ZkuxJQIkh9pQTjnN5QcYGZ9TCzf1swtsC/zax7OP8pM/s/M3vPzJaY2VFm9oQF/aM/VWEb95nZtHD9duG8983sbjP7ALjezIZa0NHZVDN7q/x2+Qrb6WVmn5rZZDO7q8KyX4bzZ5nZnTVVuPxowswyzOwvFvS//7qZTTSzc8JlS83szjD22WbWP5zfLKzrZAs6/jojnD/QzL6woL/7WWbWt5K3Xhf9HP5Sft929y//DwvvhAI2AmVR61Ssw6Vmtij8DA+Pmt/OzF4M45tsZodXtn5U+YPNbFJYl0m2+67V6DIjw+9nfPie95jZT8L6zjaz3mG50yzotG66mb1rZh3MrCdwFXBD+NmMCP9+7jez94B7o49ownVeMrOZ4eOwcP5Poz7fv5pZZsU4w64XJgF3AHcD11RXd6mFRPS5rUfDPYDtBN3zLiXonvdG4I5w2WsEd5hC0JPoy+H0U8A4gk6tzgC2AvsT/HCYCgwOyzlBf0MAtwEPh9PvA38Jp7MJ/nnbha9/DDxRSZyvAheF09cA28PpEwgSmYXv/zpwZCXrL2X3+AHl654DTAzX6whsAs6JKn9tOP1z4LFw+m7gp+F0a2AR0Ax4KKquOUBeDJ/9SGALQR8wGcCnwBExrNeJ4A7SduF7fRL12Y4t3wbQHZhfxfu+Hk63ZPf4EccBL1ZRfnP4vk2AlcCd4bLrgQfC6Tbsvgn1cuC+cPoO4Mao7T0Vfk+Z4etLouJ/HvivcDqT4G9yX4K/xexw/l/K/xYqiXUvYAfw20T/b6XSQ11MpAF332pmfweuA3ZGLRoO/DCcfoZgAIxyr7m7m9lsYI27zwYws7lAT2AGQcdlz4fln2V3p19Ezd8H2I/dXVJkAt9UEubhwNlRsdwbTp8QPqaHr5sT9KnyYQ3VBjgCeMHdI8Dq8BdqtPJ4p7L7czgBON3Mbgxf5xLscD8Ffm1mXQk6b/uS2Hzh7gUAZjaD4LP7uIZ1DiEYWKX8yOJ5gm6kIdiZD9h9YEFLM2vh7tuq2FYr4OnwCMYJEnNlJnvY1bGZfUUw0AvAbIJBbSBIaM+HR3Q5wNfV1OEFr7y7g2OAi+Db5p4tZnYhMBSYHNYrj90dq1V0JEFy3a+a95ZaUiJIHw8A04AnqykT3d/IrvA5EjVd/rqqv5vo9XeEzwbMdffhMcRYWX8nBvzO3f8aw/qVrVud8nqVsbtOBpzt7gsrlJ1vZp8DPwDeMrPL3f0/McQQ/dlFv09Nqur7JQMY7u47q1he0V3Ae+5+VtiM834V5Sp+x9Hff3nMDwH3u/urZjaS4EigKjuqWVaRAU+7+y3VFgo6mfs9QTJ5wsxOcfeJtXgfqYLOEaQJd98IjAcui5o9CTgvnP4JNf9SrSiDoPkF4IIq1l8ItDOz4QAWnOgbWEm5TyrEUu4t4Gdm1jxcv4sF/dDH4mPg7PBcQQeCJpCavAVcW96Wb2ZDwue9gSXu/meCZqxBMcZQF58DI82srZllE3YyFnobGF3+wswG17CtVgRNPRA00eyJ6G1dHDV/G9Aixm38m2C0tPIrgVqG884p/14tGKO3RyXr3gaMd/cFBM15fzKz3NpXQypSIkgv9xH0QV/uOuBSM5sFXEjQHlwbO4CBZjaV4Ffa/1Ys4O7FBMniXjObSdCkVNkVJtcD11jQe2urqPXfJmgX/zRspvonse90XiToTnoO8FeCHeyWGta5i6D5ZJYFPXWWn7j+MUE33DOA/sDfY4yh1sImmjsImqPeJTiSK3cdMCw8YT2P4ERtdX4P/M7MPiFoltsTdwAvmNlHBF1Yl3sNOKv8ZHEN27geODr8LqcCA919HnArwahjs4B3CM5XfMvMBgBnAb8FcPcZBEn7pj2sk6DeRyXFmVlzd99uZm0J+mo/3IN+30UkpHMEkupeN7PWBCc371ISEPk+HRGIiKQ5nSMQEUlzSgQiImlOiUBEJM0pEYiIpDklAhGRNPf/AdBwVfnAXk25AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "p=784\n",
    "NB = 1000\n",
    "N_acc = np.linspace(490,600,23,dtype=int)\n",
    "succes_acc=np.zeros(len(N_acc))\n",
    "\n",
    "for i in range(len(N_acc)):\n",
    "    n = N_acc[i]\n",
    "    for j in range(NB):\n",
    "        X=np.random.randn(n,p)    \n",
    "        Xtilde=X @ I\n",
    "        Sol_SLOPE = DR_SLOPE(X,Xtilde,Lambda) \n",
    "        Arrondi_Sol=np.round(Sol_SLOPE, decimals=0)\n",
    "        if np.all(Arrondi_Sol == I):\n",
    "            succes_acc[i]=succes_acc[i]+1/NB   \n",
    "\n",
    "plt.ylabel('Probabilité')  # Titre de l'axe y\n",
    "plt.xlabel('Nombre de lignes \"n\" de la matrice X')  # Titre de l'axe x\n",
    "plt.plot(N_acc, succes_acc)\n",
    "plt.title('Condition d\\'accéssibilité') # Titre\n",
    "plt.show()            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**l'algorithme FISTA :** Cet algorithme permet de calculer l'estimateur SLOPE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def FISTA(A, y, Lambda, itmax=1e3, rtol=1e-8):\n",
    "    def f(x):\n",
    "        return 0.5*np.linalg.norm(A@x-y, ord = 2)**2\n",
    "    \n",
    "    def grad_f(x):\n",
    "        return A.T@(A@x-y)\n",
    "        \n",
    "    p,n = np.shape(A)\n",
    "    gamma = 1./(np.linalg.norm(A, ord = 2)**2) # calcul de gamma = 1/constante de Lipschitz\n",
    "    x = np.zeros(n) # vecteur initial x0 \n",
    "    x_tilde = x # initialisation de x_tilde\n",
    "    x_old = x # initialisation de x_old = x_k-1\n",
    "    t = 1. # init de t\n",
    "    t_old = t # init de t_old = t_k-1\n",
    "    test = False # initialisation de la condition d'arrêt (convergence)\n",
    "    F = [f(x)+norme(x,Lambda)] # initialisation de la liste des F(xk)\n",
    "    \n",
    "    it = 1\n",
    "    \n",
    "    while test == False  and it < itmax: \n",
    "        zk = x_tilde - gamma*grad_f(x_tilde) \n",
    "        x = prox_L1_ordonnee_rapide_depistage(zk,gamma*Lambda)\n",
    "        t = 0.5*(1.+ np.sqrt(4.* t**2 + 1.))\n",
    "        x_tilde = x + ((t_old -1.)/t)*(x-x_old)        \n",
    "        x_old = x # on fixe xk et tk pour la prochaine iteration\n",
    "        t_old = t\n",
    "        F.append(f(x)+norme(x,Lambda)) # on calcule F(xk)\n",
    "        test = np.abs(F[-1]-F[-2]) < np.abs(F[-1]*rtol) # on vérifie la condition de convergence\n",
    "        it+=1\n",
    "    return x,F"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Application de l'opérateur proximal :** Calcul de l'estimateur SLOPE pour un paramètre de régularisation donné et application de l'opérateur proximal de la norme $\\ell_1$ ordonnée à cet estimateur. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAU5klEQVR4nO3dXWxV15UH8P/CEMAmgM01xgYTQvNJBgUqi4ySEcqkmirkhfQhVfNQZaRo6EMitVIfJso8NI/RaNqqD6MqdIJKR51UVWiUPCSTIlQpqpCamIgQAhpCkEOMAZtP822w1zz4IDnEZ/1v7rk+9yr7/5OQzV0+526fe5evfddee5u7Q0S++WY1egAiUg4lu0gilOwiiVCyiyRCyS6SiNll3lmlUvFVq1blxsfHx8PjZ83K/9lkZuGxExMTheKzZ+dfKnYsGxuLX79+PYzPmTMnN3bjxo3wWHbNb7vttjDOxl4EG1tLS0sYjx6X6LlUDXZd2fmL3n+egYEBnDp1atoHpVCym9njAH4FoAXAf7n7y9HXr1q1Cu+//35ufHR0NLy/tra23FiUjABw9erVMH758uUw3t7enhu7dOlSeOy8efPCeJSsAHDixIkwvmzZstzYmTNnwmNPnz4dxleuXBnG586dG8aLOHfuXBhfvHhxGL9y5UpujP0QYz9I2HWdP39+GI+eE+wHaPRD8KGHHsqN1fzjxcxaAPwngE0A1gB42szW1Ho+EZlZRX6X2ADgsLsfcfcxAH8AsLk+wxKReiuS7MsBfDHl/4PZbV9iZlvMrN/M+kdGRgrcnYgUUSTZp/vD4itzb919q7v3uXtfZ2dngbsTkSKKJPsggN4p/18BYKjYcERkphRJ9g8A3G1md5rZbQB+AOCt+gxLROqt5tKbu98ws+cBvIvJ0ts2d/+kyGBYGScqxVQqlULnZqWYqGZ7++23h8eyzkJWT2bnj7B6bm9vbxhnxx8+fDiM33XXXWE80traGsavXbsWxqPyFyvFssckKgMD/DGP5k6wY6PnYhQrVGd397cBvF3kHCJSDk2XFUmEkl0kEUp2kUQo2UUSoWQXSYSSXSQRpfazT0xM0NpoJKqNnjp1KjyW1dFZbTNqp2T95kX6rquJR73VrB7MHo+hoXhSZJE6+sWLF8M4a1tm1yUae09PT3gsw/rZixzPWqIj0bwIvbKLJELJLpIIJbtIIpTsIolQsoskQskukohSS2+zZs0KW01ZO+WKFStqvm+2Aiwr80SlufPnz4fHLlmyJIwXXQb77NmzuTFWehseHg7jY2NjYbzIUtLd3d1hnC1jxlp/9+zZkxtjLa4HDhwI42vWxGursvJZkZVvo+dD9DzVK7tIIpTsIolQsoskQskukgglu0gilOwiiVCyiySi9BbXaLdUtotrtLQw2zVzwYIFYZwtHRy1Y7JdWFkrJ6sXs7prVNON6rkAsHbt2jDOdrfdsGFDGI927WXzE1gbaTS/AABWr14dxiOs5ZnFmei5zOaERC3TqrOLiJJdJBVKdpFEKNlFEqFkF0mEkl0kEUp2kUSU3s8e1bvZkstRbzXrhY+2ewZ4LTsaG6vhs6WmWb2ZjS363tvb28Nj2ZbNrC87qqMD8bW5cOFCeOxTTz0Vxl9//fUwXgSro7M5Iewxi+aFsPULonkd0foChZLdzAYAXAAwDuCGu/cVOZ+IzJx6vLL/o7vHOzSISMPpb3aRRBRNdgfwZzPbY2ZbpvsCM9tiZv1m1s/WFBORmVM02R9x928D2ATgOTPbeOsXuPtWd+9z977Ozs6CdycitSqU7O4+lH0cBvAGgLgFSkQapuZkN7M2M7v95ucAvgtgf70GJiL1VeTd+C4Ab2R1vdkA/sfd/zc6YGJiIqytsjr7okWLcmOsH53V4VnfdqVSyY2xnnH2fbE4G3sUv+eee8JjDx06FMaLirZNZtftlVdeCeOszv7ggw/mxj766KPwWDY3gu0zwBw9ejQ3xvZHiB7vGamzu/sRAPlXU0SaikpvIolQsoskQskukgglu0gilOwiiSi9xTVaNjkqRwBxuyYrTy1cuDAeHBFt4cvaSNlS0WxJZHb+6Htbt25deCxrz2XLYLOSZVTCYuVS1l5bZDlntlwza1Fl3/fixYvDeFRGZiXJaLvpaPltvbKLJELJLpIIJbtIIpTsIolQsoskQskukgglu0giSq2zM1HtEYhrm11dXfUezpfceeeduTG2XfTw8HAYZ8s5R22LQFzHZy2srI5+6lS8lij73qNlj4u2iTJRHX/u3Lk1HwvwWjirs0fXhc0ZifIkapfWK7tIIpTsIolQsoskQskukgglu0gilOwiiVCyiySi1Dr7xMREWNdl9cWodhlt5wwAZ86cCeNLliwJ44cPH86NrV27Njx26dKlYZzV0VlP+okTJ3Jj7LoU6QkHgGvXroXxqJ7NesLZls6tra1hPNr6mNXZWZzV0aOecyCus7Maf/R4R+sH6JVdJBFKdpFEKNlFEqFkF0mEkl0kEUp2kUQo2UUSUfq68dE65ay3OsL6rtm68ayme999933tMVUr2loYAPbu3VvzuVnP+Pnz58M42066ra3ta4/pJlYnZ/FojXQgnrfB1o2P6uAAn1/A6vDRvA+2ln93d3duLFrvnr6ym9k2Mxs2s/1Tbusws51m9mn2Md7FQEQarppf438L4PFbbnsBwC53vxvAruz/ItLEaLK7+3sAbv2dYzOA7dnn2wE8Wd9hiUi91foGXZe7HweA7GPu5G8z22Jm/WbWPzIyUuPdiUhRM/5uvLtvdfc+d+/r7Oyc6bsTkRy1JvtJM+sGgOxjvHyqiDRcrcn+FoBnss+fAfBmfYYjIjOF1tnN7DUAjwKomNkggJ8BeBnAH83sWQBHATxVzZ2Nj4+HdV1WXwzXxCa98OzcTK17YgN8PfyZxPqq2fwC1u9++vTpMB6tac/OzfZnL7LHOnv/iO1DEPXKA3yvgGiNA7b2QnRdonHRZHf3p3NC32HHikjz0HRZkUQo2UUSoWQXSYSSXSQRSnaRRJTa4jo+Ph6WeliJ6uzZs7mxZcuWhcey9llWmotaHll7bbS8L8BbOR977LEw/u677+bGWIsqa8VkJaaotAbEpb2BgYHw2HvvvTeMs9JctBw0WzqctbB2dHSEcfaciEqirCwYXZeoBK1XdpFEKNlFEqFkF0mEkl0kEUp2kUQo2UUSoWQXSUSpdfaWlpawLsvaMaN69JUrV2oeF8DbVD///PPcWE9PT6H7Zlsy79y5M4xHraKsZsvmNrB687lz58J4VMdny3MfOHAgjLM6fLQlNGuJZls2M5VKJYyPjo7mxnp7e8Nj2ZbOefTKLpIIJbtIIpTsIolQsoskQskukgglu0gilOwiiSh9y+aofsn6k2utLwJ82WLWfxzVPqMli6uxe/fuML58+fKaz83GxpYtZusEsC2ho356dizbypr12kdrGBSts7P5BWx9hGi+iZmFx0ZrBETXRK/sIolQsoskQskukgglu0gilOwiiVCyiyRCyS6SiKZaN57VF6Pa55EjR8JjV65cGcbZOuJRTXhoaCg8tr29PYxH6+EDwLFjx8J4VFtl/eqs3sz62aMtuIF4fXW2nn60Vj8ADA4OhnH2mEfY+ghtbW1hnK3NED2X2fyBaO5ElEP0ld3MtpnZsJntn3LbS2Z2zMz2Zv+eYOcRkcaq5tf43wJ4fJrbf+nu67J/b9d3WCJSbzTZ3f09APGcShFpekXeoHvezPZlv+bn/lFqZlvMrN/M+k+fPl3g7kSkiFqT/dcAvgVgHYDjAH6e94XuvtXd+9y9j70JJiIzp6Zkd/eT7j7u7hMAfgNgQ32HJSL1VlOym1n3lP9+D8D+vK8VkeZA6+xm9hqARwFUzGwQwM8APGpm6wA4gAEAP6r2DqP+ZlaznT9/fm5s9erV4bGs7sl65aOaL+s3v//++8M4q8P39fWF8ai2ymrZrN/9xIkTYZx9b1G9mu1xztY3GB4eDuPR9z42NhYey/YRYP3ubI5ANL+BzTdhY8tDk93dn57m5ldrujcRaRhNlxVJhJJdJBFKdpFEKNlFEqFkF0lEqS2us2fPDsstbDnnaAte1pLISkxsWeMiDh48WOj4tWvXhvGoVMPaJVl5ayaxkuNnn31W6PxFSpJsKWhWJo62qgbirbRZGTgq20XH6pVdJBFKdpFEKNlFEqFkF0mEkl0kEUp2kUQo2UUSUWqdHYjrvgsXLgyPjWrlbMljtmQyi0fYUtCsnswcOnQojEfb//b09BS6766urjB+8uTJms/N6uhs7kOtrZ4Ab0FlLbCsjs7mN0QtstHjyUTXTK/sIolQsoskQskukgglu0gilOwiiVCyiyRCyS6SiNLr7FGPMauVnzmTv+Ucqwezcxepu7Lte1tbW8P40qVLwzhbzvnixYu5saJbbrE6+qZNm8L4O++8U/N9szr65s2bw/ibb76ZG2O7E7HlnBk2byPqp2f3HT0X3T1/TOFZReQbQ8kukgglu0gilOwiiVCyiyRCyS6SCCW7SCJKrbNPTEzg0qVLuXHWvxzVRtla2+fOnQvjbJ3wqJa9bNmy8Fi2ve/AwEAYL4L1VTOsr5vp7OzMjbG119lW2FEdnWH7DLDnE3uuXrhwIYxHPeujo6PhsdG8jqhGT1/ZzazXzP5iZgfN7BMz+3F2e4eZ7TSzT7OPxVZoEJEZVc2v8TcA/NTd7wfw9wCeM7M1AF4AsMvd7wawK/u/iDQpmuzuftzdP8w+vwDgIIDlADYD2J592XYAT87QGEWkDr7WG3RmtgrAegB/A9Dl7seByR8IAKad4G1mW8ys38z62V5uIjJzqk52M1sAYAeAn7h7/A7CFO6+1d373L2vUqnUMkYRqYOqkt3M5mAy0X/v7n/Kbj5pZt1ZvBvA8MwMUUTqgZbebPK9/FcBHHT3X0wJvQXgGQAvZx9pHWTWrFlh2WB4OP55EZVxWKmElcdYiSpqgWUtqIODg2GctcgWwZa5Zq2YbKvrO+64I4yz8lokKndWY+PGjbmxHTt2hMey0lpUQgaKLQcdtakCtS97Xk2d/REAPwTwsZntzW57EZNJ/kczexbAUQBP1TQCESkFTXZ3/yuAvEr9d+o7HBGZKZouK5IIJbtIIpTsIolQsoskQskukohSW1yvX78e1tKLbLvMao+spZHVk6M2VdYey7Dvm7XIdnR05Mai5beB4i2wR48eDeMPP/xwbmz37t3hsaxGz2rdUbvnvn37wmMfeOCBMM7mL7DHLGpjZUtJs2XP8+iVXSQRSnaRRCjZRRKhZBdJhJJdJBFKdpFEKNlFElFqnX3OnDl0a+VItM0t6z9mdc9jx46F8ag/mZ17/vz5YZzNAWDLEs+bNy83xq4Lw7ZNbmlpCeNRHT96PKu5b3bdI+vXrw/jbN5Gb29vGGePWaS9PV6oOerzj663XtlFEqFkF0mEkl0kEUp2kUQo2UUSoWQXSYSSXSQRTbVlM1s/PerzZT3AIyMjYZzV/6N6Mqs1M2zsly9fDuNRnZ716S9atCiMs+2D2RyCKM761VkdndWyo++NrUnP1ihgcyPYOgLRevtsD4TW1tbcWLjmQ3hWEfnGULKLJELJLpIIJbtIIpTsIolQsoskQskukohq9mfvBfA7AMsATADY6u6/MrOXAPwLgJsF7Bfd/e3oXGx/9iJru7N6MKvhs5rtkiVLcmNXr14Nj436zesRj9adZ7Vqtrc829d+bGwsjEdjW7x4cXhstMcAACxdujSMR/MT2PwChj0mK1asCOPRc4adu9a1/quZVHMDwE/d/UMzux3AHjPbmcV+6e7/UdM9i0ipqtmf/TiA49nnF8zsIIDlMz0wEamvr/U3u5mtArAewN+ym543s31mts3Mpl1Lx8y2mFm/mfWzKasiMnOqTnYzWwBgB4CfuPsogF8D+BaAdZh85f/5dMe5+1Z373P3vs7OzuIjFpGaVJXsZjYHk4n+e3f/EwC4+0l3H3f3CQC/AbBh5oYpIkXRZLfJlqxXARx0919Mub17ypd9D8D++g9PROqlmnfjHwHwQwAfm9ne7LYXATxtZusAOIABAD9iJ7p+/TqGhoZy46zMw84dWbhwYRhn2yZHZRzW4srKeqw8VmQ7aVYaY6297h7GWato1I7JxlapVMI4O56VWyOsDMy2TWZLeLNlsiNRC2z0eFXzbvxfAUzXcB3W1EWkuWgGnUgilOwiiVCyiyRCyS6SCCW7SCKU7CKJKH3L5p6entz42bNnaz531IIK8Lopqxd3dHTkxtjSvwyrozPR8tys3svaSNl1Za3F0f2zuRHsvotsR/3FF1+EcbYlMxs7O3/Unsueq9FjFo1Lr+wiiVCyiyRCyS6SCCW7SCKU7CKJULKLJELJLpIIY/3Kdb0zsxEAn0+5qQLgVGkD+HqadWzNOi5AY6tVPcd2h7tPu/5bqcn+lTs363f3voYNINCsY2vWcQEaW63KGpt+jRdJhJJdJBGNTvatDb7/SLOOrVnHBWhstSplbA39m11EytPoV3YRKYmSXSQRDUl2M3vczP7PzA6b2QuNGEMeMxsws4/NbK+Z9Td4LNvMbNjM9k+5rcPMdprZp9nHaffYa9DYXjKzY9m122tmTzRobL1m9hczO2hmn5jZj7PbG3rtgnGVct1K/5vdzFoAHALwTwAGAXwA4Gl3P1DqQHKY2QCAPndv+AQMM9sI4CKA37n732W3/TuAM+7+cvaDst3d/7VJxvYSgIuN3sY7262oe+o24wCeBPDPaOC1C8b1fZRw3Rrxyr4BwGF3P+LuYwD+AGBzA8bR9Nz9PQBnbrl5M4Dt2efbMflkKV3O2JqCux939w+zzy8AuLnNeEOvXTCuUjQi2ZcDmLpmzyCaa793B/BnM9tjZlsaPZhpdLn7cWDyyQMgf32jxqDbeJfplm3Gm+ba1bL9eVGNSPbptpJqpvrfI+7+bQCbADyX/boq1alqG++yTLPNeFOodfvzohqR7IMApq7mtwJA/m6PJXP3oezjMIA30HxbUZ+8uYNu9jFeMbJEzbSN93TbjKMJrl0jtz9vRLJ/AOBuM7vTzG4D8AMAbzVgHF9hZm3ZGycwszYA30XzbUX9FoBnss+fAfBmA8fyJc2yjXfeNuNo8LVr+Pbn7l76PwBPYPId+c8A/FsjxpAzrtUAPsr+fdLosQF4DZO/1l3H5G9EzwJYAmAXgE+zjx1NNLb/BvAxgH2YTKzuBo3tHzD5p+E+AHuzf080+toF4yrlumm6rEgiNINOJBFKdpFEKNlFEqFkF0mEkl0kEUp2kUQo2UUS8f/J744D8Ij1uAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAALC0lEQVR4nO3dT4ic9R3H8c+nVi/qIWnGsMTQtZJDpdAoQyikiEUqMZfowWIOkoKwHhQUPFTsQY+hVKWHIqw1mBarCCrmEFpDEMSLcZQ0fxraqGw1ZslOyMF4stFvD/ukrHFnZ5zneeZ5st/3C5aZfWY2882Qd56d+c3M44gQgNXve00PAGAyiB1IgtiBJIgdSILYgSS+P8kbW7duXUxPT0/yJoFU5ubmdPbsWS93WanYbW+T9AdJV0j6U0TsXun609PT6vV6ZW4SwAq63e7Ay8b+Nd72FZL+KOlOSTdJ2mn7pnH/PAD1KvOYfYukDyPi44j4UtLLknZUMxaAqpWJfYOkT5d8f6rY9g22Z2z3bPf6/X6JmwNQRpnYl3sS4FuvvY2I2YjoRkS30+mUuDkAZZSJ/ZSkjUu+v17S6XLjAKhLmdjfk7TJ9g22r5J0r6R91YwFoGpjL71FxAXbD0n6uxaX3vZExPHKJgNQqVLr7BGxX9L+imYBUCNeLgskQexAEsQOJEHsQBLEDiRB7EASxA4kQexAEsQOJEHsQBLEDiRB7EASxA4kMdGPksbqc+jQoaZHGGjLli1Nj9Aq7NmBJIgdSILYgSSIHUiC2IEkiB1IgtiBJFhnT67N6+Rllfm7rcY1evbsQBLEDiRB7EASxA4kQexAEsQOJEHsQBKss69ybV5HH7aW3ebZL0elYrc9J+m8pK8kXYiIbhVDAaheFXv2X0TE2Qr+HAA14jE7kETZ2EPSm7bftz2z3BVsz9ju2e71+/2SNwdgXGVj3xoRt0i6U9KDtm+99AoRMRsR3YjodjqdkjcHYFylYo+I08XpgqTXJa2+twoBq8TYsdu+2va1F89LukPSsaoGA1CtMs/Gr5f0uu2Lf85fI+JvlUyFVaPM+8JZh6/W2LFHxMeSflrhLABqxNIbkASxA0kQO5AEsQNJEDuQBG9xXQWaXIJq8iOXV+PHPdeJPTuQBLEDSRA7kASxA0kQO5AEsQNJEDuQBOvsl4Gs6+ioFnt2IAliB5IgdiAJYgeSIHYgCWIHkiB2IAnW2VuAdXRMAnt2IAliB5IgdiAJYgeSIHYgCWIHkiB2IAliB5IYGrvtPbYXbB9bsm2t7QO2Txana+odE0BZo+zZX5C07ZJtj0k6GBGbJB0svgfQYkNjj4i3JZ27ZPMOSXuL83sl3VXtWACqNu5j9vURMS9Jxel1g65oe8Z2z3av3++PeXMAyqr9CbqImI2IbkR0O51O3TcHYIBxYz9je0qSitOF6kYCUIdxY98naVdxfpekN6oZB0Bdhr6f3fZLkm6TtM72KUlPSNot6RXb90v6RNI9dQ55ueP96miDobFHxM4BF91e8SwAasQr6IAkiB1IgtiBJIgdSILYgSSIHUiC2IEkiB1IgtiBJIgdSILYgSSIHUiC2IEkOGTzKsDbWDEK9uxAEsQOJEHsQBLEDiRB7EASxA4kQexAEsQOJEHsQBLEDiRB7EASxA4kQexAEsQOJEHsQBLEDiQxNHbbe2wv2D62ZNuTtj+zfbj42l7vmADKGmXP/oKkbctsfyYiNhdf+6sdC0DVhsYeEW9LOjeBWQDUqMxj9odsHyl+zV8z6Eq2Z2z3bPf6/X6JmwNQxrixPyvpRkmbJc1LemrQFSNiNiK6EdHtdDpj3hyAssaKPSLORMRXEfG1pOck8fGmQMuNFbvtqSXf3i3p2KDrAmiHoZ8bb/slSbdJWmf7lKQnJN1me7OkkDQn6YH6Rmy/Q4cONT1CY9r8d+fz9L9paOwRsXOZzc/XMAuAGvEKOiAJYgeSIHYgCWIHkiB2IAkO2Zxcm5fOylrp75ZxWY49O5AEsQNJEDuQBLEDSRA7kASxA0kQO5AE6+wVGLZmW/da9mpeK6/LsPtsNa7Ds2cHkiB2IAliB5IgdiAJYgeSIHYgCWIHkiB2IAliB5IgdiAJYgeSIHYgCWIHkiB2IAliB5Lg/exYUdn3dfNe+/YYume3vdH2W7ZP2D5u++Fi+1rbB2yfLE7X1D8ugHGN8mv8BUmPRsSPJf1M0oO2b5L0mKSDEbFJ0sHiewAtNTT2iJiPiA+K8+clnZC0QdIOSXuLq+2VdFdNMwKowHd6gs72tKSbJb0raX1EzEuL/yFIum7Az8zY7tnu9fv9kuMCGNfIsdu+RtKrkh6JiM9H/bmImI2IbkR0O53OODMCqMBIsdu+UouhvxgRrxWbz9ieKi6fkrRQz4gAqjB06c22JT0v6UREPL3kon2SdknaXZy+UcuEq0DTHzVdRptnK2M1flT0MKOss2+VdJ+ko7YPF9se12Lkr9i+X9Inku6pZUIAlRgae0S8I8kDLr692nEA1IWXywJJEDuQBLEDSRA7kASxA0nwFtcWuJzX4ZuUca28DPbsQBLEDiRB7EASxA4kQexAEsQOJEHsQBKss18GWE9GFdizA0kQO5AEsQNJEDuQBLEDSRA7kASxA0kQO5AEsQNJEDuQBLEDSRA7kASxA0kQO5AEsQNJDI3d9kbbb9k+Yfu47YeL7U/a/sz24eJre/3jAhjXKB9ecUHSoxHxge1rJb1v+0Bx2TMR8fv6xgNQlVGOzz4vab44f972CUkb6h4MQLW+02N229OSbpb0brHpIdtHbO+xvWbAz8zY7tnu9fv9ctMCGNvIsdu+RtKrkh6JiM8lPSvpRkmbtbjnf2q5n4uI2YjoRkS30+mUnxjAWEaK3faVWgz9xYh4TZIi4kxEfBURX0t6ThKfigi02CjPxlvS85JORMTTS7ZPLbna3ZKOVT8egKqM8mz8Vkn3STpq+3Cx7XFJO21vlhSS5iQ9UMN8ACoyyrPx70jyMhftr34cAHXhFXRAEsQOJEHsQBLEDiRB7EASxA4kQexAEsQOJEHsQBLEDiRB7EASxA4kQexAEsQOJOGImNyN2X1J/1myaZ2ksxMb4Ltp62xtnUtitnFVOdsPI2LZz3+baOzfunG7FxHdxgZYQVtna+tcErONa1Kz8Ws8kASxA0k0Hftsw7e/krbO1ta5JGYb10Rma/QxO4DJaXrPDmBCiB1IopHYbW+z/S/bH9p+rIkZBrE9Z/tocRjqXsOz7LG9YPvYkm1rbR+wfbI4XfYYew3N1orDeK9wmPFG77umD38+8cfstq+Q9G9Jv5R0StJ7knZGxD8nOsgAtuckdSOi8Rdg2L5V0heS/hwRPym2/U7SuYjYXfxHuSYiftOS2Z6U9EXTh/EujlY0tfQw45LukvRrNXjfrTDXrzSB+62JPfsWSR9GxMcR8aWklyXtaGCO1ouItyWdu2TzDkl7i/N7tfiPZeIGzNYKETEfER8U589LuniY8UbvuxXmmogmYt8g6dMl359Su473HpLetP2+7Zmmh1nG+oiYlxb/8Ui6ruF5LjX0MN6TdMlhxltz341z+POymoh9uUNJtWn9b2tE3CLpTkkPFr+uYjQjHcZ7UpY5zHgrjHv487KaiP2UpI1Lvr9e0ukG5lhWRJwuThckva72HYr6zMUj6BanCw3P839tOoz3cocZVwvuuyYPf95E7O9J2mT7BttXSbpX0r4G5vgW21cXT5zI9tWS7lD7DkW9T9Ku4vwuSW80OMs3tOUw3oMOM66G77vGD38eERP/krRdi8/IfyTpt03MMGCuH0n6R/F1vOnZJL2kxV/r/qvF34jul/QDSQclnSxO17Zotr9IOirpiBbDmmpotp9r8aHhEUmHi6/tTd93K8w1kfuNl8sCSfAKOiAJYgeSIHYgCWIHkiB2IAliB5IgdiCJ/wH4aYKwZWL9CQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(1) \n",
    "n=600\n",
    "X = np.random.randn(n,784)/np.sqrt(n)\n",
    "E = 0.05 * np.random.randn(n)\n",
    "Y = X @ I + E\n",
    "Lambda=np.sqrt(range(1,785))-np.sqrt(range(0,784))\n",
    "gmax=norme_duale(X.T @ Y, Lambda)\n",
    "G = np.exp(np.linspace(-3*np.log(gmax),np.log(gmax),100))\n",
    "EQM = np.zeros(100)\n",
    "i=0\n",
    "for g in G:\n",
    "    Sol,_ = FISTA(X,Y,g*Lambda)\n",
    "    EQM[i] = np.linalg.norm(Sol-I,ord=2)**2\n",
    "    i+=1\n",
    "i0 = np.argmin(EQM)\n",
    "Sol0,_ = FISTA(X,Y,G[i0]*Lambda)\n",
    "M=max(np.abs(Sol0))\n",
    "\n",
    "#10.527\n",
    "# 7.0465659212\n",
    "ratio = 7.0465659212/10.527\n",
    "\n",
    "tau=ratio*norme_duale(Sol0,Lambda)\n",
    "prox_Sol0 = prox_L1_ordonnee(Sol0,tau*Lambda)\n",
    "\n",
    "image_Sol0 = Sol0.reshape(28, 28)\n",
    "plt.imshow(np.abs(image_Sol0), cmap=\"Greys\",vmin=0,vmax=M)\n",
    "plt.show()\n",
    "plt.close()\n",
    "\n",
    "image_prox_Sol0 = prox_Sol0.reshape(28, 28)\n",
    "plt.imshow(np.abs(image_prox_Sol0), cmap=\"Greys\",vmin=0,vmax=M)\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}