Mandelbrot¶
In this tutorial, we will take a look at complex numbers and one of their applications. We will try to figure out on how to compute a simple mandelbrot. Lets`s start with a simple example on how python handles complex numbers.
Let's start with a basic Python Mandelbrot set. We use numpy arrays for the image and display it using pylab imshow.
from PIL import Image
import numpy as np
from numba import autojit
from timeit import default_timer as timer
from skimage import data, exposure
The following few lines just briefly demonstrate on how to use complex numbers in python.
x = 1
y = 1
z = complex(x,y)
print z
print z.real
print z.imag
The complex plane
from matplotlib import pylab as plt
plt.figure(figsize = (6,4))
plt.grid()
plt.plot(z.real, z.imag, "o")
plt.arrow(0,0, z.real-0.1, z.imag-0.1,
head_width=0.05,
head_length=0.1,
fc='b',
ec='b')
plt.xlim([-2.1,+2.1])
plt.ylim([-1.1,+1.1])
plt.xlabel("real")
plt.ylabel("imag")
plt.show()
Mandelbrot computation¶
The following function computes a mandelbrot set for a given point in the complex plane $z = x + yj$. It returns the number of iterations before the computation "escapes".
def Mandel(x,y, itter = 1, thresh = 4):
c = complex(x,y)
z = 0.0j
for i in range(itter):
z = z**2+c
if (z.real**2. + z.imag**2) >=thresh:
return i
return itter
The following function performs a mandelbrot set calculation for a given set of points within a 2D array. It also scales between pixel space and complex space. It returns a 2D numpy array. With this function we can produce a single mandlebrot set by giving it the image dimensions, a cutoff value and the number of maximum iterations.
def draw_Mandel(min_x, max_x, min_y, max_y, image, iters, thresh):
height = image.shape[0] ### frame dimension y
width = image.shape[1] ### frame dimension x
pixel_size_x = (max_x - min_x) / width ### pixel dimension x
pixel_size_y = (max_y - min_y) / height ### pixel dimesnion y
for x in range(width):
real = min_x + x * pixel_size_x
for y in range(height):
imag = min_y + y * pixel_size_y
color = Mandel(real, imag, iters, thresh)
image[y, x] = color
return image
To make the mandelbrot computation a little bit more interesting, we can generate a number of images while zooming in after the 20th iteration. We will also use the timeit module to determin the time it takes to compute the individual sets.
def print_example(Z = 1, resx = 1920, resy = 1080):
ZOOM = 1/float(Z)
ratio = float(resx)/float(resy)
X, Y = -0.74814501, -0.071673006
if Z > 20:
Z = Z/10
XMIN, XMAX = X-ZOOM*ratio, X+ZOOM*ratio
YMIN, YMAX = Y-ZOOM, Y+ZOOM
else:
XMIN, XMAX = X-1.*ratio, X+1.*ratio
YMIN, YMAX = Y-1., Y+1.
plt.figure(figsize = (16,8))
image = np.zeros((resy,resx))
image = draw_Mandel(XMIN, XMAX, YMIN, YMAX, image, Z,6)
#np.save("mandelbrot",image)
p2, p98 = np.percentile(image, (2, 98))
im = plt.imshow(exposure.rescale_intensity(image, in_range=(p2, p98)), cmap ="PRGn", extent=[XMIN, XMAX,YMIN, YMAX])
plt.colorbar(im,fraction=0.046, pad=0.005)
plt.ylabel("Im")
plt.xlabel("Re")
plt.grid()
plt.text(XMIN,YMIN, "x"+str(Z), fontsize = 25, color = "w")
plt.show()
timesCPU = []
for i in [1,2,3,4,10,20]:
start = timer()
print_example(i)
dt = timer() - start
print "time to compute "+str(i)+" iterations: "+str(dt)+"s"
timesCPU.append(dt)
Let's stop after 20 iterations and lets plot the times it takes to compute each set of iterations. This appears to take very long for long eterations. The plot below shows the time per set.
bar_width = 0.4
plt.figure(figsize = (12,4))
plt.bar(np.arange(len(timesCPU)), timesCPU,bar_width,label = "CPU")
plt.xticks(np.arange(len(timesCPU)), [1,2,3,4,10,20])
plt.legend()
plt.xlabel("iteration")
plt.ylabel("time [s]")
plt.show()
Compute on the GPU¶
To optimise and speed up the computation, we can use numba's autojit and hope that it utelises the GPU.
@autojit
def Mandel(x,y, itter = 1, thresh = 4):
c = complex(x,y)
z = 0.0j
for i in range(itter):
z = z**2+c
if (z.real**2. + z.imag**2) >=thresh:
return i
return itter
def print_example(Z = 1, resx = 1920, resy = 1080):
ZOOM = 1/float(Z)
ratio = float(resx)/float(resy)
X, Y = -0.74814501, -0.071673006
if Z > 20:
Z = Z/10
XMIN, XMAX = X-ZOOM*ratio, X+ZOOM*ratio
YMIN, YMAX = Y-ZOOM, Y+ZOOM
else:
XMIN, XMAX = X-1.*ratio, X+1.*ratio
YMIN, YMAX = Y-1., Y+1.
plt.figure(figsize = (16,8))
image = np.zeros((resy,resx))
image = draw_Mandel(XMIN, XMAX, YMIN, YMAX, image, Z,6)
#np.save("mandelbrot",image)
p2, p98 = np.percentile(image, (2, 98))
im = plt.imshow(exposure.rescale_intensity(image, in_range=(p2, p98)), cmap ="PRGn", extent=[XMIN, XMAX,YMIN, YMAX])
plt.colorbar(im,fraction=0.046, pad=0.005)
plt.ylabel("Im")
plt.xlabel("Re")
plt.grid()
plt.text(XMIN,YMIN, "x"+str(Z), fontsize = 25, color = "w")
plt.show()
timesGPU = []
for i in [1,2,3,4,10,20,100,1000,10000,100000,1000000,1000000000000]:
start = timer()
print_example(i)
dt = timer() - start
print "time to compute "+str(i)+" iterations: "+str(dt)+"s"
timesGPU.append(dt)
