{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "__Programmation vectorielle__\n",
    "\n",
    "Le but des exercices est\n",
    "\n",
    "* d'avoir un programme qui donne la bonne réponse\n",
    "* qui soit le plus rapide possible (et pour cela on utilise massivement Numpy)\n",
    "\n",
    "En règle général si vous avez des `for` imbriqués c'est mauvais signe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "np.set_printoptions(precision=10, linewidth=150, suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Méthode du pivot de Gauss partiel\n",
    "\n",
    "L'ennoncé est dans le cours. On verifiera sur le cas qui génère des\n",
    "erreurs d'arrondis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x mieux = [1.0000019  0.99999905]\n",
      "[1. 3.]\n"
     ]
    }
   ],
   "source": [
    "def solve_partial_gauss(A, b):\n",
    "    A = A.astype(np.float64)   # si A a des entiers, on va avoir des erreurs de calculs\n",
    "    for i in range(len(A)-1):\n",
    "        pivot = i + np.argmax(np.abs(A[i:,i]))\n",
    "        A[[i,pivot]] = A[[pivot, i]]\n",
    "        b[[i,pivot]] = b[[pivot, i]]\n",
    "        coefs = - A[i+1:,i] / A[i,i]  # Normalement il faut tester que A[i,i] != 0\n",
    "        A[i+1:, i:] += np.outer(coefs, A[i, i:]) # ou np.einsum('i,j -> ij', coefs, A[i, i:])\n",
    "        b[i+1:] += coefs * b[i]\n",
    "    # A est maintenant triangulaire surpérieur\n",
    "    res = np.zeros(len(b), dtype=b.dtype)\n",
    "    res[-1] = b[-1] / A[-1,-1]\n",
    "    for i in range(len(A)-1)[::-1]:\n",
    "        res[i] = (b[i] - A[i,i+1:] @ res[i+1:]) / A[i,i]\n",
    "    return res\n",
    "e = 1E-6\n",
    "A = np.array([[e, 1], [1, 2]], dtype='float32')\n",
    "b = np.array([1., 3.], dtype='float32')\n",
    "x = solve_partial_gauss(A,b)\n",
    "print(\"x mieux =\", x)\n",
    "print(A @ x )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Factorisation de Choleski\n",
    "\n",
    "Il s'agit de décomposer A en $A = B\\, B^T$ avec B une matrice triangulaire inférieure. Cela n'est possible\n",
    "que si la matrice A est symétrique et définie positive (c'est d'ailleurs un facon de vérifier qu'une\n",
    "matrice est définie positive).\n",
    "\n",
    "Écrire l'algorithme de Choleski qui prend A et retourne B (pour deviner l'algorithme, essayez de trouver les \n",
    "coefficients de B à partir des coefficients de A sur une matrice A 4x4)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Rappel : pas de boucles `for` imbriquées mais des opérations vectorielles et matricielles (sur des sous-matrices).\n",
    "\n",
    "Créer une matrice symétrique définie positive et vérifier que sonprogramme marche bien."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Amméliorer Jacobi\n",
    "\n",
    "Lorsqu'on écrit une itération de la méthode de Jacobi avec l'ensemble des coefficients, on constate que\n",
    "si on calcule la nouvelle valeur de **x** élément par élement alors lorsqu'on veut mettre à jour `x[1]`, \n",
    "on connait déjà `x[0]`. Idem lorsqu'on met à jour `x[2]` on connait déjà `x[0]` et `x[1]`, etc.\n",
    "\n",
    "L'idée de la version amméliorée de Jacobi est d'utiliser la nouvelle valeur de `x[0]` et non pas l'ancienne\n",
    "comme c'est le cas dans l'algorithme du cours. Ainsi en utilisant des valeurs mise à jour on peut espérer\n",
    "converger plus vite.\n",
    "\n",
    "Dans ce cas, chaque itération demande un calcul ligne par ligne et donc une boucle `for` dans la boucle `for` sur\n",
    "les itérations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "#### Test d'arrêt\n",
    "\n",
    "On ajoutera un argument `error` à la fonction qui indique la précision désirée du résultat. Malheureusement pour connaitre la précision de notre calcul il faut connaitre la solution. On pourrait aussi calculer \n",
    "${\\bf A \\, x}^t - {\\bf b}$ mais ce n'est pas non plus l'écart entre ${\\bf x}^t$ et la solution (en plus\n",
    "c'est un calcul en n² donc lourd).\n",
    "\n",
    "Aussi on va regarder quand la\n",
    "convergence ralentit et donc on s'arrête lorsque\n",
    "$||{\\bf x}^{t+1} - {\\bf x}^t|| < \\texttt{error}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "#### Bonus\n",
    "\n",
    "Réécrire la méthode sans la boucle `for` mais en prenant bien le nouvelles valeurs de `x`. Pour cela on va découper la matrice **A** en deux matrices triangulaires."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# les méthodes de Numpy triu et tril seront utiles\n"
   ]
  },
  {
   "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.13.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}