Implémenter fractales de Julia

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 $c$ 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 $x_1$, $x_2$, $y_1$ et $y_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 $c$ et i. "$c$" 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 $Z_n+1$ (le "$+ 1$" s’applique à $n$ 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 à $Z$ (rappelez-vous, la suite pour $Z_n$) ;
  • nous entrons dans une boucle for qui bouclera i fois ;
  • nous appliquons la formule $z = z² + c$ ;
  • affectons la valeur ($Z_n$) à la liste z_list ($Z$) ;
  • 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 $Z$ pour $c = -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.

Réponse corrigée, en effet mauvaise interprétation et erreur dans le retour de l’itération.

Bonsoir,

Mon programme fonctionne correctement grâce à vos remarques, et il est disponible ici sur mon GitHub. Pour $c = -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 $c$ et ci sa partie imaginaire. Donc quand on entre cette commande, on a:

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.

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 :

Avec 200 itérations :

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()

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 ?

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.

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.

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.

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.

Nous parlons bien de la même chose. Je pensais à la base que np.vectorise faisait de la compilation à la volée par exemple comme numba ou un binding avec un executable plus efficace.

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.