{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Numpy : manipulation des tableaux\n",
    "\n",
    "Un des principaux apports du module `numpy` (numerical python) est son type de tableaux. Ce sont des conteneurs comme les listes Python mais ils sont fait pour faire du calcul.\n",
    "\n",
    "## Pourquoi des nouveaux tableaux à la place des listes Python ?\n",
    "\n",
    "- La somme de deux listes concataine les deux listes, ce n'est pas la somme de deux vecteurs.\n",
    "- Les listes Python sont lentes et demande plus de mémoire.\n",
    "- Pour faire des matrices ou des tableaux à plusieurs dimensions, il faut faire des listes de listes de listes... Ce qui devient compliqué à utiliser et à débuguer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np   # on importe le module `numpy` que l'on renomme `np` (plus rapide)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l1 = [1, 2, 3]\n",
    "l2 = [4, 5, 6]\n",
    "l1 + l2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a1 = np.array(l1)  # transformation de la liste l1 en tableau numpy\n",
    "a2 = np.array(l2)  # pareil avec l2\n",
    "a1 + a2            # somme terme à terme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from random import random      # génération de nombre aléatoire\n",
    "from operator import truediv   # fonction division"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "N = 1_000_000\n",
    "l1 = [random() for i in range(N)]\n",
    "l2 = [random() for i in range(N)]\n",
    "%timeit s = sum(map(truediv, l1, l2))   # on mesure le temps pour faire la division terme à terme de deux liste en Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "a1 = np.array(l1)\n",
    "a2 = np.array(l2)\n",
    "%timeit s = np.sum(a1/a2)   # La même chose avec numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Les tableaux `numpy`\n",
    "\n",
    "- Il y a des différences importantes entre les tableau numpy et les listes python :\n",
    "    - les tableaux numpy ont une taille fixée à la création\n",
    "    - les éléments des tableaux numpy doivent tous être du même type\n",
    "- La plupart des logiciel scientifique en python utilisent numpy de nos jours.\n",
    "- numpy utilise du code compilé (meilleures performances)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "a = np.array([0, 1, 2, 3])  # création d'un array (tableau numpy) depuis une liste\n",
    "b = np.array((4, 5, 6, 7))  # création d'un array depuis un tuple"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(a)\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Les opérations sont effectuées élément par élément"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(a * b)\n",
    "print(a + b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(5 * a)\n",
    "print(5 + a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(a @ b)   # @ est la multiplication de matrices. \n",
    "               # Si le premier array est 1D, il est vu comme une matrice à une ligne.\n",
    "               # Si le deuxième array est 1D, il est vu comme une matrice à une colonne.\n",
    "               # Ici, cela correspond donc au produit scalaire.\n",
    "print(np.dot(a, b))   # La même chose que @."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Propriétés des array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "a = np.array([[1, 2, 3], [4, 5, 6]]) # création d'un tableau 2D à partir d'une liste de liste\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "type(a) # Le type de a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "a.dtype # Le type des éléments de a (tous les mêmes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "np.shape(a)  # Retourne le nombre d'éléments dans chaque dimension\n",
    "a.shape      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "np.size(a)   # Retourne le nombre total d'éléments\n",
    "a.size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "a.ndim  # nombre de dimensions de a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Important** :\n",
    "- Toujours utiliser `shape` ou `size` pour les tableaux numpy au lieu de `len`.\n",
    "- `len` donne les mêmes informations mais seulement pour un tableau 1D."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Functions pour allouer des tableaux numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x = np.zeros((3,))\n",
    "print(f\"np.zeros((3,))\\n{x}, {x.dtype}\\n\")\n",
    "\n",
    "x = np.zeros((3,), dtype=np.int64)\n",
    "print(f\"np.zeros((3,), dtype=np.int64)\\n{x}, {x.dtype}\\n\")\n",
    "\n",
    "x = np.ones((3, 5))\n",
    "print(f\"np.ones((3, 5))\\n{x}\\n\")\n",
    "\n",
    "x = np.eye(4)\n",
    "print(f\"np.eye(4)\\n{x}\\n\")\n",
    "\n",
    "y = np.zeros_like(x)\n",
    "print(f\"np.zeros_like(x)\\n{y}\\n\")\n",
    "\n",
    "z = np.empty((2, 3))\n",
    "print(f\"np.empty((2, 3))\\n{z}\\n\")\n",
    "\n",
    "z.fill(42)\n",
    "print(f\"z.fill(42)\\n{z}\\n\")\n",
    "\n",
    "w = np.empty((2, 3))\n",
    "print(f\"np.empty((2, 3))\\n{w}\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il est possible de créer un tableau depuis une fonction qui prend en argument des tableaux d'indices."
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "Voici deux fonctions pour créer des tableaux un peu spéciaux que l'on va beaucoup utiliser par la suite !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# `linspace` retourne des nombres équirépartis dans un intervalle.\n",
    "# Les extremitées de l'intervalle sont incluses.\n",
    "x = np.linspace(0, 1, 10)   # 10 points entre 0 et 1\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# On peut récupérer le pas entre deux points avec l'argument `retstep` mis à `True`\n",
    "x, dx = np.linspace(0, 1, 10, retstep=True)\n",
    "print(dx)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# `arange` retourne des nombres équirépartis dans un intervalle, comme `linspace`.\n",
    "# La différence est que l'on donne le pas entre deux points plutôt que le nombre de points.\n",
    "# L'extrémité gauche est incluse, celle de droite non.\n",
    "x = np.arange(0, 1, 0.1)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Depuis une **diagonale** :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.diag(range(4)))\n",
    "print()\n",
    "print(np.diag([1, 4, 2, 3], 1))\n",
    "print()\n",
    "print(np.diag(np.diag(a)))  # np.diag retourne la diagonale d'un tableau 2D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les fonctions qui génèrent des nombres aléatoirement sont dans [numpy.random](https://numpy.org/doc/stable/reference/random/index.html#module-numpy.random)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Création d'un générateur de nombre aléatoire basique\n",
    "rng = np.random.default_rng()\n",
    "\n",
    "print(dir(rng))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = rng.random((5, 5))\n",
    "print(f\"rng.random(5, 5)\\n{x}\\n\")\n",
    "\n",
    "x = rng.integers(-5, 3, 10)\n",
    "print(f\"rng.integers(-5, 3, 10)\\n{x}\\n\")\n",
    "\n",
    "x = rng.poisson(lam=10, size=(3, 5))\n",
    "print(f\"rng.poisson(lam=10, size=(3, 5))\\n{x}\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# création d'un générateur de nombres aléatoires en fixant la graine pour la reproduction des résultats\n",
    "rng = np.random.default_rng(1234)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Pour allouer un array avec la même taille et le même type d'élément qu'un autre :\n",
    "\n",
    "- `empty_like` \n",
    "- `ones_like`\n",
    "- `zeros_like`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Indexation, slicing et autres"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'indexation est similaire aux listes avec des tableaux 1D\n",
    "\n",
    "- On accède aux élément d'un array avec l'opérateur `[]`.\n",
    "- Les indices démarrent à $0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([10, 11, 12, 13, 14])\n",
    "print(a[0])\n",
    "print(a[4])\n",
    "print(a[-1])\n",
    "print(a[-5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(a[5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(a[-6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Le slicing aussi est similaire aux listes pour les tableaux 1D : x[lower:upper:step]\n",
    "\n",
    "- On extrait un sous-array en donnant un indice de départ et un indice d'arrivé.\n",
    "- L'indice de départ est inclus mais pas celui d'arrivé !\n",
    "- Le pas par défaut est 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "a = np.array([10, 11, 12, 13, 14])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(a[:2])\n",
    "print(a[-5:-3])    # indices négatifs\n",
    "print(a[0:2])\n",
    "print(a[-2:])    \n",
    "print(a[:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(a[::2])\n",
    "print(a[::-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**EXERCICES**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Approximation d'une dérivée**\n",
    "\n",
    "On calcul une approximation de la dérivée de $f(x) = \\sin(x)$ avec une méthode de différences finies.\n",
    "\n",
    "$$\n",
    "f'(x) \\approx \\frac{f(x+h)-f(x)}{h}\n",
    "$$\n",
    "\n",
    "1. Discrétiser l'intervalle $[0, 4\\pi]$ avec $100$ points (utiliser `np.linspace`). Mettre ces points dans une variable `x` et le pas d'espace dans une variable `h`.\n",
    "2. Créer le vecteur `y` contenant $y = f(x) = \\sin(x)$ (utiliser `np.sin`). On rappelle que les opérations sur les arrays sont effectuées éléments par éléments avec les fonctions `numpy`.\n",
    "3. Créer le vecteur `dfx` contenant les quantités $\\frac{f(x_{i+1}) - f(x_i)}{h}$ (utiliser des différences de sous-tableaux bien choisis)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# %load solutions/numpy/derivative.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Affichage de l'approximation**\n",
    "\n",
    "Il existe beaucoup de module python pour la visualisation de données. Le principal et le plus utilisé est `matplotlib`. C'est un très gros module qui peut faire plein de chose. Il y a un liste d'exemples dont on peut s'inspirer [ici](https://matplotlib.org/stable/gallery/)\n",
    "\n",
    "Dans l'exercice précédent, nous avons calculé une approximation de la dérivée de la fonction $f$. On pourrait afficher cette approximation et la dérivée exacte pour les comparer visuellement.\n",
    "\n",
    "Tout d'abord, il faut importer le module `matplotlib` avec la commande de la cellule suivante."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import du module `matplotlib` pour faire des représentations graphiques\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La commande pour afficher une courbe est `plt.plot`. On lui donne deux tableaux représentant les coordonnées de points en 2D (abscisses et ordonnées). La commande `plt.plot` va relier ces points entre eux par des droites pour faire la courbe correspondante.\n",
    "\n",
    "Voici des exemples !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()  # créer une nouvelle fenêtre graphique\n",
    "\n",
    "x = np.linspace(-6*np.pi, 6*np.pi, 100)     # abscisses des points. On prend 100 points dans [-6pi, 6pi]\n",
    "y = np.sin(x) / x                           # ordonnées des points.\n",
    "\n",
    "plt.plot(x, y)    # Affichage du graphique\n",
    "plt.show()        # Permet de dire que l'on veut afficher les graphiques. Cette instruction n'est pas obligatoire dans un notebook.\n",
    "                  # Sans `plt.show()`, les graphiques seront affichés à la fin de la cellule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Deux graphiques sur la même figure. Le `plt.show` est après les deux `plt.plot`.\n",
    "plt.figure()\n",
    "plt.plot(x, y)\n",
    "plt.plot(x, np.sin(2*x) / x)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# On peut ajouter un titre et une légende à la figure\n",
    "plt.figure()\n",
    "plt.plot(x, y, label=\"sin(x)/x\")                    # l'option `label` est la légende de cette courbe\n",
    "plt.plot(x, np.sin(2*x) / x, label=\"sin(2x)/x\")\n",
    "plt.legend()                                        # affichage de la légende\n",
    "plt.title(\"Un graphique avec deux fonctions et une légende\")  # titre de la figure\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Il est possible de changer le style et la couleur\n",
    "plt.figure()\n",
    "plt.plot(x, y, \"b--\", label=\"sin(x)/x\")                  # b-- = ligne discontinue bleue\n",
    "plt.plot(x, np.sin(2*x) / x, \"r-x\", label=\"sin(2x)/x\")   # r-x = ligne continue rouge et des marker en forme de croix \n",
    "plt.legend()\n",
    "plt.title(\"Un graphique avec deux fonctions et une légende\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "À votre tour maintenant ! Tracer la courbe de la dérivée exacte de $f$ et son approximation sur la même figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# %load solutions/numpy/derivative_plot.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Intégration numérique**\n",
    "\n",
    "On cherche à calculer une approximation de l'intégrale suivante :\n",
    "$$\n",
    "I = \\int_{0}^{1} x \\sin\\left(\\frac{7}{2} \\pi \\sqrt{x}\\right) dx\n",
    "$$\n",
    "\n",
    "**Méthode des rectangles**\n",
    "\n",
    "Soit $f$ la fonction définie par $f(x) = x \\sin\\left(\\frac{7}{2} \\pi \\sqrt{x}\\right)$ pour $x$ dans $[0, 1]$.\n",
    "La méthode des rectangles approche l'aire sous la courbe de $f$ par la somme de l'aire de rectangles bien définis.\n",
    "\n",
    "Etant donné $N$ points $x_i, 0 \\leq i \\leq N-1$ dans $[0, 1]$ tels que $x_i < x_{i+1}$ et $x_0 = 0$ et $x_{N-1} = 1$, le j-ième rectangle est le domaine $[x_j, x_{j+1}] \\times [0, f(x_j)]$ pour $0 \\leq j \\leq N-2$.\n",
    "\n",
    "\n",
    "1. Écrire une fonction $f$ qui prend un tableau de points $x$ en argument et retourne un tableau avec les valeurs $x \\sin\\left(\\frac{7}{2} \\pi \\sqrt{x}\\right)$.\n",
    "2. Créér un tableau `x` contenant une discrétisation de $[0, 1]$ avec $N = 100$ points.\n",
    "3. Afficher la courbe de la fonction $f$ en utilisant les points d'abscisses contenus dans `x`.\n",
    "4. Écrire une fonction `rectL` qui prend deux arguments :\n",
    "    - une fonction $f$\n",
    "    - un tableau de points $x$\n",
    "   et qui retourne l'approximation de l'intégrale de $f$ avec la méthode des rectangles.\n",
    "5. Utiliser votre fonction `rectL` pour calculer une approximation de l'intégrale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load solutions/numpy/integration_rect.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le module `scipy.integrate` contient des fonctions pour faire de l'intégration numérique. Utilisons les pour vérifier notre résultat.\n",
    "\n",
    "Regarder la documentation de la fonction `scipy.integrate.quad` pour connaître ses arguments et regarder les exemples pour voir comment calculer $I$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.integrate as spi\n",
    "spi.quad(f, 0, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On peut visualiser l'approximation en traçant la courbe de $f$ et les rectangles utilisés pour approcher l'aire sous cette courbe.\n",
    "\n",
    "La fonction `plt.fill_between` rempli une région définie par les points $(x, y1)$ et $(x, y2)$ avec la couleur de votre choix. La deuxième liste de points $(x, y2)$ n'est pas obligatoire. Si elle n'est pas spécifié, `plt.fill_between` utilise l'axe des abscisses. \n",
    "On peut l'utiliser pour remplir nos rectangles. Voici des exemples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "s = np.linspace(0, np.pi, 100)\n",
    "y1 = np.sin(s)\n",
    "y2 = np.cos(s)\n",
    "plt.fill_between(s, y1)  # on remplit le domaine entre les courbes y1 et l'axe des abscisses (pas de y2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "s = np.linspace(0, np.pi, 100)\n",
    "y1 = np.sin(s)\n",
    "y2 = np.cos(s)\n",
    "plt.fill_between(s, y1, y2)  # on remplit le domaine entre les courbes y1 et y2\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# on change la couleur, l'opacité et on ajoute les courbes y1 et y2\n",
    "s = np.linspace(0, np.pi, 100)\n",
    "y1 = np.sin(s)\n",
    "y2 = np.cos(s)\n",
    "\n",
    "plt.figure()\n",
    "plt.fill_between(s, y1, y2, color=\"g\", alpha=0.15)   # alpha = opacité.\n",
    "plt.plot(s, y1, \"r\", label=\"y1\")\n",
    "plt.plot(s, y2, \"b\", label=\"y2\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Écrire un code qui affiche les rectangles utilisés pour calculer l'approximation de l'intégrale $I$ en bleu avec une opacité de $0.2$ et la courbe de $f$ (en bleue aussi).\n",
    "\n",
    "1. Pour la visualisation, utiliser $N = 10$ points et une discrétisation de $[0, 1]$ avec $100$ points pour la courbe de $f$.\n",
    "2. Changer la valeur de $N$ pour observer la convergence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load solutions/numpy/integration_rect_plot.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Méthode des trapèzes**\n",
    "\n",
    "On a utlisé des rectangles car l'aire d'un rectangle est simple à calculer, mais ils approchent mal une courbe. Ici, on veut faire la même chose mais avec des trapèzes plutôt que des rectangles.\n",
    "\n",
    "Refaire les questions précédentes mais avec des trapèzes. La fonction s'appellera `trapz` au lieu de `rectL`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load solutions/numpy/integration_trapz.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Autres méthodes et subtilités pour extraire des sous-tableaux\n",
    "\n",
    "### Les slices sont des références\n",
    "\n",
    "- Les slices sont des références !\n",
    "- Changer la valeur dans un slice change la valeur dans le tableau.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "a = np.arange(10)\n",
    "print(a)\n",
    "b = a[3:6]\n",
    "print(b)  # `b` est une \"vue\" sur le sous-tableau de `a` et `a` est appelé la \"base\" de `b`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "b[0] = -1  # l'indice 0 de b est l'indice 3 de a\n",
    "print(a)   # si on change une vue, la base est aussi changée."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "c = a[7:8].copy() # copie explicite du sous-tableau, pas de problème\n",
    "c[0] = -1 \n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sous-tableaux avec des tableaux d'indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(10)\n",
    "print(a)\n",
    "\n",
    "indices = np.array([1, 3, 5, -1])  # Si on ne veut pas des indices équirépartis par exemple\n",
    "print(a[indices])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On peut modifier les valeurs du sous tableau :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a[indices] = -1\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mais **ce n'est pas une vue**. En effet, assigner le sous-tableau à une variable fait une copie."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = a[indices]\n",
    "b[:] = 10\n",
    "print(a)         # a n'est pas modifié"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tableau de masque"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- un masque est un tableau de booléens (`True`  ou `False`)\n",
    "- indexer un tableau avec un masque **de la même taille** retourne les éléments du tableau ou le masque contient `True`\n",
    "- un masque doit donc avoir la même taille que le tableau, il y a une erreur sinon\n",
    "- les masques peuvent être créer en utilisant les opérateurs de comparaisons (>, <, >=, <=, ==, ...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(10)\n",
    "print(a)\n",
    "mask = np.array([True if i%2==0 else False for i in range(a.size)])\n",
    "print(mask)\n",
    "print(a[mask])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a[mask] = -10 # On peut utiliser le masque pour modifier un tableau directement\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# En utilisant un opérateur de comparaison\n",
    "mask = a <= 3\n",
    "a[mask] = 0\n",
    "print(mask)\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ou bien directement sans intermédiaire\n",
    "a[a <= 3] = 1\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comme pour les tableaux d'indices, les sous-tableaux obtenus avec des masques ne sont pas des vues sur le tableau d'origine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = a[mask]\n",
    "b[:] = 42\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Tableaux multi-dimensions\n",
    "\n",
    "On a déjà vu des tableaux 2D au début !\n",
    "\n",
    "Voici quelques différences avec les listes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "a = np.arange(4*3).reshape(4, 3)                   # tableau numpy 2D\n",
    "l = [[0, 1, 2], [3, 4, 5], [6, 7, 8], [9, 10, 11]] # liste de liste en python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(f\"a = \\n{a}\\n\")\n",
    "print(f\"l = \\n{l}\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "l[-1][-1] # Accès au dernier élément"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(a[-1, -1])  # L'indexation des tableaux multi-D est différente avec les tableau numpy. On sépare les dimensions par une virgule.\n",
    "print(a[0, 0])    # retourne le premier élément\n",
    "print(a[1, :])    # retourne la deuxième ligne. : = toutes les colonnes, c´est un slice\n",
    "\n",
    "# Comme on peut le voir, on indique l'indice ou le slice pour chaque dimension que l'on sépare par une virgule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# La deuxième ligne avec un tableau 2D\n",
    "r = a[1]\n",
    "print(r, r.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# La dernière colonne (comme un tableau 1D)\n",
    "c = a[:, -1]\n",
    "print(c, c.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# La dernière colonne (comme un tableau 2d). On prend ici un tableau d'indice pour spécifier la colonne et non l'indice seul.\n",
    "c = a[:, [-1]]\n",
    "print(c, c.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(a)\n",
    "print()\n",
    "print(a[1:4, 0:3])    # avec des slices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(a)\n",
    "print()\n",
    "ind_rows = np.array([1, 2, 3])\n",
    "ind_cols = np.array([0, 1, 2])\n",
    "print(a[ind_rows, ind_cols])      # avec des tableaux d'indices\n",
    "                                  # ind_rows et ind_cols doivent avoir la même taille !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Conversion d'une liste (de liste (de liste (...))) vers un tableau numpy\n",
    "b = np.array(l)\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**EXERCICES**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Planche de Galton**\n",
    "\n",
    "On veut simuler une planche de Galton (qui illustre le théorème central limite).\n",
    "\n",
    "Une balle est lachée en haut de la pyramide :\n",
    "\n",
    "          o\n",
    "         o o\n",
    "        o o o\n",
    "       o o o o\n",
    "      o o o o o\n",
    "     o o o o o o\n",
    "     \n",
    "À chaque niveau, la balle a une probabilité $p$ d'aller à gauche et $p-1$ d'aller à droite. La pyramide a $n$ niveaux. On fait cette expérience avec $N$ balles, quel est la répartition des balles en bas de la pyramide ?\n",
    "\n",
    "1. Écrire une fonction `bernoulli(p)` qui retourne $1$ avec probabilité $p$ et $0$ avec probabilité $1-p$. On peut le faire en utilisant par exemple la fonction `rng.random` qui génère un nombre réel suivant une loi uniforme entre $0$ et $1$. Si ce nombre est plus petit que $p$, on retourne $1$, sinon on retourne $0$.\n",
    "2. Ajouter un argument $n$ à la fonction `bernoulli` et modifier la afin qu'elle retourne un tableau de $0$ et de $1$ de taille $n$ au lieu d'une seule valeur.\n",
    "3. Utiliser votre fonction `bernoulli` avec $n = 10000$ et $0 < p < 1$ et calculer la valeur moyenne du tableau obtenu. Vérifier que vous avez une valeur moyenne proche de la valeur moyenne théorique (moyenne d'une loi de Bernoulli de paramètre $p$). On pourra utiliser la fonction `np.mean`.\n",
    "4. Écrire une fonction `binomial(n, p)` qui appelle `bernoulli` et retourne le nombre de $1$ dans le tableau obtenu (qui correspond au nombre de fois où la balle est allée à gauche sur les $n$ niveaux de la pyramide).\n",
    "5. Modifier la fonction `binomial` pour qu'elle retourne un tableau de taille $N$, où $N$ est un nouvel argument de la fonction `binomial` (le nombre de balle), au lieu d'une seule valeur.\n",
    "6. Appeler votre fonction `binomial` avec les arguments $p=0.5$, $n=10000$ et $N=10000$ et mettre le résultat dans un tableau `b`. Vérifier que la valeur moyenne de `b` est proche de la valeur thérorique (moyenne d'une loi binomiale de paramètre $n, p$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load solutions/numpy/galton-board.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# En réalité, on vient de réécrire une fonction qui génère des nombres aléatoirement suivant une loi binomiale de paramètre n et p.\n",
    "# Mais elle existe déjà dans numpy\n",
    "rng = np.random.default_rng()\n",
    "b = rng.binomial(n, p, N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On aimerait voir la répartition des $N$ balles en bas de la pyramide. Pour cela, on pourrait tracer un histogramme du résultat. Voici quelques exemple d'histogramme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Un histogramme de valeurs aléatoires suivant une loi normale centrée réduite\n",
    "x = rng.normal(0, 1, 10000)\n",
    "plt.figure()\n",
    "plt.hist(x)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# On peut changer le nombre de barres\n",
    "plt.figure()\n",
    "plt.hist(x, bins=30)    # bins = \"auto\" marche et calcule un nombre adapté de façon automatique.\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Normalisation de l'histogramme (l'aire fait 1, comme une densité de probabilité)\n",
    "plt.figure()\n",
    "plt.hist(x, bins=\"auto\", density=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $\\mu = n p$ et $\\sigma = \\sqrt{n p (1 - p)}$. Afficher sur la même figure :\n",
    "\n",
    "- l'histogramme normalisé des $N$ valeurs de $(b - \\mu) / \\sigma$, où $b$ est le tableau `b` obtenu précédemment.\n",
    "- la fonction $x \\mapsto \\frac{1}{\\sqrt{2 \\pi}} \\exp(-\\frac{x^2}{2})$ sur l'intervalle $[-4, 4]$. C'est la densité de la loi normale centrée réduite (moyenne nulle, variance 1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load solutions/numpy/galton-board_plot.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Marche aléatoire**\n",
    "\n",
    "Une marche aléatoire en 1D consiste en une suite $(x_n)$ telle que $x_{n+1} = x_n + d$ où $d$ est choisi aléatoirement dans $\\{-1, 1\\}$.\n",
    "\n",
    "On suppose que l'on part toujours du point $x_0 = 0$.\n",
    "\n",
    "**Marche aléatoire 1D**\n",
    "\n",
    "1. Générer un tableau `d` contenant $N = 10000$ entiers dans $\\{-1, 1\\}$ (voir `rng.integers` ou `rng.choice`). Ce tableau représente les directions prisent à chaque étape $n$.\n",
    "2. Calculer la position $x_n$ à chaque étape $n$ avec $x_0 = 0$ à partir de `d` (voir `np.cumsum`) et les mettre dans un tableau `x`.\n",
    "3. Afficher la suite $n \\mapsto x_n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load solutions/numpy/random_walk.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Moyenne des positions**\n",
    "\n",
    "1. Générer $p = 2000$ marches aléatoires 1D de taille $N = 10000$.\n",
    "2. Pour chaque marche aléatoire, garder la **dernière position** dans un tableau `last_pos` (de taille $p$ donc).\n",
    "3. Calculer la moyenne des positions à chaque étape $\\left<x_n\\right>$ des marches aléatoires sur les $p$ realisations. Les mettre dans un tableau `mean_pos`.\n",
    "4. Calculer la moyenne quadratique des positions : $\\sqrt{\\left<x_n^2\\right>}$. Les mettres dans un tableau `mean2_pos`.\n",
    "5. Tracer les deux moyennes sur un même graphique. Dans une autre figure, tracer l'histogramme des dernières positions des marhces aléatoires."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load solutions/numpy/random_walk_stats.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Temps de retour à la position 0**\n",
    "\n",
    "1. Créer une marche aléatoire de taille $N = 10000$ dans un tableau `x`.\n",
    "2. Trouver le premier temps auquel la marche est revenue à la position $0$. On pourra utiliser les fonctions `np.nonzero` ou `np.argwhere`.\n",
    "3. Afficher ce premier temps de retour à 0 ou un message si ce n'est jamais le cas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load solutions/numpy/random_walk_return_to_0.py"
   ]
  },
  {
   "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
}