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
CPU times: user 360 ms, sys: 4 ms, total: 364 ms
Wall time: 499 ms

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


CPU times: user 384 ms, sys: 0 ns, total: 384 ms
Wall time: 1.1 s
CPU times: user 4.55 s, sys: 4 ms, total: 4.55 s
Wall time: 8.94 s
CPU times: user 356 ms, sys: 32 ms, total: 388 ms
Wall time: 598 ms

In [58]:
j = JuliaSetPlot(-2)
%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: 820 ms
CPU times: user 736 ms, sys: 0 ns, total: 736 ms
Wall time: 830 ms
CPU times: user 368 ms, sys: 40 ms, total: 408 ms
Wall time: 607 ms

In [59]:
j = JuliaSetPlot(-3)
%time j.set_spacing(0.006)
%time j.generate()
%time j.interact()


CPU times: user 372 ms, sys: 0 ns, total: 372 ms
Wall time: 388 ms
CPU times: user 596 ms, sys: 0 ns, total: 596 ms
Wall time: 824 ms
CPU times: user 352 ms, sys: 36 ms, total: 388 ms
Wall time: 595 ms

In [60]:
j = JuliaSetPlot(-4)
%time j.set_spacing(0.006)
%time j.generate()
%time j.interact()


CPU times: user 316 ms, sys: 0 ns, total: 316 ms
Wall time: 643 ms
CPU times: user 544 ms, sys: 0 ns, total: 544 ms
Wall time: 864 ms
CPU times: user 360 ms, sys: 36 ms, total: 396 ms
Wall time: 426 ms

In [61]:
j = JuliaSetPlot(-5)
%time j.set_spacing(0.006)
%time j.generate()
%time j.interact()


CPU times: user 184 ms, sys: 12 ms, total: 196 ms
Wall time: 291 ms
CPU times: user 492 ms, sys: 0 ns, total: 492 ms
Wall time: 650 ms
CPU times: user 360 ms, sys: 40 ms, total: 400 ms
Wall time: 783 ms

In [62]:
j = JuliaSetPlot(-6)
%time j.set_spacing(0.006)
%time j.generate()
%time j.interact()


CPU times: user 376 ms, sys: 0 ns, total: 376 ms
Wall time: 799 ms
CPU times: user 552 ms, sys: 0 ns, total: 552 ms
Wall time: 1.98 s
CPU times: user 368 ms, sys: 4 ms, total: 372 ms
Wall time: 917 ms

There is an apparent symmetry along the y axis.


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


CPU times: user 376 ms, sys: 4 ms, total: 380 ms
Wall time: 1.1 s
CPU times: user 4.5 s, sys: 8 ms, total: 4.51 s
Wall time: 11.7 s
CPU times: user 364 ms, sys: 24 ms, total: 388 ms
Wall time: 419 ms

In [67]:
j = JuliaSetPlot(-2)
%time j.set_spacing(0.006)
%time j.generate()
%time j.interact()


CPU times: user 364 ms, sys: 0 ns, total: 364 ms
Wall time: 387 ms
CPU times: user 744 ms, sys: 0 ns, total: 744 ms
Wall time: 814 ms
CPU times: user 368 ms, sys: 24 ms, total: 392 ms
Wall time: 534 ms

In [68]:
j = JuliaSetPlot(-3)
%time j.set_spacing(0.006)
%time j.generate()
%time j.interact()


CPU times: user 364 ms, sys: 4 ms, total: 368 ms
Wall time: 375 ms
CPU times: user 596 ms, sys: 0 ns, total: 596 ms
Wall time: 606 ms
CPU times: user 368 ms, sys: 32 ms, total: 400 ms
Wall time: 425 ms

In [69]:
j = JuliaSetPlot(-4)
%time j.set_spacing(0.006)
%time j.generate()
%time j.interact()


CPU times: user 356 ms, sys: 0 ns, total: 356 ms
Wall time: 948 ms
CPU times: user 564 ms, sys: 0 ns, total: 564 ms
Wall time: 2.11 s
CPU times: user 372 ms, sys: 40 ms, total: 412 ms
Wall time: 643 ms

Negative and positive c values return the same graphs perhaps inverted along the y axis.


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


CPU times: user 384 ms, sys: 0 ns, total: 384 ms
Wall time: 880 ms
CPU times: user 1.04 s, sys: 8 ms, total: 1.04 s
Wall time: 1.84 s
CPU times: user 356 ms, sys: 52 ms, total: 408 ms
Wall time: 1.15 s

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


CPU times: user 384 ms, sys: 0 ns, total: 384 ms
Wall time: 1.28 s
CPU times: user 692 ms, sys: 0 ns, total: 692 ms
Wall time: 1.21 s
CPU times: user 364 ms, sys: 44 ms, total: 408 ms
Wall time: 738 ms

In [ ]: