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,
imshowfait correspondre la valeur max de l'image à la valeur max de la colormap (idem pour la valeur min).
En admettant que nous désirions tracer l'ensemble de Julia pour $u_0=a+ib$ avec $-2\le a\le 1$ et $-1.5\le b\le 1.5$, dans une matrice de taille 200x200, le changement de repère adéquat est :$$\left\{\begin{array}{l}a=\frac{i}{200}\times3-2\\b=\frac{j}{200}\times3-1.5\end{array}\right.$$
Vérifiez (par le calcul ou mentalement) que le point (0,0) est bien envoyé sur (-2,-1.5) (les valeurs mini pour a et b) et que le point (200,200) est envoyé sur (1,1.5) (les valeurs maxi pour a et b).
Dans votre programme, dans un premier temps, vous utiliserez le cadrage suivant : $-1.25\le a\le 1.25$ et $-1.25\le b\le 1.25$. Écrivez la fonction
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.Vérifiez que votre fonction
==== Tracé et programme principal ====
Nous savons maintenant faire les choses suivantes :
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
juliaprenant 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.
Avant de passer à la suite, faites le point avec l'encadrant afin de vérifier que vous avez correctement rempli le contrat. Une fois le travail validé, rendez l'intégralité du code realisé.
=====Version en couleur=====
Nous allons maintenant essayer d'obtenir des images en couleur, modifier la valeur de c et modifier éventuellement le cadrage.
==== Obtenir de la couleur ====
Le principe est simple: la fonction 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à?
Modifiez la fonction
Nous disposons maintenant d'une fonction 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.
En jouant sur la colormap et les valeurs vmin, vmax de
Voici un exemple d'image obtenue par ce principe :
imshow, produisez une image couleur de l'ensemble de Julia.
L'utilisation des colormap nous épargne le calcul qui permet de fabriquer les couleurs à partir d'un nombre unique et la manière dont elles se succèdent dans le dégradé. Ce serait une part importante du travail si
====Modification de la valeur de c ====
Nous avons fixé une valeur de $c$ de manière arbitraire. Vous pouvez changer cette valeur pour obtenir d'autres images (aidez vous de la cartographie de l'ensemble de Julia indiquée en début de sujet).
==== Modification du cadrage ====
Nous avons fixé les bornes de tracé de l'ensemble de Julia à -1.25, 1.25 en abscisse et en ordonnées. Vous pouvez modifier ces bornes pour obtenir un «zoom» sur certaines parties de l'image.
==== Résultats obtenus ====
Voici quelques échantillons de résultats possible.
matplotlib ne contenait pas déjà certaines fonctionnalités.
Essayez d'obtenir une belle image, avec une colormap de votre choix, votre propre valeur de c et un cadrage éventuellement différent de celui proposé en début de TP. NE mettez pas n'importe quoi pour les paramètres
===== Optimisation du calcul =====
Comme vous l'avez remarqué, le calcul n'est pas très rapide. Si nous voulions obtenir une image de taille conséquente… ce serait très long (essayez une image 400x400).
Pour accélérer le calcul, nous allons utiliser le principe de la vectorisation. En effet, les mêmes calculs sont menés sur chaque pixel de l'image. Utiliser la vectorisation consiste à fournir l'ensemble des points sur lesquels faire le calcul et 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$.
