Exemples de langages interprétés : Python, Matlab, Scilab, Octave, R
Exemples de langages compilés : C, C++, Java
Vitesse d’exécution : interprété < compilé
Pourquoi, dans ce cas, considérer Python pour du calcul scientifique ?
Temps d’implémentation vs temps d’exécution : langage lisible et épuré ⇒ développement et maintenance rapides
Exécution en Python : rapide si les passages critiques sont exécutés avec un langage compilé : de nombreuses fonctions sont compilées et le code est interfaçable avec C/C++/FORTRAN
Tutoriel/aide dans la console : nom_fonction?
Voyons maintenant les principaux types et opérations sous Python.
Deux valeurs possibles : False, True
type(False)
x = (1<2)
x=3
type(x)
Opérateurs de comparaison : ==, !=, >, >=, <, <=
2 <=8 <15
Opérateurs logiques : not, or, and
(3 == 3) or (9 > 24)
(9 > 24) and (3 == 3)
type(2**100)
type(3.6)
type(3+2j)
2*3 # produit
#2**3 # puissance
#20/3 # division flottante
#20//3 # division entière
#20%3 # modulo
#(9+5j).real
#(9+5j).imag
#abs(3+4j)
type('abc')
c1 = 'L’eau vive'
c2 = ' est "froide" !'
c1+c2 # concatenation
#c1*2 #repetition
#c1[2]
#c1[-2]
#c1[2:5]
#len("abc") # longueur
#"abcde".split("c") # scinde
#"a−ha".replace('a','o') # 'o-ho'
#'-'.join(['ci', 'joint']) # 'ci-joint'
#'abracadabra'.count('bra') # 2
#'PETIT'.lower() # 'petit'
#'grand'.upper() # 'GRAND'
list : collection hétérogène, ordonnée et modifiable d’éléments séparés par des virgules, entourée de crochets.
my_list=[4,7,3.7,'E',5,7]
#type(my_list) # list
#len(my_list) # longueur
#my_list[0] # premier terme
#my_list[1:3] # [7, 3.7]
#[0,1] + [2,4] # [0,1,2,4] (concatenation)
#l = [2,4,5,9,1,6,4]
#k = [x for x in l if x<6] #(extract under condition)
#range(6) # attention différent en Python 2 (liste) et 3 (boucle)
#range(2, 9, 2)
#[x for x in range(3,6)]
nombres = [17, 38, 10, 25, 72]
nombres.sort()
#nombres.append(12)
#nombres.reverse()
#nombres.remove(38)
#print(nombres.index(17))
#nombres[0] = 11
#nombres[1:3] = [14, 17, 2]
#nombres.count(17) # 2
nombres
ATTENTION, par défaut, en Python, les listes ne sont pas copiées.
x = [4, 2, 10, 9, "toto"]
print(x)
y = x # y: seulement un alias de x, pas une copie
y[2] = "tintin" # change x(2) en "tintin"
print(x)
#x = [4, 2, 10, 9, "toto"]
#y = x[:] # On demande une copie
#y[2] = "tintin" # modifie y mais pas x
#print(x) # [4, 2, 10, 9, "toto"]
#print(y) # [4, 2, "tintin", 9, "toto"]
tuple : collection hétérogène, ordonnée, immuable.
t=(5,7)
type(t) # tuple
t[0] # 5
#t[0]=2 # error: item assignment for tuple
x = 11
if x < 0:
print("x est negatif")
elif x % 2:
print("x est positif et impair")
else:
print("x n'est pas negatif et est pair")
print("Eh oui !")
N = 0
x = 687687.567476
while (x > 0):
x//=2
N+=1
print("Approx. de log_2(x) : " + str(N-1))
for lettre in "ciao":
print(lettre)
for x in ["\n",2,'a', 3.14,"\n"]:
print(x)
Ecrire une fonction f qui prend $n$ en entrée et calcule $$\sum_{k=1}^{n-1} k^2$$
def f(n):
print(sum([k**2 for k in range(n)]))
f(3)
Ecrire une fonction f qui calcule le nième terme de la suite de Fibonacci initialisée par a et b, dont la relation de récurrence est $$ x_{n+2} = x_{n+1} + x_n.$$
def fib(n,a=0,b=1): #0,1: val. par defaut
'''n-th term of Fibonacci sequence starting at a,b.'''#tutoriel de fib
for i in range(n):
z = a+b
a=b
b=z
return a
fib(7,2,2)
Calcul numérique = manipulations de nombres décimaux $\neq$ Calcul symbolique = manipulation d’expressions symboliques
Exemple : racines de $x^2 − x − 1 = 0$
import ma_bibliotheque
ma_bibliotheque.la_fonction(...)
import ma_bibliotheque as bibli # raccourci
bibli.la_fonction(...)
Moins précis (car la bibliothèque d’origine des fonctions n’est pas précisée à leur appel) : from ma_bibliotheque import la_fonction
from ma_bibliotheque import la_fonction
la_fonction (...)
from ma_bibliotheque import ∗
la_fonction (...)
Attention si on importe plusieurs bibliothèques ayant les mêmes fonctions...
Numpy
¶import numpy as np # toujours commencer par cette commande quand on veut utiliser la bibliothèque numpy
b = np.array([[8,3,2,4],[5,1,6,0],[9,7,4,1]])
print(b)
#type(b) # numpy.ndarray
#b.dtype # datatype: int
b.shape # (3,4)
#c = np.array([[8,2],[5,6],[9,7]], dtype=complex)
#c.dtype # datatype: complex
#c[0,0] # 8+0j
#More than 2 dimensions:
#d = np.array([[[8,3],[1,2]],[[5,1],[4,5]],[[9,7],[4,5]]])
#d.shape # (3, 2, 2)
##Reshaping:
#x = d.reshape(4,3) # tableau de taille (4,3)
#d.reshape(12,1) # tableau de taille (12,1)
#d.reshape(12,) # tableau unidimensionel de taille 12
#np.insert(np.arange(4,9),3,17) # 4,5,6,17,7,8
Attention, il existe aussi une classe numpy.matrix mais il est recommandé d'utiliser numpy.array.
X = np.arange(start=5,step=3,stop=16) # 5, 8,11,14
#A = np.ones((2,3)) # matrix filled with ones
#B = X.reshape(2,2)
#C = np.zeros((3,2)) # matrix filled with zeros
#D = np.eye(2) # identity matrix
#np.diag([1,2]) # diagonal matrix
#E = C+np.ones(C.shape) # addition: same as C+1
#F = B*D # entry-wise multiplication
#J=np.dot(B,D) # linear algebra product
#G = F.T # transpose matrix
#H = np.exp(G) #as most functions, exp is entry-wise(else use np.vectorize(my_function))
#x = np.array([4, 2, 1, 5, 1, 10])
#y=np.logical_and(x>=3, x<= 9, x!=1) # [T,F,F,T,F,F]
#x[y]
#print(np.mean(np.random.randn(1000)>1.96))
A = np.array([[2, 1, 1], [4, 3, 0]])
B = np.array([[1, 2], [12, 0]])
C = np.array([[1, 2], [12, 0], [-1, 2]])
D = np.array([[1, 2, -4], [2, 0, 0], [1, 2, 3]])
E = np.concatenate((A,B), axis=1)
F = np.concatenate((C,D), axis=1)
G = np.concatenate((E,F),axis=0)
H = np.random.randn(5,5)
I = H*G # produit terme à terme
B5 = B**5 # puissance terme à terme
B5 =np.linalg.matrix_power(B, 5)
Bm1 = np.linalg.inv(B)
dB = np.linalg.det(B)
x = np.linalg.solve(B,[3,12]) #résout B*x=[[3],[12]]
print(E)
A = np.array([[1, 2], [12, 3]])
x = np.linalg.eigvals(A) # eigenvalues
# eigenvalues and eigenvectors:
valp, vectp = np.linalg.eig(A)
#Hermitian matrices methods:
S = np.array([[1, 2], [2, 3]])
y = np.linalg.eigvalsh(S) # eigenvalues
# eigenvalues and eigenvectors:
valp, vectp = np.linalg.eigh(S)
#Singular Value Decomposition
U,s,V=np.linalg.svd(A)
Ap = np.matrix(U)*np.diag(s)*V
print(A-Ap)
Attention encore une fois !!!
x = np.array([[8,3,2],[5,1,0],[9,7,1]])
y = x
x[0,0]+=1
x[0,0]-y[0,0] # 0
z=x.copy()
x[0,0]+=1
x[0,0]-z[0,0] # 1
import numpy.random as npr
my_sample = npr.ma_loi(paramètres, taille_du_tableau)
npr.randint(low=a,high=b,size=n) : v.a. unif. sur [ a, b[
npr.choice([a1,...,an],p=[p1,...,pn],size=n) : tirages indép. dans [a1,...,an] de loi [p1,...,pn]
npr.permutation(mon_urne) : permutation de mon_urne
npr.binomial(N,p,size=n)
npr.geometric(p,size=n)
npr.multinomial(n,tableau_des_probas,size=n)
npr.poisson(alpha,size=n)
Beaucoup d’autres exemples sur http://docs.scipy.org/doc/numpy/reference/routines.random.html
import numpy.random as npr
my_sample = npr.ma_loi(paramètres, taille_du_tableau)
npr.rand(d1,d2,...) : tableau d1 x d2 x... de v.a.i. unif. sur [0, 1]
npr.uniform(low=a,high=b,size=n) : v.a.i. unif. sur [a, b[
(size=n peut être remplacé par size=(d1,d2,...), comme partout dans ce qui suit)
Beaucoup d’autres exemples sur http://docs.scipy.org/doc/numpy/reference/routines.random.html
** -> Eviter si possible les boucles en Python !!!**
from time import time
n = int(1e7)
# Methode 1. Boucle for
t1 = time()
gamma1=sum([1./i for i in range(1,n+1)]) - np.log(n)
t2 = time()
temps1 = t2 - t1
# Methode 2. Numpy
t1 = time()
gamma2=np.sum(1. / np.arange(1,n+1)) - np.log(n)
t2 = time()
temps2 = t2 - t1
print("Facteur de gain: ", temps1/temps2)
from timeit import timeit
N = 100
setup = """""
import numpy as np
n = int(1e5)
"""""
code_boucle = """
np.sum([1. / i for i in range(1, n)]) - np.log(n)
"""
time_boucle=timeit(code_boucle,setup=setup,number=N)
code_numpy = """
np.sum(1. / np.arange(1, n)) - np.log(n)
"""
time_numpy=timeit(code_numpy,setup=setup,number=N)
print("Facteur : {}".format(time_boucle/time_numpy))
matplotlib
¶import matplotlib.pyplot as plt
plt.plot(x,y) affiche la courbe affine par morceaux reliant les points d’abscisses x et d’ordonnées y (nombreuses options) pour x, y vecteurs de même dimension,
plt.hist trace un histogramme. Deux options pour les colonnes : bins= nombre de colonnes ou bins= abscisses des séparations des colonnes
plt.bar trace un diagramme en bâtons
plt.scatter(x,y) affiche le nuage de points d’abs. x et d’ord. y
plt.stem(x,y) affiche des barres verticales d’abs. x et hauteur y
plt.axis([xmin,xmax,ymin,ymax]) définit les intervales couverts par la figure
plt.axis(’scaled’) impose que les échelles en x et en y soient les mêmes
plt.show() affiche les fenêtres créées dans le script
plt.figure() crée une nouvelle fenêtre graphique
plt.title("mon titre") donne un titre à une figure
plt.legend(loc=’best’) affiche la légende d’un graphique (en position optimale)
plt.subplot subdivise la fenêtre graphique de façon à y afficher plusieurs graphiques
import scipy.stats as sps
import matplotlib.pyplot as plt
%matplotlib inline
n, p, N = 20, 0.3, int(1e4)
B = np.random.binomial(n, p, N)
f = sps.binom.pmf(np.arange(n+1), n, p)
plt.hist(B,bins=n+1,density = 1,range=(-.5,n+.5),color="white",label="loi empirique")
plt.stem(np.arange(n+1),f,"r",label="loi theorique")
plt.legend()
E = np.random.randn(int(1e5))#echantillon
x = np.linspace(-4,4,1000)
f_x = sps.norm.pdf(x) #Densite gaussienne
plt.plot(x,f_x,"r",label="Theory")
#Affichage histo:
plt.hist(E,bins=50,density=1,label="Data")
plt.legend(loc='best')
Soient $(X_i)_{i\ge 1}$ copies indépendantes d'une même variable aléatoire $X$:\begin{eqnarray*} \mathbb{E}[|X|] <\infty&\implies& S_n:=\frac{X_1+\dots+X_n}{n}%\underset{n\to\infty}{\stackrel{p.s.}{\longrightarrow}} \underset{n\to\infty}{\longrightarrow} \mathbb{E}[X],\\ \mathbb{E}[|X|] =\infty&\implies& S_n:=\frac{X_1+\dots+X_n}{n}\;\textrm{ diverge}.\end{eqnarray*}
wn=10000
S=np.cumsum(np.random.rand(n))/np.arange(1,n+1)
plt.plot(range(1,n+1),S,'r',label="S_n")
plt.plot((1,n),(.5,.5),"b--",label="Esperance")
plt.ylabel('S_n')
plt.xlabel("n")
plt.legend(loc='best')
plt.title("LGN")
alpha, n = 0.9, 1000
U = np.random.rand(n)
X = (2*np.random.randint(0,2,n)-1)*U**(-1/alpha)
S=np.cumsum(X)/np.arange(1,n+1)
plt.title("LGN: convergence ou pas selon alpha")
plt.plot((0,n),(0,0),"b--")
plt.plot(S,"r",label="Moyenne empirique")
plt.xlabel("n")
plt.ylabel("S_n")
plt.legend(loc='best')
tel que ${n^\beta(S_n-\mu)\longrightarrow \ell\ne 0}$, i.e. $${S_n\approx \mu +\frac{\ell}{n^\beta}}.$$ Pour $\beta$ fixé, afficher sur un même graphique plusieurs courbes $(n, n^\beta(S_n-\mu))$ à l'aide d'une boucle for. Essayez par exemple $\beta =0.2$, $\beta = 0.5$, $\beta = 0.7$.
beta = 0.8
n=1000
for k in range(10):
S=np.cumsum(np.random.rand(n))/np.arange(1,n+1)
T = (S - 0.5)*np.arange(1,n+1)**beta
plt.plot(range(1,n+1),T)
plt.plot((1,n),(0,0),"b--")
plt.ylabel('n^beta(S_n-mu)')
plt.xlabel("n")
plt.legend(loc='best')
plt.title("n^beta (erreur dans la LGN)")
$(X_i)_{i\ge 1}$ v.a. indépendantes de même loi d'espérance $\mu$ et d'écart-type $\sigma$: $$T_n = \frac{n^{1/2}}{\sigma}\left(\frac{X_1+\dots+X_n}{n}-\mu\right)\underset{n\to\infty}{\stackrel{loi}{\longrightarrow}} \mathcal{N}(0,1). $$
En d'autres termes, la moyenne empirique des $X_i$ $\approx$ moyenne théorique, avec une erreur aléatoire gaussienne d'ordre ${1/\sqrt{n}}$.
n,m,sigma=1000,10000,3**(-1/2)
X=2*np.random.rand(m,n)-1
S=np.sum(X,axis=1)/(np.sqrt(n)*sigma)
M=max(np.abs(S))
x=np.linspace(-M,M,1000)
y=sps.norm.pdf(x)
plt.plot(x,y,'r',label="densite")
plt.hist(S,bins=int(round(m**(1./3)*M*.5)),density=1,histtype='step',label="Histogramme")
plt.legend(loc='best')
plt.title("TCL")
Généralement, une grande somme de petits aléas peu corrélés fluctue autour de sa moyenne selon une distribution gaussienne.
Exemples : anatomie (ex : taille des individus de sexe donné), QI, nombreuses mesures physiques, données économiques, etc...
Ecriture a la fin d'un fichier existant:
Lecture:
Voir aussi np.save, np.load, np.savetxt, csv.reader...
f=open("PopLynxRegionCanada_1821_1934.dat","r")
ytxt=f.readlines() # list de str
y=[int(row) for row in ytxt] # convertit str en int
plt.plot(range(1821,1935),y,"r")
plt.title("Population de lynx")
plt.tight_layout() # pratique pour l’export
plt.show()
import numpy as np
import scipy
import matplotlib.pyplot as plt
def f(x):
return np.exp(x)+x
#zeros of f, computed starting at -0.2:
a=scipy.optimize.fsolve(f,-.2)
print(a,f(a))
#integral of f from 0 to 1:
b=scipy.integrate.quad(f,0,1)
print(b)
def g(y,t):
return y
T=np.arange(start=0,stop=1,step=.001)
##solution, at T, of y’=g(y,t), y(T[0])=1:
y=scipy.integrate.odeint(g,1,T)
plt.plot(T,np.log(y),"r")
plt.show()