.. _2020numpyrst:
==================================
Tech - calcul matriciel avec numpy
==================================
.. only:: html
**Links:** :download:`notebook <2020_numpy.ipynb>`, :downloadlink:`html <2020_numpy2html.html>`, :download:`python <2020_numpy.py>`, :downloadlink:`slides <2020_numpy.slides.html>`, :githublink:`GitHub|_doc/notebooks/td1a_home/2020_numpy.ipynb|*`
`numpy `__ est la librairie incontournable pour
faire des calculs en Python. Ces fonctionnalités sont disponibles dans
tous les langages et utilisent les optimisations processeurs. Il est
hautement improbable d’écrire un code aussi rapide sans l’utiliser.
`numpy `__ implémente ce qu’on appelle les
opérations matricielles basiques ou plus communément appelées
`BLAS `__.
Quelque soit le langage, l’implémentation est réalisée en langage bas
niveau (C, fortran, assembleur) et a été peaufinée depuis 50 ans au gré
des améliorations matérielles.
.. code:: ipython3
from jyquickhelper import add_notebook_menu
add_notebook_menu()
.. contents::
:local:
.. code:: ipython3
%matplotlib inline
Enoncé
------
La librairie `numpy `__ propose principalement deux
types :
`array `__
et
`matrix `__.
Pour faire simple, prenez toujours le premier. Ca évite les erreurs. Les
`array `__
sont des tableaux à plusieurs dimensions.
La maîtrise du slice
~~~~~~~~~~~~~~~~~~~~
Le slice est l’opérateur ``:`` (décrit sur la page
`indexing `__).
Il permet de récupérer une ligne, une colonne, un intervalle de valeurs.
.. code:: ipython3
import numpy
mat = numpy.array([[0, 5, 6, -3],
[6, 7, -4, 8],
[-5, 8, -4, 9]])
mat
.. parsed-literal::
array([[ 0, 5, 6, -3],
[ 6, 7, -4, 8],
[-5, 8, -4, 9]])
.. code:: ipython3
mat[:2], mat[:, :2], mat[0, 3], mat[0:2, 0:2]
.. parsed-literal::
(array([[ 0, 5, 6, -3],
[ 6, 7, -4, 8]]),
array([[ 0, 5],
[ 6, 7],
[-5, 8]]),
-3,
array([[0, 5],
[6, 7]]))
La maîtrise du nan
~~~~~~~~~~~~~~~~~~
`nan `__ est une
convention pour désigner une valeur manquante. Elle réagit de façon un
peut particulière. Elle n’est égale à aucune autre y compris elle-même.
.. code:: ipython3
numpy.nan == numpy.nan
.. parsed-literal::
False
.. code:: ipython3
numpy.nan == 4
.. parsed-literal::
False
Il faut donc utiliser une fonction spéciale
`isnan `__.
.. code:: ipython3
numpy.isnan(numpy.nan)
.. parsed-literal::
True
La maîtrise des types
~~~~~~~~~~~~~~~~~~~~~
Un tableau est défini par ses dimensions et le type unique des éléments
qu’il contient.
.. code:: ipython3
matint = numpy.array([0, 1, 2])
matint.shape, matint.dtype
.. parsed-literal::
((3,), dtype('int32'))
C’est le même type pour toute la matrice. Il existe plusieurs type
d’entiers et des réels pour des questions de performances.
.. code:: ipython3
%timeit matint * matint
.. parsed-literal::
508 ns ± 44.7 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
.. code:: ipython3
matintf = matint.astype(numpy.float64)
matintf.shape, matintf.dtype
.. parsed-literal::
((3,), dtype('float64'))
.. code:: ipython3
%timeit matintf * matintf
.. parsed-literal::
484 ns ± 30.4 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
.. code:: ipython3
%timeit matintf * matint
.. parsed-literal::
894 ns ± 78.8 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
Un changement de type et le calcul est plus long.
La maîtrise du broadcasting
~~~~~~~~~~~~~~~~~~~~~~~~~~~
Le
`broadcasting `__
signifie que certaines opérations ont un sens même si les dimensions des
tableaux ne sont pas tout à fait égales.
.. code:: ipython3
mat
.. parsed-literal::
array([[ 0, 5, 6, -3],
[ 6, 7, -4, 8],
[-5, 8, -4, 9]])
.. code:: ipython3
mat + 1000
.. parsed-literal::
array([[1000, 1005, 1006, 997],
[1006, 1007, 996, 1008],
[ 995, 1008, 996, 1009]])
.. code:: ipython3
mat + numpy.array([0, 10, 100, 1000])
.. parsed-literal::
array([[ 0, 15, 106, 997],
[ 6, 17, 96, 1008],
[ -5, 18, 96, 1009]])
.. code:: ipython3
mat + numpy.array([[0, 10, 100]]).T
.. parsed-literal::
array([[ 0, 5, 6, -3],
[ 16, 17, 6, 18],
[ 95, 108, 96, 109]])
La maîtrise des index
~~~~~~~~~~~~~~~~~~~~~
.. code:: ipython3
mat = numpy.array([[0, 5, 6, -3],
[6, 7, -4, 8],
[-5, 8, -4, 9]])
mat
.. parsed-literal::
array([[ 0, 5, 6, -3],
[ 6, 7, -4, 8],
[-5, 8, -4, 9]])
.. code:: ipython3
mat == 5
.. parsed-literal::
array([[False, True, False, False],
[False, False, False, False],
[False, False, False, False]])
.. code:: ipython3
mat == numpy.array([[0, -4, 9]]).T
.. parsed-literal::
array([[ True, False, False, False],
[False, False, True, False],
[False, False, False, True]])
.. code:: ipython3
(mat == numpy.array([[0, -4, 9]]).T).astype(numpy.int64)
.. parsed-literal::
array([[1, 0, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]], dtype=int64)
.. code:: ipython3
mat * (mat == numpy.array([[0, -4, 9]]).T).astype(numpy.int64)
.. parsed-literal::
array([[ 0, 0, 0, 0],
[ 0, 0, -4, 0],
[ 0, 0, 0, 9]], dtype=int64)
La maîtrise des fonctions
~~~~~~~~~~~~~~~~~~~~~~~~~
On peut regrouper les opérations que `numpy `__
propose en différents thèmes. Mais avant il
- L’\ **initialisation** :
`array `__,
`empty `__,
`zeros `__,
`ones `__,
`full `__,
`identity `__,
`rand `__,
`randn `__,
`randint `__
- Les **opérations basiques** : ``+``, ``-``, ``*``, ``/``, ``@``,
`dot `__
- Les **transformations** :
`transpose `__,
`hstack `__,
`vstack `__,
`reshape `__,
`squeeze `__,
`expend_dims `__
- Les **opérations de réduction** :
`minimum `__,
`maximum `__,
`argmin `__,
`argmax `__,
`sum `__,
`mean `__,
`prod `__,
`var `__,
`std `__
- Tout le reste comme la génération de matrices aléatoires, le calcul
des valeurs, vecteurs propres, des fonctions commme
`take `__,
…
Q1 : calculer la valeur du :math:`\chi_2` d’un tableau de contingence
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
La formule est
`là `__.
Et il faut le faire sans boucle. Vous pouvez comparer avec la fonction
`chisquare `__
de la librairie `scipy `__ qui est une extension
de `numpy `__.
.. math:: \chi_2 = N \sum_{i,j} p_{i.} p_{.j} \left( \frac{\frac{O_{ij}}{N} - p_{i.} p_{.j}}{p_{i.} p_{.j}}\right)^2
Q2 : calculer une distribution un peu particulière
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
La fonction
`histogram `__
permet de calculer la distribution empirique de variables. Pour cette
question, on tire un vecteur aléatoire de taille 10 avec la fonction
`rand `__,
on les trie par ordre croissant, on recommence plein de fois, on calcule
la distribution du plus grand nombre, du second plus grand nombre, …, du
plus petit nombre.
Q3 : on veut créer une matrice identité un million par un million
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Vous pouvez essayer sans réfléchir ou lire cette page d’abord :
`csr_matrix `__.
Q4 : vous devez créer l’application StopCovid
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Il existe une machine qui reçoit la position de 3 millions de téléphones
portable. On veut identifier les cas contacts (rapidement).
Réponses
--------
Q1 : calculer la valeur du :math:`\chi_2` d’un tableau de contingence
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
La formule est
`là `__.
Et il faut le faire sans boucle. Vous pouvez comparer avec la fonction
`chisquare `__
de la librairie `scipy `__ qui est une extension
de `numpy `__.
.. math:: \chi_2 = N \sum_{i,j} p_{i.} p_{.j} \left( \frac{\frac{O_{ij}}{N} - p_{i.} p_{.j}}{p_{i.} p_{.j}}\right)^2
.. code:: ipython3
import numpy
O = numpy.array([[15., 20., 13.], [4., 9., 5.]])
O
.. parsed-literal::
array([[15., 20., 13.],
[ 4., 9., 5.]])
.. code:: ipython3
def chi_square(O):
N = numpy.sum(O)
pis = numpy.sum(O, axis=1, keepdims=True) / N
pjs = numpy.sum(O, axis=0, keepdims=True) / N
pispjs = pis @ pjs
chi = pispjs * ((O / N - pispjs) / pispjs) ** 2
return numpy.sum(chi) * N
chi_square(O)
.. parsed-literal::
0.5798254016266716
Q2 : calculer une distribution un peu particulière
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
La fonction
`histogram `__
permet de calculer la distribution empirique de variables. Pour cette
question, on tire un vecteur aléatoire de taille 10 avec la fonction
`rand `__,
on les trie par ordre croissant, on recommence plein de fois, on calcule
la distribution du plus grand nombre, du second plus grand nombre, …, du
plus petit nombre.
.. code:: ipython3
rnd = numpy.random.rand(10)
rnd
.. parsed-literal::
array([0.98556467, 0.47377301, 0.77148185, 0.26135908, 0.27373018,
0.0240458 , 0.55360714, 0.3575369 , 0.71740732, 0.3260206 ])
.. code:: ipython3
numpy.sort(rnd)
.. parsed-literal::
array([0.0240458 , 0.26135908, 0.27373018, 0.3260206 , 0.3575369 ,
0.47377301, 0.55360714, 0.71740732, 0.77148185, 0.98556467])
.. code:: ipython3
def tirage(n):
rnd = numpy.random.rand(n)
trie = numpy.sort(rnd)
return trie[-1]
tirage(10)
.. parsed-literal::
0.876020129318981
.. code:: ipython3
def plusieurs_tirages(N, n):
rnd = numpy.random.rand(N, n)
return numpy.max(rnd, axis=1)
plusieurs_tirages(5, 10)
.. parsed-literal::
array([0.99594032, 0.67914189, 0.98105965, 0.93181536, 0.86827764])
.. code:: ipython3
t = plusieurs_tirages(5000, 10)
hist = numpy.histogram(t)
hist
.. parsed-literal::
(array([ 5, 8, 20, 35, 111, 221, 407, 785, 1273, 2135],
dtype=int64),
array([0.437878 , 0.49408914, 0.55030028, 0.60651142, 0.66272256,
0.7189337 , 0.77514485, 0.83135599, 0.88756713, 0.94377827,
0.99998941]))
.. code:: ipython3
import matplotlib.pyplot as plt
plt.plot(hist[1][1:], hist[0] / hist[0].sum());
.. image:: 2020_numpy_50_0.png
Q3 : on veut créer une matrice identité un million par un million
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Vous pouvez essayer sans réfléchir ou lire cette page d’abord :
`csr_matrix `__.
:math:`(10^6)^2=10^{12}`>10 Go, bref ça ne tient pas en mémoire sauf si
on a une grosse machine. Les matrices creuses (ou sparses en anglais),
sont adéquates pour représenter des matrices dont la grande majorité des
coefficients sont nuls car ceux-ci ne sont pas stockés. Concrètement, la
matrice enregistre uniquement les coordonnées des coefficients et les
valeurs non nuls.
.. code:: ipython3
import numpy
from scipy.sparse import csr_matrix
ide = csr_matrix((1000000, 1000000), dtype=numpy.float64)
ide.setdiag(1.)
.. parsed-literal::
c:\python395_x64\lib\site-packages\scipy\sparse\_index.py:125: SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient.
self._set_arrayXarray(i, j, x)
Q4 : vous devez créer l’application StopCovid
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Il existe une machine qui reçoit la position de 3 millions de téléphones
portable. On veut identifier les cas contacts (rapidement).
Si on devait calculer toutes les paires de distance, cela prendrait un
temps fou. Il faut ruser. Le plus simple est de construire une grille
sur le territoire français puis d’associer à chaque téléphone portable
la grille dans laquelle il se trouve. Dans une cellule de la grille, le
nombre de paires est beaucoup plus réduit. Ce n’est pas la seule astuce
qu’il faudra utiliser. Mais c’est un bon début.