Notes de cours Python scientifique

dtype décrit un simple objet dans un tableau :

>>> np.dtype(int).type
<type 'numpy.int64'>
>>> np.dtype(int).itemsize
8
>>> np.dtype(int).byteorder
'='

L'entête du fichier .wav :

• bloc de données brutes à 44 octets (au début du fichier)
• … suivi par  data_size d'octets de données réelles audio

L'entête du fichier .wav en tant que type de données structurées Numpy :

>>> wav_header_dtype = np.dtype([
...     ("chunk_id", (bytes, 4)),   # type scalaire à taille flexible, taille de l'objet 4
...     ("chunk_size", "<u4"),    # Nombre entier non signé à 32 bits au format little-endian
...     ("format", "S4"),         # chaîne de caractères à 4 octets
...     ("fmt_id", "S4"),
...     ("fmt_size", "<u4"),
...     ("audio_fmt", "<u2"),     #
...     ("num_channels", "<u2"),  # .. plus encore...
...     ("sample_rate", "<u4"),   #
...     ("byte_rate", "<u4"),
...     ("block_align", "<u2"),
...     ("bits_per_sample", "<u2"),
...     ("data_id", ("S1", (2, 2))), # sous-tableau, juste pour rire !
...     ("data_size", "u4"),
...     #
...     # Les données du son lui-même ne peuvent pas être représentées ici:
...     # ça n'a pas de taille fixe
...    ])

>>> wav_header_dtype['format']
dtype('|S4')
>>> wav_header_dtype.fields     
<dictproxy object at ...>
>>> wav_header_dtype.fields['format']
(dtype('|S4'), 8)

• Le premier élément est le « sub-dtype » dans les données structurées, correspondant au nom format
• Le deuxième élément est son décalage ( en octets) à partir du début de l'objet

Exercice

Faire un dtype « creux » en utilisant les décalages, et seulement quelques champs :

>>> wav_header_dtype = np.dtype(dict(
...   names=['format', 'sample_rate', 'data_id'],
...   offsets=[offset_1, offset_2, offset_3], # compté à partir du début de la structure en octets
...   formats=list of dtypes for each of the fields,
... ))

et l'utiliser pour lire la fréquence d'échantillonnage et data_id (en tant que sous-tableau).

>>> f = open('data/test.wav', 'r')
>>> wav_header = np.fromfile(f, dtype=wav_header_dtype, count=1)
>>> f.close()
>>> print(wav_header)
[ ('RIFF', 17402L, 'WAVE', 'fmt ', 16L, 1, 1, 16000L, 32000L, 2, 16, [['d', 'a'], ['t', 'a']], 17366L)]
>>> wav_header['sample_rate']
array([16000], dtype=uint32)

Essayons d'accéder au sous-tableau :

>>> wav_header['data_id']
array([[['d', 'a'],
        ['t', 'a']]],
      dtype='|S1')
>>> wav_header.shape
(1,)
>>> wav_header['data_id'].shape
(1, 2, 2)

Quand on accède aux sous-tableaux, les dimensions sont ajoutées à la fin!

Conversion de type 

• sur attribution
• sur la construction d'un tableau
• sur l'arithmétique
• etc.
• et manuellement: .astype(dtype)

données réinterprétées

• manuellement: .view(dtype)

• En résumé, voici ce qu'il faut retenir du cast en arithmétique:

  • seul le type (not value!) des opérandes nous intéresse
  • le plus grand type « certifié » en mesure de représenter l'ensemble est considéré
  • Les scalaires peuvent « perdre » face aux tableaux dans certaines situations

• Le cast copie les données en général:

   

  Sélectionnez

  1.
  2.
  3.
  4.
  5.
  6.
  7.
  8.
  9.
  10.
  11.
  12.
  13.
  14.
  

  >>> x = np.array([1, 2, 3, 4], dtype=np.float)
  >>> x
  array([ 1.,  2.,  3.,  4.])
  >>> y = x.astype(np.int8)
  >>> y
  array([1, 2, 3, 4], dtype=int8)
  >>> y + 1
  array([2, 3, 4, 5], dtype=int8)
  >>> y + 256
  array([257, 258, 259, 260], dtype=int16)
  >>> y + 256.0
  array([ 257.,  258.,  259.,  260.])
  >>> y + np.array([256], dtype=np.int32)
  array([257, 258, 259, 260], dtype=int32)
  

• Cast sur un objet setté : le « dtype » du tableau n'est pas modifié sur l'affectation d'un objet :

