Implémenter fractales de Julia

Comment afficher l'une des fractales de Julia tout en, si possible, modifiant le code de Mandelbrot

a marqué ce sujet comme résolu.

Bonjour,

J’ai hier développé un simple programme qui me permettait d’afficher l’ensemble de Mandelbrot en Python. Aujourd’hui, j’ai envie de reprendre ce programme afin de le transformer en quelque chose qui m’afficherai selon cc une fractale de Julia. Alors, pour premier jet, même si ça ne marche pas, voilà ce que ça donne:

import numpy as np
import matplotlib.pyplot as plt


def create_complex_matrix(x1, x2, y1, y2, quality):
    x = np.linspace(x1, x2, int((x2 - x1) * quality))
    y = np.linspace(y1, y2, int((y2 - y1) * quality))
    return x[np.newaxis, :] + y[:, np.newaxis] * 1j

def julia(c, i):
    z = c
    z_list = []
    for _ in range(i):
        z = z ** 2 + c
        z_list.append(z)
    return z_list

plan = create_complex_matrix(-2, 0.5, -1.5, 1.5, quality=1000)
mask = julia(complex(-0.1, 0.651), i=30)

plt.imshow(plan[mask], cmap="binary")
plt.show()

Ce programme, fonctionne de la sorte:

  • Nous importons la librairie numpyainsi que matplotlib.pyplot qui nous servirons respectivement à créer des matrices et à afficher ces matrices.
  • Nous définissons une fonction create_complex_matrix qui prend en entrée x1x_1, x2x_2, y1y_1 et y2y_2. A partir de ceux-ci il créé deux matrice unidimensionnel et ressort la "fusion" de ces dernières en une matrice bidimensionnel qui nous permet donc d’avoir une matrice ayant des nombres flottants, et négatif (ce qui n’est pas permis nativement).
  • Nous définissons la fonction julia qui prend en argument cc et i. "cc" correspond au nombre complexe à partir duquel nous voulons que la fractale s’organise. i n’est qu’une simple valeur qui déterminera le nombre de fois que nous ferons Zn+1Z_n+1 (le "+1+ 1" s’applique à nn mais je n’ai pas réussi à le mettre à l’échelle), donc, déterminera la complexité de la fractale rendue (puisque il détermine en soit le nombre de pixel de la fractale). Le fonctionnement interne de la fonction est relativement simple:
    • nous affectons c à z et non pas 0 pour gagner un tour ;
    • nous définissons la liste z_list qui correspond à ZZ (rappelez-vous, la suite pour ZnZ_n) ;
    • nous entrons dans une boucle for qui bouclera i fois ;
    • nous appliquons la formule z=z²+cz = z² + c ;
    • affectons la valeur (ZnZ_n) à la liste z_list (ZZ) ;
    • puis retournons la liste z_list.

Ensuite, nous créons une matrice bidimensionnelle qui correspond aux normes de la surface nécessaires nativement à Mandelbrot, créons une variable mask qui contient la suite ZZ pour c=0.1+0.651ic = -0.1 + 0.651i (pris arbitrairement) et enfin, appliquons se masque à la matrice bidimensionnelle plan afin de n’afficher uniquement les nombres de ce dernier qui matcheraient avec le masque.

Or ce programme ne fonctionne pas et renvoie d’ailleurs une erreur:

    plt.imshow(plan[mask], cmap="binary")
IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

En plus de ça, je me doute qu’il à un pépin dans la logique de mon code, et j’attend vos retours quand au remaniement de ce programme.

Cordialement.

Update

L’algorithme de Julia diffère de celui de Mandelbrot. Dans mandelbrot, la constante c est la position dans le plan complexe. Dans Julia c’est le z0 initial qui est la position dans le plan complexe et c est une constante à part entière. D’autres part, on veut retourner l’itération maximale qui garde la norme de z dans la borne qu’on a fixée.

@np.vectorize # pour itérer sur chaque élément du plan complexe z0
def julia(z0, c, i):
    z = z0
    for it in range(i):
        z = z*z + c
        if (abs(z) > 2) : break # on stop l'algorithme si on sort de notre borne
    return it # on retourne l'itération (plus it est grand plus ça diverge lentement)

Comme alternative à la méthode des masques, on peut appliquer la fonction julia (grâce à @np.vectorize) à tout les point du plan complexe (la matrice plan). Ici on aura pas de résultats binaires mais un gradient d’intensité et donc on pourra colorer l’image en fonction de la vitesse de divergence.

Pour créer votre plan complexe vous pouvez aussi utiliser np.meshgrid() qui permet de faire de l’algèbre multidimensionnel et retourner des matrices :

X,Y = np.meshgrid(x,y)
M = X + 1j*Y

Vous trouverez différentes implémentations ici

+1 -0

C’est dommage de cracher une solution toute faite, qui plus est fausse. :-°

@np.vectorize # pour itérer sur chaque élément du plan complexe z0
def julia(z0, c, i):
    z = z0
    for _ in range(i):
        z = z*z + c
        if (abs(z) > 2) : return 0 # on stop l'algorithme si ça commence à diverger
    return i # on retourne l'itération si ça converge (la norme sera la vitesse de convergence)

Retourner le nombre d’itérations (soit tout le temps 30, ou 0 sur les points qui divergent) ne donnera pas une mesure de la vitesse de convergence. Ce qu’on veut plotter est une vitesse de divergence d’ailleurs…


A partir de ceux-ci il créé deux matrice unidimensionnel et ressort la "fusion" de ces dernières en une matrice bidimensionnel qui nous permet donc d’avoir une matrice ayant des nombres flottants, et négatif (ce qui n’est pas permis nativement).

Par "pas nativement", tu parles des indices flottants ? La matrice que tu as construit ne contourne aucunement ceci, elle est remplie de sorte que z[i, k] = x(k) + y(i) * 1j. Aucun langage (ni bibliothèque) un minimum sain ne va proposer l’indexage par des flottants. C’est ultra casse-gueule à cause des problèmes des calculs flottants qui ont une précision limitée.

et enfin, appliquons se masque à la matrice bidimensionnelle plan afin de n’afficher uniquement les nombres de ce dernier qui matcheraient avec le masque.

Les masques numpy ne fonctionnent pas comme ça, ce sont des tableaux de booleans qui disent si l’élément doit être montré ou non. Ce que tu voudrais faire serait équivalent à plan[plan==mask] (avec mask ne contenant que la dernière itération, pas toute la liste). Ce serait une façon logique d’afficher les points fixes (donc l’ensemble de Julia lui même), mais en fait ça ne peut pas marcher en pratique à cause de la précision finie des nombres flottants et de l’échantillonage discret du plan qui est fait. C’est pour ça qu’à la place, soit on fait un plot binaire des points qui ont divergés ou non (ce que tu as fait pour Mandelbrot), soit la vitesse de divergence des points qui ne sont pas dans l’ensemble, pour avoir une idée de la "distance" à laquelle on est de l’ensemble de Julia (et puis en bonus c’est joli comme ça).

+1 -0

Bonsoir,

Mon programme fonctionne correctement grâce à vos remarques, et il est disponible ici sur mon GitHub. Pour c=0.4+0.6ic = -0.4 + 0.6i, voici la commande à entrer pour générer la fractale associé, car oui, j’ai armé mon programme d’une CLI:

python /path/to/julia/run.py 1000 27 -0.4 0.6 -c "inferno"

Ou plus explicitement:

... <density-pixel> <detail-density> <cr> <ci> [-c THEME]

<cr> représente la partie réelle de cc et ci sa partie imaginaire. Donc quand on entre cette commande, on a:

Fractale de Julia pour c = -0.4 + 0.6i
Fractale de Julia pour c = -0.4 + 0.6i

Le problème, comme vous pouvez le voir, c’est qu’il n’y a aucun détail intérieur à la figure: celui-ci est juste remplie d’une couleur. Les seuls détails se trouvent en extérieurs, et avec n’importe quelle autre fractale, ces détails ne se montrent pas, et la cause m’échappe.

Pouvez-vous m’aider à résoudre ce problème ?

Cordialement.

+0 -0

Ton problème n’est qu’un effet du choix de la constante c. Suivant la constante c que tu choisis choisi, la fractale de Julia associé va beaucoup varier. Pour tous les valeurs de c qui correspondent à des points situés à l’intérieur de l’ensemble de Mandelbrot, les fractales de Julia associés sont entièrement connecté et tu vas te retrouver avec des zones unis. À l’inverse, pour tous les autres valeurs de c, les fractales de Julia associés ne sont composés que de points isolés et tu n’auras donc aucune zones unis.

Quand tu regardes l’ensemble de Mandelbrot, chaque pixel correspond à un c particulier, donc tu peux générer un ensemble de Julia pour chaque point d’une image d’un ensemble de Mandelbrot. En faisant le lien entre les deux, tu pourra mieux comprendre comment choisir un c qui te permettra d’avoir un ensemble de Julia qui va plus corresponde à ce que tu cherches.

En l’occurrence, c’est juste une question du nombre d’itérations que tu fais pour tester si un point diverge ou non. Avec 50 itérations :

image.png
image.png

Avec 200 itérations :

image.png
image.png

Le code, avec une approche un peu différente, pour ceux que ça intéresse

from dataclasses import dataclass

from numpy.typing import NDArray
import numpy as np
from PIL import Image


@dataclass
class ComplexPlane:
    xleft: float
    xright: float
    yleft: float
    yright: float

    def xspan(self) -> float:
        return self.xright - self.xleft

    def yspan(self) -> float:
        return self.yright - self.yleft

    def create(self, resolution: int) -> NDArray[np.complexfloating]:
        x_s = np.linspace(self.xleft, self.xright,
                          int(self.xspan() * resolution))
        y_s = np.linspace(self.yleft, self.yright,
                          int(self.yspan() * resolution))
        return x_s[:, np.newaxis] + y_s[np.newaxis, :] * 1j


@dataclass
class JuliaSetDivergence:
    c_0: complex
    threshold: float = 2.0
    n_iterations: int = 50
    resolution: int = 1000

    def over(self, plane: ComplexPlane) -> NDArray[np.floating]:
        z_s = plane.create(self.resolution)
        div = np.full_like(z_s, self.n_iterations, dtype=np.uint32)
        for i in range(self.n_iterations):
            z_s = np.square(z_s) + self.c_0
            diverged = np.abs(z_s) >= self.threshold
            div[diverged] = np.uint32(i)
            z_s[diverged] = np.nan
        return div / self.n_iterations


def main() -> None:
    plane = ComplexPlane(xleft=-2.0, xright=2.0, yleft=-2.0, yright=2.0)
    julia_div = JuliaSetDivergence(
        c_0=complex(-0.4, 0.6),
        resolution=500,
        n_iterations=200,
    )
    div_map = julia_div.over(plane)

    im = Image.fromarray((div_map * 255).astype(np.uint8))
    im.save("plot.tiff")


if __name__ == "__main__":
    main()
+1 -0

Mon programme fonctionne correctement grâce à vos remarques, et il est disponible ici sur mon GitHub.

b4b4

En lisant ton code je me demande à quoi sert la classe Julia : tu ne l’instancies jamais et elle ne contient donc que deux méthodes statiques (qui ne sont pas définies comme telle, ce qui pourrait poser problème).
Pourquoi ne pas simplement définir les fonctions à la racine du module julia ?

Alors, premièrement j’instencie belle et bien cette classe dans le fichier run.py à la racine du projet. Et oui en effet j´aurai très bien pu ne pas mettre de classe, mais je trouve que ça organise les choses d’une certaine manière.

Et il y a quelque chose que je ne comprend pas bien: les images générés par @adri1, proviennent elles de mon programme ou bien d’ailleurs ? Auquel cas, si il ne provient pas de mon code, le soucie vien-il du nombre d’itération, que j’ai intentionellement limité à 30 à cause de la puissance de mon ordinateur ?

+0 -0

Non tu ne l’instancies pas (la classe n’est jamais appelée et tu n’as donc pas d’instance). Si tu le faisais tes appels de méthodes échoueraient (puisqu’elles ne sont pas définies comme statiques) :

>>> julia = Julia()
>>> julia.create_complex_matrix(-4, 4, -4, 4, 1000)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: Julia.create_complex_matrix() takes 5 positional arguments but 6 were given

En Python les modules sont déjà une bonne manière d’organiser les choses. Il n’est donc pas utile d’ajouter une classe coquille vide, à part si tu as bien sûr un état à conserver (ce qui se fait dans une instance de la classe donc).

Aussi tu gagnerais à lister tes dépendances (numpy et matplotlib en l’occurrence) dans un fichier requriements.txt par exemple pour faciliter l’installation.

Et il y a quelque chose que je ne comprend pas bien: les images générés par @adri1, proviennent elles de mon programme ou bien d’ailleurs ? Auquel cas, si il ne provient pas de mon code, le soucie vien-il du nombre d’itération, que j’ai intentionellement limité à 30 à cause de la puissance de mon ordinateur ?

b4b4

Mon programme, que j’ai mis en fin de message. Oui, le soucis vient du nombre d’itérations, c’est ce que je dis et montre dans mon message… Si ton ordi est un peu lent, tu pourrais diminuer le nombre de pixels de l’image pour vérifier qu’augmenter le nombre d’itérations permet effectivement de résoudre plus de détails.

Salut,

Par curiosité, j’ai fait un petit test de performances pour voir si utiliser np.vectorize est pertinent par rapport à une solution qui passe son temps dans numpy quitte a faire un peu de travail pour rien. Le résultat est sans appel (pun intended), utiliser np.vectorize est hyper-lent. C’est normal, numpy se contente de faire une boucle sur la fonction Python, donc on perd beaucoup de temps à construire les objets Python, les manipuler, et appeler la fonction elle-même.

Pour l’histoire, je génère une image de 4000x4000 pixels avec 100 itérations. Je compare en gros cette implémentation :

def divergence(z_0: complex, c_0: complex, threshold: float, itermax: int):
    for i in range(itermax):
        z_0 = z_0**2 + c_0
        if abs(z_0) >= threshold:
            return i
    return itermax

div_vect = np.vectorize(divergence, otypes=[np.uint32])
div = div_vect(z_s, c_0, threshold, itermax)

Avec celle-ci:

div = np.full_like(z_s, itermax, dtype=np.uint32)
for i in range(itermax):
    z_s = np.square(z_s) + self.c_0
    diverged = np.abs(z_s) >= self.threshold
    div[diverged] = np.uint32(i)
    z_s[diverged] = np.nan

Le code à base de np.vectorize met 139 secondes à faire les calculs. Le code effectivement factorisé (mais qui fait des calculs pour rien) met seulement 19 secondes ! C’est un exemple de plus que faire un truc idiot point de vue algorithmique (plein de calculs sur des np.nan sur les points qui ont déjà divergés) peut être largement plus rapide à cause des préfacteurs immenses côté Python (et d’un autre côté, numpy fait probablement des trucs pas trop idiots pour bien vectoriser les calculs).

Un autre truc qui prend du temps dans le code de @b4b4 est l’utilisation de pcolormesh (quelques secondes) plutôt que construire une image directement avec PIL (coût négligeable).

Intéressant ! Je viens de lire ça, np.vectorize ne semble être qu’une macro qui ne fait rien d’autre qu’une boucle for. Je n’ai jamais essayé numba mais c’est peut-être une piste pour gagner en efficacité. @b4b4 comme l’algorithme ne dépend que du pixel considéré, tu peux paralléliser le code (avec numba justement ou en faisant du multi-threading), ça peut être intéressant pour voir à quel point tu peux augmenter la résolution de l’image pour un même temps de traitement. Et c’est aussi un bon exercice si tu n’as jamais fait de multi-threading.

+0 -0

Intéressant ! Je suppose que la vectorisation devient effectivement efficace à partir d’un rang et/ou d’une taille de tableau donnée.

Quand tu dis vectorisation, tu parles de np.vectorize, ou de la version qui fait juste des manipulations sur les np.ndarray ? Parce qu’en fait, ironiquement, vectorize ne vectorise pas les opérations au sens où on l’entend classiquement en HPC (et qui idéalement donne accès aux instructions CPU pour faire ça efficacement).

np.vectorize se contente de faire des boucles sur les tableaux passés en arguments et appeler la fonction Python sur chaque élément. D’un autre côté, les opérations sur les tableaux numpy sont effectivements des opérations sur des tableaux contigus qui peuvent exploiter le cache et les instructions processeurs qui vont bien (et on pourrait même faire encore mieux en compilant pour la machine sur laquelle on tourne plutôt qu’un dénominateur commun). C’est pour ça que np.vectorize est hyper-lent par rapport à des opérations sur un array numpy, et à opérations identiques ce sera toujours plus lent. Vu le facteur 7 que se prend np.vectorize dans le cas présent, sa seule chance pour concurrencer avec une implémentation sur les arrays est d’avoir beaucoup plus d’itérations (ce qui en pratique n’est utile que si on zoome à mort) pour que le gâchis d’opérations commence à faire sentir.

Sur une fonction aussi simple, j’imagine que numba s’en sortirait pas trop mal en effet, il arriverait à vectoriser et virer tout appel à l’API Python.

Dernièrement j’ai pas mal joué avec Rust + PyO3 + ndarray et c’est vraiment pas mal pour exposer une API Python au-dessus d’un code vectorisé proprement.

Connectez-vous pour pouvoir poster un message.
Connexion

Pas encore membre ?

Créez un compte en une minute pour profiter pleinement de toutes les fonctionnalités de Zeste de Savoir. Ici, tout est gratuit et sans publicité.
Créer un compte