Vecteurs, matrices et tableaux
==============================

Les bases
------------------------------


Numpy ajoute le type ``array`` qui est similaire à une liste (``list``) avec
la condition supplémentaire que tous les éléments sont du même type. 
Nous concernant ce sera donc un tableau d'entiers, de flottants voire 
de booléens.

Une première méthode consiste à convertir une liste en un tableau via la commande
``array``. Le deuxième argument est optionnel et spécifie le type des éléments du tableau.

>>> a = np.array([1, 4, 5, 8], float)
>>> a
array([ 1., 4., 5., 8.])
>>> type(a)
<type 'numpy.ndarray'>


Un tableau peut être multidimensionnel ; ici dimension 3

.. code-block:: python

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

          [[1, 2],
           [1, 2]]])
   >>> type(a)
   <type 'numpy.ndarray'>
   >>> a[1,1,1]
   2

Nous limiterons notre étude aux tableaux {uni,bi}-dimensionnels.
En Python  modifier une donnée d'une extraction d'un tableau entraîne
aussi une modification du tableau initial ! Si nécessaire la fonction ``np.copy(a)`` ou ``a.copy()``
permet de faire une copie d'un tableau ``a``.

>>> a=np.array([1, 2, 3, 4, 5])
>>> c=np.array([1, 2, 3, 4, 5])
>>> b=a[1:3]
>>> b
array([2, 3])
>>> b[1]=0
>>> a
array([1, 2, 0, 4, 5])   # a  modifié
>>> b=c[1:3].copy()      # première méthode
>>> b[1]=0
>>> bb=np.copy(c[1:3])   # seconde méthode
>>> bb[1]=0
>>> c
array([1, 2, 3, 4, 5])   # c non modifié


Comme pour les listes le `slicing` extrait les tableaux

>>> a = np.arange(10)
>>> a
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> a[2:9:3] # [début:fin:pas]
array([2, 5, 8])
>>> a[2:8:3] # le dernier élément n'est pas inclus
array([2, 5])
>>> a[:5]    # le dernier élément n'est pas inclus
array([0, 1, 2, 3, 4])

Dans l'instruction ``[début:fin:pas]`` deux des arguments peuvent être omis : par défaut l'indice de 
début vaut 0 (le 1er élément), l'indice de fin est celui du dernier élément et le pas vaut 1. 
Il faut tout de même noter une différence entre ``[2:9]`` et ``[2:]``

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

Un ``pas`` négatif inversera l'ordre du tableau et les tableaux étant considérés cycliques
la commande ``a[-2::1]`` extrait l'avant dernier et dernier élément. Pour un tableau bi-dimensionnel
on peut bien sûr jouer avec les deux indices.

>>> a=np.eye(5,5)
>>> print a
[[ 1.  0.  0.  0.  0.]
 [ 0.  1.  0.  0.  0.]
 [ 0.  0.  1.  0.  0.]
 [ 0.  0.  0.  1.  0.]
 [ 0.  0.  0.  0.  1.]]
>>> a[:3,:4]=4
>>> a[::3,::4]=5
>>> a[:,3]=6
>>> print a
[[ 5.  4.  4.  6.  5.]
 [ 4.  4.  4.  6.  0.]
 [ 4.  4.  4.  6.  0.]
 [ 5.  0.  0.  6.  5.]
 [ 0.  0.  0.  6.  1.]]

Une autre méthode d'extraction consiste à écrire ``a[b]`` où ``b``
est un tableau d'entiers qui correspondra aux indices à extraire
et ``a`` un vecteur. Pour les matrices c'est ``a[b,c]`` et on prend les
indices de ligne dans ``b``, les indices de colonnes dans ``c``, successivement.

>>> a=np.array([2,4,6,8],float)
>>> b=np.array([0,0,1,3,2,1],int)
>>> a[b]
array([ 2.,  2.,  4.,  8.,  6.,  4.])
>>> a=a.reshape(2,2)
>>> b=np.array([0,0,1,1,0],int)
>>> c=np.array([0,1,1,1,1],int)
>>> a[b,c]
array([ 2.,  4.,  8.,  8.,  4.])
>>> 



Numpy/Scipy proposent aussi le type ``mat`` comme matrice, exclusivement un tableau bi-dimensionnel.
Comme les fonctions telles que ``ones``, ``eye`` retournent un objet de type ``array`` nous n'utiliserons pas dans la suite le type ``mat``. Il permet cependant de saisir les matrices à-la-Matlab et de
faire le produit matriciel par le simple symbole ``*``

>>> a=np.mat('[1 2 4 ; 3 4 5.]')
>>> b=np.mat('[2. ; 4 ; 6]')
>>> print a
[[ 1.  2.  4.]
 [ 3.  4.  5.]]
>>> print b
[[ 2.]
 [ 4.]
 [ 6.]]
>>> print a*b
[[ 34.]
 [ 52.]]

Tout objet ``array`` est convertible en type ``mat`` et réciproquement (sous la condition que le
tableau (``array``) soit {uni,bi}-dimensionnel)

>>> a=np.array([1, 2, 4])
>>> np.mat(a)
matrix([[1, 2, 4]])



Construire
----------


Nous avons déjà vu comment accéder aux éléments d'un tableau soit par ``a[i,j]`` soit par le `slicing`.
Pour définir certaines matrices (au sens ``array`` {uni,bi}-dimensionnel) il est possible
d'utiliser les boucles, des matrices prédéfinies de Numpy (Scipy) et les opérations élémentaires
(multiplication par un scalaire, produit matriciel, transposition, produit de Kronecker, 
jeu avec les booléens, etc).

Numpy propose ``arange`` équivalent de ``range`` mais de type ``array``, ce qui évite une conversion.

>>> np.arange(10)
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> np.arange(2, 10, dtype=np.float)
array([ 2., 3., 4., 5., 6., 7., 8., 9.])
>>> np.arange(2, 3, 0.1)
array([ 2. , 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9])
>>> np.arange(3)
array([0, 1, 2])
>>> np.arange(3.0)
array([ 0.,  1.,  2.])

Comme il y quelques subtilités dans la fonction ``range`` et ``arange`` quant au dernier élément, Pour
éviter tout désagrément la fonction ``linspace(premier,dernier,n)`` renvoie un ``array`` commençant 
par ``premier``, se terminant par ``dernier`` avec ``n`` éléments régulièrement espacés. 

>>> np.linspace(1., 4., 6)
array([ 1. ,  1.6,  2.2,  2.8,  3.4,  4. ])


Numpy propose le redimensionnement d'un tableau avec la fonction ``reshape``. Il faut tout de
même respecter une condition : le nombre d'éléments doit être le même !

>>> a=np.arange(16)
>>> a
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15])
>>> a.reshape(4,4)
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11],
       [12, 13, 14, 15]])
>>> a.reshape(2,8)
array([[ 0,  1,  2,  3,  4,  5,  6,  7],
       [ 8,  9, 10, 11, 12, 13, 14, 15]])


Il est possible aussi de transformer un ``array`` bi-dimensionnel en un ``array`` uni-dimensionnel 
avec ``reshape`` ou ``flatten``. L'instruction ``b.flatten()`` renvoie 
une copie de ``b``, ce qui n'est pas le cas de ``reshape``.
Ici un exemple avec la syntaxe ``np.reshape(array,dimension)``

>>> b=np.ones((4,4),float)+np.eye(4,4)
>>> print b
[[ 2.  1.  1.  1.]
 [ 1.  2.  1.  1.]
 [ 1.  1.  2.  1.]
 [ 1.  1.  1.  2.]]
>>> np.reshape(b,16)
array([ 2.,  1.,  1.,  1.,  1.,  2.,  1.,  1.,  1.,  1.,  2.,  1.,  1.,
        1.,  1.,  2.])