>>> y[:] = y + 1.5
>>> y
array([2, 3, 4, 5], dtype=int8)

II-C-4-a. Réinterprétation / Observation▲

>>> x = np.array([1, 2, 3, 4], dtype=np.uint8)
>>> x.dtype = "<i2"
>>> x
array([ 513, 1027], dtype=int16)
>>> 0x0201, 0x0403
(513, 1027)

>>> y = x.view("<i4")
>>> y
array([67305985], dtype=int32)
>>> 0x04030201
67305985

>>> x[1] = 5
>>> y
array([328193], dtype=int32)
>>> y.base is x
True

>>> x = np.zeros((10, 10, 4), dtype=np.int8)
>>> x[:, :, 0] = 1
>>> x[:, :, 1] = 2
>>> x[:, :, 2] = 3
>>> x[:, :, 3] = 4

>>> y = ...                     

>>> assert (y['r'] == 1).all()  
>>> assert (y['g'] == 2).all()  
>>> assert (y['b'] == 3).all()  
>>> assert (y['a'] == 4).all()

>>> y = x.view([('r', 'i1'),
... ('g', 'i1'),
... ('b', 'i1'),
... ('a', 'i1')]
... )[:, :, 0]

>>> y = np.array([[1, 3], [2, 4]], dtype=np.uint8).transpose()
>>> x = y.copy()
>>> x
array([[1, 2],
       [3, 4]], dtype=uint8)
>>> y
array([[1, 2],
       [3, 4]], dtype=uint8)
>>> x.view(np.int16)
array([[ 513],
       [1027]], dtype=int16)
>>> 0x0201, 0x0403
(513, 1027)
>>> y.view(np.int16)
array([[ 769, 1026]], dtype=int16)

>>> 0x0301, 0x0402
(769, 1026)

La question :

>>> x = np.array([[1, 2, 3],
...              [4, 5, 6],
...              [7, 8, 9]], dtype=np.int8)
>>> str(x.data)
'\x01\x02\x03\x04\x05\x06\x07\x08\t'

À quel octet débute l'objet x[1,2] dans x.data ?

La réponse (dans Numpy) :

• strides: le nombre d'octets à franchir pour trouver le prochain élément
• 1 stride par dimension

>>> x.strides
(3, 1)
>>> byte_offset = 3*1 + 1*2   # to find x[1,2]
>>> x.flat[byte_offset]
6
>>> x[1, 2]
6

• Simple, flexible

II-D-1-a. Programmation en C et Fortran▲

 

Sélectionnez

1.
2.
3.
4.
5.
6.

>>> x = np.array([[1, 2, 3],
...               [4, 5, 6]], dtype=np.int16, order='C')
>>> x.strides
(6, 2)
>>> str(x.data)
'\x01\x00\x02\x00\x03\x00\x04\x00\x05\x00\x06\x00'

• Nécessité d'augmenter de 6 octets pour trouver la ligne suivante
• Nécessité d'augmenter de 2 octets pour trouver la colonne suivante

 

Sélectionnez

1.
2.
3.
4.
5.

>>> y = np.array(x, order='F')
>>> y.strides
(2, 4)
>>> str(y.data)
'\x01\x00\x04\x00\x07\x00\x02\x00\x05\x00\x08\x00\x03\x00\x06\x00'

• Nécessité d'augmenter 2 octets pour trouver la ligne suivante
• Nécessité d'augmenter 4 octets pour trouver la colonne suivante
• De la même façon avec des dimensions supérieures
• C : les dernières dimensions varient plus vite (= strides plus petits)
• F : les premières dimensions varient plus vite

Maintenant nous pouvons comprendre le comportement de .view() :

 

Sélectionnez

>>> y = np.array([[1, 3], [2, 4]], dtype=np.uint8).transpose()
>>> x = y.copy()

La transposition n'affecte pas la mise en mémoire des données, seulement les strides

 

Sélectionnez

1.
2.
3.
4.
5.
6.
7.
8.
9.

>>> x.strides
(2, 1)
>>> y.strides
(1, 2)

>>> str(x.data)
'\x01\x02\x03\x04'
>>> str(y.data)
'\x01\x03\x02\x04'

• les résultats sont différents lorsque c'est interprété comme deux « int16 »
• .copy() crée de nouveaux tableaux comme en C (par défaut)

II-D-1-b. Partage avec des nombres entiers▲

>>> x = np.array([1, 2, 3, 4, 5, 6], dtype=np.int32)
>>> y = x[::-1]
>>> y
array([6, 5, 4, 3, 2, 1], dtype=int32)
>>> y.strides
(-4,)
>>> y = x[2:]
>>> y.__array_interface__['data'][0] - x.__array_interface__['data'][0]
8
>>> x = np.zeros((10, 10, 10), dtype=np.float)
>>> x.strides
(800, 80, 8)
>>> x[::2,::3,::4].strides
(1600, 240, 32)

>>> x = np.zeros((10, 10, 10), dtype=np.float)
>>> x.strides
(800, 80, 8)
>>> x.T.strides
(8, 80, 800)

>>> a = np.arange(6, dtype=np.int8).reshape(3, 2)
>>> b = a.T
>>> b.strides
(1, 2)

>>> str(a.data)  
'\x00\x01\x02\x03\x04\x05'
>>> b
array([[0, 2, 4],
       [1, 3, 5]], dtype=int8)
>>> c = b.reshape(3*2)
>>> c
array([0, 2, 4, 1, 3, 5], dtype=int8)

Manipulation d'un stride

>>> from numpy.lib.stride_tricks import as_strided
>>> help(as_strided)    
as_strided(x, shape=None, strides=None)
   Fare un ndarray à partir du tableau avec les paramètres shape et strides fournis

>>> x = np.array([1, 2, 3, 4], dtype=np.int16)
>>> as_strided(x, strides=(2*2, ), shape=(2, ))
array([1, 3], dtype=int16)
>>> x[::2]
array([1, 3], dtype=int16)

Exercice

array([1, 2, 3, 4], dtype=np.int8)
-> array([[1, 2, 3, 4],
          [1, 2, 3, 4],
          [1, 2, 3, 4]], dtype=np.int8)

en utilisant seulement as_strided :

Hint: byte_offset = stride[0]*index[0] + stride[1]*index[1] + ...

Solution

Stride peut aussi être 0

>>> x = np.array([1, 2, 3, 4], dtype=np.int8)
>>> y = as_strided(x, strides=(0, 1), shape=(3, 4))
>>> y
array([[1, 2, 3, 4],
       [1, 2, 3, 4],
       [1, 2, 3, 4]], dtype=int8)
>>> y.base.base is x
True

• Faisons quelque chose d'utile avec ça : produit externe de [1, 2, 3, 4] et [5, 6, 7]

>>> x = np.array([1, 2, 3, 4], dtype=np.int16)
>>> x2 = as_strided(x, strides=(0, 1*2), shape=(3, 4))
>>> x2
array([[1, 2, 3, 4],
       [1, 2, 3, 4],
       [1, 2, 3, 4]], dtype=int16)

>>> y = np.array([5, 6, 7], dtype=np.int16)
>>> y2 = as_strided(y, strides=(1*2, 0), shape=(3, 4))
>>> y2
array([[5, 5, 5, 5],
       [6, 6, 6, 6],
       [7, 7, 7, 7]], dtype=int16)

 >>> x2 * y2
array([[ 5, 10, 15, 20],
       [ 6, 12, 18, 24],
       [ 7, 14, 21, 28]], dtype=int16)

… semble quelque peu familier…

>>> x = np.array([1, 2, 3, 4], dtype=np.int16)
>>> y = np.array([5, 6, 7], dtype=np.int16)
>>> x[np.newaxis,:] * y[:,np.newaxis]
array([[ 5, 10, 15, 20],
       [ 6, 12, 18, 24],
       [ 7, 14, 21, 28]], dtype=int16)

• En interne, le principe de diffusion de tableau est en effet implémenté en utilisant des 0-strides.

Défi

• Extraire la diagonale principale de la matrice : (en supposant que la mémoire est ordonnée comme en C):

   

  Sélectionnez

  1.
  2.
  3.
  4.
  5.
  

  >>> x = np.array([[1, 2, 3],
  ...               [4, 5, 6],
  ...               [7, 8, 9]], dtype=np.int32)
  
  >>> x_diag = as_strided(x, shape=(3,), strides=(???,))
  

• Extraire les éléments de la première super-diagonale [2, 6].

• Et les sous-diagonales?

(Astuce pour les deux derniers points : le premier découpage déplace le point où le « striding » a commencé.)

Solution

Sélectionnez les diagonales:

>>> x_diag = as_strided(x, shape=(3, ), strides=((3+1)*x.itemsize, ))
>>> x_diag
array([1, 5, 9], dtype=int32)

