Corrigé TD 7 à 10 ¶
Les entiers et les réels¶
eps=1.0
while not (1.0+0.5*eps) == 1.0:
eps *= 0.5
print (eps, 2**(-52))
Rappelons que la partie fractionnaire d'un réel est codée en base 2 selon $\displaystyle \sum_{i=1}^{M} a_i 2^{-i}$ avec $a_i=0$ ou 1. Dans les lignes ci-dessus, eps commence à 1 et est divisé par deux à chaque itération de la boucle. Il va donc parcourir les fractions $2^{-i}$. Il est important dans le test de la boucle while de comparer 1.0+0.5eps à 1.0 et non 0.5eps à 0 pour que l'exposant ne joue pas (voir ci-après). La valeur de eps lorsque le test s'arrête est telle que $1.0+0.5*eps \equiv\underbrace{1.0000000000000_2}_{M bits} = 1.0$. Pour la valeur précédente, il correspond à une représentation pour laquelle on a en gros le dernier bit de la mantisse qui est non nul, soit environ $\underbrace{1.0000000000001_2}_{M bits}$ mais il peut être plus grand. D'après le site baseconvert, il semble que ce soit plutôt $\underbrace{1.00000000000011_2}_{M bits}$. On obtient donc une bonne estimation de la précision machine, qui de l'ordre de $2^{-M}$ avec ici $M=52$ sur une machine 64bits qui correspond à l'erreur d'arrondi entrainé par la troncation de la partie fractionnaire sur un nombre de bits fini. Sur le site baseconvert.com, on peut voir que 1.000000000000001 est le dernier décimal dont la représentation binaire sur 64bits est différente de 1.
eps=1.0
while not (0.0+0.5*eps) == 0.0:
eps *= 0.5
print (eps, 2**(-1022), 2**(-1022-52))
Dans ce cas, l'exposant n'est plus 0 dans le test et eps peut descendre à des valeurs beaucoup plus petites au fur et à mesure que l'on divise par deux. L'exposant le plus petit que l'on peut faire est -1022. En effet, l'exposant est codé sur 11bits avec un décalage de sorte que $E=0\equiv -1022$, $E=1022\equiv 0$ et $E=2047=2^{11}-1\equiv +1023$. L'exposant le plus petit possible est donc -1022 et le plus grand 1023. Cela correspond à l'ordre de grandeur des nombres les plus grands dont la mantisse est de l'ordre de 1 (c'est l'ordre de grandeur donné dans le cours $2^{-1022}\simeq 2.22507385851\times 10^{-308}$). Maintenant, la mantisse peut représenter un nombre très petit, de l'ordre de $2^{-52}$ si tout les bits sont mis à 0 sauf le dernier. Au total, le nombre le plus petit accessible est donc de l'ordre de $2^{-1022-52} \simeq 4.94065645841\times 10^{-324}$ qui est bien celui affiché par le code.
0.125+0.25
0.1+0.2
Cela est surprenant, on s'attendrait à ce que l'ordinateur pour calculer exactement, sans erreur d'arrondi cette addition toute simple.
from decimal import Decimal
print( Decimal("0.1") ) # affiche la représentation décimale à partir d'une string
print( Decimal(0.1) )# affiche la représentation décimale du réel 0.1 codé sur un mot mémoire
print( Decimal(0.1+0.2), Decimal(0.3) )
print( Decimal("0.1")+Decimal("0.2") )
print( Decimal(0.25), Decimal(0.125) )
print( Decimal(0.35-0.1) )
On voit que 0.1 comme float représenté sur la machine ne correspond en fait pas à la valeur attendue 0.1 en base décimale mais l'approximation 0.1000000000000000055511151231257827021181583404541015625... Cela est dû au fait que 0.1 a une représentation infinie en base 2 de la même manière que 1/3 a une représentation infinie en base 10 (base décimale) que l'on écrit souvent 0.33333333333.... Le site baseconvert nous permet de faire afficher le début tronqué de la représentation en binaire de 0.1: il s'agit de: 0.0001 1001 1001 1001 1001 1001 1001 1001 1001 1001 1001 1001 101. Du coup, dans les opérations suivantes, cette inexactitude de représentation due à la troncation se propage (un peu heureusement mais pas trop), et le résultat 0.1+0.2 n'est pas représenté par la représentation correspondant à 0.3 d'où l'erreur d'arrondi apparaissant dans l'affichage. En revanche, pour $0.25=2^{-2}$ et $0.125=2^{-3}$ qui sont représentable exactement, le résultat de leur addition est aussi exact. Par contre, pour $0.35-0.1$, bien que le résultat soit représentable exactement, les intermédiaires de calcul ne le sont pas et le résultat ne tombe pas sur celui attendu.
import numpy as np
print( np.uint8(1000), 1000%256 )
print( np.uint8(230)+np.uint8(120) )
print( np.uint8(-20), np.uint8(3.4*41) )
On observe qu'il faudra se méfier avec les conversions modulo 256 des valeurs intermédiaires ainsi que celles qui ont des résultats intermédiaires négatifs. La valeur maximale 256 est assez vite atteinte et peut conduire à un overflow.
Algèbre linéaire avec Numpy¶
import numpy as np
np.set_printoptions(precision = 4, suppress=True)
L = 10
x = np.arange(1,L+1)
V = np.sqrt(2/(L+1))*np.sin(np.pi*x/(L+1))
print( V, np.linalg.norm(V))
W = np.diag([-2.]*L,k=0) + np.diag([1.]*(L-1),k=1) + np.diag([1.]*(L-1),k=-1)
print (W)
A = np.dot(W,V)
print (A / V)
(val,U) = np.linalg.eigh(W)
print( "valeurs propres:", val)
print( "vecteurs propres:")
for i in range(L): print( i,":",U[:,i])
Um1 = np.linalg.inv(U)
D = np.diag(val)
print( np.dot(U, np.dot(D,Um1)))
print( np.dot(Um1, np.dot(W,U)))
Application à la dynamique d'une chaîne d'oscillateurs couplés¶
Omega = np.sqrt(-val); print( Omega)
Conditions initiales $Y^0 = U^{-1}X^0$
X0 = np.zeros( (L,) ); X0[0] = 0.2; X0[1] = 0.1; print( X0)
# en partant d'un état propre
# X0 = U[:,1]
Y0 = np.dot(Um1,X0); print( Y0)
Solution des équations: $X(t) = UY(t)$ avec $Y_j(t) = Y^0_j\cos(\omega_j t)$
t = 10
X = np.dot(U,Y0*np.cos(Omega*t))
print (X)
t = np.linspace(0.0,10.0,101)
allX = np.zeros( shape=(L,len(t)) )
for i in range(len(t)):
allX[:,i] = np.dot(U,Y0*np.cos(Omega*t[i]))
from pylab import *
%matplotlib inline
for x in range(L):
plot(t, allX[x,:] )
xlabel('temps')
ylabel('Position')
title("Dynamique d'oscillateurs couples")
grid(True)
show()
import matplotlib.pyplot as pl
import matplotlib.image as mpim
%matplotlib inline
img = mpim.imread("grenouille.jpg")
type(img)
print( img.shape )
L'image jpeg est importée sous la forme d'un tableau type ndarray de numpy qui a 3 dimensions.
pl.imshow(img)
On voir que la première dimension est la dimension verticale (lignes) et la seconde la dimension horizontale (colonnes), comme une matrice $N\times M$. Selon la numérotation des axes qui correspondent aux indices, les pixels sont numéroter de haut en bas et de gauche à droite avec donc le premier (indice 0,0) en haut à gauche et le dernier (indice $N-1,M-1$ en bas à droite.
pl.axis('off')
pl.imshow(img)
print( img[100,100,:])
Les deux premières dimensions correspondent à la position du pixel et la troisième au codage RVB de la couleur.
imgcopy = np.copy(img)
mpim.imsave("grenouille_copy.jpg",imgcopy)
On a $N\times M$ pixels contenant chacun 3 octets pour les couleurs RVB, on s'attend donc à une taille en Ko de
print( 250*400*3/1024,"Ko")
La taille est plus petite car le format jpeg est compressé.
Agir sur la structure spatiale¶
pl.imshow(img[::,::-1]) # réflection par rapport à la verticale
pl.imshow(img[::-1,::]) # réflection par rapport à l'horizontale
pl.imshow(img[::-1,::-1]) # rotation de 180°
pl.imshow(img[::2,::2]) # réduction de la taille par deux en sélectionnant les indices pairs
pl.imshow(img[:img.shape[0]//2:,:img.shape[1]//2:]) # extraction du quadrant supérieur gauche
nl, nc, s = img.shape # on récupère les dimensions dans nl et nc, s=3
rayon = 120 # rayon voulu en pixel
Y, X = np.ogrid[:nl,:nc] # crée le tableau 2D des indices de la photo
masque = (X-nc/2)**2 + (Y-nl/2)**2 > rayon**2 # définition de l'extérieur du cercle
# centré sur le milieu de l'image
print( masque[0,0] ) # le masque est une matrice de booléen
pl.imshow(masque, cmap=pl.cm.gray) # on affiche le masque, cmap est la palette de couleur
img_coupee = np.copy(img)
pixel_noir = np.array([0,0,0],dtype=np.uint8) # définition du pixel noir
img_coupee[masque] = pixel_noir # les indices du masque renvoyant True sont mis en noir
pl.imshow(img_coupee)
a,b = 150.,100. # demi grand axe et demi petit axe
masque = ((X-nc/2)/a)**2 + ((Y-nl/2)/b)**2 > 1 # définition de l'extérieur de l'ellipse
# centrée sur le milieu de l'image
pl.imshow(masque, cmap=pl.cm.gray) # on affiche le masque, cmap est la palette de couleur
img_coupee = np.copy(img)
pixel_noir = np.array([255,255,255],dtype=np.uint8) # définition du pixel blanc
img_coupee[masque] = pixel_noir # les indices du masque renvoyant True sont mis en blanc
pl.imshow(img_coupee)
Agir sur les couleurs¶
neg = 255*np.ones(img.shape,dtype=np.uint8)
neg -= img
pl.imshow(neg)
La première ligne crée une image blanche de même taille que la photo. La seconde soustrait les couleurs de l'image au blanc et renvoie donc la couleur complémentaire.
def linear(source,v):
res = np.zeros(source.shape,dtype=np.uint8)
ny, nx, s = source.shape
for i in range(ny):
for j in range(nx):
res[i,j] = np.array(source[i,j]*v,dtype=np.uint8)
return res
rouge = linear(img,np.array([1,0,0]))
vert = linear(img,np.array([0,1,0]))
bleu = linear(img,np.array([0,0,1]))
pl.subplot(1, 3, 1); pl.axis('off'); pl.imshow(rouge)
pl.subplot(1, 3, 2); pl.axis('off'); pl.imshow(vert)
pl.subplot(1, 3, 3); pl.axis('off'); pl.imshow(bleu)
pl.show()
pl.axis('off'); pl.imshow(rouge+vert+bleu)
blanc = 255*np.ones(img.shape,dtype=np.uint8)
magenta = blanc - vert
pl.axis('off'); pl.imshow(magenta)
# réécriture en une ligne de la fonction linear
def linear(source,v):
return np.array(source[:,:]*v,dtype=np.uint8)
def contraste(source,f):
return np.array(f(source[:,:]),dtype=np.uint8)
s = 5.0
a, b = np.tanh(s*(-127.5)/255.), np.tanh(s*(255.-127.5)/255.)
f = lambda x: 255*(np.tanh(s*(x-127.5)/255.)-a)/(b-a)
x = np.arange(256)
print( np.max(f(x)), np.min(f(x)), f(127.5))
pl.plot(x,f(x))
pl.imshow(contraste(img,f))
Compression par décomposition en valeurs singulières¶
construction de la matrice en niveau de gris.
def rgb2gray(rgb):
return np.dot(rgb, [0.299, 0.587, 0.144])
A = rgb2gray(img) # matrice en niveau de gris
pl.imshow(A, cmap=pl.cm.gray)
pl.clim(0,255)
U, s, V = np.linalg.svd(A)
pl.semilogy()
pl.plot(s/np.max(s))
kept = 25
S = np.diag(s[:kept])
R = np.dot(U[:,:kept],np.dot(S,V[:kept,:]))
pl.subplot(1, 2, 1)
pl.imshow(R,cmap=pl.cm.gray)
pl.clim(0,255); pl.axis('off')
pl.subplot(1, 2, 2)
pl.imshow(A,cmap=pl.cm.gray)
pl.clim(0,255); pl.axis('off')
taille_orig = A.size
taille_comp = U[:,:kept].size+S.size+V[:kept,:].size
print( "Taux de compression:",taille_orig/float(taille_comp))
mpim.imsave("PetiteGrenouille.jpg",R,cmap=pl.cm.gray)
n, bins, patches = pl.hist(A.flatten(), 256, range=(0.0,255), fc='k', ec='k')