>>> np.reshape(b,(2,8))
array([[ 2.,  1.,  1.,  1.,  1.,  2.,  1.,  1.],
       [ 1.,  1.,  2.,  1.,  1.,  1.,  1.,  2.]])
>>> b.flatten()
array([ 2.,  1.,  1.,  1.,  1.,  2.,  1.,  1.,  1.,  1.,  2.,  1.,  1.,
        1.,  1.,  2.])

Numpy permet d'assembler les vecteurs et les matrices, de les concaténer. Par défaut
l'assemblage se fait selon la 1ère dimension (les lignes, donc assemblage vertical). L'option ``axis=1`` assemble "horizontalement".

>>> a=np.arange(4).reshape(2,2)
>>> b=4+np.arange(4).reshape(2,2)
>>> b
array([[4, 5],
       [6, 7]])
>>> c=(np.arange(4)+6)[newaxis,:]
>>> np.concatenate((a,b))
array([[0, 1],
       [2, 3],
       [4, 5],
       [6, 7]])
>>> np.concatenate((a,b),axis=0)
array([[0, 1],
       [2, 3],
       [4, 5],
       [6, 7]])
>>> np.concatenate((a,b),axis=1)
array([[0, 1, 4, 5],
       [2, 3, 6, 7]])
>>> np.concatenate((c,c,np.concatenate((a,b),axis=1)))
array([[6, 7, 8, 9],
       [6, 7, 8, 9],
       [0, 1, 4, 5],
       [2, 3, 6, 7]])

Pour concaténer un vecteur et une matrice 
il est nécessaire de convertir le vecteur en matrice à 1 ligne (ou 1 colonne).
Par rapport à Matlab/Scilab c'est bien dommage. Au passage la commande ``np.newaxis`` permet 
d'ajouter une dimension à un tableau (si ``b`` est un tableau unidimensionnel
``b[np.newaxis,:]`` retournera une matrice à 1 ligne tandis que ``b[:,np.newaxis]`` 
retournera une matrice à 1 colonne).

>>> c=(np.arange(4)+6)[np.newaxis,:]
>>> d=(np.arange(4)+6)
>>> np.concatenate((c,d))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: arrays must have same number of dimensions
>>> np.concatenate((c,d),axis=1)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: arrays must have same number of dimensions
>>> np.concatenate((c,d[np.newaxis,:]),axis=0)
array([[6, 7, 8, 9],
       [6, 7, 8, 9]])
>>> np.concatenate((c.transpose(),d[:,np.newaxis]),axis=0)
array([[6],
       [7],
       [8],
       [9],
       [6],
       [7],
       [8],
       [9]])


Le "slicing" permet d'extraire des éléments d'un vecteur ou matrice. Il existe aussi
la commande ``np.take`` qui extrait lignes/colonnes d'un ``array``. L'option ``axis`` détermine
une extraction selon les lignes/colonnes. La syntaxe se comprend aisément.

>>> a=np.arange(16).reshape(4,4)
>>> a.take([0,3],axis=1)
array([[ 0,  3],
       [ 4,  7],
       [ 8, 11],
       [12, 15]])
>>> a.take([0,3])
array([0, 3])
>>> a.take([0,3],axis=0)
array([[ 0,  1,  2,  3],
       [12, 13, 14, 15]])
>>> a.take([0,3,2,2,2],axis=0)
array([[ 0,  1,  2,  3],
       [12, 13, 14, 15],
       [ 8,  9, 10, 11],
       [ 8,  9, 10, 11],
       [ 8,  9, 10, 11]])

La réciproque de ``np.take`` est ``np.put``

>>> a=np.array([0, 1, 2, 3, 4, 5],float)
>>> b=np.array([9,8,7],float)
>>> a.put([0,3],b)
>>> a
array([ 9.,  1.,  2.,  8.,  4.,  5.])

