L
Lazaro / Aug 20 2019

Aggregating values to the Mandelbrot and Julia sets

by Lazaro Alonso, find me on twitter as @LazarusAlon

Since the mandelbrot and julia sets are well known fractal examples in the literature, here only the computational side will be posted.

The traditional approach for plotting this sets is by pre-pixeling the coordinates( on the complex plane) and then computing their values(iteration process). But now, thanks to datashader we can directly perform the iteration test on the complex plane coordinates and then aggregate the iteration outputs into the corresponding pixel.

Getting ready

First let's see the python version installed and conda install some packages. 

import sys; sys.version.split()[0]
0.4s
Python
Python 3
'3.6.8'
conda install datashader
47.4s
Bash in Python
Python 3
conda install numba
9.0s
Bash in Python
Python 3
conda install colorcet
8.8s
Bash in Python
Python 3

Now, let's work... we are ready to import all that we need.

import numpy as np, pandas as pd, colorcet as cc
import datashader as ds, datashader.transfer_functions as tf
from numba import jit
from colorcet import palette
from functools import partial
from datashader.colors import colormap_select
from numpy import log, arange
from datashader.colors import inferno, viridis
background = "black"
cm = partial(colormap_select, reverse=(background!="black"))
4.5s
Python
Python 3

Mandelbrot set

The Mandelbrot set, M, is the set of values of c in the complex plane for which the orbit of 0 under iteration of the quadratic map

zn+1=zn2+cz_{n + 1} = z_{n}^2 + c

remains bounded. That is, a complex number c is part of the Mandelbrot set if, when starting with z_0 = 0 and applying the iteration repeatedly, the absolute value of z_n remains bounded however large n gets.

c∈M→lim⁡n→∞sup ∣zn+1∣≤2.c\in M \rightarrow \lim_{n\rightarrow\infty} sup\,|z_{n + 1}| \leq 2.

In code, the iteration process for a complex number c is simply as follows

@jit(nopython=True)
def mandelbrot(c, b, itmap):
    z = c
    for n in arange(b):
        if z.real * z.real + z.imag * z.imag > 4:
            return itmap(n, z, b)
        z = z*z + c
    return 0.0
0.4s
Python
Python 3

where b known as the bailout radius. And itemap will assign a real value to be used later for colouring the whole space of points. Actually, for itemap we will use the following normalized iteration count function:

normIterations(n,zn,b)=n−log⁡(log(∣zn∣))−log⁡(log⁡(b))log⁡(2)normIterations(n, z_{n}, b) = n - \frac{\log(log(|z_{n}|)) - \log(\log(b))}{\log(2)}

and in code as follows, (if someone knows the original reference for this normalization please let me know).

@jit(nopython=True)
def normIterations(n, zn, b):
    return n - (log(log(abs(zn))) - log(log(b)))/log(2.0)
0.3s
Python
Python 3

Then, a function that works for a whole section of the complex plane will be

@jit(nopython=True)
def mandelbrot_set(xi, xf, yi, yf, npoints, b):
    xycolor = np.zeros((npoints*npoints, 3))
    s = 0
    for x in np.linspace(xi, xf, npoints):
        for y in np.linspace(yi, yf, npoints):
            color = mandelbrot(x + 1j*y, b, normIterations)
            xycolor[s,:] = [x, y, color]
            s += 1
    return xycolor
0.2s
Python
Python 3

Disclaimer: Although this code is pretty fast, probably there are clever ways to do it. Nevertheless, I think this versions are easy to follow and understand.

Testing... datashader magic

xycm = mandelbrot_set(-2.0,0.5,-1.25,1.25,1000,1024)
3.0s
Python
Python 3

Please note that the previous array contains one million points.

def get_img(xyc, sx=400, sy = 400, pick_color = "fire", fhow = "eq_hist"):
    df = pd.DataFrame(xyc, columns =["x","y","color"]) # because ds takes DataFrames
    cvs = ds.Canvas(plot_width = sx, plot_height = sy)
    agg = cvs.points(df, 'x', 'y', ds.mean('color')) # ds.mean, ds.min, ds.max
    return tf.shade(agg, cmap=cm(palette[pick_color], 0.0), how= fhow)  # eq_hist, linear, cbrt, log
0.3s
Python
Python 3

We can set different kinds of colors from colorcet as well as scales for colouring.

escalas = ["linear", "cbrt", "log", "eq_hist"]
0.2s
Python
Python 3
tf.Images(*[get_img(xycm, fhow = scale) for scale in escalas])
2.5s
Python
Python 3

Let's do some more with different colours and scales.

xycm1 = mandelbrot_set(-0.74877,-0.74872,0.06505,0.06510,1000,2500)
xycm2 = mandelbrot_set(0.32642717997233067, 0.3265289052327473,
                       -0.05451000956418885, -0.05443371561887635,1000,3000)
xycm3 = mandelbrot_set(-1.9963806954442953, -1.996380695443582, 
                       2.628704923646517e-7, 2.62871027270105e-7,1000,1000)
xycm4 = mandelbrot_set(0.3476108223238668926295, 0.3476108223245338122665, 
                       0.0846794087369283253550, 0.0846794087374285150830,1000,10000)
