Une très grande partie des exemples donnés dans ce tutoriel vient directement de https://cs231n.github.io/python-numpy-tutorial/ (sous licence MIT) :
The MIT License (MIT)
Copyright (c) 2015 Andrej Karpathy
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
import numpy as np
import matplotlib.pyplot as plt
Pour vous lancer dans le premier devoir, voici un code Python calculant l'interpolation poynomiale par la technique de la matrice de Vandermonde. Il ne s'agit donc pas de la solution du devoir :-), mais peut être une bonne source d'inspiration...
On observe que le graphe ci-dessous est exactement le même que celui obtenu à la toute fin de ce notebook avec polyval-polyfit (voir point 2).
Notez qu'on peut aussi remplir la matrice de Vandermonde en 1 seule ligne de code en utilisant les compréhensions de liste de Python (voir lignes en commentaire).
Dans la suite de ce notebook se trouve plusieurs exemples sur les fonctions utilisées dans ce petit script. Pour en citer quelques unes :
################################################################
def interpol(X, U):
"""
Calcule les coefficients a = [a_0, a_1, ..., a_n] du polynome d'interpolation
passant par les points (X_i, U_i).
"""
N = len(X)
V = np.zeros((N, N))
for i in range(N):
V[:, i] = X**i
# ou
# V = np.array([X**i for i in range(len(X))]).T
a = np.linalg.solve(V, U) # Résout le système linéaire.
return a
def eval_interpol(a, x):
"""
Evalue le polynome p(x) = a_0 + a_1 x + a_2 x^2 + ... + a_n x^n
"""
N = len(a)
p = np.zeros((N, len(x)))
for i in range(N):
p[i] = a[i]*x**i
return np.sum(p, axis=0)
# ou
# return np.sum([a[i]*x**i for i in range(a.size)], axis=0)
#################################################################
u = lambda x: np.cos(np.pi/2*x)
n = 5
X = np.linspace(-np.pi, np.pi, n)
U = u(X)
x = np.linspace(-np.pi, np.pi, 100*n)
uh = eval_interpol(interpol(X, U), x) # Calcule + évaluation du polynome d'interpolation
plt.figure()
plt.title(r"Interpolate $u(x) = \cos\left(\dfrac{\pi}{2} x\right)$")
plt.plot(x, u(x), label="original function")
plt.plot(x, uh, label="interpolation")
plt.plot(X, U, "or", label="interpolation points")
plt.grid()
plt.legend()
plt.show()
plt.close()
La suite de ce notebook donne des exemples sur les différents outils offerts par Numpy.
En cas d'incompréhension, manque d'explications, ... ne pas hésiter à se référer à la documentation officielle de numpy (et utiliser la barre de recherche): https://numpy.org/doc/stable/reference/generated/numpy.reshape.html#numpy-reshape
N'hésitez pas non plus à recopier les exemples et les modifier pour voir ce qu'il se passe, voir quand est-ce que le code ne marche plus, ...
a = np.array([1, 2, 3]) # Create a rank 1 array
print(type(a)) # Prints "<class 'numpy.ndarray'>"
print(a.shape) # Prints "(3,)". Corresponds to the dimension of the array as a tuple (m,n,...,p)
print(a[0], a[1], a[2]) # Prints "1 2 3"
<class 'numpy.ndarray'> (3,) 1 2 3
a[0] = 5 # Change an element of the array
print(f"{a = }\n") # Prints "[5, 2, 3]"
print("a =", a)
a = array([5, 2, 3]) a = [5 2 3]
b = np.array([[1,2,3],[4,5,6]]) # Create a rank 2 array
print(b.shape) # Prints "(2, 3)"
print(b[0, 0], b[0, 1], b[1, 0]) # Prints "1 2 4"
(2, 3) 1 2 4
a = np.zeros((2,2)) # Create an array of all zeros
print(f"{a = }\n") # Prints "[[ 0. 0.]
# [ 0. 0.]]"
b = np.ones((1,2)) # Create an array of all ones
print(f"{b = }\n") # Prints "[[ 1. 1.]]"
c = np.full((2,2), 7) # Create a constant array
print(f"{c = }\n") # Prints "[[ 7. 7.]
# [ 7. 7.]]"
d = np.eye(2) # Create a 2x2 identity matrix
print(f"{d = }\n") # Prints "[[ 1. 0.]
# [ 0. 1.]]"
e = np.random.random((2,2)) # Create an array filled with random values
print(f"{e = }\n") # Might print "[[ 0.91940167 0.08143941]
# [ 0.68744134 0.87236687]]"
a = array([[0., 0.], [0., 0.]]) b = array([[1., 1.]]) c = array([[7, 7], [7, 7]]) d = array([[1., 0.], [0., 1.]]) e = array([[0.53565381, 0.71381543], [0.61510839, 0.86926919]])
Numpy offre plusieurs façons d'indexer sur les arrays (bien plus puissant que ce qui est permis avec les listes de Python de base).
# Create the following rank 2 array with shape (3, 4)
# [[ 1 2 3 4]
# [ 5 6 7 8]
# [ 9 10 11 12]]
a = np.array([[1,2,3,4], [5,6,7,8], [9,10,11,12]])
print("a (original):\n\n", a)
print()
# Use slicing to pull out the subarray consisting of the first 2 rows
# and columns 1 and 2; b is the following array of shape (2, 2):
# [[2 3]
# [6 7]]
b = a[:2, 1:3]
# Remark:
# A slice of an array (b) is a view into the same data, so modifying it
# will modify the original array (a).
print(a[0, 1]) # Prints "2"
b[0, 0] = 77 # b[0, 0] is the same piece of data as a[0, 1]
print(a[0, 1]) # Prints "77"
print()
print("a (modified):\n\n", a)
print()
print("b:\n\n", b)
a (original): [[ 1 2 3 4] [ 5 6 7 8] [ 9 10 11 12]] 2 77 a (modified): [[ 1 77 3 4] [ 5 6 7 8] [ 9 10 11 12]] b: [[77 3] [ 6 7]]
# Create the following rank 2 array with shape (3, 4)
# [[ 1 2 3 4]
# [ 5 6 7 8]
# [ 9 10 11 12]]
a = np.array([[1,2,3,4], [5,6,7,8], [9,10,11,12]])
# Two ways of accessing the data in the middle row of the array.
# Mixing integer indexing with slices yields an array of lower rank,
# while using only slices yields an array of the same rank as the
# original array:
row_r1 = a[1, :] # Rank 1 view of the second row of a
row_r2 = a[1:2, :] # Rank 2 view of the second row of a
print(row_r1, row_r1.shape) # Prints "[5 6 7 8] (4,)"
print(row_r2, row_r2.shape) # Prints "[[5 6 7 8]] (1, 4)"
print()
# We can make the same distinction when accessing columns of an array:
col_r1 = a[:, 1]
col_r2 = a[:, 1:2]
print(col_r1, col_r1.shape) # Prints "[ 2 6 10] (3,)"
print(col_r2, col_r2.shape) # Prints "[[ 2]
# [ 6]
# [10]] (3, 1)"
[5 6 7 8] (4,) [[5 6 7 8]] (1, 4) [ 2 6 10] (3,) [[ 2] [ 6] [10]] (3, 1)
On se rappelle que la commande $\texttt{range}$ de Python suit la synthaxe suivante :
$\texttt{range(begin, end, step)}$
Il est aussi possible de définir le pas (step) lorsqu'on effectue un "slicing" sur un numpy array :
$\texttt{a[begin:end:step]}$
Ceci est illustré sur plusieurs exemples :
a = np.arange(1, 5.5, 0.5) # idem range but create a numpy array and allows float numbers for "begin", "end" and "step".
print(a)
print()
print(a[::2]) # start from 0 and finish in len(a)-1, but only select elements of a with even indexes (step=2).
print(a[1::3]) # start from 1 and finish in len(a)-1, with a step of 3.
print(a[:5:2]) # start from 0 and finish in 5, with step of 2
print(a[1:5:4]) # steart from 1 and finish in 5, with a step of 4
[1. 1.5 2. 2.5 3. 3.5 4. 4.5 5. ] [1. 2. 3. 4. 5.] [1.5 3. 4.5] [1. 2. 3.] [1.5]
Utiliser la synthaxe Python des listes pour indexer sur un numpy array fonctionne aussi... Mais attention, parfois cela ne produit pas du tout le même résultat !!!
a = np.array([[1,2,3,4], [5,6,7,8], [9,10,11,12]])
# List syntax :
print(a[:2][0]) # a[:2] gets the first two lists in a, that is [[1,2,3,4], [5,6,7,8]]
# a[:2][0] gets the first list in a[:2], that is [1,2,3,4]
# Numpy array syntax :
print(a[:2, 0]) # a[:2, :] gets the first two rows of a, that is the matrix [[1 2 3 4]
# [5 6 7 8]]
# a[:2, 0] gets the first column of the previous matrix, that is the vector [[1]
# [5]]
[1 2 3 4] [1 5]
Si $\texttt{a}$ est un numpy array de taille (n,m) et $\texttt{A,B}$ sont des numpy arrays d'entiers (ou des listes d'entiers) de taille $p$, alors
$\texttt{a[A,B]}$ est équivalent à $\texttt{[a[A[0],B[0]] a[A[1],B[1]] ... a[A[p-1],B[p-1]]]}$.
Exemple :
a = np.array([[1,2], [3, 4], [5, 6]])
# An example of integer array indexing.
# The returned array will have shape (3,) and
print(a[[0, 1, 2], [0, 1, 0]]) # Prints "[1 4 5]"
# The above example of integer array indexing is equivalent to this:
print(np.array([a[0, 0], a[1, 1], a[2, 0]])) # Prints "[1 4 5]"
# When using integer array indexing, you can reuse the same
# element from the source array:
print(a[[0, 0], [1, 1]]) # Prints "[2 2]"
# Equivalent to the previous integer array indexing example
print(np.array([a[0, 1], a[0, 1]])) # Prints "[2 2]"
[1 4 5] [1 4 5] [2 2] [2 2]
# Create a new array from which we will select elements
a = np.array([[1,2,3], [4,5,6], [7,8,9], [10, 11, 12]])
print(a) # prints "array([[ 1, 2, 3],
# [ 4, 5, 6],
# [ 7, 8, 9],
# [10, 11, 12]])"
print()
# Create an array of indices
b = np.array([0, 2, 0, 1])
# Select one element from each row of a using the indices in b
print(a[np.arange(4), b]) # Prints "[ 1 6 7 11]", that is : [a[0, 0] a[1, 2] a[2, 0] a[3, 1]]
print()
# Mutate one element from each row of a using the indices in b
a[np.arange(4), b] += 10
print(a) # prints "array([[11, 2, 3],
# [ 4, 5, 16],
# [17, 8, 9],
# [10, 21, 12]])
[[ 1 2 3] [ 4 5 6] [ 7 8 9] [10 11 12]] [ 1 6 7 11] [[11 2 3] [ 4 5 16] [17 8 9] [10 21 12]]
Numpy autorise aussi d'indexer sur un numpy array via un array de booléens de même taille. Cela permet de sélectionner les éléments d'un array satisfaisant une certaine condition, et ce en très peu de lignes de code :
a = np.array([[1,2], [3, 4], [5, 6]])
bool_idx = (a > 2) # Find the elements of a that are bigger than 2;
# this returns a numpy array of Booleans of the same
# shape as a, where each slot of bool_idx tells
# whether that element of a is > 2.
print(bool_idx) # Prints "[[False False]
# [ True True]
# [ True True]]"
print()
# We use boolean array indexing to construct a rank 1 array
# consisting of the elements of a corresponding to the True values
# of bool_idx
print(a[bool_idx]) # Prints "[3 4 5 6]"
print()
# We can do all of the above in a single concise statement:
print(a[a > 2]) # Prints "[3 4 5 6]"
[[False False] [ True True] [ True True]] [3 4 5 6] [3 4 5 6]
Pour des questions de performances ou de "robustesse", il est parfois préférable d'indiquer explicitement à Numpy quels types de données contiendra un array à son initialisation. Cela se fait via le paramètre $\texttt{dtype}$.
En particulier, si on initialise un array comme :
$\texttt{x = array([1, 2])}$,
numpy définira par défaut $\texttt{x}$ comme un array d'entiers. Cependant, si par la suite on souhaite que des éléments de $\texttt{x}$ deviennent des floats, cela peut poser quelques soucis (du type arrondis automatiques, ...)
x = np.array([1, 2]) # Let numpy choose the datatype
print(x, " ", x.dtype) # Prints "int64"
x = np.array([1.0, 2.0]) # Let numpy choose the datatype
print(x, " ", x.dtype) # Prints "float64"
x = np.array([1.0, 2.0], dtype=np.int64) # Force a particular datatype
print(x, " ", x.dtype) # Prints "int64"
x = np.array([1, 2], dtype=np.float64) # Force a particular datatype
print(x, " ", x.dtype) # Prints "float64"
[1 2] int32 [1. 2.] float64 [1 2] int64 [1. 2.] float64
Pour modifier la dimension d'un numpy array (i.e. la "shape"), on peut par exemple utiliser :
$\texttt{myarray.reshape((m, n))}$.
Cela fonctionne du moment que le nombre d'éléments de $\texttt{myarray}$ reste inchangé avant et après reshape. Dans l'exemple suivant, le array possède 6 éléments. Les seules reshape possibles sont donc : (1, 6), (6, 1), (2, 3) ou (3, 2).
a = np.arange(6).reshape((3, 2)) # array([[0, 1],
# [2, 3],
# [4, 5]])
print(a)
[[0 1] [2 3] [4 5]]
Dans le cas où l'on veut convertir une matrice en vecteur (comme fait classiquement en math), on peut utiliser $\texttt{flatten}$ ou $\texttt{ravel}$. Cela revient à faire $\texttt{myarray.reshape((m*n,))}$, si $\texttt{myarray}$ est de dimension (m,n) à l'origine.
La différence entre $\texttt{flatten}$ ou $\texttt{ravel}$ est assez subtile et est expliquée en détail à l'adresse suivante : https://www.geeksforgeeks.org/differences-flatten-ravel-numpy/.
print(a, " ", a.shape)
print()
b = a.flatten()
print(b, " ", b.shape)
c = a.ravel()
print(c, " ", c.shape)
d = a.reshape((np.prod(a.shape),)) # np.prod(x) = x[0]*x[1]*...*x[n] (voir plus bas pour plus d'exemples)
print(d, " ", d.shape)
[[0 1] [2 3] [4 5]] (3, 2) [0 1 2 3 4 5] (6,) [0 1 2 3 4 5] (6,) [0 1 2 3 4 5] (6,)
Les numpy arrays peuvent être vu comme des vecteurs et des matrices surlequels les opérations mathématiques suivantes sont autorisées :
Soit $A, B \in \mathbb{R}^{n\times m}$, deux matrices chacune représentée par un numpy array de taille (m,n) et $C \in \mathbb{R}^{m\times p}$. Soit $x \in \mathbb{R}^n$, un vecteur représenté par un numpy array de taille (n,).
1) Somme/Différence/Produit/Division composante par composante (uniquement sur des numpy arrays de même shape):
2) Somme/Différence/Produit/Division par un scalaire ($\alpha \in \mathbb{R}$) sur chaque composante :
3) Produit matriciel, produit scalaire et outer-product (uniquement sur des arrays compatibles en dimensions : cf. Algèbre)
Remarque : Lorsqu'on utilise $@$ sur des numpy arrays de shape (n,) (c'est-à-dire des vecteurs), numpy calcule toujours le produit scalaire classique, i.e. $x^T y$. Si jamais on souhaite calculer un "outer-product" : $x y^T$, il faut alors utiliser la fonction $\texttt{outer(x,y)}$ de numpy.
x = np.array([[1,2],[3,4]], dtype=np.float64)
y = np.array([[5,6],[7,8]], dtype=np.float64)
# Elementwise sum; both produce the array
# [[ 6.0 8.0]
# [10.0 12.0]]
print(x + y)
print(np.add(x, y))
print()
# Elementwise difference; both produce the array
# [[-4.0 -4.0]
# [-4.0 -4.0]]
print(x - y)
print(np.subtract(x, y))
print()
# Elementwise product; both produce the array
# [[ 5.0 12.0]
# [21.0 32.0]]
print(x * y)
print(np.multiply(x, y))
print()
# Elementwise division; both produce the array
# [[ 0.2 0.33333333]
# [ 0.42857143 0.5 ]]
print(x / y)
print(np.divide(x, y))
print()
# Matrix product; both produce the array
# [[19.0 22.0]
# [43.0 50.0]]
print(x @ y)
print(np.dot(x,y))
print()
# Elementwise square root; produces the array
# [[ 1. 1.41421356]
# [ 1.73205081 2. ]]
print(np.sqrt(x))
# All the usual math functions are implemented in numpy
print(np.cos(x))
print(np.exp(x))
[[ 6. 8.] [10. 12.]] [[ 6. 8.] [10. 12.]] [[-4. -4.] [-4. -4.]] [[-4. -4.] [-4. -4.]] [[ 5. 12.] [21. 32.]] [[ 5. 12.] [21. 32.]] [[0.2 0.33333333] [0.42857143 0.5 ]] [[0.2 0.33333333] [0.42857143 0.5 ]] [[19. 22.] [43. 50.]] [[19. 22.] [43. 50.]] [[1. 1.41421356] [1.73205081 2. ]] [[ 0.54030231 -0.41614684] [-0.9899925 -0.65364362]] [[ 2.71828183 7.3890561 ] [20.08553692 54.59815003]]
D'autres fonctions utiles permettent d'effectuer diverses opérations mathématiques sur les arrays :
# Sum of elements of an array :
x = np.array([[1,2],[3,4]])
print(np.sum(x)) # Compute sum of all elements; prints "10"
print(np.sum(x, axis=0)) # Compute sum of each column; prints "[4 6]"
print(np.sum(x, axis=1)) # Compute sum of each row; prints "[3 7]"
print()
# Product of elements of an array :
print(np.prod(x)) # Compute product of all elements; prints "24"
print(np.prod(x, axis=0)) # Compute product of each column; prints "[3 8]"
print(np.prod(x, axis=1)) # Compute product of each row; prints "[2 12]"
print()
# Transpose of a matrix :
x = np.array([[1,2], [3,4]])
print(x) # Prints "[[1 2]
# [3 4]]"
print(x.T) # Prints "[[1 3]
# [2 4]]"
print()
# Note that taking the transpose of a rank 1 array does nothing:
v = np.array([1,2,3])
print(v) # Prints "[1 2 3]"
print(v.T) # Prints "[1 2 3]"
10 [4 6] [3 7] 24 [3 8] [ 2 12] [[1 2] [3 4]] [[1 3] [2 4]] [1 2 3] [1 2 3]
a = np.array([[1, 2], [3, 4]])
b = np.array([[5, 6]])
conc1 = np.concatenate((a, b), axis=0) # array([[1, 2],
# [3, 4],
# [5, 6]])
print(conc1,"\n")
conc2 = np.concatenate((a, b.T), axis=1) # array([[1, 2, 5],
# [3, 4, 6]])
print(conc2,"\n")
conc3 = np.concatenate((a, b), axis=None) # array([1, 2, 3, 4, 5, 6])
print(conc3,"\n")
[[1 2] [3 4] [5 6]] [[1 2 5] [3 4 6]] [1 2 3 4 5 6]
a = np.linspace(0, 1, 5)
print(a)
[0. 0.25 0.5 0.75 1. ]
En particulier, si $n = \texttt{len(X)}-1$, les coefficients calculés correspondent à ceux du polynome d'interpolation passant par les points $\texttt{(X[i], U[i])}$ (i.e. $R = 0$).
Polyval et polyfit combinés ensemble permettent de calculer en 1 seule ligne de code le polynome d'interpolation passant par les points (X_i, U_i) :
u = lambda x: np.cos(np.pi/2*x)
n = 5
X = np.linspace(-np.pi, np.pi, n)
U = u(X)
x = np.linspace(-np.pi, np.pi, 100*n)
uh = np.polyval(np.polyfit(X, U, n-1), x) # Calcule + évaluation du polynome d'interpolation
import matplotlib.pyplot as plt
plt.figure()
plt.title(r"Interpolate $u(x) = \cos\left(\dfrac{\pi}{2} x\right)$")
plt.plot(x, u(x), label="original function")
plt.plot(x, uh, label="interpolation")
plt.plot(X, U, "or", label="interpolation points")
plt.grid()
plt.legend()
plt.show()
plt.close()