Homework 5: JuliaSet optimization with numpy


In [37]:
import numpy as np

class JuliaSet(object):

    def __init__(self, c, n=100):
        self.c = c
        self.n = n
        self._d = 0.001
        self._complexplane= np.array([])
        self.set= np.array([])

    def juliamap(self, z):
        return z**2 + self.c

    def iterate(self, z):
        m = 0
        while True:
            z = self.juliamap(z)
            m +=1
            if abs(z)>2:
                return m
            elif m>=self.n:
                return 0
    """
    def drange(self, start, stop, step):
        r = start
        while r <= stop:
            yield r
            r += step
    """
    
    # Generate evenly spaced values over x and y planes
    def get_complexplane(self):
        r=np.arange(-2, 2, self._d)
        x, y = np.meshgrid(r,r)
        self._complexplane = x + y*1j
        return self._complexplane
        
    
    def set_spacing(self, d):
        self._d = d
        self.get_complexplane()
        
    def generate(self):
        for x in self._complexplane:
            self.set (self.iterate(x), self.iterate(x))
        return self.set

Julia Set Plotting Extension

Load module for a JuliaSet that conforms to the specified interface.

It is wise to run the test suite in test_juliaset.py with nosetests prior to attempting to plot here.


In [38]:
from juliaset import JuliaSet
# Math libraries
import numpy as np
from math import sqrt

# Matplotlib plotting libraries
import matplotlib.pyplot as plt
%matplotlib inline

# Bokeh plotting libraries
import bokeh.plotting as blt
blt.output_notebook()

class JuliaSetPlot(JuliaSet):
    """Extend JuliaSet to add plotting functionality"""
    
    def __init__(self, *args, **kwargs):
        # Invoke constructor for JuliaSet first, unaltered
        JuliaSet.__init__(self, *args, **kwargs)
        # Add another attribute: a rendered image array
        self.img = np.array([])
    
    def get_dim(self):
        """Return linear number of points in axis"""
        return int(4.0 / self._d)
    
    def render(self):
        """Render image as square array of ints"""
        if not self.set: self.generate()
        # Convert inefficient list to efficient numpy array
        self.img = np.array(self.set)
        # Reshape array into a 2d complex plane
        dim = int(sqrt(len(self.img)))
        self.img = np.reshape(self.img, (dim,dim)).T
        
    def show(self):
        """Use matplotlib to plot image as an efficient mesh"""
        if not self.img.size: self.render()
        plt.figure(1, figsize=(12,9))
        xy = np.linspace(-2,2,self.get_dim())
        plt.pcolormesh(xy, xy, self.img, cmap=plt.cm.hot)
        plt.colorbar()
        plt.show()
        
    def interact(self):
        """Use bokeh to plot an interactive image"""
        from matplotlib.colors import rgb2hex
        if not self.img.size: self.render()
        # Mimic matplotlib "hot" color palette
        colormap = plt.cm.get_cmap("hot")
        bokehpalette = [rgb2hex(m) for m in colormap(np.arange(colormap.N))]
        # Create bokeh figure
        f = blt.figure(x_range=(-2,2), y_range=(-2,2), plot_width=600, plot_height=600)
        f.image(image=[self.img], x=[-2], y=[-2], dw=[4], dh=[4], palette=bokehpalette, dilate=True)
        blt.show(f)


BokehJS successfully loaded.

Optimization of Juliaset using Numpy

The first change to the previous Juliaset algorithm was to use the numpy array type for self._complexplane and self.set arrays.

The second change was to use np.arange to create evenly spaced vectors for the complexeplane array instead of using the for loop to iterate over values between -2 and 2 for x and y.

The third change was using np.meshgrid to create the coordinate matrices from the vectors for the complexplane.

Overall the Juliasets3 algorithm performs at around 3 times faster Wall-clock time for the both the generate() and set_spacing() functions than the algorithm of juliasets2.

Benchmark info: To run the function: j = JuliaSetPlot( -0.835 - 0.2321j) %time j.set_spacing(0.006) %time j.generate() %time j.interact()

Juliasets2 results: CPU times: user 388 ms, sys: 16 ms, total: 404 ms Wall time: 1.06 s CPU times: user 2 s, sys: 36 ms, total: 2.04 s Wall time: 6.63 s

Juliasets3 results: CPU times: user 372 ms, sys: 8 ms, total: 380 ms Wall time: 383 ms CPU times: user 1.89 s, sys: 0 ns, total: 1.89 s Wall time: 2.02 s


In [40]:
j = JuliaSetPlot( -0.835 - 0.2321j)
%time j.set_spacing(0.006)
%time j.generate()
%time j.interact()


CPU times: user 372 ms, sys: 8 ms, total: 380 ms
Wall time: 383 ms
CPU times: user 1.89 s, sys: 0 ns, total: 1.89 s
Wall time: 2.02 s
CPU times: user 364 ms, sys: 24 ms, total: 388 ms
Wall time: 507 ms

In [39]:
j = JuliaSetPlot(1)
%time j.set_spacing(0.006)
%time j.generate()
%time j.interact()


CPU times: user 368 ms, sys: 4 ms, total: 372 ms
Wall time: 399 ms
CPU times: user 828 ms, sys: 0 ns, total: 828 ms
Wall time: 1e+03 ms
CPU times: user 372 ms, sys: 56 ms, total: 428 ms
Wall time: 2.03 s

Testing of juliasets settings

Setting a very high c value returns an empty plot as seen with the value 100.


In [42]:
j = JuliaSetPlot(100)
%time j.set_spacing(0.006)
%time j.generate()
%time j.interact()


CPU times: user 388 ms, sys: 0 ns, total: 388 ms
Wall time: 395 ms
CPU times: user 548 ms, sys: 0 ns, total: 548 ms
Wall time: 706 ms