Tranchez d'abord, pour ajuster le pointeur de données:

>>> as_strided(x[0, 1:], shape=(2, ), strides=((3+1)*x.itemsize, ))
array([2, 6], dtype=int32)

>>> as_strided(x[1:, 0], shape=(2, ), strides=((3+1)*x.itemsize, ))
array([4, 8], dtype=int32)

>>> y = np.diag(x, k=1)
>>> y
array([2, 6], dtype=int32)

>>> y.flags.owndata  
False

Défi

Calculer la trace du tenseur :

>>> x = np.arange(5*5*5*5).reshape(5,5,5,5)
>>> s = 0
>>> for i in xrange(5):
...    for j in xrange(5):
...       s += x[j,i,j,i]

en utilisant les strides et sum() sur le résultat.

>>> y = as_strided(x, shape=(5, 5), strides=(TODO, TODO))   
>>> s2 = ...   
>>> assert s == s2

Solution

>>> y = as_strided(x, shape=(5, 5), strides=((5*5*5 + 5)*x.itemsize,
... (5*5 + 1)*x.itemsize))
>>> s2 = y.sum()

La mise en mémoire peut affecter les performances:

In [1]: x = np.zeros((20000,))

In [2]: y = np.zeros((20000*67,))[::67]

In [3]: x.shape, y.shape
((20000,), (20000,))

In [4]: %timeit x.sum()
100000 loops, best of 3: 0.180 ms per loop

In [5]: %timeit y.sum()
100000 loops, best of 3: 2.34 ms per loop

In [6]: x.strides, y.strides
((8,), (536,))

Les strides plus petits sont-ils plus rapides ?

• Le CPU extrait des données de la mémoire principale à son cache en blocs

• Si de nombreux éléments du tableau sont ajustés en un seul bloc (petit stride):

  • moins de transferts nécessaires
  • plus rapide

• Parfois,

>>> a -= b

n'est pas pareil que :

>>> a -= b.copy()

>>> x = np.array([[1, 2], [3, 4]])
>>> x -= x.transpose()
>>> x
array([[ 0, -1],
       [ 4,  0]])

>>> y = np.array([[1, 2], [3, 4]])
>>> y -= y.T.copy()
>>> y
array([[ 0, -1],
       [ 1,  0]])

• x et x.transpose() partagent les données
• x -= x.transpose() modifie les données élément par élément
• parce que x and x.transpose() ont un striding différent, des données modifiées réapparaissent sur le RHS

1. Fournies par l'utilisateur

    

   Sélectionnez

   1.
   2.
   3.
   4.
   5.
   6.
   7.
   8.
   9.
   10.
   11.
   12.
   13.
   14.
   15.
   16.
   17.
   18.
   19.
   20.
   

   void ufunc_loop(void **args, int *dimensions, int *steps, void *data)
   {
       /*
        * int8 output = elementwise_function(int8 input_1, int8 input_2)
        *
   * Cette fonction doit calculer la ufunc pour plusieurs valeurs à la fois
        * de la manière suivante.
        */
       char *input_1 = (char*)args[0];
       char *input_2 = (char*)args[1];
       char *output = (char*)args[2];
       int i;
   
       for (i = 0; i < dimensions[0]; ++i) {
           *output = elementwise_function(*input_1, *input_2);
           input_1 += steps[0];
           input_2 += steps[1];
           output += steps[2];
       }
   }
   

2. La partie Numpy intégrée

    

   Sélectionnez

   1.
   2.
   3.
   4.
   5.
   6.
   7.
   8.
   9.
   10.
   11.
   12.
   13.
   14.
   15.
   16.
   17.
   

   char types[3]
   
   types[0] = NPY_BYTE   /* type du premier argument en entrée */
   types[1] = NPY_BYTE   /* type du second argument en entrée */
   types[2] = NPY_BYTE   /* type du troisième argument en entrée */
   
   PyObject *python_ufunc = PyUFunc_FromFuncAndData(
       ufunc_loop,
       NULL,
       types,
       1, /* ntypes */
       2, /* num_inputs */
       1, /* num_outputs */
       identity_element,
       name,
       docstring,
       unused)
   

   • Une ufunc peut également prendre en charge plusieurs combinaisons différentes de type entrée-sortie.

ufunc_loop est de forme très générique, et Numpy en fournit des préfabriquées

• Seul elementwise_func doit être fourni
• … sauf si votre fonction par élément n'est pas dans l'une des formes ci-dessus.

Vous pouvez modifier le masque par assignation :