Comme deux positions seulement ont été données, seules la valeurs 0 et 3 ont été
remplacées.


Matrices particulières
----------------------

Nous avons déjà utilisé certaines commandes



.. tabularcolumns:: |L|L|
   
.. table::

    ==================================  ==================================================================
      Commande                             Description (rapide)
    ==================================  ==================================================================
     ``np.zeros(n) ((n,p))``             vecteur nul   de taille ``n``, matrice nulle de taille ``n,p``
      ``np.eye(n) ((n,p))``              matrice de taille ``n`` (``n,p``) avec des 1 sur la diagonale
                                         et des zéros ailleurs
      ``np.ones(n) ((n,p))``             vecteur de taille ``n``, matrice de taille ``n,p`` remplie de 1
      ``np.diag(v)``                     matrice diagonale dont la diagonale est le vecteur ``v``
      ``np.diag(v,k)``                   matrice dont la 'diagonale' décalée de ``k`` est le vecteur ``v``
                                         (``k`` est un entier relatif)
     ``np.random.rand(n) ((n,p))``       vecteur (taille ``n``), matrice (taille ``n,p``) à coefficients
                                         aléatoires uniformes sur [0,1] 
    ==================================  ==================================================================

>>> v=np.ones(5)
>>> v
array([ 1.,  1.,  1.,  1.,  1.])
>>> np.diag(v,2)
array([[ 0.,  0.,  1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  1.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]])
>>> np.diag(v,-1)
array([[ 0.,  0.,  0.,  0.,  0.,  0.],
       [ 1.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  1.,  0.]])



Opérations
----------

L'opérateur ``+`` additionne terme à terme deux tableaux de même dimension. La commande ``2.+a`` renvoie
le vecteur/matrice dont tous les éléments sont ceux de ``a`` plus 2. Multplier un vecteur/matrice
par un scalaire se fait selon le même principe : la commande ``2*a`` renvoie le vecteur/matrice de même dimension dont tous les éléments ont été multipliés par 2.

L'opérateur ``*`` ne fait que multiplier terme à terme deux tableaux de même dimension:

>>> 4*np.ones((4,4))*np.diag([2, 3, 4., 5])
array([[  8.,   0.,   0.,   0.],
       [  0.,  12.,   0.,   0.],
       [  0.,   0.,  16.,   0.],
       [  0.,   0.,   0.,  20.]])

Dans le même ordre ``a**3`` renvoie le vecteur/matrice de même dimension dont tous les
éléments sont ceux de ``a`` élévés à la puissance 3. Devinez ce que fait ``1./a`` ?

Pour le produit matriciel (le vrai) la commande ``np.dot`` est là

>>> np.dot(4*np.ones((4,4)),np.diag([2, 3, 4., 5]))
array([[  8.,  12.,  16.,  20.],
       [  8.,  12.,  16.,  20.],
       [  8.,  12.,  16.,  20.],
       [  8.,  12.,  16.,  20.]])

Si ``v`` est un vecteur (``array`` uni-dimensionnel) et ``A`` une matrice (``array`` bi-dimensionnel)
alors ``np.dot(A,v)`` renvoie le produit :math:`Av` tandis que ``np.dot(v,A)`` renvoie :math:`v^tA`.
Si le produit n'est pas possible, Python vous le signale.


>>> a=np.arange(4).reshape(2,2)
>>> v=np.array([-3,2])
>>> np.dot(a,v)
array([2, 0])
>>> np.dot(v,a)
array([4, 3])
>>> w=np.concatenate((v,v))
>>> w
array([-3,  2, -3,  2])
>>> dot(a,w)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: objects are not aligned


``np.vdot(v,w)`` : produit scalaire des vecteurs ``v`` et ``w``.

