Notes de cours Python scientifique

I-A-1. Les tableaux NumPy▲

>>> import numpy as np
>>> a = np.array([0, 1, 2, 3])
>>> a
array([0, 1, 2, 3])

In [1]: L = range(1000)

In [2]: %timeit [i**2 for i in L]
1000 loops, best of 3: 403 us per loop

In [3]: a = np.arange(1000)

In [4]: %timeit a**2
100000 loops, best of 3: 12.7 us per loop

I-A-2. La documentation de référence de NumPy▲

>>> np.lookfor('create array') 
Search results for 'create array'
---------------------------------
numpy.array
    Create an array.
numpy.memmap
    Create a memory-map to an array stored in a *binary* file on disk.

In [6]: np.con*?
np.concatenate
np.conj
np.conjugate
np.convolve

I-A-3. Les conventions d'import▲

>>> import numpy as np

I-B-1. Construction manuelle de tableaux▲

>>> a = np.array([0, 1, 2, 3])
>>> a
array([0, 1, 2, 3])
>>> a.ndim
1
>>> a.shape
(4,)
>>> len(a)
4

>>> b = np.array([[0, 1, 2], [3, 4, 5]])    # 2 x 3 array
>>> b
array([[0, 1, 2],
       [3, 4, 5]])
>>> b.ndim
2
>>> b.shape
(2, 3)
>>> len(b)     # returns the size of the first dimension
2

>>> c = np.array([[[1], [2]], [[3], [4]]])
>>> c
array([[[1],
        [2]],

       [[3],
        [4]]])
>>> c.shape
(2, 2, 1)

I-B-2. Fonctions pour créer des tableaux▲

>>> a = np.random.rand(4)       # uniform in [0, 1]
>>> a  
array([ 0.95799151,  0.14222247,  0.08777354,  0.51887998])

>>> b = np.random.randn(4)      # Gaussian
>>> b  
array([ 0.37544699, -0.11425369, -0.47616538,  1.79664113])

>>> np.random.seed(1234)        # Setting the random seed

Vous avez peut-être remarqué que dans certains cas, les éléments d'un tableau sont affichés au format float (2. au lieu de 2). Cela est dû à la gestion du type de données :

>>> a = np.array([1, 2, 3])
>>> a.dtype
dtype('int64')

>>> b = np.array([1., 2., 3.])
>>> b.dtype
dtype('float64')

Les différents types de données existants nous permettent de mieux gérer l'occupation mémoire des données, mais la plupart du temps, nous travaillerons simplement avec des nombres à virgule flottante. Veuillez noter que dans l'exemple ci-dessus, NumPy détecte automatiquement le type de données en entrée.

Vous pouvez explicitement spécifier le type à utiliser :

>>> c = np.array([1, 2, 3], dtype=float)
>>> c.dtype
dtype('float64')

Par défaut, le type utilisé est le float :

>>> a = np.ones((3, 3))
>>> a.dtype
dtype('float64')

Mais il existe d'autres types :

>>> d = np.array([1+2j, 3+4j, 5+6*1j])
>>> d.dtype
dtype('complex128')

>>> e = np.array([True, False, False, True])
>>> e.dtype
dtype('bool')

>>> f = np.array(['Bonjour', 'Hello', 'Hallo',])
>>> f.dtype     # <--- strings containing max. 7 letters
dtype('S7')

Maintenant que nous avons nos premiers tableaux de données, nous allons visualiser leur contenu.

Commençons par lancer IPython en mode pylab.

$ ipython --pylab

Ou le notebook :

$ ipython notebook --pylab=inline

Sinon, si IPython est déjà lancé :

>>> %pylab

Ou, depuis le notebook :

>>> %pylab inline

Le inline est important dans le notebook, car il permet l'affichage dans le notebook même et non dans une nouvelle fenêtre.

Matplotlib est un outil de tracés de graphiques 2D. Nous pouvons importer ses fonctions ainsi :

>>> import matplotlib.pyplot as plt  # the tidy way

Puis on les utilise comme suit (notez que vous devez utiliser show de manière explicite) :

>>> plt.plot(x, y)       # line plot    
>>> plt.show()           # <-- shows the plot (not needed with pylab)

Ou encore, si vous utilisez pylab :

>>> plot(x, y)       # line plot

Utiliser import matplotlib.pyplot as plt est recommandé dans les scripts. Tandis que pylab est plus recommandé pour le travail interactif.

• Tracé de courbe 1D

>>> x = np.linspace(0, 3, 20)
>>> y = np.linspace(0, 9, 20)
>>> plt.plot(x, y)       # line plot    
[<matplotlib.lines.Line2D object at ...>]
>>> plt.plot(x, y, 'o')  # dot plot    
[<matplotlib.lines.Line2D object at ...>]

• Tableaux 2D (par exemple, les pixels d'une image) :

>>> image = np.random.rand(30, 30)
>>> plt.imshow(image, cmap=plt.cm.hot)    
>>> plt.colorbar()    
<matplotlib.colorbar.Colorbar instance at ...>

Les éléments d'un tableau peuvent être accédés et assignés de la même façon que les listes en Python :

>>> a = np.arange(10)
>>> a
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> a[0], a[2], a[-1]
(0, 2, 9)

L'idiome usuel de Python pour lire une séquence à l'envers est également supporté :

>>> a[::-1]
array([9, 8, 7, 6, 5, 4, 3, 2, 1, 0])

Pour des tableaux multidimensionnels, les index sont des tuples d'entiers :

>>> a = np.diag(np.arange(3))
>>> a
array([[0, 0, 0],
       [0, 1, 0],
       [0, 0, 2]])
>>> a[1, 1]
1
>>> a[2, 1] = 10 # third line, second column
>>> a
array([[ 0,  0,  0],
       [ 0,  1,  0],
       [ 0, 10,  2]])
>>> a[1]
array([0, 1, 0])

Échantillonnage (ou slicing) : les tableaux, comme toute séquence Python, peuvent être échantillonnés :

>>> a = np.arange(10)
>>> a
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> a[2:9:3] # [start:end:step]
array([2, 5, 8])

Notez que le dernier index n'est pas inclus :

>>> a[:4]
array([0, 1, 2, 3])

Les trois arguments start, end et step du slicing ne sont pas forcément tous requis : par défaut, l'index de départ start vaut 0, l'index d'arrivée end est celui du dernier élément et le pas step vaut 1 :

>>> a[1:3]
array([1, 2])
>>> a[::2]
array([0, 2, 4, 6, 8])
>>> a[3:]
array([3, 4, 5, 6, 7, 8, 9])

La petite illustration suivante résume bien la façon dont NumPy gère l'indexation et l'échantillonnage :

Vous pouvez également combiner les deux techniques :

>>> a = np.arange(10)
>>> a[5:] = 10
>>> a
array([ 0,  1,  2,  3,  4, 10, 10, 10, 10, 10])
>>> b = np.arange(5)
>>> a[5:] = b[::-1]
>>> a
array([0, 1, 2, 3, 4, 4, 3, 2, 1, 0])

>>> np.arange(6) + np.arange(0, 51, 10)[:, np.newaxis]
array([[ 0,  1,  2,  3,  4,  5],
       [10, 11, 12, 13, 14, 15],
       [20, 21, 22, 23, 24, 25],
       [30, 31, 32, 33, 34, 35],
       [40, 41, 42, 43, 44, 45],
       [50, 51, 52, 53, 54, 55]])

[[1, 1, 1, 1],
 [1, 1, 1, 1],
 [1, 1, 1, 2],
 [1, 6, 1, 1]]

[[0., 0., 0., 0., 0.],
 [2., 0., 0., 0., 0.],
 [0., 3., 0., 0., 0.],
 [0., 0., 4., 0., 0.],
 [0., 0., 0., 5., 0.],
 [0., 0., 0., 0., 6.]]

[[4, 3, 4, 3, 4, 3],
 [2, 1, 2, 1, 2, 1],
 [4, 3, 4, 3, 4, 3],
 [2, 1, 2, 1, 2, 1]]

Toute opération d'échantillonnage crée une vue du tableau original, qui consiste juste en une façon différente d'accéder aux données du tableau. Aussi, ce tableau n'est pas recopié en mémoire. Vous pouvez utiliser np.may_share_memory() pour vérifier si deux tableaux partagent le même espace mémoire. Notez également que cela utilise des heuristiques et peut vous renvoyer de faux positifs.

Quand vous modifiez une vue, le tableau originel est également modifié :

>>> a = np.arange(10)
>>> a
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> b = a[::2]
>>> b
array([0, 2, 4, 6, 8])
>>> np.may_share_memory(a, b)
True
>>> b[0] = 12
>>> b
array([12,  2,  4,  6,  8])
>>> a   # (!)
array([12,  1,  2,  3,  4,  5,  6,  7,  8,  9])

>>> a = np.arange(10)
>>> c = a[::2].copy()  # force a copy
>>> c[0] = 12
>>> a
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

>>> np.may_share_memory(a, c)
False

Cela peut paraître surprenant de prime abord, mais cela permet de gagner en mémoire et en temps de traitement.

Exemple : filtre pour nombres premiers

Calculer les nombres premiers entre 0 et 99, à l'aide d'un masque.

• Construisez un tableau booléen is_prime de type (100,), rempli avec True

 

Sélectionnez

>>> is_prime = np.ones((100,), dtype=bool)

• Annotez à False 0 et 1 qui ne sont pas des nombres premiers :

 

Sélectionnez

>>> is_prime[:2] = 0

• Pour chaque entier j à partir de 2, mettez à False tous ses multiples :

 

Sélectionnez

>>> N_max = int(np.sqrt(len(is_prime)))
>>> for j in range(2, N_max):
...     is_prime[2*j::j] = False

• À l'aide de help(no.nonzero), affichez les nombres premiers.

• Effectuez les tâches suivantes :

  • placez les codes précédents dans un module et nommez-le prime_sieve.py ;
  • lancez-le pour vous assurer qu'il fonctionne ;

  • implémentez la version optimisée du Crible d'Eratosthèneshttp://en.wikipedia.org/wiki/Sieve_of_Eratosthenes

    • ne tenez pas compte de j qui n'est pas un nombre premier,
    • le premier nombre à supprimer est j².

Les tableaux Numpy peuvent être indexés grâce à l'échantillonnage, mais aussi avec des tableaux de booléens ou d'entiers (masques). Cette méthode est appelée : indexation avancée (fancy indexing). Cela crée des copies et non des vues.

I-G-1. Utilisation de masque booléen▲

>>> np.random.seed(3)
>>> a = np.random.random_integers(0, 20, 15)
>>> a
array([10,  3,  8,  0, 19, 10, 11,  9, 10,  6,  0, 20, 12,  7, 14])
>>> (a % 3 == 0)
array([False,  True, False,  True, False, False, False,  True, False,
        True,  True, False,  True, False, False], dtype=bool)
>>> mask = (a % 3 == 0)
>>> extract_from_a = a[mask] # or,  a[a%3==0]
>>> extract_from_a           # extract a sub-array with the mask
array([ 3,  0,  9,  6,  0, 12])

>>> a[a % 3 == 0] = -1
>>> a
array([10, -1,  8, -1, 19, 10, 11, -1, 10, -1, -1, 20, -1,  7, 14])

I-G-2. Indexation à l'aide d'un tableau d'entiers▲

 

Sélectionnez

>>> a = np.arange(0, 100, 10)
>>> a
array([ 0, 10, 20, 30, 40, 50, 60, 70, 80, 90])

L'indexation peut être effectuée avec des tableaux d'entiers, où le même index est répété à plusieurs reprises :

 

Sélectionnez

>>> a[[2, 3, 2, 4, 2]]  # note: [2, 3, 2, 4, 2] is a Python list
array([20, 30, 20, 40, 20])

De nouvelles valeurs peuvent être assignées comme suit :

 

Sélectionnez

>>> a[[9, 7]] = -100
>>> a
array([   0,   10,   20,   30,   40,   50,   60, -100,   80, -100])

Quand un nouveau tableau est créé en indexant par l'intermédiaire d'un tableau d'entiers, le nouveau tableau a la même forme que le tableau d'origine :

 

Sélectionnez

>>> a = np.arange(10)
>>> idx = np.array([[3, 4], [9, 7]])
>>> idx.shape
(2, 2)
>>> a[idx]
array([[3, 4],
       [9, 7]])

L'image ci-dessous montre diverses utilisations d'indexations avancées :

Exercice : indexation avancée

Une fois de plus, reproduisez l'indexation avancée du diagramme ci-dessus.

Utilisez l'indexation avancée d'un côté et la création de tableaux de l'autre pour assigner les valeurs d'un tableau, par exemple en modifiant les éléments du tableau du diagramme qui sont en dessous de zéro.

II-A-1. Opérations basiques▲

>>> a = np.array([1, 2, 3, 4])
>>> a + 1
array([2, 3, 4, 5])
>>> 2**a
array([ 2,  4,  8, 16])

>>> b = np.ones(4) + 1
>>> a - b
array([-1.,  0.,  1.,  2.])
>>> a * b
array([ 2.,  4.,  6.,  8.])

>>> j = np.arange(5)
>>> 2**(j + 1) - j
array([ 2,  3,  6, 13, 28])

>>> a = np.arange(10000)
>>> %timeit a + 1  
10000 loops, best of 3: 24.3 us per loop
>>> l = range(10000)
>>> %timeit [i+1 for i in l] 
1000 loops, best of 3: 861 us per loop

>>> c = np.ones((3, 3))
>>> c * c                   # NOT matrix multiplication!
array([[ 1.,  1.,  1.],
       [ 1.,  1.,  1.],
       [ 1.,  1.,  1.]])

>>> c.dot(c)
array([[ 3.,  3.,  3.],
       [ 3.,  3.,  3.],
       [ 3.,  3.,  3.]])

II-A-2. Opérations complémentaires▲

>>> a = np.array([1, 2, 3, 4])
>>> b = np.array([4, 2, 2, 4])
>>> a == b
array([False,  True, False,  True], dtype=bool)
>>> a > b
array([False, False,  True, False], dtype=bool)

>>> a = np.array([1, 2, 3, 4])
>>> b = np.array([4, 2, 2, 4])
>>> c = np.array([1, 2, 3, 4])
>>> np.array_equal(a, b)
False
>>> np.array_equal(a, c)
True

>>> a = np.array([1, 1, 0, 0], dtype=bool)
>>> b = np.array([1, 0, 1, 0], dtype=bool)
>>> np.logical_or(a, b)
array([ True,  True,  True, False], dtype=bool)
>>> np.logical_and(a, b)
array([ True, False, False, False], dtype=bool)

>>> a = np.arange(5)
>>> np.sin(a)
array([ 0.        ,  0.84147098,  0.90929743,  0.14112001, -0.7568025])
>>> np.log(a)
array([       -inf,  0.        ,  0.69314718,  1.09861229,  1.38629436])
>>> np.exp(a)
array([  1.00000000e+00,   2.71828183e+00,   7.38905610e+00,
         2.00855369e+01,   5.45981500e+01])

>>> a = np.arange(4)
>>> a + np.array([1, 2])  
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: operands could not be broadcast together with shapes (4) (2)

>>> a = np.triu(np.ones((3, 3)), 1)   # see help(np.triu)
>>> a
array([[ 0.,  1.,  1.],
       [ 0.,  0.,  1.],
       [ 0.,  0.,  0.]])
>>> a.T
array([[ 0.,  0.,  0.],
       [ 1.,  0.,  0.],
       [ 1.,  1.,  0.]])

>>> a += a.T

II-B-1. Calculer des sommes▲

 

Sélectionnez

>>> x = np.array([1, 2, 3, 4])
>>> np.sum(x)
10
>>> x.sum()
10

Somme de lignes et de colonnes :

 

Sélectionnez

>>> x = np.array([[1, 1], [2, 2]])
>>> x
array([[1, 1],
       [2, 2]])
>>> x.sum(axis=0)   # columns (first dimension)
array([3, 3])
>>> x[:, 0].sum(), x[:, 1].sum()
(3, 3)
>>> x.sum(axis=1)   # rows (second dimension)
array([2, 4])
>>> x[0, :].sum(), x[1, :].sum()
(2, 4)

Dans le même genre, pour des tableaux de dimensions plus élevées :

 

Sélectionnez

>>> x = np.random.rand(2, 2, 2)
>>> x.sum(axis=2)[0, 1]     
1.14764...
>>> x[0, 1, :].sum()     
1.14764...

II-B-2. Autres opérations de réduction▲

Extrema :

 

Sélectionnez

>>> x = np.array([1, 3, 2])
>>> x.min()
1
>>> x.max()
3

>>> x.argmin()  # index of minimum
0
>>> x.argmax()  # index of maximum
1

Opérations logiques :

 

Sélectionnez

>>> np.all([True, True, False])
False
>>> np.any([True, True, False])
True

Comparaison de tableaux :

 

Sélectionnez

>>> a = np.zeros((100, 100))
>>> np.any(a != 0)
False
>>> np.all(a == a)
True

>>> a = np.array([1, 2, 3, 2])
>>> b = np.array([2, 2, 3, 2])
>>> c = np.array([6, 4, 4, 5])
>>> ((a <= b) & (b <= c)).all()
True

Statistiques :

 

Sélectionnez

>>> x = np.array([1, 2, 3, 1])
>>> y = np.array([[1, 2, 3], [5, 6, 1]])
>>> x.mean()
1.75
>>> np.median(x)
1.5
>>> np.median(y, axis=-1) # last axis
array([ 2.,  5.])

>>> x.std()          # full population standard dev.
0.82915619758884995

Et beaucoup plus.

Exercice : réductions

• En dehors de sum, quelle autre fonction peut-on utiliser ?
• Quelle est la différence entre sum et cumsum ?

Exemple pratique : données statistiques

Les données contenues dans le fichier populations.txt décrivent les populations de lièvres et de lynx (et de carottes) dans le nord du Canada sur 20 ans.

Vous pouvez visualiser ces données dans un éditeur ou dans Ipython :

 

Sélectionnez

In [1]: !cat data/populations.txt

Pour commencer, chargez les données dans un tableau NumPy :

 

Sélectionnez

>>> data = np.loadtxt('data/populations.txt')
>>> year, hares, lynxes, carrots = data.T  # trick: columns to variables

Puis affichez-les :

 

Sélectionnez

>>> from matplotlib import pyplot as plt
>>> plt.axes([0.2, 0.1, 0.5, 0.8]) 
>>> plt.plot(year, hares, year, lynxes, year, carrots) 
>>> plt.legend(('Hare', 'Lynx', 'Carrot'), loc=(1.05, 0.5))

La moyenne des populations dans le temps est :

 

Sélectionnez

>>> populations = data[:, 1:]
>>> populations.mean(axis=0)
array([ 34080.95238095,  20166.66666667,  42400.        ])

La variation moyenne est :

 

Sélectionnez

>>> populations.std(axis=0)
array([ 20897.90645809,  16254.59153691,   3322.50622558])

Quelle espèce a la population la plus importante chaque année ?

 

Sélectionnez

>>> np.argmax(populations, axis=1)
array([2, 2, 0, 0, 1, 1, 2, 2, 2, 2, 2, 2, 0, 0, 0, 1, 2, 2, 2, 2, 2])

Exemple pratique : broadcasting avec algorithme de marche aléatoire

Considérons un simple algorithme de marche aléatoire unidimensionnel : à chaque pas de temps, nous avançons ou reculons, avec une probabilité similaire.

Nous allons ici tâcher de trouver la distance typique parcourue depuis le départ, après t sauts en avant ou en arrière.

Nous allons pour cela simuler plusieurs marcheurs pour trouver cette loi, et nous allons procéder à quelques simulations : nous allons créer un tableau 2D avec des histoires (chaque marcheur a sa propre histoire) avec sur une dimension le temps, et sur l'autre l'histoire.

 

Sélectionnez

>>> n_stories = 1000 # number of walkers
>>> t_max = 200      # time during which we follow the walker

Nous allons choisir de manière totalement aléatoire si nous avançons ou reculons :

 

Sélectionnez

>>> t = np.arange(t_max)
>>> steps = 2 * np.random.random_integers(0, 1, (n_stories, t_max)) - 1
>>> np.unique(steps) # Verification: all steps are 1 or -1
array([-1,  1])

Nous faisons avancer le marcheur en additionnant les pas dans le temps :

 

Sélectionnez

>>> positions = np.cumsum(steps, axis=1) # axis = 1: dimension of time
>>> sq_distance = positions**2

Nous obtenons alors le fil conducteur de l'histoire :

 

Sélectionnez

>>> mean_sq_distance = np.mean(sq_distance, axis=0)

On affiche le résultat :

 

Sélectionnez

>>> plt.figure(figsize=(4, 3)) 
<matplotlib.figure.Figure object at ...>
>>> plt.plot(t, np.sqrt(mean_sq_distance), 'g.', t, np.sqrt(t), 'y-') 
[<matplotlib.lines.Line2D object at ...>, <matplotlib.lines.Line2D object at ...>]
>>> plt.xlabel(r"$t$") 
<matplotlib.text.Text object at ...>
>>> plt.ylabel(r"$\sqrt{\langle (\delta x)^2 \rangle}$") 
<matplotlib.text.Text object at ...>

Nous venons d'établir un fait physique reconnu : la distance moyenne augmente proportionnellement avec la racine carrée du temps.

• Les opérations basiques entre tableaux NumPy se font élément par élément.
• Cela ne fonctionne que sur des tableaux de dimensions identiques.
  Dans tous les cas, il est possible d'effectuer les opérations entre des tableaux de tailles différentes si NumPy peut transformer ces tableaux de façon à ce qu'ils aient les mêmes dimensions. Cette conversion est alors appelée broadcasting.

L'image suivante donne un exemple de broadcasting :

Vérifions cela :

>>> a = np.tile(np.arange(0, 40, 10), (3, 1)).T
>>> a
array([[ 0,  0,  0],
       [10, 10, 10],
       [20, 20, 20],
       [30, 30, 30]])
>>> b = np.array([0, 1, 2])
>>> a + b
array([[ 0,  1,  2],
       [10, 11, 12],
       [20, 21, 22],
       [30, 31, 32]])

Nous avons déjà fait du broadcasting sans le savoir :

>>> a = np.ones((4, 5))
>>> a[0] = 2  # we assign an array of dimension 0 to an array of dimension 1
>>> a
array([[ 2.,  2.,  2.,  2.,  2.],
       [ 1.,  1.,  1.,  1.,  1.],
       [ 1.,  1.,  1.,  1.,  1.],
       [ 1.,  1.,  1.,  1.,  1.]])

Une petite astuce fort utile :

>>> a = np.arange(0, 40, 10)
>>> a.shape
(4,)
>>> a = a[:, np.newaxis]  # adds a new axis -> 2D array
>>> a.shape
(4, 1)
>>> a
array([[ 0],
       [10],
       [20],
       [30]])
>>> a + b
array([[ 0,  1,  2],
       [10, 11, 12],
       [20, 21, 22],
       [30, 31, 32]])

Le broadcasting semble être un peu magique, mais il est en réalité tout naturel quand nous devons résoudre un problème dont les données de sortie sont un tableau avec plus de dimensions qu'en entrée.

Exemple pratique : broadcasting

Nous allons construire un tableau des distances (en miles) entre chaque ville de la Route 66 : Chicago, Springfield, Saint-Louis, Tulsa, Oklahoma City, Amarillo, Santa Fe, Albuquerque, Flagstaff et Los Angeles.

 

Sélectionnez

>>> mileposts = np.array([0, 198, 303, 736, 871, 1175, 1475, 1544,
...        1913, 2448])
>>> distance_array = np.abs(mileposts - mileposts[:, np.newaxis])
>>> distance_array
array([[   0,  198,  303,  736,  871, 1175, 1475, 1544, 1913, 2448],
       [ 198,    0,  105,  538,  673,  977, 1277, 1346, 1715, 2250],
       [ 303,  105,    0,  433,  568,  872, 1172, 1241, 1610, 2145],
       [ 736,  538,  433,    0,  135,  439,  739,  808, 1177, 1712],
       [ 871,  673,  568,  135,    0,  304,  604,  673, 1042, 1577],
       [1175,  977,  872,  439,  304,    0,  300,  369,  738, 1273],
       [1475, 1277, 1172,  739,  604,  300,    0,   69,  438,  973],
       [1544, 1346, 1241,  808,  673,  369,   69,    0,  369,  904],
       [1913, 1715, 1610, 1177, 1042,  738,  438,  369,    0,  535],
       [2448, 2250, 2145, 1712, 1577, 1273,  973,  904,  535,    0]])

Beaucoup de problèmes de calcul sur des grilles ou des réseaux peuvent être résolus par broadcasting. Par exemple, si nous désirons calculer la distance depuis l'origine jusqu'à chacun des points d'une grille de 10 x 10, nous ferons :

>>> x, y = np.arange(5), np.arange(5)[:, np.newaxis]
>>> distance = np.sqrt(x ** 2 + y ** 2)
>>> distance
array([[ 0.        ,  1.        ,  2.        ,  3.        ,  4.        ],
       [ 1.        ,  1.41421356,  2.23606798,  3.16227766,  4.12310563],
       [ 2.        ,  2.23606798,  2.82842712,  3.60555128,  4.47213595],
       [ 3.        ,  3.16227766,  3.60555128,  4.24264069,  5.        ],
       [ 4.        ,  4.12310563,  4.47213595,  5.        ,  5.65685425]])

Ou en couleur :

>>> plt.pcolor(distance)    
>>> plt.colorbar()

Remarque : la fonction numpy.ogrid permet de créer directement les vecteurs x et y de l'exemple précédent, avec deux « dimensions significatives ».

>>> x, y = np.ogrid[0:5, 0:5]
>>> x, y
(array([[0],
       [1],
       [2],
       [3],
       [4]]), array([[0, 1, 2, 3, 4]]))
>>> x.shape, y.shape
((5, 1), (1, 5))
>>> distance = np.sqrt(x ** 2 + y ** 2)

np.ogrid est très pratique dès que l'on désire effectuer des calculs sur une grille. D'un autre côté, np.mgrid fournit des matrices préremplies pour les cas où nous ne pouvons (où ne voulons) utiliser le broadcasting.

>>> x, y = np.mgrid[0:4, 0:4]
>>> x
array([[0, 0, 0, 0],
       [1, 1, 1, 1],
       [2, 2, 2, 2],
       [3, 3, 3, 3]])
>>> y
array([[0, 1, 2, 3],
       [0, 1, 2, 3],
       [0, 1, 2, 3],
       [0, 1, 2, 3]])

II-D-1. Aplatissement▲

>>> a = np.array([[1, 2, 3], [4, 5, 6]])
>>> a.ravel()
array([1, 2, 3, 4, 5, 6])
>>> a.T
array([[1, 4],
       [2, 5],
       [3, 6]])
>>> a.T.ravel()
array([1, 4, 2, 5, 3, 6])

II-D-2. Redimensionnement▲

>>> a.shape
(2, 3)
>>> b = a.ravel()
>>> b = b.reshape((2, 3))
>>> b
array([[1, 2, 3],
       [4, 5, 6]])

>>> a.reshape((2, -1))    # unspecified (-1) value is inferred
array([[1, 2, 3],
       [4, 5, 6]])

>>> b[0, 0] = 99
>>> a
array([[99,  2,  3],
       [ 4,  5,  6]])

>>> a = np.zeros((3, 2))
>>> b = a.T.reshape(3*2)
>>> b[0] = 9
>>> a
array([[ 0.,  0.],
       [ 0.,  0.],
       [ 0.,  0.]])

II-D-3. Ajout d'une dimension▲

>>> z = np.array([1, 2, 3])
>>> z
array([1, 2, 3])

>>> z[:, np.newaxis]
array([[1],
       [2],
       [3]])

>>> z[np.newaxis, :]
array([[1, 2, 3]])

II-D-4. Réarrangement▲

>>> a = np.arange(4*3*2).reshape(4, 3, 2)
>>> a.shape
(4, 3, 2)
>>> a[0, 2, 1]
5
>>> b = a.transpose(1, 2, 0)
>>> b.shape
(3, 2, 4)
>>> b[2, 1, 0]
5

>>> b[2, 1, 0] = -1
>>> a[0, 2, 1]
-1

II-D-5. Redimensionnement avec resize▲

>>> a = np.arange(4)
>>> a.resize((8,))
>>> a
array([0, 1, 2, 3, 0, 0, 0, 0])

>>> b = a
>>> a.resize((4,))   
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: cannot resize an array that has been referenced or is
referencing another array in this way.  Use the resize function

Tri le long d'un axe :

>>> a = np.array([[4, 3, 5], [1, 2, 1]])
>>> b = np.sort(a, axis=1)
>>> b
array([[3, 4, 5],
       [1, 1, 2]])

Tri in situ :

>>> a.sort(axis=1)
>>> a
array([[3, 4, 5],
       [1, 1, 2]])

Tri par indexation avancée (fancy indexing) :

>>> a = np.array([4, 3, 1, 2])
>>> j = np.argsort(a)
>>> j
array([2, 3, 1, 0])
>>> a[j]
array([1, 2, 3, 4])

Trouver le minimum et le maximum :

>>> a = np.array([4, 3, 1, 2])
>>> j_max = np.argmax(a)
>>> j_min = np.argmin(a)
>>> j_max, j_min
(0, 2)

II-F-1. Qu'avez-vous besoin de savoir pour débuter ?▲

>>> a[a < 0] = 0

III-A-1. Typage▲

>>> np.array([1, 2, 3]) + 1.5
array([ 2.5,  3.5,  4.5])

>>> a = np.array([1, 2, 3])
>>> a.dtype
dtype('int64')
>>> a[0] = 1.9     # <-- float is truncated to integer
>>> a
array([1, 2, 3])

>>> a = np.array([1.7, 1.2, 1.6])
>>> b = a.astype(int)  # <-- truncates to integer
>>> b
array([1, 1, 1])

>>> a = np.array([1.2, 1.5, 1.6, 2.5, 3.5, 4.5])
>>> b = np.around(a)
>>> b                    # still floating-point
array([ 1.,  2.,  2.,  2.,  4.,  4.])
>>> c = np.around(a).astype(int)
>>> c
array([1, 2, 2, 2, 4, 4])

III-A-2. Les différentes tailles de données▲

>>> np.array([1], dtype=int).dtype
dtype('int64')
>>> np.iinfo(np.int32).max, 2**31 - 1
(2147483647, 2147483647)
>>> np.iinfo(np.int64).max, 2**63 - 1
(9223372036854775807, 9223372036854775807L)

>>> np.iinfo(np.uint32).max, 2**32 - 1
(4294967295, 4294967295)
>>> np.iinfo(np.uint64).max, 2**64 - 1
(18446744073709551615L, 18446744073709551615L)

>>> np.finfo(np.float32).eps
1.1920929e-07
>>> np.finfo(np.float64).eps
2.2204460492503131e-16

>>> np.float32(1e-8) + np.float32(1) == 1
True
>>> np.float64(1e-8) + np.float64(1) == 1
False

>>> samples = np.zeros((6,), dtype=[('sensor_code', 'S4'),
...                                 ('position', float), ('value', float)])
>>> samples.ndim
1
>>> samples.shape
(6,)
>>> samples.dtype.names
('sensor_code', 'position', 'value')

>>> samples[:] = [('ALFA',   1, 0.37), ('BETA', 1, 0.11), ('TAU', 1,   0.13),
...               ('ALFA', 1.5, 0.37), ('ALFA', 3, 0.11), ('TAU', 1.2, 0.13)]
>>> samples
array([('ALFA', 1.0, 0.37), ('BETA', 1.0, 0.11), ('TAU', 1.0, 0.13),
       ('ALFA', 1.5, 0.37), ('ALFA', 3.0, 0.11), ('TAU', 1.2, 0.13)],
      dtype=[('sensor_code', 'S4'), ('position', '<f8'), ('value', '<f8')])

Les accès aux différents champs s'effectuent par indexation avec le nom du champ :

>>> samples['sensor_code']
array(['ALFA', 'BETA', 'TAU', 'ALFA', 'ALFA', 'TAU'],
      dtype='|S4')
>>> samples['value']
array([ 0.37,  0.11,  0.13,  0.37,  0.11,  0.13])
>>> samples[0]
('ALFA', 1.0, 0.37)

>>> samples[0]['sensor_code'] = 'TAU'
>>> samples[0]
('TAU', 1.0, 0.37)

Multiples accès à la fois :

>>> samples[['position', 'value']]
array([(1.0, 0.37), (1.0, 0.11), (1.0, 0.13), (1.5, 0.37), (3.0, 0.11),
       (1.2, 0.13)],
      dtype=[('position', '<f8'), ('value', '<f8')])

L'indexation avancée fonctionne également :

>>> samples[samples['sensor_code'] == 'ALFA']
array([('ALFA', 1.5, 0.37), ('ALFA', 3.0, 0.11)],
      dtype=[('sensor_code', 'S4'), ('position', '<f8'), ('value', '<f8')])

Les masques fonctionnent pour tous les types :

>>> x = np.ma.array([1, 2, 3, 4], mask=[0, 1, 0, 1])
>>> x
masked_array(data = [1 -- 3 --],
             mask = [False  True False  True],
       fill_value = 999999)


>>> y = np.ma.array([1, 2, 3, 4], mask=[0, 1, 1, 1])
>>> x + y
masked_array(data = [2 -- -- --],
             mask = [False  True  True  True],
       fill_value = 999999)

Sur les fonctions usuelles :

>>> np.ma.sqrt([1, -1, 2, -2])
masked_array(data = [1.0 -- 1.41421356237 --],
             mask = [False  True False  True],
       fill_value = 1e+20)

III-C-1. Bonnes pratiques▲

NumPy est également capable d'effectuer des calculs sur des polynômes :

Par exemple kitxmlcodeinlinelatexdvp3x^2+2x-1finkitxmlcodeinlinelatexdvp :

>>> p = np.poly1d([3, 2, -1])
>>> p(0)
-1
>>> p.roots
array([-1.        ,  0.33333333])
>>> p.order
2

>>> x = np.linspace(0, 1, 20)
>>> y = np.cos(x) + 0.3*np.random.rand(20)
>>> p = np.poly1d(np.polyfit(x, y, 3))

>>> t = np.linspace(0, 1, 200)
>>> plt.plot(x, y, 'o', t, p(t), '-')   
[<matplotlib.lines.Line2D object at ...>, <matplotlib.lines.Line2D object at ...>]

IV-A-1. Allons un peu plus loin▲

NumPy est également capable de calcul polynomial plus complexe, basé par exemple sur les polynômes de Chebychev. kitxmlcodeinlinelatexdvp3x^2+2x-1finkitxmlcodeinlinelatexdvp :

 

Sélectionnez

>>> p = np.polynomial.Polynomial([-1, 2, 3]) # coefs in different order!
>>> p(0)
-1.0
>>> p.roots()
array([-1.        ,  0.33333333])
>>> p.degree()  # In general polynomials do not always expose 'order'
2

Voici un exemple utilisant les polynômes de Chebyshev, sur l'intervalle [-1,1] :

 

Sélectionnez

>>> x = np.linspace(-1, 1, 2000)
>>> y = np.cos(x) + 0.3*np.random.rand(2000)
>>> p = np.polynomial.Chebyshev.fit(x, y, 90)

>>> t = np.linspace(-1, 1, 200)
>>> plt.plot(x, y, 'r.')   
[<matplotlib.lines.Line2D object at ...>]
>>> plt.plot(t, p(t), 'k-', lw=3)   
[<matplotlib.lines.Line2D object at ...>]

Les polynômes de Chebychev sont très intéressants pour les calculs d'interpolations.

IV-B-1. Fichiers texte▲

# year  hare    lynx    carrot
1900    30e3    4e3     48300
1901    47.2e3  6.1e3   48200
1902    70.2e3  9.8e3   41500
1903    77.4e3  35.2e3  38200

>>> data = np.loadtxt('data/populations.txt')
>>> data    
array([[  1900.,  30000.,   4000.,  48300.],
       [  1901.,  47200.,   6100.,  48200.],
       [  1902.,  70200.,   9800.,  41500.],
...

>>> np.savetxt('pop2.txt', data)
>>> data2 = np.loadtxt('pop2.txt')

In [1]: pwd      # show current directory
'/home/user/stuff/2011-numpy-tutorial'
In [2]: cd ex
'/home/user/stuff/2011-numpy-tutorial/ex'
In [3]: ls
populations.txt  species.txt

IV-B-2. Images▲

Avec Matplotlib :

 

Sélectionnez

>>> img = plt.imread('data/elephant.png')
>>> img.shape, img.dtype
((200, 300, 3), dtype('float32'))
>>> plt.imshow(img)     
<matplotlib.image.AxesImage object at ...>
>>> plt.savefig('plot.png')

>>> plt.imsave('red_elephant', img[:,:,0], cmap=plt.cm.gray)

En sauvegardant uniquement un canal (RGB) :

 

Sélectionnez

>>> plt.imshow(plt.imread('red_elephant.png'))  
<matplotlib.image.AxesImage object at ...>

Avec une autre bibliothèque :

 

Sélectionnez

>>> from scipy.misc import imsave
>>> imsave('tiny_elephant.png', img[::6,::6])
>>> plt.imshow(plt.imread('tiny_elephant.png'), interpolation='nearest')  
<matplotlib.image.AxesImage object at ...>

IV-B-3. Fichier au format NumPy▲

>>> data = np.ones((3, 3))
>>> np.save('pop.npy', data)
>>> data3 = np.load('pop.npy')

IV-B-4. Autres formats courants▲

1- Créez le tableau suivant (sans le saisir explicitement) et générez un nouveau tableau contenant sa 2^e et 4^e ligne.

[[1,  6, 11],
 [2,  7, 12],
 [3,  8, 13],
 [4,  9, 14],
 [5, 10, 15]]

2- Divisez chaque colonne du tableau a, élément par élément, avec le tableau b. (astuce : np.newaxis).

>>> a = np.arange(25).reshape(5, 5)
>>> b = np.array([1., 5, 10, 15, 20])

3- Un peu plus difficile : générez un tableau de 10 lignes de 3 colonnes avec des nombres aléatoires compris entre 0 et 1. Pour chaque ligne, extrayez le nombre le plus proche de 0,5.

• Utilisez abs et argsort pour trouver la colonne j la plus proche pour chaque ligne.
• Utilisez l'indexation avancée pour extraire les nombres. (astuce : a[i, j] ; le tableau i doit contenir le numéro de la colonne correspondante dans j).

Effectuons quelques manipulations sur les tableaux NumPy en commençant par la célèbre image de Lena (http://www.cs.cmu.edu/~chuck/lennapg/). Scipy fournit un tableau 2D de cette image avec la fonction scipy.lena :

>>> from scipy import misc
>>> lena = misc.lena()

Voici quelques images montrant ce que l'on peut obtenir avec quelques manipulations : utilisations de palettes de couleurs, rognage ou encore modification partielle de l'image.

• Commençons par utiliser la fonction imshow pour afficher l'image

   

  Sélectionnez

  In [3]: import pylab as plt
  In [4]: lena = misc.lena()
  In [5]: plt.imshow(lena)

• Lena est affichée en fausse couleur. Une palette de couleurs doit être précisée pour l'afficher en niveaux de gris.

   

  Sélectionnez

  In [6]: plt.imshow(lena, cmap=plt.cm.gray)

• Créez un tableau à partir de l'image : par exemple, retirez 30 pixels en bordure de l'image. Pour visualiser le résultat, utilisez à nouveau imshow.

   

  Sélectionnez

  In [9]: crop_lena = lena[30:-30,30:-30]

• Nous allons maintenant styliser le visage de Lena avec une marquise noire. Pour cela, nous avons besoin de créer un masque correspondant aux pixels que nous désirons colorer en noir. Ce masque est défini par la condition (y-256)**2 + (x-256)**2.

   

  Sélectionnez

  In [15]: y, x = np.ogrid[0:512,0:512] # x and y indices of pixels
  In [16]: y.shape, x.shape
  Out[16]: ((512, 1), (1, 512))
  In [17]: centerx, centery = (256, 256) # center of the image
  In [18]: mask = ((y - centery)**2 + (x - centerx)**2) > 230**2 # circle

• Ensuite, nous assignons la valeur 0 aux pixels de l'image correspondant à ce masque. La syntaxe est extrêmement simple et intuitive.

   

  Sélectionnez

  In [19]: lena[mask] = 0
  In [20]: plt.imshow(lena)
  Out[20]: <matplotlib.image.AxesImage object at 0xa36534c>

  

• Copiez toutes les instructions de cet exercice dans un script que vous nommerez lena_locket.py, puis exécutez-le.

• Remplacez le cercle par une ellipse.

Les données contenues dans le fichier populations.txt décrivent les populations de lièvres et de lynx (et de carottes) dans le nord du Canada sur 20 ans.

>>> data = np.loadtxt('data/populations.txt')
>>> year, hares, lynxes, carrots = data.T  # trick: columns to variables

>>> plt.axes([0.2, 0.1, 0.5, 0.8]) 
<matplotlib.axes.Axes object at ...>
>>> plt.plot(year, hares, year, lynxes, year, carrots) 
[<matplotlib.lines.Line2D object at ...>, ...]
>>> plt.legend(('Hare', 'Lynx', 'Carrot'), loc=(1.05, 0.5)) 
<matplotlib.legend.Legend object at ...>

Calculez et affichez à partir du jeu de données fourni :

• la moyenne de la population pour chaque espèce, pour chaque année de la période ;
• l'année pour laquelle chaque espèce a atteint son maximum de population ;
• quelle espèce a la population la plus importante pour chaque année (astuce : argsort et indexation avancée ) ;
• quelles années aucune des espèces citées ne dépasse les 50 000 (astuce : comparaisons et np.any) ;
• les deux années pour chaque espèce, où elles disposent de la population la plus basse (astuce : argsort, indexation avancée) ;
• comparez (graphique) l'évolution de la population de lièvres (voir help(np.gradient)) et de lynx. Vérifier la corrélation (voir herlp(np.corrcoef)).

…Le tout sans faire la moindre boucle for.

import numpy as np

data = np.loadtxt('../../../data/populations.txt')
year, hares, lynxes, carrots = data.T
populations = data[:,1:]

print "       Hares, Lynxes, Carrots"
print "Mean:", populations.mean(axis=0)
print "Std:", populations.std(axis=0)

j_max_years = np.argmax(populations, axis=0)
print "Max. year:", year[j_max_years]

max_species = np.argmax(populations, axis=1)
species = np.array(['Hare', 'Lynx', 'Carrot'])
print "Max species:"
print year
print species[max_species]

above_50000 = np.any(populations > 50000, axis=1)
print "Any above 50000:", year[above_50000]

j_top_2 = np.argsort(populations, axis=0)[:2]
print "Top 2 years with lowest populations for each:"
print year[j_top_2]

hare_grad = np.gradient(hares, 1.0)
print "diff(Hares) vs. Lynxes correlation", np.corrcoef(hare_grad, lynxes)[0,1]

import matplotlib.pyplot as plt
plt.plot(year, hare_grad, year, -lynxes)
plt.savefig('plot.png')

Écrivez une fonction f(a,b,c) qui retourne a^b-c. Créez un tableau de 24x12x6 contenant les valeurs de la fonction calculées dans le domaine [0,1]x[0,1]x[0,1].

Approchez la valeur de l'intégrale triple :

Le résultat exact est . Quelle est l'erreur relative ?

(Astuces : utilisez les opérations élément par élément et le broadcasting). Vous pouvez utiliser np.grid qui fournit un nombre de points dans un intervalle avec np.ogrid[0:1:20j]).

Rappelez-vous de la fonction Python :

def f(a, b, c):
    return some_result

import numpy as np
from numpy import newaxis

def f(a, b, c):
    return a**b - c

a = np.linspace(0, 1, 24)
b = np.linspace(0, 1, 12)
c = np.linspace(0, 1, 6)

samples = f(a[:,newaxis,newaxis],
            b[newaxis,:,newaxis],
            c[newaxis,newaxis,:])

# or,
#
# a, b, c = np.ogrid[0:1:24j, 0:1:12j, 0:1:6j]
# samples = f(a, b, c)

integral = samples.mean()

print "Approximation:", integral
print "Exact:", np.log(2) - 0.5

Écrivez un script qui calcule la fractale de Mandelbrot. L'itération de Mandelbrot est la suivante :

N_max = 50
some_threshold = 50

c = x + 1j*y

for j in xrange(N_max):
    z = z**2 + c

Les points (x, y) appartiennent à l'ensemble de Mandelbrot si |C| < some_threshold.

• Effectuez ce calcul en suivant les étapes :

  • construisez une grille de c=x+1j*y valeurs dans l'intervalle [-2,1]x[-1.5,1.5] ;
  • effectuez l'itération ;
  • créez le masque booléen 2D indiquant quels points appartiennent à l'ensemble ;
  • sauvegardez le résultat dans une image avec :

>>> import matplotlib.pyplot as plt
>>> plt.imshow(mask.T, extent=[-2, 1, -1.5, 1.5]) 
<matplotlib.image.AxesImage object at ...>
>>> plt.gray()
>>> plt.savefig('mandelbrot.png')

"""
Compute the Mandelbrot fractal
"""
import numpy as np
import matplotlib.pyplot as plt
from numpy import newaxis

def compute_mandelbrot(N_max, some_threshold, nx, ny):
    # A grid of c-values
    x = np.linspace(-2, 1, nx)
    y = np.linspace(-1.5, 1.5, ny)

    c = x[:,newaxis] + 1j*y[newaxis,:]

    # Mandelbrot iteration

    z = c
    for j in xrange(N_max):
        z = z**2 + c

    mandelbrot_set = (abs(z) < some_threshold)

    return mandelbrot_set

# Save

mandelbrot_set = compute_mandelbrot(50, 50., 601, 401)

plt.imshow(mandelbrot_set.T, extent=[-2, 1, -1.5, 1.5])
plt.gray()
plt.savefig('mandelbrot.png')

Soit P la matrice de transition d'une chaîne de Markov, et p la distribution des états :

1. 0 <= P[i,j] <= 1 : la probabilité de passer de l'état i à l'état j ;
2.  : la règle de transition ;
3. all(sum(P, axis=1) == 1), p.sum() == 1 : leur normalisation.

Écrivez un script pour une chaîne comportant cinq états qui :

• construit une matrice aléatoire et normalise chacune de ses lignes afin d'en faire une matrice de transition ;
• débute par une distribution (normalisée) aléatoire p et effectue 50 pas => p_50 ;
• calcule la distribution stationnaire : le vecteur propre de P.T associé à la valeur propre 1 (numériquement proche de 1) => p_stationary.

N'oubliez pas de normaliser le vecteur propre — Je ne l'ai pas fait…

• Vérifiez si p_50 et p_stationary sont égaux avec une tolérance de 1e^-5.

Boite à outils : np.random.rand, .dot(), np.linlg.eig, réductions, abs(), argmin, comparaisons, all, np.linalg.norm, etc.

import numpy as np

np.random.seed(1234)

n_states = 5
n_steps = 50
tolerance = 1e-5

# Random transition matrix and state vector
P = np.random.rand(n_states, n_states)
p = np.random.rand(n_states)

# Normalize rows in P
P /= P.sum(axis=1)[:,np.newaxis]

# Normalize p
p /= p.sum()

# Take steps
for k in xrange(n_steps):
    p = P.T.dot(p)

p_50 = p
print p_50

# Compute stationary state
w, v = np.linalg.eig(P.T)

j_stationary = np.argmin(abs(w - 1.0))
p_stationary = v[:,j_stationary].real
p_stationary /= p_stationary.sum()
print p_stationary

# Compare
if all(abs(p_50 - p_stationary) < tolerance):
    print "Tolerance satisfied in infty-norm"

if np.linalg.norm(p_50 - p_stationary) < tolerance:
    print "Tolerance satisfied in 2-norm"