Une introduction à Numpy

  • Tableaux Numpy
  • Création de tableaux
  • Opérations basiques
  • Indexation et slicing
  • Algèbre linéaire
  • Méthodes sur les tableaux

Contenu sous licence CC BY-SA 4.0, largement inspiré de https://github.com/pnavaro/python-notebooks

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 ufuncs : 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

In [1]:
import numpy as np
print(np.__version__)

1.26.2

Utilisez l'auto-complétion pour lister les objets numpy :

In [2]:
#np.nd<TAB>

La rubrique d'aide est également précieuse

In [3]:
?np.ndarray

Tableaux Numpy

Une question de performance

  • 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
In [4]:
from random import random
from operator import truediv
In [5]:
L1 = [random() for i in range(100000)]
L2 = [random() for i in range(100000)]
%timeit s = sum(map(truediv, L1, L2))

3.33 ms ± 66.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [6]:
a1 = np.array(L1)
a2 = np.array(L2)
%timeit s = np.sum(a1/a2)

124 µs ± 485 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

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)
In [7]:
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)
In [8]:
print(a, b, c)

[0 1 2 3] [4 5 6 7] [[8 9 0 1]]

Propriétés

In [9]:
a = np.array([1, 2, 3, 4, 5])  # On crée un tableau
In [10]:
type(a)  # On vérifie son type
Out[10]:
numpy.ndarray
In [11]:
a.dtype  # On affiche le type de ses éléments
Out[11]:
dtype('int64')
In [12]:
a.itemsize  # On affiche la taille mémoire (en octets) de chaque élément
Out[12]:
8
In [13]:
a.shape  # On retourne un tuple indiquant la longueur de chaque dimension
Out[13]:
(5,)
In [14]:
a.size  #  On retourne le nombre total d'éléments
Out[14]:
5
In [15]:
a.ndim   # On retourne le nombre de dimensions
Out[15]:
1
In [16]:
a.nbytes  # On retourne l'occupation mémoire
Out[16]:
40
  • Toujours utiliser shape ou size pour les tableaux numpy plutôt que len
  • len est réservé aux tableaux 1D

Création de tableaux Numpy

Avec des valeur constantes

In [17]:
x = np.zeros((5, 3))  # On ne précise pas le type : on crée des flottants
print(x)
print(x.dtype)

