Trying to get Julia fractal to run fast in python. Mainly using numba is the standard way possible.
Soit la suite $$ \left(z^2_i +c\right)_i,n\in \mathbb{N}. $$ l'ensemble de Julia qualifie les points définissant la frontière dans l'espace des \(z_0\) entre la zone où elle est bornée et celle où elle diverge. Cette frontière est fonction de la variable complexe \(c\).
Afin de visualiser l'ensemble fractal, nous allons effectuer une cartographie du domaine des $z_0$ et tester la convergence de la suite en tout ces points. C'est à dire rechercher les ensembles $$ Q_c^{k} = \left\{z_0 : |z_k| \le R(c) \right\} $$ pour chaque niveau de convergence \(k\), chacun sera associé à un niveau de gris dans l'image finale. Pour calculer le niveau de convergence \(k\) d'un \(z_0\), nous pouvons écrire une fonction du style:
def suite(z0, c=0.0, kmax=25):
zk = z0
Rc = max(2, abs(c))
k = 0
while abs(zk) <= Rc and k < kmax:
zk = zk**2 + c
k += 1
return k
où $k_{\text{max}}$ est un garde-fou pour éviter de rester coincé dans le cas où la suite n'est pas bornée.
Une fois cette suite définie, il nous reste à cartographier tous les $z_0$ potentiellement intéressants
import numpy as np
resolution = 256
realpart = np.linspace(-2, 2, resolution)
imagpart = np.linspace(-2, 2, resolution)
z0 = np.vstack(np.meshgrid(realpart, imagpart))
L'axe des abscisses représente donc les parties réelles de $z_0$ croissantes et l'axe des ordonnées, les parties imaginaires.
Nous pouvons enfin en finir en calculant le $k$ de chaque pixel de cette carte des $z_0$.
Ce résultat n'est pas exactement ce à quoi on s'attend quand on parle de fractal, la clé est de modififier la valeur de \(c\).
Il faudra aussi modifier la résolution et les valeurs de \(k_{max}\) pour bien visualiser le contraste et les détails.
Numba is a python package that allows you to compile just in time your functions. This has the benefit of increasing speed of execution for a rather small startup fee (the compilation time). I really like numba for its ease of integration. If you have a small algorithm to write (that is not in numpy, scipy) and want to make it decently fast: give numba a try. The next step would be Cython.
@njit
def suite(z0, c=0.0, kmax=25):
zk = z0
Rc = max(2, abs(c))
k = 0
while abs(zk) <= Rc and k < kmax:
zk = zk**2 + c
k += 1
return k
This is not a great improvement. To understand this we need to know that the passage from a python function to a compiled numba function has a cost. This overhead can be nullified by integrating the loop inside the numba space of our code. To do so, we create a function that takes in both the table of $z_0$ and the $k$ table.
@njit
def compute1(z0s, ks):
c = 0.375 + 0.2j
for i in range(len(z0s)):
z0 = z0s[i,0] + 1j * z0s[i,1]
ks[i] = suite(z0, c=c)
Try this and you will see drastically faster speeds! The question is now: can we go even faster? Hopefully yes, even a lot faster! We can for example get rid of the `meshgrid` and get rid of the function overhead. The main gain here comes from getting rid of the `meshgrid`.
@njit
def compute3(ks, xmin, xmax, resolution, kmax=25):
dx = (xmax-xmin) / resolution
c = 0.375 + 0.2j
Rc = max(2, abs(c))
for i in range(resolution):
realpart = -2 + i * dx
for j in range(resolution):
imagpart = -2 + j * dx
zk = realpart+1j*imagpart
k = 0
while abs(zk) <= Rc and k < kmax:
zk = zk**2 + c
k += 1
ks[i + resolution + j] = k
Numba also allows for multithreading, we can use this by addind `parallel=True` to the decorator and turning the outer loop in a `prange`. For `compute3` and `compute2` this leads to $\%$ and $\%$ improvements. To see what is happening, you can write
ks[i + resolution * j] = get_thread_id()
so that the fractal is ... the id of the thread that computes the specific pixel. This leads to this map where each color patch represents the work done by each of your thread. `get_thread_id()` is a numba function that can be imported like so
from numba.np.ufunc.parallel import get_thread_id
home