=================================== Numpy : calcul vectoriel en Python =================================== .. topic:: Conventions :: >>> import numpy as np >>> import scipy as sp >>> import pylab as pl Le calcul avec des tableaux ============================ +--------------------------+-------------------------------------+ | Python | numpy | +--------------------------+-------------------------------------+ | List: `a = [1, 2, 3]` | Tableau: `a = np.array([1, 2, 3])` | +--------------------------+-------------------------------------+ Faire des opérations sur beaucoup de nombres --------------------------------------------- * Calcul numérique classique = boucles :: def square(data): for i in range(len(data)): data[i] = data[i]**2 return data .. sourcecode:: ipython In [1]: %timeit data = range(1000) ; square(data) 1000 loops, best of 3: 314 us per loop * Calcul vectoriel: remplacer les boucles par des opérations sur des tableaux :: def square(data): return data**2 .. sourcecode:: ipython In [2]: %timeit data=np.arange(1000) ; square(data) 100000 loops, best of 3: 10.6 us per loop Bénéfice du polymorphisme: même code pour un seul nombre, ou un tableau. Des objets multi-dimensionels ----------------------------- :: >>> a = np.arange(10) >>> a array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) >>> b = np.reshape(a, (2, 5)) >>> b array([[0, 1, 2, 3, 4], [5, 6, 7, 8, 9]]) >>> b[:, 1] array([1, 6]) Création de tableaux --------------------- * Création avec des constantes:: >>> np.ones((2, 3)) array([[ 1., 1., 1.], [ 1., 1., 1.]]) * Les tableaux contienent des entrées typées uniformement:: >>> np.ones(3, dtype=np.int) array([1, 1, 1]) * Création d'une grille:: >>> x, y = np.indices((2, 2)) >>> x array([[0, 0], [1, 1]]) >>> y array([[0, 1], [0, 1]]) >>> x+1j*y array([[ 0.+0.j, 0.+1.j], [ 1.+0.j, 1.+1.j]]) Slicing --------- Parcourt multidimensionel des tableaux .. image:: numpy_indexing_fr.png :align: center Un exemple d'application: calcul du laplacien ---------------------------------------------- .. image:: laplacien.jpg :align: center :: image[1:-1, 1:-1] = (image[:-2, 1:-1] - image[2:, 1:-1] + image[1:-1, :-2] - image[1:-1, 2:])*0.25 .. sourcecode:: ipython In [3]: import pylab as pl In [4]: l = sp.lena() In [5]: pl.imshow(l, cmap=pl.cm.gray()) In [6]: e = l[:-2, 1:-1] - l[2:, 1:-1] + l[1:-1, :-2] - l[1:-1, 2:] In [7]: pl.imshow(e, pl.cm.gray()) .. plot:: lena_laplacien.py :hide-links: :align: center ____ .. figure:: timings.jpg :align: center **Gains en temps** Indexage avancé ================ Avec des entiers ou des masques .. image:: numpy_fancy_indexing_fr.png :align: center Utilisation de tableaux d'entiers ------------------------------------ .. >>> np.random.seed(4) * Trier un vecteur avec un autre:: >>> a, b = np.random.random_integers(10, size=(2, 4)) >>> a array([8, 6, 2, 9]) >>> b array([ 8, 9, 3, 10]) >>> a_order = np.argsort(a) >>> a_order array([2, 1, 0, 3]) >>> b[a_order] array([ 3, 9, 8, 10]) Utilisation de masques ------------------------- * Mettre tous les éléments pair d'un tableau à zéro:: >>> a = np.arange(10) >>> a array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) >>> a[a % 2] = 0 >>> a array([0, 0, 2, 3, 4, 5, 6, 7, 8, 9]) * Masque à partir d'une grille pour selectionner le centre de l'image: .. sourcecode:: ipython In [8]: n, m = l.shape In [9]: x, y = np.indices((n, m)) In [10]: distances = np.sqrt((x - 0.5*n)**2 + (y - 0.5*m)**2) In [11]: l[distance > 200] = 255 In [12]: pl.imshow(l, cmap=pl.cm.gray) .. plot:: lena_mask.py :hide-links: :scale: 75 :align: center Le broadcasting ================ Opérations multidimensionelles ------------------------------- * Ajouter un nombre et un tableau marche:: >>> a = np.ones((3, )) >>> a array([ 1., 1., 1.]) >>> a + 1 array([ 2., 2., 2.]) * Et si on ajoute deux tableaux de forme différentes? :: >>> b = 2*np.ones((2, 1)) >>> b array([[ 2.], [ 2.]]) >>> a + b array([[ 3., 3., 3.], [ 3., 3., 3.]]) * Les dimensions se complètent:: .. image:: broadcasting.jpg :align: center Pour la performance ------------------- * Creation d'une grille 3D .. image:: 3d_radius.jpg :align: center :: np.sqrt(x**2 + y**2 + z**2) ____ .. image:: 3d_radius_non_broadcasting.jpg :align: right :Sans broadcasting: .. raw:: latex \rule{0pt}{2em} :: >>> x, y, z = np.mgrid[-100:100, -100:100, -100:100] >>> print x.shape, y.shape, z.shape (200, 200, 200) (200, 200, 200) (200, 200, 200) >>> r = np.sqrt(x**2 + y**2 + z**2) * Temps : **2.3s**: création de `x`, `y`, `z`: 0.5s, calcul de `r`: 1.8s * Mémoire : 64Mo par tableaux, 6 tableaux, (`x`, `y`, `z`, `r`) et 2 temporaires => **400Mb** * 200^3 opérations élémentaires par opération de tableaux: **48 million d'opérations**. ____ .. image:: 3d_radius_broadcasting.jpg :align: right :Avec broadcasting: .. raw:: latex \rule{0pt}{2em} :: >>> x, y, z = np.ogrid[-100:100, -100:100, -100:100] >>> print x.shape, y.shape, z.shape (200, 1, 1) (1, 200, 1) (1, 1, 200) >>> r = np.sqrt(x**2 + y**2 + z**2) * Temps : **1.1s**: création de `x`, `y`, `z`: 6ms * Mémoire: `x`, `y`, `z` : 1.6Kb. `r` : 64Mo, et 1 temporaire de 64Mo => **120Mb** * **16 million d'opérations** .. raw:: html   .. topic:: `numpy`: une view structurée sur la mémoire avec des opérations * données identiques (`dtype`) * indexage rapide * vue/copies * reshape pas couteux * opérations comprenant la forme des tableaux .. :vim:spell: :vim:spelllang=fr: