Outils pour utilisateurs

Outils du site


tp:python:julia2a

📓 Ensemble de Julia et vectorisation

sol

Nous proposons ici de tracer des ensembles de Julia en noir et blanc, puis en couleur. Dans un premier temps nous ferons un calcul non optimisé pour la vitesse, puis nous verrons comment vectoriser le calcul en utilisant numpy

Références

Version monochrome

Rappels sur les ensembles de Julia

Notons tout d'abord qu'il existe plusieurs ensembles de Julia. Pour chaque nombre complexe $c$ choisi, il y a un ensemble. Nous noterons $J©$ l'ensemble des points du plan appartenant à l'ensemble de Julia paramétré par $c$.

L'ensemble de Julia $J©$ étant un ensemble de points du plan, nous pouvons le voir comme un ensemble de nombres complexes. Les nombres complexes appartenant à $J©$ sont les nombres $u_0$ tels que la suite : $$\forall n>0, u_n=u_{n-1}^2+c$$ a un module qui reste fini.

Dans un premier temps, nous tracerons l'ensemble de Julia correspondant à $c=-0.85 + 0.2 i$. Il faut que cette valeur soit facile à modifier par la suite\ : gardez ceci à l'esprit tout au long de la construction du programme (car nous testerons d'autres valeurs).

On ne peut pas déterminer sur une machine si les termes de la suite précédente ont un module qui reste fini ou non en calculant simplement un grand nombre de termes. Il faudra donc faire quelques concessions à l'exactitude. Une bonne approximation consiste à calculer les 100 premiers termes de la suite et à utiliser cette propriété (valable si le module de c est strictement inférieur à 2) :

S'il existe $N$, entier, tel que le module de $u_N$ soit strictement plus grand que 2, alors la suite n'a pas un module qui reste fini.

Ainsi, si nous trouvons un terme de module supérieur à deux, nous pourrons conclure que le point n'est pas dans l'ensemble de Julia. Si au bout du 100ème terme, aucun n'a dépassé 2 en module, on conclura (parfois faussement) que le point est dans l'ensemble de Julia.

Première fonction

En vous aidant des indications données précédement, et de la manière de manipuler les nombres complexes en Python (voir [[stu:python:python_notes|Aide mémoire / Notes sur Python 3]]), écrivez une fonction julia(u,c) qui indique en renvoyant 0 ou 1 si u est dans ''J(c). <code python > def julia(u,c) : “”“ Indique si le complexe u appartient (1) ou non (0) à J© ”“” pass </code>

Testez votre fonction en comparant vos résultats avec ceux-ci, obtenus pour $c=-0.85+0.2i$ :

  • les 3 points : 0, 0.3j et 1 sont dans l'ensemble de Julia
  • les 2 points 0.3 et 1+j ne sont pas dans l'ensemble de Julia

Donnez un copié/collé des tests effectués.

==== Représentation de l'image ==== La première chose à décider est le nombre de points que nous allons calculer. Restons modestes dans un premier temps, et calculons 200×200 pixels. Les valeurs de ces pixels (1 ou 0 selon qu'on est ou pas dans l'ensemble de Julia) seront stockés dans une matrice numpy : <code python> import numpy as np image = np.zeros1)

1)
200,200), dtype=np.uint8) </code> Pour tester l'affichage, nous allons mettre aléatoirement quelques pixels à 1 dans cette image :
for i in range(200*200 // 4):
    x = random.randint(0, 199)
    y = random.randint(0, 199)
    image[x, y] = 1
Puis 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).
==== Changement de repère ==== Ce qui précède nous permet, étant donné un nombre complexe $u_0$, de savoir s'il est ou non dans l'ensemble de Julia. De plus, nous savons représenter une matrice de valeurs (dont les coordonnées sont entières) à l'écran. Il reste maintenant à faire correspondre les nombres complexes manipulés aux coordonnées des valeurs dans la matrice numpy. Ceci s'apparente à un simple changement de repère. Nous devons décider quelles sont les limites (parties réelles et imaginaires) de $u_0$, et en déduire le changement de repère (mise à l'échelle et translation).
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 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.
==== Tracé et programme principal ==== Nous savons maintenant faire les choses suivantes :
  • 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.
Nous pouvons maintenant proposer un algorithme pour le calcul de l'image, qui crée une matrice en prenant en paramètre la valeur de c (paramètre de l'ensemble de Julia) :
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 img
Attention à 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 imshow
Il reste à écrire la fonction calcul_img et à tester le programme. Vous devriez obtenir une image similaire à celle-ci :  Ensemble de Julia 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 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…)
Nous disposons maintenant d'une fonction 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 imshow, produisez une image couleur de l'ensemble de Julia.
Voici un exemple d'image obtenue par ce principe :
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 matplotlib ne contenait pas déjà certaines fonctionnalités.
====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.
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 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.
===== 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 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$.
tp/python/julia2a.txt · Dernière modification: 2021/05/04 15:12 (modification externe)