{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Quelques opérations sur les tableaux numpy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mélange" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "a = np.fromfunction(lambda i, j: (i+1)*10 + j, (4, 5), dtype=int)\n", "a" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "rng = np.random.default_rng()\n", "rng.shuffle(a) # Mélange le talbeau a \"en-place\" (modifie le tableau a). Seul le premier axe est mélangé (les lignes ici donc)\n", "a" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "rng.shuffle(a.ravel()) # a.ravel() retourne une vue 1D sur a. Tous les éléments sont mélangés du coup.\n", "a" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Tri" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "a = np.array([5, 3, 6, 1, 6, 7, 9, 0, 8])\n", "np.sort(a) # retourne une copie triée de a" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "print(a)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "a.sort() # Tri le tableau en-place\n", "print(a)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Opération de transposition" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "a = np.array([5, 3, 6, 1, 6, 7, 9, 0, 8, 2])\n", "print(a)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "a = a.reshape((2, 5))\n", "print(a)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "c = a.T # .T retourne une vue quand c'est possible\n", "print(c, c.shape)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "c[1, 0] = -1\n", "print(a)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "c # c est une vue de la transposée de a" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# La transposition d'un tableau 1D ne fait rien. Un tableau 1D n'est pas un vecteur ligne ou colonne.\n", "# Il faut le voir simplement comme un tableau 1D, il n'y a donc pas d'axe à inverser comme pour un tableau multi-D\n", "a = np.array([5, 3, 6, 1, 6, 7, 9, 0, 8, 2])\n", "print(a)\n", "print(a.T)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Opérations sur un axe donné d'un tableau numpy\n", "\n", "La plupart des fonctions numpy ont un argument `axis` qui permet d'effectuer l'opération non pas sur tous les éléments du tableau mais sur un axe du tableau en particulier." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "a = np.arange(20).reshape(5, 4)\n", "rng.shuffle(a.ravel())\n", "print(a)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "s = a.sum(axis=0) # Somme sur les lignes\n", "print(s)\n", "print(s.shape)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "s = a.sum(axis=1) # Somme sur les colonnes\n", "print(s)\n", "print(s.shape)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "m = a.mean(axis=0) # Moyenne sur les lignes\n", "print(m)\n", "print(m.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Changer la taille d'un tableau\n", "\n", "`numpy.reshape` et `ndarray.reshape` change les tailles d'un tableau. Attention, le nombre total d'élément doit rester le même !\n", "Cette fonctoin retourne une vue lorsque c'est possible." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A = rng.integers(-5, 6, (4, 4))\n", "print(A)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(np.reshape(A, (2, 8)))\n", "print(np.reshape(A, (2, -1))) # -1 permet de laisser numpy déterminer le nombre d'éléments à mettre" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Le nombre d'élément total doit rester le même mais il est possible de modifier le nombre de dimensions\n", "print(np.reshape(A, (2, 2, 4))) # tableau 3D\n", "print()\n", "print(np.reshape(A, -1)) # tableau 1D" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "A.reshape(2, 8) # version ndarray.reshape" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# A n'a pas été changé\n", "print(A.shape)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Règles de broadcasting\n", "\n", "Le terme \"broadcasting\" décrit comment numpy gère les opérations arithmétiques entre deux tableaux de tailles différentes. Ces opérations sont faîtes élément par élément, il faudrait donc que les deux tableaux aient la même \"shape\" pour que cela marche. Le broadcasting permet à numpy de faire ces opérations lorsque les shape vérifient certaines règles sans être nécessairement égales.\n", "\n", "Deux dimensions sont compatibles si :\n", "- elles ont le même nombre d'élément ;\n", "- une des deux vaut $1$ et dans ce cas le reste du tableau est dupliqué dans cette dimension pour obtenir le même nombre d'élément que l'autre.\n", "\n", "Lors de l'opération entre deux tableaux ont fait donc, lorsque toutes les dimensions sont compatibles, l'opération sur deux tableaux dont les tailles sont le maximum des tailles des deux tableaux.\n", "\n", "**Example**: on considère deux tableaux avec les shape suivantes :\n", "```\n", "A (4d array): 8 x 1 x 6 x 1\n", "\n", "B (3d array): 7 x 1 x 5\n", "\n", "expanded B (4d): 1 x 7 x 1 x 5\n", "\n", "Result (4d array): 8 x 7 x 6 x 5\n", "```\n", "\n", "Les deux tableaux n'ont pas les mêmes shape, on ne devrait donc pas pouvoir faire d'opérations élément par élément. Mais les règles de brocasting s'appliquent :\n", "\n", "1. On compare la première dimension, celle de droite dans le résultat de `A.shape` et `B.shape`. On a $1$ pour `A` et $5$ pour `B`. Le nombre d'élément est différent mais l'un d'eux vaut 1. Le tableau `A` est donc dupliqué dans cette dimension pour obtenir $5$ éléments. Techniquement, on aurait donc `A[:, :, :, 0] = A[:, :, :, 1] = ... = A[:, :, :, 4]`.\n", "2. On passe à la dimension suivante. On a $6$ pour `A` et $1$ pour `B`. Encore une fois, le nombre d'élément n'est pas égal mais on peut \"dupliquer\" `B` dans cette dimension pour obtenir $6$ élément, comme `A`.\n", "3. Pareil pour la troisième dimension.\n", "4. `A` possède une quatrième dimension mais pas `B`. Dans ce cas, numpy ajoute une dimension à `B` avec un élément. Les règles de broadcasting s'appliquent alors comme précédemment.\n", "5. Le résultat de l'opération (+, -, *, /, par exemple) est donc un tableau de taille `8 x 7 x 6 x 5`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = np.arange(4).reshape((4, 1)) # shape 4, 1\n", "print(x.shape)\n", "print(x)\n", "print()\n", "\n", "y = np.ones(5) # shape 5,\n", "print(y.shape)\n", "print(y)\n", "print()\n", "\n", "\n", "z = x + y # ici, numpy fait du broadcasting car les tailles de ne sont pas les mêmes\n", "print(z.shape) # numpy ajoute une dimension à y pour avoir 1, 5 puis les règles de broadcasting s'appliquent\n", "print(z) # On obtient un tableau de taille 4, 5" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A = rng.random((4, 4))\n", "x = rng.random(4)\n", "\n", "print(A)\n", "print(x)\n", "\n", "u = A * x # multiplication terme à terme (broadcasting)\n", "v = A @ x # produit matrice / vecteur\n", "\n", "print(\"Produit A * x (élément par élément)\")\n", "print(u.shape)\n", "print(u)\n", "\n", "print(\"Produit matrice / vecteur\")\n", "print(v.shape)\n", "print(v)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Le broadcasting ne duplique pas vraiment les tableaux dans les dimensions qui ne sont pas égale.\n", "# numpy les voit comme ça mais ce n'est pas fait comme ça.\n", "\n", "nx = 5000\n", "x = np.arange(nx).reshape((nx, 1))\n", "\n", "ny = 10000\n", "y = np.ones(ny)\n", "\n", "%timeit np.tile(x, (1, ny)) + np.tile(y, (nx, 1)) # l'équivalent du broadcasting en dupliquant réellement les tableaux\n", "%timeit x + y # avec du broadcasting\n", "\n", "print(np.all(np.tile(x, (1, ny)) + np.tile(y, (nx, 1)) == x + y)) # On vérifie que les résultats sont les mêmes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Ajout de dimensions\n", "\n", "Vu la façon dont le broadcasting fonctionne, il pourrait être intéressant de pouvoir ajouter des dimensions de taille $1$ comme on veut.\n", "C'est possible, voici des exemples." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "B = A[:, :, np.newaxis] # np.newaxis ajoute une dimension de taille 1\n", "print(f\"A = \\n{A} with size {A.shape}\\n\")\n", "print(f\"B = \\n{B} with size {B.shape}\\n\")\n", "\n", "C = B[:, :, 0]\n", "print(f\"C = \\n{C} with size {C.shape}\\n\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# B et C sont des vues de A\n", "B[1, 2, 0] = 42\n", "print(A)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "B = A[:, np.newaxis, :]\n", "print(f\"A = \\n{A} with size {A.shape}\\n\")\n", "print(f\"B = \\n{B} with size {B.shape}\\n\")\n", "\n", "C = B[:, 0, :]\n", "print(f\"C = \\n{C} with size {C.shape}\\n\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Il est aussi possible d'utiliser `None` au lieu de `np.newaxis`\n", "B = A[:, None, :, None]\n", "print(f\"A = \\n{A} with size {A.shape}\\n\")\n", "print(f\"B = \\n{B} with size {B.shape}\\n\")\n", "\n", "C = B[:, 0, :, 0]\n", "print(f\"C = \\n{C} with size {C.shape}\\n\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "u = np.array([1, 2, 3])\n", "v = np.array([4, 5, 6])\n", "u[None, :] + v[:, None] # (1, 3) + (3, 1) -> (3, 3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercices\n", "\n", "## Matrice de distance\n", "\n", "1. Générer un tableau 1D `x` contenant les entiers de $0$ à $4$.\n", "2. Calculer la matrice `D` telle que $D_{i, j} = |x_i - x_j|$. Essayer de le faire avec du broadcasting.\n", "3. Générer un tableau 2D `X` contenant des valeurs aléatoires dans l'intervalle $[-5, 5]$ (`rng.uniform`) et de taille $(5, 2)$. Ce tableau contient les coordonnées de $5$ points en 2D.\n", "4. Calculer la matrice `D` telle que $D_{i, j} = \\|x_i - x_j\\|$ (voir `np.linalg.norm` pour la norme)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load solutions/numpy/distance-matrix.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Estimation d'aire par méthode de Monte-Carlo\n", "\n", "On cherche à estimer l'aire d'un forme compliquée. L'idée de la méthode de Monte-Carlo est de plonger cette forme complexe dans un forme plus simple dont on sait calculer l'aire (un carré par exemple). On tire ensuite des points aléatoirement dans la forme simple. On compte le nombre de point qui tombe à l'intérieur de la forme complexe (il faut pouvoir le faire). L'aire de la forme complexe peut être approché par le rapport entre le nombre de points à l'intérieur de la forme complexe et le nombre de point tiré, le tout multiplié par l'aire de la forme simple.\n", "\n", "On peut par exemple estimer l'aire du disque unité en 2D (de rayon $1$ et centré en $(0, 0)$) en le plongeant dans le carré $[-1, 1] \\times [-1, 1]$.\n", "\n", "1. Écrire une fonction `disk_area(n)` qui approche l'aire du disque en tirant `n` points aléatoirement dans le carré (utiliser `rng.uniform` et `np.count_nonzero`.\n", "2. Essayer votre fonction avec $10000$ points et comparer le résultat à la valeur exacte.\n", "3. Ajouter un argument `plot` à votre fonction `disk_area` avec la valeur par défaut `False`. Si `plot` est `True`, afficher le carré, le disque, et tous les points dans le disque. Le cercle unité peut se représenter avec la courbe $(\\cos(t), \\sin(t))$ pour $t \\in [0, 2 \\pi[$.\n", "4. Essayer avec $n = 500$ points." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load solutions/numpy/area_disk.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On cherche à connaître la vitesse de convergence de cette méthode. Au lieu de retourner l'estimation de l'aire on veut retourner un tableau `areas` de taille `n` tel que `areas[i]` est l'estimation de l'aire après avoir tiré les `i` premiers points aléatoirement.\n", "\n", "1. Modifier la fonction `disk_area` pour retourner le tableau `areas` au lieu de l'aire (utiliser `np.cumsum`).\n", "2. Appeler cette fonction avec `n` egal à $10000000$ et mettre le résultat dans un tableau `areas`.\n", "3. Calculer l'erreur entre la valeur exacte et tous les éléments de `areas`. On appelera ce tableau `err`.\n", "4. Afficher sur une même figure l'évolution du logarithme de l'erreur en fonction du logarithme du nombre de points aléatoires tirés et une ligne de pente $-0.5$. Que peut-on dire du résultat ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load solutions/numpy/area_disk_conv.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Répéter les premières questions avec le trifolium. Le trifolium est défini par l'équation paramétrique $t \\mapsto (\\cos(3t) \\cos(t), \\cos(3t) \\sin(t))$ pour $t \\in [0, \\pi]$.\n", "Un point $(x, y)$ est à l'intérieur du trifolium si $(x^2 + y^2)^2 - x(x^2 -3y^2) < 0$.\n", "Comme pour le disque unité, le trifolium est aussi inclus dans le carré $[-1, 1] \\times [-1, 1]$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load solutions/numpy/area_trifolium.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Interpolation de Lagrange\n", "\n", "Étant donné $N$ points $(x_i, y_i), 0 \\leq i \\leq N$, on veut trouver le polynôme $P$ de degré $N$ tel que $P(x_i) = y_i$ pour tout $0 \\leq i \\leq N$. Il existe un unique polynôme de degré $N$ avec cette propriété ($N+1$ équations, $N+1$ inconnues).\n", "\n", "Si on écrit $P(X) = \\sum_{i = 0}^N a_i X^i$, il faut résoudre le système linéaire suivant pour trouver les coefficients $a_i$ : $\\sum_{i=0}^N a_i x_j^i = y_j$, pour $0 \\leq j \\leq N$.\n", "\n", "Au lieu de procédé de cette façon, on va construire le polynôme $P$ comme une somme de $N+1$ polynômes $P_i$ tels que $P_i(x_i) = 1$ et $P_i(x_j) = 0$ pour tout $j \\neq i$.\n", "On peut écrire $P_i(X) = \\prod_{j \\neq i} \\frac{X - x_j}{x_i - x_j}$ et $P$ est alors défini par $P(X) = \\sum_{i=0}^N y_i P_i(X)$.\n", "\n", "On veut suivre cet algorithme pour interpoler la fonction $f(x) = -x \\sin(\\frac{7}{2} \\pi \\sqrt{x})$.\n", "\n", "Avant de commencer, voici quelques commandes pour manipuler les polynômes avec numpy." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "P = np.polynomial.Polynomial([1, 2, 3]) # Création d'un polynôme à partir de ses coefficient\n", "print(P)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Q = np.polynomial.Polynomial.fromroots([0, 1, 2]) # Création d'un polynôme à partir de ses racines\n", "print(Q)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(P+Q)\n", "print(P*Q)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(P.deriv())\n", "print(P.integ())" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "X = np.polynomial.Polynomial.identity() # Pour récupérer le polynôme identité\n", "R = 3 * X + 5 * X**2 - 1 # à partir duquel on peut construire d'autres polynômes\n", "print(R)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# On peut évaluer un polynôme simplement\n", "print(R(1.0))\n", "print(R(np.arange(5))) # élément par élément si on passe un tableau" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Faire un \"fit\" par un polynôme\n", "x = np.linspace(-2, 2, 100)\n", "y = 1 / (1 + x**2)\n", "P = np.polynomial.Polynomial.fit(x, y, deg=5) # On cherche le polynôme de degré 5 qui passe au plus proche des points (x, y)\n", "\n", "plt.figure()\n", "plt.plot(x, y, label=\"Fonction 1 / (1 + x^2)\")\n", "plt.plot(x, P(x), label=\"Fit\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Revenons à l'interpolation de Lagrange\n", "\n", "1. Définir la fonction $f$\n", "3. Définir une fonction `lagrange(x, y)` qui calcule et retourne le polynôme $P$ à partir des points $(x_i, y_i)$ contenus dans les deux tableaux en argument `x` et `y`.\n", "4. Définir un tableau `x` contenant $6$ points régulièrement espacés dans $[0, 1]$ (`np.linspace`).\n", "5. Interpoler la fonction $f$ aux points `x` (`y` est le tableau avec les valeurs $f(x_i)$).\n", "6. Afficher la fonction $f$, le polynôme $P$, et les points d'interpolation $(x_i, f(x_i))$ dans $[0, 1]$. Ajouter une légende et un titre à la figure.\n", "7. Vérifier que le polynôme $P$ a les mêmes valeurs que $f$ aux points d'interpolations\n", "\n", "Le module [scipy.interpolate](https://docs.scipy.org/doc/scipy/reference/interpolate.html) contient beaucoup de fonction pour l'interpolation." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load solutions/numpy/lagrange.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Théorème d'approximation de Weierstrass\n", "\n", "Il existe une preuve probabiliste du théorème d'approximation de Weierstrass par Bernstein qui utilise les polynômes de Bernstein. La base des polynômes de Bernstein de degré $n$ est définie par $b_{\\nu, n}(x) = C_n^\\nu x^\\nu (1 - x)^{n - \\nu}$, où $C_n^\\nu$ est le coefficient binomiale.\n", "\n", "Étant donné une fonction continue $f : [0, 1] \\longrightarrow \\mathbb{R}$, le polynôme de Bernstein $B_{n, f}(x) = \\sum_{\\nu = 0}^n f\\left(\\frac{\\nu}{n}\\right) b_{\\nu, n}(x)$ converge uniformément vers $f$ quand $n$ tend vers l'infini.\n", "\n", "Si $X$ est une variable aléatoire suivant une loi binomiale de paramètres $n \\in \\mathbb{N}$ et $x \\in [0, 1]$, alors $P(X = \\nu) = b_{\\nu, n}(x)$. Ainsi, $B_{n, f}(x)$ est la valeur moyenne de $f(X / n)$.\n", "\n", "1. Définir une fonction $f$ qui retourne $f(x) = -x \\sin(\\frac{7}{2} \\pi \\sqrt{x})$.\n", "2. Définir une fonction `bernstein(x, f, n, N)`, qui retourne le tableau contenant les valeurs $B_{n, f}(x_i), x_i \\in x$ (voir `rng.binomial` et `np.mean`) et où :\n", " - `x` est un tableau de points dans $[0, 1]$\n", " - `f` est la fonction à approcher\n", " - `n` est le degré du polynôme de Berstein\n", " - `N` est le nombre de nombres aléatoires à tirer suivant la loi binomiale de paramètre $n$ et $x_i$ pour chaque élément $x_i$ de `x`.\n", "4. Définir le tableau `x` contenant des valeurs régulièrement espacées dans $[0, 1]$.\n", "5. Appeler votre fonction `bernstein` avec ce tableau `x`, la fonction $f$, `n = 50`, et `N = 100000`.\n", "6. Afficher la fonction $f$ et son approximation en utilisant le tableau `x` défini plus haut et ajouter un légende.\n", "7. Pour $n$ dans $\\{10, 20, 40, 80, 160\\}$, afficher l'approximation de $f$ par le polynôme de Bernstein correspondant en bleu et changer l'opacité des courbes de telle manière que plus $n$ est petit, plus la courbe est transparente (utiliser le paramètre `alpha` de `plt.plot`). Ajouter la courbe de $f$ en rouge." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load solutions/numpy/bernstein.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Méthode d'Euler explicite\n", "\n", "On veut résoudre numériquement l'équation différentielle $y'(t) = -2 y(t) + 2 - e^{-4t}$ avec $y(0) = 1$.\n", "Cette équation peut se réécrire sous la forme plus générale $y'(t) = f(t, y(t))$ où dans notre cas $f(t, y) = -2 y + 2 - e^{-4t}$.\n", "\n", "La solution exacte de cette équation différentielle est $y_{ex}(t) = \\frac{1}{2} (2 - e^{-2t} + e^{-4t})$.\n", "\n", "On peut résoudre cette équation numériquement en utilisant la méthode d'Euler explicite. L'algorithme est le suivant :\n", "\n", "- Étant donné un pas d'espace $h > 0$ et un nombre d'itération $N$, on pose $y_0 = y(0)$, $t_0 = 0$ et $k = 0$.\n", "- Tant que $k < N$: \n", " - calculer $y_{k+1}$ comme étant la valeur obtenue en partant de $y_k$ et en avançant d'un pas $h$ dans la direction $f(t_k, y_k)$. On a donc $y_{k+1} = y_k + h f(t_k, y_k)$.\n", " - poser $t_{k+1} = t_k + h$.\n", "- Retourner les valeurs $y_k$ dans un tableau qui sont les approximations des $y_{ex}(t_k)$.\n", "\n", "1. Définir une fonction `sol(t)` qui retourne la solution exacte au temps `t`.\n", "2. Définir la fonction `f(t, y)` de l'énoncé. En particulier, cette fonction retourne $y_{ex}'(t)$ si on prend `y` $= y_{ex}(t)$.\n", "3. Écrire une fonction `forward_euler(f, y0, h, N)` qui retourne un tableau `t` et un tableau `y` contenant respectivement les temps $t_k$ et les approximations $y_k$.\n", "4. Appeler la fonction `forward_euler` avec la fonction `f` définie précédemment, $y_0 = 1$, $h = 0.05$, et $N = 40$.\n", "5. Tracer la solution exacte aux temps `t` et son approximation par la méthode d'Euler explicite et ajouter une légende." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load solutions/numpy/forward_euler.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Méthodes de gradient\n", "\n", "La méthode de gradient à pas fixe est une méthode simple utilisée pour trouver des minimum de fonctions lorsqu'on sait calculer leur dérivée.\n", "Étant donné un pas $h$ et un point de départ $x_0$, l'idée est de suivre la direction de plus forte pente sur une distance $h$. La direction de plus forte pente est donnée par l'opposée de la dérivée de la fonction. On construit donc une suite de point $x_n$ qui vérifient $x_{n+1} = x_n - h f'(x_n)$.\n", "\n", "On veut appliquer cette méthode à la fonction $f$ définie par $f(x) = -x \\sin(\\frac{7}{2} \\pi \\sqrt{x})$.\n", "\n", "1. Écrire une fonction `f(x)` qui retourne la valeur de $f(x)$.\n", "2. Afficher la courbe de $f$ pour $x \\in [0, 1]$.\n", "3. Écrire une fonction `df(x)` qui retourne la valeur de $f'(x)$.\n", "4. Écrire une fonction `steepest_descent(df, x0, h, eps)` qui implémente l'algorithme de gradient à pas fixe. L'algorithme s'arrête lorsque la distance entre deux itérations consécutives est plus petite que la tolérance `eps` donné en argument.\n", "5. Appeler la fonction `steepest_descent` avec `x0 = 0.9`, `h = 0.01`, et `eps = 1e-10`.\n", "6. Dans le module `scipy.optimize` il existe beaucoup de fonction d'optimisation. Essayer d'utiliser la fonction `minimize_scalar` de ce module et comparer le résultat avec celui donnée par votre fonction `steepest_descent`. On pourra importer le module avec la commande `import scipy.optimize as spo` et utiliser la fonction `spo.minimize_scalar`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load solutions/numpy/steepest_descent.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On peut modifier la méthode de gradient à pas fixe pour ajouter de \"l'inertie\".\n", "\n", "En partant d'un point $x_0$ et d'une \"vitesse\" $v_0 = 0$:\n", "\n", "- mettre à jour la vitesse : $v_{k+1} = \\mu v_k - h f'(x_k)$, où $\\mu$ est un paramètre donné (qui pourrait ressembler à de la \"friction\").\n", "- mettre à jour la position : $x_{k+1} = x_k - \\mu v_k + (1 - \\mu) v_{k+1}$.\n", "\n", "L'algorithme s'arrête quand deux positions consécutives sont suffisament proche (plus proche qu'une tolérance $\\epsilon$).\n", "\n", "C'est la méthode de Nesterov. Une variante de la méthode de Nesterov appelée méthode d'Adam est souvent utilisé en machine learning.\n", "\n", "1. Implémenter la méthode de Nesterov en écrivant une fonction `nesterov(df, x0, h, mu, eps)`.\n", "2. Appeler cette méthode avec les mêmes paramètres que pour la méthode de gradient à pas fixe et différentes valeurs de $\\mu \\in [0, 1]$. En particulier, si $\\mu = 0$, on retrouve l'algorithme de gradient à pas fixe." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load solutions/numpy/nesterov.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On veut maintenant comparer les suites $x_n$ construites avec ces deux méthodes.\n", "\n", "1. Si ce n'est pas déjà le cas, modifier les fonctions `steepest_descent` et `nesterov` afin qu'elles retournent toute la suite $(x_k)$ construite au lieu du dernier point uniquement.\n", "2. Afficher les approximations avec des markers différents et ajouter une légende." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load solutions/numpy/optimisation_plot.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Régression linéaire\n", "\n", "L'objectif est d'approcher un nuage de points par une fonction affine.\n", "\n", "Soit $(x_i, y_i), 0 \\leq i \\leq N$ un nueage de points.\n", "On veut trouver deux paramètres $a$ et $b$ tels que la fonction $f(x) = ax + b$ est aussi \"proche\" que possible des points.\n", "Plus précisément, on veut minimiser la fonction $g(a, b) = \\sum_{i = 0}^N (a x_i + b - y_i)^2$.\n", "\n", "Au minimum, le gradient de $g$, $\\nabla g (a, b) = (\\frac{\\partial g}{\\partial a} (a, b), \\frac{\\partial g}{\\partial b} g(a, b)$, est nul. Ainsi, on doit résoudre un système linéaire de deux équations (les deux composantes du gradient sont égales à $0$ en même temps) avec deux inconnues ($a$ et $b$).\n", "\n", "1. Écrire une fonction affine `f(x)` (de votre choix).\n", "2. Créer un tableau `x` de $20$ points régulièrement espacé dans $[0, 1]$.\n", "3. Créer un tableau `y` les valeurs de `f(x)` et ajouter y une petite perturbation (utiliser des nombres aléatoires suivant une loi normal centrée réduite par exemple : `rng.normal`).\n", "4. Afficher les points $(x_i, y_i)$ et la droite $f$. La droite devrait être proche des points mais sans passer exactement dessus à cause de la perturbation ajoutée.\n", "5. Question de maths, sans python : calculer le gradient de $g$ et écrire l'équation $g(a, b) = (0, 0)$ sous la forme d'un système linéaire $A (a, b) = (c1, c2)$ où $A$ est une matrice $2 \\times 2$ et $(c1, c2)$ un vecteur qui ne dépend ni de $a$ ni de $b$.\n", "6. Dans une fonction `reglin(x, y)` résoudre le système linéaire et retourner les deux coefficients $a$ et $b$ trouvés. Vous pouvez ou bien résoudre le système linéaire vous même ou bien utiliser la fonction `np.linalg.solve` par exemple.\n", "7. Appeler la fonction `reglin` avec les tableaux `x` et `y` créés au tout début.\n", "8. Quelles sont les valeurs de $a$ et $b$ obtenue ? Sont-elles proches de la droite affine $f$ ?\n", "9. Afficher la fonction $f$, les points $(x_i, y_i)$, et la droite affine obtenue avec les coefficients $a$ et $b$ trouvés.\n", "10. Essayer d'utiliser la fonction `np.polynomial.Polynomial.fit` pour faire cette régression linéaire et comparer le résultat." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load solutions/numpy/linear_regression.py" ] }, { "cell_type": "markdown", "metadata": { "scrolled": true }, "source": [ "Parce qu'on minimise une fonction, on peut réutiliser la méthode de Nesterov de l'exercice précédent par exemple. Pour cela, il faut créer une fonction qui prend un seul argument sous forme de tuple $(a, b)$ et qui retourne le gradient de $g$ au point $(a, b)$.\n", "\n", "Le problème est que le gradient de $g$ dépend aussi des points $(x_i, y_i)$. Une manière de résoudre ce problème est de créer une fonction qui retourne une fonction.\n", "\n", "- On va créé une fonction `make_dg(x, y)` qui définie une autre fonction (en utilisant `def`) `dg(ab)` d'un seul argument `ab` qui désignera le tuple `(a, b)`.\n", "- La fonction `dg(ab)` calcule et retourne le gradient de $g$ au point $(a, b)$ en utilisant le `x` et le `y` en argument de `make_dg`.\n", "- La fonction `make_dg` retourne la fonction `dg`.\n", "\n", "1. Écrire la fonction `make_dg` décrite au dessus.\n", "2. L'appeler avec les tableaux `x` et `y`. Elle retourne la fonction `dg` à un seul argument.\n", "3. Appeler la fonction `nesterov` avec la fonction `dg` en argument.\n", "4. Vérifier le résultat." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load solutions/numpy/linear_regression_nesterov.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Affichons la fonction $g$ ! C'est une fonction de deux variables à valeurs réelles, sa représentation est une surface.\n", "\n", "1. Écrire une fonction `make_g(x, y)` qui définie une fonction `g(a, b)` et la retourne. La fonction `g(a, b)` calcule le tableau contenant les valeurs $g(a_i, b_i)$ où $a_i$ sont les éléments de `a` et $b_i$ les éléments de `b`.\n", "2. Appeler la fonction `make_g` avec les points `x` et `y` pour construire la fonction `g`.\n", "3. Créer un tableau `a` contenant des points régulièrement espacés et centrés autours de la valeur du coefficient directeur de la droite affine $f$. Faire la même chose pour `b`.\n", "4. La fonction `np.meshgrid` prend deux tableaux 1D en argument (`u` et `v` par exemple) et retourne deux tableaux 2D (`U` et `V`) tels que `U[i, j] = u[j]` and `V[i, j] = v[i]`. Utiliser cette fonction avec les tableaux `a` et `b` afin d'obtenir deux tableaux `A` et `B` de même taille avec les valeurs de `a` sur les lignes de `A` et les valeurs de `b` sur les colonnes de `B`.\n", "5. Appliquer la fonction `g` aux tableaux `A` et `B` et mettre le résultat dans une variable `Z`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load solutions/numpy/linear_regression_plot.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`plt.contour` et `plt.contourf` affiche les lignes de niveaux d'une fonction. Elles prennent au moins trois arguments :\n", "\n", "- un tableau 2D X contenant les abscisses des points d'évaluation.\n", "- un tableau 2D Y contenant les ordonnées a 2D array Y containing the coordinates on the axe y.\n", "- un tableau 2D Z contenant les valeurs de la fonction dont on veut les contours aux points (X, Y)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure()\n", "plt.contour(A, B, Z, levels=40) # 40 contours\n", "plt.colorbar()\n", "\n", "plt.figure()\n", "plt.contourf(A, B, Z, levels=40) # 40 contours\n", "plt.colorbar()\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Un bout de code pour afficher des surfaces en 3D" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from mpl_toolkits.mplot3d import Axes3D\n", "from matplotlib import cm\n", "\n", "fig = plt.figure(figsize=(10,10))\n", "ax = fig.add_subplot(111, projection='3d')\n", "ax.plot_surface(A, B, Z, cmap=cm.coolwarm)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.13" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": { "height": "calc(100% - 180px)", "left": "10px", "top": "150px", "width": "336px" }, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }