{
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4,
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Programmer le gradient conjugué\n",
    "\n",
    "A partir de ce [cours sur le gradient conjugué](http://perso.unifr.ch/ales.janka/numeroptim/07_conjgrad.pdf) programmez en Python + Numpy le gradient conjugué en exploitant les astuces mathématiques indiquées pour optimiser\n",
    "votre code.\n",
    "\n",
    "* Effectuez des tests pour valider votre code. \n",
    "* Comparez la vitesse de convergence à celle du gradient avec μ optimal. Tracez des courbes de convergence (cf la feuille qui en parle)\n",
    "* Comparez les temps de calcul.\n",
    "\n",
    "\n",
    "Note : Veuillez écrire des fonctions les plus propres possibles, en particulier qui n'utilisent pas des variables globales comme c'est le cas dans ma correction du gradient (ma33)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as lin\n",
    "import matplotlib.pylab as plt\n",
    "\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_system(N):\n",
    "    A = 1.0 * np.random.randint(-10, 10, size=(N,N))\n",
    "    A[np.diag_indices(N)] = 0.1 + np.abs(A).sum(axis=0)    # diag dominante\n",
    "    A = A + A.T                                            # symétrique\n",
    "    A = A / np.abs(A).sum(axis=0).mean()\n",
    "    b = 1.0 * np.random.randint(-10,10,size=(N))\n",
    "    return A, b\n",
    "\n",
    "A, b = make_system(2)\n",
    "print(A, \"\\n\\n\", b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradient_conjugué(A, x0, b, error=1E-9, convergence=False):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# vérfions que ca marche\n",
    "\n",
    "# combien faut-il d'itérations\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_error(N, method=gradient_conjugué):\n",
    "    A, b = make_system(N)\n",
    "    x = method(A, np.zeros(N), b)\n",
    "    err = A @ x - b\n",
    "    return np.sqrt(err @ err)\n",
    "\n",
    "compute_error(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparons avec le gradient simple"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradient(A, x0, b, e = 1E-9, convergence=False, max_iterations=1000):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# vérifions que ca marche\n",
    "\n",
    "compute_error(10, method=gradient)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perfs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# comparons les performances\n",
    "\n",
    "seed = 123\n",
    "np.random.seed(seed)\n",
    "\n",
    "%timeit compute_error(1000, method=gradient)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed = 123\n",
    "np.random.seed(seed)\n",
    "\n",
    "%timeit compute_error(1000, method=gradient_conjugué)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Nombre d'iteration dans les 2 cas\n",
    "\n",
    "On teste avec N=1000 puis N=10000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1000\n",
    "A,b = make_system(N)\n",
    "x0 = np.zeros(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Pour le gradient simple"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "err = gradient(A, x0, b, convergence=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Pour le gradient conjugué"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "err = gradient_conjugué(A, x0, b, convergence=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Un cas réel\n",
    "\n",
    "Logiquement vous devriez être décu aussi on va tester avec un problème réel qui correspond à cet exemple\n",
    "de l'[équation de Poisson](https://doc.freefem.org/tutorials/poisson.html). Le système matriciel de ce problème est téléchargeable \n",
    "[ici](https://www.lrde.epita.fr/~ricou/cama/data/Ab.npz). Une fois le fichier sauvé, pour récupérer\n",
    "A et b faites :\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npz = np.load('/tmp/Ab.npz',allow_pickle=True)\n",
    "A = npz['A']\n",
    "b = npz['b']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Faites une étude rapide de A, indiquez quel pourcentage des valeurs de A est différent de 0. Afficher l'image de la matrice avec `plt.imshow(A)` (faire une grande image pour voir quelque chose).\n",
    "* Refaites la comparaison entre les deux méthodes avec ce système matriciel.\n",
    "* Regardez la documentation de `lin.solve` (en particulier les options) et comparer `lin.solve` à vos deux algorithmes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparaison gradient simple et conjugué"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time err = gradient_conjugué(A, np.zeros(len(A)), b, convergence=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# le gradient simple\n",
    "\n",
    "%time err = gradient(A, np.zeros(len(A)), b, convergence=True, max_iterations=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On voit la supériorité du gradient conjugué tant en nombre d'itérations (175 contre 7800) qu'en temps de calcul (environ 40 fois plus rapide)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparaison avec `lin.solve` de Scipy\n",
    "\n",
    "Le solveur de Scipy a des option que celui de Numpy n'a pas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# regarder la doc pour avoir les options optimales\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculer le résidu\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On note aussi `lin.solve` est plus rapide et sa solution est nettement meilleure... `lin.solve` utilise une méthode directe ici. Cela est dû au fait que Scipy utilise la bibliothèque Lapack (qui est imbatable)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Le gradient conjugué de Scipy (avec Lapack)\n",
    "\n",
    "Le gradient conjugué à tout son sens pour les matrices creuses aussi il est dans la partie \"sparse\" de Scipy.\n",
    "On a vu que notre matrice à plus de 99 % de valeur nulles ce qui en fait bien une matrice creuse. Aussi je\n",
    "la charge dans le format COO (téléchargez __[Acoo.npz](https://www.lrde.epita.fr/~ricou/cama/data/Acoo.npz)__) qui ne stocke que les valeurs non nulles et"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.sparse as sparse\n",
    "from scipy.sparse.linalg import cg\n",
    "\n",
    "Ac = sparse.load_npz('/tmp/Acoo.npz')\n",
    "%time x = cg(Ac, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On gagne un ordre de grandeur."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ]
}