[[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]
float64
In [18]:
x = np.zeros((5, 3), dtype=int)  # On explicite le type
print(x)
print(x.dtype)

[[0 0 0]
 [0 0 0]
 [0 0 0]
 [0 0 0]
 [0 0 0]]
int64

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

Création à partir d'une séquence

arange

C'est l'équivalent de range pour les listes.

In [19]:
np.arange(5)  # entiers de 0 à 4
Out[19]:
array([0, 1, 2, 3, 4])
In [20]:
np.arange(5, dtype=np.double)  # flottants de 0. à 4.
Out[20]:
array([0., 1., 2., 3., 4.])
In [21]:
np.arange(2, 7)  # entiers de 2 à 6.
Out[21]:
array([2, 3, 4, 5, 6])
In [22]:
np.arange(2, 7, 0.5)  # flottants avec incrément de 0.5.
Out[22]:
array([2. , 2.5, 3. , 3.5, 4. , 4.5, 5. , 5.5, 6. , 6.5])

linspace et logspace

In [23]:
# 5 éléments régulièrement espacés entre 1 et 4, 1 et 4 inclus
np.linspace(1, 4, 5)
Out[23]:
array([1.  , 1.75, 2.5 , 3.25, 4.  ])
In [24]:
# 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)
Out[24]:
array([   10.        ,    56.23413252,   316.22776602,  1778.27941004,
       10000.        ])

Consulter l'aide contextuelle pour plus de fonctionnalités

In [25]:
?np.logspace

Exercice : créer les tableaux suivants

[100 101 102 103 104 105 106 107 108 109]

Astuce : np.arange()

In [26]:
# Votre code ici
In [27]:
# Solution
np.arange(100, 110)
Out[27]:
array([100, 101, 102, 103, 104, 105, 106, 107, 108, 109])
[-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()

In [28]:
# Votre code ici
In [29]:
# Solution
np.linspace(-2., 2., num=20, endpoint=False)
Out[29]:
array([-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])
[ 0.001	0.00129155 0.0016681 0.00215443 0.00278256 
  0.00359381 0.00464159 0.00599484 0.00774264 0.01]

Astuce : np.logspace()

In [30]:
# Votre code ici
In [31]:
# Solution
np.logspace(-3, -2, num=10)
Out[31]:
array([0.001     , 0.00129155, 0.0016681 , 0.00215443, 0.00278256,
       0.00359381, 0.00464159, 0.00599484, 0.00774264, 0.01      ])
[[ 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()

In [32]:
# Votre code ici
In [33]:
# Solution
np.tri(7, 5, k=1) - np.ones((7, 5))
Out[33]:
array([[ 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.]])

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.

In [34]:
x = np.arange(0.0, 5.0, 1.0)
y = x*10.
z = x*100.
In [35]:
np.savetxt('test.out', (x, y, z))
%pycat test.out
In [36]:
np.savetxt('test.out', (x, y, z), fmt='%1.4e')   # Notation exponentielle
%pycat test.out
In [37]:
np.loadtxt('test.out')
Out[37]:
array([[  0.,   1.,   2.,   3.,   4.],
       [  0.,  10.,  20.,  30.,  40.],
       [  0., 100., 200., 300., 400.]])

savetxt et loadtxt ont leurs correspondants binaires :

  • save : enregistrer un tableau dans un fichier binaire au format .npy
  • load: créer un tableau numpy à partir d'un fichier binaire numpy

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

In [38]:
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)

In [39]:
with h5.File('test.h5', 'r') as f:
    for field in f.keys():
        print(f"{field}: {f[field][()]}")

x: [0. 1. 2. 3. 4.]
y: [ 0. 10. 20. 30. 40.]
z: [  0. 100. 200. 300. 400.]

Opérations basiques sur les tableaux

Par défaut, Numpy réalise les opérations arithmétiques éléments par éléments

In [40]:
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)

[ 0  5 12 21]
[ 4  6  8 10]
In [41]:
print(a**2)

[0 1 4 9]
In [42]:
print(5 * a)
print(5 + a)

[ 0  5 10 15]
[5 6 7 8]
In [43]:
print(a < 2)

[ True  True False False]
In [44]:
print(np.cos(a*np.pi))  # Utilisation de ufunc

[ 1. -1.  1. -1.]

De nombreuses ufunc dans la doc officielle.

Indexation et slicing

Indexation

Les règles différent légèrement des listes pour les tableaux multi-dimensionnels

In [45]:
# 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])

[[0 1 2]
 [3 4 5]
 [6 7 8]]
5
In [46]:
# Indexation de la liste équivalente
liste = a.tolist()
print(liste)
print(liste[1][2])

[[0, 1, 2], [3, 4, 5], [6, 7, 8]]
5

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.

In [47]:
print(a)
print(a[2, :])  # Retourne la 3ème ligne
print(a[:, 1])  # Retourne la deuxième colonne

[[0 1 2]
 [3 4 5]
 [6 7 8]]
[6 7 8]
[1 4 7]

Attention : contrairement aux listes, le slicing de tableaux ne renvoie pas une copie mais constitue une vue du tableau.

In [48]:
b = a[:, 1]
b[0] = -1
print(a)

[[ 0 -1  2]
 [ 3  4  5]
 [ 6  7  8]]

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 :

In [49]:
b.base is a.base  # b.base retourne le tableau dont b est la vue
Out[49]:
True

Si on veut réaliser une copie d'un tableau, il faut utiliser la fonction copy()

In [50]:
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é.

In [51]:
print(b.base)
b[0] = 200
print(a)

None
[[ 0 -1  2]
 [ 3  4  5]
 [ 6  7  8]]

Exercice

Approcher la dérivée de $f(x) = \sin(x)$ par la méthode des différences finies :

$$\frac{\partial f}{\partial x} \approx \frac{f(x+\Delta x)-f(x)}{\Delta x}$$

Les valeurs seront calculées au milieu de deux abscisses discrètes successives.

In [52]:
# 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)
In [53]:
%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');
No description has been provided for this image
In [54]:
# 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)$");
No description has been provided for this image

Quelques opérations d'algèbre linéaire

In [55]:
A = np.array([[1, 1],
              [0, 1]])
B = np.array([[2, 0],
              [3, 4]])
print(A*B)           # produit élément par élément
print(A.dot(B))      # produit matriciel
print(np.dot(A, B))  # une autre syntaxe de produit matriciel
print(A@B)           # encore un autre produit matriciel

[[2 0]
 [0 4]]
[[5 4]
 [3 4]]
[[5 4]
 [3 4]]
[[5 4]
 [3 4]]
In [56]:
a = np.array([[1.0, 2.0],
              [3.0, 4.0]])
print(a)

[[1. 2.]
 [3. 4.]]
In [57]:
# Deux syntaxes équivalentes pour la transposition
print(a.transpose())
print(a.T)

[[1. 3.]
 [2. 4.]]
[[1. 3.]
 [2. 4.]]
In [58]:
print(np.linalg.inv(a))  # inversion de matrice

[[-2.   1. ]
 [ 1.5 -0.5]]
In [59]:
print(np.trace(a))  # trace

5.0
In [60]:
print(np.eye(2))     # "eye" représente "I", la matrice identité

[[1. 0.]
 [0. 1.]]
In [61]:
y = np.array([[5.], [7.]])  # Résoudre A*x = y
x = np.linalg.solve(a, y)
print(x)
print(a@x == y)

[[-3.]
 [ 4.]]
[[ True]
 [ True]]
In [62]:
j = np.array([[0.0, -1.0], [1.0, 0.0]])
print(np.dot(j, j))     # produit matriciel
print(np.linalg.eig(j))  # Extraction des valeurs propres

[[-1.  0.]
 [ 0. -1.]]
EigResult(eigenvalues=array([0.+1.j, 0.-1.j]), eigenvectors=array([[0.70710678+0.j        , 0.70710678-0.j        ],
       [0.        -0.70710678j, 0.        +0.70710678j]]))

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.

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.)

Références