# A small intro on the rationale

So I’m working on a Symbolic Regression Machine written in C/C++ called Shine, which is intended to be a JIT for Genetic Programming libraries (like Pyevolve for instance). The main rationale behind Shine is that we have today a lot of research on speeding Genetic Programming using GPUs (the GPU fever !) or any other special hardware, etc, however we don’t have many papers talking about optimizing GP using the state of art compilers optimizations like we have on clang, gcc, etc.

The “hot spot” or the component that consumes a lot of CPU resources today on Genetic Programming is the evaluation of each individual in order to calculate the fitness of the program tree. This evaluation is often executed on each set of parameters of the “training” set. Suppose you want to make a symbolic regression of a single expression like the Pythagoras Theorem and you have a linear space of parameters from 1.0 to 1000.0 with a step of 0.1 you have 10.000 evaluations for each individual (program tree) of your population !

What Shine does is described on the image below:

It takes the individual of the Genetic Programming engine and then converts it to LLVM Intermediate Representation (LLVM assembly language), after that it runs the transformation passes of the LLVM (here is where the true power of modern compilers enter on the GP context) and then the LLVM JIT converts the optimized LLVM IR into native code for the specified target (X86, PowerPC, etc).

You can see below the Shine architecture:

This architecture brings a lot of flexibility for Genetic Programming, you can for instance write functions that could be used later on your individuals on any language supported by the LLVM, what matters to Shine is the LLVM IR, you can use any language that LLVM supports and then use the IR generated by LLVM, you can mix code from C, C++, Ada, Fortran, D, etc and use your functions as non-terminal nodes of your Genetic Programming trees.

Shine is still on its earlier development, it looks a simple idea but I still have a lot of problems to solve, things like how to JIT the evaluation process itself instead of doing calls from Python using ctypes bindings of the JITed trees.

# Doing Genetic Programming on the Python AST itself

During the development of Shine, an idea happened to me, that I could use a restricted Python Abstract Syntax Tree (AST) as the representation of individuals on a Genetic Programming engine, the main advantage of this is the flexibility and the possibility to reuse a lot of things. Of course that a shared library written in C/C++ would be useful for a lot of Genetic Programming engines that doesn’t uses Python, but since my spare time to work on this is becoming more and more rare I started to rethink the approach and use Python and the LLVM bindings for LLVM (LLVMPY) and I just discovered that is pretty easy to JIT a restricted set of the Python AST to native code using LLVM, and this is what this post is going to show.

# JIT’ing a restricted Python AST

The most amazing part of LLVM is obviously the amount of transformation passes, the JIT and of course the ability to use the entire framework through a simple API (ok, not so simple sometimes). To simplify this example, I’m going to use an arbitrary restricted AST set of the Python AST that supports only subtraction (-), addition (+), multiplication (*) and division (/).

To understand the Python AST, you can use the Python parser that converts source into AST:

 12345678910 >>> import ast >>> astp = ast.parse("2*7") >>> ast.dump(astp) 'Module(     body=[Expr(         value=BinOp(             left=Num(n=2), op=Mult(), right=Num(n=7)         )     )] )'

What the parse created was an Abstract Syntax Tree containing the BinOp (Binary Operation) with the left operator as the number 2, the right operator as the number 7 and the operation itself as Multiplication(Mult), very easy to understand. What we are going to do to create the LLVM IR is to create a visitor that is going to visit each node of the tree. To do that, we can subclass the Python NodeVisitor class from the ast module. What the NodeVisitor does is to visit each node of the tree and then call the method ‘visit_OPERATOR’ if it exists, when the NodeVisitor is going to visit the node for the BinOp for example, it will call the method ‘visit_BinOp’ passing as parameter the BinOp node itself.

The structure of the class for for the JIT visitor will look like the code below:

 12345678910 # Import the ast and the llvm Python bindings import ast from llvm import * from llvm.core import * from llvm.ee import * import llvm.passes as lp class AstJit(ast.NodeVisitor):     def __init__(self):         pass