Il y a enfin le produit de Kronecker, ``np.kron``, assez pratique pour créer des matrices/vecteurs
particuliers. Si :math:`A` et :math:`B` sont deux matrices non nécessairement de même taille, le produit
de Kronecker de :math:`A` :math:`B` est la matrice bloc (:math:`B` de taille :math:`n\times p`) 

.. math::
   \begin{pmatrix}
     b_{11}A & b_{12} A & \cdots & b_{1p} A \\
     \vdots &           &        & \vdots  \\
     b_{n1}A & b_{n2} A & \cdots & b_{np}A 
   \end{pmatrix}

>>> a=np.arange(4,dtype=float).reshape(2,2)
>>> b=np.arange(5,dtype=float)
>>> print a,b
[[ 0.  1.]
 [ 2.  3.]] [ 0.  1.  2.  3.  4.]
>>> np.kron(a,b)
array([[  0.,   0.,   0.,   0.,   0.,   0.,   1.,   2.,   3.,   4.],
       [  0.,   2.,   4.,   6.,   8.,   0.,   3.,   6.,   9.,  12.]])
>>> np.kron(b,a)
array([[  0.,   0.,   0.,   1.,   0.,   2.,   0.,   3.,   0.,   4.],
       [  0.,   0.,   2.,   3.,   4.,   6.,   6.,   9.,   8.,  12.]])

Un autre exemple qui peut être utilie pour certains exercices.

>>> a=np.arange(5,dtype=float)
>>> b=np.ones((1,5))
>>> np.kron(a,b.transpose())
array([[ 0.,  1.,  2.,  3.,  4.],
       [ 0.,  1.,  2.,  3.,  4.],
       [ 0.,  1.,  2.,  3.,  4.],
       [ 0.,  1.,  2.,  3.,  4.],
       [ 0.,  1.,  2.,  3.,  4.]])


Booléens
--------

On rappelle les opérateurs logiques en Python.

.. tabularcolumns:: |L|L|L|L|
   
.. table::

    =================  ======================= ==================== =======================
      Symbole             Signification             Symbole             Signification
    =================  ======================= ==================== =======================
      ``a and b``         a et b                   ``a>b``              a>b 
     ``a or b``           a ou b                       ``a>=b``            a≥b
     ``a==b``             égalité de a et b        ``a != bb``         a différent de b
     ``a<b``              a<b                       ``a<=b``            a≤b
    =================  ======================= ==================== =======================

Comme dans Matlab/Octave/Scilab, Numpy étend ``==``, ``<``, ``>``, ``<=``, ``>=`` et ``!=``
aux vecteurs/matrices (type ``array``) et ajoute le "ou exclusif" avec ``np.logical_xor(a,b)``. 
Les opérateurs de comparaisons ``==``, ``<``, ``>``, ``<=``, ``>=`` et ``!=`` existent aussi respectivement sous les noms  de
``equal(a,b)``, ``np.less(a,b)``, ``np.greater(a,b)``, ``np.less_equal(a,b)``, ``np.greater_equal(a,b)`` 
et ``np.not_equal(a,b)``. 

Ce qui est formidable avec Numpy (si si) : ``a`` et ``b`` étant deux ``array`` de
tailles identiques, la commande ``a<b`` renvoie un ``array`` de taille égale à celle de ``a`` dont
les éléments sont des booléens et correspondent au test ``<`` terme à terme. Voici
une illustration plus compréhensible

>>> b=6*np.eye(4,4)-np.diag(4+np.arange(3),1)+2*np.ones((4,4))
>>> a=np.arange(16).reshape(4,4)
>>> (a<b)
array([[ True, False, False, False],
       [False,  True, False, False],
       [False, False, False, False],
       [False, False, False, False]], dtype=bool)
>>> (a==b)
array([[False, False,  True, False],
       [False, False, False, False],
       [False, False, False, False],
       [False, False, False, False]], dtype=bool)

Si les deux arguments sont de tailles différentes, Python vous insulte. Mais si
``a`` ou ``b`` est de type ``float`` ou ``int``,  la commande ``a<2`` renvoie le tableau
de taille égale à celle de ``a`` correspondant à comparaison des éléments
de ``a`` avec 2.

