In [1]:
%pylab inline
%load_ext cythonmagic
from scipy.integrate import simps


Populating the interactive namespace from numpy and matplotlib

In [2]:
import pyximport; pyximport.install(setup_args={"include_dirs":numpy.get_include()})
import paramless as pmain
import paramless_cython as pm

Using Cython to improve the runtime

This notebook will show a simple example of improving the runtime of an evolution run using Cython functions. To do this we'll build on the Seasonal Flowering example.

The functions used there are given below:


In [3]:
class FloweringFitness(pm.FitnessFunction):
    def __init__(self, domain, alpha, beta):
        self.domain = domain
        self.alpha = alpha
        self.beta = beta
        
    def get(self, resident, mutant = None):
        pass

In [4]:
def p_x(self, a, x, k):
    #takes a, and functions x and k
    return math.e**(-x[a]/k(a))

FloweringFitness.p_x = p_x

In [5]:
def k(self, a):
    return 100*(2 + math.sin(2*math.pi*((a-1.0)**2.0 - 0.25)))

FloweringFitness.k_function = k

In [6]:
def c(self, e, alpha):
    return 2.0/(1 + math.e**(-alpha*e))

FloweringFitness.c = c

In [7]:
def birth(self, x_mutant, x_resident, domain, alpha, beta, k_function):
    x_resident_dict = dict(zip(domain, x_resident))
    x_mutant_dict = dict(zip(domain, x_mutant))
    integrand = [x_mutant_dict[a]*self.c((x_mutant_dict[a]-x_resident_dict[a])/(x_mutant_dict[a]**beta), alpha)*self.p_x(a, x_resident_dict, k_function)
                 if x_mutant_dict[a] != 0.0
                 else 0.0
                 for a in domain]
    return simps(integrand, domain)

FloweringFitness.birth = birth

In [8]:
def flowering_fitness(self, resident, mutant = None):
    fitness_resident = self.birth(resident, resident, self.domain, self.alpha, self.beta, self.k_function)
    if not mutant is None:    
        fitness_mutant = self.birth(mutant, resident, self.domain, self.alpha, self.beta, self.k_function)
        return fitness_resident, fitness_mutant
    else:
        return fitness_resident

FloweringFitness.get = flowering_fitness

These can all be redefined using Cython definitions and types:


In [9]:
%%cython
import numpy as np
cimport numpy as np
import math
from scipy.integrate import simps
from libc.math cimport exp, sin

cdef class cyFloweringFitness:
    #Class members need to be declared before assignment
    cdef double [:] domain
    cdef double alpha, beta
    def __init__(self, np.ndarray[double, ndim=1] domain, double alpha, double beta):
        self.domain = domain
        self.alpha = alpha
        self.beta = beta
        
    cdef double p_x(self, double a, x):
        return exp(-x[a]/self.k(a))

    cdef double k(self, double a):
        return 100*(2 + sin(2*math.pi*((a-1.0)**2.0 - 0.25)))

    cdef double c(self, e):
        return 2.0/(1 + exp(-self.alpha*e))

    cdef double birth(self, double [:] x_mutant, double [:] x_resident):
        x_resident_dict = dict(zip(np.asarray(self.domain), np.asarray(x_resident)))
        x_mutant_dict = dict(zip(np.asarray(self.domain), np.asarray(x_mutant)))
        integrand = [x_mutant_dict[a]*self.c((x_mutant_dict[a]-x_resident_dict[a])/(x_mutant_dict[a]**self.beta))*self.p_x(a, x_resident_dict) for a in self.domain]
        return simps(integrand, self.domain)
        
        
    def get(self, double [:] resident, double [:] mutant = None):
        fitness_resident = self.birth(resident, resident)
        if not mutant is None:    
            fitness_mutant = self.birth(mutant, resident)
            return fitness_resident, fitness_mutant
        else:
            return fitness_resident

An alternative method for utilising Cython is through the "Pure Python Mode", described at http://docs.cython.org/src/tutorial/pure.html

Set up timing


In [10]:
import timeit
times = []
names = []

In [11]:
domain = np.arange(0.001, 1.0, 0.01)

In [12]:
initial_surface = np.ones_like(domain)

In [13]:
number_of_generations = 5000

In [16]:
%%capture
seed = 777

fitness_function = FloweringFitness(domain, 1.0, 0.9)
mutator = pm.GaussianDistributionMutator(0.01, domain, 0.05)
evolver = pm.StandardEvolver(fitness_function, mutator)
pm.setup(seed)

ans, series = pmain.evolve(initial_surface, evolver,
                     iterations=number_of_generations,
                     return_time_series=True, seed=seed)

times.append(min(timeit.Timer('pmain.evolve(initial_surface=initial_surface, evolver=evolver, iterations=number_of_generations, return_time_series=True, seed=seed)', 
            'from __main__ import initial_surface, evolver, number_of_generations, seed;import pyximport; pyximport.install();import paramless_cython as pm;import paramless as pmain').repeat(repeat=3, number=3)))
names.append('Original')

In [18]:
%%capture
cy_fitness = cyFloweringFitness(domain, 1.0, 0.9)
mutator = pm.GaussianDistributionMutator(0.01, domain, 0.05)
evolver = pm.StandardEvolver(cy_fitness, mutator)
pm.setup(seed)

cyAns, series = pmain.evolve(initial_surface, evolver,
                     iterations=number_of_generations,
                     return_time_series=True, seed=seed)

times.append(min(timeit.Timer('pmain.evolve(initial_surface=initial_surface, evolver=evolver, iterations=number_of_generations, return_time_series=True, seed=seed)', 
            'from __main__ import initial_surface, evolver, number_of_generations, seed;import paramless as pmain').repeat(repeat=3, number=3)))
names.append('Cython')

Ensure that the code is equivalent:


In [19]:
plot(domain, ans)
plot(domain,cyAns)


Out[19]:
[<matplotlib.lines.Line2D at 0x7f1cbfe3bc90>]

In [20]:
#Adapted from http://nbviewer.ipython.org/github/rasbt/python_reference/blob/master/tutorials/running_cython.ipynb
from matplotlib import pyplot as plt
import numpy as np

fig = plt.figure(figsize=(10,8))

# plot bars
y_pos = np.arange(len(times)) + .5
plt.yticks(y_pos, names, fontsize=14)
bars = plt.barh(y_pos, times,
         align='center', alpha=0.4, color='g')

# annotation and labels

for b,d in zip(bars, times):
    plt.text(max(times)+1.5, b.get_y() + b.get_height()/2.5, 
             '{:.3}'.format(d),
             ha='center', va='bottom', fontsize=12)

plt.ylim([-1,len(times)+0.5])
plt.grid()
plt.ylim([0,len(times)])

plt.show()



In [ ]: