1- Manipulation des tableaux et calculs matriciels

Le but de ce TP est que vous vous fassiez la main sur la manipulation de tableaux sous numpy. Pour écrire du code efficace, il est absolument essentiel de manipuler les opérations vectorisées et de limiter le recours à des boucles pour itérer sur les valeurs des tableaux. Rappelez vous, les opérations vectorisées exécutent directement du code compilé alors que si vous utilisez des boucles, vous ajoutez une sur-couche interprétée qui est très pénalisante sur les temps de calcul.

En terme de resources, je vous conseille :

Pour réaliser les exercices ci-dessous, on pourra se créer un environnement virtuel, par exemple en utilisant uv :

Pour installer uv :

curl -LsSf https://astral.sh/uv/install.sh | sh

Pour créer un environnement virtuel, par exemple avec python 3.11:

uv venv --python 3.11  /tmp/venv
source /tmp/venv/bin/activate

Pour installer les dépendances de ce TP

uv pip install matplotlib numpy Pillow

Clustering avec K-Means

Fonctions/concepts utiles : broadcast, numpy.linalg.norm, numpy.argmin,

Le clustering consiste à résumer un nombre potentiellement infini d'observations en un ensemble discret de représentants ou prototypes.

Supposons que $x_i \in \mathbb{R}^d$ soit tiré d'une distribution de probabilité inconnue. Prenons $K$ prototypes $v_k \in \mathbb{R}^d, k \in [0,K-1]$ initialisés aléatoirement. L'algorithme de K-moyennes en ligne itère les étapes suivantes :

Trouver le prototype $k^\star$ le plus proche de $x_i$ : $k^\star = argmin_{k \in [0, K-1]}(|v_k - x_i|_2^2)$. On l'appelle l'Unité de Meilleure Correspondance (Best Matching Unit en anglais).
Mettre à jour la BMU selon $v_{k^\star} \leftarrow v_{k^\star} + \epsilon (x_i - v_{k^\star})$

Un code de base vous est fourni kmeans.py et il vous est demandé d’implémenter l’algorithme ci-dessus de manière efficace avec Numpy.

Le code implémenté produira la sortie suivante.

Extension Un exercice intéressant est de coder la version "batch" de l’algorithme. Cet algorithme fonctionne ainsi. Étant donné une collection d’échantillons $x_i$ :

Pour chaque prototype $k$, trouver tous les échantillons pour lesquels $k$ est le plus proche
Calculer la nouvelle position $v_k$ comme la moyenne des échantillons précédemment identifiés

Dans ce cas, pour être très efficace, vous devriez envisager d’utiliser la fonction cdist pour calculer les distances entre deux collections de points (prototypes et échantillons). Sur mon ordinateur portable, la version batch est $10$ fois plus rapide que la version en ligne.

Evaluation d'une fonction sur une grille régulière

Fonctions/concepts utiles : linspace, meshgrid et opérations vectorisées cos, ...

J’ai écrit un script Python pour tracer la fonction de Griewank. Mais il est assez lent. Pourriez-vous m’aider à améliorer ses performances ?

La fonction de Griewank s’écrit :

\[ f(x, y) = \frac{x^2 + y^2}{4000} - \cos(x) \cos\left(\frac{y}{\sqrt{2}}\right) + 1 \]

Sur mon ordinateur portable, pour une grille de $1000 \times 1000$, cela prend 1 seconde pour générer le tracé. Et en passant à une grille de $10,000 \times 10,000$, cela prend 48 secondes.

Voici le script : regular_grid.py.

Le résultat attendu est illustré ci-dessous.

Sur la grille $1000 \times 1000$, vous pouvez facilement diviser le temps d’exécution par deux.

Traitement d'image

Fonctions/concepts utiles : - indexation, indexation avancée, slicing, - numpy.pad, - numpy.logical_not

Une image RGB peut être représentée par un tableau numpy tridimensionnel. Par exemple, l’image bunny.jpg est une image couleur de taille $1500 \times 2122$ et est représentée par un tableau numpy de taille $1500 \times 2122 times 3$.

Nous allons travailler sur l’implémentation de quelques outils de traitement d’image. Étant donnée une image couleur RGB, je souhaite calculer sa luminance puis son Laplacien discret 2D. Les étapes sont les suivantes :

Pour le calcul de la luminance

Lorsqu'une image est chargée, par exemple avec PIL, vous pouvez la convertir en tableau numpy 3D avec un type de données unsigned int 8, c’est-à-dire des valeurs dans l’intervalle $[0, 255]$. Vous devez d’abord la convertir en flottant dans l’intervalle $[0, 1]$.
Ensuite, vous pouvez calculer la luminance $Y$ à partir des valeurs normalisées de couleur $R$, $G$ et $B$, selon les formules suivantes :

