Marc BUFFAT

Marc BUFFAT

Professeur au département de Mécanique, Lyon 1 e-mail

Blog scientifique et pédagogique utilisant des notebooks IPython et Linux

cours en ligne INPROS: chapitre 6

Ipython notebook : cours INPROS LyonHPC

Auteur: Marc BUFFAT, Pr dpt de Mécanique, UCB Lyon 1

Contributeurs: Violaine Louvet, Michel Kern, Loic Gouarin, Laurence Viry </h5>

Licence Creative Commons
Mise à disposition selon les termes de la Licence Creative Commons
Attribution - Pas d’Utilisation Commerciale - Partage dans les Mêmes Conditions 2.0 France.
In [ ]:
%matplotlib inline
%autosave 300
import numpy as np
from IPython.display import HTML,display
css_file = 'style.css'
try:
    display(HTML(open(css_file, "r").read()))
    print "using ",css_file
except :
    print "using default css"

Les tableaux (vecteur, matrice)

LyonHPC LyonHPC

Contenu de la leçonLyonHPC

  1. Les tableaux 1D
    1. Tableaux numériques en Python
    2. Manipulation de tableaux 1D en Python
  2. Les tableaux 2D
    1. Matrices en Python
    2. Exemple: produit matrice-vecteur
  3. Boucle en Python sur les tableaux
  4. Notation matricielle à la Matlab

Tableau 1D

structure de données pour représenter un vecteur

tableau 1D : liste de n éléments du même type (entier ou réel) (n = dimension du tableau)

  • en mathématique: $\mathbf{X} \in \mathcal{R}^n $
  • en algorithmique: X tableau de n réels (contenant les composantes de X)

élément d’un tableau : X[i] $\equiv\mathbf{X}_{i+1}$

attention on compte à partir de 0, donc i varie de 0 à n-1

Exemple: point dans l’espace

le point P dans l’espace $\mathcal{R}^3$ de coordonnées cartesiennes: $$ \mathbf{OP}=\left[\begin{array}{c} 1.\ 3.\ -0.5 \end{array}\right] $$ est représenté par un tableau de réels de dimension 3, tel que:

P=[1.,3.,5.] (i.e. P[0]=1., P[1]=3.0, P[2]=5.0)

Tableau numérique sous Python

  • plusieurs structures de données possible (list, array, ..)
  • mais on utilisera essentiellement la structure array de la bibliothéque numpy
In [ ]:
import numpy as np
P = np.array([1., 3., 5.])
print P
print type(P)
print P.size
print P[1]
print type(P[1])

visualisation de l’execution

In [ ]:
HTML('<iframe width="800" height="350" frameborder="0" src="http://pythontutor.com/iframe-embed.html#code=P%3D%5B1.,3.,5.%5D%0Ax%3DP%5B2%5D%0Ay%3D(P%5B0%5D%2BP%5B1%5D%2BP%5B2%5D)/3%0Az%3DP%5B3%5D&cumulative=false&heapPrimitives=false&drawParentPointers=false&textReferences=false&showOnlyOutputs=false&py=2&curInstr=4&codeDivWidth=350&codeDivHeight=400"> </iframe>')

Indice d’un tableau

attention pour un tableau de dimension n : indice $0 \le i \lt n$

In [ ]:
print P[2]
print P[3]

Manipulation de tableaux 1D en Python

création d’un tableau

  • à partir d’une liste (entre [ ]):  array(liste)
  • dimension = nbre d’éléments de la liste
  • type = type des élèments de la liste
In [ ]:
Q = np.array([1., 2.])
print Q[0], type(Q[0])
I = np.array([1, 2])
print I[0], type(I[0])

tableaux particuliers

  • tableau de 0.0 ou de 1.0 : zeros(n) ou ones(n) ( n dimension)
  • tableau de n valeurs non initialisées (attention !!!): empty(n)
In [ ]:
X = np.zeros(5)
print X
Y = np.ones(3)
print Y
Z = np.empty(3)
print Z
In [ ]:
print np.arange(0, 10, 1)
print np.array(range(10))
print np.arange(0., 10., 1.5)
  • intervalle: linspace(debut, fin, nbre d’éléments)
In [ ]:
print np.linspace(0, 10, 11)
print np.linspace(0, 10, 3)
  • valeurs aléatoires: random.rand(nbre)
In [ ]:
print np.random.rand(5)
In [ ]:
X = np.array([i*i for i in range(5)])
print X
  • spécification du type des éléments du tableau: argument dtype (int,float,..)
In [ ]:
I = np.zeros(3, dtype=int)
print I
print type(I), type(I[0])

taille et élément d’un tableau

  • nbre d’éléments du tableau X: X.size ou len(X)
  • taille du tableau X: X.shape
In [ ]:
n = 5
X = np.zeros(n)
print X.size, len(X), X.shape
  • accès aux éléments: premier X[0] et dernier X[-1]
In [ ]:
X = np.arange(1., 2., 0.2)
print X
print X[0], X[-1]
In [ ]:
X = np.arange(0, 10, 1)
print X[:]
print X[1:3]
print X[::2]
print X[1:-1]

aliasing et copie de tableau

  • aliasing par défaut: Y=X ne crée pas un nouveau tableau mais un alias
    • de même Y=X[:] (contrairement aux listes)
  • copie explicite: Y[:]=X[:]
    • ou Y = X.copy()
In [ ]:
X=np.zeros(2)
Y=X
Z=np.zeros(2)
Z[:]=X[:]
Y[0]=1
Z[0]=2
print "X=",X,"Y=",Y,"Z=",Z

boucle sur les éléments d’un tableau

boucle for : for i in range(n):

In [ ]:
n = 5
X = np.zeros(n)
for i in range(n):
    X[i] = i*i
print X

Exemple: produit scalaire

Calcul du produit scalaire scal de 2 vecteurs X,Y: $$ \mathbf{X}\in\mathcal{R}^n , \mathbf{Y}\in\mathcal{R}^n \ \ scal =\sum_{i=1}^{n} X_i Y_i$$ Remarque si X=Y $\rightarrow$ norme $$ ||\mathbf{X}|| = \sqrt{\sum_{i=1}^{n} X_i^2}$$

Algorithme produit scalaire

Algorithme ProdScal(X,Y)
    n = dimension(X) (et Y)
    scal=0
    pour i de 0 a n-1
        scal = scal + X[i]*Y[i]
    fin pour
    retour scal

Programme Python produit scalaire

In [ ]:
def ProdScal(X,Y):
    """ calcul le produit scalaire de 2 vecteurs """
    n=X.size
    scal = 0
    for i in range(n):
        scal = scal + X[i]*Y[i]
    return scal
#
X = np.random.rand(2)
Y = np.random.rand(2)
print "X =", X
print "Y =", Y
print "X.Y =", ProdScal(X,Y), ' ou', np.dot(X, Y)

Exercices

  1. Test de colinéarité de 2 vecteurs

  2. Normalisation d’un vecteur

Test de colinéarité de 2 vecteurs

$$\mathbf{U} \in \mathcal{R}^n , \mathbf{V} \in \mathcal{R}^n \mbox{ si } \left|\mathbf{U}.\mathbf{V}\right| = \left\Vert\mathbf{U}\right\Vert \left\Vert\mathbf{V}\right\Vert \mbox{ alors } \mathbf{U} \mbox{ et } \mathbf{V} \mbox{ sont colinéaires} $$

Ecrire une fonction qui teste la colinéarité

    colineaire(U,V)

Attention à la précision

Normalisation d’un vecteur X

$$\mathbf{U} \in \mathcal{R}^n \mbox{ calcul } \mathbf{U} \leftarrow \frac{\mathbf{U}}{\left\Vert\mathbf{U} \right\Vert}$$

Ecrire une fonction python qui normalise un vecteur X

Xn = normalisation(X)

attention au cas X=0

Tableau 2D

structure de données pour représenter une matrice

tableau 2D : liste de $n*p$ éléments du même type (entier ou réel) (n,p = dimensions de la matrice)

  • en mathématique: $\mathbf{A} \in \mathcal{R}^{n,p} $
  • en algorithmique: A tableau de $n*p$ réels (contenant les composantes de A)

élément i,i du tableau : A[i,j] $\equiv\mathbf{A}_{i+1,j+1}$

Attention on compte à partir de 0, donc i varie de 0 à n-1 et j de 0 à p-1

Exemple: matrice de rotation dans $\mathcal{R}^2$

Rotation $\theta$ d’un repère : $$ \mathbf{R}=\left[\begin{array}{cc} \cos(\theta)& -\sin(\theta)\ \sin(\theta)& \cos(\theta) \end{array}\right] $$ représenté par un tableau de réels de dimension $2*2$, tel que:

R=[[cos(theta),-sin(theta)],[sin(theta),cos(theta)]]

(i.e. R[0,0]=cos(theta), R[0,1]=-sin(theta), R[1,0]=sin(theta), R[1,1]=cos(theta))

remarque theta est une variable contenant la valeur de l’angle

Matrices en Python

bibliothéque numpy: array, dimension couple (tuple) (n,p)

In [ ]:
import numpy as np
theta = np.pi / 3
R = np.array([[np.cos(theta), -np.sin(theta)],
              [np.sin(theta), np.cos(theta)]])
print R
print type(R)
print R.size, R.shape

création de tableaux 2D

  • explicite A = array([[1.,2.],[3.,4.]])
  • fonctions: zeros, ones,.. : identiques en 1D en remplaçant n par (n,p) (voir doc numpy)
In [ ]:
A = np.array([[1.,2.],[3.,4.]])
Z = np.zeros((3, 2))
X = np.reshape(np.arange(6),(2,3))
print A.shape,Z.shape,X.shape

taille et éléments d’un tableau 2D

  • nbre total d’éléments du tableau A: A.size ou len(A)
  • par dimension  A.shape
  • élément: A[i, j] (à partir de 0)
  • sous-matrices (convention identique en 1D)
    • ligne i : A[i, :]
    • colonne j : A[:, j]
In [ ]:
A = np.array([[1, 2, 3], [4, 5, 6]])
print "A=\n",A, "\nA11 =",A[1,1] ,"dim A=", A.shape
B = A[:, 0:2]
print "B=\n",B, "\nB11 =",B[1,1] ,"dim B=", B.shape
C = A[1, :]
print "C=\n",C, "\nC1 =",C[1], "dim C=",C.shape
A[1,1]=0
print "B=\n",B
print "C=\n",C

équivalence avec les tableaux 1D

$A$ tableau 2D (array 2D) $n*p$ réels $\equiv$ $B$ tableau 1D de même dimension globale $np$

les données sont identiques, seul l’accès est différent:

reshape(tableau,dimension) (attention aliasing !!)

  • $B = A$
  • $A_{i,j}$ pour $i=0,n-1$ et $j=0,p-1$
  • $B_k$ pour $k=0..np-1$
  • $A_{i,j}\equiv B_k$ avec $k=i*p+j$ (stocage par ligne C/C++)
In [ ]:
A = np.array([[1., 2., 3.], [4., 5., 6.]])
B = np.reshape(A ,6)
C = np.reshape(A, (3, 2))
print "dim A =", A.shape," A=\n",A
print "dim B =", B.shape," B=\n",B
B[2] = 0
print "B=\n",B
print "A=\n",A
print "C=\n",C

Copie et aliasing

aliasing = noms de variable différents, mais données partagées

  • aliasing par défaut sous python:  B=A
  • sinon copie explicite: C[:,:]=A[:,:]
In [ ]:
A=np.zeros((2,2))
B=A
B[0,0]=2
C=A.copy()
C[0,:]=[3,4]
print "A=",A,"\nB=",B,"\nC=",C

Exemple: produit matrice-vecteur

Calcul du produit matriciel Y=A X: $$ \mathbf{A}\in\mathcal{R}^n\times\mathcal{R}^p , \mathbf{X}\in\mathcal{R}^p , \mathbf{Y}\in\mathcal{R}^n \ \ Y_i =\sum_{j=1}^{p} A_{i,j} X_j$$

Algorithme produit matrice-vecteur Y = A X

Algorithme MatVec(A,X)
    n,p = dimension de A
    Y = tableau(n)
    pour i de 0 a n-1
        Y[i]=0.
        pour j de 0 a p-1
            Y[i]=Y[i]+A[i,j]*X[j]
        fin pour j
    fin pour i
    retour Y

Programme Python produit scalaire

In [ ]:
def MatVec(A,X):
    """ calcul le produit Y=A.X de A par X """ 
    n,p = A.shape
    Y = np.zeros(n)
    for i in range(n):
        S = 0.
        for j in range(p):
            S +=  A[i, j] * X[j]
        Y[i] = S
    return Y
#
n, p = 3, 2
A = np.reshape(np.arange(n*p), (n, p))
X = np.array(np.arange(p))
print "A =\n",A
print "X =",X
print "Y =",MatVec(A,X)

Exercices

  1. Calcul du produit de 2 matrices
  2. Calcul de la matrice de Vandermonde d’ordre n

