Une introduction à Numpy¶
- Tableaux Numpy
- Création de tableaux
- Opérations basiques
- Indexation et slicing
- Algèbre linéaire
- Méthodes sur les tableaux
- Broadcasting de tableaux
Contenu sous licence CC BY-SA 4.0, largement inspiré de https://github.com/pnavaro/python-notebooks et https://github.com/jakevdp/PythonDataScienceHandbook
Numpy¶
Le Python pur est peu performant pour le calcul. Les listes ne sont pas des objets efficaces pour représenter les tableaux numériques de grande taille. Numpy a été créé à l'initiative de développeurs qui souhaitaient combiner la souplesse du langage python et des outils de calcul algébrique performants.
Numpy se base sur :
- le
ndarray
: tableau multidimensionnel - des objets dérivés comme les masked arrays et les matrices
- les
ufunc
s : opérations mathématiques optimisées pour les tableaux - des méthodes pour réaliser des opération rapides sur les tableaux :
- manipulation des shapes
- tri
- entrées/sorties
- FFT
- algébre linéaire
- statistiques
- calcul aléatoire
- et bien plus !
Numpy permet de calculer à la Matlab en Python. Il est un élément de base de l'écosystème SciPy
Démarrer avec Numpy¶
import numpy as np
print(np.__version__)
Utilisez l'auto-complétion pour lister les objets numpy :
#np.nd<TAB>
La rubrique d'aide est également précieuse
np.ndarray?
- Les listes Python sont trop lentes pour le calcul et utilisent beaucoup de mémoire
- Représenter des tableaux multidimensionnels avec des listes de listes devient vite brouillon à programmer
from random import random
from operator import truediv
L1 = [random() for i in range(100000)]
L2 = [random() for i in range(100000)]
%timeit s = sum(map(truediv, L1, L2))
a1 = np.array(L1)
a2 = np.array(L2)
%timeit s = (a1/a2).sum()
La différence avec les listes¶
Les différences entre tableaux Numpy et listes Python :
- un
ndarray
a une taille fixée à la création - un
ndarray
est composé d'éléments du même type - les opérations sur les tableaux sont réalisées par des routines C pré-compilées et optimisées (éventuellement parallèles)
a = np.array([0, 1, 2, 3]) # list
b = np.array((4, 5, 6, 7)) # tuple
c = np.matrix('8 9 0 1') # string (syntaxe matlab)
print(a, b, c)
Propriétés¶
a = np.array([1, 2, 3, 4, 5]) # On crée un tableau
type(a) # On vérifie son type
a.dtype # On affiche le type de ses éléments
a.itemsize # On affiche la taille mémoire (en octets) de chaque élément
a.shape # On retourne un tuple indiquant la longueur de chaque dimension
a.size # On retourne le nombre total d'éléments
a.ndim # On retourne le nombre de dimensions
a.nbytes # On retourne l'occupation mémoire
- Toujours utiliser
shape
ousize
pour les tableaux numpy plutôt quelen
len
est réservé aux tableaux 1D
x = np.zeros((5, 3)) # On ne précise pas le type : on crée des flottants
print(x)
print(x.dtype)
x = np.zeros((5, 3), dtype=int) # On explicite le type
print(x)
print(x.dtype)
On dispose d'une panoplie de fonctions pour allouer des tableaux avec des valeurs constantes ou non initialisées (empty
) :
empty, empty_like, ones, ones_like, zeros, zeros_like, full, full_like
np.arange(5) # entiers de 0 à 4
np.arange(5, dtype=np.double) # flottants de 0. à 4.
np.arange(2, 7) # entiers de 2 à 6.
np.arange(2, 7, 0.5) # flottants avec incrément de 0.5.
linspace
et logspace
¶
# 5 éléments régulièrement espacés entre 1 et 4, 1 et 4 inclus
np.linspace(1, 4, 5)
# 5 éléments régulièrement espacés selon une progression géométrique entre 10^1 et 10^4
np.logspace(1, 4, 5)
Consulter l'aide contextuelle pour plus de fonctionnalités
np.logspace?
Exercice : créer les tableaux suivants¶
[100 101 102 103 104 105 106 107 108 109]
Astuce :
np.arange()
# Votre code ici
# Solution
np.arange(100, 110)
[-2. -1.8 -1.6 -1.4 -1.2 -1. -0.8 -0.6 -0.4 -0.2 0.
0.2 0.4 0.6 0.8 1. 1.2 1.4 1.6 1.8]
Astuce :
np.linspace()
# Votre code ici
# Solution
np.linspace(-2., 2, num=20, endpoint=False)
[ 0.001 0.00129155 0.0016681 0.00215443 0.00278256
0.00359381 0.00464159 0.00599484 0.00774264 0.01]
Astuce :
np.logspace()
# Votre code ici
# Solution
np.logspace(-3, -2, num=10)
[[ 0. 0. -1. -1. -1.]
[ 0. 0. 0. -1. -1.]
[ 0. 0. 0. 0. -1.]
[ 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0.]]
Astuce :
np.tri()
,np.ones()
# Votre code ici
# Solution
np.tri(7, 5, k=1) - np.ones((7, 5))
Création de tableaux à partir de fichiers¶
Afin d'illustrer la création d'un tableau numpy à partir de données lues dans un fichier texte, on commence par sauvegarder un tableau dans un fichier.
x = np.arange(0.0, 5.0, 1.0)
y = x*10.
z = x*100.
np.savetxt('test.out', (x, y, z))
%pycat test.out
np.savetxt('test.out', (x, y, z), fmt='%1.4e') # Notation exponentielle
%pycat test.out
np.loadtxt('test.out')
Format HDF5 avec H5py¶
Le format .npy
n'est lisible que par Numpy.
À l'inverse, le format HDF5 est partagé par un grand nombre d'applications.
De plus, il permet de structurer des données binaires :
- en les nommants
- en ajoutant des métadonnées
- en assurant une portabilité indépendante de la plateforme
voir le manuel utilisateur
H5py est une interface Python pour HDF5.
On écrit dans le fichier test.h5
import h5py as h5
with h5.File('test.h5', 'w') as f:
f['x'] = x
f['y'] = y
f['z'] = z
On lit le fichier test.h5
(qui pourrait provenir d'une autre application)
with h5.File('test.h5', 'r') as f:
for field in f.keys():
print(f"{field}: {f[field][()]}")
Opérations basiques sur les tableaux¶
Par défaut, Numpy réalise les opérations arithmétiques éléments par éléments
a = np.array([0, 1, 2, 3])
b = np.array((4, 5, 6, 7))
print(a * b) # Produit éléments par éléments : pas le produit matriciel !
print(a + b)
print(a**2)
print(5 * a)
print(5 + a)
print(a < 2)
print(np.cos(a*np.pi)) # Utilisation de ufunc
De nombreuses ufunc dans la doc officielle.
# Indexation d'un tableau numpy
a = np.arange(9).reshape(3, 3) # Notez l'effet de la méthode reshape()
print(a)
print(a[1, 2])
# Indexation de la liste équivalente
liste = a.tolist()
print(liste)
print(liste[1][2])
Slicing¶
Pour les tableaux unidimensionnels, les règles de slicing sont les mêmes que pour les séquences. Pour les tableaux multidimensionnels numpy, le slicing permet d'extraire des séquences dans n'importe quelle direction.
print(a)
print(a[2, :]) # Retourne la 3ème ligne
print(a[:, 1]) # Retourne la deuxième colonne
Attention : contrairement aux listes, le slicing de tableaux ne renvoie pas une copie mais constitue une vue du tableau.
b = a[:, 1]
b[0] = -1
print(a)
a
est aussi une vue du tableau np.arange(9)
obtenue avec la méthode reshape()
donc a
et b
sont deux vues différentes du même tableau :
b.base is a.base # b.base retourne le tableau dont b est la vue
Si on veut réaliser une copie d'un tableau, il faut utiliser la fonction copy()
b = a[:, 1].copy()
ici b
n'est pas une vue mais une copie donc b.base
vaut None
et a
n'est pas modifié.
print(b.base)
b[0] = 200
print(a)
Exercice¶
Approcher la dérivée de $f(x) = \sin(x)$ par la méthode des différences finies centrées :
$$\frac{\partial f}{\partial x} \approx \frac{f \left(x + \frac{\Delta x}{2}\right)-f \left(x - \frac{\Delta x}{2} \right)}{\Delta x}$$
Les valeurs seront calculées au milieu de deux abscisses discrètes successives.
# On crée un tableau de 40 points d'abscisse
x, dx = np.linspace(0., 4.*np.pi, 40, retstep=True)
y = np.sin(x)
%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(x, np.cos(x)); # la dérivée analytique de sin() est cos()
# Votre code ici
# Décommentez la dernière ligne pour tracer la dérivée numérique dy_dx
# en fonction du milieu de 2 abscisses x_mil
# plt.plot(x_mil, dy_dx, 'rx');
# Solution
%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(x, np.cos(x)) # la dérivée analytique de sin() est cos()
x_mil = 0.5*(x[1:] + x[:-1])
dy = y[1:] - y[:-1]
dy_dx = dy/dx
plt.plot(x_mil, dy_dx, 'rx') # x_mil est le milieu de deux abscisses
plt.title(r"$\rm{Derivative\ of}\ \sin(x)$")
plt.xlabel("x")
plt.ylabel(r"$\frac{\partial f}{\partial x}$");
Quelques opérations d'algèbre linéaire¶
A = np.array([[1, 1],
[0, 1]])
B = np.array([[2, 0],
[3, 4]])
# Produit élément par élément
print(A*B)
# Produit matriciel (3 syntaxes équivalentes)
print(A.dot(B))
print(np.dot(A, B))
print(A@B)
a = np.array([[1., 2.],
[3., 4.]])
print(a)
# Deux syntaxes équivalentes pour la transposition
print(a.transpose())
print(a.T)
print(np.linalg.inv(a)) # inversion de matrice
print(np.trace(a)) # trace
print(np.eye(2)) # "eye" représente "I", la matrice identité
y = np.array([[5.], [7.]]) # Résoudre A*x = y
x = np.linalg.solve(a, y)
print(x)
print(a@x == y)
j = np.array([[0., -1.], [1., 0.]])
print(np.dot(j, j)) # produit matriciel
print(np.linalg.eig(j)) # Extraction des valeurs propres
Les méthodes associées aux tableaux¶
Elles sont très nombreuses : impossible de toutes les lister dans le cadre de ce cours. Citons brièvement :
min
,max
,sum
sort
,argmin
,argmax
- statistiques basiques :
cov
,mean
,std
,var
À chercher dans la doc officielle.
Broadcasting de tableaux¶
Le broadcasting est un mécanisme qui permet de faire des opérations sur des tableaux de différentes tailles ou shapes.
Pour rappel, les opérations sur les tableaux de même taille sont réalisées élément par élément :
a = np.array([0, 1, 2])
b = np.array([1, 1, 1])
a + b
Si on remplace b
par le scalaire 1
, numpy déclenche le mécanisme de broadcasting où ce scalaire est virtuellement transformé en tableau de même taille que a
en répétant la valeur, ce qui aboutit au même résultat :
a + 1
Si on additionne maintenant un tableau 1D et un tableau 2D :
A = np.ones((3, 3), dtype=int)
A
a
A + a
a
est étendu suivant la 2e dimension jusqu'à couvrir la shape de A
.
Considérons maintenant un tableau ligne et un tableau colonne :
a = np.arange(3)
a
b = np.arange(3).reshape(3, 1)
b
a + b
a
et b
sont étendus simultanément pour couvrir la shape (3, 3)
.
Règles de broadcasting¶
- Si les deux tableaux n'ont pas la même dimension, la shape de celui avec la plus faible dimension est complétée avec des $1$ sur son bord gauche.
- Si les shapes des tableaux n'ont aucune valeur commune, le tableau dont la shape vaut $1$ est étendue jusqu'à la valeur de l'autre tableau.
- Si les tailles sont différentes et non-unitaires dans toutes les dimensions, alors le broadcasting est impossible et on lève une erreur.
Exemple 1 : on étend un tableau¶
A = np.ones((2, 3), dtype=int)
a = np.arange(3)
print(f"{A.shape = }")
print(f"{a.shape = }")
- Selon la règle 1., la shape de
a
devient(1, 3)
. - Selon la règle 2., la première dimension de
a
est étendue à(2, 3)
.
Maintenant que les shapes correspondent, l'addition terme à terme est possible :
A + a
Exemple 2 : on étend les deux tableaux¶
a = np.arange(3).reshape((3, 1))
print(a)
print(f"{a.shape = }")
b = np.arange(3)
print(b)
print(f"{b.shape = }")
Dans ce cas, les deux tableaux doivent être étendus.
- Selon la règle 1., la shape de
b
doit être étendue en(1, 3)
. - Selon la règle 2., on étend la shape de
a
etb
en prenant le maximum de la taille dans chaque dimension soit(3, 3)
.
a + b
Exemple 3 : broadcasting impossible¶
a = np.ones((3, 2), dtype=int)
print(a)
print(f"{a.shape = }")
b = np.arange(3)
print(b)
print(f"{b.shape = }")
- Selon la règle 1., la shape de
b
devient(1, 3)
. - Selon la règle 2., la shape de
b
devient(3, 3)
. - Selon la règle 3., le broadcasting est alors impossible :
a + b
Exemple d'utilisation¶
Un cas classique d'utilisation du broadcasting consiste à retirer la moyenne dans certaines dimensions :
rng = np.random.default_rng(seed=0)
X = rng.random((6, 3))
X
On calcule la moyenne selon la première dimension :
Xmean = X.mean(0)
Xmean
On peut maintenant calculer la valeur centrée de X
en retirant la moyenne :
X_centered = X - Xmean
X_centered
On vérifie que X_centered
est à moyenne nulle suivant la première dimension :
X_centered.mean(0)
Programmation avec Numpy¶
- Les opérations sur les tableaux sont rapides, les boucles python sont lentes => Éviter les boucles !
- C'est une gymnastique qui nécessite de l'entraînement
- Le résultat peut devenir difficile à lire et à débugguer, par exemple dans le cas de boucles contenant de multiples conditions
- D'autres options sont alors envisageables (Cython, Pythran, Numba, etc.)