Open Source Your Knowledge, Become a Contributor
Technology knowledge has to be shared and made accessible for free. Join the movement.
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 .
Résumons la définition de l'ensemble de Mandelbrot : Un complexe
Pour déterminer son aire, nous allons utiliser une méthode dite de Monte Carlo : Prendre des points au hasard dans le rectangle
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 :
Dérouler pour voir le programme
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 :
Dérouler pour voir les modifications
Donnons quelques explications :
-
Pour la fonction
mandelbrot(c)
, elle va agir directement sur un vecteur contenant tous lesc
choisis aléatoirement. Ainsi, au lieu de renvoyer juste 0 ou 1 selon que lec
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 lesc
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 commez*z
sont faites coordonnées par coordonnées tout commenp.abs(z)
appliqueabs
à chaque coordonnée). -
Pour la fonction
monte_carlo()
, on ne fait plus une boucle pour créer à chaque fois unc
mais on crée un vecteur de tous lesc
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 fonctionsmandelbrot
etmonte_carlo
- Admirez : vous obtenez le même résultat en moins de 2 sec...
En voici la preuve :
Dérouler pour voir les modifications
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 :
Dérouler pour voir les modifications
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
on calcule en réalitéz 2 etx 2 − y 2 (si on note2 x y ). De plus on vérifie à chaque fois siz = x + y i or ce calcul est très gourmand car il faut calculera b s ( z ) ≥ 2 . Il est donc plus sage de vérifier simplement quex 2 + y 2 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.x 2 + y 2 ≥ 4 -
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 calculerMAX_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 appartient à cette cardioide ou ce cercle avant de calculer les termes de la suite. On pourra trouver les équations de ces ensembles ici.c
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 :
Dérouler pour voir les modifications
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).