NumPy 2 – vectoriser, agréger, diffuser, masquer

Au programme aujourd'hui

Quiz 1 (10 min, papier) sur la séance 1, puis correction flash.

Notebooks :

• 1-04 vectorisation — la suite : ufunc, out=
• 1-07 agrégation et axis
• 1-08 broadcasting
• 1-09 tests et masques sur tableaux

TP : 1-10 TP images (2/2) — patchwork, sépia, somme & overflow

Devoir : D2 « Monte-Carlo & Mandelbrot », 2 semaines (pour le dimanche 27 septembre, 23 h 59)

Retour sur le quiz 1

Les trois confusions les plus fréquentes — à corriger tout de suite :

• shape est un tuple, pas un nombre : (5,) pour un vecteur de 5,
  (2, 5) pour 2 lignes × 5 colonnes. len(a) ne donne que la 1re dimension.
• reshape et le slicing ne copient pas : modifier la vue modifie l'original.
  « b = a[1:4] puis b[0] = 0 » change a.
• Le dtype ne s'élargit pas tout seul : int8 × 100 déborde, uint8 + uint8
  boucle à 256. C'est à vous de faire astype.

Ces trois pièges ne disparaissent jamais : ils peuvent revenir au quiz 2, qui porte sur la séance d'aujourd'hui.

Vectoriser plutôt que boucler

Pour appliquer une fonction à tous les éléments, on ne boucle jamais en Python :
on applique la fonction au tableau entier — c'est une ufunc (universal function).

x = np.linspace(0, 2 * np.pi, 1_000_000)

y = np.sin(x)                      # OUI : une ufunc, la boucle est en C
# y = [math.sin(e) for e in x]     # NON : 10 à 100 fois plus lent

+ - * / ** < == % abs sqrt exp … sont toutes des ufunc : élément par élément, même forme en sortie.
np.add(a, 1, out=a) écrit le résultat dans un tableau existant, sans allouer.

En effet, la boucle Python paie un objet par élément, quand la ufunc parcourt un buffer contigu.

Agréger : combiner les éléments

tab = np.arange(12).reshape(3, 4)
tab.sum()     # 66 — sans axis, on agrège TOUT : un scalaire
tab.max()     # 11

Toutes existent aussi en méthode : tab.sum(), tab.argmax()…

axis : la dimension qui disparaît

tab = np.arange(12).reshape(3, 4)
#  [[ 0  1  2  3]
#   [ 4  5  6  7]
#   [ 8  9 10 11]]
tab.sum(axis=0)   # [12 15 18 21]  shape (4,) : on a consommé les LIGNES
tab.sum(axis=1)   # [ 6 22 38]     shape (3,) : on a consommé les COLONNES
tab.argmax(axis=1)   # [3 3 3]  indice du max de chaque ligne

La règle : axis=k est la dimension le long de laquelle on agrège,
et elle disparaît du résultat : shape (3, 4) avec axis=0 donne (4,).

D'où vient l'erreur classique ? On pense « axis=0 = les lignes » et on attend 3 valeurs ;
or on écrase les lignes les unes sur les autres : il reste une valeur par colonne.

Le broadcasting : l'idée

Pour incrémenter tous les éléments, pas besoin de fabriquer un tableau de 1 de même forme :
NumPy élargit virtuellement le plus petit opérande, sans jamais le construire en mémoire.

mat = np.arange(15).reshape(3, 5)
mat + 1                          # le scalaire est diffusé partout
mat + [10, 20, 30, 40, 50]       # une ligne (5,) ajoutée à chacune des 3 lignes
mat + [10, 20, 30]               # ValueError: operands could not be
                                 #   broadcast together with shapes (3,5) (3,)

Attention, ce n'est pas le + des listes : [1, 2] + [10] vaut [1, 2, 10] en Python !

L'intérêt : « ajouter un vecteur à chaque ligne » est l'opération la plus courante du calcul numérique ; elle mérite une syntaxe sans boucle ni copie.

La règle : aligner par la droite

NumPy aligne les deux formes par la droite, puis compare dimension par dimension :

1. dimensions égales → conservées ;
2. l'une vaut 1 (ou est absente) → elle est étirée à l'autre ;
3. sinon → ValueError.

np.broadcast_shapes((3, 5), (3, 1)) vous donne la réponse sans rien calculer.

Ligne + colonne = matrice

Pour ajouter une colonne, il faut la forme (n, 1) : on crée la dimension manquante avec None (alias np.newaxis).

ligne = np.arange(5)              # (5,)    [0 1 2 3 4]
col = 10 * np.arange(3)           # (3,)    [0 10 20]
col[:, None].shape                # (3, 1)

ligne + col[:, None]              # (3, 5) : une table d'addition
# [[ 0  1  2  3  4]
#  [10 11 12 13 14]
#  [20 21 22 23 24]]

col.reshape(3, 1) fait la même chose. Retenez bien ce geste : c'est lui qui permettra,
en S6, de calculer toutes les distances entre N corps sans double boucle.