What we need to do now is to create an initialization method to keep the last state of the JIT visitor, this is needed because we are going to JIT the content of the Python AST into a function and the last instruction of the function needs to return what was the result of the last instruction visited by the JIT. We also need to receive a LLVM Module object in which our function will be created as well the closure type, for the sake of simplicity I’m not type any object, I’m just assuming that all numbers from the expression are integers, so the closure type will be the LLVM integer type.

 1234567891011121314151617181920212223242526272829 def __init__(self, module, parameters):         self.last_state = None         self.module = module         # Parameters that will be created on the IR function         self.parameters = parameters         self.closure_type = Type.int()         # An attribute to hold a link to the created function         # so we can use it to JIT later         self.func_obj = None         self._create_builder()     def _create_builder(self):         # How many parameters of integer type         params = [self.closure_type] * len(self.parameters)         # The prototype of the function, returning a integer         # and receiving the integer parameters         ty_func = Type.function(self.closure_type, params)         # Add the function to the module with the name 'func_ast_jit'         self.func_obj = self.module.add_function(ty_func, 'func_ast_jit')         # Create an argument in the function for each parameter specified         for index, pname in enumerate(self.parameters):             self.func_obj.args[index].name = pname         # Create a basic block and the builder         bb = self.func_obj.append_basic_block("entry")         self.builder = Builder.new(bb)

Now what we need to implement on our visitor is the ‘visit_OPERATOR’ methods for the BinOp and for the Numand Name operators. We will also implement the method to create the return instruction that will return the last state.

 12345678910111213141516171819202122232425262728293031323334353637383940 # A 'Name' is a node produced in the AST when you     # access a variable, like '2+x+y', 'x' and 'y' are     # the two names created on the AST for the expression.     def visit_Name(self, node):         # This variable is what function argument ?         index = self.parameters.index(node.id)         self.last_state = self.func_obj.args[index]         return self.last_state     # Here we create a LLVM IR integer constant using the     # Num node, on the expression '2+3' you'll have two     # Num nodes, the Num(n=2) and the Num(n=3).     def visit_Num(self, node):         self.last_state = Constant.int(self.closure_type, node.n)         return self.last_state       # The visitor for the binary operation     def visit_BinOp(self, node):         # Get the operation, left and right arguments         lhs = self.visit(node.left)         rhs = self.visit(node.right)         op = node.op         # Convert each operation (Sub, Add, Mult, Div) to their         # LLVM IR integer instruction equivalent         if isinstance(op, ast.Sub):             op = self.builder.sub(lhs, rhs, 'sub_t')         elif isinstance(op, ast.Add):             op = self.builder.add(lhs, rhs, 'add_t')         elif isinstance(op, ast.Mult):             op = self.builder.mul(lhs, rhs, 'mul_t')         elif isinstance(op, ast.Div):             op = self.builder.sdiv(lhs, rhs, 'sdiv_t')                 self.last_state = op         return self.last_state     # Build the return (ret) statement with the last state     def build_return(self):         self.builder.ret(self.last_state)

And that is it, our visitor is ready to convert a Python AST to a LLVM IR assembly language, to run it we’ll first create a LLVM module and an expression:

 12345 module = Module.new('ast_jit_module') # Note that I'm using two variables 'a' and 'b' expr = "(2+3*b+33*(10/2)+1+3/3+a)/2" node = ast.parse(expr) print ast.dump(node)

Will output:

