In [ ]:

    

%load_ext cythonmagic


    

In [ ]:

    

import numpy as np
from time import time

def kernel(z, c, lim, cutoff=1e6):
    ''' Computes the number, `n`, of iterations necessary such that 
    |z_n| > `lim`, where `z_n = z_{n-1}**2 + c`.
    '''
    count = 0
    while abs(z) < lim and count < cutoff:
        z = z * z + c
        count += 1
    return count


def compute_julia(c, N, bound=2, lim=1000., kernel=kernel):
    ''' Pure Python calculation of the Julia set for a given `c`.  No NumPy
    array operations are used.
    '''
    julia = np.empty((N, N), dtype=np.uint32)
    grid_x = np.linspace(-bound, bound, N)
    grid_y = grid_x * 1j
    c = complex(c)
    t0 = time()
    for i, x in enumerate(grid_x):
        for j, y in enumerate(grid_y):
            julia[i,j] = kernel(x+y, c, lim)
    return julia, time() - t0


    

In [ ]:

    

import pylab as pl

def plot_julia(kwargs, compute_julia):
    ''' Given parameters dict in `kwargs` and a function to compute the Julia
    set (`compute_julia`), plots the resulting Julia set with appropriately
    labeled axes.
    '''
    kwargs = kwargs.copy()

    def _plotter(kwargs):
        bound = kwargs['bound']
        julia, time = compute_julia(**kwargs)
        print "execution time (s):", time
        julia = np.log(julia)
        pl.imshow(julia, 
                  interpolation='nearest',
                  extent=(-bound, bound)*2)
        pl.colorbar()
        title = r"Julia set for $C={0.real:5.3f}+{0.imag:5.3f}i$ $[{1}\times{1}]$"
        pl.title(title.format(kwargs['c'], kwargs['N']))
        pl.xlabel("$Re(z)$")
        pl.ylabel("$Im(z)$")

    pl.figure(figsize=(14, 12))

    cvals = [0.285+0.01j, -0.1+0.651j, -0.4+0.6j, -0.8+0.156j]
    subplots = ['221',    '222',       '223',     '224'      ]

    for c, sp in zip(cvals, subplots):
        kwargs.update(c=c)
        pl.subplot(sp)
        _plotter(kwargs)

    pl.show()


    

In [ ]:

    

kwargs = dict(N=100, bound=1.5)
plot_julia(kwargs, compute_julia)


    

In [ ]:

    

%%cython

from time import time

import numpy as np
cimport cython
cimport numpy as cnp

ctypedef double complex cpx_t
ctypedef double         real_t

cdef inline real_t cabs_sq(cpx_t z) nogil:
    ''' Helper inline function, computes the square of the abs. value of the
    complex number `z`.
    '''
    return z.real * z.real + z.imag * z.imag
    
cpdef unsigned int kernel(cpx_t z, 
                                 cpx_t c,
                                 real_t lim,
                                 real_t cutoff=1e6) nogil:
    ''' Cython implementation of the kernel computation.

    This is implemented so that no C-API calls are made inside the function
    body.  Even still, there is some overhead as compared with a pure C
    implementation.
    '''
    cdef unsigned int count = 0
    cdef real_t lim_sq = lim * lim
    while cabs_sq(z) < lim_sq and count < cutoff:
        z = z * z + c
        count += 1
    return count

@cython.boundscheck(False)
@cython.wraparound(False)
def compute_julia_opt(cpx_t c,
                       unsigned int N,
                       real_t bound=1.5,
                       real_t lim=1000.):
    '''
    Cython `compute_julia()` implementation with Numpy array buffer
    declarations and appropriate compiler directives.  The body of this
    function is nearly identical to the `compute_julia_no_opt()` function.

    '''

    cdef cnp.ndarray[cnp.uint32_t, ndim=2, mode='c'] julia 
    cdef cnp.ndarray[real_t, ndim=1, mode='c'] grid
    cdef unsigned int i, j
    cdef real_t x, y

    julia = np.empty((N, N), dtype=np.uint32)
    grid = np.linspace(-bound, bound, N)
    t0 = time()
    for i in range(N):
        x = grid[i]
        for j in range(N):
            y = grid[j]
            julia[i,j] = kernel(x+y*1j, c, lim)
    return julia, time() - t0


    

In [ ]:

    

kwargs = dict(N=1000, bound=1.5)
plot_julia(kwargs, compute_julia_opt)


    

In [ ]: