Apprendre Python dans le secondaire

Exemple de vectorisation : Aire de l'ensemble de Mandelbrot

Le but de cette page est de présenter sur un exemple un aspect propre à python : La vectorisation. En effet, python a de nombreux avantages mais en contrepartie, c'est un langage qui est très lent (sur les boucles par exemple). Cependant de nombreux modules, pour gagner en rapidité, sont programmés directement dans un autre langage (le C). C'est le cas du module numpy qui permet, entre autres, des calculs sur des vecteurs/tableaux très rapides. On appelle donc vectorisation le fait d'utiliser des vecteurs (arrays) Numpy pour éviter un maximum de boucles python.

Une fois cette vectorisation faite, on utilisera le module Numba qui permet d'accélérer encore davantage les calculs.

Cette page est une simple présentation sans exercice. Il suffira de valider pour voir le résultat. Vous pouvez bien sûr modifier ces exemples pour voir ce que donnent vos modifications. Les constantes seront toujours les mêmes pour bien voir le gain en rapidité. vous pourrez donc les modifier pour affiner la précision du résultat.

Présentation de l'exemple : Calcul de l'aire de l'ensemble de Mandelbrot

On peut trouver une présentation de l'ensemble de Mandelbrot sur Wikipédia très complète. Une remarque : une page est dédiée dans ce recueil sur les ensembles de Mandelbrot et Julia mais avec une définition un peu différente et moins classique. Nous revenons ici à la définition à partir de la fonction $$ f(z)=z^2+c$$.

Résumons la définition de l'ensemble de Mandelbrot : Un complexe $$ c$$ est dans l'ensemble de Mandelbrot si et seulement si la suite définie par $$ z_0=0$$ et $$ z_{n+1}=z_n^2+c$$ ne tend pas vers l'infini. Voici sa représentation :

Pour déterminer son aire, nous allons utiliser une méthode dite de Monte Carlo : Prendre des points au hasard dans le rectangle $$ -2\le x\le 1$$ et $$ -1\le y\le 1$$ et compter le nombre de points qui font partie de l'ensemble de Mandelbrot. On peut trouver un exemple d'utilisation de cette méthode pour le calcul d'intégrales ici.

Dans la pratique, si un terme de la suite a un module plus grand que 2, cela suffit pour savoir si la suite tend vers l'infini. Voici à quoi ressemblerait un programme directement inspiré des définitions pour calculer cette aire :

Utilisation de Numpy pour vectoriser

L'idée de la vectorisation est de remplacer une (ou plusieurs) boucle par des calculs sur un tableau/ vecteur numpy. Par exemple ici, au lieu de choisir une par une les valeurs aleatoires de c, on va créer directement un tableau contenant tous les c choisis aléatoirement et faire les calculs directement en utilisant le tableau. C'est a priori plus compliqué et demande plus de ressources mais les boucles en python sont tellement lentes que faire ainsi permet de gagner en vitesse. Voici le programme :

Donnons quelques explications :

• Pour la fonction mandelbrot(c), elle va agir directement sur un vecteur contenant tous les c choisis aléatoirement. Ainsi, au lieu de renvoyer juste 0 ou 1 selon que le c est dans l'ensemble ou pas, la fonction va renvoyer un vecteur de 0 ou de 1 correspondant à la même chose mais pour tous les c d'un coup.
  Toute la fonction est donc repensée pour faire tous les calculs d'un coup sur le vecteur et en utilisant les propriétés des tableaux numpy (par exemple les opérations comme z*z sont faites coordonnées par coordonnées tout comme np.abs(z) applique abs à chaque coordonnée).

• Pour la fonction monte_carlo(), on ne fait plus une boucle pour créer à chaque fois un c mais on crée un vecteur de tous les c qu'on va considérer.

Interlude : Utilisation de numba

Attention, ce qui suit va ressembler à un tour de magie :

• Prenez le premier code (celui qui met presque 30 sec)
• Rajoutez from numba import jit
• Rajoutez @jit sur la ligne précédent les fonctions mandelbrot et monte_carlo
• Admirez : vous obtenez le même résultat en moins de 2 sec...

En voici la preuve :

Si j'ai bien compris, Numba traduit plus ou moins le langage python pour directement traduire notre programme en code machine. Ce qui le rend, dans notre cas, aussi rapide que des langages comme le C. Malheureusement la magie n'opère pas toujours car certaines fonctions sont propres au langage python (les listes en compréhension par exemple) mais en tout cas cela marche très bien avec les fonctions de bases et numpy pour peu que cela soit fait proprement.

Vectorisation en utilisant numpy et numba

On peut encore gagner un peu en vitesse en combinant l'utilisation des tableaux numpy et numba. Pour cela il va falloir vectoriser notre fonction mandelbrot pour qu'elle s'applique à chaque élément d'un tableau numpy directement. De plus, cela nous permet de préciser les types des entrées et des sorties ce qui permet d'économiser encore du temps. Voici à quoi cela peut ressembler :

Peaufinage

On peut améliorer encore un peu notre code en utilisant deux remarques :

• La première est classique lorsqu'on veut optimiser des calculs : Eviter les calculs inutiles. On ne s'en rend pas vraiment compte en regardant le programme précédent mais on choisit 1 million de points au hasard et pour chacun, il faut calculer dans le pire des cas 500 termes de la suite. Cela fait beaucoup de calculs mais en plus, certains sont faits plusieurs fois. En effet, lorqu'on calcule $$ z^2$$ on calcule en réalité $$ x^2-y^2$$ et $$ 2xy$$ (si on note $$ z=x+yi$$). De plus on vérifie à chaque fois si $$ abs(z)\ge 2$$ or ce calcul est très gourmand car il faut calculer $$ \sqrt{x^2+y^2}$$. Il est donc plus sage de vérifier simplement que $$ x^2+y^2\ge 4$$ et on voit ainsi qu'on a intérêt à ne pas passer par les complexes mais plutôt faire soi-même les calculs des parties réelles et imaginaires.

• La seconde est que ce qui coûte cher, ce sont les valeurs de c qui sont dans l'ensemble de Mandelbrot car il faut calculer MAX_ITER termes de la suite pour constater qu'elle ne diverge pas. Or, si on observe l'ensemble de Mandelbrot, on constate qu'un grosse partie de sa surface provient d'une cardioide et du cercle situé à sa gauche. L'idée est donc de vérifier d'abord si $$ c$$ appartient à cette cardioide ou ce cercle avant de calculer les termes de la suite. On pourra trouver les équations de ces ensembles ici.

Si on modifie le programme précédent pour y ajouter ces deux remarques, on obtient le programme suivant qui fait tomber à moins de 0.2 sec le calcul qui à l'origine prenait autour de 30 sec :

Conclusion

Vous pouvez modifier les valeurs de MAX_ITER et NB_POINTS pour affiner la précision prendre de trop grandes valeurs (le temps de calcul maximal permis sur ce site est 30 sec). Par exempe MAX_ITER = 10 000 et NB_POINTS = 20 000 000. On n'est finalement pas si loin de la valeur de l'aire de l'ensemble de Mandelbrot donnée sur Wikipédia est de 1,50659177 ± 0,00000008.

Pour aller plus loin en précision, il faudrait non seulement donner des valeurs plus grandes à MAX_ITER et NB_POINTS mais aussi prendre en compte davantage de chiffres après la virgule pour éviter les erreurs d'arrondi. On peut trouver sur la page anglaise de wikipédia d'autres optimisations possibles (comme l'utilisation d'outils de la théorie des pertubations).