Module(body=[Expr(value=BinOp(left=BinOp(left=BinOp(left=BinOp(


Now we can finally run our visitor on that generated AST the check the LLVM IR output:

 1234 visitor = AstJit(module, ['a', 'b']) visitor.visit(node) visitor.build_return() print module

Will output the LLVM IR:

; ModuleID = 'ast_jit_module'

define i32 @func_ast_jit(i32 %a, i32 %b) {
entry:
%mul_t = mul i32 3, %b
%sdiv_t = sdiv i32 %add_t4, 2
ret i32 %sdiv_t
}


Now is when the real fun begins, we want to run LLVM optimization passes to optimize our code with an equivalent GCC -O2 optimization level, to do that we create a PassManagerBuilder and a PassManager, the PassManagerBuilder is the component that adds the passes to the PassManager, you can also manually add arbitrary transformations like dead code elimination, function inlining, etc:

 12345678910 pmb = lp.PassManagerBuilder.new() # Optimization level pmb.opt_level = 2 pm = lp.PassManager.new() pmb.populate(pm) # Run the passes into the module pm.run(module) print module

Will output:

; ModuleID = 'ast_jit_module'

define i32 @func_ast_jit(i32 %a, i32 %b) nounwind readnone {
entry:
%mul_t = mul i32 %b, 3
%sdiv_t = sdiv i32 %add_t4, 2
ret i32 %sdiv_t
}


And here we have the optimized LLVM IR of the Python AST expression. The next step is to JIT that IR into native code and then execute it with some parameters:

 123456 ee = ExecutionEngine.new(module)     arg_a = GenericValue.int(Type.int(), 100)     arg_b = GenericValue.int(Type.int(), 42)         retval = ee.run_function(visitor.func_obj, [arg_a, arg_b])     print "Return: %d" % retval.as_int()

Will output:

Return: 197


And that’s it, you have created a AST->LLVM IR converter, optimized the LLVM IR with the transformation passes and then converted it to native code using the LLVM execution engine. I hope you liked =)

# Darwin on the track

From The Economist article:

WHILE watching the finale of the Formula One grand-prix season on television last weekend, your correspondent could not help thinking how Darwinian motor racing has become. Each year, the FIA, the international motor sport’s governing body, sets new design rules in a bid to slow the cars down, so as to increase the amount of overtaking during a race—and thereby make the event more interesting to spectators and television viewers alike. The aim, of course, is to keep the admission and television fees rolling in. Over the course of a season, Formula One racing attracts a bigger audience around the world than any other sport.

# New SIGEVOlution Volume 5 Issue 1

From Pier Luca Lanzi site, the issue features:

• EC Testing of Embedded Systems by Peter M. Kruse, Joachim Wegener
and Stefan Wappler
• Competitions @ GECCO-2010 by Christian Gagné
• Events reports: NICSO-2010 by David Pelta
• Dissertation corner
• New issues of journals
• Calls & calendar

# Pyevolve 0.6rc1 released !

I’m proud to announce the Pyevolve 0.6rc1 ! This is the first release candidate, but it’s pretty stable for production use (as from this 0.6 version we are very closer to a stable codebase, thanks to community).

See the documentation site and the What’s New.

I would like to thank the people who directly or indirectly have contributed with this release: Boris Gorelik, Amit Saha, Jelle Feringa, Henrik Rudstrom, Matteo de Felice, SIGEVOlution Board, Mike Benoit, Ryan Campbell, Jonas A. Gustavsson, Soham Sadhu, Ernesto Costa, Ido Ben-Tsion, Frank Goodman, Vishal Belsare, Benjamin Bush; and a lot of people who gave us feedback of the experience with the use of Pyevolve on their applications.

If you see something wrong with this release candidate, please create a ticket reporting your problem at Pyevolve Trac, so we can provide fixes to the official release.

Happy coding !

– Christian S. Perone

# New issue of SIGEVOlution (Volume 4 Issue 3)

The new issue of SIGEVOlution (the newsletter of ACM Special Interest Group on Genetic and Evolutionary Computation) was released:

Issues in applying computational intelligence
By Arthur Kordon

JavaXCSF
By Patrick O. Stalph, Martin V. Butz

And a lot of information about new PhD theses, new journal issues and about events to come.

# Successful pyevolve multiprocessing speedup for Genetic Programming

As we know, Genetic Programming usually requires intensive processing power for the fitness functions and tree manipulations (in crossover operations), and this fact can be a huge problem when using a pure Python approach like Pyevolve. So, to overcome this situation, I’ve used the Python multiprocessing features to implement a parallel fitness evaluation approach in Pyevolve and I was surprised by the super linear speedup I got for a cpu bound fitness function used to do the symbolic regression of the Pythagoras theorem: $c = \sqrt{a^2 + b^2}$. I’ve used the same seed for the GP, so it has consumed nearly the same cpu resources for both test categories. Here are the results I obtained:

The first fitness landscape I’ve used had 2.500 points and the later had a fitness landscape of 6.400 points, here is the source code I’ve used (you just need to turn on the multiprocessing option using the setMultiProcessing method, so Pyevolve will use multiprocessing when you have more than one single core, you can enable the logging feature to check what’s going on behind the scenes):

from pyevolve import *
import math

rmse_accum = Util.ErrorAccumulator()

def gp_sub(a, b): return a-b
def gp_mul(a, b): return a*b
def gp_sqrt(a):   return math.sqrt(abs(a))

def eval_func(chromosome):
global rmse_accum
rmse_accum.reset()
code_comp = chromosome.getCompiledCode()

for a in xrange(0, 80):
for b in xrange(0, 80):
evaluated     = eval(code_comp)
target        = math.sqrt((a*a)+(b*b))
rmse_accum   += (target, evaluated)
return rmse_accum.getRMSE()

def main_run():
genome = GTree.GTreeGP()
genome.setParams(max_depth=4, method="ramped")
genome.evaluator += eval_func
genome.mutator.set(Mutators.GTreeGPMutatorSubtree)

ga = GSimpleGA.GSimpleGA(genome, seed=666)
ga.setParams(gp_terminals       = ['a', 'b'],
gp_function_prefix = "gp")

ga.setMinimax(Consts.minimaxType["minimize"])
ga.setGenerations(20)
ga.setCrossoverRate(1.0)
ga.setMutationRate(0.08)
ga.setPopulationSize(800)
ga.setMultiProcessing(True)

ga(freq_stats=5)
best = ga.bestIndividual()

if __name__ == "__main__":
main_run()

As you can see, the population size was 800 individuals with a 8% mutation rate and a 100% crossover rate for a simple 20 generations evolution. Of course you don’t need so many points in the fitness landscape, I’ve used 2.500+ points to create a cpu intensive fitness function, otherwise, the speedup can be less than 1.0 due the communication overhead between the processes. For the first case (2.500 points fitness landscape) I’ve got a 3.33x speedup and for the last case (6.400 points fitness landscape) I’ve got a 3.28x speedup. The tests were executed in a 2 cores pc (Intel Core 2 Duo).

# New SIGEvolution issue – Volume 3 Issue 4!

Great news, the new issue of SIGEvolution was released, click in the cover above to go to the SIGEvolution site; this issue has an interview is with Hans-Paul Schwefel.

# Approximating Pi number using Genetic Programming

As many (or very few in the real life haha) people know, today is the Pi Approximation Day ! So it’s time to make a contribution to celebrate this funny day =)

My contribution is to use Python and Pyevolve to approximate Pi number using Genetic Programming approach. I’ve created the functions gp_add(+), gp_sub(-), gp_div(/), gp_mul(*) and gp_sqrt (square root) to use as non-terminals of the GP. The fitness function is very simple too, it simple returns the absolute difference between the Python math.pi and the evaluated individual. I’ve used also a population size of 1k individuals with max tree depth of 8 and the random ephemeral constants as random integers. The best approximation I’ve got while running the GP for about 8 minutes (40 generations) was 3.1416185511, best for 3 digits, you can improve it and let it run for more time to get better approximations.

Here is the formulae I’ve got with the GP (click to enlarge):

And here is the output of the script:

Best (0): 3.1577998365
Error: 0.0162071829
Best (10): 3.1417973679
Error: 0.0002047143
Best (20): 3.1417973679
Error: 0.0002047143
Best (30): 3.1417973679
Error: 0.0002047143
Best (40): 3.1416185511
Error: 0.0000258975

- GenomeBase
Score:                   0.000026
Fitness:                 15751.020831

Params:          {'max_depth': 8, 'method': 'ramped'}

Slot [Evaluator] (Count: 1)
Slot [Initializator] (Count: 1)
Name: GTreeGPInitializator - Weight: 0.50
Doc: This initializator accepts the follow parameters:

*max_depth*
The max depth of the tree

*method*
The method, accepts "grow" or "full"

The *GTreeGPInitializator* function.

Slot [Mutator] (Count: 1)
Name: GTreeGPMutatorSubtree - Weight: 0.50
Doc:  The mutator of GTreeGP, Subtree Mutator

The *GTreeGPMutatorSubtree* function

Slot [Crossover] (Count: 1)
Name: GTreeGPCrossoverSinglePoint - Weight: 0.50

- GTree
Height:                 8
Nodes:                  21

GTreeNodeBase [Childs=1] - [gp_sqrt]
GTreeNodeBase [Childs=2] - [gp_div]
GTreeNodeBase [Childs=0] - [26]
GTreeNodeBase [Childs=2] - [gp_div]
GTreeNodeBase [Childs=2] - [gp_mul]
GTreeNodeBase [Childs=2] - [gp_sub]
GTreeNodeBase [Childs=0] - [34]
GTreeNodeBase [Childs=2] - [gp_sub]
GTreeNodeBase [Childs=0] - [44]
GTreeNodeBase [Childs=0] - [1]
GTreeNodeBase [Childs=2] - [gp_mul]
GTreeNodeBase [Childs=0] - [49]
GTreeNodeBase [Childs=0] - [43]
GTreeNodeBase [Childs=1] - [gp_sqrt]
GTreeNodeBase [Childs=0] - [18]
GTreeNodeBase [Childs=0] - [16]
GTreeNodeBase [Childs=0] - [24]
GTreeNodeBase [Childs=0] - [35]

- GTreeGP
gp_sub(44, 1)), gp_mul(49, 43)), gp_sqrt(18)),
16)), gp_add(24, 35)))

And finally, here is the source code:

from __future__ import division
from pyevolve import *
import math

def gp_sub(a, b): return a-b
def gp_div(a, b): return 1 if b==0 else a/b
def gp_mul(a, b): return a*b
def gp_sqrt(a):   return math.sqrt(abs(a))

def eval_func(chromosome):
code_comp = chromosome.getCompiledCode()
ret = eval(code_comp)
return abs(math.pi - ret)

def step_callback(engine):
gen = engine.getCurrentGeneration()
if gen % 10 == 0:
best = engine.bestIndividual()
best_pi = eval(best.getCompiledCode())
print "Best (%d): %.10f" % (gen, best_pi)
print "\tError: %.10f" % (abs(math.pi - best_pi))

return False

def main_run():
genome = GTree.GTreeGP()

genome.setParams(max_depth=8, method="ramped")
genome.evaluator += eval_func

ga = GSimpleGA.GSimpleGA(genome)
ga.setParams(gp_terminals       = ['ephemeral:random.randint(1, 50)'],
gp_function_prefix = "gp")

ga.setMinimax(Consts.minimaxType["minimize"])
ga.setGenerations(50000)
ga.setCrossoverRate(1.0)
ga.setMutationRate(0.09)
ga.setPopulationSize(1000)
ga.stepCallback.set(step_callback)

ga.evolve()
best = ga.bestIndividual()
best.writeDotImage("tree_pi.png")

print best

if __name__ == "__main__":
main_run()

If you are interested why today is the Pi Approximation day, see some resources:

Little Cartoon

Some Background History

Some Pi Approximations

# Genetic Programming and Flex layouts

To show how Genetic Programming of Pyevolve can be flexible, I’ve done a simple example using Adobe Flex and Pyevolve, the example is just to show how to evolve some kind of Flex layouts, I’ve not implemented the fitness function, this example will just create a random Flex layout using MXML. So, here is the Pyevolve code of the example:

import random
from pyevolve import *

def gp_hbox(x, y):
return "%s %s" % (x,y)

def gp_vbox(x, y):
return "%s %s" % (x,y)

def gp_panel(x, y):
return "%s %s" % (x,y)

def eval_func(chromosome):
code_comp = chromosome.getCompiledCode()

for a in xrange(0, 5):
for b in xrange(0, 5):
evaluated     = eval(code_comp)
return random.randint(1,100)

def main_run():
genome = GTree.GTreeGP()
genome.setParams(max_depth=5, method="ramped")
genome.evaluator += eval_func

ga = GSimpleGA.GSimpleGA(genome)

button     = repr("<mx:Button label='Button'/>")
label      = repr("<mx:Label text='Label'/>")
text_input = repr("<mx:TextInput width='50'/>")

ga.setParams(gp_terminals       = [button, label, text_input],
gp_function_prefix = "gp")
ga.setMinimax(Consts.minimaxType["minimize"])
ga.evolve(freq_stats=5)
print ga.bestIndividual()

if __name__ == "__main__":
main_run()

As you can see, I’ve created the layout tags like HBox, VBox and Panel as functions of GP and the Button, Labe, TextInput as terminals of the GP, the result is very funny, it’s just a random layout, but you can use your imagination to create some nice and interesting fitness functions.

Here is the SWF generated from a random individual of the population:

I hope you enjoyed =)

# Genetic Programming meets Python

I’m proud to announce that the new versions of Pyevolve will have Genetic Programming support; after some time fighting with these evil syntax trees, I think I have a very easy and flexible implementation of GP in Python. I was tired to see people giving up and trying to learn how to implement a simple GP using the hermetic libraries for C/C++ and Java (unfortunatelly I’m a Java web developer hehe).

The implementation is still under some tests and optimization, but it’s working nice, here is some details about it:

The implementation has been done in pure Python, so we still have many bonus from this, but unfortunatelly we lost some performance.

The GP core is very very flexible, because it compiles the GP Trees in Python bytecodes to speed the execution of the function. So, you can use even Python objects as terminals, or any possible Python expression. Any Python function can be used too, and you can use all power of Python to create those functions, which will be automatic detected by the framework using the name prefix =)

As you can see in the source-code, you don’t need to bind variables when calling the syntax tree of the individual, you simple use the “getCompiledCode” method which returns the Python compiled function ready to be executed.

Here is a source-code example:

from pyevolve import *
import math

error_accum = Util.ErrorAccumulator()

# This is the functions used by the GP core,
# Pyevolve will automatically detect them
# and the they number of arguments
def gp_sub(a, b): return a-b
def gp_mul(a, b): return a*b
def gp_sqrt(a):   return math.sqrt(abs(a))

def eval_func(chromosome):
global error_accum
error_accum.reset()
code_comp = chromosome.getCompiledCode()

for a in xrange(0, 5):
for b in xrange(0, 5):
# The eval will execute a pre-compiled syntax tree
# as a Python expression, and will automatically use
# the "a" and "b" variables (the terminals defined)
evaluated     = eval(code_comp)
target        = math.sqrt((a*a)+(b*b))
error_accum += (target, evaluated)
return error_accum.getRMSE()

def main_run():
genome = GTree.GTreeGP()
genome.setParams(max_depth=5, method="ramped")
genome.evaluator.set(eval_func)

ga = GSimpleGA.GSimpleGA(genome)
# This method will catch and use every function that
# begins with "gp", but you can also add them manually.
# The terminals are Python variables, you can use the
# ephemeral random consts too, using ephemeral:random.randint(0,2)
# for example.
ga.setParams(gp_terminals       = ['a', 'b'],
gp_function_prefix = "gp")
# You can even use a function call as terminal, like "func()"
# and Pyevolve will use the result of the call as terminal
ga.setMinimax(Consts.minimaxType["minimize"])
ga.setGenerations(1000)
ga.setMutationRate(0.08)
ga.setCrossoverRate(1.0)
ga.setPopulationSize(2000)
ga.evolve(freq_stats=5)

print ga.bestIndividual()

if __name__ == "__main__":
main_run()

I’m very happy and testing the possibilities of this GP implementation in Python.

And of course, everything in Pyevolve can be visualized any time you want (click to enlarge):

The visualization is very flexible too, if you use Python decorators to set how functions will be graphical represented, you can have many interesting visualization patterns. If I change the function “gp_add” to:

@GTree.gpdec(representation="+", color="red")
def gp_add(a, b): return a+b

We’ll got the follow visualization (click to enlarge):

I hope you enjoyed it, I’m currently fixing some bugs, implementing new features, docs and preparing the next release of Pyevolve, which will take some time yet =)