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]
conda install datashader
conda install numba
conda install colorcet

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

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.

cMlimnsupzn+12.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

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

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)=nlog(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).

def normIterations(n, zn, b):
    return n - (log(log(abs(zn))) - log(log(b)))/log(2.0)

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

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

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)

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

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

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

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)
coordxyc = [xycm1, xycm2, xycm3, xycm4]
colores = ['bmw', 'bmy','fire', 'gray', 'kbc']
escalas = ["linear", "eq_hist"]
tf.Images(*[get_img(xycv, pick_color = c, fhow = scale) for xycv in coordxyc 
            for scale in escalas for c in colores]).cols(5)

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)
coordxycDeeper = [xycm5, xycm6, xycm7, xycm8, xycm9, xycm10]
colores = ['bmw', 'fire', 'kbc']
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)

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


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

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

and for a region in the complex plane

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

And the plots are

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

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

One last colormap

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

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

That was fun !

Find me on twitter as @LazarusAlon