for i in range(200*200 // 4): x = random.randint(0, 199) y = random.randint(0, 199) image[x, y] = 1Puis nous l'affichons avec matplotlib :
import matplotlib.pyplot as plt plt.imshow(image, cmap=plt.get_cmap('gray'), interpolation='nearest', vmin=0.0, vmax=1.0, aspect='equal') plt.show()La fonction
imshow
permet :- de choisir une carte de couleurs (la liste est ici : Matplotlib Colormaps)
- de ne pas interpoler entre les pixels (avec
nearest
, sans cette option, il y aura de l'interpolation) - de préciser les valeurs min et max de l'image et leur correspondance dans la colormap. Si on ne précise pas,
imshow
fait correspondre la valeur max de l'image à la valeur max de la colormap (idem pour la valeur min).
conversion
qui convertit des coordonnées de pixels entières (coordonnées dans la matrice) en nombres complexes en utilisant le cadrage $-1.25\le a\le 1.25$ et $-1.25\le b\le 1.25$. Votre fonction prendra donc deux entiers en paramètres et renverra un nombre complexe.conversion
est correcte (pour des entrées dont vous pouvez facilement calculer la sortie…). Donnez un copié/collé des tests (que vous choisirez bien) que vous avez réalisés.- associer un nombre complexe $u_0$ à des coordonnées (i,j) de pixels sur l'écran : fonction
conversion
, prenant en paramètres deux entiers et renvoyant un nombre complexe - décider si un nombre complexe $u_0$ appartient à $J©$ ou pas : fonction
julia
prenant en paramètres deux complexes $u_0$ et $c$ et décidant si $u_0$ appartient à $J©$ ou non. Cette fonction renvoie un 0 ou 1.
def calcul_img(c) | img <- création de la matrice à la bonne taille | Pour chaque coordonnée (i,j) de la matrice : | | u <- conversion(i,j) | | img[i,j] <- julia(u,c) | retourner imgAttention à la boucle de calcul qui en cache en fait 2 : une pour les abscisses et une pour les ordonnées. Le programme principal, se résumera donc à :
c <- -0.85+0.2i img <- calcul_img(c) afficher l'image avec imshowIl reste à écrire la fonction
calcul_img
et à tester le programme.
Vous devriez obtenir une image similaire à celle-ci :
La suite du travail va être réalisée à partir du programme traçant l'image en noir et blanc.
Il est prudent d'en faire une copie (du programme donnant l'image en monochrome) pour conserver une version fonctionnelle.
julia
renvoyait jusqu'à présent 0 ou 1, ce qui correspond aux deux couleurs utilisées. Nous allons donc modifier la fonction julia
de telle manière qu'elle puisse renvoyer un plus grand nombre de valeurs différentes.
Considérez le fait suivant : Pour un point donné $u_0$ qui conduit à une suite dont le module ne reste pas fini, nous nous en apercevons si le module d'un terme de la suite dépasse 2. Mais nous pouvons nous en apercevoir dans les premiers termes ou au contraire au bout du centième etc… Nous pouvons donc savoir si une suite a rapidement eu un terme dont le module dépasse 2. Autrement dit, ce qui nous intéresse lorsqu'un terme $u_i$ de la suite a un module qui dépasse 2, c'est : combien vaut $i$ à ce moment là?
julia
de telle manière que :- si au bout de 100 termes, tous les modules sont restés inférieurs à 2, elle renvoie la valeur 0
- si lors du calcul du terme numéro i, le module dépasse 2, la fonction renvoie un nombre strictement positif (peut être i ou i+1…)
julia
qui renvoie 0 si le point est dans l'ensemble, et qui renvoie un entier strictement positif sinon.
imshow
, produisez une image couleur de l'ensemble de Julia.matplotlib
ne contenait pas déjà certaines fonctionnalités.vmin
et vmax
de imshow
. Essayez de comprendre à quoi ils correspondent.Avant de passer à la suite, faites à nouveau un point avec l'encadrant de TP.Une fois le programme validé, vous rendrez le code correspondant.numpy
les fera en parallèle (ils sont faits en parallèle du point de vue de l'utilisateur, mais la machine, finalement, les fait séquentiellement, par contre toutes les optimisations possibles sont utilisées (raccourcis de calcul, utilisation de bibliothèques déjà compilées,…).
L'idée générale est de fournir une matrice contenant toutes les valeurs de $u_0$ pour lesquels on veut faire le calcul. Si cette matrice est M, le calcul se limite à : ''M <- M ** 2 + C (et les boucles, très lentes en Python, seront supprimées).
Voici les différentes étapes illustrées par une petite matrice 5×5.
Il est probable que ce qui suit soit nouveau pour vous, lisez attentivement et testez les exemples. Assurez-vous de les comprendre parfaitement.
Le premier objectif est d'obtenir une matrice de ce type, chaque case correspond à un pixel de l'image et contient la valeur de $u_0$ utilisée :
<code>
np.set_printoptions(precision=3)
-1.250-1.25j -0.625-1.25j 0.000-1.25j 0.625-1.25j 1.250-1.25j ]
[-1.250-0.625j -0.625-0.625j 0.000-0.625j 0.625-0.625j 1.250-0.625j]
[-1.250+0.j -0.625+0.j 0.000+0.j 0.625+0.j 1.250+0.j ]
[-1.250+0.625j -0.625+0.625j 0.000+0.625j 0.625+0.625j 1.250+0.625j]
[-1.250+1.25j -0.625+1.25j 0.000+1.25j 0.625+1.25j 1.250+1
</code>
Cette matrice est une combinaison linéaire de ces deux matrices :
<code>
-1.25 -0.625 0. 0.625 1.25 ]
[-1.25 -0.625 0. 0.625 1.25 ]
[-1.25 -0.625 0. 0.625 1.25 ]
[-1.25 -0.625 0. 0.625 1.25 ]
[-1.25 -0.625 0. 0.625 1.25
ET
-1.25 -1.25 -1.25 -1.25 -1.25 ]
[-0.625 -0.625 -0.625 -0.625 -0.625]
[ 0. 0. 0. 0. 0. ]
[ 0.625 0.625 0.625 0.625 0.625]
[ 1.25 1.25 1.25 1.25 1.25
</code>
La première de ces deux matrices est obtenue par produit (
np.dot) de :
<code>
| 1. |
| 1. |
| 1. |
| 1. |
| 1. |
PAR
| -1.25 -0.625 0. 0.625 1.25 |
</code>
et enfin, la première de ces deux matrices est obtenue avec la fonction
np.ones et la seconde en combinant
np.linspace, qui permet d'obtenir des valeurs réparties linéairement et
np.reshape qui permet de modifier la dimension d'une matrice (on pourrait utilise
np.transpose aussi).
Une fois cette matrice obtenue (supposons qu'elle s'appelle u), nous pouvons réaliser une itération sur
toutes les valeurs de la matrice ainsi;
<code>
u = u 2 + c
print(u)
-0.850+3.325j -2.022+1.762j -2.413+0.2j -2.022-1.363j -0.850-2.925j]
[ 0.322+1.762j -0.850+0.981j -1.241+0.2j -0.850-0.581j 0.322-1.363j]
[ 0.713+0.2j -0.459+0.2j -0.850+0.2j -0.459+0.2j 0.713+0.2j ]
[ 0.322-1.363j -0.850-0.581j -1.241+0.2j -0.850+0.981j 0.322+1.762j]
[-0.850-2.925j -2.022-1.363j -2.413+0.2j -2.022+1.762j -0.850+3.325j
</code>
Dans ce qui précède, si c est un scalaire (complexe), il est ajouté à chaque élément de la matrice
après la mise au carré.
Le test de sortie du disque de rayon 2 peut aussi être vectorisé très simplement :
<code>
mask = (abs(u) < 2)
print(mask)
False False False False False]
[ True True True True True]
[ True True True True True]
[ True True True True True]
[False False False False False
</code>
Encore mieux, cette matrice de booléens peut être utilisée pour filtrer l'accès aux cases d'une matrice.
Par exemple, si nous avons une matrice 5×5 contenant des valeurs aléatoires :
<code>
alea = np.random.randint(10, size=(5,5))
print(alea)
9 5 1 8 6]
[6 0 1 5 1]
[8 8 6 7 1]
[4 4 7 0 0]
[1 9 7 5 7
</code>
Nous pouvons ajouter 1 à tous les éléments de cette matrice qui correspondent à une valeur True dans le masque :
<code>
alea[mask] = alea[mask] + 1
print(alea)
9 5 1 8 6]
[7 1 2 6 2]
[9 9 7 8 2]
[5 5 8 1 1]
[1 9 7 5 7
</code>
Cette fonctionnalité est très puissante. Elle va nous permettre de faire le compteur d'itérations (la valeur de ce compteur donne la couleur d'un pixel), et aussi de n'itérer que sur les valeurs de la matrice u dont le module n'a pas encore dépassé 2.
Votre nouveau programme n'aura plus de fonction julia
ou conversion
. La fonction calcul_img
contiendra la plupart des opérations vectorisées.
Avant de passer à la suite, faite à nouveau un point avec l'encadrant.Si tout est fonctionnel (attendez qu'il soit validé), rendez le programme obtenu.
===== Lapin de Douady =====
Parmi les nombreux ensembles de Julia, certains ont des noms. C'est le cas du lapin de Douady qui vérifie la propriété suivante : la trajectoire de l'origine (la suite obtenue pour $u_0=0$) a pour longueur 3. Autrement dit, si $u_0=0$, $u_3=0$.
Calculez la valeur de $c$ qui permet de tracer le lapin de Douady. Mis en équation, vous obtenez un polynome de degré 3 ou 4 en $c$ dont vous pourrez chercher les racines avec NumPy. Tracez ensuite l'image correspondante.
Racines d'un polynôme avec NumPy import numpy
# création du polynôme x^2-3x+2
p=numpy.poly1d([1,-3,2])
# calcul des racines
numpy.roots(p)
# ou bien
p.roots
Faites valider l'image par l'encadrant.Une fois le résultat validé, rendez le code qui vous a permis de calculer la valeur de $c$.