24.4s
Python
Python 3
coordxyc = [xycm1, xycm2, xycm3, xycm4]
colores = ['bmw', 'bmy','fire', 'gray', 'kbc']
escalas = ["linear", "eq_hist"]
0.2s
Python
Python 3
tf.Images(*[get_img(xycv, pick_color = c, fhow = scale) for xycv in coordxyc 
            for scale in escalas for c in colores]).cols(5)
9.1s
Python
Python 3

Going deeper

# this looks ugly (a lot of numbers) but the pics are nice. 
xycm5 = mandelbrot_set(-0.650790400000000001192,-0.648844800000000001192, 
                       0.44539837880859369792, 0.44685757880859369792,1000,10_000)
xycm6 = mandelbrot_set(-0.9548899408372031, -0.9548896813770819,
                       0.2525416487455764, 0.2525418433406673, 1000,10_000)
xycm7 = mandelbrot_set(0.254828857465066226270, 0.254828889245416226270,
                       -0.000605561881950000235, -0.000605538046687500235,1000,50_000)
xycm8 = mandelbrot_set(-0.882297664710767940063, -0.882297662380940440063,
                       0.235365461981556923486,0.235365463728927548486, 1000, 10_000)
xycm9 = mandelbrot_set(-0.6534376561891502063520, -0.6534376520406489856480,
                       0.3635691455538367401120, 0.3635691486652126556400,1000,10_000)
xycm10 = mandelbrot_set(0.25740289813988496306, 0.25740289814296891154, 
                        0.49283797442018109421, 0.49283797442249278541,1000,50_000)
136.8s
Python
Python 3
coordxycDeeper = [xycm5, xycm6, xycm7, xycm8, xycm9, xycm10]
colores = ['bmw', 'fire', 'kbc']
0.2s
Python
Python 3
tf.Images(*[get_img(xycv, sx = 800, sy = 600, 
                    pick_color = c, fhow = "eq_hist") 
            for xycv in coordxycDeeper 
            for c in colores]).cols(3)
11.4s
Python
Python 3

everything works beautiful ! Now, let's see the julia set.

Julia set

The Julia Set associated to the function

fc(z)=z2+cf_c (z) = z^2 + c

for a constant c the complex plane is the set

J(fc)\mathcal{J}(f_c)

of points z in the complex plane, such that

∣fcn(z)∣=∣fc(fc(…fc(z)))∣≤2  ∀n|f_{c}^{n}(z)| = |f_{c} (f_{c}(\dots f_{c}(z)))| \leq 2\; \forall_{n}

The code for this condition is

@jit(nopython=True)
def julia_map(x, y, c, b, itmap):
    z = (x + 1j*y)**2 + c
    for k in arange(b):
        if z.real * z.real + z.imag * z.imag  > 4:
            return itmap(k, z, b)
        z = z*z + c
    return 0.0
0.3s
Python
Python 3

and for a region in the complex plane

@jit(nopython=True)
def create_julia_dots(xi, xf, yi, yf, c, b, npoints):
    xycolor = np.zeros((npoints**2, 3))
    s = 0
    for x in np.linspace(xi, xf, npoints):
        for y in np.linspace(yi, yf, npoints):
            color = julia_map(x, y, c, b, normIterations)
            xycolor[s,:] = [x, y, color]
            s +=1
    return xycolor
0.2s
Python
Python 3

And the plots are 

xyc_julia = create_julia_dots(-1.7, 1.7, -1.7, 1.7,
                              -0.835 - 0.2321*1j, 500, 2000)
3.0s
Python
Python 3
colores = ['bmw', 'bmy','fire', 'gray', 'kbc']
tf.Images(*[get_img(xyc_julia, pick_color = c, fhow ="log") 
            for c in colores])
2.4s
Python
Python 3

More examples... 

cvalues = [0.274 - 0.008 * 1j, 0.285 + 0.01*1j,
           -0.70176 - 0.3842*1j, -0.8 + 0.156*1j,
          -0.7269 + 0.1889*1j, -0.1 + 0.65*1j,
           -0.382 + 0.618*1j, -0.449 + 0.571*1j]
0.2s
Python
Python 3
def get_imgJulia(c, escala, color):
    xyc = create_julia_dots(-1.5, 1.5, -1.5, 1.5, c, 500, 1000)
    df = pd.DataFrame(xyc, columns =["x","y","c"])
    cvs = ds.Canvas(plot_width=400, plot_height=400)
    agg = cvs.points(df, 'x', 'y', ds.mean('c'))
    return tf.shade(agg, cmap=cm(palette[color], 0.0), how=escala)
0.3s
Python
Python 3
colores = ['bmw', 'bmy','fire', 'gray', 'kbc', 'bmw', 
           'bmy','fire', 'gray', 'kbc']
tf.Images(*[get_imgJulia(cvalue, "linear", colores[indx]) 
            for indx, cvalue in enumerate(cvalues)]).cols(4)
7.6s
Python
Python 3

One last colormap

tf.Images(*[get_imgJulia(cvalue, "linear", "CET_L17") 
            for indx, cvalue in enumerate(cvalues)]).cols(4)
7.0s
Python
Python 3

For more colors go to colorcet or simply load matplotlib colormaps.

That was fun !

Find me on twitter as @LazarusAlon