>>> (a<2.)
array([[ True,  True, False, False],
       [False, False, False, False],
       [False, False, False, False],
       [False, False, False, False]], dtype=bool)
>>> (2<b)
array([[ True, False, False, False],
       [False,  True, False, False],
       [False, False,  True, False],
       [False, False, False,  True]], dtype=bool)

Encore plus formidable, si ``d`` est un tableau de booléens de taille
égale à celle de ``a``, ``a[d]`` extrait les éléments de ``a``
qui correspondent à la valeur ``True`` dans ``d``. La commande
``b[b<4]`` renvoie les éléments de ``b`` strictement inférieurs à 4.

>>> b=6*np.eye(4,4)-np.diag(4+np.arange(3),1)+2*np.ones((4,4))
>>> b
array([[ 8., -2.,  2.,  2.],
       [ 2.,  8., -3.,  2.],
       [ 2.,  2.,  8., -4.],
       [ 2.,  2.,  2.,  8.]])
>>> b[b<4]
array([-2.,  2.,  2.,  2., -3.,  2.,  2.,  2., -4.,  2.,  2.,  2.])

(à faire : ``np.all()`` ``np.any()``, `np.where``)

Autres commandes
-----------------

Numpy permet de "vectoriser", i.e. appliquer une fonction à un vecteur/matrice et
éviter les boucles. Comme
nous avons choisi d'utiliser Numpy à travers ``import numpy as np``, il faut
choisir les fonctions usuelles définies dans Numpy

>>> from math import cos
>>> a=np.arange(4, dtype=float)
>>> a
array([ 0.,  1.,  2.,  3.])
>>> cos(a)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: only length-1 arrays can be converted to Python scalars
>>> np.cos(a)
array([ 1.        ,  0.54030231, -0.41614684, -0.9899925 ])

Minimum/maximum global ou selon une direction : ``np.amin`` et ``np.amax``

>>> a = np.arange(4).reshape((2,2))
>>> a
array([[0, 1],
       [2, 3]])
>>> np.amin(a)           # élément minimum
0
>>> np.amin(a, axis=0)         # minimum selon ligne
array([0, 1])
>>> np.amin(a, axis=1)         # minimum selon colonne
array([0, 2])

Les commandes ``np.argmin`` et ``np.argmax`` retournent l'indice ou le tableau d'indices des éléments réalisant les extrema. Si plusieurs éléments réalisent l'extrema seul le premier indice est retourné.

>>> b=6*np.eye(4,4)-np.diag(4+np.arange(3),1)+2*np.ones((4,4))
>>> b
array([[ 8., -2.,  2.,  2.],
       [ 2.,  8., -3.,  2.],
       [ 2.,  2.,  8., -4.],
       [ 2.,  2.,  2.,  8.]])
>>> np.argmin(b)
11
>>> np.argmax(b)
0
>>> np.argmin(b,axis=1)
array([1, 2, 3, 0])
>>> np.argmin(b,axis=0)
array([1, 0, 1, 2])

Les sommes, produits de tous les éléments ou selon les lignes/colonnes se font
aisément avec ``np.sum`` et ``np.prod``. À vous de vérifier que les résultats
sont justes.

>>> b=6*np.eye(4,4)-np.diag(4+np.arange(3),1)+2*np.ones((4,4))
>>> np.sum(b)
41.0
>>> np.prod(b)
-50331648.0
>>> np.sum(b,axis=0)
array([ 14.,  10.,   9.,   8.])
>>> np.sum(b,axis=1)
array([ 10.,   9.,   8.,  14.])
>>> np.prod(b,axis=1)
array([ -64.,  -96., -128.,   64.])

Pour faire les moyennes de tous les éléments ou selon les lignes/colonnes, ``np.mean`` a le fonctionnement attendu.