{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "a308ef0e",
   "metadata": {},
   "source": [
    "# TP ACP 2\n",
    "\n",
    "Dans ce TP, nous allons étendre notre résumé à une dimension fait au TP précédent à plusieurs dimensions. Nous allons également reconstruire les données compressées, et appliquer cette reconstruction à des images. \n",
    "\n",
    "### Compétences associées :\n",
    "\n",
    "| Numéro        | Compétence           | \n",
    "|:------------- |:--------------------:|\n",
    "|AS405|\tSavoir déterminer la matrice de projection au sens de l’ACP|\n",
    "|AS406|\tSavoir déterminer la projection des données sur un sous espace vectoriel au sens de l’ACP|\n",
    "|AS407|\tSavoir interpréter les valeurs propres en fonction de chaque axe factoriel|\n",
    "|AS408|\tSavoir expliciter la transformation inverse pour reconstruire les données|\n",
    "|AS409|\tConnaître les limites de l’ACP|\n",
    "|PY306|\tSavoir calculer la variance correspondant à k premiers axes factoriels|\n",
    "|PY307|\tSavoir calculer une ACP en d dimensions sur un jeu de données|\n",
    "|PY308|\tÊtre capable de reconstruire une image après ACP|\n",
    "\n",
    "\n",
    "\n",
    "## Ex1 : Information expliquée \n",
    "\n",
    "1) Chargez le dataset Iris et normalisez les données. Gardez les 'targets', c'est à dire les noms d'espèces de chaque plante dans un vecteur `y`. Vous pouvez utiliser l'option `return_X_y=True` de la fonction `load_iris`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25d6dfdf",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_iris\n",
    "import numpy as np\n",
    "X,y = load_iris(return_X_y=True)\n",
    "n,p = X.shape\n",
    "mu,std = np.mean(X,axis=0), np.std(X,axis=0)\n",
    "X = (X- mu)/std"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "811efb69",
   "metadata": {},
   "source": [
    "2. Plottez les données `X` en 2D (les variables au choix) et affichez chaque point avec une couleur par espèce. Vous pouvez utiliser le paramètre `c=y` pour y arriver facilement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a8b78b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.scatter(X[:,0],X[:,1],c=y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "caba152d",
   "metadata": {},
   "source": [
    "Nous allons calculer l'ACP de ces données. Vous pouvez aller vite en utilisant `sklearn.decomposition.PCA` mais il est également important de savoir coder la méthode pour bien la comprendre. Dans les questions suivantes, nous allons détailler pas à pas commment calculer la PCA avec `numpy`. La version `sklearn` est donnée plus loin. \n",
    "\n",
    "3.Calculez la décomposition valeurs propres/vecteurs propres de la matrice de variance/covariance associées aux données normalisées."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19f073ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "cov = (X.T@X)/n # La division par n n'est pas obligatoire. Juste ne pas oublier dans l'interprétation des variances \n",
    "w,v = np.linalg.eigh(cov)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d13d66ef",
   "metadata": {},
   "source": [
    "4. Triez les valeurs propres dans l'ordre décroissant et appliquez le même ordre de tri aux colonnes de la matrice des vecteurs propres.\n",
    "\n",
    "**Aide** : \n",
    "\n",
    "    1. On pourra utiliser la fonction `argsort` de `numpy`. [doc](https://numpy.org/doc/stable/reference/generated/numpy.argsort.html)\n",
    "    2. L'utilisation des slices avec `-1` en dernière position permet de parcourir un élément de la fin au début. cf exemple."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68b6b2d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([3,1,5,2,4]) # un vecteur\n",
    "print(a)\n",
    "print(a[::-1]) # un vecteur affiché \"à l'envers\"\n",
    "print(np.argsort(a)) # affiche les positions de chaque élement dans l'ordre croissant\n",
    "sorted_indexes = np.argsort(a)\n",
    "print(a[sorted_indexes]) #applique le tri \"sorted_indexes\" à \"a\"\n",
    "print(np.argsort(a)[::-1]) # affiche les positions de chaque élement dans l'ordre décroissant\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9a690ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "indexes = np.argsort(w)[::-1]\n",
    "w=w[indexes]\n",
    "v=v[:,indexes]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09984b50",
   "metadata": {},
   "source": [
    "5. Calculez la matrice de projection $P$ de $R^4$ vers $R^2$, en prenant les deux premiers axes de l'ACP. Projetez les données $X$ avec $P$. Maintenant que nous sommes en 2D, plottez les données avec une couleur par espèce. Comparez avec le premier plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0daa31b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "P = v[:,:2]\n",
    "Xproj = X@P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ba9d780",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.subplot(121)\n",
    "plt.scatter(X[:,0],X[:,1],c=y)\n",
    "plt.subplot(122)\n",
    "plt.scatter(Xproj[:,0],Xproj[:,1],c=y)\n",
    "#meilleure séparation même si pas parfaite."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2f993ef",
   "metadata": {},
   "source": [
    "6. Comme vu en cours, la variance expliquée par chaque axe est donnée par la  valeur propre associée à ce même axe. Calculez le pourcentage de variance expliqué par chaque axe. On pourra afficher cette information avec `plt.bar(np.arange(4), ...)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b70b2f20",
   "metadata": {},
   "outputs": [],
   "source": [
    "pourcent_var = (w/w.sum())*100\n",
    "print(pourcent_var)\n",
    "plt.bar(np.arange(4),height=pourcent_var)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d3cb7e9",
   "metadata": {},
   "source": [
    "7. Enfin, affichez le pourcentage cumulé de variance expliqué en prenant les $p$ premiers axes. Que pouvez vous en conclure ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "781c4050",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "plt.bar(np.arange(4),height=np.cumsum(pourcent_var))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd42f471",
   "metadata": {},
   "source": [
    "**Bonus** : Faites la même analyse sur Diabetes ou d'autres datasets."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "506b269a",
   "metadata": {},
   "source": [
    "**Bonus 2**  : la même chose avec sklearn :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4aa32d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_iris\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "X,y = load_iris(return_X_y=True)\n",
    "scaler = StandardScaler()\n",
    "#on normalise nos données\n",
    "X_norm = scaler.fit_transform(X) \n",
    "#on calcule la PCA\n",
    "pca = PCA() #ici on peut spécifier le nb de composants. Si on a pas besoin de tout, cela peut nous faire gagner du temps de calcul\n",
    "pca.fit(X_norm)\n",
    "X_proj =pca.transform(X_norm)\n",
    "plt.scatter(X_proj[:,0],X_proj[:,1],c=y)\n",
    "\n",
    "plt.figure()\n",
    "plt.bar(np.arange(4),pca.explained_variance_ratio_*100)\n",
    "plt.figure()\n",
    "plt.bar(np.arange(4),np.cumsum(pca.explained_variance_ratio_)*100)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0e232c4",
   "metadata": {},
   "source": [
    "## Ex2 : Reconstruction\n",
    "Ici nous allons étudier la reconstruction d'une image suite à une ACP, et après avoir gardé qu'une partie des composantes principales. Cet exercice nous permettra de mettre en évidence les capacités de compression des méthodes basées sur l'ACP.\n",
    "\n",
    "1) Charger l'image du cameraman, une image bien connue pour le traitement d'images. Cette image étant en niveaux de gris, nous allons pouvoir la représenter sous forme d'une matrice `X` $\\in \\mathbb{R}^{256 \\times 256}$. Ici, nous faisons l'hypothèse que chaque ligne de l'image est une observation, et que chaque colonne est une variable. Vous pouvez afficher l'image avec `plt.imshow(X,cmap=\"gray\")`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9bb37cf5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import imageio.v2 as imageio\n",
    "import matplotlib.pyplot as plt\n",
    "X = imageio.imread('cameraman.png')\n",
    "plt.imshow(X,cmap=\"gray\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac5d07a1",
   "metadata": {},
   "source": [
    "2. Calculez l'ACP (en utilisant votre code ou `sklearn`) de l'image codée par `X`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb7aa9fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "mu,std = np.mean(X,axis=0), np.std(X,axis=0)\n",
    "Xnorm = (X - mu ) / std\n",
    "cov = Xnorm.T@Xnorm\n",
    "w,v = np.linalg.eigh(cov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6ac2bcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "indexes = w.argsort()[::-1]\n",
    "w=w[indexes]\n",
    "v=v[:,indexes]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bdbacbd",
   "metadata": {},
   "source": [
    "3. Calculez une projection de l'image en dimension 100 via votre ACP. Vous devriez pobtenir une matrice de taille $256 \\times 100$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b779530",
   "metadata": {},
   "outputs": [],
   "source": [
    "P = v[:,:100]\n",
    "Xproj = Xnorm@P\n",
    "plt.matshow(Xproj) #on voit rien !"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d015861",
   "metadata": {},
   "source": [
    "4. Quelle est le pourcentage de variance expliquée par votre projection ? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd353f65",
   "metadata": {},
   "outputs": [],
   "source": [
    "print((np.sum(w[:100]) / np.sum(w))*100)\n",
    "#on reduisant la taille par 2.5, on a 99.5% de la variance expliquée"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60c86b7d",
   "metadata": {},
   "source": [
    "5. Reconstruisons l'image pour revenir dans le domaine d'origine. Pensez à bien rajouter la moyenne et multiplier par l'écart type pour retrouver le même ordre de grandeur des données. \n",
    "\n",
    "** NB** : La projection inverse $P^{-1}$ est donnée par $P^\\top$ étant donné que $P$ est orthogonale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f396fa3",
   "metadata": {},
   "outputs": [],
   "source": [
    "Xreconstruit = Xproj@P.T\n",
    "Xreconstruit = Xreconstruit *  std + mu\n",
    "plt.imshow(Xreconstruit,cmap=\"gray\")\n",
    "plt.figure()\n",
    "plt.imshow(X,cmap=\"gray\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8728ddc8",
   "metadata": {},
   "source": [
    "6. Essayez différentes valeurs de dimensions intermédiaires jusqu'à trouver un bon équilibre en taille et qualité de reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c751cc0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "d= 20\n",
    "P = v[:,:d]\n",
    "Xproj = Xnorm@P\n",
    "Xreconstruit = Xproj@P.T\n",
    "Xreconstruit = Xreconstruit *  std + mu\n",
    "plt.subplot(121)\n",
    "plt.imshow(Xreconstruit,cmap=\"gray\")\n",
    "plt.subplot(122)\n",
    "plt.imshow(X,cmap=\"gray\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2aa2b32e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mlenv",
   "language": "python",
   "name": "mlenv"
  },
  "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}