{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1cea692b",
   "metadata": {},
   "source": [
    "# VIII.2 - Travaux pratiques\n",
    "\n",
    "Dans cette séance de travaux pratiques, l'objectif est de vous faire implémenter par vous-même l'apprentissage d'un premier réseau de neurones simple. En pratique, il existe de nombreuses bibliothèques qui implémentent à notre place les routines qui réalisent cet apprentissage. Toutefois, il est très formateur de manipuler les concepts de base au moins une fois pour comprendre comment fonctionnent les outils que nous utiliseront plus tard, comme Keras.\n",
    "\n",
    "## Jeu de données MNIST\n",
    "\n",
    "Pour cette séance, nous travaillerons avec la base de données MNIST. Elle est constituée de 70000 images (60000 images en apprentissage, 10000 en test) de chiffres manuscrits en noir et blanc. L'objectif est d'identifier automatiquement le chiffre à partir de l'image.\n",
    "\n",
    "Commençons par importer les données. Comme il s'agit d'un jeu de données très classique et courant, il est disponible de base dans de nombreuses bibliothèques, comme Keras. Cela nous permet de l'importer directement en une ligne de code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e683a2f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from keras.datasets import mnist\n",
    "# Import de MNIST depuis Keras\n",
    "(X_train, y_train), (X_test, y_test) = mnist.load_data()\n",
    "# Transformation des images 28x28 en vecteur de dimension 784\n",
    "X_train = X_train.reshape(60000, 784).astype('float32')\n",
    "X_test = X_test.reshape(10000, 784).astype('float32')\n",
    "# Normalisation entre 0 et 1\n",
    "X_train /= 255\n",
    "X_test /= 255\n",
    "\n",
    "# Affichage du nombre de'exemples\n",
    "print(f\"{X_train.shape[0]} exemples d'apprentissage\")\n",
    "print(f\"{X_test.shape[0]} exemples de test\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ab16c3a",
   "metadata": {},
   "source": [
    "<div class=\"admonition-question admonition\">\n",
    "\n",
    "**Question**\n",
    "\n",
    "Afficher à l'aide de matplotlib les premières images du jeu d'apprentissage. La fonction `plt.imshow()` (cf. [sa documentation](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.imshow.html)) peut vous être utile.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0373bbc8",
   "metadata": {},
   "source": [
    "<div class=\"admonition-question admonition\">\n",
    "\n",
    "**Question**\n",
    "\n",
    "Quel est l'espace dans lequel se trouvent les images ? Quelle est sa dimension ?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "931ef081",
   "metadata": {},
   "source": [
    "## Régression logistique\n",
    "\n",
    "### Modèle de prédiction\n",
    "\n",
    "Nous allons implémenter un modèle de classification linéaire simple : la régression logistique. Concrètement, la régression logistique est équivalente à un réseau de neurones à une seule couche. Il s'agit d'une projection du vecteur d'entrée $\\mathbf{x_i}$ par un vecteur de paramètres $\\mathbf{w_{c}}$, plus un biais sclaaire $b_c$, pour chaque classe. Le schéma ci-dessous illustre le modèle de régression logistique avec un réseau de neurones.\n",
    "\n",
    "![](figures/LR.png)\n",
    "\n",
    "En l'occurrence, pour MNIST $\\mathbf{x}_i$ est de dimension 784 et il y a dix chiffres possibles, donc 10 classes différentes. Dans notre cas, on considère que l'image d'entrée est représentée sous sa forme « aplatie », c'est-à-dire un vecteur (1, 784).\n",
    "\n",
    "Pour simplifier les notations, on regroupe l'ensemble des jeux de paramètres $\\mathbf{w_{c}}$ pour les 10 classes possibles dans une unique matrice $\\mathbf{W}$ de dimensions $784\\times 10$. De la même façon, les biais sont regroupés dans un vecteur $\\mathbf{b}$ de longueur 10. La sortie de la régression logistique est un vecteur contenant une activation pour chaque classe, c'est-à-dire $\\mathbf{\\hat{s_i}} =\\mathbf{x_i} \\mathbf{W} + \\mathbf{b}$ de dimensions (1, 10).\n",
    "\n",
    "Afin de transformer les activations en de sortie en probabilités pour une distribution catégorielle, on ajoute une fonction d'activation de *softmax* sur $\\mathbf{\\hat{y_i}} = \\sigma(\\mathbf{s_i})$. Cela nous permet d'obtenir en sortie un vecteur de prédictions $\\mathbf{\\hat{y_i}}$, de dimensions (1, 10), qui représente la probabilité *a posteriori* $p(\\mathbf{\\hat{y_i}} | \\mathbf{x_i})$ pour chacune des 10 classes :\n",
    "\n",
    "$$\n",
    "p(\\hat{y}_{c,i} | \\mathbf{x_i}) ) = \\frac{e^{\\langle \\mathbf{x_i} ; \\mathbf{w_{c}}\\rangle + b_{c}}}{\\sum_{c'=1}^{10} e^{\\langle \\mathbf{x_i} ; \\mathbf{w_{c'}}\\rangle + b_{c'}}}\n",
    "$$\n",
    "\n",
    "<div class=\"admonition-question admonition\">\n",
    "\n",
    "**Question**\n",
    "\n",
    "Quel est le nombre de paramètres du modèle utilisé ? Justifier le calcul.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a099941",
   "metadata": {},
   "source": [
    "### Formulation du problème d'apprentissage\n",
    "\n",
    "Pour entraîner le réseau de neurones, c'est-à-dire déterminer les valeurs optimales des paramètres $\\mathbf{W}$ et $\\mathbf{b}$, on va comparer pour chaque exemple d'apprentissage la sortie prédite $\\mathbf{\\hat{y_i}}$ (équation (softmax)) pour l'image $\\mathbf{x_i}$ à la sortie réelle $\\mathbf{y_i^*}$ (vérité terrain issue de la supervision). Dans notre cas, on choisit d'encoder la catégorie de l'image $\\mathbf{x_i}$ sous forme *one-hot*, c'est-à-dire :\n",
    "\n",
    "$$\n",
    "y_{c,i}^* =\n",
    "\\begin{cases}\n",
    "  1 & \\text{si c correspond à l'indice de la classe de } \\mathbf{x_i}  \\\\\n",
    "  0 & \\text{sinon}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Générons les étiquettes (*labels*) au format *one-hot* ((one-hot)) à l'aide de la fonction `to_categorical` (cf. [documentation de Keras](https://www.tensorflow.org/api_docs/python/tf/keras/utils/to_categorical))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4414da51",
   "metadata": {},
   "outputs": [],
   "source": [
    "from keras.utils import to_categorical\n",
    "n_classes = 10\n",
    "# Conversion des étiquettes (int) au format vectoriel one-hot\n",
    "Y_train = to_categorical(y_train, n_classes)\n",
    "Y_test = to_categorical(y_test, n_classes)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abdc7ac8",
   "metadata": {},
   "source": [
    "L'erreur de prédiction sera définie à l'aide de l'entropie croisée (*cross-entropy*). Cette fonction de coût s'applique entre $\\mathbf{\\hat{y_i}}$ et $\\mathbf{y_i^*}$ par la formule :\n",
    "$𝓛(\\mathbf{\\hat{y_i}}, \\mathbf{y_i^*}) = -\\sum_{c=1}^{10} y_{c,i}^* \\log(\\hat{y}_{c,i}) = - \\log(\\hat{y}_{c^*,i})$, où $c^*$ correspond à l'indice de la classe donnée par la supervision pour l'image $\\mathbf{x_i}$.\n",
    "\n",
    "<div class=\"admonition note\">\n",
    "\n",
    "**Note**\n",
    "\n",
    "L'entropie croisée correspond en réalité à la divergence de Kullback-Leiber pour des distributions catégorielles. La divergence KL est une mesure de dissimilarité entre distributions de probabilité. Autrement dit, l'erreur que l'on mesure vise à réduire l'écart entre la distribution réelle des catégories et la distribution prédite.\n",
    "\n",
    "</div>\n",
    "\n",
    "La fonction de coût finale correspond à l'erreur d'apprentissage, c'est-à-dire la moyenne l'entropie croisée sur l'ensemble de la base d'apprentissage $𝓓$ constituée des $N=60000$ images :\n",
    "\n",
    "$$\n",
    "𝓛_{\\mathbf{W},\\mathbf{b}}(𝓓)  = - \\frac{1}{N}\\sum_{i=1}^{N} \\log(\\hat{y}_{c^*,i})\n",
    "$$\n",
    "\n",
    "<div class=\"admonition-question admonition\">\n",
    "\n",
    "**Question**\n",
    "\n",
    "La fonction de coût de (CE) est-elle convexe par rapport aux paramètres $\\mathbf{W}$, $\\mathbf{b}$ du modèle ? Avec un pas de gradient bien choisi, peut-on assurer la convergence vers le minimum global de la solution ?\n",
    "\n",
    "</div>\n",
    "\n",
    "### Optimisation du modèle\n",
    "\n",
    "Nous allons minimiser la fonction de coût à l'aide de l'algorithme de descente de gradient appliqué sur les paramètres $\\mathbf{W}$ et $\\mathbf{b}$ du modèle de régression logistique. Pour ce faire, nous allons avoir besoin des gradients de l'entropie croisée par rapport à $\\mathbf{W}$ ainsi que $\\mathbf{b}$. Nous pouvons nous appuyer sur la des dérivées chaînées (*chain rule*, ou théorème de dérivation des fonctions composées) :\n",
    "\n",
    "$$\n",
    "\\frac{\\partial 𝓛}{\\partial \\mathbf{W}} =  \\frac{1}{N}\\sum_{i=1}^{N} \\frac{\\partial 𝓛}{\\partial \\mathbf{\\hat{y_i}}}  \\frac{\\partial \\mathbf{\\hat{y_i}}}{\\partial \\mathbf{s_i}} \\frac{\\partial \\mathbf{s_i}}{\\partial \\mathbf{W}}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\frac{\\partial 𝓛}{\\partial \\mathbf{b}} =  \\frac{1}{N}\\sum_{i=1}^{N} \\frac{\\partial 𝓛}{\\partial \\mathbf{\\hat{y_i}}}  \\frac{\\partial \\mathbf{\\hat{y_i}}}{\\partial \\mathbf{s_i}} \\frac{\\partial \\mathbf{s_i}}{\\partial \\mathbf{b}}\n",
    "$$\n",
    "\n",
    "<div class=\"admonition-question admonition\">\n",
    "\n",
    "**Question**\n",
    "\n",
    "Montrer que :\n",
    "\n",
    "$$\n",
    "\\frac{\\partial 𝓛}{\\partial \\mathbf{s_i}} = \\mathbf{\\delta^y_i} = \\frac{\\partial 𝓛}{\\partial \\mathbf{\\hat{y_i}}}  \\frac{\\partial \\mathbf{\\hat{y_i}}}{\\partial \\mathbf{s_i}} = \\mathbf{\\hat{y_i}} - \\mathbf{y_i^*}\n",
    "$$\n",
    "\n",
    "</div>\n",
    "\n",
    "<div class=\"admonition-question admonition\">\n",
    "\n",
    "**Question**\n",
    "\n",
    "En déduire que les gradients de $𝓛$ par rapport aux paramètres du modèle s'écrivent :\n",
    "\n",
    "$$\n",
    "\\frac{\\partial 𝓛}{\\partial \\mathbf{W}} = \\frac{1}{N} \\mathbf{X}^T (\\mathbf{\\hat{Y}} - \\mathbf{Y^*}) = \\frac{1}{N} \\mathbf{X}^T \\mathbf{\\Delta^y}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\frac{\\partial 𝓛}{\\partial \\mathbf{b}} = \\frac{1}{N}\\sum_{i=1}^{N}(\\mathbf{\\hat{y_i}} - \\mathbf{y_i^*})\n",
    "$$\n",
    "\n",
    "où $\\mathbf{X}$ est la matrice des données (taille $60000\\times 784$), $\\mathbf{\\hat{Y}}$ est la matrice des labels prédits sur l'ensemble de la base d'apprentissage (taille $60000\\times 10$) et $\\mathbf{Y^*}$ est la matrice des labels issue de la supervision (*ground truth*, taille $60000\\times 10$), et $\\mathbf{\\Delta^y}=\\mathbf{\\hat{Y}}-\\mathbf{Y^*}$.\n",
    "\n",
    "</div>\n",
    "\n",
    "### Implémentation de l'apprentissage\n",
    "\n",
    "Les gradients obtenus par les équations (gradientW) et (gradientB) s'écrivent sous forme « vectorielle », ce qui rend les calculs efficaces avec des bibliothèques de calcul scientifique telles que `numpy`. Après calcul du gradient, les paramètres sont mis à jour de la façon suivante :\n",
    "\n",
    "$$\n",
    "\\mathbf{W}^{(t+1)} = \\mathbf{W}^{(t)} - \\eta \\frac{\\partial 𝓛}{\\partial \\mathbf{W}}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbf{b}^{(t+1)} = \\mathbf{b}^{(t)} - \\eta \\frac{\\partial 𝓛}{\\partial \\mathbf{b}}\n",
    "$$\n",
    "\n",
    "où $\\eta$ est le pas de gradient (*learning rate*).\n",
    "\n",
    "En théorie, la descente de gradient nécessite de calculer les gradients de la fonction de coût sur tout le jeu de données d'apprentissage. Toutefois, ce jeu de données est assez grand et les gradients peuvent être longs à calculer. En pratique, on implémente plutôt une descente de gradient *stochastique*, c'est à dire que les gradients aux équations (gradientW) et (gradientB) ne seront pas calculés sur l'ensemble des $N=60000$ images d'apprentissage, mais sur un sous-ensemble de $n$ images appelé *batch* ou *lot*. Cette technique permet une mise à jour des paramètres plus fréquente qu'avec une descente de gradient classique, un temps de calcul réduit et une convergence plus rapide, au détriment d'une approximation du gradient.\n",
    "\n",
    "Le code ci-dessous décrit le squelette de l'algorithme de descente de gradient qui va permettre l'optimisation des paramètres du modèle :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7cbcf78d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "N, d = X_train.shape # N exemples, dimension d\n",
    "W = np.zeros((d, n_classes)) # initialisation de poids\n",
    "b = np.zeros((1, n_classes)) # initialisation des biais\n",
    "\n",
    "n_epochs = 20 # Nombre d'epochs de la descente de gradient\n",
    "eta = 1e-1 # Learning rate (pas d'apprentissage)\n",
    "batch_size = 100 # Taille du lot\n",
    "n_batches = int(float(N) / batch_size)\n",
    "\n",
    "# On alloue deux matrices pour stocker les valeurs des gradients\n",
    "gradW = np.zeros((d, n_classes))\n",
    "gradb = np.zeros((1, n_classes))\n",
    "\n",
    "for epoch in range(n_epochs):\n",
    "    for batch_idx in range(n_batches):\n",
    "        # ********* À compléter **********\n",
    "        # FORWARD PASS : calculer la prédiction y à partir des paramètres courants pour les images du batch\n",
    "        # BACKWARD PASS :\n",
    "        # 1) calculer les gradients de l'erreur par rapport à W et b\n",
    "        # 2) mettre à jour les paramètres W et b selon la descente de gradient\n",
    "        ..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12a782a8",
   "metadata": {},
   "source": [
    "<div class=\"admonition-question admonition\">\n",
    "\n",
    "**Question**\n",
    "\n",
    "Compléter ce code. Vous devez notamment :\n",
    "\n",
    "- Écrire une fonction `forward(batch, W, b)` qui calcule la prédiction (vecteur de sortie $\\hat{\\mathbf{y}}$) pour chaque exemple d'un batch de données. Si on considère un batch des données de taille $tb\\times 784$, les paramètres $\\mathbf{W}$ (taille $784\\times 10$) et $\\mathbf{b}$ (taille $1\\times 10$), la fonction `forward` renvoie la prédiction $\\mathbf{\\hat{Y}}$ sur le batch (taille $tb\\times 10$). La fonction `forward` sera appelée pour chaque itération de la double boucle précédente.\n",
    "- Utiliser la fonction `softmax` ci-dessous pour calculer le résultat du passage du softmax sur chaque élément de la matrice de la projection linéaire (taille $tb\\times 10$) :\n",
    "\n",
    "  ```python\n",
    "  def softmax(X):\n",
    "      # Entrée: matrice X de dimensions batch x d\n",
    "      # Sortie: matrice de mêmes dimensions\n",
    "      E = np.exp(X)\n",
    "      return (E.T / np.sum(E,axis=1)).T\n",
    "  ```\n",
    "\n",
    "- Compléter le code pour la passe *backward*, consistant à :\n",
    "  - Calculer les gradients comme indiqué dans les équations (gradientW) et (gradientB).\n",
    "  - Mettre à jour les paramètres par descente de gradient comme indiqué dans les équations (gradientupdateW) et (gradientupdateB).\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcc4938d",
   "metadata": {},
   "source": [
    "<div class=\"admonition-question admonition\">\n",
    "\n",
    "**Question**\n",
    "\n",
    "Évaluer les performances du modèle de régression logistique entraîné sur MNIST. On utilisera le taux de bonne classification (*accuracy*) comme métrique. Commencer par mesurer l'évolution des performances du modèle au cours de l'apprentissage (calcul de l'*accuracy* à chaque époque), puis évaluer sur le modèle sur la base de test. Vous pouvez utiliser la [fonction de scikit-learn](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.accuracy_score.html) ou la fonction `accuracy` ci-dessous (qui effectue également la phase de prédiction).\n",
    "\n",
    "**Vous devriez obtenir un score de l'ordre de 92% sur la base de test pour ce modèle de régression logistique.**\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "341dca3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def accuracy(W, b, images, labels):\n",
    "    \"\"\" W: matrice de paramètres\n",
    "        b: vecteur de biais\n",
    "        images: images de MNIST\n",
    "        labels: étiquettes de MNIST pour les images\n",
    "\n",
    "        Renvoie l'accuracy du modèle (W, b) sur les images par rapport aux labels\n",
    "    \"\"\"\n",
    "    pred = forward(images, W, b)\n",
    "    return np.where(pred.argmax(axis=1) != labels.argmax(axis=1), 0.,1.).mean()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