Masques booléens

Une comparaison est une ufunc : elle s'applique à tout le tableau et renvoie un tableau de booléens.

tab = np.arange(10)
tab % 2 == 0          # [ True False  True False ...]   un MASQUE
tab[tab % 2 == 0]     # [0 2 4 6 8]   sélection → une COPIE, à plat
tab[tab < 5] = 0      # à gauche d'un = : modifie EN PLACE

Un masque est aussi un tableau de 0/1 déguisé : on peut donc le compter.

m = tab > 5
m.sum()      # 4     combien de vrais
m.mean()     # 0.4   quelle proportion
m.any(), m.all()

Composer : &, |, ~ — et np.where

and, or et not ne sont pas vectorisés : Python demande « ce tableau est-il vrai ? » et NumPy refuse tout net.

(tab >= 3) and (tab < 7)   # ValueError: The truth value of an array ... is ambiguous
3 <= tab < 7               # même erreur : c'est un `and` caché

(tab >= 3) & (tab < 7)     # OUI — chaque sous-condition ENTRE PARENTHÈSES
(tab < 3) | (tab >= 7)
~(tab % 2 == 0)

np.where choisit élément par élément :

np.where(tab % 2 == 0, tab, -tab)   # [ 0 -1  2 -3  4 -5  6 -7  8 -9]
np.where(tab > 5)                   # (array([6, 7, 8, 9]),)  les INDICES des vrais

Reprise : vue ou copie ? Le piège n° 1

a = np.arange(12)
b = a.reshape(3, 4)   # b est une VUE de a
b[0, 0] = 99
a[0]                  # 99 — la mémoire est partagée !

reshape, le slicing de base (a[1:5], a[:, 0]), la transposition a.T renvoient des vues :
aucune donnée copiée, juste une autre façon de lire le même buffer.

np.shares_memory(a, b)        # True
np.shares_memory(a, a[1:5])   # True

Ce qui, au contraire, copie

c = a.copy()          # copie explicite
d = a[a > 5]          # indexation par MASQUE → copie
e = a[[0, 2, 4]]      # indexation par LISTE d'indices → copie
f = a.astype(float)   # changement de dtype → copie

d[0] = -1
a[6]                  # 6 : a n'a pas bougé

La règle : slicing de base et reshape → vue ; indexation fancy (masque, liste) et changement de type → copie.

Mais a[a > 5] = 0 modifie bien a : à gauche d'un =, NumPy écrit directement dans a.

En cas de doute : b = a.copy() avant de modifier, et vous êtes tranquille.

Les 3 pièges du jour

1. axis désigne ce qui disparaît. tab.sum(axis=0) sur (3, 4) donne 4 valeurs,
   pas 3. Vérifiez le shape du résultat avant de lui faire confiance.

2. Le broadcasting aligne par la droite. (3, 5) + (3,) échoue ; il faut (3, 1).
   Réflexe : col[:, None].

3. and / or / not ne marchent pas sur des tableaux. &, |, ~,
   et des parenthèses autour de chaque condition.

Bonus : tab[mask] est une copie, tab[mask] = v écrit dans tab. Les deux sont vrais.

TP du jour : images (2/2)

Notebook 1-10 sur https://numerique.info-mines.paris — on repart du 1-06.

1. patchwork : une grille de tuiles colorées, construite par slicing et broadcasting ;
2. sépia : un produit matriciel img @ M sur les canaux, puis np.clip ;
3. somme dans une image & overflow : additionner deux images uint8… et comprendre
   pourquoi ça noircit — astype(float) avant, astype(np.uint8) après.

À rendre : rien aujourd'hui (c'est D2 qui compte) — mais finissez le sépia, il tombe souvent en quiz.

Pour les rapides : 1-11 TP dés (avancé) · 1-14 indexation avancée ·
1-13 algèbre linéaire · exercices « numpy » sur https://numerique-exos.info-mines.paris.

Pour les très rapides : second niveau secret, le sujet bonus B2 – Le jeu de la vie.
Sur sa page ou celle des devoirs, tapez IDKFA — les joueurs de Doom connaissent.

Devoir pour la prochaine séance

D2 — Monte-Carlo & Mandelbrot

Pour le dimanche 27 septembre, 23 h 59 (deux semaines : pas de séance le 22).

• Monte-Carlo : estimer π en tirant n points au hasard
  (rng = np.random.default_rng(seed), rng.random((n, 2))), compter ceux du quart de disque
  avec un masque, tracer la convergence ;
• Mandelbrot : une grille complexe np.linspace + np.meshgrid, l'itération z = z*z + c
  sur toute la grille à la fois, un masque pour marquer l'échappement, une image HD.

Ce qu'on rend : le notebook exécuté + mandelbrot.png.

Ce qui est noté : aucune boucle sur les points ou les pixels (la boucle sur les itérations est
normale), graine fixée = résultats reproductibles, figures titrées.