>>> mx[1] = np.ma.masked
>>> mx
masked_array(data = [1 -- 3 -- 5],
             mask = [False  True False  True False],
       fill_value = 999999)

Le masque est ainsi effacé :

>>> mx[1] = 9
>>> mx
masked_array(data = [1 9 3 -- 5],
             mask = [False False False  True False],
       fill_value = 999999)

Le masque est également disponible directement :

>>> mx.mask
array([False, False, False,  True, False], dtype=bool)

Les entrées masquées peuvent être remplies avec une valeur déterminée pour obtenir un tableau en retour :

>>> x2 = mx.filled(-1)
>>> x2
array([ 1,  9,  3, -1,  5])

Le masque peut aussi être effacé

>>> mx.mask = np.ma.nomask
>>> mx
masked_array(data = [1 9 3 -99 5],
             mask = [False False False False False],
       fill_value = 999999)

Le package de tableau masqué contient également des fonctions de domaine-courant :

>>> np.ma.log(np.array([1, 2, -1, -2, 3, -5]))
masked_array(data = [0.0 0.69314718056 -- -- 1.09861228867 --],
             mask = [False False  True  True False  True],
       fill_value = 1e+20)

Exemple : Statistiques masquées

Des rangers canadiens ont été distraits lors du comptage de lièvres et de lynx en 1903-1910 et 1917-1918, et ont obtenu des chiffres incorrects. (Des cultivateurs de carottes sont restés en alerte, cependant.) Calculez les populations moyennes au fil du temps, en ignorant les chiffres invalides.

>>> data = np.loadtxt('data/populations.txt')
>>> populations = np.ma.masked_array(data[:,1:])
>>> year = data[:, 0]

>>> bad_years = (((year >= 1903) & (year <= 1910))
...            | ((year >= 1917) & (year <= 1918)))
>>> # '&' signifie 'and' et '|' signifie 'or'
>>> populations[bad_years, 0] = np.ma.masked
>>> populations[bad_years, 1] = np.ma.masked

>>> populations.mean(axis=0)
masked_array(data = [40472.7272727 18627.2727273 42400.0],
             mask = [False False False],
       fill_value = 1e+20)

>>> populations.std(axis=0)
masked_array(data = [21087.656489 15625.7998142 3322.50622558],
             mask = [False False False],
       fill_value = 1e+20)

Notez que Matplotlib connaît les tableaux masqués.

>>> plt.plot(year, populations, 'o-')   
[<matplotlib.lines.Line2D object at ...>, ...]

Titre : numpy.random.permutations échoue pour des arguments de type non entier

Je suis en train de générer des permutations aléatoires, en utilisant numpy.random.permutations

Lors de l'appel de numpy.random.permutation avec des arguments non entiers.

Cela échoue avec un message d'erreur sibyllin :

 >>> np.random.permutation(12)
    array([ 6, 11,  4, 10,  2,  8,  1,  7,  9,  3,  0,  5])
    >>> np.random.permutation(12.)
    Traceback (most recent call last):
      File "<stdin>", line 1, in <module>
      File "mtrand.pyx", line 3311, in mtrand.RandomState.permutation
      File "mtrand.pyx", line 3254, in mtrand.RandomState.shuffle
    TypeError: len() of unsized object

Cela se produit aussi avec des arguments de type long, et donc

np.random.permutation(X.shape[0]) où X est un tableau qui échoue sur Windows 64-bits(où la forme est un tuple de longs).

Ce serait génial si ça pouvait caster en entier ou au moins lever une erreur spécifique pour des types non entiers.

Je suis en train d'utiliser Numpy 1.4.1, construit à partir de l'archive officielle, sur Windows

64 avec Visual studio 2008, sur Python.org 64-bit Python.

0. Qu'êtes-vous en train de faire?

1. Extrait de code reproduisant l'erreur (si possible)

   • Que se passe-t-il réellement ?
   • Qu'est-ce que vous attendez ?
2. Plate-forme(Windows / Linux / OSX, 32/64 bits, x86/PPC…)

3. Version de Numpy/Scipy

    

   Sélectionnez

   >>> print np.__version__ 
   2...

4. Vérifiez que ce qui suit est ce que vous aviez prévu :

    

   Sélectionnez

   >>> print np.__file__ 
   /...

5. Dans le cas où vous traîniez de vieilles installations de Numpy ou endommagées.
   En cas de doute, essayez de supprimer les installations existantes de Numpy, puis réinstallez…