Produit matrice-matrice

$\mathbf{C} = \mathbf{A}\ \mathbf{B} $ pour $\mathbf{A} \in \mathcal{R}^n\times\mathcal{R}^p$ , $\mathbf{B} \in \mathcal{R}^p\times\mathcal{R}^q$ et $\mathbf{C} \in \mathcal{R}^n\times\mathcal{R}^q$ $$ C_{i,j} = \sum_{k=1}^{p} A_{i,k} B_{k,j}$$

Matrice de Vandermonde

$$ V = \left[ \begin{array}{ccccc} 1 & x_0 & x_0^2 & \ldots & x_0^{n} \ 1 & x_1 & x_1^2 & \ldots & x_1^{n}\ \ldots & \ldots & \ldots & \ldots & \ldots \ 1 & x_n & x_n^2 & \ldots & x_n^{n}\ \end{array} \right] $$

Intervient dans le problème de la détermination du polynôme d’interpolation de degré $n$ passant par $n+1$ points $\{x_i,y_i\}_{i=0,n}$

Boucle en Python sur les tableaux

itérations sur les éléments val d’un tableau Tab (ou d’une liste)

for val in Tab:
In [ ]:
# calcul du carré de la norme d'un vecteur
X = np.random.rand(100)
norm2 = 0.0
for x in X:
    norm2 = norm2 + x * x
print norm2

itérations avec indice et valeur (enumerate)

for i,val in enumerate(Tab):
In [ ]:
# determination du max d'un tableau
X = np.random.rand(10)
print X
xmax = X[0]
imax = 0.
for i, x in enumerate(X):
    if x > xmax:
        xmax = x
        imax = i
print "xmax =", xmax, " pour i =", imax

boucle avec un indice décroissant (reversed)

for i in reversed(range(n)):

exemple : permutation circulaire vers la droite

In [ ]:
# permutation des valeurs d'un tableau vers la droite: X[i+1]=X[i]
X = np.random.rand(4)
print X
last = X[-1]
for i in reversed(range(1, X.size)):
    X[i] = X[i-1]
X[0] = last
print X

boucles de très grande taille: xrange

pour des boucles de très grandes tailles (n > 100 000) on peut utiliser

for i in xrange(n)

au lieu de

for i in range(n)

Exercice: modifier l’algorithme pour effectuer une permutation circulaire vers la gauche

Notation matricielle à la Matlab

Efficacité -> boucles dans la bibliothéque numpy

  • opérations (terme à terme) sur les vecteurs +,-,*,/
In [ ]:
X = np.random.rand(3)
Y = np.random.rand(3)
print "X=",X, "\nY=",Y, "\nXY=",X*Y
  • fonctions mathématiques (terme à terme)
In [ ]:
print np.cos(X)

Utilisez les fonctions de la bibliothéque numpy

  • opérations matricielle: fonction dot

    produit scalaire, produit matrice-vecteur, produit matriciel $$ C= dot(A,B)\ \rightarrow\ C_{ij} = \sum_{k} A_{ik} B_{kj}$$

In [ ]:
print "X.Y=",np.dot(X,Y)
print "|X|=",np.sqrt(np.dot(X,X))
A = np.random.rand(3,3)
print "A*X=",np.dot(A,X)
print "A*A=\n",np.dot(A,A)
  • produit tensoriel interne (somme sur la dernière dimension): fonction inner $$ C= inner(A,B)\ \rightarrow\ C_{ij} = \sum_{k} A_{ik} B_{jk}$$
In [ ]:
print "inner(X,Y)=",np.inner(X,Y)
print "inner(A,X)=",np.inner(A,X)
print "inner(A,A)=\n",np.inner(A,A)
  • produit tensoriel externe: fonction outer
$$Z = outer(X,Y)\ \rightarrow\ Z_{ij} = X_i Y_j $$
In [ ]:
print "outer(X,X)=\n",np.outer(X,X)
  • somme et produit : fonction sum et prod
  • ainsi que de nombreuses autres fonctions (inverse, déterminant, ..)

voir la documentation numpy: http://docs.scipy.org/doc/numpy/reference

Aide sur les fonctions : help() ou TAB

help(fonction) ou fonction?
In [ ]:
help(np.sum)

Bibliographie

transciption du célébre article de Edsger W. Dijkstra.

Fin de la leçon

LyonHPC LyonHPC