\[ \begin{eqnarray} f(c) &=& \begin{cases} \frac{c}{12.92} & \mbox{ si } c \leq 0.04045 \\ \left(\frac{c + 0.055}{1.055}\right)^{2.4} & \mbox{ sinon } \end{cases}\\ R_{lin} &=& f(R) \\ G_{lin} &=& f(G) \\ B_{lin} &=& f(B) \\ Y &=& 0.2126 R_{lin} + 0.7152 G_{lin} + 0.0722 B_{lin} \end{eqnarray} \]

Pour effectuer ce calcul efficacement, il est nécessaire d’utiliser des opérations vectorisées (comme multiplier un tableau par un scalaire), et l’indexation d’un tableau numpy avec un masque booléen (via l’indexation avancée).

Pour le calcul du laplacien

Étant donné la luminance $Y$ qui est un tableau numpy bidimensionnel, le Laplacien discret 2D est un opérateur de détection de contours et se calcule comme suit :

\[ \forall i,j , L_{i,j} = 4 Y_{i,j} - (Y_{i-1,j} + Y_{i+1,j} + Y_{i, j-1} + Y_{i, j+1}) \]

Pour effectuer ce calcul efficacement, il faut :

utiliser le slicing du tableau $Y$ pour calculer les tableaux « décalés » tels que $Y_{i-1,j}$, $Y_{i,j+1}$, etc. ;
être attentif aux bords. Pour ce dernier point, le padding peut être une solution. En fait, une solution efficace consiste à construire 4 tableaux : $T_{i,j} = Y_{i-1, j}$, $B_{i,j} = Y_{i+1, j}$, $L_{i,j} = Y_{i, j-1}$ et $R_{i,j} = Y_{i, j+1}$, de sorte que votre calcul s’écrive :

\[ \forall i, j \quad L_{i,j} = 4 Y_{i,j} - (T_{i,j} + B_{i,j} + L_{i,j} + R_{i,j}) \]

La difficulté réside dans la gestion des bords, mais une solution simple consiste à appliquer un padding de zéros sur le tableau $Y$ avant d’effectuer le slicing.

Un code de base vous est fourni : laplacian.py, et il vous est demandé d’implémenter de manière efficace avec Numpy l’algorithme ci-dessus. Il est strictement interdit d’utiliser une seule boucle for Python.

Equation d'onde

Nous voulons simuler l’équation des ondes. L’équation des ondes calcule le déplacement $u(x,y,t)$ à la position $(x,y)$ et au temps $t$ en utilisant l’équation différentielle suivante :

\[ \frac{\partial^2 u}{\partial t^2}(x, y, t) = c^2 \left(\frac{\partial^2 u}{\partial x^2}(x,y,t) + \frac{\partial^2 u}{\partial y^2}(x,y,t)\right) \]

Pour simuler ce système dans l’espace et dans le temps, nous allons discrétiser à la fois l’espace et le temps en utilisant les différences finies. Le système discrétisé s’écrit :

\[ \frac{u(x,y, t+\Delta t) - 2u(x,y,t) + u(x,y,t-\Delta t)}{\Delta t^2} = c^2 \left( \frac{u(x + \Delta x,y, t) - 2u(x,y,t) + u(x-\Delta x,y,t)}{\Delta x^2} + \frac{u(x,y + \Delta y, t) - 2u(x,y,t) + u(x,y- \Delta y,t)}{\Delta y^2} \right) \]

Nous ajoutons des conditions aux bords pour conserver toute l’énergie dans notre volume simulé. Les conditions aux limites stipulent qu’il n'y a pas de mouvement sortant. Cela est mis en œuvre en fixant explicitement les dérivées spatiales du premier ordre à $0$, en utilisant des "cellules fantômes" sur les bords qui reflètent les valeurs calculées.

En Numpy, nous devons stocker les valeurs du déplacement aux instants $t-\Delta, t, t+\Delta$ pour chaque position discrète de la zone simulée. Cela peut être représenté par un tableau de taille $N_x \times N_y \times 3$ et mis à jour selon les équations ci-dessus.

Un code de base vous est fourni : waves.py, et il vous est demandé d’implémenter de manière efficace avec Numpy l’algorithme ci-dessus. Vous devez implémenter les dérivées spatiales ainsi que les mises à jour temporelles.

Les dérivées spatiales peuvent être calculées comme dans l’exemple de traitement d’image vu précédemment, en définissant des sous-tableaux qui sont des slices du signal d’origine paddé. Notez que, comme indiqué dans le code, l’entrée des fonctions de différences finies est déjà paddée. La fonction step de l’objet WaveEquation vous demande d’implémenter le décalage temporel (assez direct) ainsi que la dérivée seconde temporelle. La seule difficulté liée à cette dérivée temporelle est de noter que la variable displacement est paddée, tandis que d2x et d2y ne le sont pas. Il faut donc uniquement mettre à jour la partie interne (sans la marge) du tableau displacement.

$Équation des ondes simulée sur une grille $100 \times 100$$

Conclusion

La morale de l'histoire c'est de bien penser ses algorithmes en terme d'opérations vectorisées. Dès que vous introduisez une boucle for en python, vous devez vous inquiéter de son impact néfaste sur les performances de votre implémentation.

Les solutions sont disponibles ci-dessous :