Baby Steps into Genetic Programming
Table of Contents
1 Introduction
While my final ranking in the Google AI Contest was disappointing (280th), it was a very educational experience and was totally offset by Gábor Melis' dominating win using Common Lisp as well.
One of things that piqued my interest during the contest was a post on the AI Challenge forums about a bot written using genetic programming. Genetic programming (GP) and genetic algorithms have always held my interest but seeing the bot in action really motivated me to dive into the matter.
Genetic programming is inspired by biological evolution and is a way of solving problems by setting up an environment (tuned to the problem at hand!) and allowing computer programs to evolve towards a possible solution in that environment.
This article shows my initial exploration into GP using Common Lisp and should be an example of a typical REPL session (my session was a couple of hours divided over two evenings). The code has been reviewed, made a more readable and lispier but still looks very much like what I wrote initially. From the REPL session useful output has been cut and pasted into this article.
This is a really basic introduction into GP and is meant for people interested in GP and/or interested in (Common) Lisp. It is heavy on code and examples and light on theory. My intention is for the reader to play around with the functions on the REPL as they come by in the article.
The code in this article should be portable Common Lisp. If you need to know what a function or macro does please consult the Common Lisp HyperSpec.
1.1 Toy Project
I found A Field Guide to Genetic Programming online and decided to duplicate Bill Clementson's toy project of discovering the formula for calculating the area of a circle. Except for educational purposes we will be starting from scratch.
1.2 Environment
Many Common Lispers use Emacs, Slime and usually Unix or OS X as development environment. However, do not get the impression this is the only way to write Lisp. The development page on CLiki lists many other alternatives in various states of usability.
My recommendation is to use your favorite text editor and CLISP since
it comes with command line history and tab completion built in. Write
or paste the code in your editor and save the file after every
addition as for example "gp.lisp
". Then you only need to type:

(load "/path/to/gp.lisp")
once and can use the command line history to reload it.
2 Generating Random Code
As you might have read, the space.invaders bot was written in PHP. Our weapon of choice is Common Lisp and the first matter of business is generating random code.
2.1 Operators
(defparameter *operators* '(+  * /))
To keep things simple and already knowing our target function (πr²) we opt for using four operators: addition, subtraction, multiplication and division. In the name of simplicity (again) they will have an arity of 2, so they will all take two arguments.
*OPERATORS* is the name of the variable. The *earmuffs* convention is used to signify global variables.
2.2 RANDOMELT
(defun randomelt (sequence) (let ((seqlength (length sequence))) (when (> seqlength 0) (elt sequence (random seqlength)))))
The function RANDOMFORM in the next section uses RANDOMELT to select a random element from the list of operators.
2.3 Generating Random Function Forms
(defun randomform (operators) (append (list (randomelt operators)) (loop repeat 2 ; arity collect (let ((randomnr (random 100))) (cond ((< randomnr 50) (randomform operators)) ((< randomnr 75) (random 10.0)) (t '=input=))))))
The function works as follows: it creates a list of three items. The first item is a random operator and the second and third items are the arguments to that operator:

(operator argument1 argument2)
The arguments are either a random number between 0.0 and 10.0
(arbitrary limit), a variable called =INPUT=
(more about this later)
or, and this is important, another list of three items. For
example:

(* =INPUT= 2.345)
or

(+ 6.789 ( =INPUT= 0.123))
So the RANDOMFORM function will call itself recursively if needed.
The (arbitrary) probabilities for the ARGUMENT1 and ARGUMENT2 are:
50% of the time another random form, 25% of the time a random number
and 25% of the time the =INPUT=
variable.
The process of generating random forms described above is similar to the "grow" method described in Chapter 2.2 of the Field Guide.
2.3.1 Limiting RANDOMFORM
(defun randomform (operators &optional (maxdepth 4)) (append (list (randomelt operators)) (loop repeat 2 ; arity collect (let ((randomnr (random 100))) (if (> maxdepth 0) (cond ((< randomnr 50) (randomform operators ( maxdepth 1))) ((< randomnr 75) (random 10.0)) (t '=input=)) (cond ((< randomnr 50) (random 10.0)) (t '=input=)))))))
Since RANDOMFORM can call itself recursively it can potentially create huge amounts of random code and even exhaust the stack of the Lisp implementation. Hence we need to place some limits on RANDOMFORM. The addition of the MAXDEPTH parameter and the deduction of MAXDEPTH when RANDOMFORM is called again prevents nasty things from happening.
2.4 =INPUT=
Since we want our generated code to calculate the area of a circle
from an input value, the radius, we've added the possibility to
generate the =INPUT=
variable to RANDOMFORM. This variable is the
same as the input argument of the RUNFORM function that is discussed
later.
2.5 Testing RANDOMELT and RANDOMFORM
Pasting RANDOMELT and RANDOMFORM into the REPL and running the latter a couple of times will yield some randomly generated code:

CLUSER> (randomform *operators*) (/ 2.8173013 5.378826) CLUSER> (randomform *operators*) (* =INPUT= =INPUT=) CLUSER> (randomform *operators*) (+ 7.9595613 ( ( =INPUT= (* =INPUT= 0.57189345)) (* ( (/ 6.9174767 9.723027) =INPUT=) (* =INPUT= =INPUT=)))) CLUSER> (randomform *operators*) ( ( =INPUT= ( (/ 2.002045 (* =INPUT= 9.829036)) 4.531122)) =INPUT=)
3 Running Generated Code
(defun runform (form input) (funcall (eval `(lambda (=input=) ,form)) ; note: backquote! input))
To run the generated forms we wrap them in a LAMBDA with an =INPUT=
argument and funcall the evaluated lambda with an input value. The
following REPL interaction shows line by line what happens in
RANDOMFORM:

CLUSER> (defparameter randomform (randomform *operators*)) RANDOMFORM CLUSER> randomform (* 6.341989 =INPUT=) CLUSER> `(lambda (=input=) ,randomform) (LAMBDA (=INPUT=) (* 6.341989 =INPUT=)) CLUSER> (eval `(lambda (=input=) ,randomform)) #<FUNCTION (LAMBDA (=INPUT=)) {AF2F08D}> CLUSER> (funcall (eval `(lambda (=input=) ,randomform)) 2) 12.683978 CLUSER> (* 6.341989 2) 12.683978
However, the generated code can be illegal and cause errors. Since we'll be running a great many pieces of generated code multiple times we don't want to see error messages and drop into the debugger. To keep things simple we let illegal forms return NIL so we can kill them off later. Adding an error handler to RUNFORM will do this:
(defun runform (form input) (let ((*erroroutput* (makebroadcaststream))) (handlercase (funcall (eval `(lambda (=input=) ,form)) input) (error () nil))))
*ERROROUTPUT*
is redirected so we won't be bothered by all kinds of
compilation notices. They only serve to distract in this situation.
3.1 SBCL note
SBCL is a very good, opensource CL implementation but it doesn't compile code very fast. The way we currently handle RUNFORM by recompiling the form every time causes the advancing of generations later in the article to be very slow . Compared to SBCL, which compiles to native code, CLISP, which compiles to bytecode, does the advancing of generations much faster.
Since we're in exploration mode we will leave RUNFORM as it is.
4 Population
(defun createinitialpopulation (operators &optional (size 100)) (loop repeat size collect (randomform operators)))
If you're not familiar with CL: the COLLECT statement in LOOP accumulates a random form on each iteration and returns them as a list once iteration has finished. For example:

CLUSER> (createinitialpopulation *operators* 10) ((* =INPUT= 6.8947406) (* (+ =INPUT= 2.1140647) 6.7930226) (+ (/ (+ (* =INPUT= (/ 8.296512 =INPUT=)) 1.1989254) ( =INPUT= (/ =INPUT= ( 5.1625586 1.4763731)))) =INPUT=) ( 1.5379852 2.1223724) (* =INPUT= ( 3.0632048 (* (+ ( 5.8116364 =INPUT=) 0.1915878) =INPUT=))) (/ 8.237898 =INPUT=) (* (/ (/ ( =INPUT= =INPUT=) 0.56726336) =INPUT=) =INPUT=) (+ (+ 7.147339 =INPUT=) =INPUT=) ( ( =INPUT= =INPUT=) =INPUT=) (+ (/ =INPUT= ( (* =INPUT= =INPUT=) (* ( =INPUT= =INPUT=) =INPUT=))) 6.520609))
The way we're creating the initial population is not ideal. This is because RANDOMFORM only creates one type of syntax tree (using the "grow" method mentioned earlier). For a more diverse initial population we'd need more code generation functions that would output different kinds of syntax trees.
The "full" method is easy to add and is left as an exercise to the reader. See the Field Guide for a description.
5 Fitness
(defun fitness (form fitnessfn testinput) (loop for input in testinput for output = (runform form input) for target = (funcall fitnessfn input) for difference = (when output (abs ( target output))) for fitness = (when output (/ 1.0 (+ 1 difference))) when (null output) do (returnfrom fitness nil) collect fitness into fitnessvalues finally (return (reduce #'* fitnessvalues))))
We need a way to test the output of a generated function against our desired outcome and represent this as a number: the fitness. This is commonly a number between 0.0 and 1.0. The closer to 1.0 the better the fitness.
We also need to check whether the return value of RUNFORM is NIL in which case it executed illegal code. FITNESS will return NIL in that case as well. (Poor man's exception handling but suffices for now.)
So the function needs: 1) the form to check, 2) the fitness function to check against and 3) test input. We also want to check against multiple input values.
The TESTINPUT argument is a list of input values so the most direct approach will be to iterate over these values and running both the form and the fitness functions against them.
To get a fitness value between 0.0 and 1.0 we take the absolute difference between the output of the generated form (OUTPUT) and the output of the fitness function (TARGET). We then add 1 to this difference and use that to divide 1.0. To illustrate:

CLUSER> (defun fit (x) (/ 1.0 (+ 1 x))) FIT CLUSER> (fit 0) ; no difference so the desired output 1.0 CLUSER> (fit 0.1) 0.9090909 CLUSER> (fit 5) 0.16666667 CLUSER> (fit 500) 0.001996008
We're not testing negative values since DIFFERENCE in FITNESS will never be negative.
REPL test of FITNESS:

CLUSER> (defparameter randomform (randomform *operators*)) RANDOMFORM CLUSER> randomform (/ 8.552494 =INPUT=) CLUSER> (fitness randomform (lambda (r) (* pi r r)) '(0 1 2)) NIL ; this is correct since we divided by zero CLUSER> (setf randomform (randomform *operators*)) ( =INPUT= 5.246996) CLUSER> (fitness randomform (lambda (r) (* pi r r)) '(0 1 2)) 9.168484389517398d4 CLUSER> (setf randomform (randomform *operators*)) (* =INPUT= (+ (* =INPUT= 1.8322039) ( =INPUT= 6.812643))) CLUSER> (fitness randomform (lambda (r) (* pi r r)) '(0 1 2)) 0.009196622001630076d0 CLUSER> (fitness '(* pi =input= =input=) (lambda (r) (* pi r r)) '(0 1 2)) 1.0d0 CLUSER> (fitness '(* ( pi 0.1) =input= =input=) (lambda (r) (* pi r r)) '(0 1 2)) 0.649350645706412d0
6 Generation Functions
Now everything is in place to start thinking about creating new generations from the initial population. Currently we will only support crossovers and mutations.
6.1 Traversing Nodes
Before we get to crossovers and mutations we need to look at three important functions: NNODES, RANDOMNODE and REPLACENODE. They all share a common, but for each slightly different, function: TRAVERSENODES.
It is a basic recursive function. It iterates through each element of the FORM list and if that element is a list as well it calls TRAVERSENODES again with that list element as argument.
(defun traversenodesexample (form) (labels ((traversenodes (subform &optional (indent "")) (loop for node in subform do (format t "~D:~A ~S~%" (/ (length indent) 2) indent node) (when (listp node) (traversenodes node (concatenate 'string indent " ")))))) (traversenodes form)))
TRAVERSENODESEXAMPLE goes through FORM exactly how TRAVERSENODES does and for each node it prints its nesting level and the node itself:

CLUSER> (traversenodesexample '(a (b c) (d (e f) g) h)) 0: A 0: (B C) 1: B 1: C 0: (D (E F) G) 1: D 1: (E F) 2: E 2: F 1: G 0: H
6.1.1 NNODES
(defun nnodes (form) (let ((nodes 1)) (labels ((traversenodes (subform) (loop for node in subform do (incf nodes) (when (listp node) (traversenodes node))))) (traversenodes form)) nodes))
Helper function for RANDOMNODE. Returns the number of nodes in FORM
including the root node. Note that "(B C)
" as well as B and C are
counted as nodes:

CLUSER> (nnodes '(b c)) 3 CLUSER> (nnodes '(a (b c) (d (e f) g) h)) 12 CLUSER> (nnodes '()) 1
6.1.2 Picking Random Nodes
We want to be able to pick a random node from a form to perform operations on. RANDOMNODE does this:
(defun randomnode (form) (let* ((index 1) (nodes1 ( (nnodes form) 1)) (randomnodeindex (+ (random nodes1) 1))) (labels ((traversenodes (subform) (loop for node in subform do (when (= index randomnodeindex) (returnfrom randomnode (list :index index :node node))) (incf index) (when (listp node) (traversenodes node))))) (traversenodes form))))
It picks a RANDOMNODEINDEX and starts traversing the nodes of FORM. When it arrives at the index it returns both the node and the index as a property list:

(:index randomnodeindex :node randomnode)
Will not return the root node at index 0.
6.1.3 Replacing Nodes
We need to be able to replace a node in a form with another node. The avoid bugs and confusion we'll let REPLACENODE return this result as a new form.
(defun replacenode (form nodeindex newnode) (let ((index 0)) (labels ((traversenodes (subform) (loop for node in subform do (incf index) when (= index nodeindex) collect newnode when (and (/= index nodeindex) (not (listp node))) collect node when (and (/= index nodeindex) (listp node)) collect (traversenodes node)))) (traversenodes form))))
Traverses the nodes of FORM and collects them to return as a new form. When its INDEX counter is equal to NODEINDEX it collects NEWNODE instead.
The function does not replace the root node (index 0).
6.2 Crossovers
There are different kinds of crossovers but we will be doing the most straightforward one: replace a random node in form A with a random node from form B.
Crossovers are the most used genetic operation in GP and some people have even suggested to never use mutations. They can be compared to the procreation of humans in which attributes of the male and female are represented in their offspring.
(defun crossover (form1 form2 &key (debug nil)) (let ((rnode1 (randomnode form1)) (rnode2 (randomnode form2))) (when debug (format t "form1: ~S~%form2: ~S~%rnode1: ~S~%rnode2: ~S~%" form1 form2 rnode1 rnode2)) (replacenode form1 (getf rnode1 :index) (getf rnode2 :node))))
CROSSOVER takes two forms as arguments and returns a new form. The new form is largely similar to FORM1 but one random node is replaced by a random node from FORM2. The function RANDOMNODE is used to pick these random nodes.

CLUSER> (crossover '(1 (2 3) (4 (5 6) 7) 8) '(a (b c) (d (e f) g) h) :debug t) form1: (1 (2 3) (4 (5 6) 7) 8) form2: (A (B C) (D (E F) G) H) rnode1: (:INDEX 1 :NODE 1) rnode2: (:INDEX 3 :NODE B) (B (2 3) (4 (5 6) 7) 8) CLUSER> (crossover '(1 (2 3) (4 (5 6) 7) 8) '(a (b c) (d (e f) g) h) :debug t) form1: (1 (2 3) (4 (5 6) 7) 8) form2: (A (B C) (D (E F) G) H) rnode1: (:INDEX 3 :NODE 2) rnode2: (:INDEX 7 :NODE (E F)) (1 ((E F) 3) (4 (5 6) 7) 8) CLUSER> (crossover '(1 (2 3) (4 (5 6) 7) 8) '(a (b c) (d (e f) g) h) :debug t) form1: (1 (2 3) (4 (5 6) 7) 8) form2: (A (B C) (D (E F) G) H) rnode1: (:INDEX 2 :NODE (2 3)) rnode2: (:INDEX 6 :NODE D) (1 D (4 (5 6) 7) 8) CLUSER> (crossover '(1 (2 3) (4 (5 6) 7) 8) '(a (b c) (d (e f) g) h) :debug t) form1: (1 (2 3) (4 (5 6) 7) 8) form2: (A (B C) (D (E F) G) H) rnode1: (:INDEX 5 :NODE (4 (5 6) 7)) rnode2: (:INDEX 5 :NODE (D (E F) G)) (1 (2 3) (D (E F) G) 8)
6.3 Mutation
Mutations change a part of a form without needing another form for the operation. It is perhaps analogue to a cosmic ray changing a bit of genetic information in a biological entity.
We will be using subtree mutation since it is straightforward to write and has a potentially large effect.
(defun mutate (form operators &key (debug nil)) (let ((rform (randomform operators)) (rnode (randomnode form))) (when debug (format t "form: ~S~%rform: ~S~%rnode: ~S~%" form rform rnode)) (replacenode form (getf rnode :index) rform)))
Mutation replaces a random node of FORM with a form created by RANDOMFORM.

CLUSER> (mutate '(a (b c) (d (e f) g) h) *operators* :debug t) form: (A (B C) (D (E F) G) H) rform: (+ =INPUT= ( 4.4699216 (+ 9.623513 =INPUT=))) rnode: (:INDEX 3 :NODE B) (A ((+ =INPUT= ( 4.4699216 (+ 9.623513 =INPUT=))) C) (D (E F) G) H) CLUSER> (mutate '(a (b c) (d (e f) g) h) *operators* :debug t) form: (A (B C) (D (E F) G) H) rform: (+ 4.5209084 ( 8.943897 (+ 6.657296 =INPUT=))) rnode: (:INDEX 2 :NODE (B C)) (A (+ 4.5209084 ( 8.943897 (+ 6.657296 =INPUT=))) (D (E F) G) H)
7 Advancing a Generation
7.1 Evaluating a Population
(defun evaluatepopulation (population fitnessfn testinput) (loop for form in population for fitness = (fitness form fitnessfn testinput) when fitness collect (list :fitness fitness :form form) into result finally (return (sort result (lambda (a b) (> (getf a :fitness) (getf b :fitness)))))))
To advance a generation we need to do crossovers and mutations to the forms in the population. The literature suggests we give forms with a higher fitness more chance to be a candidate for a crossover or mutation.
EVALUATEPOPULATION evaluates a population and returns a list of the forms in that population and their fitness (minus any forms that gave errors). For ease of testing on the REPL and for use in our next function the output is sorted from best fitness to worst:

CLUSER> (defparameter population (createinitialpopulation *operators* 10)) POPULATION CLUSER> (evaluatepopulation population (lambda (r) (* pi r r)) '(0 1 2)) ((:FITNESS 0.016815323605722716111d0 :FORM ( ( =INPUT= =INPUT=) =INPUT=)) (:FITNESS 0.009437720627636827027d0 :FORM ( 1.5379852 2.1223724)) (:FITNESS 0.007690745276452599039d0 :FORM (* =INPUT= 6.8947406)) (:FITNESS 0.0031801166742824690382d0 :FORM (* =INPUT= ( 3.0632048 (* (+ ( 5.8116364 =INPUT=) 0.1915878) =INPUT=)))) (:FITNESS 0.0016815216288940106286d0 :FORM (+ (+ 7.147339 =INPUT=) =INPUT=)) (:FITNESS 2.6768631466526024146d4 :FORM (* (+ =INPUT= 2.1140647) 6.7930226)))
This is the same population as created in the Population chapter. Note that four of the ten forms have been eliminated since they executed illegal code.
7.2 HEAD
(defun head (sequence &optional (amount 1)) (if (<= amount 0) nil (if (< (length sequence) amount) sequence (subseq sequence 0 amount))))
Utility function used by ADVANCEGENERATION. Returns AMOUNT elements from the start of SEQUENCE. If SEQUENCE is shorter than AMOUNT it will return the whole SEQUENCE.
7.3 Running the Advancement
(defun advancegeneration (population fitnessfn operators testinput &optional (maxpopulation 100)) (let ((epop (evaluatepopulation population fitnessfn testinput))) (format t "Best fitness of current population: ~S~%" (getf (first epop) :fitness)) (loop for plist in (head epop maxpopulation) for i from 0 for fitness = (getf plist :fitness) for form = (getf plist :form) collect form when (<= (random 1.0d0) fitness) collect (if (<= (random 100) 90) (crossover form (getf (randomelt epop) :form)) (mutate form operators)) ;; Add a new random form to the population now and then. when (<= (random 100) 2) collect (randomform operators))))
Using the list of forms and their fitness from EVALUATEPOPULATION as input we only loop over the first MAXPOPULATION items to keep the population size in check. We then check the form's fitness against a random number between 0.0 and 1.0 and if the number is lower than the fitness the form will be selected for either a crossover (90% chance) or a mutation (10% chance).
So the form's fitness is its chance to be selected. Since the fitness of the forms in the initial population is usually very low it will take quite a few generations before something starts happening.
We collect the form and if it's been selected we also collect the result from either CROSSOVER or MUTATE. These will all be accumulated and eventually returned as the new generation.
There's also a small chance a new random form will be added to the population since the way we're currently handling things it is very possible for a successful form to take over the entire population and make it too homogeneous for good results. This is a quick and easy hack to introduce some chaos into populations.
7.4 Finding a Solution
Lets run ADVANCEGENERATION a hundred times over a new population (we also update the population on each iteration with SETF):

CLUSER> (defparameter population (createinitialpopulation *operators* 100)) POPULATION CLUSER> (loop repeat 100 for i from 0 do (format t "[~S] " i) (setf population (advancegeneration population (lambda (r) (* pi r r)) *operators* '(0 1 2)))) [0] Best fitness of current population: 0.08546627574466652d0 [...] [43] Best fitness of current population: 0.08626213422506805d0 [...] [66] Best fitness of current population: 0.10050107335294403d0 [...]
And again:

CLUSER> (loop repeat 100 for i from 0 do (format t "[~S] " i) (setf population (advancegeneration population (lambda (r) (* pi r r)) *operators* '(0 1 2)))) [...] [26] Best fitness of current population: 0.3865141184760903d0 [...] [93] Best fitness of current population: 0.49884414719455095d0 [...]
Lets do another 300 runs so we've done 500 in total and see what the best form looks like in the end:

CLUSER> (loop repeat 300 for i from 0 do (format t "[~S] " i) (setf population (advancegeneration population (lambda (r) (* pi r r)) *operators* '(0 1 2)))) [...] [35] Best fitness of current population: 0.5397727425319018d0 [...] [62] Best fitness of current population: 0.559234953025742d0 [...] [117] Best fitness of current population: 0.6436990901489741d0 [...] [179] Best fitness of current population: 0.9657338407311192d0 [...] [201] Best fitness of current population: 0.9705968409735506d0 [202] Best fitness of current population: 0.9755729636966523d0 [203] Best fitness of current population: 0.9994359703065686d0 [...] CLUSER> (defparameter bestform (first (evaluatepopulation population (lambda (r) (* pi r r)) '(0 1 2)))) BESTFORM CLUSER> bestform (:FITNESS 0.9994359703065686d0 :FORM (* (+ (+ =INPUT= =INPUT=) (+ =INPUT= (/ =INPUT= (+ 3.5050452 3.5518444)))) =INPUT=)) CLUSER> (runform (getf bestform :form) 0) 0.0 CLUSER> (runform (getf bestform :form) 1) 3.1417055 CLUSER> (runform (getf bestform :form) 2) 12.566822 CLUSER> (runform (getf bestform :form) 3) 28.275349
Here's the output from the real function to calculate the area of a circle for comparison:

CLUSER> (runform '(* pi =input= =input=) 0) 0.0d0 CLUSER> (runform '(* pi =input= =input=) 1) 3.141592653589793d0 CLUSER> (runform '(* pi =input= =input=) 2) 12.566370614359172d0 CLUSER> (runform '(* pi =input= =input=) 3) 28.274333882308138d0
Not quite the impressive result of Bill Clementson's attempt but not half bad either!
8 Conclusion
I hope to have shown how one can start with a few simple functions (RANDOMFORM and its RANDOMELT helper) and explore a topic which might interest you.
If you read the literature you'll notice that there are a lot of things still wrong with the current setup. As you will find out if you start experimenting with the code in this article bloat is a big problem, as well as a successful form taking over the entire population.
If the topic of GP interests you I would suggest you read up and look into:
 eliminating bloat
 improving the creation of the initial population
 improvements to the crossover and mutation functions
for the code in this article.
9 Thanks
Thanks to the denizens of #lisp at freenode for Common Lisp help in general.
Thanks to the following people for proofreading and comments: Marijn Haverbeke, Gábor Melis, David O'Toole and Matthias of the space.invaders team.
Keep in mind that any mistakes are mine.
Labels: ai, artificialintelligence, commonlisp, geneticprogramming, lisp, machinelearning